如何关联哪个奇异值对应于哪个条目

Question

我正在使用numpy-linalg例程lstsq来求解方程组。我的A矩阵的大小为（11046504），而我的B矩阵的大小是（110461），并且确定的秩是249，所以对于x数组求解的大约一半不是特别有用。我想使用奇异值的s数组来将对应于奇异值的参数的求解值归零，但似乎s数组是按统计显著性递减的顺序排序的。有没有一种方法可以算出我的x中的哪一个对应于每个奇异值s？

Answer 1

要获得由numpy.linalg.lstsq给出的方程Mb = x的最小二乘解，还可以使用numpy.linalg.svd，它计算奇异值分解M= U S V*。然后给出最佳解x为x = V Sp U* b，其中Sp是S的伪逆。给定矩阵U和V*（包含矩阵M的左奇异向量和右奇异向量）以及奇异值s，可以计算向量z=V*x。现在，具有i > rank(M)的z的所有分量z_i可以在不改变解决方案的情况下任意选择，因此i <= rank(M)的z_i中不包含的所有分量x_j也可以。

以下示例演示了如何使用维基百科中关于singluar值分解的条目中的样本数据来获得x的重要组件：

import numpy as np
M = np.array([[1,0,0,0,2],[0,0,3,0,0],[0,0,0,0,0],[0,4,0,0,0]])
#We perform singular-value decomposition of M
U, s, V = np.linalg.svd(M)
S = np.zeros(M.shape,dtype = np.float64)
b = np.array([1,2,3,4])
m = min(M.shape)
#We generate the matrix S (Sigma) from the singular values s
S[:m,:m] = np.diag(s)
#We calculate the pseudo-inverse of S
Sp = S.copy()
for m in range(0,m):
  Sp[m,m] = 1.0/Sp[m,m] if Sp[m,m] != 0 else 0
Sp = np.transpose(Sp)
Us = np.matrix(U).getH()
Vs = np.matrix(V).getH()
print "U:n",U
print "V:n",V
print "S:n",S
print "U*:n",Us
print "V*:n",Vs
print "Sp:n",Sp
#We obtain the solution to M*x = b using the singular-value decomposition of the matrix
print "numpy.linalg.svd solution:",np.dot(np.dot(np.dot(Vs,Sp),Us),b)
#This will print:
#numpy.linalg.svd solution: [[ 0.2         1.          0.66666667  0.          0.4       ]]
#We compare the solution to np.linalg.lstsq
x,residuals,rank,s =  np.linalg.lstsq(M,b)
print "numpy.linalg.lstsq solution:",x
#This will print:
#numpy.linalg.lstsq solution: [ 0.2         1.          0.66666667  0.          0.4       ]
#We determine the significant (i.e. non-arbitrary) components of x
Vs_significant = Vs[np.nonzero(s)]
print "Significant variables:",np.nonzero(np.sum(np.abs(Vs_significant),axis = 0))[1]
#This will print:
#Significant variables: [[0 1 2 4]]
#(i.e. x_3 can be chosen arbitrarily without altering the result)