Skip to content

Gaussian Elimination

Introduction

Gaussian elimination (Gauss–Jordan elimination) is a classic algorithm for solving systems of linear equations. It holds an important position and value in contemporary mathematics and is a key component in linear algebra courses.

Besides solving linear systems, Gaussian elimination can also be used for computing determinants, finding matrix inverses, and various applications in computer science and engineering.

The Idea of Elimination and Gaussian Elimination

Definition

Elimination involves expressing the unknown of one equation in terms of another unknown and substituting 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 achieve the elimination of an unknown. Elimination is primarily used for solving systems of two linear equations in two unknowns.

Explanation

Example 1: Use elimination to solve the following system of two linear equations in two unknowns:

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

Solution: Adding the two equations eliminates \(y\):

\[ 5x = 200 \]

Solving gives:

\[ x = 40 \]

Substituting \(x = 40\) into the second equation yields:

\[ y = -60 \]

Core Principles of Elimination

The core principles of elimination are as follows:

  • Swapping two equations does not change the solution.
  • Multiplying an equation by a nonzero number \(k\) does not change the solution.
  • Adding a multiple of one equation to another does not change the solution.

The Concept of Gaussian Elimination

The German mathematician Gauss analyzed elimination and concluded:

  • In elimination, only the coefficients of the variables are involved in computation and change.
  • The variables themselves do not participate in computation and remain unchanged.
  • The positions of coefficients can be used to represent variables, allowing variables to be omitted.
  • Simplifying by omitting variables does not change the solution of the equations.

Based on these conclusions, Gauss proposed Gaussian elimination. First, the augmented matrix of the system is transformed into reduced row echelon form using elementary row operations. Then, free variables are assigned values based on linear independence criteria. Finally, the general solution of the system is expressed.

The Five-Step Gaussian Elimination Method

Explanation

After transforming the augmented matrix into reduced row echelon form, Gaussian elimination requires assigning values to free variables, which involves knowledge of linear algebra and some human judgment, making it challenging to learn. Dividing Gaussian elimination into five steps, we propose the five-step method as follows:

  1. Transform the augmented matrix into reduced row echelon form using elementary row operations.
  2. Reconstruct the system of linear equations.
  3. Solve for the first variable.
  4. Supplement the free variables.
  5. Express the general solution in column vector form.

An example is used to illustrate the operation of this algorithm.

Process

Example 2: Use the five-step Gaussian elimination method to solve the linear system:

\[ \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

The augmented matrix is the combination of the coefficient matrix \(A\) and the constant column \(b\), denoted as \((A | b)\). Transforming the augmented matrix into reduced row echelon form using elementary row operations applies the idea of Gaussian elimination: omitting variables and representing them by their coefficient positions. The vertical line in the augmented matrix separates the coefficient matrix and the constant column, representing the equals 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) \]

Transformed to 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) \]

Transformed to reduced row echelon form.

Reconstruct the Linear System

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

Reconstructing the linear system means rewriting the reduced row echelon form back into a system of linear equations. That is, reassigning the coefficients at each position to variables and replacing the vertical line with equals signs.

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 reconstructed linear system, express the first variable of each equation in terms of the other quantities. For example, the first variables \(x_1\) and \(x_3\) in the two equations.

Supplement the Free Variables

\[ \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, variables \(x_1\) and \(x_3\) were solved. This indicates that the remaining variables \(x_2\) and \(x_4\) are not constrained by the system; they are free variables and can take any value. Therefore, we need to supplement the solution obtained in step 3. The supplementation method is \(x_2 = x_2, x_4 = x_4\), which aligns with the definition of free variables and is easy to understand: since they are free and unconstrained, they can only equal themselves.

Express the General Solution in Column Vector Form

\[ \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} \]

where \(C_1\) and \(C_2\) are arbitrary constants.

Explanation

Based on step 4, express the solution as a combination of column vectors. Since \(x_2\) and \(x_4\) are free variables and can take any value, we set them as arbitrary constants \(C_1\) and \(C_2\) on the right side of the solution, thus completing the solution of the system.

Determinant Calculation

Explanation

The determinant of an \(N \times N\) square matrix can be understood as the signed volume of the geometric body spanned by 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\), \(\sigma\) is a permutation, and \(\operatorname{sgn}(\sigma)=1\) if the number of inversions in \(\sigma\) is even, otherwise \(\operatorname{sgn}(\sigma)=−1\).

Understanding determinant invariance through volume is straightforward:

  • Transposing a matrix does not change its determinant.
  • Swapping rows (or columns) of a matrix changes the sign of the determinant.
  • Adding or subtracting rows (or columns) does not change the determinant.
  • Multiplying all elements of a row (or column) by a number \(k\) scales the determinant by \(k\).

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

Note: if at any point we cannot find a nonzero element in the current column, the algorithm should stop and return 0.

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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 said to be invertible, and \(A^{-1}\) is called its inverse matrix.

Given an \(n \times n\) matrix \(A\), the method to find its inverse is as follows:

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

The proof of correctness for this method requires more linear algebra knowledge and is omitted here due to space constraints. Interested readers can consult relevant materials.

Solving XOR Systems with Gaussian Elimination

XOR systems are systems 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" (i.e., xor or ^ in C++), and all coefficients/constants (\(a_{i,j}\) and \(b_i\)) are either \(0\) or \(1\).

Since XOR is commutative and associative, we can apply Gaussian elimination step by step. Note that we should use "XOR elimination" instead of "addition/subtraction elimination," and there is no need to multiply or divide to change coefficients (since coefficients are only \(0\) and \(1\)).

Note that the augmented matrix of an XOR system is a \(01\) matrix (containing only \(0\)s and \(1\)s). Thus, we can optimize using C++'s std::bitset, reducing the time complexity to \(O(\dfrac{n^2m}{\omega})\), where \(n\) is the number of unknowns, \(m\) is the number of equations, and \(\omega\) is typically \(32\) (machine-dependent).

Reference implementation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
std::bitset<1010> matrix[2010];  // matrix[1~n]: augmented matrix, position 0 is the constant

std::vector<bool> GaussElimination(
    int n, int m)  // n: number of unknowns, m: number of equations, returns the solution
                   // (returns an empty vector for multiple/no solutions)
{
  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;
}

Practice Problems