Gaussian Process Eigenmodes for Statistical and Systematic Uncertainties in Template Fits

This paper proposes replacing traditional per-bin Barlow-Beeston factors and interpolation modifiers with a unified eigenmode basis derived from log-Gaussian Cox process posteriors to efficiently model both statistical and systematic uncertainties in LHC template fits, thereby reducing dimensionality while preserving or bounding the original variance.

The Gist

Imagine you are trying to find a tiny, rare gem (a new particle) hidden inside a massive, noisy pile of sand (background data) at a giant particle collider. To do this, physicists use a "template"—a map of what the sand pile should look like if no gem is there. They compare their actual observations to this map. If the real pile has a weird bump that the map doesn't predict, that might be the gem.

The problem is that making this map is tricky. The map is built from computer simulations (Monte Carlo), which are like taking a limited number of photos of the sand. If you don't have enough photos, the map gets grainy and full of "static" (statistical noise). If you try to make the map too detailed to see the gem clearly, the static gets so loud you can't trust the map at all.

This paper proposes a new way to build that map using Gaussian Processes (GPs), which is a fancy mathematical way of saying "smooth, intelligent guessing."

Here is the breakdown of the paper's ideas using simple analogies:

1. The Old Way: The "Pixelated" Map

Traditionally, physicists build their map by dividing the data into tiny boxes (bins) and counting the sand in each box.

• The Problem: If you have a limited number of simulation photos, some boxes will be empty or have very few grains. To handle the uncertainty of these empty boxes, the old method adds a "wobble factor" (a nuisance parameter) to every single box.
• The Consequence: If you have a 3D map with millions of boxes, you end up with millions of wobble factors. It's like trying to steer a ship by adjusting a separate rudder for every single plank of wood. It's computationally heavy, and when the data is scarce, the map becomes so shaky that it might hide the gem or create fake ones.

2. The New Way: The "Smooth River" Map

The authors suggest replacing the pixelated boxes with a smooth, flowing river (a mathematical function). Instead of counting grains in boxes, they use a Gaussian Process to draw a smooth curve that fits the sand data.

• The Magic: Because the curve is smooth, it "knows" that if one part of the river is high, the neighbors are likely high too. It borrows strength from its neighbors.
• The Result: Even with very few photos (low statistics), the map stays smooth and reliable. It doesn't get grainy. The paper proves mathematically that this smooth map is always more precise (has less uncertainty) than the old pixelated map, never worse.

3. The "Eigenmode" Trick: Compressing the Noise

The paper also tackles "systematic uncertainties"—these are like known flaws in the camera lens (e.g., the lens might be slightly blurry or shifted).

• The Old Way: You add a separate knob for every possible way the lens could be wrong, for every single box.
• The New Way: The authors use a technique called Eigenmode decomposition. Imagine the map has a few "fundamental shapes" (like a wave, a hill, or a dip) that represent the most common ways the data can wiggle due to noise or lens flaws.
• The Benefit: Instead of adjusting millions of knobs, you only need to adjust a handful of these "fundamental shape" knobs. It's like compressing a huge, high-definition video file into a small MP3; you keep the most important information (the shape of the signal) and throw away the redundant noise. This makes the math much faster and easier to solve.

4. The Trade-off: The "Two-Step" vs. "One-Pass"

The paper is honest about a limitation.

• The Old Method (Barlow-Beeston): This is like a "joint profile." It looks at the data and the map simultaneously, adjusting the map's wobbles in real-time as it searches for the gem. It is mathematically perfect for finding the gem when data is scarce.
• The New Method (GP Eigenmode): This is a "two-step" process. First, it builds the smooth map from the simulation. Second, it uses that fixed map to find the gem.
• The Catch: Because the map is fixed in the first step, it can't adapt perfectly to the specific noise in the final data. The paper shows that if you have very little data (scarce photos), the old method is slightly better at finding the gem because it adapts better. However, if you have lots of data (which is common in modern experiments), the difference is tiny, and the new method's speed and simplicity win out.

Summary of the Paper's Claims

• What they did: They replaced the standard "pixelated" histogram maps with smooth "Gaussian Process" maps and compressed the uncertainty into a few "eigenmodes" (fundamental shapes).
• What they proved:
  1. The new smooth maps are mathematically guaranteed to be more precise than the old pixelated maps when data is scarce.
  2. The new method can reduce the number of "wobble knobs" (parameters) from thousands to just a few dozen, making complex 3D analyses possible.
  3. The old method is still the "gold standard" for pure statistical efficiency when data is extremely rare, but the new method is practically superior for modern, complex experiments where systematic errors (like lens flaws) dominate.
• The Tool: They built this into a free software package called Histimator so other physicists can use it immediately.

In short, the paper offers a way to turn a grainy, shaky, and computationally heavy map into a smooth, stable, and efficient one, allowing physicists to search for new particles in higher dimensions without getting lost in the math.

Technical Summary

Technical Summary: Gaussian Process Eigenmodes for Statistical and Systematic Uncertainties in Template Fits

Problem Statement
Statistical inference at the Large Hadron Collider (LHC) relies on the HistFactory framework, which utilizes template histograms to model observable distributions. Uncertainties in these templates are traditionally handled via two mechanisms: per-bin Barlow–Beeston (BB) gamma factors for Monte Carlo (MC) statistical errors, and interpolation-based modifiers (e.g., histosys) for systematic shape variations. Both mechanisms scale linearly with the number of bins. This scaling becomes computationally and conceptually prohibitive for multi-dimensional analyses or when MC samples are limited. Furthermore, the BB approach treats bins as independent Poisson counts, discarding the physical smoothness of the underlying distributions. This independence leads to a proliferation of weakly constrained nuisance parameters, causing systematic undercoverage of profile likelihoods when MC statistics are poor.

Methodology
The authors propose replacing discrete histogram templates with smooth functional representations derived from Log-Gaussian Cox Process (LGCP) posteriors fitted to MC data. The methodology proceeds in three stages:

1. LGCP Modeling: MC counts are modeled as a Poisson process where the log-intensity is drawn from a Gaussian Process (GP). The posterior mode provides a smooth template, while the posterior covariance encodes correlated statistical uncertainty across bins.
2. Systematic Integration: Systematic shape variations are incorporated by generating GP fits for $\pm 1\sigma$ variation points. The difference in log-rates defines a systematic direction, which is added to the statistical covariance as a rank-1 update.
3. Eigenmode Decomposition: The combined covariance matrix (statistical + systematic) is eigendecomposed. The resulting eigenmodes form a compact basis. Truncating this basis to the leading $k$ modes replaces the full set of per-bin gamma factors and interpolation parameters with a small number of Gaussian-constrained amplitudes ($z_i$).

The authors prove that this construction contains the Barlow–Beeston formalism as a limiting case (when the GP lengthscale $\ell \to 0$) and that the GP posterior variance is strictly bounded above by the BB variance at every bin. Additionally, in the limit of negligible statistical uncertainty, the framework recovers the HistFactory InterpCode 4 interpolation.

Key Contributions

• Unified Uncertainty Basis: The paper introduces a single eigenmode basis that simultaneously encodes statistical and systematic template uncertainties, reducing the parameter space dimensionality significantly compared to the histogram approach.
• Theoretical Bounds: It is proven that the GP posterior variance is bounded by the BB variance, ensuring the method does not underestimate uncertainty. The framework is shown to recover both BB and standard HistFactory interpolation as limiting cases.
• Implementation: The method is implemented in the open-source Python package Histimator, providing an imperative API for constructing these likelihoods without dependence on the ROOT framework.
• Diagnostic Tools: The paper demonstrates how to project eigenmode pulls back to the bin level, allowing analysts to interpret results using familiar per-bin diagnostic tools.

Results
The method was validated against two benchmark experiments:

1. Experiment A (Statistically Limited): A rare-resonance search with limited MC statistics ($N_{MC}$ down to 100 events).

   • Binning Dilemma: The GP template resolved the tension between coarse binning (smearing signals) and fine binning (noisy templates). It maintained stable uncertainty quantification (8–15% posterior uncertainty) across the spectrum even when histogram bins contained fewer than 5 events.
   • Coverage: While the joint-profile BB method achieved better asymptotic efficiency in the low-statistics regime (due to adapting to data), the GP method provided continuous, usable estimates where histograms failed (empty bins). The GP method exhibited a bias-variance trade-off characteristic of two-step plug-in estimators.

2. Experiment B (Systematically Limited): A precision cross-section measurement with multiple backgrounds and four systematic sources.

   • Compression: The combined covariance required only 6–11 eigenmodes to capture 95–99% of the variance, compared to 44 nuisance parameters (40 gammas + 4 systematics) in the histogram approach. This represents a compression ratio of approximately 7:1.
   • Performance: The GP eigenmode method achieved equivalent linearity, pull width (0.96–0.99), and interval coverage (67.7–70.5% for 68% intervals) to the standard histogram approach.
   • Robustness: The reduced dimensionality led to a six-fold reduction in non-convergent fits compared to the BB method.

Significance and Claims
The paper claims that the eigenmode framework offers a principled alternative to histogram-based templates, particularly in regimes dominated by systematic uncertainties or high-dimensional phase spaces.

• Efficiency vs. Robustness: The authors explicitly acknowledge a theoretical limitation: the GP method is a "two-step plug-in" estimator, whereas Barlow–Beeston performs a "joint profile" that achieves the semiparametric efficiency bound. Consequently, in statistically limited, single-channel regimes (low MC-to-data luminosity ratio $\tau$), the BB method is structurally superior for signal extraction. However, in systematically limited regimes (high $\tau$), the efficiency loss is negligible ($<9\%$ for $\tau=10$), making the parameter compression and stability of the GP method the dominant operational advantage.
• Scalability: The method scales with the effective dimensionality of the GP kernel rather than the number of bins. For a 3D template with $20^3$ bins, the GP method requires $\sim 30$ amplitudes versus 8,000 BB gammas.
• Look-Elsewhere Effect: The smooth GP background provides an analytic covariance structure for the test-statistic field, enabling the calculation of look-elsewhere trial factors without additional Monte Carlo simulations, a capability absent in the histogram approach.

The work positions the GP eigenmode method not as a replacement for the joint-profile approach in all scenarios, but as a superior tool for managing high-dimensional systematic uncertainties and stabilizing fits in data-limited regimes where traditional histograms break down.