function [f,gradf]=objfun(x)
f=(x(1)^2+x(2)^2)^2-x(1)^2-x(2)+x(3)^2;
gradf=[4*x(1)*(x(1)^2+x(2)^2)-2*x(1);4*x(2)*(x(1)^2+x(2)^2)-1;2*x(3)];
options=optimset('GradObj','on');
x0=[1;1;1];[x,fval]=fminunc('objfun',x0,options)
Python code for the "red code"
{{{id=6| import numpy from scipy.optimize import fmin_bfgs def f(x): return (x[0]**2+x[1]**2)**2-x[0]**2-x[1]+x[2]**2 def fprime(x): grad = numpy.zeros((len(x),), float) grad[0] = 4*x[0]*(x[0]**2+x[1]**2)-2*x[0] grad[1] = 4*x[1]*(x[0]**2+x[1]**2)-1 grad[2] = 2*x[2] return grad x0=[1,1,1] print f(x0), fprime(x0) #fmin with gradient provided time xopt = fmin_bfgs(f,x0,fprime) print xopt #fmin with gradient estimated by fmin module time xopt1 = fmin_bfgs(f,x0) print xopt1 /// 3 [ 6. 7. 2.] Optimization terminated successfully. Current function value: -0.500000 Iterations: 11 Function evaluations: 15 Gradient evaluations: 15 Time: CPU 0.01 s, Wall: 0.01 s [ -4.99999944e-01 5.00000405e-01 -3.84579279e-07] Optimization terminated successfully. Current function value: -0.500000 Iterations: 11 Function evaluations: 75 Gradient evaluations: 15 Time: CPU 0.01 s, Wall: 0.01 s [ -4.99999979e-01 5.00000380e-01 -3.92121138e-07] }}}Finds the extreme point analytically and computes the eigenvalues of the Hessian (second derivative in one dimension) to determine
if the function has a minimum or maximum at this point.
{{{id=7| f(x,y)=x^2*y+y^2+y f.diff() # gradient res=solve(list(f.diff()),[x,y]) print res H=f.diff(2) # Hessian matrix print H H1=H(x=0,y=-1/2) print H1 H(x=0,y=-1/2).eigenvalues() /// [ [x == -I, y == 0], [x == I, y == 0], [x == 0, y == (-1/2)] ] [(x, y) |--> 2*y (x, y) |--> 2*x] [(x, y) |--> 2*x (x, y) |--> 2] [-1 0] [ 0 2] [-1, 2] }}}Sage implementation of the "red code"
1) it computes analytically the extreme point and 2) numerically using the method of BFGS which is a modified Newton
{{{id=8| sage: var('x , y, z') sage: f(x,y,z)=(x^2+y^2)^2-x^2-y+z^2 sage: f.diff() # gradient sage: res=solve(list(f.diff()),[x,y,z]) sage: print 'roots of gradient', res sage: x0=[1, 1, 1] sage: xopt=minimize(f,x0,algorithm="bfgs") sage: print 'values f at local minimum', f(xopt[0],xopt[1],xopt[2]) sage: H=f.diff(2) # Hessian matrix sage: print 'Hessian', H sage: print 'value of Hessian at xopt', H(xopt[0],xopt[1],xopt[2]) sage: H(xopt[0],xopt[1],xopt[2]).eigenvalues() /// roots of gradient [ [x == (-2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0], [x == (2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0], [x == (2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0], [x == (-2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0], [x == 0.0, y == 0.629960578187, z == 0.0], [x == (-1/2), y == (1/2), z == 0], [x == (1/2), y == (1/2), z == 0] ] Optimization terminated successfully. Current function value: -0.500000 Iterations: 11 Function evaluations: 15 Gradient evaluations: 15 values f at local minimum -0.5 Hessian [(x, y, z) |--> 12*x^2 + 4*y^2 - 2 (x, y, z) |--> 8*x*y (x, y, z) |--> 0] [ (x, y, z) |--> 8*x*y (x, y, z) |--> 4*x^2 + 12*y^2 (x, y, z) |--> 0] [ (x, y, z) |--> 0 (x, y, z) |--> 0 (x, y, z) |--> 2] value of Hessian at xopt [ 2.00000094466 -2.00000139411 0] [-2.00000139411 4.00000463179 0] [ 0 0 2] [-13/1383966880035651535063*sqrt(56667688569868626270919004239662009317529) + 345565067138652/115188248656069, 13/1383966880035651535063*sqrt(56667688569868626270919004239662009317529) + 345565067138652/115188248656069, 2] }}}3d plotting example
{{{id=10| var('x,y') p3d = plot3d(x^2-y^2, (x, -10, 10), (y, -10, 10), mesh=true,frame=true) p3d.show() /// }}} {{{id=11| help(minimize) ///|
Click to open help window
|
Implementation of the gradient method for one variable objective function
{{{id=12| # From calculation, we expect that the local minimum occurs at x=9/4 x_old = 0 x_new = 6 # The algorithm starts at x=6 eps = 0.01 # step size precision = 0.00001 def f_prime(x): return 4 * x**3 - 9 * x**2 while abs(x_new - x_old) > precision: x_old = x_new x_new = x_old - eps * f_prime(x_old) print "Local minimum occurs at ", x_new /// Local minimum occurs at 2.24996460742785 }}}vector fields of a parametric function
{{{id=13| var('x,y') plot_vector_field((sin(x),cos(y)),(x,-pi,pi),(y,-pi,pi)) ///
}}}
{{{id=14|
f = e^x + e^y*sin(x)
f.gradient()
///
(e^y*cos(x) + e^x, e^y*sin(x))
}}}
{{{id=15|
g = x^2 + x*y
plot_vector_field(g.gradient(),(x,-2,2),(y,-2,2))
///
}}}
{{{id=16|
var('x,y,z')
G = (x+y,x*y,z)
plot_vector_field3d(G,(x,0,1),(y,0,1),(z,0,1))
///
}}}
{{{id=17|
h = x*y*z
plot_vector_field3d(h.gradient(),(x,-2,2),(y,-2,2),(z,-2,2))
///
}}}
{{{id=18|
var('x,y,z,t')
H = (z,y,x)
r = (cos(6*t),sin(6*t),t)
vfield = plot_vector_field3d(H,(x,-1,1),(y,-1,1),(z,0,pi),colors='orange')
curve = parametric_plot3d(r,(t,0,pi),color='red')
vfield + curve
///
}}}
{{{id=19|
var('theta,phi')
parametric_plot3d((2*sin(phi)*cos(theta),2*sin(phi)*sin(theta),\
2*cos(phi)),(theta,0,2*pi),(phi,0,pi))
///
}}}
{{{id=20|
var('s,t')
parametric_plot3d((s,t,s^2+t^2),(s,-1,1),(t,-1,1))
///
}}}
Gradient method for two variables objective function
Notice that plots the sequence of extreme points on the contour plot of the function
{{{id=22| var('x,y') import numpy as np xold=vector([1,1]) xnew=vector([0.5, 0.5]) precision=0.00001 pt=[] # defines the list of extreme points it computes in each iteration eps = 0.01 # step size f(x,y)=x^2+y^2 fprime=f.gradient() print(f) print(fprime) print "analytic solution of fprime",solve(list(fprime),[x,y]) temp=xnew-xold nnorm=norm(temp) plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) #the implementation of the gradient algorithm iter=1 while nnorm > precision: xold=xnew xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1])) nnorm=norm(xnew-xold) pt.append(xnew) # sets the list of the sequence of extreme points it computes iter=iter+1 print "local minimum at ",xnew print "number of iterations:",iter lpt=point(pt) show(plt+lpt) /// (x, y) |--> x^2 + y^2 (x, y) |--> (2*x, 2*y) analytic solution of fprime [ [x == 0, y == 0] ] local minimum at (0.000340081752004604, 0.000340081752004604) number of iterations: 362
}}}
Computes the exteme points of f(x,y,z) using the method BFGS (a modified Newton)
{{{id=24| var('x , y, z') f(x,y,z)=(x^2+y^2)^2-x^2-y+z^2 fprime=f.diff() # gradient print(fprime) res=solve(list(f.diff()),[x,y,z]) print res x0=[1, 1, 1] minimize(f,x0,algorithm="bfgs") /// (x, y, z) |--> (4*(x^2 + y^2)*x - 2*x, 4*(x^2 + y^2)*y - 1, 2*z) [ [x == (-2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0], [x == (2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0], [x == (2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0], [x == (-2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0], [x == 0.0, y == 0.629960578187, z == 0.0], [x == (-1/2), y == (1/2), z == 0], [x == (1/2), y == (1/2), z == 0] ] Optimization terminated successfully. Current function value: -0.500000 Iterations: 11 Function evaluations: 15 Gradient evaluations: 15 (-0.499999943818, 0.50000040471, -3.84579279067e-07) }}}Gradient method for f(x,y,z)
{{{id=26| var('x,y,z') xold=vector([0.5, 0.5, 0.5]) xnew=vector([0.2, 0.2, 0.2]) precision=0.00000001 eps = 0.001 # step size f(x,y,z)=x^2+y^2+z^2 fprime=f.gradient() print(f) print(fprime) print "analytic solution of fprime",solve(list(fprime),[x,y,z]) temp=xnew-xold nnorm=norm(temp) while nnorm > precision: xold=xnew xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1],xold[2])) nnorm=norm(xnew-xold) print "local minimum at ",xnew /// (x, y, z) |--> x^2 + y^2 + z^2 (x, y, z) |--> (2*x, 2*y, 2*z) analytic solution of fprime [ [x == 0, y == 0, z == 0] ] local minimum at (2.87749398873009e-6, 2.87749398873009e-6, 2.87749398873009e-6) }}}Computes and plots the extreme point on the contour plot of the function
{{{id=60| /// }}} {{{id=29| var('x, y') f = 100 * (y - x^2)^2 + (1 - x)^2 + 100 * (2 - y^2)^2 + (1 - y)^2 min = minimize(f, [-.1,0.5], disp=1) print "extreme point",min plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) pt = point(min, color='red', size=50) show(plt+pt, aspect_ratio=1, figsize=(4, 4)) pf = plot3d(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) show(pf) /// Optimization terminated successfully. Current function value: 4.953484 Iterations: 8 Function evaluations: 17 Gradient evaluations: 17 extreme point (-1.184618747, 1.41254232561)
}}}
{{{id=28|
help(minimize)
///
|
Click to open help window
|
Displays matrices graphically
{{{id=27| sage: pi = float(pi) sage: m = matrix(RDF, 6, [sin(i^2 + j^2) for i in [0,pi/5,..,pi] for j in [0,pi/5,..,pi]]) sage: list_plot3d(m, texture='yellow', frame_aspect_ratio=[1,1,1/3]) /// }}} {{{id=30| sage: l=[] sage: for i in range(6): ... for j in range(6): ... l.append((float(i*pi/5),float(j*pi/5),m[i,j])) sage: list_plot3d(l,texture='yellow') /// }}}Displays lists graphically
{{{id=31| sage: l=[] sage: for i in range(-5,5): ... for j in range(-5,5): ... l.append((normalvariate(0,1),normalvariate(0,1),normalvariate(0,1))) sage: list_plot3d(l,interpolation_type='nn',texture='yellow',num_points=100) /// }}}Applies gradient method and displays the sequence of approximations of the extreme point in the contour graph of
the objective function
{{{id=32| var('x,y') import numpy as np xold=vector([-.25,.25]) xnew=vector([-0.1, 0.5]) precision=0.00001 pt=[] eps = 0.1 # step size f(x,y)=sin(x^2/2-y^2/4+3)*cos(2*x+1-e^y) fprime=f.gradient() print(f) print(fprime) amin=solve(list(fprime),[x,y]) print "analytic solution of fprime",amin temp=xnew-xold nnorm=norm(temp) plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) iter=1 while nnorm > precision: xold=xnew xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1])) nnorm=norm(xnew-xold) pt.append(xnew) iter=iter+1 print "local minimum at ",xnew print "number of iterations:",iter #print pt lpt=point(pt,color='red') show(plt+lpt) pf = plot3d(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) show(pf) /// (x, y) |--> sin(1/2*x^2 - 1/4*y^2 + 3)*cos(2*x - e^y + 1) (x, y) |--> (x*cos(2*x - e^y + 1)*cos(1/2*x^2 - 1/4*y^2 + 3) - 2*sin(2*x - e^y + 1)*sin(1/2*x^2 - 1/4*y^2 + 3), -1/2*y*cos(2*x - e^y + 1)*cos(1/2*x^2 - 1/4*y^2 + 3) + e^y*sin(2*x - e^y + 1)*sin(1/2*x^2 - 1/4*y^2 + 3)) analytic solution of fprime [ x*cos(-2*x + e^y - 1)*cos(-1/2*x^2 + 1/4*y^2 - 3) - 2*sin(-2*x + e^y - 1)*sin(-1/2*x^2 + 1/4*y^2 - 3), -1/2*y*cos(-2*x + e^y - 1)*cos(-1/2*x^2 + 1/4*y^2 - 3) + e^y*sin(-2*x + e^y - 1)*sin(-1/2*x^2 + 1/4*y^2 - 3) ] local minimum at (0.322615850969378, 1.60235597131258) number of iterations: 145
}}}
Application of the gradient method for the obective function $(x-y)/(x^2+y^2+1)$
Notice that it displays the sequence of the extreme point approximations on the contour graph of the function
{{{id=35| var('x,y') import numpy as np xold=vector([0,0]) xnew=vector([-0.1, 0.5]) precision=0.00001 pt=[] eps = 0.1 # step size f(x,y)=(x-y)/(x^2+y^2+1) # function to find maximun fprime=f.gradient() print(f) print(fprime) amin=solve(list(fprime),[x,y]) print "analytic solution of fprime",amin temp=xnew-xold nnorm=norm(temp) plt = contour_plot(f, (x, -3, 3), (y, -3, 3),fill=False,cmap='hsv',labels=True) iter=1 while nnorm > precision: xold=xnew xnew=vector(xold)+eps*vector(fprime(xold[0],xold[1])) nnorm=norm(xnew-xold) pt.append(xnew) iter=iter+1 print "local minimum at ",xnew print "number of iterations:",iter #print pt lpt=point(pt,color='red') show(plt+lpt) pf = plot3d(f, (x, -3, 3), (y, -3, 3),fill=False,cmap='hsv',labels=True) show(pf) /// (x, y) |--> (x - y)/(x^2 + y^2 + 1) (x, y) |--> (-2*(x - y)*x/(x^2 + y^2 + 1)^2 + 1/(x^2 + y^2 + 1), -2*(x - y)*y/(x^2 + y^2 + 1)^2 - 1/(x^2 + y^2 + 1)) analytic solution of fprime [ [x == -1/2*sqrt(2), y == 1/2*sqrt(2)], [x == 1/2*sqrt(2), y == -1/2*sqrt(2)], [x == -1/2*I*sqrt(2), y == -1/2*I*sqrt(2)], [x == 1/2*I*sqrt(2), y == 1/2*I*sqrt(2)] ] local minimum at (0.707100597445379, -0.706979236153001) number of iterations: 117
}}}
{{{id=55|
y=function('y',x)
var('x,yy,z')
f=x*y+x^2-y^3*x
ff=f.diff(x).subs({y.diff(x):z,y:yy})
implicit_plot3d(ff,(x,-3,3),(yy,-3,3),(z,-4,4))
///
}}}
{{{id=56|
var('x,y')
f(x,y)=4 * x**3 - 9 * x**2-y
f.gradient()
#lt=plot3d(f,(x,0,2),(y,0,2))
#show(lt)
#plt = contour_plot(f, (x,0,2),(y,0,2),fill=False,cmap='hsv',labels=True)
#show(plt)
fp=plot_vector_field(f.gradient(),(x,0,2),(y,0,4))
fpl=plot(4*x**3-9*x**2,(x,0,2))
show(fpl+fp)
///
}}}