“申请”家庭是否真的没有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) {" 

那么到目前为止,但看一下lapplyvapply实际上揭示了一个完全不同的画面:

 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 

换句话说,是否说实际上是向lapplyvapply 实际上是向量化的 (相对于apply这个也称为lapply循环),帕特里克·伯恩斯究竟意味着什么呢?

4 Solutions collect form web for ““申请”家庭是否真的没有vector化?”

首先,在你的例子中,你对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指令集,当执行相同的操作(例如,如上所示添加两个向量)时,该指令集能够在同一指令期间对数据向量进行操作。 根据处理器/体系结构的不同,也许可以在同一指令下添加ab四个数字,而不是一次一个。

我们希望利用单指令多数据并执行数据级并行 ,例如,一次加载4个东西,每次加4个东西,一次存储4个东西。 这是代码vector化

请注意,这与代码并行化不同,在这种并行化中同时执行多个计算。

如果编译器识别出这样的代码块并自动引导它们,这将是一件非常困难的任务。 自动代码vector化是计算机科学中具有挑战性的研究课题。 但是随着时间的推移,编译器已经变得更好了。 你可以在这里查看 GNU-gcc自动vector化function。 同样的, 在这里 LLVM-clang 。 你还可以在最后一个链接中find一些比较gccICC (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


参考文献:

  1. 英特尔公司James Reinders谈到 (这个答案大多是试图总结这个优秀的演讲)

因此,将最好的答案/评论汇总成一些一般的答案,并提供一些背景:R有4种types的循环( 从非vector化到vector化的顺序

  1. R for循环,在每次迭代中重复调用R函数( 未vector化
  2. C循环,在每次迭代中重复调用R函数( 未vector化
  3. 只有一次调用R函数的C循环( 有点向量化
  4. 一个简单的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循环

  1. 在循环中访问和分配似乎在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 
  2. 正如@Roland所提到的,它运行一个编译的C循环而不是一个解释的R循环


虽然在对代码进行vector化时,有些事情需要考虑。

  1. 如果您的数据集(我们称之为df )属于类data.frame ,则某些vector化函数(如colMeanscolSumsrowSums等)必须rowSums其转换为matrix,因为这是他们的devise。 这意味着对于一个大的df这可能会造成巨大的开销。 虽然lapply不需要这样做,因为它从df提取实际的向量(因为data.frame只是一个向量列表),因此,如果没有那么多的列但行很多, lapply(df, mean)有时可能比colMeans(df)更好的select。
  2. 另外要记住的是R有很多不同的函数types,比如.Primitive和generic( S3S4 ) 在这里可以看到一些额外的信息。 通用函数必须做一个方法调度,有时是一个代价高昂的操作。 例如, mean是普通的S3函数,而sumPrimitive 。 因此lapply(df, sum)从上面列出的原因lapply(df, sum)有些时候lapply(df, sum)可以非常有效地比较colSums
  • 如果线程启动Executor,则无法从Future <?>和SwingWorker获取ArrayIndexOutOfBoundsException
  • 用硬编码元素初始化std :: vector最简单的方法是什么?
  • vector在C ++中的<int> :: size_type
  • 从相同的向量中push_back元素是否安全?
  • 启用C ++ 11时,std :: vector的性能回归
  • VBO与std :: vector
  • 载体在Arduino
  • C ++vectorstring
  • Androidselect器可用VectorDrawables srcCompat绘制
  • 计算2Dvector的交叉产品
  • 擦除从一个std ::向量,而为每个?