我正在使用 numpy linalg 例程 lstsq 来求解方程组。我的 A 矩阵的大小为 (11046, 504),而我的 B 矩阵的大小为 (11046, 1),确定的秩为 249,因此大约一半的 x 数组求解并不是特别有用。我想使用 s 奇异值数组将求解的与奇异值相对应的参数归零,但似乎 s 数组是按统计显着性递减的顺序排序的。有没有办法可以找出我的哪个 x 对应于每个奇异值 s?
问问题
1547 次
1 回答
3
Mb = x
要获得由 给出的方程的最小二乘解numpy.linalg.lstsq
,您还可以使用numpy.linalg.svd
,它计算奇异值分解M= U S V*
。x
然后给出最佳解x = V Sp U* b
,其中Sp
是 的伪逆S
。给定矩阵U
and V*
(包含矩阵的左奇异向量和右奇异向量M
)和奇异值s
,您可以计算向量z=V*x
。现在,可以在不改变解的情况下任意选择 with 的所有组件,z_i
所有不包含在for中的组件也可以如此。z
i > rank(M)
x_j
z_i
i <= rank(M)
这是一个示例,演示如何使用 Wikipedia 条目中关于奇异值分解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)
于 2013-08-26T20:44:12.517 回答