`outer` 对我的双重求和是否足够快？

Question

我正在尝试评估 r 中的以下双倍和：

我知道outer 是一种快速的方法。我已经尝试了以下

sum(f(outer(X,X,function(x,y) (x-y)/c)))

虽然它似乎有效，但我不确定它与某些替代方案相比有多快？首先执行outer 然后执行我的功能是否会在速度方面有所不同，反之亦然？有更好的方法吗？

Answer 1

我想首先指出您可以将代码编写为

sum(f(outer(x, x, "-") / c))

这减少了函数调用开销，因为 R 中的减法已经是一个函数。试试"-"(5, 2)。

---

outer 对于您的应用程序来说已经足够快了。唯一次优的情况是您的函数f 围绕 0 对称，即f(-u) = f(u)。在这种情况下，最优计算仅对组合矩阵outer(x, x, "-") 的下三角求和，并将和乘以 2 以对非对角线求和。最后，添加对角线结果。

以下函数执行此操作。我们为组合矩阵的下三角部分（不包括对角线）生成(i, j) 索引，那么outer(x, x, "-") / c 的下三角部分将是dx <- (x[i] - x[j]) / c。现在，

• 如果f是对称的，结果是2 * sum(f(dx)) + n * f(0)，这比outer快；
• 如果f 是不对称的，我们必须做sum(f(dx)) + sum(f(-dx)) + n * f(0)，这不会比outer 有任何优势。

---

## `x` is the vector, `f` is your function of interest, `c` is a constant
## `symmetric` is a switch; only set `TRUE` when `f` is symmetric around 0
g <- function (x, f, c, symmetric = FALSE) {
  n <- length(x)
  j <- rep.int(1:(n-1), (n-1):1)
  i <- sequence((n-1):1) + j
  dx <- (x[i] - x[j]) / c
  if (symmetric) 2 * sum(f(dx)) + n * f(0)
  else sum(f(dx)) + sum(f(-dx)) + n * f(0)
  }

---

在这里考虑一个小例子。让我们假设c = 2 和一个向量x <- 1:500。我们还考虑了一个对称函数f1 <- cos 和一个非对称函数f2 <- sin。让我们做一个基准测试：

x <- 1:500
library(microbenchmark)

我们首先考虑f1 的对称情况。记得将symmetric = TRUE 设置为g。

microbenchmark(sum(f1(outer(x,x,"-")/2)), g(x, f1, 2, TRUE))

#Unit: milliseconds
#                        expr      min       lq     mean   median       uq
# sum(f2(outer(x, x, "-")/2)) 32.79472 35.35316 46.91560 36.78152 37.63580
#           g(x, f2, 2, TRUE) 20.24940 23.34324 29.97313 24.45638 25.33352
#      max neval cld
# 133.5494   100   b
# 120.3278   100  a 

在这里我们看到g 更快。

现在考虑f2 的不对称情况。

microbenchmark(sum(f2(outer(x,x,"-")/2)), g(x, f2, 2))

#Unit: milliseconds
#                        expr      min       lq     mean   median       uq
# sum(f2(outer(x, x, "-")/2)) 32.84412 35.55520 44.33684 36.95336 37.89508
#                 g(x, f2, 2) 36.71572 39.11832 50.54516 40.25590 41.75060
#      max neval cld
# 134.2991   100   a
# 142.5143   100   a

果然，这里没有优势。

---

是的，我们还想检查g 是否进行了正确的计算。考虑一个小例子就足够了，x <- 1:5。

x <- 1:5

#### symmetric case ####

sum(f1(outer(x, x, "-") / 2))
# [1] 14.71313

g(x, f1, 2, TRUE)
# [1] 14.71313

#### asymmetric case ####

sum(f2(outer(x, x, "-") / 2))
# [1] 0

g(x, f2, 2)
# [1] 0

所以g 是正确的。