Riemannian Proximal Sampler for High-accuracy Sampling on Manifolds
Authors: Yunrui Guan
Krishna Balasubramanian
Shiqian Ma
Venue: NeurIPS 2025 poster ()
Keywords: High-accuracy sampling
Riemannian manifolds
Proximal Sampler
Heat-kernels
Primary area: theory

Abstract: We introduce the \textit{Riemannian Proximal Sampler}, a method for sampling from densities defined on Riemannian manifolds. The performance of this sampler critically depends on two key oracles: the \textit{Manifold Brownian Increments (MBI)} oracle and the \textit{Riemannian Heat-kernel (RHK)} oracle. We establish high-accuracy sampling guarantees for the Riemannian Proximal Sampler, showing that generating samples with \(\varepsilon\)-accuracy requires \(\mathcal{O}(\log(1/\varepsilon))\) iterations in Kullback-Leibler divergence assuming access to exact oracles and \(\mathcal{O}(\log^2(1/\varepsilon))\) iterations in the total variation metric assuming access to sufficiently accurate inexact oracles. Furthermore, we present practical implementations of these oracles by leveraging heat-kernel truncation and Varadhan’s asymptotics. In the latter case, we interpret the Riemannian Proximal Sampler as a discretization of the entropy-regularized Riemannian Proximal Point Method on the associated Wasserstein space. We provide preliminary numerical results that illustrate the effectiveness of the proposed methodology.

Riemannian Proximal Sampler for High-accuracy Sampling on Manifolds
Yunrui Guan
Department of Computational Applied Mathematics and Operations Research Rice University Houston, TX 77005 yg83@rice.edi
Krishnakumar Balasubramanian
Department of Statistics University of California, Davis
Davis, CA 95616
kbala@ucdavis.. edu
Shiqian Ma
Department of Computational Applied Mathematics and Operations Research Rice University Houston, TX 77005
sqma@rice.edd
Abstract
We introduce the Riemannian Proximal Sampler, a method for sampling from den- sities defined on Riemannian manifolds. The performance of this sampler critically depends on two key oracles: the Manifold Brownian Increments (MBI) oracle and the Riemannian Heat-kernel (RHK) oracle. We establish high-accuracy sampling guarantees for the Riemannian Proximal Sampler, showing that generating samples with g-accuracy requires (loo(1/E)) iterations in Kullback-Leible divergence assuming access to exact oracles and Oloog" (1/e)) metric assuming access to sufficiently accurate inexact oracles. Furthermore, we present practical implementations of these oracles by leveraging heat-kernel trun- cation and Varadhan's 's asymptotics. In the latter case, we interpret the Riemannian Proximal Sampler as a discretization of the cetropy-reeularing mal Point Method on the associated Wasserstein space. We provide preliminary numerical results that illustrate the effectiveness of the proposed methodology.
(1.e)())e) iterations in the total variation
Riemannian Proxi-
Introduction
We study the problem sampling from a density x x e on a Riemannian manifold ( (M,g). where g is the metric and x is defined with respect to the volume measure dV, The normalization constant e -5 dV, is unknown. Riemannian sampling arises in Bayesian inference (e.g., hierarchical models and Bayesian deep learning) [Girolami and Calderhead 2011 Byrne and Girolami 2013
Patterson and Teh. 2013 Liu and Zhu 2018 Arnaudon et al. 2019 Liu al.) [016 Piggott and Solo, 2016 2 Muniz et 2002 Lie et al.] 2023], statistical physics (e.g., constrained molecular dynamics) Leimkuhler and Matthews] 2016]. and manifold optimization problems such as eigenvalue
computation, low-rank approximation, and diffusion models [Goyal and Shetty 2019 Li and Erdogdu
2023
Yu et al. 2023 Bonet et al. 2023 De Bortoli et al. 2022 Huang et al. 2022
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Table I: Comparison of iteration complexity with Li and Erdogdu [2023] and Cheng et al. 2022] Here, e the target accuracy. K2 the gradient Lipschitz constant, and d' is the manifold dimension; for product of n spheres of dimension m, d = mn. All works assume a uniform lower bound on
Ricci curvature and include the hypersphere as a common example. [Cheng al. [022] does not provide explicit dependence on d or other problem parameters.
Assumptions
Source
Complexity
Metric
LSI, M STh 8th x I is Ky-smooth
8m
Li and Erdogdu [2033
O
dR2 a
log
(po)
KL
Distant- Distant-dissipativig dissipativity, f
M has bounded sectional curvature
Cheng et al. [2022]
o ( Poly (d, K2))
W1
LSI, M =5 f L1 Lipschiii Corollary 10
1.3d Ö log2? a
TV
On a Riemannian manifold, Langevin dynamics takes the form dX,- --,=-grad -gad f(X;) di + v2d where grad is the Riemannian gradient and B, is Brownian motion on the manifold. This formu- lation extends Euclidean Langevin dynamics by incorporating geometric information through the Riemannian metric, but discretizing manifold Brownian motion generally intractable. Li and Erdogdu (023) considered product manifolds 5 X ... X Sth showing convergence for scheme that discretizes only the drift while assuming exact Brownian motion-feasibl on spheres. Gaatmiyy and Vempala [022 extended this to general Hessian manifolds, however, requiring exact Brownian motion. Both assume log-Sobolev inequality and achieve poly(1/a iteration complexity in KL divergence, but the need for exact Brownian motion in Gatmiry and Vempala [2022] limits practical applicability.
Cheng et al. [2022] analyzed a practical discretization of Riemannian Langevin dynamics, where both drift and noise are discretized, They established Ô(1 (s) iteration complexity in the I-Wasserstein distance under general assumptions, and in the 2-Wasserstein distance under stronger, log-concavity- like condition. A key challenge is proving Wasserstein contractivity without convexity (c.g., on compact manifolds), addressed via second-order expansion of the Jacobi equation [Cheng al. 2022, Lemma 29]. Kong and Tao [024] proposed Lic-group MCMC sampler for densities on Lie groups, achieving polynomial ((poly(1/E)) complexity in 2- Wasserstein distance.
In comparison to the above works for Riemannian sampling, for the Euclidean case, high-accuracy algorithmss, i.e., algorithms with iteration complexity of (ppolylog((/e) are available under various assumptions (that are essentially based on (strong) log-concavity or isoperimetry): see for example
ett al. [2021], Chen et al. [2022). |aan et al. [2023] He et al. [2024] for such results for the Euclidean
proximal sampler and Dwivedi et al. 2019]]
Chen et al. 22020]. Chewi et al. (2021), Lee et al. 2020).
Wu et al. 2002 Chen and Gatmiry 2023 Andrieu al. [2024], , Altschuler and Chewi [0024] for
various Metropolized algorithms including Metropolis Random Walk (MRW), Metropolis Adjusted
Langevin Algorithm (MALA) and Metropolis Hamiltonian Monte Carlo (MHMC).
High-accuracy samplers for constrained Euclidean sampling-ias -ic.. from densities supported on convex sets K. C have been developed using Hit-and-Run and Ball Walk under various con-
ditions [Lovász] 1999 Kannan et al., 2006 1997]: J; see Kook and Zhang [2025] Section 1.3] for survey. Kook et al. [2022) introduced Constrained Riemannian HMC (CRHMC) with Implicit
Midpoint integrator and proved high-accuracy guarantees. Noble et al.] (2023) proposed Barrier
HMC (BHMC) and its discretizations with asymptotic guarantees, Kook al. (2024) developed the
"In-and-Out" algorithm for uniform sampling on convex bodies. Kook and Vempala [0024] achieved state-of-the-arr accuracy for log-concave sampling via proximal method, [Srinivasan al. (2024a] b)
showed that Metropolized Mirror and preconditioned Langevin samplers also achieve high-accuracy under suitable assumptions.
Given the above, the following natural question arises:
Can one develop high-accuracy algorithms for sampling Riemannian manifolds?

---

To the best of our knowledge, no prior work exists on providing an affirmative answer to this question. In this work, we develop the Riemannian Proximal Sampler which generalizes the Euclidean Proximal Sampler from Lee et al. [021]. In contrast to the Fuclidean case, the algorithm is based on the availability of two oracles: the Manifold Brownian Increment (MBI) oracle and the Riemannian Heat Kernel (RHK) oracle. We show in Theorem ( and Theorem g that the algorithm achieves high-accuracy guarantees under functional inequality assumptions when exact oracles are available, and under Assumption [ when inexact oracles are available, respectively. We further develop practical implementations of the aforementioned oracles that satisfy the conditions in Assumption || (Section 5). and that are connected spaces (Appendix A). A comparison with the existing results for the case of sampling on spheres is provided in Table We also demonstrate the numerical performance of the algorithms via simulations in Appendix B
proximal point method on Wasserstein
2 Preliminaries
Throughout the paper, unless otherwise specified, we use O to suppress dependency on other parameters except for and only keep leading factor. For example, log( { (log log () = Ölog ! We first recall certain preliminaries on Riemannian manifolds; additional preliminaries are provided in Appendix C
We refer the readers to Lee 2018 for more details.
Let M be a Riemannian manifold of dimension d equipped with metric g. The manifold M is assumed to be complete, connected Riemannian manifold without boundary. For point v E M, T. denotes the tangent space at z. For any v, w E T,M, we can write the metric as 92 w) = (u, w) For E M and E T.M M, exp,(t) (v) denotes the exponential map. We use grad and dV, y to represent the Riemannian gradient and the Riemannian volume form respectively.
For æ E M, Cut (ar) denotes the cut locus of x. For x. y E M, we use d(x, y) to denote the geodesic distance between and y. Let div denotes the Riemannian divergence, and Laplace-Beltrani operator A (*(()
(M) -> (M) is defined as the Riemannian divergence of Riemannian gradient: Au = div (grad u). We use vtt, x;, y) denote the density of manifold Brownian motion with time t, starting at x,, evaluated at
Let F) be measurable space. Note that the Riemannian volume form dV, is measure. probability measure and its corresponding probability density function are related through dp = pdV g. Given measurable set A E F, (A) denotes the probability assigned to the set A by P. We have (A) = JA PV=) = SA dp(x).
Definition 1 (TV distance). Let P1-P2 be probability measures defined on the measurable space (M.F). The total variation distance between P1 and p2 is defined as l| P2llTV SUPAEF Ipt (A) -P2(A)]
Definition 2 (KL divergence and x² divergence). Let Pl, P2 be probability measures the measur- able space (M,F), with full support. The Kullback-Leibler (KL) divergence and x2 divergence of p with respect to P2 are defined (respectively)
H (pz( = [ log
dp dpı,
201( [ dpl 1)2
dp2
dp
where de1 is the Radon-Nikodym derivative. dpz
is known that H (p1) 2 0 with equality if and only if p1 = p2- Although the KL divergence is not symmetric, it serve as "distance" function between two probability measures. For instance, the well known Pinsker inequality states that ||p2 pilliv- < (p))
Definition (Log-Sobolev Incquality (LSD). A probability measure P2 satisfies Log-Soboler Inequality with parameter
0 (a-LSI) if (p1) I (p1),, where (pr(Pr)
JA || grad log e1 ||= is the relative Fisher information.
p1
We also recall the definition of Poincaré inequality which is a generalization of LSL. Definition 4 (Poincaré Inequality (PI). A probability measure p satisfies Poincaré Inequality with parameter > 0 (a-Pl) if (g²) - Epigl < _ grad 9 Vg E 8 (M)
3

---

For more technical details on LSI and Pl, see Appendix H.2 and H3 In Euclidean space, conditions
like LSI and PI can be viewed as relaxation of strong convexity assumption on f. and is used to establish convergence of sampling algorithms in KI. divergence. See, for example, Vempala and Wibisono [219] (for the Langevin Monte Carlo Algorithm) and [Chen et al. [202] (for the
Euclidean proximal sampler). For a Riemannian manifold, the Bakry-Emery condition can be used to establish LSI. Informally speaking, when the potential f satisfies certain convexity, the corresponding
probability measure satisfies LSI. For more details see for example Bakry et al. [2014] and (Li and Erdogdu 2023 Appendix BJ. When the manifold is compact, it is well known that the only convex function is the constant function, and therefore the Bakry-Emery condition does not yield useful information; but for non-compact manifolds, such condition may serve as useful tool to establish LSI. Moreover, recent works translated LSI/PI conditions on I to the Dolyak--ojasiewwiz (PL) condition on f. For example, considering e -fa) in the low-temperature regime, i.e., + 0 limit, Chewi and Stromme [0224] related LSI constant and the PL constant. Similarly [ong et al.) [2024] related the PI constant and a local PL constant. Chen and Sridharan (2024) considered an "optimizability" condition and analyzed LSI/PI constant for (informally) t I O(1/d). It interesting
future work to establish similar relationships in the manifold setting.
2.1 Curvature
We also need notions of curvature on manifolds to present our main results. Let (M) denote the set of all smooth vector fields M. Define map called Riemann curvature endomorphism by R: : (MM) X (MM) x {(M) + x(M) by R(X, Y)Z = Vx V - V x Vixxy Z. While such definition is very abstract, we provide an intuitive explanation of what curvature is. Intuitively, on a manifold of positive curvature (say. a 2-dimensional sphere), geodesics tend to "contract". More precisely, given 2, y E M and e E M,M, we can parallel transport to
PV E T,M. It 22 K)d(a, y) for some well-known result that (ignore higher order terms) d(xxp, Lu, exp,, lu) < (1 2 K (which is actually the sectional curvature). From this, see that for positive curvature, which means K > O, the distance between geodesics would decrease.
Formally, given v T.MM being linearly independent, the sectional curvature of the plane spanned by v and w can be computed through K (v, w) = [=] see Lee [2018) Proposition 8.29]. On the other hand, Ricci curvature can be viewed as the average of sectional curvatures. The Ricci curvature at x E M along direction v is denoted as (v), which is equal the sum of the sectional curvatures of the 2-planes spanned by (v, 6,) where v, b2, by is an orthonormal basis Proposition 8.32].
(R(e,w))m..))
for T:M; see Lee 2018
We remark that the Ricci curvature is actually symmetric 2-tensor field defined as the trace of the curvature endomorphism on its first and last indices [Lee] 2018). which sometimes is written as Ric v) for u, E M. The previous notation is shorthand of Ric (v) Ric (v, v). When
we say Ricei curvature is lower bounded by K, we mean Ric (v, v) 2 K, Ve E T,M, |oll = 1. We end this subsection through some concrete examples.
1. The hypersphere Sd has constant sectional curvature equal to 1, and constant Ricci curvature Ric (d 1)g, Vx E M (so that Ric = d for all unit tangent vector E T.M
2. For Pm C m the manifold of positive definite matrices (with affine-invariant metric), its sectional curvatures are in the interval [- 2. 0]: see, for example, Criscitiello and Boumal 2023). Hence its Ricci curvature is lower bounded by m(m+1)-
2.2 Brownian motion on manifolds
Now we briefly discuss Brownian motion Riemannian manifold. Recall that in Euclidean space, Brownian motion is described by the Wiener process. Given x E I2d and l. >>>>>>>> 0, the Brownian motion starting at & with time / has (a Gaussian) density function v(l, y) = (2m))4/7 Ir=wil? It solves 2 the heat equation 3े v x, v) = A xt, y) with initial condition r(0, x,,) 5y)
On a Riemannian manifold, we can describe the density of Brownian motion (heat kernel) through heat equation. Let Br,, be a random variable denoting manifold Brownian motion starting at x with time and let vtt,, x, y) be the density of Br,t- The Brownian motion density v(t, x, y) is then defined

---

Algorithm Riemaanninn proximal sampler
for k 0, 2,,,, do
Y x (., wx) which is manifold Brownian increment. Step (MBI); From 2k. sample yi on X (.)) ee-f(=))(( Vs).
Step 2 (RHK): From ys. sample -ę-1 ~ n end for
as the minimal solution of the following heat equation:
Ot vt)
A x.
with v 200. = (3)
2
More details can be found in Hsu 2002 Chapter 4]. Unlike the Euclidean case, on Riemannian manifold, the heat kernel does not have a closed-form solution in general. However, some properties of the Euclidean heat kernel is preserved on a Riemannian manifold. One such property is the
- As / O, we get lim;-o / log (t, x. v) = | | On a Riemannian manifold, we have the following result.
following: Consider M Rd we have / log (1,, x) y) log (2mt)277
Fact 5 (Varadhan's 's asymptotic relation [Hsu 2002]). For all x, E M with y & Cut (z), we have
lim log (1, [-0
d( and lim grad y log vtt, x, y) exp, (x).
2
When evaluation of the heat kernel is required for practical applications, the Varadhan asymptotics aforementioned is used DDe Bortoli et al. 2022].
Yet another numerical method for evaluating the heat kernel on manifold is truncation method; see, for example, Corstanje al. [024] Section 5.1] and De Bortoli et al. [2022]. In many cases, the heat-kernel has an infinite series expansion. For example, power series expansion of heat kernel on hypersphere is given in Zhao and Song [2018] Theorem IJ, and more examples can be found in Eltzner et al. 202 Example 1-5). Similar results are also available for more general manifolds; see, for example, Azangulov et al. [2022) for compact Lie groups and their homogeneous space, and Azangulov et al. [2024] for non-compact symmetric spaces. Hence, natural approach is to truncate this infinite series at an appropriate level. For example, on S C RR, the heat kernel and its truncation up to the 1-th term (denoted as 1/) can be written respectively
31 2i + 1
Si 2i + 1
v(t, v, 3)
((2,))))) and (t,
Σ 1=0
P2((2,
4T
4x
where P. are Legendre polynomials.
3 The Riemannian proximal sampler
We now describe the Riemannian Proximal Sampler, introduced in Algorithm T Similar to the Fuclidean proximal sampler [Lee al. 2021 the algorithm has two steps. The first step is sampling from the Manifold Brownian Increment (MBI) oracle. The second step is called the Riemannian Heat- Kernel (RHK) Oracle. Recall that r11), x, y) denotes the density of manifold Brownian motion with time y). Define joint distribution T,,(,, y) xx <-=()) x. 1). Then, step 2 consists of sampling from the aforementioned distribution. When there is no ambiguity, we omit the step size r and simply write T (xx, y) x e (())((n x. y). Algorithm 1 is an idealized algorithm, in the sense that we assume exact access MBI and RHK oracles. Following Chen et al. (2022). next we provide an intuitive explanation for the algorithm from a diffusion process perspective.
Step 1:
For fixed x, we see that Ty (.,ar) v(n, -) which is the density of Brownian motion starting from æ for time y7. From this we see that the first step of the algorithm is running forward manifold heat flow: dZ, = dB,
Y X
Step We will illustrate that the second step of the algorithm running the time-reversed process of the forward process. Consider stochastic process Z, t 2 0. When we have observations
5

---

of 2n ~ Zn, we can compute the conditional probability of Zo conditioned on endpoint Zn We denote pl (T)) as the posterior. Bayes Theorem says {()))=,)) x p(zro)L(x,,2))))
where p(zn)
is the prior guess and the likelihood I. depends on the model. We consider the following model
(forward heat flow): dZz = dB, with Zo ~ π & e ee(((),, Then p( = π (xo) and (o) =
(17, 2o. Zny). Thus we get (z,)) x e -((o) Fn, Fny). and we observe that p ( is exactly x> x (xo (yy). For the forward heat flow dit = dB, with initialization Zo- T x e fe=((r), there is well-defined time reversed process 2, which satisfies (Zo, 2yy) _ (2 ) See Appendix D.2 for more details. Based on this, for the time-reversed process 2, the law of Z,, conditioned on Z is the same the posterior p=(=(2) discussed previously, i.e., π =(x)) x e e-f((=))vv v H, 2). Thus we see that the RHK oracle is, from a diffusion perspective, running the time-reversed process.
Implementing Step 1 and Step 2 is non-trivial on Riemannian manifolds. In Section and Appendix [A] respectively, we discuss two approaches based on heat-kernol truncation and Varadhan's 's asymptotics. Furthermore, geodesic random walk [Mangoubi and Smith 2018 Schwarz al. [2023] is a popular approach to simulate Manifold Brownian Increments (see Appendix B.I). however to the best of our knowledge (in various metrics of interest) is known only under strong assumptions [Cheng et al.
2022 Mangoubi and Smith 2018
4 High-accuracy convergence rates
In this section, provide the convergence rates for the Riemannian Proximal Sampler (Algorithm D assuming that the target density satisfies the LSI assumption. Firstly, note that in (Lee et al. 2021 the analysis of Euclidean Proximal Sampler is done assuming the potential function is strongly convex. However, it is known that on a compact manifold, if a function geodesically convex, then it has to be constant. Hence assuming the potential f being geodesically convex is not much meaningful. Recently, Srveettly,hhnngeg et al. [022]] discussed an analog of log-concave distribution on manifolds. Although their setting works for compact manifolds, it requires the Riemannian Hessian of the potential be lower bounded by some curvature-related value, which is still restrictive. Hence, we adopt the setting as in Chen et al. [022]. assuming that the target distribution satisfies the LSI.
we consider the case where both steps of Algorithm ! are implemented exactly, and consider the when MBI and RHK oracles inexact. Regarding notation, we let x (x:) Pi (v) denote the law of x and generated by Algorithm 1 at k-th iteration, assuming
In Section
in Section 4.2
exact MBI and exact RHK oracles. When the oracles are inexact, we let P (x),, PX (y) to denote the law of and y generated by Algorithm I at k-th iteration.
4.1 Rates with exact oracles
Our first result is follows, with the proof provided in Appendix D
Theorem 6. Let M be a Riemannian manifold without boundary, i.e., 8M 0. Denote the distribution for the k-th iteration of Algoritim ! as Ig ~ X Let k. denote the lower bound of Ricci curvature. For any initial distribution xx we have
1. Assume Tx satisfies a-LSI.
For non-negative curvature we have H₂x (px) S H₂ (xx /11 + 20 An >
For negative curvature, we have H (px) < H (p ) ( VO < n $ 1/|x].
2. Assume x satisfies a--P..
For non non-negative curvature we have x2x 1 X2x (px)) na) vn > 0
For negative curvature, we have X2x (px < X2x )/1 [11++11 ()))+)))² YO < n - 1/|m]
Note that the resulting contraction rate depends on the curvature. the curvature is non-negative, then we can recover the rate in Fuclidean space. But in the case of negative curvature, the rate becomes more complicated, and order to get the contraction rate as in Fuclidean space, we need the step size to be bounded above by some -urvature-dependent constant.

---

The above result provides high- accuracy guarantee for the Riemannian Proximal Sampler in KL- divergence and x divergence. To see that, consider for example the case when the Ricci curvature is
H (mg)
non-negative. Note that to achieve & accuracy in KL divergence, we need (I+na)ak E. Taking log
on both sides, get k = O( log(H x (or /+)) ). For small step size 17, we have log(1-na) 0 1 )- (1+
no
Hence l: = O L log H (xx)) O 1 log t) As l7 does not depend on E, that need 1pas
Öllo !) number of iterations.
There are several challenges in obtaining the aforementioned result for the Riemannian Proximal Sampler. In Fuclidean space, when a probability distribution x satisfies a-LSI, its propagation along heat flow x (0, tTa) satisfies og-LSI, with 02 I+at" This fact is very important and leveraged in Chen et al. [2022] for proving their convergence rates. A quantitative generalization of such a fact for Riemannian manifolds is not immediate and we establish the required results in Appendix H.2 following Collet and Malrieu 20088, under the required Ricci curvature assumptions.
4.2 Rates with inexact oracles
Recall that Algorithm 1 is an idealized algorithm, where we assumed the availability of the MBI and RHK oracles. Note that given I E M, exact MBI oracle generate samples ~ Th VIX (-1x). And given y E M, exact RHK generate samples
Fn X|X (-ly). In practice, exactly implementing these oracles could be computationally expensive or even impossible. For the Euclidean case, we emphasize that, as the heat kernel has an explicit closed form density (which is the Gaussian), prior works, for example, FFan et al. (2023), only consider inexact Restricted Gaussian Oracles and control the propagated error along iterations.
In this section, we derive rates of convergence in the setting where both the MBI and RHK oracles implemented inexactly. Specifically, assume we are able approximately implement the MBI oracle by generating y ~ at) |X (-Iz), and approximately implement the RKH oracle by generating in X (-ly), see Assumption I below.
Assumption 1. Denote the output of exact RHK oracle as Ti
X(Y (:lv) and inexact RHK oracle _ x ((2x) and inexact MBI as in. X1 (-ly). Similarly, denote the output of exact MBI oracle as TT.
YIX (-1z). Let GRHK and Smei be the desired accuracy. We assume that, for inverse oracle as ān
step size n-' = O(log }), the RHK and MBI oracle implementations can achieve respectively | XIY (-lv) Ty X x (-1v)llvv S SRHK, Vy, and || Y (-1z) - xn X1X (-Ix)llvv S <.MBB, Vx. We then
let Ç = max(maaaaaa (MBI).
The need for assuming the step size satisfies 1y 1/1 Ölog )) for the approximation quality is as follows. Recall from the discussion below Theorem 6 that the complexity of Riemannian Proximal Sampler depends on the step size as ({ Thus if 17 became too small, for example 1y = 0 (t).
then the overall complexity would be Poly (?), which is not high-accuracy guarantee.
We also briefly explain the intuition in assuming total variation distance error bound in oracle quality, and postpone the detailed discussion to Section 5 To guarantee high quality oracle, we need high quality approximation of heat kernel. As mentioned previously, a popular method is through truncation of infinite series. Theoretically, the Le truncation error can be bounded for compact manifold [Azangulov et al. 2002)). which says that the difference between the heat kernel and the approximation of heat kernel are close. This naturally imply an error bound in total variation distance. which motivates us to consider the propagated error in total variation distance.
We first start with result quantifying the error propagated along iterations, under the availability of inexact oracles. The proof of the following result is provided in Appendix
Lemma 7. Les pxz denote the law of X through exact oracle implementation of Algorithm ! and P denote the law of x through inexact oracle implementation of Algorithm 7 Under Assumption 7 we
have 11px (x:) 2% (x)llrv & k(GRHK + (MBI).
Based on this result, we next obtain the following result analogues to Theorem 6 the proof provided in Appendix E
7

---

Theorem 8. Similar to Theorem 6 let M be Riemannian manifold without boundary. Assume Assumption holds. For any initial distribution pú to reach (e)) total variation distance with oracle accuracy Ç {RHK [MBI log² and step size 7 = Olog ;) (for negative curvature, we additionally require S 1/
1. if n satisfies L-LSI, we need k O(log* !) iterations.
2. if satisfies a-Pl, we need k = Ö log _ log
xx((00) e²
(log2 !) iterations.
Remarks By Villani 2008, Thm. 6.15], Wi (i(,,)) v) S Dil vIITv, where D is the manifold's diameter. For compact manifolds, when constant, our TV bound directly implies a Wi bound.
Implementation of inexact oracles via heat kernel trucation
Theorem g shows that long we have sufficient accuracy of MBI and RHK oracles satisfying
Assumption I we can have a high-accuracy Riemannian sampling algorithm. In this section, we
introduce approximate implementation, based on heat kernel truncation (as introduced in Section 2 and rejection sampling. Numerical simulations for this approach are provided in Appendix B.2
First note that for rejection sampling method (in general) there are two key ingredients: proposal distribution and an acceptance rate. Assume we want to generate samples from through rejection
sampling. We choose suitable proposal distribution denoted μl, and suitable scaling constant K
such that the acceptance rate K p(a)
su(r) < 1, Vx. We generate random proposal x ~ 1l and u. E [0, 1]
being a uniform random number. Then we compute K p(2)
k(x)) and accept x if u < p(x) p(=r)*
We also introduce the following definition of Riemannian Gaussian distribution, as defined next, which will be used the proposal distribution in rejection sampling. A Riemannian Gaussian distribution centered at x" with variable t is p(t, u(fxx",x) .x) x pu(t, x², := exp
d(xr "x)²" where
28
Nh denote an unnormalized version of μu and denotes the geodesic distance. proposal distribution to implement rejection sampling, exact sampling from such distribution is well-studied for certain specific manifolds; see, for example, (Said et al.] [207]] for symmetric spaces and CCharraborty and Vemuri [2019] for Stiefel manifolds. Furthermore, this notion of Riemannian Gaussian distribution is also used in the study of differential privacy on Riemannian manifolds due to their practical feasibility [Reimherr et al., 2021 Jiang et al. 2023 In section L2 we provide an explicit algorithm for sampling from the Riemannian Gaussian distribution on the sphere via rejection sampling.
use this as our
5.1 Implementation of RHK
We first recall the rejection sampling implementation of Restricted Gaussian Oracle (RGO) in the Euclidean setting. Note that, we have log Mu (7, x, yx) 27 || yis 12 where v/, exp(- 1
2n
yxll?) is an unnormalized heat kernel (or the Gaussian density) in Euclidean space. Then we have X Tn (..k) e (f ===l=ll? Then, the RGO is implemented through rejection sampling. Specifically, we can first find the minimizer x= E arg min,, ((x + 2n ||= 11.l Note that the
minimizer represents the mode of an X|Y (.)))). We then sample Gaussian proposal Tp Nx tTa for suitable centered at the mode xt and perform rejection sampling. For more details, see, for example, Chewi 20233
On Riemannian manifold with и denoting the heat kernel, to sample from an (c, yk)) xx v-f(#))(n x. yv) through rejection sampling, we need evaluations of f(x) log v(n, x. y1). But in general, we cannot evaluate the heat kernel exactly, hence we seek for certain heat kernel approxima- tions. Hence, we use the truncated heat kernel in to replace v, and perform rejection sampling, see Algorithm 2 In the rejection sampling algorithm, as mentioned previously, we use a Riemannian Gaussian distribution as the proposal for rejection sampling. When the minimizer of g is available, we can set x² to be the minimizer; otherwise, we can simply set x² yk. We choose suitable step size
X|V
8

---

Algorithm 2 RHK through rejection sampling
Set x² yk and denote g(z) f(x) log (1 x, ys).
Set suitable / and constant CRHK s.t. VRHK(() = expl g(1)) CRHH)) $1,42 E M
exp( a )2)
for 0, 2,....... do
Generate proposal x ~ u(), Generate uniformly on [0, 1].
Return if u < VRHK (1)
end for
Algorithm 3 MBI through rejection sampling
Set suitable Ê and CMBI so that (yM == exp(log (n...y) log in (n.2..) +Cmm)
$ 1, E M
exp( Mara)
for = 0, 1, 2,.. do
Generate proposal y ~ put, x..
Generate uniformly on [0, 1].
Return y if u S (v)
end for
and that depends 7) s.t. g(x) - g(x") CRHK 2 d(x,,=")? Such inequality can guarantee that the acceptance rate (with Riemannian Gaussian distribution j(t, x", 2) as proposal) would not that the output of rejection sampling would follow #n xXXI Vv)) x (x) - log (1 2;, yk)). Similarly, to implement the MBI oracle, we also use rejection sampling to get high-accuracy approximation. Specifically, Algorithm generates inexact Brownian motion starting from x with time 7).
exceed one, i.e., VRHK (ar) 1, Vx. Then
5.2 Verification of Assumption]
I
We now show that Assumption iis satisfied for the aforementioned inexact implementation of the Riemannian Proximal Sampler. To do so, we specifically consider the case when the manifold M is compact and is a homogeneous space. Recall that denote the truncated heat kernel with truncation level I. Roughly speaking, a homogeneous space is a manifold that has certain symmetry. including Stiefel manifold, Grassmann manifold, hypersphere, and manifold of positive definite matrices.
Proposition 9. Let M be a compact manifold. Assume further that M is a homogeneous space. With truncation implementation of inexact oracles, in order for Assumption 7 to be satisfied with we need truncation level to be of order polylog(1/e) < I.
Sketch of proof: We briefly mention the idea of proof. Azangulov et al. 2022 Proposition 21] provided an L2 bound on the truncation error, and by Jensen's 's inequality we get an L1 bound as desired. With truncation level | to be of order Poly (log ) we can achieve Sar |2 xx y) vn(n, x, I)l = Ö()) See Proposition
19 and Proposition 22 for a complete proof.
Remark. Proposition 9 concerns Algorithms and which is way to implement the RHK and MBI oracles. Rejection sampling is a part of the implementation in Algorithm 2 and 3 which are both based on truncated heat kernels. The cost of rejection sampling comes in the number of steps of the for loop both the algorithms. Intuitively, one can expect that we have a highly-accurate evaluation of the heat kernel, the cost of rejection sampling should be the same that of rejection sampling with the exact heat kernel.
Remark. On the Euclidean space, for which the exact heat kernel is known, the cost of rejection sampling can be proved be ((1))Chen al. [2022]. Hypothetically, even if we have the exact heat kernel on Riemannian manifold, the cost for rejection sampling is actually unknown for general Riemannian manifolds. For the case sphere, we provide an end-to-end result (including the cost of rejection sampling) in Corollary I0 In proving this result, we first showed that when the acceptance rate in rejection sampling would possibly exceed some unimportant regions, Assumption [ still holds, via explicit computations (see Appendix [.1) Then, we show that the cost of rejection sampling (even with the inexact heat-kernel based on truncation level as stated in Proposition 9. (()) similar to the Euclidean case.

---

When M is not homogeneous space, to the best of our knowledge, it is unknown how to implement the truncation method. Exploring this direction further extend the above result an interesting direction for future work.
5.3 A concrete example on hyperspheres
We provide more specific computational complexity result that consider the dimension dependency as well as cost for rejection sampling; the proof is provided in Appendix []
S and let π satisfies a-LSI with the potential function being L1 LipLipschitz on M. Assume without the loss of generality that L1 2 vd. Consider heat
Corollary 10. Let
function f additionally
kernel truncation implementation fi.e., Algorithm 2 and 3). without minimization (i.e., start rejection sampling from yk directly), and with step size
LŽdlog 1 and truncation level
(d Poly (log 2)). To get e-accurate sample in TV distance, the iteration complexity, is k
O Lid log }), where use Ô to keep only the leading factors.
6 Additional results
In Setion A we design another practical implementations of the oracles based on Varadhan's 's asymptotics. While showing that this implementation satisfies Assumption i is left as future work, we show their connection to entropy-regularized proximal point methods on Wasserstein spaces (see Theorem 12p.
We evaluate the empirical performance of both implementations through simulation studies presented in Appendix
Concluding remarks
We introduced the Riemannian Proximal Sampler for sampling from densities Riemannian manifolds. By leveraging the Manifold Brownian Increments (MBI) and the Riemannian Heat-kernll (RHK) oracles, we established high-accuracy sampling guarantees. demonstrating logarithmic dependence on the inverse accuracy parameter (i.e., polylog( [/))) in the Kullback-Leible divergence (for exact oracles) and total variation metric (for inexact oracles). Additionally, we proposed practical implementations of these oracles using heat-kernel truncation and Varadhan's 's asymptotics, providing a connection between our sampling method and the Riemannian Proximal Point Method.
Future works include: (i) characterizing the precise dependency on other problem parameters apart from é, (ii) improving oracle approximations for enhanced computational efficiency and (iii) extending these techniques to broader classes of manifolds (and other metric-measure spaces).
Acknowledgments and Disclosure of Funding
Krishnakumar Balasubramanian was supported in part by NSF grant DMS-2413426. Shiqian Ma was supported in part by ONR grant N00014-24- 1-2705, NSF grants CCF-2311275 and ECCS 2326591.
10

---

References
Jason M Altschuler and Sinho Chewi. Faster high-accuracy log-concave sampling via algorithmic warm starts. Journal of the ACM, 71(3):1-55, 2024.
Christophe Andrieu, Anthony Lee, Sam Power, and Andi Q Wang. Explicit convergence bounds for Metropolis Markov chains: Isoperimetry, spectral gaps and profiles. The Annals of Applied Probability, 34(4):4022 (012-44022-071,, 2024.
Alexis Arnaddon, Alessandro Barp, and So Takao. Irreversible Langevin MCMC lie groups. In Geometric Science of Information: 4th International Conference, GSI 2019, Toulouse, France, August 27-29, 2019, Proceedings 4, pages 171-179. Springer, 2019.
Iskander Azangulov, Andrei Smolensky, Alexander Terenin, and Viacheslav Borovitskiy, Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces I: the compact case. arXiv e-prints, pages arXiv 2208, 2022.
Iskander Azangulov, Andrei Smolensky, Alexander Terenin, and Viacheslav Borovitskiy. Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces II: non-compact symmetric spaces, Journal of Machine Learning Research, (5(281):1-514 2024.
Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and geometry of Markov diffusion operators, volume 103. Springer, 2014.
Karthik Bharath, Alexander Lewis, Akash Sharma, and Michael V Tretyakov. Sampling and estima- tion on manifolds using the Langevin diffusion. arXiv preprint arXiv:2312.14882 2023.
Clément Bonet, Paul Berg, Nicolas Courty, François Septier, Lucas Drumetz, and Minh Tan Pham. Spherical Sliced-Wassersee Wasserstein. In The Eleventh International Conference Learning Representa- fions, 2023.
Simon Byrne and Mark Girolami. Geodesic Monte Carlo on embedded manifolds. Scandinavian Journal of Statistics, 40(4):825-845, 2013.
Djalil Chafai. Entropies, convexity, and functional inequalities: On phi-entropies and phi-sobolev inequalities. Journal of Mathematics of Kyoto University, 44(2):325-363 2004.
Rudrasis Chakraborty and Baba C Vemuri. Statistics on the Stiefel manifold: Theory and applications.
The Annals of Statistics, 47, 2019.
August Y Chen and Karthik Sridharan. Optimization, isoperimetric inequalities, and sampling via Зyаpunov potentials. arXiv preprini urXiv:2410.02973 2024.
Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono. Improved analysis for a proximal algorithm for sampling. Conference on Learning Theory, pages 2984-3014. PMLR, 2022.
Yuansi Chen and Khashayar Gatmiry. When does Metropolized Hamiltonian Monte Carlo provably outperform Metropolis-adjuster Langevin algorithm? arXiv preprint arXiv:2304.04724 2023.
Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. Journal of Machine Learning Research, 21 2020.
Xiang Cheng, Jingzhao Zhang, and Suvrit Sra. Efficient sampling on Riemannian manifolds via Langevin MCMC. Advances in Neural Information Processing Systems, 35:5995-6006, 2022.
Sinho Chewi. Log-concave sampling. Book draft available at https: // chwwisihho. tit.:eiiiiiioooo.o...ou.. github. do 2023.
Sinho Chewi and Austin Stromme. The ballistic limit of the log-sobolev constant equals the Polyak- lojasiewicz constant. arXiv preprint 1rXiv:2411.114155 2024.
Sinho Chewi, Chen Lu, Kwangjun Ahn, Xiang Cheng, Thibaut Le Gouic, and Philippe Rigollet. Optimal dimension dependence of the Metropolis-adjustee Langevin algorithm. In Conference on Learning Theory, pages 1260-1300. PMLR, 2021.

---

Jean-François Collet and Florent Malrieu. Logarithmic Sobolev inequalities for inhomogeneous Markov semigroups. ESAIM: Probability and Statistics, 12:492-504, 2008.
Mare Corstanje, Frank van der Meulen, Moritz Schauer, and Stefan Sommer. Simulating conditioned diffusions manifolds. arXiv preprint arXiv:2403.0542 2024.
Christopher Criscitiello and Nicolas Boumal. An accelerated first-order method for non-convex optimization on manifolds. Foundations of Computational Mathematics, 23(4):1433-1509 2023.
Valentin De Bortoli, Emile Mathieu, Michael Hutchinson, James Thornton, Yee Whye Teh, and Arnaud Doucet. Riemannian score-based generative modelling. Advances in Neural Information Processing Systems, 35:2406-2422, 2022.
Paromita Dubey and Hans-Georg Müller. Fréchet analysis of variance for random objects. Biometrika, 106(4))803-821, 2019.
Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: Metropolis- Hastings algorithms are fast. Journal of Machine Learning Research, 20(183):1-45 2019.
Benjamin Eltzner, Pernille Hansen, Stephan Huckemann, and Stefan Sommer. Diffusion means in geometric spaces. arXiv preprint arXiv:2105.12061 2021.
Jiaojiao Fan, Bo Yuan, and Yongxin Chen. Improved dimension dependence of proximal algorithm for sampling. In The Thirty Sixth Annual Conference on Learning Theory, pages 1473-1521. PMLR, 2023.
Maurice Fréchet. Les éléments aléatoires de nature quelconque dans un espace distancié. In Annales de I' l'institut Henri Poincaré, volume 10, pages 215-310, 1948.
Khashayar Gatmiry and Santosh S Vempala. Convergence of the Riemannian Langevin Algorithm. arXiv preprint 2r0:v:2004. 1022004////4 10818, 2022.
Mark Girolami and Ben Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society Series Statistical Methodology. 1332)1123-214 2011.
Yun Gong, Niao He, and Zebang Shen. Poincare inequality for local log- log-polyak-V) \1 ojasiewicz measures: con-asymptotic analysis in low-temperature regime. arXiv preprint arXiv:2501.00429, 2024.
Navin Goyal and Abhishek Sheuty. Sampling and optimization on convex sets in Riemannian manifolds of non-negative curvature. In Conference on Learning Theory, pages 1519-1561.
PMLR, 2019.
Ye He, Alireza Mousavi-Hosseini, Krishnakumar Balasubramanian, and Murat A Erdogdu. A Separation in Heavy-Taile Tailed Sampling: Gaussian vs. Stable Oracles for Proximal Samplers. arXiv preprint arXiv:24055. 16736, 2024.
Elton P Hsu. Logarithmic Sobolev inequalities on path spaces over Riemannian manifolds, Commu- nications mathematical physics, 189(1):9-16, 1997.
Elton Hsu. Stochastic analysis on manifolds. Number in Graduate Studies in Mathematics.
American Mathematical Soc., 2002.
Chin-We Huang, Milad Aghajohari, Joey Bose, Prakash Panangaden, and Aaron C Courville. Riemannian diffusion models. Advances in Neural Information Processing Systems, 35:2750-2761, 2022.
Yangdi Jiang, Xiaotian Chang. Yi Liu, Lei Ding, Linglong Kong, and Bei Jiang, Gaussian differential privacy on Riemannian manifolds. Advances in Neural Information Processing Systems, 36: 14665-14684, 2023.
Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the Fokker-Planck equation. SIAM journal mathematical analysis, 29(1):1117 1998.
12

---

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(1):1-50, 1997.
Ravi Kannan, László Lovász, and Ravi Montenegro. Blocking conductance and mixing in random walks. Combinatorics, Probability and Computing, 15(4):541-5700 2006.
Lingkai Kong and Molei Tao. Convergence of kinetic Langevin Monte Carlo on lie groups. arXiv preprint arXiv:2403. 2024.
Yunbum Kook and Santosh S Vempala. Sampling and integration of logconcave functions by algorithmic diffusion. arXiv preprint 1rXiv:2411.134m0 2024.
Yunbum Kook and Matthew S Zhang. Rényi-infinity constrained sampling with d membership queries. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithmss (SODA), pages 5278-5306. SIAM, 2025.
Yunbum Kook, Yin- Tat Lee, Ruoqi Shen, and Santosh Vempala. Sampling with Riemannian Hamil- tonian Monte Carlo in a constrained space. Advances in Neural Information Processing Systems, 35:31684-31666, 2022.
Yunbum Kook, Santosh S Vempala, and Matthew Zhang. In-and-Out: Algorithmic Diffusion for Sampling Convex Bodies. arXiv preprint arXiv:2405.01425, 2024.
John M Lee. Introduction to Riemannian manifolds, volume 2. Springer, 2018.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Logsmooth gradient concentration and tighter runtimes for Metropolized Hamiltonian Monte Carlo. In Conference on learning theory, pages 2565-2597. PMLR, 2020.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Structured logconcave sampling with a restricted Gaussian oracle. In Conference on Learning Theory, pages 2993-3050. PMLR, 2021.
Benedict Leimkuhler and Charles Matthews. Efficient molecular dynamics using geodesic integration and solvent-solute splitting. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(1189)20160133 2016.
Mufan Li and Murat A Erdogdu. Riemannian Langevin algorithm for solving semidefinite programs. Bernoulli, 29(4):3093-311 13, 2023.
Han Cheng Lie, Daniel Rudolf, Björ Sprungk, and Timothy Sullivan. Dimeeniion-independent Markov chain Monte Carlo on the sphere. Scandinavian Journal of Statistics, 50(4):1818-1888 2023.
Chang Liu and Jun Zhu. Riemannian Stein variational gradient descent for Bayesian inference. In Proceedings of the AAAI Conference Artificial Intelligence, volume 32, 2018.
Chang Liu, Jun Zhu, and Yang Song. Stochastic gradient geodesic MCMC methods. Advances in neural information processing systems, 2016.
László Lovász. Hit-and-run mixes fast. Mathematical programming, 86:443-411, 1999.
Oren Mangoubi and Aaron Smith. Rapid mixing of geodesic walks manifolds with positive curvature. The Annals of Applied Probability, 28(4):2501-2543, 2018.
Michelle Muniz, Matthias Ehrhardt, Michael Günther, and Renate Winkler. Higher strong order methods for linear Itô SDEs on matrix Lie groups. BIT Numerical Mathematics, 6294):1095-1119, 2022.
Tomohiro Nishiyama and Igal Sason. On relations between the relative entropy and x 2-divergence, generalizations and applications. Entropy, 22(5):563, 2020.
Maxence Noble, Valentin De Bortoli, and Alain Durmus. Unbiased constrained sampling with self-concordant barrier Hamiltonian Monte Carlo. Advances Neural Information Processing Systems, 36:32672-327111 2023.
13

---

Adam Nowak. Personal Communication, 2025.
Adam Nowak, Peter Sjögren, and Tomasz Z Szarek. Sharp estimates of the spherical heat kernel. Journal de Mathématiques Pures el Appliquées, 129:23-33, 2019.
Sam Patterson and Yee Whye Teh. Stochastic gradient Riemannian Langevin dynamics on the probability simplex. Advances in neural information processing systems, 26, 2013.
Gabriel Peyré. Entropic approximation of Wasserstein gradient flows. SIAM Journal on Imaging Sciences, 8(4):2323-2351, 2015.
Mare J Piggott and Victor Solo. Geometric Muler-Maruyama Schemes for Stochastic Differential Equations SO (n) and SE (n). SIAM Journal on Numerical Analysis, 54(4):2490-2516 2016.
Matthew Reimherr, Karthik Bharath, and Carlos Soto. Differential privacy over Riemannian manifolds. Advances in Neural Information Processing Systems, 34: 12292-12300,, 2021.
Salem Said, Hatem Hajri, Lionel Bombrun, and Baba C Vemuri. Gaussian distributions Rieman- nian symmetric spaces: statistical learning with structured covariance matrices. IEEE Transactions Information Theory, (4(2):752-772 2017.
Simon Schwarz, Michael Herrmann, Anja Sturm, and Max Wardetzky, Efficient random walks Riemannian manifolds. Foundations of Computational Mathematics, pages 1-17, 2023.
Vishwak Srinivasan, Andre Wibisono, and Ashia Wilson. Fast sampling from constrained spaces using the Metropolli-adjusted Mirror Langevin algorithm. In The Thirty Seventh Annual Conference on Learning Theory, pages 4593-4635. PMLR, 2024a.
Vishwak Srinivasan, Andre Wibisono, and Ashia Wilson. High-accuracy sampling from con- strained spaces with the Metropolis-adjusted Preconditioned Langevin Algorithm. arXiv preprint arXiv:2412.18701, 2024b.
Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. Advances in neural information processing systems, 32, 2019.
Cédric Villani. Topics in optimal transportation, volume 58. American Mathematical Soc., 2003.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2008
Andre Wibisono. Sampling optimization in the space of measures: The Langevin dynamics composite optimization problem. In Conference on Learning Theory, pages 2093-3027. PMLR, 2018.
Keru Wu, Scott Schmidler, and Yuansi Chen. Minimax mixing time of the Metropolis-adjusted Langevin algorithm for log-concave sampling. Journal of Machine Learning Research, 23(270): 1-63, 2022.
Xiangjin Xu. Heat kernel gaussian bounds on manifolds i: manifolds with non-negative ricci curvature. arXiv preprint arXiv:1912.12755 2019.
Tianmin Yu, Shixin Zheng, Jianfeng Lu, Govind Menon, and Xiangxiong Zhang. Riemannian Langevin Monte Carlo schemes for sampling PSD matrices with fixed rank. arXiv preprint arXiv:2309. 04072, 2023.
Chenchao Zhao and Jun S Song. Exact heat kernel hypersphere and its applications in kernel SVM. Frontiers in Applied Mathematics and Statistics, 4:1, 2018.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Yes, the end of Section for summary of our main contributions, well as references on the theorem number.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: See section for discussion on future works.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

We did state the assumptions in each theorem, and the proof can be found in appendix.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification:
The information needed to reproduce the toy experiments are provided.
Guidelines:
The answer NA means that the paper does include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
16

---

Answer: [Yes]
Justification:
Codes were provided in supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/
sttpiiisssooooiocyyy. for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The details that are necessary to understand the results are provided in this paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [NA]
Justification: Our work is theoretical and includes some toy examples. The plots are average over 1000 number of trails.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
17

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer:
[Yes]
Justification: The toy examples are run on personal laptop using Matlab, and only CPU (AMD Ryzen 7 PRO 5850U).
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics htttp::/neuurpp..... http:///wwuurp..ccccaaaaaaaaaaaa
Answer:
[Yes]
Justification:
The authors have read the NeurIPS Code of Ethics and followed it the paper.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
Our work theoretical.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
18

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work is theoretical and poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: The paper does not use existing assets.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
19

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: Our work is theoretical.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details of the dataset/code/mmodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Our work is theoretical and does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our work is theoretical, and does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
20

---

16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy buctps://neure wwwparippps..coooeeeeeeeee for what should or should be described.
21

---

A Implementation via Varadhan's 's asymptotics and connection to entropy-regularrized JKO scheme
In this section, we consider yet another approximation scheme for implementing Algorithm [ motivated by its connection with the proximal point method in optimization, where the latter is in the sense of optimization over Wasserstein space (Jrdann et al. 1998. Wibisono 2018 Chen et al.
2022]. Note that the proximal point method is usually called as the JKO scheme after the authors of Jordan et al. [1998].
Specifically, we consider approximating the heat kernel through Varadhan's asymptoties. Let 0(1,2,8) Xy exp(- " =: Pw, x, y) be an inexact evaluation of heat kernel. According to Varadhan's 's asymptotics, lim,-- (1,,, x) v(n, 2, y). Hence when 1j is small, & is good approximation of the heat kernel. Note that (1), 2, .) in Varadhan's 's asymptotic is exactly the Rieman- nian Gaussian distribution μ(n, 2).. Denote (z. 8) exp(-f(x) d(a.v)3) ). With inexact MBI 27
implemented through Riemannian Gaussian distribution and inexact RHK implemented through rejec- tion sampling (Algorithm) 2 to generate x-XIX (xrly) x exp(-f(x))
d(z.9) 2q ), we obtain Algorithm
4
For the case when M Sd, we prove in Appendix [ that sample from x (xly) through rejection sampling, with suitable parameters, the cost is ((1) in both dimension d and step size similar results for more general manifolds seems non-trivial. Numerical simulations for this approach are provided in Appendix B.3 Verifying Assumption [ for this implementation is open.
1) Obtaining
Algorithm 4 Inexact manifold proximal sampler with Varadhan's 's asymptotics
for 0, 2, do From Zs. sample y/ ~ # xXIX (-, mι) which is a Riemannian Gaussian distribution.
From yk. sample Tb+1 ~ xX|Y (., mk4 x e- -f())) d(x.zf) end for
using Algorithm
A.1 RHK as a proximal operator on Wasserstein space
We first show that the inexact RHK output in Algorithm 4 can be viewed as a proximal operator
on Wasserstein space, generalizing the Euclidean result in Chen et al. (2020) to the Riemannian
setting. Recall that with function f and d being a distance function, prox nf (v) arg min, f(x)
L d( 2) The (approximated) joint distribution is f v) exp(-f(x) dlz.z)², By direct 27 computation we have the following Lemma (proved in Appendix G..
2n
Lemma 11. We have that
X
arg min Hęx (p) + (p. prox nHax (8y); pEPa(Mf)
27
which shows that the ineact RHK implementation is proximal operator, ie., #x\Y=y prox x (6
A.2 Connection to entropy-regularized JKO scheme
Observe that in Algorithm 4 the Riemannian Gaussian involves distance square, which naturally relates to Wasserstein distance. Now, recall that for function F in the Wasserstein space, its Wasser- stein gradient flow can be approximated through the following discrete time JKO scheme [Jordan]
al. 1998
1
Pk+1 arg min F(p)
W3
(EP(R")
2n
If' M is a smooth compact Riemannian manifold then the Wasserstein space P2(M) is the space of Borel probability measures on M. equipped with the Wasserstein metric Wz. We refer the reader to Villani 2003 for background on Wasserstein spaces.
22

---

It was proved that as 1) + O. the discrete time sequence {pk} converge to the Wasserstein gradient flow of F. Later, Peyré 2015] proposed an approximation scheme through entropic smoothing of Wasserstein distance:
arg min F(p) + W2 (p, Pa); PEP(R*)
Pk+1
21
where W2, is the entropy-regullarized 2- Wasserstein distance defined by (here H is the negative entropy)
/ (),,,))27((; tT)
W (P1,
inf rEC(OL.P))
In Euclidean space, Chen 2022 showed that the proximal sampler can be viewed an entropy- regularized JKO scheme. We extend such an interpretation to Riemannian manifolds. Specifically, we show that Algorithm [ which is an approximation of the exact proximal sampler (Algorithm IT.. can be viewed as an entropy-regularyy that on a Riemannian manifold the negative entropy is (7) == JM ? log())dV g())
JKO as stated in Theorem 12 (proved in Appendix G Note
Theorem 12. Recall that x X xe Let zks yk I be generated by Algorithm PKtz be the distribution of Tk. Uk: Th+1, respectively. Then
Let pk Pk and
/ fdx + 2n
P = arg min
W2 (P and Me+1 = arg min XEP2(M)
W 22 (P).
XMP 27
B Simulation results
B.1 Brownian motion approximation via geodesic random walk
In our experiments, to compare against the Riemannian Langevin Monte Carlo Algorithm, we used the geodesic random walk algorithm to simulate the MBI oracle following Cheng al. [2022). De Bortoli et al. 2022], Schwarz et al. [203]; see Algorithm More efficient implementation is topic of great interest in the literature; for example, [Schwarz et al. 2023
Algorithm 5 Approximation of manifold Brownian motion using geodesic random walk
Input x: E M,, > O.
Sample & being a Euclidean Brownian increment with time in the tangent space M,
Output = (C).
While it well-known that geodesic random walks converge asymptotically the Brownian motion on the manifold, non-asymptotic rates of convergence in various metrics of interest is largely unknown. basic non-asymptotic error bound for geodesic random walk is available Wasserstein distance (see CCheng et al. [0022] Lemma 7J). Mixing time results are provided in Mangoubi and Smith [2018). However, such result is not immediately applicable to establish high-accuracy guarantees for the Riemannian proximal sampler, when the MBI oracle is implemented via geodesic random walk. An important and interesting future work is establishing rates convergence for geodesic random walk in various metrics of interest so that those results could be leveraged to obtain high-accuracy guarantees for the Riemannian proximal sampler.
B.2 Numerical experiments for Algorithms and 3: von Mises-Fisher distribution on hyperspheres
In this experiment, we test the performance of Algorithms 2 and 3 for sampling from the von Mises- Fisher distribution on hyperspheres and compare it with the Riemannian LMC method. In this case, we have f(r = -kp²² z. Note that this f(x) has a unique minimizer on Sd. This implies that LSI is satisfied, see [Li and Erdogdul 2023. Theorem 3.4]. We demonstrate the performance of our Algorithm on S C R with μ (10, 0.1, and = 10, and on 5 with μ (5, and K = 10. For the purpose of numerical demonstration, we sample the Riemannian Gaussian distribution through rejection sampling.
23

---

LACT CUB Aoaittin SaTOer
spp..
(a) S f(x) = exp(K)
(b) f(x) exp(xpp x)
(c) Ps. f(X) expl
2
Figure 1: Frechét variance (i.e., E (( x=" )²1 versus number of iterations. Left and Middle figure correspond to the implementation via Algorithm 2 and Right figure corresponds to implementation via Algorithm
To evaluate the performance, we estimate E[d(x, {. where x² is the minimizer of f. representing the mode of the distribution, and plot it as function of iterations. Note that the quantity Eld(x, x²)x is referred to as Fréchet variance
19988. Dubey and Müller (209)]. For this, generate 1000
samples (by generating samples independently via different runs) and compute Toyo E x(x;,==)²²
We use rejection sampling to generate unbiased samples and get an estimation of the true value. Due the biased nature of the Riemannian LMC method, achieve high accuracy we need small step size. Contrary to the Riemanian LMC method, the proximal sampler is unbiased. and it can achieve an accuracy while using large step size and smaller number of iterations; see Figures (a) and (b). For both algorithms, we use uniform distribution on the hypersphere as initialization.
B.3 Numerical experiments for Algorithm 4 manifold of positive definite matrices
In this subsection we illustrate the performance of Algorithm for sampling on the manifold of positive definite matrices. Let P... (X E GT(m) X X and y Xy > Vy E [R'" be the set of symmetric positive definite matrices. According to Bharath et al. [0033 Section 6.2), we can choose g(U, V) tr ----- and make (P'm.) .) Riemannian manifold. it non-compact manifold with non non-positive sectional curvature, geodesically complete and is a homogeneous space of general linear group (m). Additional details are provided in Appendix C.3
We test the performance of Algorithm 4 when the potential function f(X) 2n Im), m 3, a 0.03, following Bharath et al. [2023] Note that f is not gradient Lipschitz. In the Figure we estimate E [d((2, x²" [ and plot it as function of iterations, where x" = L3 is the minimizer of f. representing the mode of the distribution. For baseline comparison, we run Riemannian Langevin Monte Carlo for 200 iterations with decreasing step size to get a reference value of E[d(x, x )²) which serves as the true E,d(x, xx" )). Similar to the previous experiment, we generate 1000 samples from independent run, and compute 1000 Σ 10000 d(x,, )² for each method. For the Riemannian I Langevin Monte Carlo method, we find that if we set step size to 0.001 instead of 0.0001, after few iterations the algorithm diverges (potentially due to lack of gradient Lipschitz condition). But for the proximal sampler (which is unbiased algorithm), even with large step size as illustrated in the plots, the approximation scheme still works well and can achieve a higher accuracy than the Riemannian LMC algorithm. For both algorithms we initialize the algorithm at Xo = 2/3. Note that if we use a random initialization, the Riemannian LMC algorithm (prior work) might diverge (potentially due to lack of gradient Lipschitz condition in our potential).
B.4 Numerical experiments for Algorithm cost of rejection sampling for higer dimensional case
We also demonstrate that the cost for rejection sampling is not exploded by dimension. Following the same setup as Appendix B.2 we consider sampling from a von Mises-Fisher distribution on $100 with μı = (10, 0.1, 1, I,... 1) and x = 10. We test the performance of Algorithm 4 in terms of rejection sampling cost. We choose 7 0.0001, and compute the average number of iterations executed by rejection sampling. With a total number of rejection sampling oracle being x 10 the
average rejection sampling cost is found to be 2.661.
24

---

C Additional preliminaries
C.1 Divergence
We will briefly discuss divergence for the manifold setting, More details can be found in Lee 2018 Recall that in Euclidean space, for vector field F ...,... Fr. in R" divergence of F is defined as V Di=1 0f: has natural generalization the manifold setting using interior multiplication 0z: and exterior derivative.
The Riemannian divergence defined as the function such that d(ix (dV)) ) = (div X X)V,, where X is any smooth vector field on M, i denotes interior multiplication and d denotes exterior derivative. See for example Lee [0188) Appendix B] for more details. On a Riemannian manifold, recall the volumn form is dV, (det(gi))dz A... dx Let Y
Ei= Y a We can compute the
Ow:"
interior multiplication as
(dV,)) = Vdct(grs)
Σ ((-1) dx))))dd A... A die A... A...A dr
j=1
=ett((gs))
> ((-1))+tyy)d A... A dil A... dr
=
j=1
We can then compute its exterior derivative as
1- 0(Y' V /de((gs)) dx )xx A dād dx"
((i)))V)))
dx,
==1
0(Y' /del(gi)) Ox; O(Y' /det(gis)) dVg.
/ /det(gs)dd A A dx
det(gi))
dzj
i=1
0( /det(n)))
Hence we get div (Y)
In Fuclidean space, this reduces to
(dct(gg))
S2,
div (Y) = Eド= For E (M) and X E X(M), the divergence operator satisfies the following product rule div (uX) = u div (X) + (grad u, X).
ays Ozy
Furthermore, we have the "integration by parts" formula (with ğ denote the induced Riemannian metric on iM)
J (grad u, X) = f( u(X, u(X,NNydV - -f u div XdV.-
When M does not have a boundary, OM = 0. So we have
S,A. (grad u, X) - f div XdV.
C.2 Normal coordinates
Riemannian normal coordinates. Let M. There exist a neighborhood V the origin in T, M and neighborhood U of æ in M such that the exponential map exp. V U is diffcomorphism. æ. Given an orthonormal basis (z,) of T,M, there
The set U called a normal neighborhood basis isomorphism from T.M to Rd. The exponential map can combined with the basis isomorphism to get smooth coordinate map &p U + Rd Such coordinates are called normal coordinates at æ. Under normal coordinates, coordinates of x is TRd. For more details see for example Lee 20188 Chapter 5]
25

---

Cut locus and injectivity radius. Consider E T.M and let be the maximal geodesic starting at with initial velocity v. Denote (z
sup(t > 0 the restriction of Yoo to [0, 1] is minimizing) The cut point of a along To is Yollout (I, w)) pro- vided cut (z. v)
{q E M q is the cut point of z along some geodesic.). The injectivity radius at x is the distance from to its cut locus the cut locus is nonempty, and infinite otherwise [Lee] 2018 Proposition 10.36].
00. The cut locus of I is denoted as Cut (x)
When M is compact, the injectivity radius is positive [Lee] 2018 Lemma 6.16].
Theorem 13. (Lee 2018 Theorem 10.34] Let M be a complete, connected Riemannian manifold and x E M. Then
1. The cut locus of z is closed subset of M of measure zero.
2. The restriction of exp= to ID (I) is surjective.
3. The restriction of exp, to ID (x)) is a diffeomorphism onto M\ Cut (ar).
Here ID (x) {u T.M: |u| < leut (x; [e)) }} is the injectivity domain of x.
Then for any E M, under normal coordinates, for all well-behaved f, we have
J MM fdv,
=J -Sn,
.f( | Cut (p))CRd
fdV,
(00 (g)d.
JM,.Cuu Cut (p)
C.3 Additional details for manifold of positive definite matrices
We briefly mention some properties of Pm. The inverse of the exponential map is globally defined and the cut locus of every point is empty. For symmetric matrix S E mmmxm.
expx (1S) x1/2 Exp (IX - 11////XX/xxxxx111 1/2)x1/2
x Exp (t Log (X, 1/2 X2XT11/3 xx//2² is geodesic that connect X1, X2, (t)
m
d()
(log(r:))² with ri being eigenvalues of X, X2.
expx (X2) = (00) = x1 Log (X1 1/2 X²x, 1/2
1x2xx111xx0 1/1/2
We have the following fact. Lemma 14. Let (6xx)= = d(x, u) with y M being fixed. We have grad Φ(x) -2exp (y).
D Proof of main Theorems
For given o define the d-divergence to (p) E [o))(. Define the following dissipation functional
(x) =- Ep (grad (Ø' f). log 6
π
We can now compute the time derivative of the di-iivergence along certain flow.
Let μu be the law of the continuous-time Langevin diffusion with target distribution x fxx), That is, have the following SDE, dX,, grad f(X,) v2dB, Then, nx satisfies the following Fokker-Planck equation (see Lemma 25 for a proof).
0 u div (Nr grad f(X)) + grad ) div (grad X grad π ur x-
div (μ. grad log 1xx xxx)
We now show that D2 (u+K) -0,4,x (4x
26

---

Lemma 15. We have that
1x xx ex).grm grad log 12 3, xx (1/)
Dax (xx = px (grad (ø'
xx
Proof. [Proof of Lemma
15
By using the fact that & u = div (ut² grad log #x), we have
0
= 0
I
Öt xxx (1x²)
0t JM
x* x0 x
(x)
= S/x o
σt a (V)(())= Ju (8' + xxx (
a
Mr div (ul grad log ) dV
me (grad ø' 1x x
Mr x (x))
grad log
M
where in the last equality we used integration by parts.
To get more intuition on the notion of d-divergence and dissipation functional, consider ((x)
r log(z). We get KL. divergence and fisher information:
Th (p) F log( f] π = [lgg L) = Hx(p);
grad (p) = EEl (ø' o L), grad log - = Epll grad log L = 18x grad L L = (p)
π
π
π
Our proof is now based on generalizing the proof in Chen et al.
2022 to the Riemannian setting. In Section I.TT we analyze the first step of proximal sampler by viewing it as simultaneous (forward) heat flow. In Section D.2 we analyze the second step of proximal sampler by viewing it as simultaneous backward flow. Combining the two steps together, we prove convergence of proximal sampler under LSI in Section D.3
D.1 Forward step: simultaneous heat flow
We can first compute the time derivative of the <--iivergence along simultaneous heat flow. Lemma 16. Define Q describe the forward heat flow. Lel px Q2 and x Qe evolve according the simultaneous heat flow, satisfying
1 0,p Q 2/2 A(0 Qi) and Oxp 1 A(Txx(
Then 0,0 x (pxQ0) = - Dr (p* (Q))
Proof. [Proof of Lemma
16
Denote p:² == p= =x and क. xxxQr Then, we have
② (p,
/
xxx 0l pt (dV)(z)
0t n,
Øt M
= 2 / o
+ 0'(
=
pr
x
(
-
p dV,
M
Ot
at
x
0t
Recall that by construction,
Dspp" 7 Ap
1
1
div (grad x div (pf grad log )
2
and O,T Am
x
1 div (Tº grad log n.² Hence, we get
2
H x (xx)) -2 ( Sl T
i $
8' 15
2
0
x
Pt dV,
(x))
n
01
at
òl
div ( grad log m.
M
27

---

Pz
div (pt grad log pi
(div (nz grad log n, ))
d2)(2))
Tp
-1 grad of ) ot grad log m.
M
Ti
X
Pi
P. )) o. grad log nH dV)
(grad d'( ) Pi grad log pt¹) ) grad 0
Pi
xxx At
-( grad di
क. grad log mt²)
grad ø'
Pf pt grad log oggpe
Ti
grad
Pi
a
grad log x
+
Ph
grad @' Pt
()
xx
) o.
grad log o, dV
Now, notice that
grad of Pi
X
) ae
x
grad log m. = (d'(
px
grad
Pi uxx
क.. grad log o..
x
"t
So get
1 H
" (pt)) = / < Pi
grad o ) grad log=y" ) - grad 8'( ) re grad log log (dV)(
M ae
T
- - pi¹ grad d'( grad log
(2)
-Bpx grad (#) mt
grad log -Dxx (1:))
D.2 Backward step:
simultaneous backward flow
We leverage the following result.
Theorem 17 (Theorem 3.1 in De Bortoli al. (2022). For SDE dX, b(X.)dt dBr, let pe denote the distribution of X,- Denote Y. = XT-1, E [0,T) to be the time-reversed diffusion. We have that dY, (-b(Y)) + grad log (og.rr-(YY)) + dB,.
Note that the time reversal can be understood as (Yr, Yo) has the same distribution as (Xo, Xx).
Recall that vtt, x, 3) the density of manifold Brownian motion starting from x with time and evaluated at y, and that 7 y) = xx (I)) y). We denote x² >xxQ Qn to be the Y -marginal of π (I, 2). Let me = x Q.. Consider the forward process d.,, = dB, with Xo ~ =x We know that the time-reversed process satisfies dY, = grad log T
dB,.
Define Q as follows. Given pk, set 0 Qi to be the law at time t, of the solution of the time- reversed SDE (with T = 17). Thus Yo p", we get XT ~ p¹. By Bayes theorem Xo ~ JM XIY (xly)dp² (y). hence Yr ~ Sar XIY (xly)dp" (v). For the channel Q we have
1. Q7 is the identity channel.
2. Given input pK, the output at time n is p² 0 = SMr x XIX (xle)ddp (1)
3. x Qi = x x -
Thus see that the RHK step of proximal sampler can be viewed going along the time reversed process. We now have the following result.
Lemma 18. For the time-reversed process, we have
0,0 0.0+r0 (p'Q7)) Q7) -1)==11-p--00 Dzvo (p' Qi).
o.
28

---

Proof. [Proof of Lemma 18 Denote T, = =T' Qr and p. = =p"
1 Qr. =p"Qu- The Fokker- Planckk equation is
1
0,m, div (म, grad T, + Am A
2
2
1
(r)
Oipi div (pi grad log क. ) Apī div (pi grad log Api
2
t 2
Hence
20,4 n, (pi =2
Φ(
Pi i
+
d'(-=)
01 Pt
0
Te
0
Pt Tt
dV, dV (z)
=
AAT
M
Pi
div (pr grad log
Pe
Pe ATT ) [VV Apr
Tt
Pr
Pr )dV((1)
=2
div (Pr grad log
M
T1
+ J/
Pt A 60( -0'
)App Api + 00
Pt Pt
Ar, dV (z)
Tt
Te
2 8' ((_)
div (pr grad log Pi [VV
Dx; (pi
M
1
2D, x, (pt + D DI (pi = -D, n. (p: )
where we used the same steps as as in the proof Lemma [6 obtain
/
8' Pt ApT
+ 0²( Pt
Pc
A., dVy(2)) Dm (pT);
M
Tt.
Tt
and used integration by parts, to obtain
ef 8'( Pt M
div (PT grad log JdV, = -2 S
Pt
(grad
0'( o'(2²), ). grad log Pe JdV (x)
Tt
nt
nt
2E [[grad 6'( Pc 6'(2t
l = -2D x;- (pr)
grad log
T
Kt
D.3 Convergence under LSI
Now prove the main theorem.
[Proof of Theorem 6 iiemmm] [] We first prove the theorem assuming curvature non-negative. For the general case, we only need to replace the LSI constant a1, 0,,0
Proof.
1. The forward step. We know xx satisfies LSI with ag H Using Lemma 16 we have
1
O.H Q1) = - (pt Qi) Jr
- (pt Qc).
This implies H x, (0) Q:) Ae H₂x (vr)) where, A. a,,s log(1 + lxx. We also have A, (1+ ta)-1111 As result,
1xxq, (P0 Qnn) < H (pr))
+ 1ja
29

---

H (pt)) 2. The backward step. Using Lemma 18
we have
0,H Q, (pă Q,) = - (P -0; H Qi (pó Qr).
Since n ='Q = ==xx, Qu- to we know the 1.SI constant for xx QT oj == (m-1)+
Same
as in step I, we get A = So at dt = log(1 + an). As a result,
Q (P0 (1)))) - H (PO) H Y
1 ta
3. Putting together. We have x = P/ =p² 1n Qn Denote pp Pb Qūj we get
Hg (pt) H 0r (0) Q S H (pt) H Hxxq,, (p Qn) (1 ta)
-
t1+
(1 + ta)²
4. Negative curvature. For negative curvature, use a, as in Proposition a == do). We compute the integral
31 (the value to be
x integrated is I-e knge where 1
a(c 1) + 16
(1: 2e- + dx log(a(e" - 1) + K) log(k) = log(
Hence we have HG (pk)) S H (pO ) Ma(enn-)rr 1)+x ) Observe that in general, for [0, 1] we have that
e Thus for nj [R]
we have |m|n < 1, hence 2 2 e -1log This implies -ee-iein hand, we have 1 - S e which implies 1-e- + 1-e-lkin @(1 an). [x]
[r] 2 2 On the other
S 17. So we have a(e"" 11
[x]
To summarize, have that when O,
K
1
a 1) K
1 + 2
18
which implies Hx (ot) S Hax (pt) )( 11
26
D.4 Convergence under Poincaré inequality
In this subsection, we extend Theorem item from LSI to PI, under x² divergence. We follow exactly the same strategy, but we first discuss how to modify the proof from the LSI setting to P setting. Recall that æ satisfies Poincaré inequality with parameter a, we have
1 Ell grad dı
du
Var
(12)
a
Note that x² divergence is d-divergence with Φ(z) = (x 1)2, which allows us have (ar E, dx ) E, [du x (n). Also, recall that we defined in Appendix D
the dissipation of d-divergence. By definition, the dissipation of x² divergence is D= HEpKgrad (#t 1). grad log MM 2E, ll grad du de1]?]. Hence we can interpret the Poincaré in- equality as
x2 < D=(p)
2a
From this perspective, to follow the proof for LSI setting, we need to know whether the Poincaré inequality constant is also preserved along heat propagation, in the same way as LSI constant. The
answer is yes, see Appendix H.3 Proposition 32
With this, can prove Theorem
6 item
following exactly the same procedure Theorem item
Proof. [Proof of Theorem item We first assume curvature is non-negative, then discuss the situation that curvature is negative.
30

---

1. The forward step. We know π rxQu satisfies PI with == it.. Using Lemma 16 and
Proposition 32 we have
Qu ( Q1) = - - Qe) S -uxxax (0b Qc). (1)
This implies x xxx (po Q1) S Se-A. H (pt)) where, A: So a,ds log(1 ta). We
(1)
also have A, = (1+ta) -1 As a result,
X (p (xx
X Qn (p 2n) -
1 17
2. The backward step. Using Lemma 18 and Proposition 32 we have
0rx Q, (N QT) - - xx (00 Q aa X 07 (.)
Since 7 xx=== = xxx0 = we know the I.I constant for x NT is o7 == (1-1)+
Same
as in step I, we get A = -0" a dt = log(1 + an) As result,
X2x(p)) 1 ta
x Q, (Pd (0)))) -
3. Putting together. We have nx n Qu v5 =nē =x²Qn, Qt Denote n³ = pă QTj we get
Xax (p%) x (pă (p) Qn)
X² (pb) (1 ta)
X2x (pt) 2 e (g)
(1 + ta:)
(1+
4. Negative curvature. For negative curvature, use ar in Proposition integrated is -- where _ == do). We compute the integral
31 (the value be
R
dx log(a(e"t - 1) + x) log(x) log( L (C2 -2- + + t)
a 1) + K
k
Hence we have ZS (px) S Hxx x (po a((enn-1)- 24
Observe that in general, for x E [0. 1] we have that 1
2 2 c= Thus for n .
we have |x/n 1, hence lel 7 2 1elo This implies 1-e-
[n] 2 17 On the other
hand, we have 1 - x < - which implies 1-e-lein I So we have a1e"n-1)+x + -he O(1 + an).
x To summarize, we have that when O,
K
1
a<exnn 1) +
1 + a I-
1 2n
which implies x3 (px) X²x x (px)) (1+ 22
E
Proof of Theorem g
Recall that pi (x), PK (1) denote the distribution generated by Algorithm I assuming exact Brownian motion and exact RHK. This notation is applied for all k. For practical implementation, using inexact RHK and inexact Brownian motion through all the iterations, we denote the corresponding distribution by 2 (x), PX (y) respectively.
1, we are at distribution (x). Denote (x) to be the distribution (z). obtained from /xx (x) using exact Brownian motion. (Note that (v) denote the distribution obtained from x (z) using inexact Brownian motion).
Note that at iteration
31

---

We now prove Lemma] 7
Proof. [Proof of Lemma]
Using triangle inequality, have
I (z) Ph (x)]l - / pk (v) XIV
( (xlu)dy / (u)* ( (alw)doll
<1/ (1 x (xlp))) I + #11 [1 -Z XIV (x
The first part can be bounded by CRHK:
/ (2) (xly) * (xlw)) TV - sf f | - x (alv)
dy
TV
SSRHK.
For the second part, we have
/ (ot (pk-1(1) pt_ (v)) ((())))) TV
1
2 ///
(Dk /1/0/11000
1 -o% (v)) (al)))))(( dx 1/2 po-)(1) P0-0001 / (x1x(x)))
2
=11p% 1 (2)1()) S 11% Pk liv + 111k-1()) 1
<(BBI + 1 (x) - (x)llv)
Here, the last inequality follows from Lemma 37 Together, we have
l (z) Px (x) S CRHK + <MBI + 1 pk (x)]ITV.
Iteratively applying this inequality and noting that 11 (x) - pă (x)llTV = O. we obtain II (x) 2 (x)llrv S k((SHH + (MB)).
Recall that Pinsker's inequality states || u/ITV (2HL(...
Proof. [Proof of Theorem item
For simplicity, we assume non-negative curvature. The negative curvature case follows from the same proof strategy. Using Pinsker's inequality, we have
l1 n 1TV S
v 1xx s
V 2 + n
H Hxx (pt))
S 5c.
We want to bound l| - x 11v S fec. It suffices to have
H (oo (1+najak
S 10a
Hence we need
log " (xd") 2k log(1 na), i.e., k log(1+7a))
log( 211 x (pp) )
For small step size 77, we have Tog[1+na] 0 Hence k: = O
log M Mxx(ox)
-6(H log 4
na
Recall that by assumption,
Ō(log 2) We pick Ç
log" and consequently ,
Ölog log' } = Ölog _ + 2 log log ) = Ö(log 4). follows that k Ō log 1 1 = Ōllg )
The result then follows from triangle inequality:
1 T lrv S 11 pk lrv + lllt - xx lrr S k(CRHK + (MBI) + _ = Ö(s)
where kç- = =Ō I log² 2) = Ō(E).
32

---

Proof. [Proof of Theorem 8 item 2 For simplicity, we assume non-negative curvature. The negative curvature case follows from the same proof strategy, We know that for two probability measures T. P we have log(1 + (p)) 2 (0) 2| - p where the first inequality is Nishiyama and Sason 2020, equation 16] and the second inequality is Pinsker's inequality.
Hence, we have
X² (p V2 log 11
l x 11V S
V V2
Hxx (px) 1 log(1 + X2x (pk -
(1 77 )2k
We want to bound l - T* || S fec It suffices to have log (1 xax (px) 1 (exp( - 1)(1 ) that is,
X=x (p0)
S 1 Jea
We need
(1+72)
log Xxx (pt)
Xxx (pg)
k 2 2 2og 1 + +
(exp(}e2)-1 - = 00 log(11 + t72
log
E²
Hence O log x xxx(
(pt))
For small step size 77, we have log(1+na/2) 0 1 naa
Recall
that by assumption, = (log t). We pick = log2! and consequently 1 = Ö(log log 1 Öllog ! + 2 2og log ! = Ö(log ) It follows that
k Ô log x²xx (po))
(og 1 log xax x
(pb)
=
17
82
g²
The result then follows from triangle inequality:
li n lrv s 1 pk lrv + llek - x x llr S k(GRHK + (MBB) + 2e = Ö(s)
where kÇ ((E)
F Verification of Assumption 1
In this section, we consider implementing inexact oracles through the truncation method. Recall that we assume M a compact manifold, which a homogeneous space.
We use + X to denote the output of MBI oracle and RHK when rejection sampling is exact. More precisely, since we use the truncated series to approximate heat kernel, we have ,YIX x (T(1), x, y) and XIY x </=))) x, y). When rejection sampling is not exact, i.e., there exists E M s.t. (z) > 1, we denote the output to be # 7 for inexact Brownian motion and inexact RHK, respectively.
In subsection F.1
we prove Proposition 9 i.e. || XIX 1lv = Ō(C) and | π XIX
1rv
Ó()) with
log - I Then we know 4 and x satisfy Assumption
general setting, where the acceptance rate is allowed to exceed at some unimportant region. We show that on Sá, for certain choices of parameters, || x l| Õ(() and || - =VIX || = Ô(C). This means that allowing the acceptance rate exceed in unimportant regions would cause significant bias for rejection sampling. It then follows from triangle inequality that म x|x and म X|Y satisfy Assumption
In subsection I.I we consider
F.1
Exact rejection sampling
We prove Proposition 9 i.e., verify that Assumption satisfied with Ç log" y required in
Theorem g
33

---

F111 Analysis in total variation distance
The first step is to bound the total variation distance, under the assumption that heat kernel evaluation of high accuracy. We consider the following characterization of total variation distance (see Lemma 36
1
f (a(())>p((x)lV((CC))
|li -22llV
2
Proposition 19. Let M be compact manifold. Let be the desired accuracy. Assume for all y E M have Jas (n x. y) vi (n, x. y) v))dV( = (()) and for all E M we have JM (1, vr(n, x, y) VVV (u) = O().. Then || n x ~ = O(s) and || T
x = O()).
Proof. [Proof of Proposition 19
Step 1. Note that A₁ == supzEM e [(x, A2 == infaem e e- f(x) are positive constants independent
tf t. Denote Z = JM e (1 I, y)VV and Zn = SAr e (() x, (z), We know
IZ2 Z11 = -((x) v) c x(1(z) (1,
- A1 / [p(n,,,)) -1
Hence, we have
(17,))()111))
1 XX
x18 1rv S
(I ( (v)
x
J JM E
Zze 11) x.
Zie 1() v) (x)
MM
2, m -z} ----(()) -1(()) (1,1,, v() (z)())
27,
+ /
[Z2 Zi .maxfe -5()) x, y). ((m) (10 dv)
M
2Z
< Ō() + 01/ 12: Z max x, v) - Ö())
where by Lemma 21 we obtain minn(2,22)
Ō(1) and
27,122
max ((1,,, (00))))) vz vi} dV, dv) OL)))
22,
M
Step 2. Denote Zi = JM (1 x, v)V ( (v) be the normalizaing constant for vi. Since v the heat
kernel, we simply have JM v() x. y vVV (v) 1. It holds that
1I (1, v. v) Sa vn (V)((
and -IXX v(1. v).
Then,
1
x
VIX
<
JN MI (1, v, WdVV - v(1, x, (3)))
2 M5
1
[1(1,, F 8) Ziv(1, xr R.1 dVy(y)
2
M
min 2i, 13 |r x) v (1,,,,,)) dV, (v)
22,
M
Zrl-
max [u(1, (1, x, vn aV,(u)
+
22:
M
34

---

O()).
Theorem 20 (Theorem 5.3.4 in Hsu 2002) ), Let M be compact Riemannian manifold, There exist
positive constants C1 C2 such that for all (1,z, (0, 1) M M.
C1
# #t
C2
Su(t,z, - (1dd-1)
21
14/2
Lemma 21.
We have 1/ JM e g)) = (1).
Proof. [Proof of Lemma
Using lower bound of heat kernel from Theorem 20 we have
e (())
v(,,,,
M
/
C1
d(z, 3)
A2C1 17012
/ exp(
d(ar y)
A2
exp(
)dV
vVy(2))
M 17012
2n
21
A2C1 17//2
A2C1
1/122 C2
C4
Hence
1/ /
e --(4) -1,0),,,,, Ö())
M
F.1.2 Analysis of truncation error
Now we discuss the truncation level needed to guarantee high accuracy evaluation of heat kernel as required in Proposition 19
Proposition 22. Let M be compact manifold, and assume M is a homogeneous space. With 1 = Öloog 2) and Ç To to reach [[[[,, x. v) ve(n, x. y) Ō()) we need I Poly (log !
Consequently, achieve
f ((n,,,y) - vn = Ö(c)
and
M
/ (,(,, m(n,,,,))dv Ō(C).
we need Poly (log ))-
Proof. [Proof of Proposition 22 Following Azangulov et al. 2022
Proof of Proposition 21] we have
2/3/
[1 x. v) - 21,,, x 2,,1, S C'I
1
22/2/4 Requiring C' O c Ö(C)
Take = log J. Recall that in Theorem & we require ( To equivalent to
2//
1 -
Log
C'l
<
log
12/a
Take log on both sides, we get
log2
S log
This further implies
C'ilog"
E
1
1
2/d 2- -log
C log
(4log log
+ log
+ logl C")
C'llog
g
=
35

---

It suffices to take I = Poly (log ) We verify that / = Poly (log ) can guarantee the bound:
12/ Co 17 = Cog -e cloxt = Poly (log _ )e Poly (log = =0(
R/a
= O()).
logg
a homogeneous space, both and in are stationary AAznguulvv et al. 2022]. Hence
JM 2,(,,,,,) (n, a, y))dV, (z) does not depend y, and Jas [Mn 2, y) 1n, 2, y) (v) does not depend on æ. Therefore using Jensen's 's inequality,
J a,, x1 ( 0(1I((1,, in m1 v, v)
20(1u(,,,,) (1, 1))
Note that the same holds for JM [v(n. x. v) (n, v) v)) Hence get the desired bound, i.e.,
JM [yn, x, y) vi(n, z; (V))) = (() and JM [1(1, x. y) vn 2; 2, ()) = Ō()).
G
Proofs for entropy-regularized JKO scheme
Proof. [Proof of Lemma
Note that, we have
x x (p) = / p(x) log p
(v) / M
p(x)
p(x;) log
dV,
exp(-f(x)
d(x.y)²
exp(
dxx.))²
)
27
2n
J/x (log Cexp(-S(x)) d(2.y)2
p(x)
1
+ +og
)dV )V))(z)
d(x.y)² 2n
exp(
2n
J. p(x) log C'(x)xxXY dV
p(x)
2 vdv
p(x)
2n
1
-
2n JM d(z, dd + CCy),
where C=
J/ To
Cyyy and Cly) are some constants that only depends on y. The above
computation implies
2n [ d(x, v) C(y)
x arg min H HXXIY arg min Hą (p) NEPz(M)
PEPa(M)
arg min HĄ (p) + 2n PEP2(M)
kg, = prox 111 (6y))
Lemma 23. The minimization problem
min JCPx(M M xx J
1
d(x, y)²²1();)) ()VV, H()),
MXM 27
where the constraint means JM 1(a, vVVV = xxxx), has solution of the form (x. pK (x) Y (u/x))
Proof. [Proof of Lemma 23 Since JM (z, vVV px (z), we have
S. M
(// 7(,, ((x,))dV - /² (x)) B(x)dVV(x) 1
0,
we can construct the following Lagrangian
M ( 11, ((VV))VV(3)) -xx(x))))))))))x)(()))) (())dV())
M 2m d(x, v)2 v) (z)dVV)) + H(7) -
36

---

Recall that (x) .MXM ( We have,
(1y + t(p) H(7)
lim (1 +7) - -0gg()))
lim
1
lim
4 plog(7 + (4) + p
t
MXM
J MxM
p(og(()) + 1)dV
For any function f. denote 1,(7) M 2lz,
We then have
1 + tp) - 15(7)
Sa pt)f 2f
lim h-0
lim
S M pfdV (z)
Thus the variation of Lagrangian is given by
JAYMI x (_, w) + log(7) 1 ((x) ( (vd
We want the above to be zero for all sp Thus we need equivalent to
1 2n dxx y)² + log(n) + 1 -x)) = 0 which is
7(. y) = ((x) d(2
This implies (xx v) ((x) ) Integrating with respect to the y variable, we get
p (x(x)= ((x)=/[x000 (I, ((x) d(=,)))- -Vv(VVV)))))))) (a) f e - dV)
then follows that
e +o d(x,))
((,,1) x (x))))
p² (x)* (yylo))
SM t(.( dVg(y)
Lemma 24. The minimization problem
(M =or J. (f(x)+ + - d(x, u) (a, (z)V, H()),
JMM
2n
where the constraint means Sar 1la, vV (z) px (y), has solution of the form
(2,1) p² xIV
Proof. [Proof of Lemma The proof follows similarly to that of Lemma 23
Since
JM (1, y)VV = px (v), we have
1l., (x, J
(())) = 0, AB.
We first constructing the following Lagrangian:
J (f(x) + d(r,, (,,,,))))((x,))d)V(( + H(w) H()) S (S (L,,,,, - -2a
1
MXM
2m
Recall that H(x) JM ? log(l))())))))))) l Then, have
H(> + t+p) H())
J (7 + txp) log(7 + +4p) - (og())
lim
lim 1+
61>0
37

---

lim
plog(> + (1)) + (x)ddV
MxM
+ + (x))
MXM
For any function f, denote 15(7) MMX v)
We have p
lim 1 + (4)) - 15(7)
lim
+00
S
pfdV,
=
MxM
Thus the variation of Lagrangian is
JNXM p. (f(x)
(-(())) + d(x, v)² + log(n) -
1011-(()) B(v) aV (x)
2m
We want the above to be zero for all y. Thus we need f(x) d 1)² + log(?) + - B(v) 0 which is equivalent to
(1, v) (())))((()-- -(1
This implies y(2, v) elv) , d(2,y)² Hence integrate with respect the variable and get
(v(2) [ (((6)))-(0)) (x) ((v)) ( e -((())- d(x,v)? dVz
A
Therefore, we obtain
- [(x) d(2
(2,0)) x p²(2)) (v)
= =px (1)=xxIx (xlv).
JM e f(x) 2, d(x.)) dV)
Proof. [Proof of Theorem 12
By definition we have
1
arg min XEPa(M))
W
2.2n (Px x)
20
arg min
d(v, u)) H(())
XEPG(M))ECU x) MMM 27
By Lemma 24
we know the solution of
S MMM 2n
min
d((z, d(z, gd) H())
(MMM) x
(xx u) x ph (x)e 24 d(x)))) Hence the -marginal of inexact proximal sampler satisfies
x(v) / xxx (z)e d(,,)) dV2(x) = S,, o (I) (ylz)dV, -x)).)
M
Similarly,
xEPa(M) 27 2,2
arg min
W 22, (pk, + fdx
1
arg min
(f(x) _d(v, w) ())),( ((d +H(7),
MM .x) MxM
2n
and its solution is x(x) SM 1x (x)))dv = (x))
38

---

Auxiliary results
H.1
Diffusion process on manifold
well known that the law of the following SDE dX, = -b(,)d + dB, related to the Fokker- Planck equation 0ep: - div (pxb(X,) grad Pt). Here we provide proof for completeness.
Lemma 25. Let B, denote Brownian motion on Riemannian manifold M. For SDE dX. -(())dd + dB, the corresponding Fokker-Planck equation is
3ep: div (pxb(X,) + grad pt).
Proof. The infinitesimal generator of the SDE is Lf f.b) + 1 f Cheng et al. [2022] We compute the adjoint of T. which is defined by JM hT. SMs h..fVV, By divergence theorem, we have
2 / div (grad f)hdV, 1 2 (grad h, grad f)dV, / div (grad fdV, M
1
(grad S. bh)dVV / div (bh)fdV,
Hence
I hLfdV
1 h(gradd f,b) + -hasdV,
S.d. f div (bh) + 2 fAhdV, / MM
f(div (bh) Ah)dVg.
Thus obtained L*h div (bh) + Ah. By Kolmogorov forward equation Bakry al. 2014 Equation 1.5.2], we get
1 AAt = -div (b(X))px 2 grad p.). O.p: = Lt = div (b(X))p.) +
We briefly mention some properties of Markov semigroup. The following results are from Bakry et al. 2014 Section 1.2].
Definition 26.
1. Given a markov process, the assoicated markov semigroup (Pi) defined (for suitable f)
P.f(z) = = =
2. Let be the law of X then Pt is the law of Xx, We have
f P.fdp
[ Jas sd(P..).
3. Markov operators (Pt)22 can be represented by kernels corresponding to the transition probabilities of the associated Markov process:
P.f(x) J (()) 0, M.
Thus by definition, we have
Bff(X)) = = =P][(f( = (1),,))
(v)
39

---

H.2 Log-Sobolev inequality and heat flow
In the sampling literature, log-Sobolev inequality usually written the following form:
J log fog J&l Rdulog [.A f S
P
2
2/a,
ll grad f11 Js
1 Jolp).
Hu(p)
2a
In the Euclidean setting, know if Un #2 satisfy a1 LSI respectively, then their convolution
Ji P2 satisfies LSI with constant
see Chewi (2023 Proposition 2.3.7]. In particular,
we take one of la to be vtt, 2, y) (which is Gaussian in the Euclidean setting), since the Gaussian density satisfies LSI, we have the following result.
Fact 27. Consider Euclidean space. Let μ be a probability measure that satisfies a-LSI. Then its propogation along heat flow, denoted by pe 4 also satisfies LSI with constant I+t 14ta
Here Vt denote the probability measure corresponding to heat flow for time t.
On a Riemannian manifold, the density for Brownian motion satisfies LSI.
Theorem 28. [Hsu 1997 Theorem 3.1] Suppose M a complete, connected manifold with Ric M Z -c. Here e 2 0. Then for any smooth function on M, we have
[ M
log log fldv
grad flle... + 11fll log Il Ieo..
With = -C, we know the Brownian motion density for time satisfies LSI with constant I- RT.
As special case M = Rd, we have 0. Hence, the LSI constant became lim od L. That is, (with representing the measure for Brownian motion with time t) (p) S (p),,, Vp. So the LSI constant for Brownian motion 0g =
In the following, we prove that on a Riemannian manifold, such fact still true. We follow the idea by Collet and Malrieu 2008 Theorem 4.1]. For notations, we denote I'(f) = I(f.f) grad f fll
also require the tollowing intermediate result.
Lemma 29 (Theorem 5.5.2 in Bakry et al. 2014 For Markov triple with semigroup (Pt)t2o, the followings equivalent:
1. sc
2. log - log (Pff) S c(s) (P) where c(s) 1-e 28.
25
Corollary 30. With Pt denote manifold Brownian motion, have
1. T(P.f) - (PP/)))
2. Po(flogf) - P.flog(P.ff s (t)P((TL)) where c(t) = 1-e -ne
2x
Proof. [Proof of Corollary 30 For the second item, can replace by g2 for some g.
J 3² log. 2 - / dd log( / gddvz) <
23g
23s
(2g)² grad gll? du
23
g²
M
4
e
20s
ggaa gl du..
2,3
Now already know the manifold Brownian motion density va satisfies
I-c at -LSI, i.c.
f f2 log Pd
/ dovelog log
/ fdv, 2 1 1-e-m [
||grad flld Vf.
M
M
40

---

So we know, with P. representing manifold Brownian motion, Poglog f) Pif Poflog(P.f) 1 c()).) where
1 nt
1 e 222
2
23
Hence we know / can be taken as K, 8 corresponds to Ltt So we get SEPP)CE U
Proposition 31. Let M be a Riemannian manifold with Ricci curvature bounded below by R. Let Po be any initial distribution. Assume pu satisfies LSI with constant do'
J 9 log dpo - J, 9 pd log J. gd S 2do S |rad gll- dpo. Vg.
Then the propagation of Po along heat flow, denoted Pt. satisfies LSI with constant
ke
2c(t) + doe e + ndge-at ext - 1 + +Kdo
where c(t)
1-e -M
If 2 have c(t) S st.
1
2c(t) + doe nt 2 + dne at. t do
Proof. [Proof of Proposition 31 Since po satisfies LSI with constant +. equivalently with replace
g², we get
/ M
/
fdpo S 2do / I|rrd dd
do 2 JM f
/ grad S
flog fdpo
M
ídpo log
Jss
dpo-
For > O, using Corollary 30 we know the manifold Brownian motion (here represented by P.) satisfies
Pilg log 9) Pig log(()) < c(t) [(g)
g
where c(t)
Using property of markov semigroup, we have
2x
1/ Pilg gogg sc(t) / P. p.(
P ddpo + S P.g Pog(()))
g log gdpt
M
Hence
1 M
1 M
/ gdp:
g log loggdp:
gdpe log
[(
J Polg)dpoloy log / P(
Sc(1)
P.
ydpo +
Pig Peglog(Pig)de
JM
M
[(g)
do
{ I(Pig)
<c(t)
P.
)dpo
dpo
g
JM
Pig
l(g)
do
(P:VT(g))
Sc(()
M
Pa
9
ddo
+
2
E
-et
M
Pig
dpo
Po [(g)
do
T(g)
Sc(t)
+ )dpo
e
P
ydpa
M
2
M
g
= =(c(t)
+ do
E kt)
/
P.
[(g)
ddp
(c(t)
do e e-m))
/
AM
[(g)
dp..
2
where the third inequality is due to Corollary 30 and in the last inequality we used Cauchy-Schwarz inequality:
r(
grad fl
Igradd f|1 ./² = -Ell gradd S[]
PA PPPC
E
E 2 E
f
41

---

= (Piv (PiVT(f))²
Hence we know or satisfies I.SI with constant
2c(t) + due
1 - c + ridge
Note that we have limx-00 be(()+dae Euclidean space in the limit.
This means that, we can recover the result for
1406
H.3 Extension of Proposition 31 to Phi-Sobolev incquality
In the last subsection, we showed in Proposition 31 that Riemannian manifold, LSI is preserved along the propagation of heat flow. In this subsection, we extend the result to the setting of Phi- Sobolev inequality. We will prove the following proposition. This result will be useful for proving convergence of the Riemannian proximal sampler under x2 divergence and Poincaré inequality Appendix (.4).
(Theorem
6 ditem [ the proof
Proposition 32. Let M be Riemannian manifold with Ricci curvature bounded below by Assume 6: 1- R. is a C+ function and convex, where I is an interval of Further assume z" is concave on Let Po be any initial distribution. Assume Po satisfies o-entropy inequality [Bakrr et al. 2014 Section 7.6.1] with constant to:
S.L (g))do $1 J gdpo) S do /
(")))T( Vg.
Then the propagation of along heat flow, denoted as Pt. satisfies o- e-entropy inequality with constant
ke
2c(t) + doe nt 1 + nnde-at
ext - 1 + kno"
where c(t) I-c
2x
We first present some useful lemmas.
Lemma 33. Let & be such that (ي( is convex. Then
P((T(())'())) ("(Pig) 2 (P.(VT()))
In particular, when σ(z) (x (0-1)- 1)?, we know (2() = (P)(V/((g)))
is a constant, and therefore P; (()))))) 2
Proof. Using Jensen's 's inequality, since
is concave, we have (F
(a))
P((T(()))(O)) E(T())" (')
2E(r() M B(
(P:(VT(G))
where in the last inequality used Schwaa-Schwarz.
Lemma 34. For diffusion Markov triple with semigroup (Pt)izo, the followings equivalent:
1. The curvature dimension CD oo) holds for some B E R.
2. Se
3. Po(flogf) - P.f (P.( c(s)
4. P pcss))
42

---

where c(s) 1-c
Furthermore, the first item implies P:(o(g)) (oa(i))))oPPss))) S (()
2Ø
Proof.
See Theorem 4.7.2 and Theorem 5.5.2 in Bakry et al. 2014 and Theorem 2.1 in Chafai
2004).
Corollary 35. With Pt denote manifold Brownian motion, we have
1. (((P)) Se (
2. Po(Slog,]] - P.flog(PPI)) & (PPP.
3. P - (P.())² & (0()).((())
4. P - ((P.)) & (()R(((())((
where c(t) 1-E-
Proof. [Proof of Corollary
35
For the second item, replace by g² for some g-
Jn g² du, J. M
dv, log( / g² gddv.)
1 e 28s
(2g)²l grad gll dv.
2,3
g²
1 25a
grad dv.
2
Now already know the manifold Brownian motion density is satisfies
T-e at i.c.
J/ log - - f2 / fdu, S 2 1- nt
/ || f12 Vf.
M
So we know, with P. representing manifold Brownian motion, Pog f) Pef log(Pif) s (()P..(F))) where
4 1 e -25s
2
2,3
Hence we know 3 can be taken as corresponds to Lt So we get item 1, and 4.
Now we are ready to prove Proposition 32
Proof. [Proof of Proposition 32
For g O, using Corollary satisfies
35
we know the manifold Brownian motion (here represented by P.)
Pe(o(g)) ((P.g) s (()))
where c(t) Ie 2x
Using property of markov semigroup, we have
[ blg)dp: S Polo(g))doo <(t / Pi(o" (1) + / ((Pig)dpo.
Hence
1 olg)dpr oi / gdpi)
M
<d() / ((())())(T))dpp)
M (P:g))do (ocillooo//p 0 Prlg)dpo)
sc(t) / ((())())((o))doo
do
+
(PP):))
2 M
M
do
<c(t) P((o"())(())dd)))
+
c
e
d'(Pa)(PPTT())))
2
M
43

---

sdt) / Pi(o"()))(o))dd + do _ / (((()))"(()))d JM
2
=(c(1) + do 2
(c(1) + do xxt
2" (e)amm)amm)))
I 6" (g)
where the third inequality due to Corollary 30 and in the last inequality we used Lemma 33 Hence we know satisfies d-entropy inequality with constant
1
2c(1)= + doe
1 e- nt + kdge
We remark that the case <(x) = (a 1)² recovers Poincaré inequality, and <(Φ) = x log x recovers log-Sobolev inequality.
H.4
Total variation distance
Lemma 36. For TV distance, we have
[ dp()) dpl2)
11 - 1200
|TV = sup
(A) - (()) 2 3.M dV, dV
dV g1
ACM
and
/ sappa - ldp 10 - 22) TV 1-
sup
Denote the set at which supremum is achieved to be A, for M : (x) 2 (x)). Denote p! p(2) to be the measure, or corresponding probability density function with respect to the Riemannian volumn form, when appropriate.
Proof. [Proof of Lemma 36
(x) - (x))VV((()) = / (x) (x))dv + / (x) - (z)(dV\)))
M
(x) n (z)
(T)) (z)ldV())
= f (x) (z) (x) (3(
1
=p(² (A9 (1) (A.) p (A.2 (A,) +1
(9 (A.) (A.)) = 2) 1 I|r..
Now we prove the second equation.
1 [dpl [ sappil $ J ((() (o)))(()(()
S sup [f(x)] / (x(2 (x)(V)) ZEM
11 (x) (2)(e)))) = 21 - "p ||V.
When f = 1A. - 1As+
1 fdp / fdp(t)/ = ((2)(A.)) - (11)(AA)) - (A.) + (1))
= -20 - (A.)) = 2)) 1rv.
44

---

(2) Lemma 37. Let ,00) p22)) be probability measures. Let pi , denote propagation of p( p(2) along heat flow on M, with Po
= p p0)) - p1 We have l (1) - p,-112v I =1 lIv..
Proof. [Proof of Lemma 37 By definition we have that for all f.
S.a. (1))
E Eff(X) =i = Pif(z)
Assuming Xu p(), we get
I M
P1f(z)dp((
I I.d
S. g(z)p
Where we denote JM (( - / pr DV = 1.
g(x). Note that
|g(x) < /
a5
Hence
(2)
1 sup
f [dpp" - fdpt"1
Pi
TV
2
f:M- [-1,1
1
sup
2 J:M 1-111
S I (V
1
g(x)p (0 - g(r) (x)dV
sup
M 2 1:M-1-1,11
1
g(z)
((z) (2))))))))))11)
sup
2 9:
= =11 (2) lTV-
I Concrete example: hypersphere
L. Truncation method on hypersphere
Let M be hypersphere. Previously, we discussed some existing results which provided a bound on the La norm of # - v. For hypersphere, we can derive a bound in Leo norm, see subsection [.1.2] We also consider the situation that the acceptance rate in rejection sampling might exceed 1, and show that for such a situation, rejection sampling can still produce a high-accuracy sample.
Let ViBI(y), VRHK(2) denote the acceptance rate in rejection sampling. Recall that #11x vn(1,,25 2r, y) and X x e en xr, y). In the actual rejection sampling implementation, if for example in Brownian motion implementation, it happens that there exists E M, s.t. V(y) > 1, then the output for rejection sampling will not be equal to # x For such situations, denote V (1) = min(1, VMM)(y)) and V VRHK(2t (x) = min(1, VRHK(2)) Note that V MBI(y) and V (a) the actual acceptance rate in rejection sampling. we denote the corresponding rejection sampling output by Y XX respectively.
45

---

Intuitively, the region (r) near I carries most of the probability for both Riemannian Gaussian distribution jutt, u, y) well as Brownian motion v'tt, æ. y). when the variable is suitably small. Thus instead of choosing parameter to guarantee VRHK(*), (yM < 1, Vx, y E M, it suffices to guarantee VRHK((*) S 1, Va E By(r) and (y) S 1, Vy E (r) for some r.
Let L be the Lipschitz constant of f. In the rejection sampling algorithm, will use Lj d(z.p)² as proposal Riemannian Gaussian
(i.e.. 77 = Lidlog ) exp( d(z.z) exp( 21 (2(0-2)1or?
2(en)
distribution.
Define
(y) = (x log vi (7, + CMBI) $
Vy E Br(r).
exp( d(z.y)²
2(mm)
VRHK (r) exp(-f(xx) + log (1, 2, f()) log (n. v. + CRHK) 11, Ve E By(r).
cxp( 2c dx. u)
1
Lemma 38 (Nowak al. 2019]] Nowak 2025 }). Let d denote the dimension, and denote the heat
kernel (that corresponds to dt v 1,, v) Stil, v, u)). For ip represent the size of neighborhood. we have
d)
C(v,d)
1 exp( - (d((,,))) vd(x, v) S 4p
exp( 21 d(,, < vll,
14
t
Fix O44). When & (2),, holds that Clo.d) and c d),
(1). this case, we denote C C(p,d)
c(p,d)
Proposition 39. Let M = S be hypersphere. We set = 1,d the step size of proximal sampler, ; = (d - 2) log & as the parameter for proposal (Riemannian Gaussian) distribution of rejection sampling, and truncation level 1 = Poly (d2 Poly log }). There exists parameters
d
CMBI - log(22) log(C + 1),
CRHK 2 log ?
log (2T) log(C 1)
S.1. क X|Y #t X|X satisfy Assumption 7
Proof. See Proposition 42 and Proposition
Lemma 40. Xu 2019 Equation 17] Let (M4, .g) be complete manifold with Ricci curvature being non-negative. Then we have
exp( d(x,z)2
(t,,
(2t))
L11.1 Proof of Proposition 39
It's important to establish the order of constant in algorithm first.
Lemma 41. There exists some CMBI log(2rr) log(C' + 1) s.t.
- log vn (n 9) log (n, 2 d(I, y)
2(sT)) + CMBII
Consequently, for we have
VMBI(y) exp(log vn (n x. log (1) 2, CMBI)
< <, (r(r).
exp(
d(a..))
2(sn)
Proof. Recall that we require 1 3 (log )
Write = C, log & where Cn Ljd. Then can write पेक S.
46

---

1. Step 1 Consider neighborhood (ro) with rō = 2.
We have
C..
C
d(z, y)² ) + e क.
C
exp( - d(z,
11),, - v) +
exp(
)11 + 8(x, 1))
211
1)
2n
nen
with 6(x, == e expl-
For this ö, we can see that
29
exp( d(a
2.))
c aon
1
8(a, n) C
s c vy E BB((o))
exp(
d(u.y12) 27
C
n
Hence C(1 + 8(x, 7)) S C +++.. 1. which further implies
(n. 2 - log C' + d
log + d(x,v)² log (1 + 6( n))
- llg u llog
2
2n
log n + d(x, y)
h - l (C' + 1) + 2
2sy
Now consider a slightly larger neighborhood where s will be set later: C, S d(a,y)
5
we have
expl
1
5(x, 7) = C exp(
S 17
d(a.y)) 2n
C
C²-- -!¹
that when ( is small, we have
1
log(1 + 6(x, n7)) = log 6(2, n) S log
+ (1
log
C
1-1 3 < log
1
1 d(r, v)
d(x.y)²
2 1
= =log C +
+
V Cn
-
C
C
n
2
On
which further implies
d
d(x,
- log (1) v. v) 2 - log C +
log +
log(1 + 6(2r, 1))
2
21
d
d(x, k2x,
d(s. y)²
2 - log C +
log +
log C
2
2m
n
2
d
log 17 + d(I,y)?
2 2
2(sT))
log(( +1) 1)+ + log 11 + d(z,y)²
2(sT)
Together, we conclude that for all y E (r(r),
- log vi x,,) v) 2 log(C + 1) + d
2 log + y2
2(s))
2. Step 2
We have 2,) 2 1
and consequently
(2 of
1
(2mm) 1
(2TT)
Thus for all y E (r)r), for some constant C' we have
- log in (1), x. + (1 x, 2 log(C + 1) + d
d(x, y))² d
d
2 log
log n 2 1lg((T))
2(sm)
2
d(x, y)
d log(2m)
(2m) - log(C (C + 1)
2(s
2
47

---

Therefore there exists some CMBI = log(2x)) - log ( + 1) s.t.
VMBI(U) = exp (7, log v (11, x, x) + CMBI) 1 E Br(r).
exp(
d(w.y)²) 2(an)
Proposition 42. Let M be hypersphere S so that the truncation error bound can be proved in Lec-
Consider Algorithm 3 with == 1j8 where -ータ > is constant that does not depend on v): Ç. For || = O())
small S, the error for inexact rejection sampling with vs is order C, i.e., 1 V x
Hence by triangle inequality, || I// (c).
Proof. [Proof of Proposition 42
Now we have for all y E (r)
d(x, 1) 2(sT)
log 2, 21,
+ CMBI-
This implies,
exp(log log + CMBI) exp( d(z, v)
)
2(s77)
Hence have
(MM = exp 14 (17 x, v) (1, 2, + C7)
1, (r).).
exp(
d(a.y)2 2(s7)
Recall that y denote the density for Riemannian Gaussian distribution. We compute
VMB
VM
p(sn, 3)
E9 Epp(s VMMI(y) JM (MMI()))u(T), x, y)dV (y) (,,,,I,))
VMBI (y (v) y)
JM on M (n,
22, VV
exp(
2tam)
exp(log in (n. v) log V (1, C7)
JM exp(log (1 x, 8) log u2 x) + C7 CV)
11(11, x,, JA (n) x, vVVV
x
Thus the desired rejection sampling output can be written
(1sn, x. 9) VMBI(y)
# x|X
By(sna MM)
On the other hand we denote VMBi = min { 1, (M) and the actual rejection sampling output
YYIX = p(sn, v, 8) Ep(an.=..) Vieerty)
Viei(y)
Following Fan et al. 2023
Proof of Theorem 6], we get
VIBI
MBI
| Y X TT < E
d(c.y)2 21am)
exp(
E
[MBBB]
E [MB
2TF [MMBI
VMBIl
[[MBBi
We aim to derive an upper bound for
2E Vini
Vini E[
Note that
exp( d(a.yy?
B|| |MBBl] 2
exp(log (1 x, log (n. x,) x) + CMBI)
2(sT)
)dV,(u) dV,(v)
(r(r))
exp(
d(x.u)²
J& exp(
d(n.y)
2(an)
2(an)
48

---

/ exp(log log (1, + (V
JAs expl
d(z.0)²
vdV(()) (r) 2(
and similarly
exp(log log (n,, + CMBI)
BM[B VMBIl =Bp(sn...9)
1M
exp(
d(z.0)²
2(em)
J / exp(logg((n,2,0)) 2, - log 10,, + CMBildV, (v).
JM exp( d(z.y)² 2(mm)
JvV (V)21)
Hence we only need bound
2 VMBIl
2 (VMBI(W))> ( - log CMBi) vV,
[B[Ve
JB. exp(log (1) 2. y) log vn CMBId
(nin.ay (V)(v)
{ [(imm(u)>1) Pr y) JV,
2 (Vmm(()>11 miln,.))
[B
mi(.,.y)) dV)()) (1(,),,))
[B.(r) U/ (1 2,
J (1 C Poly (2) expl d(w.y))
)dV,
2n
< 2
[B.00) u( x. yy)
J Ki 2-1 S Poly ( exp( d(x.y)) (xp(-day)) (dVz(()
2n
2
O())
JBAr) vr(n, x, y)V (w)
We used (n, S + + un,, v, y) In the last equality, note that when d(x,y)²
2-
it holds that
On
exp( d(x, y) < exp(
( -
2-1 1 C,Cr
2n
and then we know 2² Poly (2) = O()) because for small S, the term 5² Poly (2) = O(1).
We used lower bound S (r) 2n (1 (v)V
Ja v(n, x. VVV (v) (.(r) (2mr)
exp( d(x, y)
(dVg(1)
B.
2n
and for the choice of r, we know this is lower bounded by a constant.
Proposition 43. Let M be hypersphere Sd, Consider Algorithm 2 wwith " = 12d & and =
(d 2) log 4. For small ç, the error for inexact rejection sampling with is is of order Ç, i.e. || x 1zzv Ō())
Proof. [Proof of Proposition
1. Step 1 Set d-1 in Lemma 41 Let " LRd ). C, Lfd. Note that have /2/2 2-d-1 d-1 d+1 Lid L7 and 1 [ 2) log _. Note that / = d-2") We know, for all
d
x E By(y), we have
d(z,y)² y)²
- log (n, + log un v. (..
CMBI-
2(sT)
We want to find CRuK so that the previously defined for proposal distribution, i.c., we need f(x) log v1(n, x, v) f(y) + logui(n..vv) 2
ب can be the variable
xe d({ y2² CRHK to hold for all æ E (y) Hence require
d(x, y)²
+ CMBI - d(I, Lid(z, v) 2 CRHK
2( ( 1)
21
49

---

Note that
d(r,y)
((x,y)²
Lad(x,)) + CMBI
2( 2<111)
2t
d( y)²
1
(L} log 2) Lid(a, y) + CMBI 2 log {
+ CMBI-
2
where in the last inequality we take d(a, y) = Ta log
Also note that when Ç is small,
1 is small. We can just take ORHK 2
CMBI CMBI-
2
Hence there exists constant ORHK
2 - s.t. for all E By(r).
1
f(z) log - fly) + log v) CRHK 2 28 ((2,))²
2. Step 2 Denote
VRHK(*)= exp(-f(xx) log un x. v) + f(y) log vz v. + CRHK)
exp( 1/4 ())
and V RRHK( (HK())) min (, VRHK (z)). Recall that the desired rejection sampling output can be
written as A X plt, x. 8) Euce. (a) On the other hand the actual rejection sampling
V(HK(()
output is > 11(1, 2r 4) Bpre.ams) (x)
Venec(() (w)
Following Fan al. 2023
Proof of Theorem
6]. we get
VRHK
VRHK E [[HHK]
|| x
S E exp( dim
E [R[K
[ VRHK] |
<
[[[[VHK]]
we only need to bound
2E
VRHK-
VRHK]
2 11
log v. + CRHK (
|B [[[]HH]]
In + (())
un v) + (RHK)dVV (x)
2 (ann())>11 exp(-f(x) + +
[B exp(-f(x) + fyy) =og
< 2 [(Vnn(())211 exp(Lid(x; 8) (1),,,)) dV)
JBLUS (xl
So it suffices to upper bound
S (Venn(()>11)
We need sharper bound for distant points. With 1 [] 0.5) log t.- have
1
-(() + f(y) + logu(n,
v. v) + d(x,v) y)²
2T
d(z,v)2
d(x, d(x, (I 0.5) /
SLnd(r, v) log Poly
21
d( v) 2 1 log
LLid(I, v) log Poly
C
+ log Poly (
log {
50

---

where in the last inequality we set d(x,y) = L, log Hence
2
I (ann
exp(Ld(a, u)) v, (2)
1V (Vando))))
1
+ log Poly |
exp( log ?
d( (V,
27
< Poly ) f
cxp(
(1 v) (z)
(Venx())=1
2T
Poly 7 (2)xpp(- ( d
d
0.5 log
1
d
=O(poly 1 (=)
10000000 0.3
O(S)
Here we used the fact that
- 1.
d²
Proposition 44. Let R v/2 Lvd The cost for rejection sampling as in Proposition 39 O(1) number
of rejections, in Ç and dimension.
Proof. [Proof of Proposition
Step 1: We first show that in a local neighborhood we have
d(x,
log (1, + log vn (1. 2, 2,2) S
+ Const
2n
We have that
v( 2 2,,,,, - -
-
c
expl
) e m0
11²
21
d(z, 2
c
expl
)1 1(),
172
21y
where C, = Lid so that 1 Cr log /
We have that for d(x, y) <
2
C,
exp(2 ( d(a,y) (
8(x, 1) = exp(
d(z.,) 17
2
1
c
7
2n
Thus we can assume W.L.O.G. that 8(2, 11) < Cs is small enough s.t. log(1 8(x, n1) is of
constant order, for all d(x, y)² < 2.
C. Then we have
d
d(x,z)²
(g. < - log c log 1) +
log(1 - 5(zz n))
2
2n
d
Mr,U)²
2
< 2
log 17 +
loge, vd(z, < C
2m
On the other hand, when x = y, have v 21))
C n
C
C
bo. S C
(1 x, x) <
+
+ e
17
n
51

---

observing that nिक S S Cy log = 1. Thus we have
log xI) S log C'
log 1).
2
Hence for C, we have
-log (n(1), J. v) + log (1 x.
d
d(x,y)²
2nj loge + log log = d(x,y)²
d
C c
2 C
2 lo +
log vd(x, -
21)
Step This step similar Proposition 48
With the first step, we see that
(()) [()) log in (1), 2r, 2. 2) log en (1) x) CRHK S Lid(a, v) + d(z, y)
C log
21j
where CRHK - log(2x) - log (C + 1) 2lo ? and c can be taken as
(2m)
With Lfdlog and T 1o(d-1)log viewing left hand side as quadratic function of d(x, v)
we get
Lıd(x,y) + d(r,y)² 21
d(x, y)²
S
2T
2 log
Thus have
((x) [(y) log (1, 2, + log (1, (1,,,,,)) xx CRHK
SLid(2, v) d(z, v)2 + log
C <
1
d(x,v)² ?
C
+ log
2n
2log
2T
The neighborhood for which the bound valid is d(x, u) S 2 i.e., a ball with radius R C, We have that R dd-1 O()). This allows us to do as exactly in Proposition 49 (with from Proposition 39
/2
L.dd
sin /
to show the rejection sampling procedure finishes with ((1) number of rejections. G
Proof. [Proof of Corollary
10 Using Pinsker's 's inequality, we have
1 I (xx 10))) V 2 ([++na)22
1 Hxx(p5)) -
ll x
11v xx112v
2
We want to bound l|px - x lrv I f.c. It suffices to have H.x(e)) (I+na)ax
Hence we need
log x Tog(1 (+)())))))))
2H x (o6))
log(
< 2k log(1 + na), i.e., = o
O L 1pa ). Hence k = O
log H (o))
For small step size we have log(1+na)
na
Ō om
log }).
Using Proposition 44 by setting 1 lld f the expected number of rejections rejection sampling is (1). We pick follows that
lod and consequently n Lıdlog Lidlog" + OILdlog !
ar
Ō Lid log"
1
k
o (mr)
log
The result then follows from triangle inequality:
1
X 1lvv < ao p lrv + lll x 1lev S k(SRHK + (MBB) + _ Ö(s)
where from Proposition 39 we can set the e to be S. so that KÇ = Ö ( 1.dd (log² !)) Tedlog
((c).
52

---

1.1.2 Heat kernel truncation: hypersphere
In this subsection, we show that on hyperspheres 8 the truncation error bound 1 VL - (c) can be achieved with truncation level I. ((Poly (log })) proved Zhao and Song [2018] the
heat kernel on St can be written as the following uniformly convergent series (with yp =- (2r, ())+))
expl k(k + d - 1)t
2k + d - 1 (dd-1)/ 11/2 (cos(s))),
vnn x. v)
ch
2
(d. 1)
where C? are the Gegenbauer polynomials. Define
1( M +
T( I(d- (Td-1)r((1) r( +1)
d(
Mr
Such M2 constructed to be an upper bound for Gegenbauer polynomials: see Zhao and Song [2018] Proof of Theorem 1]. The following proposition is directly implied by Zhao and Song] [0018 Theorem I]. and we provide proof for completeness.
Lemma 45. For I = Odd)), have
M+11 - M1
(1 (1+ -- 1)d-2
(1 11
and consequently M Mi = (1).
Proof. For 1 = (dd), we have (by definition of M and Gamma function) Mr = ((+d-2) and
(d-2)11
M4+1 ([+1+d-2)! (d-2)1((+1)** Hence Mj (1+d-2)!
((+1)=-2 and M411 (d-2)((+1) (d-2)
(++++d22)!
(++d-1)4-2
(d-2)!!)
(d-2).
Then we have (note that / ))d²))
M411 1++d-1)- (1+
d d- 2 2-2 =0((1 210==
MI 11
1+1
Proposition 46. Let M = S be a hypersphere. For truncation level L = OQd2 (log !)), we achieve [v(n, x. y) Vz pn. x, v)1 (()), E Sd.
Proof. [Proof of Proposition 46
Throughout the proof, we denote yp = (a., (r..)) The parameters
Mi satisfies |C, (x)1 < Mi according to Zhao and Song 2018 Proof of Theorem 1]. Hence, we
have
p1 v) v) - 1L 2, 2,8)1 8)1
2
((++d)
1)m
21 + d (d- 13/2 (d-1)// (cos(1)))
exp(
C
2
(d - 1)Ase
1=L+1
00
1(1+
1)m,
(21 + d 1)
exp(
2 (d 1)Asa
[=
Observe that for all > L + 1, since for large L (that depends on dimension) we have
Mr
O()):
Mi
see, also, Zhao and Song 2018 Proof of Theorem 1]. Hence,
Q exp(- (L+1 11m (21 + d + M)M441
21 + d+11
(d + 20), (d+2()m
M1+1
2
21 + - exp(
exp(
1(I+d-1mm (21 + 1)M
2
Mi
(21 + +d) M++1
=O(exp(-
Mi
(2L 2 + dyn O(exp(-Lx))) 0())
=O(exp(-
2
53

---

Algorithm 6 Riemannian Gaussian on hypersphere through rejection sampling
Input xx* E S.. Require S (d+1)2 for i = 0, 2, do
Generate proposal e x e क [eell (in Euclidean space), repeat until ||a|| S
Generate u uniformly on [0, 1].
|vl ("To peillely y
Return e uf I end for
Generate E' to be an orthonormal basis for Ty5", , set v + vo E T,S4. Output sample y = expx- (v)
For the last line, note that with L = Poly (log t) and rj = Clogg I. we have that Lrj = Poly (log t This implies exp(-Ln) ((exp(log ()) O()). Now we compute the truncation error.
[m x. v) (2 2.))
111 d 1)m, 21 + d -1
C (dd-1)/2 (cos((())
exp(
2
(d 1)
1=+11
1(1 + d - 1)m
exp(
(21 d 1)M,
(d 1)Ags
2
1=L+1
(L 1) + 1
1)17
(L + -
exp(
(2L + d
( 1)AS
2
(d (d2))! 1-Q
(1+ 1)(L, d)M
1(2 (M
(1 1++ 11)
0(epp(
d+
(d 22!
1 {)) =(()). =O(exp(- Poly (L)) = Ō(exp(- Poly (log 2) Poly (log
Poly
I.2 Sampling from the Riemannian Gaussian distribution hypersphere
We show that Riemannian Gaussian distribution on hypersphere can be generated efficiently through rejection sampling (Algorithm 6).
Note that Algorithm 6 first generate a Euclidean Gaussian v in the tangent space, which is guaranteed to satisfy ||인|| I т. Its density can be computed exactly under the normal coordinates. Then we perform rejection sampling to get samples v s.t. expa- (v) follow the Riemannian Gaussian distribution μ(l,, x²," exactly.
Proposition 47. On hypersphere Sd, when t. = O Algorithm output a sample following Riemannian Gaussian distribution p(t,a", ,u) x expl- 2 "1" .))²). with tteration complexity O(1).
Proof. First note that when generate tangent space Gaussian restricted to ||a||| under normal coordinates the corresponding density would be x e - lml x E Bb(T).
Recall that the Riemannian metric of Sá under normal coordinates satisfies / g = ( For Riemannian Gaussian distribution, we change it from Riemannian volume measure to the measure in local coordinates, and under the local coordinates, it has density
1 2 |x g = exp( bo 21
sin lar
x exp(
dd
x
Therefore we can compute the number of expected rejections as S e- |=| dx
SPo(a)) 6 1x12 siu 1z
do
54

---

Recall that with R
V
6 I+d
we know ( ) d-1 = =()). Then we have
e =
sin
d- dx
sin R
e to dr
Bo(x)
[
R
Bo(R)
Then we get
[Bo(m) e- | dx
d-1 S e- 1u1 dx
R
(BO(T) e -z |ulp sin Id= dx
sin R
Sno(Rs E- -dll da
(2xt) $
-01 (2Tt) exp( - 1
)
0(1)
R2 d)))
As long as 8 2 i,e, 5 (T+-)z the last equality holds.
1.3 Varadhan's 's asymptotics
We consider the approximation scheme introduced in Section A.2
using Varadhan's 's asymptotics. Let <p()) = 2 d(s 2) Intuitively, we want to see how the function yp can improve the convexity of ++
On manifold with positive curvature, we consider the situation that we cannot compute the minimizer of g(x) fxx + 2n (xx 3) and instead use y as the approximation of it. Notice that when n is small, since fx) is uniformly bounded, the function g(x) is dominated by 1 d(r,, y)?. thus the 2n minimizer of will be close to y. Therefore it is reasonable use approximation of the mode of e 9(2), Then in rejection sampling, we use jut, y, x) as the proposal.
Let L be the Lipschitz constant of In the next proposition, show that for some constant C.. with certain choices of 17 and f, it holds that
Ce ((z,))²
f(x) + d(zz, v)² f(v) +
2n
2 2t
Consequently, the acceptance rate defined by
V(x) exp(- f(x) + f(v) 1 d(x C2
2n
exp(
dd((z))))
2:
is guaranteed to be bounded by 1. Then, in Proposition
49
show that the expected number of
rejections is (1) in dimension d and step size
Proposition 48. Let be Li-Lisschitz and Cg be some constant. Take and T.. it holds that
C. T.zd
With T
C,
[.(d+1]
2T + +
1 2n
dla, d(x f(v) + C > d(x, 1)
2 21
Consequently, the acceptance rate is bounded by 1, i.e., V(x) < 1, Vr M.
Proof. [Proof of Proposition 48) Since f is L Li-Lsschhitz we have grad f(x)] < L,. Then we
have Lid(x,y) y) 2 fz) f(v) 2 L1 g).
1. The lower bound: The goal is to find some t and constant C such that
1 _dxx y) - f(v) + C d(x, u)
f(x) +
21
21
It suffices to find t, C such that
1
1
d(x, v²²
d(x, v)² L,d(æ, v) + C 2 0.
2n
2t
55

---

The left hand side can be viewed as a quadratic function of c(z,, y). When d(z, v) = -
the left hand side is minimized, and the mimimum is -
- 0 C. Hence can take
43
L2
Take r = C. and 1) Then we have -
1
Z
2. The upper bound: For an upper bound, we want some T S rj for which we want to show that
Ce S Lid(a,)) + d(x,u)² 2
Ce
1
f(x:) d(x,v)² f(v)
- d(x,v)²
2n
2n
2
2T
Similar before, suffices to show
_ - )(,, -
1
Ce
+ Lid(r,)) v)
SO.
2nj
2
The left hand side is maximized d( v)
b with maximum
L
C. 2
Take
T TT(44T) We can then verify that
L
C.
1 L
C. 2
0.
2
T2//c
3. Combining the two steps: From the above two steps, we get
27 d( C. f(x) d(a,, w)² d(z,
1 (w))
Cs
_ 21
1
d)x
21
In the following proposition, we show that on hypersphere (where the Riemannian metric normal coordinates is well studied), the expected number of rejections which equals to
SM exp( 1 dx, v) (v)
JA exp(-f(x) + f(y) 2n d(x,, y)² Co (Vy(2))
which is independent of dimension and accuracy.
Proposition 49. Let M be hypersphere. Set C2 log Assume without loss of generality that L1 2 max(1, 1, (7d) ). Then with 17 = T2a rejections is ((1) in both dimension and
1y for small e, the expected number of C. and t
Proof.
C₂ We try to bound the expected number of rejections. We compute it as Let T 13(d+1)'
follows:
SM expp d(x 8) (VVV (x)
JN exp( 1 d( )VV)(())))))))))))))))))))) 21
SM exp(pp ((x) + (wy) 27
I dx y² T (dV)(x)
JM exp(- 2/ d(. (V)())
Using Li and Erdogdu 2023 Lemma 8.2] and Li and Erdogdu 2023 Lemma C.5], when , 2 d using Kiemannian normal coordinates we have the tollowing lower bound on the integral:
R²
/ M
3 d(x,, d(x,("
B exp( d(z, y) (VV
expl
>
B,)R)
2
sin R d-1
2T y (1 (1-exp)
(1- expl 1 (3R² d))).
R
where BR(R) denote the geodesic ball centered at y with radius R. On the other hand, we have
/ expl exp(-2_dxx d(x, v)
(V)
/
1xp(-51101)
( )
2x
56

---

We next find a suitably small R which only depends on dimension, for which we have
R -
sin R
1 Using Taylor series for sin(R), R)), we have f
R- Hence for R
Rt
S
6
have
14d
(approximately)
R sin R
S 1++ 2. Consequently we set R = 1+d
and we know
R )d--
(1).
sin R
Combining the bounds discussed previously, have
JAr exp(- dz v) )dV,
JA exp( 27 v) Ce Cz)
R
(2T1)
dd-
sin R
(2TT) (1 (1-exp(- - 2/2
[i(++)
R²a - dd))
C.
Se C +1
t
1
e Cctl d
1
1 --exp(- 2 T
[1(d+1)
R² dd)
d-
exp( - 2 C.
[1(d+1)
R² d))
Co
For small s, we have C2 5 1. Since we assumed L1
2
1 and L3 >
d+1 5 -
we have
1- e [ia+1l) R² d)) 2 exp(- (d+1)
6
d))
2
expl
As a result,
d+T we see that the expect number of rejections is of order O(1):
SM exp( d(x, y) (
c²
d
++1
d S 20
d + d
SM exp( d(z, v) 1) VV (x) S
1 - exp(
2)
1
d
d
Observe that (d+1)* ) = (1+ (d- ( - +1) = 0(1)
57


Reviews:
[{"review_id": "EvS78iBmGx", "reviewer": "Reviewer_2wfT", "summary": "High-accuracy sampling on manifold is proposed by extending the work in Lee et al. [2021]. This algorithm is firstly discussed assuming the exact oracles are available, and then discussed under the inexact case.", "questions": "* What exactly $d$ in Line 292 refers to?\n* In Proposition~9, the truncation level is given in order for the inexact sample to achieve certain accuracy. I am wondering whether this should also relies on how fast the reject sampling converges?\n* In Line 144, what does the ``minimal solution'' exactly mean?\n* Overall, what properties should $f$ satisfy in order to satisfy the LSI and PI properties?\n* In Line 584, should it be the  Riemannian Langevin Monte Carlo Algorithm?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strength: The paper provides a high accuracy sampler on manifold, which is the first attempt towards this direction. The convergence of the sampler is firstly established under the exact MBI and RHK oracles, and then discussed based on the inexact oracles via heat kernel truncation. \n\nWeaknesses: \n* LSI somehow corresponds the  KL condition in optimization. In optimization, there is a well-established framework for the convergence analysis of the algorithms under the KL condition, even for Riemannian optimization. Similarly, after assuming LSI,  it seems the analysis follows a well-established route similar to the Euclidean case, even though there is a key property that needs to be carefully checked.\n * The paper also considers analysis under the inexact oracle via kernel truncation, it seems there are no experiments to validate the effectiveness of this scheme. Overall, despite the theoretically high accurate result, I guess how to do the sampling efficiently to satisfy the two oracles approximately is still very challenging for the manifold case.", "was_revised": true, "final_justification": "Extension from the Euclidean space to manifold is novel, and future work on further improving the practicality of the algorithm is interesting. I will maintain the score.", "rebuttal": "We thank the reviewer for their thoughtful questions and positive evaluation.\n\nResponse to weaknesses:\n\n\"LSI somehow corresponds the KL condition in optimization...\"\n\n - We would like to clarify that the extension is more involved than it might appear on surface. The proximal sampling algorithm was interpreted alternatively as running heat flow and its time reversal in [1], where a key property is that certain \"convexity\" (e.g., LSI) is preserved along the heat flow. This is easy to prove in the Euclidean space setting, as one can interpret heat flow as convolution with Gaussian. But on a Riemannian manifold, such a strategy doesn’t work. Please refer to our Proposition 31 for a theoretical justification, where heat flow was interpreted as a Markov operator, which is quite different from the Gaussian convolution treatment. It appears that this subtle yet important point may not have been fully considered in your assessment.\n - Also, as opposed to the Euclidean setting, due to the lack of closed-form solution for the heat kernel, both oracles are to be approximated to make the method practical. In our work, we provide high-accuracy rates under this simultaneous inexact oracles setup, which is also different (and relatively non-trivial) from the Euclidean setup.\n\n\"The paper also considers analysis under the inexact oracle via kernel truncation...\"\n\n - Thanks for the question. Actually, we did provide experiments for both practical approaches; the results were presented in Appendix B, page 16. The first two examples are for hyperspheres, where we implemented the truncation approach. The third one is for the manifold of positive definite matrices, and we implemented the Varadhan approach. Experimental results suggest that our proximal sampler (using both approaches) achieves a higher accuracy compared with the Riemannian Langevin Monte Carlo method.\n\nResponse to questions:\n\n\"What exactly $d$ in Line 292 refers to? \"\n - $d$ denotes the geodesic distance on Riemannian manifold. The notation was mentioned in line 81-82. Informally speaking, consider $S^{2} \\subseteq \\mathbb{R}^{3}$, like our earth. The geodesic distance between two cities is the length of a path that people need to travel. Here the path is \"curved\" (because the earth itself is curved), but not a Euclidean straight line.\n\n\"In Proposition 9, the truncation level is given in order for the inexact sample to achieve certain accuracy. I am wondering whether this should also relies on how fast the reject sampling converges?\"\n\n - Yes, in Proposition 9 we mainly consider the truncation level needed for a desired accuracy. On the other hand, the cost for rejection sampling (i.e., how fast rejection sampling converges), is an open problem in general. But we are able to give an theoretical analysis for the cost of rejection sampling when $M = \\mathcal{S}^{d}$, see Remarks at line 335 and Corollary 10.\n\n\"In Line 144, what does the \"minimal solution'' exactly mean? \"\n\n- Roughly speaking, if there is another function that also solves the heat equation, then it must be greater than or equal to the heat kernel. See [2](Theorem 4.1.5, Proposition 4.1.6) for more details.\n\n\"Overall, what properties should $f$ satisfy in order to satisfy the LSI and PI properties? \"\n\n - Thanks for the intriguing question. For the Euclidean setting, recent works translated LSI/PI conditions on $\\pi$ to conditions on $f$ like the Polyak-Lojasiewicz condition. These works considered the LSI/PI constant for $e^{-f(x)/t}$. For example, in the $t \\to 0$ limit, i.e., \"low temperature\" setting, [3] considered LSI constant and [4] considered PI constant. Another work [5] considered a \"Optimizability\" condition and analyzed LSI/PI constant for (informally) $t \\le O(1/d)$. We believe that similar conditions could be derived in the manifold setting, although there might be more sophisticated technical challenges to handle (and we plan to consider this in our future work, thanks to your question). We also emphasize that even in the Euclidean setting this question is not fully resolved.\n\n\"In Line 584, should it be the Riemannian Langevin Monte Carlo Algorithm? \"\n\n - Yes. Sorry for the typo, we will correct it.\n\nWe thank you again for your insightful feedback and questions. If our response above satisfies your concerns, we would greatly appreciate an increase of score to reflect your increased confidence in our work.\n\n[1] Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono. Improved analysis for a proximal algorithm for sampling. In Conference on Learning Theory, pages 2984–3014. PMLR, 2022.\n\n[2] Elton P Hsu. Stochastic analysis on manifolds. Number 38 in Graduate Studies in Mathematics. American Mathematical Soc., 2002.\n\n[3] Sinho Chewi and Austin J Stromme. The ballistic limit of the log-sobolev constant equals the Polyak–Łojasiewicz constant. arXiv preprint arXiv:2411.11415, 2024.\n\n[4] Yun Gong, Niao He, and Zebang Shen. Poincare inequality for local log-Polyak-Łojasiewicz measures: Non-asymptotic analysis in low-temperature regime. arXiv preprint arXiv:2501.00429,2024.\n\n[5] August Y Chen and Karthik Sridharan. Optimization, isoperimetric inequalities, and sampling via Lyapunov potentials. arXiv preprint arXiv:2410.02979, 2024."}, {"review_id": "l5FC9GVnuk", "reviewer": "Reviewer_mAJz", "summary": "The paper proposes a Riemannian extension of the proximal sampler [1, 2]. Instead of the Gaussian Oracle and Restricted Gaussian Oracle used in the standard proximal sampler, the paper introduces the Manifold Brownian Increment (MBI) oracle and the Riemannian Heat Kernel (RHK) oracle (Section 3), and derives convergence guarantees (Section 4). Unlike in the Euclidean case, these oracles may not be exact. Therefore, the paper proposes approximate implementations of these oracles (Section 5), which still guarantee convergence (Section 4).\n\n\n[1] Yin Tat Lee, Ruoqi Shen, Kevin Tian. Structured Logconcave Sampling with a Restricted Gaussian Oracle. Proceedings of Thirty Fourth Conference on Learning Theory, 2021.\n\n[2] Yongxin Chen, Sinho Chewi, Adil Salim, Andre Wibisono, Improved analysis for a proximal algorithm for sampling. Proceedings of Thirty Fifth Conference on Learning Theory, 2022.", "questions": "* I am not sure whether Lines 601–602 are correct. Does the function satisfy Assumption 3.3 of [5]? At a glance, $\\lambda_*$ appears to be possibly negative. If so, Theorem 3.4 of [5] would not apply. It is unclear to me whether $f$ actually satisfies the LSI.\n\n* Related to the previous question: as mentioned in the paper, there exists no globally supported non-constant strongly log-concave density on compact manifolds. More generally, it is conjectured that the convexity and smoothness of the function are closely related to the curvature of the manifold (see the paragraph before Section 5 in [4]). In this regard, I suspect that the existence of densities satisfying the LSI is not always guaranteed and may require some curvature condition on the manifold. Do the authors have any thoughts on this? Theorems 28 and 31 provide some partial answers, but I wonder whether there are any other known examples.\n\n* In Theorems 28 or 31, if $\\kappa < 0$, then $\\alpha < 0$. In this case, Theorem 6 appears to diverge, even if the step size is controlled, since the denominator becomes less than 1. Is this divergence inevitable?\n\n[5] Mufan Bill Li, Murat A. Erdogdu. Riemannian Langevin Algorithm for Solving Semidefinite Programs. Bernoulli 29(4), 2023.", "limitations": "The authors expressed limitations of the work. I have stated additional limitations that I found in the weaknesses section.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n* To the best of my knowledge, this is the first paper to introduce a Riemannian version of the proximal sampler. In particular, the method works under a more general setting (e.g., less smoothness, LSI) than previous Riemannian samplers [3]. The fact that the proposed method does not require $L$-smoothness of the potential is favorable, especially since verifying such a condition on Riemannian manifolds is less straightforward (and is, in fact, conjecturally related to the curvature of the manifold [4]).\n\n* While the main proof follows a similar structure to the Euclidean case (Theorem 6), extending the canonical assumptions in the Euclidean space to the Riemannian setting appears to be nontrivial and involve additional technical challenges (e.g., Proposition 31).\n\n* The fact that the proposed approximate algorithm satisfies the required assumptions (Section 5.2) and still guarantees convergence is a favorable result. \n\n### Weaknesses \n\n* Minor comment: In Section 2.1, Line 133, it would be helpful to specify the metric being used. Assigning a different metric (e.g., the Bures–Wasserstein metric) on the same space can result in a positively curved geometry.\n\n* The analysis relies on the LSI condition of the distribution on the manifold. However, it is unclear how strict or mild this condition is in the context of Riemannian manifolds. It would have been helpful to illustrate how strict or mild this assumption is by providing examples or criteria that can build the intuition (see also the below question section).\n\n[3] Xiang Cheng, Jingzhao Zhang, Suvrit Sra. Efficient Sampling on Riemannian Manifolds via Langevin MCMC. Proceedings of the 36th International Conference on Neural Information Processing System, 2022.\n\n[4] Foivos Alimisis, Antonio Orvieto, Gary Becigneul, Aurelien Lucchi. A Continuous-time Perspective for Modeling Acceleration in Riemannian Optimization. Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics. 2020.", "was_revised": true, "final_justification": "As discussed in the initial review, this paper has some interesting results.\n\n1. To the best of my knowledge, this is the first paper to introduce a Riemannian version of the proximal sampler. In particular, the method works under a more general setting (e.g., less smoothness, LSI) than previous Riemannian samplers [3]. The fact that the proposed method does not require $L$-smoothness of the potential is favorable, especially since verifying such a condition on Riemannian manifolds is less straightforward (and is, in fact, conjecturally related to the curvature of the manifold [4]).\n\n2. While the main proof follows a similar structure to the Euclidean case (Theorem 6), extending the canonical assumptions in the Euclidean space to the Riemannian setting appears to be nontrivial and involve additional technical challenges (e.g., Proposition 31).\n\n3. The fact that the proposed approximate algorithm satisfies the required assumptions (Section 5.2) and still guarantees convergence is a favorable result. \n\nThe main weakness of the paper is the LSI condition on Riemannian manifolds, which is not direct and has to be analyzed case-specific. While there is a criterion, such criterion has a limited usage. I believe this is an open problem needs not be addressed directly in this work and require the further study.\n\nIn conclusion, I will keep my initial rating (accept).", "rebuttal": "We thank the reviewer for their thoughtful questions and positive evaluation.\n\nResponse to weaknesses: \n\n\"Minor comment: In Section 2.1, Line 133, it would be helpful to specify the metric being used. Assigning a different metric (e.g., the Bures–Wasserstein metric) on the same space can result in a positively curved geometry.\"\n\n - Thanks for the comment, we will specify the metric. We are using the affine invariant metric.\n\n\"The analysis relies on the LSI condition of the distribution on the manifold. However, it is unclear how strict or mild this condition is in the context of Riemannian manifolds. It would have been helpful to illustrate how strict or mild this assumption is by providing examples or criteria that can build the intuition (see also the below question section).\"\n\n - Thanks for the suggestion, we will add more examples to clarify. In general, for convex potentials, the celebrated Bakry-Emery condition provides general results for a density to satisfy LSI based on the spectrum of the Hessian of the potential $f$ and the Ricci curvature of the manifold. The question is wide open in the case when the potential is not convex.\n - For the Euclidean setting, recent works like [2, 3, 4] for some recent works of translating LSI/PI conditions on $\\pi$ to conditions on $f$ like the Polyak-Lojasiewicz condition. We believe that similar conditions could be derived in the manifold setting, although there might be more sophisticated technical challenges to handle. We also emphasize that even in the Euclidean setting this question is not fully resolved.\n\nResponse to questions:\n\n\"I am not sure whether Lines 601–602 are correct...\"\n\n - Let $f(x) = -\\kappa \\mu^{T} x$. Clearly $x = \\frac{\\mu}{||\\mu||}$ is the minimizer, and $x = -\\frac{\\mu}{||\\mu||}$ is the maximizer. Note that they are the only two critical points. To compute the Riemannian Hessian, we construct a smooth extension of $\\text{grad} f$. Let $G(x) = - (I - x x^{T}) \\kappa \\mu = - \\kappa \\mu + x x^{T} \\kappa \\mu$ be a vector field defined on $\\mathbb{R}^{d+1}$.\nThe Euclidean directional derivative of $G$ is \n\\begin{align*}\n    DG(x)[u] = u x^{T} \\kappa \\mu + x u^{T} \\kappa \\mu\n    = \\kappa (x^{T} \\mu) u + \\kappa (u^{T} \\mu) x\n\\end{align*}\nHence when $u \\in T_{x}\\mathcal{S}^{d}$, we have \n$\\text{Hess} f(x)[u] = \\kappa (x^{T} \\mu) u$. \nRecall that the only critical point that is not a global minimum is $x = - \\frac{\\mu}{||\\mu||}$. (Actually, it is the maximizer.) We can get (for all unit tangent vector $u$), the Riemannian Hessian satisfies $\\text{Hess} f(x)[u, u] = \\kappa (x^{T} \\mu) u^{T} u = \\kappa (x^{T} \\mu) = - \\kappa ||\\mu||$. We see that the Hessian has negative eigenvalue, and hence $\\lambda_{*}$ is positive.\n\n\"Related to the previous question: as mentioned in the paper, there exists no globally supported non-constant strongly log-concave density on compact manifolds...\"\n\n - Indeed, as mentioned before, this is related to the curvature-dimension conditions (e.g., Bakry-Emery), which is discussed for example in [1](Proposition 5.7.1). For convex potentials, the LSI is satisfied baesd on the Ricci curvature of the manifold. More generally, it is open with recent explorations in the Euclidean setting, as mentioned previously.\n\n\"In Theorems 28 or 31, if $\\kappa < 0$, then $\\alpha < 0$...\"\n\n - In our notation, $\\alpha = \\frac{\\kappa}{1-e^{-\\kappa t}}$. When $\\kappa < 0$, we have $e^{\\kappa t} > 1$, hence $\\alpha > 0$. For example, if $\\kappa = -1$, we have $\\alpha = \\frac{-1}{1-e^{t}} > 0$.\n\nWe thank you again for your intriguing questions and careful evaluation.\n\n[1] Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and geometry of Markov diffusion operators, volume 103. Springer, 2014.\n\n[2] Sinho Chewi and Austin J Stromme. The ballistic limit of the log-sobolev constant equals the Polyak–Łojasiewicz constant. arXiv preprint arXiv:2411.11415, 2024.\n\n[3] Yun Gong, Niao He, and Zebang Shen. Poincare inequality for local log-Polyak-Łojasiewicz measures: Non-asymptotic analysis in low-temperature regime. arXiv preprint arXiv:2501.00429,2024.\n\n[4] August Y Chen and Karthik Sridharan. Optimization, isoperimetric inequalities, and sampling via Lyapunov potentials. arXiv preprint arXiv:2410.02979, 2024."}, {"review_id": "rWcBKCBwyK", "reviewer": "Reviewer_dci1", "summary": "This article generalizes an existing Euclidean proximal sampler to Riemannian manifolds, where the performance of the sampled, measured by its convergence rate, is established. The proposed sampler is based on two key oracles MBI and RHK, aiming to sampling the exact Brownian motion on manifolds. The practical implementations of the oracles are discussed, and analyzed. Preliminary numerical results, show the effectiveness of the proposed method.", "questions": "-\tWhy in Theorem 8, there is no distinction between positive and negative curvature situations, as in Theorem 6?\n-\tAlgorithm 2 and 3, when it is possible to find a suitable t, i.e. to ensure the existence of t and C_RHK or C_MBI such that … it is unclear whether C_RHK and C_MBI are negative. \n-\tThe idea of Varadnan asymptotics makes sense when eta -> 0, according to Fact 5. However, the rates of inexact oracles rely on large eta in Assumption 1, therefore it is conceptually contradictory to use Varadnan asymptotics in a large eta regime. A clearer explanation about this is needed to ensure that these results are correct. \n\n-\tClarify the space of g in Definition 4. For all g in which space?\n-\tThe notations of pi^Y|X and pi^X|Y in Algorithm 1 are not consistent with Section 4. Are they the same as pi_eta^Y|X and pi_eta^X|Y in Assumption 1?\n-\tClarify the notation tilde O (difference to O) used in several places (e.g. thereom 8).", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strength: The article is mostly well written and easy to follow. The main results in Theorem 6 and 8 are easy to understand, though they seem to be highly technical (and novel) to prove. The nature of the results is of theoretical interest, notably Corollary 10 gives a concert case to show the implication of the main results. \n\nWeakness: The idea of Varadhan’s asymptotics mentioned in Section 6 is not very intuitive. A better explanation is needed (see questions below). Some notations should be unified to better understand the results. Due to the lack of space, it is a pity that all the numerical results are given in Appendix.", "was_revised": true, "final_justification": "The rebuttal by Authors has clarified all my questions. I maintain my score .", "rebuttal": "We thank the reviewer for their thoughtful questions and positive evaluation.\n\nResponse to weaknesses:\n\n\"Some notations should be unified to better understand the results. Due to the lack of space, it is a pity that all the numerical results are given in Appendix.\"\n\n - Thanks for the suggestion. We will reorganize and try to move some experimental results to the main draft as much as the space allows.\n\nResponse to questions:\n\n\"Why in Theorem 8, there is no distinction between positive and negative curvature situations, as in Theorem 6? \"\n\n - Thanks for pointing out this. Yes there should be a slight distinction, where in case of negative curvature we still need an curvature dependent upper bound for the step size $\\eta$. But the iteration complexity will remains to be the same. \n    The revised theorem look like this:\n - [Theorem 8] Similar to Theorem 6, let $M$ be a Riemannian manifold without boundary. Assume Assumption 1 holds. For any initial distribution $\\rho_{0}^{X}$, to reach $\\tilde{\\mathcal{O}}(\\varepsilon)$ total variation distance with oracle accuracy $\\zeta = \\zeta_{\\mathsf{RHK}} = \\zeta_{\\mathsf{MBI}} = \\frac{\\varepsilon}{\\log^{2} \\frac{1}{\\varepsilon}}$\n    and step size $\\frac{1}{\\eta} = \\tilde{\\mathcal{O}}(\\log \\frac{1}{\\varepsilon})$ (for negative curvature, we additionally require $\\eta \\le 1/|\\kappa|$), \n   -  if $\\pi^{X}$ satisfies $\\alpha$-$\\mathsf{LSI}$, we need $k = \\tilde{\\mathcal{O}}(\\log^{2} \\frac{1}{\\varepsilon})$ iterations.\n   -  if $\\pi^{X}$ satisfies $\\alpha$-$\\mathsf{PI}$, we need $k = \\tilde{\\mathcal{O}}\\left(\\log \\frac{1}{\\varepsilon} \\log \\frac{\\chi_{\\pi^{X}}^{2} (\\rho_{0}^{X})}{\\varepsilon^{2}} \\right) = \\tilde{\\mathcal{O}}\\left(\\log^{2} \\frac{1}{\\varepsilon}\\right)$ iterations.\n\n\"Algorithm 2 and 3, when it is possible to find a suitable t, i.e. to ensure the existence of t and $C_{RHK}$ or $C_{MBI}$ such that … it is unclear whether $C_{RHK}$ and $C_{MBI}$ are negative. \"\n\n - We analyzed this theoretically in the appendix, please see for example Proposition 39.\n\n\"The idea of Varadnan asymptotics makes sense when $\\eta \\to 0$, according to Fact 5. However, the rates of inexact oracles rely on large eta in Assumption 1, therefore it is conceptually contradictory to use Varadnan asymptotics in a large eta regime. A clearer explanation about this is needed to ensure that these results are correct.\"\n\n - Yes, intuitively, when $\\eta$ is large, Varadhan's approach could lead to a larger error. An error analysis for Varadhan approach is still open. Our theoretical high-accuracy results are for truncation approach. But we are able to experimentally verify the effectiveness of this approach.\n\n\"Clarify the space of g in Definition 4. For all g in which space?\"\n\n - $g$ can be any smooth function on $M$.\n\n\"The notations of $\\pi^{Y|X}$ and $\\pi^{X|Y}$ in Algorithm 1 are not consistent with Section 4. Are they the same as $\\pi_{\\eta}^{Y|X}$ and $\\pi_{\\eta}^{X|Y}$ in Assumption 1?\"\n\n - Yes. See line 170: When there is no ambiguity, we omit the step size $\\eta$ ...\n\n\"Clarify the notation tilde O (difference to O) used in several places (e.g. thereom 8).\"\n\n - Yes, we will clarify the tilde O notation in our paper. In Theorem 8, we used the tilde O notation to emphasize that we are only considering the $\\varepsilon$ dependency.\n\nWe thank you again for your insightful feedback and questions. If our response above satisfies your concerns, we would greatly appreciate an increase of score to reflect your increased confidence in our work."}, {"review_id": "O2G5jWLrQe", "reviewer": "Reviewer_m8ed", "summary": "This paper considers sampling from distributions of the form $\\exp(-f)$ over spherical manifolds. Under the high-accuracy regime, the authors provide a Riemannian proximal sampler that produces a sample whose law is $\\epsilon$-close to the target distribution in TV distance within an iteration complexity polylogarithmic in $\\epsilon$. It is well-known that the major difficulty in manifold Langevin-based sampling is the intractability of simulating manifold Brownian motion, as the kernel generally lacks a closed-form expression. The authors overcome the issue by truncating the heat kernel and show that an ideal proximal sampler achieves polylog complexity with an inexact proximal operator and the truncated manifold Brownian motion. They further instantiate the algorithm using rejection sampling and show that the implementation preserves the same iteration complexity.", "questions": "Futher questios are already mentioned in the weakness part.", "limitations": "The paper is majorly theoretical and does not have negative societal impact.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n\nOverall, this paper is well-written and easy to follow. It should be the first to show that an ideal algorithm can achieve linear convergence by incorporating the inexact proximal operator on a Riemannian manifold. The authors demonstrate strong mathematical maturity and solid knowledge of both differential geometry and numerical methods. The theoretical analysis appears to be rigorous, and the proofs are generally convincing.\n\n**Weaknesses**\n\nHowever, I remain dubious about the novelty and the broader significance of this work. First, the results seem to be a relatively easy and straightforward extension of proximal sampling techniques from the Euclidean setting. The adaption is perhaps expected, and might limit the originality of the results.\n\nA more serious concern arises from a practical perspective. The implementation of the Riemannian proximal sampler relies crucially on rejection sampling to correct the additional bias introduced by the inexact oracles. But, it is folklore that rejection sampling can be computationally intractable in real-world and high-dimensional applications as it might require an exponential number of proposals in dimension. This brings serious issues and may cause the \"practical implementations\" to be totally infeasible in practice. It might be possible that the inexact method can be made more practical by considering other common techniques like MH used in high-accucary sampling, but the authors do not consider this or offer an explanation. Moreover, it is okay to be without an empirical section for a purely theoretical paper. However, the total lack of experiments is problematic when the practicality of the algorithm is not fully justified by theoretical results alone. Therefore, I would recommend rejection unless these concerns are addressed.", "was_revised": true, "final_justification": "I changed my score to positive after the discussion with the authors, as I acknowledge that they have addressed the issue of impracticality to an acceptable extent.", "rebuttal": "We thank the reviewer for taking the time to engage with our work and for their comments. While we respectfully disagree with some aspects of the assessment (main one being ``total lack of experiments''), we appreciate the opportunity to clarify the points of misunderstanding.\n\nResponse to weaknesses:\n\n\"However, I remain dubious about the novelty and the broader significance of this work. First, the results seem to be a relatively easy and straightforward extension of proximal sampling techniques from the Euclidean setting. The adaption is perhaps expected, and might limit the originality of the results.\"\n\n - Actually, the extension is not as easy as it looks. The proximal sampling techniques was interpreted as alternatively run heat flow and its time reversal [1], where a key property is that certain \"convexity\" (e.g., LSI) is preserved along the heat flow. This is easy to prove in the Euclidean space setting, as one can interpret heat flow as convolution with Gaussian, but on a Riemannian manifold, such strategy doesn't work. See our Proposition 31 for a theoretical justification, where heat flow was interpreted as a Markov operator, which is quite different from the Gaussian convolution treatment.\n\n\"A more serious concern arises from a practical perspective. The implementation of the Riemannian proximal sampler relies crucially on rejection sampling to correct the additional bias introduced by the inexact oracles. But, it is folklore that rejection sampling can be computationally intractable in real-world and high-dimensional applications as it might require an exponential number of proposals in dimension. \"\n\n - In general, a rejection sampling algorithm would result in an exponential number of proposals in dimension. However, this can be avoided. In Euclidean space, the proximal sampler also requires rejection sampling to implement the RGO step, see [1]. But by carefully choosing a step size $\\eta$ and constructing a suitable proposal, the number of proposals is $O(1)$ in dimension. We did theoretically justify this, see Proposition 49 and Proposition 44. Roughly speaking, we showed that for $M = \\mathcal{S}^{d}$, the cost for rejection sampling is of constant order, considering the dimension.\n\n\"This brings serious issues and may cause the \"practical implementations\" to be totally infeasible in practice. It might be possible that the inexact method can be made more practical by considering other common techniques like MH used in high-accucary sampling, but the authors do not consider this or offer an explanation. Moreover, it is okay to be without an empirical section for a purely theoretical paper. However, the total lack of experiments is problematic when the practicality of the algorithm is not fully justified by theoretical results alone. Therefore, I would recommend rejection unless these concerns are addressed.\"\n\n - Please note that we had provided simulation experiments (for both implementations) in Appendix B, see pages 15 and 16.\n\n - Furthermore, based on your questions, we tested the practical applicability by sampling from $\\mathcal{S}^{100}$, and the experiment result shows that the (average) number of proposals is not exploded by dimension. We will add this new experimental result to the revised draft as well.\n\n Thank you again for your feedback and questions. We sincerely hope that our response has addressed your concerns and may lead you to reconsider your evaluation. Should our clarifications increase your confidence in the work, we would be truly grateful if you felt it appropriate to adjust the score accordingly.\n\n[1] Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono. Improved analysis for a proximal algorithm for sampling. In Conference on Learning Theory, pages 2984–3014. PMLR, 2022."}]

Author rebuttal:

Meta-review: Reviewers unanimously find merit in this work, praising its presentation quality (Reviewer dci1, Reviewer m8ed), technical soundness (Reviewer dci1, Reviewer m8ed), novelty in the sense of the setting being minimally explored (Reviewer 2wfT). The main weakness raised is by Reviewer m8ed, who worries that the theory is too similar to known Euclidean proximal sampling theory. Reviewer m8ed, based on their score, does not view this concern as too significant, and neither do I: non-Euclidean theory developed in other domains such as optimization (for instance of geodesic convexity) is interesting, and has resulted in new phenomena not necessarily present in the Euclidean case, which may well be the case here as the issues become better understood.

On basis of reviewer consensus, I recommend this work is accepted.
Decision: Accept (poster)