Vectorizing multiple for loops in Matlab

I am trying to code a mathematical model, and it involves computing a particular quantity over a grid of values thousands of times, with some changing model parameters. Currently, this is far too slow and I am looking for advice on vectorizing the most intensive part of my model.

The equation I am evaluating is:

I've currently got a basic implementation of it for ease of reading, but now want to vectorize the entire code segment below if possible. A minimal example of the code segment is:

% 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)
    temp_result = 0;
    for historic_counter = 1:size(historic_weights,2)
        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()

On my computer, this took just over 60 seconds to evaluate. I do not wish to use the parallel toolbox, as I am already using that elsewhere, and the shown segment of code is called on every process.

Answer 1

[darova]

I changed your code a bit

See attached script

[Kieran Kalair]

Thanks, there was a typo that you spotted so your solution gives the correct results. RAM is not too much of an issuie so it even scales to quite large arrays for my avaliable hardware.

[Rik]

For those not wanting to download the script to see the solution:

clc,clear

% Setup grid to evaluate and results vector
T_max = 5;
t = linspace(0, T_max, 100);
res = zeros(size(t));
% 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
hw = rand(1,10);
ht = rand(1,10);
% Fixed single parameter omega
omega = 0.5;
% Time evaluation
tic()
for i = 1:size(t,2)
    for j = 1:size(hw,2)
        tres = 0;
        for k = 0:1:T_max
            x = t(i) - ht(j) + 1440*floor(ht(j)/1440) - 1440*k;
            tres = tres + Z_func( x, omega );
        end % End of looping over k
        res(i) = res(i) + hw(j)*tres;
    end % End of looping over weights
end % End of looping over evaluation points
toc()

% VECTORIZE VERSION
tic
[T, HT, K] = ndgrid(t,ht, 0:T_max);
[~, HW] = ndgrid(t,hw);
X = T - HT + 1440*floor(HT/1440) - 1440*K;
tres = Z_func(X, omega);
tres = HW.*sum(tres,3);
res1 = sum(tres,2);
toc

res' - res1