---
layout: post
title: "Numerical Gradient Checking: How to Debug Your Autograd Engine Before Training GPT"
date: 2026-06-07
tags: ["Deep Learning", "Backpropagation", "Autograd", "Python", "From Scratch", "Testing"]
category: "ai"
description: "What is numerical gradient checking? Learn the central difference formula, why you never train with finite differences, and how to build a grad checker that catches bugs in your autograd before LayerNorm or attention."
keywords: "numerical gradient checking, gradient check autograd, finite difference gradient, central difference formula, debug backpropagation, gradcheck from scratch, test autograd engine"
image: /assets/images/og-gradient-checking.png
order: 3
---

Most people learn numerical gradients as a curiosity — a slow fallback when you can't do calculus. That's backwards. I used to think the same until I started building an autograd engine from scratch for [AL-1](/blog/build-autograd-from-scratch/). The point of numerical gradient checking isn't to compute gradients during training. It's to catch bugs in your `backward()` functions before you waste weeks debugging attention, LayerNorm, or AdamW when the real problem is a single wrong partial derivative three layers down.

If your autograd is wrong, your entire GPT trains incorrectly. Numerical checking is the unit test framework for that engine. Below I'll derive the central difference formula step by step, work a complete numeric example on $f(w) = w^2$, compare analytical vs finite-difference gradients on a tiny neural net, and show the Python `gradcheck` I run after every op in my [autograd build](/blog/build-autograd-from-scratch/).

---

## What is a numerical gradient?

A gradient is a slope. For a function of one variable, it's "how much does $f$ change when I nudge the input a tiny bit?"

Take the simplest function I know:

$$
f(w) = w^2
$$

At $w = 3$, calculus gives the **analytical** derivative:

$$
\frac{df}{dw} = 2w
$$

$$
\frac{df}{dw}\bigg|_{w=3} = 2 \times 3 = 6
$$

I don't need calculus to *estimate* that slope. I can measure it by evaluating $f$ at two nearby points and asking: how much did $f$ rise per unit of horizontal step?

Pick a small step size:

$$
h = 0.0001
$$

Evaluate one step to the right:

$$
f(3 + h) = f(3.0001) = (3.0001)^2
$$

$$
(3.0001)^2 = 9.00060001
$$

Evaluate one step to the left:

$$
f(3 - h) = f(2.9999) = (2.9999)^2
$$

$$
(2.9999)^2 = 8.99940001
$$

The rise in $f$ across the interval:

$$
f(3 + h) - f(3 - h) = 9.00060001 - 8.99940001
$$

$$
f(3 + h) - f(3 - h) = 0.0012
$$

The horizontal span of that interval:

$$
(3 + h) - (3 - h) = 2h = 0.0002
$$

Divide rise by run:

$$
\frac{f(3 + h) - f(3 - h)}{2h} = \frac{0.0012}{0.0002}
$$

$$
\frac{f(3 + h) - f(3 - h)}{2h} = 6
$$

That matches the analytical answer exactly. This estimate — using points on *both* sides of $w$ — is the **central difference formula**.

```mermaid
graph LR
  subgraph "Central difference at w=3"
    L["f(3−h)\n8.9994"] --- C["w = 3\nslope ≈ 6"]
    C --- R["f(3+h)\n9.0006"]
  end
  style C fill:#d4edda,stroke:#009900,color:#111
```

---

## How do you derive the central difference formula?

Start from the definition of the derivative:

$$
f'(w) = \lim_{h \to 0} \frac{f(w + h) - f(w)}{h}
$$

That's the **forward difference** — one point ahead, one at $w$. A symmetric version uses a point on each side:

$$
f'(w) = \lim_{h \to 0} \frac{f(w + h) - f(w - h)}{2h}
$$

Why is the symmetric version better? Taylor-expand $f(w + h)$ around $w$:

$$
f(w + h) = f(w) + h f'(w) + \frac{h^2}{2} f''(w) + O(h^3)
$$

Taylor-expand $f(w - h)$:

$$
f(w - h) = f(w) - h f'(w) + \frac{h^2}{2} f''(w) + O(h^3)
$$

Subtract the second from the first:

$$
f(w + h) - f(w - h) = 2h f'(w) + O(h^3)
$$

The $f(w)$ terms cancel.

The $\frac{h^2}{2} f''(w)$ terms also cancel — one positive, one negative.

What remains:

$$
f(w + h) - f(w - h) = 2h f'(w) + \text{(terms of order } h^3 \text{ and higher)}
$$

Divide both sides by $2h$:

$$
\frac{f(w + h) - f(w - h)}{2h} = f'(w) + O(h^2)
$$

The error is $O(h^2)$ — proportional to $h^2$, not $h$.

For forward difference, the analogous expansion leaves an $O(h)$ error term.

That means central difference gets accurate answers with a *larger* $h$, which matters when floating-point rounding noise starts fighting you.

Almost every gradient checker uses:

$$
\boxed{\frac{f(w + h) - f(w - h)}{2h}}
$$

---

## How does analytical vs numerical gradient compare on $f(w) = w^2$?

I'll run the full comparison at $w = 3$ with $h = 0.0001$ and write every number on its own line.

**Analytical path (calculus):**

$$
f(w) = w^2
$$

$$
\frac{df}{dw} = 2w
$$

$$
w = 3
$$

$$
g_{\text{analytical}} = 2 \times 3 = 6
$$

**Numerical path (central difference):**

$$
h = 0.0001
$$

$$
w + h = 3.0001
$$

$$
w - h = 2.9999
$$

$$
f(w + h) = (3.0001)^2 = 9.00060001
$$

$$
f(w - h) = (2.9999)^2 = 8.99940001
$$

$$
f(w + h) - f(w - h) = 9.00060001 - 8.99940001 = 0.0012
$$

$$
2h = 2 \times 0.0001 = 0.0002
$$

$$
g_{\text{numerical}} = \frac{0.0012}{0.0002} = 6
$$

**Comparison:**

$$
|g_{\text{analytical}} - g_{\text{numerical}}| = |6 - 6| = 0
$$

On this simple function, the two paths agree to machine precision. That's the pass condition I'd expect from a correct `backward()`.

---

## What does numerical gradient checking look like on a tiny neural net?

Calculus on $f(w) = w^2$ is a sanity check for the formula. The real use case is verifying autograd on an actual computation graph.

Consider one multiply — the smallest "network" I can write:

```python
w = 2.0
x = 3.0

y = w * x
loss = y
```

### What should autograd return? (analytical)

Forward pass:

$$
y = w \times x
$$

$$
y = 2 \times 3 = 6
$$

$$
\text{loss} = y = 6
$$

Backprop from the loss:

$$
\frac{\partial \text{loss}}{\partial y} = 1
$$

Local derivative of multiply with respect to $w$:

$$
\frac{\partial y}{\partial w} = x = 3
$$

Chain rule:

$$
\frac{\partial \text{loss}}{\partial w} = \frac{\partial \text{loss}}{\partial y} \times \frac{\partial y}{\partial w}
$$

$$
\frac{\partial \text{loss}}{\partial w} = 1 \times 3 = 3
$$

Autograd should set:

$$
\texttt{w.grad} = 3
$$

### What does central difference return? (numerical)

Treat loss as a function of $w$ alone, holding $x$ fixed:

$$
L(w) = w \times x = 3w
$$

At $w = 2$ with $h = 0.0001$:

$$
L(w + h) = L(2.0001) = 2.0001 \times 3
$$

$$
L(2.0001) = 6.0003
$$

$$
L(w - h) = L(1.9999) = 1.9999 \times 3
$$

$$
L(1.9999) = 5.9997
$$

$$
L(w + h) - L(w - h) = 6.0003 - 5.9997
$$

$$
L(w + h) - L(w - h) = 0.0006
$$

$$
2h = 0.0002
$$

$$
\frac{L(w + h) - L(w - h)}{2h} = \frac{0.0006}{0.0002}
$$

$$
\frac{L(w + h) - L(w - h)}{2h} = 3
$$

Analytical: $3$. Numerical: $3$. Pass.

---

## What are we actually comparing after `loss.backward()`?

After `loss.backward()`, autograd gives me $\frac{\partial L}{\partial w}$ — say:

$$
\texttt{w.grad} = 3.000001
$$

Numerical checking gives the same derivative from finite differences — say:

$$
\texttt{w\_numeric} = 2.999999
$$

Absolute error:

$$
|3.000001 - 2.999999| = 0.000002
$$

Tiny. The backward implementation passes.

```mermaid
sequenceDiagram
  participant A as Autograd backward()
  participant N as Numerical checker
  participant C as Compare
  A->>C: w.grad = 3.000001
  N->>C: w_numeric = 2.999999
  C->>C: |error| < tolerance → pass
```

In practice I compare **relative** error, not just absolute — but the idea is the same: two independent paths to the same derivative should agree.

---

## Why can't you train with numerical gradients?

Because they're insanely expensive.

Take a weight matrix `W.shape = (768, 768)`. Parameter count:

$$
768 \times 768 = 589{,}824
$$

For each parameter I need $L(\theta_i + h)$ and $L(\theta_i - h)$ — two forward passes per parameter.

$$
2 \times 589{,}824 = 1{,}179{,}648 \text{ forward passes}
$$

For **one** gradient computation. On one layer.

Backprop does it in one forward + one backward. The gap isn't a constant factor — it's proportional to the number of parameters. I'd never finish training a GPT this way.

| Method | Passes per gradient step | Exact? | Use case |
|---|---|---|---|
| Numerical (central diff) | $2 \times N_{\text{params}}$ | Approximate | Debugging |
| Reverse-mode autograd | 1 forward + 1 backward | Exact (mod float) | Training |

That's why every framework trains with autograd and only reaches for finite differences when something smells wrong.

---

## How do you build a `gradcheck` in Python?

Suppose I have a `Tensor` with `.data` and `.grad`. After `loss.backward()`, I check each element. Here is the loop I actually run — broken into the six steps I think through every time.

### Step 1 — store the original value

```python
original = w.data[i]
```

### Step 2 — perturb up, run forward

```python
w.data[i] = original + h
loss_plus = forward()
```

### Step 3 — perturb down, run forward

```python
w.data[i] = original - h
loss_minus = forward()
```

### Step 4 — restore

```python
w.data[i] = original
```

### Step 5 — compute numerical gradient

$$
\text{grad\_numeric} = \frac{\text{loss\_plus} - \text{loss\_minus}}{2h}
$$

```python
grad_numeric = (loss_plus - loss_minus) / (2 * h)
```

### Step 6 — compare

```python
grad_auto = w.grad[i]
error = abs(grad_auto - grad_numeric)
```

For parameter $\theta_i$, the full comparison is:

$$
\text{numerical:} \quad \boxed{\frac{L(\theta_i + h) - L(\theta_i - h)}{2h}} \qquad \text{vs} \qquad \text{autograd:} \quad \boxed{\frac{\partial L}{\partial \theta_i}}
$$

Wrapped into a reusable function:

```python
def gradcheck(forward_fn, params, h=1e-5, tol=1e-5):
    """Compare autograd gradients to central-difference estimates."""
    # Run forward + backward once to populate .grad
    loss = forward_fn()
    loss.backward()

    for param in params:
        flat_data = param.data.ravel()
        flat_grad = param.grad.ravel()
        for i in range(flat_data.size):
            original = flat_data[i]

            flat_data[i] = original + h
            loss_plus = forward_fn().data

            flat_data[i] = original - h
            loss_minus = forward_fn().data

            flat_data[i] = original  # always restore

            grad_numeric = (loss_plus - loss_minus) / (2 * h)
            grad_auto = flat_grad[i]

            denom = max(1.0, abs(grad_auto), abs(grad_numeric))
            rel_error = abs(grad_auto - grad_numeric) / denom

            if rel_error > tol:
                raise AssertionError(
                    f"gradcheck failed at index {i}: "
                    f"auto={grad_auto}, numeric={grad_numeric}, "
                    f"rel_error={rel_error}"
                )
    return True
```
{: data-lineno="true"}

```mermaid
flowchart TD
  S["Store θᵢ"] --> P["θᵢ ← θᵢ + h\nforward → L₊"]
  P --> M["θᵢ ← θᵢ − h\nforward → L₋"]
  M --> R["Restore θᵢ"]
  R --> N["(L₊ − L₋) / 2h"]
  N --> C{"|auto − numeric|\n< tol?"}
  C -->|yes| OK["✓ op passes"]
  C -->|no| FAIL["✗ fix backward()"]
  style OK fill:#d4edda,stroke:#009900,color:#111
  style FAIL fill:#f8d7da,stroke:#990000,color:#111
```

Loop over every element of every parameter tensor. Yes, it's slow. That's the point — I run this in tests, not in the training loop.

---

## What should you test before you build GPT?

I'm working through AL-1 right now — a from-scratch GPT in pure NumPy. Before I touch LayerNorm, attention, softmax, cross-entropy, or AdamW, every primitive op needs to pass `gradcheck(forward_fn, params)`:

| Operation | Forward | Expected $\partial z / \partial x$ |
|---|---|---|
| Add | $z = x + y$ | $1$ |
| Multiply | $z = xy$ | $y$ (and $\partial z / \partial y = x$) |
| ReLU | $\max(0, x)$ | $0$ or $1$ |
| Exp | $e^x$ | $e^x$ |
| Log | $\log(x)$ | $1/x$ |
| Sum | $\sum x_i$ | $1$ (broadcast back) |
| MatMul | $AB$ | Most common source of bugs — check thoroughly |

MatMul is where I expect pain. Transpose conventions, batch dimensions, broadcasting in the backward pass — one wrong axis swap and every downstream gradient is garbage. I'd rather find that here with a 2×2 matrix than inside a 12-layer transformer.

If the checker fails on `add`, every later component inherits the bug. Fix the foundation first.

---

## How do you choose epsilon and measure error?

For $f(x) = x^3$ at $x = 2$, the exact derivative is:

$$
\frac{df}{dx} = 3x^2
$$

$$
3 \times 2^2 = 3 \times 4 = 12
$$

With $\epsilon = 0.01$:

$$
x + \epsilon = 2.01
$$

$$
x - \epsilon = 1.99
$$

$$
f(2.01) = (2.01)^3 = 8.120601
$$

$$
f(1.99) = (1.99)^3 = 7.880599
$$

$$
f(2.01) - f(1.99) = 8.120601 - 7.880599 = 0.240002
$$

$$
2\epsilon = 2 \times 0.01 = 0.02
$$

$$
\frac{f(2.01) - f(1.99)}{2(0.01)} = \frac{0.240002}{0.02}
$$

$$
\frac{f(2.01) - f(1.99)}{0.02} = 12.0001
$$

Absolute error:

$$
|12 - 12.0001| = 0.0001
$$

Relative error (the comparison I actually use in `gradcheck`):

$$
\frac{|g_a - g_n|}{\max(1, |g_a|, |g_n|)} = \frac{0.0001}{\max(1, 12, 12.0001)}
$$

$$
\frac{0.0001}{12.0001} \approx 8.3 \times 10^{-6}
$$

Too-large $\epsilon$ causes truncation error because the curve is not locally linear across the interval. Too-small $\epsilon$ causes floating-point cancellation because $f(x + \epsilon)$ and $f(x - \epsilon)$ become almost equal. Around $10^{-5}$ to $10^{-4}$ is a useful starting range for float64; float32 often needs a larger step and looser tolerance.

I run checks in deterministic float64 mode, disable dropout, avoid points where ReLU or max is nondifferentiable, sample a few coordinates of huge tensors, and always restore perturbed parameters. Gradient checking validates derivative code; it does not prove that the forward equation, data, or objective is conceptually correct.

---

## Where does this fit in the build?

This post is the sanity check that sits between understanding [backprop from scratch](/blog/build-autograd-from-scratch/) and actually trusting my engine on real model code. I run `gradcheck` after implementing each op's `backward()`, before moving to the next one. When all seven pass, I stop second-guessing whether my gradients are right and start worrying about things that actually matter — data pipelines, tokenizer bugs, learning rate schedules.

Numerical gradient checking won't make my model smarter. But it will save me from the worst kind of debugging session: staring at a loss curve that looks almost right, tweaking hyperparameters that were never the problem, when a single botched `backward()` has been lying to me since step zero.

Build the checker first. Trust nothing until it passes.