Avoid evaluating the function with same input multiple times

Question

I am trying to utilise scipy.optimise.fsolve for solving a function. I noticed that the function is evaluated with the same value multiple times in the beginning and the end of the iteration steps. For example when the following code is evaluated:

from scipy.optimize import fsolve

def yy(x):
    print(x)
    return x**2+9*x+20

y = fsolve(yy,22.)

print(y)

The following output is obtained:

[ 22.]
[ 22.]
[ 22.]
[ 22.00000033]
[ 8.75471707]
[ 4.34171812]
[ 0.81508685]
[-1.16277103]
[-2.42105811]
[-3.17288066]
[-3.61657372]
[-3.85653348]
[-3.96397335]
[-3.99561793]
[-3.99984826]
[-3.99999934]
[-4.]
[-4.]
[-4.]

Therefore the function is evaluated with 22. three times, which is unnecessary.

This is especially annoying when the function requires substantial evaluation time. Could anyone please explain this and suggest how to avoid this issue?

Answer 1

The first evaluation is done only to check the shape and data type of the output of the function. Specifically, fsolve calls _root_hybr which contains the line

shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))

Naturally, _check_func calls the function:

res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))

Since only the shape and data type are retained from this evaluation, the solver will be calling the function with the value x0 again when actual root finding process begins.

The above accounts for one extraneous call (out of two). I did not track down the other one, but it's conceivable that the FORTRAN code does some kind of preliminary check of its own. This sort of thing happens when algorithms written long ago get wrapped over and over again.

If you really want to save these two evaluations of expensive function yy, one way is to compute the value yy(x0) separately and store it. For example:

def yy(x):
    if x == x0 and y0 is not None:
       return y0
    print(x)
    return x**2+9*x+20

x0 = 22.
y0 = None
y0 = yy(x0)
y = fsolve(yy, x0)