Search code examples
pythonnumpyscipysparse-matrix

How to use a sparse matrix in numpy.linalg.solve

I want to solve the following linear system for x 

Ax = b

Where A is sparse and b is just regular column matrix. However when I plug into the usual np.linalg.solve(A,b) routine it gives me an error. However when I do np.linalg.solve(A.todense(),b) it works fine.

Question.

How can I use this linear solve still preserving the sparseness of A?. The reason is A is quite large about 150 x 150 and there are about 50 such matrices and so keeping it sparse for as long as possible is the way I'd prefer it.

I hope my question makes sense. How should I go about achieving this?

Solution

  • np.linalg.solve only works for array-like objects. For example it would work on a np.ndarray or np.matrix (Example from the numpy documentation):

    import numpy as np

a = np.array([[3,1], [1,2]])
b = np.array([9,8])
x = np.linalg.solve(a, b)

    or

    import numpy as np

a = np.matrix([[3,1], [1,2]])
b = np.array([9,8])
x = np.linalg.solve(a, b)

    or on A.todense() where A=scipy.sparse.csr_matrix(np.matrix([[3,1], [1,2]])) as this returns a np.matrix object.

    To work with a sparse matrix, you have to use scipy.sparse.linalg.spsolve (as already pointed out by rakesh)

    import numpy as np
import scipy.sparse
import scipy.sparse.linalg

a = scipy.sparse.csr_matrix(np.matrix([[3,1], [1,2]]))
b = np.array([9,8])
x = scipy.sparse.linalg.spsolve(a, b)

    Note that x is still a np.ndarray and not a sparse matrix. A sparse matrix will only be returned if you solve Ax=b, with b being a matrix and not a vector.