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_{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*} \end{split}\]

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[1]

\[\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}\]

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}\]

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 \vec{0}_{n \times n}\)

  • For each column \(j\) in \(A\) do

    • For each row \(i\) in \(A\) do

      • If \(i \leq j\)

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

      • Else

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

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

Example 6.1

Determine the LU decomposition of the following matrix

\[\begin{split} A = \begin{pmatrix} 1 & 3 & 0\\ 2 & -4 & -1\\ -3 & 1 & 2 \end{pmatrix}.\end{split}\]
Solution (click to show)

Stepping through the columns of \(A\)

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

Therefore

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

Checking that \(LU=A\)

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

6.1.1. Code#

The code below defines a function called lu() which calculates the LU decomposition of a matrix.

def lu(A):
    n = A.shape[0]
    L, U = np.eye(n), np.zeros((n, n))
    for j in range(n):
        for i in range(n):
            sum_ = 0
            if i <= j:
                for k in range(i):
                    sum_ += L[i,k] * U[k,j]
        
                U[i,j] = A[i,j] - sum_   
            
            else:         
                for k in range(j):
                    sum_ += L[i,k] * U[k,j]
                    
                L[i,j] = (A[i,j] - sum_) / U[j,j]
    
    return L, U

function [L, U] = lu(A)

n = size(A, 1);
L = eye(n);
U = zeros(n);

for j = 1 : n
    for i = 1 : n
        sum_ = 0;
        if i <= j
            for k = 1 : i - 1
                sum_ = sum_ + L(i,k) * U(k,j);
            end
            U(i,j) = A(i,j) - sum_;
        else
            for k = 1 : j - 1
                sum_ = sum_ + L(i,k) * U(k,j);
            end
            L(i,j) = (A(i,j) - sum_) / U(j,j);
        end
    end
end

end