SIAM REVIEW c° 2003 Society for Industrial and Applied Mathematics
Vol. 45, No. 1, pp. 3–000
Nineteen Dubious Ways to
Compute the Exponential of a
Matrix, Twenty-Five Years
Later¤
Cleve Molery
Charles Van Loanz
Abstract. In principle, the exponential of a matrix could be computed in many ways. Methods involv-
ing approximation theory, differential equations, the matrix eigenvalues, and the matrix
characteristic polynomial have been proposed. In practice, consideration of computational
stability and efficiency indicates that some of the methods are preferable to others, but
that none are completely satisfactory.
Most of this paper was originally published in 1978. An update, with a separate bibliog-
raphy, describes a few recent developments.
Key words. matrix, exponential, roundoff error, truncation error, condition
AMS subject classifications. 15A15, 65F15, 65F30, 65L99
PII. S0036144502418010
1. Introduction. Mathematical models of many physical, biological, and eco-
nomic processes involve systems of linear, constant coefficient ordinary differential
equations
x˙(t) = Ax(t):
Here A is a given, fixed, real or complex n-by-n matrix. A solution vector x(t) is
sought which satisfies an initial condition
x(0) = x0:
In control theory, A is known as the state companion matrix and x(t) is the system
response.
In principle, the solution is given by x(t) = etAx0 where etA can be formally
defined by the convergent power series
etA = I + tA+ t
2A2
2! + ¢ ¢ ¢ :
¤Published electronically February 3, 2003. A portion of this paper originally appeared in SIAM
Review, Volume 20, Number 4, 1978, pages 801–836.
http://www.siam.org/journals/sirev/45-1/41801.html
yThe MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 01760-2098 (moler@mathworks.com).
zDepartment of Computer Science, Cornell University, 4130 Upson Hall, Ithaca, NY 14853-7501
(cv@cs.cornell.edu).
1

2 CLEVE MOLER AND CHARLES VAN LOAN
The effective computation of this matrix function is the main topic of this survey.
We will primarily be concerned with matrices whose order n is less than a few
hundred, so that all the elements can be stored in the main memory of a contemporary
computer. Our discussion will be less germane to the type of large, sparse matrices
which occur in the method of lines for partial differential equations.
Dozens of methods for computing etA can be obtained from more or less classical
results in analysis, approximation theory, and matrix theory. Some of the methods
have been proposed as specific algorithms, while others are based on less constructive
characterizations. Our bibliography concentrates on recent papers with strong algo-
rithmic content, although we have included a fair number of references which possess
historical or theoretical interest.
In this survey we try to describe all the methods that appear to be practical, clas-
sify them into five broad categories, and assess their relative effectiveness. Actually,
each of the “methods” when completely implemented might lead to many different
computer programs which differ in various details. Moreover, these details might have
more influence on the actual performance than our gross assessment indicates. Thus,
our comments may not directly apply to particular subroutines.
In assessing the effectiveness of various algorithms we will be concerned with the
following attributes, listed in decreasing order of importance: generality, reliability,
stability, accuracy, efficiency, storage requirements, ease of use, and simplicity. We
would consider an algorithm completely satisfactory if it could be used as the basis
for a general purpose subroutine which meets the standards of quality software now
available for linear algebraic equations, matrix eigenvalues, and initial value problems
for nonlinear ordinary differential equations. By these standards, none of the algo-
rithms we know of are completely satisfactory, although some are much better than
others.
Generality means that the method is applicable to wide classes of matrices. For
example, a method which works only on matrices with distinct eigenvalues will not
be highly regarded.
When defining terms like reliability, stability and accuracy, it is important to
distinguish between the inherent sensitivity of the underlying problem and the error
properties of a particular algorithm for solving that problem. Trying to find the
inverse of a nearly singular matrix, for example, is an inherently sensitive problem.
Such problems are said to be poorly posed or badly conditioned. No algorithm working
with finite precision arithmetic can be expected to obtain a computed inverse that is
not contaminated by large errors.
An algorithm is said to be reliable if it gives some warning whenever it introduces
excessive errors. For example, Gaussian elimination without some form of pivoting is
an unreliable algorithm for inverting a matrix. Roundoff errors can be magnified by
large multipliers to the point where they can make the computed result completely
erroneous, but there is no indication of the difficulty.
An algorithm is stable if it does not introduce any more sensitivity to perturbation
than is inherent in the underlying problem. A stable algorithm produces an answer
which is exact for a problem close to the given one. A method can be stable and
still not produce accurate results if small changes in the data cause large changes in
the answer. A method can be unstable and still be reliable if the instability can be
detected. For example, Gaussian elimination with either partial or complete pivoting
must be regarded as a mildly unstable algorithm because there is a possibility that
the matrix elements will grow during the elimination and the resulting roundoff errors