Methods for the approximation of the matrix exponential in a Lie-algebraic setting

Elena Celledoni^{∗}, MSRI, 1000 Centennial Drive, Berkeley CA 94720, celledon@msri.org,

Arieh Iserles, DAMTP, Cambridge University, Silver Street, CB3 9EW, Cambridge, England ai@damtp.cam.ac.uk.

February 25, 1999

Abstract

Discretization methods for ordinary differential equations based on the use of matrix exponentials have been known for decades. This set of ideas has come off age and acquired greater urgency recently, within the context of geometric integration and discretization methods on manifolds based on the use of Lie-group actions.

In the present paper we study the approximation of the matrix exponential in a particular context: given a Lie group G and its Lie algebra 𝔤, we seek approximants F (tB) of exp(tB) such that F (tB) ∈ G if B ∈ 𝔤. Having fixed a basis V_{1}, …, V_{d} of 𝔤, we write F (tB) as a composition of exponentials of the type exp(α_{i}(t)V_{i}), where α_{i} for i = 1, 2, …, d are scalar functions. In this manner it becomes possible to increase the order of the approximation without increasing the number of exponentials to evaluate and multiply together. We study order conditions and implementation details and conclude the paper with some numerical experiments.

1 Introduction

Although numerical methods for the integration of ordinary differential equations (ODEs) based on the use of the matrix exponential have long history, the subject has acquired new relevance recently with two developments. The first, which is irrelevant to the theme of this paper, is the introduction of Krylov subspace techniques and their application to large stiff systems of differential equations (Hochbruck, Lubich & Selhofer 1998). The other development is motivated by the philosophy of geometric integration and its purpose is to recover under discretization important qualitative and geometric features of the underlying dynamical system. Examples of such methods can be found inter alia in (Casas 1996, Crouch & Grossman 1993). An important technique in geometric integration is the use of Lie-group actions, which lend themselves to the design of very effective time-stepping methods for ODEs evolving on homogeneous manifolds. Such methods have been recently studied in (Munthe-Kaas 1997) and (Engø 1998). Methods based on the use of the classical Magnus and Fer expansions for integrating ODEs on Lie-groups can be brought into this formalism (Iserles & Nørsett 1997, Iserles & Nørsett 1999, Zanna 1997). All such methods require a repeated evaluation of a matrix exponential, often of large matrices. Inasmuch as typically one can expect the replacement of the exact exponential by a suitable approximant (a rational function, say, a Krylov subspace approximant or a Schur factorization), the context of Lie-group methods imposes a crucial extra requirement. The approximant in question, applied to an arbitrary element of the Lie algebra 𝔤, must produce an outcome in the Lie group G, otherwise the whole purpose of the calculation, dicretizing within G, will be null and void. This can be done is some, but by no means, all Lie algebras of interest and we refer the reader to (Celledoni & Iserles 1998) for a more substantive discussion of this issue.

∗Research at MSRI is supported in part by NSF grant DMS-9701755.

Let G be a finite-dimensional Lie group. For all practical purposes, we may assume that G is a subgroup of the general linear group GL(n), the set of all nonsingular n × n matrices. We denote by 𝔤 the Lie algebra corresponding to G, observing that it is a subalgebra of 𝔤𝔩(n), the Lie algebra of all n×n matrices. Our concern in this paper is with differential equations that evolve on a manifold M subject to the action of G. For simplicity we can assume that M coincides with G and the action is of G on itself. The numerical solution of such differential equations can be obtained considering the pull-back on 𝔤, by means of the exponential map, of the vectorfield defining the equation. We can compute the corresponding flow by a Lie-algebra discretization method and recover the approximation of the original problem via exponentiation. Given an integration method of order p, we consider order-p approximants F (tB) for exp(tB), where B ∈ 𝔤 and t ≥ 0. We require that F (tB) ∈ G, whence it is easy to prove that important qualitative features of the original equation and the order of the discretization are retained. In (Celledoni & Iserles 1998) we have introduced low-rank splitting methods for the construction of the approximant F, as the first attempt to provide a comprehensive treatment of this issue.

Although the constraint F (tB) ∈ G represents remarkable advantage in many applications, such as problems in which the conservation of invariants is at issue in numerical modelling (volume conservation in meteorology, invariance under rotations in the theory of mechanical systems and in robotics), it should not be interpreted as the sole purpose of our analysis. Our methods are relevant also for the approximation of exp(tB) in the more general setting B ∈ 𝔤𝔩(n). Suppose in fact that B ∈ 𝔤𝔩(n) and we want to approximate exp(tB). It is always possible to write B as a sum of a matrix B_{s} ∈ 𝔰𝔩(n) (the special linear algebra of n × n matrices with zero trace) and a diagonal matrix B_{d} whose nonzero entries are equal to δ = tr (B)/n. Then B_{s} = B − B_{d} and [B_{s}, B_{d}] = 0 so that exp(tB) = exp(tB_{s}) exp(tδ). This fact is a particular case of what is known in Lie theory as the Levi decomposition (Humphreys 1972, Varadarajan 1984). Using this decomposition of the matrix B, if necessary in tandem with some scaling and squaring technique, the approximation of exp(tB) can be always reduced to the approximation of exp(tB_{s}) with B_{s} ∈ 𝔰𝔩(n). As long as we can assure that our approximation of exp(tB_{s}) resides in SL(n), the outcome is an approximant F (tB) of exp(tB) that shares with the exact exponential the feature that det F (tB) = exp(tr B).

It is possible to prove that, given a splitting B = ∑_{i=1}^{k} B_{i}, the function

e^{\frac{1}{2} tB_{1}}e^{\frac{1}{2} tB_{2}} ⋯ e^{\frac{1}{2} tB_{k−1}}e^{tB_{k}}e^{\frac{1}{2} tB_{k−1}} ⋯ e^{\frac{1}{2} tB_{2}}e^{\frac{1}{2} tB_{1}},

known as the generalized Strang splitting, approximates exp(tB) to order 2. As long as B_{1}, B_{2}, …, B_{k} ∈ 𝔤, it follows at once from the definition of a Lie group that the approximant resides in G. Moreover, 2k − 1 is the least number of exponentials that render such a splitting into a second-order approximant (Celledoni & Iserles 1998). The Strang splitting is time reversible, hence it follows readily from classical theory that the order can be raised from 2 to 4 by composing three Strang splittings with different time steps (Yoshida 1990). In that case we need to evaluate 3k exponentials and multiply 6k matrices. In the case of low-rank splittings which have been considered by Celledoni & Iserles (1998) this results in the following count of flops: 4n^{3} for order 2, 12n^{3} for order 4.

In this paper we present composition methods in which the number of exponentials k equals the dimension d of the Lie algebra. Our construction allows us to increase the order of the approximation without increasing the number of exponentials to evaluate and multiply together. Letting {V_{1}, …, V_{d}} be a basis of 𝔤, we write F (tB) as a composition of exponentials of the type exp(α_{i}(t)V_{i}), where each α_{i}(t) for i = 1, …, d is a scalar function. In general d = O(n^{2}), however, with an appropriate choice of the basis elements, the computation of each exponential exp(α_{i}(t)V_{i}) requires O(n) flops, while the formation of their product adds just 2n^{3} + O(n^{2}) flops. The challenging part of the computation is the construction of the functions α_{1}(t), …, α_{d}(t), and the cost of their calculation depends on the desired order of the approximation. Naive complexity analysis might have indicated that the total cost is growing exponentially in d as the order increases. Yet, the cost remains relatively modest for small orders and the method lends itself very well to the exploitation of sparsity in the matrix B. In the sequel we show how this approach can be turned into an efficient numerical method and we obtain algorithms of order up to 4 with a cost of O(n^{3}) for dense matrices.

Our approach can be interpreted as representing the solution using canonical coordinates of the second kind, an approach that has been pioneered by Owren & Marthinsen (1998) in the context of general Lie-group methods. Having said this, the more restrictive framework of exponential approximants possesses a very great deal of special structure. This can be exploited so as to produce efficient and competitive algorithms that approximate exp(tB), B ∈ 𝔤, in the Lie group G.

2 The technique of coordinates of second kind for the approximation of the exponential matrix

Let G be a Lie group and 𝔤 its corresponding d-dimensional Lie algebra. We choose a basis {V_{1}, …, V_{d}} of 𝔤, whence every element Y ∈ G sufficiently close to the identity can be represented in a unique fashion as

Y = exp(γ_{1}V_{1}) exp(γ_{2}V_{2}) ⋯ exp(γ_{d}V_{d}),

where exp : 𝔤 → G is the exponential map. This representation is known as representation in canonical coordinates of the second kind (Varadarajan 1984). This representation is global in the case of solvable Lie algebras. We restrict ourselves to the case 𝔤 ⊆ 𝔤𝔩(n), G ⊆ GL(n), when exp is the usual matrix exponential.

Given B ∈ 𝔤, we can represent it in a unique fashion as

B = ∑_{i=1}^{d} β_{i}V_{i}.

It is possible then to write exp(tB) in the form

U(t) = exp(tB) = exp(g_{1}(t)V_{1}) exp(g_{2}(t)V_{2}) ⋯ exp(g_{d}(t)V_{d}).

Letting g = [g_{1}, …, g_{d}]^{T}, β = [β_{1}, …, β_{d}]^{T}, it can be proved that the vector function g obeys a differential equation of the form

\frac{dg}{dt} = f(β, g), g(0) = 0,

where f is a suitable function of β and g, for sufficiently small t (Wei & Norman 1963). Given a solvable Lie algebra 𝔤 Wei & Norman (1963) prove results on the global representation of U. However, an explicit form of f is known only for very simple examples of low-dimensional Lie algebras.

In this paper we seek polynomials α_{1} ≈ g_{1}, …, α_{d} ≈ g_{d} of a suitable degree so that

exp(tB) ≈ exp(α_{1}(t)V_{1}) exp(α_{3}(t)V_{2}) ⋯ exp(α_{d}(t)V_{d}).

Differentiation yields

∑_{i=1}^{d} β_{i}V_{i} = ∑_{i=1}^{d} g′_{i}(t) ∏_{j=1}^{i−1} e^{g_{j}V_{j}}V_{i} ∏_{j=i−1}^{1} e^{−g_{j}V_{j}}.  (2.1)

Evaluating this expression at the origin gives the first-order condition

g′_{i}(0) = β_{i},  i = 1, 2, …, d.  (2.2)

Further differentiations of (2.1) lead to higher-order conditions. Let us define the functions

P_{i}(g) = exp(ad_{g_{1}V_{1}}) ∘ ⋯ ∘ exp(ad_{g_{i−1}V_{i−1}})(V_{i}),  i = 1, 2, …, d,  (2.3)

where the adjoint operator ad_{x} : 𝔤 → 𝔤 is defined as ad_{x}(y) = [x, y] for any x, y ∈ 𝔤, [x, y] = xy − yx being the matrix commutator. Note that P_{i}(g(0)) = V_{i}, i = 1, 2, …, d. Moreover, the right-hand side of (2.1) can be written in the simplified form The function

T (g) = ∑_{i=1}^{d} g′_{i}(t)P_{i}(g).

Since the derivatives of the left hand side of (2.1) vanish, the conditions for order p ≥ 1, can be obtained by solving the equations

\frac{d^{r}}{dt^{r}} T (g)|_{t=0} = 0,  r = 1, 2, …, p − 1, p ≥ 1,  (2.4)

where

\frac{d^{r}}{dt^{r}} T (g) = ∑_{i=1}^{d} ∑_{k=1}^{r} \binom{r}{k}\frac{d^{r−k+1}g_{i}}{dt^{r−k+1}} \frac{d^{k}P_{i}}{dt^{k}} .  (2.5)

In particular,

\frac{d}{dt}T (g) = ∑_{i=1}^{d} (g′′_{i}P_{i} + g′_{i} \frac{d}{dt}P_{i}) ,

\frac{d^{2}}{dt^{2}} T (g) = ∑_{i=1}^{d} (g′′′_{i}P_{i} + 2g′′_{i} \frac{d}{dt}P_{i} + g′_{i} \frac{d^{2}}{dt^{2}} P_{i}) ,

\frac{d^{3}}{dt^{3}} T (g) = ∑_{i=1}^{d} (g_{i}^{IV} P_{i} + 3g′′′_{i} \frac{d}{dt}P_{i} + 3g′′_{i} \frac{d^{2}}{dt^{2}} P_{i} + g′_{i} \frac{d^{3}}{dt^{3}} P_{i}) .

Solving (2.4) for r = 1 results in the values of g′′_{i}(0) for i = 1, 2, …, d that allow us to construct an order-2 approximant. Substituting such values in (2.4) for r = 2 yields g′′′_{i}(0) for i = 1, 2, …, d and consequently an approximant of order 3. Similar procedure can be used to construct recursively approximants of arbitrarily high order.

The main part of the computation is the evaluation of the k-th derivative of P_{i}(g) at t = 0. Expanding the exponentials in (2.3) we obtain

P_{i}(g) = ∏_{k=1}^{i−1}(I + ad_{g_{k}V_{k}} + \frac{1}{2}ad_{g_{k}V_{k}}^{2} + \frac{1}{6}ad_{g_{k}V_{k}}^{3} + …)(V_{i})

and, after further algebra,

P_{i}(g) = {I + ∑_{k=1}^{i−1} ad_{g_{k}V_{k}} + ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} ad_{g_{l}V_{l}}ad_{g_{k}V_{k}}

+ \frac{1}{2} ∑_{k=1}^{i−1} ad_{g_{k}V_{k}}^{2} + ∑_{k=3}^{i−1} ∑_{l=2}^{k−1} ∑_{j=1}^{l−1} ad_{g_{j}V_{j}}ad_{g_{l}V_{l}}ad_{g_{k}V_{k}}

+ \frac{1}{2} ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} (ad_{g_{l}V_{l}}ad_{g_{k}V_{k}}^{2} + ad_{g_{l}V_{l}}^{2}ad_{g_{k}V_{k}}) + \frac{1}{6} ∑_{k=1}^{i−1} ad_{g_{k}V_{k}}^{3} + …\} (V_{i}) .

Similarly to (Owren & Marthinsen 1997), we write P_{i}(g) in the form

P_{i}(g) = I + ∑_{r=1}^{∞} ∑_{j_{1}=1}^{i−1} ∑_{j_{2}=j_{1}}^{i−1} ⋯ ∑_{j_{r}=j_{r−1}}^{i−1} \frac{1}{j!}g_{j_{1}} … g_{j_{r}}ad_{V_{j_{1}}} ∘ ⋯ ∘ ad_{V_{jr}} (V_{i}).

Here j = (j_{1}, … j_{r}) is a multi-index of integer elements with 1 ≤ j_{r} ≤ i − 1 and j! := q_{1}!q_{2}! … q_{i−1}! where q_{k} is the number of occurrences of k in (j_{1}, j_{2} … j_{r}).

A general expression for the k-th derivative of P_{i} is given as follows: since g_{i}(0) = 0, we may let f_{i}(t) = g_{i}(t)/t, i = 1, 2, …, d. We can then rewrite P_{i} in the form

P_{i}(g) = I + ∑_{r=1}^{∞} t^{r} ∑_{j_{1}=1}^{i−1} ∑_{j_{2}=j_{1}}^{i−1} ⋯ ∑_{j_{r}=j_{r−1}}^{i−1} \frac{1}{j!}f_{j_{1}} ⋯ f_{j_{r}}ad_{V_{j_{1}}} ∘ ⋯ ∘ ad_{V_{jr}} (V_{i}).

By following the construction in(Owren & Marthinsen 1997) we obtain

\begin{aligned}\frac{d^{k}P_{i}}{dt^{k}} |_{t=0} & = & ∑_{r=1}^{k} ∑_{δ_{1}+...+δ_{r}=k} \frac{k!}{∏_{ν=1}^{μ}(δ_{ν} − 1)!} ∑_{1≤j_{1}≤...≤j_{μ}≤i−1} \frac{1}{j!} \\ & & × f_{j_{1}}^{(δ_{1}−1)} ⋯ f_{j_{μ}}^{(δ_{μ}−1)}|_{t=0} ad_{V_{j_{1}}} ∘ ⋯ ∘ ad_{V_{jμ}} (V_{i}). \end{aligned}  (2.6)

Substituting (2.6) in (2.4) and (2.5) we obtain the conditions for arbitrary order p. In particular we obtain the following formulae for the derivatives of P_{i}(g) at t = 0,

\frac{dP_{i}}{dt} |_{t=0} = ∑_{k=1}^{i−1} ad_{V_{k}} (V_{i}) g′_{k}(0),

\frac{d^{2}P_{i}}{dt^{2}} |_{t=0} = ∑_{k=1}^{i−1} (∑_{l=1}^{k−1} 2ad_{V_{l}}ad_{V_{k}}(V_{i})g′_{k}(0)g′_{l}(0) + ad_{V_{k}}^{2}(V_{i})[g′_{k}(0)]^{2} + ad_{V_{k}}(V_{i})g′′_{k}(0)) ,

\frac{d^{3}P_{i}}{dt^{3}} |_{t=0} = 6 ∑_{k=3}^{i−1} ∑_{l=2}^{k−1} ∑_{j=1}^{l−1} ad_{V_{j}}ad_{V_{l}}ad_{V_{k}}(V_{i})g′_{j}(0)g′_{l}(0)g′_{k}(0)

+ 3 ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} (ad_{V_{l}}ad_{V_{k}}^{2}(V_{i})g′_{l}(0)[g′_{k}(0)]^{2} + ad_{V_{l}}^{2}ad_{V_{k}}(V_{i})[g′_{l}(0)]^{2}g′_{k}(0))

+ ∑_{k=1}^{i−1} ad_{V_{k}}^{3}(V_{i})[g′_{k}(0)]^{3}

+ 3 ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} ad_{V_{l}}ad_{V_{k}}(V_{i}) (g′′_{k}(0)g′_{l}(0) + g′_{k}(0)g′′_{l}(0))

+ 3 ∑_{k=1}^{i−1} ad_{V_{k}}^{2}(V_{i})g′′_{k}(0)g′_{k}(0)

+ ∑_{k=1}^{i−1} ad_{V_{k}}(V_{i})g′′′_{k}(0).

Substitution readily produces order conditions. Specifically,

∑_{i=1}^{d} g′′_{i}(0)V_{i} = − ∑_{i=1}^{d} g′_{i}(0) ∑_{k=1}^{i−1} ad_{V_{k}}(V_{i})g′_{k}(0),  (2.7)

are conditions for order 2, while

∑_{i=1}^{d} g′′′_{i}(0)V_{i} = − ∑_{i=1}^{d} {2g′′_{i}(0) ∑_{k=1}^{i−1} ad_{V_{k}}(V_{i})g′_{k}(0)

+ g′_{i}(0) ∑_{k=1}^{i−1} [2 ∑_{l=1}^{k−1} ad_{V_{l}}ad_{V_{k}}(V_{i})g′_{k}(0)g′_{l}(0)  (2.8)
+ ad_{V_{k}}^{2}(V_{i})(g′_{k}(0))^{2} + ad_{V_{k}}(V_{i})g′′_{k}(0)]} ,

are the order-3 conditions. Finally, conditions for order 4 are

∑_{i=1}^{d} g_{i}^{IV}(0)V_{i} = − ∑_{i=1}^{d} {3g′′′_{i}(0) ∑_{k=1}^{i−1} ad_{V_{k}}(V_{i})g′_{k}(0)+

+ 3g′′_{i}(0) [∑_{k=1}^{i−1} (∑_{l=1}^{k−1} 2ad_{V_{l}}ad_{V_{k}}(V_{i})g′_{k}(0)g′_{l}(0)

+ ad_{V_{k}}^{2}(V_{i})(g′_{k}(0))^{2} + ad_{V_{k}}(V_{i})g′′_{k}(0))]

+ g′_{i}(0) [6 ∑_{k=3}^{i−1} ∑_{l=2}^{k−1} ∑_{j=1}^{l−1} ad_{V_{j}}ad_{V_{l}}ad_{V_{k}}(V_{i})g′_{j}(0)g′_{l}(0)g′_{k}(0)

+ 3 ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} (ad_{V_{l}}ad_{V_{k}}^{2}(V_{i})g′_{l}(0)(g′_{k}(0))^{2} + ad_{V_{l}}^{2}ad_{V_{k}}(V_{i})(g′_{l}(0))^{2}g′_{k}(0)) (2.9)

+ ∑_{k=1}^{i−1} ad_{V_{k}}^{3}(V_{i})(g′_{k}(0))^{3}

+ 3 ∑_{k=2}^{i−1} ∑_{l=1}^{k−1} ad_{V_{l}}ad_{V_{k}}(V_{i}) (g′′_{k}(0)g′_{l}(0) + g′_{k}(0)g′′_{l}(0))

+ 3 ∑_{k=1}^{i−1} ad_{V_{k}}^{2}(V_{i})g′′_{k}(0)g′_{k}(0)

+ ∑_{k=1}^{i−1} ad_{V_{k}}(V_{i})g′′′_{k}(0)]} .

In Figure 1 we have plotted along the y axis the 2-norm of the error of the approximation of exp(tB) with the second-kind coordinates (SKC) methods of order ranging from 1 to 4. The values of the error are plotted against time, (along the x-axis), to logarithmic scale for matrices of 𝔰𝔩(5). The methods have been implemented using the standard basis defined in section 3.

The computation of g′′(0), g′′′(0) and g^{IV}(0) is obtained directly implementing the formulas (2.7), (2.8),(2.9) respectively. This implementation does not depend on the choice of the particular basis of 𝔰𝔩(5), but the number of commutators that must be computed with this approach is O(d^{p}) for p = 2, 3, 4. Even if we assume that the V_{i}s are very sparse matrices and that the cost of computing each commutator is O(1) operations, the total cost exceeds O(n^{2p}) flops for p = 2, 3, 4 where n is the dimension of the matrix. Such expense is not acceptable for a competitive method of approximation of exp(tB). Fortunately, it can be decreased a very great deal by an appropriate choice of the basis {V_{1}, V_{2}, …, V_{d}}. This is the theme of the next section.

Figure 1: Error in the approximation of the exponential with WN technique.

3 Choosing a basis

The choice of the right basis and sparse representation of commutators are critical to the implementation of the SKC methods. Recalling the order conditions (2.7)–(2.9), our aim is to choose a basis so that terms of the form ad_{V_{i_{1}}} ad_{V_{i_{2}}} ⋯ ad_{V_{is}} V_{j} can be represented in the most economical manner. We recall that, given the basis {V_{1}, V_{2}, …, V_{d}} of a d-dimensional Lie algebra 𝔤, the structure constants are the numbers c_{k,l}^{i}, k, l, i = 1, 2, …, d, such that

[V_{k}, V_{l}] = ∑_{i=1}^{d} c_{k,l}^{i}V_{i}

(Humphreys 1972). Let

B = ∑_{k=1}^{d} β_{k}V_{k}.

Then an order-1 condition is always

g′_{k}(0) = β_{k},  k = 1, 2, …, d.  (3.1)

To obtain the order-2 condition we substitute (3.1) in (2.7) and express commutators in terms of structure constants,

∑_{k=1}^{d} g′′_{k}(0)V_{k} = − ∑_{l=1}^{d} β_{l} ∑_{j=1}^{l−1}[V_{j}, V_{l}]β_{j} = − ∑_{l=1}^{d} β_{l} ∑_{j=1}^{l−1} ∑_{k=1}^{d} c_{j,l}^{k}V_{k}

= − ∑_{k=1}^{d} (∑_{l=1}^{d} ∑_{j=1}^{l−1} β_{l}c_{j,l}^{k}β_{j}) V_{k}.

Since c_{j,l}^{k} = −c_{l,j}^{k}, we thus deduce that

g′′_{k}(0) = ∑_{l=1}^{d} ∑_{j=1}^{l−1} β_{l}c_{l,j}^{k}β_{j},  k = 1, 2, …, d.  (3.2)

Likewise, substituting in (2.8),

∑_{k=1}^{d} g′′′_{k}(0)V_{k} = − ∑_{i=1}^{d} {2g′′_{i}(0) ∑_{l=1}^{i−1}[V_{l}, V_{i}]β_{l} + β_{i} ∑_{l=1}^{i−1} (2 ∑_{j=1}^{l−1}[V_{j}, [V_{l}, V_{i}]]β_{l}β_{j}
+ [V_{l}, [V_{l}, V_{i}]]β_{l}^{2} + [V_{l}, V_{i}]g′′_{l}(0))} .

Note that

[V_{j}, [V_{l}, V_{i}]] = ∑_{s=1}^{d} c_{l,i}^{s}[V_{j}, V_{s}] = ∑_{k=1}^{d} ∑_{s=1}^{d} c_{l,i}^{s}c_{j,s}^{k}V_{k}.

Therefore

∑_{k=1}^{d} g′′′_{k}(0)V_{k} = −2 ∑_{i=1}^{d} g′′_{i}(0) ∑_{l=1}^{i−1} ∑_{k=1}^{d} c_{l,i}^{k}β_{l}V_{k} − 2 ∑_{i=1}^{d} β_{i} ∑_{l=1}^{i−1} ∑_{j=1}^{l−1} ∑_{k=1}^{d} ∑_{s=1}^{d} c_{l,i}^{s}c_{j,s}^{k}β_{l}β_{j}V_{k}

− ∑_{i=1}^{d} β_{i} ∑_{l=1}^{i−1} ∑_{k=1}^{d} ∑_{s=1}^{d} c_{l,i}^{s}c_{l,s}^{k}β_{l}^{2}V_{k} − ∑_{i=1}^{d} ∑_{l=1}^{i−1} β_{i} ∑_{k=1}^{d} c_{l,i}^{k}g′′_{l}(0)V_{k}

and we deduce that

g′′′_{k}(0) = ∑_{i=1}^{d} ∑_{l=1}^{i−1} c_{i,l}^{k}[2g′′_{i}(0)β_{l} + β_{i}g′′_{l}(0)] + 2 ∑_{i=1}^{d} ∑_{l=1}^{i−1} ∑_{j=1}^{l−1} ∑_{s=1}^{d} c_{i,l}^{s}c_{j,s}^{k}β_{i}β_{l}β_{j}

− ∑_{i=1}^{d} ∑_{i=1}^{d} ∑_{l=1}^{i−1} ∑_{s=1}^{d} c_{i,l}^{s}c_{l,s}^{k}β_{i}β_{l}^{2}, k = 1, 2, …, d.

Bearing in mind that for order p we require

α_{k}(t) = ∑_{r=1}^{p} \frac{1}{r!}g_{k}^{(r)}(0)t^{r}, k = 1, 2, …, d,

we observe that the sheer volume of calculations required for the evaluation of the functions α_{1}, α_{2}, …, α_{d} is prohibitive for, say, order 3, unless most of the structure constants vanish. Fortunately, bases of finite-dimensional Lie algebras which are ‘sparse’ (in the sense that a very high proportion of structure constants vanish) are known. They are associated with root space decompositions of Lie algebras (Humphreys 1972) and, in the case of semisimple algebras, are known as Chevalley bases (Carter, Segal & Macdonald 1995). Wishing to avoid too much Lie-algebraic terminology in a numerical analysis paper, we reserve our exposition to just three examples which are the most important in a range applications.

The orthogonal group Let 𝔤 = 𝔰𝔬(n), the Lie algebra of n × n skew-symmetric matrices. It corresponds to two important Lie groups: the orthogonal group O(n) of n × n orthogonal matrices and its subgroup, the special orthogonal group SO(n) of matrices with unit determinant. Its dimension is d = \frac{1}{2}n(n − 1). We let

F_{i,j} = e_{i}e_{j}^{T} − e_{j}e_{i}^{T}, i = 1, 2, …, n, j = i + 1, i + 2, …, n,

where e_{i} is the i-th canonical vector of ℝ^{n}. In other words, F_{i,j} is a matrix whose (i, j)-th element is 1, the (j, i)-th element equals −1 and zero otherwise. We can trivially expand each B ∈ 𝔰𝔬(n) as B = ∑_{i=1}^{n} ∑_{j=i+1}^{n} b_{i,j}F_{i,j}. U(t) := exp(tF_{i,j}) is simply an Euler rotation in the (i, j) plane: it is identity matrix, except that

[ \begin{matrix}U_{i,i} & U_{i,j} \\ U_{j,i} & U_{j,j} \end{matrix} ] = [ \begin{matrix}cos(tF_{i,j}) & sin(tF_{i,j}) \\ − sin(tF_{i,j}) & cos(tF_{i,j}) \end{matrix} ] .

Noting that

[F_{i,j}, F_{l,k}] = \{ \begin{matrix}−F_{j,k}, & i = l, j ≠ k, \\ −F_{i,l}, & i ≠ l, j = k, \\ F_{i,k}, & i ≠ k, j = l, \\ F_{j,l}, & i = k, j ≠ l, \\ O, & otherwise, \end{matrix}

the order conditions are simplified as follows,

p ≥ 1 :  g′_{i,j}(0) = b_{i,j},
p ≥ 2 :  g′′_{i,j}(0) = ∑_{s=j+1}^{n} b_{j,s}b_{i,s} − ∑_{r=i+1}^{j−1} b_{r,j}b_{i,r}

+ ∑_{r=1}^{i−1} b_{r,j}b_{r,i}, 1 ≤ i < j ≤ n,

and similarly for higher-order terms. Thus, the cost of computing the coefficients for the secondorder method is just \frac{1}{2}(n − 2)(n − 1)n ≈ \frac{1}{2}n^{3} flops. In comparison, a naive computation of (3.2), without exploiting sparsity of structure constants, requires \frac{1}{8}(n^{2} − n − 2)(n − 1)^{2}n^{2} ≈ \frac{1}{8}n^{6} flops.

A more classical composition method for B ∈ 𝔰𝔬(n) is the Strang splitting which we can write in the form

e^{tb_{1,2}F_{1,2}/2} ⋯ e^{tb_{n−2,n}F_{n−2,n}/2}e^{tb_{n−1,n}F n−1,n}e^{tb_{n−2,n}F_{n−2,n}/2} ⋯ e^{tb_{1,2}F_{1,2}/2}

(Celledoni & Iserles 1998). It gives a second-order approximant to exp(tB) whose calculation requires ≈ 4n^{3} flops, in comparison with ≈ 3n^{3} for the second-order CSK method.

We note that, in the specific case of 𝔰𝔬(n), diagonal Padé approximants to the exponential provide an alternative to our method, since they map the algebra to O(n). Having said this, for dense matrices B the cost of evaluating the second-order approximant (I − \frac{1}{2}tB)^{−1}(I + \frac{1}{2}tB) with, say, LU factorization is O(n^{3}), comparative with our method.

The special linear group Let 𝔤 = 𝔰𝔩(n), the set of n × n matrices with zero trace, whence d = n^{2} − 1. We split the algebra in the first instance into diagonal and off-diagonal parts: in the terminology of Lie algebras, the subspace spanned by the diagonal elements is a Cartan subalgebra (Carter et al. 1995) or maximal toral algebra (Humphreys 1972) of 𝔰𝔩(n). Specifically, our basis is

{E_{i,j} : i, j = 1, 2, …, n i ≠ j} ∪ {D_{i} : i = 1, 2, …, n − 1}.

where

E_{i,j} = e_{i}e_{j}^{T},  i, j = 1, 2, …, n, i ≠ j,
D_{i} = e_{i}e_{i}^{T} − e_{i+1}e_{i+1}^{T},  i = 1, 2, …, n − 1.

The exponentials of E_{i,j} and D_{i} are trivial,

e^{tE_{i,j}} = I + tE_{i,j}, e^{tD_{i}} = e^{t}D_{i}.

We order the elements by taking first E_{i,j}, i ≠ j, in lexicographic order, followed by D_{1}, D_{2}, …, D_{n−1}. The commutator table is

[E_{i,j}, E_{r,s}] = \{ \begin{matrix}E_{i,s}, & i ≠ s, j = r, \\ −E_{r,j}, & i = s, j ≠ r, \\ ∑_{l=s}^{r−1} D_{l}, & i = s < j = r, \\ − ∑_{l=r}^{s−1} D_{l}, & i = s > j = r, \\ O, & otherwise, \end{matrix}  \begin{matrix}i, j, r, s = 1, 2, …, n, \\ i ≠ j, r ≠ s, \end{matrix}

[E_{i,j}, D_{r}] = \{ \begin{matrix}−E_{r,j}, & i = r, j ≠ r + 1, \\ −E_{i,r+1}, & i ≠ r, j = r + 1, \\ −2E_{r,r+1}, & i = r, j = r + 1, \\ E_{r+1,j}, & i = r + 1, j ≠ r, \\ E_{i,r}, & i ≠ r + 1, j = r, \\ 2E_{r+1,r}, & i = r + 1, j = r, \\ O, & otherwise, \end{matrix}  \begin{matrix}i, j = 1, …, n, i ≠ j, \\ r = 1, …, n − 1, \end{matrix}

[D_{i}, D_{j}] = O,  i, j = 1, 2, …, n − 1.

In general, for a d-dimensional Lie algebra there are (d − 1)d^{2} structure constants. In the case of 𝔰𝔩(n) this means that up to ≈ n^{6} structure constants may be nonzero. Yet, using the above basis results in just 2(n − 1)n^{2} + 4(n − 2)(n − 1) + \frac{2}{3}(n^{2} − 1)n ≈ \frac{7}{3}n^{3} nonzero structure constants and substantial saving in the implementation of the SKC technique.

Letting

[E_{i,j}, E_{r,s}] = ∑_{(k,l)} c_{(i,j),(k,l)}^{(k,l)}E_{k,l} + ∑_{k} c_{(i,j),(r,s)}^{(k)}D_{k},
[E_{i,j}, D_{r}] = ∑_{(k,l)} c_{(i,j),r}E_{k,l}

(note that [D_{r}, E_{i,j}] = −[E_{i,j}, D_{r}] and [D_{i}, D_{j}] = O) we thus have

c_{(i,j),(r,s)}^{(k,l)} = \{ \begin{matrix}+1, & k = i, l = s, r = j, s ≠ i, \\ −1, & k = r, l = j, r ≠ j, s = i, \\ 0, & otherwise, \end{matrix}

c_{(i,j),(r,s)}^{k} = \{ \begin{matrix}+1, & i = s < j = r, k = s, s + 1, …, r − 1, \\ −1, & i = s > j = r, k = r, r + 1, …, s − 1, \\ 0, & otherwise, \end{matrix}

c_{(i,j),r}^{(k,l)} = \{ \begin{matrix}+1, & k = i = r + 1, l = j, j ≠ r or k = i, l = j = r, i ≠ r + 1, \\ −1, & k = i = r, l = j, j ≠ r + 1 or k = i, l = j = r + 1, i ≠ r, \\ +2, & k = i = r + 1, l = j = r, \\ −2, & k = i = r, l = j = r + 1, \\ 0, & otherwise, \end{matrix}

c_{(i,j),r}^{k} = c_{r,s}^{(k,l)} = c_{r,s}^{k} = 0.

Letting

B = ∑_{k≠l} β_{k,l}E_{k,l} + ∑_{k} γ_{k}D_{k},

and ordering the pairs (k, l), k ≠ l, in lexicographic order, we thus have

g′_{k,l}(0) = β_{k,l},
g′_{k}(0) = γ_{k},
g′′_{k,l}(0) = ∑_{(i,j)≻(r,s)} β_{i,j}c_{(i,j),(r,s)}^{(k,l)}β_{r,s} + ∑_{(i,j),r} β_{i,j}c_{(i,j),r}^{k}γ_{r}

= ∑_{i=1}^{k−1} β_{k,i}β_{i,l} − ∑_{i=k+1}^{n} β_{k,i}β_{i,l} + β_{k,l}(γ_{k−1} + γ_{l} − γ_{k} − γ_{l−1}),

g′′_{k}(0) = ∑_{(i,j)≻(r,s)} β_{i,j}c_{(i,j),(r,s)}^{k}β_{r,s} = − ∑_{i=1}^{k} ∑_{j=k+1}^{n} β_{i,j}β_{j,i},

where γ_{0} = γ_{n} = 0.

The Lorenz group This is the 6-dimensional group SO(3, 1) of 4×4 matrices A such that AJA^{T} = J, where J = diag(1, 1, 1, −1) (Carter et al. 1995). It has important applications in special relativity theory. he corresponding Lorenz algebra 𝔰𝔬(3, 1) consists of all matrices B such that BJ +JB^{T} = O. It is easy to verify that each element of 𝔰𝔬(3, 1) can be written in the form

B = [ \begin{matrix}0 & b_{1} & b_{2} & b_{3} \\ −b_{1} & 0 & b_{4} & b_{5} \\ −b_{2} & −b_{4} & 0 & b_{6} \\ b_{3} & b_{5} & b_{6} & 0 \end{matrix} ] , b_{1}, b_{2}, …, b_{6} ∈ ℝ.

Choosing the basis

\{[ \begin{matrix}0 & 1 & 0 & 0 \\ −1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} ] , [ \begin{matrix}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ −1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} ] , [ \begin{matrix}0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & −1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} ] , [ \begin{matrix}0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{matrix} ] , [ \begin{matrix}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{matrix} ] , [ \begin{matrix}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{matrix} ]\} ,

we obtain the commutator table

\begin{matrix}[V_{1}, V_{2}] = −V_{3}, & [V_{1}, V_{3}] = V_{2}, & [V_{1}, V_{4}] = −V_{5}, & [V_{1}, V_{5}] = V_{4}, & [V_{1}, V_{6}] = O, \\ [V_{2}, V_{1}] = V_{3}, & [V_{2}, V_{3}] = −V_{1}, & [V_{2}, V_{4}] = −V_{6}, & [V_{2}, V_{5}] = O, & [V_{2}, V_{6}] = V_{4}, \\ [V_{3}, V_{1}] = −V_{2}, & [V_{3}, V_{2}] = V_{1}, & [V_{3}, V_{4}] = O, & [V_{3}, V_{5}] = −V_{6}, & [V_{3}, V_{6}] = V_{5}, \\ [V_{4}, V_{1}] = V_{5}, & [V_{4}, V_{2}] = V_{6}, & [V_{4}, V_{3}] = O, & [V_{4}, V_{5}] = V_{1}, & [V_{4}, V_{6}] = V_{2}, \\ [V_{5}, V_{1}] = −V_{4}, & [V_{5}, V_{2}] = O, & [V_{5}, V_{3}] = V_{6}, & [V_{5}, V_{4}] = −V_{1}, & [V_{5}, V_{6}] = V_{3}, \\ [V_{6}, V_{1}] = O, & [V_{6}, V_{2}] = −V_{4}, & [V_{6}, V_{3}] = −V_{5}, & [V_{6}, V_{4}] = −V_{2}, & [V_{6}, V_{5}] = −V_{3}. \end{matrix}

Thus, out of 180 structure constants, just 24 are nonzero – and they all equal ±1. After brief claculation, we drive for example the polynomials α_{k} that yield an order-2 CSK approximant,

α_{1}(t) = β_{1}t + \frac{1}{2}(β_{2}β_{3} − β_{4}β_{5})t^{2},
α_{2}(t) = β_{2}t − \frac{1}{2}(β_{1}β_{3} + β_{4}β_{6})t^{2},
α_{3}(t) = β_{3}t + \frac{1}{2}(β_{1}β_{2} − β_{5}β_{6})t^{2},
α_{4}(t) = β_{4}t − \frac{1}{2}(β_{1}β_{5} + β_{2}β_{6})t^{2},
α_{5}(t) = β_{5}t + \frac{1}{2}(β_{1}β_{4} − β_{3}β_{6})t^{2},
α_{6}(t) = β_{6}t + \frac{1}{2}(β_{3}β_{5} − β_{2}β_{4})t^{2},

where B = ∑_{k=1}^{6} β_{k}V_{k}.

4 Time symmetry

An approximant F (tB) ≈ exp(tB) is said to be time symmetric if F (tB)F (−tB) = I, t ≥ 0. Time symmetric approximants are important for a number of reasons, not least being that they lend themselves to the Yošida technique, which allows their order to be increased (Yoshida 1990). The techniques of the last section are not time symmetric. Here we describe their modification, which results in a time-symmetric approximant.

Me mention in passing that it is possible to envisage two distinct techniques to obtain high-order algorithms based on canonical coordinates of the second kind. The first, implicit in the work of the previous section, consists of evaluating the numbers g_{k}^{(l)}(0) for l = 1, 2, …, p, where p is the order of the method. The alternative, the subject matter of the present section, consists in combining a second-order or a fourth-order approximant across a number of steps to obtain a higher-order method.

Given the splitting

B = ∑_{l=1}^{s} C_{l},

it is well known that the Strang splitting The approximation

F (tB) = e^{tC_{1}/2} ⋯ e^{tC_{s−1}/2}e^{tC_{s}}e^{tC_{s−1}/2} ⋯ e^{tC_{1}/2}  (4.1)

is of order 2 and time symmetric. Note that, as a consequence of time symmetry, for sufficiently small t ≥ 0 we can represent F (tB) = e^{F(t)} where the matrix function F(t) is odd. It is precisely this feature that allows the application of the Yošida technique.

The clear reason for (4.1) being time symmetric is that it is palindromic in the alphabet {C_{1}, C_{2}, …, C_{s}}. This provides a clue how to modify techniques based on canonical coordinates of the second kind so as to render them time symmetric. Given a basis {V_{1}, V_{2}, …, V_{d}} of the Lie algebra 𝔤, we approximate e^{tB} by the product

exp[α_{1}(t)V_{1}] ⋯ exp[α_{d−1}(t)V_{d−1}] exp[α_{d}(t)V_{d}] exp[α_{d−1}(t)V_{d−1}] ⋯ exp[α_{1}(t)V_{1}],  (4.2)

where α_{1}, α_{2}, …, α_{d} are odd polynomials.

Taking α_{l} = \frac{1}{2}β_{l}t, l = 1, 2, …, d − 1 and α_{d} = β_{d}t yeilds the second-order Strang splitting. In the sequel we seek higher-order methods of this kind.

Using the Baker–Campbell–Hausdorff (BCH) formula it is possible to express the product of exponentials at the right hand side of (4.2) as a single exponential (Varadarajan 1984, p. 141). Due to the symmetric arrangement of the exponentials in (4.2), the BCH formula is an expansion in odd powers of t. If this expansion converges, which is always the case for sufficently small t, it makes sense to write the equation

tB = 2 ∑_{i=1}^{d−1} α_{i}(t)V_{i} + α_{d}(t)V_{d} + ∑_{k=1}^{∞} Q^{2k}(α).  (4.3)

Here we denote by Q^{2k}(α) the terms of order O(t^{2k+1}) in the BCH formula applied to (4.3). Moreover, we let α_{i}^{2k}(t) be the polynomial obtained by truncating the expansion of α_{i}(t) after the first k terms, and we denote the remainder by r_{i}^{2k}(t). In other words,

α_{i}(t) = α_{i}^{2k}(t) + r_{i}^{2k}(t),  r_{i}^{2k}(t) = O(t^{2k+1}) ,  i = 1, 2, …, d.

From (4.3) we deduce

2 ∑_{i=1}^{d−1} α_{i}^{2(k)}(t)V_{i} + α_{d}^{2(k)}(t)V_{d} = tB − ∑_{r=1}^{k−1} Q^{2r}(α) + O(t^{2k+1}) .

Noting that

Q^{2r}(α) = Q^{2r}(α^{2(k−1)} + r^{2(k−1)}) = Q^{2r}(α^{2(k−1)}) + O(t^{2k+r}) ,

we obtain

2 ∑_{i=1}^{d−1} α_{i}^{2(k)}(t)V_{i} + α_{d}^{2(k)}(t)V_{d} = tB − ∑_{r=1}^{k−1} Q^{2r}(α^{2(k−1)}) + O(t^{2k+1}) .  (4.4)

Dropping the O(t^{2k+1}) terms in (4.4), it is possible to compute α^{2k} from α^{2(k−1)}. This gives a procedure to derive a sequence of successively increaing-order approximants of exp(tB). It is easy to see that the approximants

F^{2k}(tB) = exp(α_{1}^{2k}(t)V_{1}) ⋯ exp(α_{d}^{2k}(t)V_{d}) ⋯ exp(α_{1}^{2k}(t)V_{1}),

of exp(tB) are such that F^{2k}(tB)F^{2k}(−tB) = I, hence time symmetry, the reason being the symmetric arrangements of the exponentials in F^{2k}(tB) and the odd-power expansion of the functions α_{i}^{2k}.

The BCH and symmetric BCH formulae for k-terms have an exceedingly complicated expansion, which can be obtained recursively. In what follows we will make use just of the term Q^{2}(α), demonstrating how it is possible to compute it explicitely for particular choices of the basis.

In the remainder of this section we consider the implementation of time-symmetric CSK methods. We split B as before and commence by considering the Strang splitting (4.1) except that, to simplify notation, we arrange the terms in reverse ordering,

e^{tC_{s}/2} ⋯ e^{tC_{2}/2}e^{tC_{1}}e^{tC_{2}/2} ⋯ e^{tC_{s}/2}.  (4.5)

Lemma 1 The term Q^{2} of the BCH formula applied to (4.5) is

Q^{2} = \frac{t^{3}}{12} ∑_{l=2}^{s}[C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, [C_{1} + ⋯ + C_{l−1}, C_{l}]].  (4.6)

Proof See the appendix. □

Let us next consider the case 𝔤 = 𝔰𝔬(n), choosing the same sparse basis as in Section 2. Therefore, according to (4.6), we have

\begin{aligned}Q^{2} & = & \frac{t^{3}}{12} ∑_{i=1}^{n−1} ∑_{j=i+1}^{n}[b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, [b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, \\ & & b_{i,j}F_{i,j}]] + \frac{1}{24} ∑_{i=1}^{n−1} ∑_{j=i+1}^{n} b_{i,j}[F_{i,j}, [b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, b_{i,j}F_{i,j}]]. \end{aligned}  (4.7)

We compute separately each part of this sum. Exploiting the commutator table of our basis we have

[b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, b_{i,j}F_{i,j}]
= b_{i,j} (− ∑_{s=i+1}^{j−1} b_{i,s}F_{s,j} − ∑_{r=1}^{i−1} b_{r,j}F_{r,i} + ∑_{t=1}^{i−1} b_{t,i}F_{t,j}) ,

and noting that b_{i,s} = −b_{s,i}, we deduce that

[b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, b_{i,j}F_{i,j}] = b_{i,j} (∑_{{t=1 \atop t≠i}}^{j−1} b_{t,i}F_{t,j} − ∑_{r=1}^{i−1} b_{r,j}F_{r,i}) .  (4.8)

Commuting the right-hand side with F_{i,j} gives

[F_{i,j}, ∑_{{t=1 \atop t≠i}}^{j−1} b_{t,i}F_{t,j} − ∑_{r=1}^{i−1} b_{r,j}F_{r,i}] = ∑_{{t=1 \atop t≠i}}^{j−1} b_{t,i}F_{t,i} + ∑_{r=1}^{i−1} b_{r,j}F_{r,j}.  (4.9)

Let b_{1}, b_{2}, …, b_{n} be the columns of B and denote

b_{l}^{s} = ∑_{k=1}^{s−1} b_{k,l}e_{k}, k = 1, 2, …, n.

Then (4.8) yields

[b_{1,2}F_{1,2} + … + b_{i,j−1}F_{i,j−1}, b_{i,j}F_{i,j}] = b_{i,j}(b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T}) − b_{i,j}(b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T}),

while (4.9) gives

[F_{i,j}, [b_{1,2}F_{1,2} + ⋯ + b_{i,j−1}F_{i,j−1}, b_{i,j}F_{i,j}]]
= b_{i,j}(b_{i}^{j}e_{i}^{T} − e_{i}b_{i}^{j}^{T}) − b_{i,j}(b_{j}^{i}e_{j}^{T} − e_{j}b_{j}^{i}^{T}).

Multiplying the latter by b_{i,j} and summing in i and j we can evaluate (4.7) in n^{3} operations. Note that we count separately multiplications and additions, for example, we assume that the cost of Euclidean inner product of two vectors of length n is 2n operations.

We now assemble together our results to calculate (4.7). We proceed by splitting the sum b_{1,2}F_{1,2} + ⋯ + b_{i,j−1}F_{i,j−1} in three parts, whereby

[b_{1,2}F_{1,2} + ⋯ + b_{i,j−1}F_{i,j−1}, b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T} − (b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T})]

= ∑_{l=1}^{i} ∑_{k=1}^{l−1}[(b_{l}^{l}e_{l}^{T} − e_{l}b_{l}^{l}^{T}), b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T} − (b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T})]

+ ∑_{l=i+1}^{j−1} ∑_{k=1}^{i}[(b_{l}^{i+1}e_{l}^{T} − e_{l}b_{l}^{i+1}^{T}), b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T} − (b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T})]

+ ∑_{l=j}^{m} ∑_{k=1}^{i−1}[(b_{l}^{i}e_{l}^{T} − e_{l}b_{l}^{i}^{T}), b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T} − (b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T})].

Finally,

[b_{1,2}F_{1,2} + ⋯ + b_{i,j−1}F_{i,j−1}, b_{i}^{j}e_{j}^{T} − e_{j}b_{i}^{j}^{T} − (b_{j}^{i}e_{i}^{T} − e_{i}b_{j}^{i}^{T})]

= −b_{i}^{i}^{T}b_{i}^{j}F_{i,j} + b_{i}^{i}b_{j}^{i}^{T} − b_{j}^{i}b_{i}^{i}^{T} + b_{j}^{i}^{T}b_{j}^{i}F_{j,i} − (b_{j}^{i}b_{i}^{j}^{T} − b_{i}^{j}b_{j}^{i}^{T})

∑_{l=1}^{i−1} b_{l,i}(b_{l}^{l}e_{j}^{T} − e_{j}b_{l}^{l}^{T}) − b_{l}^{l}^{T}b_{i}^{j}F_{l,j} − b_{l,j}(b_{l}^{l}e_{i}^{T} − e_{i}b_{l}^{l}^{T}) + b_{l}^{l}^{T}b_{j}^{i}F_{l,i}

∑_{l=i+1}^{j−1} b_{l,i}(b_{l}^{i+1}e_{j}^{T} − e_{j}b_{l}^{i+1}^{T}) − b_{l}^{i+1}^{T}b_{i}^{j}F_{l,j} + b_{i,l}(b_{j}^{i}e_{l}^{T} − e_{l}b_{j}^{i}^{T}) + b_{l}^{i+1}^{T}b_{j}^{i}F_{l,i}

+ ∑_{l=j+1}^{n} b_{l}^{i}^{T}b_{j}^{i}F_{l,i} − b_{l}^{i}^{T}b_{i}^{j}F_{l,j}.

We analyse the computational costs of the previous formula, summing over i and j and showing that (4.7) can be computed in O(n^{3}) operations. Note that, since ∑_{l=i+1}^{j−1} b_{i,l}e_{l} = b_{i}^{i} − b_{i}^{j}, we have

b_{i}^{i}b_{j}^{i}^{T} − b_{j}^{i}b_{i}^{i}^{T} + ∑_{l=i+1}^{j−1} b_{i,l}(b_{j}^{i}e_{l}^{T} − e_{l}b_{j}^{i}^{T}) − (b_{j}^{i}b_{i}^{j}^{T} − b_{i}^{j}b_{j}^{i}^{T}) = −2(b_{j}^{i}b_{i}^{j}^{T} − b_{i}^{j}b_{j}^{i}^{T}).

It is more convenient to write the previous expression in the form

−2 ∑_{i=1}^{n−1} ∑_{j=i+1}^{n} b_{i,j}(b_{j}^{i}b_{i}^{j}^{T} − b_{i}^{j}b_{j}^{i}^{T})

= − ∑_{i=1}^{n} ( ∑_{j=i+1}^{n} 2b_{i,j}b_{j}^{i}) b_{i}^{i}^{T} − b_{i}^{i} ( ∑_{j=i+1}^{n} 2b_{i,j}b_{j}^{i})^{T}

− ∑_{i=1}^{n−1} ∑_{j=i+1}^{n} 2b_{i,j} (b_{j}^{i}(b_{i}^{j} − b_{i}^{i})^{T} − (b_{i}^{j} − b_{i}^{i})b_{j}^{i}^{T}) .

The first part of this sum is computed in about \frac{2}{3}n^{3} operations and the second part, exploiting the equality

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} 2b_{i,j}b_{j}^{i}(b_{i}^{j} − b_{i}^{i})^{T} = 2 ∑_{i=1}^{n−1} ∑_{k=n−1}^{i+2} b_{i,k} (∑_{l=n}^{k+1} b_{i,l}b_{l}^{i}) e_{k}^{T},

can also be computed in \frac{2}{3}n^{3} operations.

Adding terms of the type αF_{l,j} and βF_{l,i} leads to

− ∑_{l=1}^{i} b_{l}^{l}^{T}b_{i}^{j}F_{l,j} − ∑_{l=i+1}^{j−1} b_{l}^{i+1}^{T}b_{i}^{j}F_{l,j} − ∑_{l=j+1}^{n} b_{l}^{i}^{T}b_{i}^{j}F_{l,j} = −L_{i}^{i}b_{j}^{i}e_{j}^{T} − e_{j}(−L_{i}^{i}b_{j}^{j})^{T},

and

∑_{l=1}^{i−1} b_{l}^{l}^{T}b_{j}^{i}F_{l,i} + ∑_{l=i+1}^{j−1} b_{l}^{i+1}^{T}b_{j}^{i}F_{l,i} + ∑_{l=j+1}^{n} b_{l}^{i}^{T}b_{j}^{i}F_{l,i} = L_{i}^{j}b_{i}^{j}e_{i}^{T} − e_{i}(L_{i}^{j}b_{i}^{j})^{T},

where the matrix L_{i} is the lower triangular part of b_{1,2}F_{1,2} + … + b_{i−1,n}F_{i−1,n} and we denote by L_{i}^{s}, s = i, j, the matrix L_{i} with zeros along its s-th row. Summing up with respect to i and j, we obtain

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} b_{i,j}(−L_{i}^{i}b_{j}^{i})e_{i}^{T} = ∑_{i=1}^{n−1}(−L_{i}^{i}) ( ∑_{j=i+1}^{n} b_{i,j}b_{j}^{i}) e_{i}^{T}

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} b_{i,j}L_{i}^{j}b_{i}^{j}e_{j}^{T} = ∑_{j=2}^{n} (∑_{i=1}^{j−1} b_{i,j}c_{i}) e_{j}^{T}

where we have used the notation c_{i} := L_{i}b_{i}^{i} for i = 1, …, n − 1. The cost of computing the first sum is \frac{4}{3}n^{3}, while the cost of computing the second is 2n^{3} operations.

Finally the terms

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} ∑_{l=1}^{i} b_{i,j}b_{l,i}(b_{l}^{l}e_{j}^{T} − e_{j}b_{l}^{l}^{T}) = ∑_{l=1}^{n−1} b_{l}^{l} ( ∑_{j=l+2}^{n} c_{j,l}e_{j}^{T}) − ( ∑_{j=l+2}^{n} c_{j,l}e_{j}) b_{l}^{l}^{T}

with c_{j,l} = ∑_{i=l}^{j−1} b_{i,j}b_{l,i},

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} ∑_{l=1}^{i} b_{i,j}b_{l,j}(b_{l}e_{i}^{T} − e_{i}b_{l}^{T}) = ∑_{l=1}^{n−1} b_{l}^{l} ( ∑_{i=l+1}^{n−1} d_{i,l}e_{i}^{T}) − ( ∑_{i=l+1}^{n−1} d_{i,l}e_{i}) b_{l}^{l}^{T}

with d_{i,l} = ∑_{j=i+1}^{n} b_{i,j}b_{l,j}, and

∑_{i=1}^{n−1} ∑_{j=i+1}^{n} ∑_{l=i+1}^{j−1} b_{i,j}b_{l,i}(b_{l}^{i+1}e_{j}^{T} − e_{j}b_{l}^{i+1}^{T})

= ∑_{i=1}^{n−1} ∑_{l=i+1}^{n−1} b_{l,i} (b_{l}^{i+1}(b_{i}^{l+1} − b_{i})^{T} − (b_{i}^{l+1} − b_{i})b_{l}^{i+1}^{T});

can be computed in about \frac{2}{3}n^{3}, \frac{2}{3}n^{3} and \frac{1}{2}n^{3} operations respectively. Collecting the contributions of all the terms in the sum we obtain a total count of 7\frac{1}{2}n^{3} operations.

At the present time it is not clear that this method of computation of Q^{2} in the 𝔰𝔬(n) case is optimal form the point of view of complexity theory. We did not try any other ordering of the basis elements and it is not at all certain that different orderings could give better constants in front of the term n^{3}.

Given that the construction of the (second-order) Strang splitting carries a cost of 4n^{3} operations, the total flop count for constructing a symmetric fourth order SKC approximation of an exponential in 𝔰𝔬(n) by our algorithm is 11\frac{1}{2}n^{3}. This is marginally better than obtaining an order-4 approximation by the Yos̆ida technique from three Strang splittings which, as pointed out in (Celledoni & Iserles 1998), requires 12n^{3} flops.

5 Sparse matrices

In a naive formulation, the method of canonical coordinates of the second kind is considerably too expensive for practical computation. This, however, can be alleviated by the use of a sufficiently ‘sparse’ basis of the underlying Lie algebra 𝔤. As explained in Section 2, choosing a basis so that an overwhelming majority of structure constants vanish renders the algorithm strikingly more effective. It is important to emphasize that this has nothing to do with the structure of the matrix B ∈ 𝔤, which need not be sparse. Yet, in most practical computations (in particular when n is large) one can expect B to be sparse and structured. Good algorithms should be able to exploit this phenomenon.

In the case of SKC methods we identify two mechanisms that allow us to exploit sparsity. Although this aspect of our methods is still a matter for active investigation, the interim results are substantive enough to warrant publication. For simplicity, we describe the first mechanism just in the case of a tridiagonal B ∈ 𝔰𝔬(n), hence

B = [ \begin{matrix}0 & β_{1} & 0 & ⋯ & 0 \\ −β_{1} & 0 & ... & ... & ⋮ \\ 0 & ... & ... & ... & 0 \\ ⋮ & ... & ... & 0 & β_{n−1} \\ 0 & ⋯ & 0 & −β_{n−1} & 0 \end{matrix} ] = ∑_{k=1}^{n−1} β_{k}F_{k,k+1},

where the matrices F_{k,l} = e_{k}e_{l}^{T} − e_{l}e_{k}^{T} have been introduced in Section 2. Since

b_{k,l} = \{ \begin{matrix}β_{k}, & l = k + 1, \\ 0, & otherwise, \end{matrix}

it is easy to substitute in the general formulae for the order-2 method:

g′_{i,j}(0) = \{ \begin{matrix}β_{i}, & j = i + 1, \\ 0, & otherwise, \end{matrix} g′′_{i,j}(0) = \{ \begin{matrix}β_{i−1}β_{i}, & j = i − 2, \\ 0, & otherwise, \end{matrix} 1 ≤ i < j ≤ n,

Arranging the elements of the basis in lexicographic order, we thus obtain the second-order approximant

e^{β_{n−1}tF_{n−1,n}}e^{β_{n−2}tF_{n−2,n−1}}e^{\frac{1}{2} β_{n−1}β_{n−1}t^{2}F_{n−2,n}} ⋯ e^{β_{2}tF_{2,3}}e^{\frac{1}{2} β_{2}β_{3}t^{2}F_{2,4}}e^{β_{1}tF_{1,2}}e^{\frac{1}{2} β_{1}β_{2}t^{2}F_{1,3}}.

In other words, the cost of the approximation is just O(n) flops.

Similar situation pertains to

B = [ \begin{matrix}γ_{1} & η_{1} & 0 & ⋯ & 0 \\ μ_{1} & γ_{2} & ... & & ⋮ \\ 0 & ... & ... & ... & 0 \\ ⋮ & & ... & γ_{n−1} & η_{n−1} \\ 0 & ⋯ & 0 & μ_{n−1} & γ_{n} \end{matrix} ] ∈ 𝔰𝔩(n).

Choosing the same basis and terminology as in Section 2 we can readily ascertain that

g′′_{k,k−2}(0) = μ_{k−2}μ_{k−1}  k = 3, 4, …, n
g′′_{k,k−1}(0) = −(γ_{k−2} − 2γ_{k−1} + γ_{k})μ_{k−1},  k = 3, 4, …, n,
g′′_{k,k+1}(0) = (γ_{k−1} − 2γ_{k} + γ_{k+1})η_{k},  k = 2, 3, …, n − 1,
g′′_{k,k+2}(0) = −η_{k}η_{k+1},  k = 1, 2, …, n − 2,
g′′_{k,l}(0) = 0,  |k − l| ≥ 3

and g′′_{k}(0) = −η_{k}μ_{k}, k = 1, 2, …, n−1. Thus, a second-order approximant to a tridiagonal B ∈ 𝔰𝔩(n) is itself quindiagonal and its computation requires just O(n) flops.

Higher-order approximants and matrices with greater bandwidth lend themselves to similar treatment, although the savings are less striking. In a sense, the situation is parallel to that of approximating exp(tB) by a rational approximant, when savings accrue from sparse matrix-inversion methods, except that in our case the result is assured to belong to the right Lie group.

Another observation which is highly pertinent to the approximation of exponentials of sparse matrices has been made in (Iserles 1999). Suppose that B is a banded matrix of bandwidth s ≥ 3. In general, F (t) = exp(tB) is a dense matrix. Yet, as is easy to illustrate by computer experiments, F (t) is very near to a banded matrix. Specifically, given ε > 0, there exists r = r(t, ε) ≥ s such that all the elements of F (t) outside a band of width r are less than ε in magnitude. Moreover, tight upper bounds on r can be derived with relative ease. The idea thus is to set to zero all the elements outside bandwidth r. The outcome is a banded approximant to the exponential. Moreover, with an appropriate choice of basis elements, this means that the functions α_{i} are set to zero for elements that possess terms exclusively outside the band. Consequently corresponding exponentials equal identity and need not be included in the product. Thus, the cost scales with the size r of the bandwidth. Similar phenomenon has been already encountered in the context of 𝔰𝔬(n) and 𝔰𝔩(n), when our choice of basis and order has implied a banded structure of the exponential. The present mechanism is different, even if the net outcome is similar.

6 Numerical experiments

Our numerical experiments are organized as follows. We fist consider a test on random matrices in 𝔰𝔬(50), illustrating the performance of methods based on the use of second kind coordinates techniques for full and sparse matrices. The third and last example is the solution of a third-order ODE using Runge–Kutta/Munthe-Kaas (RK/MK) methods described in (Munthe-Kaas 1997). We use the Matlab toolbox DiffMan for the integration of ODEs on manifolds, comparing the usual implementation of RK/MK methods, whereby the the exponential is approximated to machine accuracy, with a version of the methods obtained using the time-symmetric fourth-order approximation from Section 4.

All experiments have been performed in Matlab and we have computed the error while comparing the results with the built-in function expm which calculates the exponential to nearly machine accuracy.

Figure 2: Error versus time in the 𝔰𝔬(50) (full case).

We evaluated the the error computing ∥e^{−tB}F (tB)−I∥_{F} where F (tB) is the SKC approximation of exp(tB) and ∥ · ∥_{F} denotes the Frobenius norm. The matrices have been generated randomly using the Matlab function rand and scaling the Frobenius norm so that ∥B∥_{F} = 1.

We approximate exp(tB) with a single step of the methods for different values of t, (t = 1/2^{k} and k = 1, …, 5).

In both the first two figures the norm of the error is plotted (along the y-axis) to a logarithmic scale with respect to t. Figure 2 reports the results of our first test, where we have considered a full matrix in 𝔰𝔬(50). In the plots the error norm is indicated with the symbols ‘∗’ (SKC, time symmetric, order 4) and ‘∘’ (Strang splitting, order 2).

In the next example, illustrated in Figure 3, the same methods have been applied to a sparse matrix in 𝔰𝔬(50), with four non-zero diagonals (i.e., bandwidth 5). In both the examples the methods give the correct order. In the second case, however, the count of flops is drastically reduced. We counted the number of flops using the Matlab function flops. In the first case the cost for constructing Q^{2} amounts to 9.62n^{3} while in the second we counted 0.95n^{3} flops. As it is easy to understand, the described implementation of the methods allows to take advantage immediately of the sparsity structure of the matrix B, working directly on the nonzeros entries of B à la Section 4.

The last example is concerned with the use of the techniques described in this paper in substituting the exponentials computed to machine accuracy in the integration methods of (Munthe-Kaas 1997). The experiments have been performed using the Matlab toolbox DiffMan. We use a RK/MK method of order four.

The example is a problem whose solution is the soliton originating in the Korteweg–de Vries (KdV) equation. It is a third-order ODE obtained performing a symmetry reduction on the KdV equation. The resulting ODE can be written as a three-dimensional system,

y′ = [ \begin{matrix}0 & 1 & 0 \\ 0 & 0 & 1 \\ −9y(2) & 3 & 0 \end{matrix} ] y

with y(0) = [1, 0, −1.5]^{T} and t ∈ [0, 5]. The solution of the ODE f = y_{1}(t) can be easily derived explicitely and it is f(t) = αsech (tβ), α = 1, β = 1/2\sqrt{(}3).

In Figure 4 we plot the analytic solution (solid line) on a grid of 161 points. The dotted line is the numerical solution obtained with the Matlab routine ode45 with absolute and relative tolerance 1.0e − 4. The method produced this solution in 69 steps, and it was implemented with step-size

Figure 3: Error versus time in the 𝔰𝔬(50) sparse case.

control procedure. The dashed-dotted line is the numerical solution obtained with the RK/MK method using SKC symmetric tecnhniques for the approximation of the exponential, with fixed step-size h.

In Figure 5 we plot the error (along the y-axis) with respect to the numerical solution obtained with the Matlab routine ode45 to a logarithmic scale, versus the stepsize h = 1/2^{k} and k = 1, …, 5 for the cases of the implementation of RK/MK with the expm function of Matlab (marked with +) and approximating exp to order four with a SKC technique (∘). The line marked with ∗ representes the error of the numerical solution given by the RK/MK method implemented with SKC technique for the approximation of the exponential, measured with respect to the numerical solution obtained by the same method with the use of the exact exponential (expm routine of Matlab).

It is interesting to note in this case that substituting the exact exponential with suitable fourthorder approximant does not lead to a significant deterioration in the quality of the RK/MK method and the overall error does not change much. Note that in the present case the primary variable is a vector, rather than a matrix. In general, if the underlying ODE can be written in a vector form, i.e. as an action of a Lie group on ℝ^{n}, we need to approximate exp(tB)v, where v ∈ ℝ^{n}, rather than the matrix exp(tB). This leads to obvious savings in the SKC techniques, similarly, say, to the approach of rational functions. In particular, the cost of composing exponentials is O(n^{2}), rather than O(n^{3}), operations.

Acknowledgments

The authors are grateful to Brynjulf Owren for many fruitful discussions, to Per Christian Moan for bringing the reference (Wei & Norman 1963) to their attention, to the Numerical Analysis group of DAMTP Cambridge, and to the Geometric Integration members during the fall semester 1998 at MSRI Berkeley. Research at MSRI is supported in part by NSF grant DMS-9701755.

References

Carter, R., Segal, G. & Macdonald, I. (1995), Lectures on Lie Groups and Lie Algebras, LMS Student Texts, Cambridge University Press, Cambridge.

Figure 4: The soliton originating in the KdV equation. Figure 5: RK/MK: global error at t = 5 with expm and SKC Casas, F. (1996), ‘Fer’s factorization as a symplectic integrator’, Numer. Math. 74(3), 283–303.

Celledoni, E. & Iserles, A. (1998), Approximating the matrix exponential from a Lie algebra to a Lie group, Technical Report 1998/NA03, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, England.

Chacon, R. & Fomenko, A. (1991), ‘Recursion formulas for the Lie integral’, Advances in Math. 88, 200–257.

Crouch, P. E. & Grossman, R. (1993), ‘Numerical integration of ordinary differential equations on manifolds’, J. Nonlinear Sci. 3, 1–33.

Engø, K. (1998), On the construction of geometric integrators in the RKMK class, Technical Report 158, Department of Informatics, University of Bergen, Norway.

Hochbruck, M., Lubich, C. & Selhofer, H. (1998), ‘Exponential integrators for large systems of differential equations’, SIAM J. Sci. Comput. 19, 1552–1574.

Humphreys, J. E. (1972), Introduction to Lie algebras and Representation Theory, First edn, Springer.

Iserles, A. (1999), How large is the exponential of a banded matrix?, Technical Report DAMTP 1999/1, University of Cambridge.

Iserles, A. & Nørsett, S. P. (1997), Linear ODEs in Lie groups, Technical Report 1997/NA9, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, England.

Iserles, A. & Nørsett, S. P. (1999), ‘On the solution of linear differential equations in Lie groups’, Phil. Trans. Royal Soc. A. To appear.

Munthe-Kaas, H. (1997), High order Runge-Kutta methods on manifolds, Technical Report NA14, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, England.

Owren, B. & Marthinsen, A. (1997), Runge-Kutta methods adapted to manifolds and based on rigid frames, Technical Report Numerics No. 1/1997, Department of Mathematical Sciences, The Norwegian University of Science and Technology. To appear in BIT.

Owren, B. & Marthinsen, A. (1998), Integration methods based on canonical coordinates of the second kind, In preparation.

Varadarajan, V. S. (1984), Lie Groups, Lie Algebras, and their Representation, GTM 102, Springer-Verlag.

Wei, J. & Norman, E. (1963), ‘On global representations of the solutions of linear differential equations as a product of exponentials’, Advances in Mathematics.

Yoshida, H. (1990), ‘Construction of higher order symplectic integrators’, Physics Letters A 150, 262–268.

Zanna, A. (1997), Collocation and relaxed collocation for the Fer and the Magnus expansions, Technical Report 1997/NA17, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, England.

A Appendix

For completeness, we present a proof of Lemma 1. Note that a comprehensive treatment of this subject matter, inclusive of the non-symmetric case, has been presented in a different context by Chacon & Fomenko (1991). For our purposes, however, it is sufficient to derive the first term of the expansion.

Lemma 1 The leading error term in the Strang splitting is

\frac{1}{12} ∑_{l=2}^{s}[C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, [C_{1} + ⋯ + C_{l−1}, C_{l}]].  (4.6)

Proof Letting

F_{1}(t) = e^{tC_{1}},
F_{l}(t) = e^{tC_{l}/2}F_{l−1}e^{tC_{l}/2},  l = 2, 3, …, s,

we can verify at once that F_{s} is precisely the Strang splitting. We assume that

F_{l}(t) = exp[t(C_{1} + ⋯ + C_{l}) + \frac{1}{12}Q_{l}t^{3} + O(t^{4})], l = 1, 2, …, s.

We use the BCH formula:

F_{l}(t) = exp(\frac{1}{2}tC_{l}) exp[t(C_{1} + ⋯ + C_{l−1}) + \frac{1}{12}Q_{l−1}t^{3} + O(t^{4})] exp(\frac{1}{2}tC_{l})
= exp\{t(C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}) + \frac{1}{4}t^{2}[C_{l}, C_{1} + ⋯ + C_{l−1}]
+ \frac{1}{24}t^{3}[\frac{1}{2}C_{l} − (C_{1} + ⋯ + C_{l−1}), [C_{l}, C_{1} + ⋯ + C_{l−1}]]
+ \frac{1}{12}t^{3}Q_{l−1} + O(t^{4})\} exp(\frac{1}{2}tC_{l})
= exp\{t(C_{1} + ⋯ + C_{l}) + \frac{1}{4}t^{2}[C_{l}, C_{1} + ⋯ + C_{l−1}]
+ \frac{1}{24}t^{3}[\frac{1}{2}C_{l} − (C_{1} + ⋯ + C_{l−1}), [C_{l}, C_{1} + ⋯ + C_{l−1}]]
+ \frac{1}{4}t^{2}[C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, C_{l}] + \frac{1}{16}t^{3}[[C_{l}, C_{1} + ⋯ + C_{l−1}], C_{l}]
+ \frac{1}{24}t^{3}[C_{1} + ⋯ + C_{l−1}, [C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, C_{l}]] + \frac{1}{12}t^{3}Q_{l−1} + O(t^{4})\}.

However,

[C_{l}, C_{1} + ⋯ + C_{l−1}] + [C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, C_{l}] = O,

thereby annihilating the t^{2} term, and

\frac{1}{24}[\frac{1}{2}C_{l} − (C_{1} + ⋯ + C_{l−1}), [C_{l}, C_{1} + ⋯ + C_{l−1}]] + \frac{1}{16}[[C_{l}, C_{1} + ⋯ + C_{l−1}], C_{l}]
+ \frac{1}{24}[C_{1} + ⋯ + C_{l−1}, [C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, C_{l}]]
= \frac{1}{12}[C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, [C_{1} + ⋯ + C_{l−1}, C_{l}]].

Therefore

Q_{l} = [C_{1} + ⋯ + C_{l−1} + \frac{1}{2}C_{l}, [C_{1} + ⋯ + C_{l−1}, C_{l}]] + Q_{l−1}.

Since Q_{1} = O, the expression (4.6) follows by summing the above formula for l = 2, 3, …, s. □