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 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
## rmvn(N, mu, mcov, ncores = 2) 2.972677 3.139757 3.357267 3.218336
## rmvn(N, mu, mcov) 5.112280 5.528061 5.914561 5.720748
## rmvnorm(N, mu, mcov) 19.393903 20.552228 21.758993 21.370413
## mvrnorm(N, mu, mcov) 19.126469 20.359832 22.074881 21.257687
## uq max neval cld
## 3.372834 5.371075 100 a
## 6.072789 8.758817 100 b
## 22.809613 26.410541 100 c
## 22.400681 78.709154 100 c
In 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) and available here. 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
## mvnfast::rmvt(N, mu, mcov, df = 3, ncores = 2) 4.361904 4.648730
## mvnfast::rmvt(N, mu, mcov, df = 3) 7.975746 8.418861
## mvtnorm::rmvt(N, delta = mu, sigma = mcov, df = 3) 23.151408 24.642707
## mean median uq max neval cld
## 5.049300 4.774214 5.028242 10.13336 100 a
## 8.849472 8.696331 8.983552 11.33661 100 b
## 27.498763 26.223508 27.553432 137.67502 100 c
When 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:
A[1:2, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 10.036161 -1.5868616 -3.54730 9.340636 12.903309
## [2,] 7.651365 0.4993065 12.10851 -1.691894 3.742399
Pre-creating the matrix of random variables saves some more time:
microwrapper(rmvn(N, mu, mcov, ncores = 2, A = A),
rmvn(N, mu, mcov, ncores = 2),
times = 200)
## Unit: milliseconds
## expr min lq mean median
## rmvn(N, mu, mcov, ncores = 2, A = A) 2.879411 3.136307 3.224661 3.238136
## rmvn(N, mu, mcov, ncores = 2) 3.129851 3.593006 4.191768 4.035298
## uq max neval cld
## 3.293530 3.896491 200 a
## 4.478869 7.830507 200 b
Don’t look at the median time here, the mean is much more affected by memory re-allocation.
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: milliseconds
## expr min lq mean median
## dmvn(X, mu, mcov, ncores = 2, log = T) 1.463944 1.771848 1.882509 1.857530
## dmvn(X, mu, mcov, log = T) 2.659815 3.143604 3.256164 3.247075
## dmvnorm(X, mu, mcov, log = T) 3.344282 4.102889 4.820034 4.218540
## uq max neval cld
## 1.946045 4.845624 500 a
## 3.350250 5.921703 500 b
## 4.433580 114.567752 500 c
Again, 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: milliseconds
## expr min lq
## mvnfast::dmvt(X, mu, mcov, df = 4, ncores = 2, log = T) 1.826992 2.023366
## mvnfast::dmvt(X, mu, mcov, df = 4, log = T) 3.183175 3.406718
## mvtnorm::dmvt(X, delta = mu, sigma = mcov, df = 4, log = T) 4.069030 4.340188
## mean median uq max neval cld
## 2.348839 2.108286 2.206562 106.677059 500 a
## 3.509517 3.509490 3.598573 4.435387 500 b
## 5.131188 4.471048 4.662687 110.718618 500 c
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: milliseconds
## expr min lq mean median uq
## maha(X, mu, mcov, ncores = 2) 1.490052 1.619832 1.759116 1.700614 1.842814
## maha(X, mu, mcov) 2.881382 3.034995 3.243939 3.153680 3.296028
## mahalanobis(X, mu, mcov) 5.567054 5.927210 6.804350 6.126685 7.900315
## max neval cld
## 2.512757 100 a
## 6.441206 100 b
## 9.882802 100 c
The acceleration is similar to that obtained in the previous sections.
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.
Dirk Eddelbuettel and Romain Francois (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. URL http://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 http://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.