Added use of second order derivative to Newton-Rhapson Method #20

Closed
wants to merge 1 commit into
from

Projects

None yet

6 participants

@seppwerk

As described in http://projects.scipy.org/scipy/ticket/1433 now via github.

For my own work, I implemented the use of the second derivative in the Newton-Raphson-Method. This reduced the number of NR-steps on average more than one third.

This might be pretty useful, if the calculation of the function is very expensive (i.e. in case of having a PDE system) and the derivatives are very cheap (i.e. because of AD-Tools are used anyway).

In a diagram attached to the ticket, a brief explanation can be found.

@chrisjordansquire chrisjordansquire commented on the diff Jul 21, 2011
scipy/optimize/zeros.py
@@ -42,12 +42,13 @@ def results_c(full_output, r):
# Newton-Raphson method
-def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50):
+def newton(func, x0, fprime=None, fsec=None, args=(), tol=1.48e-8, maxiter=50):
"""Find a zero using the Newton-Raphson or secant method.
@chrisjordansquire
chrisjordansquire Jul 21, 2011

I think this should have a bit more of an explanation of the algorithm in the notes section, since 2nd order Newton-Raphson solvers don't seem to be very common.

@scopatz
scopatz Jul 25, 2011

A reference would probably suffice.

@scopatz
SciPy member

Also, this pull request could use a unit test or two... I would be nice to put this option in.

@pv pv commented on the diff Aug 29, 2011
scipy/optimize/zeros.py
@@ -105,7 +111,17 @@ def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50):
msg = "derivative was zero."
warnings.warn(msg, RuntimeWarning)
return p0
- p = p0 - func(*myargs) / fder
+ fval = func(*myargs)
+ if fsec is not None:
+ fder2 = fsec(*myargs)
+ if fder2 == 0:
+ p = p0 - fval / fder
+ else:
+ root = fder ** 2 - 2 * fval * fder2
+ if root < 0:
+ p = p0 - fder / fder2
+ else:
+ p = p0 - (fder - sign(fder) * (root ** 0.5)) / fder2
@pv
pv Aug 29, 2011

This loses precision when fder2 -> 0: you end up doing (a - b)/c in floating point with b ~ a + c and |c| << |a|, |b|. The formula can be written with a + sign by using a different expression for the root.

@pv pv commented on the diff Aug 29, 2011
scipy/optimize/zeros.py
@@ -105,7 +111,17 @@ def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50):
msg = "derivative was zero."
warnings.warn(msg, RuntimeWarning)
return p0
- p = p0 - func(*myargs) / fder
+ fval = func(*myargs)
+ if fsec is not None:
+ fder2 = fsec(*myargs)
+ if fder2 == 0:
+ p = p0 - fval / fder
+ else:
+ root = fder ** 2 - 2 * fval * fder2
+ if root < 0:
+ p = p0 - fder / fder2
@pv
pv Aug 29, 2011

This could use a reference (it's not immediately obvious why do this). Would it be better to fall back to the 1st order method?

@pv
pv Aug 29, 2011

Ah ok, the point seems to be to jump to the minimum/maximum.

@pv
SciPy member
pv commented Aug 29, 2011

This should go in Scipy 0.10, just needs a simple test case (and checking the loss of precision issue).

@rgommers
SciPy member

Please move fsec=None to the end of the list of parameters before committing, otherwise it breaks backwards compatibility.

@WarrenWeckesser

How close is this algorithm to Halley's method (http://en.wikipedia.org/wiki/Halley%27s_method)? A quick look shows that Halley's method avoids taking a square root.

@seppwerk

I also realized some days ago, that the proposed method is pretty similar to the rational Halley’s method, or, more precise, it is exactly the same as the originally proposed parabolical Halley’s method (http://de.wikipedia.org/wiki/Halley-Verfahren has some more explanations, but in German).

The general difference is, that the rational Halley’s method uses a hyperbolic approximation of the function, instead of a parabolic (quadratic). The convergence of both is cubic.

I want to examine the differences between both, since the rational Halley’s method does sometimes show strange convergence behaviour.

Halley’s method is a special case of the Householder’s method (http://en.wikipedia.org/wiki/Householder%27s_method), so maybe, it makes more sense to implement this method in a more general way with explicit special cases, instead of filling the newton-method with so much more than newton’s method (it already implements the secant method).

@pv
SciPy member
pv commented May 19, 2012

Ok, I changed the formula to address my comments, added the reference to the parabolic Halley, added some necessary tests, and merged this. We can as well have one cubic order method in newton.

The Householder family of methods might be a useful addition, if it turns out that someone wants to implement them.

@pv pv closed this May 19, 2012
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment