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

Question

目的是在 R 中实现正交投影非负矩阵分解 (opnmf) 的快速版本。我正在翻译此处提供的 matlab 代码。

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

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

这是我的 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);
}

这个建议的问题证实了matlab相当快，那么使用R/c++时就没有希望了吗？

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

编辑

在下面的答案中出现有趣的评论之后。我正在发布更多详细信息。我在装有 R 3.5.3（微软提供的）的 Windows 10 机器上进行了测试，比较表明带有微软 R 的 RcppArmadillo 速度最快。

右

   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 

微软与 RcppArmadillo 的公开赛

    user  system elapsed 
  85.47    4.66   23.56 

Answer 1

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

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

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

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

最后，许多年前，我“继承”了一个在（当时相当大的）Windows 计算机上在 Matlab 中运行的快速程序。我使用 RcppArmadillo 用 C++ 重写了 Matlab 部分（小心而缓慢地，这是努力），这带来了一些改进因素——并且因为我们可以在同一台计算机上从 R 并行运行该代码（现在是开源的），所以还有一些因素。总而言之，将长达一天的模拟变成了运行几分钟的模拟。所以“是的，你可以”。

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

编辑 2：从Josef 的推文中进行比较，我也是维护者的 Docker r-base容器（作为Rocker 项目的一部分）可以使用 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:/#

这是很重要的。