Conceptio › Archive › arXiv CS

arXiv CSopen access

Manifold Sampling via Entropy Maximization

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Manifold Sampling via Entropy Maximization

arXiv:2605.12338v1 [cs.LG] 12 May 2026

Cornelius V. Braun1∗ Tilman Burghoff 1∗ Marc Toussaint1,2 1 Technische Universität Berlin 2 Robotics Institute Germany {braun,t.burghoff,toussaint}@tu-berlin.de

Abstract Sampling from constrained distributions has a wide range of applications, including in Bayesian optimization and robotics. Prior work establishes convergence and feasibility guarantees for constrained sampling, but assumes that the feasible set is connected. However, in practice, the feasible set often decomposes into multiple disconnected components, which makes efficient sampling under constraints challenging. In this paper, we propose MAnifold Sampling via Entropy Maximization (MASEM) for sampling on a manifold with an unknown number of disconnected components, implicitly defined by smooth equality and inequality constraints. The presented method uses a resampling scheme to maximize the entropy of the empirical distribution based on k-nearest neighbor density estimation. We show that, in the mean field, MASEM decreases the KL-divergence between the empirical distribution and the maximum-entropy target exponentially in the number of resampling steps. We instantiate MASEM with multiple local samplers and demonstrate its versatility and efficiency on synthetic and robotics-based benchmarks. MASEM enables fast and scalable mixing across a range of constrained sampling problems, improving over alternatives by an order of magnitude in Sinkhorn distance with competitive runtime.

Introduction

Sampling from constrained distributions is a fundamental problem in machine learning, with applications including Bayesian inference [1] and molecular design [2] as well as robotics and trajectory optimization [3]. In particular, many tasks in robotics require sampling kinematically feasible states or trajectories, for instance, to generate data for behavior cloning or to sample reset states for reinforcement learning [4–7]. The feasible set is usually given implicitly by constraint functions, which can only be evaluated point-wise. This makes sampling from it a challenging problem, as we have no prior knowledge about many of its properties, like the number of connected components. When no constraints are present, the dominant sampling paradigm is built around Markov Chain Monte Carlo (MCMC) methods such as Metropolis–Hastings, Langevin dynamics, and Hamiltonian Monte Carlo [8]. To handle constraints, these approaches are typically combined with projections or landing mechanisms that ensure samples remain feasible [1, 9–12]. However, all such methods rely on local MCMC kernels, which limits them to moves within a single connected region of the feasible set, and as a result they struggle with infeasibility barriers. As a result, even if such kernels mix well locally, they cannot determine how much probability mass to assign to each component (as illustrated in Figure 1). This makes them ill-suited in many practical scenarios, where such disconnected feasible sets arise naturally [13–15]. We tackle this problem by maximizing sample entropy across the feasible set. This provides empirical coverage of the feasible set, as the samples are distributed approximately uniformly. It enables estimation of component volumes and may serve as unbiased initialization for sampling from ∗ Equal contribution; Authors in alphabetical order.

Preprint.

OLLA

SCMC

MASEM-NHR (Ours)

Ground Truth

Disconnected Disks

Sine

NHR

Figure 1: Scatter plots of 100 samples from different samplers on synthetic benchmarks. Black lines show equality constraints, and green shaded areas mark the feasible region(s) by inequality constraints. The bottom row depicts 3d samples on spheres projected into 2d. NHR and OLLA fail to allocate correct sample masses: the larger disk and the leftmost peak of the sine are undersampled. arbitrary distributions across components. The key idea is to combine local k-nearest-neighbor density estimates [16] with Sequential Monte Carlo (SMC)-inspired resampling: particles in underrepresented components receive higher importance weights, encouraging redistribution of mass toward a uniform allocation. Crucially, the resampling step is modular and can be composed with generic constrained MCMC samplers. Adding this method to a constrained sampler is minimally invasive in the sense that it only incurs a small runtime overhead and does not degrade performance on compact and connected manifolds. Our main contributions are as follows: • We introduce Manifold Sampling via Entropy Maximization (MASEM), an algorithm for uniform sampling on manifolds with disconnected components. • We show that, under mild assumptions, the induced resampling operator contracts the KLdivergence to the uniform distribution at a geometric rate of (1 − τ /p) per iteration, where τ is a temperature parameter and p is the intrinsic dimension of the manifold. • We validate the approach empirically on synthetic and robotics benchmarks and compare performance across different local samplers.

Related Work

Sampling under constraints is typically approached by adapting classical MCMC methods to ensure feasibility, for example by projection, reflection, or annealing techniques. Soft-constraint and annealing approaches such as sequentially constrained Monte Carlo (SCMC) [17] improve global exploration by relaxing the constraints, but may produce infeasible samples (see Figure 1) and can still collapse to local constraint minima as the penalty is tightened. In contrast, our method enables mixing across infeasibility barriers and exact constraint feasibility at convergence. For manifold-constrained targets, constrained MCMC [10, 18–21], barrier-based [22] and projected Hamiltonian methods [11], landing-based approaches [9, 1], constrained SVGD [1, 23] and nonlinear Hit-and-Run (NHR) [24] operate directly on or near the feasible set and can maintain exact or approximate feasibility. Despite their differences, these methods share the common structural limitation that they rely on local transition kernels and therefore cannot reliably allocate probability mass across disconnected components. Our approach is complementary to this line of work: rather than introducing a new local kernel, it adds a resampling mechanism that can be combined with existing constrained samplers.

Preliminaries & Notation

We consider sampling uniformly from a feasible set Σ := {x ∈ Rd | h(x) = 0, g(x) ≤ 0} ⊆ Rd , 2

(1)

where h : Rd → Rm and g : Rd → Rℓ are differentiable equality and inequality constraint functions, which can only be evaluated point-wise and have Lipschitz-continuous derivatives. We assume that the feasible set Σ is bounded, for example by introducing the additional constraint ∥x∥ ≤ R for some large R > 0. Under this assumption, the feasible set and its connected components are compact. Furthermore, we impose the widely used Linear Independent Constraint Qualifier (LICQ). Assumption 1 (Linear Independent Constraint Qualifier [25]). For each point x ∈ Σ, the set of gradients of active constraints {∇hj (x)} ∪ {∇gi (x) | gi (x) = 0} is linearly independent. This assumption ensures that the feasible set Σ is a Riemannian manifold with corners and intrinsic dimension p = d − m [26, Lemma 3.1.12]. Let µΣ denote the induced Hausdorff measure and S = µΣ (Σ) its volume. We consider the general setting that the manifold decomposes into a finite number of disconnected components Σ = Σ1 ⊔ · · · ⊔ ΣC . Importantly, however, we do not assume knowledge of the number C of components or any of their structural properties, as they are implicitly defined via h and g. Due to LICQ, each component has measure σc = µΣ (Σc ) > 0. Note that if p = d, this problem reduces to sampling from the unnormalized density ρ̂(x) = 1Σ (x), where 1Σ (x) is the indicator function of Σ. Our target is to maximize the entropy of samples across the constrained set. Lemma 1 (Informally). The uniform density is the unique maximizer of the entropy on Σ. We formalize this in Appendix A. Our method maximizes the entropy by estimating the local density by k-nearest neighbor distances [16]. Let samples x1 , . . . , xN ∈ Σ be samples drawn from an underlying distribution parameterized by a density ρ. We define εi,k as the distance from xi to its k-th nearest neighbor among {xj }j̸=i . The k-nearest-neighbor density estimator at xi is ρ̂(xi ) =

k , N Vp εpi,k

(2)

where Vp is the volume of the unit ball in Rp . Consistency of this estimator, meaning ρ̂(xi ) → ρ(xi ) in probability as N → ∞ with k → ∞, k/N → 0, is a classical result for interior points of the support [27].

Manifold Sampling via Entropy Maximization

We now introduce Manifold Sampling via Entropy Maximization (MASEM), a framework for uniform sampling on a constrained set Σ as defined in Section 3. Our goal is to construct a particle system whose empirical distribution approximates the uniform density uΣ on Σ. Given an existing strong local sampler that preserves feasibility, the main challenge for constrained uniform sampling is global mass allocation. As Σ is only defined implicitly via constraints, this is a hard problem, since the measures of components or even their number are not known a priori, so local samplers cannot calibrate probability mass across disconnected components efficiently. Thus, MASEM introduces a mechanism that redistributes particles across components to ensure that each component receives mass proportional to its measure. The design of MASEM is intentionally lightweight. Its key mechanism is an importance-resampling step [8] that moves particles towards the maximum entropy target. This target is the uniform density p = uΣ (Lemma 1). Let q denote the current particle density on Σ. Since the target density is uniform, importance-resampling with weights wi ∝ p(xi )/q(xi ) reduces to resampling according to wi ∝ q(x)−1 . While the current density itself is unknown, we can estimate it using the k-nearest neighbor estimator (2). This leads to q(xi ) ∝ ε−p i,k in probability in the large particle limit, where εi,k denotes the distance of sample xi to its kth neighbor. Replacing the potentially unknown dimension p by a temperature parameter τ yields τ w̄i = εi,k ,

w̄i wi = P . j w̄j 3

(3)

Algorithm 1 Entropy-based uniform sampling on disconnected manifolds Input: Initialization scale σ 2 , number of samples N , number of iterations T , neighborhood parameter k, temperature τ , rejuvenation kernel K, rejuvenation step number M , constraint functions h, g. 1: // Initialization 2: Sample N particles xi0 ∼ N (0, σ 2 ) 3: Project the particles onto Σ (e.g., with Gauss-Newton method) 4: Apply a manifold-constrained rejuvenation kernel K 5: // Iterative entropy maximization 6: for t = 1, . . . , T do 7: Compute weights wti from the entropy-based rule in (3) 8: Resample particles according to the normalized weights wti 9: Apply a manifold-constrained rejuvenation kernel K for M steps 10: end for 11: return final particle set X (T ) Thus, particles in low-density regions have larger k-NN radii and receive larger resampling weights, encouraging redistribution toward underrepresented connected components. The temperature τ controls the aggressiveness of the resampling step. When τ is small, the update is conservative and moves mass more slowly. Larger values of τ accelerate mass redistribution but can amplify errors in the k-NN density estimate when the number of particles or mixing steps is insufficient. MASEM, described in pseudocode in Algorithm 1, works by first initializing feasible particles and then iterating resampling followed by manifold sampling rejuvenation. We initialize the particles by random i.i.d. sampling followed by projection onto the feasible set using Gauss-Newton steps on the squared slack 1 Slack(x) = ∥g(x)+ ∥2 + ∥h(x)∥2 (4) 2 to ensure approximate constraint feasibility. After initialization, we apply a feasibility-preserving kernel K for rejuvenation. This ensures that the particles are approximately uniformly distributed within components. Then, at each iteration, particles are resampled using the entropy-based resampling weights (3). After resampling, we apply K again to obtain mixed samples within each component. As local kernels are only approximately feasible in practice, we use a slack-penalized variant of the resampling weights to avoid over-replicating particles with large constraint violations in our experiments. These practical modifications are described in Appendix D.5.

Theoretical Analysis

Our theoretical guarantees address the global mass-redistribution problem, which cannot be solved by existing local samplers. We analyze Algorithm 1 in the mean-field regime, where the N → ∞ limit reduces the particle system to a deterministic evolution on the simplex of component weights. This parallels the idealized-limit viewpoint taken in related work on particle-based constrained sampling [1, 9]. Our main result (Theorem 1) shows that a single resampling step contracts the KL divergence to the uniform distribution on Σ contracts at rate (1 − τ /p), leading to exponential convergence. All proofs for this section can be found in Appendix A. 5.1

Preliminaries

We need two prerequisites for our main result: Each component must be initially populated, and the sampler has to sample accurately within components. We begin by showing that the initialization step of Algorithm 1 populates all components for N large enough. Lemma 2 (Component coverage at initialization). Suppose h, g have Lipschitz-continuous derivatives and satisfy LICQ, and that the feasible set Σ is bounded and decomposes into a finite number of connected components. We draw N points xi ∼ N (0, σ 2 ) for arbitrary σ > 0 and apply any descent method on the squared slack (4) satisfying sufficient decrease. Then P(every component is hit) → 1 as N → ∞. 4

Proof sketch. Take x∗ ∈ Σc arbitrary. The local error bound d(x, Σ) ≤ κ Slack(x)1/2 holds for some κ > 0 due to LICQ [28]. Since Σ is compact and has finitely many components, the inter-component distance δ := mini̸=j d(Σi , Σj ) is strictly positive [29, Ch. XI 4.4]. Taking L := {Slack(x) ≤ (δ/(3κ))2 } ensures that all x ∈ L ∩ Uδ/3 (Σc ) converge to Σc under any descent method satisfying sufficient decrease [30]. Therefore, the basin of attraction has positive Lebesgue measure, and standard coverage arguments yield the claim. To consider component-level convergence properties, we assume that the local kernel perfectly mixes chains within each component. This assumption decouples the analysis of the resampling in Algorithm 1 from any specific local sampler. Assumption 2 (Mixing Rejuvenation). The rejuvenation kernel K is uc -invariant on each Σc (i.e., uc K = uc , where uc = σc−1 1Σc ), and mixes sufficiently, such that after application the particles within Σc are approximately i.i.d. draws from uc . This assumption ensures that chains stay within their component during mixing. In practice, if components are close together, the chain might be able to jump between components. For the sake of our analysis, we treat these clusters of components as one “meta-component”, since, given enough time, the chains will mix uniformly within this cluster. This assumption allows us to consider component-level weights instead of individual particles. Let αc := Nc /N denote the fraction of particles in Σc . Then the empirical density after rejuvenation takes the mixture form pα (x) = αc(x) , where c(x) denotes the component containing x. Let α = (α1 , . . . , αC ) be the vector of component weights that induces the distribution pα . The sampling problem thus reduces to steering α ∈ ∆C−1 towards the maximum entropy distribution α∗ := (σ1 /S, . . . , σC /S). 5.2

Convergence Guarantees

We first derive the deterministic map that governs the evolution of α under resampling, then show that its iterates contract in KL divergence. Proposition 1 (Mean-field resampling map). Under Assumption 2, the k-NN radius of a particle xi ∈ Σc concentrates at εi,k ≍ (kσc /(N Vp αc ))1/p in the limit N → ∞ with kN → ∞, kN /N → 0 [27, Ch. 2]. Consequently, the resampling weights wi ∝ ετi,k induce the map α1−β (α∗ )β , (Φα)c := PC c 1−βc (αj∗ )β j=1 αj

β := τ /p.

(5)

The map (5) is a geometric mean of the current distribution α and the target α∗ . Its unique fixed point on ∆C−1 is α∗ . Thus, iterating Φ contracts α towards α∗ , with the parameter β = τ /p controlling the aggressiveness of rebalancing. We show that iterating Φ drives α to α∗ at a geometric rate. Theorem 1. Consider twice differentiable constraints h, g fulfilling LICQ, which define a bounded constrained set Σ that decomposes into a finite number of connected components. Suppose all components have positive measure with respect to the induced Hausdorff measure on Σ. Under Assumption 2 for any τ ∈ (0, p) and any schedule kN → ∞, kN /N → 0, the mean-field iterates αt+1 = Φ(αt ) of Algorithm 1 satisfy the following bound Σ DKL (pαt | pα∗ ) ≤ C0

τ 1− p

2t (6)

with 1 C0 = 4

(0)

αc αc max log ∗ − min log ∗ c c αc αc

!2 (7)

(0)

for the initial component weights αc pαt → uΣ in KL as t → ∞.

after the projection step of the algorithm. In particular,

α Proof sketch. We define the relative estimation error yt,c := αt,c ∗ as well as and zc := log y0,c and c t (1−β)t P (1−β) show using induction that yt,c = y0,c / j αj∗ y0,j . Writing at := (1 − β)t , this admits the exponential closed form X α∗ eat zc αt,c := c , Z(a) := αj∗ eazj = Eα∗ [eaz ], (8) Z(at ) j

from which we can see αt is an exponential tilt of α∗ , with inverse temperature at , which shrinks geometrically to zero. Substituting (8) into the definition of KL gives the decomposition X Σ DKL (pαt | pα∗ ) = at αt,c log y0,c − log Z(at ), c

so it suffices to upper-bound z̄t and lower-bound log Z(at ). Both bounds follow from properties ∗ of the cumulant generating function (CGF) Λ(a) := log Z(a) of z under αP . For the lower bound, Jensen’s inequality applied to the convex function eaz gives log Z(at ) ≥ at c α∗ zc . For the upper bound, we use basic properties of the CGF and Popoviciu’s variance inequality [31] to obtain a bound (0) (0) of Λ′′ (a) = Vνa [z] ≤ R2 /4, where R := maxc log(αc /αc∗ ) − minc log(αc /αc∗ ). Integrating then 2 Σ gives z̄t ≤ z̄α∗ +at R /4. As the linear terms cancel, we end with DKL (pαt | pα∗ ) ≤ (R2 /4)(1−β)2t , which is the bound in (6). Using this result we derive an asymptotic bound on the worst case number of iterations. Corollary 1. Assume the conditions of Theorem 1 hold. In the worst case, the number of iterations Σ t to reach a KL-divergence DKL (pαt | pα∗ ) ≤ ε is log(1/ε) + log log N t∈O . τ Proof sketch. We show C0 ∈ O (log N )2 . Bounding (6) from above by ε lets us derive log(ε/C0 ) log(1/ε) + log C0 log(1/ε) + log log N = t≈ ∈O . −2 log(1 − τp ) τ 2 log 1 − τp Component loss We note that Theorem 1 is only stated in the mean field. It might happen that no particles of a particular component are sampled, resulting in that component vanishing. The probability of this extinction is Pc = (1 − (Φα)c )N ≤ exp −N (Φα)c . Using Proposition 1 we get (Φα)c ∝ αc1−β (αc∗ )β . By Lemma 2, each component contains at least one sample, bounding αc = Nc /N ≥ N −1 . Since αc∗ is constant, we can derive an asymptotic upper bound on the extinction-probability as Pc ∈ O exp(−N αc1−β ) = O exp(−N β ) . We investigated this effect and found that, in practice, just four chains per component is enough (see Appendix C.1).

Experimental Analysis

To demonstrate how MASEM improves existing manifold samplers, we consider Non-Linear Hit&-Run (NHR) [24] and OLLA [9] with- and without the proposed resampling logic. In addition, we compare against Sequentially Constrained MC (SCMC) [17] which anneals soft constraints. To understand how MASEM compares against explicit component discovery and importance sampling, we implement Cluster-NHR, which clusters samples and computes resampling weights based on the estimated cluster volumes. We report the squared Sinkhorn distance and the averaged maximum slack violation. Since SCMC cannot guarantee exact constraint feasibility, we perform an additional projection step before reporting the metrics. Implementation details are listed in Appendix D and additional figures can be found in Appendix B. 6

Table 1: Final W22 (Sinkhorn) distance performance across five seeds. Entries are mean ± 95% confidence interval. Bold indicates better than counterpart for t-test with p < 0.01 and HolmBonferroni correction. MASEM-Sampler outperform their counterpart by an order of magnitude on all non-connected problems.

Problem Connect. Disks (3d) Disconn. Disks (3d) Seven Lobes (2d) Sine (2d) Swiss Roll (2d) 6.1

NHR

MASEM-NHR

OLLA

MASEM-OLLA

.06 ± .01 1.98 ± .12 .37 ± .02 100.21 ±4.09 6.08 ± .75

.00 ± .00 .01 ± .01 .05 ± .02 .13 ± .06 .08 ± .02

.04 ± .00 1.95 ± .11 .89 ± .08 128.38 ±3.79 18.96 ±1.15

.00 ± .00 .02 ± .00 .05 ± .02 .13 ± .07 .29 ± .09

Synthetic Benchmarks

We first evaluate sampling on lowdimensional manifolds: (1)-(2) two disks that are embedded on a sphere in R3 and we vary whether the disks are connected or not, (3) the seven lobes density from Jeon et al. [9] with constant f , (4) a sine equality constraint with reducing amplitude cut into disconnected components inequality constraint, and (5) a product of an Archimedean spiral manifold and randomly sampled circles combined with a nonlinear inequality we label swiss roll. The benchmarks are visualized in Figure 6, and we report the specific constraints and ground truth sampling method in Appendix D.2. For each problem, we run 2 000 independent chains in parallel for 5 000 steps and collect the last sample per chain.

Table 2: Final W22 (Sinkhorn) distance and slack across five seeds. Entries are mean ± 95% confidence interval. Bold indicates best for t-test with p < 0.01 and HolmBonferroni correction. W22

SCMC

Cluster-NHR

MASEM-NHR

.12 ± .01 .59 ± .23 .12 ± .06 126.79 ±10.85 1.10 ± .42

.05 ± .03 .19 ± .14 .14 ± .04 .91 ±1.04 .42 ± .26

.00 ± .00 .01 ± .01 .05 ± .02 .13 ± .06 .08 ± .02

Slack ×10−3

SCMC

Cluster-NHR

MASEM-NHR

C. Disks D. Disks Seven Lobes Sine Swiss Roll

20.46 ± 1.12 29.97 ±13.61 3.00 ± 1.89 9.38 ± 1.48 2.77 ± .71

.02 ± .03 .02 ± .03 .00 ± .00 .00 ± .00 .00 ± .00

.01 ± .01 .01 ± .01 .00 ± .00 .01 ± .01 .00 ± .00

C. Disks D. Disks Seven Lobes Sine Swiss Roll

Local Kernels Fail Under Disconnected Components. In Figure 1 we identify a common failure mode of standard manifold samplers: In the presence of disconnected components, these local methods fail to correctly allocate mass; both methods over sample the lower disk on the sphere and under sample the left-most arc of the sine. While the connected disks problem permits the chains to mix well in the single feasible level set, separating the two spheres on the manifold leads to an increase in sampling error, as shown in Table 1.

Temperature τ

Sampling Accuracy & Constraint Violation of MASEM. As shown in Table 1, MASEM-NHR and MASEM-OLLA outper22 ↓ form their counterparts across all problems with a disconnected 0.1 feasible set by an order of magnitude in the Sinkhorn distance 2 W2 . Further, MASEM minimizes the constraints, matching (or 0.5 even outperforming) their local counterparts (Figure 5). We 1.0 also see in Table 2 that MASEM yields superior performance compared to the global baselines. In particular, we observe that 1.5 while SCMC effectively yields constraint satisfying samples 2.0 on the seven lobes problem, it fails to do so in particular on 10 50 100 150 200 problems with highly non-linear constraints such as the swiss Rejuvenation steps M roll problem. Cluster-NHR, in contrast, achieves low constraint violations but higher sampling error than MASEM on all prob- Figure 2: Influence of τ and M hyperparameters for MASEM-NHR lems. on the 7 lobes problem. We mean Effect of Hyperparameters τ and M . We further examine W 2 distance across 5 seeds with 2 the influence of the weight scaling hyperparameter τ and the 95% CI. number of rejuvenation steps M on the seven lobes problem. As shown in Figure 2, our method is robust to variations of both 0.07 ±0.03

0.10 ±0.04

0.12 ±0.03

0.16 ±0.03

0.20 ±0.05

0.04 ±0.01

0.05 ±0.01

0.06 ±0.01

0.06 ±0.02

0.08 ±0.02

0.06 ±0.02

0.04 ±0.01

0.05 ±0.02

0.03 ±0.01

0.04 ±0.02

0.76 ±0.62

0.03 ±0.01

0.04 ±0.03

0.07 ±0.03

0.06 ±0.02

3.57 ±3.16

3.84 ±2.57

0.94 ±0.71

0.11 ±0.03

0.08 ±0.03

1.0

0.8

0.6

0.4

0.2

0.0

NHR

OLLA

22 ↓ 10

SCMC

Cluster-NHR

MASEM-NHR (Ours)

22 ↓

Wallclock (s) ↓ 102

MASEM-OLLA (Ours)

Wallclock (s) ↓

102

103 102

100

101

32 64 128 d (dimension)

256

101

10−1

32 64 128 d (dimension)

256

(a) Increasing ambient and fixed manifold dimension with increasing m = d − 3.

32 64 128 d (dimension)

256

32 64 128 d (dimension)

256

(b) Increasing ambient and manifold dimension with fixed m = 5.

Figure 3: Sampling performance and wallclock time as the dimension d increases (with l = |C| = 5). MASEM achieves lowest W22 (Sinkhorn) distances. hyperparameters. However, excessively large τ values paired with short mixing times lead to high Sinkhorn distances to ground truth samples. Additionally, τ = 0.1 might lead to component loss, as demonstrated in Appendix C.1. 6.2

Scaling Under High-Dimensionality and Large Number of Constraints

We assess the robustness and scalability of MASEM-NHR and MASEM-OLLA using a synthetic stress-test problem that enables explicit control of ambient dimension d as well as the number of equality and inequality constraints and number of disconnected components (m, l, |C|). In our case, we guarantee different components by sampling a set of |C| disjoint spheres, which we embed into a subspace of Rd by random linear projections defined by the m equality constraints. Scaling Under Increasing Manifold Dimension. With fixed numbers of constraints m = l = 5 we study the effect of an increasing manifold dimension as d increases. Figure 3b shows that MASEM consistently outperforms NHR and OLLA as well as Cluster-NHR by an order of magnitude. In addition, we note that the penalty-based approach of SCMC deteriorates strongly as d increases. Further results in Figure 7 show a similar trend for the KL divergence while the maximum constraint violation stays below 0.1 for MASEM. Scaling Under Number of Constraints. Figure 3a illustrates the performances as the ambient dimension d increases while the manifold dimension is fixed to 2. We see that the gap between MASEM and its counterparts remains constant across dimensions while the W22 distance only slightly increases with d. The wallclock times are competitive compared to the baselines that evaluate the constraint Jacobian (OLLA, NHR, Cluster-NHR). This shows that our method adds only little overhead in practice. Additional plots for the KL divergence and constraint violation in Figure 7 further support these results. 6.3

Robotics Applications

Motion Planning. We consider trajectory sampling, where we sample constraint satisfying trajectories through an obstacle course. We sample the trajectories of a 2d pointmass across two different obstacle courses, one that is a regular 4 × 4 grid and another one with 20 randomly placed obstacles (Figure 4). Each trajectory is parameterized by a spline defined by 3 waypoints, resulting in a 9 dimensional problem. Random obstacle motion planning has 934 constraints in total, while grid based planning has 774 inequality constraints. Grasp Sampling. Secondly, we evaluate MASEM for grasping. We optimize grasps on a capsule geometry (a common shape for collision checking) with three fingers. A band around the capsule shall not be touched, mirroring grasping in the real world, where certain objects cannot be grasped at arbitrary points.2 Each finger applies a force, such that in the end the fingers counteract gravity pulling on the object. This results in a 18-dimensional problem, with 9 equality and 48 inequality constraints. Results. Due to the lack of ground truth samples, we report the feasible entropy instead of the Sinkhorn distance. We compute it by multiplying the entropy of the feasible samples and the fraction 2 For example, a robot grasping a plate should not put its finger into the food.

Table 3: Performance on robotics problems across seeds. Entries are mean ± 95% confidence interval. Column-wise bold highlights indicate better than counterpart for t-test with p < 0.01 and Holm-Bonferroni correction. Method NHR Cluster-NHR MASEM-NHR OLLA MASEM-OLLA SCMC

Planning Random Obst.

Planning Grid Obst.

Grasping

Feas. Ent. ↑

Slack ↓

Feas. Ent. ↑

Slack ↓

Feas. Ent. ↑

Slack ↓

6.85 ± .37 6.49 ±1.43 7.93 ± .22 .86 ± .15 6.11 ± .47 1.43 ± .69

.00 ± .00 .00 ± .00 .00 ± .00 .07 ± .01 .01 ± .00 .05 ± .04

6.50 ± .14 4.41 ±1.27 7.49 ± .16 1.47 ± .24 5.35 ± .38 2.26 ± .26

.00 ± .00 .00 ± .00 .00 ± .00 .03 ± .01 .01 ± .00 .02 ± .01

3.61 ± .23 3.60 ±1.05 4.79 ± .16 1.61 ± .09 4.09 ± .18 .22 ± .04

.00 ± .00 .00 ± .00 .00 ± .00 .37 ± .01 .01 ± .00 .08 ± .00

of feasible samples (made precise in Appendix D.1). Table 3 shows an increase of sample entropy when applying MASEM, regardless of sampler. This is illustrated in Figure 4 for random obstacles. We observe that MASEM-NHR samples trajectories through the narrow gaps at the center of the course, which are missed by standard NHR. For OLLA the performance gap is even larger, as the default sampler only generates a locally distributed sample sets but fails to sample the outer regions of the domain. Moreover, OLLA generates many paths crossing over obstacles. While MASEM-NHR outperforms MASEM-OLLA on all problems, OLLA benefits more from the resampling applied by MASEM. This suggests two things: first, NHR-based samplers are more suited for robotics than OLLA-based samplers (at least within this evaluation set), and secondly, MASEM greatly improves previously unsuited samplers and makes them competitive in this domain. We believe these are valuable for future research on sampling in robotics. Further, our approach outperforms alternative global sample allocation strategies: SCMC performs poorly on all tasks; MASEM-NHR outperforms its clustering counterpart on all tasks as well, with a large gap on the regular grid design. The performance of MASEM on the motion planning problems demonstrates its ability to scale to problems with a high number of constraints. NHR

OLLA

SCMC

MASEM-NHR

MASEM-OLLA

Figure 4: Plots of 50 samples on the random motion planning problem. MASEM-based approaches yield the highest sample entropy across all methods.

Conclusion & Future Work

We presented Manifold Sampling via Entropy Maximization (MASEM), a resampling-based approach for uniform sampling on disconnected manifolds implicitly defined by constraints. We proved that MASEM minimizes the KL divergence to the uniform target distribution exponentially in the number of resampling steps and provided worst-case bounds in the mean-field. Building on this, we analyzed two instantiations of MASEM on several synthetic and real world problems. While the method introduces new parameters to tune and a small runtime overhead (which we discuss together with theoretical limitations in Appendix E), it results in a significantly lower Sinkhorn distance to the ground truth distribution compared to the baselines. Moreover, it is designed to work with any constrained MCMC-Sampler. Future research could investigate which properties make a sampler work well with MASEM. In this work we focus on entropy maximization, which leads to uniform sampling across the manifold. Future work could explore generalizations to non-uniform distributions. We see potential in applying MASEM as a data generation method in robotics and other application domains. An intriguing avenue of research is investigating practical design choices and comparing them to current heuristic approaches in practical applications.

Acknowledgments and Disclosure of Funding The authors sincerely thank Paula Cordero Encinar for discussions during early stages of the project, and Eckart Cobo-Briesewitz for his feedback on the manuscript. This research was funded by the Amazon Fulfillment Technologies and Robotics team. This work has been supported by the German Federal Ministry of Research, Technology and Space (BMFTR) under the Robotics Institute Germany (RIG).

References [1] Ruqi Zhang, Qiang Liu, and Xin Tong. Sampling in Constrained Domains with OrthogonalSpace Variational Gradient Descent. Advances in Neural Information Processing Systems, 35: 37108–37120, 2022. [2] D. C. Rapaport. The Art of Molecular Dynamics Simulation. Cambridge University Press, Cambridge, 2nd edition, 2004. [3] Zachary Kingston, Mark Moll, and Lydia E. Kavraki. Sampling-Based Methods for Motion Planning with Constraints. Annual Review of Control, Robotics, and Autonomous Systems, 1(1): 159–185, 2018. [4] Thomas Power, Rana Soltani-Zarrin, Soshi Iba, and Dmitry Berenson. Sampling Constrained Trajectories Using Composable Diffusion Models. In IROS 2023 Workshop on Differentiable Probabilistic Robotics: Emerging Perspectives on Robot Learning, 2023. [5] Carlos Florensa, David Held, Markus Wulfmeier, Michael Zhang, and Pieter Abbeel. Reverse Curriculum Generation for Reinforcement Learning. In Conference on Robot Learning, pages 482–495. PMLR, 2017. [6] Marc Toussaint, Cornelius V Braun, Eckart Cobo-Briesewitz, Sayantan Auddy, Armand Jordana, and Justin Carpentier. Constrained Sampling to Guide Universal Manipulation RL. arXiv preprint arXiv:2602.08557, 2026. [7] Eckart Cobo-Briesewitz, Tilman Burghoff, Denis Shcherba, Armand Jordana, and Marc Toussaint. Stability-Guided Exploration for Diverse Motion Generation, 2026. [8] Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I Jordan. An Introduction to MCMC for Machine Learning. Machine learning, 50(1):5–43, 2003. [9] Kijung Jeon, Michael Muehlebach, and Molei Tao. Fast Non-Log-Concave Sampling under Nonconvex Equality and Inequality Constraints with Landing. In The Thirty-Ninth Annual Conference on Neural Information Processing Systems, 2026. [10] Emilio Zappa, Miranda Holmes-Cerfon, and Jonathan Goodman. Monte Carlo on Manifolds: Sampling Densities and Integrating Functions. Communications on Pure and Applied Mathematics, 71(12):2609–2647, 2018. [11] Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz. Hybrid Monte Carlo methods for sampling probability measures on submanifolds. Numerische Mathematik, 143(2):379–421, 2019. [12] David E Kaufman and Robert L Smith. Direction Choice for Accelerated Convergence in Hit-and-Run Sampling. Operations Research, 46(1):84–95, 1998. [13] An T Le, Georgia Chalvatzaki, Armin Biess, and Jan R Peters. Accelerating Motion Planning via Optimal Transport. Advances in Neural Information Processing Systems, 36:78453–78482, 2023. [14] Zechu Li, Rickmer Krohn, Tao Chen, Anurag Ajay, Pulkit Agrawal, and Georgia Chalvatzaki. Learning Multimodal Behaviors from Scratch with Diffusion Policy Gradient. Advances in Neural Information Processing Systems, 37:38456–38479, 2024. [15] Natalie Maus, Kaiwen Wu, David Eriksson, and Jacob Gardner. Discovering Many Diverse Solutions with Bayesian Optimization. In International Conference on Artificial Intelligence and Statistics, pages 1779–1798. PMLR, 2023. 10

[16] Don O Loftsgaarden and Charles P Quesenberry. A Nonparametric Estimate of a Multivariate Density Function. The Annals of Mathematical Statistics, 36(3):1049–1051, 1965. [17] Shirin Golchi and David A Campbell. Sequentially Constrained Monte Carlo. Computational Statistics & Data Analysis, 97:98–113, 2016. [18] Tony Lelievre, Mathias Rousset, and Gabriel Stoltz. Langevin dynamics with constraints and computation of free energy differences. Mathematics of computation, 81(280):2071–2125, 2012. [19] Marcus Brubaker, Mathieu Salzmann, and Raquel Urtasun. A Family of MCMC Methods on Implicitly Defined Manifolds. In Artificial Intelligence and Statistics, pages 161–172. PMLR, 2012. [20] Simon Byrne and Mark Girolami. Geodesic Monte Carlo on Embedded Manifolds. Scandinavian Journal of Statistics, 40(4):825–845, 2013. [21] Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from a Log-Concave Distribution with Projected Langevin Monte Carlo. Discrete & Computational Geometry, 59(4):757–783, 2018. [22] Yunbum Kook, Yin-Tat Lee, Ruoqi Shen, and Santosh Vempala. Sampling with Riemannian Hamiltonian Monte Carlo in a Constrained Space. Advances in neural information processing systems, 35:31684–31696, 2022. [23] Thomas Power and Dmitry Berenson. Constrained Stein Variational Trajectory Optimization. IEEE Transactions on Robotics, 40:3602–3619, 2024. [24] Marc Toussaint, Cornelius V Braun, and Joaquim Ortiz-Haro. NLP Sampling: Combining MCMC and NLP Methods for Diverse Constrained Sampling. arXiv preprint arXiv:2407.03035, 2024. [25] Francisco J Aragón, Miguel A Goberna, Marco A López, and Margarita ML Rodríguez. Nonlinear Optimization. Springer, 2019. [26] Hubertus Th. Jongen, P. Jonker, and F. Twilt. Nonlinear Optimization in Finite Dimensions: Morse Theory, Chebyshev Approximation, Transversality, Flows, Parametric Aspects. Nonconvex Optimization and Its Applications. Springer New York, 2000. [27] Gérard Biau and Luc Devroye. Lectures on the Nearest Neighbor Method, volume 246. Springer, 2015. [28] R. Andreani, G. Haeser, R.W. Prado, and L.D. Secchin. Primal-dual global convergence of an augmented Lagrangian method under the error bound condition. Technical report, University of São Paulo, 2025. [29] James Dugundji. Topology. Allyn and Beacon, Boston, 1966. [30] Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear Convergence of Gradient and ProximalGradient Methods Under the Polyak-Łojasiewicz Condition. In Paolo Frasconi, Niels Landwehr, Giuseppe Manco, and Jilles Vreeken, editors, Machine Learning and Knowledge Discovery in Databases, pages 795–811. Springer International Publishing, 2016. [31] Tiberiu Popoviciu. Sur les équations algébriques ayant toutes leurs racines réelles. Mathematica, 9(129-145):20, 1935. [32] Martin J Wainwright and Michael I Jordan. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends® in Machine Learning, 1(1-2):1–305, 2008. [33] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Yash Katariya, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye WandermanMilne, and Qiao Zhang. JAX: Composable transformations of Python+NumPy programs, 2018. 11

[34] L. F. Kozachenko and N. N. Leonenko. Sample Estimate of the Entropy of a Random Vector. Problems of Information Transmission, 23(1-2):95–101, 1987. [35] Richard M. Murray, Zexiang Li, and S. Shankar Sastry. A Mathematical Introduction to Robotic Manipulation. CRC Press, Boca Raton, 1994. ISBN 978-1-315-13637-0. [36] Michael Laskin, Hao Liu, Xue Bin Peng, Denis Yarats, Aravind Rajeswaran, and Pieter Abbeel. CIC: Contrastive Intrinsic Control for Unsupervised Skill Discovery. arXiv preprint arXiv:2202.00161, 2022. [37] Thomas Beret, Richard Samworth, and Ming Yuan. Efficient multivariate entropy estimation via k-nearest neighbour distances. The Annals of Statistics, 47(1):288–318, 2019. [38] Jasper Snoek, Hugo Larochelle, and Ryan Adams. Practical Bayesian Optimization of Machine Learning Algorithms. In F. Pereira, C.J. Burges, L. Bottou, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012. [39] Ben Harwood, Amir Dezfouli, Iadine Chades, and Conrad Sanderson. Approximate Nearest Neighbour Search on Dynamic Datasets: An Investigation. In Mingming Gong, Yiliao Song, Yun Sing Koh, Wei Xiang, and Derui Wang, editors, AI 2024: Advances in Artificial Intelligence, volume 15443, pages 95–106. Springer Nature Singapore, 2025. [40] Qiang Liu and Dilin Wang. Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.

Formal Lemmas and Proofs

We begin by formalizing Lemma 1 introduced in Section 3. Lemma 1. Let Σ be measurable with 0 < σΣ (Σ) < ∞, and let uΣ (x) :=

1 S

for S = σΣ (Σ)

denote the uniform density on Σ with respect to σΣ . Then it holds that log S = H(uΣ ) ≥ H(ρΣ ) for all densities ρΣ supported on Σ with equality iff ρ = uΣ almost everywhere. Proof. We compute the entropy of the uniform density as Z 1 1 H(uΣ ) = − log dσΣ = log S S S Σ Next, consider the relative entropy of ρ w.r.t. uΣ : Z Z Z ρ dσΣ = DKL (ρ | uΣ ) = ρ log S dσΣ + ρ log ρ log ρ dσΣ uΣ Σ Σ Σ = log S − H(ρ) = H(uΣ ) − H(ρ). Now, by Gibbs’ inequality 0 ≤ DKL (ρ | uΣ ) with equality iff ρ = uΣ almost everywhere, proving the statement. Next, we restate the results of Section 5 and provide the complete proofs. Lemma 2 (Component coverage at initialization). Suppose h, g have Lipschitz-continuous derivatives and satisfy LICQ, and that the feasible set Σ is bounded and decomposes into a finite number of connected components. We draw N points xi ∼ N (0, σ 2 ) for arbitrary σ > 0 and apply any descent method on the squared slack (4) satisfying sufficient decrease. Then P(every component is hit) → 1 as N → ∞. Proof. Consider gradient descent on the squared constraint violations 1 ∥g(x)+ ∥2 + ∥h∥2 . Slack(x) = 2 We will show that each component Σc admits a basin of attraction Ac under gradient descent on Slack(x) with µ(Ac ) > 0. Then, drawing i.i.d. from an uncorrelated Gaussian results in each set of positive measure having a positive probability of being hit. Let the probability of hitting Ac be pc > 0. The probability that no particle among N i.i.d. draws lands in Ac is (1 − pc )N → 0, where M0 is the volume of the sampling region. A union bound over the finitely many components gives the claim. It is left to show that each component admits a basin of attraction of positive measure. To that end, we show that for each component there exists a nonempty open set Uc , such that gradient descent on points in this set ends up in Σc . Then Uc ⊆ Ac and therefore µ(Ac ) > 0. Fix Σc and x∗ ∈ Σc . By LICQ, the squared slack Slack(x) fulfills the error bound [28]. That means there exists ε, κ > 0 such that d(x, Σ) ≤ κ Slack(x)1/2 for all x ∈ Uε (x∗ ) . (9) Since the components are compact, there exists δ := min d(Σi , Σj ) > 0 i̸=j √ by standard topological arguments [29, Ch. XI 4.4]. We choose α > 0 such that α ≤ δ/(3κ), and ∗ define L := {x ∈ Uε (x ) | Slack(x) < α}. For any x ∈ L, √ d(x, Σ) ≤ κ Slack(x)1/2 < κ α ≤ δ/3. By Karimi et al. [30], gradient descent from any x0 ∈ L converges to Σ, and, since Slack(xk+1 ) ≤ Slack(xk ), iterates remain in L. Therefore any x0 ∈ Uc := L ∩ Uδ/3 (Σc ) converges into Σc . The set Uc is open and nonempty, hence µ(Ac ) ≥ µ(Uc ) > 0, which completes the proof. 13

Proposition 1 (Mean-field resampling map). Under Assumption 2, the k-NN radius of a particle xi ∈ Σc concentrates at εi,k ≍ (kσc /(N Vp αc ))1/p in the limit N → ∞ with kN → ∞, kN /N → 0 [27, Ch. 2]. Consequently, the resampling weights wi ∝ ετi,k induce the map α1−β (α∗ )β , (Φα)c := PC c 1−βc (αj∗ )β j=1 αj

β := τ /p.

(5)

Proof. The mass assigned to component c after normalization is P α1−β σ β Nc w c αc (σc /αc )β xi ∈Σc wi = P c 1−βc β , =P (Φα)c = P P =P β σj j xk ∈Σj wk j Nj wj j αj (σj /αj ) j αj where we used Nc = αc N and the concentration wc ∝ (σc /αc )β . Substituting σc = Sαc∗ and cancelling S β gives (5). Theorem 1. Consider twice differentiable constraints h, g fulfilling LICQ, which define a bounded constrained set Σ that decomposes into a finite number of connected components. Suppose all components have positive measure with respect to the induced Hausdorff measure on Σ. Under Assumption 2 for any τ ∈ (0, p) and any schedule kN → ∞, kN /N → 0, the mean-field iterates αt+1 = Φ(αt ) of Algorithm 1 satisfy the following bound 2t τ Σ DKL (pαt | pα∗ ) ≤ C0 1 − p with 1 C0 = 4

(0)

αc αc max log ∗ − min log ∗ c c αc αc (0)

for the initial component weights αc pαt → uΣ in KL as t → ∞.

(6) !2 (7)

after the projection step of the algorithm. In particular,

Proof. We prove the claim in three steps: (i) derive a closed form for the iterates, (ii) decompose the KL, and (iii) bound each term. Step 1: Closed form of the iterates. We recall that under Assumption 2, the k-NN radius of a particle xi ∈ Σc concentrates at εi,k ≍ (kσc /(N Vp αc ))1/p in the limit N → ∞ with kN → ∞, kN /N → 0 [27, Ch. 2]. Hence, our weights concentrate as τ /p τ /p k σc σc wit = εk (x)τ /p ≍ ∝ . N Vp αc αc Define yt,c := αt,c /αc∗ . Substituting αt,c = yt,c αc∗ into (5) then yields: 1−β αc∗ yt,c (yt,c αc∗ )1−β (αc∗ )β . (Φαt )c = P = P ∗ 1−β (α∗ )β ∗ 1−β j j (yt,j αj ) j αj yt,j

Using this expression, we can characterize how the gap yt,c changes over time: 1−β

yt+1,c =

αc∗ yt,c αt+1,c 1−β = yt,c /Zt = P 1−β αc∗ αc∗ j αj∗ yt,j

(10)

P 1−β with Zt := j αj∗ yt,j independent of component c. This implies the following recurrence which we prove by induction: yt,c =

at y0,c

Z̃t

at := (1 − β)t ,

Z̃t :=

X j

at αj∗ y0,j .

(11)

Base case t = 0. We have a0 = 1 so Z̃0 = which is trivially true.

(0) j αj = 1. This reduces the statement to y0,c = y0,c

Inductive step. We have 1−β yt,c ∗ 1−β j αj yt,j yat 1−β

yt+1,c = P

By (10)

0,c

∗ j αj

By (11)

yat 1−β 0,j

Zt at (1−β) y0,c Zt1−β

P 1 ∗ at (1−β) 1−β j αj y0,j Z t at+1 y0,c = P ∗ at+1 , j αj y0,j

By at (1 − β) = (1 − β)t+1 = at+1

which completes the induction. P Now we can express this in tilted exponential form by writing zc := log y0,c and Z(a) := j αj∗ eazj . Thus (11) becomes α∗ eat zc . (12) αt,c = c Z(at ) This expression quantifies how much component c was initially over- or under-represented relative to the target and still is at time t. Step 2: KL decomposition. Using (12), we decompose the KL divergence as follows: X X α∗ eat zc Σ DKL (pαt | pα∗ ) = αt,c log c ∗ = αt,c (at zc − log Z(at )) = at z̄t − log Z(at ), (13) αc c c P P where z̄t := c αt,c zc is the first moment of z under αt and Z(a) := j αj∗ eazj is the normalizing constant. Step 3: Bounds on each term. Note that αt is itself a tilted distribution from the exponential form αct = αc∗ eat zc /Z(at ) and that Z(a) is its moment generating P ∗ function. We can therefore study how its mean z̄t deviates from the untilted mean a P z̄α∗ = c αc zc by treating the tilt parameter as continuous and tracking the mean µ(a) := c νa,c zc as a varies from 0 (where ν0 = α∗ ) to at (where νat = αt ). Upper Bound via Popoviciu’s Inequality. Define the log-partition function Λ(a) := log Z(a). Differentiation of this cumulant generating function yields the first moment of the tilted distribution: P ∗ azc X α zc e ′ Λ (a) = c c = νc(a) zc = Eνa [z] = µ(a). (14) Z(a) c Differentiating once more yields the variance thereof [32]: !2 µ′ (a) = Λ′′ (a) =

νa,c zc2 −

νa,c zc

= Vνa [z].

(15)

Since z only takes values in [minc zc , maxc zc ], we can apply Popoviciu’s variance inequal(0) ity [31] yields Vνa (z) ≤ R2 /4 for all a. In the last step we define R := maxc log(αc /αc∗ ) − (0) minc log(αc /αc∗ ) for clarity. Ra Integrating this bound over time from 0 to at yields a bound on the gap µ(at )−µ(0) = 0 t µ′ (a) da ≤ 2 at R4 . This establishes a bound for the first term in (13) since at z̄t = at µ(at ) ≤ at µ(0) + a2t R2 /4. 15

(16)

Lower Bound via Jensen’s Inequality. We apply Jensen’s inequality directly to the convex exponential function in the log-partition function: ! X X ∗ azc ∗ Z(a) = αc e ≥ exp a αc zc = eaz̄α∗ . c

Hence, we get log Z(at ) ≥ at z̄α∗ .

(17)

The KL Bound. Substituting (16) and (17) into (13): 2 Σ DKL (pαt | pα∗ ) ≤ at z̄α∗ + at R2 /4 − at z̄α∗ = R4 (1 − β)2t . Since pαt and uΣ differ only in their component weights and Assumption 2 ensures within-component Σ Σ uniformity, DKL (pαt | uΣ ) = DKL (pαt | pα∗ ), yielding (6). Corollary 1. Assume the conditions of Theorem 1 hold. In the worst case, the number of iterations Σ t to reach a KL-divergence DKL (pαt | pα∗ ) ≤ ε is log(1/ε) + log log N t∈O . τ Proof. We start by providing an upper bound on the initialization constant C0 (Eq. (7)) by considering worst-case initialization. We assume using Lemma 2 that each component contains at least one (0) sample, therefore αc ∈ N1 , 1 . This lets us derive (0)

max log c

1 αc ≤ max log ∗ = − min log αc∗ c c αc∗ αc

and similarly (0)

min log c

αc 1/N ≥ min log ∗ = − log(N ) − max log αc∗ . ∗ c αc αc

Therefore (0)

Z = max log c

(0)

αc αc − min log ∗ ≤ log N + max log αc∗ − min log αc∗ ∗ c c c αc αc

and, since Z > 0, we know that C0 fulfills 2 1 1 log N + max log αc∗ − min log αc∗ C0 = Z 2 ≤ c 4 4 2 ∗ and therefore C0 ∈ O (log N ) since αc = σc is independent of N . By solving C0 (1 − τp )2t = ε for t we then get log(ε/C0 ) log(1/ε) + log C0 = t≈ ∈O −2 log(1 − τp ) 2 log 1 − τp

log(1/ε) + log log N τ

since β = τp ∈ (0, 1) and therefore β ≤ − log(1 − β).

Additional Figures

B.1

Synthetic Benchmarks Convergence and Qualitative Results

We list the full learning curves in Figure 5. These curves illustrate the evolution of Sinkhorn and KL divergence as well as constraint violation. More details on the metrics can be found in Appendix D.1. Overall we observe strong performance and fast convergence of both MASEM variants on all problems. While KL divergence is low for all methods on the connected disk problem, separating the disks as depicted in Figure 6 leads to a performance deterioration of all baselines as shown in Figure 5. Even for the connected disks, MASEM leads to faster convergence. 16

NHR

Disconnected Disks

Connected Disks

101

OLLA

SCMC

22 ↓

10−1

10−2 10

10−2

MASEM-OLLA (Ours)

Slack ↓ 10

100

10−3

MASEM-NHR (Ours)

KL ↓

100

−1

10−3

−3

10−5

10−4 0

101

2000

Iteration 22 ↓

4000

2000

Iteration

4000

10−1 10−2

2000

Iteration

4000

Slack ↓ 10−1

0.75

100

KL ↓

1.00

0.50

10−3

0.25

10−5

0.00 0

2000

Iteration 22 ↓

4000

2.0

100

Seven lobes

Cluster-NHR

×10−1

2000

Iteration

4000

2000

Iteration

4000

Slack ↓

KL ↓ 100

1.5 1.0

10−3

0.5

10−6

10−1

0.0 0

102

2000

Iteration 22 ↓

4000

2000

Iteration

4000

Iteration

4000

10−1

Sine

2000

Slack ↓

KL ↓

101

10−3

100

−2

10−5 10−7

10−1 0

Swiss roll

2000

Iteration 22 ↓

4000

2000

Iteration

4000

KL ↓

101

10−1

100

−2

2000

Iteration

4000

Slack ↓

100 10−2

10−4 10−6

10−1 10−3 10−2

2000

Iteration

4000

10−8 0

2000

Iteration

4000

2000

Iteration

4000

Figure 5: Convergence over 5000 iterations on synthetic benchmarks. From left to right: (1) squared Sinkhorn distance to ground truth samples, (2) KL divergence to pairwise distances of ground truth samples, and (3) mean maximum constraint violation Ex [max{|h(x)|, g + (x)}]. Solid lines and shaded bands show the mean and 95% CI over five independent runs. Both MASEM-NHR and MASEM-OLLA quickly decrease constraint violation and maintain it there, while achieving the lowest W22 and KL values.

Disconn. Disks

Seven Lobes

Sine

Swiss Roll

Ground Truth

MASEM-NHR

NHR

Conn. Disks

Figure 6: Qualitative results for NHR and MASEM-NHR on all 2d benchmark problems. Black lines show equality constraints, and green shaded areas mark the feasible region(s) by inequality constraints. NHR struggles with moving across components, as can be seen in the disconnected disk and the sine benchmark. Moreover, on manifolds with high curvature, NHR also struggles to distribute particles correctly, since it works by locally linearizing the equality constraints. This can be seen with the seven lobes and the swiss roll benchmark. In both cases, MASEM is able to improve the distribution of particles. B.2

Scaling Analysis

Scaling d (m = 5)

NHR

OLLA

22 ↓

Cluster-NHR

101

102

100

101

10−1

Slack ↓

103

10−1

32 64 128 d (dimension)

256

32 64 128 d (dimension)

256

KL ↓

10−3 16

32 64 128 d (dimension)

10−1

32 64 128 d (dimension)

256

Slack ↓

102

10−4

10−3 16

32 64 128 d (dimension)

10−1

−1

256

Comp. Volume TV ↓

100 10

MASEM-OLLA (Ours)

Comp. Volume TV ↓ 0

101

22 ↓

MASEM-NHR (Ours)

KL ↓

103

Scaling d and m=d−3

SCMC

10−7 16

32 64 128 d (dimension)

256

32 64 128 d (dimension)

256

32 64 128 d (dimension)

256

Figure 7: Stress test for scaling ambient and manifold dimension (first row) and scaling the ambient dimension while the manifold dimension stays fixed. (second row). Curves depict mean and 95% CI across 5 seeds. We see that MASEM outperforms alternative approaches by an order of magnitude. Performances are stable across ambient dimensions, indicating that manifold dimension is the main difficulty for samplers. Soft constraint approaches suffer in particular with high-dimensional manifolds. We report the full scaling results in Figure 7. The plots report Sinkhorn distance and pairwise distance KL divergence w.r.t. to ground truth samples. In addition, we estimate the total variation distance (TV) to the ground truth component volume estimates of the resulting empirical distributions as well as expected maximum constraint violations. We observe that soft constraint approaches scale poorly 18

on problems with high ambient dimension, likely due to the high number of constraints that must be annealed. Across dimensions MASEM outperforms its NHR and OLLA counterparts for both problems. B.3

Robotics Problems Convergence Curves

We provide the convergence curves for the robotics problems. In Figure 8 we report the convergence for motion planning. In addition to the feasible entropy and the slack, we report the entropy of path homotopy classes. This corresponds approximately to the entropy of the component weights (for more details see Appendix D.1). In Figure 9 we report the convergence curves for the grasping problem. NHR

OLLA

SCMC

Cluster-NHR

Feas. Entropy ↑

Grid Obstacles

8 6 4 2

3.0

100

2.5

10−1

2.0

10−2

1.5

10−3

1.0

10−4

0.5

10−5

0.0 0

1000

2000

3000

Iteration

4000

5000

Feas. Entropy ↑

Random Obstacles

MASEM-NHR (Ours)

Homotopy Entropy ↑

1000

2000

3000

4000

Iteration

5000

Homotopy Entropy ↑

8 6 4 2

3.0

100

2.5

10−1

2.0

10−2

1.5

10−3

1.0

10−4

0.5

10−5

1000

2000

3000

Iteration

4000

5000

Slack ↓

1000

2000

3000

Iteration

4000

5000

4000

5000

Slack ↓

10−6

0.0 0

10−6

MASEM-OLLA (Ours)

1000

2000

3000

4000

Iteration

5000

1000

2000

3000

Iteration

Figure 8: Convergence curves on the motion planning problems with grid obstacles (first row) and random obstacles (second row). Curves depict mean and 95% CI across 5 seeds.

NHR

OLLA

SCMC

Cluster-NHR

MASEM-NHR (Ours)

Feas. Entropy ↑

MASEM-OLLA (Ours)

Slack ↓ 100

10−1

10−2

10−3 10−4

10−5

0 0

1000

2000

3000

4000

5000

1000

2000

3000

Iteration

4000

5000

Figure 9: Convergence on the grasping problem. Curves depict mean and 95% CI across 5 seeds.

Additional Experiments

C.1

Component Loss

We investigate the probability that after resampling no chains live within a specific component, which we call component loss. We analyze this phenomenon using a synthetic 2-dimensional benchmark, containing 100 disks in a 10 × 10 grid. Two neighboring disk-centers are 5 units apart; the disks have radius 1. We run MASEM-NHR for T = 10 iterations and count the number of components with at least one chain, averaged over ten runs. Importantly, the maximum step size of NHR is set 19

Uniform Initialization

1.0

0.5

0.1

± 3.3

± 6.6

± 2.2

± 0.8

± 0.6

± 1.3

± 4.1

± 3.1

± 2.1

± 2.2

42.6 77.2 96.4 99.9

100

99.9 99.9

100

± 4.1

± 0.0

± 0.6

± 0.0

43 ± 3.5

43 ± 6.4

± 5.4

± 3.1

± 0.6

75.8 93.1 98.8

100

± 4.4

± 0.0

± 4.9

± 2.3

72.5 90.3 96.9 99.7 99.9

100

± 4.2

± 0.6

± 0.0

99.5

100

99.8

100

± 1.8

± 1.3

± 0.0

± 0.8

± 0.0

± 2.2

± 4.7

± 0.9

62.7 89.8 95.3 97.8 98.2 ± 7.7

± 0.6

± 3.1

± 5.1

± 2.7

± 2.5

± 4.6

1.5

± 8.1

± 1.6

± 2.5

± 1.0

± 0.6

100

± 0.0

± 2.7

± 4.0

± 1.2

± 0.6

99.9

100

± 2.5

± 0.6

± 0.0

42.8 75.9 94.1 99.2

100

± 5.9

± 0.0

± 2.4

± 3.5

± 1.5

60.3 61.4 74.4 80.9 ± 10.8

100 200 300 400 500 600 700 800 900 1000 Number of Chains

± 1.3

± 5.0

± 4.7

0.5

± 2.0

42.4 80.7 98.8 99.9 99.9

42.7 77.2

1.0

0.1

100

42.1 82.9 97.5 98.5 99.3 99.5 99.6 99.6 99.9 99.9

2.0

± 4.1

± 5.5

± 6.5

± 5.0

± 5.9

± 3.9

90.3 91.4 93.4 ± 3.1

± 3.8

Covered Components

Temperature

1.5

41.9 77.7 97.3 99.8 99.9 99.4 96.8 97.5 97.9 98.7

Temperature

2.0

Worst-case Initialization

± 3.9

100 200 300 400 500 600 700 800 900 1000 Number of Chains

Figure 10: Number of covered components (out of 100), depending on temperature τ and number of chains N . Number include 95% confidence interval. to 1 to ensure that the chains only jump between components due to resampling. To initialize, we either distribute the chains uniformly across all components (uniform initialization), or populate every component except for one chosen at random starts with a single chain (worst case initialization). Figure 10 shows the results of this experiment. We see that for τ ∈ [0.5, 1.5], even just 4 times the amount of components reduces the chance of component loss over ten iterations to 0.1%. C.2

The Pairwise Distance KL-Divergence as Metric

In our reports, we use the (to our knowledge) novel metric of the KL-divergence computed on the histogram of pairwise distances. We define this metric in Appendix D.1. Here, we demonstrate the ability of this metric to distinguish between distributions with slightly different densities. To that end we sample from a disk (ρGT ), a disk with a higher weight in the center (ρcenter ) and a disk with higher weight at the border (ρedge ). For a divergence D we pose the null-hypothesis that the divergence does not distinguish between two instances of the ground truth and a non ground truth density and ground D truth. More precisely, our hypothesis Hcenter|edge is that D(ρGT , ρGT ) ∼ D(ρGT , ρcenter|edge ) holds.

Table 4: The table shows the probability that the null hypothesis “the metric cannot distinguish between a ground truth and a non-ground truth distribution” holds. It is computed with a Welsh t-test over 1000 runs. Note that only KL rejects the null-hypothesis with p < 0.01. D

D P(Hcenter )

D P(Hedge )

ED W22 KL

.068 .68 .013 .36 4.4 × 10−34 2.6 × 10−7

To measure this, we generate two sets SGT1 , SGT2 as well as Scenter and Sedge . Each set contains 2000 samples. For Scenter , we sample 50 points of those within a radius of 0.5 of the center instead of the ground truth. For Sedge , we instead sample 50 points with radius r ∈ [0.5, 1.0]. With these samples we calculate D(SGT1 , SGT2 ), D(SGT1 , Scenter ), and D(SGT1 , Sedge ). We repeat this process 1000 times and use the results to perform a Welsh t-test to D D compute the probability p of rejecting Hcenter and Hedge respectively. We compare the energy distance  12 X X X 2 1 1 ED(X, Y ) :=  ∥xi − yj ∥ − 2 ∥xi − xj ∥ − 2 ∥yi − yj ∥ , nm i,j n i,j m i,j 

used, for example, by Jeon et al. [9], the Sinkhorn distance W22 and our metric (KL). We report the results in Table 4. The KL divergence is the only metric which rejects the hypothesis that it cannot distinguish between ground truth and the other densities with p < 0.01. 20

Experimental Details

Each experiment is repeated across 5 different seeds. Where applicable, we report the mean across all 5 runs and 95% confidence intervals in our plots. Where reported we test statistical significance using t-tests and Holm-Bonferroni using SciPy. In the following, we provide details about the benchmarks and implementation that we used in this work. All the experiments are performed on a single workstation with a NVIDIA RTX PRO 6000 Blackwell GPU, a AMD Ryzen Threadripper PRO 7955WX CPU with 16 cores, and 128 gigabyte RAM. All algorithms, benchmarks and metrics are implemented in JAX [33].3 The 2d problems are depicted in Figure 6. We will provide mathematical and implementation details on the metrics, the benchmarks and their ground truth samplers (if available) below. After that, we provide details and hyperparameters for the method. D.1

Metrics

The main challenge in comparing different sampling methods is measuring the sampling quality, that is, how well the samples represent the actual distribution. We try to tackle this problem with multiple metrics. For all benchmarks where a ground truth is available, we report the Sinkhorn distance W22 . Additionally, we report the average maximum slack violation Ex [max{|h(x)|, g + (x)}] for all problems. For benchmarks with ground truth, we additionally report the KL divergence of pairwise distance histograms. If no ground truth is available, we instead report the feasible entropy. For the scaling benchmarks we calculate the TV distance between the empirical and the uniform component masses. In the path planning problem, we additionally report the entropy of the homotopy-class membership distribution. All these metrics are defined in more detail below. KL divergence. If a ground truth distribution is available, we report the KL divergence on the histogram of pairwise distances. This is well suited to detect small changes in component mass allocation (Appendix C.2). For this metric, we compute the pairwise distances of the ground truth distribution as well as the approximate distribution. With this, we compute the normalized histogram with 50 bins for values clipped between [0, b], where b = maxx,y∈Σ ∥x − y∥∞ , i. e. the maximum distance of feasible points along any dimension. We then calculate the forward KL divergence between the ground truth and the approximate distribution. Feasible Entropy If no ground truth is available, we report the feasible entropy, which is the estimated entropy of the feasible samples multiplied with the fraction of feasible samples. This metric serves as a proxy to measure whether we achieve our goal of entropy maximization. We compute P X 1[Slack(x) < tol] X Ĥfeas = Ĥ(Ci , k) , N k=2,4,8 Ci ⊆Cfeas

where Cfeas is the set of feasible samples, Ci contains 100 samples drawn without replacement from Cfeas for i = 1, . . . 10 and Ĥ(Ci , k) denotes the Kochzachenko-Leonenko Entropy estimate of Ci , estimated using the k-th nearest neighbor [34]. If less then 100 feasible samples are available, we do not compute this metric. We choose to subsample 100 points to compare entropy across methods with potentially different numbers of feasible samples, since the estimator depends on the number of samples. TV-Distance We introduce a new metric for the scaling problems, since the other metrics might struggle with high-dimensional samples. For this metric we count the number of samples contained in each component (up to a small tolerance of 10−5 ) and normalize it. This yields the component weights αc as introduced in Section 5. We then compute the total variation distance to the ground truth αc∗ as TV = max |αc − αc∗ | c=1,...,C

Homotopy-Entropy. For the motion-planning problems, we additionally measure the entropy of the distribution of paths among their homotopy-classes. Two paths belong to different homotopy 3We will provide the code upon publication.

classes if they cross an obstacle on different sides. This directly corresponds to the components of S the manifold. More formally, we consider the space W = ([−4, 4]2 \ O)/ ∼, where O = i Oi is the union of obstacles and (xi , yi ) ∼ (xj , yj ) iff xi = xj = 3.6, which is the vertical goal line. We then consider homotopy classes of paths on this space. Note that, since we constrain paths to have increasing x-coordinates, we only observe a finite subset of the hompotopy classes. In practice, we compute a signature for each path by testing whether it passes above or below for each obstacle. This results in a |obs| dimensional vector. We count the occurrences of each signature and compute the entropy of their distribution. Note that while this is a useful metric, maximizing it does not directly correspond to maximizing the entropy of paths. By computing the entropy over the signatures, we implicitly assume that each component has the same weight. While not true in practice, this is still a useful proxy, since the actual component weights are unknown a priori. D.2

Synthetic Problems

For each low-dimensional example, we run 2 000 independent chains for 5 000 steps. From each chain, we retain only the last state, yielding 2 000 samples. Connected & Disconnected Disks We consider two randomly rotated disks embedded on a sphere in R3 . h(x) = ∥x∥2 − R, R = 2.5; g(x) = min {cos(ρi ) − R⟨x, µi ⟩} . i∈{1,2}

We use two different parameters ρ1 = 0.2, ρ2 = 0.6 such that both disks have different areas. In the connected case the disks have a arc-length of δ = 0.6 between their centers; in the disconnected case we use δ = 1.35. We choose the disk centers " # T √ 1 0 0 µ1 = Rx cos − 2δ , sin − 2δ , 0 2 where Rx = 0 t t ∈ SO(3), t = T 2 µ2 = Rx cos 2δ , sin 2δ , 0 0 t t is the rotation of 45 degrees around the x-axis. In addition we impose bound constraints x ∈ [−5, 5]3 . Ground-truth samples are generated by exact surface-uniform sampling from the union of the two disks, with a correction for possible overlap. We first choose disk i ∈ {1, 2} with probability i) proportional to its spherical area, P(i) = P2 1−cos(ρ . Given disk i, we sample uniformly on its (1−cos(ρ )) j

j=1

area by drawing samples from ϕ ∼ Unif(0, 2π),

cos α ∼ Unif(cos ρi , 1),

and setting x = R [cos α µi + sin α (cos ϕ ei,1 + sin ϕ ei,2 )]. To correct for overlap, proposals are accepted with probability 1/m(x), where m(x) is the number of disks containing x. Seven Lobes We adapt this problem from Jeon et al. [9] and replace the Gaussian mixture target with constant target f (x) = 0 to sample uniformly. q h(x) = x21 + x22 − (3 + cos(7θ)) , θ = atan2(x2 , x1 ); 1 g(x) = (x1 − 2)2 − 5x1 x32 + x52 − 40. 2 Additionally, we impose bound constraints such that x ∈ [−4.1, 4.1]2 . We use rejection sampling from the polar parametrization of the equality constraint to generate ground truth samples. Specifically, we sample θ ∼ Unif(0, 2π),

r(θ) = 3 + cos(7θ),

x(θ) = r(θ)(cos θ, sin θ).

We correct the acceptance rate with a Jacobian correction that accepts samples with probability proportional to s 2 ds dr dr = r(θ)2 + , = −7 sin(7θ). dθ dθ dθ 22

Sine We define a uniform target supported on a one-dimensional sine curve in R2 . The constraints are h(x) = x2 − exp(−0.15x1 ) sin(x1 ); x2 ≥ 0. We additionally impose bound constraints such that x ∈ [−20, 20]2 . We generate ground truth samples using rejection sampling from the parametrization of the equality constraint similar to the seven lobes procedure above. Specifically, we sample t ∼ Unif(−20, 20),

x(t) = (t, exp(−0.15t) sin t) .

Since uniform sampling in t is not uniform w.r.t. arclength, we correct the acceptance probability with the Jacobian. We retain only accepted proposals satisfying the halfspace and bound constraints. Swiss Roll We define a uniform target on a disconnected one-dimensional manifold in R2 consisting of six circles and one Archimedean spiral segment. The circle components are hk (x) = ∥x − ck ∥2 − rk ,

k = 1, . . . , 6,

where the centers ck and radii rk are sampled once with layout seed 0, using non-overlapping circles with radii in [0.45, 1.15]. The spiral is parametrized by γ(t) = csp + (a + bt)(cos t, sin t),

t ∈ [tmin , tmax ],

with csp = (4.8, −0.4),

a = 0.45,

b = 0.33,

tmin = 0.9,

tmax = 3.9π.

The equality constraint is the signed residual of the closest circle or valid spiral branch. We additionally impose g(x) = softplus(−3 − x1 ) softplus(−x2 ) − 0.1, as well as box constraints x ∈ [−10, 10]2 . Ground-truth samples are generated by first choosing a circle or the spiral with probability proportional to its parameter-domain length, i.e. 2π for each circle and tmax − tmin for the spiral. For circles we sample ϕ ∼ Unif(0, 2π) and set x = ck + rk (cos ϕ, sin ϕ), correcting by the arclength p Jacobian. For the spiral we sample t ∼ Unif(tmin , tmax ) and set x = γ(t), correcting by ds = (a + bt)2 + b2 . dt Scaling Stress Test In the synthetic stress-test benchmark we target an uniform disconnected manifold embedded in Rd . We use m to denote the total number of equality constraints. Of these, m − 1 are linear equality constraints, which define a random subspace via Ax = b, where A ∈ R(m−1)×d has orthonormal rows and is sampled randomly once. In all stress-test experiments, we set b = 0. Let p = d − (m − 1) = d − m + 1 denote the dimension of the nullspace of A, and let N ∈ Rd×p be an orthonormal basis of this nullspace. We write the corresponding latent coordinates as z = N ⊤ x. Within this subspace we then randomly sample |C| = 5 disconnected hyperspheres which define the manifold on which we sample via equality constraints. For centers ci ∈ Rp and radii ri , we use 2

h(x) = ∥z − ci⋆ ∥2 − ri2⋆ ,

i⋆ = arg min ∥z − ci ∥2 − ri2 . i∈{1,...,5}

Since each sphere has dimension p − 1, the intrinsic manifold dimension is p − 1 = d − m. The base radius is set to 0.25, and the component radii are multiplied by independent jitter factors in [0.5, 1.5]. The component centers are placed randomly in latent space using rejection sampling to ensure separation between components. 23

We additionally include exactly l = 5 spherical inequality constraints, which act as ball cutouts on the sphere components. Each cutout center is placed on the surface of a randomly selected component, and its radius is set to half of the corresponding component radius. For cutout center oj and radius ρj , the inequality is 2 gj (x) = ρ2j − ∥z − oj ∥2 , j = 1, . . . , 5, so that feasible samples lie outside all cutout balls. Finally, we impose box constraints x ∈ [−28, 28]d . We compare two main scaling regimes: increasing ambient dimension with fixed manifold dimension for m = d − 3 and increasing ambient dimension with increasing manifold dimension for fixed m = 5. In this case, the intrinsic dimension grows with the ambient dimension and is given by d − m = d − 5. Ground-truth samples are generated directly from the latent representation. We first sample a component index with probability proportional to its surface area, rp−1 P(i) = P5 i p−1 . j=1 rj Given the component, we sample a direction u uniformly on the unit sphere Sp−1 by normalizing a standard Gaussian vector, and set z = ci + ri u,

x = N z.

Since all stress-test experiments include ball cutouts, we use rejection sampling from this proposal and retain only feasible samples. D.3

Motion Planning

We consider 2D motion planning problems in which a point robot must navigate from a fixed start configuration q0 ∈ R2 to a goal region while avoiding circular obstacles. We parametrize each path through nw waypoints X = (x1 , . . . , xnw ) ∈ Rnw ×2 , which form the sampling domain. Following standard robotics convention, we recover smooth trajectories as cubic B-splines whose control points are obtained from the start and waypoints via C = Φ−1 [q0 ; X], using the pseudoinverse of the interpolation matrix Φ ∈ RT ×(nw +1) . Evaluating the resulting curve at T equally spaced parameter values yields a discretized path p ∈ RT ×2 . The goal is a vertical line on the right side of the workspace, which we impose as an equality constraint on the first coordinate of the final waypoint: hgoal (X) = xnw ,1 − 3.6, leaving the second coordinate free so any point on the goal line is admissible. The main challenge of this domain is navigating around obstacles, where multiple homotopically distinct paths are admissible. Collision constraints are enforced at each path sample via coll gj,t (p) = rj2 − ∥pt − cj ∥22 ,

t = 1, . . . , T,

for each obstacle j with center cj and radius rj . We additionally impose monotone progress along the x axis, gtmon (p) = pt−1,1 − pt,1 , t = 1, . . . , T, with p0 ≡ q0 , which removes looping and backwards-going trajectories. We further impose discrete velocity and acceleration limits to suppress paths with large gaps or sharp kinks, both of which are known failure modes when the path is reconstructed from sparse waypoints: gtvel (p) = ∥pt − pt−1 ∥22 − s2max ,

smax = 1,

gtacc (p) = ∥pt+1 − 2pt + pt−1 ∥22 − a2max ,

amax = 0.65.

Finally, all waypoint coordinates are bound-constrained to [−4, 4]2 . Since each obstacle splits the configuration space into a homotopy class passing above and one passing below, the feasible set decomposes into up to 2|obs| connected components (modulo monotonicity and feasibility pruning), which we use to evaluate homotopy-class entropy. For more details see Appendix D.1. 24

Grid Layout. In the grid problem, obstacles are placed on a regular 4 × 4 grid, that is centered within the domain. All obstacles share a common radius r = 0.5. In total, this induces m = 3 equality and l = 771 inequality constraints. This setup yields a structured benchmark with a known and combinatorially large number of homotopy classes, which is well suited for testing whether a sampler covers all feasible passages between obstacle columns or collapses onto a few preferred routes. Random Layout. In the random problem, we place 20 disk obstacles with radii drawn independently from Unif[0.2, 0.5] and center coordinates drawn uniformly over the domain. We employ rejection sampling to guarantee minimal separation between obstacle pairs to increase the number of feasible paths. In total, this induces m = 3 equality and l = 931 inequality constraints. The random layout induces an irregular and often narrower set of feasible corridors than the grid, reflecting real world robot navigation challenges. D.4

Grasping

We want to grasp a capsule, which is a cylinder with a halfsphere with the same radius as the cylinder on each end. This shape is commonly used for collision checking. The capsule is placed at the origin (0, 0, 0). Its orientation is given by a vector along its long axis, pointing towards o = (1, 1, 1). It has the length l = 1 and the radius r = 0.25. Gravity impacts this capsule with the force-vector fg = (0, 0, 1). Its friction coefficient is µ = 1. Figure 11 shows the capsule and a grasp schematically. We grasp the capsule with three fingers. The fingers as hard pointfingers with friction [35]. This means each finger i is parametrized by a point of attack pi ∈ R3 as well as a force-vector fi ∈ R3 . To ensure the capsule is grasped and stable, the sum of forces of the fingers has to counter gravity, and the sum of torques at it’s center τi = pi × fi has to be zero. This yields the constraints f1 + f2 + f3 = fg τ1 + τ2 + τ3 = 0 . For each finger i, we calculate a vector pointing from the middle axis of the capsule towards the finger as Figure 11: Schematic diagram visualizing the capsule in 2d. di = clamp ⟨o, pi ⟩, −l/2, l/2 o − pi . The red band in the middle repThis allows us to ensure that the point of attack lies on the capsules resents the infeasible region. surface by constraining ∥di ∥ = r. Furthermore, each finger can In the top left, a finger conexert a force deviating from the normal of the capsules at the point tact is shown. The orange dot of attack ni = di /r as far as friction allows. This is expressed in the is the point of attack on the force cone constraint capsule surface. The orange cone shows the friction cone ∥fi − ni ⟨ni , fi ⟩∥ ≤ µ⟨−ni , fi ⟩ . and the black arrow the surWe constrain the force each finger can apply to ∥fi ∥ ∈ [0.1, 1]. face normal. The dotted black Lastly, we introduce an infeasible region on the capsule. This is arrow in the middle symbolmeant to model objects that can’t be grasped at arbitrary points. We izes gravity. choose a band across the middle of the cylinder with width w = 0.25. This results in the constraint |⟨o, pi ⟩|≥ w/2. For each finger, we have 6 dimensions, 1 equality, and 4 inequality constraints. We enforce 6 equalities in the force and torgue constraints. Additionally, we constrain each entry of the points of attacks and forces to be between −1 and 1. Overall, this results in 18 dimensions, 9 equalities and 48 inequalities. D.5

Implementation Details

All methods are implemented in JAX [33]. We implement all baselines based on their official implementation: NHR is based on the code of Toussaint et al. [24], OLLA based on Jeon et al. [9] and SCMC on Golchi and Campbell [17]. For OLLA, we use the OLLA-H variant which employs Hutchinson’s trace estimator. 25

Initialization. We initialize all sampling runs from a Gaussian in ambient space which we center at the origin and scale by a quarter of the distance between minimum and maximum bound for that problem. Following Toussaint et al. [24], we implement the initial projection for NHR and MASEM-OLLA using the Gauss-Newton method for 500 steps with additional standard Gaussian noise that we scale by ϵ = 0.01. While this is inconsistent with theory, we find that it improves mode coverage of the initialization and yields feasible samples on every problem that we tested on. MASEM Details. In Section 5, we assume that the local sampling kernel preserves sample feasibility. Since practical samplers tend to break this property resampling naively can increase the number of infeasible samples as those tend to have largest distance to the samples on the manifold. To prevent this from happening, we apply a penalty from the resampling weights: wi ← wi

1 exp µ Slack(xi )

We choose µ = 1000. as default and find that this offers a robust choice. Further, we ensemble multiple values for k to compute the weights [36, 37]. Specifically for k = 3 we would compute all Pk distances εi,1 , εi,2 , εi,3 and take the average over their resulting kNN estimators w̄i = k1 j=1 εi,j . In practice, we find that this approach tends to have lower variance and is less sensitive to the exact choice of k despite introducing bias. Hyperparameters We tune the main hyperparameters of each method and baseline using Bayesian Optimization [38], implemented in the wandb api.4 For each method we allocate the same budget of 20 tuning runs which we tune against the Sinkhorn distance on problems with ground truth samples and against the feasible entropy on practical problems. The MASEM hyperparameters for all problems are listed in Table 5. Table 5: Hyperparameters for MASEM-NHR and MASEM-OLLA. Step sizes and other samplerrelated hyperparameters are taken from the corresponding base sampler. Hyperparameters are the temperature τ , the number of mixing steps M and the considered neighbour k. MASEM-NHR MASEM-OLLA Problem

Connected Disks Disconnected Disks Seven Lobes Sine Swiss Roll Stress Test (m = d − 3) Stress Test (m = 5)

1.0 50 1.0 50 .81 5 .75 5 .65 5 .3 5 .74 20

4 4 8 16 16 16 20

1.0 1.0 .67 .98 .5 .3 .3

50 50 5 10 50 5 5

4 4 8 16 8 16 16

Motion Planning Grid Motion Planning Random Grasping

.65 50 1.0 50 .41 5

8 8 19

.65 50 .935 50 .94 10

8 4 16

Limitations

E.1

Theoretical Limitations

Some results, namely Lemma 2, Theorem 1, and Corollary 1, are valid only in the large-scale particle limit, i.e. for N → ∞. Furthermore, the main result in Theorem 1 is only stated in the mean-field. It might happen that no particles belonging to a particular component are sampled, resulting in that component vanishing. However, the probability of this occurrence decreases exponentially with the number of samples and our experiments show this effect is negligible even for just 4 chains per component, as we demonstrated in Appendix C.1. 4 https://docs.wandb.ai/models/sweeps

E.2

Runtime Overhead

While small, our method introduces some runtime overhead compared to only using constrained samplers like NHR: every M steps, we compute the k-nearest neighbor distances and resample. This is also reflected in the asymptotic runtime. Assuming each iteration of the constrained sampler is constant in the number of samples, we increase the asymptotic runtime from O(N T M ) to O(N T (M + log N ) when using the proposed resampling.5 In practice however, we discovered that on our benchmarks this overhead is small, while applying MASEM greatly enhances sample quality. E.3

Non-Uniform Sampling

This work focuses on entropy maximization, which results in uniformly sampling the feasible set. This limits it’s usefulness when non-uniform samples are needed. If that is the case, we suggest either applying MASEM and using rejection sampling afterwards, or using the output of MASEM as an unbiased initialization for a Sequential Monte Carlo sampler [8]. E.4

Parameter Selection

As shown in Theorem 1 and Corollary 1, the convergence of the algorithm depends on the temperature τ , which should be chosen to be smaller than the intrinsic dimension of the manifold. For many applications, the intrinsic dimension is not known a priori, making it harder to select a suitable parameter. We suggest a few mitigation strategies. First, if the available compute resources allow, one could test multiple different parameters and compare the resulting entropy, as we did for tuning the robotics benchmarks. Secondly, especially in robotics, it is feasible to obtain a rough estimate of the intrinsic dimension by reasoning over the degrees of freedom of the problem. This could then be used as an upper bound. Lastly, if both strategies are infeasible, we suggest using a conservative temperature (for example, τ = 0.9) and increasing the number of iterations T .

Broader Impact

This work is foundational, and we expect no direct societal impact. Nevertheless, as described in Section 1, constrained sampling has many applications in robotics. An important concern when deploying robots in the real world is safety, especially when they are controlled by a neural network. We believe safe operation can be encouraged during training by sampling from a constrained distribution, where the constraints ensure, for example, that the robot does not collide with the environment. The method presented here should be only one of many components of a holistic system for safe robot operation. While we hope it can play a small part in improving robotic operations, we believe that its direct broader impact is limited.

List of Acronyms

CGF LICQ MASEM MCMC NHR OLLA SCMC SMC SVGD

Cumulant Generating Function Linear Independent Constraint Qualifier (see assumption 1) Manifold Sampling via Entropy Maximization (introduced in this paper) Markov Chain Monte Carlo Non-linear hit & run (constrained sampler, introduced in [24]) Overdamped Langevin with Landing (constrained sampler, introduced in [9]) Sequentially Constrained Monte Carlo (constrained sampler, introduced in [17]) Sequential Monte Carlo Stein Variational Gradient Descent (introduced in [40])

5 Using approximate nearest neighbors, one can determine the k-nearest neighbor distance in O(N log N ) [39]. Further-

more, with binary search, the resampling step can be executed in O(N log N ) as well

Record · ID 178899 · SHA-256 6153d78d95d597d7
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.