6G6Z3017 Computational Methods in ODEs
Matrix decomposition, also known as matrix factorisation, is the process of breaking down a matrix into simpler components that can be more easily analyzed or manipulated.
For example, given the matrix \(\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}\) we can decompose this as the product of two \(2 \times 2\) matrices
\[ \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 3 & 1 \end{pmatrix} \begin{pmatrix} 1 & 2 \\ 0 & -2 \end{pmatrix} \]
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.
The example given on the previous slide is an example of LU decomposition
\[ \begin{align*} \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} &= \begin{pmatrix} 1 & 0 \\ 3 & 1 \end{pmatrix} \begin{pmatrix} 1 & 2 \\ 0 & -2 \end{pmatrix} \\ A \quad \,\, &= \quad\,\, L \qquad\quad U. \end{align*} \]
Consider the LU decomposition of a \(3\times 3\) matrix
\[ \begin{align*} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{32} \\ 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*} \]
which gives a system of 9 equations (one for each element in \(A\)) in 12 unknowns which has an infinite number of solutions.
If we use the condition \(\ell_{ii} = 1\) then \[ \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*} \]
which gives a system of 9 equations in 9 unknowns.
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*} \]
Note
Determine the LU decomposition of the following matrix
\[ A = \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix}.\]
Solution
We need to calculate the values of \(\ell_{ij}\) and \(u_{ij}\) such that
\[ \begin{align*} \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \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} & \ell_{21} u_{12} + u_{22} & \ell_{21} u_{13} + u_{23} \\ \ell_{31} u_{11} & \ell_{31} u_{12} + \ell_{32} u_{22} & \ell_{31} u_{13} + \ell_{32} u_{23} + u_{33} \end{pmatrix} \end{align*} \]
\[ \begin{align*} \begin{pmatrix} 2 & 3 & 1 \\ -4 & -7 & 0 \\ 6 & 7 & 10 \end{pmatrix} = \begin{pmatrix} u_{11} & u_{12} & u_{13} \\ \ell_{21} u_{11} & \ell_{21} u_{12} + u_{22} & \ell_{21} u_{13} + u_{23} \\ \ell_{31} u_{11} & \ell_{31} u_{12} + \ell_{32} u_{22} & \ell_{31} u_{13} + \ell_{32} u_{23} + u_{33} \end{pmatrix} \end{align*} \]
Solving for \(\ell_{ij}\) and \(u_{ij}\) column-by-column
\[ \begin{align*} j &= 1: & 2 &= u_{11} & \therefore u_{11} &= 2, \\ &&-4 &= \ell_{21} u_{11} = 2\ell_{21} & \therefore \ell_{21} &= -2, \\ && 6 &= \ell_{31} u_{11} = 2\ell_{31} & \therefore \ell_{31} &= 3, \\ j &= 2: & 3 &= u_{12} & \therefore u_{12} &= 3, \\ && -7 &= \ell_{21}u_{12} + u_{22} = -6 + u_{22} & \therefore u_{22} &= -1, \\ && 7 &= \ell_{31}u_{12} + \ell_{32}u_{22}= 9 - \ell_{32} & \therefore \ell_{32} &= 2, \\ j &= 3: & 1 &= u_{13} & \therefore u_{13} &= 1, \\ && 0 &= \ell_{21}u_{13} + u_{23} = -2 + u_{23} & \therefore u_{23} &= 2, \\ && 10 &= \ell_{31}u_{13} + \ell_{32}u_{23} + u_{33} = 3 + 4 + u_{33} & \therefore u_{33} &= 3. \end{align*} \] Therefore \[ \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*} \]
Checking that \(LU=A\)
\[ \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}. \qquad \checkmark \end{align*} \]
A system of linear equations with \(m\) equations and \(n\) unknowns is expressed as
\[ \begin{align*} a_{11} x_1 +a_{12} x_2 +\cdots +a_{1n} x_n &=b_1 ,\\ a_{21} x_1 +a_{22} x_2 +\cdots +a_{2n} x_n &=b_2 ,\\ &\vdots \\ a_{m1} x_1 +a_{m2} x_2 +\cdots +a_{mn} x_n &=b_n , \end{align*} \]
where \(x_i\) are the unknown variables, \(a_{ij}\) are coefficients and \(b_i\) are constant terms. It is often more convenient to express a system of linear equations as a matrix equation
\[ \begin{align*} A \mathbf{x} = \mathbf{b}, \end{align*} \]
where \(A\) is the coefficient matrix, \(\mathbf{x}\) is the variable vector and \(b\) is the constant vector
The solution (if it exists) is the vector \(\mathbf{x}\) for which the equation \(A\mathbf{x} = \mathbf{b}\) is satisfied.
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
\[ LU \mathbf{x} = \mathbf{b}.\]
Let \(\mathbf{y} = U \mathbf{x}\) then
\[ L \mathbf{y} = \mathbf{b}. \]
\(L\) is lower-triangular then we can solve \(L \mathbf{y} = \mathbf{b}\) using forward-substitution
\(U\) is upper-triangular then we can solve \(U \mathbf{x} = \mathbf{y}\) using back-substitution
The advantage of using LU decomposition is that once \(L\) and \(U\) are calculated, it can be used for any values of \(\mathbf{b}\), e.g, Newton’s method for calculating IRK stage values.
Solve the following system using LU decomposition
\[ \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*} \]
Solution
We saw in the previous example that the LU decomposition of the coefficient matrix is
\[ \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*} \]
Solving \(L \mathbf{y} = \mathbf{b}\) using forward substitution
\[ \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*} \]
gives
\[ \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*} \]
Solving \(U \mathbf{x} = \mathbf{y}\) using back substitution
\[ \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*} \]
gives
\[ \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*} \]
So the solution is \(\mathbf{x}=(1, -2, 3)\).
If \(u_{jj}\) is zero or some small number it will mean that LU decomposition is prone to computational rounding errors
For example, consider the solution of the following system using LU decomposition
\[ \begin{align*} \begin{pmatrix} 0.0001 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 1 \\ 2 \end{pmatrix}. \end{align*}\]
Calculating the LU decomposition of the coefficient matrix we have
\[ \begin{align*} L &= \begin{pmatrix} 1 & 0 \\ 10000 & 1 \end{pmatrix}, & U &= \begin{pmatrix} 1 & 1 \\ 0 & -9999 \end{pmatrix}, \end{align*} \]
and solving using forward and back substitution and working to four decimal place accuracy
\[ \begin{align*} \begin{pmatrix} 1 & 0 \\ 10000 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &= \begin{pmatrix} 1 \\ 2 \end{pmatrix} & \therefore y_1 &= 1, y_2 = -9998, \\ \begin{pmatrix} 1 & 1 \\ 0 & -9999 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 1 \\ -9998 \end{pmatrix} & \therefore x_1 &= 0, x_2 = 0.9998\dot{9} \approx 1. \end{align*} \]
Checking the solution by substituting into the original system gives
\[ \begin{align*} \mathbf{b} &= \begin{pmatrix} 0.0001 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \neq \begin{pmatrix} 1 \\ 2 \end{pmatrix}. \end{align*} \]
However, if we swap the rows of the system we have
\[ \begin{align*} \begin{pmatrix} 1 & 1 \\ 0.0001 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 2 \\ 1 \end{pmatrix}. \end{align*}\]
which has the LU decomposition
\[ \begin{align*} L &= \begin{pmatrix} 1 & 0 \\ 0.0001 & 1 \end{pmatrix}, & U &= \begin{pmatrix} 1 & 1 \\ 0 & 0.9999 \end{pmatrix}. \end{align*} \]
Solving using forward and back substitution with the same decimal place accuracy used before
\[ \begin{align*} \begin{pmatrix} 1 & 0 \\ 0.0001 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} &= \begin{pmatrix} 2 \\ 1 \end{pmatrix} & \therefore y_1 &= 2, y_2 = 0.9998, \\ \begin{pmatrix} 1 & 1 \\ 0 & 0.9999 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} &= \begin{pmatrix} 2 \\ 0.9998 \end{pmatrix} & \therefore x_1 &= 1, x_2 = 0.9998\dot{9} \approx 1. \end{align*} \]
Checking the solution by substituting into the original system gives
\[ \begin{align*} \mathbf{b} &= \begin{pmatrix} 0.0001 & 1 \\ 1 & 1 \\ \end{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1.0001 \\ 2 \end{pmatrix} \approx \begin{pmatrix} 1 \\ 2 \end{pmatrix}. \qquad \checkmark \end{align*} \]
The permutations applied to the coefficient matrix are recorded in a permutation matrix \(P\) which is determined by applying the same permutations to the identity matrix.
For example, applying partial pivoting the matrix
\[ \begin{align*} &\begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} \begin{array}{l} R_1 \leftrightarrow R_3 \\ \phantom{x} \\ \phantom{x} \end{array} \longrightarrow \begin{pmatrix} 3 & -2 & 2\\ 1 & 0 & 2\\ 0 & 1 & -2 \end{pmatrix} \begin{array}{l} \phantom{x} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array} & \longrightarrow &\begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}. \end{align*} \]
Apply the same row operations to the identity matrix
\[ \begin{align*} &\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{array}{l} R_1 \leftrightarrow R_3 \\ \phantom{x} \\ \phantom{x} \end{array} \longrightarrow \begin{pmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{pmatrix} \begin{array}{l} \phantom{x} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{array} & \longrightarrow &\begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix} = P. \end{align*} \]
Note that \(A\) by \(P\) on the left gives the matrix \(A\) after partial pivoting has been applied
\[ \begin{align*} \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} = \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}. \end{align*} \]
To implement LU decomposition with partial pivoting (LUP decomposition) we apply partial pivoting to determine a permutation matrix \(P\) before calculating the LU decomposition of \(PA\), i.e.,
\[ LU = PA. \]
Since we have applied row swaps to the coefficient matrix we must also apply the same row swaps to the constant vector. So to solve the system \(A \mathbf{x} = \mathbf{b}\) using LUP decomposition we have
\[ \begin{align*} PA \mathbf{x} &= P \mathbf{b} \\ \therefore LU \mathbf{x} &= P \mathbf{b}, \end{align*} \]
and we solve \(L \mathbf{y} = P \mathbf{b}\) for \(\mathbf{y}\) using forward substitution and \(U \mathbf{x} = \mathbf{y}\) for \(\mathbf{x}\) using back substitution.
Solve the following system of linear equations using LU decomposition with partial pivoting
\[ \begin{align*} \begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 10 \\ -4 \\ -8 \end{pmatrix}. \end{align*} \]
Solution
We have seen from the previous example that applying partial pivoting to the coefficient matrix results in
\[ \begin{align*} P &= \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix}, \\ \therefore PA &= \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & -2\\ 1 & 0 & 2\\ 3 & -2 & 2 \end{pmatrix} = \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix}, & \end{align*} \]
Calculating the LU decomposition of \(PA\)
\[ \begin{align*} \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 1 & 0 & 2 \end{pmatrix} = \begin{pmatrix} u_{11} & u_{12} & u_{13} \\ \ell_{21} u_{11} & \ell_{21} u_{12} + u_{22} & \ell_{21} u_{13} + u_{23} \\ \ell_{31} u_{11} & \ell_{31} u_{12} + \ell_{31} u_{22} & \ell_{31} u_{13} + \ell_{32} u_{23} + u_{33} \end{pmatrix} \end{align*} \]
\[ \begin{align*} j &= 1: & 3 &= u_{11} & \therefore u_{11} &= 3, \\ && 0 &= \ell_{21} u_{11} = 3\ell_{21} & \therefore \ell_{21} &= 0, \\ && 1 &= \ell_{31} u_{11} = 3\ell_{31} & \therefore \ell_{31} &= \frac{1}{3}, \\ j &= 2: & -2 &= u_{12} & \therefore u_{12} &= -2, \\ && 1 &= \ell_{21} u_{12} + u_{22} = 0 + u_{22} & \therefore u_{22} &= 1, \\ && 0 &= \ell_{31} u_{12} + \ell_{32} u_{22} = -\frac{2}{3} + \ell_{32} & \therefore \ell_{32} &= \frac{2}{3}, \\ j &= 3: & 2 &= u_{13} & \therefore u_{13} &= 2, \\ && -2 &= \ell_{21} u_{13} + u_{23} = 0 + u_{23} & \therefore u_{23} &= -2, \\ && 2 &= \ell_{31} u_{13} + \ell_{32} u_{23} + u_{33} = \frac{2}{3} - \frac{4}{3} + u_{33} & \therefore u_{33} &= \frac{8}{3} \end{align*} \]
Therefore
\[ \begin{align*} L &= \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ \frac{1}{3} & \frac{2}{3} & 1 \end{pmatrix}, & U &= \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 0 & 0 & \frac{8}{3} \end{pmatrix}. \end{align*} \]
Solving \(L \mathbf{y} = P \mathbf{b}\) using forward substitution
\[ \begin{align*} \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ \frac{1}{3} & \frac{2}{3} & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 10 \\ -4 \\ -8 \end{pmatrix} = \begin{pmatrix} -8 \\ 10 \\ -4 \end{pmatrix}, \end{align*} \]
gives
\[ \begin{align*} y_1 &= -8, \\ y_2 &= 10 - 0(8)=10, \\ y_3 &= -4 - \frac{1}{3}(8) + \frac{2}{3}(10) = -8. \end{align*} \]
Solving \(U \mathbf{x} = \mathbf{y}\) using back substitution
\[ \begin{align*} \begin{pmatrix} 3 & -2 & 2\\ 0 & 1 & -2\\ 0 & 0 & 8/3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -8 \\ 10 \\ -8 \end{pmatrix}, \end{align*} \]
gives
\[ \begin{align*} x_3 &= \frac{3}{8}(-8) = -3, \\ x_2 &= \frac{1}{1}(10+2(-3)) = 4, \\ x_1 &= \frac{1}{3}(-8+2(4)-2(-3)) = 2. \end{align*} \]
A square matrix \(A\) is said to be positive definite if
\[ \mathbf{x}^\mathsf{T} A \mathbf{x} > 0, \]
for all \(\mathbf{x} \in \mathbb{R}^n \backslash \{\mathbf{0}\}\).
A square matrix \(A\) is positive definite if it is symmetric and the determinants of all \(k \times k\) upper-left sub-matrices are all positive.
Consider the matrix
\[ \begin{align*} A = \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix}. \end{align*} \]
Matrix \(A\) is symmetric since \(A=A^\mathsf{T}\). Checking the determinants of the upper-left sub-matrices
\[ \begin{align*} \det(2) &= 2 > 0, \\ \det\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} &= 4 - 1 = 3 > 0, \\ \det\begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix} &= 6 - 2 + 0 = 4 > 0 . \end{align*} \]
Since all of the upper-left square submatrices have a positive determinant then \(A\) is a positive-definite matrix.
Given a positive definite matrix \(A\) the Cholesky decomposition decomposes \(A\) into the product of a lower triangular matrix \(L\) and its transpose \(L^\mathsf{T}\), i.e.,
\[ \begin{align*} A = LL^\mathsf{T}. \end{align*} \]
Consider the Cholesky decomposition of a \(3\times 3\) matrix
\[ \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} \ell_{11} & \ell_{21} & \ell_{31} \\ 0 & \ell_{22} & \ell_{32} \\ 0 & 0 & \ell_{33} \end{pmatrix}\\ &= \begin{pmatrix} \ell_{11}^2 & \ell_{11} \ell_{21} & \ell_{11} \ell_{31} \\ \ell_{11} \ell_{21} & \ell_{21}^2 +\ell_{22}^2 & \ell_{21} \ell_{31} +\ell_{22} \ell_{33} \\ \ell_{11} \ell_{31} & \ell_{21} \ell_{31} +\ell_{22} \ell_{33} & \ell_{31}^2 +\ell_{32}^2 +\ell_{33}^2 \end{pmatrix}. \end{align*} \]
The elements on the main diagonal are
\[ \begin{align*} a_{jj} =\ell_{jj}^2 +\sum_{k=1}^{j-1} \ell_{jk}^2 , \end{align*} \]
and the other elements are
\[ \begin{align*} a_{ij} =\sum_{k=1}^j \ell_{ik} \ell_{jk} = \ell_{jj} \ell_{ij} +\sum_{k=1}^{j-1} \ell_{ik} \ell_{jk}. \end{align*} \]
Rearranging these expressions gives
\[ \begin{align} \ell_{jj} &= \left( a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2 \right)^{\frac{1}{2}}, & j &= 1, \ldots, n,\\ \ell_{ij} &= \dfrac{1}{\ell_{jj} }\left(a_{ij} -\displaystyle \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), & i &= j + 1,\ldots ,n. \end{align} \]
Calculate the Cholesky decomposition of the following matrix
\[ \begin{align*} A = \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix}. \end{align*} \]
Solution
We need to solve \(A = LL^\mathsf{T}\) for \(\ell_{ij}\), i.e.,
\[ \begin{align*} \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix} &= \begin{pmatrix} \ell_{11} & 0 & 0 \\ \ell_{21} & \ell_{22} & 0 \\ \ell_{31} & \ell_{32} & \ell_{33} \end{pmatrix} \begin{pmatrix} \ell_{11} & \ell_{21} & \ell_{31} \\ 0 & \ell_{22} & \ell_{32} \\ 0 & 0 & \ell_{33} \end{pmatrix} \\ &= \begin{pmatrix} \ell_{11}^2 & \ell_{11} \ell_{21} & \ell_{11} \ell_{31} \\ \ell_{21} \ell_{11} & \ell_{21}^2 + \ell_{22}^2 & \ell_{21} \ell_{31} + \ell_{22} \ell_{32} \\ \ell_{31}\ell_{11} & \ell_{31} \ell_{21} + \ell_{32} \ell_{22} & \ell_{31}^2 + \ell_{32}^2 + \ell_{33}^2 \end{pmatrix}. \end{align*} \]
Note that \(LL^\mathsf{T}\) is symmetric so we only need to calculate the lower triangular elements.
\[ \begin{align*} \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix} = \begin{pmatrix} \ell_{11}^2 & \ell_{11} \ell_{21} & \ell_{11} \ell_{31} \\ \ell_{21} \ell_{11} & \ell_{21}^2 + \ell_{22}^2 & \ell_{21} \ell_{31} + \ell_{22} \ell_{32} \\ \ell_{31}\ell_{11} & \ell_{31} \ell_{21} + \ell_{32} \ell_{22} & \ell_{31}^2 + \ell_{32}^2 + \ell_{33}^2 \end{pmatrix}. \end{align*} \]
\[ \begin{align*} j &= 1: & 4 &= \ell_{11}^2 & \therefore \ell_{11} &= 2, \\ && -2 &= \ell_{21} \ell_{11} = 2\ell_{21} & \therefore \ell_{21} &= -1, \\ && -4 &= \ell_{31} \ell_{11} = 2\ell_{31} & \therefore \ell_{31} &= -2, \\ j &= 2: & 10 &= \ell_{21}^2 + \ell_{22}^2 = 1 + \ell_{22}^2 & \therefore \ell_{22} &= 3, \\ && 5 &= \ell_{31} \ell_{21} + \ell_{32} \ell_{22} = 2 + 3\ell_{32} & \therefore \ell_{32} &= 1, \\ j &= 3: & 14 &= \ell_{31}^2 + \ell_{32}^2 + \ell_{33}^2 = 4 + 1 + \ell_{33}^2 & \therefore \ell_{33} &= 3, \end{align*} \]
therefore \(L=\begin{pmatrix} 2 & 0 & 0 \\ -1 & 3 & 0 \\ -2 & 1 & 3 \end{pmatrix}\).
Checking that \(A = LL^\mathsf{T}\)
\[ \begin{align*} \begin{pmatrix} 2 & 0 & 0\\ -1 & 3 & 0\\ -2 & 1 & 3 \end{pmatrix} \begin{pmatrix} 2 & -1 & -2\\ 0 & 3 & 1\\ 0 & 0 & 3 \end{pmatrix} = \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix}. \end{align*} \]
Given a system of linear equations of the form \(A \mathbf{x} = \mathbf{b}\) where Cholesky decomposition has been applied to the coefficient matrix then since \(A = LL^\mathsf{T}\) we have
\[ LL^\mathsf{T} \mathbf{x} = \mathbf{b}.\]
Let \(\mathbf{y} = L^\mathsf{T} \mathbf{x}\) then
\[ L \mathbf{y} = \mathbf{b}. \]
Similar to solving systems using LU decomposition we solve \(L \mathbf{y} = \mathbf{b}\) for \(\mathbf{y}\) using forward substitution and then \(L^\mathsf{T} \mathbf{x} = \mathbf{y}\) using back substitution.
Solve the following system of linear equations using the Cholesky decomposition.
\[ \begin{align*} \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -2 \\ 49 \\ 27 \end{pmatrix}. \end{align*} \]
Solution
We saw in the previous example that the Cholesky decomposition of the matrix \(A\) is
\[ L=\begin{pmatrix} 2 & 0 & 0 \\ -1 & 3 & 0 \\ -2 & 1 & 3 \end{pmatrix}. \]
Solving \(L\mathbf{y}=\mathbf{b}\) using forward substitution
\[ \begin{align*} \begin{pmatrix} 2 & 0 & 0\\ -1 & 3 & 0\\ -2 & 1 & 3 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} -2 \\ 49 \\ 27 \end{pmatrix} && \therefore \mathbf{y} = \begin{pmatrix} -1 \\ 16 \\ 3 \end{pmatrix}. \end{align*} \]
Solving \(L^\mathsf{T} \mathbf{x}=\mathbf{y}\) using back substitution
\[ \begin{align*} \begin{pmatrix} 2 & -1 & -2\\ 0 & 3 & 1\\ 0 & 0 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -1 \\ 16 \\ 3 \end{pmatrix}, && \therefore \mathbf{x} = \begin{pmatrix} 3 \\ 5\\ 1 \end{pmatrix}. \end{align*} \]