## Vectorizing multiple for loops in Matlab

### Kieran Kalair (view profile)

on 20 Sep 2019
Latest activity Commented on by Rik

on 20 Sep 2019

### darova (view profile)

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.

R2019a

### darova (view profile)

on 20 Sep 2019

I changed your code a bit See attached script

Kieran Kalair

### Kieran Kalair (view profile)

on 20 Sep 2019
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

### Rik (view profile)

on 20 Sep 2019
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