Flow-Based Surrogates for High-Dimensional Likelihoods in Experimental Neutrino Physics

Technical Summary: Flow-Based Surrogates for High-Dimensional Likelihoods in Experimental Neutrino Physics

Problem Statement
Precision long-baseline neutrino experiments rely on near-detector data to constrain systematic uncertainties (nuisance parameters) affecting the unoscillated neutrino flux. These constraints are essential for accurate oscillation-parameter measurements at far detectors and for joint analyses across multiple experiments. However, the constrained likelihood functions in these scenarios are often high-dimensional (e.g., $\sim 100$ parameters) and non-Gaussian due to non-linear parameter dependencies, physical boundaries, and limited statistics.

Current propagation methods face a methodological gap:

1. Post-fit Gaussian Approximations: These are portable and closed-form but assume a quadratic log-likelihood, failing to capture asymmetric tails, curved correlations, and non-Gaussian modes, potentially introducing bias.
2. Markov Chain Monte Carlo (MCMC): While asymptotically exact and faithful to the likelihood structure, MCMC samples are not closed-form continuous functions, making them difficult to insert directly into downstream likelihoods for gradient-based optimization or efficient reuse.
3. Full Likelihood Implementations: These are accurate but require access to specific software, data, and infrastructure, rendering them impractical for portable constraint propagation.

Methodology
The authors propose using Normalizing Flows (NF) as explicit, portable surrogate models for high-dimensional systematic likelihoods. The approach involves learning an invertible map $f_\theta$ that transforms a simple base distribution into the target systematic likelihood density $p(\eta) \propto \mathcal{L}(\eta)$.

• Hybrid Architecture: To balance expressiveness and computational cost in high dimensions, the authors employ a hybrid architecture that partitions the nuisance parameters $\eta$ into two blocks:

  • $\eta_B$ (High-dimensional, approximately linear): Modeled using coupling layers (affine coupling blocks). These provide fast, analytic one-pass inverses and density evaluations, suitable for the bulk of parameters that exhibit linear behavior.
  • $\eta_A$ (Low-dimensional, highly non-linear): Modeled using conditional autoregressive Rational Quadratic Neural Spline Flows (RQ-NSF). This block captures complex non-Gaussian features. The autoregressive structure is conditioned on $\eta_B$ to preserve cross-correlations while restricting the expensive autoregressive computation to the smaller, difficult block.
  • Initialization: The flow is initialized using the post-fit Gaussian approximation (best-fit point and covariance matrix from a standard fit) to provide a strong starting proposal distribution and training samples.

• Training Strategy: The model is trained by minimizing a symmetric Kullback-Leibler (KL) divergence between the target likelihood and the flow approximation. To handle the high dimensionality and the cost of likelihood evaluations, an iterative asynchronous training scheme is used:

  • CPU workers sample from a "lagged" checkpoint of the flow to generate proposal samples.
  • These samples are evaluated against the expensive target likelihood.
  • The resulting dataset is pushed into a FIFO buffer.
  • A GPU trainer continuously optimizes the flow parameters on the aggregated dataset.
  • This process iteratively refines the proposal distribution, reducing the variance of importance weights.

Key Contributions and Results
The method was benchmarked on a realistic replica of a T2K near-detector analysis involving 110 systematic parameters, 10 of which introduce explicit non-Gaussianities via non-linear event weights.

• Fidelity in Parameter Space: The trained NF surrogate reproduces the target likelihood with high fidelity.
  • Marginal Distributions: The NF accurately captures asymmetric tails and modes displaced from the best-fit value for the 10 non-linear parameters, closely matching a reference MCMC chain (4.2 million samples). In contrast, the Gaussian approximation fails to capture these features.
  • Correlations: The NF correctly models curved correlations and distorted confidence regions in 2D projections, which the Gaussian model misses.
• Flux Predictions: When propagating systematic uncertainties to flux predictions:
  • The NF-derived uncertainty bands agree with MCMC results within $\approx \pm 1\%$.
  • The Gaussian approximation overestimates uncertainty widths by $\approx 3.5\%$ and underestimates the median event rate by $\approx 0.2\%$, indicating a bias in the Gaussian surrogate.
• Sampling Efficiency:
  • The trained NF achieves a relative effective sample size (rESS) of 98%, compared to only 5% for the initial Gaussian proposal.
  • This represents a $\sim 20$-fold improvement in sampling efficiency for importance sampling.
  • Once trained, the surrogate is computationally cheap: it can generate $\sim 38$ million samples per GPU-hour, compared to $\sim 10,000$ evaluations per CPU-hour for the original likelihood.
• Scalability: The computational cost scales primarily with the size of the non-linear block ($D_A$), not the total dimension. Even with $D_A=100$, the surrogate remains significantly faster than the original likelihood evaluation.

Significance and Claims
The paper claims that normalizing flows provide a "faithful and portable" solution for high-dimensional likelihoods in neutrino physics. By converting the fitted systematic likelihood into an explicit density model, the NF acts as a closed-form, samplable, and pointwise evaluable surrogate.

The significance lies in bridging the gap between Gaussian summaries and MCMC samples:

• Like Gaussian approximations, the NF is lightweight, portable, and does not require the original data or fitting machinery for downstream use.
• Like MCMC, it preserves the full non-Gaussian structure (skewness, curved correlations, non-Gaussian tails), thereby avoiding the biases inherent in quadratic approximations.

The authors position this method as essential for future precision oscillation measurements, particularly for joint experiment fits (e.g., T2K-NOvA) and global analyses where constraints from one stage must be rigorously propagated to another without loss of statistical information. The trained flow allows for gradient-based optimization and uncertainty propagation in downstream analyses without rerunning the upstream fit.