algorithm - 最小化到加权网格的距离

Question

假设您有一个 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 左右的操作中或附近执行此操作？

Answer 1

理论

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

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

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

令 W 表示嵌入大小为 3n x 3m 的零网格中的原始权重网格(居中)。

然后是过滤后的(cross-correlation)结果

 R = F o W

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

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

我应该注意，虽然我在上面明确地使用了互相关，但你会得到与 convolution 相同的结果。只是因为你的过滤器 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's FFT library您可以在 CPU 上花费的时间的一小部分内计算 2D FFT。关键是，FFT 方法会随着您投入更多硬件来进行计算而扩展，而蛮力方法的扩展性会更差。

观察

在执行此操作时，我观察到几乎每次 best_x,bext_y 值都是 floor(n/2)+-1 之一。这意味着距离项很可能主导整个计算，因此，您可以只计算 4 个值的 total_dist 的值，从而使该算法变得微不足道!