5

假设您有一个 1000x1000 的正整数权重 W 网格。

我们想要找到最小化平均加权距离的单元格。到每个单元格。

这样做的蛮力方法是遍历每个候选单元格并计算距离：

int best_x, best_y, best_dist;

for x0 = 1:1000,
    for y0 = 1:1000,

        int total_dist = 0;

        for x1 = 1:1000,
            for y1 = 1:1000,
                total_dist += W[x1,y1] * sqrt((x0-x1)^2 + (y0-y1)^2);

        if (total_dist < best_dist)
            best_x = x0;
            best_y = y0;
            best_dist = total_dist;

这需要大约 10^12 次操作，这太长了。

有没有办法在~10^8 左右的操作中或附近做到这一点？

4

1 回答 1

5

理论

这可以在 O(nm log nm ) 时间内使用过滤器，其中 n，m 是网格尺寸。

您需要定义一个 size 的过滤器2n + 1 x 2m + 1，并且您需要（居中）将原始权重网格嵌入到 size 的零网格中3n x 3m。过滤器需要是距原点的距离权重(n,m)：

 F(i,j) = sqrt((n-i)^2 + (m-j)^2)

让W表示嵌入在大小为 0 的网格中的原始权重网格（居中）3n x 3m。

然后过滤（互相关）结果

 R = F o W

将为您提供total_dist网格，只需使用min R（忽略您放入 W 中的额外嵌入零）即可找到最佳x0, y0位置。

图像（即网格）过滤非常标准，可以在各种不同的现有软件（如 matlab）中使用imfilter命令完成。

我应该注意，虽然我在上面明确使用了互相关，但你会得到与卷积相同的结果，只是因为你的滤波器F是对称的。一般来说，图像过滤器是互相关的，而不是卷积的，尽管这两种操作非常相似。

O(nm log nm) 运行时间的原因是因为可以使用 2D FFT 完成图像过滤。

实施

这是 Matlab 中的两种实现，两种方法的最终结果是相同的，并且它们以非常简单的方式进行基准测试：

m=100;
n=100;
W0=abs(randn(m,n))+.001;

tic;

%The following padding is not necessary in the matlab code because
%matlab implements it in the imfilter function, from the imfilter
%documentation:
%  - Boundary options
% 
%        X            Input array values outside the bounds of the array
%                     are implicitly assumed to have the value X.  When no
%                     boundary option is specified, imfilter uses X = 0.

%W=padarray(W0,[m n]);

W=W0;
F=zeros(2*m+1,2*n+1);

for i=1:size(F,1)
    for j=1:size(F,2)
        %This is matlab where indices start from 1, hence the need
        %for m-1 and n-1 in the equations
        F(i,j)=sqrt((i-m-1)^2 + (j-n-1)^2);
    end
end
R=imfilter(W,F);
[mr mc] = ind2sub(size(R),find(R == min(R(:))));
[mr, mc]
toc;

tic;
T=zeros([m n]);
best_x=-1;
best_y=-1;
best_val=inf;
for y0=1:m
    for x0=1:n

        total_dist = 0;

        for y1=1:m
            for x1=1:n
                total_dist = total_dist + W0(y1,x1) * sqrt((x0-x1)^2 + (y0-y1)^2);
            end
        end

        T(y0,x0) = total_dist;
        if ( total_dist < best_val ) 
            best_x = x0;
            best_y = y0;
            best_val = total_dist;
        end

    end
end
[best_y best_x]
toc;

diff=abs(T-R);
max_diff=max(diff(:));
fprintf('The max difference between the two computations: %g\n', max_diff);

表现

对于 800x800 网格，在我的 PC 上肯定不是最快的，FFT 方法在 700 多秒内进行评估。蛮力方法在几个小时后没有完成，我必须杀死它。

就进一步的性能提升而言，您可以通过迁移到 GPU 等硬件平台来实现它们。例如，使用CUDA 的 FFT 库，您可以在 CPU 上花费的时间的一小部分内计算 2D FFT。关键是，FFT 方法会随着您投入更多硬件来进行计算而扩展，而蛮力方法的扩展性会更差。

观察

在实现这一点时，我观察到几乎每次，这些best_x,bext_y值都是floor(n/2)+-1. 这意味着距离项很可能在整个计算中占主导地位，因此，您可以total_dist只计算 4 个值的值，从而使该算法变得微不足道！

于 2012-06-29T21:11:18.280 回答