Gaussian Elimination: the algorithm hiding inside NumPy that I was doing by hand

There's a specific moment in studying math that hits different as an engineer: when you realize the "academic exercise" you're grinding through is literally running inside production software you've used for years.

That moment happened to me recently. I've been pivoting from backend engineering (TypeScript, NestJS, distributed systems) into AI Engineering, and I decided I wasn't going to fake my way through the math. No skipping the foundations. So I went back to Gilbert Strang's MIT 18.06 and started solving linear systems by hand. And then it clicked.

The Setup

I was working through a 3×3 system:

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

Which becomes an augmented matrix:

[ 1  2 -1 | 3 ]
[ 2  1  1 | 7 ]
[ 3 -1  2 | 8 ]

The goal: zero out everything below the diagonal. Pivot by pivot.

First pivot (column 1):

L2 ← L2 - 2·L1  →  [ 0  -3   3 |  1 ]
L3 ← L3 - 3·L1  →  [ 0  -7   5 | -1 ]

Second pivot (column 2):

L3 ← 3·L3 - 7·L2  →  [ 0  0  -6 | -10 ]

Upper triangular form:

[ 1  2  -1 |   3 ]
[ 0 -3   3 |   1 ]
[ 0  0  -6 | -10 ]

Back-substitution from bottom to top gives:

z = 5/3,  y = 4/3,  x = 2

Standard stuff. Nothing fancy. Or so I thought.

The multiplier m

Every elimination step computes a multiplier:

m = element_to_zero / pivot

So when zeroing out L2[0] using L1 as pivot row:

m = 2/1 = 2  →  L2 ← L2 - 2·L1

For L3[0]:

m = 3/1 = 3  →  L3 ← L3 - 3·L1

I was doing this mechanically, column by column, treating each operation as:

A[i][j] = A[i][j] - m * A[pivot][j]

That's it. That's the whole thing. And that's when I looked at an algorithm and went quiet for a second.

This is literally the code

for pivot in range(n):
    for row in range(pivot + 1, n):
        m = A[row][pivot] / A[pivot][pivot]
        A[row] = A[row] - m * A[pivot]

The exact sequence I was doing by hand, pivot selection, multiplier computation, row update, is the algorithm. Not a simplification of it. Not a conceptual analogy. The actual algorithm.

And when you call np.linalg.solve(A, b), you're running a production-grade and optimized version of this. The math is the same. The performance engineering around it is what makes it fast.

Where it goes from here

NumPy doesn't literally run Gaussian Elimination in the naive textbook form. What it actually computes under the hood is LU decomposition, a factorization of the matrix into two triangular pieces, where U is essentially what we produced with elimination, and L stores the multipliers m along the way.

I haven't gone deep into LU yet. But understanding that the elimination I was doing by hand is the entry point to that decomposition changed how I see the abstraction. It's not magic. It's the same loop, formalized.

What this study session actually changed

I came in thinking I was filling a gap in my math background. I came out understanding something structural: the linear algebra I'm studying isn't background knowledge for ML, it is the substrate of ML.

Backprop is the chain rule applied to matrix operations. Attention in transformers is matrix multiplication with a softmax. Embeddings live in vector spaces where distance and similarity are defined by inner products. The gradient descent step is a vector subtraction.

When Gilbert Strang says "the key ideas of linear algebra" he's not being poetic. Those ideas are load-bearing walls in almost every ML system.

I'm still early in this path, backend engineer moving into AI Engineering, currently building and studying simultaneously. But I'm increasingly convinced that the engineers who understand what's happening inside np.linalg.solve will make better decisions than the ones who only know how to call it.

I'm documenting this pivot publicly. My RAG project is live at buscadegeloefogo.vercel.app, the Linear Algebra visualizer I built as a study tool is at github.com/FelipeAraujoBS/LA-Canva-Playground. More posts incoming as I go deeper.