我有一个带有5个嵌套的matlab程序
for
环路和
if
条件如下：
for x0=1:N
for y0=1:N
for k=1:N
for x1=1:N
for y1=1:N
if ~((y1-x1>N/2)||(x1-y1>N/2)) && ~((y0-x0>N/2)||(x0-y0>N/2))
A(x0,y0)=A(x0,y0)+2^(k*((x0-y0)+(x1-y1)))*B(x1,y1)
end
end
end
end
end
end
其中A和B是两个矩阵。如何使此程序运行得更快？
我试过使用meshgrid，但似乎不起作用，因为有
if
条件。
让我们先聪明地处理循环和条件，因为您使用循环索引作为条件变量。
我们从~(y1-x1>N/2)||(x1-y1>N/2)，或方式更清晰，abs(y1-x1)<N/2。
与其有if条件，为什么不强制y1始终在范围内？
最后一个循环可以写为y1=max(x1-N/2,1):min(x1+N/2,N)，因此不需要if条件的第一部分的全部。当然，我们可以对其他变量做同样的事情：
for x0=1:N
for y0=max(x0-N/2,1):min(x0+N/2,N)
for k=1:N
for x1=1:N
for y1=max(x1-N/2,1):min(x1+N/2,N)
A(x0,y0)=A(x0,y0)+2^(k*((x0-y0)+(x1-y1)))*B(x1,y1)
end
end
end
end
end
现在，为了清楚起见，让我们对k进行重新排列和矢量化。它没有必要是中间循环，事实上，它作为中间循环的唯一特征是混淆阅读代码的人。但除此之外，它也没有必要成为一个循环。
k=1:N;
for x0=1:N
for y0=max(x0-N/2,1):min(x0+N/2,N)
for x1=1:N
for y1=max(x1-N/2,1):min(x1+N/2,N)
A(x0,y0)=A(x0,y0)+sum(2.^(k*((x0-y0)+(x1-y1))))*B(x1,y1)
end
end
end
end
现在，速度快吗？
没有MATLAB非常擅长优化代码，所以它不会更快。但至少它的方式更清晰，所以我想你已经做到了。但你需要更快！好我不确定你能不能。你有5个嵌套的循环，这太慢了。我不认为你可以meshgrid这个，即使没有条件，因为你混合了所有的循环。meshgrid在网格网格上执行操作时很好，但在您的情况下，您对每个x0,y0使用所有x1,y1，因此它不是网格操作。
这里有一个矢量化的解决方案：
x0 = (1:N).';
y0 = 1:N;
x1 = (1:N).';
y1 = 1:N;
k = reshape(1:N, 1, 1, N);
conditiona = ~((y0-x0 > N/2) | (x0-y0 > N/2));
conditionb = ~((y1-x1 > N/2) | (x1-y1 > N/2));
a = 2 .^ (k .* (x0-y0)) .* conditiona;
b = 2 .^ (k .* (x1-y1)) .* B .* conditionb;
bsum = squeeze(sum(sum(b, 1) ,2));
A = A + reshape(reshape(a, [] , N) * bsum ,N ,N);
注意，创建可能/可能不需要大量存储器的两个3D阵列a和b。在这种情况下，您需要在k上循环。例如在第一次迭代中将CCD_ 15设置为CCD_。在第二次迭代中，将其设置为6:10，依此类推。您需要将每次迭代的结果与上一次迭代相加，以形成最终的A。
解释
该函数可以通过隐式展开(比使用meshgrid更有效)和使用像.^和.*这样的元素运算符来矢量化，而不是使用^和*运算符。结果创建了5D阵列(因为我们有5个循环变量)，该阵列应该在第3-5维上求和以产生最终的2D矩阵。然而，这可能需要大量内存。另一点是，包含乘积和的函数通常可以写成有效的矩阵乘法
表达式：
2^(k*((x0-y0)+(x1-y1)))*B(x1,y1);
可以写成：
2 .^ (k .* (x0-y0)) .* 2 .^ (k .* (x1-y1)) .* B(x1, y1)
------- a -------- ------------- b ---------------
这是两个子表达式的乘积，每个子表达式都有3个维度，因为每个子表达式只包含3个循环变量。因此，5D问题被简化为3D问题。if条件还有两个子表达式，每个子表达式都可以乘以a和b子表达式：
conditiona = ~((y0-x0 > N/2) | (x0-y0 > N/2));
conditionb = ~((y1-x1 > N/2) | (x1-y1 > N/2));
a = 2 .^ (k .* (x0-y0)) .* conditiona;
b = 2 .^ (k .* (x1-y1)) .* B .* conditionb;
只需使用两个循环变量x0和y0:就可以形成for循环
for x0=1:N
for y0=1:N
A(x0,y0)=A(x0,y0)+ sum(sum(sum(a(x0,x0, :) .* b, 3), 2), 1);
%or simply A(x0,y0)=A(x0,y0)+ sum(a(x0,x0, :) .* b, "all");
end
end
通过预先计算b:的和，可以将其简化为以下循环
bsum = sum(sum(b, 1) ,2);
% bsum = sum(b ,[1, 2]);
for x0=1:N
for y0=1:N
A(x0,y0)=A(x0,y0)+ sum(a(x0,x0, :) .* bsum, 3);
% or as vector x vector multiplication
% A(x0,y0)=A(x0,y0)+ squeeze(a(x0,x0, :)).' * squeeze(bsum);
end
end
这里可以通过使用矩阵x向量乘法来防止循环：
A = A + reshape(reshape(a, [] , N) * bsum ,N ,N);
更新在Matlab下，这个解决方案可能不会更快，因为执行引擎可以优化原始代码中的循环。(它确实在倍频程下提供了加速。)
在循环中处理if语句的一个技巧是将if语句(或其一部分)转换为逻辑矩阵。然后，您可以将逻辑矩阵按元素乘以您在每个步骤中添加的值的矩阵。false值将计算为零，并且不会更改结果。
只有当每个元素都可以独立于其他元素进行计算时，这才有效。它通常会使实际计算线变慢，但在Matlab中，这通常被去除for循环带来的巨大速度提高所抵消。
在您的示例中，我们可以将这个想法与网格一起使用，以删除x0和y0上的循环。计算线需要成为元素矩阵计算，因此元素运算符.*、.^和|取代了*、^和|。
% Warning: Y0 and X0 are swapped in this example
[Y0, X0] = meshgrid(1:N,1:N);
% Create a logical matrix which represents part of the if statement
C = ~((Y0-X0>N/2) | (X0-Y0>N/2));
for k=1:N
for x1=1:N
for y1=1:N
if ~((y1-x1>N/2)||(x1-y1>N/2))
% Multiply by C elementwise
A = A + C.*2.^(k*((X0-Y0)+(x1-y1)))*B(x1,y1);
end
end
end
end
你甚至可以进一步理解这个概念，并使用多校正ndgrid来删除更多的循环，但它会变得更复杂(你必须开始在维度上求和)，如果N很大，多维数组就会变得难以处理。
注意：注意索引顺序。meshgrid将y定义为行，将x定义为列，因此矩阵索引为A(y,x)，但您使用的是A(x,y)。为了使你的例子有效，我在meshgrid的输出中切换了x和y。