Matrix input in fzero solver

Hello,

I am stuck with my code. I am trying to vectorize loops for performance issues. However, I have a problem when the vectorization involves fzero. Here I post an example of what I want to get. First code line shows the script with the loops, and the second one the function file.

MAIN SCRIPT

clear all
clc

grida=10;
gridz=60;
w=1.296;
mu1=0.85;
theta=0.1;
s = rng;

bbar_AQ = zeros(grida,gridz);
pistar_n_f = zeros(1,gridz);
k_star_f = zeros(1,gridz);
auxiliar = zeros(grida,gridz);
debt = zeros(grida,gridz);
A=randn(10,60);

aa=0.001:(w-0.002)/(grida-1):w-0.001;
z = linspace(1, 12, gridz);

 opt=optimset('display', 'off'); 

 for j=1:grida
     for t=1:gridz
          k_star_f(t) =(z(t))^(1/(1-mu1-theta));
          pistar_n_f(t) = (1-mu1)*(mu1/w)^(mu1/(1-mu1))*(z(t))^(1/(1-mu1))*(k_star_f(t)^(theta/(1-mu1))); 
          bbar_AQ(j,t) = (pistar_n_f(t) + aa(j))/0.75;
          auxiliar(j,t) =k_star_f(t);

          if bbar_AQ(j,t)>=auxiliar(j,t)
              debt(j,t)=A(j,t);
          else 
                   if j==1
                       bbar_0 = 0.1;
                   else
                       bbar_0 = 0.11;
                   end
                         debt(j,t)=fzero(@(x)b_max(x,z(t),aa(j),theta,mu1,w),bbar_0,opt);
          end     
     end
 end

 

FUNCTION FILE

 function f=b_max(x,z,a,theta,mu,w)

f=(((1-mu)*(mu/w)^(mu/(1-mu))*(z*(x+a)^theta)^(1/(1-mu))))-x;

The above code runs fine. Howeve, I am trying to vectorize it but when I arrive to the fzero line of code, I am stuck about how to proceed. This is what I did:

clear all
clc

grida = 10;
gridz = 60;
w = 1.296;
mu1 = 0.85;
theta = 0.1;


aa = 0.001:(w-0.002)/(grida-1):w-0.001;
z = linspace(1, 12, gridz);


rng(s)
A=randn(10,60);

% Compute k_star_f and pistar_n_f using vectorized operations
k_star_f = z.^(1./(1-mu1-theta));
pistar_n_f = (1-mu1)*(mu1/w).^(mu1./(1-mu1)).*(z.^(1./(1-mu1))).*(k_star_f.^(theta./(1-mu1)));

% Compute bbar_AQ using vectorized operations
expanded_aa = repmat(aa(:), 1, gridz);
bbar_AQ = (pistar_n_f + expanded_aa) / 0.75;

% Compute auxiliar using vectorized operations
auxiliar = k_star_f;
auxiliar=repmat(auxiliar,grida,1);

% Compute debt using vectorized operations based on the condition
debt2 = A2 .* (bbar_AQ2 >= auxiliar2) +  fzero(@(x)b_max(x,z(t),aa(j),theta,mu1,w),bbar_0,opt).* (bbar_AQ2 < auxiliar2);

Any help would be really appreciated.

Thank in advanced!

[Torsten]

If you want to do the calls to "fzero" in a parfor loop - ok, this might save time. But writing the calls in the vectorized form as you do is nonsense. The time is spent in "fzero", not in the loop from which "fzero" is called.

[Ibai Ostolozaga Falcon]

I see. So, do you recommend me to run a loop over the values for fzero? And what about cellfun? Thank you!

[Torsten]

Yes.

Maybe you can save time by choosing bbar_0 as the solution debt of the last call to fzero. Good initial guesses for the nonlinear solver are the best means for fast convergence.

[Ibai Ostolozaga Falcon]

Yes, thank so much! I will try it with parfor loop using the solution debt of the last call to fzero as bbar_0.

[Torsten]

To use the solution of the last step within a parfor loop will not be possible. Either the loop is processed sequentially (with prescribing bbar_0 from the last call) or the loop is processed in parallel (without the possibility of prescribing bbar_0 from the last call).

[Ibai Ostolozaga Falcon]

I see. So, this is what I tryed:

clear all
clc
grida = 10;
gridz = 60;
w = 1.296;
mu1 = 0.85;
theta = 0.1;
aa = 0.001:(w-0.002)/(grida-1):w-0.001;
z = linspace(1, 12, gridz);
rng(s)
A=randn(10,60);
% Compute k_star_f and pistar_n_f using vectorized operations
k_star_f = z.^(1./(1-mu1-theta));
pistar_n_f = (1-mu1)*(mu1/w).^(mu1./(1-mu1)).*(z.^(1./(1-mu1))).*(k_star_f.^(theta./(1-mu1)));
% Compute bbar_AQ using vectorized operations
expanded_aa = repmat(aa(:), 1, gridz);
bbar_AQ = (pistar_n_f + expanded_aa) / 0.75;
% Compute auxiliar using vectorized operations
auxiliar = k_star_f;
auxiliar=repmat(auxiliar,grida,1);

bbar_constrained = zeros(grida,gridz);

      for j=1:grida
           
        for t=1:gridz
          

         if j==1 & t==1
            bbar_0 = 0.1;
         else
            bbar_0 = bbar_0;
         end

            bbar_constrained(j,t)=fzero(@(x)b_max_2_OnlyPi(x,z(t),eta,aa(j),theta,mu1,w),bbar_0,opt);  
            bbar_0=bbar_constrained(j,t)
        end    
                 
      end
      
debt = A*(bbar_AQ>=auxiliar)+bbar_constrained*(bbar_AQ<auxiliar) ;
      

So, it will not possible to use parfor loop?

Thank you in advanced!

[Torsten]

So, it will not possible to use parfor loop?

No. The loop as written has to be processed sequentially because you made the initial values of the i-th call to "fzero" depend on the result of the (i-1)th call.

You might try to set bbar_0 outside and only leave the line

bbar_constrained(j,t)=fzero(@(x)b_max_2_OnlyPi(x,z(t),eta,aa(j),theta,mu1,w),bbar_0,opt);

within the double loop.

This loop should then be possible to be processed in parallel.

[Ibai Ostolozaga Falcon]

Thank you so much for your help Torsten!