Add LDL decomposition #1515

ajefweiss · 2025-05-06T18:21:43Z

See issue #865, this decomposition is also mentioned in #762 but in the end only the UDU decomposition was implemented (PR #766). I need this for my own use case where I basically want a Cholesky decomposition of a semi-definite covariance matrix to generate random correlated numbers. Currently I was applying the Cholesky decomposition on a modified matrix with all zero rows/columns removed, which are then later re-added to the lower triangular matrix.

I used the code from #766 as a template, but changed the code to implement the LDL algorithm instead, everything else is basically the same just re-named. This is not the standard algorithm, but a variant that should also work with semi-definite matrices (i.e. with eigenvalues equal to zero), which I need. I haven't added any tests for this yet though. I do not know why the existing UDU decomposition also does not include this capability as its just an extra 1-2 lines of code.

I also added a cholesky_l function that reconstructs the lower triangular matrix from the Cholesky decomposition.

alexhroom · 2025-05-14T09:46:15Z

hi, looking forward to this being added - I was wondering whether it's necessary for it to insist on the scalar type having trait RealField? In the UDU factorisation PR it's set to RealField because nobody was sure whether the UDU decomposition made sense for complex Hermitian matrices, but LDL is pretty widely used for complex matrices.

I tested it out on a few complex Hermitian matrices and got the results I was expecting. Maybe the functionality could be extended to ComplexField and you add some tests for equality with Cholesky on positive-definite complex Hermitian examples?

ajefweiss · 2025-05-16T13:53:43Z

I added functionality for Hermitian matrices, and rewrote the tests.

Currently, the algorithm checks for diagonal values that are exactly zero, to zero out the corresponding columns/rows. One could add a threshold value eps, similar to the way its done for the pseudo_inverse() function, so that values below eps are set to zero. Only problem here is that the diagonal values have no real relation to the eigenvalues, so the interpretation is not clear.

alexhroom · 2025-05-16T17:45:11Z

@ajefweiss cool! for clarity i'm not a maintainer/reviewer of nalgebra i'm just interested in having this in the package :)

geo-ant · 2025-11-13T08:23:34Z

Hi @ajefweiss, I'd be willing to review this but I need your help 🙂 I skimmed the paper linked in #762, but I couldn't find an algorithm for the LDL decomposition. I see there's a screenshot in the issue for UDU decomposition. Could you point me to the reference that tells me how you implemented the LDL decomposition, so I can read up on it?

EDIT: well, I scrolled down in the issue and there I see this. Is that the reference you used?

ajefweiss · 2025-11-13T16:45:43Z

Hi @geo-ant!

Thank you for taking a look at this. The actual reference I used is just the wikipedia entry which has imo a nicer description of the algorithm (https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2) which cites "Fundamentals of Matrix Computations" from Watkins (PDF versions are easily obtainable online). The relevant pages here are 82-84 although the algorithm here is just explained as a variant of the LU decomposition without being explicitly written out.

The only problem is that I have not been able to find an explicit reference for an algorithm that allows for semi-positive definite matrices. Both the reference you linked, and mine, require positive-definite matrices in their algorithms due to the division by v_i / d_i. Thus when one of the diagonal entries is 0 the algorithm fails. From my understanding, the A=LL/LDL factorizations are unique if, and only if, the matrix A is positive-definite. The theorems from both references state that the LDL factorization does exist only under the requirement that A is symmetric (and not necessarily positive definite).

I think the factorization is no longer unique because you can fill out the relevant column in the L matrix with any values as they will be subsequently zero'd out in the multiplication. In my code this column is filled with zeros.

PS: I might try to think a better name for the function "cholesky_ldl" as I actually don't think that name is appropriate and might be a bit confusing.

geo-ant

thanks for the great work, see some comments in the review.

You said that this algorithm should also exist for positive semi definite matrices, but your implementation only covers positive definite ones. As long as we document this, I don't see a problem here.

geo-ant · 2025-11-19T14:41:05Z

src/linalg/ldl.rs

+use crate::dimension::Dim;
+use simba::scalar::ComplexField;
+
+/// LDL factorization.


could you expand this to include more about the pre-and-postconditions? For example:

a sentence like "LDL factorization of a (hermitian?) matrix A into a lower unit-triangular (??) matrix L and a diagonal matrix D such that A = LDL^T".

what assumptions are made about the matrix A. I think it must be Hermitian?

It's fine to keep it concise, but I find it helpful to have all the important info right there from a user's perspective.

Thinking about this some more, do you think this decomposition should be called LDL or LDLT? I'm unsure...

just want to add to make sure we clarify that D is strictly diagonal - people also use LDL to refer to a different algorithm which decomposes a Hermitian matrix into a lower unit-triangular matrix L and a block-diagonal matrix D composed of 1x1 and 2x2 blocks.

This is the algorithm implemented in the *HETRF family of algorithms in LAPACK, which are wrapped by SciPy's ldl function in Python as well as MATLAB's ldl function, so users may be expecting the block-diagonal algorithm and be confused. (the idea of the block-diagonal algorithm is that there are bounds on the size of the elements of D which do not exist for the strictly diagonal algorithm)

The Rust package faer has both versions - they call the strictly diagonal one LDLT and the block-diagonal one LBLT.

Thinking about this some more, do you think this decomposition should be called LDL or LDLT? I'm unsure...

for some comparisons summarising info above, both Eigen (C++) and faer call it LDLT, SciPy calls it LDL, as does MATLAB.

Thanks for the comments.

Correction: My implementation of the algorithm covers semi-positive (and semi-negative) definite matrices. The referenced algorithms in the textbooks DO NOT.

Regarding the naming convection LDL/LDLT. As @alexhroom pointed out, both conventions are used. Given the fact that we already have the UDU factorization implemented as UDU (and not UDU^T), I think it would make sense to keep it as LDL.

geo-ant · 2025-11-19T14:44:29Z

src/linalg/ldl.rs

+    ///
+    /// The input matrix `p` is assumed to be hermitian c and this decomposition will only read
+    /// the lower-triangular part of `p`.
+    pub fn new(p: OMatrix<T, D, D>) -> Option<Self> {


this is absolutely a me-problem, but is there a reason why the matrix is called p, rather than mat or a? p always makes me think of a permutation matrix like in AP = QR... 😅

We can change this, currently the naming is just taken from the code from the UDU factorization.

geo-ant · 2025-11-19T14:50:11Z

src/linalg/ldl.rs

+    }
+
+    /// Returns the lower triangular matrix as if generated by the Cholesky decomposition.
+    pub fn cholesky_l(&self) -> OMatrix<T, D, D> {


I'd rather not have this in there since it's easy to google how to obtain an LL^T decomposition from an LDL^T decomposition. I don't have super strong feelings about this, but the naming is a bit awkward and the use case seems pretty fringe. But if you tell me this absolutely needs to be in there, I'll concede.

Maybe a better name for this function would be lsqrtd?

One of the main uses cases of this factorization is to get the same result as from Cholesky but without the same constraints.

geo-ant · 2025-11-20T06:44:29Z

tests/linalg/ldl.rs

+                proptest! {
+                    #[test]
+                    fn ldl(m in dmatrix_($scalar)) {
+                        let m = &m * m.adjoint();


If I understood correctly what you said in the comment, then this implementation of the LDL^T decomposition requires the matrix to be symmetric positive definite (not semi definite), right? However, using M M^T will only give us pos. semidefinite, so we should expect some spurious failures here. I've actually run into the same problem in nalgebra-lapack and you can see my solution here. I've added a proptest function that gives a spd matrix in a bit of a hacky, but sound way: The idea is to use M' = M M^T + alpha * Id, where alpha is chosen suitably large to make we don't get into numerical trouble. You should be able to easily adapt the functions to return complex matrices.

just wanted to weigh in - i'm a little confused actually... the 'point' of LDL is that it's a generalisation of the Cholesky algorithm to all symmetric/Hermitian matrices, with no requirements on positive definiteness/semi-definiteness. If it only works for definite/semidefinite matrices then there's no benefit to using LDL over Cholesky, which is a simpler decomposition and a more efficient algorithm.

Sorry if the above text is slightly confusing. But no, the LDL/LDL^T decomposition only requires the matrix to be symmetric. It does NOT require it to be positive definite. This is stated in both textbooks referenced here.

The problem is that, in both cases, the algorithms that are shown DO require the matrix to be positive OR negative definite. I'd have to look deeper into this, but I think this is because in the semi-positive/negative definite cases the solution is not unique. I think there is a choice for the particular column where the corresponding value in the D matrix is zero.

geo-ant · 2025-11-20T07:03:20Z

src/linalg/ldl.rs

+                if d[j] == T::zero() {
+                    l[(i, j)] = T::zero();
+                } else {
+                    l[(i, j)] = l_ij / d[j].clone();
+                }


is it actually correct to do that or should we return None on division by zero?

See comment above about the choice for semi-positive/negative definite matrices.

sorry for being a bit dense here (no pun intended), but are you saying that your algorithm works for all symmetric positive or negative semidefinite matrices? What if you changed the algorithm to use the in-place algorithm as suggested below? That algorithm requires A to have an LU decomp, which I believe precludes semi-definiteness...

geo-ant · 2025-11-20T07:32:32Z

src/linalg/ldl.rs

+        for j in 0..n {
+            let mut d_j = p[(j, j)].clone();
+
+            if j > 0 {
+                for k in 0..j {
+                    d_j -= l[(j, k)].clone() * l[(j, k)].clone().conjugate() * d[k].clone();
+                }
+            }
+


I've thought about the algorithm a bit and I'm wondering if you'd mind implementing the algorithm in #762 (comment). The reason is, that algorithm should be a more efficient in terms of memory usage (and that might make it more cache friendly, but I'm not sure) than your implementation. The reason is that the linked algorithm works in place, while your algorithm works out of place. Your interface wouldn't even have to change because you're taking the OMatrix by value already. Your original matrix A would then contain the LDL decomposition where L is in the stricly lower triangular part and d is on the diagonal.

The only problem I see then, is that exposing access to L and D would require new allocations. Is that why you implemented the algorithm out of place in the first place?

What I mean is this: is the primary use of LDL^T decomposition to get explicit access to L and D, or would it also be used to solve linear systems. In that case we could expose a solve method that uses the triangular solvers internally and deals with D through scaling.

We could do this, but we should then change the structure a bit. LDL could then just be a new type of a single matrix without exposing the underlying field as L*D by itself does not have a use for anything. The current implementation for UDU uses the exactly same set up.

I'd prefer that, unless you say the main value of the LDL(T) decomposition is access to the L and D matrix. I'd still expose access to those matrices via dedicated methods that allocate, because who know what people might use this for. I'd also add a solve method for solving linear systems in a least squares sense, that takes advantage of the internal structure.

Code has ben re-factored to make the calculations in-place, I'm still using the same algorithm just not allocating new matrices.

ajefweiss added 2 commits May 6, 2025 13:51

basic ldl implementation, based on udu

5592c9f

test for cholesky eq

b875935

hermitian ldl, better tests

00ca97b

ajefweiss mentioned this pull request Nov 12, 2025

Add covariance matrices? #1575

Open

geo-ant self-requested a review November 13, 2025 08:23

geo-ant added the enhancement label Nov 13, 2025

geo-ant reviewed Nov 20, 2025

View reviewed changes

ajefweiss added 2 commits November 25, 2025 11:29

update docs, renamed cholesky_l to lsqrt and added extra checks

d562fbd

in place ldl factorization

43fd280

Uh oh!

Add LDL decomposition #1515

Are you sure you want to change the base?

Add LDL decomposition #1515

Conversation

ajefweiss commented May 6, 2025

Uh oh!

alexhroom commented May 14, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ajefweiss commented May 16, 2025

Uh oh!

alexhroom commented May 16, 2025

Uh oh!

geo-ant commented Nov 13, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ajefweiss commented Nov 13, 2025

Uh oh!

geo-ant left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

alexhroom Nov 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

geo-ant Nov 20, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

alexhroom commented May 14, 2025 •

edited

Loading

geo-ant commented Nov 13, 2025 •

edited

Loading

alexhroom Nov 23, 2025 •

edited

Loading

geo-ant Nov 20, 2025 •

edited

Loading