Gaussian Elimination

Introduction

Gaussian elimination (Gauss–Jordan elimination) is a classic algorithm for solving systems of linear equations. It has an important position and value in modern mathematics and is an essential part of teaching linear algebra.

In addition to solving systems of linear equations, Gaussian elimination can also be used to compute determinants, find matrix inverses, and in other computing and engineering applications.

The elimination method and the idea behind Gaussian elimination

Definition

The elimination method expresses the unknown in one equation of a system as an algebraic expression involving another unknown and substitutes it into another equation, thereby eliminating one unknown and obtaining a solution; alternatively, multiplying one equation by a constant and adding it to another equation can also eliminate one unknown. The elimination method is mainly used to solve systems of two linear equations in two variables.

Explanation

Example 1: Solve the following system of two linear equations using the elimination method:

$$\begin{cases} 4x+y&=100 \\ x-y&=100 \end{cases}$$

Solution: Add the two equations in the system to eliminate $y$:

$$5x = 200$$

Thus:

$$x = 40$$

Substituting $x = 40$ into the second equation of the system gives:

$$y = -60$$

The core of the elimination method

The core principles of the elimination method are as follows:

• Interchanging two equations does not change the solution.

• Multiplying one equation by a nonzero number $k$ does not change the solution.

• Multiplying one equation by $k$ and adding it to another equation does not change the solution.

The idea of Gaussian elimination

The German mathematician Gauss analyzed the elimination method and reached the following conclusions:

• In the elimination method, the quantities involved in the calculations and changed are the coefficients of the variables in the equations.

• The variables themselves do not participate in the calculations and are not changed.

• The positions of the coefficients can be used to represent the variables, so the variables can be omitted.

• Omitting the variables in the calculations does not change the solution of the system.

Based on these conclusions, Gauss proposed Gaussian elimination: first transform the augmented matrix of the system into reduced row echelon form using elementary row operations; then assign values to the free unknowns according to linear independence; finally write down the general solution of the system.

The five-step method of Gaussian elimination

Explanation

After Gaussian elimination transforms the augmented matrix into reduced row echelon form, assigning values to the free unknowns requires knowledge of linear algebra and involves human judgment, which can make the method difficult to learn. Gaussian elimination can therefore be divided into five steps, giving the five-step method:

1. Transform the augmented matrix into reduced row echelon form using elementary row operations.

2. Restore the system of linear equations.

3. Solve for the first variable.

4. Add the free unknowns.

5. Express the general solution of the system.

The operation of this algorithm is illustrated further with an example.

Process

Example 2: Solve the following system of linear equations using the five-step method of Gaussian elimination:

$$\begin{cases} 2x_1+5x_3+6x_4&=9 \\ x_3+x_4&=-4 \\ 2x_3+2x_4&=-8 \end{cases}$$

Transform the augmented matrix into reduced row echelon form using elementary row operations

The augmented matrix is the new matrix formed by adjoining the coefficient matrix $A$ and the constant column $b$, namely $(A | b)$. Transforming the augmented matrix into reduced row echelon form by elementary row operations applies the idea of Gaussian elimination: variables are omitted and represented by the positions of their coefficients, while the vertical line in the augmented matrix separates the coefficient matrix from the constant column and represents the equality sign.

$$\left(\begin{matrix} 2 & 0 & 5 & 6 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 2 & 2 \end{matrix} \middle| \begin{matrix} 9 \\ -4 \\ -8 \end{matrix} \right)$$ $$\xrightarrow{r_3-2r_2} \left(\begin{matrix} 2 & 0 & 5 & 6 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 9 \\ -4 \\ 0 \end{matrix} \right)$$

Transform it into row echelon form

$$\xrightarrow{\frac{r_1}{2}} \left(\begin{matrix} 1 & 0 & 2.5 & 3 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 4.5 \\ -4 \\ 0 \end{matrix} \right)$$ $$\xrightarrow{r_1-r_2 \times 2.5} \left(\begin{matrix} 1 & 0 & 0 & 0.5 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 14.5 \\ -4 \\ 0 \end{matrix} \right)$$

Transform it into reduced row echelon form

Restore the system of linear equations

$$\begin{cases} x_1+0.5x_4 &= 14.5\\ x_3+x_4 &= -4 \\ \end{cases}$$

Explanation

Restoring the system of linear equations means rewriting the reduced row echelon form as a system of linear equations: reassign variables to the coefficients at the corresponding positions and restore the vertical line as an equals sign.

Solve for the first variable

$$\begin{cases} x_1 = -0.5x_4+14.5\notag \\ x_3 = -x_4-4\notag \end{cases}$$

Explanation

For the restored system of linear equations, express the first variable in each equation in terms of the other quantities. In this system, these are the first variables $x_1$ and $x_3$ in the two equations.

Add the free unknowns

$$\begin{cases} x_1 = -0.5x_4+14.5 \\ x_2 = x_2 \\ x_3 = -x_4-4 \\ x_4 = x_4 \end{cases}$$

Explanation

In Step 3, we solved for $x_1$ and $x_3$, showing that the remaining variables $x_2$ and $x_4$ are not constrained by the system and are free unknowns that can take arbitrary values. Therefore, we supplement the result from Step 3 by setting $x_2 = x_2,x_4 = x_4$. This way of supplementing the solution conforms to the definition of free unknowns and is easy to understand: because a free unknown is unconstrained, it can only be equal to itself.

Express the general solution of the system

$$\begin{aligned} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} &= \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} x_2+ \begin{pmatrix} -0.5 \\ 0 \\ -1 \\ 1 \end{pmatrix} x_4 + \begin{pmatrix} 14.5 \\ 0 \\ -4 \\ 0 \end{pmatrix} \\ &= \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} C_1+ \begin{pmatrix} -0.5 \\ 0 \\ -1 \\ 1 \end{pmatrix} C_2 + \begin{pmatrix} 14.5 \\ 0 \\ -4 \\ 0 \end{pmatrix} \end{aligned}$$

Here $C_1$ and $C_2$ are arbitrary constants.

Explanation

Based on Step 4, express the solution as a linear combination of column vectors. Since $x_2$ and $x_4$ are free unknowns and can take arbitrary values, set them on the right-hand side of the solution to arbitrary constants $C_1$ and $C_2$, respectively. This solves the system.

Computing determinants

Explanation

The determinant of an $N \times N$ square matrix can be understood as the oriented volume of the geometric object spanned by all its column vectors.

For example:

$$\begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix} = 1$$ $$\begin{vmatrix} 1 & 2 \\ 2 & 1 \end{vmatrix} = -3$$

The determinant has the formula

$$\operatorname{det}(A)=\sum_{\sigma \in S_{n}} \operatorname{sgn}(\sigma) \prod_{i=1}^{n} a_{i, \sigma(i)}$$

where $S_n$ is the set of all permutations of length $n$, and $\sigma$ is a permutation. If the number of inversions of $\sigma$ is even, then $\operatorname{sgn}(\sigma)=1$; otherwise, $\operatorname{sgn}(\sigma)=−1$.

The following properties provide a simple way to understand the invariance of determinants geometrically:

• Transposing a matrix leaves its determinant unchanged.

• Swapping two rows (or columns) negates the determinant.

• Adding or subtracting one row (or column) from another leaves the determinant unchanged.

• Multiplying every element of a row (or column) by $k$ scales the determinant by the same factor.

Therefore, after applying Gaussian elimination to a matrix, we obtain a diagonal matrix whose determinant is determined by the product of its diagonal elements. Its sign is determined by the number of row swaps (if this number is odd, the sign of the determinant must be reversed). Thus, Gaussian elimination can compute the determinant in $O(n^3)$ time.

Note that if, at some point, no nonzero element can be found in the current column, the algorithm should stop and return 0.

Implementation

constexpr double EPS = 1E-9;
int n;
vector<vector<double>> a(n, vector<double>(n));

double det = 1;
for (int i = 0; i < n; ++i) {
  int k = i;
  for (int j = i + 1; j < n; ++j)
    if (abs(a[j][i]) > abs(a[k][i])) k = j;
  if (abs(a[k][i]) < EPS) {
    det = 0;
    break;
  }
  swap(a[i], a[k]);
  if (i != k) det = -det;
  det *= a[i][i];
  for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
  for (int j = 0; j < n; ++j)
    if (j != i && abs(a[j][i]) > EPS)
      for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
}

cout << det;

Matrix inversion

For a square matrix $A$, if there exists a square matrix $A^{-1}$ such that $A \times A^{-1} = A^{-1} \times A = I$, then $A$ is called invertible, and $A^{-1}$ is called its inverse matrix.

Given an $n$-th order square matrix $A$, its inverse can be found as follows:

1. Construct the $n \times 2n$ matrix $(A, I_n)$.
2. Use Gaussian elimination to reduce it to reduced row echelon form $(I_n, A^{-1})$; this yields $A$'s inverse matrix $A^{-1}$. If the left half of the final reduced row echelon form is not the identity matrix $I_n$, then $A$ is not invertible.

A proof of this method's correctness requires substantial linear algebra knowledge and is omitted here for brevity. Interested readers can consult the relevant references.

Solving XOR systems of equations with Gaussian elimination

An XOR system of equations is a system of the form

$$\begin{cases} a_{1,1}x_1 \oplus a_{1,2}x_2 \oplus \cdots \oplus a_{1,n}x_n &= b_1\\ a_{2,1}x_1 \oplus a_{2,2}x_2 \oplus \cdots \oplus a_{2,n}x_n &= b_2\\ \cdots &\cdots \\ a_{m,1}x_1 \oplus a_{m,2}x_2 \oplus \cdots \oplus a_{m,n}x_n &= b_m \end{cases}$$

where $\oplus$ denotes bitwise XOR (that is, xor or ^ in C++), and all coefficients/constants in the expression (namely $a_{i,j}$ and $b_i$) are either $0$ or $1$.

Because XOR satisfies commutativity and associativity, the system can be solved by eliminating variables step by step as in Gaussian elimination. It is important to use XOR elimination rather than addition/subtraction elimination, and multiplication or division to change coefficients is unnecessary because all coefficients are $0$ or $1$.

The augmented matrix of an XOR system is a $01$ matrix (a matrix containing only $0$ and $1$), so we can optimize the implementation with C++'s std::bitset, reducing the time complexity to $O(\dfrac{n^2m}{\omega})$, where $n$ is the number of variables, $m$ is the number of equations, and $\omega$ is generally $32$ (depending on the machine).

Reference implementation:

std::bitset<1010> matrix[2010];  // matrix[1~n]: the augmented matrix; position 0 is the constant

std::vector<bool> GaussElimination(
    int n, int m)  // n is the number of unknowns, m is the number of equations, and the return value is the solution of the system
                   // (multiple solutions / no solution returns an empty vector)
{
  for (int i = 1; i <= n; i++) {
    int cur = i;
    while (cur <= m && !matrix[cur].test(i)) cur++;
    if (cur > m) return std::vector<bool>(0);
    if (cur != i) swap(matrix[cur], matrix[i]);
    for (int j = 1; j <= m; j++)
      if (i != j && matrix[j].test(i)) matrix[j] ^= matrix[i];
  }
  std::vector<bool> ans(n + 1);
  for (int i = 1; i <= n; i++) ans[i] = matrix[i].test(0);
  return ans;
}

Exercises

• Codeforces - Wizards and Bets
• luogu - SDOI2010 Alien Centipede