Sampling with Riemannian Hamiltonian Monte Carlo in a Constrained Space
Authors: Yunbum Kook
YinTat Lee
Ruoqi Shen
Santosh Vempala
Venue: NeurIPS 2022 Accept (2022)
Keywords: Sampling
Hamiltonian Monte Carlo

Abstract: We demonstrate for the first time that ill-conditioned, non-smooth, constrained distributions in very high dimension, upwards of 100,000, can be sampled efficiently \emph{in practice}. Our algorithm incorporates constraints into the Riemannian version of Hamiltonian Monte Carlo and maintains sparsity. This allows us to achieve a mixing rate independent of smoothness and condition numbers. On benchmark data sets in systems biology and linear programming, our algorithm outperforms existing packages by orders of magnitude. In particular, we achieve a 1,000-fold speed-up for sampling from the largest published human metabolic network (RECON3D). Our package has been incorporated into a popular Bioinformatics library.

Sampling with Riemannian Hamiltonian Monte Carlo in a Constrained Space
Yunbum Kook
Yin Tat Lee Microsoft Research, University of Washington yintat@uw edu
Georgia Tech yb. wwo.ooattecc.....o edu
Ruoqi Shen University of Washington shenr3@es www.ingtosssinnnon..noooooo edu
Santosh S. Vempala Georgia Tech vempala@gatech..om edu
Abstract
We demonstrate for the first time that ill-conditioned, non-smooth, constrained dis- tributions in very high dimension, upwards of 100,000, can be sampled efficiently in pracfice. Our algorithm incorporates constraints into the Riemannian version of Hamiltonian Monte Carlo and maintains sparsity. This allows us to achieve a mixing rate independent of condition numbers.
On benchmark data from systems biology and linear programming, al- gorithm outperforms existing packages by orders of magnitude. In particular, we achieve 1,000-fold speed-up for sampling from the largest published human metabolic network (RECON3D). Our package has been incorporated into the COBRA toolbox.
Introduction
Sampling Fundamental. Sampling algorithms arise naturally in models of statistical physics, e.g., Ising, Potts models for magnetism, Gibbs model for gases, etc. These models directly suggest Markov chain algorithms for sampling the corresponding configurations. the vertices of a graph are assigned spin, i.e., +1, in each step, we pick vertex at random and flip its spin with some probability. The probability chosen so that the distribution of the vector of all spins approaches target distribution where the probability exponentially decays with the number of agreements in spin for pairs corresponding to edges of the graph. the Gibbs model, particles move randomly with collisions and their motion is often modeled reflecting Brownian motion. Sampling with Markov chains today the primary algorithmic approach for high-dimensional sampling. For some fundamental problems, sampling with Markov chains is the only known efficient approach or the only approach to have guarantees of efficiency. Two notable examples are sampling perfect matchings of bipartite graph and sampling points from convex body. These are the core subroutines for estimating the permanent of a nonnegative matrix and estimating the volume of convex body, respectively. The solution space for these problems scales exponentially with the dimension. In spite of this, polynomial-time algorithms have been discovered for both problems. The current best permanent algorithm scales as n² (time) 23], while the current best volume algorithm scales as n (number of membership tests) [24]. For the latter, the first polynomial-trn algorithm had a complexity 2? [16], and the current best complexity is the result of many breakthrough discoveries, including general- purpose algorithms and analysis tools.
the Ising model where
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Sampling is Ubiquitous. The need for efficient high- -igh-siimnsional sampling arises in many fields. notable setting is metabolic networks in systems biology. A constraint- based model of metabolic network consists of metabolites and n reactions, and a set of equalities and inequalities that define a set of feasible steady state reaction rates (fluxes):
Ω = {u E R" Se - 0, S a
where a stoichiometric matrix with coefficients for each metabolite and reaction. The linear equalities ensure that the fluxes into and out of every node are balanced. The incqualities arise from thermodynamical and environmental constraints. Sampling constraint- based models is a powerful tool for evaluating the metabolic capabilities of biochemical networks [33, 46]. While the most common distribution used is uniform over the feasible region, researchers have also argued for sampling from the Gaussian density restricted the feasible region; the latter has the advantage that the feasible set does not have to be bounded. A previous approach to sampling, using hit-and-run with rounding [20], has been incorporated into the COBRA package [21] for metaholic systems analysis (Bioinformatics).
second example of mathematical interest is the problem of computing the volume of the Birkhoff polytope. For a given dimension n, the Birkhoff polytope is the set of all doubly stochastic matrices (or the convex hull of all permutation matrices). This object plays a prominent role alge- braic geometry, probability, and other fields. Computing its volume has been pursued using algebraic representations; however exact computations become intractable even for n = 11, requiring years of computation time. Hit-and-run has been used to show that sampling based volume computation can go to higher dimension [11], with small error of estimation. However, with existing sampling implementations, going beyond n = 20 seems prohibitively expensive.
third example is from machine learning, a field that increasingly turning to sampling models of data according to their performance in some objective. One such commonly used criterion is the logistic regression function. The popularity of logistic regression has led to sampling being incorporated into widely used packages such STAN [44], PyMC3 [41], and Pyro [3]. However, those packages in general do not run on the constraint-basee models we are interested in.
Problem Description. In this paper, we consider the problem of sampling from distributions whose densities are of the form
- subject to Az x E K
(1.1)
where f is a convex function and K convex body. We assume that a self-concordant barrier for K is given. Note that any convex body has self-concordant barrier [32] and there are explicit barriers for convex bodies that come up in practical applications [37], so this is a mild assumption. We introduce an efficient algorithm for the problem when K is product of convex bodies K;, each with small dimension. Many practical instances can be written this form. As special case, the algorithm can handle K in the form of {z E R" : 4, S x; S u; for all i E [n]] with 1 E RR U {-oo} and in E TR U (+20). which is the common model structure in systems biology. Moreover, any generalized linear model exp(- Eiaa, x bx)), e.g., the logistic model, can be rewritten in the form
exp(- -po0-yu ki) subject to Au = b + 8, (8, () E K
(1.2)
where K = ITK, and each K; ((s,, 14:: fi(::) S 4 is a two-dimensional convex body.
The Challenges of Practical Sampling. High dimensional sampling has been widely studied in both the theoretical computer science and the statisties communities. Many popular samplers are first-order methods, such as MALA [40], basic HMC [36, 14] and NUTS [22]. which update the Markov chain based on the gradient information of f. The runtime of such methods can depend on the condition number of the function f [15, 30, 7, 8, 42]. However, the condition number of real-world applications can be very large. For example, RECONI [27], a reconstruction of the human metabolic network, can have condition number as large as 10t due to the dramatically different orders of different chemicals' concentrations. Motivated by sampling from ill-conditioned distributions, another class of samplers use higher-order information such as Hessian of take into account the local structure of the problems [43, 9]. However, such samplers cannot handle non-smooth distributions, such as hinge-loss, lasso, uniform densities over polytopes.
For non-smooth distributions, the best polytime methods are based on discretizations of Brownian motion, e.g., the Ball walk [25] (and its
cousin, the Dikin walk [26]), which takes

---

random step in ball of a fixed size around the current point. Hit-and-Run Run [34] builds on these by avoiding an explicit step size and going to a random point along random line through the current in polytope that contains unit ball, the step size
point. Both approaches hit the same bottleneck should be (1)/ to avoid stepping out of the body with large probability. This leads to quadratic bounds (in dimension) on the number of steps to "mix".
Due the reduction mentioned in (1.2), non-smooth distributions can be translated to the form in (I.1) with constraint K. Both the first and higher-order sampler and the polytime non-smooth samplers have their limitations in handling distributions with non- non-smooth objective function or constraint K. Given the limitations of all previous samplers, natural question we want to ask is the following.
Question. Can we develop a practically efficient sampler that can handle the constrained problem in (1.1) and preserve sparsity' with mixing time independent of the condition number?
In some applications, smoothness and condition number can be controlled with tailor-made models. Our goal here is propose general solver that can sample from any non-smooth distributions as given. For traditional samplers such as the Ball walk and Hit-and-Rum Run, as mentioned earlier, the step size needs to be small so that the process does not step out. An approach that gets around this bottleneck is Hamiltonian Monte Carlo (HMC), where the next step is given by a point along a Hamiltonian- -preeerving curve according to suitably chosen Hamiltonian. It has two advantages. First, the steps are no longer straight lines in Euclidean space, and we no longer have the concern of "stepping out". Second, the process is symplectic (so measure-presse preeeeriipr.sr..gg....... and hence the filtering step is easy to compute. It was shown in [31] that significantly longer steps can be taken and the process with convergence analysis in the setting of Hessian manifolds, leading subquadratic convergence for uniformly sampling polytopes.
To make this practical, however, is formidable challenge. There are two high-level difficulties. One is that many real-world instances are highly skewed (far from isotropic) and hence it is important to use the local geometry of the density function. This means efficiently computing or maintaining second-order information such as Hessian of the logarithm of the density. This can be done in the Riemannian HMC (RHMC) framework [17, 31 but the computation of the next step requires solving the Hamiltonian ODE to high accuracy, which in turn needs the computation of leverage scores, a procedure that takes least matrix caatri-muutiplicatton time in the worst case. Another important difficulty is maintaining hard linear constraints. Existing high-dimensional packages do not allow for constraints (they must be somehow incorporated into the target density). ) and RHMC is usually considered with a dile-dimensional feasible region such as full-dimensional polytope. This can also be done in the presence of linear equalities by working in the affine subspace defined by the equalities, but this has the effect of losing any sparsity inherent in the problem and turning all coefficient matrices and objective coefficients into dense objects, thereby potentially incurring a quadratic blow-up.
Our Solution: Constrained Riemannian Hamiltonian Monte Carlo (CRHMC). We develop a constrained version of RHMC, maintaining both sparsity and constraints. Our refinement of RHMC ensures that the process satisfies the given constraints throughout, without incurring a significant overhead in time sparsity. It works even if the resulting feasible region is poorly conditioned. Since many instances in practice are ill-conditioned and have degeneracies, we believe this is a crucial aspect. Our algorithm outperforms existing packages by orders of magnitude.
we give the main ingredients of the algorithm and discuss how we overcome the challenges that prevent us from sampling efficiently in practice. Following that, in Section 3, we present empirical results on several benchmark datasets, showing that CRHMC successfully samples much larger models than previously known to be possible, and is significantly faster in terms of rate of convergence ("number of steps") and total sampling time. Our complete package is available on GitHub. We refer the reader to Appendix for theory, notations, and definitions.
In Section 2,
When A sparse. preserving the sparsity of A can greatly enhance both the runtime and the space efficiency.
3

---

2 Algorithm: Constrained RHMC
In this section, we propose a constrained Riemannian Hamiltonian Monte Carlo (CRHMC)) algorithm to sample from distributions of the form
(() subject c(a:) and K for some convex body K,
where the constraint function c: R + _ satisfies the property that the Jacobian De(a:) has full rank for x such that c(x) 0. useful keep in mind the case when c(x) = is an affine subspace Ar = b, in which case De(a) A, and the full-rank condition simply says that the rows of are independent.
We refer readers to [1. 4, 39] for preliminary versions of CRHMC called the constrained Hamiltonian Monte Carlo (CHMC). In particular, framework [4] can be extended CRHMC when K R" and in fact they mention CRHMC as possible variant. However, their algorithm for CRHMC requires eigenvalue decomposition and is not efficient for large problems, which takes n time and n² space per MCMC step in practice. In this section, we propose an algorithm that overcomes those limitations and satisfies the additional constraint K by using local metric induced by the Hessian of self-concordant barriers, leading to n re..5 time and space in practice.
2.1 Basics of CRHMC
To introduce our algorithm, we first recall the RHMC algorithm (Algorithm 1). RHMC, extend the space x the pair v), where denotes the velocity. Instead of sampling from e e-((x), RHMC samples from the distribution e-H(2,)) where H (ar, v) is the Hamiltonian, and then outputs To make sure the distribution is correct, choose the Hamiltonian such that the marginal of e (a.,)) along is proportional to e e-f()), One common choice of (ar, v) is
H v) = f(x) + 2" ((x) log det (x),,
(2.1)
2
where (x) is position-dependent positive definite matrix defined on R".
Algorithm 1: Rimaniaan Hamiltonian Monte Carlo (RHMC)
Input: Initial point (0), step size
for k 1, do // Step 1: resample
Sample ((k +) ~ N'0 M(a(k-11)) and set (N- (k-11)
// Step 2: Hamiltonian dynamics
Solve the ODE
dr OH(x,)) du
OH (x v)
(2.2)
dt dv dt Ox with H defined in (2.1) and the initial point given by (I u 21)
Set (R)) < (h) and (k) + v(h).
end Output: (A)
To extend RHMC to the constrained case, we need to make sure both Step and Step satisfy constraints, the Hamiltonian dynamics has to maintain c(x: 0 throughout Step 2. Note that
d
dr: De(xx)) (H
di c((,) De(xx) di
(2.3)
Ov:
where De(œ) is the Jacobian of at x. With H defined in (2.1). Condition (2.3) be- comes (e(=) 0. However, for full rank De(x), if M(x) is invertible. then Range(u) (OO g" immediately violates this condition due to
"rooooucced "crumch".

---

n m. To get around this issue, we use non-invertible matrix (ax with its pseudo-inverse Mx) to satisfy (x) = 0 for any v E Range( (MggeM((a) (x)). Since we want the step to be able to move in all directions satisfying c(x) we impose the following condition with Range(M (x)))(M(x))) in mind:
Range
Null(Dc(z)) for all R",
(2.4)
which can be achieved by (x) proposed soon.
Under the condition (2.4), we sample from ( Step 1, which is equivalent to sampling from (z,() subject to v & Range(M(a)) = Null( Also, the stationary distribution of CRHMC should be proportional to
e --H(z,w) subject to c(x) and v E Null(De(2)).
Here, to maintain e E Null(De(z)) during Step 2 we add a Lagrangian term to H. Without the Lagrangian term, would escape from Null(De(2t.)
Range( M (zx)) in Step 2 as seen in the proof of Lemma 1, which contradicts Range(U,) = Range(N M (x))) = Range(M The constrained Hamiltonian we propose is (See its rigorous derivation in Lemma
1
H(x, = = A(x, T c(r)
with H(r,v)= f()) + z
M(x) log (dee(M((
(2.5)
where (x (Dc(=)Dd(=))>- (Doc())u,, (DD))D #1 - DD(xx) 0H
0= )/.")) Here, pdet denotes
pseudo-determinant and A(x, v) is picked so that E Null(De(x)). An algorithmic description of CRHMC is the same as Algorithm with the constrained H in place of the unconstrained H. We show the convergence of CRHMC to the correct distribution
) in Appendix B.3.
Choice of . via Self-concordant Barriers. The construction of the Hamiltonian (2.5) relies on having a family of positive semi-definite matrix M (I) satisfying the condition (2.4) (i.e.. Range(M(x)) = Null Nul((D((r))). One natural choice is the orthogonal projection to Null(Dc(z)):
Q(x) I De(=) (Dd))
(2.6)
which is similar to the choice in [4].
For the problem we care about, there are additional constraints on x other than (c(x) = 0}. In the standard HMC algorithm, we have dx ~ (0, (x)))... For example, for a simple constraint de
K 10. to ensure every direction is moving towards/away from x 0 multiplicatively. a natural choice of M is M(x) = diag(= 2)). For general convex body K, we can use a celf-concordant barrier, function defined on K such that (xx) is self-concordant and Φ(x) + +00 as æ ŌK. Using the barrier d, we can define the local metric based on g(x) V²ő(x). Intuitively, as the sampler approaches 0K, the local metric stretches accordingly that the Hamiltonian dynamics never passes the barrier, respecting x E K throughout.
In summary, we need (x) to have its range match the null space of De(x) and agree with g(x) in its range. We can verify that M(r)= Q(x) g(x))(x)), where Q(x) is the symmetric matrix defined in (2.6), satisfies these two constraints.
2.2 Efficient Computation of aH /rr and 0H/Ov
With M Q(x) ((x)((x), we have all the pieces of the algorithm. However, using this naive algorithm to compute 0H/Oz dz and 0H/Oc, we face several challenges.
1. The algorithm involves computing the pscudo-inverse and its derivatives, which takes O(z²) except for very special matrices.
2. The Lagrangian term in the constrained Hamiltonian dynamics requires additional computation such as the second-orda derivative of c(z).
naive approach to computing leverage scores 0H results in very dense matrix.
Those challenges make the algorithm hard to implement and inefficient, especially when the dimension high. In the following paragraphs, we give an overview of how we overcome each the challenges above. We defer a more detailed discussion of our approaches and the proofs to Appendix B.2.
5

---

Avoiding Pseudo-inverse and Pseudo-determinant.. We are able to show equivalent formulas for (x) and log pdet.M (z) that can take advantage of sparse linear system solvers. In particular, we show that Mxx² = g(x) $ - - P(x)) g(x) &, where
P(æ) - g s(x) Dc(x) g(z)- De(x)))-\D((a)) g(x))=
(2.7)
As mentioned earlier, a majority of convex bodies appearing in practice are of the form K= = II: K, where K; are constant dimensional convex bodies. In this case, we will choose g(x:) to be block diagonal matrix with each block of size (11). Hence, the bottleneck of applying P(x) to a vector is simply solving a linear system the form (De 9 g- De )u = b for some b. The existing sparse linear system solvers can solve large classes of sparse linear system much faster than (nn)) time [13]. For log pdet. (z), we show
log pdee (M)) log det ((x)) Hog det (De(x) g(x))- -D(x))) det (De(s). De(a) T)_
(2.8)
This simplification allows to take advantage of sparse Cholesky decomposition. We prove (2.7) and (2.8) in Lemma 2 and Lemma 3 in Appendix B.2.1. The formulas (2.7) and (2.8) avoid the expensive pseudo-inverse and pseudo-determinant computations, and significantly improve the practical performance of our algorithm.
Simplification for Subspace Constraints. For the case c(x) = Ax - b, the Hamiltonian is now
1 T g & (I - p)g t uet (log clet g + log det Ag- A - loog det AA AA) + = ((x)))=((a++
where p - The key observation is that the algorithm only needs to know zzhh) in the HMC dynamics, and not v(h). Thus, we can replace H by any other that produces the same zh)) We show in Lemma 4 (Appendix B.2.2) that the dynamics corresponding to H above equivalent the dynamics that corresponds to much simpler Hamiltonian:
H(r, u) = ((x 2" g (I p)g to + (log det g + log det Ag -AT-
Furthermore, we have
dr
dt = 3/2 (T P) u. du
dr dr
(V(x) + Dg
Tr(g (T- P)g Dg).
dt
dt
dt
Efficient Computation Leverage Score. Even after simplifying the Hamiltonian as above, we still have term for the leverage scores, Tr(g (1) - Dg) to compute di Since (Ag- - can be extremely dense even when A is very sparse, a naive approach such as direct computation of the inverse can lead to dense-matrix multiplication. To avoid dense-matrix multiplication, our approach is based on the fact that certain entries of (Ag A A sparse Cholesky decomposition of Ag-1A [45, 5J, which can be ((n) time faster than computing obtain sparse
in du ar so that we need to compute de
the diagonal entries of P 9 AAT (Ag "AT)-'Ag
can be computed as fast as computing
(Ag (Ag- 'AT)- in many settings. We first compute the Cholesky decomposition triangular matrix L such that LL AT Then, we show that only entries of Ag- AT in sp(L) U Sp (L.)
matter in computing diag(A) (A ( (Ag A )11A), where sp(L. is the sparsity pattern of L. We give the details of our approach in Appendix B.2.3.
2.3 Diseretization
Explicit integrators such leapfrog integrator, which are commonly used for Hamiltonian Monte Carlo, are no longer symplectic on general Riemannian manifolds there have been some attempts [38] to make explicit integrators work in the Riemannian setting, its variants do not work for ill-conditioned problems.
(see Appendix C.1). Even though
Our algorithm uses the implicit midpoint method (Algorithm 3) to discretize the Hamiltonian process into steps of step size h and run the process for T iterations. This integrator reversible and symplectic (so cewasure-preserving [19], which allows to Metropolis filter to ensure the

---

distribution is correct so that we no longer need to solve ODE to accuracy to maintain the correct stationary distribution, We write H(x, v) = H (I, u) + (z, v), where
H1 0) [(x) +
2 (log det g(z) + log del Ag(x) -
1 T g (I p(x)))(x)) H2(2, 0)) 2
Starting from (ro, vo), in the first step of the integrator, we run the process on the Hamiltonian H1 with step size h with step size by solving
2 to get (m1/3, 01/3). In the second step of the integrator, we run the process on H₂
HH2 dr
21 +
2
HH2 Or
r + xg
= z +
02 =
2
2
2
iteratively using the Newton's 's method. This step involves computing the Cholesky decomposition of (Ag-1 1 A )- using the Cholesky decomposition of Ag- In the third step, we run the process h on the Hamiltonian H1 with step size 2 again to get (x1,
We state the complete algorithm (Algorithm and Algorithm 3) with details on the step size in Appendix C.I and give the theoretical guarantees in Appendix C.2 (convergence of implicit midpoint method) and Appendix D (independence of condition number).
Experiments
In this section, we demonstrate the efficiency of our sampler using experiments on real-world datasets and compare our sampler with existing samplers. We demonstrate that CRHMC is able to sample larger models than previously known be possible, and is significantly faster in terms of rate of convergence and sampling time in Section 3.2, along with convergence test in Section 3.4. We examine its behavior on benchmark instances such as simplices and Birkhoff polytopes in Section 3.3.
3.1 Experimental Setting
Settings. We performed experiments on the Standard DS12 v2 model from MS Azure cloud, which has a 2.1 Intel Xeon Platinum 8171M CPU and 28GB memory. In the experiments, we used our MATLAB and C++ implementation of CRHMC', which is available here and has been integrated into the COBRA toolbox.
We used twelve constraint-basee metabolic models from molecular systems biology in the COBRA Toolbox v3.0 [21] and ten real-world LP examples randomly chosen from NETLIB LP test sets. A polytope from each model is defined by {u E F" : Ar b, 1 < < J for A E Rmxn b b R and I, u e R,," which is input CRHMC for uniform sampling. We describe in Appendix how we preprocessed these dataset, along with full information about the datasets in Table 2.
Comparison. We used as a baseline the Coordinate Hit-and- (CHAR) implemented in two different languages. The former is Coordinate Hit-and-Rum with Rounding (CHRR) written in MATLAB [11, 20] and the latter is the same algorithm (CDHR) with an interface and a C++ library, VolEsti [6]. We refer readers to Appendix A for the details of these algorithms and our comparison setup. We note that popular sampling packages such as STAN and Pyro were not included in the experiments as they do not support constrained--ric based models. Even after transforming their formats, the transformed dataset were too ill-conditioned for those algorithms to run. CHMC in [4] works only for manifolds implicitly defined by {c(x) 0} for continuously differentiable c(x) with Dc(x) full-rank everywhere, so we could not use it for comparison.
dataset to
To evaluate the quality of sampling methods, measured quantities, the number steps per efffettive sample (i.e., mixing rate) and the sampling time per effective sample, T,. The effective sample size (ESS)* can be thought of as the number of actual independent samples,
Measurements.
e Ourr package can be run to sample from general logconeave densities and has a feature for parallelization.
WWe use the minimum of the ESS of each coordinate.
7

---

Sampling Time
Mixing Rate
103
16a
10
CHRR: dim
CRNNC: dim
Time/Sampk (s)
1oa
10a
CHRR: dim CDHR: 3.16
CRHMC: dim
Step/Sample
1o.
CDHR: 2:12
102
20
10a
160
30
10-
4
10:2
1
10
10
10
10
10a
10
10
10
10a
Dimension
Dimension
Sampling Time
Mixing Rate
CRHMC: n550
10a
CRHMC: 0.53
10
Step/S ample
Time/Sample
10a
10
10
2
10a4
10a
10:3
10a
10a
10
10a
10a
10a
NNZ
10a
Figure 3.1: Mixing rate of CRHMC and the com- Figure 3.2: Sampling time of CRHMC and the competi- petitors. Mixing rate of CRHMC was sub-linear tors. The sampling time per effective sample of CRHMC dimension and nnz preprocessed matrix was sub-quadratic
dimension and the nnz of prepro- in a model, whereas the others needed quadratically cessed matrix A in a model, while the others indicates at many steps to converge uniform distribution. In least cubic dependency on dimension. In particular for particular for dataset, CRHMC mixed up 6 dataset, CRHMC was ahle to obtain statistically inde- orders magnitude earlier than the others. Note pendent sample up to orders of magnitude faster than the mixing rate of CHAR was very close to quadratic
others. This benefit of speed-up was actually straightfor- growth when using the full-dimensional scale (the first column in Table 2).
ward from the figure, since CHRR could not obtain enough samples from instances with more than 5000 variables until it ran out of time.
taking into account correlation of samples from a target distribution. Thus the number of steps per effective sample is estimated by the total number of steps divided by the ESS, and the sampling time T, estimated as the total sampling time until termination divided by the ESS.
Each algorithm attempted to draw 1000 uniform samples, with limits on running time set to 1 day (3 days for the largest instance ken ken_18) and memory usage to 6GB. If an algorithm passes either the time or the memory limit, we stop the algorithm and measure the quantities of interest based on samples drawn until that moment. After getting uniform samples, we thinned the samples twice to ensure independence of samples; first we computed the ESS of the samples, only kept ESS many samples, and repeated this again. We estimated the above quantities only if the ESS more than 10 and algorithm does not run into any eπror while running..
3.2 Mixing Rate and Sampling Time
Sub-linear Mixing Rate. We examined how the number of steps per effective sample grows with the number of nonzeros (nnz) of matrix A (after preprocessing) and the number of variables (dimension in the plots). To this end, we counted the total number of steps taken until termination of algorithms and divided it by the effective sample size of drawn samples. Note that we thinned twice to ensure independence of samples used.
The mixing rate of CRHMC was sub-linear in both dimension and nnz, whereas previous implemen- tations based on CHAR required at least n² steps per sample as seen in Figure 3. On the dataset,
'When running CDHR from the VolEsti package some instances, we got error message "R session aborted and R encountered a fatal error".

---

Bio Model
Vars (o)
maZ
CRHMC
CHRR
CDHR
ecoli
95
291
0,0098
0.0365
0.0022
LP Model
Vars (e)
noz
CRHMC
CHRR
CDHR
cardiac_ mit
220
228
0.0100
0.0059
0.0005
israel
316
2519
0.1186
1.2224
0.4426
Aci D21
851
1758
0.4257
0.6884
0.2974
giid_poc
1160
2393
0.2199
40.988
18.468
Aci MRL5
994
2859
0.9624
2,0668
0.5237
251v47
1876
10566
0.8159
199.9
49176
1069
2951
0.9608
1.9395
0.9622
pilot in
2267
11886
1.3490
5059*
Ac7_20731
1090
2946
0.1540
2.3014
1,1086
sctop2
2500
7334
0.6752
520.2
Aci PHEA
1561
4640
0.3701
12.06
ship081
4363
9434
0.6258
6512
IAF1260
2382
6368
4,4355
3687.2
a
7248
17368
2.2205
30455*
301366
2583
7284
4,1608
70.5
35.556
woodw
8418
231.88
2.0689
30307
Reconl
3742
8717
0.7184
208.5
80hau
1206
22341
11.881
47432*
Recon2
7440
19791
2,61166
10445*
kon_18
154699
295946
1616.3
Recon3
13543
48187
31.144
29211
Table 1: Sampling time effective sample of CHRR and CRIMC. We note that CRHMC is 1000 times faster than CHRR on the latest metabolic network (Recon3). Sampling time with asterisk (*) indicates that the effective sample size less than 10.
mixing rate attained was up to orders of magnitude faster for CRHMC compared to CHAR, imply- ing that CRHMC converged to uniform distribution substantially faster than the other competitors. This gap in mixing rate increased super-lnearlly in dimension, enabling CRHMC to run on large instances of dimension up to 100000.
Sub-quadratic Sampling Time. We next examined the sampling time T, in terms of both the unz of A and the dimension of the instance. We computed the runtime of algorithms until their termination divided by the effective sample size of drawn samples, where we ignored the time it takes for preprocessing. Note that the sampling time T, is essentially multiplication of the mixing rate and the per-step complexity (i.e., how much time each step takes).
As shown in Figure 3.2 and Table I, we found that the per-step complexity of CRHMC was small enough to make the sampling time sub-quadratic in both dimension and nnz, whereas CHAR had at least a cubic dependency sampling time of CRHMC was up to orders of magnitude less than that of CHRR and CDHR. While CHRR can be used on dimension only up to a few thousands, increasing benefits of sampling time in higher dimension allows CRHMC to run on dimension up to 0.1 million.
dimension, despite of a low per-step complexity. On our dataset, the
3.3 CRHMC on Structured Instances
To see the behavior of CRHMC on very large instances, we ran the algorithm on three families of structured polytopes hypercube, simplex, and Birkhoff polytope - up dimension half-million. We attempted to draw 500 uniform samples with I day time limit (except for 2 days for half-million- dimensional Birkhoff polytope). The definitions these polytopes are shown Appendix A.I.
Mixing Rate
Sampling Time
10
10a
Cube: dim 0.34
Cube: dim 0.76
dim 0.33
Time/Sample (S)
Simplex: dim 0.83
Birthoff: dim 0.43
Birkhoff: dim 1.08
Step/Sample
20a
10
10
100
20
10
20
10
10:3
10
102
Dimension
10
20
10
10a
Dimension
30
10a
6
Figure 3.3: Mixing rate and sampling time on structured polytopes including hybercubes. simplices. and Birkhoff polytopes. CRHMC is scalable up to million dimension on hypercubes and simplices and up to million dimension on Birkhoff polytopes. We note that on the 0.5 million dimensional Birkhoff polytope the ESS is only 16, which is not reliable compared to the ESS on the other instances.
To the best of our knowledge, this the first demonstration that it is possible to sample such a large model. As seen in Figure 3.3, CRHMC can scale smoothly up to half-million dimension on hypercubes and simplices and up to dimension 10 for Birkhoff polytopes (we could not obtain reliable estimate of mixing rate and sampling time on the half-million dimensional Birkhoff

---

polytope, as the ESS is only 16 after 2 days). However, we believe that one can find room for further improvement of CRHMC by tuning parameters or leveraging engineering techniques. We also expect that CRHMC enables us to estimate the volume of B, for 2 20, going well beyond the previously best possible dimension.
3.4 Uniformity Test
We used the following uniformity test to check whether samples from CRHMC form the uniform distribution over a polytope P: check that the fraction of the samples in the scaled set x P is proportional to z As in Figure 3.4, the empirical CDFs of the radial distribution to the power of (1/dim) are close to the CDFs of the uniform distribution over those polytopes.
Empirical CDF of the Radial Distribution
Empirical CDF of the Radial Distribution
Empitial CDF Uniform
Emprical GDF
Unilom
0.8
0.8
0.6
0.6
F06
0.4
0.4
0.2
0.2
0.2 0.4
0.6 0.0
0 0.2 0.4
0.5
0.8
x
Figure 3.4: We plot the empirical cumulative distribution function of the radial distribution to the power of (1/dim) with 1000 ESS obtained by running CRHMC A7CC-49176 (952 x 1069, left) and Aci-PHEA (1319 x 1561, right), and in the plot z-axis is the scaling factor. We can observe the CDFs are very close to the CDFs of the uniform distribution over the polytopes defined by two instances.
The authors grateful to Ben Cousins for helpful discussions, and to Ronan Fleming, Ines Thiele and their research groups for advice on metabolic models. This work was supported in part by NSF awards DMS-1839116, DMS-1839323, CCF-1909756, CCF-2007443 and CCF-2134105.
Acknowledgement.
References
[1] 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] Ivona Bezáková, Daniel Štefankovië, Vijay V Vazirani, and Eric Vigoda. Accelerating simulated annealing for the permanent and combinatorial counting problems. SIAM Journal on Computing (SICOMP), 37(5): 1425:)::422-11 1454, 2008
[3] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul A. Szerlip. Paul Horsfall, and Noah D. Goodman. Pyro: Deep Universal Probabilistic Programming. Journal of Machine Learning Research (JMLR), 20:28:1- 28:6, 2019.
[4] Marcus Brubaker, Mathieu Salzmann, and Raquel Urtasun. A family of MCMC methods on implicitly defined manifolds. In Artificial intelligence and statistics (AISTATS), pages 161-172, 2012.
[5] Yogin E Campbell and Timothy A Davis. Computing the sparse inverse subset: inverse multifrontal approach. University of Florida, Technical Report TR-95-021, 1995.
[6] Apostolos Chalkis and Vissarion Fisikopoulos. volEsti: Volume approximation and sampling for convex polytopes in R. arXiv preprint arXiv:2007.01578 2020.
[7] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. Journal of Muchine Learning Research (JMLR), 21:92-1,
10

---

[8] Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan. Underdamped Langevin MCMC: non-asympiotic analysis. Conference on Learning Theory (COLT), pages 300-323. PMLR, 2018.
[9] Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, Philippe Rigollet, and Austin Stromme. Exponential ergodicity of marror-Langevin diffusions. Advances in Neural Information Process- ing Systems (NeurIPS), 13:19573-19585, 2020.
[10] Adam D Cobb, Atilm Güney Baydin, Andrew Markham, and Stephen J Roberts. Introducing an explicit symplectic integration scheme for Riemannian manifold Hamiltonian Monte Carlo. arXiv preprini arXiv: 1910.06444 2019.
[11] Ben Cousins and Santosh Vempala. A practical volume algorithm. Mathematical Programming Computation, 8(2):133-160, 2016.
[12] Timothy Davis. Direct methods for spurse linear systems. SIAM, 2006.
[13] James Demmel. Applied numerical linear algebra. SIAM, 1997.
[14] Simon Duane, Anthony Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 195(2):216-22) 1987.
[15] Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: Metropolis-Hastings algorithms are fast! In Conference on Learning Theory (COLT), pages 793-797. PMLR, 2018.
[16] Martin Dyer, Alan Frieze, and Ravi Kannan. A random polynomial-time algorithm for approxi- mating the volume of convex bodies. Journal of the ACM (JACM). 38(1):1-17, 1991.
[17] 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.
[18] Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. SIAM, 2008
[19] Ernst Hairer, Marlis Hochbruck, Arieh Iserles, and Christian Lubich. Geometric numerical integration. Oberwolfach Reports, 3(1):805-88) -82, 2006.
[20] Hulda S Haraldsdóttir, Ben Cousins, Ines Thiele, Ronan MT Fleming, and Santosh Vempala. Chrr: coordinate hit-and-run with rounding for uniform sampling of constraint-based models. Bioinformatics, 33(11):774111743 2017.
[21] Laurent Heirendt, Sylvain Arreckx, Thomas Pfau, Sehastián N Mendoza, Anne Richelle, Almut Heinken, Hulda S Haraldsdóttir, Jacek Wachowiak, Sarah M Keating, Vanja Vlasov, et al. Creation and analysis of biochemical constraint-baser models using the COBRA Toolbox 3.0. Nature protocols, 14(3):639-702, 2019.
[22] Matthew D Hoffman, Andrew Gelman, et al. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research (JMIR), 15(1):1593-1623 2014.
[23] Mark Jerrum, Alistair Sinclair, and Eric Vigoda. A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries. Journal of the ACM (14C)),,5(4))661-66))) 2004.
[24] He Jia, Aditi Laddha, Yin Tat Lee, and Santosh Vempala. Reducing isotropy and volume to KLS: an O* (n²abb)) volume algorithm. Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing (STOC), pages 961-974, 2021.
[25] Ravi Kannan, László Lovász, and Miklós Simonovits. Random walks and an O*(r²) volume algorithm for convex bodies. Random Structures & Algorithms, 11(1):1-50, 1997.
[26] Ravindran Kannan and Hariharan Narayanan. Random walks on polytopes and an affine interior point method for linear programming. Mathematics of Operations Research, 37(1):1-20, 2012.
11

---

[27] Zachary A King, Justin Lu, Andreas Dräger, Philip Miller, Stephen Federowicz, Joshua A Lerman, Ali Ebrahim, Bernhard Palsson, and Nathan E Lewis. BiGG Models: A platform for integrating, standardizing and sharing genome-scale models. Nucleic acids research, 44(D1))D51--D52, 2016.
[28] Yunbum Kook, Yin Tat Lee. Ruoqi Shen, and Santosh S. Vempala. Condition-number- independent Convergence Rate of Riemannian Hamiltonian Monte Carlo with Numerical Integrators. arXiv preprint arXiv:2110002219 2022.
[29] Aditi Laddha and Santosh Vempala. Convergence of Gibbs sampling: Coordinate Hit-and-Run mixes fast. The 37th International Symposium Computational Geometry (SoCG), 2021.
[30] Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Logsmooth gradient concentration and tighter runtimes for metropolized Hamiltonian Monte Carlo. In Conference on Learning Theory PMLR, 2020.
(COLT), pages
[31] Yin Tat Lee and Santosh Vempala. Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation. In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing (STOC), pages 1 1115-1121, 2018.
[32] Yin Tat Lee and Man-Chung Yue. Universal barrier is ceseof-concordant. Mathematics of
Operations Research, 2021.
[33] Nathan E Lewis, Harish Nagarajan, and Bernhard O Palsson. Constraining the metabolic genotype-phenootpe relationship using phylogeny of in silico methods. Nature Reviews Microbiology, 10(4):291-305 2012.
[34] László Lovász and Santosh Vempala. Hit-and-run from a corner SIAM Journal on Computing (SICOMP), 35(4);985-1005, 1005, 2006.
[35] Hariharan Narayanan and Piyush Srivastava. On the mixing time of coordinate Hit-and-Rum Combinatorics, Probability and Computing, pages 1-13, 2021.
[36] Radford M Neal et al. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11):2, 2011.
[37] Yurii Nesterov and Arkadii Nemirovskii. Interior-point polynomial algorithms in convex programming. SIAM, 1994.
[38] Pauli Pihajoki. Explicit methods in extended phase space for inseparable hamiltonian problems. Celestial Mechanics and Dynamical Astronomy, 121(3):211 11-2231, 2015.
[39] Sebastian Reich. Symplectic integration of constrained Hamiltonian systems by Runge-Kuuu -Kutta methods. University of British Columbia, Department of Computer Science, 1993.
[40] Gareth O Roberts and Richard L Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pages 341-363, 1996.
[41] John Salvatier. Thomas V Wiecki, and Christopher Fonnesbeck. Probabilistic programming in Python using PyMC3. PeerJ Computer Science, 2xe55, 2016.
[42] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances in Neural Information Processing Systems (NeurIPS). 32, 2019.
[43] Umut Simsekli, Roland Badeau, Taylan Cemgil, and Gaël Richard. Stochastic quasi-newton Langevin Monte Carlo. In International Conference on Machine Learning (ICML), pages 642-651. PMLR, 2016.
[44] Stan Development Team. RStan: the R interface to Stan, 2020. R package version 2.21.2.
[45] Kazuhiro Takahashi. Formation of sparse bus impedance matrix and its application to short circuit study. In Proceeding of PICA Conference, June, 1973, 1973.
[46] Ines Thiele, Neil Swainston, Ronan MT Fleming, Andreas Hoppe, Swagatika Sahoo, Maike Aurich, Hulda Haraldsdottir, Monica L. Mo, Ottar Rolfsson, Miranda D Stobbe, et al. community-driven global reconstruction of human metabolism. Nature biotechnology. 3 31(5):419-425, 2013.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you descrihe the limitations of your work? [Yes] Section 3.3, we mention that
we could not obtain a reliable estimates of measures in dimension higher than 10 and we believe there is room for further improvement of our algorithm.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] Our work has no negative societal impacts.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] In Appendix (to which we defer our proofs), we clearly state any assumption if needed.
(b) Did you include complete proofs of all theoretical results? [Yes] We include all proofs in Appendix.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes] See Section and Appendix A.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [N/A] Our algorithm does not require training.
(c) Did you report error bars (c.g., with respect to the random seed after running experi- ments multiple times)? [N/A]
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] We specified our experimental setting in Section 3.
4. If you using existing assets (e.g., code, data, models) curating/releasing assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] We made clear reference existing packages in Section 3.
(b) Did you mention the license of the assets? [Yes] See Section 3.
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes] We provide anonymized link
algorithm package in Section 3.1.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] All existing packages we used in the experiment are publicly available for academic purpose.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A] Our datasets consist of metabolic models from molecular system biology and L. examples, which have no personally identifiable information or offensive content.
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "kujBXZvWJ_", "reviewer": "Reviewer_Cyv8", "summary": "The authors propose a method to sample from a distribution using Riemannian Hamiltonian Monte Carlo while respecting a set of constraints. The efficiency of the method relies on a combination of tricks from linear algebra, while the proposed Riemannian metric is based on the Jacobian of the constraints and a barrier function. In the experiments it is demonstrated the behavior of the sampler in different scenarios.", "questions": "In particular, it is not clear to me: \n1. why the metric should not be invertible\n2. how the associated null spaces are related to the problem\n3. why the metric is then constructed as in Eq. 2.6 \n4. how the function $\\phi$ is defined in practice.\n\nIn general, I find hard to access the Section 2.1. I think some discussion about the geometrical intuition/properties of these steps together with some figures could make the paper much more accessible. I believe that this part can be improved.", "limitations": "The authors discuss some limitations of the proposed method.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- Originality: As the authors mention their approach is based on a previous work, which they extend such that to make it much more efficient computationaly. I think that the related work is cited and discussed very briefly in the intro. I am not an expert in the topic of sampling, but the proposed approach seems like a sufficient contribution.\n\n- Quality: The technical part of the paper seems ok, but I have not checked it carefully. The claims are supported by theoretical analysis and experimentally validated.\n\n- Clarity: I think that the paper is ok written, but the Section 2 seems a bit too technical without the necessary explanations to help the reader. For example, a lot of linear algebra results are used, but I have the feeling that the gist and intuition are missing. I believe that this makes the paper hard to access from non-expert readers like me. (Please see questions)\n\n- Significance: From the claims and the empirical results it seems that the method indeed solves the problems better (more efficient) than previous approaches. It is quite likely that this work would be useful when sampling under constraints.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "geyjbu0gLzP", "reviewer": "Reviewer_y2k1", "summary": "This paper gives a method for sampling log-concave distributions supported on compact convex sets, based on Riemannian Hamiltonian Monte Carlo (RHMC). Equality constraints are enforced via a Lagrange multiplier and inequality constraints are enforced via a carefully chosen Riemannian metric constructed from a self-concordant barrier function. Algorithms for exploiting sparsity for matrix inversion enable sampling from high-dimensional distributions. Tests show substantial improvement in sampling times in comparison to existing constrained samplers. Formal guarantees are provided for convergence of the associated integration scheme the convergence to the stationary distribution is bounded.", "questions": "Would it be possible to extend this to noncompact sets? How about non-convex sets and non-log-concave distributions? What are the challenges?", "limitations": "The authors address many technical limitations, but as described above, some fairly serious technical limitations are not adequately discussed without digging into the technicalities. \n\nThe authors indicate not potential negative societal impact, and I agree.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\nThis is a nice paper. It gives strong practical results for an important, but largely overlooked problem and gives a thorough theoretical analysis. The paper is very well-organized, with helpful guides to the structure of the appendix. \n\nWeaknesses:\nProbably the biggest weakness is that some of the assumptions are not clear from the problem setup. In particular, it is just stated that the constraint set, $\\mathcal{K}$, must be convex. However, it really must be compact. Many of the arguments implicitly depend on this. These restrictions should definitely be clarified in the problem setup, and ideally also in the abstract.\n\nIf $\\mathcal{K}$ is not compact, the cross-ratio distance from Definition 5 is not well-defined. It might be possible to extend this definition to the case in which $\\mathcal{K}$ does not contain a line by using appropriate limits. However, in this case the constant $G$ from Lemma 17 could be zero. (Consider the case of the non-negative orthant in $\\mathbb{R}^2$.)\n\nIf $\\mathcal{K}$ is not compact and the $f$ is linear, then the corresponding density may not exist, rendering arguments about stationarity from B.3 irrelevant.\n\nThere are also some small typos and unclear arguments. Unfortunately, I did not have the time to step through every  argument in detail. But, I did catch a few issues:\n\n- In line 575, it should read $MNM=M$\n- In many places, \"dynamic\" is used where \"dynamics\" would be more common. For example, in 603 \"the HMC dynamic\" would normally be \"the HMC dynamics\"\n- Similarly, \"the dynamic on x\" would normally be \"the dynamics of x\", in Lemma 4.\n- $P$ gets redefined in line 601, and clashes with the definition from Lemma 2 in a confusing way.\n- In line 798, we should have $\\|x\\|_A = \\sqrt{x^\\top A x}$\n- Definition 5 implies that $\\mathcal{K}$ must be compact. This should be stated more explicitly in the problem setup.   \n- I think there is a typo in the first inequality after line 967. I feel like this is based on a lower bound on the right side of the equation after line 964. So, the second side should have an integrand of $\\mathcal{T}_x(S^c)$.\n- I think the factor on the second inequality below line 967 should be $\\frac{1}{40}$ and not $2/5$. This is because the bound seems to follow from lower bounding the integrands by $0.05=\\frac{1}{20}$ for both terms. I don't think this really matters for the rest of the proof. \n- The last inequality below line 968 should be strict.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Kh4C7YQg-AI", "reviewer": "Reviewer_CjP3", "summary": "The paper develops a constrained version of Riemannian Hamiltonian Monte Carlo, maintaining both sparsity and constraints. The paper shows that ill-conditioned, non-smooth, constrained distributions in very high dimensions, upwards of 100,000, can be sampled efficiently. The work is supported by several theoretical results on the correctness and computational efficiency, as well as bounds on the mixing time of the CRHMC and its discretizations.", "questions": "(see above)", "limitations": "None noted.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n\n- Significance:\n    - The paper considers an important problem in machine learning that is of great interest to the broader machine learning community.\n    - The proposed framework can address sampling problems that are deemed infeasible previously. The paper shows that ill-conditioned, non-smooth, constrained distributions in very high dimensions, upwards of 100,000, can be sampled efficiently.\n    - The experiments also demonstrate the approach's superiority in comparison to existing methods in terms of rate of convergence and total sampling time.\n- Quality and Clarity:\n    - The paper is well-written and well-organized.\n    - The theoretical analyses of the work are rigorous and of very high quality. The experiments are extensive and support the findings of the theoretical part.\n\n### Weaknesses\n\n- Quality:\n    - I believe the proof of Theorem 8 (Appendix, on the correctness of the discretized version of CRHMC) is incorrect. In principle, the detailed-balance conditions need to be verified on two arbitrary measurable sets. This could be done by verifying the condition pointwise (as done in Theorem 8) if the state space is discrete or the proposal distribution is absolutely continuous. HMC proposals, however, contain point masses, and the detailed-balance condition can be verified pointwise only if the proposal is volume-preserving. This can be seen in the proof of Theorem 6 for the ideal version of RCHMC but is absent in Theorem 8. \n\n      It is worth noting that this a not a minor point since the rest of the theoretical analyses (e.g., bounds on mixing time) relies on the convergence of discretized CRHMC obtained through this Theorem.\n    - This brings a broader point about the choice of integrator for RCHMC, namely, why the implicit midpoint is preferred over the leapfrog integrator (which is known to be volume-preserving). The implicit midpoint is described in the paper with no context or references and very few details about its properties are spelled out. I think such an important component of the algorithm deserves some details treatment on its background and properties.\n    - A minor point: the paper seems to be reluctant to state that RCHMC is ergodic or converges into the formal statement of the theorems: Theorem 6 states that RCHMC satisfies detailed balance; Theorem 8 states that RCHMC has a stationary distribution. Mathematically, those two statements are much weaker than ergodicity and are not strong enough to support later discussions such as Line 816.\n- Presentation:\n    - To continue on the point above, while the development of CRHMC follows well-founded rationales, I find it very hard to separate the contributions of the current construction from pre-existing methods in the field, simply because the paper doesn’t provide enough contexts. The construction of the constraint term in the Hamiltonian $\\lambda(x, v)^T c(x)$ is rigorous, but is it an original concept, or does it draw inspiration from pre-existing works? The barrier function is a well-known tool in constrained HMC, but was introduced with no proper references. Similarly, several aspects of efficient computation and simplifications in Section 3 seems to base on are pre-existing ideas. I believe a more thorough review of the key ideas would strengthen the manuscript greatly.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "K2NKERBJKX", "reviewer": "Reviewer_kCFo", "summary": "This paper provides a variation of the Constraint Riemannian Hamiltonian Monte Carlo (RHMC). Constraint variations of HMC deal with  sampling from a target distribution that by means of some constraints, is (implicitly) defined on a manifold. In the present work it is assumed that the constraint has the form c(x) = Ax − b. Such assumptions has led to fast computations and as it is mentioned, the proposed algorithm can be run on target spaces with dimensionality as high as 100K. ", "questions": "* The proposed work is only compared against (2 implementations of) CHRR algorithm. I suggest that the comparison be with other variations of constraint HMC (e.g. Brubaker's constraint baseline HMC, constraint Metropolis-Hastings) and/or more recent implementations of CHMC algorithms such as Graham's \"Manifold Markov chain Monte Carlo methods for Bayesian inference in diffusion models\".\n\n* I think CHRR is for sampling from (constraint) uniform distributions. Are the presented experimental target densities uniform? If no, how CHRR is used and if yes, what is the usage of Riemannian HMC (since gradients are always 0) ?\n\n", "limitations": "\nNot applicable. This work which is theoretical. ", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "* It seems to me that the contribution of the presented work is incremental. The idea of constraint HMC is already proposed in a more general setting by: \n    Brubaker et al, \"A Family of MCMC Methods on Implicitly Defined Manifolds\" \nwhere Constrained Riemann Manifold HMC (CRHMC) is mentioned as a specific case. The present work correctly gives this credit to Brubaker et al. (in Line 135) but then claims Brubaker's CRHMC \"requires eigenvalue decomposition and is not efficient for large problems\". As far as I understand, Brubaker's work does not necessarily require an eigenvalue decomposition and the present work is indeed an instance of Brubaker's work where extra assumptions/restrictions mentioned in Lines 192 to 220, has led to efficient computations. \n\n* If I have understood correctly and the contribution is efficient computations under extra assumptions/restrictions, then the robustness of this work could possibly be demonstrated by extensive empirical experiments. Nonetheless, the provided experimental results seem minimal to me (the proposed algorithm is only compared against one benchmark) which is suitable for sampling from (constraint) uniform distributions only. I will return to this in the \"Questions\" section. \n\n* In general the paper's writing could be improved and notations/terminology could be defined more rigorously.   \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "8oklbBzPW1q", "reviewer": "Reviewer_jRWB", "summary": "The manuscript proposes a MCMC sampler based on the Constrained HMC method for sampling from a log log-concave density $\\pi(x) \\propto \\exp[-f(x)]$ constrained on a convex subset $x \\in K \\subset \\mathbb{R}^d$ and a constraint $c(x)=0$. The main innovations are computational tricks to scale the method to high-dimensional settings and exploit the sparsity pattern of constraints of the type $Ax=b$ for a sparse matrix $A$.", "questions": "The manuscript proposes a MCMC sampler based on the Constrained-HMC method for sampling from a log-concave density $\\pi(x) \\propto \\exp[-f(x)]$ and constrained to stay within a convex subset $x \\in K \\subset \\mathbb{R}^d$ and and an implicit manifold defined through the constraint $c(x)=0$. The main innovations are computational tricks to scale the method to high-dimensional settings and exploit the sparsity pattern of constraints of the type $Ax=b$ for a sparse matrix $A$.\n\nThe proposed method seems to deliver state-of-the-art results when evaluated in a range of challenging settings (the numerical studies are well executed). Furthermore, the authors have incorporated their method within a popular Bioinformatics library.\n\n\n1. I am not finding the discussion of sparsity well-motivated or explained. The authors mention that sparsity is preserved, but it takes quite a bit of time for the reader to realize that (I think) it is related to the sparsity pattern of A within the constraint Ax=b.\n\n2. The authors are not clear within the main paper (although it is somehow clarified in the Supp Material) whether they are considering MC methods with or without an accept reject. For example, several of the cited papers on (variation of) the Langevin dynamics do not have an accept-reject step (and are consequently only approximately sampling from the posterior).\n\n3. It would be best to make it clearer whether the paper concentrates on constraints of the type Ax=b, or more general implicit manifolds. It seems like all the experiments are with linear constraints.\n\n4. It was not clear to me whether the idea of using barrier functions within an HMC scheme is new or not. If it is not, it would be best to clarify. Something that is not very well explained is the following. Suppose that one runs a standard HMC algorithm (with explicit leapfrog integrator) in $\\mathbb{R}^d$ constrained to a convex set $K \\subset \\mathbb{R}^d$ (i.e. without an equality constraint $c(x)=0$). For a barrier function $\\phi$ and positive definite mass matrix $M(x) = \\nabla^2 \\phi(x)$, there does not seem anything that prevents the Markov chain to step out of the convex set $K$ because of discretization error. Is it different when using the implicit midpoint method for integrating the Hamiltonian dynamics? Can the authors clarify this? \n\n5. Do the authors use the assumption that the target density is log-concave? where?\n\n6. For general constraint $c(x)=0$, the CHMC requires a \"projection\" step to bring back the Markov chain back to the manifold $c(x)=0$. In practice (i.e. without further assumptions on the function $c$), this projection step may fail. It is the reason why the \"reversibility checks\"  have been developed (eg. see [1]).\n\n7. There is a discussion that describes how having a position-dependent mass matrix can help take the geometry of the distribution into account for improved mixing properties (i.e. Riemaniann-HMC/etc..). But, if I am not mistaken, the authors are simply using as mass matrix the projection of $\\nabla^2 \\phi(x)$ on the tangent space to the manifold $c(x)=0$, which is *not* taking the geometry of the target distribution (eg. anisotropy of scales) into account (i.e. the target distribution $\\pi(dx)$ is not used to design the mass matrix)\n\n8. No discussion of how to tune the leapfrog step size and momentum refresh rate (i.e. coefficient $\\beta$)?\n\n9. The authors reports the ESS, but do not specify the ESS of what functional of the trajectory is considered\n\n\n\n**Minor**:\n1. if space is an issue, I do not think that it is necessary to spend so much time explaining to NeurIPS readers that sampling is important and difficult. The introduction can be shortened, if needed.\n\n\n**[Answers by the authors acknowledged 4/8/2022]**\nThank you for the careful answers that have clarified my misunderstandings. It is indeed a very solid piece of work!\n\n**References:**\n1. \"Monte Carlo on manifolds: sampling densities and integrating functions\" by Emilio Zappa, Miranda Holmes-Cerfon, Jonathan Goodman", "limitations": "OK.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The proposed method seems to deliver state-of-the-art results when evaluated on a range of challenging settings (the numerical studies are well executed). Furthermore, the authors have incorporated their method within a popular Bioinformatics library.\n\nI have a number of questions, as well as remarks on the presentation of the paper, that are listed in the next section.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept