问题描述
我已经用R和Rcpp编写了一个函数,该函数基本上只接受数据集x,比例参数gamma和矢量参数beta,因此返回拟合的概率。下面是我在R中的代码:
tp_predict<-function(x,gamma,beta){
n<-nrows(x)
x<-as.matrix(cbind(1,x))
vec<-x%*%beta
In<-vec<0
Ip<-!vec<0
vecN<-vec[In]
vecP<-vec[Ip]
fitted<-vector(,n)
t<-tanh(gamma)
b<-t+1
a<-1-t
fitted[In]<-b/(exp(-vecN/b)+1)
fitted[Ip]<-a/(exp(-vecP/a)+1)+t
fitted<-as.vector(fitted)
return(fitted)
}
在Rcpp中:
arma::vec tp_predict_c(arma::mat x,double gamma,arma::vec beta){
int n = x.n_rows;
arma::vec one(n,fill::ones);
arma::mat xx=join_rows(one,x);
arma::vec v = xx * beta;
arma::uvec In = find (v<0);
arma::uvec Ip = find (v>=0);
arma::vec vecN = v.elem(In);
arma::vec vecP = v.elem(Ip);
double t = tanh(gamma);
double a = 1-t;
double b = 1+t;
arma::vec fitted(n);
fitted.elem(In)=b/(exp(-vecN/b)+1);
fitted.elem(Ip)=a/(exp(-vecP/a)+1)+t;
return fitted;
}
使用rmvnorm()从多重正态分布生成测试数据x。在低尺寸(10列)和小数据量(1000行)中,Rcpp代码比R代码更好地工作
Unit: microseconds
expr min lq mean median uq max neval
tp_predict(x = X,0.3,Beta) 62.1 63.85 173.640 65.5 66.85 10868.9 100
tp_predict_c(x = X,Beta) 30.7 31.15 40.714 31.9 32.75 781.6 100
但是如果我将大小增加到100000(尺寸仍为10),则Rcpp代码会比R慢:
Unit: milliseconds
expr min lq mean median uq max neval
tp_predict(x = X,Beta) 5.8967 6.15060 8.477195 6.21560 6.7766 131.9571 100
tp_predict_c(x = X,Beta) 8.6547 8.82975 9.102348 9.01935 9.2005 12.6010 100
如果我将尺寸增加到100,尺寸仍然保持1000,Rcpp函数的效果就不如低尺寸:
Unit: microseconds
expr min lq mean median uq max neval
tp_predict(x = X,Beta) 188.7 197.0 204.745 201.85 206.5 280.5 100
tp_predict_c(x = X,Beta) 166.1 171.1 189.729 181.20 189.9 511.0 100
当尺寸和尺寸增加时,Rcpp函数的cpu时间缩放似乎比R更快。我是Rcpp的新手,我只是将R代码直接翻译为Rcpp。我不确定我是否做了一些愚蠢的事情,也不知道为什么会这样。有什么我可以解决的吗?
解决方法
您的代码正在将给定的向量和矩阵在入口处复制到tp_predict_c()
函数,这效率低下并且不是必需的。
默认情况下,C ++会按“按值”传递对象,这意味着“复制此对象”。您需要通过“按引用”传递每个对象来告诉编译器避免复制。
要通过引用传递,只需将函数的定义更改为:
arma::vec tp_predict_c(const arma::mat& x,double gamma,const arma::vec& beta)
请注意const
和&
中的const arma::mat&
关键字和const arma::vec&
字符。
const
告诉编译器该对象不会被修改,而&
告诉编译器通过引用传递。两者都允许进行额外的优化。