R 的 crossprod() 在我的机器上运行速度较慢

Question

为什么

t(mat1) %*% mat2

比

crossprod(mat1, mat2)

工作得更快。后者的全部意义不就是它调用了更高效的低级例程吗？

r$> mat1 <- array(rnorm(100 * 600), dim = c(100, 600))

r$> mat2 <- array(rnorm(100 * 800), dim = c(100, 800))

r$> microbenchmark::microbenchmark(crossprod(mat1, mat2), t(mat1) %*% mat2)
Unit: milliseconds
                  expr     min       lq     mean   median       uq     max neval
 crossprod(mat1, mat2) 37.1905 44.35975 48.49262 47.61035 50.99725 76.1250   100
      t(mat1) %*% mat2 25.0813 31.68405 35.46811 35.18380 39.16565 49.8131   100

继续奔跑

r$> version
               _
platform       x86_64-w64-mingw32
arch           x86_64
os             mingw32
crt            ucrt
system         x86_64, mingw32
status
major          4
minor          4.1
year           2024
month          06
day            14
svn rev        86737
language       R
version.string R version 4.4.1 (2024-06-14 ucrt)
nickname       Race for Your Life

编辑：本地机器环境

r$> sessionInfo()
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default

r$> Sys.getenv("R_BLAS_LIBS")
    Sys.getenv("R_LAPACK_LIBS")

[1] ""
[1] ""

r$> getOption("matprod")
[1] "default"

r$> extSoftVersion()
                     zlib                     bzlib                        xz 
                  "1.3.1"      "1.0.8, 13-Jul-2019"                   "5.4.6" 
               libdeflate                      PCRE                       ICU 
                   "1.19"        "10.43 2024-02-16"                    "74.2" 
                      TRE                     iconv                  readline 
"TRE 0.8.0 R_fixes (BSD)"               "win_iconv"                        "" 
                     BLAS 
                       "" 

Answer 1

这不是答案，但评论太长了：

crossprod

和简单矩阵乘法都可以执行简单的强力三重循环，或使用适当的参数调用BLAS例程

gemm

[“通用矩阵-矩阵（乘法）”]，具体取决于（1）R如何已配置以及 (2) 参数中是否有任何

NA

值。因此，鉴于

t(x) %*% y

必然涉及分配空间来保存转置矩阵和进行实际值转置的中间步骤，人们会期望

crossprod

永远不会比转置+矩阵乘法慢；但是，如果转置步骤相对于乘法中涉及的所有浮点运算非常快，则差异可能会在噪声中丢失。

如果你想深入研究代码，你可以找到例如

• “简单”crossprod（无 BLAS）的代码如下所示：

   R_xlen_t NRX = nrx, NRY = nry, NCX = ncx;
    for (int i = 0; i < ncx; i++)
    for (int k = 0; k < ncy; k++) {
        sum = 0.0;                  
        for (int j = 0; j < nrx; j++)   
        sum += x[j + i * NRX] * y[j + k * NRY];
        z[i + k * NCX] = (double) sum;  
    }

• 类似的矩阵乘法简单循环
• 决定是否对mmult使用BLAS的逻辑（

  crossprod

  有类似的逻辑）。