元素矩阵乘法:R与Rcpp(如何加速这段代码?)
Elementwise matrix multiplication: R versus Rcpp (How to speed this code up?)
我是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
-第一行的每个元素与向量testvec
的i
-第元素相乘。我的基准测试告诉我,我的实现非常慢,我希望得到如何加快速度的建议。下面是我的代码:
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
解决方案执行得非常糟糕。有什么更好的方法吗?
如果你想加快你的计算速度,你必须小心一点,不要复制。这通常意味着牺牲可读性。这是一个版本,它不复制并修改矩阵X。
// [[Rcpp::export]]
NumericMatrix Rcpp_matvecprod_elwise(NumericMatrix & X, NumericVector & y){
unsigned int ncol = X.ncol();
unsigned int nrow = X.nrow();
int counter = 0;
for (unsigned int j=0; j<ncol; j++) {
for (unsigned int i=0; i<nrow; i++) {
X[counter++] *= y[i];
}
}
return X;
}
这是我在我的机器上得到的
> library(microbenchmark)
> microbenchmark(R=R_matvecprod_elwise(X, e), Arma=A_matvecprod_elwise(X, e), Rcpp=Rcpp_matvecprod_elwise(X, e))
Unit: milliseconds
expr min lq median uq max neval
R 8.262845 9.386214 10.542599 11.53498 12.77650 100
Arma 18.852685 19.872929 22.782958 26.35522 83.93213 100
Rcpp 6.391219 6.640780 6.940111 7.32773 7.72021 100
> all.equal(R_matvecprod_elwise(X, e), Rcpp_matvecprod_elwise(X, e))
[1] TRUE
对于初学者,我将Armadillo版本(接口)编写为
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
arama::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;
}
,因为您正在进行额外的输入和输出转换(尽管wrap()
是由粘合代码添加的)。const &
是名义上的(正如您通过上一个问题了解到的那样,SEXP
是一个易于复制的指针对象),但样式更好。
你没有显示你的基准测试结果,所以我不能评论矩阵大小等的影响。我怀疑你可能会得到更好的答案在rpp -devel比这里。你选择。
编辑:如果你真的想要便宜又快的东西,我会这样做:
// [[Rcpp::export]]
mat cheapHadamard(mat X, vec y) {
// should row dim of X versus length of Y here
for (unsigned int i=0; i<y.n_elem; i++) X.row(i) *= y(i);
return X;
}
不分配新的内存,因此会更快,可能会与r竞争。
测试输出:R> cheapHadamard(testmat, testvec)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 4 10 16
[3,] 9 18 27
R>
对于c++问题给出一个本质上是C的答案,我很抱歉,但是正如已经建议的那样,解决方案通常在于有效的BLAS实现。不幸的是,BLAS本身缺少Hadamard乘法,因此您必须自己实现。
这是一个纯粹的Rcpp实现,基本上调用C代码。如果你想让它适合c++, worker函数可以被模板化,但对于大多数使用R的应用程序来说,这不是一个问题。注意,这也是"就地"操作,这意味着它修改X而不复制它。
// it may be necessary on your system to uncomment one of the following
//#define restrict __restrict__ // gcc/clang
//#define restrict __restrict // MS Visual Studio
//#define restrict // remove it completely
#include <Rcpp.h>
using namespace Rcpp;
#include <cstdlib>
using std::size_t;
void hadamardMultiplyMatrixByVectorInPlace(double* restrict x,
size_t numRows, size_t numCols,
const double* restrict y)
{
if (numRows == 0 || numCols == 0) return;
for (size_t col = 0; col < numCols; ++col) {
double* restrict x_col = x + col * numRows;
for (size_t row = 0; row < numRows; ++row) {
x_col[row] *= y[row];
}
}
}
// [[Rcpp::export]]
NumericMatrix C_matvecprod_elwise_inplace(NumericMatrix& X,
const NumericVector& y)
{
// do some dimension checking here
hadamardMultiplyMatrixByVectorInPlace(X.begin(), X.nrow(), X.ncol(),
y.begin());
return X;
}
这是一个先生成副本的版本。我不太了解Rcpp,无法在不导致性能大幅下降的情况下进行本机操作。在堆栈上创建并返回一个NumericMatrix(numRows, numCols)
会导致代码运行速度慢30%。
#include <Rcpp.h>
using namespace Rcpp;
#include <cstdlib>
using std::size_t;
#include <R.h>
#include <Rdefines.h>
void hadamardMultiplyMatrixByVector(const double* restrict x,
size_t numRows, size_t numCols,
const double* restrict y,
double* restrict z)
{
if (numRows == 0 || numCols == 0) return;
for (size_t col = 0; col < numCols; ++col) {
const double* restrict x_col = x + col * numRows;
double* restrict z_col = z + col * numRows;
for (size_t row = 0; row < numRows; ++row) {
z_col[row] = x_col[row] * y[row];
}
}
}
// [[Rcpp::export]]
SEXP C_matvecprod_elwise(const NumericMatrix& X, const NumericVector& y)
{
size_t numRows = X.nrow();
size_t numCols = X.ncol();
// do some dimension checking here
SEXP Z = PROTECT(Rf_allocVector(REALSXP, (int) (numRows * numCols)));
SEXP dimsExpr = PROTECT(Rf_allocVector(INTSXP, 2));
int* dims = INTEGER(dimsExpr);
dims[0] = (int) numRows;
dims[1] = (int) numCols;
Rf_setAttrib(Z, R_DimSymbol, dimsExpr);
hadamardMultiplyMatrixByVector(X.begin(), X.nrow(), X.ncol(), y.begin(), REAL(Z));
UNPROTECT(2);
return Z;
}
如果你对restrict
的用法感到好奇,这意味着你作为程序员与编译器签订了一个协议,不同的内存位不重叠,允许编译器进行某些优化。restrict
关键字是c++ 11(和C99)的一部分,但许多编译器为早期标准添加了c++扩展。
一些需要基准测试的R代码:
require(rbenchmark)
n <- 50000
k <- 50
X <- matrix(rnorm(n*k), nrow=n)
e <- rnorm(n)
R_matvecprod_elwise <- function(mat, vec) mat*vec
all.equal(R_matvecprod_elwise(X, e), C_matvecprod_elwise(X, e))
X_dup <- X + 0
all.equal(R_matvecprod_elwise(X, e), C_matvecprod_elwise_inplace(X_dup, e))
benchmark(R_matvecprod_elwise(X, e),
C_matvecprod_elwise(X, e),
C_matvecprod_elwise_inplace(X, e),
columns = c("test", "replications", "elapsed", "relative"),
order = "relative", replications = 1000)
和结果:
test replications elapsed relative
3 C_matvecprod_elwise_inplace(X, e) 1000 3.317 1.000
2 C_matvecprod_elwise(X, e) 1000 7.174 2.163
1 R_matvecprod_elwise(X, e) 1000 10.670 3.217
最后,就地版本实际上可能更快,因为对同一矩阵的重复乘法可能会导致一些溢出混乱。
编辑:删除循环展开,因为它没有提供任何好处,而且会分散注意力。
- 在java中解决这段代码时面临循环中的问题
- 我是如何在这段代码中出现分段错误的
- 我不明白这段代码是如何对这个pythonlist()进行排序的,也不明白如何用C++中的向量来重现它
- 为什么这段代码不起作用,我该如何解决?
- 为什么这段代码给我错误? 有没有自错?
- 有人可以解释一下这段代码如何能够反转字符串
- 可能我知道为什么这段代码没有给出任何输出吗?
- 这段代码的最后一行在做什么?
- 我不知道为什么这段代码会让核心被转储?
- 我试图用这段代码找到二叉树的高度,但它一直返回 0,有人可以告诉我为什么吗?
- 有人可以向我解释一下这段代码的作用吗?
- 存储在哪个内存段(代码/数据段)类(员工)中?
- 为什么这段代码会导致无限循环?
- 任何人都可以弄清楚这段代码如何显示运行错误?它打印无限时间 -1 以及正确答案
- 如何为一段代码启用 -permissive
- 谁能告诉我为什么这段代码没有产生正确的输出?
- 我想反转我的阵列.为什么这段代码给出垃圾值?
- 我怎样才能加速这段代码(MWE!),例如使用限制
- 如何加速这段c++代码(特别是阶乘和幂部分)
- 元素矩阵乘法:R与Rcpp(如何加速这段代码?)