# Rényi-infinity constrained sampling with $d^{3}$ membership queries

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

Uniform sampling over a convex body is a fundamental algorithmic problem, yet
the convergence in KL or Rényi divergence of most samplers remains poorly
understood. In this work, we propose a constrained proximal sampler, a
principled and simple algorithm that possesses elegant convergence guarantees.
Leveraging the uniform ergodicity of this sampler, we show that it converges
in the Rényi-infinity divergence ($\mathcal{R}_{\infty}$) with no query
complexity overhead when starting from a warm start. This is the strongest of
commonly considered performance metrics, implying rates in
$\\{\mathcal{R}_{q},\mathsf{KL}\\}$ convergence as special cases.

By applying this sampler within an annealing scheme, we propose an algorithm
which can approximately sample $\varepsilon$-close to the uniform distribution
on convex bodies in $\mathcal{R}_{\infty}$-divergence with
$\widetilde{\mathcal{O}}(d^{3}\,\text{polylog}\frac{1}{\varepsilon})$ query
complexity. This improves on all prior results in
$\\{\mathcal{R}_{q},\mathsf{KL}\\}$-divergences, without resorting to any
algorithmic modifications or post-processing of the sample. It also matches
the prior best known complexity in total variation distance.

##  1 Introduction

Uniform sampling from convex bodies is a fundamental question in computer
science. Its applications have spanned from differential privacy [MT07, HT10,
Mir17] to scientific computing [CV16, HCT+17, KLSV22] and machine learning
[BCJ+19, Sta20].

The standard computational model assumes that, given a convex body
$\mathcal{K}\subset\mathbb{R}^{d}$ which satisfies
$B_{1}(0)\subseteq\mathcal{K}\subseteq B_{D}(0)$, the algorithm can access a
_membership oracle_ [GLS88] which, given a point $x\in\mathbb{R}^{d}$, answers
_Yes_ or _No_ to the query “Is $x\in\mathcal{K}$?”. This framework possesses
two advantages: (i) it is the most general framework under which one can
conduct analysis, subsuming other computational models as particular cases,
(ii) it has been thoroughly studied both in optimization and sampling.

Under this computational model, what is the complexity in terms of membership
oracle queries to generate a random sample from a convex body $\mathcal{K}$
(i.e., from
$\pi=\frac{1}{\operatorname{vol}(\mathcal{K})}\mathds{1}_{\mathcal{K}}$)? As
many of the desired applications require the dimension $d$ to be large, it is
typically impractical to generate samples exactly from the uniform
distribution. Among other difficulties, this is due to the intractability of
deterministically computing the normalizing factor
$\operatorname{vol}(\mathcal{K})$. Instead, one resorts to sampling from an
approximate distribution which is $\varepsilon$-close to the uniform in an
appropriate sense. Subsequent work has aimed to develop algorithms and
analysis with optimal asymptotic dependence both on the dimension $d$ and
error tolerance $\varepsilon$.

The seminal work of [DFK91] proposed a randomized polynomial time algorithm
with oracle complexity $\text{poly}(d,\log\frac{D}{\varepsilon})$. Since then,
a proliferation of improvements to both the algorithm and analysis [LS90,
Lov91, AK91, LS93, KLS97, LV06b, CV18] have led to a refined understanding of
the complexity of this problem. Such algorithms are typically instances of a
random walk, which first samples a point $x_{0}\sim\mu_{0}$ for a tractable
initial measure $\mu_{0}$, and then produces $x_{1},x_{2},\ldots,x_{k_{*}}$
following a well-specified, possibly random procedure for some number of
iterations $k_{*}$. Under mild assumptions, it is possible to make
$\mu_{k_{*}}$, the distribution of $x_{k_{*}}$, $\varepsilon$-close to $\pi$
in some desired sense. These iterative schemes are called often Markov chain
Monte Carlo (MCMC) methods.

The notion of closeness can take many different forms. Natural choices of
metric111In this text, we will often refer to any of these divergences as a
“metric” in the sense of performance metric, although only the total variation
is a true metric on the space of probability measures. such as the
$\mathsf{KL}$ divergence and total variation ($\mathsf{TV}$) fit into the
following hierarchy (see §2 for definitions):

$\mathsf{TV}^{2}\leq\mathsf{KL}\text{ divergence}\leq\mathcal{R}_{q}\text{
divergence}\leq\mathcal{R}_{\infty}\text{ divergence}.$

While the complexity of uniformly sampling convex bodies was known to be
$\widetilde{\mathcal{O}}(d^{3}\log^{2}\frac{1}{\varepsilon})$ for
$\varepsilon$-$\mathsf{TV}$ [CV18, JLLV21], the best known rates in any other
metric had much worse dependence on key parameters. Additionally, they relied
on algorithmic modifications which deviated from practical implementations
[BDMP17, GHZ22, MV22]. On the other hand, the literature on unconstrained
sampling has provided stronger guarantees in $\mathsf{KL}$ or $\chi^{2}$
_without_ overhead in complexity, through refined analyses of Langevin-based
algorithms [VW19, CEL+22, CCSW22]. Adapting those guarantees to the convex
body setting is not trivial, since restricting the distribution to
$\mathcal{K}$ places a strong constraint on the candidate algorithm. The
discrepancies in complexity when going beyond $\mathsf{TV}$ leads us to a
fundamental question on the complexity of uniform sampling in a stronger
metric:

_Is it possible to achieve a query complexity of
$\widetilde{\mathcal{O}}(d^{3}\operatorname{polylog}\frac{1}{\varepsilon})$
for uniform sampling in
$\\{\mathsf{KL},\mathcal{R}_{q},\mathcal{R}_{\infty}\\}$-divergence?_

In general, almost all algorithms take two phases: ($\mathsf{P1}$) _warm-start
generation_ and ($\mathsf{P2}$) _faster sampling under warm-start_.
$\mathsf{P1}$ finds a tractable initial distribution $\mu_{0}$ towards the
target $\pi$, whose closeness to $\pi$ is measured by the
$\mathcal{R}_{\infty}$ divergence at initialization,
$M=\sup\frac{\mathrm{d}\mu_{0}}{\mathrm{d}\pi}=\exp(\mathcal{R}_{\infty}(\mu_{0}\mathbin{\|}\pi))$.
$\mathsf{P2}$ then samples from the target by leveraging the initial warmness.
This approach is also taken by [CV18], the best known uniform-sampling
algorithm in $\mathsf{TV}$. However, it fails to go beyond $\mathsf{TV}$, as a
warm-start generated by the algorithm holds only in a _weak sense_ , causing
“$\mathsf{TV}$-collapse” (we return to this shortly). The resulting sampler
can only have guarantees in $\mathsf{TV}$.

Ignoring $\mathsf{P1}$ for the moment, one hopes that, towards stronger
guarantees, there should be a sampler which converges in $\mathcal{R}_{q}$ in
$\mathsf{P2}$. Recently, [KVZ24] proposed an algorithm with query complexity
$\widetilde{\mathcal{O}}(qMd^{2}\|\Sigma\|_{\mathsf{op}}\operatorname{polylog}\frac{1}{\varepsilon})$
for $\Sigma$ covariance of $\pi$, which seamlessly extends the best complexity
of
$\widetilde{\mathcal{O}}(Md^{2}\|\Sigma\|_{\mathsf{op}}\operatorname{polylog}\frac{1}{\varepsilon})$
in $\mathsf{TV}$ that is achieved by $\mathsf{Ball\ walk}$. However, this
immediately raises a problem: it requires a $\mathcal{R}_{\infty}$ guarantee
at initialization in order to achieve $\mathcal{R}_{q}$-convergence. One might
want to sample from $\pi$ directly from some simple feasible start, such as a
uniform distribution on the unit ball. Unfortunately, such choices potentially
incur exponential warmness in dimension (i.e., $M\gtrsim\exp(d)$). This is
where $\mathsf{P1}$ demonstrates its importance, allowing us to address this
issue through an _annealing scheme_.

##### Annealing algorithms

The strategy of annealing for volume computation and uniform sampling appeared
in [DFK91]. It proceeds with a sequence of distributions $\\{\mu_{k}\\}_{k\leq
k_{*}}$ which steadily evolve from $\mu_{0}\propto\mathds{1}_{B_{1}(0)}$, the
uniform distribution on a unit ball, towards $\mu_{k_{*}}=\pi$ the uniform
distribution on the entire body. This is done by fattening $\mu_{k}$, steadily
increasing its variance, until it is close enough to $\pi$ to be sampled
directly with MCMC.222In the original paper, this is done by a sequence of
uniform distributions, but this strategy has since been superseded by
considering a sequence of truncated Gaussians with appropriate variance. This
general schema, wherein a tractable, well-understood distribution is gradually
converted to a generic one, has been called _annealing_.

This scheme has been refined by annealing via exponential distributions
[LV06c] and via Gaussians. The latter, whose use of
$\\{\mathcal{N}(0,\sigma^{2}_{k}I_{d})|_{\mathcal{K}}\\}_{k\leq k_{*}}$ was
pioneered by the $\mathsf{Gaussian\ cooling}$ algorithm [CV18], turns out to
be a robust choice for $\\{\mu_{k}\\}_{k\leq k_{*}}$ under an appropriate
update schedule for the variance $\sigma_{k}^{2}$. These intermediate
distributions are “stepping stones” towards the desired uniform distribution,
in the sense that each distribution is $\mathcal{O}(1)$ close to the
subsequent one in $\mathcal{R}_{\infty}$, and no more than
$k_{*}=\widetilde{\mathcal{O}}(d)$ annealing stages are required in total.

$\mathsf{Gaussian\ cooling}$ can guarantee $M=\mathcal{O}(1)$ in a weak sense
with the total query complexity
$\widetilde{\mathcal{O}}(C^{2}d^{3}\operatorname{polylog}\frac{1}{\varepsilon})$
under a canonical set-up where $\mathcal{K}$ is _well-rounded_ (i.e.,
$\mathbb{E}_{X\sim\pi}[\|X\|^{2}]\leq C^{2}d$ for some constant $C$) [KLS97,
LV06c, LV06a, CV18]. Later work [JLLV21] claims that one can find an affine
transformation which maps a skewed convex body to a well-rounded one with
$C=\mathcal{O}(1)$, via a randomized algorithm with query complexity
$\widetilde{\mathcal{O}}(d^{3})$.

##### TV-collapse

Such annealing schemes are interwoven with an _approximate_ sampler, typically
an MCMC subroutine such as $\mathsf{Ball\ walk}$ and
$\mathsf{Hit\text{-}and\text{-}Run}$, requiring the sampler to approximately
sample from evolving measures. Hence, the warmness for a next target
$\mu_{k+1}$ is actually measured with respect to a distribution
$\tilde{\mu}_{k}$, approximately close to $\mu_{k}$. Standard constrained
samplers like $\mathsf{Ball\ walk}$ or $\mathsf{Hit\text{-}and\text{-}Run}$
are only able to provide guarantees (between $\tilde{\mu}_{k}$ and $\mu_{k}$)
in $\mathsf{TV}$ for truncated Gaussian. However, due to failure of the
triangle inequality through $\mathsf{TV}$ and $\mathcal{R}_{q}$ in general,
the warmness of $\tilde{\mu}_{k}$ with respect to $\mu_{k+1}$ cannot be
bounded. Even worse, it is essential that the warm-start for these samplers is
in some Rényi divergence, preferably the $\mathcal{R}_{\infty}$ divergence. As
a result, the annealing schemes fail to carry a strong enough notion of
warmness across the evolving measures. Instead, it passes a _weakened form of
warmness_ to the last step via a coupling argument, causing any sampling
guarantees based on this warmness to degrade to $\mathsf{TV}$ distance.

Given this challenge, we state a refinement of ($\mathsf{P1}$):

######  Question 1.

_Can we generate a true warm-start $M=\mathcal{O}(1)$ with
$\widetilde{\mathcal{O}}(d^{3})$ queries?_

The answer to this hinges upon the existence of a sampler which converges in
$\mathcal{R}_{\infty}$. As the triangle inequality holds for
$\mathcal{R}_{\infty}$, this would enable us to anneal the distribution
_without compromising_ $\mathcal{R}_{\infty}$-closeness throughout! Examining
this sends us back to a more specific form of ($\mathsf{P2}$):

######  Question 2.

Given $M=\mathcal{O}(1)$ and uniform or truncated Gaussian target, is there a
sampler with $\mathcal{R}_{\infty}$ guarantees without any overhead in query
complexity?

##### Rényi divergence: theoretical gap and privacy

From our earlier discussion, we cannot help but circle back to the study of
Rényi divergence in constrained sampling. This is theoretically important in
its own right, since it leads us to examine a current gap between
unconstrained and constrained sampling. More specifically, the improvement in
metric from $\mathsf{TV}$ to $\\{\mathsf{KL},\mathcal{R}_{q}\\}$ nicely
parallels the work in unconstrained sampling, where convergence in weaker
metrics (Wasserstein, $\mathsf{TV}$) [DT12, Dal17] paved the way for a
characterization of rates in $\\{\mathsf{KL},\mathcal{R}_{q}\\}$ [CB18,
MFWB22, DMM19, VW19, MCC+21, CEL+22, CCSW22, ZCL+23, AC23]. This research
program led to many fruitful insights into algorithmic properties and
mathematical techniques, and inspires us to attempt the same in the
constrained setting. See §1.3 for detail.

In purely practical terms, the $q$-Rényi divergence, particularly
$\mathcal{R}_{\infty}$, has appeared in the literature most often in the
context of differential privacy (DP) [Mir17]. $\mathcal{R}_{q}$ can be used to
establish $(\epsilon,\delta)$-DP algorithms [HT10], while
$\mathcal{R}_{\infty}$ can directly establish $\epsilon$-DP algorithms. This
connection has been well-explored in other sampling applications [GT20, GLL22,
GLL+23]. However, convergence in $\mathcal{R}_{\infty}$ has remained largely
out of reach, and the existing algorithms are both suboptimal in complexity
and limited in scope [HT10, MV22].

###  1.1 Contributions

We close a gap in our current understanding of uniform sampling via a
principled approach. We first leverage recent results [KVZ24] for sampling on
convex bodies, proposing a generalization of [KVZ24] that we call the
$\mathsf{Proximal\ sampler}$. Using uniform ergodicity results from the study
of Markov chains [DMLM03], we show that it is possible to achieve convergence
in the _very strong_ notion of $\mathcal{R}_{\infty}$ for both uniform and
truncated Gaussian distributions on $\mathcal{K}$ from an $\mathcal{O}(1)$
warm-start without suffering any loss in complexity. By using the
$\mathsf{Proximal\ sampler}$ inside an annealing scheme we call
$\mathsf{Proximal\ Gaussian\ cooling}$, we show that it is possible to sample
uniformly from a well-rounded convex body (i.e.,
$\mathbb{E}_{X\sim\pi}[\|X\|^{2}]\leq C^{2}d$) with
$\widetilde{\mathcal{O}}(C^{2}d^{3}\operatorname{polylog}\frac{1}{\varepsilon})$
queries in expectation. This convergence takes place in
$\mathcal{R}_{\infty}$, much stronger than the prior best results in
$\mathsf{TV}$.

This framework is principled, and the essential benefits of the annealing
scheme and sampler are made clearly evident. Furthermore, the approach is
simple, with no need to resort to algorithmic modifications as seen in other
work [HT10, BDMP17, MV22]. We present our improvements in further detail
below.

##### Result 1: $\mathcal{R}_{\infty}$ convergence from a warm-start

We show that, given an $\mathcal{O}(1)$ warm-start, there exists an algorithm
which samples approximately with $\varepsilon$-accuracy in
$\mathcal{R}_{\infty}$ for any $\varepsilon>0$. The target distribution is
either the uniform distribution on $\mathcal{K}$ or a truncated Gaussian of
the form $\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$. The oracle
complexity is given by the following theorem.

######  Theorem 3 (Complexity of $\mathsf{Proximal\ sampler}$ from a warm-
start, informal version of Theorems 18 and 26).

Consider a target of either of the forms $\pi\propto\mathds{1}_{\mathcal{K}}$,
$\pi_{\sigma^{2}}=\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$, where
$B_{1}(0)\subseteq\mathcal{K}\subseteq B_{D}(0)$. Then, given $\pi_{0}$ with
$\mathcal{R}_{\infty}(\pi_{0}\mathbin{\|}\pi)=\mathcal{O}(1)$, the
$\mathsf{Proximal\ sampler}$ algorithm with probability $1-\eta$ succeeds in
sampling from $\nu$,
$\mathcal{R}_{\infty}(\nu\mathbin{\|}\pi)\leq\varepsilon$, requiring no more
than

$\mathcal{O}\Bigl{(}d^{2}\mathsf{D}^{2}\log^{\mathcal{O}(1)}\frac{D}{\eta\epsilon}\Bigr{)}\,,\qquad\mathsf{D}^{2}:=\begin{cases}\frac{1}{d}\vee\sigma^{2}&\text{if
the target is $\pi_{\sigma^{2}}$}\,,\\\ D^{2}&\text{if the target is
$\pi$}\,,\end{cases}$

membership queries in expectation.

The algorithm is a generalization of the $\mathsf{In\text{-}and\text{-}Out}$
sampler [KVZ24] (Algorithm 1). We call this $\mathsf{Proximal\ sampler}$
(Algorithm 1, 2), in analogy with algorithms in unconstrained sampling
[LST21]. The chief benefit of this scheme is its analytic simplicity, which
allows us to extend known guarantees with little effort. More specifically, it
is naturally connected to isoperimetric properties such as the log-Sobolev
constant of the target. This fact will be key in establishing convergence in
$\mathcal{R}_{\infty}$.

Input: initial point $x_{0}\sim\mu_{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{Prox\text{-}uniform}$

Input: initial point $x_{0}\sim\mu_{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}\bigl{(}\frac{1}{1+h\sigma^{-2}}y,\frac{h}{1+h\sigma^{-2}}I_{d}\bigr{)}\cdot\mathds{1}_{\mathcal{K}}$
until $x_{i+1}\in\mathcal{K}$ or #attemptsi $\geq N$ (declare Failure).

4: end for

Algorithm 2 $\mathsf{Prox\text{-}Gaussian}$

The surprising fact is that this convergence takes place essentially in the
strongest possible sense of $\mathcal{R}_{\infty}$ without any overhead. To
show this, we also demonstrate a strong mixing property of $\mathsf{Proximal\
sampler}$ given any point mass within $\mathcal{K}$ as the initial
distribution. While we are not the only work to obtain rates in this
divergence [HT10, MV22], we are the first to show $\mathcal{R}_{\infty}$
guarantees with efficient rates and without resorting to algorithmic
modifications.

Furthermore, the results for the truncated Gaussian are the first guarantees
in $\mathcal{R}_{q}$ and $\mathcal{R}_{\infty}$ for the non-uniform,
constrained case with $d^{2}\sigma^{2}$ query complexity. While nowhere close
to a comprehensive guarantee for sampling from constrained distributions
$\tilde{\pi}\propto e^{-f}|_{\mathcal{K}}$, our results hint that this generic
problem can be tackled using the same analytic framework.

##### Result 2: Warm-start generation and $\widetilde{\mathcal{O}}(d^{3})$
uniform sampler in $\mathcal{R}_{\infty}$

As a major application of the previous result, we show that, for any well-
rounded convex body $\mathcal{K}$, there is an annealing algorithm with
$\widetilde{\mathcal{O}}(d^{3}\operatorname{polylog}\frac{1}{\varepsilon})$
complexity that achieves $\varepsilon$-accuracy in $\mathcal{R}_{\infty}$.
Compared to $\mathsf{Gaussian\ cooling}$, this algorithm uses
$\mathsf{Proximal\ sampler}$ as a subroutine for sampling at each stage. As a
result, we call it $\mathsf{Proximal\ Gaussian\ cooling}$ (Algorithm 3). This
allows us to strengthen the previous guarantees from [CV18] to
$\mathcal{R}_{\infty}$-divergence, which encompasses all other divergences.

######  Theorem 4 ($\mathcal{R}_{\infty}$ guarantees for uniform sampling,
informal version of Theorem 28).

Let $\pi\propto\mathds{1}_{\mathcal{K}}$ for $\mathcal{K}$ well-rounded, i.e.
$\mathbb{E}_{X\sim\pi}[\|X\|^{2}]\leq C^{2}d$ while
$B_{1}(0)\subseteq\mathcal{K}$. Then, $\mathsf{Proximal\ Gaussian\ cooling}$
succeeds in sampling from $\nu$ with probability $1-\eta$, where
$\mathcal{R}_{\infty}(\nu\mathbin{\|}\pi)\leq\varepsilon$, requiring no more
than
$\mathcal{O}\bigl{(}C^{2}d^{3}\log^{\mathcal{O}(1)}\frac{1}{\eta\varepsilon}\bigr{)}$
expected number of membership queries.

We emphasize that our rate of $\widetilde{\mathcal{O}}(d^{3})$ matches that of
[CV18] while being in a much stronger metric ($\mathcal{R}_{\infty}$). To
compare, the expected rate of the proximal sampler for unconstrained
distributions $\nu$ is $\widetilde{\mathcal{O}}(\beta
C_{\mathsf{PI}}(\pi)d^{3/2})$ from a feasible start [AC23, Remark 5.5], where
$\beta$ is the Lipschitz constant of $\nabla\log\pi$. By contrast, Theorem 4
states that, in exchange for letting $\beta=\infty$, we need to pay only an
extra multiplicative factor of $d^{3/2}C_{\mathsf{PI}}^{-1}$ when sampling
from a well-rounded convex body.

###  1.2 Techniques and challenges

####  1.2.1 Constrained samplers

##### Uniform sampling from a warm-start

Before worrying about the problem of obtaining a warm start, it is necessary
to establish the best possible oracle complexity for the sampling routine in
$\mathcal{R}_{\infty}$ when given an $\mathcal{O}(1)$ warm-start.

Previous works such as [LV06c, CV18] relied on
$\mathsf{Hit\text{-}and\text{-}Run},\mathsf{Ball\ walk}$ as underlying
samplers with their annealing schemes. These possess several disadvantages.
Firstly, the method of analysis is rather complicated. As for $\mathsf{Ball\
walk}$, the best analysis (in terms of complexity) makes statements about
$\mathsf{Speedy\ walk}$, which is a biased variant of $\mathsf{Ball\ walk}$
[KLS97, KLM06]. To implement $\mathsf{Ball\ walk}$ with $\mathsf{Speedy\
walk}$, an additional correction step for debiasing the stationary measure is
required. This degrades rate estimates for $\mathsf{Speedy\ walk}$ to
$\mathsf{TV}$, making the rate estimates for $\mathsf{Ball\ walk}$ only in
$\mathsf{TV}$ as a result. $\mathsf{Hit\text{-}and\text{-}Run}$, on the other
hand, is only known to mix in $\\{\mathsf{TV},\chi^{2}\\}$ for uniform and
exponential distributions [LV06b] and in $\mathsf{TV}$ for general log-concave
distributions [LV06a].

By contrast, the recent work of [KVZ24] demonstrated that
$\mathsf{In\text{-}and\text{-}Out}$ (Algorithm 1) could achieve
$\varepsilon$-mixing in any $\mathcal{R}_{q}$, requiring
$\widetilde{\mathcal{O}}(qMd^{2}\|\Sigma\|_{\mathsf{op}}\operatorname{polylog}{\textstyle\frac{1}{\varepsilon}})$
queries in expectation from an $M$-warm start.

#####  $\mathsf{In\text{-}and\text{-}Out}$ and the $\mathsf{Proximal\
sampler}$

The algorithm of $\mathsf{In\text{-}and\text{-}Out}$ (or
$\mathsf{Prox\text{-}uniform}$) is remarkably simple. It alternately samples
from two distributions, first getting $y$ from a Gaussian centred at the
current point, then $z$ from a Gaussian centred at $y$ but truncated to
$\mathcal{K}$ (Lines 2 and 3 in Algorithm 1). This turns out to correspond to
simulating the heat equation applied to the current measure, and then
simulating a notion of time reversal of the heat equation (Lemma 20). These
dynamics cause any $\mathcal{R}_{q}$ divergence to decay exponentially, with
the speed of decay dependent on the isoperimetry of the target distribution.

The only remaining challenge is to bound the wasted queries in Line 3. This
follows a local conductance argument, where one quantifies the probability of
$y$ landing in a set that has a hard time returning to $\mathcal{K}$. This
step requires both that $\mu_{0}$ be warm (for reasons to be explained later)
and that the variance $h=\widetilde{\mathcal{O}}(d^{-2})$ is appropriately
small.

##### Gaussian sampling with $\mathsf{Proximal\ sampler}$

Since we will end up annealing with Gaussians, we need to demonstrate that the
guarantees for $\mathsf{In\text{-}and\text{-}Out}$ are compatible with
(truncated) Gaussian distributions
$\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$. This is not necessarily
trivial. While the rate estimates for $\mathsf{Ball\ walk}$ can be adapted for
Gaussians [CV18], those for $\mathsf{Hit\text{-}and\text{-}Run}$ cannot be
[LV06a]. To use $\mathsf{Hit\text{-}and\text{-}Run}$ in an annealing scheme,
[LV06c] resorts to exponential distributions, which have worse complexity.

In this work, we generalize the $\mathsf{In\text{-}and\text{-}Out}$ algorithm
for truncated Gaussians on $\mathcal{K}$. This can be done by adapting the
forward and backward heat flow interpretations to incorporate a potential
function in the target, in analogy with works in unconstrained sampling
[CCSW22]. By the same analogy, we call our algorithm the (constrained)
$\mathsf{Proximal\ sampler}$.

While the calculations for the Gaussian case are somewhat involved, the
approach to bounding the query complexity does not require any algorithmic
modifications instead. The number of iterations needed to mix in
$\mathcal{R}_{q}$ for the $\mathsf{Proximal\ sampler}$,
$\widetilde{\mathcal{O}}(q(d\vee d^{2}\sigma^{2}))$, already follows as a
consequence of a generic lemma shown in [KVZ24] (restated in Lemma 21). As for
the query complexity, a careful bound on the local conductance following the
same approach as [KVZ24] (Lemma 23) shows that $\widetilde{\mathcal{O}}(M)$
membership queries suffices in each iteration. This makes
$\widetilde{\mathcal{O}}(qM(d\vee d^{2}\sigma^{2}))$ queries in total. The
complexity matches that of [CV18] for Gaussians in $\mathsf{TV}$, and the ease
of derivation underscores another strength of the $\mathsf{Proximal\
sampler}$/$\mathsf{In\text{-}and\text{-}Out}$ framework.

We note that this is the first time that $d^{3}$ guarantees have been shown in
$\mathcal{R}_{q}$ (and, as we shall see, in $\mathcal{R}_{\infty}$ as well)
for truncated Gaussians, even if just for the specific subclass of those that
are centered and isotropic. The techniques here suggest that, with more
effort, the analysis could be conducted more generally for other distributions
of the form $\pi=e^{-f}|_{\mathcal{K}}$.

####  1.2.2 Warm starts and Rényi divergence

##### Importance of warm starts

We briefly digress to explain why warm starts are necessary in each of these
contexts. The above samplers repeat a subroutine wherein a point is proposed,
either from uniform distribution in a ball for $\mathsf{Ball\ walk}$, and a
Gaussian for $\mathsf{In\text{-}and\text{-}Out}$. However, naïvely using each
of these points will bias the sampling algorithm. Thus, the subroutines must
contain rejection steps which correct the distribution of the proposal. To
preserve our complexity bounds, we must ensure that the number of rejections
is small. Otherwise, the algorithm will waste the majority of its queries
without moving.

How to control the number of rejections? If
$\mathsf{In\text{-}and\text{-}Out}$ is initialized at the target
$\pi\propto\mathds{1}_{\mathcal{K}}$ with suitable algorithmic parameters,
then one can bound the expected number of rejections by
$\widetilde{\mathcal{O}}(1)$. Instead, with $\exp(\mathcal{R}_{\infty})\leq M$
at initialization, a straightforward calculation bounds the number of
rejections by that experienced by $\pi$, multiplied by a factor of $M$. This
explicitly reveals the importance of starting warm when bounding the number of
rejections.

Given this context, we denote by $\mathsf{M}\to\mathsf{N}$ when a sampler
starting with warmness in a metric $\mathsf{M}$ has convergence guarantees in
a possibly weaker metric $\mathsf{N}$. The guarantees of the previous samplers
can be stated for the truncated Gaussian as $\mathsf{Ball\
walk},\mathsf{Hit\text{-}and\text{-}Run}:\mathcal{R}_{\infty}\to\mathsf{TV}$
and $\mathsf{Proximal\ sampler}:\mathcal{R}_{\infty}\to\mathcal{R}_{q}$. The
previous discussion implies that we would require
$\mathcal{R}_{\infty}\to\mathcal{R}_{\infty}$ to relay
$\mathcal{R}_{\infty}$-guarantees across the annealing routine. Of the
aforementioned works, only $\mathsf{Proximal\ sampler}$ is close to achieving
this. Nonetheless, closing this gap from $\mathcal{R}_{q}$ to
$\mathcal{R}_{\infty}$ is far from trivial and requires the introduction of
new analytical techniques, as we subsequently explain.

##### The difficulties of $\mathcal{R}_{\infty}$

One would expect that convergence in $\mathcal{R}_{\infty}$ is extremely
difficult to show. By comparison, even for a continuous-time process such as
the Langevin dynamics, convergence in $\mathcal{R}_{\infty}$ is not known
outside a few specific cases. Previous sampling guarantees in
$\mathcal{R}_{q}$ usually involve a complexity at least linear in $q$ [CEL+22,
CCSW22, ZCL+23, AC23, KVZ24], which render them useless for
$\mathcal{R}_{\infty}$.

The techniques for analyzing $q$-Rényi divergence usually involve constructing
a Markov semigroup whose interpolation in discrete time yields the sampler of
interest. The time evolution equations for this semigroup can then be used to
derive exponential convergence of $\mathcal{R}_{q}$, with a constant depending
on the isoperimetry of the stationary measure for the semigroup. This approach
fails miserably for $\mathcal{R}_{\infty}$. In its standard representation,
$\mathcal{R}_{\infty}$, being the supremum of the density ratio, cannot be
written as the expectation of a continuous quantity, and its decay properties
are difficult to establish. Instead, one must look for proof strategies that
use more information about the semigroup.

##### Any-start implies $\mathcal{R}_{\infty}$

It is well known in the Markov chain literature that, if a Markov chain mixes
rapidly in $\mathsf{TV}$ from any deterministic starting point (a property
known as _uniform ergodicity_), then the Markov chain causes the density to
contract towards that of $\pi$ in any arbitrarily strong norm. In particular,
choosing the $L^{\infty}$-norm, it is possible to translate statements in a
weak norm like “the algorithm converges to within $\varepsilon$ in
$\mathsf{TV}$ distance of $\pi$ within $K$ iterations, given _any starting
point_ in the support of $\pi$” to an unconditional statement in an _extremely
strong_ norm — “the algorithm mixes in $L^{\infty}$ distance within $K$
iterations” (Theorem 14). We call this analytical technique boosting for its
remarkable strengthening of the performance metric.

This result had thus far been difficult to apply in most sampling settings.
Consider for instance sampling from an unconstrained standard Gaussian, via
some MCMC method like the Langevin diffusion. The $\mathsf{TV}$ can always be
bounded by the $\mathsf{KL}$-divergence, which contracts exponentially.
However, because $\pi$ is supported on $\mathbb{R}^{n}$, it is always possible
to pick a starting point arbitrarily far away from the Gaussian’s mode. This
causes $\mathsf{KL}$ at initialization to be arbitrarily large. The number of
iterations that are needed to converge from a point-mass at $x$ is roughly
$\gtrsim\log\|x\|$, which is not bounded for all starting points
$x\in\mathbb{R}^{d}$. Thus, no finite number of iterations of the Markov chain
will be sufficient to show uniform ergodicity across the entirety of
$\mathbb{R}^{d}$, even though the dynamics converge exponentially quickly in
any $\mathcal{R}_{q}$.

In passing, we also mention another notion of boosting from $\mathsf{TV}$ to
$\mathcal{R}_{\infty}$. This is done via a post-processing step that can be
found in the work of [MV22]. Examining the guarantees of this approach,
however, shows that the query complexity of the original algorithm unavoidably
increases by a factor of $\operatorname{poly}(d)$. Thus, even if one could
obtain warm-start at every step, this approach would give at minimum
$\widetilde{\mathcal{O}}(d^{5})$ query complexity333One can use the
$\mathsf{Gaussian\ cooling}$ algorithm followed by the post-processing of
[MV22]. The claimed complexity follows from the $\mathsf{Gaussian\ cooling}$’s
complexity of $\widetilde{\mathcal{O}}(d^{3}\log^{2}\frac{1}{\epsilon})$ along
with $\epsilon\lesssim\exp(-d)\varepsilon$. See §1.3 for detail. for the
entire algorithm. Our approach, in contrast, does not require any algorithmic
modifications, but rather extracts a convergence guarantee that is already
latent in the algorithm.

##### Proximal sampler mixes from any start

Thus, it suffices to establish convergence guarantees given any deterministic
starting point within $\mathcal{K}$. When attempting to do this, a major issue
arises in the analysis of the classic $\mathsf{Gaussian\ cooling}$, which uses
$\mathsf{Ball\ walk}$. The mixing of $\mathsf{Ball\ walk}$ is established
using a conductance-oriented proof. If the $\mathsf{Ball\ walk}$ is analyzed
directly via $s$-conductance, however, then the warmness of the initial
distribution shows up directly in the number of $\mathsf{Ball\ walk}$ steps
(rather than just the query complexity) [LV07]. This is problematic, since if
we started at a bad initial point (for example, the tip of a cone), the
resulting number of $\mathsf{Ball\ walk}$ steps can be exponentially large in
dimension $d$, making it impossible to invoke uniform ergodicity without
degrading the rate. The other way of implementing $\mathsf{Ball\ walk}$
through $\mathsf{Speedy\ walk}$ [KLS97, KLM06], requires a non-Markovian
correction step, which rules out the boosting technique. Another candidate,
$\mathsf{Hit\text{-}and\text{-}Run}$, uses a weaker isoperimetric property
(more precisely, a Cheeger inequality), so it depends poly-logarithmically on
the warmness [LV06b, LV06a]. Hence, this disqualifies
$\mathsf{Hit\text{-}and\text{-}Run}$ as a sampler due to additional overhead
of $\text{poly}(d)$ from any feasible start.

For these discrete-time samplers, one could try to directly prove a _modified
logarithmic-Sobolev inequality_ (MLSI) of their corresponding Markov kernel
(as opposed to for the stationary distribution $\pi$). While this approach can
potentially bound the complexity overhead of these samplers from any feasible
start by $\operatorname{polylog}(d)$ (or worse), they require an involved
analysis of the kernel as well as the target $\pi$. In contrast,
$\mathsf{Proximal\ sampler}$ can be interpreted in terms of continuous forward
and backward heat flows, so its mixing can be characterized by functional
inequalities of just the target $\pi$. Thus, $\mathsf{Proximal\ sampler}$
enjoys substantially simpler analysis, with richer connections to established
functional inequalities for log-concave measures.

Using $\mathsf{Proximal\ sampler}$ (Algorithm 1 and 2), we can indeed achieve
any-start results without overhead, thereby obtaining $\mathcal{R}_{\infty}$
guarantees via the boosting. Leveraging the log-Sobolev inequality for the
target, the number of iterations needs only depend doubly-logarithmically on
the warmness. From a point-mass, it can be established that the warmness
(after one iteration of $\mathsf{Proximal\ sampler}$) is exponential in
dimension, so this allows one to pay only a logarithmic factor when
establishing mixing from any point-mass within $\mathcal{K}$. These factors
combined allow $\mathsf{Proximal\ sampler}$ to do what $\mathsf{Ball\
walk},\mathsf{Hit\text{-}and\text{-}Run}$ could not: obtain
$\mathcal{R}_{\infty}$ guarantees without any computational overhead or post-
processing. To our knowledge, our paper is the first work in MCMC to
accomplish this, and we hope that our methods would serve as a blueprint for
obtaining similar guarantees, at least in the constrained sampling setting.

####  1.2.3 Proximal Gaussian cooling

At its core, our approach in the annealing part does not deviate from [CV18].
The core idea is that, for a family of distributions
$\\{\mu_{i}\\}_{i\in\\{0,\ldots,k_{*}\\}}$, one must balance two factors: (i)
Each distribution must be $\mathcal{O}(1)$-warm with respect to the succeeding
distribution, and (ii) the total number of distributions should not be larger
than $\widetilde{\mathcal{O}}(d)$. The choice of Gaussian distributions gives
more flexibility in this respect, since one is able to update the covariance
with accelerating speed. The $\mathsf{Proximal\ Gaussian\ cooling}$ algorithm,
including the choices of covariance, is given in Algorithm 3. We note that
$\mathsf{Gaussian\ cooling}$ with $\mathsf{Ball\ walk}$ sidestepped the issue
of warmness by coupling together approximate samplers. This strategy, however,
reduces all the bounds down to $\mathsf{TV}$, and cannot be applied when
examining a stronger metric.

With our approach, by combining all the aforementioned technical ingredients,
we obtain a $\widetilde{\mathcal{O}}(d^{3})$ algorithm in
$\mathcal{R}_{\infty}$ for a well-rounded convex body. One additional benefit
of our scheme is its high-level simplicity: we only need to implement the
$\mathsf{Proximal\ sampler}$ for different targets, which can be viewed as
instances of a unified algorithm. Using this, we then follow the Gaussian-
based annealing strategy.

Input: convex body $\mathcal{K}\subset\mathbb{R}^{d}$ with $0\in\mathcal{K}$
and $\mathbb{E}_{\mathcal{K}}[\|X\|^{2}]\leq C^{2}d$.

Output: approximately uniform sample $z$.

1: Let $\bar{\mathcal{K}}=\mathcal{K}\cap B_{L\sqrt{d}}(0)$ with
$L=C\log\frac{1}{\varepsilon}$, and $\sigma^{2}:=1/d$. Denote
$\pi_{\sigma^{2}}:=\mathcal{N}(0,\sigma^{2}I_{d})|_{\bar{\mathcal{K}}}$.

2: Sample $z\sim\textsf{Unif}\,(B_{1}(0))$.

3: Get $z\sim\mathsf{Prox\text{-}Gaussian}$ with initial point $z$, target
dist. $\pi_{\sigma^{2}}$, target accuracy $\log 2$ (in
$\mathcal{R}_{\infty}$), and success probability
$\eta/\widetilde{\mathcal{O}}(L^{2}d)$.

4: while $1/d\leq\sigma^{2}\leq d$ do

5: Get $z$ from $\mathsf{Prox\text{-}Gaussian}$ with the same setup. Then
update
$\sigma^{2}\leftarrow\begin{cases}\sigma^{2}\bigl{(}1+\frac{1}{d}\bigr{)}&\text{if
}\frac{1}{d}\leq\sigma^{2}\leq 1\,.\\\
\sigma^{2}\bigl{(}1+\frac{\sigma^{2}}{L^{2}d}\bigr{)}&\text{if
}1\leq\sigma^{2}\leq L^{2}d\,.\end{cases}$

6: end while

7: Return $z\sim\mathsf{Prox\text{-}uniform}$ with initial point $z$, target
dist. $\pi^{\bar{\mathcal{K}}}$, target accuracy $\varepsilon$, and success
probability $\eta/\widetilde{\mathcal{O}}(L^{2}d)$

Algorithm 3 $\mathsf{Proximal\ Gaussian\ cooling}$

###  1.3 Related work

##### Constrained samplers

The previous approaches to constrained sampling via random walks date to
[DFK91]. The most well-studied algorithms are $\mathsf{Ball\ walk}$ [LS93,
KLS97] and $\mathsf{Hit\text{-}and\text{-}Run}$ [Smi84, Lov99]. $\mathsf{Ball\
walk}$ is a simple scheme which samples $y$ from a ball around the current
point, accepting $y$ if $y\in\mathcal{K}$ and otherwise just remaining at the
current point. Although it is possible to analyze $\mathsf{Ball\ walk}$
directly [LV07], usually $\mathsf{Ball\ walk}$ refers to the the algorithm
given by composing $\mathsf{Speedy\ walk}$ with rejection sampling. The best
known guarantee in [KLM06] gives a complexity of
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log\frac{1}{\varepsilon})$ with
$M=\exp\mathcal{R}_{\infty}(\mu_{0}\mathbin{\|}\pi)$.

Another standard algorithm is $\mathsf{Hit\text{-}and\text{-}Run}$, which
samples uniformly at random from a chord $\mathcal{K}\cap\ell$, where $\ell$
is a random line passing through the current point. Its query complexity for
uniform and exponential distributions is
$\widetilde{\mathcal{O}}(d^{2}\operatorname{tr}\Sigma\log\frac{M}{\varepsilon})$,
where $\Sigma=\mathbb{E}_{\pi}[(X-\mathbb{E}X)^{\otimes 2}]$. See §1.2 for the
benefits and drawbacks of these two approaches.

Both $\mathsf{Ball\ walk}$ and $\mathsf{Hit\text{-}and\text{-}Run}$ can assume
that $\mathcal{K}$ is well-rounded, in the sense that
$\mathbb{E}_{\pi}[\|X\|^{2}]\leq C^{2}d$ for a known constant $C$. A general
convex body $\mathcal{K}$ can be converted to a well-rounded body via an
affine transformation. Finding this transformation requires an additional
algorithmic ingredient called _rounding_. It has been claimed in the
literature [JLLV21] that a randomized rounding algorithm exists with
$\widetilde{\mathcal{O}}(d^{3})$ complexity. Since this adds no computational
overhead, we make the same assumption in the present work.

Apart from these general-purpose random walks, there also exist several
samplers which exploit geometric information about $\mathcal{K}$, or which use
stronger oracle models. The $\mathsf{Dikin\ walk}$ algorithm makes use of a
self-concordant barrier function $\phi$ to draw samples from an anisotropic
Gaussian (instead of isotropic as in $\mathsf{Ball\ walk},\mathsf{Proximal\
sampler}$), and converges in $\widetilde{\mathcal{O}}(md)$ iterations for a
convex body specified by $m$ linear constraints [KN12, KV24]. Apart from this,
Riemannian algorithms such as $\mathsf{Riemannian\ Hamiltonian\ Monte\ Carlo}$
or $\mathsf{Riemannian\ Langevin\ Monte\ Carlo}$ equip $\mathcal{K}$ with a
Riemannian metric, and then run a random walk using this geometry [GC11, LV18,
LE23, KLSV23, CZS22, GKV23]. Likewise, $\mathsf{Mirror\ Langevin}$ [ZPFP20,
Jia21, AC21, LTVW22, SWW24] alters the geometry so that one can apply methods
from unconstrained sampling, and can obtain $\widetilde{\mathcal{O}}(d)$
complexity with some additional assumptions.

Constrained samplers like [BEL18, Leh23] borrow algorithms from the
unconstrained setting, and interweave them with projection steps onto
$\mathcal{K}$. The best known overall complexity is
$\widetilde{\mathcal{O}}(\frac{d^{2}D^{3}}{\varepsilon^{4}})$ in terms of
projection oracle queries. Finally, [BDMP17, GHZ22] use “soft” penalties
rather than projections. All of the aforementioned techniques, however,
require either additional assumptions on $\mathcal{K}$ or the oracle model, or
are otherwise inefficient terms of complexity.

##### Annealing strategies

Annealing as a computational tool dates back to [KGJV83]. Its original
application was to combinatorial optimization, where the landscape of
solutions is highly non-convex, and it has served as an efficient strategy for
diverse problems [HJJ03, KV06, DCM19].

The history of annealing for volume computation / uniform sampling dates back
to the inception of constrained sampling, from the earliest works of [DFK91]
through a long line of further improvements [Lov91, LS90, AK91, LS93, KLS95,
KLS97]. The original incarnation of this algorithm, which remained unchanged
through the listed references, samples uniformly from a sequence of bodies
$\\{\mathcal{K}_{k}\\}_{1\leq k\leq k_{*}}$, with
$\mathcal{K}_{k}=\mathcal{K}\cap 2^{k/d}B_{1}(0)$, with
$k_{*}=\widetilde{\mathcal{O}}(d)$ iterations. The distribution at each prior
phase is used as a warm start for the subsequent phase. Combined with a
coupling argument for the approximate sampler, the law of the sample at
iteration $k_{*}$ can be viewed as a warm start for the uniform distribution
on $\mathcal{K}$. While the number of phases is not severe, the bodies become
increasingly difficult to sample as $k$ increases. For instance, the
complexity with $\mathsf{Proximal\ sampler}$ moves from
$\widetilde{\mathcal{O}}(d^{2})$ in the earliest phases to
$\widetilde{\mathcal{O}}(d^{3})$ in the latest phases, with the total
complexity being $\widetilde{\mathcal{O}}(d^{4})$ as a result.

To rectify this, [LV06c] proposed a similar multi-phase sampling strategy, but
successively sampling with the exponential distributions $\mu_{k}\propto
e^{-a_{k}^{\mathsf{T}}x}|_{\mathcal{K}}$ where $\|a_{k}\|=2d(1-d^{-1/2})^{k}$.
This needs only $k_{*}=\widetilde{\mathcal{O}}(\sqrt{d})$ phases to converge.
However, these exponential distributions become increasingly ill-conditioned
as $k$ increases, and the warmness between $\mu_{k},\mu_{k+1}$ is only in
$\chi^{2}$. As a result, the best known complexity for obtaining a single
sample with this approach is $\widetilde{\mathcal{O}}(d^{7/2})$.

[CV18] proposes $\mathsf{Gaussian\ cooling}$, which avoids the pitfalls of
both earlier approaches. Each stage $\mu_{k}$ consists of a truncated Gaussian
with variance $\sigma_{k}^{2}$. As opposed to the previous annealing schemes,
it is possible to accelerate the schedule (in the sense that
$\sigma_{k}/\sigma_{k-1}$ is increasing in $k$) while still maintaining
$\mathcal{R}_{\infty}$-warmness between $\mu_{k},\mu_{k+1}$. As a result,
while the total number of phases is $\widetilde{\mathcal{O}}(d)$, the overall
complexity remains at $\widetilde{\mathcal{O}}(d^{3})$.

##### Warm starts and Rényi divergence

The convergence of samplers in Rényi divergences has been well studied in
unconstrained sampling, beginning with results in continuous time [CLL19,
VW19]. These were followed by Rényi guarantees of discrete time algorithms
under log-Sobolev inequalities [GT20, EHZ22, CEL+22], Poincaré inequalities
[Leh23, CEL+22] and beyond [MHFH+23]. These algorithms are vital for warm
starts in the unconstrained field as well, particularly for Metropolis
adjusted algorithms [DCWY18, CDWY18, CLA+21, CG23]. In fact, the study of
$\mathcal{R}_{q}$ divergences and their connection to warm starts directly led
to the present state-of-the-art algorithms in unconstrained sampling [AC23].

Within constrained sampling, the work of [HT10] proposed one method for
obtaining guarantees in $\mathcal{R}_{\infty}$ on transformations of the unit
ball. This approach is based on grid walk, and unrelated to that in the
current work. It has been extended by [BST14] to more general distributions.
Secondly, [MV22] show a boosting scheme which mollifies the convex body by
convolving it with a ball. Such a technique uses a sampler with
$\mathcal{O}(e^{-d}\varepsilon)$-guarantees in $\mathsf{TV}$ as an inner loop
of the routine, which adds at least $\operatorname{poly}(d)$ overhead when
applied to a sampler with $\operatorname{polylog}\frac{1}{\varepsilon}$
dependence.

The technique of proving $L^{\infty}$ bounds using uniform ergodicity was well
known in the study of Markov chains. Its history dates perhaps back to the
works of Markov himself [Mar06], among other venerable works [YK41, Doo53,
MT12]. It is connected with Doeblin’s minorization condition and other
fundamental properties of Markov chains [MT12, §16], and has perhaps been most
succinctly stated in [Rud11, Proposition 3.23]. The uniform ergodicity
property, however, is difficult to establish, barring some exceptions such as
exact Gibbs sampling or samplers on the lattice [LDM15, WL17]. As far as we
know, we are the first to apply it in the constrained sampling setting.

###  1.4 Organization

The remainder of this paper is organized as follows. In §2, we introduce key
notions that will be used in the paper, and in §3 elucidate the boosting
technique, which converts any-start guarantees in $\mathsf{TV}$ to a
contraction in $\mathcal{R}_{\infty}$. We then illustrate how it can be
applied for uniform distributions on $\mathcal{K}$, as a teaser for §4, where
we develop guarantees for truncated Gaussians. Finally, we put these all
together in §5 to obtain our final, $\widetilde{\mathcal{O}}(d^{3})$ uniform
sampling guarantees before concluding.

##  2 Preliminaries

###  2.1 Notation

Unless otherwise specified, $\|\cdot\|$ denotes the $2$-norm on
$\mathbb{R}^{d}$ and the operator norm on $\mathbb{R}^{d\times d}$. The
notation $a=\mathcal{O}(b),a\lesssim b$ will be used to signify that $a\leq
cb$ for $c>0$ a universal constant, while $a\gtrsim b,a=\Omega(b)$ will denote
$a\geq cb$. $a=\Theta(b)$ is used when $a=\Omega(b),a=\mathcal{O}(b)$
simultaneously. The notation $\widetilde{\mathcal{O}}(b)$ will mean
$a=\mathcal{O}(b\,\operatorname{polylog}(b))$, and likewise for
$\widetilde{\Theta},\widetilde{\Omega}$. Finally, we conflate a measure with
its density where there is no confusion.

###  2.2 Basic notions

Before proceeding, we reiterate our computational model.

######  Definition 5 (Membership oracle).

We are given a convex body $\mathcal{K}$ which has
$B_{1}(0)\subseteq\mathcal{K}\subseteq B_{D}(0)\subset\mathbb{R}^{d}$ for some
$D>0$. We assume access to a _membership oracle_ , which, given a point
$x\in\mathbb{R}^{d}$, answers _Yes_ or _No_ to the query “Is
$x\in\mathcal{K}$?”

Where not otherwise specified, we will write
$\pi=\frac{1}{\operatorname{vol}(\mathcal{K})}\mathds{1}_{\mathcal{K}}$ for
the uniform distribution on $\mathcal{K}$. We introduce the following metrics
between probability measures.

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

Let $\mu,\nu$ be two probability measures on $\mathbb{R}^{d}$. Their _total
variation_ distance is given by

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

where $\mathcal{B}(\mathbb{R}^{d})$ is the collection of all Borel measurable
subsets of $\mathbb{R}^{d}$. The _$q$ -Rényi divergence_ is defined as

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

if $\mu\ll\nu$, and $\mathcal{R}_{q}(\mu\mathbin{\|}\nu):=\infty$ otherwise.
In the limit $q\to 1$, it converges to the $\mathsf{KL}$ divergence,

$\mathsf{KL}(\mu\mathbin{\|}\nu):=\int\log\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\,\mathrm{d}\mu\,,$

if $\mu\ll\nu$, and again $\mathsf{KL}(\mu\mathbin{\|}\nu):=\infty$ otherwise.
The $\chi^{2}$-divergence is defined as
$\chi^{2}(\mu\mathbin{\|}\pi):=\exp(\mathcal{R}_{2}(\mu\mathbin{\|}\nu))-1$.
Finally, the limit $\lim_{q\to\infty}\mathcal{R}_{q}(\mu\mathbin{\|}\nu)$ can
be written as

$\mathcal{R}_{\infty}(\mu\mathbin{\|}\nu):=\operatorname{ess\,sup}_{\nu}\log\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\,.$

The $\mathcal{R}_{\infty}$ distance is especially important for us, through
the concept of warmness.

######  Definition 7 (Warmness).

For $\mu\ll\nu$, we denote the warmness $M$ of $\mu$ with respect to $\nu$ as

$M:=\operatorname{ess\,sup}_{\nu}\frac{\mathrm{d}\mu}{\mathrm{d}\nu}=\exp\bigl{(}\mathcal{R}_{\infty}(\mu\mathbin{\|}\nu)\bigr{)}\,.$

Alternatively, if the above holds, then we also say that $\mu$ is an $M$-warm
start for $\nu$.

The following inequality for Rényi divergences will also be useful.

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

Let $\mu,\nu$ be probability measures, $P$ a Markov kernel, and $q\geq 1$.
Then,

$\mathcal{R}_{q}(\mu P\mathbin{\|}\nu
P)\leq\mathcal{R}_{q}(\mu\mathbin{\|}\nu)\,.$

###  2.3 Functional inequalities for constrained distributions

The following functional inequalities, also known as isoperimetry
inequalities, on the target distribution are vital for the analysis of
sampling algorithms, being equivalent to the exponential mixing of a broad
class of Markov chains.

######  Definition 9 (Log-Sobolev inequality).

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

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

with $\operatorname{Ent}_{\nu}(f^{2}):=\mathbb{E}_{\nu}[f^{2}\log
f^{2}]-\mathbb{E}_{\nu}[f^{2}]\log\mathbb{E}_{\nu}[f^{2}]$.

A weaker form of isoperimetry that is implied by the above is the Poincaré
inequality.

######  Definition 10 (Poincaré inequality).

A probability measure $\nu$ satisfies a _Poincaré inequality_ ($\mathsf{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}]\,,$ (PI)

with $\mathrm{Var}_{\nu}f:=\mathbb{E}_{\nu}[(f-\mathbb{E}_{\nu}f)^{2}]$.

The following lemma summarizes the functional inequalities which are satisfied
by (truncated) Gaussians and uniform distributions on $\mathcal{K}$, due to
[Kla23, BGL14]. We refer readers to [KVZ24, Appendix C].

######  Lemma 11.

Let $\mathcal{K}\subset\mathbb{R}^{d}$ be a convex body with diameter $D$.
Then, if $\pi\propto\mathds{1}_{\mathcal{K}}$ is the uniform distribution on
$\mathcal{K}$,

$C_{\mathsf{LSI}}(\pi)\lesssim D^{2}\qquad\text{and}\qquad
C_{\mathsf{PI}}(\pi)\lesssim\|\Sigma\|\log d\,,$

where
$\Sigma=\mathbb{E}_{X\sim\pi}[(X-\mathbb{E}_{\pi}X)(X-\mathbb{E}_{\pi}X)^{\mathsf{T}}]$
is the covariance of $\pi$. For a Gaussian
$\pi_{\sigma^{2}}=\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$, we have

$C_{\mathsf{PI}}(\pi_{\sigma^{2}})\leq
C_{\mathsf{LSI}}(\pi_{\sigma^{2}})\leq\sigma^{2}\,.$

##  3 Total variation to Rényi infinity via LSI

In this section, we first review the work of [Rud11], which allows us to bound
the distance between the iterates of a Markov chain and their stationary
measure in $L^{\infty}$ (and so in $\mathcal{R}_{\infty}$) by the worst-case
distance in $\mathsf{TV}$ from any start. Then we reveal a connection between
(LSI) and convergence from any start (uniform ergodicity), demonstrating a way
to leverage this powerful boosting technique. As a concrete example, we
provide a $\mathcal{R}_{\infty}$-guarantee of $\mathsf{Proximal\ sampler}$ for
uniform sampling over a convex body without incurring additional factors to
the state-of-the-art complexity.

###  3.1 Strong contraction of a Markov chain

We recall standard notion in the theory of Markov semigroups.

######  Definition 12 (Markov kernel).

Let $(\Omega,\mathcal{F})$ be a measurable space. A _Markov kernel_
$P:\Omega\times\mathcal{F}\to[0,1]$ satisfies

  1. 1.

for each $x\in\Omega$, the map $P(\cdot|x):=P(x,\cdot)$ from $\mathcal{F}$ to
$[0,1]$ is a probability measure on $(\Omega,\mathcal{F})$.

  2. 2.

for each $E\in\mathcal{F}$, the map $P(\cdot,E)$ from $\Omega$ to $[0,1]$ is
$\mathcal{F}$-measurable.

Given a probability measure $\mu$ over $\Omega$ and function
$f:\Omega\to\mathbb{R}$, a Markov kernel $P$ acts on $\mu$ and $f$ to yield
objects $\mu P$ and $Pf$ defined by

$\displaystyle\mu P(\cdot)$
$\displaystyle:=\int_{\Omega}P(\cdot|x)\,\mu(\mathrm{d}x)\,,$ $\displaystyle
Pf(x)$ $\displaystyle:=\mathbb{E}_{Y\sim
P(\cdot|x)}[f(Y)]=\int_{\Omega}f(y)\,P(\mathrm{d}y|x)\,.$

A probability measure $\pi$ is called _stationary_ for $P$ if $\pi P=\pi$.

In Markov semigroup theory, it is of major importance to study the
_contractivity_ of a Markov kernel, since it quantifies the convergence rate
of its corresponding Markov chain toward the stationary distribution $\pi$.
This contractivity can be captured via the contraction coefficient of $P$
defined by

$\|P\|_{L^{p}\to L^{p}}:=\sup_{0\neq f\in
L^{p}_{0}}\frac{\|Pf\|_{L^{p}}}{\|f\|_{L^{p}}}\,,$

where $\|f\|_{L^{p}}:=\|f\|_{L^{p}(\pi)}=(\mathbb{E}_{\pi}[|f|^{p}])^{1/p}$
and $L^{p}_{0}:=\\{f:\mathbb{E}_{\pi}[|f|^{p}]<\infty,\mathbb{E}_{\pi}f=0\\}$.

The most classical setting that has been studied is the $L^{2}(\pi)$ space,
whose contraction coefficient is given by $\gamma$ if the _spectral gap_ of
the Markov kernel $P$ is $1-\gamma$. By substituting
$f=\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1$, it is also possible to quantify the
convergence rate of a Markov chain in $\chi^{2}$-divergence with reference to
the same constant $\gamma$.

This classical setting, however, is not sufficient for understanding
$\mathcal{R}_{\infty}$, since we only have the one-sided bound
$\mathcal{R}_{2}=\log(1+\chi^{2})\leq\mathcal{R}_{\infty}=\log\operatorname{ess\,sup}|\frac{\mathrm{d}\mu}{\mathrm{d}\pi}|=\log\|\frac{\mathrm{d}\mu}{\mathrm{d}\pi}\|_{L^{\infty}}$.
Instead, it is natural to study the contraction coefficient in
$L^{\infty}(\pi)$, defined by

$\|P\|_{L^{\infty}\to L^{\infty}}:=\sup_{0\neq f\in
L^{\infty}_{0}}\frac{\|Pf\|_{L^{\infty}}}{\|f\|_{L^{\infty}}}\,,$

where
$\|f\|_{L^{\infty}}:=\|f\|_{L^{\infty}(\pi)}=\operatorname{ess\,sup}_{\pi}\lvert
f\rvert$ and $L^{\infty}_{0}:=\\{f:\operatorname{ess\,sup}_{\pi}\lvert
f\rvert<\infty,\mathbb{E}_{\pi}f=0\\}$. One observes that $L^{\infty}\to
L^{\infty}$ contractivity implies the _uniform ergodicity_ of a Markov chain,
and that the opposite inequality also holds due to [Rud11, Proposition 3.23].

######  Proposition 13 ([Rud11, Proposition 3.23]).

Let $P$ be a Markov kernel that is reversible with respect to the stationary
distribution $\pi$. Then,

$\|P^{n}-1_{\pi}\|_{L^{\infty}\to L^{\infty}}\leq
2\operatorname{ess\,sup}_{x}\|P^{n}(\cdot|x)-\pi\|_{\mathsf{TV}}\,,$

where $1_{\pi}$ is the operator defined by $1_{\pi}(f)=\mathbb{E}_{\pi}[f]$.

We can now deduce an explicit $L^{\infty}$-convergence result for a reversible
Markov chain as follows.

######  Theorem 14.

Consider a Markov chain with kernel $P$, initial distribution $\mu$ and
stationary distribution $\pi$. Then,

$\bigl{\|}\frac{\mathrm{d}\mu
P^{n}}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\leq\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\cdot
2\operatorname{ess\,sup}_{x}\|P^{n}(\cdot|x)-\pi\|_{\mathsf{TV}}\,$

In particular, $\mathcal{R}_{\infty}(\mu
P^{n}\mathbin{\|}\pi)\leq\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\cdot
2\operatorname{ess\,sup}_{x}\|P^{n}(\cdot|x)-\pi\|_{\mathsf{TV}}$.

###### Proof.

Setting $f=\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1$, we have

$\|P^{n}-1_{\pi}\|_{L^{\infty}\to
L^{\infty}}\geq\frac{\|P^{n}f-1_{\pi}(f)\|_{L^{\infty}}}{\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}}=\frac{\|P^{n}f\|_{L^{\infty}}}{\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}}=\frac{\bigl{\|}\frac{\mathrm{d}\mu
P^{n}}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}}{\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}}\,,$

where the last equality follows from
$P(\frac{\mathrm{d}\mu}{\mathrm{d}\pi})=\frac{\mathrm{d}(\mu
P)}{\mathrm{d}\pi}$ due to the reversibility of $P$. Therefore,

$\bigl{\|}\frac{\mathrm{d}\mu
P^{n}}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\leq\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\cdot
2\operatorname{ess\,sup}_{x}\|P^{n}(\cdot|x)-\pi\|_{\mathsf{TV}}\,.$

The Rényi-infinity bound immediately follows from $\log(1+x)\leq x$. ∎

###  3.2 LSI to uniform ergodicity without overhead

We now need to control the $\mathsf{TV}$-distance uniformly over any initial
point $x\in\Omega$. That is, one should find the iteration number $n$ of the
Markov chain such that
$\|\delta_{x}P^{n}-\pi\|_{\mathsf{TV}}\lesssim\varepsilon$ for almost every
$x\in\Omega$, where $\delta_{x}$ denotes the Dirac measure at $x$. One expects
that it is impossible in general to bound this quantity for arbitrary Markov
chains and stationary measures, but we can get around this in our current
setting, which only considers probability measures with compact support.

Mixing rates of many convex bodies samplers have been studied when $\pi$
satisfies a Cheeger isoperimetric inequality (which is equivalent to a
Poincaré inequality for log-concave distributions). For comparison, standard
choices of Markov kernel in the unconstrained setting (such as the Langevin or
underdamped Langevin dynamics) relate Poincaré inequalities for $\pi$ to the
convergence of the sampler in $\chi^{2}$, and have theoretical guarantees
typically given by

$\chi^{2}(\delta_{x}P^{n}\mathbin{\|}\pi)\lesssim\exp\bigl{(}-\frac{n}{C_{\mathsf{PI}}(\pi)}\bigr{)}\,\chi^{2}(\delta_{x}P^{1}\mathbin{\|}\pi)\,.$

As $2\,\mathsf{TV}^{2}\leq\chi^{2}$, one can achieve
$\varepsilon$-$\mathsf{TV}$ in
$C_{\mathsf{PI}}(\pi)\log(\chi^{2}(\delta_{x}P^{1}\mathbin{\|}\pi)/\varepsilon^{2})\lesssim
C_{\mathsf{PI}}(\pi)\log(\|\frac{\mathrm{d}(\delta_{x}P^{1})}{\mathrm{d}\pi}\|_{L_{\infty}}/\varepsilon)$
iterations. However, the initial $\chi^{2}$ is typically $\exp(\Omega(d))$, so
under (PI) this approach ends up incurring an additional overhead of
$\Omega(d)$ to the mixing rate.

While it would already be impressive to boost a mixing rate from the weakest
metric ($\mathsf{TV}$) to the strongest metric ($\mathcal{R}_{\infty}$) under
(PI) at the cost of additional $\operatorname{poly}(d)$, one can achieve this
with only additional $\operatorname{polylog}(d)$ factors under (LSI). Just as
a Poincaré inequality for $\pi$ is normally sufficient to imply the
exponential convergence of a corresponding Markov process to $\pi$ in
$\chi^{2}$, the log-Sobolev inequality is equivalent to exponential
convergence of many process in entropy (or in $\mathsf{KL}$). Hence, under
(LSI), theoretical guarantees of Markov chains are typically of the form

$\mathsf{KL}(\delta_{x}P^{n}\mathbin{\|}\pi)\lesssim\exp\bigl{(}-\frac{n}{C_{\mathsf{LSI}}(\pi)}\bigr{)}\,\mathsf{KL}(\delta_{x}P^{1}\mathbin{\|}\pi)\,.$

Since we know that $2\,\mathsf{TV}^{2}\leq\mathsf{KL}$ (CKP inequality),
$\varepsilon$-$\mathsf{TV}$ can be achieved after
$C_{\mathsf{LSI}}(\pi)\log(\mathsf{KL}(\delta_{x}P^{1}\mathbin{\|}\pi)/\varepsilon)\lesssim
C_{\mathsf{LSI}}(\pi)\log(\frac{1}{\varepsilon}\log\|\frac{\mathrm{d}(\delta_{x}P^{1})}{\mathrm{d}\pi}\|_{L_{\infty}})$
iterations. As noted earlier, this initial distance is at worst
$\exp(d^{\mathcal{O}(1)})$, which only results in additional
$\operatorname{polylog}(d)$ factors after evaluating the double logarithm.
Therefore, one can boost from $\mathsf{TV}$ to $\mathcal{R}_{\infty}$ with
only logarithmic overhead through (LSI).

####  3.2.1 Rényi-infinity guarantees for uniform sampling under warm start

We demonstrate the effectiveness of this boosting technique using (LSI) for
the uniform distribution. Thereby we achieve a
$\mathcal{R}_{\infty}$-guarantee of uniform sampling without compromising the
well-known mixing rate
$\widetilde{\mathcal{O}}(d^{2}D^{2}\log\frac{1}{\varepsilon})$ of
$\mathsf{Ball\ walk}$ and $\mathsf{Hit\text{-}and\text{-}Run}$ (given in terms
of $\mathsf{TV}$ and $\chi^{2}$ respectively).

To this end, it is useful to work with a sampler whose mixing is well
understood under several functional inequalities such as (PI) and (LSI).
Recently, [KVZ24] studies $\mathsf{Proximal\ sampler}$444In the original work,
it is called the $\mathsf{In\text{-}and\text{-}Out}$, inspired by the
geometric behavior of the sampler. We call it $\mathsf{Proximal\ sampler}$
instead, since this geometric behaviour is not as clear when working with
arbitrary target distributions. for uniform distributions under those
functional inequalities, via calculations following the heat flow and its
$\pi$-adjoint. In particular, [KVZ24] already establishes
$\mathcal{R}_{q}$-guarantees (with $q<\infty$) of the $\mathsf{Proximal\
sampler}$ with respect to the uniform distribution $\pi$, given any starting
measure.

######  Lemma 15 ([KVZ24, Corollary 28]).

For any $\varepsilon\in(0,1)$ and convex body
$\mathcal{K}\subset\mathbb{R}^{d}$ with diameter $D$, let $P$ be the Markov
kernel of the $\mathsf{Proximal\ sampler}$ with variance $h$. For given
$x\in\mathcal{K}$, let $\mu_{x}^{n}:=\delta_{x}P^{n}$ be the law of the $n$-th
iterate of the $\mathsf{Proximal\ sampler}$, and $\pi$ be the uniform
distribution over $\mathcal{K}$. Then,
$\mathcal{R}_{q}(\mu_{x}^{n}\mathbin{\|}\pi)\leq\varepsilon$ for
$n=\widetilde{\mathcal{O}}\bigl{(}qh^{-1}C_{\mathsf{LSI}}(\pi)\log\frac{d+h^{-1}D^{2}}{\varepsilon}\bigr{)}$.

Combining this with the boosting technique, we obtain a
$\mathcal{R}_{\infty}$-guarantee for uniform sampling:

######  Lemma 16.

For a convex body $\mathcal{K}\subset\mathbb{R}^{d}$ with diameter $D$, assume
that an initial distribution $\mu$ is $M$-warm with respect to the uniform
distribution $\pi$ over $\mathcal{K}$. For any $\varepsilon\in(0,1)$, the
$\mathsf{Proximal\ sampler}$ with variance $h$ achieves
$\mathcal{R}_{\infty}(\mu P^{n}\mathbin{\|}\pi)\leq\varepsilon$ (or
$1-\varepsilon\leq\frac{\mathrm{d}\mu P^{n}}{\mathrm{d}\pi}\leq 1+\varepsilon$
on $\mathcal{K}$) after
$n=\widetilde{\mathcal{O}}\bigl{(}h^{-1}D^{2}\log\frac{M(d+h^{-1}D^{2})}{\varepsilon}\bigr{)}$
iterations.

###### Proof.

Using
$2\|\cdot-\pi\|_{\mathsf{TV}}^{2}\leq\mathsf{KL}(\cdot\mathbin{\|}\pi)=\lim_{q\downarrow
1}\mathcal{R}_{q}(\cdot\mathbin{\|}\pi)$ and
$C_{\mathsf{LSI}}(\pi)=\mathcal{O}(D^{2})$, one obtains that after $n\gtrsim
h^{-1}D^{2}\log\frac{M(d+h^{-1}D^{2})}{\varepsilon}$ iterations,

$\sup_{x\in\mathcal{K}}\|\mu_{x}^{n}-\pi\|_{\mathsf{TV}}\leq\frac{\varepsilon}{M}\,.$

By Theorem 14 with
$\bigl{\|}\frac{\mathrm{d}\mu}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\leq M$,
we have $\bigl{\|}\frac{\mathrm{d}\mu
P^{n}}{\mathrm{d}\pi}-1\bigr{\|}_{L^{\infty}}\leq\varepsilon$ and
$\mathcal{R}_{\infty}(\mu P^{n}\mathbin{\|}\pi)\leq\varepsilon$. ∎

######  Remark 17.

We note that the above guarantee bounds the _iteration number_ for mixing in
$\mathcal{R}_{\infty}$, not the query complexity, through any-start guarantees
in $\mathsf{TV}$. The query complexity needed to attain $\varepsilon$
$\mathsf{TV}$-distance from any start in our implementation can potentially be
exponential in $d$. This, however, will not be relevant to our results since
the kernel $P$ only captures the accepted proposals. This way, we can view the
result above as merely extracting a latent property of the algorithm, which is
not dependent on details of its implementation. For Metropolized algorithms
such as $\mathsf{Ball\ walk}$ however, the kernel $P$ will need to take the
number of rejections into account, and the dependence on $d$ for any-start
guarantees can potentially scale poorly as a result.

Lastly, combined with the query complexity of implementing each step of the
$\mathsf{Proximal\ sampler}$ for uniform distributions, we obtain a guarantee
on the query complexity for achieving $\mathcal{R}_{\infty}$-mixing for
uniform sampling.

######  Theorem 18.

For any $\eta,\varepsilon\in(0,1)$, $n\in\mathbb{N}$ defined below, and convex
body $\mathcal{K}$ given by a well-defined membership oracle, the
$\mathsf{Proximal\ sampler}$ (Algorithm 1) with
$h=(2d^{2}\log\frac{9nM}{\eta})^{-1}$,
$N=\widetilde{\mathcal{O}}(\frac{nM}{\eta})$, and initial distribution
$\mu_{0}$ $M$-warm with respect to $\pi$ the uniform distribution over
$\mathcal{K}$ achieves
$\mathcal{R}_{\infty}(\mu_{n}\mathbin{\|}\pi)\leq\varepsilon$ after
$n=\widetilde{\mathcal{O}}(d^{2}D^{2}\log^{2}\frac{M}{\eta\varepsilon})$
iterations, where $\mu_{n}$ is the law of the $n$-th iterate. With probability
$1-\eta$, the algorithm iterates this many times without failure, using
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log^{6}\frac{1}{\eta\varepsilon})$
expected number of membership queries in total.

###### Proof.

By [KVZ24, Lemma 14], if one takes variance
$h=(2d^{2}\log\frac{9nM}{\eta})^{-1}$ and threshold
$N=\widetilde{\mathcal{O}}(\frac{nM}{\eta})$, then for each iteration the
expected number of membership queries is
$\widetilde{\mathcal{O}}(M\log\frac{1}{\eta\varepsilon})$, and the failure
probability is at most $\eta/n$. By Lemma 16, the $\mathsf{Proximal\ sampler}$
should iterate
$n=\widetilde{\mathcal{O}}(d^{2}D^{2}\log\frac{M}{\eta\varepsilon})$ times to
output a sample whose law is $\varepsilon$-close to $\pi$ in
$\mathcal{R}_{\infty}$. Therefore, throughout the $n$ outer iterations, the
failure probability of the $\mathsf{Proximal\ sampler}$ is at most $\eta$ by a
union bound, and the total expected number of queries is
$\widetilde{\mathcal{O}}(Md^{2}D^{2}\log\frac{M}{\eta\varepsilon})$. ∎

##  4 Uniform ergodicity of proximal sampler for truncated Gaussians

In the $\mathsf{Gaussian\ cooling}$ algorithm [CV18], a truncated Gaussian
$\pi_{\sigma^{2}}:=\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$ for a convex
body $\mathcal{K}$ serves as an annealing distribution, where the variance
$\sigma^{2}$ is gradually increased across phases. In particular,
$\mathsf{Gaussian\ cooling}$ uses a Metropolized $\mathsf{Ball\ walk}$ to
sample such a truncated Gaussian. Its convergence rate is quantified in
$\mathsf{TV}$ through a Cheeger isoperimetry of the truncated Gaussian
(namely, a Poincaré inequality).

In order for us to properly carry the warmness across phases, we must need a
$\mathcal{R}_{\infty}$-guarantee for sampling truncated Gaussians. Just as we
established $\mathcal{R}_{\infty}$ guarantees for uniform sampling via Theorem
14 and (LSI) for the uniform distribution, it is natural to propose a sampler
for a truncated Gaussian whose convergence rate can be concisely related to
$C_{\mathsf{LSI}}(\pi_{\sigma^{2}})$. Thus, in §4.1 we analyze the
$\mathsf{Proximal\ sampler}$ for $\pi_{\sigma^{2}}$, providing its convergence
rate in terms of $C_{\mathsf{LSI}}(\pi_{\sigma^{2}})$ through a heat flow
perspective just as in [KVZ24]. Then in §4.2 we establish the uniform
ergodicity of the $\mathsf{Proximal\ sampler}$ for $\pi_{\sigma^{2}}$,
deducing the $\mathcal{R}_{\infty}$ guarantee in Theorem 18.

###  4.1 Convergence analysis

Compared with the $\mathsf{Proximal\ sampler}$ for uniform distributions, the
$\mathsf{Proximal\ sampler}$ for truncated Gaussians (Algorithm 2) requires a
different backward step (i.e., the implementation of RGO) while using the same
forward step (i.e., Gaussian step).

One iteration of the $\mathsf{Proximal\ sampler}$ for
$\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$ consists of two steps:

  * •

(Forward step) $y\sim\pi^{Y|X}(\cdot|x)=\mathcal{N}(x,hI_{d})$.

  * •

(Backward step)

$x\sim\pi^{X|Y}(\cdot|y)\propto\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}-\frac{1}{2h}\|y-x\|^{2}\bigr{)}\cdot\mathds{1}_{\mathcal{K}}(x)\propto\mathcal{N}\Bigl{(}\frac{1}{1+h\sigma^{-2}}y,\frac{h}{1+h\sigma^{-2}}I_{d}\Bigr{)}\Big{|}_{\mathcal{K}}\,.$

To implement the backward step, we use rejection sampling with the proposal
$\mathcal{N}\bigl{(}\frac{1}{1+h\sigma^{-2}}y,\frac{h}{1+h\sigma^{-2}}I_{d}\bigr{)}$,
accepting if the proposal lies inside of $\mathcal{K}$. Then the expected
number of trials for the first success is

$\frac{1}{\ell(y)}:=\frac{\bigl{(}2\pi(h^{-1}+\sigma^{-2})^{-1}\bigr{)}^{d/2}}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{\sigma^{-2}+h^{-1}}{2}\|x-h^{-1}(\sigma^{-2}+h^{-1})^{-1}y\|^{2}\bigr{)}\,\mathrm{d}x}\,.$

We can write down the density of $\pi^{Y}$ as follows:

$\displaystyle\pi^{Y}(y)$
$\displaystyle=\frac{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}-\frac{1}{2h}\|x-y\|^{2}\bigr{)}\,\mathrm{d}x}{(2\pi
h)^{d/2}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}$
$\displaystyle=\frac{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2}(\sigma^{-2}+h^{-1})\|x-h^{-1}(\sigma^{-2}+h^{-1})^{-1}y\|^{2}\bigr{)}\,\mathrm{d}x}{(2\pi
h)^{d/2}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\,\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}$
(4.1)
$\displaystyle=\frac{(1+h\sigma^{-2})^{-d/2}\,\ell(y)}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\,\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}\,.$
(4.2)

####  4.1.1 Mixing analysis

The two steps within $\mathsf{Proximal\ sampler}$ have continuous
interpolation via the forward and backward heat flow, so their mixing
guarantees can be naturally associated with functional inequalities (e.g. (PI)
and (LSI)) for a target distribution. Such a mixing result for the
$\mathsf{Proximal\ sampler}$ has already been established for unconstrained
distributions $\pi^{X}\propto\exp(-V)$ [CCSW22], for which the
$\mathsf{Proximal\ sampler}$ can be generalized as follows: for
$\pi(x,y)\propto\exp\bigl{(}-V(x)-\frac{1}{2h}\|x-y\|^{2}\bigr{)}$, repeat (i)
$y_{i+1}\sim\pi^{Y|X=x_{i}}=\mathcal{N}(x_{i},hI_{d})$ and (ii)
$x_{i+1}\sim\pi^{X|Y=y_{i+1}}$.

######  Proposition 19 ([CCSW22, Theorem 3]).

Let $\mu_{k}^{X}$ be the law of the $k$-th output of $\mathsf{Proximal\
sampler}$ with initial distribution $\mu_{0}^{X}$. 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})}{\bigl{(}1+h/C_{\mathsf{LSI}}(\pi^{X})\bigr{)}^{2k/q}}\,.$

This has been further extended to constrained distributions, including the
truncated Gaussian, under only mild additional assumptions.

######  Lemma 20 ([KVZ24, Lemma 22]).

Let $\nu$ be a measure, absolutely continuous with respect to the uniform
distribution over $\mathcal{K}$, and $\mu_{0}$ an arbitrary measure. The
forward and backward heat flow solutions 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 $q$-Rényi divergence with
$q\in(1,\infty)$,

$\mathcal{R}_{q}(\mu_{h}^{\leftarrow}\mathbin{\|}\nu)\leq\lim_{t\downarrow
0}\mathcal{R}_{q}(\mu_{h-t}^{\leftarrow}\mathbin{\|}\nu_{t})\,.$

We will set $\nu=\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$ in above. It
turns out that the solutions to the two equations above give precisely the
laws of $\mu_{k}^{Y},\mu_{k+1}^{X}$ when starting at $\mu_{k}^{X}$. Secondly,
it is well-known in [BGL14] that a truncated Gaussian $\nu$ has
$C_{\mathsf{LSI}}(\nu)\leq\sigma^{2}$, as truncation to a convex set only
improves the log-Sobolev constant. Then we can derive a contraction result of
the $\mathsf{Proximal\ sampler}$ for $\nu$ in $\mathcal{R}_{q}$ for any $q\geq
1$, emulating the proof for uniform distributions in [KVZ24, Lemma 23]. First,
convolve the truncated Gaussian with Gaussian with very small variance
$\epsilon$, which generates a smooth and unconstrained distribution. Then
invoke the contraction result in Proposition 19, and use lower semicontinuity
of the $\mathcal{R}_{q}$ divergence to conclude when sending $\varepsilon\to
0$.

######  Lemma 21.

Let $\mu^{X}_{k}$ be the law of the $k$-th output of the $\mathsf{Proximal\
sampler}$ with initial distribution $\mu_{0}^{X}$ and target
$\pi^{X}=\mathcal{N}(0,\sigma^{2}I_{d})|_{\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\sigma^{-2})^{2k/q}}\,.$

###### Proof.

For small $\epsilon>0$, as
$\mu_{\epsilon}=(\mu_{0}^{X})_{\epsilon}=\mu_{0}^{X}*\mathcal{N}(0,\epsilon
I_{d})$ and $\pi_{\epsilon}=\pi^{X}*\mathcal{N}(0,\epsilon I_{d})$ are
$C^{\infty}$-smooth, we can invoke the decay result with step size
$h-\epsilon$ in Proposition 19. Thus, for contraction constants
$C_{\epsilon}=(1+\frac{h-\epsilon}{\sigma^{2}+\epsilon})^{-2/q}$ (since
$C_{\mathsf{LSI}}\bigl{(}\nu*\mathcal{N}(0,\epsilon I_{d})\bigr{)}\leq
C_{\mathsf{LSI}}(\nu)+\epsilon$ in general), it follows that

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

where we used the data-processing inequality (Lemma 8) for the last
inequality. By the lower semicontinuity of $\mathcal{R}_{q}$ as noted earlier,
sending $\epsilon\to 0$ leads to

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

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

####  4.1.2 Per-step guarantees

We can find an _effective domain_ under $\pi^{Y}$, from which
$\mathsf{Proximal\ sampler}$ escapes with only negligible probability. More
precisely, we denote the $\delta$-blowup of $\mathcal{K}$ by
$\mathcal{K}_{\delta}=\\{x\in\mathbb{R}^{d}:d(x,\mathcal{K})\leq\delta\\}$.

######  Lemma 22.

Let $R=(1+h\sigma^{-2})\,\mathcal{K}_{\delta}$ with $\delta=t/d$ and
$h=\sigma^{2}c^{2}(d^{-1}\wedge(d^{2}\sigma^{2}-c^{2})^{-1})$, where $c$ is
any constant smaller than $d^{2}\sigma^{2}$ and $t\geq 2c(c+1)$. Then,

$\pi^{Y}(R^{c})\leq\exp\bigl{(}c^{2}-\frac{t^{2}}{8c^{2}}\bigr{)}\,.$

###### Proof.

Using the density formula for $\pi^{Y}$ in (4.1),

$\displaystyle\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x\cdot\pi^{Y}(R^{c})$
$\displaystyle=\int_{[(1+h\sigma^{-2})\,\mathcal{K}_{\delta}]^{c}}\frac{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2}(\sigma^{-2}+h^{-1})\|x-(1+h\sigma^{-2})^{-1}y\|^{2}\bigr{)}\,\mathrm{d}x}{(2\pi
h)^{d/2}}\,\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}\,\mathrm{d}y$
$\displaystyle\underset{(i)}{=}\frac{(1+h\sigma^{-2})^{d}}{(2\pi
h)^{d/2}}\int_{\mathcal{K}_{\delta}^{c}}\int_{\mathcal{K}}\exp\bigl{(}-\frac{\sigma^{-2}+h^{-1}}{2}\|x-z\|^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}(1+h\sigma^{-2})^{2}\|z\|^{2}\bigr{)}\,\mathrm{d}x\mathrm{d}z$
$\displaystyle=\frac{(1+h\sigma^{-2})^{d}}{(2\pi
h)^{d/2}}\int_{\mathcal{K}_{\delta}^{c}}\int_{\mathcal{K}}\exp\bigl{(}-\frac{\sigma^{-2}+h^{-1}}{2}\|x-z\|^{2}\bigr{)}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}x\mathrm{d}z$
$\displaystyle\underset{(ii)}{\leq}\frac{(1+h\sigma^{-2})^{d}}{(2\pi
h)^{d/2}}\int_{\mathcal{K}_{\delta}^{c}}\int_{\mathcal{H}(z)}\exp\bigl{(}-\frac{\sigma^{-2}+h^{-1}}{2}\|x-z\|^{2}\bigr{)}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}x\mathrm{d}z$
$\displaystyle=(1+h\sigma^{-2})^{d/2}\int_{\mathcal{K}_{\delta}^{c}}\int_{d(z,\mathcal{K})}^{\infty}\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\bigl{(}-\frac{(\sigma^{-2}+h^{-1})y^{2}}{2}\bigr{)}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}y\mathrm{d}z\,,$

where $(i)$ follows from the change of variables $z=(1+h\sigma^{-2})^{-1}y$,
and in $(ii)$ $\mathcal{H}(z)$ denotes the supporting half-space at
$\mathsf{proj}_{\mathcal{K}}(z)$ containing $\mathcal{K}$ for given
$z\in\partial\mathcal{K}_{\delta}$, given as

$\mathcal{H}(z)=\\{x\in\mathbb{R}^{d}:\langle\mathsf{proj}_{\mathcal{K}}(z)-z,x-\mathsf{proj}_{\mathcal{K}}(z)\rangle\geq
0\\}\,,$

when $z\not\in\mathcal{K}$.

We define the one dimensional Gaussian integral

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

By the co-area formula and integration by parts, for $\mathcal{H}^{d-1}$ the
$(d-1)$-dimensional Hausdorff measure

$\displaystyle\int_{\mathcal{K}_{\delta}^{c}}\int_{d(z,\mathcal{K})}^{\infty}\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\Bigl{(}-\frac{(\sigma^{-2}+h^{-1})y^{2}}{2}\Bigr{)}\exp\Bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\Bigr{)}\,\mathrm{d}y\mathrm{d}z$
$\displaystyle=\int_{\delta}^{\infty}\mathsf{T}(s)\int_{\partial\mathcal{K}_{s}}\exp\Bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\Bigr{)}\,\mathcal{H}^{d-1}(\mathrm{d}z)\,\mathrm{d}s$
$\displaystyle=\Bigl{[}\underbrace{\mathsf{T}(s)\int_{0}^{s}\int_{\partial\mathcal{K}_{r}}\exp\Bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\Bigr{)}\,\mathcal{H}^{d-1}(\mathrm{d}z)\,\mathrm{d}r}_{\eqqcolon\mathsf{I}}\Bigr{]}_{s=\delta}^{\infty}$
$\displaystyle\qquad+\int_{\delta}^{\infty}\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\Bigl{(}-\frac{(\sigma^{-2}+h^{-1})s^{2}}{2}\Bigr{)}\int_{0}^{s}\int_{\partial\mathcal{K}_{r}}\exp\Bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\Bigr{)}\,\mathcal{H}^{d-1}(\mathrm{d}z)\,\mathrm{d}r\mathrm{d}s\,.$

The double integral above can be bounded as follows:

$\displaystyle\int_{0}^{s}\int_{\partial\mathcal{K}_{c}}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathcal{H}^{d-1}(\mathrm{d}z)\,\mathrm{d}c$
$\displaystyle=\int_{\mathcal{K}_{s}\backslash\mathcal{K}}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z$
$\displaystyle\underset{(i)}{\leq}\int_{(1+s)\,\mathcal{K}}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z$
$\displaystyle=(1+s)^{d}\int_{\mathcal{K}}\exp\bigl{(}-\frac{(1+s)^{2}(1+h\sigma^{-2})}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z$
$\displaystyle\leq(1+s)^{d}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z\,.$

where we used $\mathcal{K}_{s}\subset(1+s)\mathcal{K}$ in $(i)$, which follows
from $B_{1}(0)\subset\mathcal{K}$. Hence, the double integral is bounded by
$(1+s)^{d}\operatorname{vol}(\mathcal{K})$. Recall a standard tail bound for a
Gaussian distribution:

$\mathsf{T}(s)\leq\frac{1}{2}\exp\Bigl{(}-\frac{1}{2}\bigl{(}s(\sigma^{-2}+h^{-1})^{1/2}-1\bigr{)}^{2}\Bigr{)}\,.$

Combining these two bounds, it holds that $\mathsf{I}$ vanishes at $s=\infty$.

Upper bounding $\mathsf{I}$ by $0$, we have just derived that

$\displaystyle\int_{\mathcal{K}}$
$\displaystyle\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x\cdot\pi^{Y}(R^{c})$
$\displaystyle\qquad\leq(1+h\sigma^{-2})^{d/2}\int_{\delta}^{\infty}(1+s)^{d}\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\bigl{(}-\frac{(\sigma^{-2}+h^{-1})s^{2}}{2}\bigr{)}\,\mathrm{d}s\cdot\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z\,.$

Dividing both sides by the factor
$\int_{\mathcal{K}}\exp(-\frac{1}{2\sigma^{2}}\|x\|^{2})\,\mathrm{d}x$, we
obtain the following bound,

$\displaystyle\pi^{Y}(R^{c})$
$\displaystyle\leq(1+h\sigma^{-2})^{d/2}\int_{\delta}^{\infty}(1+s)^{d}\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\bigl{(}-\frac{(\sigma^{-2}+h^{-1})s^{2}}{2}\bigr{)}\,\mathrm{d}s$
$\displaystyle\leq(1+h\sigma^{-2})^{d/2}\int_{\delta}^{\infty}\exp(sd)\sqrt{\frac{\sigma^{-2}+h^{-1}}{2\pi}}\exp\bigl{(}-\frac{(\sigma^{-2}+h^{-1})s^{2}}{2}\bigr{)}\,\mathrm{d}s$
$\displaystyle\underset{(i)}{\leq}\frac{1}{2}(1+h\sigma^{-2})^{d/2}\exp\bigl{(}\frac{h^{\prime}d^{2}}{2}\bigr{)}\exp\Bigl{(}-\frac{1}{2}\bigl{(}\frac{\delta}{\sqrt{h^{\prime}}}-d\sqrt{h^{\prime}}-1\bigr{)}^{2}\Bigr{)}\,,$

where in $(i)$ we again use the tail bound for a Gaussian distribution. Above,
we introduced a new variable
$h^{\prime}:=(\sigma^{-2}+h^{-1})^{-1}=\frac{h}{1+h\sigma^{-2}}$. Taking
$\delta=t/d$ and $h^{\prime}=c^{2}/d^{2}$ subject to $t\geq 2c(c+1)$, we can
make

$\exp\bigl{(}\frac{h^{\prime}d^{2}}{2}\bigr{)}\exp\Bigl{(}-\frac{1}{2}\bigl{(}\frac{\delta}{\sqrt{h^{\prime}}}-d\sqrt{h^{\prime}}-1\bigr{)}^{2}\Bigr{)}\leq\exp\bigl{(}\frac{c^{2}}{2}-\frac{t^{2}}{8c^{2}}\bigr{)}\,.$

Since $h\leq\sigma^{2}c^{2}/d$, we also have
$(1+h\sigma^{-2})^{d/2}\leq\exp(c^{2}/2)$, and the claim follows. ∎

We can now provide the per-step complexity of the $\mathsf{Proximal\ sampler}$
under a warm start.

######  Lemma 23.

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

###### Proof.

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}]$.

We now bound the expectation. For the effective domain
$R=(1+h\sigma^{-2})\mathcal{K}_{\delta}$,

$\displaystyle\int_{\mathbb{R}^{d}}\underbrace{(1-\ell)^{N}\,\mathrm{d}\pi^{Y}}_{\eqqcolon\mathsf{A}}$
$\displaystyle=\int_{R^{c}}\mathsf{A}+\int_{R\cap[\ell\geq
N^{-1}\log(3nM/\eta)]}\mathsf{A}+\int_{R\cap[\ell<N^{-1}\log(3nM/\eta)]}\mathsf{A}$
$\displaystyle\underset{(i)}{\leq}\pi^{Y}(R^{c})+\int_{[\ell\geq
N^{-1}\log(3nM/\eta)]}\exp(-\ell N)\,\mathrm{d}\pi^{Y}$
$\displaystyle\qquad+\int_{R\cap[\ell<N^{-1}\log(3nM/\eta)]}\frac{(1+h\sigma^{-2})^{-d/2}\ell(y)}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}\,\mathrm{d}y$
$\displaystyle\underset{(ii)}{\leq}e^{c_{1}^{2}-\frac{t^{2}}{8c_{1}^{2}}}+\frac{\eta}{3nM}+\frac{\log(3nM/\eta)}{N}\,\frac{\int_{R}(1+h\sigma^{-2})^{-d/2}\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}\,\mathrm{d}y}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}$
$\displaystyle\underset{(iii)}{\leq}e^{c_{1}^{2}}\exp(-\frac{t^{2}}{8c_{1}^{2}})+\frac{\eta}{3nM}+\frac{\log(3nM/\eta)}{N}\cdot(1+h\sigma^{-2})^{d/2}(1+\delta)^{d}$
$\displaystyle\underset{(iv)}{\leq}e^{c_{1}^{2}}\exp(-\frac{t^{2}}{8c_{1}^{2}})+\frac{\eta}{3nM}+\frac{\exp(t+c^{2})}{N}\,\log\frac{3nM}{\eta}\,,$

where in $(i)$ we bounded the $(1-\ell)^{N}\leq 1$ in the first term,
$(1-\ell)^{N}\leq\exp(-\ell N)$ in the second term, and again
$(1-\ell)^{N}\leq 1$ in the third term, and used the density formula (4.2) of
$\pi^{Y}$ in third term. In $(ii)$, the first bound follows from Lemma 22,
while the second and third uses the condition on $\ell$ over each domain.
$(iii)$ follows from the change of variables. Lastly, $(iv)$ follows from the
setup $\delta=t/d$ and
$h=\sigma^{2}c^{2}(d^{-1}\wedge(d^{2}\sigma^{2}-c^{2})^{-1})\leq\sigma^{2}c^{2}d^{-1}$
in Lemma 22.

With $c^{2}=\frac{\log\log Z}{2\log Z}$, $t=\sqrt{8}\log\log Z$, and $N=Z(\log
Z)^{4}$, the last line is bounded by $\frac{\eta}{nM}$. Therefore,

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

We now bound the expected number of trials per iteration. Let $S$ 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}}[S]$ 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_{R}\frac{1}{\ell}\,\mathrm{d}\pi^{Y}+N\pi^{Y}(R^{c})=\frac{\int_{R}(1+h\sigma^{-2})^{-d/2}\exp\bigl{(}-\frac{1}{2(h+\sigma^{2})}\|y\|^{2}\bigr{)}\,\mathrm{d}y}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}+N\pi^{Y}(R^{c})$
$\displaystyle\leq\exp(t+c_{1}^{2})+N\exp(-\Omega(t^{2}))\leq e(\log
Z)^{3}+3(\log
Z)^{4}=\mathcal{O}\Bigl{(}\bigl{(}\log\frac{nM}{\eta}\bigr{)}^{4}\Bigr{)}\,.$

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

####  4.1.3 Query complexity under a warm start

We now put together the mixing and per-step guarantees established above.

######  Proposition 24.

For any $\eta,\varepsilon\in(0,1)$, $q\geq 1$, $n\in\mathbb{N}$ defined below
and convex body $\mathcal{K}$ given by a well-defined membership oracle, the
$\mathsf{Proximal\ sampler}$ (Algorithm 2) with
$h=\sigma^{2}(\frac{1}{d}\wedge\frac{1}{d^{2}\sigma^{2}-1})(\log\frac{9nM}{\eta})^{-1}$,
$N=\widetilde{\mathcal{O}}(\frac{nM}{\eta})$, and initial distribution
$\mu_{0}^{X}$ which is $M$-warm with respect to
$\pi^{X}=\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$ achieves
$\mathcal{R}_{q}(\mu^{X}_{n}\mathbin{\|}\pi^{X})\leq\varepsilon$ after
$n=\widetilde{\mathcal{O}}\bigl{(}q(d\vee
d^{2}\sigma^{2})\log\frac{M}{\eta\varepsilon}\bigr{)}$ iterations, where
$\mu_{n}^{X}$ is the law of the $n$-the iterate. With probability $1-\eta$,
the algorithm iterates $n$ times without failure, using
$\widetilde{\mathcal{O}}\bigl{(}qM(d\vee
d^{2}\sigma^{2})(\log\nicefrac{{1}}{{\eta\varepsilon}})^{5}\bigr{)}$ expected
number of membership queries in total. In particular, the query complexity is
$\widetilde{\mathcal{O}}(qMd^{2}\sigma^{2}(\log\nicefrac{{1}}{{\eta\varepsilon}})^{5})$
when $d^{-1}\lesssim\sigma^{2}$.

###### Proof.

By Lemma 21, the $\mathsf{Proximal\ sampler}$ should iterate
$\mathcal{O}\bigl{(}q\sigma^{2}h^{-1}\log\frac{\log M}{\varepsilon}\bigr{)}$
times to achieve $\varepsilon$-distance in $\mathcal{R}_{q}$. To ensure that
the query complexity is bounded, we choose

$h=\sigma^{2}\frac{\log\log Z}{\log
Z}\bigl{(}\frac{1}{d}\wedge\frac{1}{d^{2}\sigma^{2}-1}\bigr{)}\quad\text{and}\quad
N=\frac{9nM}{\eta}\bigl{(}\log\frac{9nM}{\eta}\bigr{)}^{4}\,.$

By Lemma 23, we need the following total number of membership queries in
expectation:

$\widetilde{\mathcal{O}}\Bigl{(}qM(d\vee
d^{2}\sigma^{2})\bigl{(}\log\frac{1}{\eta\varepsilon}\bigr{)}^{5}\Bigr{)}\,.$

Hence, if $d^{-1}\lesssim\sigma^{2}$, then the query complexity is simply
$\widetilde{\mathcal{O}}(qMd^{2}\sigma^{2}\log^{5}\nicefrac{{1}}{{\eta\varepsilon}})$.
∎

###  4.2 Rényi-infinity guarantees for truncated Gaussians under warm start

To use the boosting technique, we need the uniform ergodicity of the
$\mathsf{Proximal\ sampler}$, bounding
$\|\delta_{x}P^{n}-\pi^{X}\|_{\mathsf{TV}}$ uniformly over $\mathcal{K}$,
where $P$ is the Markov kernel of the $\mathsf{Proximal\ sampler}$ for a
truncated Gaussian. To this end, for any $x\in\mathcal{K}$, we bound
$\mathcal{R}_{\infty}(\delta_{x}P^{1}\mathbin{\|}\pi^{X})$ _uniformly_ by
$\exp(\text{poly}(D,d))$, so that
$\log\log\mathcal{R}_{\infty}(\delta_{x}P^{1}\mathbin{\|}\pi^{X})$ does not
add more than a polylogarithmic factors to the complexity. In other words, we
establish a Gaussian analogue of Lemma 15.

######  Lemma 25.

For any given $\varepsilon\in(0,1)$, the $\mathsf{Proximal\ sampler}$ for a
truncated Gaussian $\pi^{X}$ with variance $h$ and any feasible start
$x_{0}\in\mathcal{K}$ achieves
$\mathcal{R}_{\infty}(\mu_{n}^{X}\mathbin{\|}\pi^{X})\leq\varepsilon$ for
$n=\widetilde{\mathcal{O}}(h^{-1}\sigma^{2}\log\frac{d+h^{-1}D^{2}}{\varepsilon})$
iterations.

###### Proof.

We bound the warmness of $\mu_{1}^{X}$ towards $\pi^{X}$ when
$\mu_{0}^{X}=\delta_{x_{0}}$. One can readily check that

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

and we should compare this with

$\pi^{X}(x)=\frac{\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\cdot\mathds{1}_{\mathcal{K}}(x)}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z}\,.$

For $D=\text{diam}(\mathcal{K})$,

$\displaystyle\exp\Bigl{(}-\frac{1+h\sigma^{-2}}{2h}\bigl{\|}z-\frac{1}{1+h\sigma^{-2}}y\bigr{\|}^{2}\Bigr{)}$
$\displaystyle=\exp\bigl{(}-\frac{1}{2h}\|z\|^{2}-\frac{1}{2\sigma^{2}}\|z\|^{2}-\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}+\frac{z^{\mathsf{T}}y}{h}\bigr{)}$
$\displaystyle\geq\exp\bigl{(}-\frac{D^{2}}{2h}\bigr{)}\cdot\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}-\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}+\frac{z^{\mathsf{T}}y}{h}\bigr{)}\,.$

As $|z^{\mathsf{T}}y|\leq\frac{1}{2}(2\|z\|^{2}+\frac{1}{2}\|y\|^{2})$ due to
Young’s inequality,

$\displaystyle\exp\bigl{(}-\frac{D^{2}}{2h}\bigr{)}\cdot\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}-\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}+\frac{z^{\mathsf{T}}y}{h}\bigr{)}$
$\displaystyle\geq$
$\displaystyle\exp\bigl{(}-\frac{D^{2}}{2h}\bigr{)}\cdot\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}-\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}-\frac{\|z\|^{2}}{h}-\frac{\|y\|^{2}}{4h}\bigr{)}$
$\displaystyle\geq$
$\displaystyle\exp\bigl{(}-\frac{3D^{2}}{2h}\bigr{)}\cdot\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}-\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}-\frac{\|y\|^{2}}{4h}\bigr{)}\,.$

Hence,

$\displaystyle\int\frac{\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2h}\bigl{\|}x-\frac{1}{1+h\sigma^{-2}}y\bigr{\|}^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2h}\|y-x_{0}\|^{2}\bigr{)}}{(2\pi
h)^{d/2}\int_{\mathcal{K}}\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2h}\bigl{\|}z-\frac{1}{1+h\sigma^{-2}}y\bigr{\|}^{2}\bigr{)}\,\mathrm{d}z}\,\mathrm{d}y$
$\displaystyle\leq$
$\displaystyle\underbrace{\frac{\exp\bigl{(}\frac{3D^{2}}{2h}\bigr{)}}{(2\pi
h)^{d/2}}\int\frac{\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2h}\bigl{\|}x-\frac{1}{1+h\sigma^{-2}}y\bigr{\|}^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2h}\|y-x_{0}\|^{2}\bigr{)}\exp\bigl{(}\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}+\frac{\|y\|^{2}}{4h}\bigr{)}}{\int_{\mathcal{K}}\exp\bigl{(}-\frac{1}{2\sigma^{2}}\|z\|^{2}\bigr{)}\,\mathrm{d}z}\,\mathrm{d}y}_{\eqqcolon(\\#)}\,.$

The numerator of the integrand can be bounded as follows:

$\displaystyle\exp\bigl{(}-\frac{1+h\sigma^{-2}}{2h}\bigl{\|}x-\frac{1}{1+h\sigma^{-2}}y\bigr{\|}^{2}\bigr{)}\exp\bigl{(}-\frac{1}{2h}\|y-x_{0}\|^{2}\bigr{)}\exp\bigl{(}\frac{1}{2h(1+h\sigma^{-2})}\|y\|^{2}+\frac{\|y\|^{2}}{4h}\bigr{)}$
$\displaystyle=$
$\displaystyle\exp\Bigl{(}-\frac{1}{4h}\|y-2(x_{0}+x)\|^{2}-\frac{1}{2h}\bigl{(}\|x_{0}\|^{2}+\|x\|^{2}+\frac{h}{\sigma^{2}}\|x\|^{2}-2\|x_{0}+x\|^{2}\bigr{)}\Bigr{)}$
$\displaystyle\leq$
$\displaystyle\exp\bigl{(}\frac{D^{2}}{h}\bigr{)}\exp\bigl{(}-\frac{1}{4h}\|y-2(x_{0}+x)\|^{2}-\frac{1}{2\sigma^{2}}\|x\|^{2}\bigr{)}\,.$

Putting this bound back to the integral above,

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

Therefore, the ratio can be bounded by

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

so $M=\operatorname{ess\,sup}\frac{\mu_{1}^{X}}{\pi^{X}}\leq
2^{d/2}\exp(5h^{-1}D^{2})$, and
$\mathcal{R}_{q}(\mu_{1}^{X}\mathbin{\|}\pi^{X})\lesssim\frac{q}{q-1}\log
M\leq\frac{q}{q-1}\bigl{(}d+h^{-1}D^{2}\bigr{)}$. By Lemma 21, one can achieve
$\mathcal{R}_{q}(\mu_{n}^{X}\mathbin{\|}\pi^{X})\leq\varepsilon$ for
$n=\widetilde{\mathcal{O}}(qh^{-1}\sigma^{2}\log\frac{d+h^{-1}D^{2}}{\varepsilon})$.

Using $2\|\cdot\|_{\mathsf{TV}}^{2}\leq\mathsf{KL}=\lim_{q\downarrow
1}\mathcal{R}_{q}$ and $C_{\mathsf{LSI}}(\pi^{X})\leq\sigma^{2}$, one can
achieve that after $n\gtrsim
h^{-1}\sigma^{2}\log\frac{M(d+h^{-1}D^{2})}{\varepsilon}$ iterations,

$\sup_{x\in\mathcal{K}}\|\mu^{X}_{n}-\pi^{X}\|_{\mathsf{TV}}\leq\frac{\varepsilon}{M}\,.$

By Theorem 14 with
$\bigl{\|}\frac{\mathrm{d}\mu^{X}_{1}}{\mathrm{d}\pi^{X}}-1\bigr{\|}_{L_{\infty}}\leq
M$, we have
$\bigl{\|}\frac{\mathrm{d}\mu^{X}_{n}}{\mathrm{d}\pi^{X}}-1\bigr{\|}_{L_{\infty}}\leq\varepsilon$
and $\mathcal{R}_{\infty}(\mu^{X}_{n}\mathbin{\|}\pi)\leq\varepsilon$. ∎

Using the uniform ergodicity above, we can now obtain a guarantee in
$\mathcal{R}_{\infty}$ of the $\mathsf{Proximal\ sampler}$ for a truncated
Gaussian over $\mathcal{K}$, boosting the metric in Proposition 24.

######  Theorem 26.

For any $\eta,\varepsilon\in(0,1)$, $n\in\mathbb{N}$ defined below, and convex
body $\mathcal{K}$ given by a well-defined membership oracle, the
$\mathsf{Proximal\ sampler}$ (Algorithm 2) with
$h=\sigma^{2}\bigl{(}\frac{1}{d}\wedge\frac{1}{d^{2}\sigma^{2}-1}\bigr{)}\bigl{(}\log\frac{9nM}{\eta}\bigr{)}^{-1}$,
$N=\widetilde{\mathcal{O}}(\frac{nM}{\eta})$, and initial distribution
$\mu_{0}^{X}$ that is $M$-warm with respect to $\pi^{X}$ the truncated
Gaussian $\mathcal{N}(0,\sigma^{2}I_{d})|_{\mathcal{K}}$ achieves
$\mathcal{R}_{\infty}(\mu^{X}_{n}\mathbin{\|}\pi^{X})\leq\varepsilon$ after
$n=\widetilde{\mathcal{O}}\bigl{(}(d\vee
d^{2}\sigma^{2})\bigl{(}\log\frac{MD}{\eta\varepsilon}\bigr{)}^{2}\bigr{)}$
iterations, where $\mu_{n}^{X}$ is the law of the $n$-the iterate. With
probability $1-\eta$, the algorithm iterates $n$ times without failure, using
$\widetilde{\mathcal{O}}\bigl{(}M(d\vee
d^{2}\sigma^{2})(\log\nicefrac{{D}}{{\eta\varepsilon}})^{6}\bigr{)}$ expected
number of membership queries in total. In particular, the query complexity is
$\widetilde{\mathcal{O}}(Md^{2}\sigma^{2}(\log\nicefrac{{D}}{{\eta\varepsilon}})^{6})$
when $d^{-1}\lesssim\sigma^{2}$.

###### Proof.

By Lemma 25, the $\mathsf{Proximal\ sampler}$ should iterate
$n=\widetilde{\mathcal{O}}\bigl{(}(d\vee
d^{2}\sigma^{2})(\log\frac{MD}{\eta\varepsilon})^{2}\bigr{)}$ times to output
a sample whose law is $\varepsilon$-close to $\pi^{X}$ in
$\mathcal{R}_{\infty}$. As argued in Proposition 24, we can conclude that
throughout $n$ outer iterations, the failure probability of the
$\mathsf{Proximal\ sampler}$ is at most $\eta$ by the union bound, and the
total expected number of queries is $\widetilde{\mathcal{O}}\bigl{(}M(d\vee
d^{2}\sigma^{2})(\log\frac{D}{\eta\varepsilon})^{6}\bigr{)}$. ∎

######  Remark 27.

It can be checked that, conditioning on the event that the algorithm has not
failed, the distribution of the iterates remains the correct distribution
$\mu_{n}^{X}$. Thus, in practical terms, the $1-\eta$ probability of failure
will not be a significant obstacle, since one can just use the samples of the
successful trials without compromising any of the guarantees.

##  5 Successive proximal scheme: Rényi infinity guarantee for uniform
sampling

We put together the ingredients prepared in previous sections, namely the
$\mathsf{Proximal\ sampler}$ for uniform distributions (Theorem 18) and
truncated Gaussian (Theorem 26), along with the annealing scheme introduced in
[CV18], and obtain $\mathsf{Proximal\ Gaussian\ cooling}$ (Algorithm 3).

In §5.1, we outline $\mathsf{Proximal\ Gaussian\ cooling}$, together with our
main result showing that it can sample a point approximately uniformly from a
well-rounded convex body with query complexity
$\widetilde{\mathcal{O}}(d^{3}\operatorname{polylog}(1/\varepsilon))$. In
§5.2, we then provide a proof for the main claim.

###  5.1 Rényi infinity guarantee with cubic complexity

As mentioned earlier in §1, using the rounding algorithm in [JLLV21], we can
assume that a given convex body $\mathcal{K}\subset\mathbb{R}^{d}$, presented
by a well-defined membership oracle, is _well-rounded_. This means that
$\mathcal{K}$ satisfies $\mathbb{E}_{X\sim\mathcal{K}}[\|X\|^{2}]\leq C^{2}d$
with a known constant $C>0$, where $X\sim\mathcal{K}$ indicates that $X$ is
drawn from the uniform distribution over $\mathcal{K}$. We state the main
result of this paper below.

######  Theorem 28.

Assume that a well-rounded convex body $\mathcal{K}$ with
$\mathbb{E}_{\mathcal{K}}[\|X\|^{2}]\leq C^{2}d$ and $0\in\mathcal{K}$ is
presented by a well-defined membership oracle. Let $\pi$ be the uniform
distribution over $\mathcal{K}$. For given $\eta,\varepsilon\in(0,1)$,
$\mathsf{Proximal\ Gaussian\ cooling}$ is a randomized algorithm, which
succeeds with probability $1-\eta$ in outputting a sample $X\sim\nu$ such that
$\mathcal{R}_{\infty}(\nu\mathbin{\|}\pi)\leq\varepsilon$. Conditioned on its
success, $\mathsf{Proximal\ Gaussian\ cooling}$ uses
$\widetilde{\mathcal{O}}(C^{2}d^{3}\log^{8}\frac{1}{\eta\varepsilon})$
membership queries in expectation.

##### Preliminaries.

We first collect a series of observations that simplify our setup and
arguments. For $\bar{X}:=\mathbb{E}_{\mathcal{K}}[X]$, by Jensen’s inequality,

$\|\bar{X}\|\leq\mathbb{E}_{\mathcal{K}}[\|X\|]\leq\sqrt{\mathbb{E}_{\mathcal{K}}[\|X\|^{2}]}\leq
C\sqrt{d}\,.$

Also, for the covariance matrix $\Sigma_{\mathcal{K}}$ of the uniform
distribution over $\mathcal{K}$,

$\operatorname{tr}(\Sigma_{\mathcal{K}})=\mathbb{E}_{\mathcal{K}}[\|X-\bar{X}\|^{2}]\leq\mathbb{E}_{\mathcal{K}}[\|X\|^{2}]\leq
C^{2}d\,.$

Lastly, it is known from [LV07, Theorem 5.17] that
$\mathbb{P}_{\mathcal{K}}\bigl{(}\|X-\bar{X}\|\geq
t\cdot\bigl{(}\operatorname{tr}(\Sigma_{\mathcal{K}})\bigr{)}^{1/2}\bigr{)}\leq\exp(-t+1)$.
Hence,

$\mathbb{P}_{\mathcal{K}}\bigl{(}\|X-\bar{X}\|\geq t\cdot
C\sqrt{d}\bigr{)}\leq\exp(-t+1)\,.$

Therefore, we can actually work with a ‘truncated’ convex body instead of the
full convex body $\mathcal{K}$:

######  Proposition 29.

There exists a constant $L=C\log\frac{3e}{\varepsilon}$ such that the volume
of $\bar{\mathcal{K}}:=\mathcal{K}\cap B_{L\sqrt{d}}(0)$ is at least
$(1-\varepsilon/3)\operatorname{vol}(\mathcal{K})$.

##### Sketch of $\mathsf{Proximal\ Gaussian\ cooling}$.

Let $\mu_{i}$ be a truncated Gaussian
$\mathcal{N}(0,\sigma_{i}^{2}I_{d})|_{\mathcal{K}}$, and $\bar{\mu}_{i}$ be an
approximate measure close to $\mu_{i}$ produced by the algorithm. Then,
$\mathsf{Proximal\ Gaussian\ cooling}$ consists of four phases:

  * •

Phase I $(\sigma^{2}=1/d$)

    * –

Initial distribution: Uniform measure over $B_{1}(0)$ denoted by
$\mu_{0}=\bar{\mu}_{0}$.

    * –

Target distribution:
$\mathcal{N}(0,\sigma^{2}I_{d})|_{\bar{\mathcal{K}}}=\mu_{1}$.

  * •

Phase II ($1/d\leq\sigma^{2}\leq 1$)

    * –

Run $\mathsf{Proximal\ sampler}$ with initial dist. $\bar{\mu}_{i}$ (not
$\mu_{i}$) and target dist. $\mu_{i+1}$.

    * –

Update $\sigma_{i+1}^{2}=\sigma_{i}^{2}\bigl{(}1+\frac{1}{d}\bigr{)}$.

  * •

Phase III ($1\leq\sigma^{2}\leq L^{2}d$)

    * –

Run $\mathsf{Proximal\ sampler}$ with initial dist. $\bar{\mu}_{i}$ (not
$\mu_{i}$) and target dist. $\mu_{i+1}$.

    * –

Update
$\sigma_{i+1}^{2}=\sigma_{i}^{2}\bigl{(}1+\frac{\sigma_{i}^{2}}{L^{2}d}\bigr{)}$.

  * •

Phase IV ($\sigma^{2}=L^{2}d$)

    * –

Run $\mathsf{Proximal\ sampler}$ with initial distribution
$\mathcal{N}(0,\sigma^{2}I)|_{\bar{\mathcal{K}}}$ and target distribution
$\pi_{\mathcal{\bar{K}}}$, the uniform distribution over $\bar{\mathcal{K}}$.

######  Remark 30.

Phases I-III can be viewed as “preprocessing steps” whose purpose is to
generate a warm start for Phase IV. We also note that, while
$\mathcal{R}_{\infty}(\nu\mathbin{\|}\pi)\leq\varepsilon$, we do not have the
(slightly stronger) property that, for $\varepsilon<1$,

$\operatorname{ess\,sup}_{\pi}\Big{|}\frac{\mathrm{d}\nu}{\mathrm{d}\pi}-1\Bigr{|}\lesssim\varepsilon\,.$

This is because the truncation of $\pi=\pi_{\mathcal{K}}$ to
$\bar{\pi}=\pi_{\bar{\mathcal{K}}}$ causes there to be
$\mathcal{A}=\mathcal{K}\backslash\bar{\mathcal{K}}$ where
$\bar{\pi}(\mathcal{A})=0<\pi(\mathcal{A})$. Since $\nu\ll\bar{\pi}$, we
cannot establish a lower bound on $\frac{\mathrm{d}\nu}{\mathrm{d}\pi}-1$
better than $0$. On the other hand, Corollary 31 shows that we can bound

$\Bigl{\|}\frac{\mathrm{d}\nu}{\mathrm{d}\pi}-1\Bigr{\|}_{L^{p}(\pi)}^{p}\leq\varepsilon\,.$

######  Corollary 31.

Under the same assumptions as Theorem 28, $\mathsf{Proximal\ Gaussian\
cooling}$ succeeds with probability $1-\eta$ in outputting a sample $X\sim\nu$
such that
$\|\frac{\mathrm{d}\nu}{\mathrm{d}\pi}-1\|_{L^{p}(\pi)}^{p}\leq\varepsilon$,
so long as $\varepsilon<1$, $p\geq 1$. Conditioned on its success,
$\mathsf{Proximal\ Gaussian\ cooling}$ uses
$\widetilde{\mathcal{O}}(C^{2}d^{3}\log^{8}\frac{1}{\eta\varepsilon})$
membership queries in expectation.

###### Proof.

As an intermediate step in the proof of Theorem 28, we obtain that
$\operatorname{ess\,sup}_{\bar{\pi}}\bigl{|}\frac{\mathrm{d}\nu}{\mathrm{d}\bar{\pi}}-1\bigr{|}\leq\frac{\varepsilon}{3}$,
so that, replacing $\bar{\mathcal{K}}$ with $\mathcal{K}$,

$\displaystyle\mathbb{E}_{\pi}\Bigl{[}\Bigl{|}\frac{\mathrm{d}\nu}{\mathrm{d}\pi}-1\Bigr{|}^{p}\Bigr{]}$
$\displaystyle\leq\frac{\operatorname{vol}(\mathcal{K}\backslash\bar{\mathcal{K}})}{\operatorname{vol}(\mathcal{K})}+\frac{\operatorname{vol}(\bar{\mathcal{K}})}{\operatorname{vol}(\mathcal{K})}\,\Bigl{|}\frac{\mathrm{d}\nu}{\mathrm{d}\bar{\pi}}\,\frac{\operatorname{vol}(\mathcal{K})}{\operatorname{vol}(\bar{\mathcal{K}})}-1\Bigr{|}^{p}$
$\displaystyle\leq\frac{\varepsilon}{3}+\bigl{(}1-\frac{\varepsilon}{3}\bigr{)}\Bigl{(}\frac{2\varepsilon}{3(1-\frac{\varepsilon}{3})}\Bigr{)}^{p}\leq\varepsilon\,,$

so long as $\varepsilon<1$ and $p\geq 1$. ∎

###  5.2 Proof details

The following subsections state a series of lemmas, each one of which proves a
guarantee on each phase of the algorithm.

##### Failure probability.

For a target failure probability $\eta$ of the entire algorithm, we can
achieve this by setting the failure probability at each phase to be on the
order of

$\hat{\eta}:=\widetilde{\mathcal{O}}\Bigl{(}\frac{\eta}{C^{2}d\log^{2}\frac{1}{\varepsilon}}\Bigr{)}\,,$
(5.1)

since the total number of inner phases is
$\widetilde{\mathcal{O}}(C^{2}d\log^{2}\frac{1}{\varepsilon})$. This will be
sufficient for the failure probability to be at most $\eta$, where we apply a
union bound over all the phases. Since the dependence on $\hat{\eta}$ of the
$\mathsf{Proximal\ sampler}$ is polylogarithmic, this will not impact the
resulting oracle complexity bounds by more than polylogarithmic factors.

##### Structure of lemmas.

We provide a lemma for each outer phase, where each lemma establishes three
facts: (i) the number of updates (to $\sigma$), (ii) a quantitative guarantee
on the warmness on each update, and (iii) the final query-complexity bound
given by Theorem 26.

Briefly, this theorem says that, starting from an $M$-warm distribution
towards a truncated Gaussian
$\mathcal{N}(0,\sigma^{2}I)|_{\bar{\mathcal{K}}}$, the oracle complexity to
achieve $\varepsilon$ error in $\mathcal{R}_{\infty}$, with success
probability $1-\hat{\eta}$, is bounded by a quantity of asymptotic order

$\widetilde{\mathcal{O}}\bigl{(}Md^{2}\sigma^{2}\log^{6}\frac{L}{\varepsilon\hat{\eta}}\bigr{)}\,.$

Since $\sigma^{2}>d^{-1}$ at all times, it suffices to choose
$h^{-1}=\widetilde{\Theta}(d^{2}\log\frac{M\log L}{\varepsilon\hat{\eta}})$
and the number of $\mathsf{Proximal\ sampler}$ iterations to be
$n=\widetilde{\mathcal{O}}(d^{2}\sigma^{2}\log^{2}\frac{L}{\varepsilon\hat{\eta}})$.

Except for the last phase, we take $\varepsilon=\log 2$, in which case the
query complexity is simply

$\widetilde{\mathcal{O}}\bigl{(}Md^{2}\sigma^{2}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,,$
(5.2)

with $h^{-1}=\widetilde{\Theta}(d^{2}\log\frac{M\log L}{\hat{\eta}})$ and
number of $\mathsf{Proximal\ sampler}$ iterations being
$n=\widetilde{\mathcal{O}}(d^{2}\sigma^{2}\log^{2}\frac{ML}{\hat{\eta}})$.
Note that $\log 2$ error in $\mathcal{R}_{\infty}$ implies that the law of the
resultant sample is at least $2$-warm with respect to the target.

####  5.2.1 Phase I

######  Lemma 32 (Phase I).

With probability at least $1-\hat{\eta}$, initial distribution
$\mu_{0}=\mathsf{Unif}(B_{1}(0))$, and target distribution
$\mu_{1}=\mathcal{N}(0,\frac{1}{d}I_{d})|_{\bar{\mathcal{K}}}$, _Phase I_
outputs a sample with law $\bar{\mu}_{1}$ satisfying

$\operatorname{ess\,sup}_{\mu_{1}}\Bigl{|}\frac{\mathrm{d}\bar{\mu}_{1}}{\mathrm{d}\mu_{1}}-1\Bigr{|}\leq
1\,.$

Its oracle complexity is

$\widetilde{\mathcal{O}}\bigl{(}d^{3/2}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,,$

using $\mathsf{Proximal\ sampler}$ with step size
$h^{-1}=\widetilde{\Theta}(d^{2}\log\frac{\log L}{\hat{\eta}})$ and
$\widetilde{\mathcal{O}}(d\log^{2}\frac{L}{\hat{\eta}})$ iterations.

###### Proof.

# Inner phases: The number of inner phases is clearly $1$.

Warmness: The initial measure $\mu_{0}$ has warmness given by

$\displaystyle\frac{\mathrm{d}\mu_{0}(x)}{\mathrm{d}\mu_{1}(x)}$
$\displaystyle=\frac{1/\operatorname{vol}(B_{1}(0))}{\exp(-\frac{d}{2}\|x\|^{2})/\int_{\bar{\mathcal{K}}}\exp(-\frac{d}{2}\|z\|^{2})\,\mathrm{d}z}$
$\displaystyle\leq\frac{(2\pi/d)^{d/2}\exp(\frac{d}{2})}{\operatorname{vol}(B_{1})}\cdot\underbrace{\frac{1}{(2\pi/d)^{d/2}}\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{d}{2}\|z\|^{2}\bigr{)}\,\mathrm{d}z}_{\leq
1}$
$\displaystyle\underset{(i)}{\leq}(2/d)^{d/2}\exp\bigl{(}\frac{d}{2}\bigr{)}\Gamma\bigl{(}\frac{d}{2}+1\bigr{)}\underset{(ii)}{\leq}\gamma\sqrt{d}\exp\bigl{(}-\frac{d}{2}\log
d+\frac{d}{2}+\frac{d}{2}\log d-\frac{d}{2}\bigr{)}\leq\gamma\sqrt{d}\,,$

where in $(i)$ we used
$\operatorname{vol}(B_{1})=\pi^{d/2}/\Gamma({\textstyle\frac{d}{2}}+1)$, and
$(ii)$ follows from Stirling’s approximation that

$\Gamma\bigl{(}\frac{d}{2}+1\bigr{)}\leq\gamma
ed^{\frac{d}{2}+1/2}(2e)^{-\frac{d}{2}}\quad\text{for some universal constant
}\gamma>0\,.$

Final complexity: It follows from substituting $M=\mathcal{O}(\sqrt{d})$ and
$\sigma^{2}=d^{-1}$ into (5.2). ∎

####  5.2.2 Phase II

######  Lemma 33 (Phase II).

With probability at least $1-\hat{\eta}i_{*}$ for
$i_{*}=\widetilde{\mathcal{O}}(d)$, initial distribution $\bar{\mu}_{1}$
(given by from Lemma 32), and target distribution
$\mu_{i_{*}}=\mathcal{N}(0,I_{d})|_{\bar{\mathcal{K}}}$, _Phase II_ outputs a
sample with law $\bar{\mu}_{i_{*}}$ satisfying

$\operatorname{ess\,sup}_{\mu_{i_{*}}}\Bigl{|}\frac{\mathrm{d}\bar{\mu}_{i_{*}}}{\mathrm{d}\mu_{i_{*}}}-1\Bigr{|}\leq
1\,,$

Its oracle complexity is

$\widetilde{\mathcal{O}}\bigl{(}d^{3}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,,$

using $\mathsf{Proximal\ sampler}$. The step size scheme for each inner phase
is presented in the proof.

###### Proof.

Note that we slightly modify the scheme, which is to take
$\sigma_{i+1}^{2}=\min(1,\sigma_{i}^{2}(1+\frac{1}{d}))$.

# Inner phases: As $(1+d^{-1})^{d}\geq 2$ for any $d\geq 1$, it takes at most
$d$ iterations to double $\sigma$. Thus, since we need to double $\sigma$ on
the order of $\mathcal{O}(\log d)$ many times, the number of inner phases
within Phase II is at most $i_{*}=\widetilde{\mathcal{O}}(d)$.

Warmness: By the construction of the algorithm, the following bound holds for
all $x\in\bar{\mathcal{K}}$:

$\displaystyle\frac{\mathrm{d}\mu_{i}(x)}{\mathrm{d}\mu_{i+1}(x)}$
$\displaystyle=\frac{\exp\bigl{(}-\frac{1}{2\sigma_{i}^{2}}\|x\|^{2}\bigr{)}}{\exp\bigl{(}-\frac{1}{2\sigma_{i+1}^{2}}\|x\|^{2}\bigr{)}}\,\frac{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i+1}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\leq\frac{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i+1}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}$
$\displaystyle\leq\frac{\bigl{(}1+\frac{1}{d}\bigr{)}^{d/2}\int_{(1+\frac{1}{d})^{-1/2}\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{i}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\leq\sqrt{e}\,.$

Furthermore, if at each previous phase, we have

$\operatorname{ess\,sup}_{\mu_{i}}\Bigl{|}\frac{\mathrm{d}\bar{\mu_{i}}}{\mathrm{d}\mu_{i}}-1\Bigr{|}\leq
1\,,$ (5.3)

then this implies warmness of constant order between $\bar{\mu}_{i}$ and
$\mu_{i+1}$:

$\operatorname{ess\,sup}_{\mu_{i+1}}\Bigl{|}\frac{\mathrm{d}\bar{\mu_{i}}}{\mathrm{d}\mu_{i+1}}\Bigr{|}\leq
2\sqrt{e}\,.$

Final complexity: For any fixed $i$, the complexity in each inner phase is
given by substituting $M\leq 2\sqrt{e}$ and $\sigma^{2}_{i}\leq 1$ into (5.2),
to obtain a complexity of

$\widetilde{\mathcal{O}}\bigl{(}d^{2}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,.$

Here, we run the $\mathsf{Proximal\ sampler}$ with
$h_{i}^{-1}=\widetilde{\Theta}(d^{2}\log\frac{\log L}{\hat{\eta}})$ and
$n_{i}=\widetilde{\mathcal{O}}(d^{2}\sigma^{2}\log^{2}\frac{L}{\hat{\eta}})$,
to get $\mathcal{R}_{\infty}$ bounded by $\log 2$. Implicitly, this verifies
the condition (5.3).

Secondly, as $i_{*}=\widetilde{\mathcal{O}}(d)$, the total query complexity is
bounded by the product of the worst-case complexity in each inner phase with
the total number of inner phases, which is given by

$\widetilde{\mathcal{O}}\bigl{(}d^{3}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,.\qed$

####  5.2.3 Phase III

######  Lemma 34 (Phase III).

With probability at least $1-\hat{\eta}j_{*}$ for
$j_{*}=\widetilde{\mathcal{O}}(d)$, initial distribution $\bar{\mu}_{i_{*}}$
(given by from Lemma 33), and target distribution
$\mu_{i_{*}+j_{*}}=\mathcal{N}(0,L^{2}dI_{d})|_{\bar{\mathcal{K}}}$, _Phase
II_ outputs a sample with law $\bar{\mu}_{i_{*}+j_{*}}$ satisfying

$\operatorname{ess\,sup}_{\mu_{i_{*}+j_{*}}}\Bigl{|}\frac{\mathrm{d}\bar{\mu}_{i_{*}+j_{*}}}{\mathrm{d}\mu_{i_{*}+j_{*}}}-1\Bigr{|}\leq
1\,,$

Its oracle complexity is

$\widetilde{\mathcal{O}}\bigl{(}L^{2}d^{3}\log^{6}\frac{1}{\hat{\eta}}\bigr{)}\,,$

using the $\mathsf{Proximal\ sampler}$. The step size scheme for each inner
phase is presented in the proof.

###### Proof.

In this phase, we will first perform the analysis over each doubling and then
aggregate over the doublings.

# Inner phases: We first partition $[1,L^{2}d]$ by a sequence of doubling
parts, where the terminal $\sigma^{2}$ is at least double of the initial
$\sigma^{2}$ in each doubling part. Clearly, the number of doubling parts is
$\log_{2}(L^{2}d)$.

Let $\sigma^{2}$ be an initial variance of a given doubling part. Since we
have

$\Bigl{(}1+\frac{\sigma^{2}}{L^{2}d}\Bigr{)}^{\frac{L^{2}d}{\sigma^{2}}}\geq
2\quad\text{for all }d\geq 1\text{ and }\sigma\leq L\sqrt{d}\,,$

the number of inner phases within the doubling part is at most
$L^{2}d/\sigma^{2}\leq L^{2}d$. Therefore, the total number of inner phases
during Phase III is $\widetilde{\mathcal{O}}(L^{2}d)$.

Warmness: Let $j\geq i_{*}$. By the construction of the algorithm, for all
$x\in\bar{\mathcal{K}}$,

$\frac{\mathrm{d}\mu_{j}(x)}{\mathrm{d}\mu_{j+1}(x)}=\frac{\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}}{\exp\bigl{(}-\frac{1}{2\sigma_{j+1}^{2}}\|x\|^{2}\bigr{)}}\,\frac{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j+1}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\leq\frac{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j+1}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\,.$

As $\|x\|^{2}\leq L^{2}d$ on $\bar{\mathcal{K}}$, we have

$\exp\bigl{(}-\frac{1}{2\sigma_{j+1}^{2}}\|x\|^{2}\bigr{)}=\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\cdot\exp\bigl{(}\frac{1}{2(L^{2}d+\sigma_{j}^{2})}\|x\|^{2}\bigr{)}\leq\sqrt{e}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,.$

As a result,

$\frac{\mathrm{d}\mu_{j}(x)}{\mathrm{d}\mu_{j+1}(x)}\leq\frac{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j+1}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\leq\frac{\sqrt{e}\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2\sigma_{j}^{2}}\|x\|^{2}\bigr{)}\,\mathrm{d}x}=\sqrt{e}\,.$

Furthermore, if at each previous phase, we have

$\displaystyle\operatorname{ess\,sup}_{\mu_{j}}\Bigl{|}\frac{\mathrm{d}\bar{\mu_{j}}}{\mathrm{d}\mu_{j}}-1\Bigr{|}\leq
1\,,$ (5.4)

then this implies warmness of constant order

$\operatorname{ess\,sup}_{\mu_{j+1}}\Bigl{|}\frac{\mathrm{d}\bar{\mu_{j}}}{\mathrm{d}\mu_{j+1}}\Bigr{|}\leq
2\sqrt{e}\,.$

Final complexity: We first bound the total query complexity of one doubling
part. Let $\sigma^{2}$ be an initial variance of a given doubling part. For
any intermediate variance $\sigma_{j}^{2}$ within the part, the query
complexity of the $\mathsf{Proximal\ sampler}$ with
$h^{-1}=\widetilde{\Theta}(d^{2}\log^{2}\frac{L}{\eta})$, iteration number
$n=\widetilde{\mathcal{O}}(d^{2}\sigma_{j}^{2}\log^{2}\frac{L}{\hat{\eta}})$,
and warmness $M\geq 2\sqrt{e}$ can be obtained from (5.2) as

$\widetilde{\mathcal{O}}\bigl{(}d^{2}\sigma_{j}^{2}\log^{6}\frac{L}{\hat{\eta}}\bigr{)}\,.$

This achieves $\log 2$-accuracy in $\mathcal{R}_{\infty}$ in each update,
which also implicitly verifies (5.4). Since within the doubling part,
$\sigma_{j}^{2}\leq 4\sigma^{2}$ (where $\sigma^{2}$ is the initial variance
within the doubling part) and the number of inner phases is at most
$L^{2}d/\sigma^{2}$, the total query complexity during Phase III, aggregating
over all the doubling parts, is

$\log_{2}(L^{2}d)\cdot\frac{L^{2}d}{\sigma^{2}}\cdot\widetilde{\mathcal{O}}\Bigl{(}d^{2}(4\sigma^{2})\log^{6}\frac{L}{\hat{\eta}}\Bigr{)}=\widetilde{\mathcal{O}}\bigl{(}L^{2}d^{3}\log^{6}\frac{1}{\hat{\eta}}\bigr{)}\,,$

achieving $\log 2$-accuracy in $\mathcal{R}_{\infty}$ for the law of the final
sample. ∎

####  5.2.4 Phase IV

######  Lemma 35 (Phase IV).

With probability at least $1-\hat{\eta}$, initial distribution
$\bar{\mu}_{i_{*}+j_{*}}$ (given by from Lemma 34), and target distribution
$\bar{\pi}\propto\mathds{1}_{\bar{\mathcal{K}}}$, _Phase IV_ outputs a sample
with law $\nu$ satisfying

$\operatorname{ess\,sup}_{\bar{\pi}}\Bigl{|}\frac{\mathrm{d}\nu}{\mathrm{d}\bar{\pi}}-1\Bigr{|}\leq\frac{\varepsilon}{3}\,,$

Its oracle complexity is

$\widetilde{\mathcal{O}}\bigl{(}L^{2}d^{3}\log^{6}\frac{1}{\hat{\eta}\varepsilon}\bigr{)}\,,$

using the $\mathsf{Proximal\ sampler}$ with
$h^{-1}=\widetilde{\Theta}(d^{2}\log\frac{\log(L/\varepsilon)}{\hat{\eta}})$
and
$\widetilde{\mathcal{O}}(L^{2}d^{3}\log^{2}\frac{1}{\hat{\eta}\varepsilon})$
iterations.

###### Proof.

# Inner phases: There is only a single update in this phase.

Warmness: Let $\bar{\pi}$ be the uniform distribution over
$\bar{\mathcal{K}}$. First, bounding $\exp(-\frac{1}{2L^{2}d}\|x\|^{2})\leq 1$
and, since $\bar{\mathcal{K}}\subseteq B_{L\sqrt{d}}(0)$,

$\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2L^{2}d}\|x\|^{2}\bigr{)}\,\mathrm{d}x\geq\operatorname{vol}(\bar{\mathcal{K}})\exp(-\frac{L^{2}d}{2L^{2}d})=\operatorname{vol}(\bar{\mathcal{K}})/\sqrt{e}\,.$

Then, we can deduce

$\frac{\mathrm{d}\mu_{i_{*}+j_{*}}(x)}{\mathrm{d}\bar{\pi}(x)}=\frac{\mathrm{d}\mu_{i_{*}+j_{*}}(x)}{1/\operatorname{vol}(\bar{\mathcal{K}})}=\operatorname{vol}(\bar{\mathcal{K}})\cdot\frac{\exp\bigl{(}-\frac{1}{2L^{2}d}\|x\|^{2}\bigr{)}}{\int_{\bar{\mathcal{K}}}\exp\bigl{(}-\frac{1}{2L^{2}d}\|x\|^{2}\bigr{)}\,\mathrm{d}x}\leq\sqrt{e}\,.$

As before, since we have

$\operatorname{ess\,sup}_{\mu_{i_{*}+j_{*}}}\Bigl{|}\frac{\mathrm{d}\bar{\mu}_{i_{*}+j_{*}}}{\mathrm{d}\mu_{i_{*}+j_{*}}}-1\Bigr{|}\leq
1\,,$

then this implies warmness of constant order

$\operatorname{ess\,sup}_{\mu_{i_{*}+j_{*}}}\mathopen{}\mathclose{{}\left\lvert\frac{\mathrm{d}\bar{\mu}_{i_{*}+j_{*}}}{\mathrm{d}\bar{\pi}}}\right\rvert\leq
2\sqrt{e}\,.$

Final complexity: Here, using Theorem 18 (instead of the Gaussian results)
with $M=2\sqrt{e}$ and $D=2L\sqrt{d}$, we obtain the query complexity of

$\widetilde{\mathcal{O}}\bigl{(}L^{2}d^{3}\log^{6}\frac{1}{\hat{\eta}\varepsilon}\bigr{)}\,.\qed$

We now put together these four lemmas to prove Theorem 28.

###### Proof of Theorem 28.

Given the start $\mu_{0}$, we first apply Lemma 32 to produce $\bar{\mu}_{1}$.
Then apply Lemma 33 to produce $\bar{\mu}_{i_{*}}$, Lemma 34 to produce
$\bar{\mu}_{i_{*}+j_{*}}$, and finally Lemma 35 to produce $\nu$. Their
complexity is dominated by that of Lemma 35, which is

$\widetilde{\mathcal{O}}\bigl{(}L^{2}d^{3}\log^{6}\frac{1}{\hat{\eta}\varepsilon}\bigr{)}=\widetilde{\mathcal{O}}\bigl{(}C^{2}d^{3}\log^{8}\frac{1}{\hat{\eta}\varepsilon}\bigr{)}\,,$

with the guarantee
$\mathcal{R}_{\infty}(\nu\mathbin{\|}\bar{\pi})\leq\frac{\varepsilon}{3}$.

From our choice of $L$, it follows that
$\mathcal{R}_{\infty}(\bar{\pi}\mathbin{\|}\pi)\leq\log(1+\frac{2\varepsilon}{3})$
due to

$\sup_{x\in\mathcal{K}}\frac{\bar{\pi}(x)}{\pi(x)}=\frac{\operatorname{vol}(\mathcal{K})}{\operatorname{vol}(\bar{\mathcal{K}})}\leq
1+\frac{2\varepsilon}{3}\,.$

Therefore,

$\displaystyle\mathcal{R}_{\infty}(\nu\mathbin{\|}\pi)$
$\displaystyle=\operatorname{ess\,sup}_{\pi}\log\frac{\mathrm{d}\nu}{\mathrm{d}\pi}\leq\mathcal{R}_{\infty}(\nu\mathbin{\|}\bar{\pi})+\mathcal{R}_{\infty}(\bar{\pi}\mathbin{\|}\pi)$
$\displaystyle\leq\frac{\varepsilon}{3}+\log\bigl{(}1+\frac{2\varepsilon}{3}\bigr{)}\leq\varepsilon\,.$

Secondly, we need to quantity the probability of successfully generating a
sample. This is given, by a union bound, as $1-\hat{\eta}(2+i_{*}+j_{*})$.
With $\hat{\eta}$ chosen according to (5.1), $\mathsf{Proximal\ Gaussian\
cooling}$ has at least $1-\eta$ probability of success. Putting all these
together concludes the proof. ∎

##  6 Conclusion

We have presented $\mathsf{Proximal\ Gaussian\ cooling}$, an algorithm which
achieves $\varepsilon$-closeness in $\mathcal{R}_{\infty}$ with
$\mathcal{O}(d^{3}\operatorname{polylog}\frac{1}{\varepsilon})$ query
complexity. Here we note some possible future directions for research.

  * •

In this work, we showed that it is possible to extend the guarantees for
$\mathsf{In\text{-}and\text{-}Out}$ for Gaussians of the form
$\mathcal{N}(0,\sigma^{2}I_{d})$. Naturally, one asks whether it is possible
to sample from general distributions of the type $e^{-f}|_{\mathcal{K}}$, when
$f$ has some nice properties? The applications of this problem are
innumerable, and the techniques required are far from trivial. It is likely
that one would need to devise a different annealing scheme in order to control
the initial warmness parameter, or find a new strategy which bypasses the need
for annealing entirely.

  * •

This work assumed the membership oracle model for $\mathcal{K}$, but it is
interesting to consider alternative oracle models. For instance, one can use
more features of the geometry, by assuming access to a self-concordant barrier
or a Riemannian metric on $\mathcal{K}$. Alternatively, one can assume access
to a stronger oracle such as a separation oracle, which returns not just a
binary response to if $x\in\mathcal{K}$, but when $x\not\in\mathcal{K}$ also
gives a hyperplane separating $x$ and $\mathcal{K}$. It is as yet unknown if
such a model can improve the query complexity, even for uniform sampling. It
is also possible that a different annealing strategy would be beneficial in
this context.

## Acknowledgements

We thank Sinho Chewi for helpful feedback. YK was supported in part by NSF
awards CCF-2007443 and CCF-2134105. MSZ was supported by NSERC through the
CGS-D program.

## References

  * [AC21] Kwangjun Ahn and Sinho Chewi.  Efficient constrained sampling via the mirror-Langevin algorithm.  Advances in Neural Information Processing Systems (NeurIPS), 34:28405–28418, 2021.
  * [AC23] Jason M Altschuler and Sinho Chewi.  Faster high-accuracy log-concave sampling via algorithmic warm starts.  In Symposium on Foundations of Computer Science (FOCS), pages 2169–2176. IEEE, 2023.
  * [AK91] David Applegate and Ravi Kannan.  Sampling and integration of near log-concave functions.  In Symposium on Theory of Computing (STOC), pages 156–163, 1991.
  * [BCJ+19] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul A. Szerlip, Paul Horsfall, and Noah D. Goodman.  Pyro: Deep universal probabilistic programming.  The Journal of Machine Learning Research (JMLR), 20:28:1–28:6, 2019.
  * [BDMP17] 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), volume 65, pages 319–342. PMLR, 2017.
  * [BEL18] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec.  Sampling from a log-concave distribution with projected Langevin Monte Carlo.  Discrete & Computational Geometry (DCG), 59:757–783, 2018.
  * [BGL14] Dominique Bakry, Ivan Gentil, and Michel Ledoux.  Analysis and geometry of Markov diffusion operators, volume 103.  Springer, 2014.
  * [BST14] Raef Bassily, Adam Smith, and Abhradeep Thakurta.  Private empirical risk minimization: Efficient algorithms and tight error bounds.  In Symposium on Foundations of Computer Science (FOCS), pages 464–473. IEEE, 2014.
  * [CB18] Xiang Cheng and Peter Bartlett.  Convergence of Langevin MCMC in KL-divergence.  In Algorithmic Learning Theory (ALT), pages 186–211. PMLR, 2018.
  * [CCSW22] Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono.  Improved analysis for a proximal algorithm for sampling.  In Conference on Learning Theory (COLT), pages 2984–3014. PMLR, 2022.
  * [CDWY18] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu.  Fast MCMC sampling algorithms on polytopes.  The Journal of Machine Learning Research (JMLR), 19(1):2146–2231, 2018.
  * [CEL+22] Sinho Chewi, Murat A Erdogdu, Mufan Li, Ruoqi Shen, and Shunshi Zhang.  Analysis of Langevin Monte Carlo from Poincaré to Log-Sobolev.  In Conference on Learning Theory (COLT), pages 1–2. PMLR, 2022.
  * [CG23] Yuansi Chen and Khashayar Gatmiry.  When does Metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjusted Langevin algorithm?  arXiv preprint arXiv:2304.04724, 2023.
  * [CLA+21] Sinho Chewi, Chen Lu, Kwangjun Ahn, Xiang Cheng, Thibaut Le Gouic, and Philippe Rigollet.  Optimal dimension dependence of the Metropolis-adjusted Langevin algorithm.  In Conference on Learning Theory (COLT), pages 1260–1300. PMLR, 2021.
  * [CLL19] Yu Cao, Jianfeng Lu, and Yulong Lu.  Exponential decay of Rényi divergence under Fokker–Planck equations.  Journal of Statistical Physics, 176:1172–1184, 2019.
  * [CV16] Ben Cousins and Santosh S Vempala.  A practical volume algorithm.  Mathematical Programming Computation, 8(2):133–160, 2016.
  * [CV18] Ben Cousins and Santosh S Vempala.  Gaussian Cooling and $O^{*}(n^{3})$ algorithms for volume and Gaussian volume.  SIAM Journal on Computing (SICOMP), 47(3):1237–1273, 2018.
  * [CZS22] Xiang Cheng, Jingzhao Zhang, and Suvrit Sra.  Theory and algorithms for diffusion processes on Riemannian manifolds.  arXiv preprint arXiv:2204.13665, 2022.
  * [Dal17] Arnak Dalalyan.  Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent.  In Conference on Learning Theory (COLT), pages 678–689. PMLR, 2017.
  * [DCM19] Daniel Delahaye, Supatcha Chaimatanan, and Marcel Mongeau.  Simulated annealing: From basics to applications.  Handbook of metaheuristics, pages 1–35, 2019.
  * [DCWY18] Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu.  Log-concave sampling: Metropolis-Hastings algorithms are fast!  In Conference on Learning Theory (COLT), pages 793–797. PMLR, 2018.
  * [DFK91] Martin Dyer, Alan Frieze, and Ravi Kannan.  A random polynomial-time algorithm for approximating the volume of convex bodies.  Journal of the ACM (JACM), 38(1):1–17, 1991.
  * [DMLM03] Pierre Del Moral, Michel Ledoux, and Laurent Miclo.  On contraction properties of Markov kernels.  Probability Theory and Related Fields, 126(3):395–420, 2003.
  * [DMM19] Alain Durmus, Szymon Majewski, and Blazej Miasojedow.  Analysis of Langevin Monte Carlo via convex optimization.  The Journal of Machine Learning Research (JMLR), 20(1):2666–2711, 2019.
  * [Doo53] Joseph L Doob.  Stochastic Processes.  John Wiley and Sons,, 1953.
  * [DT12] Arnak S Dalalyan and Alexandre B Tsybakov.  Sparse regression learning by aggregation and Langevin Monte-Carlo.  Journal of Computer and System Sciences (JCSS), 78(5):1423–1443, 2012.
  * [EHZ22] Murat A Erdogdu, Rasa Hosseinzadeh, and Shunshi Zhang.  Convergence of Langevin Monte Carlo in chi-squared and Rényi divergence.  In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 8151–8175. PMLR, 2022.
  * [GC11] Mark Girolami and Ben Calderhead.  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214, 2011.
  * [GHZ22] Mert Gürbüzbalaban, Yuanhan Hu, and Lingjiong Zhu.  Penalized Langevin and Hamiltonian Monte Carlo algorithms for constrained sampling.  arXiv preprint arXiv:2212.00570, 2022.
  * [GKV23] Khashayar Gatmiry, Jonathan Kelner, and Santosh S Vempala.  Sampling with barriers: Faster mixing via Lewis weights.  arXiv preprint arXiv:2303.00480, 2023.
  * [GLL22] Sivakanth Gopi, Yin Tat Lee, and Daogao Liu.  Private convex optimization via exponential mechanism.  In Conference on Learning Theory (COLT), pages 1948–1989. PMLR, 2022.
  * [GLL+23] Sivakanth Gopi, Yin Tat Lee, Daogao Liu, Ruoqi Shen, and Kevin Tian.  Private convex optimization in general norms.  In Symposium on Discrete Algorithms (SODA), pages 5068–5089. SIAM, 2023.