寻找Algo here并试图实现this code,我得到了不同的l2-norms
作为线性方程的参数向量。我在尝试采用代码时犯了什么错误?
import numpy as np
from scipy import linalg
np.random.seed(123)
v = np.random.rand(4)
A = v[:,None] * v[None,:]
b = np.random.randn(4)
x = linalg.inv(A.T.dot(A)).dot(A.T).dot(b) #Usually not recommended because of Numerical Instability of the Normal Equations https://johnwlambert.github.io/least-squares/
l2_0= linalg.norm(A.dot(x) - b)
print("manually: ", l2_0)
x = linalg.lstsq(A, b)[0]
l2_1= linalg.norm(A.dot(x) - b)
print("scipy.linalg.lstsq: ", l2_1)
# 2-norm of two calculations compared
print(np.allclose(l2_0, l2_1, rtol=1.3e-1))
def direct_ls_svd(x,y):
# append a columns of 1s (these are the biases)
x = np.column_stack([np.ones(x.shape[0]), x])
# calculate the economy SVD for the data matrix x
U,S,Vt = linalg.svd(x, full_matrices=False)
# solve Ax = b for the best possible approximate solution in terms of least squares
x_hat = Vt.T @ linalg.inv(np.diag(S)) @ U.T @ y
#print(x_hat)
# perform train and test inference
#y_pred = x @ x_hat
return y-x @ x_hat #x_hat
x= direct_ls_svd(A, b)
l2_svd= linalg.norm(A.dot(x) - b)
print("svd: ", l2_svd)
# LU
x= linalg.solve(A.T@A, A.T@b)
l2_solve= linalg.norm(A.dot(x) - b)
print("scipy.linalg.solve: ", l2_solve)
# manually: 2.9751344995811313
# scipy.linalg.lstsq: 2.9286130558050654
# True
# svd: 6.830550019041984
# scipy.linalg.solve: 2.928613055805065
字符串
如果我的错误是在求解最小二乘问题的SVD分解算法实现中,或者可能是在Numpy相对Scipy rounding或精度差异中?如何纠正SVD算法,使最小二乘变得与Scipy的可比较?这种算法比迭代最小二乘方法更快,更少的内存消耗?
P.S. SVD应用程序或在这里,SVD for PCA和PLS-SVD是我的最终目标-算法与最小二乘近似相同吗?我对这样的问题感到困惑(即使是代码示例)。有人能为像我这样的新手添加一些清晰度吗?
P.P.S.应用such implementation-我得到了更糟糕的结果:svd: 10.031259300735462
对于l2范数
P.P.P.S.我也缺乏对svd in singular spectral decomposition的一些理解,如果存在ref,作为第一个无监督的降维和第二个非参数TS分析,for practice
P.P.P.P.S.!如果存在多重共线性,则使用PCA进行初步估计,否则可能会得到奇怪的结果(有偏差等).(如果无共线性=>对估计误差不敏感,也可在OLS分析的小条件数中显示,- vc.vs巨大条件数==多变量/多维回归的共线性)
1条答案
按热度按时间jv2fixgn1#
这里最重要的部分是过滤掉
0
的奇异值。对于示例数据S
是[9.22354602e-01 3.92705914e-17 1.10667017e-17 5.55744006e-18]
,请注意,有一个奇异值接近~9.22
,而其他三个奇异值很小(< 1e-16
)。如果试图使用Vt
和U
的一些元素的小误差来重构解,应该是大约相同幅度的误差将被除以这些小的值,并且将对结果增加显著的误差。在这种情况下可以做的是,假设任何足够小的奇异值都是零。在下面的函数修改版本中,我假设所有小于
rcond
乘以最大奇异值的奇异值都应该是零。然后我计算一个掩码m
,并删除U
和Vt
的相应行和列矩阵中。字符串
另一种解决方案是设置
S[m]=0
,这样在最坏的情况下就可以避免额外的副本,但在秩非常低的情况下,您将失去乘法次数减少所带来的潜在节省。