矩阵与向量相乘 - R 与 Matlab 的比较

Question

我注意到在R（版本3.6.1）中，使用较大矩阵进行矩阵乘向量的操作速度明显比Matlab（版本2019b）慢，尽管两种语言都依赖于（同一个？）BLAS库。以下是一个最简示例：

在Matlab中：

n=900; 
p=900; 
A=reshape(1:(n*p),[n,p]); 
x=ones(p,1); 
tic()
for id = 1:1000
  x = A*x; 
end
toc()

R语言中：

n=900
p=900
A=matrix(c(1:(n*p)),nrow=n,ncol=p)
x=rep(1,ncol(A))
t0 <- Sys.time()
for(iter in 1:1000){
  x = A%*%x
}
t1 <- Sys.time()
print(t1-t0)

我在使用相同计算机时，在Matlab中运行时间大约为0.05秒，而在R中为3.5秒。有什么原因导致这样的差异吗？

谢谢。

[编辑]: 我添加了下面在C语言中使用CBLAS库进行的类似计算（使用gcc cblas_dgemv.c -lblas -o cblas_dgemv进行编译，其中cblas_dgemv.c表示下面的源文件）。我得到的运行时间大约为0.08秒，非常接近使用Matlab（0.05秒）获得的运行时间。我仍在尝试找出R中巨大运行时间的原因（3.5秒）。

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include "cblas.h"

int main(int argc, char **argv)
{
  int m=900,adr;
  double *A,*x,*y;
  struct timeval t0,t1;

  /* memory allocation and initialization */
  A = (double*)malloc(m*m*sizeof(double)); 
  x = (double*)malloc(m*sizeof(double));  
  y = (double*)malloc(m*sizeof(double));  
  for(adr=0;adr<m;adr++) x[adr] = 1.; 
  for(adr=0;adr<m*m;adr++) A[adr] = adr;

  /* main loop */
  gettimeofday(&t0, NULL);
  for(adr=0;adr<1000;adr++)
    cblas_dgemv(CblasColMajor,CblasNoTrans,m,m,1.,A,m,x,1,0.,y,1);
  gettimeofday(&t1, NULL);
  printf("elapsed time = %.2e seconds\n",(double)(t1.tv_usec-t0.tv_usec)/1000000. + (double)(t1.tv_sec-t0.tv_sec));

  /* free memory */
  free(A);
  free(x);
  free(y); 

  return EXIT_SUCCESS;
}

请注意，在cblas_dgemv例程中，我无法设置y=x。因此，这个C语言计算与上述R和Matlab代码略有不同。但是编译时没有使用优化标志（没有-O3选项），并且我检查了矩阵向量积确实在循环的每次迭代中被调用（执行10倍更多的迭代会导致运行时间增加10倍）。

- RémyA

---

Matlab是高度优化的。曾经我使用GNU Octave工作时，同样的代码但是Matlab运行速度要快得多。 - paul-shuvo

据我所知，无论是哪种语言，线性代数运算都依赖于BLAS/LAPACK库。 - RémyA

如果在R中运行sessionInfo()，它会显示“矩阵乘积”的内容是什么？ - MrFlick

可能相关：https://dev59.com/0Kzka4cB1Zd3GeqP4SkK - Andrew Janke

@AndrewJanke：我认为我们不会面临mayHaveNaNOrInf测试的误报响应，因为在这种情况下，矩阵-向量乘积应该使用内部未优化的3循环算法来完成，这将更慢。实际上，通过设置options(matprod = "internal")，我得到了超过1分钟的运行时间！此外，根据(https://stat.ethz.ch/R-manual/R-devel/library/base/html/options.html)，设置`options(matprod = "blas")应该可以避免这个mayHaveNaNOrInf`测试。这真是令人沮丧！ - RémyA

Answer 1

这里 有一些相当令人震惊的内容：

从CRAN下载的预编译R发行版使用参考BLAS/LAPACK实现进行线性代数操作

"参考BLAS"是未经优化、未加速的BLAS，不像OpenBLAS或Intel MKL。Matlab使用MKL，它是加速的。

这似乎在我的macOS上的R 3.6.0中通过sessionInfo()得到了确认：

> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

如果我理解正确，这意味着默认情况下，R使用缓慢的BLAS。如果您想提高速度，需要进行一些配置以改用快速的BLAS。

对我来说，这有点令人惊讶。据我所知，参考BLAS通常主要用于测试和开发，而不是用于“实际工作”。

在macOS 10.14上，我在R 3.6与Matlab R2019b中得到大约相同的时间：Matlab为0.04秒，而R为4.5秒。我认为这表明R正在使用未加速的BLAS。

- Andrew Janke

---

我链接的基准测试显示，在矩阵叉积乘法方面，MKL比参考BLAS快约40倍，这应该是我们在这里所做的。听起来差不多。对我来说，进行微不足道的矩阵乘法可以将其降至0.001秒。 - Andrew Janke

@Andrew：确实，使用MKL加速40倍可以解释我们的观察结果。我会朝这个方向进行调查，谢谢。 - RémyA

哇，参考版本和优化版本之间的差异令人惊讶。我现在闭嘴 :) - Cris Luengo

现在我有点困惑：我切换到了一个使用OpenBLAS构建的Homebrewed R，而且sessionInfo看起来也在使用它：BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.7/lib/libopenblasp-r0.3.7.dylib 但是速度仍然很慢。这与你所看到的情况类似。这需要进一步调查。 - Andrew Janke

R是否以不同的方式存储矩阵？也许需要复制才能将其转换为可以与BLAS一起使用的形式？ - Cris Luengo

Answer 2

我使用Rcpp将C++实现的矩阵-向量乘法与R进行接口交互。

源C++文件（'cblas_dgemv2.cpp'）：执行'niter'次乘积y = A*x

#include <cblas-openblas.h>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::NumericVector cblas_dgemv2(Rcpp::NumericMatrix A, Rcpp::NumericVector x, int niter)
{
  int m=A.ncol(),iter;
  Rcpp::NumericVector y(m);
  for(iter=0;iter<niter;iter++)
    cblas_dgemv(CblasColMajor,CblasNoTrans,m,m,1.,A.begin(),m,x.begin(),1,0.,y.begin(),1);
  return y; 
}

然后我使用以下R代码进行了两个实验：

• 实验1：从R中调用y=cblas_dgmev2(A,x,1000)，在C++中执行1000次计算y=Ax*的操作；
• 实验2：从R中调用1000次y=cblas_dgemv2(A,x,1)，每次调用在C++中执行y=A*x的操作。

# compile cblas_dgemv2 (you may need to update the path to the library)
PKG_LIBS <- '/usr/lib/x86_64-linux-gnu/libopenblas.a' 
PKG_CPPFLAGS <- '-I/usr/include/x86_64-linux-gnu'
Sys.setenv(PKG_LIBS = PKG_LIBS , PKG_CPPFLAGS = PKG_CPPFLAGS) 
Rcpp::sourceCpp('cblas_dgemv2.cpp', rebuild = TRUE)

# create A and x 
n=900
A=matrix(c(1:(n*n)),nrow=n,ncol=n)
x=rep(1,n)

# Experiment 1: call 1 time cblas_dgemv2 (with 1000 iterations in C++)
t0 <- Sys.time()
y=cblas_dgemv2(A,x,1000) # perform 1000 times the computation y = A%*%x 
t1 <- Sys.time()
print(t1-t0)

# Experiment 2: call 1000 times cblas_dgemv2  
t0 <- Sys.time()
for(iter in 1:1000){
  y=cblas_dgemv2(A,x,1) # perform 1 times the computation y = A%*%x 
}
t1 <- Sys.time()
print(t1-t0)

第一个实验的运行时间为0.08秒，而第二个实验的运行时间为4.8秒。

我的结论是：在运行时间方面的瓶颈来自于R和C++之间的数据传输，而不是矩阵向量乘积本身的计算。令人惊讶，不是吗？

- RémyA

---

这是一个非常大的开销。R真的在复制所有数据吗？这毫无意义！（顺便说一句，实验做得很好！） - Cris Luengo