📂Matrix Algebra

Decomposition of Symmetric Matrices into LDU

Overview

$L^{T}$ is an upper triangular matrix, so we can consider replacing $U$ in $A = LU$ with $U:= DL^{T}$. As general conditions become stricter than for the LU decomposition, the amount of computation significantly decreases.

Theorem

For the invertible symmetric matrix $A$, there exist a lower triangular matrix $L$ and a diagonal matrix $D$ satisfying $A = LDL^{T}$.

Proof

Since $A$ is an invertible matrix, there exist a lower triangular matrix $L$ and an upper triangular matrix $U$ that satisfy $A = L U$. On the other hand, as $A = A^{T}$ holds,

$$ L U = A = A^{T} = U^{T} L^{T} $$

In the LU decomposition, $L$ is obtained as a lower triangular matrix with all diagonal elements as $1$, so $\det L = 1$, indicating that the inverse matrix $L^{-1}$ exists. Since $L^{-1}$ is a lower triangular matrix and $$ L U = U^{T} L^{T} $$ is rearranged for $U$, $$ U = L^{-1} U^{T} L^{T} $$ multiplying both sides by $(L^{T})^{-1}$ on the right, $$ U (L^{T})^{-1} = L^{-1} U^{T}$$ is obtained. Since the left side of $U (L^{T})^{-1} = L^{-1} U^{T}$ is an upper triangular matrix and the right side is a lower triangular matrix, $$ D: = U (L^{T})^{-1} = L^{-1} U^{T} $$ is a diagonal matrix. Rearranging this for $U$, it becomes $U = DL^{T}$, and we obtain the following: $$ A = LU = LDL^{T} $$

■

Calculation

Now, to specifically calculate $L$ and $D$, let’s define as follows:

$$ L : = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & 1 \end{bmatrix},\quad D : = \begin{bmatrix} d_{11} & 0 & 0 & \cdots & 0 \\ 0 & d_{22} & 0 & \cdots & 0 \\ 0 & 0 & d_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & d_{nn} \end{bmatrix} $$

Let’s directly calculate $A = (a_{ij}) = LDL^{T}$.

Since $l_{ij} > 0$ and $l_{ii} = 1$,

$$ \begin{align*} a_{ij} =& \sum_{\nu = 1}^{n} l_{i \nu} \sum_{\mu = 1}^{n} d_{\nu \mu} l_{\mu j}^{T} \\ =& \sum_{\nu = 1}^{n} l_{i \nu} \sum_{\mu = 1}^{n} \delta_{\nu \mu} d_{\nu \mu} l_{ j \mu} \\ =& \sum_{\nu = 1}^{n} l_{i \nu} d_{\nu \nu} l_{ j \nu} \\ =& \sum_{\nu = 1}^{ \min (i,j) } l_{i \nu} d_{\nu \nu} l_{ j \nu} \end{align*} $$

First, let’s assume the case where $j \le i$.

$$ a_{ij} = \sum_{\nu = 1}^{j} l_{i \nu} d_{\nu \nu} l_{j \nu} = \sum_{\nu = 1}^{j-1} l_{i \nu} d_{\nu \nu} l_{j \nu} + l_{ij} d_{jj} l_{jj} = \sum_{\nu = 1}^{j-1} l_{i \nu} d_{\nu \nu} l_{j \nu} + l_{ij} d_{jj} $$

Here, in the case where $i = j = k$,

$$ a_{kk} = \sum_{\nu = 1}^{k-1} l_{k \nu} d_{\nu \nu} l_{k \nu} + d_{kk} d_{kk} $$

rearranging it,

$$ d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2 $$

Conversely, assuming $j = k$ and $i = k+1, \cdots , n$,

$$ a_{ik} = \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} + l_{ik} d_{kk} $$

rearranging $l_{ik}$,

$$ l_{ik} = {{1} \over {d_{kk}}} \left\{ a_{ik} - \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} \right\} $$

Algorithm

Suppose $(a_{ij}) \in \mathbb{R}^{n \times n}$ is an invertible symmetric matrix.

**Step 1. $k = 1$

Substitute $d_{11} = a_{11}$ and calculate $\displaystyle l_{i1} = {{1} \over {d_{11}}} a_{i1}$.

---

**Step 2. $k = 2, 3, \cdots , n-1$

• Step 2-1. Calculate the following: $$ d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2 $$

• Step 2-2. For $i = k+1 , k+2 , \cdots , n-1$, calculate the following: $$ l_{ik} = {{1} \over {d_{kk}}} \left\{ a_{ik} - \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} \right\} $$

---

Step 3. For $k = n$, calculate the following: $$ d_{nn} = a_{nn} - \sum_{\nu = 1}^{n-1} d_{\nu \nu} l_{n \nu}^2 $$

---

• Although calculations seemed complex when done manually, organizing them into an algorithm significantly reduces the computational workload. Always remember to first solve for the diagonal components of $D$ and then compute the values of $L$ column by column.

Code

Now, let’s take an example similar to the one used in LU decomposition and perform LDU decomposition.

The following code implements the above algorithm in R.

B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
 
LDU<-function(A)
{
  n=length(A[,1])
  L<-diag(1,n)
  D<-diag(0,n)
  
  k=1
  D[1,1] = A[1,1]
  L[,1] = A[,1]/D[1,1]
 
  for(k in 2:(n-1))
  {
    D[k,k] = A[k,k] - sum(diag(D)[1:(k-1)]*L[k,1:(k-1)]^2)
    for(i in (k+1):n)
    {
      L[i,k] = ( A[i,k] - sum(L[i,1:(k-1)]*diag(D)[1:(k-1)]*L[k,1:(k-1)]) )/D[k,k]
    }
  }
  
  D[n,n] = A[n,n] - sum(diag(D)[1:(n-1)]*L[n,1:(n-1)]^2)
 
  return(list(L=L, D=D))
}
 
LDU(B)
LDU(B)$L%*%LDU(B)$D%*%t(LDU(B)$L)

The results are as follows:

Since LDU decomposition isn’t separately implemented in R, I verified only the decomposition, and it was correct.

Because the size of $B$ is not very large and the computation is manageable, if in doubt, you might want to try manual calculation directly.