turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 07:04 AM

Hello,

Can you help me with an answer?

Given e.g. the following statement:

proc iml;

A=I(2);

B=expmatrix(A);

print A B;

quit;

Then the printed matrix A is I-halved instead of I. Why?

Many thanks,

Dan

Accepted Solutions

Solution

02-19-2016
04:11 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 08:24 AM

This does not happen in PROC IML, but it does happen in SAS/IML Studio.

I'm guessing that you are running this in SAS/IML Studio or got a copy of the module that is distributed with SAS/IML Studio.

You are right: The ExpMatrix function is incorrectly changing your input matrix A.

There is a simple fix: Use the following (correct) module to store a new module in SAS/IML Studio. As long as the new module is stored in a location earlier in the search path, you will pick up the new copy. This should happen automatically, unless you have changed the default search path.

In summary, run the following in SAS/IML Studio one time. From then on, you will pick up the correct behavior.

```
start ExpMatrix( A );
/* Compute a matrix exponential. Algorithm from
Golub/Van Loan, 2nd ed, p. 558. Based on
Moler & Van Loan (1978) "Nineteen Dubious Ways
to Compute the Exponential of a Matrix,"
SIAM Review 20, 801-836.
Given delta>0 and an nxn matrix A, the algorithm
computes F = exp(A+E) where
InfinityNorm(E)<= delta*InfinityNorm(A).
We use a delta of approx 1e-15. */
AX = A;
n = nrow(AX);
InfinityNorm = max( abs(AX)[+,] );
j = max( 0, 1+floor( Log2(InfinityNorm) ) );
AX = AX / 2##j;
/* Compute Pade' approximation to exp(A) */
X = AX;
c = 0.5;
M = I(n) + c*AX;
D = I(n) - c*AX;
q = 6; /* try to approximate with 6 terms (implies delta approx 1e-15)*/
do k = 2 to q;
c = c*(q-k+1)/( (2*q-k+1)*k );
X = AX*X;
M = M + c*X;
D = D + (-1)##k * c*X;
end;
F = solve(D,M);
do k = 1 to j;
F = F*F;
end;
return ( F );
finish ExpMatrix;
store module=ExpMatrix;
```

All Replies

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 07:28 AM

This works for me. Try using commas to make the output more clear.

print A, B;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 07:40 AM

thaks for the answer. it doesn't work; same effect... it seems that I should put the "print" option for A before exponetiating the matrix... quite strange...

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 08:03 AM

expmatrix is not an in-built function, it is an IML module that is read from the library. So it is possible we have different versions of the library, and that your versions modifies the matrix A. You can examine the source code of the module by entering

COPY SASHELP.IML.EXPMATRIX.SOURCE into the SAS command box.

Solution

02-19-2016
04:11 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 08:24 AM

This does not happen in PROC IML, but it does happen in SAS/IML Studio.

I'm guessing that you are running this in SAS/IML Studio or got a copy of the module that is distributed with SAS/IML Studio.

You are right: The ExpMatrix function is incorrectly changing your input matrix A.

There is a simple fix: Use the following (correct) module to store a new module in SAS/IML Studio. As long as the new module is stored in a location earlier in the search path, you will pick up the new copy. This should happen automatically, unless you have changed the default search path.

In summary, run the following in SAS/IML Studio one time. From then on, you will pick up the correct behavior.

```
start ExpMatrix( A );
/* Compute a matrix exponential. Algorithm from
Golub/Van Loan, 2nd ed, p. 558. Based on
Moler & Van Loan (1978) "Nineteen Dubious Ways
to Compute the Exponential of a Matrix,"
SIAM Review 20, 801-836.
Given delta>0 and an nxn matrix A, the algorithm
computes F = exp(A+E) where
InfinityNorm(E)<= delta*InfinityNorm(A).
We use a delta of approx 1e-15. */
AX = A;
n = nrow(AX);
InfinityNorm = max( abs(AX)[+,] );
j = max( 0, 1+floor( Log2(InfinityNorm) ) );
AX = AX / 2##j;
/* Compute Pade' approximation to exp(A) */
X = AX;
c = 0.5;
M = I(n) + c*AX;
D = I(n) - c*AX;
q = 6; /* try to approximate with 6 terms (implies delta approx 1e-15)*/
do k = 2 to q;
c = c*(q-k+1)/( (2*q-k+1)*k );
X = AX*X;
M = M + c*X;
D = D + (-1)##k * c*X;
end;
F = solve(D,M);
do k = 1 to j;
F = F*F;
end;
return ( F );
finish ExpMatrix;
store module=ExpMatrix;
```

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 08:43 AM

Yep

Thanks for your help Rick; a brilliant answer; it works !!

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 08:54 AM

In order to avoid the expmatrix function, I could tried to approximate the exponential matrix by writing in IML the Taylor expansion (summation up to e.g.100), but I would prefer your option. Thanks Rick!

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 09:01 AM

I should have mentioned an even simpler solution, pass in a COPY of your matrix instead of the matrix itself!

start MyExpMatrix(A);

A2 = A; /* copy it; A2 gets changed, but A does not */

return( ExpMatrix(A2) );

finish;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

10-21-2015 09:17 AM

Many thanks Rick.