在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

Mem消耗解决方案

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