スパース三角システム Au*x = b を scipy スパースで効率的に解決する方法を見つけようとしています。
たとえば、スパース上三角行列 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 の solve メソッドと比較することで実証できます。(ここでは 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 は既に上三角であるため、splu の呼び出しは実際には何も行うべきではありませんが、n が大きくなると、この呼び出しは (spsolve の使用と同様に) 非常にコストがかかりますが、解決時間は短いままです。
最初に splu を呼び出さずに superLU の三角ソルバーを使用する方法はありますか? または、これを完全に行うより良い方法はありますか？