我的目标是解决：

Kc=y

c=K^{+}y

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) ) )

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

# 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

最佳答案

我的研究领域涉及一种压缩算法，本质上称为傅立叶扩展。最准确的是什么？它高度依赖于向量，我相信这是由于平滑的特性。在夏天，我用了一种叫做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)

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

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

     0

关于python - python中用于计算最小范数解或从伪逆获得的解的最准确方法是什么？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/46879411/