HMC harmonic.py

From Werner KRAUTH

1 Context
2 Python program
3 Further information
4 References

Context

This page is part of my public lectures on Algorithms and Computations in theoretical physics at the University of Oxford (see the announcement for a syllabus). For more information, see the main page of the public lectures. It presents the hamiltonian Monte Carlo algorithm invented by Duane et al. (1988). For an in-depth discussion, see my below reference. Other programs in the same series are the Metropolis algorithm and the Lévy construction, both for the harmonic chain.

[edit]

Python program

import math, random

def U(x):
    U = 0.0
    for k in range(N):
        k_minus = (k - 1) % N
        x_minus = x[k_minus]
        if k == 0: x_minus -= L
        U += (x[k] - x_minus) ** 2 / 2.0
    return U

def grad_U(x, k):
    k_plus = (k + 1) % N
    k_minus = (k - 1) % N
    x_plus = x[k_plus]
    if k == N - 1: x_plus += L
    x_minus = x[k_minus]
    if k == 0: x_minus -= L
    return 2.0 * x[k] - x_minus - x_plus

def LeapFrog(x, p, DeltaTau, I):
    p = [p[k] - DeltaTau * grad_U(x, k) / 2.0 for k in range(N)]
    for iter in range(I):
        x = [x[k] + DeltaTau * p[k] for k in range(N)]
        if iter != I - 1: p = [p[k] - DeltaTau * grad_U(x, k) for k in range(N)]
    p = [p[k] - DeltaTau * grad_U(x, k) / 2.0 for k in range(N)]
    return x, p

N = 8
L = 2 * N
DeltaTau = 0.1
I = 16
Iter = 100000
x = [L / N  * k for k in range(N)]
U_old = U(x)
for step in range(Iter):
#
#       Choose momenta and compute kinetic energy
#
    p = [random.gauss(0.0, 1.0) for k in range(N)]
    K_old = sum([y ** 2 / 2.0 for y in p])
#
#       Do I sweeps of Leapfrog. CPU time is proportional to (2 * I + 1) * N
#
    x_new, p_new = LeapFrog(x, p, DeltaTau, I)
#
#       Do the Metropolis step
#
    U_new = U(x_new)
    K_new = sum([y ** 2 / 2.0 for y in p_new])
    if random.uniform(0.0, 1.0)  < math.exp(- U_new + U_old - K_new + K_old):
        x = x_new[:]
        U_old = U_new

[edit]