Gradient Checking by Finite Differences

An analytically derived gradient can be wrong in ways that raise no error and still train a model, slowly and wrongly: a transposed weight gradient in a square layer, a sign error that a symmetric initialization hides, a missing fan-out accumulation of the kind described under reverse-mode differentiation, a reduction factor applied to the biases and not the weights. Comparing the analytic derivative against a numerical estimate of the same derivative is the oldest and most direct way to catch all of these, and its accuracy is bounded by two competing errors that determine the step size, the comparison metric, and the precision the check has to run in.

Truncation against round-off

The forward difference $\big(f(x+h) - f(x)\big)/h$ equals $f’(x) + \frac{h}{2}f”(x) + O(h^{2})$: the neglected terms are linear in the step, so the truncation error falls only as $h$ shrinks. The central difference $\big(f(x+h) - f(x-h)\big)/2h$ equals $f’(x) + \frac{h^{2}}{6}f'''(x) + O(h^{4})$, with the first-order term cancelling between the two evaluations; its truncation error falls quadratically, at the price of exactly twice the function evaluations.

Round-off works in the opposite direction. Evaluating $f$ at a point returns a value with a relative error of order $\epsilon$, the unit round-off of the format, so the two function values in the central difference carry absolute errors of order $\epsilon|f|$. Dividing their difference by $2h$ divides that absolute error by $2h$ while the numerator itself is about $2hf’(x)$, so the relative error contributed by cancellation behaves like $\epsilon|f|/(h|f’|)$. Summing the two contributions gives an error of the form $c_1 h^{2} + c_2 \epsilon/h$, minimized where the two terms balance, at a step of order $\epsilon^{1/3}$ for the central difference. In 64-bit arithmetic, with $\epsilon = 2.2\times10^{-16}$, that step is around $6\times10^{-6}$ and the achievable relative accuracy is a few times $10^{-11}$. In 32-bit arithmetic, with $\epsilon = 1.2\times10^{-7}$, the optimal step is around $5\times10^{-4}$ and the achievable accuracy is a few times $10^{-5}$. For the forward difference the optimum is at $\epsilon^{1/2}$, and the best possible agreement is around $10^{-8}$ in double precision.

The practical consequence is that a check in 32-bit arithmetic cannot distinguish a correct gradient from one with an error of $10^{-4}$ relative, which is why gradient checks are run in double precision with every parameter and every intermediate promoted. A network whose normal training precision is 32-bit has to be able to run an evaluation pass in 64-bit for this purpose.

The metric and the step, chosen together

The comparison should be scale-free and symmetric, so that a gradient entry of $10^{-12}$ does not dominate the statistic. The standard form is

$$\text{rel} = \frac{|g_{\text{analytic}} - g_{\text{numeric}}|}{\max(|g_{\text{analytic}}| + |g_{\text{numeric}}|,\ \delta)}$$

with a floor $\delta$ that prevents a division by a vanishing denominator from turning an absolute error of $10^{-16}$ into a huge ratio. With the floor at $10^{-8}$ and both gradients of order $10^{-3}$, an absolute disagreement of $10^{-7}$ gives a relative error of about $5\times10^{-5}$, which is an order of magnitude above what a correct double-precision check produces and therefore a genuine suspicion rather than noise. The commonly quoted guidance is that failures below about $10^{-7}$ are consistent with a correct implementation in double precision, values near $10^{-4}$ warrant inspection, and values above $10^{-3}$ indicate a defect, but the threshold depends on the precision, on whether the loss is a sum over a large batch, and on how close the parameters sit to a kink.

The step must also be scaled to the parameter. Adding $h$ to a parameter whose magnitude is $10^{-3}$ is a large perturbation of that parameter, and for a piecewise-linear network it may cross a kink that sits between the two evaluation points, in which case the central difference returns the average of the two one-sided slopes rather than any derivative of the function at $x$. A relative step, applied to each parameter in proportion to its magnitude, keeps the perturbation small in the same sense for all of them.

Every parameter, or one random direction

Central differences cost two evaluations of the loss per parameter, so checking a model with $10^{6}$ parameters exhaustively costs two million forward passes. The cheap alternative is to project onto a random direction: with a random vector $u$ of unit norm, the analytic quantity $\sum_i u_i g_i$ is one vector inner product, and the numerical quantity is $(L(\theta + hu) - L(\theta - hu))/2h$, two evaluations regardless of the parameter count.

The projection’s weakness is that it turns a vector comparison into a scalar one. An error vector $e$ contributes to the comparison in proportion to its inner product with $u$, so a sparse error confined to one coordinate of a million-parameter model produces a signal smaller than the norm of the error by a factor near $10^{3}$, which hides it inside the noise floor of the check. Inverting the argument gives the useful diagnosis: a projection check that fails is conclusive, since the error must have a component along the direction, while a projection check that passes only rules out errors aligned with that particular direction. Checking one layer’s parameters exhaustively, where the count is small enough, combines the two: it is affordable for a layer with a few thousand parameters and it cannot miss a sparse error inside that block. A single pass over the whole model is best treated as a smoke test.

Where the check reports errors that are not errors

A gradient check compares the analytic derivative against the derivative of the function that was actually evaluated, and that symmetry is also its blind spot. If the forward pass is wrong, both sides use the same wrong function and agree. If the loss reduction sums a batch where the analytic gradient assumes a mean, both sides see the same scaling as long as the same code path is used for both. A mislabeled target, a data pipeline that shuffles labels inconsistently, or a loss that ignores part of its input are all invisible.

Several features of real networks make the numerical side inaccurate rather than the analytic one. Kernels at exactly zero sit at the kink of a rectifier, where the derivative is undefined and the two evaluations straddle it; the check will disagree with the analytic subgradient convention by a large relative amount at those parameters. Pooling and absolute value have the same structure. A masking layer that redraws its mask on each call makes the two evaluation points see different functions, so the comparison measures the masks; every stochastic layer must be disabled, or driven from a fixed seed shared by all three evaluations. And a network whose parameters are all zero produces zero gradients everywhere, which agrees trivially and verifies nothing, so the check should be run at several randomly initialized parameter settings rather than one.

Two stronger tests that cost less

The complex-step derivative avoids the cancellation term entirely. If $f$ is analytic, then $\operatorname{Im}\big(f(x + ih)\big)/h$ equals $f’(x) + O(h^{2})$ with no subtraction of nearly equal quantities in the numerator, so no round-off amplification occurs and the step can be taken as small as $10^{-20}$ with the truncation error going to zero. The price is that the function must be evaluated with complex arguments, which excludes every operation that is not analytic: comparisons and branches, maximum and absolute value, and the rectifier itself. It is nonetheless the method of choice for differentiable primitives, and it is the reason complex-step differentiation persists in sensitivity analysis of simulation codes. The one implementation trap is that the imaginary part must be read from the result rather than from a norm or a squared magnitude, since those reintroduce conjugation and destroy the imaginary component carrying the derivative.

The second test needs no step size at all. For any primitive with both a forward and a reverse rule, the adjoint identity $\langle \bar{v}, J\dot{x}\rangle = \langle J^{\top}\bar{v}, \dot{x}\rangle$ must hold for random directions, and both sides are computed from the rules themselves rather than from a differencing formula. It catches transposition and shape errors at machine precision in two evaluations per random pair, at the cost of requiring a working implementation of both products, which the machinery described in the products section provides.

import numpy as np

def numeric_grad(f, x, h, metric_floor=1e-8):
    g = np.zeros_like(x)
    for i in range(x.size):
        step = h * max(abs(x.flat[i]), 1.0)      # relative step per parameter
        xp, xm = x.copy(), x.copy()
        xp.flat[i] += step
        xm.flat[i] -= step
        g.flat[i] = (f(xp) - f(xm)) / (2.0 * step)
    return g

def rel_error(analytic, numeric, floor=1e-8):
    denom = np.maximum(np.abs(analytic) + np.abs(numeric), floor)
    return np.max(np.abs(analytic - numeric) / denom)

f = lambda t: np.sum(3.0 * t ** 2 + t)           # d/dt = 6 t + 1
x = np.array([0.5, -1.25])
analytic = 6.0 * x + 1.0
numeric = numeric_grad(f, x, h=6e-6)             # ~ eps**(1/3) in float64
assert rel_error(analytic, numeric) < 1e-9

The check also needs its own negative control. Subtracting a deliberate constant from one entry of the analytic gradient should push the relative error into the region the workflow treats as failure; a check whose thresholds cannot be broken by a known error is not measuring anything. That control, run once on a small model, establishes what a real defect looks like in the particular metric and precision being used, which is the only way to read the numbers in between.