Parallel-izing, speeding up simulation

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 1

[Peter O]

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]

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.

Answer 2

[Jan]

• 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]

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]

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.