Metropolis-Hastings for Classical Distributions

Hello! Welcome to your next lesson on numerical methods for quantum systems.

In our last lesson, you mastered Exact Diagonalization (ED), a powerful "brute-force" method for finding the exact energy spectrum and eigenstates of small quantum many-body systems. We also confronted its primary weakness: the curse of dimensionality, where the computational cost grows exponentially, making it intractable for all but the smallest systems.

To tackle larger, more realistic problems—a key part of your goal to study many-body physics—we must turn to methods that don't require storing the entire Hamiltonian. This brings us to the world of Monte Carlo methods.

Today, we will learn the engine that drives many of these powerful techniques. Your learning outcome for this lesson is to:

• Implement the Metropolis-Hastings algorithm for sampling a classical probability distribution.

This might seem like a step away from quantum mechanics, but it's a crucial one. The Metropolis-Hastings algorithm is a general and powerful tool for statistical sampling in high-dimensional spaces. We'll master its mechanics in a classical context today, and in our next lesson, we will apply it directly to solve quantum problems using Variational Monte Carlo.

1. The Challenge of High-Dimensional Integrals

At the heart of many problems in physics, especially statistical mechanics and quantum mechanics, lies the need to compute high-dimensional integrals. For instance, the expectation value of an observable $A$ in statistical mechanics is given by an integral over the entire phase space of a system:

$⟨A⟩=\frac{1}{Z}\int A(\mathbf{q},\mathbf{p})e^{-E(\mathbf{q},\mathbf{p})/k_{B}T}d\mathbf{q}d\mathbf{p}$

For $N$ particles in 3D, this is an integral over $6N$ dimensions! Standard numerical integration methods like the trapezoidal or Simpson's rule, which rely on evaluating the function on a grid, fail catastrophically here. If you use just 10 grid points per dimension, you'd need $10^{6N}$ points—an impossible number. This is another manifestation of the "curse of dimensionality."

Monte Carlo methods solve this by replacing a deterministic evaluation on a grid with a statistical sampling of the space.

The simplest Monte Carlo approach is to sample points uniformly and average them. However, this is inefficient if the function you're integrating is peaked in a small region of the space. A much better strategy is importance sampling, where you preferentially draw samples from regions where the function's value is high.

To do this, you need to be able to draw random numbers from a distribution that mimics your function. But what if that distribution is complex and non-standard? This is where the Metropolis-Hastings algorithm comes in. It's a method to generate samples from any probability distribution $P(\mathbf{x})$, even if we only know $P(\mathbf{x})$ up to a normalization constant.

2. The Metropolis-Hastings Algorithm: A "Smart" Random Walk

The core idea of the Metropolis-Hastings algorithm is to construct a "smart" random walk in the space of possible states (or configurations) $\mathbf{x}$. This walk, known as a Markov chain, is designed so that the points it visits, after an initial "burn-in" period, are distributed according to the target probability distribution $P(\mathbf{x})$.

Let's get a visual intuition for how this works.

Now let's formalize the algorithm's steps.

The Algorithm:

1. Initialize: Start at an arbitrary state ${\mathbf{x}}_0$.
2. Iterate for a large number of steps, $t=0,1,2,\dots$:
   a. Propose: Generate a candidate state ${\mathbf{x}}^{\prime }$ from a proposal distribution $q({\mathbf{x}}^{\prime }|{\mathbf{x}}_t)$. A common choice is a Gaussian distribution centered at the current state ${\mathbf{x}}_t$.
   b. Calculate Acceptance Ratio: Compute the ratio:
   $\alpha =\frac{P({\mathbf{x}}^{\prime })}{P({\mathbf{x}}_t)}\times \frac{q({\mathbf{x}}_t|{\mathbf{x}}^{\prime })}{q({\mathbf{x}}^{\prime }|{\mathbf{x}}_t)}$
   The first term is the ratio of target probabilities. The second term, the Hastings ratio, corrects for any asymmetry in the proposal distribution. If the proposal distribution is symmetric (e.g., a Gaussian step), then $q({\mathbf{x}}_t|{\mathbf{x}}^{\prime })=q({\mathbf{x}}^{\prime }|{\mathbf{x}}_t)$, and the ratio simplifies to the original Metropolis algorithm ratio:
   $\alpha =\frac{P({\mathbf{x}}^{\prime })}{P({\mathbf{x}}_t)}$
   c. Accept or Reject:
   • Generate a uniform random number $u\in [0,1]$.
   • If $u<\alpha$, accept the new state: ${\mathbf{x}}_{t+1}={\mathbf{x}}^{\prime }$.
   • Otherwise, reject the new state and stay put: ${\mathbf{x}}_{t+1}={\mathbf{x}}_t$.

A crucial feature of this algorithm is that it only depends on the ratio of probabilities, $P({\mathbf{x}}^{\prime })/P({\mathbf{x}}_t)$. This means you don't need to know the normalization constant of your probability distribution! This is immensely powerful, as the normalization constant (like the partition function $Z$) is often the very intractable integral we wanted to avoid in the first place.

The underlying principle that ensures this works is the condition of detailed balance. In equilibrium, the probability flow from any state $\mu$ to any other state $\nu$ must equal the flow from $\nu$ to $\mu$. The acceptance rule is constructed precisely to satisfy this.

3. Implementation and Diagnostics

Now, let's translate the algorithm into code. You will implement a simple Metropolis-Hastings sampler in Python.

As you implement and run your sampler, you will quickly discover that its performance depends critically on a few practical details.

Diagnostics and Tuning

How do you know if your sampler is working correctly and efficiently?

1. Burn-in: The chain needs some number of steps to forget its initial starting position and converge to the target distribution. These initial steps are called the "burn-in" period and should be discarded from your final set of samples.
2. Trace Plots: A plot of the sample value at each iteration (the "trace") is the most basic diagnostic. It should look like a fuzzy caterpillar exploring the high-probability region, without long-term drifts or getting stuck.
3. Acceptance Rate: The fraction of proposed moves that are accepted is a key indicator of efficiency.
   • If the proposal step size is too small: Most moves will be accepted, but the walker explores the space very slowly. The trace plot will look highly correlated.
   • If the proposal step size is too large: Most moves will land in low-probability regions and be rejected. The walker will get stuck at the same point for many steps.
   • A rule of thumb is to tune the proposal step size to achieve an acceptance rate of around 20-50%.

4. Application: The 2D Ising Model

Let's apply the Metropolis algorithm to a problem central to condensed matter physics: simulating the 2D Ising model. Here, the "state" is a configuration of spins on a lattice, and the target probability distribution is the Boltzmann distribution, $P(\{\sigma \})\propto e^{-E(\{\sigma \})/k_{B}T}$.

The Metropolis algorithm for the Ising model is:

1. Start with a random lattice of N spins, ${\sigma }_i=\pm 1$.
2. Repeat for many steps:
   a. Pick a spin at a random site $i$.
   b. Propose flipping it: ${\sigma }_i\to -{\sigma }_i$.
   c. Calculate the change in energy, $\Delta E$, that this flip would cause.
   d. If $\Delta E<0$, accept the flip.
   e. If $\Delta E>0$, accept the flip with probability $e^{-\Delta E/k_{B}T}$.
3. After a burn-in period, collect measurements of quantities like total magnetization or energy to find their thermal averages.

This allows you to simulate the thermal equilibrium properties of the magnet at a given temperature $T$, including the phase transition from a disordered paramagnet to an ordered ferromagnet.

Conclusion

In this lesson, you have learned the fundamentals of the Metropolis-Hastings algorithm, one of the most important algorithms in computational science.

Key Takeaways:

• Monte Carlo methods use statistical sampling to solve high-dimensional problems where grid-based methods fail due to the curse of dimensionality.
• The Metropolis-Hastings algorithm generates a sequence of samples (a Markov chain) that are distributed according to a target probability distribution $P(\mathbf{x})$.
• Its core is a propose-accept/reject mechanism that satisfies detailed balance, ensuring convergence to the correct distribution.
• A key strength is its reliance only on the ratio of probabilities, making it applicable even when the distribution's normalization constant is unknown.
• The performance of the sampler depends critically on tuning the proposal distribution, with the acceptance rate being a key diagnostic.

Preview of the Next Lesson:

We have successfully built a tool for sampling classical probability distributions. Now, we will bring this tool back into the quantum realm. In the next lesson, Variational Monte Carlo (VMC), we will use the Metropolis algorithm to sample the probability distribution defined by a quantum wavefunction itself, $|\Psi (\mathbf{R}){|}^2$. This will allow us to calculate the ground state energy and other properties for many-body quantum systems far larger than what is possible with Exact Diagonalization.