SparseArrays.jl

SparseArrays.jl copied to clipboard

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