# The data augmentation algorithm

Vivekananda Roy, Kshitij Khare and James P. Hobert

###### Abstract

The data augmentation (DA) algorithms are popular Markov chain Monte Carlo
(MCMC) algorithms often used for sampling from intractable probability
distributions. This review article comprehensively surveys DA MCMC algorithms,
highlighting their theoretical foundations, methodological implementations,
and diverse applications in frequentist and Bayesian statistics. The article
discusses tools for studying the convergence properties of DA algorithms.
Furthermore, it contains various strategies for accelerating the speed of
convergence of the DA algorithms, different extensions of DA algorithms and
outlines promising directions for future research. This paper aims to serve as
a resource for researchers and practitioners seeking to leverage data
augmentation techniques in MCMC algorithms by providing key insights and
synthesizing recent developments.

##  1 Introduction

The data augmentation (DA) algorithm (Tanner and Wong,, 1987; Swendsen and
Wang,, 1987) is a widely used class of Markov chain Monte Carlo (MCMC)
algorithms. Suppose $f_{X}$ is the target probability density function (pdf)
on $\mathbb{R}^{p}$. Generally, the goal is to evaluate the mean of some
function $h:\mathbb{R}^{p}\rightarrow\mathbb{R}$ with respect to $f_{X}$, that
is, to find $\mbox{E}_{f_{X}}h=\int_{\mathbb{R}^{p}}h(x)f_{X}(x)\,dx$. Since
$f_{X}$ is usually available only up to a normalizing constant,
$\mbox{E}_{f_{X}}h$ cannot be computed analytically. Furthermore, these days,
the target densities arising in machine learning, physics, statistics, and
other areas are so complex that direct simulation from $f_{X}$ is impossible,
making the classical Monte Carlo approximation of $\mbox{E}_{f_{X}}h$ based on
independent and identically distributed (iid) samples from $f_{X}$ infeasible.
In such situations, it is often possible, as described later in a section, to
construct a DA Markov chain $\\{X_{n}\\}_{n\geq 1}$, which has $f_{X}$ as its
stationary density. If this DA chain $\\{X_{n}\\}_{n\geq 1}$ is well behaved
then $\mbox{E}_{f_{X}}h$ can be consistently estimated by the sample average
$\bar{h}_{n}:=\sum_{i=1}^{n}h(X_{i})/n$.

For constructing a DA algorithm with stationary density $f_{X}$, one needs to
find a joint density $f(x,y)$ on $\mathbb{R}^{p}\times\mathbb{R}^{q}$ with
‘augmented variables’ $y$ satisfying the following two properties:

  1. (i)

the $x-$marginal of the joint density $f(x,y)$ is the target density $f_{X}$,
that is, $\int_{\mathbb{R}^{q}}f(x,y)\,dy=f_{X}(x)$, and

  2. (ii)

sampling from the two corresponding conditional pdfs, $f_{X|Y}(x|y)$ and
$f_{Y|X}(y|x)$ is straightforward.

The first property makes sure that the DA algorithm presented below has
$f_{X}$ as its stationary density, while the second property allows for
straightforward simulation of this DA Markov chain. As mentioned in Hobert,
(2011), the popularity of the DA algorithm is due in part to the fact that,
given an intractable pdf $f_{X}$, there are some general techniques available
for constructing a potentially useful joint density $f(x,y)$. For example, the
augmented variables $y$ generally correspond to the missing data used in the
EM algorithm (Dempster et al.,, 1977) for defining the complete data density,
which in turn is used for finding the maximum likelihood estimates of some
parameters. In Sections 2–5 we will provide several examples of widely used
nonlinear and/or high-dimensional statistical models with complex target
densities where using appropriate augmented data, popular and efficient DA
algorithms are constructed.

Each iteration of the DA algorithm consists of two steps — a draw from
$f_{Y|X}$ followed by a draw from $f_{X|Y}$. Indeed, if the current state of
the DA Markov chain is $X_{n}=x$, we simulate $X_{n+1}$ as follows.

Algorithm The $n$th iteration for the DA algorithm

1: Given $x$, draw $Y\sim f_{Y|X}(\cdot|x)$, and call the observed value $y$.

2: Draw $X_{n+1}\sim f_{X|Y}(\cdot|y)$.

From the two steps in each iteration of the DA algorithm, it follows that the
Markov transition density (Mtd) of the DA Markov chain $\\{X_{n}\\}_{n\geq 1}$
is given by

$k(x^{\prime}|x)=\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y)f_{Y|X}(y|x)\,dy\;.$
(1.1)

Let $f_{Y}(y)=\int_{\mathbb{R}^{p}}f(x,y)\,dx$ be the marginal density of the
augmented variable $y$. Since

$k(x^{\prime}|x)f_{X}(x)=f_{X}(x)\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y)f_{Y|X}(y|x)\,dy=\int_{\mathbb{R}^{q}}\frac{f(x^{\prime},y)f(x,y)}{f_{Y}(y)}\,dy=k(x|x^{\prime})f_{X}(x^{\prime}),$
(1.2)

for all $x,x^{\prime}$, that is, the Mtd $k$ satisfies the detailed balance
condition, we have

$\int_{\mathbb{R}^{p}}k(x^{\prime}|x)f_{X}(x)\,dx=f_{X}(x^{\prime})\;.$ (1.3)

Thus, $f_{X}$ is stationary for the DA Markov chain $\\{X_{n}\\}_{n\geq 1}$.
Moreover, if $\\{X_{n}\\}_{n\geq 1}$ is appropriately irreducibile, then
$\bar{h}_{n}$ is a strongly consistent estimator of $\mbox{E}_{f_{X}}h$, that
is, $\bar{h}_{n}$ converges to $\mbox{E}_{f_{X}}h$ almost surely as
$n\rightarrow\infty$ (Asmussen and Glynn,, 2011, Theorem 2). Since the DA
chain has the Mtd (1.1), irreducibility also implies Harris recurrence of the
DA Markov chain. Furthermore, if $\\{X_{n}\\}_{n\geq 1}$ is also aperiodic
then it is Harris ergodic and the marginal density of $X_{n}$ converges to the
stationary density $f_{X}$, no matter what the initial distribution of $X_{1}$
is (Meyn and Tweedie,, 1993, chapter 13). A simple sufficient condition which
guarantees both irreducibility and aperiodicity of $\\{X_{n}\\}_{n\geq 1}$ is
that the Mtd $k$ is strictly positive everywhere. Note that, if $f(x,y)$ is
strictly positive on $\mathbb{R}^{p}\times\mathbb{R}^{q}$, then $k$ is
strictly positive everywhere and the DA Markov chain $\\{X_{n}\\}_{n\geq 1}$
is Harris ergodic.

In order to keep the notations and discussions simple, here we consider the
target and augmented state space as Euclidean. Similarly, the densities
$f_{X}$ as well as $f(x,y)$ are considered with respect to the Lebesgue
measure. However, all methods and theoretical results discussed in this
chapter extend to discrete as well as other general settings. Indeed, there
are popular DA algorithms where the target and/or the augmented state space is
not continuous, for example consider the widely used DA algorithm for the
mixture models (Diebolt and Robert,, 1994).

##  2 Examples

In this section, we provide examples of some widely used high dimensional
linear models, generalized linear models, and generalized linear mixed models
where the target distributions are intractable. Then we describe some DA
algorithms that are used to fit these models.

###  2.1 Data augmentation for Bayesian variable selection

In modern datasets arising from diverse applications such as genetics, medical
science, and other scientific disciplines, the number of covariates are often
larger than the number of observed samples. The ordinary least squares method
for estimating the regression coefficients in a linear regression is not
applicable in such situations. On the other hand, penalized regression methods
such as the least absolute shrinkage and selection operator (lasso)
(Tibshirani,, 1996) are applicable as they can simultaneously induce shrinkage
and sparsity in the estimation of the regression coefficients.

Consider the standard linear model

$Z|\mu,\beta,\sigma^{2}\sim N_{m}(\mu 1_{m}+W\beta,\sigma^{2}I_{m}),$

where $Z\in\mathbb{R}^{m}$ is the vector of responses, $\mu\in\mathbb{R}$ is
the intercept, $1_{m}$ is the $m\times 1$ vector of 1’s, $W$ is the $m\times
p$ (standardized) covariate matrix,
$\beta=(\beta_{1},\dots,\beta_{p})\in\mathbb{R}^{p}$ is the vector of
regression coefficients, and $\sigma^{2}$ is the variance parameter. The lasso
is based on $L_{1}$ norm regularization, and it estimates $\beta$ by solving

$\underset{\beta}{\mbox{min}}\;(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+\lambda\sum_{j=1}^{p}|\beta_{j}|,$
(2.1)

for some shrinkage parameter $\lambda\in\mathbb{R}$, where $z$ is the vector
of observed responses, and $\tilde{z}=z-{\bar{z}}1_{m}$. The lasso estimate
can be interpreted as the posterior mode when, conditional on $\sigma^{2}$,
the regression parameters $\beta_{j}$’s have independent and identical Laplace
priors (Tibshirani,, 1996). Indeed, following Park and Casella, (2008), we
consider the hierarchical model

$\displaystyle Z|\mu,\beta,\sigma^{2}$ $\displaystyle\sim$ $\displaystyle
N_{m}(\mu 1_{m}+W\beta,\sigma^{2}I_{m}),$ $\displaystyle\pi(\mu)\propto
1;\beta|\sigma^{2}$ $\displaystyle\sim$
$\displaystyle\prod_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^{2}}}e^{-\lambda|\beta_{j}|/\sqrt{\sigma^{2}}}$
$\displaystyle\sigma^{2}$ $\displaystyle\sim$ $\displaystyle\pi(\sigma^{2}),$
(2.2)

where $\pi(\sigma^{2})$ is the prior density of $\sigma^{2}$. Since the
columns of $W$ are centered, standard calculations show that

$f(z|\beta,\sigma^{2})\equiv\int_{\mathbb{R}}f(z|\mu,\beta,\sigma^{2})\pi(\mu)d\mu=\frac{1}{(2\pi)^{(n-1)/2}\sigma^{n-1}}\exp\Big{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)}{2\sigma^{2}}\Big{]}.$
(2.3)

Thus, from (2.1) and (2.3), the joint posterior density of
$(\beta,\sigma^{2})$ is

$\pi(\beta,\sigma^{2}|z)\propto\frac{1}{(\sigma^{2})^{(n-1+p)/2}}\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+2\lambda\sigma\sum_{j=1}^{p}|\beta_{j}|}{2\sigma^{2}}\Bigg{]}\pi(\sigma^{2}).$
(2.4)

From (2.1) and (2.4) it follows that the mode of the (conditional on
$\sigma^{2}$) posterior density of $\beta$ is the lasso estimate.

The posterior density (2.4) is intractable. Introducing augmented variables
$y=(y_{1},\dots,y_{p})$ with $y_{i}>0$ for all $i$, Park and Casella, (2008)
consider the following joint density of $(\beta,\sigma^{2},y)$

$\displaystyle\pi(\beta,\sigma^{2},y|z)$
$\displaystyle\propto\frac{1}{(\sigma^{2})^{(n-1+p)/2}}\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)}{2\sigma^{2}}\Bigg{]}$
$\displaystyle\hskip
72.26999pt\times\Bigg{[}\prod_{j=1}^{p}\frac{1}{\sqrt{y_{j}}}\exp\bigg{\\{}-\frac{\beta_{j}^{2}}{2\sigma^{2}y_{j}}\bigg{\\}}\exp\bigg{\\{}-\frac{\lambda^{2}y_{j}}{2}\bigg{\\}}\Bigg{]}\pi(\sigma^{2}).$
(2.5)

Replacing $t$ with $\beta_{j}/\sigma$, $s$ with $y_{j}$, and $a$ with
$\lambda$ in the following representation of the Laplace density as a scale
mixture of normals (Andrews and Mallows,, 1974)

$\frac{a}{2}\exp(-a|t|)=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi
s}}\exp\Big{(}-\frac{t^{2}}{2s}\Big{)}\frac{a^{2}}{2}\exp\Big{(}-\frac{a^{2}}{2}s\Big{)}ds,\;a>0$
(2.6)

we see that

$\int_{\mathbb{R}^{p}}\pi(\beta,\sigma^{2},y|z)dy=\pi(\beta,\sigma^{2}|z),$

that is, the $(\beta,\sigma^{2})-$ marginal density of the joint posterior
density $\pi(\beta,\sigma^{2},y|z)$ given in (2.1) is the target posterior
density (2.4). Thus, from Section 1, if sampling from the two conditional
densities $\pi(\beta,\sigma^{2}|y,z)$ and $\pi(y|\beta,\sigma^{2},z)$ is
straightforward, then we can construct a valid DA algorithm for (2.4).

From (2.1), we have

$\pi(y|\beta,\sigma^{2},z)\propto\prod_{j=1}^{p}\frac{1}{\sqrt{y_{j}}}\exp\bigg{\\{}-\frac{\beta_{j}^{2}}{2\sigma^{2}y_{j}}\bigg{\\}}\exp\Big{\\{}-\frac{\lambda^{2}y_{j}}{2}\Big{\\}}.$
(2.7)

Recall that the Inverse-Gaussian $(\kappa,\psi)$ density is given by

$f(u)=\sqrt{\frac{\psi}{2\pi}}u^{-\frac{3}{2}}\exp\Big{\\{}-\frac{\psi(u-\kappa)^{2}}{2\kappa^{2}u}\Big{\\}}.$

Thus from (2.7), it follows that

$\frac{1}{y_{j}}\Big{|}\beta,\sigma^{2},z\stackrel{{\scriptstyle
ind}}{{\sim}}\mbox{Inverse-
Gaussian}\Big{(}\sqrt{\frac{\lambda^{2}\sigma^{2}}{\beta^{2}_{j}}},\lambda^{2}\Big{)}\;\mbox{for}\;j=1,2,\dots,p.$
(2.8)

In order to derive the conditional density $\pi(\beta,\sigma^{2}|y,z)$ we
assume that apriori $\sigma^{2}\sim\mbox{Inverse-Gamma}(\alpha,\xi)$ for some
$\alpha\geq 0$ and $\xi\geq 0$. The improper prior
$\pi(\sigma^{2})=1/\sigma^{2}$ used in Park and Casella, (2008) is obtained by
replacing $\alpha=0,\xi=0$. Note that a draw from $\pi(\beta,\sigma^{2}|y,z)$
can be made by first drawing from $\pi(\sigma^{2}|y,z)$ followed by a draw
from $\pi(\beta|\sigma^{2},y,z)$ (Rajaratnam et al.,, 2019). Denoting the
$p\times p$ diagonal matrix with diagonal elements $(y_{1},\dots,y_{p})$ by
$D_{y}$, we have

$\displaystyle\pi(\beta|\sigma^{2},y,z)$
$\displaystyle\propto\pi(\beta,\sigma^{2},y|z)\propto\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+\beta^{\top}D_{y}^{-1}\beta}{2\sigma^{2}}\Bigg{]}$
$\displaystyle\propto\exp\Bigg{[}-\frac{\beta^{\top}(W^{\top}W+D_{y}^{-1})\beta-2\beta^{\top}W^{\top}\tilde{z}}{2\sigma^{2}}\Bigg{]}.$

Thus, $\beta|\sigma^{2},y,z\sim
N_{p}((W^{\top}W+D_{y}^{-1})^{-1}W^{\top}\tilde{z},\sigma^{2}(W^{\top}W+D_{y}^{-1})^{-1})$.
Also, from (2.1), we have

$\displaystyle\pi(\sigma^{2}|y,z)$
$\displaystyle\propto\pi(\sigma^{2},y|z)=\int_{\mathbb{R}^{p}}\pi(\beta,\sigma^{2},y|z)d\beta$
$\displaystyle\propto\int_{\mathbb{R}^{p}}\frac{1}{(\sigma^{2})^{(n-1+p)/2}}\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+\beta^{\top}D_{y}^{-1}\beta}{2\sigma^{2}}\Bigg{]}\pi(\sigma^{2})d\beta$
$\displaystyle\propto\frac{\pi(\sigma^{2})}{(\sigma^{2})^{(n-1)/2}}\exp\Big{(}-\frac{\tilde{z}^{\top}\tilde{z}}{2\sigma^{2}}\Big{)}\int_{\mathbb{R}^{p}}(\sigma^{2})^{-p/2}\exp\Bigg{[}-\frac{\beta^{\top}(W^{\top}W+D_{y}^{-1})\beta-2\beta^{\top}W^{\top}\tilde{z}}{2\sigma^{2}}\Bigg{]}d\beta$
$\displaystyle\propto(\sigma^{2})^{-(n-1+2\alpha)/2-1}\exp\Bigg{[}-\frac{\tilde{z}^{\top}(I-W(W^{\top}W+D_{y}^{-1})^{-1}W^{\top})\tilde{z}+2\xi}{2\sigma^{2}}\Bigg{]}.$

Thus, $\sigma^{2}|y,z\sim\mbox{Inverse-
Gamma}((n-1+2\alpha)/2,(\tilde{z}^{\top}(I-W(W^{\top}W+D_{y}^{-1})^{-1}W^{\top})\tilde{z}+2\xi)/2).$
Hence, a single iteration of the DA algorithm uses the following two steps to
move from $(\beta,\sigma^{2})$ to $(\beta^{\prime},{\sigma^{2}}^{\prime})$.

Algorithm One iteration of the DA algorithm for Bayesian lasso

1: Given $(\beta,\sigma^{2})$, draw $y_{1},\dots,y_{p}$ independently with
$\frac{1}{y_{j}}\Big{|}\beta,\sigma^{2},z\sim\mbox{Inverse-
Gaussian}\bigg{(}\sqrt{\frac{\lambda^{2}\sigma^{2}}{\beta^{2}_{j}}},\lambda^{2}\bigg{)}\;\mbox{for}\;j=1,2,\dots,p.$

2: Draw $(\beta^{\prime},{\sigma^{2}}^{\prime})$ by first drawing
${\sigma^{2}}^{\prime}|y,z\sim\mbox{Inverse-
Gamma}\Big{(}\frac{n-1}{2}+\alpha,\frac{\tilde{z}^{\top}(I-W(W^{\top}W+D_{y}^{-1})^{-1}W^{\top})\tilde{z}}{2}+\xi\Big{)},$
and then drawing  ${\beta}^{\prime}|{\sigma^{2}}^{\prime},y,z\sim
N_{p}((W^{\top}W+D_{y}^{-1})^{-1}W^{\top}\tilde{z},{\sigma^{2}}^{\prime}(W^{\top}W+D_{y}^{-1})^{-1}).$

Although the lasso estimator has been extensively used in applications as
diverse as agriculture, genetics, and finance, it may perform unsatisfactorily
if the predictors are highly correlated. For example, if there is a group
structure among the variables, lasso tends to select only one variable from
each group. Zou and Hastie, (2005) proposed the Elastic Net (EN) to achieve
better performance in such situations. The EN estimator is obtained by solving

$\underset{\beta}{\mbox{min}}\;(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+\lambda_{1}\sum_{j=1}^{p}|\beta_{j}|+\lambda_{2}\sum_{j=1}^{p}|\beta_{j}|^{2},$
(2.9)

where $\lambda_{1}$ and $\lambda_{2}$ are tuning parameters. From (2.9) we see
that the elastic net uses both an $L_{1}$ penalty as in lasso and an $L_{2}$
penalty as in the ridge regression. Following Kyung et al., (2010) and Roy and
Chakraborty, (2017), we consider the hierarchical Bayesian EN model

$\displaystyle Z|\mu,\beta,\sigma^{2}$ $\displaystyle\sim$ $\displaystyle
N_{m}(\mu 1_{m}+W\beta,\sigma^{2}I_{m}),$ $\displaystyle\pi(\mu)\propto
1;\pi(\beta|\sigma^{2})$ $\displaystyle\propto$
$\displaystyle\prod_{j=1}^{p}\frac{1}{\sigma}\exp\Big{\\{}-\frac{\lambda_{1}|\beta_{j}|}{\sigma}-\frac{\lambda_{2}\beta^{2}_{j}}{2\sigma^{2}}\Big{\\}},$
$\displaystyle\sigma^{2}$ $\displaystyle\sim$ $\displaystyle\pi(\sigma^{2}),$
(2.10)

where $\pi(\sigma^{2})$ is the prior density of $\sigma^{2}$. From (2.1) and
(2.3) it follows that the joint posterior density of $(\beta,\sigma^{2})$ is

$\pi_{\text{EN}}(\beta,\sigma^{2}|z)\propto\frac{\pi(\sigma^{2})}{(\sigma^{2})^{(n-1+p)/2}}\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+2\lambda_{1}\sigma\sum_{j=1}^{p}|\beta_{j}|+\lambda_{2}\sum_{j=1}^{p}\beta_{j}^{2}}{2\sigma^{2}}\Bigg{]}.$
(2.11)

From (2.9) and (2.11) we see that the mode of the (conditional on
$\sigma^{2}$) posterior density of $\beta$ is the EN estimate. Since the
density (2.11) is intractable, using augmented variables
$y=(y_{1},\dots,y_{p})$, following Kyung et al., (2010) and Roy and
Chakraborty, (2017), we consider the joint density of $(\beta,\sigma^{2},y)$
given by

$\displaystyle\pi_{\text{EN}}(\beta,\sigma^{2},y|z)$
$\displaystyle\propto\frac{\pi(\sigma^{2})}{(\sigma^{2})^{(n-1+p)/2}}\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)}{2\sigma^{2}}\Bigg{]}$
$\displaystyle\hskip
36.135pt\times\Bigg{[}\prod_{j=1}^{p}\frac{1}{\sqrt{y_{j}}}\exp\bigg{\\{}-\frac{(\lambda_{2}+1/y_{j})\beta_{j}^{2}}{2\sigma^{2}}\bigg{\\}}\exp\bigg{\\{}-\frac{\lambda_{1}^{2}y_{j}}{2}\bigg{\\}}\Bigg{]}.$
(2.12)

From (2.6) we see that

$\int_{\mathbb{R}^{p}}\pi_{\text{EN}}(\beta,\sigma^{2},y|z)dy=\pi_{\text{EN}}(\beta,\sigma^{2}|z).$

From (2.1) using similar calculations as before, we see that the conditional
distributions of $1/y_{j}$’s are the same as (2.8) with $\lambda=\lambda_{1}$.
As before, we assume that apriori $\sigma^{2}\sim\mbox{Inverse-
Gamma}(\alpha,\xi)$ for some $\alpha\geq 0$ and $\xi\geq 0$. Denoting the
$p\times p$ diagonal matrix with diagonal elements
$(1/(\lambda_{2}+1/y_{1}),\dots,1/(\lambda_{2}+1/y_{p}))$ by $\tilde{D}_{y}$,
we have

$\displaystyle\pi_{\text{EN}}(\beta|\sigma^{2},y,z)\propto\pi_{\text{EN}}(\beta,\sigma^{2},y|z)\propto\exp\Bigg{[}-\frac{(\tilde{z}-W\beta)^{\top}(\tilde{z}-W\beta)+\beta^{\top}\tilde{D}_{y}^{-1}\beta}{2\sigma^{2}}\Bigg{]}.$

Thus, $\beta|\sigma^{2},y,z\sim
N_{p}((W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top}\tilde{z},\sigma^{2}(W^{\top}W+\tilde{D}_{y}^{-1})^{-1})$.
Next, doing similar calculations as for the Bayesian lasso, we see that
$\sigma^{2}|y,z\sim\mbox{Inverse-
Gamma}((n-1+2\alpha)/2,(\tilde{z}^{\top}(I-W(W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top})\tilde{z}+2\xi)/2).$
Hence, a single iteration of the DA algorithm for the Bayesian EN uses the
following two steps to move from $(\beta,\sigma^{2})$ to
$(\beta^{\prime},{\sigma^{2}}^{\prime})$.

Algorithm One iteration of the DA algorithm for Bayesian elastic net

1: Given $(\beta,\sigma^{2})$, draw $y_{1},\dots,y_{p}$ independently with
$\frac{1}{y_{j}}\Big{|}\beta,\sigma^{2},z\sim\mbox{Inverse-
Gaussian}\bigg{(}\sqrt{\frac{\lambda_{1}^{2}\sigma^{2}}{\beta^{2}_{j}}},\lambda_{1}^{2}\bigg{)}\;\mbox{for}\;j=1,2,\dots,p.$

2: Draw $(\beta^{\prime},{\sigma^{2}}^{\prime})$ by first drawing
${\sigma^{2}}^{\prime}|y,z\sim\mbox{Inverse-
Gamma}\Big{(}\frac{n-1}{2}+\alpha,\frac{\tilde{z}^{\top}(I-W(W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top})\tilde{z}}{2}+\xi\Big{)},$
and then drawing  ${\beta}^{\prime}|{\sigma^{2}}^{\prime},y,z\sim
N_{p}((W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top}\tilde{z},{\sigma^{2}}^{\prime}(W^{\top}W+\tilde{D}_{y}^{-1})^{-1}).$

Different Bayesian generalized lasso methods, such as the Bayesian group
lasso, the Bayesian sparse group lasso, and the Bayesian fused lasso models
have been proposed in the literature to handle the situations where the
covariates are known to have some structures, for example, when they form
groups or are ordered in some way (Kyung et al.,, 2010; Xu and Ghosh,, 2015).
Introducing appropriate augmented variables $y$, DA algorithms for these
models can also be constructed (Kyung et al.,, 2010; Jin and Tan,, 2021).

###  2.2 Data augmentation for Bayesian logistic models

Logistic regression model is likely the most popular model for analyzing
binomial data. Let $(Z_{1},Z_{2},\dots,Z_{m})$ be a vector of independent
Binomial random variables, and $w_{i}$ be the $p\times 1$ vector of known
covariates associated with $Z_{i}$ for $i=1,\dots,m$. Let
$\beta\in\mathbb{R}^{p}$ be the vector of unknown regression coefficients.
Assume that $Z_{i}\sim$ Binomial $(\ell_{i},F(w_{i}^{\top}\beta))$ where
$F(t)=e^{t}/(1+e^{t})$ is the cumulative distribution function of the standard
logistic random variable. Denoting the observed responses by
$z=(z_{1},z_{2},\dots,z_{m})$, the likelihood function for the logistic
regression model is

$L(\beta|z)=\prod_{i=1}^{m}{\ell_{i}\choose
z_{i}}\frac{\left[\exp\left(w_{i}^{\top}\beta\right)\right]^{z_{i}}}{[1+\exp\left(w_{i}^{\top}\beta\right)]^{\ell_{i}}}.$

Consider a Bayesian analysis that employs the following Gaussian prior for
$\beta$

$\displaystyle\pi(\beta)\propto\exp\Big{[}-\frac{1}{2}(\beta-\mu_{0})^{\top}Q(\beta-\mu_{0})\Big{]},$
(2.13)

where $\mu_{0}\in\mathbb{R}^{p}$ and $Q$ is a $p\times p$ positive definite
matrix. Then the intractable posterior density of $\beta$ is given by

$\pi(\beta\,|\,z)\propto\frac{L(\beta|z)\pi(\beta)}{c(z)}=\frac{1}{c(z)}\prod_{i=1}^{m}{\ell_{i}\choose
w_{i}}\frac{\left[\exp\left(w_{i}^{\top}\beta\right)\right]^{z_{i}}}{[1+\exp\left(w_{i}^{\top}\beta\right)]^{\ell_{i}}}\pi(\beta)\;,$
(2.14)

where

$c(z)=\int_{\mathbb{R}^{p}}\prod_{i=1}^{m}{\ell_{i}\choose
w_{i}}\frac{\left[\exp\left(w_{i}^{\top}\beta\right)\right]^{z_{i}}}{[1+\exp\left(w_{i}^{\top}\beta\right)]^{\ell_{i}}}\pi(\beta)d\beta$
(2.15)

is the marginal density of $z$.

There have been many attempts (Holmes and Held,, 2006; Frühwirth-Schnatter and
Frühwirth,, 2010) to produce an efficient DA algorithm for the logistic
regression model, that is (2.14), without much success until recently, when
Polson et al., (2013) (denoted as PS&W hereafter) proposed a new DA algorithm
based on the Pólya-Gamma (PG) distribution. A random variable $\theta$ follows
a PG distribution with parameters $a>0,\,b\geq 0$, if

$\theta\overset{d}{=}(1/(2\pi^{2}))\sum_{i=1}^{\infty}g_{i}/[(i-1/2)^{2}+b^{2}/(4\pi^{2})],$

where $g_{i}\overset{iid}{\sim}$Gamma$(a,1)$. From Wang and Roy, 2018c , the
pdf for PG$(a,b),a>0,\,b\geq 0$ is

$p(\theta\mid
a,b)=\Big{[}\cosh\big{(}\frac{b}{2}\big{)}\Big{]}^{a}\frac{2^{a-1}}{\Gamma{(a)}}\sum_{r=0}^{\infty}(-1)^{r}\frac{\Gamma{(r+a)}}{\Gamma{(r+1)}}\frac{(2r+a)}{\sqrt{2\pi\theta^{3}}}\exp\Big{(}-\frac{(2r+a)^{2}}{8\theta}-\frac{\theta
b^{2}}{2}\Big{)},$ (2.16)

for $\theta>0,$ where the hyperbolic cosine function
$\cosh(t)=(e^{t}+e^{-t})/2$. We denote this as $\theta\sim$PG$(a,b)$.

Let $y=(y_{1},y_{2},...y_{n})$, $k_{i}=z_{i}-\ell_{i}/2,\,i=1,...,n$ and
$p(y_{i}\mid\ell_{i},0)$ be the pdf of $y_{i}\sim$PG$(\ell_{i},0)$. Define the
joint posterior density of $\beta$ and $y$ given $z$ as

$\pi(\beta,y\mid
z)=\frac{1}{c(z)}\bigg{[}\prod_{i=1}^{m}\frac{\exp\\{k_{i}w_{i}^{\top}\beta-
y_{i}(w_{i}^{\top}\beta)^{2}/2\\}}{2}p(y_{i}|\ell_{i},0)\bigg{]}\pi(\beta).$
(2.17)

By Theorem 1 in Polson et al., (2013), it follows that the $\beta$-marginal
density of the joint posterior density $\pi(\beta,y\mid z)$ is the target
density $\pi(\beta\,|\,z)$ given in (2.14), that is,

$\displaystyle\pi(\beta\mid
z)=\int_{\mathbb{R}_{+}^{m}}\bigg{[}\prod_{i=1}^{m}\frac{\exp\\{k_{i}w_{i}^{\top}\beta-
y_{i}(w_{i}^{\top}\beta)^{2}/2\\}}{2}p(y_{i}\mid\ell_{i},0)\bigg{]}dy\times\frac{\pi(\beta)}{c(z)}.$

PS&W’s DA algorithm for the logistic regression model is based on the joint
density $\pi(\beta,y\mid z)$. We now derive the conditional densities,
$\pi(\beta|y,z)$ and $\pi(y|\beta,z)$. From (2.17) we see that

$\pi(y_{i}\mid\beta,z)\propto\exp(-y_{i}(w_{i}^{\top}\beta)^{2}/2)p(y_{i}|\ell_{i},0),$
(2.18)

and thus from (2.16) we have
$y_{i}|\beta,z\overset{ind}{\sim}\text{PG}\left(\ell_{i},\left|w_{i}^{\top}\beta\right|\right),$
for $i=1,\dots,m$. Let $W$ denote the $m\times p$ design matrix with $i$th row
$w_{i}^{\top}$ and $Y$ be the $m\times m$ diagonal matrix with $i$th diagonal
element $y_{i}$. From (2.17) and (2.13), it follows that the conditional
density of $\beta$ is

$\pi\left(\beta|y,z\right)\propto\exp\left[-\frac{1}{2}\beta^{\top}(W^{\top}YW+Q)\beta+\beta^{\top}(W^{\top}\kappa+Q\mu_{0})\right],$

where $\kappa=\left(\kappa_{1},\dots,\kappa_{m}\right)^{\top}$. Thus the
conditional distribution of $\beta$ is multivariate normal. In particular,

$\beta|y,z\sim
N\left(\left(W^{\top}YW+Q\right)^{-1}(W^{\top}\kappa+Q\mu_{0}),\left(W^{\top}YW+Q\right)^{-1}\right).$
(2.19)

Thus, a single iteration of PS&W’s DA algorithm uses the following two steps
to move from $\beta$ to $\beta^{\prime}$.

Algorithm One iteration of the PS&W’s DA algorithm

1: Given $\beta$, draw $y_{1},\dots,y_{m}$ independently with
$y_{i}\sim\text{PG}\left(\ell_{i},\left|w_{i}^{\top}\beta\right|\right)$.

2: Draw $\beta^{\prime}\sim
N\left(\left(W^{\top}YW+Q\right)^{-1}(W^{\top}\kappa+Q\mu_{0}),\left(W^{\top}YW+Q\right)^{-1}\right)$.

PS&W’s DA algorithm is also applicable to the situation when $Q=0$ in (2.13),
that is, $\pi(\beta)\propto 1$, the improper uniform prior on $\beta$ provided
the posterior density $\pi(\beta|z)$ is proper, that is $c(z)$ defined in
(2.15) is finite. Polson et al., (2013) demonstrate superior empirical
performance of their DA algorithm over some other DA and MCMC algorithms for
the logistic regression model.

###  2.3 Data augmentation for probit mixed models

Generalized linear mixed models (GLMMs) are often used for analyzing
correlated binary observations. The random effects in the linear predictor of
a GLMM can accommodate for overdispersion as well as dependence among
correlated observations arising from longitudinal or repeated measures
studies. Let $(Z_{1},Z_{2},\dots,Z_{m})$ denote the vector of Bernoulli random
variables. Let $w_{i}$ and $v_{i}$ be the $p\times 1$ and $q\times 1$ known
covariates and random effect design vectors, respectively associated with the
$i$th observation $Z_{i}$ for $i=1,\dots,m$. Let $\beta\in\mathbb{R}^{p}$ be
the unknown vector of regression coefficients and $u\in\mathbb{R}^{q}$ be the
random effects vector. The probit GLMM connects the expectation of $Z_{i}$
with $w_{i}$ and $v_{i}$ using the probit link function, $\Phi^{-1}(\cdot)$ as

$\Phi^{-1}(P(Z_{i}=1))=w_{i}^{\top}\beta+v_{i}^{\top}u,$ (2.20)

where $\Phi(\cdot)$ is the cumulative distribution function of the standard
normal random variable.

Assume that we have $r$ random effects with
$u=(u_{1}^{\top},\dots,u_{r}^{\top})^{\top}$, where $u_{j}$ is a $q_{j}\times
1$ vector with $q_{j}>0$, $q_{1}+\cdots+q_{r}=q$, and
$u_{j}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}N(0,\Lambda_{j}\otimes
R_{j})$ where the low-dimensional covariance matrix $\Lambda_{j}$ is unknown
and must be estimated, and the structured matrix $R_{j}$ is usually known.
Here, $\otimes$ indicates the Kronecker product. Denoting
$\Lambda=(\Lambda_{1},\dots,\Lambda_{r})$, the probit GLMM is given by

$\displaystyle Z_{i}|\beta,u,\Lambda$
$\displaystyle\overset{\text{ind}}{\sim}$
$\displaystyle\text{Bern}(\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u))\text{ for
}i=1,\dots,m\;\text{ with}$ $\displaystyle u_{j}|\beta,\Lambda$
$\displaystyle\overset{\text{ind}}{\sim}$ $\displaystyle
N(0,\Lambda_{j}\otimes R_{j}),\,j=1,\dots,r.$

Let $z=(z_{1},z_{2},\dots,z_{m})^{\top}$ be the observed Bernoulli response
variables. Note that, the likelihood function for $(\beta,\Lambda)$ is

$\displaystyle
L(\beta,\Lambda|z)=\int_{\mathbb{R}^{q}}\prod_{i=1}^{m}\left[\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u)\right]^{z_{i}}\left[1-\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u)\right]^{1-z_{i}}\phi_{q}(u;0,A(\Lambda))du,$
(2.21)

which is not available in closed form. Here,
$A(\Lambda)=\oplus_{j=1}^{r}\Lambda_{j}\otimes R_{j}$, with $\oplus$
indicating the direct sum, and $\phi_{q}(s;a,B)$ denotes the probability
density function of the $q-$dimensional normal distribution with mean vector
$a$, covariance matrix $B$ and evaluated at $s$.

There are two widely used Monte Carlo approaches for approximating the
likelihood function (2.21) and making inference on $(\beta,\Lambda)$, namely
the Monte Carlo EM algorithm (Booth and Hobert,, 1999) and the Monte Carlo
maximum likelihood based on importance sampling (Geyer and Thompson,, 1992;
Geyer,, 1994). As explained in Roy, (2022), both the Monte Carlo EM and the
Monte Carlo maximum likelihood methods for making inference on
$(\beta,\Lambda)$ require effective methods for sampling from the conditional
density of the random effect $u$

$f(u|\beta,\Lambda,z)=\frac{f(z,u|\beta,\Lambda)}{L(\beta,\Lambda|z)},$ (2.22)

where $f(z,u|\beta,\Lambda)$ is the joint density of $(z,u)$ given by

$f(z,u|\beta,\Lambda)=\Bigg{[}\prod_{i=1}^{m}\left[\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u)\right]^{z_{i}}\left[1-\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u)\right]^{1-z_{i}}\Bigg{]}\phi_{q}(u;0,A(\Lambda)).$
(2.23)

Since the likelihood function $L(\beta,\Lambda|z)$ is not available in closed
form, neither is the density (2.22).

Albert and Chib,’s (1993) DA algorithm for the probit regression model is one
of the most widely used DA algorithms and it can be extended to construct DA
samplers for (2.22). Following Albert and Chib, (1993), let
$y_{i}\in\mathbb{R}$ be the latent continuous normal variable corresponding to
the $i$th binary observation $z_{i}$, that is $z_{i}=I(y_{i}>0)$, where
$y_{i}|\beta,u,\tau\overset{\text{ind}}{\sim}N(w_{i}^{\top}\beta+v_{i}^{\top}u$,
1) for $i=1,\dots,m$. Then

$P(Z_{i}=1)=P(Y_{i}>0)=\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u),$ (2.24)

that is,
$Z_{i}|\beta,u,\tau\overset{\text{ind}}{\sim}\text{Bern}(\Phi(x_{i}^{\top}\beta+v_{i}^{\top}u))$
as in (2.20). Using the latent variables $y=(y_{1},y_{2},\dots,y_{m})$, we now
introduce the joint density

$\displaystyle f(u,y|\beta,\Lambda,z)$
$\displaystyle=\Bigg{[}\prod_{i=1}^{m}\exp\left\\{-\frac{1}{2}\left(y_{i}-w_{i}^{\top}\beta-
v_{i}^{\top}u\right)^{2}\right\\}\times\prod_{i=1}^{m}\left[1_{\text{(0,$\infty$)}}\left(y_{i}\right)\right]^{z_{i}}\left[1_{\left(-\infty,0\right]}\left(y_{i}\right)\right]^{1-z_{i}}\Bigg{]}$
$\displaystyle\hskip
72.26999pt\times\frac{\phi_{q}(u;0,A(\Lambda))}{L(\beta,\Lambda|z)}.$ (2.25)

From (2.21) and (2.24) it follows that

$\int_{\mathbb{R}^{m}}f(u,y|\beta,\Lambda,z)dy=f(u|\beta,\Lambda,z),$ (2.26)

that is, the $u-$ marginal of the joint density $f(u,y|\beta,\Lambda,z)$ is
the target density $f(u|\beta,\Lambda,z)$ given in (2.22). Thus, the augmented
data $y$ and the joint density (2.3) satisfy the first property for
constructing a valid DA algorithm for $f(u|\beta,\Lambda,z)$.

If the two conditionals of the joint density (2.3) are easy to sample from,
then a DA algorithm can be constructed. From (2.3), we see that

$y_{i}|u,\beta,\Lambda,z\overset{\text{ind}}{\sim}\text{TN}(w_{i}^{\top}\beta+v_{i}^{\top}u,1,z_{i}),\,i=1,\dots,m,$
(2.27)

where $\text{TN}(\mu,\sigma^{2},e)$ denotes the distribution of the normal
random variable with mean $\mu$ and variance $\sigma^{2}$, that is truncated
to have only positive values if $e=1$, and only negative values if $e=0$. Let
$W$ and $V$ denote the $m\times p$ and $m\times q$ matrices whose $i$th rows
are $w^{\top}_{i}$ and $v^{\top}_{i}$, respectively. Next, using standard
linear models-type calculations, Roy, (2022) shows that the conditional
distribution of $u$ is

$u|y,\beta,\Lambda,z\sim
N_{q}\left((V^{\top}V+A(\Lambda)^{-1})^{-1}V^{\top}(y-W\beta),(V^{\top}V+A(\Lambda)^{-1})^{-1}\right).$
(2.28)

Thus, a single iteration of the DA algorithm uses the following two steps to
move from $u$ to $u^{\prime}$.

Algorithm One iteration of the DA algorithm for the probit GLMM

1: Given $u$, draw
$y_{i}\overset{\text{ind}}{\sim}\text{TN}(w_{i}^{\top}\beta+v_{i}^{\top}u,1,z_{i})$
for $i=1,\dots,m$.

2: Draw $u^{\prime}\sim
N_{q}\left((V^{\top}V+A(\Lambda)^{-1})^{-1}V^{\top}(y-W\beta),(V^{\top}V+A(\Lambda)^{-1})^{-1}\right)$.

One can consider a Bayesian analysis of the probit GLMM with appropriate
priors on $(\beta,\Lambda)$ and the DA and Haar PX-DA algorithms discussed in
this section can be extended to sample from the corresponding posterior
densities (see Wang and Roy, 2018b, ; Roy,, 2022). Similarly, PS&W’s DA
algorithm for the logistic model discussed in Section 2.2 can be extended to
fit logistic mixed models (see Wang and Roy, 2018a, ; Rao and Roy,, 2021;
Roy,, 2022).

##  3 Improving the DA algorithm

DA algorithms, although popular MCMC schemes, can suffer from slow convergence
to stationarity (Roy and Hobert,, 2007; Hobert et al.,, 2011; Duan et al.,,
2018). In the literature, a lot of effort has gone into developing methods for
accelerating the speed of convergence of the DA algorithms. These methods are
mainly based on techniques like appropriate reparameterizations and parameter
expansions such as the centered and noncentered parameterizations
(Papaspiliopoulos and Roberts,, 2008), interweaving the two parameterizations
(Yu and Meng,, 2011), parameter expanded-data augmentation (PX-DA) (Liu and
Wu,, 1999), the marginal augmentation (Meng and van Dyk,, 1999) or the
calibrated DA (Duan et al.,, 2018). Unlike the other methods mentioned here,
Duan et al.,’s (2018) calibrated DA requires an additional Metropolis-Hastings
(MH) step to maintain the stationarity of the algorithm with respect to the
target density $f_{X}$. The PX-DA and marginal augmentation are among the most
popular techniques for speeding up DA algorithms, and we describe those in
more details later in this section.

Graphically, the two steps of the DA algorithm can be viewed as $x\rightarrow
y\rightarrow x^{\prime}$. It has been observed that the convergence behavior
of a DA algorithm can be significantly improved by inserting a properly chosen
extra step in between the two steps of the DA algorithm. The idea of
introducing an extra step using appropriate auxiliary variables was developed
independently by Liu and Wu, (1999), who called the resulting MCMC algorithm
the PX-DA, and Meng and van Dyk, (1999), who called it the marginal
augmentation. Following these works, Hobert and Marchev, (2008) developed
general versions of PX-DA type algorithms and Yu and Meng, (2011) referred to
these generalized methods as sandwich algorithms since the algorithms involve
an additional step that is sandwiched between the two steps of the original DA
algorithm. Graphically, a sandwich algorithm can be represented as
$x\rightarrow y\rightarrow y^{\prime}\rightarrow x^{\prime}$, where the first
and the third steps are the two steps of the DA algorithm. Suppose the middle
step $y\rightarrow y^{\prime}$ is implemented by making a draw according to a
Mtd $r(y^{\prime}|y)$. Thus, denoting the current state by $x$, a generic
sandwich algorithm uses the following three steps to move to the new state
$x^{\prime}$.

Algorithm An iteration of a generic sandwich algorithm

1: Given $x$, draw $Y\sim f_{Y|X}(\cdot|x)$, and call the observed value $y$.

2: Draw $Y^{\prime}\sim r(y^{\prime}|y)$, and call the result $y^{\prime}$.

3: Draw $X^{\prime}\sim f_{X|Y}(\cdot|y^{\prime})$.

We now describe how to construct an Mtd $r$ for a valid middle step
$y\rightarrow y^{\prime}$ for sandwich algorithms. Recall that $f_{Y}(y)$ is
the marginal density of $y$. Suppose $r(y^{\prime}|y)$ leaves $f_{Y}$
invariant, that is,

$\int_{\mathbb{R}^{q}}r(y^{\prime}|y)f_{Y}(y)\,dy=f_{Y}(y^{\prime})\;.$ (3.1)

It turns out that any $r$ satisfying this invariance condition implies that
$f_{X}$ is stationary for the corresponding sandwich algorithm and
consequently can be used to obtain approximate samples from $f_{X}$. To see
it, note that the Mtd of the sandwich algorithm is given by

$k_{\text{SA}}(x^{\prime}|x)=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y^{\prime}|y)\,f_{Y|X}(y|x)\,dy\,dy^{\prime}\;,$
(3.2)

and

$\displaystyle\int_{\mathbb{R}^{p}}k_{\text{SA}}(x^{\prime}|x)f_{X}(x)\,dx$
$\displaystyle=\int_{\mathbb{R}^{p}}\bigg{[}\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y^{\prime}|y)\,f_{Y|X}(y|x)\,dy\,dy^{\prime}\bigg{]}f_{X}(x)\,dx$
$\displaystyle=\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\bigg{[}\int_{\mathbb{R}^{q}}r(y^{\prime}|y)\,f_{Y}(y)\,dy\bigg{]}\,dy^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,f_{Y}(y^{\prime})\,dy^{\prime}$
$\displaystyle=f_{X}(x^{\prime}),$ (3.3)

where the third equality follows from (3.1). Furthermore, if $r(y^{\prime}|y)$
satisfies the detailed balance condition with respect to $f_{Y}(y)$, that is,
$r(y^{\prime}|y)f_{Y}(y)$ is symmetric in $(y,y^{\prime})$, then

$\displaystyle k_{\text{SA}}(x^{\prime}|x)f_{X}(x)$
$\displaystyle=f_{X}(x)\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y^{\prime}|y)\,f_{Y|X}(y|x)\,dy\,dy^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y^{\prime}|y)\,f(x,y)\,dy\,dy^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y^{\prime}|y)f_{Y}(y)\,f_{X|Y}(x|y)\,\,dy\,dy^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x^{\prime}|y^{\prime})\,r(y|y^{\prime})f_{Y}(y^{\prime})\,f_{X|Y}(x|y)\,\,dy\,dy^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f(x^{\prime},y^{\prime})\,r(y|y^{\prime})\,f_{X|Y}(x|y)\,\,dy\,dy^{\prime}$
$\displaystyle=f_{X}(x^{\prime})\int_{\mathbb{R}^{q}}\int_{\mathbb{R}^{q}}f_{X|Y}(x|y)\,r(y|y^{\prime})\,f_{Y|X}(y^{\prime}|x^{\prime})\,\,dy\,dy^{\prime}$
$\displaystyle=k_{\text{SA}}(x|x^{\prime})f_{X}(x^{\prime})\;,$ (3.4)

that is, the corresponding sandwich algorithm is reversible with respect to
the target density $f_{X}$.

Let us denote the Markov chain underlying the sandwich algorithm by
$\\{X_{n}^{*}\\}_{n\geq 1}$. Since from (3) or (3) it follows that the target
density $f_{X}$ is stationary for $\\{X_{n}^{*}\\}_{n\geq 1}$, the sample
averages $\bar{h}^{*}_{n}:=\sum_{i=1}^{n}h(X^{*}_{i})/n$ based on the SA
Markov chain can also be used to estimate $E_{f_{X}}h$. Indeed, as mentioned
in Section 1, if $\\{X_{n}^{*}\\}_{n\geq 1}$ is appropriately irreducible,
$\bar{h}^{*}_{n}\rightarrow\mbox{E}_{f_{X}}h$ almost surely as
$n\rightarrow\infty$. Now, $\bar{h}^{*}_{n}$ will be preferred over
$\bar{h}_{n}$ as an estimate of $\mbox{E}_{f_{X}}h$ if the former leads to
smaller (Monte Carlo) errors. In particular, if there is a Markov chain
central limit theorem (CLT) for $\bar{h}_{n}$, that is, if
$\sqrt{n}(\overline{h}_{n}-E_{f_{X}}h)\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,\sigma_{h}^{2})$ for some finite asymptotic
variance $\sigma_{h}^{2}$, the standard error of $\overline{h}_{n}$ is
$\hat{\sigma}_{h}/\sqrt{n}$ where $\hat{\sigma}_{h}$ is it consistent
estimator of $\sigma_{h}$ (see Flegal and Jones,, 2010; Vats et al.,, 2019,
for methods of constructing $\hat{\sigma}_{h}$). Thus, in practice,
establishing CLTs for $\overline{h}_{n}$ and $\overline{h}^{*}_{n}$ is
important for ascertaining and comparing the quality of these estimators. A
sufficient condition extensively used in the literature for guaranteeing the
existence of such a CLT for every square integrable function $h$
($\int_{\mathbb{R}^{p}}h^{2}(x)\,f_{X}(x)\,dx<\infty$) is that the Markov
chain converges to its stationary distribution at a geometric rate (Roberts
and Rosenthal,, 1997; Chan and Geyer,, 1994).

Formally, the chain $\\{X_{n}\\}_{n\geq 1}$ is called geometrically ergodic if
there exist a function $M:\mathbb{R}^{p}\rightarrow[0,\infty)$ and a constant
$\rho\in[0,1)$ such that, for all $x\in\mathbb{R}^{p}$ and all $n=1,2,\dots$,

$\int_{\mathbb{R}^{p}}|k^{n}(x^{\prime}|x)-f_{X}(x)|dx\leq M(x)\,\rho^{n}\;,$

where $k^{n}(\cdot|\cdot)$ is the $n-$step Mtd of the chain
$\\{X_{n}\\}_{n\geq 1}$. Unfortunately, the simple irreducibility and
aperiodicity conditions mentioned in Section 1 does not imply geometric
ergodicity. The most common method of proving geometric ergodicity of
$\\{X_{n}\\}_{n\geq 1}$ is by establishing drift and minorization conditions
(Rosenthal,, 1995; Jones and Hobert,, 2001). Indeed, several DA algorithm
examples considered in this chapter have been shown to be geometrically
ergodic by establishing appropriate drift and minorization conditions. For
example, for the Bayesian lasso model discussed in Section 2.1, Rajaratnam et
al., (2019) established geometric ergodicity of the DA Markov chain in the
special case when $\alpha=\xi=0$. The geometric rate of convergence of PS&W’s
DA Markov chain for the binary logistic regression model discussed in Section
2.2 has been established by Choi and Hobert, (2013) and Wang and Roy, 2018c
under proper normal and improper uniform prior on $\beta$, respectively.

If the SA Markov chain $\\{X_{n}^{*}\\}_{n\geq 1}$ is reversible, denoting the
asymptotic variance for $\overline{h}^{*}_{n}$ by $\sigma_{h,*}^{2}$, it turns
out that $\sigma_{h,*}^{2}\leq\sigma_{h}^{2}$ for every square integrable
(with respect to $f_{X}$) function $h$ (Hobert and Marchev,, 2008). Thus,
sandwich algorithms are asymptotically more efficient than the corresponding
DA algorithms. So if both the DA and sandwich algorithms are run for the same
number of iterations, then the errors in $\bar{h}^{*}_{n}$ are expected to be
smaller than that of $\bar{h}_{n}$. In other words, to achieve the same level
of precision, the sandwich algorithm is expected to require fewer iterations
than the DA algorithm. In Section 4, we discuss other results established in
the literature showing the superiority of the sandwich algorithms over the
corresponding DA algorithms.

Intuitively, the extra step $R$ reduces the correlation between $x$ and
$x^{\prime}$ and thus improves the mixing of the DA algorithm. On the other
hand, since the extra step in a sandwich algorithm involves more computational
effort compared to the DA algorithm, any gain in mixing should be weighed
against this increased computational burden. Fortunately, there are several
examples where the sandwich algorithm provides a “free lunch” in that it
converges much faster than the underlying DA algorithm while requiring a
similar amount of computational effort per iteration (see e.g. Laha et al.,,
2016; Pal et al.,, 2015; Roy and Hobert,, 2007; Hobert et al.,, 2011; Roy,,
2014). Following Liu and Wu, (1999) and Meng and van Dyk, (1999), the Mtd $r$
for these efficient sandwich algorithms are often constructed by introducing
extra parameters $g\in{\mathsf{G}}\subset\mathbb{R}^{d}$ and a class of
functions $t_{g}:\mathbb{R}^{q}\rightarrow\mathbb{R}^{q}$ indexed by $g$. The
middle step $y\rightarrow y^{\prime}$ is implemented by drawing $g$ from a
density depending on $y$ and then setting $y^{\prime}=t_{g}(y)$ maintaining
invariance with respect to $f_{Y}$. Typically, $d$ is small, say 1 or 2, and
$\\{t_{g}(y):g\in G\\}$ is a small subspace of $\mathbb{R}^{q}$. Thus, drawing
from such an $r$ is usually much less computationally expensive than the two
steps of the DA Algorithm. However, as mentioned before, even when $d=1$, this
perturbation $y\rightarrow y^{\prime}$ often greatly improves the mixing of
the DA algorithm.

We now describe Liu and Wu,’s (1999) Haar PX-DA algorithm where the set
${\mathsf{G}}$ is assumed to possess a certain group structure (see also
Hobert and Marchev,, 2008). In particular, ${\mathsf{G}}$ is assumed to be a
topological group; that is, a group such that the functions
$(g_{1},g_{2})\mapsto g_{1}g_{2}$ and the inverse map $g\mapsto g^{-1}$ are
continuous. An example of such a ${\mathsf{G}}$ is the multiplicative group,
$\mathbb{R}_{+}$, where the binary operation defining the group is
multiplication, the identity element is 1, and $g^{-1}=1/g$. Suppose,
$t_{g}(y)$ represents ${\mathsf{G}}$ acting topologically on the left of
$\mathbb{R}^{q}$ (Eaton,, 1989, Chapter 2), that is, $t_{e}(y)=y$ for all
$y\in\mathbb{R}^{q}$ where $e$ is the identity element in ${\mathsf{G}}$ and
$t_{g_{1}g_{2}}(y)=t_{g_{1}}\big{(}t_{g_{2}}(y)\big{)}$ for all
$g_{1},g_{2}\in{\mathsf{G}}$ and all $y\in\mathbb{R}^{q}$.

We assume the existence of a multiplier
$\chi:{\mathsf{G}}\rightarrow\mathbb{R}_{+}$, that is, $\chi$ is continuous
and $\chi(g_{1}g_{2})=\chi(g_{1})\chi(g_{2})$ for all
$g_{1},g_{2}\in{\mathsf{G}}$. Assume that the Lebesgue measure on
$\mathbb{R}^{q}$ is relatively (left) invariant with respect to the multiplier
$\chi$; that is, assume that for any $g\in{\mathsf{G}}$ and any integrable
function $\zeta:\mathbb{R}^{q}\rightarrow\mathbb{R}$, we have

$\chi(g)\int_{\mathbb{R}^{q}}\zeta\big{(}t_{g}(y)\big{)}\,dy=\int_{\mathbb{R}^{q}}\zeta(y)\,dy\;.$

Next, suppose the group ${\mathsf{G}}$ has a left-Haar measure of the form
$\nu_{l}(g)\,dg$ where $dg$ denotes Lebesgue measure on ${\mathsf{G}}$. It is
known that the left-Haar measure satisfies

$\int_{\mathsf{G}}h\big{(}\tilde{g}g\big{)}\,\nu_{l}(g)\,dg=\int_{\mathsf{G}}h(g)\,\nu_{l}(g)\,dg$
(3.5)

for all $\tilde{g}\in{\mathsf{G}}$ and all integrable functions
$h:{\mathsf{G}}\rightarrow\mathbb{R}$. Finally, assume that

$q(y):=\int_{\mathsf{G}}f_{Y}\big{(}t_{g}(y)\big{)}\,\chi(g)\,\nu_{l}(g)\,dg$

is strictly positive for all $y\in\mathbb{R}^{q}$ and finite for (almost) all
$y\in\mathbb{R}^{q}$.

We now state how the middle step $y\rightarrow y^{\prime}$ is implemented in
Liu and Wu,’s (1999) Haar PX-DA algorithm. If the current state of the Haar
PX-DA chain is $x$, an iteration of this chain uses the following steps to
move to $x^{\prime}$.

Algorithm An iteration for the Haar PX-DA Algorithm

1: Given $x$, draw $Y\sim f_{Y|X}(\cdot|x)$, and call the observed value $y$.

2: Draw $G$ from the density proportional to
$f_{Y}\big{(}t_{g}(y)\big{)}\,\chi(g)\,\nu_{l}(g)$, call the result $g$, and
set $y^{\prime}=t_{g}(y)$.

3: Draw $X^{\prime}\sim f_{X|Y}(\cdot|y^{\prime})$.

######  Example 1 (continues=sec:probit).

Improving the DA algorithm for the probit mixed model presented in Section
2.3, Roy, (2022) constructs a Haar PX-DA algorithm for (2.22). For
constructing the sandwich step $y\rightarrow y^{\prime}$, we need the marginal
density of $Y$, $f_{Y}(y|\beta,\Lambda,z)$ from the joint density
$f(u,y|\beta,\Lambda,z)$ given in (2.3). It turns out that

$\displaystyle
f_{Y}(y|\beta,\Lambda,z)\propto\prod_{i=1}^{m}\left[1_{\text{(0,$\infty$)}}\left(y_{i}\right)\right]^{z_{i}}\left[1_{\left(-\infty,0\right]}\left(y_{i}\right)\right]^{1-z_{i}}\exp\left\\{-\frac{1}{2}\left[y^{\top}V_{1}y-2y^{\top}V_{1}W\beta\right]\right\\},$
(3.6)

where

$V_{1}=\left[I_{m}-V\left(V^{\top}V+(A(\Lambda))^{-1}\right)^{-1}V^{\top}\right].$

Let $\mathcal{Y}$ denote the subset of $\mathbb{R}^{m}$ where $y$ lives, that
is, $\mathcal{Y}$ is the Cartesian product of $m$ half (positive or negative)
lines, where the $i$th component is $(0,\infty)$ (if $z_{i}=1$) or
$(-\infty,0]$ (if $z_{i}=0$).

As mentioned in this Section, take ${\mathsf{G}}=\mathbb{R}_{+}$ and let the
multiplicative group ${\mathsf{G}}$ act on $\mathcal{Y}$ through
$t_{g}(y)=gy=(gy_{1},\dots,gy_{m})$. Note that, for any
$\tilde{g}\in{\mathsf{G}}$, we have

$\int_{0}^{\infty}h\big{(}\tilde{g}g\big{)}\,\frac{1}{g}\,dg=\int_{0}^{\infty}h(g)\,\frac{1}{g}\,dg\;,$

which shows from (3.5) that $\nu(dg)=dg/g$ is a left-Haar measure for the
multiplicative group, where $dg$ is Lebesgue measure on $\mathbb{R}_{+}$.
Also, with the group action defined this way, it is known that the Lebesgue
measure on $\mathcal{Y}$ is relatively left invariant with the multiplier
$\chi(g)=g^{m}$ (Hobert and Marchev,, 2008). Thus,

$f_{Y}\left(gy|\beta,\Lambda,z\right)\chi\left(g\right)\nu(dg)\propto
g^{m-1}\exp\left\\{-\frac{1}{2}\left[g^{2}y^{\top}V_{1}y-2gy^{\top}V_{1}W\beta\right]\right\\}dg,$
(3.7)

and since $V_{1}$ is a positive definite matrix,

$q(y)=\int_{0}^{\infty}g^{m-1}\exp\left\\{-\frac{1}{2}\left[g^{2}y^{\top}V_{1}y-2gy^{\top}V_{1}W\beta\right]\right\\}dg$

is strictly positive for all $y\in\mathcal{Y}$ and finite for (almost) all
$y\in\mathcal{Y}$. Thus, a single iteration of the Haar PX-DA algorithm for
the probit mixed model uses the following three steps to move from $u$ to
$u^{\prime}$.

Algorithm One iteration of the Haar PX-DA algorithm for the probit GLMM

1: Given $u$, draw
$y_{i}\overset{\text{ind}}{\sim}\text{TN}(w_{i}^{\top}\beta+v_{i}^{\top}u,1,z_{i})$
for $i=1,\dots,m$.

2: Draw $g$ from (3.7).

3: Calculate $y_{i}^{\prime}=gy_{i}$ for $i=1,\dots,m$, and draw $u^{\prime}$
from (2.28) conditional on
$y^{\prime}=(y_{1}^{\prime},\dots,y_{m}^{\prime})^{\top}$, that is, draw
$u^{\prime}\sim
N_{q}\left((V^{\top}V+(A(\Lambda))^{-1})^{-1}V^{\top}(y^{\prime}-W\beta),(V^{\top}V+A(\Lambda)^{-1})^{-1}\right).$

The density (3.7) is log-concave, and Roy, (2022) suggests using Gilks and
Wild,’s (1992) adaptive rejection sampling algorithm to sample from it. Since
a single draw from (3.7) is the only difference between each iteartion of the
DA and the Haar PX-DA algorithms, the computational burden for the two
algorithms is similar.

Given the group ${\mathsf{G}}$, the Haar PX-DA algorithm is the best sandwich
algorithm in terms of efficiency and operator norm (Hobert and Marchev,,
2008). Khare and Hobert, (2011) established necessary and sufficient
conditions for the Haar PX-DA algorithm to be strictly better than the
corresponding DA algorithm (see Section 4 for details). Also, the Haar PX-DA
method with appropriate choices of ${\mathsf{G}}$ and $t_{g}$ has been
empirically demonstrated to result in a huge improvement in the mixing of DA
algorithms albeit with roughly the same computational cost in various examples
(see e.g. Roy and Hobert,, 2007; Liu and Wu,, 1999; Meng and van Dyk,, 1999;
Roy,, 2014; Pal et al.,, 2015). Finally, there are other strategies proposed
in the literature to improve DA algorithms without introducing extra auxiliary
variables as in the sandwich algorithms. For example, Roy, (2016) replaces one
of the two original steps of the DA algorithm with a draw from the MH step
given in Liu, (1996). One advantage of Roy,’s (2016) modified DA algorithm
over the sandwich algorithms is that the Markov transition matrix of the
modified algorithm can have negative eigenvalues, whereas the DA algorithms
and also, generally, the sandwich algorithms used in practice are positive
Markov chains leading to positive eigenvalues (Khare and Hobert,, 2011).
Hence, Roy,’s (2016) modified DA algorithm can lead to superior performance
over the DA and sandwich algorithms in terms of asymptotic variance.

##  4 Spectral properties of the DA Markov chain

The convergence of the Markov chains associated with MCMC algorithms to the
desired stationary distribution is at the heart of the popularity of MCMC
technology. In Section 1 we discussed the convergence of the sample averages
$\bar{h}_{n}$ and the marginal distributions of $X_{n}$. For a serious
practitioner, understanding the quality of both the convergences is important.
As shown in Section 3, proving geometric ergodicity of the Markov chain is key
for providing standard errors for $\bar{h}_{n}$. Establishing geometric
ergodicity for general state space Markov chains arising in statistical
practice is in general quite challenging. The drift and minorization technique
that is often used to prove geometric ergodicity has two flavors/tracks. The
first track leads to establishing qualitative geometric ergodicity, without
providing quantitative convergence bounds for distance to stationarity after
finitely many Markov chain iterations. The other track leads to explicit
convergence bounds, but these bounds are often too conservative to be useful,
especially in modern high-dimensional settings Jones and Hobert, (2004); Qin
and Hobert, (2021). While several alternate and promising methods for
obtaining convergence bounds have been developed in recent years (see Qin and
Hobert, (2022) for a useful summary), it is safe to say that a successful
geometric ergodicity analysis remains elusive for an overwhelming majority of
MCMC algorithms used in statistical practice. However, as we will see below,
the special structure/construction of the DA Markov chain often allows us in
many examples to establish stronger properties which imply geometric
ergodicity in particular, and lead to a richer and deeper understanding of the
structure and convergence properties of the relevant DA Markov chain.

###  4.1 Leveraging operator theory for comparing DA and sandwich chains

Operator theory and associated tools play a key role in convergence analysis
and comparison of the DA Markov chain (with Mtd $k$) and the sandwich Markov
chain (with Mtd $k_{SA}$). Let $L_{0}^{2}(f_{X})$ denote the space of real-
valued functions (with domain $\mathbb{R}^{p}$) that have mean zero and are
square integrable with respect to $f_{X}$. In particular,

$L_{0}^{2}(f_{X})=\left\\{h:\mathbb{R}^{p}\rightarrow\mathbb{R}:\;\int_{\mathbb{R}^{p}}h(x)^{2}f_{X}(x)dx<\infty\mbox{
and}\int_{\mathbb{R}^{p}}h(x)f_{X}(x)dx=0\right\\}.$

Let $K:L_{0}^{2}(f_{X})\rightarrow L_{0}^{2}(f_{X})$ and
$K_{SA}:L_{0}^{2}(f_{X})\rightarrow L_{0}^{2}(f_{X})$ denote the Markov
operators corresponding to the Markov transition densities $k$ and $k_{SA}$
respectively. Then for any $h\in L_{0}^{2}(f_{X})$, we have

$(Kh)(x)=\int_{\mathbb{R}^{p}}h(x^{\prime})k(x^{\prime}\mid
x)dx^{\prime}\mbox{ and
}(K_{SA}h)(x)=\int_{\mathbb{R}^{p}}h(x^{\prime})k_{SA}(x^{\prime}\mid
x)dx^{\prime}.$

Note that $L_{0}^{2}(f_{X})$ is a Hilbert space equipped with the inner
product

$<h_{1},h_{2}>=\int_{\mathbb{R}^{p}}h_{1}(x)h_{2}(x)f_{X}(x)dx$

for every pair $h_{1},h_{2}\in L_{0}^{2}(f_{X})$, and the corresponding norm
defined by $\|h\|=\sqrt{<h,h>}$ for every $h\in L_{0}^{2}(f_{X})$. The
operator norms of $K$ and $K_{SA}$, denoted respectively by $\|K\|$ and
$\|K_{SA}\|$, are defined as

$\|K\|=\sup_{h\in L_{0}^{2}(f_{X}),\|h\|=1}\|Kh\|\mbox{ and
}\|K_{SA}\|=\sup_{h\in L_{0}^{2}(f_{X}),\|h\|=1}\|K_{SA}h\|.$

The convergence properties of the DA and sandwich Markov chains are
intricately linked to the behavior of the corresponding Markov operators $K$
and $K_{SA}$. Using the structure of $k$ and $k_{SA}$, it can be shown easily
that the operator $K$ is a non-negative operator (an operator $P$ is non-
negative if $<Pg,g>$ is non-negative $\forall g$), $\|K\|\in[0,1]$ and
$\|K_{SA}\|\in[0,1]$. Also, $K$ is a self-adjoint operator, and if
$r(y^{\prime}\mid y)$ satisfies detailed balance with respect to $f_{Y}(y)$
(which will be assumed henceforth, unless specified otherwise) then so is
$K_{SA}$ (an operator $P$ is self-adjoint if
$<Ph_{1},h_{2}>=<Ph_{2},h_{1}>\;\forall h_{1},h_{2}$). Roberts and Rosenthal,
(1997) show that for a self-adjoint Markov operator $P$ corresponding to a
Harris ergodic chain, $\|P\|<1$ if and only if the corresponding Markov chain
is geometrically ergodic. In this case, $\|P\|$ in fact corresponds to the
asymptotic rate of convergence of the relevant Markov chain to its stationary
distribution. The above concepts enable us to state the first result that
rigorously shows that the sandwich chain is at least as good as the DA chain
in some key aspects.

######  Theorem 4.1 (Hobert and Marchev, (2008)).

If the DA chain with Mtd $k$ is Harris ergodic, and the sandwich chain (with
Mtd $k_{SA}$) can itself be interpreted as a DA chain, then
$\|K_{SA}\|\leq\|K\|$. Hence, geometric Harris ergodicity of the DA chain
implies geometric Harris ergodicity of the sandwich chain, and in such
settings the asymptotic rate of convergence to stationarity for the sandwich
chain is at least as good as that of the DA chain.

The structure and analysis of sandwich chains is more complicated (sometimes
significantly so) than the DA chains due to the additional sandwich step.
However, the above result allows us to completely avoid analyzing the sandwich
chain, if geometric ergodicity can be established for the corresponding DA
chain.

Intuition says that given the extra step, the sandwich chain should have
better mixing and faster convergence than the original DA chain under suitable
regularity conditions. The above result provides rigorous justification to
this intuition. See Section 4.4 for further discussion.

###  4.2 The trace class property for Markov operators

An operator $P$ on $L^{2}_{0}(f_{X})$ is defined to be Hilbert-Schmidt if for
any orthonormal sequence $\\{h_{n}\\}_{n\geq 1}$ in $L^{2}_{0}(f_{X})$, we
have

$\sum_{n=1}^{\infty}\|Ph_{n}\|^{2}<\infty.$ (4.1)

Additionally, a positive self-adjoint operator (such as the DA operator $K$)
is defined to be trace class if

$\sum_{n=1}^{\infty}<Ph_{n},h_{n}><\infty$ (4.2)

for any orthonormal sequence $\\{h_{n}\\}_{n\geq 1}$ in $L^{2}_{0}(f_{X})$.
Both the above properties imply that the operator $P$ is compact, and has
countably many singular values. The trace class property implies that these
singular values are summable, while the Hilbert-Schmidt property implies that
these singular values are square-summable. If the associated Markov chain is
Harris ergodic, then compactness of a Markov operator implies that its
spectral radius is less than $1$. Existing results (see for example
Proposition 2.1 and Remark 2.1 in Roberts and Rosenthal, (1997)) can now be
leveraged to establish geometric ergodicity. To summarize, for Harris ergodic
Markov chains, we have

$\mbox{Hilbert Schmidt }\Longrightarrow\mbox{Compact
}\Longrightarrow\mbox{Geometrically ergodic},$

and if the corresponding Markov operator is also positive and self-adjoint,
then

$\mbox{Trace Class }\Longrightarrow\mbox{Hilbert Schmidt
}\Longrightarrow\mbox{Compact }\Longrightarrow\mbox{Geometrically ergodic}.$

The above implications potentially offer an alternate mechanism/route for
establishing geometric ergodicity. However, the conditions in (4.1) and (4.2)
seem hard to verify especially for intractable Markov transition densities
arising in statistical applications. Fortunately, there are equivalent
characterizations for (4.1) and (4.2) that are easier to state and verify. In
particular, results in Davies, (1983) can be leveraged to show that a Markov
operator $P$ (with corresponding Mtd $p(\cdot,\cdot)$) on $L^{2}_{0}(f_{X})$
is Hilbert-Schmidt if and only if

$\int_{\mathbb{R}^{p}}\int_{\mathbb{R}^{p}}\frac{p(x,y)^{2}f_{X}(x)}{f_{X}(y)}dxdy=\int_{\mathbb{R}^{p}}\int_{\mathbb{R}^{p}}\left(\frac{p(x,y)}{f_{X}(y)}\right)^{2}f_{X}(x)f_{X}(y)dxdy<\infty,$
(4.3)

and in case $P$ is self-adjoint and positive, it is trace class if and only if

$\int_{\mathbb{R}^{p}}p(x,x)dx<\infty.$ (4.4)

The conditions in (4.3) and (4.4) while simpler than those in (4.1) and (4.2)
can still be quite challenging to verify for Markov chains arising in
statistical practice. The corresponding Markov transition densities are
themselves often expressed as intractable integrals, and the analysis of the
expressions in (4.3) and (4.4) can be quite lengthy, laborious and intricate.
As the authors state in Liu et al., (1995), the condition in (4.3) “is
standard, but not easy to check and understand”. Yet, there have been many
success stories in recent years, see Khare and Hobert, (2011); Pal et al.,
(2017); Chakraborty and Khare, (2017, 2019); Qin et al., (2019); Rajaratnam et
al., (2019); Jin and Tan, (2021). Markov chains used by statistical
practitioners overwhelmingly employ either the Gibbs sampling algorithm, or
the Metropolis-Hastings algorithm, or a combination of both. Any continuous
state space chain with a non-trivial Metropolis component cannot be compact
(Chan and Geyer,, 1994) and hence cannot be trace class. The same is true for
random scan Gibbs chains. (Systematic scan) Gibbs sampling Markov chains are
in general not reversible. However, the convergence of a two-block (systematic
scan) Gibbs chain is completely characterized by its “marginal chains” which
correspond to positive self- adjoint Markov operators (Liu et al.,, 1994;
Diaconis et al.,, 2008). In particular, the construction in Section 1 shows
that the DA chain is a marginal chain of a two-block Gibbs sampler for
$f_{X,Y}$. Hence, it is not surprising that all Markov chains which have been
shown to be trace class in these papers are DA Markov chains.

###  4.3 Trace class property as a tool to establish geometric ergodicity

We have already discussed how establishing the trace class property for a
positive self-adjoint Markov operator corresponding to a Harris ergodic Markov
chain establishes geometric ergodicity of that chain. In this section, we
compare the potential advantages and disadvantages of this alternative
approach compared to the drift and minorization approach for establishing
geometric ergodicity. A key advantage is that the trace class approach based
on (4.4) is more straightforward and streamlined than the drift and
minorization approach. Indeed, construction of effective drift functions
depends heavily on the Markov chain at hand and has been described as “a
matter of art” (Diaconis et al., (2008)). On the other hand, the drift and
minorization approach has broader applicability, and has been successfully
applied to non-reversible and non-compact chains.

Suppose that for a DA Markov chain one is able to establish geometric
ergodicity through a drift and minorization analysis, and is also able to
establish the trace class condition in (4.4). Since the trace class property
is much stronger than geometric ergodicity, one would expect that assumptions
needed for establishing the trace class property are stronger than those
needed for the drift and minorization analysis. Surprisingly, as demonstrated
in Mukherjee et al., (2023), this is not always true. To understand why this
might happen, we take a brief but careful look at the Markov chain analyzed in
Mukherjee et al., (2023) below.

Similar to Sections 2.2 and 2.3, consider a binary regression setting with $n$
independent binary responses $Z_{1},Z_{2},\cdots,Z_{n}$ and corresponding
predictor vectors ${\bf w}_{1},{\bf w}_{2},\cdots,{\bf
w}_{n}\in\mathbb{R}^{p}$, such that

$P(Z_{i}=1\mid{\boldsymbol{\beta}})=F_{\nu}({\bf
w}_{i}^{T}{\boldsymbol{\beta}}),$ (4.5)

for $1\leq i\leq n$. The goal is to estimate
${\boldsymbol{\beta}}\in\mathbb{R}^{p}$. Here $F_{\nu}$ is the CDF of the
Student’s $t$-distribution with $\nu$ degrees of freedom, and the
corresponding model is referred to as the robit regression model Liu, (2004).
In the binary regression setting, outliers are observations with unexpectedly
large predictor values and a misclassified response. The robit model can more
effectively down-weight outliers and produce a better fit than probit or
logistic models Pregibon, (1982).

Following Roy, 2012a ; Albert and Chib, (1993), consider a Bayesian model with
a multivariate normal prior for ${\boldsymbol{\beta}}$ with mean
${\boldsymbol{\beta}}_{a}$ and covariance matrix $\Sigma_{a}^{-1}$. Let ${\bf
z}\in\\{0,1\\}^{n}$ denote the vector of observed values for the response
variables $Z_{1},Z_{2},\cdots,Z_{n}$. As demonstrated in Roy, 2012a , the
posterior density $\pi({\boldsymbol{\beta}}\mid{\bf z})$ is intractable, and
Roy, 2012a  develops a Data DA approach to construct a computationally
tractable Markov chain with $\pi({\boldsymbol{\beta}}\mid{\bf y})$ as its
stationary density by introducing unobserved latent variables
$\\{(U_{i},\lambda_{i})\\}_{i=1}^{n}$ that are mutually independent and
satisfy $U_{i}\mid\lambda_{i}\sim\mathcal{N}({\bf
w}_{i}^{T}{\boldsymbol{\beta}},1/{\bf\lambda}_{i})$ and
$\lambda_{i}\sim\mbox{Gamma}(\nu/2,\nu/2)$. Straightforward calculations now
show that $U_{i}\sim t_{\nu}({\bf w}_{i}^{T}{\boldsymbol{\beta}},1)$, where
$t_{\nu}(\mu,\sigma)$ denotes the Student’s $t$-distribution with $\nu$
degrees of freedom, location $\mu$ and scale $\sigma$. Furthermore, if
$Y_{i}=1_{\\{Z_{i}>0\\}}$, then
$P(Y_{i}=1\mid{\boldsymbol{\beta}})=P(U_{i}>0)=F_{\nu}({\bf
w}_{i}^{T}{\boldsymbol{\beta}})$, which precisely corresponds to the robit
regression model specified above.

One can map this setting to the DA framework laid out in Section 1, by viewing
${\boldsymbol{\beta}}$ as “$X$”, $({\bf U},{\boldsymbol{\lambda}})$ as “$Y$”,
and the joint posterior density of ${\boldsymbol{\beta}},{\bf
U},{\boldsymbol{\lambda}}$ as “$f_{X,Y}$”. Straightforward calculations (see
Roy, 2012a ) show that samples from $f_{X\mid Y}$ and $f_{Y\mid X}$ are easy
to generate, as they involve sampling from standard distributions such as
multivariate normal, truncated-$t$ and Gamma. These observations allow Roy,
2012a  to use the corresponding DA Markov chain on $\mathbb{R}^{p}$ to
generate samples from the $f_{X}$ (the marginal posterior density of
${\boldsymbol{\beta}}$). We will refer to this Markov chain as the robit DA
Markov chain.

Roy, 2012a  established Harris ergodicity of the robit DA chain, and
investigated and established geometric ergodicity using a drift and
minorization analysis. However, this analysis requires the following
assumptions.

  * •

The design matrix $W$ is full rank (which implies $n\geq p$ and rules out
high-dimensional settings)

  * •

$\Sigma_{a}=g^{-1}W^{T}W$

  * •

$n\leq\frac{g^{-1}\nu}{(\nu+1)(1+2\sqrt{{\boldsymbol{\beta}}_{a}^{T}W^{T}W{\boldsymbol{\beta}}_{a}})}.$

The last upper bound on $n$ involving the design matrix, the prior mean and
covariance and the degrees of freedom $\nu$ is in particular very restrictive.
Through a tighter drift and minorization analysis, Mukherjee et al., (2023)
relax the above restrictions to some extent, but not substantially (see
Theorem S1 and Theorem S2 in Mukherjee et al., (2023)). The related probit DA
chain is obtained by (a) using the standard normal CDF $\Phi$ instead of
$F_{\nu}$ in (4.5), (b) using again a multivariate normal prior for
${\boldsymbol{\beta}}$, (c) introducing latent variables
$U_{i}\sim\mathcal{N}({\bf w}_{i}^{T}{\boldsymbol{\beta}},1)$ with
$Y_{i}=1_{\\{Z_{i}>0\\}}$ for $1\leq i\leq n$, and (d) employing the DA
machinery in Section 1 with $({\boldsymbol{\beta}},{\bf U})$ as “$(X,Y)$”. The
geometric ergodicity analysis of the probit DA chain in Roy and Hobert,
(2007); Chakraborty and Khare, (2017) uses similar drift functions as in Roy,
2012a ; Mukherjee et al., (2023) but requires minimal assumptions. Note that
the latent variables $\\{\lambda_{i}\\}_{i=1}^{n}$ used in the robit setting
are not needed in the probit setting. While having this additional layer of
latent variables definitely complicates the analysis of the robit DA chain, it
is not clear if the restrictive conditions listed above for geometric
ergodicity of the robit DA chain are really necessary or if they are an
artifact of using the drift and minorization technique in face of this added
complication.

A trace class analysis for the probit and robit DA chains is available in
Chakraborty and Khare, (2017) and Mukherjee et al., (2023) respectively. The
trace class property for the probit DA chain was established in (Chakraborty
and Khare,, 2017, Theorem 2) under some constraints on $W$ and the prior
covariance matrix $\Sigma_{a}$. This is consistent with the expectation that
weaker assumptions should be needed for establishing drift and minorization
based geometric ergodicity, as compared to establishing the trace class
property. However, for the robit DA chain the reverse phenomenon holds:
Mukherjee et al., (2023) establish the trace class property for the robit DA
chain for any $n,\nu>2,W,{\bf z},{\boldsymbol{\beta}}_{a},\Sigma_{a}$. Hence,
drift and minorization approach, with the drift functions chosen in Roy, 2012a
; Mukherjee et al., (2023), needs stronger conditions to succeed than the
trace class approach. Essentially, the additional layer of latent variables
$\lambda_{i}$ introduced in the robit setting severely hampers the drift and
minorization analysis, but does not cause additional complications for
establishing the trace class condition.

The key takeaway from the above discussion is that even if a successful drift
and minorization based analysis is available, it is worth investigating the
finiteness of the trace class integral in (4.4), it might require
comparatively weaker assumptions contrary to expectations. Also, establishing
the trace class property provides additional insights which are discussed in
the next subsection.

###  4.4 Trace class property as a tool to compare DA and sandwich chains

The results discussed previously show that, under certain regularity
conditions, the sandwich chain is at least as good as the DA chain in certain
respects. These results, while very useful, are not completely satisfactory as
they only establish that the sandwich chain is “as good as” the DA chain. The
question is - are there any conditions under which the sandwich chain can be
shown to be “strictly better” than the DA chain (in an appropriate sense)?
This question has been investigated in Khare and Hobert, (2011); Roy, 2012b ,
and we will discuss their results below.

First, we provide a brief review of some ideas and results used in these
analyses. Both papers leverage the fact (see Buja, (1990); Diaconis et al.,
(2008)) that the DA operator $K$ can be written as a product of two simple
‘projection’ operators. In particular, let $P_{X}$ denotes the operator from
$L^{2}_{0}(f_{Y})$ to $L^{2}_{0}(f_{X})$ defined by

$(P_{X}h)(x)=\int_{\mathbb{R}^{q}}h(y)f_{Y\mid X}(y\mid x)dy$

for every $h\in L^{2}_{0}(f_{Y})$, and $P_{Y}$ denotes the operator from
$L^{2}_{0}(f_{X})$ to $L^{2}_{0}(f_{Y})$ defined by

$(P_{Y}g)(y)=\int_{\mathbb{R}^{p}}g(x)f_{X\mid Y}(x\mid y)dx$

for every $g\in L^{2}_{0}(f_{X})$. Also, let $R:L^{2}_{0}(f_{Y})\rightarrow
L^{2}_{0}(f_{Y})$ denotes the operator corresponding to the Mtd $r$ (the
sandwich step Mtd). Then, it can be shown that

$K=P_{X}P_{Y}\quad\mbox{ and }\quad K_{SA}=P_{X}RP_{Y}.$

For any $g\in L^{2}_{0}(f_{X})$ and any $h\in L^{2}_{0}(f_{Y})$, simple
calculations show that $<P_{X}h,g>=<h,P_{Y}g>$ (the first inner product
corresponds to $L^{2}_{0}(f_{X})$, while the second one corresponds to
$L^{2}_{0}(f_{Y})$). Hence, the operators $P_{X}$ and $P_{Y}$ are adjoints of
each other. Using elementary operator theoretic arguments, this fact can be
used to immediately conclude that $K=P_{X}P_{Y}$ is a positive operator (as
previously mentioned) and $\|K\|=\|P_{X}\|^{2}=\|P_{Y}\|^{2}$. It follows that

$\|K_{SA}\|=\|P_{X}RP_{Y}\|\leq\|P_{X}\|\|R\|\|P_{Y}\|=\|R\|\|K\|.$

However, in almost all applications of interest, $\|R\|=1$, since the sandwich
step corresponds typically to a univariate movement in $\mathbb{R}^{q}$.
Hence, the above argument does not allow us to prove that the sandwich
algorithm is strictly better (in the norm sense).

On the other hand, if the trace class condition in (4.4) is satisfied, then a
very useful singular value decomposition of $f_{X,Y}$ becomes available (see
Buja, (1990)). In particular, it can be shown that

$\frac{f_{X,Y}(x,y)}{f_{X}(x)f_{Y}(y)}=\sum_{i=0}^{\infty}\beta_{i}g_{i}(x)h_{i}(y),$

where $\\{g_{i}\\}_{i=1}^{\infty}$ and $\\{h_{j}\\}_{j=1}^{\infty}$ form
orthonormal bases of $L^{2}_{0}(f_{X})$ and $L^{2}_{0}(f_{Y})$ respectively
(with $g_{0}\equiv 1,h_{0}\equiv 1$), $\\{\beta_{i}\\}_{i=0}^{\infty}$ is a
decreasing sequence of numbers in the unit interval with $\beta_{0}=1$, and
$g_{i}$ and $h_{j}$ are orthogonal for every $i\neq j$ in the sense that

$\int_{\mathbb{R}^{p}}\int_{\mathbb{R}^{q}}g_{i}(x)h_{j}(y)f_{X,Y}(x,y)dydx=0.$

The singular value decomposition can be used to establish that

$P_{Y}g_{i}=\beta_{i}h_{i}\quad\mbox{ and }P_{X}h_{i}=\beta_{i}g_{i}$ (4.6)

for every $i\geq 0$. In particular this implies
$\\{\beta_{i}^{2}\\}_{i=1}^{\infty}$ is the spectrum of the DA operator $K$
and $\\{g_{i}\\}_{i=1}^{\infty}$ are the corresponding eigenfunctions. Given
that $K_{SA}=P_{X}RP_{Y}$ and (4.6), a comparative spectral analysis for the
DA and sandwich algorithms now hinges on how the operator $R$ interacts with
the functions $\\{h_{i}\\}_{i=1}^{\infty}$. By investigating this carefully,
Khare and Hobert, (2011) establish the following result.

######  Theorem 4.2.

Assume that (4.4) holds, and that the sandwich operator $R$ is idempotent
($R^{2}=R$) with $\|R\|=1$. Then the sandwich operator $K_{SA}$ is positive
and trace class. If $\\{\lambda_{SA,i}\\}_{i=1}^{\infty}$ denotes the ordered
sequence of eigenvalues for $K_{SA}$, then $\lambda_{SA,i}\leq\beta_{i}^{2}$
for every $i\geq 1$. Also, for every $i$ such that $\beta_{i}>0$,
$\lambda_{SA,i}=\beta_{i}^{2}$ if and only if $Rh_{i}=h_{i}$.

In other words, the spectrum of the sandwich operator $K_{SA}$ is dominated
pointwise by the spectrum of the DA operator $K$, and strict domination for a
specific pair of eigenvalues depends exclusively on whether or not the
operator $R$ leaves the corresponding $h$-function invariant. A stronger
version of this result, which only requires $K$ to be compact is established
in Roy, 2012b .

For a more specific comparison of the operator norms of the two operators, as
expected the key factor is the interaction of the the operator $R$ with the
linear span of $\\{h_{i}\\}_{i=1}^{\ell}$, where $\ell$ is the multiplicity of
the largest eigenvalue $\beta_{1}$, i.e., number of $\beta_{i}$’s that are
exactly equal to $\beta_{1}$. The next result from Khare and Hobert, (2011)
leverages this idea to provide a necessary and sufficient condition for the
operator norm of $K_{SA}$ to be strictly smaller than the operator norm of
$K$.

######  Theorem 4.3.

Under the setup in Theorem 4.2, $\|K_{SA}\|<\|K\|$ if and only if the only
function in the linear span of $\\{h_{i}\\}_{i=1}^{\ell}$ which is left
invariant by $R$ is the (identically) zero function.

While the two results above provide necessary and sufficient conditions for
the equality of the norm (and other eigenvalues) for $K$ and $K_{SA}$, their
practical utility is somewhat limited by the requirement to identify the
functions $\\{h_{i}\\}_{i=1}^{\infty}$. However, if we restrict to the group
action based Haar PX-DA chain, then the following result can be established.

######  Theorem 4.4.

Assume that (4.4) holds, and that the operator $K_{SA}$ corresponds to the
Haar PX-DA sandwich chain outlined in Section 3. Then
$\lambda_{i,SA}=\lambda_{i}$ for all $i\geq 1$ if and only if

$f_{X\mid Y}(x\mid y)=f_{X\mid Y}(x\mid gy)\;\forall g\in
G,x\in\mathbb{R}^{p},y\in\mathbb{R}^{q}.$

If the condition in (4.4) is satisfied, then it can be shown that the Mtds $k$
and $k_{SA}$ are exactly the same. Hence, outside of this triviality, the Haar
PX-DA sandwich chain is strictly better than the DA chain in the sense that
its eigenvalues are uniformly (pointwise) dominated by those of the DA chain
with at least one strict domination.

##  5 The “two-block” DA algorithm and its sandwich variants

The applicability of the DA algorithm depends on the ease of sampling from the
conditional densities $f_{Y\mid X}$ and $f_{X\mid Y}$. While samples from both
of these densities can be generated in a straightforward way in various
settings, there are many others where this is not true. In particular, in
these settings the density $f_{X\mid Y}$ is often intractable and cannot be
directly sampled from. However, it is often possible to partition $X$ into two
components $X=(U,V)$ (with $U\in\mathbb{R}^{u}$ and $V\in\mathbb{R}^{p-u}$)
such that samples from $f_{U\mid V,Y}$ and $f_{V\mid U,Y}$ can be easily
generated. In such settings, a modified two-block DA algorithm is used to
generate samples from $f_{X}=f_{U,V}$. Each iteration of the DA algorithm
consists of three steps — a draw from $f_{Y|X}$ followed by a draw from
$f_{U|V,Y}$, and finally a draw from $f_{V|U,Y}$. The transition of this two-
block DA Markov chain from the current state $x=(u,v)$ to the next state
$x^{\prime}=(u^{\prime},v^{\prime})$ can be described as follows.

Algorithm One-step transition of the two-block DA Markov chain

1: Given $x=(u,v)$, generate a draw from $f_{Y|X}(\cdot|x)$, and call the
observed value $y$.

2: Generate a draw from $f_{U|V,Y}(\cdot|v,y)$, and call the observed value
$u^{\prime}$.

3: Generate a draw from $f_{V|U,Y}(\cdot|u^{\prime},y)$, and call the observed
value $v^{\prime}$.

From the three steps in each iteration of the two-block DA algorithm, it
follows that the Markov transition density (Mtd) of the DA Markov chain
$\\{\tilde{X}_{n}\\}_{n\geq 1}$ (with $\tilde{X}_{n}=(U_{n},V_{n})$) is given
by

$k_{TB}(x^{\prime}|x)=\int_{\mathbb{R}^{q}}f_{V|U,Y}(v^{\prime}|u^{\prime},y)f_{U|V,Y}(u^{\prime}|v,y)f_{Y|X}(y|x)\,dy\;.$
(5.1)

Its easy to show that the two-block DA chain has $f_{X}$ as its stationary
density. However, a key difference as compared to the “single-block” Markov
chain $\\{X_{n}\\}_{n\geq 0}$ in Section 1 is that the two-block DA chain does
not satisfy the detailed balance condition.

Example (Bayesian quantile regression) Consider the linear model

$Z_{i}={\bf w}_{i}^{T}{\boldsymbol{\beta}}+\epsilon_{i};i=1,2,...,n$

subject to the constraint that the $\alpha^{th}$ quantile of the error
distribution of is $0$. Recall that $\\{{\bf w}_{i}\\}_{i=1}^{n}$ is the
collection of $p$\- dimensional covariate vectors, and
${\boldsymbol{\beta}}\in\mathbb{R}^{p}$ are the regression coefficients. The
standard (frequentist) method for estimating the regression parameter
${\boldsymbol{\beta}}$ (see Yu and Moyeed, (2001)) minimizes an objective
function which is equivalent to the negative log-likelihood for
${\boldsymbol{\beta}}$ if the errors are assumed to be i.i.d. with the
asymmetric Laplace density given by

$g_{\alpha}(\epsilon)=\alpha(1-\alpha)e^{-\rho_{\alpha}(\epsilon)}.$

where $\rho_{\alpha}(x)=x\left(\alpha-I(x<0)\right)$ with $I(\cdot)$ denoting
the standard indicator function. Note that $g_{\alpha}$ has $\alpha^{th}$
quantile equal to zero, and corresponds to the standard Laplace density with
location and scale equal to $0$ and $1/2$, respectively, when $\alpha=1/2$.
Leveraging this analogy, Yu and Moyeed, (2001); Kozumi and Kobayashi, (2011)
pursue a Bayesian approach where $\\{\epsilon_{i}\\}_{i=1}^{n}$ are assumed to
be i.i.d. with common density $g_{\alpha}$, with $\sigma>0$ is an unknown
scale parameter. The following independent priors are assigned to
${\boldsymbol{\beta}}$ and $\sigma$:
${\boldsymbol{\beta}}\sim\mathcal{N}_{p}(\beta_{0},B_{0})$ and $\sigma\sim
IG\left(\frac{n_{0}}{2},\frac{t_{0}}{2}\right)$. The posterior density of
$({\boldsymbol{\beta}},\sigma)$ is intractable, and Kozumi and Kobayashi,
(2011) propose the following DA approach which exploits a normal/exponential
mixture representation of the asymmetric Laplace distribution.

Define $\theta:=\theta(\alpha)=\frac{1-2\alpha}{\alpha(1-\alpha)}$ and
$\tau^{2}:=\tau^{2}(\alpha)=\frac{2}{\alpha(1-\alpha)}$, and consider random
pairs $\left\\{(Z_{i},R_{i})\right\\}_{i=1}^{n}$ such that $Z_{i}\mid
R_{i}=r_{i},{\boldsymbol{\beta}},\sigma\sim\mathcal{N}({\bf
w}_{i}^{T}{\boldsymbol{\beta}}+\theta r_{i},r_{i}\sigma\tau^{2})$ and
$R_{i}\mid{\boldsymbol{\beta}},\sigma\sim Exp(\sigma)$. It can be easily
verified that the marginal density of $\frac{Z_{i}-{\bf
w}_{i}^{T}{\boldsymbol{\beta}}}{\sigma}$ given $({\boldsymbol{\beta}},\sigma)$
is indeed $g_{\alpha}$. However, direct sampling or closed form computations
for the joint posterior density of ${\boldsymbol{\beta}},{\bf R},\sigma$ given
${\bf Z}={\bf z}$ are not feasible. However, it can be shown that (a) the full
posterior conditional distribution of ${\boldsymbol{\beta}}$ is multivariate
Gaussian, (b) the elements of ${\bf R}$ are conditionally independent given
${\boldsymbol{\beta}},\sigma,{\bf Z}={\bf z}$, and follow a generalized
inverse Gaussian distribution, (c) the full posterior conditional distribution
of $\sigma$ is inverse gamma.

A two-block DA algorithm with ${\boldsymbol{\beta}}$ as “$U$”, ${\bf R}$ as
“$V$” and $\sigma$ as “$Y$” can hence be used to generate approximate samples
from the joint posterior density of $({\boldsymbol{\beta}},{\bf Z})$. This in
particular yields samples from the marginal posterior density of
${\boldsymbol{\beta}}$. Since the posterior density of $\sigma$ given only
${\boldsymbol{\beta}}$ can be shown to be an inverse gamma density, this
provides a mechanism to sample from the target posterior density of
$(\beta,\sigma)$.

Coming back to the general setting, a direct application/insertion of the
single-block Haar PX-DA step (Step 2) as a sandwich step in the two-block
setting is not feasible, as the resulting Markov chain does not in general
have $f_{X}$ as its stationary density. Pal et al., (2015) develop feasible
two-block adaptations of the single-block Haar PX-DA sandwich algorithm in
Section 3, which we briefly describe below.

Consider a group $G$ acting topologically on the left of $\mathbb{R}^{q}$, as
in the discussion prior to the single-block Haar PX-DA algorithm. Given any
$x\in\mathbb{R}^{p},y\in\mathbb{R}^{q}$, define the densities

$f_{x,y}^{1}(g)=\frac{f_{V,Y}(v,gy)\chi(g)}{C_{1}(x,y)}\;\forall g\in G$

and

$f_{x,y}^{2}(g)=\frac{f_{U,V,Y}(u,v,gy)\chi(g)}{C_{2}(x,y)}\;\forall g\in G.$

Here $C_{1}(x,y)$ and $C_{2}(x,y)$ (assumed to be finite) are normalizing
constants. Either of these two densities can be used to construct a ’valid’
Haar PX-DA sandwich step in the current setting. The transition of this two-
block sandwich Markov chain from the current state $x=(u,v)$ to the next state
$x^{\prime}=(u^{\prime},v^{\prime})$ can be described as follows.

Algorithm One-step transition of the two-block Haar PX-DA sandwich Markov
chain

1: Given $x=(u,v)$, generate a draw from $f_{Y|X}(\cdot|x)$, and call the
observed value $y$.

2: Generate a draw from $f_{x,y}^{j}$ (for $j=1$ or $j=2$) and call the
observed value $g$.

3: Generate a draw from $f_{U|V,Y}(\cdot|v,gy)$, and call the observed value
$u^{\prime}$.

4: Generate a draw from $f_{V|U,Y}(\cdot|u^{\prime},gy)$, and call the
observed value $v^{\prime}$.

Let $k_{TB,SA1}$ and $k_{TB,SA2}$ denote the Mtds corresponding to the two-
block sandwich chain described in two-block Haar PX-DA algorithm (with $j=1$
and $j=2$ respectively). Pal et al., (2015) show that both $k_{TB,SA1}$ and
$k_{TB,SA2}$ have $f_{X}$ as their stationary distribution.

Returning to the Bayesian quantile regression example, recall that the the
scale parameter $\sigma$ plays the role of ”$Y$”. Consider the action of the
group $G=\mathbb{R}_{+}$ on $\mathbb{R}_{+}$, the sample space of $\sigma$,
through scalar multiplication. The left Haar measure for $G$ is given by
$\nu_{l}(dg)=\frac{dg}{g}$, and it can be shown that the Lebesgue measure on
$\mathbb{R}_{+}$ is relatively invariant with respect to the multiplier
$\chi(g)=g$. Using the above construction of two-block sandwich densities, it
can be shown that $f^{1}_{x,y}$ is a non-standard univariate density. However,
samples can be generated from this density using a rejection sampler with a
dominating inverse gamma density. On the other hand, $f^{2}_{x,y}$ can be
shown to be an inverse gamma density, and the corresponding sandwich chain is
more suitable for practical use.

We conclude with a quick summary of associated theoretical results for the
two-block DA setting. Using results in Asmussen and Glynn, (2011), it follows
that the two-block DA chain and the sandwich chain are Harris ergodic if the
corresponding Mtd is strictly positive everywhere (this is the case with most
chains encountered in statistical applications). However, things get much more
complicated (compared to the single-block setting) for a comparative study of
deeper issues such as operator norms, geometric ergodicity etc. A key reason
for these challenges is that the two-block DA operator, based on the three-
stage transition step, is a product of three relevant projection operators,
and consequently loses its self-adjointness and positivity. A theoretical
comparison of closely related “lifted” versions of the two-block DA and
sandwich chains is available in Pal et al., (2015), however several questions
involving a direct comparison of the two-block DA and sandwich chains remain
open.

  * •

Under what conditions does geometric ergodicity of the two block DA chain
imply geometric ergodicity of the corresponding sandwich chain?

  * •

Although reversibility is lost, one could still aim to establish that the two-
block DA operator is Hilbert-Schmidt by showing that the condition in (4.3) is
satisfied. In such a setting, under what conditions does the Hilbert-Schmidt
property of the two-block DA chain imply the same for the corresponding
sandwich chain?

  * •

How do answers to the above questions depend on the choice of $j=1$ vs. $j=2$
in the sandwich step in the two-block Haar PX-DA algorithm?

##  6 Distributed DA adaptations for massive data settings

The sandwich DA algorithm can serve as a useful and computationally
inexpensive tool to speed up the convergence of the DA algorithm. However, its
utility is limited in increasingly common “massive” datasets, which contain
hundreds of thousands (if not millions) of samples and/or variables. For a
Bayesian statistical analysis of such datasets, a key goal is to generate
samples from the intractable posterior density $f_{X}$ of the parameters
“$X$”. For instances where the DA algorithm is applicable, the dimensionality
of “$Y$” (the latent variables) is often equal to the number of samples or the
number of variables. For example, in the Bayesian variable selection example
described in Section 2.1, the dimensionality of $Y$ is $p$, the number of
variables. The corresponding methodology is meant to be used for high-
dimensional settings where $p$ is very large. In such settings, the
computational cost of sampling the latent variables in each iteration of the
DA and sandwich chains can be exorbitant. Alternatively, consider the Bayesian
logistic regression example in Section 2.2, where the number of latent
variables is equal to the number of samples $m$. For datasets such as the
MovieLens Data Perry, (2016); Srivastava et al., (2019), where the number of
samples is in the millions, the computational cost of sampling the latent
variables can again be prohibitive.

For settings where a massive number of samples is an issue, the following
divide-and-conquer strategy can be considered: (a) Divide the data into a
number of smaller subsets of reasonable size (b) run DA Markov chains on all
subsets in parallel, and (c) suitably combine the different subset-wise Markov
chains for inference. See Scott et al., (2016); Jordan et al., (2019); Wang
and Srivastava, (2021) and the references therein. The combined draws from the
subset-wise Markov chains do not constitute a Markov chain, and quantification
of the Monte Carlo error becomes challenging. While asymptotic statistical
guarantees based on a normal approximation of the target posterior are
available for some of these methods (as the subset sample size tends to
infinity), no rigorous bounds for the distance between the distribution of the
combined parameter draws and the target posterior distribution (based on the
entire data) are available. In other words, the lack of a Markov chain
structure for the combined draws deprives the user of the rich diagnostic and
theoretical tools that are available in the literature for Markov chains.

To address the above issues, the authors in Zhou et al., (2023) develop an
asynchronous and distributed version of the DA algorithm, called ADDA, which
outputs a Markov chain with the desired posterior as its stationary
distribution. Two key assumptions that are needed for ADDA to be useful are -
(a) the latent variables can be partitioned into blocks which are
conditionally independent (given the original parameters “$X$” and the data),
and (b) the conditional posterior distribution of each latent sub-block
exclusively depends on a relevant subset of the observed data. This is for
example true for the Bayesian variable selection example in Section 2.1, where
the augmented variables $\\{\tilde{Y}_{j}\\}_{j=1}^{p}$ are conditionally
independent given the regression coefficients and the data. Also, the
conditional posterior distribution of $\tilde{Y}_{j}$ does not depend on the
observed data (hence, requirement (b) above is trivially satisfied). Again,
for the Bayesian logistic regression example in Section 2.2, the augmented
variables $\\{\tilde{Y}_{i}\\}_{i=1}^{m}$ are conditionally independent given
the regression coefficients and the data, and the conditional posterior
distribution of $\tilde{Y}_{i}$ depends only on $(\ell_{i},w_{i})$ (quantities
related to the $i^{th}$ observation).

Consider the general DA setup of Section 1, where $X$ represents the parameter
of interest, and $Y$ represents the block of latent parameters. The goal is to
sample from $f_{X}$, which represents the marginal posterior density of $X$
given the observed data (the conditioning on the observed data will be left
implicit for continuity and ease of exposition). We assume that $Y$ can be
divided into $k$ sub-blocks $Y^{1},Y^{2},\cdots,Y^{k}$ which are conditionally
independent given $X$. Also, as mentioned above, the observed data can be
partitioned into $k$ subsets such that $f_{Y^{j}\mid X}$ only depends on the
$j^{th}$ data subset for $j=1,2,\cdots,k$.

The ADDA algorithm employs $1$ manager process and $k$ worker processes. At
the level of sampling, the job of the manager process is to sample new
parameter values from $f_{X\mid Y}$ when appropriate, and to send this sample
to all the workers. The job of the $j^{th}$ worker process is to sample new
values for the $j^{th}$ latent block from $f_{Y^{j}\mid X}$ when appropriate,
and to send these samples to the manager process. Hence, the job of sampling
the latent parameter $Y$ is distributed to the $k$ worker processes. Another
key feature of the ADDA algorithm is asynchronous sampling of the $k$ latent
parameter blocks. The degree of asynchrony is controlled by user-specified
parameters $r\in(0,1]$ and $\epsilon\in(0,1]$ as follows. To make its next
conditional draw of the parameter $X$, the manager process, with probability
$\epsilon$, waits to receive updated draws of the relevant latent parameter
blocks from all the workers. But with probability $1-\epsilon$, it only waits
to receive an $r$-fraction of updated draws, and proceeds to sample $X$ given
the most recent draws for all latent parameter blocks. If a worker is in the
midst of sampling its assigned latent parameter block, but receives a new
$X$-draw from the manager, it terminates its current sampling, and begins
afresh using the latest $X$-draw received from the manager. This entire
process is summarized below.

  1. (1)

At time $t=0$, the manager starts with initial values
$(\tilde{X}_{0},\tilde{Y}_{0})$ at $t=0$ and sends $X_{0}$ to the workers.

  2. (2)

For $t=0,1,\ldots,\infty$, the manager

    1. (M-a)

waits to receive only an $r$-fraction of updated $\tilde{Y}_{t+1}^{j}$s (see
below) from the workers with probability $1-\epsilon$, and with probability
$\epsilon$, waits to receive all the updated $\tilde{Y}_{(t+1)}^{j}$s from the
workers;

    2. (M-b)

creates $\tilde{Y}_{t+1}$ by replacing the relevant $\tilde{Y}^{j}_{t}$s with
the newly received $\tilde{Y}^{j}_{t+1}$;

    3. (M-c)

draws $X_{t+1}$ from $f(X\mid\tilde{Y}_{t+1})$; and

    4. (M-d)

sends $X_{t+1}$ to all the worker processes and resets $t=t+1$.

  3. (3)

For $t=0,\ldots,\infty$, the worker $j$ ($j=1,\ldots,k$)

    1. (W-a)

waits to receive $X_{t}$ from the manager process;

    2. (W-b)

draws $\tilde{Y}^{j}_{t+1}$ from $p(Y^{j}\mid\tilde{X}_{t+1})$; and

    3. (W-c)

sends $\tilde{Y}^{j}_{t+1}$ to the manager process, resets $t=t+1$, and goes
to (W-a) if $\tilde{X}_{t+1}$ is not received from the manager before the draw
is complete; otherwise, it truncates the sampling process, resets $t=t+1$, and
goes to (W-b).

It is clear that when $\epsilon=1$ or $r=1$, the ADDA algorithm is identical
to the DA algorithm. However, things get interesting when $\epsilon<1$ and
$r<1$. In this setting, at each iteration of the ADDA algorithm, only an
$r$-fraction of the latent parameter blocks are updated (with probability
$1-\epsilon$). The ADDA algorithm can then be viewed as a mix of systematic
and random subset scan Gibbs sampler. The systematic part comes from always
updating $X$ in any given iteration, and the random subset part comes from
only updating a random subset of $[kr]$ sub-blocks of $Y$ at each iteration.
Such algorithms have been previously considered in the literature. However,
what gives the ADDA a novel and interesting flavor is that unlike existing
approaches, the comparative speed of the $k$ workers for sampling the
respective $Y$ sub-blocks can depend on the current value of $X$. In other
words, the choice of the $[kr]$ sub-blocks of $Y$ that are updated at any
given iteration can depend on the current value of the parameter $X$. Given
this dependence, it turns out that the marginal $X$-process
$\\{\tilde{X}_{t}\\}_{t\geq 0}$ is not Markov, but the joint $(X,Y)$-process
$\\{\tilde{X}_{t},\tilde{Y}_{t}\\}_{t\geq 0}$ is Markov.

In fact, it can be shown that the Markov chain
$\\{\tilde{X}_{t},\tilde{Y}_{t}\\}_{t\geq 0}$ has $f_{X,Y}$ as its stationary
(invariant) density. For ease of exposition, we provide a proof of this
assertion when $Y$ and $X$ are supported on a discrete space and $k$, the
number of latent variable blocks, is equal to $2$. We assume that $r=0.5$ and
$\epsilon=0$. The arguments below can be extended in a straightforward way to
a general setting with more than two latent variable blocks, and to non-
discrete settings with arbitrary $r$ and $\epsilon$ values in $[0,1]$.

Let $(Y_{0},X_{0})$ denote the starting value for the ADDA chain, and assume
that it is drawn from the desired posterior $f_{Y,X}$. Let $(Y_{1},X_{1})$
denote the next iterate generated by the ADDA chain. Our goal is to show that
$(Y_{1},X_{1})\sim f_{Y,X}$. As we will see below, a key assumption which
enables this is the conditional independence assumption which implies
$f_{Y\mid X}(\tilde{y}\mid\tilde{x})=f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid\tilde{x})f_{Y^{2}\mid X}(\tilde{y}_{2}\mid\tilde{x})$
(recall that $\tilde{y}_{1}$ and $\tilde{y}_{2}$ denote the two blocks of
$\tilde{y}$).

Since $X_{1}$ given $Y_{1}$ is a draw from $f_{X\mid Y}$, it follows that

$\displaystyle P((Y_{1},X_{1})=(\tilde{y},\tilde{x}))$ $\displaystyle=$
$\displaystyle P(X_{1}=\tilde{x}\mid Y_{1}=\tilde{y})P(Y_{1}=\tilde{y})$
$\displaystyle=$ $\displaystyle f_{X\mid
Y}(\tilde{x}\mid\tilde{y})P(Y_{1}=\tilde{y}).$

Hence, to prove the desired result, it is enough to show that
$P(Y_{1}=\tilde{y})=f_{Y}(\tilde{y})$. Note that

$P(Y_{1}=\tilde{y})=\sum_{y^{\prime},x^{\prime}}P(Y_{1}=\tilde{y}\mid(Y_{0},X_{0})=(y^{\prime},x^{\prime}))f_{Y,X}(y^{\prime},x^{\prime}).$

Let us recall how $Y_{1}$ is sampled given
$(Y_{0},X_{0})=(y^{\prime},x^{\prime})$. With probability say
$c_{1}(x^{\prime})$, only the first latent variable block $Y^{1}_{1}$ is
obtained using a draw from $f_{Y^{1}\mid X}(\cdot\mid x^{\prime})$ and the
second block is left unchanged at $y_{2}^{\prime}$, and with probability
$c_{2}(x^{\prime})=1-c_{1}(x^{\prime})$, only the second latent variable block
$Y^{2}_{1}$ is obtained using $f_{Y^{2}\mid X}(\cdot\mid x^{\prime})$ and the
first block is left unchanged at $y_{1}^{\prime}$. It follows that

$\displaystyle P(Y_{1}=\tilde{y})$ $\displaystyle=$
$\displaystyle\sum_{y^{\prime},x^{\prime}}c_{1}(x^{\prime})f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid
x^{\prime})1_{\\{y_{2}^{\prime}=\tilde{y}_{2}\\}}f_{Y,X}(y^{\prime},x^{\prime})+\sum_{y^{\prime},x^{\prime}}c_{2}(x^{\prime})f_{Y^{2}\mid
X}(\tilde{y}_{2}\mid
x^{\prime})1_{\\{y_{1}^{\prime}=\tilde{y}_{1}\\}}f_{Y,X}(y^{\prime},x^{\prime})$

Using $f_{Y,X}(y^{\prime},x^{\prime})=f_{Y^{1}\mid X}(y_{1}^{\prime}\mid
x^{\prime})f_{Y^{2}\mid X}(y_{2}^{\prime}\mid x^{\prime})f_{X}(x^{\prime})$
(by conditional independence of $Y^{1}$ and $Y^{2}$ given $X$), we get

$\displaystyle P(Y_{1}=\tilde{y})$ $\displaystyle=$
$\displaystyle\sum_{x^{\prime}}\sum_{y_{1}^{\prime}}\sum_{y_{2}^{\prime}}c_{1}(x^{\prime})f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid
x^{\prime})1_{\\{y_{2}^{\prime}=\tilde{y}_{2}\\}}f_{Y^{1}\mid
X}(y_{1}^{\prime}\mid x^{\prime})f_{Y^{2}\mid X}(y_{2}^{\prime}\mid
x^{\prime})f_{X}(x^{\prime})+$
$\displaystyle\sum_{x^{\prime}}\sum_{y_{1}^{\prime}}\sum_{y_{2}^{\prime}}c_{2}(x^{\prime})f_{Y^{2}\mid
X}(\tilde{y}_{2}\mid
x^{\prime})1_{\\{y_{1}^{\prime}=\tilde{y}_{1}\\}}f_{Y^{1}\mid
X}(y_{1}^{\prime}\mid x^{\prime})f_{Y^{2}\mid X}(y_{2}^{\prime}\mid
x^{\prime})f_{X}(x^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{x^{\prime}}c_{1}(x^{\prime})f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid
x^{\prime})f_{X}(x^{\prime})\left(\sum_{y_{1}^{\prime}}f_{Y^{1}\mid
X}(y_{1}^{\prime}\mid
x^{\prime})\right)\left(\sum_{y_{2}^{\prime}}f_{Y^{2}\mid
X}(y_{2}^{\prime}\mid
x^{\prime})1_{\\{y_{2}^{\prime}=\tilde{y}_{2}\\}}\right)+$
$\displaystyle\sum_{x^{\prime}}c_{2}(x^{\prime})f_{Y^{2}\mid
X}(\tilde{y}_{2}\mid
x^{\prime})f_{X}(x^{\prime})\left(\sum_{y_{2}^{\prime}}f_{Y^{2}\mid
X}(y_{2}^{\prime}\mid
x^{\prime})\right)\left(\sum_{y_{1}^{\prime}}f_{Y^{1}\mid
X}(y_{1}^{\prime}\mid
x^{\prime})1_{\\{y_{1}^{\prime}=\tilde{y}_{1}\\}}\right)$ $\displaystyle=$
$\displaystyle\sum_{x^{\prime}}c_{1}(x^{\prime})f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid x^{\prime})f_{X}(x^{\prime})f_{Y^{2}\mid
X}(\tilde{y}_{2}\mid
x^{\prime})+\sum_{x^{\prime}}c_{2}(x^{\prime})f_{Y^{2}\mid
X}(\tilde{y}_{2}\mid x^{\prime})f_{X}(x^{\prime})f_{Y^{1}\mid
X}(\tilde{y}_{1}\mid x^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{x^{\prime}}(c_{1}(x^{\prime})+c_{2}(x^{\prime}))f_{Y,X}(\tilde{y},x^{\prime}).$

The last step again uses conditional independence of $Y^{1}$ and $Y^{2}$ given
$X$. Since $c_{1}(x^{\prime})+c_{2}(x^{\prime})=1$, it follows that

$P(Y_{1}=\tilde{y})=\sum_{x^{\prime}}f_{Y,X}(\tilde{y},x^{\prime})=f_{Y}(\tilde{y}).$

Note that the above argument considers the pure asynchronous ADDA
($\epsilon=0$). But the ADDA kernel with positive $\epsilon$ is a mixture of
the DA and pure asynchronous ADDA kernels, so the result immediately follows
for such settings as well. For a setting with more than two blocks and a
general value of $r\in(0,1)$, we will have $J={{K}\choose{\lceil Kr\rceil}}$
terms in the derivation above instead of two terms. The $j^{th}$ term, with
essentially the same manipulations as above, will simplify to
$c_{j}(x^{\prime})f_{Y,X}(\tilde{y},x^{\prime})$, where $c_{j}(x^{\prime})$
denotes the probability of choosing the relevant subset of latent variable
blocks. Since $\sum_{j=1}^{J}c_{j}(x^{\prime})=1$, the invariance result will
follow. The authors in Zhou et al., (2023) also establish the geometric
ergodicity of the ADDA chain in specific settings, including the Bayesian
variable selection example considered in Section 2.1.

The parameter $r$ controls degree of asynchrony. If $r$ is small, then the
computational cost of each ADDA iteration is lower, but the ADDA Markov chain
mixes at a slower pace. This tradeoff between slow mixing and lower
computational cost per iteration is key to the choice of $r$ and the
performance of the ADDA algorithm. The hope is that for a range of values of
$r$ where the joint effect of these two competing factors leads to a
significant decrease in the overall wall clock time required for the ADDA
chain (as compared to a pure distributed implementation of the original DA
chain). This is indeed seen in the extensive experimental evaluations
performed in Zhou et al., (2023), where the ADDA is shown to have a remarkable
overall computational gain with a comparatively small loss of accuracy. For
example, when analyzing the MovieLens data Perry, (2016); Srivastava et al.,
(2019) (with millions of samples) using Bayesian logistic regression (Section
2.2), the ADDA algorithm is three to five times faster with only 2% less
accuracy compared to the (pure distributed) DA algorithm after 10,000
iterations.

##  7 Summary and discussion

DA algorithms introduce appropriate latent variables to facilitate sampling
from intractable target distributions. These algorithms explore the parameter
space by iteratively updating the parameters and the latent variables, usually
drawing from some standard distributions. Indeed, DA algorithms are a popular
choice for MCMC computing due to its simplicity and broad applicability. The
simple structure of the transition density of the DA Markov chain allows
tractable detailed spectral analyses of its convergence behavior in addition
to drift and minorization based analyses. This chapter discusses the tools for
analyzing and comparing the spectrum of DA Markov chains. The applicability of
DA methods has been demonstrated by their successful implementation in various
widely used examples from statistics and machine learning.

Despite its widespread use, the DA algorithms can suffer from slow mixing. A
great deal of effort has gone into developing techniques to tweak the DA
chain’s transition density to improve its convergence speed. Among the
different approaches, parameter expansion based methods, the PX-DA strategies,
have been the most successful. As described in this chapter, several
theoretical results have been established showing the superiority of the PX-DA
algorithms over the DA methods. For a given reparameterization, the results in
the literature identify the ‘best’ PX-DA algorithm as the Haar PX-DA based on
the Haar measure. On the other hand, there is not much study available on
comparison among different parameter expansion strategies. One exception is
Roy, (2014), which through an empirical study involving a Bayesian robit
model, shows that a partially reparameterized Haar PX-DA algorithm can
outperform a fully reparameterized Haar PX-DA algorithm. Using the same robit
example, Roy, (2014) shows that some reparameterization may provide little
improvement over the DA algorithm. An important area of future research is to
provide theoretical and methodological results for constructing and comparing
different reparameterization strategies for improving DA algorithms.

In various fields and applications, it is common to observe large data sets
that need to be analyzed. For these data sets, the number of observations or
variables is large. The dimension of the latent vector $y$ in a DA scheme is
usually the same as the number of observations or the number of variables. On
the other hand, generally, as seen in all examples considered in this chapter,
the latent variables are conditionally independent given the observed data and
the current parameter value. Thus, the draw from $y$, the first of the two
steps in every iteration of the DA algorithm, can be parallelized,
distributing the computational burden across multiple processors or nodes.
Indeed, conditional independence of the latent blocks
$Y^{1},Y^{2},\cdots,Y^{k}$ is key for ensuring that the ADDA-process
$\\{(\tilde{X}_{t},\tilde{Y}_{t})\\}_{t=0}^{\infty}$ described in Section 6 is
a Markov chain with $f_{X,Y}$ as its stationary distribution. However, such
conditional independence is lost in many applications, such as hidden Markov
models and time series models. An important open research direction is the
extension of ADDA-type ideas that adapt to this loss of conditional
independence without compromising on some underlying Markov structure. Another
open direction is a useful integration of the sandwich technique to develop
parameter-expanded versions of the ADDA algorithm.

## Acknowledgments

The first author’s work was partially supported by USDA NIFA Grant
2023-70412-41087.

## References

  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  Journal of the American Statistical Association, 88:669–679.
  * Andrews and Mallows, (1974) Andrews, D. F. and Mallows, C. F. (1974).  Scale mixtures of normal distributions.  Journal of the Royal Statistical Society, Series B, 36:99–102.
  * Asmussen and Glynn, (2011) Asmussen, S. and Glynn, P. W. (2011).  A new proof of convergence of MCMC via the ergodic theorem.  Statistics and Probability Letters, 81:1482–1485.
  * Booth and Hobert, (1999) Booth, J. G. and Hobert, J. P. (1999).  Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm.  Journal of the Royal Statistical Society, Series B, 61:265–285.
  * Buja, (1990) Buja, A. (1990).  Remarks on functional canonical variates, alternating least squares methods and ACE.  The Annals of Statistics, 18:1032–1069.
  * Chakraborty and Khare, (2017) Chakraborty, S. and Khare, K. (2017).  Convergence properties of Gibbs samplers for Bayesian probit regression with proper priors.  Electronic Journal of Statistics, 11(1):177–210.
  * Chakraborty and Khare, (2019) Chakraborty, S. and Khare, K. (2019).  Consistent estimation of the spectrum of trace class data augmentation algorithms.  Bernoulli, 25:3832–3863.
  * Chan and Geyer, (1994) Chan, K. S. and Geyer, C. J. (1994).  Comment on “Markov chains for exploring posterior distributions” by L. Tierney.  The Annals of Statistics, 22:1747–1758.
  * Choi and Hobert, (2013) Choi, H. M. and Hobert, J. P. (2013).  The Polya-Gamma Gibbs sampler for Bayesian logistic regression is uniformly ergodic.  Electronic Journal of Statistics, 7:2054–2064.
  * Davies, (1983) Davies, E. (1983).  LINEAR INTEGRAL OPERATORS: Surveys and Reference Works in Mathematics, 7.  Wiley Online Library.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the EM algorithm (with discussion).  Journal of the Royal Statistical Society, Series B, 39:1–38.
  * Diaconis et al., (2008) Diaconis, P., Khare, K., and Saloff-Coste, L. (2008).  Gibbs sampling, exponential families and orthogonal polynomials (with discussion).  Statistical Science, 23:151–200.
  * Diebolt and Robert, (1994) Diebolt, J. and Robert, C. P. (1994).  Estimation of finite mixture distributions by Bayesian sampling.  Journal of the Royal Statistical Society, Series B, 56:363–375.
  * Duan et al., (2018) Duan, L. L., Johndrow, J. E., and Dunson, D. B. (2018).  Scaling up data augmentation MCMC via calibration.  The Journal of Machine Learning Research, 19(1):2575–2608.
  * Eaton, (1989) Eaton, M. L. (1989).  Group Invariance Applications in Statistics.  Institute of Mathematical Statistics and the American Statistical Association, Hayward, California and Alexandria, Virginia.
  * Flegal and Jones, (2010) Flegal, J. M. and Jones, G. L. (2010).  Batch means and spectral variance estimators in Markov chain Monte Carlo.  The Annals of Statistics, 38:1034–1070.
  * Frühwirth-Schnatter and Frühwirth, (2010) Frühwirth-Schnatter, S. and Frühwirth, R. (2010).  Data augmentation and MCMC for binary and multinomial logit models.  In Statistical Modelling and Regression Structures, pages 111–132. Springer.
  * Geyer, (1994) Geyer, C. J. (1994).  On the convergence of Monte Carlo maximum likelihood calculations.  Journal of the Royal Statistical Society, Series B, 56:261–274.
  * Geyer and Thompson, (1992) Geyer, C. J. and Thompson, E. A. (1992).  Constrained Monte Carlo maximum likelihood for dependent data.  Journal of the Royal Statistical Society, Series B, 54:657–699.
  * Gilks and Wild, (1992) Gilks, W. R. and Wild, P. (1992).  Adaptive rejection sampling for Gibbs sampling.  Applied Statistics, 41(2):337–348.
  * Hobert, (2011) Hobert, J. P. (2011).  Handbook of Markov chain Monte Carlo, chapter The data augmentation algorithm: theory and methodology, pages 253–293.  CRC Press, Boca Raton, FL.
  * Hobert and Marchev, (2008) Hobert, J. P. and Marchev, D. (2008).  A theoretical comparison of the data augmentation, marginal augmentation and PX-DA algorithms.  The Annals of Statistics, 36:532–554.
  * Hobert et al., (2011) Hobert, J. P., Roy, V., and Robert, C. P. (2011).  Improving the convergence properties of the data augmentation algorithm with an application to Bayesian mixture modelling.  Statistical Science, 26:332–351.
  * Holmes and Held, (2006) Holmes, C. C. and Held, L. (2006).  Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian Analysis, 1:145–168.
  * Jin and Tan, (2021) Jin, R. and Tan, A. (2021).  Fast Markov chain Monte Carlo for high-dimensional Bayesian regression models with shrinkage priors.  Journal of Computational and Graphical Statistics, 30(3):632–646.
  * Jones and Hobert, (2001) Jones, G. L. and Hobert, J. P. (2001).  Honest exploration of intractable probability distributions via Markov chain Monte Carlo.  Statistical Science, 16:312–34.
  * Jones and Hobert, (2004) Jones, G. L. and Hobert, J. P. (2004).  Sufficient burn-in for Gibbs samplers for a hierarchical random effects model.  The Annals of Statistics, 32:784–817.
  * Jordan et al., (2019) Jordan, M. I., Lee, J. D., and Yang, Y. (2019).  Communication-efficient distributed statistical inference.  Journal of the American Statistical Association, 114(526):668–681.
  * Khare and Hobert, (2011) Khare, K. and Hobert, J. P. (2011).  A spectral analytic comparison of trace-class data augmentation algorithms and their sandwich variants.  The Annals of Statistics, 39:2585–2606.
  * Kozumi and Kobayashi, (2011) Kozumi, H. and Kobayashi, G. (2011).  Gibbs sampling methods for Bayesian quantile regression.  Journal of statistical computation and simulation, 81(11):1565–1578.
  * Kyung et al., (2010) Kyung, M., Gill, J., Ghosh, M., and Casella, G. (2010).  Penalized regression, standard errors, and Bayesian Lassos.  Bayesian Analysis, 5:369–412.
  * Laha et al., (2016) Laha, A., Dutta, S., and Roy, V. (2016).  A novel sandwich algorithm for empirical Bayes analysis of rank data.  Statistics and its Interface, 10:543–556.
  * Liu, (2004) Liu, C. (2004).  Robit regression: A simple robust alternative to logistic and probit regression.  In Gelman, A. and Meng, X. L., editors, Applied Bayesian Modeling and Casual Inference from Incomplete-Data Perspectives, pages 227–238. Wiley, London.
  * Liu, (1996) Liu, J. S. (1996).  Peskun’s theorem and a modified discrete-state Gibbs sampler.  Biometrika, 83:681–682.
  * Liu et al., (1994) Liu, J. S., Wong, W. H., and Kong, A. (1994).  Covariance structure of the Gibbs sampler with applications to comparisons of estimators and augmentation schemes.  Biometrika, 81:27–40.
  * Liu et al., (1995) Liu, J. S., Wong, W. H., and Kong, A. (1995).  Covariance structure and convergence rate of the Gibbs sampler with various scans.  Journal of the Royal Statistical Society, Series B, 57:157–169.
  * Liu and Wu, (1999) Liu, J. S. and Wu, Y. N. (1999).  Parameter expansion for data augmentation.  Journal of the American Statistical Association, 94:1264–1274.
  * Meng and van Dyk, (1999) Meng, X.-L. and van Dyk, D. A. (1999).  Seeking efficient data augmentation schemes via conditional and marginal augmentation.  Biometrika, 86:301–320.
  * Meyn and Tweedie, (1993) Meyn, S. P. and Tweedie, R. L. (1993).  Markov Chains and Stochastic Stability.  Springer-Verlag, London.
  * Mukherjee et al., (2023) Mukherjee, S., Khare, K., and Chakraborty, S. (2023).  Convergence properties of data augmentation algorithms for high-dimensional robit regression.  Electronic Journal of Statistics, 17(1):19–69.
  * Pal et al., (2015) Pal, S., Khare, K., and Hobert, J. P. (2015).  Improving the data augmentation algorithm in the two-block setup.  Journal of Computational and Graphical Statistics, 24(4):1114–1133.
  * Pal et al., (2017) Pal, S., Khare, K., and Hobert, J. P. (2017).  Trace class Markov chains for Bayesian inference with generalized double Pareto shrinkage priors.  Scandinavian Journal of Statistics, 44(2):307–323.
  * Papaspiliopoulos and Roberts, (2008) Papaspiliopoulos, O. and Roberts, G. (2008).  Stability of the Gibbs sampler for Bayesian hierarchical models.  The Annals of Statistics, 38:95–117.
  * Park and Casella, (2008) Park, T. and Casella, G. (2008).  The Bayesian Lasso.  Journal of the American Statistical Association, 103:681–686.
  * Perry, (2016) Perry, P. O. (2016).  Fast moment-based estimation for hierarchical models.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(1):267–291.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya-Gamma latent variables.  Journal of the American statistical Association, 108(504):1339–1349.
  * Pregibon, (1982) Pregibon, D. (1982).  Resistant fits for some commonly used logistic models with medical applications.  Biometrics, 38:485–498.
  * Qin and Hobert, (2021) Qin, Q. and Hobert, J. P. (2021).  On the limitations of single-step drift and minorization in Markov chain convergence analysis.  The Annals of Applied Probability, 31(4):1633–1659.
  * Qin and Hobert, (2022) Qin, Q. and Hobert, J. P. (2022).  Geometric convergence bounds for Markov chains in Wasserstein distance based on generalized drift and contraction conditions.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 58(2):872–889.
  * Qin et al., (2019) Qin, Q., Hobert, J. P., and Khare, K. (2019).  Estimating the spectral gap of a trace-class Markov operator.  Electronic Journal of Statistics, 13:1790 – 1822.
  * Rajaratnam et al., (2019) Rajaratnam, B., Sparks, D., Khare, K., and Zhang, L. (2019).  Scalable Bayesian shrinkage and uncertainty quantification in high-dimensional regression.  Journal of Computational and Graphical Statistics, 28(1):174–184.
  * Rao and Roy, (2021) Rao, Y. and Roy, V. (2021).  Block Gibbs samplers for logistic mixed models: convergence properties and a comparison with full Gibbs samplers.  Electronic Journal of Statistics, 15:5598–5625.
  * Roberts and Rosenthal, (1997) Roberts, G. O. and Rosenthal, J. S. (1997).  Geometric ergodicity and hybrid Markov chains.  Electronic Communications in Probability, 2:13–25.
  * Rosenthal, (1995) Rosenthal, J. S. (1995).  Minorization conditions and convergence rates for Markov chain Monte Carlo.  Journal of the American Statistical Association, 90:558–566.
  * (55) Roy, V. (2012a).  Convergence rates for MCMC algorithms for a robust Bayesian binary regression model.  Electronic Journal of Statistics, 6:2463–2485.
  * (56) Roy, V. (2012b).  Spectral analytic comparisons for data augmentation.  Statistics & Probability Letters, 82:103–108.
  * Roy, (2014) Roy, V. (2014).  Efficient estimation of the link function parameter in a robust Bayesian binary regression model.  Computational Statistics and Data Analysis, 73:87–102.
  * Roy, (2016) Roy, V. (2016).  Improving efficiency of data augmentation algorithms using Peskun’s theorem.  Computational Statistics, 31:709–728.
  * Roy, (2022) Roy, V. (2022).  MCMC for GLMMs.  In C.R. Rao, A. S. R. and Young, A., editors, Handbook of Statistics Volume 47: Advancements in Bayesian Methods and Implementation, pages 135–159. Elsevier.
  * Roy and Chakraborty, (2017) Roy, V. and Chakraborty, S. (2017).  Selection of tuning parameters, solution paths and standard errors for Bayesian lassos.  Bayesian Analysis, 12:753–778.
  * Roy and Hobert, (2007) Roy, V. and Hobert, J. P. (2007).  Convergence rates and asymptotic standard errors for Markov chain Monte Carlo algorithms for Bayesian probit regression.  Journal of the Royal Statistical Society, Series B, 69:607–623.
  * Scott et al., (2016) Scott, S. L., Blocker, A. W., Bonassi, F. V., Chipman, H. A., I., G. E., and McCulloch, R. E. (2016).  Bayes and big data: the consensus Monte Carlo algorithm.  International Journal of Management Science and Engineering Management, 11(2):78–88.
  * Srivastava et al., (2019) Srivastava, S., DePalma, G., and Liu, C. (2019).  An asynchronous distributed expectation maximization algorithm for massive data: The dem algorithm.  Journal of Computational and Graphical Statistics, 28(2):233–243.
  * Swendsen and Wang, (1987) Swendsen, R. H. and Wang, J.-S. (1987).  Nonuniversal critical dynamics in Monte Carlo simulations.  Physical Review Letters, 58:86–88.
  * Tanner and Wong, (1987) Tanner, M. A. and Wong, W. H. (1987).  The calculation of posterior distributions by data augmentation (with discussion).  Journal of the American Statistical Association, 82:528–550.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society, Series B, 58:267–288.
  * Vats et al., (2019) Vats, D., Flegal, J. M., and Jones, G. L. (2019).  Multivariate output analysis for Markov chain Monte Carlo.  Biometrka, 106:321–337.
  * Wang and Srivastava, (2021) Wang, C. and Srivastava, S. (2021).  Divide-and-conquer Bayesian inference in hidden Markov models.  Electronic Journal of Statistics, 17:895–947.
  * (69) Wang, X. and Roy, V. (2018a).  Analysis of the Pólya-Gamma block Gibbs sampler for Bayesian logistic linear mixed models.  Statistics & Probability Letters, 137:251–256.
  * (70) Wang, X. and Roy, V. (2018b).  Convergence analysis of the block Gibbs sampler for Bayesian probit linear mixed models with improper priors.  Electronic Journal of Statistics, 12:4412–4439.
  * (71) Wang, X. and Roy, V. (2018c).  Geometric ergodicity of Pólya-Gamma Gibbs sampler for Bayesian logistic regression with a flat prior.  Electronic Journal of Statistics, 12:3295–3311.
  * Xu and Ghosh, (2015) Xu, X. and Ghosh, M. (2015).  Bayesian variable selection and estimation for group lasso.  Bayesian Analysis, 10(4):909–936.
  * Yu and Moyeed, (2001) Yu, K. and Moyeed, R. A. (2001).  Bayesian quantile regression.  Statistics & Probability Letters, 54(4):437–447.
  * Yu and Meng, (2011) Yu, Y. and Meng, X.-L. (2011).  To center or not to center: that is not the question - An ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC efficiency (with discussion).  Journal of Computational and Graphical Statistics, 20:531–570.
  * Zhou et al., (2023) Zhou, J., Khare, K., and Srivastava, S. (2023).  Asynchronous and distributed data augmentation for massive data settings.  Journal of Computational and Graphical Statistics, 32(3):895–907.
  * Zou and Hastie, (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the Elastic Net.  Journal of the Royal Statistical Society, Series B, 67:301–320.