Asked by Ivan Abraham
on 21 Feb 2017

Hello,

I am trying to simulate a stochastic equation and empirically determine it's pdf. I have the following code,

NN=100;

ind = @(x) (x>0);

x_vect=zeros(NN,1001);

for j=1:NN

ts = 0.01;

x0=0;

x_vect(j,1)=x0;

if (mod(j,1000)==0)

rng('shuffle')

end

for i=1:100000

w = rand(1,1);

if ~(w < ts)

if (x0-2*ind(x0)*ts>0)

x_vect(j,i+1)=x0-2*ind(x0)*ts;

else

x_vect(j,i+1)=0;

end

else

x_vect(j,i+1)=x0+1;

end

x0=x_vect(j,i+1);

end

end

Each sample path should last for 1000 seconds (so the inner loop runs till 10000) and I need to do a total of NN=100000 sample paths. Currently this would take me about 7 hours to run the simulation. Does anyone see an obvious way to speed up the simulation and/or parallelize this code?

Answer by Peter
on 21 Feb 2017

Accepted Answer

Taking a quick look. There's probably more you can do.

1. Eliminating the anonymous function call on the conditional inside the inner loop should save you some time since it has a healthy amount of overhead (unless Mathworks did some magic in 2016b's JIT). I'd see what inlining it right under the ~(w < ts) does for speed.

2. You're going against MATLAB's memory indexing, which is column major order. Since your inner loop i is striding columns you're paging RAM like nobody's business. Can you rearrange the data so that the inner loop is the one looping on the rows? You want to go down rows on your inner loop. I think you can just do a transpose on how you're storing your data in this case.

3. If you're only altering the seed every 1000 rows, why not grab your rand() call then rather than call it inside the inner loop 100k times? If you've got the RAM to spare. rand is fairly quick, so I'm not positive this will outweigh the memory hit. Depends on the size of the vector, I'd wager.

4. It doesn't seem that your j's are interacting with previous ones (which makes sense since you're doing this stochastically). So, MATLAB should no problem treating this as a parfor loop. The outer loop should be the parfor one.

Ivan Abraham
on 22 Feb 2017

Sign in to comment.

Answer by Jan
on 21 Feb 2017

Edited by Jan
on 21 Feb 2017

- x_vect is pre-allocated to the final size.
- anonymous function ind is replaced by direct comparison.
- x0 - 2 * (x0 > 0) * ts is calculated once only.
- x_vect accessed in columnwise order.

NN = 100;

ni = 100000;

x_vect = zeros(ni+1, NN);

for j=1:NN

ts = 0.01;

x0 = 0;

x_vect(j,1) = x0;

if mod(j,1000) == 0

rng('shuffle')

end

for i = 2:ni+1

w = rand;

if w >= ts

v = x0 - 2 * (x0 > 0) * ts;

if v > 0

x_vect(i, j) = v;

end

else

x_vect(i, j) = x0 + 1;

end

x0 = x_vect(i, j);

end

end

This give a boost from 35.1 sec to 0.60 sec on my machine. Now start experiments with parfor.

Ivan Abraham
on 22 Feb 2017

Jan
on 22 Feb 2017

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.