Julia: ottimizzazione della simulazione del sistema dinamico semplice

Question

Sto cercando di ottimizzare la simulazione di un semplice sistema dinamico in cui la risposta di una rete e dei suoi parametri (pesi) si evolvono secondo semplici equazioni lineari. La simulazione deve essere eseguita per decine di milioni di passaggi temporali, ma le dimensioni della rete saranno in genere ridotte. Quindi, le prestazioni sono limitate meno dai prodotti a matrice vettoriale, ma piuttosto da matrici temporanee, controlli vincolati e altri fattori meno visibili. Dato che sono nuovo di Julia, apprezzerei qualsiasi suggerimento per ottimizzare ulteriormente le prestazioni.

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 rete (2.000.000 passaggi) con

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

produce i tempi

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

Aggiornamento 1

Seguendo i consigli da David Sanders ho ottenuto liberare la macro di devec e scrivere i loop. Questo riduce effettivamente allocazioni di matrice e migliora le prestazioni di circa il 25%, ecco i nuovi numeri:

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

Più piccola è la dimensione della rete, maggiore è la spinta. Un esempio del codice di simulazione aggiornato può essere trovato here.

Aggiorna 2

Gran parte delle allocazioni di memoria sono dovuti ai prodotti matrice-vettore. Così, al fine di sbarazzarsi di quelli che ho sostituito i prodotti da un'operazione BLAS sul posto, BLAS.genv !, che riduce tempi da un'altra allocazione del 25% e la memoria del 90%,

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

codice here Aggiornato .

Update 3

Il più grande aggiornamento rango-1 può anche essere sostituita da due chiamate in-place funzioni BLAS, cioè BLAS.scal! per ridimensionare e BLAS.ger! per un aggiornamento di livello 1. L'avvertenza è che entrambe le chiamate sono abbastanza lento se si usa più di un filo (problema con OpenBLAS?), Quindi è meglio impostare

blas_set_num_threads(1) 

C'è un guadagno 15% in tempi per una dimensione di rete di 20 , e un guadagno del 50% per una rete di dimensioni 50. non ci sono più allocazioni di memoria, ei nuovi tempi sono

1.638287 seconds (11 allocations: 1.266 KB) 

Anche in questo caso, il codice aggiornato può essere trovato here.

Update 4

ho scritto una base Cython script per confrontare finora i risultati. La differenza principale è che non uso nessuna chiamata a BLAS ma ho loop: l'iniezione di chiamate BLAS a basso livello è un problema in Cython e le chiamate a numpy dot hanno troppe spese generali per le piccole dimensioni di rete (ho provato ...).Tempi sono

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

che è approssimativamente uguale alla versione originale (da cui, finora, il 50% viene rasato).

Answer 1

Sebbene si stia utilizzando il pacchetto Devectorize.jl, suggerisco di scrivere esplicitamente tutte quelle operazioni vettorizzate come loop semplici. Mi aspetto che questo ti darà un significativo incremento delle prestazioni.

pacchetto

Devectorize è certamente un grande contributo, ma per vedere i cerchi che salta attraverso di fare il lavoro sporco per voi, si può fare qualcosa di simile (un esempio dal README pacchetto):

using Devectorize 

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

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

Qui, macroexpand è una funzione che ti dice il codice a cui la macro @devec espande il suo argomento (il codice sul resto della riga). Non voglio rovinare la sorpresa mostrando l'output qui, ma non è solo il semplice ciclo for che si dovrebbe scrivere a mano.

Inoltre, il fatto di disporre di un'allocazione enorme suggerisce che non tutte le operazioni vettoriali vengono gestite correttamente.

A proposito, non dimenticare di eseguire prima una piccola corsa in modo da non cronometrare il passo della compilazione.

[Nota tangenziale: qui, exp è la funzione che applica la normale funzione esponenziale a ciascun elemento di una matrice, equivalente a map(exp, a+b). expm fornisce l'esponenziale di una matrice. Si è parlato di deprecare tali usi di exp.]