# In-and-Out: Algorithmic Diffusion for Sampling Convex Bodies

 Yunbum Kook Santosh S. Vempala Matthew S. Zhang  School of Computer
Science, Georgia Institute of Technology<EMAIL_ADDRESS>School of
Computer Science, Georgia Institute of Technology<EMAIL_ADDRESS>Department of Computer Science, University of Toronto, and Vector Institute,
<EMAIL_ADDRESS>
###### Abstract

We present a new random walk for uniformly sampling high-dimensional convex
bodies. It achieves state-of-the-art runtime complexity with stronger
guarantees on the output than previously known, namely in Rényi divergence
(which implies TV, $\mathcal{W}_{2}$, KL, $\chi^{2}$). The proof departs from
known approaches for polytime algorithms for the problem — we utilize a
stochastic diffusion perspective to show contraction to the target
distribution with the rate of convergence determined by functional
isoperimetric constants of the stationary density.

##  1 Introduction

Generating random samples from a high-dimensional convex body is a basic
algorithmic problem with myriad connections and applications. The core of the
celebrated result of [DFK91], giving a randomized polynomial-time algorithm
for computing the volume of a convex body, was the first polynomial-time
algorithm for uniformly sampling convex bodies. In the decades since, the
study of sampling has led to a long series of improvements in its algorithmic
complexity [LS90, LS93, KLS97, LV06, CV18], often based on uncovering new
mathematical/geometric structure, establishing connections to other fields
(e.g., functional analysis, matrix concentration) and developing new tools for
proving isoperimetric inequalities and analyzing Markov chains. With the
proliferation of data and the increasing importance of machine learning,
sampling has also become an essential algorithmic tool, with applications
needing samplers in very high dimension, e.g., scientific computing [CV16,
Har+17, Koo+22], systems biology [LNP12, Thi+13], differential privacy [MT07,
Mir17] and machine learning [Bin+19, Sta20].

Samplers for convex bodies are based on Markov chains (see Appendix A for a
summary). Their analysis is based on bounding the conductance of the
associated Markov chain, which in turn bounds the mixing rate. Analyzing the
conductance requires combining delicate geometric arguments with (Cheeger)
isoperimetric inequalities for convex bodies. An archetypal example of the
latter is the following: for any measurable partition $S_{1},S_{2},S_{3}$ of a
convex body $\mathcal{K}\subset\mathbb{R}^{d}$, we have

$\operatorname{vol}(S_{3})\geq\frac{d(S_{1},S_{2})}{C_{\mathcal{K}}}\min\\{\operatorname{vol}(S_{1}),\operatorname{vol}(S_{2})\\}\,,$

where $d(\cdot,\cdot)$ is the (minimum) Euclidean distance, and
$C_{\mathcal{K}}$ is an isoperimetric constant of the uniform distribution
over $\mathcal{K}$. (The KLS conjecture posits that
$C_{\mathcal{K}}=\mathcal{O}(1)$ for any convex body $\mathcal{K}$ in
isotropic position, i.e., under the normalization that a random point from
$\mathcal{K}$ has identity covariance). The coefficient $C_{\mathcal{K}}^{2}$
is bounded by the Poincaré constant of the uniform distribution over
$\mathcal{K}$ (and they are in fact asymptotically equal). The classical proof
of conductance uses geometric properties of the random walk at hand to reduce
the analysis to a suitable isoperimetric inequality (see e.g., [LS93, Vem05]).
The end result is a guarantee on the number of steps after which the total
variation distance (TV distance) between the current distribution and the
target is bounded by a desired error parameter. This framework has been widely
used and effective in analyzing an array of candidate samplers, e.g.,
$\mathsf{Ball\ walk}$ [KLS97], $\mathsf{Hit\text{-}and\text{-}Run}$ [Lov99,
LV06], $\mathsf{Riemannian\ Hamiltonian\ Monte\ Carlo}$ [LV18] etc.

One successful approach, studied intensively over the past decade, is based on
diffusion. The basic idea is to first analyze a continuous-time diffusion
process, typically modeled by a _stochastic differential equation_ (SDE), and
then show that a suitable time-discretization of the process, sometimes
together with a Metropolis filter, converges to the desired distribution
efficiently. A major success along this line is the $\mathsf{Unadjusted\
Langevin\ Algorithm}$ and its variants [Bes+95, DT12, Dal17, DMM19, VW19].
These algorithms have strong guarantees for sampling “nice” distributions,
such as ones that are strongly log-concave, or more generally distributions
satisfying isoperimetric inequalities, while also obeying some smoothness
conditions. The analysis of these algorithms is markedly different from the
conductance approach, and typically yields guarantees in stronger metrics such
as the $\mathsf{KL}$-divergence.

Our starting point is the following question:

Can diffusion-based approaches be used for the problem of sampling convex
bodies?

Despite remarkable progress, thus far, constrained sampling problems have
evaded the diffusion approach, except as a high-level analogy (e.g., the
$\mathsf{Ball\ walk}$ can be viewed as a discretization of Brownian motion,
but this alone does not suggest a route for analysis) or with significantly
worse convergence rates (e.g., [Bro+17, BEL18]).

Our main finding is a simple diffusion-based algorithm that can be mapped to a
stochastic process (and, importantly, to a pair of forward and backward
processes), such that the rate of convergence is bounded directly by an
appropriate functional inequality for the target distribution. As a
consequence, for the first time, we obtain clean end-to-end guarantees in the
Rényi divergence (which implies guarantees in other well known quantities such
as $\mathcal{W}_{2},\mathsf{TV},\mathsf{KL},\chi^{2}$ etc.), while giving
state-of-the-art runtime complexity for sampling convex bodies (e.g.,
$\mathsf{Ball\ walk}$ or $\mathsf{Speedy\ walk}$ [LS93, KLS97]). Besides being
a stronger guarantee on the output, Rényi divergence is of particular interest
for differential privacy [Mir17]. Perhaps most interesting is that our proof
approach is completely different from prior work on convex body sampling. In
summary,

  * •

The guarantees hold for the $q$-Rényi divergences while matching the rates of
previous work (prior work only had guarantees in the TV distance).

  * •

The analysis is simple, modular, and easily extendable to several other
settings.

Figure 1.1: Description of uniform samplers: (i) $\mathsf{Ball\ walk}$:
proposes a uniform random point $z$ from $B_{\delta}(x_{1})$, but
$z\notin\mathcal{K}$ so it stays at $x_{1}=x_{2}$. (ii) $\mathsf{Speedy\
walk}$: moves to $x_{2}$ drawn uniformly at random from $\mathcal{K}\cap
B_{\delta}(x_{1})$. (iii) $\mathsf{In\text{-}and\text{-}Out}$: first moves to
$y_{2}$ obtained by taking a Gaussian step from $x_{1}$, and then to $x_{2}$
obtained by sampling the truncated Gaussian
$\mathcal{N}(y_{2},hI_{d})|_{\mathcal{K}}$.

###  1.1 Diffusion for uniform sampling

We propose the following $\mathsf{In\text{-}and\text{-}Out}$111This name
reflects the “geometry” of how the iterates are moving. As we elaborate in
Remark 1, the name ‘proximal sampler’ may be more familiar to those from an
optimization background. sampler for uniformly sampling from $\mathcal{K}$.
Each iteration consists of two steps, one that might leave the body and the
second accepted only if it is (back) in the body.

Input: initial point $x_{0}\sim\pi_{0}$, convex body
$\mathcal{K}\subset\mathbb{R}^{d}$, iterations $T$, threshold $N$, and $h>0$.

Output: $x_{T+1}$.

1: for $i=0,\dotsc,T$ do

2: Sample $y_{i+1}\sim\mathcal{N}(x_{i},hI_{d})$.

3: Repeat: Sample $x_{i+1}\sim\mathcal{N}(y_{i+1},hI_{d})$ until
$x_{i+1}\in\mathcal{K}$ or #attemptsi $\geq N$ (declare Failure).

4: end for

Algorithm 1 $\mathsf{In\text{-}and\text{-}Out}$

It might be illuminating for the reader to compare this algorithm to the well-
studied $\mathsf{Ball\ walk}$ (Algorithm 2); each proposed step is a uniform
random point in a fixed-radius ball around the current point, and is accepted
only if the proposed point is in the body $\mathcal{K}$. In contrast, each
iteration of $\mathsf{In\text{-}and\text{-}Out}$ is a two-step process, where
the first step (Line 2) ignores the boundary of the body, and the second step
(Line 3) is accepted only if a proposal $x_{i+1}$ is a feasible point in
$\mathcal{K}$. We will presently elaborate on the benefits of this variation.

Each successful iteration of the algorithm, i.e., one that is not declared
“Failure”, can be called a proper step. We will see that the number of proper
steps is directly bounded by isoperimetric constants (such as Poincaré and
log-Sobolev) of the target distribution. In fact, this holds quite generally
without assuming the convexity of $\mathcal{K}$. The implementation of an
iteration is based on rejection sampling (Line 3), and our analysis of the
efficiency of this step relies crucially on the convexity of $\mathcal{K}$.
This is reminiscent of the $\mathsf{Speedy\ walk}$ in the literature on convex
body sampling (Algorithm 3), which is used as a tool to analyze proper steps
of the $\mathsf{Ball\ walk}$. We refer the reader to Appendix B for a brief
survey on these and related walks.

This simple algorithm can be interpreted as a composition of “flows” in the
space of measures. This view will allow us to use tools from stochastic
analysis. In particular, we shall demonstrate how to interpret the two steps
of one iteration of $\mathsf{In\text{-}and\text{-}Out}$ as alternating
_forward_ and _backward_ heat flows.

We begin by defining an augmented probability measure on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ by

$\pi(x,y)\propto\exp\bigl{(}-\frac{1}{2h}\|x-y\|^{2}\bigr{)}\,\mathds{1}_{\mathcal{K}}(x)\,.$

We denote by $\pi^{X},\pi^{X|Y}(\cdot|y)$ the marginal distribution of its
first component (resp. conditional distribution given the second component),
and similarly denote by $\pi^{Y},\pi^{Y|X}(\cdot|x)$ for the second component.
In particular, the marginal in the first component $\pi^{X}$ is the uniform
distribution over $\mathcal{K}$. Sampling from such a joint distribution to
obtain the marginal on $X$ (say), can be more efficient than directly working
only with $\pi^{X}$. This idea was utilized in Gaussian Cooling [CV18] and
later as the restricted Gaussian Oracle (RGO) [LST21, Che+22].

Under this notation, Algorithm 1 corresponds to a Gibbs sampling scheme from
the two marginals of $\pi(x,y)$. To be precise, Line 2 and Line 3 correspond
to sampling from

$y_{i+1}\sim\pi^{Y|X}(\cdot|x_{i})=\mathcal{N}(x_{i},hI_{d})\qquad\mbox{ and
}\qquad
x_{i+1}\sim\pi^{X|Y}(\cdot|y_{i+1})=\mathcal{N}(y_{i+1},hI_{d})|_{\mathcal{K}}\,.$

We implement the latter step through rejection sampling; if the number of
trials in Line 3 hits the threshold $N$, then we halt and declare _failure_ of
the algorithm. It is well known that such a Gibbs sampling procedure will
ensure the desired stationarity of $\pi(x,y)$.

#### Stochastic perspective: forward and backward heat flows.

Our algorithm can be viewed through the lens of stochastic analysis, due to an
improved analysis for the proximal sampling [Che+22]. This view provides an
interpolation in continuous-time, which is simple and powerful. To make this
concrete, we borrow an exposition from [Che24, Chapter 8.3]. Let us denote the
successive laws of $x_{i}$ and $y_{i}$ by $\mu^{X}_{i}$ and $\mu^{Y}_{i}$,
respectively. Recall that the first step of sampling from
$\pi^{Y|X}(\cdot|x_{i})$ (Line 2) yields
$\mu^{Y}_{i+1}=\int\pi^{Y|X=x}\,\mathrm{d}\mu^{X}_{i}(x)$. This is the result
of evolving a probability measure under _(forward) heat flow_ of $\mu_{i}^{X}$
for some time $h$, given precisely by the following stochastic differential
equation: for $Z_{0}\sim\mu_{i}^{X}$,

$\mathrm{d}Z_{t}=\mathrm{d}B_{t}$ ($\mathsf{FH}$)

where $B_{t}$ is the standard Brownian process. We write
$\mathsf{law}(Z_{t})=\mu_{i}^{X}P_{t}$. In particular,
$Z_{h}=Z_{0}+\zeta\sim\mu_{i}^{X}*\mathcal{N}(0,hI_{d})=\mu_{i+1}^{Y}$ for
$\zeta\sim\mathcal{N}(0,hI_{d})$. When $\mu_{i}^{X}=\pi^{X}$, the first step
of Algorithm 1 gives

$\pi^{Y}(y)=[\pi^{X}*\mathcal{N}(0,hI_{d})](y)=\frac{1}{\operatorname{vol}(\mathcal{K})\,(2\pi
h)^{d/2}}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\,\mathrm{d}x\,.$
(1.1)

The second step of sampling from $\pi^{X|Y}(\cdot|y_{i+1})$ can be represented
by $\mu_{i+1}^{X}=\int\pi^{X|Y=y}\,\mathrm{d}\mu_{i+1}^{Y}(y)$ (Line 3). The
continuous-time process corresponding to this step might not be obvious.
However, let us consider ($\mathsf{FH}$) with $Z_{0}\sim\pi^{X}$. Then,
$Z_{h}\sim\pi^{Y}$, so the joint distribution of $(Z_{0},Z_{h})$ is simply
$\pi$. This implies that $(Z_{0}|Z_{h}=y)\sim\pi^{X|Y=y}$. Imagine there is an
SDE _reversing_ the forward heat flow in a sense that if we are initialized
deterministically at $Z_{h}\sim\delta_{y}$ at time $0$, then the law of the
SDE at time $h$ would be $\mathsf{law}(Z_{0}|Z_{h}=y)=\pi^{X|Y=y}$. Then, this
SDE would serve as a continuous-time interpolation of the second step.

|  Forward flow | Backward flow
---|---|---
SDE | $\mathrm{d}Z_{t}=\mathrm{d}B_{t}$ | $\mathrm{d}Z_{t}^{\leftarrow}=\nabla\log(\pi^{X}P_{h-t})(Z_{t}^{\leftarrow})\,\mathrm{d}t+\mathrm{d}B_{t}$
Fokker-Planck | $\partial_{t}\mu_{t}=\frac{1}{2}\Delta\mu_{t}$ | $\partial_{t}\mu_{t}^{\leftarrow}=-\operatorname{div}\bigl{(}\mu_{t}^{\leftarrow}\nabla\log(\pi^{X}P_{h-t})\bigr{)}+\frac{1}{2}\Delta\mu_{t}^{\leftarrow}$
Table 1: The Fokker-Planck equations for the forward and backward heat flow
describe how the laws of $Z_{t}$ and $Z_{t}^{\leftarrow}$ in ($\mathsf{FH}$)
and ($\mathsf{BH}$) evolve over time. See Section 3.2 for details.

Such a time reversal SDE is indeed possible! The following SDE on
$(Z_{t}^{\leftarrow})_{t\in[0,h]}$ initialized at
$Z_{0}^{\leftarrow}\sim\pi^{Y}=\pi^{X}P_{h}$ ensures
$Z_{h-t}\sim\mathsf{law}(Z_{t}^{\leftarrow})=\pi^{X}P_{h-t}$:

$\mathrm{d}Z_{t}^{\leftarrow}=\nabla\log(\pi^{X}P_{h-t})(Z_{t}^{\leftarrow})\,\mathrm{d}t+\mathrm{d}B_{t}\quad\text{for
}t\in[0,h]\,.$ ($\mathsf{BH}$)

Although this is designed to reverse ($\mathsf{FH}$) initialized by
$Z_{0}\sim\pi^{X}$ (so $Z_{h}=Z_{0}^{\leftarrow}\sim\pi^{Y})$, its
construction also ensures that if $Z_{0}^{\leftarrow}\sim\delta_{y}$, a point
mass, then $Z_{h}^{\leftarrow}\sim\mathsf{law}(Z_{0}|Z_{h}=y)=\pi^{X|Y=y}$.
Thus, if we initialize ($\mathsf{BH}$) with
$Z_{0}^{\leftarrow}\sim\mu_{i+1}^{Y}$, then the law of $Z_{h}^{\leftarrow}$
corresponds to $\int\pi^{X|Y=y}\,\mathrm{d}\mu_{i+1}^{Y}(y)=\mu_{i+1}^{X}$.

######  Remark 1.

We note that $\mathsf{In\text{-}and\text{-}Out}$ is exactly the proximal
sampling scheme [LST21, Che+22, FYC23] for uniform distributions. The proximal
sampler with a target density proportional to $\exp(-V(x))$ considers an
augmented distribution $\pi(x,y)\propto\exp(-V(x)-\frac{1}{2h}\|x-y\|^{2})$
and then repeats the following two steps: (1)
$y_{i+1}\sim\pi^{Y|X=x_{i}}=\mathcal{N}(x_{i},hI_{d})$ and then (2)
$x_{i+1}\sim\pi^{X|Y=y_{i+1}}$.

###  1.2 Results

Our model of computation is the classical general model for convex bodies
[GLS88]. We assume $\operatorname{vol}(\mathcal{K})>0$ throughout this paper.
Below, $B_{r}(x)$ denotes the $d$-dimensional ball of radius $r$ centered at
$x$.

######  Definition 1 (Convex body oracle).

A _well-defined membership oracle_ for a convex body
$\mathcal{K}\subset\mathbb{R}^{d}$ is given by a point $x_{0}\in\mathcal{K}$,
a number $D>0$, with the guarantee that
$B_{1}(x_{0})\subseteq\mathcal{K}\subseteq B_{D}(x_{0})$, and an oracle that
correctly answers _YES_ or _NO_ to any query of the form “$x\in\mathcal{K}$?”

######  Definition 2 (Warmness).

A distribution $\mu$ is _$M$ -warm_ with respect to another distribution
$\pi$, i.e., for every $x$ in the support of $\pi$, we have
$\mathrm{d}\mu(x)\leq M\,\mathrm{d}\pi(x)$.

We now summarize our main result which is further elaborated in Section 3.4
(Theorem 27). Below, $\pi^{\mathcal{K}}$ is the uniform distribution over
$\mathcal{K}$, and $\mathcal{R}_{q}$ is the Rényi-divergence of order $q$ (see
Definition 6).

######  Theorem 3.

For any given $\eta,\varepsilon\in(0,1)$, $q\geq 1$, and any convex body
$\mathcal{K}$ given by a well-defined membership oracle, there exist choices
of parameters $h,N$ such that $\mathsf{In\text{-}and\text{-}Out}$, starting
from an $M$-warm distribution, with probability at least $1-\eta$, returns
$X\sim\mu$ such that
$\mathcal{R}_{q}(\mu\mathbin{\|}\pi^{\mathcal{K}})\leq\varepsilon$. The number
of proper steps is
$\widetilde{\mathcal{O}}(qd^{2}\Lambda\log\frac{M}{\eta\varepsilon})$, and the
expected total number of membership queries is
$\widetilde{\mathcal{O}}(qMd^{2}\Lambda\log^{5}(1/\eta\varepsilon))$, where
$\Lambda$ is the largest eigenvalue of the covariance of $\pi^{\mathcal{K}}$.

We note that for $X\sim\pi^{\mathcal{K}}$,

$\|\mathrm{Cov}(\pi^{\mathcal{K}})\|_{\mathsf{op}}\leq\operatorname{tr}(\mathrm{Cov}(\pi^{\mathcal{K}}))=\mathbb{E}\bigl{[}\|X-\mathbb{E}[X]\|^{2}\bigr{]}\leq
D^{2}\,.$

The above guarantee in the Rényi divergence immediately provides
$\mathcal{W}_{2},\mathsf{TV},\mathsf{KL}$, and $\chi^{2}$ guarantees as
special cases. Previous guarantees for uniformly sampling convex bodies were
only in the $\mathsf{TV}$-distance. For two distributions $\mu$ and $\pi$, we
have

  1. 1.

$\mathsf{KL}(\mu\mathbin{\|}\pi)=\lim_{q\to
1}\mathcal{R}_{q}(\mu\mathbin{\|}\pi)\leq\mathcal{R}_{q}(\mu\mathbin{\|}\pi)\leq\mathcal{R}_{q^{\prime}}(\mu\mathbin{\|}\pi)\leq\mathcal{R}_{\infty}(\mu\mathbin{\|}\pi)=\log\sup\frac{\mathrm{d}\mu}{\mathrm{d}\pi}$
for $1<q\leq q^{\prime}$.

  2. 2.

$2\,\|\mu-\pi\|_{\mathsf{TV}}^{2}\leq\mathsf{KL}(\mu\mathbin{\|}\pi)\leq\log(\chi^{2}(\mu\mathbin{\|}\pi)+1)=\mathcal{R}_{2}(\mu\mathbin{\|}\pi)$.

  3. 3.

$\mathcal{W}_{2}^{2}(\mu,\pi)\leq
2C_{\mathsf{LSI}}(\pi)\,\mathsf{KL}(\mu\mathbin{\|}\pi)$ (Talagrand’s
$\mathsf{T}_{2}$-inequality) and $C_{\mathsf{LSI}}(\pi^{\mathcal{K}})\lesssim
D^{2}$.

  4. 4.

$\mathcal{W}_{2}^{2}(\mu,\pi)\leq
2C_{\mathsf{PI}}(\pi)\,\chi^{2}(\mu\mathbin{\|}\pi)$ [Liu20] and
$C_{\mathsf{PI}}(\pi^{\mathcal{K}})\lesssim\|\mathrm{Cov}(\pi^{\mathcal{K}})\|_{\mathsf{op}}\log
d$.

The query complexity is better if the convex body is (near-)_isotropic_ ,
i.e., the uniform distribution over the body has (near-)identity covariance.
This relies on recent estimates of the worst-case Poincaré constant for
isotropic log-concave distributions [KLS95, Kla23].

######  Corollary 4.

Assume that $\pi^{\mathcal{K}}$ is isotropic. Under the same setting as above,
$\mathsf{In\text{-}and\text{-}Out}$ succeeds with probability $1-\eta$,
returning $X\sim\mu$ such that
$\mathcal{R}_{q}(\mu\mathbin{\|}\pi^{\mathcal{K}})\leq\varepsilon$. The number
of proper steps is
$\widetilde{\mathcal{O}}(qd^{2}\log\frac{M}{\eta\varepsilon})$, and the
expected total number of membership queries is
$\widetilde{\mathcal{O}}(qMd^{2}\log^{5}(1/\eta\varepsilon))$.

Our analysis will in fact show that the bound on the number of proper steps
holds for general _non-convex bodies_ and _any feasible start_ in
$\mathcal{K}$. This is deduced under an $M$-warm start in Corollaries 28 and
29. We remark that such a bound for non-convex uniform sampling is not known
for the $\mathsf{Ball\ walk}$ or the $\mathsf{Speedy\ walk}$.

######  Theorem 5.

For any given $\varepsilon\in(0,1)$ and set $\mathcal{K}\subset B_{D}(0)$ with
$\operatorname{vol}(\mathcal{K})>0$, $\mathsf{In\text{-}and\text{-}Out}$ with
variance $h$ and $M$-warm initial distribution achieves
$\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$ after the
following number of iterations:

$m=\min\begin{cases}\mathcal{O}\bigl{(}qh^{-1}C_{\mathsf{PI}}(\pi^{X})\log\frac{M}{\varepsilon}\bigr{)}&\text{for
}q\geq 2\,,\\\
\mathcal{O}\bigl{(}qh^{-1}C_{\mathsf{LSI}}(\pi^{X})\log\frac{\log
M}{\varepsilon}\bigr{)}&\text{for }q\geq 1\,.\end{cases}$

We have two different convergence results above under (8) and (7). Under (8)
we have a _doubly-logarithmic_ dependence on the warmness parameter $M$. On
the other hand, using (7), which is weaker than (8) (in general,
$C_{\mathsf{PI}}\leq C_{\mathsf{LSI}}$), the dependence on $M$ is logarithmic.
We discuss our results further in Section 1.3.

#### Outline of analysis.

We summarize our proof strategy below, which consists of two steps: (i) The
current distribution should converge to the uniform distribution, (ii) within
each iteration of the algorithm, the failure probability and the expected
number of rejections should be small enough.

  * •

We need to demonstrate that the corresponding Markov chain is rapidly mixing.
Here, we use the heat flow perspective to derive mixing rates under any
suitable divergence measure (such as $\mathsf{KL}$, $\chi^{2}$, or
$\mathcal{R}_{q}$), extending known results for the unconstrained setting
[Che+22]. As a result, the mixing analysis reduces to a suitable functional
inequality of the target distribution alone.

  * •

We show that the number of rejections in Line 3 over the entire execution of
the algorithm is bounded with high probability. To do this, we apply a
detailed argument involving local conductance and the convexity of
$\mathcal{K}$, which relies on techniques from [BNN06]. For this step, we show
that with the appropriate choice of variance $h=\widetilde{\Theta}(d^{-2})$
and threshold $N=\widetilde{\Theta}(T\eta^{-1})$, the entire algorithm
succeeds with probability $1-\eta$. The expected number of rejections is
polylogarithmic.

While each individual component resembles pre-existing work in the literature,
in their synthesis we will demonstrate how to interleave past developments in
theoretical computer science, optimal transport, and functional analysis. The
combination of these in this domain yields elegant and surprisingly simple
proofs, as well as stronger results.

###  1.3 Discussion

Here we make a few remarks, contrasting our results with known ones.

#### No need to be lazy.

Previous uniform samplers like the $\mathsf{Ball\ walk}$ are made _lazy_
(i.e., with probability $1/2$, the Markov chain does nothing), to ensure
convergence to the target stationary distribution. However, our algorithm does
not need this, since we directly show that our sampler contracts towards a
uniform distribution.

#### Unified framework.

We remark that these two different bounds also place the previously known
mixing guarantees for $\mathsf{Ball\ walk},\mathsf{Speedy\ walk}$ in a unified
framework. Existing tight guarantees for $\mathsf{Speedy\ walk}$ are in TV
distance and based on the log-Sobolev constant, assuming an oracle for
implementing each step [LV17]. The known convergence guarantees of
$\mathsf{Ball\ walk}$ (see Appendix B for details), namely the mixing time of
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log\frac{1}{\varepsilon})$ for TV
distance, are for the composite algorithm [$\mathsf{Speedy\ walk}+$rejection
sampling]. Here $\mathsf{Speedy\ walk}$ records only the accepted steps of
$\mathsf{Ball\ walk}$, so its stationary distribution differs slightly from
the uniform distribution (and can be corrected with a post-processing step).
On the other hand, $\mathsf{In\text{-}and\text{-}Out}$ actually converges to
$\pi^{\mathcal{K}}$ without any adjustments and achieves stronger Rényi
divergence bounds in the same asymptotic complexity. Our analysis shows that
the mixing guarantee is determined by isoperimetric constants of the target
(Poincaré or log-Sobolev).

#### Effective step size.

The $\mathsf{Ball\ walk}$’s largest possible step size is of order
$1/\sqrt{d}$ (see Appendix B) to keep the rejection probability bounded by a
constant. This bound could also be viewed as an “effective” step size of
$\mathsf{In\text{-}and\text{-}Out}$. This follows from the fact that the
$\ell_{2}$-norm of the Gaussian $\mathcal{N}(0,hI)$ is concentrated around
$\sqrt{hd}$, and we will set the variance $h$ of
$\mathsf{In\text{-}and\text{-}Out}$ to $\widetilde{\mathcal{O}}(1/d^{2})$, so
we have $\sqrt{hd}\approx 1/\sqrt{d}$.

#### What has really changed?

$\mathsf{In\text{-}and\text{-}Out}$ has clear similarities to both
$\mathsf{Ball\ walk}$ and $\mathsf{Speedy\ walk}$. What then are the changes
that allow us to use continuous-time interpolation? One step of $\mathsf{Ball\
walk}$ is [random step ($y\in B_{\delta}(x)$) $+$ Metropolis-filter (accept if
$y\in\mathcal{K}$)]. This filtering is an abrupt discrete step, and it is
unclear how to control contraction. It could be replaced by a step of
$\mathsf{Speedy\ walk}$ ($x\sim\text{Unif}(B_{\delta}(y)\cap\mathcal{K})$).
Then, each iteration of $\mathsf{In\text{-}and\text{-}Out}$ can be viewed as a
Gaussian version of a $\mathsf{Ball\ walk^{\prime}s\ proposal}+\mathsf{Speedy\
walk}$ algorithm.

How can we compare $\mathsf{In\text{-}and\text{-}Out}$ with $\mathsf{Speedy\
walk}?$ Iterating speedy steps leads to a biased distribution (one that is
proportional to the local conductance). As clarified in Remark 2, one step of
(a Gaussian version of) $\mathsf{Speedy\ walk}$ can be understood as a step of
backward heat flow. Therefore, if one can control the isoperimetric constants
of the biased distribution along the trajectory of the backward flow, then
contraction of $\mathsf{Speedy\ walk}$ toward the biased distribution will
follow from the simultaneous backward analysis.

##  2 Preliminaries

Unless otherwise specified, we will use $\|\cdot\|$ for the $2$-norm on
$\mathbb{R}^{d}$, and operator norm on $\mathbb{R}^{d\times d}$. We write
$a=\mathcal{O}(b)$, $a\lesssim b$ to mean that $a\leq cb$ for some universal
constant $c>0$. Similarly, $a\gtrsim b,a=\Omega(b)$ for $a\geq cb$, while
$a=\Theta(b)$ means that $a\lesssim b,b\lesssim a$ simultaneously. We will
also use $a=\widetilde{\mathcal{O}}(b)$ to denote
$a=\mathcal{O}(b\,\mathrm{polylog}(b))$. Lastly, we will use measure and
density interchangeably when there is no confusion.

To quantify the convergence rate, we introduce some common divergences between
distributions.

######  Definition 6 (Distance and divergence).

For two measures $\mu,\nu$ on $\mathbb{R}^{d}$, the _total variation_ distance
between them is defined by

$\|\mu-\nu\|_{\mathsf{TV}}:=\sup_{B\in\mathcal{F}}\lvert\mu(B)-\nu(B)\rvert\,,$

where $\mathcal{F}$ is the collection of all measurable subsets of
$\mathbb{R}^{d}\,$. The $2$-_Wasserstein distance_ is given by

$\mathcal{W}_{2}^{2}(\mu,\nu):=\inf_{\gamma\in\Gamma(\mu,\nu)}\mathbb{E}_{(X,Y)\sim\gamma}[\|X-Y\|^{2}]\,,$

where $\Gamma$ is the set of all couplings between $\mu,\nu$. Next, we define
the _$f$ -divergence_ of $\mu$ towards $\nu$ with $\mu\ll\nu$ (i.e., $\mu$ is
absolutely continuous with respect to $\nu$) as, for some convex function
$f:\mathbb{R}_{+}\to\mathbb{R}$ with $f(1)=0$ and $f^{\prime}(\infty)=\infty$,

$D_{f}(\mu\mathbin{\|}\nu):=\int
f\Bigl{(}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\Bigr{)}\,\mathrm{d}\nu\,.$

The _$\mathsf{KL}$ -divergence_ arises when taking $f(x)=x\log x$, the
_$\chi^{q}$ -divergence_ when taking $f(x)=x^{q}-1$, and the _$q$ -Rényi
divergence_ is given by

$\mathcal{R}_{q}(\mu\mathbin{\|}\nu):=\frac{1}{q-1}\log\bigl{(}\chi^{q}(\mu\mathbin{\|}\nu)+1\bigr{)}\,.$

We recall two important functional inequalities of a distribution. We use
$\nu$ to denote a probability measure over $\mathbb{R}^{d}$.

######  Definition 7.

We say that $\nu$ satisfies a _Poincaré inequality_ (PI) with parameter
$C_{\mathsf{PI}}(\nu)$ if for all smooth functions
$f:\mathbb{R}^{d}\to\mathbb{R}$,

$\mathrm{Var}_{\nu}f\leq C_{\mathsf{PI}}(\nu)\,\mathbb{E}_{\nu}[\|\nabla
f\|^{2}]\,,$

where $\mathrm{Var}_{\nu}f\coloneqq\mathbb{E}_{\nu}\lvert
f-\mathbb{E}_{\nu}f\rvert^{2}$.

The Poincaré inequality is implied by the log-Sobolev inequality.

######  Definition 8.

We say that $\nu$ satisfies a _log-Sobolev inequality_ (LSI) with parameter
$C_{\mathsf{LSI}}(\nu)$ if for all smooth functions
$f:\mathbb{R}^{d}\to\mathbb{R}$,

$\mathsf{Ent}_{\nu}(f^{2})\leq
2C_{\mathsf{LSI}}(\nu)\,\mathbb{E}_{\nu}[\|\nabla f\|^{2}]\,,$

where $\mathsf{Ent}_{\nu}(f^{2}):=\mathbb{E}_{\nu}[f^{2}\log
f^{2}]-\mathbb{E}_{\nu}[f^{2}]\log(\mathbb{E}_{\nu}[f^{2}])$. Equivalently,
for any probability measure $\mu$ over $\mathbb{R}^{d}$ with $\mu\ll\nu$,

$\mathsf{KL}(\mu\mathbin{\|}\nu)\leq\frac{C_{\mathsf{LSI}}(\nu)}{2}\,\mathsf{FI}(\mu\mathbin{\|}\nu)\,,$

where
$\mathsf{FI}(\mu\mathbin{\|}\nu):=\mathbb{E}_{\mu}[\|\nabla\log\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\|^{2}]$
is the _Fisher information_ of $\mu$ with respect to $\nu$.

We recall the data-processing inequality for Rényi divergence and
$f$-divergence.

######  Lemma 9 (Data-processing inequality).

For measures $\mu,\nu$, Markov kernel $P$, $f$-divergence $D_{f}$, and $q\geq
1$, it holds that

$D_{f}(\mu P\mathbin{\|}\nu P)\leq
D_{f}(\mu\mathbin{\|}\nu)\,,\quad\text{and}\quad\mathcal{R}_{q}(\mu
P\mathbin{\|}\nu P)\leq\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\,.$

The aforementioned functional inequalities allow us to show exponential
contraction of various divergences, through the following helpful inequality.

######  Lemma 10 (Grönwall).

Suppose that $u,g:[0,T]\to\mathbb{R}$ are two continuous functions, with $u$
being differentiable on $[0,T]$ and satisfying

$u^{\prime}(t)\leq g(t)\,u(t)\qquad\text{for all }t\in[0,T]\,.$

Then,

$u(t)\leq\exp\Bigl{(}\int_{0}^{t}g(s)\,\mathrm{d}s\Bigr{)}\,u(0)\qquad\text{for
all }t\in[0,T]\,.$

##  3 Analysis

We begin this section by proving some introductory lemmas, which will
streamline our later exposition.

We first record two fundamental lemmas, which introduce the mathematical
formalism for our analysis. The first is the existence of forward and backward
heat flows (Lemma 22), which will interpolate each line in Algorithm 1. These
flow equations describe how the laws of $Z_{t}$ and $Z_{t}^{\leftarrow}$ in
($\mathsf{FH}$) and ($\mathsf{BH}$) evolve respectively over time.

######  Lemma 11.

The forward heat flow equation with initial distribution $\mu_{0}$ is given by

$\partial_{t}\mu_{t}=\frac{1}{2}\Delta\mu_{t}\,,$

and its backward heat flow equation is given by

$\partial_{t}\mu_{t}^{\leftarrow}=-\operatorname{div}\bigl{(}\mu_{t}^{\leftarrow}\nabla\log(\pi^{X}P_{h-t})\bigr{)}+\frac{1}{2}\Delta\mu_{t}^{\leftarrow}\quad\text{with
}\mu_{0}^{\leftarrow}=\mu_{h}\,.$

These admit (weak) solutions on $[0,h]$ for any initial distribution $\mu_{0}$
with $\frac{\mathrm{d}\mu_{0}}{\mathrm{d}\pi^{X}}\leq M<\infty$.

One successful iteration of $\mathsf{In\text{-}and\text{-}Out}$ is exactly the
same as the composition of running the forward heat flow and then backward
heat flow, both for time $h$. Its validity is justified in Section 1.1, and we
give its proof below.

######  Lemma 12.

Let $\mu_{k}^{X}$ be the law of the $k$-th iterate $x_{k}$ of
$\mathsf{In\text{-}and\text{-}Out}$. If ($\mathsf{FH}$) is initialized with
$\mathsf{law}(Z_{0})=\mu_{k}^{X}$, then $\mathsf{law}(Z_{h})=\mu_{k+1}^{Y}$.
If ($\mathsf{BH}$) is initialized with
$\mathsf{law}(Z_{0}^{\leftarrow})=\mu_{k+1}^{Y}$, then
$\mathsf{law}(Z_{h}^{\leftarrow})=\mu_{k+1}^{X}$.

###### Proof.

We explicitly show that the forward and backward heat flows indeed interpolate
the two discrete steps given in Algorithm 1. For the forward part, we have
$Z_{h}=Z_{0}+\zeta$ for $\zeta\sim\mathcal{N}(0,hI_{d})$, so

$\mathsf{law}(Z_{h})=\mathsf{law}(Z_{0})*\mathcal{N}(0,hI_{d})=\mu_{k}^{X}*\mathcal{N}(0,hI_{d})=\mu_{k+1}^{Y}\,.$

Regarding the backward part, it is known from [Che+22, Lemma 14] that the
construction of the time-reversal SDE ensures that
$(Z_{h}^{\leftarrow},Z_{0}^{\leftarrow})$ and $(Z_{0},Z_{h})$ have the same
joint distribution, when $Z_{0}\sim\pi^{X}$ (and so $Z_{h}\sim\pi^{Y}$).
Hence,
$\mathsf{law}(Z_{h}^{\leftarrow}|Z_{0}^{\leftarrow}=y)=\mathsf{law}(Z_{0}|Z_{h}=y)=\pi^{X|Y=y}$,
where the last equality follows from $(Z_{0},Z_{h})\sim\pi$. Since we
initialize ($\mathsf{BH}$) with $Z_{0}^{\leftarrow}=y\sim\mu_{k+1}^{Y}$, we
have

$\mathsf{law}(Z_{h}^{\leftarrow})=\int\mathsf{law}(Z_{h}^{\leftarrow}|Z_{0}^{\leftarrow}=y)\,\mu_{k+1}^{Y}(\mathrm{d}y)=\int\pi^{X|Y}(\cdot|y)\,\mu_{k+1}^{Y}(\mathrm{d}y)=\mu_{k+1}^{X}\,,$

where the last follows from the definition of Line 3. ∎

Our analysis of these objects consists of two parts: (1) demonstrating the
mixing of $\mathsf{In\text{-}and\text{-}Out}$, i.e., how many _outer_
iterations are needed to be sufficiently close to the uniform distribution,
and (2) quantifying the failure probability and wasted steps in Line 3.

For (1), we collect in Section 3.1 some important implications of functional
inequalities, e.g. the Poincaré and log-Sobolev inequalities, for the uniform
distribution. Then in Section 3.2, we exploit the flow perspective of the
algorithm to obtain the mixing guarantees. To this end, we revisit the proofs
for the contraction results of forward and backward heat flows, and then
extend them to our _constrained_ setting:

######  Theorem 13 (Adapted for uniform distributions, [Che+22]).

Let $\mu^{X}_{k}$ be the law of the $k$-th output of
$\mathsf{In\text{-}and\text{-}Out}$ with initial distribution $\mu_{0}^{X}$.
Let $C_{\mathsf{LSI}}$ be the (8) constant of the uniform distribution
$\pi^{X}$ over $\mathcal{K}$. Then, for any $q\geq 1$,

$\mathcal{R}_{q}(\mu_{k}^{X}\mathbin{\|}\pi^{X})\leq\frac{\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})}{(1+h/C_{\mathsf{LSI}})^{2k/q}}\,.$

For $C_{\mathsf{PI}}$ the (7) constant of $\pi^{X}$ and any $q\geq 2$,

$\mathcal{R}_{q}(\mu_{k}^{X}\mathbin{\|}\pi^{X})\leq\begin{cases}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-\frac{2k\log(1+h/C_{\mathsf{PI}})}{q}&\text{if
}k\leq\frac{q}{2\log(1+h/C_{\mathsf{PI}})}\,\bigl{(}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-1\bigr{)}\,,\\\
(1+h/C_{\mathsf{PI}})^{-2(k-k_{0})/q}&\text{if }k\geq
k_{0}:=\lceil\frac{q}{2\log(1+h/C_{\mathsf{PI}})}\,\bigl{(}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-1\bigr{)}\rceil\,.\end{cases}$

This reduces the problem of obtaining a mixing guarantee to that of
demonstrating a functional inequality on the target distribution. For this, it
is not strictly necessary that $\mathcal{K}$ be convex.

On the other hand, convexity of $\mathcal{K}$ is crucial for the proof of (2).
We show in Section 3.3 that the failure probability remains under control by
taking a suitable variance $h$ and threshold $N$. We then show that the
expected number of trials per iteration is of order $\log N$, not $N$:

######  Lemma 14 (Per-iteration guarantees).

Let $\mathcal{K}$ be any convex body in $\mathbb{R}^{d}$ presented by a well-
defined membership oracle, $\pi^{X}$ the uniform distribution over
$\mathcal{K}$, and $\mu$ an $M$-warm initial distribution with respect to
$\pi^{X}$. For any given $m\in\mathbb{N}$ and $\eta\in(0,1)$, set
$Z=\frac{9mM}{\eta}(\geq 9)$, $h=\frac{\log\log Z}{2d^{2}\log Z}$ and
$N=Z(\log Z)^{4}=\widetilde{\mathcal{O}}(\frac{mM}{\eta})$. Then, the failure
probability of one iteration of $\mathsf{In\text{-}and\text{-}Out}$ is at most
$\eta/m$. Moreover, the expected number of membership queries needed per
iteration is $\mathcal{O}\bigl{(}M(\log\frac{mM}{\eta})^{4}\bigr{)}$.

###  3.1 Functional inequalities

The contraction of an outer loop of our algorithm is controlled by
isoperimetry of the uniform distribution $\pi^{X}$, which is described
precisely by a functional inequality. The most natural ones to consider in
this setting are the Poincaré inequality (7) and log-Sobolev inequality (8).
In Appendix C, we provide a more detailed discussion of how these are related
to other important notions of isoperimetry, such as the _Cheeger_ and _log-
Cheeger_ inequalities.

Below, we use $\mu,\nu$ to denote two arbitrary probability measures over
$\mathbb{R}^{d}$. The relationship between a Poincaré inequality and the
$\chi^{2}$-divergence is derived by substituting
$f=\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$ into (7).

######  Lemma 15.

Assume that $\nu$ satisfies (7) with parameter $C_{\mathsf{PI}}(\nu)$. For any
probability measure $\mu$ over $\mathbb{R}^{d}$ with $\mu\ll\nu$, it holds
that

$\chi^{2}(\mu\mathbin{\|}\nu)\leq\frac{C_{\mathsf{PI}}(\nu)}{2}\,\mathbb{E}_{\nu}\bigl{[}\bigl{\|}\nabla\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\bigr{\|}^{2}\bigr{]}\,.$

The Poincaré inequality implies functional inequalities for the Rényi
divergence.

######  Lemma 16 ([VW19, Lemma 9]).

Assume that $\nu$ satisfies (7) with parameter $C_{\mathsf{PI}}(\nu)$. For any
$q\geq 2$ and probability measure $\mu$ over $\mathbb{R}^{d}$, it holds that

$1-\exp(-\mathcal{R}_{q}(\mu\mathbin{\|}\nu))\leq\frac{q\,C_{\mathsf{PI}}(\nu)}{4}\,\mathsf{RF}_{q}(\mu\mathbin{\|}\nu)\,,$

where
$\mathsf{RF}_{q}(\mu\mathbin{\|}\nu):=q\,\mathbb{E}_{\nu}\bigl{[}\bigl{(}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\bigr{)}^{q}\|\nabla\log\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\|^{2}\bigr{]}\,/\,\mathbb{E}_{\nu}\bigl{[}\bigl{(}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\bigr{)}^{q}\bigr{]}$
is the _Rényi Fisher information_ of order $q$ of $\mu$ with respect to $\nu$.

The log-Sobolev inequality paired with the KL-divergence (8) can be understood
as a special case of the following inequality222Such inequalities are often
called Polyak-Łojasiewicz inequalities, which say for
$f:\mathbb{R}^{d}\to\mathbb{R}$, and all $y\in\mathbb{R}^{d}$ that $f(y)\leq
c\,\|\nabla f(y)\|^{2}$ for some constant $c$, if $\min f(x)=0$. paired with
the $q$-Rényi divergence for $q\geq 1$.

######  Lemma 17 ([VW19, Lemma 5]).

Assume that $\nu$ satisfies (8) with parameter $C_{\mathsf{LSI}}(\nu)$. For
any $q\geq 1$ and probability measure $\mu$ over $\mathbb{R}^{d}$, it holds
that

$\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\leq\frac{q\,C_{\mathsf{LSI}}(\nu)}{2}\,\mathsf{RF}_{q}(\mu\mathbin{\|}\nu)\,.$

Note that $\lim_{q\to 1}\mathcal{R}_{q}=\mathsf{KL}$ and
$\mathsf{RF}_{1}=\mathsf{FI}$.

We have collected below the functional inequalities used to establish the
mixing of our algorithm (see Appendix C for a detailed presentation).

######  Lemma 18.

Let $\mathcal{K}\subset\mathbb{R}^{d}$ be a convex body with diameter $D$, and
$\pi$ be the uniform distribution over $\mathcal{K}$. Then,
$C_{\mathsf{PI}}(\pi)\lesssim\|\mathrm{Cov}(\pi)\|_{\mathsf{op}}\,\log d$ and
$C_{\mathsf{LSI}}(\pi)\lesssim D^{2}$. If $\pi$ is isotropic, then
$C_{\mathsf{PI}}(\pi)\lesssim\log d$ and $C_{\mathsf{LSI}}(\pi)\lesssim D$.

###  3.2 Contraction and mixing

We start by analyzing how many outer iterations of
$\mathsf{In\text{-}and\text{-}Out}$ are required to be $\varepsilon$-close to
$\pi^{X}$, the uniform distribution over $\mathcal{K}$. The contraction of
Algorithm 1 comes from analyzing Lines 2 and 3 through the perspective of heat
flows (see Section 1.1). To exploit this view, we first revisit the previous
contraction analysis in [Che+22], which is carried out for distributions with
_smooth_ densities. Although the uniform distribution is not even continuous,
we prove a technical lemma (Lemma 22) that enables us to extend previously
known results to the uniform distribution. Lastly, combining the previous
results with our technical lemma, we obtain clean contraction results of
Algorithm 1 toward the uniform distribution $\pi^{X}$ in Theorem 23.

#### Part I: Contraction analysis for smooth distributions.

In this part, we review the contraction results for heat flow and its time-
reversal [Che+22], which are intimately connected with our algorithm. We also
provide key technical ingredients needed for its proof, such as the
computations for measures evolving under simultaneous forward/backward heat
flows. We refer interested readers to Appendix D for additional details. Only
in Part I, we assume that $\nu$ denotes a probability measure with smooth
density.

##### Forward heat flow.

We begin by introducing the “heat flow” equation (or also known as the
_Fokker-Planck_ equation), which describes the evolution of the law of $Z_{t}$
under ($\mathsf{FH}$),

$\partial_{t}\mu_{t}=\frac{1}{2}\,\Delta\mu_{t}=\frac{1}{2}\operatorname{div}(\mu_{t}\nabla\log\mu_{t})\,.$
($\mathsf{FP\text{-}FH}$)

It is well known that one can realize this equation in discrete time through a
Gaussian transition density, in the sense that, for $\mu_{h}$ (the solution at
time $h>0$ to ($\mathsf{FP\text{-}FH}$) with initial condition $\mu_{0}$), and
for any smooth function $f:\mathbb{R}^{d}\to\mathbb{R}$,

$\mathbb{E}_{\mu_{h}}[f(x)]=\mathbb{E}_{\mu_{0}}[P_{h}f(x)]\,,$

where $P_{h}f(x)=\mathbb{E}_{\mathcal{N}(x,hI_{d})}[f]$.333$\\{P_{h}\\}_{h\geq
0}$ is often called the heat semigroup. By this we can formally identify
$\mu_{h}=\mu_{0}P_{h}$, and also write $\mu_{h}$ for the law of $Z_{h}$, where
$\\{Z_{h}\\}_{h\geq 0}$ solves ($\mathsf{FH}$).

##### Backward heat flow.

Although there are many ways to define a “reversal” of $P_{h}$, we will use
the notion of _adjoint_ introduced by [KP21], which is the most immediately
useful.

Given some initial measure $\nu$ and some time horizon $h$, the adjoint
corresponds to reversing ($\mathsf{FH}$) for times in $[0,h]$ when the initial
distribution under consideration is $Z_{0}\sim\nu$. For other measures, it
must be interpreted more carefully, and is given by the following partial
differential equation starting from some measure $\mu_{0}^{\leftarrow}$ (see
(D.1) and its derivation):

$\partial_{t}\mu_{t}^{\leftarrow}=-\operatorname{div}\bigl{(}\mu_{t}^{\leftarrow}\nabla\log(\nu
P_{h-t})\bigr{)}+\frac{1}{2}\Delta\mu_{t}^{\leftarrow}\quad\text{for
}t\in[0,h]\,.$ ($\mathsf{FP\text{-}BH}$)

Write $\mu_{t}^{\leftarrow}=\mu_{0}^{\leftarrow}Q_{t}^{\nu,h}$, where
$\\{Q_{t}^{\nu,h}\\}_{t\in[0,h]}$ is a family of transition densities. Write
$\mathbf{P}_{0,h}$ for the joint distribution of the $(Z_{0},Z_{h})$-marginals
of ($\mathsf{FH}$), when $Z_{0}\sim\nu$, and $\mathbf{P}_{0|h}$ for the
conditional. Note that $\mathbf{P}_{h|0}(\cdot|x)=\mathcal{N}(x,hI_{d})$. It
is also known that ($\mathsf{FP\text{-}BH}$) gives a time-reversal of the heat
equation at the SDE level, in the sense that we can interpret
$\delta_{x}Q_{h}^{\nu,h}=\mathbf{P}_{0|h}(\cdot|Z_{h}=x)$. Thus
$\mu_{0}^{\leftarrow}Q_{h}^{\nu,h}=\int\mathbf{P}_{0|h}(\cdot|Z_{h}=x)\,\mu_{0}^{\leftarrow}(\mathrm{d}x)$,
and $\nu P_{h}Q_{t}^{\nu,h}=\nu P_{h-t}$ for all $t\in[0,h]$.

The ultimate purpose of this machinery is to affirm our earlier description of
the Gibbs sampling procedure as alternating forward and backward heat flows.
Indeed, notice that, if $\mu_{i}^{X}$ is the law of the iterate at some
iteration $i$, then $\mu_{i}^{X}P_{h}$ is precisely $\mu_{i+1}^{Y}$ under our
scheme, while $(\mu_{i}^{X}P_{h})Q_{h}^{\pi^{X},h}$ is $\mu_{i+1}^{X}$,
assuming $Q_{h}^{\pi^{X},h}$ is well defined for non-smooth measures
$\pi^{X}$. Thus, while Algorithm 1 is implemented via discrete steps, it can
be exactly analyzed through arguments in continuous time. We shall see the
benefits of this shortly.

Instead of considering the change in metrics along the evolution of $\mu
P_{t}$ with respect to “fixed” $\nu$, it will be useful to consider the
_simultaneous_ evolution of $\mu P_{t},\nu P_{t}$ (and similarly $\mu
Q_{t},(\nu P_{h})Q_{t}$). This type of computation was carried out for
specific metrics in earlier work [VW19, Che+22]. The following is a more
generalized form of one appearing in [Yua+23, Lemma 2]. In the lemma below, we
consider an arbitrary diffusion equation with corresponding Fokker-Planck
equation:

$\mathrm{d}X_{t}=b_{t}(X_{t})\,\mathrm{d}t+\,\mathrm{d}B_{t}\,\quad\text{and}\quad\partial_{t}\mu_{t}=-\nabla\cdot(b_{t}\mu_{t})+\frac{1}{2}\Delta\mu_{t}\,$
(3.1)

where $b_{t}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is smooth,
$X_{t}\in\mathbb{R}^{d}$, and $\mu_{t}=\mathsf{Law}(X_{t})$ if
$X_{0}\sim\mu_{0}$.

######  Lemma 19 (Decay along forward/backward heat flows).

Let $(\mu_{t})_{t\geq 0},(\nu_{t})_{t\geq 0}$ denote the laws of the solutions
to (3.1) starting at $\mu_{0},\nu_{0}$ respectively. Then, for any
differentiable function $g$,

$\partial_{t}g\bigl{(}D_{f}(\mu_{t}\mathbin{\|}\nu_{t})\bigr{)}=-\frac{1}{2}\,g^{\prime}\bigl{(}D_{f}(\mu_{t}\mathbin{\|}\nu_{t})\bigr{)}\times\mathbb{E}_{\mu_{t}}\Bigl{\langle}\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)},\nabla\log\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,.$

###### Proof.

The case where $g\neq\mathsf{id}$ is an application of the chain rule, so it
suffices to take $g=\mathsf{id}$ and simply differentiate an $f$-divergence.

For brevity, we drop the variable $x$ of functions involved, and proceed as
follows:

$\displaystyle\partial_{t}D_{f}(\mu_{t}\mathbin{\|}\nu_{t})$
$\displaystyle=\int\Bigl{\\{}\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\partial_{t}\nu_{t}+\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\bigl{(}\frac{\mu_{t}}{\nu_{t}}\bigr{)}^{\prime}\nu_{t}\Bigr{\\}}\,\mathrm{d}x$
$\displaystyle=\int\Bigl{\\{}\partial_{t}\nu_{t}\Bigl{(}\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}-\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\Bigr{)}+\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\partial_{t}\mu_{t}\Bigr{\\}}\,\mathrm{d}x$
$\displaystyle\underset{(i)}{=}\int\bigl{[}-\nabla\cdot(b_{t}\nu_{t})+\frac{1}{2}\Delta\nu_{t}\bigr{]}\Bigl{(}\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}-\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\Bigr{)}\,\mathrm{d}x$
$\displaystyle\qquad+\int\bigl{[}-\nabla\cdot(b_{t}\mu_{t})+\frac{1}{2}\Delta\mu_{t}\bigr{]}\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\,\mathrm{d}x\,,$

where in $(i)$ we substitute the F-P equation from (3.1). Integrating by parts
(i.e., $\int f\operatorname{div}(\mathbf{G})=-\int\langle\nabla
f,\mathbf{G}\rangle$ for a real-valued function $f$ and vector-valued function
$\mathbf{G}$), we have that

$\int\bigl{[}-\nabla\cdot(b_{t}\nu_{t})\bigr{]}\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\,\mathrm{d}x=\int\Bigl{\langle}b_{t}\nu_{t},\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\nabla\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,\mathrm{d}x\,.$
(3.2)

On the other hand, we have that

$-\int\bigl{[}-\nabla\cdot(b_{t}\nu_{t})\bigr{]}\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\,\mathrm{d}x=-\int\Bigl{\langle}b_{t}\nu_{t},\frac{\mu_{t}}{\nu_{t}}\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}+\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\nabla\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,\mathrm{d}x\,.$

The second term cancels with the RHS of (3.2). We have a similar cancellation
for the $\frac{1}{2}\Delta\nu_{t}$ term:

$\int\frac{1}{2}\Delta\nu_{t}\,\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\,\mathrm{d}x=-\int\frac{1}{2}\Bigl{\langle}\nabla\nu_{t},\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\nabla\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,\mathrm{d}x\,,$

and

$-\int\frac{1}{2}\Delta\nu_{t}\,\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\,\mathrm{d}x=\int\frac{1}{2}\Bigl{\langle}\nabla\nu_{t},\frac{\mu_{t}}{\nu_{t}}\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}+\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\nabla\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,\mathrm{d}x\,.$

Combining these, we are left with

$\displaystyle\int\bigl{[}-\nabla\cdot(b_{t}\nu_{t})+\frac{1}{2}\Delta\nu_{t}\bigr{]}\Bigl{(}\bigl{(}f\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}-\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\Bigr{)}\,\mathrm{d}x$
$\displaystyle=-\int\Bigl{\langle}b_{t}\nu_{t}-\frac{1}{2}\nabla\nu_{t},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\Bigr{\rangle}\,\mathrm{d}x$
$\displaystyle=-\int\Bigl{\langle}b_{t}\mu_{t}-\frac{1}{2}\mu_{t}\nabla\log\nu_{t},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\Bigr{\rangle}\,\mathrm{d}x\,.$

Finally, we note that

$\displaystyle\int\bigl{[}-\nabla\cdot(b_{t}\mu_{t})+\frac{1}{2}\Delta\mu_{t}\bigr{]}\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\,\mathrm{d}x$
$\displaystyle=\int\Bigl{\langle}b_{t}\mu_{t}-\frac{1}{2}\nabla\mu_{t},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\Bigr{\rangle}\,\mathrm{d}x$
$\displaystyle=\int\Bigl{\langle}b_{t}\mu_{t}-\frac{1}{2}\mu_{t}\nabla\log\mu_{t},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\Bigr{\rangle}\,\mathrm{d}x\,.$

Putting it all together, noticing that the drift terms cancel, we are left
with

$\partial_{t}D_{f}(\mu_{t}\mathbin{\|}\nu_{t})=-\int\frac{1}{2}\Bigl{\langle}\mu_{t}\nabla\log\frac{\mu_{t}}{\nu_{t}},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\Bigr{\rangle}\,\mathrm{d}x=-\frac{1}{2}\mathbb{E}_{\mu_{t}}\Bigl{\langle}\nabla\log\frac{\mu_{t}}{\nu_{t}},\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\Bigr{\rangle}\,,$

which completes the proof. ∎

To recover the decay result for the $q$-Rényi divergence, one can substitute
$g(x)=\frac{1}{q-1}\log x$ and $f(x)=x^{q}-1$. For the $\chi^{2}$-divergence,
instead substitute $g(x)=x$ and $f(x)=x^{2}-1$. From this, we can obtain a
single step of decay for the Rényi and $\chi^{2}$-divergences under different
functional inequalities.

Before proceeding, we need a standard lemma on functional inequalities under
($\mathsf{FH}$).

######  Lemma 20 (Functional inequalities under Gaussian convolutions, [Cha04,
Corollary 13]).

The following inequality holds for any $\pi$ with finite log-Sobolev and
Poincaré constants,

$C_{\mathsf{PI}}(\pi P_{t})\leq
C_{\mathsf{PI}}(\pi)+t\,,\qquad\text{and}\qquad C_{\mathsf{LSI}}(\pi
P_{t})\leq C_{\mathsf{LSI}}(\pi)+t\,.$

Combining the previous two lemmas, we can establish contraction between $\mu
P_{h}Q_{h}$ and $\nu$ after one forward/backward iteration.

######  Theorem 21 ([Che+22, Theorem 3 and 4]).

Assume $\nu$, a measure with smooth density, satisfies (8) with constant
$C_{\mathsf{LSI}}$. For any $q\geq 1$ and initial distribution $\mu$ with a
smooth density, denoting again $Q_{h}:=Q_{h}^{\nu,h}$,

$\mathcal{R}_{q}(\mu
P_{h}Q_{h}\mathbin{\|}\nu)\leq\frac{\mathcal{R}_{q}(\mu\mathbin{\|}\nu)}{(1+h/C_{\mathsf{LSI}})^{2/q}}\,.$

If $\nu$ satisfies (7) with constant $C_{\mathsf{PI}}$, then it follows that

$\chi^{2}(\mu
P_{h}Q_{h}\mathbin{\|}\nu)\leq\frac{\chi^{2}(\mu\mathbin{\|}\nu)}{(1+h/C_{\mathsf{PI}})^{2}}\,.$

Moreover, for all $q\geq 2$,

$\mathcal{R}_{q}(\mu
P_{h}Q_{h}\mathbin{\|}\nu)\leq\begin{cases}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)-\frac{2\log(1+h/C_{\mathsf{PI}})}{q}&\text{if
}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\geq 1\,,\\\
\frac{\mathcal{R}_{q}(\mu\mathbin{\|}\nu)}{(1+h/C_{\mathsf{PI}})^{2}}&\text{if
}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\leq 1\,.\end{cases}$

###### Proof.

Since the SDE in (3.1) captures the forward heat flow ($\mathsf{FH}$), we set
$\mu_{0}$ and $\nu_{0}$ in Lemma 19 to $\mu$ and $\nu$, respectively,
obtaining contraction along the forward heat flow as follows: Substituting the
$q$-Rényi into Lemma 19, we have, from the definition of the Rényi divergence
as
$\mathcal{R}_{q}(\mu\mathbin{\|}\nu):=\frac{1}{q-1}\log\mathopen{}\mathclose{{}\left(D_{f}(\mu\mathbin{\|}\nu)+1}\right),$
with $f(x)=x^{q}-1$ and $g(x)=\frac{1}{q-1}\log(x+1)$,

$\displaystyle\partial_{t}\mathcal{R}_{q}(\mu P_{t}\mathbin{\|}\nu P_{t})$
$\displaystyle=-\frac{q}{2}\,\frac{\mathbb{E}_{\mu
P_{t}}\Bigl{[}\Bigl{\langle}\nabla\Bigl{(}\frac{\mu P_{t}}{\nu
P_{t}}\Bigr{)}^{q-1},\nabla\log\frac{\mu P_{t}}{\nu
P_{t}}\Bigr{\rangle}\Bigr{]}}{(q-1)\,\mathbb{E}_{\nu
P_{t}}\Bigl{[}\Bigl{(}\frac{\mu P_{t}}{\nu
P_{t}}\Bigr{)}^{q}\Bigr{]}}\underset{(i)}{=}-\frac{q}{2}\,\frac{\mathbb{E}_{\mu
P_{t}}\Bigl{[}\Bigl{(}\frac{\mu P_{t}}{\nu
P_{t}}\Bigr{)}^{q-2}\langle\nabla\frac{\mu P_{t}}{\nu
P_{t}},\nabla\log\frac{\mu P_{t}}{\nu P_{t}}\rangle\Bigr{]}}{\mathbb{E}_{\nu
P_{t}}\Bigl{[}\Bigl{(}\frac{\mu P_{t}}{\nu P_{t}}\Bigr{)}^{q}\Bigr{]}}$
$\displaystyle\underset{(ii)}{=}-\frac{q}{2}\,\frac{\mathbb{E}_{\nu
P_{t}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\frac{\mu
P_{t}}{\nu P_{t}}}\right)^{q}\|\nabla\log\frac{\mu P_{t}}{\nu
P_{t}}\|^{2}}\right]}{\mathbb{E}_{\nu
P_{t}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\frac{\mu
P_{t}}{\nu P_{t}}}\right)^{q}}\right]}=-\frac{1}{2}\,\mathsf{RF}_{q}(\mu
P_{t}\mathbin{\|}\nu P_{t})\,,$

where in $(i)$, we use again that
$\nabla\bigl{[}f^{\prime}\bigl{(}\frac{\mu_{t}}{\nu_{t}}\bigr{)}\frac{\mu_{t}}{\nu_{t}}\bigr{]}=\nabla\bigl{(}f^{\prime}\circ\frac{\mu_{t}}{\nu_{t}}\bigr{)}\cdot\frac{\mu_{t}}{\nu_{t}}+f^{\prime}\bigl{(}\frac{\mu_{t}}{\nu_{t}}\bigr{)}\nabla\frac{\mu_{t}}{\nu_{t}}$,
and $(ii)$ uses that $\nabla\frac{\mu P_{t}}{\nu P_{t}}=\frac{\mu P_{t}}{\nu
P_{t}}\nabla\log\frac{\mu P_{t}}{\nu P_{t}}$, and the last equality recalls
the definition of the Rényi Fisher information. This yields

$\partial_{t}\mathcal{R}_{q}(\mu P_{t}\mathbin{\|}\nu
P_{t})=-\frac{1}{2}\mathsf{RF}_{q}(\mu P_{t}\mathbin{\|}\nu
P_{t})\underset{(i)}{\leq}-\frac{1}{q}\frac{\mathcal{R}_{q}(\mu
P_{t}\mathbin{\|}\nu P_{t})}{C_{\mathsf{LSI}}(\nu
P_{t})}\underset{(ii)}{\leq}-\frac{1}{q}\frac{\mathcal{R}_{q}(\mu
P_{t}\mathbin{\|}\nu P_{t})}{C_{\mathsf{LSI}}+t}\,,$

where we used Lemma 17 in $(i)$ and Lemma 20 in $(ii)$. Applying Grönwall’s
inequality (Lemma 10),

$\mathcal{R}_{q}(\mu P_{h}\mathbin{\|}\nu
P_{h})\leq\exp\Bigl{(}-\frac{1}{q}\int_{0}^{h}\frac{1}{C_{\mathsf{LSI}}+t}\,\mathrm{d}t\Bigr{)}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\leq\frac{\mathcal{R}_{q}(\mu\mathbin{\|}\nu)}{(1+h/C_{\mathsf{LSI}})^{1/q}}\,.$

Since the SDE (3.1) also captures the backward equation ($\mathsf{BH}$), we
set $\mu_{0}$ and $\nu_{0}$ in Lemma 19 to $\mu P_{h}$ and $\tilde{\nu}:=\nu
P_{h}$ respectively, obtaining contraction along the backward heat flow:

$\partial_{t}\mathcal{R}_{q}(\mu
P_{h}Q_{t}\mathbin{\|}\tilde{\nu}Q_{t})=-\frac{1}{2}\,\mathsf{RF}_{q}(\mu
P_{h}Q_{t}\mathbin{\|}\tilde{\nu}Q_{t})\leq-\frac{1}{q}\,\frac{\mathcal{R}_{q}(\mu
P_{h}Q_{t}\mathbin{\|}\tilde{\nu}Q_{t})}{C_{\mathsf{LSI}}(\tilde{\nu}Q_{t})}\underset{(i)}{\leq}-\frac{1}{q}\,\frac{\mathcal{R}_{q}(\mu
P_{h}Q_{t}\mathbin{\|}\tilde{\nu}Q_{t})}{C_{\mathsf{LSI}}+{h-t}}\,,$

where $(i)$ follows from that $\tilde{\nu}Q_{t}=\nu P_{h}Q_{t}=\nu P_{h-t}$
and $C_{\mathsf{LSI}}(\tilde{\nu}Q_{t})\leq C_{\mathsf{LSI}}+h-t$ due to Lemma
20. Applying Lemma 10 again yields

$\mathcal{R}_{q}(\mu P_{h}Q_{h}\mathbin{\|}\nu)\leq\frac{\mathcal{R}_{q}(\mu
P_{h}\mathbin{\|}\tilde{\nu})}{(1+h/C_{\mathsf{LSI}})^{1/q}}\,.$

Composing these two inequalities leads to the decay rate claimed in the
theorem.

The result in the $\chi^{2}$-divergence can be derived entirely analogously.
For instance, the decay from the forward part can be shown as follows:

$\partial_{t}\chi^{2}(\mu P_{t}\mathbin{\|}\nu
P_{t})=-\frac{1}{2}\,\mathbb{E}_{\nu P_{t}}\Bigl{[}\Bigl{\|}\nabla\frac{\mu
P_{t}}{\nu P_{t}}\Bigr{\|}^{2}\Bigr{]}\underset{(i)}{\leq}-\frac{\chi^{2}(\mu
P_{t}\mathbin{\|}\nu P_{t})}{C_{\mathsf{PI}}(\nu
P_{t})}\leq-\frac{\chi^{2}(\mu P_{t}\mathbin{\|}\nu
P_{t})}{C_{\mathsf{PI}}+t}\,,$

where $(i)$ follows from Lemma 15. Applying Grönwall’s inequality then gives

$\chi^{2}(\mu P_{h}\mathbin{\|}\nu
P_{h})\leq\exp\Bigl{(}-\int_{0}^{h}\frac{1}{C_{\mathsf{PI}}+t}\,\mathrm{d}t\Bigr{)}\,\chi^{2}(\mu\mathbin{\|}\nu)\leq\frac{\chi^{2}(\mu\mathbin{\|}\nu)}{1+h/C_{\mathsf{PI}}}\,.$

The decay along the backward heat flow in $\chi^{2}$ is entirely analogous to
the Rényi case. Then we combine two contraction results from the forward and
backward flows, completing the proof.

The result in the $\mathcal{R}_{q}$ under (7) can be shown in a similar
manner. Only difference is that in forward and backward computations, one
should use the functional inequality in Lemma 16 and the following standard
inequalities:

$1-\exp\bigl{(}-\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\bigr{)}\geq\begin{cases}\frac{1}{2}&\text{if
}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\geq 1\,,\\\
\frac{1}{2}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)&\text{if
}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\leq 1\,.\end{cases}$

∎

#### Part II: Extension to constrained distributions.

We now prove a technical lemma that extends the contraction results to
constrained distributions. This lemma guarantees the existence of weak
solutions to two stochastic processes that describe the evolution of
distributions involved in Line 2 and 3 in $\mathsf{In\text{-}and\text{-}Out}$,
in addition to lower-semicontinuity of $f$-divergence. We shall prove it for
any measure that is absolutely continuous with respect to $\pi^{X}$, since
this imposes no additional technical hurdles.

######  Lemma 22.

Let $\nu$ be a measure, absolutely continuous with respect to the uniform
measure $\pi^{X}$. The forward and backward heat flow equations given by

$\displaystyle\partial_{t}\mu_{t}$ $\displaystyle=\frac{1}{2}\Delta\mu_{t}\,,$
$\displaystyle\partial_{t}\mu_{t}^{\leftarrow}$
$\displaystyle=-\operatorname{div}\bigl{(}\mu_{t}^{\leftarrow}\nabla\log(\nu
P_{h-t})\bigr{)}+\frac{1}{2}\Delta\mu_{t}^{\leftarrow}\quad\text{with
}\mu_{0}^{\leftarrow}=\mu_{h}\,,$

admit solutions on $(0,h]$, and the weak limit $\lim_{t\to
h}\mu_{t}^{\leftarrow}=\mu_{h}^{\leftarrow}$ exists for any initial measure
$\mu_{0}$ with bounded support. Moreover, for any $f$-divergence with $f$
lower semi-continuous,

$D_{f}(\mu_{h}^{\leftarrow}\mathbin{\|}\nu)\leq\lim_{t\downarrow
0}D_{f}(\mu_{h-t}^{\leftarrow}\mathbin{\|}\nu_{t})\,.$

###### Proof.

The existence of weak solutions for the forward equation is well-known, since
$\mu_{0}$ can be weakly approximated by measures with continuous density, for
which the heat equation admits a unique solution for all time. In particular,
the weak solution is $C^{\infty}$ for $t>0$.

The reverse SDE is more subtle, since $\nabla\log\nu P_{t}$ will in general
cease to be Lipschitz as $t\to 0$. On the other hand, for any $h>0$, we can
write explicitly

$\mu_{h}(x)=\frac{1}{(2\pi
h)^{d/2}}\int\exp\bigl{(}-\frac{\|x-y\|^{2}}{2h}\bigr{)}\,\mathrm{d}\mu_{0}(y)\,.$

If one considers the system started at $\tilde{\mu}_{0}=\mu_{\epsilon}=\nu
P_{\epsilon}$ and solve the forward-backward Fokker-Planck equations on times
$[0,h-\epsilon]$, then
$\tilde{\mu}_{h-\epsilon}=\mu_{h}=\mu_{0}^{\leftarrow}=\tilde{\mu}_{0}^{\leftarrow}$
and

$\mu_{h-\epsilon}^{\leftarrow}(x)=\tilde{\mu}_{h-\epsilon}^{\leftarrow}(x)=\int\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(x)}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(z)\,\mathrm{d}z}\,\mathrm{d}\mu_{h}(y)\,.$

This follows from that if we consider system started at time $\epsilon>0$,
with initial distribution $\mu_{\epsilon}$, then we obtain the above through
the Bayesian perspective on the forward and reverse heat semigroups,
elaborated in Appendix D.

We now show that the following integral is indeed integrable, so
$\tilde{\mu}_{h}^{\leftarrow}$ is well-defined:

$\tilde{\mu}_{h}^{\leftarrow}(x):=\int\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2h}\bigr{)}\nu(x)}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2h}\bigr{)}\nu(z)\,\mathrm{d}z}\,\mathrm{d}\mu_{h}(y)\,.$

For fixed $x$ and $\epsilon<h/2$,

$\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu(z)\,\mathrm{d}z\gtrsim\exp\bigl{(}-\frac{(\|y-x_{0}\|+D)^{2}}{2(h-\epsilon)}\bigr{)}\,,$

as the support of $\nu$ is constrained to $\mathcal{K}\subset B_{D}(x_{0})$.
Since $\mu_{0}$ has bounded support,
$\mu_{h}(y)\lesssim\exp(-\frac{\|y\|^{2}}{a})$ for some constant $a>0$. Thus,

$\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2(h-\epsilon)}\bigr{)}\,\mu_{h}(y)}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(z)\,\mathrm{d}z}\lesssim\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2(h-\epsilon)}\bigr{)}\,\mu_{h}(y)}{\exp(-\frac{(\|y-x_{0}\|+D)^{2}}{2(h-\epsilon)})}\lesssim\exp\Bigl{(}\frac{\langle
2(x-x_{0}),y\rangle+2D\|y-x_{0}\|}{h}-\frac{\|y\|^{2}}{a}\Bigr{)}\,,$

and the last bound is integrable in $y$.

We then show the pointwise convergence of
$\tilde{\mu}_{h-\epsilon}^{\leftarrow}$ to $\tilde{\mu}_{h}^{\leftarrow}$ as
$\epsilon\to 0$. Note that $\nu P_{\epsilon}\to\nu$, as $\nu$ has a bounded
support. Also, the denominator is independent of $\epsilon$ due to

$\frac{1}{(2\pi(h-\epsilon))^{d/2}}\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(z)\,\mathrm{d}z=\mathcal{N}(0,(h-\epsilon)I)*\nu
P_{\epsilon}=\nu*\mathcal{N}(0,hI)\,.$

Hence, for $\epsilon\leq d^{-1}$,

$\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2(h-\epsilon)}\bigr{)}}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(z)\,\mathrm{d}z}\leq\Bigl{(}\frac{h}{h-\epsilon}\Bigr{)}^{d/2}\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2h}\bigr{)}}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2h}\bigr{)}\nu(z)\,\mathrm{d}z}\lesssim\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2h}\bigr{)}}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2h}\bigr{)}\nu(z)\,\mathrm{d}z}\,.$

As shown above, the last bound is integrable with respect to $\mu_{h}$, so the
dominated convergence theorem implies

$\lim_{\epsilon\to
0}\int\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2(h-\epsilon)}\bigr{)}}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2(h-\epsilon)}\bigr{)}\nu
P_{\epsilon}(z)\,\mathrm{d}z}\,\mathrm{d}\mu_{h}(y)=\int\frac{\exp\bigl{(}-\frac{\|x-y\|^{2}}{2h}\bigr{)}}{\int\exp\bigl{(}-\frac{\|z-y\|^{2}}{2h}\bigr{)}\nu(z)\,\mathrm{d}z}\,\mathrm{d}\mu_{h}(y)\,,$

Thus, the pointwise convergence follows. Note that if we take
$\nu(x)=\pi^{X}(x)=\frac{\mathds{1}_{\mathcal{K}}(x)}{\operatorname{vol}(\mathcal{K})}$,
then $\tilde{\mu}_{h}^{\leftarrow}$ is the distribution of the backwards step
of our algorithm. In particular, this corresponds to first sampling
$x\sim\mu_{h}$, then $y\sim Q_{h}^{\nu,h}(\cdot|x)$, which is precisely the
law of $\mu_{h}^{\leftarrow}$ given by ($\mathsf{FP\text{-}BH}$).

As for the second statement, it follows from Scheffé’s lemma [Bil95, Theorem
16.12] that the pointwise convergence of
$\mu_{h-\varepsilon}^{\leftarrow}\to\mu_{h}^{\leftarrow}$ leads to its
$\mathsf{TV}$-convergence, which in turn implies the weak convergence. It
follows from lower semicontinuity of $D_{f}$ [AFP00, Theorem 2.34] that the
weak convergence ensures
$D_{f}(\mu_{h}^{\leftarrow}\mathbin{\|}\nu)\leq\lim_{t\downarrow
0}D_{f}(\mu_{h-t}^{\leftarrow}\mathbin{\|}\nu^{\leftarrow}_{h-t})$. ∎

In the sequel, we will only consider $\nu=\pi^{X}$. Since the Rényi divergence
is a continuous function of the $\chi^{q}$ divergence (see Definition 6),
which itself is an $f$-divergence, it enjoys the same lower-semicontinuity
properties. Using this lower-semicontinuity together with the decay results in
Theorem 21, we can easily derive the contraction results of
$\mathsf{In\text{-}and\text{-}Out}$ in $\mathcal{R}_{q}$ and $\chi^{q}$ for
any $q\geq 1$. We remark that this result does not require convexity of
$\mathcal{K}$.

######  Theorem 23.

Let $\mu^{X}_{k}$ be the law of the $k$-th output of
$\mathsf{In\text{-}and\text{-}Out}$ with initial distribution $\mu_{0}^{X}$.
Let $C_{\mathsf{LSI}}$ be the (8) constant of the uniform distribution
$\pi^{X}$ over $\mathcal{K}$. Then, for any $q\geq 1$,

$\mathcal{R}_{q}(\mu_{k}^{X}\mathbin{\|}\pi^{X})\leq\frac{\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})}{(1+h/C_{\mathsf{LSI}})^{2k/q}}\,.$

For $C_{\mathsf{PI}}$ the (7) constant of $\pi^{X}$,

$\chi^{2}(\mu^{X}_{k}\mathbin{\|}\pi^{X})\leq\frac{\chi^{2}(\mu^{X}_{0}\mathbin{\|}\pi^{X})}{(1+h/C_{\mathsf{PI}})^{2k}}\,.$

Furthermore, for any $q\geq 2$,

$\mathcal{R}_{q}(\mu_{k}^{X}\mathbin{\|}\pi^{X})\leq\begin{cases}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-\frac{2k\log(1+h/C_{\mathsf{PI}})}{q}&\text{if
}k\leq\frac{q}{2\log(1+h/C_{\mathsf{PI}})}\,\bigl{(}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-1\bigr{)}\,,\\\
(1+h/C_{\mathsf{PI}})^{-2(k-k_{0})/q}&\text{if }k\geq
k_{0}:=\lceil\frac{q}{2\log(1+h/C_{\mathsf{PI}})}\,\bigl{(}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})-1\bigr{)}\rceil\,.\end{cases}$

###### Proof.

Let us set $\mu_{0}=\mu_{0}^{X}$ and $\pi_{0}=\pi^{X}$. Then,
$\mu_{h}=\mu_{0}^{\leftarrow}=\mu_{1}^{Y}$,
$\pi_{h}=\pi_{0}^{\leftarrow}=\pi^{Y}$, and
$\mu_{h}^{\leftarrow}=\mu_{1}^{X}$, $\pi_{h}^{\leftarrow}=\pi^{X}$. For small
$\epsilon>0$, as
$\mu_{\epsilon}=(\mu_{0}^{X})_{\epsilon}=\mu_{0}^{X}*\mathcal{N}(0,\epsilon
I_{d})$ is $C^{\infty}$-smooth, we can now invoke the decay results with step
size $h-\epsilon$ in Theorem 21. Thus, for contraction constants
$C_{\epsilon}=(1+\frac{h-\epsilon}{C_{\mathsf{LSI}}+\epsilon})^{-2/q}$ and
$C_{\epsilon}=(1+\frac{h-\epsilon}{C_{\mathsf{PI}}+\epsilon})^{-2}$
respectively when $\Phi=\mathcal{R}_{q}$ and $\Phi=\chi^{2}$,

$\Phi(\mu_{h-\epsilon}^{\leftarrow}\mathbin{\|}\pi_{\epsilon})\leq
C_{\epsilon}\cdot\Phi(\mu_{\epsilon}\mathbin{\|}\pi_{\epsilon})\leq
C_{\epsilon}\cdot\Phi(\mu_{0}\mathbin{\|}\pi_{0})\,,$

where we used the data-processing inequality for the last inequality. By the
second result of Lemma 22, sending $\epsilon\to 0$ leads to

$\Phi(\mu_{1}^{X}\mathbin{\|}\pi^{X})=\Phi(\mu_{h}^{\leftarrow}\mathbin{\|}\pi_{0})\leq
C\cdot\Phi(\mu_{0}\mathbin{\|}\pi_{0})=C\cdot\Phi(\mu_{0}^{X}\mathbin{\|}\pi^{X})\,.$

Repeating this argument $k$ times completes the proof. ∎

###  3.3 Failure probability and wasted steps

We begin by defining a suitable version of _local conductance_ [KLS97].

######  Definition 24 (Local conductance).

The local conductance $\ell$ on $\mathbb{R}^{d}$ is defined by

$\ell(x)\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\frac{\int_{\mathcal{K}}\exp(-\frac{1}{2h}\|x-y\|^{2})\,\mathrm{d}y}{\int_{\mathbb{R}^{d}}\exp(-\frac{1}{2h}\|x-y\|^{2})\,\mathrm{d}y}=\frac{\int_{\mathcal{K}}\exp(-\frac{1}{2h}\|x-y\|^{2})\,\mathrm{d}y}{(2\pi
h)^{d/2}}\,.$

The local conductance at $y$ quantifies the success probability of the
proposal at $y$ in Line 3. Then the expected number of trials until the first
success of Line 3 is $1/\ell(y)$. Revisiting (1.1), we can notice
$\pi^{Y}(y)=\ell(y)/\operatorname{vol}(\mathcal{K})$.

#### Naïve analysis for expected number of trials.

Starting from $\pi^{X}$, when we just naïvely sample from $\pi^{Y|X}(\cdot|x)$
for all $x$ without imposing any _failure_ condition, the expected number of
trials for one iteration is that for the probability density $p_{x}$ of
$\mathcal{N}(x,hI_{d})$,

$\int_{\mathcal{K}}\int_{\mathbb{R}^{d}}\frac{1}{\ell(y)}\,p_{x}(\mathrm{d}y)\pi^{X}(\mathrm{d}x)=\int_{\mathbb{R}^{d}}\frac{1}{\ell(y)}\,\pi^{Y}(\mathrm{d}y)=\int_{\mathbb{R}^{d}}\frac{1}{\ell(y)}\,\frac{\ell(y)}{\operatorname{vol}(\mathcal{K})}\,\mathrm{d}y=\infty\,.$

This suggests that one should consider the algorithm as having “failed” if the
number of trials exceeds some threshold.

#### Refined analysis under a failure condition.

Going forward, we assume an $M$-warm start as in previous work for uniform
sampling algorithms. By induction we have
$\frac{\mathrm{d}\mu^{X}_{i}}{\mathrm{d}\pi^{X}}\leq M$ for all $i$.

######  Lemma 25 (Propagation of warm-start).

From an $M$-warm start, we have
$\nicefrac{{\mathrm{d}\mu_{i}^{X}}}{{\mathrm{d}\pi^{X}}}\leq M$ for all $i$.

###### Proof.

Assume that $\mu_{i}^{X}$ satisfies the $M$-warm start. Then, for any
measurable $S$ and the transition kernel $T_{x}$ of Algorithm 1 at $x$,

$\mu_{i+1}^{X}(S)=\int_{\mathcal{K}}T_{x}(S)\,\mathrm{d}\mu_{i}^{X}(x)\leq
M\int_{\mathcal{K}}T_{x}(S)\,\mathrm{d}\pi^{X}(x)=M\pi^{X}(S)\,,$

where the last equality follows from the stationarity of $\pi$. Hence,
$\mathrm{d}\mu_{i+1}^{X}/\mathrm{d}\pi^{X}\leq M$. ∎

We now establish a lemma that comes in handy when analyzing the failure
probability of the algorithm. In essence, this lemma bounds the probability
that taking a Gaussian step from $\pi^{X}$ in Line 2 gets $\delta$-distance
away from $\mathcal{K}$. Let us denote the $\delta$-blowup of $\mathcal{K}$ by
$\mathcal{K}_{\delta}:=\\{x\in\mathbb{R}^{d}:d(x,\mathcal{K})\leq\delta\\}$.

######  Lemma 26.

$\pi^{Y}(\mathcal{K}_{\delta}^{c})\leq
e^{c_{1}^{2}}\exp(-\nicefrac{{t^{2}}}{{8c_{1}^{2}}})$ for $\delta=t/d$ and
variance $h=c_{1}^{2}/d^{2}$, where $c_{1}$ is some constant and $t\geq
2c_{1}(c_{1}+1)$.

###### Proof.

For $y\in\partial\mathcal{K}_{\delta}$, we can take the supporting half-space
$H$ at $\mathsf{proj}_{\mathcal{K}}(y)$ containing $\mathcal{K}$, due to
convexity of $\mathcal{K}$. Then,

$\displaystyle\pi^{Y}(\mathcal{K}_{\delta}^{c})$
$\displaystyle=\frac{1}{\operatorname{vol}(\mathcal{K})}\int_{\mathcal{K}_{\delta}^{c}}\int_{\mathcal{K}}\frac{\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}}{(2\pi
h)^{d/2}}\,\mathrm{d}x\,\mathrm{d}y\leq\frac{1}{\operatorname{vol}(\mathcal{K})}\int_{\mathcal{K}_{\delta}^{c}}\int_{H}\frac{\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}}{(2\pi
h)^{d/2}}\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=\frac{1}{\operatorname{vol}(\mathcal{K})}\int_{\mathcal{K}_{\delta}^{c}}\int_{d(y,\mathcal{K})}^{\infty}\frac{\exp\bigl{(}-\frac{1}{2h}z^{2}\bigr{)}}{\sqrt{2\pi
h}}\,\mathrm{d}z\,\mathrm{d}y\,.$ (3.3)

Let us denote the tail probability of the $1$-dimensional Gaussian with
variance $h$ by

$\mathsf{T}(s):=\mathbb{P}_{Z\sim\mathcal{N}(0,h)}(Z\geq
s)=\frac{1}{\sqrt{2\pi
h}}\int_{s}^{\infty}\exp\bigl{(}-\frac{1}{2h}z^{2}\bigr{)}\,\mathrm{d}z\,.$

By the co-area formula and integration by parts,

$\displaystyle\int_{\mathcal{K}_{\delta}^{c}}\int_{d(y,\mathcal{K})}^{\infty}\frac{\exp\bigl{(}-\frac{1}{2h}z^{2}\bigr{)}}{\sqrt{2\pi
h}}\,\mathrm{d}z\mathrm{d}y=\int_{\delta}^{\infty}\mathsf{T}(s)\operatorname{vol}(\partial\mathcal{K}_{s})\,\mathrm{d}s$
$\displaystyle=\Bigl{[}\underbrace{\mathsf{T}(s)\int_{0}^{s}\operatorname{vol}(\partial\mathcal{K}_{z})\,\mathrm{d}z}_{\eqqcolon\mathsf{F}}\Bigr{]}_{\delta}^{\infty}+\int_{\delta}^{\infty}\frac{1}{\sqrt{2\pi
h}}\exp\bigl{(}-\frac{s^{2}}{2h}\bigr{)}\int_{0}^{s}\operatorname{vol}(\partial\mathcal{K}_{z})\,\mathrm{d}z\,\mathrm{d}s\,.$
(3.4)

Recall that $\mathsf{T}(s)\leq\frac{1}{2}\exp(-\frac{1}{2}(sh^{-1/2})^{2})$
for $s\geq\delta=t/d$ due to a standard tail bound on a Gaussian distribution.
This tail bound, combined with

$\int_{0}^{s}\operatorname{vol}(\partial\mathcal{K}_{z})\,\mathrm{d}z=\operatorname{vol}(\mathcal{K}_{s})-\operatorname{vol}(\mathcal{K})\leq\operatorname{vol}\bigl{(}(1+s)\,\mathcal{K}\bigr{)}-\operatorname{vol}(\mathcal{K})=\bigl{(}(1+s)^{d}-1\bigr{)}\operatorname{vol}(\mathcal{K})\,,$

ensures that $\mathsf{F}$ vanishes at $s=\infty$. Hence, bounding the first
term in (3.4) by $0$ results in

$\displaystyle\int_{\mathcal{K}_{\delta}^{c}}\int_{d(y,\mathcal{K})}^{\infty}\frac{\exp\bigl{(}-\frac{1}{2h}z^{2}\bigr{)}}{\sqrt{2\pi
h}}\,\mathrm{d}z\,\mathrm{d}y$ $\displaystyle\leq\frac{1}{\sqrt{2\pi
h}}\int_{\delta}^{\infty}\exp\bigl{(}-\frac{s^{2}}{2h}\bigr{)}\bigl{(}\underbrace{(1+s)^{d}}_{\leq\exp(sd)}-1\bigr{)}\operatorname{vol}(\mathcal{K})\,\mathrm{d}s$
$\displaystyle\leq\frac{\operatorname{vol}(\mathcal{K})}{\sqrt{2\pi
h}}\exp(hd^{2}/2)\int_{\delta}^{\infty}\exp\Bigl{(}-\frac{1}{2h}\,(s-hd)^{2}\Bigr{)}\,\mathrm{d}s$
$\displaystyle\underset{(i)}{\leq}\operatorname{vol}(\mathcal{K})\exp(hd^{2}/2)\exp\bigl{(}-\frac{1}{2}\,(\delta
h^{-1/2}-dh^{1/2})^{2}\bigr{)}$
$\displaystyle\underset{(ii)}{\leq}\operatorname{vol}(\mathcal{K})\exp(c_{1}^{2}/2)\exp\bigl{(}-\frac{t^{2}}{8c_{1}^{2}}\bigr{)}\,,$

where in $(i)$ we used the tail bound of a Gaussian, and $(ii)$ follows from
that in the regime of $h=c_{1}^{2}d^{-2}$, we have $\delta
h^{-1/2}-dh^{1/2}-1\geq t/c_{1}-(c_{1}+1)\geq\nicefrac{{t}}{{2c_{1}}}$ for
$t\geq 2c_{1}(c_{1}+1)$. Putting the last bound into (3.3) completes the proof
of the claim. ∎

Now we choose a suitable threshold $N$ for bounding the failure probability.
Following (3.3) in the proof, one can notice that for
$y\in\mathcal{K}^{c}_{\delta}$, $\delta=\Omega(1/d)$, and $h=\Theta(d^{-2})$,

$\ell(y)\leq\int_{d(y,\mathcal{K})}^{\infty}\frac{\exp\bigl{(}-\frac{1}{2h}z^{2}\bigr{)}}{\sqrt{2\pi
h}}\,\mathrm{d}z=\mathbb{P}_{Z\sim\mathcal{N}(0,h)}(Z\geq\delta)\leq\exp(-\Omega(t^{2}))\,.$

Thus, the expected number of trials from $\mathcal{K}^{c}_{\delta}$ for the
rejection sampling in Line 3 is $\ell(y)^{-1}\geq\exp(\Omega(t^{2}))$.
Intuitively, one can ignore whatever happens in $\mathcal{K}^{c}_{\delta}$,
since $\mathcal{K}_{\delta}$ takes up most of measure of $\pi^{Y}$. As the
number of trials from $\mathcal{K}_{\delta}^{c}$ is at least
$\exp(\Omega(t^{2}))$ in expectation, the most straightforward way to ignore
algorithmic behaviors from $\mathcal{K}_{\delta}^{c}$ is simply to set the
threshold to $N=\widetilde{\mathcal{O}}(\exp(t^{2}))$. Even though the
threshold is $N$, the expected number of trials is much lower.

Lemma 14 bounds the failure probability and expected number of trials per
iteration.

###### Proof of Lemma 14.

For $\mu_{h}:=\mu*\mathcal{N}(0,hI_{d})$, the failure probability is
$\mathbb{E}_{\mu_{h}}[(1-\ell)^{N}]$. Since
$\mathrm{d}\mu/\mathrm{d}\pi^{X}\leq M$ implies
$\mathrm{d}\mu_{h}/\mathrm{d}(\pi^{X})_{h}=\mathrm{d}\mu_{h}/\mathrm{d}\pi^{Y}\leq
M$, it follows that

$\mathbb{E}_{\mu_{h}}[(1-\ell)^{N}]\leq
M\,\mathbb{E}_{\pi^{Y}}[(1-\ell)^{N}]\,.$

Then,

$\displaystyle\int_{\mathbb{R}^{d}}\underbrace{(1-\ell)^{N}\,\mathrm{d}\pi^{Y}}_{\eqqcolon\mathsf{A}}$
$\displaystyle=\int_{\mathcal{K}_{\delta}^{c}}\mathsf{A}+\int_{\mathcal{K}_{\delta}\cap[\ell\geq
N^{-1}\log(3mM/\eta)]}\mathsf{A}+\int_{\mathcal{K}_{\delta}\cap[\ell<N^{-1}\log(3mM/\eta)]}\mathsf{A}$
$\displaystyle\leq\pi^{Y}(\mathcal{K}_{\delta}^{c})+\int_{[\ell\geq
N^{-1}\log(3mM/\eta)]}\exp(-\ell
N)\,\mathrm{d}\pi^{Y}+\int_{\mathcal{K}_{\delta}\cap[\ell<N^{-1}\log(3mM/\eta)]}\frac{\ell(y)}{\operatorname{vol}(\mathcal{K})}\,\mathrm{d}y$
$\displaystyle\leq
e^{c_{1}^{2}}\exp(-\frac{t^{2}}{8c_{1}^{2}})+\frac{\eta}{3mM}+\frac{\log(3mM/\eta)}{N}\,\frac{\operatorname{vol}(\mathcal{K}_{\delta})}{\operatorname{vol}(\mathcal{K})}$
$\displaystyle\leq
e^{c_{1}^{2}}\exp(-\frac{t^{2}}{8c_{1}^{2}})+\frac{\eta}{3mM}+\frac{e^{t}}{N}\,\log\frac{3mM}{\eta}\,,$

where we used
$\operatorname{vol}(\mathcal{K}_{\delta})\subset\operatorname{vol}\bigl{(}(1+\delta)\mathcal{K}\bigr{)}=(1+\delta)^{d}\operatorname{vol}(\mathcal{K})\leq
e^{t}\operatorname{vol}(\mathcal{K})$. Taking $c_{1}^{2}=\frac{\log\log
Z}{2\log Z}$, $t=\sqrt{8}\log\log Z$, and $N=Z(\log Z)^{4}$, we can bound the
last line by $\frac{\eta}{mM}$. Therefore,

$\mathbb{E}_{\mu_{h}}[(1-\ell(\cdot))^{N}]\leq
M\,\mathbb{E}_{\pi^{Y}}[(1-\ell(\cdot))^{N}]\leq\frac{\eta}{m}\,.$

We now bound the expected number of trials per iteration. Let $X$ be the
minimum of the threshold $N$ and the number of trials until the first success.
Then the expected number of trials per step is bounded by
$M\mathbb{E}_{\pi^{Y}}[X]$ due to $\mathrm{d}\mu_{h}/\mathrm{d}\pi^{Y}\leq M$.
Thus,

$\displaystyle\int_{\mathbb{R}^{d}}\Bigl{(}\frac{1}{\ell}\wedge
N\Bigr{)}\,\mathrm{d}\pi^{Y}$
$\displaystyle\leq\int_{\mathcal{K}_{\delta}}\frac{1}{\ell}\,\mathrm{d}\pi^{Y}+N\pi^{Y}(\mathcal{K}_{\delta}^{c})=\frac{\operatorname{vol}(\mathcal{K}_{\delta})}{\operatorname{vol}(\mathcal{K})}+N\pi^{Y}(\mathcal{K}_{\delta}^{c})$
$\displaystyle\leq e^{t}+N\exp(-\Omega(t^{2}))\leq(\log Z)^{3}+3(\log
Z)^{4}=\mathcal{O}\Bigl{(}\bigl{(}\log\frac{mM}{\eta}\bigr{)}^{4}\Bigr{)}\,.$

Therefore, the expected number of trials per step is
$\mathcal{O}\bigl{(}M(\log\frac{mM}{\eta})^{4}\bigr{)}$, and the claim follows
since each trial uses one query to the membership oracle of $\mathcal{K}$. ∎

###  3.4 Putting it together

We can now show that $\mathsf{In\text{-}and\text{-}Out}$ subsumes previous
results on uniform sampling from convex bodies (such as $\mathsf{Ball\ walk}$
and $\mathsf{Speedy\ walk}$), providing detailed versions of the main results
in Section 1.2.

We first establish that the query complexity of
$\mathsf{In\text{-}and\text{-}Out}$ matches that of the $\mathsf{Ball\ walk}$
under stronger divergences. Recall that
$2\|\cdot\|_{\mathsf{TV}}^{2}\leq\mathsf{KL}\leq\log(1+\chi^{2})\leq\chi^{2}$.

######  Theorem 27.

For any given $\eta,\varepsilon\in(0,1)$, $q\geq 1$, $m\in\mathbb{N}$ defined
below and any convex body $\mathcal{K}$ given by a well-defined membership
oracle, consider $\mathsf{In\text{-}and\text{-}Out}$ (Algorithm 1) with an
$M$-warm initial distribution $\mu_{0}^{X}$,
$h=(2d^{2}\log\frac{9mM}{\eta})^{-1}$, and
$N=\widetilde{\mathcal{O}}(\frac{mM}{\eta})$. For $\pi^{X}$ the uniform
distribution over $\mathcal{K}$,

  * •

It achieves $\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$
after
$m=\widetilde{\mathcal{O}}(qd^{2}\|\mathrm{Cov}(\pi^{X})\|_{\mathsf{op}}\log\frac{M}{\eta\varepsilon})$
iterations. With probability $1-\eta$, the algorithm iterates this many times
without failure, using
$\widetilde{\mathcal{O}}(qMd^{2}\|\mathrm{Cov}(\pi^{X})\|_{\mathsf{op}}(\log\frac{1}{\eta\varepsilon})^{5})$
expected number of membership queries in total.

  * •

For isotropic $\pi^{X}$, with probability $1-\eta$, the algorithm achieves
$\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$ with
$m=\widetilde{\mathcal{O}}(qd^{2}\log\frac{M}{\eta\varepsilon})$ iterations,
using $\widetilde{\mathcal{O}}(qMd^{2}(\log\frac{1}{\eta\varepsilon})^{5})$
membership queries in expectation.

###### Proof.

We just put together Lemma 14 and Theorem 23. For target accuracy
$\varepsilon>0$, we use the $\mathcal{R}_{q}$-decay under (7) for $q\geq 2$ in
Theorem 23. The $M$-warm start assumption guarantees
$\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})\lesssim\log M$. Due to
$C_{\mathsf{PI}}(\pi^{X})=\mathcal{O}(\|\mathrm{Cov}(\pi^{X})\|_{\mathsf{op}}\,\log
d)$ (Lemma 18), $\mathsf{In\text{-}and\text{-}Out}$ can achieve
$\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$ after
$m=\widetilde{\mathcal{O}}(qd^{2}\|\mathrm{Cov}(\pi^{X})\|_{\mathsf{op}}\log\frac{M}{\eta\varepsilon})$
iterations. Since each iteration has $\eta/m$-failure probability by Lemma 14,
the union bound ensures that the total failure probability is at most $\eta$
throughout $m$ iterations. Lastly, each iteration requires
$\widetilde{\mathcal{O}}(M(\log\frac{1}{\eta\varepsilon})^{4})$ membership
queries in expectation by Lemma 14. Therefore,
$\mathsf{In\text{-}and\text{-}Out}$ uses
$\widetilde{\mathcal{O}}(qMd^{2}\min(D^{2},\|\mathrm{Cov}(\pi^{X})\|_{\mathsf{op}})(\log\frac{1}{\eta\varepsilon})^{5})$
expected number of membership queries over $m$ iterations. Since
$\mathcal{R}_{q}$ is non-decreasing in $q$, we can obtain the desired bound on
$\mathcal{R}_{q}$ for $q\in[1,2)$.

For isotropic $\pi^{X}$, we have $\mathrm{Cov}(\pi^{X})=I_{d}$, so the claim
immediately follows from $C_{\mathsf{PI}}(\pi^{X})=\mathcal{O}(\log d)$ (see
Lemma 18). ∎

We now show that the number of proper steps is bounded as claimed for general
_non-convex bodies_ and _any feasible start_ in $\mathcal{K}$. We first
establish this result under an $M$-warm start (Theorem 5).

###### Proof of Theorem 5..

By the Rényi-decay under (8) in Theorem 23,
$\mathsf{In\text{-}and\text{-}Out}$ can achieve $\varepsilon$-distance to
$\pi^{X}$ after
$\mathcal{O}\bigl{(}qh^{-1}C_{\mathsf{LSI}}(\pi^{X})\log\frac{\mathcal{R}_{q}(\mu^{X}_{1}\mathbin{\|}\pi^{X})}{\varepsilon}\bigr{)}$
iterations for $q\geq 1$.

For $q\geq 2$, we use the decay result under (7). In this case,
$\mathsf{In\text{-}and\text{-}Out}$ decays under two different rates depending
on the value of $\mathcal{R}_{q}(\cdot\mathbin{\|}\pi^{X})$. It first needs
$\mathcal{O}(qh^{-1}C_{\mathsf{PI}}(\pi^{X})\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X}))$
iterations until $\mathcal{R}_{q}(\cdot\mathbin{\|}\pi^{X})$ reaches $1$.
Then, $\mathsf{In\text{-}and\text{-}Out}$ additionally needs
$\mathcal{O}(qh^{-1}C_{\mathsf{PI}}(\pi^{X})\log\frac{1}{\varepsilon})$
iterations, and thus it needs
$\mathcal{O}(qh^{-1}C_{\mathsf{PI}}(\pi^{X})\bigl{(}\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})+\log\frac{1}{\varepsilon}\bigr{)})$
iterations in total. By substituting
$\mathcal{R}_{q}(\mu_{0}^{X}\mathbin{\|}\pi^{X})\lesssim\log M$, we complete
the proof. ∎

Next, we show that $\mathsf{In\text{-}and\text{-}Out}$ mixes from any start.

######  Corollary 28.

For any given $\varepsilon\in(0,1)$ and set $\mathcal{K}\subset B_{D}(0)$,
$\mathsf{In\text{-}and\text{-}Out}$ with variance $h$ and any feasible start
$x_{0}\in\mathcal{K}$ achieves
$\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$ after
$m=\widetilde{\mathcal{O}}(qh^{-1}C_{\mathsf{LSI}}(\pi^{X})\log\frac{d+D^{2}/h}{\varepsilon})$
iterations.

###### Proof.

We first bound the warmness of $\mu^{X}_{1}$ w.r.t. $\pi^{X}$ when
$\mu^{X}_{0}=\delta_{x_{0}}$. One can readily check that

$\mu_{1}^{X}(x)=\mathds{1}_{\mathcal{K}}(x)\cdot\int\frac{\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2h}\|y-x_{0}\|^{2}\bigr{)}}{(2\pi
h)^{d/2}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\,\mathrm{d}x}\,\mathrm{d}y\,.$

By Young’s inequality,
$\|y-x\|^{2}\leq(\|y\|+D)^{2}\leq\frac{3}{2}\|y\|^{2}+3D^{2}$ for
$x\in\mathcal{K}$. Hence,

$\displaystyle\int\frac{\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2h}\|y-x_{0}\|^{2}\bigr{)}}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\,\mathrm{d}x}\,\mathrm{d}y$
$\displaystyle\leq$
$\displaystyle\frac{\exp(2h^{-1}D^{2})}{\operatorname{vol}(\mathcal{K})}\int\exp\Bigl{(}-\frac{1}{2h}\bigl{(}\|y-x\|^{2}+\|y-x_{0}\|^{2}-\frac{3}{2}\|y\|^{2}\bigr{)}\Bigr{)}\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\frac{\exp(2h^{-1}D^{2})}{\operatorname{vol}(\mathcal{K})}\int\exp\Bigl{(}-\frac{1}{2h}\bigl{(}\frac{1}{2}\|y-2(x+x_{0})\|^{2}+(\|x\|^{2}+\|x_{0}\|^{2}-2\|x+x_{0}\|^{2})\bigr{)}\Bigr{)}\,\mathrm{d}y$
$\displaystyle\leq$
$\displaystyle\frac{\exp(5h^{-1}D^{2})}{\operatorname{vol}(\mathcal{K})}\int\exp\Bigl{(}-\frac{1}{4h}\|y-2(x+x_{0})\|^{2}\Bigr{)}\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\frac{\exp(5h^{-1}D^{2})}{\operatorname{vol}(\mathcal{K})}(4\pi
h)^{d/2}\,.$

Therefore, $M=\operatorname{ess\,sup}\frac{\mu_{1}^{X}}{\pi^{X}}\leq
2^{d/2}\exp(5h^{-1}D^{2})$. By Theorem 5 under (8),
$\mathsf{In\text{-}and\text{-}Out}$ needs
$\widetilde{\mathcal{O}}(qh^{-1}C_{\mathsf{LSI}}(\pi^{X})\log\frac{d+D^{2}/h}{\varepsilon})$
iterations. ∎

We then obtain the following corollary for a convex body $\mathcal{K}$.

######  Corollary 29.

For any given $\varepsilon\in(0,1)$ and convex body $\mathcal{K}\subset
B_{D}(0)$, $\mathsf{In\text{-}and\text{-}Out}$ with variance $h$ and an
$M$-warm initial distribution achieves
$\mathcal{R}_{q}(\mu^{X}_{m}\mathbin{\|}\pi^{X})\leq\varepsilon$ after
$m=\widetilde{\mathcal{O}}(qh^{-1}D^{2}\log\frac{1}{\varepsilon})$ iterations.
If $\pi^{X}$ is isotropic, then $\mathsf{In\text{-}and\text{-}Out}$ only needs
$\widetilde{\mathcal{O}}(qh^{-1}D\log\frac{d+d^{2}/h}{\varepsilon})$
iterations.

###### Proof.

For convex $\mathcal{K}$, it follows from Lemma 18 that
$C_{\mathsf{LSI}}(\pi^{X})=\mathcal{O}(D^{2})$ and
$C_{\mathsf{LSI}}(\pi^{X})=\mathcal{O}(D)$ for isotropic $\mathcal{K}$. The
rest of the proof can be completed in a similar way. ∎

For $h=\tilde{\Theta}(d^{-2})$, $\mathsf{In\text{-}and\text{-}Out}$ requires
$\widetilde{\mathcal{O}}(qd^{2}D^{2})$ iterations and in particular
$\widetilde{\mathcal{O}}(qd^{2}D)$ iteration for isotropic uniform
distributions. These results match those of $\mathsf{Speedy\ walk}$ [KLM06,
LV17] (see Theorem 31).

#### Acknowledgements.

We are deeply grateful to Andre Wibisono and Sinho Chewi for helpful comments
and pointers to the literature for Lemma 19. This work was supported in part
by NSF award 210644, NSERC through the CGS-D award, and a Simons Investigator
award.

range

pages14 rangepages10 rangepages39 rangepages-1 rangepages24 rangepages27
rangepages39 rangepages1 rangepages31 rangepages28 rangepages28 rangepages37
rangepages12 rangepages21 rangepages46 rangepages17 rangepages38 rangepages49
rangepages92 rangepages41 rangepages3 rangepages14 rangepages11 rangepages30
rangepages19 rangepages50 rangepages20 rangepages26 rangepages13 rangepages66
rangepages9 rangepages58 rangepages7 rangepages10 rangepages17 rangepages15
rangepages25 rangepages16 rangepages19 rangepages9 rangepages54 rangepages21
rangepages10 rangepages43 rangepages13 rangepages13 rangepages7 rangepages1
rangepages27 rangepages28

## References

  * [AC21] Kwangjun Ahn and Sinho Chewi  “Efficient constrained sampling via the mirror-Langevin algorithm”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 34, 2021, pp. 28405–28418
  * [AFP00] Luigi Ambrosio, Nicola Fusco and Diego Pallara  “Functions of bounded variation and free discontinuity problems”  Oxford University Press, 2000
  * [AGS05] Luigi Ambrosio, Nicola Gigli and Giuseppe Savaré  “Gradient flows: in metric spaces and in the space of probability measures”  Springer Science & Business Media, 2005
  * [BNN06] Mikhail Belkin, Hariharan Narayanan and Partha Niyogi  “Heat flow and a faster algorithm to compute the surface area of a convex body”  In _Foundations of Computer Science (FOCS)_ , 2006, pp. 47–56  IEEE
  * [Bes+95] Julian Besag, Peter Green, David Higdon and Kerrie Mengersen  “Bayesian computation and stochastic systems”  In _Statistical Science_ JSTOR, 1995, pp. 3–41
  * [Bil95] Patrick Billingsley  “Probability and measure”  John Wiley & Sons, 1995
  * [Bin+19] Eli Bingham et al.  “Pyro: Deep universal probabilistic programming”  In _The Journal of Machine Learning Research (JMLR)_ 20, 2019, pp. 28:1–28:6
  * [Bro+17] Nicolas Brosse, Alain Durmus, Ãric Moulines and Marcelo Pereyra  “Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo”  In _Conference on Learning Theory (COLT)_ 65 PMLR, 2017, pp. 319–342
  * [BEL15] Sebastien Bubeck, Ronen Eldan and Joseph Lehec  “Finite-time analysis of projected Langevin Monte Carlo”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 28, 2015
  * [BEL18] Sébastien Bubeck, Ronen Eldan and Joseph Lehec  “Sampling from a log-concave distribution with projected Langevin Monte Carlo”  In _Discrete & Computational Geometry (DCG)_ 59 Springer, 2018, pp. 757–783
  * [Cha04] Djalil Chafaï  “Entropies, convexity, and functional inequalities, On Phi-entropies and Phi-Sobolev inequalities”  In _Journal of Mathematics of Kyoto University_ 44.2 Duke University Press, 2004, pp. 325–363
  * [Che70] Jeff Cheeger  “A lower bound for the smallest eigenvalue of the Laplacian”  In _Problems in analysis_ 625.195-199 Princeton UP, 1970, pp. 110
  * [Che+22] Yongxin Chen, Sinho Chewi, Adil Salim and Andre Wibisono  “Improved analysis for a proximal algorithm for sampling”  In _Conference on Learning Theory (COLT)_ , 2022, pp. 2984–3014  PMLR
  * [Che21] Yuansi Chen  “An almost constant lower bound of the isoperimetric coefficient in the KLS conjecture”  In _Geometric and Functional Analysis (GAFA)_ 31 Springer, 2021, pp. 34–61
  * [Che24] Sinho Chewi  “Log-concave sampling”  Book draft available at https://chewisinho.github.io, 2024
  * [CV16] Ben Cousins and Santosh S Vempala  “A practical volume algorithm”  In _Mathematical Programming Computation_ 8.2 Springer, 2016, pp. 133–160
  * [CV18] Ben Cousins and Santosh S Vempala  “Gaussian Cooling and $\mathcal{O}^{*}(n^{3})$ algorithms for volume and Gaussian volume”  In _SIAM Journal on Computing (SICOMP)_ 47.3 SIAM, 2018, pp. 1237–1273
  * [Dal17] Arnak Dalalyan  “Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent”  In _Conference on Learning Theory (COLT)_ , 2017, pp. 678–689  PMLR
  * [DT12] Arnak S Dalalyan and Alexandre B Tsybakov  “Sparse regression learning by aggregation and Langevin Monte-Carlo”  In _Journal of Computer and System Sciences_ 78.5 Elsevier, 2012, pp. 1423–1443
  * [DMM19] Alain Durmus, Szymon Majewski and Blazej Miasojedow  “Analysis of Langevin Monte Carlo via convex optimization”  In _The Journal of Machine Learning Research (JMLR)_ 20.1 JMLR. org, 2019, pp. 2666–2711
  * [DFK91] Martin Dyer, Alan Frieze and Ravi Kannan  “A random polynomial-time algorithm for approximating the volume of convex bodies”  In _Journal of the ACM (JACM)_ 38.1 ACM New York, NY, USA, 1991, pp. 1–17
  * [Eld13] Ronen Eldan  “Thin shell implies spectral gap up to polylog via a stochastic localization scheme”  In _Geometric and Functional Analysis (GAFA)_ 23.2 Springer, 2013, pp. 532–569
  * [FYC23] Jiaojiao Fan, Bo Yuan and Yongxin Chen  “Improved dimension dependence of a proximal algorithm for sampling”  In _Conference on Learning Theory (COLT)_ , 2023, pp. 1473–1521  PMLR
  * [GKV23] Khashayar Gatmiry, Jonathan Kelner and Santosh S Vempala  “Sampling with barriers: Faster mixing via Lewis weights”  In _arXiv preprint arXiv:2303.00480_ , 2023
  * [GC11] Mark Girolami and Ben Calderhead  “Riemann manifold Langevin and Hamiltonian Monte Carlo methods”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73.2 Wiley Online Library, 2011, pp. 123–214
  * [Gop+23] Sivakanth Gopi et al.  “Algorithmic aspects of the Log-Laplace transform and a non-Euclidean proximal sampler”  In _Conference on Learning Theory (COLT)_ 195 PMLR, 2023, pp. 2399–2439
  * [GLS88] Martin Grötschel, László Lovász and Alexander Schrijver  “Geometric algorithms and combinatorial optimization”  Springer, 1988
  * [GHZ22] Mert Gürbüzbalaban, Yuanhan Hu and Lingjiong Zhu  “Penalized Langevin and Hamiltonian Monte Carlo algorithms for constrained sampling”  In _arXiv preprint arXiv:2212.00570_ , 2022
  * [Har+17] Hulda S Haraldsdóttir et al.  “CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models”  In _Bioinformatics_ 33.11 Oxford University Press, 2017, pp. 1741–1743
  * [Jia+21] He Jia, Aditi Laddha, Yin Tat Lee and Santosh S Vempala  “Reducing isotropy and volume to KLS: an $\mathcal{O}^{*}(n^{3}\psi^{2})$ volume algorithm”  In _Symposium on Theory of Computing (STOC)_ , 2021, pp. 961–974
  * [Jia21] Qijia Jiang  “Mirror Langevin Monte Carlo: the case under isoperimetry”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 34, 2021, pp. 715–725
  * [KLM06] Ravi Kannan, László Lovász and Ravi Montenegro  “Blocking conductance and mixing in random walks”  In _Combinatorics, Probability and Computing_ 15.4 Cambridge University Press, 2006, pp. 541–570
  * [KLS95] Ravi Kannan, László Lovász and Miklós Simonovits  “Isoperimetric problems for convex bodies and a localization lemma”  In _Discrete & Computational Geometry_ 13.3 Springer, 1995, pp. 541–559
  * [KLS97] Ravi Kannan, László Lovász and Miklós Simonovits  “Random walks and an $\mathcal{O}^{*}(n^{5})$ volume algorithm for convex bodies”  In _Random Structures & Algorithms (RS&A)_ 11.1 Wiley Online Library, 1997, pp. 1–50
  * [KN12] Ravindran Kannan and Hariharan Narayanan  “Random walks on polytopes and an affine interior point method for linear programming”  In _Mathematics of Operations Research_ 37.1 INFORMS, 2012, pp. 1–20
  * [KP21] Bo’az Klartag and Eli Putterman  “Spectral monotonicity under Gaussian convolution”  In _arXiv preprint arXiv:2107.09496_ , 2021
  * [Kla23] Boáz Klartag  “Logarithmic bounds for isoperimetry and slices of convex sets”  In _Ars Inveniendi Analytica_ Ars Inveniendi, 2023
  * [KL22] Boáz Klartag and Joseph Lehec  “Bourgain’s slicing problem and KLS isoperimetry up to polylog”  In _Geometric and Functional Analysis (GAFA)_ 32.5 Springer, 2022, pp. 1134–1159
  * [Koo+22] Yunbum Kook, Yin Tat Lee, Ruoqi Shen and Santosh S Vempala  “Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 35, 2022, pp. 31684–31696
  * [Koo+23] Yunbum Kook, Yin Tat Lee, Ruoqi Shen and Santosh S Vempala  “Condition-number-independent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators”  In _Conference on Learning Theory (COLT)_ 195 PMLR, 2023, pp. 4504–4569
  * [KV23] Yunbum Kook and Santosh S Vempala  “Gaussian Cooling and Dikin walks: The Interior-Point Method for logconcave sampling”  In _arXiv preprint arXiv:2307.12943_ , 2023
  * [Led94] Michel Ledoux  “A simple analytic proof of an inequality by P. Buser”  In _Proceedings of the American Mathematical Society_ 121.3, 1994, pp. 951–959
  * [LST21] Yin Tat Lee, Ruoqi Shen and Kevin Tian  “Structured logconcave sampling with a restricted Gaussian oracle”  In _Conference on Learning Theory (COLT)_ , 2021, pp. 2993–3050  PMLR
  * [LV18] Yin Tat Lee and Santosh S Vempala  “Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation”  In _Symposium on Theory of Computing (STOC)_ , 2018, pp. 1115–1121
  * [LV17] Yin Tat Lee and Santosh Srinivas Vempala  “Eldan’s stochastic localization and the KLS hyperplane conjecture: An improved lower bound for expansion”  In _Foundations of Computer Science (FOCS)_ , 2017, pp. 998–1007  IEEE
  * [Leh23] Joseph Lehec  “The Langevin Monte Carlo algorithm in the non-smooth log-concave case”  In _The Annals of Applied Probability_ 33.6A Institute of Mathematical Statistics, 2023, pp. 4858–4874
  * [LNP12] Nathan E Lewis, Harish Nagarajan and Bernhard O Palsson  “Constraining the metabolic genotype–phenotype relationship using a phylogeny of in silico methods”  In _Nature Reviews Microbiology_ 10.4 Nature Publishing Group, 2012, pp. 291–305
  * [Li+22] Ruilin Li, Molei Tao, Santosh S Vempala and Andre Wibisono  “The mirror Langevin algorithm converges with vanishing bias”  In _International Conference on Algorithmic Learning Theory (ALT)_ , 2022, pp. 718–742  PMLR
  * [Liu20] Yuan Liu  “The Poincaré inequality and quadratic transportation-variance inequalities”  In _Electronic Journal of Probability_ 25.1, 2020, pp. 1–16
  * [Lov99] László Lovász  “Hit-and-run mixes fast”  In _Mathematical Programming_ 86 Springer, 1999, pp. 443–461
  * [LS90] László Lovász and Miklós Simonovits  “The mixing rate of Markov chains, an isoperimetric inequality, and computing the volume”  In _Foundations of Computer Science (FOCS)_ , 1990, pp. 346–354  IEEE
  * [LS93] László Lovász and Miklós Simonovits  “Random walks in a convex body and an improved volume algorithm”  In _Random Structures & Algorithms (RS&A)_ 4.4 Wiley Online Library, 1993, pp. 359–412
  * [LV03] László Lovász and Santosh S Vempala  “Hit-and-run is fast and fun”  In _preprint, Microsoft Research_ , 2003
  * [LV06] László Lovász and Santosh S Vempala  “Hit-and-run from a corner”  In _SIAM Journal on Computing (SICOMP)_ 35.4 SIAM, 2006, pp. 985–1005
  * [MT07] Frank McSherry and Kunal Talwar  “Mechanism design via differential privacy”  In _Foundations of Computer Science (FOCS)_ , 2007, pp. 94–103  IEEE
  * [Mil09] Emanuel Milman  “On the role of convexity in isoperimetry, spectral gap and concentration”  In _Inventiones Mathematicae_ 177.1 Springer, 2009, pp. 1–43
  * [Mir17] Ilya Mironov  “Rényi differential privacy”  In _Computer Security Foundations Symposium (CSF)_ , 2017, pp. 263–275  IEEE
  * [Smi84] Robert L Smith  “Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions”  In _Operations Research_ 32.6 INFORMS, 1984, pp. 1296–1308
  * [SWW23] Vishwak Srinivasan, Andre Wibisono and Ashia Wilson  “Fast sampling from constrained spaces using the Metropolis-adjusted mirror Langevin algorithm”  In _arXiv preprint arXiv:2312.08823_ , 2023
  * [Sta20] Stan Development Team  “RStan: the R interface to Stan” R package version 2.21.2, 2020  URL: http://mc-stan.org/
  * [Thi+13] Ines Thiele et al.  “A community-driven global reconstruction of human metabolism”  In _Nature Biotechnology_ 31.5 Nature Publishing Group, 2013, pp. 419–425
  * [Vem05] Santosh S Vempala  “Geometric random walks: a survey”  In _Combinatorial and Computational Geometry_ 52.573-612 Cambridge University Press Cambridge, 2005, pp. 2
  * [VW19] Santosh S Vempala and Andre Wibisono  “Rapid convergence of the unadjusted Langevin algorithm: Isoperimetry suffices”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 32, 2019
  * [Vil09] Cédric Villani  “Optimal transport: old and new”  Springer, 2009
  * [Yua+23] Bo Yuan et al.  “On a class of Gibbs sampling over networks”  In _Conference on Learning Theory (COLT)_ , 2023, pp. 5754–5780  PMLR
  * [Zha+20] Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili and Marcelo Pereyra  “Wasserstein control of mirror Langevin Monte Carlo”  In _Conference on Learning Theory (COLT)_ , 2020, pp. 3814–3841  PMLR

##  Appendix A Related work

Sampling from constrained log-concave distributions is a fundamental task
arising in many fields. Uniform sampling with convex constraints is its
simplest manifestation, which was first studied as a core subroutine for a
randomized volume-computation algorithm [DFK91]. Since then, this fundamental
problem has been studied for over three decades [LS90, LS93, KLS97, LV03,
LV06, BEL18, Bro+17]. We review these algorithms, grouping them under three
categories — geometric random walks, structured samplers, and diffusion-type
samplers. Below, $\mathcal{K}$ is convex.

#### Geometric random walk.

We discuss two geometric random walks – $\mathsf{Ball\ walk}$ [LS93, KLS97]
and $\mathsf{Hit\text{-}and\text{-}Run}$ [Smi84, Lov99].

$\mathsf{Ball\ walk}$ is a simple metropolized random walk; it draws $y$
uniformly at random from a ball of radius $\delta$ centered at a current point
$x$, and moves to $y$ if $y\in\mathcal{K}$ and stays at $x$ otherwise. In the
literature, $\mathsf{Ball\ walk}$ actually refers to a composite algorithm
consisting of [$\mathsf{Speedy\ walk}+$ rejection sampling], where
$\mathsf{Speedy\ walk}$ records only the accepted steps of $\mathsf{Ball\
walk}$ (see Appendix B for details). The step size $\delta$ should be set to
$\mathcal{O}(d^{-1/2})$ to avoid stepping outside of $\mathcal{K}$. [KLS97]
showed that $\mathsf{Ball\ walk}$ needs
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log\frac{1}{\varepsilon})$ membership
queries to be $\varepsilon$-close to $\pi^{\mathcal{K}}$ in $\mathsf{TV}$,
where $D$ is the diameter of $\mathcal{K}$, and the warmness parameter $M$
measures the closeness of the initial distribution to the target uniform
distribution $\pi^{\mathcal{K}}$.

$\mathsf{Hit\text{-}and\text{-}Run}$ is another zeroth-order algorithm that
needs _no step size_ ; it picks a uniform random line $\ell$ passing a current
point, and move to a uniform random point on $\ell\cap\mathcal{K}$. [LV06]
shows that, if we define the second moment as
$R^{2}:=\mathbb{E}_{X\sim\pi^{\mathcal{K}}}[\|X-\mathbb{E}[X]\|^{2}]$, then
$\mathsf{Hit\text{-}and\text{-}Run}$ requires
$\mathcal{O}(d^{2}R^{2}\log\frac{M}{\varepsilon})$ queries. Notably, this
algorithm has a poly-logarithmic dependence on $M$ as opposed to
$\mathsf{Ball\ walk}$.

Both algorithm are affected by skewed shape of $\mathcal{K}$ (i.e., large $D$
or $R$), so these samplers are combined with pre-processing step called
_rounding_. This procedure finds a linear transformation that makes the
geometry of $\mathcal{K}$ less skewed and so more amenable to sampling. In
literature, there exists a randomized algorithm [Jia+21] that rounds
$\mathcal{K}$ and generates a good warm start (i.e., $M=\mathcal{O}(1)$), with
$\mathsf{Ball\ walk}$ used as a core subroutine. This algorithm takes up
$\widetilde{\mathcal{O}}(d^{3})$ queries in total, and in such position with
the good warm start, $\mathsf{Ball\ walk}$ only needs
$\widetilde{\mathcal{O}}(d^{2}\log\frac{1}{\varepsilon})$ queries to sample
from $\pi^{\mathcal{K}}$.

#### Structured samplers.

The aforementioned samplers based on geometric random walks require _only_
access to the membership oracle of the convex body _without_ any additional
structural assumptions. The alternate paradigm of _geometry-aware sampling_
attempts to exploit the _structure_ of convex constraints, with the aim of
expediting the convergence of the resultant sampling schemes. One common
assumption is to make available a _self-concordant barrier function_ $\phi$
which has regularity on its high-order derivatives and blows up when
approaching the boundary $\partial\mathcal{K}$. The Hessian of $\phi$ encodes
the local geometry of the constraint, and the samplers often work directly
with $\nabla^{2}\phi$.

The first canonical example of such a zeroth-order sampler is $\mathsf{Dikin\
walk}$ used when $\mathcal{K}$ is given by $m$ linear constraints [KN12]; it
draws a uniform sample from an ellipsoid (characterized by $\nabla^{2}\phi$)
of fixed radius around a current point, and is often combined with a
Metropolis adjustment. [KN12] shows that $\mathsf{Dikin\ walk}$ mixes in
$\mathcal{O}(md\log\frac{M}{\varepsilon})$ steps, although each iteration is
slightly more expensive than one membership query. This algorithm requires no
rounding, but still needs a good warm-start, which can be achieved by an
annealing-type algorithm using $\widetilde{\mathcal{O}}(md)$ iterations of
$\mathsf{Dikin\ walk}$ [KV23].

$\mathsf{Riemannian\ Hamiltonian\ Monte\ Carlo}$ is a structured sampler that
exploits the first-order information of the potential (i.e.,
$\nabla\log(1/\pi)$) [GC11]; its proposal is given as the solution to the
Hamilton’s ODE equation, followed by the Metropolis-filter. In the linear-
constraint setting above, this sampler requires
$\mathcal{O}(md^{2/3}\log\frac{M}{\varepsilon})$ many iterations to achieve
$\varepsilon$-close distance to $\pi^{\mathcal{K}}$ [LV18]. This sampler is
further analyzed for practical ODE solvers [Koo+23] and for more sophisticated
self-concordant barriers [GKV23].

Similarly, $\mathsf{Mirror\ Langevin}$ [Zha+20, Jia21, AC21, Li+22] is a class
of algorithms which converts the constrained problem into an unconstrained one
obtained by considering the pushforward of the constrained space by
$\nabla\phi$. The algorithm can also be metropolized [SWW23]. The best known
rates for this algorithm are
$\widetilde{\mathcal{O}}(d\log\frac{1}{\varepsilon})$ under some strong
assumptions on $\phi$.

#### Diffusion-type samplers.

Samplers based on discretizations of Itô diffusions, stochastic processes
which rapidly mix to $\pi$ in continuous time, have long been used for
sampling without constraints [Bes+95, DT12, Dal17, Che24]. While the
underlying stochastic processes generalize easily to constrained settings, the
discretization analysis relies crucially on the smoothness of the target
distribution. This is clearly impossible to achieve in the constrained
setting, and so some techniques are required to circumvent this difficulty.
These algorithms, however, generalize easily to the more general problem of
sampling from distributions of the form $\tilde{\pi}^{X}\propto
e^{-f}\mathds{1}_{\mathcal{K}}$, by naturally incorporating first order
information from $f$.

The first approach for adapting diffusion-based samplers [BEL15, BEL18, Leh23]
iterates a two-step procedure. First, a random step is taken, with
$x_{k+1/2}\sim\mathcal{N}(x_{k},2hI_{d})$ for some appropriately chosen step
$h$,444A gradient step can be added in the more general case, for sampling
from $\tilde{\pi}^{X}$ and then project it to $\mathcal{K}$, i.e.,
$x_{k+1}=\mathsf{proj}_{\mathcal{K}}(x_{k+1/2})$. The complexity is given in
terms of queries to a _projection oracle_ , each call to which can be
implemented with a polynomial number of membership oracle queries; a total of
$\tilde{\mathcal{O}}(\frac{d^{2}D^{3}}{\varepsilon^{4}})$ queries are needed
to be $\varepsilon$-close in $\mathcal{W}_{2}$ to $\pi^{X}$. Another approach,
which uses an algorithmically designed “soft” penalty instead of a projection,
was proposed in [GHZ22], and achieves a rate estimate of
$\tilde{\mathcal{O}}(d/\varepsilon^{10})$.

A second approach, suggested by [Bro+17], considers a different proximal
scheme, which performs a “soft projection” onto $\mathcal{K}$, by taking steps
like
$\mathcal{N}((1-h\lambda^{-1})x_{k}+h\mathsf{proj}_{\mathcal{K}}(x_{k}),2hI_{d})$.
It is called Moreau-Yosida regularized Langevin, named after an analogous
regularization scheme for constrained optimization. This scheme also relies on
access to a projection oracle for $\mathcal{K}$, and quantifies their query
complexity accordingly. Their final rate estimate is
$\tilde{\mathcal{O}}(\frac{d^{5}}{\varepsilon^{6}})$ to be $\varepsilon$-close
in $\mathsf{TV}$ distance to $\pi^{X}$.

Observing the prior work integrating diffusion-based sampling with convex
constraints, the dependence on the key parameters $d,\varepsilon$, while
polynomial, are many orders worse than the rates for zeroth-order samplers
such as $\mathsf{Ball\ walk},\mathsf{Hit\text{-}and\text{-}Run}$. In contrast,
our analysis not only recovers but in some sense surpasses the known rates for
$\mathsf{Ball\ walk},\mathsf{Hit\text{-}and\text{-}Run}$, while harmonizing
well with the continuous-time perspective of diffusions.

#### Proximal schemes for sampling.

The Gibbs sampling scheme used in this paper was inspired by the restricted
Gaussian oracle introduced in [LST21] (in turn inspired by Gaussian Cooling
[CV18]), which alternately iterates between a pure Gaussian step, and a
“proximal” step (which we elaborate in our exposition). This scheme was given
novel interpretations by [Che+22], which showed that it interpolates the
forward and backward heat flows, in the sense defined by [KP21]. The backward
heat flow itself is intimately related to stochastic localization schemes,
invented and popularized in [Eld13, Che21].

This formulation proved surprisingly powerful, allowing many existing rates in
unconstrained sampling to be recovered from a relatively simple analysis. This
was further extended by [FYC23] to achieve the current state-of-the-art rate
in unconstrained sampling. Finally, [Gop+23] suggest that this could be
applied to tackle some constrained problems. However, the assumptions in this
final mentioned work are not compatible with the uniform sampling problem on
general convex bodies.

##  Appendix B $\mathsf{Ball\ walk}$ and $\mathsf{Speedy\ walk}$

We restate the previously known guarantees for uniform sampling by
$\mathsf{Ball\ walk}$ and $\mathsf{Speedy\ walk}$. Below, let $B_{r}(x)$
denote the $d$-dimensional ball of radius $r$ centered at $x$.

Input: initial distribution $\pi_{0}$, convex body
$\mathcal{K}\subset\mathbb{R}^{d}$, iterations $T$, step size $\delta>0$.

1: Sample $x_{0}\sim\pi_{0}$.

2: for $i=1,\dotsc,T$ do

3: Sample $y\sim\mathcal{\textnormal{Unif}}(B_{\delta}(x_{i-1}))$.

4: If $y\in\mathcal{K}$, then $x_{i}\leftarrow y$. Else, $x_{i}\leftarrow
x_{i-1}$.

5: end for

Algorithm 2 $\mathsf{Ball\ walk}$

$\mathsf{Ball\ walk}$ is particularly simple; draw a uniform random point from
$B_{\delta}$ around the current point, and go there if the drawn point is
inside of $\mathcal{K}$ and stay at the current point otherwise. Its
stationary distribution can be easily seen to be
$\pi\propto\mathds{1}_{\mathcal{K}}$, the uniform distribution over
$\mathcal{K}$.

In the literature, there are two approaches to analyzing the convergence rate
of this sampler: (i) a direct analysis via the $s$-conductance of
$\mathsf{Ball\ walk}$ and (ii) an indirect approach which first passes through
$\mathsf{Speedy\ walk}$.

#### Direct analysis.

The following $\mathsf{TV}$-guarantee is obtained by lower bounding the
$s$-conductance of $\mathsf{Ball\ walk}$, which requires a one-step coupling
argument and the Cheeger inequality for $\pi$. We refer interested readers to
[Vem05, Section 5].

######  Theorem 30 (Convergence of $\mathsf{Ball\ walk}$).

For any $\varepsilon\in(0,1)$ and convex body
$\mathcal{K}\subset\mathbb{R}^{d}$ presented by a well-defined membership
oracle, let $\pi_{t}$ be the distribution after $t$ steps of $\mathsf{Ball\
walk}$ with an $M$-warm initial distribution $\pi_{0}$. Then, $\mathsf{Ball\
walk}$ with step size $\delta=\Theta(\frac{\varepsilon}{M\sqrt{d}})$ achieves
$\|\pi_{t}-\pi\|_{\mathsf{TV}}\leq\varepsilon$ for $t\gtrsim
d^{2}D^{2}\frac{M^{2}}{\varepsilon^{2}}\log\frac{M}{\varepsilon}$. If $\pi$ is
isotropic, then $\mathsf{Ball\ walk}$ only needs $\mathcal{O}(d^{2}\log
d\frac{M^{2}}{\varepsilon^{2}}\log\frac{M}{\varepsilon})$ iterations.

The mixing time of $\mathsf{Ball\ walk}$ under this approach has a polynomial
dependence on $1/\varepsilon$, rather than a polylogarithmic dependence.

#### Indirect analysis through $\mathsf{Speedy\ walk}$.

[KLS97] introduced $\mathsf{Speedy\ walk}$, which could be viewed as a version
of $\mathsf{Ball\ walk}$ and converges to a _speedy distribution_ (see
Proposition 1), which is slightly biased from $\pi$. Then, $\mathsf{Speedy\
walk}$ is used together with another algorithmic component (rejection
sampling) [KLS97, Algorithm 4.15] that converts the speedy distribution to the
uniform distribution. In the literature, $\mathsf{Ball\ walk}$ often refers to
‘$\mathsf{Speedy\ walk}$ combined with the conversion step’, rather than a
direct implementation of Algorithm 2. Strictly speaking, a mixing guarantee of
this combined algorithm should not be referred to as a provable guarantee of
$\mathsf{Ball\ walk}$.

Input: initial distribution $\pi_{0}$, convex body
$\mathcal{K}\subset\mathbb{R}^{d}$, iterations $T$, step size $\delta>0$.

1: Sample $x_{0}\sim\pi_{0}$.

2: for $i=1,\dotsc,T$ do

3: Sample
$x_{i}\sim\textnormal{Unif}\mathopen{}\mathclose{{}\left(\mathcal{K}\cap
B_{\delta}(x_{i-1})}\right)$.

4: end for

Algorithm 3 $\mathsf{Speedy\ walk}$

As opposed to $\mathsf{Ball\ walk}$, $\mathsf{Speedy\ walk}$ _always_ takes
some step at each iteration. However, the problem of sampling from
$x_{i}\sim\textnormal{Unif}\mathopen{}\mathclose{{}\left(\mathcal{K}\cap
B_{\delta}(x_{i-1})}\right)$ in Line 3 is not straightforward. This step
admits the following implementation based on rejection sampling, via a
procedure denoted by $(\ast)$:

  * •

Propose $y\sim\textnormal{Unif}(B_{\delta}(x_{i-1}))$.

  * •

Set $x_{i+1}\leftarrow y$ if $y\in\mathcal{K}$. Otherwise, repeat the
proposal.

Each actual step (indexed by $i$) in $\mathsf{Speedy\ walk}$ is called a
_proper step_ , and rejected steps during $(\ast)$ are called _improper
steps_. For example, if $x_{1},x_{1},x_{2},x_{3},x_{3},x_{3},x_{4},\dotsc$ are
the positions produced by $\mathsf{Ball\ walk}$, then only proper steps
$x_{1},x_{2},x_{3},x_{4},\dotsc$ are recorded by $\mathsf{Speedy\ walk}$.

To describe the theoretical guarantees of $\mathsf{Speedy\ walk}$, we define
the _local conductance_ $\ell(x)$ at $x\in\mathcal{K}$, which measures the
success probability of the rejection sampling scheme in $(\ast)$:

$\ell(x):=\frac{\operatorname{vol}(\mathcal{K}\cap
B_{\delta}(x))}{\operatorname{vol}(B_{\delta}(x))}\,,$

and define the _average conductance_ :

$\lambda:=\mathbb{E}_{\pi}[\ell]=\frac{1}{\operatorname{vol}(\mathcal{K})}\int_{\mathcal{K}}\ell(x)\,\mathrm{d}x\,.$

######  Proposition 1 ([KLS97]).

The stationary distribution $\nu$ of $\mathsf{Speedy\ walk}$ has density

$\nu(x)=\frac{\ell(x)\,\mathds{1}_{\mathcal{K}}(x)}{\int_{\mathcal{K}}\ell(x)\,\mathrm{d}x}\,.$

The speedy distribution $\nu$ is indeed different from the uniform
distribution $\pi$, and this discrepancy is quantified in terms of the average
conductance.

######  Proposition 2 ([KLS97, Page 22]).

$\|\nu-\pi\|_{\mathsf{TV}}\leq\frac{1-\lambda}{\lambda}$.

One can relate the step size $\delta$ to the average conductance.

######  Proposition 3 (Bound on average conductance, [KLS97, Corollary 4.5]).

$\lambda\geq 1-\frac{\delta\sqrt{d}}{2}$.

The best known result for $\mathsf{Speedy\ walk}$’s mixing is due to [KLM06]
devising the _blocking conductance_ and using the _log-Cheeger_ inequality.
When $\nu$ is isotropic (i.e., it has covariance proportional to the identity
matrix), [LV17] improves the mixing bound via the _log-Cheeger_ constant.

######  Theorem 31 (Mixing of $\mathsf{Speedy\ walk}$).

For any $\varepsilon\in(0,1)$ and convex body
$\mathcal{K}\subset\mathbb{R}^{d}$ presented by a well-defined membership
oracle, let $\nu_{t}$ be the distribution after $t$ proper steps of
$\mathsf{Speedy\ walk}$ started at any feasible point $x_{0}\in\mathcal{K}$.
Then, $\mathsf{Speedy\ walk}$ with step size $\delta=\Theta(d^{-1/2})$
achieves $\|\nu_{t}-\nu\|_{\mathsf{TV}}\leq\varepsilon$ for
$t\gtrsim(D^{2}+\log(D\sqrt{d}))\,d^{2}\log\frac{1}{\varepsilon}$. From an
$M$-warm start, the expected number of improper steps during $t$ iterations is
$\widetilde{\mathcal{O}}(tM)$. When $\nu$ is isotropic, $\mathsf{Speedy\
walk}$ needs $\mathcal{O}(d^{2}D\log\frac{1}{\varepsilon}\log\log D)$ proper
steps to achieve $\varepsilon$-$\mathsf{TV}$ distance to $\nu$.

Then, [KLS97] uses the following post-processing step to obtain an
approximately uniform distribution on $\mathcal{K}$, with a provable
guarantee.

$\mathcal{A}$: Call $\mathsf{Speedy\ walk}$ to obtain a sample $X\sim\nu_{t}$
until $\frac{2d}{2d-1}\,X\in\mathcal{K}$. If so, return
$\bar{X}=\frac{2d}{2d-1}\,X$.

######  Proposition 4 ([KLS97, Theorem 4.16]).

Under the same setting above, assume
$\|\nu_{t}-\nu\|_{\mathsf{TV}}\leq\varepsilon$ for step size
$\delta\leq(8d\log\frac{d}{\varepsilon})^{-1/2}$ and fixed $t\in\mathbb{N}$.
For $\bar{\nu}=\mathsf{law}(\bar{X})$ given by $\mathcal{A}$, it holds that
$\|\bar{\nu}-\pi\|_{\mathsf{TV}}\leq\varepsilon$, and the expected number of
calls on the conversion algorithm is at most $2$.

Combining the previous two results, we conclude that the total expected number
of membership queries to obtain a sample $\varepsilon$-close to $\pi$ in
$\mathsf{TV}$ is
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log\frac{1}{\varepsilon})$, which now has
a poly-logarithmic dependence on $1/\varepsilon$.

######  Remark 2 (Backward heat flow analysis of $\mathsf{Speedy\ walk}$).

Consider a Gaussian version of $\mathsf{Speedy\ walk}$, whose one-step
corresponds to $x_{i+1}\sim\mathcal{N}(x_{i},hI_{d})|_{\mathcal{K}}$, and this
transition kernel exactly matches integrating ($\mathsf{BH}$) for time $h$.
Thus, $\nu Q^{\pi^{X},h}_{h}=\nu$ due to the stationarity of $\nu$ under
$\mathsf{Speedy\ walk}$, where $Q^{\pi^{X},h}_{h}$ is the transition kernel
defined by the backward heat flow for time $h$ that reverses
$\pi^{X}*\mathcal{N}(0,hI_{d})$ to $\pi^{X}$. Hence, if we can control the
LSI/PI constants of $\nu$ along the backward heat-flow’s trajectory, then we
could directly analyze $\mathsf{Speedy\ walk}$ by emulating computations in
Lemma 21.

##  Appendix C Functional inequalities

We provide full details on functional inequalities omitted in Section 3.1. We
use $\mu$ and $\mu_{\mathsf{LC}}$ to denote a probability measure and log-
concave probability measure over $\mathbb{R}^{d}$, respectively.

#### Cheeger and PI constants.

The _Cheeger isoperimetric constant_ $C_{\mathsf{Ch}}(\mu)$ measures how large
surface area a measurable subset with larger volume has, defined by

$C_{\mathsf{Ch}}(\mu):=\inf_{S\subset\mathbb{R}^{d}}\frac{\mu^{+}(S)}{\min(\mu(S),\mu(S^{c}))}\,,$

where the infimum is taken over all measurable subsets $S$, and $\mu^{+}(S)$
is the Minkowski content of $S$ under $\mu$ defined as, for
$S^{\varepsilon}:=\\{x\in X:d(x,S)<\varepsilon\\}$,

$\mu^{+}(S):=\liminf_{\varepsilon\to
0}\frac{\mu(S^{\varepsilon})-\mu(S)}{\varepsilon}\,.$

[Che70] established $C_{\mathsf{PI}}(\mu)\lesssim
C_{\mathsf{Ch}}^{-2}(\mu)$555The opposite direction
$C_{\mathsf{PI}}(\mu_{\mathsf{LC}})\gtrsim
C_{\mathsf{Ch}}^{-2}(\mu_{\mathsf{LC}})$ also holds for log-concave
distributions due to [Mil09], while $C_{\mathsf{PI}}(\mu)\gtrsim
C_{\mathsf{Ch}}^{-2}(\mu)/d$ for general distributions due to [Led94]., and
then [KLS95] showed that for covariance matrix
$\Sigma_{\mu}:=\mathbb{E}_{\mu}[(\cdot-\mathbb{E}_{\mu}X)(\cdot-\mathbb{E}_{\mu}X)^{\mathsf{T}}]$,

$C_{\mathsf{Ch}}(\mu_{\mathsf{LC}})\gtrsim\frac{1}{(\mathbb{E}_{\mu_{\mathsf{LC}}}[\|X-\mathbb{E}_{\mu_{\mathsf{LC}}}X\|^{2}])^{1/2}}=\frac{1}{(\operatorname{tr}(\Sigma_{\mu_{\mathsf{LC}}}))^{1/2}}\,.$
(C.1)

This immediately leads to
$C_{\mathsf{PI}}(\pi)\lesssim(\mathbb{E}_{\pi}[\|X-\mathbb{E}_{\pi}X\|^{2}])^{1/2}\leq
D^{2}$ for the uniform distribution $\pi$ over a convex body $\mathcal{K}$
with diameter $D>0$.

Kannan et al. proposed the _KLS conjecture_ in the same paper, which says that
for the spectral norm $\|\cdot\|_{2}$,

$C_{\mathsf{Ch}}(\mu_{\mathsf{LS}})\gtrsim\frac{1}{\|\Sigma_{\mu_{\mathsf{LS}}}\|_{2}^{1/2}}\,.$

While the original result in [KLS95] ensures $C_{\mathsf{Ch}}\gtrsim d^{-1/2}$
for an isotropic log-concave distribution (due to $\Sigma=I_{d}$), this
conjecture indeed claims $C_{\mathsf{Ch}}\gtrsim 1$ for such case. Following a
line of work [LV17, Che21, KL22, Kla23], the current bound is

$C_{\mathsf{Ch}}(\mu_{\mathsf{LS}})\gtrsim\frac{(\log
d)^{-1/2}}{\|\Sigma_{\mu_{\mathsf{LS}}}\|_{2}^{1/2}}\,,$

which implies that $C_{\mathsf{PI}}(\pi)\lesssim\log d$ when $\pi$ is
isotropic for convex $\mathcal{K}$.

#### Log-Cheeger and LSI constants.

Just as the Cheeger and PI constants are related above, there are known
connections between LSI and _log-Cheeger_ constants. The log-Cheeger constant
$C_{\mathsf{logCh}}(\mu)$ of a distribution
$\mu\in\mathcal{P}(\mathbb{R}^{d})$ is defined as

$C_{\mathsf{logCh}}(\mu):=\inf_{S\subset\mathbb{R}^{d}:\mu(S)\leq\frac{1}{2}}\frac{\mu^{+}(S)}{\mu(S)\sqrt{\log\frac{1}{\mu(S)}}}\,.$

[Led94] established that $C_{\mathsf{LSI}}(\mu)\lesssim
C_{\mathsf{logCh}}^{-2}(\mu)$666The opposite direction holds under dimension-
scaling due to [Led94]: $C_{\mathsf{LSI}}(\mu)\gtrsim
C_{\mathsf{logCh}}^{-2}(\mu)/d$., and [KLM06] showed that any log-concave
distributions with support of diameter $D>0$ satisfy
$C_{\mathsf{logCh}}(\mu_{\mathsf{LS}})\gtrsim D^{-1}$. Later in 2016, [LV17]
improved this to $C_{\mathsf{logCh}}(\mu_{\mathsf{LS}})\gtrsim D^{-1/2}$ under
isotropy. Therefore, for convex $\mathcal{K}$, it follows that
$C_{\mathsf{LSI}}(\pi)\lesssim D^{2}$ and that $C_{\mathsf{LSI}}(\pi)\lesssim
D$ if $\pi$ is isotropic.

##  Appendix D The Wasserstein geometry

We present additional technical background on the Wasserstein geometry and
Markov semigroup theory. Interested readers can refer to [Vil09, AGS05, Che24]
for standard references on Wasserstein spaces and applications to sampling.

#### Wasserstein gradient.

Let $\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ be the space of probability
measures admitting densities on $\mathbb{R}^{d}$ with finite second moment.
Although there are many ways to metrize $\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})$, the geometry induced by the Wasserstein-$2$ distance
$\mathcal{W}_{2}$ is a particularly useful structure for analysis.

Under the $\mathcal{W}_{2}$-geometry, one can define a “gradient” of a
functional defined over $\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$.
Specifically, for a functional $\mathcal{F}:\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})\to\mathbb{R}\cup\\{\infty\\}$, the _Wasserstein gradient_
of $\mathcal{F}$ at $\mu\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ is defined
as $\nabla_{\mathcal{W}_{2}}\mathcal{F}(\mu)=\nabla(\delta\mathcal{F})(\mu)\in
L^{2}(\mu)$, where $\nabla$ is the standard gradient and $\delta\mathcal{F}$
is the first variation of $\mathcal{F}$777The first variation can be defined,
for any measures $\nu_{0},\nu_{1}\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$,
as $\lim_{t\downarrow
0}\frac{\mathcal{F}((1-t)\nu_{0}+t\nu_{1})-\mathcal{F}(\nu_{0})}{t}=\langle\delta\mathcal{F}(\nu_{0}),\nu_{1}-\nu_{0}\rangle$.
This definition is unique up to an additive constant, which is irrelevant as
we are only concerned with its gradient.. Equipped with this
$\mathcal{W}_{2}$-gradient, one can define the _Wasserstein gradient flow_ of
$\mathcal{F}$ that describes the evolution of a measure $\\{\mu_{t}\\}_{t\geq
0}$, from some initial measure $\mu_{0}$, as follows:

$\partial_{t}\mu_{t}=\operatorname{div}\bigl{(}\mu_{t}\nabla_{\mathcal{W}_{2}}\mathcal{F}(\mu_{t})\bigr{)}\,.$

More generally, we can identify the Wasserstein “velocity” for some measure
$\mu_{t}$ as $v_{t}$ if the time derivative of $\mu_{t}$ can be written in the
form

$\partial_{t}\mu_{t}=-\operatorname{div}(\mu_{t}v_{t})\,.$

Under this identification, the time derivative of a functional $\mathcal{F}$
on $\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ with smooth Wasserstein gradient
under these dynamics can be written as

$\partial_{t}\mathcal{F}(\mu_{t})=\mathbb{E}_{\mu_{t}}\langle\nabla_{\mathcal{W}_{2}}\mathcal{F}(\mu_{t}),v_{t}\rangle\,,$

when $v_{t}\in\overline{\\{\nabla\psi:\psi\in
C_{c}^{\infty}(\mathbb{R}^{d})\\}}^{L^{2}(\mu_{t})}$, where
$\overline{\\{\cdot\\}}^{L^{2}(\mu_{t})}$ denotes the closure of a set with
respect to $L^{2}(\mu_{t})$. This is the appropriate notion of tangent space
in this geometry.

For instance, when we take the functional to be the entropy of the measure,
$\mathcal{H}(\mu)\coloneqq\frac{1}{2}\int\mu\log\mu$, one can verify
$\nabla_{\mathcal{W}_{2}}\mathcal{H}(\mu)=\frac{1}{2}\nabla\log\mu$. The heat
flow equation can be written as
$\partial_{t}\mu_{t}=\frac{1}{2}\Delta\mu_{t}=\frac{1}{2}\operatorname{div}(\nabla\mu_{t})=\frac{1}{2}\operatorname{div}(\mu_{t}\nabla\log\mu_{t})$,
which indicates that the velocity of measures $\mu_{t}$ under the heat flow is
$v_{t}=-\frac{1}{2}\nabla\log\mu_{t}$. Hence, we can notice that
$\nabla_{\mathcal{W}_{2}}\mathcal{H}(\mu_{t})=-v_{t}$, and thus recover the
heat flow as the Wasserstein gradient flow of the entropy of the measure.

#### Fokker-Planck equation and time-reversal of SDE.

Consider a stochastic differential equation $(X_{t})$ given by

$\mathrm{d}X_{t}=-a_{t}(X_{t})\,\mathrm{d}t+\,\mathrm{d}B_{t}\qquad\text{with
}X_{0}\sim\mu_{0}\,.$ (D.1)

It is well known that measures $\mu_{t}$ described by

$\partial_{t}\mu_{t}=\mathsf{div}(\mu_{t}a_{t})+\frac{1}{2}\Delta\mu_{t}\,,$
(D.2)

correspond to $\mathsf{law}(X_{t})$. In this context, (D.2) is referred to as
the _Fokker-Planck equation_ corresponding to (D.1).

From this equation, one can deduce the Fokker-Planck equation of the time
reversal $\mu_{t}^{\leftarrow}:=\mu_{T-t}$:

$\partial_{t}\mu_{t}^{\leftarrow}=-\operatorname{div}(\mu_{t}^{\leftarrow}a_{T-t})-\frac{1}{2}\Delta\mu_{t}^{\leftarrow}=-\operatorname{div}\bigl{(}\mu_{t}^{\leftarrow}(a_{T-t}+\nabla\log\mu_{T-t})\bigr{)}+\frac{1}{2}\Delta\mu_{t}^{\leftarrow}$

In particular, this describes the evolution of $\mathsf{law}(X_{t})$ of the
stochastic differential equation:

$\mathrm{d}X_{t}=\bigl{(}a_{T-t}(X_{t})+\nabla\log\mu_{T-t}(X_{t})\bigr{)}\,\mathrm{d}t+\,\mathrm{d}B_{t}\qquad\text{with
}X_{0}\sim\mu_{0}^{\leftarrow}=\mu_{T}\,.$ (D.3)

While the law of this process will give $\mu_{T}^{\leftarrow}=\mu_{0}$ at time
$T$, it is also true that it will give $\mu_{0|T}(\cdot|z)$ if one starts
(D.3) at $X_{0}=z$. This is a subtle fact, whose justification requires the
introduction of a tool called _Doob’s $h$-transform_. The presentation of this
subject is beyond the scope of this paper, and we refer interested readers to
[KP21] as a reference to its application in this context.