Reality check: How useful exactly are SuiteSparse and LDLFactorizations' LDLT?

https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl/issues/101

I realize that calling for a native implementation of a multi-frontal or supernodal factorization
is a tall order: people spend decades of their life on these efforts. But I think it would be a useful
addition, truly empowering for Julia devotees who have to deal with large sparse matrices.
Perhaps a good target for some funding to make it happen!

2 Likes

Sorry, I don’t understand: what exactly is the problem interfacing fast external codes (we are doing it all the time AFAIU)?

Codes such as HLS or SuiteSparse cannot work directly with CSC (or other common format) matrices. So, the matrix is converted into internal representation recognized by the library. Hence, suddenly one has doubled the amount of storage because both the original sparse matrix and the internal rep are held in memory at the same time.

1 Like

OK, I can see how this is a problem memory wise, but is it dominating the computation (that would be surprising for me)?

As long as you have double the amount of memory that you actually need, the computation can be unaffected by the increased amount of memory requested. However,
it is quite easy to get an “Out of memory” error when using SuiteSparse. Then the speed
of the computation is immaterial. :wink:

1 Like

I see. But this sounds more like a call to support more memory layouts in the involved packages then (I wondered about the missing support for CSR and COO already)?

When dealing with sparse factorizations, it is typically the case that the factors require a lot (as in, a LOT) more memory than the initial matrix.
Thus, while it is correct that calling SuiteSparse’s LDLt or Cholesky will duplicate the input data, this increase in memory is often negligible when compared to how much memory you use for the factorization itself.

I use LDLFactorizations inside the LP solver Tulip.jl to solve KKT systems, when the arithmetic is anything other than Float64. It’s great to be able to do that out of the box.

That being said: most of the problems out there are in Float64. There are amazing libraries (like HSL, Pardiso, Mumps, etc…) with decades of development that we can already use, at the price of a minimal overhead. Unless a large share of user start solving linear systems in “exotic” arithmetic, I don’t see the need for a julia-native sparse factorization library.
More specifically: I would not advise anyone to start such an endeavor, but it someone does I’ll happily use it.

1 Like

This is a good point. Nevertheless, memory not used is memory available.

Also, the keyword is oftern “we can already use”. Most of these libraries have heavy dependencies which make installation on some architectures a serious undertaking. PETSC, MUMPS, HLS, … Try to install these on Windows… With Julia one can get used to the idea of write once, run anywhere, but these packages bust that wide open.

4 Likes

I think one issue is that efficient sparse direct factorisations actually rely on CSR format while other sparse operations (matrix-vector product) rely on CSC. So, even if one codes up a (faster?) Cholesky factorisation than Cholmod in Julia, I would expect that conversions steps (CSC to CSR, or even COO to CSC to CSR) will still be needed. I’m far from being an expert but I do use Suitesparse. I usually sparse from COO (triplets), this converts it automatically to CSC. Then (I guess) it 's internally converted to CSR prior to factorization. It would be nice to avoid some steps but are they avoidable?

Row-based storage of the sparse matrix is better for parallel execution of matrix vector products. However, I believe that Suitesparse use a column-based storage. CSparse definitely does.

OK, I got all mixed up :smiley:
Looks like UMFPACK and Cholmod do use CSC, it’s Pardiso that uses CSR (and variants in intel Pardiso).