python - python中用于计算最小范数解或从伪逆获得的解的最准确方法是什么？

Question

我的目标是解决：

Kc=y

对于伪逆（即最小范数解）：

c=K^{+}y

这样的模型（希望）是高次多项式模型。我特别感兴趣的是未确定的情况，在这种情况下，我们有比数据更多的多项式特征（很少有方程太多变量/未知量） f(x) = sum_i c_i x^i。注： columns = deg+1 > N = rows是多项式特征的范德模式矩阵。
我最初使用的是python函数 np.linalg.pinv，但后来我注意到了一些奇怪的事情正在发生，正如我在这里注意到的那样： Why do different methods for solving Xc=y in python give different solution when they should not?。在这个问题中，我用一个方阵来学习一个高阶多项式区间上的函数。这里的答案建议我降低多项式的阶数和/或增加区间大小。主要的问题是，在事情变得不可靠之前，我不清楚如何选择间隔或最大程度。我认为我的主要问题是，选择这样一个数值稳定的范围取决于我可以使用的方法。最后我真正关心的是
对于这个多项式拟合问题，我使用的方法与伪逆精确（或非常接近）
它的数值稳定
理想情况下，我想尝试一个大的多项式，但这可能会受到我的机器精度的限制。是否可以使用比浮点数更精确的数值来提高机器的数值精度？
另外，我真的很关心我使用的python中的任何函数，它提供了最接近伪逆的答案（并且希望它的数值稳定，这样我就可以实际使用它）。为了检查伪反转的答案，我编写了以下脚本：

import numpy as np
from sklearn.preprocessing import PolynomialFeatures

def l2_loss(y,y_):
    N = y.shape[0]
    return (1/N)*np.linalg.norm(y-y_)

## some parameters
lb,ub = -200,200
N=100
D0=1
degree_mdl = 120
## target function
freq_cos = 2
f_target = lambda x: np.cos(freq_cos*2*np.pi*x)
## evaluate target_f on x_points
X = np.linspace(lb,ub,N) # [N,]
Y = f_target(X) # [N,]
# get pinv solution
poly_feat = PolynomialFeatures(degree=degree_mdl)
Kern = poly_feat.fit_transform( X.reshape(N,D0) ) # low degrees first [1,x,x**2,...]
c_pinv = np.dot(np.linalg.pinv( Kern ), Y)
## get polyfit solution
c_polyfit = np.polyfit(X,Y,degree_mdl)[::-1] # need to reverse to get low degrees first [1,x,x**2,...]
##
c_lstsq,_,_,_ = np.linalg.lstsq(Kern,Y.reshape(N,1))
##
print('lb,ub = {} '.format((lb,ub)))
print('differences with c_pinv')
print( '||c_pinv-c_pinv||^2 = {}'.format( np.linalg.norm(c_pinv-c_pinv) ))
print( '||c_pinv-c_polyfit||^2 = {}'.format( np.linalg.norm(c_pinv-c_polyfit) ))
print( '||c_pinv-c_lstsq||^2 = {}'.format( np.linalg.norm(c_pinv-c_lstsq) ))
##
print('differences with c_polyfit')
print( '||c_polyfit-c_pinv||^2 = {}'.format( np.linalg.norm(c_polyfit-c_pinv) ))
print( '||c_polyfit-c_polyfit||^2 = {}'.format( np.linalg.norm(c_polyfit-c_polyfit) ))
print( '||c_polyfit-c_lstsq||^2 = {}'.format( np.linalg.norm(c_polyfit-c_lstsq) ))
##
print('differences with c_lstsq')
print( '||c_lstsq-c_pinv||^2 = {}'.format( np.linalg.norm(c_lstsq-c_pinv) ))
print( '||c_lstsq-c_polyfit||^2 = {}'.format( np.linalg.norm(c_lstsq-c_polyfit) ))
print( '||c_lstsq-c_lstsq||^2 = {}'.format( np.linalg.norm(c_lstsq-c_lstsq) ))
##
print('Data set errors')
y_polyfit = np.dot(Kern,c_polyfit)
print( 'J_data(c_polyfit) = {}'.format( l2_loss(y_polyfit,Y) ) )
y_pinv = np.dot(Kern,c_pinv)
print( 'J_data(c_pinv) = {}'.format( l2_loss(y_pinv,Y) ) )
y_lstsq = np.dot(Kern,c_lstsq)
print( 'J_data(c_lstsq) = {}'.format( l2_loss(y_lstsq,Y) ) )

我注意到，很少有 K与 [-1.+1]使用的参数匹配。我知道pinv最终返回伪逆，所以我猜如果我的主要目标是“确保我使用伪逆”，那么使用 polyfit是个好主意。然而，我也从数学上知道，伪逆总是最小化最小二乘误差 pinv无论什么（证明 here定理11.1.2，第446页）。因此，也许我的目标应该是只使用返回最小最小二乘误差的python函数。因此，我（在不确定的情况下）比较了这三种方法
多元土， np.pinv
伪逆， J(c) = || Kc - y ||^2
最小二乘法， J
并比较了他们在数据上给我的误差最小二乘误差：

然后我检查了函数似乎经历的奇怪的下降（如果算法不是随机的，那么这似乎是一个完全的谜，为什么会有下降），对于polyfit，数字通常较小，例如：

lb,ub = (-100, 100)
Data set errors
J_data(c_polyfit) = 5.329753025633029e-12
J_data(c_pinv) = 0.06670557822873546
J_data(c_lstsq) = 0.7479733306782645

考虑到这些结果以及伪逆是最小二乘的极小值，最好的方法似乎是忽略 np.polyfit。这是最好的办法吗？或者我错过了一些明显的东西？
作为一个额外的说明，我进入了 polyfit code来看看究竟是什么使得它有更好的最小二乘误差（现在我用它来表示伪逆的最佳近似值），它似乎有一些奇怪的条件/数值稳定性代码：

# scale lhs to improve condition number and solve
scale = NX.sqrt((lhs*lhs).sum(axis=0))
lhs /= scale
c, resids, rank, s = lstsq(lhs, rhs, rcond)
c = (c.T/scale).T  # broadcast scale coefficients

我假设，是什么给 np.linalg.pinv没有的polyfit带来了额外的稳定性？
对于我的高阶多项式线性系统逼近任务，这是正确的决定吗？
另外，在这一点上，我愿意使用其他软件，如matlab，如果它提供了正确的伪逆和更多的数值稳定性（对于大多数度和任何边界）。
我刚才还有一个随机的想法，那就是也许有一个很好的方法来对函数进行采样，这样伪逆函数的稳定性就很好了。我的猜测是，用多项式近似一个余弦需要一些类型的样本或样本之间的距离（就像尼奎斯特-香农抽样定理所说的如果基函数是正弦函数…）
应该指出的是，可能是颠倒（或伪ivnerting），然后乘以是一个坏主意。见：
https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
那个只讲逆，但我假设它也延伸到了伪逆。
现在我的困惑是，通常我们不想显式地计算伪逆，而要得到最小范数解。然而，我本以为 np.linalg.lstsq会给出我想要的答案，但它的错误也与其他的大不相同。我发现这非常令人困惑……让我觉得我使用了错误的方法来获得Python中的最小范数解决方案。
我不想得到一个正规化的解决方案。我试图得到最小范数解，而没有其他的，尽可能精确的数值计算。

Answer 1

我的研究领域涉及一种压缩算法，本质上称为傅立叶扩展。最准确的是什么？它高度依赖于向量，我相信这是由于平滑的特性。在夏天，我用了一种叫做Savitsky Golay的东西。有相当稳定的数值和精确的过滤方法。然而，我的顾问有一个相对快速和数值稳定的方法。这个区域称为傅立叶延拓或延拓。怎样？我不知道是否允许我发布它，这里是article.如果我相信我可能已经在夏天在这里发布了部分python。
这与python无关，因为python使用的底层库与大多数数字编码方案（blas和lapack）相同。Netlib在线。
有许多其他类似的快速和数字稳定的想法，可能是合适的，我会推荐。有一整本书专门论述这一点，第6章和第7章都是关于这一点的。它是关于在我想象的信号中，由于潜在的噪声，随正则化的总变化。
其他方面。由于条件不好，您可能需要调整SVD。通常都有专门的书。简单地回答你的问题，什么是最好的算法。算法有多个维度，您还没有真正说明问题的性质。如果你不知道y Boyd.这是使用高度多项式是不利的。
有整整一类厄米多项式来处理吉布斯现象和其他滤波技术，但这并不是很好的提出。您使用的是通用函数。我建议你去买汉森和鲍。有时他们会做切比切夫式的回绝。
K的条件数是多少。另外，当拟合称为梯级现象的多项式时，会发生一些事情。你应该考虑一下。如果条件值太高，请使用需要进行低阶近似的通用函数进行正则化。实际上我刚刚读过。你用的是范德蒙矩阵。我将很容易地证明这一点。范德蒙矩阵不好。不要使用它们。Runge's phenomenon.

v = (1:.5:6);

V = vander(v);

c1 = cond(V)

v2 = (1:.5:12);
c2 = cond(vander(v2));
display(c1)
display(c2)

C1=
6.0469E+12型
C2=
9.3987E+32
我尝试了一个低阶近似，但范德蒙矩阵并不好。看。

function B = lowapprox(A)
% Takes a matrix A
% Returns a low rank approx of it
% utilizing the SVD
chi = 1e-10;
[U,S,V]  = svd(A,0);

DS = diag(S);
aa = find(DS > chi);
s= S(aa,aa);
k = length(aa);
Un = U(:,1:k);
Vn = V(:,1:k)';

B = Un*s*Vn;

end


V2  = vander(v2);
r2 = rank(V2);
c2=cond(V2);
B = lowapprox(V2);
c3 = cond(B);
display(c3)
c2 =

   9.3987e+32


c3 =

   3.7837e+32

真的没有什么……如果你不知道当你得到这个倒数时发生了什么，条件数等于最大奇异值超过最小值，所以在机器精度上你有一些非常小的奇异值。
另外，我认为你对最小范数和正则化有些混淆。你说你想要一个最小范数在最小二乘意义上。SVD给出 They have knots.
它的属性9，A是由SVD根据基础构造的。这在Trefetten中有介绍，但vandermonde矩阵是病态的。
即使是构造不良的范德蒙矩阵也会丢失它。关于你的解决方案。不要使用范德蒙矩阵。否则构造多项式。一个更好的想法是重心拉格朗日插值。图书馆 least squares.
下面是matlab中的一个例子。

t= (0:.01:pi)';
f = cos(t);
data = [t,f];
f1 = barylag(data,t)
display(err =norm(f1-f1))
err =

     0

Barylag取自Matlab网站。由于我不能对你的差异做出真正的评论，你应该认识到LSQR的实际方式。LSQR算法是Krylov方法。这是在特雷费顿。SVD is here我在我的Quora页面上有一个关于数值稳定性的例子，它是你实际构造这些算法的方法。