We are going to explain Gaussian elimination on an example. Suppose we have the following linear system comprising three equations.
We can solve this system by using the first equation to eliminate x1 from the remaining equations. This yields a new system of two equations with two unknowns. By applying the procedure recursively we end up with one equation and one unknown. The procedure is referred to as Gaussian elimination. In our example we start it by multiplying the first equation with -2/1 (i.e. the negative of the coefficient corresponding to x1 in the second equation divided by the coefficient corresponding to x1 in the first equation) and adding it to the second equation. This eliminates x1 from the second equation resulting in
Next, we multiply the first equation by -(-6)/1=6 and add it to the third equation.
At this point we eliminated x1 from the second and third equation. These two equations now contain only two unknowns. Next, we repeat the procedure and eliminate x2 from the third equation by multiplying the second equation with -4/(-1)=4 and adding it to the third equation.
A more convenient way for representing a system of equations is to use the form of an extended matrix where rows correspond to the equations, the first 3 columns correspond to the coefficients in front of the unknowns, and the fourth column corresponds to the right-hand sides of the equations.
The matrix obtained after Gaussian elimination has zeros below the main diagonal and is also referred to as the row echelon form. Suppose our system of equations has n equations. Gaussian elimination requires n(n-1)/2 divisions, (n-1)n(n+1)/3 multiplications, and (n-1)n(n+1)/3 additions. Roughly speaking the computational burden of a Gaussian elimination grows proportionally to n3.
From row echelon form we can quickly obtain the solution of the linear system by applying back-substitution. First we solve the last equation which has only one unknown (x3). This results in x3=-72/25. We substitute this into the second equation to eliminate x3 which leaves us with only one unknown (x2).
From here we obtain x2=217/25. Finally, we substitute x2 and x3 into the first equation which yields.
resulting in x1=-393/25. Now we are finished and the system is solved.
Assuming we have n equations back-substitution requires n(n-1)/2 multiplications, n(n-1)/2 subtractions, and n divisions. The computational complexity of back-substitution grows proportionally to n2.
Gaussian elimination can handle systems of several thousands of equations. Larger systems of equations can be solved with iterative methods. Also note that large systems require a lot of memory for storing the matrix of coefficients. Take for instance n=10000. Matrices of this size have 108 elements which (if double precision is used) require 800MB of memory. On the other hand, solving such a large system requires roughly n3/3=0.33 1012 multiplications. With a processor running at 3.3GHz and performing one multiplication per clock cycle the multiplications alone require over 100 seconds of CPU time.
Rounding errors can significantly affect the result obtained with Gaussian elimination and back-substitution. Take for instance the system of equations represented by the following extended matrix
The solution of this system is x1=1, x2=2, and x3=-3. Due to finite precision the result of every operation is rounded. For double floating point precision this rounding takes place at the 16-th significant decimal digit. To illustrate the problem we perform Gaussian elimination and round every obtained matrix entry to 3 significant digits. First we eliminate nonzero entries below the first diagonal element.
Next, we eliminate the nonzero entry below the second diagonal element.
From here we obtain x3 with the first step of back-substitution which yields x3=0.00316. This is way off from what we expected to get (-3). What went wrong? Obviously numerical error accumulated throughout Gaussian elimination and ruined the result.
Suppose we are solving a linear system of equations
Numerical errors caused by rounding actually make us solve a perturbed system of equations. We assume the right-hand side (b) is not perturbed. A complete error analysis would also include the latter but the final result would be similar to the one we state here.
The following bound can be obtained for the relative error of the solution
All norms in this equation are L2 norms (i.e. for vectors this is the Euclidean norm). For matrices the L2 norm can be obtained as
The condition number of matrix A is defined as
The relative error of the solution depends on the condition number of matrix A and the matrix perturbation caused by rounding. We cannot do anything to change the condition number of matrix A. But we can make sure the matrix perturbation δA remains small.
When we are eliminating element aki in the k-th row by adding the scaled i-th row to the k-th row the scaling factor is determined by -aki/aii. When the absolute value of this factor is large the added entries of the scaled i-th row "drown" the much smaller entries of the k-th row in numerical noise. We can avoid this problem if we keep the absolute value of aki/aii as low as possible by making sure the absolute value of aii is as large as possible.
This can be achieved if we introduce row swapping. Swapping two rows is equivalent to swapping two equations. Clearly the final result is not affected. We swap the i-th row with the one that has the largest absolute entry in the i-th column between elements aii and ani. By swapping rows we can minimize the absolute value of the scaling factor -aki/aii. The entry that replaces aii after we swap rows is also referred to as the pivot. The aforementioned procedure is called partial pivoting.
Let us apply Gaussian elimination with partial pivoting to see if we do any better. Again, we are going to round obtained matrix entries to 3 significant digits. The pivot in the first column is found in the third row. Therefore we swap the first and the third row and start Gaussian elimination on the following extended matrix
After eliminating entries below the first diagonal element we get.
Before eliminating entries below the second diagonal element we swap the second and the third row so that a32 becomes the new pivot in the second column.
The next step completes the Gaussian elimination and yields
The first step of back-substitution now gives us x3=-2.97 which is close to the correct value (-3) and much better than the result obtained without pivoting.
Note that there also exists a more elaborate procedure called complete pivoting where the pivot is chosen between matrix entries in the rectangle between aii and ann. In this case we not only swap rows, but also swaps columns. Swapping columns corresponds to swapping unknowns. Therefore the values of the unknowns obtained after the back-substitution must be swapped in the opposite manner as columns were swapped during complete pivoting to obtain the correct result.
Often we need to solve multiple systems of equations with identical coefficient matrices which differ only in the value of the RHS vector. In such cases we use matrix decomposition. A matrix can be decomposed in a product of two matrices in many ways. A commonly used decomposition is the LU decomposition where we express matrix A as
Matrix L is a lower triangular matrix with zeros above its main diagonal and U is an upper triangular matrix with zeros below its main diagonal. The diagonal entries of matrix L are equal to 1. Once the LU decomposition of a matrix is known solving a linear system is fairly cheap in computational sense. Suppose we are solving a linear system expressed in matrix form as Ax=b. We can write
where z is a vector. To solve the linear system we first solve Lz=b. Because L is lower triangular we can apply a procedure referred to as forward substitution which is similar to back-substitution, except that now we start by expressing the first component of z from the first equation. We substitute this into the second equation which gives us the second component of z. By repeating this procedure we obtain all components of z in the same number of operations as are required for one back-substitution (i.e. the computational complexity grows proportionally to n2). In fact we don't even have to perform division because the diagonal entries of L are all equal to 1.
Because z=Ux by definition and U is upper triangular we can solve for x by applying back-substitution to Ux=z. The computational complexity of this step also grows with n2. To conclude, the computational complexity of solving a linear system of equations when the LU decomposition is known in advance is proportional to n2. This is, of course, much cheaper than performing a Gaussian elimination.
Take, for instance, the following linear system Ax=b where the LU decomposition of A is
and
We start by solving Lz=b which yields
and
Finally, we solve Ux=z.
The result is
The LU decomposition of a matrix can be computed as a byproduct of Gaussian elimination. The upper triangular part of the matrix (including the diagonal entries) obtained after Gaussian elimination is the U matrix. The L matrix is composed of ones on the main diagonal while subdiagonal entries (lij) are the negatives of the factors we used for multiplying the row containing the pivot (aii) when we were eliminating element aij. To demonstrate the algorithm let us perform Gaussian elimination side by side with LU decomposition. Note that this time we are going to work with matrix A instead of the extended matrix so there will only be n columns.
Take, for instance, the following matrix A. At this point we know little of L and U.
To eliminate the entries in the first column below the main diagonal we multiply the first row with -2/1 and add it to the second row. This eliminates the entry in the first column of the second row. The entry in the first column of the second row of matrix L is equal to the negative of the multiplier, i.e. 2/1. Similarly we multiply the first row with -(-6)/1 and add it to the third row to eliminate the entry in the first column of the third row. The element of L in the first column of the third row is therefore -6/1=-6. After we eliminate all subdiagonal entries in the first column we can copy the diagonal element of the first row along with the ones lying to its right into matrix U.
Next, we eliminate the subdiagonal entries in the second column. There is only one such entry. We eliminate it by multiplying the second row with -4/(-1)=4 and add the result to the third row. The corresponding element of matrix L is therefore -4. Because at this point the elimination of subdiagonal entries in the second column is complete we can copy the diagonal element of the second row along with all elements to its right into matrix U.
Finally, there are no subdiagonal elements to eliminate in the last (third) column. Therefore we copy the diagonal element from the third row into matrix U. The LU decomposition is complete.
Note how the matrix we obtained after Gaussian elimination is actually the U matrix.
In our last example we did not use partial pivoting. In order to avoid large numerical errors partial pivoting is necessary. When performing LU decomposition with partial pivoting we must record the row exchanges that took place during LU decomposition. When solving for z from LU=b' we use b' instead of b where b' is obtained by exchanging the elements of b in the same way we exchanged rows during LU decomposition.