Factorize a sparse symmetric indefinite matrix

[sloisel]

(I was told to open this issue here.) This seems weird:

sebastienloisel@macbook-pro julia % julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.9.3 (2023-08-24)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using LinearAlgebra, SparseArrays

julia> factorize(sparse([0.0 1.0;1.0 0.0]))
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] #ldlt!#11
   @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/solvers/cholmod.jl:1311 [inlined]
 [2] ldlt!
   @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/solvers/cholmod.jl:1303 [inlined]
 [3] ldlt!(F::SparseArrays.CHOLMOD.Factor{Float64}, A::Hermitian{Float64, SparseMatrixCSC{Float64, Int64}}; shift::Float64, check::Bool)
   @ SparseArrays.CHOLMOD /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/solvers/cholmod.jl:1336
 [4] ldlt!
   @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/solvers/cholmod.jl:1336 [inlined]
 [5] factorize
   @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/linalg.jl:1655 [inlined]
 [6] factorize(A::SparseMatrixCSC{Float64, Int64})
   @ SparseArrays /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/linalg.jl:1633
 [7] top-level scope
   @ REPL[2]:1

I didn't put more details like versioninfo() because the cause is clear. factorize doesn't pivot when doing an ldlt decomposition, which means that factorize has a weird failure mode, because if you perturb the matrix slightly so that it does lu decomposition, that one does pivot and will most likely succeed.

Connected to this, the documentation to LinearAlgebra.ldlt is pretty bad. Follow me along if you please.

A = sparse([1 2;2 -4])
F = ldlt(A)
sparse(F.Lt)

SparseArrays.CHOLMOD.CHOLMODException("Lt not supported for sparse LDLt matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")
...

This is the only way I found to get this highly informative error message printed. The documentation also says something about this but it's very confusing.

sparse(F.L)

SparseArrays.CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations")
...

sparse(F.U)

CanonicalIndexError: getindex not defined for SparseArrays.CHOLMOD.FactorComponent{Float64, :U}
...

The only one I found that works is LD = sparse(F.LD). However, only people who are expert in sparse linear algebra would understand $LD$, because this matrix is not the product of $L$ and $D$, instead it is the matrix $L$ with its diagonal overwritten with $D$.

The following explains some of the weirdness but maybe is also part of the problem:

fieldnames(typeof(F))

(:ptr,)

To make things worse, I found absolutely no documentation to the permutation option for the ldlt function.
The entirely undocumented "fieldname" F.p does something but I have no way of knowing if it works because I don't know how to ask ldlt to use pivoting.

[djturizo]

I encountered the same issues, and yeah, the lack of documentation is pretty frustrating. After some searching I found that the sparse solvers used by Julia are just wrappers of the SuiteSparse package. It is not clear from the SuiteSparse documentation either, but the package doesn't have a generic LDL' solver. In the author's own words:

I don't have a true LDL' solver. I can compute and LDL' factorization only for a diagonal D with nonzero diagonal entries, and it's not suitable for all symmetric indefinite matrices. In particular, for my simple LDL' solver inside CHOLMOD, and also my LDL package, the matrix must have the property that all leading submatrices are well conditioned. I only have it because it's the format required for the update/downdate in CHOLMOD.

If Julia wants to use my LDL' factorization directly, it should be checking the condition number estimate. It doesn't look like it's doing that.

So basically the LDL' solver provided by SuiteSparse does not do pivoting, which is a key ingredient to successfully factor symmetric indefinite matrices (like in the Bunch-Kaufman and Aasen factorization methods for dense matrices).

I wonder if there are any efficient, open source LDL' solvers for sparse matrices that could be wrapped for use in Julia...