r - mvtnorm::pmvnorm 的 Rcpp 实现比原始 R 函数慢
问题描述
我试图让 pmvnorm 的 Rcpp 版本至少与 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
所以在我看来,前面的代码中一定有什么事情发生了。要么是我如何处理犰狳的事情,要么是其他事情是如何联系起来的。我假设与最后一个实现相比应该有性能提升。
解决方案
我不会尝试为此使用其他库,而是尝试使用由导出的 C API mvtnorm
,参见https://github.com/cran/mvtnorm/blob/master/inst/NEWS#L44-L48。这样做时,我发现了结果不同的三个原因。其中之一还负责性能差异:
mvtnorm
使用 R 的 RNG,虽然这已从您正在使用的库中删除,请参见https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/randomF77.c。你的
triangl
函数不正确。它以列优先顺序返回下三角矩阵。但是,底层的 fortran 代码期望它以行为主,参见https://github.com/cran/mvtnorm/blob/master/src/mvt.f#L36-L39和https://github.com/zhanxw /libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L60libMvtnorm
使用1e-6
而不是1e-3
相对精度,参见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
文件中。
推荐阅读
- php - 发送一封电子邮件而不是两封收据?
- aop - Specman/e 列表列表(多维数组)
- tensorflow - 为 Go 构建 TensorFlow 时缺少 framework_go_proto
- ruby - 如何查看“更新!”的源代码 方法
- c++ - 导出从另一个模块导入的类
- r - 在 R 中创建一个 24 小时向量,时间间隔为 60 分钟和 1 分钟
- c++ - 编辑控件无法正常工作
- android - Android Studio 编译 Gradle 错误
- algorithm - 如何将值插入此二项式堆?
- python - 将 methodcaller 和 attrgetter 与 sorted 相结合