Implementing Gaussian Elimination#
We start with an implementation of Gaussian elimination as described previously.
Caution
A feature of this implementation is that the input A and b are changed by this routine, and on output they reflect the row-echelon form. This is done to save memory.
We will print the steps along the way
import numpy as np
def gauss_elim(A, b, *, quiet=False, pivot=True):
""" perform gaussian elimination with pivoting, solving A x = b.
A is an NxN matrix, x and b are an N-element vectors. Note: A
and b are changed upon exit to be in upper triangular (row
echelon) form """
assert b.ndim == 1, "ERROR: b should be a vector"
N = len(b)
assert A.shape == (N, N), "ERROR: A should be square with each dim of same length as b"
x = np.zeros((N), dtype=A.dtype)
if not quiet:
print_Ab(A, b)
# main loop over rows
for k in range(N):
if not quiet:
print(f"working on row {k}")
if pivot:
# find the pivot row based on the size of column k -- only consider
# the rows >= k (then add k to the index so it is 0-based)
row_max = np.argmax(np.abs(A[k:, k])) + k
if row_max != k:
# swap the row with the largest element in the current column
# with the current row (pivot) -- do this with b too!
A[[k, row_max], :] = A[[row_max, k], :]
b[[k, row_max]] = b[[row_max, k]]
if not quiet:
print("pivoted, updated system:")
print_Ab(A, b)
# do the forward-elimination for all rows below the current
for i in range(k+1, N):
coeff = A[i, k] / A[k, k]
for j in range(k+1, N):
A[i, j] += -A[k, j] * coeff
A[i, k] = 0.0
b[i] += -coeff * b[k]
# check if the row is all zeros -- singular
if np.abs(A[i, :]).max() == 0:
if not quiet:
print("singular")
print_Ab(A, b)
raise ValueError("matrix is singular")
if not quiet:
print_Ab(A, b)
# back-substitution
# last solution is easy
x[N-1] = b[N-1] / A[N-1, N-1]
for i in reversed(range(N-1)):
bsum = b[i]
for j in range(i+1, N):
bsum += -A[i, j] * x[j]
x[i] = bsum / A[i, i]
return x
This routine prints the augmented matrix \(({\bf A} | {\bf b})\).
def print_Ab(A, b):
"""printout the matrix A and vector b in a pretty fashion."""
N = len(b)
space = 8*" "
top_str = "⎧" + N*" {:>7.03f} " + "⎫" + space + "⎧" + " {:6.3f} " + "⎫"
sid_str = "⎪" + N*" {:>7.03f} " + "⎪" + space + "⎪" + " {:6.3f} " + "⎪"
bot_str = "⎩" + N*" {:>7.03f} " + "⎭" + space + "⎩" + " {:6.3f} " + "⎭"
for i in range(N):
if i == 0:
pstr = top_str
elif i == N-1:
pstr = bot_str
else:
pstr = sid_str
out = tuple(A[i, :]) + (b[i],)
print(pstr.format(*out))
print(" ")
Now we can test this out.
Important
Be sure to create your arrays as floating point arrays. If we don’t specify, and just include integers in the initialization, then the array will be an integer array, and we will truncate all divisions and won’t get the right answer.
A = np.array([[1, 1, 1],
[-1, 2, 0],
[2, 0, 1]], dtype=np.float64)
b = np.array([6, 3, 5], dtype=np.float64)
# since our method changes A and b as it works, we'll send in copies
x = gauss_elim(A.copy(), b.copy())
⎧ 1.000 1.000 1.000 ⎫ ⎧ 6.000 ⎫
⎪ -1.000 2.000 0.000 ⎪ ⎪ 3.000 ⎪
⎩ 2.000 0.000 1.000 ⎭ ⎩ 5.000 ⎭
working on row 0
pivoted, updated system:
⎧ 2.000 0.000 1.000 ⎫ ⎧ 5.000 ⎫
⎪ -1.000 2.000 0.000 ⎪ ⎪ 3.000 ⎪
⎩ 1.000 1.000 1.000 ⎭ ⎩ 6.000 ⎭
⎧ 2.000 0.000 1.000 ⎫ ⎧ 5.000 ⎫
⎪ 0.000 2.000 0.500 ⎪ ⎪ 5.500 ⎪
⎩ 0.000 1.000 0.500 ⎭ ⎩ 3.500 ⎭
working on row 1
⎧ 2.000 0.000 1.000 ⎫ ⎧ 5.000 ⎫
⎪ 0.000 2.000 0.500 ⎪ ⎪ 5.500 ⎪
⎩ 0.000 0.000 0.250 ⎭ ⎩ 0.750 ⎭
working on row 2
⎧ 2.000 0.000 1.000 ⎫ ⎧ 5.000 ⎫
⎪ 0.000 2.000 0.500 ⎪ ⎪ 5.500 ⎪
⎩ 0.000 0.000 0.250 ⎭ ⎩ 0.750 ⎭
print(x)
[1. 2. 3.]
To ensure that we got the right answer, we can compare \({\bf A}{\bf x}\) and \({\bf b}\). We can use the python matrix multiplication operator, @:
print(A @ x - b)
[0. 0. 0.]
\(4\times 4\) example#
Let’s try another, larger example
A = np.array([[1., 2., 3., 4.],
[5., 1., 6., -1.],
[10., 2., 13., -2.],
[4., 10., -2., -5.]])
b = np.array([3., 2., 4., 1.])
x = gauss_elim(A.copy(), b.copy())
⎧ 1.000 2.000 3.000 4.000 ⎫ ⎧ 3.000 ⎫
⎪ 5.000 1.000 6.000 -1.000 ⎪ ⎪ 2.000 ⎪
⎪ 10.000 2.000 13.000 -2.000 ⎪ ⎪ 4.000 ⎪
⎩ 4.000 10.000 -2.000 -5.000 ⎭ ⎩ 1.000 ⎭
working on row 0
pivoted, updated system:
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 5.000 1.000 6.000 -1.000 ⎪ ⎪ 2.000 ⎪
⎪ 1.000 2.000 3.000 4.000 ⎪ ⎪ 3.000 ⎪
⎩ 4.000 10.000 -2.000 -5.000 ⎭ ⎩ 1.000 ⎭
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 0.000 -0.500 0.000 ⎪ ⎪ 0.000 ⎪
⎪ 0.000 1.800 1.700 4.200 ⎪ ⎪ 2.600 ⎪
⎩ 0.000 9.200 -7.200 -4.200 ⎭ ⎩ -0.600 ⎭
working on row 1
pivoted, updated system:
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 9.200 -7.200 -4.200 ⎪ ⎪ -0.600 ⎪
⎪ 0.000 1.800 1.700 4.200 ⎪ ⎪ 2.600 ⎪
⎩ 0.000 0.000 -0.500 0.000 ⎭ ⎩ 0.000 ⎭
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 9.200 -7.200 -4.200 ⎪ ⎪ -0.600 ⎪
⎪ 0.000 0.000 3.109 5.022 ⎪ ⎪ 2.717 ⎪
⎩ 0.000 0.000 -0.500 0.000 ⎭ ⎩ 0.000 ⎭
working on row 2
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 9.200 -7.200 -4.200 ⎪ ⎪ -0.600 ⎪
⎪ 0.000 0.000 3.109 5.022 ⎪ ⎪ 2.717 ⎪
⎩ 0.000 0.000 0.000 0.808 ⎭ ⎩ 0.437 ⎭
working on row 3
⎧ 10.000 2.000 13.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 9.200 -7.200 -4.200 ⎪ ⎪ -0.600 ⎪
⎪ 0.000 0.000 3.109 5.022 ⎪ ⎪ 2.717 ⎪
⎩ 0.000 0.000 0.000 0.808 ⎭ ⎩ 0.437 ⎭
x
array([0.47186147, 0.18181818, 0. , 0.54112554])
A @ x - b
array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 6.66133815e-16])
Singular example#
changing A[2, 2] from 13 to 12 makes the system singular, since row 2 is just double row 1.
A = np.array([[1., 2., 3., 4.],
[5., 1., 6., -1.],
[10., 2., 12., -2.],
[4., 10., -2., -5.]])
b = np.array([3., 2., 4., 1.])
x = gauss_elim(A, b)
⎧ 1.000 2.000 3.000 4.000 ⎫ ⎧ 3.000 ⎫
⎪ 5.000 1.000 6.000 -1.000 ⎪ ⎪ 2.000 ⎪
⎪ 10.000 2.000 12.000 -2.000 ⎪ ⎪ 4.000 ⎪
⎩ 4.000 10.000 -2.000 -5.000 ⎭ ⎩ 1.000 ⎭
working on row 0
pivoted, updated system:
⎧ 10.000 2.000 12.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 5.000 1.000 6.000 -1.000 ⎪ ⎪ 2.000 ⎪
⎪ 1.000 2.000 3.000 4.000 ⎪ ⎪ 3.000 ⎪
⎩ 4.000 10.000 -2.000 -5.000 ⎭ ⎩ 1.000 ⎭
singular
⎧ 10.000 2.000 12.000 -2.000 ⎫ ⎧ 4.000 ⎫
⎪ 0.000 0.000 0.000 0.000 ⎪ ⎪ 0.000 ⎪
⎪ 1.000 2.000 3.000 4.000 ⎪ ⎪ 3.000 ⎪
⎩ 4.000 10.000 -2.000 -5.000 ⎭ ⎩ 1.000 ⎭
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[11], line 7
1 A = np.array([[1., 2., 3., 4.],
2 [5., 1., 6., -1.],
3 [10., 2., 12., -2.],
4 [4., 10., -2., -5.]])
5 b = np.array([3., 2., 4., 1.])
----> 7 x = gauss_elim(A, b)
Cell In[2], line 53, in gauss_elim(A, b, quiet, pivot)
51 print("singular")
52 print_Ab(A, b)
---> 53 raise ValueError("matrix is singular")
55 if not quiet:
56 print_Ab(A, b)
ValueError: matrix is singular
Roundoff example#
Let’s try our example with roundoff both with and without pivoting
eps = 1.e-15
A = np.array([[eps, 1.0, 1.0],
[-1.0, 2.0, 0.0],
[2.0, 0.0, 1.0]])
b = np.array([6.0, 3.0, 5.0])
x = gauss_elim(A.copy(), b.copy(), quiet=True)
x
array([0.33333333, 1.66666667, 4.33333333])
x = gauss_elim(A.copy(), b.copy(), quiet=True, pivot=False)
x
array([1.77635684, 2. , 4. ])
we see that with \(\epsilon \rightarrow 0\), without pivoting, we get a very poor solution.