LU Decomposition of Invertible Matrices
Buildup
When the matrix Matrix $A \in \mathbb{R}^{m \times m}$ is multiplied on the left side such that the component $(i, j)$ becomes $0$, the matrix $E_{ij}$ is defined as the $A$ Elimination Operator for $ij$.
Specifically, for the square matrix Square Matrix $(a_{ij}) \in \mathbb{R}^{m \times m}$, $E_{ij}$ has the diagonal components as $1$, the $(i,j)$ component as $\displaystyle -m_{ij} = -{{a_{ij}} \over {a_{jj}}}$, and the remaining components as $0$. This represents the operation of eliminating the same variables to match the coefficients in the solution of simultaneous equations in matrix form. Such elimination is meaningful only when there exists only one solution to the simultaneous equations, so the $A$ we deal with is assumed to be an Invertible Matrix.
Example
For a simple example, $$ B := \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix} $$ then we get $E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$. And since $B = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$, the $(2,1)$ elimination operator for $B$ can be calculated as $$ E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ -m_{ij} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -{{1} \over {2}} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ Meanwhile, repeated operations of this type on $$ U_{B} : = E_{32} E_{31} E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 2 \end{bmatrix} $$ results in an Upper Triangular Matrix, and all $E_{21} , E_{31} , E_{32}$ are Lower Triangular Matrices. Since the product of lower triangular matrices is a lower triangular matrix, $L_{B} : = E_{21} E_{31} E_{32}$ is also a lower triangular matrix and can be represented as $L^{-1} B = U$. Rearranging for $B$ yields $B = L_{B} U_{B}$, decomposing into the product of a lower and upper triangular matrix.
This decomposition method can be generalized without issue even for the case of $A \in \mathbb{R}^{n \times n}$ due to its simplicity. Thinking about the solution method for the simultaneous equations $A \mathbf{x} = \mathbf{b}$, since $A$ can be decomposed into $A = LU$, $$ \begin{align*} LU \mathbf{x} =& \mathbf{b} \\ \mathbf{y} :=& U \mathbf{x} \end{align*} $$ can be stated. Here, since $U$ is an upper triangular matrix, $\mathbf{y}$ can be easily found using back substitution, and $$ \begin{align*} L (U \mathbf{x} ) =& L \\ \mathbf{y} =& \mathbf{b} \end{align*} $$ since $L$ is a lower triangular matrix, $\mathbf{b}$ can be easily found using forward substitution. Therefore, being able to actually calculate LU decomposition means that a wide variety of problems can be solved quite easily and quickly.
Derivation
Now, to directly find $L$ and $U$, let’s define the lower triangular matrix $L$ that satisfies $A := L U$ and the Upper Triangular Matrix $U$ as follows. $$ \begin{align*} 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} \\ U : =& \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{bmatrix}\end{align*} $$ Here, it can easily be verified that all the diagonal components of $L$ are $1$. According to the above definition, calculating the component $a_{ij}$ of $A = (a_{ij})$ yields $$ a_{ij} = \sum_{s = 1}^{n} l_{is} u_{sj} = \sum_{s = 1}^{min(i , j)} l_{is} u_{sj} $$ Now, let’s consider $a_{ij}$ depending on whether it is on the diagonal, above the diagonal, or below the diagonal.
If $i = j = k$,
$$ \begin{align*} a_{kk} =& \sum_{s = 1}^{\min(i, j)} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + l_{kk} u_{kk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + u_{kk} \end{align*} $$
Arranging this for $u_{kk}$,
$$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$
can be obtained. Similarly, if $k < j$, $$ a_{kj} = \sum_{s = 1}^{k-1} l_{ks} u_{sj} + u_{kj} $$ and if $k < i$, $$ a_{ik} = \sum_{s = 1}^{k-1} l_{is} u_{sk} + l_{ik} u_{kk} $$ can be obtained. Arranging these for $u_{kj}$ and $l_{ik}$ respectively yields: $$ \begin{align*} u_{kj} =& a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} \\ l_{ik} =& {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} \end{align*} $$
■
Algorithm
Assuming $(a_{ij}) \in \mathbb{R}^{n \times n}$ is an invertible matrix.
Step 1. $k = 1$
Substitute $u_{1j} = a_{1j}$ and calculate $\displaystyle l_{i1} = {{1} \over {u_{11}}} a_{i1}$.
Step 2. $k = 2, 3, \cdots , n-1$
- Step 2-1. Calculate the following. $$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$
- Step 2-2. For $j = k+1 , k+2 , \cdots , n-1$, calculate the following: $$ u_{kj} = a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} $$
- Step 2-3. For $i = k+1 , k+2 , \cdots , n-1$, calculate the following: $$ l_{ik} = {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} $$
Step 3. For $k = n$, calculate the following: $$ u_{nn} = a_{nn} - \sum_{s = 1}^{n-1} l_{ns} u_{sn} $$
- Careful examination of the algorithm reveals that as $k$ increases, $L$ is calculated from left to right and $U$ from top to bottom. The necessary values at each step are ensured to have been calculated in the previous steps, making the calculation unit at each step either one row or one column.
Code
Let’s actually decompose the matrix $B$ used in the example above.
Below is the actual implementation of the above algorithm in R.
B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
LU<-function(A)
{
n=length(A[,1])
L<-diag(1,n)
U<-diag(0,n)
k=1
j=1
i=1
U[1,]<-A[1,]
L[,1]<-A[,1]/A[1,1]
for(k in 2:(n-1))
{
U[k,k] = A[k,k] - sum(L[k,1:(k-1)]*U[1:(k-1),k])
for(j in (k+1):n)
{
U[k,j] = A[k,j] - sum(L[k,1:(k-1)]*U[1:(k-1),j])
}
for(i in (k+1):n)
{
L[i,k] = ( A[i,k] - sum(L[i,1:(k-1)]*U[1:(k-1),k]) ) / U[k,k]
}
}
U[n,n] = A[n,n] - sum(L[n,1:(n-1)]*U[1:(n-1),n])
return(list(L=L,U=U))
}
LU(B)
require(Matrix)
expand(lu(B))
The result is as follows.
It can be confirmed that the LU decomposition implemented in R itself matches exactly.