python - 求解多个线性稀疏矩阵方程:“numpy.linalg.solve”与“scipy.sparse.linalg.spsolve”
问题描述
我必须为 x 求解大量“Ax = B”类型的线性矩阵方程,其中 A 是一个稀疏矩阵,主要填充了主对角线,B 是一个向量。
我的第一种方法是通过 numpy.linalg.solve 使用密集的 numpy 数组来实现此目的,并且它适用于 (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 的稀疏矩阵的解决方案,我会很高兴。有人知道实现这一目标的方法吗?或者也许还有另一种方法可以以不同的方式解决问题?我很乐意提出建议。
解决方案
一个小演示从我上面的评论中概述了这个想法:
""" 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
推荐阅读
- flutter - 通过 kafka mqtt-proxy 发送 Protobuf 消息
- go - golang如何自动将logrus错误写入supervisord stderr_logfile
- blockchain - 如何检查以太坊地址是否有效?
- python - SQLALCHAM 关系
- java - 使用 JPA 的 varchar 中的内容长度
- aggregation-framework - 在 $slice 中使用表达式
- kotlin - (保存设置失败,请重启IntelliJ IDEA)下载intellij idea 2020.2并导入2020.1制作的项目后报错
- c# - SMTP 需要安全连接网络邮件
- swift - 如何获取管理 API 访问令牌以在 Auth0 /api/v2/users 端点上使用
- replace - Notepad++ - 增量替换所有数组索引