The inverse problem and parametrizations
LH toy, MSHT20 Chebyshev form in full, NNPDF nets, grids — with the MSHT20 basis written out.
A PDF fit asks a question that is, taken literally, impossible: reconstruct seven unknown functions of $x$ from 648 numbers. This chapter is about the move that makes it possible — choosing a parametrization — and about every parametrization this lab has actually run: the Les Houches toy that validated the machinery (PPDF-1…7), the MSHT20 Chebyshev form that finally described real data (PPDF-12), neural networks, and grids. At the end, the Bayesian evidence referees them all.
The inverse problem, formalized
Chapter 3.3 will show that every DIS prediction is a linear functional of the PDFs at the parametrization scale: $T_I[f] = \sum_{j,k} \mathrm{FK}_{Ijk}\, f_j(x_k)$. A fit is then the inverse problem
with $D$ a vector of 648 measurements and $f$ a point in an infinite-dimensional function space. This is ill-posed in exactly Hadamard's sense. Uniqueness fails: the map $f \mapsto T[f]$ has an enormous null space — any function wiggle orthogonal to all 648 functionals (wiggles between grid nodes, flavour combinations photons don't couple to, everything at $x$ outside the data's reach) changes the prediction by nothing, so the data cannot distinguish $f$ from $f + \delta f$. Stability fails: inverting a smoothing operation (the FK convolution averages $f$ over $x$) amplifies noise; nudge one datum within its error bar and the "exact" inverse lurches wildly. Some regularization — extra information beyond the data — is mandatory, not optional.
Every mainstream fit regularizes the same way: restrict $f$ to a finite-dimensional family $f(x; \theta)$, $\theta \in \mathbb{R}^d$, and infer $\theta$. The null space is tamed because the family cannot wiggle arbitrarily; stability comes from the family's smoothness. The price is parametrization bias: if the true proton is not in the family, no amount of data brings you to it. Fits differ only in $d$ and in the geometry of the family — everything from $d = 13$ (Les Houches) through 52 (MSHT20) to $\sim$800 (NNPDF's networks).
The universal backbone — and why physics picks it
Almost every fixed-form parametrization since the 1980s shares one skeleton:
power-law endpoints, with a flexible interpolating factor $\mathcal{P}$ in between. The two powers are not guesses — each comes from a real argument.
At small $x$ the virtual photon–proton system has invariant mass $W^2 \approx Q^2/x \to \infty$: DIS becomes a high-energy $\gamma^* p$ total cross-section. Regge theory says high-energy total cross-sections are governed by $t$-channel exchanges of whole trajectories, $\sigma_{\rm tot} \sim (W^2)^{\alpha(0) - 1}$, where $\alpha(0)$ is the trajectory's intercept. Since $F_2 \propto Q^2 \sigma_{\gamma^* p}$ at fixed $Q^2$,
$$ F_2(x) \;\sim\; (W^2)^{\alpha(0)-1} \;\sim\; x^{\,1 - \alpha(0)} \qquad (x \to 0) $$
so $xf \sim x^{\delta}$ with $\delta = 1 - \alpha(0)$. Flavour-singlet combinations (sea, gluon) couple to the Pomeron, $\alpha_P(0) \approx 1.08$, predicting $\delta \approx -0.08$: a slow small-$x$ rise. Valence (non-singlet) combinations couple to the $\rho/\omega/f/a$ Reggeons, $\alpha_R(0) \approx 0.5$, predicting $\delta \approx 0.5$: valence vanishes at small $x$. Compare MSHT20's fitted values: $\delta_{u_V} = 0.34$, $\delta_S = -0.010$ — the right ordering, with the exact powers left to the data (DGLAP evolution modifies the naive intercepts).
At $x \to 1$ one parton must carry all the proton's momentum, and the other constituents — the $n_s$ "spectators" — must be squeezed to zero. Squeezing a spectator costs a hard gluon exchange; each exchange brings a far-off-shell propagator $\sim 1/(1-x)$ into the amplitude. Counting powers (Drell–Yan–West, Brodsky–Farrar) gives
$$ f(x) \;\sim\; (1-x)^{2 n_s - 1} \qquad (x \to 1) $$
A valence quark in $|uud\rangle$ has $n_s = 2$ spectators: $\eta = 3$. A gluon (one extra unit from the helicity mismatch) gets $\eta \approx 5$; a sea quark, which needs its $q\bar q$ partner squeezed too ($n_s = 4$), gets $\eta \approx 7$. Now look at the MSHT20 publication fit our port carries: $\eta_{u_V} = 3.75$, $\eta_{g} = 5.25$, $\eta_S = 7.04$. Half-century-old counting rules, sitting inside a modern 52-parameter fit, to within a unit.
The Les Houches toy, written in full
Colibri's validation model — the one behind all our closure tests and the first real-data fit — dresses the backbone minimally. Its 13 free parameters are $(\alpha_g, \beta_g)$, $(\alpha_u, \beta_u, \epsilon_u, \gamma_u)$, $(\alpha_d, \beta_d, \epsilon_d, \gamma_d)$ and $(N_\Sigma, \alpha_\Sigma, \beta_\Sigma)$:
with $\mathcal{P} = 1 + \epsilon\sqrt{x} + \gamma x$ the classic MRST-era interpolator. Three normalizations are not free: $A_{u_V}, A_{d_V}$ come from the valence counting rules and $A_g$ from momentum (16 slots $\to$ 13 — the same sum-rule elimination §3.2 does for MSHT20, at training-wheels scale). The remaining evolution basis is assumed: $V = u_V + d_V$, $V_3 = u_V - d_V$, $V_8 = V_{15} = V_{24} = V_{35} = V$, $T_8 = 0.4\,\Sigma$, $T_{15} = T_{24} = T_{35} = \Sigma$, and — note it for §3.3 — $T_3 = 0$. This model earned its keep: L0/L1 closure (PPDF-1–3), the three-method cross-checks (PPDF-4–6). But on real data it is a straitjacket: $\chi^2/N = 5.60$ (PPDF-7).
MSHT20: the Chebyshev parametrization, completely
The form Maria pointed us to — and the one our JAX port (PPDF-10) runs — takes $\mathcal{P}$ to be a Chebyshev series in a transformed variable:
with the Chebyshev polynomials of the first kind, written out to the order you'll meet in every formula of §3.2:
Why Chebyshevs? On $y \in [-1,1]$ they oscillate between $\pm 1$ — a bounded, orthogonal basis, so coefficients stay $\mathcal{O}(1)$ and truncation is graceful (MSHT checked that going beyond $i = 6$ moves nothing). Why the $\sqrt{x}$ argument? Because data span decades of $x$, and a polynomial in $x$ itself wastes all its resolution at large $x$: under $y = 1 - 2x$, the entire small-$x$ world $x \le 10^{-2}$ crams into $y \in [0.98, 1]$ — 1% of the basis's domain. The square root stretches it (worked example (a) makes this quantitative). The fitted basis is $u_V,\ d_V,\ S,\ s_+ = s + \bar s,\ s_-= s - \bar s,\ \rho = \bar d/\bar u$, and the gluon — with two special cases:
The gluon's second term (fitted $A_{g'} = -0.37$, $\delta_{g'} = -0.28$) lets $xg$ turn over and even go negative at very small $x$ at the low input scale — one power law cannot do that. The $s_-$ factor $(1 - x/x_0)$ builds in the sign change that zero net strangeness forces (§3.2 computes $x_0$). The $\rho = \bar d/\bar u$ combination is fitted as a ratio — no $x^\delta$ factor at all. In the FPPDF bookkeeping our port mirrors, this is a 73-slot vector; sum rules and ties cut it to 52 free parameters (worked example (b) does the audit line by line).
Neural networks, and grids
NNPDF replaces $\mathcal{P}$ by a neural network: $x f_k(x) = A_k\, x^{1-\alpha_k} (1-x)^{\beta_k}\, \mathrm{NN}_k(x, \ln x)$, one small feed-forward net per evolution-basis flavour, $\sim$800 weights in total. The preprocessing powers $\alpha_k, \beta_k$ still carry the Regge/counting endpoints (the net cannot learn behaviour where there is no data); they are sampled within self-consistently determined ranges rather than fixed. The philosophy: make parametrization bias negligible and let the closure-test machinery certify the uncertainties. The cost: 800 weights need regularization by early stopping and ensembling, and "the fit" becomes a distribution over networks.
Grid models go maximally local: fit the PDF values at a handful of $x$-nodes and interpolate. We ran this too, and it taught us the pitfalls the hard way. Locality means node values far from data are set purely by the prior — and our uniform prior box (NNPDF40 central $\pm 5\sigma$) is secretly informative: if the SLAC+BCDMS-only truth lies outside the global-fit band anywhere, the model can't reach it. Node resolution binds (6 nodes per flavour in PPDF-9); sum rules aren't automatic; and the flavours you don't parametrize are silently zeroed in the FK contraction — the PPDF-8 lesson, retold technically in §3.3.
Every choice on this page injects information the data doesn't contain — powers, truncation orders, network architectures, prior boxes. That is not a flaw; it is the regularization that makes an ill-posed problem solvable. The scientific obligations are (1) to declare the choice, (2) to test whether conclusions survive changing it, and (3) where possible, to let a principled statistic — not aesthetics — arbitrate between choices. In this lab that statistic is the Bayesian evidence.
Bias, variance, and the evidence as referee
Rigid families have low variance and high bias: the 13-parameter toy produced beautifully stable posteriors — around a fit with $\chi^2/N = 5.60$ that no honest statistician would call a description of the data. Flexible families flip the tradeoff: more of the data's noise becomes fit-able, posteriors widen, and with the wrong inductive bias you can add parameters and get worse. The evidence $Z = \int \mathcal{L}(\theta)\,\pi(\theta)\, d\theta$ scores this tradeoff automatically — better likelihood helps, wasted prior volume hurts. Our own chain is the cleanest illustration we know: 36-parameter grid vs 13-parameter toy, $\Delta \log Z = -108.8$ (PPDF-9: flexibility lost — it didn't buy likelihood where it mattered); 52-parameter MSHT20 vs the same toy, $\Delta \log Z = +1453$ (PPDF-12: the right flexibility, overwhelming win, $\chi^2/N = 1.002$). Parameter count is not the axis that matters; where the family can flex is.
Map data-region $x$ values through both candidate arguments, $y = 1 - 2\sqrt{x}$ vs the naive $y = 1 - 2x$:
$$ \begin{array}{c|cc|cc} x & 10^{-4} & 10^{-2} & 0.1 & 0.5 \\ \hline y = 1-2\sqrt{x} & 0.980 & 0.800 & 0.368 & -0.414 \\ y = 1-2x & 0.9998 & 0.980 & 0.800 & 0.000 \end{array} $$
Under the naive map, two full decades $10^{-4} \le x \le 10^{-2}$ occupy $\Delta y = 0.02$ — 1% of the Chebyshev domain, where even $T_6$ is locally indistinguishable from a straight line: the basis is blind there. Under $\sqrt{x}$, the same decades get $\Delta y = 0.18$, and the well-measured region $x \in [10^{-2}, 0.75]$ spreads across $y \in [-0.73, 0.80]$ — over three-quarters of the domain. Equivalently, walk down $y$ in equal steps of 0.4 from 0.8 and invert $x = \left(\tfrac{1-y}{2}\right)^2$: you land at $x = 0.01,\ 0.09,\ 0.25,\ 0.49,\ 0.81$ — a graceful sweep from sea to valence territory, one Chebyshev oscillation per physics regime. That is the design: equal polynomial resolution per decade-ish of $x$, achieved by a change of variable instead of more parameters.
Our port stores the parametrization as a 73-slot vector (block = flavour, slots = $A, \delta, \eta$, then Chebyshev coefficients): $u_V$ 9 $+$ $d_V$ 9 $+$ $S$ 9 $+$ $s_+$ 9 $+$ $g$ 10 (two terms) $+$ $s_-$ 10 (includes $x_0$) $+$ $\rho$ 8 (no $\delta$ slot) $+$ $c_+$ 9 $\;=\; \mathbf{73}$. Now subtract, with reasons:
−9 charm block: perturbative charm, all zeros (charm is generated by evolution, not fitted).
−4 sum-rule outputs, computed not sampled: $A_{u_V}$, $A_{d_V}$ (valence counting), $A_g$ (momentum),
$x_0$ (zero strangeness) — the four analytic eliminations of §3.2.
−1 tie: $\delta_{s_+} \equiv \delta_S$ (the sea and total-strangeness small-$x$ powers are locked
together — F2 data cannot separate them).
−7 frozen by MSHT: $\delta_{s_-}$ (fixed at 0.0015) and all six $s_-$ Chebyshev coefficients (fixed at
zero — $s - \bar s$ is barely constrained even in MSHT's global fit; with our DIS-only data it would be pure
prior).
$73 - 9 - 4 - 1 - 7 = \mathbf{52}$ — exactly the dimension of the prior box in our PPDF-12 runcard: 8 ($u_V$) + 8 ($d_V$) + 9 ($S$) + 8 ($s_+$) + 9 ($g$) + 2 ($s_-$) + 8 ($\rho$). The same audit explains the sampler's 52-dimensional space, the "52" in every MSHT20 paper, and why closure found only $\sim$11 of these directions actually constrained by our dataset (PPDF-11).