MATLAB矢量化:计算邻域矩阵

MGA*_*MGA 3 performance matlab vector matrix vectorization

给定两个向量XY长度n,表示平面上的点和邻域半径rad,是否有一种矢量化方法来计算点的邻域矩阵?

换句话说,可以对以下(对于大的痛苦缓慢n)循环进行矢量化:

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
Run Code Online (Sandbox Code Playgroud)

Div*_*kar 5

方法#1

bsxfun 基于方法 -

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

方法#2

Distance matrix calculation using matrix-multiplication 基于方法(可能更快) -

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
Run Code Online (Sandbox Code Playgroud)

方法#3

随着pdistsquareform-

A = [X(:) Y(:)]
out = squareform(pdist(A))<radius
out(1:n+1:end)= 0
Run Code Online (Sandbox Code Playgroud)

方法#4

你可以pdist像以前的方法一样使用,但是避免squareform使用一些逻辑索引来获得邻域矩阵的最终输出,如下所示 -

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)
Run Code Online (Sandbox Code Playgroud)

注意:可以看出,对于上述所有方法,我们使用预分配输出.预先分配的快速方法将是out(n,n) = 0并且基于this wonderful blog on undocumented MATLAB.这应该真的加快这些方法!