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

Question

我的目标是解决:

Kc=y

使用伪逆(即最小范数解):

c=K^{+}y

这样的模型是(希望)高次多项式模型f(x) = sum_i c_i x^i.我特别感兴趣的是欠定的情况,我们有更多的多项式特征而不是数据(几个方程太多的变量/未知数)columns = deg+1 > N = rows.注意K是多项式特征的vandermode矩阵.

我最初使用的是python函数np.linalg.pinv,但后来我注意到有一些时髦的事情正如我在这里所说:为什么在python中解决Xc = y的不同方法会给出不同的解决方案呢？.在那个问题中,我使用方阵来学习[-1.+1]具有高次多项式的区间函数.那里的答案建议我降低多项式的次数和/或增加区间大小.主要问题是我不清楚如何在事物变得不可靠之前选择间隔或最大程度.我认为我的主要问题是选择这样一个数值稳定的范围取决于我可能使用的方法.最后我真正关心的是那个

1. 我使用的方法与该多项式拟合问题的伪逆完全(或非常接近)
2. 它的数值稳定

理想情况下,我想尝试一个大程度的多项式,但这可能受到我的机器精度的限制.是否可以通过使用比浮子更精确的东西来提高机器的数值精度？

另外,我真的关心我使用的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页).因此,也许我的目标应该是使用返回最小最小二乘误差的python函数J.因此,我(在未确定的情况下)运行了三种方法的比较

1. Polygit, np.polyfit
2. 伪逆, np.linalg.pinv
3. 最小二乘, np.linalg.lstsq

并比较了他们给我的数据误差最小二乘误差:

然后我检查了函数似乎经历的奇怪的下降(如果算法不是随机的那么,为什么有逢低,这似乎是一个完全的谜)并且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.pinv.这是最好的事吗？还是我错过了一些明显的东西？

---

作为一个额外的注释,我进入了polyfit代码,看看究竟是什么使它具有更好的最小二乘误差(现在我用它来表示伪逆的最佳近似)并且它似乎有些奇怪条件/数值稳定性代码:

# 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带来了额外的稳定性？

这是否polyfit适用于我的高次多项式线性系统近似任务？

---

在这一点上,我愿意使用像matlab这样的其他软件,如果它为我提供了正确的伪逆和更多的数值稳定性(对于大多数度数和任何边界).

---

我刚才有一个随意的想法是,可能有一种很好的方法来对函数进行采样,这样伪逆的稳定性就很好.我的猜测是用多项式逼近余弦需要某种类型的样本或它们之间的距离(就像奈奎斯特香农采样定理所说,如果基函数是正弦曲线......)

---

应该指出的是,可能反转(或伪陀螺)然后乘法是一个坏主意.看到:

https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/

那个人只讨论逆,但我认为它也延伸到伪逆.

---

现在我的困惑是,通常我们不想明确地计算伪逆,并且做到A^+y=x_min_norm最小的规范解决方案.但是,我原以为这np.lstsq会产生我想要的答案,但它的错误也与其他错误有很大不同.我觉得非常混乱......让我觉得我用错误的方法在python中获得最小的规范解决方案.

---

我不是想要一个正规化的解决方案.我试图获得最小的规范解决方案而不是其他任何东西,尽可能在数值上准确.

Answer 1

我的研究领域涉及一种基本上称为傅立叶扩展的压缩算法.什么是最准确的？由于平滑性的特性,它高度依赖于我认为的载体.在夏天,我使用了一种名为Savitsky Golay的东西.有相当数量稳定和准确的过滤方法.但是,我的顾问有一种相对快速且数值稳定的方法.该区域称为傅立叶扩展或延续.怎么样？我不知道我是否被允许发布,这是文章.如果我相信我可能已经在夏天在这里部分地发布了python.

这与python无关,因为python使用与大多数数字编码方案相同的底层库,即BLAS和LAPACK.Netlib在线.

还有许多其他类似的快速和数字稳定的想法可能适合我推荐.博伊德有一本专门写这本书的书.第6章和第7章就是这样.它是关于正则化的总变化,因为你可能在我想象的信号中产生潜在的噪声.

其他方面.由于病情恶化,您可能需要规范SVD.通常都会有专门的书籍.只需回答您的问题什么是最佳算法.算法有多个维度,您还没有真正说明问题的属性.如果你不知道龙格的现象.那就是使用高次多项式是不利的.

有一整类Hermite多项式来处理Gibbs现象和其他滤波技术,但这并不是很好.您正在使用通用功能.我建议买Trefthen和Bau.有时他们会做Chebychev重投影.

K的条件数是多少.另外,在拟合称为Runges现象的多项式时会发生一些事情.你应该考虑这个.如果条件数太高,请使用您需要的通用函数进行低阶近似以进行正则化.我其实只是读它.您正在使用Vandermonde矩阵.我将很容易地证明这一点.范德蒙德矩阵很糟糕.不要使用它们.他们有结.

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

我尝试了低等级近似,但vandermonde矩阵并不好.看到.

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给出最小二乘法. 它是属性9,A是由SVD构建的.这包括在trefethen中,但vandermonde矩阵是病态的.

即使是小病的vandermonde矩阵也会失去它.现在关于你的解决方案.不要使用vandermonde矩阵.否则构造多项式.更好的想法是重心拉格朗日插值.图书馆在这里

这是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方法.这在Trefethen中有所涉及.SVD 也在 我的quora页面上有一个关于QR数值稳定性的例子,这就是你实际构建这些算法的方法