在 R 中计算两个瘦高矩阵之间的逐行点积的最快方法答案

【问题标题】：Fastest way to compute row-wise dot products between two skinny tall matrices in R在 R 中计算两个瘦高矩阵之间的逐行点积的最快方法
【发布时间】：2017-03-06 19:57:07
【问题描述】：

考虑 A 和 B 是两个尺寸为 10^8 X 5 的高瘦矩阵。即;

r=10^8
c=5
A=matrix(runif(r*c,0,1),r,c)
B=matrix(runif(r*c,0,1),r,c)

我想计算A的每一行与B的对应行的点积，即；

rowSums(A * B)

但它很慢，我想知道是否有更快的方法。

【问题讨论】：

标签： r

【解决方案1】：

这可能会让您失望，但在 R 级别，这已经是您无需自己编写一些 C 代码即可获得的最佳效果。问题是通过rowSums(A * B)，你实际上是在做

C <- A * B
rowSums(C)

第一行对三个大而细的矩阵进行全扫描；而第二行则对 1 个大而细的矩阵进行全扫描。所以总的来说，我们等效地扫描了一个又高又薄的矩阵 4 次（内存密集型）。

实际上，对于这样的操作，最优算法只需要扫描n * p高瘦矩阵两次，通过rowwise cross product：

rowsum <- numeric(n)
for j = 1, 2, ... p
  rowsum += A[,i] * B[,i]

这样，我们也避免了生成矩阵C。请注意，以上只是一个假代码，而不是有效的 R 代码甚至 C 代码。但是这个想法很明确，我们想用 C 来编程。

与您的情况类似的是sum(x * y) 和crossprod(x, y) 之间的速度差异，假设x 和y 是相同长度的大向量。

x <- runif(1e+7)
y <- runif(1e+7)

system.time(sum(x * y))
#   user  system elapsed 
#  0.124   0.032   0.158 

system.time(crossprod(x, y))
#   user  system elapsed 
#  0.036   0.000   0.056

在第一种情况下，我们扫描了一个长向量 4 次，而在第二种情况下，我们只扫描了两次。

统计计算的相关性

rowSums(A * B) 实际上是对diag(tcrossprod(A, B)) 的有效评估，常见于与逐点预测方差相关的回归计算中。例如，在使用来自模型矩阵的 QR 分解的细Q 矩阵的普通线性平方回归中，拟合值的逐点方差为diag(tcrossprod(Q))，由rowSums(Q ^ 2) 计算更有效。但是，这仍然不是最快的评估，原因已经解释过。

【讨论】：