Matrix Factorization: LU Decomposition

LU (short for lower-upper) decomposition is a form of matrix factorization that decomposes a matrix into a product of lower and upper triangular matrices. This matrix factorization approach was introduced by the Polish astronomer Tadeusz Banachiewicz in 1938. In a way, LU decomposition can be seen as a Gaussian elimination (row reduction) method executed in matrix form. The LU decomposition is one of the most widely used matrix factorization methods for solving a system of linear equations (realize that such problems arise in many different fields including engineering, physics, economics, computer science, and applied mathematics), especially when such a system results in an equation $\mathbf{A}\mathbf{x}=\mathbf{b}$ where $\mathbf{A}\in\mathbb{R}^{n\times n}$ is square and known, $\mathbf{b}\in\mathbb{R}^n$ is given, and $\mathbf{x}\in\mathbb{R}^n$ is the variable to be solved.

In this article, we will discuss the steps to execute an LU decomposition algorithm on a square matrix $\mathbf{A}$ such that it can be written as $\mathbf{A} = \mathbf{L} \mathbf{U}$ for a lower triangular matrix $\mathbf{L}\in\mathbb{R}^{n\times n}$ and an upper triangular matrix $\mathbf{U}\in\mathbb{R}^{n\times n}$ . However, note that, when $\mathbf{A}$ has at least a zero on the main diagonal, then one needs to interchange either some of its row or column, and this will involve a row permutation matrix $\mathbf{P}\in\mathbb{R}^{n\times n}$ or a column permutation matrix $\mathbf{Q}\in\mathbb{R}^{n\times n}$ that are not equal to the identity matrix. However, for the sake of simplicity, let us assume for now that the main diagonal of $\mathbf{A}$ has no zero entries.

The Wikipedia page on LU decomposition is an excellent starting point for understanding the mechanics of this factorization method. Following some of the notations, let us assume that we have a square matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ having no zero on the diagonal. As such, denote $\mathbf{A}^{(i)}$ as the $i-$ th version of $\mathbf{A}$ where all of the elements below its main diagonal have been eliminated to zero through the row reduction method for the first $i-$ th columns. Then, at step $i$ , we have $\mathbf{A}^{(i)} = \mathbf{L}_{i}^{-1}\mathbf{A}^{(i-1)}$ where $\mathbf{L}_{i}^{-1}$ is an identity matrix with the exception that its entries below the diagonal at the i-$th column are (most likely) nonzero. Starting from $\mathbf{A}^{(0)} = \mathbf{A}$ , then we carried out the following steps until $n-1$ :

$\mathbf{A} = \mathbf{L}_{1}\mathbf{L}_{1}^{-1}\mathbf{A}^{(0)}=\mathbf{L}_{1}\mathbf{A}^{(1)}$ $\Leftrightarrow \mathbf{A} = \mathbf{L}_{1}\mathbf{L}_{2}\mathbf{L}_{2}^{-1}\mathbf{A}^{(1)}=\mathbf{L}_{1}\mathbf{L}_{2}\mathbf{A}^{(2)}$ $\Leftrightarrow \cdots \Leftrightarrow$ $\mathbf{A} =\mathbf{L}_{1}\mathbf{L}_{2}\hdots\mathbf{L}_{n-1}\mathbf{A}^{(n-1)},$

which suggest that $\mathbf{L} \triangleq \mathbf{L}_{1}\mathbf{L}_{2}\hdots\mathbf{L}_{n-1}$ and $\mathbf{U} \triangleq\mathbf{A}^{(n-1)}.$

To illustrate how the above makes sense in the LU decomposition algorithm, let us consider a simple example where $\mathbf{A}\in\mathbb{R}^{3\times 3}$ is given as:

$\mathbf{A}= \begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0.5&2&2\end{bmatrix}.$

Firstly, we do the row reduction on the first column. This is performed as follows:

$\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&2.5&2\end{bmatrix} = \begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0.5&2&2\end{bmatrix} - \begin{bmatrix}0& 0 & 0\\\frac{0}{1}(1) & \frac{0}{1}(-1)& \frac{0}{1}(0)\\\frac{1}{2}(1) & \frac{1}{2}(-1)& \frac{1}{2}(0)\end{bmatrix},$

which is equal to:

$\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&2.5&2\end{bmatrix}}_{\mathbf{A}^{(1)}}= \underbrace{\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\-0.5&0&1\end{bmatrix}}_{\mathbf{L}^{-1}_1}\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0.5&2&2\end{bmatrix}}_{\mathbf{A}^{(0)}}.$

Thanks to the lower triangular structure of $\mathbf{L}^{-1}_1$ , then it is easy to verify that its inverse, $\mathbf{L}_1$ , is equal to:

$\mathbf{L}_1=\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0.5&0&1\end{bmatrix}.$

Secondly, we repeat the same step but on the second column. This results in:

$\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&0&1.1667\end{bmatrix} = \begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&2.5&2\end{bmatrix} - \begin{bmatrix}0& 0 & 0\\0 & 0& 0\\\frac{5}{6}(0) & \frac{5}{6}(3)& \frac{5}{6}(1)\end{bmatrix},$

which is equal to:

$\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&0&1.1667\end{bmatrix}}_{\mathbf{A}^{(2)}}= \underbrace{\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0&-\frac{5}{6}&1\end{bmatrix}}_{\mathbf{L}^{-1}_2}\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&2.5&2\end{bmatrix}}_{\mathbf{A}^{(1)}}.$

Again, it is straightforward to verify that $\mathbf{L}_2$ is simply equal to:

$\mathbf{L}_2=\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0&\frac{5}{6}&1\end{bmatrix}.$

That’s it! Now we are ready to compute the LU decomposition. The matrix $\mathbf{L}$ is equal to:

$\mathbf{L}= \mathbf{L}_1\mathbf{L}_2= \begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0.5&\frac{5}{6}&1\end{bmatrix}= \begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0.5&0.8333&1\end{bmatrix},$

whereas $\mathbf{U}$ is nothing but:

$\mathbf{U}= \mathbf{A}^{(2)}= \begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&0&1.1667\end{bmatrix}.$

At last, we finally have:

$\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0.5&2&2\end{bmatrix}}_{\mathbf{A}}= \underbrace{\begin{bmatrix}1 & 0 & 0\\0 & 1 & 0\\0.5&0.8333&1\end{bmatrix}}_{\mathbf{L}}\underbrace{\begin{bmatrix}1 & -1 & 0\\0 & 3 & 1\\0&0&1.1667\end{bmatrix}}_{\mathbf{U}}.$

Pretty easy, right? Now I want to mention that the process above, which is called the pure LU decomposition, will only work when there is no zero on the main diagonal, or in more technical terms, all its leading principal minors are nonzero. A generalization of LU decomposition is the LUP decomposition, where the letter P here denotes the row permutation matrix. LUP is simply an LU decomposition but with partial pivoting (row). That is, LU decomposition is equal to LUP decomposition in which the permutation matrix is equal to the identity matrix. The aforementioned Wikipedia page on the LU decomposition neatly covers an example of how the LUP decomposition is performed.

Hopefully, this short article can help you in understanding the LU decomposition. For the next articles, I am planning on explaining how the QR decomposition works and comparing the performance between LU and QR decomposition to solve some systems of linear equations on an embedded system. Stay tuned!

Author: Sebastian Nugroho

I am an engineer and researcher working on the field of control theory and optimization. View all posts by Sebastian Nugroho

Share this:

Related

Author: Sebastian Nugroho

Leave a comment Cancel reply