∎

11institutetext: Changye WU 22institutetext: CEREMADE, Université Paris
Dauphine, PSL Research University, France
22email<EMAIL_ADDRESS>33institutetext: Christian P. ROBERT
44institutetext: CEREMADE, Université Paris Dauphine, PSL Research University,
France,
Department of Statistics, University of Warwick, UK,
Università Ca’ Foscari Venezia, Italy
44email<EMAIL_ADDRESS>
# Coordinate Sampler: A Non-Reversible Gibbs-like MCMC Sampler

Changye WU  Christian P. ROBERT

(Received: date / Accepted: date)

###### Abstract

We derive a novel non-reversible, continuous-time Markov chain Monte Carlo
(MCMC) sampler, called Coordinate Sampler, based on a piecewise deterministic
Markov process (PDMP), which is a variant of the Zigzag sampler of Bierkens et
al., (2016). In addition to providing a theoretical validation for this new
simulation algorithm, we show that the Markov chain it induces exhibits
geometrical ergodicity convergence, for distributions whose tails decay at
least as fast as an exponential distribution and at most as fast as a Gaussian
distribution. Several numerical examples highlight that our coordinate sampler
is more efficient than the Zigzag sampler, in terms of effective sample size.

###### Keywords:

Markov chain Monte Carlo Piecewise deterministic Markov processes Zigzag
sampling Gibbs sampling

##  1 Introduction

A powerful and generic sampling technique, the Markov chain Monte Carlo (MCMC)
method, (see, e.g., Robert and Casella,, 2004) has been widely exploited in
computational statistics to become a standard tool in Bayesian inference,
where posterior distributions are often analytically intractable and at best
known up to a normalizing constant. However, almost all existing MCMC
algorithms, such as the Metropolis-Hastings algorithm (MH), the Hamiltonian
Monte Carlo (HMC) (Neal et al.,, 2011) and Metropolis adjusted Langevin
algorithm (MALA), satisfy detailed balance conditions, dating back to
Metropolis et al., (1953) and Hastings, (1970). Recently, a different
technology of MCMC sampling – piecewise deterministic Markov process (PDMP) –
was introduced in computational statistics, towards removing reversibility
constraints. The basic theory of PDMP was developed in Davis, (1984) and
Davis, (1993), while an application to computational statistics was
implemented by, e.g., Peters et al., (2012), Bierkens et al., (2016), and
Bouchard-Côté et al., (2018).

Since piecewise deterministic Markov processes for sampling from distributions
was introduced by Peters et al., (2012), PDMP-based, continuous-time, non-
reversible, MCMC algorithms have become relevant tools, from applied
probability (Bierkens et al.,, 2017; Fontbona et al.,, 2016) to physics
(Peters et al.,, 2012; Harland et al.,, 2017; Michel et al.,, 2014), to
statistics (Bierkens et al.,, 2016; Fearnhead et al.,, 2018; Bierkens et al.,,
2018; Bouchard-Côté et al.,, 2018; Michel and Sénécal,, 2017; Vanetti et al.,,
2017; Pakman et al.,, 2016). However, almost all existing PDMP-based MCMC
samplers are based on two original versions: the Bouncy Particle Sampler (BPS)
of Bouchard-Côté et al., (2018) and and the Zigzag Sampler of Bierkens et al.,
(2016). Bouchard-Côté et al., (2018) exhibit that BPS can provide state-of-
the-art performance compared with the reference HMC for high dimensional
distributions, while Bierkens et al., (2016) show that the PDMP-based sampler
is easier to scale in big data settings, without introducing bias. Bierkens et
al., (2018) considers the application of PDMP for distributions on restricted
domains. Fearnhead et al., (2018) unify BPS and Zigzag samplers within the
framework of PDMPs: they propose a choice of the process velocity, at event
times, over the unit sphere, based on the angle between this velocity and the
gradient of the potential function. (This perspective relates to the
transition dynamics used here.) To overcome the main difficulty met by PDMP-
based samplers, namely the simulation of time-inhomogeneous Poisson processes,
Sherlock and Thiery, (2017) and Vanetti et al., (2017) resort to a
discretization of such continuous-time samplers. Furthermore, a pre-
conditioning of the velocity set is shown to accelerate the algorithms, see
Pakman et al., (2016).

In this article, we propose the Coordinate Sampler (CS), a novel PDMP-based
MCMC sampler that is a variant of the Zigzag sampler (ZS) of Bierkens et al.,
(2016). However, it differs from ZS in three significant aspects. First, the
velocity set considered in the coordinate sampler consists of an orthonormal
basis of the Euclidean space $\mathbb{R}^{d}$, while the one in the Zigzag
sampler is restricted to $\\{-1,1\\}^{d}$, if $d$ denotes the dimension of the
target distribution. Second, the event rate function in the Zigzag sampler is
much larger than the one for the coordinate sampler, especially for high
dimensional targets. This means that events occur more frequently in the
Zigzag sampler and hence this lowers its efficiency compared with our
approach. Thirdly, the coordinate sampler targets only one component at a time
when exploring the target space, and it keeps the other components unchanged,
while the Zigzag sampler modifies all components at the same time.
The outline of this article is as follows. Section 2 introduces the necessary
background of PDMP-based MCMC samplers, the techniques used in its
implementation, and two specified samplers, BPS and ZS. Section 3 describes
the methodology behind the coordinate sampler, provides some theoretical
validation along with a proof of geometrical ergodicity, obtained under quite
mild conditions, and compares this proposal with the Zigzag sampler in an
informal analysis. Section 4 further compares the efficiency of both
approaches on banana-shaped distributions, multivariate Gaussian distributions
and a Bayesian logistic model, when effective sample size is measuring
efficiency. Section 5 concludes by pointing out further research directions
about this special MCMC sampler.

##  2 Piecewise deterministic Markov process

In this section, we briefly introduce piecewise deterministic Markov processes
(PDMP) and describe how to apply this methodology into statistical computing
problems. We describe two specified PDMP-based MCMC samplers: the bouncy
particle sampler (BPS) and the Zigzag sampler (ZS).

###  2.1 PDMP-based Sampler

Let $\pi$ be the continuous target distribution over $\mathbb{R}^{d}$ and for
convenience sake, denote $\pi(\mathbb{x})$ for the probability density
function of $\pi$, when $\mathbb{x}\in\mathbb{R}^{d}$. We define
$U(\mathbb{x})$ as the potential function of $\pi(\mathbb{x})$, that is,
$\pi(\mathbb{x})\propto\exp\\{-U(\mathbb{x})\\}$, with $U$ positive. In the
PDMP framework, an auxiliary variable, $\mathbb{V}\in\mathcal{V}$, is
introduced and a PDMP-based sampler explores the augmented state space
$\mathbb{R}^{d}\times\mathcal{V}$, targeting a variable
$\mathbb{Z}=(\mathbb{X},\mathbb{V})$ with distribution
$\rho(d\mathbb{x},d\mathbb{v})$ over $\mathbb{R}^{d}\times\mathcal{V}$ as its
invariant distribution. By construction, the distribution $\rho$ enjoys $\pi$
as its marginal distribution in $\mathbb{x}$. In practice, the existing PDMP-
based samplers choose $\mathcal{V}$ to be the Euclidean space
$\mathbb{R}^{d}$, the sphere $\mathbb{S}^{d-1}$, or the discrete set
$\mathcal{V}\mathcal{V}==\\{\mathbb{v}=(v_{1},\cdots,v_{d})|v_{i}\in\\{-1,1\\},i=1,\cdots,d\\}$.
Following Fearnhead et al., (2018), a piecewise deterministic Markov process
$\mathbb{Z}_{t}=(\mathbb{X}_{t},\mathbb{V}_{t})$ consists of three distinct
components: a deterministic dynamic between events, an event occurrence rate,
and a transition dynamic at event times. Specifically,

  1. 1.

Deterministic dynamic: between two events, the Markov process evolves
deterministically, according to some ordinary differential equation:
$\frac{d\mathbb{Z}_{t}}{dt}=\Psi(\mathbb{Z}_{t})$.

  2. 2.

Event occurrence rate: an event occurs at time $t$ with rate
$\lambda(\mathbb{Z}_{t})$.

  3. 3.

Transition dynamic: At an event time, $\tau$, the state prior to $\tau$ is
denoted by $\mathbb{Z}_{\tau-}$, with the new state being generated by
$\mathbb{Z}_{\tau}\sim Q(\cdot|\mathbb{Z}_{\tau-})$.

Here, an “event” refers to an occurrence of a time-inhomogeneous Poisson
process with rate $\lambda(\cdot)$ (Kingman,, 1992). Following (Davis,, 1993,
Theorem 26.14), this Markov process had an extended generator equal to

$\mathcal{L}f(\mathbb{z})=\nabla
f(\mathbb{z})\cdot\Psi(\mathbb{z})+\lambda(\mathbb{z})\int_{\mathbb{z}^{\prime}}\left[f(\mathbb{z}^{\prime})-f(\mathbb{z})\right]Q(d\mathbb{z}^{\prime}|\mathbb{z})$
(1)

In order to guarantee invariance with respect to $\rho(d\mathbb{z})$, the
extended generator need satisfy
$\displaystyle{\int\mathcal{L}f(\mathbb{z})\rho(d\mathbb{z})}=0$ for all $f$
in an appropriate function class on $\mathbb{R}^{d}\times\mathcal{V}$ (Davis,,
1993, Theorem 34.7).

###  2.2 Implementation of a PDMP-based Sampler

In practice, choosing an appropriate deterministic dynamic, an event rate and
a transition dynamic, produces a Markov chain with invariant distribution
$\rho(d\mathbb{z})$. As for regular MCMC, generating such a Markov chain for a
duration $T$, leads to an estimator,
$\displaystyle{\frac{1}{T}\int_{t=0}^{T}h(\mathbb{X}_{t})dt}$, converging to
the integral of interest, $\displaystyle{I=\int
h(\mathbb{x})\pi(d\mathbb{x})}$, by the Law of Large Numbers for Markov
processes (Glynn and Haas,, 2006), under appropriate assumptions. More
specifically,

$\frac{1}{T}\int_{t=0}^{T}g(\mathbb{Z}_{t})dt\longrightarrow\int
g(\mathbb{z})\rho(d\mathbb{z}),\quad\text{as $T\rightarrow\infty$}$

and defining $g(\mathbb{z})=g(\mathbb{x},\mathbb{v}):=h(\mathbb{x})$ induces,
as $T\rightarrow\infty$,

$\frac{1}{T}\int_{t=0}^{T}h(\mathbb{X}_{t})dt=\frac{1}{T}\int_{t=0}^{T}g(\mathbb{Z}_{t})dt\rightarrow\int
g(\mathbb{z})\rho(d\mathbb{z})=\int h(\mathbb{x})\pi(d\mathbb{x}),$

where $p(d\mathbb{v}|\mathbb{x})$ is the conditional distribution of the
variable $\mathbb{V}$, given $\mathbb{X}=\mathbb{x}$. Algorithm 1 contains a
pseudo-code reproducing the simulation of a PDMP in practice:

Algorithm 1 General PDMP-based sampler

1:Input: start at position $\mathbb{x}_{0}$, velocity $\mathbb{v}_{0}$ and
simulation time threshold $T^{total}$.

2:Generate a set of event times of the PDMP
$\\{\tau_{0},\tau_{1},\cdots,\tau_{M}\\}$ and their associated states
$\\{\mathbb{Z}_{\tau_{0}},\mathbb{Z}_{\tau_{1}},\cdots,\mathbb{z}_{\tau_{M}}\\}$,
where $\tau_{0}=0$, $\tau_{M-1}<T^{total}$, $\tau_{M}\geq T^{total}$. Set
$\mathbb{Z}_{0}=(\mathbb{X}_{0},\mathbb{V}_{0})$

3:Set $t\leftarrow 0$, $T\leftarrow 0$, $m\leftarrow 0$, $\tau_{m}\leftarrow
0$

4:while $T<T^{total}$ do

5: $m\leftarrow m+1$

6: $u\leftarrow\text{Uniform}(0,1)$

7: Solve the equation  $\int_{0}^{\eta_{m}}\lambda_{m}(t)dt=-\log(u),$ (2) to
obtain $\eta_{m}$, where
$\lambda_{m}(t)=\lambda\left(\Phi_{t}(\mathbb{X}_{\tau_{m-1}},\mathbb{V}_{\tau_{m-1}})\right)$,
and $\Phi_{t}(\cdot,\cdot)$ is the flow of the deterministic dynamic.

8: $\tau_{m}\leftarrow\tau_{m-1}+\eta_{m}$, $T\leftarrow\tau_{m}$,
$\mathbb{Z}_{\tau_{m}}=\Phi_{\eta_{m}}(\mathbb{X}_{\tau_{m-1}},\mathbb{V}_{\tau_{m-1}})$,
$\mathbf{Z}_{\tau_{m}}\sim Q(\mathbf{Z}_{\tau_{m}},\cdot)$.

9:Output: A trajectory of the Markov chain up to time $\tau_{M}$,
$\left\\{\mathbb{Z}_{t}\right\\}_{t=0}^{\tau_{M}}$, where
$\mathbb{Z}_{t}=\Phi_{t-\tau_{m}}(\mathbb{Z}_{\tau_{m}})$ for $\tau_{m}\leq
t<\tau_{m+1}$.

In many cases, evaluating the path integral
$\int_{t=0}^{T}h(\mathbb{Z}_{t})dt$ may however be expensive, or even
impossible, and a discretization of the simulated trajectory is a feasible
alternative. This means estimating the quantity of interest, $I$, by the
following estimator

$\hat{I}=\frac{1}{N}\sum_{n=1}^{N}h(\mathbb{X}_{\frac{nT}{N}}).$

In practice, the main difficulty in implementing a PDMP-based sampler is the
generation of the occurrence times of the associated time-inhomogeneous
Poisson process with event rate $\lambda(\cdot)$. Fortunately, the following
two theorems alleviate this difficulty.

######  Theorem 2.1 (Superposition Theorem)

(Kingman,, 1992) Let $\Pi_{1},\Pi_{2},\cdots,$ be a countable collection of
independent Poisson processes on state space $\mathbb{R}^{+}$ and let
$\Pi_{n}$ have rate $\lambda_{n}(\cdot)$ for each $n$. If
${\sum_{n=1}^{\infty}}\lambda_{n}(t)<\infty$ for all $t$, then the
superposition

$\Pi=\bigcup_{n=1}^{\infty}\Pi_{n}$

is a Poisson process with rate

$\lambda(t)=\sum_{n=1}^{\infty}\lambda_{n}(t)$

######  Theorem 2.2 (Thinning Theorem)

(Lewis and Shedler,, 1979) Let
$\lambda:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ and
$\Lambda:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ be continuous functions such
that $\lambda(t)\leq\Lambda(t)$ for all $t\geq 0$. Let
$\tau^{1},\tau^{2},\cdots,$ be the increasing finite or infinite sequence of a
Poisson process with rate $\Lambda(\cdot)$. For all $i$, if the point
$\tau^{i}$ is removed from the sequence with probability
$1-\lambda(t)/\Lambda(t)$, then the remaining points
$\tilde{\tau}^{1},\tilde{\tau}^{2},\cdots$ form a non-homogeneous Poisson
process with rate $\lambda(\cdot)$.

In practice, according to Theorem 1, we can split the event rate function into
the summation of several event sub-rate functions and take the minimum of the
first arrival times of the Poisson processes, induced by these sub-rate
functions, as the desired event time interval. In addition, in order to
generate the first arrival times of the sub-Poisson processes, we can choose
an upper bound function, whose induced Poisson process is easy to simulate,
for each sub-rate function, and resort to Theorem 2.

###  2.3 Two reference PDMD-based samplers

Almost all existing PDMD-based samplers are based on two specific versions,
both of which rely on linearly deterministic dynamics, a feature that
facilitates the determination of the state of the Markov chain between Poisson
events. Vanetti et al., (2017) uses Hamiltonian dynamics over an approximation
of the target distribution to accelerate the bouncy particle sampler, but the
efficiency of that modification depends on the quality of the approximation
and it only transfers the difficulty from setting the deterministic dynamics
to computing the event rate function.

####  2.3.1 Bouncy Particle sampler

For the Bouncy Particle sampler, as described by Bouchard-Côté et al., (2018),
the velocity set $\mathcal{V}$ is either the Euclidean space $\mathbb{R}^{d}$,
or the unit sphere $\mathbb{S}^{d-1}$. The associated augmented target
distribution is either
$\rho(d\mathbb{x},d\mathbb{v})=\pi(d\mathbb{x})\mathcal{N}(d\mathbb{v}|0,I_{d})$,
or
$\rho(d\mathbb{x},d\mathbb{v})=\pi(d\mathbb{x})\mathcal{U}_{\mathbb{S}^{d-1}}(d\mathbb{v})$,
where $\mathcal{N}(\cdot|0,I_{d})$ represents the standard $d$-dimensional
Gaussian distribution and $\mathcal{U}_{\mathbb{S}^{d-1}}(d\mathbb{v})$
denotes the uniform distribution over $\mathbb{S}^{d-1}$, respectively. The
corresponding deterministic dynamic is

$\frac{d\mathbb{X}_{t}}{dt}=\mathbb{V}_{t},\quad\frac{d\mathbb{V}_{t}}{dt}=\mathbb{0},$

the event rate satisfies
$\lambda(\mathbb{z})=\lambda(\mathbb{x},\mathbb{v})=\langle\mathbb{v},\nabla
U(\mathbb{x})\rangle_{+}+\lambda^{\text{ref}}$, where $\lambda^{\text{ref}}$
is a user-chosen non-negative constant and the transition dynamic is as

$Q((d\mathbb{x}^{\prime},d\mathbb{v}^{\prime})|(\mathbb{x},\mathbb{v}))=\frac{\langle\mathbb{v},\nabla
U(\mathbb{x})\rangle_{+}}{\lambda(\mathbb{x},\mathbb{v})}\delta_{\mathbb{x}}(d\mathbb{x}^{\prime})\delta_{R_{\nabla
U(\mathbb{x})}\mathbb{v}}(d\mathbb{v}^{\prime})+\frac{\lambda^{\text{ref}}}{\lambda(\mathbb{x},\mathbb{v})}\delta_{\mathbb{x}}(d\mathbb{x}^{\prime})\varphi(d\mathbb{v}^{\prime})$

where $\varphi(d\mathbb{v})=\mathcal{U}_{\mathbb{S}^{d-1}}(d\mathbb{v})$ or
$\varphi(d\mathbb{v})=\mathcal{N}(d\mathbb{v}|0,I_{d})$, depending on the
choice of the velocity set, and the operator $R_{\mathbb{w}}$, for any non-
zero vector $\mathbb{w}\in\mathbb{R}^{d}-\\{\mathbb{0}\\}$, is
$R_{\mathbb{w}}\mathbb{v}=\mathbb{v}-2\frac{\langle\mathbb{w},\mathbb{v}\rangle}{\langle\mathbb{w},\mathbb{w}\rangle}\mathbb{w}$.

####  2.3.2 Zigzag sampler

For the Zigzag sampler (Bierkens et al.,, 2016), the velocity set,
$\mathcal{V}$, is the discrete set
$\\{\mathbb{v}=(v_{1},\cdots,v_{d})|v_{i}\in\\{-1,1\\},i=1,\cdots,d\\}$ and
$\rho(d\mathbb{x},d\mathbb{v})=\pi(d\mathbb{x})\varphi(d\mathbb{v})$, where
$\varphi$ is the uniform distribution over $\mathcal{V}$. ZS uses the same
linear deterministic dynamics as BPS. Its event rate is
$\lambda(\mathbb{z})=\sum_{i=1}^{d}\lambda_{i}(\mathbb{x},\mathbb{v})=\sum_{i=1}^{d}\left[\left\\{v_{i}\nabla_{i}U(\mathbb{x})\right\\}_{+}+\lambda^{\text{ref}}_{i}\right]$,
where the $\lambda^{\text{ref}}_{i}$’s are user-chosen non-negative constants.
The transition dynamics is

$Q((d\mathbb{x}^{\prime},d\mathbb{v}^{\prime})|(\mathbb{x},\mathbb{v}))=\sum_{i=1}^{d}\frac{\lambda_{i}(\mathbb{x},\mathbb{v})}{\lambda(\mathbb{x},\mathbb{v})}\delta_{\mathbb{x}}(d\mathbb{x}^{\prime})\delta_{F_{i}\mathbb{v}}(d\mathbb{v}^{\prime}),$

where $F_{i}$ denotes the operator that flips the $i$-th component of
$\mathbb{v}$ and keeps the others unchanged. In practice, ZS relies on the
Superposition Theorem: At each event time, ZS simulates $d$ Poisson processes,
with rates $\lambda_{i}(\mathbb{x}+t\mathbb{v},\mathbb{v})$, computes their
first occurrence time, and takes their minimum, e.g., the $i$-th, for the
duration between current and next events, and flips the $i$-th component of
the velocity $\mathbb{v}$.

##  3 Coordinate sampler

We now describe the coordinate sampler (CS), in which only one component of
$\mathbb{x}$ evolves and the others remain inactive between event times. For
CS, the velocity set $\mathcal{V}$ is chosen to be $\\{\pm
e_{i},i=1,\cdots,d\\}$, where $e_{i}$ is the vector with $i$-th component
equal to one and the others set to zero. The augmented target distribution is
$\rho(d\mathbb{x},d\mathbb{v})=\pi(d\mathbb{x})\varphi(d\mathbb{v})$, with
$\varphi(d\mathbb{v})$ the uniform distribution over $\mathcal{V}$. The PDMP
characteristics of CS are thus

  1. 1.

Deterministic dynamic:

$\frac{d\mathbb{X}_{t}}{dt}=\mathbb{V}_{t},\quad\frac{d\mathbb{V}_{t}}{dt}=\mathbb{0}\,.$

  2. 2.

Event occurrence rate: $\lambda(\mathbb{z})=\langle\mathbb{v},\nabla
U(\mathbb{x})\rangle_{+}+\lambda^{\text{ref}}$, where $\lambda^{\text{ref}}$
is a user-chosen positive constant.

  3. 3.

Transition dynamic:

$Q((d\mathbb{x}^{\prime},d\mathbb{v}^{\prime})|(\mathbb{x},\mathbb{v}))=\sum_{\mathbb{v}^{*}\in\mathcal{V}}\frac{\lambda(\mathbb{x},-\mathbb{v}^{*})}{\lambda(\mathbb{x})}\delta_{\mathbb{x}}(d\mathbb{x}^{\prime})\delta_{\mathbb{v}^{*}}(d\mathbb{v}^{\prime})$

where
$\lambda(\mathbb{x})=\sum_{\mathbb{v}\in\mathcal{V}}\lambda(\mathbb{x},\mathbb{v})=2d\lambda^{\text{ref}}+\sum_{i=1}^{d}\big{|}\frac{\partial
U(\mathbb{x})}{\partial x_{i}}\big{|}$,

which translates into the pseudo-code

Algorithm 2 Coordinate Sampler

1:Input: Start with position $\mathbb{X}_{0}$, velocity $\mathbb{V}_{0}$ and
set simulation time threshold $T^{total}$.

2:Generate a set of event times $\\{\tau_{0},\tau_{1},\cdots,\tau_{M}\\}$ and
their associated states
$\\{\mathbb{Z}_{\tau_{0}},\mathbb{Z}_{\tau_{1}},\cdots,\mathbb{Z}_{\tau_{M}}\\}$,
where $\tau_{0}=0$, $\tau_{M}\geq T^{total}$ and $\tau_{M-1}<T^{total}$,
$\mathbb{Z}_{0}=(\mathbb{X}_{0},\mathbb{V}_{0})$

3:Set $t\leftarrow 0$, $T\leftarrow 0$, $m\leftarrow 0$, $\tau_{m}\leftarrow
0$

4:while $T<T^{total}$ do

5: $m\leftarrow m+1$

6: $u\leftarrow\text{Uniform}(0,1)$

7: Solve the equation  $\int_{0}^{\eta_{m}}\lambda_{m}(t)dt=-\log(u)$ with
respect to $\eta_{m}$, where
$\lambda_{m}(t)=\lambda\left(\mathbb{X}_{\tau_{m-1}}+t\mathbb{V}_{\tau_{m-1}},\mathbb{v}_{\tau_{m-1}}\right)$.

8: $\tau_{m}\leftarrow\tau_{m-1}+\eta_{m}$, $T\leftarrow\tau_{m}$,
$\mathbb{Z}_{\tau_{m}}=(\mathbb{X}_{\tau_{m-1}}+\eta_{m}\mathbb{V}_{\tau_{m-1}},\mathbb{v}_{\tau_{m-1}})$,
$\mathbf{Z}_{\tau_{m}}\sim Q(\mathbf{Z}_{\tau_{m}},\cdot)$.

9:Output: A trajectory of the Markov chain over $[0,\tau_{M}]$,
$\left\\{\mathbb{Z}_{t}\right\\}_{t=0}^{\tau_{M}}$, where
$\mathbb{Z}_{t}=(\mathbb{X}_{\tau_{m}}+(t-\tau_{m})\mathbb{V}_{\tau_{m}},\mathbb{V}_{\tau_{m}})$
for $\tau_{m}\leq t<\tau_{m+1}$.

###  3.1 Theoretical properties of the coordinate sampler

We now establish that CS is associated with the augmented target distribution,
$\rho(d\mathbb{x},d\mathbb{v})$, as its invariant distribution under the
condition that $U:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ is $C^{1}$.
Furthermore, under the following assumptions, the Markov process induced by CS
is $V$-uniformly ergodic for the Lyapunov function

$V(\mathbb{x},\mathbb{v})={e^{U(\mathbb{x})/2}}/{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle_{+}}},$

which was also used in Deligiannidis et al., (2017).

###### Theorem 3.1

For any positive $\lambda^{\text{ref}}>0$, the PDMP produced by CS enjoys
$\rho(d\mathbb{x},d\mathbb{v})$ as its unique invariant distribution, provided
the potential $U$ is $C^{1}$.

It is easy to check that the generator of CS, $\mathcal{L}$, satisfies

$\int\mathcal{L}f(\mathbb{z})\rho(d\mathbb{z})=0,$

for all functions $f$ in its extended generator, which means that $\rho$ is an
invariant distribution of CS (Davis,, 1993, Theorem 34.7). Uniqueness follows
from the positivity of $\lambda^{\text{ref}}$, which enables the Markov
process to reach any state $(\mathbb{x}^{*},\mathbb{v}^{*})$ from any starting
state $(\mathbb{x}_{0},\mathbb{v}_{0})$, in finite time. (Details of the proof
are provided as supplementary material.) In practice, it appears that the
constraint $\lambda^{\text{ref}}>0$ is unnecessary for convergence in many
examples.

Assumptions: Assume $U:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ satisfy the
following conditions, reproduced from Deligiannidis et al., (2017),

  1. A.1

$\frac{\partial^{2}U(\mathbb{x})}{\partial x_{i}x_{j}}$ is locally Lipschitz
continuous for all $i,j$,

  2. A.2

$\int\big{|}\nabla U(\mathbb{x})\big{|}\pi(d\mathbb{x})<\infty$,

  3. A.3

$\underline{\lim}_{|\mathbb{x}|\rightarrow\infty}{e^{U(\mathbb{x})/2}}/\sqrt{\big{|}\nabla
U(\mathbb{x})\big{|}}>0$

  4. A.4

$V\geq c_{0}$ for some positive constant $c_{0}$.

Conditions: We set conditions

  1. C.1

$\displaystyle{\underline{\lim}_{|x|\to\infty}}\big{|}\nabla
U(x)\big{|}=\infty$, $\displaystyle{\overline{\lim}_{|x|\to\infty}}\|\Delta
U(x)\|\leq\alpha_{1}<\infty$ and $\lambda^{\text{ref}}>\sqrt{8\alpha_{1}}$.

  2. C.2

$\underline{\lim}_{|\mathbb{x}|\rightarrow\infty}\big{|}\nabla
U(\mathbb{x})\big{|}=2\alpha_{2}>0$,
$\displaystyle{\overline{\lim}_{|x|\to\infty}}\|\Delta U(x)\|=0$ and
$\lambda^{\text{ref}}<\frac{\alpha_{2}}{14d}$.

where C.1 corresponds to distributions whose tails decay at rate
$\mathcal{O}(|\mathbb{x}|^{\beta})$, where $1<\beta\leq 2$, and C.2 to
distributions with tails of order $\mathcal{O}(|\mathbb{x}|^{1})$.

###### Theorem 3.2

Suppose assumptions $A.1-A.4$ hold, $\lambda^{\text{ref}}>0$, and one of the
conditions C.1 or C.2 holds, then CS is $V$-uniformly ergodic: There exist
constants $\Gamma<\infty$ and $0<\gamma<1$, such that

$\|P^{t}(\mathbb{z},\cdot)-\rho\|_{V}\leq V(\mathbb{z})\Gamma\gamma^{t}\,,$

where $P^{t}(\mathbb{z},\cdot)$ is the distribution of the Markov chain with
starting state $\mathbb{z}$ at time $t$, and the norm $\|\cdot\|_{V}$ is
defined by

$\|\mu\|_{V}=\sup_{|f|<V}\int f(\mathbb{z})\mu(d\mathbb{z})\,.$

The proof appears in the supplementary material, based on techniques quite
similar to those in Deligiannidis et al., (2017).

###  3.2 An informal comparison between Zigzag and coordinate samplers

For CS, each event time sees a change of a single component of $\mathbb{X}$,
in contrast with ZS, which modifies all components at the same time. At first
this gives the impression that CS is less efficient than ZS in its exploration
of the target space, because of this restriction. However, this intuition is
misleading: Suppose that the $\lambda_{i}$’s, $i=1,\cdots,d$ in ZS and
$\lambda$ in CS are of a similar scale, for instance taking the expected
duration between two Poisson events to be the same value $\ell$. Assume
further that computing an occurrence time have the same computation cost, $c$,
for all Poisson processes. In ZS, the event rate is the summation of the rates
$\lambda_{1},\cdots,\lambda_{d}$. Therefore, the time duration between two
events is $\frac{\ell}{d}$ and the induced computation cost is $dc$. Thus,
that each component of $\mathbb{X}$ evolves for a time duration $\ell$ costs
$d^{2}c$ for ZS. By contrast, in CS, a $dc$ computation cost will result from
the Markov chain moving for a duration time $d\ell$. Hence, the computation
cost for monitoring each component for a time duration $\ell$ is also $dc$. As
a result, CS is $\mathcal{O}(d)$ times more efficient than ZS in terms of the
evolution of a given component of $\mathbb{X}$.

##  4 Numerical experiments

In this section, we compare the efficiency of both samplers over benchmarks (a
banana-shaped distribution, two multivariate Gaussian distributions, and a
Bayesian logistic model). In each model, we run both samplers for the same
computer time or the same number of calls of the event rate functions and we
compare their efficiency in terms of effective sample size (ESS) (Liu,, 2008)
per second or per call of an event rate function. The models are reproduced
fourty times to produce an averaged efficiency ratio, namely the ratio of ESS
per second for CS over the one for ZS. We use the function ess of package
mcmcse in R to compute ESS of samples. In the first three experiments, we use
canonical ZS and canonical CS, meaning that
$\lambda_{i}^{\text{ref}}=0,i=1,\cdots,d$ in ZS and $\lambda^{\text{ref}}=0$
in CS, since such settings guarantee ergodicity. For the Bayesian logistic
model, we set $\lambda_{i}^{\text{ref}}=1,i=1,\cdots,d$ and
$\lambda^{\text{ref}}=1$. For the log-Gaussian Cox point process, we set
$\lambda_{i}^{\text{ref}}=0.1,i=1,\cdots,d$ and $\lambda^{\text{ref}}=0.1$ to
achieve a 10$\%$ complete refreshment of velocity.

Banana-Shaped Distribution: The target distribution is a $2$-dimensional
banana-shaped distribution with density

$\pi(\mathbb{x})\propto\exp\left\\{-(x_{1}-1)^{2}-\kappa(x_{2}-x_{1}^{2})^{2}\right\\}$

where $\kappa$ controls the similarity between $x_{2}$ and $x_{1}^{2}$. A high
$\kappa$ enforces the approximate constraint $x_{2}\simeq x_{1}^{2}$. The
comparison between Zigzag and coordinate samplers runs over the configurations
$2^{-2}\leq\kappa\leq 2^{5}$. With an increase in $\kappa$, the distribution
becomes more difficult to simulate and the event rate functions in CS and ZS
make the generations of time durations more costly. Figure 1 shows that CS is
more efficient than ZS across a large range of $\kappa$ in this model.

Strongly Correlated Multivariate Gaussian Distribution (MVN1): Here, the
target is a multivariate Gaussian distribution with zero mean and covariance
matrix equal to $A\in\mathbb{R}^{d\times d}$, where $A_{ii}=1$ and
$A_{ij}=0.9,i\neq j$. We consider the values $d=10,20,\ldots,100$ in our
comparison of the sampling methods.

Correlated Multivariate Gaussian Distribution (MVN2): In this scenario, the
target distribution is again a multivariate Gaussian distribution with zero
mean and covariance matrix such that $A_{ii}=1$ and $A_{ij}=0.9^{|{i-j}|}$.
Once again, the comparison runs for $d=10,\ldots,100$.

Figure 2 presents the comparison between CS and ZS for both models MVN1 and
MVN2 in terms of the minimal ESS, mean ESS, median ESS and maximal ESS taken
across all $d$ components per generation of occurrence time induced by event
(sub-) rate function. In both models, the efficiency ratio and thus the
improvement brought by CS over ZS increases with the dimension $d$.

In Table 1, we further compare CS with several standard MCMC algorithms for a
20-dimensional MVN2 model in terms of Kolmogorov-Smirnov statistic (KS) to the
target. Since it is infeasible to compute such quantities for multivariate
distributions, we compute marginal distances between samples from each
algorithm and from the target, across coordinates, and take the minimum, mean,
median and maximum of these as a summary of the efficiency of each algorithm,
for identical computation times about $155$ seconds. In this experiment, HMC
performs best in terms of Kolmogorov-Smirnov statistic. However, among the
PDMP-based MCMC algorithms, CS exceeds ZS and BPS.

Table 1: Comparison for a 20-dimensional MVN2 model based on 40 repetitions, in terms of minimum, mean, median and maximum of the marginal distances across the components for each criterion. The smaller the numerical value, the better the algorithm performs. Sampler | Min KS | Mean KS | Median KS | Max KS
---|---|---|---|---
CS | $4.02\times 10^{-3}$ | $7.00\times 10^{-3}$ | $7.04\times 10^{-3}$ | $10.02\times 10^{-3}$
ZS | $9.17\times 10^{-3}$ | $16.53\times 10^{-3}$ | $16.07\times 10^{-3}$ | $24.34\times 10^{-3}$
BPS | $4.94\times 10^{-3}$ | $8.98\times 10^{-3}$ | $9.04\times 10^{-3}$ | $12.58\times 10^{-3}$
HMC | $1.26\times 10^{-3}$ | $2.12\times 10^{-3}$ | $2.04\times 10^{-3}$ | $3.31\times 10^{-3}$

Bayesian Logistic Model: In this example, the target is the posterior of a
Bayesian logistic model under a flat prior, with no intercept. The simulated
dataset contains $N$ observations $\\{(\mathbb{r}_{n},t_{n})\\}_{n=1}^{N}$,
where each $r_{n,i}$, $n=1,\cdots,N,i=1,\cdots,d$, is drawn from a standard
normal distribution and $t_{n}$ is drawn from $\\{0,1\\}$ uniformly. The
targeted density function is thus

$\pi(\mathbb{x})\propto\prod_{n=1}^{N}\frac{\exp(t_{n}\mathbb{x}^{T}\mathbb{r}_{n})}{1+\exp(\mathbb{x}^{T}\mathbb{r}_{n})}$

In the simulations, we set $N=40,d=10$, and $\lambda^{\text{ref}}=1$ for CS,
and $\lambda_{i}^{\text{ref}}=1,i=1,\cdots,d$ for ZS. Figure 3 presents the
comparison between the two samplers, with a massive improvement brought by our
proposal.

Log-Gaussian Cox Point Process In this example, already implemented by
Galbraith, (2016), the observations $\mathbf{Y}=\\{y_{ij}\\}$ are Poisson
distributed and conditionally independent given a latent intensity process
$\mathbf{\Lambda}=\\{\lambda_{ij}\\}$ with means
$s\lambda_{ij}=s\exp(x_{ij})$, where $s=1/d^{2}$. The underlying process
$\mathbf{X}=\\{x_{ij}\\}$ is a Gaussian process with mean function
$m(x_{ij})=\mu\mathbf{1}$ and covariance function
$\Sigma(x_{i,j},x_{i^{\prime},j^{\prime}})=\sigma^{2}\exp(-\delta(i,i^{\prime},j,j^{\prime})/(\beta
d))$, where
$\delta(i,i^{\prime},j,j^{\prime})^{2}={(i-i^{\prime})^{2}+(j-j^{\prime})^{2}}$.
In our experiment, we set $d=20$ and choose
$\sigma^{2}=1.91,\mu=\log(126)-\sigma^{2}/2$ and $\beta=1/6$. The target is
conditional on the observations $\mathbf{Y}$,

$\pi(\mathbf{X}|\mathbf{Y},\mu,\sigma,\beta)\propto\exp\left\\{\sum_{i,j=1}^{d}(y_{ij}x_{ij}-s\exp(x_{ij}))-\frac{1}{2}(\mathbf{X}-\mu\mathbf{1})^{T}\Sigma^{-1}(\mathbf{X}-\mu\mathbf{1})\right\\}$

We run CS and ZS for 160 seconds each and obtain about $1,000$ draws from each
sampler. Figure 4 shows the first two components of the samples generated by
both samplers. In Figure 5, left, the values of the log-densities explored by
CS (red) are more diverse than those visited by the ZZ (blue), while the right
graph shows a similar pattern for the last component of the generated samples.
As also shown in the raw plots of Figure 4, CS is thus more efficient than ZS
in exploring the target distribution.

##  5 Conclusion

We have introduced and studied the coordinate sampler as an alternative to the
Zigzag sampler of Bierkens et al., (2016) and compared the efficiencies of the
two samplers in terms of effective sample size over several simulation
experiments. In all examples, CS exhibits a higher efficiency, which gain
increases with the dimension of the target distribution, while enjoying the
same ergodicity guarantees. While our intuition about the advantage of a
component-wise implementation led to our proposal, exhibiting a theoretical
reason for this improvement requires further investigation.

We also stress that, among PDMP-based MCMC samplers, CS is quite easy to scale
for big data problems, as is the Zigzag sampler. In addition, taking advantage
of the techniques exposed in Bierkens et al., (2018), CS can also be
implemented for distributions defined on restricted domains. In such settings,
since only one component of the target distribution is active between Poisson
events, the efficiency of CS relatively to ZS may suffer, especially in cases
when the variances across the components are of different magnitudes. An
appropriate reparametrization of the target distribution should however
alleviate this problem, and accelerate CS, which amounts to a pre-conditioning
of the velocity set. An interesting extension that needs further investigation
is to build CS that take advantage of the curvature of the target by Riemann
manifold techniques as in Girolami and Calderhead, (2011).

Figure 1: Banana-Shaped Distribution: the $x$-axis is indexed by
$\log(\kappa)$, the $y$-axis corresponds to the ratio of ESS per second of
coordinate versus Zigzag samplers. The red line shows the efficiency ratio for
the first component, the blue line for the second component, and the green
line for the log probability.

Figure 2: The upper plot shows the results for MVN1 and the lower for MVN2.
The $x$-axis indexes the dimension $d$ of the distribution, and the $y$-axis
the efficiency ratios of CS over ZS in terms of minimum, mean, median and
maximum of ESS across the components over the number of recall event rate
function. Figure 3: Comparison of CS versus ZS for the Bayesian logistic
model: the $y$-axis represents the ESS per second. Figure 4: Samples generated
by CS and ZS samplers for the same computation time, when targeting a log-
Gaussian Cox point process. Only the first two components are represented
here. Figure 5: The plots of log-density and the final component of the
samples generate by CS (red) and ZS (blue) for the same computation time, when
targeting a log-Gaussian Cox point process.

## References

  * Bierkens et al., (2018) Bierkens, J., Bouchard-Côté, A., Doucet, A., Duncan, A. B., Fearnhead, P., Lienart, T., Roberts, G., and Vollmer, S. J. (2018).  Piecewise deterministic Markov processes for scalable Monte Carlo on restricted domains.  Statistics & Probability Letters, 136:148–154.
  * Bierkens et al., (2016) Bierkens, J., Fearnhead, P., and Roberts, G. (2016).  The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data.  arXiv preprint arXiv:1607.03188.
  * Bierkens et al., (2017) Bierkens, J., Roberts, G., et al. (2017).  A piecewise deterministic scaling limit of lifted Metropolis–Hastings in the Curie–Weiss model.  The Annals of Applied Probability, 27(2):846–882.
  * Bouchard-Côté et al., (2018) Bouchard-Côté, A., Vollmer, S. J., and Doucet, A. (2018).  The bouncy particle sampler: A nonreversible rejection-free Markov chain Monte Carlo method.  Journal of the American Statistical Association, pages 1–13.
  * Davis, (1984) Davis, M. H. (1984).  Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models.  Journal of the Royal Statistical Society. Series B (Methodological), pages 353–388.
  * Davis, (1993) Davis, M. H. (1993).  Markov Models & Optimization, volume 49.  CRC Press.
  * Deligiannidis et al., (2017) Deligiannidis, G., Bouchard-Côté, A., and Doucet, A. (2017).  Exponential Ergodicity of the Bouncy Particle Sampler.  arXiv preprint arXiv:1705.04579.
  * Down et al., (1995) Down, D., Meyn, S. P., and Tweedie, R. L. (1995).  Exponential and uniform ergodicity of markov processes.  The Annals of Probability, 23(4):1671–1691.
  * Fearnhead et al., (2018) Fearnhead, P., Bierkens, J., Pollock, M., Roberts, G. O., et al. (2018).  Piecewise deterministic Markov processes for continuous-time Monte Carlo.  Statistical Science, 33(3):386–412.
  * Fontbona et al., (2016) Fontbona, J., Guérin, H., and Malrieu, F. (2016).  Long time behavior of telegraph processes under convex potentials.  Stochastic Processes and their Applications, 126(10):3077–3101.
  * Galbraith, (2016) Galbraith, N. (2016).  On event-chain Monte Carlo methods.  Master’s thesis, Department of Statistics, Oxford University, 9\.
  * Girolami and Calderhead, (2011) Girolami, M. and Calderhead, B. (2011).  Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214.
  * Glynn and Haas, (2006) Glynn, P. W. and Haas, P. J. (2006).  Laws of large numbers and functional central limit theorems for generalized semi-Markov processes.  Stochastic Models, 22(2):201–231.
  * Harland et al., (2017) Harland, J., Michel, M., Kampmann, T. A., and Kierfeld, J. (2017).  Event-chain Monte Carlo algorithms for three-and many-particle interactions.  EPL (Europhysics Letters), 117(3):30001.
  * Hastings, (1970) Hastings, W. K. (1970).  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57(1):97–109.
  * Kingman, (1992) Kingman, J. F. C. (1992).  Poisson processes, volume 3.  Clarendon Press.
  * Lewis and Shedler, (1979) Lewis, P. A. and Shedler, G. S. (1979).  Simulation of nonhomogeneous Poisson processes by thinning.  Naval Research Logistics (NRL), 26(3):403–413.
  * Liu, (2008) Liu, J. S. (2008).  Monte Carlo strategies in scientific computing.  Springer Science & Business Media.
  * 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.  The journal of Chemical Physics, 21(6):1087–1092.
  * Michel et al., (2014) Michel, M., Kapfer, S. C., and Krauth, W. (2014).  Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps.  The Journal of Chemical Physics, 140(5):054116.
  * Michel and Sénécal, (2017) Michel, M. and Sénécal, S. (2017).  Forward Event-Chain Monte Carlo: a general rejection-free and irreversible Markov chain simulation method.  arXiv preprint arXiv:1702.08397.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  MCMC using Hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11).
  * Pakman et al., (2016) Pakman, A., Gilboa, D., Carlson, D., and Paninski, L. (2016).  Stochastic bouncy particle sampler.  arXiv preprint arXiv:1609.00770.
  * Peters et al., (2012) Peters, E. A. et al. (2012).  Rejection-free Monte Carlo sampling for general potentials.  Physical Review E, 85(2):026703.
  * Robert and Casella, (2004) Robert, C. P. and Casella, G. (2004).  Monte Carlo statistical methods.  Springer New York.
  * Sherlock and Thiery, (2017) Sherlock, C. and Thiery, A. H. (2017).  A Discrete Bouncy Particle Sampler.  arXiv preprint arXiv:1707.05200.
  * Vanetti et al., (2017) Vanetti, P., Bouchard-Côté, A., Deligiannidis, G., and Doucet, A. (2017).  Piecewise Deterministic Markov Chain Monte Carlo.  arXiv preprint arXiv:1707.05296.

Supplementary material

##  6 Proof of Theorem 3

To establish this result, we need both following lemmas:

###### Lemma 1

$\rho(d\mathbb{x},\mathbb{v})=\pi(d\mathbb{x})\varphi(d\mathbb{v})$ is the
invariant distribution of the Markov process induced by the Coordinate
Sampler.

Proof of Lemma 1: The generator of the Markov process induced by the
Coordinate Sampler is, by Davis (1993, Theorem 26.14),

$\begin{split}\mathcal{L}f(\mathbb{x},\mathbb{v})&=\left\langle\frac{\partial
f(\mathbb{x},\mathbb{v})}{\partial\mathbb{x}},\mathbb{v}\right\rangle+\lambda(\mathbb{x},\mathbb{v})\int_{\mathbb{x}^{\prime}}\int_{\mathbb{v}^{\prime}}\left(f(\mathbb{x}^{\prime},\mathbb{v}^{\prime})-f(\mathbb{x},\mathbb{v})\right)Q((\mathbb{x},\mathbb{v}),(d\mathbb{x}^{\prime},d\mathbb{v}^{\prime}))\\\
&=\left\langle\frac{\partial
f(\mathbb{x},\mathbb{v})}{\partial\mathbb{x}},\mathbb{v}\right\rangle+\lambda(\mathbb{x},\mathbb{v})\sum_{i=1}^{2d}\frac{\lambda(\mathbb{x},-\mathbb{v}_{i})}{\lambda(\mathbb{x})}f(\mathbb{x},\mathbb{v}_{i})-\lambda(\mathbb{x},\mathbb{v})f(\mathbb{x},\mathbb{v})\end{split}$

Since we have

$\begin{split}\lambda(\mathbb{x},-\mathbb{v})-\lambda(\mathbb{x},\mathbb{v})&=\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle_{+}-\lambda^{\text{ref}}-\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle_{+}\\\ &=\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle\end{split}$

as a result,
$\int_{\mathbb{x}}\int_{\mathbb{v}}\mathcal{L}f(\mathbb{x},\mathbb{v})\pi(\mathbb{x})\varphi(\mathbb{v})d\mathbb{x}d\mathbb{v}=0$,
for all $f\in\mathcal{D}(\mathcal{L})$ ($\mathcal{D}(\mathcal{L})$ is defined
in Deligiannidis et al.,, 2017, Section 2.1). That is,

$\displaystyle\int_{\mathbb{x}}\int_{\mathbb{v}_{i}}\mathcal{L}f(\mathbb{x},\mathbb{v})\pi(\mathbb{x})\varphi(\mathbb{v})d\mathbb{x}d\mathbb{v}$
$\displaystyle=\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\mathcal{A}f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle=\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\left\langle\frac{\partial
f(\mathbb{x},\mathbb{v}_{i})}{\partial\mathbb{x}},\mathbb{v}_{i}\right\rangle\pi(\mathbb{x})d\mathbb{x}+\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\lambda(\mathbb{x},\mathbb{v}_{i})\sum_{j=1}^{2d}\frac{\lambda(\mathbb{x},-\mathbb{v}_{j})}{\lambda(\mathbb{x})}f(\mathbb{x},\mathbb{v}_{j})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle\quad-\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\lambda(\mathbb{x},\mathbb{v}_{i})f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle=\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\left\langle-\nabla
U(\mathbb{x}),\mathbb{v}_{i}\right\rangle
f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}+\frac{1}{2d}\sum_{j=1}^{2d}\int_{\mathbb{x}}\left(\sum_{i=1}^{2d}\lambda(\mathbb{x},\mathbb{v}_{i})\right)\frac{\lambda(\mathbb{x},-\mathbb{v}_{j})}{\lambda(\mathbb{x})}f(\mathbb{x},\mathbb{v}_{j})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle\quad-\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\lambda(\mathbb{x},\mathbb{v}_{i})f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle=\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\left\langle-\nabla
U(\mathbb{x}),\mathbb{v}_{i}\right\rangle
f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}+\frac{1}{2d}\sum_{j=1}^{2d}\int_{\mathbb{x}}\lambda(\mathbb{x},-\mathbb{v}_{j})f(\mathbb{x},\mathbb{v}_{j})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle\quad-\frac{1}{2d}\sum_{i=1}^{2d}\int_{\mathbb{x}}\lambda(\mathbb{x},\mathbb{v}_{i})f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{2d}\int_{\mathbb{x}}\left\\{\left\langle-\nabla
U(\mathbb{x}),\mathbb{v}_{i}\right\rangle+\lambda(\mathbb{x},-\mathbb{v}_{i})-\lambda(\mathbb{x},\mathbb{v}_{i})\right\\}f(\mathbb{x},\mathbb{v}_{i})\pi(\mathbb{x})d\mathbb{x}=0$

Following Davis (1993, Theorem 34.7), $\rho$ is the invariant distribution of
the Markov chain induced by Coordinate Sampler.
The following lemma is the same as Lemma 2 of Deligiannidis et al., (2017) and
we include some details of the proof of irreducibility for the Markov process
induced by the Coordinate Sampler since some details are not identical in the
original proof. In Deligiannidis et al., (2017), two events and transition
dynamics will ensure the path reach any desired state, say $\mathbb{z}^{*}$.
However, in our case, $d+2$ times are required, which makes the proof more
complicated and we further resort to a Dirichlet distribution, instead of a
uniform distribution. For simplicity, we represent the expectation over
$\varphi$ in the form of an integral, instead of a summation.

###### Lemma 2

For all $T>0$,
$\lambda^{\text{ref}}>0$,$\mathbb{z}_{0}=(\mathbb{x}_{0},\mathbb{v}_{0})\in
B(0,\frac{T}{6})\times\mathcal{V}$, and a Borel set $A\subset
B(0,\frac{T}{6})\times\mathcal{V}$

$\mathbb{P}(\mathbb{z}_{0},\mathbb{Z}_{T}\in A)\geq
C(T,d,\lambda^{\text{ref}})\int\int_{A}\varphi(d\mathbb{v})d\mathbb{x}$

for some constant $C>0$ depending only on $T,d,\lambda^{\text{ref}}$. Hence,
all compact sets are small and the Markov process induced by Coordinate
Sampler is irreducible.

Proof Let $E$ be the event that there are exactly $d+2$ events during the time
interval $[0,T]$. Suppose $f:B(0,\frac{T}{6})\times\mathcal{V}\to[0,\infty)$
be a bounded, positive function. Then

$\displaystyle\mathbb{E}^{\mathbb{z}}[f(\mathbb{Z}_{T})]\geq\mathbb{E}^{\mathbb{z}}[f(\mathbb{Z}_{T}\mathbb{I}_{E})]$
$\displaystyle=\int_{\mathbb{v}_{1}}\varphi(d\mathbb{v}_{1})\int_{t_{1}=0}^{T}dt_{1}\left(\frac{\lambda(\mathbb{x}_{0}+t_{1}\mathbb{v}_{0},-\mathbb{v}_{1})}{\lambda(\mathbb{x}_{0}+t_{1}\mathbb{v}_{0})}\exp\left\\{-\int_{u_{1}=0}^{t_{1}}\lambda(\mathbb{x}_{0}+u_{1}\mathbb{v}_{0},\mathbb{v}_{0})du_{1}\right\\}\lambda(\mathbb{x}_{0}+t_{1}\mathbb{v}_{0},\mathbb{v}_{0})\right)$
$\displaystyle\times\int_{\mathbb{v}_{2}}\varphi(d\mathbb{v}_{2})\int_{t_{2}=0}^{T-T_{1}}dt_{2}\left(\frac{\lambda(\mathbb{x}_{T_{1}}+t_{2}\mathbb{v}_{1},-\mathbb{v}_{2})}{\lambda(\mathbb{x}_{T_{1}}+t_{2}\mathbb{v}_{1})}\exp\left\\{-\int_{u_{2}=0}^{t_{2}}\lambda(\mathbb{x}_{T_{1}}+u_{2}\mathbb{v}_{1},\mathbb{v}_{1})du_{2}\right\\}\lambda(\mathbb{x}_{T_{1}}+t_{2}\mathbb{v}_{1},\mathbb{v}_{1})\right)$
$\displaystyle\times\cdots$
$\displaystyle\times\int_{\mathbb{v}_{d+2}}\varphi(d\mathbb{v}_{d+2})\int_{t_{d+2}=0}^{T-T_{d+1}}dt_{d+2}\Big{\\{}\Big{\\{}\frac{\lambda(\mathbb{x}_{T_{d+1}}+t_{d+2}\mathbb{v}_{d+1},-\mathbb{v}_{d+2})}{\lambda(\mathbb{x}_{T_{d+1}}+t_{d+2}\mathbb{v}_{d+1})}$
$\displaystyle\times\exp\left\\{-\int_{u_{d+2}=0}^{t_{d+2}}\lambda(\mathbb{x}_{T_{d+1}}+u_{d+2}\mathbb{v}_{d+1},\mathbb{v}_{d+1})du_{d+2}\right\\}\lambda(\mathbb{x}_{T_{d+1}}+t_{d+2}\mathbb{v}_{d+1},\mathbb{v}_{d+1})\Big{\\}}$
$\displaystyle\times\exp\left\\{-\int_{u_{d+3}=0}^{T-T_{d+2}}\lambda(\mathbb{x}_{T_{d+2}}+u_{d+3}\mathbb{v}_{d+2},\mathbb{v}_{d+2})du_{d+3}\right\\}f\left(\mathbb{x}_{T_{d+2}}+(T-T_{d+2})\mathbb{v}_{d+2},\mathbb{v}_{d+2}\right)\Big{\\}}$

where $T_{i}=\sum_{k=1}^{i}t_{k}$ and
$\mathbb{x}_{T_{i}}=\mathbb{x}_{0}+\sum_{k=1}^{i}t_{k}\mathbb{v}_{k-1}$. Since
$\mathbb{x}_{0}\in B(0,T)$, then $\mathbb{x}_{t}\in B(0,2T)$ for all
$t\in[0,T]$. As a result, there exists a constant $K<\infty$, such that

$K\geq\sup_{\mathbb{x}\in B(0,2T)}\Big{|}\nabla U(\mathbb{x})\Big{|}$

Since $\lambda(\mathbb{x},\mathbb{v})\geq\lambda^{\text{ref}}$ for all
$(\mathbb{x},\mathbb{v})\in\mathbb{R}^{d}\times\mathcal{V}$, then

$\displaystyle\mathbb{E}^{\mathbb{z}_{0}}[f(\mathbb{Z}_{T})]$
$\displaystyle\geq\int_{\mathbb{v}_{1}}\varphi(d\mathbb{v}_{1})\int_{t_{1}=0}^{T}dt_{1}\left\\{\frac{\lambda^{\text{ref}}}{\lambda^{\text{ref}}+K}\exp\left\\{-\int_{u_{1}=0}^{t_{1}}\left(\lambda^{\text{ref}}+K\right)du_{1}\right\\}\lambda^{\text{ref}}\right\\}$
$\displaystyle\times\int_{\mathbb{v}_{2}}\varphi(d\mathbb{v}_{2})\int_{t_{2}=0}^{T-T_{1}}dt_{2}\left\\{\frac{\lambda^{\text{ref}}}{\lambda^{\text{ref}}+K}\exp\left\\{-\int_{u_{2}=0}^{t_{2}}\left(\lambda^{\text{ref}}+K\right)du_{2}\right\\}\lambda^{\text{ref}}\right\\}$
$\displaystyle\times\cdots$
$\displaystyle\times\int_{\mathbb{v}_{d+2}}\varphi(d\mathbb{v}_{d+2})\int_{t_{d+2}=0}^{T-T_{d+1}}dt_{d+2}\Big{\\{}\left\\{\frac{\lambda^{\text{ref}}}{\lambda^{\text{ref}}+K}\exp\left\\{-\int_{u_{d+2}=0}^{t_{d+2}}\left(\lambda^{\text{ref}}+K\right)du_{d+2}\right\\}\lambda^{\text{ref}}\right\\}$
$\displaystyle\times\exp\left\\{-\int_{u_{d+3}=0}^{T-T_{d+2}}\left(\lambda^{\text{ref}}+K\right)du_{d+3}\right\\}f\left(\mathbb{x}_{T_{d+2}}+(T-T_{d+2})\mathbb{v}_{d+2},\mathbb{v}_{d+2}\right)\Big{\\}}$
$\displaystyle=\left(\frac{(\lambda^{\text{ref}})^{2}}{\lambda^{\text{ref}}+K}\right)^{d+2}\exp\left\\{-T\left(\lambda^{\text{ref}}+K\right)\right\\}\int_{\mathbb{v}_{1}}\varphi(d\mathbb{v}_{1})\int_{\mathbb{v}_{2}}\varphi(d\mathbb{v}_{2})\cdots\int_{\mathbb{v}_{d+2}}\varphi(d\mathbb{v}_{d+2})$
$\displaystyle\times\int_{t_{1}=0}^{T}dt_{1}\int_{t_{2}=0}^{T-T_{1}}dt_{2}\cdots\int_{t_{d+2}=0}^{T-T_{d+1}}dt_{d+2}f\left(\mathbb{x}_{T_{d+2}}+(T-T_{d+2})\mathbb{v}_{d+2},\mathbb{v}_{d+2}\right)$
$\displaystyle=\left(\frac{(\lambda^{\text{ref}})^{2}}{\lambda^{\text{ref}}+K}\right)^{d+2}\exp\left\\{-T\left(\lambda^{\text{ref}}+K\right)\right\\}\int_{\mathbb{v}_{1}}\varphi(d\mathbb{v}_{1})\int_{\mathbb{v}_{2}}\varphi(d\mathbb{v}_{2})\cdots\int_{\mathbb{v}_{d+2}}\varphi(d\mathbb{v}_{d+2})$
$\displaystyle\times\int_{t=\frac{5T}{6}}^{T}dt\times\int_{r_{1}=0}^{1}dr_{1}\int_{r_{2}=0}^{1-r_{1}}dr_{2}\cdots\int_{r_{d+1}=0}^{1-r_{1}-\cdots-
r_{d}}dr_{d+1}$ $\displaystyle\times
f(\mathbb{x}_{0}+t\sum_{k=1}^{d+1}r_{k}\mathbb{v}_{k-1}+t(1-\sum_{k=1}^{d+1}r_{k})\mathbb{v}_{d+1}+(T-t)\mathbb{v}_{d+2})$

Set $t>\frac{5T}{6}$ and $\mathbb{v}_{d+2}$, then
$\mathbb{x}^{\prime}=\mathbb{x}_{0}+(T-t)\mathbb{v}_{d+2}$ is also set. Since
$\mathbb{x}_{0}\in B(0,\frac{T}{6})$, then $\mathbb{x}^{\prime}\in
B(0,\frac{T}{3})$. For any $\mathbb{x}^{\prime\prime}\in B(0,\frac{T}{6})$,
$\|\mathbb{x}^{\prime}-\mathbb{x}^{\prime\prime}\|<\frac{T}{2}$. There exist
$r_{1}^{*},\cdots,r_{d+1}^{*}$ and
$\mathbb{v}_{1}^{*}\cdots,\mathbb{v}_{d+1}^{*}$ such that

$t\left(r_{1}^{*}\mathbb{v}_{0}+r_{2}^{*}\mathbb{v}_{1}^{*}+\cdots+r_{d+1}^{*}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}r_{k}^{*})\mathbb{v}_{d+1}^{*}\right)=\mathbb{x}^{\prime\prime}-\mathbb{x}^{\prime},$
$r_{k}^{*}\in[0,1],\quad\mathbb{v}_{k}^{*}\in\mathcal{V},\text{ for
}k=1,\cdots,d+1,\text{ and }\sum_{k=1}^{d+1}r_{k}^{*}\leq 1$

Let
$\mathbb{R}=(R_{1},\cdots,R_{d+1},R_{d+2})\sim\text{Dirichlet}(1,1,\cdots,1)$
and $\mathbb{V}=(\mathbb{v}_{1},\cdots,\mathbb{v}_{d+1})\sim\varphi^{d+1}$ be
independent, then for a small enough $\delta>0$,

$\displaystyle\int_{\mathcal{V}}\varphi(d\mathbb{v}_{1})\cdots\int_{\mathcal{V}}\varphi(d\mathbb{v}_{d+1})\int_{r_{1}=0}^{1}dr_{1}\int_{r_{2}=0}^{1-r_{1}}dr_{2}\cdots\int_{r_{d+1}=0}^{1-\sum_{k=1}^{d}r_{k}}dr_{d+1}$
$\displaystyle\times\mathbb{I}_{B(\mathbb{x}^{\prime\prime},\delta)}\left(\mathbb{x}^{\prime}+t\left(r_{1}\mathbb{v}_{0}+r_{2}\mathbb{v}_{1}+\cdots+r_{d+1}\mathbb{v}_{d}+(1-\sum_{k=1}^{d+1}r_{k})\mathbb{v}_{d+1}\right)\right)$
$\displaystyle=\Gamma(d+2)\mathbb{P}\left\\{\Big{|}\mathbb{x}^{\prime}+t\left(R_{1}\mathbb{v}_{0}+R_{2}\mathbb{V}_{1}+\cdots+R_{d+1}\mathbb{V}_{d}+(1-\sum_{k=1}^{d+1}R_{k})\mathbb{V}_{d+1}\right)-\mathbb{x}^{\prime\prime}\Big{|}\leq\delta\right\\}$
$\displaystyle=\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}(R_{1}\mathbb{v}_{0}+R_{2}\mathbb{V}_{1}+\cdots+R_{d+1}\mathbb{V}_{d}+(1-\sum_{k=1}^{d+1}R_{k})\mathbb{V}_{d+1})$
$\displaystyle\qquad-(r_{1}^{*}\mathbb{v}_{0}+r_{2}^{*}\mathbb{v}_{1}^{*}+\cdots+r_{d+1}^{*}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}r_{k}^{*})\mathbb{v}_{d+1}^{*})\Big{|}\leq\frac{\delta}{t}\Big{\\}}$
$\displaystyle\geq\Gamma(d+2)\mathbb{P}\Big{\\{}\\{\Big{|}(R_{1}\mathbb{v}_{0}+R_{2}\mathbb{V}_{1}+\cdots+R_{d+1}\mathbb{V}_{d}+(1-\sum_{k=1}^{d+1}R_{k})\mathbb{V}_{d+1})$
$\displaystyle\qquad-(r_{1}^{*}\mathbb{v}_{0}+r_{2}^{*}\mathbb{v}_{1}^{*}+\cdots+r_{d+1}^{*}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}r_{k}^{*})\mathbb{v}_{d+1}^{*})\Big{|}\leq\frac{\delta}{t}\\}$
$\displaystyle\qquad\cap\left\\{\mathbb{V_{k}}=\mathbb{v}^{*}_{k},k=1,\cdots,d+1\right\\}\Big{\\}}$
$\displaystyle=\Gamma(d+2)\mathbb{P}\Big{\\{}\\{\Big{|}(R_{1}\mathbb{v}_{0}+R_{2}\mathbb{v}_{1}^{*}+\cdots+R_{d+1}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}R_{k})\mathbb{v}_{d+1}^{*})$
$\displaystyle\qquad-(r_{1}^{*}\mathbb{v}_{0}+r_{2}^{*}\mathbb{v}_{1}^{*}+\cdots+r_{d+1}^{*}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}r_{k}^{*})\mathbb{v}_{d+1}^{*})\Big{|}\leq\frac{\delta}{t}\\}$
$\displaystyle\qquad\cap\left\\{\mathbb{V_{k}}=\mathbb{v}^{*}_{k},k=1,\cdots,d+1\right\\}\Big{\\}}$
$\displaystyle=\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}(R_{1}\mathbb{v}_{0}+R_{2}\mathbb{v}_{1}^{*}+\cdots+R_{d+1}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}R_{k})\mathbb{v}_{d+1}^{*})$
$\displaystyle\qquad-(r_{1}^{*}\mathbb{v}_{0}+r_{2}^{*}\mathbb{v}_{1}^{*}+\cdots+r_{d+1}^{*}\mathbb{v}_{d}^{*}+(1-\sum_{k=1}^{d+1}r_{k}^{*})\mathbb{v}_{d+1}^{*})\Big{|}\leq\frac{\delta}{t}\Big{\\}}$
$\displaystyle\qquad\times\mathbb{P}\Big{\\{}\left\\{\mathbb{V_{k}}=\mathbb{v}^{*}_{k},k=1,\cdots,d+1\right\\}\Big{\\}}$
$\displaystyle\geq\left(\frac{\lambda^{\text{ref}}}{2d\lambda^{\text{ref}}+K}\right)^{d+1}\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}(R_{1}-r_{1}^{*})\mathbb{v}_{0}+\sum_{k=2}^{d+1}(R_{k}-r_{k}^{*})\mathbb{v}_{k-1}^{*}-\sum_{k=1}^{d+1}(R_{k}-r_{k}^{*})\mathbb{v}^{*}_{d+1}\Big{|}\leq\frac{\delta}{t}\Big{\\}}$
$\displaystyle\geq\left(\frac{\lambda^{\text{ref}}}{2d\lambda^{\text{ref}}+K}\right)^{d+1}\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}(R_{1}-r_{1}^{*})\mathbb{v}_{0}\Big{|}+\sum_{k=2}^{d+1}\Big{|}(R_{k}-r_{k}^{*})\mathbb{v}_{k-1}^{*}\Big{|}+\sum_{k=1}^{d+1}\Big{|}(R_{k}-r_{k}^{*})\mathbb{v}^{*}_{d+1}\Big{|}\leq\frac{\delta}{t}\Big{\\}}$
$\displaystyle\geq\left(\frac{\lambda^{\text{ref}}}{2d\lambda^{\text{ref}}+K}\right)^{d+1}\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}(R_{1}-r_{1}^{*})\mathbb{v}_{0}\Big{|}\leq\frac{\delta}{2(d+1)t},$
$\displaystyle\qquad\Big{|}(R_{k}-r_{k}^{*})\mathbb{v}_{k-1}^{*}\Big{|}\leq\frac{\delta}{2(d+1)t},k=2,\cdots,d+1,$
$\displaystyle\qquad\Big{|}(R_{k}-r_{k}^{*})\mathbb{v}^{*}_{d+1}\Big{|}\leq\frac{\delta}{2(d+1)t},k=1,\cdots,d+1\Big{\\}}$
$\displaystyle=\left(\frac{\lambda^{\text{ref}}}{2d\lambda^{\text{ref}}+K}\right)^{d+1}\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}R_{k}-r_{k}^{*}\Big{|}\leq\frac{\delta}{2(d+1)t},k=1,\cdots,d+1\Big{\\}}$
$\displaystyle\geq\left(\frac{\lambda^{\text{ref}}}{2d\lambda^{\text{ref}}+K}\right)^{d+1}\Gamma(d+2)\mathbb{P}\Big{\\{}\Big{|}R_{k}-r_{k}^{*}\Big{|}\leq\frac{\delta}{2(d+1)T},k=1,\cdots,d+1\Big{\\}}\geq
C_{1}$

If we set an arbitrary value $\delta=\delta(d,T)>0$, such that $\delta$ is
small enough, we define $h:\mathbb{R}^{d+1}\rightarrow\mathbb{R}$ as

$h(r_{1},\cdots,r_{d+1})=\mathbb{P}\Big{\\{}\Big{|}R_{k}-r_{k}\Big{|}\leq\frac{\delta}{2(d+1)T},k=1,\cdots,d+1\Big{\\}}$

where $\mathbf{R}=(R_{1},\cdots,R_{d+2})\sim\text{Dirichlet}(1,1,\cdots,1)$.
Then $h(r_{1},\cdots,r_{d+1})>0$ for any vector $(r_{1},\cdots,r_{d+1})\in F$,
where $F=\\{(r_{1},\cdots,r_{d+1})\in[0,1]^{d+1}:\sum_{k=1}^{d+1}r_{k}\leq
1\\}$. Since $F$ is compact, as a result, there exists a constant
$\eta_{0}>0$, such that $\min_{\mathbb{r}\in F}h(\mathbb{r})>\eta_{0}$. Here
$\eta_{0}$ is fixed and only depends on $T,d$. Hence, $C_{1}$ only depends on
$d,T,\lambda^{\text{ref}}$

As a result, we have, for any $t>\frac{5T}{6}$,

$\begin{split}&\int_{\mathcal{V}}\varphi(d\mathbb{v}_{1})\cdots\int_{\mathcal{V}}\varphi(d\mathbb{v}_{d+1})\int_{r_{1}=0}^{1}dr_{1}\int_{r_{2}=0}^{1-r_{1}}dr_{2}\cdots\int_{r_{d+1}=0}^{1-\sum_{k=1}^{d}r_{k}}dr_{d+1}\\\
&\times
f(\mathbb{x}_{0}+t\sum_{k=1}^{d+1}r_{k}\mathbb{v}_{k-1}+t(1-\sum_{k=1}^{d+1}r_{k})\mathbb{v}_{d+1}+(T-t)\mathbb{v}_{d+2})\\\
&\geq
C_{2}(T,d,\lambda^{\text{ref}})\int_{B(0,\frac{T}{6})}f(\mathbb{x}^{\prime\prime},\mathbb{v}_{d+2})d\mathbb{x}^{\prime\prime}\end{split}$

and

$\begin{split}&\mathbb{E}^{\mathbb{z}_{0}}[f(\mathbb{Z}_{T})]\\\
&\geq\left(\frac{(\lambda^{\text{ref}})^{2}}{\lambda^{\text{ref}}+K}\right)^{d+2}\exp\left\\{-T\left(\lambda^{\text{ref}}+K\right)\right\\}\int_{\mathcal{V}}\varphi(d\mathbb{v}_{d+2})\times\int_{t=\frac{5T}{6}}^{T}dt\int_{B(0,\frac{T}{6})}f(\mathbb{x}^{\prime\prime},\mathbb{v}_{d+2})d\mathbb{x}^{\prime\prime}\\\
&\geq
C_{3}(T,d,\lambda^{\text{ref}})\int_{\mathcal{V}}\varphi(d\mathbb{v}_{d+2})\int_{t=\frac{5T}{6}}^{T}dt\int_{B(0,\frac{T}{6})}f(\mathbb{x}^{\prime\prime},\mathbb{v}_{d+2})d\mathbb{x}^{\prime\prime}\\\
&\geq
C_{4}(T,d,\lambda^{\text{ref}})\int_{\mathcal{V}}\varphi(d\mathbb{v}_{d+2})\int_{B(0,\frac{T}{6})}f(\mathbb{x}^{\prime\prime},\mathbb{v}_{d+2})d\mathbb{x}^{\prime\prime}\\\
\end{split}$

Hence, for any Borel set $A\subset\mathbb{R}^{d}\times\mathcal{V}$ and
$\mathbb{z}_{0}=(\mathbb{x}_{0},\mathbb{v}_{0})\in\mathbb{R}^{d}\times\mathcal{V}$,
setting $f=\mathbb{I}_{A}$ and using above arguments, we have

$\mathbb{P}(\mathbb{z}_{0},\mathbb{Z}_{T}\in A)\geq
C_{4}(T,d,\lambda^{\text{ref}})\int\int_{A}\varphi(d\mathbb{v})d\mathbb{x}$

Consequently, for any $R>0$, the set $B(0,R)\times\mathcal{V}$ is petite.
Hence, any compact set is petite and irreducibility follows.

Proof of Theorem 3: Using the same arguments as in Lemma 3 of Deligiannidis et
al., (2017) and the above Lemma 2, the Markov process induced by our
Coordinate Sampler is ergodic, hence its invariant distribution is unique. By
Lemma 1, $\rho(d\mathbb{z})$ is the unique invariant distribution of the
coordinate sampler.

##  7 Proof of Theorem 4

In this section, we use the techniques developed in Deligiannidis et al.,
(2017). We will again detail the proof that $V$ is the desired Lyapunov
function, since there are some differences between our proof and the original
one.

###### Lemma 3

(Down et al.,, 1995, Theorem 5.2) Let $\left\\{\mathbb{Z}_{t}:t\geq
0\right\\}$ be a Borel right Markov process taking values in a locally
compact, separable metric space $\mathcal{Z}$ and assume it is non-explosive,
irreducible and aperiodic. Let $(\mathcal{L},\mathcal{D}(\mathcal{L}))$ be its
extended generator. Suppose that there exists a measurable function
$V:\mathcal{Z}\rightarrow[1,\infty)$ such that $V\in\mathcal{D}(\mathcal{L})$,
and that for a petite set $C\in\mathcal{B}(\mathcal{Z})$ and constants
$b,c>0$, we have

$\mathcal{L}V\leq-cV+b\mathbb{I}_{C},$

Then $\\{\mathbb{Z}_{t}:t\geq 0\\}$ is V-uniformly ergodic.

Proof of Theorem 4: In Section 5.1. of Deligiannidis et al., (2017), $V$
defined in the paper belongs to the extended generator
$\mathcal{D}(\mathcal{L})$, given Assumptions $A.1-A.4$. We next show that $V$
is a Lyapunov function.

Case 1: $\langle\nabla U(\mathbb{x}),\mathbb{v}\rangle>0$.
$V(\mathbb{x},\mathbb{v})=\frac{e^{U(\mathbb{x})/2}}{\sqrt{\lambda^{\text{ref}}}}$
and
$\nabla_{\mathbb{x}}V(\mathbb{x},\mathbb{v})=\frac{1}{2}V(\mathbb{x},\mathbb{v})\nabla
U(\mathbb{x})$.

$\displaystyle\mathcal{L}V(\mathbb{x},\mathbb{v})$
$\displaystyle=\frac{1}{2}V(\mathbb{x},\mathbb{v})\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle+\left(\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle_{+}\right)\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}_{i}\rangle_{+}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}$ $\displaystyle\times
V(\mathbb{x},\mathbb{v})\sqrt{\frac{\lambda^{\text{ref}}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}_{i}\rangle}_{+}}-\left(\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle_{+}\right)V(\mathbb{x},\mathbb{v})$
$\displaystyle=V(\mathbb{x},\mathbb{v})\left\\{-\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\left(\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle\right)\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x}\right)\big{|})}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}$
$\displaystyle=V(\mathbb{x},\mathbb{v})\Bigg{\\{}\lambda^{\text{ref}}\left[\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x}\right)\big{|})}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}-1\right]$ $\displaystyle+\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle\left[\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x}\right)\big{|})}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}-\frac{1}{2}\right]\Bigg{\\}}$

Case 2: $\langle\nabla U(\mathbb{x}),\mathbb{v}\rangle<0$.

$V(\mathbb{x},\mathbb{v})=\frac{e^{U(\mathbb{x})/2}}{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}}$ $\frac{\partial
V(\mathbb{x},\mathbb{v})}{\partial\mathbb{x}}=\frac{1}{2}V(\mathbb{x},\mathbb{v})\nabla
U(\mathbb{x})+\frac{1}{2}V(\mathbb{x},\mathbb{v})\frac{\Delta
U(\mathbb{x})\mathbb{v}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}$
$\displaystyle\mathcal{L}V(\mathbb{x},\mathbb{v})$
$\displaystyle=\frac{1}{2}V(\mathbb{x},\mathbb{v})\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle+\frac{1}{2}V(\mathbb{x},\mathbb{v})\frac{\langle\mathbb{v},\Delta
U(\mathbb{x})\mathbb{v}\rangle}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}-\lambda^{\text{ref}}V(\mathbb{x},\mathbb{v})$
$\displaystyle+\lambda^{\text{ref}}\sum_{i=1}^{2d}\frac{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}_{i}\rangle_{+}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}V(\mathbb{x},\mathbb{v})\sqrt{\frac{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}_{i}\rangle}_{+}}$
$\displaystyle=V(\mathbb{x},\mathbb{v})\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{1}{2}\frac{\langle\mathbb{v},\Delta
U(\mathbb{x})\mathbb{v}\rangle}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+V(\mathbb{x},\mathbb{v})\left\\{\lambda^{\text{ref}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\sum_{i=1}^{d}\frac{\sqrt{\lambda^{\text{ref}}}+\sqrt{\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x})\big{|}}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}$

Case 3: $\langle U(\mathbb{x}),\mathbb{v}\rangle=0$. The generator is defined
as

$\mathcal{L}V(\mathbb{x},\mathbb{v})=\frac{dV(\mathbb{x}+t\mathbb{v},\mathbb{v})}{dt}\Big{|}_{t=0+}+\lambda(\mathbb{x},\mathbb{v})\sum_{\mathbb{v}^{\prime}\in\mathcal{V}}\frac{\lambda(\mathbb{x},-\mathbb{v}^{\prime})}{\lambda(\mathbb{x})}\left(f(\mathbb{x},\mathbb{v}^{\prime})-f(\mathbb{x},\mathbb{v})\right)$

$(i)$: If $\langle\mathbb{v},\Delta U(\mathbb{x})\mathbb{v}\rangle>0$, then

$\begin{split}&\frac{dV(\mathbb{x}+t\mathbb{v},\mathbb{v})}{dt}\Big{|}_{t=0+}=\lim_{t\to
0+}\frac{1}{t}\left\\{\frac{e^{U(\mathbb{x}+t\mathbb{v})/2}}{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}+t\mathbb{v}),-\mathbb{v}\rangle_{+}}}-\frac{e^{U(\mathbb{x})/2}}{\sqrt{\lambda^{\text{ref}}}}\right\\}\\\
&=\frac{1}{2}\frac{e^{U(\mathbb{x}+t\mathbb{v})/2}}{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}+t\mathbb{v}),-\mathbb{v}\rangle_{+}}}\langle\nabla
U(\mathbb{x}+t\mathbb{v}),\mathbb{v}\rangle\Big{|}_{t=0+}=0\end{split}$

$(ii)$: $\langle\mathbb{v},\Delta U(\mathbb{x})\mathbb{v}\rangle<0$, then

$\begin{split}&\frac{dV(\mathbb{x}+t\mathbb{v},\mathbb{v})}{dt}\Big{|}_{t=0+}\\\
&=\lim_{t\to
0+}\frac{1}{t}\left\\{\frac{e^{U(\mathbb{x}+t\mathbb{v})/2}}{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}+t\mathbb{v}),-\mathbb{v}\rangle_{+}}}-\frac{e^{U(\mathbb{x})/2}}{\sqrt{\lambda^{\text{ref}}}}\right\\}\\\
&=\frac{1}{2}\frac{e^{U(\mathbb{x}+t\mathbb{v})/2}}{\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}+t\mathbb{v}),-\mathbb{v}\rangle_{+}}}\langle\nabla
U(\mathbb{x}+t\mathbb{v}),\mathbb{v}\rangle\Big{|}_{t=0+}-\frac{1}{2}\frac{e^{U(\mathbb{x}+t\mathbb{v})/2}\langle\mathbb{v},-\Delta
U(\mathbb{x})\mathbb{v}\rangle}{\left(\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}+t\mathbb{v}),-\mathbb{v}\rangle_{+}\right)^{3/2}}\Big{|}_{t=0+}\\\
&=0+\frac{1}{2}\frac{e^{U(\mathbb{x})/2}\langle\mathbb{v},\Delta
U(\mathbb{x})\mathbb{v}\rangle}{(\sqrt{\lambda^{\text{ref}}})^{3}}\end{split}$

As a result,

$\begin{split}\quad\frac{dV(\mathbb{x}+t\mathbb{v},\mathbb{v})}{dt}\Big{|}_{t=0+}=-\frac{1}{2}\frac{e^{U(\mathbb{x})/2}\langle\mathbb{v},-\Delta
U(\mathbb{x})\mathbb{v}\rangle_{+}}{(\sqrt{\lambda^{\text{ref}}})^{3}}\end{split}$
$\begin{split}&\mathcal{L}V(\mathbb{x},\mathbb{v})\\\
&=V(\mathbb{x},\mathbb{v})\left\\{-\frac{1}{2}\frac{\langle\mathbb{v},-\Delta
U(\mathbb{x})\mathbb{v}\rangle_{+}}{\lambda^{\text{ref}}}+\lambda^{\text{ref}}\left[\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x}\right)\big{|})}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}-1\right]\right\\}\end{split}$

Condition 1: $\displaystyle{\overline{\lim}_{|x|\to\infty}}\|\Delta
U(x)\|\leq\alpha_{1}<\infty$,
$\displaystyle{\underline{\lim}_{|x|\to\infty}}|\nabla U(x)|_{1}=\infty$ and
$\lambda^{\text{ref}}>\sqrt{8\alpha_{1}}$.

$\begin{split}&\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}<\frac{1}{2}\\\
&\iff\frac{d\lambda^{\text{ref}}+d\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}<\frac{1}{2}\\\ &\iff
2\sqrt{d}\sqrt{d(\lambda^{\text{ref}})^{2}+\lambda^{\text{ref}}\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}<\big{|}\nabla U(\mathbb{x})\big{|}_{1}\\\ &\iff
4d^{2}(\lambda^{\text{ref}})^{2}+4d\lambda^{\text{ref}}\big{|}\nabla
U(\mathbb{x})\big{|}_{1}<\big{|}\nabla U(\mathbb{x})\big{|}_{1}^{2}\\\
&\iff\big{|}\nabla
U(\mathbb{x})\big{|}_{1}>2(\sqrt{2}+1)d\lambda^{\text{ref}}\end{split}$

Denote $K=K_{1}\cup K_{2}\cup K_{3}$, where $K_{1}=\\{\mathbb{x}:\big{|}\nabla
U(\mathbb{x})\big{|}_{1}\leq 2(\sqrt{2}+1)d\lambda^{\text{ref}}\\}$,
$K_{1}=\\{\mathbb{x}:\big{|}\nabla
U(\mathbb{x})\big{|}_{1}<16d\lambda^{\text{ref}}\\}$ and
$K_{3}=\\{\mathbb{x}:\|\Delta U(\mathbb{x})\|\leq 2\alpha_{1}\\}$. On $K^{c}$,
we have

$\begin{split}&\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x})\big{|}\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\\\
&\leq\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\text{ (Jensen's inequality)}\\\
&<\frac{1}{2}\end{split}$

As a result, for case $1$ and case $3$, we have

$\frac{\mathcal{L}V(\mathbb{x},\mathbb{v})}{V(\mathbb{x},\mathbb{v})}\leq-\frac{1}{2}\lambda^{\text{ref}},\quad\text{On
$K^{c}\times\mathcal{V}$}$

For case $2$, we have

$\displaystyle\frac{\mathcal{L}V(\mathbb{x},\mathbb{v})}{V(\mathbb{x},\mathbb{v})}$
$\displaystyle=\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{1}{2}\frac{\langle\mathbb{v},\Delta
U(\mathbb{x})\mathbb{v}\rangle}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\left\\{\lambda^{\text{ref}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\sum_{i=1}^{d}\frac{\sqrt{\lambda^{\text{ref}}}+\sqrt{\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x})\big{|}}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\left\\{\lambda^{\text{ref}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\frac{d\sqrt{\lambda^{\text{ref}}}+d\sqrt{\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}\text{ (Jensen's Inequality)}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\times\left\\{\frac{\sqrt{\lambda^{\text{ref}}}}{2}+\frac{\sqrt{d\lambda^{\text{ref}}}}{\sqrt{2d+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/\lambda^{\text{ref}}}}\right\\}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\times\left\\{\frac{\sqrt{\lambda^{\text{ref}}}}{2}+\frac{\sqrt{\lambda^{\text{ref}}}}{4}\right\\}\text{(
since $\big{|}\nabla U(\mathbb{x})\big{|}_{1}>16d\lambda^{\text{ref}}$)}$
$\displaystyle=\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}+\frac{3}{4}\sqrt{\lambda^{\text{ref}}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}$

Denote $w_{0}=\langle\nabla U(\mathbb{x}),-\mathbb{v}\rangle$ and define
$g:\mathbb{R}^{+}\rightarrow\mathbb{R}$ as

$g(w)=-\frac{1}{2}w-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}+w}+\frac{3}{4}\sqrt{\lambda^{\text{ref}}}\sqrt{\lambda^{\text{ref}}+w}$
$g^{\prime}(w)=-\frac{1}{2}-\frac{\alpha_{1}}{\lambda^{(\text{ref}}+w)^{2}}+\frac{3}{8}\frac{\sqrt{\lambda^{\text{ref}}}}{\sqrt{\lambda^{\text{ref}}+w}}<-\frac{1}{8}<0$

As a result, $g$ is a decreasing function on $[0,\infty)$, and

$g(w_{0})\leq
g(0)=-\lambda^{\text{ref}}+\frac{\alpha_{1}}{\lambda^{\text{ref}}}+\frac{3}{4}\lambda^{\text{ref}}\leq-\frac{1}{8}\lambda^{\text{ref}},\quad\text{(since
$\lambda^{\text{ref}}>\sqrt{8\alpha_{1}}$)}$

As a result, on $K^{c}\times\mathcal{V}$, we have

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq-\frac{1}{8}\lambda^{\text{ref}}V(\mathbb{x},\mathbb{v})$

Since $K$ is compact and $\mathcal{V}$ is finite, hence, $K\times\mathcal{V}$
is compact. As a result, there exists $b>0$ such that

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq b,\quad\text{ for all
$(\mathbb{x},\mathbb{v})\in K\times\mathcal{V}$}$

Hence, under condition 1, there exist constants $b>0$, $c>0$ and a function
$V$ such that

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq-
cV(\mathbb{x},\mathbb{v})+b\mathbb{I}_{K\times\mathcal{V}}$

Condition 2: $\underline{\lim}_{|\mathbb{x}|\rightarrow\infty}\big{|}\nabla
U(\mathbb{x})\big{|}_{1}=2\alpha_{2}>0$,
$\displaystyle{\overline{\lim}_{|x|\to\infty}}\|\Delta U(x)\|=0$ and
$\lambda^{\text{ref}}<\frac{\alpha_{2}}{14d}$. Denote
$K_{1}=\\{\mathbb{x}:\big{|}\nabla U(\mathbb{x})\big{|}_{1}\leq\alpha_{2}\\}$,
by the same arguments as above, we have, on $(\mathbb{x},\mathbb{v})\in
K_{1}^{c}\times\mathcal{V}$,

$\begin{split}&\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x})\big{|}\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\\\
&\leq\sum_{i=1}^{d}\frac{\lambda^{\text{ref}}+\sqrt{\lambda^{\text{ref}}\left(\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d\right)}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\text{ (Jensen's inequality)}\\\ &<\frac{1}{2}\text{
\quad(since $\big{|}\nabla
U(\mathbb{x})\big{|}_{1}>\alpha_{2}>14d\lambda^{\text{ref}}>2(\sqrt{2}+1)d\lambda^{\text{ref}}$)}\end{split}$

In cases $1$ and $3$, we have, for all $(\mathbb{x},\mathbb{v})\in
K_{1}^{c}\times\mathcal{V}$

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq-\frac{1}{2}V(\mathbb{x},\mathbb{v})$

Since ${\overline{\lim}_{|x|\to\infty}}\|\Delta U(x)\|=0$, there exist $M>0$
and $\epsilon<(\lambda^{\text{ref}})^{2}/8$, such that $\|\Delta
U(x)\|<\epsilon$, for all $x\in
K_{2}^{c}=\\{\mathbb{x}:\big{|}\mathbb{x}\big{|}>M\\}$. Define $K=K_{1}\cup
K_{2}$, then, in case $2$, for all $(\mathbb{x},\mathbb{v})\in
K^{c}\times\mathcal{V}$

$\displaystyle\frac{\mathcal{L}V(\mathbb{x},\mathbb{v})}{V(\mathbb{x},\mathbb{v})}$
$\displaystyle=\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{1}{2}\frac{\langle\mathbb{v},\Delta
U(\mathbb{x})\mathbb{v}\rangle}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\left\\{\lambda^{\text{ref}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\sum_{i=1}^{d}\frac{\sqrt{\lambda^{\text{ref}}}+\sqrt{\lambda^{\text{ref}}+\big{|}\nabla_{i}U(\mathbb{x})\big{|}}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\left\\{\lambda^{\text{ref}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\frac{d\sqrt{\lambda^{\text{ref}}}+d\sqrt{\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}/d}}{2d\lambda^{\text{ref}}+\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}\right\\}\text{ (Jensen's Inequality)}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\times\left\\{\frac{\sqrt{\lambda^{\text{ref}}}}{2}+\frac{\sqrt{\lambda^{\text{ref}}}}{\sqrt{2+\frac{\big{|}\nabla
U(\mathbb{x})\big{|}_{1}}{d\lambda^{\text{ref}}}}}\right\\}$
$\displaystyle\leq\left\\{\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\right\\}$
$\displaystyle+\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}\times\left\\{\frac{\sqrt{\lambda^{\text{ref}}}}{2}+\frac{\sqrt{\lambda^{\text{ref}}}}{4}\right\\}\text{(
since $\big{|}\nabla U(\mathbb{x})\big{|}_{1}>\alpha_{2}$ and
$\frac{\alpha_{2}}{d\lambda^{\text{ref}}}>14$)}$
$\displaystyle=\frac{1}{2}\langle\nabla
U(\mathbb{x}),\mathbb{v}\rangle-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}+\frac{3}{4}\sqrt{\lambda^{\text{ref}}}\sqrt{\lambda^{\text{ref}}+\langle\nabla
U(\mathbb{x}),-\mathbb{v}\rangle}$

Denote $w_{0}=\langle\nabla U(\mathbb{x}),-\mathbb{v}\rangle$, and define
$g:\mathbb{R}^{+}\rightarrow\mathbb{R}$ as

$g(w)=-\frac{1}{2}w-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}+w}+\frac{3}{4}\sqrt{\lambda^{\text{ref}}}\sqrt{\lambda^{\text{ref}}+w}$
$g^{\prime}(w)=-\frac{1}{2}-\frac{\epsilon}{\lambda^{(\text{ref}}+w)^{2}}+\frac{3}{8}\frac{\sqrt{\lambda^{\text{ref}}}}{\sqrt{\lambda^{\text{ref}}+w}}<-\frac{1}{8}<0$

As a result, $g$ is a decreasing function on $[0,\infty)$, and

$g(w_{0})\leq
g(0)=-\lambda^{\text{ref}}+\frac{\epsilon}{\lambda^{\text{ref}}}+\frac{3}{4}\lambda^{\text{ref}}\leq-\frac{1}{8}\lambda^{\text{ref}},\quad\text{(since
$\lambda^{\text{ref}}>\sqrt{8\epsilon}$)}$

As a result, on $K^{c}\times\mathcal{V}$, for all three cases, we have

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq-\frac{1}{8}\lambda^{\text{ref}}V(\mathbb{x},\mathbb{v})$

Since $K_{1}\subset K$ and $K$ is compact, $\mathcal{V}$ is finite, therefore,
$K\times\mathcal{V}$ is compact. As a result, there exists $b>0$ such that

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq b,\quad\text{ for all
$(\mathbb{x},\mathbb{v})\in K\times\mathcal{V}$}$

Hence, under condition 2, there exist constants $b>0$, $c>0$ and a function
$V$ such that

$\mathcal{L}V(\mathbb{x},\mathbb{v})\leq-
cV(\mathbb{x},\mathbb{v})+b\mathbb{I}_{K\times\mathcal{V}}$

Since each compact set is a petite set, $V$ is hence a Lyapunov function. As a
result, by Lemma 3, Theorem 4 in the paper is proved.

##  8 Event rates for all experiments

In this section, we produce the form of the event rate
$\lambda(\mathbb{x},\mathbb{v})$ of each model for both Zigzag and Coordinate
samplers.

###  8.1 Banana-shaped Distribution

In this example, the potential function is

$U(\mathbb{x})=U(x_{1},x_{2})=(x_{1}-1)^{2}+\kappa(x_{2}-x_{1}^{2})^{2},$

and its gradient is as follows,

$\frac{\partial U(\mathbb{x})}{\partial
x_{1}}=2(x_{1}-1)+4\kappa(x_{1}^{2}-x_{2})x_{1},\quad\frac{\partial
U(\mathbb{x})}{\partial x_{2}}=2\kappa(x_{2}-x_{1}^{2})$

For the Zigzag Sampler, recall that we set $\lambda_{i}^{\text{ref}}=0$,
$i=1,2$,

$\displaystyle\lambda_{1}(x+tv,v)=\left\\{v_{1}\frac{\partial
U(x+tv)}{\partial x_{1}}\right\\}_{+}$
$\displaystyle=\left\\{v_{1}\left\\{2(x_{1}+tv_{1}-1)+4\kappa(x_{1}+tv_{1})((x_{1}+tv_{1})^{2}-(x_{2}+tv_{2}))\right\\}\right\\}_{+}$
$\displaystyle=\left\\{2x_{1}v_{1}+2tv_{1}^{2}-2v_{1}+4\kappa
v_{1}(x_{1}+tv_{1})^{3}-4\kappa
v_{1}(x_{1}+tv_{1})(x_{2}+tv_{2})\right\\}_{+}$
$\displaystyle=\left\\{2x_{1}v_{1}+2tv_{1}^{2}-2v_{1}+4\kappa
v_{1}(x_{1}^{3}+3x_{1}^{2}v_{1}t+3x_{1}v_{1}^{2}t^{2}+v_{1}^{3}t^{3})-4\kappa
v_{1}(x_{1}x_{2}+(x_{1}v_{2}+x_{2}v_{1})t+v_{1}v_{2}t^{2})\right\\}_{+}$
$\displaystyle=\Big{\\{}(4\kappa v_{1}^{4})t^{3}+(12\kappa
x_{1}v_{1}^{3}-4\kappa v_{1}^{2}v_{2})t^{2}+(2v_{1}^{2}+12\kappa
x_{1}^{2}v_{1}^{2}-4\kappa x_{1}v_{1}v_{2}-4\kappa x_{2}v_{1}^{2})t$
$\displaystyle+(2x_{1}v_{1}-2v_{1}+4\kappa v_{1}x_{1}^{3}-4\kappa
v_{1}x_{1}x_{2})\Big{\\}}_{+}$
$\displaystyle=\left\\{a_{1,3}t^{3}+a_{1,2}t^{2}+a_{1,1}t+a_{1,0}\right\\}_{+}$

where $a_{1,3}=4\kappa v_{1}^{4}$, $a_{1,2}=12\kappa x_{1}v_{1}^{3}-4\kappa
v_{1}^{2}v_{2}$, $a_{1,1}=2v_{1}^{2}+12\kappa x_{1}^{2}v_{1}^{2}-4\kappa
x_{1}v_{1}v_{2}-4\kappa x_{2}v_{1}^{2}$, and
$a_{1,0}=2x_{1}v_{1}-2v_{1}+4\kappa v_{1}x_{1}^{3}-4\kappa v_{1}x_{1}x_{2}$.

$\begin{split}&\lambda_{2}(x+tv,v)=\left\\{v_{2}\frac{\partial
U(x+tv)}{\partial x_{2}}\right\\}_{+}\\\
&=\left\\{v_{2}2\kappa\left(x_{2}+tv_{2}-(x_{1}+tv_{1})^{2}\right)\right\\}_{+}\\\
&=\left\\{(-2\kappa v_{1}^{2}v_{2})t^{2}+\left(2\kappa
v_{2}(v_{2}-2x_{1}v_{1})\right)t+\left(2\kappa
v_{2}(x_{2}-x_{1}^{2})\right)\right\\}_{+}\\\
&=\left\\{a_{2,2}t^{2}+a_{2,1}t+a_{2,0}\right\\}_{+}\end{split}$

where $a_{2,2}=-2\kappa v_{1}^{2}v_{2}$, $a_{2,1}=2\kappa
v_{2}(v_{2}-2x_{1}v_{1})$ and $a_{2,0}=2\kappa v_{2}(x_{2}-x_{1}^{2})$. As a
result, we have the following upper bounds for $\lambda_{1}$ and
$\lambda_{2}$.

$\lambda_{1}(\mathbb{x}+t\mathbb{v},\mathbb{v})\leq\bar{\lambda}_{1}(t):=\left(\sum_{i=0}^{3}\max\\{0,a_{1,i}\\}\right)\max\\{1,t^{3}\\}$
$\lambda_{2}(\mathbb{x}+t\mathbb{v},\mathbb{v})\leq\bar{\lambda}_{2}(t):=\left(\sum_{i=0}^{2}\max\\{0,a_{2,i}\\}\right)\max\\{1,t^{2}\\}$

First, we use the Superposition Theorem: we set $T_{1}=0$ and generate a time
duration, $\tau$, from the Poisson process with rate $\bar{\lambda}_{1}(t)$,
then compute
$p=\lambda_{1}(\mathbb{x}+t\mathbb{v},\mathbb{v})/\bar{\lambda}_{1}(t)$ and
accept $\tau$ with probability $p$. If it is rejected, we update
$\mathbb{x}\leftarrow\mathbb{x}+\tau\mathbb{v},\mathbb{v}\leftarrow\mathbb{v},T_{1}\leftarrow
T_{1}+\tau$ and repeat the above process, until we obtain one $\tau$ and set
$T_{1}=T_{1}+\tau$. Apply this procedure on $\lambda_{2}$ and get $T_{2}$. By
the Thinning Theorem, $\min\\{T_{1},T_{2}\\}$ follows the Poisson process with
rate $\lambda_{1}+\lambda_{2}$.

For the Coordinate Sampler, if $\mathbb{v}=(v_{1},0)$, where
$v_{1}\in\\{-1,1\\}$, then

$\begin{split}&\lambda(x+tv,v)=\left\\{v_{1}\frac{\partial U(x+tv)}{\partial
x_{1}}\right\\}_{+}\\\
&=\left\\{v_{1}\left\\{2(x_{1}+tv_{1}-1)+4\kappa(x_{1}+tv_{1})((x_{1}+tv_{1})^{2}-x_{2}))\right\\}\right\\}_{+}\\\
&=\left\\{2x_{1}v_{1}+2tv_{1}^{2}-2v_{1}+4\kappa
v_{1}(x_{1}+tv_{1})^{3}-4\kappa v_{1}(x_{1}+tv_{1})x_{2}\right\\}_{+}\\\
&=\left\\{(4\kappa v_{1}^{4})t^{3}+(12\kappa
x_{1}v_{1}^{3})t^{2}+(2v_{1}^{2}+12\kappa x_{1}^{2}v_{1}^{2}-4\kappa
x_{2}v_{1}^{2})t+(2x_{1}v_{1}-2v_{1}+4\kappa x_{1}^{3}v_{1}-4\kappa
x_{1}x_{2}v_{1})\right\\}_{+}\\\
&=\left\\{b_{1,3}t^{3}+b_{1,2}t^{2}+b_{1,1}t+b_{1,0}\right\\}_{+}\end{split}$

if $\mathbb{v}=(0,v_{2})$, where $v_{2}\in\\{-1,1\\}$, then

$\begin{split}&\lambda_{1}(x+tv,v)=\left\\{v_{1}\frac{\partial
U(x+tv)}{\partial x_{1}}\right\\}_{+}\\\
&=\left\\{v_{2}\left\\{2\kappa(x_{2}+tv_{2}-x_{1}^{2})\right\\}\right\\}_{+}\\\
&=\left\\{(2\kappa v_{2}^{2})t+\left(2\kappa
v_{2}(x_{2}-x_{1}^{2})\right)\right\\}_{+}\\\
&=\left\\{b_{2,1}t+b_{2,0}\right\\}_{+}\end{split}$

At current event time, if $\mathbb{v}=(v_{1},0)$, we generate the event
duration as above via the Superposition Theorem. If $\mathbb{v}=(0,v_{2})$, we
generate the event duration directly. That is, $U\sim\text{Uniform}[0,1]$, if
$b_{2,0}>0$, then the time duration is
$\left(\sqrt{-2\log(U)b_{2,1}+b_{2,0}^{2}}-b_{2,0}\right)\Big{/}b_{2,1}$.
Otherwise, the time duration is
$\left(\sqrt{-2\log(U)b_{2,1}}-b_{2,0}\right)\Big{/}b_{2,1}$.

###  8.2 Multivariate Gaussian Distribution

In this model,

$U(\mathbb{x})=\frac{1}{2}\mathbb{x}^{T}A^{-1}\mathbb{x},\text{ and }\nabla
U(\mathbb{x})=A^{-1}\mathbb{x}$

For simplicity, we denote $B=A^{-1}$ and $B=(b_{ij})_{i,j=1,\cdots,d}$.

In Zigzag Sampler,

$\lambda_{i}(\mathbb{x}+t\mathbb{v},\mathbb{v})=\left\\{v_{i}\sum_{j=1}^{d}b_{ij}(x_{j}+tv_{j})\right\\}_{+}=\left\\{\left(v_{i}\sum_{j=1}^{d}b_{ij}v_{j}\right)t+\left(v_{i}\sum_{j=1}^{d}b_{ij}x_{j}\right)\right\\}_{+}$

In Coordinate Sampler, if only $v_{i}\neq 0$ for $\mathbb{v}$, then

$\lambda(\mathbb{x}+t\mathbb{v},\mathbb{v})=\left\\{v_{i}\sum_{j=1}^{d}b_{ij}x_{j}+tb_{ii}v_{i}^{2}\right\\}_{+}$

###  8.3 Bayesian Logistic Model

In this example, we sample from the posterior of a Bayesian logistic model
with flat prior and without intercept. Let $\mathbb{r}\in\mathbb{R}^{d}$ be
independent variable, $s\in\\{0,1\\}$ be the response variable. The model is

$\mathbb{P}(s=1|\mathbb{r})=\frac{\exp(\mathbb{r}^{T}\mathbb{x})}{1+\exp(\mathbb{r}^{T}\mathbb{x})}$

where $\mathbb{x}$ denotes the set of parameters. For a sample of observations
$\\{(\mathbb{r}_{n},s_{n})\\}_{n=1}^{N}$, the likelihood function is

$L(\mathbb{x})=\prod_{n=1}^{N}\left(\frac{\exp(\mathbb{r}_{n}^{T}\mathbb{x})}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}\right)^{s_{n}}\left(\frac{1}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}\right)^{1-s_{n}}=\prod_{n=1}^{N}\frac{\exp(s_{n}r_{n}^{T}\mathbb{x})}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}$

The potential function is

$U(\mathbb{x})=\sum_{n=1}^{N}\left\\{\log\left(1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})\right)-s_{n}\mathbb{r}_{n}^{T}\mathbb{x}\right\\}$
$\begin{split}\nabla
U(\mathbb{x})&=\sum_{n=1}^{N}\left\\{\frac{\exp(\mathbb{r}_{n}^{T}\mathbb{x})}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}-s_{n}\right\\}\mathbb{r}_{n}\\\
\end{split}$
$\begin{split}\nabla_{i}U(\mathbb{x})&=\sum_{n=1}^{N}\left\\{\frac{\exp(\mathbb{r}_{n}^{T}\mathbb{x})}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}-s_{n}\right\\}r_{n,i}\\\
\end{split}$
$\begin{split}\big{|}\nabla_{i}U(\mathbb{x})\big{|}&=\Big{|}\sum_{n=1}^{N}\left\\{\frac{\exp(\mathbb{r}_{n}^{T}\mathbb{x})}{1+\exp(\mathbb{r}_{n}^{T}\mathbb{x})}-s_{n}\right\\}r_{n,i}\Big{|}\leq\sum_{n=1}^{d}\big{|}r_{n,i}\big{|}\end{split}$

For the Zigzag Sampler, the event rate is

$\lambda_{i}(\mathbb{x}+t\mathbb{v},\mathbb{v})\leq\sum_{n=1}^{N}\big{|}r_{n,i}\big{|}+\lambda_{i}^{\text{ref}}$

For the Coordinate Sampler, if $v_{i}\neq 0$ of $\mathbb{v}$, then

$\lambda(\mathbb{x}+t\mathbb{v},\mathbb{v})\leq\sum_{n=1}^{N}\big{|}r_{n,i}\big{|}+\lambda^{\text{ref}}$

###  8.4 Log-Cox Gaussian Model

The energy function is

$U(\mathbf{x})=\sum_{i,j}\left(-y_{i,j}x_{i,j}+s\exp\\{x_{ij}\\}\right)+\frac{1}{2}(\mathbf{x}-\mu\mathbf{1})^{T}\Sigma^{-1}(\mathbf{x}-\mu\mathbf{1})$
$\nabla_{ij}U(\mathbf{x})=-y_{ij}+s\exp\\{x_{ij}\\}+(\Sigma^{-1})_{ij\cdot}(\mathbf{x}-\mu\mathbf{1})$

Denote $x_{k}=x_{30*(i-1)+j}$, $B=\Sigma^{-1}$,

$\nabla_{k}U(\mathbf{x})=-y_{k}+s\exp\\{x_{k}\\}+(\Sigma^{-1})_{k\cdot}(\mathbf{x}-\mu\mathbf{1})$

For Coordinate sampler:

$\displaystyle\lambda_{k}(\mathbf{x}+t\mathbf{v},\mathbf{v})=\left\\{v_{k}\nabla_{k}U(\mathbf{x}+t\mathbf{v})\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{-y_{k}v_{k}+sv_{k}\exp\\{x_{k}+tv_{k}\\}+v_{k}\sum_{\ell\neq
k}b_{k\ell}(x_{\ell}-\mu)+b_{kk}(x_{k}+tv_{k}-\mu)v_{k}\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{-y_{k}v_{k}+sv_{k}\exp\\{x_{k}+tv_{k}\\}+v_{k}\sum_{\ell}b_{k\ell}(x_{\ell}-\mu)+tb_{kk}v_{k}^{2}\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{v_{k}\left(\sum_{\ell}b_{k\ell}(x_{\ell}-\mu)-y_{k}\right)+tb_{kk}v_{k}^{2}+sv_{k}\exp\\{x_{k}\\}\exp\\{tv_{k}\\}\right\\}_{+}+\lambda_{0}$
$\displaystyle\leq\left\\{v_{k}\left(\sum_{\ell}b_{k\ell}(x_{\ell}-\mu)-y_{k}\right)+tb_{kk}v_{k}^{2}\right\\}_{+}+se^{x_{k}}e^{tv_{k}}\left\\{v_{k}\right\\}_{+}+\lambda_{0}$

For Zig-Zag sampler:

$\displaystyle\lambda_{k}(\mathbf{x}+t\mathbf{v},\mathbf{v})=\left\\{v_{k}\nabla_{k}U(\mathbf{x}+t\mathbf{v})\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{-y_{k}v_{k}+sv_{k}\exp\\{x_{k}+tv_{k}\\}+v_{k}(\Sigma^{-1})_{k\cdot}(\mathbf{x}+t\mathbf{v}-\mu\mathbf{1})\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{-y_{k}v_{k}+sv_{k}\exp\\{x_{k}+tv_{k}\\}+v_{k}\sum_{\ell=1}^{d}b_{k\ell}(x_{\ell}+tv_{\ell}-\mu)\right\\}_{+}+\lambda_{0}$
$\displaystyle=\left\\{v_{k}\left(\sum_{\ell=1}^{d}b_{k\ell}(x_{\ell}-\mu)-y_{k}\right)+\left(v_{k}\sum_{\ell=1}^{d}b_{k\ell}v_{\ell}\right)t+sv_{k}\exp\\{x_{k}+tv_{k}\\}\right\\}_{+}+\lambda_{0}$
$\displaystyle\leq\left\\{v_{k}\left(\sum_{\ell}b_{k\ell}(x_{\ell}-\mu)-y_{k}\right)+\left(v_{k}\sum_{\ell}b_{k\ell}v_{\ell}\right)t\right\\}_{+}+se^{x_{k}}e^{tv_{k}}\left\\{v_{k}\right\\}_{+}+\lambda_{0}$