11
$\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
    $\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 '15 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$ – J. M. is in limbo May 15 '15 at 1:08
15
$\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$
7
$\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$
  • 4
    $\begingroup$ Cholesky decomposition does not apply to undetermined matrices while LDL does. Try CholeskyDecomposition[{{1, 2}, {2, 1}}] $\endgroup$ – uranix May 14 '15 at 20:15
  • 1
    $\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$ – J. M. is in limbo May 15 '15 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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