我想加快我的Matlab代码。通常，我会找到避免循环以获得计算时间的方法，但在这种情况下，我遇到了困难。我必须计算点网格中的值，但计算值需要逻辑运算和对向量的求和，这使实现变得复杂。这个代码在我的机器上运行大约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
加速代码并不是为了去除for循环。我见过很多情况，矢量化代码比循环等效代码慢。在过去的20年里，MATLAB的循环变得越来越快，它们不再是减速的重要来源。例如，对于100万个元素的求和，以下仅比sum(x)慢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。
x0-x_input和y0-y_input，已经将计算时间缩短了50-60%。
此外，当在另一个循环中使用代码时，它帮助我区分哪些更改和哪些只能计算一次。在我的情况下，x0、x_input、y0、y_input始终保持不变，z_input随着每一次新的迭代而变化。所以我把计算分为两部分：第一，计算一次所有不变的东西；其次，计算所需的值。
为了能够用简单的矩阵乘法进行第二次计算，我将x-y值从矩阵(网格(重新排列为向量。这就是我所做的：
% 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(:))))
对结果的一些评论：
原始代码：在我的家用电脑上大约122秒(比我的办公室电脑的原始帖子慢得多(
Cris提出的代码：~50s
上面的代码：第一部分约75秒~第二部分为0.6秒。
替代计算的相对误差约为1E-17。
因此，在多次重复计算的情况下，预先计算没有变化的内容是有回报的。缺点是存储变量CCD_ 15所需的存储器使用量较大。
感谢您的反馈。