Global optimization using space filling curves

MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Global Optimization Using Space Filling Curves Michaela BAILOVA, Jiri BOUCHALA, Petr VODSTRCIL Department of Applied Mathematics, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, Czech Republic michaela.bailov[email protected], jiri.bouc[email protected], petr.vo[email protected] DOI: 10.15598/aeee.v15i2.2303 Abstract. The existence of space filling curves opens the way to reducing multivariate optimization problems to the minimization of univariate functions. In this paper, we analyze the Hoelder continuity of space filling curves and exploit this property in the solution of global optimization problems. Subsequently, an algorithm for minimizing univariate Hoelder continuous functions is presented and analyzed. It is shown that the algorithm computes the approximate minimum with the guaranteed precision. The algorithm is tested on some types of two-dimensional functions. Keywords Dimension reduction, global optimization, Hoelder continuity, Lipschitz continuity, space filling curves. 1. Introduction Many engineering problems lead to a multivariate global optimization. Such issue can be difficult to solve. An approach presented in this paper is to turn the problem into its one-dimensional equivalent. A way how to provide such simplification is to use a space filling curve. This paper deals with the problem of searching a global minimum F(y∗) = min y∈DF(y),(1) and a global minimizer y∗∈D, where Dis an Ndimensional hypercube defined as follows D={y∈RN,−1 2≤yj≤1 2,1≤j≤N}.(2) The objective function Fis assumed to satisfy the Lipschitz condition with a constant L,0< L < ∞. The main idea of the paper is to turn the multivariate optimization problem into its one-dimensional equivalent, which can be solved using some techniques evolved for univariate optimization. One way how to do so is to develop a continuous correspondence ymapping a one-dimensional interval onto the hypercube. The problem Eq. (1) turns into the following one F(y∗) = F(y(x∗)) = min x∈[0,1] F(y(x)),(3) where x∗∈[0,1]. A complete analysis is done to prove that a univariate algorithm based on Hoelder continuity can be used to find an approximation of the optimal value of Fover the domain D. The performance of the algorithm is illustrated by numerical experiments. 2. Space Filling Curves Main ideas in this section are inspired by [1] and [3]. Definition 1. A space filling curve is a single-valued continuous correspondence ymapping the unit interval [0,1] onto the hypercube Dfrom Eq. (2). If yis a space filling curve, then F(y∗) = min x∈[0,1] F(y(x)).(4) Though the concept of a space filling curve is useful for the analysis of the algorithm, the effective computation is usually based on a continuous correspondence ynmapping the unit interval only into D. The following theorem gives some information about the optimal value of F◦ynusing the quality of yn. Theorem 1. Let (yn), where yn: [0,1] →D, be a sequence of curves such that sup y∈D dist(yn([0,1]), y) =: εn→0,(5) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 251 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE for n→ ∞ and let (x∗ n)be an arbitrary sequence in [0,1] satisfying F(yn(x∗ n) |{z } =:y∗ n ) = min x∈[0,1] F(yn(x)).(6) Then: 1.(∀n∈N) : 0 ≤F(y∗ n)−min y∈DF(y)≤Lεn,(7) 2. if ey∈Dis an accumulation point of (y∗ n), then F(ey) = miny∈DF(y). Proof. 1. Let us assume that F(y∗) := min y∈DF(y)≤F(y∗ n)≤F(yn(xn)),(8) where xn∈[0,1] is chosen so that kyn(xn)−y∗k ≤ εn. Hence 0≤F(y∗ n)−F(y∗)≤F(yn(xn)) −F(y∗)≤ ≤L(kyn(xn)−y∗k)≤Lεn.(9) The last inequality is based on the Lipschitz continuity of Fon Dand on the quality of yn. This completes the proof of Eq. (7). 2. If eyis an accumulation point of (y∗ n), then there exists a subsequence (for the sake of clarity labeled the same way as the original sequence) such that y∗ n→ey. (10) Using the continuity of Fleads to F(y∗ n)→F(ey),(11) but, from Eq. (9), we also get F(y∗ n)→F(y∗).(12) It follows that F(ey) = F(y∗) = min y∈DF(y).(13) Remark 1. If y∗= arg miny∈DF(y)is unique, then y∗ n→y∗.(14) Using ynreduces the multidimensional optimization problem into its one-dimensional equivalent that can be solved by univariate algorithms. If such method provides a lower bound Mnof the onedimensional function F◦yn, then Mnis obviously a lower bound of Falong yn. Is it possible to establish a lower bound for Fover the whole D? The following theorem gives an answer. Theorem 2. Assume that the curve yn,n∈N, satisfies the assumptions of Thm. 1 and (∀x∈[0,1]) : Mn≤F(yn(x)).(15) Then the value M=Mn−Lεn(16) is a lower bound of Fover the entire D, i.e. M≤min y∈DF(y).(17) Proof. Using the result and notation of the previous theorem, we get (∀y∈D) : F(y)≥F(y∗) = (F(y∗)− −F(yn(xn))) + F(yn(xn)) ≥ −Lεn+Mn.(18) In what follows, we assume that F◦ynis Hoelder continuous with some real constants H≥0, α ∈(0,1i, i.e. (∀x0, x00 ∈[0,1]) : |F(yn(x0)) −F(yn(x00))| ≤ H|x0−x00|α.(19) Let us describe such curve for N-dimensional case. At first, divide the domain Dinto 2Nequal hypercubes. Number all of the subcubes using the index z1, 0≤z1≤2N−1. Each subcube with the index z1is designated D(z1). Moreover, D= 2N−1 [ i=0 D(zi).(20) Using the same approach, divide each of the subcubes from the previous partitioning into 2Nequal subcubes and number them with the index z2,0≤z2≤ 2N−1. Each subcube of the second partitioning is now designated D(z1, z2). Continuing the same process we get hypercubes D(z1, z2, . . . , zM)and the edge length will be 2−M. The total number of the subcubes in M-th partition will be equal to 2MN and D⊃D(z1)⊃D(z1, z2)⊃. . . ··· ⊃ D(z1, z2, . . . , zM).(21) Now cut the interval [0,1] into 2Nequal subintervals. Every single part is designated d(z1),0≤z1≤2N− 1. In the same manner, cut all the subintervals once again, etc. Continuing the same process we get 2MN subintervals with the length equal to 2−MN , which are designated d(z1, z2, . . . , zM). Moreover, [0,1] ⊃d(z1)⊃d(z1, z2)⊃. . . ··· ⊃ d(z1, z2, . . . , zM).(22) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 252 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Fig. 1: The first and the second partition of the unit interval. The process is illustrated by Fig. 1. The interval d(z1, z2, . . . , zM)can be also written as d(z1, z2, . . . , zM)=[x, x + 2−MN ],(23) and can be referred to as d(M, x). The corresponding subcube D(z1, z2, . . . , zM)is designated D(M, x). The process of partitioning has to satisfy the following condition. Condition 1. If the subintervals d(M, v0)and d(M, v00),M∈N, have a common end point, then the corresponding subcubes D(M, v0)and D(M, v00) have a common face. Now consider a space filling curve y: [0,1] onto →D such that (∀M∈N0) : y(d(z1, . . . , zM)) ⊂D(z1, . . . , zM).(24) Then F◦yis Hoelder continuous with the constants 2L√N+ 3 and 1 N, i. e. for all x0, x00 ∈[0,1] |F(y(x0)) −F(y(x00))| ≤ 2L√N+ 3(|x0−x00|)1 N.(25) The proof can be found in [1]. In what follows, we use Hilbert-type curves. The process of partitioning for two-dimensional Hilbert-type curve is illustrated by the following figure. Fig. 2: The first and the second partition of the cube Dfor twodimensional Hilbert-type curve. For the practical computation, we use only "iterates" of the Hilbert curve. Consider the result of the M-th partition illustrated in this section. The construction of the "M-th iteration" of the Hilbert curve yMis described for example in [1] and [3]. The curve yM: [0,1] into →D, (26) satisfies the following properties. Let 0 = x0< x1< ··· < x2M N = 1 be the end points of the subintervals d(z1, . . . , zM)and let yi:= yM(xi), i ∈ {0,...,2MN }.(27) Then yM(x) = yi+ (yi+1 −yi)2MN (x−xi),(28) kyi+1 −yik ≤ 2−M,(29) where x∈[xi, xi+1], i ∈ {0,...,2MN −1}, and (∀K < M) : xi∈d(z1, . . . , zK)⇒ yi∈D(z1, . . . , zK).(30) The curves yMfor M= 0,1,2and N= 2 are illustrated by Fig. 3. The construction of the curves Fig. 3: Iterations of the Hilbert curve for N= 2. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 253 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE is based on complex transformations and is completely described in [3]. The next goal is to prove that also F◦ yMsatisfies the Hoelder condition similar to Eq. (25). Theorem 3. The function F◦yM,M∈N, fulfills the Hoelder condition with the constants 2L√N+ 3 and 1 Non interval [0,1], i.e. for all x0, x00 ∈[0,1] |F(yM(x0)) −F(yM(x00))| ≤ ≤2L√N+ 3(|x0−x00|)1 N.(31) Proof. Let x0,x00 ∈[0,1],x06=x00. Then there exists n∈N0such that 2−(n+1)N≤ |x0−x00| ≤ 2−nN .(32) The proof can be divided into two parts. 1. Suppose that n≥M. If x0and x00 are from the same subinterval d(M, xi), then, according to the properties of yM, the points yM(x0)and yM(x00)belong to the same linear segment with end points yi,yi+1. From Eq. (28), Eq. (29) and Eq. (32), we get kyM(x0)−yM(x00)k≤kyi+ (yi+1 −yi)2MN · ·(x0−xi)−yi−(yi+1 −yi)2MN (x00 −xi)k= = 2MN kyi+1 −yik|x0−x00| ≤ 2MN 2−M2−nN . (33) Using the first inequality in Eq. (32), it follows that kyM(x0)−yM(x00)k ≤ ≤2(M−n)(N−1)2|x0−x00|1/N ≤2|x0−x00|1/N .(34) Hence |F(yM(x0)) −F(yM(x00))| ≤ LkyM(x0)−yM(x00)k ≤ 2L|x0−x00|1 N.(35) If x0∈d(M, xi)and x00 ∈d(M, xi+1), then the points yM(x0),yM(x00)belong to two different segments with a common end point yi+1. Using the previous result, we get kyM(x0)−yM(x00)k ≤ kyM(x0)−yi+1k+ +kyM(x00)−yi+1k ≤ 4|x0−x00|1/N .(36) Hence |F(yM(x0)) −F(yM(x00))| ≤ 4L|x0−x00|1 N.(37) 2. Now suppose that n∈N0,n<M. If x0,x00 are both from d(n, ˜xi), then, from Eq. (30), yM(x0),yM(x00)∈ D(n, ˜xi). The maximal distance can be estimated as follows kyM(x0)−yM(x00)k ≤ 2−n√N. (38) If x0∈d(n, ˜xi),x00 ∈d(n, ˜xi+1), then yM(x0)∈D(n, ˜xi)and yM(x00)∈D(n, ˜xi+1). According to the Cond. 1, D(n, ˜xi)and D(n, ˜xi+1)are contiguous. Therefore kyM(x0)−yM(x00)k ≤ 2−n√N+ 3.(39) Using Eq. (32), we can derive the final estimate kyM(x0)−yM(x00)k ≤ 2√N+ 3|x0−x00|1/N .(40) Since the function Fis Lipschitz continuous, it follows that |F(yM(x0)) −F(yM(x00))| ≤ ≤LkyM(x0)−yM(x00)k ≤ ≤2L√N+ 3|x0−x00|1/N , (41) which completes the proof. Remark 2. Let MMbe the lower bound of F◦yMand let N= 2. Then (∀y∈D) : MM−L2−M≤F(y).(42) 3. Optimization Algorithm The algorithm itroduced in this section is inspired by [1], [2] and is designed for the optimization of the univariate functions that are Hoelder continuous. Consider g: [a, b]→Rsatisfying the Hoelder condition with the constants H,α,a, b ∈R,a < b. Let x0, x1, . . . , xnbe the trial points obtained in the previous iterations. These points divide [a, b]into nintervals. Let Iibe an interval with end points xj,xk that are consecutive, j, k ∈ {0, . . . , n},xj< xk. Then the lower bounding function on the interval Iiis constructed as follows: ln i(x) := max{lL i(x), lR i(x)},(43) where lL i(x) = g(xj)−H(x−xj)α,(44) lR i(x) = g(xk)−H(xk−x)α.(45) The functions lL iand lR iare illustrated by Fig. 4. It can be shown that the function Ln(x) := ln i(x), x ∈Ii, i ∈ {0, . . . , n −1}(46) is a lower bounding of gover [a, b]. Fig. 4: The functions lL iand lR i. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 254 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE The algorithm computes a value yi∈Iias an xcoordinate of the intersection of lL iand lR i. Then, the characteristic Miis computed in the following way Mi=lL i(yi) = lR i(yi).(47) It is true that Mi≤inf{g(x), x ∈Ii}. In Fig. 4, we can see a graphical interpretation of yiand Mi. The interval Itwith the lowest value of Mtis chosen and xn+1 := ytbecomes the next trial point, that divides Itinto two subintervals. For both of them, new characteristics are computed. The algorithm goes on until diam(It)≥δ. (48) Let us denote the approximation of the minimizer generated by the algorithm by ¯x. Let εbe the desired precision of the algorithm, i.e. the goal is to find ¯x∈[a, b]such that g(¯x)−min x∈[a,b]g(x)≤ε. (49) The value δ=δεin Eq. (48) can be chosen so that δ≤2ε H1 α.(50) The steps of the algorithm can be described as follows: •first iteration: Set x0=aand x1=b and compute the values g(x0)and g(x1). Then the functions lL 0and lR 0are constructed, similarly to the Fig. 4. The point y0is found as the xcoordinate of their intersection and the characteristic M0is computed using Eq. (47). If M0=g(y0), the optimization is done and ¯x=y0. Otherwise, the algorithm sets the next trial point x2=y0, that divides I0into two subintervals, I0= [x0, x2] and I1= [x2, x1]. For both intervals, the values y0,y1and the corresponding characteristics M0, M1are computed. •n-th iteration Let x0, x1, . . . , xnbe the trial points gained from the previous iterations, not necessarily sorted. Let I0, I1, . . . , In−1be the intervals (generated in the previous steps) with the end points xi,i∈ {0, . . . , n}. These intervals are characterized by the values y0, y1, . . . , yn−1and M0, M1, . . . , Mn−1computed in a way described above. The algorithm chooses the interval Itsuch that Mt= min{M0, M1, . . . , Mn−1},(51) and sets xn+1 =yt. If Mt=g(yt), then ¯x=ytand the optimization is done. Otherwise, the point ytdivides Itinto two subintervals, Itand In. For both subintervals, the points yt,ynand the characteristics Mt,Mnare computed. If diam(It)≤δ, (52) the algorithm computes the approximation of the optimal value g∗= min {g(xi),0≤i≤n+ 1},(53) and the minimizer ¯x= arg min {g(xi),0≤i≤n+ 1}.(54) Otherwise, the next iteration is done in the same manner. For illustration, the first 6 iterations for the function g(x)=(x−0.3)2+ 1,(55) on the interval [0,1] were chosen. The red lines illustrate the curves lL iand lR iand the optimal values of Mi are marked by the red dots. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.4 0.6 0.8 1 1.2 1.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 Fig. 5: First 6 iterations for the function from Eq. (55). The algorithm can be used to compute the optimal value of Falong the curve yM. As we have shown in the previous section, the function F◦yMis Hoelder continuous on [0,1] with the constants α=1 Nand H= 2L√N+ 3. In the following observations, we assume that N= 2. We get sup y∈D dist(yM([0,1]), y) = 1 2M.(56) c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 255 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE Let us return to the problem Eq. (1). If ε > 0is a desired precision, i.e. if we want to find ¯y∈Dsuch that F(¯y)−min y∈DF(y)≤ε, (57) then we can choose the parameters Mand δso that M≥log2 2L ε,(58) δ≤2ε 2H2 .(59) To prove the statement, we use Thm. 1, the quality of yM, the inequality Eq. (50) and a triangular inequality and we get |F(¯y)−min y∈DF(y)|≤|F(¯y)−min x∈[0,1] F(yM(x))|+ +|min x∈[0,1] F(yM(x)) −min y∈DF(y)| ≤ ε 2+ε 2=ε. (60) The algorithm described above is rather theoretical. For the practical computation, some simplifications have to be made. For the two-dimensional case it is not difficult to compute the intersection of lL iand lR i, but for the high-dimensional problem it is more useful to approximate lL iand lR iby two lines, as we can see in Fig. 6. The value ˜yican be computed as the first Fig. 6: The approximations of lL iand lR i. coordinate of their intersection and ˜ Mi= min{lL i( ˜yi), lR i( ˜yi)}.(61) It is also more convenient to use an approximated value of the Hoelder constant. In [2], two possible approximations are proposed. The choice of the interval Itcan be improved. Two methods are presented in [2]. 4. Numerical Experiments For our experiments, three types of functions were chosen. Let us start with the function F1(x, y) = (x−0.3)2+ (y−0.7)2+ 1.(62) It is not difficult to realize that y∗= [0.3,0.7] and F(y∗) = 1. The algorithm was tested for M= 1,2,...,10 and the precision 10−3. The results are summarized in the Tab. 1. Tab. 1: Results summarization. MF(yM(x)) −F(y∗)yM(x)kyM(x)−y∗k 14·10−2[0.2998,0.5000] 0.2000 22.5·10−3[0.3009,0.7500] 0.0500 35.1·10−3[0.2500,0.7518] 0.0720 41.5·10−4[0.3125,0.7016] 0.0126 51.5·10−4[0.3016,0.6875] 0.0126 67.6·10−5[0.3082,0.7031] 0.0087 71.1·10−4[0.2891,0.6994] 0.0110 85.6·10−5[0.2930,0.7026] 0.0075 92.9·10−5[0.2949,0.7019] 0.0054 10 1.8·10−5[0.2959,0.7013] 0.0043 In Fig. 7, we can see the solution for M= 5. 0 00.5 0.5 1 1 0.8 1.5 0.6 2 0.4 0.2 1 0 Fig. 7: The function F1together with the fifth level of the Hilbert curve and the approximate minimizer. Consider next the multiextremal function F2(x, y) = −0.5 sin(2πx) sin(2πy)+1.(63) The algorithm found one of the minimizers and the result is illustrated by Fig. 8. The last tested function is F3(x, y) = p(x−0.6)2+ (y−0.4)2+ 0.5.(64) As we can be seen Fig. 9, the algorithm works also for the nonsmooth function. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 256 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 15 |NUMBER: 2 |2017 |JUNE 1 0 10.5 0.5 0.8 1 0.6 0.4 1.5 0.2 0 0 Fig. 8: The function F2together with the fifth level of the Hilbert curve and the approximate minimizer. 0 0.2 1 0.4 1 0.6 0.8 1 0.8 1.2 1.4 0.6 0.5 0.4 0.2 00 Fig. 9: The function F3together with the fifth level of the Hilbert curve and the approximate minimizer. 5. Conclusion In this paper, we described a minimization algorithm based on reducing the problem by means of Hilberttype space filling curve. Subsequently, complete mathematical analysis was done to show that the function F◦yMis Hoelder continuous. Thus an algorithm for minimizing univariate Hoelder continuous functions could be used to find an approximation of the optimal value of F. It was shown that the algorithm computed the approximate minimum of Fwith the guaranteed precision. The algorithm and the analysis will be extended for the dimensions N > 2and used for some practical problems of the multivariate optimizations. Acknowledgment The work is partially supported by Grant of SGS No. SP2017/122, VSB–Technical University of Ostrava, Czech Republic. References [1] SERGEYEV, Y. D, R. G. STRONGIN and D. LERA. Introduction to Global Optimization Exploiting Space-Filling Curves. New York: Springer, 2013. ISBN 978-1-4614-8041-9. [2] LERA, D. and Y. D. SERGEYEV. Lipschitz and Hoelder global optimization using space-filling curves. Applied Numerical Mathematics. 2010, vol. 60, iss. 1–2, pp. 115–129. ISSN 0168-9274. DOI: 10.1016/j.apnum.2009.10.004. [3] SAGAN, H. Space-Filling Curves. New York: Springer, 1994. ISBN 978-0-387-94265-0. [4] KUFNER, A., O. JOHN and S. FUCIK. Function spaces. Leyden: Noordhoff International Publishing, 1977. ISBN 90-286-0015-9. [5] BUTZ, A. R. Space filling curves and mathematical programming. Information and Control. 1968, vol. 12, iss. 4, pp. 314–330. ISSN 0019-9958. DOI: 10.1016/S0019-9958(68)90367-7. About Authors Michaela BAILOVA was born in Ostrava, Czech Republic. She received her M.Sc. from VSB–Technical University of Ostrava in 2015. Her research interests include functional and numerical analysis. Jiri BOUCHALA was born in Novy Jicin, Czech Republic. He received his Ph.D. from Faculty of Applied Sciences, University of West Bohemia in Pilsen in 2000. His research interests include functional analysis, boundary element methods, variational methods, partial differential equations. Petr VODSTRCIL was born in Svitavy, Czech Republic. He received his Ph.D. from Masaryk University in Brno in 2005. His research interests include calculus and functional analysis. c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 257