C'è un modo per invertire in modo efficiente un array di matrici con numpy?

Question

Normalmente invertirò una matrice di matrici 3x3 in un ciclo for come nell'esempio seguente. Purtroppo i loop for sono lenti. Esiste un modo più rapido ed efficiente per fare ciò?

import numpy as np 
A = np.random.rand(3,3,100) 
Ainv = np.zeros_like(A) 
for i in range(100): 
    Ainv[:,:,i] = np.linalg.inv(A[:,:,i]) 

Answer 1

Si scopre che si sta bruciando due livelli nel codice numpy.linalg. Se osservi numpy.linalg.inv, puoi vedere che è solo una chiamata a numpy.linalg.solve (A, inv (A.shape [0]). Questo ha l'effetto di ricreare la matrice di identità in ogni iterazione del tuo for loop. Poiché tutti i tuoi array hanno le stesse dimensioni, è una perdita di tempo. Salta questo passo pre-allocando la matrice delle identità shaves ~ 20% di sconto (fast_inverse). I miei test suggeriscono che pre-allocare l'array o allocare da un elenco di risultati non fa molta differenza

Guardare un livello più in profondità e trovare la chiamata alla routine di lapazi, ma è racchiusa in diversi controlli di integrità. Se si eliminano tutti questi e basta chiamare lapack in il tuo ciclo for (visto che conosci già le dimensioni della tua matrice e forse sai che è reale, non complessa), le cose funzionano MOLTO più velocemente (Nota che ho ingrandito il mio array):

import numpy as np 
A = np.random.rand(1000,3,3) 
def slow_inverse(A): 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.inv(A[i]) 
    return Ainv 

def fast_inverse(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.solve(A[i], identity) 
    return Ainv 

def fast_inverse2(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 

    return array([np.linalg.solve(x, identity) for x in A]) 

from numpy.linalg import lapack_lite 
lapack_routine = lapack_lite.dgesv 
# Looking one step deeper, we see that solve performs many sanity checks. 
# Stripping these, we have: 
def faster_inverse(A): 
    b = np.identity(A.shape[2], dtype=A.dtype) 

    n_eq = A.shape[1] 
    n_rhs = A.shape[2] 
    pivots = zeros(n_eq, np.intc) 
    identity = np.eye(n_eq) 
    def lapack_inverse(a): 
     b = np.copy(identity) 
     pivots = zeros(n_eq, np.intc) 
     results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) 
     if results['info'] > 0: 
      raise LinAlgError('Singular matrix') 
     return b 

    return array([lapack_inverse(a) for a in A]) 


%timeit -n 20 aI11 = slow_inverse(A) 
%timeit -n 20 aI12 = fast_inverse(A) 
%timeit -n 20 aI13 = fast_inverse2(A) 
%timeit -n 20 aI14 = faster_inverse(A) 

I risultati sono impressionanti:

20 loops, best of 3: 45.1 ms per loop 
20 loops, best of 3: 38.1 ms per loop 
20 loops, best of 3: 38.9 ms per loop 
20 loops, best of 3: 13.8 ms per loop 

EDIT: non ho guardato abbastanza da vicino a ciò che viene restituito nella risolvere. Risulta che la matrice 'b' viene sovrascritta e contiene il risultato alla fine. Questo codice ora dà risultati coerenti.

Answer 2

Per i cicli non sono infatti necessariamente molto più lento rispetto alle alternative e anche in questo caso, non vi aiuterà molto. Ma qui è un suggerimento:

import numpy as np 
A = np.random.rand(100,3,3) #this is to makes it 
          #possible to index 
          #the matrices as A[i] 
Ainv = np.array(map(np.linalg.inv, A)) 

Timing questa soluzione contro la soluzione produce un piccolo ma evidente differenza:

# The for loop: 
100 loops, best of 3: 6.38 ms per loop 
# The map: 
100 loops, best of 3: 5.81 ms per loop 

ho cercato di utilizzare la routine NumPy 'Vectorize' con la speranza di creare un soluzione ancora più pulita, ma dovrò dare una seconda occhiata a questo. Il cambiamento di ordine nell'array A è probabilmente il cambiamento più significativo, poiché utilizza il fatto che gli array di numpy sono ordinati in termini di colonne e quindi una lettura lineare dei dati è sempre leggermente più veloce in questo modo.