Increase performance in an R function using a vectorized approach
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