Open Educational Resources

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.

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.

Lecture video



Leave a Reply

Your email address will not be published. Required fields are marked *