Efficient Random Walks on Riemannian Manifolds

According to a version of Donsker's theorem, geodesic random walks on Riemannian manifolds converge to the respective Brownian motion. From a computational perspective, however, evaluating geodesics can be quite costly. We therefore introduce approximate geodesic random walks based on the concept of retractions. We show that these approximate walks converge in distribution to the correct Brownian motion as long as the geodesic equation is approximated up to second order. As a result we obtain an efficient algorithm for sampling Brownian motion on compact Riemannian manifolds.


Introduction
Probabilistic models in continuous time with geometric constraints lead to Brownian motion and stochastic differential equations (SDEs) on Riemannian manifolds (M, g).The theory of Brownian motion and SDEs on manifolds has been extensively studied (see, e.g., [16,20])and many results from the Euclidean setting can be generalized to manifold-valued SDEs.One example is Donsker's theorem.Here, one considers a geodesic random walk for which the update step at the current location x i is obtained by uniformly sampling a unit tangent vector v ∼ Unif(S xi ) with S x := {u ∈ T x M : ∥u∥ g = 1} and following a geodesic path by setting x i+1 := Exp xi (εv) for some small enough parameter ε, where v := √ mv, m denotes the dimension of the manifold, and Exp x : T x M → M is the exponential map.We will discuss the reason for scaling by √ m below.Jørgensen proved that on complete Riemannian manifolds that satisfy mild assumptions (which in particular hold for compact manifolds), the scaling limit of such a geodesic random walk is the Brownian motion on (M, g), see [21].This result has recently been generalized to the setting of Finsler manifolds, see [26].
In principle, it is thus possible to simulate Brownian motion on Riemannian manifolds by computing geodesic random walks -analogously to how one would use Donsker's theorem in the Euclidean case.Unfortunately, this procedure is not computationally efficient on manifolds since the solution of a non-linear second-order geodesic ODE has to be computed in every step.With an eye on applications, we are interested in efficient and simple, yet convergent, methods in order to approximate Brownian motion on Riemannian manifolds.
We introduce a class of algorithms that yield approximate geodesic random walks by making use of retractions that approximate the exponential map and have been introduced in the context of numerical optimization on manifolds, see [1][2][3]8].We show that retractions which approximate the geodesic ODE up to second order converge to Brownian motion in the Skorokhod topology, see Section 3. We consider two prevalent scenarios, where random walks based on such second-order retractions can be computed efficiently: The case of parameterized manifolds and the case of implicitly defined manifolds, see Section 2. Our approach generalizes the setting of approximate geodesic walks on manifolds with positive sectional curvatures presented in [27].Moreover, our method works for arbitrary dimensions, and its cost only depends on the evaluation of the respective retraction.
A retraction in our sense can be understood as both, an approximation of the geodesic ODE, or, equivalently, as a 2-jet in the sense of [5], which allows for treating coordinate-free SDEs.Armstrong and King developed a 2-jet scheme for simulating SDEs on manifolds in [6], which uses an unbounded stepsize in each iteration.Different from our retraction-based approach, which is based on bounded stepsizes, the case of unbounded stepsizes hampers efficiency (for the case of local charts) or might even prohibit an implementation (for the case of implicitly defined submanifolds when using a step-and-project approach).Moreover, our random-walk method is based on a Donsker-type theorem, while their approach improves Euler-Maruyama approximations -resulting in a different type of convergence in terms of the respective topologies.
Outside the setting of 2-jets, one typically considers SDEs on embedded manifolds M ⊂ R n and first solves them in ambient Euclidean space, followed by projecting back to M using Lagrange multipliers, see, e.g., [7,25].Projection-based approaches, however, do not immediately extend to the case of parameterized manifolds.Within our unified framework based on retractions, we cover both, the case of embedded manifolds and the case of parameterized manifolds, see Sections 2.1 & 2.2.In particular, our approach leads to an efficient and convergent numerical treatment of Brownian motion on Riemannian manifolds (including drift -see Remark 2.4).
A different scenario where geodesic random walks are commonly used is the problem of sampling the volume measure of polytopes (defined as the convex hull resulting from certain linear constraints).Originally, Euclidean random walk methods had been considered in this context, see [30].Although the Euclidean perspective is natural in this situation, small stepsizes are required in order to reach points near the boundary, which hampers efficiency in higher dimensions.To circumvent this issue, Lee and Vempala [24] introduced a Riemannian metric based on the Hessian of a log-barrier -to the effect of re-scaling space as the walk approaches the boundary.In their implementation, they used collocation schemes for solving the requisite ODEs in order to simulate geodesic random walks.We leave for future research the question whether such approaches can be reformulated in terms of retractions.
Furthermore, the problem of sampling on manifolds is commonly addressed by Markov chain Monte Carlo (MCMC) methods, see, e.g., [10][11][12]31].Different from our focus, MCMC methods are typically concerned with obtaining samples from a given probability distribution on M .Our algorithm, however, yields approximate sample paths of a Brownian motion on a manifold.We are not aware of any results concerned with proving that an MCMC method samples the correct dynamics of a given SDE.Although not the main focus of our exposition, we show that our method also correctly recovers the stationary measure of a geodesic random walk in the limit of stepsize tending to zero, see Theorem 3.2.

Retraction-based random walks
Throughout this exposition, we consider m-dimensional compact and orientable Riemannian manifolds (M, g) without boundary.We use retractions to approximate the exponential map on such manifolds.Definition 2.1.Let Ret : T M − → M be a smooth map, and denote the restriction to the tangent space T x M by Ret x for any x ∈ M .Ret is a retraction if the following two conditions are satisfied for all x ∈ M and all v ∈ T x M : (1) Ret x (0) = x, where 0 is the zero element in T x M and (2) A retraction is a second-order retraction if it additionally satisfies that for all x ∈ M and for all v ∈ T x M one has that (1) D dτ where D dτ ( d dτ γ(τ )) denotes covariant differentiation of the tangent vector field γ(τ ) = d dτ γ(τ ) along the curve γ (following standard notation, see, e.g., [13]).If M is a submanifold of a Euclidean space, Equation ( 1) is equivalent to where N x M denotes the normal bundle of M .Consequently, for the case of submanifolds, one has Clearly, the exponential map is itself a retraction.The main benefit of retractions is, however, that they can serve as computationally efficient approximations of the exponential map.Using retractions yields Algorithm 1 for simulating random walks on a Riemannian manifold.Set v := √ mv 5: x ε i := Ret x ε i−1 (εv) 6: end for 7: Return: (x ε j ) 0≤j≤N Notice that the uniform sampling in Step 3 of the above algorithm is performed with respect to the Riemannian metric restricted to the unit sphere in the attendant tangent space.We provide details for this step in Sections 2.1 and 2.2 below.Notice furthermore that ε takes the role of a spatial stepsize, while the total physical time simulated by Algorithm 1 is given by ε 2 N ; see also the time rescaling in Theorem 2.2 below.Step 4 of the algorithm is necessary since uniform sampling on the unit sphere in a space of dimension m results in a covariance matrix 1  m Id, while convergence to Brownian motion requires the identity as covariance matrix.The normalizing factor √ m in Step 4 precisely ensures the latter.Indeed, without the renormalization in Step 4 the random walk constructed in Algorithm 1 converges to a time-changed Brownian motion.
The next theorem shows convergence of the retraction-based random walk resulting from Algorithm 1 and provides an explicit expression for the generator of the resulting limit process: Theorem 2.2.Consider the sequence of random variables (X ε i ) i∈N constructed in Algorithm 1 (with N = ∞).The continuous-time process X ε := X ε ⌊ε −2 t⌋ t≥0 converges in distribution to the stochastic process with generator dv in the Skorokhod topology (see Section 3 for details).Here S x denotes the unit sphere in T x M , ω m is the volume of this sphere, and ∆ g is the Laplace-Beltrami operator on (M, g).
We defer the proof to Section 3. Notice that Theorem 2.2 immediately implies a necessary and sufficient condition for a retraction-based random walk to converge to Brownian motion on a Riemannian manifold (M, g): Corollary 2.3.Using the same notation as in Theroem 2.2, the continuous-time process X ε converges in distribution to the Brownian motion on M if and only if the Laplacian (also known as the tension) of Ret x : T x M − → M vanishes at 0 ∈ T x M for every x ∈ M .In particular, this holds for second-order retractions.
Proof.By definition, the Laplacian of Ret x : T x M − → M is defined as where ∇d Ret x is the Hessian of the C ∞ -map Ret x .Clearly, the Laplacian of Ret x vanishes at 0 ∈ T x M if and only if dv .
Applying Theorem 2.2 proves the first claim.If Ret x is a second-order retraction, the second claim immediately follows since Thus, random walks generated by second-order retractions converge to the Brownian motion on (M, g).
Remark 2.4.Drift can seamlessly be incorporated into our approach.Indeed, when modeling drift by a vector field X, then Step 4 of Algorithm 1 needs to be altered according to v := √ mv + εX| x ε i−1 .The stochastic process resulting from this altered version of Algorithm 1 converges to the process with generator Lf + df (X), where Lf is the generator given by (2).The proof of this fact follows the proof of Theorem 2.2 with obvious modifications.
In the next two subsections we consider two concrete examples of second-order retractions that are computationally efficient.and where the relatively open sets {ϕ i (U i )} cover M .In concrete applications, such parameterizations often arise naturally, e.g., for isometrically embedded manifolds M ⊂ R n .

Retractions based on local parameterizations. Our first example of computationally efficient retractions is based on local parameterizations of compact m-dimensional Riemannian manifolds (M, g). Consider an atlas consisting of finitely many charts
For a given x ∈ M , v ∈ T x M , and We assume that for some fixed ε > 0 and any x ∈ M , there exists a chart , where the (Γ ij ) denote the Christoffel symbols of the respective chart.This assumption is readily satisfied for any compact M , provided that ε is chosen small enough.In order to find such a sufficiently small ε in Algorithm 1, we propose to restart the algorithm with stepsize ε 2 if the condition is violated.Then we define the retraction Notice that Christoffel symbols can be computed numerically.For low-dimensional manifolds, symbolic differentiation accelerates their computation.
Lemma 2.5.The retraction defined by Equation (3) is of second order.
Proof.Let (γ(τ )) τ ≥0 denote the geodesic satisfying γ(0) = x and γ(0) = v.In the parameter domain, γ satisfies the ordinary differential equation The Taylor expansion of γ reads as τ → 0, which proves that p-Ret is indeed a second-order retraction.□ Another aspect of parameterization-based approaches concerns the computation of tangent vectors ṽ in the parameter domain that are uniformly distributed with respect to the Riemannian metric g.An efficient procedure for this task is presented in Algorithm 2.
Algorithm 2 Uniform sampling of ṽ w.r.t. the Riemannian metric 4: Compute the singular value decomposition (SVD) of the pullback metric where V i the ith column vector of V. 6: Return: ṽ Notice that Equation ( 5) is well defined since w k ̸ = 0 a.s.for all 1 ≤ k ≤ m and that z is a sample of a uniformly distributed random variable on S m−1 .It is straightforward to verify the following claim: Lemma 2.6.The vector ṽ constructed in Algorithm 2 is a unit tangent vector that is uniformly sampled with respect to the Riemannian metric (g ij ).
In summary, using ṽ as computed in Algorithm 2, setting v = ṽ in Step 3 of Algorithm 1, and using the retraction p-Ret from (3), provides an efficient method for computing random walks for parameterized manifolds.See Table 1 for a comparison between geodesic random walks and retraction-based random walks, showing that the latter are significantly faster than the former.Notice that the attendant algorithm might require to update the respective local chart during the simulated random walk.Results of a simulated random walk on a 2-torus (using a single chart with periodic boundary conditions) are presented in Figure 1.

2.2.
Projection-based retractions.Our second example of computationally efficient retractions is based on a step-and-project method for compact manifolds that are given as zero level sets of smooth functions.More specifically, we consider the setting where f : R n → R k is a smooth function, k < n, and the m-dimensional manifold in question, with m = (n − k), is given by ( 6) We assume that 0 is a regular value of f , i.e, that the differential df has full rank along M .In this case, a retraction can be defined via ) for i = 1, . . ., m, ḡij = 0 for i ̸ = j, and B is a random orthogonal matrix.In each case 10 000 steps were calculated in varying dimensions with stepsize ε = 0.1.All computations were performed on a standard laptop in Matlab, using Matlab's method ode45 for solving the geodesic equation.The Christoffel symbols were computed numerically.
where x ∈ M , v ∈ T x M , and π M (y) = arg min p∈M ∥p − y∥ denotes the (closest point) projection to M .Notice that π M (y) is well defined as long as the distance d(y, M ) is less then the so-called reach of M .The reach of M is the distance of M to its medial axis, which in turn is defined as the set of those points in R n that do not have a unique closest point on M .Smoothly embedded compact manifolds always have positive reach.
The following result is important in our setup.For a proof, see, e.g., Theorem 4.10 in [3].
It remains to specify the computation of a random unit direction v ∈ T x M as well as the implementation of the projection step.As for the former, a uniformly sampled unit tangent direction at some x ∈ M can be computed by first sampling a randomly distributed unit vector u ∈ S n−1 in ambient space, and denoting its projection to im(df | x ) by ũ.Then (7) v = u − ũ ∥u − ũ∥ yields the requisite uniformly distributed unit tangent vector.Finally, for computing the closest point projection of some given point y ∈ R n to M , consider the Lagrangian where λ ∈ R k denotes the (vector-valued) Lagrange multiplier.Then Newton's method offers an efficient implementation for solving the Euler-Lagrange equations resulting from L y .In our implementation, we let Newton's method run until the threshold |f (z)|/∥df | z ∥ < ε 3 is satisfied.Notice that efficiency results from the fact that the point y = x + τ v is close to M for x ∈ M and τ small enough.In summary, using v as computed in (7) for Step 3 of Algorithm 1, together with the retraction π-Ret from (6), provides an efficient method for computing random walks for implicitly given submanifolds.Results of a simulated random walk on the 2-torus are presented in Figure 2.
Remark 2.8.Our projection-based retraction is an instance of retractions induced by so-called normal foliations.The latter are defined as foliations in a local neighborhood of M = f −1 (0) whose leaves have dimension k and intersect M orthogonally.As shown in [32], retractions induced by normal foliations are always of second order and therefore provide alternative approaches for sampling Brownian motion on implicit manifolds.Examples of retractions induced by normal foliations include the gradient flow along −∇f and Newton retractions based on the update step δz = −(df ) † (z)f (z), where (df ) † is the Moore-Penrose pseudo-inverse of the Jacobian df .

Stationary measure.
For second order retractions, we show below that the stationary measure of the process in Algorithm 1 converges to the stationary measure of the Brownian motion, which is the normalized Riemannian volume form, see Theorem 3.2.Nonetheless, the empirical approximation 1 N N i=1 δ x ε i of the Riemannian volume form using random walkse.g., in order to obtain results as in Figs. 1 & 2 -requires the stepsize ε and the number of steps N in Algorithm 1 to be chosen appropriately.An indicator at which time the empirical approximation is close to the stationary measure is the so called δ-cover time, i.e., the time needed to hit any δ-ball.Hence, as a rough heuristic choice of ε and N , one may utilize [14,15].For a Brownian motion on a 2-dimensional manifold M , the δ-cover time C δ of M , satisfies ([15], Theorem 1.3) where A denotes the Riemannian area of M .Based on the observation that a second-order retraction converges to Brownian motion in the limit of vanishing stepsize ε, and using that the physical time of the walk resulting from Algorithm 1 is given by ε 2 N , one may set If δ is small enough, one may then choose N according to For dimension m > 2, the corresponding result in [14], Theorem 1.1 yields the heuristic choice  V (M ) denotes the volume of M , and ω m the volume of the unit sphere in R m .Such heuristics, however, do not quantify the error between the empirical measure obtained by simulating a random walk and the Riemannian volume measure.In the setting of Euclidean random walks (Y n ) n∈N with invariant distribution ν, it has been shown in [23] that under some conditions on Y the expectation of the 1-Wasserstein distance between the empirical and the invariant measure asymptotically decreases according to the law as N → ∞ with constants α, β, γ > 0 depending on the dimension m and ϑ ∈ (0, 1) depending on the process Y .Our simulations, see Figure 3, followed by a simple fit to the model given by the right hand side of Equation ( 8) suggest that a similar qualitative behavior might be present for retraction-based random walks on Riemannian manifolds.We leave a deeper analysis for future research.

Convergence
This section is devoted to the proof of Theorem 2.2, which yields convergence in the Skorokhod topology of the process constructed in Algorithm 1 as the stepsize ε → 0. The Skorokhod space D M [0, ∞) is the space of all càdlàg, i.e., right-continuous functions with left-hand limits from [0, ∞) to M , equipped with a metric turning D M [0, ∞) into a complete and separable metric space, and hence a Polish space (see e.g.[17], Chapter 3.5).The topology induced by this metric is the Skorokhod topology, which generalizes the topology of local uniform convergence for continuous functions [0, ∞) → M in a natural way.For local uniform convergence we have uniform convergence on any compact time interval K ⊂ [0, ∞).In the Skorokhod topology functions with discontinuities converge if the times and magnitudes of the discontinuities converge (in addition to local uniform convergence of the continuous parts of the functions).Thus, D M [0, ∞) is a natural and -due to its properties -also a convenient space for convergence of discontinuous stochastic processes.
With reference to the notation used in Theorem 2.2, notice that the process has deterministic jump times and is therefore not Markov as a process in continuous time.In order to apply the theory of Markov processes, let (η t ) t≥0 be a Poisson process with rate ε −2 and define the pseudo-continuous process Z ε = (Z ε t ) t≥0 by Z ε t := X ε ηt for any t ≥ 0, which is Markov since the Poisson process is Markov.Then, the convergence of Z ε implies convergence of X ε to the same limit by the law of large numbers, see [22], Theorem 19.28.First we restate Theorem 2.2, using an equivalent formulation based on the process Z ε : Theorem 3.1.Consider the sequence of random variables (X ε i ) i∈N obtained by the construction from Algorithm 1.The process Z ε = X ε ηt t≥0 converges in distribution in the Skorokhod topology to the L-diffusion Z, i.e., the stochastic process with generator L defined in (2).This means that for any continuous and bounded function f : )) , ε → 0 .Certainly one cannot, in general, deduce convergence of stationary measures from time-local convergence statements about paths.However, standard arguments from the theory of Feller processes indeed allow to do so.This yields the following statement regarding stationary measures.
Theorem 3.2.Let µ ε be the stationary measure of the retraction-based random walk with stepsize ε.Then the weak limit of (µ ε ) ε>0 as ε → 0 is the stationary measure of the L-diffusion, with generator L defined in (2).In the case of second-order retractions, the limit is the Riemannian volume measure.
Remark 3.3.Notice that due to the approximation of the exponential map, the process Z ε is not in general reversible and the generator of Z ε is not in general self-adjoint.This is relevant since the stationary measure corresponds to the kernel of the adjoint of the generator.
Our proof of Theorem 3.1 and Theorem 3.2 (presented below) hinges on convergence of generators that describe the infinitesimal evolution of Markov processes, see, e.g., [17], Chapter 4. We show that the generator of Z ε converges to the generator L defined in (2).The generator of the process Z ε is spelled out in Lemma 3.4, and the convergence of this generator to L is treated in Lemma 3.5.
The following result is standard for transition kernels U ε on compact state spaces, see, e.g., [17], Chapter 8.3: Lemma 3.4.The process (Z ε t ) t≥0 is Feller, and its generator is given by , the continuous real valued functions on M .Here, and ω m is the volume of the unit (m − 1)-sphere.
Notice that Lemma 3.4 additionally states that the process Z ε satisfies the Feller property, i.e., that its semigroup is a contraction semigroup that also fulfills a right-continuity property.Therefore, in particular, the process Z ε is Markov.
Proof.Let P ε be the transition kernel of the geodesic random walk, i.e., Proposition 2.3 in [21] proves that the second summand on the right hand side converges to 0 as ε → 0. Let ε := √ m ε.Then Now consider the Taylor expansions of the functions f • Ret x and f • Exp x in ε at 0. To this end, notice that for any x ∈ M and v ∈ T x M , one has where v τ = d dτ Ret x (τ v).Hence, we obtain that where Hess f = ∇df is the Riemannian Hessian of f .Using that v 0 = v, the desired Taylor expansion takes the form , where for the last equality we used the fact that D Since Hess f is Hölder-continuous by assumption and Exp and Ret are smooth, both R Exp and R Ret are Hölder-continuous as products and sums of Hölder-continuous functions.Therefore, there exists a constant C > 0 such that one has the uniform bound Hence, which proves the lemma.□ Using the above results, we can finish the proof of Theorem 3.1.
Proof of Theorem 3.1.By the proof of Theorem 2.1 and the remark on page 38 in [21], the space C (2,α) (M ) is a core for the differential operator L. For a sequence of Feller processes (Z ε ) ε>0 with generators L ε , ε > 0 and another Feller process Z with generator L and core D the following are equivalent ( [22], Theorem 19.25): In our case D = C (2,α) (M ) and f ε = f ∈ D satisfy (i) due to Lemma 3.5.This concludes the proof of Theorem 3.1 since the L-diffusion is also a Feller process on a compact manifold.□ From an analytical point of view it is perhaps not surprising that convergence of generators, see Lemma 3.5, implies convergence of stationary measures as stated in Theorem 3.2.For the proof we here follow standard arguments of the theory of Feller processes.
Proof of Theorem 3.2.The family of (unique) stationary measures (µ ε ) ε>0 is a tight family of measures since all measures are supported on the same compact manifold M .By Prokhorov's theorem, which provides equivalence of tightness and relative compactness, any subsequence of the family (µ εn ) n∈N with ε n → 0 as n → ∞ has a convergent subsubsequence.But the uniform convergence of the generators that we showed in Lemma 3.5, see also (i) in the proof of Theorem 3.1, is also equivalent to the convergence of the associated semigroups with respect to the supremum norm in space, see again [22], Theorem 19.25.This implies ( [17], Chapter 4, Theorem 9.10) that all subsequential limits must be the unique stationary measure of the L-diffusion, and therefore all subsubsequences converge to the same measure.A standard subsubsequence argument then proves the theorem.□ Retractions on sub-Riemannian manifolds.A possible extension of our work concerns random walks in sub-Riemannian geometries.Indeed, sub-Riemannian structures can be used to model low-dimensional noise that lives in a higher dimensional state space [19].A sub-Riemannian structure on a smooth manifold M consists of a smoothly varying positive definite quadratic form on a sub-bundle E of the tangent bundle T M .Similar to the Riemannian setting, so called normal sub-Riemannian geodesics arise from a first order ODE on the cotangent bundle [28].Such geodesics are uniquely determined by an initial position x ∈ M and a 1-form α ∈ T * x M , and can therefore be approximated by a second order Taylor expansion of the respective ODE on the cotangent bundle.Such an approximation would algorithmically resemble our retraction-based approach, when attempting to efficiently simulate a sub-Riemannian geodesic random walk.However, in general, there exists no canonical choice for sampling a requisite 1-form α in every step in order to construct a sub-Riemannian geodesic random walk (or an approximation thereof).Indeed, one requires an additional choice of a measure for such sampling, and different choices may lead to different limiting processes, see [4,9].Notice that once such a choice has been made, one can seamlessly adapt our algorithm to the sub-Riemannian setting.We leave for future work a respective convergence analysis.Moreover, there might exist canonical choices of such measures for special sub-Riemannian structures -e.g., when considering the frame bundle of a smooth manifold M in order to model anisotropic Brownian motion and general diffusion processes [18,29].

Figure 1 .
Figure 1.Results for a parametric surface defined as ϕ(s, t) = ((1.1 + cos(s)) cos(t), (1.1 + cos(s)) sin(t), sin(s)).Left: Contour plot for the density of the stationary measure (i.e., the volume form of the Riemannian metric) in the parameter domain.Middle & right: Result of an approximate geodesic random walk based on p-Ret (using a single chart with periodic boundary conditions) with 100.000 steps and stepsize ε = 0.5, depicted in the parameter domain and in R 3 , respectively.Using an implementation in Mathematica on a standard laptop, the computation takes a few seconds.For a discussion of the relationship between the stationary measure and the empirical measure resulting from the random walk, see Section 2.3.

Figure 2 .
Figure 2. Result of a random walk using π-Ret for a genus two surface given as the zero-level set of f : R 3 → R, f (x, y, z) = (x 2 (1 − x 2 ) − y 2 ) 2 + z 2 − 0.01 with 100.000 steps and stepsize ε = 0.1.Using an implementation in Mathematica on a standard laptop, the computation takes less than a minute.
Convergence speed in different dimensions with fixed stepsize ε = 0.25.Convergence speed of walks with different stepsizes in dimension m = 3.

Figure 3 .
Figure 3. Convergence speed of retraction-based random walks to the invariant measures on m-dimensional tori parametrized by [0, 2π] m with periodic boundary conditions and equipped with the diagonal metric g ii = 1.5 + cos(x m−(i−1) ) for i = 1, . . ., m and g ij = 0 for i ̸ = j.Results are shown for varying dimension m and fixed stepsize ε (left) as well as fixed dimension m and varying stepsize ε (right).Simulations were performed 10 times for each combination of dimension and stepsize.The points display the arithmetic mean of the total variation distances between the invariant measure and the empirical measures resulting from a box count using 20 m equally sized cubes.