Cosine similarity of one vector with many I'm keen to hear ideas for optimising R code to compute the cosine similarity of a vector x (with length l) with n other vectors (stored in any structure such as a matrix m with n rows and l columns).

Values for n will typically be much larger than values for l.

I'm currently using this custom Rcpp function to compute the similarity of a vector x to each row of a matrix m:

library(Rcpp)
cppFunction('NumericVector cosine_x_to_m(NumericVector x, NumericMatrix m) {
  int nrows = m.nrow();
  NumericVector out(nrows);
  for (int i = 0; i &lt; nrows; i++) {
    NumericVector y = m(i, _);
    out[i] = sum(x * y) / sqrt(sum(pow(x, 2.0)) * sum(pow(y, 2.0)));
  }
  return out;
}')


Varying n and l, I'm getting the following sorts of timings:



Reproducible code below.



# Function to simulate data
sim_data &lt;- function(l, n) {
  # Feature vector to be used for computing similarity
  x &lt;- runif(l)

  # Matrix of feature vectors (1 per row) to compare against x
  m &lt;- matrix(runif(n * l), nrow = n)

  list(x = x, m = m)
}

# Rcpp function to compute similarity of x to each row of m
library(Rcpp)
cppFunction('NumericVector cosine_x_to_m(NumericVector x, NumericMatrix m) {
  int nrows = m.nrow();
  NumericVector out(nrows);
  for (int i = 0; i &lt; nrows; i++) {
    NumericVector y = m(i, _);
    out[i] = sum(x * y) / sqrt(sum(pow(x, 2.0)) * sum(pow(y, 2.0)));
  }
  return out;
}')    

# Timer function
library(microbenchmark)
timer &lt;- function(l, n) {
  dat &lt;- sim_data(l, n)
  microbenchmark(cosine_x_to_m(dat$x, dat$m))
}

# Results for grid of l and n
library(tidyverse)
results &lt;- cross_d(list(l = seq(200, 1000, by = 200), n = seq(500, 4000, by = 500))) %&gt;% 
  mutate(timings = map2(l, n, timer))

# Plot results
results_plot &lt;- results %&gt;%
  unnest(timings) %&gt;% 
  mutate(time = time / 1000000) %&gt;%  # Convert time to seconds
  group_by(l, n) %&gt;% 
  summarise(mean = mean(time), ci = 1.96 * sd(time) / sqrt(n()))

pd &lt;- position_dodge(width = 20)

results_plot %&gt;% 
  ggplot(aes(n, mean, group= l)) +
  geom_line(aes(color = factor(l)), position = pd, size = 2) +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), position = pd, width = 100) +
  geom_point(position = pd, size = 2) +
  scale_color_brewer(palette = "Blues") +
  theme_minimal() +
  labs(x = "n", y = "Seconds", color = "l") +
  ggtitle("Algorithm Runtime",
          subtitle = "Error bars represent 95% confidence intervals")
I'm using Microsoft R (with Intel MKL) which makes matrix multiplications faster, but for fair comparison I set it to be single threaded.

setMKLthreads(1)


In my tests this pure R version cosine_x_to_m is twice faster than yours.

cosine_x_to_m2 = function(x,m){
  x = x / sqrt(crossprod(x));
  return(  as.vector((m %*% x) / sqrt(rowSums(m^2))) );
}


Rewriting rowSums(m^2) in C/C++ makes it even faster, about four times faster than the original.

library(ramwas)
cosine_x_to_m2 = function(x,m){
  x = x / sqrt(crossprod(x));
  return(  as.vector((m %*% x) / sqrt(rowSumsSq(m))) );
}


Initial performance:



Final version performance: