Thomas Firmin
email: thomas.firmin@univ-lille.fr

El-Ghazali Talbi
email: el-ghazali.talbi@univ-lille.fr

A fractal-based decomposition framework for continuous optimization

Keywords: Continuous optimization, Metaheuristic, High-dimensional optimization, Decomposition, Fractal, Tree search 1 Bent

In this paper, we propose a generic algorithmic framework which defines a unified view of fractal decomposition algorithms for continuous optimization. Fractals allow building a hierarchical decomposition of the decision space by using a self-similar geometrical object. The proposed generic framework is made of five distinct and independent search components: fractal geometrical object, tree search, scoring, exploration and exploitation. The genericity of the framework allowed the instantiation of popular algorithms from the optimization, machine learning and computational intelligence communities. Moreover, new optimization algorithms can be designed using various strategies of the search components. This shows the modularity of the proposed algorithmic framework. The computational experiments emphasize the behaviors of fractal-based approaches in terms of scalability, robustness, and the balance between exploitation and exploration in the search space. The obtained results show the significance of each search component of the fractal framework, and the necessity to build harder and well-defined benchmarks which can assess the performance of deterministic, axis-aligned and symmetrical decomposition-based algorithms.

Introduction

Optimizing a non-linear, non-convex, derivative free, or a black-box function in a high dimensional continuous search space is a complex task. Commonly we consider a minimization problem for an objective function f : S ⊂ R n → R:

x ∈ argmin x∈S f (x) ( 1
)
where x is the global optima, f the objective function, and S a compact set made of inequalities (e.g. upper and lower bounds of the search space).

In the past two decades, several optimization algorithms have been proposed to tackle continuous optimization problems. Population-based algorithms generally rely on evolutionary algorithms (e.g. differential evolution, evolution strategies) and swarm intelligence (e.g. particle swarm optimization) [START_REF] Bansal | Particle Swarm Optimization[END_REF]. Local search strategies are based on a single solution improvement (e.g. gradientbased algorithms, simulated annealing). Most of the state-of-the-art optimization algorithms are stochastic and require a high number of evaluations to be efficient. For optimization problems with expensive objective functions, surrogatebased optimization (e.g. Bayesian optimization) can be an alternative. A major drawback of surrogate-based strategies is that they can suffer from the curse of dimensionality, and so, they poorly scale on high dimensional optimization problems [START_REF] Garnett | Bayesian Optimization[END_REF].

In this paper, we investigate a particular class of optimization algorithms, which could be classified as divide-and-conquer strategies based on hierarchical decomposition of the search space. This class of algorithms can overcome the dimensionality problem by creating flexible, scalable and massively parallel decomposition-based algorithms for high dimensional optimization problems.

This work has been inspired by two distinct families of decomposition-based optimization algorithms. The first one is based on algorithms derived from Lipschitzian global optimization [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF], such as DIRECT and its various extensions (e.g. eDIRECT, BIRECT, HD-DIRECT) [START_REF] Jones | The DIRECT algorithm: 25 years Later[END_REF]. The second one concerns metaheuristics based on fractal decomposition, such as FRACTOP [START_REF] Demirhan | FRACTOP: A Geometric Partitioning Metaheuristic for Global Optimization[END_REF] and FDA [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF].

The main contribution of this paper is to unify fractal-based decomposition approaches into a generic and flexible algorithmic framework. The modularity of the framework allows the instantiation of popular optimization algorithms from the global optimization, machine learning and computational intelligence communities. This flexibility also allows the design of new optimization algorithms by developing new strategies of the framework's elementary building blocks, here called search components.

In our framework, a fractal is a generalized concept describing a high dimensional geometrical object structuring the search space. It is also a subset of an initial decision space or of another fractal, and a node of a tree.

A fractal is a never-ending pattern that can be used to generate any-scale fractal trees [START_REF] Mandelbrot | The fractal geometry of nature[END_REF].

The self-reproducing characteristic of fractals suggests a new way for decomposing large-scale search spaces with low time and space complexity. Indeed, a fractal bears an unlimited number of levels that can be generated using simple and efficient analytical procedures with constant O(1) complexity. More exactly, a fractal has certain properties which can be exploited to better structure the search process including recursion, scalability and self-similarity. By considering at least these three fundamental properties, it is possible to explore a search space indefinitely (recursion) in a structured manner and on all scales (scalability and self-similarity). The search space is then decomposed in a cascade of fractals organized by levels. Starting from a first fractal supposed to cover the whole search space, each fractal at a given level gives rise recursively to smaller fractals.

A fractal-based decomposition algorithm is organized around five search components: fractal, tree search, scoring, exploration and exploitation. It decomposes the search space via a given fractal geometrical object. Tree search provides a dynamic and hierarchical fractal decomposition of the search space. A scoring search component allows to balance the exploration and exploitation of the search space by giving a fitness value to fractals. A generic software framework named Zellij1 has been developed. It offers a unified programming paradigm for the instantiation of state-of-the-art fractal-based decomposition algorithms such as DIRECT (DIviding RECTangles), SOO (Simultaneous Optimistic Optimization) and FDA (Fractal Decomposition Algorithm). Moreover, one can also design new optimization algorithms according to the definition of new strategies for the search components.

The paper is organized as follows. In section 2, the Zellij generic algorithmic framework based on fractal decomposition is presented. The following sections (section 3 to 6) detail successively the search components of the algorithmic optimization framework: fractal geometrical object, tree search, scoring, exploration, and exploitation. In section 7, state-of-the art decomposition-based algorithms have been instantiated to the framework and extended by defining new strategies for their search components. Section 8 depicts the experimental setup for comparing five popular fractal-based decomposition algorithms and six new designed ones. Section 9 analyzes the obtained computational results on CEC2020 and SOCO2011 benchmarks. Finally, we summed up in section 10, the framework, future works and new research opportunities for tackling high dimensional optimization problem.

Zellij: A fractal-based decomposition algorithmic framework

This section introduces the basic concepts and search components of the Zellij framework.

Definition 1 (Search space). Let us define a continuous search space S of dimension n as a bounded subset of a metric space:

S = L × U = n i=1 [l i , u i ] (2)
with L, U ∈ R n , the lower and upper bounds, such

that ∀i ∈ [1, ..., n], l i < u i and ∀x ∈ S, ∀i ∈ [1, ..., n], l i ≤ x i ≤ u i
In our framework, a fractal is a self-similar geometrical object which does not depend on any information besides the nature of its parent. As a simple example, decomposing a hypercube (i.e ncube) into smaller hypercubes of equal size, relies solely on the boundaries of the parent. To generalize our approach to other geometrical objects than fractals, we can also consider geometrical objects that depend on additional knowledge. For instance, in the DIRECT algorithm, the decomposition process in smaller hyperrectangles depends on the sampling and the evaluation of their centers. Selecting a particular geometrical object (e.g. hypercube, hypersphere) has a major impact on the decomposition process. Five properties characterize a fractal: partition size, building complexity, coverage, overlap, and memory complexity (see Table 1 and Figure 5). These properties are discussed in Section 3.

The decomposition of the search space can be defined by a k-ary rooted partition tree. The root represents the initial complete search space, and nodes correspond to the generated fractals. The design of a tree search algorithm is crucial to efficiently explore and exploit the fractals. By introducing pruning strategies, one can reduce the search space and tackle memory issues, by eliminating some fractals. Popular tree search algorithms are Best First Search [START_REF] Dechter | Generalized best-first search strategies and the optimality of a*[END_REF], Beam Search [START_REF] Frohner | Parallel Beam Search for Combinatorial Optimization (Extended Abstract)[END_REF] or Epsilon Greedy Search [START_REF] Valenzano | On the completeness of best-first search variants that use random exploration[END_REF].

In the design of an optimization algorithm, high performance requires a trade-off between the exploration of the search space and the exploitation of the acquired knowledge. One has to find the best strategy for this dilemma [START_REF] Morales-Castañeda | A better balance in metaheuristic algorithms: Does it exist?[END_REF]. Exploration (i.e. diversification) of the search space allows one to obtain new knowledge. Non-explored fractals must be visited to ensure that all fractals are evenly explored and that search is not confined to a reduced number of fractals. Exploitation (i.e. intensification) into a reduced region of the search space uses that knowledge (e.g. best found fractals) to improve it. The promising fractals are searched more thoroughly in the hope to find better solutions.

Sampling in high dimensional search spaces is a critical task. Moreover, one does not sample in the same way in a hypercube or in a hypersphere.

In DIRECT [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF] and SOO [START_REF] Munos | Optimistic Optimization of a Deterministic Function without the Knowledge of its Smoothness[END_REF], only the center of the hyperrectangle is used in sampling. In our framework, we consider both deterministic and stochastic sampling. We also consider unique or multiple points methods to sample inside a fractal. The sampling process is essential for both exploration and exploitation search components.

The exploration could be done in a passive way (e.g. Markov chain Monte Carlo (MCMC) sampling, low discrepancy sequences [START_REF] Drmota | Sequences, Discrepancies and Applications[END_REF]) or in an active way (e.g. metaheuristics [START_REF] Talbi | Metaheuristics[END_REF], surrogatebased optimization [START_REF] Garnett | Bayesian Optimization[END_REF]). The challenge here is to find efficient sampling algorithms for diverse and complex fractals such as polytopes.

When a fractal reaches the maximum depth of the partition tree, an exploitation algorithm can be applied to it. The exploitation phase has not to be constrained within a fractal. The only bounds will be ones from the initial search space so that the exploitation can move freely toward a local or global optimum. One can use local search strategies such as gradient-based algorithms or simulated annealing [START_REF] Talbi | Metaheuristics[END_REF].

The role of the scoring search component is to assign a quality, a fitness value, for a given fractal, by using acquired knowledge following an exploration phase. It can be seen as the quality value obtained by an acquisition function in Bayesian optimization. Many scoring methods can be used, such as minimum, mean, median, distance-to-thebest [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF][15], or belief [START_REF] Demirhan | FRACTOP: A Geometric Partitioning Metaheuristic for Global Optimization[END_REF].

The Zellij workflow is described in Figure 1.

One can identify the five search components, their interactions and their algorithmic behaviors. The two For each instructions correspond to line 14 and line 16 of Algorithm To sum up, three main properties characterize Zellij :

• Generalization: our goal is to build a framework which unifies and generalizes various popular fractal decomposition-based algorithms from different communities (e.g. global optimization, reinforcement learning, computational intelligence).

• Modularity: the framework must be as modular as possible, so one can easily develop new optimization algorithms using the fractal decomposition-based approach.

• Massively parallel: a transparent and efficient parallel implementation of the algorithms on various architectures (e.g. multicores, GPUs) is carried out. The main challenge is the parallelization of the tree search component of the framework. Many parallel tree search algorithms can be considered.

In the following sections, we will further describe the five search components, and introduce some theoretical background to fractal-based decomposition algorithms.

Geometrical fractal object

The fractal search component within Zellij framework allows structuring high dimensional search spaces to better explore and exploit it. Several types of fractals can be used in the decomposition of the search space, such as hyperspheres, hypercubes or Voronoï cells (see Figure 5). This geometrical object has a great impact on the behaviors of fractal-based decomposition algorithms.

Definition 2 (k-partition). A k-partition of a subset α ⊆ S can be defined by the union of k nonempty and disjoint subspaces. These subsets can be obtained by a fractal building function

F : α → {S i } i∈[1,...k] ,
where S i ⊂ α, and:

α = F (α) = k i=1 S i (3)
Here ∀i ∈ [1, .., k], S i ̸ = ∅, and k i=1 S i = ∅. This definition does not consider overlapping subspaces, nor partial partitions.

Using a given fractal, one can build a recursive decomposition of S or a decomposition of a subset of S. Building a k-partition of an existing subset is called a refinement. Thus, when refining a fractal, the definition of S does not change. The refinement is a self-similar object, it uses the same definition as S to build children of a fractal; all S i are of the same nature.

Definition 3 (Hierarchical k-refinement). Let D be the number of successive refinements of a subset S l,j,i ⊆ S. We write j ∈ [1, ..., E l-1 ], with E l the number of sets that have been refined l times, and i ∈ [1, ..., k] the ith set of a k-partition of a superset j. A hierarchical k-refinement of S l,j,i of maximum level D (i.e. depth) is written as:

Here, a subset identified by (l, j, i) corresponds to the fractal at level l, child number i of the fractal j at level l -1. An appropriate data structure used to model Definition 3 is a k-ary rooted tree. Hence, one can rewrite the initial search space S as the root of this tree: S (1,0,1) = S, where E 0 = ∅.

A fractal is now considered as a node of a k-ary rooted tree T . Figure 2 shows an example of a 2-refinement of depth 4, of a 2D square using a bisection, drawn as a red dotted line, along the longest side.

For a fixed l ∈ [1, ..., D-1], one can write the set P l of supersets (parents) at level l that have children at level l + 1. Let us denote a leaf at the level l of a tree T as •[l, j, i]. An interesting property of using a k-ary tree on a k-refinement that fully covers S is that the initial search space will be equal to the union of all the leaves.

Theorem 1. Let S be a search space and the root of a k-ary partition tree T of maximum depth D. For all 1 < l ≤ D, the union of all leaves •[l, j, i], children of nodes numbered by [1, ..., E l-1 ] as defined in Definition 3, is equal to S:

S = S (1,0,1) = D l=2 E l-1 j=1 k i=1 •[l, j, i] (5)
Proof. We consider a k-ary rooted tree T (d) of depth d, with 1 < d ≤ D. Considering P l the set of fractals at level l having children at level l + 1:

P d-1 = E d-1 j=1 k i=1 •[d, j, i]
If we remove all •[d, j, i] from T d , we obtain a tree 1) . So, the leaves of

T (d-
T (d-1) at level d -1 can be written as {•[d -1, j, i], P d-1 }. Thus, P d-2 = E d-2 j=1 k i=1 •[d -1, j, i] P d-1
and so on, until T (1) :

S = P 1 = k i=1 •[2, 1, i] P 2
For simplicity, one can write all leaves •[l, j, i] of T (d) , a tree T of depth d, as L T (d) .

As mentioned before, Theorem 1 is valid for an exhaustive and mutually exclusive partition. This is not the case for some fractal-based algorithms such as FDA, which uses hyperspheres [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF]. Indeed, fractals can overlap and do not necessarily cover S. Additionally, it is important to mention that, the higher the dimension of the problem, the smaller the covering of the hyperspheres. In FDA, the covering is maximized according to a parameter, called the inflation ratio, no matter the overlap between fractals. However, the inflation ratio can make a part of the hypersphere going outside S. In our algorithmic framework, we consider the intersection between fractals and S. A fractal, or a part of a fractal, cannot be outside the initial search space S. If so, then the fractal is trimmed.

We previously described the basic principles of fully covering fractal-based decomposition algorithms. To extend our framework to partial partitions, we have to measure what the overlap and coverage are for fractal-based decomposition algorithms. For a tree T (d) , let Σ = P(S) be a σ-algebra of S. So (S, Σ) is a measurable set, where

L T (d) ⊆ Σ and L T (d) ∈ Σ. We can define a measure, µ : Σ → [0, +∞]. A stricter condition is applied to µ, µ(A) = 0 ⇐⇒ A = ∅.
Hence, the result of a k-partition cannot be made of null sets. Because F cannot produce null sets, and because D, the maximum tree depth, is finite, we can define a measure of the coverage of S by L T (d) and their overlap.

Definition 4 (Coverage). Let (S, Σ, µ) a measure space, and A, B ⊆ S, A, B ∈ Σ and B ⊆ A. A measure of the coverage of A by B can be written as:

C(A, B) = µ ∁ A B (6)
However if B ⊈ A such as in FDA:

C(A, B) = µ A \ A B (7)
From Equation 7we can infer B \ A, the part of B outside A.

One can distinguish two types of coverage: C(S, L T (d) ) measures the coverage of the initial fractal by L T (d) , and C(S l,j,i , F (S l,j,i )) measures the coverage of a fractal by its children.

Definition 5 (Overlap). Let (S, Σ, µ) a measure space, and A ̸ = ∅ a countable set such that ∀n ∈ N, ∀a n ∈ A : a n ∈ Σ and n∈N a n ∈ Σ. A measure of the overlap between all a n can be written as:

O(A) = µ   |A|-1 i=1 |A| j=i+1 a i a j   (8)
Coverage and overlap properties are illustrated in Figure 3 as a 2-dimensional FDA scheme. One can see, in hatched-magenta, the uncovered space and in solid-green the overlap between fractals. Moreover, because of the inflation of hyperspheres, the gridded-gray space represents parts of the fractals outside S. As the monotonic property of measures is not strict, A ⊆ B =⇒ µ(A) ≤ µ(B), one cannot say that L T (d) covers S by only looking at C(S, L T (d) ). Thus, several assumptions are made:

• For continuous dimensions, we consider that

L T (d) covers S if C(S, L T (d) ) = 0, even if L T (d) does not include S boundaries.
• By construction, we consider that F , the decomposition function, produces fractals that are strictly smaller than the parent:

∀A i ∈ F (A), µ(A i ) < µ(A)
Hence, we can modify Theorem 1 for overlap and partial k-partition.

Theorem 2. Let S be a search space, Σ = P(S) a σ-algebra on S, and (S, Σ, µ) a measure space. S is the root of a k-ary partition tree T (d) of depth d and L T (d) its leaves. C(., .) is a measure of the coverage according to Definition 4. We suppose that L T (d) is not a null set and µ(A) = 0 ⇐⇒ (A = ∅) ∨ ( A = ∅). Then, a given decomposition-based algorithm is said to be:

1. Preservative: if C(S, L T (d) ) = 0. Then, µ(S) = µ L T (d) (9) 2. Low-sacrifice based: if C(S, L T (d) ) > 0 and C(S, F (S)) > 0 and if ∀l > 1, C(S l,j,i , F (S l,j,i )) = 0. We have, µ(S) > µ L T (d) (10) with ∀d > 1, µ ( L T (d) ) ≤ µ ( L T (d+1) )
3. High-sacrifice based: if C(S, L T (d) ) > 0 and C(S, F (S)) > 0 and if ∀l > 1, C(S l,j,i , F (S l,j,i )) > 0 and F (S l,j,i )/S l,j,i ̸ = ∅. Then, ∀d > 1, we have,

µ(S) > µ L T (d) > µ L T (d+1) (11)
Proof. We consider a k-ary rooted tree T of depth d, with 1 < d ≤ D. This tree is denoted as T (d) , and the leaves of T (d) as L T (d) . One can consider a decomposition function F defined by Equation 3.

Preservative:

If ∀d > 1, C(S, L T (d) ) = 0. Then S = L T (d) , because by construction L T (d) ⊆ S and C(A, B) = 0 ⇐⇒ µ (A \ (A B)) = 0 ⇐⇒ A = B.
So, we have:

S = L T (d) =⇒ µ(S) = µ L T (d)
Low-sacrifice based:

If ∀d > 1, C(S, L T (d) ) > 0 and C(S d,j,i , F (S d,j,i )) = 0. Then, S ⊃ L T (d) =⇒ µ(S) > µ ( L T (d) ) and L T (d) ⊆ L T (d+1) . So we have, ∀a ∈ L T (d) , F (a) ∈ L T (d+1) =⇒ µ(a) ≤ µ F (a) =⇒ µ L T (d) ≤ µ L T (d+1)
High-sacrifice based:

If ∀d > 1, C(S, L T (d) ) > 0, C(S d,j,i , F (S d,j,i )) > 0 and F (A)\A = ∅ (Subsets cannot have solutions outside their parents). Then, S ⊃ L T (d) =⇒ µ(S) > µ ( L T (d) ) and L T (d) ⊃ L T (d+1) .
So we have,

∀a ∈ L T (d) , F (a) ∈ L T (d+1) =⇒ µ(a) > µ F (a) =⇒ µ L T (d) > µ L T (d+1)
This theorem describes behaviors of a partial, k-refinement where leaves might not fully cover the initial search space. Preservative algorithms describe a partition that covers S no matter whether there is overlap between leaves or not, as we consider L T (d) .

However, low-sacrifice based algorithms describe leaves of the tree that do not fully cover S, but where children of the fractal cover its parent.

In this case, one accepts to lose a part of the search space at the first refinement of S. Then future refinements should not sacrifice solutions anymore. Moreover, sometimes the refinement of a fractal can result in children, which when unified form a bigger object than the parent. This is the case for fractals that can cover solutions outside the parent space.

According to Theorem 2, let us consider an initial search space S defined by a hypercube. In FDA [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF], the algorithm considers the inscribed hypersphere as the initial search space. The algorithm accepts to "sacrifice" points between the hypercube and the inscribed hypersphere to perform a k-refinement using hyperspheres. To overcome this, the exploitation phase of FDA can ignore fractals' borders, so it can reach uncovered solutions. The algorithm is based on a high sacrifice rate when such a sacrifice happens at each refinement. Hence, some solutions are lost at each refinement. So, the deeper is the tree, the fewer solutions are reachable by the algorithm.

For high-sacrifice based algorithms, one only considers fractals that cannot expand outside their parents. If so, it becomes hard to define the notion of sacrifice, as we can imagine a decomposition where successive children do not fully cover their parents and shift outside an ancestor fractal. One should be aware of this property, and determine if it is acceptable to lose some solutions to carefully choose the right decomposition function F .

These three levels of sacrifice are described in Figure 4. Preservative fractals in Figure 4a, based on SOO in this example, preserve the initial search space S, no space is lost. While low sacrifice fractals, on Figure 4b, here based on Sierpinski triangles, sacrifice space of S to use a geometry based on triangles. And finally, the example on Figure 4c, based on an FDA scheme, sacrifices space after each refinement.

The main challenge when selecting a type of fractal is its scalability in dimension. Scalability can be measured by three properties defined in Table 1: the partition size k, the building complexity O(F ) and the data structure. According to the partition size property, the hypercube [START_REF] Demirhan | FRACTOP: A Geometric Partitioning Metaheuristic for Global Optimization[END_REF] and simplices [START_REF] Paulavičius | Simplicial Lipschitz optimization without the Lipschitz constant[END_REF] are not scalable for high dimensional search space. This non-scalability impacts the complexity of the partitioning function F .

According to the building complexity property O(F ), the Voronoï cell needs a particular close attention. Indeed, most of the algorithms used for computing the Voronoï diagram are well studied for 2D and 3D (e.g. QuickHull [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF], Fortune's algorithm [START_REF] Fortune | A sweepline algorithm for Voronoi diagrams[END_REF], Bowyer-Watson algorithm [START_REF] Watson | Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes*[END_REF]). However, they suffer from the curse of dimensionality. Hence, one has to use heuristics to create such fractals. Indeed, one cannot build the exact diagram in high dimensions [START_REF] Liu | A global optimization algorithm for simulationbased problems via the extended DIRECT scheme[END_REF][START_REF] Rushdi | VPS: VORONOI PIECEWISE SURROGATE MODELS FOR HIGH-DIMENSIONAL DATA FITTING[END_REF][START_REF] Mitchell | Spoke-Darts for High-Dimensional Blue-Noise Sampling[END_REF][START_REF] Villagran | Nonparametric Sampling Approximation via Voronoi Tessellations[END_REF][START_REF] Khachiyan | Generating All Vertices of a Polyhedron Is Hard[END_REF], and then one has to use stochastic methods to build an approximation, such as SpokeDart with hyperplane sampling [START_REF] Mitchell | Spoke-Darts for High-Dimensional Blue-Noise Sampling[END_REF]. We distinguish two approaches, the fixed (Figure 5f) and dynamic (Figure 5e) Voronoï diagrams. In the fixed Voronoï diagram, a child cell will inherit its parent's boundaries, whereas in the dynamic diagram, when building a child, the whole diagram is re-computed to consider children's boundaries. This approach implies variations of the borders of fractals. In a dynamic approach, a Voronoï cell is a very complex type of fractal, as its definition depends on varying information on other fractals. For example, when recomputing the diagram, points sampled in a dynamic Voronoï cell, can be relocated to another ones. Whereas in the fixed approach, once a fractal, i.e., a Voronoï cell is built, then its borders cannot be changed anymore. Other hypervolumes showed in Figure 5 are easier to compute, and further details can be found in [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF][5] [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF].

Furthermore, if we look at the coverage and overlap properties, one can clearly see the major drawback of hyperspheres compared to other hypervolumes. C(S, L T (d) ) and C(S l,j,i , F (S l,j,i )) measure the uncovered space. If these measures are strictly superior to 0, then the partition does not fully cover the search space or a parent fractal.

According to [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF], the inflation ratio applied on n-spheres can reduce the substantial lack of coverage by increasing the overlap. However, we have to mention the impact of the curse of dimensionality on such objects. Indeed, if we look at the Hausdorff measure of a unit n-ball, the n-volume and the n-surface tend to zero as the dimension tends to infinity. This can explain the empirical results obtained in this paper, indicating that the deeper is the decomposition tree, the lower FDA performs. According to Figure 5, 2D representations are misleading. Even if they are intuitive and allow to better understand basic principles, we cannot infer what will happen in n-dimension by only looking at 2D or 3D drawings. Thus, the choice of a fractal can be informed by carefully examining the six properties defined in Table 1.

Tree search

One can define a tree search algorithm as a function τ which selects Q unique fractals among all L T (d) according to a vector of size s = |L T (d) |:

τ : L T (d) , R s → {e1, ..., e Q } ( 12
)
where ∀q ∈ [1, ..., Q] : e q ∈ L T (d) , 1 ≤ Q ≤ s. The vector R s resulting from the scoring search component describes the quality value of each leaf.

One can instantiate tree search algorithms using an OPEN and CLOSED lists [START_REF] Dechter | Generalized best-first search strategies and the optimality of a*[END_REF]. The OPEN list contains non-explored or non-exploited fractals. The CLOSED list contains expanded fractals. Thus, the tree search component is a rule defining how to append non-expanded fractals to these lists, and how to select them. The CLOSED list is generally used to prevent the algorithm from selecting expanded fractals. a Number of centroids defined by the user.

O(F ) O(2 n ) O(n) O(n) O(n) O(2 n ) b O(n! ) Coverage of S C(S, L T (d) ) 0 0 0 > 0 0 0 Children coverage C(A, F (A)) 0 0 0 > 0 0 0 Overlap O(L T (d) ) 0 0 0 > 0 0 0
b Valid for usual algorithms, we can reduce this complexity by approximating the Voronoï diagram in high dimensions. This complexity, also depends on c, the number of centroids. But here we consider the complexity depending on the dimension n.

c It depends on the algorithm used to compute the diagram. It can be a set of vertices for the QuickHull algorithm or a set of hyperplanes for sampling methods.

When building a k-refinement, tree search algorithms are crucial and have an impact on the exploration and exploitation tradeoff. In the DIRECT algorithm, this issue is tackled by the Potentially Optimal Rectangle (POR), and many variations of DIRECT are based on this selection strategy [START_REF] Gablonsky | A Locally-Biased form of the DIRECT Algorithm[END_REF] [START_REF] Finkel | An Adaptive Restart Implementation of DIRECT[END_REF]. The FDA algorithm uses a sorted depth first search, called Move-up, which favors the selection of deep fractals. The FDA algorithm focuses on the exploitation phase applied at the last level D. Such strategy impacts the capacity of the algorithm to efficiently explore the search space.

Moreover, Breadth First Search (BFS) is totally useless in our case. Indeed, all fractals at a certain level will be decomposed before selecting fractals of the next level. Thus, Depth First Search (DFS) can be seen as a greedy exploitation only, whereas Breadth First Search (BFS) can be qualified as a greedy exploration only algorithm [START_REF] Even | Graph Algorithms[END_REF]. Both strategies are uneffective, as the scores given to fractals do not impact the selection. We lose the notion of hierarchy between fractals.

A tree search algorithm can be characterized by its balance between exploration and exploitation of the tree. Thus, one efficient and easy to use algorithm, which can replace DFS or BFS, is the Best First Search [START_REF] Dechter | Generalized best-first search strategies and the optimality of a*[END_REF]. Some of these algorithms are stochastic (e.g. Epsilon Greedy Search, Diverse Best First Search [START_REF] Imai | A novel technique for avoiding plateaus of greedy best-first search in satisficing planning[END_REF]) or, others, allow having the hand on the exploration-exploitation tradeoff (e.g. Cyclic best First Search [START_REF] Morrison | Cyclic best first search: Using contours to guide branch-and-bound algorithms[END_REF]). The notion of sacrifice defined in Theorem 2 also concerns tree search algorithms and pruning techniques. For time and memory complexity reasons, it can be necessary to prune some leaves. For instance, in the Beam Search algorithm [START_REF] Frohner | Parallel Beam Search for Combinatorial Optimization (Extended Abstract)[END_REF], only a given number of leaves are stored. explored the entirety of each level before tackling the next one, no matter the hierarchy. Whereas DFS is focused on the deepest fractal.

Scoring fractals

The tree search component needs a heuristic value determining how promising a fractal is. We define a scoring method γ which takes a set of solutions restricted to a fractal and assigns a quality value:

γ : •[l, j, i], P, f | •[l,j,i] (P ) → R (13)
The scoring method γ takes a leaf •[l, j, i], a sample of solutions P restricted to this leaf, and their corresponding objective values f | •[l,j,i] (P ). It returns a score within R which defines the quality of •[l, j, i].

For DIRECT and SOO, γ returns f (x), with x the center of the fractal. In FRACTOP, the Belief is used. It depends on a fuzzy measure computed with the best evaluation found so far, and sampled points within the fractal. This score has the particularity to inherit a part of the parent's score.

In FDA, the algorithm maximizes the distance-tothe-best solution found so far among all sampled solutions.

We have noticed that the combination between τ and γ is essential, and has different purposes concerning the exploration and exploitation tradeoff.

Finally, τ and γ can use additional information, such as a measure of the size of a fractal. For example, the σ function in DIRECT, with σ 2 or σ ∞ , measures the size of hyperrectangles according to their level and the length of their longest side [START_REF] Gablonsky | A Locally-Biased form of the DIRECT Algorithm[END_REF]. We can notice that for a fractal for which its children are smaller and of identical size (regularity), the level of a fractal can be considered as a measure of its size. This is the case for FDA, but not in DIRECT and for Voronoï fractals.

Exploration and exploitation strategies

The exploration Explore and exploitation Exploit search components can be defined by the following functions applied on a fractal α:

Explore : α → P, f | α (P )

Exploit : α → P, f | α (P ) (14)
To assess the potential of a fractal, one has to sample relevant solutions from it. This is the purpose of the Explore function. For instance, FRACTOP uses a genetic algorithm, whereas in FDA, three fixed solutions are computed.

Sampling in hyperrectangles or hypercubes is a simple procedure. One can apply all sampling methods or metaheuristics using upper and lower bounds. For example, low discrepancy sequences such as Sobol, Halton, and Kronecker methods can be used [START_REF] Drmota | Sequences, Discrepancies and Applications[END_REF]. Sampling in a hypersphere requires a few tricks to satisfy the equation of a n-ball. The Box-Muller method can be a solution [START_REF] Voelker | Efficiently sampling vectors and coordinates from the n-sphere and n-ball[END_REF][START_REF] Muller | A note on a method for generating points uniformly on n-dimensional spheres[END_REF][START_REF] Harman | On decompositional algorithms for uniform sampling from nspheres and n-balls[END_REF]. However, sampling inside a Voronoï cell (i.e. polytope) is a complex procedure. One could use methods to approximate the Lebesgue measure [START_REF] Ge | A fast and practical method to estimate volumes of convex polytopes[END_REF], hit-and-run sampling [START_REF] Corte | Novel matrix hit and run for sampling polytopes and its GPU implementation[END_REF], hyperplanes sampling [START_REF] Mitchell | Spoke-Darts for High-Dimensional Blue-Noise Sampling[END_REF], or MCMC sampling [START_REF] Chen | Fast mcmc sampling algorithms on polytopes[END_REF].

For active algorithms, such as metaheuristics (e.g. local search, evolutionary algorithms, swarm optimization), one must adapt the search operators (e.g. neighborhood, mutation, crossover, velocity update) to the type of fractal.

Finally, if the algorithm lacks of exploitation, one could apply an exploitation algorithm within leaves. In our framework, this algorithm is named Exploit. Such a function starts from a solution within a leaf fractal and if needed can ignore its boundaries to converge towards an optimum. This is the case for FRACTOP and FDA algorithms.

It is a thorough task to select an exploration and an exploitation search components. In our framework, the exploitation search component is exclusively applied to fractals of maximum depth to emphasize the search around promising areas. Consequently, it is appropriate to allocate larger budget to the exploitation compared to the exploration strategy within each fractal. However, according to the cost of the objective function, the maximum depth of the tree or the partition size (i.e. number of children per fractal), one should carefully choose the budget of the exploration and exploitation. Indeed, a budget that is too high for the exploration search component, can result in a low exploitation phase. Whereas, a budget that is too low for the exploration phase, can result in expensive exploitation in low confidence areas.

Extension of popular decomposition-based algorithms

This section instantiates some popular fractalbased optimization algorithms to our algorithmic framework. Moreover, it shows how one can extend these algorithms using various search strategies for the different components.

FRACTOP, FDA and PolyFRAC: Fractal-based optimization algorithms

FRACTOP is one of the first algorithms based on fractal decomposition [START_REF] Demirhan | FRACTOP: A Geometric Partitioning Metaheuristic for Global Optimization[END_REF]. It uses hypercubes to decompose the search space, a genetic algorithm to explore each fractal, and simulated annealing to exploit a promising leaf fractal. The algorithm implements a fuzzy measure, called belief, as a scoring method. One major drawback of this algorithm is its poor scalability in high dimension (see Table 1).

The FDA algorithm partly solves the curse of dimensionality problem of FRACTOP [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF]. Instead of a hypercubes-based decomposition, it uses if min(values) < x then 27:

x ← min(values) current ← τ (OPEN, scores) return x hyperspheres. By using such fractals, the decomposition has a lower complexity, but at the cost of overlapping fractals due to an inflation ratio (see Table 1). This ratio partially reduces the lack of space coverage implied by hyperspheres decomposition. The exploration component, called Promising Hypersphere Selection (PHS), computes three points: the center of the hypersphere and two opposite points equidistant to the center. The heuristic, used to score a fractal, is the distanceto-the-best solution found so far. Finally, a leaf at the maximum depth level of the tree, is exploited with an Intensive Local Search (ILS), which is a coordinate descent algorithm with adaptive step size.

The polyFRAC algorithm is a modification of FDA [START_REF] Khodabandelou | H -polytope decomposition-based algorithm for continuous optimization[END_REF]. Rather than using hyperspheres, polyFRAC takes advantage of H-polytope fractals. Since, finding the vertices (i.e. n-faces) of a Hpolytope is a complex procedure, thus polyFRAC approximates them. The exploration, exploitation, and scoring strategies are similar to those of FDA.

Thanks to our generic and flexible framework, we have extended and improved those algorithms by replacing the tree search component of FDA (Move-up procedure), by a beam search2 , we called this extension FDABs. We also tested FDA and FDABs using deep trees 3 , we named them FDAD and FDADBs. Additionally, we replaced the distance-to-the-best (DTTB) scoring:

f (x)
∥x-Best∥ by a centered version of it (C-DTTB).

The measure is centered on the best solution found so far, f (x)-f (Best) ∥x-Best∥ , where x is a sampled point and Best the best solution found so far. As we select the highest ratio, on the first version, the ratio associated to a better solution can be considered as non-promising if the point is far from Best, whereas a bad solution can be considered as good if it is close to Best. We named this version, FDAC.

Direct: Dividing Rectangles

The DIRECT algorithm is initially a modification of the Schubert's algorithm [START_REF] Shubert | A Sequential Method Seeking the Global Maximum of a Function[END_REF]. It assumes that the objective function is Lipschitz-continuous, with a positive constant K:

|f (x) -f (y)| ≤ K∥x -y∥, ∀x, y ∈ S
Thus, there are several advantages of such an assumption. Indeed, it allows to easily prove convergence towards a global optimum. Moreover, few hyperparameters have to be set, and the algorithm is deterministic.

However, Lipschitzian optimization has several drawbacks. The constant K must be known; there is a poor balance between exploration and exploitation, and the computation complexity is subject to the curse of dimensionality. The algorithm computes a solution at the center of a hyperrectangle before decomposing. At each iteration, the algorithm uses a series of trisections on the longest sides to subdivide the search space into smaller hyperrectangles. DIRECT introduces the concept of Potentially Optimal Rectangle (POR), which is a selection criterion that determines if a fractal is promising [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF]. Many extended versions of DIRECT have been proposed in the literature. These can be considered within our five search components:

• Fractal component: BIRECT (bisection) [START_REF] Paulavičius | Global optimization based on bisection of rectangles, function values at diagonals, and a set of lipschitz constants[END_REF],

eDIRECT (Voronoï cell) [START_REF] Liu | A global optimization algorithm for simulationbased problems via the extended DIRECT scheme[END_REF], DISIMPL (simplices) [START_REF] Paulavičius | Simplicial Lipschitz optimization without the Lipschitz constant[END_REF].

• Tree search component: Pareto-Lipschitzian optimization [START_REF] Mockus | On the Pareto Optimality in the Context of Lipschitzian Optimization[END_REF], DIRECT Restart [START_REF] Finkel | An Adaptive Restart Implementation of DIRECT[END_REF].

• Exploitation component: DIRMIN [START_REF] Liuzzi | A DIRECT-based approach exploiting local minimizations for the solution of largescale global optimization problems[END_REF].

These modifications try to overcome some DIRECT drawbacks, such as lower performances in high dimensions, low convergence rate when trapped by local optima, and a lack of a local optimizer. In this paper, we introduce another extension of DIRECT, called DIRECTBs, where we replace the POR tree search, by a Beam search. Furthermore, we have implemented under our framework two well-known versions of DIRECT, which are locally biased DIRECT (DIRECTL) [START_REF] Gablonsky | A Locally-Biased form of the DIRECT Algorithm[END_REF], and DIRECT Restart (DIRECTR) [START_REF] Finkel | An Adaptive Restart Implementation of DIRECT[END_REF].

SOO: Simultaneous Optimistic Optimization

DOO and SOO claim to be a generalization of the DIRECT algorithm. These algorithms make a strong assumption on the existence of a semimetric l. This assumption simplifies the Lipschitzcontinuous property by assuming a local smoothness around the global optimum x [START_REF] Munos | Optimistic Optimization of a Deterministic Function without the Knowledge of its Smoothness[END_REF] :

f (x) -f (x) ≤ l(x, x), ∀x ∈ S
DOO is used when l is known; otherwise, SOO is more adapted. The strength of these algorithms, is their low number of parameters and the proof of a convergence bound. Both algorithms are deterministic. At each iteration and at each level of the partition tree, the best fractal is selected according to the evaluation of a representative solution inside it (e.g. center). Here, the balance between exploration and exploitation relies on the tree search algorithm. In addition, a stochastic version called Sto-SOO has been designed for noisy loss function, where each fractal has to be evaluated multiple times [START_REF] Valko | Stochastic Simultaneous Optimistic Optimization[END_REF]. In the same way as DIRECTBs, we replaced the SOO tree search algorithm, by a Beam Search component, denoted by SOOBs algorithm.

8 Experimental setup

Algorithms selection

We have experimented more than 20 extensions of the presented popular algorithms. We tried different tree search algorithms, such as Best First Search and Cyclic Best First Search, different beam lengths for Beam Search, different maximum tree depth and scoring methods. In this paper, we have selected the most relevant ones. We focus on deterministic fractal-based decomposition algorithms. Hence, only one run on each function is necessary to evaluate their performances. We have excluded some popular stochastic algorithms such as FRACTOP, Sto-SOO and some versions of DIRECT, such as eDIRECT [START_REF] Liu | A global optimization algorithm for simulationbased problems via the extended DIRECT scheme[END_REF] or glcCluster [START_REF] Rios | Derivativefree optimization: a review of algorithms and comparison of software implementations[END_REF]. Even that some designed algorithms have been parallelized, a focus is made on a sequential version of the algorithms. Indeed, some parallel algorithms may have different behaviors in terms of convergence and search time. For a fair comparison, all the studied algorithms have been implemented under the Zellij framework.

The first aim of this comparison is to show that we can instantiate several popular algorithms within the Zellij framework. One can also evaluate their scalability and their sensitivity according to the different search components: fractal, tree search, scoring, exploration and exploitation components. Thus, we have instantiated FDA [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF], DIRECT [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF], DIRECTL [START_REF] Gablonsky | A Locally-Biased form of the DIRECT Algorithm[END_REF], DIRECTR [START_REF] Finkel | An Adaptive Restart Implementation of DIRECT[END_REF] and SOO [START_REF] Munos | Optimistic Optimization of a Deterministic Function without the Knowledge of its Smoothness[END_REF]. We resumed all 11 implemented algorithms and their different search components in Table 2. An unsought property of our selection is that all algorithms sample points in an axis-aligned fashion. This is a behavior that we must keep in mind when selecting benchmark functions.

Benchmark selection

For the selection of benchmarks, five essential properties have been considered. Functions must be shifted or rotated, or at least not have their global optimum located at the center of the search space. Indeed, for some well-known benchmark functions (e.g. Rastrigin, Ackley, Sphere), some selected algorithms (e.g. DIRECT, FDA, SOO) can directly place a point, or the center of a fractal onto the global optimum. Moreover, some benchmark functions must be multimodal, so to trap the algorithms into local optima, and thus evaluate their exploration capability. The benchmark functions must be non-separable, meaning that all dimensions have to be optimized to find the global optimum. This prevents algorithms from taking advantage of their axis align property to optimize only one or a few dimensions. Moreover, functions must be non-symmetrical, as some algorithms use symmetry to sample points. Therefore, we decided to select the CEC2020 benchmark [START_REF] Liang | Problem Definitions and Evaluation Criteria for the CEC 2020 Special Session on Multimodal Multiobjective Optimization[END_REF]. We consider all functions as black box, and we reduce the maximum number of calls to the objective function at 5000n, as in [6][15]. We also tried to reproduce the computational results obtained in [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF], on the SOCO2011 single objective benchmark. However, we select the first six functions, as the others are not shifted nor rotated, and their global optimum are located on (0, 0, ..., 0).

In our results' analysis, we use the two-sided Wilcoxon signed-rank test and corresponding mean ranks for each of the 11 algorithms. We performed this test on each benchmark. Thus, we were able to see their reliability. We take an error rate α = 0.05 for the statistical test. Tested dimensions are 10, 15, 20, 30, 50, 100, which are the dimensions that both benchmarks have in common. For the Wilcoxon test, the two hypotheses are:

• H0: The two samples come from the same distribution.

• H1: The two samples come from different distributions. a The maximum depth is set to 600, as the maxdeep variable in the original FORTRAN implementation. In the original code, a maxdiv variable, set to 3000, limits the number of successive decomposition of a fractal. In Zellij, when the difference between a lower and an upper bound is inferior to a value ϵ set to 1e -13, the fractal cannot be decomposed anymore.

b The beam length was set to 600.

Results analysis and discussion

First, all results are based on the gaps between the best solution found by each algorithm and the known global optimum of each function. Figure 7 shows the mean of the ranking of all 11 algorithms on all CEC2020 benchmark functions for each dimension. Summing these means and taking the lowest results is not enough to determine if an algorithm is better than another. Indeed, giving a rank to an algorithm for a function is arbitrary. An algorithm ranked second does not mean that the first one is significantly better. We need a corresponding statistical test and p-value for each algorithm. We resumed the comparison between ranks and p-values in Figure .8. These figures can be read column by column. Each column represents comparisons of an algorithm (column label) with all other 10 algorithms (row labels). As an . The color indicates that there is no statistical evidence that FDA is worse or better than SOO. In Figure 8, Figure 9 and Figure 10, there are three color codes:

• Solid-Grey: α > 0.05, we cannot reject the null hypothesis.

• Gridded-Green: α ≤ 0.05 and rank(column i ) < rank(row j ). We can reject the null hypothesis, and the algorithm with the label of the column i has a lower rank (is better) than the algorithm of the row j.

• Dotted-Red: α ≤ 0.05 and rank(column i ) > rank(row j ). We can reject the null hypothesis, and the algorithm with the label of the column i has a higher rank (is worse) than the algorithm of the row j.

The following analysis is based on the CEC2020, as it satisfies the properties we have previously defined.

Sensitivity of dimensionality

The initial observations reveal that for dimensions from 10 to 20, FDAD and FDADBs are the worst algorithms. Moreover, for all dimensions, we can see that these two algorithms are the worst among other FDA versions. This confirms what authors in [START_REF] Nakib | Deterministic metaheuristic based on fractal decomposition for largescale optimization[END_REF] noticed, the deeper the tree, the worse the results. Which appears to be true even if we replace the tree search (FDADBs).

For dimensions from 10 to 50, it is not certain that FDA is better than other algorithms, except FDAD and FDADBs, which are always worse. But for dimensions 100, FDA is better than DIRECT, DIRECTL and DIRECTR, which was expected. Indeed, DIRECT does not scale well in dimensions. In dimension 50, one can note that FDA, DIRECT, DIRECTL, DIRECTR, SOO, DIRECTBs, FDABs are equivalent.

Concerning SOO, we can observe that it is unclear if it is working better than other algorithms besides FDAD and FDADBs. Because SOO must generate many fractals to significantly reduce the search space, SOOBs quickly fill the beam, which can explain a significant decrease in terms of performances as n increases. SOO becomes significantly less effective than FDA, DIRECTL, DIRECTBs, FDABs and FDAC in dimension n = 100.

In low dimensions, DIRECTBs does not perform better compared to other DIRECT variations. However, it scales better when dimension increases. Indeed for n = 100, DIRECTBs, DIRECTR and DIRECTL are better than the original DIRECT. It shows the effect of the tree search algorithm replacement. We can say that depending on which dimension DIRECT is applied, we should adapt some of its features to make it scale better.

Finally, the best algorithms for high dimensional continuous problems, are, FDA, FDAC, FDABs and DIRECTBs.

Sensitivity of tree search

The selected algorithms are sensitive to different parameters. For FDA, increasing D, the maximum tree depth, reduces drastically the performances. It is interesting to compare the maximum tree depth of each algorithm. DIRECT and SOO require a deeper tree to perform well, conversely to FDA. However, in FDA, modifying Move-up by a Beam search does not have a significant impact. For DIRECTL we can see that modifying the tree search, All POR by 1 POR per level and the measure of the size of a fractal (from σ 2 to σ ∞ ), makes DIRECTL better than DIRECT for higher dimensions. Nonetheless, it does not scale well. For SOO, the replacement of the tree search by a Beam Search, worsens it.

Finally, this work, empirically shows that by modifying some parts of fractal decomposition-based algorithms, we can obtain different behaviors and scalability. The most prominent example is the comparison between DIRECT, DIRECTL and DIRECTBs. DIRECT does not scale well for dimensions 100, while other versions can perform better.

Sensitivity of scoring

FDAC and FDA are similar for all dimensions. We noticed that the modification of the scoring method did not have a significant impact on the algorithm's behaviors. It is not clear which part of the algorithm makes it work. We have changed the scoring and tree search components, but we did not notice any change in terms of efficacy, except for FDAD and FDADBs. The step size of ILS depends on the radius of the hypersphere in which it is applied, this can explain bad performances when we increase the maximum depth. Future works might focus on the exploitation phase, and its consequences on the performances of FDA.

Sensitivity of instances

The goal here, is to compare results obtained on the CEC2020 with the SOCO2011 benchmarks, in terms of sensitivity to the benchmark properties. The following analysis is based on Figure 9 and Figure 10.

With the 6 functions SOCO2011, the obtained results do not permit to clearly demonstrate which of FDA, FDAC, FDABs, DIRECT, DIRECTL, DIRECTR, DIRECTBs, SOO is the best. Evenif DIRECTBs seems to have better performances for n = 100 than the other DIRECT based versions. As for CEC2020, FDAD seems to be the worst algorithm for all dimensions, followed by FDADBs. Moreover, for low dimensions from 10 to 20, it is unclear, we cannot conclude that any of the algorithms performs better than the others (except FDAD and FDADBs). For dimensions 30 to 100, DIRECTBs, is better than all other DIRECT versions.

We even obtain different scalability. According to SOCO2011, the best algorithm to use in very high dimension (n = 100) is DIRECTBs, whereas CEC2020 shows that FDA, FDAC, FDABs and DIRECTBs are suitable solutions for this dimensionality.

Selecting all 10 functions of the SOCO2011 is not relevant, as F 7 to F 10 have their global optimum located at (0, 0, ..., 0). But here we want to show that when selecting a benchmark we have to study behaviors of the selected algorithms and properties of the benchmark to avoid worthless comparisons.

We can notice that by adding 4 new functions that are obviously easy to solve as, DIRECT, DIRECTL, DIRECTR, SOO, DIRECTBs, Fig. 7: Mean ranks on CEC2020 benchmark. SOOBs find instantly the global optimum, does not impact the results. Except FDAC which becomes better than FDADBs and FDAD in dimensions n = 20, n = 30 and n = 50. By adding these worthless functions to the benchmark and seeing no significant difference, we can suppose that SOCO2011 is not sufficiently robust to assess performances of selected algorithms. Therefore, we must question the relevance of this benchmark to compare decomposition-based algorithms.

Real-world problem: Sharpe ratio maximization

Portfolio optimization is a problem coming from the financial and investment management fields. It aims to find a portfolio of assets that offers the best possible trade-off between expected return and risk. The primary goal of portfolio optimization is to help investors make informed decisions about how to allocate their investments among different assets or securities. A portfolio is made of various assets such as bonds or real estate. These are defined by different levels of risks and returns. The expected return is the average return one can expect from a given portfolio over a specific period.

Let's consider w the weight vector of assets and µ the expected return for each one. Then the portfolio expected return is equal to R p = w T µ. Weights w are subject to the following constraints: n i=0 (w i ) = 1 and w i ≥ 0. Here the goal is not to tackle constrained optimization, so we do (w i ) = 1. However, when weights w i are passed to the objective functions, these are normalized to follow the constraint. So the search space is a unit hypercube of dimension n, where n is the number of assets, and the bounds are [0, 1] n .

One way to select the best portfolio, is to maximize the Sharpe ratio:

S p = Rp-R f σp
With, S p the Sharpe ratio, R p the expected return of a portfolio p, R f the risk-free rate, which is the expected minimum return on a zero risk investment, and σ p the standard error of p's expected return. A negative So Sharpe ratio optimization is a high dimensional and scalable problem, where some dimensions have to be set at 0 (no investment for a given asset). The dataset contains stocks from the New York Stock Exchange and Nasdaq Stock Market (SP500) from January 7, 2000, to December 29, 2017 [START_REF] Leung | Datasets for Portfolio Optimization[END_REF]. The base Sharpe ratio for comparison Results in Table 5 show that FDA-based algorithms can find the same optimum as the Efficient Frontier algorithm. Even FDAD and FDADBs were able to optimize this problem, whereas they were poorly ranked by the CEC2020 experiments.

Other fractal-based decomposition algorithms are significantly less efficient. We can state that in very high dimension, here n = 356, a local optimizer appears to be needed. Relevant computational results have been obtained comparing popular and extended deterministic algorithms according to the search components and their parameters. Their sensitivities to the dimension of the problem and the benchmark instances have also been analyzed.

The proposed framework opens a door for developing and analyzing various search components. We need larger and well-defined benchmarks with more functions so that we can have a sufficiently large sample to compute significant statistical analysis.

Future works will focus on the extension of the fractal-based decomposition algorithmic framework for combinatorial and mixed optimization problems, containing various types of variables (e.g. discrete, continuous, categorical). We will also investigate other geometrical fractal objects such as Voronoï fractals and Latin hypercubes, which implies stochasticity. We will have the opportunity to experiment the algorithms for noisy functions and multi-objective optimization problems.

From a long-term perspective, we will tackle the massive parallelization of fractal-based decomposition algorithms on large scale heterogeneous architectures including multi-cores and accelerators such as GPU, towards Exascale architectures.

Because our algorithmic framework is made of five distinct and independent components, one could imagine an search for the best combinations of them. Thus, we could figure out autonomous fractal decomposition-based algorithms dedicated to specific problems. In terms of application perspective, some computational experiments are under development on solving high-impact optimization problems such as hyperparameter optimization and automated design of deep neural networks.

Fig. 1 :

 1 Fig. 1: The workflow of the Zellij framework. In blue, the five search components. In orange, the tests made at each iteration.

Fig. 2 :

 2 Fig. 2: Example of a 2-refinement of depth 4, with a bisection along the longest side.

Fig. 3 :

 3 Fig. 3: 2-dimensional illustration of covering and overlap applied on a FDA scheme.

  Fig. 4: 2-dimensional illustrations of the three levels of sacrifice of the fractal search component. The hatched-magenta area corresponds to uncovered space resulting from successive decomposition: C(S l,j,i , F (S l,j,i ))

Fig. 5 :

 5 Fig. 5: Various examples of fractals in 2 dimensions: (a) hypercubes, (b) bisections, (c) trisections, (d) hyperspheres, (e) dynamic Voronoï, (f) fixed Voronoï.

Figure 6 Fig. 6 :

 66 Figure 6 illustrates how the tree search algorithm impacts the behaviors of fractal decompositionbased algorithms. One can clearly see the difference between depth and breadth first search. BFS

28 : 30 : 31 :

 283031 Append leaf to CLOSED 29: index ← Index of leaf in OPEN Remove element at index from OPEN Remove element at index from scores 32:

(a) n = 10 (Fig. 8 :

 108 Fig. 8: Ranks and Pair Wise Wilcoxon test comparison on CEC2020. Solid-Grey: Statistically unsignificative (α > 0.05). Gridded-Green: Better. Dotted-Red: Worse

Fig. 9 :

 9 Fig. 9: Ranks and Pair Wise Wilcoxon test comparison on 6 functions SOCO2011. Solid-Grey: Statistically unsignificative (α > 0.05). Gridded-Green: Better. Dotted-Red: Worse

(a) n = 10 (Fig. 10 :

 1010 Fig. 10: Ranks and Pair Wise Wilcoxon test comparison on 10 functions SOCO2011. Solid-Grey: Statistically unsignificative (α > 0.05). Gridded-Green: Better. Dotted-Red: Worse

  10 Conclusions and perspectivesIn this paper, we propose a unified and flexible algorithmic framework for fractal-based decomposition algorithms. This algorithmic framework allows modeling decomposition-based algorithms according to five search components: geometrical fractal, tree search, scoring, exploration, and exploitation. A Python package named Zellij has been developed and is available on GitHub 4 . Thanks to this framework, we can model various decomposition-based algorithms such as DIRECT, SOO, FRACTOP, FDA and much more. Furthermore, the modular programming standard used in Zellij, allows to prototype new algorithms quickly.

Table 1 :

 1 Example of fractals and their properties

		S l,j,i	n-cube	Trisection Bisection	n-sphere	Voronoï	Simplex
	Partition size	k	2 n	3	2	2n	c a	n!
	Partition
	building
	complex-
	ity

Table 2 :

 2 Instantiated algorithms using Zellij

	Algorithm	Fractal	Tree Search	Explor.	Exploit.	Scoring	Depth	Source
	FRACTOP	n-cube	Best First Search	GA	SA	Belief	4	[5]
	FDA	n-sphere	Sorted DFS	PHS	ILS	DTTB	5	[6]
	DIRECT	Trisection	All POR	Center	∅	∅	600 a	[3]
	DIRECTL	Trisection	1 POR per level	Center	∅	∅	600 a	[26]
	DIRECTR Trisection	Adaptive POR	Center	∅	∅	600 a	[27]
			Best
	SOO	Trisection	fractal at	Center	∅	∅	600	[12]
			each level
	FDABs	n-sphere	Beam Search b	PHS	ILS	DTTB	5	This work
	DIRECTBs Trisection	Beam Search b	Center	∅	∅	600 a	This work
	SOOBs	Trisection	Beam Search b	Center	∅	∅	600	This work
	FDAC	n-sphere	Sorted DFS	PHS	ILS	C-DTTB	5	This work
	FDAD	n-sphere	Sorted DFS	PHS	ILS	DTTB	10	[6]
	FDADBs	n-sphere	Beam Search b	PHS	ILS	DTTB	10	This work

Table 3 :

 3 Functions of the CEC2020 benchmark[START_REF] Liang | Problem Definitions and Evaluation Criteria for the CEC 2020 Special Session on Multimodal Multiobjective Optimization[END_REF]

	#	Function	Shifted	Rotated	Unimodality	Separability	Symmetrical

Table 5 :

 5 Portfolio optimization of SP500

	Algorithm	Best Sharpe ratio
	FDA	6.14E+00
	DIRECT	2.45E+00
	DIRECTL	2.64E+00
	DIRECTR	3.78E+00
	SOO	2.89E+00
	DIRECTBs	3.94E+00
	FDABs	6.14E+00
	FDADBs	6.14E+00
	FDAC	6.14E+00
	FDAD	6.14E+00
	SOOBs	2.90E+00
	Efficient Frontier a	6.14E+00

a The problem is converted into a convex one within the Mean-Variance framework.

The Zellij software is available under GitHub https://github.com/ThomasFirmin/zellij.

The beam length was set to

[START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF] The maximum depth was set to 5 for FDA and FDABs, and to 10 for FDAD and FDADBs

https://github.com/ThomasFirmin/zellij

Acknowledgment

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr). This work has been supported by the University of Lille, the ANR-20-THIA-0014 program AI PhDLille and the ANR PEPR Numpex.

Data availability statement

The datasets generated during the current study, containing raw data, summaries, and statistics of all experiments, are available free of charge at https://doi.org/10.57745/0JEUEK. The code used for experiments is available on GitHub at https://github.com/ThomasFirmin/fdb zellij exp.

Yes Yes

No No No


A Brun

M Greusard

J M Reynes

M Grenier

J Bamoulid

P Giraudoux

Q Lepiller

C Chirouze

K Bouiller

B Bailly

Dr Bailly Benoît
email: benoit.bailly@hopitaux-jura.fr

Description of an outbreak of hemorrhagic fever with renal syndrome in the southern Jura Mountains, France, in 2021

Keywords: Hantavirus, Puumala virus, hemorrhagic fever diuretics, renin angiotensin aldosterone system inhibitors, and nephrotoxic antimicrobial drugs (aminoglycosides, glycopeptides)

  L'archive ouverte pluridisciplinaire

Introduction

Hemorrhagic fever with renal syndrome (HFRS) is a zoonotic disease caused by various species of hantaviruses. Puumala orthohantavirus (PUUV) is the most common species detected in France, transmitted by the rodent Myodes glareolus (bank vole) [START_REF] Marc Reynes | SURVEILLANCE OF HUMAN HANTAVIRUS INFECTIONS IN METROPOLITAN FRANCE, 2012-2016[END_REF]. Clinical presentation is often benign, although half of cases require hospitalization [START_REF] Makary | Disease burden of Puumala virus infections, 1995-2008[END_REF]. An epidemic of HFRS occurred in the southern Jura Mountains (France) in 2021. The aim of our study was to accurately describe this outbreak.

Method

We included all laboratory-confirmed HFRS cases reported during the epidemic period from April to September 2021 in the three hospitals with emergency services located in the cities of Lons-le-Saunier, Pontarlier, and Saint Claude.

Laboratory-confirmed HFRS cases were suspected HFRS cases firstly tested positive using the rapid tests Reascan and Reagena POC Puumala IgM (Reagena, Finland; performed at the hospital), and secondarily tested positive for both IgM and IgG using enzyme-linked immunosorbent assays (ELISA) and immunofluorescence assay (IF), and/or tested positive using real-time RT-PCR detecting PUUV (performed at the French National Reference Center for Hantavirus -Hantavirus NRC, Institut Pasteur, Paris).

Patients were interviewed by telephone. Clinical and biological information were retrospectively collected from the medical record. We assessed the differences between hospitalized and non-hospitalized patient using Fisher's exact test for qualitative variables and Mann-Whitney test for continuous variables, as appropriate. All patients gave informed consent before participation. Parental consent was requested for minor patients. The study was approved by the Ethics Committee of Bourgogne Franche-Comté university (CERUBFC-2021-034).

Results

During the epidemic period, 112 patients had a positive or equivocal rapid test (performed in hospital laboratories), of which 110 were secondarily confirmed by the French NRC. Among these patients, 75 were hospitalized (hospitalization rate: 68.2%). Two hospitalized cases and 18 cases among the 35 non-hospitalized cases were excluded due to non-consent.

Lastly, 90 patients were enrolled in the study: 73 hospitalized and 17 non-hospitalized patients. Exposure factors to rodent excreta reported by patients were activities in uninhabited premises (barns, cottages, attics) for 35 patients (38.9%), activity in the forest for 31 patients (34.4%), and rodent contact for 11 patients (12.2%).

Median age of the 90 patients was 43.5 years, ranging from 3 to 89 years, and 68.9% were males. Symptoms were fever for 77 patients (85.6%), myalgia for 63 patients (70%), headaches for 49 patients (54.4%), abdominal pain for 43 patients (47.8%), and transient myopia for 18 patients (20%). Forty-one (45.6%) patients had a combination of fever, elevated serum creatinine value (>100 µmol/L), and thrombocytopenia (platelet count <150 G/L) at the time of hantavirus serology testing.

The detailed clinical and laboratory characteristics of hospitalized and non-hospitalized patients are shown in Table 1. Hospitalized patients were older than non-hospitalized patients (median age of 47.0 [35.0; 64.0] years vs. 36.0 [31.0; 41.0] years, respectively, p=0.046), and more frequently had at least one pre-existing comorbidity (31.5% vs. 11.8%, p=0.14). Transient myopia was only reported in hospitalized patients (20% vs. 0%, p=0.02).

The mean duration of hospitalization was 6.48 days (range 1-23). Indications for hospitalization were acute renal failure whether or not combined with fever for 46 patients (63.0%), isolated fever without diagnosis for 20 patients (27.4%), abdominal pain or headache refractory to analgesics for five patients (6.8%), and symptomatic bradycardia for two patients (2.7%). Some patients had multiple indications for hospitalization. Forty (44.4%) patients underwent medical imaging before hantavirus diagnosis. An abdominal computed tomography (CT) scan or abdominal ultrasound was performed for kidney injury or/and abdominal pain for 25 and 17 patients, respectively. Cerebral CT scan and cerebral magnetic resonance imaging were also carried out for severe headache in 10 and 3 patients, respectively. Moreover, lumbar puncture was performed for febrile headache in three patients. Finally, CT scan, abdominal ultrasound, cerebral magnetic resonance imaging, and lumbar puncture showed no abnormality related to the symptoms. Analgesic treatment, parenteral hydration, and biological monitoring were the only treatments performed in 69 (94.5%) of the hospitalized patients. During hospitalization, 21 patients (28.8%) presented a severe acute kidney injury (AKI) (plasma creatinine >353.6 µmol/L), but no patient developed severe hyperkalemia (serum potassium level ≥6.0 mmol/L) or severe metabolic acidosis (arterial blood pH ≤7.20). The median of the maximum plasma creatinine among patients with severe AKI was 549 (interquartile range: 397-621 µmol/L). Only one patient was dialyzed for two days. He had severe AKI with plasma creatinine at 550 μmol/L, but no major metabolic disorder or congestive heart failure.

Among hospitalized patient with thrombocytopenia the median platelet count was 77 G/L (interquartile range: 62-106). Three patients had severe thrombocytopenia (platelet count <30 G/L), and two of them received platelet transfusions for bleeding prevention. Among them, a 78-year-old patient with a history of congestive heart disease and atrial fibrillation treated with direct anticoagulation presented sudden death after 3 days of hospitalization.

Another 78-year-old patient presented lipothymia caused by atrioventricular blocks.

Cardiovascular implantable electronic device was implanted before diagnosis of HFRS. In total, 23 patients had a severe form of HFRS or presented complications: 21 patients had severe AKI, one patient was hospitalized in the cardiologic intensive care unit for bradycardia, and one patient died. Of note, 11 of the 25 hospitalized patients (44%) with a Hentzien score ≥20 at admission developed a severe disease and 33of the 45 hospitalized patients (73%) with a Hentzien score <20 developed a benign disease [START_REF] Hentzien | Bioclinical Test to Predict Nephropathia Epidemica Severity at Hospital Admission[END_REF]. At the end of the epidemic period, our infectious diseases team proposed an algorithm, which has yet to be validated in a prospective study, to assist the ambulatory management of patients with HFRS for the next epidemic (Figure 1).

Discussion

The hantavirus epidemic in the Jura Mountains led to a high number of hospitalizations. The clinical presentation and the clinical course of our patients were comparable to those of other studies: the prominent clinical findings were the combination of fever with widespread pain (abdominal pain, myalgia, headaches) [START_REF] Hentzien | Bioclinical Test to Predict Nephropathia Epidemica Severity at Hospital Admission[END_REF][START_REF] Latus | Clinical Course and Long-Term Outcome of Hantavirus-Associated Nephropathia Epidemica, Germany[END_REF][START_REF] Bermejo | Eosinophilia during Hantavirus infection: a cohort study[END_REF]. Transient myopia was present in 20% of patients and was only present in hospitalized patients. Myopia has already been reported as a risk factor for a severe form of HFRS probably because it reflected higher tissue permeability [START_REF] Hentzien | Bioclinical Test to Predict Nephropathia Epidemica Severity at Hospital Admission[END_REF][START_REF] Hautala | Central nervous system-related symptoms and findings are common in acute Puumala hantavirus infection[END_REF][START_REF] Hautala | Prospective study on ocular findings in acute Puumala hantavirus infection in hospitalised patients[END_REF]. In our study, 23% of all patients presented severe AKI but only one patient was dialyzed. No patient presented severe metabolic disorder. In literature, the fatality rate of HFRS caused by PUUV is relatively low (<1%) and 4-6% of hospitalized patients needed transient dialysis treatment [START_REF] Settergren | Clinical aspects of nephropathia epidemica (Puumala virus infection) in Europe: a review[END_REF][START_REF] Mustonen | Kidney disease in Puumala hantavirus infection[END_REF] In France, the national health insurance does not cover hantavirus laboratory diagnosis.

Moreover, the wide clinical spectrum of HFRS make diagnosis challenging. In our cohort, the triad of fever, elevation of creatinine, and thrombocytopenia was present in less than half of the cases at the time of the first consultation in hospital. This might explain the CT scan, abdominal ultrasound, and lumbar puncture performed before diagnosis. Performing CT scans with injection of contrast media may promote the development of severe AKI in this context [START_REF] Hentzien | Bioclinical Test to Predict Nephropathia Epidemica Severity at Hospital Admission[END_REF]. During an epidemic period and in high-risk areas, rapid tests detecting Puumala IgM could be performed in patients presenting fever associated with abdominal disorder or headache in the absence of other diagnoses. Excellent performances have been reported for these tests, with analytical sensitivity and specificity reaching 100% and 99%, respectively [START_REF] Hujakka | New immunochromatographic rapid test for diagnosis of acute Puumala virus infection[END_REF] . Performing this rapid test could avoid unnecessary radiological and biological investigations.

Hentzien et al. reported a score to identify patients at high risk of severe HFRS at hospital admission [START_REF] Hentzien | Bioclinical Test to Predict Nephropathia Epidemica Severity at Hospital Admission[END_REF]. In this study, nephrotoxic drug exposure, visual disorders, microscopic or macroscopic hematuria, leukocyte count >10 G /L and thrombocytopenia <90 G/L were predictive factors of severe HFRS. In hospitalized patients, the use of the bioclinical score to predict severe acute kidney injury could be useful in reducing the length of hospitalization. In our study, the positive predictive value of this score seemed low. However, the retrospective nature of our study and the sample size do not enable us to assess this score.

In our cohort, most hospitalized patients only received parenteral hydration and analgesic treatments. Most of these patients could have benefited from ambulatory monitoring alone. Ambulatory medical management could include analgesic treatment, oral hydration, discontinuation of nephrotoxic and anticoagulant treatments. Ambulatory laboratory monitoring could include bi-weekly platelet count, plasma creatinine level and a daily diuresis. As specific treatments are lacking, hospitalization could be limited to severe HFRS or to patients with severe comorbidities such as patients on diuretics, anticoagulant treatments, and/or bradycardia treatment.

Cleaning activities in uninhabited places, activities in forest, and contact with rodents were the main sources of exposure identified in our cohort, as described in the literature [START_REF] Crowcroft | Risk factors for human hantavirus infection: Franco-Belgian collaborative case-control study during 1995-6 epidemic[END_REF][START_REF] Van Loock | A case-control study after a hantavirus infection outbreak in the south of Belgium: who is at risk?[END_REF].

To date, available data on bank vole population dynamics are insufficient to explain the time pattern of epidemic peaks [START_REF] Reil | Environmental conditions in favour of a hantavirus outbreak in 2015 in Germany?[END_REF][START_REF] Khalil | Population Dynamics of Bank Voles Predicts Human Puumala Hantavirus Risk[END_REF][START_REF] Giraudoux | Pics épidémiques de fièvre hémoragique à syndrome rénal et stabilité historique des foyers : les connaissances écologiques sont encore insuffisantes pour les comprendre[END_REF]. As a consequence, for lack of well-designed wildlife observatories, there is currently no possibility of anticipating epidemics in France, which makes it difficult to implement strong preventive actions.

Conclusion

The hantavirus epidemic in the Jura Mountains led to numerous hospitalizations, one dialysis treatment, and one death. Early diagnosis by rapid test could avoid unnecessary radiological and biological investigations.

Figure 1 .

 1 Figure 1. Clinical management proposition for non-hospitalized adult HFSR patients.

Table 1 .

 1 Clinical characteristics and laboratory results for 90 patients with acute hantavirus infection, South Jura Mountains, France, 2021

	All	Hospitalized	Non	p value
	patients	patients (n=73)	hospitalized
	(n=90)		patients
			(n=17)