Markov chain Monte Carlo algorithms with sequential proposals

TL;DR

This paper introduces a flexible sequential-proposal framework for MCMC algorithms, enhancing existing methods like HMC and bouncy particle sampler, leading to improved efficiency and better exploration of complex target distributions.

Contribution

It develops new sequential-proposal MCMC methods, including automatic trajectory tuning in HMC variants and a bouncy particle sampler that improves mixing in multimodal distributions.

Findings

New methods outperform NUTS in efficiency

Sequential proposals enable better exploration of multimodal targets

Automatic trajectory tuning reduces backtracking in HMC

Abstract

We explore a general framework in Markov chain Monte Carlo (MCMC) sampling where sequential proposals are tried as a candidate for the next state of the Markov chain. This sequential-proposal framework can be applied to various existing MCMC methods, including Metropolis-Hastings algorithms using random proposals and methods that use deterministic proposals such as Hamiltonian Monte Carlo (HMC) or the bouncy particle sampler. Sequential-proposal MCMC methods construct the same Markov chains as those constructed by the delayed rejection method under certain circumstances. In the context of HMC, the sequential-proposal approach has been proposed as extra chance generalized hybrid Monte Carlo (XCGHMC). We develop two novel methods in which the trajectories leading to proposals in HMC are automatically tuned to avoid doubling back, as in the No-U-Turn sampler (NUTS). The numerical efficiency of these new methods compare favorably to the NUTS. We additionally show that the sequential-proposal bouncy particle sampler enables the constructed Markov chain to pass through regions of low target density and thus facilitates better mixing of the chain when the target density is multimodal.

Equations

$$\bar{\pi }(x):=\frac{\pi (x)}{Z}$$

$$\int 1_{A\times B}(x,y)P_j(x,dy)\pi (x)dx=\int 1_{B\times A}(x,y)P_j(x,dy)\pi (x)dx,\forall A,B\in \mathcal{X},j=1,2.$$

$$P_1(x,A\setminus \{x\})\ge P_2(x,A\setminus \{x\}),\forall x\in \mathbb{X},\forall A\in \mathcal{X}.$$

$$v(f,P_j):=\underset{M\to \infty }{\lim }M{\text{Var}}_{P_j}\left(\hat{I}(f)\right),j=1,2,$$

$$\min \left(1,\frac{\pi (y)}{\pi (x)}\right).$$

$$\mathcal{P}[X^{(i)}\in A,X^{(i+1)}\in B]=\mathcal{P}[X^{(i)}\in B,X^{(i+1)}\in A]$$

$$\begin{array}{rl}& \mathcal{P}[X^{(i)}\in A,X^{(i+1)}\in B,\text{the }n\text{-th proposal is taken as }{X}^{(i+1)}]\\ & =\int 1_A(y_0)1_B(y_n)\bar{\pi }(y_0)q(y_1|y_0)\cdots q(y_n|y_{n-1})\\ & \cdot 1\left[\Lambda \ge \frac{\pi (y_1)}{\pi (y_0)}\right]\cdots 1\left[\Lambda \ge \frac{\pi (y_{n-1})}{\pi (y_0)}\right]1\left[\Lambda <\frac{\pi (y_n)}{\pi (y_0)}\right]1[0<\Lambda <1]d\Lambda d{y}_{0}d{y}_1\dots d{y}_n,\end{array}$$

$$1\left[\Lambda \ge \frac{\pi (y_1)}{\pi (y_0)}\right]\cdots 1\left[\Lambda \ge \frac{\pi (y_{n-1})}{\pi (y_0)}\right]1\left[\Lambda <\frac{\pi (y_n)}{\pi (y_0)}\right]1[0<\Lambda <1]$$

$$\Lambda \ge \underset{k\in 1:n-1}{\max }\frac{\pi (y_k)}{\pi (y_0)}\text{and}\Lambda <\min \left(1,\frac{\pi (y_n)}{\pi (y_0)}\right).$$

$$\min \left(1,\frac{\pi (y_n)}{\pi (y_0)}\right)-\min \left(1,\frac{\pi (y_n)}{\pi (y_0)},\underset{k\in 1:n-1}{\max }\frac{\pi (y_k)}{\pi (y_0)}\right),$$

$$\begin{array}{rl}& \int 1_A(y_0)1_B(y_n)\bar{\pi }(y_0)\prod _{k=1}^{n}q(y_k|y_{k-1})\cdot \left[\min \left(1,\frac{\pi (y_n)}{\pi (y_0)}\right)-\min \left(1,\frac{\pi (y_n)}{\pi (y_0)},\underset{k\in 1:n-1}{\max }\frac{\pi (y_k)}{\pi (y_0)}\right)\right]d{y}_{0:n}\\ & =\frac{1}{Z}\int 1_A(y_0)1_B(y_n)\prod _{k=1}^{n}q(y_k|y_{k-1})\cdot \left[\min \{\pi (y_0),\pi (y_n)\}-\min \{\pi (y_0),\pi (y_n),\underset{k\in 1:n-1}{\max }\pi (y_k)\}\right]d{y}_{0:n}.\end{array}$$

$$\frac{1}{Z}\int 1_A(y_n)1_B(y_0)\prod _{k=1}^{n}q(y_k|y_{k-1})\left[\min \{\pi (y_n),\pi (y_0)\}-\min \{\pi (y_n),\pi (y_0),\underset{k\in 1:n-1}{\max }\pi (y_k)\}\right]d{y}_{0:n},$$

$$\begin{gathered} \mathcal{P}[X^{(i)}\in A,X^{(i+1)}\in B,\text{the }n\text{-th proposal is taken as }{X}^{(i+1)}] \\ =\mathcal{P}[X^{(i)}\in B,X^{(i+1)}\in A,\text{the }n\text{-th proposal is taken as }{X}^{(i+1)}]. \end{gathered}$$

$$\begin{gathered} \mathcal{P}[X^{(i)}\in A,X^{(i+1)}\in B,\text{all }N\text{ proposals are rejected}] \\ =\mathcal{P}[X^{(i)}\in A,X^{(i)}\in B,\text{all }N\text{ proposals are rejected}], \end{gathered}$$

$$\mathcal{P}[X^{(i)}\in A,X^{(i+1)}\in B]=\mathcal{P}[X^{(i)}\in B,X^{(i+1)}\in A].$$

$$\Lambda <\frac{\pi (Y_n){\prod }_{j=1}^{n}q(Y_{j-1}|Y_j)}{\pi (Y_0){\prod }_{j=1}^n{q}_j(Y_j|Y_{j-1})}.$$

$$\Lambda <\frac{\pi (Y_k){\prod }_{j=1}^k{q}_j(Y_{k-j}|Y_{k-j+1:k})}{\pi (Y_0){\prod }_{j=1}^k{q}_j(Y_j|Y_{j-1:0})}.$$

$$\Lambda <\frac{\pi (Y_k){\prod }_{j=1}^{n-k}{q}_j(Y_{k+j}|Y_{k+j-1:k}){\prod }_{j=n-k+1}^n{q}_j(Y_{n-j}|Y_{n-j+1:n})}{\pi (Y_0){\prod }_{j=1}^n{q}_j(Y_j|Y_{j-1:0})}.$$

$${\alpha }_1(y_0,y_1)=1\wedge \frac{\pi (y_1)q(y_0|y_1)}{\pi (y_0)q(y_1|y_0)},$$

$${\alpha }_2(y_0,y_1,y_2)=1\wedge \frac{\pi (y_2)q(y_1|y_2)q(y_0|y_1)\{1-{\alpha }_1(y_2,y_1)\}}{\pi (y_0)q(y_1|y_0)q(y_2|y_1)\{1-{\alpha }_1(y_0,y_1)\}}.$$

$${\alpha }_n(y_{0:n})=1\wedge \frac{\pi (y_n){\prod }_{j=1}^{n}q(y_{j-1}|y_j){\prod }_{j=1}^{n-1}\{1-{\alpha }_j(y_{n:n-j})\}}{\pi (y_0){\prod }_{j=1}^{n}q(y_j|y_{j-1}){\prod }_{j=1}^{n-1}\{1-{\alpha }_j(y_{0:j})\}}.$$

$${\mathcal{R}}_x\circ {\mathcal{R}}_x=\mathcal{I},$$

$$\frac{\psi ({\mathcal{R}}_{x}v;x)}{\psi (v;x)}\left|\frac{\partial {\mathcal{R}}_{x}v}{\partial v}\right|=1$$

$$\mathcal{T}\circ {\mathcal{S}}_{\tau }\circ \mathcal{T}\circ {\mathcal{S}}_{\tau }=\mathcal{I}.$$

$${\int }_A\Pi (x,v)dv={\int }_{{\mathcal{R}}_x(A)}\Pi (x,v)dv$$

$$\min \left(1,\frac{\Pi ({\mathcal{S}}_{\tau }(x,v))}{\Pi (x,v)}\left|\det D{\mathcal{S}}_{\tau }(x,v)\right|\right),$$

$$\Lambda <\frac{\Pi (Y_n,W_n)}{\Pi (Y_0,W_0)}\left|\det D{\mathcal{S}}_{\tau }^n(Y_0,W_0)\right|,$$

$$H(x,v):=-\log \Pi (x,v)=-\log \bar{\pi }(x)-\log \psi (v;x).$$

$$\psi (v;x)\equiv {\psi }_C(v):=\frac{1}{\sqrt{(2\pi {)}^d\left|\det C\right|}}\exp \left\{-\frac{v^T{C}^{-1}v}{2}\right\}.$$

$$\begin{array}{rl}\frac{dx}{dt}& =C\frac{\partial H}{\partial v}\\ \frac{dv}{dt}& =-C\frac{\partial H}{\partial x}.\end{array}$$

Full text

1 Introduction

Markov chain Monte Carlo (MCMC) methods are widely used to sample from distributions with analytically tractable unnormalized densities. In this paper, we explore a MCMC framework in which proposals for the next state of the Markov chain are drawn sequentially. We consider the objective of obtaining samples from a target distribution on a measurable space $(\mathbb{X},\mathcal{X})$ with density

[TABLE]

with respect to a reference measure denoted by $dx$, where $\pi(x)$ denotes an unnormalized density, and $Z$ denotes the corresponding normalizing constant. MCMC methods construct Markov chains such that, given the current state of the Markov chain $X^{(i)}$, the next state $X^{(i+1)}$ is drawn from a kernel which has the target distribution $\bar{\pi}$ as its invariant distribution. The widely used Metropolis-Hastings (MH) strategy constructs a kernel with a specified invariant distribution in the following two steps (Metropolis et al., 1953; Hastings, 1970). First, a proposal $Y$ is drawn from a proposal kernel, and second, the proposal is accepted as $X^{(i+1)}$ with a certain probability. When the proposal is not accepted, the next state of the chain is set equal to the current state $X^{(i)}$. The acceptance probability depends on the target density and the proposal kernel density at $X^{(i)}$ and $Y$ in a way that ensures that $\bar{\pi}$ is a stationary density of the constructed Markov chain.

The typical size of proposal increments and the mean acceptance probability affect the rate of mixing of the constructed Markov chain and thus the numerical efficiency of the algorithm. There is often a balance to be made between the size of proposal increments and the mean acceptance probability. Theoretical studies on this trade-off have been carried out for several widely used algorithms, such as random walk Metropolis (Roberts et al., 1997), Metropolis adjusted Langevin algorithm (MALA) (Roberts and Rosenthal, 1998), or Hamiltonian Monte Carlo (HMC) (Beskos et al., 2013), in an asymptotic scenario where the target density is given by the product of $d$ identical copies of a one dimensional density and where $d$ tends to infinity. These results suggest that the optimal balance can be made by aiming at a certain value of the mean acceptance probability which depends on the algorithm but not on the target density, provided that the marginal density satisfies some regularity conditions.

Alternative methods to the basic Metropolis-Hastings strategy have been proposed to improve the numerical efficiency beyond the optimal balance between the proposal increment size and the acceptance probability. The multiple-try Metropolis method by Liu et al. (2000) makes multiple proposals given the current state of the Markov chain and select one of them as a candidate for the next state of the Markov chain. Calderhead (2014) proposed a different algorithm that makes multiple proposals and allows more than one of them to be taken as the samples in the Markov chain. Since multiple proposals can be made independently in these methods, parallelization can increase computational efficiency. These methods make a preset number of proposals conditional on the current state in the Markov chain at each iteration.

Developments in various other directions have been made to improve the numerical efficiency of MCMC sampling. Adaptive MCMC methods use transition kernels that adapt over time using the information about the target distribution provided by the past history of the constructed chain (Haario et al., 2001; Andrieu and Thoms, 2008). The update scheme for the transition kernel is designed to induce a sequence of transition kernels that converges to one that is efficient for the target distribution. The convergence of the law of the constructed chain and the rate of convergence have been studied under certain sets of conditions (Haario et al., 2001; Atchadé and Rosenthal, 2005; Andrieu and Moulines, 2006; Andrieu and Atchade, 2007; Roberts and Rosenthal, 2007; Atchadé and Fort, 2010; Atchade and Fort, 2012). Note however that the performance of an adaptive MCMC algorithm is limited by the efficiencies of the candidate transition kernels. In a different approach, Goodman and Weare (2010) proposed using ensemble samplers that construct Markov chains that are equally efficient for all target distributions that are affine transformations of each other. These methods draw information about the shape of the target distribution from parallel chains which jointly target the product distribution given by identical copies of the target density.

There also exist a class of methods that address difficulties in sampling from multimodal distributions using local proposals. Methods in this class include parallel tempering (Geyer, 1991; Hukushima and Nemoto, 1996), simulated tempering (Marinari and Parisi, 1992), and the equi-energy sampler (Kou et al., 2006). In these methods, the mixing of the constructed Markov chain is aided by a set of other Markov chains that target alternative distributions for which the moves between separated modes happen more frequently. The equi-energy sampler bears a similarity with the approach of slice sampling, where a new sample is obtained within a randomly chosen level set of the target density (Roberts and Rosenthal, 1999; Mira et al., 2001a; Neal et al., 2003).

In this paper, we explore a novel approach where proposals are drawn sequentially conditional on the previous proposal in each iteration. The proposal draws continue until a desired number of “acceptable” proposals are made, so the total number of proposals is variable. A key element in this approach is that the decision of acceptance or rejection of proposals are coupled via a single uniform$(0,1)$ random variable drawn at the start of each iteration. This feature makes a straightforward generalization of the Metropolis-Hastings acceptance or rejection strategy. The approach is applicable to a wide range of commonly used MCMC algorithms, including ones that use proposal kernels with well defined densities and others that use deterministic proposal maps, such as Hamiltonian Monte Carlo (Duane et al., 1987) or the recently proposed bouncy particle sampler (Peters et al., 2012; Bouchard-Côté et al., 2018). We will demonstrate that the sequential-proposal approach is flexible; it is possible to make various modifications in order to develop methods that possess specific strengths.

The advantage of the sequential-proposal approach can be explained using Peskun-Tierney ordering (Peskun, 1973; Tierney, 1998; Andrieu and Livingstone, 2019). Suppose two transition kernels $P_{1}$ and $P_{2}$ defined on $(\mathbb{X},\mathcal{X})$ are reversible with respect to $\bar{\pi}$:

[TABLE]

The transition kernel $P_{1}$ is said to dominate $P_{2}$ off the diagonal if

[TABLE]

For a $\mathcal{X}$-measurable function $f$ such that $\int f^{2}(x)\bar{\pi}(x)dx<\infty$, consider an estimator of $\int f(x)\bar{\pi}(x)dx$ given by $\hat{I}(f):=\frac{1}{M}\sum_{i=1}^{M}f(X^{(i)})$. Define a scaled asymptotic variance of the estimator with respect to the kernels $P_{j}$, $j\,{=}\,1,2$ by

[TABLE]

where $X^{(0)}$ is assumed to be drawn from the stationary density $\bar{\pi}$ and $X^{(i)}$ is drawn from $P_{j}(X^{(i-1)},\cdot)$, $i\,{\geq}\,1$. Provided that $P_{1}$ dominates $P_{2}$ off the diagonal, we have $v(f,P_{1})\leq v(f,P_{2})$ (Tierney, 1998, Theorem 4). Suppose now we denote by $P_{1}(x,A\,{\setminus}\,\{x\})$ the conditional probability that $X^{(i)}$ is in $A\,{\setminus}\,\{x\}$ given $X^{(i-1)}\,{=}\,x$ when a sequential-proposal method is used and by $P_{2}(x,A\,{\setminus}\,\{x\})$ the same conditional probability when a standard MH method is used. Then $P_{1}$ dominates $P_{2}$ off the diagonal, because the sequential-proposal method tries additional proposals when the first proposal is rejected. Due to Peskun-Tierney ordering, the asymptotic variance of the estimate $\hat{I}(f)$ from the sequential-proposal method is always smaller than or equal to that from the standard MH. In fact, a similar argument also motivated the development of the delayed rejection method (Tierney and Mira, 1999; Green and Mira, 2001), which we will show to be equivalent to the sequential-proposal approach under certain circumstances in terms of the law of the constructed Markov chains.

After the initial arXiv posting of this manuscript, we became aware of the extra chance generalized hybrid Monte Carlo (XCGHMC) of Campos and Sanz-Serna (2015) that develops a method equivalent to the sequential-proposal approach in the context of HMC. Our work departs from Campos and Sanz-Serna (2015) in two important ways. First, the current paper develops sequential-proposal MCMC more broadly, and as a general recipe for improving on MCMC algorithms. Second, in applying the idea to HMC our emphasis differs from that in Campos and Sanz-Serna (2015). One of the main challenges in using HMC in practice is the tuning of the number of leapfrog jumps. Motivated by this issue, we use the sequential-proposal MCMC framework to develop two novel HMC methods that automatically tune the lengths of the leapfrog trajectories in such a way that the well-known double-backing issue in HMC is avoided, in the same spirit as the No-U-Turn sampler (NUTS) of Hoffman and Gelman (2014). Our numerical results on a multivariate normal distribution and a real data example show that the efficiencies of these new methods compare favorably to that of the NUTS.

Our paper is organized as follows. In Section 2 we explain the sequential-proposal approach in the case where proposal kernels have well defined densities. We call the resulting methods as sequential-proposal Metropolis-Hastings algorithms. An equivalence between between our approach and the delayed rejection method of Mira et al. (2001b) under certain settings is discussed in this section. In Section 3, we explain how the sequential-proposal approach can be applied to a class of MCMC algorithms that use deterministic proposal maps. Based on this formulation, we develop two variants of the NUTS algorithm in Section 4. Section 5 gives a summarizing conclusion. Proofs of most theoretical results are given in Appendices A–E. In Appendix F we propose a novel discrete time bouncy particle sampler method based on the sequential-proposal approach and demonstrate some desirable numerical properties that enable faster mixing of the Markov chain for multimodal target distributions.

2 Sequential-proposal Metropolis-Hastings algorithms

2.1 Sequential-proposal Metropolis algorithm

We will first explain the sequential-proposal approach when the proposal kernel has well defined density with respect to the reference measure of the target density $\bar{\pi}$. For a simpler presentation, we will first describe a sequential-proposal Metropolis algorithm, which uses a proposal kernel with symmetric density. Various generalizations will be introduced in Section 2.2. In standard Metropolis algorithms, given the current state $X^{(i)}\,{=}\,x$ at the $i$-th iteration of the algorithm, the proposal $Y$ is drawn from a probability kernel with conditional density $q(y{\,|\,}x)$ that is symmetric in the sense that $q(y{\,|\,}x)=q(x{\,|\,}y)$ for all $x,y\,{\in}\,\mathbb{X}$. The proposal $Y\,{=}\,y$ is accepted with the probability

[TABLE]

This is often implemented by drawing a uniform random variable $\Lambda\sim\text{unif}(0,1)$ and accepting the proposal by setting $X^{(i+1)}\leftarrow Y$ if and only if $\Lambda<\frac{\pi(y)}{\pi(x)}.$ If $Y$ is not accepted, the algorithm sets $X^{(i+1)}\leftarrow X^{(i)}$.

We will call $Y_{1}$ the first proposal drawn from $q(\cdot{\,|\,}X^{(i)})$. The proposal $Y_{1}$ is rejected if and only if a uniform random number $\Lambda\sim\text{unif}(0,1)$ is greater than or equal to $\pi(Y_{1})/\pi(X^{(i)})$. If rejected, a second proposal $Y_{2}$ is drawn from $q(\cdot{\,|\,}Y_{1})$. The second proposal is accepted if and only if $\Lambda<\pi(Y_{2})/\pi(X^{(i)})$ using the same value of $\Lambda$ used previously. If accepted, the algorithm sets $X^{(i+1)}\leftarrow Y_{2}$. In the case where $Y_{2}$ is rejected, a third proposal is drawn from $q(\cdot{\,|\,}Y_{2})$ and checked for acceptability using the same type of criterion, $\Lambda<\pi(Y_{3})/\pi(X^{(i)})$. This procedure is repeated until an acceptable proposal is found or until a preset number $N$ of proposals are all rejected, whichever is reached sooner. In the case where all $N$ proposals are rejected, the algorithm sets $X^{(i+1)}\leftarrow X^{(i)}$. A pseudocode for a sequential-proposal Metropolis algorithm is given in Algorithm 1. The algorithm reduces to a standard Metropolis algorithm if we set $N\,{=}\,1$.

We will now show that the sequential-proposal Metropolis algorithm just described constructs a reversible Markov chain with respect to the target distribution with density $\bar{\pi}$. Throughout this paper, for two integers $n$ and $m$, we will denote by $n{\hskip 1.70717pt:\hskip 1.70717pt}m$ the sequence $(n,n\,{+}\,1,\dots,m)$ if $n\leq m$ and the sequence $(n,n\,{-}\,1,\dots,m)$ if $n>m$. Also, given a sequence $(a_{n})_{n\geq 0}=(a_{0},a_{1},a_{2},\dots)$, we will denote by $a_{n:m}$ the subsequence $(a_{j})_{n\leq j\leq m}$.

Proposition 1**.**

Algorithm 1 constructs a reversible Markov chain $\left(X^{(i)}\right)$ with respect to the target density $\bar{\pi}$.

Proof.

We will show that the detailed balance equation

[TABLE]

holds for every pair of measurable subsets $A$ and $B$ of $\mathbb{X}$, provided that $X^{(i)}$ is distributed according to $\bar{\pi}$. We will write $Y_{0}:=X^{(i)}$, and the subsequent proposals as $Y_{1},Y_{2},\dots,Y_{N}$. The case where the $n$-th proposal $Y_{n}$ is taken for $X^{(i+1)}$ will be considered; then the claim of detailed balance will follow by combining the cases for $n$ in $1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and the case where all proposals are rejected. Under the assumption that $X^{(i)}$ is distributed according to $\bar{\pi}$, the probability that $X^{(i)}$ is in $A$ and the $n$-th proposal is in $B$ and taken as $X^{(i+1)}$ is given by

[TABLE]

where $\mathbb{1}_{A}$ denotes the indicator function for the set $A$ and $\mathbb{1}[\cdot]$ denotes the indicator function of the event specified between the brackets. The quantity

[TABLE]

is equal to unity if and only if

[TABLE]

It can be readily observed that for real numbers $x$, $a$, and $b$, the conditions $x\geq a$ and $x<b$ are satisfied if and only if $x\in[\min\{a,b\},\,b)$, where the interval length is given by $b\,{-}\,\min(a,b)$. Thus the interval length corresponding to the conditions (3) is given by

[TABLE]

which gives the integral of (2) over $\Lambda$. It follows that (1) is equal to

[TABLE]

If we change the notation of the dummy variables by writing $y_{0}\leftarrow y_{n}$, $y_{1}\leftarrow y_{n-1}$, $\dots$, $y_{n}\leftarrow y_{0}$, then (4) is given by

[TABLE]

where we have used the fact that the kernel density $q$ is symmetric; that is, $q(y_{k-1}{\,|\,}y_{k})=q(y_{k}{\,|\,}y_{k-1})$, for $k\in 1{\hskip 1.70717pt:\hskip 1.70717pt}n$. It is now obvious that (5) is equal to the quantity obtained by swapping the positions of $A$ and $B$ in (1). Thus we see that

[TABLE]

In the case where all $N$ proposals are rejected, the algorithms sets $X^{(i+1)}\leftarrow X^{(i)}$. Thus,

[TABLE]

which is obviously unchanged under the swap of $A$ and $B$. Thus summing (6) over all $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and adding (7) gives

[TABLE]

∎

2.2 Algorithm generalizations

The sequential-proposal Metropolis algorithm described in the previous subsection can be generalized in various ways. Firstly, the algorithm may use proposal kernels with asymmetric density. The $n$-th proposal $Y_{n}$ is drawn from a probability kernel with density $q$ which may not satisfy $q(y{\,|\,}x)=q(x{\,|\,}y)$, $\forall x,y\,{\in}\,\mathbb{X}$. A proposed value $Y_{n}$ is deemed acceptable if and only if

[TABLE]

Here $Y_{0}$ denotes the current state of the Markov chain. Clearly, (8) reduces to, if the proposal density $q$ is symmetric, the acceptance probability $\pi(Y_{n})/\pi(Y_{0})$ in Algorithm 1. We call a sequential-proposal MCMC algorithm that uses a proposal kernel that has possibly asymmetric density a sequential-proposal Metropolis-Hastings algorithm.

A sequential-proposal Metropolis-Hastings algorithm can be further generalized by taking the $L$-th acceptable proposal instead as the next state of the Markov chain for general $L\geq 1$. The algorithms previously described correspond to the case where $L\,{=}\,1$. A pseudocode for this generalized Metropolis-Hastings algorithm is given in Algorithm 2. The algorithmic parameters $N$ and $L$ may be randomly selected at each iteration, provided that they are independent of the proposals $\{Y_{n}{\,;\,}n\,{\geq}\,1\}$ and $\Lambda$. If there are less than $L$ acceptable proposals in the first $N$ proposals, the Markov chain stays at its current position. The proof that Algorithm 2 constructs a reversible Markov chain with respect to the target density $\bar{\pi}$ is given in Appendix A.

A sequential-proposal Metropolis-Hastings algorithm can also employ proposal kernels that depend on the sequence of previous proposals. Suppose that proposals are sequentially drawn in such a way that the $k$-th candidate $Y_{k}$ is drawn from a proposal kernel with density $q_{k}(\cdot{\,|\,}Y_{k-1},\dots,Y_{0})$, where $Y_{k-1}$, $\dots$, $Y_{1}$ denote the previous proposals and $Y_{0}$ denotes the current state $X^{(i)}$ in the Markov chain at the $i$-th iteration. The candidate $Y_{k}$ is deemed acceptable if

[TABLE]

Proposals are sequentially drawn until $L$ acceptable proposals are found. If there are less than $L$ acceptable proposals among the first $N$ proposals, the next state in the Markov chain is set to the current state, $X^{(i+1)}\,{\leftarrow}\,X^{(i)}$. Suppose now that the $L$-th acceptable state is obtained by the $n$-th proposal $Y_{n}$ for some $n\,{\leq}\,N$. In the case where the proposal kernel depends on the sequence of previous proposals, in order to take $Y_{n}$ as the next state of the Markov chain, an additional condition needs to be checked, namely that there are exactly $L\,{-}\,1$ numbers $k\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$ that satisfy

[TABLE]

If this additional condition is satisfied, $Y_{n}$ is taken as the next state of the Markov chain, that is $X^{(i+1)}\,{\leftarrow}\,Y_{n}$. Otherwise, the next state is set to the current state in the Markov chain, $X^{(i+1)}\,{\leftarrow}\,X^{(i)}$. A pseudocode for sequential-proposal Metropolis-Hastings algorithms that employ kernels dependent on the sequence of previous proposals are given in Appendix B. The role of the additional condition (10) is to establish detailed balance between $X^{(i)}$ and $X^{(i+1)}$ by creating a symmetry between the sequence of proposals $Y_{0}\,{\to}\,Y_{1}\,{\to}\,\cdots\,{\to}\,Y_{n}$ and the reversed sequence $Y_{n}\,{\to}\,Y_{n-1}\,{\to}\,\cdots\,{\to}\,Y_{0}$. To see this, we note that the candidate $Y_{n}$ can be taken as the next state of the Markov chain only when there are exactly $L\,{-}\,1$ acceptable proposals among $Y_{1}$, $\dots$, $Y_{n-1}$. The additional symmetry condition accounts for a mirror case where there are $L-1$ acceptable proposals among $Y_{n-1}$, $\dots$, $Y_{1}$, assuming that these proposals are sequentially drawn in the reverse order starting from $Y_{n}$. A proof of detailed balance for this algorithm is also given in Appendix B. This algorithm reduces to Algorithm 2 in the case where the proposal kernel is dependent only on the most recent proposal.

We note that sequential-proposal Metropolis-Hastings algorithms in the case where $L\,{=}\,1$ construct the same Markov chains as those constructed by delayed rejection methods (Tierney and Mira, 1999; Mira et al., 2001b; Green and Mira, 2001) when the proposal kernel depends only on the most recent proposal. A brief description of the delayed rejection method, following Mira et al. (2001b), is given as follows. Given the current state of the Markov chain $y_{0}$, the first candidate value $y_{1}$ is drawn from $q(\cdot{\,|\,}y_{0})$ and accepted with probability

[TABLE]

where $a\land b:=\min(a,b)$. If $y_{1}$ is rejected, a next candidate value $y_{2}$ is drawn from $q(\cdot{\,|\,}y_{1})$. The acceptance probability for $y_{2}$ is given by

[TABLE]

If $y_{1},\dots,y_{n-1}$ are rejected, $y_{n}$ is drawn from $q(\cdot{\,|\,}y_{n-1})$ and accepted with probability

[TABLE]

If all proposals are rejected up to a certain number $N$, the next state of the Markov chain is set to the current state $y_{0}$. We show in Appendix C the equivalence between the delayed rejection method and the sequential-proposal Metropolis-Hastings algorithm with $L\,{=}\,1$ under the case where each proposal is made depending only on the most recent proposal. The delayed rejection method can also use proposal kernels dependent on the sequence of previous proposals to construct a reversible Markov chain with respect to the target distribution. In this case, however, the law of the constructed Markov chain will be different from that by a sequential-proposal Metropolis-Hastings algorithm.

In our view, there are several advantages sequential-proposal Metropolis-Hastings algorithms have over the delayed rejection method:

Sequential-proposal Metropolis-Hastings algorithms are more straightforward to implement than the delayed rejection method. The evaluation of $\alpha_{n}(y_{0:n})$ in delayed rejection involves the evaluation of a sequence of reversed acceptance probabilities $\{\alpha_{j}(y_{n:n-j}){\,;\,}j\in 1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1\}$. This involves computation of a total of $\mathcal{O}(n^{2})$ acceptance probabilities. In comparison, sequential-proposal Metropolis-Hastings algorithms only compare the ratio in (8) to a uniform random number $\Lambda$ for the same task of checking the acceptability of $y_{n}$. The algorithmic simplicity of sequential-proposal Metropolis-Hastings facilitates the use of a large number of proposals in each iteration. Moreover, one may choose to take the $L$-th acceptable proposal for the next state of the Markov chain for a large $L\,{>}\,1$. 2. 2.

The sequential-proposal MCMC framework can be readily applied to MCMC algorithms using deterministic maps for proposals, as explained in Section 3. In particular, the sequential-proposal MCMC framework applies to Hamiltonian Monte Carlo and the bouncy particle sampler methods, leading to improved the numerical efficiency. Applications to these algorithms are discussed in Section 4 and Appendix F. We note that the delayed rejection method has been generalized to algorithms using deterministic maps in Green and Mira (2001), although only the case for the second proposal was discussed. 3. 3.

The conceptual simplicity of the sequential-proposal MCMC framework allows for various generalizations and modifications. For example, in Section 4.2, we develop sequential-proposal No-U-Turn sampler algorithms (Algorithms 6 and 7) that automatically adjust the lengths of trajectories leading to proposals in HMC, similarly to the No-U-Turn sampler algorithm proposed by Hoffman and Gelman (2014). The proofs of detailed balance for these algorithms can be obtained by making minor modifications to the proof for the sequential-proposal Metropolis-Hastings algorithms.

3 Sequential-proposal MCMC algorithms using deterministic kernels

The sequential-proposal MCMC framework can be applied to algorithms that use deterministic proposal kernels. MCMC algorithms that employ deterministic proposal kernels often target a distribution on an extended space $\mathbb{X}\,{\times}\,\mathbb{V}$ whose the marginal distribution on $\mathbb{X}$ is equal to the original target distribution $\bar{\pi}$. An additional variable $V$ drawn from a distribution on $\mathbb{V}$ serves as a parameter for the deterministic proposal kernel. In this section, we will explain a general class of MCMC algorithms using deterministic proposal kernels and show how the sequential-proposal scheme can be applied to these algorithms. Applications to specific algorithms, such as HMC or the bouncy particle sampler (BPS), are discussed in subsequent sections (Section 4 and Appendix F).

We suppose that the extended target distribution on $\mathbb{X}\times\mathbb{V}$ has density $\Pi(x,v)$ with respect to a reference measure denoted by $dx\,dv$. We further assume that the original target density $\bar{\pi}$ equals the marginal density of $\Pi$, such that $\Pi(x,v)=\bar{\pi}(x)\psi(v{\,;\,}x)$ for some $\psi(v{\,;\,}x)$, the conditional density of $v$ given $x$. We define a collection of deterministic maps $\mathcal{S}_{\tau}:\mathbb{X}\times\mathbb{V}\to\mathbb{X}\times\mathbb{V}$ for possibly various values of $\tau$. In HMC and the BPS, $\mathcal{S}_{\tau}$ has an analogy with the evolution of a particle in a physical system for a time duration $\tau$. In this analogy, the variable $x\in\mathbb{X}$ is considered as the position of a particle in the system and the variable $v\in\mathbb{V}$ as the velocity of the particle. The point $\mathcal{S}_{\tau}(x,v)$ then represents the final position-velocity pair of a particle that moves with initial position $x$ and initial velocity $v$ for time $\tau$. We suppose that the map $\mathcal{S}_{\tau}$ for each $\tau$ satisfies the following condition:

Reversibility condition. There exists a velocity reflection operator $\mathcal{R}_{x}:\mathbb{V}\to\mathbb{V}$ defined for every point $x\in\mathbb{X}$ such that

[TABLE]

holds for every $x\,{\in}\,\mathbb{X}$ and

[TABLE]

holds for almost every $(x,v)\,{\in}\,\mathbb{X}\times\mathbb{V}$ with respect to the reference measure $dx\,dv$. Furthermore, if we define a map $\mathcal{T}:\mathbb{X}\times\mathbb{V}\to\mathbb{X}\times\mathbb{V}$ as $\mathcal{T}(x,v):=(x,\mathcal{R}_{x}v)$, we have

[TABLE]

Similar sets of conditions appear routinely in the literature on MCMC (Fang et al., 2014; Vanetti et al., 2017) and on Hamiltonian dynamics (Leimkuhler and Reich, 2004, Section 4.3). In (11) and (13), $\mathcal{I}$ denotes the identity map in the corresponding space $\mathbb{V}$ or $\mathbb{X}\times\mathbb{V}$, and the symbol $\circ$ denotes function composition. In (12), $\left|\frac{\partial\mathcal{R}_{x}v}{\partial v}\right|$ denotes the absolute value of the Jacobian determinant of the map $\mathcal{R}_{x}$ at $v$. The condition (12) is equivalent to the condition that

[TABLE]

for every measurable subset $A$ of $\mathbb{V}$ and for almost every $x\,{\in}\,\mathbb{X}$, due to the change of variable formula. The condition (13) can be understood as an abstraction of a property in Hamiltonian dynamics that if we reverse the velocity of a particle and advance in time, the particle traces back its past trajectory.

Given $X^{(i)}\,{=}\,x$ and $V^{(i)}\,{=}\,v$ at the start of the $i$-th iteration, a MCMC algorithm can make a deterministic proposal $\mathcal{S}_{\tau}(x,v)$, which is accepted with probability

[TABLE]

where $D$ denotes the differential operator (i.e., $D\mathcal{S}_{\tau}(x,v)=\frac{\partial\mathcal{S}_{\tau}(x,v)}{\partial(x,v)}$). In algorithms such as HMC or the BPS, the extended target density $\Pi(x,v)$ is often taken as a product of independent densities, $\bar{\pi}(x)\psi(v)$, where a common choice for $\psi(v)$ is a multivariate normal density. The map $\mathcal{S}_{\tau}$ is often taken to preserve the reference measure, such that it has unit Jacobian determinant (i.e., $\left|\det D\mathcal{S}_{\tau}(x,v)\right|=1$, for all $(x,v)$).

The sequential-proposal framework can be used to generalize MCMC algorithms using deterministic kernels in a similar way that it is applied to Metropolis-Hastings algorithms. A pseudocode of a sequential-proposal MCMC algorithm using a deterministic kernel is shown in Algorithm 3. Proposals are obtained sequentially as $(Y_{n},W_{n})\leftarrow\mathcal{S}_{\tau}(Y_{n-1},W_{n-1})$, where we write $(Y_{0},W_{0}):=(X^{(i)},V^{(i)})$. The pair $(Y_{n},W_{n})$ is deemed acceptable if

[TABLE]

where $\mathcal{S}_{\tau}^{n}=\mathcal{S}_{\tau}\circ\cdots\circ\mathcal{S}_{\tau}$ denotes a map obtained by composing $\mathcal{S}_{\tau}$ $n$ times. If there are less than $L$ acceptable proposals in the sequence of $L$ proposals, the next state of the Markov chain is set to $(X^{(i+1)},V^{(i+1)})\leftarrow(X^{(i)},\mathcal{R}_{X^{(i)}}V^{(i)}).$ The velocity $V^{(i+1)}$ may be refreshed at the end of the iteration by drawing from $\psi(\,\cdot{\,;\,}X^{(i+1)})$ with a certain probability $p^{\text{ref}}(X^{(i+1)})$ that may depend on $X^{(i+1)}$. The parameter $\tau$ for the evolution map $\mathcal{S}_{\tau}$ can be drawn randomly. The pseudocode in Algorithm 3 shows the case where $\tau$ is drawn once per iteration and the same value is used for all $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, but $\tau$ can also be drawn separately for each $n$, provided that the draws are independent of each other and of all other random draws in the algorithm.

We state the following result for Algorithm 3. The proof is given in Appendix D.

Proposition 2**.**

The extended target distribution with density $\Pi(x,v)$ is a stationary distribution for the Markov chain $\left(X^{(i)},V^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 3. Furthermore, the Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 3, marginally for the $x$-component, is reversible with respect to the target distribution $\bar{\pi}(x)$.

4 Connection to Hamiltonian Monte Carlo methods

4.1 Sequential-proposal Hamiltonian Monte Carlo

In this section, we consider applications of the sequential-proposal approach described in Section 3 to Hamiltonian Monte Carlo algorithms and discuss the numerical efficiency. We first briefly summarize basic features of HMC algorithms. A function on $\mathbb{X}\times\mathbb{V}$, called the Hamiltonian, is defined as the negative log density of the extended target density:

[TABLE]

We assume both $\mathbb{X}$ and $\mathbb{V}$ are equal to the $d$ dimensional Euclidean space $\mathbb{R}^{d}$. The velocity distribution $\psi(v{\,;\,}x)$ is often taken as a multivariate normal density independent of $x$,

[TABLE]

An analogy with a physical Hamiltonian system is drawn by interpreting the first term $-\log\bar{\pi}(x)$ as the static potential energy of a particle and the second term $-\log\psi(v)$ as the kinetic energy. In this analogy, the covariance matrix $C$ can be interpreted as the inverse of the mass of the particle. Hamiltonian dynamics is defined as a solution to the Hamiltonian equation of motion (HEM):

[TABLE]

If we denote the solution to the HEM as $(x(t),v(t))$, the exact Hamiltonian flow $S_{\tau}^{*}$ defined by $\mathcal{S}_{\tau}^{*}(x(0),v(0)):=(x(\tau),v(\tau))$ satisfies the reversibility condition (12) and (13) when the velocity reflection operator is given by $\mathcal{R}_{x}(v)\,{=}\,{-}v$ for all $x\,{\in}\,\mathbb{X}$ and $v\,{\in}\,\mathbb{V}$. The map $\mathcal{S}_{\tau}^{*}$ preserves the Hamiltonian, that is, $H(x,v)=H(\mathcal{S}_{\tau}^{*}(x,v))$ for all $x\,{\in}\,\mathbb{X}$, $v\,{\in}\,\mathbb{V}$, and $\tau\,{\geq}\,0$. The map $\mathcal{S}_{\tau}^{*}$ also preserves the reference measure $dx\,dv$, that is, $\left|\det D\mathcal{S}_{\tau}^{*}(x,v)\right|=1$ for all $x\in\mathbb{X}$, $v\in\mathbb{V}$ and $\tau\geq 0$, which is known as Liouville’s theorem (Liouville, 1838).

A commonly used numerical approximation method for solving the HEM is called the leapfrog method (Duane et al., 1987; Leimkuhler and Reich, 2004). One iteration of the leapfrog method approximates time evolution of a Hamiltonian system for duration $\epsilon$ by alternately updating the velocity and position $(x,v)$ as follows:

[TABLE]

We call the time increment $\epsilon$ the leapfrog step size.

A standard Hamiltonian Monte Carlo algorithm is a specific instance of MCMC algorithms using deterministic kernels described in Section 3, where the extended target density $\Pi(x,v)$ is given by $\bar{\pi}(x)\psi_{C}(v)$ and the proposal map $\mathcal{S}_{\tau}$ is given by $l$ leapfrog jumps with step size $\epsilon$, such that the time duration parameter $\tau$ can be understood as the pair $(\epsilon,l)$. The reversibility condition (11)–(13) is satisfied by this $\mathcal{S}_{\tau}$ with $\mathcal{R}_{x}\,{=}\,{-}\mathcal{I}$ for all $x\,{\in}\,\mathbb{X}$. Each step in the leapfrog method (17) preserves the reference measure $dx\,dv$, so we have $|\det D\mathcal{S}_{\tau}|\equiv 1$. It is common to refresh the probability at every iteration (i.e., $p^{\text{ref}}(x)\equiv 1$).

Sequential-proposal HMC (Algorithm 4) is obtained as a specific case of sequential-proposal MCMC algorithms using deterministic kernels (Algorithm 3) under the same setting, $\Pi(x,v)=\bar{\pi}(x)\psi_{C}(v)$ and $\mathcal{S}_{\tau}=\mathcal{S}_{(\epsilon,l)}$. In other words, a proposal $(Y_{1},W_{1})$ is made by making $l$ leapfrog jumps of size $\epsilon$ starting from $(Y_{0},W_{0})$, and if the proposal is rejected, a new proposal $(Y_{2},W_{2})$ is made by making $l$ leapfrog jumps from $(Y_{1},W_{1})$. The procedure is repeated until $L$ acceptable proposals are found, or until $N$ proposals have been tried, whichever comes sooner. The leapfrog jump size $\epsilon$ and the unit number of jumps $l$ may be re-drawn at every iteration or for every new proposal. As mentioned earlier, Campos and Sanz-Serna (2015) has proposed extra chance generalized hybrid Monte Carlo (XCGHMC), which is identical to the sequential-proposal approach, except possibly in the way the velocity is refreshed at the end of each iteration. In generalized HMC (Horowitz, 1991), the velocity is partially refreshed by setting

[TABLE]

where $V$ is the velocity before refreshement, $U$ is an independent draw from $\mathcal{N}(0,C)$, and $\theta$ is an arbitrary real number. It was shown in Campos and Sanz-Serna (2015) that Markov chains constructed by XCGHMC have the same law as those constructed by Look Ahead Hamiltonian Monte Carlo (LAHMC) developed by Sohl-Dickstein et al. (2014).

A major advantage of HMC algorithms over random walk based algorithms such as random walk Metropolis or Metropolis adjusted Langevin algorithms is that HMC can make a global jump in one iteration (Neal, 2011). The leapfrog method is able to build long trajectories that are numerically stable, provided that the target distribution satisfies some regulatory conditions and the leapfrog step size is less than a certain upper bound (Leimkuhler and Reich, 2004). Since the solution to the HEM preserves the Hamiltonian, proposals obtained by a numerical approximation to the solution can be accepted with reasonably high probabilities. Given a fixed length of leapfrog trajectory, the number of leapfrog jumps is inversely proportional to the leapfrog jump size. Thus an increase in the leapfrog step size leads to a reduced number of evaluations of the gradient of the target density. On the other hand, decreasing the leapfrog step size tends to increase the mean acceptance probability. As $\epsilon\to 0$, the average increment in the Hamiltonian at the end of the leapfrog trajectory scales as $\epsilon^{4}$ (Leimkuhler and Reich, 2004). In an asymptotic scenario where the target distribution is given by a product of $d$ independent, identical low dimensional distributions and $d$ tends to infinity, the increment in the Hamiltonian converges in distribution to a normal distribution with mean $\mu\epsilon^{4}d$ and variance $2\mu\epsilon^{4}d$ for some constant $\mu\,{>}\,0$ dependent on the target density $\bar{\pi}$ (Gupta et al., 1990; Neal, 2011). Beskos et al. (2013) showed under some mild regulatory conditions on the target density that as $\epsilon=\epsilon_{0}d^{-1/4}$ and $d\to\infty$, the mean acceptance probability tends to $a(\epsilon_{0}):=2\Phi\left(-\epsilon_{0}^{2}\sqrt{\mu/2}\right)$ where $\Phi(\cdot)$ denotes the cdf of the standard normal distribution. The computational cost for obtaining an accepted proposal that is fixed distance away from the current state in HMC is approximately given by

[TABLE]

which is minimized when $a(\epsilon_{0})=0.651$ to three decimal places (Beskos et al., 2013; Neal, 2011). Empirical results also support targeting the mean acceptance probability of around 0.65 (Sexton and Weingarten, 1992; Neal, 1994). HMC using sequential proposals can improve on the numerical efficiency by increasing the probability that the constructed Markov chain makes a nonzero move at each iteration. A numerical study in Section 4.4 shows that HMC with sequential proposals leads to higher effective sample sizes per computation time compared to the standard HMC on a toy model.

4.2 Sequential-proposal No-U-Turn sampler algorithms

As previously mentioned, a key advantage of HMC over random walk based methods comes from its ability to make long moves. If the number of leapfrog jumps is too small, the Markov chain from HMC may essentially behave like a random walk because the velocity is randomly refreshed before a long leapfrog trajectory is built. Conversely, if the number of leapfrog jumps is too large, the trajectory may double back on itself, since the solution to the Hamiltonian equation of motion is confined to a level set of the Hamiltonian. However, simply stopping the leapfrog jumps when the trajectory starts doubling back on itself generally destroys the detailed balance of the Markov chain with respect to the target distribution. In order to solve this issue, Hoffman and Gelman (2014) proposed the No-U-Turn sampler (NUTS). In this section, we will briefly explain the NUTS algorithm and discuss the connection with sequential-proposal framework. In addition, we will propose two new algorithms that address the same issue of trajectory doubling.

In the No-U-Turn sampler, leapfrog trajectories are repeatedly extended twice in size in either forward or backward direction in the form of binary trees, until a “U-turn” is observed (see Figure 1). The binary tree starts from the initial node $(X^{(i)},V)$, where $V$ is the velocity drawn from the standard multivariate normal distribution at the beginning of the $i$-th iteration. The direction of binary tree expansion is determined by a sequence of $\text{unif}(\{-1,1\})$ variables denoted by $(\sigma_{j})_{j\geq 0}$. The expansion of the binary tree stops if a U-turn is observed between the two leaf nodes on opposite sides of any of the sub-binary trees of the current tree. A position-velocity pair $(x,v)$ and another pair $(x^{\prime},v^{\prime})$ that is ahead of $(x,v)$ on a leapfrog trajectory are said to satisfy the U-turn condition if either

[TABLE]

where $\cdot$ denotes the inner product in Euclidean spaces. If there is a U-turn within the lastly added half of the current binary tree, the other half without a U-turn is taken as the final binary tree. On the other hand, if a U-turn is only observed between the two opposite leaf nodes of the current binary tree but not within any of the sub-binary trees, the current binary tree is taken as the final binary tree. The next state of the Markov chain $X^{(i+1)}$ is set to one of the acceptable leaf nodes in the final binary tree. A leaf node $(x,v)$ is deemed acceptable if $\frac{\Pi(x,v)}{\Pi(X^{(i)},V)}>\Lambda$. Here $\Pi(x,v):=\bar{\pi}(x)\psi_{I_{d}}(v)$, where $\psi_{I_{d}}$ denotes the density of the $d$ dimensional standard normal distribution. Hoffman and Gelman (2014) gives two versions of the NUTS algorithm. The naive version selects the next state of the Markov chain uniformly at random among the acceptable leaf nodes in the final binary tree. The efficient version preferentially selects a random leaf node in sub-binary trees that are added later. A pseudocode for these NUTS algorithms is given in Algorithm 5. By construction, for every leaf node in the final binary tree, it is possible to build the same final binary tree starting from that leaf node using a unique sequence of directions. Since each direction is drawn from $\text{unif}(\{-1,1\})$, the probability of constructing the final binary tree is the same when started from any of its leaf nodes. This symmetric relationship ensures that the constructed Markov chain is reversible with respect to the target distribution.

The NUTS algorithm shares with the sequential-proposal MCMC framework the key feature that the decisions of acceptance or rejection of proposals are mutually coupled via a single uniform$(0,1)$ random variable drawn at the start of each iteration. Furthermore, the naive version of the algorithm (Algorithm 2 in Hoffman and Gelman (2014)) can be viewed as a specific case of the sequential-proposal MCMC algorithm as follows. At each iteration a binary tree starting from $(X^{(i)},V)$ is expanded until a U-turn is observed, as described above. Proposals are made sequentially by selecting one of the leaf nodes of the final binary tree uniformly at random. The first proposal that is acceptable is taken as the next state of the Markov chain. Since the next state of the Markov chain is then selected uniformly at random among the acceptable leaf nodes in the final binary tree, this sequential-proposal approach is equivalent to the naive NUTS.

There are two features of the NUTS algorithm that may, unfortunately, compromise the numerical efficiency. First, the point chosen for the next state of the Markov chain is generally not the farthest point on the leapfrog trajectory from the initial point. The NUTS typically constructs a leapfrog trajectory that is longer than the distance between the initial point and the point selected for the next state of the Markov chain due to the requirement of detailed balance. Second, the NUTS evaluates the log target density at every point on the constructed leapfrog trajectory to determine the acceptability. This can result in a substantial overhead if the computational cost of evaluating the log target density is at least comparable to that of evaluating the gradient of the log target density. We propose two alternative No-U-Turn sampling algorithms, which we call spNUTS1 and spNUTS2, addressing these two issues.

In spNUTS1, leapfrog trajectories are extended in one direction according to a given length schedule until a U-turn is observed, and only the endpoint of the trajectory is checked for acceptability. If the endpoint is not acceptable, a new trajectory is started from that point with a refreshed velocity. A pseudocode of spNUTS1 is given in Algorithm 6. At the start of each iteration, a velocity vector is drawn from a multivariate normal distribution $\mathcal{N}(0,C)$ where $C$ is a $d\,{\times}\,d$ positive definite matrix. A leapfrog trajectory started from the current state of the Markov chain and the drawn velocity vector, denoted by $(x_{0},v_{0})$, is repeatedly extended in units of $l$ jumps. The position-velocity pair after $lk$ leapfrog jumps is denoted by $(x_{k},v_{k})$. We note that the leapfrog updates (17) should also use the same matrix $C$ as the covariance of the velocity distribution. At preset checkpoints determined by a finite increasing sequence $(b_{j})_{j\in 1:j_{\max}}$, the algorithm calculates the angles between the displacement $x_{b_{j}}{-}x_{0}$ and the velocities $v_{0}$ and $v_{b_{j}}$. In order to take into account the given covariance structure $C$, we define a $C$-norm of a vector $x\in\mathbb{R}^{d}$ as

[TABLE]

and the cosine of the angle between two vectors $x$ and $x^{\prime}$ as

[TABLE]

The leapfrog trajectory stops at $(x_{b_{j}},v_{b_{j}})$ if either of the following inequalities hold for a given $c$:

[TABLE]

The value of $c$ can be fixed at a constant value or randomly drawn for each trajectory. Algorithm 6 describes a case where $c$ is randomly drawn from a distribution denoted by $\zeta$. If the above stopping condition (20) is not satisfied until $j\,{=}\,j_{\max}$, the trajectory stops at $(x_{b_{j_{\max}}},v_{b_{j_{\max}}})$.

The final state at the stopped trajectory $(x_{b_{j}},v_{b_{j}})$ makes the first proposal $(Y_{1},W^{\prime}_{1})$. It is taken as the next state in the Markov chain if the following two conditions are met. First, the state $(x_{b_{j}},v_{b_{j}})$ has to be acceptable by satisfying

[TABLE]

Since the Hamiltonian of the state $(x_{b_{j}},v_{b_{j}})$ is given by $-\bar{\pi}(x_{b_{j}}){-}\log\psi_{C}(v_{b_{j}})$, the acceptability criterion (21) can be interpreted as that the increase in the Hamiltonian compared to the initial state $(x_{0},v_{0})$ is at most $-\log\Lambda$. The second required condition is that

[TABLE]

Since the trajectory has been extended to $(x_{b_{j}},v_{b_{j}})$, the stopping condition (20) was not satisfied between the initial state $(x_{0},v_{0})$ and any of the previously visited states $\{(x_{b_{j^{\prime}}},v_{b_{j^{\prime}}}){\,;\,}1\,{\leq}\,j^{\prime}\,{<}\,j\}$. When the trajectory is viewed in the reverse order, an analogous situation is that the final state $(x_{b_{j}},v_{b_{j}})$ and the intermediate states $\{(x_{b_{j}-b_{j^{\prime}}},v_{b_{j}-b_{j^{\prime}}}){\,;\,}1\,{\leq}\,j^{\prime}\,{<}\,j\}$ satisfy (22). This symmetry condition is necessary to establish detailed balance of the Markov chain. If the symmetry condition is not satisfied, the next state of the Markov chain is set to the current state $x_{0}$. In the case where both the acceptability condition (21) and the symmetry condition (22) are satisfied, $x_{b_{j}}$ is taken as the next state of the Markov chain. If the symmetry condition is satisfied but the acceptability condition is not, the algorithm starts a new leapfrog trajectory from $x_{b_{j}}$ with a new velocity $W_{1}$. The new velocity $W_{1}$ is obtained by drawing a random vector from the velocity distribution $\psi_{C}$ and rescaling it such that the $C$-norm is preserved: $\|W_{1}\|_{C}=\|W^{\prime}_{1}\|_{C}$, or $\log\psi_{C}(W_{1})=\log\psi_{C}(W^{\prime}_{1})$. The new trajectory is extended until the stopping condition (20) with respect to the starting pair $(Y_{1},W_{1})$ is satisfied. This procedure is repeated for subsequent proposals $(Y_{n},W^{\prime}_{n})$, $n\,{\geq}\,2$. If all $N$ proposals $\{(Y_{n},W^{\prime}_{n}){\,;\,}n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N\}$ are unacceptable, the next state of the Markov chain is set to the current state $X^{(i)}$. The Markov chain also stays at its current state if the symmetry condition is not satisfied by at least one of the constructed leapfrog trajectories.

Algorithm 6 can be considered as a special case of sequential-proposal MCMC algorithms using deterministic kernels (Algorithm 3) where $L\,{=}\,1$ and the proposal map $\mathcal{S}_{\tau}$ is given by a function that maps the starting state $(x_{0},v_{0})$ to the final state $(x_{b_{j}},v_{b_{j}})$ of the stopped trajectory, but with a differing feature that the direction of the velocity is randomly refreshed between proposals. In practice, extending the trajectory length at an exponential rate by setting, for example, $b_{j}\,{=}\,2^{j-1}$ can lead to a high probability that the symmetry condition is satisfied. Hoffman and Gelman (2014) also remarked (based on Figure 5 in their paper) that as the size of binary trees were repeatedly doubled, the U-turn condition was satisfied most of the time only by the two opposite leaf nodes of the final binary tree but not by the opposite leaf nodes of any of the sub-binary trees, for all examples they considered including ones arising from practical applications. The choice $b_{j}\,{=}\,2^{j-1}$ also makes it easy to predict the checkpoints for the symmetry condition, $\{b_{j}\,{-}\,b_{j^{\prime}}{\,;\,}j^{\prime}\,{\leq}\,j{-}1\}$. We note that the numerical efficiency of Algorithm 6, as well as that of the original NUTS algorithm, can be improved by tuning the covariance $C$ (see the numerical results in Section 4.4).

The computational efficiency of Algorithm 6 may be compared favorably to that of the NUTS algorithm. Some numerical results are given in Section 4.4. Suppose that the average computational cost of evaluating the log target density is denoted by $c_{\pi}$ and that of evaluating the gradient of the log target density by $c_{\triangledown}$. Assume also that in the NUTS and spNUTS1 (Algorithm 6), the computational cost of checking the U-turn condition such as (20) is denoted by $c_{U}$. If a leapfrog trajectory stops after making $\xi$ sets of $l$ unit jumps on average, the NUTS algorithm evaluates the log target density $\xi$ times and checks the U-turn condition $\xi$ times on average. Thus the average computational cost for one iteration of the NUTS algorithm is given by $(lc_{\triangledown}+c_{\pi}+c_{U})\cdot\xi$. In comparison, spNUTS1 evaluates the log target density once and checks the U-turn condition $2\log_{2}\xi\,{+}\,1$ times if $b_{j}\,{=}\,2^{j-1}$ for $j\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}j_{\max}$. The average computational cost of obtaining a proposal in spNUTS1 is given by $lc_{\triangledown}\xi+c_{\pi}+c_{U}(2\log_{2}\xi\,{+}\,1)$, and the average cost of finding a new state for the Markov chain different from the current state is roughly given by \frac{1}{a\cdot\tilde{a}}\big{(}lc_{\triangledown}\xi\,{+}\,c_{\pi}\,{+}\,c_{U}(2\log_{2}\xi\,{+}\,1)\big{)}, where $a$ denotes the mean acceptance probability of a proposal and $\tilde{a}$ denotes the average probability that the symmetry condition is satisfied. Both $a$ and $\tilde{a}$ can be made close to unity in practice, so there is a computational gain in using spNUTS1 over the original NUTS if $\xi$ is large and $c_{\pi}$ is at least comparable to $c_{\triangledown}$. The number $l$ can be chosen to one unless there is an issue of numerical instability of leapfrog trajectories. We note that the increase in the cost by a factor of $\frac{1}{a}$ can be partially negated, in terms of the overall numerical efficiency, due to the fact if a proposal is deemed unacceptable, the next proposal can be further away from the initial state $Y_{0}$. The average distance between two consecutive states in the constructed Markov chain is a measure widely used to evaluate the numerical efficiency of a MCMC algorithm (Sherlock et al., 2010).

The proof of the following proposition is given in Appendix E.

Proposition 3**.**

The Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by the sequential-proposal No-U-Turn sampler of type 1 (spNUTS1, Algorithm 6) is reversible with respect to the target density $\bar{\pi}$.

Another algorithm that automatically tunes the lengths of leapfrog trajectories, called spNUTS2, is given in Algorithm 7. Unlike spNUTS1, spNUTS2 applies the sequential-proposal scheme within one trajectory. The spNUTS2 algorithm takes the endpoint of the constructed leapfrog trajectory as a candidate for the next state of the Markov chain, as in spNUTS1. However, it evaluates the log target density at every point on the trajectory like the original NUTS. Starting from the current state of the Markov chain $X^{(i)}\,{=}\,x_{0}$ and a velocity vector $v_{0}$ randomly drawn from $\psi_{C}$, the algorithm extends a leapfrog trajectory in units of $l$ leapfrog jumps. We will denote by $(x_{1},v_{1})$ the first acceptable state along the trajectory that is a multiple of $l$ leapfrog jumps away from the initial state. Here, $(x_{1},v_{1})$ is acceptable if

[TABLE]

In order to avoid indefinitely extending the trajectory when the leapfrog approximation is numerically unstable, the algorithm ends the attempt to find the next acceptable state if $N$ consecutive states at intervals of $l$ leapfrog jumps are all unacceptable. In this case, the next state of the Markov chain is set to $(x_{0},v_{0})$. For $k\,{\geq}\,2$, the state $(x_{k},v_{k})$ is likewise found as the first acceptable state along the leapfrog trajectory that is a multiple of $l$ jumps from $(x_{k-1},v_{k-1})$. If for any $k\,{\geq}\,1$ the next acceptable state is not found in $N$ consecutive states visited after $(x_{k-1},v_{k-1})$, the next state in the Markov chain is also set to $(x_{0},v_{0})$. In practice, however, this situation can be avoided by taking the leapfrog step size $\epsilon$ reasonably small to ensure numerical stability and $N$ large enough. The algorithm takes a preset increasing sequence of integers $(b_{j})_{j\in 1{\hskip 1.70717pt:\hskip 1.70717pt}j_{\max}}$ and checks if the angles between the displacement vector $x_{b_{j}}\,{-}\,x_{0}$ and the initial and the last velocity vectors $v_{0}$ and $v_{b_{j}}$ are below a certain level $c$. The trajectory is stopped at $(x_{b_{j}},v_{b_{j}})$ if either

[TABLE]

Upon reaching $(x_{b_{j_{\max}}},v_{b_{j_{\max}}})$, however, the trajectory stops regardless of whether (23) is satisfied for $j\,{=}\,j_{\max}$. As in spNUTS1, a symmetry condition is checked to ensure detailed balance. That is, the state $(x_{b_{j}},v_{b_{j}})$ is taken as the next state in the Markov chain if and only if

[TABLE]

If the symmetry condition is not satisfied, the next state of the Markov chain is set to $(x_{0},v_{0})$. As in spNUTS1, the choice of $b_{j}\,{=}\,2^{j-1}$, $j\in 1{\hskip 1.70717pt:\hskip 1.70717pt}j_{\max}$, allows the symmetry condition in spNUTS2 to be satisfied with high probability and makes the checkpoints for the symmetry condition, $\{b_{j}\,{-}\,b_{j^{\prime}}{\,;\,}j^{\prime}\,{\leq}\,j{-}1\}$, readily predictable.

When $c_{\pi}$, $c_{\triangledown}$, and $c_{U}$ denote the same average computational costs as before and $b_{j}\,{=}\,2^{j-1}$, the average computational cost of finding a distinct sample point for the Markov chain using spNUTS2 is roughly given by $\frac{1}{\tilde{a}}\{(lc_{\triangledown}+c_{\pi})\cdot\xi+c_{U}(2\log_{2}(a\xi)\,{+}\,1)\}$, where $\xi$ denotes the average length of stopped trajectories in units of $l$ leapfrog jumps, $\tilde{a}$ the average probability that the symmetry condition is satisfied, and $a$ the mean acceptance probability. Since $\tilde{a}$ can be close to unity and $c_{U}$ is often smaller than $c_{\pi}$ or $c_{\triangledown}$ in practice, the computational cost of spNUTS2 per distinct sample is comparable to that of the NUTS. However, the overall numerical efficiency of spNUTS2 can be higher because the average distance between the current and the next state of the Markov chain can be larger.

The proof of the following proposition is also given in Appendix E.

Proposition 4**.**

The Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by the sequential-proposal No-U-Turn sampler of type 2 (spNUTS2, Algorithm 7) is reversible with respect to the target density $\bar{\pi}$.

4.3 Adaptive tuning of parameters in HMC

Adaptively tuning parameters in MCMC algorithms using the history of the Markov chain can often lead to enhanced numerical efficiency (Haario et al., 2001; Andrieu and Thoms, 2008). Here we discuss adaptive tuning of some parameters in HMC algorithms. As discussed in Section 4.1, tuning the leapfrog step size $\epsilon$ is one of the critical decisions to make in running HMC. Numerical efficiency of HMC algorithms can be increased by targeting an average acceptance probability that is away from both zero and one (Beskos et al., 2013). Since the mean acceptance probability tends to increase with decreasing step size, we use the following recursive formula to update the step size,

[TABLE]

where $\epsilon_{i}$ and $a_{i}$ denote the leapfrog step size and the acceptance probability of a proposal at the $i$-th iteration, and $a^{*}$ the target mean acceptance probability. We follow a standard approach for the sequence of adaptation sizes by taking $\alpha\,{\in}\,(0,1]$ and $\lambda\,{>}\,0$ (Andrieu and Thoms, 2008).

Tuning the covariance $C$ of the velocity distribution can also increase the numerical efficiency of HMC algorithms. If the marginal distributions for the target density $\bar{\pi}$ along different directions have orders of magnitude differences in standard deviation, the size of leapfrog jumps should typically be on the order of the smallest standard deviation in order to avoid numerical instability (Neal, 2011). In this case a large number of leapfrog jumps are needed to make a global move in the direction having the largest standard deviation. Figure 4(a) shows an example of leapfrog trajectory when the target distribution is $\mathcal{N}(0,\Sigma)$, where the standard deviation of one principal component of $\Sigma$ is twenty times larger than the other. In this diagram, the leapfrog trajectory takes about 120 jumps to move across the less constrained direction from one end of a level set to the other. Choosing a covariance $C$ for the velocity distribution close to the covariance of the target distribution can substantially reduce the number of leapfrog jumps needed to explore the sample space in every direction (Neal, 2011, Section 5.4.1). Figure 4(b) shows that a leapfrog trajectory can loop around the level set with only fifteen jumps when the covariance of the velocity distribution $C$ is equal to $\Sigma$. We note that the covariance $C$ affects not only the velocity distribution and the leapfrog updates, but also the U-turn condition (23) for NUTS-type algorithms (the original NUTS, spNUTS1, and spNUTS2) via the cosAngle function.

For adaptive tuning, the covariance $C_{i}$ used at the $i$-th iteration can be set equal to the sample covariance of the Markov chain sampled up to the previous iteration. During initial iterations, a fixed covariance $C_{0}$ can be used to avoid numerical instability (Haario et al., 2001):

[TABLE]

It is possible to take $C_{i}$ as a diagonal matrix whose diagonal entries are given by the sample marginal variances of each component of the Markov chain (Haario et al., 2005). This approach is effective when the components have different scales. The computational cost can be substantially reduced by using a diagonal covariance matrix when the target distribution is high dimensional, because operations such as the Cholesky decomposition of $C_{i}$ can be avoided. Marginal sample variances can be updated with little overhead at each iteration using a recursive formula.

4.4 Numerical examples

4.4.1 Multivariate normal distribution

We used two examples to study the numerical efficiency of various algorithms discussed in this paper. We first considered a one hundred dimensional normal distribution $\mathcal{N}(0,\Sigma)$ where the covariance matrix $\Sigma$ is diagonal and the marginal standard deviations form a uniformly increasing sequence from 0.01 to 1.00.

We compared the numerical efficiency of the following five algorithms, first without adaptively tuning the covariance of the velocity distribution $C$: the standard HMC, HMC with sequential proposals (abbreviated as spHMC, which is equivalent to XCGHMC in Algorithm 4), the NUTS algorithm by Hoffman and Gelman (2014) (the efficient version), spNUTS1 (Algorithm 6), and spNUTS2 (Algorithm 7). All experiments were carried out using the implementation of the algorithms in R (R Core Team, 2018). The source codes are available at https://github.com/joonhap/spMCMC. The covariance matrix $C$ of the velocity distribution was set equal to the one hundred dimensional identity matrix. The leapfrog step size $\epsilon$ was adaptively tuned using (25) with $\alpha\,{=}\,0.7$ and target acceptance probabilities $a^{*}$ varying from 0.45 to 0.95. The adaptation started from the one hundredth iteration. The acceptance probability at the $i$-th iteration $a_{i}$ was computed using the state that was one leapfrog jump away from the current state of the Markov chain to ensure that the leapfrog jump size $\epsilon$ converges to the same value for the same target acceptance probability across the various algorithms. When running HMC, spHMC, spNUTS1, and spNUTS2, the leapfrog step size was randomly perturbed at each iteration by multiplying to $\epsilon_{i}$ a uniform random number between 0.8 and 1.2. Randomly perturbing the leapfrog step size can improve the mixing of the Markov chain constructed by HMC algorithms (Neal, 2011). For the NUTS, we found perturbing the leapfrog step size did not improve numerical efficiency and thus used $\epsilon_{i}$. In HMC and spHMC, each proposal was obtained by making fifty leapfrog jumps. In spHMC, a maximum of $N\,{=}\,10$ proposals were tried in each iteration and the first acceptable proposal was taken as the next state of the Markov chain (i.e., $L\,{=}\,1$). In spNUTS1 and spNUTS2, the stopping condition was checked according to the schedule $b_{j}\,{=}\,2^{j-1}$ for $j\in 1{\hskip 1.70717pt:\hskip 1.70717pt}j_{\max}$ with $j_{\max}\,{=}\,15$, and the unit number of leapfrog trajectories was set to one (i.e., $l\,{=}\,1$ in Algorithms 6 and 7). The value $c$ in the stopping condition (23) was randomly drawn from a uniform$(0,1)$ distribution for each trajectory, as we found randomizing $c$ yielded better numerical results than fixing it at zero. For the NUTS algorithm, randomizing $c$ did not improve the numerical efficiency, so each trajectory was stopped when the cosine angle fell below zero (i.e., $c\,{=}\,0$), as was in Hoffman and Gelman (2014). In spNUTS1, the maximum number of proposals $N$ in each iteration was set to either one or five. In spNUTS2, a maximum of $N\,{=}\,20$ consecutive states on a leapfrog trajectory were tried in each attempt to find an acceptable state. Every algorithm ran for $M\,{=}\,20{,}200$ iterations. As a measure of numerical efficiency, the effective sample size (ESS) of each component of the Markov chain was computed using an estimate of the spectral density at frequency zero via the effectiveSize function in R package coda (Plummer et al., 2006). The first two hundred states of the Markov chains were discarded when computing the effective sample sizes. Each experiment was independently repeated ten times. All computations were carried out using the Boston University Shared Computing Cluster.

Figure 5 shows both the minimum and the average effective sample size for the one hundred variables divided by the runtime in seconds when the covariance $C$ of the velocity distribution was fixed at the identity matrix. We observed that there were large variations in the effective sample sizes among the $d\,{=}\,100$ variables for the Markov chains constructed by HMC and spHMC, resulting in minimum ESSs much smaller than average ESSs. This happened due to the fact that for some variables the leapfrog trajectories with fifty jumps consistently tended to return to states close to the initial positions. The Markov chains mixed slowly in these variables. On the other hand, the leapfrog trajectories tended to reach the opposite side of the level set of Hamiltonian for some variables, for which the autocorrelation at lag one was close to $-1$. For these variables, the effective sample size was greater than the length of the Markov chain $M$. There were much variations in the effective sample size among the variables for the Markov chains constructed by spNUTS1 and spNUTS2 when the stopping cosine angle $c$ was fixed at zero, but variations diminished when $c$ was varied uniformly in the interval $(0,1)$.

The highest value of the minimum ESS per second achieved by spHMC among various values of the target acceptance probability was about fifty percent higher than that by the standard HMC. For this multivariate normal distribution, the number of leapfrog jumps $l\,{=}\,50$ for HMC and spHMC was within the range of the average number of jumps in the leapfrog trajectories constructed by the NUTS, spNUTS1, and spNUTS2 algorithms. Thus the effective sample sizes by HMC and spHMC were comparable to those by the other three algorithms, but the runtimes tended to be shorter. The highest minimum ESS per second by spNUTS1 with $N\,{=}\,5$ was 7.6 times higher than that by the NUTS and 6.9 times higher than that by spNUTS2. The runtimes of the NUTS were more than ten times longer than those of spNUTS1 and twice longer than those of spNUTS2. This happened because the evaluation of the gradient of the log target density took much less computation time than the evaluation of the log target density for this example. The highest minimum ESS per second by spNUTS1 when up to five sequential proposals were made (i.e., $N\,{=}\,5$) was twenty percent higher than when only one proposal was made ($N\,{=}\,1$).

Next we ran the NUTS, spNUTS1, and spNUTS2 algorithms for the same target distribution $\mathcal{N}(0,\Sigma)$ but with adaptively tuning the covariance $C$ of the velocity distribution. The covariance $C_{i}$ at the $i$-th iteration for $i\,{\geq}\,100$ was set to a diagonal matrix whose diagonal entries were given by the marginal sample variances of the Markov chain constructed up to that point. We did not test HMC or spHMC when $C$ was adaptively tuned, because the leapfrog step size, and thus the total length of the leapfrog trajectory with a fixed number of jumps, varied depending on the tuned values for $C$. Figure 6 shows the minimum and average ESS among $d\,{=}\,100$ variables divided by the runtime in seconds. The highest minimum ESS per second improved more than fifty times, compared to when the covariance $C$ was fixed, for the NUTS. There was more than a five-fold improvement for spNUTS1, and more than a 25-fold improvement for spNUTS2. The highest minimum ESS per second by the NUTS was $19\%$ higher than that by spNUTS1 ($N\,{=}\,5$) and $86\%$ higher than that by spNUTS2. The NUTS was relatively more efficient when $C$ was adaptively tuned because the trajectories were built using fewer leapfrog jumps. The computational advantage of spNUTS1 that the log target density is not evaluated at every leapfrog jump is relatively small when there are only few jumps per trajectory. When $C$ is close to $\Sigma$, the sampling task is essentially equivalent to sampling from the standard normal distribution, in which case larger leapfrog step sizes may be used. The trajectories were made of five to eight leapfrog jumps at the most efficient target acceptance probability when $C$ was adaptively tuned. In comparison, the number of leapfrog jumps in a trajectory was between 80 and 250 when $C$ was not adaptively tuned.

4.4.2 Bayesian logistic regression model

We also experimented the numerical efficiency of the NUTS, spNUTS1, and spNUTS2 using the posterior distribution for a Bayesian logistic regression model. The Bayesian logistic regression model and the data we used are identical to those considered by Hoffman and Gelman (2014). The German credit dataset from the UCI repository (Dua and Graff, 2017) consists of twenty four attributes of individuals and one binary variable classifying those individuals’ credit. The posterior density is proportional to

[TABLE]

where $x_{i}$ denotes the twenty four dimensional covariate vector for the $i$-th individual and $y_{i}$ denotes the classification result taking a value from $\pm 1$. We did not normalize the covariates to zero mean and unit variance as in Hoffman and Gelman (2014), because we let $C$ be adaptively tuned. The covariance $C$ was set to a diagonal matrix having as its diagonal entries the marginal sample variances of the constructed Markov chain up to the previous iteration. All algorithms were run under the same settings as those used for the multivariate normal distribution example.

Figure 7 shows the minimum and average ESS across $d\,{=}\,25$ variables per second of runtime. The minimum ESS per second by spNUTS1 at the most efficient target acceptance probability was 2.6 times higher than that by the NUTS and 1.7 times higher than that by spNUTS2. The differences in the numerical efficiency was led mostly by the differences in the runtime. The numbers of leapfrog jumps in stopped trajectories tended to be larger than those for the normal distribution example due to the correlations between the variables in this Bayesian logistic regression model; the numbers of leapfrog jumps were about fifty for the NUTS, twenty seven for spNUTS1, and twenty two for spNUTS2.

5 Conclusion

The sequential-proposal MCMC framework is readily applicable to a wide range of MCMC algorithms. The flexibility and simplicity of the framework allow for various adjustments to the algorithms and offer possibilities of developing new ones. In this paper, we showed that the numerical efficiency of MCMC algorithms can be improved by using sequential proposals. In particular, we developed two novel NUTS-type algorithms, which showed higher numerical efficiency than the original NUTS by Hoffman and Gelman (2014) on two examples we examined. In Appendix F, we apply the sequential-proposal framework to the bouncy particle sampler (BPS) and demonstrate an advantageous property that the sequential-proposal BPS can readily make jumps between multiple modes. The possibilities of other applications of the sequential-proposal MCMC framework can be explored in future research.

Acknowledgement

This work was supported by National Science Foundation grants DMS-1513040 and DMS-1308918. The authors thank Edward Ionides, Aaron King, and Stilian Stoev for comments on an earlier draft of this manuscript. The authors also thank Jesús María Sanz-Serna for informing us about related references.

Appendix A Proof of detailed balance for Algorithm 2 (sequential-proposal Metropolis-Hastings algorithm)

Here we give a proof that Algorithm 2 constructs a reversible Markov chain with respect to the target density $\bar{\pi}$. In what follows, we denote the $l$-th rank of a given finite sequence $a_{n:m}$ by $r_{l}(a_{n:m})$; that is, if we reorder the sequence $a_{n:m}$ as $a_{(1)}\geq a_{(2)}\geq\cdots\geq a_{(m-n+1)}$, then $r_{l}(a_{n:m})=a_{(l)}$. If $l$ is greater than the length of the sequence $a_{n:m}$, we define $r_{l}(a_{n:m}):=0$. We also define $r_{0}(a_{n:m}):=\infty.$

Proposition 5**.**

The Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 2 is reversible with respect to the target density $\bar{\pi}$.

Proof.

It suffices to show the claim for fixed $N$ and $L$. The general case immediately follows by considering a mixture over $N$ and $L$ according to $\nu(N,L)$.

We will show that for a given $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, the probability density of taking $y_{n}$ as the next state of the Markov chain starting from the current state $y_{0}$ after rejecting a sequence of proposals $y_{1:n-1}$ is the same as the probability density of taking $y_{0}$ starting from $y_{n}$ after going through a reversed sequence of proposals $y_{n-1:1}$. The case for $n\,{=}\,1$ coincides with a standard Metropolis-Hastings algorithm. We now fix $n\,{\geq}\,2$. Denoting the uniform$(0,1)$ random variable drawn at the beginning of the iteration by $\Lambda$, the $k$-th proposal $y_{k}$ is considered acceptable if and only if

[TABLE]

Multiplying both the numerator and the denominator by $\prod_{j=k+1}^{n}q(y_{j}{\,|\,}y_{j-1})$, we see that the above condition is equivalent to

[TABLE]

For $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$, we define the following quantities

[TABLE]

such that the condition (26) can be concisely written as

[TABLE]

In what follows, $p_{k}(y_{0},y_{1},\dots,y_{n})$ will be denoted by $p_{k}$ for brevity. The proposal $y_{n}$ is taken as the next state of the Markov chain if and only if it is the $L$-th acceptable proposal among the sequence of proposals $y_{1:n}$. This happens if and only if $\Lambda<p_{n}/p_{0}$ and there are exactly $L{-}1$ proposals among $y_{1:n-1}$ such that $\Lambda<p_{k}/p_{0}$. The latter condition is satisfied if and only if $\Lambda$ is less than the $L{-}1$-th largest number among $p_{1:n-1}/p_{0}$ but greater than or equal to the $L$-th largest number among the same sequence, that is,

[TABLE]

Under the assumption that $X^{(i)}$ is distributed according to the target density $\bar{\pi}$, the probability that the current state of the Markov chain is in a set $A\,{\in}\,\mathcal{X}$ and the $n$-th proposal, which is in a set $B\,{\in}\,\mathcal{X}$, is taken as the next state of the Markov chain is given by

[TABLE]

We will change the notation of dummy variables by writing $y_{0}\leftarrow y_{n}$, $y_{1}\leftarrow y_{n-1}$, $\dots$, $y_{n}\leftarrow y_{0}$. But note that $p_{k}(y_{n},y_{n-1},\dots,y_{0})$ can be expressed as

[TABLE]

which is the same as an expression for $p_{n-k}(y_{0},y_{1},\dots,y_{n})$. Thus under the change of notation, (27) can be re-written as

[TABLE]

where $p_{k}$ denotes $p_{k}(y_{0},y_{1},\dots,y_{n})$ for $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$. The above integral is equal to (27) with the sets $A$ and $B$ interchanged. Thus we have proved that the probability that the current state of the Markov chain is in $A$ and the $n$-th proposal, which is in $B$, is taken as the next state of the Markov chain is equal to the probability that the current state is in $B$ and the $n$-th proposal, which is in $A$, is taken as the next state. Summing the established equality over all $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and finally noting that the next state in the Markov chain is set equal to the current state in the case where there were less than $L$ proposals found in the first $N$ proposals, we reach the conclusion that under the assumption that $X^{(i)}$ is distributed according to $\bar{\pi}$, the probability that the current state of Markov chain is in $A$ and the next state is in $B$ is the same as the probability that the current state is in $B$ and the next state is in $A$. This finishes the proof of detailed balance for Algorithm 2. ∎

Appendix B Sequential-proposal Metropolis-Hastings algorithms with proposal kernels dependent on previous proposals

In Section 2.2 we presented a generalization of Algorithm 2 in which the proposal kernel can depend on previous proposals made in the same iteration. A pseudocode for this generalized version is given in Algorithm 8. A proof that Algorithm 8 constructs a reversible Markov chain with respect to the target density $\bar{\pi}$ is given below.

Proposition 6**.**

Algorithm 8 constructs a reversible Markov chain with respect to the target density $\bar{\pi}$.

Proof.

Again we consider fixed $N$ and $L$, because the general case can easily follow by considering a mixture over $N$ and $L$. Let $\Lambda$ denote the uniform$(0,1)$ number drawn at the start of the iteration. We will denote the value of current state of the Markov chain by $y_{0}$, and the values of a sequence of proposals up to the $n$-th proposal as $y_{1}$, $\dots$, $y_{n}$. The $n$-th proposal $y_{n}$ is taken as the next state of the Markov chain if and only if

[TABLE]

and there are exactly $L-1$ numbers $k$ among $1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$ that satisfy

[TABLE]

and also there exist exactly $L-1$ numbers $k^{\prime}$ among $1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$ that satisfy

[TABLE]

The inequality (28) can be expressed as

[TABLE]

We note that the numerator in the expression above is the probability density of drawing a sequence of proposals in the order $y_{k}\,{\to}\,y_{k-1}\,{\to}\,\cdots\,{\to}\,y_{0}\,{\to}\,y_{k+1}\,{\to}\,y_{k+2}\,{\to}\,\dots\,{\to}\,y_{n}$, where the value $y_{j}$ for $j\,{\geq}\,k\,{+}\,1$ is drawn from a proposal density $q_{j}(\cdot{\,|\,}y_{j-1},y_{j-2},\dots,y_{0})$. We denote this probability density by

[TABLE]

We also denote the numerator in (29) by

[TABLE]

which gives the probability density of drawing proposals in the order $y_{k}\to y_{k+1}\to\cdots y_{n}\to y_{k-1}\to\cdots\to y_{0}$, where $y_{j}$ for $j\,{\leq}\,k\,{-}\,1$ is drawn from $q_{n-j}(\cdot{\,|\,}y_{j+1},\dots,y_{n})$. One can easily check the following relations:

[TABLE]

where we remind the reader of our notation $y_{0:n}:=(y_{0},y_{1},\dots,y_{n})$ and $y_{n:0}:=(y_{n},y_{n-1},\dots,y_{0})$. Now (28) and (29) can be concisely expressed as

[TABLE]

respectively. The conditions required for taking $y_{n}$ as the next state of the Markov chain can be summarized by the following inequalities:

[TABLE]

where $r_{L}$ denotes the function returning the $L$-rank as defined in Section A, and $\frac{p_{1:n-1}}{p_{0}}(y_{0:n})$ denotes the sequence of values \big{(}\frac{p_{1}(y_{0:n})}{p_{0}(y_{0:n})},\dots,\frac{p_{n-1}(y_{0:n})}{p_{0}(y_{0:n})}\big{)}. In what follows, $p_{k}(y_{0:n})$ and $\overline{p}_{k}(y_{0:n})$ will be written as $p_{k}$ and $\overline{p}_{k}$ for brevity. Under the assumption that at the current iteration the state of the Markov chain is distributed according to $\bar{\pi}$, the probability that the current state is in $A$ and the $n$-th proposal, which is in $B$, is taken as the next state of the Markov chain is given by

[TABLE]

We now change the notation of dummy variables by writing $y_{0}\leftarrow y_{n}$, $y_{1}\leftarrow y_{n-1}$, $\dots$, $y_{n}\leftarrow y_{0}$, and noting the relations (30), we may rewrite (31) as

[TABLE]

But the above display is equal to what is obtained when the sets $A$ and $B$ are interchanged in (31). Thus we have proved that, denoting the current state of the Markov chain as $X^{(i)}$ and the next state as $X^{(i+1)}$ and assuming that $X^{(i)}$ is distributed according to $\bar{\pi}$,

[TABLE]

Summing the above equation for $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and considering that $X^{(i+1)}$ is set equal to $X^{(i)}$ for all scenarios except when a proposal among $y_{0}$, $\dots$, $y_{N}$ is taken as the next state of the Markov chain, we reach the conclusion that

[TABLE]

which shows the desired detailed balance of the Markov chain with respect to $\bar{\pi}$. ∎

Appendix C Equivalence between the sequential-proposal Metropolis-Hastings algorithm and the delayed rejection method (only in the case where proposals are not path-dependent)

Before showing the equivalence between sequential-proposal Metropolis-Hastings algorithms for $L\,{=}\,1$ and the delayed rejection method when the proposal kernel is not path-dependent, we briefly check that the target density $\bar{\pi}$ is invariant in the delayed rejection method. It suffices to check that the detailed balance equation holds for each $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$. The probability density that $y_{0}$ is drawn from $\bar{\pi}$, the proposals $y_{1},\dots,y_{n}$ are drawn sequentially, and $y_{n}$ is the first accepted value is given by

[TABLE]

Since the above quantity is symmetric with respect to reversing the order of sequence from $y_{0:n}$ to $y_{n:0}$, it also equals the probability density that starting from $y_{n}$, the proposals $y_{n-1},\dots,y_{1},y_{0}$ are drawn and $y_{0}$ becomes the first accepted value, which is given by

[TABLE]

Combining the case for $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, we see that the Markov chain constructed by the delayed rejection method is reversible with respect to the target density $\bar{\pi}$. We now prove the following proposition.

Proposition 7**.**

Sequential-proposal Metropolis-Hastings algorithm (Algorithm 2) for $L\,{=}\,1$ and fixed $N$ constructs a Markov chain that has the same law as that constructed by the delayed rejection method.

Proof.

For both Algorithm 2 and the delayed rejection method, let the state of the Markov chain at the start of a certain iteration be denoted by $y_{0}$. In both algorithms, the probability density of drawing $y_{0}$ from $\bar{\pi}$ and a sequence of proposals $y_{1:N}$ using a proposal kernel with density $q$ is given by

[TABLE]

Given the current state $y_{0}$ and a sequence of proposals $y_{1:n}$, the probability that $y_{n}$ is taken as the next state of the Markov chain is obtained by subtracting the probability that all of $y_{1:n}$ are rejected from the probability that $y_{1:n-1}$ are rejected. Thus it suffices to show that for an arbitrary $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and a given sequence of proposals $y_{1:n}$, the probability that all of the proposals are rejected is the same in both algorithms. In what follows we assume $y_{1:n}$ are drawn and fixed. We define for $k\in 1{\hskip 1.70717pt:\hskip 1.70717pt}N$,

[TABLE]

For brevity, $c_{k}(y_{0:n})$ will be denoted simply by $c_{k}$ in cases where the argument $y_{0:n}$ can be clearly understood. In Algorithm 2, all $y_{1:n}$ are rejected if and only if $\Lambda\geq c_{k}$ for all $k\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}n$ where $\Lambda\sim\text{unif}(0,1)$ is the random number drawn at the start of the current iteration, so the probability that all $y_{1:n}$ are rejected is given by

[TABLE]

For the delayed rejection method, we denote

[TABLE]

Since the probability that all $y_{1:n}$ are rejected in an implementation of the delayed rejection method is given by $\prod_{k=1}^{n}(1-\beta_{k})$, our goal is to show that

[TABLE]

which we will prove by induction. The case for $n=1$ is obvious from the definition of $\beta_{1}$. Suppose we have

[TABLE]

We denote for $k\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$ and $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$

[TABLE]

In (33), both $\beta_{1:n-1}$ and $c_{1:n-1}$ are functions of $y_{0:n}$. If we change the notation by writing $y_{0}\leftarrow y_{n}$, $y_{1}\leftarrow y_{n-1}$, $\dots$, $y_{n}\leftarrow y_{0}$, the equation (33) becomes

[TABLE]

It can also be easily checked from the definitions of $c_{k}$ and $\bar{c}_{k}^{n}$ that $c_{n}\bar{c}^{n}_{k}=c_{n-k}$ for $k\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$. Since the acceptance probability $\beta_{k}$ in the delayed rejection method is given by

[TABLE]

we observe that

[TABLE]

using the induction hypothesis (33) and (34). For real numbers $u$, $v$, and $w$, the following relation holds:

[TABLE]

We see that

[TABLE]

Thus from (35), showing (32) reduces to checking

[TABLE]

However, one can check the following relation holds:

[TABLE]

from which (36) follows by letting $u=c_{n}$, $w=\max(c_{1:n-1})$, and $v=1$. ∎

Appendix D Proof of Proposition 2

In this section, we will show that sequential-proposal MCMC algorithms using deterministic kernels (Algorithm 3) construct reversible Markov chains with respect to the target density $\bar{\pi}$. We first present the following lemma.

Lemma 1**.**

Suppose (11) and (13) hold. Define recursively $\mathcal{S}_{\tau}^{n}:=\mathcal{S}_{\tau}^{n-1}\circ\mathcal{S}_{\tau}$ where $\mathcal{S}_{\tau}^{1}=\mathcal{S}_{\tau}$. Then for any $n\geq 1$, we have $\mathcal{T}\circ\mathcal{S}_{\tau}^{n}\circ\mathcal{T}\circ\mathcal{S}_{\tau}^{n}=\mathcal{I}.$ Moreover, $\mathcal{S}_{\tau}$ is a bijective map.

Proof.

From (11), we have $\mathcal{T}\circ\mathcal{T}=\mathcal{I}$. Thus from (13), we have $\mathcal{T}=\mathcal{T}\circ\mathcal{T}\circ\mathcal{S}_{\tau}\circ\mathcal{T}\circ\mathcal{S}_{\tau}=\mathcal{S}_{\tau}\circ\mathcal{T}\circ\mathcal{S}_{\tau}$. Thus, we can see that $\mathcal{S}_{\tau}^{n}\circ\mathcal{T}$ is a self-inverse for any $n\geq 1$ from induction

[TABLE]

It also follows that $\mathcal{S}_{\tau}\circ\mathcal{T}\circ\mathcal{S}_{\tau}\circ\mathcal{T}=\mathcal{T}\circ\mathcal{T}=\mathcal{I}$. Thus, since $f\circ g=\mathcal{I}$ implies that function $f$ is surjective and $g$ is injective, the relation $\mathcal{S}_{\tau}\circ(\mathcal{T}\circ\mathcal{S}_{\tau}\circ\mathcal{T})=\mathcal{I}$ implies that $\mathcal{S}_{\tau}$ is surjective and $(\mathcal{T}\circ\mathcal{S}_{\tau}\circ\mathcal{T})\circ\mathcal{S}_{\tau}=\mathcal{I}$ implies that $\mathcal{S}_{\tau}$ is injective. ∎

Proposition 2.

The extended target distribution with density $\Pi(x,v)$ is a stationary distribution for the Markov chain $\left(X^{(i)},V^{(i)}\right)_{i\in 1:M}$ constructed by the sequential-proposal MCMC algorithm using a deterministic kernel (Algorithm 3). Furthermore, the Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 3, marginally for the $x$-component, is reversible with respect to the target distribution $\bar{\pi}(x)$.

Proof.

We will prove the claim for the case where $N$, $L$, and $\tau$ are fixed, since the general case can easily follow by considering a mixture over these parameters. In Algorithm 3, the $n$-th proposal $(Y_{n},W_{n})$ is obtained as $\mathcal{S}_{\tau}^{n}(Y_{0},W_{0})$, and if there are less than $L$ acceptable proposals in the first $N$ proposals, the next state of the Markov chain $(X^{(i+1)},V^{(i+1)})$ is set to $(X^{(i)},\mathcal{R}_{X^{(i)}}V^{(i)})$. Each iteration of Algorithm 3 can thus be understood a composition of two operations, where the first operation is simply reflecting the velocity component from $(X^{(i)},V^{(i)})$ to $(X^{(i)},\mathcal{R}_{X^{(i)}}V^{(i)})$, and the second operation proposes a sequence of proposals $(Y_{n},W_{n})=\mathcal{S}_{\tau}^{n}\circ\mathcal{T}(X^{(i)},\mathcal{R}_{X^{(i)}}V^{(i)})$ until $L$ acceptable proposals are found or until $N$ proposals have been made. The reason that we view the algorithm this way is to use the fact that both $\mathcal{T}$ and $\mathcal{S}_{\tau}^{n}\circ\mathcal{T}$ are self-inverse maps. If both the first and the second operations preserve $\Pi$ as an invariant density, the Markov chain constructed by Algorithm 3 preserves $\Pi$ as an invariant density. In fact, we will show that both the first and the second operations satisfy detailed balance with respect to $\Pi$. However, we note that this does not imply that the constructed Markov chain satisfy detailed balance with respect to $\Pi$, because carrying out the first and then the second operation is not the same as carrying out the second and then the first.

It is rather straightforward to see that velocity reflection operation $\mathcal{T}$ establishes detailed balance with respect to $\Pi$. Supposing that $(X,V)\sim\Pi$, we have for $A$, $B$ measurable in $\mathbb{X}\times\mathbb{V}$,

[TABLE]

Upon denoting $v^{\prime}:=R_{x}v$, we can express the right hand side as

[TABLE]

where we have used the condition (12). This shows that $\mathcal{T}$ establishes detailed balance with respect to $\Pi$.

Now we will show that the second operation also establishes detailed balance with respect to $\Pi$. If the current state in the Markov chain $(X^{(i)},V^{(i)})$ is denoted by $(Y_{0},W_{0})$, the second operation starts at $(Y_{0},\mathcal{R}_{Y_{0}}W_{0})$ since the velocity was reflected by the first operation. We will show that for arbitrary $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and for measurable subsets $A$, $B$ of $\mathbb{X}\times\mathbb{V}$,

[TABLE]

provided that $(Y_{0},\mathcal{R}_{Y_{0}}W_{0})$ is distributed according to $\Pi$. Then by combining the cases for $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, we can establish detailed balance for the second operation. For notational convenience, for $(y_{0},w_{0})\,{\in}\,\mathbb{X}\times\mathbb{V}$, we will write $(y_{k},w_{k})=\mathcal{S}_{\tau}^{k}(y_{0},w_{0})$ and $\bar{w}_{k}=\mathcal{R}_{y_{k}}w_{k}$ for $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$. We will write

[TABLE]

for $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$. Note that this definition leads to

[TABLE]

due to (12). Also, since $\mathcal{S}_{\tau}^{n-k}\circ\mathcal{T}(y_{k},\bar{w}_{k})=(y_{n},w_{n})$ and $\mathcal{S}_{\tau}^{n-k}\circ\mathcal{T}$ is a self-inverse map, we have $S_{\tau}^{n-k}\circ T(y_{n},w_{n})=(y_{k},\bar{w}_{k})$. This leads to

[TABLE]

The following steps are similar to corresponding steps in the proof of Proposition 5. We have

[TABLE]

where all functions $p_{k}$, $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$, in the above display take the argument $(y_{0},\bar{w}_{0})$. Using (37), the above equation is equal to

[TABLE]

We change the dummy variables by writing $(y_{0},\bar{w}_{0})\leftarrow(y_{n},w_{n})$. Since $\mathcal{S}_{\tau}^{n}\circ\mathcal{T}(y_{0},\bar{w}_{0})=(y_{n},w_{n})$, we can also write $(y_{n},w_{n})\leftarrow(y_{0},\bar{w}_{0})$. The above display can be re-written as

[TABLE]

which is equal to (38) where the sets $A$ and $B$ are interchanged. Thus we have proved

[TABLE]

By adding the cases for $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, we can conclude the proof of detailed balance for the second operation with respect to $\Pi$. Since both the first and the second operations preserves $\Pi$ as an invariant density, Algorithm 3 preserves $\Pi$ as an invariant density. Finally, refreshing the velocity $V^{(i+1)}$ from $\psi(\cdot{\,;\,}X^{(i+1)})$ at the end of the iteration with an arbitrary probability $p^{\text{ref}}(X^{(i+1)})$ clearly preserves the invariant density $\Pi(x,v)=\bar{\pi}(x)\psi(v{\,;\,}x)$.

In order to prove the claim that the marginally for the $x$-component, the Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 3 is reversible with respect to the target distribution $\bar{\pi}(x)$, we denote the position-velocity pair taken as the next state of the Markov chain at the end of the second operation by $(Y^{\prime},W^{\prime})$. We showed above that when $(Y_{0},W_{0})$ is drawn from $\Pi$, $(Y_{0},\mathcal{R}_{Y_{0}}W_{0})$ is also distributed according to $\Pi$. Due to the fact that the second operation satisfies detailed balance with respect to $\Pi$, we see that for measurable subsets $A$, $B$ of $\mathbb{X}$,

[TABLE]

This shows that the Markov chain $\left(X^{(i)})\right)_{i\in 1:M}$ constructed by Algorithm 3 is reversible with respect to $\bar{\pi}$, which is the marginal distribution of $\Pi$ for the $x$-component. ∎

Appendix E Proofs of detailed balance for sequential-proposal No-U-Turn samplers (spNUTS1 and spNUTS2)

We prove that both spNUTS1 and spNUTS2 algorithms (Algorithms 6 and 7) construct reversible Markov chains with respect to the target distribution $\bar{\pi}$.

Proposition 3.

The Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by the sequential-proposal No-U-Turn sampler of type 1 (spNUTS1, Algorithm 6) is reversible with respect to the target distribution $\bar{\pi}$.

Proof.

We assume that the cosine value $c$ at which trajectory extensions in Algorithm 6 stop is fixed, as the general case readily follows by considering a mixture over $c$. The state of the Markov chain constructed by the algorithm in the current iteration is denoted by $Y_{0}$, and assumed to be distributed according to $\bar{\pi}$. The velocity drawn from $\psi_{C}$ at the start of the iteration is denoted by $W_{0}$. For $k\,{\geq}\,1$, the leapfrog trajectory starting from $(Y_{k-1},W_{k-1})$ stops at $(Y_{k},W^{\prime}_{k})$, and the function that maps the initial position-velocity pair to the final pair will be denoted by $\mathcal{S}$, such that $(Y_{k},W^{\prime}_{k})=\mathcal{S}(Y_{k-1},W_{k-1})$. We will show that for $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$ and for measurable subsets $A$ and $B$ of $\mathbb{X}\,{\times}\,\mathbb{V}$,

[TABLE]

Then, by considering the cases $A=A_{0}\,{\times}\,\mathbb{V}$ and $B=B_{0}\,{\times}\,\mathbb{V}$ for some $A_{0},B_{0}\subset\mathbb{X}$ and summing over $n\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}N$, we show that the Markov chain constructed by Algorithm 6 is reversible with respect to $\bar{\pi}$.

When $(Y_{k},W^{\prime}_{k})$ is rejected, the velocity is refreshed by drawing $U_{k}\sim\psi_{C}$ and setting

[TABLE]

In following equations, $y_{1:n}$, $w_{1:n-1}$, and $w^{\prime}_{1:n}$ will denote functions of $y_{0}$, $w_{0}$, and $u_{1:n-1}$ defined recursively by

[TABLE]

We also denote $\bar{w}_{k}:=-w_{k}$ and $\bar{w}^{\prime}_{k}:=-w^{\prime}_{k}$. The left hand side of (39) is then given by

[TABLE]

To establish a symmetric relationship between $(Y_{0},W_{0})$ and $(Y_{n},\bar{W}^{\prime}_{n})$, we define

[TABLE]

and write $\bar{u}^{\prime}_{k}:=-u^{\prime}_{k}$ for $k\,{\in}\,1{\hskip 1.70717pt:\hskip 1.70717pt}n{-}1$. Due to the symmetric nature of the stopping condition (20), we have

[TABLE]

It is also readily observed from the definition of $u^{\prime}_{k}$ that

[TABLE]

The two relations above form a counterpart to (40) in a symmetric relationship between $(y_{k},w_{k})$ and $(y_{n-k},\bar{w}^{\prime}_{n-k})$ for $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}n$. Furthermore, since $\psi_{C}(v)$ is a function of $\|v\|_{C}$, we have $\psi_{C}(u_{k})=\psi_{C}(\bar{u}^{\prime}_{k})$ and $\psi_{C}(w^{\prime}_{k})=\psi_{C}(\bar{w}_{k})$. Finally, we use the equation

[TABLE]

which is stated as Lemma 2 and proved below. This, together with the fact that leapfrog jumps preserve the volume element, leads to

[TABLE]

Thus, (41) is equal to

[TABLE]

which, under the change of notation of dummy variables $(y_{0},w_{0})\leftarrow(y_{n},\bar{w}^{\prime}_{n})$, becomes

[TABLE]

Since the above expression equals the right hand side of (39), the claim of detailed balance is proved. ∎

Lemma 2**.**

Let $C$ be a positive definite symmetric matrix in $\mathbb{R}^{d\times d}$. Given $(w^{\prime},u)$ in $\mathbb{R}^{2d}$, define $u^{\prime}:=w^{\prime}\frac{\|u\|_{C}}{\|w^{\prime}\|_{C}}$ and $w:=u\frac{\|w^{\prime}\|_{C}}{\|u\|_{C}}$. Then we have

[TABLE]

Proof.

It is sufficient to prove the claim for $C=I$, the identity matrix. To see this, we denote $\tilde{w}:=C^{-1/2}w$, $\tilde{w}^{\prime}:=C^{-1/2}w^{\prime}$, $\tilde{u}:=C^{-1/2}u$, $\tilde{u}^{\prime}:=C^{-1/2}u^{\prime}$, and denote the Euclidean norm as $\|v\|:=\sqrt{v^{T}v}$. Then since $\|\tilde{w}\|^{2}=w^{T}C^{-1}w=\|w\|_{C}^{2}$ and the same kind of relation holds for the other three variables, we have

[TABLE]

and similarly $\tilde{w}=\tilde{u}\frac{\|\tilde{w}^{\prime}\|}{\|\tilde{u}\|}$. But we also have

[TABLE]

from which we see that it is sufficient to prove that the rightmost term is equal to unity.

Now we will assume $C=I$. Computing partial derivatives yields

[TABLE]

We carry out elementary column and row operations as follows to obtain

[TABLE]

[TABLE]

The absolute value of the determinant of the above matrix can be directly computed as

[TABLE]

∎

Proposition 4.

The Markov chain $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by the sequential-proposal No-U-Turn sampler of type 2 (spNUTS2, Algorithm 7) is reversible with respect to the target distribution $\bar{\pi}(x)$.

Proof.

We assume the cosine angle $c$ is fixed. Let $Y_{0}$ denote the current state of the Markov chain constructed by Algorithm 7 and $W_{0}$ the velocity drawn from $\psi_{C}$ at the start of the current iteration. We recursively let $(Y_{k},W_{k})$ denote the state reached after $l$ leapfrog jumps starting from $(Y_{k-1},W_{k-1})$, for $k\,{\geq}\,1$. We consider the case where the trajectory stops at the $b_{j}$-th acceptable state which is equal to $(Y_{k^{*}},W_{k^{*}})$ for some $k^{*}$. We will consider the probability

[TABLE]

for measurable subsets $A$ and $B$ of $\mathbb{X}\,{\times}\,\mathbb{V}$. Let $\mathcal{K}$ be a subset of $0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}$ and write $\mathcal{K}^{c}=(0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*})\,{\setminus}\,\mathcal{K}$. Given the values $(Y_{k},W_{k})=(y_{k},w_{k})$, $k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}$ and the uniform$(0,1)$ random variable $\Lambda$ drawn at the start of the current iteration, the states $\{(y_{k},w_{k}){\,;\,}k\,{\in}\,\mathcal{K}\}$ are deemed acceptable and $\{(y_{k},w_{k}){\,;\,}k\in\mathcal{K}^{c}\}$ not acceptable if and only if the quantity

[TABLE]

is equal to unity. Note that we consider $(y_{0},w_{0})$ as an acceptable state here. Integrating the above quantity over $\Lambda$, we see that the probability of finding $\{(y_{k},w_{k}){\,;\,}k\in\mathcal{K}\}$ acceptable and $\{(y_{k},w_{k}){\,;\,}k\in\mathcal{K}^{c}\}$ not acceptable is given by

[TABLE]

We will consider $\mathcal{K}$ that satisfies the following three conditions: $\{0,k^{*}\}\subset\mathcal{K}\subset 0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}$, $|\mathcal{K}|=b_{j}\,{+}\,1$, and $\mathcal{K}(i)-\mathcal{K}(i{\,-\,}1)\leq N$ for all $i\in 1{\hskip 1.70717pt:\hskip 1.70717pt}b_{j}$, where the elements of $\mathcal{K}$ are ordered as $0=\mathcal{K}(0)<\mathcal{K}(1)<\cdots<\mathcal{K}(b_{j})=k^{*}$. The last condition is related to the fact that spNUTS2 tries at most $N$ consecutive states to find each new acceptable state.

Given $\{(y_{k},w_{k}){\,;\,}k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}\}$ and $\mathcal{K}$, we let \Upsilon\big{(}\{(y_{k},w_{k}){\,;\,}k\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}\},\mathcal{K}\big{)} denote the indicator function that takes the value of unity if and only if the pair $(y_{k^{*}},w_{k^{*}})$ is the first position-velocity pair satisfying the U-turn condition among $\{(y_{\mathcal{K}(b_{j^{\prime}})},w_{\mathcal{K}(b_{j^{\prime}})}){\,;\,}{j^{\prime}}\,{\in}\,0{\hskip 1.70717pt:\hskip 1.70717pt}j\}$ and the stopped trajectory satisfies the symmetry condition. That is, we define

[TABLE]

where the dependence of cosAngle on $C$ is suppressed. Thus the probability of drawing $(y_{0},w_{0})$ and taking the $b_{j}$-th acceptable pair $(y_{k^{*}},w_{k^{*}})$ as the next state of the Markov chain while finding $\{(y_{k},w_{k}){\,;\,}k\,{\in}\,\mathcal{K}\}$ acceptable and $\{(y_{k},w_{k}){\,;\,}k\,{\in}\,\mathcal{K}^{c}\}$ not acceptable is given by

[TABLE]

[TABLE]

We now consider a reverse scenario where we draw $y_{k^{*}}$ from $\bar{\pi}$ and $-w_{k^{*}}$ from $\psi_{C}$ at the start of the current iteration of the algorithm. The new leapfrog trajectory exactly reverses the original trajectory and is given by $\{(y_{k^{*}-k},-w_{k^{*}-k}){\,;\,}k\in 0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}\}$. We denote $\bar{y}_{k}:=y_{k^{*}-k}$ and $\bar{w}_{k}:=-w_{k^{*}-k}$ for $k\in 0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}$. We also write $\bar{\mathcal{K}}:=\{k^{*}-k{\,;\,}k\in\mathcal{K}\}$ and $\bar{\mathcal{K}}^{c}:=(0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*})\,{\setminus}\,\bar{\mathcal{K}}$. The elements of $\bar{\mathcal{K}}$ will also be denoted as $0=\bar{\mathcal{K}}(0)<\bar{\mathcal{K}}(1)<\cdots<\bar{\mathcal{K}}(b_{j})=k^{*}$. It is easy to see that $\bar{\mathcal{K}}(i)=k^{*}-\mathcal{K}(b_{j}-i)$ for $i\in 0{\hskip 1.70717pt:\hskip 1.70717pt}b_{j}$. It follows that

[TABLE]

From this, the following equations hold:

[TABLE]

Thus we see

[TABLE]

Finally, due to the measure-preserving property of leapfrog jumps, we have $dy_{0}\,dw_{0}=dy_{k^{*}}\,dw_{k^{*}}$. From these facts, we see that the probability of drawing $\bar{y}_{0}$ from $\bar{\pi}$, drawing $\bar{w}_{0}$ from $\psi_{C}$, finding the states $\{(\bar{y}_{k},\bar{w}_{k}){\,;\,}k\in\bar{\mathcal{K}}\}$ acceptable and $\{(\bar{y}_{k},\bar{w}_{k}){\,;\,}k\in\bar{\mathcal{K}}^{c}\}$ not acceptable, and taking $(\bar{y}_{k^{*}},\bar{w}_{k^{*}})$ as the next state of the Markov chain equals

[TABLE]

[TABLE]

[TABLE]

which is the same as (42). Recall that (42) gives the probability corresponding to the case where the index set of the acceptable states is given by $\mathcal{K}$. Let $\mathbb{K}(j,k^{*})$ be the set of index sets $\mathcal{K}$ that satisfy the following three conditions: $\{0,k^{*}\}\subset\mathcal{K}\subset 0{\hskip 1.70717pt:\hskip 1.70717pt}k^{*}$, $|\mathcal{K}|=b_{j}+1$, and $\mathcal{K}(i)-\mathcal{K}(i{\,-\,}1)\leq N$ for all $i\in 1{\hskip 1.70717pt:\hskip 1.70717pt}b_{j}$. Due to the symmetric nature of these conditions, we have $\mathcal{K}\in\mathbb{K}(j,k^{*})$ if and only if $\bar{\mathcal{K}}\in\mathbb{K}(j,k^{*})$. Denoting

[TABLE]

we have from (43),

[TABLE]

Thus

[TABLE]

The third equality follows from the symmetry in $\mathbb{K}(j,k^{*})$. Summing the above equation for $j{\,\geq\,}0$ and $k^{*}{\,\geq\,}1$ and taking $A=A_{0}\,{\times}\,\mathbb{V}$ and $B=B_{0}\,{\times}\,\mathbb{V}$ for some measurable subsets $A_{0}$ and $B_{0}$ of $\mathbb{X}$, we establish detailed balance for the Markov chains $\left(X^{(i)}\right)_{i\in 1:M}$ constructed by Algorithm 7. ∎

Appendix F Connection to the bouncy particle sampler

Recently, a non-reversible, piecewise deterministic MCMC sampling method called the bouncy particle sampler (BPS) has been proposed (Peters et al., 2012; Bouchard-Côté et al., 2018). The BPS constructs a rejection free, continuous time Markov chain. A key advantage of the BPS is that it allows for local updates of the target variables, meaning that the algorithm can update one subset of the target variables at a time while the rest of the variables evolve according to a flow that is easy to compute. However, this algorithm has a limitation in terms of the target distributions it can be used for, because the user needs to be able to draw the arrival times of a non-homogeneous Poisson process which have a rate depending on the gradient of the target density. In this section, we present a new, discrete time version of BPS, which is not rejection free. This discrete time BPS is readily applicable to any target distribution with evaluable unnormalized density. We note that there exists an alternative discrete time BPS, which was given in Vanetti et al. (2017, Algorithm 4).

We describe the discrete time BPS algorithm we propose within the framework of MCMC algorithms using deterministic kernels (see Section 3). We assume that the sample space $\mathbb{X}$ and the velocity space $\mathbb{V}$ are both given by $\mathbb{R}^{d}$. As in Section 3, the density of velocity distribution is denoted by $\psi(v{\,;\,}x)$. We suppose that for each $x\,{\in}\,\mathbb{X}$, there exists a linear operator $\mathcal{R}_{x}\,{:}\,\mathbb{V}\,{\to}\,\mathbb{V}$, which satisfies $\mathcal{R}_{x}\circ\mathcal{R}_{x}=\mathcal{I}$ and $\psi(\mathcal{R}_{x}v{\,;\,}x)=\psi(v{\,;\,}x)$. The time evolution map $\mathcal{S}_{\tau}$ is defined as

[TABLE]

The self-inverse property of the linear operator $\mathcal{R}_{x}$ ensures that the absolute determinant of $\mathcal{R}_{x}$ when viewed as a matrix is equal to unity. This translates into unit Jacobian determinant of the map $\mathcal{R}_{x}$, so the condition (12), $\frac{\psi(\mathcal{R}_{x}v{\,;\,}x)}{\psi(v{\,;\,}x)}\left|\frac{\partial\mathcal{R}_{x}v}{\partial v}\right|=1$, is satisfied. It can also be readily checked that the condition (13) is satisfied for $\mathcal{S}_{\tau}$ and that

[TABLE]

A sequential-proposal discrete time BPS can be obtained as a specific case of Algorithm 3 where $\mathcal{S}_{\tau}$ and $\mathcal{R}_{x}$ are given as above. Algorithm 9 gives a pseudocode for the sequential-proposal discrete time BPS.

A convenient choice for $\psi$ is a multivariate Gaussian density

[TABLE]

where $C$ is a positive definite matrix. In this case, the conditions (11) and (12) hold if and only if

[TABLE]

for a symmetric projection matrix $P$, that is $PP{=}P$ and $P^{T}{=}P$. A matrix $P$ is a symmetric projection matrix in $\mathbb{R}^{d}$ if and only if it is a projection onto the linear span of a subset of an orthonormal basis of $\mathbb{R}^{d}$, that is $P=\sum_{j\in A}e^{j}(e^{j})^{T}$ for some $A\subseteq\{1,2,\dots,d\}$ and some orthonormal basis $(e^{1},\dots,e^{d})$.

A possible choice for $\mathcal{R}_{x}$ includes $-\mathcal{I}$, in which case the proposal map is given by $\mathcal{S}_{\tau}(x,v)=(x+v\tau,v)$. Since the map $\mathcal{S}_{\tau}$ can be readily evaluated, this choice has a computational advantage when multiple sequential proposals are made. Another sensible choice for the operator $\mathcal{R}_{x}$ is the reflection on the hyperplane perpendicular to the gradient of the log target density under the metric given by $C^{-1}$. We write $U(x):=-\log\pi(x)$ and denote the velocity reflection operator by $\mathcal{R}_{\nabla U(x)}$. This velocity reflection operator can be written as

[TABLE]

where $\langle u,w\rangle_{C}:=u^{T}C^{-1}v$. This $\mathcal{R}_{\nabla U(x)}$ is the reflection operator used by the original BPS algorithm of Peters et al. (2012). Since

[TABLE]

we have, in the case where $\tau$ is small such that $\nabla U(x-\mathcal{R}_{x}v\tau)\approx\nabla U(x)$,

[TABLE]

Therefore, repeated application of the map $\mathcal{S}_{\tau}$ has an approximate net effect of moving the particles along the gradient of the log target density.

F.1 Numerical examples

To graphically illustrate the performance of the sequential-proposal discrete time BPS (Algorithm 9), we created a target distribution defined on a unit square. The regions of high likelihood density look like four open rings, or four rotated letters of “C”, as shown in Figure 9. We applied the sequential-proposal discrete time BPS on this model with varying algorithmic parameters. In every experiment, we ran the algorithm up to 120,000 iterations, where the number of acceptable proposals $L$ was fixed at one and the jump size $\tau$ at each iteration varied uniformly between $0.08$ and $0.12$.

Figure 9(a) shows the 120,000 sampled points as black dots. The target density is represented by a color map on a green-white scale at the background. Starting from the initial point, the trajectory of every fourth point is shown by red segments. We varied the maximum number of proposals $N$ from one to ten and twenty. We used the reflection operator $\mathcal{R}_{x}\,{=}\,-\mathcal{I}$ and the velocity refreshment probability $p^{\text{ref}}\equiv 0.1$ for this experiment. In the case of $N\,{=}\,1$, where no subsequent proposals were made if the first proposal was rejected, there was no jumps between “C”s. As we increased $N$ to ten and twenty, the jump between the “C”s happened more frequently, and mixing happened faster. Figure 9(b) shows the distances between the sample means of the constructed Markov chains and the center of the target distribution as a function of runtime in seconds. The sample means clearly did not converge to the true mean when $N\,{=}\,1$, but the sample means converged to the true mean with similar rates when $N\,{=}\,10$ or $N\,{=}\,20$.

Figure 10 shows the same experiment, when the target density was inverted from the original model (i.e., the log target density was multiplied by $-1$). The four “C”s acted as barriers that were difficult for particles to pass through. The velocity reflection operator $\mathcal{R}_{x}\,{=}\,-\mathcal{I}$ was used, and $p^{\text{ref}}$ was 0.1. For $N\,{=}\,1$, there were no jumps across the barriers. For $N\,{=}\,10$ or $20$, however, the jumps across the barriers happened frequently.

Figure 11 shows the numerical results for the original four “C” model for various choices of velocity reflection operator and velocity refreshment probability. The maximum number of proposals $N$ was fixed at twenty. The left column shows the results when the reflection operator $\mathcal{R}_{x}\,{=}\,-\mathcal{I}$ was used. In the middle column, the reflection operator $\mathcal{R}_{\nabla U(x)}$ defined in (44) was used. In the right column, the reflection operator was randomly chosen between $-\mathcal{I}$ and $\mathcal{R}_{\nabla U(x)}$ with equal probability whenever the reflection operator was used by the algorithm. The reflection operator $\mathcal{R}_{x}\,{=}\,-\mathcal{I}$ resulted in a non-ergodic Markov chain when we did not refresh the velocity (top row, left). When $p^{\text{ref}}=0.1$, the choice of $\mathcal{R}_{\nabla U(x)}$ resulted in a slower convergence to the target distribution compared to the case $\mathcal{R}_{x}\,{=}\,-\mathcal{I}$. The speed of convergence was improved if the algorithm used both $-\mathcal{I}$ and $\mathcal{R}_{\nabla U(x)}$ with equal probability. From these results, we see that it is crucial to occasionally refresh the velocity for certain choices of $\mathcal{R}_{x}$ and that using a mixture of different velocity reflection operators can speed up the mixing of the Markov chain.