在python中，计算最小范数解或伪逆解的最精确方法是什么？
我的目标是解决：
Kc=y对于伪逆(即最小范数解)：
c=K^{+}y因此，该模型是(希望)高次多项式模型f(x) = sum_i c_i x^i。我特别感兴趣的是未确定的情况，我们比数据有更多的多项式特征(方程少，变量多/未知数)，columns = deg+1 > N = rows。注K是多项式特征的范德模矩阵。
我最初使用的是python函数np.linalg.pinv，但随后我注意到正在发生一些奇怪的事情，如我在这里所指出的：为什么在python中求解Xc=y的不同方法给出了不同的解决方案，而它们不应该这样做？。在这个问题中，我使用一个平方矩阵来学习区间[-1.+1]上的函数，该函数具有高次多项式。那里的答案建议我降低多项式的程度和/或增加区间大小。主要的问题是，在事情变得不可靠之前，我还不清楚如何选择区间或最大度。我认为我的主要问题是，选择这样一个数字稳定的范围取决于我可能使用的方法。最后我真正关心的是
- 对于这个多项式拟合问题，我使用的方法是精确地(或非常接近)伪逆。
- 它的数值稳定
理想情况下，我想尝试一个大度多项式，但这可能受到我的机器精度的限制。是否有可能通过使用比浮子更精确的东西来提高机器的数值精度？
而且，我真的很关心--不管我从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) ) )使用这种方法，我注意到polyfit很少与pinv使用的参数匹配。我知道pinv肯定会返回伪逆，所以我想如果我的主要目标是“确保我使用伪逆”，那么使用np.pinv是个好主意。然而，我在数学上也知道，无论发生什么，伪逆总是最小化最小二乘误差J(c) = || Kc - y ||^2 (证明这里定理11.1.2页446)。因此，也许我的目标应该是只使用返回最小最小二乘误差J的python函数。因此，我(在未确定的情况下)对三种方法进行了比较。
- Polygit，
np.polyfit - 伪逆，
np.linalg.pinv - 最小二乘，
np.linalg.lstsq
并比较了他们给我的数据上的最小二乘误差：
然后，我检查了函数似乎经历的奇怪的倾斜(顺便说一句，如果算法不是随机的，为什么会有dips，这似乎是一个完全的谜团)，而且对于多匹配，数字通常要小一些，例如：
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.pinv。这是最好的做法吗？还是我漏掉了什么明显的东西？
作为额外的注意，我进入多适合码查看究竟是什么使它具有更好的最小二乘误差(现在我用它来表示，它是伪逆的最佳近似)，并且它似乎有一些奇怪的条件/数值稳定性代码：
# 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我猜想，是什么带来了pinv没有的额外的稳定性呢？
这是正确的决定，使用polyfit的任务，我的高次多项式线性系统逼近？
此外，在这一点上，我愿意使用其他软件，如matlab，如果它提供了正确的伪逆和更多的数值稳定性(在大多数程度和任何界限)。
我刚才的另一个随机想法是，也许有一个很好的方法来采样函数，这样伪逆的稳定性是很好的。我的猜测是，用多项式逼近余弦需要某种类型的样本数或它们之间的距离(比如nyquist shannon抽样定理说，如果基函数是正弦波.)
应该指出的是，可能倒置(或假想)然后乘积是一个坏主意。请参见：
https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
它只讨论逆，但我假设它也扩展到伪逆。
现在我的困惑是，通常我们不想显式地计算伪逆，并且做A^+y=x_min_norm来得到最小范数解。然而，我认为np.lstsq会给出我想要的答案，但它的错误与其他错误也有很大的不同。我发现极端的confusing...make，我认为我使用了错误的方法来获得python中的最小范数解决方案。
我不是想得到一个正规化的解决方案。我试图得到最小范数解，而没有其他，尽可能精确的数值。
回答 1
Stack Overflow用户
发布于 2018-01-11 22:05:34
我的研究领域涉及到一种压缩算法，本质上被称为Fourier扩展。什么是最准确的？由于光滑性的性质，它高度依赖于向量。在夏天，我用了一种叫做萨维茨基的东西。有相当稳定的数字和准确的方法过滤这一点。然而，我的顾问有一个相对快速和数值稳定的方法。这一区域称为傅里叶扩张或延拓。多么?我不知道我是否被允许发布，这是文章。，如果我认为我可能已经在夏天张贴在这里的python部分。
这与python无关，因为python使用与大多数数值编码方案相同的底层库，即BLAS和LAPACK。Netlib在线。
还有一些类似的，快速和数字稳定的想法，可能是合适的，我建议。这里有一整本关于b的书，你是博伊德。第6章和第7章就是关于这一点的。这是关于正则化的总变化，因为你在我想象的信号中可能有潜在的噪声。
其他方面。您可能需要规范的SVD由于病态。通常有专门讨论这个问题的书。简单地回答你的问题什么是最好的算法。算法有多个维度，而且您还没有真正说明问题的性质。如果你不知道朗格现象。，那就是使用高次多项式是不利的。
有一整类Hermite多项式可以处理Gibbs现象和其他滤波技术，但这不是很好的估计。你在使用通用函数。我建议你去找特瑞菲和鲍。有时他们会做切比切夫的重新投射。
K的条件数是多少，另外，拟合多项式时会发生什么，叫做Runges现象。你应该考虑一下。使用泛型函数，如果条件数太高，则需要进行低秩近似来进行正则化。其实我刚读过。你用的是范德蒙矩阵。我将很容易地证明这一点。范德蒙矩阵不好。别用它们。他们有结。
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不做任何事情，really...If，你不知道当你得到条件数等于最小的最大奇异值时会发生什么，所以在机器精度上有一些非常小的奇异值。
另外，我认为你对最小范数和正则化有一些混淆。你说过你想要一个最小二乘意义上的最小范数。SVD给出了最小二乘。的性质9，A是由奇异值分解的一个基构造的。这是涵盖在trefethen，但范德蒙矩阵是病态的。
即使是构造不良的小范德蒙矩阵也会失去它。关于你的解决方案。不要使用vandermonde矩阵。否则构造多项式。一个更好的想法是重心拉格朗日插值。库是这里
下面是matlab中的一个例子。
t= (0:.01:pi)';
f = cos(t);
data = [t,f];
f1 = barylag(data,t)
display(err =norm(f1-f1))
err =
0重力赛是从matlab的网站上提取的。由于我不能真正评论您的差异，您应该认识到lsqr的实际操作方式。Lsqr算法是Krylov方法。这是Trefethen所涵盖的。SVD是也是，我在quora页面上有一个关于QR数值稳定性的例子，这就是如何构造这些算法的。
https://stackoverflow.com/questions/46879411
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号