Finding close-to-linear solution

I am using fsolve within fsolve. The most time spent in my code is on the inner fsolve function which matlab calls 1.5 million times.

The inner function that needs to be solved is as follows:

function F = solveE(E,c)

F = E - log(c'*exp(0.99*E));

Here, E is an Nx1 vector and c is an NxN matrix. N is typically around 400. Values for c are all positive and real. Values for E are negative and real. Is there a faster way of solving this?

Within the outer fsolve, I call this function using

E0 = fsolve(@(E) solveE(E,c),E0,options);

One thing I am already doing is to use my solution E0 as an initial guess for the next iteration (when matlab adjusts its guess for the variable that solves the outer fsolve, which then will change the value for c).

Thanks for any suggestions.

Answer 1

[Matt J]

Pre-transpose c before the optimization to avoid repeatng the tranpose every iteration.

ct=c';

function outer(p,ct)
  E0 = fsolve(@(E) solveE(E,ct),E0,options);
end

Also, for the inner problem, supply the Jacobian of F,

 options.SpecifyObjectiveGradient=true;
 E0 = fsolve(@(E) solveE(E,ct),E0,options);
 
 
function [F,J] = solveE(E,ct)

  N=length(ct);
  expE=exp(0.99*E)';
  tmp=ct*expE(:);

  F = E - log(tmp);

  if nargout>1
  J=eye(N)-0.99*(ct.*expE)./tmp;
  end
end

[Tintin Milou]

Thanks! The transpose works, but the effect is very tine. The Jacobian makes calculations take twice as much time... Can that be?

[Matt J]

No. That should not happen.

[Matt J]

I replaced speye(N) with eye(N) above.

[Tintin Milou]

It works if I use

options = optimoptions('fsolve','SpecifyObjectiveGradient',true);

However, the resulting values for E are not the same as before. Initially, the difference is small (1e-6), but as Matlab is looking for new values for p, the differences between the two become bigger and at some point this method with the Jacobian can no longer find a solution. I tried playing around with the tolerance thresholds for fsolve, but that didn't help either.

[Matt J]

The solver was successful in the simplified test below, both in the derivative check and finding the solution. It doesn't appear that there is anything mathematically wrong with this step, at least. Perhaps your method of choosing the initial guess is suspect. Or perhaps you are not handling the case where no solution to solveE can be found. You don't appear to check the exitflag in the code you've presented to us.

ct=eye(400)*3;
E0=ones(400,1);
options=optimoptions('fsolve','SpecifyObjectiveGradient',true,...
                              'CheckGradient',true,...
                              'Display','iter');
 [E,Fopt,exitflag] = fsolve(@(E) solveE(E,ct),E0,options);

____________________________________________________________
   CheckGradients Information

Objective function derivatives:
Maximum relative difference between supplied 
and finite-difference derivatives = 4.38836e-08.

CheckGradients successfully passed.
____________________________________________________________


                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          1         474.031                        0.0109               1
     1          2         473.595              1         0.0109               1
     2          3         472.508            2.5         0.0109             2.5
     3          4         469.795           6.25         0.0108            6.25
     4          5         463.046         15.625         0.0108            15.6
     5          6         446.387        39.0625         0.0106            39.1
     6          7         406.075        97.6562         0.0101            97.7
     7          8         313.641        244.141        0.00885             244
     8          9         134.708        610.352         0.0058             610
     9         10               0        1160.64              0        1.53e+03

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

 exitflag

exitflag = 1

range= [min(Fopt) ,max(Fopt)]

range = 1×2

0 0

function [F,J] = solveE(E,ct)

  N=length(ct);
  expE=exp(0.99*E)';
  tmp=ct*expE(:);

  F = E - log(tmp);

  if nargout>1
  J=eye(N)-0.99*(ct.*expE)./tmp;
  end
end

[Matt J]

One thing I am already doing is to use my solution E0 as an initial guess for the next iteration (when matlab adjusts its guess for the variable that solves the outer fsolve, which then will change the value for c).

I don't know if this is advisable. What contingency plan do you follow when the inner optimization of solveE() fails and produces a bad E0 for the next iteration?

Also, the outer call to fsolve is calling the inner function solveE() not only when updating its guess. It is also calling it in order to assist in finite difference calculations. This means that the last solution to solveE might in fact be further away from the current E This is even more likely at points wher the inner solution E is non-differentiable as a function of the outer parameters (as was shown could happen in this comment).

A safer way of selecting E0, I think, would be to round the 0.99 up to 1 and thereby approximate the inner equations as,

E-log(Ct*exp(E))=0

which is equivalent to,

exp(E)=Ct*exp(E)

After a change of variables z=exp(E), this becomes a linear least squares problem,

min ||z - Ct*z||,   0<=z<=1

which can be solved with lsqlin().

[Tintin Milou]

That might help. The Jacobian works now and it cuts the running time in half. So that's pretty good!

Answer 2

[Catalytic]

Using fsolve inside of fsolve is a doubtful-sounding thing to do. Why not just combine the equations from solveE with the equations in your outer problem and solve a single system?

As a simpler example, instead of having something like this as your equation function...

function F=Equations(x,z0)
 
 z=fsolve(@(z) z^3-x, z0);
 
 F(1)=x+1-z;
 
end

...it is probably better to have this

function F=Equations(xz)

 x=xz(1); z=xz(2);

 F(1)=z^3-x;
 F(2)=x+1-z;

end

[Matt J]

Yes, indeed. The first formulation is equivalent to,

function F=Equations(x,z0)
 
 F=x+1-x.^(1/3);
 
end

which is non-differentiableat x=0, thus violating the differntiability assumptions of fsolve. I speculate that it wouldn't cause too much trouble in this instance, though, since the solver is unlikely to stray too close to x=0 as long as your initial guess is decently well-informed

[Tintin Milou]

The reason for doing the nested fsolve is that the outer fsolve is costly (mostly because it evaluates the null space to solve a homogenous system of equations), whereas the inner fsolve is quite quick.

I think my function is well behaved. So non-differentiability is unlikely to be a concern.

[Matt J]

whereas the inner fsolve is quite quick.

In what way? Your original post says that that's where most of the time is spent.

I think my function is well behaved. So non-differentiability is unlikely to be a concern.

How do you know that?

[Tintin Milou]

It's quick per evaluation; it's just that I had 1.5 million evaluations. Having 1.5 million evaluations of the outer loop would be infeasible.

Re 2: The economic model behind the equations tells me that the functions are well behaved; mathematically, there might be a range for which the function is no longer well behaved, but values in that range wouldn't make any economic sense.

[Matt J]

I don't know what "well-behaved" refers to here. It seems very hard to know in advance whether the result of the inner fsolve will be differentiable with respect to the outer unknowns. It certainly isn't guaranteed by the smoothness of solveE. You can see in Catalytic's example that the inner equation function @(z) z^3-x is a highly smooth, polynomial function of both z and x. Yet, solving for z leads to a non-differentiable function of x.