Non-Reversible Parallel Tempering: a Scalable Highly Parallel MCMC Scheme

TL;DR

This paper introduces a non-reversible parallel tempering method for MCMC that outperforms traditional reversible schemes, providing theoretical insights, optimal scheduling, and practical algorithms for sampling complex distributions.

Contribution

It formalizes the distinction between reversible and non-reversible PT schemes, demonstrating the superiority of non-reversible methods and developing an optimal annealing schedule.

Findings

Non-reversible PT dominates reversible PT in performance.

Scaling limits differ: non-reversible is piecewise-deterministic, reversible is diffusive.

Numerical examples validate theoretical and methodological advances.

Abstract

Parallel tempering (PT) methods are a popular class of Markov chain Monte Carlo schemes used to sample complex high-dimensional probability distributions. They rely on a collection of $N$ interacting auxiliary chains targeting tempered versions of the target distribution to improve the exploration of the state-space. We provide here a new perspective on these highly parallel algorithms and their tuning by identifying and formalizing a sharp divide in the behaviour and performance of reversible versus non-reversible PT schemes. We show theoretically and empirically that a class of non-reversible PT methods dominates its reversible counterparts and identify distinct scaling limits for the non-reversible and reversible schemes, the former being a piecewise-deterministic Markov process and the latter a diffusion. These results are exploited to identify the optimal annealing schedule for non-reversible PT and to develop an iterative scheme approximating this schedule. We provide a wide range of numerical examples supporting our theoretical and methodological contributions. The proposed methodology is applicable to sample from a distribution $\pi$ with a density $L$ with respect to a reference distribution $\pi_0$ and compute the normalizing constant. A typical use case is when $\pi_0$ is a prior distribution, $L$ a likelihood function and $\pi$ the corresponding posterior.

Figures

Equations

$${\pi }^{(\beta )}(x)=\frac{L(x{)}^{\beta }{\pi }_0(x)}{\mathcal{Z}(\beta )}=\frac{e^{-\beta V(x)-V_0(x)}}{\mathcal{Z}(\beta )},$$

$${\text{K}}_n^{\text{PT}}(\bar{\text{x}},A)={\text{K}}_n^{\text{comm}}{\text{K}}^{\text{expl}}(\bar{\text{x}},A):=\int {\text{K}}_n^{\text{comm}}(\bar{\text{x}},\text{d}{\bar{\text{x}}}^{\prime }){\text{K}}^{\text{expl}}({\bar{\text{x}}}^{\prime },A).$$

$${\text{K}}^{\text{expl}}((\text{x},\text{i}),A_0\times A_1\times \dots A_N\times \{{\text{i}}^{\prime }\})=\prod _{i=0}^N{K}^{({\beta }_i)}(x^i,A_i){\delta }_{\text{i}}({\text{i}}^{\prime }),$$

$${\text{x}}^{(i,i+1)}$$

$${\text{K}}^{(i,i+1)}(\bar{\text{x}},\cdot )=(1-{\alpha }^{(i,i+1)}(\bar{\text{x}})){\delta }_{\bar{\text{x}}}(\cdot )+{\alpha }^{(i,i+1)}(\bar{\text{x}}){\delta }_{{\bar{\text{x}}}^{(i,i+1)}}(\cdot ).$$

$${\alpha }^{(i,i+1)}(\bar{\text{x}})$$

$${\text{K}}^{\text{even}}=\prod _{i\text{ even}}{\text{K}}^{(i,i+1)},{\text{K}}^{\text{odd}}=\prod _{i\text{ odd}}{\text{K}}^{(i,i+1)},$$

$${\text{K}}^{\text{SEO}}=\frac{1}{2}{\text{K}}^{\text{even}}+\frac{1}{2}{\text{K}}^{\text{odd}},{\text{K}}_n^{\text{DEO}}=\{\begin{array}{ll}{\text{K}}^{\text{even}}& \text{if n is even,}\\ {\text{K}}^{\text{odd}}& \text{if n is odd.}\end{array}$$

$$I_{n+1}^j$$

$$T_{\uparrow ,k}^j$$

$$T_{\downarrow ,k}^j$$

$$s^{(i,i+1)}=\mathbb{E}\left[{\alpha }^{(i,i+1)}(\bar{\text{X}})\right]=\mathbb{E}\left[\exp \left(\min \left\{0,\left({\beta }_{i+1}-{\beta }_i\right)\left(V^{(i+1)}-V^{(i)}\right)\right\}\right)\right],$$

$$K^{\mathrm{DEO}}((i,\epsilon ),(i^{\prime },{\epsilon }^{\prime }))=K^{\mathrm{DEO}}((i^{\prime },-{\epsilon }^{\prime }),(i,-\epsilon )),$$

$$\tau =\sum _{j=0}^N\underset{n\to \infty }{\lim }\frac{\mathbb{E}[{\mathcal{R}}_n^j]}{n}=\frac{N+1}{\mathbb{E}[T]}.$$

$${\mathbb{E}}_{\mathrm{SEO}}[T]$$

$${\mathbb{E}}_{\mathrm{DEO}}[T]$$

$${\tau }_{\mathrm{SEO}}({\mathcal{P}}_N)$$

$${\tau }_{\mathrm{DEO}}({\mathcal{P}}_N)$$

$$s(\beta ,{\beta }^{\prime })$$

$$r(\beta ,{\beta }^{\prime })$$

$$\lambda (\beta )=\underset{\delta \to 0}{\lim }\frac{r(\beta ,\beta +\delta )}{|\delta |}.$$

$$\lambda (\beta )=\frac{1}{2}\mathbb{E}\left[|V_1^{(\beta )}-V_2^{(\beta )}|\right],$$

$$r(\beta ,{\beta }^{\prime })$$

$$\sum _{i=1}^{N}r({\beta }_{i-1},{\beta }_i)=\Lambda +O(N\Vert {\mathcal{P}}_N{\Vert }^3).$$

$$\Lambda \le E({\mathcal{P}}_N)=\Lambda +O(\Vert {\mathcal{P}}_N\Vert ).$$

$${\tau }_{\mathrm{SEO}}({\mathcal{P}}_N)$$

$${\tau }_{\mathrm{DEO}}({\mathcal{P}}_N)⟶\bar{\tau }=\frac{1}{2+2\Lambda }>0,$$

$$\frac{N}{1+N{r}^{\ast }}\le \frac{{\tau }_{\text{DEO}}({\mathcal{P}}_N)}{{\tau }_{\text{SEO}}({\mathcal{P}}_N)}=\frac{N}{1+(N-1)r^{\ast }}\le \frac{1}{r^{\ast }}.$$

$$\frac{1}{r^{\ast }}-\frac{N}{1+N{r}^{\ast }}=\frac{1}{r^{\ast }}\left(\frac{1}{1+N{r}^{\ast }}\right)\approx \frac{1}{r^{\ast }}\left(\frac{1}{1+\Lambda }\right).$$

$$\Lambda ({\beta }_k^{\ast })\approx \frac{k}{N}\Lambda ,$$

$${\hat{r}}^{(i-1,i)}=\frac{1}{n}\sum _{k=1}^n{\alpha }^{(i-1,i)}({\bar{\text{X}}}_k)=r^{(i-1,i)}+O_p(n^{-1/2}).$$

$$\hat{\Lambda }({\beta }_k)=\sum _{i=1}^k{\hat{r}}^{(i-1,i)},$$

Full text

1 Introduction

Markov Chain Monte Carlo (MCMC) methods are widely used to approximate expectations with respect to a probability distribution with density $\pi(x)$ known up to a normalizing constant, i.e., $\pi(x)=\gamma(x)/\mathcal{Z}$ where $\gamma$ can be evaluated pointwise but the normalizing constant $\mathcal{Z}$ is unknown. When $\pi$ has multiple well-separated modes, highly varying curvature or when one is interested in sampling over combinatorial spaces, standard MCMC algorithms can perform very poorly. This work is motivated by the need for practical methods for these difficult sampling problems. A natural direction to address them is to use multiple cores and to distribute the computation.

1.1 Parallel Tempering

One popular approach for multi-core and distributed exploration of complex distributions is Parallel Tempering (PT) which was introduced independently in statistics [Geyer, 1991] and physics [Hukushima and Nemoto, 1996]; see also [Swendsen and Wang, 1986] for an earlier related proposal. Since its inception, PT remains to this day a very popular MCMC method to sample from complex multimodal target distributions arising in physics, chemistry, biology, statistics, and machine learning; see, e.g., [Issaoun et al., 2021, Ballnus et al., 2017, Chandra et al., 2019, Cho et al., 2010, Desjardins et al., 2014, Diaz et al., 2020, Dorri et al., 2020, Kamberaj, 2020, Müller and Bouckaert, 2020].

To sample from the target distribution $\pi$, PT introduces a sequence of auxiliary tempered or annealed probability distributions with densities $\pi^{(\beta_{i})}(x)\propto L(x)^{\beta_{i}}\pi_{0}(x)$ for $i=0,1,...,N$, where $\pi_{0}$ is an easy-to-sample reference distribution, $L(x)=\pi(x)/\pi_{0}(x)$ and the sequence $0=\beta_{0}<\beta_{1}<\cdots<\beta_{N}=1$ defines the annealing schedule. This bridge of auxiliary distributions is used to progressively transform samples from the reference distribution ($\beta=0$) into samples from the target distribution ($\beta=1$), for which only poorly mixing MCMC kernels may be available. For example, in the Bayesian setting where the target distribution is the posterior, we can choose the reference distribution as the prior, from which we can often obtain independent samples.

More precisely, PT algorithms are based on Markov chains in which the states are $(N+1)$-tuples, $\textbf{x}=(x^{0},x^{1},x^{2},\dots,x^{N})\in\mathcal{X}^{N+1}$, and whose stationary distribution is given by $\boldsymbol{\pi}(\textbf{x})=\prod_{i=0}^{N}\pi^{(\beta_{i})}(x^{i})$ [Geyer, 1991]. At each iteration, PT proceeds by applying in parallel $N+1$ MCMC kernels targeting $\pi^{(\beta_{i})}$ for $i=0,...,N$. We call these model-specific kernels the local exploration kernels. The chains closer to the reference chain (i.e. those with annealing parameter $\beta_{i}$ close to zero) can typically traverse regions of low probability mass under $\pi$ while the chain at $\beta=1$ ensures that asymptotically we obtain samples from the target distribution. Frequent communication between the chains at the two ends of the spectrum is therefore critical for good performance, and achieved by proposing to swap the states of chains at adjacent annealing parameters. Even swap moves (rows labelled ‘E’), respectively Odd swap moves (labelled ‘O’), propose to exchange states at chains with an even index $i$, respectively odd index $i$, and $i+1$. These proposals are accepted or rejected according to a Metropolis mechanism.

1.2 Deterministic and stochastic even-odd schemes

A key notion used to analyze the behaviour of PT is the index process. To provide intuition on this process, it is helpful to discuss briefly how PT is distributed over several machines. An important point is that instead of having pairs of machines exchanging high-dimensional states when a swap is accepted (which could be detrimental due to network latency), the machines should just exchange the annealing parameters. Suppose now we initialize machine $j$ with annealing parameter $\beta_{j}$. Then after $n$ scans, the annealing parameters are permuted among the $N+1$ machines according to the permutation $\textbf{I}_{n}=(I_{n}^{0},\dots,I_{n}^{N})$ of $\{0,\dots,N\}$, so that machine $j$ has annealing parameter $\beta_{I_{n}^{j}}$ at iteration $n$ of PT. Each index process $I_{n}^{j}$, formally introduced in Section 2.4, is initialized using $I^{j}_{0}=j$ and tracks how the state of the corresponding chain evolves over the annealing schedule thanks to the swap moves; see Figure 1. The index process $I_{n}^{j}$ thus monitors the information transfer between the reference and target on machine $j$ and as such determines partly the effectiveness of PT.

There have been many proposals made to improve this information transfer by adjusting the annealing schedule; see, e.g., [Kone and Kofke, 2005, Atchadé et al., 2011, Miasojedow et al., 2013]. These proposals are useful but do not address a crucial limitation of standard PT algorithms. In a distributed context, one can select randomly at each iteration whether to apply Even or Odd swap moves in parallel. The resulting stochastic even-odd swap (SEO) scheme, henceforth referred to as reversible PT as it admits a reversible scaling limit (see Section 6), yields index processes exhibiting a diffusive behaviour; see top row of Figure 1. Hence we can expect that when $N$ is large it takes roughly $O(N^{2})$ swap attempts for a state at $\beta_{0}=0$ to reach $\beta_{N}=1$ [Diaconis et al., 2000]. The user thus faces a trade-off. If $N$ is too large, the acceptance probabilities of the swap moves are high but it takes a time of order $O(N^{2})$ for a state at $\beta=0$ to reach $\beta=1$. If $N$ is too low, the acceptance probabilities of swap moves deteriorate resulting in poor mixing between the different chains. Informally, even in a multi-core or distributed setting, for $N$ large, the $O(N)$ gains in being able to harness more cores do not offset the $O(N^{2})$ cost of the diffusion (see Section 3.5 where we formalize this argument). As a consequence, the general consensus is that the temperatures should be chosen to allow for about a 20–40% acceptance rate to maximize the squared jump distance travelled per swap in the space of annealing parameters $[0,1]$ [Rathore et al., 2005, Kone and Kofke, 2005, Lingenheil et al., 2009, Atchadé et al., 2011]. Adding more chains past this threshold actually deteriorates the performance of reversible PT and there have even been attempts to adaptively reduce the number of additional chains [Łacki and Miasojedow, 2016]. This is a lost opportunity, since PT is otherwise particularly suitable to implementation on multi-core or distributed architectures.

An alternative to the SEO scheme is the deterministic even-odd swap (DEO) scheme introduced in [Okabe et al., 2001] where one deterministically alternates Even and Odd swap moves. We refer to DEO as non-reversible PT as it admits a non-reversible scaling limit (see Section 6). In particular, the resulting index processes do not appear to exhibit a diffusive, i.e. random walk type, behaviour, illustrated by the bottom row of Figure 1. This non-diffusive behaviour of non-reversible PT explains its excellent empirical performance when compared to alternative reversible PT schemes observed in practice [Lingenheil et al., 2009] for any given schedule. However non-reversible PT, like reversible PT, is sensitive to the choice of the schedule. All the aforementioned tuning strategies developed for PT implicitly assume a reversible framework and do not apply to non-reversible PT (as empirically verified in [Lingenheil et al., 2009]). The main contribution of this paper is to identify some of the theoretical properties of non-reversible PT so as to establish optimal tuning guidelines for this algorithm and propose a novel schedule optimization scheme to implement them. In all our experiments, the resulting iteratively optimized non-reversible PT scheme markedly outperforms the reversible and non-reversible PT schemes currently available (see Section 7).

1.3 Overview of our contributions

After introducing formally the SEO and DEO schemes in Section 2, our first contribution is a non-asymptotic result showing that the non-reversible DEO scheme is guaranteed to outperform its reversible SEO counterpart. The notion of optimality we analyze is the round trip rate, which quantifies how often information from the reference distribution percolates to the target; see Section 3. The theoretical analysis is based on a simplifying assumption called Efficient Local Exploration (ELE), which is not expected to hold exactly in real scenarios. However we show empirically that there are practical methods to approximate ELE and that even when ELE is violated the key predictions made by the theory closely match empirical behaviour. In this sense ELE can be thought of as a useful model for understanding PT algorithms.

In Section 4 we introduce the local and global communication barrier $\lambda(\beta),\Lambda$ which encode the local and global efficiency of PT. We then show that for non-reversible PT the round trip rate converges to $(2+2\Lambda)^{-1}$, in contrast to the reversible counterpart for which it decays to zero. To provide some intuition on how the small algorithmic difference between SEO and DEO can have such a profound impact, consider the scenario where PT would use the same distributions at the two end-points, $\pi=\pi_{0}$, as well as for all intermediate distributions so that $\Lambda=0$. Clearly this is not a realistic scenario, but in this simple context the index processes of DEO and SEO are easy to describe and contrast. In both DEO and SEO, $\pi=\pi_{0}$ implies that all proposed swaps will be accepted. For DEO, this makes the index process fully deterministic, performing direct trips from index [math] to $N$ and back. Such a process could be compared to a “conveyor belt” with the property that no matter what is the value of $N$, one novel trip from chain [math] reaches chain $N$ every two iterations, i.e. the round trip rate is $(2+2\Lambda)^{-1}=1/2$ as $\Lambda=0$. For SEO, even when $\pi=\pi_{0}$ the index process is still random and can be readily seen to be a simple discrete random walk. This implies that as $N$ increases, the round trip rate decreases to zero.

In practice, achieving high round trip rates requires careful tuning of the annealing schedule. In Section 5 we combine the analysis from Section 3 and Section 4 to develop a novel methodology to optimize the annealing parameters. The optimal tuning guidelines we provide are different from existing (reversible) PT guidelines and the novel methodology is highly parallel as its performance does not collapse when a very large number of chains is used. However using a large number of chains does have a diminishing return, therefore we propose a mechanism to determine the optimal trade-off between the number of chains and the number of independent PT algorithms one should use.

In Section 6 we identify the scaling limit of the index processes for both reversible and non-reversible PT as the number of parallel chains goes to infinity. We show that this scaling limit is a piecewise-deterministic Markov process for non-reversible PT whereas it is a diffusion for reversible PT as suggested by the dynamics of the bold paths in Figure 1.

Finally in Section 7, we present a variety of experiments validating our theoretical analysis and novel methodology. The method is implemented in an open source Bayesian modelling language available at https://github.com/UBC-Stat-ML/blangSDK. Our software implementation allows the user to specify the model in BUGS-like language [Lunn et al., 2000]. From this model declaration, a suitable sequence of annealed distributions is instantiated and a schedule optimized using our iterative method.

2 Setup and notation

2.1 Parallel tempering

Henceforth we will assume that the target and reference probability distributions $\pi$ and $\pi_{0}$ on $\mathcal{X}$ admit strictly positive densities with respect to a common dominating measure $\mathrm{d}x$. We will also denote these densities somewhat abusively by $\pi$ and $\pi_{0}$. It will be useful to define $V_{0}(x)=-\log\pi_{0}(x)$ and $V(x)=-\log L(x)$, where $L(x)=\pi(x)/\pi_{0}(x)$ is assumed finite for all $x\in\mathcal{X}$. Using this notation, the annealed distribution at an annealing parameter $\beta$ is given by

[TABLE]

where $\mathcal{Z}(\beta)=\int_{\mathcal{X}}L(x)^{\beta}\pi_{0}(x)\mathrm{d}x$ is the corresponding normalizing constant. We denote the annealing schedule by $0=\beta_{0}<\beta_{1}<\cdots<\beta_{N}=1$. In our analysis we will view it as a partition $\mathcal{P}_{N}=\{\beta_{0},\dots,\beta_{N}\}$ of $[0,1]$ with mesh-size $\|\mathcal{P}_{N}\|=\sup_{i}\leavevmode\nobreak\ \{\beta_{i}-\beta_{i-1}\}$.

We define the PT state $\bar{\textbf{x}}=(\textbf{x},\textbf{i})$ where $\textbf{x}=(x^{0},\dots,x^{N})\in\mathcal{X}^{N+1}$ and $\textbf{i}=(i^{0},\dots,i^{N})\in\textrm{Perm}(N)$ the group of permutations on $\{0,\dots,N\}$. PT involves constructing a Markov Chain $\bar{\textbf{X}}_{n}=(\textbf{X}_{n},\textbf{I}_{n})$ over $\mathcal{X}^{N+1}\times\textrm{Perm}(N)$ invariant with respect to $\bar{\boldsymbol{\pi}}(\bar{\textbf{x}})=(N+1)!^{-1}\prod_{i=0}^{N}\pi^{(\beta_{i})}(x^{i})$. The $i$-th component of $\textbf{X}_{n}=(X_{n}^{0},\dots,X_{n}^{N})$ tracks the states associated with annealing parameter $\beta_{i}$ at iteration $n$ while the permutation $\textbf{I}_{n}=(I_{n}^{0},\dots,I_{n}^{N})$ tracks how the annealing parameters are shuffled among machines. In particular at iteration $n$, machine $j$ stores the $I^{j}_{n}$-th component of $\textbf{X}_{n}$ and annealing parameter $\beta_{I_{n}^{j}}$. The sequence of states associated with $\beta_{i}$ is called the $i$-th chain and the state on machine $j$ is the $j$-th replica.

For both SEO and DEO, the overall $\bar{\boldsymbol{\pi}}$-invariant Markov kernel ${\textbf{K}}^{\textrm{PT}}_{n}$ describing the algorithm is obtained by the composition of a $\bar{\boldsymbol{\pi}}$-invariant local exploration kernel ${\textbf{K}}^{\textrm{expl}}$ and communication kernel ${\textbf{K}}^{\textrm{comm}}_{n}$,

[TABLE]

The difference between SEO and DEO is in the communication phase, namely ${\textbf{K}}^{\textrm{comm}}_{n}={\textbf{K}}^{\textrm{SEO}}$ in the former case and ${\textbf{K}}^{\textrm{comm}}_{n}={\textbf{K}}^{\textrm{DEO}}_{n}$ in the latter. Markov kernels corresponding to the reversible (SEO) and non-reversible (DEO) PT algorithms are described informally in the introduction and illustrated in Figure 1.

2.2 Local exploration kernels

The local exploration kernels are defined in the same way for SEO and DEO. They are also model specific, so we assume we are given one $\pi^{(\beta_{i})}$-invariant kernel $K^{(\beta_{i})}$ for each annealing parameter $\beta_{i}\in\mathcal{P}_{N}$. These can be based on $n_{\textrm{expl}}$ steps of Hamiltonian Monte Carlo, Metropolis–Hastings, Gibbs Sampling, Slice Sampling, etc. We construct the overall local exploration kernel by applying the annealing parameter specific kernels to each component independently from each other:

[TABLE]

where $\delta_{\textbf{i}}$ denotes the Dirac delta.

In our computational model, we implicitly assume that the local exploration kernel at $\beta=0$ is special in that it can provide independent exact samples from $\pi_{0}$. Mathematically, $K^{(0)}(x,A_{0})=\pi_{0}(A_{0})$. This assumption is satisfied in most Bayesian models equipped with proper prior distributions, but also in other situations such as Markov random fields (see Appendix I.6).

2.3 Communication kernels.

Before defining the communication scheme, we first construct its fundamental building block, a swap. A swap is a Metropolis–Hastings move with a deterministic proposal which consists of swapping $\beta_{i}$ and $\beta_{i+1}$ across machines. From a current state $\bar{\textbf{x}}=(\textbf{x},\textbf{i})$, a proposed swap state is denoted $\bar{\textbf{x}}^{(i,i+1)}=(\textbf{x}^{(i,i+1)},\textbf{i}^{(i,i+1)})$ where

[TABLE]

and $\textbf{i}^{(i,i+1)}\in\textrm{Perm}(N)$ is the permutation obtained by swapping $i$ and $i+1$ in i. The Metropolis–Hastings kernel $\textbf{K}^{(i,i+1)}$ corresponding to this update is given by

[TABLE]

The function $\alpha^{(i,i+1)}(\bar{\textbf{x}})$ is the corresponding acceptance probability equal to

[TABLE]

and the maximal collection of adjacent swap kernels that can be proposed in parallel without interference are

[TABLE]

which we call the even and odd kernels respectively.

For SEO, the kernel ${\textbf{K}}^{\textrm{comm}}_{n}={\textbf{K}}^{\textrm{SEO}}$ is given by a mixture of the even and odd kernels in equal proportion while for DEO the kernel ${\textbf{K}}^{\textrm{comm}}_{n}={\textbf{K}}^{\textrm{DEO}}_{n}$ is given by a deterministic alternation between even and odd kernels, that is

[TABLE]

We provide pseudo-code for the DEO scheme in Algorithm 1. The pseudo-code also estimates the average rejection probabilities $r^{(i,i+1)}$ of swap moves between chains $i$ and $i+1$ which are used to optimize the annealing schedule in Section 5. When the schedule is fixed, lines 1, 12, 17 can be omitted, and one should use “for $i\in P$” on line 10. For simplicity, the swap in line 15 is shown for a parallel computing context, where several cores have a shared memory, and hence line 15 is simply an exchange of pointers in an array, which is efficient thanks to memory sharing. For a distributed computing implementation, where several machines do not share memory and instead need to communicate over the network, it becomes advantageous to swap annealing parameters instead of states. Refer to Appendix H for pseudo-code for the distributed computing implementation.

2.4 The index process

As discussed in Section 1.2, the $j$-th component of permutation $\textbf{I}_{n}=(I^{0}_{n},\dots,I^{N}_{n})$ encodes the flow of information between the reference and target on machine $j$. We define the index process for machine $j$ as $(I_{n}^{j},\varepsilon_{n}^{j})\in\{0,\dots,N\}\times\{-1,1\}$, where $\varepsilon_{n}^{j}=1$ if the annealing parameter on machine $j$ is proposed an increase after scan $n$, and $-1$ otherwise. The concept is best understood visually: refer to the bold piecewise linear paths illustrating $I^{N-2}_{n}$ for $N=8$ and $I^{N-1}_{n}$ for $N=30$ in Figure 1.

Let $P^{(i,i+1)}\in\{0,1\}$ denote an indicator that a swap is proposed between chains $i$ and $i+1$ at iteration $n$. The realized swaps are then defined from the proposal indicators as $S^{(i,i+1)}=P^{(i,i+1)}A^{(i,i+1)}$, where $A^{(i,i+1)}|\bar{\textbf{X}}\sim\mathrm{Bern}(\alpha^{(i,i+1)}(\bar{\textbf{X}}))$ are acceptance indicator variables (Figure 2 in Appendix J). The index process satisfies the following recursive relation: initialize $I^{j}_{0}=j$ and $\varepsilon_{0}^{j}=1$ if $P_{0}^{(j,j+1)}=1$ and $-1$ otherwise. For $n>0$, we have

[TABLE]

We will use the dynamics of the index process to explain the differences between SEO and DEO communication. The only difference between the two is in the proposal indicators. Define $\textbf{P}_{n}=(P_{n}^{(0,1)},P_{n}^{(1,2)},\dots,P_{n}^{(N-1,N)})$, $\textbf{P}_{n}$ is deterministic for DEO, i.e. $\textbf{P}_{n}=\textbf{P}_{\text{even}}=(1,0,1,\dots)$ for even $n$ and $\textbf{P}_{n}=\textbf{P}_{\text{odd}}=(0,1,0,\dots)$ for odd $n$. In SEO, we have $\textbf{P}_{n}\sim\text{Unif}\{\textbf{P}_{\text{even}},\textbf{P}_{\text{odd}}\}$.

For SEO, the variables $\varepsilon_{n}\sim\mathrm{Unif}\{-1,1\}$ are i.i.d, and consequently the index process exhibits a random walk behaviour. In contrast for DEO, we have $\varepsilon^{j}_{n+1}=\varepsilon^{j}_{n}$ so long as $I^{j}_{n+1}=I^{j}_{n}+\varepsilon^{j}_{n}$ and $\varepsilon^{j}_{n+1}=-\varepsilon^{j}_{n}$ otherwise. Therefore the index process for DEO performs a more systematic exploration of the space as the direction $\varepsilon^{j}_{n+1}$ is only reversed when a swap involving machine $j$ is rejected or if the boundary is reached. The qualitative differences between the two regimes can be seen in Figure 1 (see also Figure 15 in the Supplementary Material). In particular the index process for DEO in these figures behaves very differently as $N$ increases, this will be explored formally in Section 6.

As mentioned in the introduction, we refer to the PT algorithm with SEO and DEO communication as reversible PT and non-reversible PT respectively. Our terminology is somewhat abusive but is justified by the analysis in Section 3.4 and Section 6, where it is shown that, under certain assumptions, the index process is reversible for SEO while it is non-reversible for DEO.

3 Non-asymptotic analysis of PT algorithms

3.1 Model of compute time

We start with a definition of what we model as one unit of compute time: throughout the paper, we assume a massively parallel or distributed computational setup and sampling once from each of local exploration kernel $K^{(\beta)}$ has cost $O(n_{\textrm{expl}})$ which dominates the cost of the swap kernel ${\textbf{K}}^{(i,i+1)}$. Consequently a scan of PT has cost $O(n_{\textrm{expl}})$ and for a fixed computational budget $n_{\textrm{total}}$, the total number of scans is $n_{\textrm{scan}}=O(n_{\textrm{total}}/n_{\textrm{expl}})$.

The assumption that the per-iteration cost of PT is independent of the number of chains is reasonable in GPU and parallel computing scenarios, since the communication cost for each swap does not increase with the dimension of the problem (by swapping annealing parameters instead of states). We also assume that the number of PT scans will dominate the number of parallel cores available, i.e. $n_{\textrm{scan}}\gg N$. This is reasonable when addressing challenging sampling problems. Although there are numerous empirical studies on multi-core and distributed implementation of PT [Altekar et al., 2004, Mingas and Bouganis, 2012, Fang et al., 2014], we are not aware of previous theoretical work investigating such a computational model.

3.2 Performance metrics for PT methods

The standard notion of computational efficiency of MCMC schemes is the effective sample size (ESS) per compute time. However, for PT methods, since the ESS per compute time depends on the details of the problem specific local exploration kernels ${\textbf{K}}^{\textrm{expl}}$, alternatives have been developed in the literature to assess the performance of ${\textbf{K}}^{\textrm{comm}}_{n}$ which are independent of ${\textbf{K}}^{\textrm{expl}}$ [Katzgraber et al., 2006, Lingenheil et al., 2009].

We are motivated by the Bayesian context where it is typically possible to obtain one independent sample from the reference distribution $\pi_{0}$ (i.e. from the prior distribution) at each iteration. We say that an annealed restart has occurred on machine $j$ when $I_{n}^{j}$ goes from [math] to $N$ (i.e. $\beta$ goes from [math] to $1$), which corresponds to a sample generated from $\pi_{0}$ propagating to the target $\pi$. Informally an annealed restart can be thought of as a sampling equivalent to what is known in optimization as a random restart. We say a round trip has occurred on machine $j$ when $I_{n}^{j}$ goes from [math] to $N$ and then goes back to [math] (i.e. $\beta$ goes from [math] to $1$ to [math]).

Formally, we recursively define $T_{\downarrow,0}^{j}=\inf\{n:(I_{n}^{j},\varepsilon_{n}^{j})=(0,-1)\}$ and for $k\geq 1$,

[TABLE]

The $k$-th annealed restart and round trip for machine $j$ occurs at scan $T_{\uparrow,k}^{j}$ and $T_{\downarrow,k}^{j}$ respectively. Let $\mathcal{T}_{n}$ and $\mathcal{R}_{n}$ be the total number of annealed restarts and round trips respectively during the first $n$ scans of the DEO algorithm.

We wish to optimize for the percentage of iterations that result in an annealed restart, i.e. $\tau=\lim_{n\to\infty}{\mathbb{E}}[\mathcal{T}_{n}]/n$, where we use abusively the same random variables for SEO and DEO but differentiate these schemes by using the probability measures $\mathbb{P}_{\textrm{SEO}}$ and $\mathbb{P}_{\textrm{DEO}}$ with associated expectation operators $\mathbb{E}_{\textrm{SEO}}$ and $\mathbb{E}_{\textrm{DEO}}$. We use $\mathbb{P}$ and ${\mathbb{E}}$ for statements that hold for both algorithms. If $\mathcal{T}_{n}^{j}$ and $\mathcal{R}_{n}^{j}$ are the total number of annealed restarts and round trips during the first $n$ iterations on machine $j$ respectively, then we have $\mathcal{R}_{n}^{j}\leq\mathcal{T}_{n}^{j}\leq\mathcal{R}_{n}^{j}+1$. Consequently, $\mathcal{R}_{n}\leq\mathcal{T}_{n}\leq\mathcal{R}_{n}+N+1$ and thus $\tau=\lim_{n\to\infty}{\mathbb{E}}[\mathcal{R}_{n}]/n$.

In the PT literature, $\tau$ is commonly referred to as the round trip rate and has been used to compare the effectiveness of various PT algorithms [Katzgraber et al., 2006, Lingenheil et al., 2009]. Empirically we observed that round trips per unit cost strongly correlate with ESS per unit cost as seen in Figure 3, making the round trip rate a natural objective function to compare and tune parallel tempering algorithms.

Another performance metric commonly used in the PT literature is the expected square jump distance (ESJD) [Atchadé et al., 2011, Kone and Kofke, 2005]. While this criterion is useful within the context of reversible PT for selecting the optimal number of parallel chains, the ESJD is too coarse to compare reversible to non-reversible PT methods as, for any given annealing schedule, the ESJD is identical in both cases.

3.3 Model assumptions

The analysis of the round trip times is in general intractable because the index process is not Markovian. Indeed, simulating a transition depends on the swap indicators $S_{n}^{(i,i+1)}$ (see Section 2.4), the distributions of which themselves depend on the state configuration $\bar{\textbf{X}}$. To simplify the analysis, we will make in the remainder of the paper the following simplifying assumptions:

1. (A1)

Stationarity: $\bar{\textbf{X}}_{0}\sim\bar{\boldsymbol{\pi}}$ and thus $\bar{\textbf{X}}_{n}\sim\bar{\boldsymbol{\pi}}$ for all $n$ as the kernel ${\textbf{K}}^{\textrm{PT}}_{n}$ is $\bar{\boldsymbol{\pi}}$-invariant. 2. (A2)

Efficient Local Exploration (ELE): For $X\sim\pi^{(\beta)}$ and $X^{\prime}|X\sim K^{(\beta)}(X,\cdot)$, the random variables $V(X)$ and $V(X^{\prime})$ are independent. 3. (A3)

Integrability: $V^{3}$ is integrable with respect to $\pi_{0}$ and $\pi$.

It follows from Assumptions (A1)–(A2) and (6) that the behaviour of the communication scheme only depends on the distribution of the state $\bar{\textbf{X}}_{n}$ via the $N+1$ univariate distributions of the chain-specific energies $V^{(i)}=V\left(X^{(i)}\right)$, $i\in\{0,1,2,\dots,N\}$. This allows us to build a theoretical analysis which makes no structural assumption on the state space $\mathcal{X}$ or the target $\pi$ as typically done in the literature: for example, [Atchadé et al., 2011] assume a product space $\mathcal{X}=\mathcal{X}_{0}^{d}$ for large $d$, and [Predescu et al., 2004] assume $\pi^{(\beta)}$ satisfies a constant heat capacity.

Admittedly the ELE assumption (A2) does not hold in practical applications. ELE can be approximated by increasing the number of local exploration kernels applied between consecutive swap ($n_{\textrm{expl}}$). However one may worry that to achieve a good approximation in challenging problems, $n_{\textrm{expl}}$ would have to be set to a value so large as to defy the practicality of our analysis. Surprisingly, we have observed empirically that this was not the case in the multimodal problems we considered. Figure 4 displays results in four models where a local exploration kernel alone induces good mixing of the energy chain $V(X_{n})$ (hence ELE can be approximated) yet the local exploration kernel alone is insufficient to achieve good mixing on the full state space, $\bar{\textbf{X}}_{n}$ (so that PT is justified and indeed yields efficient exploration of the configuration space). This gap is possible since $V(X)$ is 1-dimensional and potentially unimodal even when $X$ is not. This is the motivation for ELE since assuming the independence of $V(X)$ and $V(X^{\prime})$ is weaker than assuming the independence of $X$ and $X^{\prime}$ (as hypothesized e.g. in Section 5.1 of [Atchadé et al., 2011]). Obviously ELE is still expected to be a somewhat crude simplifying assumption in very complex problems; e.g. for the highly challenging high-dimensional copy number inference problem illustrated in Figure 4.

In Section 7.2, we describe additional empirical results supporting that ELE is a useful model for the purpose of designing and analyzing PT algorithms. In the severe ELE violation regime, we show that the key quantities used in our analysis are either well approximated, or approached as $N$ increases.

Assumptions (A1)-(A2) allow us to express the swap indicators as independent Bernoulli random variables $S_{n}^{(i,i+1)}\sim\text{Bern}(s^{(i,i+1)})$ where $s^{(i,i+1)}$ is given by the expectation of Equation (6),

[TABLE]

the expectation being over two independent random variables $V^{(i)},V^{(i+1)}$, satisfying $V^{(i)}\stackrel{{\scriptstyle d}}{{=}}V(X^{(\beta_{i})})$ for $X^{(\beta_{i})}\sim\pi^{(\beta_{i})}$.

3.4 Reversibility and non-reversibility of the index process

Under assumptions (A1)-(A2), each index process $(I^{j}_{n},\varepsilon^{j}_{n})$ is Markovian for $j=0,\dots,N$ with transition kernel $K^{\mathrm{SEO}}$ and $K^{\mathrm{DEO}}$ for reversible and non-reversible PT respectively. See Appendix A for an explicit representation of these kernels. We drop the superscript $j$ when the particular machine is not relevant.

The kernel $K^{\mathrm{SEO}}$ defines a reversible Markov chain on $\{0,\dots,N\}\times\{-1,1\}$ with uniform stationary distribution while $K^{\mathrm{DEO}}$ satisfies the skew-detailed balance condition with respect to the same distribution,

[TABLE]

and is thus non-reversible. It falls within the generalized Metropolis–Hastings framework, see, e.g., [Lelièvre et al., 2010].

Reversibility necessitates that the Markov chain must be allowed to backtrack its movements. This leads to inefficient exploration of the state space. As a consequence, non-reversibility is typically a favourable property for MCMC chains. A common recipe to design non-reversible sampling algorithms consists of expanding the state space to include a “lifting” parameter that allows for a more systematic exploration of the state space [Chen et al., 1999, Diaconis et al., 2000, Turitsyn et al., 2011, Vucelja, 2016].

The index process $(I_{n},\varepsilon_{n})$ for non-reversible PT can be interpreted as a “lifted” version of the index process for reversible PT with lifting parameter $\varepsilon_{n}$. Under DEO communication, $I_{n}$ travels in the direction $\varepsilon_{n}$ and only reverses direction when $I_{n}$ reaches a boundary or when a swap rejection occurs. This “lifting” construction helps explain the qualitatively different behaviour between reversible and non-reversible PT and will be further explored when identifying the scaling limit of $(I_{n},\varepsilon_{n})$ in Section 6. The lifted PT of [Wu, 2017] exploits a similar construction but only one of the $N+1$ index processes is lifted instead of all of them for DEO. A lifted version of simulated tempering has also been proposed by [Sakai and Hukushima, 2016].

3.5 Non-asymptotic domination of non-reversible PT

Assumptions (A1)-(A2) ensure that for each $j=0,\dots,N$, $\mathcal{R}_{n}^{j}$ is a delayed renewal process with round trip times $T_{k}^{j}=T^{j}_{\downarrow,k}-T^{j}_{\downarrow,k-1}$ for $k\geq 1$ and $j=0,\dots,N$. In particular, $T_{k}^{j}$ are independent and identically distributed and, for convenience, we introduce the random variable $T\stackrel{{\scriptstyle d}}{{=}}T_{k}^{j}$. By the key renewal theorem, we have

[TABLE]

An analytical expression for ${\mathbb{E}}[T]$ for reversible PT was first derived by [Nadler and Hansmann, 2007]. We extend this result to non-reversible PT in Theorem 1.

Theorem 1**.**

For any annealing schedule $\mathcal{P}_{N}=\{\beta_{0},\dots,\beta_{N}\}$,

[TABLE]

where $E(\mathcal{P}_{N})=\sum_{i=1}^{N}r^{(i-1,i)}/s^{(i-1,i)}$, and $r^{(i-1,i)}=1-s^{(i-1,i)}$ is the probability of rejecting a swap between chains $i$ and $i+1$.

The proof can be found in Appendix B.

Intuitively, Theorem 1 implies ${\mathbb{E}}[T]$ can be decomposed as the independent influence of communication scheme ${\textbf{K}}^{\textrm{comm}}_{n}$ and schedule $\mathcal{P}_{N}$ respectively. When all proposed swaps are accepted (i.e. $\pi=\pi_{0}$), the index process for reversible PT reduces to a simple random walk on $\{0,\dots,N\}$, whereas for non-reversible PT, the index processes takes a direct path from [math] to $N$ and back. Therefore, the first term in (15) and (16) represents the expected time for a round trip to occur in this idealized, rejection-free setting. The second term of (15) and (16) are identical and represent the additional time required to account for rejected swaps under schedule $\mathcal{P}_{N}$. Motivated by Theorem 1, we will refer to $E(\mathcal{P}_{N})$ as the schedule inefficiency.

By applying Theorem 1 to Equation (14), we get a non-asymptotic formula for the round trip rate in terms of $E(\mathcal{P}_{N})$.

Corollary 1**.**

For any annealing schedule $\mathcal{P}_{N}$ we have

[TABLE]

Consequently, $\tau_{\mathrm{DEO}}(\mathcal{P}_{N})>\tau_{\mathrm{SEO}}(\mathcal{P}_{N})$ for $N>1$.

Corollary 1 implies that non-reversible PT dominates reversible PT for any $N>1$ and any annealing schedule $\mathcal{P}_{N}$.

4 Asymptotic analysis of PT algorithms

4.1 The communication barrier

We begin by analyzing the behaviour of the PT swaps as $\|\mathcal{P}_{N}\|=\max_{i}|\beta_{i}-\beta_{i-1}|$ goes to zero. We define the swap and rejection functions $s,r:[0,1]^{2}\to[0,1]$ respectively as,

[TABLE]

where $V^{(\beta)}\stackrel{{\scriptstyle d}}{{=}}V(X^{(\beta)})$ for $X^{(\beta)}\sim\pi^{(\beta)}$ and $V^{(\beta)},V^{(\beta^{\prime})}$ are independent. The quantities $s(\beta,\beta^{\prime})$ and $r(\beta,\beta^{\prime})$ are symmetric in their arguments and represent the probability of swap and rejection occurring between $\beta$ and $\beta^{\prime}$ respectively under the ELE assumption (A2). Note that $s^{(i-1,i)}=s(\beta_{i-1},\beta_{i})$.

To take the limit as $\|\mathcal{P}_{N}\|\to 0$, it will be useful to understand the behaviour of $r(\beta,\beta^{\prime})$ when $\beta\approx\beta^{\prime}$. The key quantity that drives this asymptotic regime is given by a function $\lambda:[0,1]\to[0,\infty)$ defined as the instantaneous rate of rejection of a proposed swap at annealing parameter $\beta$,

[TABLE]

See Figure 7 (center) for examples of estimated $\lambda$ from various models. We define the integral of $\lambda$ by $\Lambda(\beta)=\int_{0}^{\beta}\lambda(\beta^{\prime})\mathrm{d}\beta^{\prime}$ and denote $\Lambda=\Lambda(1)$. Extending Proposition 1 in [Predescu et al., 2004] provides the following result.

Theorem 2**.**

The instantaneous rejection rate $\lambda$ is twice continuously differentiable and is equal to

[TABLE]

where $V_{1}^{(\beta)},V_{2}^{(\beta)}$ are independent random variables with law identical to the law of $V^{(\beta)}$. Moreover, we have

[TABLE]

See Appendix C for a proof and general smoothness properties of $\lambda$. In particular, the existence and smoothness of $\lambda$ is guaranteed by Assumption (A3).

Theorem 2 shows that $\lambda$ encodes up to second order the behaviour of $r$ as the annealing parameter difference between the chains goes to 0. When $\lambda(\beta)$ is high, swaps are much more likely to be rejected, implying $\lambda(\beta)$ measures the difficulty of local communication for a chain with annealing parameter $\beta$.

Notice that $\Lambda\geq 0$ with equality if and only if $\lambda(\beta)=0$ for all $\beta\in[0,1]$. It can be easily verified from (22) that $\lambda=0$ if and only if $V^{(\beta)}$ is constant $\pi^{(\beta)}$-a.s. for all $\beta\in[0,1]$ which happens precisely when $\pi_{0}=\pi$. So $\Lambda$ defines a natural symmetric divergence measuring the difficulty of communication between $\pi_{0}$ and $\pi$. In particular, for any schedule $\mathcal{P}_{N}$, the sum of the rejection rates is approximately constant and equal to $\Lambda$ as formalized in Corollary 2.

Corollary 2**.**

For any schedule $\mathcal{P}_{N}$, we have

[TABLE]

Motivated by Theorem 2 and Corollary 2 we will henceforth refer to $\lambda$, $\Lambda$ and $\Lambda(\cdot)$ as the local, global and cumulative communication barriers respectively.

4.2 Asymptotic analysis of round trip rate

The communication barrier allows us to study the asymptotic performance of PT when the number of parallel chains is large. When $N$ increases, the round trip rate for reversible PT decays to [math] while it converges to a positive constant $\bar{\tau}$ for non-reversible PT as seen in left plot in Figure 3. We formally show this in Theorem 3.

Theorem 3**.**

As $\|\mathcal{P}_{N}\|\to 0$ we have:

1. (a)

The communication inefficiency satisfies

[TABLE] 2. (b)

The round trip rate for reversible PT, $\tau_{SEO}$, goes to zero:

[TABLE] 3. (c)

The round trip rate for non-reversible PT, $\tau_{DEO}$ satisfies

[TABLE]

where the converges occurs with rate $O(\|\mathcal{P}_{N}\|)$.

See Appendix D for a proof.

In general, $\Lambda$ is large when $\pi_{0}$ deviates significantly from $\pi$. Since $\Lambda$ is problem specific, this identifies a limitation of PT present even in its non-reversible flavour, namely that adding more cores to the task will never be harmful, but does have a diminishing return. The bound $(2+2\Lambda)^{-1}$ could indeed be very small for complex problems. Moreover, it is independent of the choice of annealing schedule, hence it cannot be improved by the schedule optimization procedure described in Section 5.

When $N$ is fixed, we will see in Section 5.1 that $\tau_{\textrm{DEO}},\tau_{\textrm{SEO}}$ are both maximized when $r^{(i-1,i)}=r^{*}$ for all $i$. In this case, $E(\mathcal{P}_{N})=N(1-r^{*})^{-1}\leavevmode\nobreak\ r^{*}$ and one has

[TABLE]

This implies that for an optimally chosen schedule, non-reversible PT improves the round trip rate over reversible PT by at most a factor of $1/r^{*}$. Corollary 2 implies $Nr^{*}\approx\Lambda$, allowing us to quantify the looseness of (28) directly in terms of $\Lambda$.

[TABLE]

We can also investigate the performance of PT for multi-modal or high-dimensional targets by studying the behaviour of $\Lambda$. Assuming the target distribution factorizes [Woodard et al., 2009, Atchadé et al., 2011, Roberts and Rosenthal, 2014], we can show that $\Lambda$ is robust to the number of modes and grows at rate $O(\sqrt{d})$ as the dimension $d$ of the target increases. See Appendices E and F for details.

5 Tuning non-reversible PT

5.1 Optimal annealing schedule

We first discuss how to optimize the annealing schedule $\mathcal{P}_{N}$ to maximize $\tau_{\textrm{DEO}}$ when $N$ is fixed, which is equivalent to minimizing the schedule inefficiency, $E(\mathcal{P}_{N})=\sum_{i=1}^{N}r_{i}/(1-r_{i})$ where $r_{i}=r(\beta_{i-1},\beta_{i})$. To get a tractable approximate characterization of the feasible region of $r_{1},r_{2},\dots,r_{N}$, we use Corollary 2, which implies $\sum_{i=1}^{N}r_{i}\approx\Lambda$ for all schedules $\mathcal{P}_{N}$. Therefore assuming $\|\mathcal{P}_{N}\|$ is small enough to ignore the error term, finding the optimal schedule $\mathcal{P}_{N}^{*}$ is approximately equivalent to minimizing $\sum_{i=1}^{N}r_{i}/(1-r_{i})$, subject to the constraint $\sum_{i=1}^{N}r_{i}=\Lambda$ and $r_{i}>0$. Using Lagrange multipliers this leads to a solution where the rejection probabilities $r_{i}^{*}$ are constant in $i$, with a shared value denoted by $r^{*}$.

Consequently the optimal schedule $\mathcal{P}_{N}^{*}=\{\beta_{0}^{*},\dots,\beta_{N}^{*}\}$ satisfies $r(\beta^{*}_{i-1},\beta^{*}_{i})=r^{*}$ for all $i$. Theorem 2 and Corollary 2 imply that $r^{*}$ must satisfy, $r^{*}\approx\Lambda(\beta_{i}^{*})-\Lambda(\beta_{i-1}^{*})$ for all $i$ and $r^{*}\approx\Lambda/N$ with $O(N^{-3})$ error. By equating these two estimates for $r^{*}$ and summing from $i=0,\dots,k$ we get

[TABLE]

with an error of $O(N^{-2})$. If we ignore error terms, (30) implies that $\beta_{k}^{*}\approx G(k/N)$ where $G=F^{-1}$ and $F(\beta)=\Lambda(\beta)/\Lambda$ (see Figure 5(2)). See Appendix I.4.1 for an example where an analytic result is available, however in general $\Lambda(\cdot)$ is not known but instead estimated from the MCMC output.

The “equi-acceptance” result in (30) is consistent with other theoretical frameworks and notions of efficiency [Atchadé et al., 2011, Lingenheil et al., 2009, Kofke, 2002, Predescu et al., 2004]. However implementing this equi-acceptance recommendation in practice is non-trivial. Previous work relied on Robbins-Monro schemes [Atchadé et al., 2011, Miasojedow et al., 2013], which introduce sensitive tuning parameters. In contrast, provided $\Lambda(\cdot)$ can be estimated reliably (which we establish with Algorithm 3 in the next section), (30) provides the basis for a straightforward and effective schedule optimization scheme, described in Algorithm 2.

5.2 Estimation of the communication barrier

Assume we have access to a collection of samples $(\bar{\textbf{X}}_{1},\bar{\textbf{X}}_{2},\dots,\bar{\textbf{X}}_{n})$ from a non-reversible PT scheme based on an arbitrary annealing schedule $\mathcal{P}_{N}$. These samples may come from a short pilot run, or, as described in the next section, from the previous round of an iterative scheme. For a given schedule $\mathcal{P}_{N}$, when the central limit theorem for Markov chains holds, the Monte Carlo estimates for the rejection rates satisfy

[TABLE]

Next, using Theorem 2 we obtain $\sum_{i=1}^{k}r^{(i-1,i)}=\Lambda(\beta_{k})+O(N\|\mathcal{P}_{N}\|^{3})$. This motivates the following approximation for $\Lambda(\beta_{k})$,

[TABLE]

which has an error of order $O_{p}(\sqrt{N/n}+N\|\mathcal{P}_{N}\|^{3})$ (see Figure 5 (left)).

Given $\hat{\Lambda}(\beta_{0}),\dots,\hat{\Lambda}(\beta_{N})$, we estimate the function $\Lambda(\beta)$ via interpolation, with the constraint that the interpolated function should be monotone increasing since $\lambda(\beta)\geq 0$. Specifically, we use the Fritsch-Carlson monotone cubic spline method [Fritsch and Carlson, 1980] and denote the monotone interpolation by $\hat{\Lambda}(\beta)$. While we only use $\hat{\Lambda}(\beta)$ in our schedule optimization procedure, it is still useful to estimate $\lambda(\beta)$ for visualization purpose. We use the derivative of our interpolation, $\hat{\lambda}(\beta)=\hat{\Lambda}^{\prime}(\beta)$, which is a piecewise quadratic function.

The ideas described in this section so far are summarized in Algorithm 3, which given rejection statistics collected for a fixed annealing schedule provides an estimate of the communication barrier.

As a byproduct of Algorithm 3 we also obtain a consistent estimator $\hat{\tau}=(2+2\hat{\Lambda})^{-1}$ for the optimal round trip rate $\bar{\tau}$, where $\hat{\Lambda}=\hat{\Lambda}(1)$. We show in Figure 8 an example where (A2) is severely violated, yet $\lambda$ is still reliably estimated. This allows us to compare the empirically observed round trip rate against $\hat{\tau}$, and hence estimate how far an implementation deviates from optimal performance.

5.3 Tuning $N$

Theorem 3 shows that non-reversible PT does not deteriorate in performance as we increase $N$ unlike reversible PT, however the gains in round trip rate eventually become marginal. When $N$ is very large we expect to accumulate more round trips by running $k>1$ parallel copies of PT. As we shall see, the large $N$ asymptotic is still required however in order to determine the optimal number $k^{*}$ of PT copies.

Suppose there are $k$ copies of PT running in parallel consisting of $N+1$ chains with optimal annealing schedule $\mathcal{P}_{N}$. If $\bar{N}$ is the total number of cores available then $k$ and $N$ satisfy the constraint $k(N+1)=\bar{N}$. By Corollary 1 the total round trip rate across all $k$ copies of PT is

[TABLE]

From Section 5.1, the optimal schedule $\mathcal{P}_{N}^{*}$ has a corresponding swap rejection rate $r^{*}=\Lambda/N+O(N^{-3})$. Substituting this into (33) we get $\tau=\tau_{\Lambda}(N)+O(N^{-1})$, where

[TABLE]

Ignoring error terms, $\tau_{\Lambda}(\cdot)$ is optimized when we run $k^{*}=\bar{N}/(1+N^{*})$ copies of PT with $N^{*}+1=2\Lambda+1$ chains to achieve an optimum round trip rate $\tau^{*}=k^{*}/(2+4\Lambda)$.

Note that when $\tau$ is optimized, we have $r^{*}\approx 1/2$, which differs from from the $0.77$ optimal rejection rate from the reversible PT literature [Atchadé et al., 2011, Kofke, 2002, Predescu et al., 2004]. In Section 7.2, we show empirically that when the ELE assumption (A2) is severely violated, the recommendation in the present section is turned into a bound, $r^{*}<1/2$, as increasing $N$ appears to alleviate ELE violations.

5.4 Iterative schedule optimization

Suppose we have $\bar{N}$ cores available and a computational budget of $n_{\textrm{scan}}$ scans of PT, one scan consisting in one application of ${\textbf{K}}^{\textrm{PT}}$. The first $n_{\textrm{tune}}$ scans are used to find an accurate estimate of the communication barrier $\hat{\Lambda}(\cdot)$ and the remaining $n_{\textrm{sample}}=n_{\textrm{scan}}-n_{\textrm{tune}}$ scans to sample from the target. Algorithm 4 approximates $\hat{\Lambda}$ with error $O_{p}(\sqrt{\bar{N}/n_{\textrm{tune}}}+\bar{N}^{-2})$ by iteratively using Algorithms 2 and 3 for maxRound$=\log_{2}n_{\textrm{tune}}$ rounds until the tuning budget is depleted. Equipped with $\hat{\Lambda}(\cdot)$, we can use the remaining $n_{\textrm{sample}}$ scans to run $k^{*}=\bar{N}/(N^{*}+1)$ copies of non-reversible PT with $N^{*}=2\hat{\Lambda}(1)$ chains with optimal schedule $\mathcal{P}_{N^{*}}^{*}$ from Algorithm 2. In cases where it is suspected that (A2) may be severely violated (as evidenced for example by a large gap between $\bar{\tau}$ and the estimated $\tau$ as described in Section 5.2), it may be advantageous to also attempt $N^{*}=\bar{N}$ in line 9 of Algorithm 4. We show empirically in Figure 7 (top) that in a range of synthetic and real world problems our scheme converges using a small number of iterations, in the order of maxRound $=10$. In our experiments we used $n_{\textrm{tune}}=n_{\textrm{scan}}/2$ and discarded the samples from the first $n_{\textrm{tune}}$ scans as burn-in.

5.5 Normalizing constant computation

Algorithm 4 can be easily modified to obtain efficient estimators for quantities of the form $\int_{0}^{1}{\mathbb{E}}[f(V^{(\beta)},\beta)]\mathrm{d}\beta$ with an error of $O_{p}(\sqrt{\bar{N}/n_{\textrm{tune}}}+\bar{N}^{-2})$ for well-behaved functions $f$. Examples of such quantities include log-normalizing constants [Kirkwood, 1935, Gelman and Meng, 1998, Xie et al., 2011], KL-divergence [Dabak and Johnson, 2002], and Fisher-Rao distance [Amari, 2016] between $\pi$ and $\pi_{0}$.

For example, by taking a Riemann sum in the thermodynamic integration identity [Ogata, 1989], the log-normalizing constant $\log\mathcal{Z}(1)$ of $\pi$ can be approximated using

[TABLE]

where $\mu(\beta)={\mathbb{E}}[V^{(\beta)}]$ and $\mathcal{Z}(0)$ the normalizing constant of $\pi_{0}$. By running DEO with schedule $\mathcal{P}_{N}$ for $n$ scans we substitute a Monte Carlo estimate $\hat{\mu}(\beta_{k})$ for $\mu(\beta_{k})$ into (35) to get an estimator $\log(\hat{\mathcal{Z}}(1)/\hat{\mathcal{Z}}(0))$ with error $O_{p}(\sqrt{N/n}+N^{-2})$. See also [Gelman and Meng, 1998, Xie et al., 2011] for closely related normalization constant estimators. Figure 7 (bottom) in Section 7.1 shows how the estimate for $\log\mathcal{Z}(1)$ evolves with the number of rounds in Algorithm 4 for 16 different models.

6 Scaling limits of index process

The non-reversible communication scheme introduced in [Okabe et al., 2001] was presumably devised on algorithmic grounds (it performs the maximum number of swap attempts in parallel) since no theoretical justification was provided and the non-reversibility of the scheme was not mentioned. The arguments given in [Lingenheil et al., 2009] to explain the superiority of non-reversible communication over various PT algorithms rely on a misleading assumption, namely a diffusive scaling limit for the index process. Figure 1 (see also Figure 15 in the Supplementary Material) suggests that the index process behaves qualitatively differently as $N$ increases for reversible and non-reversible PT. The goal of this section is to investigate these differences by identifying some scaling limits. Such limits exist under assumptions (A1)-(A3) specified in Section 3.3 and Section 4.1.

Suppose $G:[0,1]\to[0,1]$ is an increasing differentiable function satisfying $G(0)=0$ and $G(1)=1$. We say that $G$ is a schedule generator for $\mathcal{P}_{N}=\{\beta_{0},\dots,\beta_{N}\}$ if $\beta_{i}=G(i/N)$. We will now assume without loss of generality that the sequence of schedules $\mathcal{P}_{N}$ are generated by some common $G$. In particular the mean value theorem implies $\|\mathcal{P}_{N}\|=O(N^{-1})$ as $N\to\infty$. This is not a strict requirement as most annealing schedules commonly used fall within this framework: the uniform schedule $\mathcal{P}_{N}^{\text{uniform}}=\{0,1/N,\dots,1\}$ is generated by $G(w)=w$, the optimal schedule $\mathcal{P}_{N}^{*}=\{\beta_{0}^{*},\dots,\beta_{N}^{*}\}$ derived in Section 5.1 is approximately generated by $G(w)=F^{-1}(w)$, where $F(\beta)=\Lambda(\beta)/\Lambda$. If $\pi_{0}(x)\propto\pi(x)^{\gamma}$ for some $\gamma\in(0,1)$, and $L(x)\propto\pi(x)^{1-\gamma}$ then $G(w)=\frac{\gamma^{1-w}-\gamma}{1-\gamma}$ corresponds to the geometric schedule commonly used by practitioners. In contrast, we show in Figure 20 examples of schedule generators $G$ corresponding to estimated optimal $\mathcal{P}_{N}^{*}$ from various models. Figure 20 shows that in real world inference scenarios, the optimal $G$ often qualitatively differs from the simple geometric schedule.

Suppose $(I_{n},\varepsilon_{n})$ is the index process for an annealing schedule $\mathcal{P}_{N}$ generated by $G$. To establish a scaling limit for $(I_{n},\varepsilon_{n})$, it will be convenient to work in a continuous time setting. To do this, we suppose the times that PT iterations occur are distributed according to a Poisson process $\{M(\cdot)\}$ with mean $\mu_{N}$. The number $M(t)$ of PT iterations that occur by time $t\geq 0$ thus satisfies $M(t)\sim\mathrm{Poisson}(\mu_{N}t)$. We define the scaled index process by $Z^{N}(t)=(W^{N}(t),\varepsilon^{N}(t))$ where $W^{N}(t)=I_{M(t)}/N$ and $\varepsilon^{N}(t)=\varepsilon_{M(t)}$.

Define the piecewise-deterministic Markov process (PDMP) [Davis, 1993] $Z(t)=(W(t),\varepsilon(t))$ on $[0,1]\times\{-1,1\}$ as follows: $W(t)$ moves in $[0,1]$ with velocity $\varepsilon(t)$ and the sign of $\varepsilon(t)$ is reversed at an inhomogeneous rate $\lambda(G(W(t))G^{\prime}(W(t))$ or when $W(t)$ hits a boundary; see [Bierkens et al., 2018] for a discussion of PDMP on restricted domains. The process $Z$ corresponds to an inhomogeneous persistent random walk with reflective boundary conditions; see [Masoliver et al., 1992] for a description of the Fokker-Plank equation and first passage times. As the optimal schedule is generated by $G=F^{-1}$ for $F(\beta)=\Lambda(\beta)/\Lambda$, we have $\lambda(G(w))G^{\prime}(w)=\Lambda$ for all $w\in[0,1]$, so $\varepsilon(t)$ changes sign at a constant rate $\Lambda$. See Figure 13 for sample trajectories of $Z$ for different values of $\Lambda$. The main result of this section is the following theorem.

Theorem 4**.**

Suppose $\mathcal{P}_{N}$ is family of schedules generated by $G$, then

1. (a)

For reversible PT if $\mu_{N}=N^{2}$ and if $W^{N}(0)$ converges weakly to $W(0)$ then $W^{N}$ converges weakly to a diffusion $W$ with initial condition $W(0)$, where $W$ is a Brownian motion on $[0,1]$ with reflective boundary conditions. The process $W$ admits $\mathrm{Unif}([0,1])$ as stationary distribution. 2. (b)

For non-reversible PT if $\mu_{N}=N$ and if $Z^{N}(0)$ converges weakly to $Z(0)$, then $Z^{N}$ converges weakly to the PDMP $Z$ with initial condition $Z(0)$. The process $Z$ admits $\mathrm{Unif}([0,1]\times\{-1,1\})$ as stationary distribution.

See Appendix G for a detailed construction of the scaled index process via their infinitesimal generators, and the proof of Theorem 4.

Theorem 4 implies for reversible PT that, if we speed time by a factor of $N^{2}$, then the index process scales to a diffusion $W$ independent of the choice of $\pi_{0},\pi$ and of the schedule. This is in contrast to non-reversible PT where, if we speed time by a factor of $N$, the index process converges to a PDMP $Z$ depending on $\pi_{0},\pi$ through $\lambda$ and the schedule through $G$.

7 Examples

7.1 Empirical behaviour of the schedule optimization method

We applied the proposed NRPT algorithm to 16 models from the statistics and physics literature to demonstrate its versatility. The models include 9 Bayesian models ranging from simple standard models such as generalized linear models and Bayesian mixtures, to complex ones such as cancer copy-number calling, ODE parameter estimation, spike-and-slab classification and two types of phylogenetic models. This is complemented by three models from statistical mechanics and four artificial models. The datasets considered include eight real datasets spanning diverse data-types and size, e.g. state-of-the-art measurements such as whole-genome single-cell sequencing data (494 individual cells from two types of cancer, triple negative breast cancer and high-grade serous ovarian [Dorri et al., 2020]) and mRNA transfection time series [Leonhardt et al., 2014], as well as more conventional ones such as mtDNA data and various feature selection/classification datasets. See Figure 6 and references therein for details. See also Appendix I for more information on the models and additional experimental results. All models and algorithms are implemented in the Blang modelling language [Bouchard-Côté et al., 2019].

In each of the 16 models and for each round of Algorithm 4 (line 4), we computed the mean swap acceptance probability across all neighbour chains. We then summarized these means across chains using a box plot. The equi-acceptance objective function of Section 5.1 can be visually understood as collapsing this box plot into a single point. We show in Figure 7 (top) the progression of these swap acceptance probabilities. In all examples considered, equi-acceptance is well approximated within 10 rounds.

Within the range of models considered, we observed a diversity of local barriers $\lambda$ estimated by Algorithm 4 (Figure 7 (middle)). Most statistical models exhibit a high but narrow peak in the neighbourhood of the reference ($\beta=0$). However, a subset of models including statistical models (mixture, ode, phylo-cancer, spike-slab) and physics models (Ising, magnetic, rotor) exhibit additional peaks away from $\beta=0$. See Figure 20 in Appendix J.1 for the corresponding schedule generators $G$ and global barriers $\Lambda$.

7.2 Robustness to ELE violation

To investigate empirically whether the NRPT methodology is robust to the violation of the ELE assumption, we ran NRPT with a range of values for $n_{\textrm{expl}}$ on the models shown in Figure 8. Let $d_{\text{var}}$ denote the number of variables in each model. We run experiments with $n_{\textrm{expl}}=0$, $(1/2)d_{\text{var}}$, $d_{\text{var}}$, $2d_{\text{var}}$, $4d_{\text{var}}$, …, $32d_{\text{var}}$ (the only exception is the reference chain ($\beta=0$), where we always use $n_{\textrm{expl}}=1$ since we can get exact samples from $\pi_{0}$). The key quantity used by the NRPT algorithm for schedule optimization is the communication barrier $\lambda$. The results shown in Figure 8 demonstrate that in all models considered the function $\lambda$ is reliably estimated even when ELE is severely violated, provided $n_{\textrm{expl}}>0$. Similarly, since the estimates of $\Lambda$ and $\bar{\tau}$ are derived from $\lambda$, these quantities can be accurately estimated even when ELE is severely violated (see Figure 19 in Appendix I.4.2). The results shown in Figure 3 support that increasing $n_{\textrm{expl}}$ reduces the difference between the theoretical and observed round trip rate. Moreover, Figure 3 also demonstrate a second strategy to alleviate ELE violation, which is to increase $N$ while fixing $n_{\textrm{expl}}$.

In the next experiment shown in Figure 9, we compare the two mechanisms available for alleviating ELE violation, namely increasing $N$ and increasing $n_{\textrm{expl}}$. First consider the black line in Figure 9 showing the regime where ELE is best approximated in these experiments. This is in close agreement with the theoretical guidelines developed in Section 5.3. At the same time, the light blue lines in Figure 9 show that in situations where the local exploration kernel is not easy to parallelize, it can be advantageous to decrease $n_{\textrm{expl}}$ and to compensate with a number of chains $N$ higher than $2\Lambda$, leading to an average swap rejection rate $r<1/2$.

7.3 Comparison with other parallel tempering schemes

We benchmarked the empirical running time of Algorithm 4 compared to previous adaptive PT methods [Atchadé et al., 2011, Miasojedow et al., 2013]. The methods we considered are: (1) the stochastic optimization adaptive method for reversible schemes proposed in [Atchadé et al., 2011]; (2), a second stochastic optimization scheme, which still selects the optimal number of chains using the $23\%$ rule but uses an improved update scheme from [Miasojedow et al., 2013], refer to Appendix I.3 for details; (3) our non reversible schedule optimization scheme (NRPT); and finally, (4) our scheme, combined with a better initialization based on a preliminary execution of a sequential Monte Carlo algorithm (more precisely, based on a “sequential change of measure,” labelled SCM, as described in [Del Moral et al., 2006]), we use this to investigate the effect on the violation of the stationarity assumption, and for fairness, we use this sophisticated initialization method for all the methods except (3). We benchmarked the methods on four models: (a) a 369-dimensional hierarchical model applied to a dataset of rocket launch failure/success indicator variables [McDowell, 2019]; (b) a 19-dimensional Spike-and-Slab variable selection model applied to the RMS Titanic Passenger Manifest dataset [Hind, 2019]; (c) A 25-dimensional Ising model from Section I.4.2 ($M=5$); (d) a 9-dimensional model for an end-point conditioned Wright-Fisher stochastic differential equation (see, e.g., [Tataru et al., 2017]).

We provide the same number of cores and parallelized all the methods to the greatest extent possible. We refer the reader to Appendix I for implementation details, experimental setup, and detailed description of the models and datasets used.

In Figure 10, each dot summarized in the box plots represents the ESS per total wall clock time in seconds including schedule optimization time for the marginal of one of the model variables. The results show that schedule optimization is more efficient with our proposed non-reversible scheme, with a speed-up in the 10–100 range for the four models considered. The results also suggest that sequential Monte Carlo based initialization may not have a large impact on the performance of our method. See Figure 21 in Appendix J.2 for other summary statistics from this benchmark.

7.4 Mixture models

Bayesian analysis of mixture models can give rise to a label-switching symmetry, leading to a multimodal posterior distribution (see Appendix I.4.4). In Figure 11 we compare the following inference methods on a label-switching posterior distribution: our proposed algorithm (NRPT), an MCMC run based on a single chain (i.e. the exploration kernel alone), the stochastic optimization method of Miasojedow-Moulines-Vihola (MMV), DEO but optimized using MMV (MMV-DEO), Annealed Importance Sampling (AIS) [Neal, 2001], and a sequential Monte Carlo based on a sequence of annealed distributions [Del Moral et al., 2006] (labelled SCM as before). For both SCM and AIS, we use the adaptive scheme of [Zhou et al., 2016] (see Figure 24 in Appendix J.3 for diagnostics of the SCM adaptation). The number of iterations are set so that the method with the smallest wall clock time is our proposed algorithm (NRPT: 1.508min, MMV: 2.019min, MMV-DEO:1.665min, AIS: 1.887min, SCM: 1.778min–5.064min). For all methods, timing includes the time spent to perform schedule optimization.

Amongst MCMC methods, only NRPT correctly captures the multimodality of the target distribution. This is confirmed by the trace plots of the three MCMC methods, shown in Figure 25 in Appendix J.3. MMV automatically selected $N=8$ by targeting an swap acceptance probability of $23\%$. Post-adaptation, the swap acceptance probability was $27\%$. For NRPT we use $N=30$ but to avoid penalizing MMV for a lack of parallelism, we limited all methods to use no more than $8$ threads. After schedule optimization, NRPT estimates the global communication barrier to a value of $\hat{\Lambda}\approx 8$, and the average swap acceptance probability was $72\%$ (see Figure 26 in Appendix J.3 for more NRPT diagnostics).

7.5 Multimodality arising from single cell, whole genome copy number inference

In this section we describe an application to copy number inference in which multimodality arises from the unknown ploidy of a cancer cell. The likelihood of this model is based on a hidden Markov model over $n=6\,206$ observations, and after analytic marginalization of the corresponding $6\,206$ hidden states, the multimodal sampling problem is defined over $d=30$ remaining latent variables. The model is described in Appendix I.4.5, in this section we summarize the results concerning the performance of the algorithms.

We compare the quality of posterior approximations from four samplers on the High-Grade Serous Ovarian cancer dataset from [Dorri et al., 2020]. The four methods compared are: the stochastic optimization adaptive PT method MMV (adaptation selected $7$ chains, average post-adaptation swap acceptance, $28.7\%$), DEO but optimized using MMV (MMV-DEO, adaptation selected $11$ chains, average post-adaptation swap acceptance, $26.9\%$) our proposed algorithm algorithm (NRPT with 25 and 50 chains, average post-optimization swap acceptance of respectively $48.8\%$ and $73.5\%$), and inference based on a single MCMC chain. As in Section 7.4, to avoid overly penalizing MMV for a lack of parallelism, we limited all methods to use no more than $8$ threads, hence obtaining the following wall clock running times including schedule optimization when applicable: MMV, 15.91h; MMV+DEO, 14.33h; 25 chains NRPT, 8.809h; 50 chains NRPT, 15.35h; 1 chain MCMC, 2.91h. The running time of MMV was dominated by adaptation (88% of the wall clock time), and is higher than the schedule optimization time used by NRPT (67%).

Refer to Figure 12 for a comparison of the trace plots between all three methods. Since the experimental protocol provides only proportionality between read counts and copy number, not absolute expected read count for a given copy number, the likelihood does not distinguish between a given copy number profile, and another one obtained by genome duplication of the same profile (refer to Appendix I.4.5 for details). The prior favours the former, but when local copy number events occur after the genome duplication it may be possible to detect recently evolved genome duplication. Challenging multimodality arises when these subsequent events involve only small regions. In such case, the tension between the prior distribution favouring low ploidy and the noisy observation favouring higher ploidy creates a multimodal posterior distribution. We use data from one such cell in the following experiments to illustrate a realistic multimodal inference problem.

The multimodality of the posterior distribution was successfully captured by NRPT but not by a single-chain MCMC nor by the stochastic optimization method MMV (Figure 12). MMV failed to achieve any round trip because the learnt schedule is highly suboptimal: while the mean swap acceptance probability across chains is $28.7\%$, the minimum across chains is nearly zero due to the noise in the optimization procedure. The difficulty of exploring this particular multimodal target is compounded by the fact that the state space of each chromosome’s copy number is random, being upper-bounded by random variables $m_{c}$ corresponding to each chromosome $c$. So in order to perform a jump doubling the cell’s ploidy the sampler has to increase a large number of discrete variables $m_{c}$ simultaneously. This is reflected in the posterior distribution of the model variables denoted $m_{c}$, shown in Figure 28 in Appendix J.4. Hence the local exploration kernel, which in this example samples the variables $m_{c}$ one at the time, is insufficient to jump mode, and only excursions through the prior can achieve mode jumping, via a regeneration based on sampling from the prior at $\beta=0$.

NRPT used 50 chains, estimated $\hat{\Lambda}=13$ (see Figure 30 in Appendix J.4), and converged quickly in this large scale example, performing drastic changes in the annealing schedule in the first 8 rounds then relatively little changes in rounds 8 through the final tenth round (Figure 30).

Comparing the performance of the two NRPT variants, we estimated a round trip rate of $9\times 10^{-4}$ for 25-chain, and $6\times 10^{-3}$ for 50-chain. Notice the increase being more than a factor two supports the use of large $N$ for exploration of complex multimodal posterior distributions. Neither MMV nor MMV-DEO achieved any round trips in the post-burn-in iterations.

8 Discussion

We have established that under the ELE assumption, non-reversible PT dominates reversible PT in terms of round trip rates for any choice of $N$ and annealing schedule $\mathcal{P}_{N}$ (Corollary 1). Moreover we have characterized the optimal round trip rate in terms of the communication barrier $\Lambda$, which we have used to identify an intrinsic limitation of both reversible and non-reversible PT (Section 4.2), and to design a novel schedule optimization method (Section 5). We have shown that this schedule optimization method outperforms previous schedule optimizers in a range of numerical simulations, and we are already aware of other research groups that have successfully applied our methodology in astronomy [Akiyama et al., 2021, Issaoun et al., 2021] and engineering [Langmore et al., 2021].

By combining (28) with the estimate $r^{*}\approx 1/2$ derived from (34), our theoretical results indicate a $1/r^{*}\approx 2$-fold increase in the round trip rate of DEO compared to SEO in an idealized setting where the ELE assumption holds and the schedule, including the number of chains, is well-tuned. The efficiency gains can be more pronounced during schedule optimization, in the course of which the number of chains can be initially larger than the optimal number of chains specified in Section 5.3 as a consequence of Theorem 3. Moreover, our experimental results show that previous methods to tune PT schedules can be fragile. Compared to these previous methods, we have observed speedups of one to two orders of magnitude thanks to our schedule optimization method.

Supplementary material

Appendix A Markov kernel for the index process

For SEO, initialize $(I_{0},\varepsilon_{0})=(i,\varepsilon)$. We then define the Markov transition kernel, $(I_{n+1},\varepsilon_{n+1})|(I_{n},\varepsilon_{n})\sim P^{\mathrm{SEO}}((I_{n},\varepsilon_{n}),\cdot)$ in two steps. In the first step, simulate

[TABLE]

where the expression “$\wedge N\vee 0$” enforces the annealing parameter boundaries. In the second step, independently sample $\varepsilon_{n+1}\sim\text{Unif}\{-1,+1\}$.

Similarly for DEO, initialize $(I_{0},\varepsilon_{0})=(i_{0},\varepsilon_{0})$. Analogous to the SEO construction, we define $(I_{n+1},\varepsilon_{n+1})|(I_{n},\varepsilon_{n})\sim P^{\mathrm{DEO}}((I_{n},\varepsilon_{n}),\cdot)$ in two steps. We first update $I_{n+1}|(I_{n},\varepsilon_{n})=(i,\varepsilon)$ as in (36), but apply the deterministic update in the second step,

[TABLE]

Appendix B Proof of Theorem 1

To simplify notation for the rest of the proof, let $T_{\uparrow}$ and $T_{\downarrow}$ be the hitting times to the target and reference defined by,

[TABLE]

We will also denote

[TABLE]

1. (a)

If we define $a_{\bullet}^{i}={\mathbb{E}}_{\mathrm{SEO}}(T_{\bullet}|I_{0}=i)$ for $i=0,\dots,N$ and $\bullet\in\{\uparrow,\downarrow\}$, then we have

[TABLE]

By the Markov property, for $i=1,\dots,N-1$, $a_{\bullet}^{i}$ satisfies the recursion

[TABLE]

For $i=1,\dots,N$, we substitute in $b_{\bullet}^{i}=a_{\bullet}^{i}-a_{\bullet}^{i-1}$ into (41). After simplification, $b^{i}_{\bullet}$ satisfies the following recursive relation

[TABLE]

The solutions to (42) are

[TABLE]

or equivalently

[TABLE]

We now deal with the case of $\uparrow$ and $\downarrow$ separately.

1. (a)

To determine $a_{\uparrow}^{0}$, we note that if $I_{0}=0$ then $I_{1}=1$ with probability $\frac{1}{2}s_{1}$ and $I_{1}=0$ otherwise. So $a_{\uparrow}^{0}$ satisfies

[TABLE]

or equivalently

[TABLE]

Substituting this into (43) implies $s_{i}b^{i}_{\uparrow}=-2i$. By summing $b^{i}_{\uparrow}=a^{i}_{\uparrow}-a^{i-1}_{\uparrow}$ from $i=1,\dots,N$ and, noting $a^{N}_{\uparrow}=0$, we get

[TABLE] 2. (b)

Similarly to determine $a_{\downarrow}^{N}$ we note that if $I_{0}=N$ then $I_{1}=N-1$ with probability $\frac{1}{2}s_{N}$ and $I_{1}=N$ otherwise. So $a_{\downarrow}^{N}$ satisfies

[TABLE]

or equivalently

[TABLE]

Substituting this into (44) implies $s_{i}b^{i}_{\downarrow}=2+2(N-i)$. By summing $b^{i}_{\downarrow}=a^{i}_{\downarrow}-a^{i-1}_{\downarrow}$ from $i=1,\dots,N$ and, noting $a^{0}_{\downarrow}=0$, we get

[TABLE]

Substituting in (47) and (50) into (40), it follows that

[TABLE] 2. (b)

If we define $a^{i,\varepsilon}_{\bullet}={\mathbb{E}}_{\mathrm{DEO}}(T_{\bullet}|I_{0}=i,\varepsilon_{0}=\varepsilon)$ for $i=0,\dots,N$, $\varepsilon\in\{+,-\}$ and $\bullet\in\{\uparrow,\downarrow\}$, then we have

[TABLE]

Note that for $i=1,\dots,N-1$ $a_{\bullet}^{i,\varepsilon}$ satisfies the recursion relations

[TABLE]

If we substitute $c_{\bullet}^{i}=a_{\bullet}^{i,+}+a_{\bullet}^{i-1,-}$, and $d_{\bullet}^{i}=a_{\bullet}^{i,+}-a_{\bullet}^{i-1,-}$ into (55) and (56) and simplify, we obtain

[TABLE]

By subtracting and adding (57) and (58), we obtain a joint recursion relation for $c_{\bullet}^{i}$ and $d_{\bullet}^{i}$ of the form

[TABLE]

Note that (60) can be rewritten as

[TABLE]

Once one has expressions for $c_{\bullet}^{i}$ and $d_{\bullet}^{i}$, then we can recover $a_{\bullet}^{i,\varepsilon}$ by using

[TABLE]

We now deal with the $\uparrow$ and $\downarrow$ cases separately.

1. (a)

Note that $a_{\uparrow}^{0,-}=a_{\uparrow}^{0,+}+1$. We can substitute this into (57) to get $s_{1}d_{\uparrow}^{1}=-2$, which combined with (61) implies

[TABLE]

Since $a_{\uparrow}^{N,+}=0$ we have $c_{\uparrow}^{N}=-d_{\uparrow}^{N}$, so by summing (59) we get

[TABLE]

After substituting in (64) into (68), we obtain

[TABLE] 2. (b)

Note that $a_{\downarrow}^{N,+}=a_{\uparrow}^{N,-}+1$. We can substitute this expression into (58) to get $s_{N}d_{\downarrow}^{N}=2$, which combined with (61) implies

[TABLE]

Since $a_{\downarrow}^{0,-}=0$ we have $c_{\downarrow}^{1}=d_{\downarrow}^{1}$, so by summing (59) we get

[TABLE]

After substituting in (70) into (74), we obtain

[TABLE]

Finally, by substituting (69) and (75) into (54), it follows that

[TABLE]

Appendix C Proof of Theorem 2

We begin by recalling the following estimate for $r(\beta,\beta^{\prime})$ from Proposition 1 in [Predescu et al., 2004]:

Proposition 1**.**

[Predescu et al., 2004]** Suppose $V^{3}$ is integrable with respect to $\pi_{0}$ and $\pi$. For $\beta\leq\beta^{\prime}$, let $\bar{\beta}=\frac{\beta+\beta^{\prime}}{2}$, then we have

[TABLE]

where $\lambda$ satisfies (22).

When $\beta<\beta^{\prime}$, $(\beta^{\prime}-\beta)\lambda(\bar{\beta})$ is the Riemann sum for $\int_{\beta}^{\beta^{\prime}}\lambda(b)db$ with a single rectangle. Let $C^{k}([0,1])$ be the set of $k$ times continuously differentiable function on $[0,1]$. If $\lambda\in C^{2}([0,1])$, then standard midpoint rule error estimates yield

[TABLE]

Therefore, if $\lambda\in C^{2}([0,1])$ then we can substitute (78) in (77) in Proposition 1, to obtain Theorem 2. This follows from Proposition 2.

Proposition 2**.**

If $V^{k}$ is integrable with respect to $\pi_{0}$ and $\pi$, then $\lambda\in C^{k-1}([0,1])$.

Proof.

Suppose $V^{k}$ is integrable with respect to $\pi_{0}$ and $\pi$, we want to show here that $\lambda:[0,1]\to{\mathbb{R}}_{+}$ given by

[TABLE]

is in $C^{k-1}([0,1])$. If we define $L(x,y)=L(x)L(y)$ and $\pi_{0}(x,y)=\pi_{0}(x)\pi_{0}(y)$, we can rewrite (79) as

[TABLE]

where $\mathcal{Z},g:[0,1]\to{\mathbb{R}}_{+}$ are defined by

[TABLE]

Since $\mathcal{Z}(\beta)>0$ on $[0,1]$, if we can show that $\mathcal{Z},g\in C^{k-1}([0,1])$ then it implies that $\lambda\in C^{k-1}([0,1])$. This is established in Lemma 2. ∎

Lemma 1**.**

If $V^{k}$ is integrable with respect to $\pi_{0}$ and $\pi$ for $k\in{\mathbb{N}}$. Then for all $\beta\in[0,1]$, $j\leq k$, $V^{j}$ is $\pi^{(\beta)}$-integrable.

Proof.

We begin by noting that for $L>0$ and $\beta\in[0,1]$, we have $L^{\beta}\leq 1+L$. This implies

[TABLE]

Therefore since $V^{k}$ is $\pi_{0}$ and $\pi$-integrable, $V^{k}$ is $\pi^{(\beta)}$-integrable. Finally by Jensen’s inequality we have for $j\geq k$,

[TABLE]

∎

Lemma 2**.**

Suppose $V^{k}$ is integrable with respect to $\pi_{0}$ and $\pi$ for some $k\in{\mathbb{N}}$ then:

1. (a)

$\mathcal{Z}\in C^{k}([0,1])$* with derivatives satisfying,*

[TABLE]

for $j\leq k$. 2. (b)

$g\in C^{k-1}([0,1])$* with derivatives satisfying,*

[TABLE]

for $j<k$.

Proof.

1. (a)

Let $h(x,\beta)=L(x)^{\beta}\pi_{0}(x)=\exp(-\beta V(x))\pi_{0}(x)$ which satisfies

[TABLE]

Note for all $\beta\in[0,1]$ and $j\leq k$,

[TABLE]

The left hand side of (93) dominates $\frac{\partial^{j}h}{\partial\beta^{j}}$ uniformly in $\beta$ and is integrable by Lemma 1. The result follows using the Leibniz integration rule. 2. (b)

Let $\tilde{h}(x,y,\beta)=|V(x)-V(y)|L(x,y)^{\beta}\pi_{0}(x,y)$. By using $\log L(x,y)=-V(x)-V(y)$, we get

[TABLE]

Similar to (a), we have for all $\beta\in[0,1]$, $j\leq k-1$,

[TABLE]

The left hand side of (95) dominates $\frac{\partial^{j}\tilde{h}}{\partial\beta^{j}}$ uniformly in $\beta$. It is integrable by Lemma 1 and using the fact that $V^{k}$ is integrable with respect to $\pi_{0}$ and $\pi$. The result follows using the Leibniz integration rule.

∎

Appendix D Proof of Theorem 3

We first note that (b) and (c) follow immediately from (a) and Corollary 1. So to prove Theorem 3 it is sufficient to show (a).

We use the fact that for all $\mathcal{P}_{N}$

[TABLE]

By Theorem 2, we have $\min_{j}s_{j}=1+O(\|\mathcal{P}_{N}\|)$ and, by Corollary 2, we have $\sum_{i=1}^{N}r^{(i-1,i)}=\Lambda+O(N\|\mathcal{P}_{N}\|^{3})$ which combined with (96) implies

[TABLE]

Therefore $E(\mathcal{P}_{N})$ converges to $\Lambda$ at a $O(\|\mathcal{P}_{N}\|)$ rate as $\|\mathcal{P}_{N}\|\to 0$.

Appendix E Multimodal decomposition of communication barrier

Since PT is often used to sample from multimodal targets, it is natural to ask how the communication barrier behaves under the presence of modes. Similar to [Woodard et al., 2009], we partition $\mathcal{X}$ into the disjoint union $\mathcal{X}=\bigcup_{k=1}^{K}\mathcal{X}_{k}$ where we have $\mathcal{X}_{k}$ represents region of the $k$-th mode of $\pi$ and the target decomposes as a mixture of its modes its modes $\pi(x)=\sum_{k=1}^{K}p_{k}\pi(x|\mathcal{X}_{k})$ where $p_{k}=\pi(\mathcal{X}_{k})$ and $\pi(x|\mathcal{X}_{k})=p_{k}^{-1}\pi(x)\mathbb{I}_{\mathcal{X}_{k}}(x)$ are the probability mass and distribution of the $k$-th mode. If we assume that the reference distribution puts the same relative mass on each mode as the target, $\pi_{0}(\mathcal{X}_{k})=\pi(\mathcal{X}_{k})$, then $\pi^{(\beta)}$ decomposes as

[TABLE]

Similarly $V(x)=\sum_{k=1}^{K}V_{k}(x)\mathbb{I}_{\mathcal{X}_{k}}(x)$ where $V_{k}(x)=-\log(\pi(x|\mathcal{X}_{k})/\pi_{0}(x|\mathcal{X}_{k}))$. Define $\lambda_{k,k^{\prime}}(\beta)$ and $\Lambda_{k,k^{\prime}}$ as the local and global communication barrier between mode $k$ and $k^{\prime}$ by

[TABLE]

where $V_{k}^{(\beta)}\stackrel{{\scriptstyle d}}{{=}}V_{k}(X_{k}^{(\beta)})$ for $X_{k}^{\beta}\sim\pi^{(\beta(}(x|\mathcal{X}_{k})$ and $V_{k}^{(\beta)},V_{k^{\prime}}^{(\beta)}$ are independent. An immediate consequence of (22) in Theorem 2 implies the following decomposition of the communication barrier.

Proposition 3** (Multimodal decomposition).**

If $\pi_{0}(\mathcal{X}_{k})=\pi(\mathcal{X}_{k})$ then,

[TABLE]

In particular (100) implies $\Lambda$ is a weighted average over the communication barriers between modes. In the particular case where the modes are exchangeable, we have $\Lambda_{k,k^{\prime}}=\Lambda$ and thus $\Lambda$ is exactly invariant to $K$. We remark that the meta-model used in this section is not realistic for practical problems as we do not know a priori the location of the modes. This meta-model is only used to obtain intuition on the multimodal scalability of the method. Section 7.1 illustrates the empirical behaviour of the method in several genuinely challenging multimodal problems.

Appendix F High-dimensional scaling of communication barrier.

We determine here the asymptotic behaviour of $\lambda$ and $\Lambda$ when the dimension of $\mathcal{X}$ is large. To make the analysis tractable, we assume that $\pi_{d}(x)=\prod_{i=1}^{d}\pi(x_{i})$ as in [Atchadé et al., 2011, Roberts and Rosenthal, 2014]. This provides a model for weakly dependent high-dimensional distributions. We only make this structural assumption on the state space and distribution to establish Proposition 4 below.

The corresponding annealed distributions are thus given by

[TABLE]

Let $\lambda_{d}$ and $\Lambda_{d}$ be the local and global communication barriers for $\pi_{d}$ respectively.

Proposition 4** (High Dimensional Scaling).**

Define $\sigma^{2}(\beta)=\mathrm{Var}(V^{(\beta)})$, for all $\beta\in[0,1]$ we have as $d\to\infty$,

[TABLE]

It follows from Proposition 4 that $\lambda_{d}$ and $\Lambda_{d}$ increase at a $O(d^{1/2})$ rate as $d\to\infty$.

Proof.

For $k=1,2$, let us define $\textbf{V}_{k}^{(\beta)}\stackrel{{\scriptstyle d}}{{=}}\textbf{V}(X_{k}^{(\beta)})$ where $X_{k}^{(\beta)}\sim\pi_{d}^{(\beta)}$ and $\textbf{V}(x)=\sum_{i=1}^{d}V(x_{i})$. The independence structure from Equation (101) tells us that $\textbf{V}_{k}^{(\beta)}$ can be decomposed as $\textbf{V}_{k}^{(\beta)}=\sum_{i=1}^{d}V_{ki}^{(\beta)}$ where $V_{ki}^{(\beta)}$ are iid with a distribution identical to $V^{(\beta)}$, and therefore

[TABLE]

The random variables $\{V_{1i}^{(\beta)}-V_{2i}^{(\beta)}\}_{i=1}^{d}$ are iid with mean zero and variance $2\sigma^{2}(\beta)$. By the central limit theorem,

[TABLE]

Thus we have

[TABLE]

The sequence of variables indexed by $d$ in the expectation in (106) is also uniformly integrable. This follows by noting that the second moment of the integrand in (106) is uniformly bounded in $d$:

[TABLE]

By $d\to\infty$ and using (104) we have

[TABLE]

which proves (102).

To obtain the dimensional scaling limit for $\Lambda_{d}$, we use Jensen’s inequality,

[TABLE]

Finally, (102), (111) along with dominated convergence yield

[TABLE]

∎

Appendix G Scaling limit of index process

G.1 Scaled index process

To establish scaling limits for $(I_{n},\varepsilon_{n})$, it will be convenient to work in a continuous time setting. To do this, we suppose the times that PT iterations occur are distributed according to a Poisson process $\{M(\cdot)\}$ with mean $\mu_{N}$. The number of PT iterations that occur by time $t\geq 0$ satisfies $M(t)\sim\mathrm{Poisson}(\mu_{N}t)$. We define the scaled index process by $Z^{N}(t)=(W^{N}(t),\varepsilon^{N}(t))$ where $W^{N}(t)=I_{M(t)}/N$ and $\varepsilon^{N}(t)=\varepsilon_{M(t)}$ (example realizations in Figure 13).

For convenience, we will denote $\beta_{w}=G(w)$ and use $z=(w,\varepsilon)\in[0,1]\times\{-1,1\}$ to be a scaled index. Define $C({\mathbb{R}}_{+},\mathcal{S})$ and $D({\mathbb{R}}_{+},\mathcal{S})$ to be set of functions $f:{\mathbb{R}}_{+}\to\mathcal{S}$ that are continuous and càdlàg respectively.

The process $Z^{N}\in D({\mathbb{R}}_{+},[0,1]\times\{-1,1\})$ takes values on the discrete set $\mathcal{P}_{N}^{\text{uniform}}\times\{-1,1\}$ and is only well-defined when $Z^{N}(0)=z_{0}\in\mathcal{P}_{N}^{\text{uniform}}\times\{-1,1\}$. To establish convergence, it is useful to extend it to a process $Z^{N}$ which can be initialized at any $z_{0}\in[0,1]\times\{-1,1\}$. Suppose $Z^{N}(0)=z_{0}\in[0,1]\times\{-1,1\}$ and let $T_{1},T_{2},\dots$ be the iteration times generated by the Poisson process $M$. We construct $Z^{N}(t)$ as follows: define $Z^{N}(t)=z_{n}$ for $t\in[T_{n},T_{n+1})$ and update $z_{n+1}|z_{n}$ via a transition kernel dependent on the communication scheme. We determine this transition kernel mirroring the construction from Section 3.5.

Before doing this, it will be useful to define the backward and forward shift operators $\Phi^{N}_{-},\Phi^{N}_{+}:[0,1]\to[0,1]$ by,

[TABLE]

and similarly,

[TABLE]

Intuitively $\Phi^{N}_{\varepsilon}(w)$ represents the location in $[0,1]$ after $w$ moves a distance $\frac{1}{N}$ in the direction of $\varepsilon$ with a reflection at [math] and $1$.

G.2 Scaled index process for reversible PT

Under the reversible communication scheme, if $z_{n}=(w_{n},\varepsilon_{n})\in\{0,1/N,\dots,1\}\times\{-1,1\}$, then we have $w_{n+1}=\Phi^{N}_{\varepsilon_{n}}(w_{n})$ if a swap successfully occurred and $w_{n+1}=w_{n}$ otherwise. In both cases, $\varepsilon_{n+1}\sim\mathrm{Unif}\{-1,+1\}$. Since $\Phi^{N}_{\varepsilon}(w)$ is not only well-defined for $w\in\mathcal{P}_{N}^{\text{uniform}}$ but for $w\in[0,1]$, we naturally extend this construction to any $w\in[0,1]$.

Formally, we generate $(w_{n+1},\varepsilon_{n+1})$ in two steps. In the first step we simulate,

[TABLE]

In the second step we simulate $\varepsilon_{n+1}\sim\mathrm{Unif}\{-1,+1\}$. This defines a continuous time Markov pure jump process $W^{N}\in D({\mathbb{R}}_{+},[0,1])$ with jumps occurring according to an exponential of rate $\mu_{N}$ and is well defined when initialized at any state $w_{0}\in[0,1]$.

From Theorem 19.2 in [Kallenberg, 2002], the infinitesimal generator for $W^{N}$ with reversible communication is

[TABLE]

where the domain $\mathcal{D}(\mathcal{L}_{W^{N}})$ is given by the set of functions such that $\mathcal{L}_{W^{N}}f$ is continuous. Since $\Phi^{N}_{+},\Phi^{N}_{-}$ are continuous, we have $\mathcal{D}(\mathcal{L}_{W^{N}})=C([0,1])$.

G.3 Scaled index process for non-reversible PT

Before defining the transition kernel for the scaled index process under non-reversible communication, it will be convenient to define the propagation function $\Phi^{N}:[0,1]\times\{-1,1\}\to[0,1]\times\{-1,1\}$ for $z=(w,\varepsilon)$,

[TABLE]

and similarly the rejection function $R:[0,1]\times\{-1,1\}\to[0,1]\times\{-1,1\}$,

[TABLE]

Under the non-reversible scheme, if $z_{n}=(w_{n},\varepsilon_{n})\in\mathcal{P}_{N}^{\text{uniform}}\times\{-1,1\}$, then we have $z_{n+1}=\Phi^{N}(z_{n})$ when a swap is accepted and $z_{n+1}=R(z_{n})$ otherwise. Since $\Phi^{N}(z)$ and $R(z)$ are well-defined for all of $z\in[0,1]\times\{-1,1\}$, we naturally extend this construction to any $z\in[0,1]\times\{-1,1\}$.

Formally, we generate $z_{n+1}$ according to the transition kernel,

[TABLE]

This defines a continuous time Markov pure jump process $Z^{N}\in D({\mathbb{R}}_{+},[0,1]\times\{-1,1\})$ with jumps occurring at an exponential of rate $\mu_{N}$. This process is well defined when initialized at any $z_{0}\in[0,1]\times\{-1,1\}$.

Analogously to the reversible case, under non-reversible communication, the infinitesimal generator for $Z^{N}$ is

[TABLE]

where $z=(w,\varepsilon)$ and $\mathcal{D}(\mathcal{L}_{Z^{N}})$ is given by the set of functions $f$ such that $\mathcal{L}_{Z^{N}}f$ is continuous. Since $\Phi^{N}$ has discontinuities at $(\frac{1}{N},-1)$ and $(1-\frac{1}{N},1)$, we can verify that $f\in\mathcal{D}(\mathcal{L}_{Z^{N}})$ if and only if $f(w_{0},-1)=f(w_{0},1)$ for $w_{0}\in\{0,1\}$.

G.4 Weak limits for scaled index processes

Define $W\in C({\mathbb{R}}_{+},[0,1])$ to be the diffusion on $[0,1]$ with generator

[TABLE]

where the domain $\mathcal{D}(\mathcal{L}_{W})$ is the set of functions $f\in C^{2}([0,1])$ such that $f^{\prime}(0)=f^{\prime}(1)=0$. $W$ is a Brownian motion on $[0,1]$ with reflective boundary conditions admitting the uniform distribution $\mathrm{Unif}([0,1])$ as stationary distribution.

Define $Z\in C({\mathbb{R}}_{+},[0,1]\times\{-1,1\})$ to be the PDMP on $[0,1]\times\{-1,1\}$ given by $Z(t)=(W(t),\varepsilon(t))$ where $W(t)$ moves in $[0,1]$ with velocity $\varepsilon(t)$ and the sign of $\varepsilon(t)$ is reversed at the arrivals times of a non-homogeneous Poisson process of rate $\lambda(G(W(t)))G^{\prime}(W(t))$ or when $W(t)$ reaches the boundary $\{(0,-1),(1,+1)\}$; see [Bierkens et al., 2018] for a discussion of PDMP on restricted domains. The infinitesimal generator of $Z$ is given by

[TABLE]

for any $f\in\mathcal{D}(\mathcal{L}_{Z})$, the set of functions $f\in C^{1}([0,1]\times\{-1,1\})$ such that $f(w_{0},-1)=f(w_{0},1)$ and $\frac{\partial f}{\partial w}(w_{0},-1)=-\frac{\partial f}{\partial w}(w_{0},1)$ for $w_{0}\in\{0,1\}$.

G.5 Proof of scaling limit for reversible PT

We will prove Theorem 4(a) by using Theorem 17.25 from [Kallenberg, 2002].

Theorem 5** (Trotter, Sova, Kurtz, Mackevic̆ius).**

Let $X,X^{1},X^{2},\dots$ be Feller processes defined on a state space $S$ with generators $\mathcal{L},\mathcal{L}_{1},\mathcal{L}_{2},\dots$ respectively. If $D$ is a core for $\mathcal{L}$, then the following statements are equivalent:

If $f\in D$, there exist $f_{N}\in\mathcal{D}(\mathcal{L}_{N})$ such that $\|f_{N}-f\|_{\infty}\to 0$ and $\|\mathcal{L}_{N}f_{N}-\mathcal{L}f\|_{\infty}\to 0$ as $N\to\infty$. 2. 2.

If $X^{N}(0)$ converges weakly to $X(0)$ in $S$, then $X^{N}$ converges weakly to $X$ in $D({\mathbb{R}}_{+},S)$.

We will be applying Theorem 5 with $\mathcal{L}=\mathcal{L}_{W}$ defined as $\mathcal{L}_{W}f=\tfrac{1}{2}f^{\prime\prime}$ for $f\in\mathcal{D}(\mathcal{L}_{W})$ where

[TABLE]

and $\mathcal{L}_{N}=\mathcal{L}_{W^{N}}$ defined in (116), which we recall here for the reader’s sake

[TABLE]

with $\Phi_{\pm}^{N}(w)$ defined in (113), (114) and $\beta_{w}=G(w)$. Recall from the discussion just before (116) that $\mathcal{L}_{W^{N}}$ defines a Feller semigroup.

First notice that in [Kallenberg, 2002], the transition semi-group and generator of a Feller process taking values in a metric space $S$ are defined on $C_{0}(S)$, the space of functions vanishing at infinity. Equivalently $f\in C_{0}(S)$ if and only for any $\delta>0$ there exists a compact set $K\subset S$ such that for $x\notin K$, $|f(x)|<\delta$. In our case since $S=[0,1]$ is compact $C_{0}(S)=C(S)$, which justifies the definition of the generator $\mathcal{L}_{W}$ given above.

G.5.1 The Feller property of $\mathcal{L}_{W}$.

Similarly $\mathcal{L}_{W}$ can be seen to define a Feller semigroup on $C([0,1])$ by the Hille-Yosida theorem; see [Kallenberg, 2002, Theorem 19.11].

Indeed the first condition is satisfied since any function $f\in C([0,1])$ can be uniformly approximated within $\epsilon>0$ by a polynomial $p_{\epsilon}$, that is a smooth function, by the Stone-Weirstrass theorem. We can further uniformly approximate $p_{\epsilon}$ within $\epsilon$ by a $C^{2}$ function $\hat{p}_{\epsilon}$ with vanishing derivatives at the endpoints. For example one can let, for a $\delta$ to be chosen later, $\hat{p}_{\epsilon}(x)=p_{\epsilon}(x)$ for $x\in(\delta,1-\delta)$ and for $x\leq\delta$ set $\hat{p}_{\epsilon}(x)=\int_{0}^{x}\rho_{\delta}(y)p_{\epsilon}^{\prime}(y)\mathrm{d}y+c$, where $\rho_{\delta}$ is a smooth, increasing transition function such that $\rho_{\delta}(x)=0$ for $x<0$, $\rho_{\delta}(x)=1$ for $x>\delta$, for example let $\rho_{\delta}=\rho(x/\delta)$, $\rho(x)=g(x)/(g(x)+g(1-x))$ and $g(x)=\exp(-1/x)\mathbf{1}_{\{x>0\}}$. We choose $c$ so that $\hat{p}_{\epsilon}(x)$ is continuous at $\delta$. A similar construction can be used for the right-endpoint. One can then check that indeed $\hat{p}_{\epsilon}\in C^{2}([0,1])$, $\hat{p}_{\epsilon}^{\prime}(0)=\hat{p}_{\epsilon}^{\prime}(1)=0$ and that for $\delta$ small enough $\|\hat{p}_{\epsilon}-p_{\epsilon}\|_{\infty}<\epsilon$.

The second condition of [Kallenberg, 2002, Theorem 19.11] requires that for some $\mu>0$, the set (\mu-\mathcal{L}_{W})(\mathcal{D}\big{(}\mathcal{L}_{W})\big{)} is dense in $C([0,1])$. Let $g\in C([0,1])$ be given. We apply [Saranen and Seikkala, 1988, Corollary 2.2], with $f(t,y)=2\mu y-2g$, which is clearly square integrable in $t$ and $2\mu$-Lipschitz in $y$. Then [Saranen and Seikkala, 1988, Corollary 2.2] implies that for small enough $\mu>0$ the two-point Neumann-boundary value problem

[TABLE]

admits a solution in the Sobolev space $H^{2}([0,1])$ of functions with square integrable first and second derivatives. This already implies that $u\in C^{1}([0,1])$, whereas the continuity of $g$ and of $u$ a priori implies the continuity of $u^{\prime\prime}$ since $u^{\prime\prime}=2\mu u-g$. Overall, for any $g\in C([0,1])$ we can find $u\in\mathcal{D}(\mathcal{L}_{W})$ such that $g=(\mu-\mathcal{L}_{W})g$ establishing the second condition of [Kallenberg, 2002, Theorem 19.11].

The third condition of [Kallenberg, 2002, Theorem 19.11] is that $(\mathcal{L}_{W},\mathcal{D}(\mathcal{L}_{W}))$ satisfies the positive maximum principle, that is if for some $f\in\mathcal{D}(\mathcal{L}_{W}))$ and $x_{0}\in[0,1]$ we have $f(x_{0})\geq f(x)\vee 0$ for all $x\in[0,1]$, then $f^{\prime\prime}(x_{0})\leq 0$. Suppose first that the maximum is attained at an interior point $x_{0}\in(0,1)$; since $f\in C^{2}([0,1])$, by definition of $\mathcal{D}(\mathcal{L}_{W}))$, $f^{\prime\prime}(x_{0})\geq 0$. If on the other hand the positive maximum is attained at $x_{0}=0$, suppose that $f^{\prime\prime}(0)>\epsilon$ for all $x\leq\epsilon$. Thus for $0<y<\epsilon$ small enough, since $f^{\prime}(0)=0$ we have

[TABLE]

thus arriving at a contradiction.

We have thus established that $\left(\mathcal{L}_{W},\mathcal{D}\left(\mathcal{L}_{W}\right)\right)$ satisfies all conditions of [Kallenberg, 2002, Theorem 19.11] and therefore generates a Feller process. Now we can apply Theorem 5 to prove Theorem 4(a). We only need to check the first condition of Theorem 5. In this direction, first note that by definition $\Phi_{\pm}^{N}(w)=w\pm 1/N$ for $w\in[1/N,1-1/N]$. Thus in this case using a Taylor expansion we have for $w^{\ast}_{-}\in[w-1/N,w]$ and $w^{\ast}_{+}\in[w,w+1/N]$ that for any $f\in\mathcal{D}\left(\mathcal{L}_{W}\right)$,

[TABLE]

Since $f^{\prime\prime}$ is uniformly continuous it follows that as $N\to\infty$,

[TABLE]

and therefore for $w\in[1/N,1-1/N]$ we have

[TABLE]

When $w\in[0,1/N)$ or $w\in(1-1/N,1]$ we instead perform a Taylor expansion around 0 or 1 respectively. We only do the calculation in the first case, the other case being similar. Let $w\in[0,1/N)$ in which case, since $f^{\prime}(0)=0$, for $w^{\ast},w^{\ast}_{-},w^{\ast}_{+}\in[0,2/N]$

[TABLE]

where the error term is uniform in $w$ and was obtained by combining the facts that $f^{\prime\prime}$ is uniformly continuous and that $|\Phi^{N}_{\pm}|,|w|/\leq 2/N$. Finally notice that since $w\in[0,1/N]$, then

[TABLE]

Finally we will need the following weaker version of Theorem 1, whose proof is postponed to the end of this section.

Lemma 3**.**

Suppose that $\pi_{0}(|V|),\pi(|V|)<\infty$. Then there exists a constant $C>0$ such that

[TABLE]

Using Lemma 3, we see that for some constant $C>0$

[TABLE]

and therefore

[TABLE]

where the error term is uniform in $w$. Thus $\mathcal{L}_{N}f\to\mathcal{L}f$ uniformly.

of Lemma 3.

Using the bound $0\leq 1-\exp(-x)\leq x$ for $x\geq 0$ we have

[TABLE]

∎

G.6 Proof of scaling limit for non-reversible PT

We will prove Theorem 4(b) in a slightly round about way. We will define the auxiliary processes $\{U^{N}(\cdot)\}$, $\{U(\cdot)\}$ living on the unit circle $\mathds{S}^{1}:=\{z\in\mathds{C}:|z|=1\}$ along with a mapping $\phi:\mathds{S}^{1}\mapsto[0,1]\times\{\pm 1\}$ such that $Z^{N}=\phi(U^{N})$ and $Z=\phi(U)$. We will first show that the law of $U^{N}$ converges weakly to $U$.

Before defining the processes we point out that we will identify $\mathds{S}^{1}$ with $[0,2\pi)$ in the usual way by working in$\mod 2\pi$ arithmetic. Notice that in this way

[TABLE]

The reason for working with these auxiliary processes is that we can now avoid working with PDMPs with boundaries, helping us to remove a layer of technicalities.

For any $N$ we define $\Sigma^{N}:\mathds{S}^{1}\mapsto\mathds{S}^{1}$ through $\Sigma^{N}(\theta)=\theta+2\pi/N$. Consider then a continuous-time process $U^{N}$ that jumps at the arrival times of a homogeneous Poisson process with rate $N$ according to the kernel

[TABLE]

where

[TABLE]

Define the map

[TABLE]

Essentially we think of the circle as comprising of two copies of $[0,1]$ glued together at the end points. The top one is traversed in an increasing direction and the bottom one in a decreasing direction. When glued together and viewed as a circle these dynamics translate in a counter-clockwise rotation with occasional reflections with respect to the $x$-axis at the time of events. With this picture in mind it should be clear that $\phi(U^{N})=Z^{N}$.

We also define the limiting process $U$ as follows. First let

[TABLE]

where $\phi^{1}(\theta)$ is the first coordinate of $\phi(\theta)$. Notice at this point that $\phi^{1}:\mathds{S}^{1}\mapsto[0,1]$ is continuous and satisfies $\phi^{1}(\theta)=\phi^{1}(-\theta)$ for any $\theta\in[0,2\pi)$, whence we obtain that $\tilde{\lambda}(-\theta)=\tilde{\lambda}(\theta)$. Given $U(0)=\theta$, let $T_{1}$ be a random variable such that

[TABLE]

and define the process as $U(s)=\theta+s\mod 2\pi$ for all $s<T_{1}$ and set $U(T_{1})=-U(T_{1}-)\mod 2\pi$. Iterating this procedure will define the $\mathds{S}^{1}$-valued PDMP $\{U(\cdot)\}$. We first need the next lemma.

Lemma 4**.**

Suppose $V$ is integrable with respect to $\pi_{0}$ and $\pi$. The process $U$ defined above is a Feller process, its infinitesimal generator is given by

[TABLE]

with domain

[TABLE]

and invariant measure $\mathrm{d}\theta/2\pi$.

Proof.

First, note that since $\mathds{S}^{1}$ is compact $C_{0}(\mathds{S}^{1})=C(\mathds{S}^{1})$ and thus to study the Feller process we consider the semi-group $\{P_{U}^{t}\}_{t}$ defined by the process $U$ as acting on $C(\mathds{S}^{1})$. To prove the Feller property we can thus use [Davis, 1993, Theorem 27.6]. Since there is no boundary in the definition of $U$ the first assumption is automatically verified, $Qf(\theta)=f(-\theta)\in C\left(\mathds{S}^{1}\right)$ for any continuous $f$. We also know that the rate $\tilde{\lambda}$ is bounded whereas by Proposition 2 and the fact that $G\in C^{1}[0,1]$ we know that $\tilde{\lambda}$ is also continuous. Therefore the third condition of [Davis, 1993, Theorem 27.6] holds and thus $U$ is Feller.

The infinitesimal generator will be defined on $\mathcal{D}(\mathcal{L}_{U})\subseteq C(\mathds{S}^{1})$. The domain is defined as the class of functions $f\in C(\mathds{S}^{1})$ such that

[TABLE]

where the limit is uniform in $\theta$. However by [Böttcher et al., 2013, Theorem 1.33], we can also consider pointwise limits without enlarging the domain. Using the definition of $U$ we then have for $\theta\in[0,2\pi)$ that

[TABLE]

Since for $x\geq 0$ we have $|\exp(-x)-1+x|\leq Cx^{2}$ for some constant $C>0$, and using the uniform continuity of $\tilde{\lambda}$ we can see that

[TABLE]

uniformly in $\theta$, and thus

[TABLE]

In addition

[TABLE]

for any $f\in C(\mathds{S}^{1})$ by strong continuity of $\{P_{U}^{t}\}$ (Feller property) and continuity of $\tilde{\lambda}$. Overall we thus have that $f\in\mathcal{D}(\mathcal{L}_{U})$ if and only if

[TABLE]

which is clearly equivalent to $f\in C^{1}(\mathds{S}^{1})$.

Finally to see that $\mathrm{d}\theta/2\pi$ is invariant, having identified the domain we can easily check that for any $f\in C(\mathds{S}^{1})$ we have

[TABLE]

Since $f\in\mathcal{D}(\mathcal{L}_{U})$ we have that $P_{U}^{s}g\in\mathcal{D}(\mathcal{L}_{U})$. Since for any $g\in\mathcal{D}(\mathcal{L}_{U})$ we have

[TABLE]

∎

Proposition 5**.**

Suppose $U^{N}(0)$ converges weakly to $U(0)$, then $U^{N}$ converges weakly to $U$ in $D({\mathbb{R}}_{+},[0,1])$.

Proof.

We will once again use Theorem 5. The generator of $U_{N}$ is given by

[TABLE]

We will consider the two terms separately. To this end notice that by (23), the boundedness of $\lambda$ and the fact that $G\in C^{1}[0,1]$

[TABLE]

for some $C>0$. Thus, using the mean value theorem, for each $\theta\in[0,2\pi)$, there exists $g_{N}(\theta)\in[\theta,\theta+1/N]$ such that

[TABLE]

where the errors are uniformly bounded and to obtain the second equality above we have used the fact that $|g_{N}(\theta)-\theta|\leq 1/N$ and that $f^{\prime}$ is uniformly continuous, being continuous on a compact set.

Overall we can see that as $N\to\infty$

[TABLE]

Next, using (23) we have that

[TABLE]

where the error is uniform in $\theta$, whence we easily conclude that

[TABLE]

uniformly in $\theta$. ∎

G.6.1 Proof of Theorem 4(b)

Now we are ready to prove the main result of this section. Notice that $Z^{N}(\cdot)=\phi\left(U^{N}(\cdot)\right)$ and $Z(\cdot)=\phi\left(U(\cdot)\right)$.

From Proposition 5 we know that the finite dimensional distributions of $U_{N}$ converge to those of $U$. If $\phi$ were continuous we could conclude using the continuous mapping theorem. Since it is not continuous at the points $\{0,\pi\}$, we will be using [Billingsley, 2013, Theorem 2.7]. We have to check that the law of the limiting process, that is the law of $\{U(\cdot)\}$ places zero mass on finite dimensional distributions that hit $\{0,\pi\}$, that is for $n\in\mathbb{N}$ and $0<t_{1}<\dots<t_{n}$ we want

[TABLE]

when $U(0)$ is initialized according to $\mathrm{d}\theta/2\pi$. But the above follows from the fact that $\mathbb{P}[U(t_{i})\in\{0,\pi\}]=0$, by stationarity when $U(0)$ is initialised uniformly on $\mathds{S}^{1}$.

Relative compactness of $\{Z_{N}(\cdot)\}_{N}$ can be easily seen to follow from the compact containment condition [Ethier and Kurtz, 2009, Remark 3.7.3]. This combined with convergence of the finite dimensional distributions of $Z^{N}$ to those of $Z$ concludes the proof.

Appendix H Distributed implementation of PT

For simplicity, in Section 2 we have shown pseudo-code for PT with a parallel computing context in mind rather than a distributed computing context. The key difference is that in the latter case, it is advantageous to swap annealing parameters instead of states to minimize network communication. For completeness, we show the pseudo-code for the distributed implementation in this section.

Appendix I Experiment supplements

I.1 Reproducibility.

To make our NRPT method easy to use we implemented it as an inference engine in the open source probabilistic programming language (PPL) Blang https://github.com/UBC-Stat-ML/blangSDK. A full description of the models used in the paper are available at https://github.com/UBC-Stat-ML/blangDemos, see in particular https://github.com/UBC-Stat-ML/blangDemos/blob/master/src/main/resources/demos/models.csv for a list of command line options and data paths used for each model. All methods use the same local exploration kernels, namely slice sampling with exponential doubling followed by shrinking [Neal, 2003]. Scripts documenting replication of our experiments are available at https://github.com/UBC-Stat-ML/ptbenchmark.

I.2 Multi-core implementation.

We use lightweight threads [Friesen, 2015] to parallelize both the local exploration and communication phases, as shown in Algorithm 1. We use the algorithm of [Leiserson et al., 2012] as implemented in [Steele and Lea, 2013] to allow each PT chain to have its own random stream. This technique avoids any blocking across threads and hence makes the inner loop of our algorithm embarrassingly parallel in $N$. Moreover, the method of [Leiserson et al., 2012] combined with the fact that we fix random seeds means that the numerical value output by the algorithm is not affected by the number of threads used. Increasing the number of threads simply makes the algorithm run faster. In all experiments unless noted otherwise we use the maximum number of threads available in the host machine, by default an Intel i5 2.7 GHz (which supports 8 threads via hyper-threading) except for Section 7.3 where we use an Amazon EC2 instance of type c4.8xlarge, which is backed by a 2.9 GHz Intel Xeon E5-2666 v3 Processor (20 threads).

I.3 Stochastic optimization methods.

All baseline methods are implemented in Blang (https://github.com/UBC-Stat-ML/blangSDK), the same probabilistic programming language used to implement our method. The code for the baseline adaption methods are available at https://github.com/UBC-Stat-ML/blangDemos. All methods therefore run on the Java Virtual Machine, so their wall clock running times are all comparable.

Both [Atchadé et al., 2011] and [Miasojedow et al., 2013] are based on reversible PT together with two different flavours of stochastic optimization to adaptively select the annealing schedule. In [Atchadé et al., 2011], the chains are added one by one, each chain targeting a swap acceptance rate of $23\%$ from the previous one. In [Miasojedow et al., 2013], this scheme is modified in two ways: first, all annealing parameters are optimized simultaneously, and second, a different update for performing the stochastic optimization is proposed. To optimize all chains simultaneously, the authors assume that both the number of chains and the equi-acceptance probability are specified. Since this information is not provided to the other methods, in order to perform a fair comparison, for the method we label as “Miasojedow, Moulines, Vihola” we implemented a method which adds the chain one at the time while targeting the swap acceptance rate of $23\%$ but based on the improved stochastic optimization update of [Miasojedow et al., 2013]. Specifically, both [Atchadé et al., 2011] and [Miasojedow et al., 2013] rely on updates of the form $\rho_{n+1}=\rho_{n}+\gamma_{n}(\alpha_{n+1}-0.23)$ where $\gamma_{n}$ is an update schedule and $\rho_{n}$ is a re-parameterization of difference in annealing parameter from the previous chain $\beta$ to the one being added $\beta^{\prime}$. The work of [Atchadé et al., 2011] uses the update $\beta^{\prime}_{n}=\beta(1+\exp(\rho_{n}))^{-1}$, whereas the work of [Miasojedow et al., 2013] specifies the explicit parameterization used for $\rho$, namely $\rho=\log(\beta^{\prime-1}-\beta^{-1})$, from which the update becomes $\beta^{\prime}_{n}=\beta(1+\beta\exp(\rho_{n}))^{-1}$. Moreover, while [Atchadé et al., 2011] use $\gamma_{n}=(n+1)^{-1}$, [Miasojedow et al., 2013] suggest to use $\gamma_{n}=(n+1)^{-0.6}$. The experiments confirm that the latter algorithm is more stable.

I.4 Description of models and datasets

In Section 7.3 we benchmark the methods on the following four models. First, a hierarchical model applied to a dataset of rocket launch failure/success indicator variables [McDowell, 2019]. We organized the data by types of launcher, obtaining $5,667$ launches for $367$ types of rockets (processed data available at https://github.com/UBC-Stat-ML/blangDemos/blob/master/data/failure_counts.csv). Each type is associated with a Beta-distributed parameter with parameters tied across rocket types, with the likelihood given by a Binomial distribution (full model specification available at https://github.com/UBC-Stat-ML/blangDemos/blob/master/src/main/java/hier/HierarchicalRockets.bl). The second model is a Spike-and-Slab variable selection model applied to the RMS Titanic Passenger Manifest dataset [Hind, 2019]. The preprocessed data is available at https://github.com/UBC-Stat-ML/blangDemos/tree/master/data/titanic. The data consist in binary classification indicators for the survival of each individual passenger as well as covariates such as age, fare paid, etc. We used a Spike-and-Slab prior with a point mass at zero and a Student-t continuous component (full model specification available at https://github.com/UBC-Stat-ML/blangDemos/blob/master/src/main/java/glms/SpikeSlabClassification.bl). Third, we used the Ising model from Section I.4.2. Finally, we also used an end-point conditioned Wright-Fisher stochastic differential equation (see, e.g., [Tataru et al., 2017]). For this last model we used synthetic data generated by the model. The specification of this last model is available at https://github.com/UBC-Stat-ML/blangSDK/blob/master/src/main/java/blang/validation/internals/fixtures/Diffusion.bl.

The model-specific command line options used for all four experiments is available at https://github.com/UBC-Stat-ML/blangDemos/blob/master/src/main/resources/demos/models.csv. The same file also documents the location of the probabilistic programming source code and precise command line arguments for the 16 models summarized in Figure 7.

I.4.1 Gaussian model

Suppose $\pi\sim N(0,\sigma^{2}\mathbb{I}_{d})$, and $\pi_{0}\sim N(0,\sigma^{2}_{0}\mathbb{I}_{d})$ with $\sigma_{0}>\sigma$. It can be shown that $\pi^{(\beta)}\sim N(0,\sigma^{2}_{\beta}\mathbb{I}_{d})$ where $\sigma_{\beta}^{-2}=(1-\beta)\sigma_{0}^{-2}+\beta\sigma^{-2}$. Theorem 1 in [Predescu et al., 2004] implies the following closed form expressions for $\lambda(\beta)$ and $\Lambda(\beta)$

[TABLE]

where $B(a,b)$ is the Beta function. As $d\to\infty$, we have $\Lambda\sim\sqrt{\frac{2d}{\pi}}\log\left(\frac{\sigma_{0}}{\sigma}\right)$, which is consistent with Proposition 4. We see from Figure 16 that the empirical approximation of $\hat{\lambda},\hat{\Lambda}$ from Algorithm 4 are consistent with the theoretical values from in (168).

Substituting (168) into (30) we have $\beta^{*}_{k}$ satisfies

[TABLE]

This is the same spacing obtained (based on a different theoretical approach) in [Atchadé et al., 2011] and [Predescu et al., 2004] for the Gaussian model. This combined with Proposition 3 also justifies why the geometric schedule works for Gaussian mixture models with well-separated modes. See Figure 15 for realizations of the index process for the reversible and non-reversible PT algorithms.

I.4.2 Ising model

We now compute numerically $\lambda$ for the Ising model on a 2-dimensional lattice of size $M\times M$ with magnetic moment $\mu$. Note that the terminology “Ising” is technically incorrect since we are using a finite $M$ instead of an infinite lattice. The terminology “autologistic” would be more appropriate but we abuse the terminology “Ising” since it is widespread in the computational statistics literature. Similarly, when we talk about “phase transition” this should be interpreted as an approximate phase transition since $M$ is finite.

Using the notation $x_{i}\sim x_{j}$ to indicate sites are nearest neighbours on the lattice, the target distribution is annealed by the inverse temperature $\beta$ and the tempered distributions are given by

[TABLE]

This is an $M^{2}$ dimensional model which undergoes a phase transition as $M\to\infty$ at some critical inverse-temperature $\beta_{c}$. When $\mu=0$ it is known that $\beta_{c}=\log(1+\sqrt{2})/2$ [Baxter, 2007]. We consider experiments with $\mu=0$ and with $\mu=0.1$, the latter denoted “Ising with magnetic field” and abbreviated “magnetic” in composite figures.

We observe that $\lambda$ exhibits very different characteristics in this scenario compared to the Gaussian model: it is not monotonic and is maximized at the critical temperature. Consequently, the optimal annealing schedule is denser near the phase transition. We also note from Figure 17 that both $\lambda$ and $\Lambda$ increases roughly linearly with respect to $M$, similarly to the conclusion of Proposition 4, however here the conditions of Proposition 4 do not apply. As $N$ increases, the round trip rate of reversible PT decays to 0 and non-reversible PT increase towards $(2+2\Lambda)^{-1}$ as seen in Figure 18. This is consistent with Theorem 3.

I.4.3 Discrete-multimodal problem

Consider a discrete state space $\mathcal{X}=\{0,\dots,2k\}$, and let $1_{\mathrm{Even}}:\Omega\to\{0,1\}$ denote the indicator function for even numbers. Define $\pi(x)\propto a^{1_{\mathrm{Even}}(x)}$ for $a>1$ and $\pi_{0}(x)\propto 1$ with $V(x)=-1_{\mathrm{Even}}(x)\log a$. The distribution $\pi$ has $k+1$ modes located where $x$ is even with low probability “barriers” located at $x$ odd. The parameter $a$ controls the relative mass put on the modes. Therefore we have

[TABLE]

where $Z(\beta)=k+(k+1)a^{\beta}$. A simple computation using (22) shows that the local communication barrier is,

[TABLE]

By integrating we obtain the global communication barrier between $\pi$ and $\pi_{0}$,

[TABLE]

I.4.4 Mixture model

We consider a Bayesian mixture model with two mixture components. The likelihood for each component is a normal distribution with a non-conjugate Uniform$(0,100)$ prior on the standard deviation and a normal prior on the means (standard deviation of $100$). We placed a uniform prior on the mixture proportion. We used simulated data generated from the model. While the mixture membership indicator latent random variables can be marginalized in this model, we sample them to make the posterior inference problem more challenging. Sampling these mixture membership random variables is representative of more complex models from the Bayesian non-parametric literature where marginalization of the latent variable is intractable; for example, this is the case for the stick-breaking representation of general completely random measures [Zhu et al., 2020].

In addition to the three MCMC methods described in the main text, we also ran baselines based on SMC and AIS, which are popular methods to explore complex posterior distributions. These methods also depend on the construction of a sequence of annealed distributions from prior to posterior. For the SMC and AIS baselines, to select the sequence of distributions we used an adaptive scheme based on relative ESS as described in [Zhou et al., 2016]. Diagnostics of the adaptation are shown in Figure 26. These methods were parallelized at the particle level. We set the number of particles to achieve a similar running time compared to NRTP, namely $2\,000$ particles for SMC and $2\,500$ particles for AIS. We found that the quality of the posterior approximation was highly dependent on performing several rounds of rejuvenations on the final particle population. The wall-clock time with 5 and 20 rounds of rejuvenation (SMC-5, SMC-20) was comparable to NRPT (1.778min, 2.302min) however the posterior approximation is markedly poor compared to NRPT. With 100 rounds of rejuvenation, the posterior approximation matches closely that of NRPT, however this brings the computation cost to 5.064min. AIS did not perform well since the weights were highly unbalanced in the last iteration, effectively resulting in an approximation putting all mass to a single particle.

I.4.5 Copy number inference

A copy number alteration is a widespread type of mutation in cancers. Whereas in healthy cells each non-sexual chromosomes comes in two copies, in certain cancer cells, chunks of chromosomes of small and large size are frequently deleted or duplicated [Beroukhim et al., 2010]—one extreme type of such event is a whole-genome duplication, in which the copy number of every genomic position is doubled. We consider the task of inferring the copy number for each location in the genome based on single-cell, whole genome data. For simplicity, we consider here a model for performing inference over the copy number profile of one individual cell at a time.

The data consists in the number of reads $y_{i,c}$ measured at different loci. Here $i$ denotes a sub-region or bin of a chromosome $c$. We use the same binning method as [Lai et al., 2020]. The goal is to infer for each bin $i$ and chromosome $c$ an integer copy number $x_{i,c}$. As previously reported in [Zahn et al., 2017], and confirmed in Figure 27, to a first order approximation the location parameter of the variable $y_{i,c}$ is determined by (1) the known GC contents of each bin, denoted $g_{i,c}$ and (2) the unknown copy number $x_{i,c}$. Here instead of performing GC-normalization as a pre-processing step, as done in [Zahn et al., 2017], we perform GC-normalization jointly with copy number inference. We describe in this section how doing so allows us to capture the uncertainty over plausible whole-genome duplication events, and induces multimodal posterior distributions.

The single cell protocol used for the dataset considered here [Zahn et al., 2017] minimizes the PCR amplification biases so that after correcting the GC-bias, the count in a bin is approximately proportional to the copy number of the bin: $\log y_{i,c}=f(\theta,\log g_{i,c})+\log(x_{i,c})+\epsilon_{i,c},$ for $x_{i,c}>0$. Based on Figure 27 in Appendix J.4, we approximate $f$ using a latent quadratic function, and hence $\theta$ is a global parameter containing the 3 coefficients specifying a quadratic polynomial. We assume $\epsilon_{i,c}$ are independent zero-mean normal random variables with standard deviation $\sigma$. To model the prior information that copy number events affect contiguous regions of the genome, for each chromosome $c$, we model the collection of random variables $x_{1,c},x_{2,c},\dots$ using a hidden Markov model (HMM). The maximum copy number in each chromosome is modelled using a geometric random variable $m_{c}$ with a shared latent parameter $p$.

We placed an exponential prior on $\sigma$ with rate $1/10$. We posit normal priors on the coefficients $\theta$ (mean zero, variance 100). We put a uniform prior on $p$. The transition probabilities of the Markov chain are determined by the marginals of a Juke Cantor model with a global rate parameter $\tau$. We put a unit exponential prior on $\tau$. Full model specification is available at https://github.com/UBC-Stat-ML/nowellpack/blob/multimodality-example/src/main/java/chromobreak/SingleCell.bl.

I.5 Multimodality of the examples considered

Figures 4 (bottom left), 22 and 23 support the multimodality of two of the examples considered in Section 7.3, namely the Ising and Spike-and-Slab examples. Multimodality of other examples considered in Section 7.1, 7.2, 7.4 and 7.5 is demonstrated in Figures 4, 27, 25, 28 and 29. Moreover, Figures 4, 25, 28 and 29 demonstrate that using standard MCMC is insufficient to explore these multimodal distributions.

I.6 Selection of $\pi_{0}$

In our experiments all Bayesian models considered have proper priors, and hence we set $\pi_{0}$ to the prior in all these situations. The local exploration kernel in this case consists in an independent draw from the prior (performed by sorting the latent random variables according to a linearization of the partial order induced by the directed graphical model, and sampling their values according to these sorted laws).

For the Ising model, we let $\pi_{0}$ denote a product of independent and identically distributed Bernoulli($1/2$) random variables, one for each node in the Ising grid. It is then straightforward to interpolate between this product distribution and the Ising model of interest using a geometric average. Independent sampling from $\pi_{0}$ is used for the local exploration kernel at $\beta=0$.

Similarly, for the rotor (XY) model, we take $\pi_{0}$ to be a product of independent and identically distributed Uniform($-\pi,\pi$) random variables, and we use a geometric interpolation between $\pi_{0}$ and the target distribution.

Appendix J Supplementary figures for Section 7

J.1 Section 7.1

J.2 Section 7.3

J.3 Section 7.4

J.4 Section 7.5