1.1. Related works

The DM sampler draws inspiration from Titsias and Dellaportas (Citation2019). In this paper, the authors optimize an objective function composed of the product of the entropy function and the average proposal acceptance rate. The proposals for the adaptive MCMC algorithm are then based off of gradient updates that aim to maximize this function, producing a wide range of proposals that maintain a balance between entropy and acceptance rate. In addition to supporting significant adaptation at early stages of a chain, Gradient-based Adaptive MCMC also allows for adaptation upon rejecting a proposal, a noteworthy feature as most adaptive algorithms do not directly consider the information offered by rejected samples.

While the entropy function is a general function applied to the entire distribution, the algorithm that we present in this paper is based on the premise that, in cases with difficult geometry, it is necessary to focus on specific local regions instead of the entire target distribution when attempting to sample from the target distribution. This is accomplished by leveraging the I-Projection of the target over the set of proposal distributions. The I-projection underestimates the support of the target distribution and will hone in on one area as opposed to the entropy function which attempts to discover a range of samples from the entire target function at once (Shannon Citation1948; Murphy Citation2012). This regional behavior helps to overcome limitations in the Gaussian function class typically used for proposal distributions by reducing the current region of interest into manageable pieces.

Within the broader adaptive MCMC literature, there are a number of proposals that use non-parametric strategies to construct adaptive proposal distributions that asymptotically converge to the target distribution. In Gilks, Best, and Tan (Citation1995), the Adaptive Rejection Metropolis Sampling (ARMS) algorithm iteratively constructs a piecewise linear proposal density in the log-space that approximates the target density. As the sampler progresses, new support points are collected that improve the approximation of the proposal to the target. The piecewise linear approximation can then be used in the Metropolis-Hastings step within a Gibbs sampler. Martino, Read, and Luengo (Citation2015) build on the ARMS algorithm by improving the mechanism by which support points are included in the construction of the proposal density such that all regions of the target support are eligible for inclusion and Meyer, Cai, and Perron (Citation2008) present an alternative specification of ARMS that uses piecewise quadratic functions to construct the proposal density. The ARMS family of algorithms have proven effective at improving the efficiency of Gibbs samplers, particularly in settings where the conditional distributions are not well-known and thus cannot be sampled from directly. To contrast with our own proposal, we note that the ARMS algorithms are designed to construct proposals that approximate univariate targets, whereas our focus is on leveraging the correlation structure found in local regions of a multivariate target distribution.

Parallel tempering is another popular method, and is related to our Scout MCMC algorithm. In parallel tempering, multiple chains are run simultaneously on the target distribution with different levels of tempering applied. The intuition behind parallel tempering is that in the highly tempered chains, it will be easier to cross low-probability boundaries which can subsequently be randomly swapped with the non-tempered chain for mixing (Swendsen and Wang Citation1986; Geyer Citation1991). However, parallel tempering does have its limitations. From a computational perspective, executing many chains but ultimately only using the samples from the non-tempered chain is burdensome.

Using parallel chains for similar purposes, in Craiu, Rosenthal, and Yang (Citation2009), the authors introduce the algorithm Inter-chain Adaptation (INCA), which uses multiple stages of sampling. The first stage involves sampling the state space with parallel chains to partition the state space, while the second stage uses these predetermined regions as a guide to sample from the target distribution. The acceptance probabilities of new proposed points then depend on the region in which the current and proposed points reside.

Finally, the Jumping Adaptive Multimodal Sampler (JAMS) algorithm addresses the challenges of multimodal sampling by front loading the computational burden of mode discovery, using optimization techniques to search for modes and subsequently incorporating this information into the sampling phase (Pompe, Holmes, and Łatuszyński Citation2019). In the sampling phase, dedicated” jump moves” are used to move between modes directly. Once in a mode, any sampler adept at unimodal sampling can be employed. Similar to how we incorporate the DM sampler into parallel tempering to form Scout MCMC, one could envision an algorithm that links the DM sampler to the JAMS algorithm to address both the irregular geometry and the multimodal sampling challenges at once.

2.1. Objective

To find a proposal distribution that minimizes the divergence with the local geometry of the target distribution, we consider using gradient updates performed at each iteration of the MCMC chain. Meanwhile, we must be cognizant of the acceptance rate. In essence, we want to have both a small I-projection so that the proposal and the target are similar, as well as a reasonably high acceptance rate so that we are able to use the samples from our proposals. As such, we propose the following as an objective function that balances both the exponential of the negative I-projection and the average acceptance rate of the proposal: $$s(x)=\hspace{0.17em}\exp \hspace{0.17em}\left[-\beta \mathcal{D}(q||p)\right]·\int \alpha (x,y;L)q(y|x)dy$$

In the above, β is a hyperparameter that balances the impact of the I-projection with the average Metropolis acceptance rate defined by: $$\alpha (x,y;L)=\min \left\{1,\frac{p(y)}{p(x)}\right\}$$ where x is the current position, y is the proposal, and L is the proposal distribution Cholesky factor (Brooks et al. Citation2011). Notice the negative inside the exponential term of s(x). As the I-projection is non-negative, the negative exponent bounds the exponential term between 0 and 1 with the maximum obtained when $\mathcal{D}(q||p)=0.$ Also note that the average acceptance rate ranges between 0 and 1. As a result of these bounds, $s\in [0,1]$ and is maximized when we have high acceptance rates with a proposal that is similar to the target. Thus, the problem of identifying a suitable proposal distribution has been reduced to maximizing s(x) where the optimal proposal distribution at any given x can be characterized by the corresponding optimal Cholesky factor L_x at the global optimum.

To make the objective function easier to manipulate, instead of optimizing s(x), we can optimize the logarithm of s(x). That is: $$\begin{array}{cc}\hspace{0.17em}\log \hspace{0.17em}s(x)\hfill & =-\beta D(q||p)+\hspace{0.17em}\log \hspace{0.17em}\int \alpha (x,y;L)q(y|x)dy\hfill \\ \hfill & =-\beta E_q\left[\log \hspace{0.17em}\frac{q(y|x)}{p(y)}\right]+\hspace{0.17em}\log \hspace{0.17em}{E}_q\left[\alpha (x,y;L)\right]\hfill \\ \hfill & =\beta E_q\left[-\hspace{0.17em}\log \hspace{0.17em}q(y|x)\right]+\beta E_q\left[\log \hspace{0.17em}p(y)\right]+\hspace{0.17em}\log \hspace{0.17em}{E}_q\left[\alpha (x,y;L)\right]\hfill \\ \hfill & =\beta H_q+\beta E_q\left[\log \hspace{0.17em}p(y)\right]+\hspace{0.17em}\log \hspace{0.17em}{E}_q\left[\alpha (x,y;L)\right]\hfill \end{array}$$

The above statement of $\hspace{0.17em}\log \hspace{0.17em}s(x)$ contains expectations entangled with both the p and q distributions that precludes a closed form solution. In particular, notice that the final term is the logarithm of an expectation. Such a term is certainly not ideal for optimization purposes. The most advisable path forward to maximize $\hspace{0.17em}\log \hspace{0.17em}s(x)$ is to instead bound it below using Jensen’s inequality. We can then optimize the lower bound instead of the objective directly. Thus we have: $$\begin{array}{cc}\hspace{0.17em}\log \hspace{0.17em}s(x)\hfill & \ge \beta H_q+\beta E_q\left[\log \hspace{0.17em}p(y)\right]+E_q\left[\log \hspace{0.17em}\alpha (x,y;L)\right]\hfill \\ \hfill & =\beta H_q+\beta E_q\left[\log \hspace{0.17em}p(y)\right]+E_q\left[\log \hspace{0.17em}\min \left\{1,\frac{p(y)}{p(x)}\right\}\right]\hfill \\ \hfill & =\beta H_q+\beta E_q\left[\log \hspace{0.17em}p(y)\right]+E_q\left[\min \left\{0,\hspace{0.17em}\log \hspace{0.17em}p(y)-\hspace{0.17em}\log \hspace{0.17em}p(x)\right\}\right]\hfill \\ \hfill & =\beta H_q+\beta E_{ϵ}\left[\log \hspace{0.17em}p(x+Lϵ)\right]+E_{ϵ}\left[\min \left\{0,\hspace{0.17em}\log \hspace{0.17em}p(x+Lϵ)-\hspace{0.17em}\log \hspace{0.17em}p(x)\right\}\right]\hfill \\ \hfill & =:\mathcal{J}(x)\hfill \end{array}$$

$\mathcal{J}(x)$ can be used as a lower bound for the objective function and for optimization. However, while $\mathcal{J}(x)$ is certainly simpler than $\hspace{0.17em}\log \hspace{0.17em}s(x),$ a general closed form solution of the maximum at each value of x is not attainable. Instead, we turn to iterative optimization methods. We choose gradient ascent as a generally accessible method to maximize $\mathcal{J}(x).$

Gradient ascent requires specifying the gradient of $\mathcal{J}(x)$ with respect to the Cholesky factor L. We leave the detailed derivation of the approximate gradient of $\mathcal{J}(x)$ to Appendix A.1 and present the final result here: $$\begin{array}{cc}{\nabla }_L\mathcal{J}(x)\hfill & =\beta \text{diag}\left(\frac{1}{L_{11}},\dots ,\frac{1}{L_{kk}}\right)+\frac{1}{J}\sum _{j=1}^J\frac{\beta }{p(x+L{ϵ}_j)}p\prime (x+L{ϵ}_j){ϵ}_j^T\hfill \\ \hfill & \hspace{1em}+\frac{1}{J}\sum _{j=1}^J{\nabla }_L\min \left\{0,\hspace{0.17em}\log \hspace{0.17em}p(x+L{ϵ}_j)-\hspace{0.17em}\log \hspace{0.17em}p(x)\right\}.\hfill \end{array}$$ where x is the current position, L is the current value of the Cholesky factor, and ϵ_j are a sample of J standard normal values used to approximate the gradients of the expectations found in $\mathcal{J}.$

Note further that the interior of the second summation in ${\nabla }_L\mathcal{J}(x)$ reduces into the following two cases depending on the value of ϵ_j, $$\begin{array}{cc}\hfill & \hspace{1em}{\nabla }_L\min \left\{0,\hspace{0.17em}\log \hspace{0.17em}p(x+L{ϵ}_j)-\hspace{0.17em}\log \hspace{0.17em}p(x)\right\}\hfill \\ \hfill & =\{\begin{array}{cc}0\hfill & \text{if}\mathrm{ }\hspace{0.17em}\log \hspace{0.17em}p(x+L{ϵ}_j)\ge \hspace{0.17em}\log \hspace{0.17em}p(x)\hfill \\ \frac{1}{p(x+L{ϵ}_j)}p\prime (x+L{ϵ}_j){ϵ}_j^T\hfill & \text{if}\mathrm{ }\hspace{0.17em}\log \hspace{0.17em}p(x+L{ϵ}_j)<\hspace{0.17em}\log \hspace{0.17em}p(x)\hfill \end{array}.\hfill \end{array}$$

The above gradient characterizes the gradient update rule $L_{t+1}=L_t+\gamma {\nabla }_L\mathcal{J}(x)$ where t is the time step of gradient ascent and γ is the step size used to maximize $\mathcal{J}(x).$ Here we make the practical note that due to the presence of the $p{(x+L{ϵ}_j)}^{-1}$ term in the gradient, ϵ_j values that result in proposals with negligible density can cause an explosion of the gradient. We thus set a large threshold value of h to catch elements in the gradient matrix with absolute values greater than h, and set the offending values to $\pm h$ respectively. This event is rare in practice but more common in tail geometries where the fraction of potentially offending proposals is higher.

Now, if we use this procedure to identify a value of L to maximize $\mathcal{J}(x)$ for each point x, call these L_x, we could then characterize a Metropolis-Hastings algorithm using these Cholesky factors.

However, we recognize that a great number of steps would be necessary to optimize $\mathcal{J}(x)$ to within some small error threshold. As gradient updates can be computationally expensive, executing a complete run of gradient ascent at every iteration of an MCMC algorithm would be untenable.

We propose that instead of fully optimizing $\mathcal{J}(x)$ at every iteration, a process that requires many expensive steps, we perform one step of gradient ascent at every MCMC iteration. This will provide approximations of the point-wise optimal sampler discussed so far with the following justifications. First, we note that the early steps of gradient ascent tend to be the most influential and thus a complete run of gradient ascent is not absolutely necessary. Secondly, in practical contexts, changes in geometry are typically gradual which implies that nearby points experience similar behavior, and by extension, similar gradients. While the proposal distribution is not fully optimized at every iteration, on aggregate, the proposal distributions become more optimal as iterations progress.

2.2. Algorithm details

We now gather the results of the above discussions into a complete algorithm summary. The divergence minimization sampler’s objective function and gradient update rule produce a series of covariance matrices for generating Gaussian proposals with an MCMC framework. Consistent with the acceptance rule in the objective function, we incorporate a Metropolis rule for proposal acceptance. At each iteration, we accept the proposal y from the current position x_t with probability: $$\alpha (x_t,y|C_t)=\min \left\{1,\frac{p(y)}{p(x_t)}\right\}$$ where t is the current MCMC iteration, and C_t is the current Cholesky factor of the proposal distribution’s covariance matrix. Note that C_t represents the partially optimized Cholesky factor as opposed to the fully optimized L_x used previously. We reserve discussion of convergence issues for Sec. 2.4.

Algorithm 1.

Divergence Minimization Sampler with Perpetual Adaptation

1: Inputs (defaults): target p(x), balancing parameter β (0.2), initial point x₀, step size γ (0.002), step threshold h (10/γ), initial scaling σ (2), iterations M

2: Initialize: C₀ := $\sigma 1$

3: for t = 0,…,M do

4:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

5:  Propose $y=x_t+C_t{ϵ}_t$

6:  Compute $G={\nabla }_L\mathcal{J}(x_t)$

7:  Accept y with probability $\alpha (x_t,y|C_t)$

8:  Update $x_{t+1}=y$ if accepted or $x_{t+1}=x_t$ if rejected.

9:  If element $\left|G_{ij}\right|>h,$ set $G_{ij}=\text{sign}(G)·h$

10:  Update Cholesky Factor: $C_{t+1}\leftarrow C_t+\gamma G$

11: end for

Algorithm 1 summarizes the DM sampler with perpetual adaptation. In its most basic form, the initial Cholesky factor is set to a diagonal matrix with equal scaling along each dimension though a more complex initialization would also be valid. Furthermore, parameters including the step size and balancing parameters are constant and supplied as inputs although they could be adapted along with the Cholesky factor.

2.3. Divergence minimization: a case study

To understand the behavior of the DM sampler, we examine a case study using a single banana distribution. The banana distribution is a unimodal distribution with non-Gaussian contours. For context, the contours of this distribution are presented in Figure 1. The banana distribution is known to be a difficult distribution to sample from with basic MCMC algorithms because of its quickly changing local geometry, especially in the two tails (Haario, Saksman, and Tamminen Citation1999; Citation2001).

Figure 1. Banana distribution contours. (Note: lighter contours indicate higher density, red dot indicates the origin).

Display full size

An intuitive way to understand the behavior of the DM sampler is to examine its samples. Figure 2 presents parallel results from an adaptive Random Walk Metropolis (aRWM) and DM sampler run. The aRWM algorithm used in this case study and subsequent examples is described by Roberts and Rosenthal (Citation2009). Each algorithm was run for 30,000 iterations with the first 1,000 removed as burn-in. Visually, we notice in the DM sampler results in Figure 2b that the interior of the contours is evenly explored whereas aRWM has blank gaps within the tails of the contours.

Figure 2. Banana distribution samples. The aRWM samples are more sparse and there are gaps in the tails whereas the DM sampler produces more samples in the tails that reach further outwards. (Note: Red dot indicates starting points and blue dots indicate samples).

Display full size

Quantitatively, we compare the algorithms using the acceptance rate and the expected squared jumping distance (ESJD). The ESJD used here balances the goals of a high acceptance rate with the increased exploration of larger steps and is defined as $\mathit{ESJD}={\sum }_{t=2}^M\left|\right|x_t-x_{t-1}|{|}_2^2$ (Roberts and Rosenthal Citation2001; Gelman and Pasarica Citation2007). The DM sampler produced an acceptance rate of 72.25% with an ESJD of 2.4 as compared to the 8.95% acceptance rate and ESJD of 7.3 of aRWM. Since we know that the DM sampler has a higher acceptance rate, this suggests that the DM sampler takes smaller steps and is perhaps less efficient in terms of exploration than aRWM. With that said, more careful steps suggest a lower risk of missing regions of interest.

Such behaviors can be explained by examining the contours of the proposal distribution at different points of the target distribution. Figure 3 presents the contours of the final proposal distribution of the aRWM run centered at the final sample. In other words, they are the contours of the covariance matrix of all samples generated. Notice that the contours have largely failed to adapt to the specific geometry of the target distribution. They have simply expanded so that all regions of meaningful density in the target are covered by the proposal distribution at any given time but have not conformed to the unique geometry of the target distribution (Bédard Citation2007). One might expect that we could achieve similar success with even a simple RWM algorithm, given a large enough proposal distribution.

Figure 3. Proposal distribution contours from the final iteration of aRWM centered at the final sample imposed on the banana distribution contours. Notice that the contours do not match the behavior of the target distribution.

Display full size

We contrast this behavior to that demonstrated by the DM sampler in Figure 4. The DM sampler delivers on the promise of adaptation to local behavior as illustrated by the contours closely matching the region of interest. The proposal distributions benefit from adapted covariance matrices that align with the current tail, resulting in a dramatically reduced likelihood of bad proposals as compared to the aRWM proposals.

Figure 4. Sample DM sampler contours from the same algorithm execution. Notice that the proposal contours in the given iterations conform to the local geometry of the target distribution.

Display full size

In summary, the contour plots demonstrate the intended behavior of the DM sampler to adapt to local regions of interest. Such behavior aligns with the objective of producing desirable adaptation. While it is common for algorithms to simply adapt to the scale of a target distribution, the DM sampler adapts to the behavior of the target distribution, a completely different and much more challenging task that is especially useful for target distributions with unique geometry. Furthermore, while in this instance it seems that aRWM outperforms in efficiency, we must question whether adapting to just the scale of a target distribution is scalable in higher dimensions given that the density will become more and more sparse.

2.4. Convergence and finite adaptation

In a standard adaptive scheme, the algorithm typically involves certain technical conditions (such as diminishing adaptation and containment, or finite adaptation) to guarantee convergence to the target distribution (Roberts and Rosenthal Citation2007; Rosenthal Citation2011). In this work, we argue that certain target distributions, such as those with unusual geometry or those with many unique modes, lend themselves to perpetual adaptation as no single Gaussian proposal distribution could hope to sample well in all regions of interest. The banana example in the previous section is a good example to illustrate this. As the banana distribution is clearly non-Gaussian, a single non-adapting Gaussian proposal cannot appropriately orient itself in the apex and in both tails. However, the consequence of embracing perpetual adaptation is that the standard convergence framework for adaptive MCMC is no longer compatible.

Our goal is to sample efficiently using region-specific proposal distributions while still fulfilling the requirements of the standard convergence framework. As such, we propose a two phase approach that limits adaptation to a finite number of iterations and subsequently transfers the lessons learned in adaptation to a Metropolis-Hastings framework. By limiting adaptation to a finite number of iterations, convergence of the non-adaptive phase to the target distribution is guaranteed (Roberts and Rosenthal Citation2007; Rosenthal Citation2011).

Recall that the basis of the DM sampler is to approximate the optimal proposal distribution characterized by the Cholesky factor L_x that maximizes the objective function s(x). As discussed in Sec. 2.1, if we knew the values of all L_x, we could produce a simple Metropolis-Hastings algorithm with defined proposal distributions. Of course, we have seen that optimizing s(x) is difficult for a single point, let alone all points in space. Fortunately, the procedure described in Algorithm 1 constructs Cholesky factors C_t at each iteration t to approximate the given location’s optimal proposal structure for sampling. If we record these Cholesky factors after each iteration, they can act as a proxy for the optimal Cholesky factor for nearby points as well, assuming some degree of continuity. Thus, after generating a collection of points and their associated Cholesky factors in the adaptive phase, in each iteration of the non-adaptive phase we select the Cholesky factor associated with the closest adaptive phase sample to construct a proposal covariance matrix for the current iteration. The algorithm thus proposes points from the distribution $q(y|x_t)\sim N(x_t,C_t{C}_t^T)$ and accept with the following rule: $${\alpha }_f(x_t,y|C_t,{\tilde{C}}_y)=\min \left\{1,\frac{p(y)q(x_t|y)}{p(x_t)q(y|x_t)}\right\}$$ where x_t is the current position, y is the proposal, $q(y|x_t)\sim N(x_t,C_t{C}_t^T),q(x_t|y)\sim N(y,{\tilde{C}}_y{\tilde{C}}_y^T),$ and C_t and ${\tilde{C}}_y$ are the Cholesky factors from the adaptive phase iterations that correspond to the points closest to x_t and y respectively. In other words, instead of calculating a new Cholesky factor for every new point, we select the point from our adaptive phase that is closest to the new point and use its corresponding (approximate) Cholesky factor. This non-adaptive phase adheres to the standard validity criteria of a non-adaptive Metropolis-Hastings algorithm. A complete algorithm summary of this scheme is presented in Algorithm 2.

We test the finite adaptation variant of the DM sampler on the banana distribution presented in Sec. 2.3. Essentially, once the adaptive phase completes, we consolidate the samples and covariance matrices and then begin the non-adaptive phase at the last adaptive phase sample. In Figure 5, we present 20,000 samples generated from the non-adaptive phase. In addition to the visual indication of the non-adaptive samples covering the relevant portions of the state space, we note that the acceptance rate is 55.93%, the proportion of samples in the left side of the distribution is 51%, and the sample mean of $\left[\begin{array}{cc}0.31& -8.04\end{array}\right]$ is approaching the true mean. These diagnostics indicate the algorithm is sampling well and is converging to the target distribution as expected.

Figure 5. Finite adaptation DM sampler variant. Samples presented only include those from the non-adaptive phase.

Display full size

We note that the finite adaptation version of the DM sampler performs similarly to the perpetually adapting version with the added benefit of adhering to established convergence criteria.

Algorithm 2.

Divergence Minimization Sampler with Finite Adaptation

1: Inputs (defaults): target p(x), balancing parameter β (0.2), initial point x₀, step size γ (0.002), step threshold h (10/γ), initial scaling σ (2), iterations M, finite adaptation threshold F (M/2), finite subsample size s (M/20)

2: Initialize: C₀: = $\sigma 1$

3: Adaptive Phase

4: for t = 0,…,F do

5:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

6:  Propose $y=x_t+C_t{ϵ}_t$

7:  Compute $G={\nabla }_L\mathcal{J}(x_t)$

8:  Accept y with probability $\alpha (x_t,y|C_t)$

9:  Update $x_{t+1}=y$ if accepted or $x_{t+1}=x_t$ if rejected.

10:  If element $\left|G_{ij}\right|>h,$ set $G_{ij}=\text{sign}(G)·h$

11:  Update Cholesky Factor: $C_{t+1}\leftarrow C_t+\gamma G$

12: end for

13: Let S be a sample of s points from $0,1,\dots F$

14: Non-Adaptive Phase

15: for t = F + 1,…,M do

16:  Select $C_t:=C_i$ where $i\in S,$ such that $d(x_i,x_t)$ is minimized.

17:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

18:  Propose $y_t=x_t+C_t{ϵ}_t$

19:  Select ${\tilde{C}}_y:=C_j$ where $j\in S,$ such that $d(x_j,y)$ is minimized.

20:  Accept y with probability ${\alpha }_f(x_t,y|C_t,{\tilde{C}}_y)$

21:  Update $x_{t+1}=y$ if accepted or $x_{t+1}=x_t$ if rejected.

22: end for

Remark.

We now comment on the choice of the Metropolis acceptance rule for the original DM sampler as well as discuss an alternative that could perhaps motivate future work. Recall that each iteration of the adaptive phase triggers the gradient update rule. Any proposal under this framework will be asymmetric which at first glance would suggest the use of a Metropolis-Hastings acceptance rule (Hastings Citation1970). Suppose for a moment that we were to consider a Metropolis-Hastings rule. In other words, we replace the acceptance rule α with the following: $${\alpha }^{*}(x_t,y|C_t)=\min \left\{1,\frac{p(y)q(x_t|y)}{p(x_t)q(y|x_t)}\right\}$$ where x_t is the current position, y is the proposal, C_t is the Cholesky factor of the proposal covariance matrix, $q(y|x_t)\sim N(x_t,C_t{C}_t^T),$ and $q(x_t|y)\sim N(y,(C_t+\gamma {\nabla }_L\mathcal{J}(x_t)){(C_t+\gamma {\nabla }_L\mathcal{J}(x_t))}^T).$ The distribution of $q(x_t|y)$ considers the gradient step made in the process of moving from x_t to y, reflecting the asymmetry involved in returning from y to x_t. In this case, reversibility would be upheld at each individual iteration without introducing any finite adaptation (Roberts and Smith Citation1994; Bai, Roberts, and Rosenthal Citation2011; Craiu et al. Citation2015). However, the proposal kernels across iterations are not necessarily identical. Concretely, visiting, leaving, and then returning to a point can result in different proposal kernels at the same point due to the use of only a single gradient step at each iteration. Thus, each individual step under the hypothetical Metropolis-Hastings setup would be reversible but, in aggregate, the entire chain may not be. This perpetual adaptation represents a departure from the established convergence theory. In this paper, we have instead decided to proceed with a finite adaptation scheme that does guarantee convergence to the target distribution.

3.1. Algorithm details

Scout MCMC generates proposals for the main chain $q(y_t|x_t)\sim N(x_t,C_t{C}_t^T),$ and accepts with probability: $$\begin{array}{cc}\alpha (x_t,y_t|C_t)\hfill & =\min \left\{1,\frac{p(y_t)}{p(x_t)}\right\}\hfill \end{array}$$ Then, it adapts the main chain Cholesky factor by $\gamma {\nabla }_L\mathcal{J}(x_t)$ as before. Next, a proposal for the scout chain is generated as $q(c_t|s_t)\sim N(s_t,{\sigma }_{s}1),$ which is accepted according to the following rule: $$\begin{array}{cc}{\alpha }_s(s_t,c_t)\hfill & =\min \left\{1,\frac{p{(c_t)}^{\tau }}{p{(s_t)}^{\tau }}\right\}\hfill \end{array}$$

Finally, Scout MCMC considers swapping $x_{t+1}$ and $s_{t+1}$ every k iterations according to the swap rule: $$\begin{array}{cc}{\alpha }_{\text{swap}}(x,s)\hfill & =\min \left\{1,\frac{p{(x)}^{\tau }p(s)}{p(x)p{(s)}^{\tau }}\right\}\hfill \end{array}$$

Algorithm 3.

Scout MCMC with Perpetual Adaptation

1: Inputs (defaults): target p(x), balancing parameter β (0.2), temperature τ (0.1), initial point x₀, step size γ (0.002), step threshold h (10/γ), initial scaling σ (2), tempered scaling σ_s (9), iterations M, swap frequency k (20)

2: Initialize: C₀: = $\sigma 1,$ s₀ = x₀

3: for t = 0,…,M do

4:  Main Chain Step

5:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

6:  Propose $y_t=x_t+C_t{ϵ}_t$

7:  Compute $G={\nabla }_L\mathcal{J}(x)$

8:  Accept y_t with probability $\alpha (x_t,y_t|C_t)$

9:  Update $x_{t+1}=y_t$ if accepted or $x_{t+1}=x_t$ if rejected.

10:  If element $\left|G_{ij}\right|>h,$ set $G_{ij}=\text{sign}(G)·h$

11:  Update Cholesky Factor: $C_{t+1}\leftarrow C_t+\gamma G$

12:  Scout Step

13:  Propose $c_t\sim N(s_t,{\sigma }_{s}1)$ and accept with probability $\min \left\{1,\frac{p{(c_t)}^{\tau }}{p{(s_t)}^{\tau }}\right\}$

14:  Update $s_{t+1}=c_t$ if accepted or $s_{t+1}=s_t$ if rejected.

15:  Swap Step

16:  if $t\equiv 0\hspace{0.5em}\mathrm{mod}k$ then

17:   Swap $x_{t+1}$ and $s_{t+1}$ with probability $\min \left\{1,\frac{p(s_{t+1})p{(x_{t+1})}^{\tau }}{p(x_{t+1})p{(s_{t+1})}^{\tau }}\right\}$

18:  end if

19: end for

Algorithm 3 provides pseudocode for the implementation of Scout MCMC. Once again, control over step size and initial scaling is determined by the user to allow flexibility between targets. For example, depending on the expected global region of interest, the tempered chain scaling can be adjusted. Additional details can be added such as adapting the scaling of the tempered chain or varying the limit on the frequency of swap moves.

3.2. Finite adaptation

Similar to the DM sampler, we present a two-phase finitely adapting variant of Scout MCMC. The first phase is the procedure presented in Algorithm 3. In the second, non-adaptive phase, the structure of the scout chain does not change. However, the main chain follows the same process as the finitely adapting DM sampler where the Cholesky factor corresponding to the nearest iteration of the adapting phase is used to construct proposal distributions in the non-adaptive phase. This reduces the non-adapting phase to a Metropolis-Hastings algorithm. We present the pseudocode associated with the finitely adapting Scout MCMC in Algorithm 4.

Algorithm 4.

Scout MCMC with Finite Adaptation

1: Inputs (defaults): target p(x), balancing parameter β (0.2), temperature τ (0.1), initial point x₀, step size γ (0.002), step threshold h (10/γ), initial scaling σ (2), tempered scaling σ_s (9), iterations M, finite adaptation threshold F (M/2), swap frequency k (20), finite subsample size s (M/20)

2: Initialize: C₀: = $\sigma 1,$ s₀ = x₀

3: Adaptive Phase

4: for t = 0,…,F do

5:  Main Chain Step

6:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

7:  Propose $y_t=x_t+C_t{ϵ}_t$

8:  Compute $G={\nabla }_L\mathcal{J}(x)$

9:  Accept y_t with probability $\alpha (x_t,y_t|C_t)$

10:  Update $x_{t+1}=y_t$ if accepted or $x_{t+1}=x_t$ if rejected.

11:  If element $\left|G_{ij}\right|>h,$ set $G_{ij}=\text{sign}(G)·h$

12:  Update Cholesky Factor: $C_{t+1}\leftarrow C_t+\gamma G$

13:  Scout Step

14:  Propose $c_t\sim N(s_t,{\sigma }_{s}1)$

15:  Accept c_t with probability $\min \left\{1,\frac{p{(c_t)}^{\tau }}{p{(s_t)}^{\tau }}\right\}$

16:  Update $s_{t+1}=c_t$ if accepted or $s_{t+1}=s_t$ if rejected.

17:  Swap Step

18:  if $t\equiv 0\hspace{0.5em}\mathrm{mod}k$ then

19:   Swap $x_{t+1}$ and $s_{t+1}$ with probability $\min \left\{1,\frac{p(s_{t+1})p{(x_{t+1})}^{\tau }}{p(x_{t+1})p{(s_{t+1})}^{\tau }}\right\}$

20:  end if

21: end for

22: Let S be a sample of s points from $0,1,\dots F$

23: Non-Adaptive Phase

24: for t = F + 1,…,M do

25:  Main Chain Step

26:  Select $C_t:=C_i$ where $i\in S,$ such that $d(x_i,x_t)$ is minimized.

27:  Generate ${ϵ}_t\sim N(\mathbf{0},\mathrm{ }1)$

28:  Propose $y_t=x_t+C_t{ϵ}_t$

29:  Select ${\tilde{C}}_t:=C_j$ where $j\in S,$ such that $d(x_j,y)$ is minimized.

30:  Accept y with probability ${\alpha }_f(x_t,y|C_t,{\tilde{C}}_t)$

31:  Update $x_{t+1}=y$ if accepted or $x_{t+1}=x_t$ if rejected.

32:  Scout Step

33:  Propose $c_t\sim N(s_t,{\sigma }_{s}1)$ and accept with probability $\min \left\{1,\frac{p{(c_t)}^{\tau }}{p{(s_t)}^{\tau }}\right\}$

34:  Update $s_{t+1}=c_t$ if accepted or $s_{t+1}=s_t$ if rejected.

35:  Swap Step

36:  if $t\equiv 0\hspace{0.5em}\mathrm{mod}k$ then

37:   Swap $x_{t+1}$ and $s_{t+1}$ with probability $\min \left\{1,\frac{p(s_{t+1})p{(x_{t+1})}^{\tau }}{p(x_{t+1})p{(s_{t+1})}^{\tau }}\right\}$

38:  end if

39: end for

In the next section, we present examples demonstrating that similar to the DM sampler, the finite adaptation version of Scout MCMC performs similarly in practice to the perpetually adapting version but has the theoretical advantage of adhering to established convergence criteria.

4.1. Double banana distribution

The first distribution we consider is an extension of the banana distribution examined in Sec. 2.3. Specifically, we consider a pair of banana distributions with overlap in the tails. This results in two primary modes along the curves of the two bananas as well as two secondary modes at the intersections. Figure 6 provides the contours of this distribution.

Figure 6. Double banana distribution contours. (Note: lighter contours indicate higher density, red dot indicates the origin).

Display full size

All of the algorithms are run for 50,000 iterations with the first 1,000 samples discarded as burn-in. Table 1 presents the results of this experiment. For this specific distribution, the target mean is $\left[\begin{array}{cc}0& -25\end{array}\right].$

Notice that the samples of aRWM and Scout MCMC are closest to the mean of the distribution. It is worth noting that there is negligible density at the mean as illustrated by Figure 6 so the ability to achieve the correct mean indicates that both bananas have been visited. In comparison, standard RWM, MALA, and parallel tempering have not achieved the level of success of the other algorithms. Even with 5 chains, parallel tempering has largely failed to converge within the 50,000 iterations. The distinction between aRWM and Scout MCMC lies in the efficiency diagnostics. We see that Scout MCMC accepts over 10 times as many proposals as aRWM though it has a smaller ESJD. A large acceptance rate is not necessarily indicative of a better algorithm but with both algorithms performing similarly, this could indicate that Scout MCMC produces higher quality proposals. Since ESJD is a measure of both the acceptance rate and the step size made with each move but Scout MCMC has a much higher acceptance rate, this would indicate that aRWM proposes moves with much greater step sizes than Scout MCMC. This is expected, however, as Scout MCMC uses a user-specified cooldown period where the main chain makes local moves and does not swap with the scout chain. Finally, we note that the finite adaptation variants of the DM sampler and Scout MCMC both perform similarly to their fully adapting counterparts though they tend to accept fewer proposals.

In addition to sample diagnostics, we also examine the samples themselves in Figure 7. A notable observation is the dramatic imbalance of the parallel tempering samples in Figure 7c, d as well as the DM samples in Figure 7f. RWM also experiences slight imbalance but more notably does not reach far into the tails within the number of iterations. Most surprising is that MALA has not managed to reach the bottom banana regardless of the choice of parameters. We attribute this largely to the gradient pulling proposals away from the tails and thus hampering exploration. The contrast between the DM sampler and Scout MCMC samples highlights the regulating abilities of the Scout chain to help the DM sampler escape from extreme regions. In this example, the left tail in Figure 7f could be considered an extreme region. Both the aRWM and Scout MCMC plots exhibit desirable sampling behavior as the samples are well dispersed over the target and seemingly balanced. However, the primary difference between the two algorithms in this example is the relative concentration of samples due to aRWM’s tendency to reject proposals and stay at the same points whereas Scout MCMC produces a larger number of unique points.

Figure 7. Double Banana Samples. PT and the DM sampler have trouble moving away from tail regions. The best performing algorithms are aRWM and Scout MCMC as they achieve the most accurate sample mean and highest ESJD values. (Note: Red dot indicates starting points and blue dots indicate samples).

Display full size

Finally, we plot the samples generated by the finite versions of the DM sampler and Scout MCMC in Figure 8. The samples presented largely match those of the fully adaptive versions. This indicates that the bank of covariance matrices generated in the adaptive phase is sufficient to produce region specific samples as intended.

Figure 8. Double Banana Samples from the finite variants of the DM Sampler and Scout MCMC. Both algorithms variants perform in line with their perpetually adapting counterparts.

Display full size

4.2. Basis vector distribution

As described briefly in the prelude to this section, the basis vector target consists of a series of normal distributions along the basis vectors in ${\mathbb{R}}^4.$ Formally, the distribution is a mixture of the following normal distributions: $$\begin{array}{c}N\left(10{\mathbf{e}}_1,{\mathbf{I}}_4\right),N\left(-10{\mathbf{e}}_1,{\mathbf{I}}_4\right),N\left(10{\mathbf{e}}_2,{\mathbf{I}}_4\right),N\left(-10{\mathbf{e}}_2,{\mathbf{I}}_4\right),\\ N\left(10{\mathbf{e}}_3,{\mathbf{I}}_4\right),N\left(-10{\mathbf{e}}_3,{\mathbf{I}}_4\right),N\left(10{\mathbf{e}}_4,{\mathbf{I}}_4\right),\text{and}\mathrm{ }N\left(-10{\mathbf{e}}_4,{\mathbf{I}}_4\right)\end{array}$$ where ${\mathbf{e}}_i$ is a basis vector in the ith direction and ${\mathbf{I}}_4$ is the 4D identity matrix.

The key feature of this distribution is that the expected value is at the origin but there is negligible density there. As such, any MCMC algorithm that hopes to be successful must be able to cross a vast low probability desert to move between modes. Each algorithm was run for 40,000 iterations with the first 2,000 iterations discarded as burn-in. The results of the 4D basis vector target are presented in Table 2.

Note that all of RWM, aRWM, MALA, and the DM sampler produced sample means that were far from the origin, the true mean. This indicates that these algorithms did not visit all of the modes in a balanced manner and likely got stuck in one or more select modes. This behavior is not unexpected, however, as these algorithms have no mechanism to cross low probability boundaries. Scout MCMC and parallel tempering, in contrast, produced sample means approaching the origin with Scout MCMC outperforming both cases of parallel tempering. In this instance, there does not appear to be any major impact from increasing the number of chains from two to five in parallel tempering aside from a larger ESJD. Finally, we note that the finite variants of the DM sampler and Scout MCMC perform in line with their respective fully adapting versions, validating their specification as non-adapting approximations.

To confirm that all eight modes were visited (as opposed to, say, two opposing modes with a mean that is equal to the origin), we examine the trace plots of parallel tempering with 5 chains and Scout MCMC in Figure 9. Each trace plot represents one of the dimensions and for this target distribution, we should see the trace plots reaching −10, 0 and 10 in all dimensions. From the plots, it is clear that both algorithms are capable of moving between modes in a frequent manner and that all modes have been visited. It is at this juncture that we turn to the point of efficiency. Notice the difference in acceptance rate and ESJD between Scout MCMC and parallel tempering in Table 2. Scout MCMC has a tendency to accept almost twice as many proposals as parallel tempering even though parallel tempering takes larger steps. This is in part due to parallel tempering having no limit on the frequency of swap moves whereas Scout MCMC is set to only be able to consider swapping every 20 iterations.

Figure 9. 4D basis vector trace plots. Notice that both PT and Scout MCMC are able to reach −10, 0, and 10 in all four dimensions.

Display full size

4.3. Banana bunch distribution

The next target consists of a mixture of 12 banana distributions in ${\mathbb{R}}^3$ arranged such that there is even less interaction than there is in the double banana example. We call this distribution the banana bunch. As this example is in ${\mathbb{R}}^3,$ we cannot simply present the contours. However, the distribution can be understood as the mixture of three groups. The projection of the target on each pair of axes (x-y, x-z, and y-z) appears as the contours in Figure 10.

Figure 10. Projections of the banana bunch distribution on each pair of axes (x-y, x-z, and y-z). (Note: lighter contours indicate higher density, red dot indicates the origin).

Display full size

Combining all three groups will result in the intersection of the apexes of two component distributions at ±40 along each axis. This results in a target with six modes, each far from the origin, and 24 tails extending from the modes toward each other. Once again, we capitalize on the symmetry of our targets and find the expected value to be at the origin. Though the origin has negligible density, scatter plots of samples projected down to planes composed of the basis vectors can be slightly misleading. Figure 11 presents a sample generated directly from the target distribution. The points seemingly at the origin are actually” above” and” below” the origin with respect to the axis missing from the respective plot.

Figure 11. True banana bunch samples. Presented are samples generated directly from the target distribution projected onto the x-y, x-z, and y-z plane.

Display full size

Given the more complex nature of this distribution, we increase the number of samples for all algorithms to 100,000 with the first 1,000 discarded as burn-in. In addition to acceptance rate, first moment, and ESJD, we also consider the second moment, the expectation of the element-wise square of the samples. Such an expectation will assess how well the tails have been explored. A sample that concentrates too heavily in the 6 modes will overshoot this expectation even if it successfully produces a mean near the origin. The true value of this expectation is $\left[\begin{array}{ccc}400& 400& 400\end{array}\right].$ Table 3 presents the results of our experiment.

Perhaps the most surprising result in Table 3 is that the sample mean produced by each algorithm is quite close to the true expected value although aRWM and Scout MCMC are clearly the best performers on that front. Regarding the squared expected value, Scout MCMC performs the best even though the results for parallel tempering with 5 chains presents a convincing case for its convergence properties. We note however that parallel tempering with 2 chains fails to match the performance of either Scout MCMC or parallel tempering with 5 chains. This suggests that 2 chains is not generally sufficient without a more nuanced strategy on the main chain such as using the DM sampler in Scout MCMC. Finally, with respect to ESJD, we note that aRWM and Scout MCMC are clearly the best performers using this efficiency metric. However, we warn the reader to recognize that aRWM accepted less than 2% of proposals, did not produce a sample second moment that was close to the true second moment, and finished the algorithm with a covariance matrix with diagonal $\left[\begin{array}{ccc}862& 1000& 209\end{array}\right].$ This is an indication that as the dimension increases, aRWM relies heavily on expanding its reach to cover the whole region of interest rather than conforming to the shape of the target distribution. The cost of this behavior is that the proposals are not always of high quality and one must hope that it produces enough proposals to extract a decent set of good samples in a limited amount of time. In this case, the enormous proposal distribution was not sufficient to fully explore the tails leading in from the modes toward the origin, thus producing an inaccurate sample second moment.

In order to be complete, however, we must also examine the plots of samples (projected to 2D) to understand whether our algorithms truly explore all modes and tails. Once again, we refer the reader to Figure 11 to view the expected behavior. Given that RWM, MALA, and the DM sampler have no mechanisms for inter-mode movement and that the numerical results reflect this fact, we focus on the remaining algorithms from here. We also focus on the 5 chain version of parallel tempering given that it outperformed the 2 chain version.

The visual results of Figure 12 reflect the behaviors noted by the diagnostics in Table 3. Adaptive RWM struggles to explore the tails as well as the z-axis which results in the poor squared expected value. Parallel tempering, in contrast, performs better in the tails but does not explore the distribution evenly within the specified number of iterations which manifests in the poorer expected value. In addition, there are a number of samples well beyond the modes that are the result of a swap from a tempered chain to the main chain. With Scout MCMC and its finite variant, we see the most appropriate distinction between mode and tail concentrations which manifests in the best squared expectation. They also do so with a higher frequency of points than aRWM or parallel tempering. The highly efficient proposals are realized even when the large distance between modes of the basis vector example is combined with the unusual geometry of the banana examples thus illustrating the ability of Scout MCMC to deliver on the promises of a multimodal sampler that excels at rapid adaptation to local geometry.

Figure 12. Banana Bunch Samples (projected onto the x-y, x-z, and y-z planes). aRWM produces much sparser plots that reflect its inability to explore the tails of the bananas. PT, despite having more samples, still results in inaccurate sample statistics because of its imbalance, with many points skewing to one side. The two Scout MCMC variants produce the most optimal results with a large and balanced set of samples as well as strong exploration of the banana tails.

Display full size

4.4. Bayesian linear regression with horseshoe priors

The final example target distribution is motivated by a more practical scenario. MCMC methods are frequently used in Bayesian inference and there is a lot of work being done to develop methods that can sample from complex posterior distributions. To test the DM sampler and Scout MCMC in such a situation, we construct a Bayesian linear regression scenario with half of the coefficients set to zero. In Bayesian settings, sparsity in regression coefficients can be induced using a horseshoe prior (Carvalho, Polson, and Scott Citation2009). The horseshoe prior is known to be difficult to sample from given the use of the Cauchy distribution. Our example uses 20 independent variables and 1000 observations generated from a standard normal distribution. The true values of the first 10 coefficients are $\left[\begin{array}{cccccccccc}10& 5& 2.5& 2& 1& 0.75& 0.5& 0.25& 0.2& 0.1\end{array}\right]$ and the second set of 10 coefficients are all set to 0. In addition to the 20 coefficients, there is one error variance parameter, 20 horseshoe coefficient parameters, and one global horseshoe parameter for a total of 42 parameters to sample.

We execute each of the algorithms for 50,000 iterations with the first 5000 discarded as burn-in. We modify some of the algorithms to accommodate the higher dimension of this example. For parallel tempering, we execute a 2 chain version as before but replace the 5 chain version with 10 chains. For the DM sampler, we adapt the step size to shrink upon rejecting a sample. Finally, for both the DM sampler and Scout MCMC, we only execute the adaptive versions. To compare across algorithms, we consider the acceptance rate, the distance between the true and estimated regression coefficients, and the ESJD. Table 4 presents the results of the experiment.

As illustrated in the results, the majority of the algorithms tested produced posterior means that are close to the true coefficient values with aRWM being a notable exception. Also of interest are the ESJD values. Of the algorithms that did converge to the true parameter values, only the DM sampler and Scout MCMC had reasonably high ESJD values, suggesting that they sampled most effectively from the posterior. This example illustrates that both algorithms introduced here can be competitive with established MCMC methods in a practical setting.

We conclude this section by repeating the notion of efficient and effective sampling. We find that aRWM may be efficient from an ESJD perspective and parallel tempering is effective as a way to explore different modes, but neither prove to be adequate across the board. Instead, Scout MCMC proves itself as an efficient and effective “smart” sampler by adopting a strategy of combining rapid regional adaptation with heavy tempering for mode swapping.