Elimination: A = LU

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 a 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 the LU-decomposition or 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 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:

It should come as no surprise that Mathematica has a built-in LU factorisation command, LUDecomposition, with returns both upper and lower triangular matrices as a single n×n matrix, and a pivot list:

{lu, p, c} = LUDecomposition[matrix]
where l is in the strictly lower‐triangular part of lu with ones assumed along the diagonal:
l = lu SparseArray[{i_, j_} /; j < i -> 1, {n, n}] + IdentityMatrix[n];
u is in the upper‐triangular part of lu:
u = lu SparseArray[{i_, j_} /; j >= i -> 1, {n, n}];
Then their product l.u reconstructs the original matrix.

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 column 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.     ▣

To see all steps in LU-factorization without row exchange, you may want to use the following subroutine:
LUfactorization[A0_, n_] := Module[{A = A0}, U = A0; L = IdentityMatrix[n];
Print[MatrixForm[L], MatrixForm[U], " = " , MatrixForm[A0]];
For[p = 1, p <= n - 1, p++,
For[i = p + 1, i <= n, i++,
m = A[[i, p]]/A[[p, p]];
L[[i, p]] = m;
A[[i]] = A[[i]] - m A[[p]]; U = A;
Print[MatrixForm[L], MatrixForm[U], " = ", MatrixForm[A0]]; ];];]
To control every column elimination, one may want to use
PivotDown[m_, {i_, j_}, oneflag_: 0] :=
Block[{k}, If[m[[i, j]] == 0, Return[m]];
Return[Table[
Which[k < i, m[[k]], k > i, m[[k]] - m[[k, j]]/m[[i, j]] m[[i]],
k == i && oneflag == 0, m[[k]] , k == 1 && oneflag == 1,
m[[k]]/m[[i, j]] ], {k, 1, Length[m]}]]]
One may want to use the following subroutine to determine LU-factorization:
LUfactor1[n_] := Module[{}, r = Table[j, {j, n}];
For[p = 1, p <= n - 1, p++,
For[j = p + 1, j <= m, j++,
If[Abs[A[[r[[j]], p]]] > Abs[A[[r[[p]], p]]],
r[[{j, p}]] = r[[{p, j}]] ];];
For[k = p + 1, k <= n, k++,
A[[r[[k]], p]] = A[[r[[k]], p]] / A[[r[[p]], p]] ;
A[[r[[k]], Range[p + 1, n] ]] =
A[[r[[k]], Range[p + 1, n] ]] - A[[r[[k]], p]] A[[r[[p]], Range[p + 1, n] ]]; ];];
L = P = IdentityMatrix[n];
Z = Table[0, {j, n}];
P = P[[r]];
U = A[[r]];
For[i = 1, i <= n, i++,
L[[i, Range[1, i - 1] ]] = A[[r[[i]], Range[1, i - 1] ]];
U[[i, Range[1, i - 1] ]] = Z[[Range[1, i - 1] ]]; ];]

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: Form the decomposition \( {\bf A} = {\bf L} {\bf U} . \)

Recall that the elementary operations on the rows of a matrix 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: For illustration, consider 4×4 elementary matrices
\[ {\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} . \]
E3 = DiagonalMatrix[{1, 1, alpha, 1}]
E42 = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0},{0,beta,0,1}}
E13 = {{0,0,1,0},{0,1,0,0}, {1,0,0,0}, {0,0,0,1}}
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} . \]
Inverse[E3]
Inverse[E42]
Inverse[E13]
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} , \]
A = Array[Subscript[a, #1, #2] &, {4, 4}]
(A) // MatrixForm
we obtain
\[ {\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} , \qquad {\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} . \]

Example: We find LU-factorization of
\[ {\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} . \]

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.

Now we use standard Mathematica commands:

A = {{2, 6, 2}, {-3, -8, 0}, {4, 9, 2}}
lu = LUDecomposition[A][[1]]
l = lu SparseArray[{i_, j_} /; j < i -> 1, {3, 3}] + IdentityMatrix[3];
u = lu SparseArray[{i_, j_} /; j >= i -> 1, {3, 3}];
l.u
{{2, 6, 2}, {-3, -8, 0}, {4, 9, 2}}
l // MatrixForm
\( \begin{pmatrix} 1&0&0 \\ -12&1&0 \\ 18&-\frac{26}{17}&1 \end{pmatrix} \)
u // MatrixForm
\( \begin{pmatrix} 1&4&16 \\ 0&34&185 \\ 0&0&-\frac{18}{17} \end{pmatrix} \)
■

we start by considering the matrix

\[ {\bf E}_1 = {\bf I} - \begin{bmatrix} 0&0&0&\cdots &0 \\ k_{21}&0&0&\cdots &0 \\ k_{31}&0&0&\cdots &0 \\ \vdots&0&0&\cdots &0 \\ k_{n1}&0&0&\cdots &0 \end{bmatrix} , \]
where \( k_{i1} = a_{i1} / a_{11} . \)
Example: Consider the system of algebraic equations
\begin{align*} 3\,x_1 + x_2 + x_3 &= -1 , \\ 2\,x_1 + x_2 + 2\, x_3 &= 4 , \\ x_1 + x_2 + 2\, x_3 &= 0 , \end{align*}
which can be written in matrix-vector form as
\[ \begin{bmatrix} 3&1&1 \\ 2&1&2 \\ 1&1&2 \end{bmatrix} = \begin{bmatrix} -1 \\ 4 \\ 0 \end{bmatrix} \qquad \Longleftrightarrow \qquad {\bf A} \, {\bf x} = {\bf b} . \]
To this system of equations corresponds an augmented matrix
\[ {\bf B} = \left[ \begin{array}{ccc|c} 3&1&1&-1 \\ 2&1&2&4 \\ 1&1&2&0 \end{array} \right] . \]
We ask Mathematica to build the augmented matrix.
A = {{3, 1, 1}, {2, 1, 2}, {1, 1, 2}}
b = {-1, 4, 0}
B = MapThread[Append, {A, b}]
{{3, 1, 1, -1}, {2, 1, 2, 4}, {1, 1, 2, 0}}
or
B = Insert[A // Transpose, b, 4] // Transpose
We apply the first elimination part to the augmented matrix.
B1 = {B[[1]], B[[2]] - (2/3)*B[[1]], B[[3]]}
{{3, 1, 1, -1}, {0, 1/3, 4/3, 14/3}, {1, 1, 2, 0}}
or, using elementary matrix
E1 = {{1, 0, 0}, {-2/3, 1, 0}, {0, 0, 1}}
E1.B
We eliminate the next entry in the first column
B2 = {B1[[1]], B1[[2]], B1[[3]] - (1/3)*B1[[1]]}
{{3, 1, 1, -1}, {0, 1/3, 4/3, 14/3}, {0, 2/3, 5/3, 1/3}}
or multiply by the elementary matrix
E2 = {{1, 0, 0}, {0, 1, 0}, {-1/3, 0, 1}}
E2.E1.B
{{3, 1, 1, 1}, {0, 1/3, 4/3, 4/3}, {0, 2/3, 5/3, 8/3}}
Now we finish the first part---forward elimination:
E3 = {{1, 0, 0}, {0, 1, 0}, {0, -2, 1}}
B3 = E3.E2.E1.B
{{3, 1, 1, 1}, {0, 1/3, 4/3, 4/3}, {0, 0, -1, -9}}
So we reduced our augmented matrix to the upper triangular one:
\[ {\bf B} \,\sim \,{\bf B}3 = \begin{bmatrix} 3&1&1&-1 \\ 0&\frac{1}{3}& \frac{4}{3}&\frac{14}{3} \\ 0&0&-1&-9 \end{bmatrix} . \]
From the last line, it follows that x3 = 9. Using this value, we immediately obtain that x2 = -22 and x1 = 4. We check the answer with Mathematica
LinearSolve[A, b]
{4, -22, 9}

 

Example: In \( \mathbb{R}^n , \) the vectors \( e_1 [1,0,0,\ldots , 0] , \quad e_2 =[0,1,0,\ldots , 0], \quad \ldots , e_n =[0,0,\ldots , 0,1] \) form a basis for n-dimensional real space, and it is called the standard basis. Its dimension is n.

 

Example: Let us consider the set of all real \( m \times n \) matrices, and let \( {\bf M}_{i,j} \) denote the matrix whose only nonzero entry is a 1 in the i-th row and j-th column. Then the set \( {\bf M}_{i,j} \ : \ 1 \le i \le m , \ 1 \le j \le n \) is a basis for the set of all such real matrices. Its dimension is mn.

 

Example: The set of monomials \( \left\{ 1, x, x^2 , \ldots , x^n \right\} \) form a basis in the set of all polynomials of degree up to n. It has dimension n+1. ■

 

Example: The infinite set of monomials \( \left\{ 1, x, x^2 , \ldots , x^n , \ldots \right\} \) form a basis in the set of all polynomials. ■

 

  1. Without row exchange, use elementary matrices to find LU-factorizations for the following matrices. \[ \mbox{(a)} \quad \begin{bmatrix} 1&-2&-2&-3 \\ 3&-9&0&-9 \\ -1&2&4&7 \\ -3&-6&26&2 \end{bmatrix} , \qquad \mbox{(b)} \quad \begin{bmatrix} 2&3&2&3 \\ 4&3&2&3 \\ 8&7&6&9 \\ 9&2&1 &3 \end{bmatrix} ; \qquad \mbox{(c)} \quad \begin{bmatrix} 8&5&7&6 \\ 2&1&3&3 \\ 6&5&7&9 \\ 4&1&0&1 \end{bmatrix} ; \] \[ \mbox{(d)} \quad \begin{bmatrix} 9&3&3&3 \\ 3&10&-2&-2 \\ 3&-2&18&10 \\ 3&-2&10&10 \end{bmatrix} ; \qquad \mbox{(e)} \quad \begin{bmatrix} 3&-9&9&3 \\ 1&2&-2&6 \\ -8&3&3&-8 \\ 9&6&6&3 \end{bmatrix} ; \qquad \mbox{(f)} \quad \begin{bmatrix} 2&4&-4&0 \\ 1&5&-5&-3 \\ 2&3&1&3 \\ 1&4&-2&2 \end{bmatrix} . \]

  1. Higham, Nicholas, Gaussian Elimination, Manchester Institute for Mathematical Sciences, School of Mathematics, The University of Manchester, 2011.