A gradient checker you can trust

From week 1 on you will derive every gradient by hand. This unit builds the instrument that tells you whether a derivation is right: compare it with a brute-force estimate that needs no calculus at all. If they agree to about seven digits, the derivation is almost certainly correct.

Needs: single-variable derivatives and Taylor’s theorem. These are Math course units C2 and C5 on the syllabus, not yet written; until then, Deisenroth, Faisal & Ong, Mathematics for Machine Learning, §5.1. Used by: every week that has a gradient in it.

The estimate

Definition W0.1Finite-difference estimates#

For $f : \mathbb{R} \to \mathbb{R}$ and a step $h > 0$, the forward and centered differences are

$$D^{+}_h f(x) := \frac{f(x + h) - f(x)}{h}, \qquad D^{0}_h f(x) := \frac{f(x + h) - f(x - h)}{2h}.$$

For $f : \mathbb{R}^n \to \mathbb{R}$, estimate each partial derivative the same way along a coordinate direction $e_i$ (the vector with $1$ in position $i$, $0$ elsewhere):

$$\big(\hat\nabla_h f(x)\big)_i := \frac{f(x + h e_i) - f(x - h e_i)}{2h}, \qquad i = 1, \dots, n.$$

Both are slopes of a chord. As $h \to 0$ they tend to $f'(x)$, if $f$ is differentiable at $x$. The centered one looks both ways and is far more accurate; for a full gradient it costs $2n$ evaluations of $f$ instead of $n + 1$, since the forward difference can share $f(x)$ across coordinates.

Theorem W0.2Truncation error#

Let $h > 0$ and let $f$ be defined on an open interval containing $[x - h, x + h]$.

1. If $f''$ exists and is continuous there, then $D^{+}_h f(x) - f'(x) = \tfrac{h}{2} f''(\xi)$ for some $\xi \in (x, x + h)$.
2. If $f'''$ exists and is continuous there, then $D^{0}_h f(x) - f'(x) = \tfrac{h^2}{6} f'''(\eta)$ for some $\eta \in (x - h, x + h)$.

So the forward error shrinks like $h$ and the centered error like $h^2$: halve the step and you halve the first error but quarter the second.

Proofof W0.2

Both parts use Taylor’s theorem with the Lagrange form of the remainder.

1. $f(x + h) = f(x) + h f'(x) + \tfrac{h^2}{2} f''(\xi)$ for some $\xi \in (x, x + h)$. Subtract $f(x)$, divide by $h$, and subtract $f'(x)$.
2. Expand one order further on each side: $f(x \pm h) = f(x) \pm h f'(x) + \tfrac{h^2}{2} f''(x) \pm \tfrac{h^3}{6} f'''(\xi_\pm)$, with $\xi_+ \in (x, x + h)$ and $\xi_- \in (x - h, x)$. Subtracting, the $f(x)$ and $f''(x)$ terms cancel: $f(x + h) - f(x - h) = 2h f'(x) + \tfrac{h^3}{6}\big(f'''(\xi_+) + f'''(\xi_-)\big)$. Divide by $2h$: $D^{0}_h f(x) - f'(x) = \tfrac{h^2}{12}\big(f'''(\xi_+) + f'''(\xi_-)\big)$. The average $\tfrac12\big(f'''(\xi_+) + f'''(\xi_-)\big)$ lies between the two values, so since $f'''$ is continuous, the intermediate value theorem gives $\eta$ between $\xi_-$ and $\xi_+$ where $f'''$ equals it. That gives $\tfrac{h^2}{6} f'''(\eta)$.

Full explanationWhy the centered difference is so much better

The cancellation in step 2 is the whole story. The forward difference compares $x$ with $x + h$, so the curvature $f''$ tilts the chord and the error is first order. The centered chord is symmetric around $x$: curvature lifts both ends equally, and its effect on the slope cancels. What survives is the next term, the third derivative, multiplied by $h^2$.

The same symmetry trick appears all over numerical computing: the trapezoid and midpoint rules for integrals, the leapfrog integrator used in physics simulations and in Hamiltonian Monte Carlo, and Richardson extrapolation builds on it, combining two step sizes to cancel the next error term as well (it works especially well with symmetric formulas, whose errors contain only even powers of $h$).

The second error: rounding

Theorem W0.2 says “make $h$ small”. A computer disagrees.

Remark W0.3Rounding error, and the best step#

A float64 number carries about 16 significant digits: its machine epsilon is $\varepsilon \approx 2.2 \times 10^{-16}$. Model each computed value of $f$ as $f(1 + \delta)$ with $|\delta|$ at most a small multiple of $\varepsilon$ (a well-implemented library function is accurate to a few units in the last place; forming $x \pm h$ adds a rounding term of the same order, which we absorb), and let $F$ bound $|f|$ near $x$. The numerator $f(x + h) - f(x - h)$ can then be off by up to $2\varepsilon F$, and dividing by $2h$ turns that into an error up to $\varepsilon F / h$, which grows as $h$ shrinks. With $M_3$ bounding $|f'''|$, the total error of the centered difference is at most

$$E(h) = \frac{M_3}{6} h^2 + \frac{\varepsilon F}{h}, \qquad \text{minimised at} \qquad h^{*} = \left(\frac{3 \varepsilon F}{M_3}\right)^{1/3}.$$

For the forward difference the same reasoning gives $\tfrac{M_2}{2} h + \tfrac{2 \varepsilon F}{h}$, minimised at $h^{*} = 2\sqrt{\varepsilon F / M_2}$. When $F$ and the derivative bounds are of order $1$, that is $h^{*} \approx 10^{-5}$ for centered and $h^{*} \approx 3 \times 10^{-8}$ for forward.

These are worst-case bounds; real rounding errors partly cancel, so the observed errors are usually smaller. The best step still lands in the predicted decade. Before reading on, test that:

Predict first. For f = sin at x = 1 in float64, roughly which step size h gives each estimate its smallest error?

Show the real curves

Read the slopes (the axis runs from large $h$ on the left to small $h$ on the right). Left of each dip, the centered curve falls twice as steeply as the forward one: on log-log axes the truncation error has slope of magnitude $2$ against $1$, which is Theorem W0.2. Right of the dips, both rise together with slope of magnitude $1$, because rounding error does not care which formula you used. This is why CS231n recommends $h \approx 10^{-5}$; PyTorch’s gradcheck defaults to $10^{-6}$.

Full explanationWhy float32 is hopeless here

Repeat the analysis with float32, where $\varepsilon \approx 1.2 \times 10^{-7}$. The best centered step becomes about $(3\varepsilon)^{1/3} \approx 7 \times 10^{-3}$, and the error bound at that step is about $\varepsilon^{2/3} \approx 2 \times 10^{-5}$ relative to the scale of $f$, before any of the extra rounding inside a real model. (For $\sin$ at $1$ the observed float32 minimum is smaller, about $6 \times 10^{-7}$, but at a step in the predicted decade.) CS231n’s authors report having sometimes seen relative errors drop from about $10^{-2}$ to $10^{-8}$ just by switching a correct implementation from single to double precision. A gradient check is a float64 tool.

Judging the agreement

Definition W0.4Relative error#

For an analytic gradient $a$ and a numerical one $b$ of the same shape, and a small floor $\epsilon_0 > 0$,

$$\operatorname{rel}(a, b) := \max_i \frac{|a_i - b_i|}{\max\big(|a_i|, |b_i|, \epsilon_0\big)}.$$

An absolute difference of $10^{-4}$ means nothing on its own: tiny next to a gradient of $1000$, huge next to one of $10^{-5}$. Dividing by the larger magnitude makes the measure scale-free and symmetric in $a$ and $b$. The floor $\epsilon_0$ (the stub uses $10^{-12}$) only prevents $0/0$ when both entries are exactly zero; it does not rescue entries whose true value is near zero (failure mode 5 in W0.6). Read the result against the CS231n scale, for a single smooth function in float64:

$\operatorname{rel}$	Verdict
$> 10^{-2}$	The gradient is probably wrong.
$10^{-4}$ to $10^{-2}$	Suspicious. Investigate before moving on.
$10^{-7}$ to $10^{-4}$	Acceptable only if the function has kinks (ReLU, max); too large for smooth ones.
$< 10^{-7}$	Good.

Errors accumulate with depth, so a check on a whole deep network can be looser than this.

Algorithm W0.5Gradient check#

Input: $f : \mathbb{R}^n \to \mathbb{R}$, a claimed gradient function $g$, a test point $x$, step $h = 10^{-5}$, tolerance $\tau = 10^{-7}$.

1. Convert $x$ to a float64 copy.
2. For each coordinate $i$: set $x_i \leftarrow x_i + h$, evaluate $f$; set $x_i \leftarrow x_i - h$ (from the original), evaluate $f$; restore $x_i$; store the centered difference.
3. Return $\operatorname{rel}\big(g(x), \hat\nabla_h f(x)\big)$ and whether it is below $\tau$.

Cost: $2n$ evaluations of $f$ and one of $g$. That is why checks run on small inputs: a handful of data points and small layer sizes.

Where it fails

Failure mode W0.6Six ways a gradient check misleads you#

1. Kinks. At points where $f$ is not differentiable (ReLU at $0$, $\max$, $|\cdot|$), the two evaluations can land on different sides of the kink, and the estimate averages two slopes. Symptom: a few entries with large error, at inputs within $h$ of a kink. Fix: use test points away from kinks, few data points, or track whether a kink was crossed.
2. Single precision. A correct gradient checked in float32 can show errors near $10^{-2}$. Fix: float64 throughout, including the copy of $x$.
3. Randomness inside $f$. Dropout, sampling or data augmentation make $f(x + h e_i)$ and $f(x - h e_i)$ use different random draws, so $f$ is not a deterministic function of $x$. Fix: fix the seed or switch them off during the check.
4. The wrong step. Too small and rounding dominates; too large and truncation does (the plot above). Fix: $h = 10^{-5}$; if a check fails, retry with $10^{-4}$ and $10^{-6}$ before blaming the gradient.
5. Negligible gradients. Where the true derivative is zero or tiny, the numerical estimate is pure noise of size about $\varepsilon F / h + M_3 h^2 / 6$, and the relative error is large and meaningless. For $g(x) = (x - 0.6)^2 (x + 2)$ at $x = 0.6$ the derivative is exactly $0$, the centered estimate is about $10^{-10}$, and the relative error is $1$. Fix: an absolute tolerance alongside the relative one, set above that noise level (this is what PyTorch’s atol is for). A floor $\epsilon_0$ only helps if it, too, sits above the noise.
6. An unlucky test point. A check passes at one $x$, not everywhere. At $x = 0$ the buggy gradient $x$ and the correct gradient $2x$ of $\|x\|^2$ agree perfectly. With square shapes, a transposed gradient has the right shape and only its values give it away; with all-equal entries, even the values can agree. A batch of one hides mistakes in summing or averaging over the batch. Fix: random inputs, non-square shapes and a batch of at least two.

Where it fits

Full explanationThree ways to get a derivative, and what the check is checking

There are three ways to differentiate a program:

Method	How	Cost for a full gradient of $f : \mathbb{R}^n \to \mathbb{R}$	Accuracy
Symbolic	Apply differentiation rules to a formula	Can blow up in expression size	Exact
Numerical (this unit)	Finite differences	$2n$ evaluations of $f$	Limited by truncation and rounding
Automatic, reverse mode (backprop)	Record the computation, apply the chain rule backwards	Time: a small constant multiple of one evaluation, for any $n$. Memory: proportional to the stored intermediate values	Exact up to floating point

Reverse-mode automatic differentiation is why deep learning is possible: a model with a billion parameters gets its full gradient for roughly the price of a few forward passes. Finite differences would need two billion. So finite differences never train anything; their job is to check the other two, because they are too simple to be wrong in the same way. You will build reverse-mode differentiation yourself in Arc 2.

A cheaper check. Instead of all $n$ coordinates, pick one random direction $u$ and compare the claimed directional derivative $g(x)^\top u$ with $\big(f(x + hu) - f(x - hu)\big)/2h$. That costs two evaluations whatever $n$ is. If $g(x)$ is wrong, the error vector $g(x) - \nabla f(x)$ is non-zero, and a random direction drawn with a density (Gaussian, say, or uniform on the sphere) is exactly orthogonal to it with probability zero. With a tolerance the picture is less clean: a small error in one entry of a large gradient is diluted by roughly $\sqrt{n}$ in the projection and can slip under the threshold. This is what JAX’s check_grads does. The coordinate-wise check, by contrast, tells you which entries are wrong.

In production

Gradient checks are a testing tool, not a training one. Frameworks ship them for people writing new operations: PyTorch’s torch.autograd.gradcheck compares analytical gradients against small finite differences, with default tolerances tuned for float64 inputs; JAX’s jax.test_util.check_grads compares forward- and reverse-mode derivatives against finite differences along a random direction. They run in unit tests and CI on tiny inputs, when someone writes a custom layer, a custom backward pass or a hand-written GPU kernel. PyTorch tests closeness with allclose, $|a - b| \le \text{atol} + \text{rtol} \cdot |b|$ with $b$ the numerical gradient, rather than the single symmetric relative error of W0.4; its fast_mode option also checks random projections, as JAX does. Different conventions for the same question.

Check yourself

Check Why does the truncation error of the centered difference shrink like h² rather than h?

Subtracting the expansions of $f(x + h)$ and $f(x - h)$ cancels the even-power terms, including the $\tfrac{h^2}{2} f''$ term that limits the forward difference. The first survivor is $\tfrac{h^3}{6} f'''$, which becomes an $h^2$ error after dividing by $2h$ (W0.2).

Check In float64 you shrink h from 1e-5 to 1e-12. What happens to the centered estimate?

It gets worse. Rounding error grows like $\varepsilon F / h$ (W0.3); in the plot, the error at $h = 10^{-12}$ is several orders of magnitude larger than at $10^{-5}$.

Check A smooth loss's analytic gradient has relative error 3e-3 against the numerical one in float64. What do you do?

Suspect a bug in the analytic gradient first. For a smooth function in float64 you expect about $10^{-7}$ or better; $3 \times 10^{-3}$ is in the suspicious band. Retry with $h = 10^{-4}$ and $10^{-6}$ to rule out the step.

Check Why divide by max(|a|, |b|) rather than by |a| alone?

The result does not depend on which gradient you call the reference, and a zero in one of them (common with ReLU) does not cause a division by zero.

Check The analytic gradient of f(x) = ‖x‖² is 2x. A buggy version returns x. At which test point does a gradient check fail to notice?

At $x = 0$, where both are the zero vector. Any point with a non-zero entry exposes the bug. This is failure mode 6 in W0.6.

Check With a random direction u, why does the directional check almost never miss a wrong gradient?

In exact arithmetic, the check misses only if $(g(x) - \nabla f(x))^\top u = 0$, that is if $u$ is orthogonal to a non-zero error vector. The set of such directions is a hyperplane, which a random $u$ drawn with a density hits with probability zero. In floating point with a tolerance, small errors can still be missed (see Where it fits).

Build it

Build mlkit.gradcheck

• mlkit/gradcheck.py
• tests/test_gradcheck.py

Implement the three functions in mlkit/gradcheck.py with NumPy only: numerical_grad (W0.1, centered, any array shape), relative_error (W0.4, vectorised, no loop) and check_grad (W0.5). The seven tests in tests/test_gradcheck.py are the specification. One of them, a preview of week 1, checks the gradient $2(Wv - y)v^\top$ of $f(W) = \|Wv - y\|^2$ with respect to a matrix.

Done when this passes: uv run pytest

Hint Iterating over an array of any shape

np.ndindex(x.shape) yields every index tuple of an array of any shape, so one loop handles vectors and matrices alike.

Hint The two tests that catch most people

test_integer_input: adding $10^{-5}$ to an integer array silently does nothing. test_does_not_mutate_input: the caller’s array must come back unchanged. One line at the top of numerical_grad fixes both: work on a float64 copy.

Break it. In weeks/w00-setup/error_vs_h.py, redraw the error curves with Matplotlib for $f(x) = e^x$ at $x = 1$ (log-log axes, $h$ from $10^{-1}$ to $10^{-15}$), and save error_vs_h.png. Then run check_grad on $f(x) = \sum_i \max(x_i, 0)$, with analytic gradient $1$ where $x_i > 0$ and $0$ elsewhere, at a point where one entry of $x$ equals $10^{-7}$. Predict the relative error first, then explain the result with W0.6.

Check Do the minima for eˣ at x = 1 sit at different step sizes than for sin at x = 1?

Barely. By W0.3, the centered $h^{*}$ depends on $F / M_3$ through a cube root. For $e^x$ at $1$, $F = M_3 = e$, a ratio of $1$; for $\sin$ at $1$ the ratio is $\sin 1 / \cos 1 \approx 1.56$, whose cube root is about $1.16$: a shift of about $0.06$ of a decade, which the noisy curves cannot resolve. The forward $h^{*}$ depends on $F / M_2$, which is $1$ in both cases ($\sin 1 / \sin 1$ and $e / e$), so the prediction does not move at all.

Derive it. Write the proof of W0.2 on paper before opening it, then compare. If you stall, note the step in math-debt.md.

Stretch. Add directional_check(f, grad_f, x, h=1e-5, rng=None) using one random direction, with a test that it catches the buggy gradient $x$ for $\|x\|^2$.

Then commit, push, and ask for a review.

Sources

• Stanford CS231n, “Neural Networks Part 3: Gradient checks”: centered differences, relative error, the verdict scale, step size, precision, kinks.
• PyTorch documentation, torch.autograd.gradcheck.
• JAX documentation, jax.test_util.check_grads.
• Baydin, Pearlmutter, Radul & Siskind, “Automatic differentiation in machine learning: a survey”, Journal of Machine Learning Research 18 (2018): symbolic, numerical and automatic differentiation compared.
• Deisenroth, Faisal & Ong, Mathematics for Machine Learning, §5.1: derivatives and Taylor series.