Distance between configurations in Markov chain Monte Carlo simulations

TL;DR

This paper introduces a universal distance measure between configurations in Markov chain Monte Carlo algorithms, demonstrating its properties and effects, especially in multimodal distributions, and revealing geometric insights.

Contribution

It defines a universal distance for local move algorithms in MCMC, explicitly calculates it for Langevin, and explores its implications for multimodal distributions and geometry.

Findings

Distance reduces significantly with tempering in multimodal distributions.

The distance exhibits expected properties for local move algorithms.

An anti-de Sitter-like geometry emerges in highly multimodal cases.

Abstract

For a given Markov chain Monte Carlo algorithm we introduce a distance between two configurations that quantifies the difficulty of transition from one configuration to the other configuration. We argue that the distance takes a universal form for the class of algorithms which generate local moves in the configuration space. We explicitly calculate the distance for the Langevin algorithm, and show that it certainly has desired and expected properties as distance. We further show that the distance for a multimodal distribution gets dramatically reduced from a large value by the introduction of a tempering method. We also argue that, when the original distribution is highly multimodal with large number of degenerate vacua, an anti-de Sitter-like geometry naturally emerges in the extended configuration space.

Equations

$$P(x|y)e^{-S(y)}=P(y|x)e^{-S(x)},$$

$$\displaystyle\hat{T}\equiv e^{S(\hat{x})/2}\,\hat{P}\,e^{-S(\hat{x})/2}=e^{S(\hat{x})/2}\,\bigl{(}\hat{P}\,e^{-S(\hat{x})}\bigr{)}\,e^{S(\hat{x})/2}$$

$$P_n(x_1|x_2)=K_n(x_1,x_2)e^{-(1/2)S(x_1)+(1/2)S(x_2)}$$

$$K_n(x_1,x_2)\equiv ⟨x_1|{\hat{T}}^n|x_2⟩.$$

$$\hat{T}=\sum _{k\ge 0}{\lambda }_k|k⟩⟨k|=|0⟩⟨0|+\sum _{k\ge 1}{\lambda }_k|k⟩⟨k|$$

$$K_n(x_1,x_2)=\sum _{k\ge 0}{c}_k(x_1,x_2){\lambda }_k^n$$

$$c_k(x_1,x_2)\equiv ⟨x_1|k⟩⟨k|x_2⟩.$$

$$f_n(x_1,x_2)\equiv \frac{|{\mathcal{X}}_n(x_1,x_2)|}{|{\mathcal{X}}_n|},$$

$$f_n(x_1,x_2)=P_n(x_1|x_2)\frac{1}{Z}{e}^{-S(x_2)}=K_n(x_1,x_2)\frac{1}{Z}{e}^{-(1/2)S(x_1)-(1/2)S(x_2)}.$$

$$F_n(x_1,x_2)$$

$$F_n(x_1,x_2)=\frac{K_n(x_1,x_2)}{\sqrt{K_n(x_1,x_1)K_n(x_2,x_2)}}=\frac{⟨x_1,n/2|x_2,n/2⟩}{|||x_1,n/2⟩|||||x_2,n/2⟩||}.$$

$$\bullet F_n(x_1,x_2)=F_n(x_2,x_1),$$

$$\bullet 0\le F_n(x_1,x_2)\le 1,$$

$$\bullet F_n(x,x)=1,$$

$$\bullet \underset{n\to \infty }{\lim }{F}_n(x_1,x_2)=1.$$

$$\displaystyle\bullet~{}\mbox{$F_{n}(x_{1},x_{2})\simeq 1$ when $x_{1}$ can be easily reached from $x_{2}$ in $n$ steps.}$$

$$\displaystyle\bullet~{}\mbox{$F_{n}(x_{1},x_{2})\ll 1$ when $x_{1}$ and $x_{2}$ are separated by high potential barriers.}$$

$${\theta }_n(x_1,x_2)\equiv \arccos (F_n(x_1,x_2)).$$

$$\bullet {\theta }_n(x_1,x_2)\ge 0,$$

$$\bullet x_1=x_2\iff {\theta }_n(x_1,x_2)=0,$$

$$\bullet {\theta }_n(x_1,x_2)={\theta }_n(x_2,x_1),$$

$$\bullet {\theta }_n(x_1,x_2)+{\theta }_n(x_2,x_3)\ge {\theta }_n(x_1,x_3).$$

$$d_n^2(x_1,x_2)$$

$$D_n^2(x_1,x_2)$$

$${\dot{x}}_t$$

$$P_t(x|x_0)=⟨x|e^{-t{\hat{H}}_{\mathrm{FP}}}|x_0⟩$$

$$\displaystyle\hat{H}_{\rm FP}={}-D\,\frac{\partial}{\partial x}\,\Bigl{[}\frac{\partial}{\partial x}+S^{\prime}(\hat{x})\Bigr{]}.$$

$$\hat{T}=e^{-ϵ\hat{H}},$$

$$\hat{H}$$

$$\displaystyle={}\Bigl{[}{}-\frac{\partial}{\partial x}+\frac{1}{2}\,S^{\prime}(\hat{x})\Bigr{]}\Bigl{[}\frac{\partial}{\partial x}+\frac{1}{2}\,S^{\prime}(\hat{x})\Bigr{]}$$

$$=-\frac{{\partial }^2}{\partial x^2}+V(\hat{x})$$

$$V(x)\equiv \frac{1}{4}[S^{\prime }(x){]}^2-\frac{1}{2}{S}^{\prime \prime }(x).$$

$$K_t(x,x_0)=⟨x|e^{-t\hat{H}}|x_0⟩$$

$$F_t(x_1,x_2)$$

$$S(x)=\frac{\omega }{2}{x}^2,$$

Full text

1. 1 Introduction

2. 2 Definition of distance

3. 2.1 Preparation

4. 2.2 Half-time overlap

5. 2.3 Definition of distance

6. 2.4 Universality of distance

7. 3 Examples

8. 3.1 Example 1: Gaussian

9. 3.2 Example 2: perturbation around the Gaussian

10. 3.3 Example 3: double-well potential

11. 4 Distance for the simulated tempering and the emergence of AdS-like geometry

12. 4.1 Distance for the simulated tempering

13. 4.2 AdS-like geometry of the extended configuration space

14. 5 Conclusion and outlook

15. A Large scale structure of the transfer matrix

1 Introduction

In Markov chain Monte Carlo (MCMC) simulations, one often encounters a situation where the equilibrium distribution is multimodal and the computation requires an extraordinarily long computer time to make the system reach global equilibrium. In order to speed up the relaxation to global equilibrium, one usually implements an additional method such as the overrelaxation [1] or the simulated/parallel tempering [2, 3, 4, 5]. However, it is not always possible to find a nice overrelaxation method. Also, the tempering method requires an adjustment of tempering parameter such that each local move has a finite amount of acceptance rate, which is usually done in a manual or adaptive way (and sometimes in an empirical way). Thus, it should be useful if there is a numerical measure that quantifies how much two configurations are separate for a given MCMC algorithm and how fast the system approaches global equilibrium. Such measure then can be used to numerically evaluate how the situation gets improved by the introduction of additional methods, which will make the above adjustment easier.

So far there have been proposed various methods for the above purpose, but they only quantify the separation between two different probability distributions (such as the $L_{1}$ distance) or measure the similarity between two data points for a given observable (such as the autocorrelation), and there has still not been a useful apparatus that directly measures an effective distance between two different configurations. In this paper, for a given MCMC algorithm we introduce a distance between two configurations that quantifies the difficulty of transition from one configuration to the other configuration. We argue that the distance takes a universal form for the class of algorithms which generate local moves in the configuration space. We make explicit calculations of distance for the Langevin algorithm, and show that it has desired and expected properties as distance which quantifies the extent of separation between two configurations. We further show that the distance for a multimodal distribution gets dramatically reduced from a large value by the introduction of a tempering method.

The introduction of such distance opens a way to investigate relaxation processes in an MCMC algorithm in terms of the geometry of the configuration space itself (which should not be confused with the geometry of the set of probability distributions). As an example, we argue that, when the original distribution is highly multimodal with large number of degenerate vacua, our distance can be regarded as a geodesic distance with respect to an anti-de Sitter-like (AdS-like) metric. Such a geometrical viewpoint enables us to adjust the tempering parameter in a purely geometrical way; the adjustment can be carried out by requiring that the resulting geodesic distance be minimized.

This paper is organized as follows. In section 2, we first introduce a positive, symmetric operator (to be called the transfer matrix) from the transition matrix of a given MCMC algorithm, and then define the distance between two configurations, which can be written only with the transfer matrix. We argue that the distance takes a universal form for the class of MCMC algorithms that generate local moves of configuration. The statement is confirmed explicitly in Appendix for the Langevin and Metropolis algorithms by using a simple model. In section 3, in order to exemplify that the distance actually has desired and expected properties, we study the distance for the Langevin algorithm in a one-dimensional configuration space with both unimodal and multimodal distributions. In section 4, we study the distance for a multimodal distribution with the implementation of the simulated tempering method, and show that the distance certainly receives a huge amount of reduction from a large value. We also investigate the local geometry of the extended configuration space of the simulated tempering, and show that it has an AdS-like geometry. Section 5 is devoted to conclusion and outlook for future work.

2 Definition of distance

In this section, we define the distance between two configurations for a given MCMC algorithm. The distance will be written only with the kernel of a positive, symmetric matrix. We argue that the distance should take a universal form for the class of MCMC algorithms that generate local moves in the configuration space.

2.1 Preparation

Let $\mathcal{M}=\{x\}$ be a configuration space, and suppose that we are given an MCMC algorithm which generates a new configuration $x$ from a configuration $y$ with the conditional probability $P(x|y)$ at each step. We assume that this yields a stochastic process which has suitable ergodic properties such that $P_{n}(x|x_{0})\equiv(P^{n})(x|x_{0})$ converges to a unique equilibrium distribution $p_{\rm eq}(x)$ in the limit $n\to\infty$, irrespectively of the initial value $x_{0}$. We further make two assumptions:

The algorithm satisfies the detailed balance condition for a given real-valued action $S(x)$,

[TABLE]

which ensures that the equilibrium distribution is given by $p_{\rm eq}(x)=e^{-S(x)}/Z$ ($Z=\int dx\,e^{-S(x)}$). 2. 2.

The eigenvalues of $P$ are all positive.111Note that the second condition is not too restrictive. In fact, when a transition matrix does not satisfy this condition (i.e., when some of the eigenvalues are negative), one can consider a new transition matrix $P_{\rm new}\equiv(P_{\rm old})^{2}$, for which the eigenvalues are all positive and the same equilibrium distribution is reached. Note also that this condition is always satisfied for the Langevin algorithm [see, e.g., (3.2)–(3.6) below].

By using the bra-ket notation $P(x|y)=\langle x|\hat{P}|y\rangle$ and the configuration operator $\hat{x}\equiv\int dx\,x\,|x\rangle\langle x|$, the above two conditions can be rephrased as a single statement that the operator $\hat{P}\,e^{-S(\hat{x})}$ is positive and symmetric. This also means that the “transfer matrix”

[TABLE]

is positive and symmetric. $\hat{T}$ shares the same set of eigenvalues as $\hat{P}$, and according to our assumptions, all the eigenvalues are positive and the largest eigenvalue is unity with no degeneracy. Note that $P_{n}(x_{1}|x_{2})=\langle x_{1}|\hat{P}^{n}|x_{2}\rangle$ can be expressed in the form

[TABLE]

with the kernel of $\hat{T}$,

[TABLE]

If we introduce the spectral decomposition of $\hat{T}$ as222 The ground state wave function is given by $\langle x|0\rangle=e^{-(1/2)S(x)}/\sqrt{Z}$.

[TABLE]

with $\lambda_{0}=1>\lambda_{1}\geq\lambda_{2}\geq\cdots>0$, then the kernel is expressed as

[TABLE]

with

[TABLE]

The relaxation of the system to global equilibrium in the stochastic process corresponds to the decoupling of higher modes from $\hat{T}^{n}=|0\rangle\langle 0|+\sum_{k\geq 1}\lambda_{k}^{n}\,|k\rangle\langle k|$ as $n$ increases. Thus, the large scale behavior of relaxation333 By “large scale behavior (or structure)” we mean the behavior (or structure) both at large step numbers and at scales larger than the increment of configuration at each Monte Carlo step. is determined by the wave functions $\langle x|k\rangle$ with small eigenvalues $\lambda_{k}$. Note that the decoupling occurs earlier for modes with larger $k$.

If we further introduce the Hamiltonian $\hat{H}$ by using a small time increment $\epsilon$ as $\hat{T}=e^{{}-\epsilon\hat{H}}$, then $\hat{T}^{n}$ is expressed as $e^{{}-t\hat{H}}$ with $t\equiv n\epsilon$. For generic MCMC algorithms, $\hat{H}$ is a nonlocal operator acting on functions over $\mathcal{M}$. However, for the class of MCMC algorithms that generate only local moves of configuration, $\hat{H}$ should become local in the limit $\epsilon\to 0$ (see subsection 2.4 for further arguments on locality).

We would like to introduce a distance $d_{n}(x_{1}|x_{2})$ for a pair of configurations, $x_{1},\,x_{2}\in\mathcal{M}$, such that it enjoys the following properties in addition to the usual axioms of distance:

• •

(P1)  The distance $d_{n}(x_{1}|x_{2})$ vanishes in the limit $n\to\infty$ for any pairs $x_{1}$ and $x_{2}$, in order to reflect the fact that any configuration can be reached from every configuration in finite steps.

• •

(P2)  If $x_{1}$ can be easily reached from $x_{2}$, then the distance is small even for finite $n$.

• •

(P3)  If the distribution $e^{-S(x)}/Z$ is multimodal, and if $x_{1}$ and $x_{2}$ belong to different modes, then the distance is very large for finite $n$.

As we see below, such a distance can be naturally introduced if we look at transitions in $\mathcal{M}$ that is already in global equilibrium.

2.2 Half-time overlap

Let the system be in global equilibrium with probability distribution $p_{\rm eq}(x)=e^{-S(x)}/Z$. We denote by $\mathcal{X}_{n}$ the set of sequences of $n$ processes in $\mathcal{M}$ and by $\mathcal{X}_{n}(x_{1},x_{2})$ the subset of $\mathcal{X}_{n}$ that consists of sequences which start from $x_{2}$ and end at $x_{1}$. We then define $f_{n}(x_{1},x_{2})$ to be the ratio of the sizes of two sets:

[TABLE]

which can be expressed as

[TABLE]

The latter expression shows that $f_{n}(x_{1},x_{2})$ is a symmetric function of $x_{1}$ and $x_{2}$. $f_{n}(x_{1},x_{2})$ gives the probability to find a sequence of $n$ processes from $x_{2}$ to $x_{1}$ (or from $x_{1}$ to $x_{2}$) out of all the possible $n$ processes, and thus expresses “mobility” between two configurations $x_{1}$ and $x_{2}$. Note that $f_{n}(x_{1},x_{2})$ should take a very small value if the equilibrium distribution is multimodal and $x_{1}$ and $x_{2}$ belong to different modes.

Since our interest is in the mobility between two different configurations, we normalize the mobility such that it takes the maximal value ($=1$) for closed loops which start from and end at the same configuration:

[TABLE]

We will call the normalized mobility $F_{n}(x_{1},x_{2})$ the half-time overlap of configurations $x_{1}$ and $x_{2}$ for $n$ steps. In fact, by using (2.9) and introducing the “half-time elapsed state” $|x,n/2\rangle\equiv\hat{T}^{n/2}\,|x\rangle$, the function $F_{n}(x_{1},x_{2})$ can actually be expressed as the overlap of two normalized half-time elapsed states:

[TABLE]

One can easily prove that $F_{n}(x_{1},x_{2})$ enjoys the following properties:

[TABLE]

Furthermore, from the properties of $f_{n}(x_{1},x_{2})$, we expect that

[TABLE]

2.3 Definition of distance

Based on the half-time overlap given above, we define a distance between $x_{1},\,x_{2}\in\mathcal{M}$ as follows:

[TABLE]

This should satisfy the properties (P1)–(P3) due to (2.15)–(2.17), as well as the following axioms of distance:444The axiom of nondegeneracy (i.e., $\theta_{n}(x_{1},x_{2})=0~{}\Rightarrow~{}x_{1}=x_{2}$) holds for all the algorithms that involve relaxations. Note that the algorithm consisting of the global updates based on the global heat bath is not in this class, because the system transits to the equilibrium state at the first Monte Carlo step (the transfer matrix is given by the projection, $T=|0\rangle\langle 0|$). In fact, for this case, the distance $\theta_{n}(x_{1},x_{2})$ vanishes for any pairs $x_{1}$ and $x_{2}$, which is also consistent with our definition of distance.

[TABLE]

It is often convenient to introduce other distances from the half-time overlap:555Note the similarity of the definition of distance between configurations to that between states in quantum information. There, $F_{n}(x_{1},x_{2})$ corresponds to the fidelity of two pure states $\rho_{1,2}=|x_{1,2}\rangle\langle x_{1,2}|/||\,|x_{1,2}\rangle\,||^{2}$, and the distances $\theta_{n}(x_{1},x_{2})$ and $D_{n}(x_{1},x_{2})$ to Bures length and Bures distance, respectively.

[TABLE]

Although these distances do not satisfy the triangle inequality, they agree with $\theta_{n}(x_{1},x_{2})$ up to higher order corrections when $\theta_{n}(x_{1},x_{2})$ is small enough. For the rest of this paper we will mainly use $d_{n}(x_{1},x_{2})$ as the definition of distance, because it gives the simplest expressions for the following examples.

2.4 Universality of distance

We close this section with a comment on the universality of our distance. As long as the chosen algorithm generates only local moves of configuration, we expect that the large scale structure of distance $d_{n}(x_{1},x_{2})$ takes a universal form, in the sense that differences of distance between two such algorithms can always be absorbed into a rescaling of $n$. In fact, our distance is totally expressed with the kernel of the transfer matrix $\hat{T}$, and, when the transition is sufficiently local, the corresponding Hamiltonian $\hat{H}$ is a local operator acting on functions over $\mathcal{M}$ in almost the same way. We thus expect that the eigenvalues $\lambda_{k}$ and the wave functions $\langle x|k\rangle$ are almost the same for small $k$’s, which ensures the same large scale structure of the kernel $K_{n}(x_{1},x_{2})$ and thus that of the distance $d_{n}(x_{1},x_{2})$. This statement is explicitly checked in Appendix for the Langevin and Metropolis algorithms using a simple one-dimensional model. Note that the argument for universality are more trustworthy when the degrees of freedom of system become larger.

This universality will not hold for algorithms that generate nonlocal moves of configuration (such as the overrelaxation algorithm).666 The tempering algorithms are nonlocal when viewed from the original configuration space $\mathcal{M}=\{x\}$, but they actually generate local moves in the extended configuration space, so that the universality of the distance is expected to hold also for these algorithms. See arguments below (4.5) in subsection 4.2 for details. The Hybrid Monte Carlo (HMC) algorithm [6] is marginal in this sense, and we leave it for future study to investigate whether the distance for the HMC algorithm exhibits the same behavior as local algorithms.

3 Examples

Expecting the universality of distance commented in the previous section, we consider the Langevin method as an MCMC algorithm, and write down the explicit form of distance for a configuration space $\mathcal{M}=\mathbb{R}$. We show that the resulting distance actually satisfies the properties (P1)–(P3).

Let $\nu_{t}$ be the Gaussian white noise with diffusion coefficient $D$, $\langle\nu_{t}\,\nu_{t^{\prime}}\rangle_{\nu}=2D\,\delta(t-t^{\prime})$, and $x_{t}=x_{t}(x_{0},[\nu])$ be the solution to the Langevin equation

[TABLE]

Then, the probability distribution777In this section we exclusively use a continuous notation; namely, the distribution $P_{n}(x|x_{0})$ will be written as $P_{t}(x|x_{0})$ with $t=n\epsilon$, where $\epsilon$ is the time increment. $P_{t}(x|x_{0})\equiv\langle\delta(x-x_{t}(x_{0},[\nu])\rangle_{\nu}$ is expressed as

[TABLE]

with the Fokker-Planck Hamiltonian (we will set $D=1$ below):

[TABLE]

The transfer matrix is expressed as

[TABLE]

where the positive, symmetric Hamiltonian $\hat{H}$ is given by

[TABLE]

with

[TABLE]

The corresponding kernel is then given by

[TABLE]

with which the half-time overlap is expressed as

[TABLE]

3.1 Example 1: Gaussian

We first consider the action

[TABLE]

for which the Hamiltonian $\hat{H}$ takes the form

[TABLE]

Note that the last term removes the zero-point energy of this harmonic oscillator. By using the kernel

[TABLE]

the distance $d_{t}(x_{1},x_{2})$ is easily obtained to be

[TABLE]

This gives a flat and translation invariant metric in the entire configuration space $\mathcal{M}=\mathbb{R}$.888Our discussion can be easily generalized to the case $\mathcal{M}=\mathbb{R}^{N}=\{{\bf x}=(x^{i})\}$ and $S({\bf x})=\sum_{i}\omega_{i}(x^{i})^{2}$. The distance is then given by $d_{t}^{2}({\bf x}_{1},{\bf x}_{2})=\sum_{i}\,(\omega_{i}/2\sinh(\omega_{i}t))\,|x_{1}^{i}-x_{2}^{i}|^{2}.$ Note that the distance decreases exponentially as $d_{t}^{2}\sim e^{-\omega t}$, from which we find that the relaxation time of the system is given by $\sim 1/\omega$.999If we take the limit $\omega\to 0$ (corresponding to the pure Brownian motion), the distance is given by $d_{t}^{2}(x_{1},x_{2})=(1/2t)\,|x_{1}-x_{2}|^{2}$, and thus the exponential damping disappears. This reflects the fact that the relaxation time becomes infinite in the limit $\omega\to 0$ in the sense that there is no normalizable equilibrium distribution for the pure Brownian motion. Since the manner of relaxation is almost the same for unimodal distributions, we see that the distance rapidly decreases to zero when the action gives a unimodal distribution. We thus confirm that the properties (P1) and (P2) hold in this example.

3.2 Example 2: perturbation around the Gaussian

We now consider the case where the action (3.9) is perturbed with a quartic term:

[TABLE]

The perturbative calculation of distance can be done easily, and we find to the first order perturbation:

[TABLE]

where $s\equiv\sinh(\omega t)$ and $c\equiv\cosh(\omega t)$. This shows that the distance generically does not take a flat or translation invariant form when the unimodal distribution is not Gaussian.

3.3 Example 3: double-well potential

Finally, in order to confirm the property (P3), we consider the case where a high potential barrier exists between local minima:

[TABLE]

We assume that $\beta$ takes a large value so that the equilibrium distribution $e^{-S(x)}/Z$ is multimodal. The potential $V(x)=(1/4)[S^{\prime}(x)]^{2}-(1/2)S^{\prime\prime}(x)$ becomes sextic (see Fig. 1),

[TABLE]

and has local minima at $x=0$ and $x={}\pm x_{+}$ with x_{+}\equiv\bigl{[}\bigl{(}2+\sqrt{1+9/\beta}\bigr{)}/3\bigr{]}^{1/2}\simeq 1.

One can easily find that $x={}\pm x_{+}$ are global minima by looking at the values of $S^{\prime\prime}(x)$ (note that $S^{\prime}(x)$ vanishes there).

The eigenvalues $E_{k}$ $(E_{0}=0<E_{1}<E_{2}\cdots)$ of $\hat{H}$ can be roughly estimated as follows. We first note that the Gaussian approximation around the global minima should be effective when $\beta\gg 1$, and thus the first two eigenstates of $\hat{H}$ can be approximated as superpositions of the ground states $|0_{+}\rangle$ and $|0_{-}\rangle$ of the approximated Gaussian potential $V(x)\simeq(1/2)\,V^{\prime\prime}(\pm x_{+})\,(x\mp x_{+})^{2}\simeq 4\beta^{2}\,(x\mp 1)^{2}$ around the right ($x=x_{+}\simeq+1$) and left ($x=-x_{+}\simeq-1$) minima:

[TABLE]

The energy difference between two states is exponentially small, $E_{1}=O(e^{-\beta/2})$, as can be estimated by an instanton calculation.

As for the second excitation of $\hat{H}$, we expect the relations $E_{0}=0\lesssim E_{1}=O(e^{-\beta/2})\ll E_{2}=O(\beta)$. In fact, there are two possible approximations for the second excitation; one is to represent the second excitation as a superposition of the first excited states of the approximated Gaussian potentials around the global minima $x={}\pm x_{+}$, and the other is to represent it as the ground state of the approximated Gaussian potential around the local minimum $x=0$. In our case, the latter approximation is applicable, because the former gives $E_{2}\simeq 4\beta$ and the latter gives $E_{2}\simeq 2\beta$ (see Fig. 3 in Appendix).

By using the expansion of the kernel $K_{t}(x_{1},x_{2})=\langle x_{1}|e^{-tH}|x_{2}\rangle$ [see (2.6)],

[TABLE]

the distance $d_{t}(x_{1},x_{2})$ is expressed as

[TABLE]

We now consider the following two cases:

[TABLE]

We first discuss Case 1 where $x_{1}$ and $x_{2}$ belong to the same mode. If we expand $d_{t}^{2}(x_{1},x_{2}=x_{1}+\delta x)$ to the second order of $\delta x$, the coefficient of $e^{-kE_{1}t}$ is given by

[TABLE]

which is vanishingly small as can be understood by considering the supports of the wave functions $\langle x|0_{\pm}\rangle$. Thus, the dominant contribution to the distance comes from the second term in (3.20), and the two configurations $x_{1}$ and $x_{2}$ can be reached from each other with a rather short time $\sim 1/E_{2}=O(1/\beta)$. This confirms that two configurations are close to each other even for a multimodal distribution if they belong to the same mode. In contrast, as for Case 2 where $x_{1}$ and $x_{2}$ belong to different modes, the coefficient of $e^{-E_{1}t}$ is not small, so that the dominant contribution to the distance comes from the first excited state, and the transition between $x_{1}$ and $x_{2}$ requires an exponentially long time $\sim 1/E_{1}=O(e^{\beta/2})$. Accordingly, the distance $d_{t}(x_{1},x_{2})$ between $x_{1}\simeq+1$ and $x_{2}\simeq-1$ has a large value $d_{t}^{2}(x_{1},x_{2})\propto\beta$,101010$F_{t}(x_{1},x_{2})=e^{-(1/2)\,d_{t}^{2}(x_{1},x_{2})}$ for large $\beta$ can be estimated to be $e^{-{\rm const.}\beta}$ by using an instanton analysis. which decreases only very slowly as $t$ elapses. We thus see that the property (P3) certainly holds in this example.

4 Distance for the simulated tempering

and the emergence of AdS-like geometry

In this section, we show that the value of distance $d_{n}(x_{1},x_{2})$ for a multimodal probability distribution gets dramatically reduced when one introduces a tempering algorithm. We further argue that the effective distance defined for the extended configuration space gives an AdS-like geometry when the original distribution is highly multimodal with large number of degenerate vacua.

4.1 Distance for the simulated tempering

Suppose that we are considering a multimodal distribution, and that we implement a simulated tempering method [2]. We take as the tempering parameter the overall coefficient $\beta$ of the action and introduce the parameter set $\mathcal{A}\equiv\{\beta_{a}\}_{a=0,1,\ldots,A}$ such that $\beta_{0}>\beta_{1}>\cdots>\beta_{A}$, where $\beta_{0}$ is the overall coefficient of the original action. We denote by $S(x;\beta_{a})$ the action with $\beta_{0}$ replaced by $\beta_{a}$.

In the simulated tempering algorithm, we extend the original configuration space $\mathcal{M}$ to $\mathcal{M}\times\mathcal{A}=\{X=(x,\beta_{a})\}$, and introduce a stochastic process such that it converges to global equilibrium with the probability distribution

[TABLE]

Ideally the weights $w_{a}$ $(a=0,1,\ldots,A)$ are chosen such that the appearance ratio of the $a$-th configuration is the same for all $a$ [i.e. $\int dx\,P_{\rm eq}(x,\beta_{a})=1/(A+1)$].111111Such weights are given by the inverse of the $a$-th partition functions, $w_{a}\propto 1/Z(\beta_{a})\equiv[\int dx\,e^{-S(x;\beta_{a})}]^{-1}$. Since it is difficult to evaluate $Z(\beta_{a})$ numerically, $w_{a}$ are usually determined in a manual or adaptive way by investigating the appearance ratio. In this paper, we assume that $w_{a}$ are already chosen in this way. A possible stochastic process that converges to $P_{\rm eq}(X)$ is given by a Markov chain which consists of the following two steps:

• •

Step 1. Generate a transition in the $x$ direction, $X=(x,\beta_{a})\to X^{\prime}=(x^{\prime},\beta_{a})$, with some proper algorithm (such as the Langevin or Metropolis algorithm).

• •

Step 2. Generate a transition in the $\beta$ direction, $X=(x,\beta_{a})\to X^{\prime}=(x,\beta_{a^{\prime}=a\pm 1})$, with the probability

[TABLE]

It is easy to see that each process satisfies the detailed balance condition with respect to $P_{\rm eq}(X)$. After one obtains a sample from the extended configuration space with probability distribution $P_{\rm eq}(X)$, one estimates the expectation values with respect to the original action $S(x;\beta_{0})$ by using only a subsample with $\beta_{a=0}$.

Since the distribution with smaller $\beta_{a}$ is less multimodal, two configurations belonging to different modes at $\beta_{0}$ can now be easily reached from each other by moving around in the extended configuration space, as long as moves in the $\beta$ direction occur frequently (as we assume here). This improves the convergence to global equilibrium, and accordingly, the distance between two configurations $X_{1}=(x_{1},\beta_{0})$ and $X_{2}=(x_{2},\beta_{0})$ will be reduced. To demonstrate that this actually happens, we calculated the distance between $X_{1}=(x_{1}=+1,\beta_{0})$ and $X_{2}=(x_{2}=-1,\beta_{0})$ for the action $S(x;\beta_{0})=(\beta_{0}/2)\,(x^{2}-1)^{2}$ with $\beta_{0}=20$, by using the simplest setting $A=1$ and $\beta_{1}=1$. As for Step 1 above, we adopted the Metropolis algorithm using Gaussian proposal distribution with variance $\sigma^{2}=0.01$, where the configuration space is restricted to the interval $[-3,3]$ and is latticized with cutoff $a=0.001$. Furthermore, by denoting the transition matrix for Step $s$ by $\hat{P}_{(s)}$ $(s=1,2)$, we set the transition matrix $\hat{P}$ in the simulated tempering algorithm to $\hat{P}=\hat{P}_{(1)}\hat{P}_{(2)}\hat{P}_{(1)}$ so that the combined transition matrix satisfies the detailed balance condition. Below is the result we obtained:121212 In order to align the scale of $n$, the transition matrix $\hat{P}=\hat{P}_{(1)}^{\,2}$ is used for the original local algorithm without tempering (central column in the table).

[TABLE]

We clearly see that the introduction of the simulated tempering method dramatically reduces the distance for such a multimodal distribution.

4.2 AdS-like geometry of the extended configuration space

When a system has very large degrees of freedom (as in the case of the large-volume or low-temperature limit), the issue related to the multimodality of probability distribution becomes serious. However, it will then become a good approximation to coarse-grain the elements of configuration space $\mathcal{M}$ by regarding configurations in the same mode as a single configuration, and to introduce a new configuration space $\bar{\mathcal{M}}$ that consists of coarse-grained configurations which are separate from each other. We assume that the way of separation between two neighboring configurations is uniform over $\bar{\mathcal{M}}$, keeping in mind a situation where the original distribution has highly degenerate vacua.

Suppose that we introduce the simulated tempering algorithm to such system. The original configuration space $\mathcal{M}$ is then extended to $\mathcal{M}\times\mathcal{A}$, as discussed in the previous subsection. When the degrees of freedom are very large, the spacing in the parameter set $\mathcal{A}=\{\beta_{a}\}$ must be sufficiently small such that two adjacent configurations $(x,\beta_{a})$ and $(x,\beta_{a+1})$ has an acceptance rate of $O(1)$. Then, we may (and we will) regard $\mathcal{A}$ as a continuous set, $\mathcal{A}=\{\beta\,|\,\beta_{A}\leq\beta\leq\beta_{0}\}$.

Now let us consider the distance in the extended, coarse-grained configuration space: $\bar{\mathcal{M}}\times\mathcal{A}=\{X=(x,\beta)\}$. We define the metric $ds^{2}=g_{MN}^{(n)}(X)\,dX^{M}dX^{N}$ [$(X^{M})=(x,\beta)$] such that its geodesic distance between two arbitrarily chosen points $X_{1}=(x_{1},\beta_{1})$ and $X_{2}=(x_{2},\beta_{2})$ gives our distance $d_{n}(X_{1},X_{2})$.

Such metric will take the following form due to the uniformity over $\bar{\mathcal{M}}$:131313We have assumed that the coefficient of $d\beta\,dx$ is relatively small and can be neglected.

[TABLE]

The functions $f(\beta)$ and $g(\beta)$ will be studied in detail in the subsequent paper [8], and we here make a simple analysis of their properties. Since two different configurations $(x_{1},\beta)$ and $(x_{2},\beta)$ in $\bar{\mathcal{M}}\times\mathcal{A}$ belong to two different modes of the distribution $e^{-S(x;\beta)}/Z(\beta)$, the transition between them becomes more difficult as $\beta$ increases. Thus, $g(\beta)$ should be an increasing function at least when $\beta$ is large, and we assume that it takes an asymptotic form $g(\beta)\propto\beta^{q}$ for large $\beta$ with some positive number $q$.141414Without a tempering method, the distance between two configurations belonging to different modes for fixed large $\beta$ can be roughly estimated by an instanton analysis to be proportional to $\beta$ [see discussions below (3.21)]. This implies that the function $g(\beta)$ increases more slowly than $\beta$ when the simulated tempering is implemented, because configurations at larger $\beta$ will get more benefit from the tempering method. We thus expect that $q$ satisfies the inequality $0<q\leq 1$. Our on-going work with numerical calculation [8] shows that $q$ is actually in this range, excluding a logarithmic growth of $g(\beta)$. On the other hand, since the transition in the $\beta$ direction has no obstacle, we expect that $f(\beta)$ has a simple behavior. Note that, when the measure in the $\beta$ direction is scale invariant [i.e., when $f(\beta)\propto 1/\beta^{2}$] for large $\beta$, the above metric has the following asymptotic form:

[TABLE]

which is nothing but the Euclidean AdS metric, as can be easily seen by the redefinition of variable, $z\equiv{\rm const.}\,\beta^{-q/2}$:

[TABLE]

We here argue that the appearance of AdS-like geometry [eq. (4.3) with $g(\beta)\sim\beta^{q}$] should be universal. In fact, as discussed above, when implementing the tempering method to an MCMC simulation for a large system, the spacing in the parameter set $\mathcal{A}=\{\beta_{a}\}$ must be sufficiently small, which allows us to regard $\mathcal{A}$ as a continuous set. Then, the stochastic processes in the $x$ and $\beta$ directions given in subsection 4.1 can be regarded as local moves in the extended configuration space $\mathcal{M}\times\mathcal{A}$, which in turn define local moves in the extended, coarse-grained configuration space $\bar{\mathcal{M}}\times\mathcal{A}$. Thus, combining with arguments in subsection 2.4, we expect that the distance takes a universal form for $\bar{\mathcal{M}}\times\mathcal{A}$ in the sense that it does not depend on the particular MCMC algorithms that generate local moves along the $x$ and $\beta$ directions.151515However, the distance should depend on the way of preparing the parameter set $\{\beta_{a}\}$ even though the set $\mathcal{A}$ can be regarded as almost continuous. This is why we expect that the emergence of AdS-like geometry is universal.

5 Conclusion and outlook

In this paper, we have defined the distance $d_{n}(x_{1},x_{2})$ between two configurations $x_{1}$ and $x_{2}$ for a given MCMC algorithm. The distance is written only with the kernel of the transfer matrix $\hat{T}$, and various properties of the distance (including the universality for local MCMC algorithms) are easily understood as the reflection of properties of the transfer matrix.

We made a detailed study of the distance for multimodal distributions, and showed that distances between two different modes get dramatically reduced by the introduction of the simulated tempering method. We also considered the effective distance in the extended configuration space by regarding configurations in the same mode as a single configuration, and argued that the distance between two configurations may then be regarded as the geodesic distance in the extended configuration space with respect to an AdS-like metric. It will be demonstrated in the subsequent paper [8] that this is actually the case.

The introduction of distance between configurations opens a way to investigate relaxation processes in an MCMC algorithm in terms of the geometry of the configuration space itself. Among possible applications of the present formalism, one interesting application is to determine the parameter set $\{\beta_{a}\}$ in the simulated tempering method by requiring that it minimize the geodesic distance in the bulk with given ends on the boundary at $\beta=\beta_{0}$. Furthermore, since the coefficient function $f(\beta)$ in (4.3) has a dependence on $\{\beta_{a}\}$, it should be interesting to investigate whether the optimized parameter set gives an AdS geometry. This point will be argued positively in [8].

As discussed in subsection 2.4, we expect that our distance takes a universal form for local MCMC algorithms. It should be interesting to study to which extent this universality holds. In particular, it is important to investigate whether the HMC algorithm has the distance similar to that for local algorithms.

In this paper, we have discussed only the case where the action $S(x)$ is real. When the action takes complex values as in QCD at finite density (see [7] for a review), we cannot directly use the present definition of distance because there can be no real-valued transfer matrix for such complex Boltzmann weight. The reweighting method gives a real-valued transfer matrix, but the reweighting method is generically inapplicable because of the sign problem. There are various approaches to the sign problem. One approach which is currently under intensive study is the complex Langevin method [9] (see also [10, 11, 12, 13, 14, 15]), and another is the Lefschetz thimble method [16] (see also [17, 18, 19, 20, 21, 22, 23, 24]).161616 Recently a new interesting method has been proposed that uses a complex path optimized with respect to a cost function [25]. It must be important to investigate whether nice distances can be introduced to these algorithms. As for the complex Langevin method, in particular, it will be very interesting if the condition [14] to be free from the wrong convergence problems [12, 13, 14, 15] can be rephrased in terms of the distance. We expect that the formalism of [26] developed for a complex Hamiltonian will be useful. As for the Lefschetz thimble algorithm, it must be interesting to investigate the geometry of the tempering algorithm that was recently introduced for integration over the Lefschetz thimbles, where the tempering parameter is set to the flow time of the antiholomorphic gradient flow [22, 23].

A study along these lines is now in progress and will be reported elsewhere.

Acknowledgments

The authors thank Yuho Sakatani and Sotaro Sugishita for useful discussions. This work was partially supported by JSPS KAKENHI (Grant Numbers 16K05321 and JP17J08709).

Appendix A Large scale structure of the transfer matrix

In this Appendix, we explicitly evaluate the transfer matrix for the Langevin and Metropolis algorithms with a one-dimensional configuration space $\mathcal{M}=\mathbb{R}$, and show that they have the same large scale structure.

As for the Langevin algorithm, if we set the infinitesimal time increment to be $\epsilon$, the transfer matrix can be written as [see (3.4)–(3.6)]

[TABLE]

with $V(x)=(1/4)[S^{\prime}(x)]^{2}-(1/2)S^{\prime\prime}(x)$. As for the Metropolis algorithm, if we use a symmetric Gaussian proposal distribution with variance $\sigma^{2}$, the off-diagonal elements of the transfer matrix are given by171717The diagonal elements can be read off from the condition that $e^{-(1/2)S(x)}$ be the eigenstate of $\hat{T}$ with unit eigenvalue.

[TABLE]

We thus see that, under the identification $\sigma^{2}\sim\epsilon$, the Hamiltonians $\hat{H}={}-(1/\epsilon)\,\ln\hat{T}$ obtained from (A.1) and (A.2) are both local in the limit $\epsilon\to 0$ and have the same tendency to enhance the values of matrix elements when $|x-y|$ and $|S(x)-S(y)|$ are small.

We numerically diagonalize the transfer matrix for the action $S(x)=(\beta/2)\,(x^{2}-1)^{2}$ with $\beta=20$ by latticizing the configuration space. Restricting the configuration space to the interval $[-2,\,2]$, we set the spatial cutoff $a$, the time increment $\epsilon$ and the variance $\sigma^{2}$ of the proposal distribution in the Metropolis to $a=0.005$, $\epsilon=a^{2}$ and $\sigma^{2}=2a^{2}$, respectively. Below are the obtained results of eigenvalues:

[TABLE]

We see that the eigenvalues for the two algorithms can be coincided by rescaling the unit. The wave functions also have the same forms as depicted in Fig. 3.