Numerical differentiation and gradient checking

Before this lesson: lesson 03 (Taylor's theorem with remainder, the Taylor test), lesson 06 (the gradients you will check), and floating point and cancellation from linear algebra lesson 13. Time: about 60 minutes, plus the notebook.

---

Why this matters #

A hand-written backward pass, a custom CUDA or Triton kernel, a new loss term, the adjoint of a PDE solver: each is a gradient that can be wrong without any error message. Training still runs; it is only slower or worse. A finite-difference check needs nothing but evaluations of $f$, so it is the standard unit test for gradients (torch.autograd.gradcheck, jax.test_util.check_grads).

It is also a small, complete lesson in floating point: a formula that is exact as $h \to 0$ gives garbage at $h = 10^{-12}$. In science codes, the same ideas appear as the gradient test of data-assimilation systems. The complex step is a standard way to get derivatives from simulation code that has no automatic differentiation.

1. Predict first #

2. The picture #

Secant slopes. The derivative $f'(x)$ is the slope of the tangent. A finite difference uses the slope of a secant through two nearby points:

• the forward difference $D_+(h) = \dfrac{f(x+h) - f(x)}{h}$ uses the points $x$ and $x + h$;
• the central difference $D_0(h) = \dfrac{f(x+h) - f(x-h)}{2h}$ uses the points $x - h$ and $x + h$.

The figure shows both secants next to the tangent.

Interactive figure: loading…

The central secant is better because it is symmetric. The curvature makes the forward secant too steep on one side; with points on both sides of $x$, the two curvature errors cancel. Only an $h^2$ term is left (Theorem 1 below).

Two errors with opposite trends.

• Truncation error. The secant is not the tangent. Its slope differs from $f'(x)$ by a term that shrinks like $h$ (forward) or $h^2$ (central). This is pure math: it exists even with exact arithmetic.
• Rounding error. In float64 every computed value of $f$ has a relative error of about $\varepsilon \approx 2.2\cdot 10^{-16}$ (machine epsilon). For small $h$, $f(x+h)$ and $f(x)$ agree in their leading digits. The subtraction cancels those digits, and what is left is mostly rounding noise. This is catastrophic cancellation (linear algebra lesson 13). Dividing by a small $h$ makes the noise large: this error grows like $1/h$.

On a log–log plot, the total error against $h$ is a V. The right arm is truncation, with slope $+1$ (forward) or $+2$ (central). The left arm is rounding, with slope $-1$. The best step sits at the bottom of the V, far from "as small as possible".

Interactive figure: loading…

The complex step leaves the real axis. It evaluates $f(x + ih)$, with $i = \sqrt{-1}$, and keeps the imaginary part. For a real-analytic $f$, $f(x + ih) \approx f(x) + ih\,f'(x)$: the derivative appears in the imaginary part, and nothing is subtracted. So there is no cancellation and no left arm.

In $n$ dimensions a finite-difference gradient moves one coordinate at a time: $n$ or $2n$ evaluations of $f$. Reverse-mode automatic differentiation gives the whole gradient for a small constant times one evaluation (lesson 05). So finite differences are a check, never a way to train.

3. Definitions and main results #

Notation. $h > 0$ is the step. $\varepsilon = 2^{-52} \approx 2.2\cdot10^{-16}$ is float64 machine epsilon. $f$ is a real function of a real $x$, smooth near $x$. We use $\varepsilon$ as the size of the relative error of one computed value of $f$. (One correctly rounded operation has error at most $u = \varepsilon/2$, the unit roundoff; a computed $f$ has several roundings, so $\varepsilon$ is a fair round number.)

Theorem 1 (truncation error) #

If $f \in C^2$, then $D_+(h) = f'(x) + \tfrac{h}{2} f''(\xi)$ for some $\xi \in (x, x+h)$. If $f \in C^3$, then $D_0(h) = f'(x) + \tfrac{h^2}{6} f'''(\eta)$ for some $\eta \in (x-h, x+h)$.

So the forward difference is first order, $O(h)$, and the central difference is second order, $O(h^2)$.

Proof.

1. Taylor with the Lagrange remainder (lesson 03): $f(x+h) = f(x) + h f'(x) + \tfrac{h^2}{2} f''(\xi)$.
2. Subtract $f(x)$ and divide by $h$: this is the forward statement.
3. One order more, on both sides:

4. Subtract the two lines. The even terms $f(x)$ and $\tfrac{h^2}{2} f''(x)$ cancel:

5. Divide by $2h$. The average of $f'''(\eta_+)$ and $f'''(\eta_-)$ lies between the two values, so it equals $f'''(\eta)$ for some $\eta$ (intermediate value theorem). $\blacksquare$

The same expansions give the second difference for the second derivative:

Richardson extrapolation

If $f$ is smooth, only even powers of $h$ appear in the central difference: $D_0(h) = f'(x) + c\,h^2 + O(h^4)$. Halving the step gives $D_0(h/2) = f'(x) + c\,h^2/4 + O(h^4)$. Take the combination that cancels the $h^2$ term:

$$R(h) = \frac{4D_0(h/2) - D_0(h)}{3} = f'(x) + O(h^4).$$

This is Richardson extrapolation. It turns two second-order values into one fourth-order value. It is the idea behind Romberg integration and many high-order finite-difference stencils. The rounding arm becomes worse, so the best step is larger ($\propto \varepsilon^{1/5}$).

Theorem 2 (the optimal step) #

Error model. Each computed value $\tilde f$ satisfies $|\tilde f - f| \le \varepsilon F$, with $F \approx |f|$ near $x$. Near $x$, $|f''| \le M_2$ and $|f'''| \le M_3$. Then the total errors are bounded by

and the bounds are smallest at

Why this is true (intuition). One term falls and one term rises as $h$ shrinks. The sum is smallest where the two terms are of the same size: the bottom of the V.

Proof.

1. Forward: the computed numerator is $f(x+h) - f(x) + \delta_1 - \delta_2$ with $|\delta_i| \le \varepsilon F$. After dividing by $h$, the rounding part is at most $2\varepsilon F/h$. Add the truncation bound $M_2 h/2$ from Theorem 1.
2. Set the derivative to 0: $\tfrac{d}{dh}E_+ = M_2/2 - 2\varepsilon F/h^2 = 0$ gives $h_+^*$. There the two terms are equal, each $\sqrt{\varepsilon F M_2}$, so $E_+ = 2\sqrt{\varepsilon F M_2}$.
3. Central: the rounding part is at most $2\varepsilon F/(2h) = \varepsilon F/h$; truncation at most $M_3 h^2/6$.
4. $\tfrac{d}{dh}E_0 = M_3 h/3 - \varepsilon F/h^2 = 0$ gives $h_0^*$. There $\varepsilon F/h_0^* = M_3 h_0^{*2}/3$, so $E_0 = M_3 h_0^{*2}/6 + M_3 h_0^{*2}/3 = M_3 h_0^{*2}/2$, which is the formula. $\blacksquare$

Numbers (float64, $F = M_2 = M_3 = 1$).

	best step $h^*$	best error	correct digits
forward	$2\sqrt\varepsilon \approx 3.0\cdot10^{-8}$	$\approx 3.0\cdot10^{-8}$	about 8 (half)
central	$(3\varepsilon)^{1/3} \approx 8.7\cdot10^{-6}$	$\approx 3.8\cdot10^{-11}$	about 10–11 (two thirds)

In general, a formula of order $m$ has $h^* \propto \varepsilon^{1/(m+1)}$ and best error $\propto \varepsilon^{m/(m+1)}$. These are worst-case bounds; measured errors are often a bit smaller.

Scale. If $f$ changes on the length scale $|x|$ (for example $x^k$ or $\log x$), then $F/M_2 \sim x^2$ and $h^* \propto \sqrt{\varepsilon}\,|x|$. So use a relative step $h_i = c \cdot \max(1, |x_i|)$, with $c = \sqrt\varepsilon \approx 1.5\cdot10^{-8}$ (forward) or $c = \varepsilon^{1/3} \approx 6\cdot10^{-6}$ (central). A fixed $h = 10^{-6}$ at $x_i = 10^6$ changes $x_i$ only in about its 13th digit.

Theorem 3 (the complex step) #

Let $f$ be analytic near $x$ and real on the real axis, so that $f(x)$ and all its derivatives at a real $x$ are real. Then

Proof.

1. Taylor series in the direction $ih$, with $i^2 = -1$ and $i^3 = -i$:

2. All derivatives at $x$ are real. So the imaginary part is $h f'(x) - \tfrac{h^3}{6} f'''(x) + O(h^5)$, and the real part is $f(x) - \tfrac{h^2}{2} f''(x) + O(h^4)$.
3. Divide the imaginary part by $h$. $\blacksquare$

Why $h = 10^{-20}$ works. No two nearly equal numbers are subtracted. The imaginary part, about $h f'(x)$, is computed with a relative error of about $\varepsilon$ for any $h$ down to underflow. The truncation term $h^2 f'''/6$ is about $10^{-40}$. So the result is accurate to machine precision, from one complex evaluation per directional derivative. A complex evaluation costs a few real ones (a complex multiply is 4 real multiplies and 2 adds).

What the code must satisfy. $f$ must be analytic, and the code must be complex-safe:

• no abs: it returns the modulus, a real number, so the derivative of that term is lost. Use np.where(z.real >= 0, z, -z) instead;
• comparisons and max only on the real part;
• no complex conjugation (np.vdot, .conj(), conjugate transposes);
• no cast to real, and complex dtypes everywhere.

The complex step is close to forward-mode automatic differentiation with dual numbers (lesson 05): dual numbers use $\epsilon^2 = 0$, the complex step uses $i^2 = -1$, and the difference is the $O(h^2)$ term.

Gradients, Jacobians and directional checks #

For $f: \R^n \to \R$, the central-difference gradient is

where $e_i$ is the $i$-th unit vector: $2n$ evaluations of $f$. For $F: \R^n \to \R^m$, the same difference of output vectors gives column $i$ of the $m \times n$ Jacobian.

A cheaper test uses one random direction $v$: the directional check

two evaluations for the whole gradient. Compare it with the claimed $g^\top v$. A few random directions catch almost every bug.

The relative error #

Compare the claimed gradient $g$ with $g_{\text{num}}$ by the relative error

Gradient sizes differ by orders of magnitude between layers. An absolute difference of $10^{-6}$ is a bug when $\|g\| = 10^{-8}$ and noise when $\|g\| = 10^4$. The relative error does not depend on the scale. It lies in $[0, 2]$. (Some libraries divide by $\|g\| + \|g_{\text{num}}\|$; the result differs at most by a factor of 2.)

Thresholds for a smooth function in float64 with a central difference: $\text{rel} < 10^{-7}$ is fine; $\text{rel} > 10^{-3}$ is almost surely a bug; in between, look closer (kinks, scale, step).

A reliable gradient check #

1. float64. In float32 ($\varepsilon \approx 1.2\cdot10^{-7}$) the best central error is about $2.5\cdot10^{-5}$: too coarse to separate a small bug from rounding.
2. Central differences with $h \approx \varepsilon^{1/3}\max(1, |x_i|)$, about $10^{-6}$ to $10^{-5}$ times the scale. If you are unsure, try two or three steps: a real bug does not depend on $h$.
3. Relative error with the thresholds above.
4. No kinks. ReLU, max, abs, hinge losses, sorting and top-$k$ are not differentiable everywhere. If $x \pm h$ crosses a kink, the numerical derivative is wrong (figure below).
5. Determinism. Turn off dropout and data augmentation, fix seeds, put batch normalization in a fixed mode: $f$ must be the same function at every evaluation.
6. The same $f$ on both sides. Include every loss term (regularization, auxiliary losses) in both $f$ and $g$.
7. A generic point. Use random inputs and weights, not zeros or a symmetric initialization, where many terms vanish and bugs hide.
8. Large models. Check a random subset of coordinates or a few random directions, layer by layer, on a tiny batch.

Interactive figure: loading…

Library defaults, as a reference. torch.autograd.gradcheck uses central differences with eps=1e-6 and passes when every entry satisfies

with atol=1e-5 and rtol=1e-3. It expects float64 inputs. SciPy's scipy.optimize.check_grad returns the absolute norm $\|g - g_{\text{num}}\|$ with forward differences; divide by the norms yourself.

The Taylor test (lesson 03) avoids a single threshold: the remainder $|f(x + td) - f(x) - t\, g^\top d|$ shrinks like $t^2$ with a correct $g$ and only like $t$ with a wrong one.

Answers to section 1. (1) $h = 10^{-8}$, close to $h^* = 2\sqrt{\varepsilon} \approx 3\cdot10^{-8}$ (here $F = M_2 = e$). The error is about $7\cdot10^{-9}$: about 8 correct digits. At $10^{-4}$ truncation dominates (error $1.4\cdot10^{-4}$); at $10^{-12}$ rounding dominates (error $4\cdot10^{-4}$). (2) No difference formula can: in float64, $1 + 10^{-20} = 1$, so the numerator is exactly 0. The complex step can: $\operatorname{Im} e^{1 + 10^{-20}i}/10^{-20}$ equals $e$ to all 16 digits.

4. Worked example #

(a) Exact truncation errors. Take $f(x) = x^3$ at $x = 1$: $f'(1) = 3$, $f''(1) = 6$, $f''' = 6$. For a cubic every truncation error can be written exactly.

Forward. Expand $(1+h)^3 = 1 + 3h + 3h^2 + h^3$. So

The leading error term $3h$ is $\tfrac{h}{2} f''(1)$. With $h = 0.1$: $D_+ = 3.31$, error $0.31$. With $h = 0.05$: $D_+ = 3.1525$, error $0.1525$, about half: first order.

Central. $(1+h)^3 - (1-h)^3 = 6h + 2h^3$. So

exactly $\tfrac{h^2}{6} f'''$. With $h = 0.1$: $3.01$, error $0.01$. With $h = 0.05$: $3.0025$, error $0.0025$, a quarter: second order.

Richardson. $\big(4 \cdot 3.0025 - 3.01\big)/3 = (12.01 - 3.01)/3 = 3$, exact, because the error was a pure $h^2$ term.

Complex step. Expand with $i^2 = -1$ and $i^3 = -i$:

So

the same size of error as the central difference, with the opposite sign (Theorem 3: $-\tfrac{h^2}{6}\cdot 6$). With $h = 0.1$: $2.99$. With $h = 10^{-20}$ the error is $10^{-40}$.

(b) Optimal steps. A cubic at $x = 1$ is a poor test of rounding (float64 computes $(1+h)^3 - 1$ almost exactly), so take $f = e^x$ at $x = 1$, where $F = M_2 = M_3 = e$.

• Forward: $h_+^* = 2\sqrt{\varepsilon e/e} = 2\sqrt{\varepsilon} \approx 3.0\cdot10^{-8}$, with error bound $2\sqrt{\varepsilon e\cdot e} = 2e\sqrt{\varepsilon} \approx 8.1\cdot10^{-8}$.
• Central: $h_0^* = (3\varepsilon)^{1/3} \approx 8.7\cdot10^{-6}$, with error bound $\tfrac12 3^{2/3}(\varepsilon e)^{2/3} e^{1/3} \approx 1.0\cdot10^{-10}$.

Section 5 measures $6.6\cdot10^{-9}$ at $h = 10^{-8}$ (forward) and $5.9\cdot10^{-11}$ at $h = 10^{-5}$ (central): below the bounds, as expected for worst-case bounds.

(c) A gradient check by hand. $f(x, y) = x^2 y$ at $(1, 2)$. The true gradient is $(2xy, x^2) = (4, 1)$. A buggy backward pass returns $g = (4, 1.01)$: a 1% error in one entry.

$f$ is quadratic in $x$ and linear in $y$, so central differences are exact here (the $f'''$ term is 0): $g_{\text{num}} = (4, 1)$. Then

far above $10^{-7}$: a bug. With the correct $g$ and $h = 10^{-6}$, float64 gives $\text{rel} \approx 10^{-11}$.

5. In code #

NumPy and SciPy check the numbers above. The notebook has these cells and the exercises.

import numpy as np
from scipy.optimize import check_grad

eps = np.finfo(float).eps                       # 2.2e-16 (float64 machine epsilon)

def forward(f, x, h):
    return (f(x + h) - f(x)) / h

def central(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

def complex_step(f, x, h):
    return f(x + 1j * h).imag / h

# 1. The V for f = exp at x = 1 (f'(1) = e). Columns: forward, central, complex step.
for h in (1e-2, 1e-5, 1e-8, 1e-12, 1e-20):
    errs = [abs(m(np.exp, 1.0, h) - np.e) for m in (forward, central, complex_step)]
    print(f"h={h:.0e}  " + "  ".join(f"{e:.1e}" for e in errs))
print(f"theory h*: {2 * np.sqrt(eps):.1e} (forward), {(3 * eps) ** (1 / 3):.1e} (central)")
# h=1e-02  1.4e-02  4.5e-05  4.5e-05
# h=1e-05  1.4e-05  5.9e-11  4.5e-11
# h=1e-08  6.6e-09  6.6e-09  0.0e+00
# h=1e-12  4.3e-04  2.1e-04  0.0e+00
# h=1e-20  2.7e+00  2.7e+00  0.0e+00     <- x + h == x: the difference is 0
# theory h*: 3.0e-08 (forward), 8.7e-06 (central)

# 2. A gradient check with a relative error
def numerical_gradient(f, x, h=1e-6):
    g = np.zeros_like(x)
    for i in range(x.size):
        e = np.zeros_like(x)
        e[i] = h * max(1.0, abs(x[i]))          # relative step
        g[i] = (f(x + e) - f(x - e)) / (2 * e[i])
    return g

def rel_error(a, b):
    den = max(np.linalg.norm(a), np.linalg.norm(b))
    return 0.0 if den == 0 else np.linalg.norm(a - b) / den

f = lambda z: z[0] ** 2 * z[1]
grad_ok = lambda z: np.array([2 * z[0] * z[1], z[0] ** 2])
grad_bug = lambda z: np.array([2 * z[0] * z[1], 1.01 * z[0] ** 2])   # 1% error in one entry
x0 = np.array([1.0, 2.0])
g_num = numerical_gradient(f, x0)
print(f"{rel_error(grad_ok(x0), g_num):.0e} {rel_error(grad_bug(x0), g_num):.1e}")   # 7e-12 2.4e-03

# 3. A kink: ReLU at 1e-7 from 0, with the step 1e-6
relu = lambda z: np.maximum(z, 0.0)
print(central(relu, 1e-7, 1e-6))               # 0.55, but the true derivative is 1

# 4. SciPy's check_grad returns an absolute error (forward differences)
print(f"{check_grad(f, grad_bug, x0):.4f}")    # 0.0100

Read the output. The forward error falls like $h$ down to about $10^{-8}$, and the central error like $h^2$ down to about $10^{-5}$; below these steps rounding takes over. The rounding arm is noisy (single steps can be lucky), so judge the trend, not one point. The complex step is exact to the last digit at every small $h$.

6. Where it shows up #

• Unit tests for gradients. torch.autograd.gradcheck (and gradgradcheck for second derivatives) and jax.test_util.check_grads, which compares forward and reverse mode with finite differences. Every custom autograd function, fused kernel or hand-written layer should pass one, including the layers of lesson 06.
• Debugging training. When a model trains badly after a change in the loss or a layer, a gradient check on a tiny batch is the fastest way to rule out a wrong gradient.
• Adjoints in science codes. Hand-written adjoints of PDE solvers and of variational data assimilation are tested with the gradient test, the ratio $\big(J(x + \alpha d) - J(x)\big)/\big(\alpha\, \nabla J(x)^\top d\big) \to 1$ as $\alpha \to 0$, and with the Taylor test. A linear operator and its adjoint (for example seismic modelling and migration) are tested with the dot-product test $\langle Au, v\rangle = \langle u, A^\top v\rangle$ for random $u, v$.
• Black-box simulators. When legacy code (for example Fortran without AD) must give sensitivities, finite differences or the complex step (recompile with complex types) are the practical options. Aerospace design uses the complex step for this.
• Stencils. The same Taylor analysis gives the order of accuracy of finite-difference stencils in PDE solvers and in grid-based physics-informed residuals.

7. Common mistakes #

• float32 for the check. A correct gradient fails, or a small bug hides in the noise.
• $h$ "as small as possible". $h = 10^{-12}$ is on the rounding arm of the V. Also wrong: one fixed absolute $h$ for inputs of very different sizes; use a relative step.
• Absolute error. $\|g - g_{\text{num}}\| = 10^{-5}$ means nothing without $\|g\|$.
• Crossing a kink. A failure at one sample of a ReLU network may be a kink, not a bug: move the point or use a smooth activation, then check again.
• A non-deterministic $f$. Dropout, random augmentation, or batch normalization in training mode with a changing batch make $f(x+h) - f(x-h)$ meaningless.
• Complex step with unsafe code. np.abs(z) returns $\sqrt{x^2 + h^2}$, a real number: part of the derivative silently becomes 0. Use np.where(z.real >= 0, z, -z) instead.
• One passing point as proof. A bug term can vanish at one point; check a few random points.

8. Check yourself #

1. Central difference of $f(x) = x^4$ at $x = 1$ with $h = 0.1$: the value, the error, and the leading-term prediction?
   Answer

   $D_0 = (1.1^4 - 0.9^4)/0.2 = (1.4641 - 0.6561)/0.2 = 4.04$, error $0.04$. Prediction $\tfrac{h^2}{6} f'''(1) = \tfrac{0.01}{6} \cdot 24 = 0.04$ (in fact $D_0 = 4 + 4h^2$ exactly).
2. In float32 ($\varepsilon \approx 1.2\cdot10^{-7}$) with $F = M_2 = M_3 = 1$: optimal steps and best error bounds?
   Answer

   Forward: $h^* = 2\sqrt\varepsilon \approx 6.9\cdot10^{-4}$, error $\approx 6.9\cdot10^{-4}$. Central: $h^* = (3\varepsilon)^{1/3} \approx 7.1\cdot10^{-3}$, error $\approx 2.5\cdot10^{-5}$. Only 3–4 digits, which is why checks run in float64.
3. A claimed gradient is $(2, -1, 0.5)$; the central difference gives $(2, -1, 0.505)$. Relative error and verdict?
   Answer

   $\|g - g_{\text{num}}\| = 0.005$ and $\max(\|g\|, \|g_{\text{num}}\|) \approx 2.2924$, so $\text{rel} \approx 2.2\cdot10^{-3}$. For a smooth function in float64 this is almost surely a bug, in the third entry.
4. $f(x) = \max(x, 0)$ at $x = 10^{-7}$, central difference with $h = 10^{-6}$: what do you get, and why?
   Answer

   $(1.1\cdot10^{-6} - 0)/(2\cdot10^{-6}) = 0.55$ instead of 1: $x - h < 0$ crosses the kink. Use a point farther than $h$ from 0, or a smaller step for this point.
5. Why can the complex step use $h = 10^{-20}$ while the central difference cannot?
   Answer

   At $x = 1$, $x \pm 10^{-20}$ rounds to $x$, so the central numerator is exactly 0. The complex step stores $h f'(x)$ in the imaginary part without any subtraction, so its relative error stays about $\varepsilon$; its truncation error $h^2 f'''/6$ is negligible.
6. Why divide by $\max(\|g\|, \|g_{\text{num}}\|)$ and not by $\|g\|$ alone?
   Answer

   It is symmetric in the two gradients, bounded (by 2), and defined when $g = 0$ but $g_{\text{num}} \ne 0$ (common with ReLUs). Only "both exactly zero" needs a rule: return 0.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the forward and central difference formulas and their orders.

   Show answer

2. Why is the central difference second order?

   Show answer

3. How big is the rounding error of a finite difference with step $h$?

   Show answer

4. What are the optimal steps and best errors of the forward and central differences in float64 (all scales 1)?

   Show answer

5. Write the complex-step formula for $f'(x)$.

   Show answer

6. Why can the complex step use $h = 10^{-20}$?

   Show answer

7. Which relative error do we use in a gradient check?

   Show answer

8. A gradient check fails at one point of a ReLU network. What should you check before you look for a bug?

   Show answer

9. Name three rules for a reliable gradient check.

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you know your backward pass is correct?" — A gradient check: float64, central differences, $h \approx 10^{-6}$ relative to the scale, relative error with thresholds, no kinks, randomness off, a few random coordinates or directions per layer.
• "Why not train with finite differences?" — $2n$ evaluations per gradient ($n$ parameters, about $10^9$ for an LLM) against a small constant times one evaluation for reverse mode, and at best about $10^{-10}$ accuracy.
• "Why is the central difference better, and which step do you use?" — The even Taylor terms cancel, so the error is $O(h^2)$. Balancing $h^2$ against $\varepsilon/h$ gives $h \propto \varepsilon^{1/3} \approx 10^{-5}$.
• "Why does a smaller step make the result worse?" — Cancellation: $f(x+h)$ and $f(x)$ share their leading digits, the difference is mostly rounding noise, and dividing by $h$ magnifies it.
• "What is the complex step?" — $\operatorname{Im} f(x+ih)/h$: error $O(h^2)$ and no cancellation, so $h = 10^{-20}$ gives machine precision; it needs analytic, complex-safe code.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. forward_diff, central_diff, complex_step; then reproduce the V plot.
2. optimal_step(method, eps, F, M): the step of Theorem 2, compared with the observed minimum of the V.
3. numerical_gradient and numerical_jacobian with relative steps.
4. gradient_check: relative error of a hand-written MLP backward pass; find the planted bug.
5. A kink experiment: how often a gradient check of a ReLU layer fails, as a function of $h$.

11. Watch and read #

Gradient Checking (C2W1L13)Andrew Ng · DeepLearning.AI, Improving Deep Neural Networks · 7 min

The practical recipe: central differences, a relative error (with a sum in the denominator) and thresholds like the ones in this lesson.

Gradient Checking Implementation Notes (C2W1L14)Andrew Ng · DeepLearning.AI, Improving Deep Neural Networks · 5 min

Regularization terms, dropout and other traps; compare with the rules in section 3.

ME564 Lecture 14: Numerical differentiation using finite differenceSteve Brunton · University of Washington, ME564 · 50 min

Forward, backward and central differences derived from Taylor series, with their orders.

Lecture 3 Part 2: Finite-Difference ApproximationsAlan Edelman and Steven G. Johnson · MIT 18.S096 Matrix Calculus for Machine Learning and Beyond · 51 min

Why and how to check derivatives with finite differences, and how accurate the approximation is.

• Driscoll and Braun, Fundamentals of Numerical Computation: §5.4 Finite differences and §5.5 Convergence of finite differences (rounding and the optimal step).
• Nick Higham, What is the complex step approximation?
• Martins, Sturdza and Alonso, "The complex-step derivative approximation", ACM Transactions on Mathematical Software 29(3), 2003: the standard reference, with the test function $e^x/\sqrt{\sin^3 x + \cos^3 x}$ used in the figure.
• CS231n, Neural Networks Part 3, section "Gradient checks": practical rules for networks.