SciPy.sparse迭代求解器：沒有稀疏的右側支持？

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

雖然文檔要求b的類型爲{array, matrix}，但如果不存在稀疏右側的支持（例如在有限元中出現），那將是非常令人驚訝的。

那麼，我在這裏做錯了什麼或爲什麼不支持？

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 asis。 的文檔表明實施LinearOperator兩種方式，但

每當有什麼事情可以通過兩種方式來完成，有人會混淆。

反正Linop下面稀疏x工程 - 用gist.github.com/denis-bz下修補lsmr.py， 。 其他稀疏迭代求解器？不知道。

---

如果你真的想要做的是：
減少| A X - B |
並且還保留| x |小，又名正規化，在L1或L2規範

那麼你一定要看看 scikit-learn。 它的目標不同角落
速度 - 準確性 - 問題 - 人（SAPP）空間 比scipy.sparse.isolve。 我發現scikit-學習固體，愉快，相當有據可查， ，但沒有比較兩個真正的問題。

你的問題有多大，多少稀疏？你能指出網絡上的測試案例嗎？

---

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