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. 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.
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:
At $w = 3$, calculus gives the analytical derivative:
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:
Evaluate one step to the right:
Evaluate one step to the left:
The rise in $f$ across the interval:
The horizontal span of that interval:
Divide rise by run:
That matches the analytical answer exactly. This estimate — using points on both sides of $w$ — is the central difference formula.
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:
That’s the forward difference — one point ahead, one at $w$. A symmetric version uses a point on each side:
Why is the symmetric version better? Taylor-expand $f(w + h)$ around $w$:
Taylor-expand $f(w - h)$:
Subtract the second from the first:
The $f(w)$ terms cancel.
The $\frac{h^2}{2} f’‘(w)$ terms also cancel — one positive, one negative.
What remains:
Divide both sides by $2h$:
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:
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):
Numerical path (central difference):
Comparison:
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:
w = 2.0
x = 3.0
y = w * x
loss = y
What should autograd return? (analytical)
Forward pass:
Backprop from the loss:
Local derivative of multiply with respect to $w$:
Chain rule:
Autograd should set:
What does central difference return? (numerical)
Treat loss as a function of $w$ alone, holding $x$ fixed:
At $w = 2$ with $h = 0.0001$:
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:
Numerical checking gives the same derivative from finite differences — say:
Absolute error:
Tiny. The backward implementation passes.
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:
For each parameter I need $L(\theta_i + h)$ and $L(\theta_i - h)$ — two forward passes per parameter.
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
original = w.data[i]
Step 2 — perturb up, run forward
w.data[i] = original + h
loss_plus = forward()
Step 3 — perturb down, run forward
w.data[i] = original - h
loss_minus = forward()
Step 4 — restore
w.data[i] = original
Step 5 — compute numerical gradient
grad_numeric = (loss_plus - loss_minus) / (2 * h)
Step 6 — compare
grad_auto = w.grad[i]
error = abs(grad_auto - grad_numeric)
For parameter $\theta_i$, the full comparison is:
Wrapped into a reusable function:
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
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:
With $\epsilon = 0.01$:
Absolute error:
Relative error (the comparison I actually use in gradcheck):
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 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.
Support the writing
If this post helped, a coffee keeps the deep dives, paper picks, and news digests coming.
Buy Me a Coffee