Гамильтоново Монте-Карло

Гамильтоново Монте-Карло (Hamiltonian Monte Carlo, HMC) — это метод марковских цепей Монте-Карло (MCMC) для генерации выборок из непрерывных распределений вероятностей, особенно эффективный в задачах байесовской статистики и машинного обучения. В отличие от классических методов случайных блужданий (например, алгоритма Метрополиса — Гастингса), HMC использует градиенты логарифма целевой плотности для построения траекторий, что позволяет значительно ускорить сходимость и уменьшить корреляцию между последовательными выборками. Метод основан на моделировании динамики физической системы с помощью уравнений Гамильтона, где целевое распределение интерпретируется как потенциальная энергия, а вспомогательные импульсные переменные — как кинетическая энергия.

История

Идея использования гамильтоновой динамики для выборки из распределений была впервые предложена физиками Саймоном Дуэйном, Аланом Кеннеди и Брайаном Пендлтоном в 1987 году в работе «Hybrid Monte Carlo» (HMC). Первоначально метод разрабатывался для задач квантовой хромодинамики на решётке, где требовалось эффективно интегрировать по многомерным конфигурациям полей. В 1990-е годы метод был адаптирован для статистического вывода, но его широкое распространение началось только в 2010-х годах с ростом популярности байесовского анализа и появлением автоматических библиотек дифференцирования (например, Stan, PyMC, TensorFlow Probability). Ключевой вклад в популяризацию HMC внёс Радул Нил, который в 2011 году опубликовал обзор «MCMC Using Hamiltonian Dynamics», систематизировавший теорию и практические рекомендации.

Физическая аналогия

Гамильтоново Монте-Карло опирается на формализм классической механики. Рассматривается система с позиционной переменной \( q \) (соответствует параметрам модели) и импульсной переменной \( p \) (вспомогательные переменные, не имеющие прямого статистического смысла). Полная энергия системы (гамильтониан) записывается как сумма потенциальной и кинетической энергий:

\[ H(q, p) = U(q) + K(p), \]

где \( U(q) = -\log \pi(q) \) — потенциальная энергия, связанная с целевым распределением \( \pi(q) \), а \( K(p) = \frac{1}{2} p^T M^{-1} p \) — кинетическая энергия с массовой матрицей \( M \) (обычно диагональной или единичной). Эволюция системы во времени описывается уравнениями Гамильтона:

\[ \frac{dq}{dt} = \frac{\partial H}{\partial p} = M^{-1} p, \quad \frac{dp}{dt} = -\frac{\partial H}{\partial q} = \nabla_q \log \pi(q). \]

Интегрирование этих уравнений на некоторый временной интервал \( \tau \) даёт новое состояние \( (q', p') \), которое затем принимается или отклоняется по правилу Метрополиса — Гастингса. Благодаря сохранению энергии в идеальной системе, вероятность принятия шага близка к 1, что обеспечивает высокую эффективность.

Алгоритм

Основные шаги

1. Инициализация: Задать начальное значение \( q_0 \).
2. Генерация импульса: Сэмплировать импульс \( p \sim \mathcal{N}(0, M) \) из многомерного нормального распределения.
3. Интегрирование: Применить численный интегратор (обычно метод «перепрыг-лягушка» — leapfrog) для приближённого решения уравнений Гамильтона на \( L \) шагов с шагом \( \epsilon \):

• Полушаг импульса: \( p \leftarrow p + \frac{\epsilon}{2} \nabla_q \log \pi(q) \)
• Полный шаг позиции: \( q \leftarrow q + \epsilon M^{-1} p \)
• Полушаг импульса: \( p \leftarrow p + \frac{\epsilon}{2} \nabla_q \log \pi(q) \)

1. Принятие/отклонение: Вычислить изменение гамильтониана \( \Delta H = H(q', p') - H(q, p) \) и принять новое состояние с вероятностью \( \min(1, \exp(-\Delta H)) \). Если отклонено, возвращается предыдущее состояние \( q \).
2. Повторение: Повторить шаги 2–4 для получения необходимого числа выборок.

Параметры настройки

• Шаг интегрирования \( \epsilon \): Малый шаг повышает точность, но увеличивает вычислительные затраты. Слишком большой шаг приводит к высокой частоте отклонений.
• Число шагов \( L \): Определяет длину траектории. Оптимальное значение зависит от геометрии целевого распределения. Для многомерных задач часто используется адаптивная настройка (например, в алгоритме No-U-Turn Sampler, NUTS).
• Массовая матрица \( M \): Влияет на масштабирование импульсов. В стандартной реализации \( M = I \), но для улучшения сходимости может быть оценена по апостериорной ковариации.

Преимущества и недостатки

Преимущества

• Высокая эффективность: HMC генерирует менее коррелированные выборки по сравнению с алгоритмами случайного блуждания, особенно в пространствах высокой размерности (десятки и сотни параметров).
• Использование градиентов: Метод эксплуатирует информацию о локальной геометрии целевого распределения, что позволяет быстро исследовать области с высокой плотностью.
• Сохранение энергии: Вероятность принятия шага близка к 1 при точном интегрировании, что минимизирует потери вычислительных ресурсов.

Недостатки

• Вычислительная сложность: Требует вычисления градиента логарифма целевой плотности на каждом шаге, что может быть дорого для моделей с большим числом параметров или сложными функциями правдоподобия.
• Чувствительность к настройке: Неправильный выбор \( \epsilon \) и \( L \) приводит к низкой эффективности или высокой частоте отклонений.
• Неприменимость к дискретным пространствам: HMC работает только с непрерывными переменными. Для дискретных параметров требуются модификации (например, использование вспомогательных переменных).
• Требования к гладкости: Целевое распределение должно быть дифференцируемым по параметрам.

Применение

Байесовская статистика

HMC является стандартным методом для апостериорного вывода в сложных байесовских моделях, таких как:

• Иерархические модели (например, в анализе клинических испытаний).
• Модели с нелинейными связями (например, в экологии и эпидемиологии).
• Модели с большим числом параметров (например, в нейросетях с байесовским обучением).

Машинное обучение

• Байесовские нейронные сети: HMC используется для оценки апостериорного распределения весов, что позволяет учитывать неопределённость прогнозов.
• Генеративные модели: В вариационных автокодировщиках и нормализующих потоках HMC может применяться для сэмплирования из скрытых распределений.
• Оптимизация: В некоторых задачах HMC используется как инструмент для глобальной оптимизации (например, в методе «Гамильтоновой оптимизации»).

Физика и инженерия

• Квантовая хромодинамика: Исходная область применения — интегрирование по конфигурациям калибровочных полей.
• Обратные задачи: В геофизике и медицинской томографии HMC применяется для восстановления параметров по зашумлённым данным.

Сравнение с другими методами MCMC

Реализации и программное обеспечение

• Stan: Вероятностный язык программирования, в котором HMC (в форме NUTS) является основным алгоритмом сэмплирования.
• PyMC: Библиотека Python для байесовского моделирования, поддерживает HMC и NUTS.
• TensorFlow Probability: Инструмент для вероятностного программирования на основе TensorFlow, включает реализацию HMC.
• NumPyro: Библиотека на базе JAX, обеспечивающая высокопроизводительный HMC с автоматическим дифференцированием.

Интересные факты

• Название «Hybrid Monte Carlo» отражает сочетание детерминированной гамильтоновой динамики и стохастического принятия решений.
• Алгоритм No-U-Turn Sampler (NUTS), предложенный Мэттью Хоффманом и Эндрю Гельманом в 2014 году, автоматически выбирает длину траектории \( L \), что устраняет необходимость ручной настройки.
• HMC лежит в основе многих современных систем байесовского вывода, используемых в научных исследованиях, от астрофизики до экономики.
• В 2020 году группа исследователей из Google предложила «Гамильтоново Монте-Карло с римановой метрикой», которое учитывает кривизну целевого распределения для ещё более эффективной выборки.

Источники

• Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. (1987). «Hybrid Monte Carlo». Physics Letters B, 195(2), 216–222.
• Neal, R. M. (2011). «MCMC Using Hamiltonian Dynamics». In Handbook of Markov Chain Monte Carlo (Chap. 5). CRC Press.
• Hoffman, M. D., & Gelman, A. (2014). «The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo». Journal of Machine Learning Research, 15(1), 1593–1623.
• Betancourt, M. (2017). «A Conceptual Introduction to Hamiltonian Monte Carlo». arXiv:1701.02434.
• Carpenter, B., et al. (2017). «Stan: A Probabilistic Programming Language». Journal of Statistical Software, 76(1).