“ * apply”家庭真的没有向量化吗？

grep("for", capture.output(getAnywhere("apply")), value = TRUE)
## [1] "        for (i in 1L:d2) {"  "    else for (i in 1L:d2) {"

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

colMeans
# function (x, na.rm = FALSE, dims = 1L) 
# {
#   if (is.data.frame(x)) 
#     x <- as.matrix(x)
#   if (!is.array(x) || length(dn <- dim(x)) < 2L) 
#     stop("'x' must be an array of at least two dimensions")
#   if (dims < 1L || dims > length(dn) - 1L) 
#     stop("invalid 'dims'")
#   n <- prod(dn[1L:dims])
#   dn <- dn[-(1L:dims)]
#   z <- 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) > 1L) {
#     dim(z) <- dn
#     dimnames(z) <- dimnames(x)[-(1L:dims)]
#   }
#   else names(z) <- dimnames(x)[[dims + 1]]
#   z
# }
# <bytecode: 0x0000000008f89d20>
#   <environment: namespace:base>

m <- 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 

#1 楼

首先，在您的示例中，您对colMeans，apply和"[.data.frame"不公平的“ data.frame”进行测试，因为它们具有开销：

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

在矩阵上，情况有些不同：

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

关于问题的主要部分，lapply / mapply / etc与简单的R循环之间的主要区别是循环的地方。正如Roland指出的那样，C和R循环都需要在每次迭代中评估R函数，这是最昂贵的。真正快速的C函数是在C中完成所有功能的函数，所以，我想，这应该是“向量化”的含义吗？

在每个“列表”元素中找到均值的示例：

（EDIT，2013年5月11日：我相信在示例中找到“均值”对于迭代评估R函数和编译后的代码之间的差异而言，这不是一个很好的设置，（1）由于R的均值算法在简单的sum(x) / length(x)上对“数字”的特殊性，（2）对它进行测试应该更有意义因此，将“均值”示例移至末尾并替换为另一个。）

作为一个简单示例，我们可以考虑找到以下每个length(x) >> lengths(x)元素的相反元素一个“列表”：

在length == 1文件中：

#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);
}

在R侧：

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

数据：

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)
})

#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

（原始均值查找示例）：

#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

#2 楼

对我来说，向量化主要是为了使您的代码更易于编写和理解。

向量化函数的目标是消除与for循环相关的簿记。例如，您可以这样写：

means <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  means[i] <- mean(mtcars[[i]])
}
sds <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  sds[i] <- sd(mtcars[[i]])
}

您可以这样写：

means <- vapply(mtcars, mean, numeric(1))
sds   <- vapply(mtcars, sd, numeric(1))

这样更容易看清什么是相同的（输入数据）和什么不同的（您要应用的功能）。

向量化的第二个优点是for循环通常用C而不是R编写。这具有显着的性能优势，但是我不认为这是矢量化的关键特性。向量化从根本上讲是在节省大脑，而不是在节省计算机工作。

#3 楼

我同意Patrick Burns的观点，即循环隐藏而不是代码矢量化。原因如下：
考虑以下C代码段：

for (int i=0; i<n; i++)
  c[i] = a[i] + b[i]

我们想要做的很清楚。但是，如何执行任务或如何执行任务并不是真的。默认情况下，for循环是串行构造。它不会告知是否或如何并行执行操作。
最明显的方法是代码以顺序方式运行。将a[i]和b[i]加载到寄存器中，将它们相加，然后将结果存储在c[i]中，然后对每个i执行此操作。
但是，现代处理器具有向量或SIMD指令集，能够在运算过程中对数据向量进行运算执行相同操作时（例如，将两个向量相加如上所述）时使用相同的指令。根据处理器/架构的不同，有可能在同一条指令下从a和b中添加四个数字，而不是一次添加一个。

我们想利用单一指令多个数据并执行数据级别的并行性，例如，一次加载4个事物，一次添加4个事物，一次存储4个事物。这就是代码矢量化。
请注意，这与代码并行化不同，后者可以同时执行多个计算。

如果编译器能够自动识别出此类代码块，那就太好了将它们矢量化，这是一项艰巨的任务。自动代码矢量化是计算机科学中一个具有挑战性的研究主题。但是随着时间的流逝，编译器在此方面做得越来越好。您可以在此处检查GNU-gcc的自动矢量化功能。同样适用于此处的LLVM-clang。而且，您还可以在最后一个链接中找到一些与gcc和ICC（英特尔C ++编译器）相比的基准测试。
例如，gcc（我在v4.9上）不会在-O2级别优化时自动对代码进行矢量化。因此，如果我们要执行上面显示的代码，它将按顺序运行。这是添加两个长度为5亿的整数矢量的时间。
我们需要添加标志-ftree-vectorize或将优化更改为级别-O3。 （请注意，-O3还会执行其他其他优化）。标志-fopt-info-vec非常有用，因为它可以告知何时成功对循环进行矢量化处理。

# compiling with -O2, -ftree-vectorize and  -fopt-info-vec
# test.c:32:5: note: loop vectorized
# test.c:32:5: note: loop versioned for vectorization because of possible aliasing
# test.c:32:5: note: loop peeled for vectorization to enhance alignment    

这告诉我们函数已被矢量化。以下是在长度为5亿的整数向量上比较非向量化版本和向量化版本的时间：

x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector

# non-vectorised, -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   1.830   0.009   1.852

# vectorised using flags shown above at -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   0.361   0.001   0.362

# both results are checked for identicalness, returns TRUE

可以安全地跳过此部分而不会失去连续性。
编译器不会总是有足够的信息可以向量化。我们可以使用OpenMP规范进行并行编程，该规范还提供了simd编译器指令，以指示编译器对代码进行向量化。重要的是要确保没有内存重叠，竞争条件等。手动对代码进行矢量化处理时，否则会导致错误的结果。

#pragma omp simd
for (i=0; i<n; i++) 
  c[i] = a[i] + b[i]

这样做，我们特别要求编译器执行以下操作：不管它如何矢量化。我们需要使用编译时标记-fopenmp激活OpenMP扩展。通过这样做：

# timing with -O2 + OpenMP with simd
x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector
system.time(.Call("Cvecsum", x, y, z))
#    user  system elapsed 
#   0.360   0.001   0.360

太好了！已使用gcc v6.2.0和llvm clang v3.9.0（均通过Homebrew，MacOS 10.12.3进行安装）进行了测试，两者均支持OpenMP 4.0。

从这个意义上讲，即使Wikipedia页面位于数组编程提到在整个数组上运行的语言通常将其称为矢量化操作，它实际上是对IMO进行循环隐藏（除非它实际上是矢量化的）。
对于R，即使C中的rowSums()或colSums()代码也不会利用代码矢量化IIUC；它只是C中的一个循环。lapply()也是如此。对于apply()，它在R中。因此，所有这些都隐藏在循环中。

简而言之，通过以下方式包装R函数：
只需在C中编写for循环！=向量化您的代码。只需在R中编写for循环！=向量化代码即可。
例如，英特尔数学内核库（MKL）实现向量形式的函数。

HTH

参考文献：


英特尔James Reinders的谈话（此答案主要是试图总结这个出色的演讲）

#4 楼

因此，将出色的答案/评论总结为一个通用的答案并提供一些背景知识：R具有4种类型的循环（从非向量化到向量化顺序）


R for循环在每次迭代中反复调用R函数（未向量化）
C循环，在每次迭代中重复调用R函数（未向量化）
仅调用R函数一次的C循环（某种程度上被向量化）
A普通的C循环，它根本不调用任何R函数，而是使用它自己的编译函数（矢量化）

，因此*apply系列是第二种类型。除了apply以外，它更是第一类

/ * .Internal（lapply（X，FUN） ））* /

/ *这是一个特殊的.internal，因此未评估的参数也是如此。它是从封闭包装中调用的，因此X和FUN是应许。 FUN必须
未经评估才能用于例如bquote。
* /


这意味着lapply的C代码接受R中的未求值函数，然后在C代码本身中对其求值。这基本上是lapply的.Internal调用之间的区别

.Internal(lapply(X, FUN))

其中有一个FUN参数包含一个R函数

和colMeans q 4312079q调用它没有.Internal参数

.Internal(colMeans(Re(x), n, prod(dn), na.rm))

FUN，不同于colMeans确切知道它需要使用什么函数，因此它在C代码内部计算均值。

您可以在lapply C代码内的每次迭代中清楚地看到R函数的求值过程。

 for(R_xlen_t i = 0; i < n; i++) {
      if (realIndx) REAL(ind)[0] = (double)(i + 1);
      else INTEGER(ind)[0] = (int)(i + 1);
      tmp = eval(R_fcall, rho);   // <----------------------------- here it is
      if (MAYBE_REFERENCED(tmp)) tmp = lazy_duplicate(tmp);
      SET_VECTOR_ELT(ans, i, tmp);
   }

总结起来， lapply没有向量化，尽管它比普通的R lapply循环有两个可能的优点。



在C中（在q4312079中，对函数进行访问），在循环中进行访问和分配似乎更快。尽管两者之间的差异似乎很大，但我们仍停留在微秒级，代价是每次迭代中R函数的估值都很高。一个简单的例子：

ffR = function(x)  {
    ans = vector("list", length(x))
    for(i in seq_along(x)) ans[[i]] = x[[i]]
    ans 
}

ffC = inline::cfunction(sig = c(R_x = "data.frame"), body = '
    SEXP ans;
    PROTECT(ans = allocVector(VECSXP, LENGTH(R_x)));
    for(int i = 0; i < LENGTH(R_x); i++) 
           SET_VECTOR_ELT(ans, i, VECTOR_ELT(R_x, i));
    UNPROTECT(1);
    return(ans); 
')

set.seed(007) 
myls = replicate(1e3, runif(1e3), simplify = FALSE)     
mydf = as.data.frame(myls)

all.equal(ffR(myls), ffC(myls))
#[1] TRUE 
all.equal(ffR(mydf), ffC(mydf))
#[1] TRUE

microbenchmark::microbenchmark(ffR(myls), ffC(myls), 
                               ffR(mydf), ffC(mydf),
                               times = 30)
#Unit: microseconds
#      expr       min        lq    median        uq       max neval
# ffR(myls)  3933.764  3975.076  4073.540  5121.045 32956.580    30
# ffC(myls)    12.553    12.934    16.695    18.210    19.481    30
# ffR(mydf) 14799.340 15095.677 15661.889 16129.689 18439.908    30
# ffC(mydf)    12.599    13.068    15.835    18.402    20.509    30

@Roland提到，它运行一个已编译的C循环，而不是一个解释的R循环


尽管在对代码进行矢量化处理时，仍需要考虑一些事项。


如果您的数据集（称为for）属于lapply类，则某些矢量化函数（例如df，data.frame，colMeans等）将不得不将其转换为矩阵首先，仅因为这是它们的设计方式。这意味着对于较大的Q4312079q，这可能会产生巨大的开销。虽然colSums不必这样做，因为它会从rowSums中提取实际的向量（因为df只是向量列表），因此，如果您没有那么多的列但有很多行，那么lapply有时会比df更好。
要记住的另一件事是R具有多种不同的函数类型，例如data.frame和通用（lapply(df, mean)，colMeans(df)），请参见此处以获取一些其他信息。泛型函数必须执行方法分派，有时这是一项昂贵的操作。例如，.Primitive是通用S3函数，而S4是mean。因此，出于上述原因，有时候S3可能比sum效率更高