sqp behavior for functions with small values

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

sqp behavior for functions with small values

Benjamin Schwab
All:

    I am optimizing a likelihood function which contains really small probabilities (on the order of 10^-35) and want to use Octave to help me find optimal parameters that maximize the likelihood function.  Using sqp on a function with such small values appears not to do anything and returns the guess that I give it.  After getting some advice on where to seek help in the IRC channel I posted this apparent issue to the bug tracker.  It was closed as not being a bug but I was advised to post the issue to the mailing list in an effort to improve Octave.  The exact inputs I can use to recreate the situation from a fresh start of Octave is in the post script.  I should mention that I'm running Octave 4.2.2 on Windows 10 home and when I installed Octave, I was warned that it might behave incorrectly on Windows 10.  I have some background with mathematics but almost none with programing.

    The suggested change of taking the log of the likelihood function works.  As I can change the tolerance for the sqp function to something that should still work with the unmodified function, this led me to believe that sqp would work.  If sqp cannot work on functions that take on such small values then maybe that should be noted in the documentation.  In any event, it should be tested to see if indeed this is an issue and if it is an issue on other operating systems before changing the documentation.  If desired, I could test sqp on Windows 10 to figure out what the true minimum tolerance is for sqp.  I don't have time at the moment to do so.

    I hope this is useful to someone else.  I don't know what more details could be requested but I would be happy to provide them or if I need to clarify my thoughts for you'all.

In Gratitude,
Benjamin A Schwab

Post Script:
>> ttt = @(x) x(1)^2+x(2)^2
ttt =

@(x) x (1) ^ 2 + x (2) ^ 2

>> sqp([70,40],ttt,[],[],[-500,-500],[500,500],200,1e-40)
ans =

   0
   0

>> ggg=@(x) -factorial(27)*prod(1/sqrt(2*pi*x(2)^2).*exp(-([1,2,3,7,20,26,35,52,58,60,64,69,72,77,84,90,90,101,104,111,121,121,132,139,141,155,15
9].-x(1)).^2./(2*x(2)^2)))
ggg =

@(x) -factorial (27) * prod (1 / sqrt (2 * pi * x (2) ^ 2) .* exp (-([1, 2, 3, 7, 20, 26, 35, 52, 58, 60, 64, 69, 72, 77, 84, 90, 90, 101, 104, 11
1, 121, 121, 132, 139, 141, 155, 159] - x (1)) .^ 2 ./ (2 * x (2) ^ 2)))

>> ggg([70,40])
ans =  -2.9010e-035
>> ggg([77,47])
ans =  -1.1955e-034
>> ggg([77.543,47.707])
ans =  -1.2047e-034
>> sqp([70,40],ggg,[],[],[-500,0],[500,500],200,1e-40)
ans =

   70
   40

>> bbb=@(x) -log(factorial(27)*prod(1/sqrt(2*pi*x(2)^2).*exp(-([1,2,3,7,20,26,35,52,58,60,64,69,72,77,84,90,90,101,104,111,121,121,132,139,141,15
5,159].-x(1)).^2./(2*x(2)^2))))
bbb =

@(x) -log (factorial (27) * prod (1 / sqrt (2 * pi * x (2) ^ 2) .* exp (-([1, 2, 3, 7, 20, 26, 35, 52, 58, 60, 64, 69, 72, 77, 84, 90, 90, 101, 10
4, 111, 121, 121, 132, 139, 141, 155, 159] - x (1)) .^ 2 ./ (2 * x (2) ^ 2))))

>> bbb([70,40])
ans =  79.525
>> bbb([77,47])
ans =  78.109
>> bbb([77.543,47.707])
ans =  78.102
>> sqp([70,40],bbb,[],[],[-500,0],[500,500])
ans =

   77.556
   47.691

>>


Reply | Threaded
Open this post in threaded view
|

Re: sqp behavior for functions with small values

Juan Pablo Carbajal-2
Hi Benjamin,

I am happy to hear that the -log solved this particular issue. The
problem you are facing is not a problem of sqp, but of qp.
If you inspect the results of line 412 of sqp, you see that the search
direction is 0,0, qp did 2 iterations and thinks it got a global
solution.
qp.m is just calling __qp__ defined in qp.cc [1]

in line 97 qp uses a fixed tolerance for the search given by sqrt(eps).

This is a bug in qp and one would need to pass the tolerance as a
parameter (or compute the correct tolerance for qp based on the
arguments of sqp). I am adding the author fo the __qp__.cc in CC to
see if he can fix this quickly. I also created a new bug report [2]
(thanks!).

[1]: http://octave.org/doxygen/4.3/da/dbd/____qp_____8cc_source.html
[2]: https://savannah.gnu.org/bugs/index.php?53506