Rcpp port of a linear algebra function in R

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.