Speed improvement of the random generator

Hi,

I've built a Monte Carlo program, very simple, and I want to improve the CPU time because as soon as I increase the number of simulations, the time exponentially exploses. I've run a profiler on it, and suprisingly, most of the time (>60%) is spent in the generation of the random numbers. If there is a way to reduce this computing time, it would improve a lot the speed of my program.

Any suggestion ?

Answer 1

[the cyclist]

It would be best if you could post a snippet of code that illustrates the problem.

My guess is that it is not really the random number generation that is slow. Could it be that, instead, you are growing an array by appending a random number to it every iteration of a loop? Preallocating that array would speed it up.

Also, could you pregenerate all your random numbers first, in a vectorized way, then access them later as you need them?

[Jan]

I agree: exponential slowdown is often a hint to a missing pre-allocation. Try this:

tic; x=[]; for i=1:1e5; x(i)=rand; end; toc

tic; x=zeros(1,1e5); for i=1:1e5; x(i)=rand; end; toc

[Andrew Newell]

If the number of iterations isn't too large, you might save even more time using

tic; x = rand(1,1e5); toc

Answer 2

[Oleg Komarov]

Even generating all the random numbers at once (moving it outside of the loops the gain is small), but the difference lies in the legacy mode:

k      = 20;
NbTraj = 10000;
NbPas  = 100;

tic
s      = RandStream('mcg16807','Seed',100);
dW     = randn(s,NbTraj,5,NbPas,20);
toc % Elapsed time is 8.288121 seconds.

% No legacy mode
tic
dW     = randn(NbTraj,5,NbPas,20);
toc % Elapsed time is 2.695166 seconds.

[Guillaume A.]

well, it makes no difference on my computer :

Elapsed time is 10.211157 seconds. (part I)

Elapsed time is 10.012667 seconds. (part II)

[Oleg Komarov]

Even if you clear the seed before executing the second part?

Answer 3

[Guillaume A.]

Ok, here is a sample code and the profiler result :

s = RandStream('mcg16807','Seed',100);
NbTraj=10000;
NbPas=100;
dW1 = zeros(NbTraj,NbPas);
dW2 = zeros(NbTraj,NbPas);
dW3 = zeros(NbTraj,NbPas);
dW4 = zeros(NbTraj,NbPas);
dW5 = zeros(NbTraj,NbPas);
for k=1:20      
    for i=1:NbPas        
      dW = randn(s,NbTraj,5);
      dW = dW-repmat(mean(dW,1),NbTraj,1);
      MatCov_rand = cov(dW);
      dW = dW*inv(chol(MatCov_rand))*A_theo;
      dW1(:,i) = dW(:,1);
      dW2(:,i) = dW(:,2);
      dW3(:,i) = dW(:,3);
      dW4(:,i) = dW(:,4);
      dW5(:,i) = dW(:,5);
    end
    %...
    %Do some calculus on dW1, dW2,...
    %...
end

And the result of the profiler :

RandStream.randn (MEX-function) 2000 10,559s

cov 2000 0,988s

repmat 2000 0,579s

mean 2021 0,245s

[Matt Fig]

Also needed is the sizes of the arrays involved. Perhaps put a call to the WHOS function after these loops then show the output.

Answer 4

[the cyclist]

Crap. I accidentally just deleted a whole answer instead of a comment I made. Sorry!

One more suggestion. You should be able to pull the random number generation outside of at least one of those loops, assuming it does not use too much memory to do so. That should give you some more speedup.

[Guillaume A.]

you're right and I had the same idea. However, this is just a sample code of just a "few" iterations. The program takes only a few seconds to execute, I can live with it. But I would like to increase the number of simulations to NbTraj=2e6 and k=50000 !!! A rapid calculus suggests me that it could take..... 37days ! and about 3Go of RAM memory.

So pull the random number generation outside one or more loop is not a solution :-(

But thanks so much for your effort lol !!

[the cyclist]

Be aware that your code, with the parameters you have put in, is generating 100,000,000 random normals (in about 10 seconds on my machine). That is by no means "slow". Fully vectorized [r=randn(1.e8,1)] takes 7 seconds. So, I would say you need to find another way!

[Guillaume A.]

It is not slow, but increase NbTraj to 2e6 and k to 50e3, and it could become an issue ;-)

Answer 5

[Guillaume A.]

Oleg, thanks for your answer, but unfortunately it does not change anything....