Unbiased constrained sampling with Self-Concordant Barrier Hamiltonian Monte Carlo
Authors: Maxence Noble
Valentin De Bortoli
Alain Durmus
Venue: NeurIPS 2023 poster (2023)
Keywords: Hamiltonian Monte Carlo
Riemannian manifold
self-concordant barrier
constrained sampling

Abstract: In this paper, we propose Barrier Hamiltonian Monte Carlo (BHMC), a version of the
  HMC algorithm which aims at sampling from a Gibbs distribution $\pi$ on a manifold
  $\mathsf{M}$, endowed with a Hessian metric $\mathfrak{g}$ derived from a self-concordant
  barrier. Our method relies on Hamiltonian
  dynamics which comprises $\mathfrak{g}$. Therefore, it incorporates the constraints defining
  $\mathsf{M}$ and is able to exploit its underlying geometry. However, 
  the corresponding Hamiltonian dynamics is defined via non separable Ordinary Differential Equations (ODEs) in contrast to the Euclidean case. It implies unavoidable bias in existing generalization of HMC to Riemannian manifolds. In this paper, we propose a new filter step, called ``involution checking step'', to address this problem. This step is implemented in two versions of BHMC, coined continuous BHMC (c-bHMC) and  numerical BHMC (n-BHMC) respectively.
  Our main results establish that these two new algorithms  generate reversible Markov
  chains with respect to $\pi$ and do not suffer from any bias in comparison to previous implementations. Our conclusions are supported by numerical experiments where
  we consider target distributions defined on polytopes.

Unbiased constrained sampling with Self-Concordant Barrier Hamiltonian Monte Carlo
Maxence Noble*
CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, 91 120 Palaiseau, France
Valentin De Bortoli
Alain Oliviero Durmus
Computer Science Department, ENS, CNRS, PSI. University
CMAP, CNRS, École polytechnique. Institut Polytechnique de Paris, 91 120 Palaiseau, France
Abstract
In this paper, we propose Barrier Hamiltonian Monte Carlo (BHMC), version of the HMC algorithm which aims at sampling from a Gibbs distribution on manifold M, endowed with Hessian metric g derived from a self-concordant barrier. Our method relies on Hamiltonian dynamics which comprises g. Therefore, it incorporates the constraints defining M and able to exploit its underlying geometry. However, the corresponding Hamiltonian dynamics is defined via non separable Ordinary Differential Equations (ODEs) in contrast to the Euclidean case. It implies unavoidable bias in existing generalization of HMC to Riemannian mani- folds. this paper, we propose a new filter step, called "involution checking step", to address this problem. This step is implemented in two versions of BHMC, coined continuous BHMC (c- (c-BHMC) and numerical BHMC (n-BHMC) respectively. Our main results establish that these chains with respect to π and do not suffer from any bias in comparison to previous implementations. Our conclusions supported by numerical experiments where we consider target distributions defined on polytopes.
algorithms generate reversible Markov
1
Introduction
Markov Chain Monte Carlo (MCMC) methods is one of the primary algorithmic approaches obtain approximate samples from a target distribution л. They have been successively applied over these past decades in large panel of practical settings, (Liu & Liu, 2001). In particular, gradient-base based MCMC methods have shown their efficiency and robustness in high-dimensional settings, and come nowadays with strong theoretical guarantees, (Dalalyan, 2017: Durmus & Moulines, 2017). However, they still struggle in facing the case where the target distribution supported on a constrained subset M of Rd (Gelfand et al., 1992; Pakman & Paninski, 2014; Lan & Shahbaba, 2015). Yet, this problem appears in various fields; see e.g., (Morris, 2002: Lewis et al., 2012; Thiele et al., 2013) with some important applications in computational statistics and biology.
Drawing samples from such distributions indeed a challenging problem that has been intensively studied in the literature, (Dyer & Frieze, 1991; Lovász & Simonovits. 1993; Lovász & Kannan, 1999; Lovász & Vempala, 2007; Brubaker et al., 2012; Cousins & Vempala, 2014; Pakman Paninski, 2014; Lan & Shahbaba, 2015; Bubeck al., 2015). In particular, some recent extensions of the
"Coresponding author. Contact
waw.eeeennbbeebbulllloo @ pol.eecniiuuu..dduu
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

popular Cetropolis-Hastting (MH) algorithm to constrained spaces consist in designing proposals based on dynamies for which invariant, (Zappa al., 2018: Lelièvre et al.. 2019, 2022). In practice, the corresponding solutions are numerically integrated based on implicit and symplectic schemes, which however come with additional difficulties. As first observed by Zappa al. (2018), an additional "involution checking step" in the usual MH filter is necessary to ensure that the resulting Markov kernel admits x as a stationary distribution.
In this paper, we aim at generalizing this family of methods by taking into account the geometry of the constrained subspace, building on the Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) algorithm introduced by Girolami & Calderhead (2011) I). Similarly to HMC, (Duane et al., 1987; Neal et al., 2011: Betancourt, 2017). RMHMC aims to target positive distribution π on E' and relies on the integration of a canonical Hamiltonian equation. In contrast, RMHMC incorporates some geometrical information in the definition of the Hamiltonian via Riemannian metric g on IR.. When considering applications of RMHMC to a convex, open and bounded subset M. a natural choice for this metric the Hessian metric associated with a self-concordant barrier on M. (Nesterov & Nemirovskii, 1994), as suggested by Kook et al. (2022a). We adopt here the same approach and now focus on a constrained subspace M which is assumed to be equipped with an appropriately designed "self-concordant" metric, (Nesterov & Nemirovskii, 1994).
Given this setting, we propose the BHMC (Barrier HMC) algorithm. To the best of our knowledge, this is the first version of RMHMC incorporating an "involution checking step" to assess the issue of asymptotic bias arising from the use of implicit integration methods solve the Hamiltonian dynamics. In this work, we focus on the numerical version of BHMC (n-BHMC), for which we provide numerical implementation. We refer to Appendix D for an introduction to the ideal version of BHMC, called continuous BHMC (c-BHMC), where it is assumed that we have access to the continuous Hamiltonian dynamics but which cannot be computed in practice. For each algorithm, we rigorously prove that the associated Markov chain preserves x.
The rest of the paper is organized as follows. In Section 2, we introduce our sampling framework and review some background on self- concordnceee and Riemannian geometry. In Section 3, we present n-BHMC, and derive theoretical results for this algorithm in Section 4. We review related works in Section 5 and provide numerical experiments for n- R-BHMC in Section 6.
Notation. For any f C3 , R), we denote by Df (and D2 the Jacobian (resp. the Hessian) of f. For any A E Raxd and any (x,y) E Rd x Rd, we denote by A D f() the vector (TT(AT(D)) (D* (())))))) E Rd, and define [x) v) as the vector x. D3 J f(x) E Ted. For any positive-defintte matrix A Baxd () stands for the scalar product induced by A on Rd. defined by (x, v) = (x,, Ay). The "momentum reversal" operator 8 Rd x R2 Rd x Rd is defined for any (x,p) E Fd x Rd by s(x,p) = (z, -p). Let E, be two sets and fa E = 2F, where 2 is the set of sets F. We say that h set-valued map. Note that any map g : + can be extended to set-valued map by identifying, for any x E E, g(z) and {g(z)}. Let f F + 2 We define where by convention Ue- = 0. For any topological space X, we denote B(X) its Borel sets. For any probability measure μ E P (X) and measurable map &p X Y, we denote # (Y) the pushforward of μ by yp In general, will equivalently denote by 2 or (f. p) a state of the Hamiltonian system.
foh: E E 2c for any x E E by f h(x) /
2 Setting and Background
In this paper, consider open subset M c Rd. and we aim at sampling from target distribution ı given for any x E M by
dx = exp]==
where V E CMM R) and Z JM e We view M embedded d-dimensional submanifold of Red equipped with a metric g. satisfying the following assumptions.
AI. M is an open convex bounded subset of Idd
A2. There exists o, a o-regular and u-self-concorlont barrier on M such that g D2
We provide in Section 2.1 basic Riemannian facts along with the definition of the Hamiltonian dynamics of RMHMC and introduce self-concordance and c-regularity in Scction 2.2.

---

2.1 Riemannian Manifold Hamiltonian dynamies
Basics on Riemannian geometry. Let M be d' d-dimensional smooth manifold, endowed with metric g. Denoting by dæ a dual coframe, (Lee, 2006, Lemma 3.2.), we recall that the Riemannian volume element corresponding to (M, g) is given in local coordinates by
dvolm(2) = = det(g))xx.
We denote by T,M the dual the langent space T is vector space naturally endowed with the scalar product (- a(xx) 1,, (Mok, 1977). We recall that the cotangent bundle T*M is defined by T*M = [a) U T.M. This set is a 2d-dimensional manifold endowed with a Riemannian metric g" given by Mok (1977), inherited from g. With this metric, the volume form on T does not depend on g anymore, and satisfies for any (x,p) E T
M. i.e. the cotangent space. For any z E M.
dvolT (z,p) = dædp
where dædp is dual coframe for T.M.
Therefore, under AI, identifying T*M with x RR the volume form on T induced by g" the Lebesgue measure of Rd x R restricted to T*M. Despite adopting at first Riemannian perspective on M, we actually recover the Euclidean setting on T*M, which motivates us to consider T*M as our sampling space. We refer to Appendix B for details on the cotangent bundle and its metric g*.
Hamiltonian dynamics. Note that with the previously introduced notation, и can be expressed as
dr//voll (x:) exp(- V(x) log(del (g(x))))
This motivates the introduction of the following Hamiltonian on T*M
(x, (=x) log (det g(z)) + 2 [
(1)
In this definition, we notably take into account the scalar product (+, ())))(()- on T,M. Finally, we consider the joint distribution on T*M
d= (1/2) expl- (oooo
where Z = --T-M exp!- ( pl) p), for which the first marginal is π. Indeed, for any E C(M, RR,, we have
JT-M ((x)dd(( (z, -) JT (2))dd(()) 0, Id)dp - JM v(x)d=(x))
where we denote by N₂ (0. Iu) the centered standard Gaussian distribution w.r.t. || -lgga))-1 The Hamiltonian dynamics for H is given by the coupled Ordinary Differential Equations (ODEs)
i = (x,
أ -= -= 2-,H (x,p),
where the derivatives of H can be computed explicitly H (x. = g(x) Pp
0, (,H(x,p) De(x) "p.. pp) I(x),
where L(x) = VV - Dg(x).
2.2 Seff-concordance and regularity
Until now, we have considered an arbitrary Riemannian metric g. In the rest of this work, we focus on metrics given by Hessian of self-concordant barriers.
Self-concordance. We first introduce self-concordant barriers, family of smooth convex functions which are well-suited for minimization by the Newton method.
Definition (Nesterov & Nemirovskii (1994)). Let U be non-empty open convex domain in Rd A function & : - R is said to be a u-self-concordant barrier (with v 2 1) on if satisfies: (a) & (U, R) and & is convex,
(b) d(x) too as
OU.
(c) (z) h, M|] < 2| alay for any M, h he (d) [h]() < w|/h for any x M, h Rd
where g(z) = (z)
Balls for |I=Il a()) called Dikin ellipsoids, key for the study of self-concordance, Appendix C.
3

---

Regularity. The property of a-regularity for some 2 1 is shared by many self celf-concordant barriers, including logarithmic and quadratic programming barriers. ensures stability for interior point polynomial-time methods. Definition and properties of o-regularity are recalled in Appendix C.
Example of the polytope. Let us assume that M is the polytope M = {x Ax < b), where E Rmndd and b E R". We endow with the Riemannian metric g(z) D2)(() where 6 + R is the logarithmic barrier given for Ei=
any M by <(x) = E log (b: - A æ). The barrier Φ is both m-selff-conoordant bar- rier and 2-regular function, (Nesterov & Ne- mirovski, 1998, page 3). Moreover, we have for any x M, ((x) A S(z)))=²AA where ((x) Diag (b, AJ 2) sEjm] We provide in Figure an illustration of this barrier.
0
Figure Á polytope M C R with three Dikin ellipsoids {y TR y "g((yy < derived from its logarithmic barrier.
3 The n-BHMC algorithm
In practice, it is not possible to exactly compute the Riemannian Hamiltonian dynamics (3). We thus introduce numerical version of BHMC (n-BHMC), in which replace the continuous ODE integration by a symplectic numerical scheme. We first define the Hamiltonian integrators used in n-BHMC in Section 3.1 and provide details on the different steps of our algorithm in Section 3.2.
3.1 Hamiltonian integrators of n-BHMC
In the same spirit as Shahbaba et al. (2014), we first rewrite the Hamiltonian given in (1) as H H1 + H2, where we highlight the non separable aspect of H in H2, i.e., for any (z,p) TMM
H1 p) = ((x) + lo))) H2(x, p) = 1/2 [pllg(a) g(a)-1-
Therefore, we have
(z,,(z (x) -=) ( Dg(z).
H =
(4)
0₂ HzH,( = - (g()) (px)
HgH2(z = -(x)
(5)
define implicit scheme to discretize the Hamiltonian dynamics (3) by splitting it into the Hamiltonian dynamics (4) and (5). In the rest of this section, we consider some step-size h E R.
Leveraging the separable structure of H1,
Explicit integrator of Hi. We first approximate the dynamics (4) on a step-size h/2 using a first-order Fuler method (Hairer et al., 2006, Theorem VI.3.3.). Since Hi is separable, this integrator simply reduces to the map Sn/2 M + T*M M defined by
S (I. p) = (z, (z.H)
We have: (i) (A5/2(TTMM C T*M, (ii Say2 is symplectic (we refer to Appendix F for more details on symplecticity), (iii) S_1/2 o S,/2 Id and (iv) S_h/2 = so Sn/2 s. Note that S1y2 is an involution on T M, which inherits from properties (i) and (ii) of S.,/2-
Implicit integrator of H2. Since H2 is not separable, a common discretization scheme such as Euler method is not symplectic anymore, which requires to use implicit method instead. We thus approximate these dynamics on a step-size with a symplectic second-order integrator denoted by focus the Störmer- Verlet scheme, (Hairer et al., 2006, Theorem
Gg. For theoretical purposes, VI.3.4.), also known as the generalized Leapfrog integrator, which is widely used in geometric integration, (Girolami & Calderhead, 2011; Betancourt, 2013; Cobb et al., 2019; Brofos & Lederman. 2021a). We refer Appendix F for discussion other numerical schemes. For any 2!!) E T (z C. T"M consists of points (1)) = (1(1), (!)) which are solution of
(1/2) ()) 20 , (12)).
(1) (0) + 110 (0)
.p (222))) + Hz(() (002))p)01))))))) p0/213)]
(()) -(112) 0 H2 (2 (1221)
(6)

---

As defined, there is no guarantee that (6) admits a (unique) solution. As matter of fact, the integrator GA, can be seen as a set-valued map. Moreover, it is easy to check that: (i) GA. (T*M) C 2 M. (ii) G_A = Gj, o s, since (s)))(((=)) -22))((2 and o₂ 0zH))(((2)) O2 H2(2), and (ii) if [G) > 0 then E (G-n (z).
Implicit integrator of H. Relying on the integrators Sp.2 and G2, previously defined, we now explain how we approximate the Hamiltonian dynamics (3) on step-size h. We first define the set-valued maps FA : T + 2T and RA T + 2 by
F = as, RA = (s S F Sn/22)
Using properties (ii) and (iii) of Gr. we have for any 2 E T such that [Fx(2]] 0, (FA (=), and for any z' E T M such that [ S/
> O, we have 2' E (R, o (z') Thus, FA. and RA are involutive in the sense of set-valued maps. Note also that s RA = S Gr, o S,/2 boils down to the Strang splitting (Strang. 1968) of the dynamics (3) and thus, is a symplectic scheme approximating the dynamics of Hamiltonian H/ a step-size h.
Numerical integrators. In practice, we do not have access to Fù but approximate it with a nunerical map n. For clarity's 's sake, we denote by doms, C. T M the domain of this integrator with 4, (domo,) ) C M. In practice, dome, corresponds to the set of points for which the numerical integration of Fú, outputs a solution. We also approximate RA, with the numerical map RØ (s (doma.)) C T*M M + T*M
R (s o Sn/2) Lb (s Sp/22)
Similarly to RA, s o R approximates the dynamics of Hamiltonian H on a step-size h. In our experiments, we design 4 using fixed-point solver with given number of iterations following Brofos & Lederman (2021a,b). We refer Appendix K for details on computations of /.--
3.2 Steps of the algorithm
For any = = (I. p) E TM,, we define on T M the norm 1 by
(a p')l== p + +Ilp' [p)))l) V(=',p) E T*M,
(7)
where || 1-11((e)-1 is the canonical norm on T induced by the Riemannian metric g (see Section 1) and || (la(e)) is the common norm used on M to study properties of self-concordance (see Section 2.2). As defined, this norm will be crucial to correct the bias of the numerical integration in n-BHMC, presented in Algorithm 1, for which we are now ready to detail the steps.
Steps 1 and 5: applying a momentum update. Assume that the current state at stage n 2 1 is 1x,-1, Pa-1) E M. Then, the momentum is first partially refreshed such that Pn <- + vBG where Gn (O, La) is independent from P This update is applied both at the beginning and the end of the iteration, similarly to Lelièvre et al. (2022).
Step 2: solving a discretized version of ODE (3). Starting from (X,, Pt) = (Xn-1, Pr), we now approximate the dynamics of H on a step-size h with the integrators presented in Section 3.1, while ensuring that the proposal state belongs to T M. We proceed as follows:
1. We first run the explicit integrator Sa,2 and compute Z (s (x.,,)
2. If 2!!0) & dom.s,. then (X,,,K)) PK)) is not updated and we directly go Step 3. Otherwise, we define z = =Z,(Z()) To ensure the reversibility of this step, we then perform an "involution checking step". (Zappa et al., 2018: Lelièvre et al., 2019), highlighted in yellow in Algorithm ie., we verify (a) z,(") E domæ, and (b) 2 $ i,, (z())). In practice, we replace condition (b) by combining explicit tolerance threshold with the norm defined in (7) and accept the proposal if
_ 4n (Z( || g0 + 0 2n 11 (2!") S 1
(8)
We discuss the choice of this norm to perform the "involution checking step" in Appendix K. If both of
these conditions are satisfied, we finally set (X,,P)) P²)) = (s o S2 (2 ), to ensure the reversibility of the update of (X,, P.).. Otherwise, (X,, P() is not updated and we go to Step 3.
Note that R is function and not a sel-valued map.
5

---

Algorithm 1: n-BHMC with Momentum Refresh
Input: (Xo, Pb) T E (0, N 0.
2% with domain dome.
Output: (X Pn)nE[N]
for n 1,.., N do
Step 1:
Pa
BBG,
Step 2: solving a discretizen zed version of ODE (3)
X,, P. < Nn-1,PP. Pni
X p.o) (s (soSa)2)(Xn:
if 2 =(xi.) )P(o) E dome, then
_ (2()))))
erto 11 4b (2x |1
z020
em | (n(
(l (2!
(if defined)
if - dome, erro er, then X P/ + (s Sp. (Zn))
end
10
Step 3: An + min(1, (xxp_leepp_(xx,PPPP Ph) H P,))) -a~ U
11
Un S An then Na. f- Pa
12
else X,, Pa + -x--1,
13
Step 4: Nn, s(Nn,,
14
Stop 5: N,,00, + VI- 3P. + VBG²
end
Step 3: computing the acceptance filter. We denote by a(x', p' (x. the acceptance probability to move from (z.p) E T to (z',p') E T*M, which is given by
alx', p'' pr,p) p) A expl-H((',))))
(xx
(9)
After Step 2, we perform simple MH filter by accepting the proposal with probability An a(X,, n. Pk X 1...1.P, P,)) similarly to classical HMC algorithm.
Step 4: applying momentum reversal. To ensure a move along the dynamics of H with step-size h, we reverse the momentum after the acceptance step. Indeed, if the acceptance filter is successful, then (X,, Pn) (s (X Pa-1)
approximates the Hamiltonian dynamics (3) step-size starting at (X,-1, Pn- as seen in Section 3.1. Otherwise, (Xn: P.) = -
3.3 About the usefulness of the "involution checking step"
In their paper, Kook et al. (2022a) also incorporate self-concordance into RMHMC to sample from distributions supported on polytopes. To integrate the Hamiltonian dynamics, they propose to use numerical version of the Implicit Midpoint integrator (Hairer al., 2006, Theorem VI.3.5.), which shares the same theoretical properties as the Störmer- Verlet scheme. This results in CRHMC (Kook et al., 2022a, Algorithm 3), whose main difference compared to n-BHMC is the lack of the "involution checking step" (Line 8 in Algorithm 1). Crucially, Kook et al. (2022a) assume that their implicit integrator admits a unique solution given any starting point and any step-size h, which then approximated by their numerical scheme. However, this statement may not be true in the self-concordant setting, as we explain below with simple example. This explains why we always refer to the implicit integrator as a set-valued map. In their setting, the numerical integrator is not guaranteed be involutive, which results in asymptotic bias that the MH filter cannot solve.
In our paper, we do not make such an assumption on the Störmer- Verlee integrator, and consider numerical solver which outputs one approximate solution to this implicit scheme. By doing so, our analysis is meant to be as close as possible to the practical implementation of RMHMC. Then, the "involution checking step"" is critical make the numerical integrator locally involutive, which enables us to derive rigorous reversibility results, see Theorem 3, and to solve the issue of asymptotic bias caused by implicit integration, see Section 6. We refer to Appendix J for detailed comparison between n- n-BHMC and CRHMC. We now present a simple setting where the assumption made by Kook al. (2022a) does not hold in the case of the Soomerr-errlel integrator, which is defined in (6).
Consider the 1 1-dimensional setting M = (-00, 0) with g x + 1/2². Assume that 3 1. Let h > 0 and let Xo E M be the (position) state of the Markov chain at the beginning of the first iteration in n- BHMC. this case, Po N(0, 1/X3) and we have p(0) Po (h/2)0,, H, (Xo, Po). where the function Hi is defined in Section 3.1. Using (4), we have (0()) N(po. 1/ 1/x3), where

---

Mo - (1/ Po) only depends on Xo- Then, the implicit equation in Po p given by
the first equation of (6), reads as p(1/2) = Po (L/2)) PX Pf Using (5), it can be rewritten a polynomial equation of degree in P/1/2)
2 ( P/ 12) P( =
(10)
Denote A 1+ P0 (0) Then, (10) may admit 0,1 or 2 solutions depending on the sign of A. If POo I - 11/(2hXx), i.e., A 2 O. then (10) admits 1 or 2 solutions. However, if h is small enough, only one solution is valid when considering the constraint on the position update (6). Whenever Po -1/(2hX))) i.e. A < 0, (10) admits no solution. Recalling that P/ ~ N(uo. 1/X0), this event occurs with positive probability, thus violating the assumption made by Kook al. (2022a).
4 Theoretical results
We now study the reversibility of n-BHMC with respect to the target distribution. To do so, we first present theoretical results the exact integrators of the discretized Hamiltonian dynamics in Section 4.1, from which we derive our assumption on the numerical integrator used in n-BHMC. Finally, we state our main result on n-BHMC in Section 4.2.
4.1 From implicit to numerical integrators
We show in Proposition 2 that even though F is a ser-valued map, can locally be identified with a c--diffeomorphism. In a manner akin Lelièvre et al. (2022), this justifies assumption A3 that the numerical map 4,, used to approximate F2, in Step 2 of Algorithm is locally a C² -ivvoluution.
Proposition 2. Assume A1, A2. Let (0) E T.M, then there exists h** > 0 (explicit in Appendix G) such that for any h E (0, h"), there exist (1) E F2(0)))), neighborhood U C T of 20) and
C -difeommophhsmm 7th U + (U) C T with
(a) 7h (2(0)) = 2/ and det Jac(Tn) 1. (b) The (2) is the only element of Fz( in Th (U) any = E U.
Proposition 2 shows that while Fn, can take multiple (or none) values on T*M, for any 2!!) E T there exists h small enough and neighborhood U of 20)) such that the set-valued map F2 is locally symplectic U. The proof of Proposition 2 is given in Appendix G. It first relies considering the Störmer- Verlet scheme introduced in (6) as the composition of maps for which we derive the existence of solutions and secondly applying the implicit function thcorem on these maps. Motivated by this result on the implicit map F,, and the fact that for any z E T M with [F [2)] > 0, z E Fn (2) (see Section 3.1), make the following assumption on the numerical map 4n.
A3. There exists X (0, s.r. for any 20) E M, there exists h4 >> O and for any h E (0, hw,
(a) B = BiilIl.co (z(0), Ar (x(( C dom4y.
(b) 4, E (B, M) and 4, 4,, ld on B,
with r": M + (0, +00) depending only on g and defined in Appendix I.
Assumption A3 can be thought as a strengthening of Proposition 2-(a), where (i) ht refers to h+, and (ii) B and 4, correspond to explicit versions of U and 7he We conjecture that the involution condition in A3-(b) could in fact be replaced by the condition that for any 2 E B, (2) e (x),,, in manner akin to Lelièvre et al. (2022). We leave this study for future work.
4.2 Reversibility results
Beyond n-BHMC. While Algorithm can be implemented practically, cannot be easily analysed under A1, A2 and A3. The main reason for this limitation is that the self-concordance properties defined locally, whereas deriving reversibility results require global control. In particular, we cannot ensure that + is locally an involution around z00) if h > ht. To circumvent this issue, we enforce a condition on h to be small enough in n- n-B Using the notation from A3, define the set
A f M h minseB (x.1) (2)) C. domp,
7

---

Let 20() E A,, By A3, we know that 4, is an involution on a neighbourhood of 20);; in particular, it comes that (LA o (2))2 = 20)) Hence, the condition (b) of the "involution checking step" in Algorithm is de facto satisfied. This naturally leads replace "Z%"" Edoms," ", i.e., the condition (a) of the "involution checking step", by the more restrictive condition " "Z"." E A,," as presented in Algorithm (see Appendix H), for which we are able derive a reversibility result.
We denote by Q T X ((T + [0, 1], the transition kernel of the (homogeneous) Markov chain (Xn, Pn NEN obtained with Algorithm 3. We now state main result on n- n-BHMC.
Theorem 3. Assume AI, A2, A3. Then, Q reversible up to momentum reversal, i.e. we have for any f E C(T'M X T R) with compact support
SI-MXI-M [2, 2 (=)=)
JI-MXT-M ((()))))
In particular, it is invariant measure for Q-
Proof. We provide here sketch of the proof, technical details being postponed to Appendix I. The reversibility (up to momentum reversal) of the momentum update is straightforward. To establish the reversibility up to momentum reversal of the numerical Hamiltonian integration step. we cover the compact support of f with a finite family of open balls with respect to the metric g. Combining A2 and A3, we show that 4 is volume- preserring C c'ddiffcomorphism on each one of these sets. We then conclude upon combining this result with the fact that R s Sp/2 = s o Sh,2 4/.
Although Algorithm 3 may be implemented. this would come with huge (and unrealistic) com- putational cost since verifying that = E A,, implies to find the minimal critical step-size h, on a neighborhood of The question of the existence of global step-size h such that Algorithm can be properly analyzed is not the topic of this paper and is left for future work.
Comparison with Theorem 8 in Kook et al. (2022a). Although Kook et al. (2022a) present theoretical result similar to Theorem 3. we claim that their statement is not correct. Indeed, the authors make critical confusion between the ideal integrator as considered in Hairer et al. (2006) and the numerical version they implement. Then, in the proof of reversibility for their scheme, they act if the two algorithms were the same while it is not true (see Appendix J for further details). On the other hand, (i) we make a clear distinction between the ideal integrator and its numerical implementation (see Section 3.1), and (ii) implement the "involution checking step" to enforce reversibility (Line 8 in Algorithm 1).
Related work
Sampling on manifolds. Traditional constrained sampling methods in Euclidean spaces include the Hit-and-Run algorithm (Lovász& Vempala, 2004), the Random Walk cetropolis-Hastings (RWMH) algorithm, also referred to as Ball Walk (Lee & Vempala, 201 7a), and HMC, (Duane et al., 1987). However, it has been empirically demonstrated that RWMH and HMC require small step-sizes in order to correctly sample from target distribution over submanifold M C Rdd thus resulting in poor mixing time (Girolami & Calderhead, 201 1, Figures 1 and 3). In the specific case where M = {x E IRd c(x) 0} for some 2d + " Brubaker et al. (2012) combine HMC with the RATTLE integrator (Leimkuhler & Skeel, 1994), incorporating the constraints of M in the Hamiltonian dynamics via Lagrange multipliers. Girolami & Calderhead (2011) adopt an original approach by endowing M with a Riemannian metric g and propose RMHMC (Riemannian Manifold HMC), a version of HMC where the Hamiltonian depends Hamiltonian dynamics of (3) on short-time steps using the generalized Leapfrog integrator (6) and the acceptance filter defined in (9). This method does not include an "involution checking step". and thus the reversibility of the algorithm cannot be ensured in practice and in theory. Similarly, Byrne & Girolami (2013) propose design numerical integrator for RMHMC using geodesics.
g in (1). It consists of integrating the
Choice of the metric in RMHMC. There are various ways to design a meaningful metric g given a submanifold M. In a Bayesian perspective, Girolami & Calderhead (2011) 1) aim at computing the a posteriori distribution of statistical model of interest and thus choose g to be the Risher-Rao metric. In this paper, we consider another approach which exploits the geometrical structure of M. For instance one can always define self-coneordant barrier e when M is a convex body, (Nesterov & Nemirovskii, 1994), and choose g D26 as done by Kook et al. (2022a).
8

---

Self-concordance in RMHMC. Elaborating on Riemannian geodesics, Lee & Vempala (2018)
provide theoretical guarantees of fast mixing time for RMHMC when (i) M is a polytope, and (ii) is the Hessian of self-concordant function Φ (as in our setting). Their result notably improves the complexity of uniform polytope sampling from algorithms relying on self-concordance such as the Dikin Walk (Kannan & Narayanan, 2009) and the Geodesic Walk (Lee & Vempala, 2017b). However, they only consider the case where the exact continuous Hamiltonian dynamics is used, as in c-BHMC (see Appendix D). On the other hand, Kook et al. (2022a) integrate the Hamiltonian dynamics via implicit schemes without any "involution checking step" in a similar self-concordant setting. In particular, they consider convex body K equipped with self-concordan barrier ø. combined with a linear equality constraint Ar = b. This is a special case of our framework (see Al and 22) by rewriting this whole set K' {A u : the
Ker(A)}², equipped with barrier u -> @(Atb + u.. Besides this, Kook et al. (2022a) provide an efficient implementation of their algorithm (CRHMC) in the case of a convex bounded manifold of the form M M- {z E Rd Az 6, l < < uj. Although CRHMC demonstrates empirical fast mixing time, we show in Section that suffers from an asymptotic bias. More recently, Kook et al. (2022b) built upon the empirical results of Kook et al. (2022a) to derive theoretical results of fast mixing time. However, they do not question the issue of asymptotic bias in CRHMC.
Enforcing reversibility. Zappa et al. (2018) propose version of RWMH including projection steps after each proposal. To our knowledge, they are the first to practice "involution checking" of the proposal, and thus enforce the reversibility of the Markov chain with respect to æ. Lelièvre et al. (2019) notably combine this method with the discretization suggested by Brubaker et al. (2012) to also enforce the constraints of the manifold. Lelièvre et al. (2022) elaborate this framework, by designing a symplectic numerical integrator with multiple possible outputs, and provide a rigorous proof of reversibility. Note that it only applies when g induced by the flat metric of Rd and therefore cannot be combined with our approach as such.
6 Numerical experiments
In our experiments", we illustrate the performance of n-BHMC (Algorithm 1) to sample from target distributions which are supported on polytopes. We compare our method with the numerical implementation of CRHMC provided by Kook al. (2022a). In all of our settings, we compute g as the Hessian of the logarithmic barrier, see Section 2.2. The algorithms are always initialized at the center of mass of the considered polytope. At each iteration of n-BHMC, we perform one step of numerical integration, using the Störmer- Verlet scheme with K = 30 fixed-point steps and keep the refresh parameter 3 equal to 1. We refer to Appendix K for more details on the setting of our experiments and additional results.
Synthetic data. We first consider the problem of sampling from the truncated Gaussian distribution
(dr/dLeb)(x) expl- ||. - --1l² /2] ( JM expl || 11l² /2]dē, Vo E Ed, where M is the hypercube or the simplex, and μ E dR is defined by
/ (10/V x (1-
+ (vd-1 - 1)e2)
where e, stands for the i-th canonical vector of R and 1 1 Σi e,. In particular, =1= = 0. #22 = 10 for any j E (3.....d d). Therefore, the mass of ι is not evenly distributed on the boundary of M, which is key to observe the impact of the reversibility condition. For this experiment, we consider the fow-dimensional setting d E (5, 10). To sample from т, we run 10 times the algorithms CRHMC and n-BHMC for N- 10s iterations, and update the step-size h such that the average acceptance rate in the MH filter is roundly equal to 0.5, following Roberts & Rosenthal (2001). We discuss the setting of the tolerance parameter 77 for B-B in Appendix K. In order to correctly assess the bias of each numerical method, we aim at accurately computing the target quantity Q = fm (z, p)2d (x). Although our choice of Q is arbitrary, it highlights well the bias in CRHMC, which is corrected when using n-BHMC.
and (4 = P2//V
3 At is the pscudo-inverse of A. 1 code:
/ http://////iiiooo/////////aa///a
connanenont-ca

---

d =
5
d
We evaluate their performance by taking as ground truth the esti- mate given by (i) the Metropolis- Adjusted Langevin Algorithm (MALA) (Roberts & Stramer, 2002) for the hypercube and (ii) the Independent MH (IMH) algo- rithm (Liu, 1996) for the simplex, which are also run 10 times, but 10 times longer than n- n-BHMC and CRHMC, i.e., for N 10 iterations. To compute the target quantity Q, we keep the whole trajectories and report the confi- dence intervals in Figure 2. For both polytopes, we notably ob- serve that our method is less bi- ased than CRHMC.
6.9
9.2
6.7
8.8
8.6
CRHMC
n-BHMC
MALA
CRHMC
n-BHMC
IMH
10
=
10
10.6
11.5
10.4
11.4
11.3
10.2
11.2
11.1
10
11
CRHMC
8-BHMC
MALA
CRHMC
n-BHMC
IMH
Figure 2: Comparison between n-BHMC and CRHMC on the hypercube (left) and the simplex (right).
Real-world data. We then consider 10 polytopes given in the COBRA Toolbox v3.0 (Heirendt et al.,
Model
d
NNZ
n n-BHMC
CRHMC
2019), which model molecular systems, and follow the method provided in (Kook et al., 2022a, Ap- pendix A) to pre-process them. Here, aim uniformly sampling from these polytopes. How- ever, we do not have access to a realistic baseline that yields an unbiased estimator, since the sam- pling dimension is too high and running MALA or IMH would be too costly. Hence, instead of assess- ing the bias of the algorithms, we rather want to highlight that the "involution checking step" does not hurt the convergence properties of BHMC com- pared CRHMC. We evaluate the efficiency of the algorithms by computing the sampling time per effective sample (in seconds), defined as the total sampling time until termination divided by the Effective Sample Size. We set the initial step-size to 0.01 for both algorithms and 77 = 10 in n- BHMC. We then follow the exact same procedure as in (Kook et al., 2022a, Table 1) by drawing 1000 uniform samples with limit rumning time set to day. Results are given in Table 1. For each molecular model, we specify by NNZ the number of non-zero entries in the matrix A defining the corresponding pre-processed polytope, which thus reflects its complexity. In particular, the sampling dimension corresponds to the number of columns of A. Note that the sampling time values are not of interest in themselves, but are only meant to compare CRHMC and n-B. While it is clear that adding an "involution checking step" implies a trade-off between accuracy and complexity of the method, our results demonstrate that it does not penalize BHMC in the considered settings, and may even make it more efficient in some cases.
ccoli
95
291
0.08822
5.719
cardiac-mit
220
228
0.1905
5.304
Aci-D21
851
1758
141.9
46.17
Aci-MR99
994
2859
325.1
61.34
Abi-49176
1069
2951
179.4
74.22
Aci-20731
1090
2946
0.7124
88.82
Aci-PHEA
1561
4640
2.343
100.8
iAF1260
2382
6368
537.9
170.5
iJO1366
2583
7284
537.9
169.4
Recon
3742
8717
3.740
3280
Table 1: world data setting: comparison
with Kook et al. (2022a).
7 Discussion
In this paper, we introduced a novel version of RMHMC, Barrier HMC (BHMC), which addresses the problem of sampling from distribution over constrained convex subset M Rdd equipped with a self-concordant barrier ạ). Our contribution highlights that the use of well known implicit schemes for ODE integration combined with these space constraints leads asymptotic bias when comes to their numerical implementation. We propose to solve it in straightforward manner via our algorithm, n-BHMC, which relies on an additional "involution checking step".. Under reasonable assumptions, our theory shows that this critical step removes the asymptotic bias. This result is supported by numerical experiments where we highlight this lack of bias compared to sftte-o--hh--ad methods. In future work, we would like to investigate the "coupled" behaviour of the hyperparameters h and in practice, the study of irreductiblity of n-BHMC and the influence of rj on the bias. Moreover, we are interested in the application of n- n-BHMC for efficient polytope volume computation. More generally, we are convinced that our study first step for future work on designing implicit integration schemes for unbiased constrained sampling methods.
10

---

Acknowledgments
AD acknowledges support from the Lagrange Mathematics and Computing Research Center. AD and MN would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme The mathematical and statistical foundation of future data-driven engineering when work on this paper was undertaken. MN acknowledges funding from the grant SCAI (ANR-19-CHIA-0000)
References
Betancourt, M. general metric for Riemannian manifold Hamiltonian Monte Carlo. In International Conference Geometric Science Information, pp. 327-334. Springer, 2013.
Betancourt, A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint
arXiv:1701.02434, 2017.
Brofos, J. and Lederman, Evaluating the implicit midpoint integrator for Riemannian Hamil- tonian Monte Carlo. In International Conference on Machine Learning, pp. 1072-1081. PMLR, 2021a.
Brofos, J. A. and Lederman, R. R. On numerical considerations for Riemannian Manifold Hamiltonian Monte Carlo. arXiv preprint arXiv:2111.099 10.9995 2021b.
Brubaker, M., Salzmann, M., and Urtasun, R. A family of MCMC methods on implicitly defined manifolds. In Artificial Intelligence and Statistics, pp. 161-172. PMLR, 2012.
Bubeck, S., Eldan, R., and Lehec, J. Finite-time analysis of projected Langevin Monte Carlo. In Proceedings the 28th International Conference on Neural Information Processing Systems, NIPS' 15, pp. 1243-1251 Cambridge, MA, USA, 2015. MIT Press. URL http: htt.://....cc..rr.//// citation.cfie 2692 22w..999....99....
Byrne S. and Girolami, M. Geodesic Monte Carlo on embedded manifolds. Scandinavian Journal of Statistics, 40(4)):25-845, 2013.
Cobb, A. D., Baydin, A. G., Markham, A., and Roberts, S. Introducing an explicit symplec- tic integration scheme for Riemannian Manifold Hamiltonian Monte Carlo. arXiv preprint 1rXiv:1910.06243, 2019.
Cousins, B. and Vempala, S. Gaussian Cooling and (*(n"3) Algorithms for Volume and Gaussian Volume. ArXiv e-prints, September 2014.
Dalalyan, A. S. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society, Series B, 79(3)x651-676, 2017.
Douc, R., Moulines, E., Priouret, P., and Soulier, P. Markov Chains. Springer, 2018.
Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. Hybrid Monte Carlo. Physics letters B, 195(2):216-2))) 2222, 1987.
Durmus, A. and Moulines, E. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. The Annals of Applied Probability, 275))1151 06 2017. doi: 12.1214 AAP1238.
Dyer, M. and Frieze, A. Computing the volume of convex bodies: where randomness provably helps. Probabilistic combinatorics and its applications, 44:123-170, 1991.
Gelfand, A. E., Smith, A. F., and Lee, T.-M. Bayesian analysis of constrained parameter and truncated data problems using Gibbs sampling. Journal of the American Statistical Association, 87(418): 523-532. 1992.
Girolami, M. and Calderhead, B. Riemann Manifold Langevin Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology). (22))123-214, 201

---

Hairer, E., Hochbruck, M., Iserles, A., and Lubich, C. Geometric numerical integration. Oberwolfach Reports, 3(1):805-8822 2006.
Heirendt, L.,, Arreckx, S., Pfau, T., Mendoza, S. N., Richelle, A., Heinken, A., Haraldsdóttir, H. S., Wachowiak, J., Keating, S. M., Vlasov, V., et al. Creation and analysis biochemical constraint- based models using the COBRA Toolbox v. 3.0. Nature protocols, 14(3):639-777 2019.
Kannan, R. and Narayanan, H. Random walks on polytopes and an affine interior point method for linear programming. In Proceedings of the forty-first annual ACM symposium on Theory of computing, pp. 561-570, 2009.
Kook, Y., Lee, Y.TT., Shen, R., and Vempala, S. Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Advances in Neural Information Processing Systems, 35:31684-31663 2022a.
Kook, Y., Lee, Y. T., Shen, R., and Vempala, S. S. Condttio--ummber-nneeeeeenn conver- gence of Riemannian Hamiltonian Monte Carlo with numerical integrators. arXiv preprint arXiv:2210.07219, 2022b.
Lan, S. and Shahbaba, B. Sampling constrained probability distributions using Spherical Augmenta- tion. ArXiv e-prints, June 2015.
Lee, J. M. Riemannian Manifolds: an introduction to curvature, volume 176. Springer Science Business Media, 2006.
Lee, Y. and Vempala, S. Eldan's stochastic localization and the KLS hyperplane conjecture: an improved lower bound for expansion. In 2017 IEEE 58th Annual Symposium Foundations of Computer Science (FOCS), pp. 1007. IEEE, 2017a.
Lee, Y. T. and Vempala, S. S. Geodesic walks in polytopes. In Proceedings of the 49th Anmual ACM SIGACT Symposium on theory of Computing, pp. 927-940, 2017b.
Lee, Y. and Vempala, S. S. 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, pp. 1115-1121, 2018.
Leimkuhler, B. J. and Skeel, R. D. Symplectic numerical integrators in constrained Hamiltonian systems. Journal of Computational Physics, 112(1):117-225, 1994.
Lelièvre, T., Rousset, M.. and Stoltz, G. Hybrid Monte Carlo methods for sampling probability measures on submanifolds. Numerische Muthematik, 143(2):379-421, 2019.
Lelièvre, T., Stoltz, G., and Zhang, W. Multiple projection Markov Chain Monte Carlo algorithms on submanifolds. IMA Journal of Numerical Analysis, 2022.
Lewis, N. E.. Nagarajan, H., and Palsson, B. O. Constraining the metabolic senotype-phenotye relationship using a phylogeny in silico methods. Nature Reviews Microbiology, 10(4):291-305, 2012.
Liu, J. S. Metropolized independent sampling with comparisons to rejection sampling and importance sampling. Statistics and computing, 6:1 6:113-119, 1996.
Liu, J. S. and Liu, J. Monte Carlo strategies in scientific computing, volume 10. Springer, 2001.
Lovász, and Kannan, R. Faster mixing via average conductance. In Proceedings of the thirty-first annual ACM symposium on Theory of computing, pp. 282-287. ACM, 1999.
Lovász, and Simonovits, M. Random walks in a convex body and an improved volume algorithm. Random structures & algorithms, 4(4):359-412, 1993.
Lovász, I.. and Vempala, Hit-and-run from a corner. In Proceedings of the thirty-sixth annual ACM symposium Theory of computing, pp. 310-314, 2004.
12

---

Lovász, L. and Vempala, S. The geometry of logconcave functions and sampling algorithms. Random Struct. Algorithms, 300))307-338, May 2007. ISSN 1042-9832. doi: 10.1002//sa.v000 URL http: wwt.dx..ddo.oo..... http://tdtt..ooo1.100000000 00t0000.0000.0...300000 3.
Mok, K.-P. Metrics and connections the cotangent bundle. In Kodai Mathematical Seminar Reports, volume 28, pp. 226-238, Department of Mathematics, Tokyo Institute Technology, 1977.
Morris, B. Improved bounds for sampling contingency tables. Random Structures Algorithms, 21(2):135-146, 2002.
Neal, R. M. et al. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2 (11):2, 2011.
Nesterov, Y. Introductory lectures on convex optimization: basic course, volume 87. Springer Science & Business Media, 2003.
Nesterov, Y. and Nemirovski, A. Multi Multi-parameter surfaces of analytic centers and long- step surface- following interior point methods. Mathematics of Operations Research, 23(1:11-38, 1998.
Nesterov, Y. and Nemirovskii, A. Interior-point polynomial algorithms in convex programming. SIAM, 1994.
Pakman, A. and Paninski, L. Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians. Journal of Computational and Graphical Statistics, 23(2):518-542, 2014.
Roberts, G. O. and Rosenthal, S. Optimal scaling for various Metropolis- algorithms. Statistical science, 16(4):351-367, 2001.
Roberts, G. O. and Stramer, O. Langevin diffusions and Metropolis-Hasting algorithms. Methodology and computing in applied probabilityy 4:337-357. 2002.
Sasaki, S. On the differential geometry of tangent bundles of Riemannian manifolds. Tohoku Mathematical Journal, Second Series, 10(3):33-334,, 1958.
Shahbaba, B., Lan, S., Johnson, W. O., and Neal, R. M. Split Hamiltonian Monte Carlo. Statistics and Computing, 24(3)339-339, 2014.
Strang, G. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506-517, 1968.
Tao, M. Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm and long time performance. Physical Review E, 94(4):043303, 2016.
Thiele, I., Swainston, N., Fleming, R. M., Hoppe, A., Sahoo, S., Aurich, M. K., Haraldsdottir, H., Mo, M. L., Rolfsson, O., Stobbe, M. D., al. A community-driven global reconstruction of human metabolism. Nature biotechmology, 3 1155441194425, 2013.
Zappa, E., Holmes-Cerfon, M., and Goodman, J. Monte Carlo on manifolds: sampling densities and integrating functions. Communications on Pure and Applied Mathematics, 71(12):2609-2 2647 2018.
13

---

Organisation of the supplementary
This appendix is organized as follows. Appendix A summarizes general facts that will be useful for proofs. Appendix B and Appendix C provide additional details Riemannian metries and technicalities on self concordcce,,, respectively. The c- c-BHMC algorithm is presented in Appendix D with some results whose proofs are given in Appendix E. Appendix presents general facts about numerical integration and specific facts on the integrators used in n- n-BHMC. Appendix G dispenses the proof of the result implicit integrators of n-BHMC, stated in Section 4.1. Appendix H presents the modified version of n- n-BHMC, incorporating a step-size condition, for which we state the reversibility with respect to the target distribution in Section 4.2. Proof of this last result given in Appendix A detailed comparison between n- n-BHMC and CRHMC (Kook et al., 2022a) is given in Appendix J. Finally, Appendix K provides more details on the experiments of Section
Notation. For any Riemannian manifold (M, g) and any Z (ar, p) T we denote by || - I1s the norm (7) defined on T'M by [1<2', || [(x) + |P'll [(xx) for any (''p') E T For any open subset U c Rd and any N, we denote by C CU, R"") the set of functions f U + Rd such that f is k times continuously differentiable.
A Useful facts and lemmas
We recall here some basic knowledge on linear algebra and probability, and state useful inequalities for our proofs.
Linear algebra reminders. For any matrices A, B E RRdxd we write A < B if for any x E g *TB A)æ 0. Any positive-definite matrix A Rdxd induces scalar product () on Rd defined by (z. (x, Ag). This scalar product induces the norm || IIA on R defined by ||^|| 2) = ||A ||A//22|| 2. For any positive-definite matrices A, B E iddd and for any a 2 3 aB is equivalent to || |a S vall- [|B. The canonical norm of Rd induces the norm || ll2 on Raxd defined by || || sup[ || ||Au||2- In partieular, for any matrices A, B R and any vector x E Rd, ||Az||2 I |||||||=/|| and ||AB|| S || Bl|2. Moreover, if A is non-negative- definite, then |||||2 = X(A), where X(A) is the largest eigenvalue of A, and ||A [|2 = X(A) for any a > 0.
Probability reminders. In this section, we consider a smooth manifold M C Rd. We first recall the definition of reversibility (Douc et al., 2018, Definition 1.5.1) before stating general results on reversibility.
Definition 4. Let Q T*M x B(T*M) + [o, 1] be a transition probability kernel and let i be a probability distribution on T M. Then,
(a) Q is said to be reversible with respect to स for any f C(T x T,M, R) with compact support
SI-MxT-M [2, (=)(() = -TMMTTM (()
(b) Q is said to be reversible up to momentum reversal with respect to for any f E C(T
T M, R) with compact support
ST-MXT-MM f(z (2," S...MMTTM (()))))
where we recall that for any (s, p) E Rdd RRd, s(x,p) = -p)-
Lemma 5. Let Q T B(T [0, 1] be transition probability kernel and let be probability distribution on T M. Assume that u मे and that Q is reversible up to momentum reversal with respect म. Then an invariant measure for Q for any f E C(T*M, R. with compact support
S f(2))(22, de' (dzz) Sr-MxI-M f(2)7(dz).
14

---

Proof. Let f C(TMM, with compact support. We have
J (())((2dd>))
=
(Definition 4)
- (2) (dz) = ST-M f(z)#(dz) (momentum reversal on 2)
which concludes the proof.
Lemma Let Q: T x B(T*)) - [0.1] be a transition probability kernel and let be a probability distribution on T Assume that s Q2 Q and that QQ is reversible with respect n. Then, Q is reversible up to momentum reversal with respect to i.
Proof. Let C(TMM, R) with compact support. We have
ST-MXT-M (2, 20= dz))
-TT-MTT-M f(2,
(momentum reversal on 2')
(Definition 4)
(momentum reversal on z')
= - f(s(2)), (2))
which concludes the proof.
Useful inequalities.
The following inequalities hold:
(a) For any u E [0, 1/5), (1 u)-² < 3u and (1 u) -3 S + 5u.
(b) For any E [0, 1/2], (1 u)- - 1 2u.
(c) For any u E [0, 1]. (1 -2) 11 3u.
(d) For any 2 0, we have 1 - (1 u) S u and (1 u)² - 1 -2u.
Details on Riemannian metrics
Let M a smooth d-dimensional manifold, endowed with a metric g. We recall that the Riemannian volume element corresponding to (M. g) is given in local coordinates by (z) / (gt(())az where dr is dual coframe, (Lee, 2006, Lemma 3.2.). For any T.M and T,M, the tangent space at æ and its dual space, i.e., the cotangent space at æ. Note that T. and M space vectors, and T is endowed with the scalar product (fz))(z)) by definition of the Riemannian metric. For clarity sake, we denote by u (resp. p) an element of the tangent (resp. cotangent) space. We recall that the tangent bundle TM and the cotangent bundle T M are respectively defined by TM = UxeM(x) U T, M and T*M = UzeM(2 U T M. These two sets are 2d-dimensional manifolds.
M, we respectively denote by
Metric on the tangent bundle TM. Sasaki (1958) originally introduced on TM a Riemannian metric which, among other properties, preserves the Fuclidean metric induced by g on each tangent space. This metric is defined by
[i + ot o g و 6=(0))
- gí
uT
where Biy and و respectively refer g and g 1 and I corresponds the Christoffel symbol. Since (TM, g) is a Riemannian manifold, the volume form on TM satisfies for any (x, v) E TM
dvolTM(z, v) det(0ā dddd det(g(z))dzdu.
15

---

Metric on the cotangent bundle T*M. Elaborating on the metric defined by Sasaki (1958). Mok (1977) showed that for any 2 M. T is naturally endowed with the scalar product (+ p0x=)=+ and proposed Riemannian metric g² on M,, which notably preserves this result on each cotangent space. This metric is closely related and defined by
Dij + p 60"
9 T
nd TT
Since (T*M, g²) is a Riemannian manifold, the volume form on T M satisfies for any (z, p) E TMM
dvolT-M p) = v det(gg prdd = dzdp.
In contrast to the tangent bundle, the volume form on the cotangent bundle does not depend on which motivates to augment on T (instead of TM) any target measure π defined on M.
C Properties of seff-concordance
We first recall the definition of self-conccrdance and derive some of its properties in Lemmas 8 and 9. Definition (Nesterov Nemirovskii (1994)), Let be a non-emye non-empty open convex domain in R. function á V + R is said to be u-seeffconcooddmt barrier (with 2 1) on U if it satisfies:
(a) Φ (U, IR) and & is convex,
(b) (0x)) + as au.
(c) [D* h M|]l < | || alay for any x M, h R4,
(d) [M]] - / /||b (ay for any x M, heRY
where g(z) = (zx)
Self-concordance be thought certain kind of regularity. Indeed, if Φ is u-seff-concordant barrier on convex body U, then D2 is 2-Lipschitz continuous on with respect to the local norm induced by (see Property (c)), and more restrictively, is v-Liischitz continuous with respect to the same norm (see Property (d)).
Lemma 8. Let U be a v-self-concordant barrier with g = D2 Assume that bounded. Then, [v(()]l + too and [la(z)]l2 + +00 as x - au.
Proof. Since @ is convex, we have, for any (1, y) E U2
Φ(x)S ((y) + V ( - y)
|0(x)<5 (o()) Vo(z)]2l| v|l
< |o(y)) + Vo(z) diam (U).
where we used Sauchy-Schwarte inequality in the second line. Using Property (b) of d. we obtain that || [[(x)]l + +00 + au. By combining this result with Property (d) of ø. we obtain that [g(~))]2 + too as x aU.
According to Nesterov (2003) (see Lemma 4.1. 4.1.2) Property (c) of self-concordance is equivalent to (D* ho, ha]l 2II [|<ll,(=))
for any M and any h2, ha) E Rdd Rd Rdd
Let u Rd + Rdd be a linear operator. For any & E Id, we define the operator norms [|u]|a(z) and ||u|] x by
[| ((x) sup(u(g(z)A)) [|]]ي(z)) = =}, [|ulla(a)> = sup(u(g(x)~1n) [[hlla(a) 1
In addition, set V c R contains no straight line if for any Dao. = (1ro 1 E R.) with xo * we have y n Dzo * 0.
16

---

Lemma 9. Let Ø: U + R be c-selffconeordant barrier with D² Assume that U contains no straight line. For any a, given by (x,
and any > we denote by W the open Dikin ellipsoid of & at { E Rd |lg xl]a(>)) r). Then, the following properties hold:
(a) For any x E U and any (h1, h2) E Rd Rd
[D' hellla()))))) S 2 2 1100ell
(b) For any x E U, (x, 1) C U, and for any E (x,
(1 - lyy x x 9(x) ( 3 (1 |Il x p(x)) *(x))
(c) For any x E U, any y E (2, 1) and any R the following hold
[Tbll(()) S (1 - lly =l(a(a) ||
[Ablb(()) $ (1 |1y ==(g(=)) || la())
[|bll(()) S 1ly 2 =)(a(m)) [|bbla((>>~
| & (1 - 1ly - x >)(s(>)) -|//l|a((a))>
Proof. The result is direct consequence of (Nesterov, 2003, Theorem 4.1.3, Theorem 4.1.5, Theorem 4.1.6).
We now introduce a-regularity, which slightly strengthens coff-concordance, by ensuring that D a(a + 1) 1)-Lipschitz continuous with respect to the local norm induced by g (see (11)). We state below the definition of a-regularity well as of its properties in Lemma 11.
Definition 10 (Nesterov & Nemirovski (1998)), Let in TR. A self-concordant function & on said a-regular, if E CA( R) and for any x E and any h Rd
and U be non-empty open conver domain
( [ hh h, (.h,M MI S ala 1) a(x) |
|
where |1bllu.z = inf(t t 0, + E U}.
Lemma 11. Let @ U + R be a-regular function with g = D2 Assume that U contains no straight line.
(a) For any E U and any h E Rdd
D'o(z)(h, h hl] s a(a 1)
(11)
(b) For any E U, any (hn, bz ha ha) Rd x Rd x Rd Rd
4 D* h2 ha ha]l S a(a 1) II (hlll
(12)
(c) For any x E U, any y E (x 1) and any h2) E Rd x Rd
[) h2] Dd (² (,, hello sala + 1)/3 ll -lyllly x=(o(z))- 3-3 1). (13)
Proof. We first prove (11). First remark that the result is true for h = 0. Consider now h + 0.
Let > O. We define t2 (1 + $) = p(=) > 0, Ve" = x + tgh and y = I tgh. We have Hence, VE U and v E U by Lemma
|| - y²" (l(())) || ve (la)..) = 1/(1 + + Therefore, |/blluur I t- and
[ h]]l S ala + 1) ((z) (1 + =)²
17

---

We conclude upon taking + 0. Then, (12) is an immediate consequence of (11) using (Nesterov & Nemirovskii, 1994, Proposition 9.1.1). We now prove (13). Using Lemma 9, we have
[ h2] - Palllar
s So |D -y- - s
s (1 t1yy - -
t(u x) 1)) h2 - x:ll ıt
s So a(a =ll |00
s aa(a + [|]|0a||@ac)|1/2211g66 fo lly /lttt/ t|yy - =11(()))
S Sa(a 1)/3 (o2))///lll -Ily --u-=/la -==l((()))) 1)
which concludes the proof.
We now prove the equivalence between the Fuclidean norm and the norms induced by g and g Lemma 12. Let (U.g) be non-empry Riemannian manifold in R.. For any Do E U, any (x.p) U x Rd, we have
(lo(z0) 112 1/21111111111122/11 |||||
| s 1|g(xo))l2 /
llo(z))l1
1/2 lpll S ( S (1(((=))~ l1lll
In addition, let
Cz. = max((1[(z))11 [1l(<o) > 0.
Then, we have
1/Czo (1|==1e + llpll2) S | + [|P1lo(00) S Czp (120112 + llpl12)
D C-BHMC: Algorithm and Results
In this section, we first state general results on the Hamiltonian dynamics (3) under our main assumptions. allows us to introduce the continuous version of BHMC (c-BHMC), for which we derive the reversibility with respect to m. Full proofs of this section provided in Appendix E.
Results on the Hamiltonian dynamics. Using tion 13 the existence and uniqueness of the Hamiltonian dynamics (3), starting from any 2o T.
's theorem, we obtain in Proposi-
Proposition 13. Assume A1, A2. Let E M. Then (3) admits unique maximal solution (Jz, 2), where (i) Jan C R. an open neighbourhood of 0, (ii) 2 E C (Jzg: T*M), (iii) 2(0) 2g and (iv) for any t E Jan H(=(t)) H(=o).
In contrast the Euclidean case (i.e., g(z) I2), the Hamiltonian dynamics (3) relies on metric derived from self-concordant barrier, and thus is defined up to explosion time. In the next result, we characterize the behaviour of the solutions when this explosion time finite.
Proposition 14. Assume A1, A2. Let 20 TMM, (Jao 2) as in Proposition 13. Assume T sup Jzp < +00, then we have
(a) lim,-T,, (lp(t)]22 +00.
(b) There exists E 0M such that lim;-TT, xx() JM-
In the rest of this section, we define the Hamiltonian flow as the map T D - T*M, where D {(h, 2o) T h < sup Jzn), such that Tu(=0) =(h) for any (h, =o) E D, being
uniquely defined from 2g in Proposition 13.
Introduction to c-BHMC. In the case where D = T*M x R, the definition of the Hamiltonian dynamics ensures that T, is involutive (up to momentum reversal) for any h > O, which key to derive reversibility in BHMC. However, in our manifold setting, there is one subtlety that needs to be checked, namely that the Hamiltonian flow does not leave the cotangent bundle in finite time. As detailed in Proposition 14, there is indeed guarantee that this flow is defined for all times, i.e., that
18

---

we have D = T*M x R. due to the ill- ill-conditioned behavior of the dynamics near the boundary of the manifold. Therefore, given a time horizon O, we restrict the cotangent bundle the points where the Hamiltonian flow is defined up to time h. Hence, we define, for any h >> 0, the sets
E { E T : < sup } E, = {=o E En h < sup ((oTs))))))))
Note that elements of E, correspond to points 2p for which:
(a) we are able to compute exactly the Hamiltonian dynamics (3) starting from 2y on the time interval
[0, hl, by considering (TT((e) hrejo,) (since EA Ex).
(b) we are also able to compute exactly the reversed dynamics up to momentum reversal starting from (x (zo) on the time interval [0, hl.
In Appendix E, prove that Tr, is symplectic (up momentum reversal) E,, for any h > which thus ensures that the Hamiltonian dynamics is reversible with respect to म. see Lemma 16.
Let h O. Our algorithm c-BHMC, whose pseudo-code is provided in Algorithm 2, then proceeds as follows at stage n 2 Assume that the current state (Xn-1, Px-1) E T M. First define a partial refreshment of the momentum Pa f- VI- 3Pn + v3G,, where G, ~ (,, 1a) is independent from Pn-1- 1- Then, verify that 1X, P.) E Ex. which ensures that the Hamiltonian flow is involutive the time interval [0, h.. starting this state. In the same spirit as in Algorithm 1, we refer to this step as an "involution checking step", which is highlighted in yellow in Algorithm 2, although the involution property directly derives from the Hamiltonian dynamics. If this step is valid, define (Xn, Px) (Xn-1, Px): otherwise, set (Xn, Pw) + s(Xn-1, P,).) Finally, update the momentum as in the beginning of the iteration to obtain the new state (Xn, Pn).
Algorithm 2: c-BHMC with Momentum Refresh
Input: (Xo, Po) E T 3 e (0, 11, N NN
Output: (Xn; Pa Pa)ne[M]
for n 1,... N do
Step 1: ~ (2,, P. -1-8pu- + VBG
Step 2: solving continuous ODE (3)
if (X) then (Xn, Pw)
else (Xn: Px) P.)
Step 3: G, ~ (2(O,Ta),, Pa + -1 + VBG
7 end
Denote by Qe T x B(T*)) - (o, 1]. the transition kernel of the (homogeneous) Markov chain (Xn, Pa)ne(N] obtained with Algorithm 2. Using the properties of the Hamiltonian dynamics, we get the following result.
Theorem 15. Assume A1, A2. Then, Qe is reversible up momentum reversal, i.e., we have for any f C(T*M x T*M, R) with compact support
S m f(z (=) dz')
In particular, क is an invariant measure for Qc-
E Proofs of Appendix D
E.1 Existence, uniqueness and explosion time of Hamiltonian dynamics
We prove below Proposition 13 and Proposition 14.
Proof of Proposition 13, Assume A1, A2. In particular, M contains no straight line, and Lemmas and 11 apply here. We rewrite the ODE (3) as a Cauchy problem defined on Banach space (R" x R I1 llE) where (a(a, b)llE = ||a||2 + [lbll for any pair (a. E Rd x IR.. We define = M x d which is an open set in E. solution (J, z) to this Cauchy problem is defined for all
t t by
2(t) F(z(t)) 20) where F(z,p) = (0, (0,H((,))) -0,1(z,p))
19

---

such that J is an open interval of R with 0 J, and x = (z,p): J - St is differentiable on J. We now prove two results on to establish existence and uniqueness of maximal solution (J,z) for this Cauchy problem:
(a) F is continuously differentiable on Ω
(b) F is locally Lipschitz-continoous on St with respect to || IlE.
Since VV, Dg E C (M, IR), we directly obtain Item (a). We now prove Item (b). Let (i. EQ We endow E with the norm || I|.. Since M is an open subset of g., there exists 0 < <r (/CCC) where we recall that C2 is defined in Lemma 12, such that (ā, (r) C M. We consider such r, and we define B B|l-ile (z, r) and Bg = B1H-1= (2, where r = Car S 1/11. In particular, we have by Lemma 9-(b) that Bg C Bll-llate) (z. gc (p. F) C 52 since 1. Moreover, using Lemma 12, we have B C Bg C 5.
Since V E C2, R), VV is Lipschitz-contimuuus on Bg. i.e., there exists Cv such that for any (2,z') E Bg x Bg.
VV(x) (V''')) S Crllæ 2' (l((z))
(14)
We first show that F is Lipschitz-continuuus on Bg with respect to || lls. Consider (2,2') E Bg x Bg denoted by z = (x. and = =<''p'). Note that (a, E (" 1) x W 1), where (I, 1) is defined in Lemma 9. We first bound [[F((()) - (()((")) = (10,H(2)
0, (z')lla))))) We have
(z) 0x(1(2)) = - -(a")""pp(z')-ppp - -0x='''0(x)=>-p->>> p²) + (a(x)~ -0(0')-)pp
then,
[z0, - (2'(((2))) S -())))))(((()----)) 12ll=///)
(15)
Using Lemma 9-(b), we have
(1 - || - x (0lo(())) 3 )())) 3 (1 - | - x x)<ا(=>) - -²0())
and thus,
(1 - | - =m(ا(=)) ) - 0(8) 0(x) 0/7) 3 (1- || ==la(a>)) - TT
(1 (1- II - < 0(3) <(1 -1-r) <11-7) -
Therefore, we have
( S (1 -P)-
(16)
In addition, we have
|| x' (b(=)
S (1
-
|| - x' 1l())) - x =1]
< (1- < 1/5,
(17)
where we used Lemma 9-(c) in the first inequality and F S 1/11 in the last inequality. Therefore, x' E W9 1) and we have using Lemma 9-(b)
{(1 - || =0 1))(=')-1 = g(x)-1 - -g(x')1
1 (11 ||a' - ||a'=xl -2- 1)0(x')-1
- x ((1-11==21lota))) -2
- 1) = 1}(1-
-g()))-
1 (11 -1 - x ==(o(=))) -2 - 1] - -1)-
Then.
(o())) (0(=)) - g(x') ()()) 121
(18)
< ma - || - =² - 11(1 -
- 7)2, {11 - || ==lo(z)) -=1( 11(1-))-2)
< (1 - (1)-))=mmx((11
- - x 22 (1 - 1|x' 1z''~xl ==(p((m))) -2 1)
Using Inequalities (a) and (c) with u =|aa- xllo(e). where - 1/5 by (17), we obtain
20

---

(a) [[1 || x <l(g(a)))
1| - 31/2 z)|a(n) < < f) 2 2' 2
(b) (1 - || - x: xllo(x)) < 3|12 - =l) <3 - F) -| r x (l)0(2))
Moreover, we have
1 -
Combining (15), (16) and (18), it comes that
(2)- (
- (1-ド) (1 (1000-0---- (19)
We define Az,1 = (1 - F) 3(1 -) (++ [1 such that (z) - -= S
Ag, || = and we aim to bound (z) (z)) [[,Hl(()) - (2))
We have
0gH(2)) - -1,()) -VV()) - -V(2a)
(20)
+ + "'' - pp,)())"''))))
+ = Dg(x) 1=00=))-+ :Do(a')
We have
())l 'Pllaaa)
s(x')²) p(lp()))))
+ [Pgg()))( +| "/''" (lll
(21)
Note that we have g(x) (2) (resp. g(x')= = (() Using Lemma 9-(c), the first upper bound in (21) is bounded by
(1 - (1=))"'>((<>)>" (3)
(3))
s 3(11=)))| (2) (2)]lla)))
+ - (D (2) (()()))))((())))))))))((())))))) + (()))())))
$ 34(1-r) "
+ + -
(Lemma 9-(a))
<4(1 - (
+2(1 (2)
(Lemma 9-(c))
+ 2(1 F) A 1 1/1 < (1-)) 11)1g)=))))))))) 44 F| 2,
4 -1
-
(see (19))
where we used that (=)) g()) S (1- ((a)-1) )) in the last inequality. Combining Lemma 9-(c) and Lemma (c((c, the second upper bound (21) is bounded by
- (+) = (x'),
- - (1-
s 1) ((1 - ||1 - -0())))))- [| (2))12
sa(a 1)(1 (1(1-7) -)) - (1- (100-) xx1 -3
sa 1)(1 -))1-r)- (F+ ) 1) (1) (2(2)((2)) (z') (1
xx
((1 - (1 - - x )lt)) -3 1}
Using Inequality (a) with (1 (1- (1- x (2)) where 5 1/5 by (17). obtain
- - 1 5(1 -- -) - x ==|((0))
(22)
Then, the second upper bound in (21) is bounded by
(5/3)a(a 1)(1 (((1-))--0 - | |0( =110(0)
21

---

We now define h E M - g(y)- Dg(y). Since Φ E CM)M), h CM)). Moreover, [x,21] x'] M by convexity of M and we can define As supye[z,a1] [D where | xllo) [ Using the mean value theorem on h, we have
[|(=)) h(x') (Iae)lllee )a()))- - Aslla' 2: 20la(a)
(23)
By combining (14), (20), (21) and (23), we have
- F (
$ Cv - 2'' ((
+ (1/2)
1) (F 1 x1
+ +1/2)44lllz x1|
We define Az,2 = Cv + (1 +, (-)
(1/00)
F) (F
|| A such that
(F) - p(2) (=)(lo))) S -2,2|l==2'11
-
Finally, we have
LF(z)-F(z2)) F(z'Ilz = (2) - - (A..1 + As,2) - 2 ||=
- F (z') (2)))l((x)) + (1) F(() (=2
i.e., F is (Aa.1 Aa.2) Lishhhippshh Bg with respect to 1|x-
We now show that F is Lipschitz-coontnnuuss on B with respect to || -|-lle- Consider (2, 2') B denoted by = (z, p) and x'' = (o', ,'). In particular, (2, E B x B Using the previous result and Lemma 12, we have
Ca(A,, (al)F(2)-F(()))
+
which proves that F is (2 + As,2) ) Lostchizz-commtnuuuss B with respect to || -llE.
Conclusion. Combining Items (a) and (b), we obtain the results (i), (ii) and (iii) of Proposition 13
upon using Cauchy-Lipsehiee 's theorem. Moreover, for any t E Jzo- (z(t)) = 0, ((z()) =(t) H, p(() = -i(1) &(1) +()) i()) = O, which proves the result (iv) of Proposition 13. O
Proof of Proposition 14. Assume AI, A2. Let 2g E T (Ja, 2) as in Proposition 13. Assume that theorem, z leaves any compact set of E at the
Tz, = sup J2. < +oo. Then, by neighbourhood of Tz. By construction of || lle, this property induced on both variables x M and p Rd, each with respect to || ||2. We directly obtain the result of Proposition 14 by continuity of (x,p).
E.2 Proof of reversibility in Algorithm 2
In (RM)HMC, one often sets the time horizon of the Hamiltonian dynamics with a hyperparameter h. However, under AI and A2, we are not able to prove that the Hamiltonian dynamics (3) is defined at time h > O, given any starting point 2g M. Actually, this dynamics explodes if the time horizon is finite (see Proposition 14). This theoretical limitation requires to implement some sort of "involution checking step" (see Line 4 in Algorithm 2), which is theoretically explained below.
In the rest of this section, for any 2g E T M, we will denote by 2(. z0) C² T TMM) the maximal solution to (3) with starting point 20, uniquely defined in Proposition 13. For any h >> O, we recall below the definition of the sets E,, and En the map Tr, introduced in Appendix D:
(a) E {2o E T h sup Je
(b) E. = {zo E E, h < sup (soTT
(c) the map TA E, T M is such that Tz(20) 2(h, 20) for any 20 En.
22

---

In particular, we have
E E,. n (s Tx,)
(24)
We first prove that the restriction of o Ta to Ex. is involutive integrator, which is crucial to derive the "reversibility" of Algorithm 2 in Theorem 15.
Lemma 16. Assume A1, A2. Then, for any h > 0. 8 Th is C' "diffeomorphism on En such that
detJaccs Tall
Proof. Assume AI, A2. Let h O. clear that TĄ (and thus s Tn) is well defined and continuously differentiable on E, (in particular Ex), by elaborating on the proof of Proposition 13 combined with Cauchy- Saanch--ippshitz's theorem. Let =g E Eg. By definition of Ex. Tn. (zo), resp. Jzo. and (T, 8 (zo( resp. ((soTsX(=)) are well defined. We now aim to prove that (Tn, 880 Tx)(=o) s(zo).
We define 2' [0, > s(2(h -t, =o)) which existence is straightforward since (h t) Jzp for any E [0, h]. In particular, z'(0) = (s (20) and is a solution of the ODE (3) on the interval [0, kj. Yet, h0 (soTh))(=))) and =0) =00, (soTT, (=o)). Then, by unicity of 21, (so (e) (see Proposition 13), 2(,,((())))) (s (z))) and a' coincide on [0, kj. In particular, we have time h
s(zo) z' z'(h) =hh (s (zu)) (T, o Tz))
Then for any 20 En. (s o Tn o (=o) zo. i.e. 8 Tr involution on Es.. Since o TA E (E,,, T*M), o Th is Ca ddiifeeomoophism on E. such that detJac|s Tx = 1. O
We recall that we denote by Qe M M(TM) - [0, 1]. the transition kernel of the (homogeneous) Markov chain (Xn, Pa)NE[N obtained with Algorithm 2. We also denote by:
(a) Qn T x B(T*)) + [0, 1]. the transition kemel referring to Step I (also Step 3) in Algorithm 2.
(b) Q₂ : T x B(T*M) + [0, 1]. the transition kernel referring to Step in Algorithm 2.
We provide below details on Markov kernels Qo. Qe, and Qe-
Kernel Qo- This kernel corresponds to the Gaussian momentum update with refreshing rate B. For any 2,,2)) E T*M T*M, we have
Qa(2,dz') (p': / Bp, (Td)dp's,, = -2T33 - d/2
(25)
- Bpll (cxa'')
Kernel Qe,1- This kernel is deterministic and corresponds to the exact integration of the Hamiltonian up until time h. For any (z, E T M x T*M, we have
Qe,(z,de = le (2.(())) (2)d)ac))(d=") LEg (=) (=)))(=)(d=').
(26)
Kernel Q.- This kernel corresponds to one step of Algorithm 2 (i.e., comprising Steps to 3). For any (z,²') E T*M T*M, have
Qe d=') ST-MXT-M Qo(2, (e,))Q.. Q2 dz')
(27)
We recall that defined in (2), admits density with respect the product Lebesgue measure given for any z = (x, E T M by
(dw (dr/(ded))) =
p( det(g(x))) expl-V(x))
Lemma 17. Assume A1. Then, the Markov kernel Qo defined (25), reversible (up to momentum reversal) with respect to T.
23

---

Proof. Assume AI. For any (2,2²) TMM x T we have
0,r3h:16666222""""""""4000077443 (cx''
(-)))- - 2 3(p g(=)) + 1pl 1
Qo(=(,dz)7(d2"))
which proves the reversibility of Qo with respect to म. Since h Qo = Qo. we conclude the proof with Lemma 6.
Lemma 18. Assume A1, A2. Then, for any h > O. the Markov kernel Qa,1 with step-size h, defined in (26), is reversible up to momentum reversal with respect to
Proof. Assume AI, A2. We recall that the set E, is defined in (24). Let f C(T*M X T* with compact support. According to Definition 4, we aim to show that
JT-MXT-M (2, Qr, (z, (dz) = ST-MXT-M / (s(=')
(28)
We denote by the left integral of (28). We have / = /1 + 12 where
1 JE, (z, (z)) 12 JEg (2)
We denote by the right integral of (28). By symmetry, we have = J + J2 where
JI je. #(=)[((s +()))
J2 = Jeg (=)s sz))ddz.
We directly have 1 = J2. Let us now prove that I = J1. By change of variable = + (8 (z)
4T we have
4 SE, (2) (z)))
-(OOTW)(E.) ((ss (=)( (T)(s),
(Lemma 16)
= =JE, ((2)f((s sz))
(Proposition 13-(iv) & (s TE,) E,)
-1.
Finally, we obtain / J and thus prove (28) for any continuous function with compact support, which concludes the proof.
We are now ready to prove Theorem 15, which states that Qe is reversible up to momentum reversal with respect to
Proof of Theorem 15. Assume A1, A2. Let f T x T + IR be a continuous function with compact support. We have
M flxx=))00((x.d)2)(()))
(see (27))
= hbtp
(Lemma 17)
(rae (momentum reversal on
= h|tpr:/eehhhhhoolloiaoooooaaoooh (Lemma 18)
hbt.://hhhhciiooooocc.....e.eeor (momentum reversal on 2)
= - (Lemma 17)
= ~/wpgy//e66hhchcc00000ososhahol-

---

Moreover, $ t. Hence, by combining Definition 4 and Lemma 5, we obtain the result of Theorem 15.
F Numerical integration in HMC
In this section, we first recall the definition of symplectic maps, and then discuss the choice of integrators in RMHMC. Finally, focus the implicit integrator defined in (6) and provide some notations, which will be notably used in appendix G.
Reminders on symplecticity. We define Ja = --1.
Ld 0
2dx2d
Definition 19 (Hairer et al. (2006)), A linear mapping A R2dd Red is called symplectic if A JaA = Ja-
Definition 20 (Hairer et al. (2006)). A differentiable map F: Redd where R2d an open set, is called symplectic if the Jacobian matrix Jac[F) J is symplectic, i.e, [[]a[[[((2] [aaca[F[e = Ja for any E particular, if F is symplectic, then det.Jac] [F] =1
Introducing Riemannian metric g in the Hamiltonian priori makes it non separable. Thus, designing an integrator which is (i) symplectic, (ii) reversible and (iii) not too computationally heavy is challenging. The generalized Leapfrog integrator (GLI), or equivalently the Störmer- Verlet integrator, combined with fixed point iterations, as chosen by Girolami & Calderhead (2011), considered as the standard scheme for RMHMC. Brofos & Lederman (2021b) analyze the impact of GLI on the ergodicity of RMHMC under metric which is designed in the same manner as Cobb et al. (2019). In particular, they show that the convergence threshold used to compute the fixed-point process only matters to an extent, and identify diminishing return to using smaller convergence thresholds. They compare the fixed-poite point approach with Newton's 's method, which gives less iterations to converge. On the other hand, the implicil midpoint integrator (IMI) enjoys the same theoretical properties as GLI, when these two integrators are combined with fixed-point iterations. As shown by Brofos & Lederman (2021a), IMI may also show numerical advantages in terms of reversibility and volume preservation for specific target distributions. However, this integrator is hard to use for reversibility proofs due to its non-separability.. Besides this, Cobb et al. (2019) combine 2-state augmented Hamiltonian, based on the work Tao (2016), with Strang splitting (Strang, 1968) to propose symplectic explicit integrator which thus has the advantage not to rely fixed-point iterations. Yet, this integrator does not satisfy reversibility in theory.
Choice of the integrators in RMHMC.
Details on the Störeer-Verlet integrator. The scheme defined in (6) can be written as G₁, G2 o Gj,, where:
(a) G T + 2T M is the set-valued map implicitly defined by the first two aligns of (6), i.e., for
any (2, 2') E T x T z' E (2) if and only if
p' p / H2H,,')
z'=x "10.
(29)
We also define the map g, T T + R" R by
5 (2,2') = (a' x ; 1 (z. p¹) H, p'),,''-p+ -
(30)
such that g, (2,2') 2²) 0 if and only if 2' E G (z).
(b) G T + M is explicitly defined by the last align of (6), ie. for any E T*M,
(2) (z.p h/2 (g,HH(Z.p))
(31)
G Proofs of Section 4.1
In this section, we state several results on the maps defined by the implicit integrators of the Hamiltonian (see Section 3.1), in order to prove Proposition 2.
25

---

Lemma 21. Let U C RR non-empty open set and 2 E U. Assume there exist neighbourhood
U, U of 2, L > O and a L-Lipschitz-contiuuus map vi : U= + RR such that vo E C. (U,,, R.). Then, for any h E (0, 1/L) and any a E RR, the map th : U2 - E defined for any x' Ug by (= + (z')²) + a is a C² diiifeeomorphism on neighbourhood U²2 C U2 of
Proof. Assume we are provided with U, U2, L and da as described in Lemma 21. Let h E (0,
and a E R We define the map th U2 + Rd by (n(2") = x' + hzb((z) + for any o' E Uz. Since 1 E CU Rd), it is clear that 1/1 E C' IR.. We have, for any 2' E U2, Duha(z)== Id +DD(22).. In particular, for any w E R. it comes that
(2)oo)) 2 ||m|| - h = (1 hL)
Since 1 1 hL > 0, Jac(en](s)) is invertible. We conclude the proof by using the local inverse function theorem.
We recall that the set-valued map GR. defined in (6), is the generalised Leapfrog integrator, or equivalently the Störmer- Verlet integrator, of the non-separable Hamiltonian H2 defined in Section 3.1. We state below a result of local smoothness of this integrator.
Lemma 22. Let (M, g) be a smooth manifold of Rd and h > 0. Then, for any (x:10, 2(0) E M x M. the vector p (0) T* 2
M such that (z("), p) E Gz(x=0), (0)) for any p()) T M is uniquely
determined by p()) = Gp, (x:(1)), where the map GR : M + T M is defined by
Gh (v) =0 (11/2(()) (v)
h1
(((x)))))((()))) ( (v), g
(32)
4
x1-)x(0)) x))). In particular, if E (M, ddd), then, for with p(1/2) (x)= = (00(0(0))- + a(v) (a)))))-(() (y -
any x: E M, G E T M).
For any a 2 1, we define ra > 0 by
ra min(1/50000, 1/ 1000a(a 1)3)
(33)
We recall that we denote by Wºr,, r) the open Dikin ellipsoid (w.rt. g) at x E M of radius r > O, given by W where a is given by A2. Then, for any x E Wo r), x E M Jac (o (x) is invertible for any h > O.
{y E Rd lly z =llat=) < r). Assume A1, A2, Let (0) E M and r E (0, ra).
Proof. Let (M.g) g) be a a smooth manifold of Rd and h > O. We first prove the existence and uniqueness of the map defined in (32). Let (I:(0) x())) E M x M. We define for any y E M
2 (1/2) (y) h (o)x + s(v) -1))))))-- (y- - x (0)),
plv) =
=00))p//1)00 (1/21 + (22(2 .p (v))
Note that these expressions are obtained by simply inverting the first two aligns of (6). Then, by definition of Gn, the vector p E T T M such that (g()),p()) E G., (x: .p) for any p()) E T M is uniquely determined by p = ()x(")) We thus obtain (32), using the expression of 02 H2 given in Section 3.1. Moreover, it is clear that Gh is continuously differentiable if g E C2 Rdxd).
We are now going to prove the second result of Lemma 22. Assume A1, A2. Let x0()) E M and r E (0, rāj, where r is defined in (33) and a is given by A2, In particular r < 1/11. For the sake of clarity, we will now denote (()( by go for the rest of the proof.
We define the open Dikin ellipsoid B, (z))(()), .r). Since r 1, By C M by Lemma 9-(b). Let h > O. We now define the following maps on Bg
{ Rd
Bo
+ +st+
Bo Rdd
My =
Kg
+ 90 Mo(y) "yy - x())
V
+ Do
y
Bp
Bp
+ Rd
Φ
y
+ Mo(y) Ko()), Ko(y)
+ G/, (3)
26

---

Note that we have, for any E Bu. o(y) = (y - (v-r(") 2 (v)).. Since is twice continuously differentiable on Bg. the maps previously defined are continuously differentiable. We now compute the derivatives of these maps. Let y E Bo and w E TRd. We have
[Mu())(u) .(v) g(y))-D((y)(wag()))
[Ko
9o
x00) Do Mo(y) -1
90 (o(y)0 + Mo) s(v) + -at- - -x 20
+ + Mo()) '
90 (5(())0' Daa
200) + 9o Mo( -
Doo())[w) =
+ Mo (Ko (DKo) Ko]b
Duo(v)(u]= = La- Doo(u)(w)
We have in particular (x: 0 and (oo(x(0)) 0. Let us first study the smoothness of Dou Bo-
Smoothness of Dōu- We prove here that DΦo is (r%) Lipschitz-cconiinuuus on Bo with respect || ll Let (y: s/) E Bo and w IR... We have
to
[V) [w] (wlll/
(34)
I || DMM (goKK)))) DMo('') Dol Kols/)lllee (w))) Ko(v")lllo
+ [oo (Koo())u) Ko(w) -
[Mo(y) (o (v)
(go(k/)) DKo(y)) (Ko()))([]lle..
Remark that the two last terms in (34) are very similar and can be bounded in the same way. The rest of the proof on the Lipschitz-continuig of DΦo is separated in steps, each of them consisting of bounding from above a term of (34).
Step I. Let us first bound the first term in (34). For any x E Bu. we denote DM (oo(Ko(()) Ko(x)] by as (z). We have
Do 1/2 (a1()) a1()))
(35)
= =(8 (DMM)()))) - DML(s)( 2o 1/2 (v)
+ (
1/2 (y)))(((())))))) /2 1/
1/2 (D (DagoKooo) Ko))) DaolKo(y Ko(r')))
We now aim to bound each one of the four terms that appear in (35).
Step 1.1, First, we have
/ (DM(())) DMo('') 90 - Do)
1/2 - (1= 1/2
0n
90
- + go g(v) (D
1/2
90 p(v) "D((y) - 00v)~ //b//
1/2
We recall that S ra < 1/11, and thus, we have for any x E 1/ the following inequalities
Ilgb 1/2 (a('') ((y) -1 1/2 S 3(1 - (0(1-))~lly w'll 1
(on the model of (18))
1 1/
1/2 ((x))(())))))))) - 1 1 < (1 (1-r)-2,
(Lemma 9-(b))
|lg 1/2 DD(x) 112 < (1 r) -1)0(()))))))))))))))))))) 1/2 D (g(x)(wlo(a)) -1/112-
27

---

In particular, we have for any 2 E {v.3/
(lx()) 1/2 (g())()((()) 1/2 =
sup
( [w,, g(x) 1/2 1/12///111 w'llo(x))
(36)
(WeRR): : 1o 13
2 (|(w|))(0))
sup
[0(x))11/2w -1/2 lo(s))
(Lemma 9-(a))
fwTER* [=]]==
< 2(1 - --- /|/w0ll0
(Lemma 9-(c))
and using again that r S r. - 1/11, we have
(ls()) 1/2 (Dgl')) [w]] Do
(37)
fw'eRd: [|=||2=1} [| g(v ((1')'1//// Do(v)lw, 100) (alt)
= sup
S Sa(a + 1)/3 wl ((1 lly v v'llo(()))3 -3-1 -1
(Lemma 1 11-(c)
S 5aa(a + 1)/3(1 - 1 - 1/1000
(on the model of (22))
Thus, by using (36) and (37), we have
llob 1/2 (DMo w] (M)()))())))) 1/2 - (12(1- - r)- 8 + 55a(a + 11/3(1 -1y - y 311l20 100llo
Step 1.2. Secondly, we have
lle 1/2 D Dgo(Ko(v)) Ko 21]
(Lemma 9-(a))
1/2 ((o))
(definition of Kg)
s 21100 1/2 Moly)
1/2 1
1/2
x (0)
l2 go
90
$ 2r ||00 1/2 Moly) 9o
-1 1/22
where we have by Lemma 9-(b)
90
Mo()) -oo
1/2 1-/2/ = (Ia + 1 1/2 g()) 1/ (1 (1- (1
(38)
Hence, we obtain
[00 1/2 (go s 21 (1 -2)))-2
(39)
Step 1.3. Thirdly, we have
| (M/(s)))) 1/2 1|2 S (1 -r) -))()()))))))))))))))))))))))) 1/2 Da(>D///// (((0')(/)(()))11/3)
(40)
< 2(1 r) (17))/111 -1lwlaa
(using (36))
Step 1.4. Fourthly, have
1/2 (Don Ko(v), Ko(w)) (w)) - Da Ko Ko
(41)
$ (gl - Ko Kn (v) | - Ko()), Kn(v) (/,()),, - Ko (Il))
$ (v) Koly') Kolo')ll (To + ((](()) Ko (y'))))loo Ko(v) Ko(y'))l oloo
(Lemma 9-(a))
s 41|Mo(u) - (0)) - Mo(30) - z())) || MM 'Iy - -x(0)) )
+ +2)lMo(v) -1 - (0)) - Mo(3/) - 11
Bo
In particular, My -y (0)) Mo(v') '(8" - =0))) Mo(y) v) (Mo())
Mo(s/)-1)( =))) and thus we have
[MM(y)) -y - x (0) Mo('') 2 x (0) M|
(42)
s 1 -1/2M((v) -,a0
(Mo(v) Mo(oo) 1/2 1/2 1/2 12lly Vo +1100 1/2 -1/2(Ma(vv)
< (1+ + (1 - -1)r)3)-1111 v'lloo + rl
- 1/2 (Mo(v') Mo(s') 1000 1/
(using (38))
28

---

We now aim to find an upper bound for ll9, (Mo(v) Mo(0')-1)o
1/22 We have
1/2 (Mo(y) -1 -1 -Mo'''' 1/2
8o
=(Ia + 90
1/2 (vv) -1000=0-6 1/2
90
(Ia B0 8(1/) 0(10)) 00/ )
1/2
(Ia + 00 - 1/2
50
+ 80
2 (a(y)) -(1')=1000
000 - (T. 90
1/2 9(10))1000
90 /--1
1// -By//--1/2((L + +By)) (0//)-1//// 1/2
90 {g(y)) 0(1/)
1/2
(o(v))+-((()) '000
1/2 (B2) BalB(y))
where B(y') = It + 90 ((/) 1/2
1/2
90 In particular, we have
(1 + (1 r) 3) 1 3 B(y)) 1/2 3 (1 + (1 - -)2))-1211
by Lemma 9-(b). Note that (17) still holds, where æ'' and F are respectively replaced by y, v' and r. Thus, on the model on (18), we have
1/2 (a(y) - {(1 - llv' - vlla(w)) - 13(1 - r2 3 go
1/2 1 {(1 - -1ly' -llaww)) 1)(1 - r) -1
and then
(1 + 1011(111b -
- -1----------111-------------1---
Ila + B(y) 11/2 1/2 / -(00 // B(y')
90
(1 + - 11l -110 (1
-))
Therefore, we have ||Bo (Mo(())))---MM) 1/2 S (1=((-1)))-==m<((((1)))))))) where
1/3
fi(r) - (1 (1-1(1111) 11 vllooo)) 1)(11 - (1+))))))++(())))---))
f2(r) - (1+ (01-((111)) -1l'->ll 13 -1) (1++1+111)))))))1 (1)))
We now aim to control the upper bound max max (ff(()))
(a) We first bound [fr(r). Since /13/ vlla()) 2 0, is clear that fi) S 0. Using Inequality (a)
with u= 1/ vllo(s). where S 1/5 by (17), we obtain
[f((r)] = (r(
< - (1+3(1 (1+31-))=((+(1-)))--- 11 - vvlo(w)))
< 3(( -r) -1 (1++11111 (1- yllp(e)
(Inequality (d)
3 - -)- + +(1 - (1-r)))-1 yll
(Lemma 9-(c))
(b) We now bound f2(r). We have
[2(r) < (1 - 2(1 (1 + (1 - 2- 11 v -1
< - (1
+ (1))))))))))) 11y' -ylo(w)
< 4(1 - (1 + (1 (1-1)-9)-11 - vlloo
(Lemma 9-(c))
where we used (i) Inequality (d) in the first line with u = l1y vllg(w)) and (ii) Inequality (b) in the
second line with u = 2(1 (1 -1- vllo(() S 4r(1 - (1 (1 (1-))+))) - 1/2
with r $ 1/11.
We recall that r E (0, 1), and thus, we have respectively (1-r)- > (1 (1-r) and (1 (1-r)2) >
(1+ (1 - 21)- Therefore, max(( (4/3)) |l., 1/2 - { and
1/ 401 ",
3و go
(43)
[Ko(())
(1 <(=+<=---0(3333-<<<<~~=)
v11.
29

---

By combining (41), (42) and (43), we finally have
11
1/2 (D)) (o
S (1 (1- (1--)))))) y'll Do
+ + - ( 3 oo
4rr(1 + (1-r) (2 (1 + (1 - ())) - - y 1l:0
Conclusion of Step 1. By combining the results of Steps 1.1 to 1.4, it comes that (Koo(( [Ko()), Ko(v)] (w Ko (y/ Il [ui
= =1000 1/2 (a1 (v) - (1()))ll
- (2 + (1 (1-r)) (12r)))(1(((11))))) 5a(a 1)/3(1 -r)- 6)
8r(1 (1 (1 (1- 4r(1 -) (1+ + "r) v el lwolo
Step 2. Let now bound the second term in (34). For any x
Bo.
denote
Mo Dgu (DKo(y)lw) Ko)) by (xx). We have
0o (d2()) - a2(3/)) (90 122 (M) Mo (y/) { (3/)))/²' Ho. 1/2 (Dgo Ko()) D (KK)('))) Ko()))))
0 1o: 1/2 Dgo (w)]), Ko(D)]
(44)
We now aim to bound each the four terms that appear (44).
Step 2.1, First, we have
1/2 (Mo(y) Mo / /1 (g(v)- g(y') 00 2
9
s 3(1 - 31y y 1l.
(on the model of (18))
Step 2.2. Secondly, we have
9o
1/2 D (KK)()))))) Ko(vill2 < [2)
- x )l -
(Lemma 9-(a))
2 (1=(+++(11))))]]>(((<0))))
(see Step 1.2.)
where
[Doo( S <1l00
1/1 (B())s5' +1) 1/2 1/2 Da(y) -122|]
+
[M) (Mo(y)-\wllg""" g0"lo:.
Using (36) and (38), comes that
[|DKo())((]lle I (2r(1 -) r3 31 + (1- (
Then, have
1/2 Dgo DKo Ko(v)]] < 2r(1 (1 - r) {1 + 27(1 31 + - --))))- |1ollo0
Step 2.3. Thirdly, have 1100
1 M, (1))) 0
1/2 11 - +Lell2 1 (1 (1-r)-2.
Step 2.4, Fourthly, using Lemma 9-(a), we have
1/2
< || DD.. (Koo DK, DK. Ko B| + |1 (DKo(v)(w) Ko(v) - ko(r')lll
+ Dol
(o()) (o(s'))l 3.
< ()) DK,
21I (v) [ oo (K) Ko(D')ll.
- - DHI (K((y)
( - Ko('')ll.
< 2[ + | - DK, (1ka [w]]l]]
+ ( (o)
(K())ll
30

---

We recall that the following meequalities hold
(
Kooeoilli 0 5 -
- ylloo (1-) ",
|| (DKK(t)) 00
s <
1|
|| (v) - Ko g0 51
-
We are now going to bound || [DKo(y) (w) [ DKo
We have
[[DKov)(w] |DDKo - DKo(s') (w]ll) 50 < || (
)wolo + Ai + A2 + As As,
where
A1 [[ + 0o 9(y)00 1/2 (La + 90 g
1/2
- 1/2
you 1 ) 1
1/2
1/2 (1a+ + 9o 1/2 g0)) 80 1/2 11 90 (y x 112
1/2
x 9o
A2 = (I + 90
1/2 g(v ((v'))0
1 90
(Dg(v)) - Da(y' [1c] 90
1/2
1/2 (Ta + -00
1/2 1/2 s(y)00
(v - 1/0"
A3 = (I + 90
-1
1/2
g 0o 1/2 9o
Dg(y) go 1/2
x {Ia 90
1/2
g(y)g
1/3-1 -La 90 1/2 g 00
(110))/011 1/2)
300
( x x1 112
1/2
A3 = 1I + 0o In particular, we have
1/2 g(('(())) 0o ) 90 1/2
1/2
Dg(y)][[l 1/2 (La + 90
1/2
8(v')) (00')o
1/2 -
1/2 (y 1/ ll2
(iMo(''---MMoll )w/l
S 4(1 (1+
y oo l 10
((43)
40
11111111111111111111111111111111 x (4(1- (4(1-r)((++111)))))
(AI, Aa) S 8r(1 (1--)))) - y' al /olll/lllo go
((36), (38), (43))
A2 S 5a(a +1/3r(1-1) (1 r) 21le 2/ ol w0ll +
(37),, (38)
A4 s 2(1 -
s /l w0loo
(36) and (38))
ll.o ( / J40p0//0J666660o///////////////1
Therefore, we have
3,0111111))011++++1--11111111111 (f0--"""0001---33-----0"0""""""1 - s +
-= rr) (1 (1-r)3)-2) 11y - o/llool10/
+ +==
+ +[ + 201++22(1-)-((11 -"---"11==""00(0-711"
x (1==(=-))001-00000>0((0100-11111
- 11l1
Conclusion of Step By combining the results of Steps 2.1 to 2.4, it comes that
[[Mo(y) ( (DKo())l)) 1
|do goooo (v))[ Ko ll
= - (a2))) -2 (22)))l
+ 2r + 1,,,
-
(6r(1- 1,6600000"""00000111111111""""""
x (4(1-
+ + + + -1)-)((11000000----))>>--------- -3(1 + (11(11-))*) (1-r) -2
+ 2 + (1 (1-r) -1--"""111111---00------11--")11
x (1+ 01------03300-------0--00-011-1-
|y - 3 80 wol.
Conclusion. Finally, we have [D [
- DD0
90 S <c(r)]y v/lloo //lllo0 where
c(rr < 24r 21 (1-r)
+ 100 + + (100(011111111111111 -10
-
+ 32r(1 - + -1) 128 (1 - --)-- 17
-
- (1-2
+ +12r(1 - --) 24r2
 + (64r (1-)) 256r (1 --00011------------------------1
+ (48r(1 -r) 9 + + - - 2,h -))-
+ -00- + + (1- 11- 2- (1-
+ +(((1- + -- (1+
31

---

by combining the results of Step (two first lines) and Step 2 (following lines). Hence, using that
r
$ 1/11, we have
dr) < max (24 x 256, 5
80a(a + - -1)
s imax(6144, 400a(a 15/3
$ max 1000a(a + 1))
< 1/raa
(see (33))
i.., for any (v.3/) E Bo Bg and E RR, we have
[[ - Doo(y)(w]laa 1 (ra) (ra))111 fly - y Illloolwwlo..
and thus
[Doulv)~ - Dou(v) S I1y y 1l.0
This last inequality finally proves that Déo is (ra) Lipp--LLppchhtzzoonninuuuu on Bo with respect to || Bo
Inequality on Déo. Elaborating on the smoothness of Dén, we prove here that [D 1/2 for any y E Bg. Let y
Bo, have
|al go = (Doo(y) + 1
1 1 Déu Doo
1 1/2 Dibo(lll
We recall that (zo = and that DΦn is (ra) -al))----in L-Lipschitz-continnnn on Bn. Thus, we obtain [|Doo())llen (IDouly) ())( (ra) x 10 < (ra)
and then ( 1/2, which proves the result.
Smoothness of do and Φo- We prove here that do (and thus Φo) is Lipschitz-coontnuuus on Bo with respect || ll Let (y, 2/)
E Bo Bo. We have
((o(v)) (o(1 pa
(M))(v) -
Ko + Ma (v') (Do (Ko().. Ka(v)] Dga (Ko('), Ko (o(y')) Il
(45)
To prove the Lipschitz-continutty of Φo. we are going to proceed two steps, each them consisting of bounding from above a of (45).
Step I. Let first bound the first term in (45). Since r S 1/11, we have
[KMo(v) -
(v)
< / (Mo(2)) - Mo(r))) 2lll Ko(v)]1
< 3( - r) ->y o /l. 2 (b)
(see (18) and Lemma 9-(a))
x 61 - r) 3 (11)))))) Il.o
(see (39))
Step 2. Let now bound the second term in (45), We have
MM (Koo - (go(Ko(v)), (1)))))))
< 11g0
Kn(v)) Dgo Ko(y'), Ko(y')) I|
(1-11)1))) (K) p.' + D Ko(y) - Ko(v)), Ko(v) - ko((')ll
(1)+ (|(])))))) - ko(s
(o -1 + 2l](n(() - Ko(v')
<11+ (1)))-))(K(())) ko(2 (2) + [|oo()) Ko(r')llao)
32

---

where we recall that the following inequalities hold
11-33311--- 4r(1-4411110000 ((())))))))((--))))))))))))>))-) (1 ())) (19))"111 - ylloo
[o S r1 + (1 - -
[)o(()) - Ko(L/)llo - (1+
Thus, we have
Ko())) Dgo Ko 2 4r + (1 r1-
- 2(1 + 4r(4 r) (1 + (1 - -))))=1}3 y'lloo
Conclusion, Since r 2 1/111, we finally have
(1do(v) doly') 1oo S (6r2(1 - r) (1
4r + (1- - (1 (1 ())))))))))((()))))))))())))))) - y/llo0
4r 32r - -r) + 641 -r) -1) 2/ so
x (11)
4864/5rlly - y'll 3/1l
Therefore, do and do are respectively (4864/5r) and (1 + Linhhttous on Bo with respect to || Il
We are now ready to prove the second result of Lemma 22.
Invertibility of Ja Let x E Bn and w E Rdd We have (h/2)G = Mo(y)
and thus,
(h/2) D (x) = Deo(x)(u) = (o/ -
+ DMo(x))
((w]o((x)
Then, we have
21 (nC (1))
< (z)) 1/2
(46)
+ +1000 / M2 1/2
1/2
DMo [w] 1/2 | (x) s0
$ (1 + (1 (1- (h//2) DG 20 (a) (x) [mo]llo_ + |1 1/2 DMo(ar) ( - 1/2
"[w]0
12211220000 (1)10
where DMo =
Mo and
|l 1/2 (=) 1/ ||00 M, (x:) (z 1/2 12
-1 1/2 DMo(x) 1/1/2
(47)
< 2(1 r) + (1- (1-) 1wllo
(see (40))
< 2 - (--)))- 1
Moreover, using the Lipschitz- of @ on Bo and 20(2 = O, we have
[10o(=)]lo < [100(2) 20(xx + (z) |l.o (1 + (2432/55)l1 x (0) 1l < (1 + 2432/5r22
(48)
Since r S 1/11, we have + (1 r1) 5/2, and by combining (47) and(48), we obtain (100] < (5h/4) (0) (x)(w]]] + 4874(1 -r) -r)
Since (x) Bo > || / /2 and (1 - r) < 2 with r S 1/11, the last inequality becomes
(5h/4)][DG,...o (x)| .. 2 {1/2 9748r) } | ll
2 1/2 - - 2 1/4.
(see (33))
In particular, Jac[Gr....oo) (xx) is invertible, which concludes the proof.
33

---

We prove below that, for any 200) E T the set (z()(() is reduced to single point, for h small enough depending on 2 Lemma 23. For any w 2 0, any r E (0, 1) and any a 2 1, we define
(1 r) h1 h1 r, a) = min ( (1+ 2(2r + 3(r + w)
1
r + 1/2(r + w
(49)
We recall that the set-valued map G2 is defined in (29). Assume AI, A2. Then, for any r (0, 1/11].
any = (x.p) E T M and any h E (0, h1 (1plll(t))2 a)), there exists a unique z' E (2, C T*M such that E (z).
Proof. Assume AI. A2. Let E (0,1/11), = (z.p) T and h E (0.l r.))). where b is defined in (49). We recall that the map 8, is defined in (30). We define B Bii-1= (2,r) and the map B,,= T*M- + R" R by (x')
&, (z, z') for any z' T. Then, EA.2 (=') = z' if and only if 2' E G2 (2). The proof is divided in two steps.
Step 1. We first prove that B (B) c B. Let z'
B, we have by Lemma 9-(a)
E (z') =||= - (1((o())1++
Jp' |l(=)) + 11
)) + ( ''. "((())""""" 1
S / (211e' f11p'llo<=>>>>>>> 1 + ||(())
o(x) o(t>
1) + (())
where the following inequalities hold
(a) Ip')) Ilp' p pll((=) S + [|
(b) ((()))))---())))((()(>>>>))))>)) I 3|1 x'll x'llo(s) S 3r, on the model of (18), using
r < 1/11.
Therefore, we have
(1,,,(()) (z') - z|1e S h((r + (|ll 1) + 3r(r + (PPl)(()) 1) + (rr + llellot>) ))²)
< (((++ + lllll ((1 + 3r) + ( (1e()
< hv /r/ (11lll <r,
which proves the statement.
Step 2.
We now prove that en... is a contraction on B. This proaly notably recovers some elements of the proof of Proposition 13 (see Appendix E). Let (21, 22) E B x B with 21 = (#1 P1) and 22 = (22.P2). Remark that 12, E (2r,,, 1) (2,,, 1). Let us first bound [gg (21) (22) We have
(1) (21) 8 Bhh (z2) 1/2 [g(x)) (p1 p2) + (((1)) 9(x2) -pa)
(1)
'22}
(o(z)) (p) p2) + a0zz)=)(( (pI p2) + (o(x))- -01=-(x))2"00 (pp3
and then, since r S 1/11, we have on the model of (18)
(z1) - (21) (1) g (22))l - \/2 (1 +
1/2
|2 1ali1--Pllacce p2ll
+ | 1/2
((x2)) )0(z)) -(x))))))0000000))00)
||
< 5 [1 + (1 (1-1) r) H|pp p2ll -/ (1 r) - || x +
Let us now bound (2)
(21)
gi. (22) || We have
(2)
(2) (21)
_ (E2) ((D((z) 1g (z) pr. (()) 'pll -
34

---

which gives using Lemma 9-(a)
!l (21) (x1) 8 (=2)lltt))- (o))) pp) p(x)-
s h)]s(()) - po)l + ))((())))))))) "(p - p-)2p (C))
- h(r + || - - P2 hrllp: p2lla(0)-
Finally, it comes that
E (21) - Br,
S + (1 - P2l1a(2)
+ 3/ - ||0l|(()))))))))))))))) )| - #2llo(x)
< (1/2)<[| - x2
+ Ilp (Illat)) (1/2)][ 221l
which proves that Sn,s is contraction B.
Conclusion, We obtain the result of Lemma 23 by applying the fixed point theorem B,,= and
B.
Elaborating Lemma 23, prove in Lemma 24 that the only element of G2 (2(0)) verifies smoothness properties if h is chosen small enough, depending on 20))
Lemma 24. For any = (x, T any 10 and any 2 1, we define h(2,r, a) min (he
1,,, h2(2), hz(2,r))
(50)
where hi is defined in (49), and h2 and ha are defined by
1 + [1P1)((a)-1))
h2(2) 3/2 31plllaa)))))--
(+ -
ha(z,r)
We recull that rà is defined in (33). We also recull that the maps Goo ER Gn and Gh are respectively defined in (29), (30), (31) and (32). Assume A1, A2, Then, for any 200) E T*M and any h ((()()()), ra a)), there exists a unique element 2/ (2() B 11 (z(0), .ra). Moreover, have
(1/2) E T such that 2/ (1/2) E
(a) Jac. [8, (x(0), 2 (1/2)) and Jac[Ga (122) are invertible,
(b) Ho (2
(1x (122) is invertible.
Proof. Assume AI. A2. Let 20)) T and h E (0, (0, ((2(0),,,,,)) a)), where defined in (50). Since r I 1/11. Lemma 23 ensures the existence and uniqueness of 2/ (1/2) = (z) (1/2)
(2/2))
Phh
T'M such that =/ (1/2) G2 (2())) n Bi--1_oo, (=(0),rz) Moreover, we have Tn
(1/2) E (0 ra),
and thus (o. (1/2) is invertible by Lemma 22. Let us prove that Jac (8, (z)(0) =1 (1/2))
and Jac G2 (1/2)) are invertible.
Invertibility of Jacz (=(00, 22/22)). We first remark that Jacz (g, (=(0), (1/22) = Jac[vo. (122) where t'o.s = ld + ui and th is defined on T by
<0o(2) (-(0(x(0)))- + g(v)
1a(= (0))
We recall that || (1/2) 200
2000 |0000 rr= < 1/11 and thus define r1 1/11 - || (1/2) - -2(0) ||_<0)
and B, B Bl (2h r1). We set 7 1/11. Note that B C Bi1-1l-0s (2 ()) 1/11). We show
(1/2)
35

---

below that dy is Lipschitz continiz-- on B1 with respect to -12.00. Let (2, B1 x B1. Similarly to the proof of Lemma 23, we have
( -
s 11p 1--r'l 1'(1o(o10))-1
+
12 p' (()))
+
g -
$ (1 + (1 - f) -) Ip p'l1
+ 3( F) (F + |p00) |a 1 )|[ x' 11a
In the same manner, we have
110% - "/0 (=))) (2) S 2(22 + |p
| 11
Therefore, recalling that F = 1/11, we obtain with the previous inequalities
[1xo(2) - vo (=()l)) I {3 + 11 ()) H|| - =²11=000)
s 2| - ='1l 200 ()))
Hence, to is (iine=="")-iisscchio-mmtnmoonn on B1 with respect to || - ol Then, since < h2(2(0)), (ac(v),,x(2)) is invertible for any 2 B1 by Lemma 21. In particular,
Jac. Jack/|s,, (2 2h (1/21) Jaclvo.<l(=) (11/2) is invertible.
Invertibility of Jac Jac(Gn](=1)))) (1/2)) Let us remark that G.. = Id + hy where ( is defined on T M by
((z) = (0, 4 Dp(=))((x)"'p (a)))((=pp)
We define r2 = 1/2 and B2 B (//a (z) r2). We show below that Ç is Lipschitz
on B2 with respect to || ||_g/>) with Lipschitz constant which does not depend on 24 (1/2) Let (z, 2 E B2 x B2. Using Lemma 9-(a) with r2 = 1/2, it is clear that
ll<(=) ((2' _|g n/2)) (1 1p/ (1/2)
a( (/2)) 2' 111_0//))
Moreover, since || (1/2 - z 20)
|xx00 < r.. we have by Lemma 12
11r 1/2) 2/21)
S (1 ---r,)² ~1pp" (23) Bolaton
< -1--ra)-10 ( + |1 (0) (1)))))))))
and thus
[I((2) - ((22)((1) (/2) S [1+ (1 - : - (r + lle 10 |
< || <' (|ايا (2 r())
Hence, < is ((((()))))))))))
(2th)-Lipchitzzconnioooss on B2 with respect to |[ Since
ha(2(0)) Jac[ Jac[Gr]( (=) is invertible for any E B2 by Lemma 21. In particular, Jac G2](2 (1/2) is
invertible, which concludes the proof.
The following lemma states the existence of a diffeomorphism between (0) E T and z(1) (z((2()) under smoothness conditions verified by 2(1)
Lemma 25. Let (1)) E T M We recall that the maps EN Gb and G₁, are respectively defined in (30), (31) and (6). Assume that E CM, Raxd) and that there exist h > 0 and (1/2) E T with the following properties:
24
(a) g (2 ,2 (1/2)= = 0 and 2 (1) G2(221/22)
36

---

(b) Jac )x (/21) and Jac Jac[Gr](( (1/2)) are invertible.
Then, there exists a neighbourhood U C T of 20 and a C' Ci-dfffomorphiis &: U + E(U) c T such that (i) [(=(0)) = = for any = E U, {(=) E Ga(z) and (ii) det = 1.
Proof. Let (0)) E T x T.M. Assume that g E ((M,R4x0 and consider h > 0 and (1/2) T'M as described in Lemma 25. Under the assumption g. 8,, and G₁, are continuously 2
differentiable respectively on T*M T*M and M. We are going prove Lemma 25, by first deriving intermediary results on 5y and Gb.
2 (1/2)) is invertible. Then, by applying the implicit function theorem on 8, (z!0), -/ (1/2)), ) we obtain the existence of neighbourhood Uo C. M of (0)) and C df-diffoomorphiss to Up - Co(U)) C T*M such that (u(2(0)) 2/h (1/2) and for any
Result on Ex* We recall that SA (2 2/ (122) and Jac (g, (=
Uo. g, (z.Eo()) 0, i.e. Eo(z) (z) Moreover, for any E Uo. the Jacobian of &o is given by
Jac[Eu] (z) Jac. [g, (2,, Eo(2)) Jac, B (c,))()))
Result Gh. We recall that 2(1) Gb (2) 2 and Jac[Gn] (1/2) is invertible. Then, we apply the inverse function theorem on Gn, at =// (1/2) and obtain the existence of neighbourhood U₁ C T (1/2) such that & = Gulu, is C² -dfffcomorphism on U1.
of 2/
Final result.
We now consider the subset U and the map respectively defined by
(a) = 5 (U, n (o()) c T*M, neighbourhood of (0))
(b) E Eo. C diffeommmrphis U such that E(2(0)) 21) and for any E U,
{(z) E (z).
Let us now prove that d 1. Let z U. We define zo Eo(z). By the chain rule, we have
[ac(E](z) Jac([ı Eu](2) Jac(Ei] (=)
= Jac({i](Eo(=)) Jac. Ig, (2,, Jac.[g,11z,50l1
where we have
H2(2 (-100_ La H0nnHH((2') 0g
La
0
(a) Jac Jac())((')==
II -"a- (2,nn
- 10p H2(, Opg OppHa((''p')) ) Ia h a H2(x111)
(b) Jac,, [s, (2,
(c) Jac, 1,, (-l - (~
0
10. -14)
Therefore, obtain
[dctJac[[[(=)]
det - ha (o(o))
h Or. - (2(2))
(Ta h H2 det(Ta + h 2 Hz. =
= dde
2
2
which concludes the proof.
Lemma 26. Let is > O. recall that the set-valued map Fú, is defined in Section 3.1. Let z M.
Assume that there exist (z',²") T x T M such that
37

---

(a) 2' E (2)
(b) z" E Fn(2), and
(c) x'= = x²"
Then, we have p' p",, and thus 2' z"
Proof. Let h > O. Let E T We consider (z', 2") E Fn(2) x (z) such that x' = x" = B. By using the last two aligns from (6) with z' and =", we have p/ = -y(1/2) + h where p(1/2) (g(x)) =()
2 (0 1101)2)))11) p",
z),, which concludes the proof.
Under the assumption A2, we define, for any (0)) E T*M, (z("))) = ((s()(0)), ra, a), where a is given by A2, and ra and h are respectively defined in (33) and (50). Note that h² appears in Proposition 2, for which we derive the proof below.
Proof of Proposition 2. recall that maps Go. 8, and G2, are respectively defined in (29), (30) and (31), Assume A1, A2. Let z0)) E T*M. We define z(0) s(2()). Let E (0, (2())) By using Lemma 24 on 2(0), we obtain the existence of 2/ (1/2) E T with the following properties:
(1/2) (a) 2h
E G (3()), ie., Sn (3(0), (2/21) O,
(b) Jac K (z(1), ( and Jac[Gn](=]] are invertible,
(c) Jac ( (0) (112) is invertible.
We then define 2/ = (2 In particular, 2f E Gz ie., 10 E Fx(2(0) and Jac[Gr,, to) (xx (2) Jac[G,, [r
2 [ccoo][(]) (2/ (1/2)) is invertible. By combining the properties of (1/2) and
with Lemma 25, we also obtain the existence of a neighbourhood U' C T of g and a C diffcomorphism EA U' + (U')))
T such that (i) (x())) = (1) (ii) for any E U'. (z) E (z) and (iii) det Jac Enl = 1. Therefore, we can define the subset and the map ? of Proposition 2 by
(a) = =(U), neighbourhood of (0) in T*M,
(1) (b) 7h = En o 8, such that (i) (2() = 2/ (ii) for any z E U, (2(2) E Fx(2) and (iii) det Jac[nk] 1.
We now prove that U can be reduced to a smaller subset such that Yn (2) the only element of F (2) in 7n (U) for any E U. Motivated by Lemma 22, we first define, for any (a, d',x") ) E M x M, F (x') as the only element p E T such that (a', .p') E Fn (x, for any p' E T M. It clear that Fa)((') GGhz (x²), where G is defined in (32). We also define 7,, M x M T by
Zn (x = (x,, F h. (2')) = (ar, G2),)
which is continuously differentiable since g C2 (M,, Rd Besides this, we obtain by Lemma 22 that Jac (o, (x) - (2)) is invertible, and then Jac [Zn] (1) is invertible. Therefore, by applying the inverse function theorem on Zn (x(0), (²)) (()) ). it comes that Zμ is C' difiimoeomoop in
neighbourhood of (z) 2z We are now going to prove the result by contradiction. Assume for now that there is no neighbourhood U of 20)) such that 7/ the only element of Fx(2) in (U)) for any 2 E U. Then,
find sequence (2))ig which converges 2 such that for any i E N. there exist two different elements 2,, E (2.)) and 21,2 E (2,). Therefore, for any e N. 2,,1 j 2222 and by Lemma 26, have
F P (21.1) Frizz (II.2) Pi
38

---

and thus
Zu (21,, Z4 (2 22,2)
(51)
Moreover, by continuity of The the sequences keNN and (Zi,2)kEN converge to 2h and therefore, (1511) DieN and (2i,2) also converge to x(1), Combined with (51), this result of convergence is in contradiction with the fact that Z₁, a diffeomorphism can reduce U to smaller subset such that V/n (2) is the only element of (z) in you (U) for any 2 E U, which concludes the proof of Proposition 2.
neighbourhood of (x:(0) Thh
(1) Therefore,
H Modification of n-BHMC algorithm with step-size conditioning
In the rest of the paper, for any 20) E M,M,, we will denote by h, (2) the value of h, given by A3.
Beyond n-BHMC. A crucial part of the proof of reversibility in BHMC, as much for c-BHMC for n-BHMC, relies local symplectic properties of the integrator of the Hamiltonian dynamics. Although Algorithm can be implemented without any practical limitation, it is hard to state such properties for its numerical integrator 4, under A1, A2 and A3, given any value of h. Indeed, we know from A3 that 2 is local involution around 20) E T when h < h. (=()))); however, cannot ensure this result when h h, (z("))... To circumvent this issue, Theoretical n- n-BHMC (Tn-Ba BHMC), presented in Algorithm 3. In this modified version of n- BHMM, actually enforce condition h to be small enough. We now get into the details of this new algorithm and assume A3 for the rest of this section.
propose to study
Theoretical motivations.
Let h > O. We recall the definition of the set A,, introduced Section 4.2
A E h minseB (a.1) (s.1) (3))
clear that A₁ C. domp,; indeed, if (0) An. we have in particular that h h, (2) and therefore z(0) E dome, by A3. Let (0) Ah. By A3, we know that 4, is an involution on neighbourhood of 2(0); in particular, it comes that (4,, o (2)) of the "involution checking step" in Algorithm I is de facto satisfied. This naturally leads to replace the condition "2 E dome. by the more restrictive condition "2 An."
z00) Hence, the condition (b)
Algorithm 3: Tn-BHMC with Momentum Refreshment
Input: (Fo, T*M, (0, N E N, 0, 0, tb with domain domp. Ab.
Output: (Xn, Pn) Pn)ne[N]
for n N do
Step 1: Gn ~ No P. -T /BG,
Step 2:
solving a discretized version of ODE (3)
X,, P. < Xo-1,Pu Pu:
p!o) (
if ( (x!") p(u) A then
7 (7(0))))
if
z E An then X,, P (so (z
end
10
end
11
Step 3: An + (x H(X..
U 00
13
if U. < An then X,, P. c- Xt,PP
14
else X., P. < Xn-1, Pr
15
Step 4: Xa, Pa + s(xn, Pr)
16
step 5: G La), Pa + -I VBG,
17 end
39

---

Implementation of Tn-BHMC. In Algorithm 3, we highlight in yellow the modifications of Tn-BHMC contrast to n-BHMC. Namely, we replace
(a) "Z(00 dome,"" "(Line 5 in Algorithm 1) by "Zi E A₁," (Line 5 in Algorithm 3).
(b) "Z(") E domo,"" "Line 8 in Algorithm 1) by "h E An" (Line in Algorithm 3).
Note that there no need to maintain the "involution checking step" of Algorithm 1, since it is automatically verified once Z As. On the whole, these new conditions more restrictive than the conditions of Algorithm since As. C dome,, moreover, they can be thought as conditions directly applied on the step-size h. The specific choice of Ay. instead of another subset of dome, is actually sufficient to derive the proof of reversibility Tn-BHMC (see Section 4.2).
Proofs of Section 4.2
I.1 Expression and properties r*
Given Riemannian manifold (M, g). we define r (I) for any E M by
(x) min()l 1/2 (|l(x=)
(52)
Note that r* is used in A3 and that (xr) = 1/C2, where C2 is defined in Lemma 12. We prove below that r* has smooth behaviour on M under our main assumptions.
Lemma 27. Assume A1, A2. Also assume that M (x) ll2 is bounded from above. As defined in (52), r" M + (0, +oo) satisfies the following properties:
(a) (²)x) + 0 as x + 0M.
(b) There exists L > such that r* is L-Lipxchitz-coniinuuus on with respect to 1|2.
(c) There exists such that r* (x) < M.for any x E M.
Proof. Assume AI, A2. Also assume that x E M + |gg (x) 2 is bounded from above. We first define M + ||g(=)]l and r2
M ||((x)=11 || ||/ such that r" min(r1, r2).
Since g E C2 CM, dxd), it is clear that 71 and ra continuously differentiable on M. We have:
(a) (1) + 0 as + 0M by Lemma 8, and
(b) (2( A 0 as x - 0M. since 1/r3 æ + [g(x)) -1 ||2 is bounded on M.
Combining the fact that (2 > for any x E M and item (a), we obtain item (a) of Lemma 27. We denote by d the distance induced by and now define, for any >> 0,
(i) inf (yEM: d(). r2(U).
(ii M, Int((x M d(z, OM) (x) S pa)). open M.
(iii) M_ M MMM, closed and bounded (and thus, compact) set in M.
Using items (a) and (b), and (ii) Me and M_ are not empty. We consider such for the rest of the proof.
ensure the existence of some (0, diam(M)) such that (i) M 0
Smoothness of rz. We define d(Ms, M..) and M, M B(0,8/4). Note that
(i) 8 > 0 since M C M,
iii M3 is closed since the ball (0,55/4 is compact, and
iiii M2 (1 M_g 0.
40

---

According to the smooth version Urysohn's lemma applied to M2 and M-2, there exists E C' "IR", [0. 1)) such that x(M_c) = and x(M;) O. We then define ñ2 R2 + (0, +00) by '2 x + (1 x)Me.. In particular, (i) there exists L2 > 0 such that 12 is Lg-Lipschizz-conimmon on M with respect to 12. since F22 E C² (12d, [0, 1)), and (ii) for any x E Mg, F2 (a) > Me.
Smoothness of r1. We now prove that 71 is 1-Lipschitz continuous on M with respect to || 1l2. Let x E M. Note that Ti (x) = (Ilo(>)1l3) 1/4 and we thus have
Vri (2) = (- (-1/4)ll(=))1 5/ h(xx) Dg(x),
where h(z) 0 اا()) 21]g())) (2)()) u(x) being normal eigenvector of g(x))
corresponding to the eigenvalue [g(())]l2 Hence,
[Vrr(>))22 - 11
Dg(x)ll2
s $ (1/2)l1(()112 3/2 2|
2|10(>))1 1
(Definition 1-(c))
which proves the result on Ti.
Smoothness of r. Let E M. We can face two cases: either, x e M_ then r2 = f2(2); or. r E M., then (z(z) > He 2 (I) and r2(2) 2 Ti (2) by definition of M.. Thus, we have for any E M, (x) = min (ra (f]((((() (2()) where 71 and f22 are respectively and L2 Lipschitz on M with respect to
|| |l2. By observing that 2 min(r1, f2) = + F22 Ir Pzl. we set L = + L2 and thus obtain item (b) of Lemma 27. Finally, item (c) of Lemma 27 directly derives from item (b),
since M is bounded.
Note that the assumptonn on g used in Lemma 27 is not directly ensured by self-concordance but can be proved when M is a polytope, as shown below.
Lemma 28, Consider a polytope M defined by m constraints with > d such that M {a Ax < b}. where A E mRmxdd is a full-rank matrix and b E TR. We endow M with the Riemannian D²o(x) where Φ M + R is the logarithmic barrier given for any x M by <(x) -m -i= In (b. ATT)). In particular, M and g verify Al and A2. Then, the function r: M > (0, +00) defined by r(x) ||g (x))22, for any E M, is bounded from above.
metric g(x)
Proof. Consider such manifold and metric g. We aim to show that the smallest eigenvalue of g(x) bounded from below for any M by constant >> O, which does not depend on x, Le., for any E IRd and any I E M. g(x:) [h, h] 2 c//1
Since A is full-ranked, A A is positive-definite In particular, for any h E Rd, (ATA) A) > Amin (A |/ where Amin AA) > 0 is the smallest eigenvalue of ATA. We recall that have for any E M, g(x) = AT ((z))((A)))))) where S(z) Diag(b: AT (m) Let i E [m]. The function r; x E M + S()),,? has the following properties: (i) r, is continuous on M, (ii) rı(x) for any x E and (iii) r(z) + +00 as -> ÖM. Thus, there exists G; > such that for any x E M, ri(x) 2 q. We define ē = minge(m) G and we have for any E M, S(x) > T.. We now define
Jimin (A A) and have for any M
g(z) h 2 8- (AT A) M| 2 c|>b
In particular, g(x)-1 = (1/c) Le. [1((x)) 1l2 (1/c), which concludes the proof.
I.2 Markov kernels of Algorithm
Based on the model of É, defined in (24), we define the set Eg C TT,, which will ensure that the maps derived from implicit integrators of n-BHMC are properly expressed
Ex (s o Sp - n t (doms,)) )) = (s (ddm) n n
n
(doma,))) (
(53)
41

---

For any (2,2') E Ex x T,, we define the acceptance probability ā(21|z) 2) to move from 2 to 2' by a(=1l=) a(z' (t) 03) a( 2) (R
where a(=' =) is the acceptance probability defined (9). We denote by Qn : T x B(T*M) +
[0, 1), the transition kernel of the (homogeneous) Markov chain (Fm:Po)nE[NY generated by Algo- rithm We also denote by:
(a) Qa T"M x M(T)) + [0, 11, the transition kernel referring to Step 1 (also Step 5) in Algorithm 1, defined in (25).
(b) QR, T M() + [0, 1], the transition kernel referring Step 2-3-4 in Algorithm 1.
We provide below details on Markov kernels Qn and Qn,1-
Kernel 1z This kernel is deterministic and corresponds to the numerical integration of the Hamiltonian up until time h. For any (2,2') E T TMM, we have
2n,1 = (m,)), de, + + (Et)
(54)
where
Q2,2(2, (2') a(R (2) 2) (dz') + [1 a(R 12)]] 22)]) (55)
Kernel Qn. This kernel corresponds to one step of Algorithm 1 (i.e., comprising Steps 1 to 5), For any (z, 2') E T T*M, we have
(n(2,dz') ST-MXT-M Qo(2, ) (= (z2 dz)).
1.3 Proof of reversibility in Algorithm
Let O. Using notation from A3, we recall definition of the set A₁ introduced in Section 4.2
A f= E T*M h min (=.1) h, (2)) domp,
(56)
We also recall that म, as defined in (2), admits density with respect to the product Lebesgue measure given for any 2 = (x, T M by
(dz/(dzdp)) p) (1/Z) 1/2 det(((())))))) exp]-V
Since Algorithm 3 can be thought Tn-BHMC are similar to the kernels from n-BHMC, defined in Appendix I.2. particular, the transition kernel corresponding to the Gaussian momentum update (Step and 5 in Algorithm 3, Step I and 5 in Algorithm 1) is the same and is reversible (up to momentum reversal) with respect to т (see Lemma 17). Namely, replace
restrictive version of Algorithm 1, the Markov kernels from
(a) the set Eg, detined in (53), by Eg. (b) the kernels Qz,1 and Qn, respectively defined in (54) and (55), by Q1 and Q2.
where
Ex = (so ( (A, n 2, 14, = (s S,/2) n (An))
(57)
Q1(z,dz') ياا (2)
dz') 1 ( (2)
(58)
Q2(2, aRR (2) 2) 2) (dz') + + a (z) =)]6
(59)
We denote by Q M M(T*M) + [0, 1). the transition kernel of the (homogeneous) Markov chain (In, Pm) NEIN generated by Algorithm 3. For any (2, z') E T x T*M, we have
Q(2. d2') -T--xTTMM Qu(2, dz dz) (#1, d22) Qu dz')
(60)
We first turn to the reversibility up to momentum reversal of Qı with respect to म. We start with the following lemma which is key establish this result.
42

---

Lemma 29. Assume AI, A2, A3. Then, for any compact set G C T for any C(TM x
T*M.R R) such that supp(g) C G. we have
JG. 9(2, (x(2))d= Jc.. (2),
where GA = Fn nnn (G), Fa An 2 (An), An being defined in (56).
Proof. Assume AI, A2, A3. Let G C T be a compact set, g E C(T*M x T such that supp(g) G x G. We first define the sets F = As n 2n (An). G. = FxnGn o (G) and the integrals l' and J'
I' -==G, s(2, (z))ddz J' J' g( ((), z)dz.
(61)
We recall the existence under Al and A2 of >> 0 and M > 0 such that r*, given by A3 and defined in (52), is L-Lischitz-coniiuuuus on M with respect to || ||2 and bounded from above by M (see
Lemma 27).
We define r T*M (0, +o0) by r(z.p) r (a)/m /m for any (z.p) E T*M, where m max (1/X, 4M, J and A is given by A3. Using the properties of r*, comes that
(a) is L-Lipschitz on T'M with respect to ||
1|2.
(b) r2, S 1/(LCC,, for any (z.p) E T where C2 is defined in Lemma 12,
(c) - 1/4, and
(d) S Art.
Note that G UzeGB rz)). Since G is compact set, there exist K N and (Zi)iE[K) T*M* such that G c U B,, where B, BiHll, (2;)).
We consider the sequence (Vi) constructed as follows: V1 = Gn and for any i E K Then, we have that, for any i E {1,..., NJ. Uj=VV if i, * iz. Therefore, we get that LK V; and $ 1 LK k (V:). In particular, for any A E B(T and any Ç C. R)
V, (Gn (u)
(n B,).
and that for any i, ig (1....,KK) K}, Vi, n Vi,
Sana ((=)(= E Su S(z)d2,
(62)
and for any À E M(T*)) and any & C,
(ann č(z)d== = Eie Sor ((z)d=.
(63)
Using (62) and (61), obtain 1' ΣK 1/ where 7i Iu.w (G)),, (x))dd for any i E [K]. We are now going to show that, for any i [K]
I = J (V.)GGFF
(64)
Let i E [K]. We proceed by making the following case disjunction:
(a) Either V; n (G) F. 0. and then 1/ = 0. We prove by contradiction in this case that nnnF, 0. Assume that there exists E 4 GV. By definition of F A 4 (z) E An. and we get that (2,) = 2 using A3. Hence, we have:
E (x)))(2)) 2 G. Un(z) E An.
(2)( An.
Therefore, we get (2) E V, n 2 (G) = 0. which absurd. Finally, (64) holds since
r = 0 = Sog !(V.)nGF,
43

---

(b) Either, there exists some Σ E V2 n & (G) such that Σ E Fh- In particular, Σ E B,, and thus || 2ella. r(2,) 1. this case, by combining A3 with Lemma 9-(c), we have for any 2' T
|| 2 $ (1 || - = < (1 - -===))
By considering 2' = 2, - 2, it comes that
|| =|| < (1 r(z)) -| - zilla < < - (1-r(z)))~r((
By combining properties (a) and (b) with Lemma 12, we have the following upper bound of r(2)) r(z) S r(3) + LI 21112 r(3) LC2, 2:112, < =(()
and thus
c - =1 (1 r(3) - 1 1/4) <1
using property (c) r. Hence, 2 E BI (2, 1). Moreover, Fa. C Ans and then it comes that (z.). Therefore, by combining property (d) r with A3, we have
(i) V2 B, C BIl (z,, (z)) doma..
(ii) the restriction of 4, to V, is de-diiffeoorphhss such that 4, 4/ Id.
This last result provides proper assumptions on 4, operate a change of variable in 1/, We now ( (F)) and G2,4 2 "(V)) G Fb. and prove that G1 G2,4 define G1 4 n (G) in two steps.
ii We first prove that GI C G2,4- Let = E G1.b- Then, there exists z' E V; (G) such that (2(2).. Since 4, is involution on V,, we have (2) 2'. Since 2' E V, V,nA,, comes that z b (V+) & (Aa). Moreover, z' 4% (G) (An). and thus, have 2 E G2 which proves this first result.
G Ai. Then, we
(ii) We now prove that G2.A C GL Let E G2 In particular, E G. Then, there exists j E [K] such that z E Vz. Therefore, E B, and thus || 2y11s, < <r(zg) 1. Since 2 E A those written above. In particular, this proves with A.3 that 4 is an involution on Vj. Then, 2 = 4 (*)),, with (=) E V, n A, and (z) E 2% (G) n 1 (An), since = E G A,. Therefore, we have E GI.,. which proves this last result.
obtain that h (zz) with the same computations
Given the fact that tb is involution on Vs, operate the change of variable Fn(2) in 7/ and obtain
= - 9 ()), So g(9)(2);
which gives (64).
Therefore, combining (62), (63) and (64), get
1= ΣK T (C)a, g(2. (z))ddd EK So (V)GGEF a(0n(2), = J which concludes the proof.
We are now ready to establish the reversibility up to momentum reversal of Q1, defined in (58), with respect to #.
Lemma 30. Assume A1, A2, A3. Then, for any > the Markov kernel Q1 with step-size defined in (58), is reversible up momentum reversal with respect to f.
Proof. Assume A1, A2, A3. Let h > We recall that the kernels Q1 and Q2 are respectively defined in (58) and (59). We define the transition kernel Q3 : T*M M x B(T M) + [0, 1] by
Q2 dz') leş
dz') + ( (d))
(65)
such that Q: = # Q3. The rest of the proof is divided into two parts. First, we prove that Q3 is reversible with respect to т. Then, we prove that Q1 is reversible up to momentum reversal with respect to म.
44

---

(a) Let f E C(T"M x T*M, R) with compact support. We consider compact set K with respect the topology induced by the set [B (2. T*M,r (0, 1)} such that supp(f) K x K. According to Definition 4, we aim to show that
ST-XTT-M
JT-MXT-M (z))) (d)))
(66)
We denote by I the left integral of (66). By combining (59) and (65), have I 11 12 13 where
11 Jep (())a(R))) (2)
T2 JE: = a (R (=) =))
1a jēgr (z)
Since supp( C K x K, we have for any = E ke (2.-) = 0 and f(:22) = 0. Note also that for
any (K)) R & K, and thus (RR) = 0 and f(; (=)) 0. By combining
these preliminary results with (9), we get
1 = Jeg Jegnknak (2))a(R((3) (+) =)) (=))dz
= (1/2), (1/Z) (ee.eeennnnn)))- minfe "(*), - (+)) Hf2,, (z))ddz
1 = Jegnk # (x) (2)
Is = Jegy- (())f(zz)ddz
We denote by J the right integral of (66). By symmetry, we have J = J1 + J2 + Ja where
J1 (Egr ((K) ((z)a(R (z) (R) z)dz
(1/2) Jegnknaee 2(K) min(e H(x) - ((
J2 JEgrK 712)[1 (=) 2))
Ja Seggr- (z,))(=
We directly have T2 J2 and I3 J3. Let us now prove that J1-
We recall that s o Sn/2 is a symplectic C'ddiffeomorphism (see Section 3.1) and Ex = (s (b, where Fr. , ). using (57). We define K₁, Fsn(soS,, (K))) ( (()))))/2)())), such = (1,)), and we operate the change of variable + (a0S),22(=) that E (R
in 11 and J1
I = (1/Z) Jk, min(e H1( ()((((((\\//)())))) H1 M((000 (f((s) (apposeeee (so Sa/2 o (n)(2))dz
J1 = (1/Z) Jk. min(e- (((ff=-)((SS,,)(())) H(( o0 [[(( o 41) (z), (z)
We now define the map g: T*M x T R and the set C T by
g(z, z') minfe ((s)) #2(so8)/2)==))) } (((((((((]))) Sn,z2 Sw)2)(=)))
G = (s o (K)
Note that G is compact set of T by continuity of Sp/2 and g continuous function by continuity of H and 8 Sa/2- Then, we obtain Jı = Jr, , by applying Lemma 29 with and G. Finally, obtain T J and thus prove (66) for any continuous function f with compact support.
(b) Let f T T + R be continuous function with compact support. We have fT-MTTMM 22,
-S-MMTTMM fls,=")
- (t
(momentum reversal on z')
(
(using (66))
= f
(momentum reversal on 2')
which concludes the proof.
45

---

We are now ready to prove Theorem 3, which states that Q is reversible up momentum reversal with respect to т.
Proof of Theorem 3. Assume A1, A2, A3. This proof is very similar to the proof of Theorem 15 (see Appendix E.2). Let f: TMM x T*M - R be continuous function with compact support. We have
M f(z, 1rw.:h//hha7r000oocccchh1aa1aaaa
(see (60))
(d2)
(Lemma 17)
(utp (dz)
(momentum reversal on 2)
-STT-My [tterh/eheeaaano///...oo.oooo:o0 (Lemma 30)
(momentum reversal on 2)
(Lemma 17)
Moreover, S Hence, by combining Definition and Lemma we obtain the result of Theorem 3.
L
J Comparison with Kook et al. (2022a).
In this section, we provide precise comparison between our work and the results stated in Kook et al. (2022a). We first dwell on the differences related to the general setting and the algorithms, and then explain how methodology may solve the limitations of Kook et al. (2022a).
J.1 General framework
Given convex body K Rd and function Bd - R" consider
= {a E K c(x) = 0}
(67)
In their paper, Kook et al. (2022a) aim at sampling from target distribution with density given for M by
(x)//xx= =xx((-
where (M, z SM ex)]-V(=))dz. They make the following assumptions:
(a) K is provided with a self-concordant barrier ø.
(b) The differential of x, De(x:,, is full-ranked for any M.
Although this setting might be general than ours (consider for instance the where K is polytope and e is a non-trivial function), Kook et al. (2022a) focus on the case
(68) where é, a e ddd with << u, m and A 1dxm Combined with assumptions (a) and (b), A must have independent rows and the Hessian of ø at x given by g(z) = D²x(z) is diagonal positive matrix for any x. Moreover, defined by (67)-(68) is a polytope. However, not any polytope can be rewritten in a manner akin to (68), and therefore, the setting chosen by Kook et al. (2022a) can actually be considered a special case of our setting, see Section 2.
k= (a Rd E S u), c(x) Ar
Then, they consider the extended state-space M = ((x,p) p) E IRd x Rd :ad:: : M. E Ker(Dc(z)) and define M constrained Hamiltonian H. given by
H(x,p) = H(x,)) p) + X(x, p) c(x))
with (xxp)) = V(x) + log pdet M(x) + 7 M(x)1,,
and M(x) = Q(()))))),
46

---

where Q is the orthogonal projection into Ker(De(z)), implying that M is and Mx, p) is Lagrangian term given by
positive,
p) = (Dod
dr/dt] (x.)))
Remark that simplifies when c is chosen as in (68), but may not be well suited for higher order constraints. Then, Kook et al. (2022a) provide simplifications of И for the setting (68) based on formulas for pdet M and M
Finally, they consider the joint distribution म given for any (z,p) E M by
dm (1/Z) (x, p)]dadp,
where ž = JM expl-H (2, drdp, for which the first marginal is æ. They naturally propose to sample from by implementing a version of RMHMC relying the constrained Hamiltonian H, which they call CRHMC.
Similarly to our approach, they consider two cases.
(a) First, Kook et al. (2022a) assume that the Hamiltonian dynamics given by H can be explicitly computed in continuous time (which is not the case practice). They present the continuous version of RMHMC which incorporates the simplified constrained Hamiltonian H, see Algorithm 1. We refer to this algorithm as C-CRHMC. This algorithm is identical to our algorithm c-BHMC provided in Appendix D. apart from the involution checking step (see Line 4 in Algorithm 2). They state in Theorem 6 that the Markov kernel corresponding to one iteration of c- C-CRHMC satisfies detailed balance with respect to म. thus ensuring that it preserves म.
(b) Then, Kook et al. (2022a) provide practical implementation of CRHMC, which is similar to n-BHMC, based on dimediiscrtiization of the Hamiltonian dynamics for a given time-step h > 0. They first split the simplified constrained Hamiltonian H into Hj and H2, by leveraging the non-separable aspect of H in H2, as we do in Section 3.1. Then, they propose to use the first-order Euler method to solve the discretized ODE associated to H1, as we also suggest, and the Implicit Midpoint Integrator (IMI) to solve the discretized ODE associated with H2, while we implement the Generalized Leupfrog Integrutor (GLI). We insist on the fact that these two second-order methods share the same theoretical properties (Hairer al., 2006). Since IMI is implicit, they propose to approximate it with a numerical integrator t,,, which is computed with a fixed-point method (as we do) detailed in Algorithm 3. Finally, they implement CRHMC which contains the same steps as n- n-BHMC apart from the involution checking step (see Line 8 in Algorithm I): (i) refreshing the momentum with a symplectic scheme (Step I in both algorithms), (ii) solving the discretized ODE associated to H in three steps according to the splitting given by Hi and Ifz (Step 2 in both algorithms), (iii) applying a Metropolis-Hastinng final momentum reversal (Step 4 in n-BHMC). They state in Theorem 8 that the Markov kernel corresponding to one iteration of CRHMC satisfies detailed balance with respect to म. thus ensuring that it preserves 7.
filter (Step 3 in both algorithms), and (iv) operating
J.2 Theoretical gaps in the reversibility of CRHMC Kook et al. (2022a)
First, we call into question the proof of the reversibility of c-CRHMC stated in (Kook et al., 2022a, Theorem 6). Indeed, Kook et al. (2022a) implicitly assume that the time-continuous Hamiltonian dynamics have a solution any time, but do not provide any proof of this result. We emphasize here that this result is not trivial and possibly may not hold due to the pathological behaviour of the barrier at the border of M, see Proposition 14. In contrast, the involution checking step implemented in c- C-BHMC verifies this condition for the forward and the reversed dynamics after momentum reversal, which then allows us to derive the reversibility of C-BMMC, see Theorem 15.
Finally, we report some gaps in the proof of reversibility of CRHMC provided in (Kook et al., 2022a, Theorem 8). Indeed, while Kook al. (2022a) claim that (Hairer et al., 2006, Theorem VI.3.5.) applies to CRHMC which contains the numerical integrator detailed in (Kook et al., 2022a, Algorithm 3), in fact (Hairer et al., 2006, Theorem VI.3.5.) only applies to the implicit integrator (IMI). This is a fundamental limitation of their theory. This theoretical confusion thus leads to numerical errors, as we detail in Section 6. contrast, we present Section theoretical results on the implicit integrator (GLI), which allows us to derive reasonable assumptions on the numerical integrator &,, to obtain reversibility of n-BHMC Theorem 3. Contrary to Kook et al. (2022a), our proof relies on the properties of 4, and highlights the critical role of the involution checking step.
47

---

K Experimental details
The numerical experiments presented in Section 6 based on the MATLAB implementation of CRHMC provided by Kook et al. (2022a). We adapt this code for n-BHMC for sake of faimesss. In partieular, our implementation differs from CRHMC by the use of the Generalized Leapfrog integrator and the "involution checking" mechanism.
Details on experiments with synthetic data. In these experiments, the hypereube refers the set (-1/2, 1/2/d where d E {5, 10). We recall that the ground truth on the quantities we estimate is given by the Metropolis Adjusted Langevin Algorithm (MALA) (Roberts Stramer, 2002) for the hypercube and the Independent Metropolis Hastings (IMH) sampler (Liu, 1996) for the simplex. We now give details on how the parameters these two algorithms are chosen. For MALA, we use constant step-size h = 0.05 for both dimensions, which results in average acceptance probability and 0.44 for 10. For IMH, we use as proposal the uniform distribution, which is simply Dirichlet distribution, resulting in average acceptance probability of order 0.36. We recall that we use an adaptive step-size h in CRHMC and n-BHMC such that we obtain an average acceptance probability of order 0.5 in the MH filter. We now discuss the setting of the tolerance parameter r) in n-BHMC. We choose so that (i) the step-size is at least greater than 10 -2 over the iterations of the algorithm and (ii) the average acceptance probability is roundly equal to 0.5. This heuristic results in defining (a) for the hypercube, 1) = 5
equal to 0.55 for
and 10 10 and (b) for
the simplex, 1 10 if d and 200 if d 10.
Details on experiments with real-world data. For these experiments, we consider the exact same setting in Kook et al. (2022a). In particular, step-size and the barrier functions.
the same adaptation strategies to update the
Influence of the norm in (8). The norm chosen for the "involution checking step" is arbitrary and many options are available. In particular, could design || Dy (2("))]22 17, and thus pick the Euclidean norm (8). However, while this approach should theoretically also reduce the bias in the method, remark that the obtained Markov chains have very poor mixing time. This is due to
fact that in practice the update on the momentum is of order [l9(>))l2 while the update on the position is of order [I((z)) 11 Hence the Euclidean norm is ill-suited for the "involution checking step" and the proposal states are rejected a lot more near the boundaries, see Figure 3.
0.5
0.5
0.0
0.0
0.5
0.5
-1.0
-1.0
Figure 3: Outputs of n BMB after 25k iterations to sample from the uniform distribution over [-1, 1]2 with 7/ = 10-3, 0.8, using in (8) the "self-concordant" norm (left) or the Euclidean norm (right). Red samples are rejected and blue samples are accepted.
48


Reviews:
[{"review_id": "LaLITqmIzK", "reviewer": "Reviewer_nmMp", "summary": "The authors resolve an existing bias with naive manifold versions of HMC by adding an involution checking step, and therefore leading the new HMC algorithm to be unbiased. \n\nI'm not super familiar with this problem, and I would like the authors to explain some of the technical details, as the proofs are quite involved for a conference review. ", "questions": "While I have many questions, I am still open to raising my score once these questions are adequately addressed. \n\n1. When you say there's no guarantee that $G_h$ has a unique solution, do you mean that (1) there are cases with multiple solutions, or (2) there are cases with no solutions, or both? If possible, can you also provide an intuitive example where this failure occurs? \n\n2. Can the authors explain to me the main intuitions behind why the naive implementations of Kook et al. (2022a) leads to a biased stationary measure? In particular, why does the Metropolis accept-reject step not fix this bias? \n\n3. To add to the previous question, on an intuitive level, how does this involution check fix this issue? \n\n4. It appears that this is not a problem when we can simulate the ideal Hamiltonian dynamics, is this related to the fact that the domain of flow is the entire cotangent bundle $T^* M$? \n\n5. The proofs are quite involved. Can you explain to me on a high level (a) how does the proof work, and (b) where do the complicated calculations mostly arise from? ", "limitations": "N/A ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors proposes a solution to an existing bias in manifold HMC, both justifying theoretically and proposing convincing empirical simulations. ", "weaknesses": "The technical details are quite involved and I believe could be explained better. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comment and thoughtful questions.\n\n### Answer to Question 1.\n\nWe refer to the general response for a discussion on the number of solutions of the implicit mapping $G_h$. We provide a one-dimensional example, where $G_h$ may admit 0,1 or 2 solutions for the implicit update of the momentum, depending on the initial condition and $h$.\n\n### Answer to Questions 2 & 3.\n\nThe main reason for the bias in the implementation of [1] is that **there is no guarantee that the numerical integrator used to solve the implicit scheme (Alg. 3 in [1]) is involutive**. This is however essential to derive theoretical results relying on the Metropolis Hastings (MH) filter in the Hamiltonian setting. Indeed, to obtain the reversibility of the Markov kernel corresponding to the Hamiltonian integration step followed by the MH step (see Lemma 30 in our paper), we need to show that for any compactly supported function $f$ \n\n$ \\int \\bar{\\pi}(z)   a(R^{\\Phi}_h(z) \\mid z)f(z,R^{\\Phi}_h(z)) \\mathrm{d} z= \\int \\bar{\\pi}(z)   a(R^{\\Phi}_h(z) \\mid z)f(R^{\\Phi}_h(z),z) \\mathrm{d} z$  ,  (a)\n\nwhere $\\bar \\pi$ is the extended target distribution, $a(z' \\mid z)=\\min(1, \\exp(-H(z')+H(z))$ is the MH acceptance filter and \n$R^{\\Phi}_h$ is the numerical integrator of the Hamiltonian $H$ with a step-size $h$. Since $R^{\\Phi}_h= T_h \\circ \\Phi_h \\circ T_h$, where $T_h$  is an involutive function, and $\\Phi_h$ is the numerical approximation of the implicit scheme $G_h \\circ s$ (see Section 3.1 in our paper), showing (a) amounts to prove that\n\n$ \\int   g(z,R^{\\Phi}_h(z)) \\mathrm{d} z = \\int g(R^{\\Phi}_h(z),z) \\mathrm{d} z$ , (b)\n\nfor a certain function $g$. **A sufficient condition to obtain (b) is thus to ensure that $\\Phi_h$ is an involution**. However, this condition is not satisfied in the CRHMC [1], and **the application of the MH filter does not solve this issue**.\n\nIn our setting, we directly enforce this condition through the involution checking step (up to some numerical relaxation). To support our point, we would like to mention the paper [2], where the authors also enforce an involution condition in a slightly different fashion than ours. Namely, given a current state $z$, they compute the whole set of solutions $\\Phi_h(z)$, choose $z’\\in \\Phi_h(z)$ according to some probability distribution defined on this set and then verify that $z \\in \\Phi_h(z’)$ (see line 7, Alg. 2 in [2]).\n\n### Answer to Question 4. \n\nIn the case where the Hamiltonian flow is well defined for all times, the domain of the flow is indeed the entire cotangent bundle but this is not sufficient to ensure the reversibility of HMC. The reversibility actually comes from the fact that the flow is involutive (up to momentum reversal) by definition of the Hamiltonian. \n\nHowever, in our manifold setting, there is one subtlety that needs to be checked, namely that the Hamiltonian flow does not leave the cotangent bundle in finite time. As detailed in Proposition 14, there is indeed no guarantee that this flow is defined for all times, due to the ill-conditioned behavior of the dynamics near the boundary of the manifold. Therefore, given a time horizon $h>0$, we restrict the cotangent bundle to the points where the flow is defined (see the definition of $E_h$ and $\\bar{E}_h$ in Appendix D) on the time interval $[0,h]$. Then, we prove that **the Hamiltonian flow is involutive up to momentum reversal on this subset (see Lemma 16)**, which is exactly the property needed to derive the reversibility of the algorithm with the ideal dynamics. Hence, we do not need to perform an *involution check* in this case. As a result, we introduce c-BHMC in Appendix D, where **the only checking step verifies that the flow is well defined on $[0,h]$, starting at the current state**.\n\n### Answer to Question 5.\n\nOur main theoretical results are Proposition 2 and Theorem 3. Proposition 2 justifies Assumption A3 made on the numerical integrator. Theorem 3 proves the reversibility of Algorithm 3. Below, we give more details on these results and provide a sketch of their proof.\n\n**Proposition 2**. In this proposition, we prove that in the self-concordant setting, for any starting point $z^0$, the implicit scheme $F_h=G_h \\circ s$ is locally an involution as long as $h$ is small enough. To prove this result, we proceed as follows. We first prove that we can build a symplectic diffeomorphism $\\gamma_h$ between $z^0$ and a solution $z^1 \\in F_h(z^0)$ on a neighborhood of $z^0$ (see Lemma 23, 24-(a) and 25). This proves the first item of Proposition 2. Next, using Lemma 22, 24-(b) and 26, we prove that $F_h$ coincides with $\\gamma_h$ on a given neighborhood of $z^0$. The main difficulty when establishing Proposition 2 comes from Lemma 22, which identifies the momentum $p^0$ used to go from $x^0$ to $x^1$ and shows that it is locally invertible as a perturbation of the identity. Proving this result requires computing higher order derivatives of the Riemannian metric and controlling their smoothness.\n\n**Theorem 3**. In this theorem, we prove the reversibility up to momentum reversal of the whole Markov kernel of Algorithm 3 (including momentum sampling, ODE integration, MH filter). The critical step of this procedure is to ensure that $\\Phi_h$ preserves the volume form on some compact sets under Assumption A3, derived from Proposition 2. This is proved in Lemma 29. To do so, we have to restrict the cotangent bundle to the set of points $A_h$, where h is small enough so that $\\Phi_h$ acts like an involution. Then, the proof consists in operating a change of variable on a partition of open subsets of a given compact set. The main difficulty in this proof is the control the compact sets on the domain of definition of $\\Phi_h$ in order to perform the change of variable in Lemma 29.\n\n\n[1] Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Kook et al., 2022a\n\n[2] Multiple projection Markov Chain Monte Carlo algorithms on submanifolds. Lelièvre et al. 2022.\n"}, {"review_id": "BtJutEaRCi", "reviewer": "Reviewer_prii", "summary": "This paper introduces a \"involution checking step\" in Riemannian Hamiltonian Monte Carlo methods that ensures time reversibility of the resulting Markov chain. The numerical experiments demonstrate the competitiveness of the proposed algorithm with regard to the method proposed in Kook et al. (2022a). \n\n---\n\nThe rebuttal addressed my quetions. I have raised the rating to 6.", "questions": "1. How does this work compare with Kook et al. (2022b)? \n1. Please address the gap between theory and practice in the weaknesses block. \n1. What do the highlights in the submitted files mean?", "limitations": "This is a paper on theoretically rigorous algorithm design. The assumptions are explicitly stated.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The algorithms are novel. \n2. The algorithms and analyses take practical numerical integrators into consideration. \n3. Appendix J points out clearly the theoretical gaps in Kook (2022a).", "weaknesses": "1. **Comparison with existing work.** The most relevant previous work seems to be Kook et al. (2022a) and Kook et al. (2022b). The second one appeared later and focused on theoretical analyses, but this submission only discusses about the gaps in the analysis of the first one. \n2. **Gap between theory and practice.** The algorithm that actually ensures time reversibility is Algorithm 3, but the numerical experiments are conducted with Algorithm 1. \n3. **Imprecise statement about existing work.** In Ln. 304, it is claimed that the convex constraint set considered by Kook et al. (2022a) takes a specific form. But that is simply a special case of the class of constraint sets in Kook et al. (2022a). \n4. **Presentation.** The presentation involves many terminologies without explanations. The paper may not be readable to a layperson. \n    1. There are some highlights with unclear meanings. \n    2. What do the vertical axes mean in Figure 2? \n5. **Typos.**\n    1. Ln. 115: Concordance -> self-concordace\n    2. Ln. 118: polynomial -> polynomial-time\n    3. What does the colon mean in (4)?\n    4. Ln. 287: a involution -> an involution\n    5. Some words in the reference list are not properly capitalized.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to review our submission and for your constructive feedback.\n\n> 1. How does this work compare with Kook et al. (2022b)?\n\nWe refer to the general response for a detailed comparison with [1]. In particular, we show that [1] suffers from the same shortcomings as [2] when it comes to the assumption made on the implicit integrator. We will add this detail in the section of the related work.\n\n> 2. Please address the gap between theory and practice in the weaknesses block.\n\nWe are aware of this slight discrepancy between theory and practice, and will highlight this in the Discussion section as a limitation of our study. We emphasize that **directly analyzing Algorithm 1 as such is a very challenging task**. Indeed, the self-concordance properties are only defined locally, whereas deriving reversibility results require global control. More precisely, given a certain current state, we have to **ensure that $h$ is small enough** (depending on this current state) to obtain that the implicit integrator (or its numerical approximation) is locally an involution (see Proposition 2 and Assumption 3). However, in practice, $h$ is fixed by the user. Thus guaranteeing that Proposition 2 holds for any element of the constrained set and a fixed step size $h$ is difficult.\n\nTo overcome this limitation, we study a modified version of Algorithm 1, given by Algorithm 3, where the involution checking step (ICS) is replaced by a stronger checking step involving the set $A_h$, *which implies ICS* (see Section 4.2) and for which we are able to state its reversibility. Although this algorithm may be implemented, this would come with a huge (and unrealistic) computational cost since verifying that $z \\in A_h$ implies to find the minimal critical step-size $h_\\star$ on a neighborhood of $z$. The question of the existence of a global step-size $h$ such that Algorithm 1 can be properly analyzed is not the topic of this paper and is left for future work.\n\n> 3. Imprecise statement about existing work\n\nIndeed, the setting presented in the introduction of [2] is a convex body $K$ provided with a self concordant barrier $\\phi$, combined with a linear equality constraint $Ax=b$ (equation 1.1). We remark that this setting is still a **special case of ours** by rewriting this whole set as $K’$={$A^\\dagger b +u \\in K, u \\in \\text{Ker} A$}, which is associated with the self-concordant barrier $u \\to \\phi(A^\\dagger b +u)$. We will clarify this result in the updated version of the paper.\n\nHowever, we emphasize that the **main case of interest in [2] is the case where $K$={$x: \\ell < x <u$}**. Indeed, the authors state in their introduction that their algorithm is ‘[efficient] when $K$ is a product of convex bodies $K_i$, each with small dimension’, i.e., $\\mathfrak{g}$ is a block-diagonal matrix. Indeed, the design of their numerical integrator heavily relies on the computation of the Cholesky decomposition of the matrix $A \\mathfrak{g}^-1 A^T$, which is efficient in the case where $\\mathfrak{g}$ is a (block-)diagonal matrix. In addition, through their paper, the authors consider the particular case where the $K_i$ are 1D, both in theory and in the experiments. The efficiency of the whole procedure in this specific case is justified in Appendix B.2.3. For these reasons, we highlighted this framework of [2] in the related work section.\n \nWe will change our discussion regarding [2] (lines 302-306) as follows: \n\n'On the other hand, Kook et al. (2022a) integrate the Hamiltonian dynamics via implicit schemes without any “involution checking step” in a similar self-concordant setting. In particular, they consider a convex body $K$ equipped with a self concordant barrier $\\phi$, combined with a linear equality constraint $Ax=b$. This setting is a special case of our framework (see A1 and A2) by rewriting this whole set as $K’$={$A^\\dagger b +u \\in K, u \\in \\text{Ker} A$}, equipped with the self concordant barrier $u \\to \\phi(A^\\dagger b +u)$. Kook et al. (2022a) provide an efficient implementation of their algorithm in the case of a convex bounded manifold of the form $\\mathsf{M}$={$x \\in \\mathbb{R}^d : Ax = b, \\ell < x < u$}'.\n\n> 4. Presentation\n\nWe refer to the general response for more details about our presentation. We hope that these explanations clarify our notation and we are ready to further update our paper if the reviewer has specific concerns regarding our terminology.\n\n\n[1] Condition-number-independent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators. Kook et al., 2022b.\n\n[2] Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Kook et al., 2022a.\n"}, {"review_id": "ZoQQxNrkb7", "reviewer": "Reviewer_iZM7", "summary": "The paper proposes a new algorithm for sampling on a constrained space via a Hamiltonian Monte Carlo algorithm, which uses a Riemannian Manifold HMC algorithm, with the Hessian metric given by a self-concordant barrier. This generalizes and improves upon existing work in constrained sampling, while removing a source of bias by adding an “involution checking step”.", "questions": "How difficult is it to establish any non-asymptotic convergence results in this context? What would be the expected rates, and would they differ (in the order of dependencies) with Kook et al.?\n\nAs noted under “weaknesses”, the experiments on “real-world data” do not seem fully convincing. Can the authors elaborate further on why this would be the case, and if there was extensive tuning of the parameters involved? It seems from the Appendix that the parameters of both algorithms were not extensively tuned, which then begs the question whether this is really a fair comparison.\n\nIt is curious that a (constrained) Gaussian distribution is chosen for the synthetic data, when Kook et al. consider a uniform distribution instead. Their choice allows for some simple tests of uniformity. Does considering a uniform distribution make any appreciable difference in the output?\n\nThe discussion in Appendix J is quite useful and I think some more of these details could be also summarized in the main text.\n\nIt is my personal preference, but I believe text in the paper should not be highlighted. Use boldface or perhaps underlining to emphasize the text in Ll. 192-193 and in Algorithm 1.\n\nL. 91 “motivates to” -> “motivates us to”\n\nL. 170 “the definition domain” -> “the domain”\n\nL. 213 “enables” -> “enables one/enables us”\n\nThe authors should be careful about the formatting of some equations in the Appendix, since they exceed the column width.\n", "limitations": "None beyond those raised in earlier sections.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The intuition for the algorithm, including the choice of self-concordant barrier, is well-informed and intuitive, although it has appeared in some form in Kook et al. (2022). Thus, I will discuss many of the strengths and weaknesses in relation to the prior work of Kook et al. (2022).\n\nThis paper possesses two generalizations of Kook et al. Firstly, the primary problem setting in Kook et al. (Eq 11 of that paper) is less general than that in this paper. The proof for this method is also cleaner, and purports to avoid some technical issues present in the proof Kook et al. (I did not carefully verify this point). Furthermore, it avoids the problems of bias presented in Kook et al., which is carefully justified by a Theorem.\n\nThe method is clearly unbiased when evaluated on the hypercube and simplex, by a margin that is usually many standard deviations.\n\nThe broad details of the proof seem to be correct and well-written, although I did not check every detail.\n", "weaknesses": "No rigorous non-asymptotic theory was presented, although this is also true of Kook et al., likely because of the complexity of these algorithms.\n\nIt is difficult to draw firm conclusions from the real-world instances in the experiments, since the method does not show improvement in all/a clear majority of cases.\n\nThe notation in this paper is a bit cumbersome, particularly in the Appendix. Nonetheless this is understandable given the complexity of the algorithm, and the authors did a good job in reducing the notational burden within the main text.\n\nTo summarize, it was difficult to review all the details/proofs in this manuscript, given its length and technical nature. Nonetheless, I feel like it presents an important development within constrained sampling and merits publication. This informs my final score.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments, we appreciate your acknowledgment of the paper’s merits.\n\n> No rigorous non-asymptotic theory was presented, although this is also true of Kook et al., likely because of the complexity of these algorithms.\n\n> How difficult is it to establish any non-asymptotic convergence results in this context? What would be the expected rates, and would they differ (in the order of dependencies) with Kook et al.?\n\nWe believe that the non-asymptotic convergence results of [1] could be extended to our framework, i.e., comprising the involution checking step. However, such a study would be very technical and is outside of the scope of the current paper. We conjecture that the convergence rate obtained in [1] would also hold in our case. In particular, we do not expect a change in the dependency of this rate with respect to the hyperparameters of the algorithm.\n\n> It is difficult to draw firm conclusions from the real-world instances in the experiments, since the method does not show improvement in all/a clear majority of cases.\n\n> It is curious that a (constrained) Gaussian distribution is chosen for the synthetic data, when Kook et al. consider a uniform distribution instead. Their choice allows for some simple tests of uniformity. Does considering a uniform distribution make any appreciable difference in the output?\n\nWe refer to the general response for more details on the setting of our experiments, including a comparison with [2] on real-world data and a discussion on the setting chosen for synthetic data. In particular, we explain that the experiments on real-world data are meant to show that the involution checking step does not hurt the complexity of BHMC, and justify the choice of a 'non-uniform' target distribution for synthetic datasets.\n\n> As noted under “weaknesses”, the experiments on “real-world data” do not seem fully convincing. Can the authors elaborate further on why this would be the case, and if there was extensive tuning of the parameters involved? It seems from the Appendix that the parameters of both algorithms were not extensively tuned, which then begs the question whether this is really a fair comparison.\n\nRegarding the choice of the hyperparameters except $\\eta$ (including $h$, $\\beta$, the number of iterations in the implicit solver..), we rely on the code provided by [2].  In particular, we do not fine-tune these hyperparameters for both CRHMC and BHMC. Since we are not provided with a truthful baseline due to the high-dimensional setting, we do not tune the parameter $\\eta$ and set $\\eta=10$ as in the experiments on toy datasets.\n\n> The notation in this paper is a bit cumbersome, particularly in the Appendix.\n\n> It is my personal preference, but I believe text in the paper should not be highlighted.\n\nWe refer to the general response for more details about our presentation. We hope that these explanations clarify our notation and we are ready to further update our paper if the reviewer has specific concerns regarding our terminology.\n\n> To summarize, it was difficult to review all the details/proofs in this manuscript, given its length and technical nature.\n\nWe are aware that our work is technical, and for completeness, we provide below a short description of our main theoretical results (Proposition 2 and Theorem 3), as well as a sketch of their proof. Proposition 2 justifies Assumption A3 made on the numerical integrator. Theorem 3 proves the reversibility of Algorithm 3. Below, we give more details on these results and provide a sketch of their proof.\n\n**Proposition 2**. In this proposition, we prove that in the self-concordant setting, for any starting point $z^0$, the implicit scheme $F_h=G_h \\circ s$ is locally an involution as long as $h$ is small enough. To prove this result, we proceed as follows. We first prove that we can build a symplectic diffeomorphism $\\gamma_h$ between $z^0$ and a solution $z^1 \\in F_h(z^0)$ on a neighborhood of $z^0$ (see Lemma 23, 24-(a) and 25). This proves the first item of Proposition 2. Next, using Lemma 22, 24-(b) and 26, we prove that $F_h$ coincides with $\\gamma_h$ on a given neighborhood of $z^0$. The main difficulty when establishing Proposition 2 comes from Lemma 22, which identifies the momentum $p^0$ used to go from $x^0$ to $x^1$ and shows that it is locally invertible as a perturbation of the identity. Proving this result requires computing higher order derivatives of the Riemannian metric and controlling their smoothness.\n\n**Theorem 3**. In this theorem, we prove the reversibility up to momentum reversal of the whole Markov kernel of Algorithm 3 (including momentum sampling, ODE integration, MH filter). The critical step of this procedure is to ensure that $\\Phi_h$ preserves the volume form on some compact sets under Assumption A3, derived from Proposition 2. This is proved in Lemma 29. To do so, we have to restrict the cotangent bundle to the set of points $A_h$, where h is small enough so that $\\Phi_h$ acts like an involution. Then, the proof consists in operating a change of variable on a partition of open subsets of a given compact set. The main difficulty in this proof is the control the compact sets on the domain of definition of $\\Phi_h$ in order to perform the change of variable in Lemma 29.\n\n\n[1] Condition-number-independent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators. Kook et al., 2022b.\n\n[2] Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Kook et al., 2022a."}, {"review_id": "UlSE86NHJZ", "reviewer": "Reviewer_SWNi", "summary": "This paper propose a computational framework for sampling via Hamiltonian Monte Carlo via self-concordant barrier functions for the purposes of constrained sampling on suitably nice Riemannian manifolds, where the Riemannian metric defines the barrier function. They provide computational and theoretical guarantees on their procedure, which hinges on an \"involution checking step\", making it more theoretically sound than other papers on this matter, in particular this makes the procedure unbiased. The authors supplement their theoretical findings with computational experiments in both real-world and synthetic data. ", "questions": "- What is this [x,y] or x : y notation in lines 101-102 (and other places)? I was not able to parse this...\n\n- I was also largely confused and unconvinced by the numerical experiments... In Table 1 it appears what the authors propose is about \"as good\" as CRHMC (as in, 50/50 chance it is better performing). I also cannot understand Fig 2 at all --- what does the y-axis mean, and what is the ground truth value the algorithms are trying to estimate? Is the goal to simply be as close to MALA as possible? Seems odd, but I'm far from an expert", "limitations": "Yes", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Strong theoretical contributions relative to prior work --- kudos to the authors for fixing the discrepancies in Kook et al (2022a), and proposing their solution (note, I didn't read the proofs in this paper either, only skimmed, but the results appear reasonable. I thought the background section, etc, contained all the relevant information, and the related works section also provided ample context.", "weaknesses": "I think the presentation is a bit suboptimal, namely the paper is very dense and hard to read in some parts. Otherwise, I really enjoyed this paper. I think the technical contributions appears to be somewhat minimal, though surely it took a lot of work to justify the \"involution checking step\", but again, not an expert in this area.\n\nI also understand that Kook et al. (2022b) provided some theoretical understanding of the algorithm from the previous paper --- would the authors comment a bit more on that? Do the theoretical results in that paper suffer from the same problems as in this one, or is everything remedied and this paper is concurrent with that one?\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback and their positive evaluation.\n\n> What is this [x,y] or x : y notation in lines 101-102 (and other places)? I was not able to parse this...\n\nWe refer to the general response for more details about our presentation. We hope that these explanations clarify our notation and we are ready to further update our paper if the reviewer has specific concerns regarding our terminology.\n\n> I think the technical contributions appears to be somewhat minimal, though surely it took a lot of work to justify the \"involution checking step\", but again, not an expert in this area.\n\nWe acknowledge that our algorithmic modification is minimal compared to the algorithm of [1], but we strongly believe that the study of such implementation is of fundamental interest for the MCMC community. Indeed, in a similar way as [2], our contribution (i) highlights that the **use of well known implicit schemes for ODE integration** (see [3] for a complete reference) **combined with space constraints leads to asymptotic bias when it comes to their numerical implementation**, and (ii) proposes to solve it in a straightforward manner. More generally, we are convinced that our work is a first step for future work on designing implicit integration schemes for constrained sampling methods. \n\n> I also understand that Kook et al. (2022b) provided some theoretical understanding of the algorithm from the previous paper --- would the authors comment a bit more on that? Do the theoretical results in that paper suffer from the same problems as in this one, or is everything remedied and this paper is concurrent with that one?\n\nWe refer to the general response for a detailed comparison with [4]. In particular, we show that [4] suffers from the same shortcomings as [1] when it comes to the assumption made on the implicit integrator.\n\n> I was also largely confused and unconvinced by the numerical experiments... In Table 1 it appears what the authors propose is about \"as good\" as CRHMC (as in, 50/50 chance it is better performing). I also cannot understand Fig 2 at all --- what does the y-axis mean, and what is the ground truth value the algorithms are trying to estimate? Is the goal to simply be as close to MALA as possible?\n\nWe refer to the general response for more details on the setting of our experiments, including a comparison with [1] on real-world data and and a discussion on the setting chosen for synthetic data. In particular, experiments on real-world data are meant to show that CRHMC [1] and BHMC have roughly the same complexity, while BHMC implements an additive involution checking step (which may have an effect on the autocorrelation of the obtained Markov chain). Moreover, we explain that we take the estimators of $\\int_{\\mathsf{M}}\\langle x, \\mu \\rangle \\mathrm{d}\\pi(x)$ given by MALA and IMH as ground truth for our experiments on synthetic datasets, where $\\pi$ refers to the Gaussian target distribution with mean $\\mu\\in \\mathbb{R}^d$. In this setting, we prove that the estimator given by BHMC is closer to this ground truth than the one given by CRHMC [1]. \n\n\n[1] Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Kook et al., 2022a.\n\n[2] Geometric numerical integration. Hairer et al., 2006.\n\n[3] Multiple projection Markov Chain Monte Carlo algorithms on submanifolds. Lelièvre et al, 2022.\n\n[4] Condition-number-independent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators. Kook et al., 2022b.\n"}]

Author rebuttal:
We thank the reviewers for their insightful comments and are encouraged by their positive feedback regarding the proposition, soundness, and clarity of our work. We provide detailed responses to each reviewer but summarize here their main feedback.

### 1. Comparison with [1] and counterexample

In [1], Kook et al. analyze the convergence rate of time-discretized Riemannian Hamiltonian Monte Carlo (RHMC) in a self-concordant setting (i.e., the Riemannian metric is the Hessian of the self-concordant barrier of the manifold), either relying on the Implicit midpoint integrator (Alg. 2 in [1]) or using the Generalized Leapfrog (GL) integrator (Alg. 3 in [1]). Both of these integrators are implicit, and in practice, must be replaced by numerical schemes such as fixed-point procedure. Crucially, **in [1], the authors assume that these implicit integrators admit a unique solution given any starting point and any step-size h**. However, proving such a statement in the self-concordant setting is highly non-trivial and might not even be true, as we explain below with a simple example. 

In our paper, we do not make the assumption from [1] on the GL integrator, and rather consider a numerical solver which outputs a unique approximate solution to this implicit scheme. By doing so, **our analysis is meant to be as close as possible to the practical implementation of RHMC**. In this case, the involution checking step (ICS) is needed to ensure the measure preservation property in RHMC. 

Therefore, **the theoretical (and practical) limitation about the reversibility of CRHMC [5] is still an issue in [1]**. To go further, we show below that the assumption made by [1] and [5] may not hold in the case of the GL integrator, defined in (6) in our paper. 

Consider the 1-dimensional setting $\mathsf{M}=(-\infty, 0)$. In this case, $\mathfrak{g}(x)=1/x^2$. Let $h>0$, $X_0\in \mathsf{M}$. Assume that $\beta=1$, and consider the first iteration of BHMC. We have $\tilde{P}_1 \sim \mathrm{N}(0, 1/X_0^2)$ and $P_1^{(0)}=\tilde{P}_1-\frac{h}{2}\partial_x H_1(X_0,\tilde{P}_1)$, where $H_1$ is defined in Section 3.1. As mentioned in our paper, $\partial_x H_1$ does not depend on the momentum variable, and thus, $P_1^{(0)} \sim \mathrm{N}(\mu, 1/X_0^2)$, where $\mu=-\frac{h}{2}\partial_x H_1(X_0,\tilde{P}_1)$. Using (5) in our paper, the implicit equation in $P_1^{(1/2)}$ in the GL integrator is actually a polynomial equation of degree 2: $\frac{h}{2}X_0 (P_1^{(1/2)})^2 + P_1^{(1/2)} -P_1^{(0)}=0$.

Denote $\Delta=1+2 h X_0 P_1^{(0)}$. Then, whenever $P_1^{(0)}> -1/(2 h X_0)$, i.e., $\Delta < 0$, this equation admits no solution. Recalling that $P_1^{(0)} \sim \mathrm{N}(\mu, 1/X_0^2)$, this event occurs with positive probability, thus violating the setting of [1] and [5].

### 2. Details on the experiments

**Synthetic data**. In this setting, the target distribution $\pi$ is a constrained Gaussian distribution with mean $\mu=(10 / \sqrt{d-1}) *(0, \sqrt{d-1}, 1,…,1)^T$, variance $1$ and the domain is the hypercube or the simplex. This choice of a ‘non-uniform’ distribution is motivated by the experimental observation that the importance of the reversibility condition is **best revealed when the mass of $\pi$ is not evenly distributed on the boundary of the domain**. The choice of the target quantity $A=\int_{\mathsf{M}} \langle x, \mu \rangle \mathrm{d}\pi(x)$  is arbitrary but is enough to highlight the bias in CRHMC [5], which is corrected when using BHMC. In our experiments, we aim at computing $A$ with the best accuracy. To do so, we consider as ground truth the estimators obtained with competitive algorithms such as MALA (for the hypercube) and IMH (for the simplex). Note that these two algorithms are run 10 times longer than CRHMC [5] and BHMC (i.e., $10^6$ iterations) to ensure accurate unbiased estimation. Finally, each algorithm  (CRHMC [5], BHMC, IMH, MALA) is run 10 times to obtain confidence intervals. We present boxplots in Figure 2 (here the vertical axis refers to the estimated target quantity). 

**Real-world data**. In this setting, **we do not have access to a realistic baseline that yields an unbiased estimator**, because the dimension of the polytope is too high and running MALA or IMH would be too costly. Hence, the goal of this experiment is instead to highlight that the ICS does not hurt the convergence properties of the algorithm. More precisely, we show that our method is comparable with the one from [5] in terms of time complexity to obtain independent samples. While it is clear that adding an ICS implies a tradeoff between accuracy and complexity of the method, our results demonstrate that this tradeoff does not penalize BHMC in the considered settings.

### 3. Precisions on the notation
We precise the following notations: (i) for any vectors $u$ and $v$, $D\mathfrak{g}(x)[u,v]$ is the vector $u^T D\mathfrak{g}(x) v$, (ii) $\mathfrak{g}(x)^{−1} : D\mathfrak{g}(x)$ stands for the vector $Tr(\mathfrak{g}(x)^{-1} D\mathfrak{g}(x))$. Regarding the yellow highlighted parts: these highlights are used in Alg. 1 to emphasize our main algorithmic contribution, compared to the implementation of [5], namely the ICS. We will add a sentence to explain this in the updated version of the manuscript. In the appendix, highlights are used in Alg. 2 to emphasize our contribution compared [5] in the case of continuous RHMC and are also used in Alg. 3 to emphasize the differences with Alg. 1.

[1] Condition-number-independent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators. Kook et al., 2022b.

[2] A family of MCMC methods on implicitly defined manifolds. Brubaker et al., 2012.

[3] Monte Carlo on manifolds: sampling densities and integrating functions. Zappa et al., 2018.

[4] Multiple projection Markov Chain Monte Carlo algorithms on submanifolds. Lelièvre et al., 2022.

[5] Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. Kook et al., 2022a.



Meta-review: The paper considers the problem of sampling on a constrained space using a variant of Hamiltonian Monte Carlo algorithm. Authors suggest using a Riemannian Manifold HMC algorithm with the Hessian metric given by a self-concordant barrier. Their framework improves the existing work in the context of constrained sampling, removing the bias of the algorithm through novel techniques.

The paper is reviewed by 4 expert reviewers with the following Rating/Confidence scores: 6/4, 6/2, 6/3, 7/3. All 4 reviewers agree that the paper is well written and has significant contributions. I also agree with the reviewers and recommend accepting this paper.
Decision: Accept (poster)