让我们首先对循环和条件进行智能化处理，因为您将循环索引用作条件变量。
我们从~(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进行循环。例如，在第一次迭代中，将k设置为1:5。在第二次迭代中，将其设置为6:10，依此类推。你需要把每次迭代的结果加到上一次迭代中，形成最终的A。

说明

此函数可通过隐式展开进行矢量化（这比使用meshgrid更有效率），并使用.^和.*等元素运算子来取代^和*运算子。（因为我们有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  ---------------

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下可能不会更快，因为执行引擎可以优化原始代码中的循环。（它确实在Octave下提供了加速。）

在循环中处理if语句的一个技巧是将if语句（或其中的一部分）转换为逻辑矩阵。然后，您可以将逻辑矩阵 * 元素方式 * 乘以您在每一步中添加的值的矩阵。false值将计算为零，并且不会改变结果。
这只在每个元素都可以独立计算的情况下才有效。这通常会使实际的计算线变慢，但在Matlab中，这往往被删除for循环所带来的巨大速度提升所抵消。
在您的示例中，我们可以使用meshgrid来消除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。