Can I speed up an ismember operation?

I have a list of 100 random numbers 1-10:

 numpoints = 100;
 gridnum = 10;
 grid_ind = randi(gridnum,[numpoints 1]); % 100 points in 10 random cells

and each of the 10 possible values has a corresponding list (of random length) of values, which are indexes to grid_ind:

 for i=1:gridnum
     neighbgrids{i} = randi(numpoints,[randi(5) 1]);
 end

So far, so good. Now, I need to find the indexes for all elements in grid_ind which have the same values as in each cell list in neighbgrids. So I use:

 for i=1:numpoints
     neighblist{i} = find(ismember(grid_ind,neighbgrids{grid_ind(i)}));
 end

But this turns out to be very slow. Is there a faster way to do this part?

[Matt J]

The ismember operation looks very strange. The expression neighbgrids{grid_ind(i)}) will always take values in 1...numpoints, but grid_ind will always take values in the much smaller set 1...gridnum. In other words, you have lots of numbers in neighbgrids that grid_ind simply can never have membership in. Are you sure the code is correct?

[Christopher]

Hi Matt. Yes I think that is right and it is correct. This is a simplification of another code (complicated codes seem less likely to get answers, so I simplified) which is the following:

 % model dimensions
 xnum = 10;
 ynum = 10;
 xsize = 1;
 ysize = 1;

 numpoints = 100; % number of particles in domain

 y = rand(numpoints,1)*ysize; % random positions
 x = rand(numpoints,1)*xsize; % random positions

 % find nearest element 
 ygridpos = ceil(ypos.*ynum/ysize); % bottom node
 xgridpos = ceil(xpos.*xnum/xsize); % right node

 % get grid indices for each point
 grid_ind = sub2ind([ynum xnum], ygridpos, xgridpos);

 % get list of points which have an index in the same and neighboring elements
 for i=1:numpoints
    neighblist{i} = find(ismember(grid_ind,neighbgrids{grid_ind(i)}));
 end

The part that is more complicated is finding the indices for the neighboring elements (the cell matrices 'neighbgrids'), which I contruct as follows:

 refs    = vec2mat(1:1:ynum*xnum,xnum)';
 refs    = padarray(refs,[1 1], 'circular');
 refs_tl = refs(1:end-2,1:end-2); refs_tl(:,1)=0; refs_tl(1,:)=0;
 refs_t  = refs(1:end-2,2:end-1); refs_t(1,:)=0;
 refs_tr = refs(1:end-2,3:end);   refs_tr(:,end)=0; refs_tr(1,:)=0;
 refs_l  = refs(2:end-1,1:end-2); refs_l(:,1)=0;
 refs_m  = refs(2:end-1,2:end-1);
 refs_r  = refs(2:end-1,3:end);   refs_r(:,end)=0;
 refs_bl = refs(3:end,1:end-2);   refs_bl(:,1)=0; refs_bl(end,:)=0;
 refs_b  = refs(3:end,2:end-1);   refs_b(end,:)=0;
 refs_br = refs(3:end,3:end);     refs_br(:,end)=0; refs_br(end,:)=0;

 for i=1:ynum
    for j=1:xnum
        node = refs(i+1,j+1);
        if node>0
            thisgrid = [refs_tl(i,j); refs_t(i,j); refs_tr(i,j); refs_l(i,j); refs_m(i,j); refs_r(i,j); refs_bl(i,j); refs_b(i,j); refs_br(i,j)];
            thisgrid(thisgrid==0)=[];
            neighbgrids{node} = thisgrid;
        end
    end
 end

The above is just a preallocation so it doesn't need to be efficient. The construction of neighblist cells, however, is repeated millions of times. Since the positions of the points are continuously changed, it needs to be efficient.

[per isakson]

"% 100 points in 10 random cells" &nbsp Is this comment misplaced? &nbsp grid_ind is a [100x1 double]

[Matt J]

@Christopher,

Your code takes .006 sec to run on my machine. What speed were you hoping for?

If these are just example numbers, what are the actual typical values of numpoints and gridnum? Is numpoints always much larger than gridnum?

[Christopher]

gridnum is ideally going to be around 100-1000 and numpoints will always be around 2-25 times that number.

Answer 1

[Matt J]

This should speed things up.

    for i=1:gridnum
     neighblist{i} = find(ismember(grid_ind,neighbgrids{i}));
    end

    neighblist=neighblist(grid_ind);

I see a 10x speed-up over your originally posted example.

[Christopher]

Thanks! This gives me about a 5x speedup the way it is implemented in my code.

Hopefully I can think of a way to avoid ismember all together, in the future, as I was hoping for significantly more speedup.

[Matt J]

It seems inefficient to be using find(). The output of ismember is enough to index into grid_inds and then you could maintain neighblist as a matrix instead of a cell, which is more efficient:

neighblist=zeros(gridnum,numpoints);

   for i=1:gridnum
     neighblist(i,:) = ismember(grid_ind,neighbgrids{i});
   end