5.1: Gaussian Elimination Revisited
\(\newcommand{\twovec}[2]{\begin{pmatrix} #1 \\ #2 \end{pmatrix} } \)
\(\newcommand{\threevec}[3]{\begin{pmatrix} #1 \\ #2 \\ #3 \end{pmatrix} } \)
\(\newcommand{\fourvec}[4]{\begin{pmatrix} #1 \\ #2 \\ #3 \\ #4 \end{pmatrix} } \)
\(\newcommand{\fivevec}[5]{\begin{pmatrix} #1 \\ #2 \\ #3 \\ #4 \\ #5 \end{pmatrix} } \)
In this section, we revisit Gaussian elimination and explore some problems with implementing it in the straightforward way that we described back in Section 1.2 . In particular, we will see how the fact that computers only approximate arithmetic operations can lead us to find solutions that are far from the actual solutions. Second, we will explore how much work is required to implement Gaussian elimination and devise a more efficient means of implementing it when we want to solve equations \(A\mathbf x = \mathbf b\) for several different vectors \(\mathbf b\text{.}\)
Preview Activity 5.1.1.
To begin, let's recall how we implemented Gaussian elimination by considering the matrix
- What is the first row operation we perform? If the resulting matrix is \(A_1\text{,}\) find a matrix \(E_1\) such that \(E_1A = A_1\text{.}\)
- What is the matrix inverse \(E_1^{-1}\text{?}\) You can find this using your favorite technique for finding a matrix inverse. However, it may be easier to think about the effect that the row operation has and how it can be undone.
-
Perform the next two steps in the Gaussian elimination algorithm to obtain \(A_3\text{.}\) Represent these steps using multiplication by matrices \(E_2\) and \(E_3\) so that
\begin{equation*} E_3E_2E_1A = A_3\text{.} \end{equation*}
- Suppose we need to scale the second row by \(-2\text{.}\) What is the \(3\times3\) matrix that perfoms this row operation by left multiplication?
- Suppose that we need to interchange the first and second rows. What is the \(3\times3\) matrix that performs this row operation by left multiplication?
Partial pivoting
The first issue that we will address is the fact that computers do not perform arithemtic operations exactly. For instance, if we ask Python to evaluate
0.1 + 0.2
, it reports
0.30000000000000004
though we know that the true value is 0.3. There are a couple of reasons for this.
First, computers perform arithmetic using base 2 numbers, which means that numbers we enter in decimal form, such as \(0.1\text{,}\) must be converted to base 2. Even though 0.1 has a simple decimal form, its representation in base 2 is the repeating decimal
To accurately represent this number inside a computer would require an infinite number of digits. Since a computer can only hold a finite number of digits, we are necessarily using an approximation when simply representing this number in a computer.
In addition, arithmetic operations, such as addition, are prone to error. To keep things simple, suppose we have a computer that represents numbers using only three decimal digits. For instance, the number 1.023 would be represented as
1.02
while 0.023421 would be
0.0234
. If we add these numbers, we have 1.023 + 0.023421 = 1.046421; the computer reports this sum as
1.02 + 0.0234 = 1.04
, whose last digit is not correctly rounded. Generally speaking, we will see this problem, which is called
round off error
, whenever we add numbers of signficantly different magnitudes.
Remember that Gaussian elimination, when applied to an \(n\times n\) matrix, requires approximately \(\frac 23 n^3\) operations. If we have a \(1000\times1000\) matrix, performing Gaussian elimination requires roughly a billion operations, and the errors introduced in each operation could accumulate. How can we have confidence in the final result? We can never completely avoid these errors, but we can take steps to mitigate them. The next activity will introduce one such technique.
Activity 5.1.2.
Suppose we have a hypothetical computer that represents numbers using only three decimal digits, as we have been discussing. We will consider the linear system
-
Show that this system has the unique solution
\begin{equation*} \begin{aligned} x & {}={} \frac{10000}{9999} = 1.00010001\ldots, \\ y & {}={} \frac{9998}{9999} = 0.99989998\ldots\text{.} \end{aligned} \end{equation*}
- If we represent this solution inside our computer that only holds 3 decimal digits, what do we find for the solution? This is the best that we can hope to find using our computer.
- Let's imagine that we use our computer to find the solution using Gaussian elimination; that is, after every arithmetic operation, we keep only three decimal digits. Our first step is to multiply the first equation by 10000 and subtract it from the second equation. If we represent numbers using only three decimal digits, what does this give for the value of \(y\text{?}\)
- By substituting our value for \(y\) into the first equation, what do we find for \(x\text{?}\)
- Compare the solution we find on our computer with the actual solution and assess the quality of the approximation.
-
Let's now modify the linear system by simplying interchanging the equations:
\begin{equation*} \begin{alignedat}{3} x & {}+{} & y & {}={} & 2 \\ 0.0001x & {}+{} & y & {}={} & 1\text{.} \\ \end{alignedat} \end{equation*}
Of course, this doesn't change the actual solution. Let's imagine we use our computer to find the solution using Gaussian elimination. Perform the first step where we multiply the first equation by 0.0001 and subtract from the second equation. What does this give for \(y\) if we represent numbers using only three decimal digits?
- Substitute the value you found for \(y\) into the first equation and solve for \(x\text{.}\) Then compare the approximate solution found with our hypothetical computer to the exact solution.
- Which approach produces the most accurate approximation?
This activity demonstrates how the practical aspects of computing differ from the theoretical. We know that the order in which we write the equations has no effect on the solution space; row interchange is one of our three allowed row operations in the Gaussian elimination algorithm. However, when we are only able to perform arithmetic operations approximately, applying row interchanges can dramatically improve the accuracy of our approximations.
If we could compute the solution exactly, we find
Since our hypothetical computer represents numbers using only three decimal digits, our computer finds
This is the best we can hope to do with our computer since it is impossible to represent the solution exactly.
When the equations are written in their original order and we multiply the first equation by 10000 and subtract from the second, we find
In fact, we find the same value for \(y\) when we interchange the equations. Here we multiply the first equation by 0.0001 and subtract from the second equation. We then find
The difference occurs when we substitute \(y\approx 1\) into the first equation. When the equations are written in their original order, we have
When the equations are written in their original order, we find the solution \(x\approx 0.00, y \approx 1.00\text{.}\)
When we write the equation in the opposite order, however, substituting \(y\approx 1\) into the first equation gives
In this case, we find the approximate solution \(x\approx 1.00, y\approx1.00\text{,}\) which is the most accurate solution that our hypothetical computer can find. Simply interchanging the order of the equation produces a much more accurate solution.
We can understand why this works graphically. Each equation represents a line in the plane, and the solution is the intersection point. Notice that the slopes of these lines differ considerably.
When the equations are written in their original order, we substitute \(y\approx1\) into the equation \(0.00001x + y = 1\text{,}\) which is a nearly horizontal line. Along this line, a small change in \(y\) leads to a large change in \(x\text{.}\) The slight difference in our approximation \(y\approx 1\) from the exact value \(y=0.9998999\ldots\) leads to a large difference in the approximation \(x\approx0\) from the exact value \(x=1.00010001\ldots\text{.}\)
If we exchange the order in which the equations are written, we substitute our approximation \(y\approx 1\) into the equation \(x+y=2\text{.}\) Notice that the slope of the associated line is \(-1\text{.}\) On this line, a small change in \(y\) leads to a relatively small change in \(x\) as well. Therefore, the difference in our approximation \(y\approx1\) from the exact value leads to only a small difference in the approximation \(x\approx1\) from the exact value.
This example motivates the technique that computers usually use to perform Gaussian elimation. In practice, we only need to perform a row interchange when a zero occurs in a pivot position, such as
Instead, we will perform a row interchange to put the entry having the largest possible absolute value entry into the pivot position. For instance, when performing Gaussian elimination on the following matrix, we begin by interchanging the first and third rows so that the upper left entry has the largest possible absolute value.
This technique is called partial pivoting , and it means that, in practice, we will perform many more row interchange operations than we typically do when computing exactly by hand.
\(LU\) factorizations
In Subsection 1.3 , we saw that the number of arithmetic operations needed to perform Gaussian elimination on an \(n\times n\) matrix is about \(\frac23 n^3\text{.}\) This means that a \(1000\times1000\) matrix, requires about two thirds of a billion operations.
Suppose that we have two equations, \(A\mathbf x = \mathbf b_1\) and \(A\mathbf x = \mathbf b_2\text{,}\) that we would like to solve. Usually, we would form augmented matrices \(\left[\begin{array}{c|c} A & \mathbf b_1 \\ \end{array}\right]\) and \(\left[\begin{array}{c|c} A & \mathbf b_2 \\ \end{array}\right]\) and apply Gaussian elimination. Of course, the steps we perform in these two computations are nearly identical. Is there a way to store some of the computation we perform in reducing \(\left[\begin{array}{c|c} A & \mathbf b_1 \\ \end{array}\right]\) and reuse it in solving subsequent equations? The next activity will point us in the right direction.
Activity 5.1.3.
We will consider the matrix
and begin performing Gaussian elimination.
-
Perform two row replacement operations to find the row equivalent matrix
\begin{equation*} A' = \left[\begin{array}{rrr} 1 & 2 & 1 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \\ \end{array}\right]\text{.} \end{equation*}
Find elementary matrices \(E_1\) and \(E_2\) that perform these two operations so that \(E_2E_1 A = A'\text{.}\)
-
Perform a third row replacement to find the upper triangular matrix
\begin{equation*} U = \left[\begin{array}{rrr} 1 & 2 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array}\right]\text{.} \end{equation*}
Find the elementary matrix \(E_3\) such that \(E_3E_2E_1A = U\text{.}\)
-
We can write \(A=E_1^{-1}E_2^{-1}E_3^{-1} U\text{.}\) Find the inverse matrices \(E_1^{-1}\text{,}\) \(E_2^{-1}\text{,}\) and \(E_3^{-1}\) and the product \(L = E_1^{-1}E_2^{-1}E_3^{-1}\text{.}\) Then verify that \(A=LU\text{.}\)
-
Suppose that we want to solve the equation \(A\mathbf x = \mathbf b = \threevec4{-7}{12}\text{.}\) We will write
\begin{equation*} A\mathbf x = LU\mathbf x = L(U\mathbf x) = \mathbf b \end{equation*}
and introduce an unknown vector \(\mathbf c\) such that \(U\mathbf x = \mathbf c\text{.}\) Find \(\mathbf c\) by noting that \(L\mathbf c = \mathbf b\) and solving this equation.
- Now that we have found \(\mathbf c\text{,}\) find \(\mathbf x\) by solving \(U\mathbf x = \mathbf c\text{.}\)
- Use the factorization \(A=LU\) and this two-step process, solve the equation \(A\mathbf x = \threevec{2}{-2}{7}\text{.}\)
This activity introduces a method for factoring a matrix \(A\) as a product of two triangular matrices, \(A=LU\text{,}\) where \(L\) is lower triangular and \(U\) is upper triangular. The key to finding this factorization is to represent the row operations that we apply in the Gaussian elimination algorithm through multiplication by elementary matrices.
Beginning with the matrix \(A = \left[\begin{array}{rrr} 1 & 2 & 1 \\ -2 & -3 & -2 \\ 3 & 7 & 4 \\ \end{array}\right] \text{,}\) we multiply the first row by \(2\) and add the result to the second row. The result is obtained by multiplying \(A\) by the matrix \(E_1 = \left[\begin{array}{rrr} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array}\right] \text{.}\) Notice that the inverse of \(E_1\) is \(E_1^{-1} = \left[\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 0 \\ \end{array}\right] \text{.}\) There is no computation required to find \(E_1^{-1}\text{;}\) to undo the effect of \(E_1\text{,}\) we simply multiply the first row by \(-2\text{,}\) rather than \(2\text{,}\) and add it to the second row. In practice, to find the inverse of an elementary matrix corresponding to a row replacement operation, we simply change the sign of the number below the diagonal. That's easy!
Likewise, the next two row replacements are performed by
From here, we find \(A = E_1^{-1}E_2^{-1}E_3^{-1} U\) so that
If we want to solve the equation \(A\mathbf x = LU\mathbf x = \mathbf b\text{,}\) we first solve for \(U\mathbf x = \mathbf c\) by solving \(L\mathbf c = \mathbf b\text{.}\) Once we have \(\mathbf c\text{,}\) we solve \(U\mathbf x = \mathbf c\) for the vector \(\mathbf x\text{,}\) which is the solution to the original equation \(A\mathbf x = \mathbf b\text{.}\) In our example,
If we want to solve \(A\mathbf x = \mathbf b\) for a different right-hand side \(\mathbf b\text{,}\) we can simply repeat this two-step process.
We have now introduced \(LU\) factorizations as a technique for solving equations \(A\mathbf x = \mathbf b\text{,}\) but we have yet to explain why they are useful. It becomes apparent when we consider the type of equations that we are solving; for instance, \(L\mathbf c = \mathbf b\) in our example has the form
Because \(L\) is a lower-triangular matrix, it contains lots of zeros. In fact, we are able to read off the first component of \(\mathbf c\) directly from the equations: \(c_1 = 4\text{.}\) We then have \(-2c_1+c_2 = -7\text{,}\) which gives \(c_2 = 1\text{,}\) and \(3c_1 + c_2 + c_3 = 12\text{,}\) which gives \(c_3=-1\text{.}\) As this shows, solving equations involving triangular matrices is quite easy because we essentially only need to perform a sequence of substitutions.
In fact, solving an equation with an \(n\times n\) triangular matrix requires approximately \(\frac 12 n^2\) operations. Once we have the factorization \(A=LU\text{,}\) we solve the equation \(A\mathbf x=\mathbf b\) by solving two equations involving triangular matrices, which requires about \(n^2\) operations. For example, if \(A\) is a \(1000\times1000\) matrix, we solve the equation \(A\mathbf x = \mathbf b\) using about one million steps. The compares with the third of a billion operations we need to perform Gaussian elimination. That represents a significant savings. Of course, we have to first find the \(LU\) factorization of \(A\) and this requires roughly the same amount of work as performing Gaussian elimination. However, once we have the \(LU\) factorization, we can use it to solve \(A\mathbf x=\mathbf b\) for different right hand sides \(\mathbf b\text{.}\)
Our discussion so far has ignored one issue, however. Remember that we sometimes have to perform row interchange operations in addition to row replacement. A typical row interchange is represented by multiplication by a matrix such as
which has the effect of interchanging the first and third rows. Notice that this matrix is not triangular so performing a row interchange will ruin the \(LU\) factorization we seek. We will not give the details here, but linear algebra software packages provide a matrix \(P\) that describes how the rows are permuted in the Gaussian elimination process. In particular, we will write \(PA = LU\text{,}\) where \(P\) is a permutation matrix, \(L\) is lower triangular, and \(U\) is upper triangular.
Therefore, to solve the equation \(A\mathbf x = \mathbf b\text{,}\) we first multiply both sides by \(P\) to obtain
That is, we multiply \(\mathbf b\) by \(P\) and then find \(\mathbf x\) using the factorization: \(L\mathbf c = P\mathbf b\) and \(U\mathbf x = \mathbf c\text{.}\)
Activity 5.1.4.
Sage will create \(LU\) factorizations; once we have a matrix
A
, we write
P, L, U = A.LU()
to obtain the matrices \(P\text{,}\) \(L\text{,}\) and \(U\) such that \(PA = LU\text{.}\)
-
In the previous activity, we found the \(LU\) factorization
\begin{equation*} A = \left[\begin{array}{rrr} 1 & 2 & 1 \\ -2 & -3 & -2 \\ 3 & 7 & 4 \\ \end{array}\right] = LU = \left[\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 3 & 1 & 1 \end{array}\right] \left[\begin{array}{rrr} 1 & 2 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array}\right]\text{.} \end{equation*}
Using Sage, define the matrix \(A\) and then ask Sage for the \(LU\) factorization. What are the matrices \(P\text{,}\) \(L\text{,}\) and \(U\text{?}\)
Notice that Sage finds a different \(LU\) factorization than we found in the previous activity. This is because Sage uses partial pivoting, as described in the previous section, when it performs Gaussian elimination. This is reflected by the fact that the permutation \(P\) is not the identity.
- Define the vector \(\mathbf b = \threevec4{-7}{12}\) in Sage and compute \(P\mathbf b\text{.}\)
-
Use the matrices
L
andU
to solve \(L\mathbf c = P\mathbf b\) and \(U\mathbf x = \mathbf c\text{.}\) You should find the same solution \(\mathbf x\) that you found in the previous activity. - Use the factorization to solve the equation \(A\mathbf x = \threevec{11}{-17}{35}\text{.}\)
- How does the factorization show us that \(A\) is invertible and that, therefore, every equation \(A\mathbf x=\mathbf b\) has a unique solution?
-
Suppose that we have the matrix
\begin{equation*} A = \left[\begin{array}{rrr} 3 & -1 & 2 \\ 2 & -1 & 1 \\ 2 & 1 & 3 \\ \end{array}\right]\text{.} \end{equation*}
Use Sage to find the \(LU\) factorization. Explain how the factorization shows that \(A\) is not invertible.
-
Consider the matrix
\begin{equation*} A = \left[\begin{array}{rrrr} -2 & 1 & 2 & -1 \\ 1 & -1 & 0 & 2 \\ 3 & 2 & -1 & 0 \\ \end{array}\right] \end{equation*}
and find its \(LU\) factorization. Explain why \(A\) and \(U\) have the same null space and use this observation to find a basis for \(Nul(A)\text{.}\)
Summary
We returned to Gaussian elimination, which we have used as a primary tool for finding solutions to linear systems, and explored its practicality as a computation tool, both in terms of numerical accuracy and computational effort.
- We saw that the accuracy of computations implemented on a computer could be improved using partial pivoting , a technique that performs row interchanges so that the entry in a pivot position has the largest possible magnitude.
- Beginning with a matrix \(A\text{,}\) we used the Gaussian elimination algorithm to write \(PA = LU\text{,}\) where \(P\) is a permutation matrix, \(L\) is lower triangular, and \(U\) is upper triangular.
- Finding this factorization involves roughly as much work as performing Gaussian elimination. However, once we have the factorization, we are able to quickly solve equations of the form \(A\mathbf x = \mathbf b\) by first solving \(L\mathbf c = P\mathbf b\) and then \(U\mathbf x = \mathbf c\text{.}\)
Exercises 5.1.4Exercises
In this section, we saw that errors made in computer arithmetic can produce approximate solutions that are far from the exact solutions. Here is another example in which this can happen. Consider the matrix
- Find the exact solution to the equation \(A\mathbf x = \twovec{2}{2}\text{.}\)
- Suppose that this linear system arises in the midst of a larger computation except that, due to some error in the computation of the right hand side of the equation, our computer thinks we want to solve \(A \mathbf{x}=\left[\begin{array}{c}2 \\ 2.0001\end{array}\right]\). Find the solution to this equation and compare it to the solution of the equation in the previous part of this exericse.
Notice how a small change in the right hand side of the equation leads to a large change in the solution. In this case, we say that the matrix \(A\) is ill-conditioned because the solutions are extremely sensitive to small changes in the right hand side of the equation. Though we will not do so here, it is possible to create a measure of the matrix that tells us when a matrix is ill-conditioned. Regrettably, there is not much we can do to remedy this problem.
In this section, we found the \(LU\) factorization of the matrix
in one of the activities, without using partial pivoting. Apply a sequence of row operations, now using partial pivoting, to find an upper triangular matrix \(U\) that is row equivalent to \(A\text{.}\)
In the following exercises, use the given \(LU\) factorizations to solve the equations \(A\mathbf x = \mathbf b\text{.}\)
-
Solve the equation
\begin{equation*} A\mathbf x = \left[\begin{array}{rr} 1 & 0 \\ -2 & 1 \\ \end{array}\right] \left[\begin{array}{rr} 3 & 1 \\ 0 & -2 \\ \end{array}\right]\mathbf x = \twovec{-3}{0}\text{.} \end{equation*}
-
Solve the equation
\begin{equation*} A\mathbf x = \left[\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 2 & 1 \\ \end{array}\right] \left[\begin{array}{rrr} 2 & 1 & 0 \\ 0 & -1 & 3 \\ 0 & 0 & 1 \\ \end{array}\right]\mathbf x = \threevec{5}{-5}{7}\text{.} \end{equation*}
Use Sage to solve the following equation by finding an \(LU\) factorization:
Here is another problem with approximate computer arithmetic that we will encounter in the next section. Consider the matrix
- Notice that this is a positive stochastic matrix. What do we know about the eigenvalues of this matrix?
-
Use Sage to define the matrix \(A\) using decimals such as
0.2
and the \(3\times3\) identity matrix \(I\text{.}\) Ask Sage to compute \(B = A-I\) and find the reduced row echelon form of \(B\text{.}\) - Why is the computation that Sage performed incorrect?
- Explain why using a computer to find the eigenvectors of a matrix \(A\) by finding a basis for \(Nul(A-\lambda I)\) is problematic.
In practice, one rarely finds the inverse of a matrix \(A\text{.}\) It requires considerable effort to compute, and we can solve any equation of the form \(A\mathbf x = \mathbf b\) using an \(LU\) factorization, which means that the inverse isn't necessary. In any case, the best way to compute an inverse is using an \(LU\) factorization, as this exericse demonstrates.
-
Suppose that \(PA = LU\text{.}\) Explain why \(A^{-1} = U^{-1}L^{-1}P\text{.}\)
Since \(L\) and \(U\) are triangular, finding their inverses is relatively efficient. That makes this an effective means of finding \(A^{-1}\text{.}\)
-
Consider the matrix
\begin{equation*} A = \left[\begin{array}{rrr} 3 & 4 & -1 \\ 2 & 4 & 1 \\ -3 & 1 & 4 \\ \end{array}\right]\text{.} \end{equation*}
Find the \(LU\) factorization of \(A\) and use it to find \(A^{-1}\text{.}\)
Consider the matrix
- Find the \(LU\) factorization of \(A\text{.}\)
- What conditions on \(a\text{,}\) \(b\text{,}\) \(c\text{,}\) and \(d\) guarantee that \(A\) is invertible?
In the \(LU\) factorization of a matrix, the diagonal entries of \(L\) are all \(1\) while the diagonal entries of \(U\) are not necessarily \(1\text{.}\) This exercise will explore that observation by considering the matrix
- Perform Gaussian elimination without partial pivoting to find \(U\text{,}\) an upper triangular matrix that is row equivalent to \(A\text{.}\)
- The diagonal entries of \(U\) are called pivots . Explain why \(\det A\) equals the product of the pivots.
- What is \(\det A\) for our matrix \(A\text{?}\)
- More generally, if we have \(PA=LU\text{,}\) explain why \(\det A\) equals plus or minus the product of the pivots.
Please provide a justification to your responses to these questions.
- In this section, our hypothetical computer could only store numbers using 3 decimal places. Most computers can store numbers using 15 or more decimal places. Why do we still need to be concerned about the accuracy of our computations when solving systems of linear equations?
- Finding the \(LU\) factorization of a matrix \(A\) is roughly the same amount of work as finding its reduced row echelon form. Why is the \(LU\) factorization useful then?
- How can we detect whether a matrix is invertible from its \(LU\) factorization?
Consider the matrix
- Find the \(LU\) factorization of \(A\text{.}\)
- Use the factorization to find a basis for \(Nul(A)\text{.}\)
- We have seen that \(Nul(A) = Nul(U)\text{.}\) Is it true that \(Col(A) = Col(L)\text{?}\)