Skip to main content
Logo image

Section 5.1 Gaussian elimination revisited

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 Ax=b for several different vectors b.

Preview Activity 5.1.1.

To begin, let’s recall how we implemented Gaussian elimination by considering the matrix
A=[12βˆ’1210βˆ’213210]
  1. What is the first row operation we perform? If the resulting matrix is A1, find a matrix E1 such that E1A=A1.
  2. What is the matrix inverse E1βˆ’1? 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.
  3. Perform the next two steps in the Gaussian elimination algorithm to obtain A3. Represent these steps using multiplication by matrices E2 and E3 so that
    E3E2E1A=A3.
  4. Suppose we need to scale the second row by βˆ’2. What is the 3Γ—3 matrix that perfoms this row operation by left multiplication?
  5. Suppose that we need to interchange the first and second rows. What is the 3Γ—3 matrix that performs this row operation by left multiplication?

Subsection 5.1.1 Partial pivoting

The first issue that we address is the fact that computers do not perform arithemtic operations exactly. For instance, Python will evaluate 0.1 + 0.2 and report 0.30000000000000004 even 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, must be converted to base 2. Even though 0.1 has a simple decimal form, its representation in base 2 is the repeating decimal
0.000110011001100110011001100110011001100110011….,
To accurately represent this number inside a computer would require infinitely many digits. Since a computer can only hold a finite number of digits, we are necessarily using an approximation just by 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Γ—n matrix, requires approximately 23n3 operations. If we have a 1000Γ—1000 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. We will consider the linear system
0.0001x+y=1x+y=2.
  1. Show that this system has the unique solution
    x=100009999=1.00010001…,y=99989999=0.99989998….
  2. 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.
  3. 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?
  4. By substituting our value for y into the first equation, what do we find for x?
  5. Compare the solution we find on our computer with the actual solution and assess the quality of the approximation.
  6. Let’s now modify the linear system by simplying interchanging the equations:
    x+y=20.0001x+y=1.
    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?
  7. Substitute the value you found for y into the first equation and solve for x. Then compare the approximate solution found with our hypothetical computer to the exact solution.
  8. 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
x=1.00010001…,y=0.99989998….
Since our hypothetical computer represents numbers using only three decimal digits, our computer finds
xβ‰ˆ1.00,yβ‰ˆ1.00.
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
(1βˆ’10000)y=2βˆ’10000βˆ’9999y=βˆ’9998βˆ’10000yβ‰ˆβˆ’10000yβ‰ˆ1.00.
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
(1βˆ’0.0001)y=2βˆ’0.0001βˆ’0.9999y=βˆ’0.9998βˆ’yβ‰ˆβˆ’1.00yβ‰ˆ1.00.
The difference occurs when we substitute yβ‰ˆ1 into the first equation. When the equations are written in their original order, we have
0.0001x+1.00β‰ˆ1.000.0001xβ‰ˆ0.00xβ‰ˆ0.00.
When the equations are written in their original order, we find the solution xβ‰ˆ0.00,yβ‰ˆ1.00.
When we write the equation in the opposite order, however, substituting yβ‰ˆ1 into the first equation gives
x+1.00β‰ˆ2.00xβ‰ˆ1.00.
In this case, we find the approximate solution xβ‰ˆ1.00,yβ‰ˆ1.00, 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β‰ˆ1 into the equation 0.00001x+y=1, which is a nearly horizontal line. Along this line, a small change in y leads to a large change in x. The slight difference in our approximation yβ‰ˆ1 from the exact value y=0.9998999… leads to a large difference in the approximation xβ‰ˆ0 from the exact value x=1.00010001….
If we exchange the order in which the equations are written, we substitute our approximation yβ‰ˆ1 into the equation x+y=2. Notice that the slope of the associated line is βˆ’1. 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β‰ˆ1 from the exact value leads to only a small difference in the approximation xβ‰ˆ1 from the exact value.
This example motivates the technique that computers usually use to perform Gaussian elimation. We only need to perform a row interchange when a zero occurs in a pivot position, such as
[1βˆ’12200βˆ’31022βˆ’3].
However, we will perform a row interchange to put the entry having the largest possible absolute value 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.
[21231βˆ’3βˆ’21βˆ’323βˆ’2]∼[βˆ’323βˆ’21βˆ’3βˆ’212123].
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.

Subsection 5.1.2 LU factorizations

In Subsection 1.3.3, we saw that the number of arithmetic operations needed to perform Gaussian elimination on an nΓ—n matrix is about 23n3. This means that a 1000Γ—1000 matrix, requires about two thirds of a billion operations.
Suppose that we have two equations, Ax=b1 and Ax=b2, that we would like to solve. Usually, we would form augmented matrices [Ab1] and [Ab2] 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 [Ab1] 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
A=[121βˆ’2βˆ’3βˆ’2374]
and begin performing Gaussian elimination without using partial pivoting.
  1. Perform two row replacement operations to find the row equivalent matrix
    Aβ€²=[121010011].
    Find elementary matrices E1 and E2 that perform these two operations so that E2E1A=Aβ€².
  2. Perform a third row replacement to find the upper triangular matrix
    U=[121010001].
    Find the elementary matrix E3 such that E3E2E1A=U.
  3. We can write A=E1βˆ’1E2βˆ’1E3βˆ’1U. Find the inverse matrices E1βˆ’1, E2βˆ’1, and E3βˆ’1 and the product L=E1βˆ’1E2βˆ’1E3βˆ’1. Then verify that A=LU.
  4. Suppose that we want to solve the equation Ax=b=[4βˆ’712]. We will write
    Ax=LUx=L(Ux)=b
    and introduce an unknown vector c such that Ux=c. Find c by noting that Lc=b and solving this equation.
  5. Now that we have found c, find x by solving Ux=c.
  6. Using the factorization A=LU and this two-step process, solve the equation Ax=[2βˆ’27].
This activity introduces a method for factoring a matrix A as a product of two triangular matrices, A=LU, 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.

Example 5.1.1.

Suppose we have the equation
[2βˆ’31βˆ’4502βˆ’22]x=[8βˆ’138],
which we write in the form Ax=b. We begin by applying the Gaussian elimination algorithm to find an LU factorization of A.
The first step is to multiply the first row of A by 2 and add it to the second row. The elementary matrix
E1=[100210000]
performs this operation so that E1A=[2βˆ’310βˆ’122βˆ’22].
We next apply matrices
E2=[100010βˆ’101],   E3=[100010011]
to obtain the upper triangular matrix U=E3E2E1A=[2βˆ’310βˆ’12003].
We can write U=(E3E2E1)A, which tells us that
A=(E3E2E1)βˆ’1U=[100βˆ’2101βˆ’11]U=LU.
That is, we have
A=LU=[100βˆ’2101βˆ’11][2βˆ’310βˆ’12003].
Notice that the matrix L is lower triangular, a result of the fact that the elementary matrices E1, E2, and E3 are lower triangular.
Now that we have factored A=LU into two triangular matrices, we can solve the equation Ax=b by solving two triangular systems. We write
Ax=L(Ux)=b
and define the unknown vector c=Ux, which is determined by the equation Lc=b. Because L is lower triangular, we find the solution using forward substitution, c=[833]. Finally, we find x, the solution to our original system Ax=b, by applying back substitution to solve Ux=c. This gives x=[2βˆ’11].
If we want to solve Ax=b for a different right-hand side b, we can simply repeat this two-step process.
An LU factorization allow us to trade in one equation Ax=b for two simpler equations
Lc=bUx=c.
For instance, the equation Lc=b in our example has the form
[100βˆ’2101βˆ’11]c=[8βˆ’138].
Because L is a lower-triangular matrix, we can read off the first component of c directly from the equations: c1=8. We then have βˆ’2c1+c2=βˆ’13, which gives c2=3, and c1βˆ’c2+c3=8, which gives c3=3. Solving a triangular system is simplified because we only need to perform a sequence of substitutions.
In fact, solving an equation with an nΓ—n triangular matrix requires approximately 12n2 operations. Once we have the factorization A=LU, we solve the equation Ax=b by solving two equations involving triangular matrices, which requires about n2 operations. For example, if A is a 1000Γ—1000 matrix, we solve the equation Ax=b using about one million steps. The compares with roughly a billion operations needed to perform Gaussian elimination, which 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 Ax=b for different right hand sides b.
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
P=[001010100],
which has the effect of interchanging the first and third rows. Notice that this matrix is not triangular so performing a row interchange will disrupt the structure of the LU factorization we seek. Without giving the details, we simply note that 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, where P is a permutation matrix, L is lower triangular, and U is upper triangular.
Therefore, to solve the equation Ax=b, we first multiply both sides by P to obtain
PAx=LUx=Pb.
That is, we multiply b by P and then find x using the factorization: Lc=Pb and Ux=c.

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, L, and U such that PA=LU.
  1. In Example 5.1.1, we found the LU factorization
    A=[2βˆ’31βˆ’4502βˆ’22]=[100βˆ’2101βˆ’11][2βˆ’310βˆ’12003]=LU.
    Using Sage, define the matrix A, and then ask Sage for the LU factorization. What are the matrices P, L, and U?
    Notice that Sage finds a different LU factorization than we found in the previous example because Sage uses partial pivoting, as described in the previous section, when it performs Gaussian elimination.
  2. Define the vector b=[8βˆ’138] in Sage and compute Pb.
  3. Use the matrices L and U to solve Lc=Pb and Ux=c. You should find the same solution x that we found in the previous example.
  4. Use the factorization to solve the equation Ax=[9βˆ’1610].
  5. How does the factorization show us that A is invertible and that, therefore, every equation Ax=b has a unique solution?
  6. Suppose that we have the matrix
    B=[3βˆ’122βˆ’11213].
    Use Sage to find the LU factorization. Explain how the factorization shows that B is not invertible.
  7. Consider the matrix
    C=[βˆ’212βˆ’11βˆ’10232βˆ’10]
    and find its LU factorization. Explain why C and U have the same null space and use this observation to find a basis for Nul(A).

Subsection 5.1.3 Summary

We returned to Gaussian elimination, which we have used as a primary tool for finding solutions to linear systems, and explored its practicality, 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, we used the Gaussian elimination algorithm to write PA=LU, 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 Ax=b by first solving Lc=Pb and then Ux=c.

Exercises 5.1.4 Exercises

1.

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
A=[1111.0001].
  1. Find the exact solution to the equation Ax=[22].
  2. 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 Ax=[22.0001]. Find the solution to this equation and compare it to the solution of the equation in the previous part of this exercise.
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.

2.

In this section, we found the LU factorization of the matrix
A=[121βˆ’2βˆ’3βˆ’2374]
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.

3.

In the following exercises, use the given LU factorizations to solve the equations Ax=b.
  1. Solve the equation
    Ax=[10βˆ’21][310βˆ’2]x=[βˆ’30].
  2. Solve the equation
    Ax=[100βˆ’210βˆ’121][2100βˆ’13001]x=[5βˆ’57].

4.

Use Sage to solve the following equation by finding an LU factorization:
[34βˆ’1241βˆ’314]x=[βˆ’3βˆ’3βˆ’4].

5.

Here is another problem with approximate computer arithmetic that we will encounter in the next section. Consider the matrix
A=[0.20.20.40.20.30.10.60.50.5].
  1. Notice that this is a positive stochastic matrix. What do we know about the eigenvalues of this matrix?
  2. Use Sage to define the matrix A using decimals such as 0.2 and the 3Γ—3 identity matrix I. Ask Sage to compute B=Aβˆ’I and find the reduced row echelon form of B.
  3. Why is the computation that Sage performed incorrect?
  4. Explain why using a computer to find the eigenvectors of a matrix A by finding a basis for Nul(Aβˆ’Ξ»I) is problematic.

6.

In practice, one rarely finds the inverse of a matrix A. It requires considerable effort to compute, and we can solve any equation of the form Ax=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 exercise demonstrates.
  1. Suppose that PA=LU. Explain why Aβˆ’1=Uβˆ’1Lβˆ’1P.
    Since L and U are triangular, finding their inverses is relatively efficient. That makes this an effective means of finding Aβˆ’1.
  2. Consider the matrix
    A=[34βˆ’1241βˆ’314].
    Find the LU factorization of A and use it to find Aβˆ’1.

7.

Consider the matrix
A=[aaaaabbbabccabcd].
  1. Find the LU factorization of A.
  2. What conditions on a, b, c, and d guarantee that A is invertible?

8.

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. This exercise will explore that observation by considering the matrix
A=[311βˆ’6βˆ’4βˆ’10βˆ’41].
  1. Perform Gaussian elimination without partial pivoting to find U, an upper triangular matrix that is row equivalent to A.
  2. The diagonal entries of U are called pivots. Explain why detA equals the product of the pivots.
  3. What is detA for our matrix A?
  4. More generally, if we have PA=LU, explain why detA equals plus or minus the product of the pivots.

9.

Please provide a justification to your responses to these questions.
  1. 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?
  2. 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?
  3. How can we detect whether a matrix is invertible from its LU factorization?

10.

Consider the matrix
A=[βˆ’11001βˆ’21001βˆ’2100βˆ’11].
  1. Find the LU factorization of A.
  2. Use the factorization to find a basis for Nul(A).
  3. We have seen that Nul(A)=Nul(U). Is it true that Col(A)=Col(L)?
You have attempted 1 of 1 activities on this page.