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?
$\begingroup$
$\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$– NasserCommented 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$– J. M.'s missing motivation ♦Commented May 15, 2015 at 1:08
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
5
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, 2015 at 6:37
J. M.'s missing motivation♦J. M.'s missing motivation
125k11 gold badges407 silver badges581 bronze badges
-
1$\begingroup$ A variation for Hermitian matrices can be found here. $\endgroup$ Commented Feb 2, 2020 at 16:39
-
$\begingroup$ Show how this algorithm works for a matrix containing symbols $\endgroup$– ayrCommented Dec 7, 2021 at 17:51
-
1$\begingroup$ You could try feeding it a symbolic matrix and study the output. $\endgroup$ Commented 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$– ayrCommented 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];,l1is the $\mathbf L$ factor, andDiagonalMatrix[d1]is the $\mathbf D$ factor. $\endgroup$ Commented Dec 7, 2021 at 18:05
$\begingroup$
$\endgroup$
2
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, 2015 at 20:14
-
5$\begingroup$ Cholesky decomposition does not apply to undetermined matrices while LDL does. Try
CholeskyDecomposition[{{1, 2}, {2, 1}}]$\endgroup$– uranixCommented 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$ Commented May 15, 2015 at 3:39