【发布时间】：2019-12-20 10:22:37
【问题描述】：
我必须为 x 求解大量“Ax=B”类型的线性矩阵方程，其中 A 是一个稀疏矩阵，主要填充了主对角线，B 是一个向量。
我的第一种方法是为此目的使用密集的 numpy 数组和 numpy.linalg.solve，它适用于 (N,n,n) 维数组，其中 N 是线性矩阵方程的数量，n 是方阵维数。我首先将它与遍历所有方程的 for 循环一起使用，这实际上相当慢。但后来意识到您也可以将 (N,n,n) 维矩阵直接传递给 numpy.linalg.solve 而无需任何 for 循环（顺便说一下，我在阅读的文档中没有找到）。这已经大大提高了计算速度（详情见下文）。
但是，因为我有稀疏矩阵，所以我还查看了 scipy.sparse.linalg.spsolve 函数，它执行类似的操作，例如相应的 numpy 函数。使用 for 循环遍历所有方程比 numpy 解决方案快得多，但似乎不可能将 (N,n,n) 维数组直接传递给 scipy 的 spsolve。
到目前为止，这是我尝试过的：
首先，我用随机数计算了一些虚构的 A 矩阵和 B 向量，用于测试目的，包括稀疏和密集：
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
number_of_systems = 100 #corresponds to N in the text
number_of_data_points = 1000 #corresponds to n in the text
#calculation of sample matrices (dense and sparse)
A_sparse = np.empty(number_of_systems,dtype=object)
A_dense = np.empty((number_of_systems,number_of_data_points,number_of_data_points))
for ii in np.arange(number_of_systems):
A_sparse[ii] = sparse.spdiags(np.random.random(number_of_data_points),0,number_of_data_points,number_of_data_points)
A_dense[ii] = A_sparse[ii].todense()
#calculation of sample vectors
B = np.random.random((number_of_systems,number_of_data_points))
1) 第一种方法：numpy.linalg.solve 与 for 循环：
def solve_dense_3D(A,B):
results = np.empty((A.shape[0],A.shape[1]))
for ii in np.arange(A.shape[0]):
results[ii] = np.linalg.solve(A[ii],B[ii])
return results
result_dense_for = solve_dense_3D(A_dense,B)
时间：
timeit(solve_dense_3D(A_dense,B))
1.25 s ± 27.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2）第二种方法：将（N，n，n）维矩阵直接传递给numpy.linalg.solve：
result_dense = np.linalg.solve(A_dense,B)
时间：
timeit(np.linalg.solve(A_dense,B))
769 ms ± 9.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
3) 第三种方法：将 scipy.sparse.linalg.spsolve 与 for 循环一起使用：
def solve_sparse_3D(A,B):
results = np.empty((A.shape[0],A[0].shape[0]))
for ii in np.arange(A.shape[0]):
results[ii] = spsolve(A[ii],B[ii])
return results
result_sparse_for = solve_sparse_3D(A_sparse,B)
时间：
timeit(solve_sparse_3D(A_sparse,B))
30.9 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
很明显，能够从方法 1 和 2 中省略 for 循环是一个优势。到目前为止，最快的替代方案是，正如可能预期的那样，使用稀疏矩阵的方法 3。
因为这个例子中的方程数量对我来说仍然相当少，而且我必须经常做这样的事情，如果有一个没有 for 循环的使用 scipy 的稀疏矩阵的解决方案，我会很高兴。有人知道实现这一目标的方法吗？或者也许还有另一种方法可以以不同的方式解决问题？我很乐意提出建议。
【问题讨论】：
-
我不太确定丢失循环是否会有所帮助（它是缓慢的 python 循环与独立性的指示；后者在稀疏操作中可能更重要）。不用想太多，感觉就像你可以为你的独立问题构建一个block-diag matrix（当然，b 同步增长）。
-
萨沙，我也不确定。我希望能够以某种方式实现类似的效果，就像使用 numpy 方法一样，在不使用 for 循环时可以节省大约 40% 的计算时间。除此之外，我会看看你的建议。
-
循环本身在这里不应该是一个问题，尽管 numpy/scipy 函数开销可能是。此外，运行您的代码我得到一个稀疏的效率警告。事实证明，稀疏代码花费了一半以上的时间转换为 csc/csr 格式。
-
这是一个有效的测试吗？矩阵都是对角的，因此求逆只是右侧向量除以对角线的元素。
spsolve可以识别这种结构并执行那么多操作。 -
LutzL，我明白你的意思。事实上，我的真实数据看起来并不那么干净。因此，我用一些真实数据（这里有点难以展示）重复了这些方法，结果质量相似。所以我会假设测试或多或少是“有效的”。
标签： python numpy scipy sparse-matrix linear-algebra