1.4 Gaussian Elimination Revisited

1Reading¶

Material related to this page, as well as additional exercises, can be found in LAA Ch. 2.5, ALA Ch 1.3, and ILA Ch. 2.6. This page is mostly based on ALA Ch 1.3.

2Learning Objectives¶

By the end of this page, you should know:

how to use Gaussian elimination to solve linear systems $A\vv x = \vv b$ when $A$ is a regular matrix
what pivots in a matrix are
algorithms to solve large systems of linear equations using Gaussian elimination and back substitution

3Gaussian Elimination: Regular Case¶

With basic matrix arithmetic operations in our toolkit, we will develop a systematic method for solving linear systems of equations. For a linear system $A\vv x = \vv b$ , with $A$ an $m\times n$ coefficient matrix, $\vv x$ an $n \times 1$ unknowns vector, and $\vv b$ an $m \times 1$ right hand side vector, we define the augmented matrix:

M = \left[\begin{array}{c|c} A & \vv b \end{array}\right] =\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \end{array}\right],

(1)

which is an $m \times (n+1)$ matrix obtained by tacking the right-hand side vector $\vv b$ onto the right of the coefficient matrix $A$ . The extra vertical line is just to remind us that the last column of this matrix plays a special role. For example, the augmented matrix for our example system is

M = \left[ \begin{array}{ccc|c} 1 & 2 & 1 & 2\\ 2 & 6 & 1 & 7 \\ 1 & 1 & 4 & 3 \end{array}\right]

(2)

Note that it is simple to go back and forth between the original linear system and the augmented matrix, but since operations on equations also affect their right-hand sides, it is convenient to keep track of everything together using the augmented matrix.

For the time being, we will concentrate our efforst on linear systems that have the same number, $n$ , of equations as unknowns. The associated coefficient matrix $A$ is square of size $n \times n$ , and the corresponding augmented matrix $M = [ A \, | \, \vv b]$ then has size $n \times (n+1)$ .

We start with a simple observation connecting Linear System Operation #1 to its equivalent matrix operation

For example, when solving our example system, our first step was to subtract two times the first equation from the second. This is equivalently done by subtracting two times the first row of the augmented matrix (2) from the second row:

-2\bm 1 & 2 & 1 & 2 \em + \bm 2 & 6 & 1 & 7 \em = \bm 0 & 2 & -1 & 3\em.

(3)

We recognize this as the second row of the modified augmented matrix

\left[ \begin{array}{ccc|c} 1 & 2 & 1 & 2\\ 0 & 2 & -1 & 3\\ 1 & 1 & 4 & 3 \end{array}\right],

(4)

that corresponds to the first equivalent example system. When elementary row operation #1 is performed, it is critical that the result replaces the row being added to and not the row being multiplied by the scalar. Notice that the elimination of a variable in an equation, in this case the first variable in the second equation, amounts to making its entry in the coefficient matrix equal to zero.

3.1Pivots¶

We will call the $(1,1)$ entry of the coefficient matrix the first pivot. The precise definition of a pivot will become clear as we continue, but one key requirement is that a pivot must always be nonzero. Eliminating the first variable $x_1$ from the second and third equations is the same as making all of the matrix entries in the column below the pivot equal to zero. We have already done this with the $(2,1)$ entry in (4). To make the $(3,1)$ entry equal to zero, we subtract the first from from the last row, resulting in the augmented matrix

\left[ \begin{array}{ccc|c} 1 & 2 & 1 & 2\\ 0 & 2 & -1 & 3\\ 0 & -1 & 3 & 1 \end{array}\right],

(5)

which we again recognize as the corresponding to the second equivalent example system. The second pivot is the $(2,2)$ entry of this matrix, which is 2, and is the coefficient of the second variable $x_2$ in the second equation. Again, the pivot must be nonzero. We use the Observation 1 of adding $1/2$ of the second row to the third row to make the entry below the second pivot equal to 0, resulting in the augmented matrix

\left[ \begin{array}{ccc|c} 1 & 2 & 1 & 2\\ 0 & 2 & -1 & 3\\ 0 & 0 & \frac{5}{2} & \frac{5}{2} \end{array}\right],

(6)

that corresponds to the triangular system equivalent to our example system. We write the final augmented matrix as

N = [U \, | \, \vv c], \quad \text{where} \quad U = \bm 1 & 2 & 1 \\ 0 & 2 & -1 \\ 0 & 0 & \frac{5}{2}\em, \quad \vv c = \bm 2 \\ 3 \\ \frac{5}{2} \em.

(7)

The corresponding linear system can be written as $U\vv x = \vv c$ . A special feature of this system is that the coefficient matrix $U$ is upper triangular^[1], which means that all entries below the main diagonal are zero, i.e., $u_{ij}=0$ whenever $i>j$ . The three nonzero entries on its diagonal, 1, 2, and $5/2$ , including the last one in the $(3,3)$ slot, are the three pivots. Once the system has been reduced to this triangular form, we can easily solve it via Back Substitution.

What we just described is an algorithm for solving a linear system of $n$ equations in $n$ unknowns, and is known as regular Gaussian Elimination. We’ll call a square matrix $A$ regular if the algorithm successfully reduces it to the upper triangular form $U$ with all nonzero pivots on the diagonal. If this fails to happen, i.e., if a pivot appearing on the diagonal is zero, then the matrix is not regular. We then use the pivot row to make all entries lying in the column below the pivot equal to zero through elementary row operations. The solution is then found by applying Back Substitution to the resulting system. We present both of these algorithms in pseudocode and python code below.

Here we use what are called in place updates, meaning that the same letter $M$ (with entries $m_{ij}$ ) denotes the current augmented matrix at each stage in the computation. We initialize with $M=[A \, | \, \vv b]$ , and output (assuming $A$ is regular) the upper triangular equivalent augmented matrix $M = [U \, | \, \vv c]$ , where $U$ is the upper triangular matrix with diagonal entries the pivots, and $\vv c$ is the resulting vector of the right-hand sides of the triangular system $U\vv x = \vv c$ .

Python Break!¶

Here’s what this looks like implemented as a function in NumPy. You should be able to map the code below to the pseudo-code above.

import numpy as np

def GaussElim(A, b):
    n = A.shape[0] # number of rows of A = number of equations
    M = np.hstack((A, b.reshape((n,1)))) # build the augmented matrix by horizontally stacking A and b.
    
    # Gaussian elimination
    for j in range(n):
        if M[j, j] == 0:
            print("A is not regular")
            break
        else:
            for i in range(j+1, n):
                scalar = M[i, j]/M[j, j]
                M[i, :] = M[i, :] - scalar*M[j, :]
    
    # return the matrix M = [U | c] in upper triangular form
    return M        

# Let's test our function
A = np.array([[1., 2., 1.],
              [2., 6., 1.],
              [1., 1., 4.]])

b = np.array([2, 7, 3])

M = GaussElim(A,b)
print(M)

[[ 1.   2.   1.   2. ]
 [ 0.   2.  -1.   3. ]
 [ 0.   0.   2.5  2.5]]

# The importance of a period!
Aint = np.array([[1, 2, 1],
              [2, 6, 1],
              [1, 1, 4]])
bint = np.array([2, 7, 3])

newM = GaussElim(Aint,bint)
print(M-newM) # these are different!

[[0.  0.  0.  0. ]
 [0.  0.  0.  0. ]
 [0.  0.  0.5 0.5]]

Next, let’s take a look at the pseudocode for Back Substitution.

Here’s what this looks like implemented as a function in NumPy. You should be able to map the code below to the pseudo-code above.

def back_substitution(M):
    # We assume M = [U|c] has U in upper triangular form
    # Better software engineering practices would call for us to check this before proceeding but we
    # focus on the math here
    
    n = M.shape[0] # number of rows of M = number of equations = number of unknowns for square A
    x = np.zeros((n,)) # set up an all zeros array x to be populated with the solution
    
    x[n-1] = M[n-1, -1]/M[n-1, n-1] # -1 is the last index of the particular axis in a python array
    for i in range(n-2, -1, -1): # range(a, b, c) returns the list [a, a+c ,a+2c, ..., b] 
        x[i] = (M[i, -1] - np.sum(M[i, i+1:-1]*x[i+1:]))/M[i, i] # * is elementwise multiplication of arrays
    
    return x

# Let's test out our code using M from above: it gives us the right answer!
x = back_substitution(M)
print(x)

[-3.  2.  1.]

Now let’s put the two functions together to define a solve function that takes in a matrix A and a right hand side b and returns a solution x such that A@x = b.

def solve(A,b):
    M = GaussElim(A,b)
    x = back_substitution(M)
    return x

x = solve(A,b)

print(f'x = {x}\n A @ x - b = {A @ x - b}')

x = [-3.  2.  1.]
 A @ x - b = [0. 0. 0.]

3.2Worked Examples¶

Exercise 1 (Regular Gaussian elimination and backsubstitution)

For the following system, write its augmented matrix. By hand, execute the regular Gaussian elimination algorithm; if the algorithm succeeds, then solve the system via backsubstitution on the upper triangular form augmented matrix.

\begin{align*} 2x + y - z &= 8\\ -3x - y + 2z &= -11\\ -2x + y + 2z &= -3 \end{align*}

(8)

Solution to Exercise 1

The augmented matrix is:

\begin{align*} \left[ \begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ -3 & -1 & 2 & -11 \\ -2 & 1 & 2 & -3 \end{array} \right] \end{align*}

(9)

Our first pivot is the element in position (1, 1), which is 2, which is nonzero. Let’s use the first row to eliminate the first elements in rows 2 and 3, by doing $R_2 \gets R_2 + (3/2)R_1$ and $R_3 \gets R_3 + R_1$ (where $R_i$ denotes the ith row of the augmented matrix). Doing this gives us the equivalent augmented matrix:

\begin{align*} \left[ \begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 1/2 & 1/2 & 1 \\ 0 & 2 & 1 & 5 \end{array} \right] \end{align*}

(10)

Our second pivot is the element in position (2, 2), which is $1/2$ , which is nonzero. Using the second row to eliminate the second element in row 3, we get:

\begin{align*} \left[ \begin{array}{ccc|c} 2 & 1 & -1 & 8 \\ 0 & 1/2 & 1/2 & 1 \\ 0 & 0 & -1 & 1 \end{array} \right] \end{align*}

(11)

The third and last pivot is -1, which is nonzezro, meaning the regular Gaussian elimination algorithm succeeded.

Now, let’s solve for $x$ , $y$ , and $z$ using backsubstitution. Starting from the last row:

The last row $\left[ \begin{array}{ccc|c} 0 & 0 & -1 & 1 \end{array} \right]$ tells us that $-z = 1$ , in other words, $z = -1$ .
The second to last row tells us $y/2 + z/2 = 1$ . Since we know $z = -1$ , we have $y/2 -1/2 = 1$ , meaning $y = 3$ .
The final row tells us that $2x + y - z = 8$ . Since we know $z = -1$ and $y = 3$ , we find that $x = 2$ .

So our final answer is $(x, y, z) = (2, 3, -1)$ .

Exercise 2 (Regular Gaussian elimination and backsubstitution (cont.))

For the following system, write its augmented matrix. By hand, execute the regular Gaussian elimination algorithm. Show how the algorithm fails for this system.

\begin{align*} x + 2y + 3z &= 6\\ 2x + 4y + 6z &= 12\\ x + y + z &= 3 \end{align*}

(12)

Solution to Exercise 2

The augmented matrix is:

\begin{align*} \left[ \begin{array}{ccc|c} 1 & 2 & 3 & 6 \\ 2 & 4 & 6 & 12 \\ 1 & 1 & 1 & 3 \end{array} \right] \end{align*}

(13)

Our first pivot is the element in position (1, 1), which is 1, which is nonzero. Let’s use the first row to eliminate the first elements in rows 2 and 3, by doing $R_2 \gets R_2 - 2R_1$ and $R_3 \gets R_3 - R_1$ (where $R_i$ denotes the ith row of the augmented matrix). Doing this gives us the equivalent augmented matrix:

\begin{align*} \left[ \begin{array}{ccc|c} 1 & 2 & 3 & 6 \\ 0 & 0 & 0 & 0 \\ 0 & -1 & -2 & -3 \end{array} \right] \end{align*}

(14)

Our second pivot is the element in position (2, 2), which is 0, which is zero. Therefore regular Gaussian elimination fails, meaning our coefficient matrix is irregular.

4Solving Big Systems of Linear Equations¶

The code that we wrote above works equally well when we have 1000 equations in 1000 unknowns (or more!). The code below generates a random A and b of size 1000 x 1000 and 1000 x 1, respectively, and solves for a 1000 x 1 solution x satisfying A @ x = b. Compare how long our code takes to solve this problem and the built in function from NumPy np.linalg.solve.

If you want to play around with this, launch the interactive notebook and see how big of a system you can solve before it takes too long or you run out of memory!

# Setup a large random system of equation Ax = b
n = 1000 # number of equations and unknowns

# Generate random A and b of size (n,n) and (n,1)
A_big = np.eye(n) + np.random.randn(n,n)
b_big = np.random.randn(n)

# solve linear system of equations Ax=b using our homebrewed solution
x = solve(A_big, b_big)

# print max absolute value of Ax - b
print(np.max(np.abs(A_big @ x - b_big)))

4.564420619246334e-09

# compare to using NumPy's built in function
x_np = np.linalg.solve(A_big, b_big)

# print max absolute value of Ax - b
print(np.max(np.abs(A_big @ x_np - b_big)))

3.06916714265526e-11

# use some python magic to compare timing: which should you use in practice?!
%timeit solve(A_big, b_big)
%timeit np.linalg.solve(A_big, b_big)

805 ms ± 4.06 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.83 ms ± 641 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)