Math for AI course · Lesson 23 of 180

Gaussian Elimination Preserves the Solution Set

Recover model parameters from a zero-pivot system, then implement safe row operations and residual checks in NumPy.

Share this article

A parameter-recovery system can contain the right information and still begin with an inconvenient zero. If the first pivot is zero, the safe response is not to invent a new equation or divide by zero. Rewrite the same constraints until their structure becomes visible.

That is the promise of Gaussian elimination: use reversible elementary row operations to simplify a linear system while preserving exactly the same solution set. By the end of this lesson, you will be able to apply the three legal row operations, solve a small parameter-recovery problem, implement the operations in NumPy, and diagnose zero pivots, redundant rows, contradictions, and floating-point residuals.

Learning outcome: Apply elementary row operations without changing solutions.

Lesson 13: No solution, one solution, or many solutions classified linear systems by their solution sets. Lesson 18: Matrices organize linear rules gave us a compact way to store the coefficients. The official Mathematics for Machine Learning companion places the next step in Section 2.3 of its current book PDF. The source explains elementary transformations and augmented matrices; the numbers, parameter scenario, code, visual, and assessments here are original.

The goal is a simpler system, not a different question

Consider a system with mm linear equations and nn unknowns. We can write it as

Aθ=b,A\boldsymbol{\theta}=\boldsymbol{b},

where A∈Rm×nA\in\mathbb{R}^{m\times n} is the coefficient matrix, θ∈Rn\boldsymbol{\theta}\in\mathbb{R}^{n} contains the unknown parameters, and b∈Rm\boldsymbol{b}\in\mathbb{R}^{m} contains the observed right-hand sides. Each row of the matrix equation represents one scalar equation:

ai1θ1+ai2θ2+⋯+ainθn=bi.a_{i1}\theta_1+a_{i2}\theta_2+\cdots+a_{in}\theta_n=b_i.

Gaussian elimination changes the appearance of those rows so that pivots and unknowns are easier to read. It must not change which parameter vectors make all the equations true. The invariant is the set

S={θ∈Rn:Aθ=b}.\mathcal{S}=\left\{\boldsymbol{\theta}\in\mathbb{R}^{n}: A\boldsymbol{\theta}=\boldsymbol{b}\right\}.

The matrix notation becomes more useful when we attach the right-hand side to the coefficient row:

[a11a12⋯a1nb1a21a22⋯a2nb2⋮⋮⋮⋮am1am2⋯amnbm].\left[\begin{array}{ccc|c} a_{11}&a_{12}&\cdots&a_{1n}&b_1\\ a_{21}&a_{22}&\cdots&a_{2n}&b_2\\ \vdots&\vdots&&\vdots&\vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn}&b_m \end{array}\right].

This is an augmented matrix. The vertical divider is a reminder that the last column is not another unknown coefficient. It is the output of the equation, and it must travel with its row whenever that row is changed.

Three reversible moves keep the solution set fixed

Let RiR_i denote one complete row of an augmented matrix. There are exactly three elementary row operations:

Legal operationNotationWhy it is reversible
Exchange two rowsRi↔RjR_i\leftrightarrow R_jExchange them again.
Multiply a row by a nonzero scalarRi←λRiR_i\leftarrow\lambda R_i, λ≠0\lambda\ne0Multiply by 1/λ1/\lambda.
Add a multiple of one row to anotherRi←Ri+λRjR_i\leftarrow R_i+\lambda R_j, i≠ji\ne jSubtract λRj\lambda R_j from RiR_i.

The third operation is often called row replacement. It includes adding one equation to another, subtracting one equation from another, and clearing a coefficient with a carefully chosen multiple.

The nonzero condition matters. Turning RiR_i into 0Ri0R_i erases a constraint; there is no inverse operation that recovers it. The same full-row rule matters too: multiplying coefficients without multiplying the corresponding right-hand-side entry creates a new equation.

Why the augmented matrix carries the whole equation

Suppose one row represents

θ1+2θ2=5.\theta_1+2\theta_2=5.

Its augmented row is [1    2∣5][1\;\;2\mid5]. Scaling the equation by 33 means

3(θ1+2θ2)=3(5)⟺3θ1+6θ2=15,3(\theta_1+2\theta_2)=3(5) \quad\Longleftrightarrow\quad 3\theta_1+6\theta_2=15,

so the legal matrix operation is

[1    2∣5]⟶[3    6∣15].[1\;\;2\mid5]\longrightarrow[3\;\;6\mid15].

The row [3    6∣5][3\;\;6\mid5] would instead describe 3θ1+6θ2=53\theta_1+6\theta_2=5. It is a different constraint. Gaussian elimination is not permission to edit the left side and right side independently; it is a controlled rewrite of the complete equation.

The same principle explains row replacement. If

Ri:aiTθ=bi,Rj:ajTθ=bj,\begin{aligned} R_i &: \boldsymbol{a}_i^{\mathsf T}\boldsymbol{\theta}=b_i,\\ R_j &: \boldsymbol{a}_j^{\mathsf T}\boldsymbol{\theta}=b_j, \end{aligned}

then replacing RiR_i with Ri+λRjR_i+\lambda R_j means

(ai+λaj)Tθ=bi+λbj.(\boldsymbol{a}_i+\lambda\boldsymbol{a}_j)^{\mathsf T} \boldsymbol{\theta} =b_i+\lambda b_j.

Every coefficient and the right-hand side receive the same combination. The new row is easier to use, but it still encodes information already present in the old pair.

Worked example: recover three scoring parameters

Imagine a toy linear score built from three measured signals:

  • f1f_1: a retrieval-match feature;
  • f2f_2: a verification feature; and
  • f3f_3: a tool-confirmation feature.

The score rule is

s=θ1f1+θ2f2+θ3f3.s=\theta_1f_1+\theta_2f_2+\theta_3f_3.

The counts and scores below are deliberately small. They are an instructional parameter-identification problem, not a claim that a production ranking system should use these features or weights.

Observationf1f_1f2f_2f3f_3observed score ss
A0215
B1114
C24515

The three rows become

0θ1+2θ2+θ3=5,θ1+θ2+θ3=4,2θ1+4θ2+5θ3=15.\begin{aligned} 0\theta_1+2\theta_2+\theta_3&=5,\\ \theta_1+\theta_2+\theta_3&=4,\\ 2\theta_1+4\theta_2+5\theta_3&=15. \end{aligned}

In matrix form,

[0215111424515].\left[\begin{array}{ccc|c} 0&2&1&5\\ 1&1&1&4\\ 2&4&5&15 \end{array}\right].

The first entry is zero, so the first row cannot be scaled into a first pivot. That is a bookkeeping problem, not evidence that the system has no solution. We can swap in a nonzero row. The static-first visual records every state, including the operation label and the complete augmented matrix:

Every row operation in the three-parameter recovery example.

Every row operation in the three-parameter recovery example

Complete static record

Every state remains visible below the optional controls. The operation changes a complete row, so each state represents an equivalent system.

  1. Step 0Start: original augmented matrix
    Original augmented matrix
     θ₁θ₂θ₃b
    R10215
    R21114
    R324515
  2. Step 1R₁ ↔ R₃
    Matrix after step 1
     θ₁θ₂θ₃b
    R124515
    R21114
    R30215
  3. Step 2R₁ ← ½R₁
    Matrix after step 2
     θ₁θ₂θ₃b
    R1125/215/2
    R21114
    R30215
  4. Step 3R₂ ← R₂ − R₁
    Matrix after step 3
     θ₁θ₂θ₃b
    R1125/215/2
    R20-1-3/2-7/2
    R30215
  5. Step 4R₂ ↔ R₃
    Matrix after step 4
     θ₁θ₂θ₃b
    R1125/215/2
    R20215
    R30-1-3/2-7/2
  6. Step 5R₂ ← ½R₂
    Matrix after step 5
     θ₁θ₂θ₃b
    R1125/215/2
    R2011/25/2
    R30-1-3/2-7/2
  7. Step 6R₁ ← R₁ − 2R₂
    Matrix after step 6
     θ₁θ₂θ₃b
    R1103/25/2
    R2011/25/2
    R30-1-3/2-7/2
  8. Step 7R₃ ← R₃ + R₂
    Matrix after step 7
     θ₁θ₂θ₃b
    R1103/25/2
    R2011/25/2
    R300-1-1
  9. Step 8R₃ ← −R₃
    Matrix after step 8
     θ₁θ₂θ₃b
    R1103/25/2
    R2011/25/2
    R30011
  10. Step 9R₁ ← R₁ − 3/2R₃
    Matrix after step 9
     θ₁θ₂θ₃b
    R11001
    R2011/25/2
    R30011
  11. Step 10R₂ ← R₂ − 1/2R₃
    Matrix after step 10
     θ₁θ₂θ₃b
    R11001
    R20102
    R30011

Step through one legal transformation at a time

Focus this panel and use the buttons or the left and right arrow keys to move between states.

Step 0 / 10

OperationStart: original augmented matrix

Current augmented matrix
 θ₁θ₂θ₃b
R10215
R21114
R324515

The first entry in the first column is zero, so scaling that row would be illegal. A nonzero row below can supply the first pivot.

The row labels, operation text, and vertical divider before b carry the meaning without relying on color. The static sequence is the fallback if JavaScript is unavailable.

Here is the same calculation in algebraic order. We use partial pivoting in the first two columns: when several nonzero candidates are available, choose the largest magnitude candidate in the current column. The lesson’s visual keeps the fractions visible so each operation can be checked by hand.

The recovered parameters pass a direct substitution check:

[021111245][121]=[0+4+11+2+12+8+5]=[5415]=b.\begin{aligned} \begin{bmatrix} 0&2&1\\ 1&1&1\\ 2&4&5 \end{bmatrix} \begin{bmatrix}1\\2\\1\end{bmatrix} &= \begin{bmatrix} 0+4+1\\ 1+2+1\\ 2+8+5 \end{bmatrix}\\ &= \begin{bmatrix}5\\4\\15\end{bmatrix} =\boldsymbol{b}. \end{aligned}

The solution did not appear because a row operation guessed new parameters. Each step changed the representation of the same three constraints. The final identity block makes the answer easy to read.

A NumPy implementation mirrors the mathematics

The following implementation makes the three legal operations explicit. It works on a copy of a finite, real-valued two-dimensional augmented array. The last column is treated as the right-hand side, and atol is the threshold below which a floating-point value is treated as zero. The implementation performs Gauss–Jordan elimination, so it clears entries both below and above each pivot.

import numpy as np


def swap_rows(matrix, first, second):
    result = np.array(matrix, dtype=float, copy=True)
    result[[first, second], :] = result[[second, first], :]
    return result


def scale_row(matrix, row, factor):
    if factor == 0:
        raise ValueError("row scaling needs a nonzero factor")
    result = np.array(matrix, dtype=float, copy=True)
    result[row, :] *= factor
    return result


def replace_row(matrix, target, source, factor):
    if target == source:
        raise ValueError("target and source rows must be different")
    result = np.array(matrix, dtype=float, copy=True)
    result[target, :] += factor * result[source, :]
    return result


def gauss_jordan(augmented, atol=1e-12):
    work = np.array(augmented, dtype=float, copy=True)
    if work.ndim != 2 or work.shape[1] < 2:
        raise ValueError("expected a 2D augmented matrix")
    if not np.all(np.isfinite(work)):
        raise ValueError("all entries must be finite")

    row_count, column_count = work.shape
    variable_count = column_count - 1
    pivot_row = 0

    for pivot_column in range(variable_count):
        if pivot_row == row_count:
            break

        candidates = np.abs(work[pivot_row:, pivot_column])
        pivot_offset = int(np.argmax(candidates))
        if candidates[pivot_offset] <= atol:
            continue

        pivot = pivot_row + pivot_offset
        if pivot != pivot_row:
            work = swap_rows(work, pivot_row, pivot)

        pivot_value = work[pivot_row, pivot_column]
        work = scale_row(work, pivot_row, 1.0 / pivot_value)

        for row in range(row_count):
            if row == pivot_row:
                continue
            coefficient = work[row, pivot_column]
            if abs(coefficient) > atol:
                work = replace_row(
                    work, row, pivot_row, -coefficient
                )

        pivot_row += 1

    work[np.abs(work) <= atol] = 0.0
    return work


A = np.array(
    [
        [0.0, 2.0, 1.0],
        [1.0, 1.0, 1.0],
        [2.0, 4.0, 5.0],
    ]
)
b = np.array([5.0, 4.0, 15.0])
augmented = np.column_stack((A, b))

rref = gauss_jordan(augmented)
theta = rref[:, -1]  # valid here because the coefficient block is full rank
residual = A @ theta - b

assert np.allclose(residual, 0.0, rtol=1e-12, atol=1e-12)
assert np.allclose(
    theta,
    np.linalg.solve(A, b),
    rtol=1e-12,
    atol=1e-12,
)

print(theta)
print(np.linalg.norm(residual, ord=np.inf))

The verified output is

[1. 2. 1.]
0.0

The final assertion compares the hand-built elimination result with numpy.linalg.solve, which the current NumPy documentation defines for a square, full-rank system and reports as raising LinAlgError for a singular or non-square coefficient matrix. The comparison is a cross-check, not a reason to hide the row operations: the point of this lesson is to make the equivalence-preserving mechanism inspectable.

The residual

r=Aθ^−b\boldsymbol{r}=A\widehat{\boldsymbol{\theta}}-\boldsymbol{b}

measures how well the computed parameters satisfy the supplied equations. The infinity norm used in the code is the largest absolute residual coordinate. For exact arithmetic, the example gives the zero vector. In floating-point arithmetic, a correct result may instead produce values on the order of machine round-off.

Boundary cases expose what the pivots mean

A zero pivot is a branching point in the algorithm, not a verdict. If the current pivot position is zero and a nonzero entry appears lower in the same column, swap rows. If every entry below is zero, move to the next column; that column may simply contain no new independent information.

After elimination, inspect rows with all-zero coefficients. There are two different outcomes:

  1. An all-zero row with a zero right-hand side, [0    0    0∣0][0\;\;0\;\;0\mid0], is redundant. It adds no constraint. If a variable also lacks a pivot, the system may have infinitely many solutions.
  2. An all-zero row with a nonzero right-hand side, such as [0    0    0∣1][0\;\;0\;\;0\mid1], says 0=10=1. The system is inconsistent and has no solution.

For example, the two rows

[123246]→R2←R2−2R1[123000]\left[\begin{array}{cc|c} 1&2&3\\ 2&4&6 \end{array}\right] \xrightarrow{R_2\leftarrow R_2-2R_1} \left[\begin{array}{cc|c} 1&2&3\\ 0&0&0 \end{array}\right]

encode the same line twice. Changing only the second right-hand side gives a different result:

[123247]→R2←R2−2R1[123001],\left[\begin{array}{cc|c} 1&2&3\\ 2&4&7 \end{array}\right] \xrightarrow{R_2\leftarrow R_2-2R_1} \left[\begin{array}{cc|c} 1&2&3\\ 0&0&1 \end{array}\right],

which is a contradiction. The coefficient rows look dependent in both examples; the right-hand side decides whether the observed constraints agree.

The shortest misconception test is to ask what happens to θ1+θ2=3\theta_1+\theta_2=3 when the left side is multiplied by 22. The legal result is 2θ1+2θ2=62\theta_1+2\theta_2=6, not 2θ1+2θ2=32\theta_1+2\theta_2=3. Row operations preserve a solution set because the complete equations are transformed reversibly.

Why this matters in an AI pipeline

The toy score shows the direct use: several observations can recover unknown weights when the equations identify one parameter vector. In a real pipeline, the matrix may describe calibration weights, a small physical model, a constraint system, or a diagnostic slice of a larger model. Elimination helps answer three practical questions:

  • Do the observations provide a pivot for every parameter?
  • Is an apparently repeated observation merely redundant?
  • Do the measured targets contradict the assumed linear rule?

The method also marks the boundary of the exact problem. Noisy measurements often make the equations inconsistent, and an overdetermined dataset may have no vector that satisfies every row exactly. In that situation, forcing a contradiction away would change the question. A least-squares method asks a different question—find a vector with small aggregate error—and belongs to the next numerical-linear-algebra decisions, not to an exact row-equivalence claim. Likewise, numpy.linalg.solve expects a square full-rank coefficient matrix; a singular or non-square problem needs a method chosen for its actual objective.

The next lesson names the staircase pattern that elimination reveals: pivots, basic variables, free variables, and the shape of a general solution. The invariant to carry forward is simple: every legal row operation changes the view of the system, not the set of answers.

Check your understanding

Question 1

Why is multiplying a row by zero not an elementary row operation, while multiplying it by 5 is legal?

Show the step-by-step solution

Start with a row equation g(θ)=0g(\boldsymbol{\theta})=0. Multiplying by 55 gives 5g(θ)=05g(\boldsymbol{\theta})=0. Because 5≠05\ne0, the original equation can be recovered by multiplying by 1/51/5. The set of solutions is unchanged.

Multiplying by zero gives 0=00=0 for every θ\boldsymbol{\theta}. The original constraint has disappeared, and no multiplication can recover it. The operation is therefore not reversible and can enlarge the solution set.

Question 2

Apply R₁ ← R₁ − 2R₂ to the augmented matrix [[1, 2 | 7], [0, 1 | 3]]. What is the new first row?

Show the step-by-step solution

Use the complete second row, including its right-hand-side entry:

R1=[1    2∣7],R2=[0    1∣3].R_1=[1\;\;2\mid7], \qquad R_2=[0\;\;1\mid3].

Multiply the second row by 22:

2R2=[0    2∣6].2R_2=[0\;\;2\mid6].

Subtract component by component:

R1−2R2=[1−0    2−2∣7−6]=[1    0∣1].R_1-2R_2 =[1-0\;\;2-2\mid7-6] =[1\;\;0\mid1].

The new first row is [1    0∣1][1\;\;0\mid1]. Omitting the right-hand-side subtraction would have produced a different equation.

Question 3

Solve the system 0x + y = 2 and x + y = 3 using a row swap before eliminating.

Show the step-by-step solution

Write the augmented matrix:

[012113].\left[\begin{array}{cc|c} 0&1&2\\ 1&1&3 \end{array}\right].

The first pivot position is zero, so swap the rows:

[113012].\left[\begin{array}{cc|c} 1&1&3\\ 0&1&2 \end{array}\right].

Clear the entry above the second pivot:

R1←R1−R2⟹[101012].R_1\leftarrow R_1-R_2 \quad\Longrightarrow\quad \left[\begin{array}{cc|c} 1&0&1\\ 0&1&2 \end{array}\right].

Reading the identity block gives x=1x=1 and y=2y=2. Substitution checks both original equations: 0(1)+2=20(1)+2=2 and 1+2=31+2=3.

Question 4

What does an all-zero row mean in each case: [0, 0 | 0] and [0, 0 | 4]?

Show the step-by-step solution

The row [0    0∣0][0\;\;0\mid0] represents 0=00=0. It is redundant and does not restrict the variables. If another variable has no pivot, that freedom can produce infinitely many solutions.

The row [0    0∣4][0\;\;0\mid4] represents 0=40=4, which is impossible. One contradictory row makes the whole system have no solution, even if other rows look consistent.

The difference is entirely in the right-hand side. Elimination must carry that entry along so the algorithm can distinguish redundancy from contradiction.

Question 5

A developer transforms [1, 1 | 3] into [2, 2 | 3] and says this is just row scaling. Is the claim correct?

Show the step-by-step solution

Row scaling by 22 must multiply every entry:

2[1    1∣3]=[2    2∣6].2[1\;\;1\mid3]=[2\;\;2\mid6].

The claimed row [2    2∣3][2\;\;2\mid3] represents 2θ1+2θ2=32\theta_1+2\theta_2=3. Dividing that equation by 22 gives θ1+θ2=3/2\theta_1+\theta_2=3/2, not the original θ1+θ2=3\theta_1+\theta_2=3. The claim is false because the right-hand side was changed independently.

Question 6

Why does the NumPy implementation use a tolerance when deciding whether a pivot or residual is zero?

Show the step-by-step solution

In exact arithmetic, zero is an exact value. Floating-point operations store many fractions approximately, so a mathematically zero result can appear as a tiny value such as 10−1610^{-16}. Treating every nonzero bit pattern as a genuine pivot can create a meaningless row operation.

The implementation therefore compares magnitudes with atol. It also checks the residual with relative and absolute tolerances. The tolerance must match the scale of the data: a threshold that is sensible for values near 11 may be too large for measurements near 10−910^{-9} or too small for very large values. A small residual is evidence of fit, not proof that an ill-conditioned parameter estimate is stable.

Question 7

A recovery system has three unknown parameters but elimination ends with only two pivots and one all-zero row. What should you conclude before choosing a numerical method?

Show the step-by-step solution

Two pivots constrain two independent parameter directions. The all-zero row adds no new information, so at least one parameter direction is free unless a different row reveals a contradiction.

If no row has the form [0    0    0∣c][0\;\;0\;\;0\mid c] with c≠0c\ne0, the system is consistent and has infinitely many solutions. The data do not identify one unique parameter vector under this model. Adding arbitrary numbers to the missing pivot would be a modeling choice, not a consequence of elimination; the next step is to report the free parameter or collect an independent measurement.

Sources

  1. Mathematics for Machine Learning companion website
  2. Mathematics for Machine Learning book PDF
  3. NumPy documentation: numpy.linalg.solve
  4. NumPy documentation: numpy.allclose
  5. NumPy documentation: numpy.finfo
  6. NumPy documentation: numpy.linalg.norm
  7. NumPy documentation: numpy.linalg.cond