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 consideration.
The purpose of this website is to introduce a recently studied application of the Weeks method for the numerical inverse Laplace transform to the computation of the matrix exponential. Much of the material here comes from this published PAPER and the dissertation by P. Kano KANO2005. For additional introductory material see recent overview talk Numerical Laplace Transform Inversion and Selected Applications. Also provided are MATLAB scripts for computing the matrix exponential with this method.

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 Laplace transform

the inverse Laplace transform

can be evaluated by complex integration along a contour . For isolated singularities, the Bromwich contour is the standard approach. The position along the real axis is known as the abscissa of convergence and is chosen by convention to be larger than the real part of all singularities of .

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

is an inherently ill-posed problem. From a functional analysis perspective, this can property can be seen as a special case of the difficulty inherent in solving an inverse problem described by a Fredholm integral equation of the first kind DAVIESB2002. In contrast, from a numerical analysis point of view, the inherent sensitivity of the inversion procedure is clear when one considers the need to multiply by a potentially increasing large exponential . Algorithmic and finite precision errors can lead to exponential divergence of numerical solutions.

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:

• Talbot's Method
• Fourier Series Method
• Post Widder Formula Method
• The Weeks Method

The Weeks Method

Introduction

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

where is a contour in the complex plane, can be expressed as the limit of an expansion in scalar Laguerre polynomials

The functions for are defined by the equation

on . The coefficients , which may be scalars, vectors, or matrices, contain the information particular to the Laplace space function and may be complex if is complex. More importantly, these coefficients are time independent so that can be evaluated at multiple times from a single set of coefficients. Details on how to compute these coefficients are provided below.

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

then equating the two expressions yields

It is known that the weighted Laguerre coefficients have the Fourier representation

Performing the appropriate substitution, assuming it is possible to interchange the sum and integral, and equating integrands thus leaves

The functions form a complete, orthogonal basis in HIGGINS1977. In principle, one could try to use directly the orthogonality of the basis to determine . However, these functions are high oscillatory and thus not amendable to numerical integration. One method of recourse to this problem is to apply a Moebius transformation to map to a new complex variable

For the Bromwich contour

so that is mapped to the unit circle . A more subtle result of this transformation is that the singularities of in the half-plane are mapped to the exterior of the unit circle in the plane. An exception to this occurs when has a singularity at infinity which then maps onto the unit circle DAVIESB2002. In the following I assume that the singularities of occur in a finite region of the complex plane. The Laplace space function corresponding to the matrix exponential belongs to this class of functions.

With the M\"obius transformation and using along the Bromwich contour, equation

can now be expressed as


In this form it is clear that the coefficients of the original expansion

are also the coefficients of a Maclaurin series. The radius of convergence is strictly greater than unity due to my selection of functions which do not have a singularity at infinity. Furthermore, within this radius, the power series converges uniformly.

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


Numerically the evaluation of the integral can be computed very accurately using the midpoint rule
, where and

Implementation of this formula for the coefficients for vectors and matrices requires consideration, particularly, with respect to memory and computing costs.

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 Algorithms

• MinMax Grid Algorithm
• Weideman's Algorithm
• Pseudospectrum Algorithm

MinMax Grid Algorithm

Weideman's Algorithm

Pseudospectrum Algorithm

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 estimate

For some functions it may be possible to determine analytically estimates for or compute the norm without directly computing on circular contours. The matrix exponential provides an excellent example where this is true.

In the case of the matrix exponential

It is a known fact of matrix norms that

where is the dimension of the square matrix . With this, I can replace the Frobenius norm

The constant is not a function of the parameters and nor the radius and is therefore not involved in the minimization of . For this reason, I drop this factor in the following discussion.

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

where is the smallest singular value of the matrix . In particular,

Computing the full singular value decomposition of to determine the norm is hardly advantageous but there are many algorithms for computing the minimum singular value. The method used here to accelerate the computation of is the ``triangularization+inverse Lanczos iteration" by Lui LUI1997 The concept behind this method is the fact that for

where is the smallest eigenvalue. For my implementation of the Weeks method, the Schur decomposition has already been performed and the number of inversions required for the power method is often less than required by doubling the number coefficients to estimate the tail of the Laguerre expansion.

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
over .

Matrix Exponential Codes

• Weeks Method Matrix Exponential

• Weeks Method Solver

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"

  

Disclaimer

Website supported by NSF Grant # ITR-0325097