I am trying to invert a large (150000,150000) sparse matrix as follows:
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)
leads to the following error message:
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: <function umfpack_di_numeric at 0x7f2c76b1d320> failed with UMFPACK_ERROR_out_of_memory
I rejigged the program to simply solve a system of linear differential equations:
import numpy as np
N=Bs.shape[0]
I=np.ones(N)
M=splu.spsolve(Bs,I)
and I encounter the same error again
I was using this code on a machine with 16 GB of RAM and then moved it onto a server with 32 GB of RAM, still to no avail.
has anyone encountered this before?
First let me say that this question should be better asked on http://scicomp.stackexchange.com where there is a great community of experts in computational science and numerical linear algebra.
Let's start from the basics: never invert a sparse matrix, it's completely meaningless. See this discussion on MATLAB central and particularly this comment by Tim Davis.
Briefly: there are no algorithms for numerically inverting a matrix. Whenever you try to numerically compute the inverse of a NxN matrix, you solve in fact N linear systems with N rhs vectors corresponding to the columns of the identity matrix.
In other words, when you compute
from scipy.sparse import eye
from scipy.sparse.linalg import (inv, spsolve)
N = Bs.shape[0]
iBs = inv(Bs)
iBs = spsolve(Bs, eye(N))
the last two statements (inv(eye) and spsolve(Bs, eye(N))) are equivalent. Please note that the identity matrix (eye(N)) is not a ones vector (np.ones(N)) as you question falsely assumes.
The point here is that matrix inverses are seldom useful in numerical linear algebra: the solution of Ax = b is not computed as inv(A)*b, but by a specialised algorithm.
Going to your specific problem, for big sparse system of equations there are no black-box solvers. You can chose the correct class of solvers only if you have a good understanding of the structure and properties of your matrix problem. The properties of your matrices in turn are a consequence of the problem you are trying to solve. E.g. when you discretise by the FEM a system of elliptic PDE, you end up with a symmetric positive sparse system of algebraic equations. Once you know the properties of your problem, you can choose the correct solving strategy.
In your case, you are trying to use a generic direct solver, without reordering the equations. It is well known that this will generate fill-ins which destroy the sparsity of the iBs matrix in the first phase of the spsolve function (which should be a factorisation.) Please note that a full double precision 150000 x 150000 matrix requires about 167 GB of memory. There are a lot of techniques for reordering equations in order to reduce the fill-in during factorisation, but you don't provide enough info for giving you a sensible hint.
I'm sorry, but you should considering reformulating your question on http://scicomp.stackexchange.com clearly stating which is the problem you are trying to solve, in order to give a clue on the matrix structure and properties.
A sparse array only fits the non-zero entries of your matrix into memory. Now suppose you do an inversion. This means that almost all entries of the matrix become non-zero. Sparse matrices are memory optimized.
There are some operations which you can apply on sparse matrices without losing the "spare" property:
- Addition, just adding a constant can keep the sparse matrix sparse.
