julia> A=[ 10.0 2.0 0.0 -2.0 0.0 -4.0 0.0 0.0;
2.0 2.0 0.0 -2.0 0.0 0.0 0.0 0.0;
0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0;
-2.0 -2.0 0.0 6.0 -2.0 -2.0 0.0 0.0;
0.0 0.0 0.0 -2.0 20.0 8.0 -8.0 -8.0;
-4.0 0.0 0.0 -2.0 8.0 8.0 -4.0 -4.0;
0.0 0.0 0.0 0.0 -8.0 -4.0 8.0 0.0;
0.0 0.0 0.0 0.0 -8.0 -4.0 0.0 8.0];
julia> Asp=sparse(A);
julia> C=cholesky(Asp)
julia> # Matrix(sparse(C.L)[C.p,:])
8×8 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.41421
0.0 0.0 0.0 -2.0 1.41421 -1.41421 1.41421 0.0
0.0 0.0 0.0 0.0 1.41421 0.0 0.0 0.0
0.0 0.0 -1.0 -1.0 -1.41421 1.41421 0.0 0.0
0.0 2.82843 0.0 0.0 0.0 0.0 0.0 0.0
-1.41421 -1.41421 8.88178e-16 2.0 0.0 0.0 0.0 0.0
2.82843 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-2.82843 -2.82843 2.0 0.0 0.0 0.0 0.0 0.0
which is not lower triangular. Then I read carefully ?cholesky which is something like
cholesky(A; shift = 0.0, check = true, perm = nothing) -> CHOLMOD.Factor
Compute the Cholesky factorization of a sparse positive definite matrix A. A
must be a SparseMatrixCSC or a Symmetric/Hermitian view of a
SparseMatrixCSC. ...
reading on, it says that what is actually computed is the Cholesky factorization LL' of a permuted matrix, PAP' . The permutation minimizes the fill-ins. Example that this is the case:
julia> L=sparse(C.L)
8×8 SparseMatrixCSC{Float64,Int64} with 19 stored entries:
[1, 1] = 2.82843
[3, 1] = -2.82843
[4, 1] = -1.41421
[2, 2] = 2.82843
...
julia> Matrix(L) # is lower triangular
8×8 Array{Float64,2}:
2.82843 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 2.82843 0.0 0.0 0.0 0.0 0.0 0.0
-2.82843 -2.82843 2.0 0.0 0.0 0.0 0.0 0.0
-1.41421 -1.41421 8.88178e-16 2.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.41421 0.0 0.0 0.0
0.0 0.0 -1.0 -1.0 -1.41421 1.41421 0.0 0.0
0.0 0.0 0.0 -2.0 1.41421 -1.41421 1.41421 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.41421
julia> # is this the cholesky decomposition of permuted A? Yes, it is
L*L' ≈ A[C.p,C.p]
true
We can explicitely recover P and construct the (not any more lower triangular) Cholesky factor of A, P'L, which multiplied by its transpose should give A:
julia> P=sparse(1:8,C.p,ones(8))
julia> Matrix(P)
8×8 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
julia> # P'L
Matrix(P'*L)
8×8 Array{Float64,2}:
0.0 0.0 0.0 -2.0 1.41421 -1.41421 1.41421 0.0
0.0 0.0 0.0 0.0 1.41421 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.41421
0.0 0.0 -1.0 -1.0 -1.41421 1.41421 0.0 0.0
-2.82843 -2.82843 2.0 0.0 0.0 0.0 0.0 0.0
-1.41421 -1.41421 8.88178e-16 2.0 0.0 0.0 0.0 0.0
0.0 2.82843 0.0 0.0 0.0 0.0 0.0 0.0
2.82843 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> #and we recover our initial A
P'*L*L'*P ≈ Asp
true
The final thing is that, contrary to this, I recover P'L using, not the permutation kept in the component p (which represents P) of the internal structure of C, but its inverted permutation (invperm) which represents P':
julia> sparse(C.L)[invperm(C.p),:]
8×8 SparseMatrixCSC{Float64,Int64} with 19 stored entries:
[5, 1] = -2.82843
[6, 1] = -1.41421
...
julia> Matrix(sparse(C.L)[invperm(C.p),:]) ≈ Matrix(P'*L)
true
