Julia: No way to extract sparse Cholesky factorization

[aterenin]

I'm not sure if this is a problem with the code or if I'm just not able to figure out how the API works from the documentation - but I cannot seem to obtain a sparse Cholesky factorization of a matrix. I am not using Cholesky to solve a linear system and need the factorization itself to use in simulating from a multivariate Gaussian distribution.

A = sparse(Float64[ 10 1 1 1; 1 10 0 0; 1 0 10 0; 1 0 0 10])
L_dense = sparse(chol(Matrix(A)))
L_sparse = sparse(cholfact(A))

Expected: L_dense == L_sparse at least approximately
Result: L_dense is correct, but L_sparse == A

In the event I am not using the API correctly, I also tried the following.

A = sparse(Float64[ 10 1 1 1; 1 10 0 0; 1 0 10 0; 1 0 0 10])
L_dense = sparse(chol(Matrix(A)))
L_sparse = sparse(cholfact(A)[:PtL])

This fails with the following error: indexing not defined for Base.SparseArrays.CHOLMOD.FactorComponent{Float64,:PtL}. Note that I really do need [:PtL] and that [:L], which works correctly, will not suffice.

[andreasnoack]

We currently don't support conversion of CHOLMOD.FactorComponent{Float64,:PtL} to SparseMatrixCSC which would be the right thing to try here since CHOLMOD.Factor is a factorized representation of the original input matrix. Hence sparse(cholfact(A)) gives something close to A. The function chol produces the triangular Cholesky factor but isn't supported in the sparse case. The many uses of "factor" here are indeed a bit confusing but we are following the naming scheme of CHOLMOD which is the C library that computes the factorization. You can try to take a look at the header files for CHOLMOD and see how the permuted triangular Cholesky factor can be extracted from a CHOLMOD.Factor.

[andreasnoack]

I realized that you can get the factor with

Lp = sparse(F[:L])[F[:p],:]

but there might be a more efficient solution.

[aterenin]

That actually gives the upper triangular part - but suffices for my purposes. Thanks!

[ViralBShah]

The syntax now (post 1.0) is: Lp = sparse(C.L)[C.p, :]

[alegarra]

The solution does not yield the same values, and they're not in the right positions, it is not triangular any more. It seems to give a decomposition of a permuted A. I am using Julia 1.1. Did I misunderstood something?

A = sparse(Float64[ 10 1 1 1; 1 10 0 0; 1 0 10 0; 1 0 0 10])
#https://github.com/JuliaLang/julia/issues/29607
aa=convert(Array,cholesky(Matrix(A)).L) #lower !!
aa
4×4 Array{Float64,2}:
 3.16228    0.0         0.0        0.0   
 0.316228   3.14643     0.0        0.0   
 0.316228  -0.0317821   3.14627    0.0   
 0.316228  -0.0317821  -0.0321048  3.1461
cc=cholesky(A)
bb=sparse(cc.L)[cc.p,:]
Matrix(bb)
4×4 Array{Float64,2}:
 0.316228  0.316228  0.316228  3.11448
 0.0       0.0       3.16228   0.0    
 0.0       3.16228   0.0       0.0    
 3.16228   0.0       0.0       0.0 

[masuday]

@alegarra In your example, cholesky(Matrix(A)) is not permuted, but cc=cholesky(A) is permuted. bb is permuted back to the original order. It is no longer lower-triangular, but you can restore the original matrix as bb*bb' ≈ A. If you want to export the factor as lower-triangular, the permutation vector should be obtained as well.

using LinearAlgebra
using Test

# slightly different example
A = sparse(Float64[10 0 1 0; 0 10 1 0; 1 1 10 1; 0 0 1 10])
fac = cholesky(A)

# equivalent
@test sparse(fac.L) ≈ cholesky(Matrix(A[fac.p,fac.p])).L

# Lower triangular factor and the permutation vector
L = LowerTriangular(sparse(fac.L))
perm = fac.p
invp = invperm(perm)

# solve a system of linear equations
b = A*Float64[1,2,3,4]
@test A\b ≈ (L' \ (L \ b[perm]))[invp]

[alegarra]

right, I put this in a blog: http://artadia.blogspot.com/2019/10/sparse-cholesky-decomposition-in-julia.html