Czy rodzina "* apply " naprawdę nie jest wektoryzowana?

Dla mnie wektoryzacja polega przede wszystkim na tym, aby Twój kod był łatwiejszy do napisania i łatwiejszy do zrozumienia.

Celem wektoryzowanej funkcji 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 napisać:

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

To ułatwia sprawdzenie, co jest takie samo (dane wejściowe), a co inne (funkcja, którą stosujesz).

Dodatkową zaletą wektoryzacji jest to, że pętla for jest często zapisywana w C, raczej niż w R. ma to istotne zalety, ale nie sądzę, że jest to kluczowa właściwość wektoryzacji. Wektoryzacja polega zasadniczo na ratowaniu mózgu, a nie na ratowaniu pracy komputera.

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

Rozważ ten C fragment kodu:

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

to, co chcielibyśmy zrobić, jest dość jasne. Ale to, jak zadanie jest wykonywane lub jak może być wykonane, nie jest tak naprawdę. A for-loop domyślnie jest konstrukcją szeregową. Nie informuje, czy i w jaki sposób można robić rzeczy równolegle.

Najbardziej oczywistym sposobem jest to, że kod jest uruchamiany w sposób sekwencyjny. Wczytaj a[i] i b[i] do rejestrów, dodaj je, zapisz wynik w c[i] i zrób to dla każdego i.

[[30]} chcielibyśmy wykorzystać pojedyncza Instrukcja wiele danych i wykonać równoległość poziomu danych, np. załadować 4 rzeczy na raz, dodać 4 rzeczy na raz, zapisać 4 rzeczy na raz, na przykład. A to jest wektoryzacja kodu .

Zauważ, że różni się to od paralelizacji kodu - gdzie wiele obliczenia są wykonywane jednocześnie.

Byłoby wspaniale, gdyby kompilator zidentyfikował takie bloki kodu i automatycznie wektoryzował je, co jest trudnym zadaniem. Automatyczna wektoryzacja kodu to trudny temat badawczy w informatyce. Ale z czasem Kompilatory stały się w tym lepsze. Możesz sprawdzić możliwości auto wektoryzacji GNU-gcc tutaj . Podobnie dla LLVM-clang tutaj . I można również znajdź kilka benchmarków w ostatnim linku w porównaniu z gcc i ICC (kompilator Intel C++).

gcc (jestem na v4.9) na przykład nie wektoryzuje kodu automatycznie przy optymalizacji poziomu -O2. Więc jeśli wykonamy kod pokazany powyżej, będzie on uruchamiany sekwencyjnie. Oto Czas dodania dwóch wektorów całkowitych o długości 500 milionów.

Musimy dodać flagę -ftree-vectorize lub zmienić optymalizację na poziom -O3. (Zauważ, że -O3 wykonuje inne dodatkowe optymalizacje również). Flaga -fopt-info-vec jest 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. Poniżej przedstawiono timingi porównujące zarówno wersje bez wektoryzowane, jak i wektoryzowane 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ące informacje do wektoryzacji. Moglibyśmy użyć Specyfikacja OpenMP do programowania równoległego , która również dostarcza dyrektywę kompilatora simd do instruowania kompilatorów do wektoryzacji kodu. Ważne jest, aby upewnić się, że nie ma nakładania się pamięci, warunków wyścigu itp.. podczas ręcznego wektoryzowania kodu, 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, zwracamy się do kompilatora o wektoryzację bez względu na wszystko. Będziemy musieli aktywować rozszerzenia OpenMP używając znacznika 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

W tym sensie, nawet jeśli Strona Wikipedii dotycząca programowania tablic wspomina, że języki operujące na całych tablicach zwykle nazywają to operacjamiwektoryzowanymi , to tak naprawdę jest toukrywanie pętli IMO (chyba że jest to faktycznie wektoryzowane).

W przypadku R, nawet rowSums() lub colSums() kod w C nie wykorzystuje kod wektoryzacji IUC; jest to tylko pętla w C. to samo dotyczy lapply(). W przypadku apply(), jest w R. wszystkie z nich są zatem ukrywaniem pętli.

W skrócie, zawijanie funkcji R przez:

Wystarczy napisać for-loop W C!= wektoryzacja kodu.
wystarczy napisać for-loop w R!= wektoryzacja kodu.

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

HTH

Bibliografia:

1. W 2007 roku, w ramach programu "Horyzont 2020", w ramach programu "Horyzont 2020", w ramach programu "Horyzont 2020" - programu ramowego w zakresie badań naukowych i innowacji (2014-2020), zrealizowano]}

Więc sumując świetne odpowiedzi/komentarze do jakiejś ogólnej odpowiedzi i podaj trochę tła: R ma 4 rodzaje pętli (W porządku nie-wektoryzowanym do wektoryzowanego )

1. R for pętla, która wielokrotnie wywołuje funkcje R w każdej iteracji (Nie wektoryzowana )
Pętla C, która wielokrotnie wywołuje funkcje R w każdej iteracji (Nie wektoryzowana )
2. pętla C wywołująca funkcję r tylko raz (nieco wektoryzowana )
3. A prosta pętla C, która nie wywołuje żadnej r funkcji w ogóle i używa jej własnych skompilowanych funkcji (Wektoryzowanych)

Rodzina *apply jest więc drugim rodzajem. Z wyjątkiem apply, który jest bardziej pierwszego typu

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

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

/* to jest coś specjalnego .Wewnętrzne, tak samo bezcenne argumenty. Jest to
wywołany z zamknięcia wrapper, więc X i FUN to obietnice. Zabawa musi być bezcenne do wykorzystania np. w bquote . */

Oznacza to, że lapplykod S C akceptuje nieocenioną funkcję z R, a następnie ocenia ją w samym kodzie C. Jest to w zasadzie różnica pomiędzy lapplys .Internal wywołaniem

.Internal(lapply(X, FUN))

Który ma FUN argument, który posiada funkcję R

I colMeans .Internal wywołanie, które nie ma } argumentu FUN

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

colMeans, W przeciwieństwie do lapply wie dokładnie jakiej funkcji potrzebuje użyć, dlatego oblicza średnią wewnątrz kodu C.

Można wyraźnie zobaczyć proces oceny funkcji R w każdej iteracji wewnątrz lapply kod 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 wszystko, lapply nie jest wektoryzowana, choć ma dwie możliwe przewagi nad zwykłą pętlą R for

1. Dostęp i przypisywanie w pętli wydaje się być szybsze w C (tj. w lapply ing a funkcja) chociaż różnica wydaje się duża, to jednak pozostajemy na poziomie mikrosekundy, a kosztowną rzeczą 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 wspomniano przez @ Roland, uruchamia skompilowaną pętlę C, a raczej zinterpretowaną pętlę R

Chociaż podczas wektoryzacji kodu, są pewne rzeczy, które musisz wziąć pod uwagę.

1. Jeśli twój zbiór danych (nazwijmy go df) należy do klasy data.frame, niektóre funkcje wektorowe (takie jakcolMeans, colSums, rowSums, itd.) będzie musiał najpierw przekształcić go w matrycę, po prostu dlatego, że tak zostały zaprojektowane. Oznacza to, że dla dużego df może to spowodować ogromne obciążenie. Chociaż lapply nie będzie musiał tego robić, ponieważ wyciąga rzeczywiste wektory z df (ponieważ data.frame jest tylko listą wektorów) i dlatego, jeśli masz nie tyle kolumn, ale wiele wierszy, lapply(df, mean) może być czasami lepszą opcją niż colMeans(df).
2. kolejną rzeczą do zapamiętania jest to, że R ma duża różnorodność różnych typów funkcji, takich jak .Primitive, oraz generyczne (S3, S4) zobacz tutaj aby uzyskać dodatkowe informacje. Funkcja generyczna musi wykonać wysyłkę metody, która czasami jest kosztowną operacją. Na przykład, mean jest ogólną funkcją S3, podczas gdy sum jest Primitive. Tak więc czasami lapply(df, sum) może być bardzo wydajne w porównaniu colSums z powodów wymienionych powyżej