"*apply"系列真的没有矢量化吗？

Question

因此,我们习惯于对每个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编写的函数.

兔子 洞的快速浏览显示了几乎相同的图片

此外,让我们以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(...我们也可以在兔子洞里找到的.那有什么不同.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)以及帕特里克·伯恩斯的真正含义是什么？

Answer 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-loop 之间的主要区别在于循环完成的地方.正如Roland所说,C和R循环都需要在每次迭代中评估R函数,这是最昂贵的.真正快速的C函数是用C做的所有事情,所以,我想,这应该是"矢量化"的意思吗？

我们在每个"list"元素中找到平均值的示例:

(编辑2016年5月11日:我认为找到"均值"的例子不是迭代评估R函数和编译代码之间差异的良好设置,(1)因为R的均值算法对"数字"的特殊性简单地说sum(x) / length(x),(2)用"list"来测试它应该更有意义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

Answer 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语言编写.这有很大的性能优势,但我认为它不是向量化的关键属性.矢量化基本上是关于拯救大脑,而不是保存计算机工作.

Answer 3

我同意Patrick Burns的观点,即它隐藏了循环,而不是代码矢量化.原因如下:

请考虑以下C代码段:

for (int i=0; i<n; i++)
  c[i] = a[i] + b[i]

什么我们希望做的是相当清楚的.但是,如何执行任务或如何执行任务并非如此.一个for循环在默认情况下是串行结构.它不会告知是否可以或如何并行完成任务.

最明显的方式是代码以顺序方式运行.加载a[i]并b[i]打开寄存器,添加它们,存储结果c[i],并为每个寄存器执行此操作i.

然而,现代处理器具有矢量或SIMD指令集,其能够在执行相同操作时在同一指令期间对数据矢量进行操作(例如,如上所示添加两个矢量).根据处理器/架构,有可能从增加,也就是说,四个数字和相同的指令下,而不是一次一个.ab

我们想利用单指令多数据并执行数据级并行,即一次加载4个东西,一次添加4个东西,例如一次存储4个东西.这是代码矢量化.

请注意,这与代码并行化不同 - 其中多个计算同时执行.

如果编译器识别出这样的代码块并自动对它们进行矢量化,那就太好了,这是一项艰巨的任务.自动代码矢量化是计算机科学中具有挑战性的研究课题.但随着时间的推移,编译器已经变得更好了.您可以在此处查看自动矢量化功能.同样在这里.您还可以在最后一个链接中找到一些基准测试与(英特尔C++编译器)进行比较.GNU-gcc LLVM-clang gccICC

gcc(我正在v4.9)例如,在-O2级别优化时不会自动向量化代码.因此,如果我们要执行上面显示的代码,它将按顺序运行.这是添加两个长度为5亿的整数向量的时间.

我们要么需要添加标志,-ftree-vectorize要么将优化更改为level -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    

这告诉我们该函数是矢量化的.以下是在长度为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规范进行并行编程,它还提供了一个simd编译器指令来指示编译器对代码进行矢量化.必须确保没有内存重叠,竞争条件等.手动向量化代码时,否则会导致错误的结果.

#pragma omp simd
for (i=0; i<n; i++) 
  c[i] = a[i] + b[i]

通过这样做,我们特别要求编译器将其矢量化,无论如何.我们需要使用编译时标志来激活OpenMP扩展-fopenmp.通过这样做:

# 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.

从这个意义上讲,尽管关于数组编程的维基百科页面提到在整个数组上运行的语言通常将其称为向量化操作,但实际上它是隐藏 IMO的循环(除非它实际上是向量化的).

在R的情况下,偶数rowSums()或colSums()C中的代码不利用代码矢量化 IIUC; 它只是C中的一个循环.同样适用lapply().如果是apply(),它在R中.因此所有这些都是循环隐藏.

简而言之,通过以下方式包装R函数:

只需在!=中编写一个for循环C代码.
只需在!=中编写一个for循环R代码.

例如,英特尔数学核心库(MKL)实现了矢量化形式的函数.

HTH

参考文献:

1. 英特尔James Reinders的演讲(这个答案主要是试图总结这个优秀的演讲)

Answer 4

所以将好的答案/评论总结成一些一般的答案并提供一些背景:R有4种类型的循环(从非矢量化到矢量化顺序)

1. for在每次迭代中重复调用R函数的R 循环(未向量化)
2. 在每次迭代中重复调用R函数的C循环(未向量化)
3. 只调用一次R函数的C循环(Somewhat vectorized)
4. 一个普通的C循环,根本不调用任何 R函数并使用它自己编译的函数(Vectorized)

所以这个*apply家庭是第二种.除了apply哪个更像是第一种类型

您可以从其源代码中的注释中理解这一点

/*.内部(lapply(X,FUN))*/

/*这是一个特殊的.Internal,所以有未评估的参数.它是
从一个闭包装中调用的,因此X和FUN是承诺.必须对FUN进行评估,以便在例如bquote中使用.*/

这意味着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代码内的平均值.

您可以在C代码中的每次迭代中清楚地看到R函数的评估过程lapply

 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);
   }

总而言之,虽然它比普通的R 循环有两个可能的优点,但它lapply没有矢量化for

1. 在循环中访问和分配似乎在C中更快(即在lapply函数中)虽然差异看起来很大,但我们仍然保持在微秒级别并且代价高昂的是每次迭代中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循环

虽然在对代码进行矢量化时,您需要考虑一些事项.

1. 如果设置你的数据(我们称之为df)是一流的data.frame,一些量化的功能(如colMeans,colSums,rowSums等),必须把它转换到一个矩阵,只是因为这是他们是如何设计的.这意味着对于一个大的df这可能会产生巨大的开销.虽然lapply它不必这样做,因为它提取实际的向量df(因为data.frame只是一个向量列表),因此,如果你没有那么多的列但很多行,lapply(df, mean)有时可能是更好的选择colMeans(df).
2. 要记住的另一件事是R有各种不同的函数类型,例如.Primitive,和泛型(S3,S4)在这里可以看到一些其他信息.通用函数必须执行方法调度,有时这是一种代价高昂的操作.例如,mean是通用的S3功能,而sum为Primitive.因此,lapply(df, sum)与colSums上面列出的原因相比,有时可能非常有效