8. Simultaneous Linear Equations, Part 1: Row Reduction/Gaussian Elimination¶
Last Revised on March 18, 2021,
adding further uses of Python slicing and “vectorization” in function rowReduction,
eliminating function zeros_below_diagonal in favor of always doing this in function row_reduce,
and fixing a few typos.
References:
Section 2.1.1 Naive Gaussian Elimination of Sauer
Section 6.1 Linear Sysems of Equations of Burden&Faires
Section 7.1 of Chenney&Kincaid
8.1. Introduction¶
The problem of solving a system of \(n\) simultaneous linear equations in \(n\) unknowns, with matrix-vector form \(A x = b\), is quite thoroughly understood as far as having a good general-purpose methods usable with any \(n \times n\) matrix \(A\): essentially, Gaussian elimination (or row-reduction) as seen in most linear algebra courses, combined with some modifications to stay well away from division by zero: partial pivoting. Also, good robust software for this general case is readily available, for example in the Python packages NumPy and SciPy, or in Matlab.
Nevertheless, this basic algorithm can be very slow when \(n\) is large – as it often is when dealing with differential equations (even more so with partial differential equations). We will see that it requires about \(n^3/3\) arithmetic operations.
Thus I will summarise the basic method of row reduction or Gaussian elimination, and then build on it with methods for doing things more robustly, and then with methods for doing it faster in some important special cases:
When one has to solve many systems \(A x^{(m)} = b^{(m)}\) with the same matrix \(A\) but different right-hand side vectors \(b^{(m)}.\)
When \(A\) is banded: most elements are zero, and all the non-zero elements \(a_{i,j}\) are near the main diagonal: \(|i - j|\) is far less than \(n\). (Aside on notation: “far less than” is sometimes denoted \(\ll\), as in \(|i-j| \ll n\).)
When \(A\) is strictly diagonally dominant: each diagonal element \(a_{i,i}\) is larger in magnitude that the sum of the magnitudes of all other elements in the same row.
We might also explore some further topics, perhaps as individual projects:
When \(A\) is positive definite: symmetric (\(a_{i,j} = a_{j,i}\)) and with all eigenvalues positive. This last condition would seem hard to verify, since computing all the eigenvalues of \(A\) is harder that solving \(Ax = b\), but there are important situations where this property is automatically guaranteed, such as with Galerkin and finite-element methods for solving boundary value problems for differential equations.
When \(A\) is sparse: most elements are zero, but not necessarily with all the non-zero elements near the main diagonal.
8.1.1. Python Note 1. Module numpy.linalg with standard nickname la¶
Python package Numpy provides a lot of useful tools for numerical linear algebra through a module numpy.linalg.
Just as package numpy is used so often that there is a conventional nickname np, so numpy.linalg is usually nicknamed la.
import numpy as np
import numpy.linalg as la
# As in recent sections, we import some items from modules individually, so they can be used by "first name only".
from numpy import array, inf, zeros_like, empty
8.2. Strategy for getting from mathematical facts to a good algorithm and then to its implentation in [Python] code¶
Here I take the opportunity to illustrate some useful strategies for getting from mathematical facts and ideas to good algorithms and working code for solving a numerical problem. The pattern we will see here, and often later, is:
8.2.1. Step 1. Get a basic algorithm:¶
Start with mathematical facts (like the equations \(\sum_{j=1}^n a_{ij}x_j = b_i\)).
Solve to get an equation for each unknown — or for an updated aproximation of each unknown — in terms of other quantitities.
Specify an order of evaluation in which all the quantities at right are evaluated earlier.
In this, it is often best to start with a verbal description before specifying the details in more precise and detailed mathematical form.
8.2.2. Step 2. Refine to get a more robust algorithm:¶
Identify cases that can lead to failure due to division by zero and such, and revise to avoid them.
Avoid inaccuracy due to problems like severe rounding error. One rule of thumb is that anywhere that a zero value is a fatal flaw (in particular, division by zero), a very small value is also a hazard when rounding error is present. So avoid very small denominators. (We will soon examine this through the phenomenon of loss of significance, and its extreme case catastrophic cancellation.)
8.2.3. Step 3. Refine to get a more efficient algorithm¶
For example,
Avoid repeated evaluation of exactly the same quantity.
Avoid redundant calculations, such as ones whose value can be determnied in advance; for example, values that can be shown in advance to be zero.
Compare and choose between alternative algorithms.
8.3. Gaussian Elimination, a.k.a. Row Reduction¶
We start by considering the most basic algorithm, based on ideas seen in a linear algebra course.
The problem is best stated as a collection of equations for individual numerical values:
Given coefficients \(a_{i,j} 1 \leq i \leq n,\, 1 \leq j \leq n\) and right-hand side values \(b_i,\, 1 \leq i \leq n\), solve for the \(n\) unknowns \(x_j,\, 1 \leq j \leq n\) in the equations $\( \sum_{j=1}^n a_{i,j} x_j = b_i,\, 1 \leq i \leq n. \)$
In verbal form, the basic strategy of row reduction or Gaussian elimination is this:
Choose one equation and use it to eliminate one chosen unknown from all the other equations, leaving that chosen equation plus \(n-1\) equations in \(n-1\) unknowns.
Repeat recursively, at each stage using one of the remaining equations to eliminate one of the remaining unknowns from all the other equations.
This gives a final equation in just one unknown, preceeded by an equation in that unknown plus one other, and so on: solve them in this order, from last to first.
8.3.1. Determining those choices, to produce a first algorithm: “Naive Gaussian Elimination”¶
A precise algorithm must include rules specifying all the choices indicated above. The simplest “naive” choice, which works in most but not all cases, is to eliminate from the top to bottom and left to right:
Use the first equation to eliminate the first unknown from all other equations.
Repeat recursively, at each stage using the first remaining equation to eliminate the first remaining unknown. Thus, at step \(k\), equation \(k\) is used to eliminate unknown \(x_k\).
This gives one equation in just the last unknown \(x_n\); another equation in the last two unknowns \(x_{n-1}\) and \(x_n\), and so on: solve them in this reverse order, evaluating the unknowns from last to first.
This usually works, but can fail because at some stage the (updated) \(k\)-th equation might not include the \(k\)-th unknown: that is, its coefficient might be zero, leading to division by zero.
We will refine the algorithm to deal with that in the later section Solving Simultaneous Linear Equations, Part 2
8.3.2. Python Note 2. Using Numpy for matrices, vectors and their products in Python 3.5 and beyond¶
As of version 3.5 of Python, vectors, matrices, and their products can be handled very elegantly using Numpy arrays, with the one quirk that the product is denoted by the at-sign @.
That is, for a matrix \(A\) and compatible matrix or vector \(b\) both stored in Numpy arrays, their product is given by A @ b.
This means that, along with my encouragement to totally ignore Python arrays in favor of Numpy arrays, and to usually avoid Python lists when working with numerical data, I also recommend that you ignore the now obsolescent Numpy matrix data type, if you happen to come across it in older material on Numpy.
Aside: Why not A * b? Because that is the more general “point-wise” array product:
c = A * b gives array c with c[i,j] equal to A[i,j] * b[i,j], which is not how matrix multiplication works.
8.4. The general case of solving \(Ax = b\), using Python and NumPy¶
The problem of solving \(Ax = b\) in general, when all you know is that \(A\) is an \(n \times n\) matrix and \(b\) is an \(n\)-vector, can in most cases be handled well by using standard software rather than by writing your own code. Here is an example in Python, solving
using the array type from package numpy and the function solve from the linear algebra module numpy.linalg.
A = array([[4., 2., 7.], [3., 5., -6.],[1., -3., 2.]])
print(f"A =\n{A}")
b = array([2., 3., 4.])
print(f"b = {b}")
A =
[[ 4. 2. 7.]
[ 3. 5. -6.]
[ 1. -3. 2.]]
b = [2. 3. 4.]
8.4.1. Python Note 3.¶
It is important to specify that the entries are real numbers (type “float”); otherwise Numpy does integer arithmetic.
One way to do this is as above: putting a decimal point in the numbers (or to be lazy, in at least one of them!)
Another is to tell the function array that the type is float:
A = array([[4, 2, 7], [3, 5, -6],[1, -3, 2]], dtype=float)
b = array([2, 3, 4], dtype=float)
x = la.solve(A, b)
print("numpy.linalg.solve says that the solution of Ax = b is")
print(f"x = {x}")
# Check the backward error, also known as the residual
r = b - A @ x
print(f"\nAs a check, the residual (or backward error) is")
print(f" r = b-Ax = {r},")
print(f"and its infinity (or 'maximum') norm is ||r|| = {la.norm(r, inf)}")
print("\nAside: another way to compute this is with max(abs(r)):")
print(f"||r|| = {max(abs(r))}")
print(f"and its 1-norm is ||r|| = {la.norm(r, 1)}")
numpy.linalg.solve says that the solution of Ax = b is
x = [ 1.81168831 -1.03246753 -0.45454545]
As a check, the residual (or backward error) is
r = b-Ax = [ 0.0000000e+00 -4.4408921e-16 8.8817842e-16],
and its infinity (or 'maximum') norm is ||r|| = 8.881784197001252e-16
Aside: another way to compute this is with max(abs(r)):
||r|| = 8.881784197001252e-16
and its 1-norm is ||r|| = 1.3322676295501878e-15
8.5. The naive Gaussian elimination algorithm, in pseudo-code¶
Here the elements of the transformed matrix and vector after step \(k\) are named \(a_{i,j}^{(k)}\) and \(b_{k}^{(k)}\), so that the original values are \(a_{i,j}^{(0)} = a_{i,j}\) and \(b_{i}^{(0)} = b_{i}\).
The name \(l_{i,k}\) is given to the multiple of row \(k\) that is subtracted from row \(i\) at step \(k\). This naming might seem redundant, but it becomes very useful later, in the section on LU factorization.
for k from 1 to n-1: \(\qquad\) Step k: get zeros in column k below row k:
\(\quad\) for i from k+1 to n:
\(\qquad\) Evaluate the multiple of row k to subtract from row i:
\(\quad\quad l_{i,k} = a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}\) \(\qquad\) If \(a_{k,k}^{(k-1)} \neq 0\)!
\(\qquad\) Subtract \((l_{i,k}\) times row k) from row i in matrix A …:
\(\quad\quad\) for j from 1 to n:
\(\quad\quad\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)} - l_{i,k} a_{k,j}^{(k-1)}\)
\(\quad\quad\) end for
\(\qquad\) … and at right, subtract \((l_{i,k}\) times \(b_k)\) from \(b_i\):
\(\quad\quad b_i^{(k)} = b_i^{(k-1)} - l_{i,k} b_{k}^{(k-1)}\)
\(\quad\) end for
The rows before \(i=k\) are unchanged, so they are ommited from the update;
however, in a situation where we need to complete the definitions of \(A^{(k)}\) and \(b^{(k)}\) we would also need the following inside the for k loop:
\(\quad\) for i from 1 to k:
\(\quad\quad\) for j from 1 to n:
\(\quad\quad\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)}\)
\(\quad\quad\) end for
\(\quad\quad b_i^{(k)} = b_i^{(k-1)}\)
\(\quad\) end for
However, the algorithm will usually be implemented by overwriting the previous values in an array with new ones, and then this part is redundant.
The next improvement in efficiency: the updates in the first \(k\) columns at step \(k\) give zero values (that is the key idea of the algorithm!), so there is no need to compute or store those zeros, and thus the only calculations needed in the above for j from 1 to n loop are covered by for j from k+1 to n.
Thus from now on we use only the latter: except when, for demonstration purposes, we need those zeros.
Thus, the standard algorithm looks like this:
for k from 1 to n-1: \(\qquad\) Step k: Get zeros in column k below row k:
\(\quad\) for i from k+1 to n: \(\qquad\) Update only the rows that change: from k+1 on:
\(\qquad\) Evaluate the multiple of row k to subtract from row i:
\(\quad\quad l_{i,k} = a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}\) \(\qquad\) If \(a_{k,k}^{(k-1)} \neq 0\)!
\(\qquad\) Subtract \((l_{i,k}\) times row k) from row i in matrix A, in the columns that are not automaticaly zero:
\(\quad\quad\) for j from k+1 to n:
\(\quad\quad\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)} - l_{i,k} a_{k,j}^{(k-1)}\)
\(\quad\quad\) end for
\(\qquad\) and at right, subtract \((l_{i,k}\) times \(b_k)\) from \(b_i\):
\(\quad\quad b_i^{(k)} = b_i^{(k-1)} - l_{i,k} b_{k}^{(k-1)}\)
\(\quad\) end for
8.5.1. Python Note 4. Syntax for for loops and 0-based array indexing¶
Since array indices in Python (and in Java, C, C++, C#, Swift, etc.) start from zero, not from one, it will be convenient to express linear algebra algorithms in a form compatible with this.
Every index is one less than in the above! Thus in an array with \(n\) elements, the index values \(i\) are \(0 \leq i < n\), excluding n, which is the half-open interval of integers \([0, n)\).
In the indexing of an array, one can refer to the part the array with indices \(a \leq i < b\), excluding b, with the slice notation
a:b.Similarly, when specifiying the range of consecutive integers \(i\), \(a \leq i < b\) in a
forloop, one can use the expressionrange(a,b).
Also, when indices are processed in order (from low to high), these notes will abuse notation slightly, refering to the values as a set — specifically, a semi-open interval of integers.
For example, the above loop
for j from k+1 to n:
first gets all indices lowered by one, to
for j from k to n-1:
and then this will sometimes be described in terms of the set of j values:
for j in [k,n):
which in Python becomes
for j in range(k, n):
This new notation needs care initially, but helps with clarity in the long run.
For one thing, it means that the indices of an \(n\)-element array, \([0,n-1)\), are described by range(0,n) and by 0:n.
In fact, the case of “starting at the beginning”, with index zero, can be abbreviated: range(n) is the same as range(0,n), and :b ia the same as 0:b.
Another advantage is that the index ranges a:b and b:c together cover the same indices as a:c, with no gap or duplication of b, and likewise range(a,b) and range(b,c) combine to cover range(a,c).
8.6. The naive Gaussian elimination algorithm, in Pythonic zero-based pseudo-code¶
Here the above notational shift is made, along with eliminating the above-noted redundant formulas for values that are either zero or are unchanged from the previous step. It is also convenient for \(k\) to be the index of the row being used to reduce subsequent rows, and so also the index of the column in which values below the main diagonal are being set to zero.
for k in [0, n-1):
\(\quad\) for i in [k+1, n):
\(\quad\quad l_{i,k} = a_{i,k}^{(k)}/a_{k,k}^{(k)}\qquad\) If \(a_{k,k}^{(k)} \neq 0\)!
\(\quad\quad\) for j in [k+1, n):
\(\quad\quad\quad a_{i,j}^{(k+1)} = a_{i,j}^{(k)} - l_{i,k} a_{k,j}^{(k)}\)
\(\quad\quad\) end for
\(\quad\quad b_i^{(k+1)} = b_i^{(k)} - l_{i,k} b_{k}^{(k)}\)
\(\quad\) end for
end for
8.7. The naive Gaussian elimination algorithm, in Python¶
Conversion to actual Python code is now quite straightforward; there is litle more to be done than:
Change the way that indices are described, from \(b_i\) to
b[i]and from \(a_{i,j}\) toA[i,j].Use case consistently in array names, since the quirk in mathematical notation of using upper-case letters for matrix names but lower case letters for their elements is gone! In these notes, matrix names will be upper-case and vector names will be lower-case (even when a vector is considered as 1-column matrix).
Rather than create a new array for each matrix \(A^{(0)}\), \(A^{(0)}\), etc. and each vector \(b^{(0)}\), \(b^{(1)}\), we overwite each in the same array.
Aside: We will see that this simplicity in translation is quite common once algorithms have been expressed with zero-based indexing. The main ugliness is with loops that count backwards; see below.
for k in range(n-1):
for i in range(k+1, n):
L[i,k] = A[i,k] / A[k,k]
for j in range(k+1, n):
A[i,j] -= L[i,k] * A[k,j]
b[i] -= L[i,k] * b[k]
To demonstrate this, some additions are needed:
Putting this algorithm into a function.
Getting the value \(n\) needed for the loop, using the fact that it is the length of vector
b.Creating the array \(L\).
Copying the input arrays
Aandbinto new ones,Uandc, so that the original arrays are not changed. That is, when the row reduction is completed,Ucontains \(A^{(n-1)}\) andccontains \(b^{(n-1)}\).
Also, for some demonstrations, the zero values below the main diagonal of U are inserted, though usually they would not be needed.
def rowReduce(A, b):
"""To avoid modifying the matrix and vector specified as input,
they are copied to new arrays, with the method .copy()
Warning: it does not work to say "U = A" and "c = b";
this makes these names synonyms, referring to the same stored data.
"""
U = A.copy()
c = b.copy()
n = len(b)
# The function zeros_like() is used to create L with the same size and shape as A,
# and with all its elements zero initially.
L = np.zeros_like(A)
for k in range(n-1):
for i in range(k+1, n):
# compute all the L values for column k:
L[i,k] = U[i,k] / U[k,k] # Beware the case where U[k,k] is 0
for j in range(k+1, n):
U[i,j] -= L[i,k] * U[k,j]
# Put in the zeros below the main diagonal in column k of U;
# this is not important for calculations, since those elements of U are not used in backward substitution,
# but it helps for displaying results and for checking the results via residuals.
U[i,k] = 0.
c[i] -= L[i,k] * c[k]
return (U, c)
def zeros_below_diagonal(U):
"""Insert the below-diagonal zero values into U, ignored in the function rowReduce.
These are needed to display U correctly,
and to demonstrate that the new system of equations Ux=c$\ has the same solution as Ax=b.
"""
for i in range(1,len(U)):
for j in range(i):
U[i,j] = 0
Note: As usual, you could omit the above defs and instead import these functions with
from numerical_methods_module import rowReduce
#from numerical_methods_module import zeros_below_diagonal
(U, c) = rowReduce(A, b)
#zeros_below_diagonal(U)
print(f"U =\n{U}")
print(f"c = {c}")
U =
[[ 4. 2. 7. ]
[ 3. 3.5 -11.25]
[ 1. -3.5 -11. ]]
c = [2. 1.5 5. ]
Let’s take advantage of the fact that we have used la.solve to get a very accurate approxiamti oof teh solution x
of \(Ax=b\); this should also solve \(Ux=c\), so chej the backward error, a.k.a. the residual:
r = c - U@x
print(f"\nThe residual (backward error) c-Ux is {r}, with maximum norm {max(abs(r))}.")
The residual (backward error) c-Ux is [ 0. -5.43506494 -5.42532468], with maximum norm 5.4350649350649345.
8.7.1. Python Note 5. Operations on a sequence of array indices, with “slicing”: vectorization¶
Python code can specify vector operations on a range of indices \([c,d)\), referred to withthe slice notaiton c:d.
For example, the slice notation A[c:d,j] refers to the array containing the \(d-c\) elements A[i,j] for \(i\) in the semi-open interval \([c,d)\).
Thus, each of the three arithmetic calculations above can be specified over a range of index values in a single command, eliminating all the inner-most for loops;
this is somtimes called vectorization.
Only for loops that contains other for loops remain.
Apart from mathematical elegance, this usually allows far faster execution.
for k in range(n-1):
L[k+1:n,k] = U[k+1:n,k] / U[k,k] # compute all the L values for column k
for i in range(k+1, n):
U[i,k+1:n] -= L[i,k] * U[k,k+1:n] # Update row i
c[k+1:n] -= L[k+1:n,k] * c[k] # update c values
I will break my usual guideline by redefining rowReduce, since this is just a different statement of exactly the same algorithm:
def rowReduce(A, b):
"""To avoid modifying the matrix and vector specified as input,
they are copied to new arrays, with the method .copy()
Warning: it does not work to say "U = A" and "c = b";
this makes these names synonyms, referring to the same stored data.
"""
U = A.copy()
c = b.copy()
n = len(b)
# The function zeros_like() is used to create L with the same size and shape as A,
# and with all its elements zero initially.
L = zeros_like(A)
for k in range(n-1):
# compute all the L values for column k:
L[k+1:n,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0
for i in range(k+1, n):
U[i,k+1:n] -= L[i,k] * U[k,k+1:n] # Update row i
# Insert the below-diagonal zeros in column k;
# this is not important for calculations, since those elements of U are not used in backward substitution,
# but it helps for displaying results and for checking the results via residuals.
U[i,k] = 0.0
c[k+1:n] -= L[k+1:n,k] * c[k] # update c values
return (U, c)
UPDATE: this was further updated on March 7 as follows, mainly to add demoMode.
def rowReduce(A, b, demoMode=False):
"""To avoid modifying the matrix and vector specified as input,
they are copied to new arrays, with the method .copy()
Warning: it does not work to say "U = A" and "c = b";
this makes these names synonyms, referring to the same stored data.
2021-03-07: added a demonstration mode.
"""
U = A.copy()
c = b.copy()
n = len(b)
# The function zeros_like() is used to create L with the same size and shape as A,
# and with all its elements zero initially.
L = np.zeros_like(A)
for k in range(n-1):
if demoMode: print(f"Step {k=}")
# compute all the L values for column k:
L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0
if demoMode:
print(f"The multipliers in column {k+1} are {L[k+1:,k]}")
for i in range(k+1, n):
U[i,k+1:n] -= L[i,k] * U[k,k+1:n] # Update row i
# Insert the below-diagonal zeros in column k;
# this is not important for calculations, since those elements of U are not used in backward substitution,
# but it helps for displaying results and for checking the results via residuals.
U[i,k] = 0.0
c[k+1:n] -= L[k+1:n,k] * c[k] # update c values
if demoMode:
# insert zeros in U:
U[k+1:, k] = 0.
print(f"The updated matrix is\n{U}")
print(f"The updated right-hand side is\n{c}")
return (U, c)
UPDATE: this was updated yet again on March 17, with further vectorization, eliminating the for loop over i.
def rowReduce(A, b, demoMode=False):
"""To avoid modifying the matrix and vector specified as input,
they are copied to new arrays, with the method .copy()
Warning: it does not work to say "U = A" and "c = b";
this makes these names synonyms, referring to the same stored data.
2021-03-07: added a demonstration mode.
2021-03-17: added further "vectorization", eliminating the loop over i
"""
U = A.copy()
c = b.copy()
n = len(b)
# The function zeros_like() is used to create L with the same size and shape as A,
# and with all its elements zero initially.
L = np.zeros_like(A)
for k in range(n-1):
if demoMode: print(f"Step {k=}")
# compute all the L values for column k:
L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0
if demoMode:
print(f"The multipliers in column {k+1} are {L[k+1:,k]}")
# In the following, the index "k:k+1" refers to the single value k,
# but is needed to make it clear that
# - L[k+1:n,k:k+1] is a matrix with one column (not a vector), and
# - U[k:k+1,k+1:n] is a matrix with one row (not a vector),
# to get the appropriate matrix multiplicaton
U[k+1:n,k+1:n] -= L[k+1:n,k:k+1] @ U[k:k+1,k+1:n]
# Insert the below-diagonal zeros in column k;
# this is not important for calculations, since those elements of U are not used in backward substitution,
# but it helps for displaying results and for checking the results via residuals.
U[k+1:n,k] = 0.0
c[k+1:n] -= L[k+1:n,k] * c[k] # update c values
if demoMode:
U[k+1:n, k] = 0. # insert zeros in couln k of U:
print(f"The updated matrix is\n{U}")
print(f"The updated right-hand side is\n{c}")
return (U, c)
def zeros_below_diagonal(U):
"""Insert the below-diagonal zero values into U, ignored in the function rowReduce.
These are needed to display U correctly,
and to demonstrate that the new system of equations Ux=c has the same solution as Ax=b.
"""
for i in range(1,len(U)):
U[i,:i] = 0.
Repeating the above testing:
(U, c) = rowReduce(A, b, demoMode=True)
#zeros_below_diagonal(U)
r = c - U@x
print(f"U =\n{U}")
print(f"c = {c}")
print(f"\nThe residual (backward error) c - Ux is {r}, with maximum norm {max(abs(r))}.")
Step k=0
The multipliers in column 1 are [0.75 0.25]
The updated matrix is
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. -3.5 0.25]]
The updated right-hand side is
[2. 1.5 3.5]
Step k=1
The multipliers in column 2 are [-1.]
The updated matrix is
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. -11. ]]
The updated right-hand side is
[2. 1.5 5. ]
U =
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. -11. ]]
c = [2. 1.5 5. ]
The residual (backward error) c - Ux is [0. 0. 0.], with maximum norm 0.0.
8.8. Backward substitution with an upper triangular matrix¶
The transformed equations have the form
and can be solved from bottom up, starting with \(x_n = c_n/u_{n,n}\).
All but the last equation can be written as
and so solved as
This procedure is backward substitution, giving the algorithm
\(x_n = c_n/u_{n,n}\)
for i from n-1 down to 1
\(\displaystyle \quad x_i = \frac{c_i - \sum_{j=i+1}^{n} u_{i,j} x_j}{u_{i,i}}\)
end for
This works so long as none of the main diagonal terms \(u_{i,i}\) is zero, because when done in this order, everything on the right hand side is known by the time it is evaluated.
For future reference, note that the elements \(u_{k,k}\) that must be non-zero here, the ones on the main diagonal of \(U\), are the same as the elements \(a_{k,k}^{(k)}\) that must be non-zero in the row reduction stage above, because after stage \(k\), the elements of row \(k\) do not change any more: \(a_{k,k}^{(k)} = a_{k,k}^{(n-1)} = u_{k,k}\).
8.9. The backward substitution algorithm in zero-based pseudo-code¶
Again, a zero-based version is more convenient for programming in Python (or Java, or C++):
\(x_{n-1} = c_{n-1}/u_{n-1,n-1}\)
for i from n-2 down to 0
\(\displaystyle \quad x_i = \frac{c_i - \sum_{j=i+1}^{n-1} u_{i,j} x_j}{u_{i,i}}\)
end for
8.9.1. Python Note 6. Indexing from the end of an array and counting backwards¶
To express the above backwards counting in Python, we have to deal with the fact that range(a,b) counts upwards and excludes the “end value” b.
The first part is easy: the extended form range(a, b, step) increments by step instead of by one, so that range(a, b, 1) is the same as range(a,b), and range(a, b, -1) counts down: \(a, a-1, \dots, b+1\).
But it still stops just before \(b\), so getting the values from \(n-1\) down to \(0\) requires using \(b= -1\), and so the slightly quirky expression range(n-1, -1, -1).
One more bit of Python: for an \(n\)-element single-index array v, the sum of its elements \(\sum_{i=0}^{n-1} v_i\) is given by sum(v).
Thus \(\sum_{i=a}^{b-1} v_i\), the sum over a subset of indices \([a,b)\), is given by sum(v[a:b]).
And remember that multiplication of Numpy arrays with * is pointwise.
8.9.2. The backward substitution algorithm in Python¶
With all the above Python details, the core code for backward substitution is:
x[n-1] = c[n-1]/U[n-1,n-1]
for i in range(n-2, -1, -1):
x[i] = (c[i] - sum(U[i,i+1:] * x[i+1:])) / U[i,i]
Aside/preview: Note that the backward substitution algorithm and its Python coding have a nice mathematical advantage over the row reduction algorithm above: the precise mathematical statement of the algorithm does not need any intermediate quantities distinguished by superscripts \({}^{(k)}\), and correspondingly, all variables in the code have fixed meanings, rather than changing at each step.
In other words, all uses of the equal sign are mathematically correct as equations!
This can be advantageous in creating algorithms and code that is more understandable and more readily verified to be correct, and is an aspect of the functional programming approach. We will soon go part way to that functional ideal, by rephrasing Gaussian elimination in a form where all variables have clear, fixed meanings, corresponding to the natural mathematical description of the process: the method of LU factorization to be seen in Section 8.1 of the text Numerical Mathematics and Computing.
8.9.3. Python Note 7. Another way to count backwards along an array¶
On the other hand, there is an elegant way access array elements “from the top down”.
Firstly (or “lastly”) x[-1] is the last element: the same as x[n-1] when n = len(x), but without needing to know that length \(n\).
More generally, x[-i] is x[n-i].
Thus, one possibly more elegant way to describe backward substitution is to count with an increasing index, the “distance from the bottom”:
from x[n-1] which is x[-1] to x[0], which is x[-n].
That is, index -i replaces index \(n - i\):
x[-1] = c[-1]/U[-1,-1]
for i in range(2, n+1):
x[-i] = (c[-i] - sum(U[-i,1-i:] * x[1-i:])) / U[-i,-i]
There is still the quirk of having to “overshoot”, referring to n+1 in range to get to final index -n.
As a final demonstration, we put this second version of the code into a complete working Python function and test it:
def backwardSubstitution(U, c):
x = np.zeros_like(c)
x[-1] = c[-1]/U[-1,-1]
for i in range(2, len(x) + 1):
x[-i] = (c[-i] - sum(U[-i,1-i:] * x[1-i:])) / U[-i,-i]
return x
UPDATE: this was updated on March 7 as follows, mainly to add demoMode.
def backwardSubstitution(U, c, demoMode=False):
"""Solve U x = c for b.
2021-03-07: aded a demonstration mode.
"""
n = len(c)
x = np.zeros(n)
x[-1] = c[-1]/U[-1,-1]
if demoMode: print(f"x_{n} = {x[-1]}")
for i in range(2, n+1):
x[-i] = (c[-i] - sum(U[-i,1-i:] * x[1-i:])) / U[-i,-i]
if demoMode: print(f"x_{n-i+1} = {x[-i]}")
return x
also available via
from numerical_methods_module import backwardSubstitution
x = backwardSubstitution(U, c)
print(f"x = {x}")
r = b - A@x
print(f"\nThe residual b - Ax = {r},")
print(f"with maximum norm {max(abs(r)):.3}.")
x = [ 1.81168831 -1.03246753 -0.45454545]
The residual b - Ax = [ 0.0000000e+00 -4.4408921e-16 8.8817842e-16],
with maximum norm 8.88e-16.
8.10. Two code testing hacks: starting from a known solution, and using randomly generated examples¶
An often useful strategy in developing and testing code is to create a test case with a known solution; another is to use random numbers to avoid accidently using a test case that in unusually easy.
# Prefered style is to have all "import" statements at the top,
# but since this is the first time we've heard of module 'random',
# I did not want it to be mentioned mysteriously above.
import random
x_random = empty(len(b)) # An array the same length as b, with no values specified yet
for i in range(len(x)):
x_random[i] = random.uniform(-1, 1) # gives random real value, from uniform distribution in [-1, 1]
print(f"x_random = {x_random}")
x_random = [-0.12411689 -0.15377245 0.42636405]
Create a right-hand side b that automatically makes x_random the correct solution:
b_random = A @ x_random
print(f"A =\n{A}")
print(f"\nb_random = {b_random}")
(U, c_random) = rowReduce(A, b_random)
#zeros_below_diagonal(U)
print(f"\nU=\n{U}")
print(f"\nResidual c_random - U@x_random = {c_random - U@x_random}")
x_computed = backwardSubstitution(U, c_random)
print(f"\nx_computed = {x_computed}")
print(f"\nResidual b_random - A@x_computed = {b_random - A@x_computed}")
print(f"\nBackward error |b_random - A@x_computed| = {max(abs(b_random - A@x_computed))}")
print(f"\nError x_random - x_computed = {x_random - x_computed}")
print(f"\nAbsolute error |x_random - x_computed| = {max(abs(x_random - x_computed))}")
A =
[[ 4. 2. 7.]
[ 3. 5. -6.]
[ 1. -3. 2.]]
b_random = [ 2.18053587 -3.69939722 1.18992855]
U=
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. -11. ]]
Residual c_random - U@x_random = [0. 0. 0.]
x_computed = [-0.12411689 -0.15377245 0.42636405]
Residual b_random - A@x_computed = [0. 0. 0.]
Backward error |b_random - A@x_computed| = 0.0
Error x_random - x_computed = [0. 0. 0.]
Absolute error |x_random - x_computed| = 0.0
8.11. What can go wrong? Three examples¶
8.11.1. Example 1¶
8.11.1.1. An obvious division by zero problem¶
Consider the system of two equations
It is easy to see that this has the solution \(x_1 = x_2 = 1\); in fact it is already in “reduced form”. However when put into matrix form
the above algorithm fails, because the fist pivot element \(a_{11}\) is zero:
A1 = array([[0., 1.], [1. , 1.]])
b1 = array([1., 1.])
(U1, c1) = rowReduce(A1, b1)
#zeros_below_diagonal(U1)
print(f"U1 = \n{U1}")
print(f"c1 = {c1}")
x1 = backwardSubstitution(U1, c1)
print(f"x1 = {x1}")
U1 =
[[ 0. 1.]
[ 0. -inf]]
c1 = [ 1. -inf]
x1 = [nan nan]
<ipython-input-12-a7d463b36315>:20: RuntimeWarning: divide by zero encountered in true_divide
L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0
/Users/lemesurierb/OneDrive - College of Charleston/UNCo-MATH375/elementary-numerical-analysis-python/python_notebooks_modules/numerical_methods_module.py:388: RuntimeWarning: invalid value encountered in double_scalars
x[-1] = c[-1]/U[-1,-1]
Python Note 8.
inf, meaning “infinity”, is a special value given as the result of operations like division by zero. Surprisingly, it can have a sign! (This is available in Python from package Numpy asnumpy.inf)nan, meaning “not a number”, is a special value given as the result of calculations like0/0. (This is available in Python from package Numpy asnumpy.nan)
8.11.2. Example 2¶
8.11.2.1. A less obvious division by zero problem¶
Next consider this system
The solution is \(x_1 = x_2 = x_3 = 1\), and this time none of th diagonal elements is zero, so it is not so obvoius that a divisin be zero probelms will occur, but:
A2 = array([[1., 1., 1.], [1., 1., 2.],[1., 2., 2.]])
b2 = array([3., 4., 5.])
(U2, c2) = rowReduce(A2, b2)
#zeros_below_diagonal(U2)
print(f"U2 = \n{U2}")
print(f"c2 = {c2}")
x2 = backwardSubstitution(U2, c2)
print(f"x2 = {x2}")
U2 =
[[ 1. 1. 1.]
[ 0. 0. 1.]
[ 0. 0. -inf]]
c2 = [ 3. 1. -inf]
x2 = [nan nan nan]
<ipython-input-12-a7d463b36315>:20: RuntimeWarning: divide by zero encountered in true_divide
L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0
What happens here is that the first stage subtracts the first row from each of the others …
A2[1,:] -= A2[0,:]
b2[1] -= b2[0]
A2[2,:] -= A2[0,:]
b2[2] -= b2[0]
… and the new matrix has the same problem as above at the next stage:
print(f"Now A2 is \n{A2}")
print(f"and b2 is {b2}")
Now A2 is
[[1. 1. 1.]
[0. 0. 1.]
[0. 1. 1.]]
and b2 is [3. 1. 2.]
Thus, the second and third equations are
with the same problem as in Example 1.
8.11.3. Example 3¶
8.11.3.1. Problems caused by inexact arithmetic¶
The equations
again have the solution \(x_1 = x_2 = 1\), and the only division that happens in the above algorithm for row reduction is by that pivot element \(a_{11} = 1, \neq 0\), so with exact arithmetic, all would be well. But:
A3 = array([[1., 1e16], [1. , 1.]])
b3 = array([1. + 1e16, 2.])
print(f"A3 = \n{A3}")
print(f"b3 = {b3}")
A3 =
[[1.e+00 1.e+16]
[1.e+00 1.e+00]]
b3 = [1.e+16 2.e+00]
(U3, c3) = rowReduce(A3, b3)
#zeros_below_diagonal(U3)
print(f"U3 = \n{U3}")
print(f"c3 = {c3}")
x3 = backwardSubstitution(U3, c3)
print(f"x3 = {x3}")
U3 =
[[ 1.e+00 1.e+16]
[ 0.e+00 -1.e+16]]
c3 = [ 1.e+16 -1.e+16]
x3 = [2. 1.]
This gets \(x_2 = 1\) correct, but \(x_1\) is completely wrong!
One hint is that \(b_1\), which should be \(1 + 10^{16} = 1000000000000001\), is instead just given as \(10^{16}\).
On the other hand, all is well with less large values, like \(10^{15}\):
A3a = array([[1., 1e15], [1. , 1.]])
b3a = array([1. + 1e15, 2.])
print(f"A3a = \n{A3a}")
print(f"b3a = {b3a}")
A3a =
[[1.e+00 1.e+15]
[1.e+00 1.e+00]]
b3a = [1.e+15 2.e+00]
(U3a, c3a) = rowReduce(A3a, b3a)
#zeros_below_diagonal(U3a)
print(f"U3a = \n{U3a}")
print(f"c3a = {c3a}")
x3a = backwardSubstitution(U3a, c3a)
print(f"x3a = {x3a}")
U3a =
[[ 1.e+00 1.e+15]
[ 0.e+00 -1.e+15]]
c3a = [ 1.e+15 -1.e+15]
x3a = [1. 1.]
8.11.4. Example 4¶
8.11.4.1. Avoiding small denominators¶
The first equation is Example 3 can be divided by \(10^{16}\) to get an equivalent system with the same problem:
Now the problem is more obvious: this system differs from the system in Example 1 just by a tiny change of \(10^{-16}\) in that pivot elements \(a_{11}\), and the problem is division by a value very close to zero.
A4 = array([[1e-16, 1.], [1. , 1.]])
b4 = array([1. + 1e-16, 2.])
print(f"A4 = \n{A4}")
print(f"b4 = {b4}")
A4 =
[[1.e-16 1.e+00]
[1.e+00 1.e+00]]
b4 = [1. 2.]
(U4, c4) = rowReduce(A4, b4)
#zeros_below_diagonal(U4)
print(f"U4 = \n{U4}")
print(f"c4 = {c4}")
x4 = backwardSubstitution(U4, c4)
print(f"x4 = {x4}")
U4 =
[[ 1.e-16 1.e+00]
[ 0.e+00 -1.e+16]]
c4 = [ 1.e+00 -1.e+16]
x4 = [2.22044605 1. ]
One might think that there is no such small denominator in Example 3, but what counts for being “small” is magnitude relative to other values — 1 is very small compared to \(10^{16}\).
To understand these problems more (and how to avoid them) it is time to explore Machine Numbers, Rounding Error and Error Propagation .
This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International