"*apply" 系列真的没有向量化吗？

首先，在您的示例中，您对“data.frame”进行测试，这对 colMeans、apply 和 "[.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 2016 年 5 月 11 日：我认为找到“平均值”的示例对于迭代评估 R 函数和编译代码之间的差异不是一个好的设置，(1) 因为R 的均值算法在简单 sum(x) / length(x) 上的“数字”上的特殊性和（2）用 length(x) >> lengths(x) 对“列表”进行测试应该更有意义。因此，“均值”示例移至末尾并替换为另一个。）

作为一个简单的例子，我们可以考虑找到“列表”的每个 length == 1 元素的对立面：

在 tmp.c 文件中：

#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

---

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

矢量化函数的目标是消除与 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 语言编写的。这具有显着的性能优势，但我不认为它是矢量化的关键属性。矢量化从根本上说是为了节省您的大脑，而不是节省计算机工作。

---

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

考虑这个 C 代码段：

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

我们想做的事情已经很清楚了。但如何执行任务或如何执行它并不是真正的。默认情况下，for 循环是一个串行结构。它没有告知是否或如何可以并行完成。

最明显的方式是代码以顺序方式运行。将 a[i] 和 b[i] 加载到寄存器，添加它们，将结果存储在 c[i] 中，并对每个 i 执行此操作。

但是，现代处理器具有 vector or SIMD 指令集，当执行相同的操作（例如，将两个向量相加为如上所示）。根据处理器/架构，可能会在同一指令下添加来自 a 和 b 的四个数字，而不是一次添加一个。

我们想利用单指令多数据并执行数据级并行性，例如，一次加载 4 个内容，一次添加 4 个内容，一次存储 4 个内容。这就是代码向量化。请注意，这与代码并行化不同——其中多个计算同时执行。

如果编译器能识别出这样的代码块并自动将它们向量化，那就太好了，这是一项艰巨的任务。 Automatic code vectorisation 是计算机科学中一个具有挑战性的研究课题。但是随着时间的推移，编译器已经变得更好了。您可以检查 GNU-gcc here 的 自动矢量化 功能。 LLVM-clang here 也是如此。您还可以在最后一个链接中找到与 gcc 和 ICC（英特尔 C++ 编译器）进行比较的一些基准测试。

例如，gcc（我在 v4.9）不会在 -O2 级别优化中自动矢量化代码。因此，如果我们要执行上面显示的代码，它将按顺序运行。这是添加两个长度为 5 亿的整数向量的时间。

我们需要添加标志 -ftree-vectorize 或将优化更改为级别 -O3。 （请注意，-O3 也执行 other additional optimisations）。标志 -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 specification for parallel programming，它还提供了一个 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（均通过自制软件安装，MacOS 10.12.3）进行测试的，它们都支持 OpenMP 4.0。

从这个意义上说，尽管 Wikipedia page on Array Programming 提到对整个数组进行操作的语言通常将其称为 矢量化操作，但它确实是 循环隐藏 IMO（除非它实际上是矢量化的）。

在 R 的情况下，即使是 C 中的 rowSums() 或 colSums() 代码也不会利用 代码矢量化 IIUC；它只是 C 中的一个循环。lapply() 也是如此。在 apply() 的情况下，它在 R 中。因此，所有这些都是循环隐藏。

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

高温高压

参考：

英特尔 James Reinders 的演讲（这个答案主要是试图总结这个精彩的演讲）

---

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

for 循环在每次迭代中重复调用 R 函数（未向量化） C 循环在每次迭代中重复调用 R 函数（未向量化） C 循环仅调用 R 函数一次（有点向量化）一个普通的 C 循环不调用任何 R 函数并使用它自己的编译函数（矢量化）

所以 *apply 系列是第二种类型。除了 apply，它更像是第一种类型

您可以从其 source code 中的评论中理解这一点

/* .Internal(lapply(X, FUN)) */ /* 这是一个特殊的 .Internal，因此有未评估的参数。它是从一个闭包包装器中调用的，所以 X 和 FUN 是 Promise。 FUN 必须未经评估才能在例如 bquote 中使用。 */

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

.Internal(lapply(X, FUN))

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

而 没有 的 colMeans .Internal 调用具有 FUN 参数

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

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

可以清楚的看到lapply C code内每次迭代中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 for 循环有两个可能的优势

在循环中访问和分配似乎在 C 中更快（即在应用函数时）虽然差异似乎很大，但我们仍然停留在微秒级别，代价高昂的是每次迭代中 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 = 复制(1e3, runif(1e3), 简化 = 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 156661.889 16129.689 16129.689 18439.908 30＃ffc（我commant in＃ffc（my lo y lo intry tody tody tody tody tody tocry in tocry in tocrand） R循环

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

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

---