[SciPy-User] solving large sparse linear system with Laplacian matrix
Fri Oct 30 02:46:48 CDT 2009
I'm looking for a memory-savvy algorithm in scipy.sparse.linalg, that I
could use to solve a very large sparse linear system. The system can be
written as AX = B, where A is a symmetric band matrix with non-zero
elements on only 4 or 6 upper (and lower) diagonals. The shape of A is
NxN where N is very large (up to 10e6, if possible...), the shape of B
and X is NxM, where M is much smaller (up to 50, say). Unfortunately, A
is symmetric but not positive-definite as the sum of each row and column
is null (A is a Laplacian matrix). B is also sparse (it is actually a
block of the original Laplacian matrix).
What kind of solver implemented in scipy would you recommend, so that I
can solve such a system with N as large as possible (on a bottom-end
computer with only 2Gb RAM)? Up to now I have used
scipy.sparse.linalg.spsolve, and the CSR scipy.sparse format for A and B:
this works fine for N as large as 300, but the memory requirements are
too high for greater values of N on my computer. Should I use an
iterative solver (I'm a complete newbie in linear algebra)? Also, spsolve
requires that the right hand side is a flat vector, so I have to solve as
many systems as the number of columns in B, which must be highly
inefficient... Any way I can solve the "whole" linear system? Any hints
will be much appreciated!
For those curious about the background, I'm trying to implement Grady's
random walker algorithm  to segment large 3-D X-ray tomography images.
N=l**3 is the number of voxels in the cubic image, and M is the number of
regions which I would like to segment. I don't require a very good
precision on the elements of the solution X.
Thanks in advance,
 L. Grady, "Random walks for image segmentation", IEEE Trans. on
pattern analysis and machine intelligence, Vol. 28, p. 1768, 2006.
More information about the SciPy-User