How to speed up 2 loops with the intersect function?

Hi

I have some code, which is very slow, because of the 2 loops. I was wondering if there was a clever way to speed this up? For example is it possible to somehow vectorize the following code? Thanks for your help!

for i=1:size(Lien,1)
   k=1;
   for j=1:size(Lien,1)
      sommetscommuns=intersect(Lien(i,:),Lien(j,:));
      if size(sommetscommuns,2)==2
         voisins(i,k)=j;
         k=k+1;
      end
   end
   if voisins(i,1)~=0 && voisins(i,2)~=0 && voisins(i,3)~=0
      Qnor(i,1) = (dot(Normales(i,:),Normales(voisins(i,1),:)) + ... 
                   dot(Normales(i,:),Normales(voisins(i,2),:)) + ... 
                   dot(Normales(i,:),Normales(voisins(i,3),:)))/3;
   elseif voisins(i,1)~=0 && voisins(i,2)~=0 && voisins(i,3)==0
      Qnor(i,1) = (dot(Normales(i,:),Normales(voisins(i,1),:)) + ...
                   dot(Normales(i,:),Normales(voisins(i,2),:)))/2;
   elseif voisins(i,1)~=0 && voisins(i,2)==0 && voisins(i,3)==0
      Qnor(i,1) = dot(Normales(i,:),Normales(voisins(i,1),:));
   elseif voisins(i,1)==0 && voisins(i,2)==0 && voisins(i,3)==0
      Qnor(i,1)=1;
   end
end

the 'Lien' matrix is around 6300 rows...

[Guillaume]

What is in Lien? Integers within some range (what range?) or something more complex?

Also, in the first part of your loop, you're looking for rows that have exactly 2 elements in common with the current row. Do you ever expect more than 2?

From the second part of the loop, it looks like there's never more than 3 rows that have two elements in common with any row. Is that true?

You say your Lien matrix has ~6300 rows, but how many columns?

[nicolas]

HI,

Lien is a matrix "double" with 3 columns. the doubles are positive from 1 to...infinity.

In the first part of the loop, I'm looking for rows that have exactly 2 elements in common with the current row. I only want that. There is no row with 3 elements in common. Rows with 1 or 2 elements in common exist.

Yes, from the second part of the loop, never more than 3 rows have two elements in common with any row.

Answer 1

[Roger Stafford]

Some possible improvements:

1) You could compute 'Qnor' in the following way:

   if k<=1
     Qnor(i,1) = 1;
   else
     Qnor(i,1) = sum(Normales(i,:).*(Normales(voisins(i,1),:)+...
                 (k>=3)*Normales(voisins(i,2),:)+...
                 (k>=4)*Normales(voisins(i,3),:)))/min(k-1,3);
   end

2) For every pair (i,j) you will get exactly the same intersection values for 'sommetscommuns' as for (j,i). Consequently the inner loop can be

   "for j=i:size(Lien,1)" instead of "for j=1:size(Lien,1)"

followed by both voisins(i,k)/Qnor(i,:) computations and voisins(j,k)/Qnor(j,:) computations. This cuts down the number of calls on 'intersect' by about half.

3) Make sure both 'voisins' and 'Qnor' arrays are initially preallocated using the 'zeros' function to the necessary sizes so that assignments within the loops will not result in new memory being allocated. That is very important!

[nicolas]

Hi,

Thanks you for your answer. 1)It doesn't work because if voisins(i,1)=0 then Normales(0,:) doesn't exist. 2)It doesn't work; it's stop after some loops...I don't know why? 3)Yes, both 'voisins' and 'Qnor' arrays are initially preallocated using the 'zeros' function.

Answer 2

[Jan]

If intersect is the bottleneck, take into account that it sorts the input every time. Then sort the columns of Lien once before the loop.

In older Matlab versions dot(x,y) was much slower than x*y'.

[nicolas]

Hi, Thanks for your answer, but when I sort the Lien matrix before loop, it doesn't work, it stop at yhe first loop... Moreover I need to keep the rows numbers in the good order.