Weeks' Method for the Matrix Exponential

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:


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

As 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 Algorithm
  • Weideman's Algorithm
  • Pseudospectrum Algorithm
Each of these is based on an error analysis of the Weeks method applied to matrix valued functions, details of which can be found in KANO2005. Briefly stated however, the algorithms work as follows.
MinMax Grid Algorithm
In 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 Algorithm
In 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 .
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

There 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"

  


Disclaimer

I 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





Patrick Kano, Moysey Brio / December 16, 2005