I'm implementing unconstrained minimization of 2D functions:
- For the search direction, I'm using both steepest descent and Newton descent.
- For the step size, I'm using the backtracking line search algorithm
The code is very simple (Mainly for future reference to anyone who might find it helpful): My function optimize(f, df, hess_f, method) looks like this:
# prepare contour plot
fig = pt.figure()
xmesh, ymesh = np.mgrid[-5:5:100j,-5:5:100j]
fmesh = f(np.array([xmesh, ymesh]))
pt.axis("equal")
pt.contour(xmesh, ymesh, fmesh)
pt.contour(xmesh, ymesh, fmesh, 250)
# initial guess + first update of search direction
guesses = [np.array([5, 0.1])]
x = guesses[-1]
s = search_direction(df, hess_f, x, method)
not_done = True
while not_done:
# calculate step size using backtracking line search:
alpha_opt = backtracking_alpha(f, df, 0.5, 0.5, x, s)
# update step
next_guess = x + alpha_opt * s
guesses.append(next_guess)
# plot current step to our updating contour graph:
it_array = np.array(guesses)
pt.plot(it_array.T[-2], it_array.T[-1], "-")
# check stopping condition
if (np.linalg.norm(guesses[-2] - guesses[-1]) < 0.0001):
not_done = False
# prepare for next guess according to search direction method
x = guesses[-1]
s = search_direction(df, hess_f, x, method)
pt.show()
print("method {2} converged to: {0}, in {1} iterations".format(x, len(guesses), method))
It works OK: For example, for the function
We get
method 0 converged to: [ 1.37484167e-04 -8.24905001e-06], in 22 iterations
method 1 converged to: [ 0. 0.], in 3 iterations
The only thing that really bothers me if the fact that my optimization functions explicitly requires not only the function I'm minimizing (which of course makes sense), but also the gradient and the Hessian. So right now I'm essentially "hard coding" these functions, like so:
def f2(x):
return 100*(x[0]-3)**2 + (x[1]-1)**2
def df2(x):
return np.array([200*(x[0]-3), 2*(x[1]-1)])
def hess_f2(x):
return np.array([200, 0, 0, 2]).reshape((2,2))
So my question is, what is a "Pythonic" way to generate functions that calculate the gradient and Hessian of the input function function, in a way that would fit my above implementation? I'm guessing it's something very simple but my main experience with Python is scripting, so I haven't yet done something like this. Thanks!
