How to add a force to Langevin equation with parfor?

The code below iterates over a simple Langevin equation, representing Brownian motion of 10^5 particles .

  tic 
             dt=0.1; 
             Particles=10^5; 
             x=zeros(1,Particles);  
             MaxTime=10000 ; 
             parfor ii=1:MaxTime 
             Gaussian=sqrt(pi*dt)*normrnd(0,1,[1,Particles]); 
             x=x+Gaussian;  
             end 
toc

The problem is that I don't know how to add a force term to the equation, in a way which still allows me to use parfor.

I want to do something like this:

tic 
             dt=0.1; 
             Particles=10^5; 
             x=zeros(1,Particles);  
             Force=zeros(1,Particles);  
             MaxTime=10000 ; 
             parfor ii=1:MaxTime 
             Gaussian=sqrt(pi*dt)*normrnd(0,1,[1,Particles]); 
             x=x+Force*dt+Gaussian;  
             Force=-x ; 
             end      
toc 

This gives the following error message: "The parfor cannot run due to the way variable 'x' is used".

With a standard 'for' loop of course this is not a problem.

It seems that, even though I know that Matlab is naturally paralellized, when one uses vector & matrices correctly, still for the Brownian motion

parfor seems to be 4 times faster than the standard 'for' loop.

How to do this WITH the force term?

Answer 1

[Raunak Gupta]

Hi Erez,

From the code I think parallelization is not possible in this case because if you see the Force variable is updated in every iteration and x is also dependent on the value of Force that is calculated in previous iteration. So, there is a dependency on the sequence of execution of loop. Also, if x is termed as a reduction variable it has to be initialized within parfor loop, which will not be correct as per the logic in code.

You may refer to the following documentations to get information about parfor loops.

Decide when to Use parfor

Ensure That parfor-Loop Iterations are Independent

Troubleshoot Variables in parfor-Loops

Temporary Variables

[Erez]

Hi Raunak, thanks. But why actually am I seeing a speed-up when I use parfor with Gaussian=sqrt(pi*dt)*normrnd(0,1,[1,Particles]); and x=x+Gaussian; as opposed to for?

I thought that this is rellated to parralelization with respect to the number of elements ('particles') in x, beyond what Matlab does in standard vercot coputations.

For that reason, I thought maybe Gaussian=sqrt(pi*dt)*normrnd(0,1,[1,Particles]); and x=x-x*dt+Gaussian;

(where the force is explicitly stated inside the equation for x), would at least work. It doesn't.

Is there any way at all to speed-up this kind of calculation in Matlab? Know how?

[Raunak Gupta]

Hi Erez,

I see your point about the first line of code

x = x + Gaussian;

is a straightforward reduction and x is not deemed as a temporary variable. But the second line of code

x = x - x*dt + Gaussian;

will require to declare x within the parfor as x is multiplied with the "dt" and the value of x has to be stored at that particular iteration to accurately calculate the output. Whereas in the first line of code even if the intermediate x is not known the final output will be correct. The parallelization is only successful when there is no dependency between the individual iteration otherwise one worker (parallel node) will have to wait for the output of another worker which eventually can make the code serial.

I tried doing the second line of code but it makes x as temporary variable. In my knowledge I don't think parallelization is possible in such cases because iterations are interdependent and there is a kind of scaling also applied to output of each iteration.

The 4x speed you got from reduction loop is probably because of 4 or 6 cores available in your system.

Hope this clarifies!

[Erez]

Ok, thank you. Yes, it solves most of the mystery. Matlab always uses as many cores as it can at the same time though, to parralelize jobs with vectors. So even without parfor I always see the calculation running on 4-6 cores. But I guess that using parfor cuts the code into blocks in a better way somehow, which works faster (related to Matlab's internal algorithm). So you are right, there's not much more to do beyond it in this case.