Vectorizing multiple for loops in Matlab/Python

Question

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.

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

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.

If this is not possible in Matlab, I'm happy to also try in python.

MrAzzaman · Accepted Answer · 2019-09-20 12:07:34Z

5

You can fairly easily vectorise the inner two loops by computing temp_result and result as matrices instead of one at a time. For instance:

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

This runs in ~9 seconds on my machine, compared to 73 seconds for your looped version.

Now, in theory you can do this without a single loop, as follows:

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

However this will use a significant amount of memory (>8 GB), and so may not be feasible for your current situation. I don't have sufficient memory on my current machine to test this, so I don't know how much faster this would run, but in theory it should be even faster, due to the lack of any for-loops in the code.

edited Sep 20, 2019 at 12:07

answered Sep 20, 2019 at 11:52

MrAzzaman

4,79815 silver badges28 bronze badges

Sign up to request clarification or add additional context in comments.

8 Comments

Kieran108 Over a year ago

The first set of code runs in ~3 seconds on my machine, a huge improvement already. I do have access to hundreds of GB of RAM, so I tried the second code, but on the result = line I get an error 'Matrix dimensions must agree.' Do you have any idea if something should be transposed and should not be?

MrAzzaman Over a year ago

Oh, sorry. I had forgotten to replace (0:1:T_max)' with eval_loop. I've fixed it, and it should work now.

HansHirse Over a year ago

Did you check for equivalent results? If I run the original script against your first solution for T_max = 10 and eval_points = linspace(0, T_max, 10), the two results differ (quite heavily). Maybe, I'm doing something wrong - that's why I'm asking if you could reproduce the original results!?

MrAzzaman Over a year ago

I will double check.

Kieran108 Over a year ago

There was actually a typo in my original code, the variable temp_result was not set to 0 before each loop over k. MrAzzaman's answer was actually the correct, both vectorized and fixing this typo. I've changed the question to reflect this now

|

Collectives™ on Stack Overflow

Vectorizing multiple for loops in Matlab/Python

1 Answer 1

8 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

8 Comments

Your Answer

Sign up or log in

Post as a guest

Related