线性方程求解的 SVD 分解

Question

在此处寻找算法并尝试实现此代码，对于线性方程的参数结果向量，我得到了不同的

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

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-algo 的最小二乘法以与 scipy 的进行比较？该算法会比迭代最小二乘法更快且消耗更少的内存吗？

附注 SVD 应用程序，PCA 的 SVD 是我的最终目标——算法与最小二乘近似相同吗？我一般对这样的问题感到困惑（即使有代码示例）。有人可以为像我这样的新手补充一些说明吗？

P.P.S。 应用这样的实现 - 我得到了更糟糕的结果：

svd:  10.031259300735462

对于l2范数

Answer 1

这里最重要的部分是过滤掉

0

的奇异值。对于您的示例数据

S

是

[9.22354602e-01 3.92705914e-17 1.10667017e-17 5.55744006e-18]

，请注意，您有一个接近

~9.22

的奇异值，而其他三个很小 (

< 1e-16

)。 如果您尝试使用

Vt

和

U

的某些元素的小误差来重建解决方案，那么它们的大小应该大致相同，将除以这些小值，并会增加结果的显着误差。

在这种情况下可以做的是，假设任何足够小的奇异值都是精确的零。在函数的以下修改版本中，我假设所有小于

rcond

倍最大奇异值的奇异值应为零。然后我计算一个掩码

m

并删除

U

和

Vt

矩阵的相应行和列。

def direct_ls_svd(A,b,rcond=1e-7):
  # calculate the economy SVD for the data matrix x
  U,S,Vt = linalg.svd(A, full_matrices=False)
  # Mask to remove the zeroes
  m = (abs(S) / np.max(abs(S))) > rcond
  U, S, Vt = U[:,m], S[m], Vt[m, :]
  assert np.allclose( U @ np.diag(S) @ Vt, A)
  # solve Ax = b for the best possible approximate solution in terms of least squares
  x_hat = Vt.T @ ((U.T @ b) / S)
    
  return x_hat

另一种解决方案是设置

S[m]=0

，这样你就可以在最坏的情况下避免额外的副本，但你会失去在极低秩情况下减少乘法次数所带来的潜在节省。