Monte Carlo Importance Sampling¶

Let us introduce the concept of importance sampling method by application to classical many-particle system (like Ising model or classical gas).

The basic idea of Monte Carlo Simulation: The simulation is performed by random walk through very large configuration space. The probability to make a move has to be such that the system gets to thermal equilibrium (and remains in thermal equilibrium) at certain temperature $T$ (is usually a parameter) after a lot of Monte Carlo steps.

The basic idea of simulated annealing: Slowly decreasing the temperature of the system leads to the ground state of the system. When using this type of cooling down by Monte Carlo, one can find a global minimum of a general minimization problem.

MC Importance Sampling by Markov chain¶

If a configuration in phase space is denoted by $X$, the probability for configuration according to Boltzman is \begin{equation} \rho(X)\propto e^{-\beta E(X)}\qquad \beta=\frac{1}{T} \end{equation}

How to sample over the whole phase space for a general problem? How to generate configurations?

  • Brute force: generate a truly random configuration $X$ and accept it with probability $e^{-\beta E(X)}$ where all $E>0$. Successive $X$ are statistically independent.

    VERY INEFFICIENT

  • Markov chain: Successive configurations $X_i$, $X_{i+1}$ are NOT statistically independent but are distributed according to the choosen disribution (such as Boltzman distribution).

    Can be made VERY EFFICIENT


What is the difference between Markov chain and uncorrelated sequence?
  • Truly random or uncorrelated sequence of configurations satisfies the identity $$P(X_1,X_2,\cdots,P_{X_N})=P_1(X_1)P_1(X_2)\cdots P_1(X_{N})$$
  • Markov chain satisfies the equation
    $$P(X_1,X_2,\cdots,P_{X_N})=P_1(X_1)T(X_1\rightarrow X_2) T(X_2\rightarrow X_3)\cdots T(X_{N-1}\rightarrow X_N)$$ where the transition probabilities $T(X\rightarrow X')$ are normalized $$\sum_{X'} T(X \rightarrow X')=1$$

We want to construct a Markov chain whose distribution approaches the target distribution $$\rho(X)$$ after sufficiently many steps, independently of the initial configuration.

For a finite-state Markov chain, the following is needed:

  • Stationarity: The transition matrix preserves the desired distribution $\rho$. Global balance is the required condition; detailed balance is a convenient sufficient condition (see below).

  • Irreducibility: Every state with nonzero target probability can eventually be reached from every other such state.

  • Aperiodicity: The chain is not trapped in a fixed temporal cycle. Together, irreducibility and aperiodicity are often summarized by saying that the chain is ergodic. They ensure that, after sufficiently many steps we get convergence: $$p_t(X) \longrightarrow \rho(X)$$

Global balance vs Detail balance¶

We need to find transition probability $T(X\rightarrow X')$ which leads to a given stationary distribution $\rho(X)$ (in this case $\rho(X)\propto e^{-\beta E(X)}$).

The probability for $X$ decreases, if system goes from $X$ to any other $X'$: $\Delta \rho(X)=-\sum_{X'}\rho(X)T(X\rightarrow X')$ and increases if $X$ configuration is visited from any other state $X'$: $\Delta \rho(X) = \sum_{X'}\rho(X')T(X'\rightarrow X)$. The step (time) difference of the probability $X$ is therefore \begin{equation} \rho(X,t+1)-\rho(X,t) = -\sum_{X'}\rho(X)T(X\rightarrow X')+\sum_{X'}\rho(X')T(X'\rightarrow X) \end{equation} We look for stationary solution, i.e., $\rho(X,t+1)-\rho(X,t)=0$ and therefore \begin{equation} \sum_{X'}\rho(X)T(X\rightarrow X')=\sum_{X'}\rho(X')T(X'\rightarrow X) \end{equation}

Solving the global-balance equation directly is generally difficult. A particular solution (a simpler and sufficient) is detailed balance solution: \begin{equation} \rho(X)T(X\rightarrow X')=\rho(X')T(X'\rightarrow X) \end{equation}

To construct algorithm, we divide transition prob. $T(X\rightarrow X')=\omega_{X\rightarrow X'}A_{X\rightarrow X'}$:

  • trial step probability $\omega_{X\rightarrow X'}$.

    Many times it is symmetric, i.e., $\omega_{X\rightarrow X'}=\omega_{X'\rightarrow X}$. (for example spin flip in ising: $\omega_{XX'}$ is $1/L^2$ if $X$ and $X'$ differ for a single spin flip and zero otherwise ).

  • acceptance probability $A_{X\rightarrow X'}$

(for example accepting or rejecting new configuration with probability proportional to $min(1,\exp(-\beta(E(X')-E(X))))$).

Detail balance condition becomes $$\frac{A_{X\rightarrow X'}}{A_{X'\rightarrow X}}=\frac{\omega_{X' \rightarrow X}\;\rho(X')} {\omega_{X\rightarrow X'}\;\rho(X)}$$

Metropolis chooses \begin{eqnarray} \begin{array}{lc} A_{X\rightarrow X'}=1 & \;\mathrm{if}\;\; \omega_{X'\rightarrow X} \rho(X')> \omega_{X\rightarrow X'} \rho(X)\\ A_{X\rightarrow X'}= \frac{\omega_{X'\rightarrow X}\,\rho(X')}{\omega_{X\rightarrow X'}\,\rho(X)} & \;\mathrm{if}\;\;\omega_{X'\rightarrow X}\rho(X')<\omega_{X\rightarrow X'}\;\rho(X). \end{array} \end{eqnarray} Obviously, this acceptance probability satisfies detail balance condition, and therefore can lead to desired Markov chain with stationary probability $p_t(X)\propto \rho(X)$ for long times, provided we also have irreducibility and aperiodicity.

Metropilis is not the only possible choice. For example, so called Heat-bath is an alternative. In Heath-bath we choose: \begin{eqnarray} {A_{X\rightarrow X'}}=\frac{\omega_{X' \rightarrow X}\;\rho(X')} {\omega_{X\rightarrow X'}\;\rho(X)+\omega_{X' \rightarrow X}\;\rho(X')} \\ {A_{X'\rightarrow X}}=\frac{\omega_{X \rightarrow X'}\;\rho(X)} {\omega_{X\rightarrow X'}\;\rho(X)+\omega_{X' \rightarrow X}\;\rho(X')} \end{eqnarray} The ratio for detailed balance is clearly satisifed.

To summarize Metropolis algorithm

  • $T(X\rightarrow X')=\omega_{X\rightarrow X'}A_{X\rightarrow X'}$ (we have to consider trial step probability)
  • $\sum_{X'}\omega_{X\rightarrow X'}=1$ and $\omega_{X\rightarrow\cdots\rightarrow X'}>0$ for all $X, X'$ after finite number of steps. The proposal and acceptance rules must allow every state with nonzero $\rho(X)$ to be reached from every other such state, so that the resulting Markov chain is irreducible. The chain must also be aperiodic.
  • $A_{X\rightarrow X'}=min(1,\frac{\rho(X')\omega_{X'\rightarrow X}}{\rho(X)\omega_{X\rightarrow X'}})$ (acceptance step probability for Metropolis)

How to accept a step with probability $A_{XX'}<1$? One can generate a random number $r\in[0,1]$ and accept the step if $r<A_{XX'}$.

Keep in mind:

  • Configurations that are generated by Markov chain are correlated. The theory guarantees that we arrive at invariant distribution $\rho$ after a long time.
  • Two configurations are statistically independent only if they are far apart in the Markov chain. This distance is called it correlation time, and is usually very long.
  • Be careful: To meassure distance in Markov chain, every step counts, not only successful steps.

Most of the time we are simulating configurations which correspond to the partition function, so that the probability to find a given state $X$ in the Markov's chain is equal to $\rho(X) = \frac{e^{-\beta E(X)}}{\sum_X e^{-\beta E(X)}}=\frac{\Delta Z(X)}{Z} $.

The average of any quantity is then calculated by $$\overline{A}\equiv \sum_X \rho(X) A(X) \rightarrow \frac{1}{n-n_0}\sum_{i>n_0}^n A_i$$ where $n_0$ steps are used to "warm-up". This is because the configurations in the Markov's chain are distributed according to $\rho(X)$.

We can also try compute the following quantity $$\sum_X \rho(X)(A(X)-\overline{A})^2 \rightarrow^? \frac{1}{n-n_0}\sum_{i>n_0}^n (A_i-\overline{A})^2.$$ This is A WRONG ESTIMATION OF THE ERROR OF MC. The error is much bigger than this estimation, because configurations $X$ are correlated! Imagine the limit of large correlations when almost all values $A_i$ are the same (very slowly changing configurations). We would estimate that standard deviation is zero regardless of the actual error!

To compute standard deviation, we need to group consecutive measurements into blocks whose length $N_0$ is larger than the correlation time. We then do statistics on the bins: \begin{eqnarray} B_l &=& \frac{1}{N_0}\sum_{i=l N_0 + 1}^{i=(l+1) N_0} A_i\\ \overline{B}&=&\overline{A}\\ \sigma^2_{\overline{A}} &=& \frac{1}{M(M-1)}\sum_{l=0}^{M-1}(B_l-\overline{B})^2\\ N &=& N_0 M \end{eqnarray} where we took into account that $\overline{A}=\overline{B}$. The correlation time is not very easy to estimate. Maybe the best algorithm is to compute $\sigma^2$ for few different $N_0$ and as long as $\sigma^2$ is increasing with $N_0$, the correlation time is still larger than $N_0$. When $\sigma^2$ stops changing with increasing $N_0$, we reached correlation time and $\sigma_{\overline{A}}$ is a good estimation the standard .

In [ ]: