title: Is the "*apply" family really not vectorized?
body: <p>So we are used to say to every R new user that "<em><code>apply</code> isn't vectorized, check out the Patrick Burns <a href="http://www.burns-stat.com/pages/Tutor/R_inferno.pdf">R Inferno</a> Circle 4</em>" which says (I quote):</p>
  
  <blockquote>
    <p>A common reflex is to use a function in the apply family. <strong>This is not</strong>
    <strong>vectorization, it is loop-hiding</strong>. The apply function has a for loop in
    its definition. The lapply function buries the loop, but execution
    times tend to be roughly equal to an explicit for loop.</p>
  </blockquote>
  
  <p>Indeed, a quick look on the <code>apply</code> source code reveals the loop:</p>
  
  <pre><code>grep("for", capture.output(getAnywhere("apply")), value = TRUE)
  ## [1] "        for (i in 1L:d2) {"  "    else for (i in 1L:d2) {"
  </code></pre>
  
  <p>Ok so far, but a look at <code>lapply</code> or <code>vapply</code> actually reveals a completely different picture:</p>
  
  <pre><code>lapply
  ## function (X, FUN, ...) 
  ## {
  ##     FUN &lt;- match.fun(FUN)
  ##     if (!is.vector(X) || is.object(X)) 
  ##        X &lt;- as.list(X)
  ##     .Internal(lapply(X, FUN))
  ## }
  ## &lt;bytecode: 0x000000000284b618&gt;
  ## &lt;environment: namespace:base&gt;
  </code></pre>
  
  <p>So apparently there is no R <code>for</code> loop hiding there, rather they are calling internal C written function. </p>
  
  <p>A quick look in the <a href="https://github.com/wch/r-source/blob/trunk/src/main/names.c#L627-L630">rabbit</a> <a href="https://github.com/wch/r-source/blob/trunk/src/main/apply.c#L34">hole</a> reveals pretty much the same picture</p>
  
  <p>Moreover, let's take the <code>colMeans</code> function for example, which was never accused in not being vectorised</p>
  
  <pre><code>colMeans
  # function (x, na.rm = FALSE, dims = 1L) 
  # {
  #   if (is.data.frame(x)) 
  #     x &lt;- as.matrix(x)
  #   if (!is.array(x) || length(dn &lt;- dim(x)) &lt; 2L) 
  #     stop("'x' must be an array of at least two dimensions")
  #   if (dims &lt; 1L || dims &gt; length(dn) - 1L) 
  #     stop("invalid 'dims'")
  #   n &lt;- prod(dn[1L:dims])
  #   dn &lt;- dn[-(1L:dims)]
  #   z &lt;- if (is.complex(x)) 
  #     .Internal(colMeans(Re(x), n, prod(dn), na.rm)) + (0+1i) * 
  #     .Internal(colMeans(Im(x), n, prod(dn), na.rm))
  #   else .Internal(colMeans(x, n, prod(dn), na.rm))
  #   if (length(dn) &gt; 1L) {
  #     dim(z) &lt;- dn
  #     dimnames(z) &lt;- dimnames(x)[-(1L:dims)]
  #   }
  #   else names(z) &lt;- dimnames(x)[[dims + 1]]
  #   z
  # }
  # &lt;bytecode: 0x0000000008f89d20&gt;
  #   &lt;environment: namespace:base&gt;
  </code></pre>
  
  <p>Huh? It also just calls <code>.Internal(colMeans(...</code> which we can also find in the <a href="https://github.com/wch/r-source/blob/trunk/src/main/names.c#L731">rabbit hole</a>. So how is this different from <code>.Internal(lapply(..</code>?</p>
  
  <p>Actually a quick benchmark reveals that <code>sapply</code> performs no worse than <code>colMeans</code> and much better than a <code>for</code> loop for a big data set</p>
  
  <pre><code>m &lt;- as.data.frame(matrix(1:1e7, ncol = 1e5))
  system.time(colMeans(m))
  # user  system elapsed 
  # 1.69    0.03    1.73 
  system.time(sapply(m, mean))
  # user  system elapsed 
  # 1.50    0.03    1.60 
  system.time(apply(m, 2, mean))
  # user  system elapsed 
  # 3.84    0.03    3.90 
  system.time(for(i in 1:ncol(m)) mean(m[, i]))
  # user  system elapsed 
  # 13.78    0.01   13.93 
  </code></pre>
  
  <p>In other words, is it correct to say that <code>lapply</code> and <code>vapply</code> <strong>are actually vectorised</strong> (compared to <code>apply</code> which is a <code>for</code> loop that also calls <code>lapply</code>) and what did Patrick Burns really mean to say?</p>