Checking a gradient with finite differences, and the step size that lies

The one test every hand-written gradient needs

Module 3 defined the derivative as the slope between two points as the gap shrinks, and noted that you can compute a usable one with two evaluations and a small number. That is also the way to check a gradient you derived by hand, or a custom backward function, or an autograd system you do not fully trust. Nudge one parameter, see how the loss moves, compare with what the gradient claimed:

numerical = ( L(θ + h) − L(θ − h) ) / (2h)

The two-sided form, with +h and −h, is worth the second evaluation: its error shrinks with h² rather than h, which the next section shows matters.

Why h cannot be tiny

The obvious instinct is to make h as small as possible, because the definition of the derivative is a limit. Two errors pull against each other.

Truncation error. The slope between two points is not the slope at one point; the curvature of the function contributes an error of order h² for the two-sided formula. Smaller h reduces it.

Rounding error. L(θ + h) − L(θ − h) is a subtraction of two nearly equal numbers, and the previous lessons showed what that does: the difference is known only to about machine epsilon ε times the size of L. Dividing by 2h amplifies that to ε/h. Smaller h increases it.

The total error is roughly h² + ε/h, smallest when the two terms balance, at h ≈ ε^(1/3):

float64: ε = 2.2e-16  →  best h ≈ 6e-6,   best achievable error ≈ 4e-11
float32: ε = 1.2e-7   →  best h ≈ 5e-3,   best achievable error ≈ 2e-5

The float32 line is the important one. In single precision the best possible numerical derivative agrees with the truth to about five digits, and with h = 1e-6, which is what everyone tries first, the rounding term is 0.1: the check is meaningless. Do gradient checks in float64. Convert the model, run the check, convert back.

What to compare, and the thresholds

A raw difference between the numerical and analytic gradients means nothing without a scale, so use the relative error:

rel = |analytic − numerical| / ( |analytic| + |numerical| + 1e-12 )

In float64 with h ≈ 1e-5:

rel < 1e-7    correct
rel ~ 1e-5    probably correct; possibly a kink nearby, see below
rel ~ 1e-3    something is slightly wrong: a missing factor, a transposed matrix
rel > 1e-2    the gradient is wrong

Check a random sample of parameters rather than all of them; a model with a million weights would need two million forward passes, and twenty randomly chosen weights find a systematic error just as well.

The kinks

A ReLU has no derivative at zero; it has a slope of 0 on one side and 1 on the other. If a parameter nudged by h moves some unit across zero, the numerical derivative averages the two slopes and disagrees with the analytic one, which picked a side. This is not a bug. It shows up as a handful of parameters with relative error near 0.5 while all the others are at 10^-8, and the tell is that the disagreeing ones change when you change h or the input. max, abs, sorting and any thresholding have the same kinks. If every parameter disagrees, the gradient is wrong; if a few do and the pattern moves, it is the kinks.

The code

import numpy as np
def grad_check(loss, theta, analytic, h=1e-5, n=20, rng=np.random.default_rng(0)):
    theta = theta.astype(np.float64)
    for i in rng.choice(theta.size, n, replace=False):
        e = np.zeros_like(theta); e.flat[i] = h
        num = (loss(theta + e) - loss(theta - e)) / (2 * h)
        ana = analytic.flat[i]
        rel = abs(ana - num) / (abs(ana) + abs(num) + 1e-12)
        print(i, f"{ana:.6e} {num:.6e} {rel:.1e}")

Run it once against a gradient you are sure of, so you know what "passes" looks like on your machine, before running it against the one you doubt.

When the numerical derivative is the derivative

Sometimes there is no analytic gradient: a loss that calls a simulator, a black-box metric, a function with a lookup table in it. Finite differences then are the gradient, at a cost of two function evaluations per parameter, which is affordable for tens of parameters and hopeless for millions. That is the honest statement of why backpropagation exists: it computes every parameter's gradient for the cost of about two forward passes in total, where finite differences cost two per parameter. Module 3 made that point from the algorithm's side. From the numerical side, the finite-difference gradient is also less accurate, by the ε^(2/3) above, so backpropagation is both a million times cheaper and a hundred thousand times more precise. Use finite differences to check it, in float64, at the h the arithmetic prescribes, and for nothing else.