0

I have an equation f(x)=y with the follwing properties:

  • we have an explicit expression for the derivative f'
  • the root is known to be between 0 and 1 (in fact f is defined only in this interval)
  • f is (weakly) increasing

I am currently solving this equation with scipy.optimize.brentq, which works well enough.

optimize.brentq(f_to_zero, 0, 1)

brentq uses the secant method, a finite-difference equivalent of Newton's method; it does not make use of any explicit f'.

I was hoping I could make the root finding faster by making use of this derivative. In other words, I want to use Newton's method augmented with bounds; this article suggests the following idea: if at any point the Newton guess falls outside the bounds, perform a bisection (between the current guess and the bound) instead. Is there any well-tested package for this? I would really rather not write my own implementation (it's probably not worth it for the performance gain).

Another point is that the function isn't even defined outside the bounds, so using

optimize.newton(f_to_zero, fprime=f_prime, x0=0.5)

is not just inefficient, it will throw an error

Somebody asked what the expression for f' is. f and f' are very long, so I have made gist here: https://gist.github.com/d3c48dde09389e8c48da0e990b57bf99

tmkadamcz
  • 119
  • 6

1 Answers1

1

If you're constraining yourself to scipy optimizers, you can transform this into a minimization problem with objective function f(x)2 and use the fmin_tnc optimizer, which uses Newton's method augmented with bounds:

fsq = lambda x : f(x) ** 2
fsqprime = lambda x : 2 * f(x) * fprime(x) 
optimize.fmin_tnc(fsq, x0 = [0.5], fprime = fsqprime, bounds = [(0,1)])
Cecil Cox
  • 376
  • 1
  • 7
  • Thanks! After trying this, to my surprise it was actually slower than brent's method. For 100 values, I get 'newton_tnc' 92.79 ms, 'brent' 6.95 ms. Looking at the RootResults, brent typically needs 15-18 function calls and as many iterations, whereas TNC needs only 3 or 4 iterations but the same number of function calls. For TNC I don't know if 'function calls' means calls to the Jacobian or the function itself, or the sum of both. – tmkadamcz Feb 14 '21 at 18:19