Skip to content

Conversation

ajefweiss
Copy link

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
Copy link
Collaborator

alexhroom commented May 14, 2025

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
Copy link
Author

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
Copy link
Collaborator

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants