mvtnorm::pmvnorm 的 Rcpp 实现比原始 R 函数慢

Question

我正在尝试让 Rcpp 版本的 pmvnorm 至少与 R 中的 mvtnorm::pmvnorm 一样快。

我找到了 https://github.com/zhanxw/libMvtnorm并使用相关源文件创建了一个 Rcpp 框架包。我添加了以下使用 Armadillo 的函数(因为我在编写的其他代码中使用了它)。

//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
  arma::mat LL = arma::trimatl(X, -1);  // omit the main diagonal
  return LL.elem(arma::find(LL != 0));
}

//[[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound, arma::vec& lowtrivec){
  double error;
  int n = bound.n_elem;
  double* boundptr = bound.memptr();
  double* lowtrivecptr = lowtrivec.memptr();
  double result = pmvnorm_P(n, boundptr, lowtrivecptr, &error);
  return result;
}

从构建包后的 R 中，这是一个可重现的示例:

set.seed(1)
covar <- rWishart(1, 10, diag(5))[,,1]
sds <- diag(covar) ^-.5
corrmat <- diag(sds) %*% covar %*% diag(sds)
triang <- triangl(corrmat)

bounds <- c(0.5, 0.9, 1, 4, -1)

rbenchmark::benchmark(pmvnorm_cpp(bounds, triang),
                      mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
                      replications=1000)

这表明 pmvnorm_cpp 比 mvtnorm::pmvnorm 慢得多。结果是不同的。

> pmvnorm_cpp(bounds, triang)
[1] 0.04300643
> mvtnorm::pmvnorm(upper=bounds, corr = corrmat)
[1] 0.04895361

这让我感到困惑，因为我认为基本的 Fortran 代码是相同的。我的代码中有什么东西让一切变慢了吗？或者我应该尝试直接移植 mvtnorm::pmvnorm 代码吗？我几乎没有使用 Fortran 的经验。

感谢建议，请原谅我的无能。

编辑:为了与替代品进行快速比较，这个:

//[[Rcpp::export]]
NumericVector pmvnorm_cpp(NumericVector bound, NumericMatrix cormat){
  Environment stats("package:mvtnorm"); 
  Function f = stats["pmvnorm"];

  NumericVector lower(bound.length(), R_NegInf);
  NumericVector mean(bound.length());
  NumericVector res = f(lower, bound, mean, cormat);
  return res;
}

具有与 R 调用基本相同的性能(以下是 40 维 mvnormal):

> rbenchmark::benchmark(pmvnorm_cpp(bounds, corrmat),
+                       mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
+                       replications=100)
                                              test replications elapsed relative user.self sys.self
2 mvtnorm::pmvnorm(upper = bounds, corr = corrmat)          100   16.86    1.032     16.60     0.00
1                     pmvnorm_cpp(bounds, corrmat)          100   16.34    1.000     16.26     0.01

所以在我看来，前面的代码中一定有什么事情发生了。无论是我如何处理 Armadillo 的事情，还是其他事情之间的联系。我认为与上一个实现相比应该有性能提升。

Answer 1

我不会为此尝试使用额外的库，而是尝试使用 mvtnorm 导出的 C API，c.f. https://github.com/cran/mvtnorm/blob/master/inst/NEWS#L44-L48 .在这样做的同时，我发现了导致结果不同的三个原因。其中之一还对性能差异负责:

1. mvtnorm 使用 R 的 RNG，而这已从您正在使用的库中删除，c.f. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/randomF77.c .

2. 您的 triangl 函数不正确。它以列优先顺序返回下三角矩阵。然而，底层的 Fortran 代码期望它按行优先顺序排列，c.f. https://github.com/cran/mvtnorm/blob/master/src/mvt.f#L36-L39和 https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L60

3. libMvtnorm 使用 1e-6 而不是 1e-3 作为相对精度，c.f. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L65 .这也是造成性能差异的原因。

我们可以使用以下代码进行测试:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::depends(mvtnorm)]]
#include <mvtnormAPI.h>

//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
  int n = X.n_cols;
  arma::vec res(n * (n-1) / 2);
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < i; ++j) {
      res(j + i * (i-1) / 2) = X(i, j);
    }
  }
  return res;
}

// [[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound,
           arma::vec& lowertrivec,
           double abseps = 1e-3){

  int n = bound.n_elem;
  int nu = 0;
  int maxpts = 25000;     // default in mvtnorm: 25000
  double releps = 0;      // default in mvtnorm: 0
  int rnd = 1;            // Get/PutRNGstate

  double* bound_ = bound.memptr();
  double* correlationMatrix = lowertrivec.memptr();
  double* lower = new double[n];
  int* infin = new int[n];
  double* delta = new double[n];

  for (int i = 0; i < n; ++i) {
    infin[i] = 0; // (-inf, bound]
    lower[i] = 0.0;
    delta[i] = 0.0;
  }

  // return values
  double error;
  double value;
  int inform;

  mvtnorm_C_mvtdst(&n, &nu, lower, bound_,
           infin, correlationMatrix, delta,
           &maxpts, &abseps, &releps,
           &error, &value, &inform, &rnd);
  delete[] (lower);
  delete[] (infin);
  delete[] (delta);

  return value;
}

/*** R
set.seed(1)
covar <- rWishart(1, 10, diag(5))[,,1]
sds <- diag(covar) ^-.5
corrmat <- diag(sds) %*% covar %*% diag(sds)
triang <- triangl(corrmat)
bounds <- c(0.5, 0.9, 1, 4, -1)
set.seed(1)
system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
 */

结果:

> system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
0.04896221 
   user  system elapsed 
  0.000   0.003   0.003 

> system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
0.04895756 
   user  system elapsed 
  0.035   0.000   0.035 

> system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
0.04896221 
   user  system elapsed 
  0.004   0.000   0.004 

在相同的 RNG(和 RNG 状态)、正确的下三角相关矩阵和相同的相对精度下，结果相同且性能相当。精度越高，性能就越差。

所有这些都是针对使用 Rcpp::sourceCpp 的独立文件。为了在包中使用它，您需要将 LinkingTo: mvtnorm 添加到您的 DESCRIPTION 文件。