解决稀疏上三角系统 [英] Solve sparse upper triangular system

问题描述

我正在尝试找出如何有效解决稀疏三角系统Au * x = b的问题.

I'm trying to figure out how to efficiently solve a sparse triangular system, Au*x = b in scipy sparse.

例如，我们可以使用以下公式构造稀疏的上三角矩阵Au和右手边b:

For example, we can construct a sparse upper triangular matrix, Au, and a right hand side b with:

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魔术)

We can get a solution to the problem using spsolve, however it is clear that the triangular structure is not being taken advantage of. This can be demonstrated by timing the solution and comparing it to the solve method in splu. (Here using iPython's %time magic)

%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一样)，而求解时间仍然很小

As Au is already upper-triangular the call to splu really shouldn't be doing much of anything, however, as n gets large this call becomes prohibitively expensive (as does the use of spsolve), while the solve time remains small.

是否可以在不先调用splu的情况下使用superLU的三角求解器?还是有更好的方法完全做到这一点?

Is there any way of using superLU's triangular solver without first calling splu? Or is there a better way of doing this altogether?

推荐答案

恐怕这不是很有启发性的，但是您是否尝试过更改列排列?当我使用"NATURAL"时，我得到了极大的提速.

I'm afraid this isn't terribly instructive, but have you tried changing the column permutation? When I use "NATURAL", I get huge speedups.

%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函数输出的Solve方法快，但它似乎已经相当接近(并避免了调用splu).也许就足够了? Scipy Docs

For me, it's not quite as fast as using the solve methods of the splu function output, but it seems to get reasonably close (and avoids calling splu). Perhaps that will suffice? Scipy Docs

这篇关于解决稀疏上三角系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！