“*apply”家族真的没有矢量化吗？
Is the "*apply" family really not vectorized?
所以我们习惯于对每个 R 新用户说“apply 没有向量化，查看 Patrick Burns R Inferno Circle 4 “上面写着（我引用）：
A common reflex is to use a function in the apply family. This is not
vectorization, it is loop-hiding. 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.
确实，快速查看 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 编写的函数。
此外，让我们以colMeans函数为例，它从未被指责没有被矢量化
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(...，我们也可以在 rabbit hole 中找到它。那么这与 .Internal(lapply(.. 有何不同？
实际上，快速基准测试表明 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 相比for 也调用 lapply) 的循环以及 Patrick Burns 真正想说的是什么？
首先，在您的示例中，您对 "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 中做所有事情的函数，所以，我猜，这应该是 "vectorised" 的意思？
我们在每个 "list"s 元素中找到平均值的示例：
(EDIT May 11 '16 ：我相信找到 "mean" 的示例对于评估R 函数迭代和编译代码，(1) 由于 R 的平均算法在 "numeric"s 上的特殊性超过简单的 sum(x) / length(x) 和 (2) 在 "list" 上测试应该更有意义s 与 length(x) >> lengths(x)。因此，"mean" 示例移至末尾并替换为另一个示例。）
作为一个简单的例子，我们可以考虑找到 "list" 的每个 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 语言编写。这具有显着的性能优势，但我认为这不是矢量化的关键 属性。矢量化从根本上说是为了拯救你的大脑，而不是拯救计算机的工作。
因此，将伟大的 answers/comments 总结为一些通用答案并提供一些背景知识：R 有 4 种类型的循环（从非矢量化到矢量化顺序 )
- R
for 在每次迭代中重复调用 R 函数的循环（未矢量化）
- 在每次迭代中重复调用 R 函数的 C 循环（未向量化）
- 仅调用 R 函数一次的 C 循环（有些向量化）
- 一个普通的 C 循环，根本不调用 任何 R 函数并使用它自己的编译函数 (Vectorized)
所以*apply一家是第二种。除了apply更多的是第一种
你可以从其source code
中的评论中了解这一点
/* .Internal(lapply(X, FUN)) */
/* This is a special .Internal, so has unevaluated arguments. It is
called from a closure wrapper, so X and FUN are promises. FUN must
be unevaluated for use in e.g. bquote .
*/
这意味着 lapplys C 代码接受来自 R 的未计算函数，然后在 C 代码本身中对其进行计算。这基本上就是lapplys .Internal call
的区别
.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 未矢量化，尽管它比普通 R for 循环
有两个可能的优势
在循环中访问和分配在 C 中似乎更快（即在 lapplying 函数中）虽然差异似乎很大，但我们仍然停留在微秒级别并且代价高昂的是每次迭代中 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 循环
尽管在向量化代码时，您需要考虑一些事项。
- 如果您的数据集（我们称之为
df）是 class data.frame，一些向量化函数（例如 colMeans、colSums、 rowSums, 等) 必须首先将其转换为矩阵，因为这是它们的设计方式。这意味着对于一个大的 df 这会产生巨大的开销。虽然 lapply 不必这样做，因为它会从 df 中提取实际向量（因为 data.frame 只是一个向量列表），因此，如果您没有那么多列但是很多行，lapply(df, mean) 有时比 colMeans(df) 更好。
- 另一件要记住的事情是 R 有多种不同的函数类型，例如
.Primitive 和泛型（S3、S4）参见 here一些额外的信息。通用函数必须进行方法分派，这有时是一项代价高昂的操作。例如，mean 是通用的 S3 函数，而 sum 是 Primitive。因此，有时 lapply(df, sum) 与 colSums 相比可能非常有效，原因如下
我同意 Patrick Burns 的观点，即 循环隐藏 而不是 代码矢量化 。原因如下：
考虑这个 C 代码片段：
for (int i=0; i<n; i++)
c[i] = a[i] + b[i]
我们想做什么很清楚。但是 如何 执行任务或如何执行任务并不是真的。默认情况下，for-loop 是一个串行结构。它不会告知事情是否可以并行完成或如何并行完成。
最明显的方式是代码 运行 以 顺序方式 。将 a[i] 和 b[i] 加载到寄存器，将它们相加，将结果存储在 c[i] 中，并对每个 i.
执行此操作
但是，现代处理器具有 vector or SIMD 指令集，能够在 同一指令 期间对 数据向量 进行操作] 在执行相同的操作时（例如，如上所示添加两个向量）。根据 processor/architecture，可能可以在同一条指令下添加来自 a 和 b 的四个数字，而不是一次一个。
We would like to exploit the Single Instruction Multiple Data and perform data level parallelism, i.e., load 4 things at a time, add 4 things at time, store 4 things at a time, for example. And this is code vectorisation.
Note that this is different from code parallelisation -- where multiple computations are performed concurrently.
如果编译器能够识别此类代码块并自动对它们进行矢量化，那就太好了，这是一项艰巨的任务。 Automatic code vectorisation is a challenging research topic in Computer Science. But over time, compilers have gotten better at it. You can check the auto vectorisation capabilities of GNU-gcc here. Similarly for LLVM-clang here。您还可以在最后 link 中找到一些与 gcc 和 ICC（Intel C++ 编译器）进行比较的基准测试。
gcc（我在 v4.9）例如不会在 -O2 级优化时自动矢量化代码。因此，如果我们要执行上面显示的代码，它将是 运行 顺序。下面是两个长度为5亿的整数向量相加的时间。
我们要么需要添加标志 -ftree-vectorize，要么将优化更改为级别 -O3。 （请注意，-O3 执行 other additional optimisations as well). The flag -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，即使是 rowSums() 或 colSums() C 中的代码也不会利用 代码矢量化 IIUC；它只是 C 中的一个循环。lapply() 也是如此。在apply()的情况下，它在R中。因此所有这些都是循环隐藏。
In short, wrapping an R function by:
just writing a for-loop in C != vectorising your code.
just writing a for-loop in R != vectorising your code.
Intel Math Kernel Library (MKL) for example implements vectorised forms of functions.
HTH
参考文献：
- Talk by James Reinders, Intel（这个回答主要是为了总结这个精彩的演讲）
所以我们习惯于对每个 R 新用户说“apply 没有向量化，查看 Patrick Burns R Inferno Circle 4 “上面写着（我引用）：
A common reflex is to use a function in the apply family. This is not vectorization, it is loop-hiding. 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.
确实，快速查看 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 编写的函数。
此外，让我们以colMeans函数为例，它从未被指责没有被矢量化
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(...，我们也可以在 rabbit hole 中找到它。那么这与 .Internal(lapply(.. 有何不同？
实际上，快速基准测试表明 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 相比for 也调用 lapply) 的循环以及 Patrick Burns 真正想说的是什么？
首先，在您的示例中，您对 "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 中做所有事情的函数，所以，我猜，这应该是 "vectorised" 的意思？
我们在每个 "list"s 元素中找到平均值的示例：
(EDIT May 11 '16 ：我相信找到 "mean" 的示例对于评估R 函数迭代和编译代码，(1) 由于 R 的平均算法在 "numeric"s 上的特殊性超过简单的 sum(x) / length(x) 和 (2) 在 "list" 上测试应该更有意义s 与 length(x) >> lengths(x)。因此，"mean" 示例移至末尾并替换为另一个示例。）
作为一个简单的例子，我们可以考虑找到 "list" 的每个 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 语言编写。这具有显着的性能优势，但我认为这不是矢量化的关键 属性。矢量化从根本上说是为了拯救你的大脑，而不是拯救计算机的工作。
因此，将伟大的 answers/comments 总结为一些通用答案并提供一些背景知识：R 有 4 种类型的循环（从非矢量化到矢量化顺序 )
- R
for在每次迭代中重复调用 R 函数的循环（未矢量化） - 在每次迭代中重复调用 R 函数的 C 循环（未向量化）
- 仅调用 R 函数一次的 C 循环（有些向量化）
- 一个普通的 C 循环，根本不调用 任何 R 函数并使用它自己的编译函数 (Vectorized)
所以*apply一家是第二种。除了apply更多的是第一种
你可以从其source code
中的评论中了解这一点/* .Internal(lapply(X, FUN)) */
/* This is a special .Internal, so has unevaluated arguments. It is
called from a closure wrapper, so X and FUN are promises. FUN must be unevaluated for use in e.g. bquote . */
这意味着 lapplys C 代码接受来自 R 的未计算函数，然后在 C 代码本身中对其进行计算。这基本上就是lapplys .Internal call
.Internal(lapply(X, FUN))
其中有一个包含 R 函数的 FUN 参数
并且 colMeans .Internal 调用 没有 有一个 FUN 参数
.Internal(colMeans(Re(x), n, prod(dn), na.rm))
colMeans，不像lapply知道确切它需要使用什么函数，因此它在C代码内部计算平均值。
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 中似乎更快（即在
lapplying 函数中）虽然差异似乎很大，但我们仍然停留在微秒级别并且代价高昂的是每次迭代中 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 循环
尽管在向量化代码时，您需要考虑一些事项。
- 如果您的数据集（我们称之为
df）是 classdata.frame，一些向量化函数（例如colMeans、colSums、rowSums, 等) 必须首先将其转换为矩阵，因为这是它们的设计方式。这意味着对于一个大的df这会产生巨大的开销。虽然lapply不必这样做，因为它会从df中提取实际向量（因为data.frame只是一个向量列表），因此，如果您没有那么多列但是很多行，lapply(df, mean)有时比colMeans(df)更好。 - 另一件要记住的事情是 R 有多种不同的函数类型，例如
.Primitive和泛型（S3、S4）参见 here一些额外的信息。通用函数必须进行方法分派，这有时是一项代价高昂的操作。例如，mean是通用的S3函数，而sum是Primitive。因此，有时lapply(df, sum)与colSums相比可能非常有效，原因如下
我同意 Patrick Burns 的观点，即 循环隐藏 而不是 代码矢量化 。原因如下：
考虑这个 C 代码片段：
for (int i=0; i<n; i++)
c[i] = a[i] + b[i]
我们想做什么很清楚。但是 如何 执行任务或如何执行任务并不是真的。默认情况下，for-loop 是一个串行结构。它不会告知事情是否可以并行完成或如何并行完成。
最明显的方式是代码 运行 以 顺序方式 。将 a[i] 和 b[i] 加载到寄存器，将它们相加，将结果存储在 c[i] 中，并对每个 i.
但是，现代处理器具有 vector or SIMD 指令集，能够在 同一指令 期间对 数据向量 进行操作] 在执行相同的操作时（例如，如上所示添加两个向量）。根据 processor/architecture，可能可以在同一条指令下添加来自 a 和 b 的四个数字，而不是一次一个。
We would like to exploit the Single Instruction Multiple Data and perform data level parallelism, i.e., load 4 things at a time, add 4 things at time, store 4 things at a time, for example. And this is code vectorisation.
Note that this is different from code parallelisation -- where multiple computations are performed concurrently.
如果编译器能够识别此类代码块并自动对它们进行矢量化，那就太好了，这是一项艰巨的任务。 Automatic code vectorisation is a challenging research topic in Computer Science. But over time, compilers have gotten better at it. You can check the auto vectorisation capabilities of GNU-gcc here. Similarly for LLVM-clang here。您还可以在最后 link 中找到一些与 gcc 和 ICC（Intel C++ 编译器）进行比较的基准测试。
gcc（我在 v4.9）例如不会在 -O2 级优化时自动矢量化代码。因此，如果我们要执行上面显示的代码，它将是 运行 顺序。下面是两个长度为5亿的整数向量相加的时间。
我们要么需要添加标志 -ftree-vectorize，要么将优化更改为级别 -O3。 （请注意，-O3 执行 other additional optimisations as well). The flag -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，即使是 rowSums() 或 colSums() C 中的代码也不会利用 代码矢量化 IIUC；它只是 C 中的一个循环。lapply() 也是如此。在apply()的情况下，它在R中。因此所有这些都是循环隐藏。
In short, wrapping an R function by:
just writing a for-loop in
C!= vectorising your code.
just writing a for-loop inR!= vectorising your code.Intel Math Kernel Library (MKL) for example implements vectorised forms of functions.
HTH
参考文献：
- Talk by James Reinders, Intel（这个回答主要是为了总结这个精彩的演讲）