Inversione di matrici sparse di grandi dimensioni con scipy

Question

Devo invertire una matrice sparsa di grandi dimensioni. Non riesco a sfuggire all'inversione della matrice, l'unica scorciatoia potrebbe essere solo avere un'idea degli elementi diagonali principali, e ignorare gli elementi fuori diagonale (preferirei di no, ma come soluzione sarebbe accettabile).

Le matrici che ho bisogno di invertire sono generalmente grandi (40000 * 40000) e hanno solo una manciata di diagonali non diverse da zero. Il mio approccio attuale è quello di costruire tutto radi, e poi

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

questo richiede chiaramente molto tempo e un sacco di memoria.

Qualche suggerimento o è solo questione di pazienza o di rendere il problema più piccolo?

Answer 1

Non penso che il modulo sparse abbia un metodo inverso esplicito, ma ha solutori sparsi. Qualcosa come questo esempio di giocattolo funziona:

>>> 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]]) 

Nota che il risultato è trasposto!

Se si prevede che la propria inversa sia sparsa e il ritorno denso dall'ultima soluzione non si adatta alla memoria, è anche possibile generarlo una riga (colonna) alla volta, estrarre i valori diversi da zero, e costruire la matrice inversa sparse da quelle:

>>> 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]