Generalized Polar Decompositions for the Approximation of the Matrix Exponential.
In this paper we describe the use of the theory of generalized polar decompositions [H. Munthe-Kaas, G. R. W. Quispel, and A. Zanna,  Found. Comput. Math., 1 (2001), pp. 297--324] to approximate a matrix exponential. The algorithms presented have the property that, if $Z \in {\frak{g}}$, a Lie algebra of matrices, then the approximation for exp(Z) resides in G, the matrix Lie group of ${\frak{g}}$. This property is very relevant when solving Lie-group ODEs and is not usually fulfilled by standard approximations to the matrix exponential.    We propose algorithms based on a splitting of Z into matrices having a very simple structure, usually one row and one column (or a few rows and a few columns), whose exponential is computed very cheaply to machine accuracy.    The proposed methods have a complexity of ${\cal O}(\kappa n^{3})$, with constant $\kappa$  small, depending on the order and the Lie algebra ${\frak{g}}$. % The algorithms are recommended in cases where it is of fundamental importance that the approximation for the exponential resides in G, and when the order of approximation needed is not too high. We present in detail algorithms up to fourth order.
Introduction
With the recent developements in the theory of Lie-group integration schemes for ordinary differential
equations (Iserles, Munthe-Kaas, N-rsett & Zanna 2000), the problem of approximating
the matrix exponential has lately received a renewed attention. Most Lie-group methods require
a number of computation of matrix exponentials from a Lie algebra g ' R n\Thetan to a Lie group
G ' GL(n; R), that usually constitutes a bottleneck in the numerical implementation of the
schemes (Celledoni, Iserles, N-rsett & Orel 1999).
The matrix exponentials need be approximated to the order of the underlying ODE method (hence
exact computation is not an issue), however, it is of fundamental importance that such approximations
resides in G. In generality, this property is not fullfilled by many standard approximations
to the exponential function (Moler & van Loan 1978) unless the exponential is evaluated exactly.
In some few cases (usually for small dimension) the exponential of a matrix can be evaluated exactly.
This happens, for instance, for three by three skew-symmetric matrices, whose exponential can be
Institutt for informatikk, University of Bergen, H-yteknologisenteret, Thorm-hlensgate 55, N-5020 Bergen,
Norway. Email: anto@ii.uib.no, hans@ii.uib.no
calculated exactly by means of the well known Euler-Rodriguez formula
sin ff
ff
where
(Marsden & Ratiu 1994). Exact formulas for skew-symmetric matrices and matrices in so(p; q) can
be derived up to dimension eight making use of the Cayley-Hamilton theorem (Horn & Johnson
1985) with significant savings with respect to approximation techniques (Barut, Zeni & Laufer
1994, Leite & Crouch 1999). However, the algorithms are not practical for larger dimensions, for
several reasons. First, they require high powers of the matrix in question (and each matrix-matrix
multiplication amounts to O
secondly, it is well known that the direct use of
the characteristic polynomial, for large scale matrices, may lead to computational instabilities.
The problem of approximating the exponential of a matrix from a Lie algebra to its corresponding
Lie group has been recently considered by (Celledoni & Iserles 2000, Celledoni & Iserles 1999). In
the first paper, the authors construct the approximation by first splitting the matrix X 2 g as the
sum of bordered matrices. Strang-type splittings of order two are considered, so that one could
apply a Yoshida technique (Yoshida 1990), based on a symmetric composition of a basic scheme
whose error locally expands in odd powers of time only, to increase the order. In the second
paper, the authors consider techniques based on canonical coordinates of the second kind (CCSK)
(Varadarajan 1984). To follow that approach, it is necessary to choose a basis of the Lie algebra g.
The choice of the basis plays a significant role in the computational complexity of the algorithms
(Owren & Marthinsen 1999), and, by choosing Chevalley bases (Carter, Segal & Macdonald 1995)
which entail a large number of zero structure constants, it is possible to reduce significantly the
cost of the methods from O
\Delta to O
In this paper we consider the problem of approximating to a given order of accuracy
F (t; Z) - exp(tZ) 2 G; Z 2
so that F (t; Z) 2 G, where g ' gl(R; n) and G ' GL(R;n). The techniques we introduce consist
in a Lie-algebra splitting of the matrix Z by means of an iterated generalized polar decomposition
induced by an appropriate involutive automorphism oe G, as discussed in (Munthe-Kaas
et al. 2000b). We introduce a general technique for approximations of arbitrary high order, and
discuss practical algorithms of order two, three and four. For large n, these algorithms are very
competitive with standard approximations of the exponential function (for example diagonal Pad'e
approximants).
The paper is organized as follows. In Section 2 we discuss the background theory of the polar
decomposition on Lie groups and its symmetric version. Such polar decomposition can be used to
induce splitting in the Lie algebra g. As long as this splitting is practical to compute, together
with the exponential of each 'splitted' part, it leads to splitting methods for the approximation of
the exponential of practical interest.
In Section 3 we use the theory developed in x2 to derive approximations of the exponential function
for some relevant matrix Lie groups as SO(R;n), and SL(R; n). Methods of order two, three and
four are discussed in greater detail, together with their computational complexity. The methods
are based on splittings in bordered matrices, whose exact exponentials are very easy to compute.
Section 4 is devoted to some numerical experiments where we illustrate the results derived in this
paper, and finally Section 5 is devoted to some concluding remarks.
Background theory
It is usual in differential geometry to denote Lie-group elements with lower case letters and Lie-
algebra elements with upper-case letters, whether they represent matrices, vectors or scalars (Hel-
gason 1978). We adopt this convention throughout this section.
Let G be a Lie group with Lie algebra g. We restrict our attention to matrix groups, i.e to the
case when G ' GL(R;n).
It is known that, provided G is an involutive automorphism of G, every element z 2 G
sufficiently close to the identity can be decomposed in the product
wg, the subgroup of elements of G fixed under oe and
is the subset of anti-fixed points of oe (Lawson 1994, Munthe-Kaas et
al. 2000b). The set G oe has the structure of a symmetric space (Helgason 1978) and is closed under
the product
as it can be easily verified by application of oe to the right-hand-side of the above relation. The
decomposition (2:1) is called the polar decomposition of z in analogy with the case of real matrices
with the choice of automorphism
g. The automorphism oe induces an involutive automorphims doe
on g in a natural manner,
d
dt
and it defines a splitting of the algebra g into the sum of two linear spaces,
Zg is a subalgebra of g, while \GammaZ g has the
structure of a Lie-triple system, a set closed under the double commutator,
To keep our presentation relevant to the argument matter of this paper, we refer the reader to
(Munthe-Kaas et al. 2000b, Munthe-Kaas, Quispel & Zanna 2000a) and references therein for a
more extensive treatement of such decompositions. However, it is of fundamental importance to
note that the sets k and p possess the following properties:
We denote by \Pi the canonical projection onto the subspace p and by \Pi k
k the
projection onto k. Then,
where
Assume that x and y in (2:1) are of the form
and they can be expanded in series
where the X i and Y i can be explicitely calculated by means of the following recurrence relations
c 2'
and
k1 ;:::;k 2k?0
k1+\Delta\Delta\Delta+k 2k =2q
(2m)! ad 2m
Z
(Zanna 2000). Note that Y (t) expands in odd powers of t only. The first terms in the expansions
of X(t) and Y (t) are
O
We also consider a symmetric-type generalized polar decomposition,
where, as above, X(t) 2 p and Y (t) 2 k. To compute X(t), we apply oe to both sides of (2:7) to
obtain
Isolating the y term in (2:8) and (2:7) and equating the result, we obtain
This leads to a differential equation for X which is very similar to the one obeyed by Y in (2:5)
(Zanna 2000). Using the recursions in (Zanna 2000) we obtain recursions for X(t) and Y (t). The
first terms are given as
and both X(t) and Y (t) expand in odd powers of t only.
Generalized polar decomposition and its symmetric version
for the approximation of the exponential
Assume now that we wish to approximate exp(tZ) for some Z 2 g, and that oe 1 is an involutive
automorphism so that the exponential of terms in as well as analytic
functions of adP , are easy to compute. Then and we can approximate
where X [1] and Y [1] obey the order conditions (2:4)-(2:6) to suitable order.
Alternatively, we can approximate
where X [1] and Y [1] obey now the order conditions (2:10)-(2:11) to given accuracy.
The same mechanism can be applied to split k 1 in p 2 \Phi k 2 by means of a suitable automorphism
oe 2 . The procedure can be iterated and, provided that the exponential of k m is easy to compute,
we have an algorithm to approximate exp(tZ) to a given order of accuracy. In this circumstance,
(3:1) will read
while the analogue of (3:2) is
both corresponding to the algebra splitting
3.1 On the choice of the automorphisms oe i
In what follows, we will consider automorphisms oe of the form
G; (3.6)
where S is an idempotent matrix, i.e. S I [Munthe-Kaas and Zanna, 2000]. Clearly,
and for simplicity, we will abuse notation writing oeZ in place of doeZ, given that all our computations
take place in the space of matrices.
all the eigenvalues of S are either +1 or \Gamma1. Thus, powers of matrices
as well as powers of adP , are easy to evaluate by means of the (+1)- and (\Gamma1)-eigenspace of S
(Munthe-Kaas & Zanna 2000).
3.2 Automorphisms that lead to banded matrices splittings
Let Z 2 gl(n; R) be a n \Theta n matrix and consider the automorphism
is the idempotent matrix
. 0
It is easy to verify that
Z =2
z
z
while
\Pi k1
z
In general, assume that, at the j-th step, the space consists of matrices of the form
. O
O w
O
Then, the obvious choice is
O ~
. 0
where I j \Gamma1 denotes the (j \Gamma 1) \Theta (j \Gamma 1) identity matrix and ~
so that the subspace p j consists of matrices of the form
O ~
Exponentials of matrices of the form (3:11) are very easy to compute: in effect,
exp
O ~
O exp( ~
where exp( ~
can be computed exactly either with a formula analogous to the Euler-Rodriguez
formula (1:1): denote a
exp( ~
I
~
~
a T
I
I
~
~
\Gammaa T
Note that
~
Another alternative for the exact exponential of ~
is the one proposed in (Celledoni & Iserles
exp( ~
where
a j
e 1 is the vector [1; finally 1)=z. The latter formula (3:13),
as we shall see in the sequel, leads to significant savings in the computation and assembly of the
exponentials.
Moreover, given that
O ~
where
~
w j;j a
Next, if Z 2 g, to obtain an approximation of the exponential in G by these automorphisms, we
shall require that oe i 's, defined by the above matrices S i , map g into g. Clearly, this is the case for
ffl so(n; R), since oe i
Z is a map from so(n) ! so(n) given that each S i is an
orthogonal matrix;
ffl sl(n; R), since oe i leaves the diagonal elements of Z (hence its trace) unchanged;
ffl quadratic Lie algebras and the
commute. This is for instance the case when J is diagonal, hence our formulas are valid
for not for the symplectic algebra sp(n; R). In the latter situation,
we consider different choices for the automorphisms oe i , discussed at a greater length in
(Munthe-Kaas & Zanna 2000).
3.3 Splittings of order two to four, their implementation and complexity
In this section we describe in more details the algorithms, the implementation and the complexity
of the splittings induced by the automorphisms described above. The cases of a polar-type repre-
sentation, xy, or a symmetric polar-type representation, z = xyx, are discussed separately.
Algorithm 1 (Polar-type splitting, order two) Based on the iterated generalized polar decomposition
(3:3).
Note that the \Pi p j
and \Pi k j
projections need not be stored in separate matrices but can be stored
in places of the rows and columns of the matrix Z. We truncate the expansions (2:6) to order two,
hence at each step only the p j -part needs correction. Taking in mind (3:3), the matrices X [j] are
low rank matrices with nonzero entries only on the j-th row, column
row are stored in place of the corresponding Z entries. The
matrix Y [n\Gamma1] is diagonal and is stored in the diagonal entries of Z.
Purpose: 2nd order approximation of the splitting (3:3)
overwritten with the nonzero elements of X [i] and Y [m] as:
a
The computation of the splitting requires at each step two matrix-vector multiplications, each
amounting to O
floating point operations (we count both multiplications and ad-
ditions), as well as two vector updates, which are O(n operations. Hence, for large n, the
cost of computing the splitting is of the order
ffl 2n 3 for so(n), taking into account that b
Note that both for so(p; q) and so(n) the matrix Y [n\Gamma1] is the zero matrix.
Algorithm 2 (Symmetric polar-type splitting, order two) Based on the iterated generalized
polar decomposition (3:4)
We truncate the expansions (2:10)-(2:11) to order two. The storing of the entries is as above.
Purpose: 2nd order approximation of the splitting (3:4)
overwritten with the nonzero elements of X [i] and Y [m] as:
% Computation of the splitting
a
This splitting costs only
ffl n(n\Gamma1)for so(n), because of skew-symmetry.
Algorithm 3 (Polar-type splitting, order three)
We truncate (2:6)-(2:7) to include O
terms. Note that the term [K; [P; K]] is of the form (3:15).
We need include also the term of the form [P; [P; K]]. We observe that
Purpose: 3rd order approximation of the splitting (3:3)
overwritten with the nonzero elements of X [i] and Y [m] as:
% Computation of the splitting
a
Analyzing the computations involved, the most costly part is constituted by the matrix-vector
products in the computations in c products
in the update of Z(j j). The computation of c
amounting to 8n 3 in the whole process. For the update of Z(j
need to compute two vector-vector products (O
operations
to uptdate the elements of the matrix. Thus, the whole cost of updating the matrix Z(j
n) is 5n 3 . The update of z j;j requires operations per step, which give a 2n 3
contribution to the total cost of the splitting.
In summary, the total cost of the splitting is
ffl 5n 3 for so(p; q) and sl(n)
ffl for so(n), note that d j need not be calculated as well as z Similarly, we take into
account that b and that only half of the elements of Z(j need be
updated. The total amounts to 2 1
It is easy to modify the splitting above to obtain order four. Note that
which requires the computation of the scalar b T
costing 2=3n 3 operations in the whole pro-
cess. However, all the other powers ad i
~
~
no further computation. Next
can be computed with just two (one) extra matrix-vector computations
for sl(n) (resp. so(n)), which contribute 4n 3 (resp. 2n 3 ) to the cost of the splitting, so that the
splitting of order four costs a total of 7n 3 operations for sl(n) (resp. 4n 3 for so(n)).
Algorithm 4 (Symmetric polar-type splitting, order four)
We truncate (2:10)-(2:11) to include O
terms. Also in this case, the term [K; [P; K]] is of the
form (3:15), while the term [P; [P; K]] is computed according to (3:16).
Purpose: 4th order approximation of the splitting (3:4)
overwritten with the nonzero elements of X [i] and Y [m] as:
a
We need to compute a total of four matrix-vector products, yielding 8
operations. The update
of the block Z(j costs 5n 3 operations, while the update of z(j; costs 2n 3
operations, for a total of
ffl 5n 3 operations for sl(n) and so(p
operations for so(n).
3.4 On higher order splittings
The costs of implementing splittings following (3:3) or (3:4) depend on the type of commutation
involved: commutators of the form [P; K] and [P contribute as an
O
\Delta term to the total complexity of the splitting, however, commutators of the form [K
for easily contribute an O
\Delta to the total complexity of the splittings if the
special structure of the terms involved is not taken into consideration. If carefully implemented,
also these terms can be computed with only matrix-vector and vector-vector products, contributing
O
\Delta operations to the total cost of the splitting. For example, let us consider the term
which appears in the O
contribution in the expansion of the Y part,
both for the polar-type and symmetric polar-type splitting. One has
denotes the matrix z j;j I \Gamma -
. The parenthesis indicate the correct order in which the
operations should be executed to obtain the right complexity (O
per iteration, hence
a total of O
\Delta for the splitting). Many of the terms are already computed for the lower order
conditions, yet the complexity arises significantly. Therefore we recommend this splitting type
techniques when a moderate order of approximation is required.
To construct higher order approximations with these splitting techniques, one could use our symmetric
polar-type splittings, together with a Yoshida-type symmetric combination.
3.5 Assembly of the approximation F (t; Z) to the exponential
For each algorithm that computes the approximation to the exponential, we distinguish two cases:
when the approximation is applied to a vector v and when instead the matrix exponential exp(Z)
is required. Since the matrices X [j] are never constructed explicitely and are stored as vectors,
computations of the exponentials exp(X [j] ) is also never performed explicitely but it is implemented
as in the case of the Householder reflections (Golub & van Loan 1989) when applied to a vector.
First, let us return to (3:13). It is easy to verify that, if we denote by ff
has the exact form
I
where I is the 2 \Theta 2 identity matrix. Similar remarks hold about the matrix D \Gamma1 . Thus, the
computation of can be done in a very few flops that `do not contribute'
to the total cost of the algorithm. Next, if v; k; the assembly of exp( ~
according to
can be computed in 6j operations. If we let j vary between 1 and n, the total cost of the
multiplications is hence 3n 2 . This is precisely the complexity of for the assembly of the exponential
for polar-type splittings, that has the form as in (3:3).
Algorithm 5 (Polar-type approximation)
Purpose: Computing the approximant (3:3) applied to a vector v
containing the nonzero elements of X [i] and Y [m] as:
a
old
and
new := [a
new .
In the case when the output needs be applied to a n \Theta n matrix B, we can apply the above
algorithm to each column of B, for a total of 3n 3 operations. This complexity can be reduced to
about 2n 3 taking into account that the vector can be calculated once and for
all, depending only on the splitting of the matrix Z and not in any manner on the columns of B.
Also can be computed once and stored for latter use.
Algorithm 6 (Symmetric polar-type approximation)
The approximation to the exponential is carried out in a manner very similar to that described
above in Algorithm 5, except that, being (3:4) based on a Strang-type splitting, the assembly is
also performed in reverse order.
Purpose: Computing the approximant (3:4) applied to a vector v
containing the nonzero elements of X [i] and Y [m] as:
a
old
old
new := [a
new .
a
old

Table

1: Complexity for a polar-type order-two approximant.
Algorithm sl(n); so(p; q) so(n)
1+5 vector matrix vector matrix
splitting 1 1n 3 1 1n 3 2n 3 2n 3
assembly exp 3n 2 2n 3 3n 2 2n 3
new := [a
new
The vectors ff; fi and fl need be calculated only once and stored for latter use in the reverse-order
multiplication. The cost of the assembly is roughly twice as the cost of the assembly in Algorithm 1,
hence it amounts to 5n 2 operations (we save n 2 operations omitting the computation of ff).
When the result is applied to a matrix B, again we apply the same algorithm to each column of
B, which yields n 3 operations. Also in this case the vector ff does not depend on B and can be
computed once and for all, reducing the cost to 4n 3 operations. The same remark holds for the
vectors fi and fl .
It is important to mention that the matrix D might be singular or close to singular (for example
when a j and b j are close to be orthogonal), hence the computation of exp( ~
according to (3:13)
may be lead to instabilities. In this case, it is recommended to use (3:12) instead of (3:13). The
latter choice is twice as expensive (5n 2 for polar-type assemblies and 9n 2 for symmetric assemblies
for F (t; Z) applied to a vector), but deals better with the case when D is nearly singular.
4 Numerical experiments
4.1 Non-symmetric polar-type approximations to the exponential
We commence comparing the polar-type order-2 splitting of Algorithm 1 combined with the assembly
of the exponential in Algorithm 5 with the (1; 1)-Pad'e approximant for matrices in sl(n) and
so(n), with corresponding groups SL(n) and SO(n). We choose diagonal Pad'e approximants as
benchmarck because they are easy to implement, are the rational approximant with highest order
of approximation at the origin and it is well known that they map quadratic Lie algebras into
quadratic Lie groups (but not necessarily other Lie algebras into the corresponding Lie groups).

Table

1 reports the complexity of the method 1+5. A (1; 1)-Pad'e approximant costs O
floating point operations when applied to a vector (essentially the cost of LU-factorising a linear
system) and O
operations when applied to n \Theta n matrices (2n 3 operations come from the
construction of the right-hand-side, 2
3 from the LU factorization and 2n 3 from the n forward
and backward solution of triangular systems).
In

Figure

4.1 we compare the number of floating point operations scaled by n 3 for matrices Z up
to size 500 as obtained in Matlab for our polar-type order-two algorithm (method 1+5) and the
both applied to a matrix. We consider the cases when Z is in sl(n) and
so(n). The costs of computing both approximations clearly converges to the theroretical estimates
(which in the plot are represented by solid lines) given in Table 1 for large n.
Flops/nmethod 1+5, sl(n)
method 1+5, so(n)

Figure

1: Floating point operations (scaled by n 3 ) versus size for the approximation of the exponential
of a matrix in sl(n) and in so(n) applied to a matrix with the order-2 polar-type algorithm
(method 1+5) and (1; 1)-Pad'e approximant.

Table

2: Complexity for a polar-type order-three approximant. The numbers in parenthesis correspond
to the coefficients for an order-four approximation.
Algorithm sl(n); so(p; q) so(n)
3+5 vector matrix vector matrix
splitting 5(7)n 3 5(7)n 3 2 1(4)n 3 2 1(4)n 3
assembly exp 3n 2 2n 3 3n 2 2n 3
total 5(7)n 3 7(9)n 3 2 1(4)n 3 4 1(6)n 3
In

Figure

2 we compare the accuracy of the two approximations (left plot) for the exponential of a
normalized so that kZk methods show a local truncation error of O
revealing that the order of approximation
to the exact exponential is two. The right plot shows the error in the determinant as a function
of the Pad'e approximant has an error that behaves like h 3 , while our method preserves the
determinant equal to one to machine accuracy.
In table 2 we report the complexity of the method 3+5, which yields an approximation to the
exponential of order three. The numbers in parenthesis refer to the cost of the algorithm with
order four corrections.
4.2 Symmetric polar-type approximations to the exponential
We commence comparing our method 2+6, yielding an approximation of order two, with the (1; 1)
Pad'e approximant. Table 3 reports the complexity of the method 2+6.
Clearly, in the matrix-vector case, our methods are one order of magnitude cheaper than the Pad'e
approximant, and are definitively to be preferred (see Figure 3, for matrices in sl(n)). Furthermore,
from
exact
exponential
method 1+5, sl(n)
|det
1|
method 1+5, sl(n)

Figure

2: Error in the approximation (left) and in the determinant (right) versus h for the approximation
of the exponential of a traceless matrix of unit norm with the order-2 polar-type
algorithm (method 1+5) and (1; 1)-Pad'e approximant.

Table

3: Complexity for a symmetric polar-type order-two approximant.
Algorithm sl(n); so(p; q) so(n)
2+6 vector matrix vector matrix
splitting
assembly exp 5n 2 4n 3 5n 2 4n 3
total
our method maps the approximation in SL(n), while the Pad'e approximant does not. When
comparing approximations of the matrix exponential applied to a vector, it is a must to consider
Krylov subspace methods (Saad 1992). We compare the method 2+6 with a Krylov subspace
method when Z is a matrix in sl(n), normalized so that kZk a vector of unit
norm. The Krylov subspaces are obtained by Arnoldi iterations, whose computational cost amounts
to circa 2mn 2 counting both multiplications and additions. Here m is the
dimension of the subspace Km j spanfv; vg. To obtain the total cost of a Krylov
method, we have to add O
computations arising from the evaluation of the exponential of
the Hessenberg matrix obtained with the Arnoldi iteration, plus 2nm operations arising from the
multiplication of the latter with the orthogonal basis. However, when n is large and m - n, these
costs are subsumed in that of the Arnoldi iteration, and the leading factor is 2mn 2 . The error,
computed as and the floating point operations of both approximations for
are given in Table 4. The Krylov method converges very fast: in all the three cases
eight-nine iterations are sufficient to obtain almost machine accuracy, while two iterations yield an
error which is of the order of method 2+6, at about two thirds (0:64; 0:68; 0:69 respectively) the cost.
On the other hand, Krylov methods do not produce an SL(n) approximation to the exponential,
unless the computation is performed to machine accuracy, which, in our particular example, is
3:30, 2:84 and 2:85, circa three times more costly than the 2+6 algorithm. For what the SO(n)
case is concerned, it should be noted that, if Z 2 so(n), then the approximation w - exp(Z)v
produced by the Krylov method has the feature that kwk independently of the number m
of iterations: in this case, the Hessenberg matrix produced by the Arnoldi iterations is tridiagonal
and skew-symmetric, hence its exponential orthogonal. Thus, Krylov methods are the method of
choice for actions of SO(n) on R n (Munthe-Kaas & Zanna 1997). One might extrapolate that,
Floating
point
operations
method 2+6, sl(n)

Figure

3: Floating point operations versus size for the approximation of the exponential of a matrix
in sl(n) applied to a vector with the order-2 symmetric polar-type algorithm (method 2+6) and

Table

4: Krylov subspace approximations versus the method 2+6 for the approximation of exp(Z)v.
Krylov 2+6
size n error m flops error flops0.74 1 21041
7.

Table

5: Complexity for a symmetric polar-type order-four approximant.
Algorithm sl(n); so(p; q) so(n)
4+6 vector matrix vector matrix
splitting 5n 3 5n 3 2 1n 3 2 1n 3
assembly exp 5n 2 4n 3 5n 2 4n 3
total 5n 3 9n 3 2 1
if we wish to compute the exponential exp(Z)Q, where Q 2 SO(n), one could perform only a
few iterations of the Krylov method to compute w being
columns of Q. Unfortunately, the approximation [w ceases to be orthogonal: although
the vectors w i cease to be linearly independent and the final approximation is not in
SO(n). Similar analysis yields for Stiefel manifolds, unless Krylov methods are implemented to
approximate the exponential to machine accuracy.
In passing, we recall that our methods based on a symmetric polar-type decomposition are time-
symmetric. Hence it is possible to compose a basic scheme in a symmetric manner, following
a technique introduced by Yoshida (Yoshida 1990), to obtain higher order approximations: two
orders of accuracy can be obtained at three times the cost of the basic method. For instance we
can use the method 2+6 as a basic algorithm to obtain an approximation of order four. Thus an
approximation of order four applied to a vector can be obtained in 17n 2 operations for sl(n) (two
splittings and three assemblies), compared to O
operations required by the method 4+6.
To conclude our gallery, we compare the method 4+6, an order-four scheme, whose complexity is
described in Table 5, with a (2; 2)-Pad'e approximant, which requires 2 2
floating point operations
when applied to vectors (2n 3 for the assembly and 2
for the LU factorization) and 6 2
matrices (since we have to resolve for multiple right-hand sides). The figures obtained by numerical
simulations for matrices in sl(n) and SO(n) clearly agree with the theoretical asymptotic values
(plotted as solid lines), as shown in Figure 4. The cost of both methods is very similar, as is
the error from the exact exponential although, in the SL(n) case, the 4+6 scheme preserves the
determinat to machine accuracy while the Pad'e scheme does not (see Figure 5).
Conclusions
In this paper we have introduced numerical algorithms for approximating the matrix exponential.
The methods discussed possess the feature that, if Z 2 g, then the output is in G, the Lie group
of g, a property that is fundamental in the integration of ODEs by means of Lie-group methods.
The proposed methods have a complexity of O
denotes the size of the matrix whose
exponential we wish to approximate. Typically, for moderate order (up to order four), the constant
- is less than 10 whereas the exact computation of a matrix exponential in Matlab (which employes
the scaling and squaring with a Pad'e approximant) generally costs between 20n 3 and 30n 3 .
Comparing methods of the same order of accuracy applied to a vector v 2 R n and to a matrix
G:
ffl For the case F (t; Z)v - exp(tZ)v, where v is a vector: Symmetric polar-type methods are
slightly cheaper than their non-symmetric variant. For the SO(n) case, the complexity of
symmetric methods is very comparable with that of diagonal Pad'e approximants of the same
Flops/nmethod 4+6, sl(n)
method 4+6, so(n)

Figure

4: Floating point operations (scaled by n 3 ) versus size for the approximation of the exponential
of a matrix in sl(n) applied to a n \Theta n matrix with the order-4 symmetric polar-type
algorithm (method 4+6) and (2; 2)-Pad'e approximant.
from
exact
exponential
methods 4+6, sl(n)
|det
1|
methods 4+6, sl(n)

Figure

5: Error in the approximation(left) and in the determinant (right) versus h for the approximation
of the exponential of a traceless matrix of unit norm with the order-4 symmetric
polar-type algorithm (method 4+6) and (2; 2)-Pad'e approximant.
order.
The complexity of the method 2+6 is O
while for the rest of our methods it is O
Krylov subspace methods do, however, have the complexity O
if the number of iterations
is independent of n. Thus, if it is important to stay on the group, we recommend Krylov
methods with iteration to machine accuracy for this kind of problems. If convergence of
Krylov methods is slow, our methods might be good alternatives. See (Hochbruck & Lubich
1997) for accurate bounds on the number m of iterations of Krylov methods.
ffl For the case F (t; Z)B - exp(tZ)B, with B an n \Theta n matrix: Non-symmetric polar-type
methods are marginally cheaper than their symmetric counterpart; however the latter should
be preferred when the underlying ODE scheme is time-symmetric. The proposed methods
have a complexity very comparable with that of diagonal Pad'e approximants of the same
order (they require slightly less operations in the SO(n) case) in addition they map sl(n) to
SL(n), a property that is not shared by Pad'e approximants. For these problems our proposed
methods seem to be the best choice.
It should also be noted that significant advantages arise when Z is a banded matrix. For instance,
the cost of method 2+6 scales as O(nr) for F (t; Z) applied to a vector and O
applied to a matrix when Z has bandwidth 2r + 1. The savings are less striking for higher order
methods since commutation usually causes fill-in in the splitting.
Our schemes have an implementation cost smaller than those proposed by (Celledoni & Iserles
1999), that also produce an output in G when Z 2 g. For the SO(n) case, Celledoni et al. propose
an order-four scheme whose complexity is 11 1
our order-four schemes (method 3+5 with
order-four corrections and method 4+6) costs 6n 3 , 6 1
operations - very comparable with the
diagonal Pad'e approximant of the same order. Furthermore, the implementation of the schemes
of Celledoni et al. requires a precise choice of a basis in g, hence the knowledge of the structure
constants of the algebra. Our approach is instead based on the inclusion relations (2:3) and is
easily expressed in very familiar linear algebra formalism.



--R


Lectures on Lie Groups and Lie Algebras
Methods for the approximation of the matrix exponential in a Lie-algebraic setting

Complexity theory for Lie-group solvers
Matrix Computations
Differential Geometry

Matrix Analysis

'Polar and Ol'shanskii decompositions'
'Closed forms for the exponential mapping on matrix Lie groups based on Putzer's method'
Introduction to Mechanics and Symmetry

Numerical integration of differential equations on homogeneous manifolds



Integration methods based on canonical coordinates of the second kind

Lie Groups

Recurrence relation for the factors in the polar decomposition on Lie groups
--TR

--CTR
Ken'Ichi Kawanishi, On the Counting Process for a Class of Markovian Arrival Processes with an Application to a Queueing System, Queueing Systems: Theory and Applications, v.49 n.2, p.93-122, February   2005
Jean-Pierre Dedieu , Dmitry Nowicki, Symplectic methods for the approximation of the exponential map and the Newton iteration on Riemannian submanifolds, Journal of Complexity, v.21 n.4, p.487-501, August 2005
lie algebra;matrix exponential;lie-group integrator