Enhancement: Add LDLt factorization to scipy

[daniel-perry]

It would be useful to be able to use an LDLt decomposition in scipy (similar to http://www.mathworks.com/help/matlab/ref/ldl.html).

This has been discussed and partially developed for numpy in numpy/numpy#3960.
From that discussion it appears a simple wrapper around an LAPACK function for LDLt would suffice, and someone has submitted some code to do that, though no pull request has been submitted.

However (also discussed there) it appears that ideally new linear algebra routines should go into scipy instead of numpy ?

Which has lead to this issue being opened.

Please advise if this sounds appropriate for scipy and if so I am happy to either move @tmwright 's changes over from numpy/numpy#3960 over to scipy (once he gives permission) or start from scratch on a new wrapper.

Referencing @seberg so he is aware of this new ticket in scipy.

[rgommers]

Yes, Scipy sounds like the right place to me. (rationale given in numpy/numpy#3960 (comment))

[srkunze]

I noticed a missing LDL decomposition in scipy as well, so I wrote one myself trying to vectorize as much as possible (only 1 for-loop left). :)

The resulting Python function performs a decomposition of a 2000x2000 matrix in 4.56 secs whereas the scipy default cholesky decomposition finishes in 1.61 secs. I know that both algorithms are not 100% comparable but I think it's enough to get an idea of how worse the performance would be.

Does a naive LDL decomposition besides wrappers for LAPACK makes sense for inclusion in scipy? If so, I would start a pull request.

[argriffing]

I've naively implemented and tested LDL decomposition, solving, and inversion in arb_mat in response to flintlib/arb#110 (comment). arb itself is LGPL but I think I could give permission for code derived from the LDL part to be used in a less restrictively licensed project.

[sturlamolden]

A "naïve" implementation makes no sence. It will be slow and error prone. Use LAPACK if you can for scipy.linalg. If you cannot find a suitable LAPACK routine, write a solver on top of BLAS. In this case (LDL' and UDU' decompositions) the LAPACK *SYTRF routines will do what we need, so preferably use those.

[sturlamolden]

The main coding task will be to parse the output of *SYTRF, which is really gross.

[tmwright]

@sturlamolden Please see my remarks over on numpy/numpy#3960.

[srkunze]

@sturlamolden I would hardly call two lines of code error-prone. ;)

I just thought it would be a good idea to give something back to SciPy as a way to say "thank you".

@tmwright Do you mean both implementations (with and without pivoting) might be useful to have in SciPy?