0

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!

karmic_mishap
  • 129
  • 5

2 Answers2

0

I think this Mathematica notebook perfectly summarize my answerenter image description here

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) enter image description here . I also can't understand how you get the minimum around 1.93

Vivian Miranda
  • 2,467
  • 1
  • 17
  • 27
  • Thank you for your assistance. I have indeed plotted the function, and, when I plot it without subtracting it from beta the way you have, I get that the global minimum is at about 1.93434. May I ask why you are plotting it the way you are? One concern I have about this is that I am actually generating these functions from a program and trying to minimize them right then and there, so I generally don't plot them and adjust things manually. I suppose there are few enough of them that I can, and you're saying it is a good idea, which makes sense. – karmic_mishap Sep 15 '13 at 02:17
  • Further, when I do try to minimize the function you present by putting a "beta" before the negative sign in my above function, I still can't find this root near zero, even if I make the interval as small as [-0.000000002,0.000000002]. Mathematica does find a very small root, but I'm writing this using GSL because I want to distribute my code, so I'd like to use open source software as much as possible to broaden my potential user-base. – karmic_mishap Sep 15 '13 at 02:27
  • I don't know how to program in Ruby, but it seems the expression in Mathematica was exact the function you sent to GSL minimization. That is why I plotted it. II - if it is an automatic function generation - you must coded in a way to avoid numerical cancellation between large numbers. – Vivian Miranda Sep 15 '13 at 02:28
  • Ah, well, I'm actually looking for the maximum of the function you've plotted above, which is why the negative sign was in the original post. The maximum of your above post still _looks_ to be about 1.93, to me. – karmic_mishap Sep 15 '13 at 02:53
  • Ah, now I understand what happened with that part. What I sent GSL to minimize was the function without the Abs[beta] part; the |beta| you see in the above Ruby code is just telling Ruby the name of the parameter I'm using. Sorry for the confusion. – karmic_mishap Sep 15 '13 at 02:57
  • I think I get the idea, though. Generally, then, I can't just set a general interval like [-6,6], even if all of my functions should have their minima in that interval, because of issues like this. Therefore, I should plot them before setting the interval for each one, as you've suggested. It's more work, and less automatic, but it should still be something I can get others to replicate later,if they're interested. Thanks again! – karmic_mishap Sep 15 '13 at 02:59
  • But you see that you can't include [-6,6] because it includes a lot of local minimum/maximum. – Vivian Miranda Sep 15 '13 at 02:59
  • If there is too much functions to plot by hand: grid the interval first to find a 10% approximation of the minimum, then apply gsl – Vivian Miranda Sep 15 '13 at 03:05
  • Sorry I can't take it to chat, since I don't have enough reputation on this particular site. May I ask you what "gridding the interval" means? Is there a good reference you could point me towards? I plan to go back to read the appropriate chapter of Numerical Recipes again to see if there's anything on it. – karmic_mishap Sep 15 '13 at 04:15
  • I meant just to make a grid (example : split your guess interval in 100 points and evaluate the function in each of this points) to have 10% estimate of the location of the minimum/maximum . Then you can set gsl interval to be ~ few times distance between 2 grid points. – Vivian Miranda Sep 15 '13 at 04:28
0

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.

Aureli Alabert
  • 11
  • 1