优化for循环RcppArmadillo中的矩阵乘法

Question

目的是在 R 中实现正交投影非负矩阵分解 (opnmf) 的快速版本。我正在运行列出可用的 matlab 代码 here。

我实现了一个普通的 R 版本，但它比我的数据 (~ 225000 x 150) 对于 20 因子解决方案的 matlab 实现要慢得多（大约慢 5.5 倍）。

所以我认为使用 c++ 可能会加快速度，但它的速度与 R 相似。我认为这可以优化但不确定如何，因为我是 c++ 的新手。 Here 是一个讨论类似问题的帖子。

这是我的 RcppArmadillo 实现。

// [[Rcpp::export]]
Rcpp::List arma_opnmf(const arma::mat & X, const arma::mat & W0, double tol=0.00001, int maxiter=10000, double eps=1e-16) {
  arma::mat W = W0;
  arma::mat Wold = W;
  arma::mat XXW = X * (X.t()*W);
  double diffW = 9999999999.9;
  
  Rcout << "The value of maxiter : " << maxiter << "\n";
  Rcout << "The value of tol : " << tol << "\n";
  
  int i;
  for (i = 0; i < maxiter; i++) {
    XXW = X * (X.t()*W);
    W = W % XXW / (W  * (W.t() * XXW));
    //W = W % (X*(X.t()*W)) / (W*((W.t()*X)*(X.t()*W)));
    
    arma::uvec idx = find(W < eps);
    W.elem(idx).fill(eps);
    W = W / norm(W,2);
    diffW = norm(Wold-W, "fro") / norm(Wold, "fro");
    if(diffW < tol) {
      break;
    } else {
      Wold = W;
    }
    
    if(i % 10 == 0) {
      Rcpp::checkUserInterrupt();
    }
    
  }
  return Rcpp::List::create(Rcpp::Named("W")=W,
                            Rcpp::Named("iter")=i,
                            Rcpp::Named("diffW")=diffW);
}

这suggested issue证实了matlab是相当快的，那么用R/c++就没有希望了吗？

测试是在 Windows 10 和 Ubuntu 16 上进行的，R 版本为 4.0.0。

编辑

在下面的答案中有趣的评论之后。我正在发布更多详细信息。我运行在装有 R 3.5.3 的 Windows 10 机器上进行测试（这是 Microsoft 提供的），比较表明 RcppArmadillo 与 Microsoft 的 R 是最快的。

R

   user  system elapsed 
 213.76    7.36  221.42

R 与 RcppArmadillo

   user  system elapsed 
 179.88    3.44  183.43

微软的 Open R

   user  system elapsed 
 167.33    9.96   45.94

Microsoft 的 Open with RcppArmadillo

    user  system elapsed 
  85.47    4.66   23.56

Answer 1

您是否知道这段代码“最终”由一对名为 LAPACK 和 BLAS 的库执行？

您知道 Matlab 附带一个高度优化的吗？您知道吗，在所有使用 R 运行的系统上，您可以更改正在使用 LAPACK/BLAS。

差异非常重要。就在今天早上，一位朋友发布了这个 tweet 对比 相同的 R 代码 运行ning 在 相同的 Windows 计算机 但在两个不同的 R 环境中。 快六倍 一个“简单地”使用并行 LAPACK/BLAS 实现。

在这里，您甚至都没有告诉我们您使用的是哪种操作系统。您可以获得适用于 R 运行的所有操作系统的 OpenBLAS（使用并行性）。您甚至可以在某些操作系统上相当轻松地获得 Intel MKL（IIRC 也是 Matlab 使用的）。对于 Ubuntu/Debian，我发布了 a script on GitHub，一步完成。

最后，很多年前我在一台（当时比较大的）Windows 计算机上“继承”了一个用 Matlab 编写的快速程序运行ning。我使用 RcppArmadillo 在 C++ 中重写了 Matlab 部分（小心而缓慢地努力），带来了一些改进因素——因为我们可以运行在同一台计算机上的另一台计算机上从 R 并行编写（现在开源）代码几个因素。将为期一天的模拟变成运行几分钟的事情，这是一个数量级的变化。所以“是的，你可以”。

编辑： 因为你可以访问 Ubuntu，你可以通过一个命令从基本的 LAPACK/BLAS 切换到 OpenBLAS，虽然我不再熟悉 Ubuntu 16.04（因为我自己运行 20.04）。

编辑2：从Josef's tweet, the Docker r-base container I also maintainer (as part of the Rocker Project)中提取比较可以使用OpenBLAS。 [1] 所以一旦我们添加它，例如通过 apt-get install libopenblas-dev 一个简单的重复矩阵叉积的时间从

移动

root@0eb44b1fcc06:/# Rscript -e 'v <- matrix(1:1e6,1e3); system.time(replicate(10, crossprod(v,v)))'
   user  system elapsed 
  9.289   0.084   9.373 
root@0eb44b1fcc06:/#

至

root@67bd334f53d4:/# Rscript -e 'v <- matrix(1:1e6,1e3); system.time(replicate(10, crossprod(v,v)))'
   user  system elapsed 
  2.259   2.370   0.447 
root@67bd334f53d4:/#

这很重要。

优化for循环RcppArmadillo中的矩阵乘法

optimize matrix multiplication in for loop RcppArmadillo

r

matrix-multiplication

nmf

rcpparmadillo