How to vectorize this code?

Hi all,

I need to vectorize the following code in order to make it faster:

ar=[    115  139  168;
        133  160  193;
        155  187  226];

at=[    803  667     552;
        3212  2667  2208;
        5621  4667  3864];       

afr(:,:,1)=[0  17  3;
            4  3208  29;
            6  127  12];

afr(:,:,2)= [0  0  0;
            0  14  2;
            0  3  2];

afr(:,:,3)=[0  0  0;
            0  4  1;
            0  1  0];

here my code with nested loop:

[Sar Iar]=sort(reshape(ar,1,prod(size(ar))));
[Sat Iat]=sort(reshape(at,1,prod(size(at))));

Output=zeros(prod(size(ar)),prod(size(at)));
for i=1:length(afr(1,1,:))
    for j=1:length(afr(1,:,1))
        for k=1:length(afr(:,1,1))
            Output(find(Sat==at(k,i)),find(Sar==ar(j,i)))=afr(k,j,i)
        end
    end
end

afr(:,:,i) is linked to the i_th column of ar and at. I need to create the a matrix of the following size [prod(size(at)),prod(size(ar))] where the elements of afr are sorted according to ar along the column and to at along the row. Here the result with the aforementioned data:

Ris=[0  0       0  0  0  0  0  0  0;
    0          0       0  0  0  0  0  0  0;
    0          17      0  3  0  0  0  0  0;
    0          0       0  0  0  0  0  4  1;
    0          0       0  0  14  0  2  0  0;
    4          3208  0  29  0  0  0  0  0;
    0          0       0  0  0  0  0  1  0;
    0          0       0  0  3  0  2  0  0;
    6          127     0  12  0  0  0  0  0];

how can I make it faster?

Thanks

[Geoff Hayes]

Hi pietro,

Why does the code need to be faster as it runs pretty quick as is. Are you anticipating larger matrices? Note that you could put the find(Sar==ar(j,i)) outside of the k (second inner) loop since ar does not depend on k,just i and j.

Geoff

[pietro]

Hi Geoff,

Thanks for your reply. The code must be faster because I work with bigger matrices around [350x430] and I need to perform this operation for more than 100 matrices and the computation is quite slow. I posted that matrices just as example.

Pietro

Answer 1

[Jan]

Some standard methods, which improve the code and make it more Matlabish. But the acceleration will not be high:

[Sar, Iar] = sort(ar(:));
[Sat, Iat] = sort(at(:));

Output = zeros(numel(ar), numel(at));
[s1, s2, s3] = size(afr);
for i = 1:s1
    for j = 1:s2
        v = (Sar==ar(j,i));
        for k = 1:s3
            Output(Sat==at(k,i), v) = afr(k,j,i);
        end
    end
end

With your tiny test data set I get for 1000 iterations under Matlab 2009a:

 Original: Elapsed time is 0.193134 seconds.
 Improved: Elapsed time is 0.079364 seconds.

The main ideas are the logical indexing (avoid the find()) and moving repeated calculations out of the loop (Sar==ar(j,i) does not depend on the inner loop).

length(afr(1,1,:)) creates a temporary vector at first. But this is a waste of time, because size replies the dimensions directly.