An Introduction to mvnfast
Introduction
The mvnfast R package provides computationally efficient
tools related to the multivariate normal and Student’s t distributions.
The tools are generally faster than those provided by other packages,
thanks to the use of C++ code through the
Rcpp\RcppArmadillo packages and
parallelization through the OpenMP API. The most important
functions are:
rmvn(): simulates multivariate normal random vectors.rmvt(): simulates Student’s t normal random vectors.dmvn(): evaluates the probability density function of a multivariate normal distribution.
dmvt(): evaluates the probability density function of a multivariate Student’s t distribution.
maha(): evaluates mahalanobis distances.
In the following sections we will benchmark each function against equivalent functions provided by other packages, while in the final section we provide an example application.
Simulating multivariate normal or Student’s t random vectors
Simulating multivariate normal random variables is an essential step
in many Monte Carlo algorithms (such as MCMC or Particle Filters), hence
this operations has to be as fast as possible. Here we compare the
rmvn function with the equivalent function
rmvnorm (from the mvtnorm package) and
mvrnorm (from the MASS package). In
particular, we simulate \(10^4\)
twenty-dimensional random vectors:
# microbenchmark does not work on all platforms, hence we need this small wrapper
microwrapper <- function(..., times = 100L){
ok <- "microbenchmark" %in% rownames(installed.packages())
if( ok ){
library("microbenchmark")
microbenchmark(list = match.call(expand.dots = FALSE)$..., times = times)
}else{
message("microbenchmark package is not installed")
return( invisible(NULL) )
}
}
library("mvtnorm")
library("mvnfast")
library("MASS")
# We might also need to turn off BLAS parallelism
library("RhpcBLASctl")
blas_set_num_threads(1)
N <- 10000
d <- 20
# Creating mean and covariance matrix
mu <- 1:d
tmp <- matrix(rnorm(d^2), d, d)
mcov <- tcrossprod(tmp, tmp)
microwrapper(rmvn(N, mu, mcov, ncores = 2),
rmvn(N, mu, mcov),
rmvnorm(N, mu, mcov),
mvrnorm(N, mu, mcov))## Unit: milliseconds
## expr min lq mean median uq
## rmvn(N, mu, mcov, ncores = 2) 1.254560 1.322246 1.560605 1.393784 1.454910
## rmvn(N, mu, mcov) 2.293801 2.323354 2.677443 2.361872 2.436003
## rmvnorm(N, mu, mcov) 5.111312 5.228837 5.859860 5.389511 5.568271
## mvrnorm(N, mu, mcov) 4.861119 4.961914 5.810931 5.088143 6.775554
## max neval
## 4.177810 100
## 6.117382 100
## 11.983339 100
## 13.722443 100In this example rmvn cuts the computational time,
relative to the alternatives, even when a single core is used. This gain
is attributable to several factors: the use of C++ code and efficient
numerical algorithms to simulate the random variables. Parallelizing the
computation on two cores gives another appreciable speed-up. To be fair,
it is necessary to point out that rmvnorm and
mvrnorm have many more safety check on the user’s input
than rmvn. This is true also for the functions described in
the next sections.
Notice that this function does not use one of the Random Number Generators (RNGs) provided by R, but one of the parallel cryptographic RNGs described in (Salmon et al., 2011). It is important to point out that this RNG can safely be used in parallel, without risk of collisions between parallel sequence of random numbers, as detailed in the above reference.
We get similar performance gains when we simulate multivariate Student’s t random variables:
# Here we have a conflict between namespaces
microwrapper(mvnfast::rmvt(N, mu, mcov, df = 3, ncores = 2),
mvnfast::rmvt(N, mu, mcov, df = 3),
mvtnorm::rmvt(N, delta = mu, sigma = mcov, df = 3))## Unit: milliseconds
## expr min lq mean
## mvnfast::rmvt(N, mu, mcov, df = 3, ncores = 2) 1.744291 1.819108 1.947058
## mvnfast::rmvt(N, mu, mcov, df = 3) 3.239691 3.351984 3.583185
## mvtnorm::rmvt(N, delta = mu, sigma = mcov, df = 3) 6.477329 6.703206 7.352407
## median uq max neval
## 1.885797 1.937800 3.736913 100
## 3.397988 3.436424 5.340133 100
## 6.806560 8.416801 10.284626 100When d and N are large, and
rmvn or rmvt are called several times with the
same arguments, it would make sense to create the matrix where to store
the simulated random variable upfront. This can be done as follows:
A <- matrix(nrow = N, ncol = d)
class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric".
rmvn(N, mu, mcov, A = A) Notice that here rmvn returns NULL, not the
simulated random vectors! These can be found in the matrix provided by
the user:
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.131527 -4.991416 -3.786537 -3.985325 0.2758604
## [2,] 2.051915 1.534035 2.477572 11.413259 -0.3046329Pre-creating the matrix of random variables saves some more time:
## Unit: milliseconds
## expr min lq mean median
## rmvn(N, mu, mcov, ncores = 2, A = A) 1.126300 1.173885 1.209128 1.191101
## rmvn(N, mu, mcov, ncores = 2) 1.230015 1.289127 1.739339 1.313038
## uq max neval
## 1.214927 1.829509 200
## 1.338932 60.565865 200Don’t look at the median time here, the mean is much more affected by memory re-allocation.
Evaluating the multivariate normal and Student’s t densities
Here we compare the dmvn function, which evaluates the
multivariate normal density, with the equivalent function
dmvtnorm (from the mvtnorm package). In
particular we evaluate the log-density of \(10^4\) twenty-dimensional random
vectors:
# Generating random vectors
N <- 10000
d <- 20
mu <- 1:d
tmp <- matrix(rnorm(d^2), d, d)
mcov <- tcrossprod(tmp, tmp)
X <- rmvn(N, mu, mcov)
microwrapper(dmvn(X, mu, mcov, ncores = 2, log = T),
dmvn(X, mu, mcov, log = T),
dmvnorm(X, mu, mcov, log = T), times = 500)## Unit: microseconds
## expr min lq mean median
## dmvn(X, mu, mcov, ncores = 2, log = T) 719.693 756.012 842.0439 772.5975
## dmvn(X, mu, mcov, log = T) 1329.563 1389.692 1492.6464 1414.2540
## dmvnorm(X, mu, mcov, log = T) 1228.782 1292.753 1809.1898 1338.5660
## uq max neval
## 811.140 2908.377 500
## 1449.842 6027.709 500
## 1444.189 57789.595 500Again, we get some speed-up using C++ code and some more from the parallelization. We get similar results if we use a multivariate Student’s t density:
# We have a namespace conflict
microwrapper(mvnfast::dmvt(X, mu, mcov, df = 4, ncores = 2, log = T),
mvnfast::dmvt(X, mu, mcov, df = 4, log = T),
mvtnorm::dmvt(X, delta = mu, sigma = mcov, df = 4, log = T), times = 500)## Unit: microseconds
## expr min lq
## mvnfast::dmvt(X, mu, mcov, df = 4, ncores = 2, log = T) 797.520 849.882
## mvnfast::dmvt(X, mu, mcov, df = 4, log = T) 1412.766 1482.366
## mvtnorm::dmvt(X, delta = mu, sigma = mcov, df = 4, log = T) 1320.679 1397.113
## mean median uq max neval
## 910.3674 881.700 917.498 1969.957 500
## 1531.0674 1504.403 1538.299 3416.154 500
## 1721.9453 1441.168 1515.219 4725.116 500Evaluating the Mahalanobis distance
Finally, we compare the maha function, which evaluates
the square mahalanobis
distance with the equivalent function mahalanobis (from
the stats package). Also in the case we use \(10^4\) twenty-dimensional random
vectors:
# Generating random vectors
N <- 10000
d <- 20
mu <- 1:d
tmp <- matrix(rnorm(d^2), d, d)
mcov <- tcrossprod(tmp, tmp)
X <- rmvn(N, mu, mcov)
microwrapper(maha(X, mu, mcov, ncores = 2),
maha(X, mu, mcov),
mahalanobis(X, mu, mcov))## Unit: microseconds
## expr min lq mean median uq
## maha(X, mu, mcov, ncores = 2) 706.243 727.7905 820.1541 745.1215 763.038
## maha(X, mu, mcov) 1331.997 1356.7430 1454.5460 1372.4065 1396.622
## mahalanobis(X, mu, mcov) 1233.680 1285.1465 2149.8402 1355.6865 1413.408
## max neval
## 2663.873 100
## 3449.083 100
## 58295.178 100The acceleration is similar to that obtained in the previous sections.
Example: mean-shift mode seeking algorithm
As an example application of the dmvn function, we
implemented the mean-shift mode
seeking algorithm. This procedure can be used to find the mode or
maxima of a kernel density function, and it can be used to set up
clustering algorithms. Here we simulate \(10^4\) d-dimensional random vectors from
mixture of normal distributions:
set.seed(5135)
N <- 10000
d <- 2
mu1 <- c(0, 0); mu2 <- c(2, 3)
Cov1 <- matrix(c(1, 0, 0, 2), 2, 2)
Cov2 <- matrix(c(1, -0.9, -0.9, 1), 2, 2)
bin <- rbinom(N, 1, 0.5)
X <- bin * rmvn(N, mu1, Cov1) + (!bin) * rmvn(N, mu2, Cov2)Finally, we plot the resulting probability density and, starting from 10 initial points, we use mean-shift to converge to the nearest mode:
# Plotting
np <- 100
xvals <- seq(min(X[ , 1]), max(X[ , 1]), length.out = np)
yvals <- seq(min(X[ , 2]), max(X[ , 2]), length.out = np)
theGrid <- expand.grid(xvals, yvals)
theGrid <- as.matrix(theGrid)
dens <- dmixn(theGrid,
mu = rbind(mu1, mu2),
sigma = list(Cov1, Cov2),
w = rep(1, 2)/2)
plot(X[ , 1], X[ , 2], pch = '.', lwd = 0.01, col = 3)
contour(x = xvals, y = yvals, z = matrix(dens, np, np),
levels = c(0.002, 0.01, 0.02, 0.04, 0.08, 0.15 ), add = TRUE, lwd = 2)
# Mean-shift
library(plyr)
inits <- matrix(c(-2, 2, 0, 3, 4, 3, 2, 5, 2, -3, 2, 2, 0, 2, 3, 0, 0, -4, -2, 6),
10, 2, byrow = TRUE)
traj <- alply(inits,
1,
function(input)
ms(X = X,
init = input,
H = 0.05 * cov(X),
ncores = 2,
store = TRUE)$traj
)
invisible( lapply(traj,
function(input){
lines(input[ , 1], input[ , 2], col = 2, lwd = 1.5)
points(tail(input[ , 1]), tail(input[ , 2]))
}))
As we can see from the plot, each initial point leads one of two points
that are very close to the true mode. Notice that the bandwidth for the
kernel density estimator was chosen by trial-and-error, and less
arbitrary choices are certainly possible in real applications.
References
Dirk Eddelbuettel and Romain Francois (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. URL https://www.jstatsoft.org/v40/i08/.
Eddelbuettel, Dirk (2013) Seamless R and C++ Integration with Rcpp. Springer, New York. ISBN 978-1-4614-6867-7.
Dirk Eddelbuettel, Conrad Sanderson (2014). RcppArmadillo: Accelerating R with high-performance C++ linear algebra. Computational Statistics and Data Analysis, Volume 71, March 2014, pages 1054-1063. URL https://dx.doi.org/10.1016/j.csda.2013.02.005
John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw (2011). Parallel Random Numbers: As Easy as 1, 2, 3. D. E. Shaw Research, New York, NY 10036, USA.