# Solving Linear Systems via Gaussian Elimination


```
%matplotlib inline
```

We remember that the problem we are trying to solve is the system of linear equations written as

\begin{equation}
  A {\bf x} = {\bf b}.
\end{equation}

This is equivalent to $n$ simultaneous linear equations for the $n$ unknowns $x_i$. 

The first class of methods we consider are **direct**: a finite number of steps are followed which, in principle, result in the correct answer.

Most direct methods use the same core approach: convert the linear system to an equivalent problem that is easy to solve. 

So the first thing we should do is consider which linear systems are easy to solve. Obviously if the matrix is diagonal it is easy to solve. More generally, it is straightforward to solve a linear system if the matrix is *triangular*.

As an example, consider

\begin{equation}
  \begin{pmatrix} 1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 6 \\ 9 \\ 6 \end{pmatrix}.
\end{equation}

Perform **back substitution** by solving from the bottom up. That is, read off the equation defined by the bottom row of the matrix:

\begin{equation}
  6 x_3 = 6 \qquad \implies \qquad x_3 = 1.
\end{equation}

Then read off the equation defined by the second row of the matrix, using the information found for $x_3$:

\begin{equation}
  4 x_2 + 5 x_3 = 9 \qquad \implies \qquad 4 x_2 = 4 \qquad \implies \qquad x_2 = 1.
\end{equation}

Finally, read off the equation defined by the first row of the matrix, using the information found for $x_2$ and $x_3$:

\begin{equation}
  x_1 + 2 x_2 + 3 x_3 = 6 \qquad \implies \qquad x_1 = 1.
\end{equation}

Back substitution will work for any *upper* triangular matrix, provided that no diagonal element is zero (guaranteed if the matrix is non-singular). If the matrix is lower triangular we would start at the top and work down, using **forward** substitution.

## Gaussian elimination

So, we need to be able to convert our original linear system to an equivalent system in triangular form - here we will aim to construct an upper triangular matrix. By equivalent, we mean that the solution is the same.

Think of each row corresponding to a single linear equation, and each column corresponding to a single unknown variable. Then we see that there are three basic operations that can be performed without changing the solution:

1. Multiply any row by a non-zero scalar (equivalent to scaling the whole equation).
2. Swap two rows (equivalent to changing the order of the equations).
3. Swap two columns (equivalent to changing the order of the unknown variables).

We will concentrate on the first two operations, as swapping the columns requires us to keep track of which entries of the final solution vector correspond to which unknown variables.

To ensure that the operations are performed consistently to both the entries of the matrix $A$ *and* the right hand side vector ${\bf b}$, it is best to construct the **augmented matrix** $\left( A | {\bf b} \right)$. All operations are applied to this matrix, which has $n+1$ columns and $n$ rows.

### Basic algorithm

The basic algorithm takes the augmented matrix with arbitrary coefficients and returns an upper triangular augemented matrix. This can then be split into the form of the original system and solved by back substitution.

The algorithm considers each row in turn. Assume 

1. We are looking at row $i$, where $1 \le i \le n$;
2. The previous steps of the algorithm mean that the coefficients in the first $i-1$ columns of row $i$, and all rows below row $i$ are zero. 
3. The $i^{\text{th}}$ entry of row $i$, $a_{i, i}$ is non-zero.

Then we can eliminate all coefficients in column $i$ in the rows below $i$ by subtracting a suitable of row $i$. Explicitly, consider row $j$ where $i < j \le n$. For all coefficients in row $j$, i.e. for $a_{j, k}$ with $1 \le k \le n$, we replace the coefficient by

\begin{equation}
  a_{j, k} \rightarrow a_{j, k} - \frac{a_{j, i}}{a_{i, i}} a_{i, k}.
\end{equation}

We see that when $k = i$ the coefficient $a_{j, i}$ will be set to zero, as required.

If we repeat this algorithm for all rows $i$ then we will have constructed a triangular matrix.

Let us write this algorithm explicitly in code.


```
import numpy as np

def GaussianEliminationReduction(A, b):
    """
    Simple Gaussian Elimination algorithm. Solves linear system A x = b, assuming well-conditioned square matrix A.
    
    Parameters
    ----------
    
    A : float, numpy array, 2d.
        Square matrix, size n x n.
    b : float, numpy array, 1d.
        RHS vector, size n (must conform with A).
    """
    
    # Store size of system
    n = len(b)
    assert(np.all(A.shape == (n, n)))
    
    # Form augmented matrix
    aug = np.hstack((A, b.reshape(n,1)));
    
    # Loop over rows
    for i in range(n):
        # Loop over rows below i
        for j in range(i+1, n):
            aug[j, :] = aug[j, :] - aug[j, i] / aug[i, i] * aug[i, :]
    
    # Return the separated, reduced, matrix and right hand side vector.
    return (aug[:,:-1], aug[:, -1])
```

Now test the algorithm on some simple examples. The first is already in triangular form as above - just checks that it doesn't do anything stupid!


```
A = np.array([[1.0, 2.0, 3.0], \
              [0.0, 4.0, 5.0], \
              [0.0, 0.0, 6.0]])
b = np.array([6.0, 9.0, 6.0])

(A1, b1) = GaussianEliminationReduction(A, b)
print("Resulting triangular matrix is \n{}.\n Right hand side vector becomes {}.\n".format(A1, b1))

C = np.array([[3.0, 0.0, 1.0], \
              [6.0, 2.0, 4.0], \
              [9.0, 2.0, 6.0]])
d = np.array([4.0, 10.0, 15.0])

(A2, b2) = GaussianEliminationReduction(C, d)
print("Resulting triangular matrix is \n{}.\n Right hand side vector becomes {}.\n".format(A2, b2))
```

    Resulting triangular matrix is 
    [[ 1.  2.  3.]
     [ 0.  4.  5.]
     [ 0.  0.  6.]].
     Right hand side vector becomes [ 6.  9.  6.].
    
    Resulting triangular matrix is 
    [[ 3.  0.  1.]
     [ 0.  2.  2.]
     [ 0.  0.  1.]].
     Right hand side vector becomes [ 4.  2.  1.].
    


### Pivoting

However, the simple algorithm above does not work robustly. As a simple example, consider the system below.


```
E = np.array([[1.0, 1.0, 1.0], \
              [0.0, 0.0, 2.0], \
              [0.0, 1.0, 1.0]])
f = np.array([1.0, 1.0, 2.0])

(A3, b3) = GaussianEliminationReduction(E, f)
print("Resulting triangular matrix is \n{}.\n Right hand side vector becomes {}.\n".format(A3, b3))
```

    Resulting triangular matrix is 
    [[  1.   1.   1.]
     [  0.   0.   2.]
     [ nan  nan -inf]].
     Right hand side vector becomes [  1.   1. -inf].
    


    -c:27: RuntimeWarning: divide by zero encountered in double_scalars
    -c:27: RuntimeWarning: invalid value encountered in multiply


The algorithm has failed as the matrix does not satisfy the assumptions necessary: in particular, the element $a_{2, 2}$ is zero.

We can simply avoid this problem by swapping rows. In this instance we would swap the second and third row and the algorithm would work (in this case, swapping the second and third rows puts the matrix in triangular form).

However, there is a similar problem if an entry is very small, but non-zero. Consider the problem

\begin{equation}
  \left( \begin{array}{c c|c} \varepsilon & 1 & 1 \\ 1 & 1 & 2 \end{array} \right) \rightarrow \left( \begin{array}{c c|c} \varepsilon & 1 & 1 \\ 0 & 1 - \varepsilon^{-1} & 2 - \varepsilon^{-1} \end{array} \right)
\end{equation}

which has been reduced to triangular form using Gaussian elimination. If $\varepsilon$ is very small then the terms $\varepsilon^{-1}$ will dominate, making the problem approximately

\begin{equation}
  \left( \begin{array}{c c|c} \varepsilon & 1 & 1 \\ 0 & - \varepsilon^{-1} & - \varepsilon^{-1} \end{array} \right) \quad \implies \quad {\bf x} = \begin{pmatrix} 0 \\ 1 \end{pmatrix}.
\end{equation}

The correct answer *should* be ${\bf x} = (1, 1)^T$. There is nothing wrong with the condition number. The problem is the Gaussian Elimination step where it divides through by $\varepsilon$, a very small number.

To avoid both these problems, the standard technique is to **pivot**, or use **pivoting**. At each stage in the algorithm we are working with row $i$ and looking to eliminate the entries in rows $j > i$. The first step is to check if there is a row $j > i$ with a larger (in magnitude) coefficient in column $i$. If so, we swap the rows. Note that we cannot consider rows above $i$ in the matrix, as these may have non-zero entries in columns $<i$.

We modify the algorithm to include pivoting.


```
def GaussianEliminationPivotingReduction(A, b):
    """
    Simple Gaussian Elimination algorithm with pivoting. Solves linear system A x = b, 
    assuming well-conditioned square matrix A.
    
    Parameters
    ----------
    
    A : float, numpy array, 2d.
        Square matrix, size n x n.
    b : float, numpy array, 1d.
        RHS vector, size n (must conform with A).
    """
    
    # Store size of system
    n = len(b)
    assert(np.all(A.shape == (n, n)))
    
    # Form augmented matrix
    aug = np.hstack((A, b.reshape(n,1)));
    
    # Loop over rows
    for i in range(n):
        # Find the row with largest magnitude, and then swap the rows.
        max_row = np.argmax(aug[i:, i])
        if (max_row): # Only swap rows if the maximum is not this row. \
                      # NOTE: the max_row is counted relative to i, so max_row = 0 => row i.
            tmp               = np.copy(aug[i, :])
            aug[i, :]         = np.copy(aug[i+max_row, :])
            aug[i+max_row, :] = np.copy(tmp)
        # Loop over rows below i
        for j in range(i+1, n):
            aug[j, :] = aug[j, :] - aug[j, i] / aug[i, i] * aug[i, :]
    
    # Return the separated, reduced, matrix and right hand side vector.
    return (aug[:,:-1], aug[:, -1])
```

Check on the previously failing case.


```
(A4, b4) = GaussianEliminationPivotingReduction(E, f)
print("Resulting triangular matrix is \n{}.\n Right hand side vector becomes {}.\n".format(A4, b4))
```

    Resulting triangular matrix is 
    [[ 1.  1.  1.]
     [ 0.  1.  1.]
     [ 0.  0.  2.]].
     Right hand side vector becomes [ 1.  2.  1.].
    


## Special case - tridiagonal systems

Tridiagonal systems turn up frequently in the solution of differential equations, particularly when using finite difference methods.
      
Consider the simple $4 \times 4$ system
      
\begin{equation}
        \begin{pmatrix}
          2 & 1 & 0 & 0 \\
          1 & 2 & 1 & 0 \\
          0 & 1 & 2 & 1 \\
          0 & 0 & 1 & 2
        \end{pmatrix}
        \begin{pmatrix}
          x_1 \\ x_2 \\ x_3 \\ x_4
        \end{pmatrix}
        =
        \begin{pmatrix}
          -1 \\ 0 \\ 0 \\ -1
        \end{pmatrix}.
\end{equation}
      
We proceed as we would with Gaussian elimination: remove the (single) coefficient below the diagonal first, rescaling as we go:
      
\begin{equation}
        \begin{pmatrix}
          1 & \tfrac{1}{2} & 0 & 0 \\
          0 & 1 & \tfrac{2}{3} & 0 \\
          0 & 0 & 1 & \tfrac{3}{4} \\
          0 & 0 & 0 & 1
        \end{pmatrix}
        \begin{pmatrix}
          x_1 \\ x_2 \\ x_3 \\ x_4
        \end{pmatrix}
        =
        \begin{pmatrix}
          -\tfrac{1}{2} \\ \tfrac{1}{3} \\ -\tfrac{1}{4} \\ -\tfrac{3}{5}
        \end{pmatrix},
\end{equation}
      
and then use back-substitution to get
      
\begin{equation}
        \left\{
          \begin{aligned}
            x_4 & =-\tfrac{3}{5} \\
            x_3 & = \tfrac{1}{5} \\
            x_2 & = \tfrac{1}{5} \\
            x_1 & =-\tfrac{3}{5}
          \end{aligned} \right. .
\end{equation}

The general formulation of this is the **Thomas algorithm**. It has the significant advantage that only the three non-trivial diagonals need to be stored and worked with. This can massively reduce the memory required.

38ba29c2a40c1be3384c583659405747c9675491
ipynb
Jupyter Notebook
Lectures/02 - Gaussian Elimination.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-02-12T19:51:52.000Z
2022-03-26T15:34:11.000Z
Lectures/02 - Gaussian Elimination.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2017-05-24T19:49:52.000Z
2018-01-23T21:40:42.000Z
Lectures/02 - Gaussian Elimination.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-01-05T13:30:47.000Z
2022-02-15T09:59:39.000Z
Qwen/Qwen-72B
1. YES
2. YES