3

概述：我对并行化（跨链）Gibbs 采样器感兴趣，以解决我已经通过 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);

}
4

1 回答 1

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。

于 2020-10-30T18:47:50.020 回答