arrays - 通过对 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 上的内部循环可以轻松地用求和替换，但外部循环......在 this question 中，建议的解决方案是sapply，但这意味着必须对所有元素进行求和。

这是我在大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

更新:

事实上，根据您的 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