🧮 Gaussian Elimination

What Gaussian elimination (row reduction) is and what it is used for
The three elementary row operations
Row echelon form vs. reduced row echelon form, and Gauss–Jordan elimination
Forward elimination, back substitution, and the matrix-decomposition view
A worked example solving a 3-variable linear system
Applications: determinants, matrix inverses, ranks and bases
Computational efficiency, the Bareiss algorithm, and numerical instability
Generalizations and pseudocode with partial pivoting

💡 Overview

Gaussian elimination, also known as row reduction, is an algorithm for solving systems of linear equations. It consists of a sequence of row-wise operations performed on the corresponding matrix of coefficients.

It can also be used to compute:
- the rank of a matrix
- the determinant of a square matrix
- the inverse of an invertible matrix
The method is named after Carl Friedrich Gauss (1777–1855).
Goal of row reduction: use a sequence of elementary row operations to modify the matrix until the lower left-hand corner of the matrix is filled with zeros, as much as possible.

🔧 Elementary Row Operations

There are three types of elementary row operations:

Type	Operation
1	Swapping (interchanging) two rows
2	Multiplying a row by a nonzero number (non-zero scalar)
3	Adding a multiple of one row to another row

If the matrix is associated to a system of linear equations, these operations do not change the solution set.
Therefore, if the goal is to solve a system of linear equations, using these row operations can make the problem easier.

📐 Echelon Forms

Leading entries (pivots)

For each row that does not consist of only zeros, the leftmost nonzero entry is called the leading entry (or pivot) of that row.
If two leading entries are in the same column, a type 3 row operation can be used to make one of those entries zero.
Using row swaps, one can always order the rows so that for every non-zero row, the leading entry is to the right of the leading entry of the row above.

Row echelon form

A matrix with the ordering above is said to be in row echelon form.

The lower left part of the matrix contains only zeros.
All of the zero rows are below the non-zero rows.
The word "echelon" is used because one can roughly think of the rows being ranked by their size, with the largest at the top and the smallest at the bottom.

Example of a matrix in row echelon form:

$\begin{bmatrix}0&2&1&-1\\0&0&3&1\\0&0&0&0\end{bmatrix}.$

It is in echelon form because the zero row is at the bottom and the leading entry of the second row (in the third column) is to the right of the leading entry of the first row (in the second column). Its leading entries are 2 and 3.

Reduced row echelon form

A matrix is in reduced row echelon form if, furthermore:

Each nonzero row is above every zero row.
All of the leading entries are equal to 1 (achievable with a type 2 operation).
In every column containing a leading entry, all of the other entries in that column are zero (achievable with type 3 operations).
The leading 1 in each nonzero row is to the right of the leading 1 in the previous row.

Using the three row operations, a matrix can always be transformed into reduced row echelon form. This final form is unique: it is independent of the sequence of row operations used.

Example sequence of row operations

Two elementary operations on different rows are done at the first and third steps. The third and fourth matrices are in row echelon form, and the final matrix is the unique reduced row echelon form:

$\begin{bmatrix}1&3&1&9\\1&1&-1&1\\3&11&5&35\end{bmatrix}\to \begin{bmatrix}1&3&1&9\\0&-2&-2&-8\\0&2&2&8\end{bmatrix}\to \begin{bmatrix}1&3&1&9\\0&-2&-2&-8\\0&0&0&0\end{bmatrix}\to \begin{bmatrix}1&0&-2&-3\\0&1&1&4\\0&0&0&0\end{bmatrix}$

Gaussian vs. Gauss–Jordan elimination

Term	Meaning
Gaussian elimination	The process until it has reached its upper triangular, (unreduced) row echelon form
Gauss–Jordan elimination	Using row operations to convert a matrix all the way into reduced row echelon form

For computational reasons, when solving systems of linear equations, it is sometimes preferable to stop row operations before the matrix is completely reduced.

⚙️ The Algorithm: Two Parts

The process of row reduction uses elementary row operations and can be divided into two parts.

Forward elimination
- Reduces a given system to row echelon form.
- From this form one can tell whether there are no solutions, a unique solution, or infinitely many solutions.
Back substitution
- Continues to use row operations until the solution is found.
- In other words, it puts the matrix into reduced row echelon form.

Matrix-decomposition view

Another point of view, very useful for analyzing the algorithm, is that row reduction produces a matrix decomposition of the original matrix.

The elementary row operations may be viewed as multiplication on the left of the original matrix by elementary matrices.
A sequence of elementary operations that reduces a single row may be viewed as multiplication by a Frobenius matrix.
Then:
- the first part of the algorithm computes an LU decomposition
- the second part writes the original matrix as the product of a uniquely determined invertible matrix and a uniquely determined reduced row echelon matrix

📝 Worked Example: Solving a System

Goal: find and describe the set of solutions of

$\begin{aligned}2x+y-z&=8 &&(L_1)\\-3x-y+2z&=-11 &&(L_2)\\-2x+y+2z&=-3 &&(L_3)\end{aligned}$

In practice, one does not usually deal with systems in terms of equations, but makes use of the augmented matrix, which is more suitable for computer manipulations.
Procedure summary:
1. Eliminate $x$ from all equations below $L_1$ .
2. Eliminate $y$ from all equations below $L_2$ .
3. This puts the system into triangular form.
4. Using back-substitution, solve for each unknown.

First steps

$x$ is eliminated from $L_2$ by adding $\tfrac{3}{2}L_1$ to $L_2$ .
Next, $x$ is eliminated from $L_3$ by adding $L_1$ to $L_3$ .
These row operations are labelled as:

$\begin{aligned}L_{2}+{\tfrac {3}{2}}L_{1}&\to L_{2},\\L_{3}+L_{1}&\to L_{3}.\end{aligned}$

Finishing

Once $y$ is also eliminated from the third row, the result is a system in triangular form, and the first part of the algorithm is complete.
From a computational point of view, it is faster to solve the variables in reverse order (back-substitution).
The solution is $z = -1$ , $y = 3$ , $x = 2$ . There is a unique solution to the original system.
Instead of stopping at row echelon form, one can continue to reduced row echelon form. This is sometimes referred to as Gauss–Jordan elimination, to distinguish it from stopping after reaching echelon form.

🚀 Applications

Historically, the first application of row reduction is solving systems of linear equations. Other important applications follow.

Computing determinants

How elementary row operations change the determinant:

Row operation	Effect on the determinant
Swapping two rows	Multiplies the determinant by $-1$
Multiplying a row by a nonzero scalar	Multiplies the determinant by the same scalar
Adding to one row a scalar multiple of another	Does not change the determinant

If Gaussian elimination applied to a square matrix $A$ produces a row echelon matrix $B$ , let $d$ be the product of the scalars by which the determinant has been multiplied. Then the determinant of $A$ is the quotient by $d$ of the product of the diagonal elements of $B$ :

$\det(A)=\frac{\prod \operatorname{diag}(B)}{d}.$

Cost comparison for an $n \times n$ matrix:

Method	Operations
Gaussian elimination	only $O(n^3)$ arithmetic operations
Leibniz formula	$(n\,n!)$ operations (number of summands times the number of multiplications in each summand)
Recursive Laplace expansion	$O(n\,2^n)$ operations if the sub-determinants are memorized and computed only once

Even on the fastest computers, the Leibniz and Laplace methods are impractical or almost impracticable for $n$ above 20.

Finding the inverse of a matrix

A variant called Gauss–Jordan elimination can be used to find the inverse of a matrix, if it exists.

For an $n \times n$ square matrix $A$ , augment the $n \times n$ identity matrix to the right of $A$ , forming an $n \times 2n$ block matrix $[A \mid I]$ .
Apply elementary row operations to find the reduced echelon form of this matrix.
$A$ is invertible if and only if it can be reduced to the identity matrix $I$ . In that case the right block of the final matrix is $A^{-1}$ .
If the algorithm cannot reduce the left block to $I$ , then $A$ is not invertible.

Example. Take

$A=\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix}.$

Row-reduce the $3 \times 6$ augmented matrix

$[A|I]=\left[\begin{array}{ccc|ccc}2&-1&0&1&0&0\\-1&2&-1&0&1&0\\0&-1&2&0&0&1\end{array}\right].$

Its reduced row echelon form is

$\left[\begin{array}{rrr|rrr}1&0&0&\frac{3}{4}&\frac{1}{2}&\frac{1}{4}\\0&1&0&\frac{1}{2}&1&\frac{1}{2}\\0&0&1&\frac{1}{4}&\frac{1}{2}&\frac{3}{4}\end{array}\right]=:[I|B].$

Why it works:

Each row operation is a left product by an elementary matrix.
On the right, the product of these elementary matrices is $B$ (since $B = BI$ ).
On the left, multiplying that product into $A$ yields the identity, so $BA = I$ .
Therefore $B = A^{-1}$ . The procedure works for square matrices of any size.

Computing ranks and bases

Gaussian elimination can be applied to any $m \times n$ matrix $A$ . For example, some $6 \times 9$ matrices can be transformed to a row echelon form like

$T=\begin{bmatrix}a&*&*&*&*&*&*&*&*\\0&0&b&*&*&*&*&*&*\\0&0&0&c&*&*&*&*&*\\0&0&0&0&0&0&d&*&*\\0&0&0&0&0&0&0&0&e\\0&0&0&0&0&0&0&0&0\end{bmatrix},$

where the stars are arbitrary entries and $a, b, c, d, e$ are nonzero entries. This echelon matrix $T$ contains a wealth of information about $A$ :

The rank of $A$ is 5, since there are 5 nonzero rows in $T$ .
The vector space spanned by the columns of $A$ has a basis consisting of its columns 1, 3, 4, 7 and 9 (the columns with $a, b, c, d, e$ in $T$ ).
The stars show how the other columns of $A$ can be written as linear combinations of the basis columns.
All of this also applies to the reduced row echelon form, which is a particular row echelon format.

⏱️ Computational Efficiency

The number of arithmetic operations is one way of measuring efficiency. To solve a system of $n$ equations in $n$ unknowns by row operations to echelon form and then solving in reverse order requires:

Operation	Count
Divisions	$n(n+1)/2$
Multiplications	$(2n^3+3n^2-5n)/6$
Subtractions	$(2n^3+3n^2-5n)/6$

Total: approximately $2n^3/3$ operations.
Thus the arithmetic complexity (time complexity, where each arithmetic operation takes one unit of time, independent of input size) is $O(n^3)$ .

When this is a good measure of time:

When the time per arithmetic operation is approximately constant, as with floating-point coefficients or coefficients in a finite field.
If coefficients are exactly represented integers or rational numbers, intermediate entries can grow exponentially large, so the bit complexity is exponential.

Scale:

Gaussian elimination and its variants can be used on computers for systems with thousands of equations and unknowns.
The cost becomes prohibitive for systems with millions of equations; these large systems are generally solved using iterative methods.
Specific methods exist for systems whose coefficients follow a regular pattern.

Bareiss algorithm

The first strongly-polynomial time algorithm for Gaussian elimination was published by Jack Edmonds in 1967.
Independently, and almost simultaneously, Erwin Bareiss discovered another algorithm, based on a remark that applies to a division-free variant of Gaussian elimination.
It keeps the same arithmetic complexity $O(n^3)$ but has a bit complexity of $O(n^5)$ , and is therefore strongly-polynomial.

How it differs from standard elimination:

Standard: subtract from each row $R_i$ below the pivot row $R_k$ a multiple of $R_k$ by $r_{i,k}/r_{k,k}$ , where $r_{i,k}$ and $r_{k,k}$ are the entries in the pivot column of $R_i$ and $R_k$ .
Bareiss: replace $R_i$ with

$\frac{r_{k,k}R_{i}-r_{i,k}R_{k}}{r_{k-1,k-1}}.$

This produces a row echelon form with the same zero entries as standard Gaussian elimination.

Key facts:

Each matrix entry generated by this variant is the determinant of a submatrix of the original matrix.
If one starts with integer entries, the divisions are exact and all intermediate and final entries are integers.
Hadamard's inequality bounds the absolute values of the intermediate and final entries, giving a bit complexity of $\tilde{O}(n^5)$ (soft O notation).
Since an upper bound on the size of final entries is known, $\tilde{O}(n^4)$ can be obtained with modular computation followed by Chinese remaindering or Hensel lifting.

As a corollary, these problems can be solved in strongly polynomial time with the same bit complexity:

Testing whether $m$ given rational vectors are linearly independent
Computing the determinant of a rational matrix
Computing a solution of a rational equation system $Ax = b$
Computing the inverse matrix of a nonsingular rational matrix
Computing the rank of a rational matrix

Numeric instability

One possible problem is numerical instability, caused by the possibility of dividing by very small numbers.
If the leading coefficient of a row is very close to zero, row reduction requires dividing by it, so any error in that number is amplified.
Gaussian elimination is numerically stable for diagonally dominant or positive-definite matrices.
For general matrices it is usually considered stable when using partial pivoting, even though there are examples of stable matrices for which it is unstable.

🌐 Generalizations

Gaussian elimination can be performed over any field, not just the real numbers.
Buchberger's algorithm is a generalization to systems of polynomial equations. It depends heavily on the notion of a monomial order. The choice of an ordering on the variables is already implicit in Gaussian elimination, as the choice to work from left to right when selecting pivot positions.
Computing the rank of a tensor of order greater than 2 is NP-hard. Therefore, if $P \neq NP$ , there cannot be a polynomial time analog of Gaussian elimination for higher-order tensors (matrices are array representations of order-2 tensors).

💻 Pseudocode

Gaussian elimination transforms a given $m \times n$ matrix $A$ into row-echelon form. Here $A[i, j]$ denotes the entry in row $i$ and column $j$ , with indices starting from 1. The transformation is done in place, so the original matrix is lost, replaced by its row-echelon form.

h := 1 /* Initialization of the pivot row */
k := 1 /* Initialization of the pivot column */

while h ≤ m and k ≤ n:
    /* Find the k-th pivot: */
    i_max := argmax (i = h ... m, abs(A[i, k]))
    if A[i_max, k] = 0:
        /* No pivot in this column, pass to next column */
        k := k + 1
    else:
        swap rows(h, i_max)
        /* Do for all rows below pivot: */
        for i = h + 1 ... m:
            f := A[i, k] / A[h, k]
            /* Fill with zeros the lower part of pivot column: */
            A[i, k] := 0
            /* Do for all remaining elements in current row: */
            for j = k + 1 ... n:
                A[i, j] := A[i, j] - A[h, j] * f
        /* Increase pivot row and column */
        h := h + 1
        k := k + 1

This algorithm differs slightly from the one discussed earlier by choosing a pivot with the largest absolute value (partial pivoting).
Partial pivoting may be required if, at the pivot place, the entry is zero.
In any case, choosing the largest possible absolute value of the pivot improves numerical stability when floating point is used.
Upon completion, the matrix is in row echelon form and the corresponding system may be solved by back substitution.

Gaussian Elimination: Linear Algebra Study Notes