2010-06-07 60 views
11

我試圖解決一組形式爲Ax = 0的方程組。A是已知的6x6矩陣,我用SVD編寫了下面的代碼以獲得向量x在一定程度上。答案大致正確,但不夠好對我有用,我該如何提高計算精度?將eps降至1.e-4以下會導致該功能失敗。計算矩陣的零空間

from numpy.linalg import * 
from numpy import * 

A = matrix([[0.624010149127497 ,0.020915658603923 ,0.838082638087629 ,62.0778180312547 ,-0.336 ,0], 
[0.669649399820597 ,0.344105317421833 ,0.0543868015800246 ,49.0194290212841 ,-0.267 ,0], 
[0.473153758252885 ,0.366893577716959 ,0.924972565581684 ,186.071352614705 ,-1 ,0], 
[0.0759305208803158 ,0.356365401030535 ,0.126682113674883 ,175.292109352674 ,0 ,-5.201], 
[0.91160934274653 ,0.32447818779582 ,0.741382053883291 ,0.11536775372698 ,0 ,-0.034], 
[0.480860406786873 ,0.903499596111067 ,0.542581424762866 ,32.782593418975 ,0 ,-1]]) 

def null(A, eps=1e-3): 
    u,s,vh = svd(A,full_matrices=1,compute_uv=1) 
    null_space = compress(s <= eps, vh, axis=0) 
    return null_space.T 

NS = null(A) 
print "Null space equals ",NS,"\n" 
print dot(A,NS) 

回答

9

A滿秩---所以x

因爲它看起來像你需要一個最小二乘解,即min ||A*x|| s.t. ||x|| = 1,做SVD這樣[U S V] = svd(A)和的最後一列V(假設列按奇異值遞減排序)爲x

也就是說,

U = 

    -0.23024  -0.23241  0.28225  -0.59968  -0.04403  -0.67213 
     -0.1818  -0.16426  0.18132  0.39639  0.83929  -0.21343 
    -0.69008  -0.59685  -0.18202  0.10908  -0.20664  0.28255 
    -0.65033  0.73984 -0.066702  -0.12447  0.088364  0.0442 
    -0.00045131 -0.043887  0.71552  -0.32745  0.1436  0.59855 
    -0.12164  0.11611  0.5813  0.59046  -0.47173  -0.25029 


S = 

     269.62   0   0   0   0   0 
      0  4.1038   0   0   0   0 
      0   0  1.656   0   0   0 
      0   0   0  0.6416   0   0 
      0   0   0   0  0.49215   0 
      0   0   0   0   0 0.00027528 


V = 

    -0.002597  -0.11341  0.68728  -0.12654  0.70622 0.0050325 
    -0.0024567  0.018021  0.4439  0.85217  -0.27644 0.0028357 
    -0.0036713  -0.1539  0.55281  -0.4961  -0.6516 0.00013067 
     -0.9999 -0.011204 -0.0068651 0.0013713 0.0014128 0.0052698 
    0.0030264  0.17515  0.02341 -0.020917 -0.0054032  0.98402 
    0.012996  -0.96557  -0.15623  0.10603  0.014754  0.17788 

所以,

x = 

    0.0050325 
    0.0028357 
    0.00013067 
    0.0052698 
     0.98402 
     0.17788 

而且,||A*x|| = 0.00027528而不是以前的解決方案x其中||A*x_old|| = 0.079442

+0

X = 0是該問題的解決方案,但一個無趣的。真正解決問題的方法是: [0.880057009282733,0.571293018023548,0.0664250041765576,1,186.758799941964,33.7579819749057] T – Ainsworth 2010-06-07 20:50:03

+0

你確定嗎?我在'A * x' ---'[[-0.056356 -0.055643 -7.3896e-013 -0.0043278 0.004483 -2.1316e-014]' – Jacob 2010-06-07 20:52:05

+0

的結果中看到一些非零元素當然,除非你不想要零空間,但最小二乘解,即'min || A * x || S.T. || X || = 1' – Jacob 2010-06-07 20:53:44

5

注意:有可能是與SVD在python混亂與matlab語法(?): in python,numpy.linalg.svd(A)返回矩陣u,s,vs (嚴格地說:點(u,點(diag(s),v)= A,因爲s是一個向量而不是numpy中的二維矩陣)。

最上面的答案在這個意義上說是正確的通常你寫u * s * vh = A並返回vh,這個答案討論v AND NOT vh。

爲了使長話短說:如果你有矩陣U,S,V,使得U * S * V = A,那麼最後行訴 V的最後colums,描述零空間。

編輯:[對我這樣的人:]每個的最後行的是一個矢量V0,使得A * V0 = 0(若對應的奇異值爲0)