求解稀疏上三角系统

Question

我正试图找出如何有效地求解一个稀疏三角形系统，在scipy稀疏中Au*x=b。

例如，我们可以构造一个稀疏的上三角矩阵Au和一个右手边b，其中：

import scipy.sparse as sp
import scipy.sparse.linalg as sla
import numpy as np
n = 2000
A = sp.rand(n, n, density=0.4) + sp.eye(n)
Au = sp.triu(A).tocsr()
b = np.random.normal(size=(n))

我们可以使用spsolve来解决这个问题，但很明显，三角形结构没有被利用。这可以通过对解决方案进行计时并将其与splu中的解决方法进行比较来证明。（此处使用iPython的%时间魔术）

%time x1 = sla.spsolve(Au,b)
CPU times: user 3.63 s, sys: 79.1 ms, total: 3.71 s
Wall time: 1.1 s
%time Au_lu = sla.splu(Au)
CPU times: user 3.61 s, sys: 62.2 ms, total: 3.67 s
Wall time: 1.08 s
%time x2 = Au_lu.solve(b)
CPU times: user 25 ms, sys: 332 µs, total: 25.4 ms
Wall time: 7.01 ms

由于Au已经是上三角的，调用splu实际上不应该做太多事情，然而，随着n变大，这个调用变得昂贵得令人望而却步（使用spsolve也是如此），而求解时间仍然很短。

有没有任何方法可以在不首先调用splu的情况下使用superLU的三角形解算器？或者有更好的方法来做到这一点吗？

Answer 1

恐怕这不是很有指导意义，但您尝试过更改列排列吗？当我使用"自然"时，我会得到巨大的加速。

%time x1 = sla.spsolve(Au, b, permc_spec="NATURAL")
CPU time: user 46.7 ms, sys: 0 ns, total: 46.7 ms
Wall time: 49 ms

对我来说，它没有使用splu函数输出的求解方法那么快，但它似乎相当接近（并避免调用splu）。也许这就足够了？Scipy Docs