How to obtain U^-1 factor from sparse Cholesky factorization?

Specific Domains Numerics
, ,
1

Given a dense positive semi-definite symmetric matrix A, we can easily obtain the U^{-1} factor from the Cholesky factorization object:

julia> using LinearAlgebra

julia> A = Symmetric([2 1; 1 2])
2×2 Symmetric{Int64, Matrix{Int64}}:
 2  1
 1  2

julia> F = cholesky(A)
Cholesky{Float64, Matrix{Float64}}
U factor:
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 1.41421  0.707107
  ⋅       1.22474

julia> inv(F.U)
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 0.707107  -0.408248
  ⋅         0.816497

How to obtain the equivalent result with sparse matrices?

julia> using SparseArrays

julia> F = cholesky(sparse(A))
SparseArrays.CHOLMOD.Factor{Float64, Int64}
type:    LLt
method:  simplicial
maxnnz:  3
nnz:     3
success: true


julia> sparse(F.U)
ERROR: CanonicalIndexError: getindex not defined for SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64}
Stacktrace:
  [1] error_if_canonical_getindex(::IndexCartesian, ::SparseArrays.CHOLMOD.FactorComponent{…}, ::Int64, ::Int64)
    @ Base ./abstractarray.jl:1327
  [2] getindex
    @ ./abstractarray.jl:1311 [inlined]
  [3] _getindex
    @ ./abstractarray.jl:1358 [inlined]
  [4] getindex
    @ ./abstractarray.jl:1312 [inlined]
  [5] iterate
    @ ./abstractarray.jl:1209 [inlined]
  [6] iterate
    @ ./abstractarray.jl:1207 [inlined]
  [7] SparseMatrixCSC{Float64, Int64}(M::SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64})
    @ SparseArrays ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/sparsematrix.jl:891
  [8] SparseMatrixCSC
    @ ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/sparsematrix.jl:885 [inlined]
  [9] convert
    @ ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/sparsematrix.jl:980 [inlined]
 [10] sparse(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64})
    @ SparseArrays ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/SparseArrays/src/sparsematrix.jl:1016
 [11] top-level scope
    @ REPL[116]:1
Some type information was truncated. Use `show(err)` to see complete types.

I tried to convert the various factors of the factorization object to sparse matrices, but only sparse(F.L) works.

I have also tried the following according to suggestions in error messages:

julia> F.U \ I
ERROR: CanonicalIndexError: getindex not defined for SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64}
Stacktrace:
 [1] error_if_canonical_getindex(::IndexCartesian, ::SparseArrays.CHOLMOD.FactorComponent{…}, ::Int64, ::Int64)
   @ Base ./abstractarray.jl:1327
 [2] getindex
   @ ./abstractarray.jl:1311 [inlined]
 [3] factorize(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/dense.jl:1436
 [4] inv(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1061
 [5] \(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :U, Int64}, J::UniformScaling{Bool})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.4+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/uniformscaling.jl:277
 [6] top-level scope
   @ REPL[138]:1
Some type information was truncated. Use `show(err)` to see complete types.
2

Can’t you do U = sparse(F.L)'? Beware that the Cholesky factorization doesn’t mean quite the same thing in the sparse case, because the matrix is permuted first (to reduce fill-in).

As for the inverse, the inverse of a sparse matrix is (generically) dense, so you might as well convert to Matrix (or UpperTriangular) before calling inv. (That’s why inv is not implemented for sparse matrices.)

3

I will go with this option. I wonder if sparse(F.U) should work though. Perhaps an issue in LinearAlgebra?

Thank you for the suggestion. Definitely an improvement over the full Matrix.