Czy rodzina „* zastosuj” naprawdę nie jest wektoryzowana?

Więc przyzwyczailiśmy się do każdego nowego użytkownika R, że „ applynie jest wektoryzowany, sprawdź Patrick Burns R Inferno Circle 4 ”, który mówi (cytuję):

Powszechnym odruchem jest użycie funkcji w rodzinie stosowanej. To nie jest wektoryzacja, to ukrywanie pętli . Funkcja zastosuj w swojej definicji pętlę for. Funkcja lapply zakopuje pętlę, ale czasy wykonywania są z grubsza równe jawnej pętli for.

Rzeczywiście, szybkie spojrzenie na applykod źródłowy ujawnia pętlę:

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

Na razie ok, ale spojrzenie na lapplylub vapplyfaktycznie ujawnia zupełnie inny obraz:

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>

Więc najwyraźniej nie ma tam ukrytej forpętli R , raczej wywołują wewnętrzną funkcję napisaną w C.

Szybkie spojrzenie w króliczą nory ujawnia prawie ten sam obraz

Ponadto weźmy colMeansna przykład funkcję, której nigdy nie zarzucano, że nie jest wektoryzowana

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>

Co? To także tylko połączenia, .Internal(colMeans(...które możemy również znaleźć w króliczej norze . Więc czym się to różni od .Internal(lapply(..?

Właściwie szybki test porównawczy ujawnia, że sapplydziała nie gorzej niż colMeansi znacznie lepiej niż forpętla dla dużego zbioru danych

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 

Innymi słowy, czy słuszne jest powiedzenie tego lapplyi vapply faktycznie są one wektoryzowane (w porównaniu z applyktórym jest forpętla, która również wywołuje lapply) i co naprawdę miał na myśli Patrick Burns?

Po pierwsze, w twoim przykładzie wykonujesz testy na „data.frame”, co nie jest sprawiedliwe colMeans, applya "[.data.frame"ponieważ mają one narzut:

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

Na matrycy obraz wygląda trochę inaczej:

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

Biorąc pod uwagę główną część pytania, główną różnicą między lapply/ mapply/ etc a prostymi pętlami R jest miejsce, w którym pętle są wykonywane. Jak zauważa Roland, obie pętle C i R muszą oceniać funkcję R w każdej iteracji, która jest najbardziej kosztowna. Naprawdę szybkie funkcje C to te, które robią wszystko w C, więc myślę, że o to właśnie chodzi w „wektoryzacji”?

Przykład, w którym znajdujemy średnią w każdym z elementów „listy”:

( EDYTUJ 11 maja 2016 : Uważam, że przykład ze znalezieniem "średniej" nie jest dobrą konfiguracją dla różnic między iteracyjną oceną funkcji R a skompilowanym kodem, (1) ze względu na specyfikę algorytmu średniej R na "numerycznej" s ponad prostym sum(x) / length(x)i (2) bardziej sensowne powinno być testowanie na „liście” z length(x) >> lengths(x). Tak więc przykład „średni” jest przenoszony na koniec i zastępowany innym).

Jako prosty przykład możemy rozważyć znalezienie przeciwieństwa każdego length == 1elementu „listy”:

W tmp.cpliku:

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

A po stronie R:

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

z danymi:

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

Benchmarking:

#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

(Zgodnie z oryginalnym przykładem ustalenia średniej):

#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

Dla mnie wektoryzacja to przede wszystkim ułatwienie pisania i zrozumienia kodu.

Celem funkcji wektoryzowanej jest wyeliminowanie księgowości związanej z pętlą for. Na przykład zamiast:

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

Możesz pisać:

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

Dzięki temu łatwiej jest zobaczyć, co jest takie samo (dane wejściowe), a co inne (funkcja, którą stosujesz).

Drugorzędną zaletą wektoryzacji jest to, że pętla for jest często zapisywana w języku C, a nie R. Ma to znaczące korzyści w zakresie wydajności, ale nie sądzę, że jest to kluczowa właściwość wektoryzacji. Wektoryzacja polega zasadniczo na ratowaniu mózgu, a nie pracy komputera.

Zgadzam się z poglądem Patricka Burnsa, że ​​jest to raczej ukrywanie w pętli, a nie wektoryzacja kodu . Dlatego:

Rozważ ten Cfragment kodu:

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

Co co chcielibyśmy zrobić, jest całkiem jasne. Ale to, jak zadanie jest wykonywane lub jak można je wykonać, nie jest tak naprawdę. Dla pętli domyślnie jest szeregowym konstrukt. Nie informuje, czy i jak można coś zrobić równolegle.

Najbardziej oczywistym sposobem jest to, że kod jest uruchamiany sekwencyjnie . Załaduja[i] i b[i]przejdź do rejestrów, dodaj je, zapisz wynik c[i]i zrób to dla każdego i.

Jednak nowoczesne procesory mają zestaw instrukcji wektorowych lub SIMD, które mogą działać na wektorze danych podczas tej samej instrukcji podczas wykonywania tej samej operacji (np. Dodawanie dwóch wektorów, jak pokazano powyżej). W zależności od procesora / architektury może być możliwe dodanie, powiedzmy, czterech liczb z ai bpod tą samą instrukcją zamiast jednej naraz.

Chcielibyśmy wykorzystać jedną instrukcję i wiele danych i wykonać równoległość na poziomie danych , czyli obciążenie 4 rzeczy na raz, dodać 4 rzeczy na raz, sklep 4 rzeczy naraz, na przykład. A to jest wektoryzacja kodu .

Zwróć uwagę, że różni się to od równoległości kodu - gdzie wiele obliczeń jest wykonywanych jednocześnie.

Byłoby wspaniale, gdyby kompilator zidentyfikował takie bloki kodu i automatycznie je wektoryzował, co jest trudnym zadaniem. Automatyczna wektoryzacja kodu to trudny temat badawczy w informatyce. Ale z biegiem czasu kompilatory stały się w tym coraz lepsze. Możesz sprawdzić możliwości autowektoryzacji GNU-gcc tutaj . Podobnie LLVM-clang tutaj . W ostatnim linku można również znaleźć testy porównawczegcc i ICC(kompilator Intel C ++).

gcc (Jestem na v4.9 ) nie wektoryzuje kodu automatycznie przy -O2optymalizacji poziomu. Więc gdybyśmy mieli wykonać kod pokazany powyżej, byłby uruchamiany sekwencyjnie. Oto czas dodawania dwóch wektorów całkowitych o długości 500 milionów.

Musimy albo dodać flagę, -ftree-vectorizealbo zmienić optymalizację na poziom -O3. (Zauważ, że -O3wykonuje również inne dodatkowe optymalizacje ). Flaga -fopt-info-vecjest przydatna, ponieważ informuje, kiedy pętla została pomyślnie wektoryzowana).

# 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    

To mówi nam, że funkcja jest wektoryzowana. Oto czasy porównania wersji niewektoryzowanych i wektoryzowanych na wektorach całkowitych o długości 500 milionów:

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

Tę część można bezpiecznie pominąć bez utraty ciągłości.

Kompilatory nie zawsze będą miały wystarczającą ilość informacji do wektoryzacji. Moglibyśmy użyć specyfikacji OpenMP do programowania równoległego , która również zapewnia rozszerzenie dyrektywę kompilatora simd, aby poinstruować kompilatory, aby wektoryzowały kod. Podczas ręcznej wektoryzacji kodu należy upewnić się, że nie ma nakładania się pamięci, warunków wyścigu itp. W przeciwnym razie spowoduje to nieprawidłowe wyniki.

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

Robiąc to, prosimy kompilator o wektoryzację bez względu na wszystko. Będziemy musieli aktywować rozszerzenia OpenMP za pomocą flagi czasu kompilacji-fopenmp . Robiąc to:

# 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

który jest świetny! Zostało to przetestowane z gcc v6.2.0 i llvm clang v3.9.0 (oba zainstalowane przez homebrew, MacOS 10.12.3), z których oba obsługują OpenMP 4.0.

W tym sensie, mimo że strona Wikipedii na temat programowania tablic wspomina, że ​​języki, które operują na całych tablicach, zwykle nazywają to operacjami wektoryzowanymi , tak naprawdę jest ukrywa pętlę IMO (chyba że jest faktycznie wektoryzowana).

W przypadku R, nawet rowSums()lub colSums()kod w C nie wykorzystuj wektoryzacji kodu IIUC; to tylko pętla w C. To samo dotyczy lapply(). W przypadku apply(), jest w R. Wszystkie te elementy są więc ukrywane w pętli .

Krótko mówiąc, zawijanie funkcji R przez:

właśnie pisanie dla pętli w C! = vectorising kodu.
właśnie pisanie dla pętli wR ! = vectorising kodu.

Na przykład Intel Math Kernel Library (MKL) implementuje wektoryzowane formy funkcji.

HTH

Bibliografia:

1. Przemówienie Jamesa Reindersa, Intel (ta odpowiedź jest głównie próbą podsumowania tego znakomitego przemówienia)

Podsumowując świetne odpowiedzi / komentarze w ogólną odpowiedź i dostarczając trochę informacji: R ma 4 typy pętli ( w kolejności od nie wektoryzowanej do zwektoryzowanej )

1. forPętla R, która wielokrotnie wywołuje funkcje języka R w każdej iteracji ( nie zwektoryzowany )
2. Pętla C, która wielokrotnie wywołuje funkcje języka R w każdej iteracji ( nie wektoryzowane )
3. Pętla C, która wywołuje funkcję języka R tylko raz ( nieco wektoryzowana )
4. Zwykła pętla C, która w ogóle nie wywołuje żadnej funkcji R i używa jej własnych skompilowanych funkcji ( wektoryzacja )

Więc *applyrodzina jest drugim typem. Z wyjątkiem tego, applyktóry jest bardziej pierwszego typu

Możesz to zrozumieć z komentarza w jego kodzie źródłowym

/ *. Wewnętrzna (lapply (X, FUN)) * /

/ * To jest specjalny .Wewnętrzny, więc ma nieocenione argumenty. Nazywa się to
z opakowania zamknięcia, więc X i FUN to obietnice. ZABAWA musi być nieoceniona do wykorzystania np. W bquote. * /

Oznacza to, że lapplykod w języku C przyjmuje nieocenioną funkcję z języka R, a następnie ocenia ją w samym kodzie C. Jest to po prostu różnica między lapplyów .Internalrozmowy

.Internal(lapply(X, FUN))

Który ma FUNargument, który przechowuje funkcję R.

I colMeans .Internalwezwanie, które nie ma FUNargumentu

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

colMeans, w przeciwieństwie do tego, lapplywie dokładnie, jakiej funkcji potrzebuje, więc oblicza średnią wewnętrznie w kodzie C.

Możesz wyraźnie zobaczyć proces oceny funkcji R w każdej iteracji w lapplykodzie 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);
   }

Podsumowując, lapplynie jest wektoryzowany , chociaż ma dwie możliwe zalety w stosunku do zwykłej forpętli R.

1. Dostęp i przypisywanie w pętli wydaje się być szybsze w języku C (tj. W lapplyfunkcji). Chociaż różnica wydaje się duża, nadal pozostajemy na poziomie mikrosekund, a kosztowna jest wycena funkcji R w każdej iteracji. Prosty przykład:

   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. Jak wspomniał @Roland, uruchamia skompilowaną pętlę C, a raczej zinterpretowaną pętlę R.

Chociaż podczas wektoryzacji kodu należy wziąć pod uwagę kilka rzeczy.

1. Jeśli zbiór danych (call Chodźmy go df) jest klasy data.frame, niektóre vectorized funkcje (takie jak colMeans, colSums, rowSums, itd.) Będą musiały przekształcić go do macierzy po pierwsze, po prostu dlatego, że jest to w jaki sposób zostały one zaprojektowane. Oznacza to, że w przypadku dużego dfmoże to spowodować ogromne obciążenie. Chociaż lapplynie będzie musiał tego robić, ponieważ wyodrębnia rzeczywiste wektory z df(jak data.framejest to tylko lista wektorów), a zatem, jeśli masz nie tyle kolumn, ale wiele wierszy, lapply(df, mean)czasami może być lepszą opcją niż colMeans(df).
2. Inną rzeczą do zapamiętania jest to, że R ma wiele różnych typów funkcji, takich jak .Primitive, i generic ( S3, S4), patrz tutaj, aby uzyskać dodatkowe informacje. Funkcja ogólna musi wykonać wysyłkę metody, która czasami jest kosztowną operacją. Na przykład meanjest S3funkcją ogólną, podczas gdy sumjest Primitive. Dlatego czasami lapply(df, sum)może być bardzo wydajne w porównaniu colSumsz przyczynami wymienionymi powyżej