使用 RcppEigen 将随机状态传递给 setRandom答案

【问题标题】：Pass random state to setRandom with RcppEigen使用 RcppEigen 将随机状态传递给 setRandom
【发布时间】：2014-05-23 13:18:37
【问题描述】：

有没有办法使用 RcppEigen 将随机状态传递给 Eigen 的 setRandom 还是我需要使用 runif？

这是一个例子：

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using  Eigen::MatrixXd;
using  Eigen::VectorXd;

// [[Rcpp::export]]
NumericVector fx() {
  RNGScope   scope;
  MatrixXd   x(3,2);
  x=x.setRandom();
  x.col(1)=as<VectorXd>(runif(3,0,1));

return wrap(x);
}

测试它：

set.seed(42); fx()
#           [,1]      [,2]
#[1,] -0.8105760 0.9148060
#[2,]  0.6498853 0.9370754
#[3,]  0.6221027 0.2861395

set.seed(42); fx()
#           [,1]      [,2]
#[1,] -0.9449154 0.9148060
#[2,]  0.8063267 0.9370754
#[3,] -0.0673205 0.2861395

请注意第 2 列（即runif）如何可重现，但第 1 列（即setRandom）则不能。

【问题讨论】：

标签： r eigen rcpp

【解决方案1】：

Eigen 中的 RNG 与 R 的 RNG 正交。我们有RNGScope 处理R 并允许您set.seed() 如您所知；你用 Eigen 的 RNG 做什么是分开的。

特别是，您必须添加胶水将 Eigen 的 RNG 注册为 R 的用户提供的 RNG。默认情况下，R 不知道 Eigen，这是有道理的，因为 R 用户期望 R 的 RNG。

您的问题似乎是一个普通的“使用 Eigen 编程”问题，您无法初始化 Eigen RNG。与Rcpp无关。

编辑：这是您的示例，已更正。原来 Eigen 使用系统 RNG，所以你需要 srand() 来播种它。

R> sourceCpp("/tmp/roland.cpp")
R> set.seed(42); fx(42)
          [,1]     [,2]
[1,] -0.933060 0.914806
[2,] -0.340072 0.937075
[3,]  0.381271 0.286140
R> set.seed(42); fx(42)
          [,1]     [,2]
[1,] -0.933060 0.914806
[2,] -0.340072 0.937075
[3,]  0.381271 0.286140
R>

它使用你的代码的这个修改版本：

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using  Eigen::MatrixXd;
using  Eigen::VectorXd;

// [[Rcpp::export]]

  RNGScope   scope;
  MatrixXd   x(3,2);
  srand(seed);
  x=x.setRandom();
  x.col(1)=as<VectorXd>(runif(3,0,1));

  return wrap(x);
}

编辑 2 以响应 OP 的 cmets：

我不认为 Rcpp::runif() 和 Eigen 发出的 srand() 调用之间的速度差异很大，所以你仍然坚持srand() 一直存在问题，并且系统之间的行为可能不同.

快速演示脚本：

#include <RcppEigen.h> 

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
Rcpp::NumericVector v1(int n) {
  return Rcpp::runif(n);
}

// [[Rcpp::export]]
Rcpp::NumericVector v2(int n) {
  Eigen::VectorXd   x(n);
  x = x.setRandom();
  return Rcpp::wrap(x);
}

/*** R
library(rbenchmark)
N <- 1e7
benchmark(v1(N), v2(N))
*/

产生

R> sourceCpp("/tmp/roland.cpp")

R> library(rbenchmark)

R> N <- 1e7

R> benchmark(v1(N), v2(N))
   test replications elapsed relative user.self sys.self user.child sys.child
1 v1(N)          100  12.633    1.000    11.356    1.261          0         0
2 v2(N)          100  17.222    1.363    13.981    3.198          0         0
R>

请注意，RcppEigen 在这里较慢，即使在仅创建向量的更简单设置中也是如此。但是我们在这里讨论的是微秒，这可能不是我担心的任何一种方式，在一个很可能有其他瓶颈的真正的应用程序中。

【讨论】：

是的，我还发现我可以使用srand，但这意味着我必须为种子使用参数。我当然可以这样做，但是我必须在 R 级别处理与 set.seed 的接口，这似乎不是最佳的。
你看到我说的正交了吗？这是两个不同的RNG，其中一个实际上并不是那么好。但简而言之，当您坚持使用两个不同的 RNG 时，您还需要播种两个不同的 RNG。没有免费的午餐，等等。
我明白这一点。我还将对as<VectorXd>(runif(3,0,1)); 替代方案的速度进行基准测试。如果速度太慢，我可以使用sample.int(2^31-1, 1) 来获取我传递给srand 的整数。
你还没有解释为什么你坚持srand()。它依赖于系统，并且（至少在某些系统上）不是很好。使用 R 的 RNG，直到你有充分的理由不这样做。你没有提供这样的理由。
感谢您的提示。 setRandom 比 as<VectorXd>(runif()) 快得多（用于创建相同长度的向量）。我现在使用runif，但由于它是在循环中使用的，因此我创建了一次迭代所需的随机数的 n 倍，并在while 循环耗尽预先计算的随机数时更新它。问题是 n 很难优化。太小意味着对runif 的调用过多，太大意味着内存需求过多并且速度也会变慢。最佳值取决于用户输入，所以我设置了一个合理的默认值，并让用户可以选择调整它...