The goal is to construct a Markov Chain $\{X_k\}$ with a stationary distribution $\pi$.

Metropolis Algorithm

First, we construct the proposal transition matrix $Q=[q(i,j)]$ where $Q$ is symmetric. Then, choose an initial state $X_0\in S$ randomly. Given $X_k=i\in S$,

1. Sample $Y\in S$ from $q(i,-)$, i.e. $\mathbb{P}(Y=y|X_k=i)=q(i,y)$ for all $y\in S$. We assume $Y=j$ is observed.
2. Let $$\alpha ={\alpha }_{ij}=\min \left(1,\frac{\pi (j)}{\pi (i)}\right)\in [0,1]$$ and generate $U\sim \mathcal{U}[0,1]$.
3. If $U\le \alpha$ then accept $Y$ and $X_{k+1}\leftarrow Y$.
4. If $U>\alpha$ then reject $Y$ and $X_{k+1}\leftarrow X_k$.

This means we accept $Y$ with probability $\alpha$. The output is a sequence $X_0,X_1,\dots$

Remarks About Computation

• $\{X_n{\}}_{n=0}^{\infty }$ is a Markov Chain with ${\lim }_{n\to \infty }{X}_n\sim \pi$.
• We can choose $Q$ to be simple, like uniform jumps, where $q(i,j)=1/N$.
• We only need to now the ratio $\pi (j)/\pi (i)$, so we don’t need to know the normalization value $Z$.
• $Y$ is a proposal state.
• In practice, $\alpha ={\alpha }_{ij}$ is easy to calculate.
• This algorithm works for non-symmetric $Q$ and $$\alpha ={\alpha }_{ij}=\min \left(1,\frac{\pi (j)q(j,i)}{\pi (i)q(i,j)}\right)$$

Example 1 (1D Ising Model)

Recall the Ising Model, where we have $m$ particles each with spin $+$ or $-$. So,

$$S=\{\sigma =({\sigma }_1,\dots ,{\sigma }_m)|{\sigma }_i=\pm 1,1\le i\le m\}=\{-1,+1{\}}^m$$

or the set of length $m$ strings of alphabet $-1,+1$. Obviously $|S|=2^m$.

Let

$$H(\sigma )=-J\sum _{i=1}^{m-1}{\sigma }_i{\sigma }_{i+1}-h\sum _{i=1}^m{\sigma }_i$$

where $h>0$¹. $H(\sigma )$ represents the “energy” of state $\sigma$.

Our target distribution² is

$${\pi }_{\beta }(\sigma )=\frac{1}{Z_{\beta }}\exp \left(-\beta H(\sigma )\right)\forall \sigma \in S$$

where $\beta >0$ is the inverse temperature $(\beta =1/k_{\beta }T)$, and $Z_{\beta }$ is normalization (partition function). This is hard to compute.

Note that a large ${\pi }_{\beta }(\sigma )⟺$ small $H(\sigma )⟺$ many ${\sigma }_j=+1$ and many

$${\sigma }_j{\sigma }_{j+1}=\{\begin{array}{ll}>0& \text{if }J>0\\ <0& \text{if }J<0\end{array}$$

Here, we can apply the Metropolis Algorithm. Choose $Q$ given by the random spin flips. For example, if $m=2$:

$"\\usepackage{tikz-cd}\n\\usetikzlibrary{arrows.meta,calc}\n\\tikzset{curve/.style={settings={#1},to path={(\\tikztostart) .. controls ($(\\tikztostart)!\\pv{pos}!(\\tikztotarget)!\\pv{height}!270:(\\tikztotarget)$) and ($(\\tikztostart)!1-\\pv{pos}!(\\tikztotarget)!\\pv{height}!270:(\\tikztotarget)$) .. (\\tikztotarget)\\tikztonodes}}, settings/.code={\\tikzset{quiver/.cd,#1}\\def\\pv##1{\\pgfkeysvalueof{/tikz/quiver/##1}}}, quiver/.cd,pos/.initial=0.35,height/.initial=0}\n\\begin{document}\n$

source code

we get the following. If $S=\{++,+-,-+,--\}$, then the transition matrix is:

$$P=\left[\begin{array}{cccc}0& 1/2& 1/2& 0\\ 1/2& 0& 0& 1/2\\ 1/2& 0& 0& 1/2\\ 0& 1/2& 1/2& 0\end{array}\right]$$

which is symmetric. Choose $X_0\in S$ given $X_k=({\sigma }_1,\dots ,{\sigma }_m)\in S$. Some notes,

• $J$ represents the interaction strength between particles. It can be treated as a constant.
• $h$ is the external magnetic field on the particles. We can treat it as a constant.

Algorithm:

1. Choose a site $i$ for $1\le i\le m$ uniformly at random among $1,2,\dots ,m$.
2. Flip the $i$th site ${\sigma }_i\to -{\sigma }_i$ to obtain the proposal state: $$\hat{\sigma }=({\sigma }_1,\dots ,{\sigma }_{i-1},-{\sigma }_i,{\sigma }_{i+1},\dots ,{\sigma }_m)\in S$$ Note that $q(\sigma ,\hat{\sigma })=q(\hat{\sigma },\sigma )=1/m$. This shows $Q$ is symmetric.
3. Compute $$\begin{array}{rl}\frac{\pi (\hat{\sigma })}{\pi (\sigma )}& =\exp \left[-\beta (H(\hat{\sigma })-H(\sigma ))\right]\\ & =\exp (-\beta \Delta H)\end{array}$$ Recall that we do this to avoid the partition function.
4. Three cases for $\Delta H$:
   1. If $1<i<m$, (we picked a particle in the middle of the chain) then $$\begin{array}{rl}H(\hat{\sigma })-H(\sigma )& =-J({\sigma }_{i-1}(-{\sigma }_i)-{\sigma }_i{\sigma }_{i+1})+h{\sigma }_i-\left[-J({\sigma }_{i-1}{\sigma }_i+{\sigma }_i{\sigma }_{i+1})-h{\sigma }_i\right]\\ & =2J{\sigma }_i({\sigma }_{i-1}+{\sigma }_{i+1})+2h{\sigma }_i\end{array}$$
   2. If $i=1$ (the left boundary), then $$\begin{array}{rl}H(\hat{\sigma })-H(\sigma )& =-J(-{\sigma }_1){\sigma }_2+h{\sigma }_1-\left[-J{\sigma }_1{\sigma }_2-h{\sigma }_1\right]\\ & =2J{\sigma }_1{\sigma }_2+2h{\sigma }_1\end{array}$$
   3. If $i=m$ (the right boundary), similarly $$\begin{array}{r}H(\hat{\sigma })-H(\sigma )=2J{\sigma }_{m-1}{\sigma }_m+2h{\sigma }_m\end{array}$$ Thus, $$\Delta H=H(\hat{\sigma })-H(\sigma )=\{\begin{array}{ll}2J{\sigma }_i({\sigma }_{i-1}+{\sigma }_{i+1})+2h{\sigma }_i& 1<i<m\\ 2J{\sigma }_1{\sigma }_2+2h{\sigma }_1& i=1\\ 2J{\sigma }_{m-1}{\sigma }_m+2h{\sigma }_m& i=m\end{array}$$
   4. We can then accept $\hat{\sigma }$ with probability $$\alpha =\min \left(1,\frac{\pi (\hat{\sigma })}{\pi (\sigma )}\right)=\min \left(1,e^{-\beta \Delta H}\right)$$
   5. Thus,
      1. If $\Delta H\le 0$ then $\alpha =1$ and set $X_{k+1}\leftarrow \hat{\sigma }$. (Here, the system lost energy, i.e. became more stable, so the second term is positive and thus greater than $1$. We will accept).
      2. If $\Delta H>0$ then $\alpha =\exp (-\beta \Delta H)$. (This is the opposite of above.)
         1. Generate $U\sim \mathcal{U}[0,1]$.
         2. If $U\le \exp (-\beta \Delta H)$, accept $\hat{\sigma }$ and $X_{k+1}\leftarrow \hat{\sigma }$
         3. If $U>\exp (-\beta \Delta H)$, reject $\hat{\sigma }$ and $X_{k+1}\leftarrow \sigma$.

Why does it work?

• In step 2, we do not calculate $H$ since it is computationally inefficient. However, we can calculate the change, $\Delta H$. Every term in $H(\sigma )$ cancels out but the changed term. This means we only need to find the sum of a few terms.
• The Ising model only involves interactions with the nearest neighbors. This is why we see changes in ${\sigma }_{i-1}$ and ${\sigma }_{i+1}$, reducing the computations needed from $O(m)$ to $O(1)$.

Theorem (Metropolis Satisfies Detailed Balance)

1. The transition matrix $P=[p(i,j)]$ of the Metropolis Chain is $$p(i,j)=\{\begin{array}{ll}q(i,j)\cdot \min \left(1,\frac{\pi (j)}{\pi (i)}\right)& i\ne j\\ 1-{\sum }_{k\in S\setminus \{i\}}q(i,k)\cdot \min \left(1,\frac{\pi (j)}{\pi (i)}\right)& i=j\end{array}$$
   1. The idea in case $1$ is that we propose moving from $i\to j$ with probability $q(i,j)$ in Metropolis, but we accept this movement $\alpha (i,j)$ times. This is just the probability of these two events happening together.
   2. In case $2$, accounts for the probability of rejection. Indeed, $p(i,i)$ is the sum of chances we tried to move but failed.
2. $\pi$ and $P$ satisfy detailed balance: $$\pi (i)\cdot p(i,j)=\pi (j)\cdot p(j,i)$$ for all $i,j\in S$. (Hence $\pi$ is a stationary distribution.)
3. If $\forall i,j\in S,q(i,j)>0$ and $\forall i\in S,p(i,i)>0$, then $\pi$ is unique, and $$\underset{n\to \infty }{\lim }{\pi }_0{P}^n=\pi$$

Proof:

$(1)$: For $i\ne j,$

$$\begin{array}{rl}p(i,j)& =\mathbb{P}(\text{propose }j\mid i)\cdot \mathbb{P}(\text{accept }j\mid \text{proposed }j)\\ & =q(i,j)\cdot \min \left(1,\frac{\pi (j)}{\pi (i)}\right)\end{array}$$

Then since $P$ is a stochastic matrix, for any state $i$, the sum ${\sum }_{j\in S}p(i,j)=1$ by law of total probability. This implies

$$\begin{array}{rl}p(i,i)& =1-\sum _{\begin{array}{c}k\in S\\ k\ne i\end{array}}p(i,k)\\ & =1-\sum _{\begin{array}{c}k\in S\\ k\ne i\end{array}}q(i,k)\cdot \min \left(1,\frac{\pi (k)}{\pi (i)}\right)\end{array}$$

which represents the “rejection probability”. The chain stays at state $i$ if either

1. The proposal $q(i,i)$ suggests staying at $i$,
2. or if the the proposal $j\ne i$ was made, but was rejected with probability $1-\alpha (i,j)$.

$(2)$: We want to show that

$$\forall i,j\in S,\pi (i)\cdot p(i,j)=\pi (j)\cdot p(j,i)$$

Trivially, this is true for $i=j$. Suppose $i\ne j$. Then

$$\begin{array}{rl}\pi (i)\cdot p(i,j)& =\pi (i)\cdot q(i,j)\cdot \min \left(1,\frac{\pi (j)}{\pi (i)}\right)\\ & =q(i,j)\cdot \min (\pi (i),\pi (j))\\ & =q(j,i)\cdot \min (\pi (j),\pi (i))\\ & =\pi (j)\cdot q(j,i)\cdot \min \left(1,\frac{\pi (i)}{\pi (j)}\right)\\ & =\pi (j)\cdot p(j,i)\end{array}$$

However, this is only true because we let $q(i,j)=q(j,i)$, i.e. the matrix $Q$ must be symmetric. Thus Metropolis satisfies detailed balance, and $\pi$ is stationary.

$(3)$: The condition $q(i,j)>0$ ensures the chain is irreducible (any state can reach any other in one step). The condition $p(i,i)>0$ ensures the chain is aperiodic. An irreducible, aperiodic chain on a finite state space is regular. By the FTMC, a regular matrix has a unique stationary distribution $\pi$ and ${\lim }_{n\to \infty }{\pi }_0{P}^n=\pi$.

Remark (Metropolis on Ising Model)

Using Metropolis on the Ising Model, we had

$$S=\{\sigma =({\sigma }_1,\dots ,{\sigma }_m),{\sigma }_j\in \{+1,-1\}\}$$

where sign flips would give us

$$({\sigma }_1,\dots ,{\sigma }_j,\dots ,{\sigma }_m)\mapsto ({\sigma }_1,\dots ,-{\sigma }_j,\dots ,{\sigma }_m)$$

which defines our proposal probability as

$$q(\sigma ,{\sigma }^{\prime })=\{\begin{array}{ll}1/m& \sigma \to {\sigma }^{\prime }\text{ is 1 flip}\\ 0& \text{otherwise}\end{array}$$

We can show that $q$ gives a regular Markov Chain. The idea is that any $2$ states in $S$ are at most $m$ flips apart.

Metropolis-Hastings Algorithm

We can run Metropolis with an asymmetric proposal matrix $Q$. The purpose of this is to relax the condition that $Q$ must be symmetric in Theorem (Metropolis Satisfies Detailed Balance). This allows us to correct the model over-sampling states frequently proposed by $Q$ regardless of the actual probability in the target distribution $\pi$.

1. Given $X_k\in S$, choose $Y$ from $q(X_k,\cdot )$.
2. Accept $Y$ with probability

$$\alpha =\min \left(1,\frac{\pi (Y)}{\pi (X_k)}\cdot \frac{q(Y,X_k)}{q(X_k,Y)}\right)$$

The output is a Markov Chain of $X_0,X_1,\cdots \in S$.

Convergence is hard to prove. Heuristically, consider that $R(n)=\text{Cov}(X_n,X_0)$. So as $n\to \infty$, $R(n)\to 0$.

Footnotes

1. This is the Hamiltonian. ↩

2. This is the Boltzmann Distribution. ↩