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)
|