如何使用 Rcpp 将 R 函数转换为 C++ 函数？
How to convert a R function into a C++ function using Rcpp?
我定义了以下函数：
pij = function(vec){
out = vec %*% t(vec)
diag(out) = NA
out = sum(out, na.rm = T)
return(out)
}
其中 vec 是向量，例如 vec = rnorm(10^4,0,1)。
我想知道如何使用 Rcpp 包在 C++ 中编写此函数。
这是我找到的解决方案：
library(Rcpp)
library(inline)
rcpp_inc = "using namespace Rcpp;
using namespace arma;"
src = "
vec vec1 = as<vec>(vecin);
mat out = vec1*trans(vec1);
out.diag().zeros();
return(wrap(accu(out)));
"
pij_rcpp = cxxfunction(signature(vecin="numeric"), src, plugin='RcppArmadillo', rcpp_inc)
但是，它比用R写的函数慢。例如，如果我运行这个例子，
set.seed(1)
x = runif(1e4)
system.time({pij_r(x)})
system.time({pij_rcpp(x)})
我知道 pij_r 的经过时间是 1.101，pij_rcpp 的经过时间是 1.323。
这是一个更好、更直接的版本，其中 C++ 最终获胜：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
double pij_cpp(const arma::vec & v) {
arma::mat m = v * v.t();
m.diag().zeros();
double s = arma::as_scalar(arma::accu(m));
return(s);
}
/*** R
library(rbenchmark)
set.seed(123)
pij <- function(vec){
out <- vec %*% t(vec)
diag(out) <- NA
out <- sum(out, na.rm = T)
}
x <- rnorm(1000)
## make sure they are the same
all.equal(pij(x), pij_cpp(x))
## benchmark
benchmark(R=pij(x), Cpp=pij_cpp(x))
*/
在我的机器上，C++ 领先：
R> sourceCpp("~/git/so-r/53105055/answer.cpp")
R> library(rbenchmark)
R> set.seed(123)
R> pij <- function(vec){
+ out <- vec %*% t(vec)
+ diag(out) <- NA
+ out <- sum(out, na.rm = T)
+ }
R> x <- rnorm(1000)
R> ## make sure they are the same
R> all.equal(pij(x), pij_cpp(x))
[1] TRUE
R> ## benchmark
R> benchmark(R=pij(x), Cpp=pij_cpp(x))
test replications elapsed relative user.self sys.self user.child sys.child
2 Cpp 100 0.127 1.000 0.283 0.356 0 0
1 R 100 0.583 4.591 2.607 4.011 0 0
R>
更大的收获是……您看错了问题。你的 R 函数是
已经高度矢量化并调出大部分已编译代码，因此没有太多收获。
我建议首先考虑问题背后的数学原理。对于向量 v，您正在尝试计算
sum_{i=1}^{N-1} sum_{j=i+1}^{N} 2 * v_i * v_j
您可以先创建矩阵 v_i * v_j 来实现，但如果 v 很大，这样做的代价可能很大。所以直接用C++实现双和比较简单：
#include <Rcpp.h>
// [[Rcpp::export]]
double pij_cpp(Rcpp::NumericVector vec) {
double out{0.0};
int N = vec.size();
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
out += 2 * vec[i] * vec[j];
}
}
return out;
}
不过，上面的公式其实可以重新排列：
2 * sum_{i=1}^{N-1} v_i * sum_{j=i+1}^{N} v_j
这允许我们通过从高端开始到低端来摆脱双重循环:
#include <Rcpp.h>
// [[Rcpp::export]]
double pij_opt(Rcpp::NumericVector vec) {
double out{0.0};
double sum{0.0};
int N = vec.size();
for (int i = N -1; i > 0; --i) {
sum += vec[i];
out += sum * vec[i-1];
}
return 2 * out;
}
我们可以将这些版本与您的 R 代码和基于 Armadillo 的版本进行比较，以获得长度为 10^4:
的向量
> bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_arma(vec))
# A tibble: 4 x 14
expression min mean median max `itr/sec` mem_alloc n_gc n_itr total_time result
<chr> <bch:t> <bch:t> <bch:t> <bch:t> <dbl> <bch:byt> <dbl> <int> <bch:tm> <list>
1 pij(vec) 716.4ms 716.4ms 716.4ms 716.4ms 1.40 1.49GB 1 1 716ms <dbl …
2 pij_cpp(v… 59.9ms 61.4ms 61.5ms 62.3ms 16.3 2.49KB 0 9 552ms <dbl …
3 pij_opt(v… 14.2µs 15.6µs 14.9µs 864.5µs 64072. 2.49KB 0 10000 156ms <dbl …
4 pij_arma(… 834.5ms 834.5ms 834.5ms 834.5ms 1.20 2.49KB 0 1 834ms <dbl …
# ... with 3 more variables: memory <list>, time <list>, gc <list>
R 和 Armadillo 差不多（并且可能受内存分配限制）。第一个 C++ 版本快了 10 倍，第二个版本快了 50000 倍！
完整代码：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
double pij_arma(arma::vec vec) {
arma::mat out = vec * vec.t();
out.diag().zeros();
return arma::accu(out);
}
// [[Rcpp::export]]
double pij_cpp(Rcpp::NumericVector vec) {
double out{0.0};
int N = vec.size();
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
out += 2 * vec[i] * vec[j];
}
}
return out;
}
// [[Rcpp::export]]
double pij_opt(Rcpp::NumericVector vec) {
double out{0.0};
double sum{0.0};
int N = vec.size();
for (int i = N -1; i > 0; --i) {
sum += vec[i];
out += sum * vec[i-1];
}
return 2 * out;
}
/*** R
pij = function(vec){
out = vec %*% t(vec)
diag(out) = NA
out = sum(out, na.rm = T)
return(out)
}
set.seed(42)
vec = rnorm(10^4,0,1)
pij(vec)
bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_arma(vec))
*/
为了完整性：这实际上是一个算法问题，所以即使 for 循环在 R 中也比 pij_cpp:
快
pij_opt_r <- function(vec) {
out <- 0
sum <- 0
N <- length(vec)
for (i in seq.int(from = N, to = 2, by = -1)) {
sum <- sum + vec[i]
out <- out + sum * vec[i-1]
}
2 * out
}
在 R 中使用矢量化函数甚至更快，但仍然不如 pij_opt:
pij_opt_r2 <- function(vec) {
N <- length(vec)
vec <- rev(vec)
sums <- cumsum(vec)
2 * sum(vec[2:N] * sums[1:N-1])
}
完整基准：
> bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_opt_r(vec), pij_opt_r2(vec), pij_arma(vec))
# A tibble: 6 x 14
expression min mean median max `itr/sec` mem_alloc n_gc n_itr total_time result
<chr> <bch:t> <bch:tm> <bch:tm> <bch:t> <dbl> <bch:byt> <dbl> <int> <bch:tm> <list>
1 pij(vec) 733.6ms 733.6ms 733.6ms 733.6ms 1.36 1.49GB 1 1 734ms <dbl …
2 pij_cpp(v… 60ms 61.41ms 60.84ms 64.2ms 16.3 2.49KB 0 9 553ms <dbl …
3 pij_opt(v… 14.2µs 15.83µs 15.35µs 750.1µs 63164. 2.49KB 0 10000 158ms <dbl …
4 pij_opt_r… 981.1µs 1.04ms 1.02ms 1.5ms 960. 119.2KB 0 480 500ms <dbl …
5 pij_opt_r… 157µs 272.95µs 241.57µs 66.3ms 3664. 547.28KB 1 1832 500ms <dbl …
6 pij_arma(… 878.4ms 878.38ms 878.38ms 878.4ms 1.14 2.49KB 0 1 878ms <dbl …
# ... with 3 more variables: memory <list>, time <list>, gc <list>
我定义了以下函数：
pij = function(vec){
out = vec %*% t(vec)
diag(out) = NA
out = sum(out, na.rm = T)
return(out)
}
其中 vec 是向量，例如 vec = rnorm(10^4,0,1)。
我想知道如何使用 Rcpp 包在 C++ 中编写此函数。
这是我找到的解决方案：
library(Rcpp)
library(inline)
rcpp_inc = "using namespace Rcpp;
using namespace arma;"
src = "
vec vec1 = as<vec>(vecin);
mat out = vec1*trans(vec1);
out.diag().zeros();
return(wrap(accu(out)));
"
pij_rcpp = cxxfunction(signature(vecin="numeric"), src, plugin='RcppArmadillo', rcpp_inc)
但是，它比用R写的函数慢。例如，如果我运行这个例子，
set.seed(1)
x = runif(1e4)
system.time({pij_r(x)})
system.time({pij_rcpp(x)})
我知道 pij_r 的经过时间是 1.101，pij_rcpp 的经过时间是 1.323。
这是一个更好、更直接的版本，其中 C++ 最终获胜：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
double pij_cpp(const arma::vec & v) {
arma::mat m = v * v.t();
m.diag().zeros();
double s = arma::as_scalar(arma::accu(m));
return(s);
}
/*** R
library(rbenchmark)
set.seed(123)
pij <- function(vec){
out <- vec %*% t(vec)
diag(out) <- NA
out <- sum(out, na.rm = T)
}
x <- rnorm(1000)
## make sure they are the same
all.equal(pij(x), pij_cpp(x))
## benchmark
benchmark(R=pij(x), Cpp=pij_cpp(x))
*/
在我的机器上，C++ 领先：
R> sourceCpp("~/git/so-r/53105055/answer.cpp")
R> library(rbenchmark)
R> set.seed(123)
R> pij <- function(vec){
+ out <- vec %*% t(vec)
+ diag(out) <- NA
+ out <- sum(out, na.rm = T)
+ }
R> x <- rnorm(1000)
R> ## make sure they are the same
R> all.equal(pij(x), pij_cpp(x))
[1] TRUE
R> ## benchmark
R> benchmark(R=pij(x), Cpp=pij_cpp(x))
test replications elapsed relative user.self sys.self user.child sys.child
2 Cpp 100 0.127 1.000 0.283 0.356 0 0
1 R 100 0.583 4.591 2.607 4.011 0 0
R>
更大的收获是……您看错了问题。你的 R 函数是 已经高度矢量化并调出大部分已编译代码，因此没有太多收获。
我建议首先考虑问题背后的数学原理。对于向量 v，您正在尝试计算
sum_{i=1}^{N-1} sum_{j=i+1}^{N} 2 * v_i * v_j
您可以先创建矩阵 v_i * v_j 来实现，但如果 v 很大，这样做的代价可能很大。所以直接用C++实现双和比较简单：
#include <Rcpp.h>
// [[Rcpp::export]]
double pij_cpp(Rcpp::NumericVector vec) {
double out{0.0};
int N = vec.size();
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
out += 2 * vec[i] * vec[j];
}
}
return out;
}
不过，上面的公式其实可以重新排列：
2 * sum_{i=1}^{N-1} v_i * sum_{j=i+1}^{N} v_j
这允许我们通过从高端开始到低端来摆脱双重循环:
#include <Rcpp.h>
// [[Rcpp::export]]
double pij_opt(Rcpp::NumericVector vec) {
double out{0.0};
double sum{0.0};
int N = vec.size();
for (int i = N -1; i > 0; --i) {
sum += vec[i];
out += sum * vec[i-1];
}
return 2 * out;
}
我们可以将这些版本与您的 R 代码和基于 Armadillo 的版本进行比较，以获得长度为 10^4:
> bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_arma(vec))
# A tibble: 4 x 14
expression min mean median max `itr/sec` mem_alloc n_gc n_itr total_time result
<chr> <bch:t> <bch:t> <bch:t> <bch:t> <dbl> <bch:byt> <dbl> <int> <bch:tm> <list>
1 pij(vec) 716.4ms 716.4ms 716.4ms 716.4ms 1.40 1.49GB 1 1 716ms <dbl …
2 pij_cpp(v… 59.9ms 61.4ms 61.5ms 62.3ms 16.3 2.49KB 0 9 552ms <dbl …
3 pij_opt(v… 14.2µs 15.6µs 14.9µs 864.5µs 64072. 2.49KB 0 10000 156ms <dbl …
4 pij_arma(… 834.5ms 834.5ms 834.5ms 834.5ms 1.20 2.49KB 0 1 834ms <dbl …
# ... with 3 more variables: memory <list>, time <list>, gc <list>
R 和 Armadillo 差不多（并且可能受内存分配限制）。第一个 C++ 版本快了 10 倍，第二个版本快了 50000 倍！
完整代码：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
double pij_arma(arma::vec vec) {
arma::mat out = vec * vec.t();
out.diag().zeros();
return arma::accu(out);
}
// [[Rcpp::export]]
double pij_cpp(Rcpp::NumericVector vec) {
double out{0.0};
int N = vec.size();
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
out += 2 * vec[i] * vec[j];
}
}
return out;
}
// [[Rcpp::export]]
double pij_opt(Rcpp::NumericVector vec) {
double out{0.0};
double sum{0.0};
int N = vec.size();
for (int i = N -1; i > 0; --i) {
sum += vec[i];
out += sum * vec[i-1];
}
return 2 * out;
}
/*** R
pij = function(vec){
out = vec %*% t(vec)
diag(out) = NA
out = sum(out, na.rm = T)
return(out)
}
set.seed(42)
vec = rnorm(10^4,0,1)
pij(vec)
bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_arma(vec))
*/
为了完整性：这实际上是一个算法问题，所以即使 for 循环在 R 中也比 pij_cpp:
pij_opt_r <- function(vec) {
out <- 0
sum <- 0
N <- length(vec)
for (i in seq.int(from = N, to = 2, by = -1)) {
sum <- sum + vec[i]
out <- out + sum * vec[i-1]
}
2 * out
}
在 R 中使用矢量化函数甚至更快，但仍然不如 pij_opt:
pij_opt_r2 <- function(vec) {
N <- length(vec)
vec <- rev(vec)
sums <- cumsum(vec)
2 * sum(vec[2:N] * sums[1:N-1])
}
完整基准：
> bench::mark(pij(vec), pij_cpp(vec), pij_opt(vec), pij_opt_r(vec), pij_opt_r2(vec), pij_arma(vec))
# A tibble: 6 x 14
expression min mean median max `itr/sec` mem_alloc n_gc n_itr total_time result
<chr> <bch:t> <bch:tm> <bch:tm> <bch:t> <dbl> <bch:byt> <dbl> <int> <bch:tm> <list>
1 pij(vec) 733.6ms 733.6ms 733.6ms 733.6ms 1.36 1.49GB 1 1 734ms <dbl …
2 pij_cpp(v… 60ms 61.41ms 60.84ms 64.2ms 16.3 2.49KB 0 9 553ms <dbl …
3 pij_opt(v… 14.2µs 15.83µs 15.35µs 750.1µs 63164. 2.49KB 0 10000 158ms <dbl …
4 pij_opt_r… 981.1µs 1.04ms 1.02ms 1.5ms 960. 119.2KB 0 480 500ms <dbl …
5 pij_opt_r… 157µs 272.95µs 241.57µs 66.3ms 3664. 547.28KB 1 1832 500ms <dbl …
6 pij_arma(… 878.4ms 878.38ms 878.38ms 878.4ms 1.14 2.49KB 0 1 878ms <dbl …
# ... with 3 more variables: memory <list>, time <list>, gc <list>