Scipy Sparse invertido o spsolve lleva a UMFPACK_ERROR_OUT_OF_MEMORY

Estoy tratando de invertir una matriz dispersa grande (150000,150000) siguiente manera:

 import scipy as sp import scipy.sparse.linalg as splu #Bs is a large sparse matrix with shape=(150000,150000) #calculating the sparse inverse iBs=splu.inv(Bs) 

lleva al siguiente mensaje de error:

 Traceback (most recent call last): iBs=splu.inv(Bs) File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 134, in spsolve autoTranspose=True) File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/umfpack/umfpack.py", line 603, in linsolve self.numeric(mtx) File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/umfpack/umfpack.py", line 450, in numeric umfStatus[status])) RuntimeError:  failed with UMFPACK_ERROR_out_of_memory 

Reorganicé el progtwig para simplemente resolver un sistema de ecuaciones diferenciales lineales:

 import numpy as np N=Bs.shape[0] I=np.ones(N) M=splu.spsolve(Bs,I) 

y me encuentro con el mismo error otra vez

Estaba usando este código en una máquina con 16 GB de RAM y luego lo moví a un servidor con 32 GB de RAM, aún sin éxito.

¿Alguien encontró esto antes?

Primero, permítame decir que esta pregunta debería formularse mejor en http://scicomp.stackexchange.com, donde hay una gran comunidad de expertos en ciencias computacionales y álgebra lineal numérica.

Empecemos por lo básico: nunca invierta una matriz dispersa, no tiene ningún significado. Vea esta discusión en MATLAB central y particularmente este comentario de Tim Davis.

Brevemente: no hay algoritmos para invertir numéricamente una matriz. Cuando intenta calcular numéricamente la inversa de una matriz NxN, soluciona de hecho N sistemas lineales con N rhs vectores correspondientes a las columnas de la matriz de identidad.

En otras palabras, cuando computas

 from scipy.sparse import eye from scipy.sparse.linalg import (inv, spsolve) N = Bs.shape[0] iBs = inv(Bs) iBs = spsolve(Bs, eye(N)) 

las dos últimas afirmaciones ( inv(eye) y spsolve(Bs, eye(N)) ) son equivalentes. Tenga en cuenta que la matriz de identidad ( eye(N) ) no es un vector de np.ones(N) ), ya que su pregunta supone falsamente.

El punto aquí es que las inversas matriciales rara vez son útiles en álgebra lineal numérica: la solución de Ax = b no se calcula como inv (A) * b, sino por un algoritmo especializado.

Dirigiéndose a su problema específico, para el gran sistema de ecuaciones dispersas no hay solucionadores de caja negra . Puede elegir la clase correcta de solucionadores solo si tiene una buena comprensión de la estructura y las propiedades de su problema de matriz. Las propiedades de sus matrices a su vez son una consecuencia del problema que está tratando de resolver. Por ejemplo, cuando discrimina con el FEM un sistema de PDE elíptico, termina con un sistema simétrico de ecuaciones algebraicas positivas simétricas. Una vez que conozca las propiedades de su problema, puede elegir la estrategia de solución correcta.

En su caso, está tratando de usar un solucionador genérico directo, sin reordenar las ecuaciones. Es bien sabido que esto generará rellenos que destruirán la poca densidad de la matriz iBs en la primera fase de la función spsolve (que debería ser una factorización). Tenga en cuenta que una matriz de doble precisión de 150000 x 150000 requiere aproximadamente 167 GB. de la memoria. Existen muchas técnicas para reordenar ecuaciones con el fin de reducir el relleno durante la factorización, pero no proporciona suficiente información para darle una pista sensata.

Lo siento, pero debería considerar reformular su pregunta en http://scicomp.stackexchange.com indicando claramente cuál es el problema que está tratando de resolver, para dar una pista sobre la estructura y las propiedades de la matriz.

Una matriz dispersa solo encaja en la memoria las entradas que no sean cero de su matriz. Ahora supongamos que haces una inversión. Esto significa que casi todas las entradas de la matriz se convierten en no cero. Las matrices dispersas son memoria optimizada.

Hay algunas operaciones que puede aplicar en matrices dispersas sin perder la propiedad “de repuesto”:

  • Además, solo agregar una constante puede mantener la matriz dispersa dispersa.