Matrix Padé-Type Method for Computing the Matrix Exponential

Matrix Padé approximation is a widely used method for computing matrix functions. In this paper, we apply matrix Padé-type approximation instead of typical Padé approximation to computing the matrix exponential. In our approach the scaling and squaring method is also used to make the approximant more accurate. We present two algorithms for computing A e and for computing At e with many 0 t  respectively. Numerical experiments comparing the proposed method with other existing methods which are MATLAB’s functions expm and funm show that our approach is also very effective and reliable for computing the matrix exponential A e . Moreover, there are two main advantages of our approach. One is that there is no inverse of a matrix required in this method. The other is that this method is more convenient when computing At e for a fixed matrix A with many 0 t  .


Introduction
The matrix exponential is the most studied and the most widely used matrix function.It plays a fundamental role in the solution of differential equations in many application problems such as nuclear magnetic resonance spectroscopy, Markov models, and Control theory.These problems often require the computation of At e for a fixed matrix and many 0 t  For example, with some suitable assumptions on the smooth of a function f , the solution to the inhomogeneous system     , , 0 , , ,

t A t s At y t e c e f s y ds
A great number of methods for computing the matrix exponential have been proposed in the past decades.Moler and Van Loan's classic papers [1,2], studied nineteen dubious methods for computing the matrix exponential.Higham [3] which improved the traditional scaling and squaring method by a new backward error analysis is by far one of the most effective methods for com-puting the matrix exponential and his method has been implemented in MATLAB's expm function.A revisited version can be seen in [4].In Higham's method and many other methods, Padé approximation is used to evaluate the matrix exponential because of its high accuracy.
In this paper we present a new approach to compute the matrix exponential.On one hand, we use the matrix Padé-type approximation given by Gu [5] instead of the typical Padé approximation to evaluate the matrix exponential.Padé-type approximation was first proposed by Brezinski [6,7] in the scalar case.Draux [8] explored this method for the case of matrix-valued function and proposed a matrix Padé-type approximation.On the other hand, the powerful scaling and squaring method is also applied in our approach.This method exploits the relation   .Thus we need to evaluate 2 s A e at first and then take squarings for s times to obtain the evaluation of A e .In the spirit of [3], the scaling number s is determined by backward error analysis.Besides, the problem of overscaling is also considered in our approach.
We briefly introduce the definition and basic property of the matrix Padé-type approximation in Section 2.
Then we develop a specific analysis for our approach in Section 3. In Section 4, two algorithms and some numerical experiments compared with other existing methods are presented.Conclusions are made in Section 5.

Matrix Padé-Type Approximation
In this section, we give a brief introduction about Matrix Padé-type approximation or MPTA for short.There are various definitions for MPTA.We are concerned with those based on Brezinski [6,7] and Gu [5].
Let ( ) f z be a given power series with n n  ma- trix coefficients, i.e.,   0 , , .
be a generalized linear function on a polynomial space P , defined by Let m v be a scalar polynomial of degree m 0 ( ) and assume that the coefficient 1 and is denoted by    .

f k m z
The approximant    .

MPTA for the Matrix Exponential
Let f be a function having the Taylor series expansion with the radius of convergence r of the form, It follows from Higham [11] that for any and z   with ( ) Az r   , where  represents the spectral radius of a matrix, According to (2.8) and (2.9) in the previous section, the [k/m] matrix Padé-type approximant to ( ) f Az can be expressed by the form where In fact, the difference between the matrix Padé-type approximants and the typical matrix Padé approximants is that the denominator of the former is a scalar polynomial, as denoted by   km q z in (3.2).In the case of the matrix exponential, At e can be defined as follows Now we give an error expression in the following theorem.
Theorem 1 Let km q and km P be expressed as (3.4) and (3.5) respectively, then ).So far, the coefficients of the denominator km q are arbitrary.These coefficients may affect the accuracy of the approximant.Sidi [12] proposed three procedures to choose these coefficients in the case of vector-valued rational approximation.We can generalize his procedures to the matrix case.We only introduce one procedure, which is called SMPE based on minimal polynomial extrapolation.In the spirit of SMPE in [12], we aim to minimize the norm of the coefficient of the first term of the numerator, i.e., the 0 i  term in the error expression (3.6).Since this term can be viewed as the major part of the error.Under the condition 1 m b  , other coefficients of km q are the solution to the following minimization problem where we choose the Frobenius norm.The corresponding inner product of two matrices A and B is defined by The next lemma can transform the minimization problem (3.7) into the solution of a linear system.
Lemma 1 The minimization problem (3.7) is equivalent to the following linear system where !l l D A l  .Proof: See Sidi [12].Different from the series form error expression in (3.6), we present another form of error bound based on Taylor series truncation error in the following theorem.
Theorem 2 Suppose f has the Taylor series expansion of the form as in ( where , 0, , j j m    are some real constants in (0,1).Proof: The result follows from a similar analysis used in [13].
R Az P Az q z  be the matrix Padé-type approximant to the matrix exponential At e , where ,  and km q and km P have the form (3.4) and (3.5) respectively and km q is a real coefficient polynomial.If where , 0, , j j m    are some real constants in (0,1).Proof: The result follows directly from Theorem 2. The result of corollary 1 shows that the MPTA to At e is feasible.In fact, we do not have to worry about the denominator of the right side in (3.11) since our numerical experiments in the next section show that all j b are much smaller than 1.Now we turn to the implementation of the scaling and squaring method of our approach.If the norm of the matrix At is not small enough, we use the scaling and squaring technique which computes the matrix exponenttial indirectly.We want to reduce the norm of At to a sufficiently small number because the function can be well approximated near the origin.Thus we need to compute   The scaling and squaring method goes back at least to Lawson [14].
The scaling integer s should be carefully determined.We choose this number based on the error analysis.Higham [3] explored a new backward error analysis of the matrix exponential.According to the analysis in [3] Proof: See the proof in [3].The result of (3.12) shows, as Higham [3] points out, that the truncation error E in the approximant to .In terms of [11] or [13], the norm of T  is bounded by the inequality Actually, according to the choice of j b according to the minimization problem (3.7), we can assume that .
We want to choose such s that the backward relative error E At does not exceed , which is the unit roundoff in IEEE double precision arithmetic.To this end, in accordance with (3.12), we require the following inequality hold, Then we define this value However, in practice, the choice of s in (3.18) may lead to overscaling, so we need to take measures to reduce the effect caused by overscaling.Therefore the method in Al-Mohy and Higham [15] can be applied.

Lemma 3 Define ( )
with the radius of convergence r and Then according to lemma 3, we have .Since if k is too small, a large s is required which is potentially dangerous and if k is too large, evaluating matrix multiplications may require excessive computations and storages.Now we need to consider the total computational cost.The computation of choosing j b consists of two parts.One part is computing the inner products in the above linear system (3.9).It needs k matrix multiplications to obtain where the powers of A have already been obtained in the previous steps.So there are no extra matrix multiplications.Finally, we should take s matrix multiplications for squaring which costs 3 2sn flops.

Algorithms and Numerical Experiments
Based on the analysis above, we present the following algorithm for computing the matrix Padé-type approximation to the matrix exponential.Algorithm 1 This algorithm evaluates the matrix exponential using Padé-type approximation together with scaling and squaring method.The algorithm is intended for IEEE double precision arithmetic.
1) Compute and Store 2 3 4) Compute the inner products in (3.9) via (3.8). 5)Solve (3.9) to obtain the coefficients of km q . 6)Compute km q via (3.4). 7) Note that the m m  linear system in (3.9) maybe ill-conditioned.To avoid this difficulty, we can use a fast and stable iterative method to solve it.The iteration will stop if it fails to converge after the maximum number of iterations.We can also modify the algorithm by resetting j b or choosing 1, 0, 0, , 1, when the system is ill-conditioned.The latter case is actually a truncated Taylor series.Now we present some numerical experiments that compare the accuracy of three methods for computing the matrix exponential.
First we took 53 test matrices obtained from MAT-LAB R2009b's gallery function in the Matrix Computation Toolbox and most of them are real and 8 8  .Figure 1 shows the relative error in the Frobenius-norm of expmpt (Algorithm 1) and MATLAB R2009b's functions expm and funm.For details about expm and funm or other algorithms for computing the matrix exponential see [12].The exact matrix exponential A e is obtained by computing a 150 order Taylor polynomial using MATLAB's Symbolic Math Toolbox.We can make the observations from Figure 1 as follows.
 MATLAB's function expm is still the best code in general.
 The proposed algorithm expmpt is also very effecttive.Precisely, expmpt is more accurate than expm in 16 examples.
 Both expm and expmpt are more reliable than funm in general.Second we concentrate on a class of 2 2   matrix of the following form , 0 where a, c is generated randomly from (-1, -0.5) and b is generated randomly from (-1000, -500).We tested 20 matrices of this form and plotted the results in Figure 2.
We can make the observations from Figure 2 as follows.
 Algorithm 1 and MATLAB's function funm perform similarly in most examples., to the function f is defined by the properties that that km p and km q are polynomials of degrees k and m, respectively, and that So Padé-type approximation requires a polynomial of higher degree to achieve the same accuracy as Padé approximation does.Another disadvantage of our approach is that each , 2, , has to be evaluated in our approach because we need them to determine the coefficients j b .In [10], Paterson-Stockmeyer (PS) method was used to evaluate the matrix polynomials.For example, a the matrix polynomial such as can be evaluated by only 7 matrix multiplications.Therefore, our approach costs more than expm if we only compute At e for one or few t.But in another point of view, our method has some advantages.Real problems usually require At e for many 0 t  , where t represents time.For example, these t are where r is not a small number.We should first divide these t into different intervals according to the criterion that points in the same interval have the same scaling integer.According to (3.20)At and take three extra matrix multiplications for each t.Therefore, when the number of t is very large, PS method may be worthless.Moreover, in our method km q is a scalar polynomial and therefore no matrix division is required except for solving the m m  linear system with m much smaller than n .Nevertheless, Padé approximation requires the matrix division     In the end of this section, we present a modified version of Algorithm 1 to compute At e for many 0 t  .Note that in the following algorithm, we preprocessed the given matrix A by the Schur decomposition to reduce the computational cost since only triangle matrices are involved after the decomposition.
Algorithm 2 This algorithm based on Algorithm e for many different t and algorithm 1 is designed to compute A e only.But we emphasize here that in each inner loop, algorithm uses the same approach as algorithm 1 to compute each At e .In addition, only Schur decomposition is added at the beginning of the algorithm to reduce the computational cost when the number t is very much.Therefore, we do not present any numerical results of algorithm 2.

Conclusions
In this paper, we develop a new approach based on matrix Padé-type approximation mixed with the scaling and squaring method to compute the matrix exponential instead of typical Padé approximation.Two numerical algorithms for computing A e and for At e with many 0 t  respectively are proposed.Our approach is established closely relative to the backward error analysis and computational considerations.Numerical results comparing the proposed algorithm with existing functions expm and funm in MATLAB have shown the proposed algorithm is effective and reliable for computing the matrix exponential.Compared with typical Padé approximation, the most significant advantages of matrix Padétype approximation lie in two aspects.One is the convenience for computing At e for a large amount of t.The other is its avoiding n n  matrix divisions, where n is the size of the given matrix A.
and (3.3), we immediately obtain the [k/m] matrix Padé-type approximant to At e as follows, At e as equivalent to a perturbation in the original matrix At.So it is a backward error analysis.If we define In order to obtain a clear error analysis for the perturbation E, we should give a practical estimation for  .Though(3.11)gives an upper bound for  , it is a rough estimation for us rather than a series truncation error of 2 s At e  do not exceed   3 O n .The other part is solving the m m  linear system (3.9) to obtain the coefficients j b .It costs 3 2m flops which is not expensive since m is much smaller than n.When we evaluate   km P At , we rewrite (3.5) as the following form

Figure 2 .
Figure 2. Normwise relative errors for MATLAB's expm, funm and Algorithm 1 on testing 2/times 2 matrices described above. MATLAB's function expm is less reliable than the other two methods.Therefore we can conclude from our experiments that the proposed method has certain value in practice.Compared with typical Padé approximation which is used in traditional methods for the matrix exponential, Padé-type approximation we used in this paper has both some advantages and disadvantages.Recall that the [k/m] Padé approximation, denoted by       km km km r x p x q x , to the function f is defined by the properties that that km p and km q are polynomials of degrees k and m, respectively, and that Note that the purposes of algorithm 1 and of algorithm 2 are different.Algorithm 2 aims to computing At