中身の詰まった dense 行列をそのまま sparse 型にしていろいろ計算すると、添字アクセスのせいで処理が重くなるはず…なんだけど、 これを前提として先へ進む前に、まあ一応、確認しておこう． Julia での計算時間の測定の練習にもなるしな．
さて、ではまず、
N = 1000
A = randn(N,N)
b = randn(N)
sA = sparse(A)
sb = sparse(b)
という感じで舞台を用意するわな． で、連立一次方程式 $A \vecbm{x} = \vecbm{b}$ を解いて、$\vecbm{x}$ を求める問題を考える． このとき、dense 型と sparse 型でそれぞれやってみて、計算時間を比べてみようっちゅうわけだ．
まず dense 型の解法だが、こちらは簡単に
x = A \ b
と書ける(計算方法は，先に確かめたとおり，LU 分解のはずだなあ)
sparse 型についても本来は同じ様に sx = sA \ sb と書けてしかるべきなんだけど 1 ，
なぜか右辺項 b を dense にとらないと動かない．
つまり，先の dense な b を使って
sx = sA \ b
とするか，持っている情報は sb のはずだという立場にたてば
sx = sA \ full(sb)
とすることになる．
どうしても右辺を sparse な sb として対象にしたい！ というなら、ごりごりと地道に
slu = lufact(sA)
sx[slu[:q]] = slu[:U] \ (slu[:L] \ (slu[:Rs] .* sb[slu[:p]]))
などとする 2 ことになるかなあ． 添字が変換されてたりして，いろんな意味であんまりスジがよくなさそうな気配だが．
さて話を戻してだ，今回は連立一次方程式のこの二つの計算 A \ b , sA \ b あたりで 3 実際にかかる時間がどれくらい違うのか、自分で確かめてみよう，というわけだな．他のももちろん試してみてもいいが．
そして計算時間の測定だけど、何回か計算して平均したほうが良さそうだ． こういうときは、@elapsed というマクロが便利．かかった時間の数字だけを出してくれるしな．
あとは、このマクロの結果を array に入れて、それを mean 関数で平均してしまえばいいので、 例えば具体的には
A \ b; # JIT コンパイルのために、一回空回ししておく．
mean( [ @elapsed A \ b for k in 1:10] )
などとすると、dense 形式での計算時間の(10回での)平均が出てくる． そうそう、JIT コンパイルのせいで妙に一回目だけ遅かったりするので、こういう感じで計算を一回空回ししておくといいな．
で、ちなみに、この結果はおおよそ 0.03878 秒だった．
次に同様に
sA \ b; # JIT コンパイルのために、一回空回ししておく．
mean( [ @elapsed sA \ b for k in 1:10] )
として sparse 型での計算時間の平均をみてみると、こちらは 0.25084 秒だった． つまり，無駄に中身の多い sparse 型行列での計算が 6,7倍遅い ので、まあ予想通り．そういう意味でメデタシメデタシだな．
あと、蛇足までに，確かに LU 分解相当だということを計算時間からも確認して見るために，
# JIT コンパイルのために、一回空回ししておく．
slu = lufact(sA);
slu \ b;
mean( [ @elapsed (
slu = lufact(sA);
slu \ b;
)
for k in 1:10] )
とやってみると、これもほぼ同じ 0.24565 秒なので、やっぱり解法としては LU 分解なんだな，ということが納得できる．
さらに、さっきの無理やり版でも
# JIT コンパイルのために、一回空回ししておく．
slu = lufact(sA);
sx[slu[:q]] = slu[:U] \ (slu[:L] \ (slu[:Rs] .* sb[slu[:p]]));
mean( [
@elapsed (
slu = lufact(sA);
sx[slu[:q]] = slu[:U] \ (slu[:L] \ (slu[:Rs] .* sb[slu[:p]]));
) for k in 1:10] )
として測定してみると、こいつは 0.46099 秒もかかる 4 ． やっぱ，あんまり良くないなあ…
というわけで、まああったりまえだけど、sparse 型の実力を見るには、ちゃんと疎行列を用意しないといけないな． それについてはまた次回ということで．