Subsampling MCMC - an Introduction for the Survey Statistician

Abstract

The rapid development of computing power and efficient Markov Chain Monte Carlo (MCMC) simulation algorithms have revolutionized Bayesian statistics, making it a highly practical inference method in applied work. However, MCMC algorithms tend to be computationally demanding, and are particularly slow for large datasets. Data subsampling has recently been suggested as a way to make MCMC methods scalable on massively large data, utilizing efficient sampling schemes and estimators from the survey sampling literature. These developments tend to be unknown by many survey statisticians who traditionally work with non-Bayesian methods, and rarely use MCMC. Our article explains the idea of data subsampling in MCMC by reviewing one strand of work, Subsampling MCMC, a so called Pseudo-Marginal MCMC approach to speeding up MCMC through data subsampling. The review is written for a survey statistician without previous knowledge of MCMC methods since our aim is to motivate survey sampling experts to contribute to the growing Subsampling MCMC literature.

Appendices

Appendix A: Algorithms

This appendix contains the main sampling algorithms discussed in the paper.

Appendix B: Details for the Poisson Regression Example

This appendix gives the details for the control variates in our illustrative Poisson regression example. Quiroz et al. (2018a) give general expressions for the gradients and hessians in the GLM class, and provide general compact expression that reduces the computational complexity of the control variates.

The Poisson Regression Model

The Poisson regression is of the form

$$y_{i} | \mathbf{x}_{i}, \boldsymbol{\theta} \overset{indep.}{\sim} \text{Pois}(\lambda_{i}), \hspace{0.5cm} \lambda_{i} = \exp(\alpha + \mathbf{x}_{i}^{T} \beta), $$

where 𝜃 = (α, β)^T.

Parameter-Expanded Control Variates

Let $\mathbf {w}_{i} = (1,\mathbf {x}_{i}^{T})^{T}$. The log-likelihood contribution from the i th observation is

$$\ell_{i}(\boldsymbol{\theta}) = y_{i} \mathbf{w}_{i}^{T} \boldsymbol{\theta} - \exp(\mathbf{w}_{i}^{T} \boldsymbol{\theta}) - \log(y_{i}!) $$

with gradient and hessian

$$\nabla_{\boldsymbol{\theta}} \ell_{i}(\boldsymbol{\theta}) = (y_{i} - \exp(\mathbf{w}_{i}^{T} \boldsymbol{\theta}))\mathbf{w}_{i} $$

$$\nabla_{\boldsymbol{\theta} \boldsymbol{\theta}^{T}}^{2} \ell_{i}(\boldsymbol{\theta}) = - \exp(\mathbf{w}_{i}^{T} \boldsymbol{\theta})\mathbf{w}_{i}\mathbf{w}_{i}^{T}. $$

Let μ(𝜃, x) = α + x^Tβ = w^T𝜃. The parameter-expanded control variate in (3.4) is then

$$\begin{array}{@{}rcl@{}} \ell_{i}(\boldsymbol{\theta}) &\approx& y_{i} \mu(\boldsymbol{\hat \theta},\mathbf{x}_{i}) - \exp(\mu(\boldsymbol{\hat \theta},\mathbf{x}_{i})) - \log(y_{i}!) \\ &&+ [y_{i} - \exp(\mu(\boldsymbol{\hat \theta},\mathbf{x}_{i}))](\mu_{i}(\theta)-\mu_{i}(\boldsymbol{\hat \theta})) \\ &&- \frac{1}{2}\exp(\mu(\boldsymbol{\hat \theta},\mathbf{x}_{i}))(\mu(\theta,\mathbf{x}_{i})-\mu(\boldsymbol{\hat \theta}, \mathbf{x}_{i}))^{2}. \end{array} $$

Data-Expanded Control Variates

The log-likelihood contribution from the i th observation is

$$\ell_{i}(\theta) = y_{i} (\alpha + \mathbf{x}_{i}^{T} \beta) - \exp(\alpha + \mathbf{x}_{i}^{T} \beta) - \log(y_{i}!) $$

with gradient and hessian

$$\nabla_{y_{i}} \ell_{i}(\theta) = \alpha + \mathbf{x}_{i}^{T} \beta - \psi_{0}(y_{i}+ 1), $$

where $\psi _{k}(z) = {\nabla _{z}^{k}} \log {\Gamma }(z)$ is the polygamma function of order k,

$$\begin{array}{*{20}l} \nabla_{\mathbf{x}_{i}} \ell_{i}(\theta) &= (y_{i} - \exp(\alpha + \mathbf{x}_{i}^{T} \beta)) \beta, \hspace{0.1cm} \nabla_{y_{i} y_{i}}^{2} \ell_{i}(\theta) = - \psi_{1}(y_{i}+ 1), \\ \nabla_{\mathbf{x}_{i} \mathbf{x}_{i}^{T}}^{2} \ell_{i}(\boldsymbol{\theta}) &= - \exp(\alpha + \mathbf{x}_{i}^{T} \beta) \beta \beta^{T}, \text{ and } \nabla_{y_{i} \mathbf{x}_{i}^{T}}^{2} \ell_{i}(\boldsymbol{\theta}) = \beta. \end{array} $$

We can write the gradients and hessian compactly by defining $\mathbf {z}_{i}=(y_{i},\mathbf {x}_{i}^{T})^{T}$,

$$\nabla_{z_{i}} \ell_{i}(\boldsymbol{\theta}) = \left[\begin{array}{l} \alpha + \mathbf{x}_{i}^{T} \beta - \psi_{0}(y_{i}+ 1) \\ (y_{i} - \exp(\alpha + \mathbf{x}_{i}^{T} \beta)) \beta \end{array}\right] $$

$$\nabla_{z_{i} {z_{i}^{T}}}^{2} \ell_{i}({\boldsymbol{\theta}}) = \left[\begin{array}{ll} - \psi_{1}(y_{i}+ 1) & \beta^{T} \\ \beta & - \exp(\alpha + \mathbf{x}_{i}^{T} \beta) \beta \beta^{T} \end{array}\right]. $$

Let μ(𝜃, x) = α + x^Tβ. The data-expanded control variate in Eq. 3.5 can after some simplifications be expressed as

$$\begin{array}{@{}rcl@{}} \ell_{i}(\boldsymbol{\theta}) &\approx& y_{c_{i}} \mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}) - \exp(\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}})) - \log(y_{c_{i}}!) \\ &&+ (y_{i} - y_{c_{i}})(\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}) - \psi_{0}(y_{c_{i}}+ 1)) -\frac{1}{2} (y_{i} - y_{c_{i}})^{2}\psi_{1}(y_{c_{i}}+ 1) \\ &&+[y_{i}-\exp(\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}))] (\mu(\boldsymbol{\theta}, \mathbf{x}_{i})-\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}) ) \\ &&-\frac{1}{2} \exp(\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}))(\mu(\boldsymbol{\theta}, \mathbf{x}_{i})-\mu(\boldsymbol{\theta}, \mathbf{x}_{c_{i}}) )^{2} . \end{array} $$