Hamiltonian Monte Carlo ทำงานอย่างไร

Question

ฉันทำกราฟิกด้านล่างเพื่ออธิบายว่าฉันเข้าใจอัลกอริทึม HMC ในปัจจุบันอย่างไร ฉันต้องการการตรวจสอบจากผู้เชี่ยวชาญเรื่องนี้หากความเข้าใจนี้ถูกหรือไม่ถูกต้อง ข้อความในสไลด์ด้านล่างถูกคัดลอกด้านล่างเพื่อความสะดวกในการเข้าถึง:

---

Hamiltonian Monte Carlo: ดาวเทียมโคจรรอบดาวเคราะห์ ยิ่งดาวเทียมอยู่ใกล้โลกมากเท่าไหร่ผลของแรงโน้มถ่วงก็จะยิ่งมากขึ้นเท่านั้น ซึ่งหมายความว่า (A) พลังงานศักย์สูงขึ้นและ (B) พลังงานจลน์ที่สูงขึ้นซึ่งจำเป็นในการรักษาวงโคจร พลังงานจลน์เดียวกันในระยะที่ไกลออกไปจากโลกจะขับดาวเทียมออกจากวงโคจร ดาวเทียมได้รับมอบหมายให้รวบรวมภาพถ่ายของพื้นที่ทางภูมิศาสตร์ที่เฉพาะเจาะจง ยิ่งดาวเทียมโคจรใกล้โลกมากเท่าไหร่ดาวเทียมก็ยิ่งเคลื่อนที่เข้าสู่วงโคจรได้เร็วขึ้นและเวลาที่โคจรผ่านพื้นที่ดังกล่าวจะยิ่งสะสมภาพถ่ายได้มากขึ้น ในทางกลับกันยิ่งดาวเทียมอยู่ห่างจากโลกมากเท่าไหร่ดาวเทียมก็ยิ่งเคลื่อนที่ช้าลงในวงโคจรเวลาที่โคจรผ่านบริเวณนั้นน้อยลงภาพถ่ายก็จะเก็บรวบรวมได้น้อยลง ในบริบทของการสุ่มตัวอย่างระยะห่างจากดาวเคราะห์หมายถึงระยะห่างจากความคาดหวังของการกระจาย พื้นที่ที่มีความเป็นไปได้ต่ำอยู่ห่างไกลจากความคาดหมาย เมื่อ "โคจรรอบความเป็นไปได้นี้" พลังงานจลน์ที่ต่ำกว่าหมายถึงเก็บตัวอย่างน้อยลงในช่วงเวลาที่กำหนดในขณะที่เมื่อโคจรรอบความเป็นไปได้ที่สูงขึ้นหมายถึงการเก็บตัวอย่างมากขึ้นในช่วงเวลาที่กำหนดเท่ากัน ในวงโคจรที่กำหนดพลังงานจลน์และศักย์ทั้งหมดจะคงที่ อย่างไรก็ตามความสัมพันธ์ระหว่างทั้งสองไม่ใช่เรื่องง่าย สมการแฮมิลตันเกี่ยวข้องกับการเปลี่ยนแปลงในสมการอื่น กล่าวคือการไล่ระดับสีของตำแหน่งเทียบกับเวลาเท่ากับโมเมนตัม และการไล่ระดับของโมเมนตัมเทียบกับเวลาเท่ากับการไล่ระดับของพลังงานศักย์เทียบกับตำแหน่ง ในการคำนวณว่าดาวเทียมจะเดินทางไปตามเส้นทางการโคจรของมันได้ไกลเพียงใดต้องใช้การรวมแบบก้าวกระโดดการอัปเดตโมเมนตัมและเวกเตอร์ตำแหน่งซ้ำ ๆ ในบริบทของการสุ่มตัวอย่างความเป็นไปได้นั้นใกล้เคียงกับระยะห่างจากโลกและการไล่ระดับของพลังงานศักย์ตามตำแหน่งคือการไล่ระดับของฟังก์ชันความหนาแน่นของความน่าจะเป็นที่เกี่ยวข้องกับพารามิเตอร์อินพุต x ข้อมูลนี้ช่วยให้สามารถสำรวจเส้นทางการโคจรรอบอินพุตต่างๆ X ซึ่งสอดคล้องกับความเป็นไปได้เดียวกัน y ที่จะสำรวจ
อย่างไรก็ตามเราไม่เพียง แต่สนใจที่จะสำรวจความเป็นไปได้เพียงอย่างเดียว แต่เราต้องสำรวจเส้นทางการโคจรหลาย ๆ ในการบรรลุเป้าหมายนี้โมเมนตัมจะต้องถูกเพิ่มแบบสุ่มโดยนำดาวเทียมเข้าใกล้หรือไกลออกไปจากโลก "โมเมนตัมคิก" แบบสุ่มเหล่านี้ทำให้มีโอกาสโคจรที่แตกต่างกัน โชคดีที่สมการแฮมิลตันช่วยให้มั่นใจได้ว่าไม่ว่าจะเป็นไปได้เท่าใดจำนวนตัวอย่างที่เก็บรวบรวมจะเป็นไปตามสัดส่วนของความเป็นไปได้ดังนั้นตัวอย่างที่เก็บรวบรวมจึงเป็นไปตามรูปร่างของการกระจายเป้าหมาย

---

คำถามของฉันคือ - นี่เป็นวิธีที่ถูกต้องในการคิดว่า Hamiltonian Monte Carlo ทำงานอย่างไร

แก้ไข:

ฉันได้ติดตั้งในโค้ดบางส่วนตามความเข้าใจของอัลกอริทึม มันใช้ได้กับ gaussian ที่มี 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()

ความคิดใด ๆ ? ฉันรู้สึกว่าฉันใกล้จะ "ได้รับ" แล้ว

Answer 1

ก่อนที่จะตอบคำถามเกี่ยวกับวิธีที่ใช้งานง่ายในการคิดเกี่ยวกับ Hamiltonian Monte Carlo ควรทำความเข้าใจกับ MCMC ปกติให้ดีที่สุด มาดูอุปมาอุปมัยของดาวเทียมกันก่อน

MCMC มีประโยชน์เมื่อคุณต้องการตัวอย่างที่ไม่มีอคติจากการแจกจ่ายที่คุณมีเฉพาะบางอย่างที่มีให้ซึ่งเป็นสัดส่วนกับ PDF แต่ไม่ใช่ PDF สิ่งนี้เกิดขึ้นใน (เช่น) การจำลองทางฟิสิกส์: PDF ถูกกำหนดโดยการแจกแจงแบบ Boltzmann, p ~ exp (-E / kT) แต่สิ่งที่คุณสามารถคำนวณได้สำหรับการกำหนดค่าใด ๆ ของระบบคือ E ไม่ใช่ p ไม่ทราบค่าคงที่ของสัดส่วนเนื่องจากอินทิกรัลของ exp (-E / kT) ในพื้นที่ทั้งหมดของการกำหนดค่าที่เป็นไปได้มักจะคำนวณยากเกินไป MCMC แก้ปัญหานั้นโดยการสุ่มเดินในลักษณะเฉพาะโดยที่ความน่าจะเป็นของการ ("ยอมรับ") แต่ละขั้นตอนจะสัมพันธ์กับอัตราส่วนของค่า p (ค่าคงที่ของสัดส่วนจะถูกยกเลิก) เมื่อเวลาผ่านไปความแตกต่างของตัวอย่างที่ได้รับการยอมรับจากการเดินสุ่มจะมาบรรจบกันเป็น PDF ที่เราต้องการโดยไม่จำเป็นต้องคำนวณ p อย่างชัดเจน

โปรดทราบว่าในข้างต้นวิธีการสุ่มขั้นตอนใด ๆ ก็ใช้ได้เท่าเทียมกันตราบใดที่ผู้เดินแบบสุ่มสามารถสำรวจพื้นที่ทั้งหมดได้ เกณฑ์การยอมรับจะรับประกันว่าตัวอย่างที่เลือกจะมาบรรจบกับ PDF จริง ในทางปฏิบัติจะใช้การแจกแจงแบบเกาส์เซียนรอบ ๆ ตัวอย่างปัจจุบัน (และซิกม่าสามารถเปลี่ยนแปลงได้เพื่อให้เศษส่วนของขั้นตอนที่ยอมรับยังคงค่อนข้างสูง) โดยหลักการแล้วจะไม่มีอะไรผิดพลาดในการทำตามขั้นตอนจากการแจกแจงแบบต่อเนื่องอื่น ๆ ("การกระจายแบบกระโดด") รอบ ๆ ตัวอย่างปัจจุบันแม้ว่าการลู่เข้าอาจช้ากว่ามาก

ตอนนี้มิล Monte Carlo ขยายอุปมาฟิสิกส์โดยเฉพาะความพยายามที่จะใช้ขั้นตอนในทิศทางซึ่งเป็นมากขึ้นมีแนวโน้มที่จะได้รับการยอมรับมากกว่าขั้นตอนของเกาส์ ขั้นตอนเหล่านี้คือสิ่งที่ผู้บูรณาการแบบก้าวกระโดดจะต้องใช้หากพยายามแก้การเคลื่อนที่ของระบบที่พลังงานศักย์เป็น E สมการการเคลื่อนที่เหล่านี้ยังรวมถึงระยะพลังงานจลน์ด้วย "มวล" (ไม่ใช่ทางกายภาพ) และ "โมเมนตัม". ขั้นตอนที่ผู้รวมระบบก้าวกระโดดใช้เวลาใน "เวลา" จะถูกส่งไปเป็นข้อเสนอไปยังอัลกอริทึม MCMC

ทำไมถึงได้ผล? Gaussian MC ก้าวไปในระยะทางเท่ากันในทุกทิศทางโดยมีความน่าจะเป็นเท่ากัน สิ่งเดียวที่ทำให้มันเอนเอียงไปสู่พื้นที่ที่มีประชากรหนาแน่นมากขึ้นของ PDF ก็คือการก้าวไปในทิศทางที่ไม่ถูกต้องมีแนวโน้มที่จะถูกปฏิเสธ พิธีกรชาวแฮมิลตันเสนอขั้นตอนทั้งในทิศทางของการไล่ระดับสี E และทิศทางของการเคลื่อนที่สะสมในขั้นตอนล่าสุด (ทิศทางและขนาดของ "โมเมนตัม") สิ่งนี้ช่วยให้สามารถสำรวจอวกาศได้เร็วขึ้นและยังมีโอกาสสูงที่จะไปถึงภูมิภาคที่มีประชากรหนาแน่นได้เร็วขึ้น

ตอนนี้อุปมาจากดาวเทียม: ฉันคิดว่านี่ไม่ใช่วิธีที่มีประโยชน์มากที่จะคิดเกี่ยวกับเรื่องนี้ ดาวเทียมเคลื่อนที่ในวงโคจรที่แน่นอน สิ่งที่คุณมีอยู่ที่นี่ค่อนข้างสุ่มเหมือนอนุภาคของก๊าซในภาชนะที่มีอนุภาคอื่น ๆ การชนกันแบบสุ่มแต่ละครั้งจะทำให้คุณมี "ก้าว"; เมื่อเวลาผ่านไปอนุภาคจะอยู่ทุกที่ในคอนเทนเนอร์โดยมีความน่าจะเป็นเท่ากัน (เนื่องจาก PDF ที่นี่มีค่าเท่ากันทุกที่ยกเว้นผนังซึ่งแสดงถึงพลังงานที่สูงมาก / เป็นศูนย์ PDF อย่างมีประสิทธิภาพ) Gaussian MCMC เปรียบเสมือนอนุภาคที่มีมวลเป็นศูนย์ที่มีประสิทธิภาพในการเดินแบบสุ่ม (หรืออนุภาคที่ไม่มีมวลเป็นศูนย์ในตัวกลางที่ค่อนข้างหนืด) มันจะไปที่นั่นผ่านการเคลื่อนที่แบบบราวเนียน แต่ไม่จำเป็นต้องเร็ว Hamiltonian MC เป็นอนุภาคที่มีมวลไม่เป็นศูนย์: มันอาจรวบรวมโมเมนตัมเพียงพอที่จะไปในทิศทางเดียวกันแม้จะมีการชนกันและบางครั้งมันอาจยิงจากปลายด้านหนึ่งไปยังอีกด้านหนึ่ง (ขึ้นอยู่กับมวลเทียบกับความถี่ / ขนาดของการชนกัน) แน่นอนว่ามันจะยังคงกระเด็นจากกำแพง แต่โดยทั่วไปแล้วมันจะมีแนวโน้มที่จะสำรวจได้เร็วขึ้น