How to make MATLAB code run faster

How do I speed this script up?

nuke;   
h = 16; %hours available
K = 15; %Capital
z = 1; %Total factor productivity
t = 0; %Taxes
G = 0; %Government Spending?
iteration = 0; % run counter
precision = 1e-3; 
distance = 2; % random value greater than precision
w = 1; %wage rate
pi_max = 1e-3; %profit?

% declare lgrid
lmin = 1e-5;
lmax = h - 1e-5;
lnum = 100000;
lgrid = linspace(lmin,lmax,lnum);

while distance > precision
    c = zeros(size(lgrid)); %init array
    u = zeros(size(lgrid)); 
    for i = 1:lnum  
        c(i) = w*(h - lgrid(i)) + pi_max - t;
        u(i) = log(c(i)) + log(lgrid(i));
    end
    [u_max, position] = max(u); 
    l_max = lgrid(position);
    N_max_supply = h - l_max;
    N = h - lgrid;
    pi = zeros(size(lgrid));
    for i = 1:lnum
      pi(i) = z*K^0.3*N(i)^0.7 - w*N(i); 
    end
    
    [pi_max, location] = max(pi);

    N_max_demand = N(location); 

  
    
    
    distance = abs(N_max_demand - N_max_supply); 
    w = 0.99999*w + 0.00001*(N_max_demand - N_max_supply);       

    if w < 0 
        stop 
    end 
    iteration = iteration + 1; % this is the number of times while loop has been executed befores topping
    s = sprintf ( ' w = %6.3f   iteration %4d    ||labor demand - labor supply|| = %8.6f      labor demand = %8.6f   labor supply = %8.6f  ', ...
        w, iteration, distance, N_max_demand, N_max_supply);
    disp(s)
end 
% now we have found the general equilibrium point, we can play with it and
% see what other values we have: 
% output:
y = z*K^0.3*N_max_demand^0.7; 
C = w*(N_max_supply) + pi_max - t;
Utility = log(C) + log(h-N_max_supply); 

[James Tursa]

Have you run the profiler? Which line(s) are the bottleneck?

[Kai Jensen]

No, im trying to parallelize the while loop I have setup.

[Rik]

It might not be possible to remove the while loop. At least I don't see an easy way to do it. If you really want to, you should try finding the direct formula. In my experience that's usually easier with pen and paper than with Matlab.

If you can't find the direct formula, you should run the profiler and look what your code spends most of the time on.

I notice the two for loops can probably be vectorized, so if that's the bottleneck, just post a comment.

[Kai Jensen]

yeah I ran the profiler, it seems like the two for loops are the bottleneck.

[Kai Jensen]

Even after vectorizing it still takes a bit of time. Any way to speed up the exponental function?

[Jan]

Of cousre you should not compute the expensive K^0.3 in each iteration, but once before the loop. Or replace the 2nd loop:

pi = z * K ^ 0.3 * N .^ 0.7 - w * N;

But do I see correctly that z * K ^0.3 * N .^ 0.7 is a constant and does not change in the while loop? Then compute this once before the loop.

By the way "pi" is a bad choice for the name of a variable, because this shadows the builtin function pi.

You calculate u and c in the first loop only to find the maximum value. Why not calcuklating the first derivative are find the formula to determine the maximum? This avoids a lot of expansive log() commands.

[Kai Jensen]

what would that look like?

Answer 1

[Jan]

I just moved some repeted calculations out of the loop:

h = 16; %hours available
K = 15; %Capital
z = 1; %Total factor productivity
t = 0; %Taxes
G = 0; %Government Spending?
iteration = 0; % run counter
precision = 1e-3; 
distance = 2; % random value greater than precision
w = 1; %wage rate
pi_max = 1e-3; %profit?
% declare lgrid
lmin = 1e-5;
lmax = h - 1e-5;
lnum = 100000;
lgrid = linspace(lmin,lmax,lnum);

logLgrid = log(lgrid);
N  = h - lgrid;
c1 = z * K ^ 0.3 * N .^ 0.7;

while distance > precision
    c = w * N + pi_max - t;
    u = log(c) + logLgrid;
    [u_max, position] = max(u); 

    l_max = lgrid(position);
    N_max_supply = h - l_max;
    
    pi = c1 - w * N;    
    [pi_max, location] = max(pi);
    N_max_demand = N(location); 
 
    distance = abs(N_max_demand - N_max_supply); 
    w = 0.99999*w + 0.00001*(N_max_demand - N_max_supply);       
    if w < 0 
        stop  % ??? What is "stop"? Do you mean: break 
    end 
    iteration = iteration + 1; % this is the number of times while loop has been executed befores topping
    s = fprintf(' w = %6.3f   iteration %4d    ||labor demand - labor supply|| = %8.6f      labor demand = %8.6f   labor supply = %8.6f\n', ...
        w, iteration, distance, N_max_demand, N_max_supply);
end 
% now we have found the general equilibrium point, we can play with it and
% see what other values we have: 
% output:
y = z*K^0.3*N_max_demand^0.7; 
C = w*(N_max_supply) + pi_max - t;
Utility = log(C) + log(h-N_max_supply); 

This version needs 7 seconds in Matlab R2016b, and 5 if the slow output to the command window is omitted, while the original took 202 seconds. All I've done was moving the repeated expensive |log and power operations out of the loop.

[Kai Jensen]

Thank you sir