python - 我如何关联哪个奇异值对应于哪个条目？

我正在使用 numpy linalg 例程 lstsq 来求解方程组。我的 A 矩阵的大小为 (11046, 504) 而我的 B 矩阵的大小为 (11046, 1)，并且确定的秩为 249，因此大约一半的 x 数组求解并不是特别有用。我想使用奇异值的 s 数组将与奇异值对应的参数的求解归零，但似乎 s 数组是按统计显着性递减的顺序排序的。有没有办法可以找出我的哪个 x 对应于每个奇异值 s？

最佳答案

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

下面的示例演示了如何使用 singluar value decomposition 上的维基百科条目中的示例数据获取 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)

关于python - 我如何关联哪个奇异值对应于哪个条目？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/18452633/