Symmetric spaces and Lie triple systems in numerical analysis of differential equations

A remarkable number of different numerical algorithms can be understood and analyzed using the concepts of symmetric spaces and Lie triple systems, which are well known in differential geometry from the study of spaces of constant curvature and their tangents. This theory can be used to unify a range of different topics, such as polar-type matrix decompositions, splitting methods for computation of the matrix exponential, composition of selfadjoint numerical integrators and dynamical systems with symmetries and reversing symmetries. The thread of this paper is the following: involutive automorphisms on groups induce a factorization at a group level, and a splitting at the algebra level. In this paper we will give an introduction to the mathematical theory behind these constructions, and review recent results. Furthermore, we present a new Yoshida-like technique, for self-adjoint numerical schemes, that allows to increase the order of preservation of symmetries by two units. Since all the time-steps are positive, the technique is particularly suited to stiff problems, where a negative time-step can cause instabilities.


Introduction
In numerical analysis there exist numerous examples of objects forming a group, i.e. objects that compose in an associative manner, have an inverse and identity element. Typical examples are the group of orthogonal matrices or the group of Runge-Kutta methods. Semigroups, sets of objects close under composition but not inversion, like for instance the set of all matrices and explicit Runge-Kutta methods, are also well studied in literature.
However, there are important examples of objects that are neither a group nor a semigroup. One important case is the class of objects closed under a 'sandwich-type' product, (a, b) → aba. For example, the collection of all symmetric positive definite matrices and all selfadjoint Runge-Kutta methods. The sandwich-type composition aba for numerical integrators was studied at length in [17] and references therein. However, if inverses are well defined, we may wish to replace the sandwich product with the algebraically nicer symmetric product (a, b) → ab −1 a. Spaces closed under such products are called symmetric spaces and are the objects of study in this paper. There is a parallel between the theory of Lie groups and that of symmetric spaces. For Lie groups, fundamental tools are the Lie algebra (tangent space at the identity, closed under commutation) and the exponential map from the Lie algebra to the Lie group. In the theory of symmetric spaces there is a similar notion of tangent space. The resulting object is called a Lie triple system (LTS), and is closed under double commutators, [X, [Y, Z]]. Also in this case, the exponential maps the LTS into the symmetric space.
An important decomposition theorem is associated with symmetric spaces and Lie triple systems: Lie algebras can be decomposed into a direct sum of a LTS and a subalgebra. The well known splitting of a matrix as a sum of a symmetric and a skew symmetric matrix is an example of such a decomposition, the skew symmetric matrices are closed under commutators, while the symmetric matrices are closed under double commutators. Similarly, at the group level, there are decompositions of Lie groups into a product of a symmetric space and a Lie subgroup. The matrix polar decomposition, where a matrix is written as the product of a symmetric positive definite matrix and an orthogonal matrix is one example.
In this paper, we are concerned with the application of such structures to the numerical analysis of differential equations of evolution. The paper is organised as follows: in §2 we review some general theory of symmetric spaces and Lie triple systems. Applications of this theory in numerical analysis of differential equations are discussed in §3, which, in turn, can be divided into two parts. In the first ( §3.1- §3.3), we review and discuss the case of differential equations on matrix groups. The properties of these decompositions and numerical algorithms based on them are studied in a series of papers. Polar-type decompositions are studied in detail in [21], with special emphasis on optimal approximation results. The paper [32] is concerned with important recurrence relations for polar-type decompositions, similar to the Baker-Campbell-Hausdorff formula for Lie groups, while [34,10,12] employ this theory to reduce the implementation costs of numerical methods on Lie groups [9,20]. We mention that polar-type decompositions are also closely related to the more special root-space decomposition employed in numerical integrators for differential equations on Lie groups in [23]. In [13] it is shown that the generalized polar decompositions can be employed in cases where the theory of [23] cannot be used.
In the second part of §3 ( §3.4 and beyond), we will consider the application of this theory to numerical methods for the solution of differential equations with symmetries and reversing symmetries. By backward error analysis, numerical methods can be thought of as exact flows of nearby vector fields. The main goal is then to remove from the numerical method the undesired part of the error (either the one destroying the symmetry or the reversing symmetry). These error terms in the modified vector field generally consist of complicated derivatives of the vector field itself and are not explicitly calculated, they are just used formally for analysis of the method. In this context, the main tools are compositions at the group level, using the flow of the numerical method and the symmetries/reversing symmetries, together with their inverses. There is a substantial difference between preserving reversing symmetries and symmetries for a numerical method: the first can be always be attained by a finite (2 steps) composition (Scovel projection [27]), the second requires in general an infinite composition. Thus symmetries are generally more difficult to retain than reversing symmetries. For the retention of symmetry, we review the Thue-Morse symmetrization technique for arbitrary methods and present a new Yoshida-like symmetry retention technique for selfadjoint methods. The latter has always positive intermediate step sizes and can be of interest in the context of stiff problems, which typically require step size restrictions. We illustrate the use of these symmetrisation methods by some numerical experiments.
Finally, Section 4 is devoted to some concluding remarks.
2 General theory of symmetric spaces and Lie triple systems In this section we present some background theory for symmetric spaces and Lie triple systems. We expect the reader to be familiar with some basic concepts of differential geometry, like manifolds, vector fields, etc. For a more detailed treatment of symmetric spaces we refer the reader to [4] and [14] which also constitute the main reference of the material presented in this section. We shall also follow (unless otherwise mentioned) the notational convention of [4]: in particular, M is a set (manifold), the letter G is reserved for groups and Lie groups, gothic letters denote Lie algebras and Lie triple systems, latin lowercase letters denote Lie-group elements and latin uppercase letters denote Lie-algebra elements. The identity element of a group will usually be denoted by e and the identity mapping by id.

Symmetric spaces
Definition 1 (See [14]) A symmetric space is a manifold M with a differentiable symmetric product · obeying the following conditions: and moreover (iv) every x has a neighbourhood U such that for all y in U x · y = y implies y = x.
The latter condition is relevant in the case of manifolds M with open set topology (as in the case of Lie groups) and can be disregarded for sets M with discrete topology: a discrete set M endowed with a multiplication obeying (i)-(iii) will be also called a symmetric space.
A pointed symmetric space is a pair (M, o) consisting of a symmetric space M and a point o called base point. Note that when M is a Lie group, it is usual to set o = e. Moreover, if the group is a matrix group with the usual matrix multiplication, it is usual to set e = I (identity matrix).
The left multiplication with an element x ∈ M is denoted by S x , and is called symmetry around x. Note that S x x = x because of (i), hence x is fixed point of S x and it is isolated because of (iv). Furthermore, (ii) and (iii) imply S x is an involutive automorphism of M, i.e. S 2 x = id. Symmetric spaces can be constructed in several different ways, the following are important examples: 1. Manifolds with an intrinsically defined symmetric product. As an example, consider the n-sphere as the set of unit vectors in R n+1 . The product x · y = S x y = (2xx T − I)y turns this into a symmetric space. The above operation is the reflection of points on a sphere. This can be generalized to m-dimensional subspaces in R n (m ≤ n) in the following manner: Assume that x = [x 1 , . . . , x m ] is a full rank matrix. Define P x the orthogonal projection operator onto the range of x, P x = x(x * x) −1 x * . Consider the reflection R x = 2P x − I. Define x · y = R x y. This operation obeys the conditions (i)-(vi) whenever x, y, z are m × n full rank matrices. In particular, note that (i) is equivalent to R 2 x = I, i.e. the reflection is an involutive matrix. 2. Subsets of a continuous (or discrete) group G that are closed under the composition x · y = xy −1 x, where xy is the usual multiplication in G. Groups themselves, continuous, as in the case of Lie groups, or discrete, are thus particular instances of symmetric spaces. As another example, consider the set of all symmetric positive definite matrices as a subset of all nonsingular matrices, which forms a symmetric space with the product 3. Symmetric elements of automorphisms on a group. An automorphism on a group G is a map σ : G → G satisfying σ (ab) = σ (a)σ (b). The symmetric elements are defined as It is easily verified that A obeys (i)-(iv) when endowed with the multiplication x · y = xy −1 x, hence it is a symmetric space. As an example, symmetric matrices are symmetric elements under the matrix automorphism σ (a) = a −T . 4. Homogeneous manifolds. Given a Lie group G and a subgroup H, a homogeneous manifold M = G/H is the set of all left cosets of H in G. Not every homogeneous manifold possesses a product turning it into a symmetric space, however, we will see in Theorem 1 that any connected symmetric space arises in a natural manner as a homogeneous manifold.
5. Jordan algebras. Let a be a finite-dimensional vector space with a bilinear multiplication 1 x * y such that x * y = y * x, x * (x 2 * y) = x 2 * (x * y) (powers defined in the usual way, x m = x * x m−1 ), with unit element e. Define L x (y) = x * y and set P x = 2L 2 x − L x 2 . Then, the set of invertible elements of a is a symmetric space with the product x · y = P x (y −1 ).
In the context of symmetric matrices, take x * y = 1 2 (xy + yx), where xy denotes the usual matrix multiplication. After some algebraic manipulations, one can verify that the product Let G be a connected Lie group and let σ be an analytic involutive automorphism, i.e. σ = id and σ 2 = id. Let G σ denote fixσ = {g ∈ G : σ (g) = g}, G σ e its connected component including the base point, in this case the identity element e and finally let K be a closed subgroup such that G σ

Theorem 1 ([14])
The homogeneous space M = G/K is a symmetric space with the product xK · yK = xσ (x) −1 σ (y)K and G σ is a symmetric space with the product x · y = xy −1 x. The space of symmetric elements G σ is isomorphic to the homogeneous space G/G σ . Moreover, every connected symmetric space is of the type G/K and also of the type G σ .
The interesting consequence of the above theorem is that every connected symmetric space is also a homogeneous space, which implies a factorization: as coset representatives for G/G σ one may choose elements of G σ , thus any x ∈ G can be decomposed in a product x = pk, where p ∈ G σ and k ∈ G σ . In other words, The matrix polar decomposition is a particular example, discussed in §3.1.
The automorphism σ on G induces an automorphism on the Lie algebra g and also a canonical decomposition of g. Let g and k denote the Lie algebras of G and K respectively and denote by dσ the differential of σ at e, σ (exp(tX)), ∀X ∈ g. (3) Note that dσ is an involutive automorphism of g and has eigenvalues ±1. Moreover, X ∈ k implies dσ (X) = X. Set p = {X ∈ g : dσ (X) = −X}. Then, [4]. It is easily verified that that is, k is a subalgebra of g while p is an ideal in k. Given X ∈ g, its canonical decomposition p ⊕ k is X = P + K, with P ∈ p and K ∈ k, X = P + K, dσ (K) = K, dσ (P) = −P, (algebra splitting).
We have already observed that there is close connection between projection matrices, reflections (involutive matrices) and hence symmetric spaces. In a linear algebra context, this statement can be formalized as follows. Recall that a matrix Π is a projection if Π 2 = Π .
Note that if S is involutive, so is −S, which corresponds to the opposite identification of the ±1 eigenspaces. A matrix K, whose columns are in the span of the +1 eigenspace is said to be block-diagonal with respect to the automorphism, while a matrix P, whose columns are in the span of the −1 eigenspace, is said to be 2-cyclic.
In the context of Lemma 1, we recognize that P ∈ p is the 2-cyclic part and K ∈ k is the block-diagonal part. Namely, if X is represented by a matrix, then where X i j = Π i S XΠ j S restricted to the appropriate subspaces. Then, X −− and X ++ corresponds to the block-diagonal part K and X −+ , X +− corresponding to the 2-cyclic part P.
In passing, we mention that the decomposition (5) is called a Cartan decomposition whenever the Cartan-Killing form B(X,Y ) = tr(ad X ad Y ) is nondegenerate, hence it can be used to introduce a positive bilinear form B dσ = −B(X, dσ (Y )). If this is the case, the linear subspaces k,p are orthogonal.
The involutive automorphism σ need not be defined at the group level G and thereafter lifted to the algebra by (3). It is possible to proceed the other way around: an involutive algebra automorphism dσ on g, which automatically produces a decomposition (4)-(5), can be used to induce a group automorphism σ by the relation and a corresponding group factorization (2). Thus, we have an "upstairs-downstairs" viewpoint: the group involutive automorphisms generate corresponding algebra automorphisms and vice versa. This "upstairs-downstairs" view is useful: in some cases, the group factorization (2) is difficult to compute starting from x and σ , while the splitting at the algebra level might be easy to compute from X and dσ . In other cases, it might be the other way around.

Lie triple systems
In Lie group theory Lie algebras are important since they describe infinitesimally the structure of the tangent space at the identity. Similarly, Lie triple systems give the structure of the tangent space of a symmetric space.
is called a Lie triple system (LTS) if the following identities are satisfied: A typical way to construct a LTS is by means of an involutive automorphism of a Lie algebra g. With the same notation as above, the set p is a LTS with the composition Vice versa, for every LTS there exists a Lie algebra G and an involutive automorphism σ such that the given LTS corresponds to p. The algebra G is called standard embedding of the LTS. In general, any subset of g that is closed under the operator is a Lie triple system. It can be shown that being closed under T X guarantees being closed under the triple commutator.

The classical polar decomposition of matrices
Let GL(N) be the group of N × N invertible matrices. Consider the map It is clear that σ is an involutive automorphism of GL(N). Then, according to Theorem 1, the set of symmetric elements G σ = {x ∈ GL(N) : σ (x) = x −1 } is a symmetric space. We observe that G σ is the set of invertible symmetric matrices. The symmetric space G σ is disconnected and particular mention deserves its connected component containing the identity matrix I, since it reduces to the set of symmetric positive definite matrices. The subgroup G σ consists of all orthogonal matrices. The decomposition (2) is the classical polar decomposition, any nonsingular matrix can be written as a product of a symmetric matrix and an orthogonal matrix. If we restrict the symmetric matrix to the symmetric positive definite (spd) matrices, then the decomposition is unique. In standard notation, p is denoted by s (spd matrix), while k is denoted by q (orthogonal matrix). 2 At the algebra level, the cor- , is the classical algebra of skew-symmetric matrices, while p = {X ∈ gl(N) : dσ (X) = −X} is the classical set of symmetric matrices. The latter is not a subalgebra of gl(N) but is closed under T X , hence is a Lie triple system. The decomposition (6) is nothing else than the canonical decomposition of a matrix into its skew-symmetric and symmetric part, It is well known that the polar decomposition x = sq can be characterized in terms of best approximation properties. The orthogonal part q is the best orthogonal approximation of x in any orthogonally invariant norm (e.g. 2-norm and Frobenius norm). Other classical polar decompositions x = sq, with s Hermitian and q unitary, or with s real and q coinvolutory (i.e. qq = I), can also be fitted in this framework [2] with the choice of automorphisms σ (x) = x − * =x −T (Hermitian adjoint), and σ (x) =x respectively (wherex denotes the complex conjugate of x).
The group decomposition x = sq can also be studied via the algebra decomposition.

Generalized polar decompositions
In [21] such decompositions are generalized to arbitrary involutive automorphisms, and best approximation properties are established for the general case.
In [32] an explicit recurrence is given, if exp(X) = x, exp(S) = s and exp(Q) = q then S and Q can be expressed in terms of commutators of P and K. The first terms in the expansions of S and Q are Clearly, also other types of automorphisms can be considered, generalizing the group factorization (2) and the algebra splitting (6). For instance, there is a large source of involutive automorphisms in the set of involutive inner automorphisms that can be applied to subgroups of G = GL(n) to obtain a number of interesting factorizations. The matrix r has to be involutive, r 2 = I, but it need not be in the group: the factorization makes sense as long as σ (x) is in the group (resp. dσ (X) is in the algebra). As an example, let G = SO(n) be the group of orthogonal matrices, and let r = [−e 1 , e 2 , . . . , e n ] = I − 2e 1 e T 1 , where e i denotes the ith canonical unit vector in R n . Obviously, r ∈ SO(n), as det r = −1; nevertheless, we have (rxr) T (rxr) = r T x T r T rxr = I, as long as x ∈ SO(n), since r = r T and r T r = I, thus σ (x) = rxr ∈ SO(n). It is straightforward to verify that the subgroup G σ of Theorem 1 consists of all orthogonal n × n matrices of the form where q n−1 ∈ SO(n−1). Thus the corresponding symmetric space is G/G σ = SO(n)/SO(n− 1). Matrices belong to the same coset if their first column coincide, thus the symmetric space can be identified with the (n − 1)-sphere S n−1 .
The corresponding splitting of a skew-symmetric matrix V ∈ g = so(n) is . , e n , −e n+1 , e n+2 . . . , e 2n ] for symplectic matrices, gives the splitting in lower-dimensional symplectic matrices forming a sub-algebra and a Lietriple system, and so on. In [34,10], such splittings are used for the efficient approximation of the exponential of skew-symmetric, symplectic and zero-trace matrices. In [12] similar ideas are used to construct computationally effective numerical integrators for differential equations on Stiefel and Grassman manifolds.

Generalized polar coordinates on Lie groups
A similar framework can be used to obtain coordinates on Lie groups [13], which are of interest when solving differential equations of the typė by reducing the problem recursively to spaces of smaller dimension. Recall that where Ω obeys the differential equationΩ = dexp −1 j! ad j A , B j being the jth Bernoulli number, see [9].
, the solution can be factorized as Note that (13) depends on P, however, it is possible to formulate it solely in terms of K, Π − X and Π + X, but the expression becomes less neat. In block form: where v = ad P . The above formula paves the road for a recursive decomposition, by recognizing that dexp −1 K is the dexp −1 function on the restricted sub-algebra k. By introducing a sequence of involutive automorphisms σ i , one induces a sequence of subalgebras, g = g 0 ⊃ g 1 ⊃ g 2 ⊃ . . ., of decreasing dimention, g i+1 = Range(Π + σ i ). Note also that the functions appearing in the above formulation are all analytic functions of the ad-operator, and are either odd or even functions, therefore they can be expressed as functions of ad 2 on p. In particular, this means that, as long as we can compute analytic functions of the ad 2 operator, the above decomposition is computable.
Thus, the problem is reduced to the computation of analytic functions of the 2-cyclic part P as well as analytic functions of ad P (trivialized tangent maps and their inverse). The following theorem addresses the computation of such functions using the same framework of Lemma 1. where A similar result holds for ad P , see [13]. It is interesting to remark that if where ψ 1 , ψ 2 as above. In particular, the problem is reduced to computing analytic functions of the principal square root of a matrix [7]: numerical methods that compute these quantities accurately and efficiently are very important for competitive numerical algorithms. Typically, Θ is a low-rank matrix, hence computations can be done using eigenvalues and eigenfunctions of the ad operator restricted to the appropriate space, see [12,13]. In particular, if A, B are vectors, then B T A is a scalar and the formulas become particularly simple. These coordinates have interesting applications in control theory. Some early use of these generalized Cartan decompositions (4)-(5) (which the author calls Z 2 -grading) to problems with nonholonomic constraints can be found in [1]. In [11], the authors embrace the formalism of symmetric spaces and use orthogonal (Cartan) decompositions with applications to NMR spectroscopy and quantum computing, using adjoint orbits as main tool. Generally, these decompositions can be found in the literature, but have been applied mostly to the cases when [k, k] = {0} or [p, p] = {0}, or both, as (12)-(13) become very simple and the infinite sums reduce to one or two terms. The main contribution of [13] is the derivation of such differential equations and the evidence that such equations can be solved efficiently using linear algebra tools. See also [33] for some applications to control theory.
we say that R is a reversing symmetry of ϕ [17]. Without further ado, we restrict ourselves to involutory symmetries, the main subject of this paper. Symmetries and reversing symmetries are very important in the context of dynamical systems and their numerical integration. For instance, nongeneric bifurcations can become generic in the presence of symmetries and vice versa. Thus, when using the integration time-step as a bifurcation parameter, it is vitally important to remain within the smallest possible class of systems. Reversing symmetries, on the other hand, give rise to the existence of invariant tori and invariant cylinders [19,25,28,29].
It is a classical result that the set of symmetries possess the structure of a group -they behave like automorphisms and fixed sets of automorphisms. The group structure, however, does not extend to reversing symmetries and fixed points of anti-automorphisms, and in the last few years the set of reversing symmetries has received the attention of numerous numerical analysts. In [17] it was observed that the set of fixed points of an involutive antiautomorphism A − was closed under the operation ϕ 1 · ϕ 2 = ϕ 1 ϕ 2 ϕ 1 ∈ fixA − , ∀ϕ 1 , ϕ 2 ∈ fixA − , that McLachlan et al. called "sandwich product". 3 Indeed, our initial goal was to understand such structures and investigate how they could be used to devise new numerical integrators for differential equations with some special geometric properties. We recognise, cfr. §2.1, that the set of fixed points of an anti-automorphism is a symmetric space. Conversely, any connected space of invertible elements closed under the "sandwich product", is the set of the fixed points of an involutive automorphism (cfr. Theorem 1) and has associated to it a LTS. To show this, consider the well known symmetric BCH formula, [26], which is used extensively in the context of splitting methods [16]. Because of the sandwich-type composition (symmetric space structure), the corresponding Z must be in the LTS space, and this explains why it can be written as powers of the double commutator operators T X = ad 2 X , T Y = ad 2 Y applied to X and Y . A natural question to ask is: what is the automorphism σ having such sandwich-type composition as anti-fixed points? As fixA − = {z|σ (z) = z −1 }, we see that, by writing z = exp Z, we have σ (exp(Z)) = (exp(Z)) −1 = exp(−Z). In the context of numerical integrators, the automorphism σ consists in changing the time t to −t. This will be proved in §3.5.
If M is a finite dimensional smooth compact manifold, it is well known that the infinite dimensional group of Diff(M) of all smooth diffeomorphisms M → M is a Lie group, with Lie algebra Vect(M) of all smooth vector fields on M, with the usual bracket and exponential map. It should be noted, however, that the exponential map is not a one-to-one map, not even in the neighbourhood the identity element, since there exist diffeomorphisms arbitrary close to the identity which are not on any one-parameter subgroup and others which are on many. However, the regions where the exponential map is not surjective become smaller and smaller the closer we approach the identity [24,22], and, for our purpose, we can disregard these regions and assume that our results are formally true.
There are two different settings that we can consider in this context. The first is to analyze the set of differentiable maps that possess a certain symmetry (or a discrete set of symmetries). The second is to consider the structure of the set of symmetries of a fixed diffeomorphism. The first has a continuous-type structure while the second is more often a discrete type symmetric space.

Proposition 1
The set of diffeomorphisms ϕ that possess R as an (involutive) reversing symmetry is a symmetric space of the type G σ .
It is clear that σ acts as an automorphism, σ (ϕ 1 ϕ 2 ) = σ (ϕ 1 )σ (ϕ 2 ), moreover, if R is an involution then so is also σ . Note that the set of diffeomorphisms ϕ that possess R as a reversing symmetry is the space of symmetric elements G σ defined by the automorphism σ (cfr. §2). Hence the result follows from Theorem 1.

Proposition 2
The set of reversing symmetries acting on a diffeomorphism ϕ is a symmetric space with the composition R 1 · R 2 = R 1 R −1 2 R 1 .
Proof If R 1 is a symmetry of ϕ then so is also R −1 , since R −1 1 ϕ −1 R 1 = ϕ and the assertion follows by taking the inverse on both sides of the equality. In particular, if R 1 is a symmetry of ϕ it is also true that R −1 1 is a reversing symmetry of ϕ −1 . Next, we observe that if R 1 and R 2 are two reversing symmetries of ϕ then so is also R 1 S −1 R 1 , since It follows that the composition R 1 · R 2 = R 1 R −1 2 R 1 is an internal operation on the set of reversing symmetries of a diffeomorphism ϕ.
With the above multiplication, the conditions i)-iii) of Definition 1 are easily verified. This proves the assertion in the case when φ has a discrete set of reversing symmetries.
Acting on ϕ = exp(tX), we have where T * is the tangent map of T . The pullback is natural with respect to the Jacobi bracket, for all vector fields X,Y . Hence the map dσ is an involutory algebra automorphism. Let k σ and p be the eigenspaces of dσ in g = diff(M). Then is the Lie algebra of vector fields that have T as a symmetry and is the Lie triple system, vector fields corresponding to maps that have T as a reversing symmetry. Thus, as is the case for matrices, every vector field X can be split into two parts, having T as a symmetry and reversing symmetry respectively. In the context of ordinary differential equation, let us consider Given an arbitrary involutive function T , the vector field F can always be canonically split into two components, having R as a symmetry and reversing symmetry respectively. However, if one of these components equals zero, then the system (16) has T as a symmetry or a reversing symmetry.

Selfadjoint numerical schemes as a symmetric space
Let us consider the ODE (16), whose exact flow will be denoted as ϕ = exp(tF). Backward error analysis for ODEs implies that a (consistent) numerical method for the integration of (16) can be interpreted as the sampling at t = h of the flow ϕ h (t) of a vector field F h (the so called modified vector field) which is close to F, where p is the order of the method (note that setting t = h, the local truncation error is of order h p+1 ). Consider next the map σ on the set of flows depending on the parameter h defined as where The map σ is involutive, since σ 2 = id, and it is easily verified by means of the BCH formula that σ (ϕ 1,h ϕ 2,h ) = σ (ϕ 1,h )σ (ϕ 2,h ), hence σ is an automorphism. Consider next , namely the method ϕ h is selfadjoint. Proposition 3 The set of one-parameter, consistent, selfadjoint numerical schemes is a symmetric space in the sense of Theorem 1, generated by σ as in (17).
Next, we perform the decomposition (4). We deduce from (17) that is the subalgebra of vector fields that are odd in h, and is the LTS of vector fields that possess only even powers of h. Thus, if F h is the modified vector field of a numerical integrator ϕ h , its canonical decomposition in k ⊕ p is the first term containing only odd powers of h and the second only even powers. Then, if the numerical method ϕ h (h) is selfadjoint, it contains only odd powers of h locally (in perfect agreement with classical results on selfadjoint methods [3]).
3.6 Connections with the generalized Scovel projection for differential equations with reversing symmetries In [21] it has been shown that it is possible to generalize the polar decomposition of matrices to Lie groups endowed with an involutive automorphism: every Lie group element z sufficiently close to the identity can be decomposed as z = xy where x ∈ G σ , the space of symmetric elements of σ , and y ∈ G σ , the subgroup of G of elements fixed under σ . Furthermore, setting z = exp(tZ) and y = exp(Y (t)), one has that Y (t) is an odd function of t and it is a best approximant to z in G σ in G σ right-invariant norms constructed by means of the Cartan-Killing form, provided that G is semisimple and that the decomposition g = p ⊕ k is a Cartan decomposition. Assume that ϕ, the exact flow of the differential equation (16), has R as a reversing symmetry (i.e. F ∈ p σ , where σ (ϕ) = RϕR −1 ), while its approximation ϕ h has not. We perform the polar decomposition i.e. ψ h has R as a reversing symmetry, while χ h has R as a symmetry. Since the original flow has R as a reversing symmetry (and not symmetry), χ h is the factor that we wish to eliminate. We have ψ 2 h = ϕ h σ (ϕ h ) −1 . Hence the method obtained composing ϕ h with σ (ϕ h ) −1 has the reversing symmetry R every other step. To obtain ψ h we need to extract the square root of the flow ϕ h σ (ϕ h ) −1 . Now, if φ (t) is a flow, then its square root is simply φ (t/2). However, if φ h (t) is the flow of a consistent numerical method (p ≥ 1), namely the numerical integrator corresponds to φ h (h), it is not possible to evaluate the square root φ h (h/2) by simple means as it is not the same as the numerical method with half the stepsize, φ h/2 (h/2). The latter, however, offers an approximation to the square root: note that an expansion which, compared with φ h (h), reveals that the error in approximating the square root with the numerical method with half the stepsize is of the order of a term that is subsumed in the local truncation error. The choicẽ as an approximation to ψ h (we stress that each flow is now evaluated at t = h/2), yields a map that has the reversing symmetry R at each step by design, since is the inverse (or adjoint) method of ϕ h/2 . If σ is given by (17), then σ (ϕ −1 h/2 ) = ϕ * h/2 (h/2) and this algorithm is precisely the Scovel projection [27] originally proposed to to generate selfadjoint numerical schemes from an arbitrary integrator, and then generalized to the context of reversing symmetries [17].

Proposition 4
The generalized Scovel projection is equivalent to choosing the G σ -factor in the polar decomposition of a flow ϕ h under the involutive automorphism σ (ϕ) = RϕR −1 , whereby square roots of flows are approximated by numerical methods with half the stepsize.

Connection with the Thue-Morse sequence and differential equation with symmetries
Another algorithm that can be related to the generalized polar decomposition of flows is the application of the Thue-Morse sequence to improve the preservation of symmetries by means of a numerical integrator [8]. Given an involutive automorphism σ and a numerical method ϕ h in a group G of numerical integrators, Iserles et al. [8] construct the sequence of methods Since , it is easily observed that the k-th method corresponds to composing σ 0 ϕ [k] and σ 1 ϕ [k] according to the k-th Thue-Morse sequence 01101001 . . ., as displayed below in Table 1 (see [30,18]). Iserles et al. ( [8]) showed, by a recursive use of the BCH formula, that each iteration improves by one order the preservation of the symmetry S by a consistent numerical method, where S is the involutive automorphism such that σ (φ ) = S φ S −1 . The main argument of the proof is that if the method φ [k] has a given symmetry error, the symmetry error of σ (φ [k] ) has the opposite sign. Hence the two leading symmetry errors cancel and φ [k+1] has a symmetry error of one order higher. In other words, if the method ϕ h preserves S to order p, then ϕ [k] preserves the symmetry S to order p + k every 2 k steps. As σ changes the sign of the symmetry error only, if a method ϕ h has a given approximation order p, so does σ (ϕ h ). Thus φ [0] = σ (ϕ h ) can be used as initial condition in (20), obtaining the conjugate Thue-Morse sequence 10010110 . . .. By a similar argument, also generate sequences with increasing order of preservation of symmetry. u 0 (y, x) = u 0 (x, y), so is the solution for all t. The symmetry is σ (u(x, y)) = u(y, x). Now, assume that we solve the equation by the method of alternating directions, where the method ϕ corresponds to solving with respect to the x variable keeping y fixed, while σ (ϕ) with respect to the y variable keeping x fixed. The symmetry will typically be broken at the first step. Nevertheless, we can get a much more symmetric solution if the sequence of the directions obeys the Thue-Morse sequence (iteration ϕ [k] ) or the equivalent sequence given by iteration χ [k] . This example is illustrated in Figures 1-2. 3.8 Connections with a Yoshida-type composition and differential equations with symmetries In a famous paper that appeared in 1990 ( [31]) Yoshida showed how to construct high order time-symmetric integrators starting from lower order time-symmetric symplectic ones. Yoshida showed that, if ϕ is a selfadjoint numerical integrator of order 2p, then is a selfadjoint numerical method of order 2p + 2 provided that the coefficients α and β satisfy the condition for the PDE u t = u x + u y + 2 × 10 −3 u 2 , on a square, with initial condition u 0 = e −x 2 −y 2 and periodic boundary conditions. The integration is performed using the first order method u k+1 = u k + h(Du k + 1 2 f (u k )) (Forward Euler), where D is a circulant divided difference discretization matrix of the differential operators ∂ x , ∂ y (second order central differences). The experiments are performed with constant stepsize h = 10 −2 . The order of the overall approximation is the same as the original method (first order only), the only difference is the order of the directions, chosen according to the Thue-Morse sequence. The different symbols correspond to different sampling rates: every step, every second step, fourth, eight, . . . . From top left to bottom right: sequences '0', '01', '0110', '01101001', etc. whose only real solution is In the formalism of this paper, time-symmetric methods correspond to G σ -type elements with σ as in (17) and it is clearly seen that the Yoshida technique can be used in general to improve the order of approximation of G σ -type elements. A similar procedure can be applied to improve the order of the retention of symmetries and not just reversing symmetries. To be more specific, let S be a symmetry of the given differential equation, namely S * F = FS , with S = id, S −1 = S , S * denoting the pullback of S to g = Vect(M) (see §3.4). Here, the involutive automorphism is given by    Fig. 2 Global error and symmetry error max x,y |u((t, x, y) − u(t, y, x)| versus step size h for the method of alternating directions for the PDE u t = u x + u y + 2 × 10 −3 u 2 , as above. The basic method ϕ (method "0") is now a symmetric composition with the Heun method: half step in the x direction, full step in the y direction, half step in the x direction, all of them performed with the Heun method (improved Euler). It is immediately observed that the symmetry error is improved by choosing the directions according to the Thue-Morse sequence. The order of the method (global error) remains unchanged (order two).
Proposition 5 Assume that ϕ h (t) is the flow of a self-adjoint numerical method of order 2p, , where E j = P j + K j , P j ∈ p,K j ∈ k and F has S as a symmetry. The composition ϕ [1] h (t) = ϕ ah (at)σ (ϕ bh (bt))ϕ ah (at), has symmetry error 2p + 2 at t = h.
Application of the symmetric BCH formula, together with the fact that dσ acts by changing the signs on the p-components only, allows us to write the relation (23) as ϕ [1] h (t) = exp((2a + b)tF + (2(at)(ah) 2p + (bt)(bh) 2p )K 2p (24) where the O th 2p+2 comes from the E 2p+2 term and the O t 3 h 2p from the commutation of the F and E 2p terms (recall that no first order commutator appears in the symmetric BCH formula). The numerical method is obtained letting t = h. We require 2a + b = 1 for consistency, and 2a 2p+1 − b 2p+1 = 0 to annihilate the coefficient of P 2p , the lowest order p-term. The resulting method ϕ [1] h (t) retains the symmetry S to order 2p + 2, as the first leading symmetry error is a O h 2p+3 term.
This procedure allows us to gain two extra degrees in the retention of symmetry per iteration, provided that the underlying method is selfadjoint, compared with the Thue-Morse sequence of [8] that yields one extra degree in symmetry per iteration but does not require selfadjointness. As for the Yoshida technique, the composition (23) can be iterated k times to obtain a time-symmetric method of order 2p that retains symmetry to order 2(p + k). The disadvantage of (23), with respect to the classical Yoshida approach is the fact that the order of the method is retained (and does not increase by 2 units as the symmetry error). The main advantage is that all the steps are positive, in particular the second step, b, whereas the β is always negative for p ≥ 3 and typically larger than α, requiring step-size restrictions for stiff methods. In the limit, when p → ∞, a, b → 1/3, i.e. the step lengths become equal. Thus, the proposed technique for improving symmetry is of particular interest in the context of stiff problems. This is illustrated by the following example and Figure 3.
with a gaussian initial condition, u(0) = e −9x 2 −9y 2 , and periodic boundary conditions. The problem is semi-discretized on a uniform and isotropic mesh with spacing δ = 0.1 and is reduced to the set of ODEsu = (D xx + D yy )u + f (u) = F(u), where D xx = I ⊗ D 2 and D yy = D 2 ⊗ I, D 2 being the circulant matrix of the standard second order divided differences, with stencil 1 δ 2 [1, −2, 1]. For the time integration, we consider the splitting F = F 1 + F 2 , where F 1 = D xx + 1 2 f and F 2 = D yy + 1 2 f and the second order self-adjoint method: 4 We display the global error and the symmetry error for time-integration step sizes h = 3h 0 × [1, 1 2 , . . . , 1 64 ], where the parameter h 0 is chosen to be the largest step size for which the basic method ϕ in (25) is stable. The factor 3 comes from the fact that both the Yoshida and our symmetrization method (23) require 3 sub-steps of the basic method. So, one step of the Yoshida and our symmetrising composition can be expected to cost the same as the basic method (25).
As we can see in Figure 3, the Yoshida technique, with α = 1/(2−2 1/3 ) and β = 1−2α, does the job of increasing the order of accuracy of the method from two to four, and so does the symmetry error. However, since α < 1/2, the β -step is negative, and, as a consequence, it is observed that the method fails to converge for the two largest values of the step size. Conversely, our symmetrization method (23) has a = 1/(2 + 2 1/3 ) and b = 1 − 2a, with b positive, and the method converges for all the time-integration steps. As expected, the order is not improved, but the symmetry order is improved by two units. The symmetry σ is applied by transposing the matrix representation u i, j of the solution at (iδ , jδ ), before and after the intermediate b-step. Otherwise, the two implementations are identical.

Conclusions and remarks
In this paper we have shown that the algebraic structure of Lie triple systems and the factorization properties of symmetric spaces can be used as a tool to: 1) understand and provide a unifying approach to the analysis of a number of different algorithms; and 2) devise new algorithms with special symmetry/reversing symmetry properties in the context of the numerical solution of differential equations. In particular, we have seen that symmetries are more difficult to retain (to date, we are not aware of methods that can retain a generic involutive symmetry in a finite number of steps), while the situation is simpler for reversing symmetries, which can be achieved in a finite number of steps using the Scovel composition. So far, we have considered the most generic setting where the only allowed operations are noncommutative compositions of maps (the map ϕ, its transformed, σ (ϕ), and their inverses). If the underlying space is linear and so is the symmetry, i.e. σ (ϕ 1 + ϕ 2 ) = σ (ϕ 1 ) + σ (ϕ 2 ), the mapφ h = ϕ h +σ (ϕ h ) 2 obviously satisfies the symmetry σ , as σ (φ h ) =φ h . Because of the linearity and the vector-space property, we can use the same operation as in the tangent space, namely we identify σ and dσ . This is in fact the most common symmetrization procedure for linear symmetries in linear spaces. For instance, in the context of the alternating-direction examples, a common way to resolve the symmetry issue is to first solve, say, the x and the y direction, solve the y and the x direction with the same initial condition, and then average the two results.
In this paper we did not mention the use and the development of similar concepts in a more strict linear algebra setting. 5 Some recent works deal with the further understanding of scalar products and structured factorizations, and more general computation of matrix functions preserving group structures, see for instance [15,5,6] and references therein. Some of these topics are covered by other contributions in the present BIT issue, which we strongly encourage the reader to read to get a more complete picture of the topic and its applications.