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

06-29-2016 10:14 AM

I am new to vectorizing generally and SAS/IML though have quickly realised the potential program and computational savings. I have a specific problem I have set out below and have searched the internet, to no avail. I have been inspired by the solution to another problem in this link - in which a series calculation is vectorized from a DO loop and read a lot of Rick Wicklin's blog which have helped massively.

**Problem: **I have a transition matrix (e.g. 4x4) and I want to forecast the volume movement between the 4 states over N months. The final result will be a Nx4 matrix with the forecasted volumes in each state. I can achieve this using a DO loop but I am hopeful that someone can help me find a vectorized solution, as I may have to run similar vectorize this problem.

The attached files contains an example of my current approach (using a DO loop).

Thanks in advance.

Andrew

Accepted Solutions

Solution

06-29-2016
11:27 AM

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

06-29-2016 11:24 AM

This is an iterative computation, so you can't avoid the DO loop. If you are only interested in the asymptotic steady-state solution, there are ways to find that directly by using eigenvalues, but if you want the forecast time series then you have to compute each step.

There is a way to make your program more efficient. Currently you are computing a matrix power (trans**i), which is an expensive computation at each stepo of the loop. At the i_th step you are repeating all the computations from the previous step, plus one more matrix multiplication.

However, you can reduce this to N matrix-vector computations, which is much faster. Inside the loop, merely update the previous computation, like this:

```
v = input1;
do i = 1 to ForecastMonths ;
v = v * trans;
results[i, ] = v;
end ;
```

It isn't clear whether you want the first row to be the 0th iteration or the result after the first iteration. I've done the latter.

Not only is this method faster, but it is numerically more stable because you never form large powers of matrices.

All Replies

Solution

06-29-2016
11:27 AM

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

06-29-2016 11:24 AM

This is an iterative computation, so you can't avoid the DO loop. If you are only interested in the asymptotic steady-state solution, there are ways to find that directly by using eigenvalues, but if you want the forecast time series then you have to compute each step.

There is a way to make your program more efficient. Currently you are computing a matrix power (trans**i), which is an expensive computation at each stepo of the loop. At the i_th step you are repeating all the computations from the previous step, plus one more matrix multiplication.

However, you can reduce this to N matrix-vector computations, which is much faster. Inside the loop, merely update the previous computation, like this:

```
v = input1;
do i = 1 to ForecastMonths ;
v = v * trans;
results[i, ] = v;
end ;
```

It isn't clear whether you want the first row to be the 0th iteration or the result after the first iteration. I've done the latter.

Not only is this method faster, but it is numerically more stable because you never form large powers of matrices.

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

06-29-2016 11:30 AM

Thank you very much for the prompt response. Reading your (and other) other blog entries I was perhaps being overly hopeful there was a way to vectorize everything in my life, so it is helpful to get confirmation where a loop is required.

Also, thank you for the advice regarding the vector/matrix multiplication, as opposed to the matrix-power operation.

Andrew

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

06-29-2016 01:09 PM

A good rule of thumb is to look for a vectorized approach when your loop is over the rows or columns of a matrix or vector.

In this case, the loop is an independent of the size of the matrix.

Another possible opportunity for vectorization is when you are performin an operation on a VECTOR within a loop, and the (i+1)st computation is INDEPENDENT of the result of the i_th computation. Then you might be able to pack all those vectors into a matrix and do one columnwise or rowwise operation on the matrix. Alas, that is not the case here: the (i+1)st computation depends on the result of the i_th computation.

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

06-30-2016 05:06 AM

Thanks Rick.

I understood from your other blog posts that the dependency of the (i+1)th observation on the ith means that vectorization wouldn't be possible. However, it was when someone showed me that I can get to the ith iteration by raising the transition matrix to the ith power that I thought there might be a way to do this in one step as I only require the initial vector and the transition matrix**i. However, I note your point about matrix power being more expensive than matrix multiplication.

That is disappointing. However, I'm on a mission as I'm working with some heavily looped code at the moment and there are some definite wins for me.

Cheers,

Andrew

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

07-08-2016 05:57 AM

I wrote down some of my thoughts on how to iterate a Markov transition matrix efficiently in SAS/IML. See

Markov transition matrices in SAS/IML

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

07-17-2016 06:54 AM

For future reference, I wrote up some of these ideas and related concepts for analyzing Markov chains:

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

07-18-2016 09:54 AM

Thanks for sharing these Rick. I have been following your blog with interest!