我的目标是解决：Kc=y

对于伪逆(即最小范数解)：

^{pr2}$

这样模型(希望)是高次多项式模型f(x) = sum_i c_i x^i。我特别感兴趣的是我们有更多的多项式特征比数据(少方程太多变量/未知数)columns = deg+1 > N = rows。注K是多项式特征的范德模矩阵。在

我最初使用的是python函数np.linalg.pinv，但后来我注意到了一些奇怪的事情，如我所说：Why do different methods for solving Xc=y in python give different solution when they should not?。在这个问题中，我用一个方阵来学习区间[-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(证明here定理11.1.2第446页)。因此，也许我的目标应该是只使用返回最小平方误差J的python函数。因此，我(在不确定的情况下)对这三种方法进行了比较一夫多妻制，np.polyfit

伪逆，np.linalg.pinv

最小二乘法，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 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

我假设，是什么给pinv不具备的额外稳定性带来了？在

用polyfit来做高次多项式线性系统逼近的任务，这是正确的决定吗？在

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

我刚才的另一个随机的想法是，也许有一个很好的方法来取样函数，这样伪逆的稳定性就很好了。我的猜测是用多项式近似余弦需要一定数量的样本或它们之间的距离(就像nyquist-shannon抽样定理所说的，如果基函数是正弦波的话…)

应该指出的是，可能反转(或伪ivnerting)然后再乘以是一个坏主意。参见：

这一个只讨论逆，但我假设它也扩展到伪逆。在

现在我的困惑是，我们通常不想显式地计算伪逆，并做A^+y=x_min_norm来得到最小范数解。然而，我本以为np.lstsq会得到我想要的答案，但它的错误也与其他答案大不相同。我发现这非常令人困惑…让我觉得我用错误的方法得到python中的最小范数解。在

我不是想得到一个正规化的解决方案。我试图得到最小范数解，而不是其他，尽可能精确的数值。在