Strange behavior of fmincon when analytical gradient is supplied

Dear all

I want to find the minimum of a log likelihood function of the form

[f,gr_f] = mylogl(theta, data)

subject to linear constraints, where theta is 3x1 vector of the log likelihood parameters and gr_f is the analytical gradient. As I see it, in order to find the minimum, I can employ three methods

1. Use fmincon with numerical derivatives
2. Use fmincon with analytical derivatives
3. Use fsolve to solve the system of equations grad_f = 0 ( I have a separate function, named grad_f that computes the gradient of the log likelihood, gr_f )

The thing is that methods (1) and (3) provide the same results, and if they dont, method (3) provides better results, however, I cannot make option (2) to work. When the analytical gradient is supplied, fmincon always gets stuck to the initial point, no matter what that is. This is the output of the procedure:

if true
Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                     step    optimality
  0           1   -1.846070e+03     0.000e+00     1.000e+00     0.000e+00     1.458e+03  
  1          63   -1.846070e+03     0.000e+00     2.489e-10     1.163e-10     1.458e+03  
Local minimum possible. Constraints satisfied.
fmincon stopped because the size of the current step is less than
the selected value of the step size tolerance and constraints are 
satisfied to within the selected value of the constraint tolerance.
<stopping criteria details>
Elapsed time is 17.091739 seconds.
Optimization stopped because the relative changes in all elements of x are
less than options.StepTolerance = 1.000000e-10, and the relative maximum constraint
violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-12.
Optimization Metric                                           Options
max(abs(delta_x./x)) =   9.21e-11                       StepTolerance =   1e-10 (selected)
relative max(constraint violation) =   0.00e+00   ConstraintTolerance =   1e-12 (selected)

The gradient vector supplied by the fmincon equals

   evalmodel.gradient = [-1.458070270701287e+03;-3.697754312033735e+02]

and my optimization options are as follows:

options = optimoptions('fmincon','Algorithm','sqp','Display','iter',...
'MaxFunctionEvaluations',1000,'SpecifyObjectiveGradient',true);
options = optimoptions(options,'ConstraintTolerance',1e-12,'StepTolerance',1e-10);

Any ideas?

Answer 1

[John D'Errico]

The very first thing I would do is use the 'CheckGradients" option for fmincon. Is it just vaguely possible that the gradient is incorrect as supplied? :) Yes, it is, and given the problem you have reported, the test is trivial to do, and a logical choice.

[Manthos Vogiatzoglou]

Dear John

Thanks for the suggestion. I avoided doing as suggested for two main reasons. First, since solving the equation "gradient = 0" provided the same or better results with numerical fmincon, doesn't that mean that my gradient is correct? And second, one of my parameters has value around .01 and the other can be as high as 200. I got worried about the gradient check, which has fixed value 10^-6.

[John D'Errico]

.01 and 200 is not too large of a range, although I'd often advise scaling your variables if it was any larger.

You could also my derivest tools (downloaded from the file exchange) to verify the gradient computation. Since it will also return an estimate of the error in the computation, you can then decide if the gradient computation is correct.

Finally, nothing stops you from doing a test yourself of the gradient. A derivative is pretty trivial to compute as simply a finite difference approximation. Pick your own choice of dx, and use a central difference for accuracy.