The Steepest Descent Algorithm for Unconstrained Optimization and a Bisection Line-Search Method

The Steepest Descent Algorithm for Unconstrained Optimization and a Bisection Line-search Method Robert M. Freund February, 2004 2004 Massachusetts Institute of Technology. 1 1 The Algorithm The problem we are interested in solving is: P : minimize f (x) x ∈ n , s.t. where f (x) is differentiable. If x = x̄ is a given point, f (x) can be approximated by its linear expansion f (¯ x + d) ≈ f (¯ x) + ∇f (¯ x)T d if d “small”, i.e., if d is small. Now notice that if the approximation in the above expression is good, then we want to choose d so that the inner product ∇f (x̄)T d is as small as possible. Let us normalize d so that d = 1. Then among all directions d with norm d = 1, the direction −∇f (x̄) d˜ = ∇f (x̄) makes the smallest inner product with the gradient ∇f (x̄). This fact follows from the following inequalities:
∇f (¯ x) d ≥ −∇f (¯ x x)d = ∇f (¯) T T −∇f (¯ x) ∇f (¯ x)  = −∇f (¯) x T d˜ . For this reason the un-normalized direction: d¯ = −∇f (x̄) is called the direction of steepest descent at the point x̄. Note that d¯ = −∇f (¯ x) is a descent direction as long as ∇f (¯ x) = 0. To T T ¯ x) = −(∇f (¯ x)) ∇f (¯ x) < 0 so long as see this, simply observe that d ∇f (¯ x) = 0. ∇f (¯ 2 A natural consequence of this is the following algorithm, called the steepest descent algorithm. Steepest Descent Algorithm: Step 0. Given x0 , set k := 0 Step 1. dk := −∇f (xk ). If dk = 0, then stop. Step 2. Solve minα f (xk + αdk ) for the stepsize αk , perhaps chosen by an exact or inexact linesearch. Step 3. Set xk+1 ← xk + αk dk , k ← k + 1. Go to Step 1. Note from Step 2 and the fact that dk = −∇f (xk ) is a descent direction, it follows that f (xk+1 ) < f (xk ). 2 Global Convergence We have the following theorem: Convergence Theorem: Suppose that f (·) : n →  is continuously differentiable on the set S = {x ∈ n | f (x) ≤ f (x0 )}, and that S is a closed and bounded set. Then every point x̄ that is a cluster point of the sequence {xk } satisfies ∇f (x̄) = 0. Proof: The proof of this theorem is by contradiction. By the Weierstrass ¯ be Theorem, at least one cluster point of the sequence {xk } must exist. Let x any such cluster point. Without loss of generality, assume that limk→∞ xk = ¯ but that ∇f (¯ x) =  0. This being the case, there is a value of α ¯ > 0 x, such that δ := f (¯ x + αd) ¯ ¯ > 0, where d¯ = −∇f (¯ x) − f (¯ x). Then also x + αd) ¯ ¯ ∈ intS. (¯ 3 Now limk→∞ dk = d¯. Then since (¯ x + αd) ¯ ¯ ∈ intS, and (xk + αd ¯ k) → (¯ x + αd), ¯ ¯ for k sufficiently large we have δ δ δ x) − δ + = f (¯ x) − . f (xk + αd ¯ k ) ≤ f (¯ x + αd) ¯ ¯ + = f (¯ 2 2 2 However, δ f (¯ x) < f (xk + αk dk ) ≤ f (xk + αd ¯ k ) ≤ f (¯ x) − , 2 which is of course a contradiction. Thus d¯ = −∇f (x̄) = 0. q.e.d. 3 The Rate of Convergence for the Case of a Quadratic Function In this section we explore answers to the question of how fast the steepest descent algorithm converges. We say that an algorithm exhibits linear convergence in the objective function values if there is a constant δ < 1 such that for all k sufficiently large, the iterates xk satisfy: f (xk+1 ) − f (x∗ ) ≤ δ, f (xk ) − f (x∗ ) where x∗ is some optimal value of the problem P . The above statement says that the optimality gap shrinks by at least δ at each iteration. Notice that if δ = 0.1, for example, then the iterates gain an extra digit of accuracy in the optimal objective function value at each iteration. If δ = 0.9, for example, then the iterates gain an extra digit of accuracy in the optimal objective function value every 22 iterations, since (0.9)22 ≈ 0.10. The quantity δ above is called the convergence constant. We would like this constant to be smaller rather than larger. We will show now that the steepest descent algorithm exhibits linear convergence, but that the convergence constant depends very much on the ratio of the largest to the smallest eigenvalue of the Hessian matrix H(x) at 4 the optimal solution x = x∗ . In order to see how this arises, we will examine the case where the objective function f (x) is itself a simple quadratic function of the form: 1 f (x) = xT Qx + q T x, 2 where Q is a positive definite symmetric matrix. We will suppose that the eigenvalues of Q are A = a1 ≥ a2 ≥ . . . ≥ an = a > 0, i.e, A and a are the largest and smallest eigenvalues of Q. The optimal solution of P is easily computed as: x∗ = −Q−1 q and direct substitution shows that the optimal objective function value is: 1 f (x∗ ) = − q T Q−1 q. 2 For convenience, let x denote the current point in the steepest descent algorithm. We have: 1 f (x) = xT Qx + q T x 2 and let d denote the current direction, which is the negative of the gradient, i.e., d = −∇f (x) = −Qx − q. Now let us compute the next iterate of the steepest descent algorithm. If α is the generic step-length, then 1 f (x + αd) = (x + αd)T Q(x + αd) + q T (x + αd) 2 5 1 1 = xT Qx + αdT Qx + α2 dT Qd + q T x + αq T d 2 2 1 = f (x) − αdT d + α2 dT Qd. 2 Optimizing the value of α in this last expression yields α= dT d , dT Qd and the next iterate of the algorithm then is
x = x + αd = x + dT d d, dT Qd and 1 1 (dT d)2
f (x ) = f (x + αd) = f (x) − αdT d + α2 dT Qd = f (x) − . 2 2 dT Qd Therefore, T 2 d)
f (x) − 12 (d − f (x∗ ) f (x ) − f (x∗ ) dT Qd = f (x) − f (x∗ ) f (x) − f (x∗ ) =1− =1− 1 T 2 x Qx 1 2 (Qx =1− 1 (dT d)2 2 dT Qd + q T x + 12 q T Q−1 q + 1 (dT d)2 2 dT Qd q)T Q−1 (Qx (dT d)2 (dT Qd)(dT Q−1 d) =1− where β= + q) 1 β (dT Qd)(dT Q−1 d) . (dT d)2 6 In order for the convergence constant to be good, which will translate to fast linear convergence, we would like the quantity β to be small. The following result provides an upper bound on the value of β. Kantorovich Inequality: Let A and a be the largest and the smallest eigenvalues of Q, respectively. Then β≤ (A + a)2 . 4Aa We will prove this inequality later. For now, let us apply this inequality to the above analysis. Continuing, we have
4Aa 1 (A − a)2 f (x ) − f (x∗ ) = 1 − ≤ 1 − = = f (x) − f (x∗ ) β (A + a)2 (A + a)2
A/a − 1 A/a + 1 2 =: δ . Note by definition that A/a is always at least 1. If A/a is small (not much bigger than 1), then the convergence constant δ will be much smaller than 1. However, if A/a is large, then the convergence constant δ will be only slightly smaller than 1. Table 1 shows some sample values. Note that the number of iterations needed to reduce the optimality gap by a factor of 10 grows linearly in the ratio A/a. 4 Examples 4.1 Example 1: Typical Behavior Consider the function f (x1 , x2 ) = 5x21 + x22 + 4x1 x2 − 14x1 − 6x2 + 20. This function has its optimal solution at x∗ = (x1∗ , x∗2 ) = (1, 1) and f (1, 1) = 10. In step 1 of the steepest descent algorithm, we need to compute  k d = −∇f (xk1 , xk2 ) = −10xk1 − 4xk2 + 14 −2xk2 − 4xk1 + 6 7   = dk1 dk2  Upper Bound on  A 1.1 3.0 10.0 100.0 200.0 400.0 a 1.0 1.0 1.0 1.0 1.0 1.0 δ=  A/a−1 2 A/a+1 0.0023 0.25 0.67 0.96 0.98 0.99 Number of Iterations to Reduce the Optimality Gap by a factor of 10 1 2 6 58 116 231 Table 1: Sensitivity of Steepest Descent Convergence Rate to the Eigenvalue Ratio and in step 2 we need to solve αk = arg minα h(α) = f (xk + αdk ). In this example we will be able to derive an analytic expression for αk . Notice that h(α) = f (xk + αdk ) = 5(xk1 + αdk1 )2 + (xk2 + αdk2 )2 + 4(xk1 + αdk1 )(xk2 + αdk2 ) − −14(xk1 + αdk1 ) − 6(xk2 + αdk2 ) + 20, and this is a simple quadratic function of the scalar α. It is minimized at αk = (dk1 )2 + (dk2 )2 2(5(dk1 )2 + (dk2 )2 + 4dk1 dk2 ) Using the steepest descent algorithm to minimize f (x) starting from x1 = (x11 , x12 ) = (0, 10), and using a tolerance of  = 10−6 , we compute the iterates shown in Table 2 and in Figure 2: For a convex quadratic function f (x) = 12 xT Qx−cT x, the contours of the function values will be shaped like ellipsoids, and the gradient vector ∇f (x) at any point x will be perpendicular to the contour line passing through x, see Figure 1. 8 k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 xk1 0.000000 −2.252782 0.755548 0.204852 0.940243 0.805625 0.985392 0.952485 0.996429 0.988385 0.999127 0.997161 0.999787 0.999306 0.999948 0.999830 0.999987 0.999959 0.999997 0.999990 0.999999 0.999998 1.000000 0.999999 xk2 10.000000 8.786963 3.200064 2.903535 1.537809 1.465322 1.131468 1.113749 1.032138 1.027806 1.007856 1.006797 1.001920 1.001662 1.000469 1.000406 1.000115 1.000099 1.000028 1.000024 1.000007 1.000006 1.000002 1.000001 dk1 −26.000000 1.379968 −6.355739 0.337335 −1.553670 0.082462 −0.379797 0.020158 −0.092842 0.004928 −0.022695 0.001205 −0.005548 0.000294 −0.001356 0.000072 −0.000332 0.000018 −0.000081 0.000004 −0.000020 0.000001 −0.000005 0.000000 dk2 −14.000000 −2.562798 −3.422321 −0.626480 −0.836592 −0.153144 −0.204506 −0.037436 −0.049992 −0.009151 −0.012221 −0.002237 −0.002987 −0.000547 −0.000730 −0.000134 −0.000179 −0.000033 −0.000044 −0.000008 −0.000011 −0.000002 −0.000003 −0.000000 ||dk ||2 29.52964612 2.91071234 7.21856659 0.71152803 1.76458951 0.17393410 0.43135657 0.04251845 0.10544577 0.01039370 0.02577638 0.00254076 0.00630107 0.00062109 0.00154031 0.00015183 0.00037653 0.00003711 0.00009204 0.00000907 0.00002250 0.00000222 0.00000550 0.00000054 Table 2: Iterations of the example in Subsection 4.1. 9 αk 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1800 0.0866 2.1803 0.0866 2.1817 0.0866 0.0000 f (xk ) 60.000000 22.222576 12.987827 10.730379 10.178542 10.043645 10.010669 10.002608 10.000638 10.000156 10.000038 10.000009 10.000002 10.000001 10.000000 10.000000 10.000000 10.000000 10.000000 10.000000 10.000000 10.000000 10.000000 10.000000 2 2 5 x +4 y +3 x y+7 x+20 1200 1000 800 600 400 200 0 10 5 0 −5 −10 y 0 −5 −10 5 10 x Figure 1: The contours of a convex quadratic function are ellipsoids. 4.2 Example 2: Hem-Stitching Consider the function 1 f (x) = xT Qx − cT x + 10 2 where Q and c are given by:  Q= 20 5 5 2  and c= 14 6 10   10 8 6 4 2 0 −2 −5 0 5 Figure 2: Behavior of the steepest descent algorithm for the example of Subsection 4.1. For this problem, A a = 30.234, and so the linear convergence rate δ of  the steepest descent algorithm is bounded above by δ =  A/a−1 2 A/a+1 = 0.8760. If we apply the steepest descent algorithm to minimize f (x) starting from x1 = (40, −100), we obtain the iteration values shown in Table 3 and shown graphically in Figure 3. 4.3 Example 3: Rapid Convergence Consider the function 1 f (x) = xT Qx − cT x + 10 2 11 xk2 ||dk ||2 40.000000 −100.000000 25.542693 −99.696700 26.277558 −64.668130 16.763512 −64.468535 17.247111 −41.416980 10.986120 −41.285630 11.304366 −26.115894 7.184142 −26.029455 7.393573 −16.046575 4.682141 −15.989692 0.460997 0.948466 −0.059980 3.038991 −0.124280 3.297005 −0.132216 3.328850 −0.133195 3.332780 −0.133316 3.333265 −0.133331 3.333325 −0.133333 3.333332 286.06293014 77.69702948 188.25191488 51.13075844 123.88457127 33.64806192 81.52579489 22.14307211 53.65038732 14.57188362 1.79847660 0.22196980 0.02739574 0.00338121 0.00041731 0.00005151 0.00000636 0.00000078 xk1 k 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 f (xk ) f (xk )−f (x∗ ) f (xk−1 )−f (x∗ ) 6050.000000 3981.695128 2620.587793 1724.872077 1135.420663 747.515255 492.242977 324.253734 213.703595 140.952906 3.066216 0.965823 0.933828 0.933341 0.933333 0.933333 0.933333 0.933333 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658079 0.658078 0.658025 0.654656 0.000000 αk 0.0506 0.4509 0.0506 0.4509 0.0506 0.4509 0.0506 0.4509 0.0506 0.4509 0.4509 0.4509 0.4509 0.4509 0.4509 0.4509 0.4509 0.0000 Table 3: Iteration values for example of Subsection 4.2, which shows hemstitching. where Q and c are given by:  Q=  and c= For this problem, 20 5 5 16 A a 14 6   = 1.8541, and so the linear convergence rate δ of
A the steepest descent algorithm is bounded above by δ = 12 −1 a A +1 a 2 = 0.0896. Figure 3: Hem-stitching in the example of Subsection 4.2. If we apply the steepest descent algorithm to minimize f (x) starting from x1 = (40, −100), we obtain the iteration values shown in Table 4 and shown graphically in Figure 4. Figure 4: Rapid convergence in the example of Subsection 4.3. 13 xk2 ||dk ||2 40.000000 −100.000000 19.867118 −1.025060 2.513241 −4.555081 1.563658 0.113150 0.745149 −0.053347 0.700361 0.166834 0.661755 0.158981 0.659643 0.169366 0.657822 0.168996 0.657722 0.169486 0.657636 0.169468 0.657632 0.169491 0.657628 0.169490 0.657627 0.169492 0.657627 0.169491 1434.79336491 385.96252652 67.67315150 18.20422450 3.19185713 0.85861649 0.15054644 0.04049732 0.00710064 0.00191009 0.00033491 0.00009009 0.00001580 0.00000425 0.00000075 xk1 k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 f (xk ) f (xk )−f (x∗ ) f (xk−1 )−f (x∗ ) 0.0704 76050.000000 0.0459 3591.615327 0.0704 174.058930 0.0459 12.867208 0.0704 5.264475 0.0459 4.905886 0.0704 4.888973 0.0459 4.888175 0.0704 4.888137 0.0459 4.888136 0.0704 4.888136 0.0459 4.888136 0.0704 4.888136 0.0459 4.888136 0.0000 4.888136 0.047166 0.047166 0.047166 0.047166 0.047166 0.047166 0.047166 0.047166 0.047166 0.047166 0.047161 0.047068 0.045002 0.000000 αk Table 4: Iteration values for the example in Subsection 4.3, which exhibits rapid convergence. 4.4 Example 4: A nonquadratic function Let f (x) = x1 − 0.6x2 + 4x3 + 0.25x4 − 4  log(xi ) − log(5 − i=1 4  xi ) . i=1 Table 5 shows values of the iterates of steepest descent applied to this function. 4.5 Example 5: An analytic example Suppose that 1 f (x) = xT Qx + q T x 2 14 k 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 110 120 121 122 xk1 xk2 xk3 xk4 ||dk ||2 f (xk ) f (xk )−f (x∗ ) f (xk−1 )−f (x∗ ) 1.000000 0.802973 0.634501 0.616707 0.545466 0.543856 0.553037 0.547600 0.526273 0.525822 0.511244 0.503892 0.501351 0.500467 0.500161 0.500056 0.500019 0.500007 0.500002 0.500001 0.500000 0.500000 0.500000 1.000000 1.118216 1.710704 1.735137 1.972851 1.986648 2.121519 2.129461 2.205845 2.212935 2.406747 2.468138 2.489042 2.496222 2.498696 2.499550 2.499845 2.499946 2.499981 2.499993 2.499998 2.499998 2.499998 1.000000 0.211893 0.332224 0.205759 0.267358 0.204195 0.241186 0.202091 0.229776 0.202342 0.200739 0.200271 0.200095 0.200033 0.200011 0.200004 0.200001 0.200000 0.200000 0.200000 0.200000 0.200000 0.200000 1.000000 0.950743 1.121308 1.108079 1.054072 1.044882 0.991524 0.983563 0.938380 0.933770 0.841546 0.813922 0.804762 0.801639 0.800565 0.800195 0.800067 0.800023 0.800008 0.800003 0.800001 0.800001 0.800001 4.17402683 1.06742574 1.88344095 0.46928947 1.15992055 0.31575186 0.80724427 0.23764416 0.58321024 0.18499125 0.06214411 0.02150306 0.00743136 0.00256658 0.00088621 0.00030597 0.00010564 0.00003647 0.00001257 0.00000448 0.00000145 0.00000455 0.00000098 4.650000 2.276855 1.939389 1.797957 1.739209 1.698071 1.674964 1.657486 1.646479 1.637766 1.612480 1.609795 1.609480 1.609443 1.609439 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 0.219505 0.494371 0.571354 0.688369 0.682999 0.739299 0.733266 0.770919 0.764765 0.803953 0.806953 0.807905 0.808228 0.808338 0.808376 0.808381 0.808271 0.806288 0.813620 0.606005 0.468218 0.000000 Table 5: Iteration values for the example of Subsection 4.4. 15
where Q= +4 −2 −2 +2
Then ∇f (x) = 
and +4 −2 −2 +2 
x1 x2 
+  +2 −2 . 
 and so ∗ x = and q= +2 −2 0 1 f (x∗ ) = −1. Direct computation shows that the eigenvalues of Q are A = 3 + √ a = 3 − 5, whereby the bound on the convergence constant is
δ= A/a − 1 A/a + 1 2 √ 5 and ≤ 0.556. Suppose that x0 = (0, 0). Then we have: x1 = (−0.4, 0.4), x2 = (0, 0.8) and the even numbered iterates satisfy x2n = (0, 1 − 0.2n ) and so and f (x2n ) = (1 − 0.2n )2 − 2 + 2(0.2)n f (x2n ) − f (x∗ ) = (0.2)2n . Therefore, starting from the point x0 = (0, 0), the optimality gap goes down by a factor of 0.04 after every two iterations of the algorithm. 5 Final Remarks • The proof of convergence is basic in that it only relies on fundamental properties of the continuity of the gradient and on the closedness and boundedness of the level sets of the function f (x). • The analysis of the rate of convergence is quite elegant. 16 • Kantorovich won the Nobel Prize in Economics in 1975 for his contributions to the application of linear programming to economic theory. • The bound on the rate of convergence is attained in practice quite often, which is too bad. • The ratio of the largest to the smallest eigenvalue of a matrix is called the condition number of the matrix. The condition number plays an extremely important role in the theory and the computation of quantities involving the matrix Q. • What about non-quadratic functions? Most functions behave as nearquadratic functions in a neighborhood of the optimal solution. The analysis of the non-quadratic case gets very involved; fortunately, the key intuition is obtained by analyzing the quadratic case. 6 A Bisection Algorithm for a Line-Search of a Convex Function Suppose that f (x) is a continuously differentiable convex function, and that we seek to solve: α ¯ := arg min f (¯ x + αd¯), α where x̄ is our current iterate, and d¯ is the current direction generated by an algorithm that seeks to minimize f (x). Suppose further that d¯ is a descent direction of f (x) at x = x ¯, namely: f (¯ x + d¯) < f (¯ x) for all  > 0 and sufficiently small . Let h(α) := f (x̄ + αd¯), whereby h(α) is a convex function in the scalar variable α, and our problem is to solve for ᾱ = arg min h(α). α We therefore seek a value ᾱ for which
h (ᾱ) = 0. 17 It is elementary to show that
¯ T d. ¯ x + αd) h (α) = ∇f (¯
Proposition 6.1 h (0) < 0. q.e.d. Because h(α) is a convex function of α, we also have:
Proposition 6.2 h (α) is a monotone increasing function of α. q.e.d. Figure 5 shows an example of convex function of two variables to be optimized. Figure 6 shows the function h(α) obtained by restricting the function of Figure 5 to the line shown in that figure. Note from Figure 6 that
h(α) is convex. Therefore its first derivative h (α) will be a monotonically increasing function. This is shown in Figure 7.
Because h (α) is a monotonically increasing function, we can approxi
α) = 0, by a suitable bisection mately compute α, ¯ the point that satisfies h (¯

method. Suppose that we know a value α ˆ that h (ˆ α) > 0. Since h (0) < 0
α) > 0, the mid-value α ˜ = 0+2 αˆ is a suitable test-point. Note the and h (ˆ following:
• If h (α̃) = 0, we are done.
α) > 0, we can now bracket α ¯ in the interval (0, α). ˜ • If h (˜
α) < 0, we can now bracket α ¯ in the interval (α, ˜ α). ˆ • If h (˜ This leads to the following bisection algorithm for minimizing h(α) = f (x̄ + ¯ by solving the equation h
(α) ≈ 0. αd) Step 0. Set k = 0. Set αl := 0 and αu := α̂. Step k. Set α ˜= αu +αl 2
and compute h (˜ α). 18 20 10 0 −10 −20 f(x ,x−30) 1 2 −40 −50 −60 −70 10 −80 0 8 0.5 6 1 x 2 4 1.5 2 2 0 2.5 3 −2 Figure 5: A convex function to be optimized. 19 x1 0 −10 −20 h(α)−30 −40 −50 −60 −0.4 −0.2 0 0.2 α 0.4 0.6 0.8 1 Figure 6: The 1-dimensional function h(α). 1 0.5 0 h′(α)−0.5 −1 −1.5 −2 −0.4 −0.2 0 0.2
α 0.4 0.6 0.8 1 Figure 7: The function h (α) is monotonically increasing. 20
• If h (˜ α) > 0, re-set αu := α. ˜ Set k ← k + 1.
• If h (˜ α) < 0, re-set αl := α. ˜ Set k ← k + 1.
• If h (α̃) = 0, stop. Proposition 6.3 After every iteration of the bisection algorithm, the cur
¯ such that h (¯ α) = 0. rent interval [αl , αu ] must contain a point α Proposition 6.4 At the k th iteration of the bisection algorithm, the length of the current interval [αl , αu ] is
k L= 1 2 (α̂). Proposition 6.5 A value of α such that |α − ᾱ| ≤  can be found in at most 
 α̂ log2  steps of the bisection algorithm. 6.1
Computing α ˆ for which h (ˆ α) > 0
Suppose that we do not have available a convenient value α ˆ for which h (ˆ α) >
α). 0. One way to proceed is to pick an initial “guess” of α ˆ and compute h (ˆ

If h (ˆ α) > 0, then proceed to the bisection algorithm; if h (ˆ α) ≤ 0, then re-set α ˆ ← 2ˆ α and repeat the process. 6.2 Stopping Criteria for the Bisection Algorithm In practice, we need to run the bisection algorithm with a stopping criterion. Some relevant stopping criteria are: ¯ • Stop after a fixed number of iterations. That is stop when k = K, ¯ where K specified by the user. • Stop when the interval becomes small. That is, stop when αu − αl ≤ , where  is specified by the user.

α)| becomes small. That is, stop when |h (˜ α)| ≤ , • Stop when |h (˜ where  is specified by the user. This third stopping criterion typically yields the best results in practice. 21 6.3 Modification of the Bisection Algorithm when the Do­ main of f (x) is Restricted The discussion and analysis of the bisection algorithm has presumed that our optimization problem is P : minimize f (x) x ∈ n . s.t. Given a point x̄ and a direction d¯, the line-search problem then is LS : minimize h(α) := f (x̄ + αd¯) α ∈ . s.t. Suppose instead that the domain of definition of f (x) is an open set X ⊂ n . Then our optimization problem is: P : minimize f (x) x ∈ X, s.t. and the line-search problem then is LS : minimize h(α) := f (x̄ + αd¯) s.t. x ¯ + αd¯ ∈ X. In this case, we must ensure that all iterate values of α in the bisection algorithm satisfy x ¯ + αd¯ ∈ X. As an example, consider the following problem: P : minimize f (x) := − s.t. m i=1 b − Ax > 0. 22 ln(bi − Ai x) Here the domain of f (x) is X = {x ∈ n | b − Ax > 0}. Given a point x̄ ∈ X and a direction d¯, the line-search problem is: LS : minimize h(α) := f (¯ x + αd¯) = − s.t. m i=1 x + αd¯)) ln(bi − Ai (¯ b − A(x̄ + αd¯) > 0. Standard arithmetic manipulation can be used to establish that x + αd¯) > 0 if and only if α ˆ b − A(¯ ˇ<α<α where α ˇ := − min ¯ Ai d<0 bi − A i x ¯ ¯ −Ai d ˆ := min and α ¯ Ai d>0 bi − Ai x ¯ , ¯ Ai d and the line-search problem then is: LS : minimize h(α) := − s.t. m i=1 α ˇ < α < α. ˆ 23 ln(bi − Ai (x̄ + αd¯)) 7 Proof of Kantorovich Inequality Kantorovich Inequality: Let A and a be the largest and the smallest eigenvalues of Q, respectively. Then β≤ (A + a)2 . 4Aa Proof: Let Q = RDRT , and then Q−1 = RD−1 RT , where R = RT is an orthonormal matrix, and the eigenvalues of Q are 0 < a = a1 ≤ a2 ≤ . . . ≤ an = A, and, ⎛ ⎞ a1 ⎜ 0 ⎜ D = ⎜ .. ⎝ . 0 a2 .. . ... ... .. . 0 0 .. . 0 0 . . . an ⎟ ⎟ ⎟. ⎠ Then β= f T Df f T D−1 f (dT RDRT d)(dT RD−1 RT d) = (dT RRT d)(dT RRT d) fT ffT f where f = RT d. Let λi = β= n  i=1 fi2 . fT f λi ai Then λi ≥ 0 and n  i=1
λi 1 ai  i=1 λi = 1, and  n ⎛ = i=1 ⎜ ⎝ n λi n i=1 1 ai 1  ⎞. ⎟ ⎠ λi ai The largest value of β is when λ1 + λn = 1, see the illustration in Figure 8. Therefore, 24 β≤ λ1 a1 + λn A1 1 λ1 a+λn A = ( 1 A + 21 a)( 12 a + 21 A) (A + a)2 (λ1 a + λn A)(λ1 A + λn a) ≤ 2 = . Aa Aa 4Aa q.e.d. y y = 1/x Σλi(1/ai) 1/(Σλiai) 0 Σλiai 0 a1 a2 a3 =a ... an =A x Figure 8: Illustration for proof of the Kantorovich Inequality. 25 8 Steepest Descent Exercises NOTE ON COMPUTATION EXERCISES: You may use any machine and any programming language. We recommend, however, that you use MATLAB on Athena or on your own computer. (The teaching assistant is prepared to help you get started in MATLAB on Athena.) Some of the specific details in the computation exercises are purposely left to your own discretion. For example, in cases where you must use a line-search routine, you must choose your own tolerance for the line-search. Also, you must choose your own stopping criteria for your computation. Please include your code or pseudo-code in your write-up. 1. (Steepest Descent) Suppose that xk and xk+1 are two consecutive points generated by the steepest descent algorithm with exact linesearch. Show that ∇f (xk )T ∇f (xk+1 ) = 0. 2. (Steepest Descent) Suppose that we seek to minimize f (x1 , x2 ) = 5x21 + 5x22 − x1 x2 − 11x1 + 11x2 + 11. (a) Find a point satisfying the first-order necessary conditions for a solution. (b) Show that this point is a global minimum of f (x). (c) What would be the worst rate of convergence for the steepest descent algorithm for this problem? (d) Starting at (x1 , x2 ) = (0, 0), at most how many steepest descent iterations would it take to reduce the function value to 10−11 ? 3. (Steepest Descent) Suppose we seek to minimize 1 f (x) = xT Hx + cT x + 13, 2  where H= 10 −9 −9 10   and c = 4 −15  . Implement the steepest descent algorithm on this problem, using the following starting points: 26 • x0 = (0, 0)T . • x0 = (−0.4, 0)T . • x0 = (10, 0)T . • x0 = (11, 0)T . As it turns out, the optimal solution to this problem is x∗ = (5, 6)T , with f (x∗ ) = −22. What linear convergence constants do you observe for each of the above starting points? 4. (Steepest Descent) Suppose we seek to minimize 1 f (x) = xT Hx + cT x, 2 where ⎛ ⎞ ⎛ ⎞ 10 −18 2 12 ⎜ ⎟ ⎜ ⎟ H = ⎝ −18 40 −1 ⎠ and c = ⎝ −47 ⎠ . 2 −1 3 −8 Implement the steepest descent algorithm on this problem, using the following starting points. • x0 = (0, 0, 0)T . • x0 = (15.09, 7.66, −6.56)T . • x0 = (11.77, 6.42, −4.28)T . • x0 = (4.46, 2.25, 1.85)T . As it turns out, the optimal solution to this problem is x∗ = (4, 3, 1)T , with f (x∗ ) = −50.5. What linear convergence constants do you observe for each of the above starting points? 5. (Steepest Descent) Suppose that we seek to minimize the following function: f (x1 , x2 ) = −9x1 −10x2 +θ(− ln(100−x1 −x2 )−ln(x1 )−ln(x2 )−ln(50−x1 +x2 )), where θ is a given parameter. Note that the domain of this function is X = {(x1 , x2 ) | x1 > 0, x2 > 0, x1 + x2 < 100, x1 − x2 < 50}. Implement the steepest descent algorithm for this problem, using the 27 bisection algorithm for your line-search, with the stopping criterion
that |h (α̃)| ≤  = 10−6 . Run your algorithm for θ = 10 and for θ = 100, using the following starting points. • x0 = (8, 90)T . • x0 = (1, 40)T . • x0 = (15, 68.69)T . • x0 = (10, 20)T . What linear convergence constants do you observe for each of the above starting points? Helpful Hint: Note that f (x1 , x2 ) is only defined on the domain X = {(x1 , x2 ) | x1 > 0, x2 > 0, x1 + x2 < 100, x1 − x2 < 50}. If you are at a point x ¯ = (¯ x1 , x ¯2 ) and you have computed a direction d¯ = (d¯1 , d¯2 ), then note that a value of the upper bound α̂ is effectively given by the largest value of α for which the following constraints are satisfied: ¯2 +αd¯2 > 0, x ¯1 +αd¯1 +¯ x2 +αd¯2 < 100, x ¯1 +αd¯1 −x ¯2 −αd¯2 < 50. x ¯1 +αd¯1 > 0, x The largest value of α satisfying these conditions can easily be computed by applying appropriate logic. Important Note: Please keep the value of θ as a specific command in your code. The reason for this is that later in the semester, we will consider an algorithm where we iteratively change the value of θ. It will then be useful for your code to have a line in it where θ can be set and re-set. 28