“申请”家庭是否真的没有vector化？
所以我们习惯于对每个R新用户说“ apply不是vector化的，查看Patrick Burns R Inferno Circle 4 ”（我引用）：
一个常见的反应是在应用系列中使用一个函数。 这不是 vector化，而是循环隐藏 。 apply函数在其定义中有一个for循环。 lapply函数隐藏循环，但执行时间往往大致等于明确的for循环。
实际上，快速查看apply源代码揭示了循环：
grep("for", capture.output(getAnywhere("apply")), value = TRUE) ## [1] " for (i in 1L:d2) {" " else for (i in 1L:d2) {"
那么到目前为止，但看一下lapply或vapply实际上揭示了一个完全不同的画面：
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>
所以显然没有R for循环隐藏在那里，而是他们调用内部C写function。
在兔子 洞快速查看揭示几乎相同的图片
而且，我们以colMeans函数为例，这个函数从未被vector化过
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>
咦？ 它也只是叫.Internal(colMeans(... ，我们也可以在兔子洞中find。那么这与内部.Internal(lapply(.. ？
实际上，一个快速的基准testing表明， sapply并不比colMeans糟糕，并且比一个大数据集的for循环要好得多
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
换句话说，是否说实际上是向lapply而vapply 实际上是向量化的 （相对于apply这个也称为lapply循环），帕特里克·伯恩斯究竟意味着什么呢？
首先，在你的例子中，你对colMeans不公平的“data.frame”进行testing， 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
在matrix上，图片有点不同：
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-loop之间的主要区别是循环完成的地方。 正如罗兰所指出的，C和R循环都需要在每次迭代中评估R函数，这是最昂贵的。 真正快速的C函数就是那些用C来做所有事情的，所以我猜这应该是“vector化”的意思？
一个例子，我们在每个“list”元素中find均值：
（ 编辑5月11日16 ：我相信find“意思”的例子不是一个好的设置，以评估一个R函数迭代和编译代码之间的差异，（1）由于R的平均algorithm的“数字” sum(x) / length(x)和（2）在length(x) >> lengths(x) “list”上进行testing更有意义，所以“mean”示例被移动到最后，换成另一个。）
作为一个简单的例子，我们可以考虑一个“列表”的每个length == 1元素的相反的发现：
在一个tmp.c文件中：
#include <Rh> #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
对我来说，vector化主要是为了让你的代码更容易编写，更易于理解。
vector化函数的目标是消除与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))
这样可以更容易地看到相同的内容（input数据）和不同的内容（您正在应用的function）。
vector化的一个次要优点是for循环通常用C语言而不是R语言编写。这样做有很大的性能优势，但我不认为它是vector化的关键属性。 Vectorisation从根本上挽救你的大脑，而不是节省电脑的工作。
我同意帕特里克·伯恩斯（Patrick Burns）的观点，认为这是循环隐藏而不是代码vector化 。 原因如下：
考虑这个C代码片段：
for (int i=0; i<n; i++) c[i] = a[i] + b[i]
我们想做的事情很清楚。 但是如何执行任务或如何执行任务并不是真的。 一个for-loop默认是一个串行结构。 它并不告知是否或如何能够并行地完成任务。
最明显的方法是代码以顺序方式运行。 a[i]和b[i]加载到寄存器，添加它们，将结果存储在c[i] ，并为每个i执行此操作。
然而，现代处理器具有向量或SIMD指令集，当执行相同的操作（例如，如上所示添加两个向量）时，该指令集能够在同一指令期间对数据向量进行操作。 根据处理器/体系结构的不同，也许可以在同一指令下添加a和b四个数字，而不是一次一个。
我们希望利用单指令多数据并执行数据级并行 ，例如，一次加载4个东西，每次加4个东西，一次存储4个东西。 这是代码vector化 。
请注意，这与代码并行化不同，在这种并行化中同时执行多个计算。
如果编译器识别出这样的代码块并自动引导它们，这将是一件非常困难的任务。 自动代码vector化是计算机科学中具有挑战性的研究课题。 但是随着时间的推移，编译器已经变得更好了。 你可以在这里查看 GNU-gcc的自动vector化function。 同样的， 在这里 LLVM-clang 。 你还可以在最后一个链接中find一些比较gcc和ICC （Intel C ++编译器）的基准。
gcc （我在v4.9 ）例如不能自动v4.9代码在-O2水平优化。 所以，如果我们要执行上面显示的代码，它将按顺序运行。 下面是添加两个长度为5亿的整数向量的时机。
我们需要添加flag -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
这告诉我们该函数是vector化的。 以下是在长度为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
这部分可以安全地跳过而不失连续性。
编译器并不总是有足够的信息进行vector化。 我们可以将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（都通过自制软件，MacOS 10.12.3安装的）testing的，两者都支持OpenMP 4.0。
从这个意义上说，尽pipe数组编程的维基百科页面提到，在整个数组上运行的语言通常将其称为vector化操作 ，但它实际上是循环隐藏 IMO（除非实际上是vector化的）。
在R的情况下，甚至在C中的rowSums()或colSums()代码也不会利用代码vector化 IIUC; 它只是一个循环在C.同样适用于lapply() 。 在apply()情况下，它在R.所有这些都是循环隐藏 。
简而言之，通过以下方式封装R函数：
只是写一个for循环在
C！=向量化你的代码。
只要在R写一个for循环就可以了。英特尔math核心函数库（MKL）例如实现向量化的函数forms。
HTH
参考文献：
- 英特尔公司James Reinders谈到 （这个答案大多是试图总结这个优秀的演讲）
因此，将最好的答案/评论汇总成一些一般的答案，并提供一些背景：R有4种types的循环（ 从非vector化到vector化的顺序 ）
- R
for循环，在每次迭代中重复调用R函数（ 未vector化 ） - C循环，在每次迭代中重复调用R函数（ 未vector化 ）
- 只有一次调用R函数的C循环（ 有点向量化 ）
- 一个简单的C循环，根本不调用任何 R函数，并使用它自己的编译函数（ Vectorized ）
所以*apply家庭是第二种types。 除了更多的是第一类应用
您可以从源代码中的评论中了解这一点
/ * .Internal（lapply（X，FUN））* /
/ *这是一个特殊的内部，所以有没有评价的论据。 它是
从封闭包装器中调用，所以X和FUN是承诺。 FUN必须是未评价的，例如在bquote中使用。 * /
这意味着lapply的C代码接受来自R的lapply评估的函数，并且稍后在C代码本身内进行评估。 这基本上是lapply s .Internal呼叫的区别
.Internal(lapply(X, FUN))
其中有一个有R函数的FUN参数
而colMeans .Internal调用没有 FUN参数
.Internal(colMeans(Re(x), n, prod(dn), na.rm))
colMeans不像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没有vector化 ，虽然它有两个可能的优势，比普通的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 = 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循环
虽然在对代码进行vector化时，有些事情需要考虑。
- 如果您的数据集（我们称之为
df）属于类data.frame，则某些vector化函数（如colMeans，colSums，rowSums等）必须rowSums其转换为matrix，因为这是他们的devise。 这意味着对于一个大的df这可能会造成巨大的开销。 虽然lapply不需要这样做，因为它从df提取实际的向量（因为data.frame只是一个向量列表），因此，如果没有那么多的列但行很多，lapply(df, mean)有时可能比colMeans(df)更好的select。 - 另外要记住的是R有很多不同的函数types，比如
.Primitive和generic（S3，S4） 在这里可以看到一些额外的信息。 通用函数必须做一个方法调度，有时是一个代价高昂的操作。 例如，mean是普通的S3函数，而sum是Primitive。 因此lapply(df, sum)从上面列出的原因lapply(df, sum)有些时候lapply(df, sum)可以非常有效地比较colSums
