title: Parallelizing a non-trivial Gibbs Sampler in R: RcppThread vs. RcppParallel
body: <p><strong>Overview</strong>:
  I'm interested in parallelizing (across chains) a Gibbs sampler for a non-trivial regression problem that I've <a href="https://github.com/apeterson91/rstapDP/blob/master/src/FDPP/FDPPSampler.inl" rel="nofollow noreferrer">already implemented</a> in serial via Rcpp/RcppEigen. I've read the documentation for <a href="https://rcppcore.github.io/RcppParallel/index.html" rel="nofollow noreferrer"><code>RcppParallel</code></a> and <a href="https://github.com/tnagler/RcppThread" rel="nofollow noreferrer"><code>RcppThread</code></a> and I want to know if my understanding of the challenges involved in parallelizing this code are accurate and if my proposed pseudocode using <code>RcppThread</code> is viable.</p>
  <p><strong>Programming Challenge</strong>:
  This regression problem requires inverting an updated design matrix each iteration of the Gibbs sampler. Consequently any new matrix (one per chain) needs to be &quot;thread safe&quot;. That is, there is no danger of one thread writing to memory that another thread might also try to access. If this is done, I can then draw and store the regression coefficient samples (beta) by giving the <code>Rcpp::parallelFor</code> a unique index with which to assign the samples. I'm wondering <em>where/how would be best to initialize these thread specific matrices</em>?. See below for my overall conceptual understanding and first guess at how I could essentially use the sample principle of assigning samples in parallel, to assign X's in parallel. <em>Note</em> This is assuming that Eigen objects are okay with concurrent index access in the same way I've seen std::vector&lt;&gt;'s memory accessed in the <code>RcppThread</code> documentation.</p>
  <pre><code>#include &quot;RcppEigen.h&gt;
  // [[Rcpp::plugins(cpp11)]]
  // [[Rcpp::depends(RcppThread)]] 
  // [[Rcpp::depends(RcppEigen)]] 
  
  // Sampler class definition
  #include &quot;Sampler.h&quot; 
  #include &quot;RcppThread.h&quot;
  
  // [[Rcpp::export]]
  Eigen::ArrayXXd fancyregression(const Eigen::VectorXd &amp;y, // outcome vector
                                  const Eigen::MatrixXd &amp;Z, // static sub-matrix of X
                                  const int &amp;num_iterations,
                                  const int &amp;num_chains_less_one,
                                  const int &amp;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,[&amp;beta, &amp;shared_X, &amp;n,&amp;sampler, &amp;dim_X, &amp;rng](unsigned int chain){
          // chain specific iteration for loop
          for(unsigned int iter_ix = 0; iter_ix &lt; 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);
  
  }
  </code></pre>