Need to accelerate an operation on a matrix
1 view (last 30 days)
Show older comments
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
4 Comments
Accepted Answer
Matt J
on 28 Jan 2013
Edited: Matt J
on 28 Jan 2013
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
2 Comments
Matt J
on 4 Feb 2013
Edited: Matt J
on 9 Feb 2013
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
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!