14

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?

denis
  • 21,378
  • 10
  • 65
  • 88
laila
  • 1,009
  • 3
  • 15
  • 27
  • I have the exact same problem on a PC running Linux Mint 17.3 with 64 Gb of RAM. I checked that maximum ram consumption was less than 20.3 Gb thus I assume there is some limit of memory consumption related to software issue... I checked "ulimit" is "unlimited"... No clue what to do there... – Alain Dec 19 '15 at 01:57
  • It might be good to add that this seems to be linux specific, I did run successfully the exact same code on Windows 7 (both systems with Python 2.7 64 bits). – Alain Dec 19 '15 at 02:06
  • I found this link mentioning a possible re-build of UMFPACK ( http://scicomp.stackexchange.com/questions/7973/memory-management-for-solving-large-sparse-systems-with-umfpack ). But I do not know what to do with it. – Alain Dec 19 '15 at 02:21
  • thanks @Alain , if I figure this out will share the solution :) – laila Dec 21 '15 at 10:29
  • I was wondering if you found a way out ? I have been looking at it on my side but I could not find any explanation for this limitation... – Alain Jan 06 '16 at 17:24
  • @Alain http://scicomp.stackexchange.com/questions/7973/memory-management-for-solving-large-sparse-systems-with-umfpack hints at 32/64 builds of UMFPACK. I do not think that this is the problem here: I would suggest excessive fill-in during factorisation. – Stefano M Jan 23 '16 at 23:36
  • @laila in case you did not solve your problem, I recently investigated the bug I had, which seems very similar to yours. I found out that my Linux distribution was not optimized for my recent i7 processors and this was the origin of the extra RAM memory consumption. You may find more details here : http://stackoverflow.com/questions/36344049/scipy-sparse-linalg-spsolve-surprising-behaviour-for-large-sparse-matrices-on-li After installing Fedora 23, everything runs smoothly. – Alain Apr 05 '16 at 20:05

2 Answers2

4

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.

Stefano M
  • 4,267
  • 2
  • 27
  • 45
0

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.
www.data-blogger.com
  • 4,076
  • 7
  • 43
  • 65