【发布时间】：2014-07-24 12:09:01
【问题描述】：
我是 C++ 编程的新手（使用 Rcpp 无缝集成到 R），我希望能得到一些关于如何加快计算速度的建议。
考虑以下示例：
testmat <- matrix(1:9, nrow=3)
testvec <- 1:3
testmat*testvec
# [,1] [,2] [,3]
#[1,] 1 4 7
#[2,] 4 10 16
#[3,] 9 18 27
在这里，R 回收了testvec，因此，粗略地说，testvec“变成”了一个与testmat 相同维度的矩阵，用于此乘法。然后返回 Hadamard 产品。我希望使用Rcpp 来实现此行为，即我希望矩阵testmat 中i-th 行的每个元素都与向量testvec 的i-th 元素相乘。我的基准测试告诉我，我的实现非常慢，我希望能就如何加快速度提出建议。这是我的代码：
首先，使用Eigen：
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
using namespace Rcpp;
using namespace Eigen;
// [[Rcpp::export]]
NumericMatrix E_matvecprod_elwise(NumericMatrix Xs, NumericVector ys){
Map<MatrixXd> X(as<Map<MatrixXd> >(Xs));
Map<VectorXd> y(as<Map<VectorXd> >(ys));
int k = X.cols();
int n = X.rows();
MatrixXd Y(n,k) ;
// here, I emulate R's recycling. I did not find an easier way of doing this. Any hint appreciated.
for(int i = 0; i < k; ++i) {
Y.col(i) = y;
}
MatrixXd out = X.cwiseProduct(Y);
return wrap(out);
}
这是我使用 Armadillo 的实现（调整以遵循 Dirk 的示例，请参阅下面的答案）：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
arma::mat A_matvecprod_elwise(const arma::mat & X, const arma::vec & y){
int k = X.n_cols ;
arma::mat Y = repmat(y, 1, k) ; //
arma::mat out = X % Y;
return out;
}
使用 R、Eigen 或 Armadillo 对这些解决方案进行基准测试表明，Eigen 和 Armadillo 都比 R 慢大约 2 倍。有没有办法加快这些计算速度或至少与 R 一样快？有没有更优雅的设置方法？任何建议表示赞赏和欢迎。 （我也鼓励对一般编程风格进行切题的评论，因为我是 Rcpp / C++ 的新手。）
这里有一些可重现的基准：
# for comparison, define R function:
R_matvecprod_elwise <- function(mat, vec) mat*vec
n <- 50000
k <- 50
X <- matrix(rnorm(n*k), nrow=n)
e <- rnorm(n)
benchmark(R_matvecprod_elwise(X, e), A2_matvecprod_elwise(X, e), E_matvecprod_elwise(X,e),
columns = c("test", "replications", "elapsed", "relative"), order = "relative", replications = 1000)
这会产生
test replications elapsed relative
1 R_matvecprod_elwise(X, e) 1000 10.89 1.000
2 A_matvecprod_elwise(X, e) 1000 26.87 2.467
3 E_matvecprod_elwise(X, e) 1000 27.73 2.546
如您所见，我的Rcpp-solutions 表现相当糟糕。有什么更好的办法吗？
【问题讨论】：
-
您可以在 Eigen 中执行
X = X.cwiseProduct(y.replicate(1, X.cols()));，然后在return Xs;中执行，但是对于 1e3*1e3 矩阵，它并没有提高我系统上的 Eigen 函数的基准。 -
在犰狳中，尝试 .each_row() 函数：arma.sourceforge.net/docs.html#each_colrow 类似以下内容：X.each_row() %= y
-
仔细查看您的代码，使用.each_col() 可能更正确：
X.each_col() %= y -
谢谢@mtall！您建议的更改大大加快了
arma-函数的速度，因此它现在超过了R，即使Rcpp 版本仍然优于arma。无论如何感谢您的评论。非常感谢您的帮助。 -
@coffeinjunky - 这是代码可读性、紧凑性和速度之间的权衡。虽然您始终可以编写针对手头特定任务进行优化的低级代码，但它几乎总是以更难以阅读和维护为代价（因此您更有可能产生错误）。通常只使用库提供的工具会更安全、更容易，例如 .each_col() 方法。
标签： c++ r eigen rcpp armadillo