求解稀疏上三角系统答案

【问题标题】：Solve sparse upper triangular system求解稀疏上三角系统
【发布时间】：2013-09-07 07:24:48
【问题描述】：

我试图弄清楚如何有效地解决稀疏三角系统，Au*x = b in scipy sparse。

例如，我们可以构造一个稀疏的上三角矩阵 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 得到问题的解决方案，但是很明显三角形结构没有被利用。这可以通过定时求解并将其与 slu 中的求解方法进行比较来证明。（这里使用 iPython 的 %time 魔法）

%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 已经是上三角函数，所以对 spl 的调用实际上不应该做任何事情，但是，随着 n 变大，这个调用变得非常昂贵（就像使用 spsolve 一样），而求解时间仍然很短.

有没有什么方法可以在不先调用 slu 的情况下使用 superLU 的三角求解器？还是有更好的方法来做到这一点？

【问题讨论】：

我在 ipython. spsolve 中使用 timeit，测量时间约为 90 毫秒。 splu 需要几秒钟。 sla.spsolve(Au,b,use_umfpack=False) 在 1-2 秒范围内。 linalg.solve_triangular(Au.toarray(),b) 比 spsolve 慢（200 毫秒）。还要比较答案。 x2 的大值与 x1 的值不接近。
您应该考虑迭代解决方案是否比直接解决方案更适合您的问题。请参阅此讨论：scicomp.stackexchange.com/questions/81/…
我需要三角求解器来实现迭代求解器的 SSOR 预处理器。我已经使用 fortran 和 f2py 编写了一个快速的解决方案，但仍然宁愿使用原生 python/主要包解决方案。

标签： python scipy linear-algebra sparse-matrix

【解决方案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

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

【讨论】：