MATLAB 矢量化:计算邻域矩阵
MATLAB vectorization: computing a neighborhood matrix
给定两个向量 X
和 Y
,长度 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
方法 #1
基于bsxfun
的方法-
out = bsxfun(@minus,X,X').^2 + bsxfun(@minus,Y,Y').^2 < radius^2
out(1:n+1:end)= 0
方法 #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
方法 #3
与pdist
and squareform
-
A = [X(:) Y(:)]
out = squareform(pdist(A))<radius
out(1:n+1:end)= 0
方法 #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)
注意: 可以看出,对于上述所有方法,我们都对输出使用预分配。一种快速的预分配方法是使用 out(n,n) = 0
并且基于 this wonderful blog on undocumented MATLAB
。这应该真的加快这些方法!
如果您的社区中的点数很少,或者您 运行 使用蛮力方法的内存不足,则以下方法非常有用:
如果你安装了统计工具箱,可以看看rangesearch
方法。
(免费替代品包括文件交换中的 k-d tree implementations of a range search。)
rangesearch
的用法很简单:
P = [X,Y];
[idx,D] = rangesearch(P, P, rad);
它 returns 一个元胞数组 idx
范围内节点的索引及其距离 D
。
根据您的数据大小,这可能会在速度和内存方面有所帮助。
该算法不是计算所有成对距离然后过滤掉那些较大的距离,而是构建一个称为 k-d tree 的数据结构来更有效地搜索接近点。
然后您可以使用它来构建一个 sparse
矩阵:
I = cell2mat(idx.').';
J = runLengthDecode(cellfun(@numel,idx));
n = size(P,1);
S = sparse(I,J,1,n,n)-speye(n);
(这使用了 中的 runLengthDecode
函数。)
如果你的数据点没有变化,而且你想多次查询你的数据,你也可以看看 KDTreeSearcher
class。
给定两个向量 X
和 Y
,长度 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
方法 #1
基于bsxfun
的方法-
out = bsxfun(@minus,X,X').^2 + bsxfun(@minus,Y,Y').^2 < radius^2
out(1:n+1:end)= 0
方法 #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
方法 #3
与pdist
and squareform
-
A = [X(:) Y(:)]
out = squareform(pdist(A))<radius
out(1:n+1:end)= 0
方法 #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)
注意: 可以看出,对于上述所有方法,我们都对输出使用预分配。一种快速的预分配方法是使用 out(n,n) = 0
并且基于 this wonderful blog on undocumented MATLAB
。这应该真的加快这些方法!
如果您的社区中的点数很少,或者您 运行 使用蛮力方法的内存不足,则以下方法非常有用:
如果你安装了统计工具箱,可以看看rangesearch
方法。
(免费替代品包括文件交换中的 k-d tree implementations of a range search。)
rangesearch
的用法很简单:
P = [X,Y];
[idx,D] = rangesearch(P, P, rad);
它 returns 一个元胞数组 idx
范围内节点的索引及其距离 D
。
根据您的数据大小,这可能会在速度和内存方面有所帮助。 该算法不是计算所有成对距离然后过滤掉那些较大的距离,而是构建一个称为 k-d tree 的数据结构来更有效地搜索接近点。
然后您可以使用它来构建一个 sparse
矩阵:
I = cell2mat(idx.').';
J = runLengthDecode(cellfun(@numel,idx));
n = size(P,1);
S = sparse(I,J,1,n,n)-speye(n);
(这使用了 runLengthDecode
函数。)
如果你的数据点没有变化,而且你想多次查询你的数据,你也可以看看 KDTreeSearcher
class。