2.4. Solving \(Ax = b\) with LU factorization#
Last revised on October 28, 2025, adding Algorithm 2.9 for forward substitution with upper triangular matrices that are not “unit”, Algorithm 2.10 for the Cholesky factorization of positive definite matrices, and expanded notes on when LU factorization works.
References:
[Chasnov, 2012] Section 3.2, LU decomposition.
[Sullivan, 2021] Section 4.4, The LU Factorization.
[Sauer, 2022] Section 2.2, The LU Factorization.
[Burden et al., 2016] Section 6.5, Matrix Factorizations[Burden et al., 2016].
[Dionne, 2023] Sections 4.2, LU Factorization and 4.3, Cholesky Factorization. (Section 4.2 also includes pivoting, as seen in the section Solving Ax = b With Both Pivoting and LU factorization.)
[Chenney and Kincaid, 2013] Section 8.1, Matrix Factorizations.
[Kincaid and Chenney, 1990] Section 4.2, LU and Cholesky Factorizations.
[Trefethen and Bau, 2022] Lecture 20, Gaussian Elimination.
[Dahlquist and Björck, 2009] Section 7.2.3, Pivoting Strategies.
# Import some items from modules individually, so they can be used by "first name only".
from numpy import array, zeros_like, identity
2.4.1. Avoiding repeated calculation, excessive rounding and messy notation: the LU factorization (a.k.a. LU decomposition)#
Putting aside pivoting for a while, there is another direction in which the algorithm for solving linear systems \(Ax=b\) can be improved. It starts with the idea of being more efficient when solving multiple systems with the same right-hand side: \(A x^{(m)} = b^{(m)}, m=1,2, \dots\).
However it has several other benefits:
allowing a strategy to reduce rounding error, and
a simpler, more elegant mathematical statement.
We will see how to merge this with partial pivoting in the section Solving Ax = b With Both Pivoting and LU factorization.
Some useful jargon:
Definition 2.7 (Triangular matrix)
A matrix is triangular if all its non-zero elements are either on the main diagonal or on only one side of it. There are two possibilities:
Matrix \(U\) is upper triangular if \(u_{ij} = 0\) for all \(i > j\).
Matrix \(L\) is lower triangular if \(l_{ij} = 0\) for all \(j > i\).
If in addition, the elements of the main diagonal of such a matrix are all 1’s, it is called unit triangular.
One important example of an upper triangular matrix is \(U\) formed by row reduction; note well that it is much quicker and easier to solve \(Ux = c\) than the original system \(Ax=b\) exactly because of its triangular form.
We will soon see that the multipliers \(l_{ij}\), \(i > j\) for row reduction that were introduced in the section Row Reduction (Gaussian Elimination) help to form a very useful unit lower triangular matrix \(L\).
The key to the LU factorization idea is finding a lower triangular matrix \(L\) and an upper triangular matrix \(U\) such that \(L U = A\), and then using the fact that it is far quicker to solve a linear system when the corresponding matrix is triangular.
Indeed it will be shown in Theorem 2.6 that if naive row reduction for \(Ax=b\) succeeds, giving row-reduced form \(Ux = c\):
The matrix \(A\) can be factorized as \(A = LU\) with \(U\) an \(n \times n\) upper triangular matrix and \(L\) an \(n \times n\) lower triangular matrix.
With the further condition that \(L\) is unit lower triangular, there is a unique such factorization: this is called the Doolittle factorization of \(A\).
In the Doolittle factorization, the matrix \(U\) is the one given by naive row reduction, and the elements of \(L\) below its main diagonal are the multipliers arising in naive row reduction. (The other elements of \(L\), on and above the main diagonal, are the ones and zeros dictated by it being unit lower triangular: the same as for those elements in the \(n \times n\) identity matrix.)
The transformed right-hand side \(c\) arising from naive row reduction is the solution of \(Lc = b\), and this is solvable by a procedure caled forward substitution, very similar to the backward subsitution used to solve \(Ux = c\).
Putting all this together: if naive row reduction works for \(A\), we can introduce the name \(c\) for \(Ux\), and note that \(Ax = (LU)x = L(Ux) = Lc = b\). Then solving of the system \(Ax = b\) can be done in three steps:
Using \(A\), find the Doolittle factors, \(L\) and \(U\).
Using \(L\) and \(b\), solve \(Lc = b\) to get \(c\). (Forward substitution)
Using \(U\) and \(c\), solve \(Ux = c\) to get \(x\). (Backward substitution)
2.4.2. The direct method for the Doolittle LU factorization#
If you believe the above claims, we already have one algorithm for finding an LU factorization; basically, do naive row reduction, but ignore the right-hand side \(b\) until later. However, there is another “direct” method, which does not rely on anything we have seen before about row reduction, and has other advantages as we will see.
(Aside: If I were teaching linear algebra, I would be tempted to start here and skip row reduction!)
This method starts by considering the apparently daunting task of solving the \(n^2\) simultaneous and nonlinear equations for the initially unknown elements of \(L\) and \(U\):
The first step is to insert the known information; the already-known values of elements of \(L\) and \(U\). For one thing, the sums above stop when either \(k=i\) or \(k=j\), whichever comes first, due to all the zeros in \(L\) nd \(U\):
Next, when \(i \leq j\) — so that the sum ends at \(k=i\) and involves \(l_{i,i}\) — we can use \(l_{i,i} = 1\).
So break up into two cases:
On and above the main diagonal (\(i \leq j\), so \(\min(i,j) = i\)):
Below the main diagonal (\(i > j\), so \(\min(i,j) = j\)):
In each equation, the last term in the sum has been separated, so that we can use them to “solve” for an unknown:
Here comes the characteristic step that gets us from valid equations to a useful algorithm: we can arrange these equations in an order such that all the values at right are determined by an earlier equation!
First look at what they say for the first row and first column.
With \(i=1\) in the first equation, there is no sum, and so:
which is the familiar fact that the first row is unchanged in naive row reduction.
Next, with \(j=1\) in the second equation, there is again no sum:
which is indeed the multipliers in the first step of naive row reduction.
Remember that one way to think of row reduction is recursively: after step \(k\), one just applies the same process recursively to the smaller \(n-k \times n-k\) matrix in the bottom-right-hand corner. We can do something similar here; at stage \(k\):
First use the first of the above equations to solve first for row \(k\) of \(U\), meaning just \(u_{k,j}, j \geq k\),
Then use the second equation to solve for column \(k\) of \(L\): \(l_{i,k}, i > k\).
Algorithm 2.7 (Doolittle factorization)
Stage \(k=1\) is handled by the simpler special equations above:
for j from 1 to n:
\(\quad u_{1,j}=a_{1,j}\)
end
for i from 1 to n:
\(\quad l_{i,1}=a_{i,1}/a_{1,1}\)
end
and for the rest:
for k from 2 to n:
\(\quad\) for j from k to n: \(\qquad\) Get the non-zero elements in row \(k\) of \(U\)
\(\quad\quad\) \(u_{k,j}=a_{k,j} - \sum_{s=1}^{k-1}l_{k,s}u_{s,j}\)
\(\quad\) end
\(\quad\) for i from k+1 to n: \(\qquad\) Get the non-zero elements in column \(k\) of \(L\) (except the 1’s on its diagonal)
\(\quad\quad\) \(l_{i,k}=\displaystyle\frac{a_{i,k}-\sum_{s=1}^{k-1}l_{i,s}u_{s,k}}{u_{k,k}}\)
\(\quad\) end
end
Remark 2.18
The special handling for \(k=1\) could be omitted, doing the main loop from 1 instead;
that would just have empty sums \sum_{s=1}^{0}.
That version is described in Algorithm 2.15 in section Operation Counts for Solving Simultaneous Linear Equations, where it is convenient for counting how much arithmetic is involved.
Note well that in the formulas to evaluate at the right,
The terms \(l_{k,s}\) are for \(s < k\), so from a column \(s\) that has already been computed for a previous \(k\) value.
The terms \(u_{s,j}\) are for \(s < k\), so from a row \(s\) that has already been computed for a previous \(k\) value.
The denominator \(u_{k,k}\) in the second inner loop is computed just in time, in the first inner loop for the same \(k\) value.
So the only thing that can go wrong is the same as before: a zero pivot element \(u_{k,k}\).
Remark 2.19 (On the Doolittle factorization algorithm)
For \(k=n\), the second inner loop is redundant, so could be eliminated. Indeed it might need to be eliminated in actual code, where “empty loops” might not be allowed. On the other hand, allowing empty loops makes the above correct also for \(k=1\); then the
for kloop encompases the entire factorization algorithm.This direct factorization algorithm avoids any intermediate modification of arrays, and thus eliminates all those superscripts like \(a_{i,j}^{(k)}\). This is not only nicer mathematically, but can help to avoid mistakes like code that inadvertently modifies the array containing the matrix \(A\) and then uses it to compute the residual, \(b - Ax\). More generally, such purely mathematical statements of algorithms can help to avoid coding errors; this is part of the philosophy of the functional programming approach.
Careful examination shows that the product \(l_{k,s} u_{s,j}\) that is part of what is subtracted at location \((k,j)\) is the same as what is subtracted there at stage \(k\) of row reduction, just with different names. More generally, every piece of arithmetic is the same as before, except arranged in a different order, so that the \(k-1\) changes made to an element in row \(k\) are done together, via those sums.
Rewrites with zero-based indexing provided in Linear algebra algorithms using 0-based indexing and semi-open intervals.
We can now verify the claim above that this is essentially naive row-reduction, just achieved in a different way and with new notation:
Theorem 2.6
If naive row reduction for \(Ax=b\) succeeds, giving row-reduced form \(Ux = c\):
The matrix \(A\) can be factorized as \(A = LU\) with \(U\) an \(n \times n\) upper triangular matrix and \(L\) an \(n \times n\) lower triangular matrix.
With the further condition that \(L\) is unit lower triangular there is a unique such factorization: This is called the Doolittle Factorization of \(A\).
In the Doolittle factorization, the matrix \(U\) is the one given by naive row reduction, and the elements of \(L\) below its main diagonal are the multipliers arising in naive row reduction. (The other elements of \(L\), on and above the main diagonal, are the ones and zeros dictated by it being unit lower triangular: the same as for those elements in the \(n \times n\) identity matrix.)
The transformed right-hand side \(c\) arising in naive row reduction is the solution of the system \(Lc = b\), and this is solvable by a procedure caled forward substitution, very similar to the backward subsitution used to solve \(Ux = c\).
Proof. First, note that in both naive row reduction (“NRR”) and the Doolittle factorization, the first row is unchanged, so \(u_{1,i}\) is the same in each case. Also, the multipliers in column one of NRR and in Doolittle are \(l_{i,1} = a_{i,1}/a_{1,1}\); see Equation (2.7). Thus at step \(k=1\), the first row of \(U\) and the first column of \(L\) are the same.
Next, we can work by induction, assuming that all is well up to stage \(k-1\), as we have checked for \(k=2\): the \(u_{i,j}\) are the same for rows \(i < k\) and the \(l_{i,j}\) are the same for columns \(j < k\).
At stage \(k\), NRR only changes rows \(i>k\), so the final form of the i-th row is that from stage \(i-1\): \(u_{i,j} = a_{i,j}^{(i-1)}\).
Thus the above formula (2.4) for elements of \(U\) can be rewritten as
and here, all the \(l_{i,k}\) are for columns \(k < i\), so are already know to be the same for both algorithms.
Look at the terms that are subtracted from \(a_{i,j}\): they are \(l_{i,k} a_{k,j}^{(k-1)}\) for \(1 \leq k < i\)
and these are the same terms subtracted in position \((i,j)\) in stages 1 through \(k-1\) of NRR: see Algorithm 2.1:
In other words, the arithmetic for getting \(u_{i,j}\) is the same, just with the Doolittle algorithm doing all the subtractions at once.
def lu_factorize(A, demo_mode=False):
"""Compute the Doolittle LU factorization of A.
Sums like $\\sum_{s=1}^{k-1} l_{k,s} u_{s,j}$ are done as matrix products;
in the above case, row matrix L[k, 1:k-1] by column matrix U[1:k-1,j] gives the sum for a give j,
and row matrix L[k, 1:k-1] by matrix U[1:k-1,k:n] gives the relevant row vector.
"""
n = len(A) # len() gives the number of rows in a 2D array.
# Initialize U as the zero matrix;
# correct below the main diagonal, with the other entries to be computed below.
U = zeros_like(A)
# Initialize L as the identity matrix;
# correct on and above the main diagonal, with the other entries to be computed below.
L = identity(n)
# Column and row 1 (i.e Python index 0) are special:
U[0,:] = A[0,:]
L[1:,0] = A[1:,0]/U[0,0]
if demo_mode:
print(f"After step k=0")
print(f"U=\n{U}")
print(f"L=\n{L}")
for k in range(1, n-1):
U[k,k:] = A[k,k:] - L[k,:k] @ U[:k,k:]
L[k+1:,k] = (A[k+1:,k] - L[k+1:,:k] @ U[:k,k])/U[k,k]
if demo_mode:
print(f"After step {k=}")
print(f"U=\n{U}")
print(f"L=\n{L}")
# The last row (index "-1") is special: nothing to do for L
U[-1,-1] = A[-1,-1] - sum(L[-1,:-1]*U[:-1,-1])
if demo_mode:
print(f"After the final step, k={n-1}")
print(f"U=\n{U}")
return (L, U)
2.4.2.1. A test case on LU factorization#
A = array([[4, 2, 7], [3, 5, -6],[1, -3, 2]], dtype=float)
print(f"A=\n{A}")
A=
[[ 4. 2. 7.]
[ 3. 5. -6.]
[ 1. -3. 2.]]
(L, U) = lu_factorize(A, demo_mode=True)
After step k=0
U=
[[4. 2. 7.]
[0. 0. 0.]
[0. 0. 0.]]
L=
[[1. 0. 0. ]
[0.75 1. 0. ]
[0.25 0. 1. ]]
After step k=1
U=
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. 0. ]]
L=
[[ 1. 0. 0. ]
[ 0.75 1. 0. ]
[ 0.25 -1. 1. ]]
After the final step, k=2
U=
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. -11. ]]
print(f"A=\n{A}")
print(f"L=\n{L}")
print(f"U=\n{U}")
print(f"L times U is \n{L@U}")
print(f"The 'residual' A - LU is \n{A - L@U}")
A=
[[ 4. 2. 7.]
[ 3. 5. -6.]
[ 1. -3. 2.]]
L=
[[ 1. 0. 0. ]
[ 0.75 1. 0. ]
[ 0.25 -1. 1. ]]
U=
[[ 4. 2. 7. ]
[ 0. 3.5 -11.25]
[ 0. 0. -11. ]]
L times U is
[[ 4. 2. 7.]
[ 3. 5. -6.]
[ 1. -3. 2.]]
The 'residual' A - LU is
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
2.4.2.2. Forward substitution: solving \(Lc = b\) for \(c\)#
This is the last piece missing. The strategy is very similar to backward substitution, but slightly simplified by the ones on the main didogonal of \(L\). The equations \(L c = b\) can be written much as above, separating off the last term in the sum:
Then solve for \(c_i\):
These are already is usable order: the right-hand side in the equation for \(c_i\) involves only the \(c_j\) values with \(j < i\), determined by earlier equations if we run through index \(i\) in increasing order.
First, \(i=1\)
Next, \(i=2\)
Next, \(i=3\)
Algorithm 2.8 (forward substitution with a unit lower triangular matrix)
\(c_1 = b_1\)
for i from 2 to n:
\(\displaystyle \quad c_i = b_i - \sum_{j=1}^{i-1} l_{i,j} c_j\)
end
Remark 2.20
As with Algorithm 2.7 above, the special handling for \(k=1\) could be omitted, doing the main loop from 1 instead; that version is described in Algorithm 2.14 in section Operation Counts for Solving Simultaneous Linear Equations.
I leave expressing this in zero-based pseudo-code to Exercise 2.9;
there is Python code for it in the module numerical_methods, described at
forward_substitution
in the appendix
numerical_methods,
so we can import it with
from numerical_methods import forward_substitution
(Note however that this version is not in the form asked for in the following exercise!)
See Exercise 2.8
2.4.2.3. A test case on forward substitution#
b = array([2., 3., 4.])
c = forward_substitution(L, b)
print(f"c = {c}")
print(f"The residual b - Lc is {b - L@c}")
print(f"\t with maximum norm {max(abs(b - L@c)):0.3}")
c = [2. 1.5 5. ]
The residual b - Lc is [0. 0. 0.]
with maximum norm 0.0
2.4.2.4. Completing the test case, with backward substitution#
As this step is unchanged, just import the version seen in Row Reduction (Gaussian Elimination).
from numerical_methods import backward_substitution
x = backward_substitution(U, c)
print(f"The residual c - Ux for the backward substitution step is {c - U@x}")
print(f"\t with maximum norm {max(abs(c - U@x)):0.3}")
print(f"The residual b - Ax for the whole solving process is {b - A@x}")
print(f"\t with maximum norm {max(abs(b - A@x)):0.3}")
The residual c - Ux for the backward substitution step is [ 0.00000000e+00 -2.22044605e-16 0.00000000e+00]
with maximum norm 2.22e-16
The residual b - Ax for the whole solving process is [0.0000000e+00 0.0000000e+00 8.8817842e-16]
with maximum norm 8.88e-16
2.4.2.4.1. Creating modules#
See the notes on creating Python modules in the Python Tutorial; for now I repeat those notes here.
A module is created simply as a Python file containing a collection of defibitions of functions, variables and such;
for example, the module numerical_methods used here is defined by the file numerical_methods.py
One way to create such a file is with Spyder (or another Python IDE, or indeed any text editor).
However, there is the alternative of working with JupyterLab and Jupyter notebooks: first put the definitions (and optionally, some comments and testing code) in a notebook, and then convert that to an eponymous Python code file with the JupyerLab menu command
File > Export Notebook As ... > Export Notebook to Executable Script
This is how the module numerical_methods is created from the notebook file linear_algebra.ipynb in the Appendices.
One advantage of this aproach is that the notebook can contain useful documentation, including mathematical descriptions of algorithms.
2.4.3. When does LU factorization work?#
Theorem 2.1 says that naive row reduction works with SDD matrices (in the sense of avoiding division by zero); the corresponding result here is:
Theorem 2.7
(A) Any matrix that is strictly diagonally dominant be either rows or columns
with the factor \(U\) also strictly diagonally dominant diagonal in the same direction (i.e. by rows of columns). In particular, the diagonal elements of \(U\) all non-zero, so backward substitution also works.
(B) Further, when the matrix is column-wise SDD matrix, this LU factorization is also “optimal”, in the sense that maximal element partial pivoting would in fact not make any row swaps.
Proof. For part (A), see [Dahlquist and Björck, 2009], Theorem 7.2.6.
Part (B) follows from the fact that, as seen in Theorem 2.2, column-wise SDD matrices have all \(|l_{i,j}| \le 1\), which is exactly what the row swaps in maximal element partial pivoting are chosen to achieve, so this pivoting is un-needed with such matrices.
This nice second property can be got for row-wise SDD matrices via a twist, or actually a transpose.
For an SDD matrix, it transpose \(B = A^T\) is column-wise SDD and so has the nice Doolitle factorization described above: \(B = L_B U_B\), with \(L_B\) being column-wise diagonally dominant and having ones on the main diagonal.
Transposing back, \(A = B^T = (L_B U_B)^T = U_B^T L_B^T\), and defining \(L = U_B^T\) and \(U = L_B^T\),
\(L\) is lower triangular
\(U\) is upper triangular, row-wise diagonally dominant and with ones on its main diagonal: it is “unit upper triangular”.
This is thus another LU factorization of \(A\), with \(U\) rather than \(L\) being the factor with ones on its main diagonal.
2.4.4. Crout factorization#
This sort of \(LU\) factorization is called the Crout factorization; as with the Doolittle version, if such a factorization exists, it is unique.
As a consequence of Theorem Theorem 2.7 above:
Theorem 2.8
Every matrix that is SDD by either rows or columns has a Crout factorization, and all the diagonal elements are non-zero, so forward- and backward-substitution work.
With the Crout factorization, the diagonal elements of \(L\) are not all ones, so a more general version of forward-substitution is needed, with divisions by the diagonal elements as in Algorithm 2.5 for backward substitution:
Algorithm 2.9 (forward substitution with a general lower triangular matrix)
\(c_1 = b_1 / l_{1,1}\)
for i from 2 to n:
\(\displaystyle \quad c_i = \frac{b_i - \sum_{j=1}^{i-1} l_{i,j} c_j}{l_{i,i}}\)
end
2.4.5. Cholesky factorization#
As was mentioned at the end of the section Row Reduction (Gaussian Elimination), naive row reduction also works for positive definite matrices, and thus so do both of these LU factorizations.
However, there is another LU factorization that works even better in that case: the Cholesky factorization, where \(U = L^T\) so that \(A = L L^T\).
The derivation is omitted (see [Dionne, 2023] Section 4.3 , Cholesky Factorization), but here is the algorithm. (Note that forward-substitution then requires the more general form of Algorithm 2.9.)
Algorithm 2.10 (Cholesky factorization)
\(l_{1,1} = \sqrt{a_{1,1}}\)
for i from 2 to n:
\(\quad l_{i,1} = a_{i,1}/a_{1,1}\)
end
for k from 2 to n-1:
\(\quad\) \(l_{k,k} = \sqrt{{a_{k,k}} - \sum_{j=1}^{k-1}l_{k,j}^2}\)
\(\quad\) for i from k+1 to n:
\(\quad\quad\) \(l_{i,k} = \displaystyle\frac{a_{i,k} - \sum_{s=1}^{k-1}l_{i,s} l_{k,s}}{l_{k,k}}\)
\(\quad\) end
\(l_{n,n} = \sqrt{{a_{n,n}} - \sum_{j=1}^{n-1}l_{n,j}^2}\)
end
2.4.6. Computing Determinants#
Theorem 2.9
For an \(n \times n\) matrix with Dolittle factorization \(A = LU\), its determinant is the product of the main diagonal elements of \(U\):
More generally, for any LU factorization (e.g. Crout’s), it is the product of all the diagonal elements of both matrices:
Proof. This follows from the fact that for any triangular matrix, the determinant is the product of its main diagonal elements (plus the fact that these are all 1’s for the unit-triangular matrix \(L\) in the Dolittle factorization).
To check that, recall the formula for the determinant in terms of sums and differences of all the \(n!\) possible products consisting of \(n\) elements, one from each row and column: for a triangular matrix, this is just the “main diagonal product” because each of the others includes elements from both above and below that main diagonal, and so had a factor of zero.
2.4.7. Computing the Inverse of a Matrix (for emergency use only!)#
It is usually neither necessary nor desirable to explicitly compute the inverse of a matrix; for example, the process of computing \(A^{-1}\) and then solving \(Ax = b\) with \(x = A^{-1} x\) requires about four times as much computational effort as either solving with row reduction or by the “LU” method seen in this section.
However, the inverse is needed in some situations, and then the LU factorization gives a relatively efficient way of finding it.
To describe the strategy, we use the superscript notation \(M^{(i)}\) for column \(i\) of a matrix \(M\). Thus for the identity matrix, \(I^{(i)}\) is the column vector whose \(i\)-th entry is 1 and all ohers are zero, better knwn as the unit basis vector \(e^{(i)}\).
The strategy is then:
Note that \(A A^{-1} = I\)
For any matrix B, column \(i\) of the product \(A B\) is the product of \(A\) with column \(i\) of \(B\); using the above superscript notation for columns, \(A (B^{(i)}) = (A B)^{(i)}\).
In particular, \(A\) times column \(i\) of its inverse is \(I^{(i)}\) a.k.a. \(e^{(i)}\): \(A (A^{-1})^{(i)} = e^{(i)}\).
Each column of the inverse can be computed by solving the above; that is, solve \(Ax = b\) successively for each choice \(b = e^{(i)}\).
This only requires computing an LU factorization of \(A\) once, and then doing only the far faster forward- and backward- substitution steps for each column of the inverse.
2.4.8. Exercises#
Exercise 2.8
Solve the equation \(A \vec{x} = \vec b\) with
(as seen in Exercise 2.1 of of Row Reduction (Gaussian Elimination) and Exercise 2.7 of Partial Pivoting), this time by computing and using the (Doolittle) \(LU\) factorization of \(A\).
As in that earlier exercise, compute all intermediate values as decimal approximations rounded to four significant digits. Note that one can evaluate the sum giving each new value \(u_{i,j}\) and \(l_i,j\) “exactly” and then round the result; this emulates the way that a computer can evaluate compound expression “exactly” within the processor and then round the final result for storing back to memory.
Call this approximate solution \(\vec{x}_3\), and compute the residual \(\vec{r}_3 = \vec{b} - A \vec{x}_3\) and its maximum norm \(\| \vec{r}_3 \|_\infty\). This should be done with high accuracy (not rounding to four decimal places), and could be done using a Python notebook as a “calculator”.
Compare these results to those of Exercise 2.1 and Exercise 2.7, and discuss.
See Exercise 2.16 for a followup on this.
Exercise 2.9
A) Express the forward substitution strategy of Forward substitution: solving Lc = b for c as pseudo-code, adjusting to Python’s zero-based indexing. Spell out all the sums explicitly rather than using “\(\Sigma\)” notation for sums or any other matrix multiplication short-cut.
B) Then implement it “directly” in a Python function, with format:
function forward_substitution(L, b)
. . .
return c
Again do this with explicit evaluation of each sum rather than using the function sum or any matrix multiplication short-cut.
C) Test it, using this often-useful “reverse-engineering” tactic:
Create suitable test arrays
Landc. (Use \(n\) at least three, and preferably larger.)Compute their product, with
b = L @ c(If this notation is unfamiliar, see Numpy Array Operations and Linear Algebra in the Python Tutorial.)Check if
c_solution = forward_substitution(L, b)gives the correct value (within rounding error.)
The following exercises all refer to the equations
The explorations here wil be continued in the exercises of section Error bounds for linear algebra, matrix norms, condition numbers, etc. where we learn more about measuring and bounding errors.
Exercise 2.10
Express the system of equations (2.8) in matrix-vector form \(Ax= b\) and compute the LU factorization of \(A\), rounding to three significant digits: call these factors \(L_a\) and \(U_a\). (Show working; this is for comparison to “exact” (meaning far more accurate) values to be calculated below.)
Use these factors to solve \(L_a c = b\) for \(c\), again rounding at each step to three significant digits: call this approximation \(c_a\).
Complete this approximate solution process by solving \(U_a x_a = c_a\) for \(x\), again rounding at each arithmetic step to three significant digits.
Compute the residual \(r = b - A x_a\) and its maximum norm \(\|r\|_\infty\). This time compute “exactly”, or at least to far more that three significant digits; for example to do this using Python as your calculator.
Aside: you could also check each stage of the working with the intermediate residuals or backward errors: \(A - L_a U_a\), \(b - L_a c_a\) and \(c_a - U_a x_a\). These backwrd error checks can useful for debugging code, as in the next exercise.
Exercise 2.11
Write Python implementations of the algorithms used above (the Doolittle factorization and forward and backward substitution) in a file linalg.py providing module linalg, and use this to compute “exact” results for \(L\), \(U\), \(c\) and \(x\).
Use the results of Exercise 2.10 to compute the absolute and relative errors in \(x_a\).
Exercise 2.12
Use the Python code from Exercise 2.11 and the LU factorization that it gives to compute the inverse \(A^{-1}\).
Check your work in the now familar way, by evalauting the “residual” \(I - A A^{-1}\).
Remark 2.21
There is a continutation of this exploration in Exercise 2.15 of Error bounds for linear algebra, matrix norms, condition numbers, etc..
Exercise 2.13
(An ongoing activity.)
Start building a Python module linear_algebra in a file linear_algebra.py, with all our linear algebra functions: for now, forward_substitution(L, b) as above and also row_eeduce(A, b) and backward_substitution(U, c) from the section
Row Reduction (Gaussian Elimination)
Include testing/demonstration at the bottom of the module definition file, in the block of the if statement
if __name__ == "__main__":
We will add the Doolittle method and such in a while, and could use this module in assignments and projects.