Is this a bug in GSL's minimum finding routine?

Question

I have run into a very strange problem using GSL recently. I'm trying to find the maximum of an ugly function by asking GSL to find the minimum of the function's negative image. This has been working fine with most of the functions I've been using it on up to now, but with certain functions, I get an error.

Specifically, the ruby GSL gem throws an exception when I feed my function to the minimum-finder, stating that the given endpoints do not enclose a minimum. However, the given endpoints DO enclose a minimum; furthermore, the results seem to depend on the initial estimate given.

When I ask GSL to find the minimum starting with a guess of 0.0, I get the error. When I ask it to find the minimum starting with a guess of 2.0, it finds a minimum. It seems strange that GSL would complain that there isn't a minimum in a given interval just depending on the initial guess.

Below is a script I've written to reproduce this error. I'm using GSL version 1.15 and the latest version of the Ruby/GSL gem wrapper, with Ruby version 1.9.3p392.

Will someone please try this script and see if they can reproduce these results? I would like to report this as a bug to the GSL maintainers, but I'd feel better doing that if I could get the error to occur on somebody else's computer too.

There are three values for the initial minimum guess provided in the script; starting with 0.0 causes GSL to throw the aforementioned error. Starting with 1.0 causes GSL to report a minimum, which happens to only be a local minimum. Starting with 2.0 causes GSL to report another minimum, which seems to be the global minimum I'm looking for.

Here is my test script:

require("gsl")
include GSL::Min

function = Function.alloc { |beta|  -(((6160558822864*(Math.exp(4*beta))+523830424923*(Math.exp(3*beta))+1415357447750*(Math.exp(5*beta))+7106224104*(Math.exp(6*beta)))/(385034926429*(Math.exp(4*beta))+58203380547*(Math.exp(3*beta))+56614297910*(Math.exp(5*beta))+197395114*(Math.exp(6*beta))))-((1540139705716*(Math.exp(4*beta))+174610141641*(Math.exp(3*beta))+283071489550*(Math.exp(5*beta))+1184370684*(Math.exp(6*beta)))*(1540139705716*(Math.exp(4*beta))+174610141641*(Math.exp(3*beta))+283071489550*(Math.exp(5*beta))+1184370684*(Math.exp(6*beta)))/(385034926429*(Math.exp(4*beta))+58203380547*(Math.exp(3*beta))+56614297910*(Math.exp(5*beta))+197395114*(Math.exp(6*beta)))**2)) }

def find_maximum(fn1)
  iter = 0;  max_iter = 500
  minimum = 0.0             # reasonable initial guess; causes GSL to crash!
  #minimum = 1.0             # another initial guess, gets a local min
  #minimum = 2.0             # this guess gets what appears to be the global min
  a = -6.0
  b = 6.0
  #pretty wide interval

  gmf = FMinimizer.alloc(FMinimizer::BRENT)

  gmf.set(fn1, minimum, a, b)
  #THIS line is failing (sometimes), complaining that the interval given doesn't contain a minimum. Which it DOES.

  begin
    iter += 1
    status = gmf.iterate
    status = gmf.test_interval(0.001, 0.0)
    # puts("Converged:") if status == GSL::SUCCESS
    a = gmf.x_lower
    b = gmf.x_upper
    minimum = gmf.x_minimum
    # printf("%5d [%.7f, %.7f] %.7f %.7f\n",
           # iter, a, b, minimum, b - a);
  end while status == GSL::CONTINUE and iter < max_iter
  minimum
end

puts find_maximum(function)

Please let me know what happens when you try this code, commenting out the different initial values for minimum. If you can see a reason why this is actually the intended behavior of GSL, I would appreciate that as well.

Thank you for your help!

Answer 1

I think this Mathematica notebook perfectly summarize my answer

Whenever you need to numerically minimize a particular function, you must understand it qualitatively (make plots) to avoid failure due to numerical errors. This plot shows a very important fact: Your function seems to require numerically difficult cancellations between very large numbers at beta ~ 6. Numerically issues can also arise at beta ~ -6.

Then, when you set the interval to be [-6, 6], your numerical errors can prevent GSL to find the minimum. A simple plot and few evaluations of the potentially dangerous terms can give a better understanding of the appropriate limits you must set in GSL routines.

Conclusion: if the minimum is around zero (always make plots of the functions you want to minimize!), then you must set a small interval around zero to avoid numerical problems when you have functions that depends on precise cancellation between exponentials.

Edit: In a comment you asked me to not consider the factor Abs[beta]. In this case, GSL will find multiple local minimum (=> you must diminish the interval) . I also can't understand how you get the minimum around 1.93

Answer 2

GSL one-dimensional minimising routines report this error whenever your initial guess has a function value greater or equal than the function value at either extreme of the interval. So, indeed, it depends on the initial guess that you find one or another local minimum, or none at all. The error message is certainly misleading.