bugGNU Octave - Bugs: bug #57113, fminunc return wrong Hessian matrix

 
 

bug #57113: fminunc return wrong Hessian matrix

Submitter:  None
Submitted:  Thu 24 Oct 2019 09:51:57 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Flynn_Cz Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * 5.1.0
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 29 Oct 2019 09:34:02 PM UTC, comment #10: 

fminunc is implemented in an m-file so experiments are pretty easy to make.  Note this code at the end of the function:


  ## When info != 1, recalculate the gradient and Hessian using the final x.
  if (nargout > 4 && (info == -1 || info == 2 || info == 3))
    grad0 = grad;
    if (has_grad)
      [fval, grad] = fcn (reshape (x, xsz));
      grad = grad(:);
    else
      grad = __fdjac__ (fcn, reshape (x, xsz), fval, typicalx, cdif)(:);
    endif

    if (nargout > 5)
      ## Use the damped BFGS formula.
      y = grad - grad0;
      sBs = sumsq (w);
      Bs = hesr' * w;
      sy = y' * s;
      theta = 0.8 / max (1 - sy / sBs, 0.8);
      r = theta * y + (1-theta) * Bs;
      hesr = cholupdate (hesr, r / sqrt (s' * r), "+");
      hesr = cholupdate (hesr, Bs / sqrt (sBs), "-");
    endif
    ## Return the gradient in the same shape as x
    grad = reshape (grad, xsz);
  endif

  ## Restore original shapes.
  x = reshape (x, xsz);

  if (nargout > 3)
    output.iterations = niter;
    output.successful = nsuciter;
    output.funcCount = nfev;
  endif

  if (nargout > 5)
    hess = hesr'*hesr;
  endif


So there is specific code that calculates the Hessian at the final value of x.  One issue, is that the gradient and Hessian aren't calculated when info == -3 which is also an error condition.  Maybe the test condition should be updated to "info != 1".

Second, the calculation uses the damped BFGS formula and the damping will mean that the Hessian only slowly changes over to the true value.  It seems like this would be the location to insert a different Hessian calculation algorithm.

Rik <rik5>
Group administrator
Tue 29 Oct 2019 09:32:36 PM UTC, comment #9: 

@Rik No wait, I said that fminunc returns the wrong hessian non welder_mead_min; I just used welder_mead_min to estimate a good starting point for fminunc!
Bye

comment #7:

> @Giuseppe: Could you file a separate bug report about the welder_mead_min function returning an incorrect result.  That is different from this report which is focussed on fminunc.

Giuseppe <flynn_cz>
Tue 29 Oct 2019 09:28:44 PM UTC, comment #8: 

I think there are some improvements to be made in Octave.  In Matlab, running


[x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), ones (2,1))


shows that "info = 1" which means the algorithm finished successfully, and the number of iterations recorded in output is equal to 1.  For this case,


hess = [NaN NaN
        NaN NaN]


If I change the initial starting point x0,


[x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), 10*ones (2,1))


the info field is still 1 which indicates a successful completion, and the number of iterations is now 2.  In this case, the calculated Hessian is


hess = [2.0000 0.0000
        0.0000 2.0000]


This seems like "good behavior" on the part of Matlab.  When it cannot calculate the Hessian properly, it returns a matrix of NaN to indicate that the value should not be trusted.  When it is possible to calculate the Hessian (> 1 iteration) then it calculates something close to the correct value.

When you run the original example in Octave, you find that the return code info = 2.  This is not the same thing as the algorithm having converged successfully at a solution.  According to the documentation for fminunc, this return code means "Last relative step size was less than 'TolX'".  Similarly, if I start with an initial position that is even farther away


[x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), 30*ones (2,1))


I end up with an info return code of -3 which means "The trust region radius became excessively small".

It may be a false trail, but investigating why Octave is arriving at a different return code could be useful.

Rik <rik5>
Group administrator
Tue 29 Oct 2019 08:02:18 PM UTC, comment #7: 

@Giuseppe: Could you file a separate bug report about the welder_mead_min function returning an incorrect result.  That is different from this report which is focussed on fminunc.

Rik <rik5>
Group administrator
Mon 28 Oct 2019 08:01:32 PM UTC, comment #6: 

Ok, but note that using welder_mead_min from optim package to get the starting point for fminunc, this return the correct parameters but a diagonal matrix for their variance!! Look at this

par =

   2.0179e+07   2.5534e+06

V =

Diagonal Matrix

   1.0000e+00            0
            0   1.0000e+00


I can assure that V, that is the inverse of fisher matrix, is totally incorrect!


comment #4:

> For minimizing a quadratic form like @(x) sum (x.^2), BFGS [1] converges in one iteration from any starting point, but the approximate Hessian is generally not the right value when it's initially set to the identity (which is the case in fminunc). This is expected. If you need the Hessian for something like estimating confidence intervals, you should compute it at the minimum in some other way.
>
> We could explore initializing the Hessian more intelligently in fminunc using, e.g., the ideas here [2]. These would give the right answer using BFGS for this particular problem, and could conceivably be what Matlab does. fminunc in Octave actually combines BFGS with a double dogleg trust region minimization (which is different from Matlab's method, based on their documentation). Although unfortunately there's no reference provided in the fminunc code to the algorithm used, there are some details in an old mailing list thread [3]. But in general, you shouldn't count on the approximate Hessian from this sort of solver to be accurate.
>
> [1] https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
> [2] https://scicomp.stackexchange.com/q/11324
> [3] https://octave.1599824.n4.nabble.com/algorithm-for-fminunc-td4648956.html

Giuseppe <flynn_cz>
Mon 28 Oct 2019 07:54:06 PM UTC, comment #5: 

Hi all, this is an example of how fminsearch diverges.

Using Octave


x0 =

   2.1686e+07   3.2452e+06

fminsearch(@(x) loglik_LEVD(x, M_max), x0)
Exceeded target...quitting
ans =

  -5.6748e+07  -7.4920e+04


whilst, using MATLAB

x0 =

   2.1686e+07   3.2452e+06

fminsearch(@(x) loglik_LEVD(x, M_max), x0)
ans =

   2.0179e+07   2.5534e+06


loglik_LEVD is just the function that have to be minimised in order to get the parameter of the LEVD (Gumbel) distribution



comment #3:

> Matlab online gets the hessian right if x0 is sufficiently large (otherwise the Hessian is NaN(2)):


> [x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), 10*ones (2,1))
>
> hess =
> [2.0000    0.0000
>  0.0000    2.0000]


>
> I'm not an expert in numerical mathematics. But maybe there is something we can do to at least come closer to Matlab's results?
>
> The Hessian returned by Octave seems to strongly depend on x0.
>
> Please, close again if you think what we have is good enough.
>
> @Giuseppe: Can you please file another bug report if you have an example for which fminsearch diverges?

Giuseppe <flynn_cz>
Mon 28 Oct 2019 07:41:08 PM UTC, comment #4: 

For minimizing a quadratic form like @(x) sum (x.^2), BFGS [1] converges in one iteration from any starting point, but the approximate Hessian is generally not the right value when it's initially set to the identity (which is the case in fminunc). This is expected. If you need the Hessian for something like estimating confidence intervals, you should compute it at the minimum in some other way.

We could explore initializing the Hessian more intelligently in fminunc using, e.g., the ideas here [2]. These would give the right answer using BFGS for this particular problem, and could conceivably be what Matlab does. fminunc in Octave actually combines BFGS with a double dogleg trust region minimization (which is different from Matlab's method, based on their documentation). Although unfortunately there's no reference provided in the fminunc code to the algorithm used, there are some details in an old mailing list thread [3]. But in general, you shouldn't count on the approximate Hessian from this sort of solver to be accurate.

[1] https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
[2] https://scicomp.stackexchange.com/q/11324
[3] https://octave.1599824.n4.nabble.com/algorithm-for-fminunc-td4648956.html

Nir Krakauer <nir_krakauer>
Mon 28 Oct 2019 03:52:04 PM UTC, comment #3: 

Matlab online gets the hessian right if x0 is sufficiently large (otherwise the Hessian is NaN(2)):

[x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), 10*ones (2,1))

hess =
[2.0000    0.0000
 0.0000    2.0000]


I'm not an expert in numerical mathematics. But maybe there is something we can do to at least come closer to Matlab's results?

The Hessian returned by Octave seems to strongly depend on x0.

Please, close again if you think what we have is good enough.

@Giuseppe: Can you please file another bug report if you have an example for which fminsearch diverges?

Markus Mützel <mmuetzel>
Group administrator
Fri 25 Oct 2019 08:37:03 AM UTC, comment #2: 

Hi! I don’t need symbolic solution. I was using the maximul likelihood method in order to estimates parameters for the pdf of experimental data and looking at the hessian, I noted there was incorrect results. So I tried with the paraboloid and I saw that the function return the totally incorrect result.
I think that there are some errors in the functions for the optimization of multi variable functions, in fact I also noted that fminsearch often diverges despite a good starting point.

Bye!
comment #1:

> Thank you for your bug report.  Please note, that Octave "primarily intended for numerical computations".  It it not a computer algebra system https://en.wikipedia.org/wiki/Computer_algebra_system.  Octave uses some adapted BFGS (a quasi-Newton) method in "fminunc" (see "edit fminunc") and only provides approximations of the Hessian matrix.
>


> >> [x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), ones(2,1))
> x =
>
>   -0.0000000058530
>   -0.0000000058530
>
> fval =    6.8516e-17
> info =  2
> output =
>
>   scalar structure containing the fields:
>
>     iterations =  16
>     successful =  7
>     funcCount =  29
>
> grad =
>
>   -0.000000026607
>   -0.000000026607
>
> hess =
>
>    1.50000   0.50000
>    0.50000   1.50000


>
> If you want to get a symbolic solution, please use the symbolic toolbox of Octave, which does a great job for this toy problem:
>


> >> pkg load symbolic
> >> syms x y
> >> f(x, y) = x^2 + y^2
> f(x, y) = (symfun)
>
>    2    2
>   x  + y
>
> >> grad(x, y) = [diff(f, x); diff(f, y)]
> grad(x, y) = (symfun)
>
>   ⎡2⋅x⎤
>   ⎢   ⎥
>   ⎣2⋅y⎦
>
> >> hess(x, y) = [diff(f, x, 2), diff(diff (f, x), y); diff(diff (f, y), x), diff(f, y, 2)]
> hess(x, y) = (symfun)
>
>   ⎡2  0⎤
>   ⎢    ⎥
>   ⎣0  2⎦


Giuseppe <flynn_cz>
Fri 25 Oct 2019 01:44:21 AM UTC, comment #1: 

Thank you for your bug report.  Please note, that Octave "primarily intended for numerical computations".  It it not a computer algebra system https://en.wikipedia.org/wiki/Computer_algebra_system.  Octave uses some adapted BFGS (a quasi-Newton) method in "fminunc" (see "edit fminunc") and only provides approximations of the Hessian matrix.


>> [x, fval, info, output, grad, hess] = fminunc (@(x) sum (x.^2), ones(2,1))
x =

  -0.0000000058530
  -0.0000000058530

fval =    6.8516e-17
info =  2
output =

  scalar structure containing the fields:

    iterations =  16
    successful =  7
    funcCount =  29

grad =

  -0.000000026607
  -0.000000026607

hess =

   1.50000   0.50000
   0.50000   1.50000


If you want to get a symbolic solution, please use the symbolic toolbox of Octave, which does a great job for this toy problem:


>> pkg load symbolic
>> syms x y
>> f(x, y) = x^2 + y^2
f(x, y) = (symfun)

   2    2
  x  + y

>> grad(x, y) = [diff(f, x); diff(f, y)]
grad(x, y) = (symfun)

  ⎡2⋅x⎤
  ⎢   ⎥
  ⎣2⋅y⎦

>> hess(x, y) = [diff(f, x, 2), diff(diff (f, x), y); diff(diff (f, y), x), diff(f, y, 2)]
hess(x, y) = (symfun)

  ⎡2  0⎤
  ⎢    ⎥
  ⎣0  2⎦


Kai Torben Ohlhus <siko1056>
Group Member
Thu 24 Oct 2019 09:51:57 PM UTC, original submission:  

Hi! Minimising f(x, y) = x^2 + y^2 the function return the obvious wrong hessian matrix

H = [1.51691   0.47148
        0.47148   1.53984]

Bye!

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by flynn_cz (Posted a comment)
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-10-28 mmuetzel StatusInvalid / Not an Octave Bug Confirmed
        Open/ClosedClosed Open
    2019-10-25 siko1056 Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        StatusNone Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
        Operating SystemMac OS Any

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code