Engineering at Alberta Courses

Direct Methods: LU Decomposition

Introduction

Another method that is comparable in efficiency and speed to the Gauss elimination methods stated above is the LU decomposition. It was implied when the Gauss elimination method was introduced that if a linear system of equations $Ax=b$ is such that the matrix $A$ is an upper triangular matrix, then the system can be solved directly using backward substitution. Similarly, if a linear system of equations $Ax=b$ is such that the matrix $A$ is a lower triangular matrix, then the system can be solved directly using forward substitution. The following two examples illustrate this.

Solving Lower Triangular Systems

Consider the linear system $Ax=b$ defined as:

$\begin{split} 5x_1+0+0&=15\\ 2x_1+3x_2+0&=6\\ 4x_1+3x_2+2x_3&=2\\ \end{split}$

which in matrix form has the form:

$\left(\begin{matrix}5&0&0\\2&3&0\\4&3&2\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}15\\6\\2\end{array}\right)$

The matrix $A$ is a lower triangular matrix. The system can be solved by forward substitution. Namely, the first equation gives $x_1=3$ . The second equation gives: $2(3)+3x_2=6\Rightarrow x_2=0$ . The third equation gives: $4(3)+0+2x_3=2\Rightarrow x_3=-5$ . Forward substitution can be programmed in a very fast and efficient manner. Consider the following system:

$\left(\begin{matrix}L_{11}&0&0&\cdots&0\\L_{21}&L_{22}&0&\cdots&0\\L_{31}&L_{32}&L_{33}&\cdots&0\\\vdots&\vdots&\vdots&&\vdots\\L_{n1}&L_{n2}&L_{n3}&\cdots&L_{nn}\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\\\vdots\\x_n\end{array}\right)= \left(\begin{array}{c}b_1\\b_2\\b_3\\\vdots\\b_n\end{array}\right)$

Therefore:

$\begin{split} x_1&=\frac{b_1}{L_{11}}\\ x_2&=\frac{b_2-L_{21}x_1}{L_{22}}\\ x_3&=\frac{b_3-L_{31}x_1-L_{32}x_2}{L_{33}}\\ &\vdots\\ x_n&=\frac{b_n-\sum_{j=1}^{n-1}L_{nj}x_j}{L_{nn}} \end{split}$

This can be written in the following simple form with $i$ starting at $i=1$ to $n$ :

$x_i=\frac{b_i - \sum_{j=1}^{i-1}L_{ij}x_j}{L_{ii}}$

Solving Upper Triangular Systems

Similarly, consider the linear system $Ax=b$ defined as:

$\begin{split} 2x_1+3x_2+4x_3&=2\\ 0+3x_2+2x_3&=6\\ 0+0+5x_3&=15 \end{split}$

which in matrix form has the form:

$\left(\begin{matrix}2&3&4\\0&3&2\\0&0&5\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}2\\6\\15\end{array}\right)$

The matrix $A$ is an upper triangular matrix. The system can be solved by backward substitution. Namely, the third equation gives $x_3=3$ . The second equation gives: $3x_2+2(3)=6\Rightarrow x_2=0$ . The first equation gives: $2x_1+0+4(3)=2\Rightarrow x_1=-5$ . Backward substitution can be easily programmed in a very fast and efficient manner. Consider the following system:

$\left(\begin{matrix}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&&\vdots\\0&0&0&\cdots&U_{nn}\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\\\vdots\\x_n\end{array}\right)= \left(\begin{array}{c}b_1\\b_2\\b_3\\\vdots\\b_n\end{array}\right)$

Therefore:

$\begin{split} x_n&=\frac{b_n}{U_{nn}}\\ &\vdots\\ x_{3}&=\frac{b_3-\sum_{j=4}^nU_{3j}x_j}{U_{33}}\\ x_{2}&=\frac{b_2-\sum_{j=3}^nU_{2j}x_j}{U_{22}}\\ x_1&=\frac{b_1-\sum_{j=2}^{n}U_{1j}x_j}{U_{11}} \end{split}$

This can be written in the following simple form with $i$ starting at $i=n$ and decreasing to $1$ :

$x_i=\frac{b_i - \sum_{j=i+1}^{n}U_{ij}x_j}{U_{ii}}$

Example

We will first introduce the LU decomposition using an example, consider the system $Ax=b$ defined as:

$\left(\begin{matrix}2&1&1\\5&2&2\\4&3&2\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}5\\6\\3\end{array}\right)$

Without explaining how at this point, we can decompose the matrix $A$ as the multiplication of two matrices $A=LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. For this particular example, $L$ and $U$ have the forms:

$L=\left(\begin{matrix}1&0&0\\2.5&1&0\\2&-2&1\end{matrix}\right)\qquad U=\left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\0&0&-1\end{matrix}\right)$

The reader should verify that indeed $A=LU$ . The linear system of equations can then be written as:

$\left(\begin{matrix}1&0&0\\2.5&1&0\\2&-2&1\end{matrix}\right)\left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\0&0&-1\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}5\\6\\3\end{array}\right)$

The above system can be solved quickly as follows. In the equation $LUx=b$ , set $y = Ux$ , then $Ly=b$ . Since $L$ is a lower triangular matrix, we can solve for $y$ easily. Then, since $Ux=y$ and $U$ is an upper triangular matrix, we can solve for $x$ easily. The equation $Ly=b$ has the form:

$\left(\begin{matrix}1&0&0\\2.5&1&0\\2&-2&1\end{matrix}\right)\left(\begin{array}{c}y_1\\y_2\\y_3\end{array}\right)=\left(\begin{array}{c}5\\6\\3\end{array}\right)$

The above system can be solved using forward substitution. Therefore $y_1=5$ . The second equation yields: $y_2=6-2.5(5)=-6.5$ . The third equation yields: $y_3=3-2(5)+2(-6.5)=-20$ . We can now solve $Ux=y$ as follows:

$\left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\0&0&-1\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}y_1\\y_2\\y_3\end{array}\right)=\left(\begin{array}{c}5\\-6.5\\-20\end{array}\right)$

This system can be solved using backward substitution. The third equation yields $x_3=20$ . The second equation yields $-0.5x_2-0.5(20)=-6.5\Rightarrow x_2=-7$ . The third equation yields $2x_1-7+20=5\Rightarrow x_1=-4$ . This example implies that if the matrix $A$ can be decomposed such that $A=LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix, then the system can be easily solved using forward and backward substitution. It turns out that such decomposition always exists provided that the rows of the matrix $A$ are properly ordered (by that we mean there are no zero pivots during the Gauss elimination step).

Formal Representation

The following is a more formal representation of the method. Let $Ax=b$ be a linear system of equations. Let $A=LU$ be an LU decomposition for $A$ with $L$ being a lower triangular matrix and $U$ being an upper triangular matrix. Then: $Ax=LUx=L(Ux)$ . Set $y = Ux$ , then $Ly=b$ . The elements of the vector $y$ can be solved using forward substitution as explained above. Then, the system $Ux=y$ can be solved using backward substitution as explained above.

LU Decomposition Algorithm

There are various algorithms to find the decomposition $LU$ . Mathematica has a built-in function to find the LU decomposition called “LUDecomposition”. Here we present the most common algorithm. In the first step, using Gauss elimination, the matrix $A$ is turned into an upper triangular matrix which is the matrix $U$ . The matrix $L$ has values of 1 in the diagonal entries. The entries $L_{ij}$ that are below the diagonals are such that $L_{ij}=$ negative the factor that was used to reduce row $i$ using the pivot row $j$ during the Gauss elimination step. The following example illustrates the process. Consider the matrix

$\left(\begin{matrix}2&1&1\\5&2&2\\4&3&2\end{matrix}\right)$

To find both $L$ and $U$ we will write the identity matrix on the left and the $A$ matrix on the right:

$\left(\begin{matrix}1&0&0\\0&1&0\\0&0&1\end{matrix}\right) \qquad \left(\begin{matrix}2&1&1\\5&2&2\\4&3&2\end{matrix}\right)$

Gauss elimination is then used to eliminate the coefficient $A_{21}$ by multiplying row 1 by $\frac{-5}{2}$ and then adding it to row 2. Then, the coefficient $L_{21}=-\frac{-5}{2}=2.5$

$\left(\begin{matrix}1&0&0\\2.5&1&0\\0&0&1\end{matrix}\right) \qquad \left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\4&3&2\end{matrix}\right)$

Gauss elimination is then used to eliminate the coefficient $A_{31}$ by multiplying row 1 by $\frac{-4}{2}$ and then adding it to row 3. Then, the coefficient $L_{31}=-\frac{-4}{2}=2$

$\left(\begin{matrix}1&0&0\\2.5&1&0\\2&0&1\end{matrix}\right) \qquad \left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\0&1&0\end{matrix}\right)$

Gauss elimination is then used to eliminate the coefficient $A_{32}$ by multiplying row 2 by $\frac{1}{0.5}$ and then adding it to row 3. Then, the coefficient $L_{32}=-\frac{1}{0.5}=-2$

$\left(\begin{matrix}1&0&0\\2.5&1&0\\2&-2&1\end{matrix}\right) \qquad \left(\begin{matrix}2&1&1\\0&-0.5&-0.5\\0&0&-1\end{matrix}\right)$

The matrix on the left hand side is $L$ and the one on the right hand side is $U$ .

Programming the Method

Programming the LU decomposition method requires three steps. In the first one, the Gauss elimination is used to create the $L$ and $U$ matrices. Then, forward substitution is used to solve for the vector $y$ in $Ly=b$ . Then, backward substitution is used to solve for the vector $x$ in $Ux=y$ . The following code has two procedures, the first one takes the matrix $A$ and produces the two matrices $L$ and $U$ . The second procedure utilizes the first procedure to produce $L$ and $U$ and then uses the backward and forward substitutions to find the vectors $y$ and then $x$ .

View Mathematica Code

LU[A_] :=
 (n = Length[A];
  L = IdentityMatrix[n];
  U = A;
  For[k = 1, k <= n, k++, 
   Do[L[[i, k]] = U[[i, k]]/U[[k, k]];U[[i]] = U[[i]] - U[[i, k]]/U[[k, k]]*U[[k]], {i, k + 1, n}]];
  {L, U})

LUAB[A_, b_] := 
(n = Length[A];
  L = LU[A][[1]];
  U = LU[A][[2]];
  y = Table[0, {i, 1, n}];
  x = Table[0, {i, 1, n}];
  For[i = 1, i <= n,i = i + 1, (BB = 0; Do[BB = L[[i, l]]*y[[l]] + BB, {l, 1, i - 1}]; y[[i]] = (b[[i]] - BB)/L[[i, i]])];
  For[i = n, i >= 1, i = i - 1, (BB = 0; Do[BB = U[[i, l]]*x[[l]] + BB, {l, i + 1, n}]; x[[i]] = (y[[i]] - BB)/U[[i, i]])];
  x)

A = {{1, 2, 3, 5}, {2, 0, 1, 4}, {1, 2, 2, 5}, {4, 3, 2, 2}};
b = {-4, 8, 0, 10};
LUAB[A, b]

View Python Code

import numpy as np
def LU(A):
  n = len(A)
  L = np.identity(n)
  U = np.array(A).astype(np.float)
  for k in range(n):
    for i in range(k + 1, n):
      L[i,k] = U[i,k]/U[k,k]
      U[i] = U[i] - U[i,k]/U[k,k]*U[k]
  return [L, U]

def LUAB(A,b):
  n = len(A)
  L = LU(A)[0]
  U = LU(A)[1]
  y = np.zeros(n)
  x = np.zeros(n)
  for i in range(n):
    BB = 0
    for l in range(0, i):
      BB = L[i,l]*y[l] + BB 
    y[i] = (b[i] - BB)/L[i,i]
  for i in range(n-1,-1,-1):
    BB = 0
    for l in range(i,n):
      BB = U[i,l]*x[l] + BB 
      x[i] = (y[i] - BB)/U[i,i]      
  return x

A = [[1, 2, 3, 5], [2, 0, 1, 4], [1, 2, 2, 5], [4, 3, 2, 2]]
b = [-4, 8, 0, 10]
LUAB(A, b)

We are providing two MATLAB files; the first calculates the decomposition, while the second solves a linear system of equation using the decomosition.

MATLAB file to calculate the LU decomposition Download

MATLAB file to solve a linear system using the LU decomposition Download

Comparison with Previous Direct Methods

In comparison to the naive Gauss elimination method, this method reduces the computations in the forward elimination step by conducting the forward elimination step on the matrix $A$ without adding the vector $b$ . However, it adds another step which is the forward substitution step described above. The advantage of this method is that no manipulations are done on the vector $b$ . So, if multiple solutions corresponding to multiple possible values of the vector $b$ are sought, then this method could be faster than the Gauss elimination methods.

Open Educational Resources