Need to accelerate an operation on a matrix

I'm writing code that involves looping through each element in arrays with appx. 10e5 elements. For each value in the array, a filter is applied to the cells around it. This filtered neighborhood is used to determine the value assigned to a value in a new matrix. This process is the bottleneck in the code, so I'd appreciate any help in making this loop faster. I tried a gpuArray-based implementation but it was significantly slower than what I'm currently doing. One thing I find curious is that a simple indexing of the large array takes longer than some of the more complex operations.

Here's the code:

for k = 1:ncells

    pos = allinds(k,:);
    yf = pos(1) + sensdist;
    xf = pos(2) + sensdist;

    % Slow operation: apply the neighborhood filter
    neighbors = subfield(yf-sensdist:yf+sensdist, xf-sensdist:xf+sensdist) .*obj.nhood;

    %Slow operation: look up next value for this index
    nextstate = trans(1,subfield(yf,xf));
    %Slow operation: check how many neighbors possess are the next value
    quorum = sum(sum(neighbors == nextstate));
    if quorum > 1;
      if ~isempty(zipintersect(quorum,obj.go))
        temp(k) = nextstate;
      end
    end
  end

[Matt J]

How big is sensdist, and hence the neighborhoods? How big is ncells as compared to numel(subfield)?

One thing I find curious is that a simple indexing of the large array takes longer than some of the more complex operations.

There's nothing "simple" about indexing arrays. Getting things in and out of memory is always one of the most troublesome operations to perform, especially when jumping around nonsequentially in the array, as you're doing here.

[Davis]

sensdist is typically 3 or 5 (this code is for generating cellular automata).

[Matt J]

What about ncells? How big is that?

[Davis]

ncells ranges from 10e4 to ~10e5

Answer 1

[Matt J]

This may help some. The key points are,

• Pre-allocate temp if you're not doing so already.
• Create a table of NextStates in advance for all positions (yf,xf). Makes it simpler to access
• Get rid of 2D indexing and replaces it by linear indexing.
• Get rid of double sum() by shaping the neighborhood as a vector all the time instead of a 2D submatrix.
• Get rid of double "if"

    sz=size(subfield);
    w=sensdist+1;
    ww=w+sensdist;

    [ii,jj]=ndgrid(1:ww);
    jumps= sub2ind(sz,ii,jj) - sub2ind(sz, w, w) ;
      jumps=jumps(:);

      obj.nhood=obj.nhood(:);  

    NextStates=reshape( trans(1,subfield(:)), sz ) ;

    allpos=sub2ind(sz,allinds(:,1)+sensdist,allinds(:,2)+sensdist);

    temp=nan(1,ncells);

    for k = 1:ncells

        pos=allpos(k);
        % Slow operation: apply the neighborhood filter
        neighbors = subfield(pos+jumps) .*obj.nhood;
        %Slow operation: look up next value for this index
        nextstate = NextStates(pos);

        %Slow operation: check how many neighbors possess are the next value
        quorum = sum(neighbors == nextstate);
        if quorum > 1 && ~isempty(zipintersect(quorum,obj.go))
            temp(k) = nextstate;
        end

      end

[Matt J]

Since ncells seems to cover almost your entire sufield array, it seems to make sense to pre-tabulate all the weighted neighbor values as well

    sz=size(subfield);
    w=sensdist+1;
    ww=w+sensdist;

    [ii,jj]=ndgrid(1:ww);
    jumps= sub2ind(sz,ii,jj) - sub2ind(sz, w, w) ;
      jumps=jumps(:);
      obj.nhood=obj.nhood(:); 

    NextStates=reshape( trans(1,subfield(:)), sz ) ;

    allpos=sub2ind(sz,allinds(:,1)+sensdist,allinds(:,2)+sensdist);

    NeighborTable=subfield(bsxfun(@plus,allpos,jumps(:).'));
      NeighborTable=bsxfun(@times,NeighborTable, obj.nhood(:).');

    temp=nan(1,ncells);

    for k = 1:ncells

        neighbors = NeighborTable(k,:);

        nextstate = NextStates(pos);

        quorum = sum(neighbors == nextstate);
        if quorum > 1 && ~isempty(zipintersect(quorum,obj.go))
            temp(k) = nextstate;
        end
      end

[Davis]

Your suggestions dropped runtime by a factor of 3. Thanks! I'm going to see how fast I can get things from this point using parfor/gpuarray.