解决多个线性稀疏矩阵方程: "numpy.linalg.solve" vs. "scipy.sparse.linalg.spsolve"

Question

解决多个线性稀疏矩阵方程: "numpy.linalg.solve" vs. "scipy.sparse.linalg.spsolve"

pythonnumpyscipysparse-matrixlinear-algebra

6

我需要解决大量形如“Ax=B”的线性矩阵方程，其中A是一个稀疏矩阵，主对角线上有值，而B是一个向量。

我的第一种方法是使用numpy中的密集数组和numpy.linalg.solve函数。使用(N,n,n)维数组（其中N是线性矩阵方程的数量，n是矩阵的维度）通过for循环迭代所有方程，虽然可以正常运行，但速度相对较慢。后来我发现其实还可以直接将(N,n,n)维数组传递给numpy.linalg.solve函数，这在我读的文档中并没有提到。这样就大大提高了计算速度（详见下文）。

然而，由于我有稀疏矩阵，因此我也查看了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)

第二种方法：直接将（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)

第三种方法：使用带有for循环的scipy.sparse.linalg.spsolve：

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。

由于本例中方程的数量对我来说仍然较低，并且我经常需要执行此类操作，如果能使用scipy的稀疏矩阵解决问题而不需要for循环，我会非常高兴。有人知道如何实现吗？或者也许还有其他不同的解决问题的方法？欢迎提出建议。

- Snijderfrey

我不确定去掉循环是否会有帮助（这是缓慢的Python循环与独立性指示之间的权衡；在稀疏操作中，后者可能更重要）。没有过多考虑，感觉可以构建一个分块对角矩阵，将独立问题放入其中（当然要同步增加b）。 - sascha

Sascha，我也不确定。希望能以类似于numpy的方法实现相似的效果，这样可以节省大约40%的计算时间，而不必使用for循环。除此之外，我会看看你的建议。 - Snijderfrey

循环本身应该不是问题，但numpy/scipy函数的开销可能会成为问题。此外，运行您的代码时，我收到了一个稀疏效率警告。结果表明，目前稀疏代码花费了一半以上的时间用于转换为csc/csr格式。 - Paul Panzer

这是一个有效的测试吗？矩阵都是对角线矩阵，因此求逆只是将右侧向量按对角线元素进行逐元素除法。spsolve 可以识别这种结构并执行相应的操作。 - Lutz Lehmann

LutzL，我明白你的意思。实际上，我的真实数据看起来并不是那么干净。因此，我用一些真实数据（这些数据有点难以在这里展示）重复了这些方法，并得到了类似的结果。因此，我认为这个测试或多或少是“有效”的。 - Snijderfrey

Paul Panzer，你说得对：将format='csr'添加到所有稀疏矩阵构造中，可以将第三种方法的时间降低到17.4毫秒。 - Snijderfrey

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- sascha · Accepted Answer

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

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