Search code examples
rquadratic

Efficient way of calculating quadratic forms: avoid for loops?


I want to calculate N (N is big) quadratic forms. I am using the command 'quad.form' from the R package 'emulator'. How can I implement this without using a for loop?

So far, I am using

library(emulator)

A = matrix(1,ncol=5,nrow=5) # A matrix

x = matrix(1:25,ncol=5,nrow=5) # The vectors of interest in the QF

# for loop
QF = vector()
for(i in 1:5){
QF[i] = quad.form(A,x[,i])
}

Is there a more direct and efficient way to calculate these quadratic forms?

Something intriguing is that

quad.form(A,x)

is (10 times) faster than the for loop, but I only need the diagonal of this outcome. So, it would still be an inefficient way of calculating the N quadratic forms of interest.


Solution

  • How about

    colSums(x * (A %*% x))
    

    ? Gets the right answer for this example at least ... and should be much faster!

    library("rbenchmark")
    A <- matrix(1, ncol=500, nrow=500)
    x <- matrix(1:25, ncol=500, nrow=500)
    
    library("emulator")
    aa <- function(A,x) apply(x, 2, function (y) quad.form(A,y))
    cs <- function(A,x) colSums(x * (A %*% x))
    dq <- function(A,x) diag(quad.form(A,x))
    all.equal(cs(A,x),dq(A,x))  ## TRUE
    all.equal(cs(A,x),aa(A,x))  ## TRUE
    benchmark(aa(A,x),
              cs(A,x),
              dq(A,x))
    ##       test replications elapsed relative user.self sys.self
    ## 1 aa(A, x)          100  13.121    1.346    13.085    0.024
    ## 2 cs(A, x)          100   9.746    1.000     9.521    0.224
    ## 3 dq(A, x)          100  26.369    2.706    25.773    0.592