.. _elimination: Elimination =========== .. index:: rref .. index:: lu The elimination procedure is used to find solutions to square matrix equations. This section reviews the manual elimination procedure 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 Python for solving systems of equations. Two elimination-based functions are introduced. The ``rref`` function from the SymPy module performs Gauss-Jordan elimination and is a valuable educational and analysis tool. The ``lu`` function from the SciPy linalg module implements a variation of elimination known as LU decomposition, which is the preferred method for solving square matrix systems of equations. .. index:: pivot .. index:: elimination .. index:: Gaussian elimination .. index:: augmented matrix Triangular Systems ------------------ A diagonal system is the only system that is easier to solve than a triangular system. Our efforts are nearly complete after changing the :math:`\bf{A}` matrix from a system of equations of the form :math:`\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. Continuing the elimination procedure to get to a diagonal system is unnecessary. .. math:: \begin{array}{ccc} \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} & \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} & \begin{bmatrix} d_{11} & 0 & 0 \\ 0 & d_{22} & 0 \\ 0 & 0 & d_{33} \end{bmatrix} \\ %\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 :math:`\mathbf{A}\,\bm{x} = \bm{b}` takes the form :math:`\mathbf{L}\,\bm{x} = \bm{c}`. .. math:: :label: eq-forward-sub \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix}\, \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix} \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 :math:`\mathbf{U}\,\bm{x} = \bm{c}`. .. math:: :label: eq-back-sub \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} \, \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix} \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. .. _gaussElim: Gaussian Elimination -------------------- Gaussian elimination and Gauss–Jordan elimination  [1]_ are procedures for changing a matrix to upper triangular form. 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. #. Reorder the rows, which is called row exchanges or partial pivoting. Row exchanges can be critical to finding an accurate solution to some systems of equations. We will address row exchanges in the context of numerical accuracy in :ref:`num-stability`. Initially, however, we will use systems of equations where row exchanges are not needed. #. Add a multiple of one row to another, replacing that row. #. Multiply a row by a nonzero constant. When a row has all zeros to the left of the main diagonal, the 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 :math:`\mathbf{A}\,\bm{x} = \bm{b}`. Any operations done on the :math:`\bf{A}` matrix must also be applied to the :math:`\bm{b}` vector. So, we form an augmented matrix. Using the numbers from our previous example, we have the following. The pivots in each step are underlined. .. math:: \left[\begin{array}{ccc|c} \underline{2} & -3 & 0 & 3 \\ 4 & -5 & 1 & 9 \\ 2 & -1 & -3 & -1 \end{array}\right] Add :math:`-1` of *row 1* to *row 3*. .. math:: \left[\begin{array}{ccc|c} \underline{2} & -3 & 0 & 3 \\ 4 & -5 & 1 & 9 \\ 0 & 2 & -3 & -4 \end{array}\right] Add :math:`-2` of *row 1* to *row 2*. The pivot then moves to *row 2*. .. math:: \left[\begin{array}{ccc|c} 2 & -3 & 0 & 3 \\ 0 & \underline{1} & 1 & 3 \\ 0 & 2 & -3 & -4 \end{array}\right] Add :math:`-2` of *row 2* to *row 3* to finish the row operations. .. math:: \left[\begin{array}{ccc|c} 2 & -3 & 0 & 3 \\ 0 & 1 & 1 & 3 \\ 0 & 0 & -5 & -10 \end{array}\right] Our matrix equation is now upper triangular. .. math:: \begin{bmatrix} 2 & -3 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & -5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 3 \\ 3 \\ -10 \end{bmatrix} The final step is to apply back substitution as shown in equation :eq:`eq-back-sub`. .. math:: \begin{aligned} &- 5\,x_3 = -10 \\ &x_2 + 2 = 3 \\ &2\,x_1 - 3 = 3 \end{aligned} .. math:: x_1 = 3, x_2 = 1, x_3 = 2 .. describe:: Practice Problem Here is another system of equations with an integer solution. Use elimination to solve for variables :math:`x_1`, :math:`x_2`, and :math:`x_3`. Then use the ``solve`` function as demonstrated in :ref:`python-solve` to verify your answer. .. math:: \left\{\begin{aligned} -3\,x_1 + 2\,x_2 - x_3 &= -1 \\ 6\,x_1 - 6\,x_2 + 7\,x_3 &= -7 \\ 3\,x_1 - 4\,x_2 + 4\,x_3 &= -6 \end{aligned}\right. .. _elim-inverse: 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 :math:`\left[\bf{A}\,|\,\bf{I}\right]` and do row operations until we have :math:`\left[\mathbf{I}\,|\,\mathbf{A}^{-1}\right]`. .. math:: \left[\begin{array}{ccc|ccc} \underline{2} & -3 & 0 & 1 & 0 & 0 \\ 4 & -5 & 1 & 0 & 1 & 0 \\ 2 & -1 & -3 & 0 & 0 & 1 \end{array}\right] Add :math:`-1` of *row 1* to *row 3*. .. math:: \left[\begin{array}{ccc|ccc} \underline{2} & -3 & 0 & 1 & 0 & 0 \\ 4 & -5 & 1 & 0 & 1 & 0 \\ 0 & 2 & -3 & -1 & 0 & 1 \end{array}\right] Add :math:`-2` of *row 1* to *row 2*. The pivot then moves to *row 2*. .. math:: \left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 0 & \underline{1} & 1 & -2 & 1 & 0 \\ 0 & 2 & -3 & -1 & 0 & 1 \end{array}\right] Add :math:`-2` of *row 2* to *row 3*. .. math:: \left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 0 & \underline{1} & 1 & -2 & 1 & 0 \\ 0 & 0 & -5 & 3 & -2 & 1 \end{array}\right] Add :math:`3` of *row 2* to *row 1*. The pivot then moves to *row 3*. .. math:: \left[\begin{array}{ccc|ccc} 2 & 0 & 3 & -5 & 3 & 0 \\ 0 & 1 & 1 & -2 & 1 & 0 \\ 0 & 0 & \underline{-5} & 3 & -2 & 1 \end{array}\right] Add :math:`3/5` of *row 3* to *row 1*; and add :math:`1/5` of *row 3* to *row 2*. Then a pivot is no longer needed. .. math:: \left[\begin{array}{ccc|ccc} 2 & 0 & 0 & -16/5 & 9/5 & 3/5 \\ 0 & 1 & 0 & -7/5 & 3/5 & 1/5 \\ 0 & 0 & -5 & 3 & -2 & 1 \end{array}\right] Divide *row 1* by :math:`2`; and divide *row 3* by :math:`-5`. .. math:: \left[\begin{array}{ccc|ccc} 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 \end{array}\right] .. math:: \mathbf{A}^{-1} = \begin{bmatrix} -1.6 & 0.9 & 0.3 \\ -1.4 & 0.6 & 0.2 \\ -0.6 & 0.4 & -0.2 \end{bmatrix} .. index:: matrix inverse .. index:: Gaussian elimination .. index:: inverse of matrix .. _rref: Reduced Row Echelon Form ------------------------ .. index:: rref The word *echelon* 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 using the Gauss-Jordan elimination procedure. 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 their column. In RREF, the matrix’s first :math:`m` (number of rows) columns should look as close as possible to an identity matrix. In RREF, an augmented full rank matrix system has the identity matrix in the first :math:`m` columns. The last column is the solution to the system of equations. The ``rref`` function from the SymPy module puts a matrix or augmented matrix into RREF. The first value returned is the RREF matrix. The second value returned is the pivot column of each row, which we ignore in the following example. :: In [58]: import numpy as np In [59]: from sympy import Matrix In [60]: A = np.array([[8, 1, 6], [3, 5, 7], [4, 9, 2]]) In [61]: b = np.array([[27, 38, 55]]).T In [62]: Ab = np.hstack((A, b)); Ab Out[62]: array([[ 8, 1, 6, 27], [ 3, 5, 7, 38], [ 4, 9, 2, 55]]) In [63]: mAb = Matrix(Ab) In [64]: RR, _ = mAb.rref(); In [65]: RR = np.array(RR) In [66]: print(RR) [[1 0 0 2] [0 1 0 5] [0 0 1 1]] As we saw from the CR factoring in :ref:`CR`, the RREF of a singular matrix shows how dependent columns are weighted sums of the independent columns. With :math:`r = \text{rank}(\mathbf{A})`, the last :math:`(n - r)` rows are all zeros. :: In [67]: A = Matrix([[1, -8, -7], [4, -5, -1], [6, 2, 8]]) In [68]: Rref, _ = A.rref() In [69]: print(np.array(Rref)) [[1 0 1] [0 1 1] [0 0 0]] .. index:: RREF .. index:: reduced row echelon form .. rubric:: Footnote .. [1] The Gauss–Jordan elimination procedure divides the pivot row by the pivot value before row operations are performed. Thus, the pivot values used are always one. The elimination process also clears elements above the diagonal. We will use the Gaussian elimination procedure in our examples, which does not take this step.