The fractal geometry of fitness landscapes at the local optima level

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 search performance, and visual analysis of LONs shows insight into this relationship.


Introduction
Fractals are patterns which contain parts resembling the whole (Mandelbrot 1972). Under this definition fractals are ubiquitous in the complex simplicity of nature, from microscopic blood vessel networks to the macroscopic pattern of the rings of Saturn. Nature and evolution seem to favour fractal design: using a pattern repeatedly allows replicability with very few instructions. The fractal dimension (Mandelbrot 1975) is a complexity index capturing how the detail in a pattern changes when one views it using a different resolution or scale. Fractal dimension analysis has been used in diagnostic imaging (detecting colon cancer (Esgiar et al. 2002); characterising images in mammography (Caldwell et al. 1990); characterising leukaemia cells (Mashiah et al. 2008), search and rescue [analysing the layout of victim location after building collapses (Saeedi and Sorensen, 2009)], and in engineering for the design of antenna (Werner and Ganguly 2003), 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 (Uahabi and Atounti 2017), which indicates high spatial complexity and convolution.
Fitness landscapes of some combinatorial optimisation problems have been viewed under a fractal lens (Weinberger and Stadler 1993). 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 (Weinberger and Stadler 1993) 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 (Zelinka et al. 2014;Locatelli 2005;Richter 2018) and some have emphasised the potential lying dormant in the largely untapped field of fractal analysis for landscapes.
A local optima network (LON) (Ochoa et al. 2008) 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 (Daolio et al. 2010;Daolio et al. 2011;Verel et al. 2011;Herrmann et al. 2016;Ochoa and Veerapen 2018;Ochoa et al. 2017).
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; Thomson et al. 2018). 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 (Burkard et al. 1997) 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 (Ochoa and Herrmann 2018) 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 crossvalidation; using intelligible predictors such as the extent of multifractality and the median fractal dimension).
The article is structured as follows: Sect. 2 contains the necessary background information to render this article self-contained; Sect. 3 details aspects of the methodology used; Sect. 4 gives the experimental setup, with Sect. 5 presenting the results; finally, Sect. 6 finishes the article with conclusions and directions for future work.

Fitness landscapes
A fitness landscape (Stadler 2002) 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 2 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 (Pitzer and Affenzeller 2012). Indeed, landscapes have been used to facilitate algorithm selection (Hoos et al. 2004), operator selection (Merz and Freisleben 2000), and parameter tuning (Hutter et al. 2007).

Local optima networks
The local optima network (LON) model (Ochoa et al. 2008) 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 (Chicano et al. 2012;Herrmann et al. 2016;McMenemy et al. 2018). 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: 8n 2 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 2 LO which are the local optima, and edges e ij 2 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 Mandelbrot (1975) 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: To understand what the fractal dimension of a shape means we can begin by revisiting the familiar shapes associated with the topological dimension: a one-dimensional line; a two-dimensional square; a three-dimensional cube. We can observe the relationship between scale and detail for a square in Fig. 1. Looking first at Fig. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 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 (Hadjileontiadis and Douka 2007); in biology for characterising the tortuosity of animal trails (Dicke and Burrough 1988); and also in medicine for characterising mammographic patterns (Caldwell et al. 1990) and detecting colon cancer (Esgiar et al. 2002).
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 (Song et al. 2005). 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 Eq. 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 Eq. 1: Fractal complexity in local optima networks has been calculated previously using box counting (Thomson et al. 2018a; Thomson et al. 2018). The box counting algorithm was altered in Thomson et al. (2018) 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 b. 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 (Mandelbrot et al. 1997). 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 (Song et al. 2015;Furuya and Yakubo 2011). 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 Sects. 3.3 and 4.3.

The quadratic assignment problem
Our analysis is conducted on the much-studied quadratic assignment problem (QAP) (Lawler 1963) which is often used in fitness landscape analysis (Merz and Freisleben 2000;Merz 2004;Daolio et al. 2011;Pitzer and Affenzeller 2012;Verel et al. 2018). 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 (a) (b) Fig. 2 The relationship between detail and scale for a fractal with topological dimension two and fractal dimension of $ 1.585

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 (Ochoa and Herrmann 2018). 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 nodethis 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. LON plateaus are collapsed in preparation for multifractal analysis in this work, and also to facilitate computation of the funnel features and the number of compressed local optima feature described in Sect. 4.5. 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 Sect.*4.2.

Fractal analysis algorithms
As stipulated in Sect. 2.3 the standard approach for calculating and defining fractal dimension of a complex network is with a box counting algorithm (Song et al. 2005(Song et al. , 2006. 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. 2018). 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 jf ðn 1 Þ À f ðn 2 Þj\ 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  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 c if either the distance dðn; cÞ = 1 or dðn; cÞ ¼ r À 1 and jðf ðnÞ À f ðcÞÞj\. Pseudocode is given in Algorithm 1.
At the end of each ''sandboxing'' iteration conducted with particular values for the parameters q, r and , the associated fractal dimension is calculated: where detail is the average ''sandbox'' size (as a proportion of the network size), q is an arbitrary real-valued value, and Fractal geometry in LONs 321 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 b. 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 Sect. 4.3.

Instances used
All instances used are from the QAPLIB benchmark library for QAP, the quadratic assignment problem library (QAPLIB) (Burkard et al. 1997). We cap the maximum problem size at 50 due to the computational expense  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 (Stützle 2006): 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 Sect. 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 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 Sect. 4.4.

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  15b, 20b, 25b, 30b, 35b, 40b, 50b} networks. The comparison between two fitness values is conducted through logarithmic returns: 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 6 f 2 , the result of Eq. 4 is negative. This value can then be compared with a set threshold, . A range of ten values is used for that algorithm: 2 f0:01; 0:19g 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 f2; diameter À 1g 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 b must be chosen. Recall that b sets the minimum edge weight between two nodes. After preliminary runs it was noted that if b was set as greater than the minimum weight in the weights distribution then little-to-no ''boxing'' occurred. For this reason b 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 Sect. 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 b, 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 (Stützle 2006). We use his ILS configured as follows: firstimprovement 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) (Taillard 1991) 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: 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ðROTSÞ 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 Sect. 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.  (2017))]; 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; 7. the LON diameter; 8. the number of compressed local optima (after connected LON nodes of the same fitness are compressed together-labelled as comp.opt in the Figures); 9. the correlation between the fitnesses of neighbours in the LON (fit.fit.corr); 10. the number of sink nodes present (sinks); 11. and the sub-optimal sink strength (that is, the total incoming edge weight to any sub-optimal sink nodes-so.strength).

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.

Distribution analysis
In Figs. 3 and 4, 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 are considered in these plots. We 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 Fig. 3a, b 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 Fig 3a, b, 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 Sect. 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 (Fig. 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 Fig. 4a, b. The range of fractal dimensions for a LON is a way to quantify the extent of multifractality present and is calculated as maximumvalueminimumvalue 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 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 lowestthe 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 Fig. 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 (Ochoa and Veerapen 2016). We begin with visualisation before moving onto correlation analysis (Sect. 5.3) and machine learning models (Sect. 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.
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 Fig. 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- Fig. 5a shows that the ''had20'' LON has many, and they appear to be densely connected to other nodes. Contrarily, the ''lipa20b'' LON in Fig. 5b has a single global optimum, which seems to be more sparsely connected within its network. Also noteworthy is the relative sizes of the nonoptimal (grey) nodes. In Fig. 5a there are many large nodes which are sub-optimal 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 (Fig. 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 Figs. 6a, b.
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 Fig. 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 Fig. 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 and 8 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ðROTSÞ on the QAP instances, alongside the proposed fractal dimension features (unique FDs, range FD, median FD, max FD and min FD)-in Fig. 7, these concern deterministic dimensions; in Fig. 8, they are probabilistic. Also shown are LON features which are not associated with fractal complexity-these were introduced in Sect. 4.5. 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 Fig. 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 Fig. 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 the approach of previous work on multifractality in LONs (Thomson et al. 2018a), although the difference is not pronounced. In the pðROTSÞ column, there are only weak positive correlations with the fractal dimension variables, and indeed there are no strong correlations with any of the fitness landscape variables included.
Next we will consider the correlation plot which includes probabilistic fractal dimension variables in Fig. 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 Fractal geometry in LONs 327 0.001. These two also have a weak negative correlation with pðROTSÞ-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ðROTSÞ 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).

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 regression produces a stronger model fit. This is likely because random forest trees are adept at considering nonlinearities 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ðROTSÞ 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ðROTSÞ model. Now let us look at the random forest predictor importance rankings, which are provided in Fig. 9. The values are averaged over 10,000 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, Fig. 7 Spearman correlations between pairs of variables including pðILSÞ and pðROTSÞ, fractal dimension metrics for the LONs, and other landscape features 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 Fig. 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, 9). 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ðROTSÞ  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 Fig. 10, we draw your attention to the two fractal dimension metrics in the iterated local search plot (Fig. 10a): they are among the most important predictors, ranking third and fourth. Their values are between 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 Fig. 10b. Instead, the strongest predictors for pðROTSÞ 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 (Ochoa and Herrmann 2018) was used to build the LONs.
(a) (b) Fig. 9 Variable importance values for random forest models; the models include deterministic fractal dimension features as part of the predictor set 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 (Bo_ zejko et al. 2018). 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 Daolio et al. (2011)). 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.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/. Fig. 10 Variable importance values for random forest models; the models include probabilistic fractal dimension features as part of the predictor set Fractal geometry in LONs 331