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

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. will be back soon May 15 '15 at 1:08

2 Answers

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


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

• Cholesky decomposition does not apply to undetermined matrices while LDL does. Try CholeskyDecomposition[{{1, 2}, {2, 1}}] – uranix May 14 '15 at 20:15
• @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. will be back soon May 15 '15 at 3:39