""" Gradient descent ================== An example demoing gradient descent by creating figures that trace the evolution of the optimizer. """ import numpy as np import pylab as pl from scipy import optimize from cost_functions import mk_quad, mk_gauss, rosenbrock,\ rosenbrock_prime, rosenbrock_hessian, LoggingFunction,\ CountingFunction x_min, x_max = -1, 2 y_min, y_max = 2.25/3*x_min - .2, 2.25/3*x_max - .2 ############################################################################### # A formatter to print values on contours def super_fmt(value): if value > 1: if np.abs(int(value) - value) < .1: out = '$10^{%.1i}$' % value else: out = '$10^{%.1f}$' % value else: value = np.exp(value - .01) if value > .1: out = '%1.1f' % value elif value > .01: out = '%.2f' % value else: out = '%.2e' % value return out ############################################################################### # A gradient descent algorithm # do not use: its a toy, use scipy's optimize.fmin_cg def gradient_descent(x0, f, f_prime, hessian=None, adaptative=False): x_i, y_i = x0 all_x_i = list() all_y_i = list() all_f_i = list() for i in range(1, 100): all_x_i.append(x_i) all_y_i.append(y_i) all_f_i.append(f([x_i, y_i])) dx_i, dy_i = f_prime(np.asarray([x_i, y_i])) if adaptative: # Compute a step size using a line_search to satisfy the Wolf # conditions step = optimize.line_search(f, f_prime, np.r_[x_i, y_i], -np.r_[dx_i, dy_i], np.r_[dx_i, dy_i], c2=.05) step = step[0] if step is None: step = 0 else: step = 1 x_i += - step*dx_i y_i += - step*dy_i if np.abs(all_f_i[-1]) < 1e-16: break return all_x_i, all_y_i, all_f_i def gradient_descent_adaptative(x0, f, f_prime, hessian=None): return gradient_descent(x0, f, f_prime, adaptative=True) def conjugate_gradient(x0, f, f_prime, hessian=None): all_x_i = [x0[0]] all_y_i = [x0[1]] all_f_i = [f(x0)] def store(X): x, y = X all_x_i.append(x) all_y_i.append(y) all_f_i.append(f(X)) optimize.fmin_cg(f, x0, f_prime, callback=store, gtol=1e-12) return all_x_i, all_y_i, all_f_i def newton_cg(x0, f, f_prime, hessian): all_x_i = [x0[0]] all_y_i = [x0[1]] all_f_i = [f(x0)] def store(X): x, y = X all_x_i.append(x) all_y_i.append(y) all_f_i.append(f(X)) optimize.fmin_ncg(f, x0, f_prime, fhess=hessian, callback=store, avextol=1e-12) return all_x_i, all_y_i, all_f_i def bfgs(x0, f, f_prime, hessian=None): all_x_i = [x0[0]] all_y_i = [x0[1]] all_f_i = [f(x0)] def store(X): x, y = X all_x_i.append(x) all_y_i.append(y) all_f_i.append(f(X)) optimize.fmin_bfgs(f, x0, f_prime, callback=store, gtol=1e-12) return all_x_i, all_y_i, all_f_i def powell(x0, f, f_prime, hessian=None): all_x_i = [x0[0]] all_y_i = [x0[1]] all_f_i = [f(x0)] def store(X): x, y = X all_x_i.append(x) all_y_i.append(y) all_f_i.append(f(X)) optimize.fmin_powell(f, x0, callback=store, ftol=1e-12) return all_x_i, all_y_i, all_f_i def nelder_mead(x0, f, f_prime, hessian=None): all_x_i = [x0[0]] all_y_i = [x0[1]] all_f_i = [f(x0)] def store(X): x, y = X all_x_i.append(x) all_y_i.append(y) all_f_i.append(f(X)) optimize.fmin(f, x0, callback=store, ftol=1e-12) return all_x_i, all_y_i, all_f_i ############################################################################### # Run different optimizers on these problems levels = dict() for index, ((f, f_prime, hessian), optimizer) in enumerate(( (mk_quad(.7), gradient_descent), (mk_quad(.7), gradient_descent_adaptative), (mk_quad(.02), gradient_descent), (mk_quad(.02), gradient_descent_adaptative), (mk_gauss(.02), gradient_descent_adaptative), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), gradient_descent_adaptative), (mk_gauss(.02), conjugate_gradient), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), conjugate_gradient), (mk_quad(.02), newton_cg), (mk_gauss(.02), newton_cg), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), newton_cg), (mk_quad(.02), bfgs), (mk_gauss(.02), bfgs), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), bfgs), (mk_quad(.02), powell), (mk_gauss(.02), powell), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), powell), (mk_gauss(.02), nelder_mead), ((rosenbrock, rosenbrock_prime, rosenbrock_hessian), nelder_mead), )): # Compute a gradient-descent x_i, y_i = 1.6, 1.1 counting_f_prime = CountingFunction(f_prime) counting_hessian = CountingFunction(hessian) logging_f = LoggingFunction(f, counter=counting_f_prime.counter) all_x_i, all_y_i, all_f_i = optimizer(np.array([x_i, y_i]), logging_f, counting_f_prime, hessian=counting_hessian) # Plot the contour plot if not max(all_y_i) < y_max: x_min *= 1.2 x_max *= 1.2 y_min *= 1.2 y_max *= 1.2 x, y = np.mgrid[x_min:x_max:100j, y_min:y_max:100j] x = x.T y = y.T pl.figure(index, figsize=(3, 2.5)) pl.clf() pl.axes([0, 0, 1, 1]) X = np.concatenate((x[np.newaxis, ...], y[np.newaxis, ...]), axis=0) z = np.apply_along_axis(f, 0, X) log_z = np.log(z + .01) pl.imshow(log_z, extent=[x_min, x_max, y_min, y_max], cmap=pl.cm.gray_r, origin='lower', vmax=log_z.min() + 1.5*log_z.ptp()) contours = pl.contour(log_z, levels=levels.get(f, None), extent=[x_min, x_max, y_min, y_max], cmap=pl.cm.gnuplot, origin='lower') levels[f] = contours.levels pl.clabel(contours, inline=1, fmt=super_fmt, fontsize=14) pl.plot(all_x_i, all_y_i, 'b-', linewidth=2) pl.plot(all_x_i, all_y_i, 'k+') pl.plot(logging_f.all_x_i, logging_f.all_y_i, 'k.', markersize=2) pl.plot([0], [0], 'rx', markersize=12) pl.xticks(()) pl.yticks(()) pl.xlim(x_min, x_max) pl.ylim(y_min, y_max) pl.draw() pl.figure(index + 100, figsize=(4, 3)) pl.clf() pl.semilogy(np.maximum(np.abs(all_f_i), 1e-30), linewidth=2, label='# iterations') pl.ylabel('Error on f(x)') pl.semilogy(logging_f.counts, np.maximum(np.abs(logging_f.all_f_i), 1e-30), linewidth=2, color='g', label='# function calls') pl.legend(loc='upper right', frameon=True, prop=dict(size=11), borderaxespad=0, handlelength=1.5, handletextpad=.5) pl.tight_layout() pl.draw()