¿Hay alguna manera de invertir eficientemente una matriz de matrices con numpy?

Question

Normalmente invertiría una matriz de matrices 3x3 en un bucle for como en el ejemplo siguiente. Desafortunadamente, los bucles for son lentos. ¿Hay una manera más rápida y más eficiente de hacer esto?

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

Resulta que se está quemando dos niveles hacia abajo en el código numpy.linalg. Si miras en numpy.linalg.inv, puedes ver que es solo una llamada a numpy.linalg.solve (A, inv (A.shape [0]). Esto tiene el efecto de recrear la matriz de identidad en cada iteración de tu Debido a que todas sus matrices tienen el mismo tamaño, eso es una pérdida de tiempo. Si omite este paso al asignar previamente la matriz de identidad, se afeitará ~ 20% de descuento en el tiempo (fast_inverse). Mi prueba sugiere que se preasigne la matriz o que se asigne de una lista de resultados no hace mucha diferencia.

Mire un nivel más profundo y encontrará la llamada a la rutina lapack, pero está envuelto en varios controles de cordura. Si los quita todos y simplemente llame a lapack en su bucle for (dado que ya conoce las dimensiones de su matriz y tal vez sepa que es real, no complejo), las cosas se ejecutan MUCHO más rápido (Tenga en cuenta que he hecho mi matriz más grande):

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) 

Los resultados son impresionantes:

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: no me veía con suficiente atención a lo que se devuelve en resolver. Resulta que la matriz 'b' se sobrescribe y contiene el resultado al final. Este código ahora da resultados consistentes.

Answer 2

Para los bucles de hecho no son necesariamente mucho más lentos que las alternativas y también en este caso, no le ayudará mucho. Pero aquí es una sugerencia:

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 esta solución frente a la solución produce una pequeña pero notable diferencia:

# 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 

He intentado utilizar la rutina numpy 'vectorize' con la esperanza de crear una incluso una solución más limpia, pero tendré que echarle un segundo vistazo a eso. El cambio de orden en la matriz A es probablemente el cambio más significativo, ya que utiliza el hecho de que las matrices numpy se ordenan en columnas y, por lo tanto, una lectura lineal de los datos es ligeramente más rápida de esta manera.