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


138

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?


8
To wszystko jest w semantyce, ale nie uważałbym ich za wektoryzowane. Rozważam podejście wektoryzowane, jeśli funkcja R jest wywoływana tylko raz i można jej przekazać wektor wartości. *applyfunkcje wielokrotnie wywołują funkcje języka R, co powoduje ich zapętlenie. Odnośnie dobrej wydajności sapply(m, mean): Prawdopodobnie kod C lapplymetody wywołuje tylko raz, a następnie wywołuje ją wielokrotnie? mean.defaultjest dość zoptymalizowany.
Roland

4
Doskonałe pytanie i dziękuję za sprawdzenie kodu źródłowego. Szukałem, czy zostało ostatnio zmienione, ale nic na ten temat nie ma w informacjach o wydaniu R od wersji 2.13.0 i nowszych.
ilir

1
W jakim stopniu wydajność zależy zarówno od platformy, jak i użytych flag kompilatora C i konsolidatora?
smci

3
@DavidArenburg Właściwie nie sądzę, że jest dobrze zdefiniowany. Przynajmniej nie znam odniesienia kanonicznego. Definicja języka wspomina o operacjach „wektoryzowanych”, ale nie definiuje wektoryzacji.
Roland

3
Bardzo powiązane: czy rodzina R's to coś więcej niż cukier syntaktyczny? (I, podobnie jak te odpowiedzi, również dobrze się czyta.)
Gregor Thomas

Odpowiedzi:


73

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

10
Świetna uwaga na temat kosztów konwersji ramki data.frame na macierz i dzięki za dostarczenie testów porównawczych.
Joshua Ulrich

To bardzo fajna odpowiedź, chociaż nie mogłem skompilować twoich funkcji all_Ci C_and_R. I również w dokumentacje compiler::cmpfunw starej wersji R lapply który zawiera rzeczywistą R forpętlę, zaczynam podejrzewać, że Burns odnosił się do tej starej wersji, która została wektoryzowane od tamtego czasu i jest to rzeczywista odpowiedź na moje pytanie .. ..
David Arenburg,

@DavidArenburg: Benchmarking la1z ?compiler::cmpfunwydaje się nadal zapewniać taką samą wydajność przy wszystkich all_Cfunkcjach. Myślę, że to - rzeczywiście - staje się kwestią definicji; czy „wektoryzowany” oznacza każdą funkcję, która akceptuje nie tylko skalary, jakąkolwiek funkcję, która ma kod w C, jakąkolwiek funkcję, która korzysta z obliczeń tylko w C?
alexis_laz

1
Wydaje mi się, że wszystkie funkcje w R mają kod C, po prostu dlatego, że wszystko w R jest funkcją (która musiała być napisana w jakimś języku). Więc w zasadzie, jeśli dobrze to rozumiem, mówisz, że lapplynie jest on wektoryzowany tylko dlatego, że ocenia funkcję R w każdej iteracji w swoim kodzie C?
David Arenburg,

5
@DavidArenburg: Gdybym musiał w jakiś sposób zdefiniować „wektoryzację”, wybrałbym podejście lingwistyczne; tj. funkcja, która akceptuje i wie, jak obsłużyć "wektor", czy to szybki, wolny, napisany w C, w R czy w czymkolwiek innym. W R, znaczenie wektoryzacji polega na tym, że wiele funkcji jest napisanych w C i obsługuje wektory, podczas gdy w innych językach użytkownicy zwykle przechodzą w pętlę nad danymi wejściowymi, aby -eg- znaleźć średnią. To sprawia, że ​​wektoryzacja wiąże się pośrednio z szybkością, wydajnością, bezpieczeństwem i solidnością.
alexis_laz

65

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.


5
Nie sądzę, aby istniała znacząca różnica w wydajności między forpętlami C i R. OK, pętla C może zostać zoptymalizowana przez kompilator, ale głównym punktem wydajności jest to, czy zawartość pętli jest wydajna. Oczywiście skompilowany kod jest zwykle szybszy niż kod interpretowany. Ale prawdopodobnie to chciałeś powiedzieć.
Roland

3
@Roland tak, to nie jest sama pętla for, to wszystko wokół niej (koszt wywołania funkcji, możliwość modyfikacji w miejscu, ...).
hadley

10
@DavidArenburg "Niepotrzebna konsekwencja jest hobgoblinem małych umysłów";)
hadley

6
Nie, nie sądzę, aby wydajność była głównym punktem wektoryzacji kodu. Przepisanie pętli na okrążenie jest korzystne, nawet jeśli nie jest szybsze. Głównym celem dplyr jest to, że ułatwia wyrażanie manipulacji danymi (i jest po prostu bardzo fajne, że jest również szybki).
hadley

12
@DavidArenburg to dlatego, że jesteś doświadczonym użytkownikiem języka R. Większość nowych użytkowników uważa, że ​​pętle są bardziej naturalne i należy ich zachęcać do wektoryzacji. Dla mnie, przy użyciu funkcji takich jak colMeans niekoniecznie jest o wektoryzacji, chodzi o ponowne szybki kod, który ktoś już napisał
Hadley

49

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)

35

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

1
Bardzo spójne podsumowanie. Kilka uwag: (1) C wie, jak obsługiwać pliki „data.frame”, ponieważ są to „listy” z atrybutami; jest to colMeansitd., które są zbudowane do obsługi tylko macierzy. (2) Jestem trochę zdezorientowany twoją trzecią kategorią; Nie mogę powiedzieć, o czym mówisz - exaclty. (3) Ponieważ odnosisz się konkretnie do tego lapply, uważam, że nie robi to różnicy między "[<-"R i C; obaj wstępnie przydzielają „listę” (SEXP) i wypełniają ją w każdej iteracji ( SET_VECTOR_ELTw C), chyba że nie rozumiem.
alexis_laz

2
Rozumiem do.call, że buduje wywołanie funkcji w środowisku C i po prostu ją ocenia; chociaż trudno mi porównać to z zapętleniem lub wektoryzacją, ponieważ robi to inaczej. Właściwie masz rację co do uzyskiwania dostępu i przypisywania różnic między C i R, chociaż oba pozostają na poziomie mikrosekund i nie wpływają znacząco na wynik wyniku, ponieważ kosztowne jest iteracyjne wywołanie funkcji R (porównaj R_loopiw R_lapplymojej odpowiedzi ). (Zmienię twój post z benchmarkiem; mam nadzieję, że nadal nie będziesz miał nic przeciwko)
alexis_laz

2
Nie próbuję się nie zgodzić - i szczerze mówiąc nie wiem, z czym się nie zgadzasz. Mój wcześniejszy komentarz mógłby być lepiej sformułowany. Próbuję doprecyzować używaną terminologię, ponieważ termin „wektoryzacja” ma dwie definicje, które są często mylone. Nie sądzę, żeby to było dyskusyjne. Burns i wydaje się, że chcesz go używać tylko w sensie implementacji, ale Hadley i wielu członków R-Core ( Vectorize()na przykładzie) używają go również w sensie interfejsu użytkownika. Myślę, że duża część nieporozumień w tym wątku jest spowodowana użyciem jednego terminu dla dwóch oddzielnych, ale powiązanych ze sobą koncepcji.
Gregor Thomas

3
@DavidArenburg i czy to nie jest wektoryzacja w sensie interfejsu użytkownika, niezależnie od tego, czy pod spodem znajduje się pętla for w R czy C?
Gregor Thomas

2
@DavidArenburg, Gregor, myślę, że istnieje pomyłka między „wektoryzacją kodu” a „wektoryzacją funkcji”. W R użycie wydaje się skłonne do tego drugiego. „Wektoryzacja kodu” opisuje działanie na wektorze o długości „k” w tej samej instrukcji. Zawijanie fn. wokół zapętlonego kodu skutkuje „wektoryzowanymi funkcjami” (tak, to nie ma sensu i jest mylące, zgadzam się, lepiej byłoby ukrywać pętlę lub wektorowe funkcje i / p ) i nie musi mieć nic wspólnego z wektoryzacją kodu . W R zastosuj byłby funkcją wektoryzowaną , ale nie wektoryzuje kodu, a raczej operuje na wektorach.
Arun
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.