Sparse Cholesky decomposition in Julia

Julia has nice built-in sparse matrices and algorithms. I spent a whole afternoon trying to extract the sparse Cholesky factorization of a sparse matrix but it did not look what I wanted. I was misunderstanding the decomposition. The documentation is not exceedingly clear but apparently it had to be done this way:

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