Weight-Preserving Simulated Tempering

08/14/2018 ∙ by Nicholas G. Tawn, et al. ∙ UNIVERSITY OF TORONTO University of Warwick 0

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).

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

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

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

(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.

1:A temperature schedule ; initialising chain value, ; a within temperature proposal mechanism, ; , the number of algorithm iterations and , the number of within-temperature proposals.
2:
3:function ST()
4:     for  do
5:         
6:          Unif
7:          +w
8:         Compute:
(2)
9:         Sample
10:         if  then
11:              () ()
12:         else
13:              () ()
14:         end if
15:         Perform updates to the -marginal according to to get .
16:     end for
17:     return
18:end function
Algorithm 1 The Simulated Tempering (ST) Algorithm

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

The PT algorithm procedure is given in Algorithm 2.

1:A temperature schedule ; initialising chain values, ; a within temperature proposal mechanism, ; , the number of algorithm iterations and , the number of within-temperature proposals.
2:
3:function PT()
4:     for  do
5:         
6:         Sample uniformly from
7:         Compute:
(3)
8:         Sample
9:         if  then
10:              
11:         else
12:              
13:         end if
14:         for  do
15:               updates to the -marginal chain according to to get .
16:         end for
17:     end for
18:     return
19:end function
Algorithm 2 The Parallel Tempering (PT) Algorithm

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:

(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

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:

(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)):

(6)

where the weights are approximated as

(7)

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.

Figure 1: Power-tempered target densities of a bimodal Gaussian mixture using inverse temperature levels respectively. At the hot state it is evident that the mode centred on 40 begins to dominate the weight as .

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

(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

(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

(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.

Figure 2: For a bimodal Gaussian mixture target, plots of the normalised tempered target densities generated by both powering the target (black) and the WSGM targets (red) at inverse temperature levels .

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

(11)

where is a normalised target density. Then set

(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

  2. [Weight-preserving tempering, version #1]
    Denoting and ; if takes the form

    then .

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

    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.

3.3 Hessian Adjusted Tempering

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

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

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

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

(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

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

(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

(15)

where with

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

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 ,

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.

Figure 3: For the same bimodal Gaussian target from Figure 2, here is a comparison of the HAT vs WSGM targets at inverse temperatures and respectively. Note they are almost identical at the colder temperature; but they do differ slightly in the interval at the hotter temperature where the “G” function is allowing the footprint of the narrow mode to expand.

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

(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:

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.

Figure 4: Top: Trace plots of the functional of the simulated tempering chains given in (18). On the left is the version using the WSGM targets, which mixes well through the temperature schedule and finds both modal regions. On the right is the version using the standard power-based targets, which fails to ever find one of the modes. Bottom: Trace plots of in each of the cases respectively.

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.

(19)

where the skew normal density is given by

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 5: Two trace plots of the first marginal component cold state chain targeting the distribution in (19) using the HAT and PT algorithms respectively. Note the HAT algorithm run illustrates a chain that is performing rapid exploration between all four modes whereas the PT algorithm exhibits significant sticky patches.

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

(20)

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.

Figure 6: Running estimate of , i.e.  given in (20), for 10 runs of the PT (left) and HAT (right) algorithms. The horizontal black line at the level 0.25 represents the true probability that one aims to target in this case. In both cases a burn-in of 10,000 iterations was removed. Observe the lack of convergence of the weight estimates for the PT runs compared to the relatively impressive estimates from the HAT runs.
PT 0.23 0.36 0.19 0.31 0.10 0.12 0.18
0.39 0.51 0
HAT 0.27 0.24 0.26 0.22 0.22 0.27 0.23
0.28 0.25 0.26
Table 1: The end point estimates, , of from the 10 runs of the PT and HAT algorithms. The true value of 0.25 appears to be well approximated by HAT but not by PT.

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.

() RT (secs)
PT 0.288 0.187 0.0593 217
HAT 0.249 0.019 0.0063 451
Table 2: Using the 10 runs of each algorithm in a batch-means approach to estimate the Monte Carlo variance of the pooled estimator i.e. . Also displayed is the average run time (RT, measured in seconds) of a single one of the 10 repaeted runs for both methods respectively.

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

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

(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

(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

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)):

(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

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

then is skew Brownian motion with reflection at

(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

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

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

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

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

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

In a Gaussian setting, up to a proportionality constant

(25)

and at any point

(26)