Skip to main content

endeav: a course on code efficiency and software development with R

4 Profiling code

Profiling is about determining which part of the code take the most time to compute (and also memory-wise). Once you have found the block of code that takes the longest time to execute, our goal is only to optimize that small part of the code.

To get a profiling of the code below, select the lines of code of interest and go to the “Profile” menu then “Profile Selected Lines”. It uses the package profvis, and in particular its profvis() function.

n <- 10e4
pdfval <- mvnpdf(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE)

OK, OK, we get it ! Concatenating a vector as you go in a loop is really not a good idea.

4.1 Comparison with an improved implementtation of mnvpdf().

Consider a new version of mvnpdf(), called mvnpdfsmart(). Download the file, and then include it in the package.

Now profile the following command:

n <- 10e4
pdfval <- mvnpdfsmart(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE)

We have indeed removed the main computational bottleneck, and we can now learn in a more detailed way what takes time in our function.

To confirm that mvnpdfsmart() is indeed much faster than mvnpdf() we can make a comparison using microbenchmark():

n <- 1000
mb <- microbenchmark(mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE),
                     mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n),
                                 Log = FALSE),
                     times=100L)
mb
## Unit: milliseconds
##                                                            expr      min
##       mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3.297466
##  mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2.297312
##        lq     mean  median       uq      max neval cld
##  3.400622 3.769656 3.50345 3.620095 7.665442   100  a 
##  2.329722 2.458756 2.35996 2.405450 6.446061   100   b

We can also check whether mvnpdfsmart() becomes competitive with dmvnorm():

n <- 1000
mb <- microbenchmark(mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2)),
                     mvnpdf(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE),
                     mvnpdfsmart(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE),
                     times=100L)
mb
## Unit: microseconds
##                                                            expr      min
##              mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2))   43.829
##       mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3252.161
##  mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2285.258
##        lq       mean    median       uq      max neval cld
##    49.815   58.54759   58.3225   63.263  104.919   100 a  
##  3481.843 3832.78045 3518.6815 3580.161 7859.126   100  b 
##  2314.019 2372.60440 2329.2305 2347.270 6178.208   100   c

There is still work to be done…

4.2 Comparison with an optimized pure implementation

After several research, tests, trials and errors, Boris arrived at an optimized version using capabilities.

Include this mvnpdfoptim() function in your package, and then profile it:

n <- 10e4
profvis::profvis(mvnpdfoptim(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE))

And the microbenchmark() that goes with it:

n <- 1000
mb <- microbenchmark(mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2)),
                     mvnpdf(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE),
                     mvnpdfsmart(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE),
                     mvnpdfoptim(x=matrix(1.96, nrow = 2, ncol = n), Log=FALSE),
                     times=100L)
mb
## Unit: microseconds
##                                                            expr      min
##              mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2))   42.558
##       mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3322.845
##  mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2309.120
##  mvnpdfoptim(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 1741.762
##         lq      mean   median        uq      max neval  cld
##    51.8445  103.3696   63.673   73.9435 4088.807   100 a   
##  3444.2050 3778.0147 3549.042 3615.4620 7582.048   100  b  
##  2333.5355 2483.4606 2354.097 2392.2475 9521.799   100   c 
##  1772.5735 1848.4469 1796.456 1834.3400 5345.211   100    d

Finally, we can profile the dmvnorm() function:

n <- 10e6
profvis::profvis(mvtnorm::dmvnorm(matrix(1.96, nrow = n, ncol = 2)))
profvis::profvis(mypkgr::my_dmvnorm(matrix(1.96, nrow = n, ncol = 2)))

You can download the my_dmvnorm() function [here](here and include it in your package, in order to have the source code available in the profile result.