Hamiltonian Monte Carlo

Hamiltonian Monte Carlo

1. Dirichlet Distribution

1. 개요

Hamiltonian Monte Carlo(이하 HMC)는 MCMC(Markov Chain Monte Carlo) 방법 중 하나로, 우리가 많이 알고 있는 Metropolis Hastings 알고리즘에 약간의 물리학의 개념이 합쳐진 알고리즘으로 생각하면 쉽다. Hamiltonian Dynamics(자세한 내용 생략)을 활용하여, 보다 효율적으로 공간을 탐색하여 샘플링을 하는 것이 이 알고리즘의 핵심이라고 할 수 있다. HMC는, 연속되는 두 샘플간의 correlation을 효과적으로 줄임으로써 MH 알고리즘에서 제안한 Gaussian Random Walk보다 더 효율적이고 빠르게 수렴한다.

2. 그림으로 이해하기

아래의 그림을 보자.

아래와 같은 아주 간단한 분포가 있다고 해보자. X축은 우리가 샘플링하고자 하는 대상인 \(\theta\)라고 하고, Y축은 posterior density인 \(P(\theta \mid X)\)라고 해보자.

우리는 위의 분포에서 가장 높은 y값(=\(P(\theta\mid X)\) )을 가지는 x값(=\(\theta\))에서 많은 샘플링이 이루어지기를 원한다 ( 반대로, 높이가 낮을수록 샘플이 적게끔 ). 이는, 반대로 말해서 위 그래프를 아래와 같이 위아래로 뒤집은 다음, 가장 낮은 지점에서 샘플링이 이루어지길 원하는 것과 같은 꼴이다.

위와 같은 분포를 골짜기로 생각하고, 썰매를 탄다고 생각해보자. 한번 타고 내려오면, 우리는 주로 낮은 부분에 위치하고, 아주 적은 확률로 높은 위치에 있게 될 것이다. 따라서 썰매가 더 자주 있게 되는 낮은 위치에서 샘플링이 많이 이루어지고, 썰매가 드문 확률로 있게 되는 높은 곳에서는 샘플링이 많이 이루어지지 않는다고 생각하면 쉽게 이해할 수 있을 것이다.

위의 뒤집어진 그래프에서 y축을 “에너지량”(=\(E_i\))이라고 생각해보자. 즉, 에너지가 높을 수록 (골짜기의 높은 부분에 위치할 수록), 해당 에너지량의 확률 \(P(E_i)\)는 낮을 것이다. 간단하게 식으로 예를 들면, 다음과 같이 나타낼 수 있을 것이다.

[식1] \(P(E_i) \propto e^{\frac{-E_i}{T}}\) ( 에너지가 증가할수록 해당 에너지값을 가질 확률이 낮아지는 형태 )

3. Hamiltonian, \(H\)

위의 이야기 예시처럼, 우리는 “에너지”를 사용하여 위 알고리즘에 대해 이해를 할 것이다.

고등학교 때 과학시간때 배웠 듯, 에너지는 위치에너지 + 운동에너지로 구성된다. 따라서, 쉽게 생각해자면 에너지라고 할 수 있는 Hamiltonian \(H\)의 식은 아래와 같이 나타낼 수 있다.

[식2] \(H(x,p) = U(x) + K(p)\)

• \(x\) : 물체의 “위치(position)”

• \(p\) : 물체의 “운동량(momentum)”

• \(U(x)\) : 위치 에너지
• \(K(p)\) : 운동 에너지

( 우리가 샘플링하고자 하는 값은 물체의 위치로 표현되는 \(x\)이지, \(p\)가 아니다. \(p\)는 \(x\)를 샘플링하는데에 보조 역할(auxiliary variable)이라고 생각하면 된다. )

여기서 위치 에너지인 \(U(x)\)는 우리가 삼은 target distribution \(f(x)\)에 마이너스 log를 씌운 값이다

[식3] \(U(x) = -lnf(x)\)

여기서 운동 에너지인 \(K(m)\)​는 운동량(momentum)에 비례하고, 물질의 질량(mass)에 반비례하는 것으로, 아래와 같은 식으로 나타낸다.

[식4] \(V(p) = \sum_{i=1}^{d} \frac{m^2}{2\cdot\text{mass}}\) ( 여기서 \(d\)는 벡터의 차원이다 )

( 위 식을 matrix형태로 간단히 나타내면 다음과 같다 : \(\frac{1}{2}p^TM^{-1}p\) , \(M\)은 mass)

[식1]에서 알아 봤듯, 특정 위치에 속할 확률은, 해당 위치에서 가지게 되는 에너지의 값에 아래와 같이 반비례(\(P(E_i) \propto e^{\frac{-E_i}{T}}\)) 한다. 방금 배운 Hamiltonian \(H\)를 이 식에 대입하여 정리하면 다음과 같다.

\[\begin{aligned} P(x,p) &\propto e^{\frac{-H(x,p)}{T}} \\ &= e^{\frac{-U(x)-K(p)}{T}}\\ &= e^{\frac{lnf(x)-\frac{1}{2}p^TM^{-1}p}{T}} \\ & \propto f(x) \times -\frac{1}{2}p^TM^{-1}p \\ &\equiv f(x) \times -N(p \mid 0,1)\end{aligned}\] \[P(x) = \int P(x,p)dp = \frac{1}{Z}f(x) \int N(p \mid 0,1)dp = \frac{1}{Z}f(x)\]

우리는 위 식을 통해, joint distribution $P(x,p)$에서 $x$와 $p$는 서로 상관 없음을 알 수 있다. 따라서, $x$를 샘플링하기 위해 \(P(x)\)에서 샘플링하는 방법 대신, joint distribution인 \(P(x,p)\)에서 \((x,p)\)를 샘플링하고 그냥 $p$를 버리면 된다.

4. Steps of HMC

1. \(p\)를 샘플링한다. \(p \sim N(0,1)\)

2. Leap frog steps

   • \(L\) : leap frog step의 횟쉬
   • \(\Delta t\)​ : step size

   아래의 과정을 총 \(L\)번 반복한다

   (2-1) \(p\)를 절반 update한다 : \(p \leftarrow p + \frac{1}{2}\epsilon\frac{d\;lnf(x)}{d x}\)

   ( 여기서 \(\epsilon\) 는 자유롭게 튜닝 가능하다)

   (2-2) \(x\)를 update한다 : \(x \leftarrow x + \epsilon M^{-1}p\)

   (2-3) \(p\)의 나머지 절반을 update한다 : \(p \leftarrow p + \frac{1}{2}\epsilon\frac{d\;lnf(x)}{d x}\)

3. Acceptance ratio를 구한다

   \[\alpha = min(\frac{f(x^{*})P(p^{*})}{f(x^{t-1})P(p^{t-1})},1)\]

4. Acceptance ratio에 따라 이동할지 말지를 결정한다.

   \[\theta^t = \left\{\begin{matrix} \begin{aligned} &x^{*} \;\;\;\; \text{with probability} \;\; \alpha \\ &x^{t-1} \;\;\;\;\;\;\;\ \text{otherwise} \end{aligned} \end{matrix}\right.\]