# Transition Matrix Monte Carlo

## 1 Introduction

This paper presents both an approach to analyzing the information contained in configurations generated by general Monte Carlo simulations, and a closely related method of simulation that provides great flexibility and surprising efficiency. The approach uses information contained in the configurations about the set of possible changes on the next Monte Carlo step, which we encode in a “Transition Matrix.”<sup>-</sup> This information is complementary to a histogram analysis<sup>,</sup>. Indeed, it was originally developed as an improvement to the histogram approach. The discovery that the transition matrix, by itself, had definite advantages over the older methods has opened up the development of several new techniques.

Combining the Transition Matrix with an N-fold way simulation technique produces an extremely flexible and efficient approach to rather general Monte Carlo simulations. Because the same information is needed for the simulation and the analysis, the two methods work extremely efficiently together. In particular, the fact that this information is updated and used for every step of the N-fold way simulation enables contributions to the transition matrix to be made at every MC step, instead of after every sweep, as in cluster simulations.

Since no use is made of cluster methods that could be limited to non-frustrated systems, the techniques described in this paper are extremely general. They can be used directly with any system that has equally spaced energy levels, and, as is the case for histograms, can be generalized to systems with continuous symmetry with the use of binning.

In the following sections, we will define the transition matrix, describe how thermodynamic information is extracted from it, describe the N-fold way (as it is used with the transition matrix), describe the creation of generalized ensembles, and discuss the advantages of using these methods as part of parallel simulations similar to Replica Monte Carlo or Parallel Simulated Tempering.

## 2 The Transition Matrix

We will use the two-dimensional Ising model to demonstrate the transition matrix approach, although generalizations to higher dimensions and more complicated models is trivial. In particular, the generalization to a spin glass is obvious. The Ising model is given by the Hamiltonian

$$H=J\underset{i,j}{}\sigma _i\sigma _j,$$
(1)   
where $`\sigma _i`$ takes on the values $`\pm 1`$.

Although transition matrices can be defined for any kind of Monte Carlo simulation that depends only of the energies of the initial and final states, the primary transition matrix is defined in terms of the standard Monte Carlo dynamics at infinite temperature. For this dynamics, each spin is chosen with equal probability and “flipped” with probability one.

For every spin in a configuration, it is easy to calculate the change in energy of the configuration when that spin is reversed by simply counting the number of neighbors with the same sign. Denoting the number of sites that will produce an energy change $`\delta E`$ by $`n_{\delta E}`$, we define the “infinite temperature transition matrix” by

$$T_{E,\delta E}n_{\delta E}_E/N,$$
(2)   
where $`N`$ is the total number of sites and the average is taken over all configurations at energy $`E`$.

The largest eigenvalue of this matrix is unity, and the corresponding eigenvector is the density of states, $`W(E)`$.

$$\underset{\delta E}{}W(E\delta E)T_{E\delta E,\delta E}=W(E)$$
(3)   
The elements of the transition matrix must satisfy the condition of detailed balance,

$$W(E)T_{E,\delta E}=W(E+\delta E)T_{E+\delta E,\delta E},$$
(4)   
which requires the “TTT-identity”,

$$T_{E,\mathrm{\Delta }}T_{E+\mathrm{\Delta },\mathrm{\Delta }}T_{E+2\mathrm{\Delta },2\mathrm{\Delta }}=T_{E,2\mathrm{\Delta }}T_{E+2\mathrm{\Delta },\mathrm{\Delta }}T_{E+\mathrm{\Delta },\mathrm{\Delta }},$$
(5)   
where $`\mathrm{\Delta }`$ is the smallest allowable energy change.

For Potts models, the Ising model in three-dimensions, or other models with a greater number of possible energy changes, there are additional TTT-identities of the form,

$$T_{E,(m1)\mathrm{\Delta }}T_{E+(m1)\mathrm{\Delta },\mathrm{\Delta }}T_{E+m\mathrm{\Delta },m\mathrm{\Delta }}=T_{E,m\mathrm{\Delta }}T_{E+m\mathrm{\Delta },\mathrm{\Delta }}T_{E+(m1)\mathrm{\Delta },(m1)\mathrm{\Delta }}.$$
(6)   
The situation is only slightly more complicated if other quantities, like the magnetization or second-neighbor interactions, are introduced. In each case, the introduction of an additional energy change creates two new elements at each energy ($`+\delta E`$ and $`\delta E`$), along with one new identity (possibly involving four matrix elements), leaving one new independent variable.

Transition matrices corresponding to general Monte Carlo dynamics, $`\stackrel{ˇ}{T}_{E,\delta E}`$, can be constructed from the acceptance probabilities,

$$\stackrel{ˇ}{T}_{E,\delta E}=\{\begin{array}{cc}T_{E,\delta E}\times a_{E,\delta E}\hfill & \mathrm{for}\delta E0\hfill \\ T_{E,\delta E}\times a_{E,\delta E}+\mathrm{\Sigma }_{\delta E0}T_{E,\delta E}\times \left(1a_{E,\delta E}\right)\hfill & \mathrm{for}\delta E=0\hfill \end{array}.$$
(7)   
The corresponding probability distribution can then be constructed from the leading eigenvector of this matrix.

For example, the standard Metropolis acceptance probabilities are given by

$$a_{E,\delta E}=\mathrm{min}[1,\mathrm{exp}\left(\delta E/kT\right)],$$
(8)   
and the leading eigenvector of the corresponding transition matrix is the usual canonical probability distribution.

## 3 The N-Fold Way

Instead of choosing spins and determining whether to accept each move, it is possible to use the set of numbers, $`\left\{n_{\delta E}\right\}`$, to determine in advance what the probability of picking a spin that will change the energy by $`\delta E`$ and calculating its probability of acceptance. A class of spins with $`\delta E`$ is then picked with probability

$$\left(\frac{n_{\delta E}}{N}\right)\left(\frac{a_{E,\delta E}}{A}\right),$$
(9)   
where

$$A=\underset{\delta E}{}\left(\frac{n_{\delta E}}{N}\right)a_{E,\delta E}.$$
(10)   
A spin from that class is chosen and flipped. The contributions to the transition matrix (and any other quantity being computed) are weighted with a factor of $`\left(1/A\right)`$.

Since the set of numbers $`\left\{n_{\delta E}\right\}`$ are updated at every step for both the simulation and the transition matrix, little extra computer time is needed to record contributions to the transition matrix at every Monte Carlo step.

It is clear that any acceptance rates may be used in defining an ensemble for Monte Carlo simulations, as long as detailed balance is satisfied. Because of the reweighting of the contributions of each configuration, it is not even necessary for the acceptance rates to be normalized.

## 4 Generalized Ensembles

To create an ensemble with a desired probability distribution, $`P(E)`$, a necessary condition is

$$P(E)\times T_{E,\delta E}\times a_{E,\delta E}=P(E+\delta E)\times T_{E+\delta E,\delta E}\times a_{E+\delta E,\delta E}.$$
(11)   
For example, the multicanonical ensemble is characterized by $`P(E)=P(E+\delta E)`$, so that the condition on the acceptance rates is

$$\frac{a_{E,\delta E}}{a_{E+\delta E,\delta E}}=\frac{T_{E+\delta E,\delta E}}{T_{E,\delta E}}.$$
(12)   
Eqn. 4 relates this to the usual condition that the ratios of acceptance rates is equal to the ratio of the densities of state at the two energies.

Eqn. 12 is not sufficient to determine the acceptance rates uniquely. In addition to the usual acceptance rates given by the minimum of the ratio in Eqn. 12 or unity, both

$$a_{E,\delta E}=T_{E+\delta E,\delta E}$$
(13)   
and

$$a_{E,\delta E}=1/T_{E,\delta E}$$
(14)   
are valid. Many other options exist.

## 5 The “Equal-Hit” Ensemble

Although the multicanonical ensemble was designed to visit every energy level with equal probability, that is not really optimal when using the N-fold way. For low energies, the acceptance ratio is small and the N-fold way achieves a uniform $`P(E)`$ by visiting an energy level very few times, but weighting it with a large factor $`A^1_{E,Nf}`$. To scan all energy levels equally, it would be more appropriate to specify that the number of visits to each energy level, $`H(E)`$, is uniform. Note that the average of $`A^1`$ over the visits to energy level, $`E`$, with the N-fold way is not the same as the microcanonical average over configurations with that energy. Denoting the microcanonical average of a quantity, $`B`$, by, $`B_E`$, we have

$$B_E=\frac{BA^1_{E,Nf}}{A^1_{E,Nf}}.$$

If we take $`B=A`$, this gives

$$A^1_{E,Nf}=\frac{1}{A_E}.$$

Since the probability, $`P(E)`$, is given by the product of the number of times an N-fold way move ends at the energy $`E`$ times the average inverse acceptance ratio, $`A^1_{E,Nf}`$, the condition for an equal-hit ensemble with $`H(E)=H(E+\delta E)`$ is

$$A^1_{E,Nf}\times T_{E,\delta E}\times a_{E,\delta E}=A^1_{E+\delta E}\times T_{E+\delta E,\delta E}\times a_{E+\delta E,\delta E}.$$
(15)   
This gives the condition on the acceptance rates.

$$\frac{a_{E,\delta E}}{a_{E+\delta E,\delta E}}=\frac{A^1_{E+\delta E,Nf}\times T_{E+\delta E,\delta E}}{A^1_{E,Nf}\times T_{E,\delta E}}$$
(16)   
This condition can also be fulfilled in many ways, including

$$a_{E,\delta E}=A^1_{E+\delta E,Nf}\times T_{E+\delta E,\delta E}$$
(17)   
and

$$a_{E,\delta E}=A^1_{E+\delta E,Nf}/T_{E,\delta E}.$$
(18)   
## 6 Equilibration

At the beginning of a Monte Carlo simulation, there is usually little information about either the density of states or the transition matrix. For a canonical simulation at a given temperature, this does not cause a problem. However, for either a multicanonical simulation or an equal-hit simulation, such information is necessary. However, it turns out that the problem of dealing with the lack of information is greatly simplified by the transition matrices.

The primary approach consists of three stages.

For the first stage, we start with a random configuration and use one of the acceptance rates if it is known, otherwise we arbitrarily set it equal to unity. As the simulation in this stage progresses, we have a changing random walk that does not satisfy detailed balance and is biased towards states that have not yet been visited.

For the second stage, we take the approximate (and biased) transition matrix from the first stage, and impose the TTT-identities. In the second stage of the simulation, we use this approximate transition matrix to define the acceptance rates. Although the acceptance rates are not those of the ensemble that we really want to simulate, they do satisfy detailed balance by their construction, so that the transition matrix calculated in the second stage is not biased.

For the third stage, we use the estimate of the transition matrix (after imposing the TTT-identities) from the second stage to determine the acceptance rates. This algorithm represents a good approximation to the ensemble we wish to simulate, as well as satisfying detailed balance. The third stage is the main part of the simulation; it is where most of the computer time is used.

A modification of the third stage would be to continually update the transition matrix as the simulation progresses. This has the advantage that the simulation provides an increasingly accurate representation of the desired ensemble. The obvious disadvantage is that, strictly speaking, the algorithm no longer satisfies detailed balance. However, we have carried out some very long simulations on small systems with this modification, and found that the errors decrease with the square root of the length of the simulation, as expected. Note that the effect of any bias from the early part of the simulation is expected to decrease directly as the length of the simulation, rendering it unimportant. Therefore, we believe that the small violation of detailed balance is harmless, and will improve the efficiency of the algorithm.

The two-stage equilibration before the main simulation in the third stage requires relatively little computer time, so that this approach turns out to be more efficient than the standard multicanonical procedure, as well as being much simpler.

## 7 Parallel Simulations

A great advantage of a multicanonical simulation, which is retained by the equal-hit simulation, is that information over the entire range of energies can be obtained. This means that a single simulation provides the thermodynamic behavior at all temperatures – even negative temperatures (also known as the antiferromagnetic version of the model).

However, these advantages come with a price. The relaxation time to equilibrate a random walk is proportional to the square of its range. Since the total energy range is proportional to the number of particles, $`N`$, the relaxation time in units of MC steps is proportional to $`N^2`$. This does, indeed, turn out to be the case for the simulations described above.

We can greatly reduce this disadvantage by introducing parallel simulations. If we consider a number, $`l`$, of parallel simulations on independent replicas of the system, in which each replica is restricted to a range proportional to $`N/l`$, the relaxation time for each replica is only proportional to $`\left(N/l\right)^2`$. Taking the increased number of simulations into account, we are left with a relaxation time in units of MC steps proportional to $`N^2/l`$. The relaxation time in units of MC steps/site is proportional to $`N/l`$.

Note that this improvement occurs even with a serial machine. Implementation on a parallel machine is trivial.

We can further improve the algorithm by exchanging replicas, as introduced in the Replica Monte Carlo method<sup>-</sup> and later rediscovered as parallel simulated tempering.

Break up the total energy range into pieces with equal width. Beginning with the lowest energies, restrict the first replica to the first two pieces of the energy spectrum. The second replica is restricted to the second and third pieces, and so on. Replica $`n`$ and replica $`n+1`$ therefore overlap for half of their range of energies. After simulating for a few relaxation times (proportional to $`N^2/l`$), all neighboring replicas that are in the overlapping half of their ranges are exchanged. Since we are either using a multicanonical simulation with single spin flips, or an equal-hit ensemble with the N-fold way, the acceptance probability for such exchanges is unity.

This combination of replica exchange with an equal-hit (or multicanonical) simulation provides improvements over both. Although this approach has been discussed in terms of transition matrix Monte Carlo, its use is much broader. For example, it could easily be applied to multicanonical simulations of biological molecules, as we intend to do in future work.