Ben*_*est 4 recursion performance r memoization rcpp
我在R中编写了一个递归函数,并使用memoise来加速它.我试图通过在Rcpp中编写它然后记忆Rcpp函数来进一步加速它,但R函数更快.为什么会这样,有什么方法可以加快我的使用速度？
require(microbenchmark)
require(Rcpp)
require(memoise)
Rcpp功能:
cppFunction('
double FunCpp (unsigned int i, double d1, double d2,
double p, double s, unsigned int imax,
double n, double k, double r,
double m, double t) {
if (i == 0) return 0;
if (i == 1) return log2(-1*d1);
if (i == 2) return log2(d2*d1 - p*s);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunCpp(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunCpp(i-2, d1, d2, p, s, imax, n, k, r, m, t);
return x + log2(1 - pow(2,y-x));
}
')
FunCpp = memoise(FunCpp)
R功能:
FunR = memoise(function(i, d1, d2, p, s, imax, n, k, r, m, t) {
if(i == 0) 0
else if(i == 1) log2(-1*d1)
else if(i == 2) log2(d2*d1 - p*s)
else {
x = log2(abs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)))
x = x + FunR(i-1, d1, d2, p, s, imax, n, k, r, m, t)
y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax))
y = y + FunR(i-2, d1, d2, p, s, imax, n, k, r, m, t)
x + log2(1 - 2^(y-x))
}
})
这种速度比较对我来说是切合实际的.递归函数用于一系列整数,但之后,不会再使用相同的输入调用它.所以为了速度比较,这里我从其他函数中调用函数,在我调用递归函数之后,我使用forget()来重置缓存.
TestFunR = function() {
x = sapply(1:31, function(i) {
FunR(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
forget(FunR)
}
TestFunCpp = function() {
x = sapply(1:31, function(i) {
FunCpp(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
forget(FunCpp)
}
microbenchmark(TestFunR(), TestFunCpp())
Unit: milliseconds
expr min lq mean median uq max neval cld
TestFunR() 9.979738 10.4910 12.83228 10.91887 11.89264 61.61513 100 a
TestFunCpp() 520.955483 528.6965 547.31103 536.73058 547.66377 729.57631 100 b
编辑:在发布之前我从Dirk的书中得到了一个方法.
includeText = '
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>
class F {
public:
F(unsigned int n = 200, double d1 = 0, double d2 = 0, double p = 0, double s = 0) {
memo.resize(n);
std::fill( memo.begin(), memo.end(), NAN );
memo[0] = 0;
memo[1] = log2(-1*d1);
memo[2] = log2(d2*d1 - p*s);
}
double FunIL(int i, double d1, double d2, double p, double s, double imax,
double n, double k, double r, double m, double t) {
if (i < 0) return((double) NAN);
if (i >= (int) memo.size()) throw std::range_error(\"i too large\");
if (!std::isnan(memo[i])) return(memo[i]);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunIL(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunIL(i-2, d1, d2, p, s, imax, n, k, r, m, t);
memo[i] = x + log2(1 - pow(2,y-x));
return(memo[i]);
}
private:
std::vector< double > memo;
};
'
bodyText = '
int is = Rcpp::as<int>(i);
double d1s = Rcpp::as<double>(d1);
double d2s = Rcpp::as<double>(d2);
double ps = Rcpp::as<double>(p);
double ss = Rcpp::as<double>(s);
double imaxs = Rcpp::as<double>(imax);
double ns = Rcpp::as<double>(n);
double ks = Rcpp::as<double>(k);
double rs = Rcpp::as<double>(r);
double ms = Rcpp::as<double>(m);
double ts = Rcpp::as<double>(t);
F f(ns, d1s, d2s, ps, ss);
return Rcpp::wrap( f.FunIL(is, d1s, d2s, ps, ss, imaxs, ns, ks, rs, ms, ts) );
'
FunInline = cxxfunction(signature(i = "integer", d1 = "numeric", d2 = "numeric", p = "numeric",
s = "numeric", imax = "numeric", n = "numeric", k = "numeric",
r = "numeric", m = "numeric", t = "numeric"),
plugin = "Rcpp",
verbose = T,
incl = includeText,
body = bodyText)
它同样有效(参见TestFunInline):
microbenchmark(TestFunR(), TestFunCpp(), TestFunCpp_Mem(), TestFunInline())
Unit: microseconds
expr min lq mean median uq max neval cld
TestFunR() 8871.251 9067.758 10301.8003 9287.5725 9593.1310 19270.081 100 b
TestFunCpp() 514415.356 517160.251 522431.2980 519321.6130 523811.7640 584812.731 100 c
TestFunCpp_Mem() 245.474 264.291 284.8908 281.6105 292.0885 526.870 100 a
TestFunInline() 279.686 295.723 378.2134 306.8425 316.0370 6621.364 100 a
但是,我无法使用doParallel.我正在使用optim和optimx包优化每个进程的目标函数,当我使用%do%时,它可以工作,但是当我使用%dopar%时,我看到的是目标函数无法在初始参数上进行评估.我从他的许多其他帖子中得到了Dirk的建议,并将Coatless的方法放入了一个包中,但我不确定如何将Dirk的书中的方法放入包中.这只是我对C++的缺乏经验.
编辑2:它最终点击了如何将Dirk的方法放在我的包中的源文件中.我知道还有其他关于将Rcpp与doParallel一起使用的讨论,但是我把这个代码放在这里因为它是解决我的问题的另一个好方法,并且通过将这个代码添加到我的包中的源文件中,它碰巧变得更容易了对于我来说,我的并行方法比内联方式更好.
class F {
public:
F(unsigned int n = 200, double d1 = 0, double d2 = 0, double p = 0, double s = 0) {
memo.resize(n);
std::fill( memo.begin(), memo.end(), NAN );
memo[0] = 0;
memo[1] = log2(-1*d1);
memo[2] = log2(d2*d1 - p*s);
}
double FunIL(int i, double d1, double d2, double p, double s, double imax,
double n, double k, double r, double m, double t) {
if (i < 0) return((double) NAN);
if (i >= (int) memo.size()) throw std::range_error("\"i too large\"");
if (!std::isnan(memo[i])) return(memo[i]);
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunIL(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunIL(i-2, d1, d2, p, s, imax, n, k, r, m, t);
memo[i] = x + log2(1 - pow(2,y-x));
return(memo[i]);
}
private:
std::vector< double > memo;
};
// [[Rcpp::export]]
double FunDirk(int i, double d1, double d2, double p, double s,
double imax, double n, double k, double r,
double m, double t) {
F f(n, d1, d2, p, s);
return f.FunIL(i, d1, d2, p, s, imax, n, k, r, m, t);
}
好吧,首先让我们考虑一下这个目的memoise.的目标memoise是缓存功能的结果和重用他们.因此,在一次计算之后,它不再需要为计算中的任何其他序列再次重新计算该值,它只能从高速缓存中检索该值.这与递归结构设置特别相关.
memoise 关于R和C++的缓存访问memoisize的设置是缓存R值函数值.在这种情况下,它正在缓存这些值.但是,C++代码无法访问缓存的值.因此,C++版本重新计算每个值.从本质上讲,您实际上是在使用:
x = sapply(1:31, function(i) {
FunCpp(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
免责声明:接下来应该有一个更正式的论点,但它已经有一段时间了.
为了理解算法,有时我们需要使用所谓的Big O表示法,它允许我们观察代码如何渐近地执行.现在,在这种情况下,Big O是O(2 ^ N),因为有两个计算调用:Fun(i-1)和FunR(i-2).但是,memoise使用哈希映射/查找表,其中可能是O(n)最坏的并且O(1)处于最佳状态的大O. 本质上,我们有恒定的指数渐近结果.
但是,这并不一定意味着C++函数是垃圾.R到Rcpp和后桥的缺点之一是在两个域之间传输值之间的滞后时间.因此,我们可以稍微降低计算时间的一种方法是将循环完全放在C++中.
例如
// [[Rcpp::export]]
Rcpp::NumericVector FunCpp_loop(unsigned int e,
double d1, double d2,
double p, double s, unsigned int imax,
double n, double k, double r,
double m, double t){
Rcpp::NumericVector o(e);
for(unsigned int i = 0; i < e; i++){
o(i) = FunCpp(31-(i+1), -152, -147.33, 150, 0.03, 30, 31, 1, 1, 2, 5);
}
return o;
}
但是,这里的基准并没有真正改善这种情况(即使预先创建了矢量1:31)
Unit: milliseconds
expr min lq mean median uq max neval
TestFunR(tv) 8.467568 9.077262 9.986837 9.449952 10.60555 14.91243 100
TestFunCpp(tv) 476.678391 482.489094 487.687811 486.351087 490.25346 579.38161 100
TestFunCpp_loop() 478.348070 482.588307 488.234200 486.211347 492.33965 521.10918 100
我们可以应用memoiseC++中给出的相同的memoziation技术.实现不是那么漂亮和漂亮,但它有助于显示相同的原则是适用的.
首先,我们将制作一个hashmap.
// Memoization structure to hold the hash map
struct mem_map{
// Initializer to create the static (presistent) map
static std::map<int, double> create_map()
{
std::map<int, double> m;
m.clear();
return m;
}
// Name of the static map for the class
static std::map<int, double> memo;
};
// Actuall instantiate the class in the global scope (I know, bad me...)
std::map<int, double> mem_map::memo = mem_map::create_map();
现在,我们可能应该使用一些访问器来处理map对象.
// Reset the map
// [[Rcpp::export]]
void clear_mem(){
mem_map::memo.clear();
}
// Get the values of the map.
// [[Rcpp::export]]
std::map<int, double> get_mem(){
return mem_map::memo;
}
最后,让我们改变你函数中的一些内部事物.
// Users function
// [[Rcpp::export]]
double FunCpp_Mem (int i, double d1, double d2,
double p, double s, unsigned int imax,
double n, double k, double r,
double m, double t) {
// We have already computed the value...
if(mem_map::memo.count(i) > 0)
return mem_map::memo[i];
// Otherwise, let us get ready to compute it!
double res = 0;
if (i <= 2){
if (i <= 0) { // i == 1
res = 0.0;
}else if (i == 1) {
res = log2(-1.0*d1);
}else { // i == 2
res = log2(d2*d1 - p*s);
}
// Store result in hashmap
mem_map::memo[i] = res;
return res;
}
// Calculate if not in special case.
double x = log2(fabs(-(((imax - (n - i))/imax)*k*r + m + (n - i)*t)));
x = x + FunCpp_Mem(i-1, d1, d2, p, s, imax, n, k, r, m, t);
double y = log2((n - i + 1)*t*k*r*((imax - ((n - i + 1) - 1))/imax));
y = y + FunCpp_Mem(i-2, d1, d2, p, s, imax, n, k, r, m, t);
res = x + log2(1 - pow(2,y-x));
// Update the hashmap for uncalculated value
mem_map::memo[i] = res;
return res;
}
呼大量的工作.让我们测试一下,看看它是否值得.
# Benchmark for Rcpp Memoization
TestFunCpp_mem = function(tv) {
x = sapply(tv, function(i) {
FunCpp_Mem(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
clear_mem()
}
TestFunR = function(tv) {
x = sapply(tv, function(i) {
FunR(i = 31-i, d1 = -152, d2 = -147.33, p = 150, s = 0.03,
imax = 30, n = 31, k = 1, r = 1, m = 2, t = 5)
})
forget(FunR)
}
# Pre-generate vector
tv = 1:31
microbenchmark(TestFunR(tv),TestFunCpp_mem(tv))
结果......
microbenchmark(TestFunR(tv),TestFunCpp_mem(tv))
Unit: microseconds
expr min lq mean median uq max neval
TestFunR(tv) 8246.324 8662.694 9345.6947 9009.868 9797.126 13001.995 100
TestFunCpp_mem(tv) 203.832 214.939 253.7931 228.898 240.906 1277.325 100
具有记忆功能的Cpp功能比R版本快40.5倍!记忆绝对是值得的!