## HILBERT CURVES FOR EFFICIENT EXPLORATORY LANDSCAPE ANALYSIS NEIGHBOURHOOD SAMPLING

## A PREPRINT

## Johannes J. Pienaar

Department of Computer Science University of Pretoria Pretoria, South Africa jpienaar85@gmail.com

## Anna S. Bosman

Department of Computer Science University of Pretoria Pretoria, South Africa anna.bosman@up.ac.za

August 2, 2024

## ABSTRACT

Landscape analysis aims to characterise optimisation problems based on their objective (or fitness) function landscape properties. The problem search space is typically sampled, and various landscape features are estimated based on the samples. One particularly salient set of features is information content , which requires the samples to be sequences of neighbouring solutions, such that the local relationships between consecutive sample points are preserved. Generating such spatially correlated samples that also provide good search space coverage is challenging. It is therefore common to first obtain an unordered sample with good search space coverage, and then apply an ordering algorithm such as the nearest neighbour to minimise the distance between consecutive points in the sample. However, the nearest neighbour algorithm becomes computationally prohibitive in higher dimensions, thus there is a need for more efficient alternatives. In this study, Hilbert space-filling curves are proposed as a method to efficiently obtain high-quality ordered samples. Hilbert curves are a special case of fractal curves, and guarantee uniform coverage of a bounded search space while providing a spatially correlated sample. We study the effectiveness of Hilbert curves as samplers, and discover that they are capable of extracting salient features at a fraction of the computational cost compared to Latin hypercube sampling with post-factum ordering. Further, we investigate the use of Hilbert curves as an ordering strategy, and find that they order the sample significantly faster than the nearest neighbour ordering, without sacrificing the saliency of the extracted features.

Keywords Fitness landscape analysis · sampling · Hilbert curve

## 1 Introduction

Search landscape analysis has established itself as a useful approach for understanding complex optimisation problems and analysing evolutionary algorithm behaviour [18]. Many landscape analysis techniques have been developed over the years, with the most widely used techniques including fitness distance correlation [11], local optima networks [26], and exploratory landscape analysis (ELA) [24]. When landscape analysis produces numeric outputs, the resulting feature vectors can be used as abstract representations of problem instances, where instances with similar feature vectors are assumed to fall into similar problem classes. If these feature vectors effectively capture the important characteristics of problems, they can be used as the feature component for automated algorithm design (AAD) - specifically, automated algorithm configuration and selection. A number of recent studies have achieved different aspects of AAD using landscape analysis in specific contexts [2, 13, 14, 16, 17, 32].

Landscape analysis approaches differ in terms of what they measure or predict (e.g. ruggedness, modality, presence of funnels, and so on), and also on what they produce (e.g. numerical results or a visualisation of a phenomenon) [21]. They can also be distinguished based on the scale of the analysis [27] - a global approach attempts to characterise the features of the search space as a whole, while a local approach will consider the features of the landscape in the

## Katherine M. Malan

Department of Decision Sciences University of South Africa Pretoria, South Africa malankm@unisa.ac.za

neighbourhood of solutions. Fitness distance correlation [11] and local optima networks [26] are both examples of global approaches, whereas the average length of an adaptive walk [37] is an example of a local landscape feature.

Many landscape analysis techniques that measure local features are based on samples that are spatially correlated , i.e. sequences of neighbouring solutions, as opposed to a sample of independent solutions from the whole search space. Techniques that require such sequences of sampled solutions include correlation length for measuring ruggedness [38], entropic profiles of ruggedness and smoothness with respect to neutrality [35, 36] - adapted as a single measure of ruggedness for continuous spaces [19], approximations of gradient [20], information content features [25], and measures of neutrality [34].

In the context of numerical optimisation, a wide range of landscape analysis approaches have been implemented in the R package flacco [12], which has also been re-implemented in Python ( pflacco 1 ). Over 300 landscape metrics are included in flacco , organised into 17 sets, which include the original ELA metrics [24] covering six of the feature sets, and the set of information content features [25] consisting of five metrics. In a study of the subset of the 'cheap" feature sets in flacco , Renau et al. [29] found that the two most salient (in terms of distinguishing between problems) and robust feature metrics were from the information content and ELA meta-model feature sets. Unlike the other feature sets in flacco , the information content metrics require spatially correlated samples.The saliency of the information content feature metrics and the additional requirement of neighbourhood ordering is the motivation for this study.

To date, the two main strategies proposed for generating spatially correlated samples in continuous search spaces are random walks [22] and post-factum ordering of global samples [25]. Desirable properties of a sampling strategy for estimating local features are that the solutions provide good coverage of the search space [15], successive points are positioned close to each other (compared to other solutions in the sample) to capture landscape changes in the neighbourhood, and that the process has low computational expense. The most commonly used approaches to neighbourhood sampling for numerical optimisation are uniform random sampling or Latin hypercube sampling (LHS) [28], followed by either random or nearest-neighbour ordering. These are the approaches implemented in flacco and pflacco .

In this paper, we investigate the use of Hilbert curves [10] as a new sampling strategy. Hilbert curves are fractal spacefilling curves with two desirable properties: (1) they guarantee uniform search space coverage, and (2) neighbouring points on the curve are located close to one another. We show that Hilbert curves are comparable to ordered LHS in terms of the saliency of the extracted features, but are significantly cheaper to compute. Additionally, we show that a Hilbert curve can be used to spatially order an LHS sample, resulting in significant computational gains compared to the commonly used nearest neighbour ordering method.

## 2 Hilbert curves

A space-filling curve is a surjective continuous function from the unit interval [0 , 1] to a unit hypercube [0 , 1] d . Surjectivity implies that every point in the hypercube maps to at least one point in the interval, and continuity ensures that no areas in [0 , 1] d are missed [5]. Space-filling curves are a special case of fractal curves, and are guaranteed to fill a continuous space in the limit.

Our interest in space-filling curves derives from two useful properties that they offer, namely (1) uniform coverage of a bounded d -dimensional space, and (2) the ability to provide a unique mapping between points in the d -dimensional space and points on the 1-dimensional curve [31]. Specifically, this study considers the Hilbert space-filling curve, first proposed by D. Hilbert in 1891 [10]. The Hilbert curve is defined recursively, and can be constructed through a limit process of iteration. Each successive iteration creates an approximation of the true Hilbert curve that passes through more points in the d -dimensional unit hypercube. For practical purposes, the number of iterations is chosen to be finite, and is further referred to as the order of the Hilbert curve.

Consider the construction of a Hilbert curve in 2 D . For the first iteration, a line on the closed interval [0 , 1] and a square [0 , 1] 2 are taken. Four equidistant points are selected on the line, where the starting point is at 0 and the end point is at 1 . The square is subdivided into four equal parts. Intervals of the line connecting each pair of points are then mapped onto the square such that the intervals adjacent on the line share a common edge on the square. This results in a simple U-shape, illustrated in Figure 1(a). At each subsequent iteration, the curve from the previous iteration is divided into four equal parts. Each part is then shrunk by a factor of 1/2, rotated, and repositioned such that the four curves connect at their endpoints in a U-shaped or reverse U-shaped pattern. Figure 1 shows 2 D Hilbert curves of order 1 to 3 .

The Hilbert curve is bijective, i.e. the mapping between points on the curve and a d -dimensional space is reversible. Bijectivity of space-filling curves has been useful in organising multi-dimensional data storage and retrieval systems,

1 https://pypi.org/project/pflacco/

Figure 1: Visualisation of a 2 D Hilbert curve of order p = 1 2 3 { , , } .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-8e84d053-1552-4f47-a2af-c1ae804b82e9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-3e18e7fc-1f8e-49aa-88b2-f1edf77910b3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-87f85d8b-b6fb-4f7b-a2ee-b17972a68abc.jpg)

since it allows a linear index to the data to be constructed and searched [1, 30]. Ideally, closely related data in the d -dimensional space would also be close together on the 1 -dimensional curve. Faloutsos and Roseman [6] investigated the clustering properties of Hilbert curves compared to other space-filling curves, and found Hilbert curves to be superior when measuring the average distance (along the curve) to the nearest neighbours of a given point. This finding suggests that Hilbert curves may serve as a viable alternative to nearest neighbour sorting of a sample.

For the Hilbert curve implementation in this study, we use the hilbertcurve Python package ( https://pypi.org/ project/hilbertcurve/ ), based on the algorithm presented by Skilling [33]. All code used for the experiments conducted in this study can be found at https://github.com/jpienaar-tuks/MIT807 .

## 3 Hilbert curves as samplers

Given a Hilbert curve of a particular order, the vertices can be used as a basis for a sample in the associated multidimensional space. Table 1 shows how the number of vertices grows with the order of the curve and the dimension of the search space. When sampling for ELA, it is common practice to use sample sizes of 10 2 × d to 10 3 × d [25]. Table 1 highlights in bold the order and dimension combinations where the number of vertices exceeds a common

Table 1: Exponential growth in Hilbert curve vertices with dimension and curve order. Values that exceed a value of 10 3 × d are formatted in bold.

|                     |                     | Dimension   | Dimension   | Dimension     | Dimension      | Dimension      |
|---------------------|---------------------|-------------|-------------|---------------|----------------|----------------|
|                     |                     | 2           | 3           | 5             | 10             | 20             |
| Order               | 1                   | 4           | 8           | 32            | 1024           | 1 . 05 × 10 6  |
| Order               | 2                   | 16          | 64          | 1 024         | 1 . 05 × 10 6  | 1 . 10 × 10 12 |
| Order               | 3                   | 64          | 512         | 32 768        | 1 . 07 × 10 9  | 1 . 15 × 10 18 |
| Order               | 4                   | 256         | 4 096       | 1 . 05 × 10 6 | 1 . 10 × 10 12 | 1 . 21 × 10 24 |
| Order               | 5                   | 1 024       | 32 768      | 3 . 36 × 10 7 | 1 . 13 × 10 15 | 1 . 27 × 10 30 |
| Order               | 6                   | 4 096       | 262 144     | 1 . 07 × 10 9 | 1 . 15 × 10 18 | 1 . 33 × 10 36 |
| Budget ( 10 3 × d ) | Budget ( 10 3 × d ) | 2 000       | 3 000       | 5 000         | 10 000         | 20 000         |

sample budget of 10 3 × d . To remain within a sampling budget, we sub-sample randomly from the set of vertices on a Hilbert curve with a minimum order of 3.

## 3.1 Stochastic sampling using Hilbert curves

Since the generation of a Hilbert curve is deterministic, the following two strategies were evaluated to introduce stochasticity:

- · Selecting random points along the edges of the Hilbert curve, illustrated in Figure 2(a): Given any two sequential vertices, P i and P i +1 , a new point is selected P j = rP i +(1 -r P ) i +1 , where r ∼ U (0 1) , .
- · Selecting points near vertices of the curve, illustrated in Figure 2(b): For each vertex P i , generate a new point drawn from a normal distribution centred on P i with a standard deviation σ (constrained by the bounds of the

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-6b3d7d05-e316-4135-87a6-e44a538a4e0b.jpg)

(a) Randomising points along the edges

Figure 2: Illustration in 2D of adding stochasticity to a Hilbert curve

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-5e0cfbf7-aae6-4377-803d-729a6bec21e6.jpg)

search space). The step size of the Hilbert curve before scaling is 1, so σ was empirically chosen to be 0.3 to prevent excessive potential overlap with points generated by neighbouring vertices.

Neighbouring vertices on a Hilbert curve are equidistant, so the sequence of vertices has a constant step size. The two randomisation techniques introduce variation into the step size. Figure 3 shows that the randomisation near vertices strategy produces step sizes with a more normal distribution. For all further experiments we have used the randomisation around vertices strategy.

## 3.2 Search space coverage

We now investigate the extent to which Hilbert curve sampling covers the search space compared to competing strategies such as Latin hypercube sampling (LHS). To investigate the search space coverage, the following three strategies were used to draw samples from search space [ -5 5] , d of sizes n ∈ { 100 d, 316 d, 1000 d } for dimensions d ∈ { 5 10 20 30 , , , } :

- 1. Hilbert curve: randomly select the required number of points from the curve (without replacement).
- 2. LHS: generate the required sample using a Latin hypercube design [23].
- 3. Random walk: use a simple random walk [22] with maximum step size of 1 to generate the sample.

For each sample size and dimension combination, a random uniform sample was drawn as a reference set. Thirty independent runs of each sampling strategy were implemented, and for each run, the Hausdorff distance [9] to the reference set was calculated. The Hausdorff distance as a measure of search space coverage was first proposed and investigated by Lang and Engelbrecht [15], and their methodology is followed in this study. Statistical significance tests were performed for each sampling strategy as proposed by Derrac et al. [3] and described by Lang and Engelbrecht [15].

Table 2 shows that the Hilbert curve and LHS have similar Hausdorff distances, implying a similar coverage of the search space, whereas the random walk provides the worst coverage. Lower values were achieved in most cases by the Hilbert curve sampler, confirmed by the statistical significance tests presented in Table 3. The null hypothesis is that there is no significant difference between the Hausdorff distances for each of the sampling strategies. Results show that the null hypothesis can be rejected, given the p -values and a significance level of 0.05. We therefore conclude that the Hilbert curve sampler achieves a more uniform search space coverage than LHS.

## 3.3 Computational cost

We now compare the cost of the Hilbert curve as a sampler to alternative sampling strategies. Generating the Hilbert curve is not computationally cheap, but we expect that this upfront investment will be justified if the sample is subsequently used to calculate landscape metrics that require ordered samples.

Figure 3: Effect of randomisation strategy on step size. Sample was generated from an 8 th order Hilbert curve in 2D space resulting in a total of 65 535 points. The x -axis represents the distance between two consecutive points.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-a1969abf-45f7-4e23-b7e2-358ad509c1ad.jpg)

Table 2: Average Hausdorff distances between samplers and a uniform random sample for different sample sizes and dimensions d . Best and worst distances are highlighted in cyan and red, respectively.

|   D | Sample size   |   Sample size | Hilbert curve       | Latin hypercube     | Random walk         |
|-----|---------------|---------------|---------------------|---------------------|---------------------|
|   5 | 100 × d       |           500 | 2.0642 ( ± 0.0300)  | 2.0866 ( ± 0.0242)  | 3.9168 ( ± 0.2967)  |
|   5 | 316 × d       |          1580 | 1.6065 ( ± 0.0134)  | 1.6218 ( ± 0.0121)  | 3.9127 ( ± 0.2731)  |
|   5 | 1000 × d      |          5000 | 1.2754 ( ± 0.0097)  | 1.2690 ( ± 0.0044)  | 3.8658 ( ± 0.2952)  |
|  10 | 100 × d       |          1000 | 4.9716 ( ± 0.0295)  | 5.0651 ( ± 0.0300)  | 7.6117 ( ± 0.3323)  |
|  10 | 316 × d       |          3160 | 4.3345 ( ± 0.0108)  | 4.4297 ( ± 0.0143)  | 7.5796 ( ± 0.3217)  |
|  10 | 1000 × d      |         10000 | 3.8045 ( ± 0.0064)  | 3.8781 ( ± 0.0055)  | 7.5310 ( ± 0.3035)  |
|  20 | 100 × d       |          2000 | 9.6778 ( ± 0.0265)  | 9.9578 ( ± 0.0254)  | 12.5364 ( ± 0.2787) |
|  20 | 316 × d       |          6320 | 8.9312 ( ± 0.0121)  | 9.2413 ( ± 0.0087)  | 12.5213 ( ± 0.2007) |
|  20 | 1000 × d      |         20000 | 8.2919 ( ± 0.0082)  | 8.5967 ( ± 0.0051)  | 12.5013 ( ± 0.2640) |
|  30 | 100 × d       |          3000 | 13.5314 ( ± 0.0263) | 13.8606 ( ± 0.0205) | 16.3195 ( ± 0.2825) |
|  30 | 316 × d       |          9480 | 12.7478 ( ± 0.0177) | 13.1485 ( ± 0.0065) | 16.2629 ( ± 0.2186) |
|  30 | 1000 × d      |         30000 | 12.0615 ( ± 0.0123) | 12.4913 ( ± 0.0057) | 16.3055 ( ± 0.2677) |

Table 3: Ranks achieved by each of the sampling strategies with the p -value of the significance tests. Best (lowest) values highlighted in bold.

| Sampler         | Friedman   | Friedman-aligned   | Quade   |
|-----------------|------------|--------------------|---------|
| Hilbert curve   | 1.09       | 289.23             | 1.04    |
| Latin Hypercube | 1.91       | 431.78             | 1.96    |
| Random walk     | 3.00       | 900.50             | 3.00    |
| p -value        | < 1e-5     | < 1e-5             | < 1e-5  |

To evaluate the computational cost, a performance counter 2 was used to keep track of the time to a) generate the sample and b) calculate information content metrics for both the Hilbert curve and LHS. Using the pflacco library, the information content ( ic ) metrics were calculated for each of the 24 BBOB functions defined as part of the COCO platform [7]. This was done for each dimension and sample size listed in Table 2. Note that for LHS to be used for information content metrics, the sample needs to be ordered. We evaluated both the random ordering and the nearest neighbour ordering strategies. Also note that the information content function provided by the pflacco library was modified to accept none as an ordering argument to use the pre-ordered Hilbert curve samples.

2 Python's time.perf\_counter

Figure 4(a) shows that the Hilbert curve sampler is slower than LHS to generate a sample of the same size. While the computational cost of the Hilbert curve sampler seems to grow at least quadratically, it remains acceptable even for large sample sizes (just over 1.2 seconds to generate a sample of 30 000 points). Figure 4(b) shows that for the information content metrics, Hilbert curve sampling strategy is substantially faster (especially at larger sample sizes) than LHS with nearest neighbour ordering.

## 3.4 Predictive performance of Hilbert curve samples

We have shown that Hilbert curve sampling provides an efficient method for generating spatially correlated sequences of solutions from continuous search spaces, and that these sequences provide good coverage of the search space. We now investigate whether the landscape features extracted from these samples provide good predictive performance for algorithm selection.

Following the approach by Muñoz et al. [25], we use the task of predicting the class of each BBOB function as a proxy for algorithm selection: if the landscape features of function instances can be used to discriminate between problem classes, then they should also be effective predictor variables for algorithm selection. Five target classes corresponding to the BBOB function groupings as outlined in [8] were used, namely (1) separable functions: { f 1 . . . f 5 } , (2) functions with low or moderate conditioning: { f 6 . . . f 9 } , (3) unimodal functions and functions with high conditioning: { f 10 . . . f 14 } , (4) multimodal functions with adequate global structure: { f 15 . . . f 19 } , (5) functions with low or moderate conditioning: { f 20 . . . f 24 } .

The experimental setup is as follows: 24 BBOB functions in dimensions d ∈ { 2 5 10 20 , , , } are used, with an evaluation budget of 1 000 × d . Each BBOB function has 15 predefined instances, and three random instances per function are held out as testing data 3 . For each combination of instance and dimension, 1 000 d samples are drawn from [ -5 5] , d , using three sampling strategies: Hilbert curve (HC), Latin hypercube (LH) and random walk (RW). Four feature sets from the pflacco package were selected that did not require further function evaluations, namely: dispersion ( disp ), ELA y-distribution ( ela\_distr ), ELA meta model ( ela\_meta ), and information content ( ic ).

Table 4 gives the testing accuracy of a decision tree, k-nearest neighbour, and random forest classifier for the task of predicting the function class from the landscape metrics. The scikit-learn version 1.1.3 default settings were used in all instances. Results show that a random forest model was the most effective at predicting the function class across all sampling strategies. Overall, the RW strategy emerges as the least competitive performer, while the Hilbert curve and Latin hypercube sampling strategies were comparable for all three classifiers.

3 This leave-one-instance out approach is used instead of the leave-one-problem out approach due to the small number of classes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-1f2f7c5a-99d7-4e4f-9893-e727f71c3c91.jpg)

(a) Time required to generate sample only

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-149e255b-de7a-4872-9400-7c39b5797f91.jpg)

(b) Time to generate sample and calculate information content metrics

Figure 4: Comparison of time (in seconds) to generate samples and calculate information content metrics for Hilbert curve sampling and LHS using nearest neighbour and random ordering. Trendlines are polynomials of order 2.

Table 4: Accuracy at predicting the function class using four feature sets based on samples from the three different sampling strategies

|    | Decision Tree        | k-Nearest Neighbour   | Random Forest        |
|----|----------------------|-----------------------|----------------------|
| HC | 93 . 50%( ± 5 . 62%) | 70 . 15%( ± 7 . 41%)  | 97 . 38%( ± 1 . 64%) |
| LH | 93 . 87%( ± 5 . 01%) | 72 . 31%( ± 7 . 34%)  | 97 . 64%( ± 2 . 36%) |
| RW | 88 . 78%( ± 7 . 06%) | 63 . 97%( ± 6 . 87%)  | 93 . 60%( ± 3 . 93%) |

Table 5: Testing accuracy for predicting the function class for different classifiers and different sampling strategies: Hilbert curve (HC), Latin hypercube (LH), random walk (RW), under different dimensions ( D ).

|           |    | Decision tree   | Decision tree   | Decision tree   | k-Nearest Neighbour   | k-Nearest Neighbour   | k-Nearest Neighbour   | Random Forest   | Random Forest   | Random Forest   |
|-----------|----|-----------------|-----------------|-----------------|-----------------------|-----------------------|-----------------------|-----------------|-----------------|-----------------|
| Features  | D  | HC              | LH              | RW              | HC                    | LH                    | RW                    | HC              | LH              | RW              |
| disp      | 2  | 66.67%          | 52.78%          | 36.11%          | 70.83%                | 72.22%                | 34.72%                | 72.22%          | 69.44%          | 40.28%          |
| disp      | 5  | 61.11%          | 55.56%          | 37.50%          | 61.11%                | 56.94%                | 36.11%                | 54.17%          | 69.44%          | 25.00%          |
| disp      | 10 | 55.56%          | 55.56%          | 25.00%          | 52.78%                | 65.28%                | 20.83%                | 66.67%          | 63.89%          | 26.39%          |
| disp      | 20 | 58.33%          | 65.28%          | 34.72%          | 65.28%                | 63.89%                | 31.94%                | 63.89%          | 66.67%          | 22.22%          |
| ela_distr | 2  | 58.33%          | 48.61%          | 34.72%          | 52.78%                | 55.56%                | 33.33%                | 51.39%          | 51.39%          | 31.94%          |
| ela_distr | 5  | 66.67%          | 66.67%          | 33.33%          | 58.33%                | 62.50%                | 30.56%                | 63.89%          | 68.06%          | 40.28%          |
| ela_distr | 10 | 73.61%          | 79.17%          | 40.28%          | 70.83%                | 68.06%                | 45.83%                | 77.78%          | 75.00%          | 36.11%          |
| ela_distr | 20 | 79.17%          | 75.00%          | 41.67%          | 69.44%                | 63.89%                | 43.06%                | 76.39%          | 76.39%          | 43.06%          |
| ela_meta  | 2  | 79.17%          | 79.17%          | 59.72%          | 62.50%                | 66.67%                | 48.61%                | 86.11%          | 81.94%          | 79.17%          |
| ela_meta  | 5  | 83.33%          | 86.11%          | 83.33%          | 75.00%                | 79.17%                | 59.72%                | 87.50%          | 93.06%          | 84.72%          |
| ela_meta  | 10 | 95.83%          | 90.28%          | 86.11%          | 76.39%                | 87.50%                | 52.78%                | 91.67%          | 94.44%          | 88.89%          |
| ela_meta  | 20 | 91.67%          | 88.89%          | 83.33%          | 81.94%                | 87.50%                | 50.00%                | 94.44%          | 91.67%          | 87.50%          |
| ic        | 2  | 90.28%          | 94.44%          | 81.94%          | 90.28%                | 88.89%                | 84.72%                | 90.28%          | 87.50%          | 87.50%          |
| ic        | 5  | 90.28%          | 90.28%          | 84.72%          | 84.72%                | 91.67%                | 86.11%                | 93.06%          | 93.06%          | 91.67%          |
| ic        | 10 | 87.50%          | 91.67%          | 91.67%          | 86.11%                | 93.06%                | 86.11%                | 91.67%          | 95.83%          | 94.44%          |
| ic        | 20 | 88.89%          | 84.72%          | 94.44%          | 73.61%                | 91.67%                | 88.89%                | 86.11%          | 93.06%          | 94.44%          |

Table 5 gives the performance for each of the three sampling strategies based on the four feature sets separately. An observable trend is that both the ela\_meta and ic features are very salient features and allow the classifiers to accurately discriminate between the function classes, confirming results of Renau et al. [29].

## 4 Hilbert curves as an ordering tool

In this section, we investigate the use of Hilbert curves as an ordering aid for samples generated from other sampling methodologies. A Hilbert curve of order p maps each point on a 1-dimensional curve [0 , . . . , 2 dp ] to a point in [0 , . . . , 2 p -1] d and vice versa, as illustrated earlier in Figure 1. The mapping is bijective, and allows for an efficient mapping from a point in the d -dimensional space to a point on the 1-dimensional Hilbert curve.

This d -D to 1-D mapping can be exploited, given that a Latin hypercube sample of n points divides each axis into n intervals. By selecting the appropriate curve order p , i.e. p = ⌈ log 2 ( n +1) ⌉ , we can ensure that the resulting Hilbert curve is of sufficient length to accurately map all LHS points to points on the Hilbert curve. The ordering of the points on the Hilbert curve can then be used to provide a spatially correlated ordering of the Latin hypercube sample.

Given that the current best practice of using a greedy nearest neighbour strategy to order LHS requires calculating at least half of the pairwise distance matrix between all points in the sample, the Hilbert curve ordering strategy has the potential to be significantly faster. Furthermore, the nearest neighbour strategy tends to start with relatively short step sizes, which increase as the number of unvisited points decrease. Hilbert curve ordering is likely to provide more consistent step sizes.

## 4.1 Step size consistency

Figure 5 shows a visual comparison of the Hilbert curve and nearest neighbour ordering strategies on a 2 D Sphere function sample. The underlying sample X for both strategies is identical, and was obtained using LHS. In the figures the orderings produced by the Hilbert curve (HC) and the nearest neighbour (NN) approach are depicted as a red line. Colour scale is used to indicate the fitness values of the sampled points. It is evident from Figure 5 that the maximum

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-7966a3c1-f5cf-4b5e-80ca-120357b8f6f8.jpg)

(a) Ordering using Hilbert curves

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-e3b1297c-702d-4789-a009-fe6dc42da115.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-91fc2e28-58aa-4c95-a625-605abd15c7ac.jpg)

(b) Ordering using nearest neighbour

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-91916af9-9f6e-47cc-9a1c-3ae0a5d0cf43.jpg)

(c) Step sizes using Hilbert curves

(d) Step sizes using nearest neighbour

Figure 5: Visualisation of sample ordering using Hilbert curve (left) and nearest neighbour (right) ordering strategies.

step size is smaller for HC than for NN. The step size distributions shown in Figures 5(c) and 5(d) confirm that the maximum step size of HC is lower, and show a greater skew in the step size distribution for NN.

## 4.2 Computational cost

To evaluate the computational cost of the Hilbert curve ordering, we compare it to the nearest neighbour and random (RND) ordering strategies applied to a sample generated using LHS. The Latin hypercube samples of sizes n ∈ { 100 d, 316 d, 1000 d } were drawn from search space [ -5 5] , d for dimensions d ∈ { 5 10 20 30 , , , } . For each configuration, 24 independent samples were drawn, and timed with a performance counter (Python's time.perf\_counter ).

Figure 6 shows the performance comparison between the various ordering strategies. Hilbert curve sampling as discussed in Section 3 is included for completeness, as it does not require additional ordering. From Figure 6, NN ordering performs the worst, especially for large sample sizes. Furthermore, RND sorting and Hilbert curve sampling are equally fast (red and blue curves overlap). Hilbert curve ordering produces an intermediate result, slightly slower than random ordering, but significantly faster than the NN ordering. Since the underlying sample is generated by LHS, Hilbert curve ordering retains all the benefits of LHS for ELA features insensitive to order (such as ela\_meta , or disp ).

## 4.3 Evaluation of features generated using Hilbert curve ordering

We now compare the information content features generated from a Latin hypercube sample using the HC, NN, and RND ordering. Since the feature values are evaluated rather than the computational cost, this experiment was restricted

Figure 6: Comparison of time (in seconds) to order the samples generated by the Latin hypercube sampling strategy using Hilbert curve, nearest neighbour, and random ordering. Trendlines are polynomials of order 2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-183c0798-754e-4c6d-953c-321bea2adccf.jpg)

to dimensions d ∈ { 2 5 , } and sample sizes n ∈ { 100 d, 1000 d } . A single sample set X was generated and evaluated on each of the BBOB functions to generate corresponding y = f ( X ) results. Information content metrics were then calculated on these ( X,y ) pairs using each of the ordering strategies. This procedure was repeated 30 times.

Figures 7 and 8 present initial partial information ( M 0 ) and maximum information content ( H max ) landscape features based on the entropy of the sample. We consider BBOB problems in 5 D (problem indices listed on the x -axis), and calculate feature values based on NN, RND, and HC ordering for n = 1000 d .

It is evident from Figure 7 that the HC ordering produces values that lie between RND and NN, but follow a pattern similar to that of NN. It is clear that M RND 0 converges to the same value regardless of the underlying function (confirming the results of [25]). M NN 0 and M HC 0 , on the other hand, produce values that allow for discrimination between functions. Since larger values for M 0 indicate a more rugged landscape (more changes in concavity over the length of the walk), the HC ordering would on average indicate a more rugged landscape than the NN ordering.

According to Figure 8, H RND max varies over a wider range than M RND 0 . However, the H RND max NN values still tend to cluster around the upper edge of the graph, closer to the maximum value of 1. In comparison, H max and H HC max vary over a larger part of the [0 , 1] range. Larger values for H max indicate a larger variety of objects in the landscape (pits, peaks, plateaus). Thus, HC ordering again indicates a more rugged landscape than the NN ordering.

To evaluate the relative saliency of the features produced by the various ordering strategies, permutation-based feature importance technique [4] is employed. We train a random forest (RF) classifier using information content features as input ( { ϵ s , ϵ max , ϵ 0 5 . , H max and M 0 } , as defined in [25]), and the function class as output (as defined in [25] and used in Section 3). Feature saliency is gauged by randomly permuting each input feature in turn across the dataset, and observing the corresponding reduction in RF accuracy. For this purpose, the features are grouped by the ordering strategy, and a RF is trained on 2 3 of the data, with the remaining 1 3 held out for testing.

After training, the base accuracy of the RF classifier is evaluated on the held out data. The accuracy of the RF classifier is then re-evaluated as each feature in the held-out dataset is randomly permuted in turn, and the difference to the base accuracy is recorded. This procedure is repeated 10 times for each sampling strategy. The results are shown in Table 6.

Table 6 shows that ϵ s was the most salient feature across all sample orderings, with the RF trained on the HC ordered features relying on ϵ s more heavily than the RF trained on the NN or RND ordered features. RND sample ordering made the M 0 feature near useless for the classifier, in line with the observations made in Figure 7. Conversely, RF trained on the NN ordering relied on M 0 more than RF trained on the HC ordered features. Overall, HC ordering yielded marginally better RF performance compared to the NN and RND ordering approaches.

Figure 7: M 0 in 5 D as calculated using nearest neighbour (NN), random (RND) and Hilbert curve (HC) ordering methodologies ( n = 1000 d ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-f1b8c208-ffbb-42b9-ad65-b797cb1aef9f.jpg)

Figure 8: H max in 5 D as calculated using nearest neighbour (NN), random (RND) and Hilbert curve (HC) ordering methodologies ( n = 1000 d ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00526v1-figure-c3ed5962-cd57-4329-ad07-378a4a9bf6c1.jpg)

## 5 Conclusion

This paper proposed the use of Hilbert space-filling curves in the context of optimisation problem landscape analysis for the purpose of (1) sampling the search space in a spatially correlated manner that also guarantees uniform coverage, and (2) spatially ordering samples generated using other sampling algorithms such as the Latin hypercube. Experiments were conducted to evaluate the relative computational efficiency of the Hilbert curves, as well as the saliency of the landscape features extracted using Hilbert curve sampling and Hilbert curve ordering. In the context of sampling, Hilbert curves were significantly faster than Latin hypercube sampling for the purpose of generating order-sensitive

Table 6: Impact of permuting information content features on the accuracy of a random forest (RF) classifier. Classification target was the function groupings.

| Sampler                  | LHS (HC)         | LHS (NN)         | LHS (RND)        |
|--------------------------|------------------|------------------|------------------|
| RF base accuracy         | 98.44% (±0.39%)  | 98.25% (±0.37%)  | 97.44% (±0.54%)  |
| ϵ max ( ic.eps_max )     | -8.17% (±2.22%)  | -8.46% (±1.98%)  | -14.36% (±2.76%) |
| ϵ 0 . 5 ( ic.eps_ratio ) | -20.62% (±3.28%) | -10.45% (±1.64%) | -27.82% (±3.23%) |
| ϵ s ( ic.eps_s )         | -52.40% (±2.47%) | -47.34% (±1.69%) | -39.27% (±2.56%) |
| H max ( ic.h_max )       | -8.56% (±0.87%)  | -6.14% (±0.65%)  | -15.08% (±1.76%) |
| M 0 ( ic.m0 )            | -12.02% (±1.02%) | -20.98% (±1.60%) | -0.09% (±0.52%)  |

features such as the information content metrics. Features extracted by the Hilbert curves were informative, and allowed for successful discrimination between problem classes. As an ordering tool, Hilbert curves performed significantly faster than the commonly used nearest neighbour ordering, and also yielded salient landscape features. Thus, Hilbert curves present a viable computationally efficient alternative to both Latin hypercube sampling, and nearest neighbour ordering of a sample for the purpose of landscape analysis.

## References

- [1] Abel, D.J., Mark, D.M.: A comparative analysis of some two-dimensional orderings. International Journal of Geographical Information Systems 4 (1), 21-31 (1990), https://dx.doi.org/10.1080/02693799008941526
- [2] Beham, A., Wagner, S., Affenzeller, M.: Algorithm selection on generalized quadratic assignment problem landscapes. In: Proceedings of the Genetic and Evolutionary Computation Conference. p. 253-260. GECCO '18, Association for Computing Machinery, New York, NY, USA (2018), https://doi.org/10.1145/3205455. 3205585
- [3] Derrac, J., García, S., Molina, D., Herrera, F.: A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm and Evolutionary Computation 1 (1), 3-18 (2011). doi:https://doi.org/10.1016/j.swevo.2011.02.002
- [4] Dwivedi, R., Dave, D., Naik, H., Singhal, S., Omer, R., Patel, P., Qian, B., Wen, Z., Shah, T., Morgan, G., Ranjan, R.: Explainable ai (xai): Core ideas, techniques, and solutions. ACM Comput. Surv. 55 (9) (jan 2023). doi:10.1145/3561048, https://doi.org/10.1145/3561048
- [5] Falconer, K.: Fractal geometry: mathematical foundations and applications. John Wiley &amp; Sons (2004)
- [6] Faloutsos, C., Roseman, S.: Fractals for secondary key retrieval. In: Proceedings of the Eighth ACM SIGACTSIGMOD-SIGART Symposium on Principles of Database Systems. p. 247-252. PODS '89, Association for Computing Machinery, New York, NY, USA (1989), https://doi.org/10.1145/73721.73746
- [7] Hansen, N., Auger, A., Ros, R., Mersmann, O., Tušar, T., Brockhoff, D.: COCO: A platform for comparing continuous optimizers in a black-box setting. Optimization Methods and Software 36 (1), 114-144 (2020), https://dx.doi.org/10.1080/10556788.2020.1808977
- [8] Hansen, N., Finck, S., Ros, R., Auger, A.: Real-Parameter Black-Box Optimization Benchmarking 2009: Noiseless Functions Definitions. Research Report RR-6829, INRIA (2009), https://inria.hal.science/ inria-00362633
- [9] Heinonen, J.: Lectures on Analysis on Metric Spaces. Springer Science &amp; Business Media (2001). doi:10.1007/9781-4613-0131-8
- [10] Hilbert, D.: Über die stetige abbildung einer linie auf ein flächenstück. In: Dritter Band: Analysis· Grundlagen der Mathematik· Physik Verschiedenes, pp. 1-2. Springer (1935)
- [11] Jones, T., Forrest, S.: Fitness Distance Correlation as a Measure of Problem Difficulty for Genetic Algorithms. In: Proceedings of the Sixth International Conference on Genetic Algorithms. pp. 184-192 (1995)
- [12] Kerschke, P., Trautmann, H.: Comprehensive feature-based landscape analysis of continuous and constrained optimization problems using the R-package flacco. In: Applications in Statistical Computing, pp. 93-123. Springer (2019)
- [13] Kostovska, A., Jankovic, A., Vermetten, D., Džeroski, S., Eftimov, T., Doerr, C.: Comparing algorithm selection approaches on black-box optimization problems. In: Proceedings of the Companion Conference on Genetic and Evolutionary Computation. ACM (Jul 2023), https://doi.org/10.1145/3583133.3590697

| [14]   | Kuk, J., Goncalves, R., Pozo, A.: Combining fitness landscape analysis and adaptive operator selection in multi and many-objective optimization. In: 2019 8th Brazilian Conference on Intelligent Systems (BRACIS). IEEE (Oct 2019). doi:10.1109/bracis.2019.00094                                                |
|--------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [15]   | Lang, R.D., Engelbrecht, A.P.: Decision space coverage of random walks. In: 2020 IEEE Congress on Evolutionary Computation (CEC). pp. 1-8. IEEE (2020), https://dx.doi.org/10.1109/CEC48606.2020.9185623                                                                                                          |
| [16]   | Liefooghe, A., Malan, K.M.: Adaptive landscape-aware constraint handling with application to binary knapsack problem. In: Proceedings of the Companion Conference on Genetic and Evolutionary Computation. ACM(Jul 2023). doi:10.1145/3583133.3596405                                                             |
| [17]   | Malan, K.M.: Landscape-aware constraint handling applied to differential evolution. In: Theory and Practice of Natural Computing. Lecture Notes in Computer Science, vol. 11324, pp. 176-187. Springer International Publishing (2018)                                                                            |
| [18]   | Malan, K.M.: A survey of advances in landscape analysis for optimisation. Algorithms 14 (2), 40 (2021), https://dx.doi.org/10.3390/a14020040                                                                                                                                                                      |
| [19]   | Malan, K.M., Engelbrecht, A.P.: Quantifying ruggedness of continuous landscapes using entropy. In: 2009 IEEE Congress on Evolutionary Computation. pp. 1440-1447 (2009), https://doi.org/10.1109/CEC.2009. 4983112                                                                                                |
| [20]   | Malan, K.M., Engelbrecht, A.P.: Ruggedness, funnels and gradients in fitness landscapes and the effect on PSO performance. In: 2013 IEEE Congress on Evolutionary Computation. IEEE (Jun 2013), https://doi.org/10. 1109/cec.2013.6557671                                                                         |
| [21]   | Malan, K.M., Engelbrecht, A.P.: A survey of techniques for characterising fitness landscapes and some possible ways forward. Information Sciences 241 , 148-163 (2013), https://dx.doi.org/10.1016/j.ins.2013.04. 015                                                                                             |
| [22]   | Malan, K.M., Engelbrecht, A.P.: A progressive random walk algorithm for sampling continuous fitness landscapes. In: 2014 IEEE Congress on Evolutionary Computation (CEC). pp. 2507-2514. IEEE (2014), https://dx.doi. org/10.1109/CEC.2014.6900576                                                                |
| [23]   | McKay, M.D., Beckman, R.J., Conover, W.J.: A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 42 (1), 55-61 (2000), https://dx.doi. org/10.1080/00401706.2000.10485979                                                           |
| [24]   | Mersmann, O., Bischl, B., Trautmann, H., Preuss, M., Weihs, C., Rudolph, G.: Exploratory landscape analysis. In: Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation. pp. 829-836 (2011), https://dx.doi.org/10.1145/2001576.2001690                                                |
| [25]   | Muñoz, M.A., Kirley, M., Halgamuge, S.K.: Exploratory landscape analysis of continuous space optimization problems using information content. IEEE Transactions on Evolutionary Computation 19 (1), 74-87 (2015), https://dx.doi.org/10.1109/TEVC.2014.2302006                                                    |
| [26]   | Ochoa, G., Tomassini, M., Vérel, S., Darabos, C.: AStudy of NKLandscapes' Basins and Local Optima Networks. In: Proceedings of Genetic and Evolutionary Computation Conference. pp. 555 - 562 (July 2008)                                                                                                         |
| [27]   | Pitzer, E., Affenzeller, M., Beham, A., Wagner, S.: Comprehensive and automatic fitness landscape analysis using HeuristicLab. In: Computer Aided Systems Theory - EUROCAST 2011, pp. 424-431. Springer Berlin Heidelberg (2012)                                                                                  |
| [28]   | Renau, Q., Doerr, C., Dreo, J., Doerr, B.: Exploratory landscape analysis is strongly sensitive to the sampling strategy. In: International Conference on Parallel Problem Solving from Nature. pp. 139-153. Springer (2020)                                                                                      |
| [29]   | Renau, Q., Dreo, J., Doerr, C., Doerr, B.: Expressiveness and robustness of landscape features. In: Proceedings of the Genetic and Evolutionary Computation Conference Companion. p. 2048-2051. GECCO '19, Association for Computing Machinery, New York, NY, USA (2019), https://doi.org/10.1145/3319619.3326913 |
| [30]   | Rivest, R.L.: Partial-match retrieval algorithms. SIAM Journal on Computing 5 (1), 19-50 (1976), https: //dx.doi.org/10.1137/0205003                                                                                                                                                                              |
| [31]   | Sagan, H.: Space-filling curves. Springer Science & Business Media (2012), https://doi.org/10.1007/ 978-1-4612-0871-6                                                                                                                                                                                             |
| [32]   | Sallam, K.M., Elsayed, S.M., Sarker, R.A., Essam, D.L.: Landscape-assisted multi-operator differential evolution for solving constrained optimization problems. Expert Systems with Applications 162 , 113033 (Dec 2020). doi:10.1016/j.eswa.2019.113033                                                          |

- [33] Skilling, J.: Programming the Hilbert curve. In: Bayesian Inference and Maximum Entropy Methods in Science and Engineering. American Institute of Physics Conference Series, vol. 707, pp. 381-387 (Apr 2004), https: //dx.doi.org/10.1063/1.1751381
- [34] van Aardt, W.A., Bosman, A.S., Malan, K.M.: Characterising neutrality in neural network error landscapes. In: 2017 IEEE Congress on Evolutionary Computation (CEC). pp. 1374-1381 (2017), https://doi.org/10. 1109/CEC.2017.7969464
- [35] Vassilev, V.K., Fogarty, T.C., Miller, J.F.: Information characteristics and the structure of landscapes. Evolutionary Computation 8 (1), 31-60 (2000), https://dx.doi.org/10.1162/106365600568095
- [36] Vassilev, V.K., Fogarty, T.C., Miller, J.F.: Smoothness, ruggedness and neutrality of fitness landscapes: from theory to application. In: Advances in evolutionary computing, pp. 3-44. Springer (2003)
- [37] Verel, S., Liefooghe, A., Jourdan, L., Dhaenens, C.: On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives. European Journal of Operational Research 227 (2), 331-342 (Jun 2013), https://doi.org/10.1016/j.ejor.2012.12.019
- [38] Weinberger, E.: Correlated and Uncorrelated Fitness Landscapes and How to Tell the Difference. Biological Cybernetics 63 (5), 325-336 (September 1990)
Johannes J. Pienaar Anna S. Bosman Katherine M. Malan
Hilbert curves for efficient exploratory landscape analysis neighbourhood sampling
Landscape analysis aims to characterise optimisation problems based on their
objective (or fitness) function landscape properties. The problem search space
is typically sampled, and various landscape features are estimated based on the
samples. One particularly salient set of features is information content, which
requires the samples to be sequences of neighbouring solutions, such that the
local relationships between consecutive sample points are preserved. Generating
such spatially correlated samples that also provide good search space coverage
is challenging. It is therefore common to first obtain an unordered sample with
good search space coverage, and then apply an ordering algorithm such as the
nearest neighbour to minimise the distance between consecutive points in the
sample. However, the nearest neighbour algorithm becomes computationally
prohibitive in higher dimensions, thus there is a need for more efficient
alternatives. In this study, Hilbert space-filling curves are proposed as a
method to efficiently obtain high-quality ordered samples. Hilbert curves are a
special case of fractal curves, and guarantee uniform coverage of a bounded
search space while providing a spatially correlated sample. We study the
effectiveness of Hilbert curves as samplers, and discover that they are capable
of extracting salient features at a fraction of the computational cost compared
to Latin hypercube sampling with post-factum ordering. Further, we investigate
the use of Hilbert curves as an ordering strategy, and find that they order the
sample significantly faster than the nearest neighbour ordering, without
sacrificing the saliency of the extracted features.