There's nothing quite fancy about what you're about reading, as the essence of this post is mainly to document the learning experience of translating the handwritten mathematical procedures to that of a computational procedure, and in this case, with numpy (python).
In handwritten math solutions, there's plenty room for flexibility. You rarely have to think about the performance of how your computational data is spread out in memory (as there's almost none in that sense). Your choice here determines the ease of the type of operations you intend to carry out on the data involved. As much as numpy provides couple of data storage types, you're still left with the responsibility of choosing the right one.
The above might not be obvious until you're faced with situation where you used a one dimensional array to represent a vector matrix, and sanely have to combine it with a coefficient matrix to form an augmented matrix. And directly combining both matrix with np.hstack produces undesired outcome because you forgot to reshape the vector matrix to an n-dimensional array of single element entry.
# Linear system Ax = b
self.A = A.astype(np.float64) # M x N shape
self.b = b.astype(np.float64) # N x 1 shape
# M x N+1 shape
self.augmented = np.hstack((self.A, self.b.reshape(-1, 1)))Alright, 'tis was the first omission I had to traced down whilst working on this library. It was fun figuring it out.
Moving on. I made the decision to separate the elimination implementation code out of that of the back substitution procedure. Good call that I only recognized later, mostly because I could inspect the Row Echelon Form matrix—the residual matrix remaining after elimination—and retrieve necessary information about the submitted linear system if necessary.
Whilst implementing the internals of the elimination procedure, my first attempt was a clever but hacky steps of discovering the pivot terms and rearranging the matrix to promote the rows to intended row positions to ease elimination. This was okay, but as luck (and google search) would have it, I found an easier solution already baked into numpy that I could use to my advantage.
for k in range(min(self.rows, self.n)):
if use_pivoting:
max_index = np.argmax(np.abs(M[k:, k])) + k
if k != max_index:
M[[k, max_index]] = M[[max_index, k]]The experience from this was to always check the necessary docs in hopes to find some existing API to implement a solution before attempting a hacky solution. Don't get me wrong, this isn't a message to let loose of your tinkering spirit.
Another thing that tripped me up wasn't even strictly about linear algebra, but how computers handle floating-point precision. Genuinely, I don't know much about its internals. Hence, in my handwritten notes, zero is just 0. It's clean. But when I started running test cases, I kept getting weird behavior where my implementation assume a number like 1.23e-17 as a valid pivot.
Forgive me to say, floating-point arithmetic is messy. Operations accumulate tiny bits of noise, so a value that should be zero mathematically ends up being a microscopic non-zero number. If you don't catch this, your solution would try to divide by it, and surely, everything explodes.
The only way ahead was to add some kind of "tolerance guard" (I settled on 1e-10). Basically, if a number is small enough, the code forces it to be treated as 0.0. It feels a bit careless and like cheating, but I got to find out later that was fine, at least to keep things stable.
# To account for float ops quirkness
if abs(M[i, k]) < 1e-10:
M[i, k] = 0.0Now it was time to implement the Back Substitution logic, and my first attempt left me with some bugs to hunt me during testing. lol! My initial implementation assumed that if I had N variables, I would naturally have N rows to iterate through. silly right!? I used the loop index i for both the column (variable) and the row.
As you can guess, that worked fine for square matrices. But as soon as I threw an underdetermined system at it (where you have fewer equations than variables, M < N), the code crashed hard.
When I found the mistake, I had to rethink the loop, indeed. So what? I needed to decouple the row tracker from the variable tracker. In the final code, you'll see I manually decrement a rows counter (rows -= 1) inside the loop. This lets me walk backwards through the variables (x_4, x_3, x_2...) while independently tracking which equation row corresponds to them.
Finally, there was the question of what to do when the system is dependent (there are infinitely many solutions). The easy route would have been to raise an exception or return NaN. But nah, I wanted something more tangible.
If the pivot is zero during back substitution, it means that specific variable is "free" (it can be pin to any numeric value). Here, I chose to randomly pick a value between 0, 1, or 2, and treat it as a particular solution for the variable in question. Hence, I have a solution vector that satisfies the equations, rather than just an error message.
Putting it all together, here is the final GaussianSolver class. It handles the augmentation, the pivoting, the floating-point noise, and that tricky back-substitution indexing. You can find the source code on GitHub.
import numpy as np
import random
from numpy.typing import NDArray
class GaussianSolver:
def __init__(self, A: NDArray[np.float64], b: NDArray[np.float64]) -> None:
"""
Initialize with coefficient matrix A and constant vector b.
Constructs the Augmented Matrix [A|b]
"""
self.A = A.astype(np.float64) # M x N shape
self.b = b.astype(np.float64) # N x 1 shape
if self.A.shape[0] != self.b.shape[0]:
raise ValueError("Dimensions of A and b do not match.")
self.augmented = np.hstack((self.A, self.b.reshape(-1, 1))) # M x N+1 shape
self.rows, self.cols = self.augmented.shape
# The number of unknowns which is cols - 1 (since last column accounts for 'b')
self.n = self.cols-1
def eliminate(self, use_pivoting: bool = True) -> NDArray[np.float64]:
"""
Performs Gaussian Elimination to transform the matrix into Row Echelon Form.
Args:
use_pivoting (bool): If True, swaps rows to put largest absolute value on diagonal.
If False, proceeds blindly
"""
M = self.augmented.copy()
# min accounts for pivot columns for undetermined and overdetermined systems
for k in range(min(self.rows, self.n)):
if use_pivoting:
max_index = np.argmax(np.abs(M[k:, k])) + k
if k != max_index:
M[[k, max_index]] = M[[max_index, k]]
pivot = M[k, k]
if np.isclose(pivot, 0):
continue
for i in range(k+1, self.rows):
factor = M[i, k] / pivot
M[i, k:] -= factor * M[k, k:]
# To account for float ops quirkness
if abs(M[i, k]) < 1e-10:
M[i, k] = 0.0
return M
def back_substitution(self, ref_matrix: NDArray[np.float64]):
"""
Perform Back Substitution on a Row Echelon Form Matrix.
Returns the solution vector or raise error on Infinite/No solution.
"""
x = np.zeros(self.n)
rows = self.rows
for i in range(self.n-1, -1, -1):
rows-= 1
if np.allclose(ref_matrix[rows, :-1], 0) and not np.isclose(ref_matrix[rows, -1], 0):
return "Inconsistent system (No Solution)"
pivot = ref_matrix[rows, i]
if np.isclose(pivot, 0):
# Dependent system (Infinitely Many Solution), meaning x_i is a free variable.
# pick a constant value for it, say 0.0 or 1.0
x[i] = random.choice([0.0, 1.0, 2.0])
print(f"free variable, x[{i}]: {x[i]}")
continue
sum_ax = np.dot(ref_matrix[rows, i+1:self.n], x[i+1:])
x[i] = (ref_matrix[rows, -1] - sum_ax) / pivot
return x