Cholesky with sparse matrices

Should the following example work?

julia> using LinearAlgebra

julia> using SparseArrays

julia> # 10x10 sparse matrix
       A = sprand(10, 10, 1.0)
10×10 SparseMatrixCSC{Float64, Int64} with 100 stored entries:
 0.904426  0.425966   0.678877  0.990534    …  0.00486748  0.625129   0.28918
 0.908373  0.0387709  0.17868   0.725458       0.245102    0.0862697  0.632291
 0.361639  0.387367   0.457368  0.369276       0.0896446   0.958725   0.170507
 0.992062  0.344254   0.546643  0.780792       0.996618    0.427766   0.328221
 0.877559  0.19974    0.123014  0.00153193     0.516185    0.372909   0.37089
 0.396713  0.279064   0.855201  0.0253617   …  0.404998    0.934266   0.935675
 0.656733  0.0671018  0.859334  0.427824       0.828616    0.920709   0.101263
 0.682545  0.709925   0.173347  0.13505        0.351449    0.509649   0.263236
 0.667382  0.247777   0.20485   0.37823        0.249513    0.711105   0.919538
 0.185091  0.271876   0.971401  0.671663       0.956294    0.766601   0.786353

julia> # factorization with dense matrix
       F = cholesky(Matrix(A'A))
Cholesky{Float64, Matrix{Float64}}
U factor:
10×10 UpperTriangular{Float64, Matrix{Float64}}:
 2.24998  0.858844  1.30242   1.43303    …   1.32864     1.63633     1.30012
  ⋅       0.689552  0.437529  0.0572014      0.132471    0.763297    0.261955
  ⋅        ⋅        1.27031   0.51108        0.79303     0.99325     0.542768
  ⋅        ⋅         ⋅        0.855286      -0.202569   -0.298744   -0.130929
  ⋅        ⋅         ⋅         ⋅             0.0814408   0.101939    0.183558
  ⋅        ⋅         ⋅         ⋅         …  -0.409493    0.350895    0.311987
  ⋅        ⋅         ⋅         ⋅            -0.0795898   0.0724291   0.901631
  ⋅        ⋅         ⋅         ⋅             0.803321   -0.0278404   0.239681
  ⋅        ⋅         ⋅         ⋅              ⋅          0.487393   -0.116216
  ⋅        ⋅         ⋅         ⋅              ⋅           ⋅          0.301848

julia> # lower triangular is a valid AbstractMatrix
       F.L * rand(10)
10-element Vector{Float64}:
 0.6138698436089721
 0.5057798279914256
 0.7973250659194446
 0.7477808021543471
 0.49600451802224116
 0.42332422693787697
 0.7360032500412224
 1.2133500519668785
 1.2063489971211445
 1.4239143574981488

julia> # factorization with sparse matrix
       F = cholesky(A'A)
SuiteSparse.CHOLMOD.Factor{Float64}
type:    LLt
method:  simplicial
maxnnz:  55
nnz:     55
success: true


julia> # doesn't work, the L field doesn't implement getindex
       F.L * rand(10)
ERROR: getindex not defined for SuiteSparse.CHOLMOD.FactorComponent{Float64, :L}
Stacktrace:
  [1] error(::String, ::Type)
    @ Base ./error.jl:42
  [2] error_if_canonical_getindex(::IndexCartesian, ::SuiteSparse.CHOLMOD.FactorComponent{Float64, :L}, ::Int64, ::Int64)
    @ Base ./abstractarray.jl:1231
  [3] getindex
    @ ./abstractarray.jl:1217 [inlined]
  [4] _getindex
    @ ./abstractarray.jl:1257 [inlined]
  [5] getindex
    @ ./abstractarray.jl:1218 [inlined]
  [6] generic_matvecmul!(C::Vector{Float64}, tA::Char, A::SuiteSparse.CHOLMOD.FactorComponent{Float64, :L}, B::Vector{Float64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra ~/Desktop/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:759
  [7] mul!
    @ ~/Desktop/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:81 [inlined]
  [8] mul!
    @ ~/Desktop/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:275 [inlined]
  [9] *(A::SuiteSparse.CHOLMOD.FactorComponent{Float64, :L}, x::Vector{Float64})
    @ LinearAlgebra ~/Desktop/julia-1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:51
 [10] top-level scope
    @ REPL[41]:2

The cholesky documentation mentions a special case for sparse matrices, but I can’t find out how to extract the lower triangular or the F.PtL as AbstractMatrix for further processing.

Calling sparse(F.L) gives you the sparse Cholesky factor as a SparseMatrixCSC as opposed to the SuiteSparse object, which does implement the methods you would expect and has the * method (as in, the SparseMatrixCSC implements those things). Is that what you’re looking for?

1 Like

(Note that a sparse Cholesky factorization is not the same as dense Cholesky, because in the sparse case it first permutes the rows/columns of the matrix before factorizing, in order to maximize the sparsity of the factors.)

1 Like

An explicit call to sparse solves the issue, thank you @cgeoga.

I wonder if this special case with sparse matrices could be made more “Julian”. Shouldn’t F.L always return a valid AbstractMatrix as a Julia object? Can multiple dispatch help here? Why the sparse case works so differently? Could this end-user API be improved somehow? How to write generic code that works with both dense and sparse matrices?

2 Likes