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
.
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]
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.