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

Question

我正在评估 this formula (from Székely and Rizzo, 2013，第 1263 页），其中 a 和 b 是对称的 n×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)
}

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

How to vectorize this mathematical formula?

math

r

matrix

vectorization