Navigation

Solving linear systems with Gaussian elimination

When it comes to machine learning (ML), especially in the early stages, understanding the core mathematical operations is critical. Linear algebra plays a foundational role in many ML algorithms, and one essential technique is Gaussian elimination.

In this blog post, we’ll explore Gaussian elimination, an algorithm for solving systems of linear equations, and how it’s implemented programmatically using Python and NumPy. This understanding will give you a more profound insight into the building blocks of machine learning, especially for algorithms involving linear regression, principal component analysis (PCA), and more.

Why is Gaussian Elimination Important in ML?

Gaussian elimination is widely used to:

  • Linear Regression: Solving normal equations to find optimal coefficients.
  • Principal Component Analysis (PCA): Calculating eigenvectors and eigenvalues.
  • Neural Networks: Assisting in weight adjustments during training.
  • Control Systems: Solving state-space representations and transfer functions.
  • Optimization Problems: Used in linear and quadratic programming.
  • Computer Graphics: Handling transformations and rendering equations.
  • Data Imputation: Filling in missing values in datasets using linear equations.

What is Gaussian Elimination?

Gaussian elimination transforms a given matrix into row echelon form (upper triangular matrix) through a series of elementary row operations. The system is then solved using back substitution. It consists of two main phases:

  1. Forward elimination: Convert the matrix into row echelon form.
  2. Back substitution: Solve for the variables once the matrix is in row echelon form.

Let’s now move on to the step-by-step implementation of this algorithm.

Python Implementation of Gaussian Elimination

We will implement Gaussian elimination using Python and NumPy to handle matrices efficiently.

Step 1: Augmenting the Matrix

An augmented matrix is a matrix that combines the matrix of coefficients A and the vector of constants B from a system of linear equations.

How to Form the Augmented Matrix

Given a system of linear equations in matrix form:

A * X = B

where:

  • A is an n × n matrix of coefficients,
  • X is an n × 1 column vector of variables,
  • B is an n × 1 column vector of constants,

the augmented matrix M is created by appending B as an additional column to A.

This can be represented as:

M = [A ∣ B]

This results in an n × (n+1) matrix, where the first n columns correspond to the coefficients of the system and the last column contains the constants.

To solve AX=B, we first augment matrix A with vector B, forming a new matrix M.

Note: For this implementation, matrix (A) is always square, accommodating scenarios with (n) equations and (n) variables.

import numpy as np

def augmented_matrix(A, B):
    # Horizontally stack A and B to form the augmented matrix
    M = np.hstack((A, B))
    return M

Step 2: Row Operations and Swapping

We need to perform row operations, including swapping rows when required. This function will help swap rows to ensure we have a valid pivot (non-zero element) during the elimination process.

It does not change the original matrix but returns a new one.

def swap(M, index_1, index_2):
    # Swap the rows at index_1 and index_2
    M[[index_1, index_2]] = M[[index_2, index_1]]
    return M

Step 3: Get Non-zero Element for Pivoting

This function becomes essential when encountering a 0 value during row operations. It determines whether a non-zero value exists below the encountered zero, allowing for potential row swaps.

To ensure numerical stability, we often look for non-zero pivot elements in the current column, checking rows below the current one.

Our aim is to get all zeros below every pivot element. So whenever this function encounters the non-zero number below the pivot it will return the index of the row to be swapped — or -1 if the entire column below the pivot is zero, which signals a singular matrix.

def get_non_zero_element_below_zero_from_column(M, starting_row, column):
    M = M.copy()
    column_array = M[starting_row:, column]
    for i, val in enumerate(column_array):
        # To check for non-zero values, you must always use np.isclose instead of doing "val == 0".
        if not np.isclose(val, 0, atol=1e-5):
            index = i + starting_row
            return index
    return -1

Step 4: Forward Elimination (Row Echelon Form)

The Row Echelon Form (REF) is a simplified matrix form used to solve systems of linear equations. A matrix is said to be in row echelon form if it satisfies the following conditions:

Conditions for Row Echelon Form

  1. Leading Entries:
    • The leftmost nonzero entry in each row is called a leading entry.
    • Each leading entry is to the right of the leading entry in the row above.
  2. Zero Rows:
    • Any rows that contain only zeros are at the bottom of the matrix.
  3. Pivot Columns:
    • In a column with a leading entry, all entries below it are zeros.

Understanding it in plain terms: row echelon form is just a “staircase” shape — each row starts its nonzero entries one column further right than the row above it, so the bottom row involves the fewest variables (often just one) and the top row involves the most. That staircase is exactly what makes back substitution possible: you solve the simplest row first, then walk upward substituting known values into rows that have more unknowns. Forward elimination’s entire job is to reshape the matrix into that staircase without changing what the system of equations actually means.

Example of Row Echelon Form

For a system of 3 equations with 3 variables:

2x +  y -  z =  8
-3x -  y + 2z = -11
-2x +  y + 2z = -3

Original (augmented) matrix:

[[ 2,  1, -1,   8],
 [-3, -1,  2, -11],
 [-2,  1,  2,  -3]]

After forward elimination (each pivot row normalized to a leading 1, as the code below does):

[[ 1,  0.5, -0.5,  4],
 [ 0,    1,    1,  2],
 [ 0,    0,    1, -1]]

Back-substituting from this form gives the full solution: x = 2, y = 3, z = -1. Every step of that reduction — including this example — is checked against the code below; there’s no illustrative shortcut here that the code wouldn’t actually produce.

We now apply Gaussian elimination to reduce the matrix into row echelon form. The function performs the elimination by making the entries below the pivot in each column zero.

To eliminate a value below a pivot, two steps are required for each row:

  • Scale the pivot row by the inverse of the pivot, so its leading entry becomes 1:

    Row → (1 / pivot) · Row

  • Then, for each row below it, subtract the appropriate multiple of the pivot row to zero out that column:

    Row(below) → Row(below) − factor · Row(pivot)

def row_echelon_form(A, B):
    A = A.astype('float64')
    B = B.astype('float64')
    M = augmented_matrix(A, B)
    num_rows = len(A)

    for row in range(num_rows):
        pivot = M[row, row]
        if np.isclose(pivot, 0, atol=1e-5):
            # Swap with a non-zero row
            swap_row = get_non_zero_element_below_zero_from_column(M, row, row)
            if swap_row == -1:
                return 'Singular Matrix'
            M = swap(M, row, swap_row)
            pivot = M[row, row]

        M[row] = M[row] / pivot  # Normalize pivot row
        for j in range(row + 1, num_rows):
            factor = M[j, row]
            M[j] -= factor * M[row]

    return M

Step 5: Back Substitution

Steps of Back Substitution

  1. Start with the Last Equation:
    • In row echelon form, the last equation will typically involve only one variable.
    • Solve for this variable directly.
  2. Substitute the Value into the Previous Rows:
    • Use the value of the last variable to solve for the second-to-last variable, then substitute both values into the third-to-last, and so on.
  3. Continue Until All Variables are Solved.

Once we have the row echelon form, we can solve for the unknown variables using back substitution. Rather than continuing explicit row operations up to reduced row echelon form, the function below computes each variable directly: working from the bottom row up, it subtracts the already-known variables’ contributions from the constant term to isolate the next one.

def back_substitution(M):
    num_rows = M.shape[0]
    solution = np.zeros(num_rows)

    for row in reversed(range(num_rows)):
        solution[row] = M[row, -1] - np.sum(M[row, :-1] * solution)

    return solution

Step 6: Gaussian Elimination

Finally, we combine everything to perform the complete Gaussian elimination process:

def gaussian_elimination(A, B):
    M = row_echelon_form(A, B)
    if isinstance(M, str):  # Check if the matrix is singular
        return M
    solution = back_substitution(M)
    return solution

Example: Solving a Linear System

Let’s use the Gaussian elimination method to solve the following system of linear equations:

x + 2y + 3z = 1
        y   = 2
       5z   = 4

Matrices:

A = [[1, 2, 3], [0, 1, 0], [0, 0, 5]]
X = [x, y, z]
B = [1, 2, 4]
A = np.array([[1, 2, 3], [0, 1, 0], [0, 0, 5]])
B = np.array([[1], [2], [4]])

solution = gaussian_elimination(A, B)
print(f"x = {solution[0]}\ny = {solution[1]}\nz = {solution[2]}")

Output:

x = -5.4
y = 2.0
z = 0.8

You can check this by substitution: with y = 2 and z = 0.8, the first equation requires x + 2(2) + 3(0.8) = 1, so x = 1 - 4 - 2.4 = -5.4.

Conclusion

Understanding Gaussian elimination deepens your understanding of the linear algebra concepts behind many machine learning algorithms. This technique not only solves systems of equations but also introduces you to matrix manipulation, a skill critical in ML. Algorithms such as linear regression or PCA heavily rely on operations with matrices, and Gaussian elimination is a foundational method for solving linear systems efficiently.

By practising and coding Gaussian elimination yourself, you’ll be well-prepared to tackle more advanced topics in machine learning confidently!

Happy coding!