It seems like the iterative solvers of scipy.sparse.linalg do not support the sparse matrix data types of scipy.sparse as a right hand side of the equation system (whereas the direct solver does). Consider the following short example:
import numpy as np
import scipy
import scipy.sparse
import scipy.sparse.linalg
# Construct a random sparse but regular matrix A
A = scipy.sparse.rand(50,50,density=0.1)+scipy.sparse.coo_matrix(np.diag(np.random.rand(50,)))
# Construct a random sparse right hand side
b = scipy.sparse.rand(50,1,density=0.1).tocsc()
ud = scipy.sparse.linalg.spsolve(A,b) # Works as indented
ui = scipy.sparse.linalg.bicg(A,b)[0] # ValueError: A and b have incompatible dimensions
The shapes seem to be correct though:
print(A.shape)
print(b.shape)
returns
(50, 50)
(50, 1)
Defining b as a dense ndarray however works
b = scipy.sparse.rand(50,1,density=0.1).todense()
It would be very surprising to me if there would be no support for sparse right hand sides (as arise in FEM for example) although the documentation requests b to be of type {array, matrix}.
So what am I doing wrong here or why would this be unsupported?
A two-part approach:
A in a scipy.sparse
scipy.sparse.linalg.LinearOperator
(easy)The idea of LinearOperator is easy and powerful: it's a virtual linear operator, or actually two:
Aop = Linop( A ) # see below
A.dot(x) -> Aop.matvec(x)
A.T.dot(y) -> Aop.rmatvec(y)
x = solver( Aop, b ... ) # uses matvec and rmatvec, not A.dot A.T.dot
Here matvec and rmatvec can do anything (within reason).
For example, they could linearize horrible nonlinear equations
near args x and y of any type whatever.
Unfortunately aslinearoperator doesn't work asis for sparse x.
The doc suggests two ways of implementing LinearOperator, but
Whenever something can be done in two ways, someone will be confused.
Anyway, Linop below works with sparse x -- with a patched lsmr.py,
under gist.github.com/denis-bz .
Other sparse iterative solvers ? Don't know.
If what you really want to do is:
minimize |A x - b|
and also keep |x| small, a.k.a. regularization, in the L1 or L2 norm
then you should definitely look at
scikit-learn .
It targets different corners of
Speed - Accuracy - Problems - People (SAPP) space
than scipy.sparse.isolve .
I've found scikit-learn solid, pleasant, pretty well-documented,
but haven't compared the two on real problems.
How big, how sparse are your problems ? Could you point to test cases on the web ?
""" Linop( A ): .matvec .rmatvec( sparse vecs )
http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html
http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsmr.html
"""
from __future__ import division
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import LinearOperator # $scipy/sparse/linalg/interface.py
__version__ = "2015-12-24 dec denis + safe_sparse_dot"
#...............................................................................
class Linop( LinearOperator ): # subclass ?
""" Aop = Linop( scipy sparse matrix A )
-> Aop.matvec(x) = A dot x, x ndarray or sparse
Aop.rmatvec(x) = A.T dot x
for scipy.sparse.linalg solvers like lsmr
"""
def __init__( self, A ):
self.A = A
def matvec( self, x ):
return safe_sparse_dot( self.A, x )
def rmatvec( self, y ):
return safe_sparse_dot( self.A.T, y )
# LinearOperator subclass should implement at least one of _matvec and _matmat.
def _matvec( self, b ):
raise NotImplementedError( "_matvec" )
# not _matvec only:
# $scipy/sparse/linalg/interface.py
# def matvec(self, x):
# x = np.asanyarray(x) <-- kills sparse x, should raise an error
def _rmatvec( self, b ):
raise NotImplementedError( "_rmatvec" )
@property
def shape( self ):
return self.A.shape
def safe_sparse_dot( a, b ):
""" -> a * b or np.dot(a, b) """
# from sklearn
if sparse.issparse(a) or sparse.issparse(b):
try:
return a * b
except ValueError: # dimension mismatch: print shapes
print "error: %s %s * %s %s" % (
type(a).__name__, a.shape,
type(b).__name__, b.shape )
raise
else:
return np.dot(a, b)
#...........................................................................
if __name__ == "__main__":
import sys
from lsmr import lsmr # patched $scipy/sparse/linalg/lsmr.py
np.set_printoptions( threshold=20, edgeitems=10, linewidth=100, suppress=True,
formatter = dict( float = lambda x: "%.2g" % x ))
# test sparse.rand A m n, x n 1, b m 1
m = 10
n = 100
density = .1
bdense = 0
seed = 0
damp = 1
# to change these params in sh or ipython, run this.py a=1 b=None c=[3] ...
for arg in sys.argv[1:]:
exec( arg )
np.random.seed(seed)
print "\n", 80 * "-"
paramstr = "%s m %d n %d density %g bdense %d seed %d damp %g " % (
__file__, m, n, density, bdense, seed, damp )
print paramstr
A = sparse.rand( m, n, density, format="csr", random_state=seed )
x = sparse.rand( n, 1, density, format="csr", random_state=seed )
b = sparse.rand( m, 1, density, format="csr", random_state=seed )
if bdense:
b = b.toarray().squeeze() # matrix (m,1) -> ndarray (m,)
#...........................................................................
Aop = Linop( A )
# aslinearoperator( A ): not for sparse x
# check Aop matvec rmatvec --
Ax = Aop.matvec( x )
bA = Aop.rmatvec( b )
for nm in "A Aop x b Ax bA ".split():
x = eval(nm)
print "%s: %s %s " % (nm, x.shape, type(x))
print ""
print "lsmr( Aop, b )"
#...........................................................................
xetc = lsmr( Aop, b, damp=damp, show=1 )
#...........................................................................
x, istop, niter, normr, normar, norma, conda, normx = xetc
x = getattr( x, "A", x ) .squeeze()
print "x:", x.shape, x
# print "lsmr( A, b ) -- Valueerror in $scipy/sparse/linalg/interface.py"
# xetc = lsmr( A, b, damp=damp, show=1 ) # Valueerror
safe_sparse_dot( A, b.T ) # ValueError: dimension mismatch