Speeding up MCMC by delayed acceptance and data subsampling

Let $\theta$ be the parameter in a model with density $p(y_{k}|\theta,x_{k})$, where $y_{k}$ is a potentially multivariate response vector and $x_{k}$ is a vector of covariates for the $k$th observation. Let $l_{k}(\theta)=\log p(y_{k}|\theta,x_{k})$ denote the $k$th observation’s log-density, $k=1,\dots,n$. Given conditionally independent observations, the likelihood function can be written

(2.1)

$$p(y|\theta)=\exp\left[l(\theta)\right],$$

where $l(\theta)=\sum_{k=\text{1}}^{n}l_{k}(\theta)$ is the log-likelihood function. To estimate $p(y|\theta)$, we first estimate $l(\theta)$ based on a sample of size $m$ from the population $\left\{l_{1}(\theta),\dots,l_{n}(\theta)\right\}$ and subsequently use (2.1). The first step corresponds to the classical survey sampling problem of estimating a population total: see Särndal et al., (2003) for an introduction to survey sampling. Note that (2.1) is more general than identically independently distributed (iid) observations, although we require that the log-likelihood can be written as a sum of terms, where each term depends on a unique piece of data information. One example is models with a random effect for each subject, with possibly multiple individual observations per subject (longitudinal data). In this case, a single term in the log-likelihood sum corresponds to the log joint density for a subject, and we therefore sample subjects (rather than individual observations) when estimating (2.1).

A main distinction of sampling schemes is whether the sample is obtained with or without replacement. It is clear that sampling with replacement gives a higher variance for any function of the sample: including the same element more than once does not provide any further information about finite population characteristics. However, it results in sample elements that are independent which facilitates the derivation of Theorem 1 for the delayed acceptance MH. It also allows us to apply the theory and methodology in Quiroz et al., (2016) for the delayed acceptance PMMH in Section 4. It should be noted that the two sampling schemes are approximately the same when $m\ll n$.

Let $u=(u_{1},\dots u_{m})$ be a vector of indices obtained by sampling $m$ indices with replacement from $\{1,\dots,n\}$ and let $\left\{l_{u_{1}}(\theta),\dots,l_{u_{m}}(\theta)\right\}$ be the sample. We will consider Simple Random Sampling (SRS) which means that

$$\Pr(u_{i}=k)=\frac{1}{n}\text{ for }k=1,\dots,n,\text{ and for all }i=1,\dots,m.$$

Payne and Mallick, (2015) use SRS (but without replacement) to form an unbiased estimate of the log-likehood. However, the elements $l_{k}(\theta)$ (for a fixed $\theta$) vary substantially across the population and estimating the total $\sum_{k=1}^{n}l_{k}(\theta)$ with SRS is well known to be very inefficient under such a scenario. Instead $\Pr(u_{i}=k)$ should be (approximately) proportional to a size measure for $l_{k}(\theta)$. Quiroz et al., (2016) argue that this so called Proportional-to-Size sampling is in many cases unlikely to be successful as knowledge of a size measure (which also depends on $\theta$) for all $k$ can often defeat the purpose of subsampling. They instead propose to use SRS but incorporate control variates in the estimator for variance reduction which we now turn to.

The idea in Quiroz et al., (2016) is to homogenize the population $\left\{l_{1}(\theta),\dots,l_{n}(\theta)\right\}$: if the resulting elements are roughly of the same size then SRS is expected to be efficient. Let $q_{k}(\theta)$ denote an approximation of $l_{k}(\theta)$ and decompose

	$\displaystyle l(\theta)$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{n}q_{k}(\theta)+\sum_{k=1}^{n}\left[l_{k}(\theta)-q_{k}(\theta)\right]$
		$\displaystyle=$	$\displaystyle q(\theta)+d(\theta),$

where

$$q(\theta)=\sum_{k=1}^{n}q_{k}(\theta),\quad d(\theta)=\sum_{k=1}^{n}d_{k}(\theta),\quad\text{and}\quad d_{k}(\theta)=l_{k}(\theta)-q_{k}(\theta).$$

We emphasize that all quantities depend on $\theta$ which, from now on, is sometimes suppressed for a compact notation.

Define the random variables $\eta_{i}=nd_{u_{i}}$ and $d_{u_{i}}=l_{u_{i}}-q_{u_{i}}$ with $\mathrm{E}[\eta_{i}]=d$ and

$$\sigma_{\eta}^{2}=\mathrm{V}[\eta_{i}]=n^{2}\mathrm{V}[d_{u_{i}}]=n\sum_{k=1}^{n}(d_{k}-\bar{d})^{2},\quad\text{with }\bar{d}=\sum_{k=1}^{n}d_{k}.$$

The difference estimator estimates $d$ in (2.3) with the Hansen-Hurwitz estimator (Hansen and Hurwitz,, 1943),

$$\hat{d}_{m}=\frac{1}{m}\sum_{i=1}^{m}\eta_{i},\quad\text{with }\mathrm{E}[\hat{d}_{m}]=d\quad\text{and}\quad\mathrm{V}[\hat{d}_{m}]=\frac{\sigma_{\eta}^{2}}{m}.$$

We can obtain an unbiased estimate $\hat{\sigma}_{\eta}^{2}=n^{2}\mathrm{\hat{V}}[d_{u_{i}}]$, where $\hat{\mathrm{V}}[d_{u_{i}}]$ is the usual unbiased sample variance estimator (of the population $\{d_{1},\dots d_{n}\}$). The estimator of the log-likelihood is thus

(2.3)

$$\hat{l}_{m}=\sum_{k=1}^{n}q_{k}(\theta)+\hat{d}_{m},\quad\text{with }\mathrm{E}[\hat{l}_{m}]=l\quad\text{and}\quad\sigma^{2}=\mathrm{V}[\hat{l}_{m}]=\frac{\sigma_{\eta}^{2}}{m}.$$

It also follows that $\hat{\sigma}_{m}^{2}=\hat{\sigma}_{\eta}^{2}/m$ is an unbiased estimator of $\sigma^{2}$.

Quiroz et al., (2016) reason that observations close in data space $(x_{k},y_{k})$ should, for a fixed $\theta$, have similar $l_{k}(\theta)$ values. They cluster the data space into $K$ clusters and approximate, within each cluster, $l_{k}(\theta)$ by a second order Taylor series expansion $q_{k}(\theta)$ around the centroid of the cluster. This allows computing $q$ using $K$ evaluations (instead of $n$), see Quiroz et al., (2016) for details. Bardenet et al., (2015) propose similar control variates, but instead expand with respect to $\theta$ around a reference value $\theta^{\star}$. While it can be shown that $q$ can be computed using a single evaluation, the approximation can be inaccurate when $\theta$ is far from $\theta^{\star}$ giving a large $\sigma^{2}$.

We note that the log-likelihood estimate in Payne and Mallick, (2015) (if with replacement sampling is used) is a special case of (2.3), namely when $q_{k}=0$ for all $k$. We will see that this results in a poor performance of the delayed acceptance algorithm and that control variates are crucial for a successful implementation.

It is clear that an unbiased log-likelihood estimator becomes biased for the likelihood when transformed to the ordinary scale by the exponential function. Since by the Central Limit Theorem (CLT)

$$\sqrt{m}\left(\hat{l}_{m}(\theta)-l(\theta)\right)\rightarrow\mathcal{N}(0,\sigma_{\eta}^{2})\quad\text{as}\quad m\rightarrow\infty,$$

Quiroz et al., (2016) (see also Ceperley and Dewing, 1999; Nicholls et al., 2012) approximately bias corrects the likelihood estimate

(2.4)

$$\hat{p}_{m}(y|\theta,u)=\exp\left(\hat{l}_{m}(\theta)-\hat{\sigma}_{\eta}^{2}(\theta)/2m\right).$$

Equation (2.4) is unbiased if $\sigma_{\eta}^{2}$ is used in place of $\hat{\sigma}_{\eta}^{2}$ (and normality holds for $\hat{l}_{m}$). In practice we need to use $\hat{\sigma}_{\eta}^{2}$ (to not use all data) and Quiroz et al., (2016) show that (2.4) is asymptotically unbiased, with the bias decreasing as $O(m^{-2})$. For the rest of the paper we refer to (2.4) as a “bias-corrected” estimator, where the quotation marks highlight that the correction is not exact.

Payne and Mallick, (2015) propose a subsampling delayed acceptance MH following Christen and Fox, (2005). The aim in delayed acceptance is to simulate a Markov chain $\{\theta^{(j)}\}_{j=1}^{N}$ which admits the posterior

$$\pi(\theta)=\frac{p(y|\theta)p(\theta)}{p(y)},\text{ where }p(y)=\int p(y|\theta)p(\theta)d\theta\text{ and }p(\theta)\text{ denotes the prior,}$$

as invariant distribution. Moreover, the likelihood $p(y|\theta)$ should only be evaluated if there is a good chance of accepting the proposed $\theta$.

The algorithm in Payne and Mallick, (2015) proceeds as follows. Let $\theta_{c}=\theta^{(j)}$ denote the current state of the Markov chain. In the first stage, propose $\theta^{\prime}\sim q_{1}(\theta|\theta_{c})$ and compute

(3.1)

$\displaystyle\alpha_{1}(\theta_{c}\rightarrow\theta^{\prime})$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{\hat{p}_{m}(y|\theta^{\prime},u)p(\theta^{\prime})/q_{1}(\theta^{\prime}|\theta_{c})}{\hat{p}_{m}(y|\theta_{c},u)p(\theta_{c})/q_{1}(\theta_{c}|\theta^{\prime})}\right\},$

where $\hat{p}_{m}(y|\theta,u)$ is the estimator in Section 2.4 without control variates for $\hat{l}_{m}$ (but not “bias-corrected”, see below). Now, propose

$$\theta_{p}=\begin{cases}\theta^{\prime}&\quad\text{w.p.}\quad\alpha_{1}(\theta_{c},\theta^{\prime})\\ \theta_{c}&\quad\text{w.p.}\quad 1-\alpha_{1}(\theta_{c},\theta^{\prime}),\end{cases}$$

and move the chain to the next state $\theta^{(j+1)}=\theta_{p}$ with probability

(3.2)

$\displaystyle\alpha_{2}(\theta_{c}\rightarrow\theta_{p})$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{p(y|\theta_{p})p(\theta_{p})/q_{2}(\theta_{p}|\theta_{c})}{p(y|\theta_{c})p(\theta_{c})/q_{2}(\theta_{c}|\theta_{p})}\right\},$

where

$\displaystyle q_{2}(\theta_{p}|\theta_{c})$

$\displaystyle=$

$\displaystyle\alpha_{1}(\theta_{c}\rightarrow\theta_{p})q_{1}(\theta_{p}|\theta_{c})+r(\theta_{c})\delta_{\theta_{c}}(\theta_{p}),\quad r(\theta_{c})=1-\int\alpha_{1}(\theta_{c}\rightarrow\theta_{p})q_{1}(\theta_{p}|\theta_{c})d\theta_{p},$

and $\delta$ is the Dirac delta function. If rejected we set $\theta^{(j+1)}=\theta_{c}$.

Note that $\alpha_{2}$ in (3.2) is equivalent to the acceptance probability of a standard MH with a proposal density $q_{2}(\theta_{p}|\theta_{c})$ that is a mixture of two proposal densities: the first proposes to move from $\theta_{c}$ to $\theta_{p}$ (from the “slab” $q_{1}$) and the second proposes to stay at $\theta_{c}$ (from the “spike”). The mixture weight for the “slab” is $\alpha_{1}$ in (3.1) (with $\theta^{\prime}=\theta_{p}$): if the likelihood estimate is higher at $\theta^{\prime}$ compared to that of $\theta_{c}$ (after correcting with $q_{1}$) we propose $\theta^{\prime}=\theta_{p}$ with probability 1. Conversely, if it is lower, we propose to move but with a probability smaller than 1 which decreases the less likely we think that $\theta^{\prime}$ will be accepted as indicated by the estimated likelihood. The estimator $\hat{p}_{m}(y|\theta^{\prime},u)$ used by Payne and Mallick, (2015) does not depend on the current state of the Markov chain and hence the mixture weights in $q_{2}$ are state-independent. Convergence to the invariant distribution therefore follows from standard MH theory. The same applies when the estimator uses the control variates in Quiroz et al., (2016), or in Bardenet et al., (2015) but only for a fixed $\theta^{\star}$. Bardenet et al., (2015) suggest setting $\theta^{\star}=\theta_{c}$ every now and then to prevent that the control variates can be poor if the chain is far from $\theta^{\star}$: the resulting approximation is clearly state-dependent and standard MH theory does not apply. Instead, convergence to the invariant $\pi(\theta)$ is proved in Christen and Fox, (2005) and the delayed algorithm is exactly as above but with

$$\hat{p}_{m}(y|\cdot,u)=\hat{p}_{m}^{(\theta_{c})}(y|\cdot,u),\quad\text{emphasizing that it depends on the current state }\theta_{c}.$$

The state-dependent algorithm is a key ingredient when developing the delayed acceptance block PMMH in Section 4.

As noted above, Payne and Mallick, (2015) do not “bias correct” the likelihood estimate as they point out that the algorithm will have the correct invariant distribution anyway. We remark that without control variates it is not a good idea to apply the correction as the variance of $\hat{\sigma}_{m}^{2}$ is huge which adversely affects $\hat{p}_{m}(y|\cdot,u)$. An efficient estimate of $\hat{\sigma}_{m}^{2}$, however, would certainly improve $\hat{p}_{m}(y|\cdot,u)$. While the control variates allow us to estimate $\hat{\sigma}_{m}^{2}$ accurately, we will not implement the “bias-correction” when comparing to Payne and Mallick, (2015) in order to make the comparison fair.

It is beneficial to also update $u$ in order to avoid the risk of having a subset of observations for which the approximation is poor. This is especially important for the estimator in Payne and Mallick, (2015), as not using control variates can result in the particular subset having highly heterogeneous elements, which is detrimental for SRS. We note that updating $u$ is still a valid MH because (i) $u$ is not a state of the Markov chain and (ii) the distribution $p(u)=1/n^{m}$ (SRS) does not depend on $\theta$. Thus, the transition kernel of an algorithm that updates $u$ is a state-independent mixture of transition kernels, where each of the kernels satisfies detailed balance either by standard MH (or Christen and Fox, (2005) if the approximation depends on $\theta_{c}$). Since the weights $1/n^{m}$ do not depend on $\theta$, it follows that the mixture also satisfies detailed balance, and thus has $\pi(\theta)$ as invariant distribution. It is unnecessary to update $u$ in every iteration, instead one can update $u$ randomly to save the overhead cost of indexing the data matrix when obtaining the subset of observations.

We remark that delayed acceptance (sometimes with names as early rejection or surrogate transition), although without data subsampling, has been considered earlier in the literature. References include Fox and Nicholls, (1997); Liu, (2008); Cui et al., (2011); Smith, (2011); Solonen et al., (2012); Golightly et al., (2015); Sherlock et al., 2015a . Each stage in Banterle et al., (2014) (see Section 1) uses a partition of the data and can thus be considered as delayed acceptance with data subsampling. The advantage of our algorithm is that, because it only has two stages and the second stage evaluates the full data likelihood, we can instead estimate this likelihood in order to never do the evaluation for the full data set, see Section 4.

When considering efficiency of MCMC algorithms with additional computational costs (e.g. estimating the target), there are two types of fundamentally different efficiencies that interplay. The first is the statistical efficiency, which we will measure by the asymptotic (as the number of MCMC iterates go to infinity) variance of an estimate based on output from the Markov chain. Consider two MCMC algorithms $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ with the same invariant distribution. Then $\mathcal{A}_{1}$ is statistically more (less) efficient than $\mathcal{A}_{2}$ if it has a lower (higher) asymptotic variance. The second is the computational efficiency, which solely concerns “execution time” to produce a given number of iterates. The measures of statistical and computational efficiency (and a combination of them) used in this article are presented later in this section.

Sherlock et al., 2015b (in a non-subsampling context) study the statistical efficiency for delayed acceptance random walk Metropolis and, moreover, an efficiency that also takes into account the computational efficiency for the case where the target is estimated (DA-PMMH in Section 4).

Christen and Fox, (2005) note that, because the transition kernels of both the MH and delayed acceptance MH are derived from the same proposal $q_{1}$, and in addition $\alpha_{2}\leq 1$, the delayed acceptance MH will be less statistically efficient than MH. The intuition is that under these conditions the chain clearly exhibits a more “sticky” behavior and an estimate based on these samples will have a larger asymptotic variance under DA-MH than MH. Notice that the closer $\alpha_{2}$ is to $1$, the more statistically efficient the delayed acceptance algorithm is, and when $\alpha_{2}=1$ it is equivalent to the standard MH which gives the upper bound of the possible statistical efficiency achieved by a DA-MH.

Result 1 in Payne and Mallick, (2015) gives the alternative formulation (for state-independent approximations)

(3.3)

$\displaystyle\alpha_{2}(\theta_{c}\rightarrow\theta_{p})$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{\hat{p}_{m}(y|\theta_{c},u)/p(y|\theta_{c})}{\hat{p}_{m}(y|\theta_{p},u)/p(y|\theta_{p})}\right\}.$

Let $l_{k}(\theta_{c},\theta_{p})=l_{k}(\theta_{c})-l_{k}(\theta_{p})$ and denote by $\hat{l}_{m}(\theta_{c},\theta_{p})$ the estimate of $l(\theta_{c},\theta_{p})=\sum_{k=1}^{n}l_{k}(\theta_{c},\theta_{p})$. Similarly to (2.3),

(3.4)

$\displaystyle\hat{l}_{m}(\theta_{c},\theta_{p})$

$\displaystyle=$

$\displaystyle q(\theta_{c},\theta_{p})+\frac{1}{m}\sum_{i=1}^{m}\zeta_{i},\quad\text{with }q(\theta_{c},\theta_{p})=\sum_{k=1}^{n}q_{k}(\theta_{c},\theta_{p}),$

where $q_{k}(\theta_{c},\theta_{p})=q_{k}(\theta_{c})-q_{k}(\theta_{p})$ and the $\zeta_{i}$’s are iid with

$$\Pr\left(\zeta_{i}=nD_{k}(\theta_{c},\theta_{p})\right)=\frac{1}{n},\text{ with }D_{k}=\left(l_{k}(\theta_{c},\theta_{p})-q_{k}(\theta_{c},\theta_{p})\right)\quad\text{for }i=1,\dots m.$$

We can also show that

(3.5)

$$\mathrm{E}[\hat{l}_{m}(\theta_{c},\theta_{p})]=l(\theta_{c},\theta_{p})\text{ and }V[\hat{l}_{m}(\theta_{c},\theta_{p})]=\frac{\sigma_{\zeta}^{2}}{m}\text{ with }\sigma_{\zeta}^{2}=n\sum_{k=1}^{n}\left(D_{k}(\theta_{c},\theta_{p})-\bar{D}_{F}(\theta_{c},\theta_{p})\right)^{2},$$

where $\bar{D}_{F}$ is the mean over the full population. When not “bias-corrected”, the ratio appearing in (3.3) becomes

(3.6)

$$R_{m}=\exp\left(\hat{l}_{m}(\theta_{c},\theta_{p})-l(\theta_{c},\theta_{p})\right).$$

We now propose a theorem that relates $\mathrm{E}\left[\alpha_{2}(\theta_{c}\rightarrow\theta_{p})\right]$ to the variance $\sigma_{R}^{2}=\mathrm{V}[\hat{l}_{m}(\theta_{c},\theta_{p})]$ under the assumption that $\hat{l}_{m}(\theta_{c},\theta_{p})\sim\mathcal{N}(l(\theta_{c},\theta_{p}),\sigma_{R}^{2})$ (equivalently $R_{m}\sim\log\mathcal{N}(0,\sigma_{R}^{2})$ in (3.6)). In turn, $\alpha_{2}$ relates to the statistical efficiency as discussed above. The assumption of normality is justified by a standard CLT for $\hat{l}_{m}(\theta_{c},\theta_{p})$ since the $\zeta_{i}$’s are iid.

Theorem 1 says that if the variance of the log of the ratio $R_{m}$ in (3.6) is lower, then the algorithm has a higher $\alpha_{2}$ on average. Moreover, it shows that $\alpha_{2}$ deteriorates quickly with $\sigma_{R}$, which illustrates the importance of achieving a small $\sigma_{R}$.

Although the delayed acceptance MH is always less statistically efficient than MH it can of course still be more generally efficient in terms of balancing computational and statistical efficiency. Clearly, a necessary condition to achieve a more general efficient DA-MH algorithm than the corresponding MH requires that its computing time is faster, i.e. it must be computationally more efficient. When comparing DA-MH algorithms with the same computing time (by for example having the same subsample size) then Theorem 1 shows that a smaller variance of the log-ratio results in a more general efficient algorithm.

To also compare algorithms of different computing times we now define a measure for the general efficiency discussed above. In the rest of the paper, whenever we claim that algorithms are more or less efficient to each other it is based on this measure. The statistical (in)efficiency part is measured by the Inefficiency Factor (IF), which quantifies the amount by which the variance of $(1/N)\sum_{j=1}^{N}\theta^{(j)}$ is inflated when $\theta^{(j)}$ is obtained by Markov chain simulation compared to that of iid simulation. It is given by

(3.7)

$\displaystyle\mathrm{IF}$

$\displaystyle=$

$\displaystyle 1+2\sum_{l=1}^{\infty}\rho_{l},$

where $\rho_{l}$ is the auto-correlation at the $l$th lag of the chain, and can be computed with the coda package in R (Plummer et al.,, 2006). To include computational efficiency in the measure we use Effective Draws (ED) per computing unit

(3.8)

$\displaystyle\mathrm{ED}$

$\displaystyle=$

$\displaystyle\frac{N}{\mathrm{IF}\times t},$

where $N$ is the number of MCMC iterations and $t$ is the computing time. The measure of interest is the effective draws per computing time of a delayed acceptance algorithm $\mathcal{A}$ relative to that of MH, i.e.

(3.9)

$\displaystyle\mathrm{RED}$

$\displaystyle=$

$\displaystyle\frac{\mathrm{ED}^{\mathcal{A}}}{\mathrm{ED}^{\mathrm{MH}}}.$

Quiroz et al., (2016) propose a pseudo-marginal approach to subsampling based on the “bias-corrected” likelihood estimator in (2.4). In the next subsection this estimate replaces the true likelihood, thereby avoiding the evaluation of the full data set. To allow for a smaller subsample size without adversely affecting the mixing of the chain Quiroz et al., (2016) develop a correlated pseudo-marginal approach based on Deligiannidis et al., (2016). Tran et al., (2016) (see also Quiroz et al.,, 2016) use an alternative approach to correlated pseudo-marginal which we use for our Delayed Acceptance Block PMMH (DA-BPMMH) in Section 4.3. Although not considered here, it is straightforward to instead correlate the subsamples as in Deligiannidis et al., (2016).

Pseudo-marginal by data subsampling targets the following posterior on an augmented space $(\theta,u)$,

(4.1)

$$\tilde{\pi}_{m}(\theta,u)=\hat{p}_{m}(y|\theta,u)p(u)p(\theta)/p_{m}(y),\text{ with }p_{m}(y)=\int p_{m}(y|\theta)p(\theta)d\theta.$$

The algorithm is similar to MH except that $\theta$ and $u$ are proposed and accepted (or rejected) jointly, with probability

(4.2)

$\displaystyle\alpha_{\mathrm{PMMH}}$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{\hat{p}_{m}(y|\theta_{p},u_{p})p(\theta_{p})/q(\theta_{p}|\theta_{c})}{\hat{p}_{m}(y|\theta_{c},u_{c})p(\theta_{c})/q(\theta_{c}|\theta_{p})}\right\},$

where the subscripts denote the current and proposed values of $\theta$ and $u$. Quiroz et al., (2016) show that the algorithm converges to a slightly perturbed target (only approximately unbiased likelihood estimator) $\pi_{m}(\theta)$ and prove that

$$\frac{\left|\pi_{m}(\theta)-\pi(\theta)\right|}{\pi(\theta)}\leq O(m^{-2}).$$

Moreover, if $h(\theta)$ is a function that is absolutely integrable with respect to $\pi(\theta)$, then $\mathrm{E}_{\pi_{m}}\left[h(\theta)\right]$ is also within $O(m^{-2})$ of its true value.

The variance of $\hat{l}_{m}$ is crucial for the general efficiency of the algorithm: $\sigma^{2}$ in the approximate interval $[1,3.3]$ is optimal (Pitt et al.,, 2012; Doucet et al.,, 2015; Sherlock et al., 2015c, ). The correlated pseudo-marginal approach induces a high positive correlation between the estimates in (4.2) by correlating $u_{c}$ and $u_{p}$. This allows the use of a less precise estimator without getting stuck: the errors in the numerator and denominator tend to cancel. As a consequence, $\sigma^{2}$ can be larger than above, hence speeding up the algorithm by taking a smaller subsample. We follow Tran et al., (2016) and divide $u$ (defined in Section 2.2) into $G$ blocks,

$$u=(u_{(1)},\dots,u_{(G)})\quad\text{with }\frac{m}{G}\text{ observations within each block,}$$

and update a single block (randomly) in each iteration. Setting $G$ large induces a high positive correlation $\rho$ between $\hat{l}_{m}(\theta_{c})$ and $\hat{l}_{m}(\theta_{p})$: Tran et al., (2016) show that $\rho=1-1/G$ under certain assumptions. We set $G=100$ for the state-dependent algorithm in our application.

Delayed acceptance PMMH uses an estimate of the target in the second stage of the algorithm. Such algorithms have recently been considered by Sherlock et al., 2015b ; Sherlock et al., 2015a and Golightly et al., (2015). We propose a state-independent data subsampling DA-PMMH (uncorrelated PMMH) and a state-dependent (block PMMH, where the correlation between the subsamples follows Tran et al.,, 2016) extension DA-Block (correlated) PMMH (DA-BPMMH) in the next subsection.

The estimated “bias-corrected” likelihood in (2.4) is approximated in a first stage screening by $\hat{s}$ which is discussed in Section 4.4. Similarly to the algorithm in Section 3.1, we desire to only evaluate the (estimated) target (on the augmented space) if the proposed state is likely to be accepted. Let $\left(\theta_{c},u_{c}\right)=\left(\theta^{(j)},u^{(j)}\right)$ denote the current state of the augmented Markov chain. In the first stage, propose $\theta^{\prime}\sim q_{1}(\theta|\theta_{c})$ and $u^{\prime}\sim p(u)$ and evaluate

(4.3)

$\displaystyle\alpha_{1}^{\mathrm{PMMH}}\left\{\left(\theta_{c},u_{c}\right)\rightarrow\left(\theta^{\prime},u^{\prime}\right)\right\}$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{\hat{s}(\theta^{\prime},u^{\prime})p(\theta^{\prime})/q_{1}(\theta^{\prime}|\theta_{c})}{\hat{s}(\theta_{c},u_{c})p(\theta_{c})/q_{1}(\theta_{c}|\theta^{\prime})}\right\},$

where $\hat{s}(\theta,u)$ is the approximation of $\hat{p}_{m}(y|\theta,u)$ (2.4). Propose

$$\left(\theta_{p},u_{p}\right)=\begin{cases}\left(\theta^{\prime},u^{\prime}\right)&\quad\text{w.p.}\quad\alpha_{1}^{\mathrm{PMMH}}\left\{\left(\theta_{c},u_{c}\right)\rightarrow\left(\theta^{\prime},u^{\prime}\right)\right\}\\ \left(\theta_{c},u_{c}\right)&\quad\text{w.p.}\quad 1-\alpha_{1}^{\mathrm{PMMH}}\left\{\left(\theta_{c},u_{c}\right)\rightarrow\left(\theta^{\prime},u^{\prime}\right)\right\},\end{cases}$$

and move to $\left(\theta^{(j+1)},u^{(j+1)}\right)=\left(\theta_{p},u_{p}\right)$ with probability

(4.4)

$\displaystyle\alpha_{2}^{\mathrm{PMMH}}\left\{\left(\theta_{c},u_{c}\right)\rightarrow\left(\theta_{p},u_{p}\right)\right\}$

$\displaystyle=$

$\displaystyle\min\left\{1,\frac{\hat{p}_{m}(y|\theta_{p},u_{p})p(\theta_{p})/q_{2}(\theta_{p}|\theta_{c})}{\hat{p}_{m}(y|\theta_{c},u_{c})p(\theta_{c})/q_{2}(\theta_{c}|\theta_{p})}\right\},$

where

$\displaystyle q_{2}(\theta_{p}|\theta_{c})$

$\displaystyle=$

$\displaystyle\alpha_{1}^{\mathrm{PMMH}}q_{1}(\theta_{p}|\theta_{c})+r(\theta_{c})\delta_{\theta_{c}}(\theta_{p}),\quad r(\theta_{c})=1-\int\alpha_{1}^{\mathrm{PMMH}}q_{1}(\theta_{p}|\theta_{c})d\theta_{p},$

and $\delta$ is the Dirac delta function. If rejected we set $\left(\theta^{(j+1)},u^{(j+1)}\right)=\left(\theta_{c},u_{c}\right)$. Similarly to the argument for the DA-MH, we recognize $\alpha_{2}^{\mathrm{PMMH}}$ as the acceptance probability for a pseudo-marginal algorithm. Since the approximation is state-independent ($\hat{s}(\theta,u)$ independent of the current state implies state-independent mixture weights for $q_{2}$), convergence follows from being a pseudo-marginal algorithm (Andrieu and Roberts,, 2009). However, as the estimate is biased the target is perturbed (Quiroz et al.,, 2016) but within $O(m^{-2})$ as discussed in Section 4.1.

As $u$ is part of the state in a pseudo-marginal algorithm, obtaining a state-dependent approximation is easily achieved by correlating $u$: we sample $u_{p}\sim p(u|u_{c})$ thus $\hat{s}(\theta,u)=\hat{s}_{\theta_{c},u_{c}}(\theta,u)$. This is the state-dependent (on the augmented space) setup in Christen and Fox, (2005) and it follows that the invariant distribution is $\tilde{\pi}_{m}(\theta,u)$ in (4.1). This is the invariant distribution in the algorithm in Quiroz et al., (2016), which we already mentioned has a marginal $\pi_{m}(\theta)$ within $O(m^{-2})$ of $\pi(\theta)$. We note that this state-dependent approximation is very convenient because it automatically allows us to have a higher variance on $\hat{l}_{m}$ because of the block PMMH mechanism (see Section 4.1).

Our approximation is inspired by adaptive delayed acceptance ideas in Cui et al., (2011) and Sherlock et al., 2015a . However, our approach is not strictly adaptive as we only learn about the proposal for a fixed training period of $N_{\mathrm{train}}$ iterations, which are discarded from the final draws.

The idea is to use a sparser set of the data to construct control variates $q_{k}^{(1)}(\theta)$ in the first stage. During the training period we learn about the discrepancy between $q^{(1)}(\theta)=\sum_{k=1}^{n}q_{k}^{(1)}(\theta)$ and $q(\theta)$ of the second stage (obtained with a denser set), i.e.

(4.5)

$$\underbrace{q(\theta)-q^{(1)}(\theta)}_{=e(\theta)}=f(\theta)+\epsilon,\quad\epsilon\text{ is the noise (assumed independent of $\theta$)}.$$

Learning about $f$ is a standard regression problem: the collection of proposed $\theta$’s during the fixed training period are the inputs and the discrepancies are the training data.

The trivial decomposition

	$\displaystyle\hat{l}_{m}(\theta)$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{n}q_{k}(\theta)+\hat{d}_{m}(\theta)-\hat{\sigma}_{m}^{2}/2$
		$\displaystyle=$	$\displaystyle\sum_{k=1}^{n}q_{k}^{(1)}(\theta)+\left(\sum_{k=1}^{n}q_{k}(\theta)-\sum_{k=1}^{n}q_{k}^{(1)}(\theta)\right)+\hat{d}_{m}(\theta)-\hat{\sigma}_{m}^{2}/2$

suggests the first stage approximation

$$\hat{s}(\theta,u)=\sum_{k=1}^{n}q_{k}^{(1)}(\theta)+\hat{e}(\theta)+\hat{d}_{m}(\theta)-\hat{\sigma}_{m}^{2}/2,$$

where $\hat{e}(\theta)$ is the prediction of the discrepancy at $\theta$. Note that computing the true discrepancy requires the denser set to be evaluated, whereas the prediction is very fast. For example, in linear (or non-linear by basis functions) regression the parameters of $f(\theta)$ are estimated once after the training data has been collected. Prediction of $e(\theta)$ for a new “data-observation” $\theta$ is then a simple dot product computation, which is typically much faster than evaluating $q(\theta)$. Note that if the $u$’s are correlated then $\hat{s}(\theta,u)$ also depends on the current state, i.e. $\hat{s}_{\theta_{c},u_{c}}(\theta,u)$.

In Section 5 we implement both a linear regression and (noise free) Gaussian process to learn $f$ in (4.5), but any regression technique can be used.

We model the probability of bankruptcy conditional on a set of covariates using a data set of $534,717$ Swedish firms for the time period 1991-2008. We have in total $n=4,748,089$ firm-year observations. The variables included are: earnings before interest and taxes, total liabilities, cash and liquid assets, tangible assets, logarithm of deflated total sales and logarithm of firm age in years. We also include the macroeconomic variables GDP-growth rate (yearly) and the interest rate set by the Swedish central bank. See Giordani et al., (2014) for a detailed description of the data set.

We consider the logistic regression model

$\displaystyle p(y_{k}|x_{k},\beta)$

$\displaystyle=$

$\displaystyle\left(\frac{1}{1+\exp(x_{k}^{T}\beta)}\right)^{y_{k}}\left(\frac{1}{1+\exp(-x_{k}^{T}\beta)}\right)^{1-y_{k}},$

where $x_{k}$ includes the variables above plus an intercept term. We set $p(\beta)\sim N(0,10I)$ for simplicity. Since the bankruptcy observations $(y_{k}=1)$ are sparse in the data, we follow Payne and Mallick, (2015) and estimate the likelihood only for the $y_{k}=0$ observations. That is, we decompose

$\displaystyle l(\beta)$

$\displaystyle=$

$\displaystyle\sum_{\{k;y_{k}=1\}}l_{k}(\beta)+\sum_{\{k;y_{k}=0\}}l_{k}(\beta),$

and evaluate the first term whereas a random sample is only taken to estimate the second term. The second term clearly follows the structure presented in Section 2.1.

The quantity of interest is the effective draws as defined in (3.8). We now outline how to compute $t$ as CPU time and present an alternative measure independent of the implementation that is based number of evaluations. We first outline a robust measure of the CPU time.

The delayed acceptance algorithms we implement have two stages, where the first stage is filtering out draws unlikely to be accepted. One can view MH as an algorithm that does not filter any proposed draws at all: any draw is subject to an accept/reject decision based on the second stage likelihood. To make total CPU time comparisons fair between a Delayed Acceptance (DA) MH and its corresponding MH, we compute the total time for the latter (MH) by the median time the former (DA) spends in the second stage and multiply by the number of MCMC iterations $N$:

$$\mathrm{CPU}_{\mathrm{MH}}=N\times\text{median time $\mathrm{DA}$ stage 2}.$$

The median is used to avoid extreme values that can arise due to external disturbances (CPU time should be independent of $\theta$ here). For the delayed acceptance algorithm the CPU time is

$$\mathrm{CPU}_{\mathrm{DA}}=N\times\text{median time $\mathrm{DA}$ stage 1}\text{ }+\mathrm{FullEval\times\text{median time $\mathrm{DA}$ stage 2},}$$

where $\mathrm{FullEval}$ is the number of second stage evaluations the delayed acceptance performs.

The (total) number of evaluations measure is straightforward for MH ($N\times n$), DA-MH without control variates ($N\times m+\mathrm{FullEval}\times n$) and with ($N\times(K+m)+\mathrm{FullEval}\times n$), and PMMH/BPMMH ($N\times(K+m)$). For the delayed versions of PMMH/BPMMH we similarly add the different evaluations, but note that it will be different for the pre- and post- training period. Moreover, there is a one time cost of learning $f(\theta)$ and also a (post-training) cost of predicting $e(\theta)$. We translate these to number of evaluations as follows. First, for learning $f(\theta)$ we measure the CPU time it takes to fit it with the training data. We compare this time to the average CPU time for the iterations during the training period. We do so because we know the measure of the latter in terms of number of evaluations: $K^{(1)}+K+m$, where $K^{(1)}$ and $K$ are, respectively, the number of centroids in the sparser and denser set of data. Therefore, if learning $f(\theta)$ is say $T$ times slower in CPU time, then this translates to $T\cdot(K^{(1)}+K+m)$ number of evaluations. Finally, the number of evaluations for predicting a single $e(\theta)$ depends on the model for $f(\theta)$. For linear regression the prediction is a dot product which we assign the same cost as computing a single log-likelihood contribution (which is typically a function of a dot product). For a Gaussian process, the prediction requires evaluating the kernel for $N_{\mathrm{train}}$ observations and we let that define the number of evaluations for a single prediction.

We consider the following two examples. Example 1 estimates the model with DA-MH implemented with our efficient control variates and compares it to the implementation in Payne and Mallick, (2015). Example 2 estimates the model with DA-BPMMH and compare to the block PMMH algorithm in Quiroz et al., (2016). To make comparisons fair for Example 1 we use without replacement sampling (as in Payne and Mallick,, 2015). This sampling scheme is typically used together with the Horvitz-Thompson estimator (Horvitz and Thompson,, 1952): see Särndal et al., (2003) on how to modify the formulas in Section 2.3 for without replacement sampling.

Two main implementations of the difference estimator are considered. The first computes $q_{k}$ with the second order term evaluated at $\beta$, which we call dynamic. The second, which we call static, fixes the second order term at the optimum $\beta^{\star}$. The dynamic approach clearly provides a better approximation but is more expensive to compute. For both the dynamic and static approaches we compare four different sparse representations of the data for computing $q$ in (2.3), each with a different number of clusters. The clusters are obtained using Algorithm 1 in Quiroz et al., (2016) on the observations for which $y=0$ ($4,706,523$ observations). We note that, as more clusters are used to represent the data, the approximation of the likelihood is more accurate, although it is more expensive to compute.

We consider a Random walk MH proposal for $\beta$ where we learn the proposal scale during the first $N_{\mathrm{train}}=5,000$ (and also train $f(\theta)$) iterations in order to reach an acceptance probability of $\approx 0.23$ for MH (Roberts et al.,, 1997) and $\approx 0.10$ for BPMMH. For the delayed acceptance algorithms we have the same targets but for $\alpha_{1}$, i.e. the first stage acceptance probability. We discard the training samples and also a subsequent burn-in period of $10$% of the remaining samples ($20,000$) when doing inference. However, the computing costs (CPU and number of evaluations) include all $N$ iterations.

Finally, the delayed acceptance algorithms are implemented with an update of $u$ with probability $0.01$.

Tables 1 and 2 summarize the results, respectively, for the difference estimator with control variates (DE) and the estimator in Payne and Mallick, (2015) (PM). It is evident that the difference estimator has a larger second stage acceptance probability $\alpha_{2}$ (for a given sample size), which is a consequence of Theorem 1 because it has a lower $\sigma_{R}^{2}=V[\log(R_{m})].$ Figure 5.1 confirms that the normality assumptions, for the smallest value of $m$ and $K$ (the worst case scenario), are adequate for both methods. We also note from Table 2 that for some sample sizes Payne and Mallick, (2015) performs more poorly than the standard Metropolis-Hastings algorithm. One possible explanation is that the applications in Payne and Mallick, (2015) have a small number of continuous covariates (one in the first application and three in the second) and the rest are binary. It is clear that the continuous covariate case results in more variation among the log-likelihood contributions which is detrimental for SRS. In this application we have eight continuous covariates which explains why SRS without covariates performs poorly for small sampling fractions. As an example, for a subsample of $0.1\%$ of the data, not a single effective sample was obtained.

Our second example explores how the state-dependent delayed acceptance BPMMH improves the BPMMH proposed in Quiroz et al., (2016). The results are presented in Table 3. Quiroz et al., (2016) in turn show that they outperform other subsampling approaches and, in particular, that many of these approaches perform more poorly than the standard MH (see also Bardenet et al.,, 2015) in terms of efficiency and/or can give a very poor approximation of the posterior. Here we find that the delayed acceptance BPMMH is 30 times more efficient than MH, which is a huge improvement considering the aforementioned facts. Moreover, we find that the approximate posterior produced is very close to the true posterior (simulated by MH), as illustrated in Figure 5.2.