We know a system $A\mathbf{x} = \mathbf{b}$ has one, none, or infinitely many solutions. This page is about the algorithm that actually finds them: Gaussian elimination. It’s the engine under np.linalg.solvenp.linalg.solve, under every regression fit, and under matrix inversion. Master it once and you’ll never be mystified by a “singular matrix” again.

The idea: simplify without changing the answer

You already do this by hand: add one equation to another, scale an equation, swap two equations. None of these changes the solution set — they just make the system easier to read. Written on the coefficients alone, they’re called elementary row operations:

1. Swap two rows.
2. Scale a row by a nonzero constant $\lambda$.
3. Add a multiple of one row to another.

To avoid rewriting the variables every step, we stack $A$ and $\mathbf{b}$ into one augmented matrix $[\,A \mid \mathbf{b}\,]$ and operate on that.

Watch it happen

Here is a real $3\times 3$ system reduced to row-echelon form (a triangular “staircase”), then back-substituted to the answer. Each frame is one elementary row operation — watch the zeros appear below the diagonal and the pivots (amber) lock in:

The triangular shape is the whole point. Once the matrix is in row-echelon form, the last row gives one variable outright, and you substitute upward for the rest. That two-phase dance — forward elimination to a staircase, then back-substitution — is Gaussian elimination.

Row-echelon form, pivots, and free variables

A matrix is in row-echelon form (REF) when:

• every all-zero row sits at the bottom, and
• each row’s first nonzero entry (its pivot) is strictly to the right of the pivot above it.

That gives the staircase. If we go further so that every pivot is 1 and is the only nonzero entry in its column, we have reduced row-echelon form (RREF) — the cleanest possible version of the system.

Variables sitting in pivot columns are basic variables; the rest are free variables — we can set them to anything, and each choice gives another solution. Free variables are exactly why a system has infinitely many solutions.

Particular + general solution

When a system has infinitely many solutions, we describe all of them in two pieces:

The recipe is always the same three steps:

  flowchart TD
  A["1. Find ONE particular solution
to Ax = b"] --> C["3. General solution =
particular + all of the homogeneous"]
  B["2. Find ALL solutions
to Ax = 0 (the null space)"] --> C

The set of solutions to the homogeneous system $A\mathbf{x} = \mathbf{0}$ is called the kernel or null space — we’ll meet it again on the Linear Mappings page.

Inversion is just elimination in disguise

Want $A^{-1}$? Solve $A X = I$ for the whole matrix $X$ at once. Augment $A$ with the identity and row-reduce until the left block becomes the identity — whatever lands on the right is the inverse:

So “find the inverse” and “solve a system” are the same computation.

NumPy: solve, invert, and least-squares

import numpy as np
 
A = np.array([[1.0, 1.0, 1.0],
              [1.0, 2.0, 3.0],
              [1.0, 1.0, 2.0]])
b = np.array([6.0, 14.0, 9.0])
 
# 1. Unique solution (square, invertible)
x = np.linalg.solve(A, b)
print("solution:", np.round(x, 4))          # [1. 2. 3.]
 
# 2. Inverse via elimination (NumPy does the row-reduction for you)
print("inverse:\n", np.round(np.linalg.inv(A), 4))
 
# 3. Overdetermined system (more equations than unknowns) -> least squares
#    There is no exact solution, so we find the best fit.
M = np.array([[1.0, 1.0],
              [1.0, 2.0],
              [1.0, 3.0]])
y = np.array([1.0, 2.0, 2.0])
coef, *_ = np.linalg.lstsq(M, y, rcond=None)
print("least-squares fit:", np.round(coef, 4))

import numpy as np
 
A = np.array([[1.0, 1.0, 1.0],
              [1.0, 2.0, 3.0],
              [1.0, 1.0, 2.0]])
b = np.array([6.0, 14.0, 9.0])
 
# 1. Unique solution (square, invertible)
x = np.linalg.solve(A, b)
print("solution:", np.round(x, 4))          # [1. 2. 3.]
 
# 2. Inverse via elimination (NumPy does the row-reduction for you)
print("inverse:\n", np.round(np.linalg.inv(A), 4))
 
# 3. Overdetermined system (more equations than unknowns) -> least squares
#    There is no exact solution, so we find the best fit.
M = np.array([[1.0, 1.0],
              [1.0, 2.0],
              [1.0, 3.0]])
y = np.array([1.0, 2.0, 2.0])
coef, *_ = np.linalg.lstsq(M, y, rcond=None)
print("least-squares fit:", np.round(coef, 4))

solution: [1. 2. 3.]
inverse:
 [[ 1.  1. -1.]
 [ 1. -1.  0.]
 [-1.  0.  1.]]
least-squares fit: [0.6667 0.5   ]

solution: [1. 2. 3.]
inverse:
 [[ 1.  1. -1.]
 [ 1. -1.  0.]
 [-1.  0.  1.]]
least-squares fit: [0.6667 0.5   ]

Why this matters for ML

• Linear regression in closed form is $\hat{\boldsymbol\theta} = (X^\top X)^{-1}X^\top \mathbf{y}$ — a pseudo-inverse, i.e. Gaussian elimination on a normal-equations system.
• Singular matrixSingular matrix errors mean a pivot went missing: redundant/collinear features. The fix (regularization, dropping features) is this theory applied.
• Numerical stability matters: for millions of variables we don’t invert; we use iterative or factorization methods. But they all rest on the row-reduction ideas here.

🧪 Try It Yourself

Exercise 1 – Solve a 3×3 system

Exercise 2 – Invert by augmenting with the identity

Exercise 3 – Best fit when there’s no exact solution

Recap

• Elementary row operations (swap, scale, add) simplify a system without changing its answer.
• Reduce $[\,A \mid \mathbf{b}\,]$ to row-echelon form (a staircase), then back-substitute.
• Pivot columns → basic variables; the rest → free variables (the source of infinite solutions).
• Every solution = particular solution + null space; inverting $A$ is reducing $[\,A\mid I\,]$.
• When no exact solution exists, the pseudo-inverse / least squares gives the best fit — this is linear regression.

Next: we zoom out from single systems to the whole arena vectors live in — Vector Spaces.