Solving Systems with Python: Discovering Gaussian Elimination in Machine Learning

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.

• Support Vector Machines (SVM): Optimizing objective functions involving linear equations.

• 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.

• Robotics: Solving kinematic and dynamic equations for motion.

• Game Development: Implementing physics simulations requiring linear 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

To solve AX=BAX = BAX=B, where A is the matrix of coefficients, X is the vector of variables, and B is the vector of constants, we first augment matrix A with vector B, forming a new matrix M.

Note: For this assignment, 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 00 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 pivot it will return the index of the row to be swapped.

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)

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.

We will consider the centre element be the pivot and the value below it as variable value. To nullify the last row (row 2), two steps are required:

• Scale 𝑅1 by the inverse of the pivot:

Row 1→1/pivot⋅Row

Resulting in the updated matrix with the pivot for row 1 set to 1.

• Next, to eliminate the value below the pivot in row 11, apply the following formula:

Row 2→Row 2−value⋅Row 1

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

Once we have the row echelon form, we can now solve for the unknown variables using back substitution.

Utilizing elementary row operations, it aims to convert every element above the pivot into zeros, ending with a matrix in reduced row echelon form. The formula employed is:

Row above→Row above−value⋅Row pivot

In this equation, value denotes the value above the pivot, which initially equals 1

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:

pythonCopy codedef 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=10+y+0z=20+0+5z=4x+2y+3z=10+y+0z=20+0+5z=4x+2y+3z0+y+0z0+0+5z​=1=2=4​

In matrix form:

A=[123010005],B=[124]A = [123010005], \quad B = [124]A=​100​210​305​​,B=​124​​

pythonCopy codeA = 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 = -0.4
y = 2.0
z = 0.8

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!