r - 在 R : RcppThread vs. RcppParallel 中并行化一个非平凡的 Gibbs 采样器

Question

概览:我对并行化(跨链)Gibbs 采样器感兴趣，以解决我遇到的非平凡回归问题 already implemented通过 Rcpp/RcppEigen 串行。我已阅读 RcppParallel 的文档和 RcppThread我想知道我对并行化此代码所涉及的挑战的理解是否准确，以及我提出的使用 RcppThread 的伪代码是否可行。

编程挑战:此回归问题需要在 Gibbs 采样器的每次迭代中反转更新的设计矩阵。因此，任何新矩阵(每个链一个)都需要是“线程安全的”。也就是说，不存在一个线程写入另一个线程也可能尝试访问的内存的危险。如果完成此操作，我便可以通过为 Rcpp::parallelFor 指定用于分配样本的唯一索引来绘制和存储回归系数样本 (beta)。我想知道在哪里/如何最好地初始化这些特定于线程的矩阵？。请参阅下文了解我的总体概念理解，并首先猜测我如何基本上使用并行分配样本的样本原则来并行分配 X。 注意 这是假设 Eigen 对象可以并发索引访问，就像我在 RcppThread 文档中看到的访问 std::vector<> 内存的方式一样。

#include "RcppEigen.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppThread)]] 
// [[Rcpp::depends(RcppEigen)]] 

// Sampler class definition
#include "Sampler.h" 
#include "RcppThread.h"

// [[Rcpp::export]]
Eigen::ArrayXXd fancyregression(const Eigen::VectorXd &y, // outcome vector
                                const Eigen::MatrixXd &Z, // static sub-matrix of X
                                const int &num_iterations,
                                const int &num_chains_less_one,
                                const int &seed,
                                ...)
{ 
   std::mt19937 rng;
   rng(seed);
   const int dim_X = get_dim_X(Z,...);
   const int n = y.rows();
   const int num_chains = num_chains_less_one + 1;

   Eigen::ArrayXXd beta_samples;
   beta_samples.setZero(num_iterations,num_chains*dim_X);

   Eigen::MatrixXd shared_X(n,dim_X*num_chains);

   // sampler object only has read access to its arguments
   SamplerClass sampler(y,Z,...);
    
   //chain for loop
    RcppThread::parallelFor(0, num_chains_less_one,[&beta, &shared_X, &n,&sampler, &dim_X, &rng](unsigned int chain){
        // chain specific iteration for loop
        for(unsigned int iter_ix = 0; iter_ix < num_iterations ; iter_ix ++){
            X.block(0,dim_X*chain,n,dim_X) = sampler.create_X(rng);
            beta_samples(iter_ix,dim_X*chain) = sampler.get_beta_sample(X,rng); 
        }
    });

    return(beta_samples);

}

Answer 1

“在哪里/如何最好地初始化这些线程特定的矩阵？”

您正在寻找特定于线程的资源。这是一个准系统示例:

#include <Rcpp.h>
#include <RcppParallel.h>
using namespace Rcpp;
using namespace RcppParallel;

// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]

struct Test : public Worker {
  tbb::enumerable_thread_specific<bool> printonce;
  Test() : printonce(false) {}
  
  void operator()(std::size_t begin, std::size_t end) {
    tbb::enumerable_thread_specific<bool>::reference p = printonce.local();
    if(!p) { // print once per thread
      std::cout << 1;
      p= true;
    }
  }
};

// [[Rcpp::export(rng = false)]]
void test() {
  Test x{};
  parallelFor(0, 10000, x);
}

RcppParallel 在底层使用 TBB(适用于大多数操作系统)，因此您可以使用和查找 TBB 中的任何内容。

请注意，由于它是一个线程，因此必须在某处分配本地资源，因此您需要使用类/仿函数而不是 lambda。