Simple Routines for Optimization

1 Outline • A Bisection Line-Search Algorithm for 1-Dimensional Optimization • The Conditional-Gradient Method for Constrained Optimization (FrankWolfe Method) • Subgradient Optimization • Application of Subgradient Optimization to the Lagrange Dual Problem 2 A Bisection Line-Search Algorithm for 1-Dimensional Optimization Consider the optimization problem: P : minimizex f (x) x ∈ n . s.t. Let us suppose that f (x) is a differentiable convex function. In a typical algorithm for solving P we have a current iterate value x̄ and we choose a direction d¯ by some suitable means. The direction d¯ is usually chosen to be a descent direction, defined by the following property: f (¯ x) for all
> 0 and sufficiently small . x +
d¯) < f (¯ We then typically also perform the 1-dimensional line-search optimization: ¯ := arg min f (¯ x + αd¯) . α α 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(α). α 2 We therefore seek a value ᾱ for which � h (ᾱ) = 0. It is elementary to show that � ¯ ¯ T d. x + αd) h (α) = ∇f (¯ � Property: If d¯ is a descent direction at x̄, then h (0) < 0. Because h(α) is a convex function of α, we also have: � Property: h (α) is a monotone increasing function of α. Figure 1 shows an example of convex function of two variables to be optimized. Figure 2 shows the function h(α) obtained by restricting the function of Figure 1 to the line shown in that figure. Note from Figure 2 that � h(α) is convex. Therefore its first derivative h (α) will be a monotonically increasing function. This is shown in Figure 3. � Because h (α) is a monotonically increasing function, we can approxi� ¯ the point that satisfies h (¯ α) = 0, by a suitable bisection mately compute α, � � method. Suppose that we know a value α ˆ that h (ˆ α) > 0. Since h (0) < 0 � and h (ˆ α) > 0, the mid-value α ˜ = 0+2 αˆ is a suitable test-point. Note the following: � • If h (α̃) = 0, we are done. � • If h (˜ α) > 0, we can now bracket α ¯ in the interval (0, α). ˜ � • If h (˜ α) < 0, we can now bracket α ¯ in the interval (α, ˆ ˜ α). 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 (˜ α). � α) > 0, re-set αu := α. ˜ Set k ← k + 1. • If h (˜ � • If h (˜ α) < 0, re-set αl := α. ˜ Set k ← k + 1. 3 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 1: A convex function to be optimized. 4 x 1 0 −10 −20 h(α)−30 −40 −50 −60 −0.4 −0.2 0 0.2 α 0.4 0.6 Figure 2: The 1-dimensional function h(α). 5 0.8 1 1 0.5 0 h′(α)−0.5 −1 −1.5 −2 −0.4 −0.2 0 0.2 � α 0.4 0.6 Figure 3: The function h (α) is monotonically increasing. 6 0.8 1 � • If h (α̃) = 0, stop. Property: After every iteration of the bisection algorithm, the current � ¯ such that h (¯ α) = 0. interval [αl , αu ] must contain a point α Property: At the k th iteration of the bisection algorithm, the length of the current interval [αl , αu ] is L=
k 1 2 (α̂). Property: A value of α such that |α − ᾱ| ≤
can be found in at most 
 α̂
log2 steps of the bisection algorithm. 2.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. 2.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 is 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. 7 2.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 : minimizex 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 : minimizex 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 : minimizex f (x) := − s.t. m  i=1 b − Ax > 0. 8 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: m  x + αd¯) = − ln(bi − Ai (¯ x + αd¯)) LS : minimizeα h(α) := f (¯ i=1 b − A(x̄ + αd¯) > 0. s.t. Standard arithmetic manipulation can be used to establish that b − A(¯ x + αd¯) > 0 if and only if α ˇ<α<α ˆ 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. 3 m  i=1 ln(bi − Ai (x̄ + αd¯)) α ˇ < α < α. ˆ The Conditional-Gradient Method for Constrained Optimization (Frank-Wolfe Method) We now consider the following optimization problem: P : minimizex f (x) x∈C . s.t. We assume that f (x) is a convex function, and that C is a convex set. Herein we describe the conditional-gradient method for solving P , also called the Frank-Wolfe method. This method is one of the cornerstones of optimization, and was one of the first successful algorithms used to solve nonlinear optimization problems. It is based on the premise that the set C 9 is well-suited for linear optimization. That means that either C is itself a system of linear inequalities C = {x | Ax ≤ b}, or more generally that the problem: LOc : minimizex cT x s.t. x∈C is easy to solve for any given objective function vector c. This being the case, suppose that we have a given iterate value x̄ ∈ C. Let us linearize the function f (x) at x = x̄. This linearization is: x) + ∇f (¯ x)T (x − x) ¯ , z1 (x) := f (¯ which is the first-order Taylor expansion of f (·) at x. ¯ Since we can easily do linear optimization on C, let us solve: x) + ∇f (¯ x)T (x − x) ¯ LP : minimizex z1 (x) = f (¯ s.t. x∈C , which simplifies to: x)T x LP : minimizex ∇f (¯ x∈C . s.t. Let x∗ denote the optimal solution to this problem. Then since C is a convex set, the line segment joining x ¯ and x∗ is also in C, and we can perform a line-search of f (x) over this segment. That is, we solve: x + α(x∗ − x)) ¯ LS : minimizeα f (¯ s.t. 0≤α≤1. Let α ¯ denote the solution to this line-search problem. We re-set x: ¯ 10 x ¯←x ¯ + α(x ¯ ∗ − x) ¯ and repeat this process. The formal description of this method, called the conditional gradient method or the Frank-Wolfe method, is given below: Step 0: Initialization. Start with a feasible solution x0 ∈ C. Set k = 0. Set LB ← −∞. Step 1: Update upper bound. Set U B ← f (xk ). Set x̄ ← xk . Step 2: Compute next iterate. – Solve the problem x) + ∇f (¯ x)T (x − x) ¯ z̄ = minx f (¯ x∈C , s.t. and let x∗ denote the solution. – Solve the line-search problem: ¯ minimizeα f (x̄ + α(x∗ − x)) 0≤α≤1, s.t. and let ᾱ denote the solution. ¯ + α(x ¯ ¯ ∗ − x) – Set xk+1 ← x Step 3: Update Lower Bound. Set LB ← max{LB, z̄}. Step 4: Check Stopping Criteria. If |U B − LB| ≤
, stop. Otherwise, set k ← k + 1 and go to Step 1. 11 3.1 Upper and Lower Bounds in the Frank-Wolfe Method, and Convergence • The upper bound values U B are simply the objective function values of the iterates f (xk ) for k = 0, . . .. This is a monotonically decreasing sequence because the line-search guarantees that each iterate is an improvement over the previous iterate. • The lower bound values LB result from the convexity of f (x) and the gradient inequality for convex functions: f (x) ≥ f (¯ x) + ∇f (¯ x)T (x − x) ¯ for any x ∈ C . Therefore min f (x) ≥ min f (¯ x) + ∇f (¯ ¯ = z¯ , x)T (x − x) x∈C x∈C and so the optimal objective function value of P is bounded below by z̄. We also have the following convergence theorem for the Frank-Wolfe method: Property: Suppose that C is a bounded set, and that there exists a constant L for which ∇f (x) − ∇f (y) ≤ L x − y for all x, y ∈ C. Then there exists a constant Ω > 0 for which the following is true: Ω . f (xk ) − min f (x) ≤ x∈C k 3.2 Illustration of the Frank-Wolfe Method Consider the following instance of P : P : minimize f (x) s.t. where 12 x∈C , f (x) = f (x1 , x2 ) = −32x1 + x41 − 8x2 + x22 and C = {(x1 , x2 ) | x1 − x2 ≤ 1, 2.2x1 + x2 ≤ 7, x1 ≥ 0, x2 ≥ 0} . Notice that the gradient of f (x1 , x2 ) is given by the formula:  ∇f (x1 , x2 ) = 4x31 − 32 2x2 − 8 . Suppose that xk = x ¯ = (0.5, 3.0) is the current iterate of the FrankWolfe method, and the current lower bound is LB = −100.0. We compute x) = f (0.5, 3.0) = −30.9375 and we compute the gradient of f (x) at x: ¯ f (¯  ∇f (0.5, 3.0) = 4x31 − 32 2x2 − 8  = −31.5 −2.0 . We then create and solve the following linear optimization problem: LP : z̄ = minx1 ,x2 s.t. −30.9375 − 31.5(x1 − 0.5) − 2.0(x2 − 3.0) x1 − x2 ≤ 1 2.2x1 + x2 ≤ 7 x1 ≥ 0 x2 ≥ 0 . The optimal solution of this problem is: x∗ = (x∗1 , x2∗ ) = (2.5, 1.5) , and the optimal objective function value is: z̄ = −50.6875 . Now we perform a line-search of the 1-dimensional function = −32(¯ ¯ x1 + α(x∗1 − x ¯1 )) + (¯ x1 + α(x∗1 − x ¯1 ))4 f (¯ x + α(x∗ − x)) −8(¯ x2 + α(x∗2 − x ¯2 )) + (¯ x2 + α(x∗2 − x ¯2 ))2 13 over α ∈ [0, 1]. This function attains its minimum at α ¯ = 0.7165 and we therefore update as follows: xk+1 ← x+ ¯ ∗ −¯ x) = (0.5, 3.0)+0.7165((2.5, 1.5)−(0.5, 3.0)) = (1.9329, 1.9253) ¯ α(x and LB ← max{LB, z̄} = max{−100, −50.6875} = −50.6875 . The new upper bound is U B = f (xk+1 ) = f (1.9329, 1.9253) = −59.5901 . This is illustrated in Figure 4. 4 4.1 Subgradient Optimization Definition Suppose that f (x) is a convex function. If f (x) is differentiable, we have the gradient inequality: f (x) ≥ f (¯ x) + ∇f (¯ x)T (x − x) ¯ for any x ∈ X , where typically we think of X = n . This inequality is illustrated in Figure 5. There are many important convex functions that are not differentiable. The notion of the gradient generalizes to the concept of a subgradient of a convex function. A vector g ∈ n is called subgradient of the convex function f (x) at x = x̄ if the following inequality is satisfied: ¯ f (x) ≥ f (¯ x) + g T (x − x) This definition is illustrated in Figure 6. 14 for all x ∈ X . 8 7 6 5 4 x2 xk 3 x^ 2 k+1 x x* 1 0 −1 −1 −0.5 0 0.5 1 1.5 2 2.5 x1 Figure 4: Illustration of an iteration of the Frank-Wolfe method. 15 3 5000 4000 3000 2000 f(u)1000 0 −1000 −2000 −3000 −100 −80 −60 −40 −20 0 20 40 60 80 u Figure 5: The gradient and the gradient inequality for a differentiable convex function. 16 100 8000 6000 4000 2000 f(u) 0 −2000 −4000 −6000 −100 −80 −60 −40 −20 0 20 40 60 80 u Figure 6: Subgradients and the subgradient inequality for a nondifferentiable convex function. 17 100 4.2 Properties of Subgradients Suppose that f (x) is a convex function. For each x, let ∂f (x) denote the set of all subgradients of f (x) at x. We call ∂f (x) the “subdifferential of f (x).” • If f (x) is convex, then ∂f (x) is always a nonempty convex set. • If f (x) is differentiable, then ∂f (x) = {∇f (x)}. • Subgradients plays the same role for convex functions as the gradient does for differentiable functions. Consider the following optimization problem: min f (x) x – If f (x) is convex and differentiable, then x is a global minimum if and only if ∇f (x) = 0. – If f (x) is convex and non-differentiable, then x is a global minimum if and only if 0 ∈ ∂f (x). 4.3 Subgradients for Concave Functions If f (x) is a concave function, then g is a subgradient of f (x) at x = x̄ if: ¯ f (x) ≤ f (¯ x) + g T (x − x) for all x ∈ X . This is illustrated in Figure 7. Figure 8 shows a piecewise-linear concave function. Figure 9 illustrates the subdifferential for a concave function. 4.4 Computing Subgradients Subgradients play a very important role in non-differentiable optimization. In most algorithms, we assume that we have a subroutine that receives as input a value x, and has output g where g is a subgradient of f (x). 18 3000 2000 1000 0 f(u) −1000 −2000 −3000 −4000 −5000 −100 −80 −60 −40 −20 0 20 40 u Figure 7: The subgradient of a concave function. 19 60 80 100 500 400 300 200 f(u) 100 0 −100 −200 0 20 40 60 80 100 120 140 u Figure 8: A piecewise linear concave function. 20 160 180 200 200 150 100 50 f(u) 0 −50 −100 −150 0 20 40 60 80 100 120 140 160 u Figure 9: The subdifferential of a concave function. 21 180 200 5 The Subgradient Method for Maximizing a Con­ cave Function Suppose that Z(u) is a concave function, and that we seek to solve: P : maximizeu Z(u) u ∈ n . s.t. If Z(u) is differentiable and d := ∇Z(¯ u) satisfies d = 0, then d is an ascent direction at ū, namely Z(¯ u +
d) > Z(¯ u) for all
> 0 and sufficiently small . This is illustrated in Figure 10. However, if Z(u) is not differentiable and g is a subgradient of Z(u) at u = ū, then g is not necessarily an ascent direction. This is illustrated in Figure 11. The following algorithm generalizes the steepest descent algorithm and can be used to maximize a nondifferentiable concave function Z(u). Step 0: Initialization. Start with any point u1 ∈ n . Choose an infinite sequence of positive stepsize values {αk }∞ k=1 . Set k = 1. Step 1: Compute a subgradient. Compute g ∈ ∂Z(uk ). Step 2: Compute stepsize. Compute stepsize αk from stepsize series. g Step 3: Update Iterate. Set uk+1 ← uk + αk g . Set k ← k + 1 and go to Step 1. As it turns out, the viability of the subgradient algorithm depends critically on the sequence of stepsizes: Property: Suppose that {αk }∞ k=1 satisfies: ∞ lim αk = 0 k→∞ and αk = ∞ . k=1 Then under very mild additional assumptions, sup Z(uk ) = maxn Z(u) . k u∈ 22 g _ f(x) + Figure 10: The gradient is an ascent direction. 23 _ f(x)T(x-x) g _ _ f(x) + gT(x-x) Figure 11: A subgradient is not necessarily an ascent direction. 24 5.1 Example of Subgradient Algorithm in One Variable Consider the following concave optimization problem: P : maximizeu Z(u) = min{0.5u + 2, −1u + 20} s.t. u∈. We illustrate various implementations of the subgradient method on this simple problem. • Choose u1 = 0 and αk = 0.k14 . Figure 12 illustrates the performance of the subgradient algorithm for this stepsize sequence. • Choose u1 = 0 and αk = 0.02. Figure 13 illustrates the performance of the subgradient algorithm for this stepsize sequence. • Choose u1 = 0 and αk = 0.k01 . Figure 14 illustrates the performance of the subgradient algorithm for this stepsize sequence. • Choose u1 = 0 and αk = 0.01 × (0.9)k . Figure 15 illustrates the performance of the subgradient algorithm for this stepsize sequence. 5.2 Example of Subgradient Algorithm in Two Variables Consider the following concave optimization problem: P : maximizeu Z(u) = min{ 2.8571u1 − 0.2857u2 − 5.7143, −u1 + u2 + 2, −0.1290u1 − 1.0323u2 + 21.1613} s.t. u ∈ 2 . We illustrate the implementation of the subgradient method on this problem with u1 = (0, 0) and αk = √1k . Figure 16 shows the function level sets and the path of iterations. Figure 17 shows the objective function values, and Figure 18 shows values of the variables u = (u1 , u2 ). 25 8 7.99 7.98 Z(u)7.97 7.96 7.95 7.94 7.93 11.85 11.9 11.95 12 12.05 12.1 u 12.08 8.06 12.06 8.04 8.02 k Z(u ) 12.04 k u 12.02 8 7.98 12 7.96 11.98 7.94 11.96 7.92 0 5 10 15 0 5 k 10 k Figure 12: Illustration of subgradient algorithm, αk = 26 0.14 k . 15 8 7.99 7.98 Z(u)7.97 7.96 7.95 7.94 7.93 11.85 11.9 11.95 12 12.05 12.1 u 8.08 8.06 12.04 8.04 12.02 k8.02 k Z(u ) u 12 8 11.98 7.98 11.96 7.96 0 5 10 15 k 0 5 10 k Figure 13: Illustration of subgradient algorithm, αk = 0.02 . 27 15 8 7.99 7.98 Z(u)7.97 7.96 7.95 7.94 7.93 11.85 11.9 11.95 12 12.05 12.1 u 8.08 12.04 8.06 12.02 8.04 Z(uk) uk 8.02 12 11.98 8 11.96 7.98 0 5 10 15 0 5 k 10 k Figure 14: Illustration of subgradient algorithm, αk = 28 0.01 k . 15 8 7.99 7.98 Z(u)7.97 7.96 7.95 7.94 7.93 11.85 11.9 11.95 12 12.05 12.1 u 8.08 12.04 8.06 12.02 8.04 Z(uk) uk 8.02 12 11.98 8 11.96 7.98 0 5 10 15 k 0 5 10 k Figure 15: Illustration of subgradient algorithm, αk = 0.01 × (0.9)k . 29 15 0 2 4 6 8 u2 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 u 1 Figure 16: Illustration of the subgradient method in two variables: level sets and path of iterations. 30 20 8 7 6 5 4 Z(uk) 3 2 1 0 −1 −2 0 10 20 30 40 50 k Figure 17: Illustration of the subgradient method in two variables: objective function values. 31 60 8 6 uk 14 2 0 0 10 20 30 40 50 60 40 50 60 k 14 12 10 uk 2 8 6 4 2 0 0 10 20 30 k Figure 18: Illustration of the subgradient method in two variables: values of variables u = (u1 , u2 ). 32 6 Solution of the Lagrangian Dual via Subgradient Optimization We start with the primal problem: OP : minimumx s.t. f (x) gi (x) ≤ 0, i = 1, . . . , m x ∈ P, We create the Lagrangian: L(x, u) := f (x) + uT g(x) and the dual function: L∗ (u) := minimumx∈P f (x) + uT g(x) The dual problem then is: D : maximumu L∗ (u) s.t. u≥0 Recall that L∗ (u) is a concave function. The premise of Lagrangian duality is that it is “easy” to compute L∗ (¯ u) for any given u. ¯ That is, it is easy to compute an optimal solution x̄ of u) := minimumx∈P L∗ (¯ f (x) + u ¯T g(x) = f (¯ x) + u ¯T g(¯ x) for any given u, ¯ where x ¯ ∈ P . It turns out that computing subgradients of L∗ (u) is then also easy. We have: 33 Property: Suppose that u ¯ is given and that x ¯ ∈ P is an optimal solution of L∗ (¯ u) = min f (x) + u ¯T g(x). Then g := g(¯ x) is a subgradient of L∗ (u) at x∈P u = ū. Proof: For any u ≥ 0 we have L∗ (u) = min f (x) + uT g(x) x∈P ≤ f (¯ x) + uT g(¯ x) = x) + (u − u) ¯ T g(¯ x) f (¯ x) + u ¯T g(¯ = min f (x) + u ¯T g(x) + g(¯ ¯ x)T (u − u) = L (¯ u) + g T (u − u) ¯ . x∈P ∗ Therefore g is a subgradient of L∗ (u) at ū. q.e.d. The subgradient method for solving the Lagrangian dual can now be stated: Step 0: Initialization. Start with any point u1 ∈ n , u1 ≥ 0. Choose an infinite sequence of positive stepsize values {αk }∞ k=1 . Set k = 1. Step 1: Compute a subgradient. Solve for an optimal solution x̄ of L∗ (uk ) = min f (x) + (uk )T g(x). Set g := g(x̄). x∈P Step 2: Compute stepsize. Compute stepsize αk from stepsize series. g . If uk+1 ≥ 0, re-set Step 3: Update Iterate. Set uk+1 ← uk + αk g ← max{uk+1 , 0}, i = 1, . . . , m. Set k ← k + 1 and go to Step 1. uk+1 i i Note that we have modified Step 3 slightly in order to ensure that the values of uk remain nonnegative. 6.1 Illustration and Exercise using the Subgradient Method for solving the Lagrangian Dual Consider the primal problem: 34 cT x OP : minimumx Ax − b ≤ 0 s.t. x ∈ {0, 1}n . Here g(x) = Ax − b and P = {0, 1}n = {x | xj = 0 or 1, j = 1, . . . , n}. We create the Lagrangian: L(x, u) := cT x + uT (Ax − b) and the dual function: L∗ (u) := minimumx∈{0,1}n cT x + uT (Ax − b) The dual problem then is: D : maximumu L∗ (u) s.t. u≥0 u) is: Now let us choose u ¯ ≥ 0. Notice that an optimal solution x ¯ of L∗ (¯ x ¯j = 0 if (c − AT ū)j ≥ 0 1 if (c − AT ū)j ≤ 0 for j = 1, . . . , n. Also, n ¯T (A¯ x − b) = −u u) = cT x ¯+u ¯T b − L∗ (¯ (c − AT u) ¯j j=1 Also g := g(¯ x−b x) = A¯ 35 − . is a subgradient of L∗ (ū). Now consider the following data instance of this problem:     7 −8 12  −2 −2      −1       5  , b =  45  A= 6      −5  20  6  3 12 42 and cT = ( −4 1 ) . Solve the Lagrange dual problem of this instance using the subgradient algorithm starting at u1 = (1, 1, 1, 1, 1)T , with the following step-size choices: • αk = 1 k • αk = √1 k for k = 1, . . .. for k = 1, . . .. • αk = 0.2 × (0.75)k for k = 1, . . .. • a stepsize rule of your own. 36