在python中有效地解决具有非零对角的三对角系统答案

【问题标题】：solving tridiagonal system with non-zero opposite corners in python efficiently在python中有效地解决具有非零对角的三对角系统
【发布时间】：2020-06-10 17:59:25
【问题描述】：

我想解决系统 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 包的内置功能有效地解决这个系统？

【问题讨论】：

标签： python numpy scipy linear-algebra

【解决方案1】：

这是一个循环系统，可以用 O(N log N) 的 FFT 求解。见scipy.linalg.solve_circulant。

我不知道大规模是什么意思，但我猜它大约是 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

加速确实很重要。

【讨论】：