CODEWITHSUNDEEP
rmathmatrixvectorization
shadowtalker
shadowtalker

Reputation: 13903

How to vectorize this mathematical formula?

I'm evaluating this formula (from Székely and Rizzo, 2013, page 1263), in which a and b are symmetric, n-by-n matrices.

How can I vectorize the section with the nested for loop? I feel like there are probably clever matrix manipulation tricks that I could use here, maybe along with the sweep function.

Here is what I have now:

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)
}

For instance, I know that u1 could be computed simply as u1 <- a * b but I wanted to unroll the entire formula just to demonstrate the underlying math. What I'm looking for are similar vectorizations for u2 and u3.

Upvotes: 0

Views: 274

Answers (1)

Khashaa
Khashaa

Reputation: 7373

It can be easily vectorized

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)
}

Upvotes: 4

Related Questions