Kubo formula

2013-12-12T14:03:20Z

138.246.2.247:

'''Multicanonical ensemble''' (also called '''multicanonical sampling''' or '''flat histogram''') is a sampling technique that uses the [[Metropolis–Hastings algorithm]] to compute integrals where the integrand has a rough landscape with multiple local minima. It samples states according to the inverse of the [[density of states]],<ref name=Berg/> which has to be known a priori or be computed using other techniques like the [[Wang and Landau algorithm]].<ref name=Landau/> Multicanonical sampling is an important technique for spin systems like the [[Ising model]] or [[spin glass]]es.<ref name=Berg/><ref name=Newmann/><ref name=Dayal/>

==Motivation==

In systems with a large number of degrees of freedom, like spin systems, [[Monte Carlo integration]] is required. In this integration, [[importance sampling]] and in particular the [[Metropolis algorithm]], is a very important technique.<ref name=Newmann/> However, the Metropolis algorithm samples states according to <math>\exp(-\beta E)</math> where beta is the inverse of the temperature. This means that an energy barrier of <math>\Delta E</math> on the energy spectrum is exponentially difficult to overcome.<ref name=Berg/> Systems with multiple local energy minima like the [[Potts model]] become hard to sample as the algorithm gets stuck in the system's local minima.<ref name=Newmann/> This motivates other approaches, namely, other sampling distributions.

==Overview==

Multicanonical ensemble uses the Metropolis–Hastings algorithm with a sampling distribution given by the inverse of the density of states of the system, contrary to the sampling distribution <math>\exp(-\beta E)</math> of the Metropolis algorithm.<ref name=Berg/> With this choice, on average, the number of states sampled at each energy is constant, i.e. it is a simulation with a "flat histogram" on energy. This leads to an algorithm for which the energy barriers are no longer difficult to overcome. Another advantage over the Metropolis algorithm is that the sampling is independent of the temperature of the system, which means that one simulation allows the estimation of thermodynamical variables for all temperatures (thus the name "multicanonical": several temperatures). This is a great improvement in the study of first order [[phase transition]]s.<ref name=Berg/>

The biggest problem in performing a multicanonical ensemble is that the density of states has to be known ''a priori''.<ref name=Landau/><ref name=Newmann/> One important contribution to multicanonical sampling was the [[Wang and Landau algorithm]], which asymptotically converges to a multicanonical ensemble while calculating the density of states during the convergence.<ref name=Landau/>

The multicanonical ensemble is not restricted to physical systems. It can be employed on abstract systems which have a cost function ''F''. By using the density of states with respect to F, the method becomes general for computing higher dimensional integrals or finding local minima.<ref name=Lee/>

==Motivation==

Consider a system and it phase-space <math>\Omega</math> characterized by a configuration <math>\boldsymbol{r}\in \Omega</math> and a "cost" function ''F'' from the system's phase-space to a one dimensional space <math>\Gamma</math>: <math>F(\Omega) = \Gamma = [\Gamma_\min, \Gamma_\max]</math>, the spectrum of ''F''.

{| class="wikitable1" width=300px
|-
|
{{show
|example:
|
The [[Ising model]] with ''N'' sites is an example of such a system; the phase-space is a discrete phase-space defined by all possible configurations of ''N'' spins <math>\boldsymbol{r} = (\sigma_1, \ldots, \sigma_i, \ldots, \sigma_N)</math> where <math>\sigma_i\in\{-1,1\}</math>. The cost function is the [[Hamiltonian (quantum mechanics)|Hamiltonian]] of the system:

: <math>H(\boldsymbol{r}) = - \sum_{\langle i~j\rangle} J_{ij} (1 - \sigma_i \sigma_j),</math>

where <math><i,j></math> is the sum over neighborhoods and <math>J_{ij}</math> is the interaction matrix.

The energy spectrum is <math>\Gamma = [E_\min, E_\max]</math> which, in this case, depends on the particular <math>J_{ij}</math> used. If all <math>J_{ij}</math> are 1 (the ferromagnetic Ising model), <math>E_\min = 0</math> (e.g. all spins are 1.) and <math>E_\max = 2 D N </math> (half spins are up, half spins are down). Also notice that in this system, <math>\Gamma \in \mathbb{Z}</math>

}}
|}

The computation of an average quantity <math>\langle Q\rangle</math> over the phase-space requires the evaluation of an integral:

: <math>\langle Q\rangle = \frac{1}{V}\int_\Omega Q(\boldsymbol{r}) p(\boldsymbol{r}) \,d\boldsymbol{r}</math>

where <math>p(\boldsymbol{r})</math> is weight of each state per volume,

: <math>p(\boldsymbol{r}) = \frac{1}{V}\int_\Omega p(\boldsymbol{r}) \,d\boldsymbol{r}.</math>

The density of states in respect with ''F'' is given by

: <math>\rho(f) = \frac{1}{V}\int_\Omega \delta(f - F(\boldsymbol{r})) \,d\boldsymbol{r}</math>

which means that if both ''Q'' and ''p'' do not depend on the particular state but only on the particular F's value of the state <math>F(\boldsymbol{r}) = F_\boldsymbol{r}</math>,

: <math>Q(\boldsymbol{r}) = Q(F_\boldsymbol{r}), p(\boldsymbol{r}) = p(F_\boldsymbol{r}),</math>

the formula for <math>\langle Q\rangle</math> can be integrated over ''f'' by adding a [[dirac delta function]],

: <math>
\begin{align}
\langle Q\rangle & = \int_{\Gamma_\min}^{\Gamma_\max} \int_\Omega Q(F_\boldsymbol{r}) p(F_\boldsymbol{r}) \delta(f - F_\boldsymbol{r}) \,d\boldsymbol{r} \,df \\
& = \int_{\Gamma_\min}^{\Gamma_\max} Q(f) p(f) \int_\Omega \delta(f - F_\boldsymbol{r}) \,d\boldsymbol{r} \, df \\
& = \int_{\Gamma_\min}^{\Gamma_\max} Q(f) p(f) \rho(f) \, df \\
\end{align}
</math>

i.e. the knowledge of the density of states allows the computation of averages over ''F'' using a one-dimensional integral instead of a multidimensional integral as it is the projection of the number of states on <math>\Omega</math> in <math>\Gamma</math>.

{| class="wikitable1" width=300px
|-
|
{{show
|example:
|
A system in contact with a heat bath at inverse temperature <math>\beta</math> is a clear example for computing this kind of integral. For instance, the mean energy of the system is weighted by the [[Boltzmann factor]]:

:<math>\langle E\rangle = \frac{1}{V}\int_\Omega H(\boldsymbol{r}) \frac{e^{-\beta H(\boldsymbol{r})}}{Z} \, d\boldsymbol{r}</math>

where

: <math>Z = \frac{1}{V}\int_\Omega e^{-\beta H(\boldsymbol{r})} \, d\boldsymbol{r}.</math>

The density of states <math>\rho(E)</math> can used to compute <math>\langle E\rangle</math>:

: <math>\langle E\rangle = \int_{E_\min}^{E_\max} E \frac{e^{-\beta E}}{Z} \rho(E) \, dE </math>

}}
|}

Because the number of states can be high for systems with high number of degrees of freedom, an analytical expression can be hard to obtain and the computation of <math>\langle Q\rangle</math> can be expensive. Typically, because the problem is a multidimensional integral, [[Monte Carlo integration]] is normally employed. On the simplest formulation, the method chooses ''N'' uniform random states <math>\boldsymbol{r}_i\in \Omega</math>, and uses the [[estimator]]

:<math>\overline{Q}_N = \frac{1}{N}\sum_{i = 0}^N Q(\boldsymbol{r}_i) p(\boldsymbol{r}_i)</math>

for computing <math>\langle Q\rangle</math>: <math>\overline{Q}_N</math> converges almost surely to <math>\langle Q\rangle</math> by the [[Law_of_large_numbers#Strong_law|strong law of large numbers]]:

:<math>\lim_{N\rightarrow\infty} \overline{Q}_N = \langle Q\rangle.</math>

One typical problem of this convergence is that the variance of ''Q'' can be very high, which leads to a high computational effort to achieve reasonable results.

{| class="wikitable1" width=300px
|-
|
{{show
|example
|
On the previous example, the states that mostly contribute to the integral are the ones with low energy. If the states are sampled uniformly, on average, the number of states which are sampled with energy ''E'' is given by the density of states. This density of states can be centered far away from the energy's minima and thus the average can be difficult to obtain.

}}
|}

To improve this convergence, the [[Metropolis–Hastings algorithm]] was proposed. Generally, Monte Carlo methods' idea is to use [[importance sampling]] to improve the convergence of the estimator <math>\overline{Q}_N</math> by sampling states according to an arbitrary distribution <math>P(\boldsymbol{r})</math> (notice the capital ''P'', different from ''p''), and use the appropriate estimator:

: <math>\overline{Q}_N = \frac{1}{X}\sum_{i = 0}^N Q(\boldsymbol{r}_i) P^{-1}(\boldsymbol{r}_i) p(\boldsymbol{r}_i)</math>

where <math>X = \sum_{i = 0}^N P^{-1}(\boldsymbol{r}_i)</math>.

Notice that when ''P'' is a uniform distribution, this estimator equals the one used on a uniform sampling, as it should.

One important choice of ''P'' is to define an arbitrary temperature, and use it equals to [[Boltzmann factor]] with the energy being the cost function:

: <math>P(\boldsymbol{r}) = p_\mathrm{Boltzmann}(\boldsymbol{r}) = \frac{e^{-\beta F(\boldsymbol{r})}}{\int_\Omega \, d\boldsymbol{r} e^{-\beta F(\boldsymbol{r})}}</math>

I.e. the lower the cost function of a particular state, the more likely it is to be sampled.

Historically, this occurred because the [[Equation of State Calculations by Fast Computing Machines|original idea]]<ref name=Metropolis/> was exactly to use [[Metropolis–Hastings algorithm]] to compute averages on a system in contact with a heat bath where the weight is given by the Boltzmann factor. On these systems, the choice of states according to it led to a considerable improvement on studying physical systems.<ref name=Newmann/>

However, it is not true that the sampling distribution must equals the weight distribution. One reason for this is that [[Metropolis–Hastings algorithm]] fails to converge when the cost function has multiple minima.<ref name=Berg/> The reason is that the algorithm uses a random walk with local steps. I.e. the random walk normally performs steps whose energy difference is of order 1. This means that the computational cost to the algorithm leave a specific region with a local minimum exponentially increases with the cost function's value of the minimum. I.e. the deeper the minimum, the more time the algorithm spends there, and harder (exponentially with deep) it will leave.

This is the motivation to introduce a multicanonic ensemble. The idea is to avoid becoming stuck on local minima of the cost function by making them "invisible" to the sampling technique.

==Multicanonic ensemble==

A multicanonic ensemble is choosing the sampling distribution used in the importance sampling to be

: <math>P(\boldsymbol{r}) = \frac{1}{\rho(F(\boldsymbol{r}))}</math>

where

: <math> \rho(f) = \frac{1}{V}\int_\Omega \delta(F(\boldsymbol{r}) - f) \, d\boldsymbol{r}</math>

is the density of states of the system with respect to the cost function (''V'' is the phase-space volume). The consequence of this choice is that the probability to sample a state and it has cost function ''f'' is

: <math> P(f) = \frac{1}{f_\max - f_\min}\int_\Omega \delta(f - F(\boldsymbol{r})) P(\boldsymbol{r})\,d\boldsymbol{r} = \frac{1}{f_\max - f_\min}\frac{1}{V}\rho(f) \int_\Omega \delta(f - F(\boldsymbol{r})) \, d\boldsymbol{r} = \frac{1}{f_\max - f_\min} = \text{constant}</math>

which motivates the name "flat histogram". I.e. all costs are equally sampled, and thus there are no barriers. For systems in contact with a heat bath, there is another important advantage: because the sampling is independent of the temperature, one simulation is enough to study all temperatures (thus the name "multicanonic": several temperatures).

{| class="wikitable1" width=300px
|-
|
{{show
|example:
|
On the ferromagnetic [[Ising model]] with ''N'' sites (exemplified on previous section), the density of states can be analytically computed. In this case, a multicanonic ensemble can be used to compute any other quantity ''Q'' by sampling the system according to <math>P(\boldsymbol{r})</math> and using the proper estimator <math>\overline{Q}</math> defined on the previous section.

}}
|}

==Tunneling time and critical slowing down==

Like in any other Monte Carlo method, there are correlations of the samples being drawn from <math>P(\boldsymbol{r})</math>. A typical measurement of the correlation is the ''tunneling time''. The tunneling time is defined by the number of Markov steps (of the Markov chain) the simulation needs to perform a round-trip between the minimum and maximum of the spectrum of ''F''. One motivation to use the tunneling time is that when it crosses the spectra, it passes through the region of the maximum of the density of states, thus de-correlating the process. On the other hand using round-trips ensures that the system visits all the spectrum.

Because the histogram is flat on the variable ''F'', a multicanonic ensemble can be seen as a diffusion process (i.e. a random walk) on the one dimensional line of ''F'' values. [[Detailed balance]] of the process dictates that there is no [[stochastic drift|drift]] on the process.<ref name=Robert/> This implies that the tunneling time, in local dynamics, should scale as a diffusion process, and thus the tunneling time should scale quadratically with the size of the spectrum, ''N'':

:<math>\tau_{tt} \propto N^2</math>

However, in some systems (the Ising model being the most paradigmatic), the scaling suffers from critical slowing down: it is <math>N^{2+z}</math> where <math>z>0</math> depends on the particular system.<ref name=Dayal/>

Non-local dynamics were developed to improve the scaling to a quadratic scaling<ref name=Wolff/> (see [[wolff algorithm]]), beating the critical slowing down. However, it is still an open question whether there is a local dynamics that does not suffer from critical slowing down in spin systems like Ising model.

==References==

{{Reflist|
refs=
<ref name=Berg>{{Cite doi|10.1103/PhysRevLett.68.9}}</ref>
<ref name=Landau>{{Cite doi|10.1103/PhysRevLett.86.2050}}</ref>
<ref name=Newmann>{{cite book
|title= Monte Carlo Methods in Statistical Physics
|last1=Newmann |first1=M E J
|last2=Barkema |first2=G T
|year= 2002
|publisher= Oxford University Press
|location=USA
|isbn=0198517971
|url=
|page=
|pages=
|ref=
}}
</ref>
<ref name=Lee>{{Cite doi|10.1103/PhysRevE.50.R651}}</ref>
<ref name=Metropolis>{{Cite doi|10.1063/1.1699114}}</ref>
<ref name=Robert>{{cite book
|title=Monte Carlo statistical methods
|last1=Robert |first1=Christian |authorlink1=
|last2=Casella |first2=George |authorlink2=
|year= 2004
|publisher=Springer
|location=
|isbn=978-0-387-21239-5
|url=http://www.springer.com/statistics/statistical+theory+and+methods/book/978-0-387-21239-5
}}
</ref>
<ref name=Dayal>{{Cite doi|10.1103/PhysRevLett.92.097201}}</ref>
<ref name=Wolff>{{Cite doi|10.1103/PhysRevLett.62.361}}</ref>
}}

[[Category:Monte Carlo methods]]
[[Category:Statistical algorithms]]
[[Category:Computational physics]]

formulasearchengine - User contributions [en]

Kubo formula