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

这是可行的，但速度非常慢，因为我的矩阵非常大，而且我认为它并没有利用数组的稀疏性和接近三对角性

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

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

我是否遗漏了一些明显的东西，是否有一个技巧可以使用python numpy或scipy包的内置函数高效地解决这个系统

Answer 1

这是一个循环系统，可以用O（N logn）中的FFT来求解。见scipy.linalg.solve_circulant

我不知道“海量”意味着什么，但我猜它大约是100000，否则很可能会耗尽内存。下面是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

加速确实很重要