A família "* apply" não é realmente vetorizada?

Então, costumamos dizer a todos os novos usuários do R que " applynão são vetorizados, confira o Círculo de Inferno Patrick Burns R 4 ", que diz (cito):

Um reflexo comum é usar uma função na família Apply. Isso não é vetorização, é uma ocultação de loop . A função de aplicação possui um loop for em sua definição. A função lapply enterra o loop, mas os tempos de execução tendem a ser aproximadamente iguais a um loop for explícito.

De fato, uma rápida olhada no applycódigo fonte revela o loop:

grep("for", capture.output(getAnywhere("apply")), value = TRUE)
## [1] "        for (i in 1L:d2) {"  "    else for (i in 1L:d2) {"

Ok, até agora, mas uma olhada lapplyou vapplyrealmente revela uma imagem completamente diferente:

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>

Então, aparentemente, não há um forloop R escondido lá, eles estão chamando a função C escrita interna.

Uma rápida olhada na toca do coelho revela praticamente a mesma imagem

Além disso, vamos pegar a colMeansfunção, por exemplo, que nunca foi acusada de não ser vetorizada

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>

Hã? Também chama apenas o .Internal(colMeans(...que também podemos encontrar na toca do coelho . Então, como isso é diferente .Internal(lapply(..?

Na verdade, um benchmark rápido revela que o sapplydesempenho não é pior colMeanse muito melhor do que um forloop para um grande conjunto de dados

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 

Em outras palavras, não é correto dizer que lapplye vapply são realmente vectorized (em comparação com applyo que é um forloop que também chama lapply) eo que Patrick Burns realmente quer dizer?

Primeiro de tudo, no seu exemplo, você faz testes em um "data.frame" que não é justo e colMeans, como eles têm uma sobrecarga: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

Em uma matriz, a imagem é um pouco diferente:

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

Analisando a parte principal da questão, a principal diferença entre lapply/ mapply/ etc e R-loops diretos é onde o loop é feito. Como observa Roland, os loops C e R precisam avaliar uma função R em cada iteração que é a mais cara. As funções C realmente rápidas são aquelas que fazem tudo em C; então, acho que deve ser isso o que "vetorizado" significa?

Um exemplo em que encontramos a média em cada um dos elementos de uma "lista":

( EDIT 11 de maio de 16 : acredito que o exemplo de encontrar a "média" não é uma boa configuração para as diferenças entre avaliar iterativamente uma função R e o código compilado, (1) devido à particularidade do algoritmo médio de R em "numérico" s sobre um simples sum(x) / length(x)e (2) deve fazer mais sentido testar em "list" s com length(x) >> lengths(x). Portanto, o exemplo "mean" é movido para o final e substituído por outro.)

Como um exemplo simples, podemos considerar a descoberta do oposto de cada length == 1elemento de uma "lista":

Em um tmp.carquivo:

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

E no lado R:

system("R CMD SHLIB /home/~/tmp.c")
dyn.load("/home/~/tmp.so")

com dados:

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

Avaliação comparativa:

#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

(Segue o exemplo original da descoberta média):

#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

Para mim, vetorização é principalmente tornar seu código mais fácil de escrever e mais fácil de entender.

O objetivo de uma função vetorizada é eliminar a contabilidade associada a um loop for. Por exemplo, em vez de:

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

Você pode escrever:

means <- vapply(mtcars, mean, numeric(1))
sds   <- vapply(mtcars, sd, numeric(1))

Isso facilita ver o que é o mesmo (os dados de entrada) e o que é diferente (a função que você está aplicando).

Uma vantagem secundária da vetorização é que o loop for geralmente é escrito em C, e não em R. Isso traz benefícios substanciais de desempenho, mas não acho que seja a propriedade chave da vetorização. A vetorização é fundamental para salvar seu cérebro, não para salvar o trabalho do computador.

Concordo com a opinião de Patrick Burns de que é um pouco oculto em loop e não uma vetorização de código . Aqui está o porquê:

Considere este Ctrecho de código:

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

O que gostaríamos de fazer é bastante claro. Mas como a tarefa é executada ou como ela pode ser executada não é realmente. Um loop for, por padrão, é uma construção serial. Não informa se ou como as coisas podem ser feitas em paralelo.

A maneira mais óbvia é que o código seja executado de maneira seqüencial . Carregue a[i]e ligue os b[i]registros, adicione-os, armazene o resultado c[i]e faça isso para cada um i.

No entanto, os processadores modernos têm um conjunto de instruções de vetor ou SIMD que é capaz de operar em um vetor de dados durante a mesma instrução ao executar a mesma operação (por exemplo, adicionando dois vetores como mostrado acima). Dependendo do processador / arquitetura, pode ser possível adicionar, digamos, quatro números de ae bsob a mesma instrução, em vez de um de cada vez.

Gostaríamos de explorar os Dados Múltiplos de Instrução Única e executar o paralelismo no nível dos dados , ou seja, carregar 4 coisas por vez, adicionar 4 coisas por vez, armazenar 4 coisas por vez, por exemplo. E isso é vetorização de código .

Observe que isso é diferente da paralelização de código - onde vários cálculos são executados simultaneamente.

Seria ótimo se o compilador identificasse esses blocos de código e os vetorizasse automaticamente , o que é uma tarefa difícil. A vetorização automática de código é um tópico de pesquisa desafiador em Ciência da Computação. Mas com o tempo, os compiladores ficaram melhores nisso. Você pode verificar os recursos de vetorização automática GNU-gcc aqui . Da mesma forma por LLVM-clang aqui . E você também pode encontrar alguns benchmarks no último link em comparação com gcce ICC(compilador Intel C ++).

gcc(Estou ativo v4.9), por exemplo, não vetoriza código automaticamente na -O2otimização de nível. Portanto, se executássemos o código mostrado acima, ele seria executado sequencialmente. Aqui está o momento para adicionar dois vetores inteiros de comprimento 500 milhões.

Precisamos adicionar a sinalização -ftree-vectorizeou alterar a otimização para o nível -O3. (Observe que também -O3realiza outras otimizações adicionais ). O sinalizador -fopt-info-vecé útil, pois informa quando um loop foi vetorizado com êxito).

# 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    

Isso nos diz que a função é vetorizada. Aqui estão os horários que comparam versões não vetorizadas e vetorizadas em vetores inteiros de comprimento 500 milhões:

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

Esta parte pode ser pulada com segurança sem perder a continuidade.

Os compiladores nem sempre têm informações suficientes para vetorizar. Poderíamos usar a especificação OpenMP para programação paralela , que também fornece uma diretiva compilador simd para instruir os compiladores a vetorizar o código. É essencial garantir que não haja sobreposição de memória, condições de corrida, etc. ao vetorizar o código manualmente, caso contrário, resultará em resultados incorretos.

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

Ao fazer isso, solicitamos especificamente ao compilador que o vetorize, não importa o quê. Precisamos ativar as extensões OpenMP usando o sinalizador de tempo de compilação -fopenmp. Ao fazer isso:

# 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

o que é ótimo! Isso foi testado com o gcc v6.2.0 e o llvm clang v3.9.0 (ambos instalados via homebrew, MacOS 10.12.3), ambos compatíveis com o OpenMP 4.0.

Nesse sentido, mesmo que a página da Wikipedia sobre Programação de matrizes mencione que os idiomas que operam em matrizes inteiras geralmente chamam isso de operações vetorizadas , na verdade, ele está ocultando o IMO (a menos que seja realmente vetorizado).

No caso de R, par rowSums()ou colSums()código em C não explora a vetorização de código IIUC; é apenas um loop em C. O mesmo vale para lapply(). No caso de apply(), é em R. Todos esses são, portanto, ocultos .

Em resumo, agrupando uma função R:

apenas escrever um loop for em C! = vectorising seu código.
apenas escrever um loop for em R! = vectorising seu código.

A Biblioteca Intel Kernel da Matemática (MKL), por exemplo, implementa formas vetorizadas de funções.

HTH

Referências:

1. Palestra de James Reinders, Intel (esta resposta é principalmente uma tentativa de resumir essa excelente palestra)

Então, para resumir as ótimas respostas / comentários em algumas respostas gerais e fornecer alguns antecedentes: R tem 4 tipos de loops ( da ordem não vetorizada à vetorizada )

1. forLoop R que chama repetidamente as funções R em cada iteração ( Não vetorizado )
2. Loop C que chama repetidamente funções R em cada iteração ( Não vetorizado )
3. Loop C que chama a função R apenas uma vez (um pouco vetorizado )
4. Um loop C simples que não chama nenhuma função R e usa suas próprias funções compiladas ( vetorizadas )

Então a *applyfamília é o segundo tipo. Exceto applyqual é mais do primeiro tipo

Você pode entender isso a partir do comentário em seu código-fonte

/ * .Internal (lapply (X, FUN)) * /

/ * Este é um .Internal especial, portanto, possui argumentos não avaliados. É
chamado de um invólucro de fechamento, portanto X e FUN são promessas. O FUN não deve ser avaliado para uso em, por exemplo, bquote. * /

Isso significa que lapplyo código s aceita uma função não avaliada de R e depois a avalia dentro do próprio código C. Esta é basicamente a diferença entre lapplya .Internalchamada s

.Internal(lapply(X, FUN))

Que possui um FUNargumento que contém uma função R

E a colMeans .Internalchamada que não tem FUNargumento

.Internal(colMeans(Re(x), n, prod(dn), na.rm))

colMeans, diferente de lapplysaber exatamente qual função ele precisa usar, calcula a média internamente no código C.

Você pode ver claramente o processo de avaliação da função R em cada iteração no lapplycódigo 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);
   }

Para resumir, lapplynão é vetorizado , embora tenha duas vantagens possíveis sobre o forloop R simples

1. Acessar e atribuir em um loop parece ser mais rápido em C (ou seja, em lapplyuma função) Embora a diferença pareça grande, ainda permanecemos no nível de microssegundos e o caro é a avaliação de uma função R em cada iteração. Um exemplo simples:

   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. Como mencionado por @Roland, ele executa um loop C compilado, em vez de um loop R interpretado

Embora ao vetorizar seu código, há algumas coisas que você precisa levar em consideração.

1. Se o seu conjunto de dados (chamado de Let It df) é de classe data.frame, algumas funções vectorized (tais como colMeans, colSums, rowSums, etc.) terá de convertê-lo em uma matriz em primeiro lugar, simplesmente porque esta é a forma como eles foram concebidos. Isso significa que, para grandes empresas, dfisso pode criar uma enorme sobrecarga. Embora lapplynão precise fazer isso, pois extrai os vetores reais de df(como data.frameé apenas uma lista de vetores) e, portanto, se você não tiver muitas colunas, mas muitas linhas, lapply(df, mean)às vezes pode ser uma opção melhor do que colMeans(df).
2. Outra coisa a lembrar é que o R possui uma grande variedade de diferentes tipos de funções, como .Primitive, e genérico ( S3, S4) , veja aqui para obter mais informações. A função genérica precisa executar um despacho de método que às vezes é uma operação cara. Por exemplo, meané genéricoS3 função while sumé Primitive. Assim, algumas vezes lapply(df, sum)pode ser muito eficiente em comparação colSumscom os motivos listados acima