Векторизация нескольких циклов в Matlab / Python

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, я также радпопробуй в питоне.

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 ГБ), и поэтому он может оказаться невозможным для вашей текущей ситуации.У меня недостаточно памяти на моей текущей машине, чтобы проверить это, поэтому я не знаю, насколько быстрее это будет работать, но теоретически это должно быть еще быстрее из-за отсутствия каких-либо циклов for в коде.