Lecture 10: Effects of finite precision arithmetic#
The need for row swapping in GE#
Consider the following linear system of equations
Problem. We cannot eliminate the first column by the diagonal by adding multiples of row 1 to rows 2 and 3 respectively.
Row swapping#
Solution. Swap the order of the equations!
Swap rows 1 and 2:
\[\begin{split} \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 1 & 2 & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 7 \\ 5 \end{pmatrix} \end{split}\]Now apply Gaussian elimination
\[\begin{split} \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 1.5 & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 7 \\ 3 \end{pmatrix} ; \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & -0.75 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 7 \\ -2.25 \end{pmatrix}. \end{split}\]
Another example#
Consider another system of equations
Apply Gaussian elimination as usual:
\[\begin{split} \begin{pmatrix} 2 & 1 & 1 \\ 0 & 0 & -1 \\ 2 & 2 & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 3 \\ -1 \\ 2 \end{pmatrix} ; \begin{pmatrix} 2 & 1 & 1 \\ 0 & 0 & -1 \\ 0 & 1 & -1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 3 \\ -1 \\ -1 \end{pmatrix} \end{split}\]Problem. We cannot eliminate the second column below the diagonal by adding a multiple of row 2 to row3.
Again this problem may be overcome simply by swapping the order of the equations - this time swapping rows 2 and 3:
\[\begin{split} \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & -1 \\ 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 3 \\ -1 \\ -1 \end{pmatrix} \end{split}\]We can now continue the Gaussian elimination process as usual.
In general. Gaussian elimination requires row swaps to avoid breaking down when there is a zero in the “pivot” position.
Problems with finite precision#
Consider using Gaussian elimination to solve the linear system of equations given by
where \(\varepsilon \neq 1\).
The true, unique solution is \((x_1, x_2)^T = (1, 2)^T\).
If \(\varepsilon \neq 0\), Gaussian elimination gives
\[\begin{split} \begin{pmatrix} \varepsilon & 1 \\ 0 & 1 - \frac{1}{\varepsilon} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 2 + \varepsilon \\ 3 - \frac{2 + \varepsilon}{\varepsilon} \end{pmatrix} \end{split}\]Problems occur not only when \(\varepsilon = 0\) but also when it is very small, i.e. when \(\frac{1}{\varepsilon}\) is very large, this will introduce very significant rounding errors into the computation.
Removing the problem#
Use Gaussian elimination to solve the linear system of equations given by
where \(\varepsilon \neq 1\).
The true solution is still \((x_1, x_2)^T = (1, 2)^T\).
Gaussian elimination now gives
\[\begin{split} \begin{pmatrix} 1 & 1 \\ 0 & 1 - \varepsilon \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 3 \\ 2 - 2\varepsilon \end{pmatrix} \end{split}\]The problems due to small values of \(\varepsilon\) have disappeared.
Notes#
Writing the equations in a different order has removed the previous problem.
The diagonal entries are now always relatively larger.
The interchange of the order of equations is a simple example of row pivoting. This strategy avoids excessive rounding errors in the computations.
Gaussian elimination with pivoting#
Before eliminating entries in column \(j\):
find the entry in column \(j\), below the diagonal, of maximum magnitude;
if that entry is larger in magnitude than the diagonal entry then swap its row with row \(j\).
Then eliminate column \(j\) as before.
Notes#
This algorithm will always work when the matrix \(A\) is invertible/non-singular.
Conversely, if all of the possible pivot values are zero this implies that the matrix is singular and a unique solution does not exist.
At each elimination step the row multiplies used are guaranteed to be at most one in magnitude…
… so any errors in the representation of the system cannot be amplified by the elimination process.
As always, solving \(A \vec{x} = \vec{b}\) requires that the entries in \(\vec{b}\) are also swapped in the appropriate way.
Pivoting can be applied in an equivalent way to LU factorisation.
The sequence of pivots is independent of the vector \(\vec{b}\) and can be recorded and reused.
The constraint imposed on the row multipliers means that for LU factorisation every entry in \(L\) satisfies \(| l_{ij} | \le 1\).
In python, the function call
P, L, U = scipy.linalg.lu(A, permute_l=0)
factorises \(A\) and returns \(L\), \(U\) and the pivot matrix \(P\).
Example#
Consider the linear system of equations given by
where \(0 \le \varepsilon \ll 1\), and solve it using
Gaussian elimination without pivoting
Gaussian elimination with pivoting.
The exact solution is \(\vec{x} = (0, -1, 2)^T\) for any \(\varepsilon\) in the given range.
1. Solve the system using Gaussian elimination with no pivoting.#
Eliminating the first column gives
and then the second column gives
which leads to
There are many divisions by \(\varepsilon\), so we will have problems if \(\varepsilon\) is small.
2. Solve the system using Gaussian elimination with pivoting.#
The first stage is identical (because \(a_{11} = 10\) is largest).
but now \(|a_{22}| = \varepsilon\) and \(|a_{32}| = 2.5\) so we swap rows 2 and 3 to give
Now we may eliminate column 2:
which leads to the exact answer:
Further reading#
Mathematics stackexchange: what are pivot numbers in LU decomposition? please explain me in an example
Trefethen, Lloyd N.; Bau, David (1997), Numerical linear algebra, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.