对于我的双重总结，“外部”足够快吗？

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)。在这种情况下，最佳计算仅在组合矩阵的下三角之和上进行求和，并将总和乘以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的对称情况。记得为g设置symmetric = TRUE。

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是正确的。