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) 4.898052
## mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3.808992
## lq mean median uq max neval cld
## 4.976044 6.081604 5.293864 6.719299 15.590539 100 a
## 3.862859 4.165044 3.925445 4.211289 6.300621 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)) 81.661
## mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 5781.571
## mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3976.651
## lq mean median uq max neval cld
## 205.3415 492.5964 261.2015 311.195 12650.10 100 a
## 13190.7880 15189.2731 13696.1925 15173.957 41892.14 100 b
## 11062.0900 12206.9655 11317.6515 12046.051 34608.32 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)) 167.605
## mvnpdf(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 4933.151
## mvnpdfsmart(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 3853.252
## mvnpdfoptim(x = matrix(1.96, nrow = 2, ncol = n), Log = FALSE) 2316.105
## lq mean median uq max neval cld
## 198.0895 395.4377 245.9185 306.656 2941.054 100 a
## 12889.4015 13346.1281 13329.9065 13728.991 27165.915 100 b
## 10818.0335 11112.3187 11159.0380 11284.658 29226.204 100 c
## 6748.6690 6634.5182 6888.6625 6984.056 16913.422 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.