CODEWITHSUNDEEP
rperformancematrixstatisticsmatrix-multiplication
Zilch
Zilch

Reputation: 23

Efficient way of computing quadratic forms for outer product matrices

I'm currently computing a quadratic form by taking a known vector and using the element wise multiplication of two outer products as the input matrix. To be specific, my code looks something like this

set.seed(42)  # for sake of reproducibility
library(emulator) 
Fun <- function(a,b) sqrt((1/(2*pi)))*exp(-0.5*(a-b)^2)
n <- 5000
x <- rnorm(n)
y <- rnorm(n)
u <- rnorm(n)
I <- quad.form(outer(x,x,Fun)*outer(y,y,Fun),u)

This is quite slow and the problem gets significantly worse as n increases. As far as I can make out, the part that causes the problem is the outer(x,x,Fun)*outer(y,y,Fun) term inside the quadratic form.

Is there any way to speed this up?

Upvotes: 2

Views: 296

Answers (1)

Roland
Roland

Reputation: 132746

You can half the timings by making use of the symmetry. I find it easiest to quickly write an Rcpp function:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix foo(NumericVector x, NumericVector y) {
  double n = x.size();
  NumericMatrix M(n, n);
  for(double i = 0; i < n; i++) {
    for(double j = 0; j <= i; j++) {
      M(i,j) = ((1/(2*M_PI))) *
        (exp(-0.5*(pow(x(i)-x(j), 2) + pow(y(i)-y(j), 2))));
      M(j,i) = M(i,j);
    }
  }
  return M;
}

Timings:

set.seed(42)  # for sake of reproducibility
library(emulator) 
Fun <- function(a,b) sqrt((1/(2*pi)))*exp(-0.5*(a-b)^2)
n <- 1e4
x <- rnorm(n)
y <- rnorm(n)
u <- rnorm(n)

system.time({
  I <- quad.form(outer(x,x,Fun)*outer(y,y,Fun),u)
})
#       user      system     elapsed 
#      4.287       1.373       5.687 

system.time({
  J <- quad.form(foo(x, y),u)
})
#       user      system     elapsed 
#      2.232       0.168       2.409 

all.equal(I, J)
#[1] TRUE

Further improvements would require parallelization (or possibly use of some maths).

Upvotes: 3

Related Questions