LDL Decomposition Calculator

What is LDL Decomposition?

LDL decomposition (also called LDLᵀ factorization) is a method of decomposing a symmetric matrix A into three factors: a unit lower triangular matrix L, a diagonal matrix D, and the transpose of L. The factorization is written as:

Here, "unit lower triangular" means L has ones on its main diagonal and all entries above the diagonal are zero. The matrix D is diagonal, meaning all off-diagonal entries are zero. Unlike Cholesky decomposition, the diagonal entries of D are allowed to be negative or zero, which makes LDL decomposition applicable to a broader class of symmetric matrices, including symmetric indefinite matrices.

Every symmetric matrix that can be factored without pivoting admits an LDLᵀ decomposition. More precisely, if A is symmetric and all leading principal minors are nonzero, then A has a unique LDLᵀ factorization. When combined with symmetric pivoting (row and column permutations of the form PAPᵀ), the factorization extends to all symmetric matrices: PAPᵀ = LDLᵀ.

How LDL Relates to Cholesky Decomposition

LDL and Cholesky decomposition are closely related. In fact, Cholesky decomposition is a special case of LDL decomposition. If A is symmetric positive definite, then all diagonal entries of D are strictly positive, and we can write:

where L̃ = L√D is the Cholesky factor. In other words, Cholesky decomposition absorbs the square root of D into L to produce a single lower triangular factor. This means LDL decomposition is the "square-root free" variant of Cholesky: it avoids computing any square roots, which can be advantageous both for numerical stability and for performance on certain hardware.

The key advantage of LDL over Cholesky is generality. Cholesky requires the matrix to be symmetric positive definite (all eigenvalues positive). If any eigenvalue is negative or zero, the square root √D is not real, and Cholesky fails. LDL decomposition, by contrast, happily accommodates negative diagonal entries in D, making it suitable for symmetric indefinite matrices that arise frequently in optimization, mechanics, and constrained systems.

The LDL Algorithm Step by Step

The algorithm computes L and D one column at a time. For an n×n symmetric matrix A, the procedure is:

1. Initialize: Set L to the identity matrix and D to a zero matrix.
2. For each column j = 1, 2, ..., n:
   1. Compute D[j][j]: The j-th diagonal entry of D is given by
      D[j][j] = A[j][j] − Σ_{k=1}^{j-1} L[j][k]² · D[k][k]
   2. Compute L[i][j] for i > j: Each entry below the diagonal in column j is
      L[i][j] = (1 / D[j][j]) · (A[i][j] − Σ_{k=1}^{j-1} L[i][k] · L[j][k] · D[k][k])
3. Result: After processing all columns, L is unit lower triangular and D is diagonal, with A = LDLᵀ.

Note that the algorithm accesses only the lower triangle of A, since A is symmetric. Each column requires referencing the previously computed columns of L and entries of D, making this an inherently sequential column-by-column process.

A Worked Example

Consider the 3×3 symmetric matrix:

Column 1: D[1][1] = 4. Then L[2][1] = 2/4 = 0.5, and L[3][1] = -2/4 = -0.5.

Column 2: D[2][2] = 10 - (0.5)²(4) = 10 - 1 = 9. Then L[3][2] = (1/9)(2 - (-0.5)(0.5)(4)) = (1/9)(2 + 1) = 1/3.

Column 3: D[3][3] = 5 - (-0.5)²(4) - (1/3)²(9) = 5 - 1 - 1 = 3.

The final factorization is:

When to Use LDL vs Cholesky vs LU

Choosing the right decomposition depends on the properties of your matrix and the problem you are solving:

• LDL decomposition: Best for symmetric matrices, especially when the matrix may be indefinite (has both positive and negative eigenvalues). Common in optimization where Hessian matrices are indefinite, in saddle-point systems, and in structural mechanics with mixed formulations. LDL avoids square roots and costs about n³/3 operations.
• Cholesky decomposition: Best for symmetric positive definite matrices. Cholesky is slightly simpler to implement and is the standard choice when you know the matrix is positive definite (e.g., covariance matrices, normal equations in least squares, kernel matrices). It also costs n³/3 operations but requires square roots.
• LU decomposition: The general-purpose factorization for any non-singular square matrix. LU does not require symmetry and uses partial pivoting for stability. It costs 2n³/3 operations, roughly twice as much as LDL or Cholesky. Use LU when the matrix has no special structure.

In summary: if your matrix is symmetric positive definite, use Cholesky. If it is symmetric but possibly indefinite, use LDL. If it has no symmetry, use LU.

Applications of LDL Decomposition

LDL decomposition appears across many areas of applied mathematics and engineering:

• Optimization and Newton's method: The Hessian matrix of a twice-differentiable function is symmetric. At a saddle point or away from a minimum, the Hessian can be indefinite. LDL decomposition factors the Hessian efficiently and reveals its inertia (the number of positive, negative, and zero eigenvalues), which is critical for trust-region and line-search methods.
• Interior point methods: Modern solvers for linear and quadratic programming use interior point algorithms that require solving symmetric indefinite systems of the form known as augmented systems or KKT systems. LDL decomposition is the natural factorization for these systems.
• Structural and finite element analysis: Stiffness matrices in structural mechanics are symmetric. When constraint equations or mixed formulations are involved, the resulting system is symmetric indefinite, making LDL the preferred factorization.
• Covariance and precision matrices: In statistics, updating a covariance or precision matrix (for instance during Kalman filtering) can be done efficiently with rank-one updates to an existing LDL factorization, avoiding repeated full factorizations.
• Solving symmetric linear systems: Any symmetric system Ax = b can be solved via LDLᵀ: first decompose, then solve Ly = b by forward substitution, then Dz = y (trivial diagonal solve), then Lᵀx = z by back substitution. This requires roughly n³/3 flops for the factorization and 2n² for each subsequent solve.

Numerical Stability and Pivoting

The basic LDL algorithm without pivoting can encounter zero or near-zero diagonal entries D[j][j], leading to division by zero or severe numerical instability. To handle this, symmetric pivoting strategies are employed.

The most widely used approach is the Bunch-Kaufman algorithm, which uses 1×1 and 2×2 pivots. At each step, the algorithm examines the current column and decides whether to use a single diagonal element as a pivot (1×1 pivot) or to combine two rows/columns into a 2×2 block pivot. The 2×2 pivots allow the algorithm to handle cases where no single diagonal element is large enough to serve as a stable pivot. The resulting factorization takes the form PAPᵀ = LDLᵀ where D is block diagonal with 1×1 and 2×2 blocks.

Another pivoting strategy is the Bunch-Parlett algorithm, which searches the entire unreduced submatrix for the best pivot at each step. This provides better numerical guarantees (bounded element growth) but is more expensive because of the full search. In practice, Bunch-Kaufman is almost always preferred for its lower overhead and strong empirical stability.

The bounded Bunch-Kaufman algorithm, implemented in LAPACK routines like dsytrf, guarantees that the elements of L are bounded by a modest constant (approximately 2.781), ensuring that round-off errors do not grow uncontrollably during the factorization.

Computational Complexity

LDL decomposition of an n×n symmetric matrix requires approximately n³/3 floating-point operations. This is the same cost as Cholesky decomposition and exactly half the cost of LU decomposition (which requires 2n³/3 operations). The savings come from exploiting symmetry: only the lower triangle of the matrix needs to be accessed and stored.

In terms of storage, LDL decomposition requires only n(n+1)/2 entries for L (the lower triangle, with the unit diagonal implicitly stored) plus n entries for D, totaling n(n+1)/2 + n = n(n+3)/2 entries. In practice, L can overwrite the lower triangle of A in place, and D can be stored separately in a vector of length n.

Once the factorization is computed, solving a linear system Ax = b costs 2n² operations (forward substitution, diagonal solve, and back substitution). This makes LDL very efficient when the same matrix A is used with multiple right-hand sides.

Inertia of a Matrix

One of the most valuable byproducts of LDL decomposition is the inertia of the matrix. The inertia of a symmetric matrix A is the triple (n⁺, n⁻, n⁰) counting the number of positive, negative, and zero eigenvalues respectively.

By Sylvester's law of inertia, the factorization A = LDLᵀ preserves inertia: the signs of the diagonal entries of D exactly match the signs of the eigenvalues of A. Specifically:

• The number of positive entries in D equals the number of positive eigenvalues of A.
• The number of negative entries in D equals the number of negative eigenvalues of A.
• The number of zero entries in D equals the number of zero eigenvalues (the nullity) of A.

This makes LDL decomposition an efficient tool for determining whether a symmetric matrix is positive definite (all D entries positive), negative definite (all negative), or indefinite (mixed signs), without computing the eigenvalues directly. This is particularly useful in optimization, where checking the definiteness of a Hessian matrix determines whether a critical point is a minimum, maximum, or saddle point.

Connection to Bunch-Kaufman Factorization

The Bunch-Kaufman factorization is a pivoted variant of LDL decomposition designed for robust numerical computation. While the standard LDL decomposition produces a strictly diagonal D, the Bunch-Kaufman factorization allows D to be block diagonal with 1×1 and 2×2 blocks.

The factorization takes the form PAPᵀ = LDLᵀ where P is a permutation matrix and D has the block structure. The 2×2 blocks appear when no suitable 1×1 pivot can be found, which occurs when the diagonal elements are too small relative to the off-diagonal elements. Each 2×2 block can be individually factored or inverted as needed during the solve phase.

Bunch-Kaufman is the algorithm implemented by LAPACK (dsytrf, ssytrf), MATLAB's ldl() function, and most production numerical libraries. It requires the same n³/3 flop count as plain LDL but provides strong backward-stability guarantees. For most practical problems, the Bunch-Kaufman factorization is the recommended way to compute an LDL-type decomposition.

Limitations

• Requires the input matrix to be symmetric (or Hermitian for complex matrices). Non-symmetric matrices should use LU decomposition instead.
• Without pivoting, the factorization can fail if any leading principal minor is zero, even if the matrix is nonsingular.
• The factorization does not directly provide eigenvalues or eigenvectors; it only reveals the inertia (eigenvalue sign counts).
• For sparse matrices, fill-in during factorization can be significant, and specialized reordering strategies (e.g., approximate minimum degree) are needed to preserve sparsity.