因此，我们习惯于对每个 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 中完成所有操作的函数，所以，我想，这应该就是“向量化”的意义所在？
我们找到每个“列表”元素的平均值的示例:
(编辑 2016 年 5 月 11 日:我相信寻找“平均值”的示例对于迭代评估 R 函数和编译代码之间的差异来说并不是一个好的设置，(1 )由于 R 的平均算法在简单 sum(x)/length(x) 上对“数字”的特殊性，并且 (2) 对“列表”进行测试应该更有意义length(x) >> lengths(x)。因此，“mean”示例被移至末尾并替换为另一个示例。)
作为一个简单的例子，我们可以考虑查找“列表”中每个 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
关于r - "*apply"家族真的没有矢量化吗？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/28983292/