我在几个统计计算问题的回答中简要介绍了这个问题。
我的这个回答 的后续部分比较全面。请注意,我在那里的讨论以及以下内容实际上假定您知道什么是BLAS,特别是什么是3级BLAS例程
dgemm
和
dsyrk
。
在这里,我的答案将提供一些证据和基准测试,以确保您对我在那里的论点有信心。此外,还将解释Douglas Bates的评论。
< h1 > 交叉积和< code >"%*%"到底是什么?
让我们检查这两个例程的源代码。 R
base包的C级源代码主要可以在以下位置找到:
R-release/src/main
特别地,矩阵运算定义在
R-release/src/main/array.c
现在,
"%*%"
匹配到一个C例程matprod
,该例程调用dgemm
并设置transa = "N"
和transb = "N"
;
crossprod
很容易被视为一个复合函数,匹配到symcrossprod
,crossprod
,symtcrossprod
和tcrossprod
(复杂矩阵的对应物,如ccrossprod
在此未列出)。
你现在应该明白crossprod
避免了所有显式的矩阵转置。 crossprod(A,B)
比t(A)%*%B
更便宜,因为后者需要对A
进行显式矩阵转置。 在接下来的内容中,我将把这称为转置开销。
在R级别分析可以展示这个开销。考虑以下例子:
A <- matrix(ruinf(1000 * 1000), 1000)
B <- matrix(ruinf(1000 * 1000), 1000)
Rprof("brutal.out")
result <- t.default(A) %*% B
Rprof(NULL)
summaryRprof("brutal.out")
Rprof("smart.out")
result <- crossprod(A, B)
Rprof(NULL)
summaryRprof("smart.out")
请注意第一个案例中的
t.default
如何出现在分析结果中。
分析还给出了执行的CPU时间。您会发现两者似乎花费了相同的时间,因为开销微不足道。现在,我将向您展示什么时候开销是一个痛点。
何时转置开销才显著?
假设
A
是一个
k * m
的矩阵,
B
是一个
k * n
的矩阵,则矩阵乘法
A'B
(
'
表示转置)具有FLOP计数(浮点加法和乘法量)
2mnk
。如果您执行
t(A) %*% B
,则转置开销为
mk
(
A
中的元素数),因此比率为:
useful computation : overhead = 2n : 1
除非n
很大,否则在实际矩阵乘法中无法摊销转置开销。
最极端的情况是n=1
,即对于B
你有一个单列矩阵(或向量)。考虑进行基准测试:
library(microbenchmark)
A <- matrix(runif(2000 * 2000), 2000)
x <- runif(2000)
microbenchmark(t.default(A) %*% x, crossprod(A, x))
你会发现
crossprod
快了好几倍!
什么情况下"%*%"
结构上不如其他方法?
正如我在链接答案中提到的(以及Bates的笔记中的基准测试结果),如果你做A'A
,crossprod
肯定会快两倍。
如果你有稀疏矩阵怎么办?
拥有稀疏矩阵并不改变上述基本结论。R包Matrix
用于设置稀疏矩阵,还具有用于"%*%"
和crossprod
的稀疏计算方法。所以你仍然应该期望crossprod
稍微快一些。
Bates关于BLAS“截至2007年夏季”的评论是什么意思?
这与稀疏矩阵无关。BLAS严格用于密集数值线性代数。
它所涉及的是在Netlib的F77参考dgemm
内部使用的循环变量的差异。对于调度矩阵-矩阵乘法op(A) * op(B)
(这里的*
表示矩阵乘法而不是元素乘积!),有两个循环变量,变量的选择完全由op(A)
的转置设置决定:
- 对于
op(A) = A'
,在最内层循环中使用“内积”版本,在这种情况下无法进行零检测;
- 对于
op(A) = A
,使用“AXPY”版本,并且可以从计算中排除op(B)
中的任何零。
现在想想R如何调用dgemm
。第一种情况对应于crossprod
,而第二种情况对应于"%*%"
(以及tcrossprod
)。
在这方面,如果你的B
矩阵有很多零,而它仍然处于密集矩阵格式中,那么t(A) %*% B
将比crossprod(A, B)
更快,因为前者的循环变量更有效率。
最有启发性的例子是当B
是对角矩阵或带状矩阵时。让我们考虑一个带状矩阵(实际上是一个对称三对角矩阵):
B_diag <- diag(runif(1000))
B_subdiag <- rbind(0, cbind(diag(runif(999)), 0))
B <- B_diag + B_subdiag + t(B_subdiag)
现在让A
仍然是一个完整的密集矩阵。
A <- matrix(runif(1000 * 1000), 1000)
然后比较。
library(microbenchmark)
microbenchmark(t.default(A) %*% B, crossprod(A, B))
你会发现 "%*%"
要快得多!
我猜一个更好的例子是矩阵 B
是一个三角矩阵。这在实践中相当常见。三角矩阵不被视为稀疏矩阵(也不应该被存储为稀疏矩阵)。
注意:如果您正在使用带有优化BLAS库(如OpenBLAS或Intel MKL)的R,您仍然会发现crossprod
更快。但是,这确实是因为任何优化的BLAS库中的阻塞和缓存策略破坏了与Netlib的F77参考BLAS中的循环变量调度模式相同的模式。因此,任何“零检测”都是不可能的。所以,您将观察到对于这个特定的示例,F77参考BLAS甚至比优化的BLAS更快。
crossprod
是 C 级别的,因此通常应该更快,而对于稀疏矩阵,R 似乎以更高效的方式处理它。 - David Arenburg%*%
也可以是.Primitive()
和C级别的。所以答案是%*%
只有在期望稀疏矩阵时才应该使用(或者为了代码可读性)? - Tim