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:

  • Gaussian elimination and Gauss--Jordan elimination both use the augmented matrix \( \left[ {\bf A} \, | \, {\bf b} \right] ,\) so b must be known. In contrast, LU-decomposition uses only matrix A, so once that factorization is complete, it can be applied to any vector b.
  • For large linear systems in which computer memory is at a premium, one can dispense with the storage of the 1's and zeroes that appear on or below the main diagonal of U, since those entries are known. The space that this opens up can then be used to store the entries of L.

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 system 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 multiplying 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 1: 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} . \]
End of Example 1

Example 2: 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} . \]

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]]
End of Example 2

Example 3: Consider the 3×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 I3 to the thrid row. The product \( {\bf E}_{31} (-1)\, {\bf A} \) is
\[ {\bf E}_{31} (-1)\,\, {\bf A} = \begin{bmatrix} 1&\phantom{-}0 &2&\phantom{-}3 \\ 2&-1&3&\phantom{-}6 \\ 0&\phantom{-}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&\phantom{-}0&\phantom{-}2&\phantom{-}3 \\ 0&-1&-1&\phantom{-}0 \\ 0&\phantom{-}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&\phantom{-}0&0 \\ 2&\phantom{-}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} \)
End of Example 3
■

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 4: 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}
End of Example 4

 

  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.