存储值后MATLAB全局搜索速度骤降的原因及优化咨询

问题背景 #

我定义了一个定义域为0≤x≤π、0≤y,z≤2π的函数，编写了步长为0.01的MATLAB代码进行全局搜索寻找最大值点，同时存储每个点的函数值。示例代码如下：

x=0;
y=0;
z=0;
i=1;

Maxdeviation= ComputeDeviationSum(0,0,0,mufT, mufD, mufC);
MaximizingArg= [x, y,z] ;
while x<= pi
    y=0;
    z=0;
    while y<= 2*pi
        z=0;
        while z<=2 *pi
            deviation= x+y+z; % 仅为示例，实际函数更复杂
            if deviation> Maxdeviation
                MaximizingArg= [x, y,z];
                Maxdeviation= deviation;
            end
            a(i,:)=[x, y,z, deviation];
            z=z+0.01;
            i=i+1;
        end
        y=y+0.01;
    end
    x=x+ 0.01;
end

运行后发现代码速度极慢，删除以下存储相关代码后，代码运行速度大幅提升（20秒内迭代次数提升至原来的10倍）：

i=1;
a(i,:)=[x, y,z, deviation];
i=i+1;

现咨询两个问题：

1. 为何删除存储值的代码后速度大幅提升？
2. 若需保留存储每个点函数值的功能，该如何优化代码以提升运行速度？

问题1：删除存储代码后速度大幅提升的原因 #

• 动态数组扩容的巨大开销：MATLAB中数组a初始未预分配内存，每次执行a(i,:)=[...]时，MATLAB需要重新分配更大的内存空间，把原有数据完整复制到新空间，再写入新元素。你的迭代次数极多（步长0.01时总点数约1.97e8），这种频繁的内存重分配和数据复制会累积成致命的性能瓶颈。
• 额外的临时操作成本：每次赋值a(i,:)=[x,y,z,deviation]需要临时创建1×4的向量，再执行写入操作，这也会增加额外的计算和内存操作损耗。

问题2：保留存储功能的优化方案 #

方案1：预分配数组内存（最核心的优化） #

提前计算总迭代次数，为数组a分配足够的固定内存，彻底避免动态扩容的开销：

% 计算各维度的总点数
x_steps = round(pi / 0.01) + 1;
y_steps = round(2*pi / 0.01) + 1;
z_steps = round(2*pi / 0.01) + 1;
total_points = x_steps * y_steps * z_steps;

% 预分配4列的数组（指定类型可进一步优化内存效率）
a = zeros(total_points, 4, 'double');

x=0;
y=0;
z=0;
i=1;

Maxdeviation= ComputeDeviationSum(0,0,0,mufT, mufD, mufC);
MaximizingArg= [x, y,z] ;
while x<= pi
    y=0;
    while y<= 2*pi
        z=0;
        while z<=2 *pi
            deviation= ComputeDeviationSum(x,y,z,mufT, mufD, mufC); % 替换为实际函数
            if deviation> Maxdeviation
                MaximizingArg= [x, y,z];
                Maxdeviation= deviation;
            end
            a(i,:)=[x, y,z, deviation];
            z=z+0.01;
            i=i+1;
        end
        y=y+0.01;
    end
    x=x+ 0.01;
end

方案2：向量化运算替代嵌套循环（进一步提速） #

MATLAB原生擅长向量化操作，把嵌套循环改成矩阵运算可大幅减少循环开销：

% 生成所有网格点
[x_grid, y_grid, z_grid] = meshgrid(0:0.01:pi, 0:0.01:2*pi, 0:0.01:2*pi);
% 将网格展平为列向量
x_vec = x_grid(:);
y_vec = y_grid(:);
z_vec = z_grid(:);

% 向量化计算所有点的函数值
% 若原函数支持向量化输入（接受数组参数、元素级运算），直接调用即可；否则用arrayfun
deviation_vec = arrayfun(@(x,y,z) ComputeDeviationSum(x,y,z,mufT, mufD, mufC), x_vec, y_vec, z_vec);

% 合并为存储数组
a = [x_vec, y_vec, z_vec, deviation_vec];

% 寻找最大值点
[Maxdeviation, max_idx] = max(deviation_vec);
MaximizingArg = [x_vec(max_idx), y_vec(max_idx), z_vec(max_idx)];

提示：如果能修改ComputeDeviationSum使其支持向量化输入，速度会比arrayfun更快。

方案3：并行计算工具箱（超大规模数据场景） #

如果总点数过大，可借助并行计算工具箱将循环分配到多个CPU核心执行：

% 开启并行池
parpool;

% 预分配数组
x_steps = round(pi / 0.01) + 1;
y_steps = round(2*pi / 0.01) + 1;
z_steps = round(2*pi / 0.01) + 1;
total_points = x_steps * y_steps * z_steps;
a = zeros(total_points, 4);

Maxdeviation= ComputeDeviationSum(0,0,0,mufT, mufD, mufC);
MaximizingArg= [0,0,0];

% 用parfor替代外层循环，分配到多核心
parfor x_idx = 1:x_steps
    x = (x_idx-1)*0.01;
    for y_idx = 1:y_steps
        y = (y_idx-1)*0.01;
        for z_idx = 1:z_steps
            z = (z_idx-1)*0.01;
            deviation= ComputeDeviationSum(x,y,z,mufT, mufD, mufC);
            % 并行环境下最大值更新需保证线程安全
            if deviation > Maxdeviation
                Maxdeviation = deviation;
                MaximizingArg = [x,y,z];
            end
            % 计算全局索引
            global_idx = (x_idx-1)*y_steps*z_steps + (y_idx-1)*z_steps + z_idx;
            a(global_idx,:) = [x,y,z,deviation];
        end
    end
end

% 关闭并行池
delete(gcp);

注意：需确保ComputeDeviationSum是线程安全的，避免全局变量读写冲突。

内容的提问来源于stack exchange，提问作者Asigan