有效求解 python 中具有非零对角的三对角系统

Question

我想求解系统 A.b=x 其中 A 是几乎 [=36= 中的一个三对角矩阵]:

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:[=14 中的 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，否则可能运行内存不足。下面是 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

加速确实很重要。

有效求解 python 中具有非零对角的三对角系统

solving tridiagonal system with non-zero opposite corners in python efficiently

python

numpy

linear-algebra

scipy