Gaussian Elimination Methods

1. Introduction

Gaussian elimination is the most fundamental and commonly used direct method for solving systems of linear equations in linear algebra. Its core idea is to transform the coefficient matrix (A) into an upper triangular form through “stepwise elimination,” and then solve for the unknowns using back substitution. To improve the algorithm’s numerical stability, two improved methods have been proposed: Gaussian elimination with partial pivoting and Gaussian elimination with full pivoting. The differences among the three methods mainly lie in the “pivot selection strategy.” This article will systematically introduce the mathematical concepts, numerical stability analysis, and pure C++ implementation code for the three methods, demonstrating the evolution logic of the algorithm step by step from theoretical derivation to program implementation.

2. Gaussian Elimination: The Most Basic Method for Solving Linear Equations

1. Mathematical Principles

Given a system of linear equations:

where A is the coefficient matrix, and b is the right-hand side vector. The goal of Gaussian elimination is to use a series of elementary row operations to transform A into an upper triangular matrix, making the system become: so that the unknowns can be easily solved through back substitution.

2. Operational Process

1. Construct the augmented matrix;
2. Starting from the first row, progressively eliminate the elements below in the same column;
3. Repeat the elimination until the matrix is in upper triangular form;
4. Finally, back substitute from the bottom up to find the unknowns.

3. Numerical Characteristics

This method is logically clear and structurally simple, but if the pivot is too small, it can lead to amplified rounding errors, affecting the precision and stability of the results. Therefore, in numerical computations, a “pivot selection” strategy is usually adopted for improvement.

3. Gaussian Elimination with Partial Pivoting: Enhancing Numerical Stability

1. Basic Idea

The partial pivoting method selects the element with the largest absolute value from the current column as the pivot before each elimination step and swaps that row with the current row. This effectively avoids small pivots, thereby reducing rounding errors and improving numerical stability.

2. Calculation Steps

1. Find the element with the largest absolute value in the current column;
2. If the maximum element is not in the current row, swap the two rows;
3. Perform the usual elimination steps;
4. Finally, solve through back substitution. This method does not change the solution of the system, it merely rearranges the order of the equations. Its core idea is a “locally optimal” choice to achieve stable results at a lower cost.

4. Gaussian Elimination with Full Pivoting: Pursuing Global Stability

1. Mathematical Concept

The full pivoting method further improves upon the partial pivoting method. It searches for the element with the largest absolute value not only in the current column but also in the entire remaining submatrix as the pivot during each elimination step. This may require simultaneous row and column exchanges, thereby minimizing rounding errors. The full pivoting method is the most numerically stable form of Gaussian elimination.

2. Implementation Points

1. Search for the element with the largest absolute value in the remaining submatrix;
2. Record the row and column exchange positions;
3. Perform the row and column exchanges;
4. Complete the elimination process;
5. After obtaining the results, restore the original variable order based on the column exchange records.

3. Numerical Significance

This “globally optimal” pivot strategy can almost eliminate the amplification effect of rounding errors, but the computational load and programming complexity increase accordingly, thus it is often used in high-precision or scientific calculations.

5. Algorithm Comparison and Application Recommendations

Method Name	Pivot Selection Method	Numerical Stability	Computational Cost	Applicable Scenarios
Gaussian Elimination	None	General	Fast	Theoretical teaching, simple problems
Partial Pivoting	Maximum value in the current column	Stable	Medium	Default method for engineering calculations
Full Pivoting	Maximum value in the submatrix	Most stable	High	High-precision scientific calculations

For general engineering calculations, it is recommended to use Gaussian elimination with partial pivoting; for extremely high precision or ill-conditioned matrices, full pivoting is the best choice.

6. Pure C++ Implementation (Complete Code)

Below is the complete C++ implementation of the three algorithms. The code is completely independent, does not rely on third-party libraries, and is based solely on user-defined <span>Matrix</span> and <span>Vector</span> classes.

7. Conclusion

This article systematically introduced the mathematical concepts and pure C++ implementations of three Gaussian elimination algorithms:

• Gaussian Elimination: Simple structure, suitable for teaching and theoretical analysis.
• Gaussian Elimination with Partial Pivoting: Enhances numerical stability, the mainstream method for engineering calculations.
• Gaussian Elimination with Full Pivoting: Most stable but with the highest computational cost, suitable for scientific and high-precision problems. The three methods reflect the evolution from “direct computation” to “numerical optimization”. At the code level, pivot selection, row and column exchanges, and back substitution form the core of the algorithm. From mathematics to programming, it is about transforming abstract logic into an operable computational process; and in numerical computation, the balance between stability and efficiency is the essence of algorithm design.

#include "Matrix.h"
#include "Vector.h"
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>

// ==========================================
// 1. Gaussian Elimination (No Pivoting)
// ==========================================
Vector<double> LinearSolver::GaussianElimination(const Matrix<double>& A, const Vector<double>& b)
{
    int n = A.getRows();
    Matrix<double> Aug(n, n + 1);

    // Construct the augmented matrix
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
            Aug[i][j] = A[i][j];
        Aug[i][n] = b[i];
    }

    // Elimination
    for (int k = 0; k < n - 1; ++k)
    {
        if (fabs(Aug[k][k]) < 1e-12)
            throw std::runtime_error("Zero pivot detected!");

        for (int i = k + 1; i < n; ++i)
        {
            double factor = Aug[i][k] / Aug[k][k];
            for (int j = k; j <= n; ++j)
                Aug[i][j] -= factor * Aug[k][j];
        }
    }

    // Back substitution
    Vector<double> x(n);
    for (int i = n - 1; i >= 0; --i)
    {
        double sum = Aug[i][n];
        for (int j = i + 1; j < n; ++j)
            sum -= Aug[i][j] * x[j];
        x[i] = sum / Aug[i][i];
    }

    return x;
}

// ==========================================
// 2. Gaussian Elimination with Partial Pivoting
// ==========================================
Vector<double> LinearSolver::GaussianEliminationWithPartialPivoting(const Matrix<double>& A, const Vector<double>& b)
{
    int n = A.getRows();
    Matrix<double> Aug(n, n + 1);

    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
            Aug[i][j] = A[i][j];
        Aug[i][n] = b[i];
    }

    for (int k = 0; k < n - 1; ++k)
    {
        // Find the maximum pivot in the current column
        int pivotRow = k;
        double maxVal = fabs(Aug[k][k]);
        for (int i = k + 1; i < n; ++i)
        {
            if (fabs(Aug[i][k]) > maxVal)
            {
                maxVal = fabs(Aug[i][k]);
                pivotRow = i;
            }
        }

        if (pivotRow != k)
            std::swap(Aug[pivotRow], Aug[k]);

        // Elimination
        for (int i = k + 1; i < n; ++i)
        {
            double factor = Aug[i][k] / Aug[k][k];
            for (int j = k; j <= n; ++j)
                Aug[i][j] -= factor * Aug[k][j];
        }
    }

    // Back substitution
    Vector<double> x(n);
    for (int i = n - 1; i >= 0; --i)
    {
        double sum = Aug[i][n];
        for (int j = i + 1; j < n; ++j)
            sum -= Aug[i][j] * x[j];
        x[i] = sum / Aug[i][i];
    }

    return x;
}

// ==========================================
// 3. Gaussian Elimination with Full Pivoting
// ==========================================
Vector<double> LinearSolver::GaussianEliminationWithFullPivoting(const Matrix<double>& A, const Vector<double>& b)
{
    int n = A.getRows();
    Matrix<double> Aug(n, n + 1);

    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
            Aug[i][j] = A[i][j];
        Aug[i][n] = b[i];
    }

    std::vector<int> rowPerm(n), colPerm(n);
    for (int i = 0; i < n; ++i)
    {
        rowPerm[i] = i;
        colPerm[i] = i;
    }

    for (int k = 0; k < n - 1; ++k)
    {
        double maxVal = 0.0;
        int pivotRow = k, pivotCol = k;

        // Find the maximum pivot in the remaining submatrix
        for (int i = k; i < n; ++i)
            for (int j = k; j < n; ++j)
                if (fabs(Aug[i][j]) > maxVal)
                {
                    maxVal = fabs(Aug[i][j]);
                    pivotRow = i;
                    pivotCol = j;
                }

        if (maxVal < 1e-12)
            throw std::runtime_error("Singular matrix detected!");

        // Row exchange
        if (pivotRow != k)
        {
            std::swap(Aug[pivotRow], Aug[k]);
            std::swap(rowPerm[pivotRow], rowPerm[k]);
        }

        // Column exchange
        if (pivotCol != k)
        {
            for (int i = 0; i < n; ++i)
                std::swap(Aug[i][pivotCol], Aug[i][k]);
            std::swap(colPerm[pivotCol], colPerm[k]);
        }

        // Elimination
        for (int i = k + 1; i < n; ++i)
        {
            double factor = Aug[i][k] / Aug[k][k];
            for (int j = k; j <= n; ++j)
                Aug[i][j] -= factor * Aug[k][j];
        }
    }

    // Back substitution
    Vector<double> x(n);
    for (int i = n - 1; i >= 0; --i)
    {
        double sum = Aug[i][n];
        for (int j = i + 1; j < n; ++j)
            sum -= Aug[i][j] * x[j];
        x[i] = sum / Aug[i][i];
    }

    // Restore column exchange order
    Vector<double> xFinal(n);
    for (int i = 0; i < n; ++i)
        xFinal[colPerm[i]] = x[i];

    return xFinal;
}