Search
Duplicate

AI/ Hamiltonian Mechanics, Hamiltonian Monte Carlo(HMC)

Hamiltonian Mechanics

Hamiltonian Mechanics는 뉴턴 역학을 재구성한 이론 체계로 에너지 보존의 관점에서 물리 시스템을 바라보며 특히 보존력이 작용하는 시스템에서 유용하게 사용된다. 해밀토니안 역학이 뉴턴 운동법칙을 대체하거나 확장하지는 않지만 시스템의 전체적인 에너지 동역학을 더 명확하게 이해하는데 도움을 준다.
에너지 환경을 굴러다니는 particle(입자)의 움직임을 θ∈RD\boldsymbol{\theta} \in \mathbb{R}^D(종종 q\bold{q}로 표기됨)와 모멘텀 v∈RD\bold{v} \in \mathbb{R}^D(종종 p\bold{p}로 표기됨)의 측면에서 특징화할 수 있다. (θ,v)(\boldsymbol{\theta},\bold{v})의 가능한 값은 phase space(위상 공간)라고 부른다.
위상 공간은 물리계의 모든 가능한 상태를 나타내는 공간으로 각 점은 시스템의 특정 상태를 나타낸다. 위상 공간의 각 점은 potential 에너지와 kinetic 에너지의 합으로 표현 가능하며, 여기서 potential 에너지는 위치 에너지, 전기, 자기, 탄성에 의한 에너지 등을 일반화한 형태이다. 실제 물리계에서 시스템의 총 에너지는 potential 에너지와 kinetic 에너지의 합으로 나타난다. 해밀토니안은 시스템의 전체 에너지(potential+kinetic)를 타나내는 함수를 말한다.
위상 공간에서 각 점에 대한 Hamiltonian 함수를 다음과 같이 정의한다.
H(θ,v)≜E(θ)+K(v)\mathcal{H}(\boldsymbol{\theta},\bold{v}) \triangleq \mathcal{E}(\boldsymbol{\theta}) + \mathcal{K}(\bold{v})
여기서 E(θ)\mathcal{E}(\boldsymbol{\theta})는 potential energy이고 K(v)\mathcal{K}(\bold{v})는 kinetic(운동) 에너지이다. 해밀토니안은 전체 에너지이다. 물리학 설정에서 포텐셜 에너지는 중력의 당김 때문이고, 모멘텀은 파티클의 움직임 때문이다. 참고로 Lagrangian mechanics은 해밀토니안과 다른 방식으로 운동 에너지를 정의하고 같은 방식으로 위치 에너지를 정의한 후에 그 운동 에너지와 위치 에너지의 차로 정의된다. 해밀토니안과 라그랑지안은 르장드르 변환으로 변환 가능하다.
통계학 설정에서는 포텐셜 에너지를 다음과 같이 취한다. 여기서 p~(θ)\tilde{p}(\boldsymbol{\theta})는 p(θ,D)p(\boldsymbol{\theta},\mathcal{D})와 같은 가능한 unnormalized 분포이다.
E(θ)=−log⁡p~(θ)\mathcal{E}(\boldsymbol{\theta}) = -\log\tilde{p}(\boldsymbol{\theta})
키네틱 에너지는 다음과 같이 정의된다. 여기서 Σ\boldsymbol{\Sigma}는 양의 정부호 행렬이고 inverse mass matrix라고 한다.
K(v)=12v⊤Σ−1v\mathcal{K}(\bold{v}) = {1\over 2} \bold{v}^\top \boldsymbol{\Sigma}^{-1}\bold{v}
안정된 궤도(orbit)은 상수 에너지를 갖는 phase 공간에서 궤적(trajectories)로 정의된다. 에너지 레벨 설정 안의 파티클의 궤적은 Hamilton’s equations이라 부르는 다음의 연속 시간 미분 방정식을 풀어서 얻을 수 있다.
dθdt=∂H∂v=∂K∂vdvdt=−∂H∂θ=−∂E∂θ\begin{aligned} {d\boldsymbol{\theta} \over dt} &= {\partial \mathcal{H} \over \partial \bold{v}} = {\partial \mathcal{K} \over \partial \bold{v}} \\ {d\bold{v}\over dt} &= -{\partial \mathcal{H} \over \partial \boldsymbol{\theta}} = -{\partial \mathcal{E} \over \partial \boldsymbol{\theta}} \end{aligned}
왜 에너지가 보존되는지 보기 위해 다음에 유의하라.
dHdt=∑i=1D[∂H∂θidθidt+∂H∂vidvidt]=∑i=1D[∂H∂θi∂H∂vi−∂H∂θi∂H∂vi]=0{d\mathcal{H} \over dt} = \sum_{i=1}^D \left[{\partial \mathcal{H} \over \partial \boldsymbol{\theta}_i}{d\boldsymbol{\theta}_i \over dt} + {\partial \mathcal{H} \over\partial \bold{v}_i}{d \bold{v}_i \over dt} \right] = \sum_{i=1}^D \left[ {\partial \mathcal{H} \over \partial \boldsymbol{\theta}_i}{\partial \mathcal{H} \over \partial \bold{v}_i} - {\partial \mathcal{H} \over \partial \boldsymbol{\theta}_i}{\partial \mathcal{H} \over \partial \bold{v}_i}\right] = 0
직관적으로 이 결과는 다음과 같이 해석할 수 있다. 행성의 궤도를 도는 위성은 모멘텀 때문에 직선으로 연속하기를 원하지만 중력 때문에 행성 쪽으로 당겨지게 되고, 이 힘이 상쇄되면 궤도는 안정적이다. 위성이 행성 쪽으로 회전하면 키네틱 에너지는 증가하지만 포텐셜 에너지는 감소한다.
어떤 시간 증가분 ss에 대한 (θ(t),v(t))(\boldsymbol{\theta}(t),\bold{v}(t))에서 (θ(t+s),v(t+s))(\boldsymbol{\theta}(t+s),\bold{v}(t+s))로 매핑하는 것은 작은 충분한 시간 단계에서 가역이다. 게다가 이 매핑은 볼륨을 보존하므로 야코비안 행렬식은 1이다. 이 사실은 나중에 이 시스템을 MCMC 알고리즘으로 전환할 때 중요하다.

Hamilton’s Equations

Euler’s method

해밀토니안 방정식의 시간 진화를 모델링하는 가장 단순한 방법은 step size η\eta라고 하는 작은 양으로 위치와 모멘텀을 동시에 업데이트 하는 것이다. 이것은 Euler’s 방법이라고 한다.
vt+1=vt+ηdvdt(θt,vt)=v(t)−η∂E(θt)∂θθt+1=θt+ηdθdt(θt,vt)=θt+η∂K(vt)∂v\begin{aligned} \bold{v}_{t+1} &= \bold{v}_t + \eta{d\bold{v} \over dt}(\boldsymbol{\theta}_t, \bold{v}_t) =\bold{v}(t) - \eta{\partial \mathcal{E}(\boldsymbol{\theta}_t) \over \partial \boldsymbol{\theta}} \\ \boldsymbol{\theta}_{t+1} &= \boldsymbol{\theta}_t + \eta{d \boldsymbol{\theta} \over dt}(\boldsymbol{\theta}_t,\bold{v}_t) = \boldsymbol{\theta}_t + \eta {\partial \mathcal{K}(\bold{v}_t) \over \partial \bold{v}} \end{aligned}

12.5.2.2 Modified Euler’s method

수정된 오일러 방법은 약간 더 정확하고 다음과 같이 작동한다. 우선 모멘텀을 업데이트하고 새로운 모멘텀을 사용하여 위치를 업데이트 한다.
vt+1=vt+ηdvdt(θt,vt)=v(t)−η∂E(θt)∂θθt+1=θt+ηdθdt(θt,vt+1)=θt+η∂K(vt+1)∂v\begin{aligned} \bold{v}_{t+1} &= \bold{v}_t + \eta{d\bold{v} \over dt}(\boldsymbol{\theta}_t, \bold{v}_t) =\bold{v}(t) - \eta{\partial \mathcal{E}(\boldsymbol{\theta}_t) \over \partial \boldsymbol{\theta}} \\ \boldsymbol{\theta}_{t+1} &= \boldsymbol{\theta}_t + \eta{d \boldsymbol{\theta} \over dt}(\boldsymbol{\theta}_t,\bold{v}_{t+1}) = \boldsymbol{\theta}_t + \eta {\partial \mathcal{K}(\bold{v}_{t+1}) \over \partial \bold{v}} \end{aligned}
불행히 이 방법의 비대칭성은 어떤 이론적인 문제의 원인이 될 수 있다.

12.5.2.3 Leapfrog integrator

Leapfrog integrator는 수정된 오일러 방법의 비대칭성을 해결하기 위한 대칭 버전으로, 우선 모멘텀을 ‘절반’ 업데이트하고 위치를 전체 업데이트하고 마지막으로 모멘텀의 나머지 ‘절반’을 업데이트 한다.
vt+1/2=vt−η2∂E(θt)∂θθt+1=θt+η∂K(vt+1/2)∂vvt+1=vt+1/2−η2∂E(θt+1)∂θ\begin{aligned} \bold{v}_{t+1/2} &= \bold{v}_t - {\eta \over 2}{\partial \mathcal{E}(\boldsymbol{\theta}_t) \over \partial \boldsymbol{\theta}} \\ \boldsymbol{\theta}_{t+1} &= \boldsymbol{\theta}_t + \eta {\partial \mathcal{K}(\bold{v}_{t+1/2}) \over \partial \bold{v}} \\ \bold{v}_{t+1} &= \bold{v}_{t+1/2} - {\eta \over 2}{\partial \mathcal{E}(\boldsymbol{\theta}_{t+1}) \over \partial \boldsymbol{\theta}} \end{aligned}
여러 leapfrog 단계를 수행하면 궤적의 시작과 끝에서 v\bold{v}의 절반 업데이트를 수행하고 그 사이에서 θ\boldsymbol{\theta}와 v\bold{v}의 전체 단계 업데이트를 번갈아 수행하는 것과 같다.

Hamiltonian Monte Carlo (HMC)

Hamiltonian Monte Carlo는 해밀토니안 역학의 개념을 적용하여 샘플링 과정에서 random walk가 아니라 목표 분포를 따라 움직이도록 한다.
해밀토니안 역학을 따라 확장된 상태 공간 (θ,v)(\boldsymbol{\theta}, \bold{v})에서 타겟 분포는 다음 형식을 갖는다.
p(θ,v)=1Zexp⁡[−H(θ,v)]=1Zexp⁡[−E(θ)−12v⊤Σv]p(\boldsymbol{\theta},\bold{v}) = {1\over Z} \exp [-\mathcal{H}(\boldsymbol{\theta},\bold{v})] = {1\over Z} \exp \left[-\mathcal{E}(\boldsymbol{\theta}) - {1\over 2} \bold{v}^\top \boldsymbol{\Sigma}\bold{v} \right]
관심 있는 잠재 변수에 대한 marginal 분포는 다음 형식을 갖는다.
p(θ)=∫p(θ,v)dv=1Zqe−E(θ)∫1Zpe−12v⊤Σvdv=1Zqe−E(θ)p(\boldsymbol{\theta}) = \int p(\boldsymbol{\theta},\bold{v})d\bold{v} = {1\over Z_q} e^{-\mathcal{E}(\boldsymbol{\theta})} \int {1\over Z_p} e^{-{1\over2}\bold{v}^\top\boldsymbol{\Sigma}\bold{v}}d\bold{v}={1\over Z_q}e^{-\mathcal{E}(\boldsymbol{\theta})}
마르코프 체인의 이전 상태가 (θt−1,vt−1)(\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})이라고 가정하자. 다음 상태를 샘플하기 위해 다음과 같이 진행한다.
1.
초기 위치를 θ0′=θt−1\boldsymbol{\theta}_0' = \boldsymbol{\theta}_{t-1}로 설정하고 새로운 랜덤 모멘텀을 샘플링한다. v0′∼N(0,Σ)\bold{v}_0' \sim \mathcal{N}(\bold{0},\boldsymbol{\Sigma}).
2.
그 다음 (θ0′,v0′)(\boldsymbol{\theta}_0',\bold{v}_0')에서 시작하여 위상 공간에서 랜덤 궤적을 초기화 하고
3.
최종 제안 상태 (θ∗,v∗)=(θL′,vL′)(\boldsymbol{\theta}^*,\bold{v}^*) = (\boldsymbol{\theta}_L',\bold{v}_L')에 도달할 때까지 LL 도약 단계에 대해 따른다.
해밀토니안 역학을 올바르게 시뮬레이션 했다면 이 절차의 시작과 끝에서 에너지는 같아야 한다. 그렇지 않으면 HMC는 diverged이고 샘플은 거절된다. 에너지가 일정하면 MH 수락 확률을 계산한다.
α=min⁡(1,p(θ∗,v∗)p(θt−1,vt−1))=min⁡(1,exp⁡[−H(θ∗,v∗)+H(θt−1,vt−1)])\begin{aligned} \alpha &= \min\left(1,{p(\boldsymbol{\theta}^*,\bold{v}^*) \over p(\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})} \right) \\&= \min(1,\exp[-\mathcal{H}(\boldsymbol{\theta}^*,\bold{v}^*) + \mathcal{H}(\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})]) \end{aligned}
마지막으로 확률 α\alpha의 (θt,vt)=(θ∗,v∗)(\boldsymbol{\theta}_t,\bold{v}_t) = (\boldsymbol{\theta}^*,\bold{v}^*) 를 설정하여 제안을 수락한다. 그렇지 않으면 (θt,vt)=(θt−1,vt−1)(\boldsymbol{\theta}_t,\bold{v}_t) = (\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})을 설정한다. 아래 pseudo code 참조.
Algorithm: Hamiltonian Monte Carlo
1.
for t=1:Tt = 1:T do
a.
랜덤 모멘텀을 생성 vt−1∼N(0,Σ)\bold{v}_{t-1} \sim \mathcal{N}(\bold{0},\boldsymbol{\Sigma})
b.
(θ0′,v0′)=(θt−1,vt−1)(\boldsymbol{\theta}_0',\bold{v}_0') = (\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})을 설정
c.
모멘텀에 대한 half step v12′=v0′−η2∇E(θ0′)\bold{v}_{1\over2}' = \bold{v}_0' -{\eta \over 2} \nabla \mathcal{E}(\boldsymbol{\theta}_0')
d.
for l=1:L−1l=1:L-1 do
i.
θl′=θl−1′+ηΣ−1vl−1/2′\boldsymbol{\theta}_l' = \boldsymbol{\theta}_{l-1}' + \eta\boldsymbol{\Sigma}^{-1}\bold{v}_{l-1/2}'
ii.
vl+1/2′=vl−1/2′−η∇E(θl′)\bold{v}_{l+1/2}' = \bold{v}_{l-1/2}' - \eta \nabla\mathcal{E}(\boldsymbol{\theta}_l')
e.
위치에 대한 full step: θL′=θL−1′+ηΣ−1vL−1/2′\boldsymbol{\theta}_L' = \boldsymbol{\theta}_{L-1}' + \eta \boldsymbol{\Sigma}^{-1}\bold{v}_{L-1/2}'
f.
모멘텀에 대한 half step: vL′=vL−1/2′−η2∇E(θL′)\bold{v}_L' = \bold{v}_{L-1/2}' - {\eta\over2}\nabla \mathcal{E}(\boldsymbol{\theta}_L')
g.
제안을 계산 (θ∗,v∗)=(θL′,vL′)(\boldsymbol{\theta}^*,\bold{v}^*) = (\boldsymbol{\theta}_L',\bold{v}_L')
h.
계산 α=min⁡(1,exp⁡[−H(θ∗,v∗)+H(θt−1,vt−1)])\alpha = \min(1,\exp[-\mathcal{H}(\boldsymbol{\theta}^*,\bold{v}^*)+\mathcal{H}(\boldsymbol{\theta}_{t-1},\bold{v}_{t-1})])
i.
확률 α\alpha로 θt=θ∗\boldsymbol{\theta}_t = \boldsymbol{\theta}^*을 설정 그렇지 않으면 θt=θt−1\boldsymbol{\theta}_t = \boldsymbol{\theta}_{t-1}

참고

Made with