$\mathbf L\mathbf D\mathbf L^\top$ Cholesky decomposition

Question

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?

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. — Nasser, May 14 '15 at 20:51
@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. — J. M.'s missing motivation, May 15 '15 at 01:08

score 19 · Accepted Answer · answered Aug 28 '15 at 06:37

19

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

answered Aug 28 '15 at 06:37

J. M.'s missing motivation

124,525
11
401
574

1

A variation for Hermitian matrices can be found here. – J. M.'s missing motivation Feb 02 '20 at 16:39
Show how this algorithm works for a matrix containing symbols – dtn Dec 07 '21 at 17:51
1

You could try feeding it a symbolic matrix and study the output. – J. M.'s missing motivation Dec 07 '21 at 17:53
Is it possible to remove the word Conjugate? And yet, is it possible to display the matrix L and the matrix D separately? – dtn Dec 07 '21 at 17:54
1

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. – J. M.'s missing motivation Dec 07 '21 at 18:05

score 10 · Answer 2 · answered May 14 '15 at 20:14

10

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.}} *)

answered May 14 '15 at 20:14

Daniel Lichtblau

58,970
2
101
199

5

Cholesky decomposition does not apply to undetermined matrices while LDL does. Try CholeskyDecomposition[{{1, 2}, {2, 1}}] – uranix May 14 '15 at 20:15
2

@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. – J. M.'s missing motivation May 15 '15 at 03:39

$\mathbf L\mathbf D\mathbf L^\top$ Cholesky decomposition

2 Answers2

Linked

Related