Stability of the Generalized Polar Decomposition Method for the Approximation of the Matrix Exponential

Generalized polar decomposition method (or brieﬂy GPD method) has been introduced by Munthe-Kaas and Zanna [5] to approximate the matrix exponential. In this paper, we investigate the numerical stability of that method with respect to roundoff propagation. The numerical GPD method includes two parts: splitting of a matrix Z ∈ g , a Lie algebra of matrices and computing exp( Z ) v for a vector v . We show that the former is stable provided that (cid:3) Z (cid:3) is not so large, while the latter is not stable in general except with some restrictions on the entries of the matrix Z and the vector v .


Introduction
Since exponential matrix naturally appears in the analytical solutions of many differential equations, it is important to present the numerical approximate methods for computing exponential matrices in applied problems. In many cases that the different equations are modeled on a matrix Lie group, we need that the approximated matrix to remain in the same Lie group. Most of the well-known standard methods [3] do not satisfy this condition except in some special cases, like 3 × 3 antisymmetric matrices [2].
Recently, Munthe-Kaas and Zanna have introduced a geometric numerical method, called generalized polar decomposition method, or briefly GPD method [5]. The main advantage of this method among others is that it preserves the approximated object in the mentioned Lie group. In the GPD method, the Lie group G and its Lie algebra g are splitted to simpler ones, such that the computation of exp Z will be easier in those subgroups and subalgebras. There are two numerical algorithms in GPD method. The first one (Algorithm 1) splits the matrix Z and in the second one (Algorithm 2), exp(Z)v is computed, where v is an arbitrary vector.
As far as we know, no proof of the stability of the GPD method, with respect to propagation of roundoff error, has appeared in the literature. In this paper, we study the stability of the method for the so-called polar-type splitting order-two algorithm that has been introduced in [5]. We show that Algorithm 1 which is a simpler algorithm is stable, provided the norm Z = max i,j |Z i,j | of the matrix Z is not very large, but in the complicated one, that is, Algorithm 2, the stability is controlled by entries of Z, exp Z and v = max i |v i | as the norm of the vector v.
We introduce a sufficient condition for this aim. More precisely, given an n × n matrix Z, let for j = 1, 2, . . . , n − 1 and let then Algorithm 2 is stable under the following conditions: Z and v are not very large, (1.1a) where u is the unit roundoff error. The content of the paper is as follows. In Section 2, some preliminaries concerning the general concept of error analysis and GPD method and also details of the above-mentioned algorithms are given. In Section 3, the error analysis of the Algorithm 1 and in Section 4, the error analysis of Algorithm 2 are addressed. Section 5 is devoted to sensibility analysis of Algorithm 2.

Error analysis
In this section, we recall some standard definitions and lemmas, which are applied in the other sections. Elementary operators like +, ·, ×, and / as well as √ , sin, and cos admit the following floating point arithmetic: where op is an elementary operation, and u is the unit roundoff error.
Lemma 1 (see [1]). If |δ i | ≤ u and ρ i = ±1 for 1 ≤ i ≤ n and nu < 1, then The following technical properties are satisfied by the sequence {γn}n: Denoting the calculated value of any variable x by x and |A| = [|a ij |] for any matrix A = [a ij ], we have the following lemma.
Also for matrix multiplication, we have For any two vectors a, b ∈ R n , let (2.6) So, using this notation, we can write Similar to Lemma 1, one can prove the following lemma.

GPD method
Suppose G is a finite dimensional Lie group and g is its Lie algebra and let σ : G → G, σ = id be an involutive automorphism, that is, one-to-one smooth map such that: It is well known that for sufficiently small t and for any z = exp(tZ) ∈ G, where Z ∈ g, we can write where σ(x) = x −1 and σ(y) = y. The decomposition (2.9) is called the Generalized Polar Decomposition (GPD) of z. On the Lie algebra g, the involutive automorphism is induced by σ: hence dσ defines a splitting of g into the direct sum of two linear spaces: In fact, k is a Lie subalgebra of g, but p only admits Lie triple system property: B] is the standard commutator on Lie algebra g. Now, let where the maps Πp : g → p and Π k : g → k are the canonical projection maps. It can be easily verified that every Z ∈ g can be uniquely decomposed into Z = P + K, where (see [5]). Also, k and p satisfy If, in the decomposition (2.9), we consider x = exp X(t) and y = exp Y (t), where X(t) ∈ p and Y (t) ∈ k, for sufficiently small t, then it can be written as follows: where the coefficients X i and Y i can be found in [5]. The first terms in the expansions of X(t) and Y (t) are (2.11) Considering σ 1 , σ 2 , . . . , σm as elements of a sequence of involutive automorphisms on G, σ 1 induces a decomposition g = p 1 ⊕ k 1 and approximation: where X [1] and Y [1] are suitable truncations of (2.10) and (2.11), respectively. Then k 2 can be decomposed by σ 2 and so we have k 1 = p 2 ⊕ k 2 and exp Y [1] (t) ≈ exp X [2] (t) exp Y [2] (t) .
Repeating this process m times, (m 1) yields the following: (2.12) The number m is chosen appropriately such that the right-hand side terms can be easily computed. If the expansions (2.10) are truncated at order two, then the two following algorithms based on the iterated generalized polar decomposition (2.12) will be obtained. The first algorithm splits the matrix Z, while the second one approximates Algorithm 1 (Polar-type splitting, order two). % Purpose: 2nd order approximation of the splitting (2.12) % In: n × n matrix Z % Out: Z overwritten with the nonzero elements of X [i] and Y [m] as: for j = 1 : n − 1 a j := Z(j + 1 : n, j) b j := Z(j, j + 1 : n) T K j := Z(j + 1 : n, j + 1 : n) In Algorithm 1, we denote by Z simultaneously the input data matrix and the splitted output matrix. In this algorithm, the calculated value of c j (1 ≤ j ≤ n − 1) is Hence, By applying the norm of the matrix Z this reduces to Now, similarly for computing d j , we have Therefore, the entries of the matrix Z change as follows: where |u i j | ≤ u1, (i = 4, 5) and So from (3.1), If we replace the value of c j from Algorithm 1 in the above relation, we will have and finally for nu < 1 and from property (2.1), the upper bound of p j will be

Journal of Generalized Lie Theory and Applications
The greatest value of the right-hand side occurs at j = 1, hence for all j, However, since the diagonal entries do not change during this algorithm, (3.3) and (3.4) show the stability of the algorithm whenever Z is not so large.

Error analysis of Algorithm 2
In the GPD method, the matrix Z is splitted in Algorithm 1, and the exp(Z)v is calculated in Algorithm 2. We denote the splitted matrix of Z with Z = (z i,j ).

Error analysis of α j
Let k j = b T j a j and α j = k j , here we have assumed k j > 0. If k j is not a positive number, by considering α j = |k j | when k j < 0 and by letting sinh αj αj = 1 when k j = 0, it is easily seen that of all the arguments and the upper bounds presented in the paper remain valid.
From Lemma 2, we have where |θ where |t| ≤ (n − j)γ n−j Z 2 . (4.1) Then by assuming Z ≤ 1 and from condition (1.1), we have Hence k j + t is nonnegative and where δ, with |δ| ≤ u, denotes the floating point arithmetic error. Therefore, where where By using the identity sinh x − sinh y = 2 sinh x−y 2 cosh x+y 2 , we obtain where the quantity ε 1 j is (4.5) From (4.1) and (4.3), we have (4.6) We now consider the condition that is clearly a consequence of which itself is obtained from condition (1.1b) for Z < 1. Hence from (4.7) and by increasing monotonicity of the map x → x 1−x on the interval (0, 1), we can write The function has the elementary property that 1 ≤ f (a) ≤ f (b) whenever |a| ≤ b, from which and (4.5) we deduce |ε 1 j | ≤ η j , where (4.8) Now by using (4.4) and (4.7) in (4.8), we get

Conclusion
As we expected, the upper bound of |error j | (4.23) strongly depends on the input data Z and v. More precisely, if k 0 u, then this error bound can be controlled by max{ue Z , v }. Finally, we have shown the stability of the GPD method under condition (1.1).