Gaussian Elimination

The Gaussian elimination method is one of the most important and ubiquitous algorithms that can help deduce important information about the given matrix�s roots/nature as well determine the solvability of linear system when it is applied to the augmented matrix. As such, it is one of the most useful numerical algorithms and plays a fundamental role in scientific computation. This method has been historically around and used by Chinese mathematicians since, 179 CE. Many people have contributed to Gaussian elimination, including Isaac Newton. However, it was named by the American mathematician George Forsythe (1917--1972) in honor of the German mathematician and physicist Carl Friedrich Gauss (1777--1855). The full history of Gaussian algorithm can be found within the Joseph F. Grcar article. For our exposition, we need a shorthand way of writing systems of equations.

For a linear system of equation Ax = b, where A is an m×n matrix and b is a known n column vector, we assign the m×(n+1) augmented matrix, denoted [A|b] or (A|b), which is A with column b tagged on.
In other words, if A is the m×n matrix
${\bf A} = \begin{bmatrix} a_{11}& a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21}& a_{22} & a_{23} & \cdots & a_{2n} \\ \vdots& \vdots& \vdots& \ddots & \vdots \\ a_{m1}& a_{m2} & a_{m3} & \cdots & a_{mn} \end{bmatrix} \qquad \mbox{and} \qquad {\bf b} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix} ,$
then the corresponding augmented matrix is
$\left[ {\bf A} \vert {\bf b} \right] = \left[ \begin{array}{ccccc|c} a_{11}& a_{12} & a_{13} & \cdots & a_{1n} & b_1 \\ a_{21}& a_{22} & a_{23} & \cdots & a_{2n} & b_2 \\ \vdots& \vdots& \vdots& \ddots & \vdots & \vdots \\ a_{m1}& a_{m2} & a_{m3} & \cdots & a_{mn} & b_n \end{array} \right] .$
Sometimes a vertical line is used to separate the coefficient entries from the constants can be dropped and then augmented matrix is written as [A b] . The augmented matrix provides a precise and concise information about a linear system of equations
$\begin{cases} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1 , \\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n &= b_2 , \\ \cdots \qquad \cdots \\ a_{m1} x_1 + a_{m2} x_2 + \cdots + a_{mn} x_n &= b_m . \\ \end{cases}$
A system of linear equations above, written in compact vector form Ax = b, is said to be consistent if it has at least one solution and inconsistent if there is no solution.

To understand the idea of Gaussian elimination algorithm, we shall consider a series of examples, starting with a two dimensional case.

Example: Consider the system of algebraic equations
\begin{align*} x -2\, y &= 1 , \\ 3\,x + 4\,y &= 13 . \end{align*}
If we multiply the first equation by -3 and add to the last one (which does not change the solution), we get
${\bf Before: } \quad \begin{split} x -2\, y &= 1 , \\ 3\,x + 4\,y &= 13 ; \end{split} \qquad {\bf After: } \quad \begin{split} x -2\, y &= 1 , \\ \qquad\qquad 10\,y &= 10 ; \end{split} \qquad \begin{split} \mbox{(multiply by 3 and subtract)} \\ (x\mbox{ has been eliminated)} \end{split}$
The first stage we accomplished is called forward elimination because it deleted one variable from consideration. Forward elimination produces an upper triangular system, which can be seen with matrix notation (which is called the augmented matrix).
${\bf Before: } \qquad \begin{bmatrix} \color{red}{1}&-2&1 \\ 3&\color{red}{4}&13 \end{bmatrix} \qquad\quad {\bf After: } \qquad \begin{bmatrix} \color{red}{1}&-2&1 \\ 0&\color{red}{10}&10 \end{bmatrix} .$
The last equation 10y = 10 reveals y = 1, and we go up the triangle to x = 3. This quick process is called back substitution. It is used for upper triangular systems of any size, after forward elimination is complete. We plot our equations with Mathematica.
Plot[{(x - 1)/2, (13 - 3*x)/4}, {x, -0.5, 4}, PlotStyle -> Thick] Plot[{(x - 1)/2, 1}, {x, -0.5, 4}, PlotStyle -> Thick]

Instead of eliminating the variable x, we can go after y. If we multiply the first equation by 2 and add to the second equation, we get

${\bf Before: } \quad \begin{split} x -2\, y &= 1 , \\ 3\,x + 4\,y &= 13 ; \end{split} \qquad {\bf After: } \quad \begin{split} x -2\, y &= 1 , \\ 5\,x \qquad &= 15 ; \end{split} \qquad \begin{split} \mbox{(multiply by 2 and add)} \\ (y\mbox{ has been eliminated)} \end{split}$
or in matrix form:
${\bf Before: } \qquad \begin{bmatrix} 1&\color{red}{-2}&1 \\ \color{red}{3}&4&13 \end{bmatrix} \qquad\quad {\bf After: } \qquad \begin{bmatrix} 1&\color{red}{-2}&1 \\ \color{red}{5}&0&15 \end{bmatrix} .$
vline = Line[{{3, -0.5}, {3, 1.5}}];
Plot[{(x - 1)/2}, {x, -0.5, 4}, Epilog -> {Directive[Thick, Orange], vline}, PlotStyle -> Thick]
Which variable should you eliminate, x or y? For a computer, it does not matter, but for humans it does because it is simpler. So mostly for educational purposes, we will follow tradition and we will eliminate variables from left to right in order to reduce the matrices into upper triangular form. However, remember that when you write code for practical calculations, it does not matter which variable you eliminate and in what order---the computer does not care!

When we used matrix form corresponding to the given system of equations, we marked with the color red for special positions in the corresponding augmented matrix because they are important to understand. These positions are usually referred to as pivots. ■

Example: Now we consider slightly different system of algebraic equations
\begin{align*} x -2\, y &= 1 , \\ 3\,x - 6\,y &= 13 . \end{align*}
Eliminating variable x by subtracting 3 times first equation from the second one, we obtain
\begin{align*} x -2\, y &= 1 , \\ 0\,x - 0\,y &= \color{red}{10} . \end{align*}
There is no solution. Remember that zero is never allowed as a pivot, hence we get an equation with a pivot at the last column.
${\bf Before: } \qquad \begin{bmatrix} \color{red}{1}&-2&1 \\ 3&-6&\color{red}{13} \end{bmatrix} \qquad\quad {\bf After: } \qquad \begin{bmatrix} \color{red}{1}&-2&1 \\ 0&0&\color{red}{10} \end{bmatrix} .$
Plot[{(x - 1)/2, (3*x-13)/6}, {x, -0.5, 4}, PlotStyle -> Thick]

Example: We transfer the system of algebraic equations to matrix form using its augmented matrix
$\begin{split} x -2\, y &= 1 , \\ 3\,x - 6\,y &= 3 . \end{split} \qquad \Longrightarrow \qquad \begin{bmatrix} \color{red}{1}&-2&1 \\ 3&-6&3 \end{bmatrix} .$
The result of subtracting 3 times first equation from the last one reveals the 2�2 matrix with only one pivot:
$\begin{bmatrix} \color{red}{1}&-2&1 \\ 0&0&0 \end{bmatrix} .$
The last line in the matrix shows that every x and y satisfy the equation 0·x + 0·y = 0. There is really only one equation x - 2y = 1. One of the variables is free when another one is expressed through the free one:
$x = 1 + 2\,y \qquad \mbox{or} \qquad y = \left( x-1 \right) /2 ,$
when y or x is freely chosen, respectively. There is no need to plot one straight line that includes both equations because every its point satisfies both equations. We have a whole line of solutions. ■

Example: Consider the system of algebraic equations
\begin{align*} 0\,x -2\, y &= -2 , \\ 3\,x + 4\,y &= 13 . \end{align*}
For this system, we cannot choose the first coefficient as a pivot because by definition the pivot cannot be zero. So exchange two equations to obtain an equivalent augmented matrix:
$\begin{bmatrix} 0&-2&-2 \\ 3&4&13 \end{bmatrix} \ \sim \ \begin{bmatrix} \color{red}{3}&4&13 \\ 0&\color{red}{-2}&-2 \end{bmatrix} .$
Plot[{(x - 1)/2, 1}, {x, -0.5, 4}, PlotStyle -> Thick]
The new system is already triangular, so one of the lines is parallel to an axis. ■

To understand Gaussian elimination, you have to go beyond 2�2 systems of equations. Therefore, we present examples of 3�3 systems of equations that will be enough to see the pattern. The other example, with rectangular matrices, is following.

Example: Consider the system of algebraic equations
\begin{align*} x -3\, y +z &= 6 , \\ -6\,x + 3\,y - 15\,z&= 3 , \\ 2\,x - 8\,y + 8\,z&= 10 . \end{align*}
We convert the given system of equations into an augmented matrix:
$\left[ \begin{array}{ccc|c} 1&-3&1&6 \\ -6&3&-15&3 \\ 2&-8&8&10 \end{array} \right] .$
This matrix contains all of the information in the system of equations without the unknown variables x, y, and z to carry around. Now we perform the process of elimination. The notation to the right of each matrix describes the row operations that were performed to get the matrix on that line. For example 2R1+R2 ↦ R2 means "replace row 2 with the sum of 2 times row 1 and row 2."
$\left[ \begin{array}{ccc|c} \color{red}{1}&-3&1&6 \\ 0&\color{red}{-15}&-9&39 \\ 0&-2&6&-2 \end{array} \right] \qquad \begin{array}{c} \\ 6\,R_1 + R_2 \mapsto R_2 \\ -2\,R_1 + R_3 \mapsto R_3 \end{array}$
$\left[ \begin{array}{ccc|c} \color{red}{1}&-3&1&6 \\ 0&\color{red}{-15}&-9&39 \\ 0&-2&6&-2 \end{array} \right] \,\sim \, \left[ \begin{array}{ccc|c} \color{red}{1}&-3&1&6 \\ 0&\color{red}{-5}&-3&13 \\ 0&1&\color{red}{-3}&1 \end{array} \right] \qquad \begin{array}{c} \\ \frac{1}{3}\,R_2 \mapsto R_2 \\ -\frac{1}{2}\,R_3 \mapsto R_3 \end{array} .$
$\left[ \begin{array}{ccc|c} \color{red}{1}&-3&1&6 \\ 0&\color{red}{-5}&-3&13 \\ 0&1&\color{red}{-3}&1 \end{array} \right] \,\sim \, \left[ \begin{array}{ccc|c} \color{red}{1}&-3&1&6 \\ 0&\color{red}{-5}&-3&13 \\ 0&0&\color{red}{-\frac{18}{5}}&\frac{18}{5} \end{array} \right] \qquad \begin{array}{c} \\ \\ \frac{1}{5}\,R_2 + R_3 \mapsto R_3 \end{array}$
If we now reverse the conversion process and turn the augmented matrix into a system of equations we have
\begin{align*} x -3\, y +z &= 6 , \\ 0\,x - 5\,y - 3\,z&= 2 , \\ 0\,x + 0\,y - \frac{18}{5}\,z&= \frac{18}{5} . \end{align*}
We can now easily solve for x, y, and z by back-substitution to obtain x = 1, y = -2, and z = -1. For a system of equations with a 3x3 matrix of coefficients, the goal of the process of Gaussian Elimination is to create (at least) a triangle of zeroes in the lower left hand corner of the matrix below the diagonal. Note that you may switch the order of the rows at any time in trying to get to this form.
a1 = ContourPlot3D[x - 3 y + z == 6, {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, ContourStyle -> Directive[Red]];

f2[x_, y_] := (-2* x + y-1)/5;
a2 = ParametricPlot3D[{x, y, f2[x, y]}, {x, -10, 10}, {y, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, PlotStyle -> Directive[Green]];

f3[x_, y_] := (5- x + 4*y)/4;
a3 = Graphics3D[{LightBlue, Polygon[Flatten[#, 1] &@{#[[1]], #[[2]],f3[#[[1]], #[[2]]]} & /@ {{-10, -10}, {-10, 10}, {10, 10}, {10, -10}}]}];
Show[a1, a2, a3]
Finally we plot planes in the echelon form:
ff2[x_, y_] := -(5*y + 2)/5;
b2 = ParametricPlot3D[{x, y, ff2[x, y]}, {x, -10, 10}, {y, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, PlotStyle -> Directive[Green]];
ff3[x_, y_] := -1;
b3 = Graphics3D[{LightBlue, Polygon[Flatten[#, 1] &@{#[[1]], #[[2]], ff3[#[[1]], #[[2]]]} & /@ {{-10, -10}, {-10, 10}, {10, 10}, {10, -10}}]}];
Show[a1, b2, b3]

Example: Consider the system of algebraic equations
\begin{align*} x -2\, y-6\,z &= 1 , \\ x -4\,y -12\,z &= 2, \\ 2\,x +4\,y + 12\,z &= 3 . \end{align*}
First, we plot planes corresponding the given system of equations:
a1 = ContourPlot3D[x - 2 y -6 z == 1, {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, ContourStyle -> Directive[Orange]];

f2[x_, y_] := (x -4 y-2)/12;
a2 = ParametricPlot3D[{x, y, f2[x, y]}, {x, -10, 10}, {y, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, PlotStyle -> Directive[Green]];

f3[x_, y_] := (3- 2*x - 4*y)/12;
a3 = Graphics3D[{LightBlue, Polygon[Flatten[#, 1] &@{#[[1]], #[[2]],f3[#[[1]], #[[2]]]} & /@ {{-10, -10}, {-10, 10}, {10, 10}, {10, -10}}]}];
Show[a1, a2, a3]
We keep x in the first equation and eliminate it from the other equations. To do so, add -1 times equation 1 to equation 2. After some practice, this type of calculation is usually performed mentally. However, it is convenient to use software:
$\begin{array}{c} -1 \cdot [\mbox{equation }1] \\ + [\mbox{equation }2] \\ \hline [\mbox{new equation }2] \end{array} \qquad \begin{array}{ccccc} -x&+2\,y& +6\,z& =&-1 \\ x&-4\,y& -12\,z& =&2 \\ \hline &-2\,y& -6\,z & =& 1 \end{array}$
To subtract row 1 from row 2, enter in Mathematica:
A = {{1, -2, -6, 1}, {1, -4, -12, 2}, {2, 4, 12, 3}}
m2 = A;
m2[[2]] -= m2[[1]]
To subtract row 1 times 2 from row 3 enter:
m3 = m2;
m3[[3]] = m3[[3]] - 2*m3[[1]]
$\begin{array}{c} -2 \cdot [\mbox{equation }1] \\ + [\mbox{equation }3] \\ \hline [\mbox{new equation }3] \end{array} \qquad \begin{array}{ccccc} -2\,x&+4\,y& +12\,z& =&-2 \\ 2\,x&+4\,y& +12\,z& =&3 \\ \hline &8\,y& +24\,z & =& 1 \end{array}$
$\begin{split} x-2\,y -6\,z &=1 , \\ -2\,y -6\,z &= 1 , \\ 8\,y +24\,z &= 1 ; \end{split} \qquad \Longleftrightarrow \qquad \left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 0&\color{red}{-2}&-6&1 \\ 0&8&24&1 \end{array} \right] .$
Now, multiply equation 2 by 4 and add to the last one. This yields
$\begin{split} x-2\,y -6\,z &=1 , \\ -2\,y -6\,z &= 1 , \\ 0\,y +0\,z &= 5 ; \end{split} \qquad \Longleftrightarrow \qquad \left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 0&\color{red}{-2}&-6&1 \\ 0&0&0&\color{red}{5} \end{array} \right] .$
Notice the last equation reads: 0=5. This is not possible. So the system has no solutions; it is not possible to find values x, y, and z that satisfy all three equations simultaneously. ■

Example: Consider the system of algebraic equations
\begin{align*} x -2\, y-6\,z &= 1 , \\ x -4\,y -12\,z &= 2, \\ 2\,x +4\,y + 12\,z &= -2 . \end{align*}
First, we plot planes corresponding the given system of equations:
a1 = ContourPlot3D[x - 2 y -6 z == 1, {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, ContourStyle -> Directive[Orange]];

f2[x_, y_] := (x -4 y-2)/12;
a2 = ParametricPlot3D[{x, y, f2[x, y]}, {x, -10, 10}, {y, -10, 10}, AxesLabel -> {x, y, z}, Mesh -> None, PlotStyle -> Directive[Green]];

f3[x_, y_] := (-2- 2*x - 4*y)/12;
a3 = Graphics3D[{LightBlue, Polygon[Flatten[#, 1] &@{#[[1]], #[[2]],f3[#[[1]], #[[2]]]} & /@ {{-10, -10}, {-10, 10}, {10, 10}, {10, -10}}]}];
Show[a1, a2, a3]
In augmented matrix form we have
$\left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 1&-4&-12&2 \\ 2&4&12&-2 \end{array} \right] .$
We now use the method of Gaussian Elimination and annihilate variable x from two last equations:
$\left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 0&\color{red}{-2}&--6&1 \\ 0&8&24&-4 \end{array} \right] \qquad \begin{array}{c} \\ -1\,R_1 + R_2 \mapsto R_2 \\ -2\,R_1 + R_3 \mapsto R_3 \end{array} .$
Next we divide every term in the last equation by 4 to obtain
$\left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 0&\color{red}{-2}&-6&1 \\ 0&2&6&-1 \end{array} \right] \qquad \begin{array}{c} \\ \\ \frac{1}{4}\, R_3 \mapsto R_3 \end{array} .$
Finally we add the second row and the third row. This leads to its reduced echelon form with only two pivots:
$\left[ \begin{array}{ccc|c} \color{red}{1}&-2&-6&1 \\ 0&\color{red}{-2}&-6&1 \\ 0&0&0&0 \end{array} \right]$
Converting back to a system of equations, we have
\begin{align*} x -2\, y-6\,z &= 1 , \\ 0\,x -2\,y -6\,z &= 1, \\ 0\,x +0\,y + 0\,z &= 0 . \end{align*}
Notice the last equation: 0=0 (this resulted from equation 3 being a linear combination of the other two equations). This is always true. And, we can solve the first two equations to get x and y as functions of z alone. Solving these equations, we get
Solve[{x - 2*y - 6*z == 1, -2*x - 6*z == 1}, {x, y}]
{{x -> 1/2 (-1 - 6 z), y -> -(3/4) (1 + 6 z)}}
So
$x= - \frac{1}{2} \left( 1+ 6\,x \right) , \qquad y = - \frac{3}{4} \left( 1+ 6\,x \right) ,$
where z is a free variable. ■

Example: Consider the system of algebraic equations
\begin{align*} x_1 -2\, x_2 + x_3 + x_4 &= -6 , \\ x_1 -4\,x_2 - x_3 + 3\,x_4 &= -12, \\ 2\,x_1 -3\,x_2 +2\,x_3 + x_4 &= -8 , \\ -3\,x_1 +2\,x_2 + x_3 + 2\,x_4 &= -4 . \end{align*}
The augmented matrix for this system becomes
$\left[ \begin{array}{cccc|c} 1&-2&1&1&-6 \\ 1&-4&-1&3&-12 \\ 2&-3&2&1&-8 \\ -3&2&1&2&-4 \end{array} \right] .$
To eliminate x1 from the last three equations, we multiply the first equation by -1, -2, and 3, respectively. By adding the results to the corresponding rows, it will introduce zeroes into positions below the pivot (which we mark with red) in the first column:
$\left[ \begin{array}{cccc|c} 1&-2&1&1&-6 \\ 1&-4&-1&3&-12 \\ 2&-3&2&1&-8 \\ -3&2&1&2&-4 \end{array} \right] \,\sim \, \left[ \begin{array}{cccc|c} \color{red}{1}&-2&1&1&-6 \\ 0&-2&-2&2&-6 \\ 0&1&0&-1&4 \\ 0&-4&4&5&-22 \end{array} \right] \qquad \begin{array}{c} \\ -1\,R_1 + R_2 \mapsto R_2 \\ -2\,R_1 + R_3 \mapsto R_3 \\ 3\,R_1 + R_4 \mapsto R_4 \end{array} .$
It is convenient to divide the second row by 2:
$\left[ \begin{array}{cccc|c} \color{red}{1}&-2&1&1&-6 \\ 0&-1&-1&1&-3 \\ 0&1&0&-1&4 \\ 0&-4&4&5&-22 \end{array} \right] \,\sim \left[ \begin{array}{cccc|c} \color{red}{1}&-2&1&1&-6 \\ 0&\color{red}{-1}&-1&1&-3 \\ 0&0&-1&0&1 \\ 0&0&8&1&-10 \end{array} \right] \qquad \begin{array}{c} \\ \\ R_2 + R_3 \mapsto R_3 \\ -4\,R_2 + R_4 \mapsto R_4 \end{array} .$
Finally, multiplying the third row by 8 and adding to the last row, we obtain, what is usually called the row echelon form for the given augmented matrix:
$\left[ \begin{array}{cccc|c} \color{red}{1}&-2&1&1&-6 \\ 0&\color{red}{-1}&-1&1&-3 \\ 0&0&\color{red}{-1}&0&1 \\ 0&0&0&\color{red}{1}&-2 \end{array} \right]$
that has pivots in every row. Return this to the familiar linear system:
\begin{align*} x_1 -2\, x_2 + x_3 + x_4 &= -6 , \\ -x_2 - x_3 + x_4 &= -3, \\ -x_3 &= 1 , \\ x_4 &= -2 . \end{align*}
Solving by back substitution, we obtain
$x_1 =1, \quad x_2 =2 , \quad x_3 =-1, \quad x_4 = -2 .$

The Gaussian elimination method is basically a series of operations carried out on a given matrix, in order to mathematically simplify it to its echelon form. When it is applied to solve a linear system Ax=b, it consists of two steps: forward elimination (also frequently called Gaussian elimination procedure) to reduce the matrix to upper triangular form, and back substitution. Therefore, solving a linear system of algebraic equations using this elimination procedure can also be called forward elimination and back substitution and abbreviated as FEBS.

A rectangular matrix is said to be in echelon form or row echelon form if it has the following three properties:
1. All nonzero rows are above any rows of all zeroes.
2. Each leading entry, called the pivot, of a row is in a column to the right of the leading entry of the row above it.
3. All entries in a column below a leading entry are zeroes.
An echelon matrix is one that is in echelon form.
A matrix is said to be in its row echelon form when it meets the following two conditions:

1. Any present zero rows are placed all the way at the bottom of the matrix.
2. The very first value of the matrix (termed as the leading entry or pivot) is a non-zero term (some texts prefer it be �1�, but it is not necessarily).
Example: The following matrices are in the echelon form:
$\begin{bmatrix} \color{red}{♦} & \ast & \ast & \ast \\ 0& \color{red}{♦} & \ast & \ast \\ 0&0&0&0 \\ 0&0&0&0 \end{bmatrix} , \quad \begin{bmatrix}\color{red}{ ♦} & \ast & \ast & \ast \\ 0&0 & \color{red}{♦} & \ast \\ 0&0&0& \color{red}{⊚} \\ 0&0&0&0 \end{bmatrix} , \quad \begin{bmatrix} \color{red}{♦} & \ast & \ast & \ast & \ast & \ast \\ 0&0&0 & \color{red}{♦} & \ast & \ast \\ 0&0&0&0& \color{red}{♦} & \ast \\ 0&0&0&0&0& \color{red}{⊚} \end{bmatrix} .$
where denotes the pivot's position (the entry cannot be zero), * denotes arbitrary element that could be zero or not, and denotes lonely nonzero entry that looks as a pivot but it indicates that the corresponding system has no solution. From theoretical point of view, a pivot could be in the last column, but when dealing with augmented matrices corresponding to the linear system of equations, we avoid this. For consistent systems, pivots cannot be in lonely position at the end of a row; otherwise, the system has no solution. ■
A pivot position in a matrix A is a location in A that corresponds to a leading term in the row echelon form of A. A pivot column is a column of A that contains a pivot position.

Any nonzero matrix may be row reduced (that is, transformed by elementary row operations) into more than one matrix in echelon form, using different sequences of row operations. However, the leading entries are always in the same positions in any echelon form obtained from a given matrix.

Theorem: A linear system has a solution if and only if the rightmost column of the associated augmented matrix is not a pivot column---that is, if and only if its echelon form does not contain a row of the form
[ 0 0 ... 0 ] with nonzero.
If a linear system is consistent, then the solution set contains either a unique solution (with no free variables) or infinitely many solutions, when there is at least one free variable.

Row echelon form states that the Gaussian elimination method has been specifically applied to the rows of the matrix. It uses only those operations that preserve the solution set of the system, known as elementary row operations:

1. Addition of a multiple of one equation to another.
Symbolically: (equation j) $$\mapsto$$ (equation j) + k (equation i).
2. Multiplication of an equation by a nonzero constant k.
Symbolically: (equation j) $$\mapsto$$ k (equation j).
3. Interchange of two equations>
Symbolically: (equation j) $$\Longleftrightarrow$$ (equation i).

Row echelon and Reduced row echelon forms are the resulting matrices of the Gaussian elimination method. By applying Gaussian elimination to any matrix, one can easily deduce the following information about the given matrix:

• the rank of the matrix;
• the determinant of the matrix;
• the inverse (invertible square matrices only);
• the kernel vector (also known as �null vector�).
The main purpose of writing a matrix in its row echelon form and/or reduced row echelon, is to make it easier to examine the matrix and to carry out further calculations, especially when solving a system of algebraic equations.

When forward elimination procedure is applied to a system of algebraic equations Ax=b, the first step is create an augmented matrix, which is obtained by appending the columns vector b from right to the matrix of the system: $${\bf B} = \left[ {\bf A} \, \vert \, {\bf b} \right] .$$ The next step is to use elementary row operations to reduce the augmented matrix B to the new augmented matrix $${\bf C} = \left[ {\bf U} \, \vert \, {\bf c} \right] ,$$ where U is upper triangular matrix. This means that the new system Ux=c is easy to solve.

The actual use of the term augmented matrix was evolved by the American mathematician Maxime B�cher (1867--1918) in his book Introduction to Higher Algebra, published in 1907. B�cher was an outstanding expositor of mathematics whose elementary textbooks were greatly appreciated by students. His achievements were documented by William F. Osgood (1919).

Variables in a linear system of equations that corresponds to pivot positions in the augmented matrix for the given system are called the leading variables. the remaining variables are called free variables.
Example: We consider the following augmented matrix:
${\bf A} = \begin{bmatrix} \color{red}{-1}&3&-2&5&4&0&7 \\ 0&0& \color{red}{3}&6&-1&0&2 \\ 0&0&0&0&0&\color{red}{7}&49 \\ 0&0&0&0&0&0&0 \end{bmatrix} .$
This matrix corresponds to four equations in six unknown variables. Since pivots are located in columns 1, 3, and 6, the leading variables for this system of equations are x1, x3, and x6. The other three variables x2, x4, and x5 are free variables. ■

Theorem: If a homogeneous linear system has n unknowns, and if the row echelon form of its augmented matrix has r nonzero rows, then the system has n-r free variables.

Row operations can be applied to any matrix, not merely to one that arises as the augmented matrix of a linear system. Two matrices are called row equivalent if there is a sequence of elementary row operations that transforms one matrix into the other. It is important to note that row operations are reversible. if two rows are interchanged, they can be returned to their original positions by another interchange. If a row scaled is scaled by a nonzero constant k, then multiplying the new row by 1/k produces the original row. Finally, consider a replacement operation involving two rows---say 1 and 3---and suppose that k times row 1 is added to row 3 to obtain a new row 3. To come back, just add -k times row 1 to new row 3 and you will get the original row 3.

Theorem: Let $${\bf B} = \left[ {\bf A} \, \vert \, {\bf b} \right]$$ be the augmented matrix corresponding to the linear equation Ax=b, and suppose B is row equivalent (using a sequence of elementary row operations) to the new augmented matrix $${\bf C} = \left[ {\bf U} \, \vert \, {\bf c} \right] ,$$ which corresponds to the linear system Ux=c. Then the two linear systems have precisely the same solution set.

Naive Gaussian elimination algorithm for Ax=b with a square matrix A (pseudocode)
 
for i=1 to n-1
for j=i+1 to n
m=a(j,i)/a(i,i)
for k=i+1 to n
a(j,k) = a(j,k) -m&a(i,k)
endfor
b(j) = b(j) - m*b(i)
endfor
endfor

The outmost loop (the i loop) ranges over the columns of the matrix; the last column is skipped because we do not need to perform any eliminations there because there are no elements below the diagonal (if we were doing elimination on a nonsquare matrix with more rows than columns, we would have to include the last column in this loop).

The middle loop (the j loop) ranges down the i-th column, below the diagonal (hence j ranges only from i+1 to n---the dimension of the matrix A). We first compute the multiplier m, for each row. This is the constant by which we multiply the i-th row in order to eliminate the aji element. Note that we overwrite the previous values with the new ones, and we do not actually carry out computation that makes aji zero. Also this loop is where the right-side vector is modified to reflect the elimination step.

The innermost loop (the k loop) ranges across the j-th row, starting after the i-th column, modifying each element appropriately to reflect the elimination of aji.

Finally, we must be aware that the algorithm does not actually create the zeroes in the lower triangular half of B; this would be wasteful of computer time since we don't need t have these zeroes in place for the algorithm to work. The algorithm works because we utilize only the upper triangular part of the matrix from this point forward, so the lower triangular elements need never be referenced.

Backward solution algorithm for A (pseudocode)

 
x(n) = b(n)/a*n,n)
for i=n-1 to 1
sum = 0
for j=i+1 to n
sum = sum + a(i,j)*x(j)
endfor
x(i) = (b(i) - sum)/a(i,i)
endfor

This algorithm simply matches backward up the diagonal, computing each xi in turn. Finally, we are computing
$x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j=i+1}^n a_{ij} x_j \right) ,$
which is what is necessary to solve a triangular system. The j loop is simply accumulating the summation term in this formula. The algorithm stops if one of the diagonal terms is zero because we cannot divide by it. This case requires a special attention that yields interchange of row.
Example: Consider system of algebraic equations
$\begin{split} 3\,x-3\,y-6\,z &=3, \\ 2\,x + 4\,y + z &= 1 ,\\ -x + y + 2\,z &= -1. \end{split}$
The idea of Gaussian elimination is to replace the above system by another system with the same solutions, it is just easier to solve. First, we build the augmented matrix corresponding to the given system of equations:
$\left[ {\bf A}\,{\bf b}\right] = \begin{bmatrix} 3&-3&-6&3 \\ 2&4& 1 & 1 \\ -1 &1&2&-1 \end{bmatrix} .$
Multiplying the first row by -2/3 and adding to the second one, we obtain an equivalent augmented matrix:
$\left[ {\bf A}\,{\bf b}\right] \sim \begin{bmatrix} 3&-3&-6&3 \\ 0&6&4&-1 \\ -1 &1&2&-1 \end{bmatrix} .$
Next, we multiply the first row by 1/3 and add to the third row:
$\left[ {\bf A}\,{\bf b}\right] \sim \begin{bmatrix} 3&-3&-6&3 \\ 0&6&4&-1 \\ 0 &0&0&0 \end{bmatrix} .$
Therefore, we obtain an equivalent augmented matrix in row echelon form. Since it contains one row of all zeroes, the given system has infinite many solutions that we obtain by solving the second equation:
$6\,y + 4\,z = -1 \qquad \Longrightarrow \qquad z = - \frac{3}{2}\, y - \frac{1}{4} .$
Using this expression, we get from the first equation
$3\,x-3\,y-6\,z = 3x + 6\,y + \frac{3}{2} = 3 \qquad \Longrightarrow \qquad x = -2\,y + \frac{3}{2} .$
So y is a free variable, and we obtain the general solution:
$x = -2\,y + \frac{3}{2} , \quad z = - \frac{3}{2}\, y - \frac{1}{4} .$

Example: Consider system of algebraic equations
$\begin{split} 3\,x-3\,y-6\,z &=3, \\ 2\,x + 4\,y + z &= 1 ,\\ -x + y + 2\,z &= 1. \end{split}$
We apply the same procedure as in the previous example: forward elimination. The procedure to be used expresses some of unknowns in terms of others by eliminating certain unknowns from all the equations except one. To begin, we eliminate x from every equation except the first one by adding -2/3 times the first equation to the second and 1/3 times the first equation to the third. The result is the following new system:
$\begin{bmatrix} 1&-1&-2&1 \\ 2&4&1&1 \\ -1&1&2&1 \end{bmatrix} \qquad \Longrightarrow \qquad \begin{bmatrix} 1&-1&-2&1 \\ 0&6&5&3 \\ 0&0&0&2 \end{bmatrix} .$
Since the pivot is situated in the last column of the augmented matrix, the given system of equations has no solution because it is impossible to satisfy the equation 0 = 2.

Example: When using the Gaussian elimination technique, you can at any time exchange rows, meaning that you can switch any two rows an unlimited number of times. This is very helpful if your matrix contains a 0 in the (1,1) position.

For example, consider the system of equations

$\begin{split} 2\,x_2 + x_3 &=-8, \\ x_1 -2\,x_2 -3\,x_3 &= 0 ,\\ -x_1 + x_2 + 2\,x_3 &= 3; \end{split} \qquad \begin{bmatrix} 0&2&1&-8 \\ 1&-2&-3&0 \\ -1&1&2&3 \end{bmatrix} .$
Because we cannot have 0 in the (1,1) position, we swap row 1 and row 2:
$\begin{bmatrix} 1&-2&-3&0 \\ 0&2&1&-8 \\ -1&1&2&3 \end{bmatrix} .$
Adding row 1 to row 3 gives
$\begin{array}{c} 1 \cdot [\mbox{equation }1] \\ + [\mbox{equation }3] \\ \hline [\mbox{new equation }3] \end{array} \qquad \begin{array}{ccccc} x_1 &-2\,x_2& -3\, x_3& =&0 \\ -x_1&+ x_2 & +2\,x_3& =&3 \\ \hline &-x_2& -x_3 & =& 3 \end{array}$
Now we swap row 2 and row 3:
$\begin{bmatrix} 1&-2&-3&0 \\ 0&-1&-1&3 \\ 0&2&1&-8 \end{bmatrix} .$
Adding 2 times row 2 to row 3 gives us
$\begin{bmatrix} 1&-2&-3&0 \\ 0&-1&-1&3 \\ 0&0&-1&-2 \end{bmatrix} .$
Now add -1 times row 3 to row 2 and add -3 times row 3 to row 1 which gives us
$\begin{bmatrix} 1&-2&0&-4 \\ 0&-1&0&5 \\ 0&0&-1&-2 \end{bmatrix} .$
Add -2 times row 2 to row 1
$\begin{bmatrix} 1&0&0&-4 \\ 0&-1&0&5 \\ 0&0&-1&-2 \end{bmatrix} .$
Finally, we multiply rows 2 and 3 by -1 to obtain Gauss-Jordan form:
$\begin{bmatrix} 1&0&0&-4 \\ 0&1&0&-5 \\ 0&0&1&2 \end{bmatrix} .$
This gives x1 = -4, x2 = -5, and x3 = 2. ■

1. For each set of matrices, determine which elementary operations can transform matrix A to matrix B.
1. ${\bf A} = \begin{bmatrix} 1 & 5 & -3 & 2\\ 2 & 2 & 3 & 5\\ 0 & 2 & 4 & -3\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 1 & 5 & -3 & 2\\ 0 & -8 & 9 & 1\\ 0 & 2 & 4 & -3\\ \end{bmatrix} ;$
2. ${\bf A} = \begin{bmatrix} 2 & 3 & 7 & 1\\ 0 & 4 & 3 & 5\\ 1 & 2 & 7 & 4\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 2 & 3 & 7 & 1\\ 1 & 2 & 7 & 4\\ 0 & 4 & 3 & 5\\ \end{bmatrix} ;$
3. ${\bf A} = \begin{bmatrix} 3 & 2 & 1 & -1\\ 4 & 6 & 5 & 2\\ 1 & 7 & -2 & 3\\ \end{bmatrix} \qquad {\bf B} = \begin{bmatrix} 0 & -19 & 7 & -10\\ 4 & 6 & 5 & 2\\ 1 & 7 & -2 & 3\\ \end{bmatrix} .$
2. For the following problems, use row operations to solve the system of equations.
1. \begin{equation*} x_{1}+2\,x_{2}=-3 \end{equation*} \begin{equation*} 2\,x_{1}+3\,x_{2}=4 \end{equation*}
2. \begin{equation*} 3\,x_{1}+x_{2}=6 \end{equation*} \begin{equation*} 3\,x_{1}+3\,x_{2}=2 \end{equation*}
3. In the following problems, solve the systems of linear equations given their augmented matrices.
1. \begin{equation*} \begin{bmatrix} 1 & -1 & 2 & -1 & | & 5\\ 0 & 5 & -1 & 3 & | & 4\\ 0 & 0 & 3 & -1 & | & 5\\ 0 & 0 & 0 & 2 & | & 8\\ \end{bmatrix} \end{equation*}
2. \begin{equation*} \begin{bmatrix} 1 & -2 & 2 & | & 4\\ 0 & 3 & -2 & | & 1\\ 0 & 0 & 0 & | & 3\\ \end{bmatrix} \end{equation*}
4. For the following problems, solve the system of equations.
1. \begin{equation*} 2\,x_{2}+x_{3}=3 \end{equation*} \begin{equation*} 2\,x_{1}+x_{2} -\,2x_{3}=-3 \end{equation*} \begin{equation*} 4\,x_{1}+4\,x_{2}-3\,x_{3}=-3 \end{equation*}
2. \begin{eqnarray*} 3\,x_{1}-2\,x_{2}-1\,x_{3} &=-1 , \\ 2\,x_{1}-1\,x_{2}+2\,x_{3} &=4 , \\ x_{2}-2\,x_{3} &=-6. \end{eqnarray*}
5. For the following problems, find the value or values of f, if any, which make the matrix an augmented matrix of a consistent linear system.
1. $\begin{bmatrix} 2 & f & 3\\ -4 & -8 & -6 \end{bmatrix} ;$
2. $\begin{bmatrix} 1 & f & -3\\ 2 & 4 & -8 \end{bmatrix} ;$
3. $\begin{bmatrix} 2 & -1 & f\\ -4 & 2 & 6\\ \end{bmatrix} .$
6. For the following problems, use row operations to solve the system of equations.
1. \begin{equation*} x_{1}+x_{2}-x_{3}=9 \end{equation*} \begin{equation*} x_{2}+3x_{3}=3 \end{equation*} \begin{equation*} -x_{1}-2x_{3}=2 \end{equation*}
2. \begin{equation*} 4x_{1}+8x_{2}-4x_{3}=4 \end{equation*} \begin{equation*} 3x_{1}+8x_{2}+5_{3}=-11 \end{equation*} \begin{equation*} -2x_{1}+x_{2}+12x_{3}=-17 \end{equation*}
3. \begin{equation*} x_{1}+2x_{2}+3x_{3}=4 \end{equation*} \begin{equation*} 2x_{1}+3x_{2}+3x_{3}=1 \end{equation*} \begin{equation*} 4x_{1}+2x_{2}+3x_{3}=1 \end{equation*}
7. Determine whether the following augmented matrices represent consistent systems of equations and how many solutions exist.
1. A= \begin{bmatrix} 1 & 4 & -2 & | & 8\\ 5 & 7 & -5 & | & 6\\ -3 & 2 & -6 & | & 6\\ \end{bmatrix}
2. A= \begin{bmatrix} 1 & 1 & 1 & | & 3\\ 2 & 0 & -1 & | & 1\\ 0 & 1 & -1 & | & 0\\ \end{bmatrix}
3. A= \begin{bmatrix} 1 & -1 & -1 & | & 1\\ -1 & 2 & 2 & | & -4\\ 3 & 2 & 1 & | & -12\\ \end{bmatrix}
4. A= $\begin{bmatrix} 3 & 1 & -3 & | & 1\\ -1 & 2 & 2 & | & -1\\ 2 & 3 & -1 & | & 1\\ \end{bmatrix}$

1. Casady, Chelsea, Gaussian Elimination and its History, Liberty University
2. Grcar, Joseph F., Mathematicians of Gaussian Elimination, Notes of American Mathematical Society, 58, No 6, 2011, 782--792.
3. William F. Osgood, (1919). "The life and services of Maxime B�cher". Bull. Amer. Math. Soc. 25 (8): 337�350. doi:10.1090/s0002-9904-1919-03198-x.
4. Gaussian elimination for dummies
5. Examples of Gaussian Elimination, Dartmouth College.
6. Higham, Nicholas, Gaussian Elimination, Manchester Institute for Mathematical Sciences, School of Mathematics, The University of Manchester, 2011.
7. Systems of Linear Equations
8. Elementary Row Operations
9. Elementary Row Operations for Matrices
10. Automatic Row Reduction Input Form, A program that performs Gaussian elimination similarly to a human working on paper.
11. Reduce Row Echelon Form on line by Kardi Teknomo.
12. An online tool for performing elementary operations with rows and columns,
13. Gauss--Jordan Elimination Calculator.