Jacobian Computation — Finite Difference and AAD

Setup

The Levenberg-Marquardt algorithm requires the Jacobian $J(\theta) \in \mathbb{R}^{N \times p}$ at every iteration. For a Heston calibration with $N = 40$ instruments and $p = 5$ parameters, $J$ has 200 entries. For an LMM calibration with $N = 100$ swaptions and $p = 50$ parameters, $J$ has 5000 entries — each requiring a pricing function evaluation or a derivative computation.

How $J$ is computed determines the speed and accuracy of the entire calibration pipeline. There are four approaches: (i) analytic differentiation, (ii) finite differences, (iii) complex step differentiation, and (iv) automatic differentiation (AAD). Each involves a different trade-off between implementation cost, accuracy, and computational load.

Assumptions:

• The pricing function $f: \mathbb{R}^p \to \mathbb{R}$ (or $f: \mathbb{R}^p \to \mathbb{R}^N$) is smooth in $\theta$ on the domain of interest.
• Floating-point arithmetic: all computations are in IEEE 754 double precision, machine epsilon $\varepsilon_{\text{mach}} \approx 2.2 \times 10^{-16}$.
• For AAD: the pricing function is implemented as a sequence of elementary operations (addition, multiplication, $\exp$, $\ln$, etc.) — the standard assumption for any AD tool.

Theory

Finite Differences: Forward and Central

Error analysis: By Taylor's theorem:

$f(x + h) = f(x) + h f'(x) + \frac{h^2}{2} f''(x) + \frac{h^3}{6} f'''(x) + \cdots$

For the forward difference: $\frac{f(x+h) - f(x)}{h} = f'(x) + \underbrace{\frac{h}{2} f''(x)}_{\text{truncation}} + O(h^2)$

Truncation error: $O(h)$. Roundoff error from floating-point evaluation of $f(x+h) - f(x)$: $O(\varepsilon_{\text{mach}} |f(x)| / h)$.

For the central difference: $\frac{f(x+h) - f(x-h)}{2h} = f'(x) + \underbrace{\frac{h^2}{6} f'''(x)}_{\text{truncation}} + O(h^4)$

Truncation error: $O(h^2)$. Roundoff error: $O(\varepsilon_{\text{mach}} |f(x)| / h)$ (same structure, since $f(x)$ cancels but floating-point errors in $f(x\pm h)$ persist).

Central FD is $\sim 10^4$ times more accurate than forward FD at their respective optimal step sizes, at the cost of twice as many function evaluations. In calibration, the reduced error is almost always worth the extra evaluation.

Cost for the Jacobian:

For the Jacobian, all $N$ rows share the same perturbations: compute $f(\theta + he_j)$ for $j=1,\ldots,p$, compute all $N$ residuals at each perturbed point. Total: $p+1$ (forward) or $2p$ (central) calls to the full pricing engine (all $N$ instruments).

The Complex Step Method

Derivation: By the Taylor expansion of $f$ in $\mathbb{C}$: $f(x + ih) = f(x) + ih\,f'(x) - \frac{h^2}{2} f''(x) - \frac{ih^3}{6} f'''(x) + \cdots$

Taking the imaginary part: $\mathrm{Im}[f(x + ih)] = h f'(x) - \frac{h^3}{6} f'''(x) + \cdots$

Dividing by $h$: $f'(x) = \mathrm{Im}[f(x+ih)]/h + O(h^2)$.

The key advantage: there is no subtraction of nearly-equal quantities. Catastrophic cancellation — the main source of roundoff error in FD — does not occur. One can use $h = 10^{-200}$ and achieve machine-precision accuracy.

Automatic Differentiation (AAD)

Finite differences and complex step compute one column of $J$ per perturbation (forward mode in the sense of perturbation direction). When $p \ll N$ (few parameters, many instruments), this is efficient. When $N \ll p$ or when computing a scalar objective gradient, the reverse mode of AAD computes all partial derivatives in a single backward pass.

For the calibration Jacobian $J \in \mathbb{R}^{N \times p}$, reverse mode computes one row $\nabla_\theta f_i$ per backward pass — $N$ backward passes in total. Forward mode computes one column per pass — $p$ forward passes in total. The efficient choice depends on the ratio $N/p$:

For Heston ($p=5$, $N=40$): forward mode (5 passes) vs. reverse mode (40 passes) — forward wins. For LMM ($p=50$, $N=100$): reverse mode (100 passes) vs. forward mode (50 passes) — forward mode still wins on pass count, but each pass is cheaper with AAD than with FD.

Application: Analytic Vega in Black-Scholes

For a European call under Black-Scholes:

$C(S, K, r, T, \sigma) = S\,N(d_+) - K e^{-rT} N(d_-)$

where $d_\pm = [\ln(S/K) + (r \pm \frac{\sigma^2}{2})T] / (\sigma\sqrt{T})$.

The analytic derivative with respect to $\sigma$ (vega) is:

$\frac{\partial C}{\partial \sigma} = S\sqrt{T}\,n(d_+)$

where $n$ is the standard normal PDF. This is the benchmark for validating FD and AAD implementations: the complex step or reverse-mode AAD result must match this to within $O(\varepsilon_{\text{mach}})$.

Validation

The companion notebook verifies the four derivative methods on the Black-Scholes vega and on a 2-parameter nonlinear function, confirming:

1. Forward FD error decreases as $O(h)$ for $h > h^*$ and increases as $O(1/h)$ for $h < h^*$.
2. Central FD error decreases as $O(h^2)$ until the roundoff floor.
3. Complex step achieves machine-precision accuracy for all tested $h$ values.
4. A simple reverse-mode AAD implementation produces the correct gradient on a small model.
5. The optimal step size predictions of Theorem 3.1 are verified empirically.

Limitations

Other limitations:

• Parallelism: Column-wise FD (one column of $J$ per perturbation) parallelises trivially over $p$ threads. Reverse-mode AAD is sequential per backward pass; parallelism requires batching.
• Memory cost of AAD: Reverse mode requires storing the tape (the entire computation graph) between the forward and backward passes. For complex pricing functions, this tape can be large. Checkpointing techniques trade recomputation for memory.
• Implicit functions: If the pricing function involves solving an implicit equation (e.g., Newton's method for implied vol inversion), the Jacobian requires the implicit function theorem: $\partial \sigma_{\text{impl}} / \partial \theta = -(\partial G/\partial \sigma)^{-1} \partial G/\partial \theta$ where $G(\sigma, \theta) = 0$ defines the implied vol. This must be handled explicitly in AAD — naively differentiating through the Newton iterations works but is less efficient.

Interview Angle