Global Optimization Using Space Filling Curves

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.


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 and a global minimizer y * ∈ D, where D is an Ndimensional hypercube defined as follows The objective function F is 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 y mapping a one-dimensional interval onto the hypercube.The problem Eq. ( 1) turns into the following one 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 F over the domain D. The performance of the algorithm is illustrated by numerical experiments.

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 y mapping the unit interval [0, 1] onto the hypercube D from Eq. ( 2).
If y is a space filling curve, then 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 y n mapping the unit interval only into D. The following theorem gives some information about the optimal value of F • y n using the quality of y n .
Theorem 1.Let (y n ), where y n : [0, 1] → D, be a sequence of curves such that for n → ∞ and let (x * n ) be an arbitrary sequence in [0, 1] satisfying Then: Proof.

Let us assume that
where The last inequality is based on the Lipschitz continuity of F on D and on the quality of y n .This completes the proof of Eq. (7).
2. If y is 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 → y. (10) Using the continuity of F leads to but, from Eq. (9), we also get It follows that Using y n reduces the multidimensional optimization problem into its one-dimensional equivalent that can be solved by univariate algorithms.
If such method provides a lower bound M n of the onedimensional function F • y n , then M n is obviously a lower bound of F along y n .Is it possible to establish a lower bound for F over the whole D? The following theorem gives an answer.
Theorem 2. Assume that the curve y n , n ∈ N, satisfies the assumptions of Thm. 1 and Then the value is a lower bound of F over the entire D, i.e.
Proof.Using the result and notation of the previous theorem, we get In what follows, we assume that F • y n is Hoelder continuous with some real constants H ≥ 0, α ∈ (0, 1 , i.e. Let us describe such curve for N -dimensional case.At first, divide the domain D into 2 N equal hypercubes.Number all of the subcubes using the index z 1 , 0 ≤ z 1 ≤ 2 N − 1.Each subcube with the index z 1 is designated D(z 1 ).Moreover, Using the same approach, divide each of the subcubes from the previous partitioning into 2 N equal subcubes and number them with the index z 2 , 0 ≤ z 2 ≤ 2 N − 1.Each subcube of the second partitioning is now designated D(z 1 , z 2 ).Continuing the same process we get hypercubes D(z 1 , z 2 , . . ., z M ) and the edge length will be 2 −M .The total number of the subcubes in M -th partition will be equal to 2 M N and In the same manner, cut all the subintervals once again, etc. Continuing the same process we get 2 M N subintervals with the length equal to 2 −M N , which are designated d(z 1 , z 2 , . . ., z M ).Moreover, The process is illustrated by Fig. 1.The interval d(z 1 , z 2 , . . ., z M ) can be also written as and can be referred to as d(M, x).The corresponding subcube The process of partitioning has to satisfy the following condition.
Then F • y is Hoelder continuous with the constants 2L √ N + 3 and 1 N , i. e. for all x , x ∈ [0, 1] 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.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 y M is described for example in [1] and [3].The curve satisfies the following properties.Let be the end points of the subintervals d(z 1 , . . ., z M ) and let Then The curves y M for M = 0, 1, 2 and N = 2 are illustrated by Fig. 3.The construction of the curves is based on complex transformations and is completely described in [3].The next goal is to prove that also F • y M satisfies the Hoelder condition similar to Eq. (25). (31) The proof can be divided into two parts.
1. Suppose that n ≥ M .If x and x are from the same subinterval d(M, x i ), then, according to the properties of y M , the points y M (x ) and y M (x ) belong to the same linear segment with end points y i , y i+1 .From Eq. (28), Eq. ( 29) and Eq.(32), we get Using the first inequality in Eq. (32), it follows that Hence ), then the points y M (x ), y M (x ) belong to two different segments with a common end point y i+1 .Using the previous result, we get Hence Using Eq. (32), we can derive the final estimate which completes the proof.
Remark 2. Let M M be the lower bound of F • y M and let N = 2. Then

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] → R satisfying the Hoelder condition with the constants H, α, a, b ∈ R, a < b.Let x 0 , x 1 , . . ., x n be the trial points obtained in the previous iterations.These points divide [a, b] into n intervals.Let I i be an interval with end points x j , x k that are consecutive, j, k ∈ {0, . . ., n}, x j < x k .Then the lower bounding function on the interval I i is constructed as follows: where The functions l L i and l R i are illustrated by Fig. 4. It can be shown that the function  The algorithm computes a value y i ∈ I i as an x coordinate of the intersection of l L i and l R i .Then, the characteristic M i is computed in the following way It is true that M i ≤ inf{g(x), x ∈ I i }.In Fig. 4, we can see a graphical interpretation of y i and M i .
The interval I t with the lowest value of M t is chosen and x n+1 := y t becomes the next trial point, that divides I t into two subintervals.For both of them, new characteristics are computed.The algorithm goes on until 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 The value δ = δ ε in Eq. ( 48) can be chosen so that The steps of the algorithm can be described as follows: • first iteration: Set x 0 = a and x 1 = b and compute the values g(x 0 ) and g(x 1 ).Then the functions l L 0 and l R 0 are constructed, similarly to the Fig. 4. The point y 0 is found as the x coordinate of their intersection and the characteristic M 0 is computed using Eq. ( 47).If M 0 = g(y 0 ), the optimization is done and x = y 0 .Otherwise, the algorithm sets the next trial point x 2 = y 0 , that divides I 0 into two subintervals, I 0 = [x 0 , x 2 ] and I 1 = [x 2 , x 1 ].For both intervals, the values y 0 , y 1 and the corresponding characteristics M 0 , M 1 are computed.
• n-th iteration Let x 0 , x 1 , . . ., x n be the trial points gained from the previous iterations, not necessarily sorted.Let I 0 , I 1 , . . ., I n−1 be the intervals (generated in the previous steps) with the end points x i , i ∈ {0, . . ., n}.These intervals are characterized by the values y 0 , y 1 , . . ., y n−1 and M 0 , M 1 , . . ., M n−1 computed in a way described above.The algorithm chooses the interval I t such that and sets x n+1 = y t .If M t = g(y t ), then x = y t and the optimization is done.Otherwise, the point y t divides I t into two subintervals, I t and I n .For both subintervals, the points y t , y n and the characteristics M t , M n are computed.If the algorithm computes the approximation of the optimal value and the minimizer Otherwise, the next iteration is done in the same manner.
For illustration, the first 6 iterations for the function on the interval [0, 1] were chosen.The red lines illustrate the curves l L i and l R i and the optimal values of M i are marked by the red dots.The algorithm can be used to compute the optimal value of F along the curve y M .As we have shown in the previous section, the function F • y M is Hoelder continuous on [0, 1] with the constants α = 1 N and H = 2L √ N + 3.In the following observations, we assume that N = 2.We get Let us return to the problem Eq. ( 1).If ε > 0 is a desired precision, i.e. if we want to find ȳ ∈ D such that then we can choose the parameters M and δ so that To prove the statement, we use Thm. 1, the quality of y M , the inequality Eq. ( 50) and a triangular inequality and we get (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 l L i and l R i , but for the high-dimensional problem it is more useful to approximate l L i and l R i by two lines, as we can see in Fig. 6.The value ỹi can be computed as the first coordinate of their intersection and 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 I t can be improved.Two methods are presented in [2].

Numerical Experiments
For our experiments, three types of functions were chosen.Let us start with the function 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.
M In Fig. 7, we can see the solution for M = 5.
As we can be seen Fig. 9, the algorithm works also for the nonsmooth function.

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 • y M is 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 F with the guaranteed precision.The algorithm and the analysis will be extended for the dimensions N > 2 and used for some practical problems of the multivariate optimizations.

cFig. 1 :
Fig. 1: The first and the second partition of the unit interval.

Condition 1 .
If the subintervals d(M, v ) and d(M, v ), M ∈ N, have a common end point, then the corresponding subcubes D(M, v ) and D(M, v ) have a common face.Now consider a space filling curve y : [0, 1] onto → D such that

Fig. 2 :
Fig. 2: The first and the second partition of the cube D for twodimensional Hilbert-type curve.

Fig. 4 :
Fig. 4: The functions l L i and l R i .

Fig. 6 :
Fig. 6: The approximations of l L i and l R i .

Fig. 7 :
Fig. 7: The function F 1 together with the fifth level of the Hilbert curve and the approximate minimizer.

Fig. 8 :
Fig. 8: The function F 2 together with the fifth level of the Hilbert curve and the approximate minimizer.

Fig. 9 :
Fig. 9: The function F 3 together with the fifth level of the Hilbert curve and the approximate minimizer.