LU Decomposition

Section 2.5 LU Decomposition

In many applications where linear systems appear, one needs to solve \(A \xb = \bb\) for many different vectors \(\bb\text{.}\) For instance, a structure must be tested under several different loads, not just one. As in the example of a truss (2.2.2), the loading in such a problem is usually represented by the vector \(\bb\text{.}\) Gaussian elimination with pivoting is the most efficient and accurate way to solve a linear system. Most of the work in this method is spent on the matrix \(A\) itself. If we need to solve several different systems with the same \(A\text{,}\) and \(A\) is big, then we would like to avoid repeating the steps of Gaussian elimination on \(A\) for every different \(\bb\text{.}\) This can be accomplished by the LU decomposition, which in effect records the steps of Gaussian elimination.

🔗

Subsection 2.5.1 LU decomposition

The main idea of the LU decomposition is to record the steps used in Gaussian elimination on A in the places where the zero is produced. Consider the matrix

\begin{equation*} A = \left(\begin{array}{ccc}1 & -2 & 3 \\ 2 & -5 & 12 \\ 0 & 2 & -10\end{array}\right)\text{.} \end{equation*}

The first step of Gaussian elimination is to subtract \(2\) times the first row from the second row. In order to record what we have done, we will put the multiplier, \(2\text{,}\) into the place it was used to make a zero, i.e. the second row, first column. In order to make it clear that it is a record of the step and not an element of \(A\text{,}\) we will put it in parentheses. This leads to

\begin{equation*} \left(\begin{array}{ccc}1 & -2 & 3 \\ (2) & -1 & 6 \\ 0 & 2 & -10\end{array}\right)\text{.} \end{equation*}

There is already a zero in the lower left corner, so we don’t need to eliminate anything there. We record this fact with a \((0)\text{.}\) To eliminate the third row, second column, we need to subtract \(-2\) times the second row from the third row. Recording the \(-2\) in the spot it was used we have

\begin{equation*} \left(\begin{array}{ccc}1 & -2 & 3 \\ (2) & -1 & 6 \\ (0) & (-2) & 2\end{array}\right)\text{.} \end{equation*}

Let \(U\) be the upper triangular matrix produced, and let \(L\) be the lower triangular matrix with the records and ones on the diagonal, i.e.

\begin{equation*} L = \left(\begin{array}{ccc}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & -2 & 1\end{array}\right) \quad\text{and}\quad U = \left(\begin{array}{ccc}1 & -2 & 3 \\ 0 & -1 & 6 \\ 0 & 0 & 2\end{array}\right)\text{,} \end{equation*}

then we have the following wonderful property:

\begin{equation*} LU = \left(\begin{array}{ccc}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & -2 & 1\end{array}\right) \left(\begin{array}{ccc}1 & -2 & 3 \\ 0 & -1 & 6 \\ 0 & 0 & 2\end{array}\right) = \left(\begin{array}{ccc}1 & -2 & 3 \\ 2 & -5 & 12 \\ 0 & 2 & -10\end{array}\right) = A\text{.} \end{equation*}

Thus we see that \(A\) is actually the product of \(L\) and \(U\text{.}\) Here \(L\) is lower triangular and \(U\) is upper triangular. When a matrix can be written as a product of simpler matrices, we call that a decomposition of \(A\) and this one we call the LU decomposition.

🔗

Subsection 2.5.2 Using LU to solve equations

If we also include pivoting, then an LU decomposition for \(A\) consists of three matrices \(P\text{,}\) \(L\) and \(U\) such that

\begin{equation*} PA = LU\text{.} \end{equation*}

The pivot matrix \(P\) is the identity matrix, with the same rows switched as the rows of \(A\) are switched in the pivoting. For instance,

\begin{equation*} P = \left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0\end{array}\right)\text{,} \end{equation*}

would be the pivot matrix if the second and third rows of \(A\) are switched by pivoting. MATLAB will produce an \(LU\) decomposition with pivoting for a matrix \(A\) with the command [L U P] = lu(A) where P is the pivot matrix. To use this information to solve \(A \xb = \bb\) we first pivot both sides by multiplying by the pivot matrix:

\begin{equation*} PA \xb = P \bb \equiv \mathbf{d}\text{.} \end{equation*}

Substituting \(LU\) for \(PA\) we get

\begin{equation*} LU \xb = \mathbf{d}. \end{equation*}

Then we define \(\yb=U \xb\text{,}\) which is unknown since \(\xb\) is unknown. Using forward-substitution, we can (easily) solve

\begin{equation*} L \yb = \mathbf{d} \end{equation*}

for \(\yb\) and then using back-substitution we can (easily) solve

\begin{equation*} U \xb = \yb \end{equation*}

for \(\xb\text{.}\) In MATLAB this would work as follows:

>> A = rand(5,5)
>> [L U P] = lu(A)
>> b = rand(5,1)
>> d = P*b
>> y = L \ d
>> x = U \ y
>> rnorm = norm(A*x - b)   % Check the result

We can then solve for any other \(\bb\) without redoing the LU step. Repeat the sequence for a new right hand side: c = randn(5,1); you can start at the third line. While this may not seem like a big savings, it would be if \(A\) were a large matrix from an actual application.

🔗

The LU decomposition is an example of Matrix Decomposition which means taking a general matrix \(A\) and breaking it down into components with simpler properties. Here \(L\) and \(U\) are simpler because they are lower and upper triangular. There are many other matrix decompositions that are useful in various contexts. Some of the most useful of these are the QR decomposition, the Singular Value decomposition and Cholesky decomposition. Often a decomposition is associated with an algorithm, e.g., finding the LU decomposition is equivalent to completing Gaussian Elimination.

🔗

Exercises 2.5.3 Exercises

1.

Solve the systems below by hand using Gaussian elimination and back substitution (exactly as above) on the augmented matrix. As a by-product, give the LU decomposition of A. Pivot wherever appropriate (the number being eliminated should be smaller than the number eliminating it). Check by hand that \(LU=PA\) and \(A\xb = \bb\) and compare with MATLAB.

\(\displaystyle \displaystyle{ A = \left(\begin{array}{cc}0.5 & 4 \\ 2 & 4\end{array} \right), \qquad \bb = \left(\begin{array}{c}0 \\ -3\end{array} \right)}\)
🔗

🔗
\(\displaystyle \displaystyle{ A = \left(\begin{array}{ccc}1 & 2 & 3 \\ 4 & 1 & 1 \\ 2 & 3 & 9\end{array} \right), \qquad \bb = \left(\begin{array}{c}2 \\ 1 \\ 1\end{array} \right)}\)
🔗

🔗

🔗

2.

Finish the following MATLAB function program:

function [x1, r1, x2, r2] = mysolve(A,b)
  % Solves linear systems using the LU decomposition with pivoting
  % and also with the built-in solve function A\b.
  % Inputs: A -- the matrix
  %         b -- the right-hand vector
  % Outputs: x1 -- the solution using the LU method
  %          r1 -- the scalar residual using the LU method
  %          x2 -- the solution using the built-in method
  %          r2 -- the scalar residual using the 
  %                built-in method

Using format long, test the program on both random matrices (randn(n,n)) and Hilbert matrices (hilb(n)) with \(n\) large (as big as you can make it and the program still run). Print your program and summarize your observations. (Do not print any random matrices or vectors.)

🔗

Prev Top Next