在 Matlab/Python 中向量化多个 for 循环
Vectorizing multiple for loops in Matlab/Python
我正在尝试编写一个数学模型,它涉及在值网格上计算特定数量数千次,并使用一些不断变化的模型参数。目前,这太慢了,我正在寻找有关矢量化模型中最密集部分的建议。
为了便于阅读,我目前已经有了它的基本实现,但现在如果可能的话,我想对下面的整个代码段进行矢量化。代码段的一个最小示例是:
% 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 中尝试。
通过将 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 循环在代码中。
我正在尝试编写一个数学模型,它涉及在值网格上计算特定数量数千次,并使用一些不断变化的模型参数。目前,这太慢了,我正在寻找有关矢量化模型中最密集部分的建议。
为了便于阅读,我目前已经有了它的基本实现,但现在如果可能的话,我想对下面的整个代码段进行矢量化。代码段的一个最小示例是:
% 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 中尝试。
通过将 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 循环在代码中。