arrayfun può essere notevolmente più lento di un ciclo esplicito in MATLAB. Perché?

Considerare il seguente test di velocità semplice per 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 

Sulla mia macchina (Matlab 2011b su Linux Mint 12), l’output di questo test è:

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

Che cosa?!? arrayfun , pur ammettendo una soluzione dall’aspetto più pulito, è un ordine di grandezza più lento. Che cosa sta succedendo qui?

Inoltre, ho fatto uno stile simile di test per cellfun e cellfun trovato che era circa 3 volte più lento di un ciclo esplicito. Ancora una volta, questo risultato è l’opposto di quello che mi aspettavo.

La mia domanda è: perché arrayfun e cellfun molto più lenti? E dato questo, ci sono dei buoni motivi per usarli (oltre a far sembrare il codice buono)?

Nota: sto parlando della versione standard di arrayfun qui, NON della versione GPU dalla toolbox di elaborazione parallela.

EDIT: Solo per essere chiari, sono consapevole che Func1 sopra può essere vettorizzato come indicato da Oli. L’ho scelto solo perché produce un semplice test di velocità ai fini della domanda reale.

EDIT: Seguendo il suggerimento di Grungetta, ho rifatto il test con la feature accel off . I risultati sono:

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

In altre parole, sembrerebbe che gran parte della differenza sia che l’acceleratore JIT fa un lavoro molto migliore di accelerare il ciclo esplicito for quello che fa arrayfun . Questo mi sembra strano, dato che arrayfun fornisce effettivamente più informazioni, cioè il suo uso rivela che l’ordine delle chiamate a Func1 non ha importanza. Inoltre, ho notato che se l’acceleratore JIT è acceso o spento, il mio sistema usa sempre solo una CPU …

Puoi ottenere l’idea eseguendo altre versioni del tuo codice. Considera esplicitamente di scrivere i calcoli, invece di usare una funzione nel tuo ciclo

 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 

Tempo di calcolo sul mio computer:

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

Ora, mentre la soluzione completamente “vettorializzata” è chiaramente la più veloce, è ansible vedere che definire una funzione da chiamare per ogni voce x è un enorme sovraccarico. Solo scrivere esplicitamente il calcolo ci ha permesso di accelerare il fattore 5. Immagino che questo dimostri che il compilatore JIT MATLAB non supporta le funzioni inline . Secondo la risposta di gnovice, in realtà è meglio scrivere una funzione normale piuttosto che anonima. Provalo.

Passaggio successivo: rimuovi (vectorizza) il ciclo interno:

 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. 

Un altro fattore di 5 accelerazione: c’è qualcosa in quelle affermazioni che dicono che dovresti evitare i loop in MATLAB … O c’è davvero? Dai un’occhiata a questo allora

 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. 

Molto più vicino alla versione “completamente” vettorizzata. Matlab memorizza le matrici in base alla colonna. Dovresti sempre (quando ansible) strutturare i tuoi calcoli per essere vettorizzati ‘in termini di colonne’.

Possiamo tornare a Soln3 ora. L’ordine del ciclo è ‘riga-saggio’. Consente di cambiarlo

 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. 

Meglio, ma ancora molto male. Ciclo singolo – buono. Ciclo doppio – cattivo. Immagino che MATLAB abbia fatto del lavoro decente per migliorare le prestazioni dei loop, ma il sovraccarico del loop è ancora lì. Se avessi un lavoro più pesante dentro, non te ne accorgeresti. Ma dal momento che questo calcolo è limitato dalla larghezza di banda della memoria, si vede il sovraccarico del ciclo. E vedrai ancora più chiaramente l’overhead di chiamare Func1 lì.

Allora, che succede con arrayfun? Nessuna funzione in linea neanche lì, quindi un sacco di spese generali. Ma perché è molto peggio di un doppio ciclo annidato? In realtà, l’argomento dell’utilizzo di cellfun / arrayfun è stato ampiamente discusso molte volte (ad es. Qui , qui , qui e qui ). Queste funzioni sono semplicemente lente, non è ansible utilizzarle per calcoli a grana fine. Puoi usarli per brevità del codice e conversioni elaborate tra celle e array. Ma la funzione deve essere più pesante di quello che hai scritto:

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

Nota che Soln7 è una cella ora .. a volte è utile. Le prestazioni del codice sono piuttosto buone ora, e se hai bisogno di una cella come output, non è necessario convertire la matrice dopo aver utilizzato la soluzione completamente vettoriale.

Quindi, perché arrayfun è più lento di una semplice struttura ad anello? Sfortunatamente, è imansible per noi dirlo con certezza, poiché non esiste un codice sorgente disponibile. Si può solo supporre che poiché arrayfun è una funzione generica, che gestisce tutti i tipi di strutture e argomenti di dati diversi, non è necessariamente molto veloce in casi semplici, che è ansible esprimere direttamente come nidi di loop. Da dove viene il sovraccarico non possiamo sapere. L’overhead potrebbe essere evitato con una migliore implementazione? Forse no. Ma sfortunatamente l’unica cosa che possiamo fare è studiare la performance per identificare i casi, in cui funziona bene, e quelli, dove non funziona.

Aggiornamento Poiché il tempo di esecuzione di questo test è breve, per ottenere risultati affidabili ho aggiunto ora un ciclo attorno ai test:

 for i=1:1000 % compute end 

Alcune volte indicato di seguito:

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

Si vede che l’arrayfun è ancora cattivo, ma almeno non tre ordini di grandezza peggiore della soluzione vettoriale. D’altra parte, un singolo ciclo con calcoli basati su colonne è veloce quanto la versione completamente vettoriale … Tutto ciò è stato fatto su una singola CPU. I risultati per Soln5 e Soln7 non cambiano se passo a 2 core – In Soln5 dovrei usare un parfor per farlo parallelizzare. Dimentica la velocità … Soln7 non funziona in parallelo perché arrayfun non funziona in parallelo. Dall’altra parte la versione vettoriale di Olis:

 Oli 5.508085 seconds. 

Questo perché!!!!

 x = randn(T, N); 

non è gpuarray tipo gpuarray ;

Tutto quello che devi fare è

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