什么时候 'crossprod' 比 '%*%' 更受欢迎，什么时候不是？
When is 'crossprod' preferred to '%*%', and when isn't?
当 X 和 Y 都是矩阵时，crossprod(X,Y) 何时比 t(X) %*% Y 更受欢迎？文档说
Given matrices x and y as arguments, return a matrix cross-product.
This is formally equivalent to (but usually slightly faster than) the
call t(x) %*% y (crossprod) or x %*% t(y) (tcrossprod).
那么什么时候不更快？在网上搜索时，我发现有几个来源表明 crossprod 通常是首选并且应该用作默认值（例如 here), or that it depends and results of comparisons are inconclusive. For example, Douglas Bates (2007) says that:
Note that in newer versions of R and the BLAS library (as of summer
2007), R’s %*% is able to detect the many zeros in mm and shortcut
many operations, and is hence much faster for such a sparse matrix
than crossprod which currently does not make use of such
optimizations [...]
那么我应该什么时候使用 %*% 什么时候使用 crossprod？
我在关于统计计算问题的几个回答中简要介绍了这个问题。 的跟进部分比较全面。请注意，我在那里的讨论以及以下内容实际上假设您知道 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 计数（浮点加法和乘法的数量）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 肯定会快 2 倍。
如果你有稀疏矩阵怎么办？
稀疏矩阵不会改变上面的基本结论。用于设置稀疏矩阵的 R 包 Matrix 也有 "%*%" 和 crossprod 的稀疏计算方法。所以你仍然应该期望 crossprod 会稍微快一些。
那么贝茨对 BLAS "as of summer 2007" 的评论是什么意思？
这与稀疏矩阵无关。 BLAS 严格用于密集数值线性代数。
它涉及的是 Netlib 的 F77 参考文献 dgemm 中使用的循环变体的差异。调度矩阵-矩阵乘法有两种循环变体op(A) * op(B)（这里的*表示矩阵乘法而不是元素乘积！），变体完全由op(A)的转置设置决定：
- for
op(A) = A',最内层循环使用"inner product"版本，这种情况下不能检测零；
- 对于
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 或英特尔 MKL）一起使用，您仍然会发现 crossprod 更快。然而，这实际上是因为任何优化的 BLAS 库中的 阻塞和缓存 策略破坏了 Netlib 的 F77 参考 BLAS 中的循环变体调度模式。结果，任何 "zero-detection" 是不可能的。因此，您将观察到 对于这个特定示例，F77 参考 BLAS 甚至比优化的 BLAS 更快。
当 X 和 Y 都是矩阵时，crossprod(X,Y) 何时比 t(X) %*% Y 更受欢迎？文档说
Given matrices
xandyas arguments, return a matrix cross-product. This is formally equivalent to (but usually slightly faster than) the callt(x) %*% y(crossprod) orx %*% t(y)(tcrossprod).
那么什么时候不更快？在网上搜索时，我发现有几个来源表明 crossprod 通常是首选并且应该用作默认值（例如 here), or that it depends and results of comparisons are inconclusive. For example, Douglas Bates (2007) says that:
Note that in newer versions of R and the BLAS library (as of summer 2007), R’s
%*%is able to detect the many zeros inmmand shortcut many operations, and is hence much faster for such a sparse matrix thancrossprodwhich currently does not make use of such optimizations [...]
那么我应该什么时候使用 %*% 什么时候使用 crossprod？
我在关于统计计算问题的几个回答中简要介绍了这个问题。 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 计数（浮点加法和乘法的数量）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 肯定会快 2 倍。
如果你有稀疏矩阵怎么办？
稀疏矩阵不会改变上面的基本结论。用于设置稀疏矩阵的 R 包 Matrix 也有 "%*%" 和 crossprod 的稀疏计算方法。所以你仍然应该期望 crossprod 会稍微快一些。
那么贝茨对 BLAS "as of summer 2007" 的评论是什么意思？
这与稀疏矩阵无关。 BLAS 严格用于密集数值线性代数。
它涉及的是 Netlib 的 F77 参考文献 dgemm 中使用的循环变体的差异。调度矩阵-矩阵乘法有两种循环变体op(A) * op(B)（这里的*表示矩阵乘法而不是元素乘积！），变体完全由op(A)的转置设置决定：
- for
op(A) = A',最内层循环使用"inner product"版本，这种情况下不能检测零； - 对于
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 或英特尔 MKL）一起使用，您仍然会发现 crossprod 更快。然而，这实际上是因为任何优化的 BLAS 库中的 阻塞和缓存 策略破坏了 Netlib 的 F77 参考 BLAS 中的循环变体调度模式。结果，任何 "zero-detection" 是不可能的。因此，您将观察到 对于这个特定示例，F77 参考 BLAS 甚至比优化的 BLAS 更快。