Python (NumPy, SciPy), znajdowanie pustej przestrzeni macierzy

Question

Próbuję znaleźć pustą przestrzeń (przestrzeń rozwiązania Ax = 0) danej macierzy. Znalazłem dwa przykłady, ale nie mogę sprawić, żeby działały. Ponadto nie mogę zrozumieć, co robią, aby się tam dostać, więc nie mogę debugować. Mam nadzieję, że ktoś mógłby mnie przez to przejść.

Strony dokumentacji (numpy.linalg.svd i numpy.compress) są dla mnie nieprzejrzyste. Nauczyłem się tego, tworząc matrycę C = [A|0], znajdując zredukowaną formę rzutu rzutu i rozwiązując dla zmiennych po rzędzie. Nie wydaje mi się, aby śledzić, jak to się dzieje w tych przykładach.

Dzięki za wszelką pomoc!

Oto moja matryca próbki, która jest taka sama jak wikipedia example:

A = matrix([ 
    [2,3,5], 
    [-4,2,3] 
    ]) 

metodą (found here i here):

import scipy 
from scipy import linalg, matrix 
def null(A, eps=1e-15): 
    u, s, vh = scipy.linalg.svd(A) 
    null_mask = (s <= eps) 
    null_space = scipy.compress(null_mask, vh, axis=0) 
    return scipy.transpose(null_space) 

Kiedy próbuję go, wrócę pusty macierz:

Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56) 
[GCC 4.4.5] on linux2 
Type "help", "copyright", "credits" or "license" for more information. 
>>> import scipy 
>>> from scipy import linalg, matrix 
>>> def null(A, eps=1e-15): 
... u, s, vh = scipy.linalg.svd(A) 
... null_mask = (s <= eps) 
... null_space = scipy.compress(null_mask, vh, axis=0) 
... return scipy.transpose(null_space) 
... 
>>> A = matrix([ 
...  [2,3,5], 
...  [-4,2,3] 
...  ]) 
>>> 
>>> null(A) 
array([], shape=(3, 0), dtype=float64) 
>>> 

Answer 1

Wydaje się działać dobrze dla mnie:

A = matrix([[2,3,5],[-4,2,3],[0,0,0]]) 
A * null(A) 
>>> [[ 4.02455846e-16] 
>>> [ 1.94289029e-16] 
>>> [ 0.00000000e+00]] 

Answer 2

Otrzymujesz dekompozycję SVD macierz A. s jest wektorem wartości własnych. Jesteś zainteresowany prawie zerowymi wartościami własnymi (patrz $ A * x = \ lambda * x $ gdzie $ \ abs (\ lambda) < \ epsilon $), który jest podany przez wektor logicznych wartości null_mask.

Następnie wyodrębnia się z listy vh wektory własne odpowiadające niemal zerowym wartościom własnym, co jest dokładnie tym, czego szukasz: sposobem na rozciągnięcie pustej przestrzeni. Zasadniczo wyodrębniamy wiersze, a następnie transponujemy wyniki, aby uzyskać macierz z wektorami własnymi jako kolumnami.

Answer 3

Twoja metoda jest prawie poprawne. Problem polega na tym, że kształt s zwracany przez funkcję scipy.linalg.svd wynosi (K,) gdzie K = min (M, N). Tak więc w twoim przykładzie s ma tylko dwa wpisy (pojedyncze liczby z pierwszych dwóch wektorów pojedynczych). Następująca korekta twojej funkcji zerowej powinna pozwolić jej działać dla macierzy o dowolnej wielkości.

import scipy 
import numpy as np 
from scipy import linalg, matrix 
def null(A, eps=1e-12): 
... u, s, vh = scipy.linalg.svd(A) 
... padding = max(0,np.shape(A)[1]-np.shape(s)[0]) 
... null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0) 
... null_space = scipy.compress(null_mask, vh, axis=0) 
... return scipy.transpose(null_space) 
A = matrix([[2,3,5],[-4,2,3]]) 
print A*null(A) 
>>>[[ 4.44089210e-16] 
>>> [ 6.66133815e-16]] 
A = matrix([[1,0,1,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]]) 
print null(A) 
>>>[[ 0.   -0.70710678] 
>>> [ 0.   0.  ] 
>>> [ 0.   0.70710678] 
>>> [ 1.   0.  ]] 
print A*null(A) 
>>>[[ 0. 0.] 
>>> [ 0. 0.] 
>>> [ 0. 0.] 
>>> [ 0. 0.]] 

Answer 4

Sympy czyni to prostym.

>>> from sympy import Matrix 
>>> A = [[2, 3, 5], [-4, 2, 3], [0, 0, 0]] 
>>> A = Matrix(A) 
>>> A * A.nullspace()[0] 
Matrix([ 
[0], 
[0], 
[0]]) 
>>> A.nullspace() 
[Matrix([ 
[-1/16], 
[-13/8], 
[ 1]])] 

Answer 5

szybsza, ale mniej numerycznie stabilny metodą jest użycie a rank-revealing QR decomposition, takich jak scipy.linalg.qr z pivoting=True:

import numpy as np 
from scipy.linalg import qr 

def qr_null(A, tol=None): 
    Q, R, P = qr(A.T, mode='full', pivoting=True) 
    tol = np.finfo(R.dtype).eps if tol is None else tol 
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol) 
    return Q[:, rnk:].conj() 

Na przykład:

A = np.array([[ 2, 3, 5], 
       [-4, 2, 3], 
       [ 0, 0, 0]]) 
Z = qr_null(A) 

print(A.dot(Z)) 
#[[ 4.44089210e-16] 
# [ 6.66133815e-16] 
# [ 0.00000000e+00]]