向量化（或加速）双循环，对 R 中的不同索引求和

Question

我正在尝试优化旨在计算两个方阵元素乘积的双倍和的代码。假设我们有两个大小为 n、W 和 V 的方阵。需要计算的对象是一个带有元素的向量B

简单来说：计算两个不同矩阵中两个不同行的逐个元素乘积并取它们的总和，然后对第二个矩阵的所有行取一个额外的总和（没有相同的索引）。

问题是，这个任务的计算复杂度貌似 O(n³) 因为我们正在创建的这个对象的长度，B，是n，每个元素需要两次求和。这是我想出的：

1. 对于给定的 i 和 j (i≠j)，从 k 上的内和开始。对所有k求和，然后减去k=i和k=j的项，再乘以j≠的指标我。
2. 由于限制 j≠i 已在内部和中得到处理，因此外部和仅用于 j=1,...,n。

如果我们表示，那么这两个步骤将如下所示 和 。

然而，编写一个循环被证明是非常低效的。 n=100 工作速度很快（0.05 秒）。但是，例如，当 n=500（我们在这里讨论的是实际应用）时，平均计算时间是 3 秒，而对于 n=1000，它跳转到 22 秒。

k 上的内层循环可以很容易地用 sum 代替，但外层循环...在所有元素上完成。

这是我试图在大n宇宙热寂之前评估的代码。

set.seed(1)
N <- 500
x1 <- rnorm(N)
x2 <- rchisq(N, df=3)

bw1 <- bw.nrd(x1)
bw2 <- bw.nrd(x2)
w <- outer(x1, x1, function(x, y) dnorm((x-y)/bw1) )
w <- w/rowSums(w)

v <- outer(x2, x2, function(x, y) dnorm((x-y)/bw2) )
v <- v/rowSums(v)

Bij <- matrix(NA, ncol=N, nrow=N)
for (i in 1:N) { # Around 22 secs for N=1000
  for (j in 1:N) {
    Bij[i, j] <- (sum(w[i, ]*v[j, ]) - w[i, i]*v[j, i] - w[i, j]*v[j, j]) * (i!=j)
  }
}
Bi <- rowSums(Bij)

专家 R 程序员如何矢量化这种循环？

Answer 1

无需查看矩阵 w 和 v 的内容，您的双 for 循环可以替换为简单的矩阵运算，使用一个矩阵乘法 (tcrossprod)、转置 (t) 和对角线提取：

Mat.ij <- tcrossprod(w, v) - 
    matrix(rep(diag(w), times = N), nrow = N) * t(v) - 
    w * matrix(rep(diag(v), each = N), nrow = N)
diag(Mat.ij) <- 0

all.equal(Bij, Mat.ij)
[1] TRUE

Answer 2

更新：

事实上，鉴于你对 B_{ij} 的表达，我们也可以这样做

diag(w) <- diag(v) <- 0
BBij <-  tcrossprod(w, v)
diag(BBij) <- 0

range(rowSums(BBij) - Bi)
# [1] -2.220446e-16  0.000000e+00
range(BBij - Bij)
# [1] -6.938894e-18  5.204170e-18

因此，虽然有些明显，但对于您的目的来说，B_{ij} 和 B_i 都不依赖于 W 和 V 的对角线可能是一个有趣的观察。

初步答案：

自从 其中 W 和 V 的对角线可以设置为零，V_{.k} 表示 V 的第 k 列之和，我们有

diag(w) <- diag(v) <- 0 
A <- w * v
rowSums(sweep(w, 2, colSums(v), `*`)) - rowSums(A) + diag(A)

在哪里

range(rowSums(sweep(w, 2, colSums(v), `*`)) - rowSums(A) + diag(A) - Bi)
# [1] -1.110223e-16  1.110223e-16