bug #57113: fminunc return wrong Hessian matrix

 Submitted by: None Submitted on: 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

You are not logged in

Discussion

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))
[fval, grad] = fcn (reshape (x, xsz));
else
grad = __fdjac__ (fcn, reshape (x, xsz), fval, typicalx, cdif)(:);
endif

if (nargout > 5)
## Use the damped BFGS formula.
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
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>
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

> @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>
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>
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!

> 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

> 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.

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>
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
>
>
>   -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:
>

> >> 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)]
>
>   ⎡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

-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:

>> 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)]

⎡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>
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

Attached Files

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

Attach Files:

Comment:

No files currently attached

Dependencies

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)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.

History

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