Stochastic Lie group integrators
STOCHASTIC LIE GROUP INTEGRATORS
SIMON J.A. MALHAM∗ AND ANKE WIESE∗
Abstract. We present Lie group integrators for nonlinear stochastic differential equations with
non-commutative vector fields whose solution evolves on a smooth finite dimensional manifold. Given
a Lie group action that generates transport along the manifold, we pull back the stochastic flow on
the manifold to the Lie group via the action, and subsequently pull back the flow to the corresponding
Lie algebra via the exponential map. We construct an approximation to the stochastic flow in the
Lie algebra via closed operations and then push back to the Lie group and then to the manifold,
thus ensuring our approximation lies in the manifold. We call such schemes stochastic Munthe-Kaas
methods after their deterministic counterparts. We also present stochastic Lie group integration
schemes based on Castell–Gaines methods. These involve using an underlying ordinary differential
integrator to approximate the flow generated by a truncated stochastic exponential Lie series. They
become stochastic Lie group integrator schemes if we use Munthe-Kaas methods as the underlying
ordinary differential integrator. Further, we show that some Castell–Gaines methods are uniformly
more accurate than the corresponding stochastic Taylor schemes. Lastly we demonstrate our methods
by simulating the dynamics of a free rigid body such as a satellite and an autonomous underwater
vehicle both perturbed by two independent multiplicative stochastic noise processes.
Key words. stochastic Lie group integrators, stochastic differential equations on manifolds
AMS subject classifications. 60H10, 60H35, 93E20
1. Introduction. We are interested in designing Lie group numerical schemes
for the strong approximation of nonlinear Stratonovich stochastic differential equa-
tions of the form
yt = y0 +
Vi(yτ , τ) dW
τ . (1.1)
HereW 1, . . . ,W d are d independent scalar Wiener processes andW 0t ≡ t. We suppose
that the solution y evolves on a smooth n-dimensional submanifold M of RN with
n ≤ N and Vi : M × R+ → TM, i = 0, 1, . . . , d, are smooth vector fields which in
local coordinates are Vi =
j=1 V
i ∂yj . The flow-map ϕt : M → M of the integral
equation (1.1) is defined as the map taking the initial data y0 to the solution yt at
time t, i.e. yt = ϕt ◦ y0.
Our goal in this paper is to show how the Lie group integration methods developed
by Munthe-Kaas and co-authors can be extended to stochastic differential equations
on smooth manifolds (see Crouch and Grossman [8] and Munthe-Kaas [40]). Suppose
we know that the exact solution of a given system of stochastic differential equations
evolves on a smooth manifold M (see Malliavin [36] or Emery [14]), but we can only
find the solution pathwise numerically. How can we ensure that our approximate
numerical solution also lies in the manifold?
Suppose we are given a finite dimensional Lie group G and Lie group action Λy0
that generates transport across the manifold M from the starting point y0 ∈ M via
elements of G. Then with any given elements ξ in the Lie algebra g corresponding to
the Lie group G, we can associate the infinitesimal action λξ using the Lie group action
Λy0 . The map ξ 7→ λξ is a Lie algebra homomorphism from g to X(M), the Lie algebra
∗Maxwell Institute for Mathematical Sciences and School of Mathematical and Computer Sciences,
Heriot-Watt University, Edinburgh EH14 4AS, UK. (S.J.Malham@ma.hw.ac.uk, A.Wiese@hw.ac.uk).
(16/10/2007)
http://arxiv.org/abs/0704.0022v2
2 Malham and Wiese
of vector fields over the manifold M. Further the Lie subalgebra {λξ ∈ X(M) : ξ ∈ g}
is isomorphic to a finite dimensional Lie algebra with the same structure constants
(see Olver [42], p. 56).
Conversely, suppose we know that the Lie algebra generated by the set of govern-
ing vector fields Vi, i = 0, 1, . . . , d, on M is finite dimensional, call this XF (M). Then
we know there exists a finite dimensional Lie group G whose Lie algebra g has the
same structure constants as XF (M) relative to some basis, and there is a Lie group
action Λy0 such that Vi = λξi , i = 0, 1, . . . , d, for some ξi ∈ g (see Olver [42], p. 56 or
Kunita [30], p. 194). The choice of group and action is not unique.
In this paper we assume that there is a finite dimensional Lie group G and action
Λy0 such that our set of governing vector fields Vi, i = 0, 1, . . . , d, are each infinitesimal
Lie group actions generated by some element in g via Λy0 , i.e. Vi = λξi for some ξi ∈ g,
i = 0, 1, . . . , d. They are said to be fundamental vector fields. This means that we can
write down the set of governing vector fields Xξi for a system of stochastic differential
equations on the Lie group G that, via the Lie group action Λy0, generates the flow
governed by the set of vector fields Vi on the manifold. The vector fields Vi on M are
simply the push forward of the vector fields Xξi on G via the Lie group action Λy0 .
Typically the flow on the Lie group also needs to be computed numerically. We thus
want the approximation to remain in the Lie group so that the Lie group action takes
us back to the manifold.
To achieve this, we pull back the set of governing vector fields Xξi on G to the set
of governing vector fields vξi on g, via the exponential map ‘exp’ from g to G. Thus
the stochastic flow generated on g by the vector fields vξi generates the stochastic flow
on G generated by the Xξi . The set of governing vector fields on g are for each σ ∈ g:
vξi ◦ σ ≡
(adσ)
k ◦ ξi , (1.2)
where Bk is the kth Bernoulli number and the adjoint operator adσ is a closed operator
on g, in fact adσ ◦ ζ = [σ, ζ], the Lie bracket on g. Now the essential point is that
ξi ∈ g and so the series on the right or any truncation of it is closed in g. Hence
if we construct an approximation to our stochastic differential equation on g using
the vector fields vξi or an approximation of them achieved by truncating the series
representation, then that approximation must reside in the Lie algebra g. We can then
push the approximation in the Lie algebra forward onto the Lie group and then onto
the manifold. Provided we compute the exponential map and action appropriately, our
approximate solution lies in the manifold (to within machine accuracy). In summary,
for a given ξ ∈ g and any y0 ∈ M we have the following commutative diagram:
∗−−−−→ X(G)
(Λy0 )∗−−−−→ X(M)
exp−−−−→ G
Λy0−−−−→ M
We have implicitly separated the governing set of vector fields Vi, i = 0, 1, . . . , d,
from the driving path process w ≡ (W 1, . . . ,W d). Together they generate the unique
solution process y ∈ M to the stochastic differential equation (1.1). When there
is only one driving Wiener process (d = 1) the Itô map w 7→ y is continuous in
the topology of uniform convergence. When there are two or more driving processes
Stochastic Lie group integrators 3
(d ≥ 2) the Universal Limit Theorem tells us that the Itô map w 7→ y is continuous
in the p-variation topology, in particular for 2 ≤ p < 3 (see Lyons [32], Lyons and
Qian [33] and Malliavin [36]). A Wiener path with d ≥ 2 has finite p-variation for
p > 2. This means that from a pathwise perspective, approximations to y constructed
using successively refined approximations to w are only guaranteed to converge to
the correct solution y, if we include information about the Lévy chordal areas of
the driving path process. Note however that the L2-norm of the 2-variation of a
Wiener process is finite. In the Lie group integration procedure prescribed above we
must solve a stochastic differential system on the Lie algebra g defined by the set
of governing vector fields vξi and the driving path process w ≡ (W 1, . . . ,W d). In
light of the Universal Limit Theorem and with stepsize adaptivity in mind in future
(see Gaines and Lyons [20]), we for instance use in our examples order 1 stochastic
numerical methods—that include the Lévy chordal area—to solve for the flow on the
Lie algebra g.
We have thus explained the idea behind Munthe-Kaas methods and how they can
be generalized to the stochastic setting. The first half of this paper formalizes this
procedure.
In the second half of this paper, we consider autonomous vector fields and con-
struct stochastic Lie group integration schemes using Castell–Gaines methods. This
approach proceeds as follows. We truncate the stochastic exponential Lie series expan-
sion corresponding to the flow ϕt of the solution process y to the stochastic differential
equation (1.1). We then approximate the driving path process w ≡ (W 1, . . . ,W d) by
replacing it by a suitable nearby piecewise smooth path in the appropriate variation
topology. An approximation to the solution yt requires the exponentiation of the
approximate truncated exponential Lie series. This can be achieved by solving the
system of ordinary differential equations driven by the vector field that is the approx-
imate truncated exponential Lie series. If we use ordinary Munthe-Kaas methods as
the underlying ordinary differential integrator the Castell–Gaines method becomes a
stochastic Lie group integrator.
Further, based on the Castell–Gaines approach we then present uniformly accurate
exponential Lie series integrators that are globally more accurate than their stochastic
Taylor counterpart schemes (these are investigated in detail in Lord, Malham and
Wiese [31] for linear stochastic differential equations). They require the assumption
that a sufficiently accurate underlying ordinary differential integrator is used; that
integrator could for example be an ordinary Lie group Munthe-Kaas method. In
the case of two driving Wiener processes we derive the order 1/2, and in the case
of one driving Wiener process the order 1 uniformly accurate exponential Lie series
integrators. As a consequence we confirm the asymptotic efficiency properties for
both these schemes proved by Castell and Gaines [8] (see Newton [41] for more details
on the concept of asymptotic efficiency). We also present in the case of one driving
Wiener process a new order 3/2 uniformly accurate exponential Lie series integrator
(also see Lord, Malham and Wiese [31]).
We present two physical applications that demonstrate the advantage of using
stochastic Munthe-Kaas methods. First we consider a free rigid body which for ex-
ample could model the dynamics of a satellite. We suppose that it is perturbed by
two independent multiplicative stochastic noise processes. The governing vector fields
are non-commutative and the corresponding exact stochastic flow evolves on the unit
sphere. We show that the stochastic Munthe-Kaas method, with an order 1 stochastic
Taylor integrator used to progress along the corresponding Lie algebra, preserves the
4 Malham and Wiese
approximate solution in the unit sphere manifold to within machine error. However
when an order 1 stochastic Taylor integrator is used directly, the solution leaves the
unit sphere. The contrast between these two methods is more emphatically demon-
strated in our second application. Here we consider an autonomous underwater vehicle
that is also perturbed by two independent multiplicative stochastic noise processes.
The exact stochastic flow evolves on the manifold which is the dual of the Euclidean
Lie algebra se(3); two independent Casimirs are conserved by the exact flow. Again
the stochastic Munthe-Kass method preserves the Casimirs to within machine error.
However the order 1 stochastic Taylor integrator is not only unstable for large step-
sizes, but the approximation drifts off the manifold and makes a dramatic excursion
off to infinity in the embedding space R6.
Preserving the approximate flow on the manifold of the exact dynamics may be a
required property for physical or financial systems driven by smooth or rough paths—
for general references see Iserles, Munthe-Kaas, Nørsett and Zanna [25], Hairer, Lubich
and Wanner [22], Elworthy [13], Lyons and Qian [33] and Milstein and Tretyakov [38].
Stochastic Lie group integrators in the form of Magnus integrators for linear stochastic
differential equations were investigated by Burrage and Burrage [5]. They were also
used in the guise of Möbius schemes (see Schiff and Shnider [43]) to solve stochastic
Riccati equations by Lord, Malham and Wiese [31] where they outperformed direct
stochastic Taylor methods. Further applications where they might be applied include:
backward stochastic Riccati equations arising in optimal stochastic linear-quadratic
control (Kohlmann and Tang [28]); jump diffusion processes on matrix Lie groups
for Bayesian inference (Srivastava, Miller and Grenander [44]); fractional Brownian
motions on Lie groups (Baudoin and Coutin [3]) and stochastic dynamics triggered
by DNA damage (Chickarmane, Ray, Sauro and Nadim [10]).
Our paper is outlined as follows. In Section 2 we present the basic geometric setup,
sans stochasticity. In particular we present a generalized right translation vector field
on a Lie group that forms the basis of our subsequent transformation from the Lie
group to the manifold. Using a Lie group action, this vector field pushes forward
to an infinitesimal Lie group action vector field that generates a flow on the smooth
manifold. In Section 3 we specialize to the case of a matrix Lie group and using
the exponential map, derive the pullback of the generalized right translation vector
field on the Lie group to the corresponding vector field on the Lie algebra. To help
give some context to our overall scheme, we provide in Section 4 illustrative examples
of manifolds and natural choices for associated Lie groups and actions that generate
flows on those manifolds. Then in Section 5 we show how a flow on a smooth manifold
corresponding to a stochastic differential equation can be generated by a stochastic
flow on a Lie algebra via a Lie algebra action. We explicitly present stochastic Munthe-
Kaas Lie group integration methods in Section 6. We start the second half of our
paper by reviewing the exponential Lie series for stochastic differential equations in
Section 7. We show in Section 8 how to construct geometric stochastic Castell–Gaines
numerical methods. In particular we also present uniformly accurate exponential Lie
series numerical schemes that not only can be used as geometric stochastic integrators,
but also are always more accurate than stochastic Taylor numerical schemes of the
corresponding order. In Section 9 we present our concrete numerical examples. Finally
in Section 10 we conclude and present some further future applications and directions.
2. Lie group actions. SupposeM is a smooth finite n-dimensional submanifold
of RN with n ≤ N . We use X(M) to denote the Lie algebra of vector fields on the
manifold M, equipped with the Lie–Jacobi bracket [U, V ] ≡ U · ∇V − V · ∇U , for all
Stochastic Lie group integrators 5
U, V ∈ X(M). Let G denote a finite dimensional Lie group.
Definition 2.1 (Lie group action). A left Lie group action of a Lie group G on a
manifold M is a smooth map Λ: G ×M → M satisfying for all y ∈ M and R,S ∈ G:
(1) Λ(id, y) = y; (2) Λ(R,Λ(S, y)) = Λ(RS, y). We denote Λy ◦ S ≡ Λ(S, y).
Hereafter we suppose y0 ∈ M is fixed and focus on the action map Λy0 : G→M.
We assume that the Lie group action Λ is transitive, i.e. transport across the manifold
from any point y0 ∈ M to any other point y ∈ M can always be achieved via a group
element S ∈ G with y = Λy0 ◦ S (Marsden and Ratiu [37], p. 310).
We define the Lie algebra g associated with the Lie group G to be the vector space
of all right invariant vector fields on G. By standard construction this is isomorphic
to the tangent space to G at the identity id ≡ idG (see Olver [42], p. 48 or Marsden
and Ratiu [37], p. 269).
Definition 2.2 (Generalized right translation vector field). Suppose we are
given a smooth map ξ : M→g. With each such map ξ we associate a vector field
Xξ : G → X(G) defined as follows
Xξ ◦ S ≡ ∂τ exp
τ ξ(Λy0 ◦ S)
for S ∈ G, where ‘exp’ is the usual local diffeomorphism exp: g → G from a neigh-
bourhood of the zero element o ∈ g to a neighbourhood of id ∈ G.
Definition 2.3 (Infinitesimal Lie group action). We associate with each vector
field Xξ : G→X(G) a vector field λξ : M→X(M) as the push forward of Xξ from G to
M by Λy0, i.e. λξ ≡
Xξ, so that if S ∈ G and y = Λy0 ◦ S ∈ M, then
λξ ◦ y ≡ ∂τΛy0 ◦ γ(τ)|τ=0 ,
where γ(t) ∈ G, γ(0) = S and ∂τγ(τ) = Xξ ◦ γ(τ) (the flow generated on G by the
vector field Xξ starting at S ∈ G). Naturally, as a vector field λξ is linear, and also
λξ ◦ y ≡ LXξ ◦ Λy0 ◦ S ,
the Lie derivative of Λy0 along Xξ at S ∈ G.
Remarks.
1. The map Λ(S) : M→M defined by y 7→ Λ(S) ◦ y ≡ Λy ◦ S represents a
flow on M. Hence if y = Λ(S) ◦ y0, the push forward of λξ by Λ(S) is given by
λξ ≡ λAdSξ (Marsden and Ratiu [37], p. 317).
2. We define the isotropy subgroup at y0 ∈ M by Gy0 ≡ {S ∈ G : Λy0◦S = y0}; it
is a closed subgroup of G (see Helgason [23], p. 121 or Warner [48], p. 123). We define
the global isotropy subgroup by GM ≡ ∩y0∈MGy0 ≡ {S ∈ G : Λy0 ◦ S = y0, ∀y0 ∈ M};
it is a normal subgroup of G (see Olver [42], p. 38).
3. A Lie group action is said to be is effective/faithful if the map S 7→ Λ(S) from
G to Diff(M), the group of diffeomorphisms on M, is one-to-one. This is equivalent
to the condition that different group elements have different actions, i.e. GM ≡ {idG}.
A Lie group action is said to be free if Gy0 = {idG} for all y0 ∈ M, i.e. Λy0 is a
diffeomorphism from G to M. For more details see Marsden and Ratiu [37], p. 310
and Olver [42], p. 38.
4. The map γ : G/Gy0→M defined by γ : S · Gy0 7→ Λy0 ◦ S is a diffeomorphism,
i.e. M ∼= G/Gy0 for any y0 ∈ M (a manifold M with a Lie group action Λ: G×M→M
defined over it is thus diffeomorphic to a homogeneous manifold ; see Warner [48],
p. 123 or Olver [42], p. 40). Further, the induced action of G/GM on M is effective.
Hence if Λ is not an effective action of G, we can replace it (without loss of generality)
by the induced action of G/GM (see Olver [42], p. 38).
6 Malham and Wiese
5. Our definition for the generalized right translation vector field Xξ on G is
motivated by the standard right translation vector field used to identify g, the vector
space of right invariant vector fields on G, with TidG, the tangent space to G at the
identity. When ξ ∈ g is constant, Xξ ∈ X(G) is right invariant and a Lie bracket
on TidG can be defined via right extension by the corresponding Lie–Jacobi bracket
for the vector fields Xξ on X(G). Unless ξ ∈ g is constant, Xξ is not in general
right invariant. For further details see Varadarajan [47], Olver [42], or Marsden and
Ratiu [37].
6. The infinitesimal generator map ξ 7→ λξ from g to X(M) is a Lie algebra
homomorphism. If we identify g as the vector space of left invariant vector fields on G
this map becomes an anti-homomorphism. The Lie–Jacobi bracket as defined above
gives the right (rather than left) Lie algebra stucture over the group of diffeomorphisms
on M. If in addition we take the Lie–Jacobi bracket to be minus that defined above—
associated with the left Lie algebra structure—then the infinitesimal generator map
becomes a homomorphism again. See for example Marsden and Ratiu [37], p. 324 or
Munthe-Kaas [40].
7. The image of g under the infinitesimal generator map ξ 7→ λξ forms a finite
dimensional Lie algebra of vector fields on M which is isomorphic to the Lie algebra of
the effectively acting quotient group G/GM (see Olver [42], p. 56). Thus the tangent
space to M at any point is g and M inherents a connection from G/GM. Connections
are necessary to define martingales on manifolds, but not for defining semimartingales
(our focus here); see Malliavin [36] and Emery [14].
8. A comprehensive study of the systematic construction of symmetry Lie groups
from given vector fields can be found in Olver [42].
9. We assumed above that the vector fieldsXξ and λξ are autonomous. However
all results in this and subsequent sections up to Section 7 can be straightforwardly
extended to non-autonomous vector fields generated by ξ : M × R→g with (y, t) 7→
ξ(y, t) for all y ∈ M and t ∈ R.
10. For full generality we want to suspend reference to embedding spaces as far as
possible. However in subsequent sections to be concise we will more explicitly reclaim
this context.
3. Pull back to the Lie algebra. For ease of presentation, we will assume in
this section that G is a matrix Lie group. Recall that the exponential map exp: g → G
is a local diffeomorphism from a neighbourhood of o ∈ g to a neighbourhood of id ∈ G.
Let vξ : g→g be the pull back of the vector field Xξ : G→X(G) from G to g via the
exponential mapping exp: g→G, i.e. vξ ◦ σ ≡ exp∗Xξ ◦ σ. If σ ∈ g then
vξ ◦ σ = dexp−1σ ◦ ξ
Λy0 ◦ expσ
. (3.1)
Here dexp−1σ : g→g is the inverse of the right-trivialized tangent map of the exponential
dexpσ : g→g defined as follows. If β(τ) is a curve in g such that β(0) = σ and
β′(0) = η ∈ g then dexp: g× g→g is the local smooth map (Varadarajan [47], p. 108)
dexpσ ◦ η ≡ ∂τ expβ(τ)|τ=0 exp(−σ)
exp(adσ)− id
◦ η .
Note that as a tangent map dexpσ : g→g is linear. The inverse operator dexp
σ is the
operator series (1.2) generated by considering the reciprocal of dexpσ.
Stochastic Lie group integrators 7
To show that (3.1) is true, if exp: g→G with σ 7→ S = expσ, and β(τ) ∈ g with
β(0) = σ and ∂τβ(τ) = vξ ◦ β(τ), then:
exp∗ vξ ◦ S = ∂τ expβ(τ)|τ=0
dexpσ ◦ vξ ◦ σ
exp(σ)
≡ Xξ ◦ S .
Since ‘exp’ is a diffeomorphism in a neighbourhood of o ∈ g, this push forward calcu-
lation establishes the pull back (3.1) for all σ ∈ g in that neighbourhood.
4. Illustrative examples. Suppose the vector field V : M× R→X(M) gener-
ates a flow solution yt ∈ M starting from y0 ∈ M. Then assume there exists a:
1. Lie group G with corresponding Lie algebra g;
2. Lie group action Λy0 : G→M for which a starting point y0 ∈ M is fixed;
3. Vector field λξ : M× R→X(M) such that: V ≡ λξ, i.e. V is a fundamental
vector field corresponding to the action Λy0 .
Let us suppose G is a matrix Lie group (or can be embedded into a matrix Lie
group, for example the Euclidean group SE(3) is naturally embedded into the special
linear group SL(4;R)). We have for all S ∈ G and t ∈ R,
Xξ(S, t) ≡ ξ
Λy0(S), t
S . (4.1)
If V = λξ for some ξ : M→g, some Lie group G and corresponding action Λy0 , then
the flow generated by Xξ on G drives the flow generated by V on M. In each of the
examples below, given the manifold M, we present a natural Lie group and action
associated with the manifold structure, and identify vector fields which generate flows
on the manifold via the Lie group.
Stiefel manifold Vn,k. Suppose M = Vn,k ≡ {y ∈ Rn×k : yTy = I}. Take
G = SO(n), the special orthogonal group, and Λy0(S) ≡ Sy0, the action of left
multiplication. The corresponding Lie algebra g = so(n). Then by direct calculation
λξ(y) = ξ(y, t) y. Hence if the given vector field V (y, t) = ξ(y, t) y, then the push
forward of the flow generated by Xξ(S, t) on G in (4.1) is the flow generated by V on
M. Note that the unit sphere S2 ∼= V3,1, i.e. S2 is just a particular Stiefel manifold.
In Section 9 as an application, we consider rigid body dynamics evolving on S2.
Isospectral manifold Sn. Suppose M = Sn = {y ∈ Rn×n : yT = y}, the set of
n× n real symmetric matrices. Take G = O(n), the orthogonal group and Λy0(S) ≡
T, which is an isospectral action (Munthe-Kaas [40]). The corresponding Lie
algebra is g = so(n). Again, by direct calculation λξ(y) = ξ(y, t) y − y ξ(y, t). Hence
if the given vector field V (y, t) = ξ(y, t) y−y ξ(y, t), then the push forward of the flow
generated by Xξ(S, t) on G in (4.1) is the flow generated by V on M.
Dual of the Euclidean algebra se(3)∗. Suppose M = se(3)∗ ∼= R3, the dual
of the Euclidean algebra se(3) of the Euclidean group SE(3) =
(s, ρ) ∈ SE(3) : s ∈
SO(3), ρ ∈ R3
. Take G = SE(3) so g = se(3) and Λ ≡ Ad∗ : G × g∗→g∗, the
coadjoint action of G on g∗. Then by direct calculation λξ(y) = −ad∗ξ(y). Since
λξ(y) in linear in ξ and −λξ(y) ≡ λ−ξ(y), it follows that if V (y) = ad∗ξ(y), then
the push forward of the flow generated by X−ξ(S, t) = −ξ
Λy0(S), t
S on G is the
flow generated by V on M. For more details see Section 9 where we investigate the
dynamics of an autonomous underwater vehicle evolving on se(3)∗.
8 Malham and Wiese
Grassmannian manifold Gr(k, n). The Grassmannian manifold M = Gr(k, n)
is the space of k-dimensional subspaces of Rn. Take G = GL(n), the general linear
matrix group, where if S ∈ GL(n), we identify
where the block matrices α, β, γ and δ are sizes k × k, k × (n − k), (n− k) × k and
(n − k) × (n − k), respectively (see Schiff and Shnider [43]; Munthe-Kaas [40]). We
choose the action of GL(n) on Gr(k, n) to be the generalized Möbius transformation
Λy0(S) = (αy0 + β)(γy0 + δ)
−1. Hence if
ξ(t) =
a(t) b(t)
c(t) d(t)
then direct calculation reveals that λξ(y) = a(t)y+ b(t)− yc(t)y− yd(t). Hence if the
given vector field V (y) = a(t)y + b(t) − yc(t)y − yd(t), then the push forward of the
flow generated by Xξ(S, t) = ξ(t)S on G is the flow generated by V on Gr(k, n).
5. Stochastic Lie group integration. We show that if a Lie group action
Λ: G ×M→M exists, then for y0 ∈ M fixed, the Lie algebra action Λy0 ◦ exp: g→M
carries a flow on g to a flow on M.
Theorem 5.1. Suppose there exists a Lie group action Λ: G ×M→M. Then if
there exists a process σ ∈ g and a stopping time T∗ such that on [0, T∗), σ satisfies
the Stratonovich stochastic differential equation
vξi ◦ στ dW iτ , (5.1)
then the process y = Λy0 ◦ expσ ∈ M satisfies the Stratonovich stochastic differential
equation on [0, T∗):
yt = y0 +
λξi ◦ yτ dW iτ . (5.2)
Proof. Using Itô’s lemma, if σt ∈ g satisfies (5.1) then Λy0 ◦ expσt satisfies
Λy0 ◦ expσt = Λy0 ◦ exp o+
Lvξi ◦ Λy0 ◦ expστ dW
Now recall that for each i = 0, 1, . . . , d, Xξi is the push forward of vξi from g to G via
the exponential map, and that λξi is the push forward of Xξi from G to M via Λy0
and so the Lie derivative
Lvξi ◦ Λy0 ◦ expσt ≡ λξi ◦ yt .
Then since yt = Λy0 ◦ expσt, we conclude that y ∈ M is a process satisfying the
stochastic differential equation (5.2).
Corollary 5.2. Suppose that for each i = 0, 1, . . . , d there exists ξi : M→g such
that the vector field Vi : M→X(M) and λξi : M→X(M) can be identified, i.e.
Vi ≡ λξi . (5.3)
Stochastic Lie group integrators 9
Then the push forward by ‘Λy0◦exp’ of the flow on the Lie algebra manifold g generated
by the stochastic differential equation (5.1) is the flow on the smooth manifold M
generated by the stochastic differential equation (5.2), whose solution can be expressed
in the form yt = Λy0 ◦ expσt.
Remark. If the action is free then ‘Λy0 ◦ exp’ is a diffeomorphism from a neigh-
bourhood of o ∈ g to a neighbourhood of y0 ∈ M.
6. Stochastic Munthe-Kaas methods. Assuming that the vector fields in our
original stochastic differential equation (1.1) are fundamental and satisfy (5.3), then
stochastic Munthe-Kaas methods are constructed as follows:
1. Subdivide the global interval of integration [0, T ] into subintervals [tn, tn+1].
2. Starting with t0 = 0, repeat the next two steps over successive intervals
[tn, tn+1] until tn+1 = T .
3. Compute an approximate solution σ̂tn,tn+1 to (5.1) across [tn, tn+1] using a
stochastic Taylor, stochastic Runge–Kutta or Castell–Gaines method.
4. Compute the approximate solution ytn+1 ≈ Λytn ◦ exp σ̂tn,tn+1 .
Note that by construction σ̂tn,tn+1 ∈ g because the stochastic differential equa-
tion (5.1) (or any stochastic Taylor or other sensible approximation) evolves the so-
lution locally on the Lie algebra g via the vector fields vξi : g→g. Suitable methods
for approximating the exponential map to ensure it maps g to G appropriately can be
found in Iserles and Zanna [26]. Then by construction ytn+1 ∈ M.
For example, with two Wiener processes and autonomous vector fields vξi ◦ σ, an
order 1 stochastic Taylor Munthe-Kaas method is based on
σ̂tn,tn+1 =
J0vξ0 +J1vξ1 +J2vξ2 +
+J12vξ1vξ2 +J21vξ2vξ1 +
◦o , (6.1)
evaluated at the zero element o ∈ g. Typically ‘dexp−1σ ’ is truncated to only include
the necessary low order terms to maintain the order of the numerical scheme.
Remark. It is natural to invoke Ado’s Theorem (see for example Olver [42] p. 54):
any finite dimensional Lie algebra is isomorphic to a Lie subalgebra of gl(n) (the
general linear algebra) for some n ∈ N. However as Munthe-Kaas [40] points out,
directly using a matrix representation for the given Lie group might not lead to the
optimal computational implementation (other data structures might do so).
7. Exponential Lie series. The stochastic Taylor series is known in different
contexts as the Neumann series, Peano–Baker series or Feynman–Dyson path ordered
exponential. If the vector fields in the stochastic differential equation (1.1) are au-
tonomous (which we assume henceforth), i.e. for all i = 0, 1, . . . , d, Vi = Vi(y) only,
then the stochastic Taylor series for the flow is
Jα1···αm(t)Vα1 · · ·Vαm .
Here Pm is the set of all combinations of multi-indices α = (α1, . . . , αm) of length m
with αi ∈ {0, 1, . . . , d} and
Jα1···αm(t) ≡
· · ·
∫ τm−1
dWα1τm · · · dW
are multiple Stratonovich integrals.
10 Malham and Wiese
The logarithm of ϕt is the exponential Lie series, Magnus expansion (Magnus [34])
or Chen–Strichartz formula (Chen [9], Strichartz [45]). In other words we can express
the flow map in the form ϕt = expψt, where
Ji(t)Vi +
j>i=0
(Jij − Jji)(t)[Vi, Vj ] + · · ·
is the exponential Lie series for our system, and [· , ·] is the Lie–Jacobi bracket on
X(M). See Yamato [49], Kunita [29], Ben Arous [1] and Castell [7] for the derivation
and convergence of the exponential Lie series expansion in the stochastic context;
Strichartz [45] for the full explicit expansion; Sussmann [46] for a related product
expansion and Lyons [32] for extensions to rough paths.
Let us denote the truncated exponential Lie series by
ψ̂t =
Jα cα , (7.1)
where Qm denotes the finite set of multi-indices α for which ‖Jα‖L2 is of order up to
and including tm, where m = 1/2, 1, 3/2, . . .. The terms cα are linear combinations
of finitely many (length α) products of the smooth vector fields Vi, i = 0, 1, . . . , d.
The following asymptotic convergence result can be established along the lines of the
proof for linear stochastic differential equations in Lord, Malham and Wiese [31]; we
provide a proof in Appendix A.
Theorem 7.1. Assume the vector fields Vi have 2m+1 uniformly bounded deriva-
tives, for all i = 0, 1, . . . , d. Then for t ≤ 1, the flow exp ψ̂t ◦ y0 is square-integrable,
where ψ̂t is the truncated Lie series (7.1). Further, if y is the solution of the stochastic
differential equation (1.1), there exists a constant C
m, ‖y0‖2
such that
∥yt − exp ψ̂t ◦ y0
m, ‖y0‖2
tm+1/2 . (7.2)
8. Geometric Castell–Gaines methods. Consider the truncated exponential
Lie series ψ̂tn,tn+1 across the interval [tn, tn+1]. We approximate higher order multiple
Stratonovich integrals across each time-step by their expectations conditioned on the
increments of the Wiener processes on suitable subdivisions (Gaines and Lyons [20]).
An approximation to the solution of the stochastic differential equation (1.1) across
the interval [tn, tn+1] is given by the flow generated by the truncated and conditioned
exponential Lie series ψ̂tn,tn+1 via
ytn+1 ≈ exp
ψ̂tn,tn+1
◦ ytn .
Hence the solution to the stochastic differential equation (1.1) can be approximately
computed by solving the ordinary differential system (see Castell and Gaines [8];
Misawa [39])
u′(τ) = ψ̂tn,tn+1 ◦ u(τ) (8.1)
across the interval τ ∈ [0, 1]. Then if u(0) = ytn we will get u(1) ≈ ytn+1. We
must choose a sufficiently accurate ordinary differential integrator to solve (8.1)—we
implicitly assume this henceforth.
Stochastic Lie group integrators 11
The set of governing vector fields Vi, i = 0, 1, . . . , d, prescribes a map from the
driving path process w ≡ (W 1, . . . ,W d) to the unique solution process y ∈ M to the
stochastic differential equation (1.1). The map w 7→ y is called the Itô map. Recall
that we assume the vector fields are smooth. When there is only one driving Wiener
process (d = 1) the Itô map is continuous in the topology of uniform convergence
(Theorem 1.1.1. in Lyons and Qian [33]). When there are two or more driving pro-
cesses (d ≥ 2) the Universal Limit Theorem (Theorem 6.2.2. in Lyons and Qian [33])
tells us that the Itô map is continuous in the p-variation topology, in particular for
2 ≤ p < 3. A Wiener path with d ≥ 2 has p-variation with p > 2, and the p-variation
metric in this case includes information about the Lévy chordal areas of the path
(Lyons [32]). Hence we must choose suitable piecewise smooth approximations to the
driving path process w. The following result follows from the corresponding result
for ordinary differential equations in Hairer, Lubich and Wanner [22] (p. 112) as well
as directly from Chapter VIII in Malliavin [36] on the Transfer Principle (see also
Emery [15]).
Lemma 8.1. A necessary and sufficient condition for the solution to the stochastic
differential equation (1.1) to evolve on a smooth n-dimensional submanifold M of
RN (n ≤ N) up to a stopping time T∗ is that Vi(y, t) ∈ TyM for all y ∈ M, i =
0, 1, . . . , d.
Hence the stochastic Taylor expansion for the flow ϕt is a diffeomorphism on M.
However a truncated version of the stochastic Taylor expansion for the flow ϕ̂t will not
in general keep you on the manifold, i.e. if y0 ∈ M then ϕ̂t ◦ y0 need not necessarily
lie in M. On the other hand, the exponential Lie series ψt, or any truncation ψ̂t of
it, lies in X(M). By Lemma 8.1 this is a necessary and sufficient condition for the
corresponding flow-map exp ψ̂t to be a diffeomorphism on M. Hence if u(0) = ytn ∈
M, then ytn+1 ≈ u(1) ∈ M. When solving the ordinary differential equation (8.1),
classical geometric integration methods, for example Lie group integrators such as
Runge–Kutta Munthe-Kaas methods, over the interval τ ∈ [0, 1] will numerically
ensure ytn+1 stays in M. Additionally, as the following result reveals, numerical
methods constructed using the Castell–Gaines Lie series approach can also be more
accurate (a proof is provided in Appendix B). We define the strong global error at
time T associated with an approximate solution ŷT as E ≡ ‖yT − ŷT ‖L2 .
Theorem 8.2. In the case of two independent Wiener processes and under the
assumptions of Theorem 7.1, for any initial condition y0 ∈ M and a sufficiently
small fixed stepsize h = tn+1 − tn, the order 1/2 Lie series integrator is globally more
accurate in L2 than the order 1/2 stochastic Taylor integrator. In addition, in the
case of one Wiener process, the order 1 and 3/2 uniformly accurate exponential Lie
series integrators generated by ψ̂
tn,tn+1
= J0V0 + J1V1 +
[V1, [V1, V0]]
(3/2)
tn,tn+1
= J0V0 + J1V1 +
(J01 − J10)[V0, V1] + h
[V1, [V1, V0]]
respectively, are globally more accurate in L2 than their corresponding stochastic Tay-
lor integrators. In other words, if E lsm denotes the global error of the exponential Lie
series integrators of order m above, and Estm is the global error of the stochastic Taylor
integrators of the corresponding order, then E lsm ≤ Estm for m = 1/2, 1, 3/2.
Remarks.
1. The result for ψ̂(3/2) is new. That the order-1/2 Lie series integrator (for two
Wiener processes) and the order 1 integrator generated by ψ̂(1) are uniformly more
accurate confirms the asymptotically efficient properties of these schemes proved by
12 Malham and Wiese
Castell and Gaines [8]. The proof follows along the lines of an analogous result for
linear stochastic systems considered in Lord, Malham and Wiese [31].
2. Consider the order 1/2 exponential Lie series with no vector field commu-
tations. Solving the ordinary differential equation (8.1) using an (ordinary) Euler
Munthe-Kaas method and approximating dexp
σ ≈ id is equivalent to the order 1/2
stochastic Taylor Munthe-Kaas method (for the same Lie group and action).
9. Numerical examples.
9.1. Rigid body. We consider the dynamics of a rigid body such as a satellite
(see Marsden and Ratiu [37]). We will suppose that the rigid body is perturbed by two
independent multiplicative stochastic processes W 1 and W 2 with the corresponding
vector fields Vi(y) ≡ ξi(y) y, for i = 0, 1, 2, with ξi ∈ so(3). If we normalize the initial
data y0 so that |y0| = 1 then the dynamics evolves on M = S2. We naturally suppose
G = SO(3), and Λy0(S) ≡ Sy0 so that λξi(y) = ξi(y) y, and we can pull back the
flow generated by V on M to the flow on G generated by Xξi(S, t) = ξi
Λy0(S)
i = 0, 1, 2. We use the following matrix representation for the ξi(y) ∈ so(3):
ξi(y) =
0 −y3/αi,3 y2/αi,2
y3/αi,3 0 −y1/αi,1
−y2/αi,2 y1/αi,1 0
where the constants αi,j for j = 1, 2, 3 are chosen so that the vector fields Vi and
matrices ξi do not commute for i = 0, 1, 2: α0,1 = 3, α0,2 = 1, α0,3 = 2, α1,1 = 1,
α1,2 = 1/2, α1,3 = 3/2, α2,1 = 1/4, α2,2 = 1, α2,3 = 1/2. The vector fields Vi satisfy
the conditions of Theorem 7.1 since the manifold is compact in this case.
We will numerically solve (1.1) using three different order 1 methods: stochastic
Taylor, stochastic Taylor Munthe-Kaas based on (6.1) and Castell–Gaines (a stan-
dard non-geometric Runge–Kutta method is used to solve the ordinary differential
equation (8.1)). The vector field compositions ViVj needed for the stochastic Taylor
and Castell–Gaines methods are readily computed. For the Munthe-Kaas method we
note that we have vξi ◦ o = ξi(y0) and
vξivξj ◦ o = Â(y0, y0;αi, αj)− 12 [ξi(y0), ξj(y0)] .
Here o ∈ so(3) is the zero element on the Lie algebra, and for all y, z ∈ R3 we define
A(y, z;α, β) ≡
− y3z2
− y1z3
− y2z1
and ˆ : R3→so(3) denotes the vector space isomorphism σ 7→ σ̂ where
0 −σ3 σ2
σ3 0 −σ1
−σ2 σ1 0
Note that ŷ z ≡ y ∧ z (see Marsden and Ratiu [37]). Note also since σ ∈ so(3),
expσ ∈ SO(3) can be conveniently and cheaply computed using Rodrigues’ formula
(see Marsden and Ratiu [37] or Iserles et al. [25]).
In Figure 9.1 we show the distance from the manifold S2 of each the three approx-
imations; we start with initial data y0 = (
2, 0)T. The stochastic Taylor Munthe-
Kaas method can be seen to preserve the solution in the unit sphere to within machine
Stochastic Lie group integrators 13
0 1 2 3 4 5 6 7 8 9 10
Stochastic Taylor
Castell−Gaines
Munthe−Kaas
Fig. 9.1. Rigid body: We show the log-distance of the approximate solution to the unit sphere
as a function of time for each of the methods. Below we show the approximate solutions as a function
of time for the stochastic Taylor (blue) and Munthe-Kaas methods (magenta). The trajectory starts
at the top right and eventually drifts over the left horizon.
14 Malham and Wiese
error. We also see that the stochastic Taylor method clearly drifts off the sphere as
the integration time progresses, as does the non-geometric Castell-Gaines method—
which does however remain markedly closer to the manifold than the stochastic Taylor
scheme.
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Stochastic Taylor
Castell−Gaines
Munthe−Kaas
−3.2 −3.1 −3 −2.9 −2.8 −2.7 −2.6 −2.5 −2.4 −2.3 −2.2
(stepsize)
Number of sampled paths=100
Stochastic Taylor
Castell−Gaines
Munthe−Kaas
Fig. 9.2. Autonomous underwater vehicle: We show the log-distance of the approximate solution
to the two Casimirs C1 = π ·p (dotted line) and C2 = |p|
2 (solid line) as a function of time for each
of the methods. Below, we also show the global error as a function of stepsize.
Stochastic Lie group integrators 15
9.2. Autonomous underwater vehicle. The dynamics of an ellipsoidal au-
tonomous underwater vehicle is prescribed by the state y = (π, p) ∈ se(3)∗ where
π ∈ so(3)∗ is its angular momentum and p ∈ (R3)∗ its linear momentum (see Holmes,
Jenkins and Leonard [24], Egeland, Dalsmo and Sørdalen [12] and Marsden and
Ratiu [37]). We suppose that the vehicle is perturbed by two independent multi-
plicative stochastic processes. The governing vector fields are for i = 0, 1, 2:
Vi(y) = ad
◦ y .
Here ξi(y) =
ωi(y), ui(y)
∈ se(3) where ωi(y) = I−1i π and ui(y) = M
i p are the
angular and linear velocity, and Ii = diag(αi,1, αi,2, αi,3) andMi = diag(βi,1, βi,2, βi,3)
are the constant moment of inertia and mass matrices, respectively. Explicitly for
ξ ∈ se(3) we have
ad∗ξ ◦ y ≡ (π ∧ ω + p ∧ u, p ∧ ω) .
The system of vector fields Vi, i = 0, 1, 2 represents the Lie–Poisson dynamics on
M = se(3)∗ (Marsden and Ratiu [37]). There are two independent Casimir functions
Ck : se(3)
∗→R, k = 1, 2, namely C1 = π · p and C2 = |p|2; these are conserved by the
flow on se(3)∗. Note that the Hamiltonian, i.e. total kinetic energy 1
(π · ω + p · u),
is also exactly conserved (and helpful for establishing the sufficiency conditions in
Theorem 7.1), but that is not our focus here.
If G = SE(3) ∼= SO(3) × R3, then the coadjoint action of SE(3) on se(3)∗,
: SE(3) × se(3)∗→se(3)∗ is defined for all S = (s, ρ) ∈ SE(3), where s ∈ SO(3)
and ρ ∈ R3, and y ∈ se(3)∗ by: Λy ◦ S = Ad∗S−1 ◦ y ≡
sπ + ρ ∧ (sp), sp
. The
corresponding infinitesimal action λ : se(3)× se(3)∗→se(3)∗ for all ξ ∈ se(3) and y ∈
se(3)∗ is given by (see Marsden and Ratiu [37], p. 477)
λξ ◦ y = −ad∗ξ ◦ y .
Since ad∗ξ(y) = −λξ(y) = λ−ξ(y) the governing set of vector fields on se(3)∗ are
Vi(y) = λ−ξi ◦ y .
We can now pull back this flow on se(3)∗ to a flow on SE(3) via Λy0 . The correspond-
ing flow on SE(3) is generated by the governing set of vector fields for i = 0, 1, 2:
X−ξi ◦ S = −
ωi(y) ∧ s, ωi(y) ∧ ρ+ ui(y)
with y = Λy0(S).
To aid implementation note that SE(3) =
(s, ρ) ∈ SE(3) : s ∈ SO(3), ρ ∈ R3
embeds into SL(4;R) via the map
S = (s, ρ) 7→
where O is the three-vector of zeros. Also se(3) is isomorphic to a Lie subalgebra of
sl(4;R) with elements of the form
σ = (θ, ζ) 7→
16 Malham and Wiese
Hence the governing vector fields on SE(3) are of the form Xξi = −ξi(y)S, where
ξi(y) =
ω̂i(π) ui(p)
The governing vector fields on se(3) are vξi(σ) = −dexpσ ◦ ξi
Λy0(expσ)
. Again the
vector field compositions ViVj needed for the stochastic Taylor and Castell–Gaines
methods can be computed straightforwardly. Direct calculation also reveals that in
block matrix form
vξivξj◦o =
Â(π0, π0;αi, αj) + Â(p0, p0;βi, αj) A(π0, p0;αi, βj)
[ξi(y0), ξj(y0)] .
Here A(y, z;α, β) is defined as for the rigid body example. Note that the exponential
map exp
se(3) : se(3)→SE(3) is defined for all σ = (θ, ζ) ∈ se(3) by
se(3) σ =
so(3) θ̂ f(θ)ζ
where exp
so(3) is the exponential map from so(3) to SO(3) which can be computed
using Rodrigues’ formula and (see Bullo and Murray [4], p. 5)
f(θ) = I3×3 + (1 − cos ‖θ‖)θ̂/‖θ‖2 +
1− (sin ‖θ‖)/‖θ‖
θ̂2/‖θ‖2 .
In Figure 9.2 we show the distance from the manifold se(3)∗ of each the three ap-
proximations; in particular how far the individual trajectories stray from the Casimirs
C1 = π · p and C2 = |p|2. We start with the initial data y0 = (
2, 0, 0,
As before the stochastic Taylor Munthe-Kaas method can be seen to preserve the
Casimirs to within machine error. We also see that the stochastic Taylor method
clearly drifts off the manifold as the integration time progresses and at a particular
time depending on the Wiener path shoots off very rapidly away from the manifold.
Note also that for large stepsizes the stochastic Taylor method is unstable. However
the non-geometric Castell–Gaines and stochastic Munthe-Kaas methods still give reli-
able results in that regime. Lastly, although the the stochastic Munthe-Kaas method
adheres to the manifold to within machine error, the error of the non-geometric
Castell–Gaines method is actually smaller.
10. Conclusions. We have established and implemented stochastic Lie group
integrators based on stochastic Munthe-Kaas methods and also derived geometric
Castell–Gaines methods. We have also revealed several aspects of these integrators
that require further investigation.
1. We could construct a stochastic nonlinear Magnus method by approximating
the solution to the stochastic differential equation (5.1) on the Lie algebra using Picard
iterations (see Casas and Iserles [6]).
2. We would like to develop a practical procedure for implementing ordinary
Munthe-Kaas methods for higher order Castell–Gaines integrators. We need to de-
termine the element ξ : M→g so that in (8.1) we have ψ̂ = λξ.
3. We need to determine the properties of the local and global errors for the
stochastic Munthe-Kaas methods. Also a thorough investigation of the stability prop-
erties of the stochastic Munthe-Kaas and Castell–Gaines methods is required. For
the autonomous underwater vehicle simulations they were both superior to the direct
stochastic Taylor method, especially for larger stepsizes. We also need to compare the
relative efficiency of the methods concerned, in particular to compare an optimally
efficient geometric Castell–Gaines method with the stochastic Munthe-Kaas method.
Stochastic Lie group integrators 17
4. Although we have chiefly confined ourselves to driving paths that are Wiener
processes, we can extend Munthe-Kaas and Castell–Gaines methods to rougher driv-
ing paths (Lyons and Qian [33], Friz [18], Friz and Victoir [19]). Further, what hap-
pens when we consider processes involving jumps? For example Srivastava, Miller and
Grenander [44] consider jump diffusion processes on matrix Lie groups for Bayesian
inference. Or what if we consider fractional Brownian driving paths; Baudoin and
Coutin [3] investigate fractional Brownian motions on Lie groups?
5. Schiff and Shnider [43] have used Lie group methods to derive Möbius schemes
for numerically integrating deterministic Riccati systems beyond finite time removable
singularities and numerical instabilities. They integrate a linear system of equations
on the general linear group GL(n) which corresponds to a Riccati flow on the Grass-
mannian manifold Gr(k, n) via the Möbius action map. Lord, Malham and Wiese [31]
implemented stochastic Möbius schemes and show that they can be more accurate and
cost effective than directly solving stochastic Riccati systems using stochastic Taylor
methods. We would like to investigate further their effectiveness for stochastic Ric-
cati equations arising in Kalman filtering (Kloeden and Platen [27]) and to backward
stochastic Riccati equations arising in optimal stochastic linear-quadratic control (see
for example Kohlmann and Tang [28] and Estrade and Pontier [16]).
6. Other areas of potential application of the methods we have presented in this
paper are for example: term-structure interest rate models evolving on finite dimen-
sional invariant manifolds (see Filipovic and Teichmann [17]); stochastic dynamics
triggered by DNA damage (Chickarmane, Ray, Sauro and Nadim [10]) and stochastic
symplectic integrators for which the gradient of the solution evolves on the symplectic
Lie group (see Milstein and Tretyakov [38]).
Acknowledgments. We thank Alex Dragt, Peter Friz, Anders Hansen, Terry
Lyons, Per-Christian Moan and Hans Munthe–Kaas for stimulating discussions. We
also thank the anonymous referees, whose suggestions and encouragement improved
the original manuscript significantly. SJAM would like to acknowledge the invalu-
able facilities of the Isaac Newton Institute where some of the final touches to this
manuscript were completed.
Appendix A. Proof of Theorem 7.1. We follow the proof for linear stochastic
differential equations in Lord, Malham and Wiese [31] (where further technical details
on estimates for multiple Stratonovich integrals can be found). Suppose ψ̂t ≡ ψ̂t(m)
is the truncated Lie series (7.1). First we show that exp ψ̂t ◦ y0 ∈ L2. We see that
for any number k,
)k ◦ y0 is a sum of |Qm|k terms, each of which is a k-multiple
product of terms Jα cα ◦ y0. It follows that
)k ◦ y0
‖cα ◦ y0‖
αi∈Qm
i=1,...,k
‖Jα1Jα2 · · · Jαk‖L2 . (A.1)
Note that the maximum of the norm of the compositions of vector fields cα◦y0 is taken
over a finite set. Repeated application of the product rule reveals that for i = 1, . . . , k,
each term ‘Jα1Jα2 · · · Jαk ’ in (A.1) is the sum of at most 22mk−1 Stratonovich integrals
Jβ , where for t ≤ 1, ‖Jβ‖L2 ≤ 24mk−1 tk/2. Since the right hand side of equation (A.1)
consists of |Qm|k 22mk−1 Stratonovich integrals Jβ , we conclude that,
)k ◦ y0
‖cα ◦ y0‖ · |Qm| · 26m · t1/2
18 Malham and Wiese
Hence exp ψ̂t ◦ y0 is square-integrable.
Second we prove (7.2). Let ŷt denote the stochastic Taylor series solution, trun-
cated to included terms of order up to and including tm. We have
∥yt − exp ψ̂t ◦ y0
∥yt − ŷt
∥ŷt − exp ψ̂t ◦ y0
We know yt ∈ L2—see Lemma III.2.1 in Gihman and Skorohod [21]. Note that
the assumptions there are fulfilled, since the uniform boundedness of the derivatives
implies uniform Lipschitz continuity of the vector fields by the mean value theorem,
and uniform Lipschitz continuity in turn implies a linear growth condition for the
vector fields since they are autonomous. Note that ŷt is a strong approximation to
yt up to and including terms of order t
m, with the remainder consisting of O(tm+1/2)
terms (see Proposition 5.9.1 in Kloeden and Platen [27]). It follows from the definition
of the exponential Lie series as the logarithm of the stochastic Taylor series, that the
terms of order up to and including tm in exp ψ̂t ◦ y0 correspond with ŷt; the error
consists of O(tm+1/2) terms.
Appendix B. Proof of Theorem 8.2. Our proof follows along the lines of that
for uniformly accurate Magnus integrators for linear constant coefficient systems (see
Lord, Malham & Wiese [31] and Malham and Wiese [35]). Let ϕtn,tn+1 and ϕ̂tn,tn+1
denote the exact and approximate flow-maps constructed on the interval [tn, tn+1] of
length h. We define the local flow remainder as
Rtn,tn+1 ≡ ϕtn,tn+1 − ϕ̂tn,tn+1 ,
and so the local remainder is Rtn,tn+1 ◦ ytn . Let Rls and Rst denote the local flow
remainders corresponding to the exponential Lie series and stochastic Taylor approx-
imations, respectively.
B.1. Order 1/2 integrator: two Wiener processes. For the global order 1/2
integrators we have to leading order Rls = 1
(J12 − J21)[V1, V2] and Rst = J12V1V2 +
J21V2V1. Note that we have included the terms J11V
1 and J22V
2 in the integrators.
A direct calculation reveals that
(Rst ◦ y0)TRst ◦ y0
(Rls ◦ y0)TRls ◦ y0
+ h2mUTBU +O
. (B.1)
Here m = 1/2 (for the order 1/2 integrators), U = (V1V2 ◦ y0, V2V1 ◦ y0)T ∈ R2n, and
B ∈ R2n×2n consists of n× n diagonal blocks of the form bijIn×n where
b = 1
and In×n is the n×n identity matrix. Since b is positive semi-definite, the matrix B =
b⊗In×n is positive semi-definite. Hence the order 1/2 exponential Lie series integrator
is locally more accurate than the corresponding stochastic Taylor integrator.
B.2. Order 1 integrator: one Wiener process. For the global order 1 in-
tegrators we have to leading order Rls = 1
(J01 − J10)[V0, V1] and Rst = J01V0V1 +
J10V1V0 + J111V
h2(V0V
1 + V
1 V0). The terms of order h
2 shown are significant
when we consider the global error in Section B.4 below. The estimate (B.1) also
applies in this case with m = 1 and U = (V0V1 ◦ y0, V1V0 ◦ y0, V 31 ◦ y0)T ∈ R3n; and
B ∈ R3n×3n consists of n× n diagonal blocks of the form bijIn×n where
b = 1
3 3 3
3 3 3
3 3 5
Stochastic Lie group integrators 19
Since b is positive semi-definite, the matrix B = b ⊗ In×n is positive semi-definite.
Hence the order 1 exponential Lie series integrator is locally more accurate than the
corresponding stochastic Taylor integrator.
B.3. Order 3/2 integrator: one Wiener process. The local flow remainders
are Rls = 1
J110−2J101+J011− 12h
[V1, [V1, V0]] and R
st = J011V0V
1 +J101V1V0V1+
J110V
1 V0 + J1111V
1 − 14h
2(V0V
1 + V
1 V0 +
V 41 ). The terms of order h
2 shown are
significant when we consider the global error—but for a different reason this time—see
Section B.4 below. Again, the estimate (B.1) applies in this case with m = 3/2 and
U = (V0V
1 ◦ y0, V1V0V1 ◦ y0, V 21 V0 ◦ y0, V 41 ◦ y0)T ∈ R4n; and B ∈ R4n×4n consists of
n× n diagonal blocks of the form bijIn×n where
b = 1
11 8 5 12
8 8 8 12
5 8 11 12
12 12 12 24
Again, B is positive semi-definite and the order 3/2 exponential Lie series integrator
is locally more accurate than the corresponding stochastic Taylor integrator.
B.4. Global error. Recall that we define the strong global error at time T
associated with an approximate solution ŷT as E ≡ ‖yT − ŷT ‖L2. The exact and
approximate solutions can be constructed by successively applying the exact and
approximate flow maps ϕtn,tn+1 and ϕ̂tn,tn+1 on the successive intervals [tn, tn+1] to
the initial data y0. A straightforward calculation shows for a small fixed stepsize h,
E2 = E (R ◦ y0)TR ◦ y0 , (B.2)
up to higher order terms, where R ≡
n=0 ϕtn+1,tN ◦Rtn,tn+1 ◦ϕt0,tn is the standard
accumulated local error contribution to the global error. The important conclusion is
that when we construct the global error (B.2), the terms of leading order in the local
flow remainders Rls or Rst with zero expectation lose only a half order of convergence
in this accumulation effect. Hence in the local flow remainders shown above, for the
terms of zero expectation, the local superior accuracy for the Lie series integrators
transfers to the corresponding global errors (see Lord, Malham and Wiese [31] for
more details). Terms of non-zero expectation however behave like deterministic er-
ror terms losing a whole order (in the local to global convergence); they contribute
to the global error through their expectations. Hence we include such terms of or-
der h2 in the order 3/2 integrators above and they appear as the terms subtracted
from the remainders shown. For the order 1 integrators we do not need to include
the order h2 terms in the integrator to obtain the correct mean-square convergence.
However to guarantee that the global error for the exponential Lie series integrator
is always smaller than that for the stochastic Taylor scheme, we include this term in
the integrator.
REFERENCES
[1] G. Ben Arous, Flots et series de Taylor stochastiques, Probab. Theory Related Fields, 81
(1989), pp. 29–77.
[2] F. Baudoin, An introduction to the geometry of stochastic flows, Imperial College Press, 2004.
[3] F.Baudoin and L. Coutin, Self-similarity and fractional Brownian motions on Lie groups,
arXiv:math.PR/0603199 v1, 2006.
http://arxiv.org/abs/math/0603199
20 Malham and Wiese
[4] F. Bullo and R. M. Murray, Proportional derivative (PD) control on the Euclidean group,
CDS Technical Report 95-010, 1995.
[5] K. Burrage and P. M. Burrage, High strong order methods for non-commutative stochas-
tic ordinary differential equation systems and the Magnus formula, Phys. D, 133 (1999),
pp. 34–48.
[6] F. Casas and A. Iserles, Explicit Magnus expansions for nonlinear equations, Cambridge NA
reports, 2005.
[7] F. Castell, Asymptotic expansion of stochastic flows, Probab. Theory Related Fields, 96
(1993), pp. 225–239.
[8] F. Castell and J. Gaines, An efficient approximation method for stochastic differential equa-
tions by means of the exponential Lie series, Math. Comp. Simulation, 38 (1995), pp. 13–19.
[9] K. T. Chen, Integration of paths, geometric invariants and a generalized Baker–Hausdorff
formula, Annals of Mathematics, 65(1) (1957), pp. 163–178.
[10] V. Chickarmane, A. Ray, H. M. Sauro and A. Nadim, A model for p53 dynamics triggered
by DNA damage, SIAM J. Applied Dynamical Systems, 6(1) (2007), pp.61–78.
[11] P. E. Crouch and R. Grossman, Numerical integration of ordinary differential equations on
manifolds, J. Nonlinear Sci., 3 (1993), pp. 1–33.
[12] O. Egeland, M. Dalsmo and O. J. Sørdalen, Feedback control of a nonholonomic under-
water vehicle with a constant desired configuration, The International Journal of Robotics
Research, 15(1) (1996), pp. 24–35.
[13] K. D. Elworthy, Stochastic differential equations on manifolds, London Mathematical Society
Lecture Note Series 70, Cambridge University Press, 1982.
[14] M. Emery, Stochastic Calculus on manifolds, Universitext, Springer–Verlag, 1989.
[15] , On two transfer principles in stochastic differential geometry, Séminaire de probabilités
(Strasbourg), 24 (1990), pp. 407–441.
[16] A. Estrade and M. Pontier, Backward stochastic differential equations in a Lie group,
Séminaire de probabilités (Strasbourg), 35 (2001), pp. 241–259.
[17] D. Filipović and J. Teichmann, On the geometry of the term structure of interest rates, Proc.
R. Soc. Lond. A, 460 (2004), pp. 129–167.
[18] P. Friz, Continuity of the Itô-map for Hölder rough paths with applications to the support
theorem in Hölder norm, arXiv:math.PR/0304501 v2, 2003.
[19] P. Friz and N. Victoir, Euler estimates for rough differential equations, Preprint, 2007.
[20] J. G. Gaines and T. J. Lyons, Variable step size control in the numerical solution of stochastic
differential equations, SIAM J. Appl. Math., 57(5) (1997), pp. 1455–1484.
[21] I. I. Gihman, and A. V. Skorohod, The theory of stochastic processes III, Springer, 1979.
[22] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration, Springer Series
in Computational Mathematics, 2002.
[23] S. Helgason, Differential geometry, Lie groups, and symmetric spaces, Academic Press, 1978.
[24] P. Holmes, J. Jenkins and N. E. Leonard, Dynamics of the Kirchoff Equations I: coincident
centers of gravity and bouyancy, Phys. D, 118 (1998), pp. 311–342.
[25] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna, Lie-group methods, Acta
Numer., (2000), pp. 215–365.
[26] A. Iserles and A. Zanna, Efficient computation of the matrix exponential by generalized polar
decompositions, SIAM J. Numer. Anal., 42(5) (2005), pp. 2218–2256.
[27] P. E. Kloeden and E. Platen, Numerical solution of stochastic differential equations,
Springer, 1999.
[28] M. Kohlmann and S. Tang, Multidimensional backward stochastic Riccati equations and ap-
plications, SIAM J. Control Optim., 41(6) (2003), pp. 1696–1721.
[29] H. Kunita, On the representation of solutions of stochastic differential equations, LNM 784,
Springer–Verlag, 1980, pp. 282–304.
[30] , Stochastic flows and stochastic differential equations, Cambridge University Press, 1990.
[31] G. Lord, S. J.A. Malham and A. Wiese, Efficient strong integrators for linear stochastic
systems, 2006, Submitted.
[32] T. Lyons, Differential equations driven by rough signals, Rev. Mat. Iberoamericana, 14(2)
(1998), pp. 215–310.
[33] T. Lyons and Z. Qian, System control and rough paths, Oxford University Press, 2002.
[34] W. Magnus, On the exponential solution of differential equations for a linear operator, Comm.
Pure Appl. Math., 7 (1954), pp. 649–673.
[35] S. J.A. Malham and A. Wiese, Universal optimal stochastic expansions, 2007, Preprint.
[36] P. Malliavin, Stochastic analysis, Grundlehren der mathematischen Wissenschaften 313,
Springer, 1997.
[37] J. E. Marsden and T. S. Ratiu, Introduction to mechanics and symmetry, Second edition,
http://arxiv.org/abs/math/0304501
Stochastic Lie group integrators 21
Springer, 1999.
[38] G. N. Milstein and M. V. Tretyakov, Stochastic numerics for mathematical physics,
Springer, 2004.
[39] T. Misawa, A Lie algebraic approach to numerical integration of stochastic differential equa-
tions, SIAM J. Sci. Comput., 23(3) (2001), pp. 866–890.
[40] H. Munthe-Kaas, High order Runge–Kutta methods on manifolds, Appl. Numer. Math., 29
(1999), pp. 115–127.
[41] N. J. Newton, Asymptotically efficient Runge–Kutta methods for a class of Itô and
Stratonovich equations, SIAM J. Appl. Math., 51 (1991), pp. 542–567.
[42] P. J. Olver, Equivalence, invariants, and symmetry, Cambridge University Press, 1995.
[43] J. Schiff and S. Shnider, A natural approach to the numerical integration of Riccati differ-
ential equations, SIAM J. Numer. Anal., 36(5) (1999), pp. 1392–1413.
[44] A. Srivastava, M. I. Miller and U. Grenander, Jump-diffusion processes on matrix Lie
groups for Bayesian inference, preprint, 2000.
[45] R. S. Strichartz, The Campbell–Baker–Hausdorff–Dynkin formula and solutions of differen-
tial equations, J. Funct. Anal., 72 (1987), pp. 320–345.
[46] H. J. Sussmann, Product expansions of exponential Lie series and the discretization of stochas-
tic differential equations, in Stochastic Differential Systems, Stochastic Control Theory,
and Applications, W. Fleming and J. Lions, eds., Springer IMA Series, Vol. 10 (1988),
pp. 563–582.
[47] V. S. Varadarajan, Lie groups, Lie algebras, and their representations, Springer, 1984.
[48] F. W. Warner, Foundations of differentiable manifolds and Lie groups, Graduate Texts in
Mathematics, Springer–Verlag, 1983.
[49] Y. Yamato, Stochastic differential equations and nilpotent Lie algebras, Z. Wahrsch. Verw.
Gebiete, 47(2) (1979), pp 213–229.
