具有 (M, N) 右手尺寸矩阵的稀疏线性方程组的迭代求解答案

【问题标题】：Iterative solving of sparse systems of linear equations with (M, N) right-hand size matrix具有 (M, N) 右手尺寸矩阵的稀疏线性方程组的迭代求解
【发布时间】：2016-02-05 22:55:46
【问题描述】：

我想求解一个稀疏线性方程组：A x = b，其中 A 是一个 (M x M) 数组，b 是一个 (M x N) 数组，x 是一个 (M x 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 的任何想法？

【问题讨论】：

虽然 A 和 b 可能是稀疏的，但通常 Ax = b 的解决方案将是密集的
你看过这些方法的代码了吗？基础数学呢？也许第 2 列的解决方案本质上独立于第 1 列的解决方案，并且需要不同数量的迭代。
我还没有进入代码。

标签： python arrays numpy scipy linear-algebra

【解决方案1】：

迭代求解器和直接求解器之间的主要区别在于，直接求解器可以通过使用因式分解（通常是 Cholesky 或 LU）更有效地求解多个右手值，而迭代求解器则不能。这意味着对于直接求解器，同时求解多个列具有计算优势。

另一方面，对于迭代求解器，在同时求解多个列时没有计算增益，这可能是cg、bicg 等的 API 本身不支持矩阵求解的原因。

因此，像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 本质上都只是这样的包装器。

【讨论】：