Vectorization of 2 for loops in MATLAB

I have just started exploring world of vectorization. I got the 1-D vectorization down but i am having trouble vectorizing the following code.

CityPairs = [7 3
3 1
3 1
1 7
7 1
3 4
5 1
4 6];
Offices = [1;3;7];
nOffices = size(Offices,1);

connection = zeros(nOffices);
for i = 1:nOffices
    for j = 1:nOffices
        connection(i,j) = sum(Offices(i) == CityPairs(:,1)...
            & CityPairs(:,2) == Offices(j));
    end
    connection(i,:) = connection(i,:);
end   
disp(connection)

In this example there are 7 cities, three of which have offices. I want a pairwise matrix for the cities with offices to capture the sum of all the one way connections between each. The answer for above problem should be:

   0     0     1
   2     0     0
   1     1     0

Any suggestions are welcome. Thanks in advance.

Keith

Answer 1

[Cedric]

It is not trivial to "vectorize" this setup, as there are operations that require some look up table, and there is accumulation (unless you find a trick). This is likely to make the approach without FOR loops more complicated (code-wise) than the basic, loop-based approach. So let's start with the trick first ;-) ..

 >> A = double( bsxfun(@eq, CityPairs(:,1), Offices.') ) ;
 >> B = double( bsxfun(@eq, CityPairs(:,2), Offices.') ) ;
 >> A.' * B
 ans =
     0     0     1
     2     0     0
     1     1     0

I let you study a bit BSXFUN and understand how this solution works. If you don't understand after having spent some time, ask me for details.

Now if I hadn't found a trick, I would have done it as follows (I display the output at each step to make it more understandable)..

First, we find valid city IDs (those which have an office):

 >> isValid = ismember(CityPairs, Offices)
 isValid =
     1     1
     1     1
     1     1
     1     1
     1     1
     1     0
     0     1
     0     0

Then we extract rows of CityPairs which have two valid nodes/cities:

 >> validPairs = CityPairs(all(isValid, 2),:)
 validPairs =
     7     3
     3     1
     3     1
     1     7
     7     1

We want to build a list of office IDs that correspond to valid pairs of cities, and for that we build a look up table:

 >> LT(Offices) = 1 : nOffices
 LT =
     1     0     2     0     0     0     3

You can see that LT(7) is 3, which is the office ID of city 7. With that, we can build the list of office IDs which correspond to valid city pairs:

 >> officeIds = LT(validPairs)
 officeIds =

     3     2
     2     1
     2     1
     1     3
     3     1

The last step is to accumulate a count of 1 for each row of the above array, which actually counts entries with same indices:

 >> accumarray(officeIds, ones(size(officeIds,1),1), [nOffices, nOffices])
 ans =
     0     0     1
     2     0     0
     1     1     0

This last step should be done using SPARSE instead of ACCUMARRAY in cases where there is a large number of offices. The trick with BSXFUN could also benefit from converting one of its inputs to sparse, e.g. the vector of offices.

PS: if it was for a homework, I doubt that either method is what was required.

[Keith]

Cedric,

I spent greater part of day chasing the latter method down, but don't think I would have gotten there. Your trick is just what I need. This is not for HW but rather a simulation I am building for my dissertation. I plan to look into bsxfun in greater detail and will get back to you with any questions. Thanks again. Keith

[Cedric]

You're welcome. Keep in mind the last comment about using sparse matrices, especially in your context where the matrix of connections counts is likely to be mostly filled with zero's. If it gets too large (see note *), dealing with sparse matrices will avoid storing any zero, and keep the approach memory-friendly.

(*) You are dealing with class double elements, stored on eight bytes, so it's enough to have around 11,000 offices (=> output is 11,000x11,000) for the output matrix to require almost 1GB RAM to be stored, without talking about intermediary arrays which are created when you perform operations with/on this matrix. Dealing with sparse matrices generates a computation overhead usually though, so don't use them is the size of you problem is reasonable).

[Keith]

Cedric,

I posted one more for you if you have the time: new problem

I am trying it on my own, perhaps I can beat you to it...

Keith

[Cedric]

Seems that Andrei beat us ;-)