如何向量化这个数学公式？

Question

我正在评估this formula（来自Székely and Rizzo, 2013，第1263 页），其中a 和b 是对称的n-by-n 矩阵。

如何使用嵌套的for 循环对该部分进行矢量化？我觉得可能有一些聪明的矩阵操作技巧可以在这里使用，也许与sweep 函数一起使用。

这是我现在拥有的：

u_squared <- function(a, b, n) {
  a. <- .colSums(a, n, n)
  b. <- .colSums(b, n, n)
  a.. <- sum(a.)
  b.. <- sum(b.)
  u1 <- 0
  u2 <- 0
  u3 <- 0
  for (i in 1:n) {
    for (j in 1:n) {
      u1 <- u1 + a[i, j] * b[i, j]
      u2 <- u2 + a[i, j] * (b.. - 2*b.[j] - 2*b.[i] + 2*b[i, j])
      u3 <- u3 + a[i, j] * (b.[i] - b[i, j])
    }
  }
  u1 <- u1 / (n * (n-1))
  u2 <- u2 / (n * (n-1) * (n-2) * (n-3))
  u3 <- u3 / (n * (n-1) * (n-2))
  return (u1 + u2 - 2 * u3)
}

例如，我知道u1 可以简单地计算为u1 <- a * b，但我想展开整个公式只是为了演示基础数学。我正在寻找u2 和u3 的类似矢量化。

Answer 1

可以很容易地向量化

u_squared1 <- function(a, b, n){
  b. <- matrix(.colSums(b, n, n), n, n) 
  B1 <- b.. - 2 * b. - 2 * t(b.) + 2*b
  B2 <- b. - b
  u1 <- sum(a*b) / (n * (n-1))
  u2 <- sum(a*B1) / (n * (n-1) * (n-2) * (n-3))
  u3 <- sum(a*B2) / (n * (n-1) * (n-2))
  return (u1 + u2 - 2 * u3)
}