“*apply”家族真的没有向量化吗？
Posted
技术标签:
【中文标题】“*apply”家族真的没有向量化吗？【英文标题】：Is the "*apply" family really not vectorized? 【发布时间】：2015-05-13 01:06:37 【问题描述】：所以我们习惯于对每个 R 新用户说“apply 没有矢量化，请查看 Patrick Burns R Inferno Circle 4”，它说（我引用）：
一个常见的反应是使用 apply 系列中的一个函数。 这不是 向量化，它是循环隐藏的。 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 编写的函数。
快速查看rabbit hole 发现几乎相同的图片
此外，我们以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）帕特里克·伯恩斯到底想说什么？
【问题讨论】：
这都是语义上的，但我不会认为它们是矢量化的。如果一个 R 函数只被调用一次并且可以传递一个值向量，我会考虑一种向量化的方法。*apply 函数重复调用 R 函数，这使得它们循环。关于sapply(m, mean) 的良好性能：可能lapply 的C 代码只分派一次方法，然后重复调用该方法？ mean.default 非常优化。
很好的问题，感谢您检查底层代码。我正在查看它最近是否已更改，但从 2.13.0 版开始的 R 发行说明中没有关于此的内容。
性能在多大程度上取决于平台以及使用的 C 编译器和链接器标志？
@DavidArenburg 实际上，我认为它的定义并不明确。至少我不知道规范的参考。语言定义提到了“向量化”操作，但没有定义向量化。
非常相关：Is R's apply family more than syntactic sugar?（而且，就像这些答案一样，也是一本好书。）
【参考方案1】：
首先，在您的示例中，您在“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-loops 之间的主要区别在于循环的完成位置。正如 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
【讨论】：
关于将 data.frame 转换为矩阵的成本非常重要，感谢您提供基准。 这是一个非常好的答案，尽管我无法编译您的all_C 和 C_and_R 函数。我还在 compiler::cmpfun 的文档中发现了一个 old R version of lapply，其中包含一个实际的 R for 循环，我开始怀疑 Burns 指的是 那个 i> 从那时起被矢量化的旧版本，这是我问题的实际答案....
@DavidArenburg ：从?compiler::cmpfun 对la1 进行基准测试似乎仍然可以在除all_C 之外的所有函数中产生相同的效率。我想，这确实是一个定义问题。 “向量化”是指任何不仅接受标量的函数、任何具有 C 代码的函数、任何仅使用 C 语言计算的函数吗？
我猜 R 中的 all 函数中都有 C 代码，仅仅是因为 R 中的 everything 是一个函数（必须用一些语言）。所以基本上，如果我理解正确，你是说lapply 不是矢量化的，因为它在每个迭代中用它的C代码评估一个R函数？
@DavidArenburg ：如果我必须以某种方式定义“矢量化”，我想，我会选择一种语言方法；即一个接受并知道如何处理“向量”的函数，无论它是快的、慢的、用 C 语言、R 语言或其他任何东西编写的。在 R 中，向量化的重要性在于许多函数是用 C 编写的并处理向量，而在其他语言中，用户通常会循环输入以 - 例如 - 找到平均值。这使得矢量化与速度、效率、安全性和稳健性间接相关。【参考方案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 语言编写的。这具有显着的性能优势，但我不认为它是矢量化的关键属性。矢量化从根本上说是为了节省您的大脑，而不是节省计算机工作。
【讨论】：
我认为 C 和 Rfor 循环之间没有有意义的性能差异。好的，编译器可能会优化 C 循环，但性能的重点是循环的内容是否有效。显然编译的代码通常比解释的代码更快。但这可能就是你想说的。
@Roland 是的，它本身不是 for 循环本身，而是它周围的所有东西（函数调用的成本、就地修改的能力……）。跨度>
@DavidArenburg “不必要的一致性是小脑袋的妖精”;)
不，我不认为性能是矢量化代码的重点。将循环重写为 lapply 是有益的，即使它没有更快。 dplyr 的主要特点是它可以更轻松地表达数据操作（而且速度也很快）。
@DavidArenburg 那是因为您是一位经验丰富的 R 用户。大多数新用户发现循环更加自然，需要鼓励向量化。对我来说，使用 colMeans 之类的函数不一定是关于矢量化，而是关于重用某人已经编写的快速代码【参考方案3】：
因此，将出色的答案/cmets 总结为一些通用答案并提供一些背景知识：R 有 4 种类型的循环（从非矢量化到矢量化顺序）
-
R
for 循环在每次迭代中重复调用 R 函数（未矢量化）
在每次迭代中重复调用 R 函数的 C 循环（未矢量化）
只调用一次 R 函数的 C 循环（有些矢量化）
一个纯 C 循环，根本不调用 任何 R 函数并使用它自己的编译函数（矢量化）
所以*apply 系列是第二种类型。除了apply 更多的是第一类
您可以从其source code中的评论中理解这一点
/* .Internal(lapply(X, FUN)) */
/* 这是一个特殊的 .Internal，因此有未评估的参数。它是 从闭包包装器调用，所以 X 和 FUN 是承诺。有趣的必须 未经评估以用于例如报价。 */
这意味着lapplys C 代码接受来自 R 的未评估函数，然后在 C 代码本身中对其进行评估。这基本上就是lapplys .Internal调用的区别
.Internal(lapply(X, FUN))
其中有一个 FUN 参数，其中包含一个 R 函数
而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 中更快（即在 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）属于data.frame 类，则某些矢量化函数（例如colMeans、colSums、rowSums 等）必须将其转换为首先是矩阵，仅仅是因为它们是这样设计的。这意味着对于一个大的df，这会产生巨大的开销。虽然lapply 不必这样做，因为它会从df 中提取实际向量（因为data.frame 只是一个向量列表），因此，如果您没有那么多列但有很多行，@987654353 @ 有时比 colMeans(df) 更好。
要记住的另一件事是 R 具有多种不同的函数类型，例如 .Primitive 和泛型（S3、S4）请参阅 here 了解更多信息。通用函数必须执行方法分派，这有时是一项代价高昂的操作。例如，mean 是通用的 S3 函数，而 sum 是 Primitive。因此，由于上面列出的原因，有时lapply(df, sum) 与colSums 相比可能非常有效
【讨论】：
非常有凝聚力的总结。只需注意几点： (1) C 知道如何处理“data.frame”，因为它们是带有属性的“列表”；就是colMeans 等是为仅处理矩阵而构建的。 （2）我对你的第三类有点困惑；我不知道你指的是什么-exaclty-。 （3）由于您专门指的是lapply，我相信R和C中的"[<-"没有区别；他们都预先分配了一个“列表”（一个 SEXP）并在每次迭代中填写它（C 中的SET_VECTOR_ELT），除非我错过了你的观点。
我明白你对do.call 的看法，因为它在C 环境中构建了一个函数调用并对其进行评估；虽然我很难将它与循环或矢量化进行比较，因为它做了不同的事情。实际上，您在访问和分配 C 和 R 之间的差异方面是正确的，尽管两者都停留在微秒级别并且不会极大地影响结果，因为代价高昂的是迭代 R 函数调用（比较 R_loop 和R_lapply 在我的回答中）。 （我会用基准编辑你的帖子；我希望你仍然不会介意）
我并不是要不同意——老实说，我对您不同意的内容感到困惑。我之前的评论本来可以更好的措辞。我正在尝试改进正在使用的术语，因为术语“矢量化”有两个经常混为一谈的定义。我不认为这是有争议的。 Burns 和您似乎只想在实现意义上使用它，但 Hadley 和许多 R-Core 成员（以 Vectorize() 为例）也在 UI 意义上使用它。我认为这个帖子中的大部分分歧是由于使用一个术语来表示两个独立但相关的概念。
@DavidArenburg，这不是 UI 意义上的矢量化，不管下面的 R 或 C 中是否有 for 循环？
@DavidArenburg，Gregor，我认为混淆是在“代码矢量化”和“矢量化函数”之间。在 R 中，用法似乎倾向于后者。 “代码向量化”描述了在同一指令中对长度为“k”的向量进行操作。包装一个 fn。围绕循环代码导致“矢量化函数”（是的，它没有意义并且令人困惑，我同意，更好的是 loop hidden 或 vector i/p functions ) 并且不需要与 代码矢量化 有任何关系。在 R 中，apply 将是一个矢量化函数，但它不会对您的代码进行矢量化，而是对矢量进行操作。【参考方案4】：
我同意 Patrick Burns 的观点，即它是 循环隐藏 而不是 代码矢量化。原因如下：
考虑这个C代码sn-p：
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 指令集，当执行相同的操作（例如，添加如上图所示的两个向量）。根据处理器/架构，可能会在同一指令下添加来自a 和b 的四个数字，而不是一次添加一个。
我们想利用Single Instruction Multiple Data 并执行数据级并行性，例如，一次加载 4 个内容，一次添加 4 个内容，一次存储 4 个内容。这就是代码矢量化。
请注意，这与代码并行化不同——并行执行多个计算。
如果编译器能识别出这样的代码块并自动将它们向量化，那就太好了，这是一项艰巨的任务。 Automatic code vectorisation 是计算机科学中一个具有挑战性的研究课题。但是随着时间的推移，编译器已经变得更好了。您可以检查GNU-gcchere 的自动矢量化 功能。 LLVM-clanghere 也是如此。您还可以在最后一个链接中找到一些与 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-loop ！= 向量化您的代码。 只需在R中编写一个 for-loop ！= 向量化您的代码。Intel Math Kernel Library (MKL) 例如实现函数的向量化形式。
HTH
参考资料：
-
Talk by James Reinders, Intel（本回答主要是对这次精彩演讲的总结）
【讨论】：
以上是关于“*apply”家族真的没有向量化吗？的主要内容，如果未能解决你的问题，请参考以下文章