matrix - $mathbf Lmathbf Dmathbf 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?

Answer

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