Sarah L Thomson
email: s.l.thomson@stir.ac.uk

Gabriela Ochoa
email: gabriela.ochoa@stir.ac.uk

Sébastien Verel
email: verel@uni-littoral.fr

Clarifying the Difference in Local Optima Network Sampling Algorithms

Keywords: Combinatorial Fitness Landscapes, Local Optima Networks, Quadratic Assignment Problem

We conduct the first ever statistical comparison between two Local Optima Network (LON) sampling algorithms. These methodologies attempt to capture the connectivity in the local optima space of a fitness landscape. One sampling algorithm is based on a random-walk snowballing procedure, while the other is centred around multiple traced runs of an Iterated Local Search. Both of these are proposed for the Quadratic Assignment Problem (QAP), making this the focus of our study. It is important to note the sampling algorithm frameworks could easily be modified for other domains. In our study descriptive statistics for the obtained search space samples are contrasted and commented on. The LON features are also used in linear mixed models and random forest regression for predicting heuristic optimisation performance of two prominent heuristics for the QAP on the underlying combinatorial problems. The model results are then used to make deductions about the sampling algorithms' utility. We also propose a specific set of LON metrics for use in future predictive models alongside previously-proposed network metrics, demonstrating the payoff in doing so.

Introduction

Local optima networks (LON ) are a partial fitness landscape of a combinatorial optimisation problem [START_REF] Ochoa | A Study of NK Landscapes' Basins and Local Optima Networks[END_REF]. Features of the networks have repeatedly been linked to heuristic search in various famous problem domains ( [START_REF] Ochoa | Mapping the Global Structure of TSP Fitness Landscapes[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF][START_REF] Mcmenemy | How perturbation strength shapes the global structure of tsp fitness landscapes[END_REF][START_REF] Iclanzan | Data-driven local optima network characterization of qaplib instances[END_REF]). LONs were originally used for fully-enumerated fitness landscapes and therefore reflected small problems. More recently, sampling algorithms have been proposed [START_REF] Iclanzan | Data-driven local optima network characterization of qaplib instances[END_REF][START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF][START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF] with the aim of using the analysis on larger problems. The nature, biases and resultant samples of the algorithms are of critical importance in moving LON analysis towards real-world systems. In this study we conduct the first comparison of sampling algorithms for this purpose. Features of the samples are contrasted, and the features are used in linear and random forest models for predicting heuristic performance on the underlying combinatorial problems. In other words, we examine which of the two most recent sampling algorithms ( [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF][START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF]) produces LONs with more predictive power, in terms of their ability to measure search difficulty. We also propose using funnel metrics as features for predictive models of algorithm performance. A funnel is a large fitness landscape feature and is essentially a basin of attraction at the level of local optima. The main contributions of this paper are:

1. The first descriptive comparison of local optima network sampling algorithms. 2. Contrast of the differences in predictive power of the LON features for explaining meta-heuristic variance. 3. The addition of a specific funnel metric set for predicting algorithm performance in LON models using linear and random forest regression.

Definitions

A fitness landscape [START_REF] Stadler | Fitness landscapes[END_REF] is a triplet (S, N, f ) where S is the set of all possible solutions, N : S -→ 2 S , a neighbourhood structure, is a function that assigns to every s ∈ S a set of neighbours N (s), and f is a fitness (objective value) function such that f : S -→ R, where the fitness value is a real number that can be viewed as the height of a given solution in the landscape.

A local optima network is a representation of the fitness landscape at the level of local optima. We now formally define the constituent parts of a local optima network, before proceeding to describe the object as a whole.

Nodes. The set of nodes, LO, is comprised of local optima, i.e. a solution lo i satisfies the condition that it has superior fitness to all other solutions in its neighbourhood: ∀n ∈ N (lo i ) : f (lo i ) ≥ f (n), where N (lo i ) is the neighbourhood and n is a single neighbour.

Edges. The set of edges, E, consists of directed and weighted links. An edge is traced if the probability of 'escape' -using perturbation and then hill-climbing -from the source node to the destination is greater than zero, and is weighted with the probability. Formally, local optima lo i and lo j form the source and destination of an edge iff w ij > 0.

Local optima network (LON). The weighted local optima network LON = (LO, E) is a graph where the nodes lo i ∈ LO are the local optima, and there exists an edge e ij ∈ E, with weight w ij , between two nodes lo i and lo j if w ij > 0. Note that w ij may be different than w ji . Thus, two weights are needed in general, and so a local optima network is an oriented transition graph.

The 'Network' Feature Set

We now describe and introduce terminologies for the complex network features we use in this study. The number of edges found in one of our sampled networks is referred to as edges. We also use the mean fitness found in the sample of local optima comprising the network nodes, and denote this as meanfitness. The network diameter (the longest possible length in number of edges separating two nodes), diam, and the average out-degree (the number of edges which are directed away from a particular node), outdegree, are also included.

The 'Funnel' Feature Set

A fitness landscape funnel is a basin of attraction at the level of local optima. Instead of a single local optimum as an attractor for surrounding solutions, every member of a funnel is a local optimum. Each of these is on a path which, if followed, terminates at a single high-quality local optimum -the bottom of the funnel. These funnel bottoms are also referred to as sinks from a graph theory perspective, and are simply the nodes with no outgoing edges. We uncover the funnels by first pruning off any LON edges which are non-improving in fitness, and then by commencing a depth-first search from a funnel bottom, in the process revealing the sequences of local optima which must terminate there.

We use four features relating to fitness landscape funnels in our analysis here. The first is simply the amount of funnels as a proportion of the total node count, funnel. We also use the number of compressed local optima. The compression process happens during the network pre-processing for funnel analysis, where the non-improving edges are removed, and then plateaus at the local optima level are joined into a single node. The metric ncoptima is the number of optima left after this processing of the network. Also considered is the incoming flow to suboptimal funnel bottoms or sinks (which can be calculated by the addition of the incoming edge-weights to the sink nodes) -this is referred to as substrength, and it comes as a proportion of the total flow going to all funnel bottoms, including optimal ones. The fourth funnel feature is sinkfitness, which is simply the mean fitness of funnel bottom nodes (sinks in the network).

Experimental Setting

Benchmark Test Problem

Our test instances are Quadratic Assignment Problems (QAP), where an optimal assignment of N facilities to N locations is sought. Each facility, location pair has an associated distance and flow between them. An instance of QAP, then, is defined with two matrices: a distances matrix, and a flow matrix. A solution to the QAP, π, is a permutation of facilities, while the objective value is calculated as follows:

C(π) = n i=1 n j=1 a ij b πiπj
where A = {a ij } and B = {b ij } are the distance and flow matrices, respectively.

We focus exclusively on the much-studied Quadratic Assignment Problem Library (QAPLIB). QAPLIB boasts a diverse selection of problems, both synthetic and real-world, from different classes. Here we use 30 moderately-sized instances, from various problem classes. The problem dimensions N are between 25 and 50, meaning 25-50 facilities to be assigned. Occurrences of the QAP can largely be placed in one of four categories: Uniform random distances and flows; random flows on grids; real-world problems; and random real-world like problems. Instances from all of these categories are considered.

Uniform random distances and flows. This class of QAP is known to be difficult for optimisation heuristics. The entries for both the distance and flow matrices are taken at random from a Gaussian distribution. The naming convention for them in the QAPLIB is tainna, with nn being the problem dimension. The nature of the points on the plane is random. We use here tai25a, tai30a, tai35a, and tai40a.

Random flows on grids. For these problems, the locations are each situated in one square on an m*n grid, which is rectangular. The flow matrix entries are generated randomly. From this category, we use nug25, nug27, nug28, nug30, sko42, tho30, tho40, and wil50.

Real-world problems. These problems arise from practical applications. We briefly describe the instances used in this article. We use two bur instances, which comprise of the average stenotypist's typing data. Also included are some kra problems: these were used to plan a German hospital. The esc instance we use presents itself in computer science (the testing of sequential circuits). In the ste problem set, the goal is to minimise the length of connections between units that have to be placed on a rectangular grid. The real-world instances we use in our analysis are bur26a, bur26b, esc32e, kra30a, kra30b, kra32, lipa30a, lipa30b, lipa40a, lipa40b, ste36a, ste36b, and ste36c.

Random real-world like problems. The naming convention is tainnb, where nn is the problem dimension. We use the problems tai25b, tai30b, tai35b, and tai40b. Miscellaneous. The chr set of problems are a special occurrence of the QAP and do not fit into any of the previous categories. The flow matrix forms a mathematical tree, with no particular specification for the definition of the distance matrix. We use chr25a from this set.

The Sampling Algorithms

Snowball (Chain-referral) Sampling. Snowball sampling is a process, originating in the social sciences, where each respondent in a survey asks a few of their friends to complete the survey as well. The social network of survey respondents is then analogous to an ever-growing rolling snowball, growing larger with every step. This was introduced as the mechanism behind a sampling algorithm for local optima networks very recently [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF].

To construct the local optima network sample, the algorithm first starts from a random solution and hill-climbs to a local optimum, x t . This local optimum is the start-point of the walk. Then the recursive snowball procedure branches from x t : the neighbouring local optima are explored, and then the neighbours of those are explored, up to a specified 'depth' given by a parameter d. The local optima are found by carrying out mutations followed by hill-climbing. The hill-climbing used to obtain local optima is best-improvement and uses a random pairwise swap of facilities as the operator. The perturbation operator is four of these mutations in a row, i.e. four swaps.

After the snowball expansion of x t , the algorithm returns to the initial local optimum x t and continues on its walk from there. Specifically, a neighbouring local optimum which is not already present in the walk is found. If all neighbouring local optima have been expanded, a new random solution is generated and hill-climbing applied to find a local optimum for inclusion in the random walk. Then, that optimum is included in the walk as x t+1 and becomes the 'centre' from which a new iteration of snowballing expansion begins. The algorithm then goes back x t+1 and finds a neighbour to be x t+2 on the walk, and so on. The walk x t...t+n is of a length given by the sampling parameter l.

LON snowball sampling is configurable with three parameters: l (length of random walk), d (depth of snowballing), and m (number of sampled edges). We use the three (l, d) combinations suggested in [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF] for obtaining our LONs: (100, 60), (400, 30), and (100, 30), alongside the suggested value for m, which is two. The parameters l and d can be used to tune the sample obtained.

A full description of LON snowballing, including pseudo-code and descriptive figures, can be found in [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF]. Algorithm 2 is particularly helpful. In the text that follows, we refer to this sampling algorithm using the term walkSample, to reflect the fact it is a sample based on a random walk.

Iterated Local Search Sampling. This sampling method is built around the competent Improved Iterated Local Search [START_REF] Stützle | Iterated local search for the quadratic assignment problem[END_REF] heuristic for the QAP, and was proposed in [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF]. The local optima network is built while the heuristic algorithm is attempting optimisation. In essence, 200 ILS runs are commenced from different random solutions; for each run, as the ILS progresses in the minimisation process, all local optima encountered are recorded. Alongside these, every transition between two local optima (using perturbation and then hill-climbing) is saved as an edge in the network (LON), or if the edge has already been recorded, the weight of the edge is incremented. Each run terminates either when the global optimum is found, or after 10000 iterations without improvement.

The local search stage uses a first improvement hill-climbing variant with a random pairwise swap operator (the same as the LS in walkSample. The perturbation operator exchanges k randomly chosen items. In the original ILS algorithm code, a few settings for k are suggested. We use a subset of those in this study, namely: n 4 , n 2 , and 3n 4 . These settings are important for the construction of the LONs, because the perturbation mechanism controls the discovery and inclusion of network edges (connections between nodes). A comprehensive description of both the original ILS and the modified version for LON sampling can be found in [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF]. We refer to this sampling process henceforth by optSample, to reflect the fact that the sampling is done alongside heurisitic optimisation.

Figure 1 compares a sample from optSample algorithm with one obtained using walkSample. The Figure serves as a precursor to a formal comparative analysis later on, but captures a remarkable amount of information. Both are sampled local optima networks, extracted from the same Quadratic Assignment Problem. The sample in Figure 1a is extracted using optSample, while walkSample was used to obtain the network in Figure 1b.

In red is the global optimum (or the global optima), while all other nodes are grey in colour. Both samples have been capped by fitness: only optima in the top 5% (Figure 1a) and 0.05% (Figure 1b) of the fitness distributions are included in the plots. The latter percentile is lower because the sample had orders of magnitude more optima and edges crowding its network at this fitness level.

Comparing the two plots, there are striking differences. Taking them in turn, we can see that the high-quality sub-LON which optSample produces (Figure 1a) is extremely sparse, and indeed the two global optima are disconnected (isolated) nodes in this graph. This could indicate that a heuristic would be unlikely to end up at one of these apparently inaccessible nodes; alternatively, it could be that the global optima are well connected nodes (or even network hubs), but their network neighbours exist in lower-quality fitness levels.

Casting our attention to the LON in Figure 1b --produced by walkSample --it is a markedly more dense network. There are significantly more nodes and edges between them. Recalling that this is only the highest-quality 0.05% of local optima which were sampled, the implication is that this sampling method exposes much of the neighbourhood (at the local optima level) surrounding the global optimum, and may be more suited to characterising this promising area of the search space. A global optimum can be seen in the centre Figure 1b.

Heuristics

To gain a view of empirical complexity of the chosen problem instances, we must collect optimisation data by competent heuristic algorithms.

Two competitive algorithms from the literature are deployed: Stützle's Improved ILS for the QAP (IILS) [START_REF] Stützle | Iterated local search for the quadratic assignment problem[END_REF], and Taillard's Robust Taboo Search for the QAP (TS) [START_REF] Taillard | Robust taboo search for the quadratic assignment problem[END_REF]. IILS comes with a wealth of viable parameter configurations. We choose here first improvement hill-climbing, in combination with a perturbation strength of 3n 4 pairwise exchanges of facilities. To quantify their performance on the thirty problems, we use the performance gap metric, which is the obtained fitness as a proportion of the desired fitness after a fixed budget of iterations (1000 here). We take the average of this value over 1000 runs. In the results that follow, this metric is referred to as IILSp and TSp (shorthand for ILS and TS performance).

Predictive Model Setup

We use both linear mixed models and random forest regression models for algorithm performance prediction. When considering our full (90 sampled LON per sampling algorithm) sets of observations, we conduct random repeated subsampling cross-validation (also known as bootstrapping) for 1000 iterations, each time shuffling the observations randomly. We do this with an 80-20 training-test split. For the smaller datasets (split by parameter set), we only have 30 observations, so we do not split into training and test data, but use all 30 for regression.

In the results that follow, we focus on the R 2 associated with each of the models, which quantifies the amount of variance in the response variable (in our case either IILSp or TSp) explainable using the set of features in question (either the full metric set, which is the union of the network and funnel features, or the funnel or network features exclusively).

For the random forest regression models, we also report the best four predictors for the models in terms of importance, computed during the model formation. We do not conduct random forest regression on the smaller sets of 30 LONs -only on our two full sets of 90 LONs -to provide enough data for the process.

Sampling Algorithm Comparison

To have two sampling algorithms produce networks which are somehow correlated would be encouraging. We rely on sampling algorithms for local optima networks on problems of any reasonable size. As a result, before the algorithms reach real-world applications, they must be refined. Table 1 provides Spearman correlation coefficients, calculated between features of LONs obtained by opt-Sample and walkSample, respectively. An indication of the p-value is given, as described in the caption. Both network features and funnel features are included as variables. These were introduced in Sections 2.1 and 2.2, respectively.

Most of the features have a fairly weak correlation, but with a p-value less than the threshold 0.05, giving evidence against the null hypothesis that the sample features are unrelated. The correlation for the out-degree of the samples is moderate, and has an encouragingly small p-value of < 0.001. The two strongest associations are for the mean fitness in the network, and the mean fitness of the funnel bottoms (sinkfitness), respectively. These two features seem to show a great deal of agreement, gaining very strong correlations indicating statistical significance.

Table 2 attempts to quantify the predictability of optSample and walkSample. This is done by looking at the ranges of values for important features of the obtained LONs. The minimum value in the sample set is represented as a proportion of the maximum value found.

Looking first at the column giving the information on optSample, we can see that the proportions are mostly much nearer to zero than to one. This means that the smallest values are completely dwarfed by the values of the largest. The deduction from this is that you cannot tell this sampling methodology exactly what which trajectory to follow. This is advantageous in the sense that it is not rigid or artificially contained, but could also have the potential to result in unreliable samples. It is, however, promising when considering the wide range of the average outdegree of the obtained networks; this hints at hub-and-spoke network structure being exposed (where it exists).

In contrast, the column showing the range information for walkSample shows much larger proportions between 0.29 and 0.99, telling us that this algorithm produces a fairly predictable and therefore tuneable number of nodes and edges. However, because the range for the outdegree is very small (the smallest outdegree is almost as big as the largest, arising due to the nature of the algorithm), important connections between local optima will likely be missed and so-called network 'hubs' may never be found. These are critical because much of the flow of the network -in our case, prospective heuristic search trajectories -may pass through them or be drawn towards them.

Prediction of Heuristic Competence on the Problems

For our predictions, we computed models based on all 90 sampled LONs (that is, 90 per sampling algorithm), which comprise LONs generated using three parameter sets, each of them used for all 30 instances. We also wanted to see the predictive power when using LONs generated using a single, fixed set of parameters; therefore, we also considered subsets of the LON datasets, where we split into sets of 30.

For optSample, the parameter difference came from the perturbation setting used in the construction of the edges, connecting optima in the sample. The three settings of this were n 4 , n 2 and 3n 4 , and the LON sets corresponding to these are referred to as set 1, set 2, and set 3 (in that order) in the models that follow. The set of all 90 LONs together is set 4.

In the case of walkSample, the parameter combinations were of the sampling depth, d and length of walk, l. The sets used were (60, 100), (30, 400), and (30, 100); these are referred to as set 1, set 2, and set 3, respectively. Again, the set of all 90 LONs combined is set 4.

Table 3 reports the R 2 for several predictive models. Each row gives information on the models in question. The columns represent, in order, the sampling algorithm used to obtain the LON samples (Sample); the type of predictive methodology (Model, one of random forest, rf, or linear mixed model, lm); the set of features used as predictors (Features); the parameter set used in the construction of the network samples (Param); and an indication of whether ILS or TS performance gap was the response variable. The higher the values for the models, the more explainable the variance in the optimisation metric is using that set of predictors. Focusing on the optSample entries, we can see that the general trend is that more variance is explained in the TS than in the ILS (compare the two columns). When considering the 90 LONs together in a linear setting, using all the features together (all ) produces the strongest models of the three predictor sets (all, funnel, and network ). For TS, using only the network set is nearly as effective as using all the predictors. In the case of ILS, funnel predictors can explain more search variance than using exclusively network features.

The poorest fits of the optSample linear mixed models come from the prediction of ILS response using network variables. If we compare the random forest models with the linear mixed models -taking the bottom three optSample rows against the set 4 (full set) optSample rows -it can be seen that for ILS, these are stronger than the linear models. Notice that in the optSample random forest setting with TS response that the models are strong. We also notice that using exclusively the funnel set produces a stronger R 2 value than using the complete set or the network one.

The walkSample models are summarised in the lower half of Table 3. Let us first consider the values obtained using the full set of 90 sampled LONs (the rows where column Param. is stated as set 4 ). Comparing these three with one another, we see that for ILS, the best predictor set is all the metrics togetheralthough none of these models fit well. For Tabu Search response, using the full feature set produces a high R 2 value of 0.968. Using exclusively the funnel set produces a model which is not much weaker. Using only the network features has a very poor R 2 .

Just like with the optSample models, the trend over the models here is that more of the variance in TS is explainable with the predictors than that of the ILS. This can be seen by locating the ILS column and checking against the neighbouring TS column entries.

The random forest models for walkSample LONs are useless in the case of ILS response. In the case of TS, however, we can see that using the full set of features produces a strong model with the highest of the three R 2 values. The funnel set performs next-best, with the network features ranking lowest within these three settings.

We now compare the random forest models for the two sampling methodologies and look at six rows in total (all rows where Model is rf ). What is clear, if going by these alone, is that optSample has more predictive power for these problems and algorithms than walkSample. This can be seen in the generally higher R 2 entries.

Doing the same comparison for the linear models for each of our two sampling algorithms (every row where Model is stated as lm), it is less of a clear-cut distinction. Neither produce strong linear models with ILS as the response variable. For TS, the funnel settings are much weaker with respect to optSample than walkSample. In contrast, the network set produce much stronger models when extracted from optSample LONs when compared with walkSample.

Table 4 reports rankings of features in terms of importance in their predictive random forest models. The column labelled '1' shows the most important feature from the model specified in that row. In the cases where the full feature set are used as predictors (the rows where Features are stated as 'all'), there were eight variables used, but only the best four are reported in the Table .  We can see from examining the '1' row that almost always, the meanfitness (mean local optimum fitness in the network) and sinkfitness (mean funnel bottom fitness in the network) are ranked as the most important. outdegree and substrength also feature more than once among the top two.

Correlation Study. Figure 2 is a correlation matrix contrasting observed features (edges, meanfitness, outdegree, diam, substrength, sinkfitness, funnel, and ncoptima), of the local optima networks sampled by optSample with heuristic performance on the underlying combinatorial problem (ILSp and TSp).

Each entry in the upper-right triangle is a Spearman correlation between the row and column variables. The diagonal shows density plots for the features, with the lower-left triangle showing pairwise scatterplots. The colour in these represent the split into sampling parameter group.

We are, of course, principally interested in the connections between our heuristic performance metrics, and the features of the sampled LONs. The easiest way to approach this is by locating the ILSp and TSp rows and looking along them, checking against the intersections for interesting correlations. For ILSp, we can see that (if we are taking p <0.05 to be indication of statistical significance), there is a moderate positive correlation with indication of statistical significance with substrength. Notice also by checking the ILSp column intersections that there are similar moderate positive correlations with p < 0.001 with ncoptima and with funnel. In the case of TSp, there are weak negative correlations with substrength and ncoptima.

Figure 3 has precisely the same layout as Figure 2, but shows results for walkSample instead of optSample.

If we assess things in the same way: the intersections of ILSp and TSp with other columns and rows. In general, there are much weaker correlations with the performance gap data here than we saw for the features of the optSample LONs and there are no correlations between ILSp and TSp and landscape features which have p < 0.05. This can be seen in the smaller values and lack of p-value indications for the feature-performance pairs.

Conclusions and Thoughts

We have conducted an empirical comparison between two sampling algorithms for local optima networks of the Quadratic Assignment Problem. It is essential to refine sampling for local optima networks before we bring them to real world systems. A descriptive statistical comparison was reported, as well as findings pertaining to the algorithms' predictive power in estimating the performance gap obtained by two prominent meta-heuristics on the underlying combinatorial problem.

We found that the two sampling methods exhibited some agreement in the networks they produced and that we could reject the null hypothesis that they produce completely independent samples. They differed from a descriptive perspective in that walkSample was tuneable and predictable, while optSample varied widely but seemed good at finding hub-and-spoke structure in the local optima space.

A correlation study for the LON features and the heuristic performance metrics was conducted. The correlations were clearer when considering the features of the LONs obtained using optSample than walkSample.

We also worked on explaining heuristic algorithm performance on the problems using linear and random forest models, and found that the sampled LON features (for both optSample and walkSample) better fit the TS response variable than the ILS one. We saw that overall, including both the funnel metric set and the network set would be advantageous in explaining search discrepancies for these two heuristics. An exceptions to this was found when including optSample LON features in random forest regression with TS response -using only the funnel set produced a slightly higher R 2 values in this case. For both optSample and walkSample, the extracted funnel metrics proved useful in models involving TS as a response variable. Going off the random forest models alone, optSample generally had more explanatory power than its competitor, for these choices of instances and heuristics.

From the random forest rankings, the most important predictors were those pertaining to fitness in the sampled networks: the fitness of funnel bottoms, and of nodes in general in the network. This hints that perhaps fitness levels in the local optima space are more pertinent to heuristic search than the subset of transition edges sampled by the LON algorithms.

  (a) optSample LON sample for instance wil50. Top 5% of optima by fitness. (b) walkSample LON sample for instance wil50. Top 0.5% of optima by fitness.

Fig. 1 .

 1 Fig. 1. A representative LON for an instance from the Quadratic Assignment Problem Library. Figures show LON samples as obtained by two different sampling methods. The global optimum or optima in the samples are shown in red.

Fig. 3 .

 3 Fig. 3. Correlation matrices of performance metrics and walkSample-produced LON features. Lower triangle: pairwise scatter plots. Diagonal: density plots. Upper triangle: pairwise Spearman's rank correlation, * * * p < 0.001, * * p < 0.01, * p < 0.05.

Table 1 .

 1 Spearman's rank correlation, * * * p < 0.001, * * p < 0.01, * p < 0.05, for pairs of features from the two sampling algorithms.

	Predictor	Corr.
	optima	0.239 *
	edges	0.149
	ncoptima 0.249 *
	outdegree 0.441 * * *
	diam	-0.212 *
	meanfitness 0.990 * * *
	sinks	0.112
	sinkfitness 0.990 * * *
	substrength 0.274 * *

Table 2 .

 2 The variance (represented as the minimum value as a proportion of the maximum value) of some key metrics for the sampled LONs.

	Predictor optSample walkSample
	optima	0.0003	0.2499
	edges	0.0003	0.2499
	outdegree 0.0508	0.9999
	sinks	0.0003	0.2499

Table 3 .

 3 R 2 values for linear and random forest models to explain heuristic performance variation on the combinatorial problems.

	Sample Model Features Param. ILS TS
	optSample lm	all	set 1 0.124 0.963
	optSample lm	all	set 2 0.135 0.963
	optSample lm	all	set 3 0.182 0.964
	optSample lm	all	set 4 0.109 0.974
	optSample lm	funnel set 1 0.114 0.125
	optSample lm	funnel set 2 0.117 0.124
	optSample lm	funnel set 3 0.162 0.127
	optSample lm	funnel set 4 0.094 0.125
	optSample lm network set 1 0.069 0.957
	optSample lm network set 2 0.065 0.957
	optSample lm network set 3 0.075 0.958
	optSample lm network set 4 0.061 0.954
	optSample	rf	all	set 4 0.307 0.820
	optSample	rf	funnel set 4 0.216 0.892
	optSample	rf network set 4 0.247 0.775
	walkSample lm	all	set 1 0.640 0.980
	walkSample lm	all	set 2 0.339 0.955
	walkSample lm	all	set 3 0.294 0.964
	walkSample lm	all	set 4 0.158 0.968
	walkSample lm	funnel set 1 0.200 0.956
	walkSample lm	funnel set 2 0.104 0.946
	walkSample lm	funnel set 3 0.086 0.954
	walkSample lm	funnel set 4 0.061 0.901
	walkSample lm network set 1 0.062 0.027
	walkSample lm network set 2 0.245 0.042
	walkSample lm network set 3 0.248 0.117
	walkSample lm network set 4 0.112 0.016
	walkSample rf	all	set 4 0.045 0.910
	walkSample rf	funnel set 4 -0.017 0.632
	walkSample rf network set 4 -0.084 0.249

Table 4 .

 4 Predictor rankings for the size random forest models.

	Sample Features Param Resp.	1	2	3	4
	optSample funnel set 4 ILS sinkfitness substrength	funnel	ncoptima
	optSample network set 4 ILS meanfitness	edges	outdegree diameter
	optSample	all	set 4 ILS meanfitness sinkfitness ncoptima	funnel
	walkSample funnel set 4 ILS sinkfitness ncoptima	funnel	substrength
	walkSample network set 4 ILS meanfitness outdegree	edges	diameter
	walkSample	all	set 4 ILS outdegree meanfitness sinkfitness ncoptima
	optSample funnel set 4 TS sinkfitness ncoptima	funnel	substrength
	optSample network set 4 TS meanfitness	edges	diameter outdegree
	optSample	all	set 4 TS sinkfitness meanfitness substrength	funnel
	walkSample funnel set 4 TS substrength	funnel	sinkfitness ncoptima
	walkSample network set 4 TS meanfitness outdegree diameter	edges
	walkSample	all	set 4 TS sinkfitness meanfitness outdegree ncoptima

  Correlation matrices of performance metrics and optSample-produced LON features. Lower triangle: pairwise scatter plots. Diagonal: density plots. Upper triangle: pairwise Spearman's rank correlation, * * * p < * * p < 0.01, * p < 0.05.

	TSp edges meanfitness outdegree diam substrength sinkfitness TSp edges meanfitness outdegree diam substrength sinkfitness	r= r= r= r= r= r= r= ncoptima -0.289** 0.470*** 0.115 -0.034 0.346*** 0.652*** 0.116 r= 0.100 r= 0.898*** r= 0.315** r= 0.395*** ncoptima r= -0.146 r= 0.286** r= 0.315**	r= r= r= r= r= r= r= funnel -0.169 -0.204 0.025 -0.365*** -0.220* 0.867*** 0.035 r= 0.094 r= 0.640*** r= 0.502*** r= 0.334** funnel r= -0.228* r= 0.295** r= 0.502***		r= r= r= r= r= r= r= ILSp -0.082 0.183 -0.117 -0.133 0.104 0.682*** -0.116 r= -0.082 r= 0.039 r= -0.094 r= -0.138 ILSp r= -0.012 r= 0.161 r= -0.094		r= r= r= r= r= r= TSp -0.253* 0.103 -0.154 -0.028 -0.218* 0.102 r= 0.031 r= 0.171 r= -0.150 TSp r= 0.071 r= 0.016 r= 0.171			r= r= r= r= r= edges -0.072 0.406*** 0.832*** 0.158 -0.070 r= -0.025 r= 0.482*** edges r= -0.025 r= 0.217* r= -0.025		r= r= r= r= meanfitness outdegree 0.116 -0.102 r= 0.355*** -0.012 r= -0.227* 0.999*** r= 0.121 r= -0.027 meanfitness outdegree r= -0.247* r= 0.369*** r= 0.249* r= 0.054 r= 1.000*** r= -0.027		r= r= diam 0.089 -0.099 diam r= -0.058 r= -0.247*		r= substrength sinkfitness -0.002 substrength sinkfitness r= 0.249*	2 3 0500000 1000000 1500000 2000000 0e+00 2e+08 4e+08 6e+08 5 10 15 2 3 40000 80000 120000 160000 0e+005e+081e+09 0.49999 0.50000 0.50001 0.50002 6 0 2500 5000 7500 10000 0.000.250.500.751.00 0e+00 2e+08 4e+08 6e+08 8 10 12 0.9992 0.9996 1.0000 0e+005e+081e+09
																																1 1
	r= r= 4e-04 Fig. 2. r= 0.440*** 0.663*** ncoptima funnel ILSp 2e-04 6e-04 0.806*** r= 0.036 ncoptima funnel ILSp 0.001 0.002 0.003 0.004	0e+00 0.000	r= 0.050 0.075 r= 0.93 0.550*** 0.025 0.010 0.90 0.96	0.000 0.87	1.04 1.04	1.02 1.02	1.00 1.00	3 3	2 2	1 1	2000000	1500000 160000	1000000 120000	500000 80000	0 40000	4e+08 1e+09 6e+08	2e+08 5e+08	0e+00 0e+00	5 0.49999 0.50000 10 0.50001 15 0.50002	10000 12	2500 8 5000 10 7500	0 6	1.00 1.0000	0.50 0.9996 0.75	0.25 0.9992	0.00	6e+08	4e+08 1e+09	2e+08 5e+08	0e+00 0e+00	0 2000 4000 6000 8000 0.000 0.025 0.050 0.075 1.00 1.02 1.04 1e+052e+053e+05 0.870.900.930.96 1.00 1.02 1.04

Acknowledgements. This work is supported by the UK's Engineering and Physical Sciences Research Council (grant number EP/J017515/1). Data generated during this research are available from the Stirling Online Repository for Research Data (http://hdl.handle.net/11667/128).


Sarah L Thomson
email: s.l.thomson@stir.ac.uk

Gabriela Ochoa
email: gabriela.ocha@cs.stir.ac.uk

Sébastien Verel
email: verel@univ-littoral.fr

The Fractal Geometry of Fitness Landscapes at the Local Optima Level

Keywords: Fitness Landscapes, Fractal Analysis, Local Optima Networks 1 Introduction

A Local Optima Network (LON) encodes local optima connectivity in the fitness landscape of a combinatorial optimisation problem. Recently, LONs have been studied for their fractal dimension. Fractal dimension is a complexity index where a non-integer dimension can be assigned to a pattern. This paper investigates the fractal nature of LONs and how that nature relates to metaheuristic performance on the underlying problem. We use visual analysis, correlation analysis, and machine learning techniques to demonstrate that relationships exist and that fractal features of LONs can contribute to explaining and predicting algorithm performance. The results show that the extent of multifractality and high fractal dimensions in the LON can contribute in this way when placed in regression models with other predictors. Features are also individually correlated with iterated local search performance, and visual analysis of LONs shows insight into this relationship.

one views it using a different resolution or scale. Fractal dimension analysis has been used in diagnostic imaging (detecting colon cancer [START_REF] Esgiar | Fractal analysis in the detection of colonic cancer images[END_REF]; characterising images in mammography [START_REF] Caldwell | Characterisation of mammographic parenchymal pattern by fractal dimension[END_REF]; characterising leukaemia cells [START_REF] Mashiah | Lymphoma and leukemia cells possess fractal dimensions that correlate with their biological features[END_REF]), search and rescue (analysing the layout of victim location after building collapses [START_REF] Saeedi | An algorithmic approach to generate after-disaster test fields for search and rescue agents[END_REF]), and in engineering for the design of antenna [START_REF] Werner | An overview of fractal antenna engineering research[END_REF], among innumerable others. Fractal geometry can also facilitate vast amounts of information being embedded in a comparatively small space. Indeed, branching structures inside human lungs fill space in a fractal way; because of this, the equivalent surface area of a tennis court is compacted within the volume of the lungs. The fractal dimension of human lungs has been measured at approximately 2.88 [START_REF] Uahabi | New approach to the calculation of fractal dimension of the lungs[END_REF], which indicates high spatial complexity and convolution.

Fitness landscapes of some combinatorial optimisation problems have been viewed under a fractal lens [START_REF] Weinberger | Why some fitness landscapes are fractal[END_REF]. Fitness landscapes are both a lucid metaphor and a mathematical object; they contain the set of solutions to an optimisation problem, the fitnesses of those solutions (these can be visualised as the heights), and a function for measuring adjacency between solutions. The study of fitness landscape architecture provides insight about reactions between metaheuristic algorithms and problems. This can serve as a springboard for more informed algorithm design or selection.

The first study to conduct fractal analysis on fitness landscapes [START_REF] Weinberger | Why some fitness landscapes are fractal[END_REF] stipulated that for certain problems, landscape ruggedness scales at different levels of abstraction and that this indicated fractal structure. Subsequent studies have reported similar findings [START_REF] Zelinka | Fractal analysis of fitness landscapes[END_REF][START_REF] Locatelli | On the multilevel structure of global optimization problems[END_REF][START_REF] Richter | Scale-invariance of ruggedness measures in fractal fitness landscapes[END_REF] and some have emphasised the potential lying dormant in the largely untapped field of fractal analysis for landscapes.

A Local Optima Network (LON) [START_REF] Ochoa | A study of nk landscapes' basins and local optima networks[END_REF] models local optima and their connectivity in a fitness landscape. That is, the nodes are local optima and the edges are metaheuristic search transitions between two local optima under a chosen search operation. There is a significant body of evidence suggesting that features of LONs can correlate to, explain, or predict metaheuristic algorithm performance on the underlying combinatorial problem [START_REF] Daolio | Local optima networks of the quadratic assignment problem[END_REF][START_REF] Daolio | Communities of minima in local optima networks of combinatorial spaces[END_REF][START_REF] Verel | Local optima networks with escape edges[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF][START_REF] Ochoa | Mapping the global structure of tsp fitness landscapes[END_REF][START_REF] Ochoa | Understanding phase transitions with local optima networks: Number partitioning as a case study[END_REF].

Little is known about the fractal complexity in LONs and how their fractal nature relates to metaheuristic algorithm performance. Preliminary work has indicated that the fractal dimension might have a connection to search (Thomson et al., 2018a,b). That being said, the latter study considers only small problem instances (size N = 18 for a binary-encoded problem, NK Landscapes). The first study mentioned is on the Quadratic Assignment Problem (QAP) and they consider some benchmark instances from QAPLIB [START_REF] Burkard | Qaplib-a quadratic assignment problem library[END_REF] up to N = 28 (Thomson et al., 2018a), although only two of the library's several instance classes for this problem size range are included; consequently, the fractal analysis is conducted on only 25 QAPLIB instances.

We intend to illuminate understanding of the relationships between fractal geometry in LONs and metaheuristic algorithm performance. The QAP serves as a testbed for the analysis and we use QAPLIB instances, increasing the number of instances considered threefold when compared to previous work (Thomson et al., 2018a) and raising the maximum problem size from 28 to 50. A recent and refined LON construction algorithm [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF] is used to intelligently build LONs for the QAPLIB instances. Features of the LONs, including fractal dimension features, are computed and the parallel between them and performance is investigated using visual tools, correlation analysis, and linear and random forest regression models.

The contributions of this article can summarised as follows:

1. We bring new insight into how multifractal geometry at the local optima level can help explain and predict algorithm performance 2. A significant expansion of the data-set used for fractal analysis in LONs (using more than 3x the previous number of QAPLIB instances and raising N ≤ 28 to N ≤ 50, as well as deploying a recent refined and tested sampling algorithm for constructing the LONs) 3. Enhanced statistical techniques for properly validating the use of LON fractal analysis for algorithm explanation and prediction (random forest to model non-linearities; random repeated subsampling cross-validation; using intelligible predictors such as the extent of multifractality and the median fractal dimension).

The article is structured as follows: Section 2 contains the necessary background information to render this article self-contained; Section 3 details aspects of the methodology used; Section 4 gives the experimental setup, with Section 5 presenting the results; finally, Section 6 finishes the article with conclusions and directions for future work.

Preliminaries

Fitness Landscapes

A fitness landscape [START_REF] Stadler | Fitness landscapes[END_REF] is composed of three parts, (S, N, f ) : S is the full solution set; N : S -→ 2 S is known as the neighbourhood function and assigns a set of adjacent solutions N (s) to every s ∈ S; and f is a fitness function f : S -→ R that provides a mapping from solution to associated fitness. That fitness can be conceptualised as the solution height within the landscape metaphor.

The analysis of fitness landscape objects can provide an intense understanding of optimisation problems and their reactions with metaheuristic algorithms [START_REF] Pitzer | A comprehensive survey on fitness landscape analysis[END_REF]. Indeed, landscapes have been used to facilitate algorithm selection [START_REF] Hoos | Search space features underlying the performance of stochastic local search algorithms for max-sat[END_REF], operator selection [START_REF] Merz | Fitness landscape analysis and memetic algorithms for the quadratic assignment problem[END_REF], and parameter tuning [START_REF] Hutter | Automatic algorithm configuration based on local search[END_REF].

Local Optima Networks

The Local Optima Network (LON) model [START_REF] Ochoa | A study of nk landscapes' basins and local optima networks[END_REF] was introduced as a tool for studying the connectivity of local optima in a fitness landscape, and has subsequently shown proficiency in helping with explaining metaheuristic search dynamics [START_REF] Chicano | Local optima networks, landscape autocorrelation and heuristic search performance[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF][START_REF] Mcmenemy | How perturbation strength shapes the global structure of tsp fitness landscapes[END_REF]. We define the components of a LON before describing the model as a whole.

Nodes. The set of nodes LO are the local optima, meaning that a node lo i has superior fitness with respect to the entire neighbourhood. Formally: ∀n ∈ N (lo i ) : f (lo i ) ≤ f (n) (assuming minimisation) where N (lo i ) is the neighbourhood and n is a neighbour.

Edges. An edge is delineated between two nodes if the probability of "escape" from the source local optimum to the destination is greater than zero. The "escape" is defined with respect to a chosen search operation (or sequence of operations). The edge is weighted with the probability as w ij . Formally local optima lo i and lo j comprise the source and destination of an edge respectively iff w ij > 0. In this work, sampling is used; as a result, nodes are not necessarily associated with their complete set of potential edges.

Local optima network (LON).

A local optima network, LON = (LO, E), consists of nodes lo i ∈ LO which are the local optima, and edges e ij ∈ E between pairs of nodes lo i and lo j with weight w ij iff w ij > 0. We remark here that w ij may be different than w ji ; it follows that two weights are needed and that a LON is an oriented and weighted graph.

Fractal Dimension

The notion of a fractal dimension for patterns was conceived by [START_REF] Mandelbrot | Stochastic models for the earth's relief, the shape and the fractal dimension of the coastlines, and the number-area rule for islands[END_REF] and is defined as a complexity index which captures how the detail in the pattern changes with resolution used to measure it. The fractal dimension can be computed as the ratio between the logarithm of the detail and the logarithm of the scale used:

fractal dimension = log(detail) log(scale) (1)
To understand what the fractal dimension of a shape means we can begin by revisiting the familiar shapes associated with the topological dimension: a onedimensional line; a two-dimensional square; a three-dimensional cube.

We can observe the relationship between scale and detail for a square in Figure 1. Looking first at Figure 1a where the length scale m used to measure is one (the length of one side of the square) the detail measured is precisely one square. Moving onto Figure 1b we observe that a length scale of m = 1 2 is used here (this is a scaling factor of two because the resolution is twice as fine). That results in the measurement of four smaller copies of the larger square. The scale is two and the detail is four. Similarly, when m is one-quarter of the length of a side of the square (scaling factor of four; see Figure 1c) this results in sixteen copies of the larger square being measured, giving a scale of four and detail of sixteen. The relationship 4 x = 16 where x is dimension can be transformed into log( 16) log(4) , i.e. the ratio between detail and scale which is two in this case. The square is two-dimensional because for any scale the detail observed will be scale 2 .

For some patterns the exponent x is not an integer but rather somewhere else on the real number line. In this case, the way detail changes with resolution cannot be captured with topological dimension. An illustrative example of this can be seen in Figure 2 with the Sierpinski Triangle. Figure 2a shows that when a scaling factor m of one is used we accordingly measure the complete pattern. If we increase the resolution twofold as in Figure 2b three smaller copies of the large triangle are now measured. Recalling that fractal dimension can be obtained by solving for x the equation scale x = detail we observe that x is not an integer here. The equation is 2 x = 3 which results in a fractal dimension of x = ∼1.585.

Fractal dimensions can have efficacy in obtaining spatial and geometric information about real-world systems. They have been used, for example, in engineering for detecting cracks in plate structures [START_REF] Hadjileontiadis | Crack detection in plates using fractal dimension[END_REF]; in biology for characterising the tortuosity of animal trails [START_REF] Dicke | Using fractal dimensions for characterizing tortuosity of animal trails[END_REF]; and also in medicine for characterising mammographic patterns [START_REF] Caldwell | Characterisation of mammographic parenchymal pattern by fractal dimension[END_REF] and detecting colon cancer [START_REF] Esgiar | Fractal analysis in the detection of colonic cancer images[END_REF]. In our study we are computing fractal dimensions on LONs to obtain spatial complexity information about fitness landscapes. A widely-used method to estimate fractal dimension for a complex network is the "box counting" algorithm [START_REF] Song | Self-similarity of complex networks[END_REF]. This "boxes" together nodes which are within m network edges of each other, aiming to describe the network using as few "boxes" as possible. The parameter m is the scale of measurement used and serves as the denominator in Equation 1 to obtain fractal dimension alongside the number of "boxes" required to cover the network, which is the amount of detail observed. In stage one of the procedure, "centre" nodes are initially identified as those which are the best connected in the network. Nodes which are at a distance of no more than m edges to the centre node are then marked as "covered" and are added to the "box" associated with the centre. The process continues until all nodes are either "covered" or they are centre nodes. That means wherever a node cannot be "covered" with respect to any of the centres, it becomes a centre itself. In stage two the central distances for all nodes are calculated; this is -for each node -the closest distance to a centre. Following that, the "box" membership identity of each non-centre node is switched to that of a neighbour which is closer to a centre node. The original node is removed.

At the end of the process the number of "boxes" needed to cover the network completely, which we refer to as mb, is the number of detail units observed when using the resolution scale m. We can obtain the fractal dimension for the network by inserting mb and m into Equation 1:

fractal dimension = log(mb) log(m) (2)
Fractal complexity in local optima networks has been calculated previously using box counting (Thomson et al., 2018a,b). The box counting algorithm was altered in Thomson et al. (2018b) to specialise to LONs. For two nodes to be "boxed" as a single "unit" of detail they must either be a single edge apart or they are within m edges of each other and they also have a fitness distance less than a set threshold .

A subsequent study proposed additional mechanisms for computing and therefore defining the fractal dimension of a LON (Thomson et al., 2018a). A box counting variant which was introduced which used LON edge weights during the process. In a LON edge weights represent the probability that a search path between the local optima will be followed. The box counting variant used as the criteria for "boxing" that two nodes have a single edge between them which is weighted with a probability greater than β. The authors referred to values obtained using this method as probabilistic fractal dimensions.

In real-world complex systems a single fractal dimension can sometimes be insufficient to capture the complexity [START_REF] Mandelbrot | A multifractal model of asset returns[END_REF]. Monofractal analysis such as the box counting described earlier is based on the assumption that fractal complexity is roughly uniform in the pattern. Some networks have been found to be multifractal [START_REF] Song | Multifractal analysis of weighted networks by a modified sandbox algorithm[END_REF][START_REF] Furuya | Multifractality of complex networks[END_REF]. A multifractal algorithm has been used on LONs in a prior study (Thomson et al., 2018a) and we deploy this in our experiments. The process produces a spectrum of fractal dimensions for a single pattern (LON in our case). Details and pseudo-code for the algorithm are provided later on in Sections 3.3 and 4.3.

Methodology

The Quadratic Assignment Problem

Our analysis is conducted on the much-studied Quadratic Assignment Problem (QAP) [START_REF] Lawler | The quadratic assignment problem[END_REF] which is often used in fitness landscape analysis [START_REF] Merz | Fitness landscape analysis and memetic algorithms for the quadratic assignment problem[END_REF][START_REF] Merz | Advanced fitness landscape analysis and the performance of memetic algorithms[END_REF][START_REF] Daolio | Communities of minima in local optima networks of combinatorial spaces[END_REF][START_REF] Pitzer | A comprehensive survey on fitness landscape analysis[END_REF][START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF]. A QAP instance is specified with a distance matrix and a flow matrix. An entry in the distance matrix, D l1,l2 is the distance between two locations: d(l 1 , l 2 ). In the flow matrix this is the flow between two items: f (i 1 , i 2 ). Solutions are encoded as a permutation of length N , and are the allocation of N items to N locations. Fitness of a solution is the product of distances and flows between the locations and items according to the permutation and the aim is minimisation. The fitness function, g, for a solution x is then g

(x) = N i=1 N j=1 D ij F ij , ∀x ∈ S.

Constructing Sampled LONs

LON sampling algorithms are generally augmented on top of an existing optimisation algorithm. We align with this trend here, opting for a recently-introduced construction algorithm which joins an ILS with LON logging for QAP [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF].

The ILS algorithm is run r times from independent random starting solutions. In the ILS process, the local optimisation is a pairwise exchange of items, with the perturbation being k pairwise exchanges. Whenever the ILS Sampling detects no improving moves from the current solution, the solution is added as a local optimum node -this is an approximation of the true structure, because the algorithm does not consider the existence of saddle points. Only improving or equal fitness local optima are accepted. It follows that local optima plateaus might be explored, although not exhaustively, during sampling. When there is a local optima plateau, this is not collapsed by default; plateaus are sometimes collapsed to facilitate the extraction of certain LON features and to assess the neutrality present. Each local optimum encountered during search is stored in the set of nodes LO alongside its fitness, and if two optima l1 and l2 are connected by an ILS cycle (local search followed by k perturbations) during the search, an edge e l1,l2 is stored in the LON edge-set, E. The nodes are edges logged during the r runs are joined to form a single local optima network for the problem instance. All parameters for the algorithm are stated later on in Section 4.2.

Fractal Analysis Algorithms

As stipulated in Section 2.3 the standard approach for calculating and defining fractal dimension of a complex network is with a box counting algorithm [START_REF] Song | Self-similarity of complex networks[END_REF][START_REF] Song | Origins of fractality in the growth of complex networks[END_REF]. This process iteratively "boxes" together nodes iff the distance d(n 1 , n 2 ) < m, i.e. the nodes are < m edges apart. The parameter m provides the scaling factor which is used to compute the fractal dimension of the network alongside the associated detail observed when using that scaling factor. Empirically the detail is defined as the number of "boxes" needed to completely cover the network, taken as a proportion of the network size.

We mentioned that the box counting algorithm has been specialised for the specific case of a local optima network previously (Thomson et al., 2018b). In that they allowed nodes to be "boxed" if either the distance

d(n 1 , n 2 ) = 1 or d(n 1 , n 2 ) < m and also |f (n 1 ) -f (n 2 )|< where f (n x )
is the fitness of node x and is the maximum fitness difference between nodes n 1 and n 2 .

Multifractal Dimension Analysis

The process for calculating multifractal dimensions is different to standard box counting and it produces a spectrum of fractal dimensions for a single pattern. One approach is called the sandbox algorithm [START_REF] Liu | Determination of multifractal dimensions of complex networks by means of the sandbox algorithm[END_REF] where several nodes are randomly selected to be sandbox "centres". Members of the sandboxes are computed as nodes are r edges apart from the centre c. After that the average sandbox size is calculated. The procedure is replicated for different values of r which is the sandbox radius. To facilitate the production of a dimension spectrum the whole process is repeated for several arbitrary real-valued numbers which supply a parameter we call q. The sandbox algorithm was specialised and modified to suit LONs in a prior study (Thomson et al., 2018a) and this is the process we use for our fractal analysis experiments. In our version of the algorithm a node n can be included in the "sandbox" of a central node end for 24: end for At the end of each "sandboxing" iteration conducted with particular values for the parameters q, r and , the associated fractal dimension is calculated:

c if either the distance d(n, c) = 1 or d(n, c) = r -1 and |(f (n) -f (c))|< . Pseudocode is given in Algorithm 1.
fractal dimension = log(detail q-1 ) (q -1) * log(scale) (3)
where detail is the average "sandbox" size (as a proportion of the network size), q is an arbitrary real-valued value, and scale is r dm , with r being the radius of the boxes and dm the diameter of the network. We use this process as the foundation for obtaining multifractal dimensions for LONs in this work. In addition we separately implement a modified version of it where the metaheuristic path probabilities encoded in LON edge weights are used in the calculations. Specifically, for a node n to be a member of the sandbox with centre c there must be one of two situations: either there is a single edge between n and c (of any probability; this is to guarantee boxing momentum), or there is an edge between a direct neighbour of n and c which is weighted with a probability greater than a specified threshold β. This element was implemented with the motivation that nodes which are in close proximity to a probable path towards the central node should be included in the sandbox. We remove the fitness distance check for this algorithm variant and instead of a set of values for "sandbox" radii the sandboxes are of a fixed width, r = 2. The rest of the algorithm remains the same and a spectrum of fractal dimensions is produced. To differentiate the results produced by this particular algorithm variant in the following Sections we refer to the fractal dimensions obtained by this method as probabilistic fractal dimensions. The algorithm which does not include the probability constraint but instead includes "sandbox" radius variation as well as fitness distance constraints produces values which we refer to as deterministic fractal dimensions. The parameters for both algorithms described are stated in Section 4.3.

Experimental Setup

Instances Used

All instances used are from the QAPLIB benchmark library for QAP, the Quadratic Assignment Problem Library (QAPLIB) [START_REF] Burkard | Qaplib-a quadratic assignment problem library[END_REF]. We cap the maximum problem size at 50 due to the computational expense [START_REF] Liu | Determination of multifractal dimensions of complex networks by means of the sandbox algorithm[END_REF] associated with multifractal analysis of large networks. It follows that further study is needed in order to confirm any findings on larger problem instances. We additionally remove the "esc" instances from the group because their LONs have very few distinct fitnesses due to large amounts of neutrality present. The resultant set consists of 85 problems, with the problem sizes ranging from 12 to 50. In all cases the global optimum is known.

The nature of QAP instances can commonly be characterised into one of four classes [START_REF] St Ützle | Iterated local search for the quadratic assignment problem[END_REF]: uniform random distances and flows; random flows based on grids; real-world; and random "real-world like", which are not real-world but mimic distance and flow patterns seen in real-world presentations of QAP. Table 1 shows the QAPLIB instances used in the experiments and present them in these four categories. Numbers which form part of the instance names indicate the problem size, i.e. number of locations and flows and the length of a permutation solution.

Construction of Local Optima Networks

For each QAP instance we construct a local optima network. As stipulated in Section 3.2, this is done by using an ILS algorithm which has been augmented with LON logging mechanisms. The LON logging amalgamates the unique nodes and nug{12, 14, 16a-b, 17, 18, 20, 21, 22, 24, 25, 27 15b, 20b, 25b, 30b, 35b, 40b, 50b} edges from 200 ILS runs into a single network. Each run terminates after 10,000 iterations without an improvement. This is a deliberately lenient condition which was chosen with the motivation that ILS runs should converge to a natural stalling point. The parameter setting, however, means that some LONs (those associated with the larger problems in the instance set) are built over a number of hours; this computational cost is a limitation to our approach. We argue, however, that the benefit of the insight gained through our method outweighs the cost. The remaining ILS parameters and setup are detailed shortly in Section 4.4.

, 30} | scr{12, 15, 20} | sko{42, 49} | tho{30,40} | wil{50} real-world bur{a-h} | chr{12a-c, 15a-c, 18a-b, 20a-c, 22a-b, 25a} | {els19} | kra{30a-b, 32} | lipa{20a-b, 30a-b, 40a-b, 50a-b} | ste36{a-c} real-world like tai{12b,

Fractal Analysis

In contrast to traditional monofractal analysis, to generate multifractal dimensions for the LONs a range of arbitrary real-valued numbers is needed. We set these as q in the range [3.00, 8.90] in step sizes of 0.1. The number of "sandbox" centres in each iteration is set at 50 and the choice of these centres is randomised. As mentioned, in our deterministic multifractal algorithm variant fitness distance is considered in order to specialise to local optima networks. The comparison between two fitness values is conducted through logarithmic returns:

fitness difference = log(f 1 /f 2 ) (4)
where f 1 and f 2 are the fitnesses of two local optima at the start and end of a LON edge. We take the absolute value of the computed fitness difference because if f 1 f 2 , the result of Equation 4 is negative. This value can then be compared with a set threshold, . A range of ten values is used for that algorithm: ∈ {0.01, 0.19} in step sizes of 0.02. Another essential element of deterministic multifractal analysis is the sizes (radii) for the sandboxes. For these we use values in the range r ∈ {2, diameter -1} where diameter is the LON diameter.

For the probabilistic multifractal algorithm variant the fitness constraint is not used and the sandboxes are of a fixed width, r = 2. The probability threshold parameter β must be chosen. Recall that β sets the minimum edge weight between two nodes. After preliminary runs it was noted that if β was set as greater than the minimum weight in the weights distribution then little-to-no "boxing" occurred. For this reason β is set as the minimum weight present in the distribution.

An important note. We note here that 32 out of the 85 LONs had only a single edge weight present throughout the network. Recall that the probability-based boxing process outlined in Section 3.3.1 stipulates that nodes can be boxed together when either there is a single edge between n and c (of any probability; this is to guarantee boxing momentum), or there is an edge between a direct neighbour of n and c which is weighted with a probability greater than a specified threshold β, which is set as the lowest weight in the network. As a consequence, when all weights are equivalent, then no boxing based on probability will occur at all -no pairs of nodes will pass the acceptance condition that their connecting edge has a weight greater than the minimum weight seen in the LON. This renders these particular networks ineligible for probabilistic fractal analysis under these conditions. Consequently, results which pertain to probabilistic dimensions consider the 53 eligible LONs and their features, while those pertaining to deterministic fractal dimensions consider all 85 LONs.

Metaheuristic Performance

To obtain algorithm performance information with which to compare the LON features we use two search algorithms for the QAP. St ützle introduced iterated local search (ILS) variants for state-of-the-art performance on the QAP [START_REF] St Ützle | Iterated local search for the quadratic assignment problem[END_REF]. We use his ILS configured as follows: first-improvement pairwise exchanges for local search; 3n 4 exchanges for perturbation; accepting only improving local optima; and terminating when the global optimum is found or after 100 iterations. Taillard's Robust Taboo Search (ROTS) [START_REF] Taillard | Robust taboo search for the quadratic assignment problem[END_REF] is also a competitive heuristic for the QAP. This a best-improvement pairwise exchange local search with a variable-length tabu list tail. For each facility-location combination, the most recent point in the search when the facility was assigned to the location is retained. A potential move is deemed to be "tabu" (not allowed) if both facilities involved have been assigned to the prospective locations within the last s cycles. The value for s is changed randomly, but is always from the range [0.9n, 1.1n]. A run terminates when the global optimum is found, or after 100 iterations.

We run the ILS and the ROTS in these configurations on each QAPLIB instance 100 times from different starting solutions. As a measure for their performance we define the performance gap p as follows:

p = f (alg) f (opt) (5)
where f (alg) is the fitness obtained by the algorithm and f (opt) is the fitness of the global optimum. In this way, a "solved" run will output "1" and lower values are closer to the optimal fitness. For each QAP instance we report the mean p over 100 runs. In the results that follow, p(ILS) is this value for iterated local search and p(ROT S) is for robust tabu search.

LON Features

Features are extracted from the local optima networks. Deterministic fractal dimension sets are calculated for each of the 85 LONs considered. In those sets there are 60 * (diameter -2) * 10 dimensions, where 60 is the number of (arbitrary) q values, diameter is the LON diameter (which differs between LONs), and 10 is the number of values for fitness distance threshold . As we recall from Section 4.3 53 of the 85 LONs are eligible for probabilistic fractal analysis. Those 53 have sets of probabilistic fractal dimensions calculated in addition to the deterministic ones. In each set there are 60 dimensions (one for each value of q). This is less than there was when using the deterministic multifractal algorithm and this is because the probabilistic variant does not consider parameter ranges for "sandbox" radius and does not include fitness distance calculations. The statistics we draw from the fractal complexity data are: the minimum fractal dimension, the median fractal dimension, the maximum fractal dimension, the range of fractal dimensions (calculated as the difference between the largest and smallest values), and the number of distinct fractal dimensions. The latter two capture the degree of multifractality present.

Also considered in the experiments are other LON features which are not related to fractal dimension values:

1. number of local optima; 2. mean fitness of sampled local optima in the LON; 3. fitness range; 4. fitness of sinks ("sinks" are nodes which have no out-going edges, i.e. the algorithm used to construct the LON became trapped there); 5. extent of meta-neutrality, which is neutrality at the local optima level, computed as meta-neutrality = number unique fitnesses number local optima ;

6. mean out-degree.

Regression Model Setup

We build algorithm performance models using LON features for predictors and the performance of competitive metaheuristic algorithms as the response variables. The aim is clarifying how LON features can contribute to explaining or predicting algorithm proficiency, paying particular attention to the fractal nature of the LON. In pursuit of that we conduct linear and random forest regressions. The number of observations we have is relatively small -85 for the deterministic dimensions and 53 for probabilistic -so we use random repeated subsampling cross-validation for obtaining model statistics. This is conducted for 10,000 iterations with a training-test split of 80-20. The random forest regression uses 500 trees. Predictors are standardised (due to different value ranges) as follows: p = (p-E(p)) sd(p) , with p being the predictor in question. The model statistics we focus on are R 2 , which captures the amount of variance in the response variable which can be explained using the predictor set, and mean squared error, which expresses the mean squared difference between the model-estimated values and the actual values.

For the random forest models, variable importance rankings and values are reported. The values are calculated as the reduction in decision tree node impurities when splitting on the variable and are averaged over all 500 trees used in the regression. Node impurities are measured with the residual sum of squares.

The non-fractal LON predictors used in the models are the mean fitness; fitness range; fitness of sinks; extent of meta-neutrality; out-degree; and the number of global optima. For the deterministic fractal dimensions, we include the minimum fractal dimension and median. In the probabilistic case, these two are replaced with the fractal dimension range and number of unique dimensions.

Results

Distribution Analysis

In Figures 3 and4, box-plots convey information about the fractal dimensions calculated on the local optima networks. Each box contains values for LONs associated with a particular QAPLIB instance class -those are indicated on the x-axis labels. Only a sub-set of the instance classes which are involved in the central experimentation of this Chapter are considered in these plots. I chose these groups because displaying their distributions alongside each other illustrates evident visual differences between these particular classes. Also provided in the Figures as accompanying text for each box is the performance of iterated local search on the QAP instances associated with those LONs; this is the performance metric p(ILS).

In Figures 3a and3b the distributions concern the median LON fractal dimension which is associated with using the deterministic and probabilistic methodologies, respectively. In the case of the deterministic fractal analysis, this is the median value computed over all of the dimensions produced under these conditions; each dimension is the output resulting from using a different combination of the fractal analysis parameters q, r and . The probabilistic median is computed from the spectrum of dimensions associated with the range of values for q.

In both Figure 3a and Figure 3b, the "lip" class of LONs seem to have the highest values and the "had" group have the lowest. On both plots, the highest value belongs to the "lip" category and the lowest to "had". Notice that in 3b the "lip" and the "nug" instances -whose LONs generally have the highest fractal dimensions in this plot -also have higher values of p(ILS). As stipulated in Section 4.4, values like these reflect that metaheuristic performance was of lower quality. With deterministic analysis, the "lip" group have the largest variation, while the "had" LONs have among the smallest; with probabilistic dimensions (Figure 3b), "had" have the largest and "lip" the smallest. Deterministic fractal dimensions appear to be higher than probabilistic fractal dimensions.

Consider now the range of fractal dimensions in the deterministic and probabilistic spectra calculated for the LONs, which are given in Figures 4a and4b.

The range of fractal dimensions for a LON is a way to quantify the extent of multifractality present and is calculated as maximum valueminimum value with respect to the complete set of fractal dimensions produced using either the deterministic or probabilistic paradigm. Also provided is the average ILS performance, p(ILS), for the QAP instances included in the classes.

Looking at the two plots and noting the different scales used for them, it seems clear that the probabilistic dimension calculation process lends to more compact ranges. This is intuitive: the conditions are stricter for measuring "boxes" during Fig. 4: Boxplots for the distributions of the range of fractal dimension. Here the range is the smallest dimension for that LON subtracted from the largest. Each box displays values for LONs extracted from one of four QAPLIB instance classes, as indicated in the legend. The black line within the boxes is the median value the dimension calculation process. Let us consider in both plots the levels of the black lines (which indicate the distribution median). The "had" group has the lowest in 4a and the "lip" group has the highest. That hints that the degree of multifractality in the "lip" group is the most pronounced among the four, and it is the least pronounced in the "had" group. The previous plots told us that "lip" LONs had the highest dimensions, and "had" showed the lowest. It follows that the degree of deterministic multifractality might be associated with lower fractal dimensions. For 4b though, "had" LONs have the highest ranges of dimension and "lip" have the lowest -the opposite trend to the deterministic dimensions. With respect to algorithm performance, we can that the "lip" LONs, associated to problems with the lowest metaheuristic performance (p(ILS)), appear to have a higher extent of deterministic multifractality and a lower extent of probabilistic multifractality. In Figure 4b, the two problem groups with the best ILS performance have the widest ranges of values for dimension (i.e. amount of multifractality) of the four categories.

Visualisation

Visual analysis of LONs provides valuable insight into algorithm performance and problem structure, and can augment more empirical or statistical findings [START_REF] Ochoa | Additional dimensions to the study of funnels in combinatorial landscapes[END_REF]. We begin with visualisation before moving onto correlation analysis (Section 5.3) and machine learning models (Section 5.4) thereafter.

Figure 5 shows two partial LONs, each for a different QAPLIB instance. Only the fittest 10% of local optima are plotted for visual clarity. Global optima are red squares and all other nodes are grey circles. The node sizes are proportional to the incoming strength to that node, which is the weighted incoming degree. These two LONs were selected from the "had" and "lip" instance classes because the former have lower fractal dimensions and also a lower degree of deterministic multifractality than the latter. These two instances chosen have the same problem size, N = 20, and similar numbers of local optima.

(a) partial LON of "had20" QAPLIB instance which has median fractal dimension 2.975 and associated tabu search performance 1.011 (b) partial LON of "lipa20b" QAPLIB instance which has median fractal dimension 4.015 and associated tabu search performance 1.154 Fig. 5: partial local optima networks for two QAPLIB instances; only local optima which are in the fittest 10% are shown. Global optima are square and red; all others are grey circles. The size of the nodes captures the incoming "strength" to the node in the LON, i.e. the weighted in-degree

In accordance with the higher fractal dimensions, the algorithm performance is lower on the "lip" group of problems. Using as a performance measure the obtained fitness (as a proportion of the global fitness), robust tabu search averaged 1.096 on the "lip" instances. For the "had" group this was 1.011. Our task in this Section of the results is to seek explanation in the networks concerning the algorithm performance differences while also paying particular attention to how their fractal nature relates to what is visually seen in the structure.

The median fractal dimension for "had20", plotted in Figure 5a, is 2.975; for "lipa20b" it is 4.015. The range of fractal dimensions for "had20" is around 63, and is around 49 for "lipa20b". An evident difference in the two Figures is the number and connectivity of global optima -Figure 5a shows that the "had20" LON has many, and they appear to be densely connected to other nodes. Contrarily, the "lipa20b" LON in Figure 5b has a single global optimum, which seems to be more sparsely connected within its network. Also noteworthy is the relative sizes of the non-optimal (grey) nodes. In Figure 5a there are many large nodes which are suboptimal and they have access to the global optima. Figure 5b is not the same; in fact, many of the nodes which are one step from the global optimum are very small indeed. That tells us that these nodes have small incoming degree which might hinder ascension through fitness levels during optimisation. These grey nodes are also not well-connected to each other. The opposite is true for the other network. In the "had20" LON (Figure 5a), connectivity is so dense in the promising local optima region that visually tracking paths is impossible.

Let us now view the Figures using an algorithm performance explanation lens. Of course, the number of global optima matters and so does the accessibility of them. The "lipa20b" global optimum has many incoming edges but most of these are sourced from nodes which have low incoming degree themselves. It follows that the global optimum is less accessible. The "had20" LON, which is highly populated with edges in this promising landscape region, is probably easily solvable in part because when an algorithm reaches one of the large grey nodes (this should be likely because they have high incoming degree) there is an abundance of paths to a global optimum. The same trends are present when comparing the two networks in Figures 6a and6b.

(a) partial LON of "had18" QAPLIB instance which has median fractal dimension 3.175 and associated tabu search performance 1.011 (b) partial LON of "nug16b" QAPLIB instance which has median fractal dimension 4.090 and associated tabu search performance 1.055 Fig. 6: partial local optima networks for two QAPLIB instances; only local optima which are in the fittest 15% are shown. Global optima are square and red; all others are grey circles. The size of the nodes captures the incoming "strength" to the node in the LON, i.e. the weighted in-degree These are the partial LONs of the "had18" and "nug16b" QAPLIB instances. Figure 6a shows "had18", which has a lower median fractal dimension (3.175) and better tabu search performance on the underlying problem (1.011) when compared with "nug16b" shown in Figure 6b, which has a median fractal dimension of 4.090 and tabu search performance of 1.055. Surveying the two figures, we can again visually account for the discrepancy in fractal dimension and algorithm performance by looking at the spatial complexity. Although the LON of "nug16b" has more global optima (in red), edges appear less uniformly distributed in their vicinity when compared to the LON of "had18". In addition we notice that some nodes which are one step from a global optimum in Figure 6b are small in size. This tells us that they have low incoming degree and that the probability of search paths reaching them is small. As a consequence potential routes towards the global optima may be missed by algorithms.

Correlation Analysis

Figures 7 and8 show pairwise correlations between variables. These are Spearman rank coefficients, which are more appropriate to use where variables are not linearly related. Included are p(ILS) and p(ROT S) on the QAP instances, alongside the proposed fractal dimension features (unique FDs, range FD, median FD, max FD and min FD) -in Figure 7, these concern deterministic dimensions; in Figure 8, they are probabilistic. Also shown are LON features which are not associated with fractal complexity -these were introduced in Section 4.5. We additionally include LON features from previous works (Thomson et al., 2018a;[START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF][START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF] for comparative purposes. These are the LON diameter; the number of compressed local optima (after connected LON nodes of the same fitness are compressed together -labelled as comp.opt in the Figures); the correlation between the fitnesses of neighbours in the LON (fit.fit.corr); the number of sink nodes present (sinks); and the sub-optimal sink strength (that is, the total incoming edge weight to any sub-optimal sink nodes -so.strength). In addition, fractal dimensions excerpted from arbitrary points on the multifractal spectra are considered as features -these are arb.dfd1 and arb.dfd2 in Figure 7 and arb.pfd1 and arb.pfd2 in Figure 8. The approach of taking an arbitrary excerpt from the spectrum and using it as a feature was taken in the previous work on multifractality in LONs -its inclusion here facilitates a comparison between previous features and the proposed ones.

In particular we are interested in the correlation between fractal features of the LON and algorithm performance variance on the associated combinatorial problem. The intersections between the p(ILS) column and the fractal feature rows in Figure 7 reveal moderate positive correlations between them in the case of the fractal dimension range, median, maximum, and minimum -as well as the two dimension excerpts, arb.dfd1 and arb.dfd2. For all of these the associated p-value is less than 0.001. We notice that the correlations are stronger than the p(ILS) correlations with other LON features such as mean.fitness, fitness.sinks, LO.neutrality, edges, diameter, comp.opt, optima, fit.fit.corr, and assortativity. They are also slightly stronger than the correlations between p(ILS) and sinks and so.strength. The correlations with dimensional summary statistics such as med FD appear slightly larger than arb.dfd1 and arb.dfd2, which are the fractal dimension features calculated using Next we will consider the correlation plot which includes probabilistic fractal dimension variables in Figure 8. Only two of the proposed fractal features show a moderate negative relationship with p(ILS) -these are those related to the extent of multifractality, i.e., unique FDs and range FD. Those correlations have associated p-value of less than 0.001. These two also have a weak negative correlation with p(ROT S) -for unique FDs, the p-value is 0.0118; for range FD, it is 0.0776. Again, the correlations between the other LON features and p(ILS) and p(ROT S) are diminutive, with the exception of out-degree and clusteringcoef. Observe that the fractal dimensions which were arbitrarily excerpted from the multifractal spectra, arb.pfd1 and arb.pfd2, are far less correlated to p(ILS) than two of the fractal features proposed in this work (that is, unique FDs and range FD).

Algorithm Performance Regression Models

Deterministic Fractal Dimension Features

Table 2 contains regression model statistics whose values are estimated over 10,000 random repeated subsampling iterations. Each row represents a particular model setup. The response variable is shown in the second column. The R 2 and mean squared error are given. We can see from the R 2 values that random forest regres- 9. The values are averaged over 10000 iterations of random repeated subsampling cross-validation. For explaining p(ILS), the fitness of sinks is most important. The median fractal dimension is second most important. These two predictors have importance values noticeably higher than the rest, although even the lowest predictors, LO neutrality and minimum fractal dimension, still have importance values around 0.009, which is around half the value of the highest. Fitness of sinks is again the most important factor in Figure 9b, which is the tabu search response model setup. This is followed by other fitness-based features in second, third, and fourth place -LO neutrality, mean fitness, and fitness range. The median fractal dimension contributes moderately well, ranking fourth out of eight features. The position of minimum fractal dimension is last, but even so, it does contribute to the model.

Probabilistic Fractal Dimension Features

In Table 3 is model statistics where the predictor set includes probabilistic fractal dimension features instead of the deterministic ones seen in Table 2. This is followed by the associated random forest predictor rankings in Figure 10. The random forest p(ILS) model setup is rather weak with respect to the R 2 estimate. Indeed, in this measure it is weaker than the equivalent model which used deterministic dimensions. It should be reiterated at this point that the data-set is composed of fewer observations here than in the previous models (Tables 2 and9). There are 53 observations here, compared with a previous 85. This might impact the formulation of a well-fitting model. Nonetheless both setups with p(ILS) as the response variable have markedly lower mean squared errors than their ROTS counterparts. This is also true in Table 2. Back in Table 3, the p(ROT S) models have higher mean squared errors but the random forest model is definitely the strongest with respect to search algorithm explanation, with around 58% being accounted for by the predictors. Although a smaller portion of variance is explained in the p(ILS) models, the low mean squared errors are encouraging in accuracy terms. In the predictor rankings, seen in Figure 10, we draw your attention to the two fractal dimension metrics in the iterated local search plot (Figure 10a): they are among the most important predictors, ranking third and fourth. Their values are betweeen 0.005 and 0.006; for comparison, the importance value of the most dominant predictor, the number of global optima, is just above 0.007. The fractal dimension predictors are among the least important in the tabu model setup seen in Figure 10b. Instead, the strongest predictors for p(ROT S) appear to be relating to the local optima level fitness distribution: mean fitness, LO neutrality, fitness of sinks, and fitness range. In fact, these four form a distinct group on the plot, far higher in importance than the remaining four (which include the fractal dimension features). Nevertheless, the lower group are not useless: their values are approximately in the range 0.04 -0.07; the more important features have values between around 0.155 and 0.171.

Conclusion

We conducted multifractal analysis on the local optima networks (LONs) associated with a benchmark combinatorial optimisation problem library, QAPLIB. The QAPLIB instance set was more than three times the size of the set used in a prior study (Thomson et al., 2018a) and raised the considered problem sizes from N ≤ 28 to N ≤ 50. A recent and refined LON construction algorithm [START_REF] Ochoa | Perturbation strength and the global structure of qap fitness landscapes[END_REF] was used to build the LONs. Relationships between fractal dimension features of LONs and algorithm performance by iterated local search (ILS) were established using correlation analysis, visual analysis tools, and linear and random forest regression with random repeated subsampling cross-validation. The results showed that the extent of multifractality and the highness of values in the dimension spectrum can contribute towards partially predicting or explaining ILS algorithm performance. Features of the fractal dimension distribution for the LONs also displayed individual pairwise correlations to ILS algorithm performance. Fractal dimension features in LONs were less important for predicting tabu search but could still contribute some information. Sampled fitness levels in the LON were more important in these models.

A limitation to our approach is that the features are computed from sampled LONs, whose characteristics can alter markedly with a different sampling effort [START_REF] Bo Żejko | Local optima networks in solving algorithm selection problem for tsp[END_REF]. Nevertheless, the LON features can contribute towards explaining algorithm performance within regression models -it follows that they are useful, even if the sample illustrates a certain version of the fitness landscape. In addition, there is no alternative to sampling when analysing QAP LONs of moderate size (greater than N = 11 according to [START_REF] Daolio | Communities of minima in local optima networks of combinatorial spaces[END_REF]). Another consideration is the random selection of box centres in the sandbox algorithm for multifractal analysis. We argue that the number of algorithm iterations -each of which contains a random selection of centres and produces its own fractal dimension -should mitigate the variation induced by random selection. In addition, the resulting fractal dimension features help to explain metaheuristic performance in statistical analysis over 85 observations. This implies that the randomness inherent to the approach does not affect the empirical usefulness of the computed fractal dimensions.

The present study could serve as a foundation for further work within this research avenue which remains untapped. In particular, we would like to expand the maximum size of the problems studied, as well as venturing to other domains and to constrained problems. Finally, we conclude with a remark concerning our interest in studying the relationship between perturbation strength used to generate the LONs, and the calculated fractal dimensions of that LON.

Fig. 1 :

 1 Fig. 1: The relationship between detail and scale for a square; m is the length of scale used for measurement

  Fig.2: The relationship between detail and scale for a fractal with topological dimension two and fractal dimension of ∼1.585

Fig. 3 :

 3 Fig. 3: Boxplots for the distributions of median fractal dimension. Each box displays values for LONs extracted from one of four QAPLIB instance classes, as indicated in the legend. The black line within the boxes is the median value

  Fig. 8: Spearman correlations between pairs of variables including p(ILS) and p(ROT S), fractal dimension metrics for the LONs, and other landscape features

Fig. 9 :

 9 Fig.9: Variable importance values for random forest models; the models include deterministic fractal dimension features as part of the predictor set.

Fig. 10 :

 10 Fig.10: Variable importance values for random forest models; the models include probabilistic fractal dimension features as part of the predictor set.

  Algorithm for Multifractal Analysis of a LON

Algorithm 1 Input: LON , q.values, radius.values, f itness.thresholds, number.centres Output: mean sandbox size 1: Initialisation: 2: centre.nodes ← ∅, noncentre.nodes ← all.nodes 3: mean.sandbox.sizes ← ∅ 4: for q in q.values do 5: for rd in radius.values do 6: for in f itness.thresholds do 7: centre.nodes ← RANDOM.SELECTION(all.nodes, number.centres) 8: sandbox.sizes ← ∅ 9: for c in centre.nodes do 10: number.boxed ← 0 11: for v in all.nodes do 12: d ← DISTANCE(c, v) 13: j ← DIFFERENCE(f (c), f (v)) 14: if ( d == 1 ) OR ( d == rd -1 and j < ) then: 15: number.boxed ← number.boxed + 1 16: end if 17: end for 18: sandbox.sizes ← sandbox.sizes ∪ {[number.boxed]} 19: end for 20: bs ← MEAN(sandbox.sizes) 21: mean.sandbox.sizes[q][rd][ ] ← bs 22: end for 23:

Table 1 :

 1 The QAPLIB instances used in the experiments

	class	instance names
	uniform random tai{12a, 15a, 17a, 20a, 25a, 30a, 35a, 40a, 50a} |
		rou{12, 15, 20}
	random grids	had{12, 14, 16, 18, 20} |

Table 2 :

 2 Summary statistics averaged with 10,000 iterations of random repeated subsampling cross-validation for explaining p(ILS) and p(ROT S). Predictors include deterministic fractal dimension LON statistics, as well as other landscape features such as number of local optima and fitness distribution measures This is likely because random forest trees are adept at considering non-linearities between variables. The amount of variance in the iterated local search and tabu search performance which can be explained using the predictors is higher in the random forest models. The mean squared error is very low in the case of the random forest regression which is explaining p(ILS). The strongest model in terms of R 2 is using random forest regression with p(ROT S) as the response, with around 61% of variance being explained using the landscape features. Less variance in p(ILS) response, around 48%, is explained using the same type of regression. This model setup does, however, have a much lower error rate than the associated p(ROT S) model. Now let us look at the random forest predictor importance rankings, which are provided in Figure

	type of regression	response variable	R 2	mean squared error
	linear	iterated local search performance 0.160	0.001
	linear	robust tabu search performance	0.161	0.027
	random forest	iterated local search performance 0.482	0.000
	random forest	robust tabu search performance	0.611	0.011
	sion produces a stronger model fit.

Table 3 :

 3 Summary statistics estimated with 10000 iterations of random repeated subsampling cross-validation for explaining p(ILS) and p(ROT S). Predictors include probabilistic fractal dimension statistics, as well as other landscape features such as number of local optima and fitness distribution measures

	type of regression	response variable	R 2	mean squared error
	linear	iterated local search performance 0.293	0.002
	linear	robust tabu search performance	0.332	0.017
	random forest	iterated local search performance 0.304	0.000
	random forest	robust tabu search performance	0.578	0.033