Asked by Kieran Kalair
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.

Answer by darova
on 20 Sep 2019

Accepted Answer

I changed your code a bit

See attached script

Kieran Kalair
on 20 Sep 2019

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.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.