Avoid evaluating a function with the same input at the same time
I am trying to use scipy.optimise.fsolve to solve a function. I noticed that the function is evaluated with the same value multiple times at the beginning and 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)
This produces the following output:
[ 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 from 22. three times, which is optional.
This is especially annoying when a feature requires significant evaluation time. Can anyone explain this and suggest how to avoid this problem?
source to share
The first evaluation is performed only to validate the shape and data type of the function's output. Specifically, fsolve
calls _root_hybr
that contain the line
shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
Naturally what _check_func
calls the function :
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
Since only the shape and data type are preserved from this estimate, the solver will call the function with the value again x0
when the actual process of finding the root begins.
The above means one outside call (out of two). I haven't tracked down the other, but it's entirely possible that the FORTRAN code is doing some pre-check. This happens when algorithms written a long time ago are completed over and over again.
If you really want to keep those two estimates of the expensive function yy, one way is to compute the value of 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)
source to share
I realized that an important reason for this question is that fsolve is not designed for such a problem: Solvers should be chosen wisely :)
multivariate: fmin, fmin_powell, fmin_cg, fmin_bfgs, fmin_ncg
nonlinear: leastsq
constrained: fmin_l_bfgs_b, fmin_tnc, fmin_cobyla
global: basinhopping, brute, differential_evolution
local: fminbound, brent, golden, bracket
n-dimensional: fsolve
one-dimensional: brenth, ridder, bisect, newton
scalar: fixed_point
source to share