Solution Methods for Quadratic Optimization 1 1 Outline • Active Set Methods for Quadratic Optimization • Pivoting Algorithms for Quadratic Optimization • Four Key Ideas that Motivate Interior Point Methods • Interior-Point Methods for Quadratic Optimization • Reduced Gradient Algorithm for Quadratic Optimization • Some Computational Results 2 Active Set Methods for Quadratic Optimization In a constrained optimization problem, some constraints will be inactive at the optimal solution, and so can be ignored, and some constraints will be active at the optimal solution. If we knew which constraints were in which category, we could throw away the inactive constraints, and we could reduce the dimension of the problem by maintaining the active constraints at equality. Of course, in practice, we typically do not know which constraints might be active or inactive at the optimal solution. Active set methods are designed to make an intelligent guess of the active set of constraints, and to modify this guess at each iteration. Herein we describe a relatively simple active-set method that can be used to solve quadratic optimization problems. Consider a quadratic optimization problem in the format: QP : 1 T 2 x Qx minimizex s.t. + cT x Ax = b x≥0. 2 (1) The KKT optimality conditons for QP are as follows: Ax = b x≥0 Qx + c − AT p − s = 0 s≥0 xj sj = 0, j = 1, . . . , n. We suppose for this problem that n is very large, but that intuition suggests that very few variables xj will be non-zero at the optimal solution. Suppose that we have a feasible solution x, and that we partition the components of x into x = (xβ , xη ), where xβ ≥ 0 and xη = 0. In this partition, the number of elements in β will be small, relative to the value of n. Using this partition, we can re-write the data for our problem as: Q= Qββ Qβη Qηβ Qηη , A = [Aβ Aη ] , c= cβ cη . Now let us “guess” that the variables xj , j ∈ η will be zero in the optimal solution. That being the case, we can eliminate the variables xη from further consideration for the moment. We will then concentrate our eﬀorts on solving the much smaller problem: QPβ : zβ∗ = minimumxβ s.t. 1 T 2 xβ Qββ xβ + cTβ xβ Aβ xβ = b xβ ≥ 0 . Let x∗β be the optimal solution of QPβ , and let p∗ be the associated KKT multipliers for the constraints “Aβ xβ = b”. Also, let sβ denote the 3 associated KKT multipliers on the nonnegativity constraints xβ ≥ 0. Then x∗β and p∗ , sβ will satisfy the KKT optimality conditions for QPβ , namely: Qββ x∗β + cβ − ATβ p∗ − sβ = 0 Aβ x∗β = b x∗β ≥ 0, sβ ≥ 0 (x∗β )T sβ = 0 . We can expand x∗β into the following feasible solution to the original problem: x = (xβ , xη ) = (x∗β , 0) by setting all of the variables xη := 0. This solution will of course be feasible for the original problem, but will it be optimal? The answer lies, of course, in checking the KKT optimality conditions for the original problem using x = (x∗β , 0) and p = p∗ , which we write as follows: Qββ Qβη Qηβ Qηη [Aβ Aη ] x∗β 0 xβ∗ 0 + cβ cη − ATβ ATη p∗ = sβ sη ≥0 =b (x∗β , 0) ≥ 0 (x∗β )T sβ = 0, 0T sη = 0 . Here we see that all of these conditions are automatically satisﬁed except for the condition that sη ≥ 0. We remark that sη is deﬁned above to be: sη := Qηβ x∗β + cη − ATη p∗ . We therefore have: 4 Proposition 2.1 If sη ≥ 0, then (x∗β , 0) is an optimal solution to (1). Proof: Since x∗β is an optimal solution of QPβ , then the KKT multipliers p∗ and sβ satisfy Qββ x∗β + cβ − ATβ p∗ − sβ = 0 Aβ x∗β = b x∗β ≥ 0, sβ ≥ 0 (x∗β )T sβ = 0 . If furthermore 0 ≤ sη := Qηβ x∗β + cη − ATη p∗ we have that the vector (xβ∗ , 0) satisﬁes Qββ Qβη Qηβ Qηη [Aβ Aη ] x∗β 0 x∗β 0 + cβ cη − ATβ ATη p∗ = sβ sη ≥0 =b (x∗β , 0) ≥ 0 (x∗β )T sβ = 0, 0T sη = 0 . which are precisely the optimality conditions for (1). Just as importantly, if sj < 0 for some j ∈ η, then the optimal solution to the expanded problem QPβ∪{j} will have a strictly better objective function value than that of QPβ . Proposition 2.2 If sj < 0 for some j ∈ η, then zβ∪{j} < zβ . We defer the proof of this proposition to the end of this section. This all suggests the following active-set algorithm for solving our original problem QP: 5 0. Deﬁne some partition of the indices, {1, . . . , n} = β ∪ η. 1. Solve the reduced problem QPβ : zβ∗ = minimumxβ 1 T 2 xβ Qββ xβ s.t. + cβT xβ Aβ xβ = b xβ ≥ 0 , obtaining the optimal solution x∗β and optimal KKT multiplier p∗ on the equations “Aβ xβ = b”. 2. Compute sη := Qηβ x∗β + cη − ATη p∗ . If sη ≥ 0, STOP. In this case, (x∗β , 0) is optimal for (1). Otherwise if sj < 0 for some j ∈ η, then set β ← β ∪ {j}, η := η \ {j}, and go to 1. In practice we might add more than one index j to the set β at each iteration. Some active set methods add all indices j for which sj < 0. Also, in order to keep the set β from becoming too large, we might remove indices i from β if xi has been zero in the optimal solutions of the reduced problems for a large number of previous iterations. There are also very many other variants of this type of method. Proof of Proposition 2.2: For ease of notation we will drop the subscript β from the subproblem solved in step 1 of the algorithm. Let us call x̄ the optimal solution of this problem, together with KKT multipliers p̄ and s̄. These all satisfy the KKT conditions: Qx̄ + c − AT p̄ − s̄ = 0 Ax̄ = b x ¯ ≥ 0, s̄ ≥ 0, x ¯T s̄ = 0 . ¯2 = 0 and We separate the indices of x into sets such that x ¯1 > 0 and x we separate the data for the problem as follows: 6 Q= Q11 Q12 Q21 Q22 , c= c1 c2 , A = [A1 A2 ] . We then write the above system as: ¯1 + c1 − AT1 p̄ = 0 Q11 x ¯1 + c2 − AT2 p̄ − s̄2 = 0 Q21 x ¯1 = b A1 x s̄2 ≥ 0 . x̄1 > 0, We now consider a new variable θ to add to this system. The data associated with this variable is q, δ, ĉ, a, such that: Q→ Q qT q δ , c→ c ĉ A → [A a] , Since this variable by supposition has a “reduced cost” that is negative, then the new data satisﬁes: ¯ + cˆ − aT p̄ = −ε < 0 . qT x We will now show that this implies that there exists a feasible descent direction. The objective function of our problem and the formulas for the gradients of the objective function are: f (x, θ) = 1 T (x , θ) 2 ∇f (x, θ) = x, 0) = ∇f (¯ Q qT q δ Q qT Q¯ x+c qT x ¯ + cˆ q δ x θ +  x θ c ĉ T + (c , ĉ)  Q11 x ¯ 1 + c1   ¯ 1 + c2  . =  Q21 x qT x ¯ + ĉ 7 x θ Invoking a nondegeneracy assumption, we assume that there exists a vector d1 for which A1 d1 + a = 0. Then, since x̄1 > 0 we have that (x1 , x2 , θ) = (x̄1 + αd1 , 0, α) ≥ 0, for α > 0 and small. Next we compute:   Q11 x ¯ 1 + c1   ¯ 1 + c2  ∇f (¯ x, 0)T (d1 , 0, 1) = (dT1 , 0, 1)  Q21 x qT x ¯ + cˆ = (AT1 p̄)T d1 + q T x ¯ + cˆ ¯ + cˆ = −ε < 0 . = −aT p̄ + q T x Therefore (d1 , 0, 1) is a feasible descent direction, and so the optimal objective function value of the expanded problem is strictly less than that of the current problem. 3 Pivoting Algorithms for Quadratic Optimization To facilitate the presentation, we will consider the following format for QP: QP : minimizex 1 T 2 x Qx s.t. + cT x Ax ≥ b (π) x≥0 (y) Let us ﬁrst look at the KKT optimality conditions for (2): 8 (2) s := Ax − b ≥ 0 x≥0 Qx + c − AT π − y = 0 π≥0 y≥0 πi si = 0 i = 1, . . . , m yj xj = 0 j = 1, . . . , n. We can write the ﬁrst ﬁve conditions in the form of a tableau of linear constraints in the nonnegative variables (y, s, x, π): y s I 0 0 I x π −Q AT −A 0 Relation = = RHS c −b and keeping in mind the nonnegativity conditions: y≥0 s≥0 x≥0 π≥0 and the complementarity conditions: yj xj = 0, j = 1, . . . , n, si πi = 0, i = 1, . . . , m. This tableau has exactly (n + m) rows and 2(n + m) columns, that is, there are twice as many columns as rows. Therefore a basic feasible solution of this tableau will have one half of all columns in the basis. We can ﬁnd a basic feasible solution to this tableau (which has no objective function row by the way!) by applying Phase I of the simplex algorithm. However, what we really want to do is to ﬁnd basic feasible solution that is also “complementary” in the following sense: 9 • Either yj is in the basis or xj is in the basis, but not both, j = 1, . . . , n, • Either si is in the basis or πi is in the basis, but not both, i = 1, . . . , m. Such a basic feasible solution will have the nice property that either yj = 0 or xj = 0, j = 1, . . . , n, and also si = 0 or πi = 0, i = 1, . . . , m. That is, a complementary basic feasible solution (CBFS) will satisfy the KKT conditions and so will solve QP. One strategy for ﬁnding a CBFS is to add a nonnegative parameter θ and solve the following tableau: y s I 0 0 I x π −Q AT −A 0 Relation = = RHS c + θe −b + θe (Here e is the vector of ones of any appropriate dimension.) Now notice that for θ suﬃciently large, that the RHS is nonnegative and so the complementary basic solution:     y c + θe  s   −b + θe   =  x   0 π 0 is a complementary basic feasible solution. The strategy then is try to start from this solution, and to drive θ to zero by a sequence of pivots that maintain a complementary basic feasible solution all the while. As it turns out, it is easy to derive pivot rules for such a strategy. The resulting algorithm is called the “complementary pivot algorithm”. It was originally developed by George Dantzig and Richard Cottle in the 1960s, and analyzed and extended by Carl Lemke, also in the 1960s. We will not go over the speciﬁc pivot rules of the complementary pivot algorithm. They are simple but technical, and are easy to implement. What can we say about convergence of the complementary pivot algorithm? Just like the simplex method for linear optimzation, one can prove that the 10 complementary pivot algorithm must converge in a ﬁnite number of pivots (whenever the matrix Q is SPSD). Regarding computational performance, notice that the tableau generally has twice the number of constraints and twice the number of variables as a regular linear optimization problem would have. This is because the primal and the dual variables and constraints are written down simultaneously in the tableau. This leads to the general rule that a complementary pivot algorithm will take roughly eight times as long to solve a quadratic program via pivoting as it would take the simplex method to solve a linear optimization problem of the same size as the primal problem alone. Up until recently, most software for quadratic optimization was based on the complementary pivot algorithm. However, nowadays most software for quadratic optimization uses interior-point methods, which we will soon discuss. 4 Four Key Ideas that Motivate Interior Point Methods We will work in this section with quadratic optimization in the following format: QP : 1 T 2 x Qx minimizex s.t. + cT x Ax = b (3) x≥0. 4.1 A barrier function The feasible region of (QP) is the set F := {x | Ax = b, x ≥ 0} . The boundary of F will be where one or more of the components of x is equal to zero. The boundary is composed of linear pieces, and is generally 11 arranged in an exponentially complex manner. The simplex method takes explicit advantage of the linearity of the boundary, but can take a long time to solve a problem if there are exponentially many extreme points and edges to traverse. For this reason, we will take the curious step of introducing a barrier function B(x). A barrier function is any (convex) function deﬁned on the interior of F and whose values go to inﬁnity as the boundary is approached. That is: • B(x) < +∞ for all x ∈ intF, and • B(x) → ∞ as xj → 0. The idea in a barrier method is to dissuade points x from ever approaching the boundary of the feasible region. For a ﬁxed value of µ > 0, we consider solving: QPµ : minimizex 1 T 2 x Qx s.t. + cT x + µB(x) Ax = b (4) x>0. Let x(µ) denote the optimal solution of QPµ . For large values of µ, the optimal solution x(µ) will be more concerned about staying away from the boundary than in optimizing the original quadratic objective function. However, for small values of µ, the optimal solution x(µ) will be ever more concerned about optimizing the original quadratic objective function while simultaneously not getting too close to the boundary. Interior-point methods use the barrier function: B(x) = − n ln(xj ) . j=1 This function is convex, and approaches +∞ as any xj approaches zero. 4.2 Staying Centered relative to the Objective Function Value Suppose that x̄ is a the current iterate of an algorithm for solving (QP). If x̄ is very far from the boundary of the feasible region F, then no matter 12 what direction we wish to go at the next iteration, we will have plenty of room to take a step, see Figure 1. However, since the optimal solution will generally lie on the boundary of the feasible region, we will want to keep the algorithm iterates as centered as possible for their objective function values. In the ﬁgure, point a is close to the boundary, point b is far from the boundary, and points c and d are far from the boundary relative to their objective function level sets. Now consider the barrier version of our quadratic problem with the logarithmic barrier: BPµ : minimizex s.t. 1 T 2 x Qx + cT x − µ n j=1 ln(xj ) Ax = b x>0. The objective function here will minimize the barrier term while also accounting for the original quadratic objective function. In this way, the solutions x(µ) of BPµ will be nicely centered relative to other feasible values x ∈ F with as good or better objective function values. Therefore the solutions x(µ) will be relatively far from the boundary for their objective function values. 4.3 Treating the KKT Conditions of the Barrier Problem as a Nonlinear Equation System For a general barrier function B(x), our parameterized barrier problem is: QPµ : minimizex 1 T 2 x Qx s.t. + cT x + µB(x) Ax = b (5) x>0. Because the positivity conditions x > 0 will never be binding, the optimality conditions for this problem are: 13 Objective d c a b Figure 1: Point a is close to the boundary, point b is far from the boundary. Points c and d are far from the boundary relative to their objective function level sets. 14 Ax = b Qx + c − µ∇B(x) − AT p = 0 x>0. This is a nonlinear equation system with n + m equations in the n + m variables (x, p). A generic algorithm for this system can be described as follows: 0. Start with some x ¯ > 0 and some p. ¯ 1. If A¯ x − b ≤ and Q¯ x + c − µ∇B(¯ x) − AT p̄ ≤ for some small tolerance > 0, STOP. Otherwise go to Step 2. 2. Use some method, perhaps Newton’s method, to generate a desired step d = (dx , dp ). 3. Compute a step-size that will maintain the positivity of the iterate value of x. Set r = 0.99, and compute: θ = min 1, r min ∆xj <0 x ¯j −∆xj . 4. Re-set (¯ x, p) ¯ ← (¯ x + θdx , p̄ + θdp ), and go to Step 1. In this method, we check to see if we satisfy our nonlinear equations (to within some tolerance) in Step 1. We generate a presumably good direction x, p) ¯ ← (¯ x + dx , p̄ + d = (dx , dp ) in Step 2. Although we would like to re-set (¯ dp ), we must ensure that our iterate values of x remain strictly positive. This is accomplished in Step 3. Note that the value r = 0.99 could be set to any positive value strictly less than 1.0 . 15 4.4 Newton’s Method The last motivating idea in interior-point methods is Newton’s method. In its various forms, Newton’s method is a method for solving nonlinear equations by linearizing the equations at each iteration. Suppose that we wish to solve the equation system: h1 (x) = 0 .. . .. . .. . (6) hn (x) = 0 where x = (x1 , . . . , xn ). We can write this succinctly as: h(x) = 0 where h(x) = (h1 (x), . . . , hn (x)). Let x ¯ denote the current iterate value. Then from Taylor’s theorem, we have: ˜ x + ∆x) := h(¯ h(¯ x + ∆x) ≈ h(¯ x) + [H(¯ x)] ∆x , (7) where H(x) is the Jacobian matrix: Hij (x) := ∂hi (x) . ∂xj ˜ x + ∆x) is the linear approximation of h(x) around the Here the function h(¯ point x = x. ¯ We can try to solve the system h(¯ x + ∆x) = 0 by instead ˜ solving h(¯ x + ∆x) = 0, which from (7) is: h(¯ x) + [H(¯ x)] ∆x = 0 . (8) Solving (8) for ∆x, we obtain: x) . d := ∆x = − [H(¯ x)]−1 h(¯ (9) The direction d is called the Newton direction or the Newton step. In a pure version of Newton’s method, we would set the next iteration value to be ¯+d . xnew = x 16 Newton’s method has remarkably fast convergence properties if the iterates lie near to the solution of the system. However, if the iterates are very far away from the solution of the system, Newton’s method typically does not have particularly good convergence properties. 5 5.1 Interior-Point Methods for Quadratic Optimization The Central Path Equation System For this section it will be most convenient to work with quadratic optimization in the following format: QP : minimizex 1 T 2 x Qx s.t. + cT x Ax = b (10) x≥0. Using the usual Lagrange duality constructions, we can derive the following dual problem: DP : maximizep,s,x bT p − 12 xT Qx AT p + s − Qx = c s.t. s≥0. The following weak duality result follows by direct algebraic manipulation of the feasibility conditions: Proposition 5.1 If x if feasible for (QP) and (p, s, x) is feasible for (DP), then the duality gap is: 1 1 T x Qx + cT x − bT p − xT Qx = sT x . 2 2 17 The KKT optimality conditons for QP are as follows: Ax = b x≥0 AT p + s − Qx = c s≥0 xj sj = 0, j = 1, . . . , n. The strong duality result is as follows: Proposition 5.2 If x if feasible for (QP) and (p, s, x) is feasible for (DP), then x if optimal for (QP) and (p, s, x) is optimal for (DP) if and only if: sT x = 0 . Before we move on, let us introduce some notation that is used in interiorpoint methods. For vectors x and s, we denote:   X= x1 0 0 .. .      S= , xn s1 0 .. . 0 Therefore, for example, Se = s,   1   , . e =  ..  . 1 sn  Xe = x,   x1 s1  ..  XSe =  .  . xn sn Now let us relax the KKT conditions above, by picking a value of µ > 0 and instead solve: Ax = b, x > 0 AT p + s − Qx = c, s>0 xj sj = µ, j = 1, . . . , n. 18 With our interior-point notation we can write this system as: Ax = b, x>0 AT p + s − Qx = c, s>0 XSe = µe . A solution (x, p, s) = (xµ , pµ , sµ ) to this system for some value µ is called a point on the central path or central trajectory of the problem. Proposition 5.3 If (x, p, s) is on the central path for some value µ > 0, then the duality gap is 1 T 1 x Qx + cT x − bT p − xT Qx = sT x = nµ . 2 2 The proof of this proposition is a straightforward extension of the algebra of the system. The value of this proposition is that it shows that if we solve the system for smaller values of µ, then our duality gap will go to zero, and we will approximately solve our original problem. We can also derive the central path by considering the following logarithmic barrier problem associated with QP for a given value of µ: BPµ : minimizex s.t. 1 T 2 x Qx + cT x − µ n j=1 ln(xj ) Ax = b x>0. The KKT conditions for this logarithmic barrier problem are: Ax = b, x>0 Qx + c − µX −1 e − AT p = 0 . 19 If we then deﬁne s := µX −1 e > 0 and rearrange the above system, we obtain: Ax = b, x>0 AT p + s − Qx = c XSe = µe . Since x > 0 implies that s > 0, this is exactly the same system of equations we had derived earlier that deﬁne the central path. 5.2 The Newton Equation System Suppose that we have a current “point” (actually a triplet) (¯ x, p, ¯ s) ¯ for which ¯ > 0 and s¯ > 0, and we are given a value of µ > 0. We want to compute a x direction (∆x, ∆p, ∆s) so that (¯ x + ∆x, p̄ + ∆p, s̄ + ∆s) approximately solves the central path equations. We set up the system: (1) A(x̄ + ∆x) = b p + ∆p) + (¯ s + ∆s) = c (2) −Q(¯ x + ∆x) + AT (¯ ¯ + ∆X)(S¯ + ∆S)e = µe . (3) (X If we linearize this system of equations and rearrange the terms, we obtain the Newton equations for the current point (¯ ¯ x, p, ¯ s): (1) A∆x = b − Ax̄ =: r1 (2) −Q∆x + AT ∆p + ∆s = c + Qx̄ − AT p̄ − s̄ =: r2 ¯ ¯ ¯ =: r3 (3) S̄∆x + X∆s = µe − XSe We refer to the solution (∆x, ∆p, ∆s) to the above system as the Newton direction at the point (¯ ¯ x, p, ¯ s). 20 5.3 The Step-sizes Given our current point (¯ x, p, ¯ s) ¯ and a given value of µ, we compute the Newton direction (∆x, ∆p, ∆s) and we update our variables by choosing a step-size θ, to obtain new values: (˜ x, p, ˜ s) ˜ ← (¯ x, p, ¯ s) ¯ + θ(∆x, ∆p, ∆s) . In order to ensure that x ˜ > 0 and s˜ > 0, we choose a value of r satisfying 0 < r < 1 (r = 0.99 is a common value in practice), and determine θ as follows: θP = min 1, r min ∆xj <0 θD = min 1, r min ∆sj <0 x ¯j −∆xj s̄j −∆sj . We then set the step-size to be: θ := min{θP , θD } . 5.4 Decreasing the Path Parameter µ We also want to shrink µ at each iteration, in order to (hopefully) shrink the duality gap. The current iterate is (¯ ¯ and the current values satisfy: x, p, ¯ s), µ≈ We then re-set µ to µ← 1 10 x̄T s¯ n T x ¯ s̄ n where “10” is user-speciﬁed. 21 , 5.5 The Stopping Criterion We typically declare the problem solved when it is “almost” primal feasible, “almost” dual feasible, and there is “almost” no duality gap. We set our tolerance to be a small positive number, typically = 10−8 , for example, and we stop when: (1) Ax̄ − b ≤ (2) x + AT p̄ + s¯ − c ≤ − Q¯ ¯≤ (3) s̄T x 5.6 The Full Interior-Point Algorithm a) Given (x0 , p0 , s0 ) satisfying x0 > 0, s0 > 0, and µ0 > 0, and r satisfying 0 < r < 1, and > 0. Set k ← 0. b) Test stopping criterion. Check if: (1) Axk − b ≤ (2) − Qxk + AT pk + sk − c ≤ (3) (sk )T xk ≤ . If so, STOP. If not, proceed. c) Set µ ← 1 10 (xk )T (sk ) n d) Solve the Newton equation system: (1) A∆x = b − Axk =: r1 (2) −Q∆x + AT ∆p + ∆s = c + Qxk − AT pk − sk =: r2 (3) S k ∆x + X k ∆s = µe − X k S k e =: r3 22 e) Determine the step-size: θP = min 1, r min ∆xj <0 θD = min 1, r min ∆sj <0 xkj −∆xj skj −∆sj θ = min{θP , θD } . f) Update values: (xk+1 , pk+1 , sk+1 ) ← (xk , pk , sk ) + θ(∆x, ∆p, ∆s) . Re-set k ← k + 1 and return to (b). 5.7 Solving the Newton Equations We now show how to manipulate the Newton equations to yield some formulas for the solution of the Newton equation system. Given our current x, p, ¯ s) ¯ > 0 and s¯ > 0, and a value of µ > 0, we wish to point (¯ ¯ for which x solve the Newton equations below for (∆x, ∆p, ∆s): (1) A∆x = b − Ax̄ =: r1 (2) −Q∆x + AT ∆p + ∆s = c + Qx̄ − AT p̄ − s̄ =: r2 ¯ ¯ ¯ =: r3 (3) S̄∆x + X∆s = µe − XSe ¯ we obtain the following string of Pre-multiplying equation (2) by X, equations that derive from one another: 23 ¯ ¯ T ∆p + X∆s ¯ ¯ 2 −XQ∆x + XA = Xr ¯ ¯ T ∆p − S̄∆x + r3 = Xr ¯ 2 −XQ∆x + XA ¯ T ∆p = Xr ¯ 2 − r3 ¯ + S¯ ∆x + XA − XQ ¯ + S¯ −∆x + XQ −1 ¯ + S¯ −r1 + A XQ ¯ + S¯ A XQ −1 ¯ T ∆p = XQ ¯ + S¯ XA −1 −1 ¯ T ∆p = A XQ ¯ + S¯ XA ¯ T ∆p = A XQ ¯ + S¯ XA −1 ¯ 2 − r3 Xr −1 ¯ 2 − r3 Xr ¯ 2 − r 3 + r1 Xr Therefore ∆p is determined by solving the last system above, or equivalently, by the formula: ¯ + S¯ ∆p = A XQ −1 ¯ T XA −1 ¯ + S¯ A XQ −1 ¯ 2 − r 3 + r1 Xr . Once ∆p has been computed, we can compute ∆x and ∆s as follows: ¯ + S¯ ∆x = XQ −1 ¯ 2 + XA ¯ T ∆p r3 − Xr ∆s = Q∆x + r2 − AT ∆p . Note that if Q is very sparse or nicely structured (for example, if Q is diagonal), then the bulk of the work above lies in forming the matrix: ¯ + S¯ A XQ −1 ¯ T XA and then solving a system of equations with this matrix. In particular when Q is diagonal, then the sparsity pattern of this matrix will be the same as the sparsity pattern of AAT . 24 5.8 Remarks on Interior-Point Methods • The algorithm just described is almost exactly what is used in commercial interior-point method software. • A typical interior-point code will solve a linear or quadratic optimization problem in 25-80 iterations, regardless of the dimension of the problem. • There is a very “deep” theory of convergence for interior-point methods, when the constant “10” is replaced by the number √ 2+4 n √ 1+4 n • Interior-point methods are still an area of very active research. • These days, interior-point methods have been extended to allow for the solution of problems with convex quadratic constraints: QCQP : aTi x + 12 xT Qi x ≤ bi , i = 1, . . . , m s.t. 6 cT x + 12 xT Qx minx Reduced Gradient Algorithm for Nonlinear Optimization The reduced gradient algorithm is designed to solve a problem of the form: NLP : minimizex s.t. f (x) Ax = b x≥0. In the case of a suitably formatted quadratic problem, 1 f (x) = xT Qx + cT x 2 25 (11) and ∇f (x) = Qx + c . We will make the following assumptions regarding bases for the system “Ax = b, x ≥ 0”: • For any basis β ⊂ {1, . . . , n} with |β| = m, Aβ is nonsingular. • For any basic feasible solution (xβ , xη ) = (A−1 β b, 0) we have xβ > 0, where here η := {1, . . . , n} \ β. The reduced gradient algorithm tries to mimic certain aspects of the simplex algorithm. At the beginning of an iteration of the reduced gradient algorithm, we have a feasible solution x ¯ for our problem NLP. Pick a basis β for which x̄β > 0. Remember that such a basis must exist by our assumptions, but that such a basis might not be uniquely determined by the current solution x̄. Set η := {1, . . . , n} \ β. Now invoke the equations of our system, which state that: Aβ xβ + Aη xη = b for any x = (xβ , xη ), and so we can write: xβ = A−1 β (b − Aη xη ) . Since the basic variables are uniquely determined by the nonbasic variables above, we can substitute these into the formula for the objective function: f (x) = f (xβ , xη ) = f (A−1 β (b − Aη xη ), xη ) . We therefore can deﬁne our reduced objective function to be: h(xη ) = f (A−1 β (b − Aη xη ), xη ) . The gradient of the reduced objective function is called the reduced gradient, and has the following formula: −1 T −T ∇h(xη ) = ∇fη (A−1 β (b − Aη xη ), xη ) − Aη Aβ ∇fβ (Aβ (b − Aη xη ), xη ) , 26 which evaluates at x̄ to be: x) − ATη A−T x) . ∇h(¯ xη ) = ∇fη (¯ β ∇fβ (¯ Let us re-name the reduced gradient r¯ = r¯η for notational convenience. The reduced gradient can then be used to deﬁne a new direction d as follows: for j ∈ η, we deﬁne: dj :=   −r̄j        −r̄      0    j if x̄j > 0 if r̄j ≤ 0 if x ¯j = 0 and r¯j > 0 For j ∈ β, we deﬁne: dβ := −A−1 β Aη dη . This direction has the following properties, which are fairly easy to prove: • d is a feasible direction: Ad = 0 and x ¯ + αd ≥ 0 for all α > 0 and suﬃciently small. • If d = 0, then x ¯ is a KKT point, with associated KKT multipliers ∇f (x̄) for the equations “Ax = b”. p̄ = A−T β β • If d = 0, then d is a descent direction, that is, ∇f (¯ x)T d < 0. This leads to the following algorithm for solving NLP: 0. The current point is x. ¯ Choose β for which xβ > 0 and |β| = m. Deﬁne η := {1, . . . , n} \ β. 1. Compute the reduced gradient: x) − ATη A−T x) . r̄η := ∇fη (¯ β ∇fβ (¯ 27 2. Compute the direction d. For j ∈ η, deﬁne: dj :=   −r̄j        −r̄      0    j if x̄j > 0 if r̄j ≤ 0 if x ¯j = 0 and r¯j > 0 For j ∈ β, deﬁne: dβ := −A−1 β Aη dη . 3. Compute the step-size. Deﬁne the following three quantities: α1 := max {α|(x̄ + αd)β ≥ 0} α2 := max {α|(x̄ + αd)η ≥ 0} α3 := arg minα f (x̄ + αd)|0 ≤ α ≤ min{α1 , α2 } 4. Update values. If d = 0, then STOP (¯ x is a KKT point). If d = 0, 3 then set x ¯←x ¯ + α d, and go to (0.). 6.1 Some Comments on the Reduced Gradient Algorithm • The reduced gradient algorithm has been very eﬀective on problems with linear equality and inequality constraints. That is because it is actually mimicking the simplex algorithm as much as possible. • One very useful modiﬁcation of the reduced gradient algorithm is to set β to correspond to the index of the set of the m largest components of x̄. • As it turns out, there is no proof of convergence for the algorithm as written above. A slightly modiﬁed rule for determining the direction d is as follows, and leads to a straightforward proof of convergence: 28 For j ∈ η, deﬁne:   −r̄j          −r̄j    dj := 7 if x ¯j > 0 and r¯j ≤ 0 if x ¯j = 0 and r¯j ≤ 0  −¯ rj x ¯j          0    if x ¯j > 0 and r¯j > 0 if x ¯j = 0 and r¯j > 0 Some Computational Results In this section we report some computational results on solving quadratic pattern classiﬁcation problems. The problem under consideration is the dual of the pattern classiﬁcation quadratic program with penalties for the case when we cannot strictly separate the points, namely D1: D1 : maximizeλ,γ s.t. k  k T   k m k m 1 λi + γj −  λi ai − γj bj   λi ai − γ j bj  2 i=1 j=1 i=1 j=1 i=1 j=1 k λi − i=1 m γj = 0 j=1 (12) λi ≤ C i = 1, . . . , k γj ≤ C j = 1, . . . , m λ ≥ 0, γ≥0. We refer to this problem as the REGULAR model. By deﬁning the extra variables x: 29  x :=  k m λi ai − i=1  γj bj  , j=1 we can re-write this model as: k D2 : maximumλ,γ,x i=1 s.t. x− λi + k i=1 k λi m j=1 ai λi − i=1 γj − 12 xT x − m j=1 m γj bj =0 γj = 0 j=1 λi ≤ C i = 1, . . . , k γj ≤ C j = 1, . . . , m λ ≥ 0, γ ≥ 0 λ ∈ k , γ ∈ m . We refer to this model as the EXPANDED model. We solved ten instances of the pattern classiﬁcation model that we created, whose main characteristics are as shown in Table 1. In solving these problem instances, we used the penalty value C = 1, 000. The two-dimensional problems are illustrated in Figure 2, Figure 3, Figure 4, and Figure 5. We solved these instances using the following two algorithms: • IPM (Interior-Point Method) We solved the problems using the software code LOQO, which is an interior-point algorithm implementation for linear, quadratic, and nonlinear optimization. 30 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Figure 2: Separable model with 100 points. 31 0.8 1 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Figure 3: Non-separable model with 100 points. 32 1 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Figure 4: Separable model with 1,000 points. 33 1 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Figure 5: Non-separable model with 1,000 points. 34 1 Separable or Not Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable Dimension n 2 2 2 2 20 20 20 20 20 20 Number of points m 100 100 1,000 1,000 300 300 1,000 1,000 10,000 10,000 Table 1: Ten instances of the pattern classiﬁcation problem. • ACTSET (Active Set Method) We solved the problems using a simple active set method, that adds one new index at each outer iteration, and does not drop any indices from consideration. The method starts with 30 indices in the set β. The interior-point algorithm LOQO was used to solve the subproblems at each iteration. Note that for this particular model, that active set methods should perform very well. This is because we would expect very few points to be binding in the primal, and so very few of the dual variables λi and/or γi will be non-zero in the dual problem. And it is the dual problem which we are going to solve. Table 2 shows the optimal objective function values for each method on each problem. Notice from these tables that both the REGULAR and the EXPANDED models achieve the same optimal objective functions on the diﬀerent algorithms, when they have enough memory to solve the problem. Table 3 shows the number of iterations of each method on each problem. Notice the following: • The non-separable problems take more iterations to solve. 35 Type Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable n 2 2 2 2 20 20 20 20 20 20 m 100 100 1,000 1,000 300 300 1,000 1,000 10,000 10,000 REGULAR MODEL IPM 19.4941 992.4144 1.1398 2,101.6155 EXPANDED MODEL IPM ACTSET 19.4941 19.4941 992.4144 992.4144 88.8831 88.8831 5,751.2161 5,751.2161 1.1398 1.1398 2,101.6155 2,101.6155 1.9110 1.9110 7,024.1791 7,024.1791 4.0165 Table 2: Optimal solution values reported by each method for each problem. Blanks indicate the method ran out of memory. • Interior-point methods use very few iterations in general. • The iteration counts for interior-point methods do not really grow much at all with the size of the problem. • The exception to these observations is the active set algorithm, for which the total iteration count can be misleading. Table 4 shows the CPU time for each method and each problem, in seconds. Notice the following: • The REGULAR model takes more time to solve than does the EXPANDED model. • As the dimensions of the problem grows, the time to solve the models grows, with one exception. • The non-separable problem instances take more time to solve. • The active set method is the clear winner. This leaves several questions: 36 Type Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable n 2 2 2 2 20 20 20 20 20 20 m 100 100 1,000 1,000 300 300 1,000 1,000 10,000 10,000 REGULAR MODEL IPM 27 32 30 40 EXPANDED MODEL IPM ACTSET 26 58 32 194 28 101 60 1,237 32 340 39 778 35 417 54 2,022 598 > 20,000 Table 3: Iteration count for each problem and each method. Blanks indicate the method ran out of memory. Type Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable Separable Non-separable n 2 2 2 2 20 20 20 20 20 20 m 100 100 1,000 1,000 300 300 1,000 1,000 10,000 10,000 REGULAR MODEL IPM 2.89 3.23 27.16 41.39 EXPANDED MODEL IPM ACTSET 0.20 0.04 0.24 0.17 240.85 0.09 278.56 3.39 8.04 0.99 4.30 3.90 175.60 1.27 289.15 12.00 1.02 > 58,759 Table 4: CPU times, in seconds. Blanks indicate the method ran out of memory. 37 • Why does the REGULAR model takes more time to solve than does the EXPANDED model? • Why do the larger models run out of memory? The answer to both of these questions probably has to do with the way LOQO handles the numerical linear algebra of forming and then factoring the Newton equation system at each interior-point iteration. 38