Fast Non-Log-Concave Sampling under Nonconvex Equality and Inequality Constraints with Landing
Authors: Kijung Jeon
Michael Muehlebach
Molei Tao
Venue: NeurIPS 2025 poster ()
Keywords: constrained sampling
nonconvex manifold
overdamped Langevin dynamics
Primary area: probabilistic_methods
TL;DR: We propose OLLA, a projection-free overdamped Langevin framework that enforces both equality and inequality constraints via a deterministic “landing” correction along the manifold normal.

Abstract: Sampling from constrained statistical distributions is a fundamental task in various fields including Bayesian statistics, computational chemistry, and statistical physics.
This article considers the cases where the constrained distribution is described by an unconstrained density, as well as additional equality and/or inequality constraints, which often make the constraint set nonconvex. Existing methods for nonconvex constraint set $\Sigma \subset \mathbb{R}^d$ defined by equality or inequality constraints commonly rely on costly projection steps. Moreover, they cannot handle equality and inequality constraints simultaneously as each method only specialized in one case. In addition, rigorous and quantitative convergence guarantee is often lacking.
In this paper, we introduce Overdamped Langevin with LAnding (OLLA), a new framework that can design overdamped Langevin dynamics accommodating both equality and inequality constraints. The proposed dynamics also deterministically corrects trajectories along the normal direction of the constraint surface, thus obviating the need for explicit projections. We show that, under suitable regularity conditions on the target density 
and $\Sigma$, OLLA converges exponentially fast in $W_2$ distance to the constrained target density $\rho_\Sigma(x) \propto \exp(-f(x))d\sigma_\Sigma$.
Lastly, through experiments, we demonstrate the efficiency of OLLA compared to projection-based constrained Langevin algorithms and their slack variable variants, highlighting its favorable computational cost and reasonable empirical mixing.

Fast Non-Log-Concave Sampling under Nonconvex Equality and Inequality Constraints with Landing
Kijung Jeon
Michael Muehlebach
Molei Tao Georgia Institute of Technology mtao@gatech. edu
Georgia Institute of Technology xjeon@gatech. edu
MPI-IS mi chaelm@tue www.imittuueeee.....mmm meg.de
Abstract
Sampling from constrained statistical distributions is a fundamental task in various fields including Bayesian statistics, computational chemistry, and statistical physics. This article considers sampling from a constrained distribution that is described by an unconstrained density, as well as additional equality and/or inequality con- straints, which often make the constraint set nonconvex. Existing methods struggle in the presence of such nonconvex constraints, they rely on projections, which are computationally expensive or intractable, are specialized to either inequality or equality constraints, and often lack rigorous quantitative convergence guaran- tees. In this paper, we introduce Overdamped Langevin with LAnding (OLLA), a new framework that can design overdamped Langevin dynamics accommodating both nonlinear equality and inequality constraints. The proposed dynamics also deterministically corrects trajectories along the normal direction of the constraint surface, thus obviating the need for explicit projections. We show that, under suit- able regularity conditions on the target density and the feasible set E C Rd OLLA converges exponentially fast in 2- 2-Wasserstein distance to the constrained target density (x) x (=p(--(=))ddr::: Lastly, through experiments, we demonstrate the efficiency of OLLA compared known constrained Langevin algorithms and their slack variable variants, highlighting its favorable computational cost and fast empirical mixing.'
1 Introduction
Sampling from complex, constrained statistical distributions is fundamental problem in machine learning, with applications ranging from Bayesian inference under structured priors to training generative models with safety or fairness constraints. When there is no constraint, a prominent class of sampling techniques centered around (overdamped) Langevin dynamics, where the drift is set to the gradient of the log-target density. These have gained significant traction due to their strong theoretical guarantees: for example under log-concave target densities [1-8] or more general densities that satisfy relaxed conditions such isoperimetric inequalities [9, 10] one obtain fast, non-asymptotic convergence rates. Langevin dynamics can even be generalized, via a Riemannian Langevin approach, to sample under convex constraints [e.g., 11, 12]. However, extending Langevin- based approaches to sample from distributions supported on nonconvex sets Σ C Id remains a major challenge. Existing techniques typically rely projection steps, which computationally expensive and require convexity to ensure convergence. Moreover, most methods offer limited quantitative convergence guarantees in this case; in fact, even if the density is log-concave a nonconvex constraint can easily make the target distribution much harder to sample from, also rendering the analysis more difficult. This is critical bottleneck for many emerging machine learning
1 All code is provided in following repository:
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

applications such as imitation learning [13] or constrained generative modeling - where the feasible set implicitly defined by complex equality and inequality constraints.
In this article, we introduce OLLA (Overdamped Langevin with LAnding), a suite of projection-free stochastic dynamics that could serve as the foundation for sampling from constrained distributions. OLLA avoids projections by combining two key ideas: (1) relying on local approximations of the feasible set to guide the sampling and (2) introducing a restitution mechanism called "landing" that guarantees convergence to the feasible set. These techniques build on ideas previously developed in the context of nonconvex optimization, where they have been shown to provide scalable and effective algorithms [14-16], and share close connection to several powerful constrained sampling approaches in the literature [17-21]. We adapt and extend the ideas to the sampling setting, which makes OLLA both computationally efficient and theoretically grounded. Our main contributions are summarized as follows:
(Unified treatment of constraints) OLLA is described by a stochastic differential equation (SDE) that, in contrast to related prior works, enforces both equality and inequality constraints. This is achieved by constructing the tangent space to the constraint manifold and projecting the overdamped Langevin drift and diffusion terms onto the tangent space resulting in a simple least- squares problem. OLLA recovers the classical cquality-only constrained Langevin dynamics special case, yet seamlessly accommodates arbitrary smooth inequality constraints without resorting to slack variables or projections.
(Exponential convergence) We prove that the continuous version of OLLA converges to the constrained target distribution 3e at exponential rate under appropriate regularity assumptions, Our convergence results are non-asymptotic and characterized by the 2- 2-sassesttiin distance in all scenarios (equality constraints only, inequality constraints only, and mixed).
(Efficient SDE discretization with trace-estimation) We introduce OLLA-H, a computationally efficient Euler-Maruyama (EM) discretization of the aforementioned SDE that features a Hutchin- son trace estimator [22] for approximating the Itô-Stratonovich correction term arising from the diffusion. As result, OLLA-H has low computational cost per iteration, even in high dimensions, and it achieves relatively accurate sampling and empirical mixing, an aspect that we demonstrate in various numerical experiments.
2 Related Works
We first focus our literature review on recent works that consider Langevin sampling under nonlinear constraints. One of the closest touching points is [17], which describes a gradient descent approach on the KL divergence. The algorithm shares some similarities to ours in that projections are avoided, but the work focused on equality constraints only, whereas OLLA covers both equality and inequality constraints. The work [23] proposes the use of slack variables to incorporate inequality constraints to change a mixed problem into an equality-only problem, which comes at the cost of additional spurious dimensions. In similar vein, [18] designs particle-based variational inference method to incorporate incquality constraints. The method is effective at sampling under inequality constraints only, suffers, however, from high computational cost in high dimensions due to the estimation of associated boundary integrals.
OLLA is inspired by recent methods in nonlinear optimization [15, 14, 16) that use a similar landing mechanism and avoid projections onto the feasible set. There has also been important work by [19-21] who introduced constrained Langevin dynamics based numerical schemes such as SHAKE, RATTLE [24-26], and including
Hitropolis-Hasttngg corrections [21]. These works mainly focus on equality-only constraints, although inequality constraints can be incorporated via including slack variables or applying reflection at the boundary. We use these algorithms as basclines and refer them to Constrained Langevin (CLangevin) [20], Constrained Hamiltonian Monte Carlo (CHMC) [20], and Constrained generalized Hybrid Monte Carlo (CGHMC) [21].
In the present work, constraints are handled through careful decomposition of the stochastic dynamics on the boundary into normal and tangential parts, thereby avoiding projections and even enabling infeasible initialization. There have also been alternative algorithm designs that, however, do not share these features, for example, [27, 28] (based on projection). [29, 30] (barrier functions), [31] (reflections), or [11, 32, 12] (mirror maps). Other works by [33, 34] introduce penalties for constraint violations or relax the notion of constraint satisfaction [35]. In addition, we note that, although closely related, constrained sampling is not the same sampling on manifolds [36-39].

---

Preliminaries & Notations
We consider sampling from a target density supported the compact and connected Riemannian
submanifold of IRd defined by Σ == {a Rd (x)) = 0, g(())) v 0} where h : RRa - I and
g R + R' are smooth functions. To guarantee that all constraint-related constructions are well- posed, we impose the Linear Independence Constraint Qualification (LICQ) condition [40].
Definition 1. The functions h, satisty LICQ (4), (V))))( are linearly independent for every E 2, where La denotes the set of active inequality constraints, i.c., 1 = {i E [4] (x(x) 2 0}.
As result of LICQ, the tangent space of Σ at x E E can be defined as
TE = {ur E Rd Vh(z)u = 0, (g - 0, 12)
and its orthogonal projector onto T& is (xx) I VJ(x)
where J(x)
(1,, mm (I), Dr, (x) + e, Sa (x) +e| + .f,,,.., = 1 denotes the stacking of con-
straint functions for some
and G(x)) VJ(x) is its associated Gram matrix.
For any smooth f Σ + R and a smooth vector field X - Rd the intrinsic gradient and divergence on Σ are given by Vzf(a) [(z)) divgX( = (I(z))VX())) where V denotes the usual Euclidean gradient Jacobian in IRd Our target density is set to be p2(x)) xx (ep(-f(x))doss on E with day: being the induced Hausdorff measure on E. We write Pt, Pr. be the law of OLLA and the projected process of OLLA onto E at time On E, the natural extension of KL. divergence and Fisher information take the form:
KL <pl(p)) = -s pln
(
P5 (0'dos doy and (Pllpps) = 1 pl In ( 12dox
Ps
We now streamline notations appearing in Section 4. A complete list of symbols and precise technical definitions are included in Appendix A. In particular, let (Ix) denote the nearest-point (Euclidean) projection onto Σ, and let λy.I be the log-Sobolev constant of (E, px.). We then assume the existence of the following constants.
Mn
sup Zo (pu
[1((r))
Mg := sup
(1)))))
za Esu
over the support of the initial law Po- The constant x (Lemma C.1, Lemma C.3) and 8 captures the regularity of Σ and Ü, denotes the tubular neighborhood of width 8 with a special "recovery" property (see Theorem C.I, Theorem C.2 for the precise definitions of 8 and Üs).
4
Main results
4.1
Construction of OLLA via Least Squares
We now derive the continuous-time dynamics OLLA by choosing the drift vector and the symmetric diffusion matrix Q to be the closest-inn a least-squares sense-to the unconstrained usual Langevin coefficients, subject to enforcing both the equality constraints {hx(z)) and the active inequality constraints {g, (=)))) jez." This is achieved by applying Itô's lemma to each constraint function h; and 9j and splitting the change in, for example, h; into martingale term Vh,((X)TTQXXX (XX,))((VVVV)))) and drift term (X,, ((X,) + (XV)
dt. By
choosing Q so that Vh,))) = Vgr) (x) = 0, the martingale piece vanishes exactly in the normal directions. Simultancously, we pick the drift vector g to satisfy the linear equation
Vh q + Tr (V"1,QQ")) ah; = 0, Vg q + 2
Tr (V QQ²) a(g; +)) = (1)
2
that h (Xx)) ha (n)ee-a and gj (X.) + (gz (Xn) eat where hyperparameters a, e denote the landing or boundary repulsion rate, respectively. This enforces g, (X,) to hit 0 in finite
time = 1T In((g)( (Xo) + t)/e), after which 9; (X,) < 0 remains forever. As a result, this approach removes any noise and drift direction in the normal of constraints and implants a pure drift normal constraints, guaranteeing exponential decay of both equality and active-inequality constraints at rate This yields the closed-fore form SDE in Proposition and Lemma

---

Proposition 1 (Construction of OLLA and its closed form SDE), Consider the following SDE:
dx, = g(x) + Q(X.)dW
(2)
where
QV/, 0, Vi E [m],
Argmin||/211 - Q11
[QVg 0, Vj La,
QERaXa
[V/q i 1r (T((2),,00)) ah; 0, vi E [m],
q argminll- + f VER"
Vg [ 9; + a(g; +c) =0, Vj E Iz.
Then, there exists a closed form (DELO)))) (2) given by:
dx, =
+ "
H(X))dt + V2I1(X))dW, (3)
where
H 7: -==-VJ (T(",) .... hm II. Tr (V 9i, II) Tr (V (V m]]"
9iuzi
(4)
is the associated mean curvature correction term of Mi, = x E Rd h(x) = 0, g1, (z) 0).
Remark 1 (Mean curvature = ttoSStratonovich correction). By technical stochastic-calcut identi- ties on manifolds (see, e.g., Rousset et al. [19], Lemma 3.19), the Itô-Stratonovich correction arising from the Stratonovich SDE
dXx = [TX))V + X
+ V2I(XX)) oW;
(5)
coincides exactly with the mean-curvature term H(x) of E,, = fx E Rd Mxx)) 0, 91, (x) 0}. Remark 2 (Relation to orthogonal direction samplers from variational KL). Zhang al. [17]
introduced an overdamped-Lange Langevin sampler for equality constraints by minimizing the constrained KL divergence via orthogonal-space variational formulation, and Zhang al. [18] also used similar approach to handle single inequality constraint. In the absence of inequality constraints, our
OLLA dynamies coincide with the equality constrained sampler of [17] up to a modified potential
f(x) = fx) + In ded(ct (G), since the mean curvature correction satisfies
H(x) = diviT(x) + [[ In (Idet G(x))
(see Rousset et al. [19], Lemma 3.21). Correspondingly, whereas their framework yields stationary measure proportional to e (1)(()))g (xxx) under the coarea formula, OLLA converges to the Riemannian volume-weighted density e e-f(=) (z) on Σ (see Rousset et al. [19]. Lemma 3.2), In addition to this, the work [18] enforces a single inequality constraint via a purely deterministic normal drift VVVg/l\Vlll without any stochastic component in the tangential direction. OLLA instead projects noise and drift vectors tangentially even during the landing phase, thereby preserving exploration on the evolving manifold E := {z E Rd h(x)= = h(Xu)e and improving mixing.
Lemma 1 (Exponential decay of constraint functions). The dynamics induced by (3) satisfy the following properties almost surely for Vi E [m], Vj Ixo:
hx(X,) = hx(o))e att
t20 0
(7)
and
1
g; +e
9,(X,) = - (X,) + -c)e-at
s in
a
e
1
(x() +
9, (x(X)) 50
ln
a
with (X.)) < 1, 2 ffr j Txo, where I == {kr [1] gk (x: 0} is the index set of active
inequality constraints.
4.2
Non-asymptotic Convergence Analysis of OLLA
In this subsection, we establish non-asymptotic convergence guarantees for OLLA in three scenarios equality only case, inequality only case, mixed case- by recognizing that OLLA has rapid landing

---

property driven by landing rate and natural mixing on the landed manifold induced by the LSI constant. For the detailed notations and proofs, we refer the Appendix A and related Appendices.
Equality-Only Scenario. When the constraints consist of only smooth equalities h((x) = 0, the continuous-time OLLA dynamics (3) can be written
(8) The drift term naturally decomposes into tangential term, which moves along the constraint surface, and a normal landing term, which forces each coordinate h2 (X,) to decay exponentially fast, as summarized in Lemma
dXx [IX + G (X.)((X))) dt + (X)( dW,
Once Xt ~ Pt. lies within a tubular neighborhood of Σ, the nearest projection map ит onto y becomes available, and the projected process Y, πT (X,)) Prr can be defined (Theorem C.1). Then, the regularity of Σ naturally implies that ||X, - Y, 12 [|A(X,)]l2 = O(e=af) holds (Lemma C.I), O of) (Lemma E.1), where Pt, Pt are the laws of Xx. Yt, respectively. Furthermore, the combination of Lipschitzness of π and ||X; - Yella = O(e enable us write the projected process Y, as overdamped Langevin dynamics on E with noisy drift vector and diffusion matrix whose norm is bounded by O(e-at) (Corollary E.1). Therefore, the effect of noisy terms can dominated by the effect of the LSI constant, which leads to the exponentially fast convergence of KL a (pellpx) (Lemma E.3). A rigorous combination of these insights gives the following theorem.
leading the contraction of W2(Pt, (t)
Theorem 1 (Convergence result for equality-constrained OLLA). Suppose assumptions (CI) to (C4) hold. Let X, be the stochastic process following the equality-constriined OLLA and let Pts P be the law of X; and its projection Y, (Xx)) on Σ for 2 tean teur := max { 1 In 6, 1 In C.). respectively. Then, for > 2A1sI for all t 2 tem. it holds that
a
a
Mr e +
2 AesI
W2(Pt, s
KL (pellpe)
K
where
KL (pell))) & exp (22sst(t - bead 211s1O8, (e (e-atte e-alar) ) (pe. (penollss) Orl
-
for some constants C3 0 -1 |1 Ma +
(( (Crghin)
Mhn ()) Or
Cae Cater
O with CLa being
a
a
a a-2
the Lipschitz constant of VT(x)[I(Z) on (E).).
Inequality-One Scenario. With smooth inequalities g(x) < O, we introduced a small boundary repulsion parameter e >> so that each initially violated constraint 9j > driven the boundary within a finite landing time teut
(n((,,(Xo) + e)/e).
From the Pokker-Plannee perspective, the normal probability flux at the boundary satisfies (n, Jo) = -oept where J is the probability current density of pr and is the outward unit normal vector of Σ. The boundary repulsion enforces the outward probability flow to become zero after the landing time feut, and the normal probability flux vanishes to zero after time feur (Lemma F.1). This enables to ignore the boundary integral appearing on the time derivative of KL(pt] (ps) guaranteeing exponential decay driven by the effect of the LSI constant. Therefore, the following theorem comes by directly analyzing the time derivative of (L p5z) after time lent-
Theorem 2 (Convergence result for cecquality-constraiie (C4) hold. Let X, be the stochastic process following the inequality-constriined OLLA and let Pt be the law of X.. Then, for i. benry Lou = a
OLLA), Suppose assumptions (CI) to
1 In
(Mate)
Ma+te
2
Wa(pt,pr) S $ Aesi
KL (pellor))
where
KL <L s 2desiitt- KL (ptalllp)))
5

---

Mixed Scenario. In the mixed setting, OLLA's dy- namics are sensitive the boundary repulsion from g. It is noteworthy that Σ remains unchanged for different inequality functions g, as long as the boundary of the feasible set, where g(z) = 0 intersects h(x)) = 0, is identical. Nevertheless, the choice of function g af- fects the landing dynamics, which in turn can alter the convergence rate of OLLA.
X
h(x)
To quantify this, we define the projected manifold Ep := T x E Rd h(x) = p. (x) (}} and as- sume the norm of boundary velocity op of 82, reg-
g(x) 0
Figure I: OLLA trajectory in mixed case
ulated by Vilpl / for some 8 > 0 (Assumption (M2)). Additionally, we assume Ex, lies inside (L) for 0 < [|l|2 < 5 to avoid stopping behavior of Y, on عK (Assumption (MI)). Under these assumptions, the trajectory and inequality-only can be seamlessly combined. The following theorem arises from a rigorous inte- gration of previous high-level ideas. We also provide theorem with relaxed version of Assumption (MI) in the appendix; see Remark 3 (when Xo 5(x)),, Remark 4, and Corollary A.1 for the details.
X, can be illustrated in Figure I and the proof ideas of equality-only
Theorem 3 (Convergence result for miixed-constrain and (MI) to (M3) hold. Define X, to be the stochastic process following OLIA dynamics (3) and Pt be the law of Y, == x (Xx) after t 2 tears teus == max In
OLLA). Suppose assumptions (CI) (C4)
M, M,+e)
1 In ( 1 In (n(C5)) J.
a
Then, for 2) and 3 2 1, the following non-asymptotic convergence rate of W2 (Pc. Ps.) can be obtained as follows
W2(pt,
Mn e + V Aesi
2
KL KL (pel
1
where
-e--aaa-eaee KL (āelle) S exp -22iss(tt teu) 21esiOs (e -at - -e-amm) 2 (pem
( 1lps) + θr
for some constants G4, Gs, Go. Go > :=: + aGGG5.M)) a 2-
Ö, 0 01-0 CL + (Crn ))
C -aBtr (GoVV
0x. 2\
and with Cr, being the Lipschitz constant of (x)))((()) on Ús(E).
4.3 Euler-Maruyama Discretization & Hutchinson Trace Estimation
To implement OLLA in practice, we discretize the continuous-tin SDE by the Euler- Maruyama (EM) update. At each iteration, we compute three components for the drift vector: the projected gradient drift, the landing drift, and the mean curvature correction H. In particular, H (4) requires forming the full Hessian of each constraint and computing traces of the form Tr (IIVhh;) and (IIV 9;), an Od grad-cost) operation that quickly becomes infeasible in high dimensions. We therefore employ the Hutchinson trace estimator [22], which gives the following approximation:
N (Tx)) ((V2)vue)) Tr ( 93)
N
Tr (IVV2)s)
Σ 1=
R
~
N
N
for each i E [m]. E In where 2/k ~ AV I2) are independent standard normal samples. Each Hessian- Vector Product (HVP) 42 (or V2 V29g0z) can be computed a cost similar to one gradient evaluation, so N probes incur only O(N grad-cost) computational cost rather than O(d grad-cost), saving significant computational cost in high-dimension circumstances. In our experiments, N suffices to achieve low variance estimates that match the performance of the full Hessian computation.
By combining these numerical schemes, we arrive at the full algorithm of OLLA in Algorithm 1.

---

Star
Two Lobes
Quadratic Poly
Mixture Gaussian
OLLA
CGHMC
Equaliy boundary
Inequality boundary Feasible region
Figure 2: Scatter plots of 200 samples from OLLA (top row) and CGHMC (bottom row) on four 2D synthetic examples. Solid lines show equality constraints, dashed lines show inequality boundaries, and green shaded areas mark feasible region by inequality constraints. OLLA closely matches the CGHMC samples in each scenario.
5 Experiments
To demonstrate the sampling accuracy and efficiency, we compare OLLA and its Hutchinson- accelerated variant (OLLA-H) against three standard constrained samplers: CI Clangevin [20]. CHMC [20], and CGHMC [21]. While the three baselines were originally designed for equality constraints, we introduce one slack variable per inequality constraint via 9j (I:) + _ 0 for each j [7] so that (y(x)) < 0 is enforced. CGHMC, meanwhile, does not rely slack variables, but instead uses Metropolis-Hasting correction step not only to reject samples based on the energy the proposed samples, but also based on the inequality constraint violation. Therefore, CGHMC allows more accurate and unbiased sampling from the constrained distribution than CHMC and CLangevin in the long run, and we use it as ground truth for measuring the accuracy of our methods.
5.1 Synthetic 2D Examples
We first evaluate sampling on two-dimensional manifolds: (1) star-shaped equality manifold, (2) two- lobe inequality manifold both with uniform density, (3) manifold defined by poolaaiccppl equality and inequality constraints with a standard Gaussian target, (4) Gaussian mixture with nine components restricted to seven lobes manifold by a nonlinear inequality (both for mixed scenario).
For each problem, we ran 200 independent chains in parallel with 5,000 steps, and collected the last 200 samples. We then computed the W2 distance and energy distance between the empirical well the constraint violation defined
distribution of the samples and the target distribution, by [z][( E max ( respectively. The results are shown in Figure and Figure Further details of example setups and additional results are included in Appendix H.l.
Sampling Accuracy Constraint Violation of OLLA. As shown in Figure 3, OLLA and OLLA-H match the per- formance of CHMC and CLangevin in both Wasserstein and energy-distance metrics, demonstrating that our landing- based approach attains sampling accuracy on par with estab- lished methods. Also, constraint violations for OLLA and OLLA-H remained at low levels without computationally expensive projection steps.
Table 1: Effect of W [[]]]]]]
W3
B[[h(=)]1
0.363
0.682
10
0.130 0..00
100
0.159
0.017:4..01
200
0.12110.00
0.008
Effect of Hyperparameters a and We further examine the influence of the landing rate and boundary repulsion e on sampling accuracy and constraint satisfaction. As shown in Tables and 2, increasing a accelerates convergence, yielding smaller W/2 values and reduced equality constraint
Table Effect of on W3. B[mmax g
W3
E|max g (x)
0.1
0.151 20.06
0.08240.017
1
0.108
0.067
violations, consistent with our theoretical prediction that larger a enhances the landing and contraction rates toward However, excessively large a leads to numerical instability,
5
0.123=0.018
0.040+0.015
10
0.112=00034
0.019+0.006
7

---

OLLA
OLLA-H
CLangevin
CHMC
Energy Distance
W3 (Sinkhom)
Mean [b(o]
0.3
0.03
0.2
0.1
0.2
0.01
0.0
2000
4000
5000
1000
4000
2000 3000
5000
Iteration
Iteration
Iteration
Figure 3: Convergence diagnostics on the Gaussian mixture of 9 components on the 7 lobes manifold
with a 100, distance to CGHMC samples, and (3) mean constraint violation [[]]]- Solid lines and shaded bands show the mean +1 SD over five independent runs. Both OLLA and OLLA-H rapidly decrease TR[ [h] down small values and maintain it there, which achieving the lowest energy and W errors.
From left to right: (1) energy distance CGHMC samples, (2) squared W?
OLLA
OLLA-H
Clargevin
CHMC
CGHMC
CPU time/ESS
Pixi > 0)
K(X) (test function)
0.7
0.6
0.4
0.4
0.2
0.0
0.1
0.0
1.00
200
100
400
700
100
200
300
400
500
(dimension)
(dimension)
d (dimension)
Figure 4: Sampling performance and accuracy as the dimension d increases (with m = left to right: (1) CPU time per ESS versus d, (2) Estimated probability P(x1 0) versus d, (3) Estimated value of ((x) versus d. Shaded bands shows 1SD over five runs.
= 5). From
causing W to rise and the sampler to collapse (Table 7). A similar trend is observed for Stronger repulsion lowers inequality violations, but beyond certain range, W remains nearly unchanged, aligning with the continuous-timg time theory that e primarily affects the finite landing time rather than the asymptotic convergence rate. Overly large values can again destabilize the dynamics and lead to numerical breakdown (Table 8).
5.2 Scaling of OLLA under High- Dimenniimessonnalty and Large Number of Constraints
We assess the robustness and scalability of OLLA and OLLA-H using a synthetic stress-test problem that enables explicit control over the ambient dimension d and the numbers of equality and inequality constraints (m, l). Samples are drawn from a uniform distribution on a constrained manifold defined by linear and quadratic constraints. Each algorithm runs for 1,000 iterations with burn-in and thinning, and we vary one of d, m, or while fixing the others to disentangle their individual effects.
We evaluate two key metrics: (1) CPU time per effective sample (CPU/ESS), and (2) the accuracy of representative test function estimates such as P( 0) and K(x) sin(21)e"² + log([æs/
1) tznh(z4) + II cos(x.). Detailed experimental setups are provided Appendix H.2.
Scaling under Dimension. Figure 4 illustrates the sampling performance of algorithms as 5). On the left, CPU time/ESS OLLA-H remains essentially flat around 0. 0.05s/sample while OLLA grows linearly (reaching = 1. 1.1s/sample at d = 700), and CHMC and CLangevin stay at 0.2 0. 3.3s/sample and 0. 0.1s/sample respectively. In the center, both OLLA-H, CHMC, and CGHMC maintain P(z1 > 0) R 0.5, whereas CLangevin collapses to 0.2 in high dimensions, indicating severe bias. the right, the estimate of nonlinear test function K (z) shows that OLLA, OLLA-H, CHMC, and CGHMC all produce virtually identical estimates
increases from 10 to 700 (with m

---

OLLA
OLLA-H
Clargevin
CHMC
CGHMC
CPU time/ESS
p(x 0)
K(X) (test function)
0.9
0.7
0.8
0.6
0.6
0.5
0.4
0.3
2
0.2
0
0.1
0.0
10
20
30
60
30
10 20
30
40
30
10
20
30
40
50
60
finequalitt count)
(inequality count)
(inequaliit count)
Figure 5: Sampling performance and accuracy as the number of inequality constraints / increases (with d 100, m = 5). From left to right: (1) CPU time per ESS versus I, (2) Estimated probability P(x1 > 0) versus I, (3) Estimated value of K (æ) versus l. Shaded bands shows + ISD over five runs. Note that OLLA and CGHMC results on CPU time/ESS overlap almost perfectly, suggesting their comparable performance on this metric.
even as d grows, while CLangevin lags behind. Overall, the results indicate that OLLA-H scales, maintaining reliable performance the dimension d' increases.
Scaling under the Number of Constraints. With d 100, we separately increased the number of equalities m (with = 5) and the number of inequalities (with m = 5). In the presence of large number of equalities, OLLA and OLLA-H may lose their edge over equality constrained specialized baselines. In these situations, the maximum that is stable in the discrete algorithm becomes limited to prevent significant discretization error. However, we observe that even, at near limit value of a, [[[h]] may remain relatively high unless the step size A/ is further reduced. We add an additional study for this on Appendix H.2. By contrast, when adding more inequalities, OLLA-H stays relatively accurate and shows very low CPU time/ESS (Figure 5).
5.3 Molecular System with Realistic Potential
We further evaluate samplers on a molecular system with realistic potentials and geometric constraints. The system incorporates equality constraints h (fixed bond lengths and angles) and inequality constraints (steric hindrance), alongside torsion and Teeck--Chandle-AAndessoo potential terms. Experiments were conducted for increasing dimensions d 3N by varying the number of atoms. Detailed experimental details are provided in Appendix H.3 and Table 3 summarizes the results.
Table Estimates of radius of gyration squared (R3) and total CPU time for sampling (in brackets) on the molecular system with realistic potential (complex Σ, small [T,], large Nnewtn). The average samplers
constraint violation for OLLA-H was below 0.007 (equality), while projection-basee maintained violations below 0.0001. Inequality violations were observed to be O for all algorithms. The (d, 1) configurations are: (15, 7, 6). (30, 17, 36), (45, 27, 91), (60, 37, 171), (90, 57, 406).
method/ / dim (d)
15
30
45
60
90
1.392 + 0.026
5.414 = 0.017
12.240 + 0.102
21.820 + 0.098
49.080 + 0.223
OLLA-H (N = 0)
(78s)
(151s)
(254s)
(334s)
(704s)
OLLA-H (N = 5)
1.370 + 0.032
5.424 = 0.015
12.200 + 0.153
21.780
+ 0.098
48,940 + 0.273
(182s)
(400s)
(656s)
(916s)
(1732s)
CLangevin
1.396 i 0.012
5.526 1 0.015
12.400 1: 0.000
22.140
1: 0.049
49.960 & 0.049
(421s)
(3620s)
(10940ss)
(22600s)
(52220s)
1.410 + 0.000
5.580 + 0.000
12.500 + 0.000
22.200
+: 0.000
49.960 + 0.049
CHMC
(1478)
(468s)
(1012s)
(1712s)
(3660s)
1.410 + 0.000
5.580 + 0.000
12.500 + 0.000
22.200 + 0.000
49.940 + 0.049
CGHMC
(135s)
(282s)
(467s)
(652s)
(1116s)
As dimension d grows, the feasible set Σ becomes increasingly complex, while the number of active inequality constraints [r,] remains small. However, the numerous equality constraints make each projection step demanding for projection-based methods (CHMC, CLangevin, CGHMC), resulting in substantial Newton iteration Naewinn and increased computational costs. Particularly, CLangevin significantly affected from this, whereas CGHMC is faster but still slower than OLLA-H (N 0).

---

In contrast, increasing the Hutchinson probe N from N = 0 to N = 5 improves the mean equality constraint violation across all dimensions d, reducing it from a range of 0.0055 down to 0.00035.
Despite this, both OLLA-H configurations (N = 0, N = 5) yield comparable test function estimation accuracy. This suggests that the constraint violation at N = 0 was already sufficiently low, rendering the improvemen''s & impact on the final sampling accuracy negligible. Consequently, the N 0 wahie HVP evaluations are completely skipped-achiewe achieves significantly lower computational cost. This behavior aligns with findings in Zhang et al. [17], suggesting that omitting the Itô- Stratonovich (mean curvature) correction term has negligible practical impact on sampling accuracy, while markedly improving efficiency.
5.4 Bayesiaa Logistic Regression with Fairness and Monotonicity Constraints
We evaluate the samplers high-dimensional Bayesian logistic regression task using two-layer neural network trained on the German Credit dataset [41]. The setup enforces fairesss through equality constraints h ensuring parity in true positive rate and false positive rate between demographic
groups, along with monotonicity constraints g on selected input data. Further details of experimental setup are provided in Appendix H.3.
This problem poses significant challenges for projection baased samplers. Step sizes effective in unconstrained scenarios led to projection failures for CI Clangevin and CHMC, and to acceptance rates below 5% for CGHMC. To maintain stability, their step sizes were substantially reduced, In contrast, the projection-free OLLA-H remained stable across all settings and, for fairness, was also evaluated with the same reduced step size as CLangevin.
As summarized in Table 4, OLLA-H (N = 0) consistently attains the lowest test negative log- likelihood (NLL) while being orders of magnitude faster than projection-based baselines. Although the projection- based methods achieve tighter feasibility. OLLA-H maintains small violations without noticeable degradation in accuracy. These results highlight OLLA-H's robustness and computational advantage in nigh-dimensional constrained Bayesian logistic regression task.
Table 4: Test NLL and total CPU time for sampling (in brackets) on the Bayesian logistic regression with fairness and monotonicity constraints (tigh-dimensional 2.). The average constraint violations for OLLA-H were below 0.005 (equality) and 0.15 (inequality). while projection-base samplers (CLangevin, CHMC, CGHMC) maintained feasibility below 0.0008 (equality) and no inequality violation.
method/ / dim (d)
1986
4994
9986
49986
100002
OLLA-H (N 0)
0.514 0.013
0.521 + 0.008
0.524 + 0.014
0.523
0.011
0.520 0.015
(63s)
(70s)
(70s)
(81s)
(82s)
OLLA-H (N = 5)
0.520 = 0.013
0.524 + 0.008
0.505 + 0.004
0.517 = 0.011
0.516 + 0.004
(159s)
(180s)
(205s)
(189s)
(197s)
CLangevin
0.573 = 0.004
0.568 + 0.013
0.564 + 0.022
0.580 + 0.005
0.570 + 0.011
(1162s)
(1176s)
(1194s)
(1428s)
(1370s)
CHMC
0.599 = 0.015
0.595 + 0.020
0.599 + 0.017
0.606 = 0.004
0.605 0.001
(526s)
(532s)
(561s)
(586s)
(611s)
CGHMC
0.600 = 0.007 (76s)
0.600 + 0.009
0.606 + 0.003
0.598 = 0.020
0.601 + 0.007
(77s)
(82s)
(83s)
(88s)
6 Conclusion & Future works
We have presented Overdamped Langevin with Landing (OLLA), a projection-free SDE sampler that enforces nonlinear equality and inequality constraints by deterministically "landing" trajectories onto the feasible set while retaining full tangential noise, and proved that its continuous dynamics converge exponentially fast in 2-Wasserstein distance under appropriate regularity. Building on this, we proposed OLLA-H, an EM discretization that uses Hutchinson trace estimator for approximating the tt--Stratonovic correction only O(N grad-cost) per step, and showed in both 2D and high- dimensional tests that it matches the accuracy of established constrained samplers while drastically reducing runtime. Future work will include non-asymptotic convergence guarantees for the discrete algorithm-closin the gap between SDE theory and implementation and developing OLLA variants that remain stable even with many equality constraints in very high dimensions, further extending its scope to large-scale constrained probabilistic inference.
10

---

Acknowledgements
KJ and MT are grateful for partial supports by NSF Grants DMS- 1847824 DMS-2513699, DOE Grants NA0004261, SC0026274, and Richard Duke Fellowship. MM is supported by the German Research Foundation. The authors thank Andre Wibisono for insightful discus- sion.
References
[1] Arnak Dalalyan. Theoretical guarantees for approximate sampling from smooth and log- concave densities. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79 (3):651-676, 2017. 1
[2] Alain Durmus and Éric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. The Annals of Applied Probahility, 275)):551- 2017.
[3] Alain Durmus and Éric Moulines. High-dimensional Bayesian inference via the unadjusted Langevin algorithm. Bernoulli, 25(4A):2854-2822 2019.
[4] Murat Erdogdu, Rasa Hosseinzadch, and Shunshi Zhang. Convergence of Langevin Monte Carlo in chi-squared and Rényi divergence. In International Conference on Artificial Intelligence and Statistics, pages 8151-8175. PMLR, 2022.
[5] Xiang Cheng and Peter Bartlett. Convergence of Langevin MCMC in KL-divergence. In International Conference on Algorithmic Learning Theory, pages 186-211. PMLR, 2018.
[6] Arnak Dalalyan. Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent. In Conference on Learning Theory, pages 678-689, PMLR, 2017.
[7] Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis Langevin Monte Carlo via convex optimization. Journal of Machine Learning Research, 20(73):1-46, 2019.
[8] Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt (d) dimension dependence of Langevin Monte Carlo. In International Conference Learning Representations, 2022.
[9] Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted Langevin algo- rithm: Isoperimetry suffices. In Advances Neural Information Processing Systems, volume 32, 2019.
[10] Sinho Chewi. Log-concave sampling. Technical report, 2024. URL MtRLhttts///chewiiinhe github. ww/mmii..i........ htt.//aai..........mmmmmmmm
[11] Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror Langevin Monte Carlo. In Conference on Learning Theory, pages 3814-3841. PMLR, 2020. 1,
[12] Ruilin Li, Molei Tao, Santosh S Vempala, and Andre Wibisono. The mirror Langevin algorithm converges with vanishing bias. In International Conference on Algorithmic Learning Theory, pages 718-742 PMLR, 2022. 1,2
[13] Maryam Zare, Parham M Kebria, Abbas Khosravi. and Saeid Nahavandi. survey of imitation learning: Algorithms, recent developments, and challenges. IEEE Transactions on Cybernetics, 2024.7
[14] Michael Muehlebach and Michael I. Jordan. On constraints in first-order optimization: A view from non-smooth dynamical systems. Journal of Machine Learning Research, 23(256):1-47, 2022. 2
[15] Michael Muehlebach and Michael Jordan. Accelerated first-order optimization under nonlinear constraints. Mathematical Programming, pages 1-46, 2025. 2

---

[16] Sholom Schechtman, Daniil Tiapkin. Michael Muehlebach, and Eric Moulines. Orthogonal directions constrained gradient method: from non-linear equality constraints Stiefel manifold. In Conference on Learning Theory, pages 1228-1258. PMLR, 2023. 2
[17] Ruqi Zhang, Qiang Liu, and Xin Tong. Sampling constrained domains with orthogonal-space variational gradient descent. In Advances in Neural Information Processing Systems, volume 35, pages 37108-37 37120, 2022. 2, 4, 10
[18] Shiyue Zhang. Longlin Yu, Ziheng Cheng, and Cheng Zhang. Functional gradient flows for constrained sampling. In Advances in Neural Information Processing Systems, volume 37, pages 112356-1 12385, 2024. 4,4,54 54
[19] Mathias Rousset, Gabriel Stoltz, and Tony Lelievre. Free energy computations: a mathematical perspective. World Scientific, 2010. 2, 4, 30, 31, 32, 38, 46
[20] Tony Lelievre, Mathias Rousset, and Gabriel Stoltz. Langevin dynamics with constraints and computation of free energy differences. Mathematics of Computation, 811220):2011-2125, 2012. 2, 7, 30, 31
[21] Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz. Hybrid Monte Carlo methods for sampling probability measures on submanifolds. Numerische Mathematik, 143(2):379-221, 2019 2, 32
[22] Michael Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in 1076, 1989. 2, 6
and Computation, 18(3):1059-
[23] Thomas Power and Dmitry Berenson. Constrained stein variational trajectory optimization.
IEEE Transactions on Robotics, 40:3602-3619, 20204.44..00........ 2
[24] Giovanni Ciccotti and Jean-Paul Ryckaert. Molecular dynamics simulation of rigid molecules. Computer Physics Reports, 466):46-399) 1986.
[25] Jean-Paul Ryckaert, Giovanni Ciccotti, and Herman JC Berendsen. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. Journal of Computational Physics, 23(3):327-341, 1977.
[26] Hans C Andersen. Rattle: A "velocity" version of the shake algorithm for molecular dynamics calculations. Journal of computational Physics, 52(1):24-34, 1983. 2
[27] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from a log-concave distribution with projected Langevin Monte Carlo. Discrete Computational Geometry, 59:757-783, 2018.
2
[28] Sebastien Bubeck, Ronen Eldan, and Joseph Lehec. Finite-time analysis of projected Langevin Monte Carlo. In Advances in Neural Information Processing Systems, volume 29, pages 1253-1251, 2015.
[29] Yunbum Kook, Yin Tat Lee, Ruoqi Shen, and Santosh S Vempala. Sampling with Riemannian Hamiltonian Monte Carlo in constrained space. In Advances Neural Information Processing Systems, volume 36, pages
2022. 2
[30] Yunbum Kook and Santosh S Vempala. Gaussian cooling and Dikin walks: The interior point method for logconcave sampling. In Conference on Learning Theory. PMLR, 2024. 2
[31] Kanji Sato, Akiko Takeda, Reiichiro Kawai, and Taiji Suzuki. Convergence error analysis of reflected gradient Langevin dynamics for non-convex constrained optimization. Japan Journal of Industrial and Applied Mathematics, 42(1):127-151, 2025. 2
[32] K Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-Langevin algorithm. Advances in Neural Information Processing Systems, volume 34, pages 28405- 28418. 2021. 2
12

---

[33] Avetik Karagulyan and Armak Dalalyan. Penalized Langevin dynamics with vanishing penalty for smooth and log-concave targets. In Advances in Neural Information Processing Systems, volume 33, pages 17594-17604, 2020. 2
[34] Mert Gurbuzbalaban, Yuanhan Hu, and Lingjiong Zhu. Penalized overdamped and underdamped Langevin Monte Carlo algorithms for constrained sampling. Journal of Machine Learning Research, 25(263):1 1-67, 2024. 2
[35] Luiz F Chamon, Mohammad R Karimi, and Anna Korba. Constrained sampling with primal- dual Langevin Monte Carlo. In Advances in Neural Information Processing Systems, volume 37, pages
2024. 2
[36] Mark Girolami and Ben Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73
(2):123-214, 2011. 2
[37] Xiang Cheng, Jingzhao Zhang, and Suvrit Sra. Efficient sampling on Riemannian manifolds via Langevin MCMC. In Advances in Neural Information Processing Systems, volume 35, pages 5995 6006, 2022.
[38] Khashayar Gatmiry and Santosh S Vempala. Convergence of the Riemannian Langevin algo- rithm. arXiv preprint arXiv:2204.10815 2022.
[39] Lingkai Kong and Molei Tao. Convergence of kinetic Langevin Monte Carlo on Lie groups. In Conference on Learning Theory. PMLR, 2024.
[40] R Tyrrell Rockafellar and Roger J-B Wets. Variational analysis, volume 317. Springer Science Business Media, 2009. 3
[41] Hans Hofmann. Statlog (German Credit Data). UCI Machine Learning Repository, 1994. DOI: https//ddb.og/1224322SCCCC 10, 61
[42] Leonard Gross. Logarithmic sobolev inequalities and contractivity properties of semigroups. In Dirichlet Forms: Lectures given at the Ist Session of the Centro Internazionale Matematico Estivo (CIME) held in Varenna, Ital,, June 8-19, 1992, pages 54-88 Springer, 2006. 24
[43] Felix Otto and Cédric Villani. Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality. Journal of Funetional Analysis, 173(2)3311-000, 2000. 46
[44] OS Rothaus. Diffusion on compact Riemannian manifolds and logarithmic Sobolev inequalities. Journal of Functional Analysis, 42(1):102-109, 1981.
[45] OS Rothaus. Hypercontraccivvtty and the bakry-emery criterion for compact Lie groups. Journal of Functional Analysis, 65(3)3358-367, 1986.
[46] Feng- Yu Wang. On estimation of the logarithmic Sobolev constant and gradient estimates of heat semigroups. Probability Theory and Related Fields, 108:87-101, 1997.
[47] Xiao Wang, Qi Lei, and loannis Panageas. Fast convergence of Langevin dynamics manifold: Geodesics meet log-Sobolev. In Advances in Neural Information Processing Systems, volume 33, pages 11894-189044 2020. 24, 29
[48] Laurent Saloff-Coste. Uniformly elliptic operators on Riemannian manifolds. Journal of Differential Geometry, 36(2):417-450, 1992. 24
[49] Michel Ledoux. Concentration of measure and logarithmic Sobolev inequalities. In Seminaire de probabilites XXXIII, pages 120-216. Springer, 2006. 29
[50] John M Lee. Introduction to Riemannian manifolds, volume 2. Springer, 2018. 33
[51] Gregory Chirikjian. Stochastic Models, Information Theory, and Lie Groups, Volume Classical Results and Geometric Methods. Springer Science & Business Media, 2009. 41
13

---

[52] Chin- Wei Huang, Milad Aghajohari, Joey Bose, Prakash Panangaden, and Aaron C Courville. Riemannian diffusion models. In Advances in Neural Information Processing Systems, vol- ume 35, pages 2750-2761, 2022. 41
[53] Sergej G Bobkov and Friedrich Götze. Exponential integrability and transportation cost related to logarithmic Sobolev inequalities. Journal of Functional Analysis, 163(1):1-28, 1999. 46
[54] A Paszke. Pytorch: An imperative style, high- peefoperforance deep learning library. In Advances in Neural Information Processing Systems, volume 32, 2019. 54
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:
Summarized contributions in Section 1.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in 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 perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how be expected to generalize other settings.
much the results
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Section 5.2 and Section illustrates a limitation.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in 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 the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might 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 address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
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

---

Justification:
Appendix A, C to G and Section 2, 4 address theoretical results.
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 referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide 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 be properly referenced.
4. Experimental result reproducibility
Question: 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 whether the code and data are provided or not)?
Answer: [Yes]
Justification: Section and Appendix H covers experiment setting reproduce results.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a 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 a 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 clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either 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 to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers 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:
Code is provided. Running ipynb files can reproduce the results.
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, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: /nips. wwtppbccpppppaecccccccccccccccee / ssisssolicyyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
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 supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: Section anb Appendix illustrates the experimental setting and details. Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All the plots and tables have + SD and mean values and trial numbers specified.
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, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are 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 еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Appendix H illustrates the computing resources for experiments.
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, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: httprrpp.....ublllcchccccttttedd
Answer: [Yes]
Justification: The research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: /ettps://nne http:://eeuuusseeeee http:///.........etttccccccceee.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- cration due to laws or regulations 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: The paper covers theoretical analysis and introduces an algorithm, which does not relate to societal impacts.
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:
The paper does not pose 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 from another parties.
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 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: The paper does not introduce new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing 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: The paper 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 in 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 the NeurIPS Code of 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:
The paper 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 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 expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
20

---

Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: LLM not involved during the development of the research.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs any important, original, non-standard components.
Please refer to our LLM policy (https: /tap:::///eurpssss http///rpps.cccooffeerrreeee//// ccc///ooferrreccee222///11/ for what should or should not be described.
21

---

TABLE OF CONTENTS
1
Introduction
1
2 Related Works
2
Preliminaries Notations
3
4
Main results
3
4.1
Construction of OLLA via Least Squares
3
4.2
Non-asymptotic Convergence Analysis of OLLA
4
4.3
Mally-Maruyyma Discretization & Hutchinson Trace Estimation
6
5
Experiments
7
5.1
Synthetic 2D Examples
7
5,2
Scaling of OLLA under High-Dimensionality and Large Number Con- straints
8
5.3
Molecular System with Realistic Potential
9
5.4
Bayesian Logistic Regression with Fairness and Monotonicity Constraints
10
6 Conclusion & Future works
10
A Table of Key Notations, Assumptions, and Remarks
23
B Algorithms OLLA, OLLA-H, and baselines
30
C Proof of Basic Properties Σ
33
C.1 Intrinsic Gradient on Σ
33
C.2 Intrinsic Divergence on Σ
33
C.3 Recoverable Tubular Neighborhood of Boundaryless Riemannian Manifold
33
C.4 Recoverable Tubular Neighborhood of Riemannian Manifold with Boundary
35
D Construction of SDE with Exponentially Fast Decaying Constraints
37
E Proof of Theoretical Results Equality-constrainn OLLA
39
E.1 Upper Bound of (2(Pt,
39
E.2 Upper Bound of KL 2 (Prllps)
40
F Proof of Theoretical Results Inequality OLLA
47
F.1 Convergence Result for Inequality-constrained OLLA
47
G Proof of Theoretical Results Mixeddconstrainee OLLA
48
G.1 Upper Bound of (2(pt, Pt)
48
G.2 Upper Bound of KL Σ (pellpz)
49
H Experiment Settings and Supplementary Results
54
H.1
Experiment Settings and Supplementary Results for Synthetic 2D Examples H.2 Experiment Settings and Supplementary Results for High-dimensional
54
Manifold with Large Number Constraints
58
H.3 Experiment Settings for Molecular system and Bayesian logistic regression
task
60
22

---

A Table of Key Notations, Assumptions, and Remarks
Table 5: Table of Key Notations
Symbol
Definition
Descriptions
h
h(x) (z)(
Am (x)]²
Equality constraints
g
g(x) = (z). (x))]
Inequality constraints
Σ
(a E I h(x) = 0, ((x) 50}
Constraint manifold
I
(i [1] (x) 2 0] (11, }
Active index set of inequalities
91.
(x) 19, (z). (I))
Active inequality constraints Constraint-correcte vector Orthogonal projector onto T
((x)
{m(z)" (x) c.... ii (x) + cc
((x))
T
( Vh(z)u 0, Vgr, (z)u 0}
Tangent space of Σ at Intrinsic gradient on Σ (C) Intrinsic divergence on Σ (C)
Vxd
[())V/(>)
div=X
Tr (XI)))VVX())
doy
Induced surface (Hausdorff) measure of D
Surface measure on 2
G(w)
V(x=)VJ(x>²
Gram matrix, full rank assumed Usual tubular neighborhood See details in (C.1, C.2)
U,(E)
Tubular neighborhood of Σ with reach
(E)
Recoverable tubular neighborhood with width
M₁
Eupppe))))( [1/xx))]
Initial bound of h(x:)
Mg
suppoe) (
Initial bound of g(x)
xxx)
arg minyez/l2 v/ا
Nearest point projection onto Σ
Boundary repulsion rate
Controls effect of repulsion.
a
Landing rate
Controls constraint decay
A
Nog-Sobolev constant of ps: on E Density of the process X, at time Density of the projected process Y, at time
Enable KL Σ < 1 1
2
Pt
Law of X, (following OLLA)
Pt
Law of Y; == (X,)
P
Target (stationary) density on Σ
Proportional to exp( - -f)dax
KL (pl(m)
fy Edoy L fo Vll In f1idd
KL.-divergence on Fisher information on E
(ell
R
Regularity constant of Ep:= (fx Rd h(x) p. (x)) s $0)
See details in (C.1, C.3)
Ep
Projected manifold for ||l|| S
Boundary velocity of OΣ,
See details in Theorem
23

---

Assumption (Assumptions for exponential convergence of OLLA), The followings are assump- tions used for proving the convergence results of OLLA:
(CI) LICQ.
The linear independence constraint qualification holds on E if, for every E E,
[V), U (V)
set of linearly independent vectors.
(C2) Manifold regularity.
The Riemannian manifold Σ is compact and connected. Also, the equality constraint function h(x) is coercive, i.e., [|<(2)]]2 - oo as ||>||e - oc. Assume Vh(z) 0, Vz E E and dim(E) d when there only inequality constraints.
(C3) Initial constraints bounds.
The initial distribution po satisfies
MA == sup Exupp(po)
[lb(20)l22 < co, Mg := sup l1g(2o)lla < oc.
Espp((p)
(C4) Log-Soboley constant.
The stationary distribution p5(x) x exp(-f(z))do3 satisfies a Log-Sobolev Incquality (LSI) with constant ÄLSI so that
1 (pllps).
KL (plpe)) S 2ALSI
Note: Compactness of L with non-negativity of Ricci curvature guarantees the LSI condition. See [42-47] Remark for detailed description.
(MI) Regularity of Ep. lies inside the interior of E for 0 < llpl|2 5, where 8 is the width of recoverable tubular neighborhood (2))..
The projected manifold Σ, M({ RC h(x) = p. gx)) 0})
(M2) Regularity of ØEp.
The boundary velocity op of 0Σ, appearing on Leibniz integral rule satisfies supseos, | l Vipll for some V 0, Also, assume Ms;
plaas axx, (0E,) 00.
(M3) Bound on Pt. PV- There exists the constants G1, G2, G3 > 0 such that Gi = suPr20,reE Pr < and G2 S ps G3 for every x Σ.
Remark 3 (Comments the assumptions (M1) and (M3)).
Although the assumption (MI) stated for all small perturbation p E I"' ||pl|2 in our anaylsis on Theorem 3, we only require:
Σ = ( E Rd h(x))= = h(Xo)e .(()) S {0} int(2))
for t 2 feuts feui := max In Mate) <
a In (M))) In(CGs) Hence, when Xo ~ 8(xo) for
some zo E IRd we replace the global requirement "for all p with 0 < ||p||2 < 8" by weaker, one-dimensional condition:
(M1') If Xo 5(10) for same xo E Rdd, let M(xg) [ and assume
Eya = = E Rdd (x) su, g(x) 0}) int(E), for all (0,
For the uniform boundedness sup Pe < oo in the assumption (M3), we recall from the proof of Lemma G.2 that satisfies the following Fokker-Planck equation:
0,p: -dvv (V ln bN)) + divs:divv: + 6z)) + &x)).
Since ||62||2 2| = O(e-at, the second-order differential operator of the Fokker-Planck equation becomes uniformly elliptic for sufficiently large In the absense of boundary conditions (equality- only case), the result from Saloff-Coste [48] would then guarantee supr20.zex P: However,

---

the non-standard boundary conditions imposed on Σ in our setting preclude a direct application of those results. We therefore retain the uniform boundedness of Pr as an explicit assumption.
Remark 4 (Relaxed assumption of (M1)). We remark that the assumption (M1) can be strong in practice. To mitigate this issue, we propose the following assumption (M1") which is a milder assumption than (M1). It assumes
(MI") Let 4 == max {a In (o ) In (1 } and define pt == (Xt) 0.) for 2 4 Suppose there exist 7p, Cp 0, 2 L such that pt. O.p. S pp & 1/2 for 2 lp-
c
We note that this assumption is a necessary condition of the previous assumption (MI), and the landing property guarantees lim,-+ Pr = 0. Also, this assumption ensures that both Pt. dept decays exponentially fast, which can be satisfied depending on how fast the sticking behavior of π (Xxe) at the bounadry is attenuated as t increases.
Corollary A.1 (Convergence result for mixed-constrained OLLA with milder assumption). Sup- pose assumptions (C1) to (C4) and (MI") to (M3) hold. Define X, to be the stochastic pro- cess following OLLA dynamics
(3) and Pt. to be the law of == n (X+) after t ' tens tee == Mu. 2 (n(((()) tt)) }. Then, the following non-asymptotic convergence max 1 In M,tc In
rate of (2 PΣ) holds for a > 2dLSI. 3 2 1, and > temi
Mn at + V AisI We(pe,pE) S
2 KL 1 (ps)) + Ope diamy (2))2
where
KL (p (1z)) sexp - - -2 - tem)) 2desiOs
(e (e-at- - -mur)- Cp (e (c-Tpt - - -c-u)
7p
1 O57" + 1
for some constants G Go
GT ,. > 0, OF = (Co + aG Gb,Ma)
ofar
a- 2A1st'
Mn ))
Css :o-0 1+ Cr Mn +
(Cr,,
85 = [G VM 6 aaSter
aB 2d1si*
and with CLa being the Lipschitz constant of (x) on Úo(E), p. being the conditional laws of Yi given Y, inr(E).
Proof. We slightly adjust the argument in the proof of Theorem 3 to show this statement, and the definitions of constants are shared with Theorem 3. We first observe that the probability measure p. can be decomposed as follows:
Pe = (1-pi) + ppptt
22 teut
where p: :- (.E) and
p:= pr bg
1 1- Pi
Pt
are the conditional laws of ýr restricted on Σ and OΣ, respectively.
Because #t is the convex combination of p and p the convexity of the 2-Wasserstein distance implies
(' < (1 W 25) + P. .2)?
Here, we note that
( PP)) <
2
KL E ( l|ps)
25

---

by Lemma E.6 and the fact that p. has its support on in) Also, we note that
DePt _D
1 Pt
(O.p. (O(pt)pt
because Pt PE
holds for E int(E). Then, the same approach used in Lemma G.2 gives
[sps) "p,"
0e
Pt
In
Pt
dos
Pz
1
1
dipi In
doy
+
Pt
In
(ut
ne)dorss.
1 Pt E,
15
Ps
Tern (1)
Term
O.P: KL
+
Pt
+
0e
doy
De
Term (3)
where the last equality holds from the Leibniz integral rule with of being the velocity vector of the boundary of E. (f2 Rd h(x) h(Xo)e .g(x)) 0}). Therefore, the expression of a, p. implies that Term (1) becomes
Term (1) = /
2 (V In bbn - Σ
divs (FF" (Sa + oz)) + -xx)
55z In
Pr
doy
Ex
Ps
k=1
Te-m(1-1)
So In bx)
div (fx) - 8&)
In
Pi
ne
do
ps.
Term Te-2
using (1 P.P
Therefore, Lemma G.2 implies Term (1) can be bounded by
Term (1) S (11 (" Ilss)) + Öoe at + aG G5
6 C2
6:00
} and GF" sup SUP of" S 2G and the
with G G3My max
In
C2
Term (2) can be bounded by
Term (2)l S G VM a.8t
with G. == G G My. Also, Term (3) upper bounded by
Term (3) S Cpe 'KL 2 (p: [pps)
Therefore, combining the bounds with the LSI condition gives the following inequality:
0 (p) ( (1 - CGe att_ Cp
KL
"0t"""" [|ps) + (Co+ aG Gr.Mi) at
2
GG VM
Hence, applying the Grönwall-type inequality recovers the following inequality:
-2 tem) 2desiO's (c- - ea-aas)
Cg
KL x "0;"" |(xx) S exp
(e- 7pt e- -Toleae)
7p
271si (c-as e -alon)
Cp
IKL (2F (pz) +
exp
21LSIls Mhss(ss-oon) tem)
(e (c-2pe- c- T-pplou))
beut
a
Pp
[CGs + aG Gp,M) e aas + G VMMe--aan ds
which can again be summarized as follows:
26

---

2ALSIC' (e-at (e-a-ee -e e-abu) Cp (e- e @-oooou)
a KL (pF (ps)) eexp -2ALSi(t tewe)
7p
KKL IKL (pp)) + 07 + c
-at where C7 (Cs + aG² GyMb)
-- and C5 GG MM)
2A1si
2Aesi
Lastly, observe that
1/4
/ ds(p. xdpp dg)
W pp)) =
PE)
inf (TEll(P! .F) ExE
- diamy (2)
where diamy: (Σ) == sup {dy(z, y) 1, J < 0o due to the compactness of Σ. Therefore, combining the above upper bounds, we recover the following inequality:
( S
KL (
Ope diamg (E)²)
ALSI
Therefore, applying the same reasoning as in Theorem 3 completes the proof.
Remark 5 (Effect of number of Hutchinson probes N the decaying speed of constraint functions). In this remark, under the single equality case, we demonstrate that exponential decay of h via OLLA-H and analyze the effects of the number of Hutchinson probes N. When there is only single equality constraint, we recall that the OLLA-H update rule is given as
Xa+1 = X (6(X2) + (x A + V2AITT(X)) & ~ N Ia)
where b(x) == IT 1 (T))/((()))) - VVh Vh(r) (T<))) and
$(x) = N
T
( (xxx) = VM(r)G(Z))+ (5(x) ( (5(x)))
kit
with S(x) = [1 and UKSs ~ NNO, Io). Also, we note that E [Ex (Xx)) Xx] = and
(Xx Xx] [V((X)]]* (ar(S(())) = VXxa]] (I
Proposition A.1 (Exponential decay of constraint function under OLLA-H), Assume the single equality construint scenario and the following conditions:
(Dissipativity of ((x))
(x,)())) < =m] + for some > 0, 2 0.
(Linear growth of b(x)) [bt))| + B|| for some A, B 2 0.
(Boundedness of tangential fressian-gradie ratio)
C₁ = supneRd
< co.
(Ly-smoothness of
For some L1 0. ||VPh(z)]| < L1 for any I E Rd.
(La-Lipschitzness of 4²h) For some L2 > [[) v s La|la vll for any
Rd.
(Step size control)
< At min {1, 7
Then, the OLLA-H dynamics with Hutchinson probes have the following decaying property of h:
E [ < (1 - adt)* [h( + ( At
d 11) ))
N.)
for K 0.
27

---

Proof. We first prove that supk20 || < oo where AX2 == Xk+1 - Xg. To show this, we - j, - N) to be the & th canonical filtration and observe that m
define Fs == (XXms Em, Um..j
E
| 2 Fic = | 2 (X+)) + ( + ArE [X<s( X&] + 2ATE
[[[ 1 Xx
because [& Xx]
E (Xx) Xx]
0 and AX (b(X²) + cx + (2A)11(Xx(Rx
Applying the disspativity and linear growth assumption, we have
20h
E [||X*+1l1 Fu] - (1 (2m BDt) || + 2c + 2(d (A+
At
N
because
[III((X)) V
20hh
E [ele((Xx)ll) X S
[Xx)]
N
holds from the tangential Hessian Hessian-gradient ratio assumption, and E [Ex Xx] d Because the step size control assumption guarantees following inequality:
(2m BāI) E (0, 1), we have the
2c + 2(d - 1) + (A+ + 2Ch )a sup [| 420
N
MN <0
2m BAt
by taking expectations and iterating the recursion. Therefore, it holds that [
E
||I
Fi.
E
Xx] ||
|x++) ArE [X+
1| Xx} + 24t(d 1)
AC
A (A
2Att(d 1)
and
sup E & 2 < At A BMN + 2Ch k00
+ 22tt(d - == MAP + MA 0o
N
with Mj == A + BMN + 2C, and M2 2(d 1).
Similarly, let us define (xx)) = (b(X)) + +xxx (t. Then, under the similar proof, it holds that
E [I 12 Fx E [B(x)) 1 Xx] ((x)) + A [|(e((x))1 Xx.
At (A B1| At
N
and
20h N
sup E [[x(())) A20
< At A BMN +
= Mr
(10)
Next, we prove the decaying property of h Rd - R. From the 2nd order Taylor expansion on h, the following holds almost surely:
1(X++1) = ((xx) + Vh(X.) (X))
- (1-aA()(X)) - Tr (U(X. A (XX At
+
- MXXXX
A
where X₁ is some point E Rd between X and X and AXe = XK+1 Xt. Therefore, by applying thet previous observations, we get the following:
B(h(X++1) F2 I (1 - (AtM + L [IB(X)) X + LAFE [|AXA| Xx]
2
which implies
LAP (M 22t) + 2M2)
E[((X++1] S < (1- (M)))E)/
28

---

By applying the telescoping sum with expectation, we recover the following formula:
[[((K)] - - ast)* E [[((X0]] L + 2 + 2M22
2xx
Finally, we note that
MN = o d+ At 1+ (1+ n)), M1 -1=0 (d+ At 1 vx) + )
(a+
M2 = (d)
and, therefore,
At
d+ 1 NN) ))
E [[((KR) (1 aAt)* [[h(X,)] 0
a
As we can see in the above proposition, the effect of N scales with (1/NN) and vanishes as oo. Therefore, the usage of the Hutchinson estimator does not affect the convergence speed of constraint functions under the single-equality scenario with sufficiently regular h and f. Remark (Geometric control of ALST). We state the following result from geometric analysis that provides a lower bound for SLTI:
Theorem A.1 (Informal, [49, 47)). Let E, be a compact Riemannian manifold with diameter D and non-negative Ricci curvature. Then, À 2 1+2D-/X or A 2 [1+2x) holds, where x is the first eigenvalue of the Laplace-Beltrami operator on Σ.
x²
This theorem implies if the constraints shrink the diameter D (or increase \ of Σ, the lower bound of Si increases, so the on-manifold decay KL (px (ps:) accelerates (dominated by its exponential rate -AAssı appearing in Theorem or Theorem 3).
29

---

B Algorithms of OLLA, OLLA-H, and baselines
Algorithm 1 Euler Mareyama discretization of OLLA & OLLA-H
initial point 20 R step size At, number of steps K, landing rate a, boundary repulsion rate potential f. constraints {4. {gz (g
Input:
Hutchinson probe numbers N, mode
E (OLLA, OLLA-H)
2: Output: sample trajectory {=*]
3: for k =
1
4: Evaluate constraints: (1x(2x)) J (1 (Ox
5: Compute: (f(x),, J(zz) (ax)
6: Compute: G(r+)) (I(aa) [((x)]
7: if Gan <
La: then G + G else G-1 t G-1
8:
if mode = OLLA:
Compute: Tr [T (II( Vhh ...
9: else
Compute: Tr ビk=1 i0 (xrk) Wiiik (I* V 9i 2/ 2x2², u NO Ta)
10: Compute: H(Tx) + VJ( (mw)
Compute: (Ixx) [I((x)) fz())) G H(*x)
Update: Ja+1 t 26 + ((xx) At + (2ArI(( where & NN(,,a)
13:
end for return {
14:
Remark 7 (Pseudo-inverse of Gram matrix when LICQ fails), A second issue of OLLA arises if the Gram matrix G = VJVJ" becomes singular. This may happen when the LICQ condition momentarily fails near the neighborhood of Σ. In that case, we replace the inverse G-1 with the Moore-Penrose pseudo-inverse G Because G still projects onto the row space of V. and annihilates its null space, it enforces the same exact orthogonality to constraint gradients, preserving the exponential landing behavior and numerical stability of OLLA.
Algorithm 2 Constrained Langevin (CLangevin) with slack variables [19, 20]
Input: initial position To E Σ, step size At. number of steps K, potential f. constraints { {g/}/=1 projection iterations tolerance т, regularization
2: Output: sample trajectory {=x 1K
3: Initialize slack variables: $0.j + max 2g, (zo), OJ for j E [4]
4: Set the extended state: 8o + (8o.1...., ...,, 80,1) E R yo + (zo, $0.1....,..) ) E Rd+
5: for k= O..... K do
6: Draw & ~ (0, la+e)
7: Compute the augmented constraint vector:
1 (xx)). + 91 (Ix) +
1
(yn, =
2 .........
2 sk
8:
Constrained update:
Uk+1 Uk Viss V2 VJ(UK)
where (ex s) = (xx and Ås is chosen such that [| S T by Newton's method with regularization À and max iterations
9:
Set : 2x yx(1 d), 8a + ma 1:
10: end for
11:
Output: (zr)k
30

---

Algorithm 3 Constrained HMC (CHMC) with slack variables [19, 20]
1: Input: initial position 20 L, step size At, number of steps K, potential f. constraints
{hr Wi (ha)m--119)4 )=))) projection iterations L, tolerance regularization X, friction
2: Output: sample trajectory {==]ā=o
Initialize slack variable: $o.j f ( (mo) 0f for j [4]
4: Set the extended state: 80 -- (50.1; 80,1) E R! yo f- (Iro s0.1, 80.1) Rd+t
5: Sample momentum po ~ (0, La+z) such that (J(o)) 0
6: for k= 0.....,, K
7: Draw Ék. 1k+1/2 NO, Ia++)
8: Compute the augmented constraint vector:
1 (m))) + 2 . ....:....911....... (......((((())) + s4,1
1 2
J(yk)
2
9: Midpoint Euler step:
AL
7P% + PR+1/A) + AtrEE + VJ((k) (
Ps+1/4 Pk
4
such that VJ(( 1/4 =
Verlet step (1):
At V fext (yk) VJ(W)) >x+1/2
Pk+1/2 Pk+1/4
Wk-1 yk + P1+1/2
such that [| S by Newton's method with regularization À and max iterations L Verlet step (2):
11:
At V Text( WJ
Pk-3/4 = Pk+1/2
such that J( Df Pk43/4 =0 Midpoint Enler step:
12:
(PK43/4 + Pk41) V AM + VJ(M4+1)\
PŁ+1 = Pk+3/4
such that (z(++))) Pk+1
13: Set : 24 f 8k [1 d), 8% + [u + 1:
end for
15:
Output: {xx) kk k=
31

---

Algorithm 4 Constrained generalized HMC (CGHMC) with MH correction 21]
1: Input: initial position 20 E E, step size At, number of steps K, potential f. equality constraints [h4]2=1 (9s 13 projection iterations tolerance 5, regularization X, friction y
2: Output: sample trajectory {ax 1K
3: Sample momentum ~ (0, Ia) such that VJ(o)) = 0
4: for k = 0 K 1
5: Draw Ex- Sk+1/2 ~ No, Ia)
6: Compute the augmented constraint vector: J(xx) = (Ix)) hm (xx))
7:
Midpoint Euler step:
P4+1/1 Ph 1(p: + P1+11/A + ANTEE + VJ(())) (
such that VJ(Ik P(+1/) =0 8: Compute the Hamiltonian: (I P f(x) 2 1lpp+1//l1B
9: Verlet step (1):
At Vs(x)) VJ(xz)" A641/2
Pk+1/2 Pk+1/4 2
5k+1 3k 1k-1/2
such that [[ ol S by Newton's method with regularization À and max iterations Verlet step (2):
10:
At
Pk+3/4 = Pk+1/2
Vf((++1) + (J
2
such that VJ PK
11:
Compute the Hamiltonian: ( PN43/4) f(*++1) + 211/443/41BB
12:
Cettopolis-Hastting Correction: With probability
min (exp (-(H(ix+1 P.43//))) H(K,,Pk41/1)
= (Tx+1 if g( 1) S 0
Otherwise, reject and flip momentum (x+3) (0k: - Pk+ ;) Midpoint Euler step:
At (Pk+3/4 + pk+1) + / At + (J(
Pk+1 Pk+3/4
such that ( Px+1 =
14: end for
15: Output: {xx 1K
k==
32

---

C Proof of Basic Properties on
C.1 Intrinsic Gradient on >
Recall that our target manifold 2 is defined by = {a E R h(x) = 0, 0,g(z) 50). For smooth function f Σ + R. the Riemannian gradient Vzf is defined by the relation
(Vzf(z),e) (o), for every T
To check (z)) = (x), observe that ITx) E T.E and [ reproduces dfr on every tangent vector. This is because E T.E. ((()) = (V)
(V)()), 1) (v) due to the symmetric nature of [I())) Therefore, the intrinsic gradient on E can be extrinsically defined by Vzf(x) = [(x))
C.2 Intrinsic Divergence on
Recall that if X x(22) where X(E is set of smooth vector fields on Σ, then the divergence of X on E is defined by divs: (x) - (m+ 1..)) ( X, Ex) where A is Levi-Civita connection and {
is an orthonormal frame of T. To check its extrinsic formula, first observe that the induced Levi-Civita connection is given by V (xx) == (VrX) = my Xx) where X, Y E (E), T indicates tangential component Σ, and V is the Levi-Civita connection (or usual directional gradient) in Rd. Then, if X and E.......((mm+ are extended to the ambient space Rd holds that
d- (-m++Ual LaB
d-(m+1Ze) Σ T=1
divex= =
Σ
X,
(IVE,, E (IIVI)
where the last equation is obtained by taking basis {E1, ... Ea. (m+ +.D: 21, +1a1 of Rd with
(i- being the orthonormal basis of (T and applying the definition of trace.
C.3 Recoverable Tubular Neighborhood of Boundaryless Riemannian Manifold
Let U.(E)= == fx E Rdd dist(z, Σ) < be tubular neighborhood of Σ with reach whose existence is guaranteed by the compactness of [50]. In the proof of Theorem we require property that enables us to recover the unique x E (E) such that n (x) and h(r) = p, given the information of
and p E R..
The following theorem indicates the existence of such a nice neighborhood of Σ such that any x in this neighborhood can be recovered from the information of y, In addition this, need regularity lemma to connect the decrease of [(x)]]2 ||2 with the decrease of dist(z, E), which is crucial property to show the convergence of the W2 distance.
Lemma C.1 (Regularity lemma). Assume h E C² and LICQ condition is satisfied on E. Then, there exist constants é, such that for all x E U+(E) C (E)
[M(x)]]2 xl|æ mx))|2
where - millyeN Jmin (Vh(y)) 0.
Proof. For each y E E, Vh(y) has full rank by the LICQ condition, so its smallest singular value Ty = Omin (h(y)) 0. Since, any x (E) be decomposed into y for some
y E E, N, (y) = E) Taylor's theorem gives
h(x) hl + v) Vh(y)e+ (v)
where the norm of the remainder term (v) is bounded above by [R,(())12 1M for M supyeu. (V < co. Hence,
[|h())]l2 |1(2 v) l2 z IVVh())ll2 - [ 2 o,l|1112 1 M
1
2
33

---

Furthermore, because h E C2 and y + min (Vh(y) ( is continuous function on Σ, the compactness of Σ implies there exists o m 0 min (Vh(y) > Also, choose f e such that M - a.
Then whenever x: E (E) so that
||0||2 < é, we get
[|((z)]l2 2 a||u||2 - 2
1 M M 2 (a 2
Melluula
a |elll = a || - オ (z) 2 2 2
Lemma C.1 shows that an upper bound on ||A(z)]] I|2 yields an upper bound on || π (2) however, this guarantee holds only once z has entered Ug. Therefore, to ensure that x indeed enters Ua whenever [/(xz)]22 is sufficiently small, we appeal to Lemma C.2.
Lemma C.2 (Entrance cutoff of U2 (E) in terms of l1) Assume h is continuous and coercive, i.e. [[(())]l2 + Do whenever | + oo. Then, there exists 8 0 such that dist(z, -) if [1(())l22 < 8 where è > 0 is the constant defined in Lemma C.I.
Proof. First, we show that h is proper map. Let C be compact set in R that it is closed and bounded. Then, (C) is closed because h is continuous. Also, suppose (C) were unbounded that there is sequence { ken h- hC) with ||æ - oo. Then, by coercivity of h. (1x Il2 - oo while every h(zx)) lies in C, which is bounded and leads to contradiction. Therefore, h (C') is bounded and closed and therefore is compact by the Heine-Ba Moonel theorem. This proves h is proper map.
Now, define a set S² == for E TRd dist(rr Σ) 2 2 and assume inf [1(())]22 = this case, there must be a sequence {ma)eeN C S2 with [| + This implies Ia E h (B(O, 1)) for Vk 2 K for some K E N. Since h h(((Bo 1)) is a compact set, there is a subsequence {T } of { converging to some H, from the Bolzano- Weieessrasssss theorem. Subsequently, the and x. E Σ. However, because Sy is closed and every 2R, S₁, it should satisfy I. S = dist(æ, x. E) 2 é, which is contradiction that I, E Σ + dist(z, M) O. Therefore, 8* = infaes, [1( > 0 holds and dist(z, 2) é if [lb(x)]l2 < 8
continuity of h implies h(xx,) limy- h(xx,)
Theorem C.1 (Recoverable tubular neighborhood). Let D= fx Rd h(x) = 0} be a compact and boundaryless Riemannian manifold with LIQQ condition. Then, there exists recoverable tubular neighborhood of E with width 5, Ú&(E) = {a E Rd [1]b(=)l2 < <6} UE(E for some
>> such that
1. The nearest-point projection map Ús(E) Ex argmin y is well-defined.
2. The following recovery map ((y.p) is well-defined
((y.p) 2 8(0,8) + Ús(2), where L ::E: x B 5) + R" is the C function such that h(((y. p)) = p and x(C(y,)) p)) = y for ((..p) E Σ
= Vhly) L(y.p)
(11)
E 8).
In this case, we refer to Úg as the recoverable tubular neighborhood of Σ.
( Then, for Proof. Define F((,, L) R (d+m)+m + RT by (), L) Vhly + L)(y)) each y Σ, F(y,, 0) h(w) - - 0. Also, V.F((,,)) implies VLF((y, 0), 0) = Vh())Vh(y)T, which is invertible due to full rank assumption of Vh(y).
Therefore, by the implicit function theorem, there exists by > O and an open set Uy = (l/ E V x RM dx(y. < 6y. IIp'll2 < d C Rd such that there exists unique C' func- tion L: Un R satisfying Lyy 0) and F((U.P), I(y.p)) = for v(y.P) E Uy-
Now, observe that UsESU, the open cover of Σ x {0), which is a compact set. Therefore, using its finite subcovers, we pick > 0 such that for Vy E Σ and ||l||2 < 8. L(y is
34

---

well-defined C' function on Σ x B() and the recovery map ((y.p) : E X B(0,8) + R
<(y.p) Vhly) Lyy p) is also well-defined and C² on Σ ((0, 8) such that ((C(y. p)) = p holds (from F((v.P) p). Lyy p)) = 0).
Finally, setting 5 = min { (where the constant 8 comes from Lemma C.2) and Ús(E)
{ E Rd [1(x))]22 < < C (E) enables the well-defined nearest-point projection map π on (x) By applying x to (3 p). we recover the formula ((((,,)) p)) = T + Vh(y)* L(y
y. V( 6).
C.4 Recoverable Tubular Neighborhood of Riemannian Manifold with Boundary
The following results generalize the proceeding lemmas and theorems by incorporating inequality constraints alongside the equality constraints.
Lemma C.3 (Regularity lemma with boundary). Let Σ - E
Rd h(x)
O, g(x) S 03. As-
sume h,g C2 and LICQ condition is safisfied on E. Then, there exist constants é, K. >> 0 such that for all E UE) C UU(E
[|b(z)]l2 + 11gr (x)ll 2 x||æ- - <<>)ll2
where a. minger miin ([vh(v)T, Vor, (v) > 0.
Proof. For every subset C [I]. consider E, = {v E Σ I, = 13. Then for y E E,, VJr(y) = [Vh(y)T. (gr()) e IR (/)+)/ has full rank due to the LICQ condition, so its smallest singular value of miin (v)) > 0. Because y of is continuous
E,, also have min inf o > 0 by Lemma C.4. Also, since any x E (E) can be decomposed into y + v for ICl y
some y 2, E (2:). Taylor's theorem gives
((x) = hly + v) = Vh(y)u + (v), 91, (x) = 91, (v + v) = Vgi, (w)v + (v(()
where the norm of the remainder term (v), (v( is bounded above by [ S M I| 12, RR (v)ll2 M|lll for M max (2),
(1)V2)=)]l2. ( 00. Now,
set Vy))), which satisfies ||lm|2 2 o||u||2 by the definition of a. Then, observe
(0)
|1((z)]l2 + 119r, (z)ll2 |/l|2 (||l)22 - [|((2) 1gr, (z)]l2) a||u||2 M M|ol
where (o) comes from the following observation:
||ll| - (1)(()) 1/ + 19.r, (x)]12) S || l2 - (1,( (x)) S I [M(=)* .1,)),,,,(()) (x)
[1) (1) S |R (x)]l + (v)
< M1
Now, choose small 0 K such that Mē < g. Then whenever (+(E) so that ||u||2 get
[l(xx)]22 + I1g (2) > _ ||lll2
- g1 1|5-=()))le
Lemma C.4. For every subset I C [1], consider E, == {v Σ I, = 1}. For E E,, de-
fine vr((2) = [Vh(y)? FFT e R(m+)4))xd and of = amin (v)( > 0. Then, inf of holds. NEEr
Proof. Suppose inf of 0, then there exists sequence yk Σ, with Jmin (VJr(yk)) ye>, Because Σ compact, the sequence has a subsequence of yı converging to y. E by Bolzano-Weierstras theorem. In this case, there are two possibilities:
35

---

1. When 1g. = I, mmin (Jr(y.)) > 0 by LICQ, contradicting the assumption inf of = 0.
yy
2. When 1g. = I K for some non-empty set K C [4, the LICQ condition implies V.JuK (y-) = [Vh(y+)". Vgr(v=)", (yx()=)* has a full rank, thereby, imposing (g)( have full rank. This again implies inf oy > a contradiction. yE>a
Note that 1g. does not deactivate already activated inequality index E I. This is because gi (y)) = for Vi E / and continuity of gi implies gi (y-) = 0. Thus, the previous argument on two possibilities completes the proof.
Lemma C.5 (Entrance cutoff of U2(E) in terms of [|((z)]l2 and g(x)), Let E = {i e Rd h(z) 0, g(xx I 03. Assume h.g are continuous and h is coercive, i.e. [1(())]l2 +
whenever | 3C. Then, there exists 8* > 0 such that dist(æ, 2) < é if [[((2)]2 < 8" and 9((x) < 8* for i E [1] where > O is defined in Lemma C.3.
Proof. Define S == { E R4 dist(r, E) 2 }} and let (x) = max (114(z)) (z), 9 (x)). Assume rns (xx) = 0. In this case, there must be a sequence (Ixheen C S² with [[(())) 12 +
and (x())) S (xx) + for Vi E [4]
This implies Iy E h (B(0,, 1)), vk 2 K and for some K E N. Since (B 1)) is a compact set due to the properness of h (Lemma C.2), there exists a subsequence {2 }of {ax} converging to some H, from the Dollzano-Weierstrn
theorem. Subsequently, the continuity of h implies M(x.) = lim h(F,,)) )= = 0, (z)) < 0, Vi E [4], therefore, x. E E.
However, because S₂ closed and every 2g Sg, it should satisfy x. E.. + dist(r,, Σ) > which is a contradiction that x. E E + dist(æ, Σ) 0. Therefore, 8 := infres, (x)) > 0 holds and (x, Σ) < é if vxxx) < 8
Theorem C.2 (Recoverable tubular neighborhood with boundary). Let Σ {i E Rd hxx) = 0, g(x) S {0} be a compact Riemannian manifold with boundary. As- sume LICQ condition on E. Then, there exists recoverable tubular neighborhood of E with width 8, Ux(E) == f Rd [[h(x)]l2 < 8,9(2) < <8} Ue(E)) for some 8 >> O such that
1. The nearest-point projection map a Úb(E) + E,T(x) = argmin, 1122 is well-defined.
2. The following recovery map C(y, p- 41 well-defined
C(y. P. q1,))= : Σx (0, 1, (0, + Úz(E), -+V
where Ly is the index set of active inequalities at y E, J(y) (h(y)) 91, (v)] e m=+,,,, and L Bm (0.6) (0,, Im is the function such that
(I(((1.P.91,)) =4. i I, 9.(S(v.P99,,)) <0. I
h(C(v.P.91,)) =p,
T(((y,P, 91,)) =
V(u.p. qr,) E x Bm (0, x Bu,, (0,6). Furthermore, when E in((E), L = L(y.p) is a C function Σ x Bm (0,
Proof. Let first define F((,, p. 91, ), L) Idd + mm+,,, F((,, P. 91,,,, (hly+ + VJ(yy) - P: P. (y + VJ(y) gr,)) E m and consider the stratum E, EVE E I, = 1} for each subset / of inequality indices.
Then for y Σ,, we observe that F((,,0,0).0 = 0 and V.F = VJ(I) which is invertible by the LICQ condition. Therefore, the implicit function theorem ensures the ex-
istence of 6, 0. U, == ((V' ...pp'gh Σn x Rm x R2 dx v') < 6y, 1|/'ll2 < Sy, 1l 6,)
36

---

such that there exists unique C' map L : Uy + R=+ satisfying ((U..P. O, qr)) E Uy-
Now, observe that UhEn, is the open cover of E, x {0)m. x (0), which is a compact set. Therefore, using the finite subcovers of Σ,, we can pick 8, 0 such that for Vy E Σx, |lpll < 61, and l|gr|12 < 8y, Lyy p. pr)r) is well-defined C function on En Bm (0, x (0, and the recovery map Ç(y, 41) Σ x Bm (0, 8) (0,) + R. ((v.P. = + VJ(()) Lyy also well-defined C map on Σn x Bm (0, 81) (0, 8,) such that h(C(y. p. ()))) = p and ge (e(CV. (:))) = 9; for E I, which comes from the property Fly. 1, (.,)))) ((y, p. 41)) =.
Furthermore, from the continuity of gs, i & I. there exists y/ > such that (2(2) - - for all e By, 71) and the compactness of Ug(2) gives G-Lipschitzness of g. (2() for all E [4] for some
G > 0. Also, we recall that L(y. (1) is C map on { (v.P. qr)
Er, llpl2 < 8, 1lrrll2 < &
with Ly,, 0, = O.
Hence, the Taylor expansion of L(y. p. q with boundedness of | V V, 112 | P. on (V.P..r) y E,, Iplla < 8,, lellll ll < 81 and boundedness of [VV on Σ gives
| P: gr) < Cr(llll2 + 1lrrll2) = [[VJ(N)" Liy.p.ar)ll s Cj + ||
for some C₁, Cj > Then, choosing 8, == min (8r, 1d44
concludes that whenever llll2 <
o, llq l2 < 8,, the inequality
(.(CV..P.4)) = gi(y +
& gl(1) + GlV 1 21 < 0 2
holds for all i & T. Finally, setting ō = min {minrau 5,, 8%} (where the constant 8* comes from Lemma C.5) and (E) == {o E Rd [Ih(z)]la < 6, (l(())]l2 < < C Ue(E) enable well-defined nearest-point projection map on Úa(E). By applying to ((y, P. 9r,),, recover the formula (C(,, [ (g,,)r - y, V.y,p.rr,)
Bm (0,0) (0,), O
D Construction of SDE with Exponentially Fast Decaying Constraints
Proposition (Construction of OLLA and its closed form SDE). Consider the following SDE:
dx, = (Xx)dt + Q(X,)dW.
(12)
where
SQ/,. = 0, Vi = [m],
Q = argmin|| 22 - Q11 s.t
[QV 0, Vj E l..
QEaaXd
[Vh 1 (T("))QQQ)) + ah, =0, Vi E [m],
9 argminllg + ffl GERd
s.1
[Vg " TT(('9; + a(gg +c) =0, Vj I
Then, there exists a closed form SDE of (12) given by:
dX, = (I(())VV + G
+ H(X,)dt + V2TT(X,), (Ito) V2T((x)) dW, (Strato.)
dX, = -IX(X.)V +
where o denotes the Stratonovich integral and
H := -==vV/ (T("), Tr (V") Tr Tr gi, 1) ...., Tr(V gizz m)]"
is the related Ito-Stratonovich correction, mean curvature of f x Rd h(x)) 91. (x) = 0}.
Proof. Define J(z) := [Mx=)T. 9T. + | E mI and VJ(x)= = V9) (x) e
R Il)x) Then, the Lagrangian function associated with the optimization problem for Q becomes
37

---

[(Q.A) = 1V21 + Tr (A"VQ) where A E R(++ 1,Dxd is Lagrangian multiplier.
Then, the stationarity condition with respect to Q gives
aL = - - AJVJ aQ
Q= 1 VVV.
0
Also, the constraint condition implies
1 ATUVV 0 -=QV--- (V21 2
= -V2V (V.J.
where VJV.J is invertible due to the LICQ condition. Therefore, we get optimal Q = =
v2 (I - V.J (V.VJ"''~VVJ)) For the part, we set
1 b: 2
Tr (v<2; ah1, Tr (V Rinz QQ") + a +0]]
and the associated Lagrangian function L(g, : being an Lagrangian multiplier. Then, the stationarity with respect to q gives
= | V
with A E 1 4,
0L
+ àq
=
q= /= A
Again, the constraint condition implies
0 = -//=b==v///-VVVVV-
A= (V
using the invertibilitty of JJV By plugging this expression to 7. we recover the optimal ==
IIV VJ'(VJVV''' b. Therefore, the Ito version of closed form SDE given by
dX, - -I(((X) - V.J(X,) -(X,,,))) +
V2I,XX)
where G VJVJ" is the associated Gram matrix. Also, some tensor-calculus computation (Equation 3,46 in [19]) gives
- -V.J" [Tr h1) ..... Tr h Tr
Tr (V m) VII()II((x)
which is the Ito fio-Stratonovchh correction term. Furthermore, applying the technique in [19] (Remark 3.17), we can recover that this expression is equal to the mean curvature term H(x) of a manifold defined by Σ,, {a E Rd hxx)) 0. (x) 0}. Therefore, we get the following closed form expression of the SDE:
dX, = - V(X, + OV/
+ H(X,) + (ili))+++2III(. VII(X.), (Ito)
+ d.X; - (V/ dt v211(Xe) dW,,
(Strato.)
where H. - (T
..... Tr (V* II) Tr (V²g,, II) Tr (V m)]] is the associated Ito-Stratonovich correction term (or mean curvature term of Σ,,).
Lemma 1 (Exponential decay of constraint functions). The dynamics induced by (12) satisfies the following properties almost surely for Vi [m], Vj Txo:
hx(X.) hx(Xo)e at 120
(13)
and
1
9 (x() +e
9x(X,) - (X)(Xo) aar
< In
a
1
(X,) +e
(x,))) 5.
1
In
with (X,) < 0, Vt 2 0. j & Ixo, where Is = {k E
[4]
9u (z) 2 0} is the index set of active
inequality constraints.
38

---

Proof. Observe that, for each & E [m], the Stratonovich chain rule implies
dhg (X,) (X.) ( - (XX)) VV2((X,) o dW:
@ Vhh (
m+ 1s
a Σ (X)]) -acodd--ahhsh,aa -ahx(Xi)dt,
where (1) holds due to the fact Vhy 11(()) = By integrating both side with respect to t, we
recover h(x,) = h(Xo)e-o tt 0 almost surely. Repeating the same calculation for 9j for k E L. we obtain
m+1+
dg (X,)
Σ
c)dt -X(xx(X)) c)dt,
which again recovers 9: (X,) = -E + (Xo) + e)e-at Furthermore, once (Xt) - 0, it is instantaneously reflected into interior of E whenever it hits the boundary Oc. Therefore, gx(), S holds for Vt 2 0, Ix..
E Proof of Theoretical Results Equality-constraint OLLA
Observe that when the constraints are only equality constraints, the equality-constraint OLLA (3) given by
dx, (IX)) (AVN)//XX)))+)XXXXXX)00()(())))) (2I(X))
(2I((X) dW,
(14)
The high-level proof idea of Theorem is to decompose the convergence analysis into two parts: (1) Convergence of W22 distance between Pt and Pr. (2) Convergence of KL between Pr and ps where Pr, Pe are the law of Xx, (; (X.)) respectively and pr; is the law of the target distribution, which satisfies dpp x exp(-f(z))dos
E.1 Upper Bound of W2(pt, Pt)
Lemma E.1 (Upper bound of W2(Pe, Pt)). Let Pt be the law of Xt which follows equality- constrained OLLA (14) and define ta 1 In (3)). For i. 2 la the law Pe Y, (X,)) is well-defined and it holds that
a
Mn aat
We(pt, s
(15)
k
Proof. For lo. observe that ||X, - Yell2 = |X, - (X,)l12 <
(1)((x+)l22 I Mae-at by Lemma C.1 and Lemma 1. Then, by integrating both sides with respect optimal coupling of pr and Pe. we get
(o) I
M e
W2(Pi, Pe) & [E[l), Y
E E
Yella|
s
at
where (o) holds by Jensen's inequality.
39

---

E.2 Upper Bound of KL 1LL (p:llps)
Lemma E.2 (SDE representation of projected process). Let Xt be the stochastic process following the SDE:
d.X, b(X,, t)dt + (X)) dW,
where X E ÜE(E), Vt 2 0. Then, the stochastic process Y, defined by Y;= x (X,) follows the SDE below:
dY, (Y)))(( ()di (2)(()) dW, (V,(X))(X.)) -xx))))....VVVVVVV....
+ (Xx (X))))((Y)) dW,
Proof. From the the Stratonovich chain rule, we observe that
dY; = dv (X,) = V=(XX)6((,, t)dt + (X, a dW. This expression can be re-written as follows
dy, VT VT(Y)) + (2V dW. () VT(Y))((X,)))
+ (2 - V=(X)) o odW
[I(Y:)6(Y, t)dt V2 o () dl
(xx - (Y(()) o where (o) holds because VT(y) = (1)).. 11²)²² = l1l (idempotent) for Vy E Σ and Y, E Σ.
Corollary E.1 (SDE representation of projected process from equality-constrained OLLA), Let X be the stochastic process following equality-construiie OLIA (14). forn,for! L 2 lo(:= In (t)) the projected process Ye := т TT.)) follows the following SDE:
dY, [-II()))V bn (Yr + VI AN dW, where (V,, 1)122 Chye-at, AV,, CANE-at for 2 almost surely for some constant Coy, CAN = CLA Mr
> with Ct, being the Lipschitz constant of (x(x)[((>) on
Ü.(
Proof.
By applying Lemma E.2 to the SDE (14), holds that dY, = (TY)) t) + 22T(())))))))))))))_)))))( o dW, + (V,(X))(XX.)) (x(())((Y
+ V2
Vx()) dW,
for b(x, 1) b(x) = -VV (=))))))))))) By using Lemma Theo- rem C.I, we can set X, = C(V,, ((x))e where Ç E x Rm + (E) is the recovery map. Now, since X,, Yı Úz(E) and the closure of Ús (E) is compact, (T(z)b(x) and VT(z) is CL., CLA Liipschizz on (s)2), respectively for some
0. Therefore, it holds that Yell2 S Cr, Cr,))M(Xo)lla eaat Cr, Mr. e-at ((D)e ((x))
(Y,(Y,t))lz S Ct, CL,,,((Y,
K where (Y,,) VT(((Y, (T(s(Y,, (oxx))--a)b((Xx (((),,,))((() -aee))a)) VT(Y)6((Y)) and the second last inequality comes from Lemma C.I. Similarly, we obtain the bound of AN,Y, = h(Xo))),(,X))eee VT((( h(Xo)e (oooee-a)))((((X at)))U((((,,)))))))))))))))))))( h(X)) -aa))) (V(()))((())) follows:
e < CL, My
(A, f)12 S CLAII((Y, e - Yellz <
Finally, we complete the proof by setting C.y = Cr,M CCAy = Cz Mr and observing that VT(Z)= [I(T(z)) for Vz e Üs(E), which implies AN(Y,,t) INY,))
The following theorem is a Fokker-Planck equation of a Stratonovich SDE defined on Riemannian manifold. We will rely on this theorem describe the time derivative of pt.
40

---

Theorem E.1 Pookker-Planck equation Riemannian manifold [51, 52]). Lel X₁ E be stochastic processfollowii the SDE:
dX, = Vodt + Σ V dBt
1=
where Vo, Vk smooth vector fields E for each & E [d] and Bj are kth components of Brownian motion BN. Then, the law Pt of the stochastic process Xt satisfies the following Fokker- Planck equation:
Dep: (p: + 2
divs(divvs(::
Lemma E.3 (Upper bound of KL (pellsss) Assume that ps: satisfies the LSI condition with constant AiSI. Let X. be the stochastic process following equadivy-consrriined OLLA (14) and pe be the law of Ye == π (X.) after 2 lears leur == max { I In(Cs)). Then, for a + 2À1st. the following non-asymptotic convergence rate of KL 1KL2(pm (p:llps) can be obtained as follows
-2Xcss(t tww) 2desiCs (e (e-at -a-tte e -eaaww)) IKL (p lpss)
KL (p:ll=)) S exp
a
L exp 22(s t((m) + 2des1Cs, (c-ae -as- - e e-ahw)) as
+ C
lms
In particular, a > 2d151. it becomes
KL (pelle)) & exp (22 ttea) 21rs1Os (E - -E-ata)) ) IKL (ptwellp)) + Crl for some constants C3 = 0(1+ + CAN + -CAy): Co. Cy = Cee allo
a
0.
a-2sg
Proof. By Theorem E.1, Corollary E.1, and the choice of Vf In pe, we know that the projected process Y; given by
dY, [-I((Y)) by (Y,, t)] dt (2I(()) -W dW,
and its associated Fokker- Planckk equation can be written as follows:
0.p: divv (Vx Ps + bN)) Σ dive (divs 6f) + 5%))
where Sk NT, 6x following equations:
IT and en is kth standard basis vector for Rdd. Now observe
aKL = 1 Oepr In
Ôt dos + d / pidos =
- Ope In
Pt
dos
dive V Impe) +
divs
In
dos
ps
1=1
Term
Pt
+
- dive(ppbbx))
In
doy
Σ
pz
Term
Pr dos: + divs:(?.d2)b0)
+
[div:: (divss((div + divs:(?:
Σ
Ps
Term (3)
41

---

Analysis of Term (1) induced by the main SDE.
From integration by parts, we obtain
Term (1) /
Pt
Pt (iv)( Vv. In
pe (vsllpp
Vn
dos
dos
x
p
ps
=
{ Pi in Vg In
Pt
dop
Σ
VEPP, V
V
Pt
dos
1=
d
Pt
(faiv Veln In
dos
pr:
k=
(A)
At
( Pt
Pi (V Vyln In
ddos
Vz
dop
Ps
Pel
Vv
In
Pt
[gdor
(p-lle)
3
where (4) holds using Lemma E.4 (the third term 0) and the fact that P: VE In Vgpt.
Analysis of Term (2) induced by the noise drift bx. Again, using the integration by parts, we observe that
Pt p=:
1/5 pe(b,, V5lIn
Pt ydos
(Term (2) =
div (p.bx) In
doy
pr
Pi
11ados S Caye-at [ p
Pe
S DellbwllalVvx.1n
1/ In
llaos
ps
Coxe at {
Cox
+
1
-
In
Pt Px:
11 dos: e aat 1= (Pellps) + -bve-at eat
Cb
Obx
where (D) incquality holds using the AM-GM inequality.
Analysis of Term (3) induced by noise diffusion Ay. To analyze Term (3), we apply integration by parts and the chain rule of the divergence, i.e., divy (prak) = (Vyr,, ax) ((vis for a vector field 0g on Σ:
Term (3) = Σ
/ d
V In
P+ dos - Σ
Pr pr.
dop
k=1*
pr:
k=1-
5.
pe dos
pe((divs:fe(0&
Vn
doy
Σ
(V Ssk
Vz
3c
k=
Σ
IVIP (6k, Vz
Pt
|dss
Σ / (V (6 V In
d
P: dozz
pz
k=1
k=
Now, note that Ekm = IIANII Edm SKOT = II II, and Ed J. = IIA II. Then it holds that
Σ k=
& In
pr
doo
=166
İr dop
pr
S CANe-al aat [ i Vs In
Pt
112ds
Pr
(X)
x CANe at (Pollp))) + Ose oot
po
V In
Pt
12dox
Pp
(0)
S Ca ( + CAN Cze
pr
C3
+
In
Pt
11
doy
4
C3
CA. C2 -nt
s 2ot +
4
42

---

where C3 = max In p=lle < oo by compactness of Σ, and (x) holds using the triangle inequality; In l2 < V In
fu
ll2 + IIVr: In Also, (o) comes from the AM-GM
inequality. By following the same reasoning, the following inequalities are obtained:
Pt
S e (Pellpp)) + CAN att
fx) Vn
doy
pr:
4
5t
An C3
Vxrte (6kr
Vyln
ddy
220 AN e 2at -2,,-22 (pellp))
2at
4
C2 C3
20 AN at - +
4
Next, we observe that the following terms decay exponentially fast:
Pt
divs (fv (
pe(divx((ft0) Vx In
dor
1
pel
Σ
es
(4)
Pt
V In
Dr
S CAN
pol 0o.1-11V2
112d0s
Pr:
{ p1
Can
1
SCA
V
11
doz
4
CAN
C2
se ae-at1= (Plles)) +
cat
where C4 max divs (I(x)() xo. (0) holds because dd k- dfvx(fx)o; = (II) ILAN
once (Iv) is given by (div == div=(fa) ) for each k E [d]. Similarly, we obtain the following bound by using Lemma E.5:
Pt
Ps Pc.
गe (divs VscIn
Pr [ In dan <
|
A=1
LemE.5
S Carve-at at S VillV In
pt P
O2 aat
1220l
re +
and
d
Σ
Pe (div)(b))0 Vz In
ydoz S / Pe | )12 Vx In
Pr
Idor
Ps
kml
Lem S Cayo 2at [ Vo In
pr |idax Se at (T.(lllez))) + Oztn -e-
(ame-an eat
4
where div is similarly defined by (divgILAN), == divs( for each k E [d].
Applying Gronwall-type inequality.
By summing all the bounds, we arrive at
1St -2ALsi(1
0,KL 2Pi
(poll(p) + Cge c +
(c²= + Cay
c3
C² C2,
0a c2 + Cg. +
Co == with Cs = (4+4CA, +2C +2C2 C.
+
Therefore, the Grönwall-type inequality gives for t > tout, tewt == max 1 Ind, In(Cs)):
-2Aysi(t tcu)) 2dcsiC6 (e (a-at_ -attaaeaa c-ttm) IKL E (pt.m llpz)
KL <L*(Pellp>) S exp
a
22Ls1(s - loar) + 2ALSiO', (e (e-as e =daa)) e casds). a
+ Co
exp
43

---

In particular, a > 2dLSI. it holds that
-2xrss(t - lem)) 2ALsIC' (e-at - e-aau)) |KL (ptca (pps)) + Cil
KL (pello)) S exp
where C7 == 222si
from the fact that
[ exp
(2Nisi(s toui) + 25 (0- - -c-aaa e ds =
00 exp(2ALsI($ tcal)e ds
lar
e -aao
OC.
a 2ALSI
Lemma E.4. Let [fk] be set of vectors defined by fx [Ix))s where IT())) is the orthogonal projector onto T, and eg is the kth standard basis vector TRd. Then, it holds that
[divsfu)fk =0.
Proof. Recalling that (Tx)) = 1
we define N(x) =
Vh(x)) (Vh()))))())) E Raxm so that ( = Im and I1(x) = / (z))
we let the columns of (x)) to be (x) ..., 7m (x) then these produce an orthonormal basis
of E,E. This is because Im(N (x))) = Im (m(hh(z)²) (from the invertibility [Vh(:)\V(>>>">>=)
implies (m) (x)) span N and ((x)"N((3
guaranteess the orthonormality.
Next, we define a vector field F(x) by F(x) (lx))div=(I(())) where (div):l((z)) ==
dive(fx(I)) for each k E [d]. With this definition, we have divyll -divs(NN") -h divg(nAn))) Now observe that for E [d].
(divs(manA)) Tr (IV
ni (men
n T mmamm IisnazO5mal
(divgn))ma + Σ II) 8 = (divgn))ma)
where nu is the Ith component of ng. From this fact, have the following result:
div divs(mun?))
(divesme)m = F ITdivs:(II = A=1
Σ
Finally, the definition of gives Ed div F, which zero by the above argument.
Lemma E,5, Let Pr be the law of the projected process Ye of X,, where X, follows equality- constrained OLLA. Define (z AN (x. t)ee and denote dIVx:(IIAx) vector in Rd such that (div>llAN); == div: (6x) for each & E [d]. Then, it holds almost surely
[IVV) (Y, t))ll2 S Caice
for 1 2 la(:= 1 In (+)) and some constant Caiv > 0.
Proof. First, for each & [d], observe that Vo,
+ AI t)ex and
div):(xxv)) Tr (II))) AN ()xx)
VAN
(16)
Term
Term (2)
44

---

where, for matrix -aatue function G(y). VG is the chird-orded tensor defined by (VGiyu)))s (G,(M) for j, k [d], and the gradient is taken over y.
Oye
For the Term (1), we know that when y = Y,, 2 to.
Term (1)) - | - [[V SKIlAA/ S KICANe
where K1 == sup - -
oc.
For the Term (2), recall that (v t) (C(((y. ((c)) Vx(
conditionally on Xo. from Corollary E. For the notational convenience, we define m(y.t) Slv. h(Xo)e Then, it follows that
() ()) t)l + (ll V²T(t)[[(() VT(y)
ViT(N(V ())(v) 1) (1
+ V'r(n(y. (n) () + t)
VI(y)
At this moment, from the recovery map < in Theorem C.I, the integral form of the remainder gives
Vn /|2 || (S(v.)) C(v, 0)]l [,,(u.p)) V,C(u, 0) S K2 K2||p|12, where == ((,o)e-at K₂ == sup (x V Vv,(Vv p)ll <
By combining these results with the previous expression of VAN (v/, t), we get
(V,, 1)12 - D4e at
Üs. the Lipschitzness of (V1, VIVII on Üg and the contraction of [[,((,,)) Y2)l2 S eat Finally, when y Y,. we get the following upper bound of Term (2) by applying the previous result:
for some D4 O, using the boundedness of || ||V |2, V 22, |T
Term (2) [[V & -d- ml V AAcell2 Vd mD4e
D5e a.s.
where (o) comes from the fact that | = Tr (II())) = rank rank(II(y)) = d - m and D3 == vd - Vd. Therefore, by combining results for Term (1) and Term (2), we obtain
ddiv (oiv(6)(Y)))) < Doe
(Il)) (Y, 1)
div (dv) Dze
for 1 2 to and D6 = KICAN + D5, Dy := VdDs. Because the final result holds without any dependency on Xo, the result holds almost surely without conditioning on Xo.
45

---

Lemma E.6 (LSI implies Talagrand inequality [43, 53, 191). For the probability measures μt,
defined on smooth complete Riemannian manifold E, define the W distance between ja and E by
3
W (1, = inf [
dx(p. (dp, du)
<Ell(.))) JEXE
where (H, v) denotes the set of coupling probability measures of μ. and dy denotes the geodesic distance on so that for p, E
{[,
C7E =1,7(1) (p)>))))) = }
dz(p. = =inff
with I1-1lg being the induced metric on E. Then, the probability v is said to satisfy the Talagrand inequality (T) with constant Dr >0 if for all probability measures j with μ < V, it holds that
W (1, 1)
KL (H)l)) 'T
2
Particularly, if v satisfies a Logarithmic Sobolev Inequality (LSII) with constant NLSO then v satisfies the Talagrand inequality with constant λLSI-
Theorem (Convergence result for gquality-constrainet OLLA). Suppose assumptions (CI) to (C4) hold. Let X, be the stochastic process following the equality-constrainet OLLA (14)
and let Pt, Pt be the law of Xt and its projection Y. ((x)) on E for 2 tews teut max 2 In In CsS, respectively. Then, for all t 2 tema it holds that
W2(pe,
Mn +
2 ALsI
KL - (pellp))
K
where
KL (Tellps) S exp
2drs1Os -22is((1 ((w))
(e (c-at- -a -calmm)) IKL " (pr.
a
+ Co L" exp 22(s tewr) + 2des/Os (e-as - - admn e aad
t..
In particular, a > 2dtsI. it holds that
KL KL (Telle)) S exp (-2desi(t - bleu)) 21eSiOo (e (e-at -e E [KL (ptn + Orl
for some constants Cs - (1+ Oza Mn + (2aMMe ) Co, Cr
Cae 2--2Ar
0 with CLa
being the Lipschitz constant of (T(())[[(2) on (.)
Proof. First, observe that d2(p. = 1pppp 4/1a S ds(p, q), Vp. Σ because Σ is the submanifold of ddd with Euclidean metric. Thus, W2 ps) I (pt, holds and we have
A A-ineg S
W2(pi, PS)
W2(pt. p.) + (2 Pz) S We(p.. + W (pt.(pp)))
Now, recall that Σ = Rd h(x)) = 0} is a smooth compact and connected Riemannian mani- fold. Therefore, it is complete by the Rinf theorem. Thus, Lemma E.6 implies
KL KL (Dtllex) = W2(Pe. S W2(Pts + V KL (L"(pillpz>)
2 ALSI
WH (pt,))) S
Hence, we conclude the proof by borrowing the results of Lemma E. and Lemma E.3
46

---

Proof of Theoretical Results Inequality-constrained OLLA
In this section, we analyze the non-asymptotic convergence rate of inequallty-coonstrained OLLA. Note that Proposition gives
dX, (X, vdW,,
if ( <
dXx, = ((x)d aVg.. + cll, dt + 22T((X)) dW,, otherwise
as the closed form SDE of inequality-cconttrained OLLAA
F.1 Convergence Result for
OLLA
Lemma F.1 (Boundary behavior of Pt meequlltyy-onsstaall OLLA). Let X2 be the stochastic process following the constrained OLLA and Pt. be the law of Xt. Also, denote Jj be the
probability current density defined by d.pt = - Jt. Then, for t Z teas tem == In M, it holds that
(n(x), Je(x)) = O, Px(x) - O, VI E O where n is the unit normal vector of Σ and daz is the surface measure of 0c.
Proof.
From the Fokker-Planck equation the inequality-constrained OLLA, we know that
2 J = qp: - V !QQ pt] = (QQ")) P. Vpt
where the last equality comes from the chain rule of the matrix divergence. Then, for each Vgj. j E Iz. observe that
VØ (V. (QQ"))
Σ (ke0)) (V (0Q"))&
E
(0ko,) am
t,k=
d
0
Σ (QQ2)inOup
(QQT)) = -T
1=1
i,k=1
holds for x E OE, where (o) holds due the Vgl Q condition of Proposition 1 and Q² Q. Therefore, for each j E Iz, we have
Voz q -r (V (V
(A) -alg + +e)p; = -aep 0,
(J,, Vggi = V J.
where (A) holds from the Vg,9 Tr QQ) + alg + c) condition of Proposition 1. Finally, we conclude the proof by observing the following:
dt / pidoy / v J.do
d
/
n" J,daax
aep, doos P:
Joz
where n is the outward unit normal vector of 0E and (+) holds because supp(p:) C E for 2 feut implies Sg pido for t 2 feur-
47

---

Theorem 2 (Convergence result for the LSI condition with constant disI. Let X be the stochastic process following the inequality-
OLLA). Assume that ps: satisfies
constrained OLLA and let pt be the law of Xt. Then, for t 2 teur tew IIn Mote the following holds
W2(pt, PF)
2
KL = (pellpE)
V Arsi
where
KL (p: -e
2Aeslll-bao) KL (ptcollss)
Proof. Observe that
8,KL (p:lle)
-f 5x
Oep: In
Pt ps)
doy + +0
pidos:
Opr In
Pt
dos
ps:
(-V J1)
ln
Pt
dos
J.V In
Pt
do
(J
In
doop
5
ps
05
Ps
J.
J.V.
In
Pt
dog (A)
-I=(pellp>)
-2Lss
pr
where (o) holds due to Lemma F.1 for t 2 feul and (A) comes from the fact that Jj V In p=
Pe V In P = - V
Pz
almost everywhere in Σ. Therefore, we recover the upper bound of
KL x (pellpx) by applying the Gronwall-type inequality follows:
KL (polle))
-2t tea) (
Also, we recall that Σ f E d g(x) < }} is smooth compact and connected Riemannian manifold, thereby it is complete by the Hopf-Rinow theorem. Thus, Lemma E.6 implies
We(Ps S W (Pc. P.) S V ALSI
KL (pellee)
Hence, we conclude the proof by applying the previous upper bound of KL 2 (pillpz))
G Proof of Theoretical Results Mixed--onsttrined OLLA
G.1 Upper Bound of W2(Pt, Pt)
Lemma G.1 (Upper bound of W2(pt,Pe)) Let Pt be the law of X, which follows mixed-
constrained OLLA (12) and define t == max { (Mate) 1 In ) For t 2 t, the
a
law Pe Y; Y; = (Xe) well-defined and holds that
Mh at
W2(Pi, Pi)
(17)
Proof. For / 2 L1, observe that
Mh at
|X, +(X)1la < [1(Xi))22 + 19 (Xe)le
by Lemma C.3 and Lemma 1. Then, by integrating both side with respect to optimal coupling of pe and Pt. get
We(P p.) < [E - r.lli)) $
Mr
E II
Ydlla <
aat
K
where (o) holds by Jensen's 's inequality.
48

---

G.2 Upper Bound KL E. (prllps))
Corollary G.1 (SDE representation of projected process from mixed ceeed-conttraieed OLLA).
Let X, be the stochastic process following mixed-construined OLLA (12). Define l ==
max { In (M+te) ) 2.1 In (())) and assume Y, (Xi) int 2 via the assumption
La
(MI). Then, the projected process Ye follows the following SDE:
dY, = ====)=1())VF((x) + (Y,, t)] dt V211(Y)(1 + AN(Y,,t) dW..
where 1|bx (Y,, (Y,,t)l2 = C-at "ayeea,,1AAA Aw f)ll = Cane-at for t 2 I almost surely for some constant Oby CAN = Gr ddn > with CLa being the Lipschitz constant of (=(x)[[(>) on Ü.(E).
Proof. By applying Lemma E.2 to the SDE (3), holds that dY, II))))(()))) t)dt + (I) dW, + (V((X))(XX ()(())((Y + V
(Y,
for b(x, 1) = b(x) == [mxx) (xx) +
Now note that T (Xx)) E int(E)
fort 2 and the recovery map < (Theorem C.2) is C' for (()) E int(E), (E) Therefore, can X, = -(Y,, (Xo)ee -at) by using Lemma and Theorem C.2. Again, since Xx. Yi Us(E) and the closure of Üs (E:) is compact, (x(x)((x) and (x(x)][(x)) is Üza CLa -Lipschitz on (2) respectively for some CL,y ČLA > O. Therefore, holds that
V2ll2 Cr.11A(Xool1 -at Org c-at
[V,x( (xx,
VI(())6((() and the last inequality comes from Lemma C.I. Similarly, we obtain the bound of AN(Y,,t) (V((Y) as follows:
where (Y,tt, = V (((
h(Xo)e "IC ((X))e
C2.ol((X)llz car CL, Mr e-t
| (Y, S CC,,ls(Y, "t) - Yoll2 5
Hence, we complete the proof by setting CoN == Ma CAN = Ô. Mn and noting that VT(z)
(<T()) for 7 (xx) int(E), x Ú8(E), which implies AN (Y,,t) IN (Y,,
49

---

Lemma G.2 (Upper bound of KL (prlle))) Assume that ps satisfies the LSI condition with
constant ALSI. Let Xt be the stochastic process following mixed-constrained OLLA (12) and Pe be the law of Yt == (Xx, after 2 Leur, henr == max In Mate
1
In
Ma 8
) 1 In (n(Cs)
Suppose
(Regularity of Ep) Ep== h(x) p. g(x) S 0}} C inr(E) for 0 llpll2 S 8.
(Regularity of PPp) The boundary velocity u of PMp satisfies sup ue || I V 1p for some V 0,3 >
O. Also, assume My: =- supl|pll <6 oas, (0E,) oc.
(Bound on Prs Pr)
G1 = suPr20,EEE Ps 30 and 0 < G2 S Ps S for E
Then, for + 2A1s,, the following non-asymptotic convergence rate of KL s (pellors) can be obtained as follows
-2ALsi(t - tew) 221sIO' (E-ot -E - x
KL - (pille)) S eexp
a
[KL (Peullox) + f" exp tow
lonu) 221s1(s - leme) + 2desiOs, (c-a0 - e e-ata)) x
[(ōs + aG,GoMMn) E e GVV ass ds}.
In particular, i 2desI. 33 2 1, the inequality becomes
KL (pelle)) S exp -2>esr(t - lcar)= 2desiCs, (e (e-at_ e e-of IKL (Ptmllll) ((pe) + Ö7 Cz.
for some constants ÖS = 0(1 + ČAN + C, Ga, Gs, Go, Co Ör > 0. and Er = (Čo +
G. -o
MG,G,MM) 2Ais and Cs == (G&V M, ag 2A-2
Proof. First, Corollary G.I, and the choice of V. - In ps: gives the following SDE of the
projected process Y, of X, for t 2 fouti
dY, (V- -
dW,
where 6N(I, t) == Vr(((z, (X VT(((Ir, ((Xo)e =t)(), h(o)e Vx(r)l((ar) for x E Σ, conditionally on Xo-
Vz(z)b(r), Ax
Hence, from Corollary G.I. its associated Fokker- Planckk equation can be written as follows:
a.p: (Pv (VyImp + bN)) divy (dv( + 6x))
+
Defining E, : T(f2: E TRdd h(xx) = m(xo)e -at .g(x)= s 0}} conditionally on Xo. we observe
0,KL (pille)) = 0 /
p
p. in
doy
E.
/ 05,
Pt pr:
fel ne)dooss + / irdos
O.p, In
dos +
Pr
In
1
5,
Term (1)
Term (2)
where the last equality holds from the Leibniz integral rule with of being the velocity vector of the boundary of Eg and being the outward unit normal vector of ΦΣ, Therefore, the expression of
50

---

O.p: implies that Term (1) becomes
Term (1) / ( piVVs + bN)
Σ (Pe (&)) + +5k)
Vy in
Pr Pr
dds
Term (1-1)
-So. es,
Pe(Vsslpx: + bN)
d
Pt
div dive (delF + +z))(fk + 5k)
In
.n
Pzz
Term (1-2)
by integration by parts, where fa = Ilec, 82 = IAAddR. Now, we observe that Term (1-1) can be bounded
Term (1-1) S -1 -11- ( + Čge-at
following the same proof of Lemma E.3 with different constants ÖN 0(1 +(2x)x)): Co (note that we ignored the integrand at Oc which is measure zero with respect to doy).
For the analysis of Term (1-2), first observe the following fact:
S
Pt ps
G
Gi
}
pe In
dooz, I G3 max
In
oas, (OΣ,)
E
G2
G2
G
S G3 max
ln
= G4
G2
from the assumptions of G1 := suP120,zeR ग < 00, < G2 S ss S G3. and the regularity of QE;
such that sup12ta. Tas, (0E,) < Ms: < co.
Next, from Corollary G.1, we note that the following holds conditionally Xo:
fz(2)) + (x())) (I) + (t(t,)))e
tangent vector on for E Rd h(x) - M(Xo)e .g(x) S 0}. (x) + (z) = VT(((x, h(xo)e (()))e="))I[(((( h(Xo)eeat ")ex becomes a tangent vector of Σ, Similarly, it also becomes tangent vector of 0E, on the boundary because II is the orthogonal projector induced by and active g. Hence, (Sk + oc, nz) = 0 holds, where n, is the outward unit normal vector of ØE,,
Because (()o)er)
Therefore, Term (1-2) becomes
Pt Ps
Term (1-2) & aG5Maeeat
Pr In
doar. S SaGaGs.Mneaat
where (I) holds because Corollary G.I gives
V px)) + (x, (CCa000-w..E.DDoo.o./e///C"""""0
with b(x) == - -aaaa
and, therefore.
(x + (as
& aG5Mpe <0.
for x E 0,, with G5 == sup, llllllll (C(((((
Also, similarly, Term (2) is bounded as follows:
Term (2) < sup ||,/1/2 / (1zz Jas.
In
Pr
+ +
do S Go &GV
rEaE,
with Gg = G4 Mys Therefore, combining the results with the LSI condition gives the following inequality:
0,KL (pille))) < -2AIsT(1 -ae-aa -LL (prlle)) (Co + aGGGGsM)) eat GGV M -38t
51

---

where the last inequality comes from the LSI condition. Finally, applying the tröp inequality recovers the following inequality:
2deSiC' -e-a---aaer)
KL (prllps) exp
-2Lssi - teur)
[KL ( lpps)
a
[ exp 221SI(8 lem) 2desi (c-as e-otam) x
(Co + +GGGGGMM) as + GGV ags ds).
Also, similarly in the last argument of Lemma E.3, observe that if a 2ALsI and B 2 1
exp 2ALSi(s tem)) + 2ALSiC'
(e (a-a - e as ds $
exp
(2ALSi(S temr)) c ds
a
c -aton
<000
22LSI
and
exp 2ALSi(s - toui) + 2dLSICo (e-as - --aatu) e e- ds 5 C Caffte
< 00.
a8 -221s1
fox
Therefore, there exists C7, Ög < oo such that
2ALSi (e (e-at at k e-aton) |KL (Prullo)) Ö7 + Cz].
KL (Telle) S exp
-22is((t - lem)
a
where c= (Co + aGaGg.M.)
atta 2Ats
and Cs = (G&VM))
a8- 2ÀLsi
52

---

Theorem 3 (Convergence result for mixed-constrained OLLA), Assume that pv: satisfies the LSI condition with constant ALSI- Let X, be the stochastic process following mixed- constrained OLLA (3) and be the law of Yt = *(x)) after 2
leat, leut
max
In
Mate
In Mu 1 In(Cs)) Suppose
(Regularity of Exp
Ep: mf h(a) (( p.g(x) 50)) C inl(Z) for 0 < ||p||2 S 6.
(Regularity of p) The boundary velocity i of S..p satisfies SuPEEdE, ||2/ S V 11 for some V > 0,3 >
O. Also, assume My = supp)p] <o oas, (0Ep)
(Bound on Pts P.)
G1 == supr20.ee Pt xo and 0 < G2 S Pr S G3 E
Then, for a + 2λLsi the following non-asymptotic convergence rate of W₂ (pt, p2)) can be obtained follows
Mr
2
W2(Pr, pr) - K
+
KL E (prllps)
AISI
where
2desrOs (e 1ea e-atan) x
KL (polle)) S exp
-2Xcsi(t temr)
[KL (pr. + I" exp 221s1(s tam) + 2dLSIOs, (E (a- E -atm)
[Coo (Co++GGGGMM, cas G&V G&VMRe -aßa] ds)
In particular, if > 2d1sI and B 2 1, the previous bound simplifies to
KL (pelles) S exp (22csr(tt tew) 2 (e- - - -e-aum))
( (pps) + Cr + Cs.
for some constants C₂ = (1 [+ Na (ELaMMn
( G1, Gs, Go, Co 0. and C7
(Co + aGGGGG aaaa 2AIs and Cs = (G6V Mp
e -altoo
with Oza being the Lipschitz constant
a88 2Aiss"
of (x) on (E).
Proof. We note that (x) = 0 on OΣ for t 2 four holds from the regularity of the Ep assumption. Therefore, by the same approach Theorem 1 and Lemma E.6, holds that
2
2 ALSI
W (pt, S
KL (palls)) + W2(Pe & Wa(pt, Pc) +
KL KL* (pr
ALSI
Therefore, we conclude the proof by combining the results of Lemma G.I and Lemma G.2.
53

---

H Experiment Settings and Supplementary Results
Settings. The first two experiments were exccuted on desktop with an AMD Ryzen 7900X CPU (12 cores) with 32 GB RAM. Runs were implemented in WSL2 (Ubuntu) environment (CPU-only), using the Python and the PyTorch [54] framework.
H.1 Experiment Settings and Supplementary Results for Synthetic 2D Examples
Experiment Settings. In this experiment, we compare four samplers (OLLA, OLLA-H, Clangevin, CHMC) on the following synthetic 2D examples:
1. (Star)
star-shaped equality manifold with uniform density:
f(x) = 0, h(x) = x1 + (1.5 + 0.3 cos()))) 0 = arctan ( T1).
2. (Two Lobes)
two-lobe inequality manifold (from [18]) with uniform density:
f(x) = 0, g(x) - - gx)) 2, q(x)) e 2( 33) +e +
2(1 3)2
3. (Quadratic Poly) quadratic curve defined by mixed polynomial equality and inequality under a standard Gaussian target:
1 f(x)
- ||
12. h(z) = 2 + 2 - g(x) = =x x22 - 1.
4. (Mixture Gaussian)
nine-Gaussian mixture restricted by seven-lobe manifold:
f(x) = In Σ
exp - - cc1b ) h(x) = 22 x3 -(3 + c70)
1=1
g(x) = (x) 2)²
5212 + 0.5z. 40, 0 aretan 2112
with (c:) 9 = (- 2,200 0,0,2) 23²
where x [x1. E RR2. For each 2D example, we run 200 independent chains for K = 5000 steps each. From each chain, we retain only the state at step K, yielding 200 samples per sampler.
To provide a fixed target distribution for distance calculation, we generate 200 samples using CGHMC. This reference is held constant across compute all comparisons. For each sampler, W2, energy distance as well as the mean con- straint violations EE| [t (z)]] and B|maxx(22) (]
Table 6: Hyperparameter settings for 2D synthetic examples (At x 10 10-4)
Method
Hyperparameters
OLLA
a = 200, 6=
over the 200 samples.
OLLA-H
200, 6= 1, N
The hyperparameter setup is provided in Table 6. To mitigate ill-conditioning in the Newton solver of CLangevin, add Tikhonov regularization with Tikhonov matrix r vXI, À 1.0 for the Mixture Gaussian example, and X = 0.1 for the other three. For the Clangevin, CHMC, CGHMC, X, is initialized exactly Σ to ensure the stability
CLangevin
L = 3, T -= 10-1, =1 (1,0.1)
CHMC
7 1, 3, r 10-4,1
CGHMC
7= 1, I. I 3, == 10-1,A=0
of algorithms while OLLA and OLLA-H have noisy initialization Xo = Yo + (O 1), Yo e E and X, progressively approaches Σ by the landing mechanism.
Remark 8. The results shown in Figure and Figure were obtained under different setup described above: larger step size A/ was used; all methods were initialized from the same initializa- tion point Xo E Σ, rather than via random sampling near E; and the regularization parameter À in CLangevin was increased for improved numerical stability under large step size.
54

---

OLLA
OLLA-H
CLangevin
CHMC
Energy Distance
W! (Sinkhorn)
Mean [h(X]]
0.08
0.05
0.04
0.06
0.03
0.04
0.02
0.02
0.01
0.00
0.00
0.00
3000
5000
Iteration
Iteration
Iteration
Figure 6: Convergence diagnostics on the Star example (equality-only case). (1) energy distance to CGHMC samples (left), (2) W distance to CGHMC samples (center), and (3) mean of [h(x) (right)
OLLA-
CLangevin
Energy Distance
W? (Sinkhorn)
Mean max g(x)
0.05
0.05
0.04
0.04
0.03
0.03
0.02
0.01
0.00
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
1000
2000 3000 Iteration
4000 5000
Iteration
Iteration
Figure 7: Convergence diagnostics on the Two Lobes example (inequality-only case). (1) energy
distance to CGHMC samples (left), (2) W? distance to CGHMC samples (center), and (3) mean of max g (right)
OLLA
OLLA-
CLangevin
CHMC
Energy Distance
W (Sinkharn)
0.08
0.10
0.06
0.06
0.04
0.02
0.00
0.00
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
Iteration
Iteration
Mean of [hOx)
Mean of max gix)*
0.03
0.03
0.02
0.02
0.01
0.01
0.00
0.00
0
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
Iteration
Itoration
Figure 8: Convergence diagnostics on the Quadratic Poly example (mixed-case). (1) energy distance
to CGHMC samples (top left) and (2) W distance to CGHMC samples (top right). (3) mean of [m(z) (bottom left) and (4) mean max ( (bottom right)
55

---

OLLA
OLLA-
CLangevin
CHMC
Energy Distance
W3 (Sinkhorn)
0.08
0.20
0.06
0.15
0.02
0.05
0.00
1000
3000
4000
5000
0.00
2000
3000
4000
5000
Iteration
Mean [h0O]]
Mean of max gix)*
0.05
0.20
0.15
0.10
0.01
0.00
0.00
2000
4000
5000
lteration
Iteration
Figure 9: Convergence diagnostics on the Mixture Gaussian example (mixed-case). (1) energy distance to CGHMC samples (top left) and (2) W distance to CGHMC samples (top right). (3) mean of [h(z) (bottom left) and (4) mean of max g(x) (bottom right)
Supplementary Results Sampling Accuracy & Constraint Violation. In addition to the Mixture Gaussian example (Figure 3), we evaluated OLLA and OLLA-H on three further 2D geometries to verify that the distributional accuracy observed in Figure of the main text generalize other settings. Across all three additional examples, OLLA and OLLA-H closely match the convergence behavior of CLangevin and CHMC in both energy distance and W? metrics, while maintaining constraint violations at low levels without requiring explicit projection steps. These trends are illustrated in Figure 6 (Star), Figure 7 (Two Lobes), Figure 8 (Quadratic Poly), and Figure (Mixture Gaussian under the hyperparameter setup of Table 6). These additional experiments confirm that our landing-based sampler provide relatively accurate, ressprinn samples across diverse manifold geometries.
Supplementary Results Effect of Hyperparameters and In Table 7, we report how varying the landing rate a (with e = 1) affects sampling performance on the Mixture Gaussian example. For both OLLA and OLLA-H, increasing from to 500 leads to consistent reductions in energy distance, W?, and the average constraint violation [=]](=)]] However, setting a too large causes sampling failures and numerical errors (e.g. 5000).
Similarly, Table examines the effect of boundary repulsion rate e (with a 100). Across full range of tested, sampling accuracy-as measured by both W? and energy distance- essentially invariant, exhibiting smaller variation than when a is changed. In contrast, inequality- constraint enforcement steadily improves increases: the average violation E[g(x)+1 declines monotonically, reflecting stronger repulsion. Only when e becomes exceedingly large does one observe degradation in equality-constraint enforcement and occasional numerical instabilities, mirroring the breakdown seen at an extreme a value.
56

---

Table 7: Effect of a on energy distance, W3. E (())) and B[g(x)+ on the Mixture Gaussian example with 1 (top: OLLA, bottom: OLLA-H)
ar
energy distance
W2
[[
Elg(x))1
1
0.121=00025
0.363+40.64
0.682+0.017
1.11340.351
10
0.066=0.019
0.200+00035
0.130
0.234400032
100
0.003:10.016
0.159
0.017:10.00
0.045
200
0.040=0.012
0.121
0.008++0.001
0.054400009
500
0.033=#0011
0.104++0.200
0.004+0.000
0.001
700
0.04410001
0.132110019
0.003:00.00
0.016+0.011
5000
NaN (results unavailable)
1
0.104+0.019
0.333+0.077
0.643++0.02
1.102+0.337
10
0.059=0.018
0,181
0. 0.120100.11
0.189+00023
100
0.052
0.156100026
0:.000
0.055+m.024
200
0.052+0.018
0.153+00037
0.009+00000
0.039
500
0.041
0.124+0.027
0.004
0.013+0.009
700
0.03710.011
0.110100019
0.00210.000
0.00710.006
5000
NaN (results unavailable)
Table 8: Effect of on energy distance, w3. E[]h(r)]],
[]. and B[g(z)*
on the Mixture Gaussian
example with 100 (top: OLLA, bottom: OLLA-H)
energy distance
W3
B[(N(=)]]
[[g(x))+
0.1
0.048+0014
0. +0.00
0.01440.01
0.082+0.017
1
0.033++.004
0.108
0.017+0002
0.067=0.027
5
0.040+40.06
0.123++0018
0.016+00001
0.040=00015
10
0.036410016
0.11240.034
0.017
0.019
50
0.038+0004
0.112++0029
0.018+0.001
0.003=0.004
200
0.041 10.200
0.119
0.018
0.00610.02
10000
0.066110022
0.137
0.033
0.0000...00
0.1
0.039400012
0.126+00014
0.013+0000
0.073
1
0.048+00018
0.142+0.029
0.016+0000
0.048=0.014
5
0.035.10011
0..00
0.018 10.001
0.0271..007
10
0.04410014
0.127
0.018
0..02710.007
50
0.04740013
0.134400019
0.018+0.002
0.010=0.016
200
0.04010007
0.117:10006
0.018
0.001
10000
0.027
1.111
0.083
0.000
57

---

H.2
Experiment Settings and Supplementary Results for High-dimensional Manifold with Large Number of Constraints
Experiment Settings. In this high-dimensional experiment, we construct synthetic "stress-test" manifold in IR by imposing m - 1 linear equality and / spherical inequality constraints inside bounding sphere. Concretely, we generate m random hyperplanes h;(x)) = af afxx (with a; NO0 Ia), b; ~ (N 0.12)) for i E [m 11. together with the sphere constraint hm (x) || R²²
(with = 5), and spherical obstacles (z(x) = ,² || - cgll2 (with r = and obstacle centers Oy ~ NN R for j E [1] All randomness is fixed via a seed across experiments.
For each choice of ambient dimension d, the number of equality constraints m, and the num- ber of obstacles I, run one single chain of each algorithm for 1000 iterations. We discard the first 200 iterations burn-in and then retain every 5th iterate, yielding 160 post- burn-in samples. Similar to 2D synthetic ex- periments, the baseline samplers (CLangevin, CHMC, and CGHMC)
Table 9:
Hyperparameter settings for high- dimensional manifold example (At = 1 x 10-
Method
Hyperparameters
OLLA
200, ==
OLLA-H
a = 200, <= 1, N =5
initialized exactly on Σ so that Xo E Σ, where as OLLA and OLLA-H start from noisy initialization Yo-
CLangevin
L 5, 10 = 0.1
CHMC
= 1, I. == 5, = 10-1 0
CGHMC
7= 1, == 5, = 10 +,1=0
Xo + A (0, La), Xo E Σ.
We measure performance along two complementary criteria. First, report the computational cost per effective sample size (ESS), defined as
CPU time ESS = total CPU runtime (s)
minigega ESS;
where ESS; is the uni-variate ESS in coordinate i. Second, we assess the estimation accuracy via the sample means of representative test functions- P( > 0) and some complicated test functions. To isolate the effect of each problem parameter, we vary one element of {d, m, 1} at a time while holding the others fixed.
Remark 9. described above: a larger step size At was used and the regularization parameter in CLangevin, CHMC, and CGHMC was increased for improved numerical stability under large step size.
The results shown in Figure 4 and Figure 5 were obtained under different setup as
Supplementary Results Scaling under Dimension and the Number of Inequality Constraints In Figure 10, we plot the performance of samplers as the ambient dimension d increases from 50 to 700 (with mn 1 = 5). Across all d, OLLA and OLLA-- produce virtually identical estimates of our benchmark test functions similar to the baselines. In contrast, CPU time per ESS of OLLA-H stays essentially flat (on the order of Rs 0. 0.05s/sample), whereas OLLA grows roughly lincarly (reaching 1.1s/sample at d 700), Also, the CPU runtime plot exhibits OLLA-H achieves the lowest wall-clock times across all dimensions.
Similarly, Figure 11 shows results with dimension fixed at d = 100 and the number of inequality constraints varying form 10 60 (with m = 5). Again, OLLA and OLLA-H approximately match baseline methods in estimating the mean of test functions, and OLLA-H demonstrates dramatically lower CPU time per ESS and CPU runtime than both OLLA and the other baselines.
Supplementary Results Scaling under the Number of Equality Constraints In Figure 12, we show how sampling performance changes as the number of equality constraints m grows from 10 to 60 (with d = 100, / = 5. fixed a 200), Although OLLA-H continues to show the lowest CPU time/ESS (and total CPU time) among the methods, its sampling accuracy gradually degrades as increases, drifting away from the baseline values. Especially, we observe that the equality constraint violation worsens, and one must compensate by increasing a, reducing At with longer chain to suppress the equality constraint violation. Each of these solutions may lead to increase of computational cost, and a cannot be driven arbitrarily high due to its induced discretization instabilities (Table 7). Equality-only baselines therefore achieve more accurate estimates-albeit at higher cost-indicating that large m regimes are particularly challenging to OLLA-H compared the high-dimensional or many inequality constraints settings.
58

---

20)
2.0
0.7
1.0
0.5
0.3
0.0
100
200
500
B00
703
100
200
300
400
750
300
400
500
600
di (dimension)
(dimension)
d (dimension)
Blx +x1 3xg
21 sinxx 16% log(is]+ + 1Jtanhrxg)
costx))
5ie sle-niccosaag) + txt x(lsintss))
1.0
0.8
0.8
0.5
0.6
OLLA
OLLAH
0.2
Clangevin
CHMC
0.0
-0.2
600
(dimennionn)
Figure 10: Scaling under the change of ambient dimension d. (1) CPU time per ESS (top left), (2) total CPU runtime (top center). (3) estimates of test functions (others) as the dimension d increases from 50 to 700 (with = I = 5).
CPU time/ESS
CPU time
P(x; > 0)
OLLAH
Clargevin
CHMC
CGHMC
50
15
50
20
30
(inequality count)
/inequality count)
[inequa fiy count)
21 [[sinix le log((ks) 1)tanhix;)+ costx))
xIx! x5x +1
ZIe 2le=/ocoosxgg) (x3 xfjsintag]l
1.0
2.5
2.0
1.0
0.5
0.0
20
30
[inequality count)
(Inequulity count)
[Bnequa Byy count)
Figure 11: Scaling under the change of inequality count L. (1) CPU time per ESS (top left), (2) total CPU runtime (top center) (3) estimates of test functions (others) as the number of inequalities [ increases from 10 to 60 (with d = 100, 5). The average inequality violation is maintained [] = 0.000 + 0.000 for all samplers.
59

---

CPU time/E55
CPU time
Pig
OLLA
8-
CLangevin
80
CHMC
60
OGHMC
40
0.4
20
0.2
0-
10
30
10
20
30
40
50
10
30
40
50
m (equaltly count)
m lequaltly count)
m lequaltly count)
21sinDr))eft log((ss 1)tanhixal Costell
[[x] xf +
#le-hicoo(ic)) (x2 -115500x)]
10.0
7.5
0-
5.0
-1
2.5
0.0
40
40
40
(equaltiy count)
(equaltiy count)
(equaltiy count)
Figure 12: Scaling under the change of equality count m. (1) CPU time per ESS (top left), (2) total CPU runtime (top center), (3) estimates of test functions (others) as the number of inequalities m increases from 10 to 60 (with d = 100, 5).
H.3 Experiment Settings for Molecular system and Bayesian logistic regression task
Settings. These two experiments were executed Linux machine equipped with Intel Xeon Gold 6226 CPU (24 cores) and 192 GB RAM.
Experiment Settings (Molecular System). atoms R so the state space dimension is 3aatomss- Let p E Rd be the flattened vector of atom positions, and let Pk E I denote the position of the k-th atom (where p is reshaped into Natoms matrix).
This experiment models polymer chain of Natoms
The equality constraints h(p) = 0 enforce fixed bond lengths (l)) between adjacent atoms and fixed angles (0,) between consecutive bonds:
1. Bond length constraints (k: O,.... Natoms 2):
Athnnd (p)) |P PR+1 =0
2. Bond angle constraints (k- 1,.... Natoms 2): Let of = PP--1
Pk and uy Pk+1 Pk.
hangle,<(p)
cos(Oa) =
12
with I1 = 1.0 and 0g 109.5..
The inequality constraints (g(p) < 0) enforce steric hindrance, ensuring minimum distance (r'in)
between non-adjacent atoms (ji - jl 2 2):
Piy (p) 2 P P,l <
for j 2 2
with mmin 1.0.
The potential function fp) models the energy of the polymer configuration and includes terms for torsion angles and nondbedded interactions based the Seatss-hhhann-AAeeeaa Andersen (WCA) potential. is calculated as f(p) = (U)rr((p) + (p)), where 3 an inverse temperature parameter (default: 1.0).
The torsion potential, Unr. depends on the dihedral angles the (for k: = 1,..., Natoms 3) formed by consecutive bonds, where dy denotes the dihedral angle between atoms Pk-1, Pk. P8+1, Pk+2-
60

---

The potential is sum over modes E M (default: M = {1.3}) with corresponding force constants km (default: 0.5, k3 0.2) and phase shifts 8n (default: 61 0.0, 63 0.0):
Namu
(p)(p)
Σ MEM
km (1 + cos(mmx 6m)).
1=1
The non-bonded WCA potential, Unh, models repulsion between atoms i and j that are not directly bonded and are separated by least two bonds (ii - أl 2 3). Let Rij I|P: P be the distance between atoms i and j. With the steric minimum distance minn its associated length scale mmin/ Tm/0 and the potential cutoff is = 21/ a. The interaction energy is given by:
U2)( 4EWCA ( Raj
12
Rai
+
EWCA
where EWCA is the energy scale (default: 1.0). The total non-bonded potential is the sum over eligible pairs (j 2 i + 3) where the distance is less than the cutoff::
Na
(p)(p)
Σ
Ui(RR)) 1(Rig <Te)
2=
J=4
where I[(-) is the indicator function.
We vary Natoms {5, 10, 15, 20, 30), corresponding to d {15, 30, 45, 60, 90). We run single chain for K = 5000 steps, discard the first 1000 as burn-in, and thin by 5. To measure the accuracy of sampling, use the radius of gyration squared (Rg) as test function, which is defined
Naaas
(p(
|
Pom l1
where Pam
Ps.
N
Natoms
N
Natoms
K=
k=0
Table 10: Hyperparameter settings for the Molecular System example (At 10 10-5)
Method
Hyperparameters
OLLA-H
a = 500, e = 1.0, N E {0,5)
CLangevin
L 30, T 10-4,
0.5
CHMC
1.0, L = 30, r = 10-4 4, A = 0.0
CGHMC
7 1.0, 30,
4 4.- 0.0
Experiment Settings (Bayesian logistic regression). This experiment involves sampling the posterior distribution of weights 0 E RRd for two-layer Bayesian neural network applied to the German Credit dataset [41]. Let (-) denote the sigmoid function and Plogia (0, x, a) be the network's output probability for input features and sensitive attribute a.
The neural network consists an input layer, two hidden layers with ReLU activation (sizes H1 32, H2 16), and final linear output layer combined with bias term by and term dependent the sensitive attribute:
h = ReLU + b,) E IR f11
h2 ReLU(W,h1 + b2) E R2
Plogit (0, xr w.3 h2 + aa + by E IR
where the parameters constituting 0 have dimensions: W1 E R Xisput.dim 6, R W2 H, by E H2 Wy E HR E R. and bn E R. So, the total dimension d H4 input_dim H1 + H2 H1 + H2 + Hz + 2 varies based on the input dimension, which itself depends on the feature hashing dimension used for categorical features. In particular, the parameter size of can be changed by adjusting hashing dimension.
The potential function is the negative fog-posterior (w) = (log P(D () + log P(O)). where D is the training data, ()) the log-likelitee ligeeiillhood using the sigmoid of the logits, and log P(0) the log-prior based isotropic Gaussian distribution with precision 1
61

---

The equality constraints (h(0) 0) enforce faireess via demographic parity on True Positive Rate (TPR) and False Positive Rate (FPR) between groups defined by the sensitive attribute (default: gender) a E {0, 1):
hTPR( (u) = 1...//= (Progi (v, 2;, 1))] = (o (Pogit (v, x, 0)) = 0
(p) = Be..//=ll.y=0l for (Pogi (v, 2;, 1))] [. (or (Progit (U, z; 0)) = 0.
These expectations are estimated using averages over the training data subsets corresponding each sensitive attribute and true label y.-
The inequality constraints (0()) 0) enforce monotonicity on selected features (default: duration, credit amount, existing credit, age) by requiring the gradient the logit with respect to these features have a specific sign (or be close to zero, within margin 8 1.0) at subset of anchor data points Danchor. Let S = {duration, credit amount, existing credits} and S = {age}. The constraints are formulated as:
{
Oprogia Die a:)
Oprogi (v, Di, a.)
})
9mono (t) = max
max jES+ E Eacchre
max KES 22.ED
dzij
The dimension d varies based on feature hashing as d E (706, 1986, 4994, 9986, 49986, 100002).
We single chain for K = 200 steps, discard the first 40 burn-in, and thin by 2. Also, the test Negative Log-Likelihood (NLL) as the metric.
Table 11: Hyperparameter settings for the Bayesian logistic regression task
Method
Hyperparameters
OLLA-H
a = 100, e = 1.0, E (0,5), 5 10
CLangevin
L 10, 1.0, À = 0.5, A 5 10-4
CHMC
7 =1..0, T. = 10, r 1.0, 0.5, At 5 x 10-3
CGHMC
7= 1.0.0, L. 10, 1.0, 0.5, AI 5 x 10-3
62


Reviews:
[{"review_id": "I73CDd7MlI", "reviewer": "Reviewer_tmLS", "summary": "This paper makes use of the landing technique -- incorporting constraints directly into the guiding stochastic differential equation (SDE) -- from non-convex optimization to handle both equality and inequality constraints in sampling. The proposed algorithm, Overdamped Langevin with LAnding (OLLA), avoids projection to the feasible set, and hence runs efficiently and is proven to converge exponentially fast. With Hutchinson trace estimation, to approximate the Hessian in the mean curvature term, OLLA could scale to high dimensions.", "questions": "Here are a few more questions and suggestions:\n\n* In the submanifold $\\Sigma$, the inequality constraint is $g(x)\\leq 0$. However, in the set of active inequality constraints, why we have $g_i(x)\\geq 0$?\n* line 97: $g_{I_x}$ $\\to$ $g_{|I_x|}$.\n* In Section 3, the constants $\\kappa$, $\\delta$ and tabular neighborhood were given very abruptly and there lacked some background.\n* line 123, $h(X_t)$ $\\to$ $h_i(X_t)$.", "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": "# Strengths:\n\n* The paper integrates the both equality and inequality constraints into Langevin dynamics and avoids projection.\n* The proposed OLLA algorithm comes with theoretic convergence guarantee.\n* OLLA-H algorithm scales to high dimensions and has SOTA performance compared with similar samplers with constraints.\n\n# Weaknesses:\n\n* It is a little disappointing that OLLA-H has comparable performance with CGHMC in the included tests. Hence it remains doubt whether practitioners choose one over the other.\n* In all the simulation tests, the highest dimensionality is 700. The paper could be strengthened by including an example of dimension bigger than 1000.\n* In all the simulation tests, the target distributions are relatively simple -- either uniform or Gaussian, though the constraints are complicated and versatile. This may not be very realistic --  in practice both distribution and constraints could be complicated. It would be more convincing to test OLLA(-H) on some more practical distributions with certain difficulty to sample.", "was_revised": true, "final_justification": "Thanks for the authors' rebuttal that resolved my concerns. I will keep my score.", "rebuttal": "We thank the reviewer for the thoughtful and constructive feedback. Below, we (1) give a brief overview of what we changed, and then (2) leave itemized responses to address each of the comments, and we hope they can clarify your questions and concerns.\n\n\n\n---\n**Overview of main changes/clarifications.**\n* We added two new experiments (Table A, B) to demonstrate the advantages of OLLA-H under practical scenarios: (A) a physics-motivated molecular system with a nontrivial potential and complex constraints, and (B) a high-dimensional Bayesian logistic regression with fairness and monotonicity constraints.\n* We clarify computational trade-offs between OLLA-H and CGHMC (and other projection-based baselines), including per-iteration complexity.\n* We clarify the notations of active inequality set $|I_x|$, and backgrounds on $\\kappa, \\delta$.\n---\n>Weakness (1): It is a little disappointing that OLLA-H has comparable performance with CGHMC in the included tests. Hence, there remains doubt whether practitioners choose one over the other.\n>\n>Weakness (2): In all the simulation tests, the highest dimensionality is 700. The paper could be strengthened by including an example of a dimension bigger than 1000.\n>\n>Weakness (3): In all the simulation tests, the target distributions are relatively simple -- either uniform or Gaussian, though the constraints are complicated and versatile. This may not be very realistic -- in practice, both distribution and constraints could be complicated. It would be more convincing to test OLLA(-H) on some more practical distributions with certain difficulty in sampling.\n\n**Answer.** OLLA-H is preferable to CGHMC (also other baselines) when (1) large portion of inequality constraints stays inactive (loose boundaries $\\Leftrightarrow$ small $|I_x|$), (2) constrained domain $\\Sigma$ and potential function $\\nabla f(x)$ are so complicated that they force small step size $\\Delta t$ for CGHMC to have reasonable acceptance rate on Metropolis-Hastings correction, or (3) Newton steps for projection become larger under geometrically complicated $\\Sigma$.\n\n**Computational picture (per step).**\n\nFollowing notations on paper, let $m$ be the number of equality constraints, $l$ the number of inequalities, and $I_x \\subset \\\\{1,...,l\\\\}$ the currently active inequality set. OLLA-H forms and inverts a Gram matrix $G$ using only equalities and active inequalities:\n* $\\text{cost(OLLA-H)} = \\mathcal O(N \\cdot (m+ |I_x|)^3)$\n\nwhere $N$ is the number of Hutchinson probes. For the slack-variable + Newton's methods (CLangevin / CHMC), they inflate the cubic of constraint block to $(m+l)$:\n* $\\text{cost(CLangevin, CHMC)} = \\mathcal O(N_{newton} \\cdot (m+l)^3)$\n\nwhere $N_{newton}$ is the number of required iterations for Newton's method to converge.\nThe CGHMC (Metropolis-Hastings + Newton's method) avoids slack variables, but samples are accepted or rejected based on Metropolis-Hastings correction:\n* $\\text{cost(CGHMC)} = \\mathcal O(N_{newton} \\cdot m^3), \\quad \\text{efficiency} \\propto \\text{acceptance rate}.$\n\nWhen the boundary of the feasible set $\\Sigma$ is tight and complex (many inequalities intermittently active), $|I_x|$ can grow as OLLA-H's cubic term reflects this. However, in the same regime, CGHMC's acceptance rate typically deteriorates due to frequent boundary interactions. Also, if the step size is too big to drastically change the energy of Hamiltonian of the system, the CGHMC's acceptance rate again drops, yielding low ESS/CPU-time. Simultaneously, projection based methods (CLangevin, CHMC, CGHMC) often require larger $N_{newton}$ and can be extremely slowed down if the Newton method stalls on irregular or high-dimensional $\\Sigma$. \n\n**Additional Empirical Results (Appendix Table A,B).**\n\nFor the experiment setup (potential, constraints), please refer to our reply to the reviewer 74bE.\n\n****Table A - Molecular system with realistic potential (complex $\\Sigma$, small $|I_x|$, large $N_{newton}$)****\n\nThe below figure shows the estimation value of the test function $\\mathbb E [|x|]^2$ with total CPU time to run (wrapped by brackets). The constraint violations for both equality constraints and inequality constraints are below 0.005 / 0.02, respectively, across all algorithms. Also, the $(n, m, l)$ configurations are given by (15, 7, 6), (30, 17, 36), (45, 27, 91), (60, 37, 171), (90, 57, 406) \n| Method \\ dim | 15 | 30 | 45 | 60 | 90 |\n| --- | --- | --- | --- | --- | --- |\n| OLLA-H | 23.02 (22.51s) | 236.6 (42.48s) | 852.9 (67.38s) | 2152 (107s) | 7561 (257.8s) |\n| CLangevin | 36.44 (294.4s) | 228.5 (1611s) | 925.5 (3724s) | 2162 (6631s) | 7572 (15230s) |\n| CHMC | 31.1 (79.34s) | 266.8 (280.1s) | 881.3 (614s) | 2235 (1090s) | 7640 (2514s) |\n| CGHMC | 29.8 (49.96s) | 256.8 (107.9s) | 895.7 (174.8s) | 2174 (234.3s) | 7607 (415.1s) |\n\n\nIn this setting, only a small subset of inequalities is active at each step (so $|I_x| \\ll l$), which makes OLLA-H's Gram matrix compact and economic to compute. Also, due to the complexity of the potential and constrained domain $\\Sigma$, the $N_{newton}$ remains relatively high $(\\approx 20)$. Across the dimensions, OLLA-H attains comparable estimates of test functions while achieving best CPU-time improvements of at least 2x to 60x relative to other baselines, further widening the gap in effectiveness. This highlights OLLA-H's advantage when $\\Sigma$ is intricate but sparsely active in inequalities.\n\n****Table B - Bayesian logistic regression with fairness + monotonicity constraints (High-dimensional $\\Sigma$)****\n\nThe below figure shows the test NLL loss on UCI german credit data with total CPU time to run (wrapped by brackets). All baselines kept equality constraint biolations at $\\leq 10^{-3}$ level. For inequality constraints, each algorithms achieve around $0.99$ (OLLA-H), $0.05$ (CLangevin), $0.02$ (CHMC), across all dimensions. \n\n| Method \\ dim | 706 | 1986 | 4994 | 9986 | 49986 |\n| --- | --- | --- | --- | --- | --- |\n| OLLA-H | 0.603 (8.25s) | 0.541 (8.65s) | 0.551 (8.30s) | 0.564 (9.85s) | 0.577 (15.85s) |\n| CLangevin | 0.604 (16.45s) | 0.546 (16.55s) | 0.550 (19.75s) | 0.544 (20.55s) | 0.557 (28.50s) |\n| CHMC | 0.605 (18.10s) | 0.602 (19.40s) | 0.599 (21.35s) | 0.602 (27.05s) | 0.601 (52.20s) |\n\n\nHere, the target density is a non-Gaussian posterior, and constraints are also nonlinear functionals, both implemented via a two-layer neural network. In this experiment, the CGHMC suffers from a low-acceptance rate problem ($\\le 5 \\\\% $) due to the complicated potential (or high step size) and frequent inequality activations. So, we dropped the result of CGHMC. While the results show the projection-based algorithms (CLangevin, CHMC) attain smaller inequality constraint violations over OLLA-H, the lower computational cost of OLLA-H is still highlighted with moderate constraint violations, which corresponds to our purpose of the algorithm. \n\n---\n\n>Question (1): In the submanifold $\\Sigma$, the inequality constraint is $g(x) \\leq 0$, However, in the set of active inequality constraints, why do we have $g_i(x) \\geq 0$?\n\n**Answer.** Thank you for pointing out. Our notation follows:\n* Feasible set: $\\Sigma:= \\\\{x \\in \\mathbb R^d  \\mid h(x) =0,\\  g(x) \\leq 0 \\\\}$\n* Set of active inequality constraints: $I_x := \\\\{i \\in [l] \\mid g_i(x) \\geq 0\\\\}$\n\nUnder this notation, $g_i(x) \\geq 0$ will be flagged as \"activated\", and we enforce them to decay exponentially fast via our OLLA dynamics. We understand our sign conventions for inequalities may not align with other literature. For example, [1] use feasibility $g_i(x) \\geq 0$ and define the active set as $I_x := \\\\{i \\in [l] \\mid g_i(x) \\leq 0\\\\}$. However, this is equivalent to ours under a sign change of inequalities.\n\n---\n\n>Suggestion (1): In section 3, the constants $\\kappa$, $\\delta$, and tubular neighborhood were given very abruptly, and there lacked some background.\n\n**Answer.** Thank you. We agree that these should be introduced earlier with some intuition. We will move the definitions to Preliminaries and add a short intuitions with cross references on the revision. Here, we leave concise definitions and roles for $\\kappa, \\delta$ and their backgrounds.\n* (Lemma C.1/C.3) $\\kappa :=  \\frac{1}{2} \\min_{y \\in \\Sigma} \\sigma_{min} ([\\nabla h(y)^T, \\nabla g_{I_y}(y)^T]^T ) >0$\n\nThis constant is the proportionality constant used when showing that reducing $h$ and $g$ decreases the normal displacement $||x-\\pi(x)||_2$. The existence of this constant is guaranteed by LICQ condition and determined by the manifold structure of $\\Sigma$.\n\n* (Theorem C.1/C.2) $\\delta$ is the reach of \"recoverable tubular neighborhood\" $\\hat{U}\\_\\delta (\\Sigma)$, which is defined by $\\hat{U}_\\delta (\\Sigma) :=\\\\{x \\in \\mathbb R^d \\mid ||h(x)||_2 <\\delta, \\ g(x) < \\delta \\\\}$.\n \nFor $x \\in \\hat{U}\\_\\delta (\\Sigma)$, the nearest-point projection $\\pi(x) \\in \\Sigma$ is well-defined and there exists a recovery map $\\zeta$ such that, writing $y =\\pi(x), p = h(x), q = g_{I_y}(x)$, we can recover $x$ by $x=\\zeta(y,p,q)$. This structure underlies our analysis of the projected process $Y_t = \\pi(X_t)$.\n\n---\n>Suggestion (2): Typos in Line 97 ($g_{I_x} \\rightarrow g_{|I_x|}$), Line 123 ($h(X_t) \\rightarrow h_i(X_t)$.\n>\n**Answer.** Thank you for pointing out. We'll fix both errata.\n\n---\nThank you again for your positive review and constructive feedback regarding the paper. We'd appreciate the opportunity to discuss further if you have any additional questions or suggestions.\n\n### Reference\n\n[1] M. Muehlebach and M. I. Jordan (2022). On constraints in first-order optimization: A view from non-smooth dynamical systems. JMLR 23(256)."}, {"review_id": "5fopNfg6NV", "reviewer": "Reviewer_cZ9Z", "summary": "This paper studies the problem of sampling from a density satisfying log-Sobolev inequality with nonlinear equality and inequality constraines, even when the feasible set is nonconvex. The authors propose a novel framework, OLLA, which modifies Langevin dynamics to enforce constraints via a landing mechanism. By operating on the tangent space of the constraint manifold, OLLA avoids costly projections. The method achieves exponential convergence under some assumptions and is efficiently discretized using Hutchinson trace estimation. Empirical results demonstrate that OLLA matches or exceeds the performance of existing constrained sampling algorithms in both accuracy and computational efficiency, particularly in high-dimensional settings.", "questions": "How should one select the landing rate and step size?\n\nSince there is no theoretical analsysis for OLLA-H, how is the Hutchinson trace estimator affect the convergence rate?\n\nCould you clarify the assumptions underlying the main theorems and illustrate concrete settings in which they hold?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\nThis paper propose a novel sampling framework combined with projection free landing mechanism in optimization and shows exponential convergence guarantees in 2-Wasserstein distance under mild LICQ and log-Sobolev assumptions. To make the algorithm more efficient they further modify the algorithm by using Hutchinson trace estimation, yielding low per-step cost even in high dimensions. \n\nWeakness:\n\nThe stated assumptions are unclear and may not hold in practical scenarios, and the algorithm is highly sensitive to its hyperparameters, yet the authors provide no guidance on how to select them.", "was_revised": true, "final_justification": "The authors' rebuttal addressed my concerns. I will keep my positive score.", "rebuttal": "We appreciate the reviewer for the constructive feedback. Below, we provide itemized responses to address each of the comments, and we hope they clarify your questions and concerns.\n\n---\n>Weakness (1): The stated assumptions are unclear and may not hold in practical scenarios.\n>Question (3): Could you clarify the assumptions underlying the main theorems and illustrate concrete settings in which they hold?\n\n**Answer.** Below we restate them concisely, provide settings where they hold, and for the mixed case, we propose a milder replacement (M1'') of the strong assumption (M1) together with a corollary.\n\n**Core assumptions (C1) - (C3)**\nFirst, we recall our notations $\\Sigma := \\\\{x \\in \\mathbb R^d \\mid h(x) =0 , \\ g(x) \\leq 0\\\\}$\n* Assumption (C1)-LICQ: For every $y \\in \\Sigma$, $\\\\{\\nabla h_i(y)\\\\}\\_{i \\in [m]}  \\cup \\\\{\\nabla g_j(y)\\\\}\\_{j \\in I_y}$ is linearly independent.\n\nThis assumption guarantees a uniformly well-conditioned normal space and underlies the projection/landing analysis.\n\n* Assumption (C2)-Regularity of $\\Sigma$: $\\Sigma$ is assumed to be compact and connected. Also, $\\nabla h(x) \\neq 0$ on $\\Sigma$ and $h$ is coercive ($||h(x)||_2 \\rightarrow \\infty$ as $||x||_2 \\rightarrow \\infty$). When there are only inequality constraints, we assume $\\text{dim}(\\Sigma) = d$.\n\nThis assumption are required for the existence of $\\lambda_{LSI}$ (from C4 in Appendix) and linkage between constraint violation and normal displacement: smaller $||h(x)||_2, \\ g(x)^+$ implies smaller $||x-\\pi(x)||_2$ within its recoverable tubular neighborhood. \n\n* Assumption (C3)-Bounded initialization: For the initial law $\\rho_0$, we assume $M_h := \\sup\\_{x_0 \\in \\text{supp}(\\rho_0)} ||h(x)||_2 < \\infty, M_g :=\\sup\\_{x_0 \\in \\text{supp}(\\rho_0)} ||g(x)||_2 < \\infty$.\n\nThis ensures a finite landing time (almost surely) into the recoverable tubular neighborhood (for equality constraints) or the constrained domain $\\Sigma$ (for inequality constraints). \n\n**Examples satisfying (C1)-(C3).** Some of classical manifolds can satisfy assumptions (C1)-(C3). For example, sphere $S^{d-1} := \\\\{x \\in \\mathbb R ^d \\mid ||x||_2 = 1\\\\}$, Stiefel manifold $St(n, p) := \\\\{X \\in \\mathbb R^{n \\times p} \\mid X^T X = I_p\\\\}$ (with $n > p$ ), Tori $T^k$.\n\n**Mixed-case only assumptions (M1) - (M3)**\nOur mixed-case analysis introduced the assumption (M1) to control how the projected manifold $\\Sigma_p$ intersects the interior of $\\Sigma$. \n* Assumption (M1)-Regularity of $\\Sigma_p$: The projected manifold $\\Sigma_p := \\pi ( \\\\{ x \\in \\mathbb R^d \\mid h(x) = p, g(x) \\leq 0 \\\\})$ lies inside the interior of $\\Sigma$ for $0 < ||p||_2 < \\delta$, where $\\delta$ is the width of recoverable tubular neighborhood $\\hat{U}\\_{\\delta}(\\Sigma)$.\n\nWe agree that (M1) can be strong in many practical geometries. To understand why this is assumed. We can think of the following examples. We set $h(x,y) = 0$ and define the equality-only projection $\\pi^h(x,y) = (x,0)$ (the nearest-point projection onto $\\\\{y=0\\\\}$). Also, we write $\\Sigma:= \\\\{(x,0) \\in \\mathbb R^2 \\mid g(x,0) \\leq 0\\\\}$, $\\Sigma_p := \\pi^h ( \\\\{(x,y) \\in \\mathbb R^2 \\mid y =p, g(x,y) \\leq 0\\\\})$.\n* Example (1)-satisfy (M1): pick $g_{circ}(x,y) = 2(x^2+y^2) - 1 \\leq 0$\n\nIn this case, $\\Sigma_p$ lies inside the $\\text{int}(\\Sigma)$ for $p \\in (-1, 1)$, and satisfy the assumption (M1). Also, we note that the projected process $Y_t = \\pi(X_t)$ fully lies on $\\text{int}(\\Sigma)$, not stuck at the boundary of $\\Sigma$ (Here, $\\pi$ is the neareast-point projection onto $\\Sigma$).\n\n* Example (2)- violates (M1): pick $g_{hyp}(x,y) = 2(x^2-y^2) -1 \\leq 0$\n\nIn this case, while $\\Sigma$ is identical to example (1), $\\Sigma_p := \\\\{(x,0) \\in \\mathbb R^2 \\mid |x| \\leq \\sqrt{\\frac{1}{2} +p^2}\\\\} \\not\\subset \\text{int}(\\Sigma)$ and fails to satisfy (M1). In this example, the projected process $Y_t$ can stuck at the boundary when $\\pi^h(X_t) \\not\\in \\text{int}(\\Sigma)$, slowing down the convergence of $Y_t$.\n\n* Example (3)- dramatric case: pick $g_{\\text{loose}}(x,y) = 2(x^2 - 1000y^2) -1$.\n\nThis inequality constraint can define the same $\\Sigma$ defined on Example (1) or (2). However, under $g_{\\text{loose}}$, the feasible region of this inequality stretches far along $\\Sigma$, so boundary hits persits and the process $\\pi^h(X_t)$ routinely overshoots $\\Sigma$. This behavior degrades the sampling efficiency when $X_t$ hovers relatively high on $\\Sigma$, even though $h(X_t)\\rightarrow 0$ at the end.\n\nTo incorporate this, we propose a milder assumption (M1'') below:\n* Assumption (M1''): Define $p_t$ = $P(\\pi(X_t) \\notin \\text{int}(\\Sigma))$. we assume that after sufficently large time $t_{bd}$, there exist $\\alpha_{bd}, C_{bd}, \\epsilon_{bd} >0$ such that $p_t \\leq C_{bd}e^{-\\alpha_{bd}t}$. \n* Corollary (informal): Under assumptions (C1)-(C3) and (M1'')-(M3), for $t \\geq t_{bd}$ and $\\alpha > 2\\lambda_{LSI}, \\beta \\geq 1$, it holds that \n\\begin{equation}\n    W_2(\\rho_t, \\rho_\\Sigma) \\leq \\frac{M_h}{\\kappa} e^{-\\alpha t} + \\sqrt{\\frac{2}{\\lambda_{LSI}} \\textsf{KL}^{\\Sigma} (\\tilde{\\rho}\\_{t}^{\\Sigma\\_{in}} || \\rho\\_{\\Sigma})  + C\\_{bd}e^{-\\alpha\\_{bd}t}\\cdot \\text{diam}\\_{\\Sigma} (\\Sigma)^2} \n\\end{equation}\nwhere $\\textsf{KL}^{\\Sigma}(\\tilde{\\rho}_{t}^{\\Sigma\\_{in}} || \\rho\\_{\\Sigma}) = \\mathcal O(\\exp(-2\\lambda\\_{LSI} t))$ as in Theorem 3.\n\nTo comment on $p_t$, if the inequality constraints behave nicely so that (M1) is achieved or the feasible set of inequality constraints contracts to $\\Sigma$ very fast, the $p_t$ term decays fast, leading to a better convergence rate. This aligns with our observation from previous examples. \n\nThe remaining assumptions (M2) and (M3) are comparably milder than (M1). For the assumption (M2), we can pick sufficiently high $V, \\beta$ accordingly to $\\Sigma$, because we only require $V > 0, \\beta>1$. Also, we left the comment of Assumption (M3) on Remark 5 (in Appendix A), and this is likely to be satisfied for large enough $t$ due to uniformly elliptic behavior of the associated Fokker-Planck equation.\n\n---\n>Question (2): Since there is no theoretical analysis for OLLA-H, how does the Hutchinson trace estimator affect the convergence rate?\n>\n**Answer.** The stochastic trace is used only in the mean-curvature term $\\mathcal H$. In our formulation, $\\mathcal H$ is multiplied by $\\nabla J^T G^{-1}$ (normal operator), hence the Hutchinson noise acts only in the normal direction and tangential mixing on $\\Sigma$ is minimally affected. What change is the second-order discretization bias in the landing dynamics, which we will show below:\n* (Decomposition) Let $Y_t := \\pi(X_t)$ and $Z_t$ be the target process ($Z_t \\sim \\rho_\\Sigma$). Then, the distance between $X_t$ and $Z_t$ is decomposed into $\\mathbb W_2(\\rho_t, \\rho_\\Sigma) \\leq E||X_t-Z_t||_2 \\leq E||X_t-\\pi(X_t)||_2 + E||Y_t-Z_t||_2$\n\nThe first term primarily depends on $\\alpha$ and the statistics of $\\mathcal H$, while the second term heavily depends on the on-manifold geometry (e.g., LSI constant) and should be minimally affected by the Hutchinson noise (from our continuous dynamic analysis in Theorem 3). Therefore, we analyze how the noise affects the first term.\n* Lemma (Informal, equality-only, $m=1$). Suppose $||X_k||_2$ is bounded for $k\\geq 1$ a.s. Then, we have\n\\begin{equation}\n    \\mathbb E [h(X_k)] \\leq (1-\\alpha \\Delta t)^K \\mathbb E[h(X_0)]  + \\mathcal{O} (\\frac{\\Delta t}{\\alpha} ( (1+ \\frac{d^2}{N})) \n\\end{equation}\nwhere $K$ is the total number of iterations. \n\nTherefore, the Hutchinson estimator may not affect the converence rate, but rather increase the discretization error, which again vanishes as $N \\rightarrow \\infty$.\n\n---\n>Weakness (2): The algorithm is highly sensitive to its hyperparameters, yet the authors provide no guidance on how to select them.\n>Question (1): How should one select the landing rate and step size?\n>\n**Answer.**  Below we give a principled rule and a practical recipe.\n\n**Landing rate $\\alpha$.** First, we observe that a monotonically decreasing upper bound of $W_2$ can be achieved whenever $\\alpha > 2\\lambda_{LSI}$ (from Theorem 3). However, $\\lambda_{LSI}$ is unknown in most cases. So, we rely on the following informal theorem (Theorem 7.3, [2]):\n* Theorem (Informal) Let $\\Sigma$ be a compact Riemannian manifold with diameter $D = \\sup_{x,y \\in \\Sigma} d_\\Sigma(x,y)$ and non-negative Ricci curvature. Then, $\\lambda_{LSI} \\geq \\frac{\\lambda_1}{1+2D\\sqrt{\\lambda_1}}$, or $\\lambda_{LSI} \\geq \\frac{\\pi^2}{(1+2\\pi)D^2}$ holds, where $\\lambda_1$ is the first eigenvalue of the Laplace-Beltrami operator.\n\nTherefore, we start searching a reasonable $\\alpha$ at a scale of $\\mathcal O(1/D^2)$, which not only decreases equality constraints fast but also maintains numerical stability of discretization.\n\n**Step size $\\Delta t$.** For choosing a reasonable step size, we note the previous analysis on the decay of $h(X_k)$. Because $\\mathbb E [h(X_k)]$ decays exponentially fast with base $(1-\\alpha \\Delta t)$, we choose a reasonable $\\Delta t > 0$ such that $\\alpha \\Delta t < 1$ with some step size searching.\n\n**Boundary repulsion rate $\\epsilon$.** These terms play a crucial role in satisfying the inequality constraints. Too small $\\epsilon$ will make trajectories sticky near the boundary of $\\Sigma$, and too large $\\epsilon$ will make aggressive repulsion on the boundary so that the sample overshoots.\nIn practice, we monitor the inequality violation during sampling and adjust $\\epsilon$ so that the sample does not overshoot and show sticky behavior at the boundary. Fortunately, it turns out that $\\epsilon$ is less sensitive than $\\alpha$ (Table 6 in Appendix), and it is often not hard to choose a good $\\epsilon$.\n\n---\nThank you again for your positive review and helpful feedback. We'd appreciate the opportunity to discuss further if you have any additional questions or suggestions.\n\n### Reference\n[2] M. Ledoux (1999). Concentration of measure and logarithmic Sobolev inequalities. Séminaire de Probabilités 33."}, {"review_id": "n0w38rMyka", "reviewer": "Reviewer_74bE", "summary": "This paper proposes OLLA, a projection-free Langevin dynamics method for sampling under nonconvex equality/inequality constraints. By combining tangential noise with normal drift (\"landing\"), OLLA enforces constraints without projections or slack variables. Theoretical analysis shows exponential convergence in Wasserstein distance, while a practical variant (OLLA-H) uses trace estimation for scalability. Experiments demonstrate superior efficiency over baselines in high dimensions.", "questions": "* How should practitioners choose $(\\alpha, \\epsilon)$ for real-world problems? Are there heuristics or adaptive schemes?\n* How does OLLA behave when constraints violate LICQ (e.g., degenerate gradients)?\n* The paper notes degraded performance for large $m$. Is there a theoretical or empirical threshold beyond which OLLA becomes impractical?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n* This paper introduces a novel SDE design with a landing term for constraint enforcement, combined with rigorous non-asymptotic convergence analysis for three constraint scenarios.\n* This paper addresses a critical gap in constrained sampling, especially for nonconvex sets where projections fail, which also has practical impact for Bayesian inference, generative modeling, and other ML applications with constraints.\n* This is the first method to unify equality/inequality constraints in Langevin dynamics without projections.\n\nWeaknesses:\n* Critical assumptions are deferred to the appendix, making it harder to evaluate the method’s applicability upfront. This obscures potential limitations, such as sensitivity to constraint degeneracy or nonsmoothness, which should be highlighted in the main text.\n* Hyperparameters $(\\alpha, \\epsilon)$ are sensitive but lack principled selection guidelines.\n* Limited to synthetic experiments; real-world benchmarks would strengthen the impact.", "was_revised": true, "final_justification": "The authors answered my questions and addressed my concern, and thus I would like to keep my positive score.", "rebuttal": "We appreciate the reviewer for the constructive feedback. Below, we provide itemized responses to address each of the comments, and we hope they clarify your questions and concerns.\n\n---\n>Weakness (3): Limited to synthetic experiments; real-world benchmarks would strengthen the impact.\n\n**Answers.** We added two real-world tests that demonstrate the advantages of OLLA-H under complicated potentials and domain constraints. For the experimental results, please refer to the results we provided to reviewer tmLS.\n\n**Experiment (A): Molecular/Polymer system.**  In this experiment, we have $N$ atoms $P = (P_0, ..., P_{N-1}) \\in (\\mathbb R^3)^N$, $(d = 3N)$ and we sample from the complicated potential $f(P) = \\sum_{i <j, |i-j|>1} \\phi_{nb}(P_i, P_j) + \\sum_{k=1}^{N-3} \\phi_{tors}(P_{i-1}, P_i, P_{i+1}, P_{i+2})$, where $\\phi_{nb}$ is 12-degree piece-wise polynomials (e.g. nonbonded pair potential), $\\phi_{tors}$ is highly nonlinear function involving periodic cosine function (e.g. Periodic dihedral torsion potential). \n\nThe equality constraints are defined by $h_k^{bond}(P) = ||P_k-P_{k+1}||_2^2 - l^2 = 0, \\ k =0, ..., N-2$, and $h_k^{ang}(P) = (P\\_{k-1} - P\\_{k}) \\cdot (P\\_{k+1} - P\\_{k}) - l^2 \\cos(\\theta), k=1,...,N-2$, and the inequality constraints are given by $g\\_{ij}(P) = r\\_{min}^2 - ||P_i-P_j||_2^2$ for $|i-j| > 1$.\n\n**Experiment (B): UCI German Credit: Bayesian 2-layer logistic regression with fairness + monotonicity constraints.** In this experiment, we use a two-layer ReLU neural network with hashed categoricals (to incorporate categorical data). The potential function $f(\\theta)$ is the posterior distribution with an isotropic Gaussian prior.\n\nThe equality constraints are $h_{TPR}(\\theta) = \\mathbb E [p_\\theta(X,A) \\mid Y=1, A=1] - \\mathbb E [p_\\theta(X,A) \\mid Y=1, A=0]$, $h_{FPR}(\\theta) = \\mathbb E [p_\\theta (X,A) \\mid Y=0, A=1] - \\mathbb E [p_\\theta(X,A) \\mid Y=0, A=0]$, where $Y$ is the binary outcome label (e.g. good/bad credit), $A$ is the sensitive attribute (e.g. gender), and $p_\\theta (X,A)$ is the predicted probability based on $X,A$.\n\nThe inequality constraints are set to be $g_+(\\theta) = \\max_{r,j \\in S_+} \\\\{ -\\partial_{x_j} l_\\theta (z_r)\\\\}$, $g_-(\\theta) = \\max_{r,j \\in S_-} \\\\{\\partial_{x_j} l_\\theta (z_r)\\\\}$, where $S_+$ is the index set of monotone-increasing features, $S_-$ is the index set of monotone-decreasing features, $l_\\theta(z)$ is logit value based on $z = (x,a)$, and $z$ is anchor points to check monotonicity.\n\n**Summarized Results.** In scenarios where the set of active inequality constraints is sparse (i.e.$|I_x|$ is small) or the manifold geometry is especially intricate so that each projection requires many Newton iterations, OLLA-H proves more computationally efficient than other baseline algorithms.\n\n---\n\n>Weakness (2) & Question (1): Hyperparameters $(\\alpha, \\epsilon)$ are sensitive but lack principled selection guidelines. How should practitioners choose for real-world problems? Are there heuristics or adaptive schemes?\n>Question (3): The paper notes degraded performance for large $m$. Is there a \n>theoretical or empirical threshold beyond which OLLA becomes impractical?\n\n\n**Answers.** Below we give a principled rule and a practical recipe.\n\n**Landing rate $\\alpha$.** First, we observe that a monotonically decreasing upper bound of $W_2$ can be achieved whenever $\\alpha > 2\\lambda_{LSI}$ (from theorem 3). However, $\\lambda_{LSI}$ is unknown in most cases. So, we rely on the following informal theorem (Theorem 7.3., [2]):\n* Theorem (Informal) Let $\\Sigma$ be a compact Riemannian manifold with diameter $D = \\sup_{x,y \\in \\Sigma} d_\\Sigma(x,y)$ and non-negative Ricci curvature. Then, $\\lambda_{LSI} \\geq \\frac{\\lambda_1}{1+2D\\sqrt{\\lambda_1}}$, or $\\lambda_{LSI} \\geq \\frac{\\pi^2}{(1+2\\pi)D^2}$ holds, where $\\lambda_1$ is the first eigenvalue of the Laplace-Beltrami operator.\n\nTherefore, we start searching reasonable $\\alpha$ at a scale of $\\mathcal O(1/D^2)$, which not only decreases equality constraints fast but also maintains numerical stability of discretization.\n\n**Step size $\\Delta t$.** For choosing a reasonable step size, we note the previous analysis on the decay of $h(X_k)$. Because $\\mathbb E [h(X_k)]$ decays exponentially fast with base $(1-\\alpha \\Delta t)$, we choose a reasonable $\\Delta t > 0$ such that $\\alpha \\Delta t < 1$ with some step size searching.\n\n**Boundary repulsion rate $\\epsilon$.** This term plays a crucial role in satisfying the inequality constraints. Too small $\\epsilon$ will make trajectories sticky near the boundary of $\\Sigma$, and too large $\\epsilon$ will make aggressive repulsion on the boundary so that the sample overshoots.\nIn practice, we monitor the inequality violation during sampling and adjust $\\epsilon$ so that the sample does not overshoot and show sticky behavior at the boundary. Fortunately, it turns out that $\\epsilon$ is less sensitive than $\\alpha$ (Table 6 in Appendix), and it is often not hard to choose a good $\\epsilon$.\n\n**Adaptive method for choosing $\\alpha.$** \nIn our high-dimension stress test, we observed that keeping $\\alpha$ fixed while increasing the number of equalities $m$ slows the convergence rate (Figure 12, Appendix H). In contrast to this, this phenomenon was improved by adaptively choosing $\\alpha$; increasing $\\alpha$ with $m$ mitigates this problem (described in Section 5, Appendix H.2). A plausible explantation is that, as $m$ grows, the feasible manifolds's diameter $D$ shrinks, so the requirement $\\alpha > 2\\lambda_{LSI}$ effectively calls for a larger $\\alpha$.\n\n**Threshold hold of $m$ beyond which OLLA-H becomes ineffective.** \nHowever, when $\\alpha$ is increased, we also need to reduce the step size to keep $\\alpha \\Delta t <1$ to ensure discretization stability. If we combine the previous two criteria to choose $\\alpha$, we obtain a two-sided bound $\\mathcal 1/D^2 \\lesssim \\alpha \\leq 1/\\Delta t$. \n\nThis implies that if the lower bound $D^{-2}$ grows with $m$ (as observed in our stress test), then a fixed $\\Delta t$ will eventually violate $\\alpha \\Delta t < 1$. Consequently, $\\Delta t$ must be reduced as $m$ increases. At that point, to cover the same physical time horizon $T$, the step count $K = T/\\Delta t$ increases, starting to make OLLA-H ineffective compared to other baselines.\n\n---\n>Question (2): How does OLLA behave when constraints violate LICQ (e.g., degenerate gradients)? \n\n**Answers.** LICQ is used in the theory to guarantee a uniformly well-conditioned normal space (also invertibility of Gram matrix $G$). In practice, OLLA/OLLA-H remains well-defined and typically stable under LICQ failures (e.g., degenerated/colinear constraints) by using a Moore-Penrose pseudoinverse $G^\\dagger$ (Remark 4, Appendix B). \n\nTo understand how this can mitigate the problem, we recall the induced definition of orthogonal projector $\\Pi(x):= I-\\nabla J(x)^T G(x)^\\dagger \\nabla J(x)$, which again becomes a idempotent projector onto $\\text{Null}(\\nabla J(x)) = T_x\\Sigma$. Hence, $\\Pi(x)$ always projects onto the tangent space $T_x\\Sigma$.\n\n* Example (redundant gradient): Let $h(x,y) = x^2+y^2 - 1$, $g(x,y) = x^2 + y^2 -\\frac{3}{2}$.\nIn this example, $\\nabla h, \\nabla g$ are colinear, making $G$ singular. However, $\\Pi(x):= I-\\nabla J(x)^T G(x)^\\dagger \\nabla J(x)$ still projects onto the circle's tangent line, which is precisely the intended behavior of our algorithm.\n---\n\n>Weakness (1): Critical assumptions are deferred to the appendix, making it harder to evaluate the method’s applicability upfront. This obscures potential limitations, such as sensitivity to constraint degeneracy or nonsmoothness, which should be highlighted in the main text.\n\n**Answers.** Thank you for pointing this out. This point was also raised by reviewer cZ9Z. We will move the key assumptions to the main text with short intuitions and reasonable examples to understand them. The following response briefly highlights the major changes; for detailed assumptions and clarifications, please see our reply to reviewer cZ9Z.\n\n**Reason to introduce new assumption (M1'').** One of our original mixed-case assumption (M1) asked the projected manifold $\\Sigma_p := \\pi( \\\\{h(x) =p , g(x) \\leq 0\\\\})$ to lie in $\\text{int} (\\Sigma)$ for small enough $||p||_2>0$, which can be strong in practice. To mitigate this, we propose the following assumption (M1'')\n* Assumption (M1''): Define $p_t$ = $P(\\pi(X_t) \\notin \\text{int}(\\Sigma))$. we assume that after sufficently large time $t_{bd}$, there exist $\\alpha_{bd}, C_{bd}, \\epsilon_{bd} >0$ such that $p_t \\leq C_{bd}e^{-\\alpha_{bd}t}$. \n* Corollary (informal): Under assumptions (C1)-(C3) and (M1'')-(M3), for $t \\geq t_{bd}$ and $\\alpha > 2\\lambda_{LSI}, \\beta \\geq 1$, it holds that \n\\begin{equation}\n    W_2(\\rho_t, \\rho_\\Sigma) \\leq \\frac{M_h}{\\kappa} e^{-\\alpha t} + \\sqrt{\\frac{2}{\\lambda_{LSI}} \\textsf{KL}^{\\Sigma} (\\tilde{\\rho}\\_{t}^{\\Sigma\\_{in}} || \\rho\\_{\\Sigma})  + C\\_{bd}e^{-\\alpha\\_{bd}t}\\cdot \\text{diam}\\_{\\Sigma} (\\Sigma)^2} \n\\end{equation}\nwhere $\\textsf{KL}^{\\Sigma}(\\tilde{\\rho}_{t}^{\\Sigma\\_{in}} || \\rho\\_{\\Sigma}) = \\mathcal O(\\exp(-2\\lambda\\_{LSI} t))$ as in Theorem 3.\n\nThis implies when inequalities and initial distributions are nice (e.g. (M1) holds or $\\Sigma_p$ contracts to $\\Sigma$ quickly), the last term (including $\\text{diam}\\_{\\Sigma} (\\Sigma)$) decays rapidy and the overall rate can be improved; otherwise that term becomes the bottleneck.\n\n---\nThank you again for your positive review and helpful feedback. We'd appreciate the opportunity to discuss further if you have any additional questions or suggestions.\n\n\n### Reference\n[2] M. Ledoux (1999). Concentration of measure and logarithmic Sobolev inequalities. Séminaire de Probabilités 33."}, {"review_id": "j9miO8sqW2", "reviewer": "Reviewer_AMyy", "summary": "This paper introduces Overdamped Langevin with Landing (OLLA), a novel projection-free sampling algorithm for distributions defined by complex, and potentially non-convex, equality and inequality constraints. By decomposing Langevin dynamics into tangential and normal components and enforcing constraint satisfaction via an exponential \"landing\" drift, OLLA provably converges exponentially fast in 2-Wasserstein distance under a log-Sobolev inequality (LSI) assumption. The authors also present a practical Euler–Maruyama discretization (OLLA-H) using a Hutchinson estimator for curvature corrections, achieving only $O(N)$ gradient costs per step.", "questions": "Do the analysis and convergence guarantees extend naturally to high dimensions? If not, could the authors discuss the challenges in making this dependence explicit?", "limitations": "There are no limitations or potential negative societal impacts of their work.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n* The authors provide a unified framework to handle arbitrary smooth equality and inequality constraints without explicit projections.\n* The authors provide rigorous non-asymptotic exponential convergence guarantees in the 2-Wasserstein distance, leveraging LSI constants.\n* Practical discretization via Hutchinson trace estimation reduces per-step cost to $O(N)$ gradient evaluations.\n* The paper is generally well-structured.\n\n**Weaknesses**\n* The theoretical convergence rate has a similar mathematical form to that of unconstrained Langevin algorithms under good conditions, depending on an LSI constant. Intuitively, one might expect that constraining the sampling domain would improve the LSI constant and lead to faster convergence, especially as the domain diameter shrinks [1-3]. The paper does not provide such an analysis or an estimate of the LSI constant in terms of the domain geometry. It is unclear how difficult such an analysis would be. Could authors comment on that, and if possible, offer a theoretical advantage in convergence rate?\n* The experimental evaluation is confined to synthetic tasks. As it is a paper with a strong theoretical focus, this is understandable.\n* Typos: line 210: \"Metropolis-Hasting correction\" -> \"Metropolis-Hastings correction\"\n\n**References**\n\n[1] Reflection Couplings and Contraction Rates for Diffusions, Probability Theory and Related Fields, 2016\n\n[2] Sampling from a Log-Concave Distribution with Projected Langevin Monte Carlo, Discrete & Computational Geometry, 2018\n\n[3] Projected Stochastic Gradient Langevin Algorithms for Constrained Sampling and Non-Convex Learning, COLT, 2021", "was_revised": true, "final_justification": "My concerns have been adequately addressed, and I have updated my score accordingly. The additional theoretical insights and experimental results provided in the rebuttal helped clarify my key concerns. I encourage the authors to include these additions in the revised version to strengthen the final paper.", "rebuttal": "We appreciate the reviewer for the constructive feedback. Below, we provide itemized responses to address each of the comments, and we hope they clarify your questions and concerns.\n\n---\n>Weakness (1): The theoretical convergence rate has a similar mathematical form to that of unconstrained Langevin algorithms under good conditions, depending on an LSI constant. Intuitively, one might expect that constraining the sampling domain would improve the LSI constant and lead to faster convergence, especially as the domain diameter shrinks [1-3]. The paper does not provide such an analysis or an estimate of the LSI constant in terms of the domain geometry. It is unclear how difficult such an analysis would be. Could authors comment on that, and if possible, offer a theoretical advantage in convergence rate?\n\n**Answer.** Thank you for providing insightful feedback. We agree that we missed this analysis in our paper. Below we (1) state a geometric lower bound for the LSI constant on $\\Sigma$, (2) explain how shrinking diameter improves the on-manifold mixing rate and how this interacts with the landing rate $\\alpha$ in our bound.\n\n**Geometric control of $\\lambda\\_{\\text{LSI}}.$** First, we bring up the following informal theorem (Theorem 7.3., [2]):\n* Theorem (Informal) Let $\\Sigma$ be a compact Riemannian manifold with diameter $D = \\sup_{x,y \\in \\Sigma} d_\\Sigma(x,y)$ and non-negative Ricci curvature. Then, $\\lambda_{LSI} \\geq \\frac{\\lambda_1}{1+2D\\sqrt{\\lambda_1}}$, or $\\lambda_{LSI} \\geq \\frac{\\pi^2}{(1+2\\pi)D^2}$ holds, where $\\lambda_1$ is the first eigenvalue of the Laplace-Beltrami operator.\n\nThis implies if the constraints shrink the diameter $D$ (or increase $\\lambda_1$), the lower bound of $\\lambda\\_{\\text{LSI}}$ increases, so the on-manifold decay of $\\textsf{KL}^\\Sigma(\\rho_t || \\rho_\\Sigma)$ acceraltes (dominated by exponential rate $-2\\lambda\\_{\\text{LSI}})$  \n\n**How this yields a theoretical advantage in our bound.** Our non-asymptotic bound reads, for sufficently large enough $t \\geq t_{cut}$,  \n\\begin{equation}\n    W_2(\\rho_t, \\rho_\\Sigma) \\leq \\frac{M_h}{\\kappa} e^{-\\alpha t} + \\sqrt{\\frac{2}{\\lambda\\_{\\text{LSI}}} \\textsf{KL}^{\\Sigma} (\\tilde{\\rho}_t || \\rho\\_{\\Sigma}) } \n\\end{equation}\nWith $\\textsf{KL}^{\\Sigma} (\\tilde{\\rho}_t || \\rho\\_{\\Sigma}) = \\mathcal O(e^{-2(\\lambda\\_{\\text{LSI}} (t-t\\_{cut}))})$, the second term decays at a scale $\\mathcal O(e^{-\\lambda\\_{\\text{LSI}} t})$. Thus the effective rate is controlled by $\\min \\\\{\\alpha, \\lambda\\_{\\text{LSI}} \\\\}$. By this analysis, we can observe that stronger constraints that reduce $D$ make $\\lambda\\_{\\text{LSI}}$ larger, so, provided $\\alpha$ is not the bottleneck, the overall convergence rate in $W_2$ improves.\n\n---\n>Question (1): Do the analysis and convergence guarantees extend naturally to high dimensions? If not, could the authors discuss the challenges in making this dependence explicit?\n\n**Answer.** Our convergence guarantees have some dependency on dimension $d$. The rates $\\alpha, \\lambda\\_{\\text{LSI}}$ are dimension-free, but the explicit $d$-dependence enters through the additive constant, which we can make explicit and a little sharper than now.\n\n**Where the dimension dependency enters.** To simplify the setup, we illustrate this using Theorem 1 (Same thing can be applied to mixed scenario). In theorem 1, there exists a constant $C_7:= \\frac{C_6 e^{-\\alpha t_{cut}}}{\\alpha - 2\\lambda\\_{\\text{LSI}}}$ (there is an errata; missing $C_6$), and $C_6$ depends on the constant $C_{div}$ by $C_6 = \\mathcal O (C_{div}^2)$, which again depends on $\\mathcal O (\\sqrt{d(d-m)})$ by the proof of Lemma E.5 in appendix. However, we revised the proof by using better linear algebra techniques and were able to conclude that we can reduce the dependence of $C_{div}$ by $\\mathcal O (\\sqrt{d-m})$. \n\nTherefore, we can say our convergence rate depends on $\\mathcal O (\\sqrt{d-m})$ (due to the square-root on $\\textsf{KL}^\\Sigma$) in the equality-only case (same analysis can be applied to mixed cases). We will incorporate this refinement and reflect this dependency in our paper.\n\n---\n>Weakness (2): The experimental evaluation is confined to synthetic tasks.\n\n**Answer.** We appreciate this concern. For the related answers, please refer to our responses to tmLS and 74bE, where we added two real-world benchmarks with full setup and results.\n\n---\n>Suggestion (1): Typos in Line 210: \"Metropolis-Hasting correction\" $\\rightarrow$ \"Metropolis-Hastings correction\".\n\n**Answer.**  Thank you for pointing this out. We will fix this.\n\n---\nThank you again for your positive review and helpful feedback. We'd appreciate the opportunity to discuss further if you have any additional questions or suggestions.\n\n\n### Reference\n[2] M. Ledoux (1999). Concentration of measure and logarithmic Sobolev inequalities. Séminaire de Probabilités 33."}]

Author rebuttal:

Meta-review: The paper proposes a new and unified projection free sampling method based on Langevin dynamics that can handle both equality and inequality constraints. The paper also shows non-asymptotic convergence analysis. Given the strong results and potential applications, I recommend it to be accepted.
Decision: Accept (poster)