let N = 1000000 a = rand(N) S = spdiagm(0 => ones(N+1), -1 => a) S = S' * S C_1 = cholesky(S) L_1 = sparse(C_1.L) C_2 = cholesky(S, perm=1:N+1) L_2 = sparse(C_2.L) @assert !(L_1 * L_1' ≈ S) @assert L_2 * L_2' ≈ S end
例:精度行列から多変量正規に従うRVを乱数生成
let μ + (sparse(cholesky(Σ_inv, perm=1:D).L)' \ rand(MvNormal(I(D)))) end