Блог

Я относительно новичок в Rcpp, недавно научившись переносить свой код из R. Однако это все еще недостаточно быстро, и я пытаюсь реализовать RcppParallel без особого успеха. По сути, у меня есть сэмплер Гиббса, обновляющий набор переменных через Метрополис-Гастингс, поскольку большинство моих переменных не имеют апостериоров закрытой формы. Наконец-то я смог запустить свой код без ошибок, но, к сожалению, мой код, использующий RcppParallel, работает с той же скоростью, что и мой старый код, и я понятия не имею, почему.

Чтобы упростить ситуацию, предположим, я хочу обновить переменную U, которая представляет собой матрицу p на r. У меня есть функция logFC_U для вычисления задней части журнала, а затем функция MH_U для выполнения шага Metropolis-Hastings и возврата обновленного столбца U.

 #include lt;RcppDist.hgt; #include lt;RcppArmadillo.hgt; #include lt;RcppParallel.hgt;  #include lt;cmathgt; #include lt;randomgt; #include lt;stringgt; #include lt;vectorgt;  using namespace Rcpp; using namespace arma; using namespace RcppParallel;  // [[Rcpp::depends(RcppArmadillo, RcppDist, RcppParallel)]]  double logFC_U(arma::vec amp;var_val, std::size_t j,   const arma::mat amp;Y, const arma::mat amp;X,   const arma::mat amp;eta, const arma::mat amp;lambda,  const arma::mat amp;U0, const arma::mat amp;V,   const arma::vec amp;delta, double tau){  double lp = 0;   mat U = U0;  U.col(j) = var_val;  mat B = U*diagmat(delta)*V.t();  double SS = arma::as_scalar(sum(sum(log(1 exp(X*B eta*lambda.t())))));  lp  = arma::as_scalar(-0.5*tau*U.col(j).t()*U.col(j)     delta(j)*U.col(j).t()*X.t()*Y*V.col(j)-SS);  return lp; }  arma::vec MH_U(std::size_t i, std::size_t j,   const arma::mat amp;Y, const arma::mat amp;X,   const arma::mat amp;eta, const arma::mat amp;lambda,  const arma::mat amp;U, const arma::mat amp;V,   const arma::vec amp;delta, double tau, double eps) {  arma::mat U_new = U;  arma::vec var_current, var_new;  arma::mat B = U*diagmat(delta)*V.t();    int p = U.n_rows;  var_current = U.col(j);   var_new = U.col(j);  var_new(i) = rnorm(1, var_current(i), eps)(0);  U_new.col(j) = var_new;    arma::vec res = vec(p 1, arma::fill::zeros); // first element will store acceptance int    // get log posterior for var_current, var_new  double l0 = logFC_U(var_current, j, Y, X, eta, lambda, U, V, delta, tau);  double l1 = logFC_U(var_new, j, Y, X, eta, lambda, U_new, V, delta, tau);    // acceptance probability  double prob = std::min(1.0, arma::as_scalar(exp(l1-l0)));    // sample binomial(1,prob)  res(0) = R::rbinom(1, prob);    // return var_current if acc=0 and var_new if acc=1  res.subvec(1, p) = res(0)*var_new   (1-res(0))*var_current;  return res; }  

Затем у меня есть функция, которая перебирает все элементы U, обновляя по одному элементу за раз и возвращая другую матрицу, которая в основном представляет собой две матрицы, расположенные рядом друг с другом: первая показывает, какие элементы были обновлены или не использовались 1 и 0, а следующая-обновленная матрица U.

 // [[Rcpp::export]] arma::mat RegUpdateU(const arma::mat amp;Y, const arma::mat amp;X,   const arma::mat amp;eta, const arma::mat amp;lambda,  const arma::mat amp;U, const arma::mat amp;V,   const arma::vec amp;delta, double tau, double eps){    arma::mat res(U.n_rows, 2*U.n_cols);    for(int i=0; ilt;U.n_rows; i  ){  for(int j=0; jlt;U.n_cols; j  ){  arma::vec prop = MH_U(i, j, Y, X, eta, lambda, U, V, delta, tau, eps);  res(i,j) = prop(0);  res(i,j U.n_cols) = prop(i 1);  }  }    return(res); }  

Это шаг, который я попытался распараллелить. Поскольку мне нужно использовать RcppArmadillo (я думаю), чтобы выполнить все умножение матриц, я использовал метод, который нашел в другом посте, для преобразования входных данных в arma::mat и arma::vec. Я также пытался сохранить входные данные как NumericMatrix и NumericVector, а затем преобразовать их в arma::mat и arma::vec внутри функции MH_U, но получал сообщение об ошибке, что не было соответствующей функции для вызова «MH_U», и я не смог это решить.

 struct UpdateU : public Worker {    // inputs  const double eps;  const double tau;  const RMatrixlt;doublegt; Y;  const RMatrixlt;doublegt; X;  const RMatrixlt;doublegt; Eta;  const RMatrixlt;doublegt; Lambda;  const RMatrixlt;doublegt; U;  const RMatrixlt;doublegt; V;  const RVectorlt;doublegt; Delta;  std::string proposal;  std::size_t n, d, p, r, q;    // output matrix to write to  RMatrixlt;doublegt; res;    UpdateU(const NumericMatrix Y, const NumericMatrix X,  const NumericMatrix Eta, const NumericMatrix Lambda,  const NumericMatrix U, const NumericMatrix V,  const NumericVector Delta, const double tau, const double eps,  std::size_t n, std::size_t d, std::size_t p,   std::size_t q, std::size_t r, NumericMatrix res)  : Y(Y), X(X), Eta(Eta), Lambda(Lambda), U(U), V(V), Delta(Delta),   tau(tau), eps(eps), n(n), d(d), p(p), q(q), r(r), res(res) {}    // convert inputs to arma::mat  arma::mat convertY(){  RMatrixlt;doublegt; y = Y;  arma::mat MAT(y.begin(), n, d, false);  return MAT;  }  arma::mat convertX(){  RMatrixlt;doublegt; x = X;  arma::mat MAT(x.begin(), n, p, false);  return MAT;  }  arma::mat convertEta(){  RMatrixlt;doublegt; eta = Eta;  arma::mat MAT(eta.begin(), n, q, false);  return MAT;  }  arma::mat convertLambda(){  RMatrixlt;doublegt; lambda = Lambda;  arma::mat MAT(lambda.begin(), d, q, false);  return MAT;  }  arma::mat convertU(){  RMatrixlt;doublegt; u = U;  arma::mat MAT(u.begin(), p, r, false);  return MAT;  }  arma::mat convertV(){  RMatrixlt;doublegt; v = V;  arma::mat MAT(v.begin(), d, r, false);  return MAT;  }  arma::vec convertDelta(){  RVectorlt;doublegt; delta = Delta;  arma::vec VEC(delta.begin(), r, false);  return VEC;  }    // function call operator that work for the specified range (begin/end)  void operator()(std::size_t begin, std::size_t end) {  for (std::size_t i = begin; i lt; end; i  ) {  for(std::size_t j=0; jlt;U.ncol(); j  ){  arma::mat y = convertY();  arma::mat x = convertX();  arma::mat eta = convertEta();  arma::mat lambda = convertLambda();  arma::mat u = convertU();  arma::mat v = convertV();  arma::vec delta = convertDelta();    arma::vec prop = MH_U(i, j, y, x, eta, lambda, u, v, delta, tau, eps);    // write to output matrix  res(i,j) = prop(0);  res(i,j U.ncol()) = prop(i 1);  }  }  } };  NumericMatrix ParallelUpdateU(NumericMatrix Y, NumericMatrix X, NumericMatrix eta, NumericMatrix lambda,   NumericMatrix U, NumericMatrix V, NumericVector delta, double tau, double eps) {    // allocate the matrix we will return  NumericMatrix res(U.nrow(), 2*U.ncol());    std::size_t n = Y.nrow(), d = Y.ncol(), p = X.ncol(), q = eta.ncol(), r = U.ncol();    // create the worker  UpdateU updateU(Y, X, eta, lambda, U, V, delta, tau, eps, n, d, p, q, r, res);    // call it with parallelFor  std::size_t numThreads=8;  parallelFor(0, U.nrow(), updateU, numThreads);    return res; }   

Кто-нибудь знает, почему моя исходная функция RegUpdateU и распараллеленная версия ParallelUpdateU имеют одинаковую скорость?

 library(rbenchmark) n=1000 d=5 p=3 q=2 r=3 U=matrix(0.1, nrow=p, ncol=r) V=matrix(0.1, nrow=d, ncol=r) delta=rep(0.1, length=r) Y=matrix(rbinom(n*d, 1, 0.5), nrow=n) X=matrix(rbinom(n*p, 1, 0.5), nrow=n) b=matrix(rnorm(n*q), nrow=n) lambda=matrix(1, nrow=d, ncol=q) tauU=1 epsU=1  benchmark(ParallelUpdateU(Y, X, b, lambda, U, V, delta, tauU, epsU),  RegUpdateU(Y, X, b, lambda, U, V, delta, tauU, epsU))  

Я был бы ТАК признателен, если бы кто-нибудь мог мне помочь. Я потратил на это несколько недель и так и не смог в этом разобраться. Прямо сейчас весь мой пробоотборник Гиббса работает более 90 часов с данными, которые далеко не соответствуют размеру реальных данных, которые я хотел бы использовать для своей проблемы. Это мой первый пост в stackoverflow, поэтому, пожалуйста, дайте мне знать, если мне понадобится предоставить какую-либо дополнительную информацию. Спасибо!!!!!