Where ROOT Gets It Wrong

1.Introduction

The HistFactory model is the standard statistical model for binned analyses at the LHC. Its profile likelihood ratio test statistic q̃(μ) is the foundation of exclusion limits and discovery significances at ATLAS and CMS.

Three independent implementations are widely used:

• ›ROOT/RooFit — C++, Minuit2. De facto standard for I/O and visualization at CERN.
• ›pyhf — Pure Python, SLSQP. Formal reference for the HistFactory specification (ATLAS ATL-PHYS-PUB-2019-029).
• ›NextStat — Rust, L-BFGS-B. Reverse-mode AD, optional CUDA/Metal GPU acceleration.

Do these three implementations produce the same physics results? If not, which are correct? We answer with a rigorous 3-way comparison covering 31-point profile scans on canonical fixtures exercising all principal modifier types.

2.Methodology

2.1 Validation Pipeline

HistFactory XML + ROOT histograms
        |
        +---> hist2workspace -> RooFit    -> ROOT profile scan   (C++/PyROOT)
        +---> pyhf.readxml   -> pyhf      -> pyhf profile scan   (Python)
        +---> NextStat import -> Rust     -> NextStat profile scan (Rust)

The profile scan computes q̃(μ) at 31 evenly spaced points in μ = [0, 3]. All three pipelines use bounded optimization with the POI constrained to μ ≥ 0.

2.2 Test Statistic

              ⎧ 2 · [ NLL(μ, θ̂_μ) − NLL(μ̂, θ̂) ]   if μ̂ ≤ μ
q̃(μ)   =     ⎨
              ⎩ 0                                      if μ̂ > μ

As defined in Cowan et al. [2], implemented identically in all three tools.

2.3 Fixtures

We additionally test real-world TRExFitter exports with up to 249 nuisance parameters.

2.4 Interpolation Codes

• ›NormSys (OverallSys): Code 4 (polynomial, 6th order)
• ›HistoSys: Code 4p (piecewise polynomial)

This eliminates interpolation choice as a source of discrepancy.

2.5 Versions and Environment

Raw scan JSON is snapshotted under docs/blog/assets/numerical-accuracy/data/ and rendered into SVGs by scripts/blog/generate_numerical_accuracy_plots.py.

3.Results

3.1 The Good: ShapeSys — Near-Perfect 3-Way Agreement

The multichannel fixture with ShapeSys modifiers achieves the gold standard: all three implementations agree to better than 1e-6 on q(μ):

When all three optimizers converge, the implementations produce identical physics results.

3.2 The Bad: OverallSys — ROOT Overestimates at High μ

The xmlimport fixture reveals a systematic bias in ROOT at high signal strengths:

NextStat and pyhf agree to 1e-7. ROOT overestimates q(μ) by up to 0.051 at μ = 3.0 (0.6% relative), growing linearly with distance from μ̂.

3.3 The Ugly: Coupled HistoSys — ROOT Fails Completely

The coupled_histosys fixture produces the most dramatic divergence:

ROOT's free fit reports status = -1 (Minuit2 could not determine a positive-definite covariance matrix). The conditional fits at μ ≥ 1.0 diverge catastrophically.

3.4 Real-World Analyses

The EWK analysis is striking: ROOT's unconditional fit diverged so severely that Minuit2 returned μ̂ = 4.9e23 with NLL = 1.8e23. NextStat converges normally to μ̂ = 6.57. Since μ̂ > all scan points, all q̃ values happen to be zero — masking the catastrophic failure.

4.Optimizer Quality on Large Models

Beyond the 3-way comparison, we investigate optimizer behavior on models with >100 nuisance parameters. NextStat (L-BFGS-B) and pyhf (SLSQP) are compared via cross-evaluation: each tool evaluates its objective at the other's best-fit parameters.

4.1 Cross-Evaluation Results

On 184+ parameter models, NextStat finds NLL values 0.01–0.08 lower. Cross-evaluation confirms objective parity at ~1e-13: the difference is purely optimizer quality.

4.2 Gradient Norm Analysis

The projected gradient norm at the best-fit point is a direct measure of optimizer quality. At pyhf's optimum on the tHu workspace: ‖grad‖ = 4.63 (not stationary). At NextStat's optimum: ‖grad‖ = 0.020 (stationary to numerical precision).

4.3 Can pyhf Recover the Better Minimum?

When warm-started at NextStat's parameters, pyhf recovers NLL = 179.404 — confirming the minimum is genuine. However, pyhf's multi-start search (20 random inits) only converges in 5/20 runs, with the best result still 0.004 above NextStat's minimum.

4.4 Technical Explanation

L-BFGS-B maintains a limited-memory inverse Hessian approximation using the most recent m = 10 (s, y) pairs — O(mn) per iteration with effective curvature information. SLSQP can struggle on high-dimensional HistFactory likelihoods with tight bounds and correlated nuisance parameters, where the NLL surface contains narrow valleys.

5.Performance

Profile scan timing (31 μ-points, including free fit):

Speed comes from: compiled Rust (vs Python overhead), reverse-mode AD (vs finite-difference gradients), and zero-allocation hot path with pre-compiled modifier evaluation. CPU: Apple M5 (arm64). GPU toy benchmark: NVIDIA RTX 4000, additional 6.4× over CPU.

6.Validation Methodology

6.1 Seven-Tier Tolerance Contract

All tiers are CI-gated against pyhf's NumPy backend.

6.2 Why pyhf is the Reference, Not ROOT

• ›Published in JOSS [5] with peer review
• ›Used for published ATLAS results
• ›Maintained by ATLAS physicists (Heinrich, Feickert, Stark)
• ›Origin of the JSON workspace interchange format
• ›Tested against ROOT's output in its own CI

The HistFactory specification is defined by the mathematical model [1], not by any software implementation. Validating against pyhf IS validating against the specification.

6.3 Reproducibility

# 3-way profile scan comparison (ROOT + NextStat + optional pyhf scan)
PYTHONPATH=bindings/ns-py/python ./.venv/bin/python tests/validate_root_profile_scan.py \
  --histfactory-xml tests/fixtures/pyhf_xmlimport/config/example.xml \
  --rootdir tests/fixtures/pyhf_xmlimport \
  --include-pyhf --keep

# Cross-evaluation / optimizer diagnostics
PYTHONPATH=bindings/ns-py/python ./.venv/bin/python tests/diagnose_optimizer.py \
  --workspace tests/fixtures/workspace_tHu.json \
  --workspace tests/fixtures/tttt-prod_workspace.json \
  --multi-start 20

# Regenerate figures (reads snapshot data)
python3 scripts/blog/generate_numerical_accuracy_plots.py

7.Discussion

7.1 Implications for Published Results

• ›OverallSys bias: 0.6% overestimation at μ = 3.0 translates to a slightly tighter exclusion limit. For most analyses negligible, but it is a systematic bias growing with distance from μ̂.
• ›Coupled HistoSys failure: q(μ) = 41.6 vs 19.0 at μ = 3.0. In asymptotic approximation: √41.6 − √19.0 = 6.4σ − 4.4σ = 2.0σ overestimate.

This affects only analyses using coupled HistoSys modifiers at extreme parameter values. The majority of analyses use OverallSys, where ROOT's bias is small.

7.2 The Optimizer Matters

Cross-evaluation (§4) demonstrates unambiguously that the NLL landscape is identical — all differences arise from the optimizer. For <50 parameters, all three optimizers converge to the same minimum. For larger models, the choice of optimizer is a physics-relevant decision.

7.3 What NextStat Does Differently

• ›Exact gradients: Reverse-mode AD provides machine-precision gradients, eliminating O(h²) truncation error of finite differences.
• ›Zero-allocation optimization: NLL hot path allocates zero heap memory, reducing cache pressure.
• ›Specification-first validation: Every CI run gates against pyhf, not ROOT. If ROOT disagrees, we investigate — but we don't break CI to match ROOT's bugs.

7.4 Practical Safeguards

• ›Treat non-zero fit status / failed covariance as a hard failure
• ›Cross-check with an independent implementation (pyhf) or cross-evaluate NLL at best-fit points
• ›Track convergence metric (projected gradient norm / KKT residuals) on high-dimensional models
• ›Warm-start conditional fits along μ scan from neighboring points
• ›Record tool versions and fit tolerances with every published limit

8.Conclusion

We have demonstrated through systematic comparison that:

• ›NextStat and pyhf agree to < 1e-5 on q(μ) across all canonical fixtures and real-world analyses with up to 249 parameters.
• ›ROOT/RooFit exhibits three failure modes: optimizer non-convergence (coupled HistoSys), catastrophic divergence (μ̂ = 4.9e23), and systematic q(μ) overestimation (OverallSys, up to 0.6%).
• ›On large models (>100 params), L-BFGS-B finds NLL values 0.01–0.08 lower than SLSQP, verified at ~1e-13 objective parity.
• ›NextStat is 37×–880× faster than ROOT and 37×–75× faster than pyhf, with additional 6.4× GPU acceleration.

ROOT is the standard for I/O and visualization, but for statistical inference, independent implementations with modern optimizers and rigorous validation against the mathematical specification provide more reliable results.