Multiplicando matrices dispersas y densas de Numpy / Scipy de manera eficiente

Estoy trabajando para implementar la siguiente ecuación:

X =(YT * Y + YT * C * Y) ^ -1 

Y es una matriz (nxf) y C es (nxn) una diagonal; n es aproximadamente 300k y f variará entre 100 y 200. Como parte de un proceso de optimización, esta ecuación se utilizará casi 100 millones de veces, por lo que debe procesarse muy rápido.

Y se inicializa al azar y C es una matriz muy dispersa con solo unos pocos números de los 300 k en la diagonal será diferente a 0. Dado que las funciones diagonales de Numpy crean matrices densas, creé C como una matriz csr dispersa. Pero al intentar resolver la primera parte de la ecuación:

 r = dot(C, Y) 

La computadora se bloquea debido a los límites de memoria. Entonces decidí intentar convertir Y en csr_matrix y hacer la misma operación:

 r = dot(C, Ysparse) 

Y este enfoque tomó 1.38 ms . Pero esta solución es un tanto “complicada” ya que estoy usando una matriz dispersa para almacenar una densa, me pregunto cuán eficiente es esto realmente.

Entonces, mi pregunta es si hay alguna forma de multiplicar la C dispersa y la D densa sin tener que convertir la Y en escasa y mejorar el rendimiento. Si de alguna manera C pudiera representarse como densa diagonal sin consumir toneladas de memoria, tal vez esto llevaría a un rendimiento muy eficiente, pero no sé si esto es posible.

¡Aprecio tu ayuda!

La razón por la que el producto punto se encuentra con problemas de memoria al calcular r = punto (C, Y) es porque la función punto de numpy no tiene soporte nativo para manejar matrices dispersas. Lo que está sucediendo es que Number piensa en la matriz dispersa C como un objeto python, y no en una matriz numpy. Si inspeccionas a pequeña escala, puedes ver el problema de primera mano:

 >>> from numpy import dot, array >>> from scipy import sparse >>> Y = array([[1,2],[3,4]]) >>> C = sparse.csr_matrix(array([[1,0], [0,2]])) >>> dot(C,Y) array([[ (0, 0) 1 (1, 1) 2, (0, 0) 2 (1, 1) 4], [ (0, 0) 3 (1, 1) 6, (0, 0) 4 (1, 1) 8]], dtype=object) 

Claramente, lo anterior no es el resultado que le interesa. En cambio, lo que desea hacer es calcular utilizando la función sparse.csr_matrix.dot de scipy:

 r = sparse.csr_matrix.dot(C, Y) 

o mas compacto

 r = C.dot(Y) 

Tratar:

 import numpy as np from scipy import sparse f = 100 n = 300000 Y = np.random.rand(n, f) Cdiag = np.random.rand(n) # diagonal of C Cdiag[np.random.rand(n) < 0.99] = 0 # Compute YT * C * Y, skipping zero elements mask = np.flatnonzero(Cdiag) Cskip = Cdiag[mask] def ytcy_fast(Y): Yskip = Y[mask,:] CY = Cskip[:,None] * Yskip # broadcasting return Yskip.T.dot(CY) %timeit ytcy_fast(Y) # For comparison: all-sparse matrices C_sparse = sparse.spdiags([Cdiag], [0], n, n) Y_sparse = sparse.csr_matrix(Y) %timeit Y_sparse.T.dot(C_sparse * Y_sparse) 

Mis tiempos:

 In [59]: %timeit ytcy_fast(Y) 100 loops, best of 3: 16.1 ms per loop In [18]: %timeit Y_sparse.T.dot(C_sparse * Y_sparse) 1 loops, best of 3: 282 ms per loop 

Primero, ¿está realmente seguro de que necesita realizar una inversión de matriz completa en su problema? La mayoría de las veces, solo se necesita calcular x = A ^ -1 y que es un problema mucho más fácil de resolver.

Si esto es realmente así, consideraría calcular una aproximación de la matriz inversa en lugar de la inversión de la matriz completa. Dado que la inversión de la matriz es realmente costosa. Vea, por ejemplo, el algoritmo de Lanczos para una aproximación eficiente de la matriz inversa. La aproximación se puede almacenar escasamente como un bono. Además, solo requiere operaciones de vector de matriz, por lo que ni siquiera tiene que almacenar la matriz completa para invertir.

Como alternativa, al utilizar pioneros, también puede usar el método .todense para calcular la matriz para invertir utilizando operaciones de vector de matriz eficientes. Hay un contenedor especial disperso para matrices diagonales.

Para una implementación del algoritmo de Lanczos, puede echar un vistazo a pyoperators (descargo de responsabilidad: soy uno de los coautores de este software).