2

我想在 12*12 矩阵上执行 SVD。numpy.linalg.svd作品很好。但是当我尝试通过执行 u*s*v 来取回 12*12 矩阵 A 时，我没有取回它。

import cv2

import numpy as np

import scipy as sp

from scipy import linalg, matrix

a_matrix=np.zeros((12,12))

with open('/home/koustav/Documents/ComputerVision/A2/codes/Points0.txt','r') as f:

    for (j,line) in enumerate(f):
         i=2*j
         if(i%2==0):
        values=np.array(map(np.double,line.strip('\n').split(' ')))
        a_matrix[i,4]=-values[2]
        a_matrix[i,5]=-values[3]
        a_matrix[i,6]=-values[4]
        a_matrix[i,7]=-1
        a_matrix[i,8]=values[1]*values[2]
        a_matrix[i,9]=values[1]*values[3]
        a_matrix[i,10]=values[1]*values[4]
        a_matrix[i,11]=values[1]*1

        a_matrix[i+1,0]=values[2]
        a_matrix[i+1,1]=values[3]
        a_matrix[i+1,2]=values[4]
        a_matrix[i+1,3]=1
        a_matrix[i+1,8]=-values[0]*values[2]
        a_matrix[i+1,9]=-values[0]*values[3]
        a_matrix[i+1,10]=-values[0]*values[4]
        a_matrix[i+1,11]=-values[0]*1


s_matrix=np.zeros((12,12))

u, s, v = np.linalg.svd(a_matrix,full_matrices=1)

k=0

while (k<12):

    s_matrix[k,k]=s[k]

    k+=1
print u

print '\n'

print s_matrix

print '\n'

print (u*s_matrix*v)

这些是我使用的要点：

285.12 14.91 2.06655 -0.807071 -6.06083

243.92 100.51 2.23268 -0.100774 -5.63975

234.7 176.3 2.40898 0.230613 -5.10977

-126.59 -152.59 -1.72487 4.96296 -10.4564

-173.32 -164.64 -2.51852 4.95202 -10.3569

264.81 28.03 2.07303 -0.554853 -6.05747

请提出一些建议...

4

1 回答 1

3

除了通过使用像这样的内置函数来节省一些代码和时间之外numpy.diag，您的问题似乎是*操作员。在 numpy 中，您必须numpy.dot用于矩阵乘法。请参阅下面的代码以获取工作示例...

In [16]: import numpy as np

In [17]: A = np.arange(15).reshape(5,3)

In [18]: A
Out[18]: 
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])


In [19]: u, s, v = np.linalg.svd(A)

In [20]: S = np.diag(s)

In [21]: S = np.vstack([S, np.zeros((2,3)) ])

In [22]: #fill in zeros to get the right shape

In [23]: np.allclose(A, np.dot(u, np.dot(S,v)))
Out[23]: True

numpy.allclose检查两个数组在数值上是否接近...

于 2013-02-16T11:17:25.057 回答