Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

newton_raphson_iterate bisection fallback does not work #808

Open
wthrowe opened this issue Jul 26, 2022 · 5 comments · May be fixed by #1000
Open

newton_raphson_iterate bisection fallback does not work #808

wthrowe opened this issue Jul 26, 2022 · 5 comments · May be fixed by #1000

Comments

@wthrowe
Copy link

wthrowe commented Jul 26, 2022

The documentation for newton_raphson_iterate states that "Out-of-bounds steps revert to bisection of the current bounds", which should be robust, but the implementation fails to bisect correctly in some cases. The following program requests the root of a quartic with the sole root 1 in the bounding interval, but with an initial guess on the wrong side of an extremum that causes a step out of bounds. The returned result is approximately the endpoint of the interval farthest from the actual root, which should have been eliminated by a single bisection.

#include <boost/math/tools/roots.hpp>
#include <iostream>
#include <utility>

int main() {
  const auto quartic = [](const double x) {
    return std::make_pair((x * x - 1.0) * (x * x + 0.1),
                          2.0 * x * ((x * x - 1.0) + (x * x + 0.1)));
  };

  const double result =
      boost::math::tools::newton_raphson_iterate(quartic, 0.5, 0.1, 1.1, 5);
  std::cout << result << "\t" << quartic(result).first << "\n";
  return result;
}

Output:

0.10625 -0.110033
@jzmaddock
Copy link
Collaborator

Unfortunately, I don't think this is solvable in all cases.

In this situation, the initial NR step goes out of bounds to the left, we therefore assume (erroneously in this case) that there is a root in the interval [0.1, 0.5] and bisect that interval. Subsequent NR steps succeed, but only lead us to the left endpoint 0.1.

Now, we could have bisected the whole interval [0.1, 1.1] and in this case this would have led to the correct root, but consider a curve that looks more like:

nr

With P1 as the initial guess. This time we go out of bounds to the right, and bisection of [p1, b] leaves us in the correct domain, where as bisecting [a,b] would have led to a dead end (and another bug).

Basically, if we are given an initial guess we the right to assume that it is "downhill" from that guess to the root. For example, one could easily construct curves that look a lot like yours, but which have one or more roots between [0.1, 0.5] in addition to the root on [0.5,1.1], it would be a reasonable expectation in that case that we would follow the direction indicated by the derivative at the initial guess.

I realise this sounds like a cop-out, and I really am committed to making this routine as robust as possible, not least because NR routines have a reputation for exploding in all sorts of situations. But I just don't think it's possible to make them completely fool-proof in all situations.

Suggestions welcome!

@NAThompson
Copy link
Collaborator

@wthrowe : Does the Halley iterate succeed? In general I find the Halley iterate has a more sane basin of convergence and moreover it can degrade gracefully to a newton iterate.

@jzmaddock
Copy link
Collaborator

Does the Halley iterate succeed? In general I find the Halley iterate has a more sane basin of convergence and moreover it can degrade gracefully to a newton iterate.

It doesn't go out-of-bounds with the first step, but since we're starting the wrong side of the maxima, it just goes downhill away from the root until it hits the search endpoint. Which is what I would expect to be honest.

I guess the question is to what lengths we're prepared to go to protect against bad initial guesses?

@NAThompson
Copy link
Collaborator

I guess the question is to what lengths we're prepared to go to protect against bad initial guesses?

Yeah, I think the worry is that if you add too many branches for these cases, you either create bugs, or create performance issues.

Though maybe "if it goes out the right bound, start the guess at the left" might not be insane?

@wthrowe
Copy link
Author

wthrowe commented Jul 28, 2022

There are simple robust algorithms for this (see, for example, rtsafe from Numerical Recipes), but after converting our project to use one of them I don't think you can switch because it would break too much downstream code. Since they can fall back to an actual bisection algorithm, they require that the passed bounds are a valid bracket for a bisection search, which is currently not required for your newton_raphson. (Infinite bounds can probably be handled, but I've already hit cases where evaluating f(bound) crashes.)

wthrowe added a commit to wthrowe/spectre that referenced this issue Aug 1, 2022
The Boost version's bisection fallback does not work properly.

boostorg/math#808
wthrowe added a commit to wthrowe/spectre that referenced this issue Aug 3, 2022
The Boost version's bisection fallback does not work properly.

boostorg/math#808
@ryanelandt ryanelandt linked a pull request Jul 12, 2023 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants