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)Flame Graph
Data
| mypkgr/R/mvnpdf.R | Memory | Time | ||||||
|---|---|---|---|---|---|---|---|---|
| #' Multivariate-Normal probability density function | ||||||||
| #' | ||||||||
| #' This is a concise description of what the function does. | ||||||||
| #' | ||||||||
| #' This part gives more details on the function. | ||||||||
| #' | ||||||||
| #' @param x a p x n data matrix with n the number of observations and | ||||||||
| #'p the number of dimensions | ||||||||
| #' @param mean mean vector | ||||||||
| #' @param varcovM variance-covariance matrix | ||||||||
| #' @param Log logical flag for returning the log of the probability density | ||||||||
| #'function. Default is \code{TRUE}. | ||||||||
| #' | ||||||||
| #' @return a list containing the input matrix x and y the multivariate-Normal probability density function | ||||||||
| #' computed at x | ||||||||
| #' | ||||||||
| #' @export | ||||||||
| #' | ||||||||
| #' @examples | ||||||||
| #' mvnpdf(x=matrix(1.96), Log=FALSE) | ||||||||
| #'dnorm(1.96) | ||||||||
| #' | ||||||||
| #'mvnpdf(x=matrix(rep(1.96, 2), nrow=2, ncol=1), Log=FALSE) | ||||||||
| mvnpdf <- function(x, mean = rep(0, nrow(x)), | ||||||||
| varcovM = diag(nrow(x)), Log = TRUE) { | ||||||||
| n <- ncol(x) | ||||||||
| p <- nrow(x) | ||||||||
| x0 <- x - mean | ||||||||
| Rinv <- solve(varcovM) | ||||||||
| LogDetvarcovM <- log(det(varcovM)) | ||||||||
| y <- NULL | ||||||||
| for (j in 1:n) { | ||||||||
| yj <- - p/2 * log(2*pi) - 0.5 * LogDetvarcovM - | ||||||||
| 0.5 * t(x0[, j]) %*% Rinv %*% x0[, j] | ||||||||
| y <- c(y, yj) | ||||||||
| } | ||||||||
| if (!Log) { | ||||||||
| y <- exp(y) | ||||||||
| } | ||||||||
| res <- list(x = x, y = y) | ||||||||
| class(res) <- "mvnpdf" | ||||||||
| return(res) | ||||||||
| } | ||||||||
| #' Plot of the mvnpdf function on the first dimension | ||||||||
| #' | ||||||||
| #' @param x an object of class \code{mvnpdf} resulting from a call of | ||||||||
| #' \code{mnvpdf()} function. | ||||||||
| #' @param ... graphical parameters passed to \code{plot()} function. | ||||||||
| #' | ||||||||
| #' @return Nothing is returned, only a plot is given. | ||||||||
| #' @export | ||||||||
| #' @importFrom graphics plot | ||||||||
| #' @examples | ||||||||
| #' pdfvalues <- mvnpdf(x=matrix(seq(-3, 3, by = 0.1), nrow = 1), Log=FALSE) | ||||||||
| #' plot(pdfvalues) | ||||||||
| plot.mvnpdf <- function(x, ...) { | ||||||||
| graphics::plot(x$x[1, ], x$y, type = "l", ylab="mvnpdf", xlab="Obs (1st dim)", ...) | ||||||||
| } | ||||||||
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)Flame Graph
Data
| mypkgr/R/mvnpdfsmart.R | Memory | Time | ||||||
|---|---|---|---|---|---|---|---|---|
| #' @rdname mvnpdf | ||||||||
| #' @export | ||||||||
| mvnpdfsmart <- function(x, mean = rep(0, nrow(x)), | ||||||||
| varcovM = diag(nrow(x)), Log = TRUE) { | ||||||||
| n <- ncol(x) | ||||||||
| p <- nrow(x) | ||||||||
| x0 <- x - mean | ||||||||
| Rinv <- solve(varcovM) | ||||||||
| LogDetvarcovM <- log(det(varcovM)) | ||||||||
| y <- rep(NA, n) | ||||||||
| for (j in 1:n) { | ||||||||
| yj <- - p/2 * log(2*pi) - 0.5 * LogDetvarcovM - | ||||||||
| 0.5 * t(x0[, j]) %*% Rinv %*% x0[, j] | ||||||||
| y[j] <- yj | ||||||||
| } | ||||||||
| if (!Log) { | ||||||||
| y <- exp(y) | ||||||||
| } | ||||||||
| res <- list(x = x, y = y) | ||||||||
| class(res) <- "mvnpdf" | ||||||||
| return(res) | ||||||||
| } | ||||||||
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:
Flame Graph
Data
| mypkgr/R/mvnpdfoptim.R | Memory | Time | ||||||
|---|---|---|---|---|---|---|---|---|
| #' @rdname mvnpdf | ||||||||
| #' @export | ||||||||
| mvnpdfoptim <- function(x, mean = rep(0, nrow(x)), | ||||||||
| varcovM = diag(nrow(x)), Log=TRUE){ | ||||||||
| if(!is.matrix(x)){ | ||||||||
| x <- matrix(x, ncol=1) | ||||||||
| } | ||||||||
| n <- ncol(x) | ||||||||
| p <- nrow(x) | ||||||||
| x0 <- x-mean | ||||||||
| Rinv <- backsolve(chol(varcovM), x=diag(p)) | ||||||||
| xRinv <- apply(X=x0, MARGIN=2, FUN=crossprod, y=Rinv) | ||||||||
| logSqrtDetvarcovM <- sum(log(diag(Rinv))) | ||||||||
| quadform <- apply(X=xRinv, MARGIN=2, FUN=crossprod) | ||||||||
| y <- (-0.5*quadform + logSqrtDetvarcovM -p*log(2*pi)/2) | ||||||||
| if(!Log){ | ||||||||
| y <- exp(y) | ||||||||
| } | ||||||||
| res <- list(x = x, y = y) | ||||||||
| class(res) <- "mvnpdf" | ||||||||
| return(res) | ||||||||
| } | ||||||||
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:
Flame Graph
Data
profvis::profvis(mypkgr::my_dmvnorm(matrix(1.96, nrow = n, ncol = 2)))Flame Graph
Data
| mypkgr/R/my_dmvnorm.R | Memory | Time | ||||||
|---|---|---|---|---|---|---|---|---|
| #'Copy of mvtnorm::dmvnorm for progiling and reference | ||||||||
| #' | ||||||||
| #'@param x matrix | ||||||||
| #'@param mean vector | ||||||||
| #'@param sigma covariance matrix | ||||||||
| #'@param log logical. Default is \code{FALSE} | ||||||||
| #'@param checkSymmetry logical. Default is \code{TRUE} | ||||||||
| #' | ||||||||
| #'@export | ||||||||
| #' | ||||||||
| my_dmvnorm <- function(x, mean = rep(0, p), sigma = diag(p), log = FALSE, | ||||||||
| checkSymmetry = TRUE){ | ||||||||
| if (is.vector(x)) | ||||||||
| x <- matrix(x, ncol = length(x)) | ||||||||
| p <- ncol(x) | ||||||||
| if (!missing(mean)) { | ||||||||
| if (!is.null(dim(mean))) | ||||||||
| dim(mean) <- NULL | ||||||||
| if (length(mean) != p) | ||||||||
| stop("x and mean have non-conforming size") | ||||||||
| } | ||||||||
| if (!missing(sigma)) { | ||||||||
| if (p != ncol(sigma)) | ||||||||
| stop("x and sigma have non-conforming size") | ||||||||
| if (checkSymmetry && !isSymmetric(sigma, tol = sqrt(.Machine$double.eps), | ||||||||
| check.attributes = FALSE)) | ||||||||
| stop("sigma must be a symmetric matrix") | ||||||||
| } | ||||||||
| dec <- tryCatch(base::chol(sigma), error = function(e) e) | ||||||||
| if (inherits(dec, "error")) { | ||||||||
| x.is.mu <- colSums(t(x) != mean) == 0 | ||||||||
| logretval <- rep.int(-Inf, nrow(x)) | ||||||||
| logretval[x.is.mu] <- Inf | ||||||||
| } | ||||||||
| else { | ||||||||
| tmp <- backsolve(dec, t(x) - mean, transpose = TRUE) | ||||||||
| rss <- colSums(tmp^2) | ||||||||
| logretval <- -sum(log(diag(dec))) - 0.5 * p * log(2 * | ||||||||
| pi) - 0.5 * rss | ||||||||
| } | ||||||||
| names(logretval) <- rownames(x) | ||||||||
| if (log) | ||||||||
| logretval | ||||||||
| else exp(logretval) | ||||||||
| } | ||||||||
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.