How to Know a Function or an Operation in R Is Vectorized

How do I know a function or an operation in R is vectorized?

Vectorization in R basically means that any looping is moved to a faster, compiled language such as C or FORTRAN. For that to occur the vector(s) in question should be "atomic" - i.e. it should be "flat" and homogeneous - and the vector type, which you can check with typeof(), should make sense for the operation(s) being performed. If it is atomic then it is vectorized.

You can check if a vector is atomic using is.atomic(). Another type of vector that is not vectorized is called "recursive", which you can check using is.recursive(). Recursive objects can contain other objects of any type, i.e. they can be heterogeneous. Lists and data frames are recursive.

Try something like the following to gain some insight into atomic vs. recursive:

# Atomic:
1
1:3
c("a", "b", "c")
c(T, F, T)

# Recursive:
list(nums = 1:3, letts = c("a", "b", "c"), logics = c(T, F, T))
data.frame(nums = 1:3, letts = c("a", "b", "c"), logics = c(T, F, T))

# Vectors can be atomic or recursive:
is.vector(1:9) # TRUE
is.atomic(1:9) # TRUE
is.recursive(1:9) # FALSE

is.vector(list(nums = 1:9, chars = "x")) # TRUE
is.atomic(list(1:9)) # FALSE
is.recursive(list(1:9)) # TRUE

# Matrices are atomic, data frames are recursive:
is.vector(matrix(1:9, 3)) # FALSE
is.atomic(matrix(1:9, 3)) # TRUE
is.recursive(matrix(1:9, 3)) # FALSE

is.vector(as.data.frame(matrix(1:9, 3))) # FALSE
is.atomic(as.data.frame(matrix(1:9, 3))) # FALSE
is.recursive(as.data.frame(matrix(1:9, 3))) # TRUE

I think you can assume that many, if not most, of the R functions that you use most frequently are vectorized. I don't think there is any way to check this other than by looking at the documentation or the function internals. Whenever you think about writing a for loop to do simple element-wise operations, think about how to do it using vectorization. With enough practice it will become second nature to you. For more details I can recommend this blog post from Noam Ross.

Vectorize my thinking: Vector Operations in R

Here's what seems like another very R-type way to generate the sums. Generate a vector that is as long as your input vector, containing nothing but the repeated sum of n elements. Then, subtract your original vector from the sums vector. The result: a vector (isums) where each entry is your original vector less the ith element.

> (my.data$item[my.data$fixed==0])
[1] 1 1 3 5 7
> sums <- rep(sum(my.data$item[my.data$fixed==0]),length(my.data$item[my.data$fixed==0]))
> sums
[1] 17 17 17 17 17
> isums <- sums - (my.data$item[my.data$fixed==0])
> isums
[1] 16 16 14 12 10

Vectorized IF statement in R?

x <- seq(0.1,10,0.1)

> x
  [1]  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5
 [16]  1.6  1.7  1.8  1.9  2.0  2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0
 [31]  3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9  4.0  4.1  4.2  4.3  4.4  4.5
 [46]  4.6  4.7  4.8  4.9  5.0  5.1  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9  6.0
 [61]  6.1  6.2  6.3  6.4  6.5  6.6  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4  7.5
 [76]  7.6  7.7  7.8  7.9  8.0  8.1  8.2  8.3  8.4  8.5  8.6  8.7  8.8  8.9  9.0
 [91]  9.1  9.2  9.3  9.4  9.5  9.6  9.7  9.8  9.9 10.0

> ifelse(x < 5, 1, 2)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

How to vectorize an operation on a series of vectors in R

Up front, your myv might be organized as a series of vectors, one column each; it is better for many tools to convert it into a list of vectors.

asplit(myv, 2)
# [[1]]
#  [1]  1  2  3  4  5  6  7  8  9 10
# [[2]]
#  [1] 11 12 13 14 15 16 17 18 19 20
# [[3]]
#  [1] 21 22 23 24 25 26 27 28 29 30

base R

sapply/lapply are to a single vector/list as mapply/Map are to n of them.

Map(myf, mya, asplit(myv , 2))
# [[1]]
# [1] 7.980123
# [[2]]
# [1] 17.64959
# [[3]]
# [1] 26.80944
mapply(myf, mya, asplit(myv , 2))
# [1]  7.980123 17.649590 26.809440

tidyverse

The order of arguments is different, and instead of individual arguments it needs all of them in a list itself.

purrr::pmap(list(mya, asplit(myv , 2)), myf)
# [[1]]
# [1] 7.980123
# [[2]]
# [1] 17.64959
# [[3]]
# [1] 26.80944
purrr::pmap_dbl(list(mya, asplit(myv , 2)), myf)
# [1]  7.980123 17.649590 26.809440

Alternative approach, given the comments.

This approach truly is vectorized, but has deconstructed the function a little.

colSums(t(mya / (mya / t(myv) + 1)))
# [1]  7.980123 17.649590 26.809440

To get to this point, one needs to recognize where transpose and such is necessary. I'll start with some known points:

mya[1] / myv[,1] + 1
#  [1] 2.000000 1.500000 1.333333 1.250000 1.200000 1.166667 1.142857 1.125000 1.111111 1.100000

In order to mimic that with matrices (and not just vectors), we might try

(mya / myv + 1)
#           [,1]     [,2]     [,3]
#  [1,] 2.000000 1.181818 1.142857
#  [2,] 2.000000 1.250000 1.045455
#  [3,] 2.000000 1.076923 1.086957
#  [4,] 1.250000 1.142857 1.125000
#  [5,] 1.400000 1.200000 1.040000
#  [6,] 1.500000 1.062500 1.076923
#  [7,] 1.142857 1.117647 1.111111
#  [8,] 1.250000 1.166667 1.035714
#  [9,] 1.333333 1.052632 1.068966
# [10,] 1.100000 1.100000 1.100000

But if you notice, the division of mya over myv is column-wise, so it is expanding to

c(mya[1] / myv[1,1], mya[2] / myv[2,1], mya[3] / myv[3,1], mya[1] / myv[4,1], ...)

where we would prefer it to be transposed. Okay, so we transpose it so that the rows of myv are vertical for the division.

(mya / t(myv) + 1)[1,]
#  [1] 2.000000 1.500000 1.333333 1.250000 1.200000 1.166667 1.142857 1.125000 1.111111 1.100000

That's better. Now we need to do the same for the next step. That brings us to

t(mya / (mya / t(myv) + 1))
#            [,1]     [,2]     [,3]
#  [1,] 0.5000000 1.692308 2.625000
#  [2,] 0.6666667 1.714286 2.640000
#  [3,] 0.7500000 1.733333 2.653846
#  [4,] 0.8000000 1.750000 2.666667
#  [5,] 0.8333333 1.764706 2.678571
#  [6,] 0.8571429 1.777778 2.689655
#  [7,] 0.8750000 1.789474 2.700000
#  [8,] 0.8888889 1.800000 2.709677
#  [9,] 0.9000000 1.809524 2.718750
# [10,] 0.9090909 1.818182 2.727273

Since you wanted to sum across each of the mya values. Knowing that we have three in mya and we see three columns, one might infer we need to sum each column. We can prove that empirically:

sum(mya[1] / (mya[1] / myv[,1] + 1))
# [1] 7.980123
colSums(t(mya / (mya / t(myv) + 1)))
# [1]  7.980123 17.649590 26.809440

But really, we don't need to tranpose then sum columns when we can not-transpose and sum the rows :-)

rowSums(mya / (mya / t(myv) + 1))
# [1]  7.980123 17.649590 26.809440

Using Apply or Vectorize to apply custom function to a dataframe

For the particular task requested it could be

celebrities$newcol <- with(celebrities, age + income)

The + function is inherently vectorized. Using apply with sum is inefficient. Using apply could have been greatly simplified by omitting the first column because that would avoid the coercion to a character matrix caused by the first column.

 celebrities$newcol <- apply(celebrities[-1], function(x) sum(x) )

That way you would avoid coercing the vectors to "character" and then needing to coerce back the formerly-numeric columns to numeric. Using sum inside apply does get around the fact that sum is not vectorized, but it's an example of inefficient R coding.

You get automatic vectorization if the "inner" algorithm can be constructed completely from vectorized functions: the Math and Ops groups being the usual components. See ?Ops. Otherwise, you may need to use mapply or Vectorize.

Vectorize a product calculation which depends on previous elements?

ifelse is vectorized and there's a bit of a penalty if you're using it on one element at a time in a for-loop. In your example, you can get a pretty good speedup by using if instead of ifelse.

fun1 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- ifelse(z[i-1]==1, 1, 0)
  }
  z
}

fun2 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- if(z[i-1]==1) 1 else 0
  }
  z
}

z <- c(1,1,0,0,0,0)
identical(fun1(z),fun2(z))
# [1] TRUE
system.time(replicate(10000, fun1(z)))
#   user  system elapsed 
#   1.13    0.00    1.32
system.time(replicate(10000, fun2(z)))
#   user  system elapsed 
#   0.27    0.00    0.26

You can get some additional speed gains out of fun2 by compiling it.

library(compiler)
cfun2 <- cmpfun(fun2)
system.time(replicate(10000, cfun2(z)))
#   user  system elapsed 
#   0.11    0.00    0.11

So there's a 10x speedup without vectorization. As others have said (and some have illustrated) there are ways to vectorize your example, but that may not translate to your actual problem. Hopefully this is general enough to be applicable.

The filter function may be useful to you as well if you can figure out how to express your problem in terms of a autoregressive or moving average process.

Is the *apply family really not vectorized?

First of all, in your example you make tests on a "data.frame" which is not fair for colMeans, apply and "[.data.frame" since they have an overhead:

system.time(as.matrix(m))  #called by `colMeans` and `apply`
#   user  system elapsed 
#   1.03    0.00    1.05
system.time(for(i in 1:ncol(m)) m[, i])  #in the `for` loop
#   user  system elapsed 
#  12.93    0.01   13.07

On a matrix, the picture is a bit different:

mm = as.matrix(m)
system.time(colMeans(mm))
#   user  system elapsed 
#   0.01    0.00    0.01 
system.time(apply(mm, 2, mean))
#   user  system elapsed 
#   1.48    0.03    1.53 
system.time(for(i in 1:ncol(mm)) mean(mm[, i]))
#   user  system elapsed 
#   1.22    0.00    1.21

Regading the main part of the question, the main difference between lapply/mapply/etc and straightforward R-loops is where the looping is done. As Roland notes, both C and R loops need to evaluate an R function in each iteration which is the most costly. The really fast C functions are those that do everything in C, so, I guess, this should be what "vectorised" is about?

An example where we find the mean in each of a "list"s elements:

(EDIT May 11 '16 : I believe the example with finding the "mean" is not a good setup for the differences between evaluating an R function iteratively and compiled code, (1) because of the particularity of R's mean algorithm on "numeric"s over a simple sum(x) / length(x) and (2) it should make more sense to test on "list"s with length(x) >> lengths(x). So, the "mean" example is moved to the end and replaced with another.)

As a simple example we could consider the finding of the opposite of each length == 1 element of a "list":

In a tmp.c file:

#include <R.h>
#define USE_RINTERNALS 
#include <Rinternals.h>
#include <Rdefines.h>

/* call a C function inside another */
double oppC(double x) { return(ISNAN(x) ? NA_REAL : -x); }
SEXP sapply_oppC(SEXP x)
{
    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));
    for(int i = 0; i < LENGTH(x); i++) 
        REAL(ans)[i] = oppC(REAL(VECTOR_ELT(x, i))[0]);

    UNPROTECT(1);
    return(ans);
}

/* call an R function inside a C function;
 * will be used with 'f' as a closure and as a builtin */    
SEXP sapply_oppR(SEXP x, SEXP f)
{
    SEXP call = PROTECT(allocVector(LANGSXP, 2));
    SETCAR(call, install(CHAR(STRING_ELT(f, 0))));

    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));     
    for(int i = 0; i < LENGTH(x); i++) { 
        SETCADR(call, VECTOR_ELT(x, i));
        REAL(ans)[i] = REAL(eval(call, R_GlobalEnv))[0];
    }

    UNPROTECT(2);
    return(ans);
}

And in R side:

system("R CMD SHLIB /home/~/tmp.c")
dyn.load("/home/~/tmp.so")

with data:

set.seed(007)
myls = rep_len(as.list(c(NA, runif(3))), 1e7)

#a closure wrapper of `-`
oppR = function(x) -x

for_oppR = compiler::cmpfun(function(x, f)
{
    f = match.fun(f)  
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[[i]] = f(x[[i]])
    return(ans)
})

Benchmarking:

#call a C function iteratively
system.time({ sapplyC =  .Call("sapply_oppC", myls) }) 
#   user  system elapsed 
#  0.048   0.000   0.047 

#evaluate an R closure iteratively
system.time({ sapplyRC =  .Call("sapply_oppR", myls, "oppR") }) 
#   user  system elapsed 
#  3.348   0.000   3.358 

#evaluate an R builtin iteratively
system.time({ sapplyRCprim =  .Call("sapply_oppR", myls, "-") }) 
#   user  system elapsed 
#  0.652   0.000   0.653 

#loop with a R closure
system.time({ forR = for_oppR(myls, "oppR") })
#   user  system elapsed 
#  4.396   0.000   4.409 

#loop with an R builtin
system.time({ forRprim = for_oppR(myls, "-") })
#   user  system elapsed 
#  1.908   0.000   1.913 

#for reference and testing 
system.time({ sapplyR = unlist(lapply(myls, oppR)) })
#   user  system elapsed 
#  7.080   0.068   7.170 
system.time({ sapplyRprim = unlist(lapply(myls, `-`)) }) 
#   user  system elapsed 
#  3.524   0.064   3.598 

all.equal(sapplyR, sapplyRprim)
#[1] TRUE 
all.equal(sapplyR, sapplyC)
#[1] TRUE
all.equal(sapplyR, sapplyRC)
#[1] TRUE
all.equal(sapplyR, sapplyRCprim)
#[1] TRUE
all.equal(sapplyR, forR)
#[1] TRUE
all.equal(sapplyR, forRprim)
#[1] TRUE

(Follows the original example of mean finding):

#all computations in C
all_C = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP tmp, ans;
    PROTECT(ans = allocVector(REALSXP, LENGTH(R_ls)));

    double *ptmp, *pans = REAL(ans);

    for(int i = 0; i < LENGTH(R_ls); i++) {
        pans[i] = 0.0;

        PROTECT(tmp = coerceVector(VECTOR_ELT(R_ls, i), REALSXP));
        ptmp = REAL(tmp);

        for(int j = 0; j < LENGTH(tmp); j++) pans[i] += ptmp[j];

        pans[i] /= LENGTH(tmp);

        UNPROTECT(1);
    }

    UNPROTECT(1);
    return(ans);
')

#a very simple `lapply(x, mean)`
C_and_R = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP call, ans, ret;

    PROTECT(call = allocList(2));
    SET_TYPEOF(call, LANGSXP);
    SETCAR(call, install("mean"));

    PROTECT(ans = allocVector(VECSXP, LENGTH(R_ls)));
    PROTECT(ret = allocVector(REALSXP, LENGTH(ans)));

    for(int i = 0; i < LENGTH(R_ls); i++) {
        SETCADR(call, VECTOR_ELT(R_ls, i));
        SET_VECTOR_ELT(ans, i, eval(call, R_GlobalEnv));
    }

    double *pret = REAL(ret);
    for(int i = 0; i < LENGTH(ans); i++) pret[i] = REAL(VECTOR_ELT(ans, i))[0];

    UNPROTECT(3);
    return(ret);
')                    

R_lapply = function(x) unlist(lapply(x, mean))                       

R_loop = function(x) 
{
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[i] = mean(x[[i]])
    return(ans)
} 

R_loopcmp = compiler::cmpfun(R_loop)

set.seed(007); myls = replicate(1e4, runif(1e3), simplify = FALSE)
all.equal(all_C(myls), C_and_R(myls))
#[1] TRUE
all.equal(all_C(myls), R_lapply(myls))
#[1] TRUE
all.equal(all_C(myls), R_loop(myls))
#[1] TRUE
all.equal(all_C(myls), R_loopcmp(myls))
#[1] TRUE

microbenchmark::microbenchmark(all_C(myls), 
                               C_and_R(myls), 
                               R_lapply(myls), 
                               R_loop(myls), 
                               R_loopcmp(myls), 
                               times = 15)
#Unit: milliseconds
#            expr       min        lq    median        uq      max neval
#     all_C(myls)  37.29183  38.19107  38.69359  39.58083  41.3861    15
#   C_and_R(myls) 117.21457 123.22044 124.58148 130.85513 169.6822    15
#  R_lapply(myls)  98.48009 103.80717 106.55519 109.54890 116.3150    15
#    R_loop(myls) 122.40367 130.85061 132.61378 138.53664 178.5128    15
# R_loopcmp(myls) 105.63228 111.38340 112.16781 115.68909 128.1976    15

R: Monitor progress of vectorized operation

In my comment, I asked what your code is and how you achieve vectorization. I think this matters. Although generally speaking, vectorization is achieved by using loops in compiled code, I am not entirely sure of this. Therefore, I would like to be less confident in saying "absolutely no".

However, if you want to track progress at R level, you must be able to get an index, like i used in an R level for loop. Now, let's check what most R vectorized functions look like:

> grep
function (pattern, x, ignore.case = FALSE, perl = FALSE, value = FALSE, 
    fixed = FALSE, useBytes = FALSE, invert = FALSE) 
{
    if (!is.character(x)) 
        x <- structure(as.character(x), names = names(x))
    .Internal(grep(as.character(pattern), x, ignore.case, value, 
        perl, fixed, useBytes, invert))
}
<bytecode: 0xa34dfe0>
<environment: namespace:base>

> gsub
function (pattern, replacement, x, ignore.case = FALSE, perl = FALSE, 
    fixed = FALSE, useBytes = FALSE) 
{
    if (!is.character(x)) 
        x <- as.character(x)
    .Internal(gsub(as.character(pattern), as.character(replacement), 
        x, ignore.case, perl, fixed, useBytes))
}

In above examples, we see that those vectorized R functions are merely a thin wrapper of compiled code (see the .Internal()). There are no explicit loop index for you to refer to. Hence for those example functions, tracking progress is not possible.

I suggest you have a look at the particular function you used. That is the best way to convince yourself.

follow up

Originally, I put lapply in my examples:

> lapply
function (X, FUN, ...) 
{
    FUN <- match.fun(FUN)
    if (!is.vector(X) || is.object(X)) 
        X <- as.list(X)
    .Internal(lapply(X, FUN))
}
<bytecode: 0x9c5c464>
<environment: namespace:base>

Then @RichardScriven expressed his view of *apply family. On stack overflow, these two posts/answers are extremely useful to understanding vectorization issues in R:

is the "*apply" family really not vectorized
is R's apply family more than syntactic sugar?

Truly, though lapply calls C code to do the loop, it has to evaluate R function FUN along the loop. Hence:

if FUN dominates execution time, then lapply will not have noticeable advantage over R's for loop.
if FUN does so little work, that the loop overhead dominates the execution, then lapply will have noticeable advantage over R's for loop, because for loop in C is more "light weighted".

Discussing the performance of lapply is off-topic in this post, so I will not attach examples for demonstration.