我在几个关于统计计算问题的回答中简要介绍了这个问题。 my this answer的跟进部分比较全面。请注意，我在那里的讨论以及以下内容确实假设您知道 BLAS 是什么，尤其是 3 级 BLAS 例程 dgemm 和 dsyrk 是什么。
我在这里的回答将提供一些证据和基准，以向您保证我在那里的论点。此外，将对 Douglas Bates 的评论进行解释。
crossprod 和 "%*%" 到底是做什么的？
让我们检查一下这两个例程的源代码。 R base 包的 C 级源代码主要位于
R-release/src/main
特别是，矩阵运算定义在
R-release/src/main/array.c
现在，
-
"%*%" 与 C 例程 matprod 匹配，该例程使用 transa = "N" 和 transb = "N" 调用 dgemm；
-
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计数（浮点加法和乘法的数量）@ 987654351@。如果你做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 肯定会快 2 倍。
如果你有稀疏矩阵怎么办？
拥有稀疏矩阵不会改变上述基本结论。 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 是一个三角矩阵。这在实践中相当普遍。三角矩阵不被视为稀疏矩阵（也不应该存储为稀疏矩阵）。
警告：如果您将 R 与优化的 BLAS 库（如 OpenBLAS 或 Intel MKL）一起使用，您仍然会看到 crossprod 更快。然而，这实际上是因为任何优化的 BLAS 库中的 阻塞和缓存 策略都会破坏循环变体调度模式，就像 Netlib 的 F77 参考 BLAS 中一样。因此，任何“零检测”都是不可能的。因此，您将观察到对于这个特定示例，F77 参考 BLAS 甚至比优化后的 BLAS 还要快。