Weeks' Method for the Matrix Exponential
The accurate and reliable computation of the matrix exponential
function is a long standing problem of numerical analysis. One
approach is to compute the matrix function from a numerical inverse
Laplace transform. The numerical inverse Laplace transform is however
an inherently sensitive procedure and thus requires special
Analytic Inversion of the Laplace Transform
To understand numerical inversion of the Laplace transform, it is first necessary to understand the analytic
inversion. The latter is a well-known application of the theory of complex
variables. In general, for a time domain scalar function with
The inversion of the Laplace transform is not a unique operation. Two time domain functions which differ at a single point in time for example will have the same transform. However, due to Lerch's theorem, one can in practice assume that the analytic inverse is well-defined. Furthermore, since in applications to physical problems the time domain functions are generally continuous, the selection of a particular solution is not a fundamental problem.
Numerical Inversion Methods
The numerical inversion of the Laplace transform
One interpretation of the various numerical approaches to the inversion integral are that they are different regularization techniques. From experimentation and review, four main algorithms for numerical Laplace transform inversion have proved to be of use, the Post-Widder formula ABATE2004, VALKO2001, VALKO2004, Fourier Series Expansion DEHOOG1982, Talbot's method TALBOT1979, and the Weeks method ABATE1996, WEEKS1966, WEIDEMAN1999. Each has found a domain of application corresponding to the ability of the algorithm to invert certain classes of Laplace space functions. Additional information about the first three numerical inversion algorithms can be found in the following links:
The Weeks Method
Weeks method is one of the most well known algorithms for the numerical inversion of a scalar Laplace space function. It is also a relatively straight forward approach, going back to ideas from Lanczos in 1956 and William T. Weeks at IBM in 1966 ABATE1996, DAVIESB2002, GIUNTA1988, WEIDEMAN1999. Its popularity is due primarily to the fact that it returns an explicit expression for the time domain function.
In particular, Weeks method assumes that a smooth function of bounded exponential growth , given by the inverse Laplace transform
The two free scaling parameters and in the expansion must be selected according to the constraints that and , where is the abscissa of convergence. The restriction of to positive values ensures that the weighted Laguerre polynomials are well behaved for large t. More importantly, this condition implies that . Algorithms for the automated selection of these two tuning parameters is the main difficulty in applying the method. This issue is addressed below and in more detail in the works by WEIDEMAN1999 for scalar functions and KANO2005 for the matrix exponential.
Computation of the Coefficients
The computation of the coefficients for a general vector or matrix function
begins with an integration in the complex plane
If one chooses the Bromwich contour with ,
and assumes the expansion
With the M\"obius transformation and using along the Bromwich contour, equation
Cauchy's integral theorem ARFKEN2001 provides a method to compute . Since the function
is analytic inside the radius of convergence , the integration can be performed along the unit circle
If there is sufficient memory to store the function evaluated along the contour then a large enhancement in speed can be obtained by using the fast Fourier transform WEIDEMAN1999. The summation with can be performed in operations as opposed to required by the direct application of the midpoint rule. The coefficients corresponding to negative indices which result from the FFT algorithm are not needed for the summation and thus can be neglected. For large matrices however, it may be too demanding on memory to store a series of evaluations of . In this case, it is more efficient to overwrite the Laplace space function after each evaluation at a point in the contour. Because the midpoint rule converges quickly for a periodic function, a reasonable choice of is .
Parameter Selection AlgorithmsAs stated previously, the primary practical difficulty of applying the Weeks method is the selection of appropriate values for the methods two tuning parameters and . Through experimentation, three main automated algorithms for their selection have proven to be of use. These are
MinMax Grid AlgorithmIn this algorithm, roundoff and discretization errors are neglected. Instead, considering the expression for the truncation error from the Cauchy estimate GIUNTA1988
it is surmised the maximizing the radius of convergence
where are the singularities, will lead to a minimization of the truncation error. is therefore discretized on a grid and then maximized over the computed values. This can be performed in MATLAB using only a few lines of vectorized code.
Weideman's AlgorithmIn a recent work, J.A.C. Weideman WEIDEMAN1999 has shown that the numerical inversion of a scalar function using Weeks method can be performed robustly if one minimizes an estimate for the total error. This direct approach is clearly the most expensive and a full two parameter search is probably prohibitively costly for most matrix-valued functions . This is particularly true since the truncation error is estimated by doubling the number of expansions coefficients. Fortunately, again for the matrix exponential, I can utilize the spectrum to define as the value which maximizes the radius of convergence for a given . This allows a one parameter minimization of the total error estimate over .
My last approach is a compromise between minimization of the total error used in Weideman's algorithm and the blind
maximization of the radius of convergence. I propose to relate and by a
maximization of the radius of convergence and then perform a minimization of the Cauchy estimate for the truncation error
as a function of . I thus ignore the roundoff and discretization errors. This in turn requires one to compute
In the case of the matrix exponential
The motivation for the change of norm is the efficiency with which the two-norm of the resolvent matrix
can be estimated. The quantity is well known in numerical methods and is related to
the concept of the pseudospectrum of a matrix TREFETHEN2000
One last important detail is a reasonable choice of . I have tested two approaches. In one I
choose if , else . This selection is arbitrary but leads to an overestimate of error rather than an underestimate.
A second algorithm, which avoids choosing a value for is to maximize
Matrix Exponential CodesThere are two codes; one computes the full matrix exponential , the other . Codes for the scalar Weeks method can be found at the sight by J.A.C. Weideman Scalar Case. To run the codes, download them to a directory, unzip, and start Matlab. Examples of how to use them are provided below. Additional assistance and a description of the input parameters can be obtained by "help MWeeksExpMainPub" or "help WeeksSolverMainPub".
The Matrix exponential code....
< M A T L A B > Copyright 1984-2010 The MathWorks, Inc. Using Toolbox Path Cache. Type "help toolbox_path_cache" for more info. To get started, type one of these: helpwin, helpdesk, or demo. For product information, visit www.mathworks.com. >> Amatrix = gallery('pei',3) Amatrix = 2 1 1 1 2 1 1 1 2 >> Matlabans = expmdemo1(Amatrix) Matlabans = 20.0116 17.2933 17.2933 17.2933 20.0116 17.2933 17.2933 17.2933 20.0116 >> [WeeksAns,Sigma,Beta] = MWeeksExpMainPub(Amatrix,16,1,10,0.1,0.1,10,0.1,0,1) WeeksAns = 20.01157123002201 + 0.00000000000007i 17.29328940156454 + 0.00000000000008i 17.29328940156454 + 0.00000000000008i 17.29328940156453 + 0.00000000000008i 20.01157123002201 + 0.00000000000007i 17.29328940156454 + 0.00000000000008i 17.29328940156455 + 0.00000000000008i 17.29328940156455 + 0.00000000000008i 20.01157123002204 + 0.00000000000007i Sigma = 12.60000000000000 Beta = 10 >> AbsError = max(max(abs(WeeksAns-Matlabans))) AbsError = 2.829185933991094e-12 >> RelError = max(max(abs(WeeksAns-Matlabans)./abs(Matlabans))) RelError = 1.636002190384667e-13
The Weeks Method Solver computes . It requires less memory than the full matrix exponential by utilizing vector coefficients .
>> Amatrix = gallery('pei',5); >> bvector = rand(5,1); >> Matlabans = expmdemo1(Amatrix)*bvector Matlabans = 256.2639 254.3095 255.3307 255.0022 256.1040 >> [WeeksAns,Sigma,Beta] = WeeksSolverMainPub(Amatrix,bvector,16,1,10,0.1,0.1,10,0.1,0,1) WeeksAns = 1.0e+02 * 2.56263874102069 + 0.00000000000000i 2.54309454552695 + 0.00000000000000i 2.55330724099062 + 0.00000000000000i 2.55002192245100 + 0.00000000000000i 2.56103956714960 + 0.00000000000000i Sigma = 13.80000000000000 Beta = 10 >> AbsError = max(abs(WeeksAns-Matlabans)) AbsError = 2.361844091057682e-10 >> RelError = max(abs(WeeksAns-Matlabans)./abs(Matlabans)) RelError = 9.216453545532294e-13 Please note that in latest versions of Matlab The demos expm1, expm2, and expm3 have been renamed expmdemo1, expmdemo2, and expmdemo3, to avoid a name conflict with the new function expm1 http://www.mathworks.com/access/helpdesk/help/techdoc/ref/expm1.html For more details, see http://www.mathworks.com/access/helpdesk/help/techdoc/rn/f7-998197.html#f7-101304"
DisclaimerI can not guarantee that all the information in this website or any linked website is accurate or reliable. Also, although I have made every effort to provide accurate Matlab scripts, I can not guarantee that they are free of errors. Use at your own risk. You are free to use them or modify them as you like, please however simply provide a reference.
Website supported by NSF Grant # ITR-0325097