r - R中线性代数函数的Rcpp端口

Question

以下是我试图移植到 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)

这是我在 Rcpp 翻译中的刺探（肯定有很多错误）：

// [[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);
}

什么是正确（和最快）的单核翻译？

理论背景：这个函数的迭代驱动收敛到一个H近似HH^T等于的W。有关背景，请参阅此问题（以及我的答案）。