# Cholesky with sparse matrices

**URL:** <https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494>\
**Category:** Numerics\
**Tags:** question, linearalgebra, factorization\
**Created:** [December 2, 2021, 9:47pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494 "2021-12-02T21:47:13Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 2, 2021, 9:47pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/1 "2021-12-02T21:47:13Z")

</div>

Should the following example work?

```julia
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](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.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.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [December 2, 2021, 10:56pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/2 "2021-12-02T22:56:02Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 3, 2021, 1:07am UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/3 "2021-12-03T01:07:47Z")

</div>

(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.)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 3, 2021, 1:27am UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/4 "2021-12-03T01:27:15Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![jarapo](https://avatars.discourse-cdn.com/v4/letter/j/6bbea6/32.png) [@jarapo](https://discourse.julialang.org/u/jarapo)\
**Post date:** [November 2, 2022, 6:28pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/5 "2022-11-02T18:28:23Z")

</div>

I got stuck with a similar problem and wish to add something.  
Indeed  
`sparse(F.L) * rand(10)`  
works.

Also, following the OP,

```julia
F = cholesky(Matrix(A'A))
F.U * rand(10)

```

works but  
`sparse(F.U) * rand(10)`  
does not:

```julia
ERROR: getindex not defined for SuiteSparse.CHOLMOD.FactorComponent{Float64, :U}
Stacktrace:
  [1] error(::String, ::Type)
    @ Base ./error.jl:42
  [2] error_if_canonical_getindex(::IndexCartesian, ::SuiteSparse.CHOLMOD.FactorComponent{Float64, :U}, ::Int64, ::Int64)
    @ Base ./abstractarray.jl:1183
  [3] getindex
    @ ./abstractarray.jl:1169 [inlined]
  [4] _getindex
    @ ./abstractarray.jl:1214 [inlined]
  [5] getindex
    @ ./abstractarray.jl:1170 [inlined]
  [6] iterate
    @ ./abstractarray.jl:1096 [inlined]
  [7] iterate
    @ ./abstractarray.jl:1094 [inlined]
  [8] SparseMatrixCSC{Float64, Int64}(M::SuiteSparse.CHOLMOD.FactorComponent{Float64, :U})
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:616
  [9] convert
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:667 [inlined]
 [10] sparse(A::SuiteSparse.CHOLMOD.FactorComponent{Float64, :U})
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:692
 [11] top-level scope
    @ REPL[96]:1

```

I was working with a non-sparse matrix and the code was working (I chose U because that saved me mentally transposing a matrix in my application).

Changed it to a sparse matrix and took so long where it went wrong (now I am just using L’ instead of U).

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [September 28, 2023, 8:27am UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/6 "2023-09-28T08:27:51Z")

</div>

```julia
using SparseArrays
using LinearAlgebra
A = sprand(10, 10, 1.0)
FFF = cholesky(sparse(A'A))
@show FFF
FFF.L \ rand(10) # ok
FFF.PtL \ rand(10) #ok
sparse(FFF.L) *rand(10) # ok on 1.8.5& 1.9.3 ; ko on 1.7
sparse(FFF.PtL) * rand(10) #ko on all versions 1.7, 1.8.5, 1.9.3 

```

The “sparse” trick seems to work only with the .L factor not the more interesting .PtL factor, for which the error is

```julia
CanonicalIndexError: getindex not defined for SuiteSparse.CHOLMOD.FactorComponent{Float64, :PtL}

```

Is there a way to be able to work with the .PtL components as with normal matrix ?

---

<div class="post-metadata">

**Author:** ![Alexander-Barth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alexander-barth/32/3692_2.png) [@Alexander-Barth](https://discourse.julialang.org/u/Alexander-Barth)\
**Post date:** [September 28, 2023, 2:34pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/7 "2023-09-28T14:34:34Z")

</div>

I had a look at the cholmod.jl code and it seems that one can get the permutation matrix via `FFF.p`

> <https://github.com/JuliaSparse/SparseArrays.jl/blob/ac5c8ed5e87af863d62892693f683543b122bb76/src/solvers/cholmod.jl#L1189>

With `FFF.L` and the permutation `FFF.p`, I guess you can construct PtL ?

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [September 28, 2023, 2:46pm UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/8 "2023-09-28T14:46:34Z")

</div>

Thanks, yes that gives access to the permutation. I guess I do not even need to form PtL but just use the permutations whenever I do a matrix multiplications using L

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [September 29, 2023, 6:31am UTC](https://discourse.julialang.org/t/cholesky-with-sparse-matrices/72494/9 "2023-09-29T06:31:02Z")

</div>

For completness, how to avoid needing PtL:

```julia
using SparseArrays
using LinearAlgebra
A = sprand(10, 10, 1.0)
FFF = cholesky(sparse(A'A))
x=rand(10)
L= sparse(FFF.L) 
p=FFF.p

work=L*(L'*x[p])
work=work[invperm(p)]

norm(x-FFF\work)/norm(x)

```