[SciPy-User] solving large sparse linear system with Laplacian matrix

David Warde-Farley dwf@cs.toronto....
Fri Oct 30 04:28:13 CDT 2009

On 30-Oct-09, at 3:46 AM, Emmanuelle Gouillart wrote:

> 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).

I'm not sure they exist in scipy.sparse.linalg, but looking into your  
problem a little further it appears that several options exist in  
LAPACK (but are not wrapped in SciPy, yet). Conveniently, they have  
built-in support for multiple right-hand-sides.

Solve a general banded system of linear equations:

Convert a symmetric banded matrix to symmetric tridiagonal form:

The tridiagonal form is useful because you can then use this  
specialized routine for solving symmetric tridiagonal systems:

So, if both memory and performance are at issue, I'd recommend looking  
into wrapping dsbtrd and dptsv with f2py (these are for doubles, hence  
the 'd' - if you're interested in single precision use the 's'  
versions). Note that they require you pass in the band matrix in a  
particular compressed format, which is kid of a pain, but on the other  
hand _will_ save you space over CSR or anything like that. The layout  
you need to use for symmetric band matrices is near the bottom of this  


Hope this helps,


More information about the SciPy-User mailing list