# Weight-Preserving Simulated Tempering

Simulated tempering is popular method of allowing MCMC algorithms to move between modes of a multimodal target density π. One problem with simulated tempering for multimodal targets is that the weights of the various modes change for different inverse-temperature values, sometimes dramatically so. In this paper, we provide a fix to overcome this problem, by adjusting the mode weights to be preserved (i.e., constant) over different inverse-temperature settings. We then apply simulated tempering algorithms to multimodal targets using our mode weight correction. We present simulations in which our weight-preserving algorithm mixes between modes much more successfully than traditional tempering algorithms. We also prove a diffusion limit for an version of our algorithm, which shows that under appropriate assumptions, our algorithm mixes in time O(d [log d]^2).

## Authors

• 4 publications
• 18 publications
• 14 publications
• ### Optimal Temperature Spacing for Regionally Weight-preserving Tempering

Parallel tempering is popular method for allowing MCMC algorithms to pro...
10/13/2018 ∙ by Nicholas G. Tawn, et al. ∙ 0

• ### Skew Brownian Motion and Complexity of the ALPS Algorithm

Simulated tempering is a popular method of allowing MCMC algorithms to m...
09/25/2020 ∙ by Gareth O. Roberts, et al. ∙ 0

• ### A Framework for Adaptive MCMC Targeting Multimodal Distributions

We propose a new Monte Carlo method for sampling from multimodal distrib...
12/06/2018 ∙ by Emilia Pompe, et al. ∙ 0

• ### Simulated Tempering Method in the Infinite Switch Limit with Adaptive Weight Learning

We investigate the theoretical foundations of the simulated tempering me...
09/13/2018 ∙ by Anton Martinsson, et al. ∙ 0

• ### Does Hamiltonian Monte Carlo mix faster than a random walk on multimodal densities?

Hamiltonian Monte Carlo (HMC) is a very popular and generic collection o...
08/09/2018 ∙ by Oren Mangoubi, et al. ∙ 0

• ### Accelerating Parallel Tempering: Quantile Tempering Algorithm (QuanTA)

Using MCMC to sample from a target distribution, π(x) on a d-dimensional...
08/30/2018 ∙ by Nicholas G. Tawn, et al. ∙ 0

• ### Deep Multimodal Transfer-Learned Regression in Data-Poor Domains

In many real-world applications of deep learning, estimation of a target...
06/16/2020 ∙ by Levi McClenny, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Consider the problem of drawing samples from a target distribution, on a -dimensional state space where

is only known up to a scaling constant. A popular approach is to use Markov chain Monte Carlo (MCMC) which uses a Markov chain that is designed in such a way that the invariant distribution of the chain is

.

However, if exhibits multimodality, then the majority of MCMC algorithms which use tuned localised proposal mechanisms, e.g. Roberts et al. (1997) and Roberts and Rosenthal (2001), fail to explore the state space, which leads to biased samples. Two approaches to overcome this multimodality issue are the simulated and parallel tempering algorithms. These methods augment the state space with auxiliary target distributions that enable the chain to rapidly traverse the entire state space.

The major problem with these auxiliary targets is that in general they don’t preserve regional mass, see Woodard et al. (2009a), Woodard et al. (2009b) and Bhatnagar and Randall (2016). This problem can result in the required run-time of the simulated and parallel tempering algorithms growing exponentially with the dimensionality of the problem.

In this paper, we provide a fix to overcome this problem, by adjusting the mode weights to be preserved (i.e., constant) over different inverse-temperatures. We apply our mode weight correction to produce new simulated and parallel tempering algorithms for multimodal target distributions. We show that, assuming the chain mixes at the hottest temperature, our mode-preserving algorithm will mix well for the original target as well.

This paper is organised as follows. Section 2 reviews the simulated and parallel tempering algorithms and the existing literature for their optimal setup. Section 3 describes the problems with modal weight preservation that are inherent with the traditional approaches to tempering, and introduces a prototype solution called the HAT algorithm that is similar to the parallel tempering algorithm but uses novel auxiliary targets. Section 4 presents some simulation studies of the new algorithms. Section 5 provides a theoretical analysis of a diffusion limit and the resulting computational complexity of the HAT algorithm in high dimensions. Section 6 concludes and provides a discussion of further work.

## 2 Tempering Algorithms

There is an array of methodology available to overcome the issues of multimodality in MCMC, the majority of which use state space augmentation e.g. Wang and Swendsen (1990), Geyer (1991), Marinari and Parisi (1992), Neal (1996), Kou et al. (2006), Nemeth et al. (2017). Auxiliary distributions that allow a Markov chain to explore the entirety of the state space are targeted, and their mixing information is then passed on to aid inter-modal mixing in the desired target. A convenient approach for augmentation methods, such as the popular simulated tempering (ST) and parallel tempering (PT) algorithms introduced in Geyer (1991) and Marinari and Parisi (1992), is to use power-tempered target distributions, for which the target distribution at inverse temperature level is defined as

 πβ(x)∝[π(x)]β

for . For each algorithm one needs to choose a sequence of “inverse temperatures” such that where , so that equals the original target density , and hopefully the hottest distribution is sufficiently flat that it can be easily sampled.

The ST algorithm augments the original state space with a single dimensional component indicating the current inverse temperature level creating a - dimensional chain, , defined on the extended state space that targets

 π(β,x)∝K(β)π(x)β (1)

where ideally , resulting in a uniform marginal distribution over the temperature component of the chain. Techniques to learn exist in the literature, e.g. Wang and Landau (2001) and Atchadé and Liu (2004), but these techniques can be misleading unless used with care. The ST algorithm procedure is given in Algorithm 1.

The PT approach is designed to overcome the issues arising due to the typically unknown marginal normalisation constants. The PT algorithm runs a Markov chain on an augmented state space with target distribution defined as

 πn(x0,x1,…,xn)∝πβ0(x0)πβ1(x1)…πβn(xn).

The PT algorithm procedure is given in Algorithm 2.

### 2.1 Optimal Scaling for the ST and PT Algorithms

Atchadé et al. (2011) and Roberts and Rosenthal (2014) investigated the problem of selecting optimal inverse-temperature spacings for the ST and PT algorithms. Specifically, if a move between two consecutive temperature levels, and , is to be proposed, then what is the optimal choice of

? Too large, and the move will probably be rejected; too small, and the move will accomplish little (similar to the situation for the Metropolis algorithm, cf.

Roberts et al. (1997) and Roberts and Rosenthal (2001)).

For ease of analysis, Atchadé et al. (2011) and Roberts and Rosenthal (2014) restrict to -dimensional target distributions of the iid form:

 π(x)∝d∏i=1f(xi). (4)

They assume that the process mixes immediately (i.e., infinitely quickly) within each temperature, to allow them to concentrate solely on the mixing of the inverse-temperature process itself. To achieve non-degeneracy of the limiting behaviour of the inverse-temperature process as , the spacings are scaled as , i.e.  where a positive value to be chosen optimally.

Under these assumptions, Atchadé et al. (2011) and Roberts and Rosenthal (2014) prove that the inverse-temperature processes of the ST and PT algorithms converge, when speeded up by a factor of , to a specific diffusion limit, independent of dimension, which thus mixes in time , implying that the original ST and PT algorithms mix in time as . They also prove that the mixing times of the ST and PT algorithms are optimised when the value of is chosen to maximise the quantity

 ℓ2×2Φ(−ℓ√I(β)2)

where . Furthermore, this optimal choice of corresponds to an acceptance rate of inverse-temperature moves of 0.234 (to three decimal places), similar to the earlier Metropolis algorithm results of Roberts et al. (1997) and Roberts and Rosenthal (2001).

From a practical perspective, setting up the temperature levels to achieve optimality can be done via a stochastic approximation approach (Robbins and Monro (1951)), similarly to Miasojedow et al. (2013) who use an adaptive MCMC framework (see e.g. Roberts and Rosenthal (2009)).

### 2.2 Torpid Mixing of ST and PT Algorithms

The above optimal scaling results suggest that the mixing time of the ST and PT algorithms through the temperature schedule is , i.e. grows only linearly with the dimension of the problem, which is very promising. However, this optimal, non-degenerate scaling was derived under the assumption of immediate, infinitely fast within-temperature mixing, which is almost certainly violated in any real application. Indeed, this assumption appears to be overly strong once one considers the contrasting results regarding the scalability of the ST algorithm from Woodard et al. (2009a) and Woodard et al. (2009b). Their approach instead relies on a detailed analysis of the spectral gap of the ST Markov chain and how it behaves asymptotically in dimension. They show that in cases where the different modal structures/scalings are distinct, this can lead to mixing times that grow exponentially in dimension, and one can only hope to attain polynomial mixing times in special cases where the modes are all symmetric.

The fundamental issue with the ST/PT approaches are that in cases where the modes are not symmetric, the tempered targets do not preserve the regional/modal weights. That motivates the current work, which is designed to preserve the modal weights even when performing tempering transformations, as we discuss next.

Interestingly, a lack of modal symmetry in the multimodal target will affect essentially all the standard multimodal focused methods: the Equi-Energy Sampler of Kou et al. (2006), the Tempered Transitions of Neal (1996), and the Mode Jumping Proposals of Tjelmeland and Hegstad (2001), all suffer in this setting. Hence, the work in this paper is applicable beyond the immediate setting of the ST/PT approaches.

## 3 Weight Stabilised Tempering

In this section, we present our modifications which preserve the weights of the different modes when performing tempering transformations. We first motivate our algorithm by considering mixtures of Gaussian distributions.

Consider a -dimensional bimodal Gaussian target distribution with means, covariance matrices and weights given by for respectively. Hence the target density is given by:

 π(x)=w1ϕ(x,μ1,Σ1)+w2ϕ(x,μ2,Σ2), (5)

where is the pdf of a multivariate Gaussian with mean and covariance matrix . Assuming the modes are well-separated then the power tempered target from (1) can be approximated by a bimodal Gaussian mixture (cf. Woodard et al. (2009b), Tawn (2017)):

 π(x)=W(1,β)ϕ(x,μ1,Σ1β)+W(2,β)ϕ(x,μ2,Σ2β), (6)

where the weights are approximated as

 W(i,β) = wβi|Σi|1−β2wβ1|Σ1|1−β2+wβ2|Σ2|1−β2 (7) ∝ wβi|Σi|1−β2.

A one-dimensional example of this is illustrated in Figure 1, which shows plots of a bimodal Gaussian mixture distribution as . Clearly the second mode, which was originally wide but very short and hence of low weight, takes on larger and larger weight as , thus distorting the problem and making it very difficult for a tempering algorithm to move from the second mode to the first when is small.

Higher dimensionality makes this weight-distorting issue exponentially worse. Consider the situation with but and where is the identity matrix. Then

 W(2,β)W(1,β)≈σd(1−β), (8)

so the ratio of the weights degenerates exponentially fast in the dimensionality of the problem for a fixed

. This heuristic result in (

8

) shows that between levels there can be a “phase-transition” in the location of probability mass, which becomes critical as dimensionality increases.

### 3.1 Weight Stabilised Gaussian Mixture Targets

Consider targeting a Gaussian mixture given by

 π(x)∝J∑j=1wjϕ(x,μj,Σj) (9)

in the (practically unrealistic) setting where the target is a Gaussian mixture with known parameters, including the weights. By only tempering the variance component of the modes, one can spread out the modes whilst preserving the component weights. We formalise this notion as follows:

###### Definition 1 (Weight-Stabilised Gaussian Mixture (WSGM))

For a Gaussian mixture target distribution , as in (9), the weight-stabilised Gaussian mixture (WSGM) target at inverse temperature level is defined as

 πWSβ(x)∝J∑j=1wjϕ(x,μj,Σjβ). (10)

Figure 2 shows the comparison between the target distributions used when using power-based targets vs the WSGM targets for the example from Figure 1.

Using these WSGM targets in the PT scheme can give substantially better performance than when using the standard power based targets. This is very clearly illustrated in the simulation study section below in Section 4.1. Henceforth, when the term “WSGM ST/PT Algorithm” is used it refers to the implementation of the standard ST/PT algorithm but now using the WSGM targets from (10).

### 3.2 Approximating the WSGM Targets

In practice, the actual target distribution would be non-Gaussian, and only approximated by a Gaussian mixture target. On the other hand, due to the improved performance gained from using the WSGM over just targeting the respective power-tempered mixture, there is motivation to approximate the WSGM in the practical setting where parameters are unknown. To this end, we present a theorem establishing useful equivalent forms of the WSGM; these alternative equivalent forms give rise to a practically applicable approximation to the WSGM.

To establish the notation, let the target be a mixture distribution given by

 π(x)∝J∑j=1hj(x)=J∑j=1wjgj(x) (11)

where is a normalised target density. Then set

 πβ(x)∝J∑j=1fj(x,β)=J∑j=1W(j,β)[gj(x)]β∫[gj(x)]βdx. (12)

Then we have the following result, proved in the Appendix.

###### Theorem 3.1 (WSGM Equivalences)

Consider the setting of (11) and (12) where the mixture components consist of multivariate Gaussian distributions i.e. . Then

1. [Standard, non-weight-preserving tempering]
If then

 W(j,β)∝wβj|Σj|1−β2.
2. [Weight-preserving tempering, version #1]
Denoting and ; if takes the form

 hj(x)exp{(1−β2)(∇j(x))T[∇2j(x)]−1∇j(x)}.

then .

3. [Weight-preserving tempering, version #2]
If

 fj(x,β)=hj(x)βhj(μj)(1−β)

then .

Remark 1: In Theorem 3.1, statement (b) shows that second order gradient information of the ’s can be used to preserve the component weight in this setting.

Remark 2: Statement (c) extends statement (b) to no longer require the gradient information about the but simply the mode/mean point . Essentially this shows that by appropriately rescaling according to the height of the component as the components are “powered up” then component weights are preserved in this setting.

Remark 3: A simple calculation shows that statement (c) holds for a more general mixture setting when all components of the mixture share a common distribution but different location and scale parameters.

The results of Theorem 3.1 are derived under the impractical setting that the components are all known and that is indeed a mixture target. One would like to exploit the results of (b) and (c) from Theorem 3.1 to aid mixing in a practical setting where the target form is unknown but may be well approximated by a mixture.

Suppose herein that the modes of the multimodal target of interest, , are well separated. Thus an approximating mixture of the form given in (11) would approximately satisfy

 π(x)∝hM(x)

where . Hence it is tempting to apply a version of Theorem 3.1(b) to directly as opposed to the . So at inverse temperature , the point-wise target would be proportional to

 π(x)exp{(1−β2)(∇π(x))T[∇2π(x)]−1∇π(x)}.

where and . This only uses point-wise gradient information up to second order. At many locations in the state space, provided that is at a temperature level that is sufficiently cool that the tail overlap is insignificant, and if the target was indeed a Gaussian mixture then this approach would give almost exactly the same evaluations as from (12) in the setting of (b). However, at locations between modes when the Hessian of is positive semi-definite, this target behaves very badly, with explosion points that make it improper.

Instead, under the setting of well separated modes then one can appeal instead to the weight preserving characterisation in Theorem 3.1(c). Assume that one can assign each location in the state space to a “mode point” via some function , with a corresponding tempered target given by

 πβ(x)∝π(x)βπ(μx,β)1−β.

Note the mode assignment function’s dependence on . This can be understood to be necessary by appealing to Figure 2 where it is clear that the narrow mode in the WSGM target has a “basin of attraction” that expands as the temperature increases.

###### Definition 2 (Basic Hessian Adjusted Tempering (BHAT) Target)

For a target distribution on with a corresponding “mode point assigning function” ; the BHAT target at inverse temperature level is defined as

 πBHβ(x)∝π(x)βπ(μx,β)1−β. (13)

However, in this basic form there is an issue with this target distribution at hot temperatures when . The problem is that it leaves discontinuities that can grow exponentially large and this can make the hot state temperature level mixing exponentially slow if using standard MCMC methods for the within temperature moves.

This problem can be mitigated if one assumes more knowledge about the target distribution. Suppose that the mode points are known and so there is a collection of mode points . This assumption seems quite strong but in general if one cannot find mode points then this is essentially saying that one cannot find the basins of attraction and thus the desire to obtain the modal relative masses (as MCMC is trying to do) must be relatively impossible. Indeed, being able to find mode points either prior to or online in the run of the algorithm is possible e.g. Tjelmeland and Hegstad (2001), Behrens (2008) and Tawn et al. (2018). Furthermore, assume that the target, , is in a neighbourhood of the mode locations and so there is an associated collection of positive definite covariance matrices where . From this and knowing the evaluations of at the mode points then one can approximate the weights in the regions to attain a collection where

 ^wj=π(μj)|Σj|1/2∑Kk=1π(μk)|Σk|1/2

With denoting the pdf of a then we define the following modal assignment function motivated by the WSGM:

###### Definition 3 (WSGM mode assignment function)

With collections , and specified above then for a location and inverse temperature define the WSGM mode assignment function as

 A(x,β)=argmaxj{^wjϕ(x|μj,Σjβ)}. (14)

Under the assumption that there are collections , and that have either been found through prior optimisation or through an adaptive online approach we define the following:

###### Definition 4 (Hessian Adjusted Tempering (HAT) Target)

For a target distribution on with collections , and defined above along with the associated mode assignment function given in (14), then the Hessian adjusted tempering (HAT) target is defined as

 πHβ(x)∝{π(x)βπ(μA(x,β))1−βif A(x,β)=A(x,1)G(x,β)if A(x,β)≠A(x,1) (15)

where with

 G(x,β)=π(μ^A)((2π)dΣ^A)1/2ϕ(x|μ^A,Σ^Aβ)βd/2.

Essentially the function “” specifies the target distribution when the chain’s location, , is in a part of the state space where the narrower modes expand their basins of attraction as the temperature gets hotter. Both the choice of and the mode assignment function used in Definition 4 are somewhat canonical to the Gaussian mixture setting. With the same assignment function specified in Definition 3, an alternative and seemingly robust “” that one could use is given by

 G(x,β)=π(x,1,A) +(2P(A(x,β))P(A(x,β))+P(A(x,1))−1)[π(x,β,A)−π(x,1,A)]

where and .

With either of the suggested forms of the function then under the assumption that the target is continuous and bounded on , and that for all ,

 ∫Xπβ(x)dx<∞,

then is a well defined probability density, i.e. Definition 4 makes sense.

For a bimodal Gaussian mixture example Figure 3 compares the HAT target relative to the WSGM target; showing that the HAT targets are a very good approximation to the WSGM targets, even at the hotter temperature levels.

We propose to use the HAT targets in place of the power-based targets for the tempering algorithms given in Section 2. We thus define the following algorithms, which are explored in the following sections.

###### Definition 5 (Hessian Adjusted Simulated Tempering (HAST) Algorithm)

The HAST algorithm is an implementation of the ST algorithm (Section 2, Algorithm 1) where the target distribution at inverse temperature level is given by from Definition 4.

###### Definition 6 (Hessian Adjusted (Parallel) Tempering (HAT) Algorithm)

The HAT algorithm is an implementation of the PT algorithm (Section 2, Algorithm 2) where the target distribution at inverse temperature level is given by from Definition 4.

## 4 Simulation Studies

### 4.1 WSGM Algorithm Simulation Study

We begin by comparing the performances of a ST algorithm targeting both the power-based and WSGM targets for a simple but challenging bimodal Gaussian mixture target example. The example will illustrate that the traditional ST algorithm, using power-based targets, struggles to mix effectively through the temperature levels due to a bottleneck effect caused by the lack of regional weight preservation.

The example considered is the 10-dimensional target distribution given by the bi-modal Gaussian mixture

 π(x)=w1ϕ(μ1,Σ1)(x)+w2ϕ(μ2,Σ2)(x) (16)

where , , , , and . When power based tempering is used, then mode 1 accounts for only 20% of the mass at the cold level, but at the hotter temperature levels becomes the dominant mode containing almost all the mass.

For both runs the same geometric temperature schedule was used:

 Δ={1,0.32,0.322,…,0.326}.

This geometric schedule is justified by Corollary 1 of Tawn and Roberts (2018), which suggests this is an optimal setup in the case of a regionally weight preserved PT setting. Indeed, this schedule induces a swap move acceptance rates around 0.22 for the WSGM algorithm; close to the suggested 0.234 optimal value.

This temperature schedule gave swap acceptance rates of approximately 0.23 between all levels of the power-based ST algorithm except for the coldest level swap where this degenerated to 0.17. That shows that the power-based ST algorithm was set up essentially optimally according to the results in Atchadé et al. (2011).

In order to ensure that the within-mode mixing isn’t influencing the temperature space mixing, a local modal independence sampler was used for the within-mode moves. This essentially means that once a mode has been found, the mixing is infinitely fast. We use the modal assignment function which specifies that the location is in mode 1 if and in mode 2 otherwise. Then the within-move proposal distribution for a move at inverse temperature level is given by

 (17)

where

is the density function of a Gaussian random variable with mean

and variance matrix .

Figure 4 plots a functional of the inverse temperature at each iteration of the algorithm. The functional is

 (18)

where is the usual sign function and is the minimum of the inverse temperatures. The significance of this functional will become evident from the results of the core theoretical contributions made in this paper in Theorems 5.1 and 5.2 in Section 5. Essentially it is taking a transformation of the current inverse temperature at iteration of the Markov chain, such that when the magnitude of is 1 and when the temperature is at its hottest level, i.e. , then is zero. Furthermore, in this example the sign of is a reasonable proxy to identify the mode that the chain is contained in with a negative value suggesting the chain is in the mode centred on and otherwise.

Figure 4 clearly illustrates that the hot state modal weight inconsistency leads the chain down a poor trajectory since at hot temperatures nearly all the mass is in modal region 1. This results in the chain never reaching the other mode in the entire (finite) run of the algorithm. Indeed, the trace plots in Figure 4 show that the chain is effectively trapped in mode 1, which although it only has 20% of the mass in the cold state, is completely dominant at the hotter states.

### 4.2 Simulation study for HAT

To demonstrate the capabilities of the HAT algorithm in a non-Gaussian setting where the modes exhibit skew then a five-dimensional four-mode skew-normal mixture target example is presented. Albeit a mixture, this example can be seen to give similar target distribution geometries to non-mixture settings due to the skew nature of the modes.

 π(x)∝4∑k=1wk5∏i=1f(xi|μk,σk,α) (19)

where the skew normal density is given by

 f(z|μ,σ,α)=2σϕ(z−μσ)Φ(α(z−μ)σ)

and where , , , , , , , , and .

As will be seen in the forthcoming simulation results the imbalance of scales within each modal region ensures that this is a very challenging problem for the PT algorithm.

Since this target fits into the setting of Corollary 1 of Tawn and Roberts (2018) then a geometric inverse temperature schedule is approximately optimal for the HAT target in this setting. Indeed, Tawn and Roberts (2018) suggest that the geometric ratio should be tuned so that the acceptance rate for swap moves between consecutive temperatures is approximately 0.234. In this case, eight tempering levels were used to obtain effective mixing; these were geometrically spaced and given by , was found to be approximately optimal and gave an average of 0.22 for the swaps between consecutive levels for the HAT algorithm.

Using this temperature schedule along with appropriately tuned RWM proposals for the within temperature moves, 10 runs of both the PT and HAT algorithms were performed. In each individual run, each temperature marginal was updated with RWM proposals followed by a temperature swap move proposal and this was repeated with sweeps. This results in a sample output of 600,001 of the cold state chain prior to any burn-in removal. Herein for this example denote .

As expected, the scale imbalance between the modes resulted in the PT algorithm performing poorly and with significant bias in the sample output. In contrast, the HAT approach was highly successful in converging relatively rapidly to the target distribution, exhibiting far more frequent inter-modal jumps at the cold state.

Figure 5 shows one representative example of a run of the PT and HAT algorithms by plotting the first component of the five-dimensional marginal chain at the coldest target state. It illustrates the impressive inter-modal mixing of HAT across all 4 modal regions as opposed to the very sticky mixing exhibited by the PT algorithm.

Figure 6 shows the running approximation of (which is approximately the weight of the first mode i.e. ) after the iteration of the cold state chains, after removing a burn-in period of 10,000 initial iterations, for the ten runs of the PT and HAT runs respectively. The approximation after iteration is given by

 ^Wk1:=1k−10000k∑i=100011(−30

where is the location of the first component of the five-dimensional chain at the coldest temperature level after the

iteration. This figure indicates that the PT algorithm fails to provide a stable estimate for

with the running weight approximations far from stable at the end of the runs; in stark contrast the HAT algorithm exhibits very stable performance in this case. In fact the final estimates for the are given in Table 1.

Table 2 presents the results of using the 10 runs of each algorithm in a batch-means approach to estimate the Monte Carlo variance of the estimator of . The results in Table 2 show that the Monte Carlo error is approximately a factor of 10 higher for the PT algorithm than the HAT approach.

However, it is also important to analyse this inferential gain jointly with the increase in computational cost. Table 2 also shows that the average run time for the 10 HAT runs was 451 seconds which is a little more than 2 times slower than the average run time of the PT algorithm (217 seconds) in this example. The major extra expense is due to the cost of computing the WSGM mode assignment function in (14) at both the cold and current temperature of interest at each evaluation of the HAT target. Anyhow, this is very promising since for a little more than twice the computational cost the inferential accuracy appears to be ten times better in this instance.

## 5 Diffusion Limit and Computational Complexity

In this section, we provide some theoretical analysis for our algorithm. We shall prove in Theorems 5.1 and 5.2 below that as the dimension goes to infinity, a simplified and speeded-up version of our weight-preserving simulated tempering algorithm (i.e., the HAST Algorithm from Definition 5, equivalent to the ST Algorithm 1 with the adjusted target from Definition 4) converges to a certain specific diffusion limit. This limit will allow us to draw some conclusions about the computational complexity of our algorithm.

### 5.1 Assumptions

We assume for simplicity (though see below) that our target density is a mixture of the form (11) with just modes, of weights and respectively, with each mixture component a special i.i.d. product as in (4). We further assume that a weight-preserving transformation (perhaps inspired by Theorem 3.1(b) or (c)) has already been done, so that

 πβ(x) ∝ p[g1(x)]β∫[g1(x)]βdx+(1−p)[g2(x)]β∫[g2(x)]βdx ≡ pgβ1(x)+(1−p)gβ2(x)

for each . We consider a simplified version of the weight-preserving process, in which the chain always mixes immediately within each mode, but the chain can only jump between modes when at the hottest temperature , at which point it jumps to one of the two modes with probabilities and respectively. Let denote the indicator of which mode the process is in, taking value or .

We shall sometimes concentrate on the Exponential Power Family special case in which each of the two mixture component factors is of the form for some . This includes the Gaussian case for which and . (Note that the HAT target in (15) requires the existence of second derivatives about the mode points, corresponding to .)

As in Atchadé et al. (2011) and Roberts and Rosenthal (2014), following Predescu et al. (2004) and Kone and Kofke (2005), we assume that the inverse temperatures are given by , with

 βi=βi−1−ℓ(βi−1)/d1/2 (21)

for some fixed function . In many cases, including the Exponential Power Family case, the optimal choice of is for a constant .

We let be the inverse temperature at time for the -dimensional process. To study weak convergence, we let be a continuous-time version of the process, speeded up by a factor of , where is an independent standard rate 1 Poisson process. To combine the two modes into one single process, we further augment this process by multiplying it by when the algorithm’s state is closer to the second mode, while leaving it positive (unchanged) when state is closer to the first mode. Thus define

 X(d)t=(3−2I)β(d)N(dt). (22)

### 5.2 Main Results

Our first diffusion limit result (proved in the Appendix), following Roberts and Rosenthal (2014), states that when we are at an inverse temperature greater than , the inverse temperature process behaves identically to the case where there is only one mode (i.e. ).

###### Theorem 5.1

Assume the target is of the form (11), with modes of weights and , with inverse weights chosen as in (21). Then up until the first time the process hits , as , converges weakly to a fixed diffusion process given by (22).

Theorem 5.1 described what happens away from . However it says nothing about what happens at . Moreover its statespace is not connected, and we have not even properly defined at . To resolve these issues we define

 h(x)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫xβmin1ℓ(u)du, when % x>0−∫−xβmin1ℓ(u)du, when x<00, when x=0

and set , thus making the process continuous at 0.

###### Remark 1

The process

leaves constant densities locally invariant, for all where is the adjoint of the infinitesimal generator of , as will be shown in the Appendix. This suggests that the density of the invariant distribution of (if it exists) should be piecewise uniform, i.e. it should be constant for and also constant for though these two constants might not be equal.

To make further progress, we require a proportionality condition. Namely, we assume that the quantities corresponding to are proportional to each other in the two modes. More precisely, we extend the definition of to for (corresponding to the first mode), and for (corresponding to the second mode), and assume there is a fixed function and positive constants and such that we have for (in the first mode), while for (in the second mode). For example, it follows from Section 2.4 of Atchadé et al. (2011) that in the Exponential Power Family case, for and for , so that this proportionality condition holds in that case.

Corresponding to this, we choose the inverse temperature spacing function as follows (following Atchadé et al. (2011) and Roberts and Rosenthal (2014)):

 ℓ(β) = I−1/20(β)ℓ0 (23)

for some fixed constant .

To state our next result, we require the notion of skew Brownian motion, a generalisation of usual Brownian motion. Informally, this is a process that behaves just like a Brownian motion, except that the sign of each excursion from 0 is chosen using an independent Bernoulli random variable; for further details and constructions and discussion see e.g. Lejay (2006). We also require the function

 z(h) = h[2Φ(−ℓ02√r(h))]−1/2.

where for and for . We then have the following result (also proved in the Appendix).

###### Theorem 5.2

Under the set-up and assumptions of Theorem 5.1, assuming the above proportionality condition and the choice (23), then as , the process converges weakly in the Skorokhod topology to a limit process . Furthermore, the limit process has the property that if

 Zt=z(h(Xt)),

then is skew Brownian motion with reflection at

 (3−2i)[2Φ(−ℓ02√ri)]−1/2∫1βmin1ℓ(u)du,   i=1,2 . (24)
###### Remark 2

It follows from the proof of Theorem 5.2 that the specific version of skew Brownian motion that arises in the limit is one with excursion weights proportional to

 a=p[2Φ(−ℓ02√r1)]1/2% and\leavevmode\nobreak\ b=(1−p)[2Φ(−ℓ02√r2)]1/2.

That means that the stationary density for on the positive and negative values is proportional to and respectively. This might seem surprising since the limiting weights of the modes should be equal to and , not proportional to and (unless ). The explanation is that the lengths of the positive and negative parts of the domain are given by and respectively. Hence, the total stationary mass of the positive and negative parts – and hence also the limiting modes weights – are still and as they should be.

### 5.3 Complexity Order

Theorems 5.1 and 5.2 have implications for the computational complexity of our algorithm.

In Theorem 5.1, the limiting diffusion process is a fixed process, not depending on dimension except through the value of . It follows that if is kept fixed, then reaches 0 (and hence mixes modes) in time . Since is derived (via ) from the process speeded up by a factor of , it thus follows that for fixed , reaches (and hence mixes modes) in time . So, if is kept fixed, then the mixing time of the weight-preserving tempering algorithm is , which is very fast. However, this does not take into account the dependence on , which might also change as a function of .

Theorem 5.2 allows for control of the dependence of mixing time on the values of . The limiting skew Brownian motion process is a fixed process, not depending on dimension nor on , with range given by the reflection points in (24). It follows that reaches 0 (and hence mixes modes) in time of order the square of the total length of the interval, i.e. of order

 ⎛⎝2∑i=1[2Φ(−ℓ02√ri)]−1/2∫1βmin1ℓ(u)du⎞⎠2

In the Exponential Power Family case, this is easily computed to be .

This raises the question of how large needs to be, as a function of dimension . If the proposal scaling is optimal for within each mode at the cold temperature, then the proposal scaling is . Then, at an inverse temperature , the proposal scaling is . Hence, at an inverse temperature , the probability of jumping from one mode to the other (a distance away) is roughly of order . This is exponentially small unless . This indicates that, for our algorithm to perform well, we need to choose . With this choice, the mixing time order becomes

 ⎛⎝2∑i=1[2Φ(−ℓ02√ri)]−1/2∫11/d21ℓ(u)du⎞⎠2

In the Exponential Power Family case, this corresponds to . That is, for the inverse temperature process to hit and hence mix modes, takes iterations. This is a fairly modest complexity order, and compares very favourably to the exponentially large convergence times which arise for traditional simulated tempering as discussed in Subsection 2.2.

### 5.4 More than Two Modes

Finally, we note that for simplicity, the above analysis was all done for just two modes. However, a similar analysis works more generally. Indeed, suppose now that we have modes, of general weights with . Then when gets to , the process chooses one of the modes with probability . This corresponds to being replaced by a Brownian motion not on , but rather on a “star” shape with different length-1 line segments all meeting at the origin (corresponding, in the original scaling, to ), where each time the Brownian motion hits the origin it chooses one of the line segments with probability each. This process is called Walsh’s Brownian motion, see e.g. Barlow et al. (1989). (The case but corresponds to skew Brownian motion as above.) For this generalised process, a theorem similar to Theorem 5.1 can be then stated and proved by similar methods, leading to the same complexity bound of iterations in the multimodal case as well.

## 6 Conclusion and Further Work

This article has introduced the HAT algorithm to mitigate the lack of regional weight preservation in standard power-based tempered targets. Our simulation studies show promising mixing results, and our theorems indicate the mixing times can become polynomial rather than exponential functions of the dimension , and indeed of time under appropriate assumptions.

Various questions remain to make our HAT approach more practically applicable. The “modal assignment function” needs to be specified in an appropriate way, and more exploration into the robustness of the current assignment mechanism is needed to understand its performance on heavier and lighter tailed distributions. The suggested HAT target assumes knowledge of the mode points which typically one will not have to begin with and one would rely on effective optimisation methods to seek these out either during or prior to the run of the algorithm. Indeed, this has been partially explored by the authors in Tawn et al. (2018). The performance of HAT is heavily reliant on the mixing at the hottest temperature level; the use of RWM here can be problematic for HAT where the mode heights of the disperse modes can be far lower than the narrower modes. As such more advanced sampling schemes such as discretised tempered Langevin could give accelerated mixing at the hot state; the effects of which would be transferred to an improvement in the mixing at the coldest state.

In the theoretical analysis of Section 5, the spacing between consecutive inverse-temperature levels was taken to be to induce a non trivial diffusion limit. However, this result required strong assumptions. Accompanying work in Tawn and Roberts (2018) suggests that for the HAT algorithm under more general conditions, the consecutive optimal spacing should still be , with an associated optimal acceptance rate in the interval .

## 7 Appendix

In this Appendix, we prove the theorems stated in the paper.

### 7.1 Proof of Theorem 3.1

Herein, assume the mixture distribution setting of (11) and (12) where the mixture components consist of multivariate Gaussian distributions i.e. . We prove each of the three parts of Theorem 3.1 in turn.

###### Proof (Proof of Theorem 3.1(a))

Recall that where such that . Hence, taking gives

 fj(x,β) = wβjϕ(x;μj,Σj)β = wβj⎡⎢ ⎢⎣((2π)d|Σj|)1−β2βd/2⎤⎥ ⎥⎦ϕ(x;μj,Σjβ) ∝ wβj|Σj|1−β2ϕ(x;μj,Σjβ)
###### Proof (Proof of Theorem 3.1(b))

Recall the result of Theorem 3.1(a). To adjust for the weight discrepancy from the cold state target a multiplicative adjustment factor, is used such that

 fj(x,β)=hj(x)βαj(x,β)

where . An identical argument to Theorem 3.1(a) shows that this immediately gives .

In a Gaussian setting, up to a proportionality constant

 wj∝hj(x)[(2π)d2|Σj|12exp{12(x−μj)TΣ−1j(x−μj)}] (25)

and at any point

 ∇loghj(x) = −Σ−1j(x−μj) (26) ∇2loghj(x) = −Σ−1j.