improved speed through adding multiple parallel operations

R/RcppExports.R

@@ -25,6 +25,14 @@ cost_sparse_rcpp <- function(X, A, B, e, poisson) {
     .Call('_fastTopics_cost_sparse_rcpp', PACKAGE = 'fastTopics', X, A, B, e, poisson)
 }
 
+cost_parallel_rcpp <- function(X, A, B, e, poisson) {
+    .Call('_fastTopics_cost_parallel_rcpp', PACKAGE = 'fastTopics', X, A, B, e, poisson)
+}
+
+cost_sparse_parallel_rcpp <- function(X, A, B, e, poisson) {
+    .Call('_fastTopics_cost_sparse_parallel_rcpp', PACKAGE = 'fastTopics', X, A, B, e, poisson)
+}
+
 le_diff_rcpp <- function(X) {
     .Call('_fastTopics_le_diff_rcpp', PACKAGE = 'fastTopics', X)
 }
@@ -33,6 +41,14 @@ x_over_crossprod_rcpp <- function(i, j, x, A, B, e) {
     .Call('_fastTopics_x_over_crossprod_rcpp', PACKAGE = 'fastTopics', i, j, x, A, B, e)
 }
 
+x_over_crossprod_parallel_rcpp <- function(i, j, x, A, B, e) {
+    .Call('_fastTopics_x_over_crossprod_parallel_rcpp', PACKAGE = 'fastTopics', i, j, x, A, B, e)
+}
+
+poisson_nmf_kkt_sparse_parallel_rcpp <- function(X, L, F, e) {
+    .Call('_fastTopics_poisson_nmf_kkt_sparse_parallel_rcpp', PACKAGE = 'fastTopics', X, L, F, e)
+}
+
 mixem_rcpp <- function(L, w, x0, numiter) {
     .Call('_fastTopics_mixem_rcpp', PACKAGE = 'fastTopics', L, w, x0, numiter)
 }

R/likelihood.R

@@ -93,15 +93,15 @@ cost <- function (X, A, B, e = 1e-8, family = c("poisson","multinom"),
     A <- matrix(A)
   if (!is.matrix(B))
     B <- matrix(B,1,ncol(X))
-    
+
   # Check and process "model" and "version" input arguments.
   family  <- match.arg(family)
   poisson <- family == "poisson"
   if (missing(version)) {
     if (is.matrix(X))
       version <- "R"
     else
-      version <- "Rcpp"
+      version <- "Rcpp_parallel"
   }
 
   # Compute the terms in the log-likelihoods that depend on A or B.
@@ -113,6 +113,11 @@ cost <- function (X, A, B, e = 1e-8, family = c("poisson","multinom"),
       f <- drop(cost_rcpp(X,A,B,e,poisson))
     else
       f <- drop(cost_sparse_rcpp(X,A,B,e,poisson))
+  } else if (version == "Rcpp_parallel") {
+    if (is.matrix(X))
+      f <- drop(cost_parallel_rcpp(X,A,B,e,poisson))
+    else
+      f <- drop(cost_sparse_parallel_rcpp(X,A,B,e,poisson))
   }
   names(f) <- rownames(X)
   return(f)
@@ -179,15 +184,24 @@ deviance_poisson_const <- function (X) {
 # Compute the residuals of the first-order Karush-Kuhn-Tucker (KKT)
 # conditions for Poisson non-negative matrix factorization at solution
 # estimate (F,L).
-poisson_nmf_kkt <- function (X, F, L, e = 1e-8) {
+#' @importFrom Matrix sparseMatrix
+poisson_nmf_kkt <- function (X, F, L, e = 1e-8,
+                             version = c("Rcpp_parallel","Rcpp")) {
+  version <- match.arg(version)
   if (is.matrix(X)) {
     A <- X/(tcrossprod(L,F) + e)
     return(list(F = F*(t(1 - A) %*% L),
                 L = L*((1 - A) %*% F)))
-  } else {
-    A <- x_over_tcrossprod(X,t(L),t(F),e)
+  } else if (version == "Rcpp") {
+    d <- summary(X)
+    y <- drop(x_over_crossprod_rcpp(d$i - 1,d$j - 1,d$x,t(L),t(F),e))
+    A <- sparseMatrix(i = d$i,j = d$j,x = y,dims = dim(X))
     return(list(F = F*(repmat(colSums(L),ncol(X)) - as.matrix(t(A) %*% L)),
                 L = L*(repmat(colSums(F),nrow(X)) - as.matrix(A %*% F))))
+  } else {
+    result <- poisson_nmf_kkt_sparse_parallel_rcpp(X,L,F,e)
+    return(list(F = F*(repmat(colSums(L),ncol(X)) - result$tAL),
+                L = L*(repmat(colSums(F),nrow(X)) - result$AF)))
   }
 }
 

R/misc.R

@@ -55,7 +55,7 @@ apply.nonzeros <- function (X, f) {
 #' 
 x_over_tcrossprod <- function (X, A, B, e) {
   d <- summary(X)
-  y <- drop(x_over_crossprod_rcpp(d$i - 1,d$j - 1,d$x,A,B,e))
+  y <- drop(x_over_crossprod_parallel_rcpp(d$i - 1,d$j - 1,d$x,A,B,e))
   return(sparseMatrix(i = d$i,j = d$j,x = y,dims = dim(X)))
 }
 

src/RcppExports.cpp

@@ -97,6 +97,36 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cost_parallel_rcpp
+arma::vec cost_parallel_rcpp(const arma::mat& X, const arma::mat& A, const arma::mat& B, double e, bool poisson);
+RcppExport SEXP _fastTopics_cost_parallel_rcpp(SEXP XSEXP, SEXP ASEXP, SEXP BSEXP, SEXP eSEXP, SEXP poissonSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< double >::type e(eSEXP);
+    Rcpp::traits::input_parameter< bool >::type poisson(poissonSEXP);
+    rcpp_result_gen = Rcpp::wrap(cost_parallel_rcpp(X, A, B, e, poisson));
+    return rcpp_result_gen;
+END_RCPP
+}
+// cost_sparse_parallel_rcpp
+arma::vec cost_sparse_parallel_rcpp(const arma::sp_mat& X, const arma::mat& A, const arma::mat& B, double e, bool poisson);
+RcppExport SEXP _fastTopics_cost_sparse_parallel_rcpp(SEXP XSEXP, SEXP ASEXP, SEXP BSEXP, SEXP eSEXP, SEXP poissonSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::sp_mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< double >::type e(eSEXP);
+    Rcpp::traits::input_parameter< bool >::type poisson(poissonSEXP);
+    rcpp_result_gen = Rcpp::wrap(cost_sparse_parallel_rcpp(X, A, B, e, poisson));
+    return rcpp_result_gen;
+END_RCPP
+}
 // le_diff_rcpp
 arma::mat le_diff_rcpp(const arma::mat& X);
 RcppExport SEXP _fastTopics_le_diff_rcpp(SEXP XSEXP) {
@@ -124,6 +154,36 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// x_over_crossprod_parallel_rcpp
+arma::vec x_over_crossprod_parallel_rcpp(const arma::vec& i, const arma::vec& j, const arma::vec& x, const arma::mat& A, const arma::mat& B, double e);
+RcppExport SEXP _fastTopics_x_over_crossprod_parallel_rcpp(SEXP iSEXP, SEXP jSEXP, SEXP xSEXP, SEXP ASEXP, SEXP BSEXP, SEXP eSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::vec& >::type i(iSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type j(jSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type A(ASEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+    Rcpp::traits::input_parameter< double >::type e(eSEXP);
+    rcpp_result_gen = Rcpp::wrap(x_over_crossprod_parallel_rcpp(i, j, x, A, B, e));
+    return rcpp_result_gen;
+END_RCPP
+}
+// poisson_nmf_kkt_sparse_parallel_rcpp
+Rcpp::List poisson_nmf_kkt_sparse_parallel_rcpp(const arma::sp_mat& X, const arma::mat& L, const arma::mat& F, double e);
+RcppExport SEXP _fastTopics_poisson_nmf_kkt_sparse_parallel_rcpp(SEXP XSEXP, SEXP LSEXP, SEXP FSEXP, SEXP eSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::sp_mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type L(LSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type F(FSEXP);
+    Rcpp::traits::input_parameter< double >::type e(eSEXP);
+    rcpp_result_gen = Rcpp::wrap(poisson_nmf_kkt_sparse_parallel_rcpp(X, L, F, e));
+    return rcpp_result_gen;
+END_RCPP
+}
 // mixem_rcpp
 arma::vec mixem_rcpp(const arma::mat& L, const arma::vec& w, const arma::vec& x0, unsigned int numiter);
 RcppExport SEXP _fastTopics_mixem_rcpp(SEXP LSEXP, SEXP wSEXP, SEXP x0SEXP, SEXP numiterSEXP) {
@@ -414,8 +474,12 @@ static const R_CallMethodDef CallEntries[] = {
     {"_fastTopics_ccd_update_factors_sparse_parallel_rcpp", (DL_FUNC) &_fastTopics_ccd_update_factors_sparse_parallel_rcpp, 4},
     {"_fastTopics_cost_rcpp", (DL_FUNC) &_fastTopics_cost_rcpp, 5},
     {"_fastTopics_cost_sparse_rcpp", (DL_FUNC) &_fastTopics_cost_sparse_rcpp, 5},
+    {"_fastTopics_cost_parallel_rcpp", (DL_FUNC) &_fastTopics_cost_parallel_rcpp, 5},
+    {"_fastTopics_cost_sparse_parallel_rcpp", (DL_FUNC) &_fastTopics_cost_sparse_parallel_rcpp, 5},
     {"_fastTopics_le_diff_rcpp", (DL_FUNC) &_fastTopics_le_diff_rcpp, 1},
     {"_fastTopics_x_over_crossprod_rcpp", (DL_FUNC) &_fastTopics_x_over_crossprod_rcpp, 6},
+    {"_fastTopics_x_over_crossprod_parallel_rcpp", (DL_FUNC) &_fastTopics_x_over_crossprod_parallel_rcpp, 6},
+    {"_fastTopics_poisson_nmf_kkt_sparse_parallel_rcpp", (DL_FUNC) &_fastTopics_poisson_nmf_kkt_sparse_parallel_rcpp, 4},
     {"_fastTopics_mixem_rcpp", (DL_FUNC) &_fastTopics_mixem_rcpp, 4},
     {"_fastTopics_pnmfem_update_factors_rcpp", (DL_FUNC) &_fastTopics_pnmfem_update_factors_rcpp, 5},
     {"_fastTopics_pnmfem_update_factors_sparse_rcpp", (DL_FUNC) &_fastTopics_pnmfem_update_factors_sparse_rcpp, 5},

src/cost.cpp

@@ -1,8 +1,81 @@
+#include <RcppParallel.h>
 #include "misc.h"
 #include "cost.h"
 
 using namespace arma;
 
+// CLASS DEFINITIONS
+// -----------------
+// Worker for parallel computation of the dense cost function using
+// RcppParallel reduce. Each thread accumulates its own partial sum
+// over a range of columns of X, and the join() method combines them.
+struct cost_worker : public RcppParallel::Worker {
+  const mat&  X;
+  const mat&  A;
+  const mat&  B;
+  double      e;
+  bool        poisson;
+  vec         f;
+
+  cost_worker (const mat& X, const mat& A, const mat& B, double e,
+               bool poisson) :
+    X(X), A(A), B(B), e(e), poisson(poisson), f(X.n_rows, fill::zeros) { };
+
+  // Split constructor: creates a new worker with a zeroed accumulator.
+  cost_worker (const cost_worker& other, RcppParallel::Split) :
+    X(other.X), A(other.A), B(other.B), e(other.e), poisson(other.poisson),
+    f(other.X.n_rows, fill::zeros) { };
+
+  // Accumulate contributions from columns [begin, end).
+  void operator() (std::size_t begin, std::size_t end) {
+    vec y(A.n_rows);
+    for (unsigned int j = begin; j < end; j++) {
+      y  = A * B.col(j);
+      f -= X.col(j) % log(y + e);
+      if (poisson)
+        f += y;
+    }
+  }
+
+  void join (const cost_worker& other) { f += other.f; }
+};
+
+// Worker for parallel computation of the sparse cost function.
+struct cost_sparse_worker : public RcppParallel::Worker {
+  const sp_mat& X;
+  const mat&    A;
+  const mat&    B;
+  double        e;
+  bool          poisson;
+  vec           f;
+
+  cost_sparse_worker (const sp_mat& X, const mat& A, const mat& B, double e,
+                      bool poisson) :
+    X(X), A(A), B(B), e(e), poisson(poisson), f(X.n_rows, fill::zeros) { };
+
+  cost_sparse_worker (const cost_sparse_worker& other, RcppParallel::Split) :
+    X(other.X), A(other.A), B(other.B), e(other.e), poisson(other.poisson),
+    f(other.X.n_rows, fill::zeros) { };
+
+  void operator() (std::size_t begin, std::size_t end) {
+    unsigned int i;
+    vec y(A.n_rows);
+    for (unsigned int j = begin; j < end; j++) {
+      sp_mat::const_col_iterator xj = X.begin_col(j);
+      sp_mat::const_col_iterator xm = X.end_col(j);
+      y = A * B.col(j);
+      for (; xj != xm; ++xj) {
+        i     = xj.row();
+        f(i) -= (*xj) * log(y(i) + e);
+      }
+      if (poisson)
+        f += y;
+    }
+  }
+
+  void join (const cost_sparse_worker& other) { f += other.f; }
+};
+
 // FUNCTION DEFINITIONS
 // --------------------
 // Compute negative log-likelihoods for assessing a topic model fit or
@@ -52,6 +125,34 @@ arma::vec cost (const mat& X, const mat& A, const mat& B, double e,
   return f;
 }
 
+// This is the same as cost_rcpp, except that Intel Threading Building
+// Blocks (TBB) are used to compute the contributions from each column
+// of X in parallel.
+//
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::depends(RcppParallel)]]
+// [[Rcpp::export]]
+arma::vec cost_parallel_rcpp (const arma::mat& X, const arma::mat& A,
+			      const arma::mat& B, double e, bool poisson) {
+  cost_worker worker(X,A,B,e,poisson);
+  parallelReduce(0,X.n_cols,worker);
+  return worker.f;
+}
+
+// This is the same as cost_sparse_rcpp, except that Intel Threading
+// Building Blocks (TBB) are used to compute the contributions from
+// each column of X in parallel.
+//
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::depends(RcppParallel)]]
+// [[Rcpp::export]]
+arma::vec cost_sparse_parallel_rcpp (const arma::sp_mat& X, const arma::mat& A,
+				     const arma::mat& B, double e, bool poisson) {
+  cost_sparse_worker worker(X,A,B,e,poisson);
+  parallelReduce(0,X.n_cols,worker);
+  return worker.f;
+}
+
 // Helper function for cost_sparse_rcpp.
 arma::vec cost_sparse (const sp_mat& X, const mat& A, const mat& B,
 		       double e, bool poisson) {