Cholesky/LDL-decomposition for semidefinite matrices in python

Question

I am looking for Cholesky/LDL-decomposition for semi-definite matrices in python.

Search-results:

• Both numpy.linalg.cholesky and sympy.Matrix.LDLdecomposition only work for positive-definite.
• There is a request on the scipy-mailing list from 2011 https://mail.scipy.org/pipermail/scipy-user/2011-July/029926.html, but it seems that was not implemented.
• There is a rough description on https://mathoverflow.net/questions/155147/cholesky-decomposition-of-a-positive-semi-de%EF%AC%81nite, but no code. Is there a straightforward implementation?

My instances are rather small, about $100\times100$, so a symbolic solution is fine (the bottleneck is in some other place).

Answer 1

With SciPy v1.1.0 and later, you can use scipy.linalg.ldl for the factorization of indefinite matrices.

example : Creating a positive definite matrix arr

>>> import numpy as np
>>> import scipy.linalg as la

>>> arr = np.random.rand(100,97)
>>> M = arr @ arr.T
>>> l, d, p = la.ldl(M)
>>> np.allclose(l.dot(d).dot(l.T) - M, np.zeros([100, 100]))
True

Answer 2

You can just use an LU decomposition. For symmetric or hermitian matrices they are equivalent up to some sign ambiguities.