Cholesky with sparse matrices

Specific Domains Numerics
, ,
1

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.

2

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
3

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

2 Likes
4

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
5

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

Also, following the OP,

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

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

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

6 
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

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 ?

7

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

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

1 Like
8

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

1 Like
9

For completness, how to avoid needing PtL:

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)