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