{'fraction_non_alphanumeric': 0.03428192947195322, 'fraction_numerical': 0.008838845240270418, 'mean_word_length': 4.600537221795856, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A hybrid projection algorithm for stochastic differential equations on manifolds

Stochastic differential equations projected onto manifolds occur widely in physics, chemistry, biology, engineering, nanotechnology and optimization theory. In some problems one can use an intrinsic coordinate system on the manifold, but this is often computationally impractical. Numerical projections are preferable in many cases. We derive an algorithm to solve these, using adiabatic elimination and a constraining potential. We also review earlier proposed algorithms. Our hybrid midpoint projection algorithm uses a midpoint projection on a tangent manifold, combined with a normal projection to satisfy the constraints. We show from numerical examples on spheroidal and hyperboloidal surfaces that this has greatly reduced errors compared to earlier methods using either a hybrid Euler with tangential and normal projections or purely tangential derivative methods. Our technique can handle multiple constraints. This allows, for example, the treatment of manifolds that embody several conserved quantities. The resulting algorithm is accurate, relatively simple to implement and efficient.


Introduction
Diffusion processes and stochastic trajectories on curved manifolds inside Cartesian spaces have many applications. The diffusion of atoms or molecules on many types of curved surface occurs in surface physics, biophysics, catalysis, biochemistry, cell biology, and nanotechnology [1,2,3,4,5,6,7,8]. Curved-space diffusion also arises in more abstract problems in general relativity [9], imaginary-time path-integrals [10] and in quantum field theory [11].
Other related examples include systems of stochastic equations with conserved quantities [12,13]. These are used for sampling purposes in studying finite temperature physical systems [14]. Multidimensional constraints are found in engineering systems [15,16], such as in robotics [17,18,19], where joints may have restricted mobility due to their mechanical design, and experience stochastic forces.
Phase-space stochastic processes [20,21,22] as well as quantum phase-space methods, for example representations [23,24] of Majorana physics [25,26] can also lead to curved space Fokker-Planck [27,28,29] and stochastic equations [30]. These are just a few of the many applications of such embedded stochastic processes occurring in technology and in science.
Projected diffusion processes of this type can sometimes be treated using stochastic equations for a lower dimensional set of coordinates on the manifold, for example by using spherical polar coordinates on a spherical surface. These can be termed implicit projections: the manifold is implicit in the definition of the coordinates themselves. The embedding theorems [31,32,33] of Riemannian geometry show that surfaces of arbitrary smooth Riemannian curvature are equivalent to an appropriate projection. The resulting stochastic equations usually involve multiplicative noise, because the noise coefficients depend on the state of the system. Several numerical algorithms are known for treating such multiplicative stochastic differential equations (SDE)s [34,35,36].
There are problems with implicit projections [37], especially for closed surfaces. The projected equations may have singularities in the noise and drift terms, due to the analyticity properties of maps [38] and the Poincare-Hopf theorem [39,40]. Stochastic numerical algorithms may not converge near singularities. Also, the definition of an intrinsic coordinate system [41,42] on curved surfaces can be a complex problem. There is generally no simple reduced coordinate system for a biological cell wall or other low-dimensional structures in biology and engineering, due to their complex geometry.
The approach used here is to define a stochastic process in the original coordinates, whose projection satisfies the restrictions of an embedded manifold. This overcomes the singularity and variable change issues. Even if the stochastic equation is known on the manifold, it can be useful to employ projections in order to eliminate singularities. The drawback is that discretization can lead to errors from finite step-sizes which may take the stochastic trajectory off the manifold. Such errors usually grow in time as the trajectory moves progressively away from the manifold. These additional constraint errors are not present in implicit approaches where the constraints are automatically satisfied.
In this paper a hybrid midpoint projection algorithm is introduced which combines both tangential and normal projections. A similar, but less accurate method has been used for pure diffusion on a two-dimensional manifold embedded in a three-dimensional space. This used an Euler method for each tangential step [43]. Here we generalize this approach to a more accurate midpoint projection algorithm for curved surface stochastic equations, which is well-defined for multiple constraints and arbitrary dimensions. The resulting algorithm treats drift as well as diffusion. Both numerical applications and examples are given.
We first review the theory of manifolds embedded in a Euclidean space. An explanation of how projections can occur physically is obtained through adiabatic elimination with constraint potentials. This is then applied to projections of stochastic processes. The continuous time limit of a constrained stochastic equation is derived as a Stratonovich equation, which is shown to be valid under conditions of finite bandwidth noise and fast projection timescales. Algorithms for multi-dimensional tangential and normal projections are derived. These methods are used to generalize some previous proposals for solving projected stochastic equations.
In order to compare the methods, numerical implementations are carried out. For simplicity, we treat diffusion on spheroidal and hyperboloidal surfaces in a three-dimensional Cartesian space. We use a stable implicit projected method with very small time-step to generate reference results. Simulations using hybrid Euler [43] and tangential methods [44] are also carried out, with errors calculated both for measures of distance diffused, and for conservation of constraints. We show that the hybrid midpoint projection algorithm has much lower errors in both categories than either of these earlier methods.
An outline of the paper is as follows. In Section 2, we review the geometric background of manifolds, constraints and projections. In Section 3 an adiabatic elimination approach is used to justify the interpretation of the projected equation as a Stratonovich equation. Section 4 gives a discussion of numerical algorithms and the different types of projections. Section 5 gives numerical examples and comparisons with intrinsic diffusion results for spheroidal and hyperboloidal surfaces embedded in Euclidean space. Section 6 summarizes the results. In the Appendix we explain the details of the intrinsic SDE variable changes and derivations.

Manifolds, constraints and projections
This paper treats numerical algorithms for projected stochastic equations. We utilize Stratonovich equations [45,30], which often occur in physics applications as the broad-band limit of a finite bandwidth random process. These follow standard variable change rules. Projections using the Ito stochastic calculus are also known [44], but are not treated here, and have different properties.
We start with a Stratonovich SDE on a Euclidean space R n , so that for a real vector trajectory x = x 1 , . . . x n , Here a (x, t) is a real vector function of x of dimension n, B (x, t) is a real n × s dimensional matrix function, where we omit the arguments when there is no ambiguity. The real s-dimensional Gaussian noise vector ξ (t) is assumed to have a finite correlation time T c , and is delta-correlated in the broad-band limit of small T c , so that: For a small time interval ∆t = t 1 − t 0 T c , we define a stochastic integral ∆w i as: which has the property that ∆w σ ∆w κ = δ σκ ∆t.
A single discrete step in the stochastic trajectory is ∆x = x (t 1 ) − x (t 0 ) . From the Stratonovich interpretation of stochastic calculus, the midpointx is the location where the derivatives are evaluated, wherē In the absence of a projection, the broad-band limit of Eq (1) is defined [45] as the limiting behavior of the following implicit stochastic difference equation In many applications it is useful to project this equation onto a manifold. The broad-band limit is only taken here after adiabatic elimination [46], and we assume that the time-scales of the constraint process are much faster than the noise correlation time, T c . Although we only treat real spaces here, these results are applicable to complex spaces, with minor changes. In this section we analyze how such manifolds are obtained, and the corresponding geometrical definitions.
Since we wish to derive algorithms numerically applicable to this problem, it is preferable to use a coordinatedependent notation rather than a coordinate-free approach. To avoid complex variable changes, the results will be expressed in the coordinates of the enclosing (Euclidean) space. In the Appendix, comparisons are made to intrinsic coordinate methods.

Manifolds
In the section we review the geometrical notation to define the manifold and its properties [47,48]. We equip R n with the Euclidean metric. We use standard contravariant notation for coordinates x = {x i }, and covariant notation for derivatives and normals n = {n i } to allow the results to be extended to other metrics. Repeated induces will be summed over in the Einstein convention. However, in the Euclidean case x i = x i so if the reader prefers, the index placement can be ignored.
An m-dimensional manifold is the set of points y ∈ M ⊂ R n for which where f is a p = n − m dimensional vector constraint function that defines the projection. We assume that the constraint equations are linearly independent.
Here f = f 1 , . . . f p denotes the p constraints. As examples, the constraints could be quadratic functions, where f i 0 is a real constant, h i is a real vector, G i is an N × N real matrix and: This covers the classical quadratic geometries, including circles, parabolas, spheres, hyper-spheres, ellipsoids, and hyperboloids, as well as intersections of these surfaces when there are multiple constraints. We use these for examples, and choose cases such that f 0 = −1 and h = 0. The general algorithms discussed are not limited to these quadratic constraints, which are introduced because these manifolds have well-known intrinsic coordinates and diffusion equations for comparison purposes.. In some cases it is possible to define φ ∈ R m as a new intrinsic coordinate system for the set of points restricted to the manifold. For example, on the surface of a sphere one may use polar coordinates. In such cases, a vector mapping function Φ exists where φ = Φ (y) for every Euclidean coordinate y in the original manifold M.
Such mapping functions lead to singularities for a compact manifold, and may not be able to be expressed in a closed form. As a result, we use algorithms that can be expressed in the original Euclidean coordinates. Intrinsic coordinates are treated in the Appendices as special cases. These are used for numerical examples, although such comparisons are only available in special cases where analytic transforms exist.

Tangent and normal spaces
Near to a point y in the manifold, from a Taylor expansion of the j-th constraint equation, where ∆ i = x i − y i , ∂ i ≡ ∂/∂x i , and the constraint derivative is From this condition, v j ⊥i = ∂ i f j is called the j-th normal vector, v j ⊥ , where the lower index indicates a covariant vector, and the upper index labels the constraint, f j , so that j = 1, . . . p. Any orthogonal vector v such that v · v j ⊥ = 0 locally satisfies the j−th constraint, although it may not satisfy the other constraints. In general the normal vectors v j ⊥ may not be orthogonal to each other. They are assumed to be linearly independent at all locations in order to generate a differentiable manifold locally isomorphic to R m without a singular point [47]. Singular points may occur with some choices of multiple projections, but these are unphysical and can lead to convergence issues. One can use the gradient vectors of the constraints to define a set of orthonormal vectors n j that span the p-dimensional space of vectors normal to each constraint. This requires orthogonalization of the set v j ⊥ , either by Gram-Schmidt or other methods [49]. As a result, the orthonormal vectors n j have the property that, for i, j = 1, . . . p: Similarly, there is a set of orthonormal tangent vectors m i , where i = 1, . . . n. These define the tangent space of the manifold, and give an orthonormal basis of m contravariant vectors orthogonal to all of the normal vectors v j ⊥ , so that, for j = 1, . . . p and i = 1, . . . m: One can define a vector space consisting of all vectors tangent to M at any point y ∈ M. This is the tangent space, T y M, and it is spanned by the basis m j for j = 1, . . . m. The orthogonal complement to the tangent space is the normal space, spanned by the basis n j for j = 1, . . . p.
For every vector v ∈ R N and manifold coordinate y ∈ R m , a vector decomposition can therefore be written as follows [50]: provided v , v ⊥ are the tangential and normal projection of v ∈ R m respectively and T ⊥ y M is the orthogonal complement of T y M. If we consider the sphere as an example, the vectors perpendicular to the radii are tangents, while the radial vectors are in the normal space.
This can be expressed in terms of projection operators P ⊥ and P y . These satisfy the fundamental property of mathematical projections on sets that, once projected, a coordinate doesn't change under further projections, so P (P (x)) = P (x).

Tangential projections
Eq. (12) guarantees that any tangent vector of form y + m i is on the tangent manifold. On this manifold, in a neighborhood of a point y ∈ M , all constraints are satisfied up to terms of order 2 since for tangent vectors m i and a small coefficient : A general tangential projection of any time-evolution equation is obtained by projecting an arbitrary derivative vector δ at each location onto the tangential plane of the manifold. This type of projection depends on the location, so we use the notation of δ = Π y (δ), for a tangential projection of a derivative vector δ at a point on the manifold.
One can obtain the tangential projection δ in a simple case of a one-dimensional constraint, by removing the normal component from the vector δ using Eq. (13), so that: Here n ⊗ n is an outer vector product defined as the matrix [n ⊗ n] ij = n i n j , n is the normal vector at y, I is the identity matrix, and (I − n ⊗ n) is called the tangential operator [51]. More generally, for multiple constraints, the projection must remove multiple normal components n j , so one must define:

Normal projections
In a normal projection, points in a neighborhood of M are projected onto the nearest point on M. The projection y = P ⊥ (x) therefore takes an initial x ∈ R n , and maps it onto the nearest manifold coordinate y.
For example, in the spherical case, a path in Euclidean coordinates x in n-space is projected onto the n − 1 dimensional hyper-spherical surface of a ball of radius r around the origin, at every step on the path. An exact projection function in this case is: which is defined for all x = 0. Exact normal projections are only obtainable in special cases. Hence, we consider approximate normal projections valid near a manifold point.
To obtain the general approach which solves this, we introduce a Lagrange multiplier vector c. Given an initial point x and a projected point y on the manifold, one must find Lagrange multipliers c j (x, y) such that: Solving for y corresponds to minimizing the functional k(y) = |y − x| 2 , subject to the constraint that f i (y) = 0. This is generally a hard computational problem, but can be solved approximately with a Taylor expansion, valid if |y − x| is small. Since f j (y) = 0, therefore to first order in a Taylor series about y, one finds that This follows since the v ⊥ j vectors are all derivatives. From Eq (18), One must therefore solve the following nonlinear implicit equation: In this equation, the normal projection matrix M is defined so that its components are the the inner products of v i ⊥ (y) with the normal vectors n j (y), i.e, Solving for c j , one obtains an implicit equation: Given an initial estimate x, a normal projection is obtained by solving Eq. (23) for c, and hence obtaining y using iteration, where the matrix M (x) is used instead of M (y) to lowest order. In the algorithms described here, a single iterative step is used because the trajectories are already close to the manifold, so that: In the case of a one-dimensional constraint, since n ⊥ (x) = 1, it follows that v ⊥ (x) · n (x) = v ⊥ (x) . For small |y − x|, one therefore obtains:

Adiabatic elimination for projected equations
To understand the projection of a stochastic equation, one must ask: what physical or mathematical process restricts the path to a manifold? To prevent ambiguities, we treat a constraining potential, as a common situation found in physics, and also used for quantum applications where there are similar issues [48]. While this is not universally applicable for all projected stochastic equations, it is well-defined as a method for explaining our choice of stochastic calculus.
We introduce a local mapping function Φ inside a local patch where φ = Φ (y) for Euclidean coordinates y in the original manifold M. For purposes of adiabatic elimination, we extended this locally to an invertible mapping x = x (φ) including nearby off-manifold points, with a complete set of coordinates Specifically, we define φ ⊥ to correspond to points in the normal space, where if y is a manifold point, then for nearby points x in Euclidean space, φ ⊥ is the distance along a normal vector: To return to the spherical example, the coordinate φ ⊥ for overall dimension n = 3 is the radial coordinate in that case.

Constraining potentials
Constraining potentials are essential for constraining particles on a manifold, and provide a model for the origin of a projected stochastic equation. We assume that the constraint conditions originate in a scalar constraint potential The projected stochastic equation is the adiabatic limit for λ → ∞, of a continuous stochastic process in which the initial SDE in Eq. (1) includes the constraint as a potential, so that: Here i is the contravariant derivative, which equals the covariant derivative ∂ i u for our metric choice. The drift for motion including the confining potential is then a λ (x) = a − ∇u (x). The gradient of the potential is, In terms of the local coordinates φ, the constrained equation is: where the Jacobian J is defined as: From the definition of the normal coordinates in Eq (26), this implies that for i > m, where m is the manifold dimension, and x ≈ y, the Jacobian is equal to the corresponding normal vector: The reduction of the number of independent degrees of freedom that result from a minimization of the potential u is equivalent to adiabatic elimination of the fast transverse variables in a stochastic equation, treated next.

Adiabatic elimination
Adiabatic elimination will be used in order to analyze how projections are applied to stochastic equations. As has been recently pointed out [44], there can be ambiguity in defining the which type of stochastic equation is obtained after projection. This depends on how the projection occurs. We assume that the original SDE in Eq. (1) has noise with a finite bandwidth. Hence, it follows the Stratonovich calculus in the broad-band limit [45].
The projection is assumed to occur through an adiabatic elimination process, due to the constraint potential. The resulting diffusion or noise is generally state-dependent, even if not originally. To analyze this, we follow a similar method to the work of Gardiner [46]. In general, a direct approach to adiabatic elimination is obtained through dividing up the variables into a "fast" and "slow" set. In the present context, these are locally in the direction of the normal and tangential coordinates.
We take a given pointφ on the manifold such that f x φ = 0. We now consider a trajectory φ (t) nearφ, governed by the constrained stochastic equation in Eq. (30), so that Changes in the coordinate, ∆, are divided into "slow" coordinates ∆ for motion inside the manifold, and "fast" or normal coordinates ∆ ⊥ for motion outside the manifold. For a small displacement we obtain: In a neighborhood of the manifold where x ≈ y, these have equations given by the chain rule for variable changes, noting that we only consider smooth differentiable functions here:

Fast variables
Expanding the displacement ∆ according to Eq. (34), we define Here a i ⊥ is the i-th component of the normal drift, and the inner-product notation refers to inner products in the original Euclidean space. Each term b i ⊥ is a covariant vector in the s-dimensional space of noise terms: This gives a locally valid stochastic equation for the fast variables, where the inner product b i Here we define the constraint linear response matrix R as One can rewrite R in a factored form as R = M T M , where M ij = v i ⊥ · n j is the normal projection matrix of (22). From Eq (11), the diagonal elements of R are non-vanishing. As a result, R is a symmetric, positive definite matrix, provided the constraints are non-singular. Since it must have positive eigenvalues, the adiabatic limit is found by taking the limit of λ → ∞. This is equivalent to setting∆ ⊥ j = 0, due to the rapid equilibration that occurs.
We also assume that λT c 0 , with the result that the "fast" set ∆ ⊥ experience a large rate of change defined by λ, and relaxes to equilibrium on a short time scale of 1/λ. Therefore, in the adiabatic limit, ∆ ⊥ = 0, and the system is constrained to the tangent manifold.

Slow variables
By contrast, the slow set ∆ does not relax at the same rate. The constraint terms vanish, since the Jacobian for these variables is orthogonal to the constraint terms. Substituting the resulting values into the "slow" equations therefore leads to a simpler projected equation with fewer independent variables. We define, for µ = 1, . . . m : This equation can be written as an equation purely in terms of the intrinsic manifold coordinates, as: with the definitions that: This leads to an important question. Is the projected SDE an Ito or Stratonovich equation? Our approach is to assume that the original equation prior to projection has a finite bandwidth. This means that ordinary calculus rules are applicable at each stage. The broadband limit is then taken after adiabatic elimination, subject to the restriction that λT c 0, as explained above. This leads to a Stratonovich interpretation. The intrinsic equation can also be rewritten in a shorthand form as: Gardiner [46] demonstrates that in two-dimensional cases, then even if the original equation is in the broad-band limit, adiabatic elimination will lead to the Stratonovich interpretation of the resulting SDE, if there are no fast variables in the noise coefficients. More generally, for broad-band noise equations, one should make a case-by-case analysis to determine if there are stochastic correction terms. The constraint equations are now automatically satisfied, provided they are satisfied initially, sincė While these equations are correct in the continuous limit, it is important to take additional steps to ensure the solutions remain on the manifold. To explain this, if we use the midpoint definition of Eq. (5), the corresponding discrete midpoint algorithm in Cartesian coordinates can be written as: However, due to discretization error, the path will not in general remain on the manifold with this discretization, as we show in numerical examples below. Since our derivation requires that the constraints are satisfied exactly, it is necessary to make an additional normal projection, to satisfy the constraints, so that the final algorithm reads:

Spherical example
As an illustration, consider a hyper-spherical manifold with |x| 2 = 1 and isotropic diffusion with B = I in the original Euclidean space. Writing Eq. (1) as a difference equation as in Eq. (6) leads to ∆x = ∆w.
where the Gaussian noise ∆w integrated over an interval ∆t is correlated according to Eq. (4. Simply projecting this according to the constraint naively leads to a projected equation: where P ⊥ (x) projects the new vector y = x + δ normally onto the sphere, or more generally onto a manifold M. Instead of this approach, using the adiabatic elimination combined with a potential, one might choose that: On the projected manifold, the restoring drift towards the surface is normal, since: For a tangential projection P x on the manifold, the resulting stochastic equation is then: From the results in the previous subsection, we regard this as a Stratonovich equation [45], equivalent to evaluating the drift and projections at the midpoint, and to maintain the path on the manifold the algorithm is modified to:

Stochastic projection algorithms
We first summarize the numerical algorithms for tangential and normal projections, and then explain how these are used in projected stochastic equations. In some cases [44], a tangential projection of the coefficients of a Stratonovich SDE is used. This doesn't guarantee the solution is on the manifold. The difficulty with purely tangential methods is that since errors can accumulate at finite step-size, the final result can move arbitrarily far way from the desired manifold. Other proposed methods that use restrictions to a manifold are often limited to two-dimensional surfaces in three-dimensional spaces [52].
Another approach applied to physics problems [43], is to employ a tangentially projected Ito-Euler step followed by normal projections to remain on the manifold. Since this is a forward Euler method, the algorithm can have large step size errors. It has the advantage that the final normal projection keeps the solution on the manifold. However, Ito stochastic equations are not generally equivalent to Stratonovich stochastic equations, unless appropriate corrections are employed. Therefore, this algorithm may not be generally applicable for an arbitrary projection. We give comparative results for simple cases in the next section.
The hybrid midpoint projection algorithm described here combines a midpoint method [34] instead of an Euler step, together with normal projection. The manifold is given by the constraint equations f i (x) = 0, wherex is defined as a midpoint for each tangential step. We assume the existence of a stable orthogonalization algorithm O (S) that generates an orthonormal set of vectors , S ⊥ = n 1 , . . . n p , from a set of gradient vectors S = v j ⊥ , . . . v p ⊥ , using one of the known standard techniques [49,53].
The three algorithms compared here all use the notation that x 0 , x 1 are the initial and final locations of a step in time, ∆t is the step-size, and ∆w are random Gaussian variables. They have correlations that correspond to a discretized delta-function in time, as described in Eq (4).
In the following section, we give examples of the use of these algorithms by comparing them with well-converged, high-accuracy simulations using intrinsic coordinates, for three-dimensional projections onto two-dimensional ellipsoidal and hyperboloidal surfaces.
We focus on low order, high performance methods, as they are often most useful in applications of stochastic equations. The reason for this is that the total error combines discretization and sampling errors. Reducing the discretization error using high-order methods is not as useful as with ordinary differential equations, since the resulting time penalty may lead to fewer samples, increased sampling error and little improvement.

Hybrid Euler projection algorithm (hEP)
Holyst [43], describes a hybrid Euler projected method with pure diffusion. The original proposal had no drift vector, a = 0, and was restricted to dimensions n = 3 and m = 2. Their approach combined an Ito-Euler step in time, with a tangential projection of the noise at the initial point, and an approximate normal projection at the end of the step in time. Here we describe this approach in a slightly more general form by including a drift term and a constant noise matrix B, since an Ito-type algorithm would require additional corrections if the diffusion was space-dependent.
In the notation of Eq ( 6), this is defined for a step starting at x 0 , assuming that a is the drift and B is the noise matrix, as: It is described in detail as an algorithm in the box below, 1. Evaluate a step, δ = a (x 0 ) ∆t + B · ∆w.
2. Obtain an orthonormal set n j at x 0 . 3. Tangentially project: δ = δ − m j=1 n j δ · n j 4. Estimate the next point: x = x 0 + δ 5. Obtain normal vectors v j ⊥ and an orthonormal set n j at x .
Although approximately correct in the original scenario, this algorithm is not recommended here, since it has much greater errors than hybrid midpoint methods, as shown in the next section. It also may not converge accurately in general, since we are fundamentally considering projection as a Stratonovich process, and Ito corrections may be necessary.

Midpoint projection algorithm (MP)
We next consider purely tangential projections using a Stratonovich midpoint method, valid for small changes in x, although this can lead to error-propagation issues. The midpoint projection algorithm can be expressed symbolically as: They do not guarantee that a path stays on the manifold except in the limit. Hence it should be supplemented by another check, explained next.
In some cases this has lower errors than the Euler method, but it is not constrained to say on the manifold. As shown in the numerical examples in the next section, this leads to an increasing constraint error with time.

Hybrid midpoint projection algorithm (hMP)
One may combine these two methods to obtain a hybrid midpoint projection algorithm for the projected stochastic differential equation which has both good accuracy and long-term stability. It can be written symbolically as: This is described in the box below, where we assume that initially i = 1, and x (1) = x 0 .
This uses a tangential projection for each midpoint iteration, together with a final normal projection to ensure that the resulting step remains on the manifold. The reason for this choice is that the stochastic departure from the manifold can be relatively significant if allowed to propagate. The final normal projection ensures that this additional error is removed.
It is important to ensure that the projected path remains on the manifold, since error-propagation may cause a long term drift off the manifold. Hence, there are two parts to the calculation, which in principle would require different numbers of iterations. In the examples we use three iterations for the first part, and find it is sufficient to choose a single nonlinear iteration for the final projection.

Numerical examples
In this section we treat numerical examples corresponding to the different SDE projection algorithms. These tests used a public domain stochastic toolbox, xSPDE [54] to compare three different curved space diffusion SDE projection methods with an intrinsic method, using spheroidal and hyperboloidal surfaces embedded in a threedimensional manifold as simple examples for the comparative studies. (4.1,4.2&4.3).

Computational details
We consider the diffusion on spheroidal and hyperboloidal surfaces for testing purposes. These are useful as a test case, since they have well-studied intrinsic SDEs that can be used for comparisons. These examples have positive and negative curvatures which vary in space, so they cover a wide variety of conditions.
The plots below correspond to: • the average great circle distance, Θ = cos −1 x T Gx 0 , • the mean Euclidean distance, d = |x − x 0 | , • great circle distance error, Θ , • Euclidean distance error, d , and • the error in the manifold constraint, |f | .
Errors are calculated by comparison with intrinsic simulations. Details of the intrinsic Stratonovich equations used for the calculations are given in the Appendix. The intrinsic data has sampling errors and step size errors. Errors in these reference calculations are reduced as follows: • The sampling errors are reduced by using 2 × 10 6 trajectories in both cases, which leads to a typical sampling error of ∼ 10 −3 in the distances. • Step-sizes of 10 −3 were used for the intrinsic method. A stable implicit midpoint algorithm was employed [34].
Step-size errors were estimated by comparison with a second calculation with a step-sizes of 0.5 × 10 −3 , showing that step-size errors were negligible.
The steps in the cases used for the comparative study of projected algorithms were larger in order to investigate the size of the discretization error, with a step-size of 10 −2 . In all cases, there were 2 × 10 6 random trajectories,

Ellipsoid
An ellipsoid is defined by the constraint equation Setting a = b = 1 yields a spheroid, which is a distortion of the unit sphere along one axis. The metric tensor on the spheroid simplifies greatly compared to the general ellipsoid, and hence we restrict ourselves to this shape. In order to test the different projection algorithms outlined previously, we compare the results of the projected SDE to a direct simulation of diffusion on the spheroid. The equations used to simulate the latter are derived in the Appendix, and are: Projecting a point in R 3 on the spheroid is obtained as where G = diag 1, 1, 1/c 2 . This projection holds generally for quadratic constraints, provided the matrix G is replaced by the matrix appropriate for the manifold under consideration. However, we use the more generally applicable local normal projections in the examples plotted below. The initial point in the simulations was at θ = φ =1, so that it was a large distance from any singularity. This was allowed to diffuse following the spheroidal diffusion equations for a time-interval of T = 1. The results are shown in Figs (1, 2, 3 & 4), where c = 0.25 was chosen. The measure of distance used for these figures is The results show that the hMP algorithm gives excellent agreement between the intrinsic and the projected SDE results.

Hybrid Euler projection algorithm (hEP)
Because of the final projection employed, this algorithm shows a typical behavior of poor convergence with step-size, but good ability to maintain the trajectory on the manifold. The maximum discretization error in the mean great circle distance Θ is 0.07, and the maximum average constraint error, |f | is 0.009. While the errors are reasonably under control, convergence is rather poor for the number of integration steps used.

Midpoint projection algorithm
Because of the midpoint tangential projection employed, this algorithm shows improved convergence with stepsize, but is unable maintain the trajectory on the manifold. The step-size or discretization errors in the average great circle distance is 0.034±.001, but the average constraint error is 0.19. This demonstrates that while the distance error appears low, this is misleading. Convergence as measured by the constraint error is extremely poor for the number of integration steps used, and these trajectories have moved a substantial distance from the manifold.

Hybrid midpoint projection algorithm
Because of the hybrid midpoint tangential projection employed, this algorithm shows greatly improved convergence with step-size, and is able to maintain the trajectory on the manifold. In this case the maximum step-size or discretization error in the average great circle distance is 4 ± 0.5 × 10 −3 , and the maximum average constraint error is extremely small, at 7 × 10 −4 . This shows excellent convergence, with well-controlled distance error combined with a greatly reduced constraint error compared to either of the other two methods.

Hyperboloid
The final example we consider is that of the one-sheeted hyperboloid, defined by the constraint equation We choose a = b = 1 and c = 0.25 to simplify the metric tensor, which is required to obtain the intrinsic SDEs for diffusion on this manifold. Details are given in the Appendix.

The corresponding intrinsic equations arė
The projection used for this manifold is simply given by Eq (56), with G = diag 1, 1, −1/c 2 . Applying the three different projection algorithms for SDEs yields Figs (5, 6, 7 & 8), where we use the Euclidean distance [55] measure, The initial point in the simulations was at θ = v =1. This was allowed to diffuse following the spheroidal diffusion equations for a time-interval of T = 1. Similarly to the ellipsoid, the results show almost exact agreement between the intrinsic and hybrid midpoint projected SDE algorithms, with poor convergence in the other two cases.

Hybrid Euler projection algorithm
As in the ellipsoidal case, this algorithm shows a typical behavior of the algorithm: it has poor convergence with step-size, but good ability to maintain the trajectory on the manifold, due to the use of a normal projection. The maximum discretization error in the average Euclidean distance is 0.037, and the maximum average constraint error is 0.025. These errors are relatively large given the small step-size used here.

Midpoint projection algorithm
Again similar to the ellipsoidal case, this algorithm shows improved convergence with step-size, but is unable maintain the trajectory on the manifold. The maximum step-size or discretization error in the average Euclidean distance is 3.5 × 10 −3 , but as before, this is misleading. In fact, the trajectories have moved well off the manifold, since the maximum average constraint error is 0.15.

Hybrid midpoint projection algorithm
This algorithm repeats the trend in the ellipsoidal case, with both improved convergence with step-size, and much greater ability to maintain the trajectory on the manifold. In this case, the maximum step-size or discretization error in the mean Euclidean distance is 4.7 × 10 −3 , and the maximum average constraint error is 2.7 × 10 −3 .

Summary of results
In summary, we have shown how projected SDEs are derived using adiabatic elimination with a constraint potential, and shown that the corresponding projected equation is a Stratonovich type. We have obtained a hybrid midpoint projection algorithm for projecting a stochastic equation onto a general manifold, and compared it to earlier proposals, using numerical studies on ellipsoidal and hyperboloidal surfaces of positive and negative curvature.
A comparison of all six cases studied is shown the  The hybrid midpoint algorithm clearly gives a distinctly improved error performance compared to both the earlier proposals. The improvement is due to a combination of tangential and normal projections, together with an intrinsically accurate midpoint algorithm. The cause of this is clear. The hybrid Euler projection method proposed earlier cannot accurately track the changes in the projected diffusion matrix during a step, as it relies on an initial estimate of the diffusion, which changes in space due to the projection. The midpoint projection method has an improved treatment of diffusion, but is unable to control error growth as trajectories steadily move off the manifold due to error accumulation.
The hybrid midpoint projection method, by correcting for the changes in the diffusion matrix and maintaining the constraint in time, has much greater accuracy than either of these. As other alternatives, adaptive step-size [56] or higher order methods are also possible [57,58,35]. Higher-order algorithms have been recently proposed for projected equations [59,60] which may give, potentially, even lower time-step errors.
However, there is an important trade-off, discussed elsewhere [61,62]. The complexity of variable step-size or higher order methods means that they are often difficult to parallelize. They also are usually slower per time-step. This gives fewer parallel trajectories, with less samples in the same time. The result can be an increase in the sampling error, and an increased total error for the same computational time. This trade-off is essentially caused by the need to compute many stochastic trajectories.
By comparison, the hybrid midpoint projected method combines excellent error performance with low complexity, straightforward parallel implementation and relatively high speed. This is a great advantage when it is important to also reduce sampling errors caused by the stochastic noise.

Intrinsic SDE derivations
The projected coordinates y = P (x) in the main text are coordinates on an embedded manifold in Euclidean space R n . A further transformation is often made to obtain intrinsic coordinates on the manifold with lower dimensionality. Hence, to obtain φ, we define One can make subsequent transformations on φ to obtain other systems of coordinates on the manifold, such as defining spherical polar coordinates with rotated polar directions, but these are all isomorphic to the set of projected coordinates {P (x)}. For projected points, y = P (x), and for the examples in the main text, we take the simplest case that the original manifold has a locally Euclidean metric and Euclidean diffusion g E ij = δ ij on a tangent plane. The final mapping then defines the induced metric g, which is in general a curved metric [33], such that: We therefore consider an m-dimensional manifold equipped with a metric g µν , in an intrinsic coordinate system with coordinates φ. As an example, the Fokker-Planck equation governing non-driven Euclidean diffusion on this manifold is given by where ∆ = ∇ µ ∇ µ is the Laplace-Beltrami operator,P is the covariant probability density, and D r is the covariant diffusion coefficient. This FPE is in a covariant form, meaning that probability is conserved with respect to the measure √ gd m φ, where g = det (g µυ ) . If one scales the covariant probability density by √ g, the result is a probability density P that is conserved with respect to d m φ.
Inserting also the formula [63,65] , one finds that ∂P ∂t with g µν = [g µν ] −1 , and P = √ gP . Equivalently, using Einstein summation convention for repeated indices, ∂P ∂t = −∂ µ 1 4 g µν ∂ ν ln g P + 1 2 ∂ µ g µν ∂ ν P, where we set D r = 1. In this expression, where g µν = σ β µ σ β ν σ . the self-adjoint drift is This FPE may be transformed to the Stratonovich form, from which one easily reads off an equivalent system of SDEs, by applying the formula The corresponding projected Stratonovich SDE is theṅ However, the fact that the final diffusion matrix is equal to the induced metric is due to the form of the initial diffusion matrix, which was chosen equal to the Euclidean metric tensor for this example. This is not always the case, and the final diffusion matrix has to be worked out accordingly.

Spheroidal diffusion
As an example in the main text, a 2-spheroid may be embedded in R 3 by the constraint equation x 2 + y 2 + z 2 /c 2 = 1.

Hyperbolic diffusion
As another example used in the main text, the one-sheeted hyperboloid is defined by the constraint x 2 + y 2 − z 2 /c 2 = 1.
A convenient set of intrinsic coordinates is given by where −∞ < v < ∞ and 0 ≤ θ < π. This yields the inverse metric and making using of the above equation and the Eqs (71&72) one obtains the intrinsic Stratonovich SDE as follows: v = tanh v c 2 − 1 + (c 2 + 1) cosh 2v + ξ v