Matrix exponential

[Davidova]

Hello,
I use Proc IML function in SAS Studio for matrix exponential and the procedure is successful. But when I import the matrix in Excel and call matrix exponential function the results are different from the results from SAS. Is there anybody familiar in the method behind the procedure for matrix exponential both Excel and SAS?

[IanWakeling]

Please would you show us the syntax you are using both in IML and in EXCEL?

[Davidova]

SAS Studio

PROC IML;
A={1 1, 0 0};
t=2;
X=ExpMatrix (T*A);
ExactAnswer = ( exp (t) | | t* exp(t) ) //
( 0 | | exp (t) );
Print X, ExactAnswer;
Run;

Sorry , the excel VBA macros I can not provide now.

[IanWakeling]

I don't think I can offer much help.  You can see the source code for the function ExpMatrix from the IML library.  Look for the library catalog called SASHELP.IML, and then you can double-click on ExpMatrix to view the source, it is possible that you are using a different algorithm in VBA that would account for the differences.

[Rick_SAS · Accepted solution]

Here are a few comments:

1. Your ExactAnswer computation is wrong. It looks like you copied the documentation example but then changed the matrix. The doc says, "For the matrix used in the example...," but you are not using the matrix in the example.

2. I do not know how Excel implements the matrix exponential, but the computation is notoriously difficult. Two excellent numerical analysts (Moler and van Loan) published a famous paper titled, "Nineteen dubious ways to compute the exponential of a matrix."

3. The EXPMATRIX function in IML uses a Pade approximation (reference in doc). For a 2x2 matrix, you can check the computation by applying Putzer's 2x2 formula for the matrix exponential. See p. 866-867 of Gustafson (2022) Differential Equations and Linear Algebra

The following IML program implements the 2x2 case for real eigenvalues:

proc iml;
A={1 1, 
   0 0};
t=2;

/* call the built-in EXPMATRIX module */
expMat=ExpMatrix (t*A);
Print expMat;

/* implement Putzen's 2x2 formula for real eigenvalues */
start Putzen2x2(A, t);
   eval = eigval(A);
   /* assume eigenvalues are real */
   lambda1 = eval[1,1];
   lambda2 = eval[2,1];
   r1 = exp(lambda1*t);
   r2 = (exp(lambda1*t) - exp(lambda2*t))/(lambda1-lambda2);
   expAt = r1*I(2) + r2*(A-lambda1*I(2));
   return expAt;
finish;

Putzen = Putzen2x2(A, t);
print Putzen;

Both methods report that exp(2*A) is

{7.3890561 6.3890561,
 0                  1                  }

so there is no evidence that the IML computation is not accurate.

[Davidova]

Thanks for your comments! I appreciate your support!