Re: Matrix loop

[Jaroslav Hajek]

On Wed, May 5, 2010 at 6:27 PM, Riccardo Corradini <address@hidden> wrote:

Hi,
I would know if is there a more efficient way to perform the following operation
s_cR =21;
n=1000;
a= rand(s_cR,s_cR,n);
whos a
meant = mean(a,3);
whos meant
for tt=1:n
    b(:,:,tt) = a(:,:,tt) - meant;
    c(:,:,tt)= b(:,:,tt) * b(:,:,tt)';
endfor
result = sum(c,3)*1/n;
whos result

Thanks in advance

Lots of options. But this particular code's performance sucks primarily because you reallocate arrays b and c at each step:

octave:1> s_cR =21;
octave:2> n=1000;
octave:3> a= rand(s_cR,s_cR,n);
octave:4> tic;
octave:5> meant = mean(a,3);
octave:6> for tt=1:n
>     b(:,:,tt) = a(:,:,tt) - meant;
>     c(:,:,tt)= b(:,:,tt) * b(:,:,tt)';
> endfor
octave:7> result = sum(c,3)*1/n;
octave:8> toc
Elapsed time is 1.96773 seconds.

pre-initializing b and c to proper size gives a big speed-up:

octave:1> s_cR =21;
octave:2> n=1000;
octave:3> a= rand(s_cR,s_cR,n);
octave:4> tic;
octave:5> b = c = zeros (size (a));
octave:6> meant = mean(a,3);
octave:7> for tt=1:n
>     b(:,:,tt) = a(:,:,tt) - meant;
>     c(:,:,tt)= b(:,:,tt) * b(:,:,tt)';
> endfor
octave:8> result = sum(c,3)*1/n;
octave:9> toc
Elapsed time is 0.0619521 seconds.

yet simpler and faster:

octave:18> tic;
octave:19> b = center (a, 3);
octave:20> c = cellfun (@(x) x*x', num2cell (b, [1, 2]), "uniformoutput", false);
octave:21> result = plus (c{:}) / n;
octave:22> toc
Elapsed time is 0.027951 seconds.

hth

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz