2010-06-07 19 views
11

Je tente de résoudre un ensemble d'équations de la forme Ax = 0. A est une matrice 6x6 connue et j'ai écrit le code ci-dessous en utilisant SVD pour obtenir le vecteur x qui fonctionne pour Dans une certaine mesure. La réponse est à peu près correcte mais pas assez bonne pour être utile pour moi, comment puis-je améliorer la précision du calcul? L'abaissement d'eps en dessous de 1.e-4 entraîne l'échec de la fonction.Calcul de l'espace nul d'une matrice

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) 

Répondre

9

A est de plein rang --- si x est

Comme il semble que vous avez besoin d'une solution des moindres carrés, ce min ||A*x|| s.t. ||x|| = 1, faire la SVD telle que [U S V] = svd(A) et la dernière colonne de V (en supposant que les colonnes sont triées dans l'ordre décroissant des valeurs singulières) est x.

Ie,

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 

Ainsi,

x = 

    0.0050325 
    0.0028357 
    0.00013067 
    0.0052698 
     0.98402 
     0.17788 

Et, ||A*x|| = 0.00027528 par opposition à votre solution précédente pour x||A*x_old|| = 0.079442

+0

x = 0 est une solution au problème , mais inintéressant. La vraie solution au problème, obtenu par différents moyens est: [0.880057009282733,0.571293018023548,0.0664250041765576,1,186.758799941964,33.7579819749057] T – Ainsworth

+0

Êtes-vous sûr? Je vois des éléments non-zéro dans le résultat de 'A * x' ---' [-0.056356 -0.055643 -7.3896e-013 -0.0043278 0.004483 -2.1316e-014] ' – Jacob

+0

À moins bien sûr, vous ne voulez pas l'espace nul, mais la solution des moindres carrés, c'est-à-dire 'min || A * x || s.t. || x || = 1' – Jacob

5

Attention: Il pourrait y avoir confusion avec le SVD en python vs. matlab-syntax (?): en python, numpy.linalg.svd (A) renvoie les matrices u, s, vs uch que u * s * v = A (strictement: point (u, point (diag (s), v) = A, car s est un vecteur et non une matrice 2D en numpy).

La réponse la plus élevée est correcte en ce sens que habituellement vous écrivez u * s * vh = A et vh est retourné, et cette réponse traite de v ET NON vh.

Pour faire une longue histoire courte: si vous avez des matrices u, s, v tel que u * s * v = A, puis les dernières lignes de v, pas les derniers colums de v, décrire la nullspace.

Edit: [pour les gens comme moi:] chacune des dernières lignes est un vecteur v0 tel que A * v0 = 0 (si la valeur singulière correspondante est 0)