Memoise 一个 Rcpp 函数？

Question

我在 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

编辑：在发布这篇文章之前，我从德克的书中得到了一个工作方法。

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 包在每个进程上优化 objective 函数，当我使用 %do% 时它可以工作，但是当我使用 %dopar% 时，我看到的只是 objective 函数不能在初始参数上进行评估。我从德克的许多其他帖子中采纳了他的建议，并将 Coatless 的方法放入一个包中，但我不确定如何将德克书中的方法放入一个包中。这只是我在 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);

}

Answer 1

记住我

嗯，先想想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)
  })

大 O 算法

免责声明：接下来应该有更正式的争论，但已经有一段时间了。

为了理解算法，有时我们需要使用所谓的 Big O notation 来观察代码如何渐近执行。现在，本例中的大 O 是 O(2^N)，因为有两个计算调用：Fun(i-1) 和 FunR(i-2)。但是，memoise 使用散列映射/查找 table，在最坏的情况下可能有 O(n) 的大 O，在最好的情况下可能有 O(1)。本质上，我们有常数与指数渐近结果。

改进微基准测试 - W/O 在 C++ 中记忆

然而，这并不一定意味着 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

C++ 中的记忆

我们可以在 C++ 中应用 memoise 中给出的相同记忆技术。实现不是那么漂亮和漂亮，但它表明相同的原则是适用的。

首先，我们将制作一个哈希图。

// 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();

现在，我们可能应该制作一些访问器来处理地图对象。

// 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 倍！记忆绝对值得！

Memoise 一个 Rcpp 函数？

Memoise an Rcpp function?

r

performance

recursion

rcpp

memoization

记住我

`memoise` 关于 R 和 C++ 的缓存访问

大 O 算法

改进微基准测试 - W/O 在 C++ 中记忆

C++ 中的记忆

Memoise 一个 Rcpp 函数？

Memoise an Rcpp function?

r

performance

recursion

rcpp

memoization

记住我

memoise 关于 R 和 C++ 的缓存访问

大 O 算法

改进微基准测试 - W/O 在 C++ 中记忆

C++ 中的记忆

`memoise` 关于 R 和 C++ 的缓存访问