Class Notes for Ch 4, Nonlinear Programming

CLASS NOTES for CHAPTER 4, Nonlinear Programming 1 Space Curves and Tangent Lines Recall that space curves are defined by a set of parametric equations,   x1 (t)  x2 (t)    r(t) =  ..   .  xn (t) In Calc III, we might have written this a little differently, ~r(t) =< x(t), y(t), z(t) > but here we want to use n dimensions rather than two or three dimensions. The derivative and antiderivative of r with respect to t is done componentwise,    R  x01 (t) x (t) dt 1  x0 (t)   R x2 (t) dt    2   r(t) =  ..  , R(t) =   ..  .    . R 0 xn (t) xn (t) dt And the local linear approximation to r(t) is also done componentwise. The tangent line (in n dimensions) can be written easily- the derivative at t = a is the direction of the curve, so the tangent line is given by:     x01 (a) x1 (a)  x0 (a)   x2 (a)    2   L(t) =  ..  + t  ..   .   .  x0n (a) xn (a) In Class Exercise: Use Maple to plot the curve r(t) = [cos(t), sin(t), t] T and its tangent line at t = π2 . 2 Gradients and Tangent Planes Let f : Rn → R. In this case, we can write: y = f (x1 , x2 , x3 , . . . , xn ) 1 Note that any function that we wish to optimize must be of this form- It would not make sense to find the maximum of a function like a space curve; n dimensional coordinates are not well-ordered like the real line- so the following statement would be meaningless: (3, 5) > (1, 2). This type of function is sometimes called a surface even though we can’t graph it if the domain has a higher dimension than 2. The gradient of f , ∇f , is defined by its partial derivatives. We will denote ∂f as fxi , and so: ∂xi    ∇f (x1 , x2 , . . . , xn ) =   fx1 (x1 , x2 , . . . , xn ) fx2 (x1 , x2 , . . . , xn ) .. . fxn (x1 , x2 , . . . , xn )      NOTE: When writing the gradient, we will not distinguish between it being a row vector or a column vector. For some ideas, it may be more clear notationally to write it as a column, but for computations, it is a row vector. Normally, it will be clear from the context whether it is a row or column. The rate of change of f at x = a in the direction of the unit vector u is defined as the directional derivative of f : Du f (a) = ∇f (a) · u If we recall the useful relationship between any two vectors u, v ∈ Rn , u · v = kuk2 kvk2 cos(θ) then we see that the direction for which the rate of change of f is a maximum is in the direction of ∇f (a). In Class Example: Find the maximum rate of change of f at the given point and the direction in which it occurs: f (x, y) = x + y/z at (4, 3). The linearization of f at x = a is the tangent plane, L(x) = f (a) + ∇f (a) · (x − a) which is also the first order Taylor series approximation to f at a. 2 2.1 Surfaces and Contours In maximizing or minimizing a function, we’ll consider equations like: F (x) = k, where k is a real number. From our earlier discussion on the Implicit Function Theorem, we know that this equation could implicitly define one of the variables in terms of the others. In Calculus III, for example, we saw • F (x, y) = k defines a level curve in the plane. • F (x, y, z) = k defines a level surface in three dimensions. Now let r(t) be a curve on the surface defined by F (x) = k. Then we could write F as dependent on t, and differentiate with respect to t: ∇F (r(t)) · r0 (t) = 0 which implies that the gradient is orthogonal to any curve on the surface. In particular, we can compute the tangent line or plane at the point (a, b) or (a, b, c): • Fx (a, b) · (x − a) + Fy (a, b) · (y − b) = 0 • Fx (a, b, c) · (x − a) + Fy (a, b, c) · (y − b) + Fz (a, b, c) · (z − c) = 0 If x ∈ Rn , this can be written compactly as ∇F (a) · (x − a) = 0 In Class Example: Compute the tangent line ´ and show that the gradient ³ is orthogonal to it, if x2 + y 2 = 1 at √12 , √12 . In this case, the gradient of √ √ F (x, y) = x2 + y 2 is [2x, 2y], and at the given point, the gradient is [ 2, 2]. Notice that this vector points directly opposite from the center of the circle. On the other hand, the equation of the tangent line is: µ ¶ µ ¶ √ √ 1 1 + 2 y−√ =0 2 x− √ 2 2 √ or in the more familiar form, y = −x + 2. The “slope” of the gradient is 1, while the slope of the tangent line is −1. In two dimensions such as this, a Fy line in the direction of the gradient would have slope , while the slope of Fx Fx the tangent line would be − . Fy The important point in this section: Given a level surface of the form F (x) = k, the gradient ∇F at x = a is orthogonal. Furthermore, the gradient points in the direction where F increases the fastest (for y = F (x). 3 3 3.1 Derivatives we don’t see in Calculus III The Hessian Matrix Let f : Rn → R, so that we have something in the form of y = f (x1 , x2 , . . . , xn ). We have already seen that the derivative of f takes the form of the gradient. We now examine the second derivative of f , and we’ll see that the second derivative forms what is called the Hessian matrix. The first derivative is in the form of the gradient, where we have a vector of functions. Let fxi (x1 , x2 , . . . , xn ) = gi (x1 , x2 , . . . , xn ). Now each function g is also a function from Rn to R, and so the derivative of g is also a gradient. We could therefore write the second derivative of f in matrix form to obtain what is called the Hessian of f :     f x1 x1 f x1 x2 . . . f x1 xn ∇g1  ∇g2   fx x fx x . . . fx x  2 2 2 n    2 1  Hf =  ..  =  .. ..  ..  .   . .  . f xn x1 f xn x2 . . . f xn xn ∇gn In Class Example: Let f (x, y) = x3 + 2xy − y 2 . Compute the gradient and Hessian of f . · 2 ¸ · ¸ 3x + 2y 6x + 2 2 ∇f = , Hf = 2x − 2y 2 −2 You might note that, if the second partials of f are continuous, then Hf will always be a symmetric matrix (this is a consequence of Clairaut’s Theorem). In Class Example: Let f (x, y) = 2 + x2 + x + xy + y + y 2 . Show that f can be written as: 1 f (a) + ∇f (a)(x − a) + (x − a)T Hf (a)(x − a) 2 where a = (0, 0). First, we compute ∇f and Hf : ¸ · 2x + 1 + y , ∇f = x + 1 + 2y 4 Hf = · 2 1 1 2 ¸ Substituting (x, y) = (0, 0), the gradient is [1, 1], and now compute: ¸ ¸ ¸· · · 1 x−0 x−0 2 1 + [x − 0, y − 0] 2 + [1, 1] · y−0 y−0 1 2 2 · ¸ 1 2x + y 2 + x + y + [x, y] = 2 + x + y + x2 + xy + y 2 x + 2y 2 The Second Order Taylor Approximation If f : Rn → R, then the Taylor approximation to f at a to second order is: 1 f (x) ≈ f (a) + ∇f (a)(x − a) + (x − a)T Hf (a)(x − a) 2 3.2 Connecting the Hessian to Calculus III In Calc III, when we were finding local extrema, we had the “Second Derivatives Test”. To generalize this result requires a little linear algebra, but the end result is that if x = a is a critical point of f , and Hf (a) has all positive eigenvalues, then f has a local minimum at a. If Hf (a) has all negative eigenvalues, then f has a local maximum at a. If the eigenvalues are mixed in sign, there is a saddle point. 3.3 The Jacobian Matrix Let f : Rn → Rm . The we could write f using m coordinate functions, each depending on n variables,   f1 (x1 , . . . , xn )  f2 (x1 , . . . , xn )    f (x1 , . . . , xn ) =   ..   . fm (x1 , . . . , xn ) Definition: The derivative of f takes the form of a matrix called the Jacobian of f , constructed by taking the partial derivatives:   (f1 )x1 (f1 )x2 · · · (f1 )xn  (f2 )x (f2 )x · · · (f2 )x  n  1 2  Jf =   .. .. ..   . . . (fm )x1 (fm )x2 · · · (fm )xn 5 In Class Example: Compute the Jacobian of:  2  x + yz + z f (x, y, z) =  y 2 + xy + 3  sin(xyz)   2x z y+1  y 2y + x 0 Jf =  yz cos(xyz) xz cos(xyz) xy cos(xyz) If f : Rn → Rm , then the linearization of f at x = a is given by: L(x) = f (a) + Jf (a)(x − a) 4 The Extreme Value Theorem Let f : Rn → R be continuous on a closed and bounded domain D. Then f attains a global maximum and global minimum on D. Furthermore, this maximum and minimum value is obtained either at the critical points of f (where the gradient is zero or not defined), or on the boundary of D. In Class Example: Find the maximum value of f (x, y) = x2 − 2xy + 2y where the domain is the rectangle D = {(x, y)|x ∈ [0, 3], y ∈ [0, 2]} The gradient is [2x − 2y, −2x + 2]. Setting these to zero, we see the only solution is (1, 1). The value of the function here is f (1, 1) = 1. Since the boundary is made up of line segments, it is easy to analyze f along them: • For y = 0, 0 ≤ x ≤ 3, f (x, 0) = x2 , and the max of 9 occurs at x = 3 • For x = 0, 0 ≤ y ≤ 2, f (0, y) = 2y, and the max of 4 occurs at y = 2. • For x = 3, 0 ≤ y ≤ 2, f (3, y) = 9 − 4y, and the max of 5 occurs at y = 0. • For y = 2, 0 ≤ x ≤ 3, f (x, 2) = x2 − 4x + 4 = (x − 2)2 , and has a max of 4 at x = 0. Now we see that the global max is 9 and occurs at the boundary point (3, 0). 6 5 Newton’s Method Revisited We saw that for a problem of the form f (x) = 0, we could solve for x iteratively by using an initial guess, x0 , and iterate: xk+1 = xk − f (xk ) f 0 (xk ) We also saw that there could be problems with Newton’s Method: • If x∗ is the desired point, then f 0 (x∗ ) 6= 0. • We’re not sure which root we’re converging to. • We can only guarantee convergence if x0 is sufficiently close to x∗ . We got the formula for Newton’s Method by linearizing the function f at xk , then we set y = 0 and solved for xk+1 : y − f (xk ) = f 0 (xk )(xk+1 − xk ) ⇒ −f (xk ) = f 0 (xk )(xk+1 − xk ) 5.1 Scalar Valued Functions If f : Rn → R, then Newton’s Method is not well defined. For example, suppose f (x, y) = x2 + y 2 − 1. Then the zeros to this function are all points on the unit circle- Generally, the solution set is a level surface (or level curve). 5.2 Vector Valued Functions Let g : Rn → Rn . The derivative of g is then its Jacobian matrix (n × n), and the linearization of g about xk is given by: y − g(xk ) = Jg(xk )(xk+1 − xk ) Set y = 0 and solve for xk+1 , and we get: xk+1 = xk − (Jg(xk ))−1 g(xk ) Let’s do a simple example just to see what is at work here. Let g(x, y) = [x2 − y 2 , 2x + y − 1]. Then the Jacobian matrix at a point (x, y) is given by: · ¸ 2x −2y Jg(x, y) = 2 1 7 And since this is a 2 × 2 matrix, its inverse can be computed directly: ¸ · 1 1 2y −1 (Jg) (x, y) = 2x + 4y −2 2x Let the initial guess be (1, 1). The next iterate is then: ¸ ¸· ¸ · · ¸ · 1 1/3 0 1 1 2 = − x1 = 1/3 2 1 6 −2 2 You can double check that this is indeed a solution to g(x) = 0, but this was because we had a lucky first guess. Other initial conditions may take longer for the algorithm to converge. 5.3 Newton’s Method in Optimization In one dimension, we look for points where f 0 (x) = 0. Numerically, we can use Newton’s Method to find the x by letting g(x) = f 0 (x): xk+1 = xk − g(xk ) f 0 (xk ) = x − k g 0 (xk ) f 00 (xk ) If f : Rn → R, then we can think of ∇f : Rn → Rn . To find the zeros of the gradient, we can use Newton’s Method in a similar way. Let g(x) = ∇f . Then, from our previous computation, xk+1 = xk − (Jg(xk ))−1 g(xk ) = xk − (Hf (xk ))−1 ∇f (xk ) Note that this last product is an n × n inverse matrix of the Hessian times the gradient as a column vector. If the domain of f is unbounded, we may or may not have an optimal value- recall that we are simply finding where the gradient is zero, which may correspond to a local or global maximum, local or global minimum, or a saddle point. EXAMPLE: Find the maximum of f (x, y) = 12 + 10y − 2x2 − 8xy − y 4 . Since this is a function of two variables, we can look at the contour plot to predict what will happen, and this is given in Figure 1. After looking at the function, it is clear that there is a global maximum, since the value of the function will tend to −∞ as (x, y) become unbounded. 8 3 2 y 1 –10 –8 –4 –6 –2 0 2 4 6 x –1 –2 Figure 1: Contours of f (x, y) = 12 + 10y − 2x2 − 8xy − y 4 . We see that the contours are encircling two points and there might be a saddle in between. We compute the gradient and Hessian of f : £ ¤ ∇f = −4x − 8y 10 − 8x − 4y 3 · ¸ −4 −8 Hf = −8 −12y 2 We can look at the gradient directly- Solutions to where the gradient is zero will lie along the line y = − 21 x, from which we would have to solve: 10 − 8x + 21 x3 = 0, and we get three solutions for x, which are approximately: x ≈ −4.159, x ≈ 3.084, x ≈ 1.434 Using the Matlab script below, we can numerically solve for these using Newton’s Method: 9 x=[-8;1]; for j=1:30 gradf=[-4*x(1)-8*x(2); 10-8*x(1)-4*x(2)^3]; Hf=[-4 8; -8 -12*x(2)^2]; iHf=inv(Hf); x=x-iHf*gradf; end Running this code will give a final value of approximately (−4.519, 2.2597). We can put additional things into the code so that we can track the values of x, ∇f and Hf . We could also use a stopping criteria rather than running the algorithm 30 times (30 was chosen arbitrarily). At this point, the functional value was approximately 49.37. Starting Newton’s Method at (3, −1) gave a different ending value of approximately (3.084, −1.542), and the value of f was 9.484. Starting Newton’s Method at (1, −1) gave a different ending value of approximately (3.084, −1.542). Starting Newton’s Method at (1.5, −0.5) gave us (−4.519, 2.259). Starting very close to the third point we found by other methods, Newton’s method did not want to converge to it- it is in fact a saddle point. We see some typical behavior in this example: • If Newton’s method starts close to the optimum, it will reach that point fairly quickly. • The solution given by Newton’s method may be only a local maximum (or minimum). • Different initial conditions may give different answers. In practice, it is common to run an optimization algorithm with many randomly chosen initial values to probe the surface. 6 Lagrange Multipliers In the last section, we considered optimizing a nonlinear function over an unbounded domain. In this section, we place constraints on the domain- in particular, we assume the constraints are of a particular form: g(x) = k. 10 6 y 4 2 –2 2 4 6 x –2 Figure 2: Some contours associated with f (x, y) = x3 + y 3 + 3xy (these are 0, 3.5, 20, 100, 200, 350, 400) which are shown as dotted curves (the higher the value, the higher the curve for this function). The constraint (x − 3)2 + (y − 3)2 = 9 is shown as the solid curve. Before looking at the case in Rn , let’s consider the problem where we are optimizing a function f (x, y) under constraint that g(x, y) = k. In this case, plots of the form f (x, y) = c are level curves of f , and the plot of g(x, y) = k is a curve. For example, consider f (x, y) = x3 + y 3 − 3xy subject to the restriction that the points lie on a circle of radius 3 centered at (3, 3). In Figure 2, we show a plot of some contours of f together with the circle. We’ll recall that the gradient of f , in this case ∇f (x, y) = [3x2 + 3y, 3y 2 + 3x] is orthogonal to the level curve at a given point on the level curve. Similarly, ∇g(x, y) = [2(x − 3), 11 2(y − 3)] is orthogonal to the circle. Now, consider f restricted to the curve where we’ll assume that g(x, y) = k has been parametrized by ~r(t) = [x(t), y(t)]. In the example, this would mean x(t) = 3 cos(t) + 3, y(t) = 3 sin(t) + 3 Note that: ~r0 (t) = [−3 sin(t), 3 cos(t)] and this is the tangent vector to the curve. Now back to the general case, we can write: h(t) = f (~r(t)) = f (x(t), y(t)) and if h has an extreme point at t∗ , then h0 (t∗ ) = 0: 0 = h0 (t∗ ) = ∇f (x(t∗ ), y(t∗ )) · ~r0 (t) Thus, at an extreme point, the gradient of f is orthogonal to the tangent vector on the constraint, which is in turn orthogonal to ∇g. Method of Lagrange Multipliers If y = f (x) has an extreme point on the surface defined by g(x) = k, then candidates for where the extrema occur are where ∇f and ∇g are parallel. That is, ∇f (x) = λ∇g(x) The constant(s) λ are called the Lagrange multipliers. Let us solve our previous example: Find the optimal values of f (x, y) = x3 + y 3 + 3xy, such that (x − 3)2 + (y − 3)2 = 9. ∇f (x, y) = [3x2 + 3y, 3y 2 + 3x] ∇g(x, y) = [2(x − 3), 2(y − 3)] so that we solve the following system of equations: 3x2 + 3y = λ2(x − 3) 3y 2 + 3x = λ2(y − 3) 2 (x − 3) + (y − 3)2 = 9 12 These types of systems are generally difficult to solve by hand, so we’ll use a computer algebra system or Matlab to help us. In this example, we have: 3 x=3± √ , 2 3 y =3± √ , 2 33 λ = 21 ± √ 2 Note where these points occur in our Figure. The extrema are found by substituting these 2 values back into f . Doing this, we find that the maximum occurs at: 3 3 33 x = 3 + √ , y = 3 + √ , λ = 21 + √ 2 2 2 and the maximum value is approximately 347.32, which we expected from the contour plot. The minimum is approximately 3.67, which is also expected from the contour plot. Of course, in this example, we could actually convert f to a function of one variable by explicitly performing the parameterization. In this case, we have: f (t) = 27(cos(t) + 1)3 + 27(sin(t) + 1)3 + 27(cos(t) + 1)(sin(t) + 1) and we can use a numerical solver (like Newton’s method) to solve for where f 0 (t) = 0. In this case, we get that t ≈ .7854, 3.926 where we get a maximum value for f for the first value of t, and a minimum for f at the second value of t. This method can be employed for several constraints, g1 (x) = k1 , g2 (x) = k2 , . . . , gm (x) = km . In this case, the gradient of f is parallel to a linear combination of the gradients of gi , and we have the system of equations: ∇f g1 .. . = λ1 ∇g1 + λ2 ∇g2 + . . . + λm ∇gm = k1 .. . gm = k m Note that this is a system of n + m equations in n + m unknowns (if x ∈ Rn ). In general, we could translate this system into: G(x1 , . . . , xn , λ1 , . . . , λm ) = 0 and use Newton’s method to solve it numerically. 13 Example: Maximize 31 x31 + 13 x32 + 12 x33 subject to the conditions that x1 + x2 + x3 = 0 and x21 + x22 + x23 = 3. The system of equations we need to solve is given by: x21 x22 x23 x1 + x 2 + x 3 x21 + x22 + x23 = λ1 + 2λ2 x1 = λ1 + 2λ2 x2 = λ1 + 2λ2 x3 =0 =3 In this example, we have:    G=   x21 − 2λ2 x1 − λ1 x22 − 2λ2 x2 − λ1 x23 − 2λ2 x3 − λ1 x1 + x 2 + x 3 x21 + x22 + x23 − 3       and the Jacobian of G:   2x1 − 2λ2 0 0 −1 −2x1  0 2x2 − λ2 0 −1 −2x2     0 0 2x − 2λ −1 −2x JG =  3 2 3    1 1 1 0 0  2x1 2x2 2x3 0 0 √ Performing Newton’s method gives a maximum where λ1 = 1, λ2 = 42 , and √ triples of the form (− √12 , − √12 , 2), so any two variables can be − √12 , and √ the third is 2. 14