# Log-Concave Sampling on Compact Supports: A Versatile Proximal Framework

Lu Yu  CREST/ENSAE Paris<EMAIL_ADDRESS>
###### Abstract

In this paper, we explore sampling from strongly log-concave distributions
defined on convex and compact supports. We propose a general proximal
framework that involves projecting onto the constrained set, which is highly
flexible and supports various projection options. Specifically, we consider
the cases of Euclidean and Gauge projections, with the latter having the
advantage of being performed efficiently using a membership oracle. This
framework can be seamlessly integrated with multiple sampling methods. Our
analysis focuses on Langevin-type sampling algorithms within the context of
constrained sampling. We provide nonasymptotic upper bounds on the ${\sf
W}_{1}$ and ${\sf W}_{2}$ errors, offering a detailed comparison of the
performance of these methods in constrained sampling.

##  1 Introduction

Sampling from probability distributions plays a critical role in various
fields of science and engineering, especially when dealing with convex and
compact sets [ADFDJ03, GCSR95, Stu10]. In this context, the problem involves
sampling from a probability measure $\nu$ on such sets, characterized by its
density function

$\displaystyle\nu(x)=\frac{e^{-U(x)}}{\int_{\mathbb{R}^{p}}e^{-U(s)}{\rm
d}s}\,,$

Here, $U(x)=f(x)+\ell_{\mathcal{K}}(x)$, where $f(x)$ represents a potential
function and $\ell_{\mathcal{K}}(x)$ is an indicator function ensuring $x$
lies within the convex and compact set $\mathcal{K}\subset\mathbb{R}^{p}$.
Specifically, $\ell_{\mathcal{K}}(x)$ takes the form

$\displaystyle\ell_{\mathcal{K}}:=\begin{cases}+\infty&\text{ if
}x\notin\mathcal{K}\\\ 0&\text{ if }x\in\mathcal{K}\,.\end{cases}$

Solving this constrained sampling problem is challenging and has garnered
considerable interest across various fields, including computer science and
statistics. In the realm of computer science, a line of research initiated by
[DFK91] explored polynomial-time algorithms for uniformly sampling convex
bodies. This has been followed by seminal studies on the convergence
properties of the Ball Walk and the Hit-and-Run algorithm toward uniform
density on a convex body or, more broadly, to log-concave densities [KLS97,
KVZ24, Lov99, LS93, LV07, Smi84]. Other Markov Chain Monte Carlo (MCMC)
methods, such as Gibbs sampling [GSL92] and Hamiltonian Monte Carlo [BSU12,
GHZ22, KLSV22], have also been adapted and enhanced to sample from
distributions defined on convex and compact sets.

In recent years, leveraging optimization techniques to facilitate sampling has
become a prevalent approach. By formulating the sampling challenge as an
optimization problem, methods like projected stochastic gradient descent
[BEL15, BEL18, Lam21, Leh23], proximal approaches [BDMP17, DMP18, SR20],
particle-based algorithms [LLK+22], and mirror descent [AC21, CLGL+20, HKRC18,
ZPFP20] have proven effective in navigating the target distribution to
generate samples. Further innovations have emerged from the intersection of
deep learning and neural networks, leading to the development of novel
sampling techniques via generative adversarial networks [GPAM+14] and
variational autoencoders [KW13]. These advanced methodologies offer promising
pathways for sampling from intricate and high-dimensional distributions,
particularly those defined on convex and compact sets [OHHDT22].

In this work, we tackle the challenge posed by a non-smooth target density
$\nu$ by employing a proximal method. This method leverages a regularization
technique that involves projecting onto set $\mathcal{K}$, effectively
transforming the constrained sampling problem into an unconstrained one. Our
framework is versatile, accommodating various projection options, such as the
Euclidean projection—which corresponds to the Moreau envelope of the indicator
function $\ell_{\mathcal{K}}$ [RW09, DMP18, BDMP17, Per16]—and the Gauge
projection [LBGH22, Mha22], which can be efficiently performed using a
membership oracle. Specifically, we introduce a smooth and strongly convex
surrogate distribution that closely mimics the target density. By fine-tuning
the regularization parameter, we can significantly reduce approximation
errors. The advantageous properties of our surrogate distribution enable the
effective application of various sampling techniques. In this work, we
specifically explore the vanilla Langevin Monte Carlo [RT96, Dal17, DM17,
EH21, MHFH+23, RRT17, EMS18, MFWB22, EHZ22], kinetic Langevin Monte Carlo
[CCBJ18, DRD20, SL19, MCC+21, ZCL+23], and the parallelized versions of the
midpoint randomization method for these algorithms [YD24] within the context
of constrained sampling. We derive the convergence rates for these algorithms
in Wasserstein-1 and Wasserstein-2 distances, and provide a detailed
comparison of their performance. To this end, we make the following
contributions.

  * •

We establish both general upper and lower bounds for the distance between the
smooth approximation of the target density and its original form in
Wasserstein-$q$ distance for any $q\geqslant 1$. These bounds are detailed in
Proposition 3 and Proposition 4, respectively.

  * •

In Section 2, we demonstrate that our proposed framework can seamlessly
incorporate various projection options, including both Euclidean and Gauge
projections.

  * •

In Section 3, we incorporate several MCMC sampling methods, such as (kinetic)
Langevin Monte Carlo and the parallelized midpoint randomization method for
these algorithms, into our general framework. We present a detailed
convergence analysis of these methods in both Wasserstein-1 and Wasserstein-2
distances.

In summary, we develop a comprehensive framework specifically designed for
sampling from convex and compact sets, utilizing a regularization technique
that involves projecting onto the constrained sets. This framework is
exceptionally flexible and adaptable, compatible with various projection
operators and sampling methods. Additionally, it facilitates a lucid
comparison of the behaviors of sampling methods when deployed in constrained
sampling scenarios. Notably, our study presents an improved error bound for
Langevin Monte Carlo applied to constrained sampling settings, as well as the
first convergence analysis of kinetic Langevin Monte Carlo and the
parallelized midpoint randomization method for (kinetic) Langevin Monte Carlo
within the context of constrained sampling. This analysis offers new insights
into the dynamics of constrained sampling.

Notation. Denote the $p$-dimensional Euclidean space by $\mathbb{R}^{p}$. The
letter $\boldsymbol{\theta}$ denotes the deterministic vector and its
calligraphic counterpart $\boldsymbol{\vartheta}$ denotes the random vector.
We use $\mathbf{I}_{p}$ and $\mathbf{0}_{p}$ to denote, respectively, the
$p\times p$ identity and zero matrices. Define the relations
$\mathbf{A}\preccurlyeq\mathbf{B}$ and $\mathbf{B}\succcurlyeq\mathbf{A}$ for
two symmetric $p\times p$ matrices $\mathbf{A}$ and $\mathbf{B}$ to mean that
$\mathbf{B}-\mathbf{A}$ is positive semi-definite. The gradient and the
Hessian of a function $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ are denoted by
$\nabla f$ and $\nabla^{2}f$, respectively. Given any pair of measures $\mu$
and $\nu$ defined on $(\mathbb{R}^{p},\mathcal{B}(\mathbb{R}^{p}))$, the
Wasserstein-$q$ distance between $\mu$ and $\nu$ is defined as

$\displaystyle{\sf
W}_{q}(\mu,\nu)=\Big{(}\inf_{\varrho\in\Gamma(\mu,\nu)}\int_{\mathbb{R}^{p}\times\mathbb{R}^{p}}\\!|\\!|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}|\\!|_{2}^{q}\,{\rm
d}\varrho(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})\Big{)}^{1/q},~{}~{}~{}q\geqslant
1\,,$

where the infimum is taken over all joint distributions $\varrho$ that have
$\mu$ and $\nu$ as marginals. For a set $\mathcal{K}\subset\mathbb{R}^{p},$ we
use $\mathcal{K}^{c}$ to denote its complement. The ceiling function maps
$x\in\mathbb{R}$ to the smallest integer greater than or equal to $x$, denoted
by $\lceil{x}\rceil.$

##  2 Smooth Approximation for the Target Density

The absence of smoothness in the target distribution $\nu$ presents
significant challenges because sampling algorithms often depend on the
smoothness properties of the target distribution to effectively explore the
space and produce representative samples. Motivated by this, we consider the
approximation for $\ell_{\mathcal{K}}$ of the form

$\displaystyle\ell_{\mathcal{K}}^{\lambda}(\boldsymbol{\theta}):=\frac{1}{2\lambda^{2}}\,\\!|\\!|\boldsymbol{\theta}-\mathtt{P}_{\mathcal{K}}(\boldsymbol{\theta})|\\!|_{2}^{2}\,,$

where $\lambda>0$ is the tuning parameter and
$\mathtt{P}_{\mathcal{K}}:\mathbb{R}^{p}\to\mathcal{K}$ denotes a projection
operator which projects the vector $x\in\mathbb{R}^{p}$ onto the set
$\mathcal{K}$. Define
$U^{\lambda}(\boldsymbol{\theta}):=f(\boldsymbol{\theta})+\ell_{\mathcal{K}}^{\lambda}(\boldsymbol{\theta}),$
and we the define the corresponding surrogate target density $\nu^{\lambda}$

$\displaystyle\nu^{\lambda}(\boldsymbol{\theta})=\frac{e^{-U^{\lambda}(\boldsymbol{\theta})}}{\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(\boldsymbol{\theta}^{\prime})}{\rm
d}\boldsymbol{\theta}^{\prime}}\,.$ (2.1)

Throughout the paper, we define
$\mu_{k}(\nu)=\int\\!|\\!|\boldsymbol{\theta}|\\!|_{2}^{k}\,\nu({\rm
d}\boldsymbol{\theta})$ as the $k$-th moment of the distribution $\nu$
($k\geqslant 1$), and we assume the convex and compact set $\mathcal{K}$
satisfy the following assumption.

######  Assumption 1.

Given a positive constant $r\in(0,\infty)$, we assume the Euclidean ball
centered at the origin with radius $r$, denoted by $\mathcal{B}_{2}(r)$, is
contained in $\mathcal{K}$

$\displaystyle\mathcal{B}_{2}(r)=\\{\boldsymbol{\theta}\in\mathbb{R}^{p}:\\!|\\!|\boldsymbol{\theta}|\\!|_{2}\leqslant
r\\}\subset\mathcal{K}\,.$

This assumption has been commonly made in the work of constrained sampling
[Lam21, BEL18, BDMP17, GHZ22]. Moreover, we assume that the potential function
$f$ satisfies the following assumption.

######  Assumption 2.

The function $f:\mathbb{R}^{p}\to\mathbb{R}$ is continuously differentiable,
and its Hessian matrix $\nabla^{2}f$ satisfies

$\displaystyle
m\mathbf{I}_{p}\preccurlyeq\nabla^{2}f(\boldsymbol{\theta})\preccurlyeq
M\mathbf{I}_{p},\qquad\forall\boldsymbol{\theta}\in\mathbb{R}^{p}.$

Moreover, we assume the function $f$ is lower bounded and
$\mathbf{0}=\mathop{\mathrm{missing}}{argmin}_{\boldsymbol{\theta}\in\mathcal{K}}f(\boldsymbol{\theta}).$

The assumption that $f$ attains the minimum at the origin simplifies the
presentation of our results. All our results will still hold even if this
condition is not met. Given the minimizer of the function $f$ over
$\mathcal{K},$ denoted by $\boldsymbol{\theta}_{*}$, we simply need to shift
all the coordinates and consider the constrained set
$\mathcal{K}+\boldsymbol{\theta}_{*}$. Additionally, we note that the
approximation $U^{\lambda}$ inherits the strong convexity of $f.$ Moreover, we
require $U^{\lambda}$ to be integrable, continuously differentiable, and
smooth, even if $U$ is not.

######  Condition 2.1.

The function $U^{\lambda}$ is $m$-strongly convex, continuously
differentiable, and $M^{\lambda}$-smooth. Moreover, it holds that
$\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(\boldsymbol{\theta})}{\rm
d}\boldsymbol{\theta}<\infty$.

Below, we present two examples of the approximation
$\ell_{\mathcal{K}}^{\lambda}$, each utilizing a different, commonly used
projection operator. We will demonstrate that these examples effectively
approximate the indicator function $\ell_{\mathcal{K}}$ and the resulting
potential functions fulfill Condition 2.1.

######  Example 2.1.

[Moreau envelope] One example of a projection operator is the Euclidean
operator, which results in the Moreau envelope of the indicator function
$\ell_{\mathcal{K}}$

$\displaystyle\ell_{\mathcal{K}}^{E,\lambda}(\boldsymbol{\theta})=\inf_{\boldsymbol{\theta}^{\prime}\in\mathbb{R}^{p}}\big{(}\ell_{\mathcal{K}}(\boldsymbol{\theta}^{\prime})+\frac{1}{2\lambda^{2}}\,\\!|\\!|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}|\\!|_{2}^{2}\big{)}=\frac{1}{2\lambda^{2}}\,\\!|\\!|\boldsymbol{\theta}-\mathtt{P}^{E}_{\mathcal{K}}(\boldsymbol{\theta})|\\!|_{2}^{2}\,,$

where $\mathtt{P}_{\mathcal{K}}^{E}$ denotes the Euclidean projection on
$\mathcal{K}$.

Define the corresponding surrogate potential
$U^{E,\lambda}(\boldsymbol{\theta}):=f(\boldsymbol{\theta})+\ell_{\mathcal{K}}^{E,\lambda}(\boldsymbol{\theta}).$
This surrogate potential $U^{E,\lambda}$ satisfies the following property.

######  Lemma 1.

Assume Assumptions 1 and 2 hold, the potential $U^{E,\lambda}$ satisfies
Condition 2.1 with $M^{\lambda}=1/\lambda^{2}+M$.

######  Example 2.2.

[Gauge projection] Another example of the projection operator is the Gauge
projection $\mathtt{P}_{\mathcal{K}}^{G}:\mathbb{R}^{p}\to\mathbb{R}$, which
is defined as

$\displaystyle\mathtt{P}_{\mathcal{K}}^{G}(\boldsymbol{\theta}):=\frac{\boldsymbol{\theta}}{g_{\mathcal{K}}(\boldsymbol{\theta})}$

with a variation of Gauge function (also known as the Minkowski function)
$g_{\mathcal{K}}:\mathbb{R}^{p}\to\mathbb{R}$ of set $\mathcal{K}$

$\displaystyle g_{\mathcal{K}}(\boldsymbol{\theta}):=\inf\\{t\geqslant
1:\boldsymbol{\theta}\in t\mathcal{K}\\}\,.$ (2.2)

This function encapsulates the scaling required to project
$\boldsymbol{\theta}$ onto the set $\mathcal{K}.$

Set
$\ell_{\mathcal{K}}^{G,\lambda}(\boldsymbol{\theta})=\frac{1}{2\lambda^{2}}\,\\!|\\!|\boldsymbol{\theta}-\mathtt{P}^{G}_{\mathcal{K}}(\boldsymbol{\theta})|\\!|_{2}^{2}\,.$
The corresponding surrogate function is defined via
$U^{G,\lambda}(\boldsymbol{\theta}):=f(\boldsymbol{\theta})+\ell_{\mathcal{K}}^{G,\lambda}(\boldsymbol{\theta}).$
This surrogate potential $U^{G,\lambda}$ satisfies the following property.

######  Lemma 2.

Assume Assumptions 1 and 2 hold, the surrogate potential $U^{G,\lambda}$
satisfies Condition 2.1 with $M^{\lambda}=1/\lambda^{2}+M$.

In the following, we aim to measure the discrepancy between the surrogate
distribution and the target distribution. To this end, we establish an upper
bound in ${\sf W}_{q}$ distance between the target distribution $\nu$ and the
approximated distribution $\nu^{\lambda}$ defined in (2.1). Our aim is to
demonstrate that as $\lambda$ approaches zero, the distance between these two
distributions converges to zero.

######  Proposition 3.

Under Assumption 1, assume that the potential function $f$ is convex with 0 as
its minimizer. Then, for any $q\geqslant 1$ and any
$\lambda\in\big{(}0,\frac{r}{2(p+q)}\big{)},$ it holds that

$\displaystyle{\sf W}_{q}^{q}(\nu,\nu^{\lambda})\leqslant
3\mu_{q}(\nu){\lambda}\,(2p+q)/r\,,$

where $\mu_{q}(\nu)=\int\\!|\\!|\boldsymbol{\theta}|\\!|_{2}^{q}\,\nu({\rm
d}\boldsymbol{\theta}),q\geqslant 1.$

In this proposition, we provide a general and precise quantification of the
distance between the approximate distribution $\nu^{\lambda}$ to the target
distribution $\nu$ in ${\sf W}_{q}$ distance for any $q\geqslant 1.$ When
$q=1,$ our result aligns with the the bound established for ${\sf
W}_{1}(\nu,\nu^{\lambda})$ as stated in [BDMP17, Proposition 5].

We also derive the lower bound for the distance between the approximate
distribution $\nu^{\lambda}$ to the target distribution $\nu$ in ${\sf W}_{q}$
distance, specifically for the scenario where $\mathcal{K}=\mathcal{B}_{2}(r)$
and the projections discussed in Examples 2.1 and 2.2. In this case, the Gauge
projection and Euclidean projection coincide, yielding the following results.

######  Proposition 4.

Let $\mathcal{K}=\mathcal{B}_{2}(r)$ and assume the potential function $f$ is
convex and $M$-smooth, with $\mathbf{0}$ as its minimizer and
$f(\mathbf{0})=0.$ Set $\ell_{\mathcal{K}}=\ell_{\mathcal{K}}^{E,\lambda}.$
Under Assumptions 1, for any $q\geqslant 1$ and any
$\lambda\in\big{(}0,\frac{r}{2(p+q)}\wedge\frac{1}{\sqrt{M}}\big{)},$ it holds
that

$\displaystyle{\sf
W}_{q}^{q}(\nu,\nu^{\lambda})\geqslant\mu_{q}(\nu)\min\bigg{\\{}\Big{|}\Big{(}\sqrt{\frac{\pi}{8}}e^{-\frac{3}{4}Mr^{2}}\frac{\lambda(p+q)}{r+3\lambda
p}+\frac{r}{r+3\lambda
p}\Big{)}^{1/q}-1\Big{|}^{q},\Big{|}\Big{(}\frac{3\lambda(p+q)}{r}+1\Big{)}^{1/q}-1\Big{|}^{q}\bigg{\\}}\,.$

Before we discuss the results of the lower bound, we present the following
corollary.

######  Corollary 5.

Under the assumptions stated in Proposition 4, when $\lambda=o(r/p)$ and
$q=1$, it holds that

$\displaystyle{\sf W}_{1}(\nu,\nu^{\lambda})\geqslant
C_{M,r}\mu_{1}(\nu)\frac{p+1}{r}\lambda\,,$

where $C_{M,r}>0$ is a constant that depends on $M$ and $r$ exponentially.

Based on the results outlined in the preceding corollary, when $q=1$ and with
a sufficiently small $\lambda,$ the lower bound aligns with the rate estimates
from Proposition 3 up to a constant factor, thereby confirming the optimality
of the upper bound for the case of $q=1.$ However, when $q>1$ and $\lambda$ is
small, the lower bound follows the order
$\big{(}(1+\lambda)^{1/q}-1\big{)}^{q}$, which increases at a rate slower than
$\lambda.$ This indicates the suboptimality111We conjecture that this
suboptimality arises from the general technique used in our proof. With
specific assumptions about the density, a refined lower bound for the
Wasserstein distance between $\nu$ and $\nu^{\lambda}$. of the rate stated in
Proposition 3, highlighting a difference in behavior for higher values of $q.$

##  3 Proximal Langevin Algorithms

A central theme of this work is the approximation of the non-smooth potential
$\nu$ using a carefully crafted smooth surrogate density $\nu^{\lambda}$.
Owing to the desirable properties of the surrogate density $\nu^{\lambda}$, we
can utilize various sampling methods that are effective for distributions with
smooth and strongly convex densities.

In this section, we explore Langevin Monte Carlo (LMC) and its variant,
Kinetic Langevin Monte Carlo (KLMC), along with the parallelized randomized
midpoint method for both algorithms, denoted as pRLMC and pRKLMC,
respectively, to sample from the density $\nu^{\lambda}$. We combine the
discretization error with the approximation error analyzed in Section 2 to
evaluate errors in Wasserstein-1 and Wasserstein-2 distances between the
sample distribution and the target density $\nu$ in the context of constrained
sampling. Table 1 below offers a comparison of the results for the four
sampling methods mentioned above, using the specific projection operators
discussed in Example 2.1 and Example 2.2. For simplicity, we omit the
constants and logarithmic terms that appear in the bounds. The table
illustrates that when the number of parallel steps $R>1$, the randomized
midpoint method markedly enhances the performance of both the vanilla Langevin
Monte Carlo and kinetic Langevin Monte Carlo algorithms. When the number of
parallel steps $R=1$, the pRLMC and pRKLMC algorithms correspond exactly to
the randomized midpoint method applied to Langevin Monte Carlo (RLMC) and
kinetic Langevin Monte Carlo (RKLMC), respectively, as detailed in [YD24]. The
convergence rate for RLMC is comparable to that of LMC when assessing the
error in ${\sf W}_{2}$ distance. Further details are provided following
Corollary 9.

| LMC | KLMC | pRLMC | pRKLMC
---|---|---|---|---
${\sf W}_{1}$ | $\widetilde{\mathcal{O}}(\varepsilon^{-4})$ | $\widetilde{\mathcal{O}}(\varepsilon^{-4})$ | $\widetilde{\mathcal{O}}(R^{-1/3}\varepsilon^{-10/3})$ | $\widetilde{\mathcal{O}}(R^{-1/3}\varepsilon^{-8/3})$
${\sf W}_{2}$ | $\widetilde{\mathcal{O}}(\varepsilon^{-6})$ | $\widetilde{\mathcal{O}}(\varepsilon^{-7})$ | $\widetilde{\mathcal{O}}(R^{-1/3}\varepsilon^{-6})$ | $\widetilde{\mathcal{O}}(R^{-1/3}\varepsilon^{-5})$
Table 1: The number of iterations required by
$\\{\text{L,KL,pRL,pRKL}\\}\text{MC}$ algorithms to achieve an error bounded
by $\varepsilon\sqrt{p/m}$ in ${\sf W}_{1}$ distance and ${\sf W}_{2}$
distance. $R\geqslant 1$ denotes the number of parallel steps.

###  3.1 Langevin Monte Carlo

Let $\boldsymbol{\vartheta}_{0}$ be a random vector drawn from a distribution
$\nu^{\lambda}$ on $\mathbb{R}^{p}$ and let
$\boldsymbol{W}=(\boldsymbol{W}_{t}:t\geqslant 0)$ be a $p$-dimensional
Brownian motion that is independent of $\boldsymbol{\vartheta}_{0}$. To sample
from the approximation distribution $\nu^{\lambda}$, we consider the vanilla
Langevin diffusion, which is a strong solution to the stochastic differential
equation

$\displaystyle{\rm d}\boldsymbol{L}_{t}^{\sf LD}=-\nabla
U^{\lambda}(\boldsymbol{L}_{t}^{\sf LD})\,{\rm d}t+\sqrt{2}\,{\rm
d}\boldsymbol{W}_{t},\qquad t\geqslant 0,\qquad\boldsymbol{L}^{\sf
LD}_{0}=\boldsymbol{\vartheta}_{0}.$ (3.1)

This equation has a unique strong solution, which is a continuous-time Markov
process, termed Langevin diffusion. Under the further assumptions on the
potential $U^{\lambda}$, such as strong convexity, the Langevin diffusion is
ergodic, geometrically mixing and has $\nu^{\lambda}$ as its unique invariant
distribution [Bha78]. Moreover, we can sample from the distribution defined by
$\nu^{\lambda}$ by using a suitable discretization of the Langevin diffusion.
LMC algorithm is based on this idea, combining the considerations with the
Euler discretization. Specifically, for small values of $h\geqslant 0$ and
$\Delta_{h}\boldsymbol{W}_{t}=\boldsymbol{W}_{t+h}-\boldsymbol{W}_{t}$, the
following approximation holds

$\displaystyle\boldsymbol{L}^{\sf LD}_{t+h}$
$\displaystyle=\boldsymbol{L}^{\sf LD}_{t}-\int_{0}^{h}\nabla
U^{\lambda}(\boldsymbol{L}_{t+s}^{\sf LD})\,{\rm
d}s+\sqrt{2}\;\Delta_{h}\boldsymbol{W}_{t}\approx\boldsymbol{L}_{t}^{\sf
LD}-h\nabla U^{\lambda}(\boldsymbol{L}_{t}^{\sf
LD})+\sqrt{2}\;\Delta_{h}\boldsymbol{W}_{t}.$

By repeatedly applying this approximation with a small step-size $h$, we can
construct a Markov chain $(\boldsymbol{\vartheta}^{\sf
LMC}_{k}:k\in\mathbb{N})$ that converges to the target distribution
$\nu^{\lambda}$ as $h$ goes to zero. More precisely,
$\boldsymbol{\vartheta}^{\sf LMC}_{k}\approx\boldsymbol{L}^{\sf LD}_{kh}$, for
$k\in\mathbb{N}$, is given by

$\displaystyle\boldsymbol{\vartheta}_{k+1}^{\sf
LMC}=\boldsymbol{\vartheta}_{k}^{\sf LMC}-h\nabla
U^{\lambda}(\boldsymbol{\vartheta}_{k}^{\sf
LMC})+\sqrt{2}\,(\boldsymbol{W}_{(k+1)h}-\boldsymbol{W}_{kh}).$

We now provide explicit upper bounds for the error of the LMC algorithm in
terms of ${\sf W}_{1}$ and ${\sf W}_{2}$ distances within the context of
constrained sampling.

######  Theorem 6.

Under Assumptions 1 and 2, we further assume that the potential $U^{\lambda}$
satisfies Condition 2.1. Let the step size $h\in(0,1/M^{\lambda})$ and the
tuning parameter $\lambda\in\big{(}0,\frac{r}{2(p+2)}\big{)}$. Then, for every
$n\geqslant 1$, the distribution $\nu_{n}^{\sf LMC}$ of
$\boldsymbol{\vartheta}_{n}^{\sf LMC}$ satisfies

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu)$ $\displaystyle\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu)+2\sqrt{\frac{3\mu_{2}(\nu)\lambda(2p+2)}{r}}+\sqrt{\frac{2phM^{\lambda}}{m}}\,.$

Moreover, when the initial distribution $\nu_{0}^{\sf LMC}$ is set to be the
Dirac measure at the minimizer of $f,$ it holds that

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf LMC},\nu)$ $\displaystyle\leqslant
e^{-\frac{mnh}{2}}\sqrt{\frac{p}{m}}+\frac{3\mu_{1}(\nu)\lambda(2p+1)}{r}+\sqrt{\frac{2phM^{\lambda}}{m}}\,.$

The term $M^{\lambda}$ introduces an additional factor of $\lambda$ into the
convergence rate. To clarify this dependence on $\lambda$, we specify in the
following corollary that $M^{\lambda}=M+1/\lambda^{2}$. This specification
corresponds to the specific projection operators used in Example 2.1 and
Example 2.2. We then optimize $\lambda$ to obtain the results presented below.

######  Corollary 7.

Let $\varepsilon\in(0,1)$ be a small number, and
$M^{\lambda}=M+\frac{1}{\lambda^{2}}.$
(a) Set $\lambda=\sqrt{2/3}(ph/m)^{1/4}\sqrt{r/(\mu_{1}(\nu)(2p+1))},$ choose
$h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=2^{-10}3^{-2}(p/m)r^{2}\big{(}\mu_{1}(\nu)(2p+1)\big{)}^{-2}\varepsilon^{4}\quad\text{and}\quad
n\geqslant\frac{2}{mh}\log(2/\varepsilon)\,,$

then we have ${\sf W}_{1}(\nu_{n}^{\sf
LMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

(b) Set $\lambda=2^{-2/3}3^{-1/3}(rh)^{1/3}(\mu_{2}(\nu)m)^{-1/3}$, choose
$h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=2^{-20}3^{-2}r^{2}(\mu_{2}(\nu)m)^{-2}\varepsilon^{6}\quad\text{and}\quad
n\geqslant\frac{2}{mh}\log(2/\varepsilon)\,,$

then we have ${\sf W}_{2}(\nu_{n}^{\sf
LMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

According to this corollary, when evaluating the error in ${\sf W}_{1}$
distance, the required sample size $n$ to achieve a prespecified error level
is of order $\tilde{\mathcal{O}}(\varepsilon^{-4}).$ This represents an
improvement over the rate of $\tilde{\mathcal{O}}(\varepsilon^{-6})$ obtained
in [BDMP17]. Additionally, our general framework provides the convergence rate
in ${\sf W}_{2}$ for the proximal LMC algorithm, thus complementing the study
of the proximal LMC presented in [BDMP17].

###  3.2 Parallelized Randomized Midpoint Discretization of Vanilla Langevin
Diffusion

As an alternative to the Euler discretization method discussed in Section 3.1
for the stochastic differential equation (3.1), we explore another
discretization framework known as the randomized midpoint method, introduced
in [SL19]. This method enables the discretization and simulation of the
Langevin diffusion (3.1), ultimately converging to the distribution
$\nu^{\lambda}$. Building upon the foundations laid by [SL19, YKD23], [YD24]
developed a parallel computing scheme for this algorithm, significantly
enhancing the efficiency of the sampling process. The formal definition of the
parallelized randomized midpoint method for the Langevin Monte Carlo (pRLMC)
is defined in Algorithm 1. For simplicity, the superscript $\sf pRLMC$ is
omitted therein.

Algorithm 1 Parallelized RLMC (pRLMC)

Input: number of parallel steps $R$, number of sequential iterations $Q$, step
size $h$, number of iterations $K$, initial point
$\boldsymbol{\vartheta}_{0}.$
Output: iterate $\boldsymbol{\vartheta}_{K+1}$

1:for $k=1$ to $K$ do

2: Draw $U_{kr}$ uniformly from $\big{[}\frac{r-1}{R},\frac{r}{R}\big{]}$,
$r=1,\dots,R$.

3: Generate
$\boldsymbol{\xi}_{kr}=\boldsymbol{W}_{(k+U_{kr})h}-\boldsymbol{W}_{kh},r=1,\dots,R.$

4: Set
$\boldsymbol{\vartheta}_{k}^{(0,r)}=\boldsymbol{\vartheta}_{k},r=1,\dots,R$.

5: for $q=1$ to $Q-1$ do

6: for $r=1$ to $R$ in parallel do

7: $a_{kj}=\min\\{\frac{1}{R},U_{kr}-\frac{j-1}{R}\\},j=1,\dots,r$

8:
$\boldsymbol{\vartheta}_{k}^{(q,r)}=\boldsymbol{\vartheta}_{k}-h\sum_{j=1}^{r}a_{kj}\nabla
U^{\lambda}\big{(}\boldsymbol{\vartheta}_{k}^{(q-1,j)}\big{)}+\sqrt{2}\boldsymbol{\xi}_{kr}$.

9: end for

10: end for

11:
$\boldsymbol{\vartheta}_{k+1}=\boldsymbol{\vartheta}_{k}-\frac{h}{R}\sum_{r=1}^{R}\nabla
U^{\lambda}\big{(}\boldsymbol{\vartheta}_{k}^{(Q-1,r)}\big{)}+\sqrt{2}(\boldsymbol{W}_{(k+1)h}-\boldsymbol{W}_{kh})$.

12:end for

In parallel

We now introduce the primary findings in the work, which establish explicit
upper bounds for the error associated with the pRLMC algorithm, measured in
${\sf W}_{1}$ and ${\sf W}_{2}$ distances, specifically within the framework
of constrained sampling.

######  Theorem 8.

Under Assumptions 1 and 2, we further assume that the potential $U^{\lambda}$
satisfies Condition 2.1. Let the the tuning parameter
$\lambda\in\big{(}0,\frac{r}{2(p+2)}\big{)}$. Choose the step size $h$ such
that $M^{\lambda}h\leqslant 1/10.$ Then, for every $n\geqslant 1$, the
distribution $\nu_{n}^{\sf pRLMC}$ of $\boldsymbol{\vartheta}_{n}^{\sf pRLMC}$
satisfies

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf pRLMC},\nu)$
$\displaystyle\leqslant\Big{(}1+\sqrt{\kappa
M^{\lambda}h}\big{(}0.82(M^{\lambda}h)^{Q-1}+0.94M^{\lambda}h/R\big{)}\Big{)}e^{-mnh/2}{\sf
W}_{2}(\nu_{0}^{\sf pRLMC},\nu)$ $\displaystyle\qquad+\sqrt{\kappa
M^{\lambda}h}\big{(}3.98(M^{\lambda}h)^{Q-1}+6.91M^{\lambda}h/\sqrt{R}\big{)}\sqrt{p/m}$
$\displaystyle\qquad\Big{(}2+\sqrt{\kappa
M^{\lambda}h}\big{(}0.82(M^{\lambda}h)^{Q-1}+0.94M^{\lambda}h/R\big{)}\Big{)}\sqrt{3\mu_{2}(\nu)\lambda(2p+2)/r}\,,$

where $\kappa=M^{\lambda}/m$. Moreover, when the initial distribution
$\nu_{0}^{\sf pRLMC}$ is set to be the Dirac measure at the minimizer of $f,$
it holds that

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf pRLMC},\nu)$
$\displaystyle\leqslant\Big{(}1+\sqrt{\kappa
M^{\lambda}h}\big{(}0.82(M^{\lambda}h)^{Q-1}+0.94M^{\lambda}h/R\big{)}\Big{)}e^{-mnh/2}\sqrt{p/m}$
$\displaystyle\qquad+\sqrt{\kappa
M^{\lambda}h}\big{(}3.98(M^{\lambda}h)^{Q-1}+6.91M^{\lambda}h/\sqrt{R}\big{)}\sqrt{p/m}$
$\displaystyle\qquad+3\mu_{1}(\nu)\lambda(2p+1)/r.$

To our knowledge, this represents the first convergence rate for the pRLMC
algorithm within the context of constrained sampling. When $R=1$, these
results recover the convergence rate for the randomized midpoint method
applied to Langevin Monte Carlo algorithm (RLMC) [HBE20, YKD23, YD24] in
constrained sampling. Below, we define $M^{\lambda}=M+1/\lambda^{2}$, a
setting that corresponds to the projection choices outlined in Examples 2.1
and 2.2. We then optimize $\lambda$ to derive the following corollary.

######  Corollary 9.

Given the number of parallel steps $R,$ we set $Q=\lceil{\log R}\rceil+1.$ Let
$\varepsilon\in(0,1)$ be a small number, and
$M^{\lambda}=M+\frac{1}{\lambda^{2}}.$
(a) Set $\lambda=(\frac{2^{4}\cdot
7}{3})^{1/5}(\frac{p}{R})^{1/10}(\frac{r}{m\mu_{1}(\nu)(2p+1)})^{1/5}h^{3/10},$
choose $h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=15.6^{-10/3}p^{2/5}m^{-1}\Big{(}\frac{r}{\mu_{1}(\nu)(2p+1)}\Big{)}^{8/3}R^{1/3}\varepsilon^{10/3}\quad\text{and}\quad
n\geqslant\frac{2}{mh}\log\Big{(}\frac{2.22}{\varepsilon}\Big{)}\,,$

then we have ${\sf W}_{1}(\nu_{n}^{\sf
pRLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

(b)
$\lambda=320^{2/9}3^{-1/3}\big{(}\frac{rp}{m^{2}\mu_{2}(\nu)(2p+2)R}\big{)}^{1/9}h^{1/3}$,
choose $h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=13^{-6}m^{-7/3}\Big{(}\frac{rp}{\mu_{2}(\nu)(2p+2)}\Big{)}^{8/3}R^{1/3}\varepsilon^{6}\quad\text{and}\quad
n\geqslant\frac{2}{mh}\log\Big{(}\frac{2.22}{\varepsilon}\Big{)}\,,$

then we have ${\sf W}_{2}(\nu_{n}^{\sf
pRLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

Comparing this convergence rate with that obtained for proximal LMC in
Corollary 7, we observe superior performance of pRLMC over LMC when $R>1$.
When assessing the error in ${\sf W}_{2}$ distance for the case of $R=1$, the
sample complexity $\tilde{\mathcal{O}}(\epsilon^{-6})$ aligns with that of
LMC. This outcome results from our specific choice of $\lambda$ during the
optimization of the upper bound concerning $\lambda$, placing it within a
region where the performance of LMC is comparable to that of RLMC.

###  3.3 Kinetic Langevin Monte Carlo

In this part, we explore the application of kinetic Langevin diffusion in
constrained sampling. Recall that the kinetic Langevin process
$\boldsymbol{L}^{\textup{\sf KLD}}$ is a solution to a second-order stochastic
differential equation that can be informally written as

$\displaystyle{\textstyle\frac{1}{\gamma}}\ddot{\boldsymbol{L}}_{t}^{\sf
KLD}+\dot{\boldsymbol{L}}_{t}^{\sf KLD}=-\nabla
U^{\lambda}(\boldsymbol{L}_{t}^{\sf KLD})+\sqrt{2}\,\dot{\boldsymbol{W}}_{t},$
(3.2)

with initial conditions $\boldsymbol{L}_{0}^{\sf
KLD}=\boldsymbol{\vartheta}_{0}$ and $\dot{\boldsymbol{L}}_{0}^{\sf
KLD}=\mathbf{v}_{0}$. In (3.2), $\gamma>0$, $\boldsymbol{W}$ is a standard
$p$-dimensional Brownian motion and dots are used to designate derivatives
with respect to time $t\geqslant 0$. This can be formalized using Itô’s
calculus and introducing the velocity field $\mathbf{V}^{\sf KLD}$ so that the
joint process $(\boldsymbol{L}^{\sf KLD},\mathbf{V}^{\sf KLD})$ satisfies

$\displaystyle{\rm d}\boldsymbol{L}^{\sf KLD}_{t}=\mathbf{V}_{t}^{\sf
KLD}\,{\rm d}t;\quad\tfrac{1}{\gamma}{\rm d}\mathbf{V}^{\sf
KLD}_{t}=-\big{(}\mathbf{V}_{t}^{\sf KLD}+\nabla
U^{\lambda}(\boldsymbol{L}_{t}^{\sf KLD})\big{)}\,{\rm d}t+\sqrt{2}\,{\rm
d}\boldsymbol{W}_{t}.$ (3.3)

Similar to the vanilla Langevin diffusion (3.1), the kinetic Langevin
diffusion $(\boldsymbol{L}^{\sf KLD},\mathbf{V}^{\sf KLD})$ is a Markov
process that exhibits ergodic properties when the potential $U^{\lambda}$ is
strongly convex (see [EGZ19] and references therein). The invariant density of
this process is given by

$\displaystyle
p_{*}(\boldsymbol{\theta}_{,}\mathbf{v})\propto\exp\\{-U^{\lambda}(\boldsymbol{\theta})-{\textstyle\frac{1}{2\gamma}}\|\mathbf{v}\|^{2}\\},\qquad\text{for
all}\quad\boldsymbol{\theta},\mathbf{v}\in\mathbb{R}^{p}.$

Note that the marginal of $p_{*}$ corresponds to $\boldsymbol{\theta}$
coincides with the target density $\nu^{\lambda}$. The kinetic Langevin Monte
Carlo (KLMC) algorithm is a discretized version of KLD equation (3.3), where
the term $\nabla U^{\lambda}(\boldsymbol{L}_{t})$ is replaced by $\nabla
U^{\lambda}(\boldsymbol{L}_{kh})$ on each interval $[kh,(k+1)h)$. The
resulting error bounds are given in the following theorem.

######  Theorem 10.

Under Assumptions 1 and 2, we further assume that the potential $U^{\lambda}$
satisfies Condition 2.1. Let the tuning parameter
$\lambda\in\big{(}0,\frac{r}{2(p+2)}\big{)}$. Choose $\gamma$ and $h$ so that
$\gamma\geqslant 5M^{\lambda}$ and $\sqrt{\kappa}\,\gamma h\leqslant 0.1$,
where $\kappa=M^{\lambda}/m$. Assume that $\boldsymbol{\vartheta}_{0}^{\sf
KLMC}$ is independent of $\mathbf{v}_{0}^{\sf KLMC}$ and that
$\mathbf{v}_{0}^{\sf KLMC}\sim\mathcal{N}_{p}(0,\gamma\mathbf{I}_{p})$. Then,
for any $n\geqslant 1$, the distribution $\nu_{n}^{\textup{ \sf KLMC}}$ of
$\boldsymbol{\vartheta}_{n}^{\textup{\sf KLMC}}$ satisfies

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf KLMC},\nu)$ $\displaystyle\leqslant
2\varrho^{n}{\sf W}_{2}(\nu_{0}^{\sf
KLMC},\nu)+3\sqrt{\frac{3\mu_{2}(\nu)\lambda(2p+2)}{r}}+0.05\sqrt{\frac{\varrho^{n}\mathbb{E}[U^{\lambda}(\boldsymbol{\vartheta}_{0}^{\sf
KLMC})-f(\mathbf{0})]}{m}}+0.9\gamma h\sqrt{\frac{\kappa p}{m}}\,,$

where $\varrho=e^{-mh}$. Moreover, let the staring points
$\boldsymbol{\vartheta}_{0}^{\sf KLMC}=\mathbf{0}$, it holds that

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf KLMC},\nu)$ $\displaystyle\leqslant
2\varrho^{n}\sqrt{\frac{p}{m}}+\frac{3\mu_{1}(\nu)\lambda(2p+1)}{r}+0.9\gamma
h\sqrt{\frac{\kappa p}{m}}\,.$

To our knowledge, these results are the first reported convergence rate for
the KLMC algorithm within the context of constrained sampling. Below, we set
$M^{\lambda}=M+1/\lambda^{2}$, aligning with the projection methods detailed
in Example 2.1 and Example 2.2. We proceed by optimizing $\lambda$ to
establish the subsequent corollary.

######  Corollary 11.

Let $\varepsilon\in(0,1)$ be a small number. Set
$\gamma=5M^{\lambda},M^{\lambda}=M+\frac{1}{\lambda^{2}}$.
(a) Set
$\lambda=3^{3/4}2^{1/8}r^{1/4}\big{(}\mu_{1}(\nu)(2p+1)\big{)}^{-1/4}p^{1/8}m^{-1/4}h^{1/4},$
choose $h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=2^{-7/2}r^{3}\big{(}\mu_{1}(\nu)(2p+1)\big{)}^{-3}p^{3/2}m^{-1}\varepsilon^{4}\quad\text{and}\quad
n\geqslant\frac{1}{mh}\log(4/\varepsilon)\,,$

then we have ${\sf W}_{1}(\nu_{n}^{\sf
KLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

(b) Set $\lambda=2^{1/7}3^{3/7}(r/\mu_{2}(\nu))^{1/7}(h/m)^{2/7}$, choose
$h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=32.26^{-7}m^{-5/2}r^{3}\mu_{2}(\nu)^{-3}\varepsilon^{7}\quad\text{and}\quad
n\geqslant\frac{1}{mh}\log(4/\varepsilon)\,,$

then we have ${\sf W}_{2}(\nu_{n}^{\sf
KLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

We note that our error bounds depend in part on the synchronous coupling
between the KLMC and the KLD. However, for the vanilla Langevin algorithm,
[DMM19] have shown that the dependency of the error bound on $\kappa$ can be
improved by employing alternative couplings. We propose that similar
improvements might be achievable for the KLMC algorithm in constrained
sampling scenarios using non-synchronous coupling. For further insights,
interested readers are directed to [YKD23].

###  3.4  Parallelized Randomized Midpoint Discretization of Kinetic Langevin
Diffusion

The randomized midpoint method, introduced and studied in [SL19], aims at
providing a discretization of the kinetic Langevin process (3.2) that reduces
the bias of sampling as compared to more conventional discretizations. The
parallel computing of this algorithm is outlined in [YD24]. This algorithm is
referred to as pRKLMC, and for the convenience of the readers, we restate it
in Algorithm 2. To ease the notation, we omit the superscript $\sf pRKLMC$
therein.

Algorithm 2 Parallelized RKLMC (pRKLMC)

Input: number of parallel steps $R$, number of sequential iterations $Q$, step
size $h$, friction coefficient $\gamma$, number of iterations $K$, initial
points $\mathbf{v}_{0}$ and $\boldsymbol{\vartheta}_{0}.$
Output: iterates $\boldsymbol{\vartheta}_{K+1}$ and $\mathbf{v}_{K+1}$

1:for $k=1$ to $K$ do

2: Draw $U_{k1},\dots,U_{kR}$ uniformly from
$\big{[}0,\frac{1}{R}\big{]},\dots,\big{[}\frac{R-1}{R},1\big{]},$
respectively.

3: Generate
$\bar{\boldsymbol{W}}_{s}=\boldsymbol{W}_{kh+s}-\boldsymbol{W}_{kh}$

4: Generate
$\boldsymbol{\xi}_{kr}=\int_{0}^{U_{kr}h}\big{(}1-e^{-\gamma(U_{kr}h-s)}\big{)}d\bar{\boldsymbol{W}}_{s},r=1,\dots,R$.

5: Set
$\boldsymbol{\vartheta}_{k}^{(0,r)}=\boldsymbol{\vartheta}_{k},r=1,\dots,R$.

6: for $q=1$ to $Q-1$ do

7: for $r=1$ to $R$ in parallel do

8: $a_{kr}=\frac{1-e^{-\gamma hU_{kr}}}{\gamma}.$

9:
$b_{kj}=\int_{\frac{(j-1)h}{R}}^{h\min(\frac{j}{R},U_{kr})}(1-e^{-\gamma(U_{kr}h-s)})ds,j=1,\dots,r.$

10:
$\boldsymbol{\vartheta}_{k}^{(q,r)}=\boldsymbol{\vartheta}_{k}+a_{kr}\mathbf{v}_{k}-\sum_{j=1}^{r}b_{kj}\nabla
U^{\lambda}\big{(}\boldsymbol{\vartheta}_{k}^{(q-1,j)}\big{)}+\sqrt{2}\boldsymbol{\xi}_{kr}.$

11: end for

12: end for

13:
$\boldsymbol{\vartheta}_{k+1}=\boldsymbol{\vartheta}_{k}+\frac{1-e^{-\gamma
h}}{\gamma}\mathbf{v}_{k}-\sum_{r=1}^{R}\frac{h}{R}(1-e^{-\gamma
h(1-U_{kr})})\nabla
U^{\lambda}\big{(}\boldsymbol{\vartheta}_{k}^{(Q-1,r)}\big{)}+\sqrt{2}\int_{0}^{h}\big{(}1-e^{-\gamma(h-s)}\big{)}d\bar{\boldsymbol{W}}_{s}$.

14: $\mathbf{v}_{k+1}=e^{-\gamma
h}\mathbf{v}_{k}-\gamma\sum_{r=1}^{R}\frac{h}{R}e^{-\gamma h(1-U_{kr})}\nabla
U^{\lambda}\big{(}\boldsymbol{\vartheta}_{k}^{(Q-1,r)}\big{)}+\sqrt{2}\gamma\int_{0}^{h}e^{-\gamma(h-s)}d\bar{\boldsymbol{W}}_{s}$.

15:end for

In parallel

In the theorem below, we quantify the error bounds for this algorithm when
applied to constrained sampling.

######  Theorem 12.

Under Assumptions 1 and 2, we further assume that the potential $U^{\lambda}$
satisfies Condition 2.1. Let the tuning parameter
$\lambda\in\big{(}0,\frac{r}{2(p+2)}\big{)}$. Choose $\gamma$ and $h$ so that
$\gamma\geqslant 5M^{\lambda}$ and $\gamma h\leqslant 0.1\kappa^{-1/6}$, where
$\kappa=M^{\lambda}/m$. Assume that $\boldsymbol{\vartheta}_{0}^{\sf pRKLMC}$
is independent of $\mathbf{v}_{0}^{\sf pRKLMC}$ and that $\mathbf{v}_{0}^{\sf
pRKLMC}\sim\mathcal{N}_{p}(0,\gamma\mathbf{I}_{p})$. Then, for any $n\geqslant
1$, the distribution $\nu_{n}^{\textup{ \sf pRKLMC}}$ of
$\boldsymbol{\vartheta}_{n}^{\textup{\sf pRKLMC}}$ satisfies

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf pRKLMC},\nu)$ $\displaystyle\leqslant
1.8\varrho^{n}{\sf W}_{2}(\nu_{0}^{\sf
pRKLMC},\nu)+2.8\sqrt{\frac{3\mu_{2}(\nu)\lambda(2p+2)}{r}}+0.28\sqrt{\frac{\varrho^{n}\mathbb{E}[U^{\lambda}(\boldsymbol{\vartheta}_{0}^{\sf
pRKLMC})-f(\mathbf{0})]}{m}}$ $\displaystyle\qquad+44.78\sqrt{\frac{(\gamma
M^{\lambda}h^{2})^{3}}{R^{2}}+(\gamma
M^{\lambda}h^{2})^{2Q-1}}\sqrt{\frac{\kappa p}{m}}\,,$

where $\varrho=e^{-mh}$. Moreover, let the staring points
$\boldsymbol{\vartheta}_{0}^{\sf pRKLMC}=\mathbf{0}$, it holds that

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf pRKLMC},\nu)$ $\displaystyle\leqslant
1.8\varrho^{n}\sqrt{\frac{p}{m}}+\frac{3\mu_{1}(\nu)\lambda(2p+1)}{r}+44.78\sqrt{\frac{(\gamma
M^{\lambda}h^{2})^{3}}{R^{2}}+(\gamma
M^{\lambda}h^{2})^{2Q-1}}\sqrt{\frac{\kappa p}{m}}\,.$

To our knowledge, this represents the first convergence analysis for the
pRKLMC algorithm within the context of constrained sampling. When $R=1,$ we
recover the convergence rate for RKLMC [SL19, YKD23, YD24] when applied in the
context of constrained sampling. We define $M^{\lambda}=M+1/\lambda^{2},$ in
line with the projection options described in Examples 2.1 and 2.2. We then
select the optimal $\lambda$ to derive the following corollary.

######  Corollary 13.

Given the number of parallel steps $R,$ set $Q=\lceil{\log R}\rceil+2.$ Let
$\varepsilon\in(0,1)$ be a small number, $\gamma=5M^{\lambda}$ with
$M^{\lambda}=M+1/\lambda^{2}$.
(a) Set $\lambda=(\frac{7\cdot
8011}{3})^{1/8}(\frac{rh^{3}}{R\mu_{1}(\nu)(2p+1)})^{1/8}(\frac{p}{m^{2}})^{1/16}$,
choose $h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=23.6^{-8/3}p^{7/6}m^{-1}r^{7/3}\big{(}\mu_{1}(\nu)(2p+1)\big{)}^{-7/3}R^{1/3}\varepsilon^{8/3}\quad\text{and}\quad
n\geqslant\frac{1}{mh}\log\Big{(}\frac{3.6}{\varepsilon}\Big{)}\,,$

then we have ${\sf W}_{1}(\nu_{n}^{\sf
pRKLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

(b) Set $\lambda=(\frac{14\cdot
8011}{5})^{2/15}(\frac{rp}{m^{2}\mu_{2}(\nu)(2p+2)})^{1/15}(\frac{h^{3}}{R})^{2/15}$,
choose $h>0$ and $n\in\mathbb{N}$ so that

$\displaystyle
h=20.9^{-5}p^{7/3}m^{-13/6}r^{7/3}\big{(}\mu_{2}(\nu)(2p+2)\big{)}^{-7/3}R^{1/3}\varepsilon^{5}\quad\text{and}\quad
n\geqslant\frac{1}{mh}\log\Big{(}\frac{3.6}{\varepsilon}\Big{)}\,,$

then we have ${\sf W}_{2}(\nu_{n}^{\sf
pRKLMC},\nu)\leqslant\varepsilon\sqrt{p/m}$.

The corollary demonstrates that pRKLMC achieves the best performance in terms
of sample complexity compared to the other three algorithms considered in this
work.

##  4 Discussion

Strongly convexity In this work, we focus on a strongly convex potential
function $f$ coupled with specific choices of projection operators, resulting
in a strongly log-concave surrogate density $\nu^{\lambda}$. However, it is
important to note that the approximation error bound between the target
density $\nu$ and the surrogate density $\nu^{\lambda}$, as stated in
Proposition 3, depends solely on the convexity of $f$. By relaxing the strong
convexity requirement of $f$ and considering alternative projection operators,
one can generate a surrogate density $\nu^{\lambda}$ that is less restrictive
than strongly log-concave. This modification enables the extension of these
findings to a broader spectrum of sampling methods.

Other metrics The Wasserstein distance we employ is a natural metric for
measuring sampling errors due to its relevance to optimal transport theory.
However, recent advancements in gradient-based sampling have investigated
other metrics, including total variation distance, KL divergence, and
$\chi^{2}$ divergence. An intriguing avenue for future research would be to
establish error guarantees for constrained sampling concerning these
alternative metrics.

##  5 Acknowledgments

We are grateful to Arnak Dalalyan for inspiring discussions and constructive
comments on this manuscript. This work was supported by the center Hi! PARIS.

## References

  * [AC21] Kwangjun Ahn and Sinho Chewi, _Efficient constrained sampling via the mirror-langevin algorithm_ , Advances in Neural Information Processing Systems 34 (2021), 28405–28418.
  * [ADFDJ03] Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I Jordan, _An introduction to mcmc for machine learning_ , Machine learning 50 (2003), 5–43.
  * [BDMP17] Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra, _Sampling from a log-concave distribution with compact support with proximal langevin monte carlo_ , Conference on learning theory, PMLR, 2017, pp. 319–342.
  * [BEL15] Sebastien Bubeck, Ronen Eldan, and Joseph Lehec, _Finite-time analysis of projected langevin monte carlo_ , Advances in Neural Information Processing Systems 28 (2015).
  * [BEL18] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec, _Sampling from a log-concave distribution with projected langevin monte carlo_ , Discrete & Computational Geometry 59 (2018), no. 4, 757–783.
  * [Bha78] R. N. Bhattacharya, _Criteria for recurrence and existence of invariant measures for multidimensional diffusions_ , Ann. Probab. 6 (1978), no. 4, 541–553.
  * [BSU12] Marcus Brubaker, Mathieu Salzmann, and Raquel Urtasun, _A family of mcmc methods on implicitly defined manifolds_ , Artificial intelligence and statistics, PMLR, 2012, pp. 161–172.
  * [CCBJ18] Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan, _Underdamped Langevin MCMC: A non-asymptotic analysis_ , COLT 2018, 2018\.
  * [CLGL+20] Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, Philippe Rigollet, and Austin Stromme, _Exponential ergodicity of mirror-langevin diffusions_ , Advances in Neural Information Processing Systems 33 (2020), 19573–19585.
  * [Dal17] Arnak S. Dalalyan, _Theoretical guarantees for approximate sampling from a smooth and log-concave density_ , J. R. Stat. Soc. B 79 (2017), 651 – 676.
  * [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 (1991), no. 1, 1–17.
  * [DM17] Alain Durmus and Eric Moulines, _Nonasymptotic convergence analysis for the unadjusted Langevin algorithm_ , Ann. Appl. Probab. 27 (2017), no. 3, 1551–1587.
  * [DM19] Alain Durmus and Eric Moulines, _High-dimensional bayesian inference via the unadjusted langevin algorithm_.
  * [DMM19] Alain Durmus, Szymon Majewski, and Błażej Miasojedow, _Analysis of langevin monte carlo via convex optimization_ , The Journal of Machine Learning Research 20 (2019), no. 1, 2666–2711.
  * [DMP18] Alain Durmus, Eric Moulines, and Marcelo Pereyra, _Efficient bayesian computation by proximal markov chain monte carlo: when langevin meets moreau_ , SIAM Journal on Imaging Sciences 11 (2018), no. 1, 473–506.
  * [DRD20] Arnak S Dalalyan and Lionel Riou-Durand, _On sampling from a log-concave density using kinetic Langevin diffusions_ , Bernoulli 26 (2020), no. 3, 1956–1988.
  * [EGZ19] Andreas Eberle, Arnaud Guillin, and Raphael Zimmer, _Couplings and quantitative contraction rates for Langevin dynamics_ , Ann. Probab. 47 (2019), no. 4, 1982–2010.
  * [EH21] Murat A Erdogdu and Rasa Hosseinzadeh, _On the convergence of langevin monte carlo: The interplay between tail growth and smoothness_ , COLT 2021, 2021\.
  * [EHZ22] Murat A Erdogdu, Rasa Hosseinzadeh, and Shunshi Zhang, _Convergence of langevin monte carlo in chi-squared and rényi divergence_ , AISTATS 2022, 2022\.
  * [EMS18] Murat A Erdogdu, Lester Mackey, and Ohad Shamir, _Global non-convex optimization with discretized diffusions_ , NeurIPS 2018, 2018.
  * [GCSR95] Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin, _Bayesian data analysis_ , Chapman and Hall/CRC, 1995.
  * [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).
  * [GPAM+14] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio, _Generative adversarial nets_ , Advances in neural information processing systems 27 (2014).
  * [GSL92] Alan E Gelfand, Adrian FM Smith, and Tai-Ming Lee, _Bayesian analysis of constrained parameter and truncated data problems using gibbs sampling_ , Journal of the American Statistical Association 87 (1992), no. 418, 523–532.
  * [HBE20] Ye He, Krishnakumar Balasubramanian, and Murat A Erdogdu, _On the ergodicity, bias and asymptotic normality of randomized midpoint sampling method_ , Advances in Neural Information Processing Systems 33 (2020), 7366–7376.
  * [HKRC18] Ya-Ping Hsieh, Ali Kavis, Paul Rolland, and Volkan Cevher, _Mirrored langevin dynamics_ , Advances in Neural Information Processing Systems 31 (2018).
  * [KLS97] Ravi Kannan, László Lovász, and Miklós Simonovits, _Random walks and an o*(n5) volume algorithm for convex bodies_ , Random Structures & Algorithms 11 (1997), no. 1, 1–50.
  * [KLSV22] Yunbum Kook, Yin-Tat Lee, Ruoqi Shen, and Santosh Vempala, _Sampling with riemannian hamiltonian monte carlo in a constrained space_ , Advances in Neural Information Processing Systems 35 (2022), 31684–31696.
  * [KVZ24] Yunbum Kook, Santosh S Vempala, and Matthew S Zhang, _In-and-out: Algorithmic diffusion for sampling convex bodies_ , arXiv preprint arXiv:2405.01425 (2024).
  * [KW13] Diederik P Kingma and Max Welling, _Auto-encoding variational bayes_ , arXiv preprint arXiv:1312.6114 (2013).
  * [Lam21] Andrew Lamperski, _Projected stochastic gradient langevin algorithms for constrained sampling and non-convex learning_ , Conference on Learning Theory, PMLR, 2021, pp. 2891–2937.
  * [LBGH22] Zhou Lu, Nataly Brukhim, Paula Gradu, and Elad Hazan, _Projection-free adaptive regret with membership oracles_ , arXiv preprint arXiv:2211.12638 (2022).
  * [Leh23] Joseph Lehec, _The langevin monte carlo algorithm in the non-smooth log-concave case_ , The Annals of Applied Probability 33 (2023), no. 6A, 4858–4874.
  * [LLK+22] Lingxiao Li, Qiang Liu, Anna Korba, Mikhail Yurochkin, and Justin Solomon, _Sampling with mollified interaction energy descent_ , arXiv preprint arXiv:2210.13400 (2022).
  * [Lov99] László Lovász, _Hit-and-run mixes fast_ , Mathematical programming 86 (1999), 443–461.
  * [LS93] László Lovász and Miklós Simonovits, _Random walks in a convex body and an improved volume algorithm_ , Random structures & algorithms 4 (1993), no. 4, 359–412.
  * [LV07] László Lovász and Santosh Vempala, _The geometry of logconcave functions and sampling algorithms_ , Random Structures & Algorithms 30 (2007), no. 3, 307–358.
  * [MCC+21] Yi-An Ma, Niladri S. Chatterji, Xiang Cheng, Nicolas Flammarion, Peter L. Bartlett, and Michael I. Jordan, _Is there an analog of Nesterov acceleration for gradient-based MCMC?_ , Bernoulli 27 (2021), no. 3, 1942 – 1992.
  * [MFWB22] Wenlong Mou, Nicolas Flammarion, Martin J Wainwright, and Peter L Bartlett, _Improved bounds for discretization of Langevin diffusions: Near-optimal rates without convexity_ , Bernoulli 28 (2022), no. 3, 1577–1601.
  * [Mha22] Zakaria Mhammedi, _Efficient projection-free online convex optimization with membership oracle_ , Conference on Learning Theory, PMLR, 2022, pp. 5314–5390.
  * [MHFH+23] Alireza Mousavi-Hosseini, Tyler Farghly, Ye He, Krishnakumar Balasubramanian, and Murat A Erdogdu, _Towards a Complete Analysis of Langevin Monte Carlo: Beyond Poincaré Inequality_ , arXiv preprint arXiv:2303.03589 (2023).
  * [OHHDT22] Joaquim Ortiz-Haro, Jung-Su Ha, Danny Driess, and Marc Toussaint, _Structured deep generative models for sampling on constraint manifolds in sequential manipulation_ , Conference on Robot Learning, PMLR, 2022, pp. 213–223.
  * [Per16] Marcelo Pereyra, _Proximal markov chain monte carlo algorithms_ , Statistics and Computing 26 (2016), 745–760.
  * [RRT17] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky, _Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis_ , COLT 2017, 2017.
  * [RT96] G. O. Roberts and R. L. Tweedie, _Exponential convergence of Langevin distributions and their discrete approximations_ , Bernoulli 2 (1996), no. 4, 341–363. MR 1440273 (98j:62014)
  * [RW09] R Tyrrell Rockafellar and Roger J-B Wets, _Variational analysis_ , vol. 317, Springer Science & Business Media, 2009.
  * [SL19] Ruoqi Shen and Yin Tat Lee, _The randomized midpoint method for log-concave sampling_ , Advances in Neural Information Processing Systems 32 (2019).
  * [Smi84] Robert L Smith, _Efficient monte carlo procedures for generating points uniformly distributed over bounded regions_ , Operations Research 32 (1984), no. 6, 1296–1308.
  * [SR20] Adil Salim and Peter Richtárik, _Primal dual interpretation of the proximal stochastic gradient langevin algorithm_ , Advances in Neural Information Processing Systems 33 (2020), 3786–3796.
  * [Stu10] Andrew M Stuart, _Inverse problems: a bayesian perspective_ , Acta numerica 19 (2010), 451–559.
  * [Vil09] Cédric Villani, _Optimal transport: old and new_ , vol. 338, Springer, 2009\.
  * [YD24] Lu Yu and Arnak Dalalyan, _Parallelized midpoint randomization for langevin monte carlo_ , arXiv preprint arXiv:2402.14434 (2024).
  * [YKD23] Lu Yu, Avetik Karagulyan, and Arnak Dalalyan, _Langevin monte carlo for strongly log-concave distributions: Randomized midpoint revisited_ , arXiv preprint arXiv:2306.08494 (2023).
  * [ZCL+23] Matthew Zhang, Sinho Chewi, Mufan Bill Li, Krishnakumar Balasubramanian, and Murat A Erdogdu, _Improved Discretization Analysis for Underdamped Langevin Monte Carlo_ , arXiv preprint arXiv:2302.08049 (2023).
  * [ZPFP20] Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra, _Wasserstein control of mirror langevin monte carlo_ , Conference on Learning Theory, PMLR, 2020, pp. 3814–3841.

##  Appendix A Proofs of Section 2

In this part, we provide the proof of Lemma 1, Lemma 2, Proposition 3, and
Proposition 4.

###  A.1 Proofs of Lemma 1 and Lemma 2

Proof. [Proof of Lemma 1] By Theorem 2.26 in [RW09], the function
$\ell_{\mathcal{K}}^{E,\lambda}(x)$ is convex and continuously differentiable,
and $\ell_{\mathcal{K}}^{E,\lambda}(x)$ is $\frac{1}{\lambda^{2}}$-smooth. The
last claim follows from Lemma 1 in [BDMP17].

∎

Proof. [Proof of Lemma 2] By the definition of
$\ell_{\mathcal{K}}^{G,\lambda},$ we have

$\displaystyle\ell_{\mathcal{K}}^{G,\lambda}=\frac{1}{2\lambda^{2}}\big{(}1-\frac{1}{g_{\mathcal{K}}(x)}\big{)}^{2}\,\\!|\\!|x|\\!|_{2}^{2}\,.$

By the definition of $U^{G,\lambda}$, the function $U^{G,\lambda}$ is
$m$-strongly convex, continuously differentiable. Moreover, it holds that

$\displaystyle\\!|\\!|\nabla\ell_{\mathcal{K}}^{G,\lambda}(x)-\nabla\ell_{\mathcal{K}}^{G,\lambda}(y)|\\!|_{2}\leqslant\frac{1}{\lambda^{2}}\,\\!|\\!|x-y|\\!|_{2}\,,$

This implies $U^{G,\lambda}$ is $\frac{1}{\lambda^{2}}+M$-smooth. Now we show
that $\int_{\mathbb{R}^{p}}e^{-U^{G,\lambda}(s)}{\rm d}s<\infty.$ Let
$R:=\max_{x\in\mathcal{K}}\\!|\\!|x|\\!|_{2},$ and there is a constant
$l\in\mathbb{R}$ such that $f(x)\geqslant l,\forall x\in\mathbb{R}.$ Define
$R:=\sup_{x\in\mathcal{K}}\\!|\\!|x|\\!|_{2}.$ Note that

$\displaystyle\int_{\mathbb{R}^{p}}e^{-U^{G,\lambda}(x)}{\rm d}x=$
$\displaystyle\int_{\mathcal{K}}e^{-f(x)}{\rm
d}x+\int_{\mathcal{B}_{2}(0,R)\cap\mathcal{K}^{c}}e^{-U^{G,\lambda}(x)}{\rm
d}x+\int_{\mathcal{B}_{2}(0,R)^{c}}e^{-U^{G,\lambda}(x)}{\rm d}x\,.$

We now derive the bound for the second and the third terms on the right-hand
side of the previous display. Since $f$ is convex, there exists
$a\in\mathbb{R}$ and $b\in\mathbb{R}^{p}$ such that

$\displaystyle f(x)\geqslant a+\langle b,\,x\rangle\,.$

We then find

$\displaystyle\int_{\mathcal{B}_{2}(0,R)\cap\mathcal{K}^{c}}e^{-U^{G,\lambda}(x)}{\rm
d}x\leqslant e^{-a+bR}\int_{\mathcal{B}_{2}(0,R)}{\rm d}x\,.$

Thus,

$\displaystyle\int_{\mathcal{B}_{2}(0,R)\cap\mathcal{K}^{c}}e^{-U^{G,\lambda}(x)}{\rm
d}x\leqslant$ $\displaystyle e^{-a+bR}\int_{\mathcal{B}_{2}(0,R)}{\rm d}x$
$\displaystyle\leqslant$ $\displaystyle
e^{-a+bR}\int_{\mathcal{B}_{2}(0,R)}{\rm d}x$ $\displaystyle=$ $\displaystyle
e^{-a+bR}\frac{\pi^{p/2}}{\Gamma(\frac{d}{2}+1)}R^{p}$ $\displaystyle<$
$\displaystyle\infty\,.$

Note that

$\displaystyle\int_{\mathcal{B}_{2}(0,R)^{c}}e^{-U^{G,\lambda}(x)}{\rm
d}x\leqslant
e^{-l}\int_{\mathcal{B}_{2}(0,R)^{c}}e^{-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-R)^{2}}{\rm
d}x$

By Fubini’s theorem and the fact that

$\displaystyle
e^{-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-R)^{2}}=\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-\frac{t^{2}}{2\lambda^{2}}}\mathds{1}_{[\,\\!|\\!|x|\\!|_{2}-R,\infty)}(t){\rm
d}t\,,$

we obtain

$\displaystyle\int_{\mathcal{B}_{2}(0,R)^{c}}e^{-U^{G,\lambda}(x)}{\rm
d}x\leqslant$ $\displaystyle
e^{-l}\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-\frac{t}{2\lambda^{2}}}\int_{x:\,\\!|\\!|x|\\!|_{2}\geq
R}\mathds{1}_{[\,\\!|\\!|x|\\!|_{2}-R,\infty)}(t){\rm d}x{\rm d}t$
$\displaystyle\leqslant$ $\displaystyle
e^{-l}\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-\frac{t}{2\lambda^{2}}}\int_{\mathcal{B}_{2}(0,R+t)}{\rm
d}x{\rm d}t$ $\displaystyle=$ $\displaystyle
e^{-l}\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-\frac{t}{2\lambda^{2}}}\frac{\pi^{p/2}}{\Gamma(\frac{p}{2}+1)}({R+t})^{d}{\rm
d}t$ $\displaystyle<$ $\displaystyle\infty\,.$

Collecting pieces gives the desired result.

∎

###  A.2 Proofs of Proposition 3

To prove Proposition 3, we need the following auxiliary lemma.

######  Lemma 14.

Assume $f$ is convex and
$\mathbf{0}=\mathop{\mathrm{missing}}{argmin}_{\boldsymbol{\theta}\in\mathcal{K}}f(\boldsymbol{\theta}).$
Under Assumption 1 , it holds for any
$\lambda\in\big{(}0,\frac{r}{2(p+k)}\big{)}$ that

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}\nu^{\lambda}(x){\rm
d}x\leqslant\frac{3\lambda\,(p+k)}{r}\times\mu_{k}(\nu)\,,$

where
$\mu_{k}(\nu)=\int_{\mathbb{R}^{p}}\\!|\\!|\boldsymbol{\theta}|\\!|_{2}^{k}\,\nu({\rm
d}\boldsymbol{\theta}).$

Proof. [Proof of Lemma 14] By Fubini’s theorem and the fact that

$\displaystyle\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-\frac{t^{2}}{2\lambda^{2}}}\mathds{1}_{[\,\\!|\\!|x-\mathtt{P}_{\mathcal{K}}(x)|\\!|_{2},+\infty)}(t){\rm
d}t=e^{-\frac{1}{2\lambda^{2}}\,\\!|\\!|x-\mathtt{P}_{\mathcal{K}}(x)|\\!|_{2}^{2}}\,,$

it holds for every $k\geqslant 0$ that

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-U^{\lambda}(x)}{\rm
d}x=\int_{0}^{\infty}\Big{(}\underbrace{\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}\mathds{1}(\,\\!|\\!|x-\mathtt{P}_{\mathcal{K}}(x)|\\!|_{2}\leqslant
t)\,{\rm
d}x}_{\psi(t)}\,\Big{)}\frac{t}{\lambda^{2}}e^{-\frac{t^{2}}{2\lambda^{2}}}\,{\rm
d}t.$ (A.1)

Define the set
$\mathcal{K}_{t}:=\\{x\in\mathbb{R}^{p}:\\!|\\!|x-\mathtt{P}_{\mathcal{K}}(x)|\\!|_{2}\leqslant
t\\},$ we then have

$\displaystyle\psi(t)=\int_{\mathcal{K}_{t}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}\,{\rm
d}x-\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}\,{\rm d}x.$ (A.2)

On the other hand, it holds for every $b>0$ that

$\displaystyle\int_{\mathcal{K}_{t}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}\,{\rm d}x$
$\displaystyle=(1+b)^{p+k}\int_{\mathcal{K}_{t}/(1+b)}\\!|\\!|y|\\!|_{2}^{k}e^{-f((1+b)y)}\,{\rm
d}y$
$\displaystyle\leqslant(1+b)^{p+k}\int_{\mathcal{K}_{t}/(1+b)}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y$

where we have used the inequality

$\displaystyle
f((1+b)y)+bf(0)\geqslant(1+b)f\Big{(}\frac{1}{1+b}\,(1+b)y+\frac{b}{b+1}\,0\Big{)}=(1+b)f(y)\geqslant
f(y)+bf(0),$

implying that $f((1+b)y)\geqslant f(y)$ for every $b>0$ and every
$y\in\mathbb{R}^{p}$. In addition, it holds for any
$x\in\mathcal{K}_{t}/(1+b)$ that

$\displaystyle\\!|\\!|(1+b)x-\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}|\\!|_{2}\leqslant
t$

which is equivalent to

$\displaystyle\frac{1+b}{b}\,\\!|\\!|x-\frac{1}{1+b}\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}|\\!|_{2}\leqslant\frac{t}{b}\,.$

Note that we can rewrite $x$ as

$\displaystyle
x=\frac{1}{1+b}\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}+\frac{b}{1+b}\frac{x-\frac{1}{1+b}\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}}{\frac{b}{1+b}}\,.$

Notice that

$\displaystyle\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}\in\mathcal{K}$

and when $b=t/r$ we have

$\displaystyle\frac{x-\frac{1}{1+b}\mathtt{P}_{\mathcal{K}}\big{(}(1+b)x\big{)}}{\frac{b}{1+b}}\in\mathcal{K}\,.$

By the convexity of the set $\mathcal{K}$, this implies $x\in\mathcal{K}.$
Therefore, we obtain

$\displaystyle\int_{\mathcal{K}_{t}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}\,{\rm d}x$
$\displaystyle\leqslant\Big{(}1+\frac{t}{r}\Big{)}^{p+k}\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y\leqslant
e^{(p+k)t/r}\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm d}y\,.$

Combining this inequality with (A.2), we get

$\displaystyle\psi(t)\leqslant\big{(}e^{(p+k)t/r}-1\big{)}\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y\,.$

Substituting this upper bound in (A.1) leads to

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-U^{\lambda}(x)}{\rm
d}x$
$\displaystyle\leqslant\int_{0}^{\infty}\frac{t}{\lambda^{2}}e^{-{t^{2}}/{2\lambda^{2}}}\big{(}e^{(p+k)t/r}-1\big{)}\,{\rm
d}t\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm d}y$
$\displaystyle=\int_{0}^{\infty}\big{(}e^{(p+k)t/r}-1\big{)}\,{\rm
d}\big{(}-e^{-{t^{2}}/{2\lambda^{2}}}\big{)}\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y$
$\displaystyle=\frac{(p+k)}{r}\int_{0}^{\infty}\exp\Big{\\{}-\frac{t^{2}}{2\lambda^{2}}+\frac{(p+k)t}{r}\Big{\\}}\,{\rm
d}t\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm d}y$
$\displaystyle=\frac{(p+k)}{r}\exp\Big{\\{}\frac{\lambda^{2}(p+k)^{2}}{2r^{2}}\Big{\\}}\int_{0}^{\infty}\exp\Big{\\{}-\frac{(t-(\lambda^{2}/r)(p+k))^{2}}{2\lambda^{2}}\Big{\\}}\,{\rm
d}t\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm d}y$
$\displaystyle\leqslant\frac{\sqrt{2\pi}\lambda\,(p+k)}{r}\exp\Big{\\{}\frac{\lambda^{2}(p+k)^{2}}{2r^{2}}\Big{\\}}\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y\,.$

When $\lambda<\frac{r}{2(p+k)}$, it holds that

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-U^{\lambda}(x)}{\rm
d}x$
$\displaystyle\leqslant\frac{3\lambda\,(p+k)}{r}\times\int_{\mathcal{K}}\\!|\\!|y|\\!|_{2}^{k}e^{-f(y)}\,{\rm
d}y.$ (A.3)

Finally, note that the normalizing constant of $\nu^{\lambda}$ can be lower
bounded as follows

$\displaystyle\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(x)}{\rm
d}x\geqslant\int_{\mathcal{K}}e^{-U^{\lambda}(x)}{\rm
d}x=\int_{\mathcal{K}}e^{-f(x)}{\rm d}x\,.$

The desired result follows readily by dividing the two sides of (A.3) by
$\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(x)}{\rm d}x$ and employing this
inequality.

∎

We are now ready to prove Proposition 3

Proof. [Proof of Proposition 3] By Theorem 6.15 in [Vil09], it holds that

$\displaystyle{\sf
W}_{q}^{q}(\nu,\nu^{\lambda})\leqslant\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}|\nu(x)-\nu^{\lambda}(x)|{\rm
d}x\,.$

We note that

$\displaystyle\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}|\nu(x)-\nu^{\lambda}(x)|\,{\rm
d}x$
$\displaystyle=\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}|\nu(x)-\nu^{\lambda}(x)|{\rm
d}x+\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}(x){\rm d}x$
$\displaystyle=\Big{(}1-\frac{\int_{\mathcal{K}}e^{-f(y)}\,{\rm
d}y}{\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(y)}\,{\rm
d}y}\Big{)}\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}\nu(x){\rm
d}x+\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}(x){\rm d}x$
$\displaystyle=\frac{\int_{\mathcal{K}^{c}}e^{-U^{\lambda}(y)}\,{\rm
d}y}{\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(y)}\,{\rm
d}y}\times\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}\nu(x){\rm
d}x+\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}(x){\rm d}x$
$\displaystyle=\int_{\mathcal{K}^{c}}\nu^{\lambda}(y)\,{\rm
d}y\times\mu_{q}(\nu)+\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}(x){\rm
d}x\,.$

When $\lambda\leqslant\frac{r}{2(p+q)}$, we apply Lemma 14 with $k=0$ and
$k=q$, which leads us to the inequality

$\displaystyle\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}|\nu(x)-\nu^{\lambda}(x)|\,{\rm
d}x\leqslant\frac{3{\lambda}\big{(}p\mu_{0}(\nu)+(p+q)\big{)}\mu_{q}(\nu)}{r}=\frac{3{\lambda}(2p+q)}{r}\times\mu_{q}(\nu).$

This completes the proof.

∎

###  A.3 Proof of Proposition 4

To prove Proposition 4, we need the following auxiliary lemma.

######  Lemma 15.

Under the assumptions stated in Proposition 4, for any $k\geqslant 1$ and any
$\lambda\in\big{(}0,\frac{r}{2(p+k)}\wedge\frac{1}{\sqrt{M}}\big{)},$ it holds
that

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}\nu^{\lambda}({\rm
d}x)/\mu^{k}(\nu)\geqslant\frac{1}{2}\sqrt{\frac{\pi}{2}}e^{-\frac{3}{4}Mr^{2}}\frac{p+k}{r+3\lambda
p}\lambda\,.$

Proof. [Proof of Lemma 15] By the definition of $\nu^{\lambda}$ and
$\mu^{k}(\nu),$ we have

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}\nu^{\lambda}({\rm
d}x)/\mu^{k}(\nu)$
$\displaystyle=\frac{\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-r)^{2}}{\rm
d}x}{\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(x)}{\rm
d}x}\cdot\frac{\int_{\mathcal{K}}e^{-f(x)}{\rm
d}x}{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}{\rm d}x}\,.$ (A.4)

Note that $f$ is $M$-smooth and $f(\mathbf{0})=0$, it then holds that

$\displaystyle f(x)\leqslant\frac{M}{2}\\!|\\!|x|\\!|_{2}^{2}\,.$

This implies

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-r)^{2}}{\rm
d}x\geqslant\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-\frac{M}{2}\\!|\\!|x|\\!|_{2}^{2}-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-r)^{2}}{\rm
d}x\,.$

Employing the integration in polar coordinates gives

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-\frac{M}{2}\\!|\\!|x|\\!|_{2}^{2}-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-r)^{2}}{\rm
d}x$
$\displaystyle=\omega_{p-1}\int_{r}^{+\infty}t^{k+p-1}e^{-\frac{M}{2}t^{2}-\frac{1}{2\lambda^{2}}(t-r)^{2}}{\rm
d}t\,,$

where
$\omega_{p-1}=\int_{-\pi}^{\pi}\int_{0}^{\pi}\cdots\int_{0}^{\pi}\prod_{j=1}^{p-2}(\sin\theta_{j})^{p-j-1}{\rm
d}\theta_{1}\cdots{\rm d}\theta_{p-1}.$ Let $t=r+\lambda z,$ it then follows
that

$\displaystyle\int_{\mathcal{K}^{c}}$
$\displaystyle\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}-r)^{2}}{\rm
d}x$ $\displaystyle\geqslant\lambda\omega_{p-1}\int_{0}^{+\infty}(r+\lambda
z)^{k+p-1}e^{-\frac{M}{2}(r+\lambda z)^{2}-\frac{1}{2}z^{2}}{\rm d}z$
$\displaystyle\geqslant\lambda\omega_{p-1}r^{k+p-1}\int_{0}^{+\infty}e^{-\frac{M}{2}(r+\lambda
z)^{2}-\frac{1}{2}z^{2}}{\rm d}z$
$\displaystyle=\lambda\omega_{p-1}r^{k+p-1}\int_{0}^{+\infty}e^{-\big{(}\sqrt{\frac{1+M\lambda^{2}}{2}}x+\frac{M\lambda
r}{\sqrt{2(1+M\lambda^{2})}}\big{)}^{2}}{\rm
d}xe^{-\frac{Mr^{2}}{2}\big{(}1+\frac{M\lambda^{2}}{1+M\lambda^{2}}\big{)}}\,.$

The last step follows from a change of variables. Using the fact that
$\lambda<1/\sqrt{M}$ and $(x+y)^{2}\leqslant 2x^{2}+2y^{2},\forall
x,y\geqslant 0,$ we then obtain

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)-\frac{1}{2\lambda^{2}}(\,\\!|\\!|x|\\!|_{2}^{2}-r)^{2}}{\rm
d}x$
$\displaystyle\geqslant\lambda\omega_{p-1}r^{k+p-1}\int_{0}^{+\infty}e^{-(1+M\lambda^{2})x^{2}}{\rm
d}xe^{-\frac{Mr^{2}}{2}(1+\frac{M\lambda^{2}}{1+M\lambda^{2}})}$
$\displaystyle=\lambda\omega_{p-1}r^{k+p-1}\frac{1}{2}\sqrt{\frac{\pi}{1+M\lambda^{2}}}e^{-\frac{Mr^{2}}{2}(1+\frac{M\lambda^{2}}{1+M\lambda^{2}})}$
$\displaystyle\geqslant\lambda\omega_{p-1}r^{k+p-1}\frac{1}{2}\sqrt{\frac{\pi}{2}}e^{-\frac{3}{4}Mr^{2}}\,.$

Moreover, by display (A.3), we find

$\displaystyle\int_{\mathbb{R}^{p}}e^{-U^{\lambda}(x)}{\rm d}x$
$\displaystyle=\int_{\mathcal{K}}e^{-U^{\lambda}(x)}{\rm
d}x+\int_{\mathcal{K}^{c}}e^{-U^{\lambda}(x)}{\rm d}x$
$\displaystyle\leqslant\Big{(}1+\frac{3\lambda
p}{r}\Big{)}\int_{\mathcal{K}}e^{-\,\\!|\\!|x|\\!|_{2}^{2}/2}{\rm d}x\,.$

Furthermore, using the integration in polar coordinates again gives

$\displaystyle\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{k}e^{-f(x)}{\rm d}x$
$\displaystyle\leqslant\omega_{p-1}\int_{0}^{r}t^{k+p-1}{\rm d}t$
$\displaystyle=\omega_{p-1}\frac{r^{k+p}}{k+p}\,.$

Collecting pieces and plugging them into display (A.4) gives

$\displaystyle\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{k}\nu^{\lambda}({\rm
d}x)/\mu^{k}(\nu)\geqslant\frac{1}{2}\sqrt{\frac{\pi}{2}}e^{-\frac{3}{4}Mr^{2}}\frac{p+k}{r+3\lambda
p}\lambda$

as desired.

∎

We are now ready to prove Proposition 4.

Proof. [Proof of Proposition 4] By Proposition 7.29 in [Vil09], it holds for
any $q\geqslant 1$ that

$\displaystyle{\sf
W}_{q}^{q}(\nu^{\lambda},\nu)\geqslant\Big{|}\Big{(}\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)\Big{)}^{1/q}-\Big{(}\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu({\rm
d}x)\Big{)}^{1/q}\Big{|}^{q}\,.$

Dividing both sides by $\mu_{q}(\nu)$ gives

$\displaystyle{\sf W}_{q}^{q}(\nu^{\lambda},\nu)/\mu_{q}(\nu)$
$\displaystyle\geqslant\Big{|}\Big{(}\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)/\mu_{q}(\nu)\Big{)}^{1/q}-1\Big{|}^{q}$
$\displaystyle=\Bigg{|}\Bigg{(}\frac{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)+\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)}{\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu({\rm
d}x)}\Bigg{)}^{1/q}-1\Bigg{|}\,.$

By the definitions of $\nu$ and $\nu^{\lambda},$ we have

$\displaystyle\frac{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)}{\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu(dx)}$
$\displaystyle=\frac{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}e^{-f(x)}{\rm
d}x}{\int_{\mathbb{R}^{p}}e^{-U^{\lambda}}dx}\cdot\frac{\int_{\mathcal{K}}e^{-f(x)}{\rm
d}x}{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}e^{-f(x)}{\rm d}x}$
$\displaystyle=\frac{\int_{\mathcal{K}}e^{-f(x)}{\rm
d}x}{\int_{\mathcal{K}}e^{-f(x)}{\rm d}x+\int_{\mathcal{K}^{c}}e^{-f(x)}{\rm
d}x}\,.$

Combining this with display (A.3) then yields

$\displaystyle\frac{\int_{\mathcal{K}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)}{\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu({\rm
d}x)}\geqslant\frac{r}{r+3\lambda p}\,.$

By Lemma 15 and Lemma 14, we have

$\displaystyle\frac{1}{2}\sqrt{\frac{\pi}{2}}e^{-\frac{3}{4}Mr^{2}}\frac{p+q}{r+3\lambda
p}\lambda\leqslant\frac{\int_{\mathcal{K}^{c}}\\!|\\!|x|\\!|_{2}^{q}\nu^{\lambda}({\rm
d}x)}{\int_{\mathbb{R}^{p}}\\!|\\!|x|\\!|_{2}^{q}\nu({\rm
d}x)}\leqslant\frac{3\lambda(p+q)}{r}\,.$

Collecting pieces gives the desired result.

∎

##  Appendix B Proofs of Section 3

Proof. [Proof of Theorem 6] The triangle inequality for the Wasserstein
distance provides us with

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu)\leqslant{\sf
W}_{2}(\nu_{n}^{\sf LMC},\nu^{\lambda})+{\sf W}_{2}(\nu^{\lambda},\nu)\,.$
(B.1)

By Theorem 9 in [DMM19], it holds that

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu^{\lambda})\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu^{\lambda})+\sqrt{\frac{2M^{\lambda}ph}{m}}$ (B.2)

provided that $h\leqslant 1/M^{\lambda}.$ Using the triangle inequality for
the Wasserstein distance again gives

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu^{\lambda})\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu)+\sqrt{\frac{2M^{\lambda}ph}{m}}+{\sf W}_{2}(\nu,\nu^{\lambda})\,.$

Plugging this back to display (B.1) then provides us with

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu)\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu)+\sqrt{\frac{2M^{\lambda}ph}{m}}+2{\sf W}_{2}(\nu,\nu^{\lambda})\,.$

Combining this with Proposition 3, for any
$\lambda\in\big{(}0,\frac{r}{2(p+2)}\big{)},$ we obtain

$\displaystyle{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu)\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu)+\sqrt{\frac{2M^{\lambda}ph}{m}}+2\Big{(}\frac{3\mu_{2}(\nu)\lambda(2p+2)}{r}\Big{)}^{\frac{1}{2}}\,.$

This completes the proof of the first claim.

We then proceed to prove the second claim. By the triangle inequality and the
monotonicity of Wasserstein distance, we have

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf LMC},\nu)$ $\displaystyle\leqslant{\sf
W}_{1}(\nu_{n}^{\sf LMC},\nu^{\lambda})+{\sf W}_{1}(\nu^{\lambda},\nu)$
$\displaystyle\leqslant{\sf W}_{2}(\nu_{n}^{\sf LMC},\nu^{\lambda})+{\sf
W}_{1}(\nu^{\lambda},\nu)\,.$

When $\nu_{0}^{\sf LMC}$ is the Dirac mass at the minimizer of the function
$f$, by Proposition 1 in [DM19], it holds that

$\displaystyle{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu^{\lambda})\leqslant\sqrt{\frac{p}{m}}\,.$

Combining this with previous display, display (B.2), and Proposition 3 then
gives

$\displaystyle{\sf W}_{1}(\nu_{n}^{\sf LMC},\nu)$ $\displaystyle\leqslant
e^{-\frac{mnh}{2}}{\sf W}_{2}(\nu_{0}^{\sf
LMC},\nu^{\lambda})+\sqrt{\frac{2M^{\lambda}ph}{m}}+{\sf
W}_{1}(\nu^{\lambda},\nu)$ $\displaystyle\leqslant
e^{-\frac{mnh}{2}}\sqrt{\frac{p}{m}}+\sqrt{\frac{2M^{\lambda}ph}{m}}+\frac{3\mu_{1}(\nu)\lambda(2p+1)}{r}$

as desired.

∎

Proof. [Proof of Theorem 8] The proof employs a similar technique to that used
in the proof of Theorem 6. The key distinction lies in the utilization of
Theorem 1 from [YD24] to bound the ${\sf W}_{2}$ distance between
$\nu_{n}^{\sf pRLMC}$ and $\nu^{\lambda}$.

∎

Proof. [Proof of Theorem 10] The proof employs a similar technique to that
used in the proof of Theorem 6. The key distinction lies in the utilization of
Theorem 3 from [YKD23] to bound the ${\sf W}_{2}$ distance between
$\nu_{n}^{\sf KLMC}$ and $\nu^{\lambda}$.

∎

Proof. [Proof of Theorem 12] The proof employs a similar technique to that
used in the proof of Theorem 6. The key distinction lies in the utilization of
Theorem 2 from [YD24] to bound the ${\sf W}_{2}$ distance between
$\nu_{n}^{\sf pRKLMC}$ and $\nu^{\lambda}$.

∎