Werner KRAUTH - New pages [en]

Tartero Shiratani Krauth 2026

Werner — Thu, 02 Apr 2026 15:56:11 GMT

Summary:

'''G. Tartero, S. Shiratani, W. Krauth''' ''''' Lifting the fog - a case for non-reversible "lifted" Markov chains '''''' ''' arXiv2603.16855 (2026) '''''

'''Abstract'''
Phase transitions appear all over science, and are familiar from everyday life, as water boiling, sugar melting into caramel or as nematic molecules turning smectic in liquid-crystal displays. The dynamics of phase transitions can be extremely slow, as for example when fog in winter does not lift, that is when the coarsening takes much time from many tiny water droplets to fewer but larger rain drops that feel the pull of gravity. The dynamics of phase transitions is relevant also for the performance of computer algorithms. In the ubiquitous Metropolis Monte Carlo algorithm, the mixing dynamics towards equilibrium leads towards the solution of a sampling problem. It is governed by the same reversibility and detailed-balance principles as the overdamped physical dynamics of fog. For the phase-separated Lennard-Jones system, we describe here how the coarsening dynamics of non-reversible "lifted" variants of the Metropolis algorithm proceeds on much faster time scales, with the microscopic non-reversibility translating into large-scale relative motion of droplets that is impossible under the Ostwald-ripening condition of reversibility. A density-displacement coupling moves droplets relative to each other through a lensing effect. Efficient implementations of the long-range Metropolis algorithm and its non-reversible lifting (event-chain Monte Carlo) allow us to show that, in consequence, the coarsening growth exponent is larger under lifting. For large system sizes, the computing problem is thus solved infinitely faster than before, with the outcome strictly unchanged with respect to the Metropolis algorithm. We also discuss the larger setting of our findings, namely that "lifted" non-reversible algorithms can be set up for generic reversible sampling methods, with applications going much beyond our example of lifting fog.

[http://arxiv.org/pdf/2603.16855 Electronic version (from arXiv)]

Sathwik Krauth 2025

Werner — Tue, 25 Nov 2025 01:15:06 GMT

Summary:

''K. Sathwik, W. Krauth'' ''''' Thermodynamic integration, fermion sign problem, and real-space renormalization'''''' ''' arXiv2511.12794 (2025) '''''

'''Abstract'''
We reconsider real-space renormalization for the two-dimensional Ising model, following the path traced out by Wilson in Sect. VI of his 1975 Reviews of Modern Physics. In that reference, Wilson considerably extended the Kadanoff decimation procedure towards a possibly rigorous construction of a real-space scale-invariant hamiltonian. Wilson’s construction has, to the best of our knowledge, never been fully understood and
thus neither reproduced nor generalized. In the present work, we use Monte Carlo sampling in combination with thermodynamic integration in order to retrace Wilson’s computation for a real-space renormalization with a number of terms in the hamiltonian. We elaborate on the connection of real-space renormalization with the fermion sign problem and discuss to which extent our Monte Carlo procedure actually implements Wilson’s
program from half a century ago.

[http://arxiv.org/pdf/2511.12794 Electronic version (from arXiv)]

Massoulie et al 2025

Werner — Sun, 16 Mar 2025 11:51:53 GMT

Summary:

'''Brune Massoulié, Clément Erignoux, Cristina Toninelli, Werner Krauth''' '''''Velocity trapping in the lifted TASEP and the true self-avoiding random walk''''' ''' Physical Review Letters 135, 127102 (2025)'''

'''Abstract'''
We discuss non-reversible Markov-chain Monte Carlo algorithms that, for particle systems, rigorously sample the positional Boltzmann distribution and that have faster than physical dynamics. These algorithms all feature a non-thermal velocity distribution. They are exemplified by the lifted TASEP (totally asymmetric simple exclusion process), a one-dimensional lattice reduction of event-chain Monte Carlo. We analyze its dynamics in terms of a velocity trapping that arises from correlations between the local density and the particle velocities. This allows us to formulate a conjecture for its out-of-equilibrium mixing time scale, and to rationalize its equilibrium superdiffusive time scale. Both scales are faster than for the (unlifted) TASEP. They are further justified by our analysis of the lifted TASEP in terms of many-particle realizations of true self-avoiding random walks. We discuss velocity trapping beyond the case of one-dimensional lattice models and in more than one physical dimensions. Possible applications beyond physics are pointed out.

[https://doi.org/10.1103/mqdr-x95j Published version (subscription needed)]

[http://arxiv.org/pdf/2503.10575 Electronic version (from arXiv)]

HMC harmonic.py

Werner — Wed, 19 Feb 2025 18:18:16 GMT

Summary:

==Context==
This page is part of my public lectures on [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures Algorithms and Computations in theoretical physics] at the University of Oxford (see the announcement for a syllabus). For more information, see [[Oxford_Lectures_2025|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_harmonic.py|Metropolis algorithm]] and the [[Levy_harmonic.py|Lévy construction]], both for the harmonic chain.

==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

==Further information==

==References==
Krauth, W. [http://arxiv.org/pdf/2411.11690 Hamiltonian Monte Carlo vs. event-chain Monte Carlo: an appraisal of sampling strategies beyond the diffusive regime]. ArXiv: 2411.11690
[[Category:Python]] [[Category:Oxford2025]]

Metropolis harmonic.py

Werner — Wed, 19 Feb 2025 18:13:37 GMT

Summary:

==Context==
This page is part of my public lectures on [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures Algorithms and Computations in theoretical physics] at the University of Oxford (see the announcement for a syllabus). For more information, see [[Oxford_Lectures_2025|the main page of the public lectures]]. It presents the Metropolis Monte Carlo algorithm invented by Metropolis et al. (1953). For an in-depth discussion, see my below reference. Other programs in the same series are the [[HMC_harmonic.py|Hamiltonian Monte Carlo]] algorithm and the [[Levy_harmonic.py|Lévy construction]], both for the harmonic model.

==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

N = 8
L = 16
delta = 1.0
x = [L / N * k for k in range(N)]
Iter = 1000000
for iter in range(Iter):
k = random.randint(0, N - 1) # random slice
k_plus = (k + 1) % N
k_minus = (k - 1) % N
x_plus = x[k_plus]
if k == N - 1: x_plus += L
if k_plus == N: x_plus = x[0] + L
x_minus = x[k_minus]
if k == 0: x_minus -= L
x_new = x[k] + random.uniform(-delta, delta) # new position at slice k
U_k = ((x_plus - x[k]) ** 2 + (x[k] - x_minus) ** 2) / 2.0
U_kp = ((x_plus - x_new) ** 2 + (x_new - x_minus) ** 2) / 2.0
if random.uniform(0.0, 1.0) < math.exp(- (U_kp - U_k)): x[k] = x_new

==Further information==

==References==
Krauth, W. [http://arxiv.org/pdf/2411.11690 Hamiltonian Monte Carlo vs. event-chain Monte Carlo: an appraisal of sampling strategies beyond the diffusive regime]. ArXiv: 2411.11690
[[Category:Python]] [[Category:Oxford2025]]

Levy harmonic.py

Werner — Wed, 19 Feb 2025 18:11:05 GMT

Summary:

Oxford Lectures 2025

Werner — Tue, 21 Jan 2025 13:01:59 GMT

Summary:

My [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures Public Lectures at the University of Oxford (UK)], entitled '''Algorithms and computations in theoretical physics''', ran from 21 January 2025 through 11 March 2025 and are repeated (and adapted) from 20 January 2026 through 10 March 2026.

===Lecture 1: 20 January 2026 (21 January 2025)===
We start our parallel exploration of physics and of computing with the concept of sampling, the process of producing examples (“samples”) of a probability distribution. In week 1, we consider “direct” sampling (the examples are obtained directly) and, among the many connections to physics, will come across the Maxwell distribution. In 1859, it marked the beginning of the field of statistical physics.

[http://www.lps.ens.fr/%7Ekrauth/images/b/b4/WK_Lecture1_Oxford2025.pdf Here are the notes for the first lecture (20 January 2026 (21 January 2025)). ]

===Lecture 2: 27 January 2026 (28 January 2025)===
We start with a Special Topic A, where we supplement the first Lecture with two all-important special topics that
explore the meaning of convergence in statistics, and the fundamental usefulness
of statistical reasoning. We can discuss them here, in the direct-sampling
framework, but they are more generally relevant. The strong law of large
numbers, for example, will turn into the famous ergodic theorem for Markov
chains.

[http://www.lps.ens.fr/%7Ekrauth/images/6/69/WK_Lecture2Spec_Oxford2025.pdf Here are the notes for Special Topic A (28 January 2025). ]

Then, in this second lecture proper, we consider Markov-chain sampling, from the adults' game
on the Monte Carlo heliport to modern ideas on non-reversibility. We discuss
a lot of theory, but also six ten-line pseudo-code algorithms, none of them
approximate, and all of them as intricate as they are short.

[http://www.lps.ens.fr/%7Ekrauth/images/9/9a/WK_Lecture2_Oxford2025.pdf Here are the notes for Lecture 2.]

===Lecture 3: 03 February 2025 (04 February 2025)===
In this third lecture, we consider Markov-chain sampling in an abstract setting
in one dimension. We discuss some theory, but also seven ten-line pseudo-code
algorithms, none of them approximate, and all of them as intricate as they are short.
At the end, we discuss the foundations of statistical mechanics, as seen in a one-
dimensional example.

[http://www.lps.ens.fr/%7Ekrauth/images/8/80/WK_Lecture3_Oxford2025.pdf Here are the notes for Lecture 3.]

* Here is [[Metropolis X2X4.py|the Metropolis algorithm for a single particle in an anharmonic potential]]
* Here is [[Factor Metropolis X2X4.py|the factorized Metropolis algorithm for a single particle in an anharmonic potential]]
* Here is [[ZigZag X2X4.py|the Zig-zag algorithm for a single particle in an anharmonic potential]]

===Lecture 4: 10 February 2026 (17 February 2025)===
In this fourth lecture, we treat classical many-particle systems that interact with hard-sphere
potentials, a case of importance for the foundations of statistical mechanics but also for its applications. We move from Newtonian mechanics to Boltzmann mechanics and from classical mechanics to statistical mechanics, in a way that is again entirely
example-based. We then treat the case of one-dimensional hard spheres and discover the Asakura--Oosawa interaction, the fifth force in nature. We will end with a discussion of the double nature of pressure: kinematic and thermodynamic.

[http://www.lps.ens.fr/%7Ekrauth/images/7/72/WK_Lecture4_Oxford2025.pdf Here are the notes for Lecture 4.]

* Here is [[Event_disks_box.py| the event-driven algorithm for four hard disks in a square box]], in a few dozen lines of code.

===Lecture 5: 17 February 2026 (18 February 2025)===
In this fifth lecture, we treat classical many-particle systems that interact with general interactions, in the example of the harmonic chain.

Lecture notes for lecture 5 are forthcoming, but most of the material can be found [http://arxiv.org/pdf/2411.11690 here].

* Here is [[Levy_harmonic.py| the Lévy-construction algorithm]] for the harmonic chain, in nine lines of code. It directly samples configurations.
* Here is [[Metropolis_harmonic.py| the Metropolis algorithm]] for the harmonic chain, in two dozen lines of code. It samples configurations using the algorithm developed by Metropolis et al. in (1953).
* Here is [[HMC_harmonic.py| the Hamiltonian Monte Carlo algorithm]] for the harmonic chain, in three dozen lines of code. It implements the algorithm developed by Duane et al. (1988).

===Lecture 6: 24 February 2026 (25 February 2025)===
This sixth lecture is the first of two lectures on quantum statistical mechanics. Starting from the quantum harmonic oscillator, we introduce to density matrices and
the Feynman path integral. The computational techniques that we concentrate on during these two weeks have widespread use in statistics and physics,
and the fundamental matrix-squaring property of density matrices corresponds to the convolution of probability densities. The Lévy algorithm for the sampling of
path integrals, on the other hand, is mathematically equivalent to what we discussed, in a totally different context, in Lecture 5.

[http://www.lps.ens.fr/%7Ekrauth/images/d/d3/WK_Lecture6_Oxford2025.pdf Here are the notes for Lecture 6.]

===Lecture 7: 03 March 2026 (04 March 2025)===
This is the second of two lectures on computational quantum statistical mechanics. Building on the quantum harmonic oscillator, the density matrices and the Feynman path integral introduced in the sixth lecture, we will discuss a direct-sampling quantum Monte Carlo algorithm for bosons that illustrates the phenomenon of Bose–Einstein condensation. Before doing that, however, we must discuss quantum statistics of indistinguishable particles and set up a naive program against which to check our state-of-the-art simulation code.

[http://www.lps.ens.fr/%7Ekrauth/images/f/f5/WK_Lecture7_Oxford2025.pdf Here are the notes for Lecture 7.]

* Here is [[Two_cycles.py| Two-cycles.py]], an intricate seven-line algorithm for sampling random permutations... but only those without cycles longer than 2!
* Here is [[Direct_N_bosons.py| Direct-bosons.py,]] a three-dozen-line direct-sampling algorithm for ideal bosons in a three-dimensional harmonic trap.

===Lecture 8: 10 March 2026 (11 March 2025)===
The final lecture in the set is on the classical Ising model. We go from basic (and not so basic) enumerations to the high-temperature expansion of the Ising model, touch on its exact solution by Kac and Ward (1952) and then treat sampling methods: Metropolis, heatbath, and perfect sampling, to finish with the Wolff (1989) cluster algorithm.

[http://www.lps.ens.fr/%7Ekrauth/images/5/52/WK_Lecture8_Oxford2025.pdf Here are the notes for Lecture 8.]

* Here is [[Enumerate ising dos.py| Enumerate-ising.py]] that uses the Gray code enumeration to compute the density of states of the Ising model.
* Here is [[Markov_ising.py| Markov-ising.py]] that uses the Metropolis algorithm to sample the Boltzmann distribution of the Ising model.
* Here is [[Heatbath_ising.py| Heatbath-ising.py]] that uses the Heatbath algorithm for the above sampling problem.
* Here is [[Cluster_ising.py| CLuster-ising.py]] that uses the Wolff cluster algorithm, again for sampling the Boltzmann distribution of the Ising model.

Krauth 2024 HMC ECMC

Werner — Tue, 26 Nov 2024 00:11:18 GMT

Summary:

''W. Krauth'' '''''Hamiltonian Monte Carlo vs. event-chain Monte Carlo: an appraisal of sampling strategies beyond the diffusive regime '''''
arXiv:2411.11690 (2024)

'''Abstract'''
We discuss Hamiltonian Monte Carlo (HMC) and event-chain Monte Carlo (ECMC)
for the one-dimensional chain of particles with harmonic interactions and benchmark them
against local reversible Metropolis algorithms. While HMC achieves considerable speedup
with respect to local reversible Monte Carlo algorithms, its autocorrelation functions of
global observables such as the structure factor have slower scaling with system size than for
ECMC, a lifted non-reversible Markov chain. This can be traced to the dependence of ECMC
on a parameter of the harmonic energy, the equilibrium distance, which drops out when
energy differences or gradients are evaluated. We review the recent literature and provide
pseudocodes and Python programs. We finally discuss related models and generalizations
beyond one-dimensional particle systems.

[http://arxiv.org/pdf/2411.11690 Electronic version (from arXiv)]

Hukushima Krauth 2024

Werner — Thu, 10 Oct 2024 20:25:14 GMT

Summary:

'''K. Hukushima, W. Krauth''' '''''Damage spreading and coupling in spin glasses and hard spheres ''''' '''Frontiers in Physics 12, 1507250 (2025)'''

'''Abstract'''
We study the connection between damage spreading, a phenomenon long discussed in the physics literature, and the coupling of Markov chains, a technique used to bound the mixing time. We discuss in parallel the Edwards-Anderson spin glass model and the hard-disk system, focusing on how coupling provides insights into the region of fast coupling within the paramagnetic and liquid phases. We also work out the connection between path coupling and damage spreading. Numerically, the scaling analysis of the mean coupling time determines a critical point between fast and slow couplings. The exact relationship between fast coupling and disordered phases has not been established rigorously, but we suggest that it will ultimately enhance our understanding of phase behavior in disordered systems.

[http://arxiv.org/pdf/2410.05544 Electronic version (from arXiv)]

[https://doi.org/10.3389/fphy.2024.1507250 Published version (open access)]

Essler Krauth 2023

Werner — Thu, 18 Jul 2024 15:10:05 GMT

Summary:

'''F. H. L Essler, W. Krauth''' '''''Lifted TASEP: a Bethe ansatz integrable paradigm for non-reversible Markov chains ''''' ''' Phys. Rev. X 14, 041035'''

'''Popular Summary'''
Markov-chain Monte Carlo (MCMC) algorithms formulate the sampling problem for complex probability distributions as a simulation of fictitious physical systems in equilibrium, where all motion is diffusive and time reversible. But nonreversible algorithms can, in principle, sample distributions much more efficiently. In recent years, a class of “lifted” Markov chains has implemented this idea in practice, but the resulting algorithms are extremely difficult to analyze. In this work, we introduce an exactly solvable paradigm for nonreversible Markov chains.

Our paradigm, which we term the lifted totally asymmetric simple exclusion process (TASEP), describes a particular type of nonreversible dynamics for particles on a one-dimensional lattice. We show that this dynamics allows for polynomial speedups in particle number compared to the famous Metropolis MCMC algorithm. The lifted-TASEP dynamics is, in fact, faster than that of any other known class of models. To arrive at our conclusions, we combine exact methods from the theory of integrable models with extensive numerical simulations. In particular, we prove that the lifted TASEP is integrable and determine the scaling of its relaxation and mixing times with system size.

Our work opens the door to obtaining mathematically rigorous results for speedups of nonreversible MCMC algorithms, and more generally, of lifted Markov chains arising in interacting many-particle systems.

[http://arxiv.org/pdf/2306.13059 Electronic version (from arXiv)]

[https://doi.org/10.1103/PhysRevX.14.041035 Published version (open source)]

==Further context==
The Lifted TASEP is discussed in my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution", and a sample Python program can be found [[LiftedTASEPCompact.py|here]].

My [[Massoulie_et_al_2025|2025 paper with B. Massoulié, C. Erignoux and C. Toninelli]] provides a deeper theoretical analysis.

Hard spheres coupling.py

Werner — Wed, 12 Jun 2024 21:52:42 GMT

Summary:

Ising coupling.py

Werner — Wed, 12 Jun 2024 20:56:14 GMT

Summary:

Stopping circle.py

Werner — Wed, 12 Jun 2024 20:10:04 GMT

Summary: /* Further Information */

SSEPCoupling.py

Werner — Wed, 12 Jun 2024 19:59:51 GMT

Summary: /* Further information */

Factor ZigZag X2X4.py

Werner — Tue, 11 Jun 2024 06:50:49 GMT

Summary: /* References */

Bounded Lifted Metropolis X2X4.py

Werner — Tue, 11 Jun 2024 06:39:55 GMT

Summary: /* Python program */

Factor Metropolis X2X4 patch.py

Werner — Mon, 10 Jun 2024 23:00:54 GMT

Summary:

ZigZag X2X4.py

Werner — Mon, 10 Jun 2024 22:49:10 GMT

Summary:

Lifted Metropolis X2X4.py

Werner — Mon, 10 Jun 2024 22:43:48 GMT

Summary: /* References */

Factor Metropolis X2X4.py

Werner — Mon, 10 Jun 2024 22:42:04 GMT

Summary: /* References */

Metropolis X2X4.py

Werner — Mon, 10 Jun 2024 22:25:02 GMT

Summary: /* References */

LiftedTASEPCompact.py

Werner — Mon, 10 Jun 2024 14:39:55 GMT

Summary: /* References */

This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

My Lecture 3 is concerned with the Symmetric Simple Exclusion Process (SSEP), and its liftings, the TASEP (totally asymmetric simple exclusion process) and the lifted TASEP, treated here. All these dynamical systems carry the word "Process" in their descriptions. This is because they are usually described in continuous time. We rather use a formulation in discrete time, where at each time step, a single move is attempted. In fact, each move consists in the choice of a random particle and the choice of a random direction.

Here, we are thus concerned with the lifted TASEP. With periodic boundary conditions, we may separate the forward-and-backward moving TASEP, as discussed in Lecture 3, into two independent copies. At each time step, the TASEP samples the random particle to be moved (forward). This, as discussed, is one half of a lifting of the SSEP.

==Python program==

import math
import random
alpha = 0.8
exponent = 2.0
prefactor = 1.0
NPart = 100; NSites = 2 * NPart
NIter = int(prefactor * NPart ** exponent)
NStrob = NIter // 40
Conf = [1] * NPart + [0] * (NSites - NPart)
Active = random.randint (0, NSites - 1)
while Conf[Active] != 1: Active = random.randint(0, NSites - 1)
Text = 'Periodic lifted TASEP, N= ' + str(NPart) + ', L= ' + str(NSites) + ', alpha= ' + str(alpha)
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))
print('-' * (NSites + 2))
for iter in range(NIter):
NewActive = (Active + 1) % NSites
if Conf[NewActive] == 0:
Conf[Active], Conf[NewActive] = 0, 1
Active = NewActive
if random.uniform(0.0, 1.0) < alpha:
while True:
Active = (Active - 1) % NSites
if Conf[Active] == 1: break
if iter % NStrob == 0:
PP = '|'
for k in range(NSites):
if Conf[k] == 0: PP += ' '
elif Active == k: PP += '^'
else: PP += 'X'
print(PP + '|')
print('-' * (NSites + 2))
Text = 'Total time = ' + str(prefactor) + ' * N ^ ' + str(exponent)
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))

==Output==

Here is output of the above Python program for the lifted TASEP with, for simplicity, N=32, L=64 and only 20 configurations over the length of the simulation. The caret ^ indicates the active particle in the sample space that is lifted with respect to the SSEP.

===For alpha = 0.8 ===

Periodic lifted TASEP, N= 32, L= 64, alpha= 0.8
------------------------------------------------------------------
|XXXXXXXXXXXXXXXXXXXXXXXXXXX^XXXX |
|XXXXXXXXXXXXXXXXXXXXXXXXXX X^XXXX |
|XXXXXXXXXXXXXXXXXXXXX^ XXXXXXX XXX |
|XXXXXXXXXXXXXXXXXX X XXXXXXXX^ XXXX |
|XXXXXXXXXXXXXXXXXX X XXXXXXX^XXXXX |
|XXXXXXXX^ XXXXX XXXX X XXXXXXXXXXXX X |
|XXXXX XXXXXXX XX^XXX X XXXXXXXXXXXX X |
|XXXXX ^XXXXXX XX XXXXXXXXXXXXXXXXX X |
|XXXXX X X^XXXXXX XXXXXXXXXXXXXXXXX X |
|XXXXX X XXXXXX XXXX^XXXXXXXXXXXXXX X |
|XXXXX X XXXXXX XXXXXXXXXXXXXXXX^XX X |
|XXXXX X XXXXXX XXXX XXXXXX^XXXXXXXXX |
|XXXXX X XXXXXX XXXX XXXXXXXXXXX^ X XXX |
|XXXXX X XXXXXX XXXX X ^XX XXXXXXXXXXXX |
|XXXXX X XXXXXX XXXX XXXXXXXX^XXXXXXX |
|XXXXX X XXXXXX XXXX XXXXXXXXX XX^XXXX |
|XXXXX X XXXXXX XXXX XXXXXXXXX XXX XXX^ |
|XXXXX X XXXXXX XXXX XXXXXXX ^ XXXX X XXX |
|XXXXX X XXXXXX XXXX XXXXXX^ X XXXX XXXX |
|XXXXX X X XXXXXX^XX XXXXXXX X XXXX XXXX |
|XXXXX X X ^XXXXX XXX XXXXXXX X XXXX XXXX |
------------------------------------------------------------------
Total time = 1.0 * N ^ 2.0

Clearly, N^2 steps are not sufficient to relax to equilibrium (as the right part of the box remains empty.

===For alpha = 0.5 ===

Periodic lifted TASEP, N= 32, L= 64, alpha= 0.5
------------------------------------------------------------------
|XXXXXX^XXXXXXXXXXXXXXXXXXXXXXXXX |
|XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX^ |
|XXXXXXXXXXXXXXXXXXXXXX XXXXX X X^ XX |
|XXXXXXXXXXXXXXXXXXXXXX X^ XXX X X XX X |
|XXXXXXXXXXXXXXXXXXXXXX X X X X ^ XXXX X |
|XXXXXXXXXXXXXXXXXX ^ X XX XX XX XXXXX X |
|XXXXXXXXXXXXXXXXX XX XXXXX^ XXXXXX X |
|XXXXXXXXXXXXXXXXX XX XXXXX XXXX ^ X XX |
|XXXXXXXXXXXXXXXXX XX XXXXX XX X XX X^ X |
|XXXXXXXXXXXXXXXXX XX XXXXX XX ^ XX X X X |
|XXXXXXXXXXXX X^X X X X X X XXXX X X XXX X X X |
|XXXXXXXXXX X XXX XX X ^X X XXXX X X XXX X X X |
|XXXXXXXXXX X XXX XX X XX XX XX^ X X XXX X X X |
|XXXXXXXXXX X XXX XX X XX X^X XX XXXXX X X X |
|XXXXXXXXXX X XXX XX X XX XX XX XX ^ XX X XX X |
|XXXXXXXXXX X X XX ^ XX X X XX X X XXX XX X XX X |
|XXXXXXXX X X ^X X XXXXX X XX X X XXX XX X XX X |
| XXXX XX XX XXXX X XXXXX X XX X X XXX XX X XX ^ |
| XXXX XX XX XXXX X XXXXX X XX X X X^ X XXXXX X |
| XXXX XX XX XX X XX^X X XXX X X X X XX X XXXXX X |
| XXXX XX XX XX X XX X X XXX X X ^ XX XX X XXXXX X |
------------------------------------------------------------------
Total time = 1.0 * N ^ 2.0

At the critical value alpha = 0.5 (in general alpha = N / L), N^2 steps appear to be sufficient to relax the system to equilibrium, and we conjecture the mixing time of the lifted TASEP to equal N^2.

==Further Information==

Further analysis of the lifted TASEP is performed in Massoulié et al (2025). In that paper, we consider an ensemble average of the above picture, plotting the expected density as a function of position and of time obtained from thousands of simulations. Also

==References==
* Essler F. H. L, Krauth W., Lifted TASEP: a Bethe ansatz integrable paradigm for non-reversible Markov chains, [https://doi.org/10.1103/PhysRevX.14.041035 Phys. Rev. X 14, 041035 (2024)].
* Kapfer S. C. and Krauth W., Irreversible Local Markov Chains with Rapid Convergence towards Equilibrium, [https://doi.org/10.1103/PhysRevLett.119.240603 Phys. Rev. Lett. 119, 240603 (2017)].
* Massoulié B., Erignoux C., Toninelli C. and Krauth W., Velocity trapping in the lifted TASEP and the true self-avoiding random walk [https://doi.org/10.1103/mqdr-x95j Phys. Rev. Lett. 135, 127102 (2025)]

TASEPCompact.py

Werner — Mon, 10 Jun 2024 14:03:51 GMT

Summary:

This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

My Lecture 3 is concerned with the [[SSEPCompact.py|Symmetric Simple Exclusion Process (SSEP)]] and its liftings, the TASEP (totally asymmetric simple exclusion process) (''treated here'') and the [[LiftedTASEPCompact.py|lifted TASEP]]. All these dynamical systems carry the word "Process" in their descriptions. This is because they are usually described in continuous time. We rather use a formulation in discrete time, where at each time step, a single move is attempted. In fact, each move consists in the choice of a random particle and the choice of a random direction.

Here, as mentioned, we are concerned with the TASEP. With periodic boundary conditions, we may separate the forward-and-backward moving TASEP, as discussed in Lecture 3, into two independent copies. At each time step, the TASEP samples the random particle to be moved (forward). This, as discussed, is one half of a lifting of the SSEP.

==Python program==

import math
import random
exponent = 2.0
alpha = 0.5
prefactor = 1.0
NPart = 32; NSites = 2 * NPart
NIter = int(prefactor * NPart ** exponent)
NStrob = NIter // 20
Conf = [1] * NPart + [0] * (NSites - NPart)
Text = 'Periodic TASEP, N= ' + str(NPart) + ', L= ' + str(NSites)
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))
print('-' * (NSites + 2))
for iter in range(NIter):
Active = random.randint (0, NSites - 1)
while Conf[Active] != 1:
Active = random.randint(0, NSites - 1)
Step = 1
NewActive = (Active + Step) % NSites
if Conf[NewActive] == 0:
Conf[Active], Conf[NewActive] = 0, 1
if iter % NStrob == 0:
PP = str()
for k in range(NSites):
if Conf[k] == 0: PP += ' '
else: PP += 'X'
print('|' + PP + '|')
print('-' * (NSites + 2))
Text = 'Total time = ' + str(prefactor) + ' * N ^ ' + str(exponent)
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))

==Output==
.
Here is output of the above Python program for the TASEP with, for simplicity, N=32, L=64. Over the entire life-time of the simulations (here N^2 single moves), only 20 configurations are shown, one roughly every 50 steps.

Periodic TASEP, N= 32, L= 64
------------------------------------------------------------------
|XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX |
|XXXXXXXXXXXXXXXXXXXXXXXXXXXX XXX X |
|XXXXXXXXXXXXXXXXXXXXXXXXX XXXX XX X |
|XXXXXXXXXXXXXXXXXXXXXXXX XXXXX X X X |
|XXXXXXXXXXXXXXXXXXXXXX XXXXXX X XX X |
|XXXXXXXXXXXXXXXXXXXXX XXXXX X XXX XX |
|XXXXXXXXXXXXXXXXXXXXX XXX XX XXXX X X |
|XXXXXXXXXXXXXXXXXXX XXXXX X XX XX XXX |
|XXXXXXXXXXXXXXX XXXXXXXXX X XX X XXXX |
|XXXXXXXXXXXXXX XXXXXXXX XX XX XX XXX X |
|XXXXXXXXXXXXXX XXXXXXXX XXXX XX X X XX |
|XXXXXXXX XXXXXXXXXXXXX X XXXX XX XX XX |
|XXXX XXXXXXXXXXXXXXXXX XX XX X X XX X X X |
|XX XXXXXXXXXXXXXXXXXX X X XXX X XX X XX X |
| XXXXXXXXXXXXXXXXXX X XX XXX X XX XXX X X |
| XXXXXXXXXXXXXXX X XX XX X X X X XX X XX X X X |
| XXXXXXXXXXXXXX XX X XXX X XX X XX XXX X X X |
| XXXXXXXXXXX XX XXX XX XX XX X X XXXXX X X X |
| XXXXXXXXX XXX XXXX XX XXXX X XXX X X XX X X |
| XXXXXXXX XXXX XXXX XX XXX X XX XX XX X X X X |
| XXXXXXX XXXX XX XXX XX XX X X XX X XXX X X X X |
------------------------------------------------------------------
Total time = 1.0 * N ^ 2.0

With all its limitations, the program illustrates that the N^2 steps are not sufficient to move the compact initial state and to approach the equilibrium.

==Further Information==

* The mixing behavior of the TASEP t_mix \sim N^(5/2) has been computed by Baik & Liu (2016) (see references). This is in our units, where one move takes place per unit time.
* The relaxation time of the TASEP is likewise t_rel \sim N^(5/2). Again this is in the MCMC units of one move per time step.
* The description of the TASEP as a lifting of the SSEP is from Kapfer et al. (2017) (see references).
* An example transition matrix for the TASEP with open boundary conditions (hard walls) can be found in Essler & Krauth (2024).
* The TASEP can be solved by Bethe ansatz, and many of its properties are known.
* As discussed at length in Lecture 3, the transition matrix of the TASEP is doubly stochastic so that, in equilibrium, every allowed configuration is equally likely.

==References==
* Baik, J., and Z. Liu, Z. TASEP on a Ring in Sub-relaxation Time Scale. J. Stat. Phys., 165(6): 1051–1085, 2016.
* Dhar D. An exactly solved model for interfacial growth. Phase Transitions, 9(1):51, 1987
* Essler F. H. L., Krauth W., Lifted TASEP: a Bethe ansatz integrable paradigm for non-reversible Markov chains, [https://doi.org/10.1103/PhysRevX.14.041035 Phys. Rev. X 14, 041035 (2024)].
* Kapfer S. C. and Krauth W., Irreversible Local Markov Chains with Rapid Convergence towards Equilibrium, [https://doi.org/10.1103/PhysRevLett.119.240603 Phys. Rev. Lett. 119, 240603 (2017)].

SSEPCompact.py

Werner — Mon, 10 Jun 2024 12:31:12 GMT

Summary: /* Python program */

==Context==
This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

==Introduction==

My Lecture 3 is concerned with the ''Symmetric Simple Exclusion Process'' (SSEP), treated here, and its liftings, the TASEP (totally asymmetric simple exclusion process) and the ''lifted'' TASEP. All these dynamical systems carry the word "''Process''" in their descriptions because they are usually described in continuous time. Here, we rather use a formulation in discrete time, where at each time step t=0,1,2,..., a single move is attempted. In fact, each move consists in the choice of a random particle and the choice of a random direction. The SSEP is a local diffusive Markov chain, and it has very slow dynamics: it takes N^3 log N steps to get it from a compact initial state into equilibrium.

==Python program==
import math
import random
exponent = 3.0
alpha = 0.5
prefactor = 1.0
NPart = 100; NSites = 2 * NPart
NIter = int(prefactor * NPart ** exponent * math.log(NPart))
NStrob = NIter // 400
Conf = [1] * NPart + [0] * (NSites - NPart)
Active = random.randint (0, NSites - 1)
while Conf[Active] != 1:
Active = random.randint(0, NSites - 1)
Text = 'Periodic SSEP, N= ' + str(NPart) + ', L= ' + str(NSites)
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))
print('-' * (NSites + 2))
for iter in range(NIter):
Active = random.randint (0, NSites - 1)
while Conf[Active] != 1:
Active = random.randint(0, NSites - 1)
Step = random.choice([-1, 1])
NewActive = (Active + Step) % NSites
if Conf[NewActive] == 0:
Conf[Active], Conf[NewActive] = 0, 1
if iter % NStrob == 0:
PP = str()
for k in range(NSites):
if Conf[k] == 0: PP += ' '
else: PP += 'X'
print('|' + PP + '|')
print('-' * (NSites + 2))
Text = 'Total time = ' + str(prefactor) + ' * N ^ ' + str(exponent) + ' * log N'
print(' ' * (NSites// 2 + 1 - len(Text) // 2) + Text + ' ' * (NSites// 2 + 1 - len(Text) // 2))

This example program performs a large number of iterations of the Monte Carlo algorithm for the Symmetric Simple Exclusion Process, and plots 40 lines of output over the entire simulation time.

==Output==

Here is output of the above Python program with, for simplicity, N=32, L=64 and only 20 configurations over the length of the simulation.

The initial configuration is compact. Clearly, the simulation has not run long enough to forget its initial state, and it would be necessary to increase the simulation time from N^2 log N to N^3 log N. In our simplified setting, the logarithm is difficult to see, but better analysis tools readily extract it from the numerical data.

Periodic SSEP, N= 32, L= 64
------------------------------------------------------------------
|XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX |
|XXXXXXXXXXXXXXXXXXXXXXXXXXXXX XX X |
| XXXXXXXXXXXXXXXXXXXXXXXXXXXXX X X X|
| X XXXXXXXXXXXXXXXXXXXXXXXXXXX X X X X|
| X XXXXXXXXXXXXXXXXXXXXXX X XXX X X XX X|
|X XXXXXXXXXXXXXXXXXXXXX XXX X XX X X X X |
|X XXX XXXXXXXXXXXXXXXXXXXX XX X X X X X X |
| XXXX XXXXXXXXXXXXXXXXXXXX X X XX X X X X|
| XXXXXX XXXXXXXXXXXXXXXXXX X X X X X X X X |
| XXXXXXXX XXXXXXXXXXXXXXX X X X XX X XX X|
| XXXXXXX XXXXXXXXXXXXX XXXX X X XX XX X X |
| XXXX XXXXXXXXXXXXXX XXXXX X XXXX X X X X |
|X XX XXXXXXXXXXXX XXXXXXX X XX XX X X X XX |
|X XXXX XXXXXXXXXX XXXX XXXX XXX X X X X XX|
| XXX XXXXXXXXXXX XXXXXXXXX X X X X X X X X X|
|X XXX XXXXXXXXXXX XXXX XXXXX X X X XX X XX|
| XXXXX XXXXXXXXX X XXXXXXXXX X X XXX X XX|
|XXXX XXXXXXXXXXXX X XXXXX XX X XX X X XX X|
|XXXX XXXXXXXXX XXXXXXXX XX X XX X X XX X X|
|XXXXXX XXXXXX XXXXXXXXXXX X X X X XX X X X |
|X XXXXXX XX XXXXXXXXXXXXX XX X X X X X X X X|
------------------------------------------------------------------
Total time = 1.0 * N ^ 2.0 * log N

==Further Information==

* The mixing behavior of the SSEP t_mix \sim N^3 log N has been computed by Lacoin (2016, 2017) (see references).
* The relaxation time of the SSEP is t_rel \sim N^3 (without the logarithm). It is thus asymptotically smaller than the mixing time, leading to the cutoff phenomenon.

==References==
* Essler F. H. L, Krauth W., Lifted TASEP: a Bethe ansatz integrable paradigm for non-reversible Markov chains, [https://doi.org/10.1103/PhysRevX.14.041035 Phys. Rev. X 14, 041035 (2024)]
* Lacoin H., The cutoff profile for the simple exclusion process on the circle, Ann. Probab. 44, 3399 (2016)
* Lacoin H., The simple exclusion process on the circle has a diffusive cutoff window, Ann. Inst. H. Poincaré Probab.Statist. 53, 1402 (2017).
* Kapfer S. C. and Krauth W., Irreversible Local Markov Chains with Rapid Convergence towards Equilibrium, [https://doi.org/10.1103/PhysRevLett.119.240603 Phys. Rev. Lett. 119, 240603 (2017)].

Diffusion CFTP coupl pos.py

Werner — Thu, 06 Jun 2024 13:32:45 GMT

Summary:

Diffusion forward.py

Werner — Thu, 06 Jun 2024 13:29:20 GMT

Summary: /* Context */

Top to random TVD.py

Werner — Thu, 06 Jun 2024 13:27:39 GMT

Summary:

Top to random simul stop.py

Werner — Thu, 06 Jun 2024 13:24:03 GMT

Summary:

Sample transformation power.py

Werner — Thu, 06 Jun 2024 13:22:15 GMT

Summary:

Sample transformation exp.py

Werner — Thu, 06 Jun 2024 13:20:57 GMT

Summary:

Top to random eigenvalues.py

Werner — Thu, 06 Jun 2024 13:17:47 GMT

Summary:

Bernoulli poisson.py

Werner — Thu, 06 Jun 2024 12:51:35 GMT

Summary:

Direct surface 2d.py

Werner — Thu, 06 Jun 2024 12:48:29 GMT

Summary:

Naive surface 2d.py

Werner — Thu, 06 Jun 2024 12:47:35 GMT

Summary:

Buggy surface 2d.py

Werner — Thu, 06 Jun 2024 12:46:04 GMT

Summary:

Bernoulli two pebbles patch.py

Werner — Thu, 06 Jun 2024 05:33:42 GMT

Summary:

Bernoulli two pebbles.py

Werner — Thu, 06 Jun 2024 05:31:21 GMT

Summary:

Bernoulli.py

Werner — Thu, 06 Jun 2024 05:29:43 GMT

Summary:

Diffusion.py

Werner — Thu, 06 Jun 2024 05:26:51 GMT

Summary:

BegRohu Lectures 2024

Werner — Thu, 06 Jun 2024 05:13:44 GMT

Summary: /* Context */

==Context==
This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

A "few" years ago, as if it was yesterday, I taught at the 1996 Beg Rohu Summer School on an [https://arxiv.org/pdf/cond-mat/9612186 Introduction To Monte Carlo Algorithms]. The lecture notes of this first course were ''found readable and the drawings enjoyable'', so I made them into a 2006 book for Oxford University Press entitled "Statistical Mechanics: Algorithms and Computations".

The 2024 lectures discuss the foundations of the enormous corpus of works that have constituted the second revolution in Markov chains and Monte Carlo algorithm, namely the understanding of time scales in stochastic dynamics and the concepts of coupling, thinning and, last not least, the systematic construction of non-reversible (that is, non-equilibrium) Markov chains that nevertheless converge to the imposed equilibrium steady state. I hope that what I present here is again ''found readable and enjoyable'', ...

==List of Python program==

===Lecture 1===

[[Bernoulli.py | Bernoulli.py]]

[[Bernoulli_two_pebbles.py | Bernoulli_two_pebbles.py]]

[[Bernoulli_two_pebbles_patch.py | Bernoulli_two_pebbles_patch.py]]

[[Bernoulli_poisson.py | Bernoulli_poisson.py]]

[[Sample_transformation_exp.py | Sample_transformation_exp.py]]

[[Sample_transformation_power.py | Sample_transformation_power.py]]

[[Buggy_surface_2d.py | Buggy_surface_2d.py]]

[[Naive_surface_2d.py | Naive_surface_2d.py]]

[[Direct_surface_2d.py | Direct_surface_2d.py]]

[[Walker.py | Walker.py]]

===Lecture 2===

[[Image:One d cftp.png|left|600px|border|Coupling-from-the-past approach to sampling.]]
 

[[ Top_to_random_eigenvalues.py | Top_to_random_eigenvalues.py]] x β

[[ Top to random simul.py | Top to random simul.py]]

[[ Top_to_random_simul_stop.py | Top_to_random_simul_stop.py]]

[[Top_to_random_TVD.py |Top-to-random-TVD.py]]

[[ Top to random TVD mix rel.py | Top to random TVD mix rel.py]]

[[Diffusion.py | Diffusion.py]] x

[[Diffusion_forward.py | Diffusion_forward.py]]

[[Diffusion_CFTP_coupl_pos.py | Diffusion_CFTP_coupl_pos.py]]

[[Diffusion CFTP.py | Diffusion CFTP.py]] x β

===Lecture 3===

[[SSEPCompact.py | SSEPCompact.py]] x

[[SSEPSteady.py | SSEPSteady.py]]

[[TASEPCompact.py | TASEPCompact.py]] x

[[TASEPSteady.py | TASEPSteady.py]]

[[LiftedTASEPCompact.py | LiftedTASEPCompact.py]] x

[[LiftedTASEPSteady.py | LiftedTASEPSteady.py]]

===Lecture 4===

[[Metropolis_X2X4.py | Metropolis_X2X4.py ]]

[[Factor_Metropolis_X2X4.py | Factor_Metropolis_X2X4.py ]]

[[Factor_Metropolis_X2X4_patch.py | Factor_Metropolis_X2X4_patch.py ]]

[[Lifted_Metropolis_X2X4.py | Lifted_Metropolis_X2X4.py ]]

[[ZigZag_X2X4.py | ZigZag_X2X4.py ]]

[[Factor_ZigZag_X2X4.py | Factor_ZigZag_X2X4.py ]]

[[Bounded_Lifted_Metropolis_X2X4.py | Bounded_Lifted_Metropolis_X2X4.py ]]

[[Bounded_ZigZag_X2X4.py | Bounded_ZigZag_X2X4.py ]]

[[Bounded_Factor_ZigZag_X2X4.py | Bounded_Factor_ZigZag_X2X4.py ]]

===Lecture 5===

[[Stopping_circle.py | Stopping_circle.py ]]

[[SSEP_coupling.py | SSEP_coupling.py ]]

[[SSEP_coupling_FTP.py | SSEP_coupling_FTP.py ]]

[[Ising_coupling.py | Ising_coupling.py ]]

[[Ising_coupling_FTP.py | Ising_coupling_FTP.py ]]

[[Hard_spheres_coupling.py |Hard_spheres_coupling.py ]]

"x" means that the program "already" has some explanatory text, even though this text may not yet be complete

"β" means that it has been β tested by students at the 2024 Beg Rohu summer school

==References==

* Levin, D. A., Peres, Y. & Wilmer, E. L. Markov Chains and Mixing Times (American Mathematical Society, 2008)
* Krauth, W., Statistical Mechanics: Algorithms and Computations (Oxford University Press, 2006)
* Wasserman, L., All of Statistics (Springer Verlag, 2004)

See the individual pages of my lectures for a large number of more specific references.

Top to random TVD mix rel.py

Werner — Thu, 06 Jun 2024 00:08:16 GMT

Summary:

Top to random simul.py

Werner — Thu, 06 Jun 2024 00:05:53 GMT

Summary:

Walker.py

Werner — Thu, 06 Jun 2024 00:04:24 GMT

Summary: /* Python program */

Diffusion CFTP.py

Werner — Thu, 06 Jun 2024 00:02:29 GMT

Summary: /* Python program */

==Context==
This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

The present program illustrates the coupling-from-the-past approach to Markov-chain sampling introduced by Propp and Wilson in 1997. This paper is definitely part of the ''second Revolution''. It illustrates how Markov-chain sampling can obtain samples from the stationary distribution, without any corrections whatsoever. As in the other program in the series started by [[Diffusion.py |Diffusion.py]], we consider a particle diffusing on a path graph with five sites, with "down", "up", and "straight¨ moves that are cropped at the boundaries (there are no periodic boundary conditions). We consider the formulation of a Markov chain in terms of random maps, but run from time t=-infinity up to time t=0.

[[Image:One d cftp.png|left|600px|border|Coupling-from-the-past approach to sampling.]]
 
At time t=-infinity, as shown, the pebble starts on configuration 1 so that, by time t=0, it has run for an infinite time. Whatever the mixing time and the relaxation time, by t=0, it must sample the stationary distribution (the uniform distribution on the five sites). The below program constructs just as many sets of arrows as needed to conclude where it finishes its infinitely long journey.

==Python program==

import random
import matplotlib.pyplot as plt

N = 5
pos = []
for stat in range(100000):
all_arrows = {}
time_tot = 0
while True:
time_tot -= 1
arrows = [random.choice([-1, 0, 1]) for i in range(N)]
if arrows[0] == -1: arrows[0] = 0
if arrows[N - 1] == 1: arrows[N - 1] = 0
all_arrows[time_tot] = arrows
positions = set(range(0, N))
for t in range(time_tot, 0):
positions = set([b + all_arrows[t][b] for b in positions])
if len(positions) == 1: break
a = positions.pop()
pos.append(a)
plt.title('Backward coupling: 1-d with walls: position at t=0')
plt.hist(pos, bins=N, range=(-0.5, N - 0.5), density=True)
plt.savefig('backward_position_t0.png')
plt.show()

Note that "arrows" is a list, and all_arrows a dictionary containing the arrows for each position (0,1,...,N - 1). This dictionary contains exactly the random map of arrows in the above picture, from time time_tot to time -1.

==Output==

Here is output of the above Python program. The histogram is absolutely flat, without any corrections. But this is normal, given that the simulation has run, for each of the realizations of the random map, an infinite number of iterations.

[[Image:Backward position t0.png|left|600px|border|Coupling-from-the-past approach to sampling.]]
 

==Further Information==

* Rather than to run the simulation until t=0, one might be tempted to stop it as soon as the configurations have coupled at some time t < 0. This is tested in the following program, but it is wrong for obvious reasons: The goal was to infer the position at t=0 of a simulation that has run since time t=-infinity.
* The above program is easily adapted to a non-uniform stationary distribution pi(1),...,pi(N-1), and its output is then equivalent to what we would obtain using the direct-sampling algorithms that I discussed in BegRohu Lecture 1, first and foremost [[Walker.py| Walker's algorithm]]. However, Walker's algorithm always requires the stationary probabilities on the sample space to be evaluated and rearranged in some way, whereas the CFTP approach to Markov chains can be often applied for huge sample spaces, way larger than what can be listed. This is itself an astonishing story, that I hope to have time to relate in later lectures.
* A truly astonishing direct-sampling algorithm [[Stopping circle.py | Stopping_cycle.py]], based on the concept of strong stationary stopping times, is discussed in Lecture 5. However, this algorithm comes with a no-go theorem by Lovász and Winkler, that it is only correct for the cycle graph (the graph in the above example, but with periodic boundary conditions), and for the complete graph.

==References==

The reference is:

Heatbath ising.py

Werner — Mon, 04 Mar 2024 22:34:29 GMT

Summary:

This page presents the Python3 program heatbath_ising.py, a Markov-chain algorithm for the Ising model on an L x L square lattice in two dimensions with periodic boundary conditions. This program uses the heatbath algorithm.

__FORCETOC__
=Description=
The program is described in my 2024 Oxford Lecture No 8, and also in my book. This version of the program estimates the energy per particle, and the specific heat.
=Program=

import random, math

L = 6
N = L * L
nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N,
(i // L) * L + (i - 1) % L, (i - L) % N) \
for i in range(N)}
nsteps = 10000000
beta = 1.0
S = [random.choice([-1, 1]) for site in range(N)]
E = -0.5 * sum(S[k] * sum(S[nn] for nn in nbr[k]) \
for k in range(N))
E_tot, E2_tot = 0.0, 0.0
for step in range(nsteps):
k = random.randint(0, N - 1)
Upsilon = random.uniform(0.0, 1.0)
h = sum(S[nn] for nn in nbr[k])
Sk_old = S[k]
S[k] = -1
if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * h)):
S[k] = 1
if S[k] != Sk_old:
E -= 2.0 * h * S[k]
E_tot += E
E2_tot += E ** 2
E_av = E_tot / float(nsteps)
E2_av = E2_tot / float(nsteps)
c_V = beta ** 2 * (E2_av - E_av ** 2) / float(N)
print(E_av / N,c_V)

=Version=
See history for version information.

[[Category:Python]] [[Category:Oxford_2024]] [[Category:MOOC_SMAC]]

Two cycles.py

Werner — Tue, 27 Feb 2024 02:32:50 GMT

Summary:

This is a Python3 program to sample random permutations P of N elements subject to the constraint that the lengths of the cycles in P can only be one or two. The algorithm is based on a recursive-sampling strategy discussed in my book and presented in [[Oxford_Lectures_2025#Lecture_7:_04_March_2025|Lecture 7]] of my [[Oxford_Lectures_2025|2025 Oxford lectures]].

__FORCETOC__
=Description=

To sample permutations of N elements which only contain cycles of length 1 or 2, we derived in Lecture 7 a recursion formula for the number Y_N of such permutations (eq. (7.6)), namely
Y_N = Y_{N-1} + (N - 1) Y_{N-2}
The first term on the right-hand site contains the number of permutations such that the last element is itself in a cycle of length 1, while the second term corresponds to permutations of length 2. It thus suffices to first compute the sequence of the Y_0, Y_1, Y_2, .... Y_N, then sample a random number between 0 and Y_N, and then check whether it is smaller than Y_{N-1}. If it is, the last element is in a cycle of length 1, else in one of length 2. We can continue our sampling either with N - 1 particles left in a permutation that we know nothing about, or else with N - 2 particles. All this may sound a bit mysterious, but it certainly works out OK!

The reason we discuss this algorithm is that a very similar algorithm is used for sampling cycles in the Bose-gas problem of ideal particles (for example, in a harmonic trap in three dimensions).

=Program=

import random
from sympy.combinatorics import Permutation
Y = {-1: 0, 0:1, 1:1}
N = 4
Stats = {}
for k in range(2, N + 1):
Y[k] = Y[k - 1] + (k - 1) * Y[k - 2]
for iter in range(100000):
Q = list(range(N))
random.shuffle(Q)
M = N
P = []
while M > 0:
if random.uniform(0.0, Y[M]) < Y[M - 1]:
P.append([Q[M - 1]])
M -= 1
else:
P.append([Q[M - 1] , Q[M - 2]])
M -= 2
P = tuple(Permutation(P).array_form)
Stats[P] = Stats.get(P, 0) + 1
print(Stats)

=Output=

The output is a dictionary of permutations (in array form) together with the number of times they were sampled. Here we see, that all the permutations are (roughly) equally frequent, that is, were drawn with equal probabilities. We may check that all these permutations have cycles of at most two.

{(0, 3, 2, 1): 9801, (0, 2, 1, 3): 10281, (0, 1, 2, 3): 10145, (2, 3, 0, 1): 9940, (3, 1, 2, 0): 9869,
(0, 1, 3, 2): 9950, (2, 1, 0, 3): 9830, (3, 2, 1, 0): 10154, (1, 0, 3, 2): 10016, (1, 0, 2, 3): 10014}

=Version=
See history for version information.
For more information on this program, see [[Oxford_Lectures_2025#Lecture_7:_04_March_2025|Lecture 7]] of my [[Oxford_Lectures_2025|2025 Oxford lectures]].
[[Category:Python]] [[Category:Oxford2025]]

Oxford Lectures 2024

Werner — Mon, 26 Feb 2024 16:59:02 GMT

Summary:

My [https://www.physics.ox.ac.uk/events/algorithms-and-computations-theoretical-physics-set-lectures 2024 Public Lectures at the University of Oxford (UK)] run from 16 January 2024 through 5 March 2024.

[http://www.lps.ens.fr/%7Ekrauth/images/c/c5/WK_Lecture1_Oxford2024.pdf Here are the notes for the first lecture. ]

[http://www.lps.ens.fr/%7Ekrauth/images/5/5a/WK_Lecture2_Oxford2024.pdf Here are the notes for the second and third lectures. ]

[http://www.lps.ens.fr/%7Ekrauth/images/1/19/WK_Lecture4_Oxford2024.pdf Here are the preliminary notes for the fourth lecture.]

[http://www.lps.ens.fr/%7Ekrauth/images/6/63/WK_Lecture5_Oxford2024.pdf Here are the preliminary notes for the fifth lecture.]

[http://www.lps.ens.fr/%7Ekrauth/images/9/98/WK_Lecture6_Oxford2024.pdf Here are the notes for the sixth lecture.]

[http://www.lps.ens.fr/%7Ekrauth/images/b/b4/WK_Lecture7_Oxford2024.pdf Here are the notes for the seventh lecture.]

[[Two_cycles.py| Two_cycles.py]]

[[direct_N_bosons.py| Direct_bosons.py]]

[http://www.lps.ens.fr/%7Ekrauth/images/f/fe/WK_Lecture8_Oxford2024.pdf Here are the notes for the eighth lecture.]

[[Enumerate_ising.py|Enumerate_ising.py]]

[[Markov_ising.py|Markov_ising.py]]

[[Heatbath_ising.py|Heatbath_ising.py]]

[[Cluster_ising.py| Cluster_ising.py]]

Direct N bosons.py

Werner — Mon, 26 Feb 2024 16:49:52 GMT

Summary:

This is a Python3 program for the direct sampling of N non-interacting (that is, ideal) bosons in a three-dimensional harmonic trap. The program is discussed in my book and presented in [[Oxford_Lectures_2025#Lecture_7:_04_March_2025|Lecture 7]] of my 2025 Oxford lectures.

__FORCETOC__

=Description=

=Program=

import pylab, math, random
def z(beta, k):
sum = 1.0 / (1.0 - math.exp(-k * beta)) ** 3
return sum
def canonic_recursion(beta,N):
Z = [1.0]
for M in range(1, N + 1):
Z.append(sum(Z[k] * z(beta, M - k) for k in range(M)) / M)
return Z
def pi_list_make(Z, M):
pi_list = [0] + [z(beta, k) * Z[M - k] / Z[M] / M for k in range(1, M + 1)]
pi_sum = [0]
for k in range(1, M + 1):
pi_sum.append(pi_sum[k - 1] + pi_list[k])
return pi_sum
def tower_sample(data, Upsilon): #fully naive tower sampling
for k in range(len(data)):
if Upsilon < data[k]: break
return k
def levy_harmonic_path(Del_tau, N):
beta = N * Del_tau
xN = random.gauss(0.0, 1.0 / math.sqrt(2.0 * math.tanh(beta / 2.0)))
x = [xN]
for k in range(1, N):
Upsilon_1 = 1.0 / math.tanh(Del_tau) + 1.0 / math.tanh((N - k) * Del_tau)
Upsilon_2 = x[k - 1]/ math.sinh(Del_tau) + xN / math.sinh((N - k) * Del_tau)
x_mean = Upsilon_2 / Upsilon_1
sigma = 1.0 / math.sqrt(Upsilon_1)
x.append(random.gauss(x_mean, sigma))
return x

N = 100 # this naive program may overflow for N > 100
T_star = 0.1
T = T_star * N ** (1.0 / 3.0)
beta = 1.0 / T
Z = canonic_recursion(beta, N)
print(Z)
M = N
x_config = []
y_config = []
while M > 0:
pi_sum = pi_list_make(Z, M)
Upsilon = random.uniform(0.0, 1.0)
k = tower_sample(pi_sum, Upsilon)
x_config += levy_harmonic_path(beta, k)
y_config += levy_harmonic_path(beta, k)
M -= k
pylab.figure(1)
pylab.axis('scaled')
pylab.xlim(-5.0, 5.0)
pylab.ylim(-5.0, 5.0)
pylab.plot(x_config, y_config,'ro')
pylab.xlabel('$x$')
pylab.ylabel('$y$')
pylab.title('3d bosons in trap (2d projection):'+' $N$= '+str(N)+' T* =' + str(T_star))
pylab.show()

=Output=

[[Image:Boson trap fast.gif|100px|right|frame|Direct-sampling algorithm for ideal bosons in a trap. The output shown here was produced with a more elaborate version of the program than the one shown on this page]]

 

=Further reading=

The Path-integral Monte Carlo algorithm

=Version=
See history for version information.
For more information on this program, see Lecture 7 of my 2025 Oxford lectures.
[[Category:Python]] [[Category:Oxford2025]]

Coupling ising.py

Werner — Fri, 17 Feb 2023 17:17:32 GMT

Summary:

This page presents the program coupling_ising.py, a heat-bath algorithm for the Ising model on an LxL square lattice in two dimensions, run for two configurations at a time. The algorithm illustrates the coupling phenomenon.

__FORCETOC__
=Description=

=Program (in Python3)=
import random, math

L = 7
N = L * L
nbr = {i : ((i // L) * L + (i + 1) % L, (i + L) % N,
(i // L) * L + (i - 1) % L, (i - L) % N) \
for i in range(N)}
NIter = 100
for TT in range(20, 40):
T = TT / 10
beta = 1.0 / T
MeanCoupling = 0
for iter in range(NIter):
S1 = [-1] * N
S2 = [1] * N
step = 0
while True:
step += 1
k = random.randint(0, N - 1)
Upsilon = random.uniform(0.0, 1.0)
h1 = sum(S1[nn] for nn in nbr[k])
S1[k] = -1
if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * h1)): S1[k] = 1
h2 = sum(S2[nn] for nn in nbr[k])
S2[k] = -1
if Upsilon < 1.0 / (1.0 + math.exp(-2.0 * beta * h2)): S2[k] = 1
if S1 == S2:
MeanCoupling += step
break
print(T, MeanCoupling / NIter)
=Output=
A slightly modified graphics version of this program produces the following output:

[[Image:IsingCoupling.png|left|50px]]
 
=Version=
See history for version information.

[[Category:Python]] [[Category:Honnef_2015]] [[Category:MOOC_SMAC]]

Advanced topics MCMC 2023

Werner — Wed, 11 Jan 2023 09:27:33 GMT

Summary:

This is the homepage of my 2023 lecture course on Advanced topics in Markov-chain Monte Carlo

[http://www.lps.ens.fr/%7Ekrauth/images/f/f6/CM_2023_1.1.pdf Here is the CM1.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/2/23/CM_2023_1.2.pdf Here is the CM1.2]

[http://www.lps.ens.fr/%7Ekrauth/images/4/45/TD01_ICFP_ADV_MCMC.pdf Here is the TD1]

[http://www.lps.ens.fr/%7Ekrauth/images/4/4d/CM_2023_2.1.pdf Here is the CM2.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/1/11/CM_2023_2.2.pdf Here is the CM2.2]

[http://www.lps.ens.fr/%7Ekrauth/images/5/51/TD02_2023_ICFP_ADV_MCMC.pdf Here is the TD2]

[http://www.lps.ens.fr/%7Ekrauth/images/e/ef/CM_2023_3.1.pdf Here is the CM3.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/a/a4/CM_2023_3.2.pdf Here is the CM3.2]

[http://www.lps.ens.fr/%7Ekrauth/images/2/25/TD03_2023_ICFP_ADV_MCMC.pdf Here is the TD3]

[http://www.lps.ens.fr/%7Ekrauth/images/3/3e/CM_2023_4.1.pdf Here is the CM4.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/1/1f/CM_2023_4.2.pdf Here is the CM4.2]

[http://www.lps.ens.fr/%7Ekrauth/images/1/11/TD04_2023_ICFP_ADV_MCMC.pdf Here is the TD4]

[http://www.lps.ens.fr/%7Ekrauth/images/d/df/CM_2023_5.1.pdf Here is the CM5.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/8/87/CM_2023_5.2.pdf Here is the CM5.2]

[http://www.lps.ens.fr/%7Ekrauth/images/9/92/CM_2023_5.3.pdf Here is the CM5.3]

[http://www.lps.ens.fr/%7Ekrauth/images/2/28/TD05_2023_ICFP_ADV_MCMC.pdf Here is the TD5]

[http://www.lps.ens.fr/%7Ekrauth/images/8/82/ControleMock_2023.pdf Here is the Mock Exam]

[[Coupling ising.py| here is the program Coupling_ising.py]]

[http://www.lps.ens.fr/%7Ekrauth/images/e/e3/CM_2023_6.1.pdf Here is the CM6.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/2/2c/CM_2023_6.2.pdf Here is the CM6.2]

[http://www.lps.ens.fr/%7Ekrauth/images/e/e6/TD06_2023_ICFP_ADV_MCMC.pdf Here is the TD6]

[http://www.lps.ens.fr/%7Ekrauth/images/7/72/CM_2023_8.1.pdf Here is the CM8.1] of my 2023 lecture course Advanced topics in Markov-chain Monte Carlo.

[http://www.lps.ens.fr/%7Ekrauth/images/2/29/CM_2023_8.2.pdf Here is the CM8.2]

[http://www.lps.ens.fr/%7Ekrauth/images/a/ab/TD08_2023_ADV_MCMC.pdf Here is the TD8]