CODEWITHSUNDEEP
performancer
mrkb80
mrkb80

Reputation: 591

R nested loop slow

I have no idea why something like this should be slow:

steps=500
samples=100000
s_0=2.1
r=.02
sigma=.2
k=1.9

at<-matrix(nrow=(steps+1),ncol=samples)
at[1,]=s_0

for(j in 1:samples)
{
  for(i in 2:(steps+1))
  {
    at[i,j]=at[(i-1),j] + sigma*sqrt(.0008)*rnorm(1)
  }
}

I tried to rewrite this using sapply, but it was still awful from a performance standpoint.

Am I missing something here? This would be seconds in c++ or even the bloated c#.

Upvotes: 4

Views: 1332

Answers (3)

Kevin Ushey
Kevin Ushey

Reputation: 21285

To write fast R code, you really need to re-think how you write functions. You want to operate on entire vectors, not just single observations at a time.

If you're really deadset in writing C-style loops, you could also try out Rcpp. Could be handy if you're well accustomed to C++ and prefer writing functions that way.

library(Rcpp)
do_stuff <- cppFunction('NumericMatrix do_stuff(
  int steps,
  int samples,
  double s_0,
  double r,
  double sigma,
  double k ) {

  // Ensure RNG scope set
  RNGScope scope;

  // allocate the output matrix
  NumericMatrix at( steps+1, samples );

  // fill the first row
  for( int i=0; i < at.ncol(); i++ ) {
    at(0, i) = s_0;
  }

  // loop over the matrix and do stuff
  for( int j=0; j < samples; j++ ) {
    for( int i=1; i < steps+1; i++ ) {
      at(i, j) = at(i-1, j) + sigma * sqrt(0.0008) * R::rnorm(0, 1);
    }
  }

  return at;

}')

system.time( out <- do_stuff(500, 100000, 2.1, 0.02, 0.2, 1.9) )

gives me

   user  system elapsed 
  3.205   0.092   3.297

So, if you've already got some C++ background, consider learning how to use Rcpp to map data to and from R.

Upvotes: 1

Jan
Jan

Reputation: 11726

R can vectorize certain operations. In your case you can get rid of the outer loop by doing a following change.

for(i in 2:(steps + 1))
{
    at[i,] = at[(i - 1),] + sigma * sqrt(.0008) * rnorm(samples)
}

According to system.time the original version for samples = 1000 takes 6.83s, while the modified one 0.09s.

Upvotes: 4

Ben Bolker
Ben Bolker

Reputation: 226087

How about:

at <- s_0 + t(apply(matrix(rnorm(samples*(steps+1),sd=sigma*sqrt(8e-4)),
                   ncol=samples),
                    2,
                    cumsum))

(Haven't tested this carefully yet, but I think it should be right, and much faster.)

Upvotes: 4

Related Questions