时间戳

首页 > 解决方案 > 优化 for 循环 RcppArmadillo 中的矩阵乘法

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

问题描述

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

我实现了一个香草 R 版本，但它比我的数据（~225000 x 150）的 20 因子解决方案的 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 是最快的。

R

   user  system elapsed 
 213.76    7.36  221.42

R 与 RcppArmadillo

   user  system elapsed 
 179.88    3.44  183.43

微软的开放 R

   user  system elapsed 
 167.33    9.96   45.94

微软与 RcppArmadillo 合作

    user  system elapsed 
  85.47    4.66   23.56

标签： rmatrix-multiplicationrcpparmadillonmf

解决方案

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

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

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

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

最后，多年前，我“继承”了一个在 Matlab 中运行的快速程序，该程序在一台（当时非常大的）Windows 计算机上运行。我使用 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:/#

这是实质性的。

推荐阅读