SSEP coupling.py

From Werner KRAUTH

Contents 1 Context 2 Python program 3 Further information 4 References

Context

This page is part of my 2024 Beg Rohu Lectures on "The second Markov chain revolution" at the Summer School "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

The present program illustrates the coupling and the phenomenon of monotone coupling in Markov-chain sampling in the example of the Symmetric Simple exclusion Process (SSEP) of NPart particles on the path graph of NSites sites without periodic boundary conditions. We start, at time t=0 with two configurations, one called ConfLow and another one, called ConfHigh, and runs the same Markov chain on both of them until they resulting configurations coincide. In the output, we check that this coupling time scales as N^3 log N. NB: There are no periodic boundary conditions, but to simplify the program, I use "phantom" vertex "-1" (containing "phantom particle "-1") and phantum vertex "NSites", with phantom particle "NPart". Only particles 0 to NPart-1, which can live on vertices 0,..., NSites-1, are real.

Python program

import math
import random
Ntrials = 100

for NPart in [8, 16, 32, 64]:
    NSites = 2 * NPart
    Coupling  = []
    for Iter in range(Ntrials):
        ConfLow = {-1: -1, NPart: NSites}; ConfHigh = {-1: -1, NPart: NSites}
        for k in range(NPart):
            ConfLow[k] = k
            ConfHigh[NPart - 1 - k] = NSites - 1 - k
        iter = 0
        while True:
            iter += 1
            Active = random.randint (0, NPart - 1)
            sigma = random.choice([-1, 1])
            if ConfLow[Active + sigma] != ConfLow[Active] + sigma: ConfLow[Active] += sigma
            if ConfHigh[Active + sigma] != ConfHigh[Active] + sigma: ConfHigh[Active] += sigma
            CLow = [ConfLow[k] for k in range(NPart)]
            CHigh = [ConfHigh[k] for k in range(NPart)]
            for k in range(NPart):
                if CLow[k]> CHigh[k]: print(Error)
            if ConfLow == ConfHigh:
                Coupling.append(iter / NPart ** 3 / math.log(NPart))
                break
    print(NPart, sum(Coupling) / len(Coupling))

 N t_coup / N^3 / log N
 8 1.2698
16 1.1993
32 1.1190
64 1.1028

Further information

This program can easily be adapted to the coupling-from-the-past approach to Markov-chain sampling introduced by Propp and Wilson in 1996, and then it allows one to obtain direct samples of the Boltzmann distribution.

References

• Propp J., Wilson D., Random Struct. Algorithms 9, 223 (1996)