Preface

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0340. It is primarily for students who have some experience using Mathematica. If you have never used Mathematica before and would like to learn more of the basics for this computer algebra system, it is strongly recommended looking at the APMA 0330 tutorial. As a friendly reminder, don'r forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in regular fonts. This means that you can copy and paste all comamnds into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts to your needs for learning how to use the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0340
Return to the main page for the course APMA0330
Return to Part II of the course APMA0340

LU-factorizationΒΆ

Recall that in Gaussian Elimination, row operations are used to change the coefficient matrix to an upper triangular matrix. The solution is then found by back substitution, starting from the last equation in the reduced system. In Gauss-Jordan Reduction, row operations are used to diagonalize the coefficient matrix, and the answers are read directly.

The goal of this section is to identify Gaussian elimination with LU factorization. The original m-by-n matrix A becomes the product of two or more special matrices that are actually triangular matrices. Namely, the factorization comes from elimination: \( {\bf A} = {\bf L}\,{\bf U} , \) where L is lower triangular m-by-m matrix and U is upper m-by-n triangular matrix. Such representation is called and LU-decomposition of LU-factorization. Computers usually solve square systems of linear equations using the LU decomposition, and it is also a key step when inverting a matrix, or computing the determinant of a matrix. The LU decomposition was introduced by a Polish astronomer, mathematician, and geodesist Tadeusz Banachiewicz (1882--1954) in 1938.

For a single linear system \( {\bf A}\,{\bf x} = {\bf b} \) of n equations in n unknowns, the methods LU-decomposition and Gauss-Jordan elimination differ in bookkeeping but otherwise involve the same number of flops. However, LU-factorization has the following advantages:

Not every square matrix has an LU-factorization. However, if it is possible to reduce a square matrix A to row echelon form by Gaussian elimination without performing any row interchanges, then A will have an LU-decomposition, though it may not be unique.

LU-decomposition:
Step 1: rewrite the system of algebraic equations \( {\bf A} {\bf x} = {\bf b} \) as

\[ {\bf L}{\bf U} {\bf x} = {\bf b} . \]
Step 2: Define a new \( n\times 1 \) matrix y (which is actually a vector) by
\[ {\bf U} {\bf x} = {\bf y} . \]
Step 3: Rewrite the given equation as \( {\bf L} {\bf y} = {\bf b} \) and solve this sytem for y.
Step 4: Substitute y into the equation \( {\bf U} {\bf x} = {\bf y} \) and solve for x.

Procedure for constructing LU-decomposition:
Step 1: Reduce \( n \times n \) matrix A to a row echelon form U by Gaussian elimination without row interchanges, keeping track of the multipliers used to introduce the leading coefficients (usually 1's) and multipliers used to introduce the zeroes below the leading coefficients.
Step 2: In each position along the main diagonal L, place the reciprocal of the multiplier that introduced the leading 1 in that position of U.
Step 3: In each position below the main diagonal of L, place the negative of the multiplier used to introduce the zero in that position of U.
Step 4: For the decomposition \( {\bf A} = {\bf L} {\bf U} . \)

Recall the elementary operations on the rows of a matrix that are equivalent to premultiplying by an elementary matrix E:

(1) multiplying row i by a nonzero scalar α , denoted by \( {\bf E}_i (\alpha ) ;\)
(2) adding β times row j to row i, denoted by \( {\bf E}_{ij} (\beta ) \) (here β is any scalar), and
(3) interchanging rows i and j, denoted by \( {\bf E}_{ij} \) (here \( i \ne j \) ), called elementary row operations of types 1,2 and 3, respectively.

Correspondingly, a square matrix E is called an elementary matrix if it can be obtained from an identity matrix by performing a single elementary operation.

Theorem: Every elementary matrix is invertible, and the inverse is also an elementary matrix. ■

Example. Illustrations for m = 4:

\[ {\bf E}_3 (\alpha ) = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&\alpha &0 \\ 0&0&0&1 \end{bmatrix} , \qquad {\bf E}_{42} (\beta ) = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1 &0 \\ 0&\beta&0&1 \end{bmatrix} , \qquad {\bf E}_{13} = \begin{bmatrix} 0 &0&1&0 \\ 0&1&0&0 \\ 1&0&0 &0 \\ 0&0&0&1 \end{bmatrix} . \]
Correspondingly, we have their inverses:
\[ {\bf E}_3^{-1} (\alpha ) = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&\alpha^{-1} &0 \\ 0&0&0&1 \end{bmatrix} , \qquad {\bf E}_{42}^{-1} (\beta ) = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1 &0 \\ 0&-\beta&0&1 \end{bmatrix} , \qquad {\bf E}_{13}^{-1} = \begin{bmatrix} 0 &0&1&0 \\ 0&1&0&0 \\ 1&0&0 &0 \\ 0&0&0&1 \end{bmatrix} . \]
Multiplying by an arbitrary matrix
\[ {\bf A} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} , \]
we obtain
\begin{align*} {\bf E}_3 (\alpha )\,{\bf A} &= \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ \alpha \,a_{31} & \alpha \,a_{32} & \alpha \,a_{33} & \alpha \,a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} , \quad {\bf A} \, {\bf E}_3 (\alpha ) = \begin{bmatrix} a_{11} & a_{12} & \alpha \,a_{13} & a_{14} \\ a_{21} & a_{22} & \alpha \,a_{23} & a_{24} \\ a_{31} & a_{32} & \alpha \,a_{33} & a_{34} \\ a_{41} & a_{42} & \alpha \,a_{43} & a_{44} \end{bmatrix} , \\ {\bf E}_{42} (\beta ) \, {\bf A} &= \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ \alpha \,a_{31} & \alpha \,a_{32} & \alpha \,a_{33} & \alpha \,a_{34} \\ a_{41} + \beta \,a_{21} & a_{42} + \beta \,a_{22} & a_{43} + \beta \,a_{23}& a_{44} + \beta \,a_{24} \end{bmatrix} , \quad \\ {\bf A} \, {\bf E}_{42} (\beta ) &= \begin{bmatrix} a_{11} & a_{12} + \beta\, a_{14} & a_{13} & a_{14} \\ a_{21} & a_{22} + \beta \,a_{24}& a_{23} & a_{24} \\ a_{31} & a_{32} + \beta\, a_{34}& a_{33} & a_{34} \\ a_{41} & a_{42} + \beta \,a_{44}& a_{43} & a_{44} \end{bmatrix} , \\ {\bf E}_{13} \,{\bf A} &= \begin{bmatrix} a_{31} & a_{32} & a_{33} & a_{34} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{11} & a_{12} & a_{13} & a_{14} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} = {\bf A}\, {\bf E}_{13} . \end{align*}

Example. We find LU-factorization of the following 3-by-3 matrix:

\[ {\bf A} = \begin{bmatrix} 2&6&2 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} . \]

To achieve this, we reduce A to a row echelon form U using Gaussian elimination and then calculate L by inverting the product of elementary matrices.

Reduction Row operation Elementary matrix Inverse matrix
\( \begin{bmatrix} 2&6&2 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} \)      
Step 1: \( \frac{1}{2} \times \,\mbox{row 1} \)   \( {\bf E}_1 = \begin{bmatrix} 1/2&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} \) \( {\bf E}_1^{-1} = \begin{bmatrix} 2&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} \)      
Step 2 \( 3 \times \,\mbox{row 1} + \mbox{row 2}\) \( {\bf E}_2 = \begin{bmatrix} 1&0&0 \\ 3&1&0 \\ 0&0&1 \end{bmatrix} \)    \( {\bf E}_2^{-1} = \begin{bmatrix} 1&0&0 \\ -3&1&0 \\ 0&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 4&9&2 \end{bmatrix} \)      
Step 3 \( -4 \times \,\mbox{row 1} + \mbox{row 3}\) \( {\bf E}_3 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ -4&0&1 \end{bmatrix} \) \( {\bf E}_3^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 4&0&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&-3&-2 \end{bmatrix} \)      
Step 4 \( 3 \times \,\mbox{row 2} + \mbox{row 3}\) \( {\bf E}_4 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&3&1 \end{bmatrix} \) \( {\bf E}_4^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&-3&1 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&0&7 \end{bmatrix} \)      
Step 5 \( \frac{1}{7} \times \,\mbox{row 3} \) \( {\bf E}_5 = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&1/7 \end{bmatrix} \) \( {\bf E}_5^{-1} = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&7 \end{bmatrix} \)
\( \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&0&1 \end{bmatrix} = {\bf U}\)      
Now we find the lower triangular matrix:
\begin{align*} {\bf L} &= {\bf E}_1^{-1} {\bf E}_2^{-1} {\bf E}_3^{-1} {\bf E}_4^{-1} {\bf E}_5^{-1} \\ &= \begin{bmatrix} 2&0&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix} \begin{bmatrix} 1&0&0 \\ -3&1&0 \\ 0&0&1 \end{bmatrix} \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 4&0&1 \end{bmatrix} \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&-3&1 \end{bmatrix} \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&7 \end{bmatrix} \\ &= \begin{bmatrix} 2&0&0 \\ -3&1&0 \\ 4&-3&7 \end{bmatrix} , \end{align*}
so
\[ {\bf A} = \begin{bmatrix} 2&6&2 \\ -3&-8&0 \\ 4&9&2 \end{bmatrix} = \begin{bmatrix} 2&0&0 \\ -3&1&0 \\ 4&-3&7 \end{bmatrix}\, \begin{bmatrix} 1&3&1 \\ 0&1&3 \\ 0&0&1 \end{bmatrix} = {\bf L}\,{\bf U} . \]

We can check LU-factorization with the following Mathematica code. However, this factorization is not unique. First, we introduce subroutine:
LUdecompose[n_, A_] := Module[{AA, multiplier},
L = Array[0, {n, n}];
U = Array[0, {n, n}];
AA = Array[0, {n, n}];
For[i = 1, i <= n, i++, L[[i, i]] = 1;
For[j = 1, j <= n, j++, AA[[i, j]] = A[[i, j]];
If[j != i, L[[i, j]] = 0]]];
For[k = 1, k <= n - 1, k++,
For[i = k + 1, i <= n, i++,
multiplier = N[AA[[i, k]]/AA[[k, k]]];
L[[i, k]] = multiplier;
For[j = 1, j <= n, j++,
AA[[i, j]] = AA[[i, j]] - multiplier*AA[[k, j]]]]
For[i = 1, i <= n, i++,
For[j = 1, j <= n, j++, U[[i, j]] = AA[[i, j]];
If[j < i, U[[i, j]] = 0]]];];
Print["L", "=", L // MatrixForm];
Print["U", "=", U // MatrixForm]]
Then we apply this subroutine to our matrix:
A = {{2, 6, 2}, {-3, -8, 0}, {4, 9, 2}}
LUdecompose[3, A]
L= \( \left( \begin{array}{cccc} 1& 0& 0 \\ -1.5& 1& 0 \\ 2.& -3.& 1 \end{array} \right) \)
U= \( \left( \begin{array}{cccc} 2& 6& 2 \\ 0& 1.& 3. \\ 0& 0& 7. \end{array} \right) \)
Finally, we check the answer:
L.U
{{2., 6., 2.}, {-3., -8., 0.}, {4., 9., 2.}}
Notice that matrix \( {\bf L}_{n \times n} = [l_{i,j}] \) has only one unknown to be solved for in its first row. Once the matrix L has been formed, forward substitution step, L y = b, can be conducted, beginning with the first equation as it has only one unknown:
\[ y_1 = \frac{b_1}{l_{1,1}} . \]
Subsequent steps of forward substitution can be represented by the following formula:
\[ y_i (x) = \frac{1}{l_{i,i}} \left( b_i - \left[ \sum_{j=1}^{i-1} l_{i,j} * b_j \right]_{i=2..n} \right) . \]
The procedure below conducts forward substitution steps to solve for \( {\bf y}_{n \times 1} . \)
forwardsubstitution[n_,L_,b_]:=Module[{},
y=Array[y,n];
y[[1]]=b[[1]]/L[[1,1]];
For[i=2,i<=n,i++,summ=0;
For[j=1,j<=i-1,j++,summ=summ+L[[i,j]]*y[[j]]];
y[[i]]=(b[[i]]-summ)/L[[i,i]]];
Print["y","=",y//MatrixForm]]
Now that y has been found, it can be used in the back substitution step, U x = y, to solve for solution vector \( {\bf x}_{n \times 1} , \) where \( {\bf U}_{n \times n} = [u_{i,j}] \) is the upper triangular matrix calculated in LU-decomposition algorithm, and \( {\bf y}_{n \times 1} \) is the right hand side array.

Back substitution begins with solving the n-th equation as it has only one unknown:

\[ x_n = \frac{y_n}{u_{n,n}} . \]
The remaining unknowns are solved for working backwards from the (n-1)-th equation to the first equation using the formula below:
\[ x_i (x) = \frac{1}{u_{i,i}} \left( y_i - \left[ \sum_{j=i+1}^{n} u_{i,j} * x_j \right]_{i=n-1..1} \right) . \]
The following procedure solves for \( {\bf x}_{n \times 1} . \)
backsubstitution[n_,U_,y_]:=Module[{},
x=Array[x,n];
x[[n]]=y[[n]]/U[[n,n]];
For[i=n-1,i>=1,i--, summ=0;
For[j=i+1,j<=n,j++,summ=summ+U[[i,j]]*X[[j]]];
X[[i]]=(y[[i]]-summ)/U[[i,i]]];
Print["X","=",X//MatrixForm]]

Example. Consider the \( 3\times 4 \) matrix

\[ {\bf A} = \begin{bmatrix} 1&0&2&3 \\ 2&-1&3&6 \\ 1&4&4&0 \end{bmatrix} \]
and consider the elementary matrix
\[ {\bf E}_{31} (-1) = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ -1&0&1 \end{bmatrix} , \]
which results from adding \( -1 \) times the first row of \( {\bf I}_3 \) to the thrid row. The product \( {\bf E}_{31} (-1)\, {\bf A} \) is
\[ {\bf E}_{31} (-1)\,\, {\bf A} = \begin{bmatrix} 1&0&2&3 \\ 2&-1&3&6 \\ 0&4&2&-3 \end{bmatrix} , \]
which is precisely the matrix that results when we add \( -1 \) times the first row of A to the third row. Next, we multiply by matrix \( {\bf E}_{21} (-2) = \begin{bmatrix} 1&0&0 \\ -2&1&0 \\ 0&0&1 \end{bmatrix} . \) This yields
\[ {\bf E}_{21} (-2)\,{\bf E}_{31} (-1)\,\, {\bf A} = \begin{bmatrix} 1&0&2&3 \\ 0&-1&-1&0 \\ 0&4&2&-3 \end{bmatrix} , \]
Finally, we multiply by \( {\bf E}_{32} (4) = \begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&4&1 \end{bmatrix} \) to obtain
\[ {\bf E}_{32} (4)\,{\bf E}_{21} (-2)\,{\bf E}_{31} (-1)\,\, {\bf A} = \begin{bmatrix} 1&0&2&3 \\ 0&-1&-1&0 \\ 0&0&-2&-3 \end{bmatrix} = {\bf U} . \]
Hence \( {\bf A} = {\bf E}_{31}^{-1} \,{\bf E}_{21}^{-1} \, {\bf E}_{32}^{-1} \, {\bf U} . \) Evaluation the product of three matrices
\[ {\bf E}_{31}^{-1} \,{\bf E}_{21}^{-1} \, {\bf E}_{32}^{-1} = \begin{bmatrix} 1&0&0 \\ 2&1&0 \\ 1&-4&1 \end{bmatrix} = {\bf L} , \]
we obtain the lower triangular matrix L, so \( {\bf A} = {\bf L} \,{\bf U} \) is the required LU-factorization. ■

 

Mathematica solves linear systems

Elementary Row Operations

Row Echelon Forms

LU-factorization

PLU Factorization

Reflection

Givens rotation

Row Space and Column Space

Null Space or Kernel

Rank and Nullity

Square Matrices

Cramer's Rule

Symmetric and self-adjoint matrices

Cholesky decomposition

Projection Operators

Gram-Schmidt

QR-decomposition

Least Square Approximation

SVD Factorization

Numerical Methods

Markov chains

Inequalities

Miscellaneous

 

Return to Mathematica page

Return to the main page (APMA0330)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Equations
Return to the Part 3 Linear Systems of Ordinary Differential Equations
Return to the Part 4 Non-linear Systems of Ordinary Differential Equations
Return to the Part 5 Numerical Methods
Return to the Part 6 Fourier Series
Return to the Part 7 Partial Differential Equations