# Variable-Complexity Weighted-Tempered Gibbs Samplers for Bayesian Variable
Selection

Lan V. Truong<EMAIL_ADDRESS>Department of Engineering, University of
Cambridge, Trumpington, Cambridge, CB2 1PZ, Cambridgeshire, United Kingdom

###### Abstract

Subset weighted-Tempered Gibbs Sampler (wTGS) has been recently introduced by
Jankowiak to reduce the computation complexity per MCMC iteration in high-
dimensional applications where the exact calculation of the posterior
inclusion probabilities (PIP) is not essential. However, the Rao-Backwellized
estimator associated with this sampler has a high variance as the ratio
between the signal dimension and the number of conditional PIP estimations is
large. In this paper, we design a new subset weighted-Tempered Gibbs Sampler
(wTGS) where the expected number of computations of conditional PIPs per MCMC
iteration can be much smaller than the signal dimension. Different from the
subset wTGS and wTGS, our sampler has a variable complexity per MCMC
iteration. We provide an upper bound on the variance of an associated Rao-
Blackwellized estimator for this sampler at a finite number of iterations,
$T$, and show that the variance is
$O\big{(}\big{(}\frac{P}{S}\big{)}^{2}\frac{\log T}{T}\big{)}$ for a given
dataset where $S$ is the expected number of conditional PIP computations per
MCMC iteration. Experiments show that our Rao-Blackwellized estimator can have
a smaller variance than its counterpart associated with the subset wTGS.

###### keywords:

Bayesian Variable Selection, Tempered Gibbs Sampler, Computational Complexity,
Linear Regressions, Linear Models.

##  1 Introduction

Markov chain Monte Carlo (MCMC) methods comprise a class of algorithms for
sampling from a known function. MCMC methods are primarily used for
calculating numerical approximations of multi-dimensional integrals, for
example in Bayesian statistics, computational physics [12], computational
biology, [9], and linear models [18, 19, 20]. Monte Carlo algorithms have been
very popular over the last decade [10, 17]. Many practical problems in
statistical signal processing, machine learning and statistics, demand fast
and accurate procedures for drawing samples from probability distributions
that exhibit arbitrary, non-standard forms [1, 8, 16]. One of the most popular
Monte Carlo methods are the families of Markov chain Monte Carlo (MCMC)
algorithms [1, 17] and particle filters [5]. The MCMC techniques generate a
Markov chain with a pre-established target probability density function as
invariant density [13].

Gibbs sampler (GS) is a Markov chain Monte Carlo (MCMC) algorithm for
obtaining a sequence of observations which are approximated from a specific
multivariate probability distribution. This sequence can be used to
approximate the joint distribution, the marginal distribution of one of the
variables, or some subset of the variables. It can be also used to compute the
expected value (integral) of one of the variables [2, 3]. GS is applicable
when the joint distribution is not known explicitly or is difficult to sample
from directly, but the conditional distribution of each variable is known and
is easy (or at least, easier) to sample from. The GS algorithm generates an
instance from the distribution of each variable in turn, conditional on the
current values of the other variables. It can be shown that the sequence of
samples constitutes a Markov chain, and the stationary distribution of that
Markov chain is just the sought-after joint distribution.

GS is commonly used as a means of statistical inference, especially Bayesian
inference. However, pure Markov chain based schemes (i.e., ones which simulate
from precisely the right target distribution with no need for subsequent
important sampling correction) have been far more successful. This is because
MCMC methods are usually much more scalable to high-dimensional situations,
whereas important sampling weight variances tend to grow (often exponentially)
with dimension. [23] proposed a natural way to combine the best of MCMC and
importance sampling in a way that is robust in high-dimensional contexts and
ameliorates the slow mixing which plagues many Markov chain based schemes. The
proposed scheme is called Tempered Gibbs Sampler (TGS), involving component-
wise updating rule like Gibbs Sampling (GS), with improved mixing properties
and associated importance weights which remain stable as dimension increases.
Through an appropriately designed tempering mechanism, TGS circumvents the
main limitations of standard GS, such as the slow mixing introduced by strong
posterior correlations. It also avoids the requirement to visit all
coordinates sequentially, instead iteratively making state-informed decisions
as to which coordinate should be next updated.

TGS has been applied to Bayesian Variable Selection (BVS) problem, observing
multiple orders of magnitude improvements compared to alternative Monte Carlo
schemes [23]. Since TGS updates each coordinate with the same frequency, in a
BVS context, this may be inefficient as the resulting sampler would spend most
iterations updating variables that have low or negligible posterior inclusion
probability, especially when the signal dimension, $P$, gets large. A better
solution, called weighted Tempered Gibbs Sampling (wTGS) [23], updates more
often components with a larger inclusion probability, thus having a more
focused computational effort. However, despite the intuitive appeal of this
approach to BVS problem, approximating the resulting posterior distribution
can be computationally challenging. A principal reason for this is the
astronomical size of the model space that results whenever there more than a
few dozen covariates. To scale the high-dimensional regime, [11] has recently
introduced an efficient MCMC scheme whose cost per iteration can be
significantly reduced compared to wTGS. The main idea is to introduce an
auxiliary variable that controls which conditional posterior inclusion
probabilites (PIPs) are computed in a given MCMC iteration. However, this
scheme contains some weaknesses such as the Rao-Blackwellized estimator
associated with this sampler has very high variance as $P/S$ large at a small
number of MCMC iterations, $T$. In addition, generating an auxiliary random
set which is uniformly distributed over ${P\choose S}$ subsets in the subset
wTGS algorithm [11] requires very long running time.

In this paper, we design a new subset wTGS called variable-complexity wTGS
(VC-wTGS). The Rao-Blackwellized estimator associate with VC-wTGS is shown to
have the variance $O\big{(}\frac{\log
T}{T}\big{(}\frac{P}{S}\big{)}^{2}\big{)}$, where $P,S,T$ are the signal
dimension, the number of PIP computations per an MCMC iteration, and the
number of MCMC iterations, respectively. Experiments show that our scheme
converges to PIPs very fast for simulated datasets and that the variance of
the Rao-Blackwellized estimator can be much smaller than subset wTGS [11] when
$P/S$ is very high for MNIST dataset. More specifically, our contributions
include:

  * •

We propose a new subset wTGS, called VC-wTGS, where the expected number of
conditional PIP computations per MCMC can be much smaller than the signal
dimension.

  * •

We analyse the variance of an associated Rao-Blackwellized estimator at each
finite number of MCMC iterations. We show that this variance is
$O\big{(}\frac{\log T}{T}\big{(}\frac{P}{S}\big{)}^{2}\big{)}$ for a given
dataset.

  * •

We provide some experiments on a simulated dataset (multivariate Gaussian
dataset) and the real dataset (MNIST). Experiments show that our estimator can
have better variance than the subset wTGS-based estimator [11] at high $P/S$
for the same number of MCMC iterations $T$.

##  2 Preliminaries

###  2.1 Mathematical Backgrounds

Let a Markov chain $\\{X_{n}\\}_{n=1}^{\infty}$ on a state space $\mathcal{S}$
with transition kernel $Q(x,dy)$ and the initial state $X_{1}\sim\nu$, where
$\mathcal{S}$ is a Polish space in $\mathbb{R}$. In this paper, we consider
the Markov chains which are irreducible and positive-recurrent, so the
existence of a stationary distribution $\pi$ is guaranteed. An irreducible and
recurrent Markov chain on an infinite state-space is called Harris chain [21].
A Markov chain is called _reversible_ if the following detailed balance
condition is satisfied:

$\displaystyle\pi(dx)Q(x,dy)=\pi(dy)Q(y,dx),\qquad\forall x,y\in\mathcal{S}.$
(1)

Define

$\displaystyle
d(t)=\sup_{x\in\mathcal{S}}d_{\rm{TV}}(Q^{t}(x,\cdot),\pi),\qquad
t_{\rm{mix}}(\varepsilon):=\min\\{t:d(t)\leq\varepsilon\\},$ (2)

and

$\displaystyle\tau_{\min}:=\inf_{0\leq\varepsilon\leq
1}t_{\rm{mix}}(\varepsilon)\bigg{(}\frac{2-\varepsilon}{1-\varepsilon}\bigg{)}^{2},\qquad
t_{\rm{mix}}:=t_{\rm{mix}}(1/4).$ (3)

Let $L_{2}(\pi)$ be the Hilbert space of complex valued measurable functions
on $\mathcal{S}$ that are square integrable w.r.t. $\pi$. We endow
$L_{2}(\pi)$ with inner product $\langle f,g\rangle:=\int fg^{*}d\pi$, and
norm $\|f\|_{2,\pi}:=\langle f,f\rangle_{\pi}^{1/2}$. Let $E_{\pi}$ be the
associated averaging operator defined by $(E_{\pi})(x,y)=\pi(y),\forall
x,y\in\mathcal{S}$, and

$\displaystyle\lambda=\|Q-E_{\pi}\|_{L_{2}(\pi)\to L_{2}(\pi)},$ (4)

where $\|B\|_{L_{2}(\pi)\to
L_{2}(\pi)}=\max_{v:\|v\|_{2,\pi}=1}\|Bv\|_{2,\pi}.$ $Q$ can be viewed as a
linear operator (infinitesimal generator) on $L_{2}(\pi)$, denoted by
$\mathbf{Q}$, defined as $(\mathbf{Q}f)(x):=\mathbb{E}_{Q(x,\cdot)}(f)$, and
the reversibility is equivalent to the self-adjointness of $\mathbf{Q}$. The
operator $\mathbf{Q}$ acts on measures on the left, creating a measure
$\mu\mathbf{Q}$, that is, for every measurable subset $A$ of $\mathcal{S}$,
$\mu\mathbf{Q}(A):=\int_{x\in\mathcal{S}}Q(x,A)\mu(dx)$. For a Markov chain
with stationary distribution $\pi$, we define the _spectrum_ of the chain as

$\displaystyle
S_{2}:=\big{\\{}\xi\in\mathbb{C}:(\xi\mathbf{I}-\mathbf{Q})\enspace\mbox{is
not invertible on}\enspace L_{2}(\pi)\big{\\}}.$ (5)

It is known that $\lambda=1-\gamma^{*}$ [14], where

$\displaystyle\gamma^{*}$
$\displaystyle:=\begin{cases}1-\sup\\{|\xi|:\xi\in\mathcal{S}_{2},\xi\neq
1\\},\\\ \qquad\qquad\mbox{if eigenvalue $1$ has multiplicity $1$,}\\\
0,\qquad\qquad\mbox{otherwise}\end{cases}$

is the _the absolute spectral gap_ of the Markov chain. The absolute spectral
gap can be bounded by the mixing time $t_{\rm{mix}}$ of the Markov chain by
the following expression:

$\displaystyle\bigg{(}\frac{1}{\gamma^{*}}-1\bigg{)}\log 2\leq
t_{\rm{mix}}\leq\frac{\log(4/\pi_{*})}{\gamma_{*}},$ (6)

where $\pi_{*}=\min_{x\in\mathcal{S}}\pi_{x}$ is the _minimum stationary
probability_ , which is positive if $Q^{k}>0$ (entry-wise positive) for some
$k\geq 1$. See [22] for more detailed discussions. In [6, 22], the authors
provided algorithms to estimate $t_{\rm{mix}}$ and $\gamma^{*}$ from a single
trajectory.

Define

$\displaystyle\mathcal{M}_{2}:=\bigg{\\{}\nu\in\mathcal{M}(\mathcal{S}):\bigg{\|}\frac{dv}{d\pi}\bigg{\|}_{2}<\infty\bigg{\\}},$
(7)

where $\|\cdot\|_{2}$ is the standard $L_{2}$ norm in the Hilbert space of
complex valued measurable functions on $\mathcal{S}$.

###  2.2 Problem Set-up

Consider the linear regression with $X\in\mathbb{R}^{N\times P}$ and
$Y\in\mathbb{R}^{N}$ and define the following space of models:

  * •

inclusion variables: $\gamma_{i}\sim\texttt{Bern}(h)$

  * •

noise variance:
$\sigma_{\gamma}^{2}\in\texttt{InvGamma}\big{(}\frac{1}{2}\nu_{0},\frac{1}{2}\nu_{0}\lambda_{0}\big{)}$

  * •

coefficients:
$\beta_{\gamma}\sim\mathcal{N}(0,\sigma_{\gamma}^{2}\tau^{-1}\mathbf{I}_{|\gamma|})$

  * •

response:
$Y_{n}\sim\mathcal{N}(\beta_{\gamma}.X_{n\gamma},\sigma_{\gamma}^{2})$

where $i=1,2,\cdots,P$ and $n=1,2,\cdots,N$. Here each
$\gamma_{i}\in\\{0,1\\}$ controls whether the coefficient $\beta_{i}$ and the
$i$-th covariate are included $(\gamma_{i}=1)$ or excluded $(\gamma_{i}=0)$
from the model. In the following, we use $\gamma$ to refer to the vector
$(\gamma_{1},\gamma_{2},\cdots,\gamma_{P})$. The hyperparameter $h\in(0,1)$
controls the overall level of sparsity; in particular $hP$ is the expected
number of covariates included a priori. The $|\gamma|$ coefficients
$\beta_{\gamma}\in\mathbb{R}^{|\gamma|}$ are governed by the standard Gaussian
prior with precision proportional to $\tau>0$. Here,
$|\gamma|\in\\{0,1,2,\cdots,P\\}$ denotes the total number of included
covariates. The response $Y_{n}$ is generated from a Gaussian distribution
with variance governed by an Inverse Gamma prior. Note that we do not include
a bias term, but doing so may be desirable in practice. An attractive feature
of the model is that it explicitly reasons about variable inclusion and allows
us to define _posterior inclusion probabilities_ or PIPs, where

$\displaystyle\texttt{PIP}(i):=p(\gamma_{i}=1|\mathcal{D})\in[0,1]$ (8)

and $\mathcal{D}=\\{X,Y\\}$ is the observed dataset.

##  3 Main Results

###  3.1 Introduction to Subset wTGS

In this subsection, we review the subset wTGS which was proposed by [11].
Consider the following (unnormalized) target distribution:

$\displaystyle
f(\gamma,i,\mathcal{S}):=p(\gamma|\mathcal{D})\frac{\frac{1}{2}\eta(\gamma_{-i})}{p(\gamma_{i}|\gamma_{-i},\mathcal{D})}\mathcal{U}(\mathcal{S}|i,\mathcal{A}).$
(9)

Here, $\mathcal{S}$ ranges over all the subsets of $\\{1,2,\cdots,P\\}$ of
size $S$ that also contain a fixed ‘anchor’ set
$\mathcal{A}\subset\\{1,2,\cdots,P\\}$ if size $A<S$, and $\eta(\cdot)$ is
some weighting functions. Moreover, $U(\mathcal{S}|i,\mathcal{A})$ is the
uniform distribution over the all size $S$ subsets of $\\{1,2,\cdots,P\\}$
that contain both $i$ and $\mathcal{A}$.

In practice, the set $\mathcal{A}$ can be chosen during burn-in. Subset wTGS
proceeds by defining a sampling scheme for the target distribution (9) that
utilizes Gibbs updates w.r.t. $i$ and $\mathcal{S}$ and Metropolized-Gibbs
update w.r.t. $\gamma_{i}$.

  * •

$i$-updates: Marginalizing $i$ from (9) yields

$\displaystyle
f(\gamma,\mathcal{S})=p(\gamma|\mathcal{D})\phi(\gamma,\mathcal{S})$ (10)

where we define

$\displaystyle\phi(\gamma,\mathcal{S}):=\sum_{i\in\mathcal{S}}\frac{\frac{1}{2}\eta(\gamma_{-i})}{p(\gamma_{i}|\gamma_{-i},\mathcal{D})}\mathcal{U}(\mathcal{S}|i,\mathcal{A})$
(11)

and have leveraged that $\mathcal{U}(\mathcal{S}|i,\mathcal{A})=0$ if
$i\notin\mathcal{S}$. Crucially, computing $\phi(\gamma,\mathcal{S})$ is
$\Theta(S)$ instead of $\Theta(P)$. We can do Gibbs updates w.r.t. $i$ using
the distribution

$\displaystyle
f(i|\gamma,\mathcal{S})\sim\frac{\eta(\gamma_{-i})}{p(\gamma_{i}|\gamma_{-i},\mathcal{D})}\mathcal{U}(\mathcal{S}|i,\mathcal{A}).$
(12)

  * •

$\gamma$-updates: Just as for $wTGS$ we utilized Metropolized -Gibbs updates
w.r.t. $\gamma_{i}$ that result in deterministic flips $\gamma_{i}\to
1-\gamma_{i}$. Likewise the marginal $f(i)$ is proportional to
$\texttt{PIP}(i)+\frac{\varepsilon}{P}$ so that the sampler focuses
computational efforts on large PIP covariates.

  * •

$\mathcal{S}$-updates: $\mathcal{S}$ is updated with Gibbs moves,
$\mathcal{S}\sim\mathcal{U}(\cdot|i,\mathcal{A})$. For the full algorithm, see
the Algorithm 1.

Algorithm 1 The Subset $S$-wTGS Algorithm

Input: Dataset $\mathcal{D}=\\{X,Y\\}$ with $P$ covariates; prior inclusion
probability $h$; prior precision $\tau$; subset size $S$; anchor set size $A$;
total number of MCMC iterations $T$; number of burn-in iteration
$T_{\rm{burn}}$.

Output: Approximate weighted posterior samples
$\\{\rho^{(t)},\gamma^{(t)}\\}_{t=T_{\rm{burn}}+1}^{T}$

Initializations:
$\gamma^{(0)}=\underbrace{(0,0,\cdots,0)}_{P\enspace\text{covariates}}$ and
choose $\mathcal{A}$ be the $A$ covariate with exhibiting the largest
correlations with $Y$. Choose $i^{(0)}$ randomly from $\\{1,2,\cdots,P\\}$ and
$\mathcal{S}^{(0)}\sim\mathcal{U}(\cdot|i^{(0)},\mathcal{A})$.

for $t=1,2,\cdots,T$ do

Estimate
$f(j|\gamma^{(t-1)})\leftarrow\phi_{t-1}(\gamma)^{-1}\frac{\frac{1}{2}\eta(\gamma_{-j}^{(t-1)})}{p(\gamma_{j}^{(t-1)}|\gamma_{-j}^{(t-1)},\mathcal{D})}$
for all $j\in[P]$.

Sample $i^{(t)}\sim f(\cdot|\gamma^{(t-1)})$

$\gamma^{(t)}\leftarrow\texttt{flip}(\gamma^{(t-1)}|i^{(t)})$ where
$\texttt{flip}(\gamma|i)$ flips the $i$-th coordinate of
$\gamma:\gamma_{i}\leftarrow 1-\gamma_{i}$.

Sample $\mathcal{S}^{(t)}\sim\mathcal{U}(\cdot|i^{(t)},\mathcal{A})$

Estimate $S$ conditional PIPs
$p(\gamma_{j}^{(t)}|\gamma_{-j}^{(t)},\mathcal{D})$ for all
$j\in\mathcal{S}^{(t)}$

$\phi_{t}(\gamma)\leftarrow\sum_{j\in\mathcal{S}^{(t)}}\frac{\frac{1}{2}\eta(\gamma_{-j}^{(t)})}{p(\gamma_{j}^{(t)}|\gamma_{-j}^{(t)},\mathcal{D})}$

Compute the unnormalized weights
$\tilde{\rho}^{(t)}\leftarrow\phi^{-1}(\gamma^{(t)})$

if $t\leq T_{\rm{burn}}$ then

Adapt $\mathcal{A}$ using some adaptive scheme.

end if

end for

for $t=1,2,\cdots,T$ do

$\rho^{(t)}\leftarrow\frac{\tilde{\rho}^{(t)}}{\sum_{s>T_{\rm{burn}}}^{T}\tilde{\rho}^{(s)}}$

end for

Output: $\\{\rho^{(t)},\gamma^{(t)}\\}_{t=1}^{T}$.

The details of this algorithm is described in ALG 1. The associated estimator
for this sampler is defined as [11]:

$\displaystyle\texttt{PIP}(i)\approx\sum_{t=1}^{T}\rho^{(t)}\big{(}\mathbf{1}\\{i\in\mathcal{S}^{(t)}\\}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})+\mathbf{1}\\{i\notin\mathcal{S}^{(t)}\\}\gamma_{i}^{(t)}\big{)}.$
(13)

###  3.2 A Variable Complexity wTGS Scheme

In the subset wTGS in Subsection 3.1, the number of conditional PIP
computations per MCMC iteration is fixed, i.e., it is equal to $S$. In the
following, we propose a variable-computation complexity-based wTGS schemes
(VC-wTGS), say ALG 2, where the only requirement is that the expected number
of the conditional PIP computations per MCMC iteration is $S$. This means that
$\mathbb{E}[S_{t}]=S,$ where $S_{t}$ is the number of conditional PIP
computations at the $t$-th MCMC iteration.

Compared with ALG 1, ALG 2 allows us to use different subset sizes at MCMC
iterations. By ALG 2, the expectation of number of conditional PIP
computations in each MCMC iteration is $P\times(S/P)+0\times(1-S/P)=S$. Since
we aim to bound the variance at each finite iteration $T$, we don’t mention
about $T_{\rm{burn}}$ in ALG 2. In practice, we usually remove some initial
samples. We also use the following new version of Rao-Blackwellized estimator:

$\displaystyle\texttt{PIP}(i)\approx\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D}).$
(14)

Algorithm 2 A Variable-Complexity Based wTGS Algorithm

Input: Dataset $\mathcal{D}=\\{X,Y\\}$ with $P$ covariates; prior inclusion
probability $h$; prior precision $\tau$; total number of MCMC iterations $T$;
subset size $S$.

Output: Approximate weighted posterior samples
$\\{\rho^{(t)},\gamma^{(t)}\\}_{t=1}^{T}$

Initializations: $\gamma^{(0)}=(\gamma_{1},\gamma_{2},\cdots,\gamma_{P})$
where $\gamma_{j}\sim\texttt{Bern}(h)$ for all $j\in[P]$.

for $t=1,2,\cdots,T$ do

Set $Q^{(1)}=1$. Sample a Bernoulli random variable
$Q^{(t)}\sim\texttt{Bern}(\frac{S}{P})$ if $t\geq 2$.

if $Q^{(t)}=1$ then

Estimate
$f(j|\gamma^{(t-1)})\leftarrow\phi_{t-1}(\gamma)^{-1}\frac{\frac{1}{2}\eta(\gamma_{-j}^{(t-1)})}{p(\gamma_{j}^{(t-1)}|\gamma_{-j}^{(t-1)},\mathcal{D})}$
for all $j\in[P]$.

Sample $i^{(t)}\sim f(\cdot|\gamma^{(t-1)})$

$\gamma^{(t)}\leftarrow\texttt{flip}(\gamma^{(t-1)}|i^{(t)})$ where
$\texttt{flip}(\gamma|i)$ flips the $i$-th coordinate of
$\gamma:\gamma_{i}\leftarrow 1-\gamma_{i}$.

Estimate $P$ conditional PIPs
$p(\gamma_{j}^{(t)}|\gamma_{-j}^{(t)},\mathcal{D})$ for all $j\in[P]$

$\phi(\gamma^{(t)})\leftarrow\sum_{j\in[P]}\frac{\frac{1}{2}\eta(\gamma_{-j}^{(t)})}{p(\gamma_{j}^{(t)}|\gamma_{-j}^{(t)},\mathcal{D})}$

Compute the unnormalized weights
$\tilde{\rho}^{(t)}\leftarrow\phi^{-1}(\gamma^{(t)})$

else

$\gamma^{(t)}\leftarrow\gamma^{(t-1)}$

$\tilde{\rho}^{(t)}\leftarrow 1$

end if

end for

for $t=1,2,\cdots,T$ do

$\rho^{(t)}\leftarrow\frac{\tilde{\rho}^{(t)}Q^{(t)}}{\sum_{s=1}^{T}\tilde{\rho}^{(s)}Q^{(s)}}$

end for

Output: $\\{\rho^{(t)},\gamma^{(t)}\\}_{t=1}^{T}$.

######  Remark 1.

In ALG 2, Bernoulli random variables $\\{Q^{(t)}\\}_{t=1}^{T}$ are used to
replace for random set $\mathcal{S}$ in ALG 1. There are two main reasons for
this replacement: (1) generating a random set $\mathcal{S}$ from ${P\choose
S}$ subsets of $[P]$ takes very long running time for most pairs $(P,S)$, (2)
the associated Rao-Blackwellized estimator usually has smaller variance with
ALG 2 than ALG 1 at high $P/S$. See Section 4 for our simulation results.

###  3.3 Theoretical Bounds for Algorithm 2

First, we prove the following result. The proof can be found in Appendix A.

######  Lemma 1.

Let $U$ and $V$ be two positive random variables such that $U/V\leq M$ a.s.
for some constant $M$. In addition, assume that on a set $D$ with probability
at least $1-\alpha$, we have

$\displaystyle|U-\mathbb{E}[U]|$ $\displaystyle\leq\varepsilon\mathbb{E}[U],$
(15) $\displaystyle|V-\mathbb{E}[V]|$
$\displaystyle\leq\varepsilon\mathbb{E}[V],$ (16)

for some $0\leq\varepsilon<1$. Then, it holds that

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\frac{U}{V}-\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon^{2}}{(1-\varepsilon)^{2}}\bigg{(}\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{)}^{2}+\bigg{[}\max\bigg{(}M,\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{)}\bigg{]}^{2}\alpha.$
(17)

We also recall the following Hoeffding’s inequality for Markov chain:

######  Lemma 2.

[15, Theorem 1.1] Let $\\{Y_{i}\\}_{i=1}^{\infty}$ be a stationary Markov
chain with state space $[N]$, transition matrix $A$, stationary probability
measure $\pi$, and averaging operator $E_{\pi}$, so that $Y_{1}$ is
distributed according to $\pi$. Let $\lambda=\|A-E_{\pi}\|_{L_{2}(\pi)\to
L_{2}(\pi)}$ and let $f_{1},f_{2},\cdots,f_{n}:[N]\to\mathbb{R}$ so that
$\mathbb{E}[f_{i}(Y_{i})]=0$ for all $i$ and $|f_{i}(\nu)|\leq a_{i}$ for all
$\nu\in[N]$ and all $i$. Then for $u\geq 0$,

$\displaystyle\mathbb{P}\bigg{[}\bigg{|}\sum_{i=1}^{n}f_{i}(Y_{i})\bigg{|}\geq
u\bigg{(}\sum_{i=1}^{n}a_{i}^{2}\bigg{)}^{1/2}\bigg{]}\leq
2\exp\bigg{(}-\frac{u^{2}(1-\lambda)}{64e}\bigg{)}.$ (18)

Now, the following result can be shown.

######  Lemma 3.

Let

$\displaystyle\phi(\gamma):=\sum_{j\in[P]}\frac{\frac{1}{2}\eta(\gamma_{-j})}{p(\gamma_{j}|\gamma_{-j},\mathcal{D})}$
(19)

and define

$\displaystyle f(\gamma):=\phi(\gamma)p(\gamma|\mathcal{D}).$ (20)

Then, by ALG 2, the sequence $\\{\gamma^{(t)},Q^{(t)}\\}_{t=1}^{T}$ forms a
reversible Markov chain with the stationary distribution proportional to
$f(\gamma)q(Q)$ where $q$ is the Bernoulli $(S/P)$ distribution. This Markov
chain has transition kernel
$K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))=K^{*}(\gamma\to\gamma^{\prime})q(Q^{\prime})$
where

$\displaystyle
K^{*}(\gamma\to\gamma^{\prime})=\frac{S}{P}\sum_{j=1}^{P}f(j|\gamma)\delta(\gamma^{\prime}-\texttt{flip}(\gamma|j))+\bigg{(}1-\frac{S}{P}\bigg{)}\delta(\gamma^{\prime}-\gamma).$
(21)

In the classical wTGS [23], the Markov chain $\\{\gamma^{(t)}\\}_{t=1}^{T}$
also form a Markov chain. However, this Markov chain is different from the
Markov chain in Lemma 3. However, the two Markov chains still have the same
stationary distribution which is proportional to $f(\gamma)$. See a detailed
proof of Lemma 3 in Appendix B.

######  Lemma 4.

For the Rao-Blackwellized estimator in (14) which is applied to the output
sequence $\\{\rho^{(t)},\gamma^{(t)}\\}_{t=1}^{T}$ of ALG 2, it holds that

$\displaystyle
E_{i,T}:=\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})$
$\displaystyle\to\texttt{PIP}(i),\quad\texttt{as}\quad T\to\infty.$ (22)

###### Proof.

By Lemma 3, $\\{\gamma^{(t)},Q^{(t)}\\}_{t=1}^{T}$ forms a reversible Markov
chain with stationary distribution $f(\gamma)/Z_{f}q(Q)$ where
$Z_{f}=\sum_{\gamma}f(\gamma)$. Hence, by SLLN for Markov chain [4], for any
bounded function $h$, we have

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\phi^{-1}(\gamma^{(t)})Q^{(t)}h(\gamma^{(t)})$
$\displaystyle\to\mathbb{E}_{qf(\cdot)/Z_{f}}\big{[}\phi^{-1}(\gamma)h(\gamma)Q\big{]}$
(23)
$\displaystyle=\sum_{Q}q(Q)\sum_{\gamma}\frac{f(\gamma)}{Z_{f}}\phi^{-1}(\gamma)h(\gamma)Q$
(24)
$\displaystyle=\bigg{(}\sum_{Q}q(Q)Q\bigg{)}\bigg{(}\sum_{\gamma}\frac{f(\gamma)}{Z_{f}}\phi^{-1}(\gamma)h(\gamma)\bigg{)}$
(25)
$\displaystyle=\mathbb{E}_{q}[Q]\frac{1}{Z_{f}}\sum_{\gamma}p(\gamma|\mathcal{D})h(\gamma)$
(26)
$\displaystyle=\frac{S}{P}\frac{1}{Z_{f}}\sum_{\gamma}p(\gamma|\mathcal{D})h(\gamma),$
(27)

where (26) follows from $f(\gamma)=p(\gamma|\mathcal{D})\phi(\gamma)$.

Similarly, we have

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}Q^{(t)}\phi^{-1}(\gamma^{(t)})$
$\displaystyle\to\mathbb{E}_{qf(\cdot)/Z_{f}}\big{[}\phi^{-1}(\gamma)Q\big{]}$
(28)
$\displaystyle=\sum_{Q}q(Q)Q\sum_{\gamma}\frac{f(\gamma)}{Z_{f}}\phi^{-1}(\gamma)$
(29)
$\displaystyle=\mathbb{E}_{q}[Q]\sum_{\gamma}\frac{1}{Z_{f}}p(\gamma|\mathcal{D})$
(30) $\displaystyle=\frac{S}{P}\frac{1}{Z_{f}},$ (31)

where (30) also follows from $f(\gamma)=p(\gamma|D)\phi(\gamma)$. ∎

From (27) and (31), we obtain

$\displaystyle\frac{\frac{1}{T}\sum_{t=1}^{T}\phi^{-1}(\gamma^{(t)})Q^{(t)}h(\gamma^{(t)})}{\frac{1}{T}\sum_{t=1}^{T}Q^{(t)}\phi^{-1}(\gamma^{(t)})}\to\sum_{\gamma}p(\gamma|\mathcal{D})h(\gamma),$
(32)

or equivalently

$\displaystyle\sum_{t=1}^{T}\rho^{(t)}h(\gamma^{(t)})\to\sum_{\gamma}p(\gamma|\mathcal{D})h(\gamma)$
(33)

as $T\to\infty$.

Now, by setting $h(\gamma)=p(\gamma_{i}=1|\gamma_{-i},\mathcal{D})$, from
(33), we obtain

$\displaystyle\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})\to\texttt{PIP}(i)$
(34)

for all $i\in[P]$.

The following result bounds the variance of PIP estimator at finite $T$.

######  Lemma 5.

For any $\varepsilon\in[0,1]$, let $\nu$ and $\pi$ be the initial and
stationary distributions of the reversible Markov sequence
$\big{\\{}\big{(}\gamma^{(t)},Q^{(t)}\big{)}\big{\\}}$. Define

$\displaystyle\hat{\phi}(\gamma):=\frac{\phi^{-1}(\gamma)}{\max_{\gamma}\phi^{-1}(\gamma)},$
(35)

and

$\displaystyle\varepsilon_{0}=\frac{P}{\texttt{PIP}(i)\mathbb{E}_{\pi}[\hat{\phi}(\gamma)]S}\sqrt{\frac{64e\log
T}{(1-\lambda_{\gamma,Q})T}}.$ (36)

Then, we have

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})-\texttt{PIP}(i)\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon_{0}^{2}}{(1-\varepsilon_{0})^{2}}\texttt{PIP}^{2}(i)+\frac{4P}{S}\frac{1}{\min_{\gamma}\pi(\gamma)T}\to
0,$ (37)

as $T\to\infty$ for fixed $P,S$ and the dataset. Here, $\pi(\gamma)$ is the
marginal distribution of $\pi(\gamma,Q)$.

###### Proof.

See Appendix C. ∎

Next, we provide a lower bound for $1-\lambda_{\gamma,Q}$. First, we recall
the following Dirichlet form on spectral gap.

######  Definition 1.

Let $f,g:\Omega\to\mathbb{R}$. The Dirichlet form associated with a reversible
Markov chain $Q$ on $\Omega$ is defined by

$\displaystyle\mathcal{E}(f,g)$
$\displaystyle=\langle(\mathbf{I}-\mathbf{Q})f,g\rangle_{\pi}$ (38)
$\displaystyle=\sum_{x\in\Omega}\pi(x)[f(x)-\mathbf{Q}f(x)]g(x)$ (39)
$\displaystyle=\sum_{x\in\Omega}\pi(x)\bigg{[}\sum_{y}Q(x,y)(f(x)-f(y))\bigg{]}g(x)$
(40)
$\displaystyle=\sum_{x,y\in\Omega\times\Omega}\pi(x)Q(x,y)g(x)(f(x)-f(y)).$
(41)

######  Lemma 6.

[7] (Variational characterisation) For a reversible Markov chain $Q$ with
state space $\Omega$ and stationary distribution $\pi$, it holds that

$\displaystyle
1-\lambda=\inf_{g:\Omega\to\mathbb{R},\atop\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\mathcal{E}(g,g),$
(42)

where $\mathcal{E}(g,g):=\langle(\mathbf{I}-\mathbf{Q})g,g\rangle_{\pi}$.

######  Lemma 7.

The spectral gap $1-\lambda_{\gamma,Q}$ of the reversible Markov chain
$\\{\gamma^{(t)},Q^{(t)}\\}$ satisfies

$\displaystyle 1-\lambda_{\gamma,Q}\geq 1-\frac{S}{P}\lambda_{P}\geq
1-\frac{S}{P}$ (43)

where $1-\lambda_{P}$ is the spectral gap of the reversible Markov chain
$\\{\gamma^{(t)}\\}$ of the wTGS algorithm (i.e. $S=P$).

###### Proof.

From Lemma 6 and the fact that $\\{\gamma^{(t)},Q^{(t)}\\}$ forms a reversible
Markov chain with transition kernel
$K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))=K^{*}(\gamma\to\gamma^{\prime})q(Q^{\prime})$,
we have

$\displaystyle 1-\lambda_{\gamma,Q}$
$\displaystyle=\inf_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\langle
g,g\rangle_{\pi}-\langle\mathbf{K}g,g\rangle$ (44)
$\displaystyle=1-\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\langle\mathbf{K}g,g\rangle$
(45)
$\displaystyle=1-\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,Q}\mathbf{K}g(\gamma,Q)g(\gamma,Q)\pi(\gamma,Q)$
(46)
$\displaystyle=1-\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,Q}\sum_{\gamma^{\prime},Q^{\prime}}K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))g(\gamma^{\prime},Q^{\prime})g(\gamma,Q)\pi(\gamma,Q)$
(47)
$\displaystyle=1-\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\frac{S}{P}\sum_{\gamma,Q}\sum_{\gamma^{\prime},Q^{\prime}}K^{*}(\gamma\to\gamma^{\prime})q(Q^{\prime})g(\gamma^{\prime},Q^{\prime})g(\gamma,Q)\pi(\gamma,Q)$
(48)
$\displaystyle=1-\frac{S}{P}\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,Q}\sum_{\gamma^{\prime},Q^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\frac{f(\gamma)}{Z_{f}}q(Q)g(\gamma^{\prime},Q^{\prime})g(\gamma,Q)q(Q^{\prime})$
(49)
$\displaystyle=1-\frac{S}{P}\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,\gamma^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\frac{f(\gamma)}{Z_{f}}\sum_{Q,Q^{\prime}}g(\gamma^{\prime},Q^{\prime})g(\gamma,Q)q(Q)q(Q^{\prime})$
(50)
$\displaystyle=1-\frac{S}{P}\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,\gamma^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\pi(\gamma)\bigg{(}\sum_{Q}g(\gamma,Q)q(Q)\bigg{)}\bigg{(}\sum_{Q^{\prime}}\pi(\gamma^{\prime},Q^{\prime})q(Q^{\prime})\bigg{)}$
(51)
$\displaystyle=1-\frac{S}{P}\sup_{g(\gamma,Q):\mathbb{E}_{\pi}[g]=0,\mathbb{E}_{\pi}[g^{2}]=1}\sum_{\gamma,\gamma^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\pi(\gamma)h(\gamma)h(\gamma^{\prime})$
(52)

where

$\displaystyle\pi(\gamma)$ $\displaystyle=\frac{f(\gamma)}{Z_{f}},$ (53)
$\displaystyle Z_{f}$ $\displaystyle=\sum_{\gamma}f(\gamma),$ (54)
$\displaystyle h(\gamma)$ $\displaystyle:=\sum_{Q}g(\gamma,Q)q(Q).$ (55)

Observe that

$\displaystyle\mathbb{E}_{\pi}[h(\gamma)]$
$\displaystyle=\sum_{\gamma}h(\gamma)\pi(\gamma)$ (56)
$\displaystyle=\sum_{\gamma}\sum_{Q}g(\gamma,Q)q(Q)\pi(\gamma)$ (57)
$\displaystyle=\sum_{\gamma,Q}g(\gamma,Q)\pi(\gamma,Q)$ (58)
$\displaystyle=\mathbb{E}_{\pi}[g(\gamma,Q)]$ (59) $\displaystyle=0.$ (60)

On the other hand, we also have

$\displaystyle\mathbb{E}_{\pi}\big{[}h^{2}(\gamma)\big{]}$
$\displaystyle=\sum_{\gamma}\bigg{(}\sum_{Q}g(\gamma,Q)q(Q)\bigg{)}^{2}\pi(\gamma)$
(61)
$\displaystyle\leq\sum_{\gamma}\bigg{(}\sum_{Q}g(\gamma,Q)^{2}q(Q)\bigg{)}\pi(\gamma)$
(62) $\displaystyle=\sum_{\gamma,Q}g(\gamma,Q)^{2}\pi(\gamma,Q)$ (63)
$\displaystyle=\mathbb{E}_{\pi}\big{[}g(\gamma,Q)^{2}\big{]}$ (64)
$\displaystyle=1,$ (65)

where (62) follows from the convexity of the function $x^{2}$ on $[0,\infty)$.

From (60), (65), and (52), we obtain

$\displaystyle 1-\lambda_{\gamma,Q}$ $\displaystyle\geq
1-\frac{S}{P}\sup_{h(\gamma):\mathbb{E}_{\pi}[h]=0,\mathbb{E}_{\pi}[h^{2}]\leq
1}\sum_{\gamma,\gamma^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\pi(\gamma)h(\gamma)h(\gamma^{\prime}).$
(66)

Now, note that $\mathbb{E}_{\pi}[h]=0$ is equivalent to
$h\perp_{\pi}\mathbf{1}$. Let $|\Omega|=2^{P+1}:=n$ and
$h_{1},h_{2},\cdots,h_{n}$ are eigenfunctions of $\mathbf{K}^{*}$
corresponding to the decreasing ordered eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ and are orthogonal since
$\mathbf{K}^{*}$ is self-adjoint. Set $h_{1}=\mathbf{1}$. Since
$\|h\|_{2,\pi}=1$ and $h\perp_{\pi}\mathbf{1}$, we have
$h=\sum_{j=2}^{n}a_{j}h_{j}$ because it is perpendicular to $h_{1}$ so it can
be only represented by these eigenvectors. By taking $l_{2}$-norm on both
sizes we have $\sum_{j=2}^{n}a_{j}^{2}\leq 1$ since the form like $\langle
h_{i},h_{j}\rangle_{\pi}=0$ and $\langle
h_{i},h_{i}\rangle=\|h_{i}\|_{2,\pi}^{2}=1$. Thus,

$\displaystyle\sup_{h:\mathbb{E}_{\pi}[h]=0,\mathbb{E}_{\pi}[h^{2}]\leq
1}\sum_{\gamma,\gamma^{\prime}}K^{*}(\gamma\to\gamma^{\prime})\pi(\gamma)h(\gamma)h(\gamma^{\prime})$
$\displaystyle\leq\max_{a_{2},a_{3},\cdots,a_{n}:\sum_{j=2}^{n}a_{j}^{2}\leq
1}\sum_{j=1}^{n}a_{j}^{2}\lambda_{j}$ (67)
$\displaystyle\leq\lambda_{2}\sum_{j=2}^{n}a_{j}^{2}$ (68)
$\displaystyle=\lambda_{2},$ (69)

where $\sum_{j=2}^{n}a_{j}^{2}\leq 1$ and $\lambda_{j}\in\texttt{spec}(P)$
such that $\lambda_{2}\geq\lambda_{3}\cdots\geq\lambda_{n}$. Hence, from (69),
we obtain

$\displaystyle 1-\lambda_{\gamma,Q}$ $\displaystyle\geq
1-\frac{S}{P}\lambda_{2}$ (70) $\displaystyle=1-\frac{S}{P}\lambda_{P}$ (71)
$\displaystyle=\frac{S}{P}\big{(}1-\lambda_{P}\big{)}+1-\frac{S}{P}$ (72)
$\displaystyle\geq 1-\frac{S}{P}.$ (73)

∎

By combining Lemma 4, Lemma 5 and Lemma 7, we come up with the following
theorem.

######  Theorem 1.

For the variable-complexity subset wTGS-based estimator in (14) and given
dataset $(X,Y)$, it holds that

$\displaystyle
E_{i,T}:=\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})$
$\displaystyle\to\texttt{PIP}(i),\quad\texttt{as}\quad T\to\infty$ (74)

and

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}|\gamma_{-i}^{(t)},\mathcal{D})-\texttt{PIP}(i)\bigg{|}^{2}\bigg{]}=O\bigg{(}\frac{\log
T}{T}\bigg{(}\frac{P}{S}\bigg{)}^{2}\bigg{(}\frac{\max_{\gamma}\phi(\gamma)}{\min_{\gamma}\phi(\gamma)}\bigg{)}^{2}\bigg{)},$
(75)

where

$\displaystyle\phi(\gamma)=\frac{1}{2}\sum_{j\in[P]}\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}|\gamma_{-j},\mathcal{D})}.$
(76)

###### Proof.

First, (74) is shown in Lemma 4. Now, we show (75) by using Lemma 5 and Lemma
7.

Observe that

$\displaystyle\mathbb{E}_{\pi}[\hat{\phi}(\gamma)]$
$\displaystyle=\mathbb{E}_{\pi}\bigg{[}\frac{\phi^{-1}(\gamma)}{\max_{\gamma}\phi^{-1}(\gamma)}\bigg{]}$
$\displaystyle\geq\frac{\min_{\gamma}\phi^{-1}(\gamma)}{\max_{\gamma}\phi^{-1}(\gamma)}$
(77)
$\displaystyle=\frac{\min_{\gamma}\phi(\gamma)}{\max_{\gamma}\phi(\gamma)}.$
(78)

In addition, we have

$\displaystyle\phi(\gamma)$
$\displaystyle=\sum_{j\in[P]}\frac{\frac{1}{2}\eta(\gamma_{-j})}{p(\gamma_{j}|\gamma_{-j},\mathcal{D})}$
(79)
$\displaystyle=\frac{1}{2}\sum_{j\in[P]}\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}|\gamma_{-j},\mathcal{D})}.$
(80)

Now, note that

$\displaystyle\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}|\gamma_{-j},\mathcal{D})}=\begin{cases}1,&\qquad\gamma_{j}=1\\\
\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}=0|\gamma_{-j},\mathcal{D})},&\qquad\gamma_{j}=0.\end{cases}$
(81)

In Appendix D, we can estimate
$\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}=0|\gamma_{-j},\mathcal{D})}$
based on the dataset. More specifically, let $\tilde{\gamma}_{1}$ is given by
$\gamma_{-i}$ with $\gamma_{i}=1$, $\tilde{\gamma}_{0}$ is given by
$\gamma_{-i}$ with $\gamma_{i}=0$, then we can show that

$\displaystyle\frac{p(\gamma_{j}=1|\gamma_{-j},\mathcal{D})}{p(\gamma_{j}=0|\gamma_{-j},\mathcal{D})}=\bigg{(}\frac{h}{1-h}\bigg{)}\sqrt{\tau\frac{\det(X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I)}{\det(X_{\tilde{\gamma}_{1}}^{T}X_{\tilde{\gamma}_{1}}+\tau
I)}}\bigg{(}\frac{\|Y\|^{2}-\|\tilde{Y}_{\tilde{\gamma}_{0}}\|^{2}+\nu_{0}\lambda_{0}}{\|Y\|^{2}-\|\tilde{Y}_{\tilde{\gamma}_{1}}\|^{2}+\nu_{0}\lambda_{0}}\bigg{)}^{\frac{N+\nu_{0}}{2}}.$
(82)

Here,

$\displaystyle\|\tilde{Y}_{\gamma}\|^{2}$
$\displaystyle=\tilde{Y}_{\gamma}^{T}\tilde{Y}_{\gamma}$ (83)
$\displaystyle=Y^{T}X_{\gamma}(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1}X_{\gamma}^{T}Y.$ (84)

∎

##  4 Experiments

In this section, we show by simulation that the PIP-estimator is convergent as
$T\to\infty$. In addition, we compare the variance of associated Rao-
Blackwellized estimators for VC-wTGS and subset wTGS on simulated and real
datasets. To compute $p(\gamma_{i}|\gamma_{-i},Y)$, we use the same trick as
[23, Appendix B.1] for the new setting. See our derivations of this posterior
distribution in Appendix D.

Figure 1: VC-wTGS Rao-Blackwellized Estimators (ALG 2)

Figure 2: The variance of VC-wTGS Rao-Blackwellized Estimators

###  4.1 Simulated Datasets

First, we perform a simulated experiment. Let $X\in\mathbb{R}^{N\times P}$ be
a realization of a multivariate (random) Gaussian matrix. We consider the case
$N=100$ and $P=200$. We run $T=20000$ iterations.

Fig. 1 shows that the Rao-Blackwellized estimator in (14) converges to the
value of PIP at $T\to\infty$ for different values of $S$. Fig. 2 shows the
variance for estimating PIP(0) and PIP(1) at different values of $S$. Since
the variance of wTGS is very small $\approx 0$ at $T$ big enough, the variance
of variable-complexity wTGS is also small at $T$ big enough.

###  4.2 Real Datasets

In this simulation, we run ALG 2 on MNIST dataset. Fig. 3 and Fig. 4 plot
PIP(0) and PIP(1) and the estimated variances for the Rao-Blackwellized
estimator in (14) at different values of $S$, respectively. These plots show a
trade-off between the computational complexity and the estimated variance for
estimating PIP(0) and PIP(1). The expected number of PIP computations is only
$ST$ in ALG 2 but $TP$ in wTGS if we run $T$ MCMC iterations. However, we
suffer an increasing in variance. By Theorem 1, the variance is
$O\big{(}\big{(}\frac{P}{S}\big{)}^{2}\frac{\log T}{T}\big{)}$ for a given
dataset, i.e., increasing at most $(P/S)^{2}$ times. For many applications, we
don’t need to estimate PIPs exactly, hence VC-wTGS can be used to reduce
computational complexity especially when $P$ is very large (million
covariates). Fig. 5 shows that VC-wTGS outperforms subset wTGS [11] at high
values of $P/S$.

Figure 3: VC-wTGS Rao-Blackwellized Estimators (ALG 2)

Figure 4: The variance of VC-wTGS Rao-Blackwellized Estimators

Figure 5: Comparing the variance between subset wTGS and VC-wTGS at $S=2$.

##  Appendix A Proof of Lemma 1

Observe that with probability at least $1-\alpha$, we have

$\displaystyle(1-\varepsilon)\mathbb{E}[U]$ $\displaystyle\leq
U\leq(1+\varepsilon)\mathbb{E}[U]$ (85)
$\displaystyle(1-\varepsilon)\mathbb{E}[V]$ $\displaystyle\leq
V\leq(1+\varepsilon)\mathbb{E}[V].$ (86)

Hence, we have

$\displaystyle\bigg{(}\frac{1-\varepsilon}{1+\varepsilon}\bigg{)}\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\leq\frac{U}{V}\leq\bigg{(}\frac{1+\varepsilon}{1-\varepsilon}\bigg{)}\frac{\mathbb{E}[U]}{\mathbb{E}[V]}.$
(87)

From (87), with probability at least $1-\alpha$, we have

$\displaystyle\bigg{|}\frac{U}{V}-\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{|}\leq\frac{2\varepsilon}{1-\varepsilon}\bigg{(}\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{)}.$
(88)

It follows from (88) that

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\frac{U}{V}-\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{|}^{2}\bigg{]}$
$\displaystyle=\mathbb{E}\bigg{[}\bigg{|}\frac{U}{V}-\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{|}^{2}\bigg{|}D\bigg{]}\mathbb{P}(D)+\mathbb{E}\bigg{[}\bigg{|}\frac{U}{V}-\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{|}^{2}\bigg{|}D^{c}\bigg{]}\mathbb{P}(D^{c})$
(89)
$\displaystyle\leq\frac{4\varepsilon^{2}}{(1-\varepsilon)^{2}}\bigg{(}\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{)}^{2}+\bigg{[}\max\bigg{(}M,\frac{\mathbb{E}[U]}{\mathbb{E}[V]}\bigg{)}\bigg{]}^{2}\alpha.$
(90)

##  Appendix B Proof of Lemma 3

The transition kernel for the sequence $\\{\gamma^{(t)}\\}$ can be written as

$\displaystyle
K^{*}(\gamma\to\gamma^{\prime})=\frac{S}{P}\sum_{j=1}^{P}f(j|\gamma)\delta(\gamma^{\prime}-\texttt{flip}(\gamma|j))+\bigg{(}1-\frac{S}{P}\bigg{)}\delta(\gamma^{\prime}-\gamma).$
(91)

This implies that for any pair $(\gamma,\gamma^{\prime})$ such that
$\gamma^{\prime}=\texttt{flip}(\gamma|i)$ for some $i\in[P]$, we have

$\displaystyle K^{*}(\gamma\to\gamma^{\prime})$
$\displaystyle=\frac{S}{P}\sum_{j=1}^{P}f(j|\gamma)\delta(\gamma^{\prime}-\texttt{flip}(\gamma|j))$
(92) $\displaystyle=\frac{S}{P}f(i|\gamma).$ (93)

Now, by ALG 2, we also have

$\displaystyle f(i|\gamma)$
$\displaystyle=\phi^{-1}(\gamma)\frac{\frac{1}{2}\eta(\gamma_{-i})}{p(\gamma_{i}|\gamma_{-i},\mathcal{D})}$
(94)

and

$\displaystyle f(i|\gamma^{\prime})$
$\displaystyle=\phi^{-1}(\gamma^{\prime})\frac{\frac{1}{2}\eta(\gamma^{\prime}_{-i})}{p(\gamma^{\prime}_{i}|\gamma^{\prime}_{-i},\mathcal{D})}.$
(95)

From (94) and (95) and $\gamma_{-i}=\gamma^{\prime}_{-i}$, we obtain

$\displaystyle\frac{K^{*}(\gamma\to\gamma^{\prime})}{K^{*}(\gamma^{\prime}\to\gamma)}$
$\displaystyle=\frac{\frac{S}{P}f(i|\gamma)}{\frac{S}{P}f(i|\gamma^{\prime})}$
(96) $\displaystyle=\frac{f(i|\gamma)}{f(i|\gamma^{\prime})}$ (97)
$\displaystyle=\frac{\phi(\gamma^{\prime})p(\gamma^{\prime}|\mathcal{D})}{\phi(\gamma)p(\gamma|\mathcal{D})}$
(98) $\displaystyle=\frac{f(\gamma^{\prime})}{f(\gamma)}.$ (99)

In addition, we also have
$K^{*}(\gamma\to\gamma^{\prime})=K^{*}(\gamma^{\prime}\to\gamma)=0$ if
$\gamma^{\prime}\neq\gamma$ and $\gamma^{\prime}\neq\texttt{flip}(\gamma|i)$
for any $i\in[P]$. Furthermore,
$K^{*}(\gamma\to\gamma^{\prime})=K^{*}(\gamma^{\prime}\to\gamma)=1-\frac{S}{P}$
if $\gamma=\gamma^{\prime}$.

By combining all these cases, it holds that

$\displaystyle
f(\gamma)K^{*}(\gamma\to\gamma^{\prime})=f(\gamma^{\prime})K^{*}(\gamma^{\prime}\to\gamma)$
(100)

for all $\gamma^{\prime},\gamma$.

This means that $\\{\gamma^{(t)}\\}_{t=1}^{T}$ form a reversible Markov chain
with stationary distribution $f(\gamma)/Z_{f}$ where

$\displaystyle Z_{f}=\sum_{\gamma}f(\gamma).$ (101)

Since $\\{Q_{t}\\}_{t=1}^{T}$ is an i.i.d. Bernoulli sequence with $q(1)=S/P$
and independent of $\\{\gamma^{(t)}\\}_{t=1}^{T}$,
$\\{\gamma^{(t)},Q^{(t)}\\}_{t=1}^{T}$ forms a Markov chain with the
transition kernel satisfying:

$\displaystyle
K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))=q(Q^{\prime})K^{*}(\gamma\to\gamma^{\prime}).$
(102)

It follows from (102) that

$\displaystyle
q(Q)f(\gamma)/Z_{f}K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))=\big{[}K^{*}(\gamma\to\gamma^{\prime})f(\gamma)/Z_{f}\big{]}q(Q)q(Q^{\prime})$
(103)

for any pair $(\gamma,Q)$ and $(\gamma^{\prime},Q^{\prime})$.

Finally, from (100) and (103), we have

$\displaystyle
q(Q)f(\gamma)/Z_{f}K((\gamma,Q)\to(\gamma^{\prime},Q^{\prime}))=q(Q^{\prime})f(\gamma)/Z_{f}K((\gamma^{\prime},Q^{\prime})\to(\gamma,Q)).$
(104)

This means that $\\{\gamma_{t},Q^{(t)}\\}_{t=1}^{T}$ forms a reversible Markov
chain with stationary distribution $q(Q)f(\gamma)/Z_{f}$.

##  Appendix C Proof of Lemma 5

First, by definition of $\hat{\phi}(\gamma)$ in (35) we have

$\displaystyle\rho^{(t)}=\frac{\hat{\phi}(\gamma^{(t)})}{\sum_{t=1}^{T}\hat{\phi}\big{(}\gamma^{(t)}\big{)}}.$
(105)

In addition, observe that

$\displaystyle 0\leq\hat{\phi}(\gamma)\leq 1.$ (106)

Now, let $g:\\{0,1\\}^{P}\to\mathbb{R}_{+}$ such that $g(\gamma)\leq 1$ for
all $\gamma$. Then, by applying Lemma 2 and a change of measure, with
probability $1-2\frac{d\nu}{d\pi}\exp(-\frac{\zeta^{2}T(1-\lambda)}{64e})$, we
have

$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}\bigg{]}\bigg{|}\leq\zeta$
(107)

for any $\zeta>0$.

Similarly, by using Lemma 2, with probability at least
$1-2\frac{d\nu}{d\pi}\exp(-\frac{\zeta^{2}T(1-\lambda)}{64e})$, it holds that

$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)}-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)}\bigg{]}\bigg{|}\leq\zeta.$
(108)

By using the union bound, with probability at least
$1-4\frac{d\nu}{d\pi}\exp(-\frac{\zeta^{2}T(1-\lambda)}{64e})$, it holds that

$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}\bigg{]}\bigg{|}\leq\zeta,$
(109)
$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})\bigg{]}\bigg{|}\leq\zeta.$
(110)

Now, by setting
$\zeta=\zeta_{0}:=\frac{\varepsilon}{T}\min\big{\\{}\mathbb{E}_{\pi}\big{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}\big{]},\mathbb{E}_{\pi}\big{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})\big{]}\big{\\}}$
for some $\varepsilon>0$ (to be chosen later), with probability at least
$1-4\frac{d\nu}{d\pi}\exp(-\frac{\zeta_{0}^{2}T(1-\lambda)}{64e})$, it holds
that

$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}\bigg{]}\bigg{|}$
$\displaystyle\qquad\leq\frac{\varepsilon}{T}\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}\bigg{]},$
(111)
$\displaystyle\frac{1}{T}\bigg{|}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)}-\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)}\bigg{]}\bigg{|}\leq\frac{\varepsilon}{T}\mathbb{E}_{\pi}\bigg{[}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)}\bigg{]}.$
(112)

Furthermore, by setting

$\displaystyle U$
$\displaystyle:=\frac{1}{T}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)},$
(113) $\displaystyle V$
$\displaystyle:=\frac{1}{T}\sum_{t=1}^{T}\hat{\phi}(\gamma^{(t)})Q^{(t)},$
(114)

we have

$\displaystyle\frac{U}{V}$
$\displaystyle=\frac{\sum_{t=1}^{T}\phi^{-1}(\gamma^{(t)})g(\gamma^{(t)})Q^{(t)}}{\sum_{t=1}^{T}\phi^{-1}(\gamma^{(t)})Q^{(t)}}$
(115) $\displaystyle=\sum_{t=1}^{T}\rho^{(t)}g(\gamma^{(t)})$ (116)

and

$\displaystyle M:=\sup(U/V)\leq 1$ (117)

since $\sum_{t=1}^{T}\rho^{(t)}=1$ and $g(\gamma^{(t)})\leq 1$ for all
$\gamma^{(t)}$.

From (110)-(117), by Lemma 1, we have

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}g(\gamma^{(t)})Q^{(t)}-\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon^{2}}{(1-\varepsilon)^{2}}\bigg{(}\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{)}^{2}+\bigg{[}\max\bigg{(}1,\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{)}\bigg{]}^{2}\alpha,$
(118)

where
$\alpha:=4\frac{d\nu}{d\pi}\exp\big{(}-\frac{\varepsilon^{2}T(1-\lambda_{\gamma,Q})\min\\{\mathbb{E}_{\pi}[U],\mathbb{E}_{\pi}[V]\\}^{2}}{64e}\big{)}$,
where $\lambda_{\gamma,Q}$ is the stationary distribution of the reversible
Markov chain $\\{\gamma^{(t)},Q^{(t)}\\}$.

Now, by setting

$\displaystyle\varepsilon=\varepsilon_{0}$
$\displaystyle=\frac{1}{\min\\{\mathbb{E}_{\pi}[U],\mathbb{E}_{\pi}[V]\\}}\sqrt{\frac{64e\log
T}{(1-\lambda_{\gamma,Q})T}},$ (119)

we have $\alpha=4\frac{d\nu}{d\pi}\frac{1}{T}$. Then, we obtain

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}g(\gamma^{(t)})-\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon_{0}^{2}}{(1-\varepsilon_{0})^{2}}\bigg{(}\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{)}^{2}+\bigg{[}\max\bigg{(}1,\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{)}\bigg{]}^{2}\alpha.$
(120)

Now, observe that

$\displaystyle\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}$
$\displaystyle=\frac{\mathbb{E}_{\pi}\big{[}g(\gamma)Q\hat{\phi}(\gamma)\big{]}}{\mathbb{E}_{\pi}\big{[}\hat{\phi}(\gamma)Q\big{]}}$
(121)
$\displaystyle=\frac{\mathbb{E}_{\pi}\big{[}g(\gamma)Q\phi^{-1}(\gamma)\big{]}}{\mathbb{E}_{\pi}\big{[}\phi^{-1}(\gamma)Q\big{]}}.$
(122)

On the other hand, by Lemma 3, we have
$\pi(\gamma,Q)=\frac{q(Q)f(\gamma)}{Z_{f}}$ where
$Z_{f}:=\sum_{\gamma}f(\gamma)$ and
$f(\gamma)=p(\gamma|\mathcal{D})\phi(\gamma)$. It follows that

$\displaystyle\mathbb{E}_{\pi}\big{[}g(\gamma)Q\phi^{-1}(\gamma)\big{]}$
$\displaystyle=\mathbb{E}_{q(Q)f(\gamma)/Z_{f}}\big{[}g(\gamma)Q\phi^{-1}(\gamma)\big{]}$
(123)
$\displaystyle=\sum_{\gamma}\sum_{Q}g(\gamma)Q\phi^{-1}(\gamma)\frac{f(\gamma)}{Z_{f}}q(Q)$
(124)
$\displaystyle=\frac{1}{Z_{f}}\sum_{\gamma}\sum_{Q}g(\gamma)q(Q)Qp(\gamma|\mathcal{D})$
(125)
$\displaystyle=\frac{1}{Z_{f}}\mathbb{E}_{p(\gamma|\mathcal{D})}\big{[}g(\gamma)\big{]}\mathbb{E}_{q}[Q].$
(126)

Similarly, we have

$\displaystyle\mathbb{E}_{\pi}\big{[}\phi^{-1}(\gamma)Q\big{]}$
$\displaystyle=\mathbb{E}_{q(Q)f(\gamma)/Z_{f}}\big{[}\phi^{-1}(\gamma)Q\big{]}$
(127)
$\displaystyle=\sum_{Q}\sum_{\gamma}\phi^{-1}(\gamma)Q\frac{f(\gamma)}{Z_{f}}q(Q)$
(128)
$\displaystyle=\frac{1}{Z_{f}}\bigg{(}\sum_{\gamma}P(\gamma|\mathcal{D})\bigg{)}\mathbb{E}_{q}[Q].$
(129)

From (122), (126) and (129), we obtain

$\displaystyle\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}=\mathbb{E}_{p(\gamma|\mathcal{D})}\big{[}g(\gamma)\big{]}.$
(130)

For the given problem, by setting
$g(\gamma)=p(\gamma_{i}=1|\gamma_{-i},\mathcal{D})$, from (130), we have

$\displaystyle\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}=\texttt{PIP}(i).$
(131)

In addition, we have

$\displaystyle\mathbb{E}_{\pi}[V]$
$\displaystyle=\mathbb{E}_{\pi}\big{[}\hat{\phi}(\gamma)Q\big{]}$ (132)
$\displaystyle=\sum_{\gamma,Q}\hat{\phi}(\gamma)Q\frac{f(\gamma)}{Z_{f}}q(Q)$
(133)
$\displaystyle=\bigg{(}\sum_{\gamma}\hat{\phi}(\gamma)\frac{f(\gamma)}{Z_{f}}\bigg{)}\bigg{(}\sum_{Q}Qq(Q)\bigg{)}$
(134) $\displaystyle=\mathbb{E}_{\pi}[\hat{\phi}(\gamma)]\mathbb{E}_{Q}[Q]$
(135) $\displaystyle=\frac{S}{P}\mathbb{E}_{\pi}[\hat{\phi}(\gamma)].$ (136)

Hence, we obtain

$\displaystyle\min\\{\mathbb{E}_{\pi}[U],\mathbb{E}_{\pi}[V]\\}$
$\displaystyle=\mathbb{E}_{\pi}[V]\min\bigg{\\{}1,\frac{\mathbb{E}_{\pi}[U]}{\mathbb{E}_{\pi}[V]}\bigg{\\}}$
(137)
$\displaystyle=\mathbb{E}_{\pi}[V]\min\bigg{\\{}1,\texttt{PIP}(i)\bigg{\\}}$
(138) $\displaystyle=\mathbb{E}_{\pi}[V]\texttt{PIP}(i)$ (139)
$\displaystyle=\frac{S}{P}\mathbb{E}_{\pi}[\hat{\phi}(\gamma)]\texttt{PIP}(i).$
(140)

From (120), (131), and (140), we have

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})-\texttt{PIP}(i)\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon_{0}^{2}}{(1-\varepsilon_{0})^{2}}\texttt{PIP}^{2}(i)+4\frac{d\nu}{d\pi}\frac{1}{T},$
(141)

and

$\displaystyle\varepsilon_{0}=\frac{P}{\texttt{PIP}(i)\mathbb{E}_{\pi}[\hat{\phi}(\gamma)]S}\sqrt{\frac{64e\log
T}{(1-\lambda_{\gamma,Q})T}}.$ (142)

Now, observe that

$\displaystyle\frac{d\nu}{d\pi}(\gamma,Q)$
$\displaystyle=\frac{p_{\gamma_{1},Q_{1}}(\gamma,Q)}{\pi(\gamma,Q)}$ (143)
$\displaystyle\leq\frac{1}{\pi(\gamma,Q)}$ (144)
$\displaystyle=\frac{1}{\pi(\gamma)q(Q)}$ (145)
$\displaystyle\leq\frac{P}{S}\frac{1}{\min_{\gamma}\pi(\gamma)}.$ (146)

By combining (141) and (146), we have

$\displaystyle\mathbb{E}\bigg{[}\bigg{|}\sum_{t=1}^{T}\rho^{(t)}p(\gamma_{i}^{(t)}=1|\gamma_{-i}^{(t)},\mathcal{D})-\texttt{PIP}(i)\bigg{|}^{2}\bigg{]}\leq\frac{4\varepsilon_{0}^{2}}{(1-\varepsilon_{0})^{2}}\texttt{PIP}^{2}(i)+\frac{4P}{S}\frac{1}{\min_{\gamma}\pi(\gamma)T}.$
(147)

##  Appendix D Derive $p(\gamma_{i}|\mathcal{D},\gamma_{-i})$

Observe that

$\displaystyle p(\gamma_{i}|\mathcal{D},\gamma_{-i})$
$\displaystyle=\frac{p(\gamma_{i}|\mathcal{D},\gamma_{-i})}{p(1-\gamma_{i}|\mathcal{D},\gamma_{-i})}\bigg{(}1+\frac{p(\gamma_{i}|\mathcal{D},\gamma_{-i})}{p(1-\gamma_{i}|\mathcal{D},\gamma_{-i})}\bigg{)}^{-1}.$
(148)

In addition, we have

$\displaystyle\frac{p(\gamma_{i}=1|\mathcal{D},\gamma_{-i})}{p(\gamma_{i}=0|\mathcal{D},\gamma_{-i})}$
$\displaystyle=\frac{p(\gamma_{i}=1,\mathcal{D}|\gamma_{-i})}{p(\gamma_{i}=0,\mathcal{D}|\gamma_{-i})}$
(149)
$\displaystyle=\frac{p(\gamma_{i}=1|\gamma_{-i},X)}{p(\gamma_{i}=0|\gamma_{-i},X)}\frac{p(Y|\gamma_{i}=1,\gamma_{-i},X)}{p(Y|\gamma_{i}=0,\gamma_{-i},X)}$
(150)
$\displaystyle=\bigg{(}\frac{p(\gamma_{i}=1)}{p(\gamma_{i}=0)}\bigg{)}\bigg{(}\frac{p(Y|\gamma_{i}=1,\gamma_{-i},X)}{p(Y|\gamma_{i}=0,\gamma_{-i},X)}\bigg{)}$
(151)
$\displaystyle=\bigg{(}\frac{h}{1-h}\bigg{)}\bigg{(}\frac{p(Y|\gamma_{i}=1,\gamma_{-i},X)}{p(Y|\gamma_{i}=0,\gamma_{-i},X)}\bigg{)}.$
(152)

On the other hand, for any tuple
$\gamma=(\gamma_{1},\gamma_{2},\cdots,\gamma_{P})$ such that $\gamma_{i}=1$
(so $|\gamma|\geq 1$), we have

$\displaystyle
p(Y|\gamma_{i}=1,\gamma_{-i},\beta_{\gamma},\sigma_{\gamma}^{2},X)$
$\displaystyle=\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\|Y-X_{\gamma}\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}.$
(153)

It follows that

$\displaystyle p(Y|\gamma_{i}=1,\gamma_{-i},X\big{)}$
$\displaystyle=\int_{\beta_{\gamma}}\int_{\sigma_{\gamma}^{2}=0}^{\infty}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\|Y-X_{\gamma}\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}p(\beta_{\gamma}|\gamma_{i}=1,\gamma_{-i})p(\sigma_{\gamma}^{2}|\gamma_{i}=1,\gamma_{-i})d\beta_{\gamma}d\sigma_{\gamma}^{2}$
(154)
$\displaystyle=\int_{\sigma_{\gamma}^{2}=0}^{\infty}\texttt{InvGamma}\bigg{(}\frac{1}{2}\nu_{0},\frac{1}{2}\nu_{0}\lambda_{0}\bigg{)}\int_{\beta_{\gamma}}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\|Y-X_{\gamma}\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}$
$\displaystyle\qquad\times\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi\tau^{-1}}\big{)}^{|\gamma|}}\exp\bigg{(}-\frac{\|\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}\tau^{-1}}\bigg{)}d\beta_{\gamma}d\sigma_{\gamma}^{2}.$
(155)

Now, observe that

$\displaystyle\|Y-X_{\gamma}\beta_{\gamma}\|^{2}+\tau\|\beta_{\gamma}\|^{2}$
$\displaystyle\qquad=(Y-X_{\gamma}\beta_{\gamma})^{T}(Y-X_{\gamma}\beta_{\gamma})+\tau\beta_{\gamma}^{T}\beta_{\gamma}$
(156)
$\displaystyle\qquad=Y^{T}Y-2Y^{T}X_{\gamma}\beta_{\gamma}+\beta_{\gamma}^{T}X_{\gamma}^{T}X_{\gamma}\beta_{\gamma}+\tau\beta_{\gamma}^{T}\beta_{\gamma}$
(157)
$\displaystyle\qquad=Y^{T}Y-2Y^{T}X_{\gamma}\beta_{\gamma}+\beta_{\gamma}^{T}(X_{\gamma}^{T}X_{\gamma}+\tau
I)\beta_{\gamma}.$ (158)

Now, consider the EVD (singular value decomposition) of the positive definite
matrix $X_{\gamma}^{T}X_{\gamma}+\tau I$ (note that $\tau>0$):

$\displaystyle X_{\gamma}^{T}X_{\gamma}+\tau I=U^{T}\Lambda U$ (159)

where $\Lambda$ is the a diagonal matrix consisting of all positive eigenvalue
of $X_{\gamma}^{T}X_{\gamma}+\tau I$. Let

$\displaystyle\tilde{\beta}_{\gamma}$
$\displaystyle:=\sqrt{\Lambda}U\beta_{\gamma},$ (160)
$\displaystyle\tilde{Y}_{\gamma}$
$\displaystyle:=\sqrt{\Lambda^{-1}}UX_{\gamma}^{T}Y.$ (161)

Then, we have

$\displaystyle\|Y-X_{\gamma}\beta_{\gamma}\|^{2}+\tau\|\beta_{\gamma}\|^{2}$
$\displaystyle\qquad=Y^{T}Y-2Y^{T}X_{\gamma}\beta_{\gamma}+\beta_{\gamma}^{T}(X_{\gamma}^{T}X_{\gamma}+\tau
I)\beta_{\gamma}$ (162)
$\displaystyle\qquad=Y^{T}Y-2Y^{T}X_{\gamma}\sqrt{\Lambda^{-1}}U^{T}\tilde{\beta}_{\gamma}+\tilde{\beta}_{\gamma}^{T}\tilde{\beta}_{\gamma}$
(163)
$\displaystyle\qquad=Y^{T}Y-2\tilde{Y}_{\gamma}^{T}\tilde{\beta}_{\gamma}+\tilde{\beta}_{\gamma}^{T}\tilde{\beta}_{\gamma}$
(164)
$\displaystyle\qquad=\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}|^{2}\big{)}+\big{(}\tilde{Y}_{\gamma}^{T}\tilde{Y}_{\gamma}-2\tilde{Y}_{\gamma}^{T}\tilde{\beta}_{\gamma}+\tilde{\beta}_{\gamma}^{T}\tilde{\beta}_{\gamma}\big{)}$
(165)
$\displaystyle\qquad=\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}|^{2}\big{)}+\|\tilde{Y}_{\gamma}-\tilde{\beta}_{\gamma}\|^{2}.$
(166)

Hence, we have

$\displaystyle d\beta_{\gamma}$
$\displaystyle=\det(U^{T}\Lambda^{-1/2})d\tilde{\beta}_{\gamma}$ (167)
$\displaystyle=\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}d\tilde{\beta}_{\gamma}.$ (168)

Hence, we have

$\displaystyle\int_{\beta_{\gamma}}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\|Y-X_{\gamma}\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi\tau^{-1}}\big{)}^{|\gamma|}}\exp\bigg{(}-\frac{\|\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}\tau^{-1}}\bigg{)}d\beta_{\gamma}$
(169)
$\displaystyle\qquad=\int_{\tilde{\beta}_{\gamma}}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}|^{2}\big{)}+\|\tilde{Y}_{\gamma}-\tilde{\beta}_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}$
$\displaystyle\qquad\qquad\times\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi\tau^{-1}}\big{)}^{|\gamma|}}\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}d\tilde{\beta}_{\gamma}$ (170)
$\displaystyle\qquad=\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\tau^{|\gamma|/2}\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}|^{2}\big{)}}{2\sigma_{\gamma}^{2}}\bigg{)}\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}.$ (171)

By combining (155) and (171), we obtain

$\displaystyle p(Y|\gamma_{i}=1,\gamma_{-i},X\big{)}$
$\displaystyle=\int_{\beta_{\gamma}}\int_{\sigma_{\gamma}^{2}=0}^{\infty}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\exp\bigg{(}-\frac{\|Y-X_{\gamma}\beta_{\gamma}\|^{2}}{2\sigma_{\gamma}^{2}}\bigg{)}p(\beta_{\gamma}|\gamma_{i}=1,\gamma_{-i})p(\sigma_{\gamma}^{2}|\gamma_{i}=1,\gamma_{-i})d\beta_{\gamma}d\sigma_{\gamma}^{2}$
(172)
$\displaystyle=\int_{\sigma_{\gamma}^{2}=0}^{\infty}\texttt{InvGamma}\bigg{(}\frac{1}{2}\nu_{0},\frac{1}{2}\nu_{0}\lambda_{0}\bigg{)}\frac{1}{\big{(}\sigma_{\gamma}\sqrt{2\pi}\big{)}^{N}}\tau^{|\gamma|/2}$
$\displaystyle\qquad\times\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}|^{2}\big{)}}{2\sigma_{\gamma}^{2}}\bigg{)}\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}d\sigma_{\gamma}^{2}$ (173)
$\displaystyle\qquad=\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}\tau^{|\gamma|/2}(2\pi)^{-N/2}\int_{\sigma_{\gamma}^{2}=0}^{\infty}\texttt{InvGamma}\bigg{(}\frac{1}{2}\nu_{0},\frac{1}{2}\nu_{0}\lambda_{0}\bigg{)}(\sigma_{\gamma}^{2})^{-N/2}$
$\displaystyle\qquad\qquad\times\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}\|^{2}\big{)}}{2\sigma_{\gamma}^{2}}\bigg{)}d\sigma_{\gamma}^{2}$
(174) $\displaystyle\qquad=\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}\tau^{|\gamma|/2}(2\pi)^{-N/2}$
$\displaystyle\qquad\times\int_{\sigma_{\gamma}^{2}=0}^{\infty}\frac{(1/2\lambda_{0}\nu_{0})^{1/2\nu_{0}}}{\Gamma(1/2\nu_{0})}(1/\sigma_{\gamma}^{2})^{1/2\nu_{0}+1}\exp\bigg{(}-1/2\nu_{0}\lambda_{0}/\sigma_{\gamma}^{2}\bigg{)}(\sigma_{\gamma}^{2})^{-N/2}$
$\displaystyle\qquad\qquad\times\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}\|^{2}\big{)}}{2\sigma_{\gamma}^{2}}\bigg{)}d\sigma_{\gamma}^{2}$
(175) $\displaystyle\qquad=\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}\tau^{|\gamma|/2}(2\pi)^{-N/2}\frac{(1/2\lambda_{0}\nu_{0})^{1/2\nu_{0}}}{\Gamma(1/2\nu_{0})}$
$\displaystyle\qquad\times\int_{\sigma_{\gamma}^{2}=0}^{\infty}(1/\sigma_{\gamma}^{2})^{1/2\nu_{0}+1+N/2}\exp\bigg{(}-\frac{\big{(}\|Y\|^{2}-\|\tilde{Y}_{\gamma}\|^{2}+\nu_{0}\lambda_{0}\big{)}}{2\sigma_{\gamma}^{2}}\bigg{)}d\sigma_{\gamma}^{2}$
(176) $\displaystyle\qquad=\det(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1/2}\tau^{|\gamma|/2}(2\pi)^{-N/2}\frac{(1/2\lambda_{0}\nu_{0})^{1/2\nu_{0}}}{\Gamma(1/2\nu_{0})}$
$\displaystyle\qquad\qquad\times\Gamma\bigg{(}\frac{N+\nu_{0}}{2}\bigg{)}\bigg{(}\frac{\|Y\|^{2}-\|\tilde{Y}_{\gamma}\|^{2}+\nu_{0}\lambda_{0}}{2}\bigg{)}^{-\frac{N+\nu_{0}}{2}}.$
(177)

Let $\tilde{\gamma}_{1}$ is given by $\gamma_{-i}$ with $\gamma_{i}=1$,
$\tilde{\gamma}_{0}$ is given by $\gamma_{-i}$ with $\gamma_{i}=0$. It follows
that

$\displaystyle\frac{p(Y|\gamma_{i}=1,\gamma_{-i},X)}{p(Y|\gamma_{i}=0,\gamma_{-i},X)}=\sqrt{\tau}\sqrt{\frac{\det(X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I)}{\det(X_{\tilde{\gamma}_{1}}^{T}X_{\tilde{\gamma}_{1}}+\tau
I)}}\bigg{(}\frac{\|Y\|^{2}-\|\tilde{Y}_{\tilde{\gamma}_{0}}\|^{2}+\nu_{0}\lambda_{0}}{\|Y\|^{2}-\|\tilde{Y}_{\tilde{\gamma}_{1}}\|^{2}+\nu_{0}\lambda_{0}}\bigg{)}^{\frac{N+\nu_{0}}{2}}.$
(178)

On the other hand, we have

$\displaystyle\|\tilde{Y}_{\gamma}\|^{2}$
$\displaystyle=\tilde{Y}_{\gamma}^{T}\tilde{Y}_{\gamma}$ (179)
$\displaystyle=Y^{T}X_{\gamma}(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1}X_{\gamma}^{T}Y.$ (180)

Hence, we finally have

$\displaystyle\frac{p(Y|\gamma_{i}=1,\gamma_{-i},X)}{p(Y|\gamma_{i}=0,\gamma_{-i},X)}=\sqrt{\tau\frac{\det(X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I)}{\det(X_{\tilde{\gamma}_{1}}^{T}X_{\tilde{\gamma}_{1}}+\tau
I)}\bigg{(}\frac{S_{\tilde{\gamma}_{0}}}{S_{\tilde{\gamma}_{1}}}\bigg{)}^{N+\nu_{0}}},$
(181)

where

$\displaystyle
S_{\gamma}:=Y^{T}Y-Y^{T}X_{\gamma}(X_{\gamma}^{T}X_{\gamma}+\tau
I)^{-1}X_{\gamma}^{T}Y+\nu_{0}\lambda_{0}.$ (182)

Based on this, we can estimate

$\displaystyle p(\gamma_{i}|\mathcal{D},\gamma_{-i})$
$\displaystyle=\frac{p(\gamma_{i}|\mathcal{D},\gamma_{-i})}{p(1-\gamma_{i}|\mathcal{D},\gamma_{-i})}\bigg{(}1+\frac{p(\gamma_{i}|\mathcal{D},\gamma_{-i})}{p(1-\gamma_{i}|\mathcal{D},\gamma_{-i})}\bigg{)}^{-1}.$
(183)

Denote the set of included variables in $\tilde{\gamma}_{0}$ as
$I=\\{j:\tilde{\gamma}_{0,j}=1\\}$ . Define
$F=\big{(}X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I\big{)}^{-1}$, $\nu=X^{T}Y$ and $\nu_{\tilde{\gamma}_{0}}=(\nu_{j})_{j\in
I}$. Also define $A=X^{T}X$ and $a_{i}=(A_{ji})_{j\in I}$. Then, by using the
same arguments as [23, Appendix B1], we can show that

$\displaystyle
S(\tilde{\gamma}_{1})=S(\tilde{\gamma}_{0})-d_{i}\big{(}\nu_{\tilde{\gamma}_{0}}^{T}Fa_{i}-\nu_{i}\big{)}^{2},$
(184)

where $d_{i}=(A_{ii}+\tau-a_{i}^{T}Fa_{i})^{-1}$. In addition, we can compute
$a_{i}^{T}Fa_{i}$ by using the Cholesky decomposition of $F=LL^{T}$ and

$\displaystyle a_{i}^{T}Fa_{i}$ $\displaystyle=\|a_{i}^{T}L\|^{2}$ (185)
$\displaystyle=\sum_{j\in I}(BL)_{ij}^{2},$ (186)

where $B$ is the $p\times|\gamma|$ matrix made of the columns of $A$
corresponding to variables included in $\gamma$.

In addition, we have

$\displaystyle X_{\tilde{\gamma}_{1}}^{T}X_{\tilde{\gamma}_{1}}+\tau
I=\begin{pmatrix}X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I&a_{i}\\\ a_{i}^{T}&A_{ii}+\tau\end{pmatrix}$ (187)

Hence, by using Schur’s formula for the determinant of block matrix, we are
easy to see that

$\displaystyle\frac{\det(X_{\tilde{\gamma}_{0}}^{T}X_{\tilde{\gamma}_{0}}+\tau
I)}{\det(X_{\tilde{\gamma}_{1}}^{T}X_{\tilde{\gamma}_{1}}+\tau I)}=d_{i}.$
(188)

Using this algorithm, if pre-computing $X^{T}X$ is not possible, the
computational complexity per MCMC iteration is
$O(N|\gamma|^{2}+|\gamma|^{3}+P|\gamma|^{2})$. Otherwise, if pre-computing
$X^{T}X$ is possible, the computational complexity per MCMC iteration is
$O(|\gamma|^{3}+P|\gamma|^{2})$.

## References

  * Andrieu et al. [2004] C. Andrieu, N. de Freitas, A. Doucet, and M. I. Jordan.  An introduction to MCMC for machine learning.  _Machine Learning_ , 50:5–43, 2004.
  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Bolstad [2010] W. M. Bolstad.  _Understanding Computational Bayesian Statistics_.  John Wiley, 2010.
  * Breiman [1960] L. Breiman.  The strong law of large numbers for a class of Markov chains.  _Annals of Mathematical Statistics_ , 31:801–803, 1960\.
  * Bugallo et al. [2007] M. F. Bugallo, S. Xu, and P. M. Djurić.  Performance comparison of EKF and particle filtering methods for maneuvering targets.  _Digit. Signal Process._ , 17:774–786, 2007.
  * Combes and Touati [2019] R. Combes and M. Touati.  Computationally efficient estimation of the spectral gap of a markov chain.  _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ , 3:1 – 21, 2019.
  * Diaconis and Saloff-Coste [1993] P. Diaconis and L. Saloff-Coste.  Comparison theorems for reversible markov chains.  _Annals of Applied Probability_ , 3:696–730, 1993.
  * Fitzgerald [2001] W. J. Fitzgerald.  Markov chain Monte Carlo methods with applications to signal processing.  _Signal Process._ , 81:3–18, 2001.
  * Gupta and Rawlings [2014] A. Gupta and J. B. Rawlings.  Comparison of parameter estimation methods in stochastic chemical kinetic models: Examples in systems biology.  _AIChE journal. American Institute of Chemical Engineers_ , 60 4:1253–1268, 2014.
  * Hesterberg [2002] T. Hesterberg.  Monte carlo strategies in scientific computing.  _Technometrics_ , 44:403 – 404, 2002.
  * Jankowiak [2022] M. Jankowiak.  Bayesian variable selection in a million dimensions.  _ArXiv_ , abs/2208.01180, 2022.
  * Kasim et al. [2019] M. F. Kasim, A. F. A. Bott, P. Tzeferacos, D. Q. Lamb, G. Gregori, and S. M. Vinko.  Retrieving fields from proton radiography without source profiles.  _Physical review. E_ , 100 3-1:033208, 2019.
  * Liang et al. [2010] F. Liang, C. Liu, and R. J. Carroll.  Advanced Markov chain Monte Carlo methods: Learning from past samples.  2010\.
  * Paulin [2015] D. Paulin.  Concentration inequalities for Markov chains by Marton couplings and spectral methods.  _Electronic Journal of Probability_ , 20(79):1 – 32, 2015.
  * Rao [2018] S. Rao.  A Hoeffding inequality for Markov chains.  _Electronic Communications in Probability_ , 2018.
  * Read et al. [2012] J. Read, L. Martino, and D. Luengo.  Efficient Monte Carlo methods for multi-dimensional learning with classifier chains.  _Pattern Recognit._ , 47:1535–1546, 2012.
  * Robert and Casella [2005] C. P. Robert and G. Casella.  Monte carlo statistical methods.  _Technometrics_ , 47:243 – 243, 2005.
  * Truong [2020] L. V. Truong.  Replica analysis of the linear model with Markov or hidden Markov signal priors.  _ArXiv_ , abs/2009.13370, 2020.
  * Truong [2021] L. V. Truong.  Linear models with hidden Markov sources via replica method.  _2021 IEEE International Symposium on Information Theory (ISIT)_ , pages 396–401, 2021.
  * Truong [2022] L. V. Truong.  On linear model with markov signal priors.  In _AISTATS_ , 2022.
  * Tuominen and Tweedie [1979] P. Tuominen and R. L. Tweedie.  Markov Chains with Continuous Components.  _Proceedings of the London Mathematical Society_ , s3-38(1):89–114, 01 1979.
  * Wolfer and Kontorovich [2019] G. Wolfer and A. Kontorovich.  Estimating the mixing time of ergodic Markov chains.  In _32nd Annual Conference on Learning Theory_ , 2019.
  * Zanella and Roberts [2019] G. Zanella and G. O. Roberts.  Scalable importance tempering and Bayesian variable selection.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 81, 2019.