加速Python中的稀疏求逆、逐元素乘法和加法

Question

我正在尝试反转三对角的 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秒。

Answer 1

首先，从数值稳定性 POV 来看，求解线性系统几乎总是优于直接反演。

其次，如果您知道矩阵是三对角的，您可能需要查看 scipy.linalg 中的

_banded

和

_tridiagiaginal

函数系列。甚至来自

scipy.linalg.lapack

的裸 lapack 包装，例如 https://docs.scipy.org/doc/scipy/reference/ generated/scipy.linalg.lapack.dgtsv.html