Gaussian Elimination

This article discusses the Gaussian elimination algorithm, one of the most fundamental and important numerical algorithms of all time. It is used to solve linear equations of the form

: $\mathbf{A} \vec{x} = \vec{b},$

where $\mathbf{A}$ is a known $N\times N$ matrix, $\vec{b}$ is a known vector of length $N$ , and $\vec{x}$ is an unknown vector of length $N$ . The goal is to find $\vec{x}$ . The Gaussian elimination algorithm is implemented by Scipy's scipy.linalg.solve function.

There are many good expositions of Gaussian elimination on the web (including the Wikipedia article). Here, we will only go through the basic aspects of the algorithm; you are welcome to read up elsewhere for more details.

The basic algorithm

The best way to understand how Gaussian elimination works is to work through a concrete example. Consider the following $N=3$ problem:

: $\begin{bmatrix}1 &2 &3 \\ 3 &2 &2 \\ 2 &6 &2\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix}3 \\ 4 \\ 4\end{bmatrix}.$

The Gaussian elimination algorithm consists of two distinct phases: row reduction and back-substitution.

Row reduction

In the row reduction phase of the algorithm, we manipulate the matrix equation so that the matrix becomes upper triangular (i.e., all the entries below the diagonal are zero). To achieve this, we note that we can subtract one row from another row, without altering the solution. In fact, we can subtract any multiple of a row.

We will eliminate (zero out) the elements below the diagonal in a specific order: from top to bottom along each column, then from left to right for successive columns. For our $3 \times 3$ example, the elements that we intend to eliminate, and the order in which we will eliminate them, are indicated by the colored numbers 0, 1, and 2 in the following figure:

The first matrix element we want to eliminate is at $(1,0)$ (orange circle). To eliminate it, we subtract, from this row, a multiple of row $0$ . We will use a factor of $3/1 = 3$ :

: $(3x_0 + 2x_1 + 2x_2) - (3/1)(1x_0 + 2x_1 + 3x_2) = 4 - (3/1) 3$

The factor of $3$ we used is determined as follows: we divide the matrix element at $(1,0)$ (which is the one we intend to eliminate) by the element at $(0,0)$ (which is the one along the diagonal in the same column). As a result, the term proportional to $x_0$ disappears, and we obtain the following modified linear equations, which possess the same solution:

(Note that we have changed the entry in the vector on the right-hand side as well, not just the matrix on the left-hand side!) Next, we eliminate the element at $(2,0)$ (green circle). To do this, we subtract, from this row, a multiple of row $0$ . The factor to use is $2/1 = 2$ , which is the element at $(2,0)$ divided by the $(0,0)$ (diagonal) element:

: $(2x_0 + 6x_1 + 2x_2) - (2/1)(1x_0 + 2x_1 + 3x_2) = 4 - (2/1) 3$

The result is

Next, we eliminate the $(2,1)$ element (blue circle). This element lies in column $1$ , so we eliminate it by subtracting a multiple of row $1$ . The factor to use is $2/(-4) = -0.5$ , which is the $(2,1)$ element divided by the $(1,1)$ (diagonal) element:

: $(0x_0 + 2x_1 - 4x_2) - (2/(-4))(0x_0 - 4x_1 - 7x_2) = -5 - (2/(-4)) (-2)$

The result is

We have now completed the row reduction phase, since the matrix on the left-hand side is upper-triangular (i.e., all the entries below the diagonal have been set to zero).

Back-substitution

In the back-substitution phase, we read off the solution from the bottom-most row to the top-most row. First, we examine the bottom row:

Thanks to row reduction, all the matrix elements on this row are zero except for the last one. Hence, we can read off the solution

: $x_2 = (-4.5)/(-7.5) = 0.6.$

Next, we look at the row above:

This is an equation involving $x_1$ and $x_2$ . But from the previous back-substitution step, we know $x_2$ . Hence, we can solve for

: $x_1 = [-5 - (-7) (0.6)] / (-4) = 0.2.$

Finally, we look at the row above:

This involves all three variables $x_0$ , $x_1$ , and $x_2$ . But we already know $x_1$ and $x_2$ , so we can read off the solution for $x_0$ . The final result is

: $\begin{bmatrix}x_0 \\ x_1 \\ x_2\end{bmatrix} = \begin{bmatrix}0.8 \\ 0.2 \\ 0.6\end{bmatrix}.$

Runtime

Let's summarize the components of the Gaussian elimination algorithm, and analyze how many steps each part takes:

Row reduction
•	Step forwards through the rows. For each row $n$ ,			$N$ steps
	•	Perform pivoting (to be discussed below).		$N-n +1 \sim N$ steps
	•	Step forwards through the rows larger than $n$ . For each such row $m$ ,		$N-n \sim N$ steps
		•	Subtract $(A'_{mn}/A'_{nn})$ times row $n$ from the row $m$ (where $\mathbf{A}'$ is the current matrix). This eliminates the matrix element at $(m,n)$ .	$O(N)$ arithmetic operations
Back-substitution
•	Step backwards through the rows. For each row $n$ ,			$N$ steps
	•	Substitute in the solutions $x_m$ for $m > n$ (which are already found). Hence, find $x_n$ .		$N-n \sim O(N)$ arithmetic operations

(The "pivoting" procedure hasn't been discussed yet; we'll do that in a later section.)

We conclude that the runtime of the row reduction phase scales as $O(N^3)$ , and the runtime of the back-substitution phase scales as $O(N^2)$ . The algorithm's overall runtime therefore scales as $O(N^3)$ .

Matrix generalization

We can generalize the Gaussian elimination algorithm described in the previous section, to solve matrix problems of the form

: $\mathbf{A}\, \mathbf{x} = \mathbf{b}$ ,

where $\mathbf{x}$ and $\mathbf{b}$ are $N\times M$ matrices, not merely vectors. An example, for $M = 2$ , is

: $\begin{bmatrix}1 &2 &3 \\ 3 &2 &2 \\ 2 &6 &2\end{bmatrix} \begin{bmatrix} x_{00} & x_{01} \\ x_{10} & x_{11} \\ x_{20} & x_{21} \end{bmatrix} = \begin{bmatrix}3 & 6 \\ 4 & 8 \\ 4 & 2\end{bmatrix}.$

It can get a bit tedious to keep writing out the $x$ elements in the system of equations, particularly when $\mathbf{x}$ becomes a matrix. For this reason, we switch to a notation known as the augmented matrix:

: $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 6 \\ 3 & 2 & 2 & 4 & 8 \\ 2 & 6 & 2 & 4 & 2 \end{array}\right].$

Here, the entries to the left of the vertical separator denote the left-hand side of the system of equations, and the entries to the right of the separator denote the right-hand side of the system of equations.

The Gaussian elimination algorithm can now be performed directly on the augmented matrix. We will walk through the steps for the above example. First, row reduction:

Eliminate the element at $(1,0)$ : $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 6 \\ 0 & -4 & -7 & -5 & -10 \\ 2 & 6 & 2 & 4 & 2 \end{array}\right]$
Eliminate the element at $(2,0)$ : $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 6 \\ 0 & -4 & -7 & -5 & -10 \\ 0 & 2 & -4 & -2 & -10 \end{array}\right]$
Eliminate the element at $(2,1)$ : $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 6 \\ 0 & -4 & -7 & -5 & -10 \\ 0 & 0 & -7.5 & -4.5 & -15 \end{array}\right]$

The back-substitution step converts the left-hand portion of the augmented matrix to the identity matrix:

Solve for row $2$ : $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 3 \\ 0 & -4 & -7 & -5 & -10 \\ 0 & 0 & 1 & 0.6 & 2 \end{array}\right]$
Solve for row $1$ : $\left[\begin{array}{ccc|cc} 1 &\; 2 & \;\;\;3 \;\;& 3 & \;3 \\ 0 & \;1 & \;\;\;0\;\; & 0.2 & \;-1 \\ 0 & \;0 & \;\;\;1\;\; & 0.6 & \;2 \end{array}\right]$
Solve for row $0$ : $\left[\begin{array}{ccc|cc} 1 & \;0 & \;\;\;0\;\; & 0.8 & \;2 \\ 0 & \;1 & \;\;\;0\;\; & 0.2 & \;-1 \\ 0 & \;0 & \;\;\;1\;\; & 0.6 & \;2 \end{array}\right]$

After the algorithm finishes, the right-hand side of the augmented matrix contains the result for $\mathbf{x}$ . Analyzing the runtime using the same reasoning as before, we find that the row reduction step scales as $O\Big(N^2(N+M)\Big)$ , and the back-substitution step scales as $O\Big(N(N+M)\Big)$ .

This matrix form of the Gaussian elimination algorithm is the standard method for computing matrix inverses. If $\mathbf{b}$ is the $N\times N$ identity matrix, then the solution $\mathbf{x}$ will be the inverse of $\mathbf{A}$ . Thus, the runtime for calculating a matrix inverse scales as $O(N^3)$ .

Pivoting

In our description of the Gaussian elimination algorithm so far, you may have noticed a problem. During the row reduction process, we have to multiply rows by a factor of $A'_{mn}/A'_{nn}$ (where $\mathbf{A}'$ denotes the current matrix). If $A'_{nn}$ happens to be zero, the factor blows up, and the algorithm fails.

To bypass this difficulty, we add an extra step to the row reduction procedure. As we step forward through row numbers $n = 0, 1, \cdots, N-1$ , we do the following for each $n$ :

Pivoting: search through the matrix elements on and below the diagonal element at $(n,n)$ , and find the row $n'$ with the largest value of $|A'_{n'n}|$ . Then, swap row $n$ and row $n'$ .
Continue with the rest of the algorithm, eliminating the $A'_{mn}$ elements below the diagonal.

You should be able to convince yourself that (i) pivoting does not alter the solution, and (ii) it does not alter the runtime scaling of the row reduction phase, which remains $O(N^3)$ .

Apart from preventing the algorithm from failing unnecessarily, pivoting improves its numerical stability. If $A'_{nn}$ is non-zero but very small in magnitude, dividing by it will produce a very large result, which brings about a loss of floating-point numerical precision. Hence, it is advantageous to swap rows around to ensure that the magnitude of $A'_{nn}$ is as large as possible.

When trying to pivot, it might happen that all the values of $|A'_{n'n}|$ , on and below the diagonal, are zero (or close enough to zero within our floating-point tolerance). If this happens, it indicates that our original $\mathbf{A}$ matrix is singular, i.e., it has no inverse. Hence, the pivoting procedure has the additional benefit of helping us catch the cases where there is no valid solution to the system of equations; in such cases, the Gaussian elimination algorithm should abort.

Example

Let's work through an example of Gaussian elimination with pivoting, using the problem in the previous section:

: $\left[\begin{array}{ccc|cc} 1 & 2 & 3 & 3 & 6 \\ 3 & 2 & 2 & 4 & 8 \\ 2 & 6 & 2 & 4 & 2 \end{array}\right].$

The row reduction phase goes as follows:

( $n = 0$ ): Pivot, swapping row $0$ and row $1$ : $\left[\begin{array}{ccc|cc} 3 & 2 & 2 & 4 & 8 \\ 1 & 2 & 3 & 3 & 6 \\ 2 & 6 & 2 & 4 & 2 \end{array}\right].$
( $n = 0$ ): Eliminate the element at $(1,0)$ : $\left[\begin{array}{ccc|cc} 3 & 2 & 2 & 4 & 8 \\ 0 & 4/3 & 7/3 & 5/3 & 10/3 \\ 2 & 6 & 2 & 4 & 2 \end{array}\right].$
( $n = 0$ ): Eliminate the element at $(2,0)$ : $\left[\begin{array}{ccc|cc} 3 & 2 & 2 & 4 & 8 \\ 0 & 4/3 & 7/3 & 5/3 & 10/3 \\ 0 & 14/3 & 2/3 & 4/3 & -10/3 \end{array}\right].$
( $n = 1$ ): Pivot, swapping row $1$ and row $2$ : $\left[\begin{array}{ccc|cc} 3 & 2 & 2 & 4 & 8 \\ 0 & 14/3 & 2/3 & 4/3 & -10/3 \\0 & 4/3 & 7/3 & 5/3 & 10/3 \end{array}\right].$
( $n = 1$ ): Eliminate the element at $(2,1)$ : $\left[\begin{array}{ccc|cc} 3 & 2 & 2 & 4 & 8 \\ 0 & 14/3 & 2/3 & 4/3 & -10/3 \\0 & 0 & 15/7 & 9/7 & 30/7 \end{array}\right].$

The back-substitution phase then proceeds as usual. You can check that it gives the same results we obtained before.

LU decomposition

A variant of the Gaussian elimination algorithm can be used to compute the LU decomposition of a matrix. This procedure was invented by Alan Turing, the British mathematician considered the "father of computer science". The LU decomposition of a square matrix $\mathbf{A}$ consists of a lower-triangular matrix $\mathbf{L}$ and an upper-triangular matrix $\mathbf{U}$ , such that

: $\mathbf{A} = \mathbf{L}\, \mathbf{U}.$

In certain special circumstances, LU decompositions provide a very efficient method for solving linear equations. Suppose that we have to solve a set of linear equations $\mathbf{A} \vec{x} = \vec{b}$ many times, using the same $\mathbf{A}$ but an indefinite number of $\vec{b}$ 's which might not be known in advance. For example, the $\vec{b}$ 's might represent an endless series of measurement outcomes, with $\mathbf{A}$ representing some fixed experimental configuration. We would like to efficiently calculate $\vec{x}$ for each $\vec{b}$ that arrives. If this is done with Gaussian elimination, each calculation would take $O(N^3)$ time.

However, if we can perform an LU decomposition ahead of time, then the calculations can be performed much more quickly. The linear equations are

: $\mathbf{L}\, \mathbf{U} \vec{x} = \vec{b}.$

This can be broken up into two separate equations:

: $\mathbf{L}\, \vec{y} = \vec{b}, \quad \mathrm{and}\;\;\; \mathbf{U} \vec{x} = \vec{y}.$

Because $\mathbf{L}$ is lower-triangular, we can solve the first equation by forward-substitution (similar to back-substitution, except that it goes from the first row to last) to find $\vec{y}$ . Then we can solve the second equation by back-substitution, to find $\vec{x}$ . The whole process takes $O(N^2)$ time, which is a tremendous improvement over performing a wholesale Gaussian elimination.

However, finding the LU decomposition takes $O(N^3)$ time (we won't go into details here, but it's basically a variant of the row reduction phase of the Gaussian elimination algorithm). Therefore, if we are interested in solving the linear equations only once, or a handful of times, the LU decomposition method does not improve performance. It's useful in situations where the LU decomposition is performed ahead of time. You can think of the LU decomposition as a way of re-arranging the Gaussian elimination algorithm, so that we don't need to know $\vec{b}$ during in the first, expensive $O(N^3)$ phase.

In Python, you can perform the LU decomposition using the scipy.linalg.lu function. The forward-substitution and back-substitution steps can be performed using scipy.linalg.solve_triangular.

Computational physics notes | Prev: Numerical linear algebra Next: Eigenvalue problems