Gaussian Elimination for Solving Linear Systems

#58 · Linear Algebra · Medium

Problem

Implement Gaussian elimination with back-substitution to solve a system of linear equations Ax = b. Transform the augmented matrix to upper triangular form, then solve by back-substitution.

Solution

import numpy as np

def gaussian_elimination(A, b):
    A = np.array(A, dtype=np.float64)
    b = np.array(b, dtype=np.float64).reshape(-1, 1)
    n = A.shape[0]

    # Augmented matrix
    aug = np.hstack([A, b])

    # Forward elimination with partial pivoting
    for col in range(n):
        # Find pivot
        max_row = col
        for row in range(col + 1, n):
            if abs(aug[row, col]) > abs(aug[max_row, col]):
                max_row = row

        aug[[col, max_row]] = aug[[max_row, col]]

        if abs(aug[col, col]) < 1e-12:
            continue

        for row in range(col + 1, n):
            factor = aug[row, col] / aug[col, col]
            aug[row, col:] -= factor * aug[col, col:]

    # Back-substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = aug[i, n]
        for j in range(i + 1, n):
            x[i] -= aug[i, j] * x[j]
        x[i] /= aug[i, i]

    return x.tolist()

Explanation

Form the augmented matrix [A | b].
Forward elimination: For each column, find the row with the largest absolute value (partial pivoting) and swap it into the pivot position. Eliminate all entries below the pivot.
Back-substitution: Starting from the last row, solve for each variable by substituting already-known values.
Partial pivoting improves numerical stability by avoiding division by very small numbers.

Complexity

Time: O(n^3) for the elimination plus O(n^2) for back-substitution
Space: O(n^2) for the augmented matrix

← #57 #59 →