MATLAB vectorization: computing a neighborhood matrix

Question

Given two vectors X and Y of length n, representing points on the plane, and a neighborhood radius rad, is there a vectorized way to compute the neighborhood matrix of the points?

In other words, can the following (painfully slow for large n) loop be vectorized:

neighborhood_mat = zeros(n, n);
for i = 1 : n
    for j = 1 : i - 1
        dist = norm([X(j) - X(i), Y(j) - Y(i)]);
        if (dist < radius)
            neighborhood_mat(i, j) = 1;
            neighborhood_mat(j, i) = 1;
        end
    end
end

Just curious - Did any of the approaches listed here work for you? If it did, consider accepting it? — Divakar, Feb 26 '15 at 13:44
@Divakar Yes, all do the job nicely; I went with bsxfun in the end. Sorry for the delay in accepting, I was away for a while. — MGA, Feb 26 '15 at 17:50
Ah that's okay! Glad to see `bsxfun` make it through.. again! — Divakar, Feb 26 '15 at 17:53

score 5 · Accepted Answer · edited May 23 '17 at 12:10

Approach #1

bsxfun based approach -

out = bsxfun(@minus,X,X').^2 + bsxfun(@minus,Y,Y').^2 < radius^2
out(1:n+1:end)= 0

Approach #2

Distance matrix calculation using matrix-multiplication based approach (possibly faster) -

A = [X(:) Y(:)]
A_t = A.';  %//'
out = [-2*A A.^2 ones(n,3)]*[A_t ; ones(3,n) ; A_t.^2] < radius^2
out(1:n+1:end)= 0

Approach #3

With pdist and squareform -

A = [X(:) Y(:)]
out = squareform(pdist(A))<radius
out(1:n+1:end)= 0

Approach #4

You can use pdist as with the previous approach, but avoid squareform with some logical indexing to get the final output of neighbourhood matrix as shown below -

A = [X(:) Y(:)]
dists = pdist(A)< radius

mask_lower = bsxfun(@gt,[1:n]',1:n)  %//'
%// OR tril(true(n),-1)

mask_upper = bsxfun(@lt,[1:n]',1:n)  %//'
%// OR mask_upper = triu(true(n),1)
%// OR mask_upper = ~mask_lower; mask_upper(1:n+1:end) = false;

out = zeros(n)
out(mask_lower) = dists

out_t = out'  %//'
out(mask_upper) = out_t(mask_upper)

Note: As one can see, for the all above mentioned approaches, we are using pre-allocation for the output. A fast way to pre-allocate would be with out(n,n) = 0 and is based upon this wonderful blog on undocumented MATLAB. This should really speed up those approaches!

I bet you have these types of answers pre-prepared somewhere... Copy and paste... You are too quick. :-) — kkuilla, Feb 25 '15 at 11:14
@kkuilla Well, I should bookmark that bsxfun link somewhere! :) And yeah I guess I like these type of questions ;) — Divakar, Feb 25 '15 at 11:15

score 1 · Answer 2 · edited May 23 '17 at 10:24

The following approach is great if the number of points in your neighborhoods is small or you run low on memory using the brute-force approach:

If you have the statistics toolbox installed, you can have a look at the rangesearch method. (Free alternatives include the k-d tree implementations of a range search on the File Exchange.)

The usage of rangesearch is straightforward:

P = [X,Y];
[idx,D] = rangesearch(P, P, rad);

It returns a cell-array idx of the indices of nodes within reach and their distances D.

Depending on the size of your data, this could be beneficial in terms of speed and memory. Instead of computing all pairwise distances and then filtering out those that are large, this algorithm builds a data structure called a k-d tree to more efficiently search close points.

You can then use this to build a sparse matrix:

I = cell2mat(idx.').';
J = runLengthDecode(cellfun(@numel,idx));
n = size(P,1);
S = sparse(I,J,1,n,n)-speye(n);

(This uses the runLengthDecode function from this answer.)

You can also have a look at the KDTreeSearcher class if your data points don't change and you want to query your data lots of times.

MATLAB vectorization: computing a neighborhood matrix

2 Answers2