Julia で sparse Cholesky 分解

Julia の場合、以下のような感じで sparse Cholesky factorization ができる。

julia> A = sprand(100,100,0.1);
julia> B = A*A' + 100*speye(100);
# ここまでで B に sparse で positive definite な行列を生成

julia> C = cholfact(B);
# これで CHOLMOD.Factor が生成される

julia> L = sparse(C);
# これで普通の sparse 行列として L を取り出せる。
# ただし、AMD による行と列の順番変更がある。
# この情報は C.Perm に入っているが、CHOLMOD の情報のままのためインデックスが 0 スタートであり、Julia の１スタートのために C.Perm + 1 など１を足す必要あり。

julia> vecnorm(L*L' - B[C.Perm+1, C.Perm+1])
1.428931331049563e-13

# これで数値誤差が評価できて、数値誤差を除けば正しく計算できていることが分かる。