Где ROOT ошибается: строгая численная проверка реализаций HistFactory

1.Введение

Модель HistFactory [1] это стандартная статистическая модель для бинированных анализов на Большом адронном коллайдере (LHC). Её математическая спецификация описывает, как систематические неопределённости (мешающие параметры) модифицируют ожидаемые выходы событий через мультипликативные факторы, морфинг форм и побиновые статистические неопределённости. Выведенная из этой модели тестовая статистика отношения профильных правдоподобий q̃(μ) [2] лежит в основе пределов исключения и значимостей открытия в ATLAS и CMS.

На практике широко используются три независимые реализации этой модели:

• ›ROOT/RooFit — исходная C++ реализация. Использует hist2workspace для построения модели и Minuit2 [4] для оптимизации. Де-факто стандарт для I/O и визуализации в CERN.
• ›pyhf — чистый Python, SLSQP. Разработан физиками ATLAS как формальный эталон спецификации HistFactory и используется в задачах реинтерпретации/комбинаций (ATL-PHYS-PUB-2019-029).
• ›NextStat — Rust, L-BFGS-B. AD в обратном режиме, опциональное GPU-ускорение (CUDA f64, Metal f32). Валидируется относительно pyhf как эталона спецификации.

Естественный вопрос: дают ли эти три реализации одни и те же физические результаты? Если нет, какая из них корректна и почему? Мы отвечаем на него строгим 3-сторонним сравнением, покрывающим полный профильный скан (31 точка по μ) на канонических фикстурах, которые задействуют все основные типы модификаторов: мультипликативные (OverallSys), побиновые (ShapeSys/StatError) и морфинг форм (HistoSys).

2.Методология

2.1 Пайплайн валидации

HistFactory XML + ROOT histograms
        |
        +---> hist2workspace -> RooFit    -> ROOT профильный скан   (C++/PyROOT)
        +---> pyhf.readxml   -> pyhf      -> pyhf профильный скан   (Python)
        +---> NextStat import -> Rust     -> NextStat профильный скан (Rust)

Профильный скан вычисляет q̃(μ) в 31 равномерно распределённой точке на отрезке μ = [0, 3]. Во всех трёх пайплайнах используется оптимизация с ограничениями, причём параметр интереса (POI) ограничен как μ ≥ 0.

2.2 Тестовая статистика

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

Как определено в Cowan и др. [2]; реализовано одинаково во всех трёх инструментах.

2.3 Фикстуры

Дополнительно мы тестируем реальные экспорты TRExFitter с числом мешающих параметров до 249.

2.4 Коды интерполяции

• ›NormSys (OverallSys): Code 4 (полином 6-го порядка)
• ›HistoSys: Code 4p (кусочно-полиномиальная)

В регрессионных тестах паритета NextStat↔pyhf эти коды фиксируются явно (pyhf: modifier_settings). В 3-стороннем сравнении ROOT↔pyhf↔NextStat мы строим все три тракта из одного и того же экспорта HistFactory (один combination.xml + те же ROOT-гистограммы), что минимизирует риск «расхождения из-за разных настроек», а не из-за численной/оптимизационной природы проблемы.

2.5 Версии и окружение

Сырые JSON-данные сканов фиксируются (snapshot) в docs/blog/assets/numerical-accuracy/data/ и рендерятся в SVG скриптом scripts/blog/generate_numerical_accuracy_plots.py. Для каждого фикстура-кейса директория данных содержит артефакты root_profile_scan.json, pyhf_profile_scan.json,nextstat_profile_scan.json и summary.json.

3.Результаты

3.1 Хорошо: ShapeSys и почти идеальное 3-стороннее согласие

Фикстура multichannel с модификаторами ShapeSys показывает «золотой стандарт»: все три реализации согласуются по q(μ) лучше чем 1e-6:

Когда все три оптимизатора сходятся, реализации дают идентичные физические результаты.

3.2 Плохо: OverallSys — ROOT переоценивает при высоком μ

Фикстура xmlimport выявляет систематическое смещение в ROOT при больших значениях силы сигнала μ:

NextStat и pyhf согласуются на уровне нескольких микродолей по q(μ) на всём скане (max |Δq| ≈ 2.4e-6). ROOT завышает q(μ) до 0.051 при μ = 3.0 (0.6% по относительной величине), причём расхождение растёт примерно линейно с расстоянием от μ̂.

3.3 Ужасно: связанный HistoSys и полный отказ ROOT

Фикстура coupled_histosys даёт наиболее драматическое расхождение:

Свободный фит ROOT возвращает status_free = -1 (Minuit2 не смог определить положительно определённую ковариационную матрицу). Условные фиты при μ ≥ 1.0 затем расходятся катастрофически.

3.4 Реальные анализы

Анализ EWK показателен: безусловный фит ROOT разошёлся настолько, что Minuit2 вернул μ̂ = 4.9e23 и NLL = 1.8e23. NextStat сходится нормально и находит μ̂ = 6.57. Поскольку μ̂ > всех точек скана, все значения q̃ оказываются равными нулю, и катастрофический отказ легко «замаскировать».

4.Качество оптимизатора на больших моделях

Помимо 3-стороннего сравнения мы отдельно анализируем поведение оптимизаторов на моделях с >100 мешающими параметрами. NextStat (L-BFGS-B) и pyhf (SLSQP) сравниваются через cross-evaluation: каждый инструмент вычисляет NLL в точке параметров, найденной другим инструментом.

4.1 Результаты cross-evaluation

На моделях с 184+ параметрами NextStat находит значения NLL на 0.01–0.08 ниже. Cross-evaluation подтверждает паритет целевой функции на уровне ~1e-13: различие обусловлено только качеством оптимизатора.

4.2 Анализ нормы градиента

Норма проецированного градиента в точке наилучшего фита это прямой индикатор качества оптимизатора. В оптимуме pyhf на воркспейсе tHu: ‖grad‖ = 4.63 (точка не стационарна). В оптимуме NextStat: ‖grad‖ = 0.020 (стационарна в пределах численной точности).

4.3 Может ли pyhf найти лучший минимум?

При тёплом старте из параметров NextStat, pyhf восстанавливает NLL = 179.404, подтверждая, что минимум реален. Однако multi-start поиск pyhf (20 случайных инициализаций) сходится лишь в 5/20 запусков, а лучший результат остаётся на 0.004 выше минимума NextStat.

4.4 Техническое объяснение

L-BFGS-B поддерживает аппроксимацию обратного гессиана ограниченной памяти, используя последние m = 10 пар (s, y), что даёт O(mn) на итерацию при наличии информативной кривизны. SLSQP может испытывать трудности на высокоразмерных правдоподобиях HistFactory с жёсткими ограничениями и коррелированными мешающими параметрами, где поверхность NLL имеет узкие «долины».

5.Производительность

Тайминги профильного скана (31 точка по μ, включая свободный фит):

Скорость достигается за счёт: компилируемого Rust (вместо Python-накладных расходов), AD в обратном режиме (вместо конечных разностей) и «горячего пути» без выделений памяти (zero-allocation) с предкомпилированным вычислением модификаторов. CPU: Apple M5 (arm64). GPU-бенчмарк на псевдоэкспериментах: NVIDIA RTX 4000, дополнительное ускорение 6.4× относительно CPU.

6.Методология валидации

6.1 Семь уровней контракта допусков

Уровни 1–6 гейтятся в CI относительно pyhf (как эталона спецификации) и внутренних детерминированных режимов NextStat; сравнение с ROOT выделено в отдельные ROOT-зависимые проверки и обычно пропускается в средах без root/hist2workspace.

6.2 Почему pyhf — эталон, а не ROOT

• ›Опубликовано в JOSS [5] после рецензирования
• ›Используется для опубликованных результатов ATLAS
• ›Поддерживается физиками ATLAS (Heinrich, Feickert, Stark)
• ›Источник формата JSON для обмена воркспейсами
• ›Тестируется против вывода ROOT в собственном CI

Спецификация HistFactory задаётся математической моделью [1], а не какой-либо программной реализацией. Валидация относительно pyhf это и есть валидация относительно спецификации.

6.3 Воспроизводимость

# Команды ниже предполагают запуск из корня репозитория NextStat.

# 3-стороннее сравнение профильного скана (ROOT vs NextStat, опционально +pyhf)
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 \
  --start 0.0 --stop 3.0 --points 31 \
  --workdir tmp/root_parity \
  --include-pyhf --keep

# Если нужно отделить «проблему оптимизатора» от возможной «разницы модели»,
# добавьте --dump-root-params: скрипт пересчитает NLL NextStat в параметрах ROOT
# и запишет cross-evaluation диагностику в summary.json.

# Агрегатор для нескольких кейсов (Apex2): пишет единый JSON-репорт + per-case run dirs
PYTHONPATH=bindings/ns-py/python ./.venv/bin/python tests/apex2_root_suite_report.py \
  --cases tests/fixtures/trex_parity_pack/cases_realistic_fourtop.json \
  --dq-atol 2e-2 \
  --mu-hat-atol 5e-2 \
  --deterministic \
  --out tmp/apex2_root_suite_realistic_fourtop.json

# Cross-evaluation / диагностика оптимизатора
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

# Перегенерировать фигуры (читает snapshot-данные)
python3 scripts/blog/generate_numerical_accuracy_plots.py

Скрипт validate_root_profile_scan.py создаёт директорию видаtmp/root_parity/run_<timestamp>/, где сохраняютсяroot_profile_scan.json, nextstat_profile_scan.json (и опционально pyhf_profile_scan.json), а также summary.json с метриками расхождения. Для быстрого объяснения, где именно расходится q(μ), используйте:tests/explain_root_vs_nextstat_profile_diff.py.

7.Обсуждение

7.1 Следствия для опубликованных результатов

• ›OverallSys bias: Переоценка 0.6% при μ = 3.0 приводит к немного более жёсткому пределу исключения. Для большинства анализов незначительно, но это систематическое смещение, растущее с расстоянием от μ̂.
• ›Отказ на связанном HistoSys: q(μ) = 41.6 против 19.0 при μ = 3.0. В асимптотическом приближении: √41.6 − √19.0 = 6.4σ − 4.4σ = 2.0σ завышения.

Это влияет только на анализы, где используются связанные модификаторы HistoSys при экстремальных значениях параметров. В большинстве анализов используется OverallSys, и там смещение ROOT невелико.

7.2 Оптимизатор важен

Cross-evaluation (§4) однозначно показывает, что ландшафт NLL идентичен, а все различия возникают из-за оптимизатора. Для <50 параметров все три оптимизатора сходятся к одному и тому же минимуму. Для больших моделей выбор оптимизатора становится физически значимым решением.

7.3 Чем NextStat отличается

• ›Точные градиенты: AD в обратном режиме даёт градиенты на уровне машинной точности, устраняя O(h²) ошибку усечения конечных разностей.
• ›Оптимизация без выделений: горячий путь NLL не делает heap-аллоков, уменьшая давление на кэш.
• ›Валидация по спецификации: Каждый CI-запуск проверяется против pyhf, не ROOT. Если ROOT не согласен, мы исследуем — но не ломаем CI ради багов ROOT.

7.4 Практические меры предосторожности

• ›Считать ненулевой статус фита или провал ковариации жёстким фейлом
• ›Сверяться с независимой реализацией (pyhf) или делать cross-evaluation NLL в точках наилучшего фита
• ›Трекать метрику сходимости (норма проецированного градиента / KKT-остатки) на высокоразмерных моделях
• ›Делать тёплый старт условных фитов вдоль μ-скана из соседних точек
• ›Фиксировать версии инструментов и толерантности фита с каждым опубликованным пределом

8.Заключение

Мы продемонстрировали через систематическое сравнение, что:

• ›NextStat и pyhf согласуются до < 1e-5 по q(μ) на всех канонических фикстурах и реальных анализах с числом параметров до 249.
• ›ROOT/RooFit демонстрирует три режима отказа: несходимость оптимизатора (связанный HistoSys), катастрофическая дивергенция (μ̂ = 4.9e23) и систематическое завышение q(μ) (OverallSys, до 0.6%).
• ›На больших моделях (>100 параметров), L-BFGS-B находит минимумы NLL на 0.01–0.08 глубже, чем SLSQP; это подтверждено паритетом целевой функции на уровне ~1e-13.
• ›NextStat быстрее в 37×–880× относительно ROOT и в 37×–75× относительно pyhf, с дополнительным GPU-ускорением 6.4×.

ROOT остаётся стандартом для I/O и визуализации, но для статистического инференса независимые реализации с современными оптимизаторами и строгой валидацией относительно математической спецификации дают более надёжные результаты.