带有向量箱的 Matlab histc答案

【问题标题】：Matlab histc with vector bins带有向量箱的 Matlab histc
【发布时间】：2011-11-30 11:53:18
【问题描述】：

如果我有两个大小为 [m x n] 和 [p x n] 的矩阵 A 和 B，我想求出 B 的每一行在 A 中出现的次数的计数，例如：

>> A = rand(5,3)

A =

    0.1419    0.6557    0.7577
    0.4218    0.0357    0.7431
    0.9157    0.8491    0.3922
    0.7922    0.9340    0.6555
    0.9595    0.6787    0.1712

>> B = [A(2,:); A(1,:); A(2,:); A(3,:); A(3,:); A(4,:); A(5,:)] 

B =

    0.4218    0.0357    0.7431
    0.1419    0.6557    0.7577
    0.4218    0.0357    0.7431
    0.9157    0.8491    0.3922
    0.9157    0.8491    0.3922
    0.7922    0.9340    0.6555
    0.9595    0.6787    0.1712

在这种情况下的答案是

ans =

     1     2     2     1     1

虽然和这个例子不同，但总的来说 m >> p

如果 A 和 B 是向量，matlab 的 histc 可以完成这项工作，但如果 bin 是向量，则似乎没有等效项。

目前我是这样做的：

for i=1:length(B)
    indices(i) = find(abs(A/B(i,:)-1) < 1e-15); 
    % division requires a tolerance due to numerical issues
end
histc(indices, 1:size(A,1))

ans =

     1     2     2     1     1

但是由于我有很多这样的矩阵 B，并且 A 和 B 都很大，所以这非常慢。有什么想法可以改进吗？

编辑：

查看到目前为止的方法，我有以下数据：

A                    7871139x3                188907336  double                       
B                        902x3                    21648  double

为了让事情变得更快，我将使用前 10 行 B

B = B(1:10,:);

请注意，对于完整的应用程序，我（目前）拥有 >10^4 个此类矩阵（最终将 >10^6 ....）

我的第一种方法：

tic, C = get_vector_index(A,B); toc
Elapsed time is 36.630107 seconds.

bdecaf 的方法（可以通过删除 if 语句并使用 L1 距离而不是 L2 距离缩短到 ~25 秒）

>> tic, C1 = get_vector_index(A,B); toc
Elapsed time is 28.957243 seconds.
>> isequal(C, C1) 

ans =

     1

oli 的 pdist2 方法

>> tic, C2 = get_vector_index(A,B); toc
Elapsed time is 7.244965 seconds.

>> isequal(C,C2)

ans =

     1

oli的归一化方法

>> tic, C3 = get_vector_index(A,B); toc
Elapsed time is 3.756682 seconds.

>> isequal(C,C3)

ans =

     1

最后我想出了另一种方法，我先搜索第一列，然后在第一列的命中中搜索第二列，递归直到列用完。这是迄今为止最快的....

N = size(A,2);
loc = zeros(size(B,1),1);
for i=1:size(B,1)
    idx{1} = find(A(:,1)==B(i,1));
    for k=2:N, 
        idx{k} = idx{k-1}(find(A(idx{k-1},k)==B(i,k))); 
    end
    loc(i) = idx{end};
end
C = histc(loc, 1:size(A,1));

导致：

>> tic, C4 = get_vector_index(A,B); toc
Elapsed time is 1.314370 seconds.

>> isequal(C, C4)

ans =

     1

另请注意，使用intersect 会慢得多：

>> tic, [~,IA] = intersect(A,B,'rows'); C5 = histc(IA,1:size(A,1)); toc
Elapsed time is 44.392593 seconds.

>> isequal(C,C5)

ans = 

    1

【问题讨论】：

顺便说一下，我意识到我当前方法中的矩阵除法也是错误的——它（几乎）有可能得到一个虚假的误报，但它会产生一个错误，因为只有在我拥有的所有情况下，A 中 B 的每个向量都有一个匹配项。
只是一个技术性问题 - 您输入的代码只有在 B 中的所有行在 A 中都有一行时才有效。此外，如果 A 中有两个几乎相同的行，则会出现问题。该行为也与 histc 不同。 histc 分类到 最近的 垃圾箱 - 当它正好在垃圾箱中时，你就成功了。
还要小心使用长度——它总是返回最大的尺寸。所以如果p<n会有麻烦。
抱歉，长度(B) 应该是 size(B,1) ......当我创建示例时，它悄悄出现了。
对于前一个，是的，我的描述可能缺少，因为在我的情况下，B 中的所有行在 A 中都有一行 - 只是 histc 将解决我在一维情况下的问题

标签： matlab histogram vectorization

【解决方案1】：

也许您可以将它们标准化，以便检查它们的点积是否为1

A = rand(5,3);
B = [A(2,:); A(1,:); A(2,:); A(3,:); A(3,:); A(4,:); A(5,:)];
A2=bsxfun(@times,A,1./sqrt(sum(A.^2,2))); %%% normalize A
B2=bsxfun(@times,B,1./sqrt(sum(B.^2,2))) %%% normalize B
sum(A2*B2'>1-10e-9,2) %%% check that the dotproduct is close to 1

ans =

     1
     2
     2
     1
     1

如果您需要更快但近似的东西，我建议您使用 flann 库，该库用于快速近似最近邻：

http://www.cs.ubc.ca/~mariusm/index.php/FLANN/FLANN

【讨论】：

我很喜欢这种方法，但是矩阵乘法不适合内存（例如 A 是 7871139x3，B 是 903x3）
我也看过 FLANN，但我不认为我可以安装它，因为我没有 sudoers 帐户 - 从源代码编译需要 cmake（未安装，并且安装它对我来说简直就是炸弹，大概是因为我不是 sudoer），目前他们的预建库仅是 ubuntu/debian（服务器是 Red Hat）
如果你有 linux 64 位，我可以把我的编译版本发给你。

【解决方案2】：

其实，更简单的做法是：

sum(10e-9>pdist2(A',B'),2)

它计算所有成对距离、阈值和计数。

【讨论】：

【解决方案3】：

我会这样解决：

indices = zeros(size(A,1),1);
for i=1:size(B,1)
    distances = sum( ( repmat(B(i,:),size(A,1),1)-A ).^2 ,2);
    [md,im]=min(distances);
    if md < 1e-9
      indices(im) = indices(im)+1;
    end
end

如果你删除 if 它只会排序到最近的垃圾箱。

【讨论】：

比@oli 的方法稍慢

【解决方案4】：

我实际上将此解决方案作为问题的编辑，但为了接受答案，我也将解决方案放在这里：

N = size(A,2);
loc = zeros(size(B,1),1);
for i=1:size(B,1)
    idx{1} = find(A(:,1)==B(i,1));
    for k=2:N, 
        idx{k} = idx{k-1}(find(A(idx{k-1},k)==B(i,k))); 
    end
    loc(i) = idx{end};
end
C = histc(loc, 1:size(A,1));

导致：

>> tic, C4 = get_vector_index(A,B); toc
Elapsed time is 1.314370 seconds.

>> isequal(C, C4)

ans =

     1

【讨论】：