Sparse Symmetric LDLT

[mfalt]

I just created a SparseLDLT.jl package to be able to factorize for example symmetric quasi-definite matrices that ldltfact is currently unable to handle. The code is very simple and using SuiteSparse/LDL to do the hard work (which is already packed with julia but not compiled). Maybe it could be useful to include in a future version of julia somehow? I've seen a lot of discussion about sparse factorizations but couldn't find the relevant issue. I still don't know exactly the set of matrices that SuiteSparse/LDL is successful on.

I had to use the name spldltfact since I couldn't overload ldltfact in Base, a simple example:

using SparseLDLT
m = 400;  n = 500
A = sprandn(m, n, .2)
M = [speye(n) A'; A -speye(m)]
b = randn(m+n)
F = spldltfact(M)
x = F\b

ldltfact usually results in

julia> ldltfact(M)
ERROR: ArgumentError: matrix has one or more zero pivots

[ararslan]

I think one of the long-run goals here is to reduce the size of base Julia, and have certain packages be part of a standard library. To that end, I think it makes more sense for this to continue to live in a separate package.

[andreasnoack]

@mfalt I think you just revealed a bug in ldltfact. See https://github.com/JuliaLang/julia/pull/19045. I took a look at the documentation for LDL and, as I read it, it's not supposed to handle more situations than CHOLMOD.

[mfalt]

Interesting, although this means that I spent a couple of days on this in vain. Would you mind explaining to me which matrices these two functions are guaranteed to work for?

[andreasnoack]

This is from the CHOLMOD documentation: "It also includes a non-supernodal LDLt factorization method that can factorize symmetric indefinite matrices if all of their leading submatrices are well-conditioned". It seems that the LDL code is mainly for teaching the theory: "Its primary purpose is to illustrate much of the basic theory of sparse matrix algorithms". I believe it covers the same class of matrices as CHOLMOD's LDLt. CHOLMOD is just much more optimized. Does that answer your question?

[andreasnoack]

@mfalt Is the fix in #19045 sufficient for your application? Can this issue be closed?

[mfalt]

Yes, I've been testing it a bit and it seems to perform well. I'm not sure that it is always faster, but can't confirm that. Thanks for fixing this!