How do I speed up a loop where I need to check whether elements are less than 0

I am trying to speed-up the following loop. This loop updates a value function v_d given a guess v_ini_exp. All variables are pre-allocated. From running a profiler it seems that the major bottleneck is the line where I need to discard elements of a matrix which are less than 0. In my model, these values are not feasible. But it seems that this is taking almost 65% of the execution time. I was wondering if there is a way to speed up the loop (provided below)

I can think of two ways to try to speed it up:

(1) Speed-up the line that takes over 65% of the execution time (I tried different formulation but the current one is the fastest I could think of)

(2) Try to vectorize at least one the the loop, but this is unclear to me since the variable d in the inner loop is a matrix.

Is there any way to speed up this loop? I can use parfor to decrease the runtime by around 50% but is there another way to speed up this loop before jumping to use parallel toolbox?

% Variables defined up to this loop:
% (1) q_ini is an array (z_grid_num-by-k_grid_num-by-b_grid_num)
% (2) v_ini_exp is an array (z_grid_num-by-k_grid_num-by-b_grid_num)
% (3) wealth is an array (z_grid_num-by-k_grid_num-by-b_grid_num)
% (4) choice a matrix (k_grid_num-by-b_grid_num)

% continuation value function
for ii = 1:z_grid_num
    q_ini_ii = squeeze(q_ini(ii,:,:));
    v_ini_exp_ii = squeeze(v_ini_exp(ii,:,:));
    choice = q_ini_ii.*b_grid_choice - k_grid_choice;
    for jj = 1:k_grid_num
        for kk = 1:b_grid_num
            
            % dividends at time t
            d = wealth(ii,jj,kk) + choice;
            
            % choosing optimal consumption
            vec_d = (d + v_ini_exp_ii).*(d>0) + (-1e30).*(d<=0);  % <- THIS LINE TAKES 65% OF THE EXECUTION TIME
            [v_d(ii,jj,kk)] = max(vec_d,[],'all'); % <- THIS LINE TAKES 25% OF THE EXECUTION TIME
            
        end
    end
end

[Ameer Hamza]

I am a bit confused by your code. Why do you need this line?

vec_d = (d + v_ini_exp_ii).*(d>0) + (-1e30).*(d<=0);

which is replacing all the values less than 0 with -1e30.

However, you are using max() function, which will select the maximum value anyway. So why replace those values with -1e30.

[Walter Roberson]

You replace exact 0 with a small negative value. If it is acceptable to instead replace negative with 0 then use max(0,expression)

[Michal Szkup]

I am using the line

vec_d = (d + v_ini_exp_ii).*(d>0) + (-1e30).*(d<=0);

to make sure that I do not choose (d + v_ini_exp_ii) where d<0. It might be the case that d<0 but v_ini_exp_ii is very large and then max will pick up that number.

But your and Walter's comments, give me an idea how to speed up the look. I initialize vec_d to low number to begin with so that I only need to execute

vec_d = (d + v_ini_exp_ii).*(d>0) 

This decreased execution time by about 10-15%. Would be nice to cut it down a bit more since I will have to execute this code hunderds of time with large grids but even such small gain is useful.

[Walter Roberson]

vec_d = max(0, d + v_ini_exp_ii)

would probably be faster if it is acceptable to use 0 instead of 1e-30

[Ameer Hamza]

I guess that issue here is that the OP only wants the maximum value of this sum when d > 0. I think the OP's current solution is probably the fastest it can get in MATLAB.

[Adam Danz]

vec_d = (d + v_ini_exp_ii).*(d>0) 

is the way to go, then.

[Walter Roberson]

max((d + v_ini_exp_ii).*(d>0))

does appear to be about as fast as you are likely to get, being only a small bit slower than max(d+v_ini_exp_ii) in my tests.

[Walter Roberson]

Michal, is one of the following a real possibility ?

• v_ini_exp_ii can be quite large in a situation where d is negative -- large enough that d plus the value could end up being the maximum overall ?
• v_ini_exp_ii can be consistently enough negative where d are positive such that the maximum of the sum where d is positive could be negative ?

[Michal Szkup]

@Walter I am afraid of the first possibility. I will have to solve the problem for different parameter values so I would like to have a check against that. I cannot imagine that the second option would happen.

Apparently, writing the problematic line as

vec_d = d + v_ini_exp_ii.*(d>0) 

provides an additional very small but positive gain. This formulation is valid for my problem since the lowest value v_ini_exp_ii takes is 0.

[Michal Szkup]

I tested the code a little bit more and it seems that if I could be sure that I can execute the line

vec_d = d + v_ini_exp_ii

rather than

vec_d = d + v_ini_exp_ii.*(d>0) 

then, as @Walter suggested, this yields additional gains of up to 25% which is probably as fast as this can get without restructuring the whole code.

For the current parameters, the code returns the same solutions in both cases, but I will need to test that this holds true for a wide range of parameters.