具有(M,N)右手大小矩阵的稀疏线性方程组的迭代求解

Question

我想求解稀疏线性方程组:A x = b,其中A是(M×M)阵列,b是(M×N)阵列,x是和(M×N)阵列.

我用以下三种方式解决了这个问题:

• scipy.linalg.solve(A.toarray(), b.toarray()),
• scipy.sparse.linalg.spsolve(A, b),
• scipy.sparse.linalg.splu(A).solve(b.toarray()) #返回密集数组

我希望使用迭代scipy.sparse.linalg方法解决问题:

• scipy.sparse.linalg.cg,
• scipy.sparse.linalg.bicg,
• ...

然而,方法仅支持具有形状(M,)或(M,1)的右手侧b.关于如何将这些方法扩展到(M x N)数组b的任何想法？

Answer 1

迭代求解器和直接求解器之间的关键区别在于,直接求解器可以通过使用分解(通常是Cholesky或LU)来更有效地求解多个右手值,而迭代求解器则不能.这意味着对于直接求解器,同时求解多个列具有计算优势.

对于迭代求解,在另一方面,有一个在同时解决多个字段没有计算增益,这可能是为什么矩阵解决方案并不API中本地支持cg,bicg等等.

因此,直接解决方案scipy.sparse.linalg.spsolve可能是您的最佳选择.如果由于某种原因你仍然需要迭代解决方案,我只需创建一个简单的便利函数,如下所示:

from scipy.sparse.linalg import bicg

def bicg_solve(M, B):
    X, info = zip(*(bicg(M, b) for b in B.T))
    return np.transpose(X), info

然后,您可以创建一些数据并按如下方式调用它:

import numpy as np
from scipy.sparse import csc_matrix

# create some matrices
M = csc_matrix(np.random.rand(5, 5))
B = np.random.rand(5, 4)

X, info = bicg_solve(M, B)
print(X.shape)
# (5, 4)

任何在右侧接受矩阵的迭代求解器API基本上只是这样的包装器.