## Parallel-izing, speeding up simulation

### Ivan Abraham (view profile)

on 21 Feb 2017
Latest activity Edited by Jan

on 22 Feb 2017

### Peter (view profile)

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?

### Peter (view profile)

on 21 Feb 2017

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

### Ivan Abraham (view profile)

on 22 Feb 2017
That was one of my reasons why I posted the question here. In my opinion, each sample path can be generated on a separate worker, but MATLAB won't let me parfor it.

on 21 Feb 2017
Edited by Jan

### Jan (view profile)

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

### Ivan Abraham (view profile)

on 22 Feb 2017
Thank you for the answer, this helps a lot. I did not know about the access column-wise vs row-vise made a big difference.
Jan

### Jan (view profile)

on 22 Feb 2017
The acceleration is caused mainly by the pre-allocation and avoiding the anonymous function, while the columns-wise order pushed the run time by 0.4 seconds only.