## 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 is such that the matrix is an upper triangular matrix, then the system can be solved directly using backward substitution. Similarly, if a linear system of equations is such that the matrix 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 defined as:

which in matrix form has the form:

The matrix is a lower triangular matrix. The system can be solved by forward substitution. Namely, the first equation gives . The second equation gives: . The third equation gives: . Forward substitution can be programmed in a very fast and efficient manner. Consider the following system:

Therefore:

This can be written in the following simple form with starting at to :

#### Solving Upper Triangular Systems

Similarly, consider the linear system defined as:

which in matrix form has the form:

The matrix is an upper triangular matrix. The system can be solved by backward substitution. Namely, the third equation gives . The second equation gives: . The first equation gives: . Backward substitution can be easily programmed in a very fast and efficient manner. Consider the following system:

Therefore:

This can be written in the following simple form with starting at and decreasing to :

### Example

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

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

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

The above system can be solved quickly as follows. In the equation , set , then . Since is a lower triangular matrix, we can solve for easily. Then, since and is an upper triangular matrix, we can solve for easily. The equation has the form:

The above system can be solved using forward substitution. Therefore . The second equation yields: . The third equation yields: . We can now solve as follows:

This system can be solved using backward substitution. The third equation yields . The second equation yields . The third equation yields . This example implies that if the matrix can be decomposed such that where is a lower triangular matrix and 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 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 be a linear system of equations. Let be an LU decomposition for with being a lower triangular matrix and being an upper triangular matrix. Then: . Set , then . The elements of the vector can be solved using forward substitution as explained above. Then, the system can be solved using backward substitution as explained above.

#### LU Decomposition Algorithm

There are various algorithms to find the decomposition . 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 is turned into an upper triangular matrix which is the matrix . The matrix has values of 1 in the diagonal entries. The entries that are below the diagonals are such that negative the factor that was used to reduce row using the pivot row during the Gauss elimination step. The following example illustrates the process. Consider the matrix

To find both and we will write the identity matrix on the left and the matrix on the right:

Gauss elimination is then used to eliminate the coefficient by multiplying row 1 by and then adding it to row 2. Then, the coefficient

Gauss elimination is then used to eliminate the coefficient by multiplying row 1 by and then adding it to row 3. Then, the coefficient

Gauss elimination is then used to eliminate the coefficient by multiplying row 2 by and then adding it to row 3. Then, the coefficient

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

### Programming the Method

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

View Mathematica CodeLU[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]

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 without adding the vector . 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 . So, if multiple solutions corresponding to multiple possible values of the vector are sought, then this method could be faster than the Gauss elimination methods.