My question is yet another "Vectorize this!". Similar question appeared elsewhere (Efficient way of calculating quadratic forms: avoid for loops?), but somehow I can't seem to make it work for my case.

I want to calculate quadratic form x'Sxfor every p-dimensional observation x in the sample of size n. I couldn't figure out a nice, vectorized code, so my final resort is to do it by for loop. Following toy example is for p=2, n=100.

set.seed(123)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x <- cbind(x1,x2)
Sigma <- matrix(c(1, 2, 3, 4), ncol = 2)
z  <- rep(0, n)
for (i in 1:n) {
   z[i]  <- x[i, ] %*% solve(Sigma, x[i, ]) #quadratic form of x'S^{-1}x
}

Like many other R users who worship vectorized codes, the use of for loop caused emotional pain. In order to ease the pain, I modified my code using a couple of common vectorization technique.

ap <- function(Sigma, x) apply(x, 1, function(x) x %*% solve(Sigma, x))
lap <- function(Sigma, x) unlist(lapply(1:n, function(i) x[i, ] %*% solve(Sigma, x[i, ])))
loop <- function(Sigma, x){
  z  <- rep(0, n)
  for (i in 1:n) {
    z[i]  <- x[i, ] %*% solve(Sigma, x[i, ])
  }
  z
}

But the speed comparison shows nothing much is gained.

library(microbenchmark)
microbenchmark(lap(Sigma, x), ap(Sigma, x), loop(Sigma, x))

# Unit: milliseconds
#           expr      min       lq     mean   median       uq       max neval
#  lap(Sigma, x) 4.207434 4.444895 5.092389 4.616912 5.283504  8.440802   100
#   ap(Sigma, x) 4.360204 4.523306 5.317304 4.685396 5.412771 10.168674   100
# loop(Sigma, x) 4.518645 4.679317 6.204626 4.827831 5.438908 94.115144   100

Is there any room for improvement, or should I go to Rcpp for freeing myself from sin of using for loops?

If you store the rows of x in a list and use vapply instead of lapply, you can speed this up a bit as follows

# First, make a list of the rows of x
xl <- vector("list",nrow(x))
for (i in seq_along(xl)) xl[[i]] <- x[i, ] 

# Apply solve
solve.mat <- vapply(xl, solve, numeric(2), a = Sigma)
# Take the dot product of each pair of elements
result <- colSums(solve.mat * t(x))
all(result == lap(Sigma, x))
# [1] TRUE

Writing it in one step and comparing

library(microbenchmark)
microbenchmark(lap = lap(Sigma, x),
    csums = colSums(vapply(xl, solve, numeric(2), a = Sigma) * t(x)))
# Unit: milliseconds
#   expr      min       lq     mean   median       uq      max neval
#    lap 3.013343 3.050855 3.164558 3.097901 3.136355 4.206923   100
#  csums 2.224350 2.263772 2.354349 2.289751 2.317672 3.660294   100