해밀토니안 몬테카를로

1. 개요[편집]

해밀토니안 몬테카를로 Hamiltonian (Hybrid) Monte Carlo
약칭	HMC
원 출처	Duane·Kennedy·Pendleton·Roweth (1987), 격자 QCD
필요한 정보	$\nabla \log \tilde\pi(q)$
적분기	립프로그(Störmer–Verlet)
최적 수락률	0.651 (고차원 극한)
대표 구현	Stan, PyMC, NumPyro (전부 NUTS)

확률분포를 언덕으로 바꾸고, 그 위로 공을 굴린다. 굴러간 자리가 곧 표본이다.

해밀토니안 몬테카를로(Hamiltonian Monte Carlo, HMC. 원래 이름은 Hybrid Monte Carlo)는 목표분포를 퍼텐셜 에너지로 해석하고 가짜 운동량을 붙여 해밀토니안 역학을 수치적분한 끝점을 제안으로 쓰는 MCMC 방법이다. 무작위걷기가 고차원에서 기어 다니는 근본 원인이 “방향 정보를 안 쓴다”는 데 있으므로, 기울기를 써서 목표가 실제로 뻗어 있는 방향으로 멀리 이동하겠다는 발상이다.

마르코프 연쇄 몬테카를로 문서가 상세균형·메트로폴리스-헤이스팅스 등 MCMC 일반론과 HMC의 위치를 다룬다면, 이 문서는 립프로그 적분기의 성질이 어떻게 정확성을 만드는지, 스텝 크기와 궤적 길이를 어떻게 정하는지를 파고든다. 적분기 자체의 이론은 심플렉틱 적분기와 베를레 적분에, 표본 효율의 측정은 유효표본크기에 위임한다.

2. 증강과 정준분포[편집]

목표 $\pi(q) \propto \tilde\pi(q)$ 에 대해 퍼텐셜을 $U(q) = -\log\tilde\pi(q)$ 로 두고, 보조 운동량 $p \sim \mathcal{N}(0, M)$ 을 도입한다. 확장된 공간의 해밀토니안은

$$H(q,p) = U(q) + K(p), \qquad K(p) = \tfrac{1}{2}p^\top M^{-1} p$$

이고, 결합밀도 $\propto e^{-H(q,p)} = \tilde\pi(q)\,e^{-K(p)}$ 는 $q$ 와 $p$ 가 독립인 곱 형태다. 따라서 $p$ 를 그냥 버리면(주변화) 남는 것이 정확히 $\pi(q)$ 다. 온도 1의 정준 앙상블을 인위적으로 만들어 놓고 위치 좌표만 읽는 셈이며, 정규화 상수 $Z$ 는 어차피 어디에도 등장하지 않는다.

한 번의 반복은 두 단계다.

1. 운동량 재추출. $p \sim \mathcal{N}(0,M)$ 을 새로 뽑는다. 이건 결합분포의 조건부에서 정확히 뽑는 깁스 단계라 항상 수락된다.
2. 결정론적 궤적. $(q,p)$ 에서 출발해 립프로그로 $L$ 스텝 적분하고, 끝점 $(q^*,p^*)$ 를 제안으로 삼아 $\min(1, e^{-\Delta H})$ 로 판정한다($\Delta H = H(q^*,p^*) - H(q,p)$).

1단계가 에너지를 갈아 끼우고, 2단계가 그 에너지 등위면 위를 멀리 훑는다. 두 단계 모두 결합분포를 보존하므로 연쇄 전체가 $\pi$ 를 남긴다.

3. 왜 립프로그여야 하는가[편집]

립프로그는 반스텝 운동량 갱신 → 온스텝 위치 갱신 → 반스텝 운동량 갱신이다.

$$p_{t+\epsilon/2} = p_t - \tfrac{\epsilon}{2}\nabla U(q_t), \quad q_{t+\epsilon} = q_t + \epsilon M^{-1} p_{t+\epsilon/2}, \quad p_{t+\epsilon} = p_{t+\epsilon/2} - \tfrac{\epsilon}{2}\nabla U(q_{t+\epsilon})$$

이 도식이 선택된 이유는 정확도가 아니라 두 개의 구조적 성질 때문이다.

• 심플렉틱성 → 위상공간 부피 보존. 각 부분 갱신이 전단 사상이라 야코비안이 정확히 1이다. 덕분에 메트로폴리스 수락비에 야코비안 항이 붙지 않고 $e^{-\Delta H}$ 만 남는다. 부피를 안 지키는 적분기(예: 순진한 오일러법)를 쓰면 보정항 없이는 목표분포가 통째로 틀어진다.
• 시간 가역성. 끝점의 운동량 부호를 뒤집고 같은 스텝 수를 적분하면 정확히 출발점으로 돌아온다. 이 성질 덕에 “적분 후 운동량 반전”이라는 사상이 자기 자신의 역함수(대합)가 되고, 제안이 대칭이라 상세균형이 성립한다. 실제 구현에서 운동량 부호를 뒤집는 코드가 안 보이는 것은, 바로 다음 반복에서 $p$ 를 어차피 새로 뽑고 $K(p)$ 가 우함수라 결과가 같기 때문일 뿐이다.

여기서 가장 중요한 인식 하나. 연속시간 해밀토니안 흐름은 $H$ 를 정확히 보존하므로 $\Delta H = 0$, 즉 수락률이 1이다. 실제로 관측되는 기각은 물리도 확률도 아니고 순전히 적분 오차다. 게다가 심플렉틱 적분기는 실제 $H$ 대신 그것에 $O(\epsilon^2)$ 만큼 가까운 섀도 해밀토니안을 정확히 보존하므로, 궤적을 길게 늘여도 $\Delta H$ 가 표류하지 않고 유계로 진동한다. 오일러법으로 같은 짓을 하면 에너지가 단조 증가해 $L$ 이 커질수록 수락률이 0으로 붕괴한다.

4. 스텝 크기, 질량행렬, 안정 한계[편집]

립프로그는 무조건 안정이 아니다. 조화진동자 $U(q) = \tfrac12\omega^2 q^2$, $M = 1$ 에 적용하면 한 스텝 전달행렬의 고윳값이 단위원 위에 남을 조건이 정확히

$$\epsilon\,\omega < 2$$

다. 이 선을 넘는 순간 고윳값이 실수로 갈라져 $|p|,|q|$ 가 지수적으로 폭발하고, $\Delta H$ 가 수백~수천으로 튀면서 수락률이 서서히가 아니라 절벽처럼 0으로 떨어진다. HMC 사용자가 겪는 “스텝을 조금 키웠더니 갑자기 아무것도 수락되지 않는다”는 현상의 정체가 이것이다.¹

다차원에서는 가장 좁은 방향이 $\epsilon$ 을 지배한다. 목표 공분산 $\Sigma$ 의 최소 고윳값 방향이 $\omega_{\max}$ 를 정하므로, 조건수가 나쁘면 넓은 방향은 손도 못 대고 좁은 방향 때문에 $\epsilon$ 을 줄여야 한다. 해법이 질량행렬을 전처리기로 쓰는 것이다. $M^{-1} \approx \Sigma$ 로 잡으면(Stan이 워밍업에서 하는 일이 바로 이 추정이다) 좌표가 표준화된 것과 같아져 모든 방향의 $\omega$ 가 비슷해진다. 대각 근사가 기본값이고, 차원이 작으면 조밀 행렬도 쓴다.

질량행렬 품질의 무료 진단도 있다. 등분배 정리에 의해 정준분포에서 $\mathbb{E}[K(p)] = d/2$ 여야 하므로, 실측 운동에너지 평균이 이 값에서 크게 벗어나면 $M$ 이 목표 기하와 안 맞는다는 신호다.

5. 궤적 길이와 NUTS[편집]

$\epsilon$ 이 정해져도 몇 스텝을 갈지($L$)가 남는다. 너무 짧으면 결국 무작위걷기고, 너무 길면 궤적이 U턴을 돌아 출발점 근처로 되돌아와 계산만 버린다. 조화진동자에서 $L\epsilon$ 을 주기의 정수배로 잡으면 제자리로 돌아오는 최악의 조합이 되는데, 이걸 사람이 손으로 피하는 건 무리다.

NUTS(No-U-Turn Sampler, Hoffman & Gelman 2014)가 이 문제를 자동화했다. 궤적을 앞뒤 무작위 방향으로 배로 늘려 가며 이진 트리를 만들고, 부분궤적의 양 끝 $(q^-,p^-)$, $(q^+,p^+)$ 에 대해

$$(q^+ - q^-)\cdot p^- < 0 \quad \text{또는} \quad (q^+ - q^-)\cdot p^+ < 0$$

이 성립하면 “되돌아오기 시작했다”고 판정해 확장을 멈춘다. 이렇게 만든 궤적 위의 점들 중 하나를 다항 가중으로 뽑으면 가역성과 상세균형이 유지된다. 스텝 크기는 목표 수락률(기본 0.8)을 맞추도록 워밍업 동안 이중 평균법으로 적응시키고, 본 표본추출에서는 고정한다 — 적응을 끝까지 켜 두면 마르코프성이 깨지기 때문이다.

핵심 진단은 발산 전이(divergent transition)다. 궤적 도중 $\Delta H$ 가 임계값을 넘으면 그 반복을 발산으로 표시하는데, 이는 대개 목표 기하에 곡률이 극심한 영역(계층 모형의 깔때기)이 있어 균일한 $\epsilon$ 으로는 못 넘는다는 뜻이다. 발산이 몇 개라도 있으면 그 표본은 편향돼 있다고 봐야 하며, $\epsilon$ 을 줄이는 것보다 비중심 재매개변수화로 기하를 펴는 쪽이 근본 해결이다.²

6. 스케일링과 한계[편집]

고차원 극한에서의 비교가 HMC의 존재 이유를 요약한다. 수락률을 일정하게 유지하려면 스텝 크기를 무작위걷기 MH는 $\epsilon \propto d^{-1/2}$, MALA는 $d^{-1/6}$, HMC는 $d^{-1/4}$ 로 줄여야 하고, 그 결과 독립표본 하나를 얻는 비용이 무작위걷기는 $O(d)$, HMC는 $O(d^{1/4})$ 로 갈린다. 최적 수락률도 각각 0.234 / 0.574 / 0.651이다. 차원 1000짜리 계층 모형에서 두 방법의 실행 시간 차이가 서너 자릿수로 벌어지는 근거가 이 지수다.

그렇다고 만능은 아니다.

• 이산 파라미터를 못 다룬다. 기울기가 없으니 원리적으로 불가능하다. 실무 처방은 이산 변수를 해석적으로 주변화하는 것이며, 은닉 마르코프 모형의 전방 알고리즘으로 상태열을 합쳐 없애는 것이 전형적인 예다.
• 다봉 분포는 여전히 못 넘는다. 궤적은 에너지 등위면 위를 도는 것이라 골짜기를 뛰어넘지 못한다. 이 문제는 병렬 템퍼링이나 담금질 모사 계열의 영역이다.
• 기울기가 비싸다. 반복마다 $\nabla U$ 를 수십 번 계산하므로 자동 미분이 사실상 전제 조건이고, 미분 불가능한 솔버가 끼어 있으면 곧장 막힌다.

변형도 많다. 위치에 따라 계량을 바꾸는 리만 다양체 몬테카를로(RMHMC)는 피셔 정보를 질량행렬로 써서 깔때기 기하를 자동으로 펴 주지만, 적분기가 음함수가 되어 스텝마다 부동점 반복이 필요하다. 대규모 데이터에서 미니배치 기울기를 쓰는 확률적 기울기 HMC는 잡음이 에너지를 끌어올리므로 마찰항을 넣어 랑주뱅 동역학 꼴로 만들어야 정상분포가 유지된다.

그리고 이 알고리즘의 고향이 통계가 아니라는 점은 짚고 갈 만하다. 1987년 격자 QCD에서 게이지장 배위를 뽑기 위해 만들어진 “Hybrid Monte Carlo”가 원형이고, 통계 쪽으로 옮겨 온 것은 1990년대 닐(R. Neal)의 작업이다.³ 이징 모형 같은 격자계에서 국소 갱신이 임계점 근처에서 임계 감속에 걸리는 문제를 우회하려던 물리학자들의 도구가, 30년 뒤 베이즈 통계 소프트웨어의 기본 엔진이 됐다.⁴

7. 관련 문서[편집]

• 마르코프 연쇄 몬테카를로 · 몬테카를로 방법
• 심플렉틱 적분기 · 베를레 적분 · 해밀토니안 역학
• 유효표본크기 · 중요도 표본추출 · 기각표본추출
• 랑주뱅 동역학 · 병렬 템퍼링
• 정준 앙상블 · 등분배 정리
• 자동 미분 · 변분 추론
• 불확실성 정량화 · 통계

8. Footnotes[편집]

1. 이게 무작위걷기 메트로폴리스와 결정적으로 다른 점이다. 저쪽은 제안 폭을 키우면 수락률이 매끄럽게 내려가서 조율이 쉽지만, HMC는 안정 한계 직전까지 수락률이 멀쩡하다가 한 발 넘는 순간 낭떠러지다. 그래서 손으로 조율하지 말고 이중 평균법에 맡기라는 것. ↩

2. 발산 몇 개쯤 무시하고 결과를 쓰고 싶은 유혹이 늘 있는데, 발산은 “표본이 조금 부족하다”가 아니라 “가지 못한 영역이 있다”는 신호라서 표본을 늘려도 해결되지 않는다. 발산 지점만 따로 산점도에 찍어 보면 대개 깔때기 목 부분에 예쁘게 몰려 있다. ↩

3. 닐의 2011년 핸드북 챕터 MCMC using Hamiltonian dynamics는 지금도 이 분야 최고의 입문서로 꼽힌다. 물리 비유와 통계적 엄밀성의 배합비가 절묘하다. ↩

4. 그래서 HMC 논문을 읽다 보면 통계 용어와 물리 용어가 한 문장 안에서 섞인다. 사후분포와 퍼텐셜, 온도와 사전분포 강도, 분배함수와 주변가능도가 같은 것을 가리키는데 부르는 이름만 다르다. 두 동네 사람이 같은 세미나에 앉으면 30분쯤 서로 다른 언어로 같은 말을 하다가 뒤늦게 알아차리는 광경이 벌어진다. ↩