CODEWITHSUNDEEP
arraysrloopsmatrixvectorization
Andreï V. Kostyrka
Andreï V. Kostyrka

Reputation: 352

Vectorising (or speeding up) a double loop with summation over non-identical indices in R

I am trying to optimise the code designed to compute double sums of product of the elements of two square matrices. Let’s say we have two square matrices of size n, W and V. The object that needs to be computed is a vector B with elements

In simple terms: compute element-by-element products of two different rows in two different matrices and take their sum, then take an extra sum over all rows of the second matrix (sans identical indices).

The problem is, the computational complexity of this task seemingly O(n3) because the length of this object we are creating, B, is n, and each element requires two summations. This is what I have come up with:

  1. For given i and j (i≠j), start with the inner sum over k. Sum for all k, then subtract the terms for k=i and k=j, and multiply by the indicator of j≠i.
  2. Since the restriction j≠i has been taken care of in the inner sum, the outer sum is taken just for j=1,...,n.

If we denote , then the two steps will look like and .

However, writing a loop turned out to be very inefficient. n=100 works quickly (0.05 seconds). But, for instance, when n=500 (we are talking about real-world applications here), the average computation time is 3 seconds, and for n=1000, it jumps to 22 s.

The inner loop over k can be easily replaced by a sum, but the outer one... In this question, the suggested solution is sapply, but it implies that the summation must be done over all elements.

This is the code I am trying to evaluate before the heat death of the Universe for large n. 

set.seed(1)
N <- 500
x1 <- rnorm(N)
x2 <- rchisq(N, df=3)

bw1 <- bw.nrd(x1)
bw2 <- bw.nrd(x2)
w <- outer(x1, x1, function(x, y) dnorm((x-y)/bw1) )
w <- w/rowSums(w)

v <- outer(x2, x2, function(x, y) dnorm((x-y)/bw2) )
v <- v/rowSums(v)

Bij <- matrix(NA, ncol=N, nrow=N)
for (i in 1:N) { # Around 22 secs for N=1000
  for (j in 1:N) {
    Bij[i, j] <- (sum(w[i, ]*v[j, ]) - w[i, i]*v[j, i] - w[i, j]*v[j, j]) * (i!=j)
  }
}
Bi <- rowSums(Bij)

How would an expert R programmer vectorise such kind of loops?

Upvotes: 4

Views: 296

Answers (2)

Julius Vainora
Julius Vainora

Reputation: 48191

Update:

In fact, given your expression for B_{ij}, we may also do the following

diag(w) <- diag(v) <- 0
BBij <-  tcrossprod(w, v)
diag(BBij) <- 0

range(rowSums(BBij) - Bi)
# [1] -2.220446e-16  0.000000e+00
range(BBij - Bij)
# [1] -6.938894e-18  5.204170e-18

Hence, while somewhat obvious, it may also be an interesting observation for your purposes that neither B_{ij} nor B_i depend on the diagonals of W and V.

Initial answer:

Since enter image description here where the diagonals of W and V can be set to zero and V_{.k} denotes the sum of the k-th column of V, we have

diag(w) <- diag(v) <- 0 
A <- w * v
rowSums(sweep(w, 2, colSums(v), `*`)) - rowSums(A) + diag(A)

where

range(rowSums(sweep(w, 2, colSums(v), `*`)) - rowSums(A) + diag(A) - Bi)
# [1] -1.110223e-16  1.110223e-16

Upvotes: 4

RolandASc
RolandASc

Reputation: 3923

Without looking into the content of your matrices w and v, your double for-loop can be replaced with simple matrix operations, using one matrix multiplication (tcrossprod), transpose (t) and diagonal extraction:

Mat.ij <- tcrossprod(w, v) - 
    matrix(rep(diag(w), times = N), nrow = N) * t(v) - 
    w * matrix(rep(diag(v), each = N), nrow = N)
diag(Mat.ij) <- 0

all.equal(Bij, Mat.ij)
[1] TRUE

Upvotes: 3

Related Questions