有效地求解 Python 中具有非零对角的三对角线系统

Question

我想求解系统A.b=x，其中A几乎是python中的三对角矩阵：

A 是这样的矩阵

a b 0 0 .... 0 0 b
b a b 0 .... 0 0 0
0 b a b .... 0 0 0
.
.
0 0 0 0 .... b a b
b 0 0 0 .... 0 b a

即具有非零相对角的三对角线。

我可以使用 numpy 求解器求解和集成我的系统：

numpy.linalg.solve

这有效，但非常慢，因为我的矩阵很大，我认为它没有利用 A 数组的稀疏性和近三对角线

。如果它是一个纯三对角线系统，我知道如何使用经典的正向和后向替换算法快速有效地解决它，但我被那些非零的对角难住了。 我浏览了numpy和scipy，我唯一能想到的就是尝试将NxN矩阵转换为带状系统，并尝试使用scipy中的solve_banded：

https://docs.scipy.org/doc/scipy/reference/linalg.html

我是否错过了一些明显的东西，有没有技巧可以使用python numpy或scipy包的内置功能有效地解决这个系统？

Answer 1

这是循环系统，可以用 O(N log N 中的 FFT 求解(。请参阅scipy.linalg.solve_circulant。

我不知道 massive 是什么意思，但我想它大约是 100000，否则可能会用完 RAM。下面是 N=10000 的稍小情况的代码。

import scipy.linalg
import numpy as np
from time import time
N = 10000
a, b = 1, 2
y = np.random.uniform(size=N)
# make big matrix
M = np.zeros((N,N))
np.fill_diagonal(M, a)
np.fill_diagonal(M[1:,:], b)
np.fill_diagonal(M[:,1:], b)
M[-1, 0] = M[0, -1] = b
tic = time()
x0 = np.linalg.solve(M, y)
toc = time()
print("np.linalg.solve", toc - tic)
tic = time()
# just use first row
x1 = scipy.linalg.solve_circulant(M[0], y)
toc = time()
print("scipy.linalg.solve_circulant", toc - tic)
print(np.isclose(x0, x1).all())

结果是：

np.linalg.solve 7.422604322433472
scipy.linalg.solve_circulant 0.0010323524475097656
True

加速确实很重要。