Matrix Exponentials

For background on the computation of matrix exponentials, see "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later," SIAM Review 45, 3-49, 2003.

The Pseudospectra Gateway is also highly recommended. The web site is:

http://web.comlab.ox.ac.uk/projects/pseudospectra/

Here are three of the 19 ways to compute the exponential of a matrix.

Start with the matrix A.
Scaling and squaring
Matrix exponential via Taylor series.
Matrix exponential via eigenvalues and eigenvectors.
Compare the results.
Taylor series fails.
Eigenvalues and vectors fails.

Start with the matrix A.

A = [0 1 2; 0.5 0 1; 2 1 0]
Asave = A;

A =

         0    1.0000    2.0000
    0.5000         0    1.0000
    2.0000    1.0000         0

Scaling and squaring

expmdemo1 is an M-code implementation of the built-in algorithm used by MATLAB for the matrix exponential. See Golub and Van Loan, Matrix Computations, 3rd edition, algorithm 11.3-1.

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
   c = c * (q-k+1) / (k*(2*q-k+1));
   X = A*X;
   cX = c*X;
   E = E + cX;
   if p
     D = D + cX;
   else
     D = D - cX;
   end
   p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
for k = 1:s, E = E*E; end

E1 = E

E1 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

Matrix exponential via Taylor series.

expmdemo2 uses the classic definition for the matrix exponential. As a practical numerical method, this is slow and inaccurate if norm(A) is too large.

A = Asave;

% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
   E = E + F;
   F = A*F/k;
   k = k+1;
end

E2 = E

E2 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

Matrix exponential via eigenvalues and eigenvectors.

expmdemo3 assumes that the matrix has a full set of eigenvectors. As a practical numerical method, the accuracy is determined by the condition of the eigenvector matrix.

A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E

E3 =

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

Compare the results.

For this matrix, they all do equally well

E = expm(Asave);
err1 = E - E1
err2 = E - E2
err3 = E - E3

err1 =

     0     0     0
     0     0     0
     0     0     0


err2 =

  1.0e-014 *

   -0.3553   -0.1776   -0.2665
   -0.1332   -0.1332   -0.0444
    0.1776         0    0.3553


err3 =

  1.0e-014 *

   -0.7105   -0.5329   -0.7105
   -0.6217   -0.5329   -0.5773
   -0.5329   -0.5329   -0.5329

Taylor series fails.

Here is a matrix where the terms in the Taylor series become very large before they go to zero. Consequently, expmdemo2 fails.

A = [-147 72; -192 93];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

   -0.0996    0.0747
   -0.1991    0.1494


E2 =

  1.0e+006 *

   -1.1985   -0.5908
   -2.7438   -2.0442


E3 =

   -0.0996    0.0747
   -0.1991    0.1494

Eigenvalues and vectors fails.

Here is a matrix that does not have a full set of eigenvectors. Consequently, expmdemo3 fails.

A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E2 = expmdemo2(A)
E3 = expmdemo3(A)

E1 =

    0.3679    0.3679
         0    0.3679


E2 =

    0.3679    0.3679
         0    0.3679


E3 =

    0.3679         0
         0    0.3679