6.5. Elimination

The elimination procedure is used to find solutions to matrix equations. This section reviews the elimination procedure as it is usually performed with pencil and paper. A basic understanding of how to manually implement the elimination procedure is a prerequisite to understanding and applying the functions available in MATLAB for solving systems of equations. Two elimination based MATLAB functions are introduced. The rref function performs Gauss-Jordan elimination and is a useful educational and analysis tool. The lu function implements a variation of elimination known as LU decomposition that is the preferred method for solving square matrix systems of equations.

6.5.1. Triangular Systems

The only system easier to solve than a triangular system is a diagonal system. Our efforts are nearly complete after changing the \bf{A} matrix from a system of equations of the form \mathbf{A}\,\bm{x} = \bm{b} to lower triangular, upper triangular, or a product of the two as found with LU decomposition. Triangular matrices have zeros either above or below the main diagonal. It is not efficient to continue the elimination procedure to get to a diagonal system. Our objective in the elimination procedure is to produce triangular systems of equations.

\begin{array}{ccc} \mat{l_{11} 0 0; l_{21} l_{22} 0; l_{31} l_{32} l_{33}} & \mat{u_{11} u_{12} u_{13}; 0 u_{22} u_{23}; 0 0 u_{33}} & \mat{d_{11} 0 0; 0 d_{22} 0; 0 0 d_{33}} \\ \hfill \\[-13pt] \text{Lower} & \text{Upper} & \text{Diagonal} \\[-2pt] \text{Triangular} & \text{Triangular} & {} \end{array}

The solution to a lower triangular matrix equation uses a forward substitution procedure where we start solving for unknowns on the top row and work down from there. After elimination, the equation \mathbf{A}\,\bm{x} = \bm{b} takes the form \mathbf{L}\,\bm{x} = \bm{c}.

(6.14)\mat{l_{11} 0 0; l_{21} l_{22} 0; l_{31} l_{32} l_{33}}\, \vector{x_1; x_2; x_3} = \vector{c_1; c_2; c_3} \qquad \left | \enspace \begin{aligned} x_1 &= \frac {c_1}{l_{11}} \\ x_2 &= \frac {c_2 - l_{21}\,x_1}{l_{22}} \\ x_3 &= \frac {c_3 - l_{31}\,x_1 - l_{32}\,x_2}{l_{33}} \end{aligned} \right.

The solution to an upper triangular matrix system uses a back substitution procedure where we start solving for unknowns on the bottom row and work up from there. After elimination, the equation takes the form \mathbf{U}\,\bm{x} = \bm{c}.

(6.15)\mat{u_{11} u_{12} u_{13}; 0 u_{22} u_{23}; 0 0 u_{33}} \, \vector{x_1; x_2; x_3} = \vector{c_1; c_2; c_3} \qquad \left | \enspace \begin{aligned} x_3 &= \frac {c_3}{u_{33}} \\ x_2 &= \frac {c_2 - u_{23}\,x_3}{u_{22}} \\ x_1 &= \frac {c_1 - u_{12}\,x_2 - u_{13}\,x_3}{l_{11}} \end{aligned} \right.

6.5.2. The Gaussian Elimination Procedure

Gaussian elimination and Gauss–Jordan elimination [1] are procedures for changing a matrix to upper triangular form where the solution is simple to calculate. A sequence of linear operations is applied to the matrix’s rows to shift matrix elements below the diagonal to zero.

Three operations are allowed in elimination.

  1. Reorder the rows, which is called row exchanges or partial pivoting. This step can be critical to finding an accurate solution to some systems of equations. In Numerical Stability, we will address row exchanges in the context of numerical accuracy. Initially, we will use systems of equations where row exchanges are not needed.
  2. Add a multiple of one row to another row, replacing that row.
  3. Multiply a row by a nonzero constant.

When a row has all zeros to the left of the main diagonal, that nonzero element on the diagonal is called the pivot. The pivot is used to determine a multiple of the row to add to another row to produce a needed zero.

We begin with our classic matrix equation \mathbf{A}\,\bm{x} = \bm{b}. Any operations done on the \bf{A} matrix must also be applied to the \bm{b} vector. Note that one could continue the elimination steps until \bf{U} is a diagonal matrix, but it is faster to stop at an upper triangular matrix and use back substitution ( equation (6.15)) as illustrated in the example below.

To help carry forward the operations to the right side of the equation, we form an augmented matrix. Using the numbers from our previous example, we have the following. The pivots in each step are underlined.

\spalignaugmat{\underline{2} -3 0 3; 4 -5 1 9; 2 -1 -3 -1}

Add -1 of row 1 to row 3.

\spalignaugmat{\underline{2} -3 0 3; 4 -5 1 9; 0 2 -3 -4}

Add -2 of row 1 to row 2. The pivot then moves to row 2.

\spalignaugmat{2 -3 0 3; 0 \underline{1} 1 3; 0 2 -3 -4}

Add -2 of row 2 to row 3 to finish the row operations.

\spalignaugmat{2 -3 0 3; 0 1 1 3; 0 0 -5 -10}

Our matrix equation is now upper triangular.

\spalignmat{2 -3 0; 0 1 1; 0 0 -5} \spalignvector{x_1;x_2;x_3} = \spalignvector{3;3;-10}

The final step is to apply back substitution as shown in equation (6.15).

\spalignsys{ \. \- \. - 5\,x_3 \+ \. = -10; \. \+ x_2 \+ \. + 2 = 3; 2\,x_1 \+ \. \+ \. - 3 = 3;}

x_1 = 3, x_2 = 1, x_3 = 2

Practice Problem

Here is another system of equations with an integer solution. Use elimination to solve for variables x_1, x_2, and x_3. Then use MATLAB’s left-divide operator as demonstrated in Jumping Ahead to MATLAB to verify your answer.

\spalignsys{-3x_1 + 2x_2 - x_3 = -1; 6x_1 - 6x_2 + 7x_3 = -7; 3x_1 - 4x_2 + 4x_3 = -6}

6.5.3. Elimination to Find the Matrix Inverse

Gaussian elimination can calculate the inverse of a matrix when we start with an augmented matrix of the form \left[\bf{A}\,|\,\bf{I}\right] and do row operations until we have \left[\mathbf{I}\,|\,\mathbf{A}^{-1}\right].

\spalignaugmathalf{\underline{2} -3 0 1 0 0; 4 -5 1 0 1 0; 2 -1 -3 0 0 1;}

Add -1 of row 1 to row 3.

\spalignaugmathalf{\underline{2} -3 0 1 0 0; 4 -5 1 0 1 0; 0 2 -3 -1 0 1}

Add -2 of row 1 to row 2. The pivot then moves to row 2.

\spalignaugmathalf{2 -3 0 1 0 0; 0 \underline{1} 1 -2 1 0; 0 2 -3 -1 0 1}

Add -2 of row 2 to row 3.

\spalignaugmathalf{2 -3 0 1 0 0; 0 \underline{1} 1 -2 1 0; 0 0 -5 3 -2 1}

Add 3 of row 2 to row 1. The pivot then moves to row 3.

\spalignaugmathalf{2 0 3 -5 3 0; 0 1 1 -2 1 0; 0 0 \underline{-5} 3 -2 1}

Add 3/5 of row 3 to row 1; and add 1/5 of row 3 to row 2. Then a pivot is no longer needed.

\spalignaugmathalf{2 0 0 -16/5 9/5 3/5; 0 1 0 -7/5 3/5 1/5; 0 0 -5 3 -2 1}

Divide row 1 by 2; and divide row 3 by -5.

\spalignaugmathalf{1 0 0 -16/10 9/10 3/10; 0 1 0 -7/5 3/5 1/5; 0 0 1 -3/5 2/5 -1/5}

\mathbf{A}^{-1} = \spalignmat{-1.6 0.9 0.3; -1.4 0.6 0.2; -0.6 0.4 -0.2}

Note

This result matches what MATLAB found in Calculating a Matrix Inverse. All of the tedious row operations certainly makes one appreciate that MATLAB can perform the calculations for us.

6.5.4. Reduced Row Echelon Form

The word echelon just means resembling stair steps. A matrix is in row echelon form when:

  • All nonzero values in each column are above all zeros.
  • The leading coefficient of each row is strictly to the right of the leading coefficient in the row above it.

When Gaussian elimination is used on an augmented matrix to make an upper triangular matrix, the matrix is put into row echelon form.

Reduced row echelon form (RREF) requires that the use the Gauss-Jordan elimination procedure is used. The elimination procedure continues to address elements above and on the diagonal. The additional requirements for reduced row echelon form are:

  • Every leading coefficient must be 1.
  • The leading coefficients must be the only nonzero value in its column. In RREF, the first m (number of rows) columns of the matrix should look as close as possible to an identity matrix.

A full rank matrix system in augmented form has the identity matrix in the first m columns when in RREF. The last column is the solution to the system of equations.

MATLAB’s rref function puts a matrix or augmented matrix into RREF.

>> Ab = [A b]
Ab =
    8     1     6    27
    3     5     7    38
    4     9     2    55
>> rref(Ab)
ans =
    1     0     0     2
    0     1     0     5
    0     0     1     1

As we saw from the CR factoring in The Column and Row Factorization, the RREF of a singular matrix shows how dependent columns are weighted sums of the independent columns. With r = \text{rank}(\mathbf{A}), the last (n - r) rows are all zeros.

>> A =
      1    -8    -7
      4    -5    -1
      6     2     8
>> rref(A)
ans =
      1     0     1
      0     1     1
      0     0     0

Footnote

[1]Before row operations are performed using a pivot variable, the Gauss–Jordan elimination algorithm divides the pivot row by the pivot value. Thus the pivot values used are always one. The elimination process also clears elements above the diagonal. This step is needed when the objective is RREF ( Reduced Row Echelon Form). For the sake of simplifying the explanation, we will use the Gaussian elimination procedure in our examples, which does not take this step.