python - SciPy.sparse 迭代求解器 : No sparse right hand side support?

Question

似乎 scipy.sparse.linalg 的迭代求解器不支持 scipy.sparse 的稀疏矩阵数据类型作为方程系统的右侧(而直接求解器会)。考虑以下简短示例:

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

虽然形状似乎是正确的:

print(A.shape)
print(b.shape)

返回

(50, 50)
(50, 1)

将 b 定义为密集的 ndarray 但是有效

b = scipy.sparse.rand(50,1,density=0.1).todense() 

如果不支持稀疏右侧(例如在 FEM 中出现)，我会感到非常惊讶，尽管文档要求 b 为 {array> 类型, 矩阵}.

那么我在这里做错了什么或者为什么这不受支持？

Answer 1

两部分方法:

• 将你的 A 包装在 scipy.sparse 中 scipy.sparse.linalg.LinearOperator (简单)
• 修补 lsmr 求解器来使用它(进行了一些调试)。

LinearOperator 的思想简单而强大:它是一个虚拟的线性运算符， 或者实际上是两个:

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

这里 matvec 和 rmatvec 可以做任何事情(在合理范围内)。 例如，他们可以线性化可怕的非线性方程 靠近任何类型的参数 x 和 y。 不幸的是，aslinearoperator 不适用于稀疏 x。 该文档建议了两种实现 LinearOperator 的方法，但是

Whenever something can be done in two ways, someone will be confused.

无论如何，下面的 Linop 与稀疏 x 一起工作——带有补丁的 lsmr.py， 在 gist.github.com/denis-bz 下. 其他稀疏迭代求解器？不知道。

---

如果你真正想做的是:
最小化 |A x - b|
并保留 |x| L1 或 L2 范数中的小，也就是正则化

那你一定要看看 scikit-learn . 它针对不同的角落
速度 - 准确性 - 问题 - 人员 (SAPP) 空间 比 scipy.sparse.isolve 。 我发现 scikit-learn 可靠、令人愉快、文档齐全， 但还没有在实际问题上比较两者。

您的问题有多大，有多稀少？你能指出网络上的测试用例吗？

---

""" 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