2013-09-03 5 views
5

私は疎な三角系をどのように効率よく解くかを考えています。スパースな上三角系を解く

例えば、我々がスパース上三角行列、Au及び右辺bを構築することができます。しかし、それがあることは明らかである、我々はSPSOLVEを使用して問題の解決策を得ることができます

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

は、三角形の構造は利用されていません。これは、解をタイミング化し、spluの解法と比較することで実証できます。同じように(nが大きくなる。この呼び出しは法外に高価になるとAuなどの

%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 

は、すでに本当にspluへの呼び出しがしかし、何の多くを行うべきではない上三角である(ここでiPythonの%の時間の魔法を使用して) spsolveの使用)、解決時間は小さいままです。

最初にspluを呼び出すことなく、superLUの三角ソルバーを使用する方法はありますか?あるいは、これを全部行うより良い方法がありますか?

+1

私は 'ipythonで' timeitを使います。 'spsolve'は約90msです。 'splu'には数秒かかります。 'sla.spsolve(Au、b、use_umfpack = False)'は1-2秒の範囲です。 'linalg.solve_triangular(Au.toarray()、b)'は 'spsolve'(200ms)よりも遅いです。また、答えを比較してください。 x2の大きな値はx1の値に近くない。 – hpaulj

+0

反復的な解決策が、直接的な解決策の問題に適しているかどうかを検討する必要があります。このディスカッションを参照してください:http://scicomp.stackexchange.com/questions/81/what-guidelines-should-i-follow-when-choosing-a-sparse-linear-system-solver –

+1

三角ソルバーが必要です反復ソルバーのためのSSOR前提条件を実装します。私はfortranとf2pyを使って高速な解決法を書いたが、むしろネイティブのpython/majorパッケージの解決法である。 – dwfm

答えて

0

これはひどく有益ではありませんが、列の並べ替えを変更しようとしましたか?私が "NATURAL"を使うと、私は非常に高速化します。私にとって

%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