Gaussian Elimination

Comp.145  November 1996  Dr. A.J. Mestel

Let A be an n × n non-singular matrix with elements a_{ij}, and let \underline{b} be an n × 1 vector. The simultaneous equations A\underline{x} = \underline{b} can be solved for the n unknowns x_{1}, x_{2}, … x_{n} by a systematic process of row operations known as Gaussian elimination. For illustration, we consider the case n = 4. The equations are then:

\begin{aligned}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} + a_{14}x_{4} & = b_{1}  (1) \\ a_{21}x_{1} + a_{22}x_{2} + a_{23}x_{3} + a_{24}x_{4} & = b_{2}  (2) \\ a_{31}x_{1} + a_{32}x_{2} + a_{33}x_{3} + a_{34}x_{4} & = b_{3}  (3) \\ a_{41}x_{1} + a_{42}x_{2} + a_{43}x_{3} + a_{44}x_{4} & = b_{4}  (4) \end{aligned}

Assuming a_{11} ≠ 0, we define for i = 2, 3, 4 the numbers L_{i1} = a_{i1}/a_{11} and then subtract L_{i1} × equation (1) from equation (i). We obtain a new set of equations:

\begin{aligned}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} + a_{14}x_{4} & = b_{1}  (1) \\ a′_{22}x_{2} + a′_{23}x_{3} + a′_{24}x_{4} & = b′_{2}  (2′) \\ a′_{32}x_{2} + a′_{33}x_{3} + a′_{34}x_{4} & = b′_{3}  (3′) \\ a′_{42}x_{2} + a′_{43}x_{3} + a′_{44}x_{4} & = b′_{4}  (4′) \end{aligned}  where \{ \begin{matrix}a′_{ij} & = a_{ij} − L_{i1}a_{1j} \\ b′_{i} & = b_{i} − L_{i1}b_{1} \\ i, j & = 2, 3, 4 \end{matrix} .

The first equation is unchanged, while the others form a reduced system of dimension (n − 1). We now assume a′_{22} ≠ 0, define L_{i2} = a′_{i2}/a′_{22} and subtract L_{i2} × equation (2′) from equation (i′) for i = 3, 4. The new equations are

\begin{aligned}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} + a_{14}x_{4} & = b_{1}  (1) \\ a′_{22}x_{2} + a′_{23}x_{3} + a′_{24}x_{4} & = b′_{2}  (2′) \\ a′′_{33}x_{3} + a′′_{34}x_{4} & = b′′_{3}  (3′′) \\ a′′_{43}x_{3} + a′′_{44}x_{4} & = b′′_{4}  (4′′) \end{aligned}  where  \{ \begin{matrix}a′′_{ij} & = a′_{ij} − L_{i2}a′_{2j} \\ b′′_{i} & = b′_{i} − L_{i2}b′_{2} \\ i, j & = 3, 4 \end{matrix} .

Finally, we assume a′′_{33} ≠ 0, let L_{43} = a′′_{43}/a′′_{33}, and subtract L_{43} × (3′′) from (4′′), giving

\begin{aligned}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} + a_{14}x_{4} & = b_{1}  (1) \\ a′_{22}x_{2} + a′_{23}x_{3} + a′_{24}x_{4} & = b′_{2}  (2′) \\ a′′_{33}x_{3} + a′′_{34}x_{4} & = b′′_{3}  (3′′) \\ a′′′_{44}x_{4} & = b′′′_{4}  (4′′′) \end{aligned}  where  \{ \begin{matrix}a′′′_{44} & = a′′_{44} − L_{43}a′′_{34} \\ b′′′_{4} & = b′′_{4} − L_{43}b′′_{3} \end{matrix} .

This completes the elimination stage of the algorithm, which reduces the system to triangular form. These equations are easily solved by back-substitution. We have already assumed a_{11}, a′_{22}, a′′_{33} ≠ 0. Provided also a′′′_{44} ≠ 0, equation (4′′′) determines x_{4}, which we can substitute into (3′′) to find x_{3} and so on until all the unknowns are found.

(1) The Gaussian Elimination Algorithm is equivalent to “factorising” the matrix A into Lower and Upper triangular matrices, a process known as LU-decomposition. That is A = LU, where (with the matrix elements as defined overleaf)

L = ( \begin{matrix}1 & 0 & 0 & 0 \\ L_{21} & 1 & 0 & 0 \\ L_{31} & L_{32} & 1 & 0 \\ L_{41} & L_{42} & L_{43} & 1 \end{matrix} )  and  U = ( \begin{matrix}a_{11} & a_{12} & a_{13} & a_{14} \\ 0 & a′_{22} & a′_{23} & a′_{24} \\ 0 & 0 & a′′_{33} & a′′_{34} \\ 0 & 0 & 0 & a′′′_{44} \end{matrix} ) .

This is not obvious! Note that detA = (detL)(detU) = (1)(a_{11}a′_{22}a′′_{33}a′′′_{44}) ≠ 0.

(2) At any stage in the process, the equation being subtracted from the others is known as the pivotal equation, and its first coefficient is known as the pivot. If at any stage in the process the natural pivot (say, a′_{22}) is zero then it is necessary to reorder the reduced system of equations so as to obtain a non-zero pivot. If it is not possible to do this then necessarily detA = 0, and the matrix is singular.

(3) If one of the pivots is small, but not zero, the algorithm will work, but computational errors can become unacceptably large. In general, it is a good idea to choose a pivot of as large a modulus as possible. A strategy known as partial pivoting consists of choosing the greatest coefficient in the first column of the reduced matrix to act as pivot. Complete pivoting involves choosing the greatest coefficient in the entire reduced matrix (i.e. involves swapping columns as well as rows.)

(4) Gaussian Elimination can be used to solve the equations A\underline{x} = \underline{b} for several different RHS vectors \underline{b} simultaneously. In particular, we can calculate the inverse matrix A^{−1} by considering n different RHS. If we let \underline{b} = (1, 0, 0, … 0)^{T}, we obtain for \underline{x} the first column of A^{−1}. The other columns can be found similarly by letting \underline{b} = (0, … 1, 0 … 0)^{T}.

(5) The number of floating-point operations (flops) required to solve an n-dimensional system of equations with p different RHS is approximately \frac{1}{3}n^{3} + pn^{2}. We may therefore find A^{−1} in \frac{4}{3}n^{3} flops. In fact, we can do better than this by exploiting the fact that the n RHS-vectors contain many zeros. Matrix inversion can be achieved in n^{3} flops.

(6) The following piece of BASIC program solves the n equations A\underline{x} = \underline{b} provided no row-swapping is required. The program overwrites the arrays a(1 : n, 1 : n) and b(1 : n).

FOR k=1 TO n–1  The pivotal row is number k.
FOR i=k+1 TO n  For each remaining row
e=a(i,k)/a(k,k) : a(i,k)=0  calculate the multiplier
FOR j=k+1 TO n : a(i,j)=a(i,j)–e∗a(k,j) : NEXT j and subtract the pivotal row.
b(i)=b(i)–e∗b(k)  Then deal with the RHS.
NEXT i  Do the above for all rows and
NEXT k  commence back substitution:
x(n)=b(n)/a(n,n)  Find x(n) first,
FOR i=n–1 TO 1 STEP –1  then x(n-1) etc.
y=b(i)  as the RHS, plus
FOR j=n TO i+1 STEP –1 : y=y–a(i,j)∗x(j) : NEXT j multiples of the known x’s
x(i)=y/a(i,i)  divided by the pivots.
NEXT i  That’s all, folks!