Te sugerimos que pruebes esta solución en un ambiente controlado antes de enviarlo a producción, un saludo.
Solución:
Primero déjeme decir que esta pregunta debería hacerse mejor en http://scicomp.stackexchange.com donde hay una gran comunidad de expertos en ciencia computacional y álgebra lineal numérica.
Empecemos por lo básico: nunca invertir una matriz dispersa, no tiene ningún sentido. Vea esta discusión en MATLAB central y particularmente este comentario de Tim Davis.
Brevemente: no existen algoritmos para invertir numéricamente una matriz. Cada vez que intenta calcular numéricamente la inversa de una matriz NxN, resuelve de hecho N sistemas lineales con N vectores rhs correspondientes a las columnas de la matriz identidad.
En otras palabras, cuando calculas
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)
) es no un vector de unidades (np.ones(N)
) como usted pregunta asume falsamente.
El punto aquí es que las matrices inversas rara vez son útiles en álgebra lineal numérica: la solución de Ax = b no se calcula como inv(A)*b, sino mediante un algoritmo especializado.
Yendo a su problema específico, para un gran sistema de ecuaciones disperso no hay caja negra solucionadores Puede elegir la clase correcta de solucionadores solo si tiene una buena comprensión de la estructura y las propiedades de su problema matricial. Las propiedades de sus matrices, a su vez, son una consecuencia del problema que está tratando de resolver. Por ejemplo, cuando discretizas mediante el FEM un sistema de PDE elíptica, terminas con un sistema disperso positivo simétrico de ecuaciones algebraicas. Una vez que conoce las propiedades de su problema, puede elegir la estrategia de solución correcta.
En su caso, está tratando de usar un solucionador directo genérico, sin reordenar las ecuaciones. Es bien sabido que esto generará rellenos que destruirán la escasez de iBs
matriz en la primera fase del spsolve
(que debería ser una factorización). Tenga en cuenta que una matriz completa de doble precisión de 150000 x 150000 requiere aproximadamente 167 GB de memoria. Hay 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.
Nos puedes añadir valor a nuestra información colaborando tu veteranía en las observaciones.