Scipy sparse invert or spsolve lead to UMFPACK_ERROR_OUT_OF_MEMORY

Question

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?

Answer 1

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.