2013-02-27 16 views
5

Muszę odwrócić dużą macierz rzadką. Nie mogę uciec od inwersji matrycy, jedynym skrótem byłoby po prostu zorientowanie się w głównych przekątnych elementach i zignorowanie elementów poza diagonalnych (wolałbym nie, ale jako rozwiązanie byłoby do przyjęcia).Odwracanie dużych rzadkich macierzy za pomocą scipy

Matryce muszę odwracające są zazwyczaj duże (40000 * 40000) i które mają tylko kilka non-niezerowych przekątnych. Moje obecne podejście jest budowanie wszystko skąpe, a następnie

posterior_covar = np.linalg.inv (hessian.todense())

to wyraźnie zajmuje dużo czasu i dużo pamięci.

żadnych wskazówek, czy jest to tylko kwestia cierpliwości lub podejmowania problemu mniejsze?

+0

Opcja 'scipy' doc mówi, że nie trzeba do zagęszczenia matrycę, więc jestem nieco zdezorientowany. http://docs.scipy.org/doc/scipy-0.8.x/reference/sparse.html – BenDundee

+3

Wersja 0.12 scipy (aktualnie w fazie testów beta) ma funkcję scipy.sparse.linalg.inv. –

Odpowiedz

5

Nie sądzę, że moduł rzadki ma wyraźną metodę odwrotną, ale ma rzadkie rozwiązania. Działa coś podobnego do tego zabawkowego przykładu:

>>> a = np.random.rand(3, 3) 
>>> a 
array([[ 0.31837307, 0.11282832, 0.70878689], 
     [ 0.32481098, 0.94713997, 0.5034967 ], 
     [ 0.391264 , 0.58149983, 0.34353628]]) 
>>> np.linalg.inv(a) 
array([[-0.29964242, -3.43275347, 5.64936743], 
     [-0.78524966, 1.54400931, -0.64281108], 
     [ 1.67045482, 1.29614174, -2.43525829]]) 

>>> a_sps = scipy.sparse.csc_matrix(a) 
>>> lu_obj = scipy.sparse.linalg.splu(a_sps) 
>>> lu_obj.solve(np.eye(3)) 
array([[-0.29964242, -0.78524966, 1.67045482], 
     [-3.43275347, 1.54400931, 1.29614174], 
     [ 5.64936743, -0.64281108, -2.43525829]]) 

Pamiętaj, że wynik został transponowany!

Jeśli spodziewasz się, że twoja odwrotność będzie również rzadka, a gęsty powrót z ostatniego rozwiązania nie będzie pasował do pamięci, możesz również wygenerować go w jednym wierszu (kolumnie) na raz, wyodrębnić wartości niezerowe, rzadki i zbudować macierz odwrotną od tych:

>>> for k in xrange(3) : 
...  b = np.zeros((3,)) 
...  b[k] = 1 
...  print lu_obj.solve(b) 
... 
[-0.29964242 -0.78524966 1.67045482] 
[-3.43275347 1.54400931 1.29614174] 
[ 5.64936743 -0.64281108 -2.43525829] 
Powiązane problemy