Wrong Hessian output in fminunc

Can somebody explain weird Matlab's output for the following optimization problem? The goal is to minimize and to evaluate Hessian matrix at the solution. However, Matlab returns a matrix of NaNs for the Hessian at the solution point. Here is the code

[xOpt,~,~,~,grad,hessian] = fminunc(@(x) x(1)^2+x(2)^2,[1,1])

Answer 1

[Matt J]

Another solution would be to use standalone routines for numerical gradient and hessian estimation

https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

fun=@(x) (x(1)^2+x(2)^2);


xOpt = fminunc(fun ,[1,1]);

g=gradest(fun,xOpt)
H=hessian(fun,xOpt)

No hidden scale factors if you go that route:

xOpt =

     0     0


g =

     0     0


H =

    2.0000         0
         0    2.0000

[Alex]

Thank you! I will accept your asnwer! One question though. Is there a similar way (scripts available) when the objective function is not in analytic form?

[Matt J]

There is no requirement in the solution I've given you that the objective be in analytic form.

[Alex]

Isn't is explicitly stated in the documentation for the files that you referenced?

Answer 2

[Bruno Luong]

I guess the hessian is estimated from BFGS formula, that needs more than 1 FMINUNC iteration. In you case the minimum is reached in 1 iteration (since the gradient point directly to the minimum), the hessian estimation is not yet in the update cycle.

A small modification make # iteration > 1, and hessian starts well.

[xOpt,~,~,~,grad,hessian] = fminunc(@(x) 2*x(1)^2+x(2)^2,[1,1])

[Alex]

Thank you for the explanation. However, you modified the function to minimize and the Hessian is not the same as for the original problem. Do you know if there is any other way in Matlab to obtain a Hessian matrix at a point?

[Bruno Luong]

Not 100% satisfied answer but if you use trust-region algorithm and suppy the gradient, you 'll get back the hessian even with 1 iteration

options = optimoptions('fminunc', ...
    'Algorithm', 'trust-region', ...
    'SpecifyObjectiveGradient', true ...
    );

[xOpt,out,b,c,grad,hessian] = fminunc(@fun,[1,1],options)

function [f,g] = fun(x)
% Calculate objective f
f = x(1)^2+x(2)^2;
if nargout > 1 % gradient required
    g = 2*x;
end
end

If your function is quadratic and the -gradient points toward the minimum at any point, meaning if FMINUNC converge in 1 iteration regardless the starting point, then you must have hessian that is a scale of eye matrix. The NAN likely happens only for the trivial case as you have showed.

Answer 3

[Matt J]

Just run fmincon (twice) with no constraints. Running a second time is important, because the hessian output is not the Hessian calculated at the final point, but rather the Hessian at the point just prior to that.

fun=@(x) x(1)^2+x(2)^2;

xOpt = fmincon(fun ,[1,1]);

[xOpt,~,~,~,~,grad,hessian] = fmincon(fun,xOpt)

xOpt =

   1.0e-08 *

    0.5588    0.5588


grad =

   1.0e-07 *

    0.2608
    0.2608


hessian =

     1     0
     0     1

[Bruno Luong]

Hmm, wasn't the hessian

[2 0;
 0 2]

?

[Matt J]

Hmmm, indeed. No, it wasn't. In fact, I find that multiplying the function by any scalar leaves the hessian output completely unchanged. Some internal pre-normalization nof the cost function, I suppose..

[Matt J]

Alex's comment moved here:

Isn’t it explicitly stated in the documentation?

[Matt J]

Not as far as I can see. The documentation on the meaning of the Hessian output is here,

https://www.mathworks.com/help/optim/ug/hessian.html

but there is no mention that I can see of any pre-scaling of the objective.