在Matlab / Python中向量化多个for循环

Question

我正在尝试对数学模型进行编码，它涉及在数千个具有变化的模型参数的值网格上计算特定数量。目前，这太慢了，我正在寻找有关向量化模型中最密集部分的建议。

为了便于阅读，我目前已经有了它的基本实现，但是现在如果可能的话，希望对下面的整个代码段进行矢量化处理。该代码段的一个最小示例是：

% Setup grid to evaluate and results vector
T_max = 10000;
eval_points = linspace(0, T_max, 1000);
results = zeros(size(eval_points));
% Function that is used in computation
Z_func = @(x, omega) (1./(omega.*sqrt(2*pi))).*exp( -(x.^2)./(2.*omega.*omega) );
% Random data for now, known in full problem
historic_weights = rand(1,100);
historic_times   = rand(1,100);
% Fixed single parameter omega
omega            = 0.5;
% Time evaluation
tic()
for eval_counter = 1:size(eval_points,2)
    for historic_counter = 1:size(historic_weights,2)
    temp_result = 0;
        for k = 0:1:T_max
            temp_result = temp_result + Z_func( eval_points(eval_counter) - historic_times(historic_counter) + 1440*floor(historic_times(historic_counter)/1440) - 1440*k, omega );
        end % End of looping over k
        results(eval_counter) = results(eval_counter) + historic_weights(historic_counter)*temp_result;
    end % End of looping over weights 
end % End of looping over evaluation points
toc()

在我的计算机上，这只花了60多秒钟来评估。我不希望使用并行工具箱，因为我已经在其他地方使用了它，并且显示的代码段在每个进程中都被调用。

如果在Matlab中无法做到这一点，我很乐意也尝试使用python。

Answer 1

您可以通过计算temp_result和result将矩阵作为一次矩阵而不是一次简单地向量化内部两个循环。例如：

for eval_counter = 1:size(eval_points,2)
    temp_result = sum(Z_func( eval_points(eval_counter) - historic_times + 1440*floor(historic_times/1440) - 1440*(0:1:T_max)', omega ));
    results(eval_counter) = results(eval_counter) + sum(historic_weights.*temp_result);
end % End of looping over evaluation points

在我的机器上，此过程只需要9秒钟，而循环版本的运行时间是73秒钟。

现在，从理论上讲，您可以在没有单个循环的情况下执行此操作，如下所示：

eval_points = linspace(0,T_max,1000);
historic_weights = rand(100,1); % Note transposed from original
historic_times   = rand(100,1);
eval_loop = reshape(0:T_max,1,1,[]); % size = [1,1,10000];

result = sum(historic_weight.*sum(Z_func(eval_points - historic_times + 1440*floor(historic_times/1440) - 1440*eval_loop, omega ),3),1);

然而，这将使用一个显著的内存量（> 8 GB），因此可能不适合你现在的情况是可行的。我当前的机器上没有足够的内存来测试它，所以我不知道它的运行速度有多快，但是从理论上讲，由于代码中没有任何for循环，它甚至应该更快。