arrayfun może być znacznie wolniejsza niż jawna pętla w programie Matlab. Czemu?

Rozważ następujący prosty test szybkości dla arrayfun:

T = 4000;
N = 500;
x = randn(T, N);
Func1 = @(a) (3*a^2 + 2*a - 1);

tic
Soln1 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln1(t, n) = Func1(x(t, n));
    end
end
toc

tic
Soln2 = arrayfun(Func1, x);
toc

Na moim komputerze (Matlab 2011b na Linux Mint 12) wynik tego testu to:

Elapsed time is 1.020689 seconds.
Elapsed time is 9.248388 seconds.

Co?!? arrayfun, choć wprawdzie bardziej przejrzyste rozwiązanie, jest o rząd wielkości wolniejsze. Co tu się dzieje?

Ponadto wykonałem podobny test dla cellfuni stwierdziłem, że jest około 3 razy wolniejszy niż jawna pętla. Ponownie, ten wynik jest przeciwieństwem tego, czego się spodziewałem.

Moje pytanie brzmi: dlaczego są arrayfunio cellfunwiele wolniejsze? A biorąc pod uwagę to, czy są jakieś dobre powody, aby ich używać (poza tym, że kod wygląda dobrze)?

Uwaga: mówię tutaj o standardowej wersji arrayfun, a NIE wersji GPU z zestawu narzędzi do przetwarzania równoległego.

EDYCJA: Żeby było jasne, zdaję sobie sprawę, że Func1powyżej można wektoryzować, jak wskazał Oli. Wybrałem to tylko dlatego, że daje prosty test szybkości dla celów rzeczywistego pytania.

EDYCJA: Zgodnie z sugestią grungetty ponownie wykonałem test feature accel off. Wyniki są następujące:

Elapsed time is 28.183422 seconds.
Elapsed time is 23.525251 seconds.

Innymi słowy, wydaje się, że duża część różnicy polega na tym, że akcelerator JIT znacznie lepiej przyspiesza jawną forpętlę niż robi arrayfun. Wydaje mi się to dziwne, ponieważ w arrayfunrzeczywistości dostarcza więcej informacji, tj. Jego użycie pokazuje, że kolejność wywołań Func1nie ma znaczenia. Zauważyłem również, że niezależnie od tego, czy akcelerator JIT jest włączony, czy wyłączony, mój system używa tylko jednego procesora ...

Możesz zrozumieć pomysł, uruchamiając inne wersje swojego kodu. Rozważ jawne napisanie obliczeń, zamiast używania funkcji w pętli

tic
Soln3 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln3(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

Czas na obliczenia na moim komputerze:

Soln1  1.158446 seconds.
Soln2  10.392475 seconds.
Soln3  0.239023 seconds.
Oli    0.010672 seconds.

Teraz, podczas gdy w pełni „wektoryzowane” rozwiązanie jest zdecydowanie najszybsze, widać, że zdefiniowanie funkcji, która ma być wywoływana dla każdego wpisu x, jest ogromnym narzutem. Samo wyraźne wypisanie obliczeń dało nam 5-krotne przyspieszenie. Wydaje mi się, że to pokazuje, że kompilator MATLABs JIT nie obsługuje funkcji wbudowanych . Zgodnie z odpowiedzią udzieloną tam przez gnovice, właściwie lepiej jest napisać normalną funkcję niż anonimową. Spróbuj.

Następny krok - usunięcie (wektoryzacja) pętli wewnętrznej:

tic
Soln4 = ones(T, N);
for t = 1:T
    Soln4(t, :) = 3*x(t, :).^2 + 2*x(t, :) - 1;
end
toc

Soln4  0.053926 seconds.

Kolejne przyspieszenie o czynnik 5: w tych stwierdzeniach jest coś, co mówi, że powinieneś unikać pętli w MATLAB-ie ... A może naprawdę? Spójrz na to więc

tic
Soln5 = ones(T, N);
for n = 1:N
    Soln5(:, n) = 3*x(:, n).^2 + 2*x(:, n) - 1;
end
toc

Soln5   0.013875 seconds.

Znacznie bliżej do „w pełni” wektoryzowanej wersji. Matlab przechowuje macierze według kolumn. Zawsze (jeśli to możliwe) należy nadawać obliczeniom strukturę wektoryzacji „kolumnowej”.

Możemy teraz wrócić do Soln3. Tam kolejność pętli jest „wierszowa”. Zmieńmy to

tic
Soln6 = ones(T, N);
for n = 1:N
    for t = 1:T
        Soln6(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

Soln6  0.201661 seconds.

Lepiej, ale nadal bardzo źle. Pojedyncza pętla - dobra. Podwójna pętla - źle. Myślę, że MATLAB wykonał porządną pracę nad poprawą wydajności pętli, ale nadal istnieje obciążenie pętli. Gdybyś miał w środku cięższą pracę, nie zauważyłbyś. Ale ponieważ to obliczenie jest ograniczone pasmem pamięci, widać narzut pętli. I będzie jeszcze wyraźniej zobaczyć napowietrznej wywołanie func1 tam.

Więc o co chodzi z arrayfun? Tam też nie ma żadnej funkcji, więc dużo narzutów. Ale dlaczego jest o wiele gorszy niż podwójna zagnieżdżona pętla? Właściwie temat korzystania z cellfun / arrayfun był obszernie omawiany wiele razy (np. Tutaj , tutaj , tutaj i tutaj ). Te funkcje są po prostu powolne, nie można ich używać do tak drobnoziarnistych obliczeń. Możesz ich używać do zwięzłości kodu i fantazyjnych konwersji między komórkami i tablicami. Ale funkcja musi być cięższa niż ta, którą napisałeś:

tic
Soln7 = arrayfun(@(a)(3*x(:,a).^2 + 2*x(:,a) - 1), 1:N, 'UniformOutput', false);
toc

Soln7  0.016786 seconds.

Zauważ, że Soln7 jest teraz komórką… czasami jest to przydatne. Wydajność kodu jest teraz całkiem dobra i jeśli potrzebujesz komórki jako wyniku, nie musisz konwertować macierzy po użyciu w pełni wektoryzowanego rozwiązania.

Dlaczego więc arrayfun działa wolniej niż prosta struktura pętli? Niestety nie możemy tego powiedzieć na pewno, ponieważ nie ma dostępnego kodu źródłowego. Można się tylko domyślać, że skoro arrayfun jest funkcją ogólnego przeznaczenia, która obsługuje wszystkie rodzaje różnych struktur danych i argumentów, niekoniecznie jest bardzo szybka w prostych przypadkach, które można bezpośrednio wyrazić jako zagnieżdżenia pętli. Skąd się bierze narzut, nie wiemy. Czy można uniknąć kosztów ogólnych dzięki lepszej implementacji? Może nie. Ale niestety jedyne, co możemy zrobić, to zbadać wydajność, aby zidentyfikować przypadki, w których działa dobrze, i te, w których nie.

Aktualizacja Ponieważ czas wykonania tego testu jest krótki, aby uzyskać wiarygodne wyniki dodałem teraz pętlę wokół testów:

for i=1:1000
   % compute
end

Czasami podane poniżej:

Soln5   8.192912 seconds.
Soln7  13.419675 seconds.
Oli     8.089113 seconds.

Widzisz, że arrayfun jest nadal zły, ale przynajmniej nie o trzy rzędy wielkości gorszy niż rozwiązanie wektoryzowane. Z drugiej strony pojedyncza pętla z obliczeniami opartymi na kolumnach jest tak szybka, jak w pełni zwektoryzowana wersja ... Wszystko to zostało zrobione na jednym procesorze. Wyniki dla Soln5 i Soln7 nie zmieniają się, jeśli przełączę się na 2 rdzenie - w Soln5 musiałbym użyć parfor, aby uzyskać równoległość. Zapomnij o przyspieszeniu ... Soln7 nie działa równolegle, ponieważ arrayfun nie działa równolegle. Z drugiej strony wersja wektoryzowana Olis:

Oli  5.508085 seconds.

To dlatego, że !!!!

x = randn(T, N); 

nie jest gpuarraytypem;

Wszystko, co musisz zrobić, to

x = randn(T, N,'gpuArray');