Julia: optymalizacja symulacji prostego systemu dynamicznego

Question

Próbuję zoptymalizować symulację prostego systemu dynamicznego, w którym odpowiedź sieci, a także jej parametry (wagi) ewoluują zgodnie z prostymi równaniami liniowymi. Symulacja musi trwać dziesiątki milionów kroków, ale rozmiar sieci zwykle będzie mały. W związku z tym wydajność jest ograniczona w mniejszym stopniu przez produkty matrycowo-wektorowe, ale raczej przez tymczasowe tablice, powiązane kontrole i inne mniej widoczne czynniki. Ponieważ jestem nowy w Julia, doceniam wszelkie wskazówki, które pomogą zoptymalizować wydajność.

function train_network(A, T, Of, cs, dt) 
    N, I = size(T) 
    z = zeros(I) 
    r = zeros(N) 

    @inbounds for t in 1:size(cs, 1) 
     # precompute 
     Az = A*z 
     Ofr = Of*r 

     # compute training signal 
     @devec z += dt.*(Az + cs[t] - 0.5.*z) 
     I_teach = T*(Az + cs[t]) 
     Tz  = T*z 

     # rate updates 
     @devec r += dt.*(I_teach - Ofr - 0.1.*r) 

     # weight updates 
     for i in 1:I 
      @devec T[:, i] += dt.*1e-3.*(z[i].*r - T[:, i]) 
     end 

     for n in 1:N 
      @devec Of[:, n] += dt.*1e-3.*(Tz.*r[n] - Of[:, n])  
     end 
    end 
end 

# init parameters 
N, I = 20, 2 
dt = 1e-3 

# init weights 
T = rand(N, I)*N 
A = rand(I, I) 
Of = rand(N, N)/N 

# simulation time & input 
sim_T = 2000 
ts = 0:dt:sim_T 
cs = randn(size(ts, 1), I) 

Timing sieci (2.000.000 stopni) z

@time train_network(A, T, Of, cs, dt) 

daje taktowania

3.420486 seconds (26.12 M allocations: 2.299 GB, 6.65% gc time) 

Aktualizacja 1

Idąc za radą David Sanders dostałam pozbyć się makra devec i wypisać pętle. To istotnie zmniejsza alokacji tablicy i zwiększa wydajność o około 25%, tutaj są nowe numery:

2.648113 seconds (18.00 M allocations: 1.669 GB, 5.60% gc time) 

Im mniejszy rozmiar sieci, tym większy impuls. Istotą zaktualizowanego kodu symulacji można znaleźć here.

Aktualizacja 2

wiele alokacji pamięci, są ze względu na produkty matrycy wektora. Tak, aby pozbyć się tych Wymieniłem te produkty przez działanie w miejscu Blas, BLAS.genv !, która zmniejsza czasy przez kolejne 25% i alokacji pamięci o 90%,

1.990031 seconds (2.00 M allocations: 152.589 MB, 0.69% gc time) 

Zaktualizowany kod here .

Update 3

Największy stopień-1 Aktualizację można również zastąpić dwoma połączeniami w miejscu funkcji Blas, mianowicie BLAS.scal! do skalowania i BLAS.ger! dla aktualizacji klasy 1. Zastrzeżenie jest to, że oba połączenia są dość powolne, jeśli więcej niż jeden wątek jest używany (problem z OpenBLAS?), Więc najlepiej jest ustawić

blas_set_num_threads(1) 

Jest przyrost taktowania do wielkości sieci 20 15% oraz zysk w wysokości 50% dla sieci o rozmiarze 50. nie ma więcej alokacje pamięci, a nowe czasy są

1.638287 seconds (11 allocations: 1.266 KB) 

Znowu zaktualizowany kod można znaleźć here.

Update 4

pisałem podstawową Cython script porównanie wyników tej pory. Główną różnicą jest to, że nie używam żadnych wywołań do BLAS, ale mam pętle: Wstrzykiwanie niskopoziomowych wywołań BLAS jest bólem w Cythonie, a wywołania numpy dot mają zbyt dużo narzutów dla małych rozmiarów sieci (próbowałem ...).Parametry czasowe to

CPU times: user 3.46 s, sys: 6 ms, total: 3.47 s, Wall time: 3.47 s 

, która jest mniej więcej taka sama jak wersja oryginalna (z której, jak dotychczas, goli zostało 50%).

Answer 1

Chociaż korzystasz z pakietu Devectorize.jl, proponuję po prostu zapisać wszystkie wektoryzowane operacje jawnie jako proste pętle. Spodziewam się, że da to znaczący wzrost wydajności.

Devectorize pakiet jest z pewnością wielki wkład, ale aby zobaczyć obręcze to przeskakuje do wykonania brudnej roboty dla ciebie, możesz zrobić coś takiego (przykład z README pakietu):

using Devectorize 

a = rand(2,2); 
b = rand(2,2); 
c = rand(2,2); 

julia> macroexpand(:(@devec r = exp(a + b) .* sum(c))) 

Tutaj, macroexpand jest funkcją, która mówi ci kod, do którego makro @devec rozwija swój argument (kod na pozostałej części linii). Nie zepsuję niespodzianki, pokazując tutaj wyjście, ale to nie jest zwykła pętla for, którą napiszesz ręcznie.

Co więcej, fakt, że masz dużą alokację sugeruje, że nie wszystkie operacje wektorowe są poprawnie obsługiwane.

Przy okazji, nie zapomnij zrobić małego biegu jako pierwszego, więc nie masz czasu na krok kompilacji.

[Uwaga styczna: tutaj, exp to funkcja, która stosuje zwykłą funkcję wykładniczą do każdego elementu macierzy, równoważną map(exp, a+b). expm podaje wykładniczy macierz. Mówi się o wycofywaniu takich zastosowań z exp.]