我正在尝试反转三对角的 10000x10000 稀疏矩阵,但我遇到的问题是
scipy.sparse.linalg.inv()
太慢,所以我尝试使用 spsolve()
代替。主要问题是,每当我与 *
执行逐元素乘法时,稀疏矩阵都会改变形状。我原来的代码是:
mport numpy as np
from scipy import sparse
def generate_tri_diagonal_matrix(size):
main_diagonal = np.random.rand(size)
upper_diagonal = np.random.rand(size - 1)
lower_diagonal = np.random.rand(size - 1)
tri_diagonal_matrix = np.diag(main_diagonal) + np.diag(upper_diagonal, k=1) + np.diag(lower_diagonal, k=-1)
return sparse.csc_matrix(tri_diagonal_matrix)
M=10000
matrix_size = M
c = 1.53 # Some random float
A = generate_tri_diagonal_matrix(matrix_size)
A_0 = A*2
A_1 = A*3
A_2 = A*4
P = np.random.rand(M)
I = np.eye(M)
inv_ini_1 = sparse.csc_matrix(I - c*A_1)
inv_1 = sparse.linalg.inv(inv_ini_1)
inv_ini_2 = sparse.csc_matrix(I - c*A_2)
inv_2 = sparse.linalg.inv(inv_ini_2)
Y_0 = P + delta_t*A*P
Y_1 = inv_1 * (Y_0 - c*A_1*P)
Y_2 = inv_2 * (Y_1 - c*A_2*P)
P = Y_2
但是当我尝试使用时:
# Attempt
inv_1 = sparse.csc_matrix(I - c*A_1)
inv_2 = sparse.csc_matrix(I - c*A_2)
Y_0 = sparse.csc_matrix(P + delta_t*A*P)
Y_1 = sparse.linalg.spsolve(inv_1 , Y_0 - c*A_1*P)
Y_2 = sparse.linalg.spsolve(inv_2 , Y_1 - c*A_2*P)
P = Y_2
Y_0 - c*A_1*P
中的Y_1
从(10k,10k)矩阵变为(10k,)矩阵。为什么会发生这种情况以及如何加速这种反转?因为每个linalg.inv()
需要我的电脑15秒。
首先,从数值稳定性 POV 来看,求解线性系统几乎总是优于直接反演。
其次,如果您知道矩阵是三对角的,您可能需要查看 scipy.linalg 中的
_banded
和 _tridiagiaginal
函数系列。甚至来自 scipy.linalg.lapack
的裸 lapack 包装,例如 https://docs.scipy.org/doc/scipy/reference/ generated/scipy.linalg.lapack.dgtsv.html