1

我想加快我的 Matlab 代码。通常，我会找到避免使用 for 循环来获得计算时间的方法，但在这种情况下我遇到了障碍。我必须计算点网格中的值，但计算值需要逻辑运算和向量求和，这使实现复杂化。这段代码在我的机器上运行大约 8 秒：

clear all
% Grid
dLimsX=[-100 +100];
dLimsY=[-100 +100];
dStep=1;
[x_map, y_map]=meshgrid((dLimsX(1):dStep:dLimsX(2)),(dLimsY(1):dStep:dLimsY(2)));
nPoints_map=numel(x_map);
% Inputs
smallDistance=1e-3;
N=10e3;
scaleFactor=10;
x_input = sin(linspace(0,1,N));
y_input = cos(linspace(0,1,N));
z_input = linspace(0,1,N);
tic
A=zeros(size(x_map));
for r=1:size(x_map,1)
    y0=y_map(r,1);
    for c=1:size(x_map,2)
        x0=x_map(1,c);
        
        idxTemp = find((x0-x_input).^2+(y0-y_input).^2>smallDistance); % do not consider in the calculation the inputs too close to the point
        A(r,c) = sum( scaleFactor * z_input(idxTemp) .* (y0-y_input(idxTemp)) ./ ((x0-x_input(idxTemp)).^2 +(y0-y_input(idxTemp)).^2+eps) );
    end
end
toc
4

2 回答 2

2

加速代码并不是要删除 for 循环。我似乎有很多矢量化代码比循环等效代码慢的情况。在过去的 20 年中，MATLAB 的循环越来越快，它们不再是导致减速的重要原因。例如，以下仅比sum(x)对 100 万个元素求和慢 4 倍：

y = 0;
for ii = 1:numel(x)
   y = y+x(ii);
end

如果循环内的计算成本更高，则循环开销完全消失。

您仍然可以从矢量化中受益的原因是，在循环代码中，通常会提取矩阵的一行或一列。这涉及复制数据，这是昂贵的。另一方面，如果向量化代码需要一个大的中间矩阵，那么将该矩阵存储在内存中将成为使向量化代码显着变慢的瓶颈。内存访问通常是问题所在。

为了使您的代码更快，您应该首先关注避免重复计算。例如，(y0-y_input).^2是计算3*size(x_map,2)次数！（数据子集的 1/3 时间，但索引删除的点数很少）。

此外，您应该使用逻辑索引，并避免使用find. A(find(condition))与 相同A(condition)，但速度较慢。

您的循环在我的机器上运行约 10.5 秒，此版本运行约 5.1 秒：

tic
A = zeros(size(x_map));
for r = 1:size(x_map,1)
    y0 = y_map(r,1);
    dy2 = (y0-y_input).^2;
    for c = 1:size(x_map,2)
        x0 = x_map(1,c);
        dx2 = (x0-x_input).^2;
        idxTemp = dx2 + dy2 > smallDistance; % do not consider in the calculation the inputs too close to the point
        A(r,c) = sum(scaleFactor * z_input(idxTemp) .* (y0-y_input(idxTemp)) ./ (dx2(idxTemp) + dy2(idxTemp) + eps));
    end
end
toc

可能会有进一步的改进，例如避免y0-y_input在内循环中的重复计算。

于 2020-08-07T15:12:09.797 回答
0

Cris Luengo 的回答给了我一个很大的提示，让我思考哪些计算可以避免重复。正如 Cris 提出的那样，避免重新计算x0-x_input和已经将计算时间减少了 50-60%。y0-y_input

此外，当在另一个循环中使用代码时，它帮助我区分了哪些变化和哪些只能计算一次。就我而言，x0, x_input, y0,y_input始终保持不变；，z_input随着每次新的迭代而变化。所以我把计算分成两部分：首先，计算一次所有不变的东西；第二，计算所需值。

为了能够通过简单的矩阵乘法执行第二次计算，我将 xy 值从矩阵（网格）重新排列为向量。这就是我所做的：

% Arrange x-y values in two vectors
x = x_map(:);
y = y_map(:);
nPoints = numel(x);

z_input = linspace(0,1,N)';

tic
a = zeros(nPoints,N);
for p = 1:nPoints
    x0 = x(p);
    y0 = y(p);
    dx2 = (x0-x_input).^2;
    dy2 = (y0-y_input).^2;
    idxTemp = dx2 + dy2 > smallDistance; % do not consider in the calculation the inputs too close to the point
    a(p,:) = (scaleFactor * 1 .* dy2 ./ (dx2 + dy2 + eps)) .*idxTemp;
end
toc
tic
A2 = a*z_input;
toc

% Check that the values calculated with the alternative method are correct
mean(mean(abs((A(:)-A2)./A(:))))

对结果的一些评论：

  • 原始代码：在我的家用 PC 上约 122 秒（比我原来的办公室 PC 慢得多）

  • Cris 提出的代码：~50 s

  • 上面的代码：第一部分约 75 秒。第二部分约 0.6 秒。

  • 替代计算的相对误差约为 1E-17。

因此，在多次重复计算的情况下，预先计算不改变的内容是值得的。缺点是存储变量所需的内存使用量较大a。

我将不胜感激您的反馈。

于 2020-08-08T04:09:17.500 回答