求解多个线性稀疏矩阵方程：“numpy.linalg.solve”与“scipy.sparse.linalg.spsolve”

Question

我必须求解 x 的大量“Ax=B”类型的线性矩阵方程，其中 A 是一个稀疏矩阵，主要填充主对角线，B 是一个向量。

\n\n

我的第一种方法是通过 numpy.linalg.solve 使用密集的 numpy 数组来实现此目的，并且它可以很好地处理 (N,n,n) 维数组，其中 N 是线性矩阵方程的数量，n 是方阵维度。我首先将它与迭代所有方程的 for 循环一起使用，这实际上相当慢。但后来意识到，您也可以将 (N,n,n) 维矩阵直接传递给 numpy.linalg.solve ，而无需任何 for 循环（顺便说一下，我在阅读的文档中没有找到）。这已经大大提高了计算速度（详细信息见下文）。

\n\n

但是，因为我有稀疏矩阵，所以我还查看了 scipy.sparse.linalg.spsolve 函数，它与相应的 numpy 函数执行类似的操作。使用 for 循环迭代所有方程比 numpy 解决方案快得多，但似乎不可能将 (N,n,n) 维数组直接传递给 scipy\xc2\xb4s spsolve。

\n\n

这是我到目前为止所尝试的：

\n\n

首先，我计算一些虚构的 A 矩阵和带有随机数的 B 向量以用于测试目的，包括稀疏和密集：

\n\n

import numpy as np\nfrom scipy import sparse\nfrom scipy.sparse.linalg import spsolve\n\nnumber_of_systems = 100 #corresponds to N in the text\nnumber_of_data_points = 1000 #corresponds to n in the text\n\n#calculation of sample matrices (dense and sparse)\nA_sparse = np.empty(number_of_systems,dtype=object)\nA_dense = np.empty((number_of_systems,number_of_data_points,number_of_data_points))\n\nfor ii in np.arange(number_of_systems):\n    A_sparse[ii] = sparse.spdiags(np.random.random(number_of_data_points),0,number_of_data_points,number_of_data_points)\n    A_dense[ii] = A_sparse[ii].todense()\n\n#calculation of sample vectors\nB = np.random.random((number_of_systems,number_of_data_points))\n

\n\n

1）第一种方法：带有for循环的numpy.linalg.solve：

\n\n

def solve_dense_3D(A,B):\n    results = np.empty((A.shape[0],A.shape[1]))\n    for ii in np.arange(A.shape[0]):\n        results[ii] = np.linalg.solve(A[ii],B[ii])\n    return results\n\nresult_dense_for = solve_dense_3D(A_dense,B)\n

\n\n

定时：

\n\n

timeit(solve_dense_3D(A_dense,B))\n1.25 s \xc2\xb1 27.6 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n

\n\n

2）第二种方法：将（N，n，n）维矩阵直接传递给numpy.linalg.solve：

\n\n

result_dense = np.linalg.solve(A_dense,B)\n

\n\n

定时：

\n\n

timeit(np.linalg.solve(A_dense,B))\n769 ms \xc2\xb1 9.68 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n

\n\n

3）第三种方法：使用 scipy.sparse.linalg.spsolve 和 for 循环：

\n\n

def solve_sparse_3D(A,B):\n    results = np.empty((A.shape[0],A[0].shape[0]))\n    for ii in np.arange(A.shape[0]):\n        results[ii] = spsolve(A[ii],B[ii])\n    return results\n\nresult_sparse_for = solve_sparse_3D(A_sparse,B)\n

\n\n

定时：

\n\n

timeit(solve_sparse_3D(A_sparse,B))\n30.9 ms \xc2\xb1 132 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n

\n\n

很明显，能够省略方法 1 和 2 中的 for 循环是有好处的。正如可能预期的那样，迄今为止最快的替代方案是使用稀疏矩阵的方法 3。

\n\n

因为这个例子中的方程数量对我来说仍然相当少，而且因为我必须经常做这样的事情，所以如果有一个使用 scipy\xc2\xb4s 稀疏矩阵而不使用 for 循环的解决方案，我会很高兴。有人知道实现这一目标的方法吗？或者也许有另一种方法以甚至不同的方式解决问题？我很乐意提供建议。

\n

Answer 1

一个小演示概述了我上面评论中的想法：

""" YOUR CODE (only imports changed + deterministic randomness) """

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
from scipy.sparse import block_diag
from time import perf_counter as pc

np.random.seed(0)

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

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

start = pc()
result_sparse_for = solve_sparse_3D(A_sparse,B)
end = pc()
print(result_sparse_for)
print(end - start)

""" ALTERNATIVE APPROACH """

def solve_sparse_3D_blockdiag(A,B):
    oldN = B.shape[0]

    A_ = block_diag(A)    # huge sparse block-matrix of independent problems
    B_ = B.ravel()        # flattened vector

    results = spsolve(A_, B_)
    return results.reshape(oldN, -1)    # unflatten results

start = pc()
result_sparse_for = solve_sparse_3D_blockdiag(A_sparse,B)
end = pc()
print(result_sparse_for)
print(end - start)

哪个输出

[[ 0.97529866  1.26406276  0.83348888 ...  0.99310639  3.90781207
0.16724226]
[ 1.23398934 28.82088739  1.6955886  ...  1.85011686  0.23386882
1.17208753]
[ 0.92864777  0.22248781  0.09445412 ...  2.5080376   0.91701228
0.97266564]
...
[ 0.33087093  0.89034736  1.7523883  ...  0.2171746   4.89236164
0.31546549]
[ 1.2163625   3.0100941   0.87216264 ...  1.62105596  0.33211353
2.07929302]
[ 5.35677404  1.23830776  0.16073721 ...  0.26492506  0.53676822
3.73192617]]
0.08764066299999995

###

[[ 0.97529866  1.26406276  0.83348888 ...  0.99310639  3.90781207
0.16724226]
[ 1.23398934 28.82088739  1.6955886  ...  1.85011686  0.23386882
1.17208753]
[ 0.92864777  0.22248781  0.09445412 ...  2.5080376   0.91701228
0.97266564]
...
[ 0.33087093  0.89034736  1.7523883  ...  0.2171746   4.89236164
0.31546549]
[ 1.2163625   3.0100941   0.87216264 ...  1.62105596  0.33211353
2.07929302]
[ 5.35677404  1.23830776  0.16073721 ...  0.26492506  0.53676822
3.73192617]]
0.07241856000000013

有一些事情要做：

• 使用您原来的更理智的基准测试方法
• 以正确的稀疏格式构建 block_diag，以消除一些潜在的警告和减速：请参阅文档
• 调整 spsolve 的参数permc_spec