Tim*_*Tim 3 r matrix blas matrix-multiplication
当恰好被crossprod(X,Y)优选t(X) %*% Y，当X和Y都是矩阵？文档说
给定矩阵
x并y作为参数，返回矩阵叉积。这在形式上等同于（但通常略快于）调用t(x) %*% y(crossprod) 或x %*% t(y)(tcrossprod)。
那么什么时候不是更快呢？在网上搜索时，我发现了几个来源，这些来源要么crossprod通常是首选的，应作为默认值使用（例如此处），要么取决于比较结果并无定论。例如，Douglas Bates (2007) 说：
请注意，在较新版本的 R 和 BLAS 库（截至 2007 年夏季）中，R
%*%能够检测到许多零mm并简化许多操作，因此对于这种稀疏矩阵来说，速度比crossprod当前没有使用的要快得多这样的优化 [...]
那么我应该什么时候使用%*%，什么时候使用crossprod？
我在几个关于统计计算问题的回答中简要介绍了这个问题。我这个回答的后续部分比较全面。请注意，我在那里的讨论以及以下内容确实假设您知道BLAS是什么，尤其是 3 级 BLAS 例程和是什么。dgemmdsyrk
我在这里的回答将提供一些证据和基准，以向您保证我在那里的论点。此外，将对 Douglas Bates 的评论进行解释。
crossprod和"%*%"真正做什么？让我们检查两个例程的源代码。R基础包的C 级源代码主要位于
R-release/src/main
特别地，矩阵运算定义在
R-release/src/main/array.c
现在，
"%*%"与 C 例程匹配，该例程使用andmatprod调用；dgemmtransa = "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肯定会快 2 倍。
具有稀疏矩阵不会改变上述基本结论。ř包Matrix用于建立稀疏矩阵还具有用于稀疏计算方法"%*%"和crossprod。所以你仍然应该期望crossprod稍微快一点。
这与稀疏矩阵无关。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)，因为前者的循环变体更有效。
最有启发性的例子是 whenB是对角矩阵或带状矩阵。让我们考虑一个带状矩阵（这里实际上是一个对称的三对角矩阵）：
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 还要快。