6.1. LU decomposition

LU decomposition (also known as LU factorisation) is a procedure for decomposing a square matrix \(A\) into the product of a lower triangular matrix \(L\) and an upper triangular matrix \(U\) such that

\[ A = LU. \]

The advantage of writing a matrix as a product of \(L\) and \(U\) is that the solution to a triangular set of equations is easy to calculate using forward or back substitution. Consider the LU decomposition of a \(3\times 3\) matrix

\[\begin{split} \begin{align*} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix}= \begin{pmatrix} \ell_{11} & 0 & 0\\ \ell_{21} & \ell_{22} & 0\\ \ell_{31} & \ell_{32} & \ell_{33} \end{pmatrix}\begin{pmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{pmatrix}, \end{align*} \end{split}\]

which gives a system of 9 equations (one for each element in \(A\)) in 12 unknowns. Since there are more unknowns than equations, the factorisation is not unique. A standard convention is to set the diagonal entries of \(L\) equal to 1 (\(\ell_{ii}=1\)), which produces a unique decomposition whenever it exists.

\[\begin{split} \begin{align*} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} &= \begin{pmatrix} 1 & 0 & 0 \\ \ell_{21} & 1 & 0 \\ \ell_{31} & \ell_{32} & 1 \end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{pmatrix} \\ &= \begin{pmatrix} u_{11} & u_{12} & u_{13} \\ \ell_{21} u_{11} & u_{22} + \ell_{21} u_{12} & u_{23} + \ell_{21} u_{13} \\ \ell_{31} u_{11} & \ell_{31} u_{12} + \ell_{32} u_{22} & u_{33} + \ell_{31} u_{13} +\ell_{32} u_{23} \end{pmatrix}. \end{align*} \end{split}\]

In general,

\[ a_{ij} = \sum_{k=1}^{\min(i,j)} \ell_{ij}u_{kj}. \]

The triangular structure of \(L\) and \(U\) limits which terms are non-zero, leading to separate formulas for \(i > j\) and \(i < j\).

The elements in the lower triangular region, where \(i>j\), are

\[ \begin{align*} a_{ij} =\sum_{k=1}^j \ell_{ik} u_{kj} =\ell_{ij} u_{jj} +\sum_{k=1}^{j-1} \ell_{ik} u_{kj}, \end{align*} \]

which is rearranged to

\[ \begin{align*} \ell_{ij} =\frac{1}{u_{jj} }\left(a_{ij} -\sum_{k=1}^{j-1} \ell_{ik} u_{kj} \right). \end{align*} \]

For the elements in the upper triangular region where \(i\le j\) we have

\[ \begin{align*} a_{ij} = u_{ij} +\sum_{k=1}^{i-1} \ell_{ik} u_{kj}, \end{align*} \]

which is rearranged to

\[ \begin{align*} u_{ij} = a_{ij} -\sum_{k=1}^{i-1} \ell_{ik} u_{kj}. \end{align*} \]

So to calculate an LU decomposition we loop through each column of \(A\) and calculate the elements of \(\ell_{ij}\) and \(u_{ij}\) for that column.

Theorem 6.1 (LU decomposition)

The LU decomposition of an \(n \times n\) square matrix \(A\) results in two \(n \times n\) matrices \(L\) and \(U\) such that \(A = LU\). The elements of \(L\) and \(U\) are calculated using

(6.1)\[\begin{split} \begin{align*} u_{ij} &= a_{ij} - \sum_{k=1}^{i-1} \ell_{ik}u_{kj}, & i &= 1, \ldots, j, \\ \ell_{ij} &= \dfrac{1}{u_{jj}} \left(a_{ij} - \displaystyle \sum_{k=1}^{j-1} \ell_{ik}u_{kj}\right), & i &= j+1, \ldots, n. \end{align*} \end{split}\]

The quantity \(u_{jj}\) is called the pivot. The algorithm fails if a pivot is zero.

Note

Not every square matrix admits an LU decomposition without row exchanges. The formulas below require the pivots \(u_{jj}\) to be non-zero. If a zero pivot is encountered, row interchanges are required, leading to the more general \(PA=LU\) factorisation discussed later.

Algorithm 6.1 (LU decomposition)

Inputs: A square \(n \times n\) matrix \(A\).

Outputs: Two square \(n \times n\) matrices \(L\) and \(U\) such that \(A = LU\).

• \(L \gets I_n\)

• \(U \gets \mathbf{0}_{n \times n}\)

• for \(j = 1, \ldots, n\)

  • \(u_{ij} \gets a_{ij} - \displaystyle\sum_{k=1}^{i-1} \ell_{ik}u_{kj}, \qquad i = 1, \ldots, j\)

  • \(\ell_{ij} \gets \dfrac{1}{u_{jj}} \left( a_{ij} - \displaystyle \sum_{k=1}^{j-1} \ell_{ik}u_{kj} \right), \qquad i = j + 1, \ldots, n\)

• Return \(L\) and \(U\)

The computational cost of LU decomposition is \(O(n^3)\), comparable to Gaussian elimination. Once the factorisation has been computed, multiple systems \(A \mathbf{x} = \mathbf{b}\) can be solved efficiently using back substitution.

Example 6.1

Determine the LU decomposition of the following matrix

\[\begin{split} A = \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix}.\end{split}\]

---

Solution

Loop through the columns of \(A\) and compute the elements of \(L\) and \(U\)

\[\begin{split} \begin{align*} j &= 1: & u_{11} &= a_{11} = 2, \\ && \ell_{21} &= \frac{1}{u_{11}}a_{21} = \frac{1}{2}(-4) = -2, \\ && \ell_{31} &= \frac{1}{u_{11}}a_{31} = \frac{1}{2}(6) = 3, \\ j &= 2: & u_{12} &= a_{12} = 3, \\ && u_{22} &= a_{22} - \ell_{21}u_{12} = -7 - (-2)(3) = -1, \\ && \ell_{32} &= \frac{1}{u_{22}}(a_{32} - \ell_{31}u_{12}) = \frac{1}{-1}(7 - 3(3)) = 2, \\ j &= 3: & u_{13} &= a_{13} = 1, \\ && u_{23} &= a_{23} - \ell_{21}u_{13} = 0 - (-2)(1) = 2, \\ && u_{33} &= a_{33} - \ell_{31}u_{13} - \ell_{32}u_{23} = 10 - 3(1) - 2(2) = 3, \end{align*} \end{split}\]

therefore

\[\begin{split} \begin{align*} L &= \begin{pmatrix} 1 & 0 & 0\\ -2 & 1 & 0\\ 3 & 2 & 1 \end{pmatrix},& U &= \begin{pmatrix} 2 & 3 & 1 \\ 0 & -1 & 2 \\ 0 & 0 & 3 \end{pmatrix}. \end{align*} \end{split}\]

Checking that \(LU=A\)

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & 0\\ -2 & 1 & 0\\ 3 & 2 & 1 \end{pmatrix} \begin{pmatrix} 2 & 3 & 1 \\ 0 & -1 & 2 \\ 0 & 0 & 3 \end{pmatrix} = \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix}. \quad \checkmark \end{align*} \end{split}\]

---

6.1.1. Solving systems of linear equations using LU decomposition

Given a system of linear equations of the form \(A \mathbf{x} = \mathbf{b}\) where LU decomposition has been applied to the coefficient matrix then since \(A = LU\) we have

(6.2)\[ LU \mathbf{x} = \mathbf{b}.\]

Introduce an intermediate vector \(\mathbf{y}\) defined by

\[ L \mathbf{y} = \mathbf{b}. \]

Since \(L\) is lower triangular, the solution of \(L \mathbf{y} = \mathbf{b}\) is easily computed using forward substitution. Once \(\mathbf{y}\) has been determined, the solution \(\mathbf{x}\) is obtained from

\[ U \mathbf{x} = \mathbf{y}, \]

using back substitution.

Algorithm 6.2 (Solving systems of linear equations using LU decomposition)

Inputs: A system of linear equations expressed using the coefficient matrix \(A\) and variable vector \(\mathbf{b}\) such that \(A \mathbf{x} = \mathbf{b}\).

Outputs: The solution vector \(\mathbf{x}\).

• Compute the LU decomposition of \(A\)

• Solve \(L \mathbf{y} = \mathbf{b}\) using forward substitution

  • \(y_1 = b_1\)

  • \(y_i \gets b_i - \displaystyle\sum_{j=1}^{i-1} \ell_{ij} y_j\), \(\qquad i = 2, \ldots, n\)

• Solve \(U \mathbf{x} = \mathbf{y}\) using back substitution

  • \(x_n \gets \dfrac{y_n}{u_{nn}}\)

  • \(x_i \gets \dfrac{1}{u_{ii}} \left( y_i - \displaystyle \sum_{j=i+1}^{n} u_{ij} x_j \right)\), \(\qquad i = n-1, \ldots, 1\)

• Return \(\mathbf{x}\).

The main advantage of LU decomposition arises when solving several systems

\[ \begin{align*} A \mathbf{x}^{(1)} &= \mathbf{b}^{(1)}, & A \mathbf{x}^{(2)} &= \mathbf{b}^{(2)}, & \ldots \end{align*} \]

with the same coefficient matrix \(A\). The decomposition \(A = LU\) is only computed once, after which each new right-hand side requires only one forware-substitution step and one back-substitution step.

This is why LU decomposition is commonly used in implicit Runge-Kutta methods. During the Newton iterations, the coefficient matrix \(I−h(A \otimes J_f)\) remains unchanged while the right-hand side vector \(g(\mathbf{Y})\) varies. Consequently, the LU decomposition is computed once and reused for each linear solve, leading to a substantial reduction in computational effort.

Example 6.2

Solve the following system using LU decomposition

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}= \begin{pmatrix} -1 \\ 10 \\ 22 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

We saw in the Example 6.1 above that the LU decomposition of the coefficient matrix is

\[\begin{split} \begin{align*} L &= \begin{pmatrix} 1 & 0 & 0\\ -2 & 1 & 0\\ 3 & 2 & 1 \end{pmatrix},& U &= \begin{pmatrix} 2 & 3 & 1 \\ 0 & -1 & 2 \\ 0 & 0 & 3 \end{pmatrix}. \end{align*} \end{split}\]

Solving \(L \mathbf{y} = \mathbf{b}\) using forward substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & 0\\ -2 & 1 & 0\\ 3 & 2 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} -1\\ 10 \\ 22 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} y_1 &= -1,\\ y_2 &= 10 + 2y_1 = 10 + 2(-1) = 8,\\ y_3 &= 22 - 3y_1 - 2y_2 = 22 - 3(-1) - 2(8) = 9. \end{align*} \end{split}\]

Solving \(U \mathbf{x} = \mathbf{y}\) using back substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & 3 & 1 \\ 0 & -1 & 2 \\ 0 & 0 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -1 \\ 8 \\ 9 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} x_3 &= \frac{9}{3} = 3,\\ x_2 &= -(8 - 2x_3) = -(8 - 2(3)) = -2, \\ x_1 &= \frac{1}{2}(-1 - 3x_2 - x_3) = \frac{1}{2}(-1 - 3(-2) - 3) = 1. \end{align*} \end{split}\]

So the solution is \(\mathbf{x}=(1, -2, 3)\).

To verify the solution we check that \(A \mathbf{x} = b\)

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix} \begin{pmatrix} 1 \\ -2 \\ 3 \end{pmatrix}= \begin{pmatrix} -1 \\ 10 \\ 22 \end{pmatrix}. \quad \checkmark \end{align*} \end{split}\]