<< (newer post) (older post) >>     

連立一次方程式や固有値問題の解法(1)


Julia の計算をもっと速く?

もちろんコードによるんだけど、数値計算の世界ではプログラムの計算時間の大半(90% を越えることも…)が線形計算にとられることは珍しくない. だから、Julia でプログラムするとき、線形計算が速いかどうかは結構大事な話.

Julia が内部でどんな解法を使っているか

というわけで、まずは少しずつ、Julia の線形計算について調べようと思うんだな.

手始めに、特に指定しないで普通に線形計算を julia にさせたとき、どんなアルゴリズム/ライブラリが使われるか、を見てみたい. 基本的に Julia は線形計算については BLAS と LAPACK を使うことになっているそうなので、そのあたりとのからみを見る方法がわかれば十分、かな.

そして、こういうときはマクロ @code_lowered を使うと,実際にどういうルーチンが呼ばれるかわかりやすい. 例えば、

N = 10
A = randn(N,N)
b = randn(N)

@code_lowered A \ b

なんてやると,次のような出力が出る.

1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:A,:B], Any[Any[Any[:A,:Any,0],Any[:B,:Any,0],Any[:m,:Any,18],Any[:n,:Any,18],Any[symbol("#s347"),:Any,2]],Any[],3,Any[]], :(begin  # linalg/generic.jl, line 329:
        GenSym(0) = (Base.LinAlg.size)(A)
        #s347 = (top(start))(GenSym(0))
        GenSym(1) = (top(indexed_next))(GenSym(0),1,#s347)
        m = (top(getfield))(GenSym(1),1)
        #s347 = (top(getfield))(GenSym(1),2)
        GenSym(2) = (top(indexed_next))(GenSym(0),2,#s347)
        n = (top(getfield))(GenSym(2),1)
        #s347 = (top(getfield))(GenSym(2),2) # linalg/generic.jl, line 330:
        unless m == n goto 3 # linalg/generic.jl, line 331:
        unless (Base.LinAlg.istril)(A) goto 1 # linalg/generic.jl, line 332:
        unless (Base.LinAlg.istriu)(A) goto 0 # linalg/generic.jl, line 333:
        return (Base.LinAlg.Diagonal)(A) \ B
        goto 1
        0:  # linalg/generic.jl, line 335:
        return (Base.LinAlg.LowerTriangular)(A) \ B
        1:  # linalg/generic.jl, line 338:
        unless (Base.LinAlg.istriu)(A) goto 2 # linalg/generic.jl, line 339:
        return (Base.LinAlg.UpperTriangular)(A) \ B
        2:  # linalg/generic.jl, line 341:
        return (Base.LinAlg.lufact)(A) \ B
        3:  # linalg/generic.jl, line 343:
        return (Base.LinAlg.qrfact)(A,(top(apply_type))(Base.LinAlg.Val,true)) \ B
    end))))

これを見ると、julia では、いろいろ場合分けをした挙句、一般の場合は、連立一次方程式を Base.LinAlg.lufact で、つまり、LU 分解によって計算しているらしいことがわかる. より詳しく知りたければ、

?lufact

と打ち込めば、

という感じで解説が出てくる(jupyter を使っている場合). そしてさらに、

methods(lufact)

とすると、multi dispatch による多重定義の source へのリンクをすべて出してくれるので,おおよそのところがわかる. この場合、たどってみると、さらに場合分けがあって、Julia 自前の LU分解関数 generic_lufact! が呼び出されるか、もしくは LAPACK の getrf が呼び出されるようで、まあいずれにせよ、LU 分解に基づくことになる. 結局、デフォルトでは、連立一次方程式の解を求めるための計算速度は LAPACK と BLAS によってほぼ決定されることになる.

同様に、 固有値問題, 例えば eigvals についても調べてみると、これについてはあまり凝ったことはしていなくて、ほぼ直接 LAPACK の geevx! を呼び出しているだけだ.


LAPACK を直接呼び出せる?

さらに蛇足を書いておくと、LAPACK の多くの命令を Julia から直接呼び出せる 1

例えば、julia の連立一次方程式解法よりも LAPACK の gesvx! の方が,条件数やエラー評価なんかを返してくれるので情報が多い. 同じ LU 分解なので計算結果は本質的に同じだが、gesvx! の方が計算しているものが多い分,ちょっとだけ遅いようだけど.

具体的には、julia から

LAPACK.gesvx!(A,b)

とすれば直接 LAPACK の gesvx! を使える.簡単だなあ、おい.

返ってくる数字の見方は,

?LAPACK.gesvx!

とすればマニュアルが出るので,その後半部分を読めばいい.まあ,あんまり使うことはなさそうだが.


近代的な解法は?

現状では、少なくとも CG 法系統の近代的な解法は Julia にデフォルトでは用意されていない. パッケージは複数あって、例えば KrylovMethods パッケージbicgstb 関数や cgls 関数などがあるし、 IDRsSolver パッケージ には IDR ソルバー idrs が用意されている.

ちょっとだけ bicgstb を試してみよう. 入力行列は疎行列型でないとダメなので,元行列が密行列型ならば

using KrylovMethods

N = 10
A = randn(N,N)
b = randn(N)

bicgstb(sparse(A),b)

などいう感じで入力することになる.

うん、ごく普通の小さな密行列を変換して入れると,よく知られたとおり「前処理をしない場合の」bicgstb はあまりよくない. 収束しないこともよくあるなあ. 前処理なしの乱暴な計算の場合は、同じパッケージ内の cgls の方がずっと安定しそうだ 2 . 実際ためしてみると、

using KrylovMethods

N = 10
A = randn(N,N)
b = randn(N)

x = cgls(sparse(A),b)[1]

A * x - b

として、誤差を出力させてみると、

10-element Array{Float64,1}:
 -1.72892e-9
 -1.31056e-8
  3.0163e-9 
  7.52434e-9
  7.79114e-9
  2.21865e-9
 -4.52298e-9
 -3.2723e-9 
 -6.0486e-9 
  1.06731e-8

という感じの結果が得られる.うん、悪くないね.

同様に IDRsSolver パッケージの idrs 関数を試してみよう.

using IDRsSolver

N = 10
A = randn(N,N)
b = randn(N)

x = idrs(sparse(A),b)[1]

A * x - b

などとして誤差を出力してみると、

10-element Array{Float64,1}:
 -1.15118e-11
  5.93725e-12
  9.61897e-12
  5.97555e-12
 -3.91798e-12
 -1.44605e-11
  3.85769e-12
  1.02041e-11
 -5.36537e-12
  1.62759e-12

という感じで、こちらも前処理無しでなかなか悪くなさそうだ.

例えばこの idrs 関数でサイズを変えて試してみると、 100x100 のサイズで既に LU分解よりずっと速いように見える. もちろん,フル行列状態で大きくしていくと(疎行列のアクセスの問題で)かえって遅くなるので,きちんと「本当に」疎な行列で試してみる必要がありそうだ. というわけで、いずれ、大規模で疎な行列に対してこれらを試してみることにしますか.


  1. lapack.jl をみてみるといい.LAPACK の多くの命令用の wrapper が用意されている. [return]
  2. たぶん,もう少し「形」をもった疎行列で大規模な場合は bicgstb の方が良いんだろうなあ.いずれ試してみよう. [return]
     << (newer post) (older post) >>