Hamiltonian Monte Carlo는 어떻게 작동합니까?
현재 HMC 알고리즘을 이해하는 방법을 설명하기 위해 아래 그래픽을 만들었습니다. 이 이해가 정확하거나 잘못된 경우 주제 전문가의 확인을 받고 싶습니다. 아래 슬라이드의 텍스트는 쉽게 액세스 할 수 있도록 아래에 복사되어 있습니다.
Hamiltonian Monte Carlo : 위성이 행성을 공전합니다. 위성이 행성에 가까울수록 중력의 영향이 커집니다. 이것은 (A) 더 높은 위치 에너지와 (B) 궤도를 유지하는 데 필요한 더 높은 운동 에너지를 의미합니다. 행성에서 먼 거리에있는 동일한 운동 에너지는 위성을 궤도에서 방출 할 것입니다. 위성은 특정 지역의 사진을 수집하는 임무를 맡고 있습니다. 위성이 행성 궤도에 가까울수록 궤도를 더 빨리 이동할수록 더 많은 사진을 수집합니다. 반대로, 위성이 행성에서 멀어 질수록 궤도에서 더 느리게 이동할수록 해당 지역을 통과하는 시간이 줄어들수록 수집되는 사진이 줄어 듭니다. 샘플링의 맥락에서 행성으로부터의 거리는 분포 기대치로부터의 거리를 나타냅니다. 가능성이 낮은 영역은 예상과는 거리가 멀다. "이 우도 궤도를 도는"경우 운동 에너지가 낮 으면 고정 된 시간 간격 동안 수집 된 샘플 수가 적다는 것을 의미하는 반면, 더 높은 우도 궤도를 도는 경우 동일한 고정 시간 간격에서 더 많은 샘플이 수집됨을 의미합니다. 주어진 궤도에서 총 에너지, 운동 및 잠재력은 일정합니다. 그러나 둘 사이의 관계는 간단하지 않습니다. 해밀턴 방정식은 하나의 변화를 다른 것과 관련시킵니다. 즉, 시간에 대한 위치의 기울기는 운동량과 같습니다. 그리고 시간에 대한 운동량의 기울기는 위치에 대한 위치 에너지의 기울기와 같습니다. 위성이 궤도 경로를 따라 얼마나 멀리 이동했는지 계산하려면 도약 적분을 사용하여 운동량과 위치 벡터를 반복적으로 업데이트해야합니다. 샘플링의 맥락에서 가능성은 행성으로부터의 거리와 유사하며 위치에 대한 위치 에너지의 기울기는 입력 매개 변수 x에 대한 확률 밀도 함수의 기울기입니다. 이 정보를 통해 동일한 가능성 y에 해당하는 다양한 입력 X 주변의 궤도 경로를 탐색 할 수 있습니다.
그러나 우리는 단순히 하나의 가능성을 탐색하는 데 관심이있는 것이 아니라 여러 궤도 경로를 탐색해야합니다. 이를 달성하려면 운동량을 무작위로 증가시켜 위성을 행성에서 더 가깝게 또는 더 멀리 가져와야합니다. 이러한 임의의 "모멘텀 킥"을 사용하면 궤도를 돌릴 수있는 다양한 가능성이 있습니다. 다행히도 해밀턴 방정식은 우도에 관계없이 수집 된 샘플의 수가 우도에 비례하므로 수집 된 샘플이 목표 분포의 형태를 따릅니다.
제 질문은-이것이 Hamiltonian Monte Carlo의 작동 방식에 대해 생각하는 정확한 방법입니까?
편집하다:
알고리즘에 대한 이해를 바탕으로 일부 코드를 구현했습니다. mu = 0, sigma = 1 인 가우스에 대해 작동합니다. 그러나 시그마를 변경하면 깨집니다. 모든 통찰력을 주시면 감사하겠습니다.
import numpy as np
import random
import scipy.stats as st
import matplotlib.pyplot as plt
from autograd import grad
def normal(x,mu,sigma):
numerator = np.exp((-(x-mu)**2)/(2*sigma**2))
denominator = sigma * np.sqrt(2*np.pi)
return numerator/denominator
def neg_log_prob(x,mu,sigma):
num = np.exp(-1*((x-mu)**2)/2*sigma**2)
den = sigma*np.sqrt(np.pi*2)
return -1*np.log(num/den)
def HMC(mu=0.0,sigma=1.0,path_len=1,step_size=0.25,initial_position=0.0,epochs=1_000):
# setup
steps = int(path_len/step_size) -1 # path_len and step_size are tricky parameters to tune...
samples = [initial_position]
momentum_dist = st.norm(0, 1)
# generate samples
for e in range(epochs):
q0 = np.copy(samples[-1])
q1 = np.copy(q0)
p0 = momentum_dist.rvs()
p1 = np.copy(p0)
dVdQ = -1*(q0-mu)/(sigma**2) # gradient of PDF wrt position (q0) aka momentum wrt position
# leapfrog integration begin
for s in range(steps):
p1 += step_size*dVdQ/2 # as potential energy increases, kinetic energy decreases
q1 += step_size*p1 # position increases as function of momentum
p1 += step_size*dVdQ/2 # second half "leapfrog" update to momentum
# leapfrog integration end
p1 = -1*p1 #flip momentum for reversibility
#metropolis acceptance
q0_nlp = neg_log_prob(x=q0,mu=mu,sigma=sigma)
q1_nlp = neg_log_prob(x=q1,mu=mu,sigma=sigma)
p0_nlp = neg_log_prob(x=p0,mu=0,sigma=1)
p1_nlp = neg_log_prob(x=p1,mu=0,sigma=1)
# Account for negatives AND log(probabiltiies)...
target = q0_nlp - q1_nlp # P(q1)/P(q0)
adjustment = p1_nlp - p0_nlp # P(p1)/P(p0)
acceptance = target + adjustment
event = np.log(random.uniform(0,1))
if event <= acceptance:
samples.append(q1)
else:
samples.append(q0)
return samples
이제 여기에서 작동합니다.
mu, sigma = 0,1
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
하지만 시그마를 2로 변경하면 깨집니다.
# Generate samples
mu, sigma = 0,2
trial = HMC(mu=mu,sigma=sigma,path_len=2,step_size=0.25)
# What the dist should looks like
lines = np.linspace(-6,6,10_000)
normal_curve = [normal(x=l,mu=mu,sigma=sigma) for l in lines]
# Visualize
plt.plot(lines,normal_curve)
plt.hist(trial,density=True,bins=20)
plt.show()
어떤 아이디어? 나는 "그것을 얻는 것"에 가깝다고 느낀다.
답변
Hamiltonian Monte Carlo에 대해 생각하는 직관적 인 방법에 대한 질문에 답하기 전에 일반 MCMC를 확실히 파악하는 것이 가장 좋습니다. 지금은 위성의 비유를 제쳐두 자.
MCMC는 PDF 자체가 아닌 PDF에 비례 하는 사용 가능한 항목 만있는 배포에서 편향되지 않은 샘플을 원할 때 유용 합니다. 이것은 (예) 물리학 시뮬레이션에서 발생합니다. PDF는 Boltzmann 분포, p ~ exp (-E / kT)에 의해 제공되지만 시스템의 모든 구성에 대해 계산할 수있는 것은 p가 아니라 E입니다. 가능한 구성의 전체 공간에 대한 exp (-E / kT)의 적분은 일반적으로 계산하기가 너무 어렵 기 때문에 비례 상수는 알 수 없습니다. MCMC는 특정 방식으로 무작위 걷기를 수행하여이 문제를 해결합니다. 여기서 각 단계를 취할 ( "수락") 확률은 p 값의 비율 (비례 상수가 상쇄 됨)과 관련이 있습니다. 시간이 지남에 따라 랜덤 워크에서 허용 된 샘플의 분포는 명시 적으로 p를 계산할 필요없이 원하는 PDF로 수렴됩니다.
위의 경우 임의의 보행자가 전체 공간을 탐색 할 수있는 한 임의의 단계를 수행하는 모든 방법이 똑같이 유효합니다. 수락 기준은 선택한 샘플이 실제 PDF로 수렴되도록 보장합니다. 실제로 현재 샘플 주변의 가우스 분포가 사용됩니다 (시그마는 허용 된 단계의 비율이 상대적으로 높게 유지되도록 변할 수 있음). 수렴이 훨씬 느릴 수 있지만 현재 샘플 주변의 다른 연속 분포 ( "점핑 분포")에서 단계를 수행하는 데 원칙적으로 잘못된 것은 없습니다.
이제 Hamiltonian Monte Carlo 는 가우스 단계보다 받아 들여질 가능성 이 더 높은 방향으로 구체적으로 조치를 취함으로써 물리학 적 은유를 확장합니다 . 단계는 위치 에너지가 E 인 시스템의 운동을 풀려고 할 때 도약 적분기가 취할 것입니다. 이러한 운동 방정식에는 운동 에너지 용어도 포함되며, (문자 그대로 물리적이 아닌) "질량"및 "기세". 도약 적분기가 "시간"내에 취하는 단계는 MCMC 알고리즘에 제안으로 전달됩니다.
왜 이것이 작동합니까? 가우스 MC는 모든 방향에서 동일한 확률로 동일한 거리를 걷는다. PDF에서 더 밀집된 영역으로 편향되는 유일한 것은 잘못된 방향으로의 단계가 거부 될 가능성이 더 높다는 것입니다. Hamiltonian MC는 E 기울기의 방향과 최근 단계 ( "운동량"의 방향 및 크기)에서 누적 된 모션의 방향으로 단계를 제안합니다. 이를 통해 공간을 더 빠르게 탐색 할 수 있으며 인구 밀도가 높은 지역에 더 빨리 도달 할 가능성이 높아집니다.
자, 위성 비유 : 저는 이것이 그것에 대해 생각하는 데 매우 유용한 방법이 아니라고 생각합니다. 위성은 정확한 궤도로 이동합니다. 여기에있는 것은 매우 무작위 적입니다. 다른 입자가 들어있는 용기의 가스 입자와 비슷합니다. 각 무작위 충돌은 "단계"를 제공합니다. 시간이 지남에 따라 입자는 동일한 확률로 컨테이너의 모든 곳에있을 것입니다 (여기서 PDF는 매우 높은 에너지를 나타내는 벽을 제외하고 PDF가 거의 동일하기 때문에). Gaussian MCMC는 무작위 걷기 (또는 상대적으로 점성이있는 매체에서 0이 아닌 질량 입자)를 수행하는 사실상 제로 질량 입자와 같습니다. 갈색 운동을 통해 도달하지만 반드시 빠르지는 않습니다. Hamiltonian MC는 질량이 0이 아닌 입자입니다. 충돌에도 불구하고 같은 방향으로 계속 나아가기에 충분한 운동량을 모을 수 있으므로 때때로 컨테이너의 한쪽 끝에서 다른 쪽 끝으로 쏠 수 있습니다 (질량 대 주파수에 따라 다름). 충돌의 크기). 물론 여전히 벽에서 튀어 나올 것이지만 일반적으로 더 빨리 탐색하는 경향이 있습니다.