The Hessian method in full
Quadratic expansion, eigenvector sets, tolerances, observable propagation — and the exact failure condition.
Every PDF error band you have ever seen from CTEQ or MSHT was made by the method of this section: expand the $\chi^2$ to second order around its minimum, diagonalize the curvature, and step along each eigendirection. It is fast, elegant, and exactly correct — under assumptions we can state precisely and, with our own fits, watch fail precisely. This chapter derives the entire machine: the $\Delta\chi^2 = 1$ rule from first principles, the covariance–Hessian relation, eigenvector PDF sets, the master formulas for propagating errors to any observable, the tolerance problem that haunts global fits, the computational routes to the Hessian (including the hard lesson our own PPDF-14 run taught us at 52 dimensions), and the exact condition under which the method stops meaning anything.
The quadratic world
Let $\chi^2(\theta)$ be the figure of merit of §3.4, minimized at $\theta_0$ with value $\chi^2_{\min}$. Taylor expansion around the minimum — where the gradient vanishes by definition — gives
Everything in the Hessian method follows from taking this quadratic form seriously — not as an approximation near the minimum, but as the whole likelihood. That substitution is the method's single load-bearing assumption, and every failure we will meet is a failure of it.
For Gaussian data the likelihood is $\mathcal L(\theta) \propto e^{-\chi^2(\theta)/2}$. With one parameter, write the quadratic form as $\chi^2 = \chi^2_{\min} + (\theta - \hat\theta)^2/s^2$. Then
$$ \mathcal L(\theta) \;\propto\; \exp\!\left[-\frac{(\theta-\hat\theta)^2}{2 s^2}\right] $$
— a Gaussian in $\theta$ with standard deviation exactly $s$. The point where $\chi^2$ has risen by $\Delta\chi^2 = 1$ is $|\theta - \hat\theta| = s$: the $1\sigma$ point, enclosing 68.27% of the probability. No convention entered; the "$\Delta\chi^2 = 1$ rule" is a theorem about Gaussians, not a choice. Hold that thought for the tolerance discussion below, where global fits replace 1 by $T^2 \sim 10$–$25$.
In many dimensions the same likelihood reads $\mathcal L \propto \exp\!\big[-\tfrac14\,\delta\theta^{\mathsf T} H\, \delta\theta\big]$. Matching to the general Gaussian $\exp\!\big[-\tfrac12\,\delta\theta^{\mathsf T} C_\theta^{-1}\,\delta\theta\big]$ gives
$$ C_\theta^{-1} = \tfrac{1}{2} H \qquad\Longleftrightarrow\qquad C_\theta = 2\,H^{-1} $$
The factor 2 is the perennial trap: $H$ is the Hessian of $\chi^2$, and $\chi^2 = -2\ln\mathcal L$. Had we defined $H$ from $\ln\mathcal L$ directly, the covariance would be the plain inverse. Note what this equation demands: $H$ must be invertible and positive definite. That is the method's existence condition, and it is exactly what our real-data run violates.
Diagonalize the (symmetric) Hessian: $H v_k = \lambda_k v_k$ with orthonormal $v_k$. Define rescaled coordinates $z_k \equiv \sqrt{\lambda_k/2}\; (v_k^{\mathsf T}\delta\theta)$. Then
$$ \Delta\chi^2 = \tfrac12 \sum_k \lambda_k \,(v_k^{\mathsf T}\delta\theta)^2 = \sum_k z_k^2 $$
— in $z$-space, constant-$\Delta\chi^2$ surfaces are perfect spheres and every direction carries error 1. The eigenvector PDF sets are the $2p$ points reached by stepping a distance $T$ along each axis of this sphere:
$$ S_k^\pm:\qquad \theta = \theta_0 \pm T\sqrt{\frac{2}{\lambda_k}}\; v_k \qquad (\Delta\chi^2 = T^2) $$
Each set is a full PDF — evolved, convolved, delivered as an LHAPDF member. The $1\sigma$ excursion along direction $k$ is $\sqrt{2/\lambda_k}$: soft directions (small $\lambda_k$) step far, stiff directions barely move. Keep the scaling $\sigma_k \propto 1/\sqrt{\lambda_k}$ in view; it is the fuse on the bomb.
Linearize an observable around the minimum: $\mathcal O(\theta) \approx \mathcal O_0 + \nabla\mathcal O^{\mathsf T}\delta\theta$. Its variance under the Gaussian of Derivation 2 is $\sigma_{\mathcal O}^2 = \nabla\mathcal O^{\mathsf T} C_\theta \nabla\mathcal O = \sum_k (2/\lambda_k)(\nabla\mathcal O^{\mathsf T} v_k)^2$. Now evaluate $\mathcal O$ on the sets: $\mathcal O(S_k^\pm) - \mathcal O_0 \approx \pm T\sqrt{2/\lambda_k}\, \nabla\mathcal O^{\mathsf T} v_k$, so the finite differences across sets reconstruct exactly the gradient projections, and
$$ \Delta\mathcal O_{\rm sym} \;=\; \frac{1}{2T}\sqrt{\sum_k \big[\mathcal O(S_k^+) - \mathcal O(S_k^-)\big]^2} \;=\; \sqrt{\nabla\mathcal O^{\mathsf T}\, C_\theta\, \nabla\mathcal O} $$
— the symmetric master formula, an identity in the linear-quadratic world. When the response is visibly asymmetric one uses the one-sided version,
$$ \Delta\mathcal O_\pm \;=\; \frac{1}{T}\sqrt{\sum_k \Big[\max\big(\pm\{\mathcal O(S_k^+)-\mathcal O_0\},\, \pm\{\mathcal O(S_k^-)-\mathcal O_0\},\, 0\big)\Big]^2} $$
which keeps, for each direction, whichever member moves the observable up (for $\Delta_+$) or down (for $\Delta_-$). This is why a Hessian PDF release is just $2p$ member grids plus these two formulas: the entire covariance has been pre-packaged into finitely many evaluations.
The tolerance problem: $T \neq 1$
Derivation 1 says the textbook value is $T = 1$. Yet CTEQ and MSHT ship error sets built at $T \sim 3$–$5$ (so $\Delta\chi^2 \sim 10$–$25$). Why? Because the global-fit likelihood is not the Gaussian its $\chi^2$ pretends: (i) datasets are mutually inconsistent at the level of their quoted errors — different experiments pull shared parameters to incompatible values, so a $\Delta\chi^2=1$ band around the compromise minimum covers neither camp; (ii) the parametrization is imperfect — a rigid functional form cannot bend to every dataset, and its stiffness masquerades as spurious precision; (iii) unquantified theory error (higher orders, higher twist) lives in nobody's covariance matrix. MSHT's dynamic tolerance makes the inflation direction-by-direction: walk along eigendirection $k$ until the first individual dataset $n$ exits its own 68% confidence range (a percentile criterion on each $\chi^2_n$, rescaled to its number of points), and set $T_k$ there — typically $T_k \approx 2$–$6$, averaging near 3, with each member's binding dataset recorded. Statistically, $T \neq 1$ is a confession: it concedes that the stated likelihood is wrong, and repairs it by scaling the parameter covariance by $T^2$ — an ad-hoc but disclosed patch, tuned per direction, in place of a model of the tensions themselves. The Bayesian route of topic 4 must instead confront those tensions explicitly (in the likelihood or priors); it has no dial equivalent to $T$.
Computing $H$: three routes and one lab lesson
With $\chi^2 = r^{\mathsf T} C^{-1} r$ and residuals $r_i = T_i(\theta) - D_i$, differentiate twice ($J_{ia} = \partial T_i/\partial\theta_a$):
The first piece is the Gauss–Newton approximation $H_{\rm GN} = 2 J^{\mathsf T} C^{-1} J$. The second
vanishes identically when the residuals vanish at the minimum ($r = 0$: perfect fit) or the model is linear
($\partial^2 T = 0$) — in either case GN is the exact Hessian. Note the structural difference:
$H_{\rm GN}$ is positive semi-definite by construction (it is a Gram matrix), while the exact Hessian can go
indefinite when nonzero residuals couple to model curvature. The remaining routes are automatic differentiation
(exact, when tractable) and finite differences of gradients or of $\chi^2$ values (approximate, but assumption-free).
In PPDF-14 we needed the exact $H$ of the 52-parameter MSHT20 model and hit an
engineering wall first: three XLA autodiff routes — fused jax.hessian, Hessian-vector-product columns,
and finite differences of the jitted gradient — each spent over 40 minutes compiling at 52 dimensions
before any arithmetic. The fix was to abandon compilation entirely: central finite differences of eager-mode
likelihood values, step $5\times10^{-4}\max(1,|\theta_i|)$, symmetrized $(H + H^{\mathsf T})/2$. Method
lesson: at high parameter count, the numerically humblest route can be the only one that terminates — and it is
the one that told us the truth below.
Since $C_\theta = 2H^{-1}$, the quantity $\lambda_k/2$ is the inverse variance — the precision — the data supply about direction $v_k$. The Hessian spectrum is therefore an information audit: large $\lambda$ = well-measured combination; $\lambda \to 0$ = the data are silent about $v_k$ and the band $\sqrt{2/\lambda_k} \to \infty$ honestly reports it; $\lambda < 0$ (possible once noisy residuals meet model curvature) = the "minimum" is not even a minimum at numerical resolution, and the formalism returns negative variance — no answer, not a wide one. The field's standard mitigation is parameter freezing: fix the offending combinations at chosen values and invert the reduced Hessian. That works — MSHT20 fits 52 parameters but builds its eigenvector sets from a reduced subset of directions — but recognize it for what it is: a delta-function prior, imposed by hand, carrying zero width, and nowhere advertised in the final band. The Bayesian analysis of §5.3 will make exactly this prior explicit and give it honest width.
Worked example — a two-parameter fit with one flat direction
Given: a two-parameter fit whose curvature at the minimum is
$$ H = \begin{pmatrix} 4.00 & 3.96 \\ 3.96 & 4.00 \end{pmatrix} $$
(this arises from data that mostly constrain the sum $\theta_1 + \theta_2$ — e.g. two nearly identical predictions).
Step 1 — spectrum. Eigenvalues $\lambda_\pm = 4.00 \pm 3.96$: $\lambda_1 = 7.96$ along $v_1 = (1,1)/\sqrt2$, $\lambda_2 = 0.04$ along $v_2 = (1,-1)/\sqrt2$. Condition number 199.
Step 2 — $1\sigma$ steps. $\sigma_k = \sqrt{2/\lambda_k}$: stiff direction $\sqrt{2/7.96} = 0.501$; soft direction $\sqrt{2/0.04} = 7.07$. The $\Delta\chi^2=1$ ellipse is 14× longer than wide; the members $S_2^\pm$ sit $\pm7.07$ away along $(1,-1)/\sqrt2$ while $\chi^2$ barely notices.
Step 3 — marginal error on $\theta_1$. $C_\theta = 2H^{-1}$: with $\det H = 16 - 15.6816 = 0.3184$, $C_{11} = 2 \times 4.00/0.3184 = 25.1$, so $\sigma(\theta_1) = 5.01$ — dominated entirely by the soft direction.
Step 4 — freeze $\theta_2$. Fixing $\theta_2 \equiv \theta_2^0$ leaves the 1×1 Hessian $H_{11} = 4.00$, hence $\sigma(\theta_1) = \sqrt{2/4} = 0.71$. The quoted error just shrank by a factor 7.1 — not because data arrived, but because a delta prior silently absorbed the flat direction. And if numerical noise had pushed $\lambda_2$ from $+0.04$ to $-0.04$: negative variance, no band at all. This is the entire pathology of the next step, in miniature.
Given (from PPDF-14, real SLAC+BCDMS data, MSHT20's 52 parameters on 648 points): exact-Hessian estimate at the verified minimum ($\chi^2/N = 0.9845$, found identically 4×) has 11 of 52 eigendirections negative (each consistent with zero at finite-difference resolution), and among the positive ones, parameter variances up to $C_{kk} \sim 1.3\times10^{9}$; the framework retains 20/52 directions after thresholding.
Step 1 — the excursion scale. A variance of $1.3\times10^9$ means a $1\sigma$ member step of $\sqrt{1.3\times10^9} \approx 3.6\times10^4$ in parameters whose priors — and physically sensible ranges — have width $\sim 0.4$. Member excursions recorded at $\pm 7\times10^4$ are exactly this scale (the precise factor depends on the per-direction step normalization): the sets $S_k^\pm$ ask for MSHT20 exponents five orders of magnitude outside any meaningful region, the same runaway that PPDF-4 saw parameters ride to $10^9$ along closure-test flat directions.
Step 2 — why the bands explode by $50$–$10^7$×. Through the master formula, an observable's band is $\sum_k (2/\lambda_k)(\nabla\mathcal O^{\mathsf T}v_k)^2$: directions with $\lambda_k \to 0^+$ contribute divergently wherever $\nabla\mathcal O$ has any overlap with them. Flavour combinations nearly orthogonal to the null space (the mid-$x$ gluon) inflate "only" by $\mathcal O(50)$; combinations living along it inflate by up to $\sim 10^7$ (PPDF-15 measures 50–$3\times10^7$ against the Bayesian reference).
Step 3 — and the 11 negative directions? They cannot enter at all — $\sqrt{2/\lambda}$ is imaginary — so they are dropped: a freezing decision made not by physics judgment but by the sign of numerical noise. At L1 closure (PPDF-4) the same disease was partial and quiet: the Hessian gluon band merely under-covered, with pulls against truth reaching $\approx 6\sigma$ while Bayesian and MC pulls stayed $\le 1.1$. On real data it is total: 52 parameters against 648 $F_2$ points leave the quadratic world, and with it the method, undefined. MSHT avoid this in practice with $\sim$4600 global points plus explicit freezing — which is precisely the data-sufficiency requirement this run exposes (background on the parametrization: see the wiki, MSHT20 entry).