vectorize a simple for loop?

  • Follow


Does anyone know how to vectorize this simple for loop? 

% dum = [1 3; 1 4; 2 4]
dum =

     1     3
     1     4
     2     4

% Dij = zeros(size(dum,1),4)
Dij =

     0     0     0     0
     0     0     0     0
     0     0     0     0

for j=1:size(dum,1)
    for k = 1:size(dum,2)
        Dij(j,dum(j,k)) = 1;
    end
end

% Dij      
Dij =

     1     0     1     0
     1     0     0     1
     0     1     0     1

I think I want to use something like:
% Dij = zeros(size(dum,1),4);
% Dij(dum)=1

Dij =

     1     1     0     0
     1     0     0     0
     1     0     0     0

but obviously this does not work...

Thanks,
StevieP
0
Reply Stevie 12/11/2009 1:20:23 AM

One approach:

% Data
dum = [1 3; 1 4; 2 4];
Dij = zeros(size(dum,1),4);

% Engine
Dij(sub2ind(size(Dij),cumsum(ones(size(dum))),dum)) = 1

You could also inline sub2ind for more speed.
0
Reply Matt 12/11/2009 3:17:21 AM


Hi Steve,

You can easily verify that the (m,n)th element of a matrix corresponds to that matrix's 

(n-1)*M+m th

element.

For example, if

A = [1 2 3 4 5;
       6 7 8 9 10;
       11 12 13 14 15;
       16 17 18 19 20;
       21 22 23 24 25];

then, A(3,2) = 12, which is equal to A(k) where k = (2-1)*size(A,1)+3.

We can use the same trick. In fact, you can see that the values in dum are all column indices. The first row is for the first row of Dij, the second is for the second and so on. So we can define:

dum2 = (dum - 1)*size(dum,1)+diag(1:size(dum,1))*ones(size(dum));

The first term in the sum is obvious: (n-1)*M. The second term does nothing but determines m in (n-1)*M+m. 

So the answer to your question is

dum = [1 3; 1 4; 2 4];
dum2 = (dum - 1)*size(dum,1)+diag(1:size(dum,1))*ones(size(dum));
Dij = zeros(size(dum,1),4);
Dij(dum2(:)) = 1;

gives what you are looking for.

Best.
0
Reply Sadik 12/11/2009 3:24:19 AM

"Sadik " <sadik.hava@gmail.com> wrote in message <hfse13$fm8$1@fred.mathworks.com>...
> Hi Steve,
> 
> You can easily verify that the (m,n)th element of a matrix corresponds to that matrix's 
> (n-1)*M+m th
> element.
....

Dear Matt Fig and Sadik, 

Very clever answers! I see that your methods were essentially the same, differing essentially only by explicit or implicit use of the function sub2ind, so far as I could tell. I ran them both on my matrices and it seems that explicitly defining  
dij2 = (dij - 1)*size(dij,1)+diag(1:size(dij,1))*ones(size(dij))
does make it slightly faster. I think this was what Matt Fig was referring to when he mentioned sub2ind could be made faster by inline function. 

Strangely though, and I did not expect this for my problem which had fairly large matrices (of several hundred by several hundred), the for-loop actually tends to be faster than either vectorized method. I thought it always paid to vectorize, but apparently not... I'll have look into that more. 

In any event, thanks for the clever solutions, 

steviep
0
Reply Stevie 12/22/2009 4:50:06 AM

Another method:

clear DUM
DUM(:,:,2)=dum;
DUM(:,:,1)=ndgrid(1:size(dum,1),1:size(dum,2));
Dij = accumarray(reshape(DUM,[],2),1);

I would suggest to take a look of using sparse matrix for matrix with a lot of zeros.

Bruno
0
Reply Bruno 12/22/2009 8:07:04 AM

4 Replies
180 Views

(page loaded in 0.17 seconds)


Reply: