解多个线性稀疏矩阵方程:对比"scipy.sparse.linalg.spsolve" [英] Solving multiple linear sparse matrix equations: "numpy.linalg.solve" vs. "scipy.sparse.linalg.spsolve"
问题描述
我必须求解大量类型为x的"Ax = B"的线性矩阵方程，其中A是主要填充对角线的稀疏矩阵，B是向量.
I have to solve a large amount of linear matrix equations of the type "Ax=B" for x where A is a sparse matrix with mainly the main diagonal populated and B is a vector.
我的第一种方法是为此使用numpy.linalg.solve密集的numpy数组，并且它在(N，n，n)维数组中也可以正常工作，其中N是线性矩阵方程的个数，而n是线性矩阵方程的个数.方阵尺寸.我首先将其与for循环一起使用，遍历所有方程式，这实际上很慢.但是，然后意识到您还可以将(N，n，n)维矩阵直接传递给numpy.linalg.solve，而无需任何for循环(通过我在阅读的文档中没有找到的方式).这已经大大提高了计算速度(详细信息请参见下文).
My first approach was to use dense numpy arrays for this purpose with numpy.linalg.solve, and it works fine with a (N,n,n)-dimensional array with N being the number of linear matrix equations and n the square matrix dimension. I first used it with a for loop iterating through all equations, which in fact is rather slow. But then realized that you can also pass the (N,n,n)-dimensional matrix directly to numpy.linalg.solve without any for loop (which by the way I did not find in the documentation I read). This already gave a good increase in computation speed (details see below).
但是，由于矩阵稀疏，因此我还查看了scipy.sparse.linalg.spsolve函数，该函数执行类似相应的numpy函数的操作.使用for循环遍历所有方程式比numpy解决方案快得多，但似乎不可能直接将(N，n，n)维数组直接传递给scipy的散点.
However, because I have sparse matrices, I also had a look at the scipy.sparse.linalg.spsolve function which does similar things like the corresponding numpy function. Using a for loop iterating through all equations is much, much faster than the numpy solution, but it seems impossible to pass the (N,n,n)-dimesional array directly to scipy´s spsolve.
这是我到目前为止尝试过的:
Here is what I tried so far:
首先，我计算一些虚构的A矩阵和B向量，这些数字具有随机数以用于测试目的(稀疏和稠密):
First, I calculate some fictional A matrices and B vectors with random numbers for test purposes, both sparse and dense:
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循环:
1) First approach: numpy.linalg.solve with for loop:
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:
2) Second approach: Passing the (N,n,n)-dimensional matrix directly to 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循环配合使用:
3) Third approach: Using scipy.sparse.linalg.spsolve with a for loop:
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是可以预期的.
It is obvoius that there is an advantage being able to omit the for loop from approach 1 and 2. By far the fastest alternative is, as could probably be expected, approach 3 with sparse matrices.
因为在这个例子中，方程式的数量对我来说仍然很低，并且因为我必须经常做类似的事情，所以如果有使用scipy的稀疏矩阵而不带for循环的解决方案，我会很高兴.有人知道实现这一目标的方法吗?或者，也许还有另一种以不同的方式解决问题的方法?我很乐意提供建议.
Because the number of equations in this example is still rather low for me and because I have to do things like that very often, I would be happy if there was a solution using scipy´s sparse matrices without a for loop. Is anybody aware of a way to achieve that? Or maybe there is another way to solve the problem in an even different way? I would be happy for suggestions.
推荐答案
通过上面的评论概述了这个想法的小演示:
A small demo outlining the idea from my comment above:
""" 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
这篇关于解多个线性稀疏矩阵方程:对比"scipy.sparse.linalg.spsolve"的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
