3

The following is an objective function for symmetrical non-negative matrix factorization that I'm trying to port into Rcpp:

fit_H <- function(W,H, num.iter){
  for(i in 1:num.iter){
    H <- 0.5*(H*(1+(crossprod(W,H)/tcrossprod(H,crossprod(H)))))
  }
  H
}

H <- fit_H(W,H, num.iter = 100)

Here's my stab (with plenty of things wrong, to be sure) at an Rcpp translation:

// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd W, Eigen::MatrixXd H){
    Eigen::MatrixXd C = 0.5 * H * (1 + ((W*H)/(trans(H)*(H*H))));

    return Rcpp::wrap(C);
}

What's the correct (and fastest) single-core translation?

Theoretical background: Iteration of this function drives convergence to an H such that HH^T approximately equals W. See this question (and my answer there) for background.

zdebruine
  • 3,687
  • 6
  • 31
  • 50
  • 1
    why the for-loop? It seems to have no purpose in your first example. – Oliver Oct 31 '20 at 20:50
  • @Oliver I'd like to repeat the function `num.iter` times. With each iteration, H is being updated for use in the next iteration. I don't want to be doing the for-loop in R and passing H through Rcpp with each iteration. – zdebruine Oct 31 '20 at 20:57
  • 1
    @Oliver Ah, code corrected! I previously indicated variable `D` as what should have been `H`. One of those Stack-O copy-paste simplification jobs gone wrong :( – zdebruine Oct 31 '20 at 21:11
  • 2
    Now that you've edited your question, If the actual problem is like the one shown, I'd simply stay in R. Matrix operations are performed in optimized C code, so I wouldn't expect to see any grand (if any at all) improvements. One thing to note however, is to set `rng = false` in `Rcpp`, as the exporting of seed takes a not insignificant amount of time (see [this answer](https://stackoverflow.com/a/60288999/10782538)). It seems that you are mostly missing a for-loop and transposing some objects (which are transposed by `crossprod`). – Oliver Nov 01 '20 at 07:52
  • @Oliver thanks again! I’d mark that comment as the accepted answer if I could. I read your answer in the question you linked and it clearly shows Rcpp is not always faster than Optimized R functions. I’ll look into ways to optimize the objective function in R. – zdebruine Nov 01 '20 at 18:36

0 Answers0