在 MATLAB 中矢量化/优化此代码？

Question

我正在构建我的第一个大型 MATLAB 程序，并且我已经设法为所有内容编写原始矢量化代码，直到我开始尝试创建一个表示立体投影中矢量密度的图像。在几次尝试失败后，我访问了 Mathworks 文件交换站点，发现了一个符合我需要的开源程序，由 Malcolm Mclean 提供。使用测试矩阵，他的函数会产生如下内容：

虽然这几乎正是我想要的，但他的代码依赖于三重嵌套的 for 循环。在我的工作站上，在这部分代码中，大小为 25000x2 的测试数据矩阵花费了 65 秒。这是不可接受的，因为我将在我的项目中扩展到大小为 500000x2 的数据矩阵。

到目前为止，我已经能够矢量化最里面的循环（这是最长/最差的循环），但如果可能的话，我想继续并完全摆脱这些循环。这是我需要向量化的 Malcolm 的原始代码：

dmap = zeros(height, width); % height, width: scalar with default value = 32
for ii = 0: height - 1          % 32 iterations of this loop
    yi = limits(3) + ii * deltay + deltay/2; % limits(3) & deltay: scalars
    for jj = 0 : width - 1      % 32 iterations of this loop
        xi = limits(1) + jj * deltax + deltax/2; % limits(1) & deltax: scalars
        dd = 0;
        for kk = 1: length(x)   % up to 500,000 iterations in this loop
            dist2 = (x(kk) - xi)^2 + (y(kk) - yi)^2;
            dd = dd + 1 / ( dist2 + fudge); % fudge is a scalar
        end
        dmap(ii+1,jj+1) = dd;
    end
end

这就是我已经对最内层循环所做的更改（这是对效率的最大消耗）。这将我机器上相同测试矩阵的时间从 65 秒减少到 12 秒，这比我想要的要好，但仍然慢得多。

     dmap = zeros(height, width);
    for ii = 0: height - 1
        yi = limits(3) + ii * deltay + deltay/2;
        for jj = 0 : width - 1
            xi = limits(1) + jj * deltax + deltax/2;
            dist2 = (x - xi) .^ 2 + (y - yi) .^ 2;
            dmap(ii + 1, jj + 1) = sum(1 ./ (dist2 + fudge));
        end
    end

所以我的主要问题是，我可以做任何进一步的更改来优化此代码吗？甚至是解决问题的替代方法？我已经考虑在程序的这一部分使用 C++ 或 F# 而不是 MATLAB，如果我无法使用 MATLAB 代码达到合理的效率水平，我可能会这样做。

还请注意，此时我没有任何其他工具箱，如果有，我知道这将是微不足道的（例如使用统计工具箱中的 hist3）。

Answer 1

内存消耗解决方案

yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width  ) - .5 * deltax;
dx = bsxfun( @minus, x(:), xi ) .^ 2;
dy = bsxfun( @minus, y(:), yi ) .^ 2;
dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
dmap = sum( 1./(dist2 + fudge ) , 3 );

编辑

通过将操作分成块来处理极大的x 和y：

blockSize = 50000; % process up to XX elements at once
dmap = 0;
yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width  ) - .5 * deltax;
bi = 1;
while bi <= numel(x)
    % take a block of x and y
    bx = x( bi:min(end, bi + blockSize - 1) );
    by = y( bi:min(end, bi + blockSize - 1) );
    dx = bsxfun( @minus, bx(:), xi ) .^ 2;
    dy = bsxfun( @minus, by(:), yi ) .^ 2;
    dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
    dmap = dmap + sum( 1./(dist2 + fudge ) , 3 );
    bi = bi + blockSize;
end     

Answer 2

这是一个很好的例子，说明了为什么从 1 开始循环很重要。 ii 和 jj 在 0 处启动的唯一原因是要取消 ii * deltay 和 jj * deltax 术语，但是这会在 dmap 索引中引入顺序性，防止并行化。 p>

现在，通过重写循环，您可以在打开matlabpool 后使用parfor()：

dmap = zeros(height, width);
yi   = limits(3) + deltay*(1:height) - .5*deltay;
matlabpool 8
parfor ii = 1: height
    for jj = 1: width
        xi    = limits(1) + (jj-1) * deltax + deltax/2;
        dist2 = (x - xi) .^ 2 + (y - yi(ii)) .^ 2;
        dmap(ii, jj) = sum(1 ./ (dist2 + fudge));
    end
end
matlabpool close

请记住，打开和关闭池会产生很大的开销（在我的 Intel Core Duo T9300、vista 32 Matlab 2013a 上需要 10 秒）。

PS.我不确定内部循环而不是外部循环是否可以有意义地并行化。您可以尝试将 parfor 切换到内部并比较速度（我建议立即使用大矩阵，因为您已经在 12 秒内运行并且开销几乎一样大）。

Answer 3

或者，这个问题可以通过使用kernel density estimation techniques 来解决。这是统计工具箱的一部分，或者有this KDE implementation by Zdravko Botev（不需要工具箱）。

对于下面的示例代码，对于 N = 500000，我得到 0.3 秒，对于 N = 1000000，我得到 0.7 秒。

N = 500000;
data = [randn(N,2); rand(N,1)+3.5, randn(N,1);];  % 2 overlaid distrib
tic; [bandwidth,density,X,Y] = kde2d(data); toc;
imagesc(density);