1.87cm1.87cm1.87cm1.87cm capbtabboxtable[][]

# Lifted samplers for partially ordered discrete state-spaces

Philippe Gagnon 1, Florian Maire 1

###### Abstract

A technique called lifting is employed in practice for avoiding that the
Markov chains simulated for sampling backtrack too often. It consists in
lifting the state-space to include direction variables for guiding these
chains. Its implementation is direct when the probability mass function
targeted is defined on a totally ordered set, such as that of a univariate
random variable taking values on the integers. In this paper, we adapt this
technique to the situation where only a partial order can be established and
explore its benefits. Important applications include simulation of systems
formed from binary variables, such as those described by the Ising model, and
variable selection when the marginal model probabilities can be evaluated, up
to a normalising constant. To accommodate for the situation where one does not
have access to these marginal model probabilities, a lifted trans-dimensional
sampler for partially ordered model spaces is introduced. We show through
theoretical analyses and empirical experiments that the lifted samplers
outperform their non-lifted counterparts in some situations, and this at no
extra computational cost. The code to reproduce all experiments is available
online.111See the ArXiv page of this paper.

1Department of Mathematics and Statistics, Université de Montréal.

Keywords: Bayesian statistics; binary random variables; Ising model; Markov
chain Monte Carlo methods; Peskun-Tierney ordering; trans-dimensional
samplers; variable selection.

##  1 Introduction

###  1.1 Partially ordered state-spaces

A partially ordered set $\bm{\mathcal{X}}$ is such that there exists a
reflexive, antisymmetric, and transitive binary relation defined through a set
$\bm{\mathcal{R}}$ in which it is possible to establish that
$\mathbf{x}\leq\mathbf{y}$ and $\mathbf{y}\geq\mathbf{x}$ as a consequence of
$(\mathbf{x},\mathbf{y})\in\bm{\mathcal{R}}$ for
$\mathbf{x},\mathbf{y}\in\bm{\mathcal{X}}$. In such a set, pairs are
comparable when either $\mathbf{x}<\mathbf{y}$ or $\mathbf{y}<\mathbf{x}$, and
are incomparable when neither $\mathbf{x}<\mathbf{y}$ nor
$\mathbf{y}<\mathbf{x}$. This represents the difference with totally ordered
sets such as $\operatorname{\mathbb{R}}$ or $\mathbb{N}$ for which every pair
of different elements is comparable. We refer the reader to Trotter (1992) for
more details on partially ordered sets.

An important example of such sets is when any $\mathbf{x}\in\bm{\mathcal{X}}$
can be written as a vector $\mathbf{x}:=(x_{1},\ldots,x_{n})$ for which each
component $x_{i}$ can be of two types, say Type A or Type B, denoted by
$x_{i}\in\\{\text{A},\text{B}\\}$. In this case, $\bm{\mathcal{R}}$ can be set
to

$\bm{\mathcal{R}}:=\\{(\mathbf{x},\mathbf{y})\in\bm{\mathcal{X}}\times\bm{\mathcal{X}}:n_{\text{A}}(\mathbf{x})\leq
n_{\text{A}}(\mathbf{y})\\},$

where $n_{\text{A}}(\mathbf{x})$ is the number of components of Type A in
$\mathbf{x}$:
$n_{\text{A}}(\mathbf{x})=\sum_{i=1}^{n}\mathds{1}_{x_{i}=\text{A}}$. The
function $n_{\text{B}}(\mathbf{x})$ is defined analogously. Note that
$n=n_{\text{A}}(\mathbf{x})+n_{\text{B}}(\mathbf{x})$ for all $\mathbf{x}$ and
that the order can be symmetrically established by instead considering
$n_{\text{B}}(\mathbf{x})$ and $n_{\text{B}}(\mathbf{y})$.

Two main statistical problems fit within this framework: simulation of binary
random variables such as graphs or networks and variable selection. Indeed,
for the former, $\bm{\mathcal{X}}$ can be parameterized such that
$\bm{\mathcal{X}}=\\{-1,+1\\}^{n}$, where for example for an Ising model,
$x_{i}\in\\{-1,+1\\}$ represents the state of a spin, or, for networks,
$\bm{\mathcal{X}}=\mathcal{M}_{n}(\\{0,1\\})$ where $x_{i,j}\in\\{0,1\\}$
indicates whether nodes $i$ and $j$ are connected or not. For variable
selection, $\bm{\mathcal{X}}=\\{0,1\\}^{n}$ and $x_{i}\in\\{0,1\\}$ indicates
whether or not the $i$th covariate is included in the model characterised by
the vector $\mathbf{x}\in\bm{\mathcal{X}}$.

###  1.2 Sampling on partially ordered state-spaces

Let $\pi$ be a probability distribution defined on a measurable space
$(\bm{\mathcal{X}},\bm{\mathsf{X}})$ where $\bm{\mathcal{X}}$ is a partially
ordered set and $\bm{\mathsf{X}}$ a sigma algebra on $\bm{\mathcal{X}}$ and
consider the problem of sampling from $\pi$. We assume that it is not possible
to generate independent draws from $\pi$. Given the complex dependency
structure of modern statistical models such as graphical models and the
intractable nature of some distributions that arise, for instance, in Bayesian
statistics, this is a realistic assumption. We turn to Markov chains based
sampling methods which typically rely on an ergodic stochastic process
$\\{\mathbf{X}(m)\,:\,m\in\mathbb{N}\\}$ whose limiting distribution is $\pi$.
A typical Markov chain based sampler, such as the Glauber dynamics for
graphical models or the tie-no-tie sampler for network models, selects
uniformly at random one of the coordinates of $\mathbf{x}$, say $x_{i}$, and
proposes to flip its value from B to A (if $x_{i}=\text{A}$), and accept or
reject this move according to a prescribed probability that guarantees that
the Markov chain limiting distribution is $\pi$. Such moves are often rejected
when the mass concentrates on a manifold of the ambient space. Zanella (2019)
recently proposed a locally informed generic approach for which the
probability to select the $i$th coordinate depends on the relative mass of the
resulting proposal, i.e. $\pi(\mathbf{y})/\pi(\mathbf{x})$, aiming at
proposing less naive moves. Yet, the sampler is of reversible
Metropolis–Hastings (MH, (Metropolis et al., 1953; Hastings, 1970)) type,
implying that the chain may often go back to recently visited states. When
this is the case, the state-space is explored through a random walk, a process
exhibiting a diffusive behaviour.

The lifting technique, which can be traced back to Gustafson (1998) and even
to Horowitz (1991), aim at producing Markov chains which do not suffer from
such a behaviour. This is achieved by introducing a momentum variable
$\nu\in\\{-1,+1\\}$ which provides the process
$\\{\mathbf{X}(m)\,:\,m\in\mathbb{N}\\}$ with some memory of its past in order
to avoid backtracking. To remain Markovian, the process is thus enlarged to
$\\{(\mathbf{X},\nu)(m)\,:\,m\in\mathbb{N}\\}$ and the momentum variable is
flipped at random times according to a prescribed dynamic which guarantees
that, marginally, the process $\\{\mathbf{X}(m)\,:\,m\in\mathbb{N}\\}$ retains
its limiting distribution. Lifted Markov chains have been quantitatively
studied in Chen et al. (1999) and Diaconis et al. (2000) and have been shown
to reduce, sometimes dramatically, the mixing time of random walks on groups.
Over the past few years, there has been a renewed interest for lifted
techniques in the computational statistics community: in addition to speeding-
up the mixing time of random walk Markov chains they are also suspected to
reduce asymptotic variances of empirical averages of observables of interests,
see Andrieu and Livingstone (2019) for some precise results. We also refer to
Gagnon and Doucet (2019), Syed et al. (2019) and Neal (2020) for examples
where popular Markov chain Monte Carlo (MCMC) algorithms such as reversible
jump (Green, 1995), parallel tempering (Geyer, 1991) and slice sampling (Neal,
2003) have seen their performance improved in some situations by considering
their lifted version. Remarkably, those lifted samplers are implemented at no
additional computational cost over their non-lifted counterparts, and also
with no additional implementation difficulty.

All these successful applications of the lifting technique have been carried
out in contexts where $\bm{\mathcal{X}}$ is one-dimensional (Gustafson, 1998)
or exhibits a one-dimensional parameter which plays a central role in the
sampling scheme: the annealing parameter in Syed et al. (2019), the model
indicator reflecting the size/complexity of the nested models in Gagnon and
Doucet (2019), and the height of the level-lines
$\\{\pi(\mathbf{X}(m))\,:\,m\in\mathbb{N}\\}$ in Neal (2020). When such a one-
dimensional feature does not exist, there does not exist a straightforward way
of lifting the state-space without facing issues of reducibility or obtaining
inefficient samplers. A possibility, deemed as naive in Gustafson (1998), is
to lift each marginal component of the state-space and update them one at the
time. When the state-space is uncountable, it is possible to construct a
persistent walk by introducing bounces at random event times which change the
direction of propagation (see Vanetti et al. (2017)). However, when the state-
space is countable and partially ordered, such an approach is infeasible.

The objective of this paper is to present and analyse generic methods based on
the lifting technique to sample from a given probability mass function (PMF)
with a partially ordered countable support. In particular, we break free from
the requirement of having a one-dimensional parameter by exploiting the local
one-dimensional neighborhood structure induced by the partial order on
$\bm{\mathcal{X}}$: the neighbourhood of $\mathbf{x}$, denoted by
$\mathcal{N}(\mathbf{x})$, is separated into two parts where one comprises
states with $\mathbf{y}>\mathbf{x}$, denoted by
$\mathcal{N}_{+1}(\mathbf{x})$, and the other one comprises states with
$\mathbf{y}<\mathbf{x}$, denoted by $\mathcal{N}_{-1}(\mathbf{x})$
(considering that $\mathcal{N}(\mathbf{x})$ is only composed of states that
can be compared with $\mathbf{x}$). Looking for instance at the variable
selection example, the partial order is defined by mean of the model sizes, or
in other words, the number of covariates included in the models. If the
momentum is $\nu=+1$, the chain is forced to attempt visiting models with more
variables until a move is rejected, then $\nu$ is flipped to $\nu=-1$. As a
consequence, the momentum variable remains one-dimensional while the Markov
chain is often irreducible and efficiently explore the state-space. An
illustration showing the benefit of this approach is provided at Figure 1.
Again, we stress that the typical lifted sampler is implemented at no
additional computational cost over its non-lifted counterpart, and also with
no additional implementation difficulty.

$\begin{array}[]{cc}\textbf{Random walk behaviour}&\textbf{Persistent
movement}\cr\vspace{-1mm}\textbf{ESS = 0.12 per it.}&\textbf{ESS = 0.33 per
it.}\cr\includegraphics[width=173.44534pt]{Fig_1_a.pdf}&\includegraphics[width=173.44534pt]{Fig_1_b.pdf}\end{array}$

Figure 1: Trace plots for the statistic of number of covariates included in
the model for a MH sampler with a locally informed proposal distribution
(discussed in more details in Section 3.2) and its lifted counterpart, when
applied to solve a real variable selection problem (presented in Section 4.2);
ESS stands for effective sample size

###  1.3 Overview of our contributions

In this paper, we focus on the simulation of two-dimensional Ising models and
variable selection problems, without restricting ourselves to these examples
of applications when we present the samplers and analyse them. For these
examples, a generic sampler that we study corresponds to the discrete-time
version of a specific sampler independently developed in Power and Goldman
(2019), a paper in which the focus is rather on exploiting any structure of
$\bm{\mathcal{X}}$ identified by users. The structure identified here is, in a
sense, that $\bm{\mathcal{X}}$ exhibits a partial order. We consequently do
not claim originality for the samplers that will be presented. Our
contributions are the following:

  * •

statement of the sampling problem within the specific framework of partially
ordered discrete state-spaces (so that it becomes straightforward to implement
a sampler using the lifting technique in this framework);

  * •

identification of situations in which the lifted samplers are expected to
outperform their non-lifted counterparts, based on theoretical analyses and
numerical experiments;

  * •

introduction of a trans-dimensional lifted sampler useful, among others, for
variable selection when it is not possible to integrate out the parameters.

###  1.4 Organisation of the paper

The generic algorithm is first presented in Section 2. We next analyse in
Section 3 two important versions with uniform proposal distributions and
locally informed ones, allowing to establish that they can outperform their
non-lifted counterparts under some assumptions. In Section 4, we show the
difference in empirical performance for a Ising model (Section 4.1) and real
variable selection problem (Section 4.2). In Section 5, we consider that
$\bm{\mathcal{X}}$ is a model space and propose a lifted trans-dimensional
sampler allowing to simultaneously achieve model selection and parameter
estimation. The paper finishes in Section 6 with retrospective comments and
possible directions for future research.

##  2 Generic algorithm

The sampler that we present is a MCMC algorithm that generates proposals
belonging to a subset of $\mathcal{N}(\mathbf{x})$ chosen according to a
“direction” $\nu\in\\{-1,+1\\}$, when the current state is
$\mathbf{x}\in\bm{\mathcal{X}}$. In particular, the proposals belong to
$\mathcal{N}_{+1}(\mathbf{x}):=\\{\mathbf{y}\in\mathcal{N}(\mathbf{x}):\mathbf{y}>\mathbf{x}\\}\subseteq\mathcal{N}(\mathbf{x})$
when $\nu=+1$ or
$\mathcal{N}_{-1}(\mathbf{x}):=\\{\mathbf{y}\in\mathcal{N}(\mathbf{x}):\mathbf{y}<\mathbf{x}\\}\subseteq\mathcal{N}(\mathbf{x})$
when $\nu=-1$, where $\mathcal{N}_{-1}(\mathbf{x})$ and
$\mathcal{N}_{+1}(\mathbf{x})$ thus denote two subsets of
$\mathcal{N}(\mathbf{x})$ such that
$\mathcal{N}_{-1}(\mathbf{x})\cup\mathcal{N}_{+1}(\mathbf{x})=\mathcal{N}(\mathbf{x})$.
More formally, assuming that the Markov chain state is
$\mathbf{x}\in\bm{\mathcal{X}}$, the proposal distribution
$q_{\mathbf{x},\nu}$ has its support restricted to
$\mathcal{N}_{\nu}(\mathbf{x})$. There exist natural candidates for such
distributions, as will be explained in the following. This makes the
implementation of the proposed sampler almost straightforward; once the
neighbourhood structure has been specified. In our framework, the partial
ordering induces a natural neighbourhood structure.

The sampler is based on the well known technique of lifting: the state-space
$\bm{\mathcal{X}}$ is lifted to an extended state-space
$\bm{\mathcal{X}}\times\\{-1,+1\\}$ such that the marginal and the conditional
distributions of the direction variable $\nu$ is the uniform distribution on
$\\{-1,+1\\}$. The algorithm, which bares a strong resemblance with the guided
walk (Gustafson, 1998), is now presented in Algorithm 1.

Algorithm 1 A lifted sampler for partially ordered discrete state-spaces

  1. 1.

Generate $\mathbf{y}\sim q_{\mathbf{x},\nu}$ and $u\sim\mathcal{[}0,1]$.

  2. 2.

If

$\displaystyle
u\leq\alpha_{\nu}(\mathbf{x},\mathbf{y}):=1\wedge\frac{\pi(\mathbf{y})\,q_{\mathbf{y},-\nu}(\mathbf{x})}{\pi(\mathbf{x})\,q_{\mathbf{x},\nu}(\mathbf{y})},$
(1)

set the next state of the chain to $(\mathbf{y},\nu)$. Otherwise, set it to
$(\mathbf{x},-\nu)$.

  3. 3.

Go to Step 1.

If $\bm{\mathcal{X}}$ is finite, there exists a boundary, in the sense that
there may be no mass beyond a state $\mathbf{x}$ when the direction followed
is $\nu$. For instance in variable selection, the posterior probability of a
model with more covariates than the maximum number is zero. Algorithm 1 may
thus seem incomplete, in the sense that it is not explicitly specified what to
do on the boundary. We in fact consider that $q_{\mathbf{x},\nu}$ is defined
even on the boundary, and that it is defined to generate a point outside of
$\bm{\mathcal{X}}$. This point will be automatically rejected (because its
mass is 0) and the chain will remain at $\mathbf{x}$ and the direction will be
reversed. Note that this is a technical requirement. In practice, one can
simply skip Step 1 in this case and set the next state to $(\mathbf{x},-\nu)$.

It is possible to establish the correctness of the algorithm through that of a
more general version based on the lifted algorithm presented in Andrieu and
Livingstone (2019). Before presenting this more general version which has
interesting features, we introduce necessary notation. Let
$\rho_{\nu}:\bm{\mathcal{X}}\to[0,1]$, for $\nu\in\\{-1,+1\\}$, be a user-
defined function for which we require that:

$\displaystyle 0\leq\rho_{\nu}(\mathbf{x})\leq
1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}}),$ (2)
$\displaystyle\rho_{\nu}(\mathbf{x})-\rho_{-\nu}(\mathbf{x})=T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}}),$
(3)

where we have set, for all
$(\mathbf{x},\nu)\in\bm{\mathcal{X}}\times\\{-1,+1\\}$,

$\displaystyle
T_{\nu}(\mathbf{x},\bm{\mathcal{X}}):=\sum_{\mathbf{x}^{\prime}\in\bm{\mathcal{X}}}q_{\mathbf{x},\nu}(\mathbf{x}^{\prime})\,\alpha_{\nu}(\mathbf{x},\mathbf{x}^{\prime})=\sum_{\mathbf{x}^{\prime}\in\mathcal{N}_{\nu}(\mathbf{x})}q_{\mathbf{x},\nu}(\mathbf{x}^{\prime})\,\alpha_{\nu}(\mathbf{x},\mathbf{x}^{\prime}).$
(4)

These conditions make the algorithm valid and are thus considered satisfied in
the sequel. Finally, let $Q_{\mathbf{x},\nu}$ be a PMF such that
$Q_{\mathbf{x},\nu}(\mathbf{x}^{\prime})\propto
q_{\mathbf{x},\nu}(\mathbf{x}^{\prime})\,\alpha_{\nu}(\mathbf{x},\mathbf{x}^{\prime})$.
The more general algorithm is now presented in Algorithm 2.

Algorithm 2 A more general lifted sampler for partially ordered discrete
state-spaces

  1. 1.

Generate $u\sim\mathcal{U}[0,1]$.

    1. (i)

If $u\leq T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$, generate $\mathbf{y}\sim
Q_{\mathbf{x},\nu}$ and set the next state of the chain to $(\mathbf{y},\nu)$;

    2. (ii)

if $T_{\nu}(\mathbf{x},\bm{\mathcal{X}})<u\leq
T_{\nu}(\mathbf{x},\bm{\mathcal{X}})+\rho_{\nu}(\mathbf{x})$, set the next
state of the chain to $(\mathbf{x},-\nu)$;

    3. (iii)

if $u>T_{\nu}(\mathbf{x},\bm{\mathcal{X}})+\rho_{\nu}(\mathbf{x})$, set the
next state of the chain to $(\mathbf{x},\nu)$.

  2. 2.

Go to Step 1.

######  Proposition 1.

The transition kernel of the Markov chain
$\\{(\mathbf{X},\nu)(m):m\in\operatorname{\mathbb{N}}\\}$ simulated by
Algorithm 2 admits $\pi\otimes\mathcal{U}\\{-1,1\\}$ as invariant
distribution.

###### Proof.

See Section 7. ∎

It is interesting to notice that $T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$
represents the probability to accept a proposal when the current state is
$(\mathbf{x},\nu)$. In Algorithm 2, we thus first decide if we accept a
proposal or not, and if it is the case, in Step 1.(i), we randomly select the
value of the proposal (using the conditional distribution). It can be readily
checked that choices for $\rho_{\nu}$ include
$\rho_{\nu}(\mathbf{x})=1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$ and
$\rho_{\nu}(\mathbf{x})=\max\\{0,T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})\\}$.

If $\rho_{\nu}(\mathbf{x})=1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$, the
condition for Case (iii) of Step 1 is never satisfied, and the algorithm
either accepts the proposal and keeps the same direction, or the proposal is
rejected and the direction is reversed. In this case, one can show that
Algorithm 2 corresponds to Algorithm 1, which is what allows ensuring the
correctness of Algorithm 1 as well. Setting $\rho_{\nu}(\mathbf{x})$ otherwise
than $\rho_{\nu}(\mathbf{x})=1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$ allows in
Case (iii) of Step 1 to keep following the same direction, even when the
proposal is rejected. Intuitively, this is desirable when the rejection is due
to “bad luck”, and not because there is no mass in the direction followed. The
function $\rho_{\nu}(\mathbf{x})$ aims at incorporating this possibility in
the sampler.

In the typical MCMC framework, when one wants to sample from a probability
density function, it is not possible to directly compute
$T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$ as it requires computing an integral
with respect to this density function. In such a case, it is therefore usually
not possible to set $\rho_{\nu}(\mathbf{x})$ otherwise than
$1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$. This contrasts with our discrete
state-space framework in which it is often possible to directly compute
$T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$. A corollary of Theorem 3.15 in Andrieu
and Livingstone (2019) states that the best choice of function $\rho_{\nu}$ in
terms of asymptotic variance is

$\displaystyle\rho_{\nu}^{*}(\mathbf{x}):=\max(0,T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})),$
(5)

and that the worst choice is
$\rho_{\nu}^{\text{w}}(\mathbf{x}):=1-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$.

######  Corollary 1.

If $\bm{\mathcal{X}}$ is finite, then for any function $\rho_{\nu}$ and
$f\in\mathcal{L}_{2}^{*}(\bar{\pi}):=\left\\{f:\int
d\bar{\pi}f^{2}<\infty\text{ and }f(\mathbf{x},-1)=\right.$
$\left.f(\mathbf{x},+1)\text{ for all }\mathbf{x}\right\\}$,

$\mathrm{var}(f,P_{\rho^{*}})\leq\mathrm{var}(f,P_{\rho})\leq\mathrm{var}(f,P_{\rho^{\text{w}}}),$

where $\bar{\pi}:=\pi\otimes\mathcal{U}\\{-1,+1\\}$,
$\mathrm{var}(f,P_{\rho}):=\mathbb{V}\mathrm{ar}f(\mathbf{X},\nu)+2\sum_{k>0}\left\langle
f,P_{\rho}^{k}f\right\rangle$ and $P_{\rho}$ is the transition kernel
simulated by Algorithm 2, $\left\langle f,P_{\rho}^{k}f\right\rangle$ being
the inner product, i.e. $\left\langle f,P_{\rho}^{k}f\right\rangle:=\int
d\bar{\pi}fP_{\rho}^{k}f$.

###### Proof.

See Section 7. ∎

The price to pay for using $\rho_{\nu}^{*}$ instead of
$\rho_{\nu}^{\text{w}}$, for instance, is that the algorithm is more
complicated to implement because it is required to systematically compute
$T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$ at each iteration (it is also sometimes
required to compute $T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})$). Using
$\rho_{\nu}^{*}$ also comes with an additional computation cost (which is
discussed in Section 3.2). In our numerical experiments, it is seen that the
gain in efficiency of using Algorithm 2 with $\rho_{\nu}^{*}$ over Algorithm 1
is not significant. One may thus opt for simplicity and implement Algorithm 1.

##  3 Analysis of specific samplers

In the previous section, we presented generic algorithms with conditions
ensuring that they are valid. A necessary input to implement them is the
proposal distribution $q_{\mathbf{x},\nu}$ to be used. In this section, we
explore two natural choices and analyse their asymptotic variances. We start
in Section 3.1 with the common situation where the proposal is picked
uniformly at random. As mentioned in the introduction, this choice may lead to
poor mixing. We thus go on in Section 3.2 to a distribution incorporating
information about the target. This section finishes with a brief discussion in
Section 3.3 on the computational costs associated to these choices in regular
MH samplers and their lifted counterparts.

###  3.1 Uninformed uniform proposal

In the reversible MH sampler, it is common to set the proposal distribution,
denoted by $q_{\mathbf{x}}$ for this algorithm, to
$q_{\mathbf{x}}:=\mathcal{U}\\{\mathcal{N}(\mathbf{x})\\}$. In our framework,
the analogous proposal distribution is naturally defined as
$q_{\mathbf{x},\nu}:=\mathcal{U}\\{\mathcal{N}_{\nu}(\mathbf{x})\\}$. In this
case, the acceptance probability (1) of a proposal becomes

$\alpha_{\nu}(\mathbf{x},\mathbf{y})=1\wedge\frac{\pi(\mathbf{y})\,|\mathcal{N}_{\nu}(\mathbf{x})|}{\pi(\mathbf{x})\,|\mathcal{N}_{-\nu}(\mathbf{y})|}.$

For ease of presentation of the analysis, consider again the important example
described in Section 1.1 where each component $x_{i}$ of
$\mathbf{x}=(x_{1},\ldots,x_{n})$ can be of two types. We in this section
highlight the dependency on $n$ (the dimension) of the state-space and target
because it will be relevant in our analysis. We thus write $\pi_{n}$ for the
target and $\bm{\mathcal{X}}_{n}$ for the state-space, where each state is of
the form $\mathbf{x}:=(x_{1},\ldots,x_{n})$ with
$x_{i}\in\\{\text{A},\text{B}\\}$. For now on, consider for that $\text{A}=-1$
and $\text{B}=+1$, corresponding to the case of Ising model. We note that
there is no loss of generality of considering this special case within the
important example.

In a MH sampler, one sets
$\mathcal{N}(\mathbf{x}):=\\{\mathbf{y}\in\bm{\mathcal{X}}_{n}:\sum_{i}|x_{i}-y_{i}|=2\\}=\\{\mathbf{y}\in\bm{\mathcal{X}}_{n}:\exists
j\text{ such that }y_{j}=-x_{j}\\}$, so that the algorithm proposes to flip a
single bit at each iteration. It thus chooses uniformly at random which bit to
flip. Therefore, the size of the neighbourhoods in this sampler is constant
for any $\mathbf{x}$ and is given by $n$. This implies that the acceptance
probability in this sampler, denoted by $\alpha(\mathbf{x},\mathbf{y})$,
reduces to
$\alpha(\mathbf{x},\mathbf{y})=1\wedge\pi(\mathbf{y})/\pi(\mathbf{x})$.

In the lifted case, the acceptance probability can be rewritten as

$\displaystyle\alpha_{\nu}(\mathbf{x},\mathbf{y})=1\wedge\frac{\pi(\mathbf{y})\,n_{-\nu}(\mathbf{x})}{\pi(\mathbf{x})\,n_{\nu}(\mathbf{y})}.$
(6)

Indeed,
$\mathcal{N}_{\nu}(\mathbf{x}):=\\{\mathbf{y}\in\bm{\mathcal{X}}_{n}:\exists
j\text{ such that }y_{j}=-x_{j}=\nu\\}$, which implies that
$|\mathcal{N}_{\nu}(\mathbf{x})|=n_{-\nu}(\mathbf{x})$ (with the analogous
implication for $\mathcal{N}_{-\nu}(\mathbf{y})$). The acceptance probability
$\alpha_{\nu}$ thus depends on an additional term
$n_{-\nu}(\mathbf{x})/n_{\nu}(\mathbf{y})$ compared to that in the MH sampler.
This term may have an negative impact by decreasing the acceptance
probability. This represents in fact the price to pay for using the lifted
sampler Algorithm 2 (including Algorithm 1 as a special case): the reversible
sampler is allowed to backtrack, which makes the sizes of the neighbourhoods
constant, whereas it is the opposite for Algorithm 2. The size of the
neighbourhoods diminishes in the lifted sampler as the chain moves further in
a direction (making the neighbourhoods in the reverse direction bigger and
bigger).

The impact is alleviated when $n$ is large and the mass of $\pi_{n}$
concentrates in the interior of $\bm{\mathcal{X}}_{n}$ in an area where
$|n_{-\nu}(\mathbf{x})-n/2|\leq\kappa$, where $\kappa$ is a positive integer.
Indeed, $n_{\nu}(\mathbf{y})=n_{\nu}(\mathbf{x})+1=n-n_{-\nu}(\mathbf{x})+1$,
which implies that $\alpha_{\nu}(\mathbf{x},\mathbf{y})\approx
1\wedge\pi(\mathbf{y})/\pi(\mathbf{x})=\alpha(\mathbf{x},\mathbf{y})$. In
fact, in an ideal situation where
$|\mathcal{N}_{-1}(\mathbf{x})|=|\mathcal{N}_{+1}(\mathbf{x})|=|\mathcal{N}(\mathbf{x})|/2$
for all $\mathbf{x}\in\bm{\mathcal{X}}_{n}$ (which is impossible for the
important example with types A and B), a corollary of Theorem 3.17 in Andrieu
and Livingstone (2019) establishes that Algorithm 2 with any function
$\rho_{\nu}$ dominates the MH algorithm in terms of asymptotic variances. In
particular, Algorithm 1 dominates the MH algorithm. Before presenting this
corollary, we define the transition kernel simulated by Algorithm 2 with
$q_{\mathbf{x},\nu}:=\mathcal{U}\\{\mathcal{N}_{\nu}(\mathbf{x})\\},\mathbf{x}\in\bm{\mathcal{X}}_{n},\nu\in\\{-1,+1\\}$,
as $P_{\rho,n}$ and that simulated by the MH sampler with
$q_{\mathbf{x}}:=\mathcal{U}\\{\mathcal{N}(\mathbf{x})\\},\mathbf{x}\in\bm{\mathcal{X}}_{n}$,
as $P_{\text{MH},n}$.

######  Corollary 2.

If

(a)

$\bm{\mathcal{X}}_{n}$ is finite,

(b)

$|\mathcal{N}_{-1}(\mathbf{x})|=|\mathcal{N}_{+1}(\mathbf{x})|=|\mathcal{N}(\mathbf{x})|/2=n^{*}/2$
for all $\mathbf{x}\in\bm{\mathcal{X}}_{n}$, where $n^{*}$ is a positive
integer that does not depend on $\mathbf{x}$ (but that may depend on $n$),

then for any $f\in\mathcal{L}_{2}^{*}(\bar{\pi})$ and $n$,

$\mathrm{var}(f,P_{\rho,n})\leq\mathrm{var}(f,P_{\text{MH},n}).$

###### Proof.

See Section 7. ∎

In light of the above, one might expect the inequality to (approximately) hold
(up to an error term) when the mass is highly concentrated on the points
$\mathbf{x}$ that are not too far from the center of the domain, where here
the notion of centrality is defined in terms of the distance between
$n_{-1}(\mathbf{x})$ or $n_{+1}(\mathbf{x})$ to $n/2$, suggesting that the
lifted sampler outperforms the reversible MH algorithm in this situation. The
rest of the section is dedicated to the introduction of conditions under which
this statement is true.

The key argument in proving Corollary 2 is to show that

$\displaystyle
q_{\mathbf{x}}(\mathbf{y})\,\alpha(\mathbf{x},\mathbf{y})=(1/2)\,q_{\mathbf{x},+1}(\mathbf{y})\,\alpha_{+1}(\mathbf{x},\mathbf{y})+(1/2)\,q_{\mathbf{x},-1}(\mathbf{y})\,\alpha_{-1}(\mathbf{x},\mathbf{y}),$
(7)

for all $\mathbf{x}$ and $\mathbf{y}$. Indeed, once this is done, Theorem 3.17
in Andrieu and Livingstone (2019) can be directly applied, which yields the
result. This in fact implies that once one has designed a lifted sampler, it
is possible to identify its non-lifted counterpart through (7), and to
establish that the latter is inferior. In particular, we can establish that
the sampler that flips a coin at each iteration to next decide which PMF to
use between $q_{\mathbf{x},+1}$ and $q_{\mathbf{x},-1}$ to generate a proposal
$\mathbf{y}$ that will be subject to approval using
$\alpha_{+1}(\mathbf{x},\mathbf{y})$ or $\alpha_{-1}(\mathbf{x},\mathbf{y})$
is inferior. Denote by $P_{\text{rev.},n}$ the Markov kernel simulated by this
algorithm. Now, what we would ideally do is to show a Peskun-Tierney ordering
(Peskun, 1973; Tierney, 1998) between $P_{\text{rev.},n}$ and
$P_{\text{MH},n}$ to establish the domination of the lifted sampler over the
reversible MH algorithm. Such an ordering is difficult to obtain as one needs
to show that for any pair $(\mathbf{x},\mathbf{y})$ such that
$\mathbf{x}\neq\mathbf{y}$, $P_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq
P_{\text{MH},n}(\mathbf{x},\mathbf{y})$.

A more general ordering is presented in Zanella (2019): if
$P_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq\omega
P_{\text{MH},n}(\mathbf{x},\mathbf{y})$ for all $\mathbf{x}\neq\mathbf{y}$,
then
$\mathrm{var}(f,P_{\text{rev.},n})\leq\mathrm{var}(f,P_{\text{MH},n})/\omega+((1-\omega)/\omega)\mathbb{V}\mathrm{ar}f(\mathbf{X})$,
where $\omega$ is a positive constant. Note that $f$ is a function of $\nu$ as
well, but because of the restriction $f(\mathbf{x},-1)=f(\mathbf{x},+1)$,
$\nu$ can be treated as a constant. We show in this section that if
$P_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq\omega
P_{\text{MH},n}(\mathbf{x},\mathbf{y})$ with $\omega$ close to 1 for all
$\mathbf{x}\neq\mathbf{y}$ belonging to a specific set having a high
probability, then
$\mathrm{var}(f,P_{\text{rev.},n})\leq\mathrm{var}(f,P_{\text{MH},n})+\text{small
error term}$ (under regularity conditions).

We start by defining this set:

$\displaystyle\bm{\mathcal{X}}_{\varphi(n)}:=\\{\mathbf{x}\in\bm{\mathcal{X}}_{n}:n/2-\varphi(n)\leq
n_{-1}(\mathbf{x}),n_{+1}(\mathbf{x})\leq n/2+\varphi(n)\\},$ (8)

where $\varphi$ is a function such that $0\leq\varphi(n)<n/2$ and, for
$\epsilon>0$, there exists $N>0$ such that for all $n\geq N$, we have that
$\varphi(n)/n<\epsilon$. For $\mathbf{x}\neq\mathbf{y}$ belonging to
$\bm{\mathcal{X}}_{\varphi(n)}$, if we consider $\omega_{n}$ now a function of
$n$, it is possible to establish that: for $\epsilon>0$, there exists a $N>0$
such that for all $n\geq N$,
$P_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq\omega_{n}P_{\text{MH},n}(\mathbf{x},\mathbf{y})$
with $1-\omega_{n}<\epsilon$. The explicit form of $\omega_{n}$ is

$\displaystyle\omega_{n}:=\left(1+\frac{\varphi(n)}{n/2}\right)^{-2}\left(1-\frac{\varphi(n)}{n/2}\right).$
(9)

One can imagine that if: the mass concentrates on
$\bm{\mathcal{X}}_{\varphi(n)}$, the chains do not often leave this set, and
when they do they do not take too much time to come back, then the asymptotic
variances should not be too different to those of chains with stationary
distribution $\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$, which assigns null
mass on $\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}$ and is thus the
normalised version of $\pi_{n}$ on $\bm{\mathcal{X}}_{\varphi(n)}$:

$\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}(\mathbf{x}):=\begin{cases}\pi_{n}(\mathbf{x})\big{/}\sum_{\mathbf{x}^{\prime}\in\bm{\mathcal{X}}_{\varphi(n)}}\pi_{n}(\mathbf{x}^{\prime})&\text{if}\quad\mathbf{x}\in\bm{\mathcal{X}}_{\varphi(n)},\cr
0&\text{if}\quad\mathbf{x}\notin\bm{\mathcal{X}}_{\varphi(n)}.\end{cases}$

This is what we obtain assuming such a behaviour for the chains generated by
$P_{\text{rev.},n}$ and $P_{\text{MH},n}$. We also require that the chains
generated by the Markov kernels with the same proposal mechanisms as
$P_{\text{rev.},n}$ and $P_{\text{MH},n}$, but whose stationary distribution
is $\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$, mix sufficiently well (in a
sense that will be made precise). We denote by $\tilde{P}_{\text{rev.},n}$ and
$\tilde{P}_{\text{MH},n}$ these Markov kernels.

######  Theorem 1.

Pick $f\in\mathcal{L}_{2}^{*}(\bar{\pi})$ and consider that all states
$\mathbf{x}\in\bm{\mathcal{X}}_{n}$ are such that
$\mathbf{x}=(x_{1},\ldots,x_{n})$ with $x_{i}\in\\{-1,+1\\}$. If, for
$\epsilon>0$, there exists $N>0$ such that for all $n\geq N$ it is possible to
choose a constant that may depend on $n$, $\varrho(n)$, with

(a)

$\sum_{k=\varrho(n)+1}^{\infty}\left\langle f,P^{k}f\right\rangle<\epsilon$,
for all
$P\in\\{P_{\text{rev.},n},\tilde{P}_{\text{rev.},n},P_{\text{MH},n},\tilde{P}_{\text{MH},n}\\}$,

(b)

$\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cup_{m=0}^{k}A_{m}^{\mathsf{c}}(\bm{\mathcal{X}}_{\varphi(n)-1})}]<\epsilon$,
for all $P\in\\{P_{\text{rev.},n},P_{\text{MH},n}\\}$, where it is considered
that the chain starts at stationarity and evolves using $P$ (and
$\mathbb{E}[f(\mathbf{X}(k))]=0$ without loss of generality) and
$A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1}):=\\{\mathbf{X}(m)\in\bm{\mathcal{X}}_{\varphi(n)-1}\\}$
(see (8) for the definitions of $\bm{\mathcal{X}}_{\varphi(n)-1}$ and
$\varphi(n)$),

(c)

$(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cap_{m=0}^{k}A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1})}]<\epsilon$,
for all $P\in\\{P_{\text{rev.},n},P_{\text{MH},n}\\}$,

(d)

$(1/\omega_{n}-1)\mathrm{var}(f,\tilde{P}_{\text{MH},n})<\epsilon$ and
$((1-\omega_{n})/\omega_{n})\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}f(\mathbf{X})<\epsilon$,
where $\omega_{n}$ is defined in (9) and
$\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}$ denotes
a variance computed with respect to
$\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$,

then there exists a positive constant $\kappa$ such that

$\mathrm{var}(f,P_{\rho,n})\leq\mathrm{var}(f,P_{\text{rev.},n})\leq\mathrm{var}(f,P_{\text{MH},n})+\kappa\epsilon.$

###### Proof.

See Section 7. ∎

There are several assumptions involved in Theorem 1. But, to put this into
perspective, Assumptions (c) and (d) are automatically verified if
$\mathrm{var}(f,\tilde{P}_{\text{MH},n})$ and
$\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}f(\mathbf{X})$
are bounded by constants that do not depend on $n$. Assumptions (a) and (b)
are really the crucial ones.

###  3.2 Locally informed proposal

In this section, we discuss and analyse the use of locally-balanced proposal
distributions, as defined by Zanella (2019) in the reversible MH framework as

$\displaystyle
q_{\mathbf{x}}(\mathbf{y}):=g\left(\frac{\pi_{n}(\mathbf{y})}{\pi_{n}(\mathbf{x})}\right)\bigg{/}c_{n}(\mathbf{x}),\quad\mathbf{y}\in\mathcal{N}(\mathbf{x}),$
(10)

where $c_{n}(\mathbf{x})$ represents the normalising constant, i.e.
$c_{n}(\mathbf{x}):=\sum_{\mathbf{x}^{\prime}\in\mathcal{N}(\mathbf{x})}g\left(\pi_{n}(\mathbf{x}^{\prime})/\pi_{n}(\mathbf{x})\right)$,
and $g$ is a strictly positive continuous function such that $g(x)/g(1/x)=x$.
Note that we used the same notation as in Section 3.1 for the proposal
distribution; this is to simplify. We will use the same notation as in Section
3.1 for the proposal distribution of the lifted sampler and for the Markov
kernels as well. In this section, it will be implicit that the proposal
distributions are informed proposals and the Markov kernels are those induced
by these informed proposals. Note also that we again highlight the
dependencies on $n$ of some terms that will appear in our analysis.

Such a function $g$ defined in (10) implies that the acceptance probability in
the MH algorithm is given by

$\alpha(\mathbf{x},\mathbf{y})=1\wedge\frac{\pi_{n}(\mathbf{y})\,q_{\mathbf{y}}(\mathbf{x})}{\pi_{n}(\mathbf{x})\,q_{\mathbf{x}}(\mathbf{y})}=1\wedge\frac{c_{n}(\mathbf{x})}{c_{n}(\mathbf{y})}.$

Zanella (2019) shows that $c_{n}(\mathbf{x})/c_{n}(\mathbf{y})\longrightarrow
1$ as $n\longrightarrow\infty$ under some assumptions. More precisely, the
author defines $\mathbf{x}:=(x_{1},\ldots,x_{n})$ and considers that at any
given iteration, only a small fraction of the $n$ components is proposed to
change values. The result holds when the random variables
$(X_{1},\ldots,X_{n})$ exhibit a structure of conditional independence, which
implies that the normalising constants $c_{n}(\mathbf{x})$ and
$c_{n}(\mathbf{y})$ share a lot of terms. This is again a consequence of the
backtracking of the reversible sampler and is thus in contrast with what we
observe for the lifted algorithm.

Two choices for $g$ are $g(x)=\sqrt{x}$ and $g(x)=x/(1+x)$, the latter being
called the Barker proposal in reference to Barker (1965)’s acceptance
probability choice. The advantage of the latter choice is that it is a bounded
function of $x$, which stabilises the normalising constants and thus the
acceptance probability. This is shown in Zanella (2019) and in the continuous
random variable case in Livingstone and Zanella (2019). We use this function
in our numerical analyses.

The proposal distribution in the lifted sampler is given by

$q_{\mathbf{x,\nu}}(\mathbf{y}):=g\left(\frac{\pi_{n}(\mathbf{y})}{\pi_{n}(\mathbf{x})}\right)\bigg{/}c_{n,\nu}(\mathbf{x}),\quad\mathbf{y}\in\mathcal{N}_{\nu}(\mathbf{x}),$

where $c_{n,\nu}(\mathbf{x})$ is the normalising constant. In this case,

$\displaystyle\alpha_{\nu}(\mathbf{x},\mathbf{y}):=1\wedge\frac{\pi_{n}(\mathbf{y})\,q_{\mathbf{y},-\nu}(\mathbf{x})}{\pi_{n}(\mathbf{x})\,q_{\mathbf{x},\nu}(\mathbf{y})}=1\wedge\frac{c_{n,\nu}(\mathbf{x})}{c_{n,-\nu}(\mathbf{y})}.$
(11)

There are two main differences with the reversible sampler. Firstly, the
normalising constants $c_{n,\nu}(\mathbf{x})$ and $c_{n,-\nu}(\mathbf{y})$ are
sums with (in general) not the same number of terms. Consider again, for ease
of presentation of the analysis, the important example described in Section
1.1 where each component $x_{i}$ of $\mathbf{x}=(x_{1},\ldots,x_{n})$ can be
of two types and more specifically the special case of Ising model. We know
that in this case $c_{n,\nu}(\mathbf{x})$ is a sum of $n_{-\nu}(\mathbf{x})$
terms (see (6)), whereas $c_{n,-\nu}(\mathbf{y})$ is a sum of
$n_{\nu}(\mathbf{y})$ terms. The second main difference is that, in the MH
sampler, it is proposed to flip one of the $n_{-1}(\mathbf{x})$ components to
$+1$ or one of the $n_{+1}(\mathbf{x})$ components to $-1$, and
$c_{n}(\mathbf{x})$ is formed from these proposals. The constant
$c_{n}(\mathbf{y})$ is also formed from proposals to flip components to $+1$
or $-1$ with $n_{-1}(\mathbf{y})$ and $n_{+1}(\mathbf{y})$ close to
$n_{-1}(\mathbf{x})$ and $n_{+1}(\mathbf{x})$ (and this is why the ratio of
the two constants converges to 1 provided that there exists a structure of
conditional independence). In contrast, $c_{n,\nu}(\mathbf{x})$ is formed from
proposals to flip one of the $n_{-\nu}(\mathbf{x})$ components to $\nu$ and
$c_{n,-\nu}(\mathbf{y})$ from proposals to flip one of the
$n_{\nu}(\mathbf{y})=n_{\nu}(\mathbf{x})+1$ components to $-\nu$; the
compositions of these constants are thus fundamentally opposite. There is
therefore no guarantee that
$c_{n,\nu}(\mathbf{x})/c_{n,-\nu}(\mathbf{y})\longrightarrow 1$ even under the
conditions stated in Zanella (2019).

Nevertheless, there exists as in the previous section an ideal situation in
which the lifter sampler outperforms the reversible MH algorithm.

######  Corollary 3.

If

(a)

$\bm{\mathcal{X}}_{n}$ is finite,

(b)

$c_{n,-1}(\mathbf{x})=c_{n,+1}(\mathbf{x})=c_{n}(\mathbf{x})/2=c_{n}^{*}/2$
for all $\mathbf{x}\in\bm{\mathcal{X}}_{n}$, where $c_{n}^{*}$ is a positive
constant that does not depend on $\mathbf{x}$ (but that may depend on $n$),

then for any $f\in\mathcal{L}_{2}^{*}(\bar{\pi})$ and $n$,

$\mathrm{var}(f,P_{\rho,n})\leq\mathrm{var}(f,P_{\text{MH},n}).$

###### Proof.

Analogous to that of Corollary 2. ∎

A sufficient condition for Assumption (b) to be verified is:
$|\mathcal{N}_{-1}(\mathbf{x})|=|\mathcal{N}_{+1}(\mathbf{x})|=|\mathcal{N}(\mathbf{x})|/2=n^{*}/2$
(Assumption (b) in Corollary 2) and

$\frac{1}{n^{*}/2}\sum_{\mathbf{x}^{\prime}\in\mathcal{N}_{-1}(\mathbf{x})}g\left(\frac{\pi_{n}(\mathbf{x}^{\prime})}{\pi_{n}(\mathbf{x})}\right)=\frac{1}{n^{*}/2}\sum_{\mathbf{x}^{\prime}\in\mathcal{N}_{+1}(\mathbf{x})}g\left(\frac{\pi_{n}(\mathbf{x}^{\prime})}{\pi_{n}(\mathbf{x})}\right)=\frac{1}{n^{*}}\sum_{\mathbf{x}^{\prime}\in\mathcal{N}(\mathbf{x})}g\left(\frac{\pi_{n}(\mathbf{x}^{\prime})}{\pi_{n}(\mathbf{x})}\right)=\mu,$

for all $\mathbf{x}\in\bm{\mathcal{X}}_{n}$, where $\mu$ is a positive
constant. We thus notice that the acceptance probabilities can be directly
rewritten in terms of averages when
$|\mathcal{N}_{-1}(\mathbf{x})|=|\mathcal{N}_{+1}(\mathbf{x})|=|\mathcal{N}(\mathbf{x})|/2=n^{*}/2$
and that an additional condition to Assumption (b) in Corollary 2 is
sufficient in the locally informed case for ordering the asymptotic variances.
This thus allows establishing a connection with the uniform case.

As in the previous section, it is possible to derive conditions under which
the inequality in Corollary 3 holds approximately. They are based as before on
the definition of a set, which in this case involves states $\mathbf{x}$ that
are such that $c_{n,-1}(\mathbf{x})$ and $c_{n,+1}(\mathbf{x})$ are close to
$c_{n}^{*}/2$. These states do not have to be in an area such that
$n_{-1}(\mathbf{x})$ and $n_{+1}(\mathbf{x})$ are close to $n/2$, but in
return the mass have to be (in some sense) evenly spread out in the area to
which they belong.

We now define this set:

$\displaystyle\bm{\mathcal{X}}_{\varphi(n)}:=\\{\mathbf{x}\in\bm{\mathcal{X}}_{n}:c_{n}^{*}/2-\varphi(n)\leq
c_{n,-1}(\mathbf{x}),c_{n,+1}(\mathbf{x}),c_{n}(\mathbf{x})/2\leq
c_{n}^{*}/2+\varphi(n)\\},$ (12)

where $\varphi$ is in this section a function such that
$0\leq\varphi(n)<c_{n}^{*}/2$ and, for $\epsilon>0$, there exists $N>0$ such
that for all $n\geq N$, we have that $\varphi(n)/c_{n}^{*}<\epsilon$. We
consider to simplify that for any
$\mathbf{x},\mathbf{y}\in\bm{\mathcal{X}}_{\varphi(n)}$ there exists a
probable path from $\mathbf{x}$ to $\mathbf{y}$ generated by $P_{\text{MH},n}$
(and marginally for $P_{\rho,n}$) with all intermediate states belonging to
$\bm{\mathcal{X}}_{\varphi(n)}$ as well. We now define a restricted version of
$\bm{\mathcal{X}}_{\varphi(n)}$ for which from any state
$\mathbf{x}\in\bm{\mathcal{X}}_{\varphi(n)}$, all the possible proposals
$\mathbf{y}$ belong to $\bm{\mathcal{X}}_{\varphi(n)}$ as well; denote this
set by $\bm{\mathcal{X}}_{\varphi(n)}^{0}$. For $\mathbf{x}\neq\mathbf{y}$
belonging to $\bm{\mathcal{X}}_{\varphi(n)}$, it is possible to establish
that: for $\epsilon>0$, there exists a $N>0$ such that for all $n\geq N$,
$P_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq\omega_{n}P_{\text{MH},n}(\mathbf{x},\mathbf{y})$
with $1-\omega_{n}<\epsilon$. The explicit form of $\omega_{n}$ is

$\displaystyle\omega_{n}:=\left(1+\frac{\varphi(n)}{c_{n}^{*}/2}\right)^{-3}\left(1-\frac{\varphi(n)}{c_{n}^{*}/2}\right)^{3}.$
(13)

We are now ready to present the analogous result to Theorem 1, in which here
$\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$ is the normalised version of
$\pi_{n}$ on $\bm{\mathcal{X}}_{\varphi(n)}$ defined in (12).

######  Theorem 2.

Pick $f\in\mathcal{L}_{2}^{*}(\bar{\pi})$ and consider that all states
$\mathbf{x}\in\bm{\mathcal{X}}_{n}$ are such that
$\mathbf{x}=(x_{1},\ldots,x_{n})$ with $x_{i}\in\\{-1,+1\\}$. If, for
$\epsilon>0$, there exists $N>0$ such that for all $n\geq N$ it is possible to
choose a constant that may depend on $n$, $\varrho(n)$, with

(a)

$\sum_{k=\varrho(n)+1}^{\infty}\left\langle f,P^{k}f\right\rangle<\epsilon$,
for all
$P\in\\{P_{\text{rev.},n},\tilde{P}_{\text{rev.},n},P_{\text{MH},n},\tilde{P}_{\text{MH},n}\\}$,

(b)

$\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cup_{m=0}^{k}A_{m}^{\mathsf{c}}(\bm{\mathcal{X}}_{\varphi(n)}^{0})}]<\epsilon$,
for all $P\in\\{P_{\text{rev.},n},P_{\text{MH},n}\\}$, where it is considered
that the chain starts at stationarity and evolves using $P$ (and
$\mathbb{E}[f(\mathbf{X}(k))]=0$ without loss of generality) and
$A_{m}(\bm{\mathcal{X}}_{\varphi(n)}^{0}):=\\{\mathbf{X}(m)\in\bm{\mathcal{X}}_{\varphi(n)}^{0}\\}$,

(c)

$(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cap_{m=0}^{k}A_{m}(\bm{\mathcal{X}}_{\varphi(n)}^{0})}]<\epsilon$,
for all $P\in\\{P_{\text{rev.},n},P_{\text{MH},n}\\}$,

(d)

$(1/\omega_{n}-1)\mathrm{var}(f,\tilde{P}_{\text{MH},n})<\epsilon$ and
$((1-\omega_{n})/\omega_{n})\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}f(\mathbf{X})<\epsilon$,
where $\omega_{n}$ is defined in (13),

then there exists a positive constant $\kappa$ such that

$\mathrm{var}(f,P_{\rho,n})\leq\mathrm{var}(f,P_{\text{rev.},n})\leq\mathrm{var}(f,P_{\text{MH},n})+\kappa\epsilon.$

###### Proof.

Analogous to that of Theorem 1. ∎

It is possible to establish a connection with Theorem 1 as we did for
Corollary 3 with Corollary 2. Consider indeed that the set
$\bm{\mathcal{X}}_{\varphi(n)}$ is comprised of states $\mathbf{x}$ with
$n_{-1}(\mathbf{x})$ and $n_{+1}(\mathbf{x})$ close to $n/2$ and
$(1/n_{+1}(\mathbf{x}))\sum_{\mathbf{x}^{\prime}\in\mathcal{N}_{-1}(\mathbf{x})}g(\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x}))$,
$(1/n_{-1}(\mathbf{x}))\sum_{\mathbf{x}^{\prime}\in\mathcal{N}_{+1}(\mathbf{x})}g(\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x}))$
and
$(1/n)\sum_{\mathbf{x}^{\prime}\in\mathcal{N}(\mathbf{x})}g(\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x}))$
all close to a positive constant $\mu$. We thus notice that under a precise
sense of what we mean by close to, this special case fits within the
definition of $\bm{\mathcal{X}}_{\varphi(n)}$ (12), under an additional
condition compared to the definition of the set in the previous section (8).

###  3.3 Computational costs

We provide in this section an overview of the computational costs associated
to using the proposal distributions described in Sections 3.1 and 3.2. The
uniform distribution is the least expensive: at each iteration, one has to
generate from a uniform and then evaluate the acceptance probability which
requires the computation of a ratio $\pi(\mathbf{y})/\pi(\mathbf{x})$.
Consider that the cost of the latter is the important one, in the sense that
all the other costs are comparatively negligible. The approach of Zanella
(2019) thus costs twice as much, if we assume that the cost of computing any
ratio is the same and that the ratios
$\pi(\mathbf{x}^{\prime})/\pi(\mathbf{x})$ for all $\mathbf{x}^{\prime}$ in
the neighbourhood are all computed in parallel. Indeed, these ratios are
necessary to generate the proposal $\mathbf{y}$, but once the latter has been
generated, the process has to be repeated for the reverse move. This is true
for the reversible MH sampler and Algorithm 1. Therefore, if the informed
proposal leads to a sampler at least twice as effective (in terms of ESS for
instance), then it is beneficial. It is the case in all our numerical
experiments. Note that in light of the above, implementing the reversible MH
sampler or Algorithm 1 costs essentially the same.

Algorithm 2 is more costly. When used with a uniform distribution and
$\rho:=\rho^{*}$ (5), at each iteration, $|\mathcal{N}_{\nu}(\mathbf{x})|$
ratio evaluations are required to compute
$T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$ (4), and it is afterwards required to
compute $T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})$ from time to time. This makes
the implementation cost somewhere in between that of Algorithm 1 with a
uniform and Algorithm 1 using the approach of Zanella (2019). When Algorithm 2
is used with an informed proposal and $\rho:=\rho^{*}$ (5), the cost may
explode. Consider that parallel computing is used to compute any normalising
constant $c_{\nu}(\mathbf{x})$, but that the normalising constants are
computed sequentially, then at each iteration it is required to compute at
least $1+|\mathcal{N}_{\nu}(\mathbf{x})|$ normalising constants compared with
two in Algorithm 1. The computation time per iteration will thus roughly be at
least $(1+|\mathcal{N}_{\nu}(\mathbf{x})|)/2$ times larger.

##  4 Numerical experiments

We first consider in Section 4.1 simulation of an Ising model and use this as
a toy example for which we can control the dimension, the roughness of the
target and where the mass concentrates to show how the performance of the
lifted and non-lifted samplers varies when these parameters change. In Section
4.2, we evaluate their performance when employed to solve a real variable
selection problem.

###  4.1 Ising model

For the two-dimensional model, the ambiant space $(V_{\eta},E_{\eta})$ is a
$\eta\times\eta$ square lattice regarded here as a square matrix in which each
element takes either the value $-1$ or $+1$. We write each state as a vector
as before: $\mathbf{x}=(x_{1},\ldots,x_{n})$, where $n=\eta^{2}$. The states
can be encoded as for instance: the values of the components on the first line
are $x_{1},\ldots,x_{\eta}$, those on the second line
$x_{\eta+1},\ldots,x_{2\eta}$, and so on. The PMF is given by

$\pi(\mathbf{x})=\frac{1}{Z}\exp\left(\sum_{i}\alpha_{i}x_{i}+\lambda\sum_{\langle
ij\rangle}x_{i}x_{j}\right),$

where $\alpha_{i}\in\operatorname{\mathbb{R}}$ and $\lambda>0$ are known
parameters, $Z$ is the normalising constant and the notation $\langle
ij\rangle$ indicates that sites i and j are nearest neighbours.

We first consider a base target distribution for which $n=50^{2}$, the spatial
correlation is moderate and more precisely $\lambda:=0.5$, and which has the
external field (the values for the $\alpha_{i}$’s) presented in Figure 2.

Figure 2: External field of the base target

We generated the $\alpha_{i}$ independently as follows:
$\alpha_{i}=-\mu+\epsilon_{i}$ if the column index is smaller than or equal to
$\ell=\lfloor\eta/2\rfloor$ and $\alpha_{i}=\mu+\epsilon_{i}$ otherwise, where
$\mu=1$, the $\epsilon_{i}$ are independent uniform random variables on the
interval $(-0.1,+0.1)$ and $\lfloor\cdot\rfloor$ is the floor function. In the
simulation study, while keeping the other parameters fixed, we gradually
increase $\eta$ from 50 to 500 to observe the impact of dealing with larger
systems. Next, while keeping the other parameters fixed (with $\eta=50$), we
gradually increase the value of $\mu$ from 1 to 3. This has for effect of
increasing the contrast so that it is clearer which values the $x_{i}$ should
take, thus making the target rougher and concentrated on fewer configurations.
One could vary $\lambda$ as well, but this would also make the target rougher
and concentrated on fewer configurations.

We observe the impact on Algorithm 1 with uniform and informed proposal
distributions, and their non-lifted MH counterparts. For such a simulation
study, it would be simply too long to obtain the results for Algorithm 2 with
$\rho_{\nu}^{*}$ (5). We tried to vary the value of $\ell$ to observe what
happens with the acceptance rates for the uniform lifted sampler when the
ratios $n_{-\nu}(\mathbf{x})/n_{\nu}(\mathbf{y})$ (see (6)) are far from 1.
The impact in this example is however not the one expected: the acceptance
rates increase instead of deteriorating. To see why, consider for instance
that $\ell=5$. The ratios $n_{-\nu}(\mathbf{x})/n_{\nu}(\mathbf{y})$ with
$\nu=-1$ are on average around 9 (45 columns with $+1$’s and 5 with $-1$’s).
With $\nu=-1$, it is tried to flip a bit from $+1$ to $-1$ and
$\pi(\mathbf{y})/\pi(\mathbf{x})$ is thus multiplied by a factor of around 9.
It is likely that this bit is on the yellow side (see Figure 2). For such a
move, $\pi(\mathbf{y})/\pi(\mathbf{x})$ is often around $\exp(-2(1+0.5\times
4))=0.002$. Therefore, it is more likely to accept this move compared to in
the reversible MH sampler (in the MH sampler $\pi(\mathbf{y})/\pi(\mathbf{x})$
is not multiplied by 9). When $\nu=+1$, the multiplicative factor is thus
around $1/9$, but it is relatively likely that the proposal will be to flip a
bit from $-1$ to $+1$ on the yellow side, because there are some bits with the
value $-1$ on this larger side, and this move is often automatically accepted
(because $\pi(\mathbf{y})/\pi(\mathbf{x})$ is often around $1/0.002)$. Note
that these conclusions are not in contradiction with our analysis of Section
3.1, because what we observed here is specific to the Ising model. We do not
present the results because the graph is uninteresting: the performance is
essentially constant for the informed samplers and that of the uniform ones is
so low that we do not see the ESS vary.

We present the other results in Figure 3. They are based on 1,000 independent
runs of 100,000 iterations for each algorithm and each value of $\mu$ and
$\eta$, with burn-ins of 10,000. For each run, an ESS per iteration is
computed for $f(\mathbf{x},-1)=f(\mathbf{x},+1)=\sum_{i}x_{i}$ and then the
results are averaged out. This function is proportional to the magnetisation.
Monitoring such a statistic is relevant as a quicker variation of its value
(leading to a higher ESS) indicates that the whole state-space is explored
quicker.

For the base target (represented by the starting points on the left of the
graphs in Figure 3), the mass is concentrated on a manifold of several
configurations, which allows for persistent movement for informed samplers.
The lifted one indeed takes advantage of this; it is approximately 7 times
more efficient than its non-lifted counterpart. The gap widens as $\eta$
increases; it is approximately 20 and 70 times more efficient when $\eta$ is
3.2 and 10 times larger (i.e. when $n$ is 10 and 100 times larger),
respectively. We observed that the ratio of ESSs increases linearly with
$\eta$. The non-informed samplers both perform poorly (the lines are on top of
each other).

As $\mu$ increases, the target becomes rougher and concentrated on fewer
configurations. When the roughness and concentration level are too severe the
performance of the lifted informed sampler stagnates, whereas that of the non-
lifted MH sampler continues to improve. There are two reasons for this.
Firstly, the acceptance rates deteriorate more rapidly for the lifted than the
non-lifted sampler (as a consequence of the difference in the acceptance
probability, see (11)). Secondly, when the mass is concentrated on few
configurations, it leaves not much room for persistent movement for the lifted
sampler. The latter thus loses its avantage. Again, the non-informed samplers
both perform poorly (the lines are on top of each other).

$\begin{array}[]{cc}\includegraphics[width=216.81pt]{Fig_Ising_n.pdf}&\includegraphics[width=216.81pt]{Fig_Ising_mu.pdf}\cr\textbf{(a)}&\textbf{(b)}\end{array}$

Figure 3: ESS per iteration of
$f(\mathbf{x},-1)=f(\mathbf{x},+1)=\sum_{i}x_{i}$ for Algorithm 1 with uniform
and informed proposal distributions and their non-lifted MH counterparts when:
(a) $\eta$ increases from 50 to 500 and the other parameters are kept fixed
($\mu=1$, $\lambda=0.5$ and $\ell=25$); (b) $\mu$ increases from 1 to 3 and
the other parameters are kept fixed ($\eta=50$, $\lambda=0.5$ and $\ell=25$)

###  4.2 Variable selection: US crime data

A study of crime rate was first presented in Erhlich (1973) and then an
expended and corrected version appeared in Vandaele (1978) in which corrected
data were provided; the topic was more precisely on the connection between
crime rate and 15 covariates (some were added by Vandaele (1978)) such as
percentage of males of age between 14 and 23 and mean years of schooling in a
given state. The data were indeed aggregated by state and were from 47 U.S.
states in 1960. They were analysed in several statistics papers (see, e.g.,
Raftery et al. (1997)) and are available in the R package MASS.

The data are modelled using the usual linear regression with normal errors.
Here we set the prior distribution of the regression coefficients and scaling
of the errors to be, conditionally on a model, the non-informative Jeffreys
prior. It is proved in Gagnon (2019) that a simple modification to the uniform
prior on the model random variable (represented here by $\mathbf{X}$) prevents
the Jeffreys–Lindley (Lindley, 1957; Jeffreys, 1967) paradox from arising.
With the resulting likelihood function and prior density on the parameters,
the latter can be integrated out. It is thus possible to evaluate the exact
marginal posterior probability of any of the $2^{15}=$ 32,768 models, up to a
normalising constants. This allows us to implement the MH sampler with the
locally informed proposal distribution of Zanella (2019) and its lifted
counterparts (Algorithm 1 and Algorithm 2 with $\rho_{\nu}^{*}$ (5)). In the
previous statistical studies, it was noticed that the mass is diffused over
several models, so that we expect the lifted chains to exhibit persistent
movement (as seen in Figure 1) and to perform well. To simplify the
presentation, we do not show the performance of the uniform samplers because,
as in the previous section, it is very poor.

The performances of the algorithms are summarised in Figure 4. The results are
based on 1,000 independent runs of 10,000 iterations for each algorithm, with
burn-ins of 1,000. Each run is started from a distribution which approximates
the target. On average, Algorithm 1 and Algorithm 2 with $\rho_{\nu}^{*}$ are
$2.7$ and $3.3$ times more efficient than their non-lifted counterpart,
respectively. The benefits of persistent movement thus compensate for a
decrease in acceptance rates; the rate indeed decreases from 0.92 for the MH
sampler to 0.71 for Algorithm 1 and Algorithm 2 with $\rho_{\nu}^{*}$ (5).

Figure 4: ESS per iteration for
$f(\mathbf{x},-1)=f(\mathbf{x},+1)=\sum_{i}x_{i}$ of 1,000 independent runs
for the MH sampler with the locally informed proposal distribution and its
lifted counterparts (Algorithm 1 and Algorithm 2 with $\rho_{\nu}^{*}$)

##  5 Lifted trans-dimensional sampler

In this section, we introduce a trans-dimensional algorithm which is a non-
reversible version of the popular reversible jump (RJ) algorithms introduced
by Green (1995). We consider that $\bm{\mathcal{X}}$ is a model space and
$\mathbf{X}$ a model indicator. The latter indicates, for instance, with 0’s
and 1’s which covariates are included in the model in variable selection
contexts as in Section 4.2. Such an algorithm is useful when it is not
possible to integrate out the parameters, contrarily to the linear regression
with normal errors and suitable priors. Examples of such situations include
analyses based on linear regression with super heavy-tailed errors ensuring
whole robustness (Gagnon et al., 2018) and generalised linear models and
generalised linear mixed models (Forster et al., 2012).

The parameters of a given model $\mathbf{x}$ are denoted by
$\bm{\theta}_{\mathbf{x}}\in\bm{\Theta}_{\mathbf{x}}$. Trans-dimensional
algorithms are MCMC methods that one uses to sample from a target distribution
$\pi$ defined on a union of sets
$\cup_{\mathbf{x}\in\bm{\mathcal{X}}}\\{\mathbf{x}\\}\times\bm{\Theta}_{\mathbf{x}}$,
which corresponds in Bayesian statistics to the joint posterior of the model
indicator $\mathbf{X}$ and the parameters of Model $\mathbf{X}$,
$\bm{\theta}_{\mathbf{X}}$. Such a posterior allows to jointly infer about
$(\mathbf{X},\bm{\theta}_{\mathbf{X}})$, or in other words, simultaneously
achieve model selection and parameter estimation. In this section, we assume
for simplicity that the parameters of all models are continuous random
variables.

Given the current state of the Markov chain
$(\mathbf{x},\bm{\theta}_{\mathbf{x}})$, a trans-dimensional sampler generates
the next state by first proposing a model candidate $\mathbf{y}\sim
q_{\mathbf{x}}(\mathbf{y})$ and then a proposal for its corresponding
parameter values. When $\mathbf{y}=\mathbf{x}$, we say that a parameter update
is proposed, whereas we say that a model switch is proposed when
$\mathbf{y}\neq\mathbf{x}$. Note that $\mathbf{x}\in\mathcal{N}(\mathbf{x})$,
contrarily to the previous sections. This is to allow parameter updates. When
a parameter update is proposed, we allow any fixed-dimensional methods to be
used; we only require that the Markov kernels leave the conditional
distributions $\pi(\,\cdot\mid\mathbf{x})$ invariant. When a model switch is
proposed, a vector of auxiliary variables
$\mathbf{u}_{\mathbf{x}\mapsto\mathbf{y}}$ is typically generated and this is
followed by a proposal mechanism leading to
$(\bm{\theta}^{\prime}_{\mathbf{y}},\mathbf{u}_{\mathbf{y}\mapsto\mathbf{x}})$,
where $\bm{\theta}^{\prime}_{\mathbf{y}}$ is the proposal for the parameter
values in Model $\mathbf{y}$. We require the whole proposal mechanism for
$\bm{\theta}^{\prime}_{\mathbf{y}}$ to be valid in a RJ framework, in the
sense that the model switch transitions are reversible in this framework. The
non-reversibility in the lifted trans-dimensional sampler lies in the
transitions for the $\mathbf{x}$ variable during model switches. More
precisely, $\mathbf{y}$ is generated from $q_{\mathbf{x},\nu}(\mathbf{y})$
instead, but the proposal mechanism for $\bm{\theta}^{\prime}_{\mathbf{y}}$
during model switches is the same. We consider that the acceptance probability
of these model switches in RJ is given by

$\displaystyle\alpha_{\text{RJ}}((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}})):=1\wedge\frac{q_{\mathbf{y}}(\mathbf{x})}{q_{\mathbf{x}}(\mathbf{y})}\,r((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}})),$
(14)

where the function $r$ depends on the method and may depend on several other
variables.

The algorithm is now presented in Algorithm 3. In it, we consider that the
current model $\mathbf{x}$ always belongs to the neighbourhood
$\mathcal{N}_{\nu}(\mathbf{x})$, regardless of the current direction $\nu$,
and that $q_{\mathbf{x},-1}(\mathbf{x})=q_{\mathbf{x},+1}(\mathbf{x})$, which
will typically be the case in practice.

Algorithm 3 A lifted trans-dimensional sampler for partially ordered model
spaces

  1. 1.

Generate $\mathbf{y}\sim q_{\mathbf{x},\nu}$, a PMF with support restricted to
$\mathcal{N}_{\nu}(\mathbf{x})$.

  2. 2.(a)

If $\mathbf{y}=\mathbf{x}$, attempt a parameter update using a MCMC kernel of
invariant distribution $\pi(\,\cdot\mid\mathbf{x})$ while keeping the current
value of the model indicator $\mathbf{x}$ and direction $\nu$ fixed.

  3. 2.(b)

If $\mathbf{y}\neq\mathbf{x}$, attempt a model switch from Model $\mathbf{x}$
to Model $\mathbf{y}$. Generate $\bm{\theta}^{\prime}_{\mathbf{y}}$ using a
method which is valid in RJ and $u\sim\mathcal{U}[0,1]$. If

$\displaystyle
u\leq\alpha_{\text{NRJ}}((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}})):=1\wedge\frac{q_{\mathbf{y},-\nu}(\mathbf{x})}{q_{\mathbf{x},\nu}(\mathbf{y})}\,r((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}})),$
(15)

set the next state of the chain to
$(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}},\nu)$. Otherwise, set it to
$(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)$.

  4. 3.

Go to Step 1.

######  Proposition 2.

The transition kernel of the Markov chain
$\\{(\mathbf{X},\bm{\theta}_{\mathbf{X}},\nu)(m):m\in\operatorname{\mathbb{N}}\\}$
simulated by Algorithm 3 admits $\pi\otimes\mathcal{U}\\{-1,1\\}$ as invariant
distribution.

###### Proof.

See Section 7. ∎

The main difficulty with the implementation of trans-dimensional samplers is
the construction of efficient proposal schemes for
$(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}})$ during model switches. Gagnon
(2019) discusses this in depth in the RJ framework. The author proposes a
scheme and proves that it is possible to arbitrarily approach an asymptotic
regime in which one is able to generate $\bm{\theta}^{\prime}_{\mathbf{y}}$
from $\pi(\,\cdot\mid\mathbf{y})$ (the correct conditional distribution) and
evaluate exactly the ratios of marginal probabilities
$\pi(\mathbf{y})/\pi(\mathbf{x})$ (and is therefore able to adequately
construct $q_{\mathbf{x},\nu}$). In particular, for this scheme,
$r((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}}))$
is a consistent estimator of $\pi(\mathbf{y})/\pi(\mathbf{x})$. We refer the
reader to that paper for the details.

We thus conclude that the marginal behaviour of
$\\{(\mathbf{X},\nu)(m):m\in\operatorname{\mathbb{N}}\\}$ is the same as that
of the stochastic process generated by Algorithm 1 in the asymptotic regime
and considering only iterations for which model switches are proposed. All
conclusions previously drawn thus hold, at least approximatively. In
particular, one may analyse the same data as in Section 4.2, but using the
super heavy-tailed regression of Gagnon et al. (2018) for robust inference and
outlier detection. The results would be essentially the same because, as
mentioned in Raftery et al. (1997), “standard diagnostic checking (see, e.g.,
Draper and Smith (1981)) did not reveal any gross violations of the
assumptions underlying normal linear regression” and the robust method is
designed for leading to similar results in the absence of outliers. We thus
omit further analysis of Algorithm 3 and we do not illustrate how it performs
for brevity. We nevertheless highlight that it is important for
$r((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}^{\prime}_{\mathbf{y}}))$
to be a low variance estimator of $\pi(\mathbf{y})/\pi(\mathbf{x})$ as
persistent movement may be interrupted otherwise, as shown in Gagnon and
Doucet (2019).

##  6 Discussion

In this paper, we presented and analysed generic algorithms allowing
straightforward sampling from any PMF $\pi$ with a support $\bm{\mathcal{X}}$
on which a partial order can be established. The algorithms rely on the
technique of lifting. We showed that these are expected to perform well when
the shape of target (the level of concentration of the mass) allows for
persistent movement. This is true even when the target concentrates on a
manifold of the ambient space in the case where the lifting technique is
combined with locally informed proposal distributions (provided that the shape
of the manifold allows for persistent movement).

The samplers are in particular useful for the simulation of binary random
variables and variable selection. Algorithm 1 can be directly employed for the
latter when the parameters of the models can be integrated out. A lifted
trans-dimensional sampler for partially ordered model spaces have been
introduced in Section 5 for, among others, variable selection when it is not
possible to integrate out the parameters. We believe it would be interesting
to continue this line of research by taking steps towards automatic generic
samplers using the technique of lifting for any discrete state-space.

## References

  * Andrieu and Livingstone (2019) Andrieu, C. and Livingstone, S. (2019) Peskun-Tierney ordering for Markov chain and process Monte Carlo: beyond the reversible scenario.  arXiv:1906.06197.
  * Barker (1965) Barker, A. A. (1965) Monte Carlo calculations of the radial distribution functions for a proton-electron plasma.  Austral. J. Phys., 18, 119–134.
  * Chen et al. (1999) Chen, F., Lovász, L. and Pak, I. (1999) Lifting Markov chains to speed up mixing.  In Proceedings of the thirty-first annual ACM symposium on Theory of computing, 275–281.
  * Diaconis et al. (2000) Diaconis, P., Holmes, S. and Neal, R. M. (2000) Analysis of a nonreversible Markov chain sampler.  Ann. Appl. Probab., 726–752.
  * Draper and Smith (1981) Draper, N. R. and Smith, H. (1981) Applied regression analysis (2nd ed.).  New York: Wiley.
  * Erhlich (1973) Erhlich, I. (1973) Participation in illegitimate activities: a theoretical and empirical analysis.  J. Polit. Econ., 81, 521–567.
  * Forster et al. (2012) Forster, J. J., Gill, R. C. and Overstall, A. M. (2012) Reversible jump methods for generalised linear models and generalised linear mixed models.  Stat. Comput., 22, 107–120.
  * Gagnon (2019) Gagnon, P. (2019) A step further towards automatic and efficient reversible jump algorithms.  ArXiv:1911.02089.
  * Gagnon et al. (2018) Gagnon, P., Desgagné, A. and Bédard, M. (2018) A new Bayesian approach to robustness against outliers in linear regression.  Bayesian Anal. Advance publication.
  * Gagnon and Doucet (2019) Gagnon, P. and Doucet, A. (2019) Non-reversible jump algorithms for Bayesian nested model selection.  ArXiv:1911.01340.
  * Geyer (1991) Geyer, C. J. (1991) Markov chain monte carlo maximum likelihood.  Interface Proceedings.
  * Green (1995) Green, P. J. (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika, 82, 711–732.
  * Gustafson (1998) Gustafson, P. (1998) A guided walk Metropolis algorithm.  Stat. Comput., 8, 357–364.
  * Hastings (1970) Hastings, W. K. (1970) Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57, 97–109.
  * Horowitz (1991) Horowitz, A. M. (1991) A generalized guided Monte carlo algorithm.  Phys. Lett. B, 268, 247–252.
  * Jeffreys (1967) Jeffreys, H. (1967) Theory of Probability.  Oxford Univ. Press, London.
  * Lindley (1957) Lindley, D. V. (1957) A statistical paradox.  Biometrika, 44, 187–192.
  * Livingstone and Zanella (2019) Livingstone, S. and Zanella, G. (2019) On the robustness of gradient-based MCMC algorithms.  arXiv:1908.11812.
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E. (1953) Equation of state calculations by fast computing machines.  J. Chem. Phys., 21, 1087.
  * Neal (2003) Neal, R. M. (2003) Slice sampling.  Ann. Statist., 705–741.
  * Neal (2020) — (2020) Non-reversibly updating a uniform [0, 1] value for Metropolis accept/reject decisions.  arXiv:2001.11950.
  * Peskun (1973) Peskun, P. (1973) Optimum Monte-Carlo sampling using Markov chains.  Biometrika, 60, 607–612.
  * Power and Goldman (2019) Power, S. and Goldman, J. V. (2019) Accelerated sampling on discrete spaces with non-reversible Markov processes.  arXiv:1912.04681.
  * Raftery et al. (1997) Raftery, A. E., Madigan, D. and Hoeting, J. A. (1997) Bayesian model averaging for linear regression models.  J. Amer. Statist. Assoc., 92, 179–191.
  * Roberts and Rosenthal (2004) Roberts, G. O. and Rosenthal, J. S. (2004) General state space Markov chains and MCMC algorithms.  Probab. Surv., 1, 20–71.
  * Syed et al. (2019) Syed, S., Bouchard-Côté, A., Deligiannidis, G. and Doucet, A. (2019) Non-reversible parallel tempering: an embarassingly parallel MCMC scheme.  arXiv:1905.02939.
  * Tierney (1998) Tierney, L. (1998) A note on Metropolis-Hastings kernels for general state spaces.  Ann. Appl. Probab., 8, 1–9.
  * Trotter (1992) Trotter, W. T. (1992) Combinatorics and partially ordered sets: Dimension theory, vol. 59.  Johns Hopkins University Press Baltimore.
  * Vandaele (1978) Vandaele, W. (1978) Participation in illegitimate activities; Ehrlich revisited.  In Deterrence and incapacitation, 270–335. Washington, D.C.: National Academy of Sciences Press.
  * Vanetti et al. (2017) Vanetti, P., Bouchard-Côté, A., Deligiannidis, G. and Doucet, A. (2017) Piecewise-deterministic Markov chain Monte Carlo.  arXiv:1707.05296.
  * Zanella (2019) Zanella, G. (2019) Informed proposals for local MCMC in discrete spaces.  To appear in J. Amer. Statist. Assoc.

##  7 Proofs

###### Proof of Proposition 1.

It suffices to prove that the probability to reach the state
$(\mathbf{y},\nu^{\prime})$ in one step is equal to the probability of this
state under the target:

$\displaystyle\sum_{\mathbf{x},\nu}\pi(\mathbf{x})\,(1/2)\,P((\mathbf{x},\nu),(\mathbf{y},\nu^{\prime}))=\pi(\mathbf{y})\,(1/2).$
(16)

where $P$ is the transition kernel.

The probability to reach the state $(\mathbf{y},\nu^{\prime})$ from some
$(\mathbf{x},\nu)$ is given by:

$\displaystyle P((\mathbf{x},\nu),(\mathbf{y},\nu^{\prime}))$
$\displaystyle=T_{\nu}(\mathbf{x},\bm{\mathcal{X}})\,Q_{\mathbf{x},\nu}(\mathbf{y})\,\mathds{1}(\nu=\nu^{\prime})$
$\displaystyle\qquad+\mathds{1}(\nu=-\nu^{\prime},\mathbf{x}=\mathbf{y})\left[(\rho_{\nu}(\mathbf{x})+T_{\nu}(\mathbf{x},\bm{\mathcal{X}}))-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})\right]$
$\displaystyle\qquad+\mathds{1}(\nu=\nu^{\prime},\mathbf{x}=\mathbf{y})\left[1-\rho_{\nu}(\mathbf{x})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})\right]$
$\displaystyle=q_{\mathbf{x},\nu}(\mathbf{y})\,\alpha_{\nu}(\mathbf{x},\mathbf{y})\,\mathds{1}(\nu=\nu)$
$\displaystyle\qquad+\mathds{1}(\nu=-\nu^{\prime},\mathbf{x}=\mathbf{y})\,\rho_{\nu}(\mathbf{x})$
$\displaystyle\qquad+\mathds{1}(\nu=\nu^{\prime},\mathbf{x}=\mathbf{y})\left[1-\rho_{\nu}(\mathbf{x})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})\right].$

We have that

$\displaystyle\pi(\mathbf{x})\,(1/2)\,P((\mathbf{x},\nu),(\mathbf{y},\nu^{\prime}))$
$\displaystyle=(1/2)\,\pi(\mathbf{y})\,q_{\mathbf{y},-\nu^{\prime}}(\mathbf{x})\,\alpha_{-\nu^{\prime}}(\mathbf{y},\mathbf{x})\,\mathds{1}(-\nu^{\prime}=-\nu)$
$\displaystyle\qquad+(1/2)\,\pi(\mathbf{y})\,\mathds{1}(-\nu^{\prime}=\nu,\mathbf{y}=\mathbf{x})\,\rho_{-\nu^{\prime}}(\mathbf{y})$
$\displaystyle\qquad+(1/2)\,\pi(\mathbf{y})\mathds{1}(-\nu^{\prime}=-\nu,\mathbf{y}=\mathbf{x})\left[1-\rho_{-\nu^{\prime}}(\mathbf{y})-T_{-\nu^{\prime}}(\mathbf{y},\bm{\mathcal{X}})\right]$
$\displaystyle=(1/2)\,\pi(\mathbf{y})\,T_{-\nu^{\prime}}(\mathbf{y},\bm{\mathcal{X}})\,Q_{\mathbf{y},-\nu^{\prime}}(\mathbf{x})\mathds{1}(-\nu^{\prime}=-\nu)$
$\displaystyle\qquad+(1/2)\,\pi(\mathbf{y})\,\mathds{1}(-\nu^{\prime}=\nu,\mathbf{y}=\mathbf{x})\left[(\rho_{-\nu^{\prime}}(\mathbf{y})+T_{-\nu^{\prime}}(\mathbf{y},\mathcal{X}))-T_{-\nu^{\prime}}(\mathbf{y},\mathcal{X})\right]$
$\displaystyle\qquad+(1/2)\,\pi(\mathbf{y})\,\mathds{1}(-\nu^{\prime}=-\nu,\mathbf{y}=\mathbf{x})\left[1-\rho_{-\nu^{\prime}}(\mathbf{y})-T_{-\nu^{\prime}}(\mathbf{y},\mathcal{X})\right],$

where we used the definition of $\alpha$ for the first term and that
$\rho_{\nu}(\mathbf{x})-\rho_{-\nu}(\mathbf{x})=T_{-\nu}(\mathbf{x},\bm{\mathcal{X}})-T_{\nu}(\mathbf{x},\bm{\mathcal{X}})$
for the third term. Notice the sum on the right-hand side (RHS) is equal to
the probability to reach some $(\mathbf{x},-\nu)$, starting from
$(\mathbf{y},-\nu^{\prime})$:
$(1/2)\,\pi(\mathbf{y})\,P((\mathbf{y},-\nu^{\prime}),(\mathbf{x},-\nu))$.

Therefore,

$\displaystyle\sum_{\mathbf{x},\nu}\pi(\mathbf{x})\,(1/2)\,P((\mathbf{x},\nu),(\mathbf{y},\nu^{\prime}))$
$\displaystyle=\sum_{\mathbf{x},\nu}(1/2)\,\pi(\mathbf{y})\,P((\mathbf{y},-\nu^{\prime}),(\mathbf{x},-\nu))$
$\displaystyle=(1/2)\,\pi(\mathbf{y}),$

which concludes the proof. ∎

We now present a lemma that will be useful in the next proofs.

######  Lemma 1.

Let $Q$ be the Markov kernel of the Markov chain simulated by Algorithm 2, for
any valid switching function $\rho_{\nu}$, $\nu\in\\{-1,1\\}$. Assume that
$\bm{\mathcal{X}}$ is finite. Then, for any
$f\in\mathcal{L}_{2}^{*}(\bar{\pi})$,

$\lim_{\lambda\to
1}\sum_{k>0}\lambda^{k}\langle\,f,Q^{k}f\,\rangle=\sum_{k>0}\langle\,f,Q^{k}f\,\rangle\,.$
(17)

###### Proof.

For each $f\in\mathcal{L}_{2}^{*}(\bar{\pi})$, define the sequence of
functions $S_{n}:\lambda\mapsto\sum_{0<k\leq
n}\lambda^{k}\langle\,f,Q^{k}f\,\rangle$ defined for $\lambda\in[0,1)$ and its
limit $S(\lambda)=\sum_{k>0}\lambda^{k}\langle\,f,Q^{k}f\,\rangle$ (the
dependance of $S_{n}$ and $S$ on $f$ and $Q$ is implicit). We now show that
the partial sum $S_{n}$ converges uniformly to $S$ on $[0,1)$, and since for
each $n\in\mathbb{N}$, the function
$\lambda\to\lambda^{n}\langle\,f,Q^{n}f\,\rangle$ admits a limit when
$\lambda\to 1$, we have that $S$ admits a limit when $\lambda\to 1$, given by

$\lim_{\lambda\to 1}S(\lambda)=S(1)=\sum_{k>0}\langle\,f,Q^{k}f\,\rangle\,,$

which is Eq. (17). First, note that

$\sup_{\lambda\in[0,1)}\left|S_{n}(\lambda)-S(\lambda)\right|=\sup_{\lambda\in[0,1)}\left|\sum_{k>n}\lambda^{k}\langle\,f,Q^{k}f\,\rangle\right|\leq\sup_{\lambda\in[0,1)}\sum_{k>n}\lambda^{k}\left|\langle\,f,Q^{k}f\,\rangle\right|=\sum_{k>n}\left|\langle\,f,Q^{k}f\,\rangle\right|\,.$
(18)

Thus, to prove that
$\sup_{\lambda\in[0,1)}\left|S_{n}(\lambda)-S(\lambda)\right|\to 0$, it is
sufficient to prove that the series
$\sum_{k>0}\left|\langle\,f,Q^{k}f\,\rangle\right|$ converges.

By bilinearity of the inner product and by linearity of the iterated operators
$Q,Q^{2},\ldots$, it can be checked that for any linear mapping $\phi$ on
$\mathcal{L}_{2}^{\ast}(\bar{\pi})$

$\sum_{k=1}^{\infty}\left|\left\langle
f,Q^{k}f\right\rangle\right|<\infty\Leftrightarrow\sum_{k=1}^{\infty}\left|\left\langle\phi(f),Q^{k}\phi(f)\right\rangle\right|<\infty\,.$
(19)

Since $\bm{\mathcal{X}}$ is finite, if $f\in\mathcal{L}_{2}^{\ast}(\bar{\pi})$
then $\sup|f|<\infty$. As a consequence, we may use
$\phi(f):=(f-\bar{\pi}f)/\sup|f|$ and $\bar{\pi}f:=\int f\mathrm{d}\bar{\pi}$.
In the following we denote by $\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})$ the
subset of $\mathcal{L}_{2}^{\ast}(\bar{\pi})$ such that

$\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi}):=\left\\{f\in\mathcal{L}_{2}^{\ast}(\bar{\pi})\,:\,\bar{\pi}f=0\,,\;\sup|f|\leq
1\right\\}\,.$

By Eq. (19), we only need to check that the series
$\sum_{k>0}\left|\langle\,f,Q^{k}f\,\rangle\right|$ converges for each
$f\in\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})$. Since $\bm{\mathcal{X}}$ is
finite, $Q$ is uniformly ergodic and there exists constants $\varrho\in(0,1)$
and $C\in(0,\infty)$ such that for any $t\in\mathbb{N}$,

$\sup_{(x,\nu)\in\bm{\mathcal{X}}\times\\{-1,+1\\}}\|\delta_{x,\nu}Q^{t}-\bar{\pi}\|_{\mathrm{tv}}\leq
C\varrho^{t}\,,\ $ (20)

where for any signed measure $\mu$, $\|\mu\|_{\mathrm{tv}}$ denotes its total
variation. On the one hand, note that for each
$f\in\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})$

$\langle\,f,Q^{k}f\,\rangle=\mathbb{E}f(X,\nu)Q^{k}f(X,\nu)\leq\mathbb{E}|f(X,\nu)||Q^{k}f(X,\nu)|\,,$
(21)

and on the other hand, we have that for any
$(x,\nu)\in\bm{\mathcal{X}}\times\\{-1,1\\}$,

$|Q^{k}f(x,\nu)|=|Q^{k}f(x,\nu)-\bar{\pi}f|\leq\sup_{f\in\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})}|Q^{k}f(x,\nu)-\bar{\pi}f|\,.$
(22)

But $\|\mu\|_{\mathrm{tv}}=(1/2)\sup_{g:\bm{\mathcal{X}}\to[-1,1]}|\mu g|$,
see for instance (Roberts and Rosenthal, 2004, Proposition 3). Since
$f\in\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})$, $|f|\leq 1$ and we have by
inclusion that for all $(x,\nu)\in\bm{\mathcal{X}}\times\\{-1,1\\}$

$|Q^{k}f(x,\nu)|\leq\sup_{f\in\mathcal{L}_{2}^{\ast,0,1}(\bar{\pi})}|Q^{k}f(x,\nu)-\bar{\pi}f|\leq\sup_{g:\bm{\mathcal{X}}\times\\{-1,1\\}\to[-1,1]}|Q^{k}g(x,\nu)-\bar{\pi}g|\leq
2\|\delta_{x,\nu}Q^{k}-\bar{\pi}\|_{\mathrm{tv}}\,.$ (23)

Taking the supremum over all $(x,\nu)\in\bm{\mathcal{X}}\times\\{-1,1\\}$ in
Eq. (23), and combining with Eq. (20) yields

$\sup_{(x,\nu)\in\bm{\mathcal{X}}\times\\{-1,1\\}}|Q^{k}f(x,\nu)|\leq
2C\varrho^{k}\,.$

Plugging this into Eq. (21), we have

$\left|\langle\,f,Q^{k}f\,\rangle\right|\leq
2C\mathbb{E}|f(X,\nu)|\varrho^{k}\,.$ (24)

which is clearly summable. As a consequence, $S_{n}$ converges uniformly to
$S$ on $[0,1)$ which concludes the proof. ∎

###### Proof of Corollary 1.

The results of Theorem 3.15 in Andrieu and Livingstone (2019) holds in our
framework, implying that

$\mathrm{var}_{\lambda}(f,P_{\rho^{*}})\leq\mathrm{var}_{\lambda}(f,P_{\rho})\leq\mathrm{var}_{\lambda}(f,P_{\rho}^{\text{w}}),$

where
$\mathrm{var}_{\lambda}(f,P_{\rho}):=\mathbb{V}\mathrm{ar}f(\mathbf{X},\nu)+2\sum_{k>0}\lambda^{k}\left\langle
f,P_{\rho}^{k}f\right\rangle$ with $\lambda\in[0,1)$. Lemma 1 allows to
conclude. ∎

###### Proof of Corollary 2.

The proof is an application of Theorem 3.17 in Andrieu and Livingstone (2019)
which will allow to establish that

$\mathrm{var}_{\lambda}(f,P_{\rho})\leq\mathrm{var}_{\lambda}(f,P_{\text{MH}}).$

We will thus be able to conclude using Lemma 1.

In order to apply Theorem 3.17, we must verify that

$q_{\mathbf{x}}(\mathbf{y})\,\alpha(\mathbf{x},\mathbf{y})=(1/2)\,q_{\mathbf{x},+1}(\mathbf{y})\,\alpha_{+1}(\mathbf{x},\mathbf{y})+(1/2)\,q_{\mathbf{x},-1}(\mathbf{y})\,\alpha_{-1}(\mathbf{x},\mathbf{y}),$

for all $\mathbf{x}$ and $\mathbf{y}$. This is straightforward to verify under
the assumptions of Corollary 2:

$\displaystyle(1/2)\,q_{\mathbf{x},+1}(\mathbf{y})\,\alpha_{+1}(\mathbf{x},\mathbf{y})+(1/2)\,q_{\mathbf{x},-1}(\mathbf{y})\,\alpha_{-1}(\mathbf{x},\mathbf{y})$
$\displaystyle\qquad=\frac{1}{2}\frac{1}{(|\mathcal{N}(\mathbf{x})|/2)}\left(1\wedge\frac{\pi(\mathbf{y})}{\pi(\mathbf{x})}\right)\mathds{1}_{\mathbf{y}\in\mathcal{N}_{+1}(\mathbf{x})}+\frac{1}{2}\frac{1}{(|\mathcal{N}(\mathbf{x})|/2)}\left(1\wedge\frac{\pi(\mathbf{y})}{\pi(\mathbf{x})}\right)\mathds{1}_{\mathbf{y}\in\mathcal{N}_{-1}(\mathbf{x})}$
$\displaystyle\qquad=\frac{1}{|\mathcal{N}(\mathbf{x})|}\left(1\wedge\frac{\pi(\mathbf{y})}{\pi(\mathbf{x})}\right)\left(\mathds{1}_{\mathbf{y}\in\mathcal{N}_{+1}(\mathbf{x})}+\mathds{1}_{\mathbf{y}\in\mathcal{N}_{-1}(\mathbf{x})}\right)=q_{\mathbf{x}}(\mathbf{y})\,\alpha(\mathbf{x},\mathbf{y}).$

∎

###### Proof of Theorem 1.

We first prove that
$\mathrm{var}(f,P_{\rho,n})\leq\mathrm{var}(f,P_{\text{rev.},n})$. This is
done as in the proof of Corollary 2.

We now analyse $\mathrm{var}(f,P_{\text{rev.},n})$:

$\mathrm{var}(f,P_{\text{rev.},n})=\mathbb{V}\mathrm{ar}f(\mathbf{X})+2\sum_{k>0}\left\langle
f,P_{\text{rev.},n}^{k}f\right\rangle=\mathbb{V}\mathrm{ar}f(\mathbf{X})+2\sum_{k>0}\text{Cov}_{P_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))],$

where we omitted the dependence on $\nu$ because, as we mentioned in Section
3.1, it can be treated as a constant as a consequence of the restrictions on
$f$. In the expression above, it is considered that the chain starts at
stationarity and evolves using $P_{\text{rev.},n}$. We consider without loss
of generality that $\mathbb{E}[f(\mathbf{X}(k))]=0$ (for any $k$).

We first write

$\displaystyle\mathbb{V}\mathrm{ar}f(\mathbf{X})$
$\displaystyle=\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}]$
(25)
$\displaystyle=\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})^{2}]-\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]$
(26)
$\displaystyle\qquad+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}],$
(27)

where $\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}$ denotes an
expectation with respect to $\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$.

Also,

$\displaystyle\sum_{k>0}\text{Cov}_{P_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))]=\sum_{k=1}^{\varrho(n)}\text{Cov}_{P_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))]+\sum_{k=\varrho(n)+1}^{\infty}\text{Cov}_{P_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))],$

where $\varrho(n)$ is chosen according to the statement of Theorem 1;
therefore the second term on the RHS can be made as small as we want. For the
first term, we have

$\displaystyle\sum_{k=1}^{\varrho(n)}\text{Cov}_{P_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))]$
$\displaystyle=\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))]$
$\displaystyle=\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cap_{m=0}^{k}A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1})}]+\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cup_{m=0}^{k}A_{m}^{\mathsf{c}}(\bm{\mathcal{X}}_{\varphi(n)-1})}],$

where
$A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1}):=\\{\mathbf{X}(m)\in\bm{\mathcal{X}}_{\varphi(n)-1}\\}$.
By assumption, the second term on the RHS can be made as small as we want.

We have

$\displaystyle\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cap_{m=0}^{k}A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1})}]$
$\displaystyle=\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{P}_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))]$
$\displaystyle+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\sum_{k=1}^{\varrho(n)}\mathbb{E}_{P_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))\mathds{1}_{\cap_{m=0}^{k}A_{m}(\bm{\mathcal{X}}_{\varphi(n)-1})}],$

where $\tilde{P}_{\text{rev.},n}$ is the Markov kernel whose stationary
distribution is $\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$. This equality
holds because the paths involving transition probabilities
$P_{\text{rev.},n}(\mathbf{x},\mathbf{y})$ with
$\mathbf{x},\mathbf{y}\in\bm{\mathcal{X}}_{\varphi(n)-1}$ have the same
probabilities as those of the chain with stationary distribution
$\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}$. We just have to renormalise the
probabilities of the starting point by dividing by
$\pi(\bm{\mathcal{X}}_{\varphi(n)})$ to complete the argument. Note that, by
assumption, the second term on the RHS can be made as small as we want.

Now,

$\displaystyle\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{P}_{\text{rev.},n}}[f(\mathbf{X}(0))f(\mathbf{X}(k))]=\sum_{k=1}^{\varrho(n)}\text{Cov}_{\tilde{P}_{\text{rev.},n}}[f(\mathbf{X}(0)),f(\mathbf{X}(k))]+\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}.$

By assumption, the sum from $\varrho(n)+1$ to $\infty$ of the covariances is
small as well. If we combine this with (25), we have

$\displaystyle\mathrm{var}(f,P_{\text{rev.},n})$
$\displaystyle=\mathrm{var}(f,\tilde{P}_{\text{rev.},n})+\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]$
$\displaystyle\qquad+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}]+2\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}$
$\displaystyle\leq\mathrm{var}(f,\tilde{P}_{\text{MH},n})/\omega_{n}+((1-\omega_{n})/\omega_{n})\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}f(\mathbf{X})+\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}$
$\displaystyle\qquad+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}]+2\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}$
$\displaystyle=\mathrm{var}(f,\tilde{P}_{\text{MH},n})+\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}]$
$\displaystyle\qquad+2\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+(1/\omega_{n}-1)\mathrm{var}(f,\tilde{P}_{\text{MH},n})+((1-\omega_{n})/\omega_{n})\mathbb{V}\mathrm{ar}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}f(\mathbf{X}),$

using Theorem 2 of Zanella (2019) for the inequality and omitting the terms
that can be made as small as we want. This theorem can be used as a result of
the bound
$\tilde{P}_{\text{rev.},n}(\mathbf{x},\mathbf{y})\geq\omega_{n}\tilde{P}_{\text{MH},n}(\mathbf{x},\mathbf{y})$
with $\omega_{n}$ defined in (9) (see Section 3.1). By assumption, the last
two terms can be made as small as we want.

Now we used the previous arguments in the reverse order to show that

$\displaystyle\mathrm{var}(f,\tilde{P}_{\text{MH},n})+\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}+(1-1/\pi(\bm{\mathcal{X}}_{\varphi(n)}))\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}}]+\mathbb{E}[f(\mathbf{X})^{2}\mathds{1}_{\mathbf{X}\in\bm{\mathcal{X}}_{\varphi(n)}^{\mathsf{c}}}]$
$\displaystyle\qquad+2\sum_{k=1}^{\varrho(n)}\mathbb{E}_{\tilde{\pi}_{\bm{\mathcal{X}}_{\varphi(n)}}}[f(\mathbf{X})]^{2}=\mathrm{var}(f,P_{\text{MH},n})+\text{small
error term},$

which yields the result. ∎

###### Proof of Proposition 2.

It suffices to prove that the probability to reach the state
$\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}\in A,\nu^{\prime}$ in one step
is equal to the probability of this state under the target:

$\displaystyle\sum_{\mathbf{x},\nu}\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime}))\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
$\displaystyle=\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\,d\bm{\theta}_{\mathbf{y}}^{\prime},$
(28)

where $P$ is the transition kernel. Note that we abuse notation here by
denoting the measure $d\bm{\theta}_{\mathbf{y}}^{\prime}$ on the left-hand
side (LHS) given that we may in fact use vectors of auxiliary variables to
generate the proposal when switching models, which often do not have the same
dimension as $\bm{\theta}_{\mathbf{y}}^{\prime}$.

We consider two distinct events: a model switch is proposed, that we denote
$S$, and a parameter update is proposed (therefore denoted $S^{\mathsf{c}}$).
We know that the probabilities of these events are
$1-q_{\mathbf{x},\nu}(\mathbf{x})$ and $q_{\mathbf{x},\nu}(\mathbf{x})$,
respectively. We rewrite the LHS of (28) as

$\displaystyle\sum_{\mathbf{x},\nu}\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime}))\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
(29)
$\displaystyle\quad=\sum_{\mathbf{x},\nu}\,(1-q_{\mathbf{x},\nu}(\mathbf{x}))\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S)\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$ (30)
$\displaystyle\qquad+\sum_{\mathbf{x},\nu}\,q_{\mathbf{x},\nu}(\mathbf{x})\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S^{\mathsf{c}})\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}.$
(31)

We analyse the two terms separately. We know that

$P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S^{\mathsf{c}})=\delta_{(\mathbf{x},\nu)}(\mathbf{y},\nu^{\prime})\,P_{S^{\mathsf{c}}}(\bm{\theta}_{\mathbf{x}},\bm{\theta}_{\mathbf{y}}^{\prime}),$

where $P_{S^{\mathsf{c}}}$ is the transition kernel associated with the method
used to update the parameters. Therefore, the second term on the RHS of (29)
is equal to

$\displaystyle\sum_{k,\nu}\,q_{\mathbf{x},\nu}(\mathbf{x})\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S^{\mathsf{c}})\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
$\displaystyle\quad=q_{\mathbf{y},\nu^{\prime}}(\mathbf{y})\,\pi(\mathbf{y})\times(1/2)\int\pi(\bm{\theta}_{\mathbf{y}}\mid\mathbf{y})\left(\int_{A}P_{S^{\mathsf{c}}}(\bm{\theta}_{\mathbf{y}},\bm{\theta}_{\mathbf{y}}^{\prime})\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{y}}.$

We also know that $P_{S^{\mathsf{c}}}$ leaves the conditional distribution
$\pi(\,\cdot\mid\mathbf{y})$ invariant, implying that

$\displaystyle
q_{\mathbf{y},\nu^{\prime}}(\mathbf{y})\,\pi(\mathbf{y})\times(1/2)\int\pi(\bm{\theta}_{\mathbf{y}}\mid\mathbf{y})\left(\int_{A}P_{S^{\mathsf{c}}}(\bm{\theta}_{\mathbf{y}},\bm{\theta}_{\mathbf{y}}^{\prime})\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{y}}$
(32)
$\displaystyle\quad=q_{\mathbf{y},\nu^{\prime}}(\mathbf{y})\,\pi(\mathbf{y})\times(1/2)\int_{A}\pi(\bm{\theta}_{\mathbf{y}}^{\prime}\mid\mathbf{y})\,d\bm{\theta}_{\mathbf{y}}^{\prime}=q_{\mathbf{y},\nu^{\prime}}(\mathbf{y})\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\,d\bm{\theta}_{\mathbf{y}}^{\prime}.$
(33)

For the model switching case (the first term on the RHS of (29)), we use the
fact that there is a connection between
$P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S)$ and the kernel associated to a specific RJ. Consider that
$q_{\mathbf{x}}(\mathbf{y})=(1/2)\,q_{\mathbf{x},-1}(\mathbf{y})+(1/2)\,q_{\mathbf{x},+1}(\mathbf{y})$
for all $\mathbf{x},\mathbf{y}$ and that all other proposal distributions in
RJ are the same as in Algorithm 3 during model switches. In this case,
$\alpha_{\text{RJ}}=\alpha_{\text{NRJ}}$ in the case where at the current
iteration it is chosen to use $q_{\mathbf{x},\nu}$ (which happens with
probability $1/2$) and in the reverse move it is chosen to use
$q_{\mathbf{y},-\nu}$ (which also happens with probability $1/2$).

Consider the case where Model $\mathbf{y}$ is reached from Model
$\mathbf{x}\neq\mathbf{y}$ coming from direction $\nu^{\prime}=\nu$. Given the
reversibility of RJ, the probability to go from Model $\mathbf{x}$ with
parameters in $B$ to Model $\mathbf{y}\neq\mathbf{x}$ with parameters in $A$
is

$\displaystyle\int_{B}\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\left(\int_{A}P_{\text{RJ}}((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}))\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}=\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\left(\int_{B}P_{\text{RJ}}((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}),(\mathbf{x},\bm{\theta}_{\mathbf{x}}))\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime},$
(34)

where $P_{\text{RJ}}$ is the transition kernel of the RJ. Note that

$P_{\text{RJ}}((\mathbf{x},\bm{\theta}_{\mathbf{x}}),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}))=(1/2)\,(1-q_{\mathbf{x},\nu}(\mathbf{x}))\,P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu)\mid
S),$

given that the difference between both kernels is that in RJ, it is first
decided to use $q_{\mathbf{x},\nu}$, there is thus an additional probability
of $1/2$. Analogously,
$P_{\text{RJ}}((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}),(\mathbf{x},\bm{\theta}_{\mathbf{x}}))=(1/2)\,(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\,P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)$. Using this and taking $B$ equals the whole parameter (and auxiliary)
space in (34), we have

$\displaystyle(1-q_{\mathbf{x},\nu}(\mathbf{x}))\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu)\mid
S)\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
$\displaystyle\qquad=(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\left(\int
P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime}.$

The only other case to consider for model switches is where Model $\mathbf{y}$
is reached from Model $\mathbf{y}$ (because the proposal is rejected) and the
direction is $\nu^{\prime}=-\nu$. The probability of this transition is

$(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\left(1-\sum_{\mathbf{x}\neq\mathbf{y}}\int
P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime}.$

So, the total probability of reaching
$\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime}\in A,\nu^{\prime}$ through a
model switch is (recalling (29)):

$\displaystyle\sum_{\mathbf{x},\nu}\,(1-q_{\mathbf{x},\nu}(\mathbf{x}))\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu^{\prime})\mid
S)\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
$\displaystyle\quad=\sum_{\mathbf{x}\neq\mathbf{y}}(1-q_{\mathbf{x},\nu}(\mathbf{x}))\int\pi(\mathbf{x},\bm{\theta}_{\mathbf{x}})\times(1/2)\left(\int_{A}P((\mathbf{x},\bm{\theta}_{\mathbf{x}},\nu),(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},\nu)\mid
S)\,d\bm{\theta}_{\mathbf{y}}^{\prime}\right)\,d\bm{\theta}_{\mathbf{x}}$
$\displaystyle\qquad+(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\left(1-\sum_{\mathbf{x}\neq\mathbf{y}}\int
P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime}$
$\displaystyle\quad=\sum_{\mathbf{x}\neq\mathbf{y}}(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\left(\int
P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime}$
$\displaystyle\qquad+(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2)\left(1-\sum_{\mathbf{x}\neq\mathbf{y}}\int
P((\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime},-\nu),(\mathbf{x},\bm{\theta}_{\mathbf{x}},-\nu)\mid
S)\,d\bm{\theta}_{\mathbf{x}}\right)\,d\bm{\theta}_{\mathbf{y}}^{\prime}$
$\displaystyle\quad=(1-q_{\mathbf{y},-\nu}(\mathbf{y}))\int_{A}\pi(\mathbf{y},\bm{\theta}_{\mathbf{y}}^{\prime})\times(1/2).$

Combining this with (32) allows to conclude the proof. ∎