1 Introduction↩︎

The aim of this paper is the design of high-order sampling methods of ergodic stochastic dynamics on Riemannian manifolds. More precisely, let \((\mathcal{M},g)\) be a smooth, complete, connected \(D\)-dimensional Riemannian manifold endowed with a metric \(g\) and let \(\nabla\) be its associated Levi-Civita connection. Let \(E_1,\dots,E_D\) be a global orthonormal frame basis (for simplicity), that is, for all \(x\) in \(\mathcal{M}\), the set \(\{E_1(x),\dots,E_D(x)\}\) is an orthonormal basis of the tangent space \(T_x \mathcal{M}\). Given a smooth and Lipschitz vector field \(F : x \mapsto \underset{d=1}{\overset{D}{\sum}} f^d(x) E_d(x)\), we consider stochastic differential equations on \(\MM\) of the following form: \[\label{equation:SDE} \mathrm{d}X(t) = F(X(t))\mathrm{d}t + \sqrt{2}\sum_{d=1}^D E_d(X(t)) \circ \mathrm{d}W_d(t) , \quad X(0) = X_0 \in \mathcal{M}.\tag{1}\] Equation 1 can be understood as a SDE with additive noise, as it rewrites (see [1]) as \[\mathrm{d}X(t) = (F+\nabla_{E_n}E_n)(X(t))\mathrm{d}t + \sqrt{2} \mathrm{d}B_\mathcal{M}(t),\] where \(B_\MM(t)\) is a Brownian motion on \(\MM\) and we use the Einstein summation notation. The class of SDEs 1 includes the celebrated Riemannian Langevin dynamics \[\label{equation:Langevin} \mathrm{d}X(t) = -\nabla V(X(t))\mathrm{d}t + \sqrt{2}\mathrm{d}B_\mathcal{M}(t),\tag{2}\] when \(F\) derives from a potential, that is, when for a smooth function \(V\colon \MM\rightarrow\R\), one has \[\label{equation:correction} F = -\nabla V - \nabla_{E_n}E_n, \quad f^d = - E_d[V] - \left\langle\nabla_{E_n}E_n \vert E_d \right\rangle.\tag{3}\]

Under growth assumptions on \(F\) (see, for instance, the Bakry-Emery criterion [2] for the Langevin case 2 ), the flow \(X(t)\) of 1 is ergodic, that is, its long-time behaviour is described by a deterministic measure \(\mathrm{d}\mu_\infty=\rho_\infty \mathrm{d}\vol_\MM\), that is absolutely continuous with respect to the Riemannian volume form \(\mathrm{d}\vol_\MM\), in the sense that for a large class of smooth functions \(\varphi\) and for all initial condition \(X_0 \in \MM\), one has \[\underset{T \to \infty}{\lim} \frac{1}{T} \int_0^T \varphi(X(t))\mathrm{d}t = \displaystyle{ \int_{\mathcal{M}} } \varphi(y) \mathrm{d}\mu_\infty(y) \text{ almost surely.}\] In the case of Riemannian Langevin dynamics 2 , the invariant measure is given by the Gibbs density \(\rho_\infty\propto e^{-V}\).

In this paper, we design new numerical methods that have a high order of accuracy for sampling the invariant measure of 1 . We consider one-step numerical integrators \({(X_n)}_{n = 0,\dots,N}\), on a discretization \(t_n = nh\) of a time interval \([0,T]\) for a time step \(h=T/N\), given by a random perturbation of the identity: \[\label{equation:method} X_{n+1} = \Psi_h(X_n), \quad X_0 \in \MM.\tag{4}\] We also consider post-processed methods \[\label{equation:post-process} X_0 \in \MM, \quad X_{n+1} = \Psi_h(X_n), \quad \overline{X}_N = \overline{\Psi}_h(X_N),\tag{5}\] where the post-processor \(X \mapsto \overline{X}\) is a random perturbation of the identity map on \(\MM\), that is applied only once at the very end of the trajectory. Assuming ergodicity of the numerical method 4 for a measure \(\mathrm{d}\mu^h_\infty=\rho^h_\infty \mathrm{d}\vol_\MM\), we propose methods of high order of accuracy \(p > 1\) for the invariant measure of 1 , that is, that satisfy for \(h\) small enough, \[\abs{ \displaystyle{ \int_{\mathcal{M}} } \varphi \mathrm{d}\mu^h_\infty - \displaystyle{ \int_{\mathcal{M}} } \varphi \mathrm{d}\mu_\infty} \leq C h^p.\] A crucial difficulty of the approach is that the integrators have to evolve on \(\MM\) in order for the measure \(\mathrm{d}\mu^h_\infty\) to be absolutely continuous w.r.t.\(\mathrm{d}\vol_\MM\). Moreover, our approach focuses on intrinsic integrators, that are, numerical methods that evolve on \(\MM\) and that do not depend on an embedding of \(\MM\) in a higher-dimensional Euclidean space or on a choice of coordinates. Among the existing extrinsic approaches, we mention the popular projection methods (see [3], [4] and references therein), which are conveniently implemented by embedding \(\MM\) in a vector space of higher dimension, but often face severe timestep restrictions. We also mention the existing intrinsic approaches for SDEs on manifolds [5]–[8], that either focus on order one of accuracy or strong convergence, or are defined on specific manifolds only (i.e., Lie groups). Our approach relies on the class of frozen-flow methods recently introduced in [9], on which we add post-processors [10]. This class of methods can be seen as a stochastic extension of the Crouch-Grossman and commutator-free Lie group methods [11]–[15]. In opposition to [5], [6], our approach generalises the Lie group methods to any smooth manifold and focuses only on the approximation for the invariant measure. An important feature of the new integrators is that they do not rely on the use of embeddings or on local coordinates. Their formulation, convergence analysis, order theory, and implementation are entirely intrinsic.

To derive methods of high order for the invariant measure, it is sufficient to use methods of weak order \(p\) as they will have at least order \(p\) for the invariant measure [16]. However, there exists many efficient methods of low weak order and high order for the invariant measure in a variety of (Euclidean) contexts [17]–[21]. Among these, the celebrated Leimkuhler-Matthews method [18] is of weak order one and second order for the invariant measure for solving 1 when \(\MM=\R^D\): \[\begin{align} H_n = X_n+ \frac{\sqrt{2h}}{2} \xi^d_n E_d,\quad X_{n+1} = X_n+ h F(H_n) + \frac{\sqrt{2h}}{2}\xi^d_n E_d,\quad \overline{X_N} = X_N +\frac{\sqrt{2h}}{2} \overline{\xi}^d_n E_d. \end{align}\] Our intrinsic approach with frozen-flow methods allows us to derive in particular a straightforward generalisation of the Leimkuhler-Matthews method to the Riemannian setting,

where only one evaluation of the vector field \(F\) is needed (see Section 2.3 for the details). We recall that, in comparison, the extrinsic approach with projection methods from [4] requires four evaluations of \(F\), has a technical implementation, and faces stability issues.

The derivation of the order conditions for the invariant measure relies on exotic Lie-Butcher series and on the Weitzenböck connection algebra. The exotic Lie-Butcher series were introduced in [9] for deriving weak high-order estimates. They generalise naturally the exotic Butcher series [22]–[25] and the Lie-Butcher series [13], [15], [26]. While such formalism is necessary for dealing with the intricate calculations of order conditions, we emphasize that the algebraic objects considered here are interesting beyond their numerical use. We mention in particular the works [9], [25], [27] that study the geometric and algebraic properties of exotic forests and the works [28]–[32] that draw strong links between planar trees and computations in specific connection algebras. The numerical integrators considered here rely on elementary geometric operations (geodesic, parallel transport) for a specific curvature-free connection, called the Weitzenböck connection. This choice allows us to use planar trees for the calculations, to which we add the exotic feature for representing the Laplace-Beltrami operator in terms of forests. Thanks to this algebraic formalisation, we provide a methodology for the design of intrinsic integrators of any high-order for the invariant measure. More precisely, we define the integration by parts of planar exotic forests, thus generalising [22], [24], [25] to the manifold setting, and we describe the order conditions for the invariant measure by a character on a modified shuffle algebra of exotic forests. Contrary to the Euclidean case where such theory naturally leads to modified equations of any order and stochastic backward error analysis for the invariant measure, there are profound algebraic differences on manifolds that make intrinsic stochastic backward error analysis challenging.

The paper is organised the following way. We present in Section 2 the main results of the paper. After presenting the notation and main assumptions, we give a convenient characterisation of the high-order for the invariant measure, generalising the Euclidean works [10], [33], [34]. Then, we apply this criterion to derive new simple methods of second order for the invariant measure. The proof of the criterion for high order for the invariant measure is presented in Section 3. In Section 4, we apply the exotic Lie-Butcher series formalism for the derivation of order conditions for the invariant measure, and we extend the integration by parts on exotic forests to the planar context. We present numerical experiments on a variety of classical manifolds in Section 5 and we present outlooks and future works in Section 6.

2 Preliminaries, characterization of the invariant measure and new numerical methods↩︎

2.1 Notation and main assumptions↩︎

For \(v \in T\MM\), let \(\abs{v}=\sqrt{g(v,v)}\) be its Riemannian norm. For \(o \in \MM\). Then we denote by \(r = d(o,\cdot)\) the Riemannian distance map. It is 1-Lipschitz and \(r^2\) is smooth on \(\MM \setminus Cut_o\). For all \(\varphi \in C^\infty\), \(\mathrm{d}\varphi\) is its differential and \(\abs{\mathrm{d}\varphi}\) its Riemannian norm. We denote by \(v[\varphi](x)\) the differential \(\mathrm{d}\varphi(x) \cdot v\) of \(\varphi\) in the direction of \(v \in T_x \MM\) at the point \(x \in \MM\). Let denote the set of test function by \(C_P^\infty(\MM)\), the set of smooth functions whose derivatives of all orders have polynomial growth and satisfy estimates of the form \[\abs{E_{d_q}[\dots E_{d_1}[\varphi]\dots]}(x) \leq C(1+{r(x)}^K) , \quad \text{ for } q=0,1,\dots\]

Let \(F : x \mapsto \underset{d=1}{\overset{D}{\sum}} f^d(x) E_d(x)\) be a vector field and its unique decomposition in the orthonormal frame, then \(F\) is an element of \(\mathfrak{X}_P(\mathcal{M})\) if \(f^d \in C_P^\infty(\MM)\) and if its components are Lipschitz, that is, \(\abs{\mathrm{d}f^d} \leq C\). Equipped with the Jacobi bracket \({[-,-]}_J\), the space \((\mathfrak{X}_P(\mathcal{M}),{[-,-]}_J)\) is a Lie algebra.

In the spirit of [9], [31], [35], let the differential operators on \(\MM\) be given by the tensor algebra \(T(\mathfrak{X}_P(\mathcal{M}))\) whose product is denoted by \(\cdot\). Differential operators act on functions \(\varphi\) by: \[(X_1 \cdots X_n)\vartriangleright \varphi = x^{i_1}_1 \cdots x^{i_n}_n E_{i_1} [ \ldots E_{i_n}[\varphi] \ldots ],\] where \(X_j=\sum_{i_j=1}^D x^{i_j}_j E_{i_j}\) is the unique decomposition of the vector field \(X_j\) in the frame basis.

The product \(\vartriangleright\) extends to \(T(\mathfrak{X}_P(\mathcal{M}))\) by the so-called Guin-Oudom process [36], [37], which extends straightforwardly in a curvature-free setting (see [9]). The space of differential operators \((T(\mathfrak{X}_P(\mathcal{M})),\cdot, \Delta_{\sqcup\mathchoice{\mkern-4.8mu}{\mkern-4.5mu}{\mkern-4.5mu}{\mkern-4.5mu}\sqcup})\) is a Hopf algebra when equipped with the deshuffle coproduct [31], [35]. In the context of a connection with constant torsion, \((T(\mathfrak{X}_P(\mathcal{M})),\cdot, \Delta_{\sqcup\mathchoice{\mkern-4.8mu}{\mkern-4.5mu}{\mkern-4.5mu}{\mkern-4.5mu}\sqcup}, \vartriangleright)\) yields a post-Hopf algebra [38]. Note that taking into account the action of differential operators on functions would naturally yield post-Hopf algebroid structures [9], [39].

The generator \(\mathcal{L}\) of the SDE 1 is defined by \[\label{equation:generator} \mathcal{L} \varphi = \underset{d=1}{\overset{D}{\sum}} (f^d E_d + (E_d\cdot E_d)) \vartriangleright\varphi.\tag{6}\] In the case of the Langevin equation 2 , the generator is expressed as \[\label{equation:generator32Langevin} \mathcal{L} \varphi = - \nabla V [\varphi] + \Delta \varphi.\tag{7}\] To ensure that the equation 1 is well-posed and that the solution does not blow up, we shall assume the following.

There exist various criteria which ensure that inequality ?? is satisfied. The first criterion is compactness, which can be used for classical manifolds such as the sphere \(\mathbb{S}^n\) and the Lie group \(SO_q(\mathbb{R})\). This condition being seldom satisfied, a second classical criterion is the Bakry-Émery criterion from [2]. This criterion applies to the case of equation 1 , where the generator writes as 7 . One says that the potential \(V\) satisfies the criterion if there exists \(\kappa \in \mathbb{R}\) such that \[\Ric + \Hess (V) \geq \kappa.\] The term \(\Ric\) designs the Ricci tensor on \((\mathcal{M},g)\) and is defined as the trace of the Riemann tensor. Note that this handy criterion has been adapted in [40] for equation 1 with a lower bound of the operator \(\Ric-\nabla F\). The moment conditions from [41] implies stochastic completeness.

Under Assumption 1, the equation 1 generates a Markovian semigroup on \(C_P^\infty(\MM)\). For all \(\varphi \in C_P^\infty(\MM)\), the function \(u : (t,x) \mapsto \mathbb{E} \left[ \varphi(X(t))\vert X(0) = x \right]\) satisfies the Kolmogorov equation \[\label{equation:Kolmogorov} \partial_t u = \mathcal{L}(u), \quad u(0,x) = \varphi(x).\tag{8}\] To the contrary to [9], we are interested in the behaviour in long time of the solution of 1 .

In the spirit of [34], the ergodicity of the stochastic process is a prerequisite of our analysis.

In the case of the Riemannian Langevin equation 2 , Assumption 2 is automatically satisfied and the density \(\rho_\infty\) is explicitly given by the Gibbs measure: \[\label{equation:Gibbs} \rho_\infty = e^{-V} / Z, \quad Z = \displaystyle{ \int_{\MM} } e^{-V} \mathrm{d}\vol_\MM.\tag{9}\]

The following assumption from [9], automatically satisfied on \(\mathbb{R}^d\), ensures the regularity of the semigroup 8 .

Hence, taking the limit as \(t \to \infty\) of the time integration of Equation 8 leads to, under Assumption 4, \[\label{equation:lim95u} \underset{t \to \infty}{\lim} u(t,x) = \varphi(x) + \int_0^\infty \mathcal{L} u(s,x) \mathrm{d}s.\tag{10}\] The following assumption is an ergodicity condition concerning convergence to the invariant measure. We introduce the following norm and seminorm: \[\left\Vert\psi \right\Vert_{C^k} = \underset{\underset{0 \leq \abs{j}\leq k}{j=(j_1,\dots,j_D)}}{\sup} \underset{x \in \mathcal{M}}{\sup} \abs{\partial_j \psi(x)} \text{ and } \abs{\psi}_{C^k} = \underset{\underset{1 \leq \abs{j}\leq k}{j=(j_1,\dots,j_D)}}{\sup} \underset{x \in \mathcal{M}}{\sup} \abs{\partial_j \psi(x)}.\]

The Taylor expansion of the semigroup 8 satisfies the following.

Introduce the Taylor-Talay-Tubaro expansion [16] of the numerical scheme 4 and an assumption on its moments.

2.2 High-order characterization of the invariant measure↩︎

Let us now state one of the central results of this work. We introduce new sufficient conditions to approximate the invariant measure, adapted from [34]. This result is key for designing the high-order numerical schemes 4 presented in Subsection 2.3. Proofs of Theorems 5 and 7, adapted from [10], [33], are posponed to Section 3.

This characterization is extended to post-processed method 5 , in the spirit of [10].

2.3 Post-processed frozen-flow methods of high order for sampling ergodic dynamics↩︎

In this section, we leverage Theorem 5 to derive a system of conditions whose resolution yields three new order-2 methods for accurately sampling the invariant measure. Higher orders can be obtained by extending the Taylor expansion to the desired degree. Our analysis focuses on frozen-flow schemes, which generalize classical Lie group methods [12], [13]. See Section 4 for details of the Butcher series formalism used on the Talay-Tubaro expansion.

Following [10] and the frozen-flow methods [9], we consider the following class of methods and we define with the same way the class for the post-processor 5 \[\overline{X}_N = \overline{\Psi}_h(X_N).\]

\[\begin{align} H_n^i = & \exp \left( \left( h Z^0_{i,j,K} f^d(H_n^j) + \sqrt{h} Z_{i,K}^d \xi^{d}_n \right) E_d \right) \cdots \notag \\ & \cdots \exp \left( \left( h Z^0_{i,j,1} f^d(H_n^j) + \sqrt{h} Z_{i,1}^d \xi^{d}_n \right) E_d \right) X_n, \notag \\ X_{n+1} = & \exp \left( \left( h z^0_{i,K} f^d(H_n^i) + \sqrt{h} z^d_{K} \xi^{d}_n \right) E_d \right) \cdots \label{equation:post-processed} \\ & \cdots \exp \left( \left( h z^0_{i,1} f^d(H_n^i) + \sqrt{h} z^d_{1} \xi^{d}_n \right) E_d \right) X_n, \notag \end{align}\tag{11}\] where the \(\xi^{d}_n,\overline{\xi^{d}_N}\) are independent standard Gaussian random variables and where the \(Z^0_{i,j,k},z^0_{i,k},z^d_{k},Z^d_{i,k},\overline{z_{k}^{d}},\overline{z^0_{i,k}},\overline{Z_{i,k}^{d}},\overline{Z^0_{i,j,k}}\) are fixed real numbers.

We now express the general method as an exotic S-serie (see Definition 9), with the coefficient map of the numerical methods \(a_w\), and the algebraic formalism of Section 4.

Since we only express our conditions in terms of exotic forests, to simplify the notation, we add the coefficients of \(\left[ \mathcal{L}, \overline{A_1} \right]\), in the post-processed case, in thoses of \(A_2\).

Using the Butcher series formalism of Section 4, we obtain a set of conditions that ensures the order 2 for the invariant measure of equation 2 .

We propose the following new post-processed frozen-flow method to solve equation 2 . In the Euclidean setting, Method 1 rewrites as the Leimkuhler-Matthews method [18] with the post-processor formulation from [10].

We propose two alternative schemes without post-processors to solve equation 2 , for the sake of comparison.

3 High order analysis for the invariant measure↩︎

In this subsection, we generalize the analysis of [33] to Riemannian manifolds and prove Theorem 5 and Theorem 7.

Let \({(P_t)}_{t\geq0}\) be the semigroup associated to the Markov process \({(X(t,x))}_{t\geq0,x\in\mathcal{M}}\), in other words, \({(P_t)}_{t\geq0}\) gives to \(\psi\) a smooth function on \(\mathcal{M}\), the unique solution \(u\) of the Kolmogorov equation 8 . In the spirit of backard analysis [33], [43], let us build a modified generator \(\tilde{\mathcal{L}_\tau}\) such that the solution of the problem \[\partial_t v(t,x) = \tilde{\mathcal{L}_\tau} v , \quad v(0,x) = \varphi(x),\] at the time \(t = \tau\), coincides with the numerical flow, that is \[\exp \left( \tau \tilde{\mathcal{L}_\tau} \right) = \id + \underset{n \geq 1}{\sum} \tau^n A_n.\] Thus the expansion yields \[\tilde{\mathcal{L}_\tau} = \underset{k = 0}{\overset{n}{\sum}} \tau^k L_k + R_n , \quad R_n = \underset{k \geq n+1}{{\sum}} \tau^k L_k .\] Define \(\pi_n : P \mapsto n! {\left({\partial_\tau}^n P\right)}_{| \tau = 0}\), then \[\begin{align} A_n & = \pi_n \left( \exp \left( \tau \tilde{\mathcal{L}_\tau} \right) \right) = \underset{k \geq 0}{\sum} \frac{1}{k!} \pi_n \left( \tau^k \tilde{\mathcal{L}_\tau}^k \right) \\ & = \underset{k = 1}{\overset{n}{\sum}} \frac{1}{k!} \pi_{n-k} \left( \tilde{\mathcal{L}_\tau}^k \right) = \underset{k = 1}{\overset{n}{\sum}} \frac{1}{k!} \pi_{n-k} \left( { \left( L_0 + \tau L_1 + \cdots + \tau^n L_n + R_n \right) }^k \right) \\ & = \underset{k = 1}{\overset{n}{\sum}} \frac{1}{k!} \pi_{n-k} \left( \underset{j = 0}{\overset{nk}{\sum}} \tau^j \underset{i_1+ \cdots +i_k = j}{\sum} L_{i_1} \cdots L_{i_k} \right) \\ & = \underset{k = 1}{\overset{n}{\sum}} \frac{1}{k!} \underset{i_1+ \cdots +i_k = n-k}{\sum} L_{i_1} \cdots L_{i_k} , \end{align}\] where we use \(\pi_p (R_n) = 0\) when \(p \leq n\). Thus \(\left(L_n\right)_{n \in\mathbb{N}}\) satisfies the recurrence relation \[L_0 = A_1 = \mathcal{L} , \quad L_n = A_{n+1} - \underset{k = 2}{\overset{n+1}{\sum}} \frac{1}{k!} \underset{i_1+ \cdots +i_k = n+1-k}{\sum} L_{i_1} \cdots L_{i_k}.\]

We construct a sequence of measure which weakly converges to the invariant measure of the numerical method 4 .

To be able to prove Theorem 5, we recall a lemma from [34], that is the central result for the proof of the theorem.

With this analysis, we now prove the Theorem 5.

We conclude this section by the proof of Theorem 7.

4 Intrinsic order conditions for the invariant measure↩︎

In this section, we use the algebraic framework of exotic Lie-Butcher series introduced in [9] to derive the order conditions for sampling the invariant measure. After recalling the notion of planar exotic forests, we present a new operation \(\mathop{\mathrm{IBP}}\) on forests, derived from the Riemannian integration by parts, to express the adjoint of the operators \(A_n\) and generalise [22], [24]. The approach is successfully applied through Theorem 5 and Theorem 7 to derive the conditions of Theorem 12. We then discuss the form of the exotic Lie-Butcher series describing the adjoint of the operators \(A_n\) with irreducible forests.

4.1 Planar exotic forests↩︎

We recall the structure of planar exotic forests, in the spirit of [9] which generalize standard B-series [43]–[45] to SDEs on manifolds.

From [9], the operator \(\mathbb{F}^F: (\mathcal{EF},\cdot) \to (T(\mathcal{\mathfrak{X}}_P(\MM)),\cdot)\) is a Hopf algebra morphism. In particular, \(\mathbb{F}^F\) satisfies \[\mathbb{F}^F(\pi_1\cdot \pi_2)=\mathbb{F}^F(\pi_1)\cdot \mathbb{F}^F(\pi_2).\]

In the following, we denote by \(a_w\) the coefficient map of the numerical flow, that is the unique coefficient map \(a_w : \mathcal{EF}\to \mathbb{R}\) such that \[\mathbb{E} \left[ \varphi(X_1)\vert X_0 = x \right] = S_h(a_w) \vartriangleright \varphi .\]

4.2 Integration by parts for the invariant measure↩︎

The derivation of the order conditions for the invariant measure is performed using multiple integration by parts of the Taylor-Talay-Tubaro expansion against the invariant measure. In this context, the use of an orthonormal frame proves crucial to ensure that the operations rewrite with planar exotic forests.

Let us now identify the exotic forests whose integration by part does not write with on exotic forests.

The integration by parts of planar exotic forests is defined as the following, which generalises the non-planar map in [24].

We now derive an explicit expression of the adjoint operator \(A_j^*\) with irreducible forests. Since \(\mathop{\mathrm{IBP}}\) removes only the first numbered root, we iterate \(\mathop{\mathrm{IBP}}\) until the output is a linear combination of irreducible forests.

The algebraic characterization of Theorem 5 rewrites naturally with planar exotic forests.

In order to compute the second order condition from Theorem 12, Table 1 gives the reduction of second order exotic forests.

4.3 Exotic S-series expression of the adjoint operators↩︎

Let us distinguish the coefficient maps for the weak error \(a_w : \mathcal{EF}\to \mathbb{R}\) and for the invariant measure. Denote by \(\left\langle\cdot \vert\cdot \right\rangle\) the bilinear pairing over \(\mathcal{EF}\) such that, for \(\pi_1, \pi_2 \in EF\), \[\left\langle\pi_1 \vert\pi_2 \right\rangle = 1 \text{ if } \pi_1 = \pi_2 \neq 0 \text{ and } 0 \text{ else}.\]

It is important to mention that the conditions for the invariant measure are not independent and satisfy relations, analogously to the weak context [9]. However, the map \(a_{\mu_\infty}\) is not a character for the shuffle product as one has for instance \[{a_{\mu_\infty}( \tikz[planar forest ] { \node[b] at (0.0, 0.4) { } ; })}^2 = 2 a_{\mu_\infty}( \tikz[planar forest ] { \node[b] at (0.0, 0.4) { } ; \node[b] at (1.0, 0.4) { } ; }- \tikz[planar forest ] { \node[b] at (0.0, 0.4) { } ; \node[lc] at (1.0, 0.4) { 1 } ; \node[lc] at (2.0, 0.4) { 1 } ; }) , \quad {a_{w}( \tikz[planar forest ] { \node[b] at (0.0, 0.4) { } ; })}^2 = 2 a_{w}( \tikz[planar forest ] { \node[b] at (0.0, 0.0) { } child {node[b] at (0.0, 1.0) { } } ; }).\]

5 Numerical experiments↩︎

5.1 Ergodicity and long-time sampling on \(SO_p\)↩︎

Our first experiment focuses on \(\MM=SO_p\) the compact Lie group of special orthogonal matrices of size \(p \times p\). Let \({\left(A_d\right)}_{d=1,\ldots,D}\) be an orthonormal basis of its Lie algebra \(\mathfrak{so}_p\) of skew-symmetric matrices for the metric \(g_x(A,B) = \Trace \left( A^T B \right)\), \(A,B \in T_x \MM\). Let \(E_d(y) = A_d y\) be the associated orthonormal frame basis. As the manifold is smooth and compact, the assumptions of Theorem 7 are satisfied and our analysis applies.

Inspired by the experiments in [9], we compare the new methods of high order for the invariant measure to the frozen-flow Euler method ?? and the so-called SFF2 method from [9], which uses two random variables and is of weak order 2. We consider dynamics with the test function \[\label{eq:num95test95fct} \varphi(x) = \exp \left( - \frac{g_x (x-I_p,x-I_p)}{2p} \right).\tag{13}\] The reference solution is taken as the mean value of the SFF2 method [9] with the reference time step \(h=2^{-12}\). We compare the order of convergence in long time for two potentials of the form \(V = P( \left\Vert\cdot - I_p \right\Vert^2 )\) where \(P\) polynomial: a quadratic potential well with \(P = 10 X\) and a sextic potential with multiple local minima using the cubic polynomial \(P = \frac{X}{2} (X^2-9X+24)\). Starting from the initial condition \(X_0 = - \left(\begin{smallmatrix} 1&0&0\\0&0&1\\0&1&0 \end{smallmatrix}\right)\), we observe in Figure 5 the ergodic behaviour of the dynamic as the trajectories explore a neighbourhood of the identity matrix after an initial transient phase.

We observe the error curves for the invariant measure in Figure 6. We observe that the new methods 1, 2 and 3 exhibit the expected second-order behavior, with a reduced cost compared to SFF2, which confirms our theoretical results. Moreover, Method 1 rapidly reaches the Monte-Carlo threshold in the quadratic case, hinting that it may have a higher order of convergence for this specific potential.

5.2 Von-Mises Fisher dynamics on the sphere \(\mathbb{S}^2\)↩︎

Our second experiment focuses on the 2-dimensional sphere \(\MM=\mathbb{S}^2\). As our new methods rely on orthonormal bases, we consider the standard coordinates on \(\mathbb{S}^2\) minus the poles \(( 0 , 0 , \pm 1)\): \[y_{\theta,\phi} = \begin{pmatrix} \cos (\theta) \cos (\phi) \\ \cos (\theta) \sin (\phi) \\ \sin (\theta) \end{pmatrix} \in \mathbb{S}^2, \quad \theta \in \left(-\frac{\pi}{2},\frac{\pi}{2}\right), \quad \phi \in \left[0,2\pi\right),\] with the associated orthonormal frame: \[E_1(y_{\theta,\phi}) = \begin{pmatrix} -\sin (\theta) \cos (\phi) \\ -\sin (\theta) \sin (\phi) \\ \cos (\theta) \end{pmatrix}, \quad E_2(y_{\theta,\phi}) = \begin{pmatrix} - \cos (\theta) \sin (\phi) \\ \cos (\theta) \cos (\phi) \\ 0 \end{pmatrix}.\] Using the symmetry \(S : (x,y,z) \mapsto (z,y,x)\), we define a second frame on \(\mathbb{S}^2\) minus \(( \pm 1 , 0 , 0)\). For one step of the method, we use the first frame basis if \(\abs{\sin(\theta)}\leq 0.6\), and the second frame else. This ensures a bounded Lipschitz constant for the chosen frame in each case and the assumptions of our analysis are satisfied. Following the experiments of [4], [7], [9], we consider the potential \(V : (x,y,z) \mapsto -25z\), the associated vector field (containing the Ito correction 3 ) \[F = - \nabla V - \nabla_{E_1} E_1 - \nabla_{E_2} E_2,\] and the test function \(\varphi : (x,y,z) \mapsto z^2\). This choice of potential confines the trajectories around the north pole of \(\mathbb{S}^2\) (see Figure 7). The reference solution is chosen as the mean value of the SFF2 method [9] with the reference time step \(h=\frac{1}{ 25 \times 2^{10}}\).

We compare the new methods of high order for the invariant measure with the methods of [9]. The error curves for the invariant measure in Figure 7 display the expected orders of convergence. In particular, Methods 2 and 3 exhibit second-order behavior, while Method 1 reaches the Monte-Carlo threshold instantly.

6 Conclusion↩︎

In this paper, we designed new numerical schemes of high accuracy tailored for the long time sampling of ergodic stochastic dynamics on Riemannian manifolds. We generalized the Euclidean analysis [33], [34] to characterize the intrinsic approximation for the invariant measure on manifolds, and presented new algebraic operations on exotic forests [9] to compute the order conditions in long time. The analysis applies naturally to Riemannian Langevin dynamics and the new methods, which generalise in particular the Leimkuhler-Matthews method, show second order of convergence for a reduced cost on numerical experiments.

The present paper opens several avenues for future research. The new discretisations could naturally be combined with popular sampling techniques, such as perturbations [21], [48]–[50]. The analysis for the invariant measure, as well as the associated algebraic formalism, will be extended for the design of modified equations of arbitrarily high order for sampling SDEs on manifolds and for the challenging derivation of intrinsic stochastic backward error analysis, generalising [22], [25]. Similar to the Euclidean setting [51], one could extend the analysis to general SDEs with multiplicative noise and the creation of higher order methods. Our approach follows the one of deterministic Lie-group methods and thus relies on a curvature-free connection, unrelated to the natural Levi-Civita connection on Riemannian manifolds. We will extend the frozen-flow methods so that they rely only on Riemannian operations. The study of such new methods is already open in the deterministic setting and relies on the challenging general understanding of the connection algebra [30], [52], [53]. Finally, the new algebraic formalism of planar exotic series is interesting in itself and could be studied for its universal combinatorial, algebraic, and geometric properties [27], [54]–[57], but also for its potential applications in different fields, in the spirit of the use of Butcher series and their extensions in rough paths [58], [59], renormalisation theory [60], variational calculus [61], [62], approximation of PDEs [63], mathematical physics [64], and wave kinetic dynamics [65]. These projects will be studied in upcoming works.

The authors would like to thank Baptiste Huguet for helpful comments on a previous version of the present work. The authors acknowledge the support from the French program ANR-25-CE40–2862–01 (MaStoC — Manifolds and Stochastic Computations). Experiments presented in this paper were carried out using the Abaca testbed, supported by Inria (see https://abaca.inria.fr).

References↩︎