14
$\begingroup$

Mathematica can do a Cholesky decomposition $\mathbf A = \mathbf L\mathbf L^\top$, but how do I do a LDL decomposition $\mathbf A = \mathbf L\mathbf D\mathbf L^\top$, with $\mathbf L$ being a unit lower triangular matrix?

$\endgroup$
2
  • 2
    $\begingroup$ Here is the LAPACK Fortran implemenation using what is called Bunch-Kaufman diagonal pivoting method for LDL^T factorization for real symmetric (not necessarily positive definite) matrix. I think Mathematica should have an LDL special decomposition. Matlab has one here. $\endgroup$
    – Nasser
    May 14, 2015 at 20:51
  • $\begingroup$ @Nasser, in fact the $\mathbf D$ matrix in Bunch-Kaufman is block diagonal instead of being diagonal, so I wouldn't strictly consider it as something the OP is looking for. $\endgroup$ May 15, 2015 at 1:08

2 Answers 2

19
$\begingroup$

I needed this decomposition to answer another question, so I broke down and implemented it myself. The code is more or less a straightforward translation of the pseudocode in Golub/Van Loan:

LDLT[mat_?SymmetricMatrixQ] := 
     Module[{n = Length[mat], mt = mat, v, w},
            Do[
               If[j > 1,
                  w = mt[[j, ;; j - 1]]; v = w Take[Diagonal[mt], j - 1];
                  mt[[j, j]] -= w.v;
                  If[j < n,
                     mt[[j + 1 ;;, j]] -= mt[[j + 1 ;;, ;; j - 1]].v]];
               mt[[j + 1 ;;, j]] /= mt[[j, j]],
               {j, n}];
            {LowerTriangularize[mt, -1] + IdentityMatrix[n], Diagonal[mt]}]

A few tests:

m1 = HilbertMatrix[20];
m2 = Array[Min, {20, 20}];

{l1, d1} = LDLT[m1];
m1 == l1.DiagonalMatrix[d1].Transpose[l1]
   True

{l2, d2} = LDLT[m2];
m2 == l2.DiagonalMatrix[d2].Transpose[l2]
   True
$\endgroup$
5
  • 1
    $\begingroup$ A variation for Hermitian matrices can be found here. $\endgroup$ Feb 2, 2020 at 16:39
  • $\begingroup$ Show how this algorithm works for a matrix containing symbols $\endgroup$
    – dtn
    Dec 7, 2021 at 17:51
  • 1
    $\begingroup$ You could try feeding it a symbolic matrix and study the output. $\endgroup$ Dec 7, 2021 at 17:53
  • $\begingroup$ Is it possible to remove the word Conjugate? And yet, is it possible to display the matrix L and the matrix D separately? $\endgroup$
    – dtn
    Dec 7, 2021 at 17:54
  • 1
    $\begingroup$ The version here assumes a symmetric matrix (not a Hermitian one), so it doesn't use Conjugate[]. The $\mathbf L$ and $\mathbf D$ matrices are already separate, but are stored together in a list. If you evaluate {l1, d1} = LDLT[m1];, l1 is the $\mathbf L$ factor, and DiagonalMatrix[d1] is the $\mathbf D$ factor. $\endgroup$ Dec 7, 2021 at 18:05
10
$\begingroup$

I'll illustrate with a simple example.

mat = {{2, 1}, {1, 3.}};
ch = CholeskyDecomposition[mat]

(* Out[145]= {{1.41421356237, 0.707106781187}, {0., 1.58113883008}} *)

Pull out the diagonal. Use it to modify and get a triangular matrix with ones on the diagonal.

diag = Diagonal[ch]

(* Out[148]= {1.41421356237, 1.58113883008} *)

modch = ch*1/diag

(* Out[149]= {{1., 0.5}, {0., 1.}} *)

Since we have to account for two such factors (one on each side), the D matrix will be the square of this diagonal. We check below that this gives the correct decomposition.

Transpose[modch].DiagonalMatrix[diag^2].modch

(* Out[153]= {{2., 1.}, {1., 3.}} *)
$\endgroup$
2
  • 5
    $\begingroup$ Cholesky decomposition does not apply to undetermined matrices while LDL does. Try CholeskyDecomposition[{{1, 2}, {2, 1}}] $\endgroup$
    – uranix
    May 14, 2015 at 20:15
  • 2
    $\begingroup$ @uranix means that a symmetric matrix that is not positive definite will certainly not have a Cholesky decomposition, but it may still have an $\mathbf L\mathbf D\mathbf L^\top$ decomposition. $\endgroup$ May 15, 2015 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.