MATLAB 向量化:计算邻域矩阵

Posted

技术标签:

【中文标题】MATLAB 向量化:计算邻域矩阵【英文标题】:MATLAB vectorization: computing a neighborhood matrix 【发布时间】:2015-04-27 08:15:26 【问题描述】:

给定两个向量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

【问题讨论】:

只是好奇 - 这里列出的任何方法对您有用吗?如果有,考虑接受吗? @Divakar 是的,一切都做得很好;最后我选择了 bsxfun。很抱歉延迟接受,我离开了一段时间。 啊,没关系!很高兴看到bsxfun 再次成功! 【参考方案1】:

方法#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

使用pdistsquareform -

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。这应该真的加快这些方法!

【讨论】:

我敢打赌,您已经在某处预先准备了这些类型的答案...复制和粘贴...您太快了。 :-) @kkuilla 好吧,我应该在某处为那个 bsxfun 链接添加书签! :) 是的,我想我喜欢这类问题 ;)【参考方案2】:

如果您附近的点数很少或您使用蛮力方法内存不足,则以下方法非常有用:

如果你安装了统计工具箱,你可以看看rangesearch 方法。 (免费替代品包括 File Exchange 上的 k-d tree implementations of a range search。)

rangesearch 的用法很简单:

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

它返回一个单元数组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);

(这使用来自this answer 的runLengthDecode 函数。)


如果您的数据点没有变化并且您想多次查询数据,您也可以查看 KDTreeSearcher 类。

【讨论】:

以上是关于MATLAB 向量化:计算邻域矩阵的主要内容,如果未能解决你的问题,请参考以下文章

MATLAB中向量数组的向量范数

如何使用向量化代码从 MATLAB 中的两个向量生成所有对?

Matlab - 在邻域中应用函数

向量化 matlab 列归一化

matlab for循环 改写成 矩阵算法

电力负荷预测基于matlab粒子群算法优化支持向量机预测电力负荷含Matlab源码 1225期