MATLAB Answers

0

Vectorizing multiple for loops in Matlab

Asked by Kieran Kalair on 20 Sep 2019
Latest activity Commented on by Rik
on 20 Sep 2019
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.

  0 Comments

Sign in to comment.

Products


Release

R2019a

1 Answer

Answer by darova
on 20 Sep 2019
 Accepted Answer

I changed your code a bit
See attached script

  2 Comments

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
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

Sign in to comment.