Increase performance in an R function using a vectorized approach Goal

Write an R function to generate probabilities modeled by the following equation (IRT; 2 Parameter Logistic Model):


Data

set.seed(1)
a = runif(10, 1.5, 3)
b = rnorm(10, 0, 1)
theta = rnorm(10000000)
# the output of the implementation should result in 
# a matrix with 100 mil. rows and 10 columns


Implementation

Strategy A

# the implementation of the equation
computeProbability = function(theta, b, a) {
    x = exp(1) ^ (a * (theta - b))
    return(x / (1 + x))
}

strategy.A = function(theta, b, a) {
    n.rows = length(theta)
    n.cols = length(b)

    prob_mtx = matrix(nrow = n.rows, ncol = n.cols)

    for(i in 1:n.rows)
    {
        prob_mtx[i, ] = computeProbability(theta[i], b, a)
    }
    return(prob_mtx)
}


Strategy B

strategy.B = function(theta, b, a) {
   return(t(sapply(theta, computeProbability, b = b, a = a)))
}


Strategy C

strategy.C = function(theta, b, a) {
    return(1 / (1 + exp(-sweep(outer(theta, b, "-"), 2, a, "*"))))
}


Timings

    # Strategy A        |       # Strategy B        |       # Strategy C
                        |                           |
 user  system elapsed   |    user  system elapsed   |   user  system elapsed
64.76    0.27   65.08   |   82.01    0.91   82.93   |   7.81    0.64    8.46




Question: Strategy C is by far the most efficient way, but how can I make it even faster?
Your strategy C is probably already nearly optimal. Here is an Rcpp based solution to squeeze out some additional seconds

library('Rcpp')
cppFunction(includes = '#include &lt;math.h&gt;',
            code = 'Rcpp::NumericMatrix strategy_rcpp(const Rcpp::NumericVector&amp; theta,
                                                      const Rcpp::NumericVector&amp; b,
                                                      const Rcpp::NumericVector&amp; a) {
  const int n = theta.size();
  const int m = a.size();
  if (m != b.size())
      Rcpp::stop("a and b must have equal length");
  const double e = exp(1.0);
  Rcpp::NumericMatrix ret(n, m);
  for (int j = 0; j &lt; m; ++j) {
      for (int i = 0; i &lt; n; ++i) {
         ret(i, j) = 1.0 / (1.0 + pow(e, -a[j] * (theta[i] - b[j])));
      }
  }
  return ret;
}')


Benchmarks with your data:

library('microbenchmark')
microbenchmark(times = 1,
               strategy.C(theta, b, a),
               strategy_rcpp(theta, b, a)
)

Unit: seconds                                                                
                       expr       min        lq      mean    median        uq
    strategy.C(theta, b, a) 10.996822 10.996822 10.996822 10.996822 10.996822
 strategy_rcpp(theta, b, a)  6.921705  6.921705  6.921705  6.921705  6.921705