Come vectorize doppia sommatoria in MATLAB

Question

sto cercando di produrre codice vettorizzati per questa equazione (e molteplici altri della stessa forma):

Questo sarà valutato nel ode45 e l'unica parte che le modifiche sono Rmn(t) quindi sono pre-calcolo delle funzioni di Bessel e sin. Al momento il mio codice è simile al seguente:

for iR = 1:nR 
    p(:, iR) = sum(Rmn.*J0m(:, :, iR),1)*sinanz; 
end 

M, N sono il numero di termini nella mia sommatoria e R, Z sono il numero di r e z coordinate sto usando. Rmn è una matrice M*N. J0m è un array M*N*R. È una matrice M*R che viene ripetuta N volte. sinanz è una matrice N*Z. J0m e sinanz sono precalcolati e non cambiano.

Questo funziona ma è lento e quindi sto cercando di ottimizzarlo. Ho pensato che il primo passo sarebbe quello di ridurre J0m quindi è solo m*R ma non riesco a capire come. Sto cercando qualche suggerimento su come farlo e qualche suggerimento su come ottimizzare questo in generale.

Answer 1

Come probabilmente già sapete, si dovrebbe pre-allocare p prima del ciclo:

p = zeros(Z, nR); 

questo modo si evita l'array p di crescere ad ogni iterazione, accelerando il ciclo tremendamente.

Si potrebbe vettorizzare il tutto in virtù della bsxfun:

% C ≣ M×N×R array of all products Rmn·J0m 
C = bsxfun(@times, Rmn, J0m); 

% D ≣ M×N×R×Z array of all products C·sinanz 
D = bsxfun(@times, C, permute(sinanz, [3 1 4 2])); 

% Sum in M and N directions, and re-order 
p = permute(sum(sum(D,1),2), [4 3 1 2]); 

ma dubito che sarà più veloce; MATLAB (leggi: BLAS) è piuttosto veloce con le matrici 2D, ma di solito non è molto bello con gli array di più D.

Ti suggerisco di leggere su bsxfun; è anche il modo per ridurre l'array J0m su M × R nel modo in cui descrivi.

Naturalmente, si può sbarazzarsi dei permute s definendo correttamente le variabili, quindi cerchiamo di fare un piccolo test nelle versioni 'ideali' sia del codice di loop e vettorializzare:

%% Initialize some variables 

% Number of tests 
TESTS = 1e4; 

% Your dimensions 
M = 5; nR = 4; 
N = 2; Z = 3; 

% Some dummy data 
Rmn = rand(M,N); 
sinanz = rand(N,Z); 
J0m = rand(M,nR); % NOTE: J0m doesn't have to be 3D anymore by using bsxfun 

%% Test 1: your own version, optimized 

% Start test 
start = tic; 

% Repeat the test a couple of times to get a good average 
for ii = 1:TESTS 
    p1 = zeros(Z,nR); 
    for iR = 1:nR 
     p1(:, iR) = sum(bsxfun(@times,Rmn,J0m(:, iR)), 1)*sinanz; 
    end  
end 
stop = toc(start); 

% Average execution time 
fprintf(1, 'Average time of looped version: %f seconds.\n', stop/TESTS); 

%% Vectorized version, using 4D arrays: 

% Don't make the permutes part of the test 
J0m = permute(J0m, [1 3 2]); 
sinanz = permute(sinanz, [3 1 4 2]); 

% Start test 
start = tic; 

% Repeat the test a couple of times to get a good average 
for ii = 1:TESTS 
    C = bsxfun(@times, Rmn, J0m); 
    D = bsxfun(@times, C, sinanz); 
    p2 = sum(sum(D,1),2); 

end 
stop = toc(start); 

% Average execution time 
fprintf(1, 'Average time of vectorized version: %f seconds.\n', stop/TESTS); 

%% Check for equality 

maxDifference = max(max(p1 - squeeze(p2)')) 

risultati:

Average time of looped version : 0.000054 seconds. 
Average time of vectorized version: 0.000023 seconds. 
maxDifference = 
    4.440892098500626e-016 

che sembra abbastanza buono! Tuttavia, utilizzando

M = 50; nR = 40; 
N = 20; Z = 30; 

si ottiene

Average time of looped version : 0.000708 seconds. 
Average time of vectorized version: 0.009835 seconds. 
maxDifference = 
    8.526512829121202e-014 

così la versione vettorializzare ottenuto appena un ordine di grandezza più lento rispetto alla variante in loop.

Ovviamente, your mileage may vary, ma il messaggio da portare a casa è che dovresti aspettarti che questa differenza peggiori sempre di più per aumentare le dimensioni.

Quindi, in conclusione:

• se le tue dimensioni sono di grandi dimensioni, il bastone con il ciclo. Se sono piccoli, utilizzare anche il ciclo - è solo SO molto più facile per gli occhi rispetto alla versione vettorializzare :)
• assicurarsi di pre-allocare il p variabile
• uso bsxfun per ridurre l'occupazione di memoria (ma non aumenterà molto la velocità)