On Covering Methods for D.C. Optimization

Journal of Global Optimization 18: 265–274, 2000. © 2000 Kluwer Academic Publishers. Printed in the Netherlands. 265 On Covering Methods for D.C. Optimization? RAFAEL BLANQUERO and EMILIO CARRIZOSA Departamento de Estadística e Investigación Operativa, Universidad de Sevilla, Sevilla, Spain (e-mail: [email protected]; [email protected]) (Received 21 January 1999; accepted in revised form 19 May 2000) Abstract. Covering methods constitute a broad class of algorithms for solving multivariate Global Optimization problems. In this note we show that, when the objective fis d.c. and a d.c. decomposition for fis known, the computational burden usually suffered by multivariate covering methods is significantly reduced. With this we extend to the (non-differentiable) d.c. case the covering method of Breiman and Cutler, showing that it is a particular case of the standard outer approximation approach. Our computational experience shows that this generalization yields not only more flexibility but also faster convergence than the covering method of Breiman-Cutler. Key words: Covering methods, d.c. Functions, Outer approximation, Power diagrams 1. Introduction Given a polytope X⊂Rmand a real-valued function fdefined on X, we consider the optimization problem max x∈Xf(x) (1) Covering algorithms, [4, 7, 10, 12, 17, 19], constitute an important class of methods for solving problems of type (1) for very large classes of functions f.In essence, covering methods build a feasible sequence {xn}and a sequence of upper envelopes {E(n)}, i.e., functions satisfying E(n) ⩾ffor all n E(n) ⩾E(n+1)for all n E(n)(xj)=f(x j)for all j=1,... ,n If the detection of an ε-optimal solution is used as stopping rule, the prototype covering algorithm is given in Algorithm 1, [16]. ALGORITHM 1. Initialization Take x1∈X Set f=f(x 1),x∗=x1and construct E(1) ?This research has been supported by Grant PB96-1416-C02-02 of DGES, Spain 266 R.BLANQUERO AND E.CARRIZOSA Iteration n = 1,2, ... Set f=maxx∈XE(n)(x) and select xn+1∈argmaxx∈XE(n)(x). If f−f≤εthen Stop. Construct E(n+1)from E(n). If f(x n+1)>f then set f=f(x n+1)and x∗=xn+1 Go to iteration n+1. Whereas the theoretical convergence of this algorithm can be established for general instances from e.g. [16] (Theorem II.2) or [19], (the algorithm finishes with x∗as an ε-optimal solution and fas a lower bound differing less than εfrom the optimal value), care must be taken in order to make these procedure work in practice. Indeed, two critical aspects must be considered, namely, how the upper envelopes E(n) are defined, and then how the subproblems (Qn), max x∈XE(n)(x), (Qn) are solved. The upper envelopes E(n) are usually constructed as pointwise minimum of more simple functions ki, E(n) =min 1≤i≤nki(2) For instance, if fis Lipschitz with known Lipschitz constant L, Piyavskii algorithm, [19, 14, 16] uses kPIY igiven by kPIY i(x) =f(x i)+Lkx−xik Another example of covering method based on (2) is provided in [3, 4, 7, 10]: assuming f:X⊂Rm7−→ Ris continuously differentiable and satisfies the condition f(x)≤f(y)+∇f(y) 0(x −y) +Kkx−yk2∀x,y ∈X(3) where K≥0 is a known constant, one defines kBrC ias kBrC i(x) =f(x i)+∇f(x i)0(x −xi)+Kkx−xik2(4) This is generalized in a later paper, [4], to kBC i(x) =f(x i)+∇f(x i)0(x −xi)+q(x −xi) where q(x) =x0Hx is a quadratic form with positive definite matrix Hand kBC i satisfies, for every choice of xi∈X, the condition kBC i(x) ≥f(x) ∀x∈X ON COVERING METHOD FOR D.C. OPTIMIZATION 267 that is, f(x i)+∇f(x i)0(x −xi)+q(x −xi)≥f(x) ∀xi,x ∈X(5) Note that (3) and kBrC iare particular cases of condition (5) and kBC i, respectively, taking as Hthe diagonal matrix with value Kon the diagonal. Solving the subproblems (Qn)seems to be the bottleneck point of the procedures. In fact, each (Qn)is itself a Global-Optimization (in general non-differentiable) problem, whose resolution seems to be far from trivial. Indeed, it is mentioned in [16] (page 596) that the Piyavskii algorithm seems to be not very promising in dimensions higher than 2 due to the increasing difficulty in solving (Qn).The structure of the problems (Qn)is more simple for those coverings based on (4); however, bounds on the first and second derivatives of fare required, which, as mentioned in [4], is often an equally difficult global optimization problem as the original. In the next section we show that, when a d.c. decomposition of fis known, a new class of upper envelopes can be described. These envelopes will yield subproblems (Qn)with a structure similar to those based on (4), and, at the same time, are more versatile since they can avoid the calculus of the bounds of derivatives and can also be used for non-differentiable problems. 2. Covering algorithm for d.c. functions Let fbe a real-valued d.c. function on X, i.e., fcan be represented as the difference of two functions convex on X, and assume that a d.c. decomposition of fis known, f(x)=f1(x) +f2(x) for each x∈X, where f1is convex on Xand f2is concave on a convex open set ⊃X. We recall the reader that, although for an arbitrary d.c. function f, it is not easy to obtain such a d.c. decomposition, a series of rules enables one to obtain d.c. decompositions of functions fwhich are obtained (via summation, pointwise maximum, linear combinations, ... ) from more simple d.c. functions, [16, 21]. 2.1. THE UPPER ENVELOPE AND PARTICULAR CASES In this case, by majorizing the concave part of the decomposition, we obtain the following convex bounding function: kdc i(x) =f1(x) +(f2(xi)+ξ2(xi)0(x −xi)) (6) where ξ2(xi)is an arbitrary subgradient of f2at xi. In particular, the method of Breiman–Cutler [7] and the generalization due to Baritompa–Cutler [4] are particular cases of our approach (6). Indeed, one has 268 R.BLANQUERO AND E.CARRIZOSA PROPOSITION 2.1. Let f:X⊂Rm→Rsatisfy (5) for a quadratic form qwith symmetric positive semidefinite matrix H. Then, fis a d.c. function and f(x)=x0Hx +(f (x) −x0Hx) (7) is a d.c. decomposition for f. Moreover, if (7) is used as d.c. decomposition of f, then Algorithm 1 with envelopes kdc iyields the algorithm of Baritompa–Cutler. Proof. Since f1(x) =x0Hx is a convex function, it suffices to show that the differentiable function f2(x) =f(x)−x0Hx is concave. By (5), it follows that f(x) ⩽f(y)+∇f(y) 0(x −y) +(x −y)0H(x −y) =f(y)+∇f(y) 0(x −y) +x0Hx −y0Hy −2y0H(x −y) Rearranging terms, f(x)−x0Hx ⩽f(y)−y0Hy +(∇f(y)−2Hy)0(x −y) f2(x) ⩽f2(y) +∇f2(y)0(x −y) , from which we conclude that f2is concave. Moreover, kdc i(x) =f1(x) +f2(xi)+∇f2(xi)0(x −xi) =x0Hx +(f (xi)−x0 iHx i)+(∇f(x i)−2Hx i)0(x −xi) =f(x i)+∇f(x i)0(x −xi)+x0Hx −x0 iHx i−2x0 iH(x −xi) =f(x i)+∇f(x i)0(x −xi)+(x −xi)0H(x −xi) =kBC i(x), showing that kdc yields the envelope of Baritompa–Cutler when (7) is used as d.c. decomposition of f.2 Hence, even for functions fsatisfying (5), different envelopes can be obtained as soon as a d.c. decomposition, other than (7), is provided. This may avoid the costly process of bounding derivatives, [4] and may affect the convergence, as depicted in Example 2.2 and Section 3 (see [6] for further details). EXAMPLE 2.2. In the paper [7], Breiman and Cutler use f:I=[−0.2,1]7−→ Rdefined as f(x)=cos(5πx) −x2 This function is twice continuously differentiable, thus it admits a decomposition of the form f(x)=Kx2−(Kx2−f(x)) (8) In particular, Breiman and Cutler propose such decomposition taking K=12.5π2− 1(the smallest value for which Kx2−f(x)is convex). ON COVERING METHOD FOR D.C. OPTIMIZATION 269 However, an alternative d.c. decomposition is provided in [5]: given ˆx∈I, f(x)=f1(x) +f2(x) (9) where f1(x) =f(ˆx) +f0(ˆx)(x −ˆx) +Zx ˆx (x −t)[f00(t)]+dt f2(x) =−Zx ˆx (x −t)[f00(t)]−dt being A+=max(A, 0)and A−=min(−A, 0). This yields the d.c. decomposition f(x)=(f (x) +h(x)) −h(x), where h(x) =                    h1(x) x ∈[z1,a 1] h2(x) x ∈[a1,b 1] h3(x) x ∈[b1,a 2] h4(x) x ∈[a2,b 2] h5(x) x ∈[b2,a 3] h6(x) x ∈[a3,b 3] h7(x) x ∈[b3,z 2] h1(x) =1+5π(x −a1)sin(5πa1)−cos(5πa1)+a1(2x−a1) h2(x) =1−cos(5πx)+x2 h3(x) =1+5π(x −b1)sin(5πb1)−cos(5πb1)+b1(2x−b1) h4(x) =h2(x) +h3(x) −1−5π(x −a2)sin(5πa2)+cos(5πa2)−a2(2x−a2) h5(x) =1+h4(x) −h2(x) +5π(x −b2)sin(5πb2)−cos(5πb2)+b2(2x−b2) h6(x) =h2(x) +h5(x) −1−5π(x −a3)sin(5πa3)+cos(5πa3)−a3(2x−a3) h7(x) =1+h6(x) −h2(x) +5π(x −b3)sin(5πb3)−cos(5πb3)+b3(2x−b3) z1=−0.2ai=−1 5πarccos(−0.08 π2)+2 5(i −1)i=1,2,3 z2=1bi=1 5πarccos(−0.08 π2)+2 5(i −1)i=1,2,3 This alternative envelope yields sharper bounds than (8), as depicted in Figure 1, which shows the envelopes derived from (8) and (9) after four iterations of the algorithms using x1=0.5as starting point. For the sake of completeness, we have compared in Figure 2 the d.c. envelope from (9) and the Piyavskii envelope (using the quite sharp Lipschitz constant L= 17.51) again after four iterations using x1=0.5as starting point. 270 R.BLANQUERO AND E.CARRIZOSA Figure 1. Comparison of envelopes I: Example function (double-wide line), Breiman–Cutler envelope (dotted line), d.c. envelope (single-wide line) Figure 2. Comparison of envelopes II: Example function (double-wide line), Piyavskii envelope (dotted line), d.c. envelope (single-wide line) 2.2. SOLVING (Qn) 2.2.1. Convex maximization For the upper envelope (6), (Qn) can be formulated as a convex maximization problem with linear constrains, e.g. [8]. Indeed, (Qn)is equivalent to max (x,t)∈B(n) {f1(x) +t}, where B(n) =(x, t) :x∈X, t≤f2(xi)+ξ2(xi)0(x −xi)∀i,1≤i≤n ON COVERING METHOD FOR D.C. OPTIMIZATION 271 Hence, solving each (Qn)amounts to maximizing a convex function over the polyhedron B(n). Moreover, since we obtain B(n+1)from B(n) by adding a linear constraint, the resolution of (Qn)can be used to solve (Qn+1)if procedures for on-line enumeration of vertices, [9], are used. Moreover, the insertion of each new constraint can be done in polynomial time (the dimension mconsidered to be fixed), since, at the n-th iteration, O(nm−1)time suffices to update the geometry of the feasible region, see [1, 11]. In other words, the covering Algorithm 1 for the d.c. case (and, in particular, the Breiman–Cutler method) is simply an outer approximation algorithm. 2.2.2. Power diagrams A different though equivalent strategy for optimizing (Qn)is proposed in [7] for the particular case (4): defining for each i, n,DBrC i:nas DBrC i:n=x∈Rm:kBrC i(x) ≤kBrC j(x) ∀j,1≤j≤n, one immediately obtains that each DBrC i:nis polyhedral. Moreover, since (Qn)can be rewritten as max 1≤i≤n(max x∈X∩DBrC i:n kBrC i(x)),(10) solving (Qn)amounts to inspecting the vertices of each polytope X∩DBrC i:n,for which a geometrical method is detailed in [7]. This strategy is applicable not only to covering methods derived from kBrC,but also for the more general case in which the envelopes are derived from kdc. Indeed, kdc i(x) ≤kdc j(x) ⇔f2(xi)+ξ2(xi)0(x −xi)≤f2(xj)+ξ2(xj)0(x −xj) ⇔kx−Cik2+βi≤ x−Cj  2+βj where Ck=−1 2ξ2(xk)β k=f2(xk)−ξ2(xk)0xk−1 4kξ2(xk)k2k=i, j thus, defining Ddc i:nas Ddc i:n=x∈Rm:kdc i(x) ≤kdc j(x) ∀j,1≤j≤n, it turns out that {Ddc i:n}1≤i≤nis a power diagram, with centers C1,... ,C nand (additive) weights β1,... ,β n, [1, 2, 11, 18]. The key property of power diagrams is the fact that they induce a polyhedral subdivision of Rm(as the particular case {DBrC i:n}), and can be constructed by online procedures, [2]. In fact, the geometrical procedure described in [7] for solving (Qn)is basically the same than the on-line procedure of [2] for the (at first glance unrelated) problem of describing the geometrical structure of a power diagram. 272 R.BLANQUERO AND E.CARRIZOSA Table 1. Number of iterations for unidimensional test problems No. Piyavskii Breiman–Cutler DC covering 1 149 25 12 2 155 21 14 3 195 103 54 4 413 24 14 5 151 32 15 6 129 37 11 7 153 24 15 8 185 85 53 9 119 24 14 10 203 24 14 11 373 44 28 12 327 42 29 13 993 264 15 14 145 29 13 15 629 80 16 16 497 88 15 17 549 67 26 18 303 19 14 19 131 20 14 20 493 30 11 A geometric-based strategy is also possible for envelopes based on kPIY for Lipschitz optimization. Indeed, one can define DPIY i:n, and obtain an expression similar to (10). This is the basis of the Jaumard–Herrmann–Ribault algorithm described in [14]. The so defined sets DPIY i:nmay still be identified with geometrical objects, namely, additive Voronoi diagrams, [1, 11, 18], but (for dimension greater or equal to 2) are no longer polyhedral, thus more complicated data structures are needed in order to store and update the geometry of such diagram. 3. Computational experience In this section we present some numerical results, which show that the d.c. decomposition chosen in (6) may have a great influence on the convergence of the algorithm. See also [6] for further results. We have studied the performance of Algorithm 1 using a set of well-known test problems from the literature, collected in [7, 13]. Table 1 shows the number of iterations spent by the three covering methods addressed throughout the paper ON COVERING METHOD FOR D.C. OPTIMIZATION 273 Table 2. Computational results for multidimensional test problems Example Breiman–Cutler DC covering Iterations CPU CPU Iterations CPU CPU standard time seconds standard time seconds COS2 77 2.6404 0.0059 27 1.7696 0.0040 COS4 1392 1416.7967 3.1776 201 250.8309 0.5748 H3 2575 443.2186 1.0094 442 117.6784 0.2648 for twenty one-dimensional test problems (see Table 1 of [13]). The computational results concerning to the Piyavskii and Breiman–Cutler methods are taken from [20], where the exact value of the associated parameters Land Khave been used, as well as an accuracy of value ε=0.0001 (b −a), with [a,b]being the interval where the test function is defined. The computational implementation of Algorithm 1 has been carried out in every case using the d.c. decomposition given in Section 6 of [15] for the non-C2function 18 and the d.c. decomposition (9) for the remaining instances. The starting point for the algorithm has been always taken equal to the middle point of the interval. We just provide the number of iterations since the CPU times for our algorithm are negligible. Finally, we have considered three multidimensional test problems from [7, 10] (COS2, COS4 and H3 in Table 3 of [7]), where the d.c. representations used are based on decomposition (9) for each term in COS2 and COS4, and the decomposition suggested by Proposition 3.5 of [22] for H3. The stopping rule and the tolerance values ε1and ε2has been taken as in [7]. In Table 2 we show the number of iterations and the CPU time (measured in units of standard time and seconds) taken by our algorithm in order to solve these problems. We just need to say that the implementation of Algorithm 1 used in these examples is directly based on the program listings of the Breiman–Cutler algorithm provided in [10]. Both of them have been programmed in Fortran 90 in a personal computer with a 200 Mhz. Pentium processor and 16 Megabytes of RAM memory. 4. Acknowledgements The authors want to thank Dr. J. Fliege for his helpful remarks. References 1. Aurenhammer, F. (1987), Power diagrams: Properties, algorithms and applications, SIAM J. on Computing 16: 78–96. 2. Aurenhammer, F. (1991), Voronoi Diagrams – A Survey of a Fundamental Geometric Data Structure, ACM Computing Surveys 23, 345–495.