Eric Angel
email: eric.angel@ibisc.univ-evry.fr

Romain Campigotto
email: romain.campigotto@lamsade.dauphine.fr

Christian Laforest
email: christian.laforest@isima.fr

Implementation and Comparison of Heuristics for the Vertex Cover Problem on Huge Graphs ⋆

Keywords: implementation of algorithms, experimental analysis, huge graphs, low memory, vertex cover

Implementation and Comparison of Heuristics for the Vertex Cover Problem on Huge Graphs

Eric Angel, Romain Campigotto, Christian

Introduction

The vertex cover problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF] is a well-known classical NP-complete optimization graph problem that has received a particular attention these last decades. In particular, it occurs in many concrete situations [START_REF] Pirzada | Applications of Graph Theory[END_REF], in fields such as biology, meteorology, finance, etc. where amount of data is more and more important. This leads to the problem of designing algorithms well suited to cope with such large instances.

Notations. Graphs G = (V, E) considered throughout this paper are undirected, simple, unweighted and represent the instance to be treated here. We denote by n the number of vertices (n = |V |) and by m the number of edges (m = |E|). For any vertex u ∈ V , we denote by N (u) = {v | uv ∈ E} the set of neighbors of u and we call degree the number of neighbors of vertex u.

Definition of the Vertex Cover Problem. A cover C of G is a subset of vertices such that every edge contains (or is covered by) at least one vertex of C, that is C ⊆ V and ∀e = uv ∈ E, one has u ∈ C or v ∈ C (or both). The vertex cover problem is to find a cover of minimum size. We denote by OP T the size of an optimal cover for a given graph.

Related Work. Several studies focused on massive data sets these last decades [START_REF]Handbook of Massive Data Sets[END_REF]. In particular, for the max clique problem [START_REF] Bomze | Handbook of Combinatorial Optimization, chap. The Maximum Clique Problem[END_REF], experiments on graphs with 53 • 10 6 vertices and 170 • 10 6 edges have been performed [START_REF] Abello | External Memory Algorithms[END_REF]. But no such study has been done for the vertex cover problem. However, it has been extensively studied theoretically: many exact (exponential) algorithms, approximation algorithms, online algorithms, etc. have been proposed (due to space limitations, we do not give references about these works: some of them can be found in the introductions of [START_REF] Bar-Yehuda | Minimum Vertex Cover in Rectangle Graphs[END_REF] and [START_REF] Escoffier | Complexity and Approximation Results for the Connected Vertex Cover Problem in Graphs and Hypergraphs[END_REF]). Several experimental studies have already been made, often to compare the quality of several algorithms [START_REF] Delbot | Analytical and Experimental Comparison of Six Algorithms for the Vertex Cover[END_REF][START_REF] Gilmour | Kernelization as Heuristic Structure for the Vertex Cover Problem[END_REF] or validate specific methods [START_REF] Alber | Experimental Evaluation of a Tree Decomposition-Based Algorithm for Vertex Cover on Planar Graphs[END_REF][START_REF] Asgeirsson | Vertex Cover Approximation on Random Graphs[END_REF]. Nevertheless, no one achieved the huge graph sizes we consider in this paper (in these studies, the largest graphs has 10,000 vertices).

Our General Model of Treatment. To the intrinsic NP-completeness is added the difficulty to manipulate graphs and run algorithms with severe constraints. Indeed, with respect to their huge sizes, the processing unit (we consider a standard computer) cannot load them entirely in its memory. Moreover, the graph, which is stored on an external disk, must not be modified, since it often comes from experimentations and can be used by different users for different goals. Specifically, the important cost of graph creation forces us to preserve its integrity, in order to be able to run several algorithms on it.

Organization of the Paper. We give in the next section a general description of our experiments. We present, analyze and compare in Sect. 3 results obtained by executing the six algorithms on several instances. Finally, in Sect. 4, we conclude and give some perspectives.

General Description

In this section, we describe elements used for experiments. Programs (executables, with source code) are available at [1].

Technical Characteristics. The "Processing Unit" is a laptop computer with a Dual Core processor running at 2.8 GHz, 6 Mb cache memory and 4 Gb RAM. Graphs and Covers are stored on the same external hard drive, which is a USB 2.0 hard disk of 2 Tb, running at 7200 revolutions/minute and equipped with 8 Mb cache memory. Programs are written in language C, C99 standard, in order to use specific data type unsigned long long and associated functions to read and write binary files.

Storage and Reading of Graphs.

There exist many ways to store a graph: with an adjacency matrix, an adjacency list, etc. We use the method described in [START_REF] Ajwani | Design, Implementation and Experimental Study of External Memory BFS Algorithms[END_REF] (for more details, report to the Sect. 3.2 page 20). More precisely, our graphs are stored with two files:

.list file which contains 2m + 1 values: the number of vertices in the graph (which is needed by several algorithms to create an n bits array) and the list of the neighbors of vertices in the graph; .deg file which contains n+1 values, which are needed to access to the neighbors of a vertex and compute its degree.

The n vertices are labeled from 0 to n -1. The .list file contains first the value n, then vertices of set N (0), then vertices of set N (1), etc. (however, neighbors of each vertex can be stored in any order, not necessarily following the order of labels). The .deg file contains, for each vertex (and by increasing order of labels), the place of its first neighbor (in the .list file). It contains n + 1 values, in order to compute the degree of the last vertex (the last value of the .deg file points to the end of the .list file). Indeed, to compute degree of vertex i, we subtract the i th value from the (i + 1) th . Figure 1 shows an example on a small graph. Algorithms scan graphs by reading the two files described previously. First, the .list file is read to know the number of vertices. Then, the .deg file is read to know the place of the first neighbor of the first vertex and the place of the first neighbor of the second vertex. Thereby, there are two pointers (the first one is followed by the second one) which delimit the set of neighbors of a vertex. Once these places are known, the reading process continues in the .list file, where neighbors of the first vertex are retrieved. If the treatment unit does not need to get all of them, it can "step over" the remaining neighbors and go immediatly to the neighbors of the second vertex. It proceeds in the same way for the following vertices.

If algorithms need to know degrees of neighbors, they read independently the .deg file with another playhead, which can be moved at a precise place to get the two successive values needed to compute the degree.

All the algorithms read the .list file in a sequential way, but they can "step over" values that are not needed. Algorithms which do not need to compute degrees of neighbors read the .deg file also in a sequential way.

Algorithms Implemented. We have implemented six algorithms adapted to the treatment of huge graphs: LR, ED, S-Pitt, LL, SLL and ASLL.

LR has been proposed in [START_REF] Delbot | A Better List Heuristic for Vertex Cover[END_REF]. ED is the 2-approximative algorithm which returns vertices of a maximal matching. S-Pitt is a probabilistic algorithm inspired by the algorithm presented in [START_REF] Pitt | A Simple Probabilistic Approximation Algorithm for Vertex Cover[END_REF]: it has an expected approximation ratio equal to 2. The authors have done a theoretical study of LL, SLL and ASLL in [START_REF] Angel | Analysis and Comparison of Three Algorithms for the Vertex Cover Problem on Large Graphs with Low Memory Capacities[END_REF].

We present now a basic description of the six algorithms, by giving conditions to put vertices of the input graph into the solution.

Let G = (V, E) be a graph. Let C be the cover under construction. For each vertex u ∈ V , we have

LR: if u ∈ C, {v | uv ∈ E ∧ v ∈ C} is put in C; ED: if u ∈ C and if u has a neighbor v ∈ C, u and v are put in C; S-Pitt: if u ∈ C and if u has a neighbor v ∈ C, either u or v is put in C with
equiprobability assumption;

LL: u is put in C if it has at least one neighbor v such that v > u (their labels are compared); SLL: u is put in C if it has at least one neighbor v such that d(v) < d(u) or d(v) = d(u) and v > u; ASLL: u is put in C if it has at least one neighbor v such that d(v) > d(u) or d(v) = d(u) and v < u.
Now, we can describe how we have implemented these algorithms, in relation with the way that the graphs are stored on the external hard disk.

As described above, the algorithms scan graphs vertex by vertex and, for each current vertex u, scan its neighbors one by one. If an algorithm decides that u belongs to the solution (applying the conditions given in the descriptions of the algorithms above), u is put immediately and definitively into the cover. Then, the algorithm steps over its remaining neighbors and goes to the next vertex. Otherwise, it gets the next neighbor of u (and, at the end, requires the next vertex like in the previous case). Also, when an algorithm scans a vertex u which is already in the cover, it goes immediately to the next vertex, without scanning its neighbors.

It is worth noticing that LR, ED and S-Pitt need to allocate an n bits array to mark vertices sent to the solution (reading on the external hard drive during the execution would take too long); SLL and ASLL need to compute degrees of neighbors.

Writing the Covers on the Disk. A cover is written as a list of vertex labels into a file, which is built piece by piece: once an algorithm decides to put a vertex into the solution, it writes it into the cover file. A vertex cannot appear twice, because algorithms have been designed to produce no duplicates.

Example of Execution of Algorithm LR. We consider the graph given in Fig. 2. The execution of LR on it works as follows. At the beginning, the cover C (which is materialized by an 6 bits array in the internal memory of the computer) is empty (i.e. all the binary flags are lowered).

1. C = ∅. We consider the vertex 0. We get its neighbor, 3: we put it in C (it is written on the disk and the corresponding flag in the memory is raised). Then, we get its neighbor, 1: we put it in C. 2. C = {1, 3}. The vertex 1 is not treated since it is already in C. 3. C = {1, 3}. We treat the vertex 2. We get its neighbor, 5: it is put in C.

Then, we get its neighbor, 1: nothing is done since it is already in C. 4. C = {1, 3, 5}. The vertex 3 is not treated since it is already in C. 5. C = {1, 3, 5}. We treat the vertex 4. Its neighbors (1, 3 and 5) are retrieved but not considered, because they are already in C. 6. C = {1, 3, 5}. The vertex 5 is not treated since it is already in C.

Hence, we have just scanned seven vertices in the .list file (which contains sixteen labels) and the cover produced by LR contains three vertices: 1, 3 and 5.

Graph Families Used. We have executed our algorithms on different graphs:

on sparse graphs (where m ∈ O(n)): ButterFly graphs [START_REF] Weisstein | Butterfly graph, from MathWorld -A Wolfram Web Ressource[END_REF], de Bruijn graphs [START_REF] De Bruijn | A Combinatorial Problem[END_REF] and grid graphs; -on dense graphs (where m ∈ Θ(n 2 )): hypercubes, complete bipartite graphs and complete split graphs 4 .

We have chosen these graphs because they can be easily generated (we can produce huge size graphs with a standard computer). In that sense, we have designed generators to construct these specific graphs. For example, CPU running times for the instances generation of size 100 • 10 9 are between 5 hours than 7 hours. Also, the size of their optimal solutions is known (we can give exact sizes, excepted for de Bruijn graphs where we can only give lower bounds). Thus, it is possible to present results on the quality of algorithms.

We have also chosen to execute algorithms on random power law graphs, where degree sequences follow a power law. We have used generator described in [START_REF] Vigier | Random Generation of Large Connected Simple Graphs with Prescribed Degree Distribution[END_REF], which is based on the Molloy and Reed model [START_REF] Molloy | A Critical Point for Random Graphs With a Given Degree Sequence[END_REF].

This generator is able to produce random power law graphs, but it cannot create them in an online way: a memory space linear to the graph size is needed. However, on our computer, we can still create graphs with more than 10 • 10 6 vertices and edges, with CPU running times less than 3 hours.

Graph Sizes. The size of a graph is given by its number of vertices and edges. Hence, we denote by graph size the value n + m.

The huge size notion is relative: it depends on context considered (e.g. the size limits for algorithms with exponential complexity are lower than for algorithms with linear complexity). So, we have defined several levels for graph sizes.

1 st level. The graphs size is about 200 • 10 6 (several Gb on disk). This is the largest random power law graphs size that can be generated on our computer. 2 nd level. The graphs size is about 30 • 10 9 (more than 100 Gb on disk). The algorithms SLL and ASLL begin to reach their limits in terms of running times on our computer. 3 rd level. The graphs size is about 100 • 10 9 (around 1.5 Tb on disk). With our computer, we cannot allocate an n bits array if the number of vertices is bigger than 30 • 10 9 .

Experimentations. For the first level, we have executed algorithms S-Pitt, LL, SLL and ASLL five times on each instance (algorithms LR and ED are deterministics). From second level, we simulated a user having limited resources (time and disk space). So, for each graph, we have executed the six algorithms once.

Based on results obtained and resources already spent, we have executed again several algorithms (often one time).

A total of thirty-four executions was made: thirteen in the first level, six in the second level and two in the third level. We have also executed our algorithms on thirteen instances for which sizes are between 10 6 and 4 • 10 6 . But, due to space limitations, we do not give details on results obtained for these instances. We can however indicate that they are broadly similar to the results obtained in the first level.

Results and Observations

Evaluated Criteria. We have focused on quality of solutions produced by algorithms and complexities, expressed by the number of requests made to the instance, i.e. the number of neighbors read in the .list file. We have also considered running times. For that, we have used the UNIX command /usr/bin/time, which gives the time used by the processor during a program execution.

Results

Presentation. We give one table per criterion. For algorithms that have been executed more than once, we give values corresponding to the best solution (in terms of quality).

For each instance, the best value (among the set of values presented for the six algorithms) is in bold font. Conversely, when values are bad (e.g. an algorithm which returns almost all the vertices or performs more than m requests), numbers are in italics font.

Table 1 (resp. 2) gives quality of solutions (resp. complexity in number of requests) obtained in the first level on graphs created by our generators and on random graphs. For random power law graphs, the last digit given in the instance name (starting by rg) denotes the minimum degree of the graph. For the first level, expected running times of each algorithm are between twenty seconds and two minutes.

Tables 3 and4 give respectively quality of solutions and number of requests obtained for the second and third levels. Table 5 gives CPU running times. For the third level, we only give values for algorithms that have been executed until the end. Observations on Quality of Solutions. As we can see on Tables. 1 and 3, the algorithm LR is almost always the best. Moreover, it often returns the optimal solution. However, it can be very bad on complete bipartite and split graphs.

In general, SLL offers good performance, especially on random power law graphs (its performance is close to LR). Nevertheless, it is less efficient on regular graphs 5 .

The global performance of algorithms S-Pitt and LL is intermediate but, for LL, it fluctuates more than S-Pitt. Indeed, on one instance, LL can be the best or the worst, that is not the case for S-Pitt.

Finally, ED and ASLL are overall the worst algorithms (and ED reaches often its approximation ratio of 2). For ED, these results confirm observations made by F. Delbot et al. [START_REF] Delbot | Analytical and Experimental Comparison of Six Algorithms for the Vertex Cover[END_REF].

Observations on the Number of Requests. As we can see on Tables. 2 and 4, the algorithm ED is almost always the best. Furthermore, it always performs less than m requests.

The algorithm LR often reaches m requests (it cannot perform worse), except on instances on which it returns a bad solution (it is better on them).

The performance of S-Pitt is close to LR: it is often better on specific graphs (except on complete split graphs) but it is worse on random graphs.

Algorithms LL, SLL and ASLL can perform more than m requests. This explains the fact that LR is generally the second algorithm in terms of complexity (even if its upper bound of m requests is often reached). ASLL is better, especially on random power law graphs and complete split graphs.

Analysis of Running Times. We focus on values presented for the second level in Tab. 5 (CPU running times obtained in first level are too similar to be exploited). To obtain an estimation of real running times (observed on our computer), one can multiply by 3.2 (resp. 1.6) CPU running times given in Tab. 5 for algorithms LR, ED, S-Pitt and LL (resp. SLL and ASLL).

Algorithms SLL and ASLL are different because they have to use another playhead on .deg file to calculate degrees of neighbors. Therefore, their CPU running times are bigger. For this reason, we focus primarily on values observed for algorithms LR, ED, S-Pitt and LL.

On sparse graphs (where n and m are similar), CPU running times are close. They depend on number of requests performed by algorithms and size of covers constructed. Indeed, the size of solutions can be as huge as n, and writing on a disk is longer than reading. Moreover, there is often a trade off between the number of requests performed and the quality of solutions constructed: algorithms which produce the best solutions often perform the biggest number of requests (in any case, on one instance, an algorithm is never both the best in terms of quality of solution and complexity).

On dense graphs (where n is negligible compared to m), the analysis is less intricate because the size of covers written is tiny compared to the number of requets performed. Thus, CPU running times are mainly influenced by the number of requests done.

But these two criteria are not sufficient to explain CPU running times observed. Another technical aspects, linked to operating systems, are involved. Indeed, the access to the hard drive is indirect: the processing unit uses buffers. Also, the atomic unit of access depends on the size of disk sectors. If we read only one vertex on .list file, the system loads more vertices into its buffers. Hence, the number of physical access is lower than the number of requests performed. One overtakes in this regard practical considerations highlighted in the I/O-efficient model (see [START_REF] Vitter | Algorithms and Data Structures for External Memory[END_REF] for a survey).

Limits Encountered on our Machine. In the third level, we have generated two instances: a complete bipartite graph with 630,000 vertices and a ButterFly graph of dimension 30. On the complete bipartite graph, we were able to run algorithms LR, ED, S-Pitt and LL: their real running times do not exceed (on our computer) eight hours (executions of SLL and ASLL were stopped after twenty hours). On the ButterFly graph, we can only use LL (its real running time is about fifteen hours) because, on our computer, we cannot allocate an array of 33 • 10 9 bits (and executions of SLL and ASLL would take too long).

Therefore, LL is the only algorithm that can be run with our computer on all instances.

General Synthesis

We have implemented six algorithms for the vertex cover problem on huge graphs. We were able to run these algorithms with a standard laptop computer on instances of sizes up to 30•10 9 vertices and edges (about 300 Gb on disk). The CPU running times we obtained do not exceed eight hours (and corresponding real running times are lower than ten hours).

We have observed that SLL and ASLL are almost always the slowest, since they have to compute degrees of neighbors. In this direction, they are "less adapted". However, SLL is still interesting, because it can give good solutions.

To test limits of algorithms, we have generated two instances of sizes about 100 • 10 9 vertices and edges (at least 1 Tb on disk): a complete bipartite graph with 630,000 vertices and 95 • 10 9 edges (a dense graph), and a ButterFly graph of dimension 30, with 33 • 10 9 vertices and 64 • 10 9 edges (a sparse graph).

-On the complete bipartite graph, we were able to run LR, ED, S-Pitt and LL (executions of SLL and ASLL were stopped before the end). For the slowest (LR), its CPU running time barely reaches two hours (and its corresponding real running time is about seven hours). -On the ButterFly graph, we were only able to execute the algorithm LL:

executions of LR, ED and S-Pitt failed because we could not allocate an array of 33 • 10 9 bits (and, as for the complete bipartite graph above, SLL and ASLL were stopped before the end).

General Observations. By summarizing the set of results presented for instances we used, among the six algorithms, LR is the one which gives the best solutions. It is closely followed by SLL, while ED and ASLL gives the worst solutions. However, ED performs the smallest number of requests. Based on that, choosing an algorithm which satisfies both quality of solutions and complexity in number of requests is difficult. On sparse graphs, where running times are often similar, we should promote the quality of solutions. Therefore, LR is a good candidate. Unfortunately, it needs to allocate an n bits array to be run, that is not always possible. On dense graphs, running times can change significantly, this makes choice trickier, since the most efficient algorithms are often the slowest.

Perspectives. We could extend our work by designing efficient algorithms on large instances for other problems. Then, we could compare our treatment method with the existing ones, e.g. with the semi-external greedy randomized adaptive search procedure presented in [START_REF] Abello | External Memory Algorithms[END_REF] for the max clique problem.

15 Fig. 1 .

 151 Fig. 1. Storage of a graph with 5 vertices and 7 edges

Fig. 2 .

 2 Fig. 2. Graph with 6 vertices and 8 edges

Table 1 .

 1 Quality of solutions obtained in the first level, expressed in percentage of n (sizes of optimal covers for random power law graphs cannot be estimated)

	Instance	n	OP T	LR	ED S-Pitt LL	SLL ASLL
	butterfly-21	46,137,344 50	50	100 78.16 90.91 85.06 90.91
	debruijn-25	33,554,432 > 50 66.67 88.89 77.56 66.78 72.71 82.39
	grid-6000.9000 54,000,000 50	50 99.99 81.41 99.99 99.99 99.99
	hypercube-23	8,388,608 50	50 99.97 99.26 99.99 99.99 99.99
	compbip-7000.15000 22,000 31.82 68.18 63.64 62.98 31.82 31.82 68.18
	split-7500.12000	19,500 38.46 99.99 61.67 61.39 38.46 38.46 99.99
	rg-20m 1	20,000,000 -	9.94 19.65 13.41 49.30 10.42 99.99
	rg-20m 2	20,000,000 -	34.12 62.89 45.60 49.30 36.89 99.88
	rg-25m 1	25,000,000 -	14.19 28.09 18.91 41 14.89 99.95
	rg-25m 2	25,000,000 -	38.76 69.99 51.29 48.13 42.12 99.65
	rg-30m 1	30,000,000 -	43.48 76.61 57.17 59.02 47.44 97.57
	rg-30m 2	30,000,000 -	15.68 31.05 21.10 30.98 16.46 99.93
	rg-35m 2	35,000,000 -	43.12 76.16 56.79 53.70 46.98 99.12

Table 2 .

 2 Number of requests performed in the first level (in percentage of m)

	Instance	m	LR	ED S-Pitt LL	SLL ASLL
	butterfly-21	88,080,384 100 60.51 91.94 103.30 104.76 80.16
	debruijn-25	67,108,861 66.67 61.11 87.03 113.87 106.84 100.48
	grid-6000.9000 107,985,000 100 49.93 86.24 91.66 91.66 75.03
	hypercube-23	96,468,992 100 11.83 23.57 17.54 17.58 17.51
	compbip-7000.15000 105,000,000 100 53.34 54.31 100.01 100.01 100.02
	split-7500.12000 118,121,250 0.02 47.47 47.82 76.20 76.20 0.17
	rg-20m 1	59,624,494 49.12 47.99 55.20 34.29 59.22 57.73
	rg-20m 2	90,808,193 40.23 36.92 48.73 38.20 56.57 34.27
	rg-25m 1	70,911,180 45.44 44.57 53.13 36.23 58.19 47.08
	rg-25m 2	87,837,432	45 40.93 55.45 51.81 65.50 41.84
	rg-30m 1	82,356,722 50.09 45.04 62.78 57.89 75.78 52.29
	rg-30m 2	81,819,916 44.86 44.02 52.84 39.76 58.28 45.23
	rg-35m 2	96,555,269 50.14 45.11 62.70 62.52 75.62 53.31

Table 3 .

 3 Quality of solutions obtained in the second and third levels, expressed in percentage of n

	Instance	n	OP T	LR	ED S-Pitt LL SLL ASLL
	butterfly-28	7,784,628,224 48.28 48.28 96.55 78.76 93.24 93.10 96.55
	debruijn-33	8,589,934,592 > 50 66.67 88.89 77.56 96.55 99.99 99.99
	grid-75000.90000	6,750,000,000 50	50 99.99 81.41 99.99 99.99 99.99
	hypercube-30	1,073,741,824 50	50 99.99 99.78 99.99 99.99 99.99
	compbip-35000.500000	535,000	6.54 93.46 13.08 13.11 6.54 6.54 93.46
	split-70000.180000	250,000	28 99.99 48.19 48.01 28	28 99.99
	butterfly-30	33,285,996,544 48.28	-	-	-98.26 -	-
	compbip-250000.380000	630,000	39.68 60.32 79.37 79.46 84.39 -	-

Table 4 .

 4 Number of requests performed in the second and third levels, expressed in percentage of m

	Instance	m	LR	ED S-Pitt LL	SLL ASLL
	butterfly-28	15,032,385,540 99.99 61.91 91.59 103.50 101.79 72.62
	debruijn-33	17,179,869,183 66.67 61.11 87.03 108.91 108.33 91.67
	grid-75000.90000	13,499,835,000 100 49.98 86.24 91.67 91.67	75
	hypercube-30	16,106,127,360 100 9.10 18.22 13.42 13.33 13.33
	compbip-35000.500000 17,500,000,000 100	93 92.97 100.01 100.01 100.02
	split-70000.180000 15,049,965,000 0.002 60.24 60.45 83.72 83.72 0.02
	butterfly-30	64,424,509,440	-	-	-102.18	-	-
	compbip-250000.380000 95,000,000,000 100 34.21 34.05 25.88	-	-

Table 5 .

 5 CPU running times obtained in the second and third levels (number of executions are given in parenthesis: we give the average time here)

	Instance	LR	ED	S-Pitt	LL	SLL	ASLL
	butterfly-28	1:11:55 1:15:53 1:18:41 (2) 1:20:35 (2) 5:10:47	3:38:43
	debruijn-33	1:20:37 1:24:10 1:26:14 (2) 1:29:35 (2) 7:43:53	5:08:37
	grid-75000.90000	1:02:47 1:08:57 1:09:58 (2) 1:11:35 (2) 3:24:52	3:01:32
	hypercube-30	0:41:06 0:33:21 0:36:02 (3) 0:33:19 (3) 1:07:46 (2) 1:07:38 (2)
	compbip-35000.500000 0:23:15 0:22:23 0:22:28 (4) 0:22:13 (4) 6:11:51	6:12:03
	split-70000.180000 0:00:17 0:15:21 0:15:36 (5) 0:16:11 (5) 4:27:39 0:00:29 (8)
	butterfly-30	-	-	-	5:47:43	-	-
	compbip-250000.380000 2:02:16 1:01:19 1:01:21	0:32:17	-	-

A complete split graph is a complete bipartite graph in which the vertices subset of lowest size is changed to a clique.

A regular graph is a graph where all the vertices have the same degree.

Acknowledgements

We would like to thank the anonymous referees for their insightful comments and suggestions, which have helped to improve the presentation of this paper.

⋆ Work partially supported by the French Agency for Research under the DEFIS program TODO, ANR-09-EMER-010.


Bernard Ng
email: bernardyng@gmail.com

Gael Varoquaux

Jean Baptiste Poline

Bertrand Thirion

A Novel Sparse Graphical Approach for Multimodal Brain Connectivity Inference

Keywords: brain connectivity, diffusion MRI, fMRI, multimodal integration

Despite the clear potential benefits of combining fMRI and diffusion MRI in learning the neural pathways that underlie brain functions, little methodological progress has been made in this direction. In this paper, we propose a novel multimodal integration approach based on sparse Gaussian graphical model for estimating brain connectivity. Casting functional connectivity estimation as a sparse inverse covariance learning problem, we adapt the level of sparse penalization on each connection based on its anatomical capacity for functional interactions. Functional connections with little anatomical support are thus more heavily penalized. For validation, we showed on real data collected from a cohort of 60 subjects that additionally modeling anatomical capacity significantly increases subject consistency in the detected connection patterns. Moreover, we demonstrated that incorporating a connectivity prior learned with our multimodal connectivity estimation approach improves activation detection.

Introduction

Recent evidence suggests that the effects of many neurological diseases are manifested through abnormal changes in brain connectivity [START_REF] Fox | Clincial Applications of Resting State Functional Connectivity[END_REF]. Techniques for inferring connectivity from functional magnetic resonance imaging (fMRI) data can be largely divided into two categories. One category comprises techniques, such as the seedbased approach and independent component analysis (ICA) [START_REF] Fox | Clincial Applications of Resting State Functional Connectivity[END_REF], which group brain areas into networks. The other category includes techniques that estimate the connectivity between brain areas [START_REF] Varoquaux | Brain Covariance Selection: Better Individual Functional Connectivity Models Using Population Prior[END_REF] and apply graph-theoretic measures to characterize the estimated connection structure. We focus on connectivity estimation in this work.

The strong noise inherent in fMRI data and its high dimensionality given the typically small sample sizes pose major challenges to reliable connectivity estimation [START_REF] Varoquaux | Brain Covariance Selection: Better Individual Functional Connectivity Models Using Population Prior[END_REF]. Since neural dynamics is largely shaped by the structure of the underlying fiber pathways [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF], informing connectivity estimation with anatomical information extracted from e.g. diffusion MRI (dMRI) data should prove beneficial. Past studies that jointly examined dMRI and fMRI data have primarily focused on comparing connectivity measures estimated from each modality separately [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF]. The general finding is that brain areas with high anatomical connectivity typically exhibit high functional con-nectivity [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF], but the converse does not necessarily hold due to factors, such as noise-induced correlations in fMRI observations, indirect functional connections, and tractography errors [START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF]. Other widely-used approaches for multimodal analysis include employing fMRI to guide seed selection in tractography [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF]. The use of fiber bundle shapes to predict activated brain areas has also been explored [START_REF] Zhu | Optimization of Functional Brain ROIs via Maximization of Consistency of Structural Connectivity Profiles[END_REF]. A major limitation to the aforementioned approaches is that pooling information extracted independently from each modality does not capitalize on the complementary information that a joint analysis of the two modalities would facilitate. To exploit the joint information in dMRI and fMRI data, variants of ICA and canonical correlation analysis [START_REF] Si | Discriminating Schizophrenia and Bipolar Disorder by Fusing fMRI and DTI in a Multimodal CCA + joint ICA Model[END_REF] have been proposed for identifying brain areas that display high correlations between anatomical and functional attributes, such as fractional anisotropy and activation effects. Recently, a probabilistic model has been put forth for combining dMRI and fMRI data in detecting group differences in brain connection structure [START_REF] Venkataraman | Joint Modeling of Anatomical and Functional Connectivity for Population Studies[END_REF].

In this paper, we propose a novel multimodal integration approach based on sparse Gaussian graphical model (SGGM) for estimating intra-subject brain connectivity. Specifically, we cast connectivity estimation as a sparse inverse covariance learning problem [START_REF] Hsieh | Sparse Invers Covariance Matrix Estimation Using Quadratic Approximation[END_REF]. Since elements of the inverse covariance matrix are proportional to the partial correlations between the associated variable pairs, zero entry would indicate conditional independence [START_REF] Hsieh | Sparse Invers Covariance Matrix Estimation Using Quadratic Approximation[END_REF]. Using SGGM thus enables simultaneous estimation of connection strength and structure. To integrate anatomical information into functional connectivity estimation, we adapt the degree of sparse penalization on each functional link based on its anatomical capacity. Functional connections with less anatomical support are thus more heavily penalized, which helps reduce false detection of noiseinduced functional connections. Also, using partial correlations as a measure of functional connectivity reduces the effects of indirect interactions. Furthermore, although larger penalizations are exerted on functional connections with no anatomical support, if ample evidence from fMRI data suggest the presence of such links, these connections will not necessarily be assigned zero connectivity. Our approach hence provides some tolerances to the inconsistencies between dMRI and fMRI-derived connectivity measures that hinder integration of these modalities. On a large dataset of 60 subjects, we showed that applying our multimodal approach significantly increases subject consistency in the detected connection structure over analyzing fMRI data alone. Enhanced group activation detection was also obtained by incorporating a connectivity prior [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF] learned with our proposed approach, thus demonstrating the gain of integrating anatomical and functional information in brain connectivity estimation.

Methods

In this work, we focus on integrating dMRI with resting state (RS) fMRI for estimating brain connectivity. For this, we propose a SGGM-based approach (Section 2.1). Critical to the estimation is the choice on sparsity level, which we optimize using cross validation with a refined grid search strategy (Section 2.2). We validate our approach based on subject consistency and group activation detection (Section 2.3).

Sparse Gaussian Graphical Model

Let S be a d×d sample covariance matrix computed from data that are assumed to follow a centered multivariate Gaussian distribution. In the present context of brain connectivity estimation, S corresponds to correlations between the RS-fMRI observations of d brain areas of a given subject. To learn a well-conditioned sparse inverse covariance matrix, Λ ˆ, we minimize the negative log data likelihood over the space of positive definite (p.d.) matrices, Λ > 0, with an l 1 penalty imposed on Λ [START_REF] Hsieh | Sparse Invers Covariance Matrix Estimation Using Quadratic Approximation[END_REF]:

        d i d j ij ij tr 1 1 0 ) det( log ) ( min arg ˆΛ W Λ SΛ Λ Λ  , (1)
where λ governs the overall level of sparsity and W ij differentially weights the amount of sparse penalization on each connection based on its anatomical capacity:

  , 0 exp        j i j i ij ij  K W (2)
where K ij corresponds to some measures of anatomical capacity. We set K ij as the fiber count between brain areas i and j [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF] with the amount of penalization saturating to zero for K ij >> σ, modeling how the additional fibers may reflect redundant wiring. Other decreasing functions of K ij can also serve as W ij . We defer the selection of W ij for future work. σ is learned from data (Section 2.2). In accordance to past findings [START_REF] Honey | Can Structure Predict Function in the Human Brain[END_REF][START_REF] Damoiseaux | Greater than the Sum of its Parts: A Review of Studies Combining Structural Connectivity and Resting-state Functional Connectivity[END_REF], Λ ij associated with brain areas that have fewer connecting fibers are more strongly penalized. Note that we have explicitly set W ij to 0 for i = j, which has been theoretically proven and empirically shown to provide more accurate solutions of (1) [START_REF] Duchi | Projected Subgradient Methods for Learning Sparse Gaussions[END_REF]. To solve (1), we employ a recent second-order algorithm [START_REF] Hsieh | Sparse Invers Covariance Matrix Estimation Using Quadratic Approximation[END_REF] that facilitates efficient computation of Newton steps with iterates guaranteed to remain p.d. This algorithm provides substantially faster convergence rate than current gradientbased methods. We define convergence as having a duality gap, η, below 10 -5 [START_REF] Duchi | Projected Subgradient Methods for Learning Sparse Gaussions[END_REF]:

  . 1 1 d tr d i d j ij ij        Λ W SΛ   (3)

Parameter Selection

The estimated connection structure critically depends on λ. To learn the optimal λ and σ data-drivenly for each subject, we combine cross validation with a refined grid search strategy, as summarized in Algorithm 1. The algorithm requires the following inputs: Z = an t×d matrix containing RS-fMRI time courses of d brain areas, K = fiber count matrix, λ lb and λ ub = initial lower and upper bounds of the λ search range, σ lb and σ ub = lower and upper bounds of the σ search range, R = number of refinement levels, F = number of grid points for λ, G = number of grid points for σ, and C = number of cross validation folds. We set λ ub to max|S ij |, i ≠ j, which corresponds to the maximum λ beyond which ij Λ ˆis guaranteed to shrink to 0 [START_REF] Duchi | Projected Subgradient Methods for Learning Sparse Gaussions[END_REF]. As for λ lb , we empirically found that setting λ lb to λ ub /100 assigns non-zero value to >90% of the elements in Λ ˆin the case of uniform sparse penalization. We fix σ lb and σ ub to the 25 th and 75 th percentile of the fiber count distribution. R, F, G, and C are set to 3, 5, 5, and 3, respectively. Set λ grid to a log grid with F grid points between λ lb and λ ub in decreasing order 7:

for f = 1 to F, g = 1 to G, and c = 1 to C 10:

Estimate sample covariance, S train , with all folds of Z except the c th fold 11:

Estimate sample covariance, S test , with the c th fold of Z 12:

Solve (1) to find Λ train with λ = λ grid (f), σ = σ grid (g) 13:

Compute data likelihood, dl(f,g,c), of S test given Λ train , based on (1) 14:

without the sparse penalization term 17: end 18: Find λ opt and σ opt based on maximum of the average dl(f,g,c) over C folds 19: if λ opt == λ grid (1) 20:

Set λ lb to λ grid [START_REF] Varoquaux | Brain Covariance Selection: Better Individual Functional Connectivity Models Using Population Prior[END_REF], set λ ub to λ grid (1) 22: else if λ opt == λ grid (F) 23:

Set λ lb to λ grid (F)/10, set λ ub to λ grid (F-1) 25: else 26:

Find f opt corresponding to λ opt 27:

Set λ lb to λ grid (f opt -1), set λ lb to λ grid (f opt +1) 29: end 30: end

Validation

We base our validation on: 1) increased subject consistency in the support of Λ ˆand 2) increased group activation detection. The rationale behind the first criterion is that subjects within the same population should have similar brain connection structure [START_REF] Varoquaux | Brain Covariance Selection: Better Individual Functional Connectivity Models Using Population Prior[END_REF]. As for the second criterion, we have shown in a previous work [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF] that incorporating a RS-connectivity prior improves group activation detection. Thus, greater increase in group activation detection presumably implies that the corresponding connectivity estimates better reflect the underlying neural circuitry.

Subject Consistency.

For each pair of subjects in our dataset (Section 3), we compare the support of their Λ ˆ's using the Dice Similarity Coefficient (DSC), defined as TPR / (2TPR + FPR + FNR), with each subject alternately taken as the reference. TPR, FPR, and FNR are the true positive rate, false positive rate, and false negative rate, respectively. Both connections commonly present in each pair of subjects as well as those in one subject but not the other are thus accounted for in our chosen metric.

Group Activation Detection. Our previously proposed connectivity-informed activation model [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF] is summarized below:

                  AX Y V AX Y V V AX Y 1 1 1 1 2 1 exp 2 1 , ~2 T tr N n  , ( 4
)              AX V A X V XX XX V A 1 2 2 2 2 exp 2 , , 0 ~22 T T T T tr MN m d     , ( 5
)
where Y is a d×n matrix containing task fMRI time courses of d brain areas of a given subject. X is a m×n regressor matrix [START_REF] Friston | Statistical Parametric Maps in Functional Imaging: A General Linear Approach[END_REF]. A is a d×m activation effect matrix to be estimated. V 1 and V 2 are d×d covariance matrices of Y and A, respectively. XX T models the correlations between the m experimental conditions. MN(0,V 2 ,αXX T ) denotes the matrix normal distribution, which serves as a conjugate prior of (4) with α controlling the influence of this prior on A. We assume V 1 = I d×d as conventionally done. V 2 is set to the connectivity estimates learned with our proposed approach.

3 Materials fMRI data. Task fMRI data were collected from 60 healthy subjects at multiple imaging centers. Each subject performed 10 language, computation, and sensorimotor tasks similar to those in [START_REF] Pinel | Fast Reproducible Identification and Large-scale Databasing of Individual Functional Cognitive Networks[END_REF] over a period of ~5 min. RS-fMRI data of ~7 min were also collected. Data were acquired using 3T scanners from multiple manufacturers with TR = 2200 ms, TE = 30 ms, and flip angle = 75 o . Standard preprocessing, including slice timing correction, motion correction, temporal detrending, and spatial normalization, was performed on the task fMRI data using the SPM8 software. Similar preprocessing was performed on the RS-fMRI data except a band-pass filter with cutoff frequencies at 0.01 to 0.1 Hz was applied. White matter and cerebrospinal fluid confounds were regressed out from the gray matter voxel time courses. dMRI data. dMRI data were collected from the same 60 subjects. Acquisition sequence similar to [START_REF] Jones | Isotropic Resolution Diffusion Tensor Imaging with Whole Brain Acquisition in a Clinically Acceptable Time[END_REF] was used with TR = 15000 ms, TE = 104 ms, flip angle = 90 o , 32 gradient directions, and b-value = 1300 s/mm 2 . MedINRIA was employed for tensor estimation and fiber tractography [START_REF] Toussaint | MedINRIA: Medical Image Navigation and Research Tool by INRIA[END_REF]. We warped our brain parcel template (described below) onto each subject's B 0 volume to facilitate fiber count computation.

Brain Parcellation. We divided the brain into P parcels (P set to 500) to enable finer brain partitioning than facilitated by standard brain atlases (typically P < 150). This choice of P provides a balance between functional localization and robustness to subject variability in tractography [START_REF] Venkataraman | Joint Modeling of Anatomical and Functional Connectivity for Population Studies[END_REF]. Parcellation was performed by concatenating RS-fMRI time courses across subjects and applying Ward clustering [START_REF] Michel | A Supervised Clustering Approach for fMRI-based Inference of Brain States[END_REF]. Parcel time courses were then generated by averaging the voxel time courses within each parcel, and normalized by subtracting the mean and dividing by the standard deviation.

Results and Discussion

To explore the gain in learning connectivity jointly from dMRI and RS-fMRI data, we compared based on DSC the subject consistency in the support of ij Λ ˆas estimated using our SGGM-based approach with anatomical capacity-weighted sparse penalization (SGGM-A) versus uniform sparse penalization without anatomical information (SGGM-U). To establish a baseline for comparison, we permuted 200 times the columns and rows of ij Λ ˆlearned with SGGM-A to generate a null distribution of DSC.

The average DSC of both SGGM-A and SGGM-U (Fig. 1(a)) were significantly greater than that of the null distribution based on a Wilcoxon signed rank test (p-value < 0.05). Hence, the observed subject consistency was significantly above chance. Moreover, DSC of SGGM-A was significantly higher than that of SGGM-U (p-values < 0.05), thus demonstrating the benefits of multimodal integration for connectivity estimation. In addition, we compared the sensitivity in group activation detection using: 1) ordinary least square (OLS) [START_REF] Friston | Statistical Parametric Maps in Functional Imaging: A General Linear Approach[END_REF] without any connectivity prior, 2) ridge regression to control overfitting [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF], connectivity-informed activation model [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF] with connectivity prior learned by applying 3) oracle approximating shrinkage (OAS) on RS-fMRI data only [START_REF] Chen | Shrinkage Algorithms for MMSE Covariance Estimation[END_REF], 4) SGGM-U on RS-fMRI data only, and 5) SGGM-A on dMRI and RS-fMRI data jointly. We examined 21 contrasts between the 10 experimental conditions. To enforce strict control over FPR so that we can safely base our validation on increased group activation detection, max-t permutation test [START_REF] Nichols | Controlling the Familywise Error Rate in Functional Neuroimaging: a Comparative Review[END_REF] was used. Fig. 1(b) shows the percentage of parcels detected with significant activation averaged over contrasts. In agreement with [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF], adding a connectivity prior substantially improved activation detection over OLS and ridge regression. More importantly, SGGM-A was found to outperform both SGGM-U and OAS. To assess whether the enhanced detection was statistically significant, we used a permutation test. Specifically, for each permutation, we first randomly selected half of the parcels and exchanged the labels (active or non-active) assigned by SGGM-A and each of the contrasted methods in turn for a given p-value threshold. We then computed the difference in the average number of detected parcels, and performed this procedure 10,000 times to generate a null distribution. For all p-value thresholds in Fig. 2(b), the origi-nal difference in the number of detected parcels was found to be greater than the 95 th percentile of the corresponding null distribution. Hence, the detection improvement was statistically significant. Qualitatively, our approach additionally detected relevant areas adjacent to those found by considering functional information alone (Fig. 1(c)).

Conclusions

We proposed a novel SGGM-based approach for multimodal brain connectivity inference. We showed that integrating dMRI and RS-fMRI data significantly increases subject consistency in the learned connection structure compared to analyzing RS-fMRI data alone. Enhanced group activation detection was also demonstrated. Our results thus suggest that connectivity estimated by combining anatomical and functional information may better resemble the underlying neural pathways than solely relying on functional information. A particularly important byproduct of this work is that our multimodal connectivity estimation approach in combination with our connectivity-informed activation model [START_REF] Ng | Connectivity-informed fMRI Activation Detection[END_REF] provides a statistically-rigorous test bench for comparing different dMRI processing strategies. Quantitative evaluation, as opposed to qualitative assessment as often used in most dMRI studies, is thus facilitated.

Fig. 1 .

 1 Fig. 1. Real data results. (a) Subject consistency based on DSC. (b) Percentage of parcels with significant activation averaged over contrasts vs. p-value thresholds. (c) Parcels detected for auditory tasks. Red = detected by SGGM-A only. Purple = detected by SGGM-A and SGGM-U. Green = detected by SGGM-A, SGGM-U, and OAS. Violet = detected by all methods.

Acknowledgments

This work was supported by the ANR grant, BrainPedia ANR-10-JCJC 1408-01, and the Berkeley INRIA Stanford grant. The data were acquired within the IMAGEN project. Jean Baptiste Poline was partly funded by the IMAGEN project, which receives funding from the E.U. Community's FP6, LSHM-CT-2007-037286. This manuscript reflects only the authors' views and the Community is not liable for any use that may be made of the information contained therein.


Fabio Daolio
email: fabio.daolio@unil.ch

Sébastien Verel
email: verel@i3s.unice.fr

Gabriela Ochoa

Marco Tomassini
email: marco.tomassini@unil.ch

Local Optima Networks and the Performance of Iterated Local Search

Keywords: F.2.m [Analysis of Algorithms and Problem Complexity]: Miscellaneous, G.2.2 [Discrete Mathematics]: Graph Theory-Network problems, I.2.8 [Artificial Intelligence]: Problem Solving, Control Methods, and Search-Heuristic methods Algorithms, Measurement, Performance Combinatorial Fitness Landscape, Local Optima Network, Local Search Heuristics

Local Optima Networks (LONs) have been recently proposed as an alternative model of combinatorial fitness landscapes. The model compresses the information given by the whole search space into a smaller mathematical object that is the graph having as vertices the local optima and as edges the possible weighted transitions between them. A new set of metrics can be derived from this model that capture the distribution and connectivity of the local optima in the underlying configuration space.This paper departs from the descriptive analysis of local optima networks, and actively studies the correlation between network features and the performance of a local search heuristic. The N K family of landscapes and the Iterated Local Search metaheuristic are considered. With a statisticallysound approach based on multiple linear regression, it is shown that some LONs' features strongly influence and can even partly predict the performance of a heuristic search algorithm. This study validates the expressive power of LONs as a model of combinatorial fitness landscapes.

INTRODUCTION

One of the most conspicuous limitations of heuristic search methods in combinatorial optimization, is the ability to become trapped at a local optimum [START_REF] Glover | Future paths for integer programming and links to artificial intelligence[END_REF]. The number and distribution of local optima in a search space have, therefore, an important impact on the performance of heuristic search algorithms. A recently proposed model: Local Optima Networks [START_REF] Tomassini | Complex-network analysis of combinatorial spaces: The NK landscape case[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF], provides an intermediate level of description for combinatorial fitness landscapes. The model has a higher descriptive power than a single statistical metric or set of metrics; but it also compresses the search space into a more manageable mathematical object. Specifically, a graph having as vertices the optima configurations of the problem and as edges the possible weighted transitions between these optima. This network representation allows the application of new analytical tools and metrics to study combinatorial landscapes, namely, those of complex networks analysis (e.g. degree distribution, clustering coefficient, assortativity, and community structure, to name a few). In previous work, alternative definitions of edges have been studied, and some of these metrics have been computed on the network extracted for two combinatorial problems: the NK family of landscapes [START_REF] Tomassini | Complex-network analysis of combinatorial spaces: The NK landscape case[END_REF][START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF][START_REF] Ochoa | First-improvement vs. best-improvement local optima networks of nk landscapes[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF][START_REF] Verel | Local Optima Networks with Escape Edges[END_REF], and the Quadratic Assignment 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]. Those studies have been mainly descriptive, although distinctive correlations between some network features and previous knowledge about search difficulty in these landscapes had been found. A previous related work [START_REF] Pelikan | Nk landscapes, problem difficulty, and hybrid evolutionary algorithms[END_REF], explored the relationships between N K landscape features and the performance of a hybrid EA. The authors use standard landscape metrics, and conduct a study based mainly on scatter plots. They suggest that: "further work is necessary to gain better understanding of the escape rate on the actual problem difficulty". The present study, addresses exactly this point. The goal is to systematically explore correlations between local optima network features and the performance of a stochastic local search algorithm (Iterated Local Search) running on the underlying combinatorial optimization problem (in this study the N K family of landscapes). The 'escape rate' is a property related to the local optimum, and the LON could be considered as a new tool to better understand problem difficulty. The ultimate goal is to have predictive models of the performance of specific local search heuristics when solving a given combinatorial optimization problem. This paper proposes an initial predictive model of performance based on the most influential LON features.

METHODS

The local optima network model of combinatorial landscapes, and the iterated local search metaheuristic are considered in this study. The relevant definitions and experimental setup are given below.

Local Optima Network

A fitness landscape [START_REF] Reidys | Combinatorial landscapes[END_REF] is a triplet (S, V, f ) where S is a set of potential solutions i.e. the search space, V : S -→ 2 S is a function that assigns to every s ∈ S a set of neighbors V (s) i.e the neighborhood structure, and f : S -→ R is the evaluation of the corresponding solutions i.e. the fitness function.

The present study uses the well-known N K-landscapes [START_REF] Kauffman | The Origins of Order[END_REF] as a benchmark set. It is a problem-independent model for constructing combinatorial landscapes that are tunably rugged. In the model, N refers to the number of (binary) genes in the genotype, i.e. the string length, and K to their epistatic interaction, i.e. the number of other loci (chosen at random here) that influence the fitness contribution of a particular gene. Starting from this N -loci 2-allele additive model, by increasing the non-linearity K from 0 to N -1, the landscapes can be tuned from smooth to rugged. Hence, S is the search space of all N -bit binary strings, and its size is S = 2 N . The neighborhood is defined by the minimum possible move on it, which is the single bit-flip operation, and the neighborhood size is V (s) = N . The fitness function evaluates each genotype s as f (s) =1 N N i=1 fi(si, si1, . . . , siK ), where the values of loci contributions fi : {0, 1} K+1 → [0, 1] are drawn uniformly at random in [0, 1].

A local optimum (LO), which is taken to be a maximum here, is a solution s * such that ∀s ∈ V (s), f (s) ≤ f (s * ). All optima are determined through exhaustive search by recursively running the 1-bit-flip best-improvement hill-climber, as in Algorithm 1. Let us denote by h(s) the operator that associates to each solution s ∈ S, the solution obtained after applying that algorithm until convergence to a LO. Since S is of finite size and there is no neutrality in f (s) values, this produces a partition of the landscape in a finite number of basins of attraction, which we can denote by LO1, LO2, LO3 . . . , LOnv, the local optima.

Algorithm 1: Best-Improvement Hill-Climber Choose initial solution s ∈ S ;

repeat choose s ∈ V (s), such that f (s ) = max x∈V (s) f (x); if f (s) < f (s ) then s ← s ; until s is a Local Optimum;
The connections among them account for the chances of escaping from a LO and jumping into another one with a controlled move [START_REF] Verel | Local Optima Networks with Escape Edges[END_REF]. There exists a directed transition eij from LOi to LOj if it exists a solution s such that d(s, LOi) ≤ D and h(s) = LOj, where the distance d(si, sj) can be measured in "number of moves" (i.e Hamming distance in the bit-flip operator case). The distance-threshold D ∈ N can be chosen accordingly to the applied perturbation; in this work, it is set to D = 2. The weight wij of such a transition is then: wij = {s ∈ S | d(s, LOi) ≤ D and h(s) = LOj}, i.e. the number of paths at distance D starting at LOi and reaching the basin of LOj. This can be normalized by the number of solutions within reach w.r.t. the given distance threshold, i.e. {s ∈ S | d(s, LOi) ≤ D}.

The weighted and directed graph G = (V, E) having the set of vertices V = {LO1, . . . , LOnv} and the set of edges E = {eij|wij > 0}, is the Local Optima Network (LON) [START_REF] Verel | Local Optima Networks with Escape Edges[END_REF].

Iterated Local Search

Iterated local search is a relatively simple but successful algorithm. It operates by iteratively alternating between applying a move operator to the incumbent solution and restarting local search from the perturbed solution. This search principle has been rediscovered multiple times, within different research communities and with different names [START_REF] Baxter | Local optima avoidance in depot location[END_REF][START_REF] Martin | Large-step Markov chains for the TSP incorporating local search heuristics[END_REF]. The term iterated local search (ILS) was proposed in [START_REF] Lourenço | Iterated local search[END_REF]. Algorithm 2 outlines the procedure.

Algorithm 2: Iterated Local Search s0 ← GenerateInitialSolution; s * ← LocalSearch(s0); repeat s ← Perturbation(s * ); s * ← LocalSearch(s ); s * ← AcceptanceCriterion(s * , s * ); until termination condition met;
In the present study, the base LocalSearch heuristic is the same best-improvement hill-climber of Algorithm 1, which stops on a LO. This heuristic uses the single bit-flip move operator. Therefore, a 2-bit-flip mutation is chosen as a Perturbation operator. When a different LO is found after that, the search process accepts the move if its fitness is higher (we are assuming maximization). With these settings, ILS is performing a first-improvement hill-climbing in the configuration space of the LON with the escape-edges at distance D = 2, as defined in section 2.1.

The search terminates at the global optimum, which for benchmark problems is known a priori, or when reaching a pre-set limit of fitness evaluations F Emax.

Performance Evaluation

As the performance criterion, we selected the expected number of function evaluations to reach the global optimum (success) after independent restarts of the ILS algorithm (Algorithm 2) [START_REF] Auger | Performance evaluation of an advanced local search evolutionary algorithm[END_REF]. This measure accounts for both the success rate (ps ∈ (0, 1]) and the convergence speed. In theory, after (N -1) unsuccessful runs stopped at Tus-steps and the final successful one running for Tssteps, the total running time would be T = N -1 k=1 (Tus) k + Ts. Taking the expectation and considering that N follows a geometric distribution 1 with parameter ps, it gives:

E(T ) = 1 -ps ps E(Tus) + E(Ts)
where in the present case E(Tus) = F Emax, the ratio of successful to total runs is an estimator for ps, and E(Ts) can be estimated by the average running time of those successful runs.

The ILS variant detailed in Sec. 2.2 is essentially incomplete, i.e. there are soluble problem instances for which the success probability is < 1 even in the limit of an infinite running time [START_REF] Hoos | Stochastic local search: Foundations and applications[END_REF]. Given the chosen acceptance criterion, the search will eventually get stuck. Indeed, out of the test runs, 1 instance with K = 16 and 3 instances with K = 17 were not solved. This theoretical limitation could be overcome by performing as many random restart as to cover the whole search space, but such a solution is of limited practical interest for large problems. In the present study, the success performance has been estimated on all the instances that had been solved at least once.

The benchmark set consists of N K-landscapes with N = 18 and K ∈ {2, 4, 6, 8, 10, 12, 14, 16, 17}. Those are the largest possible parameter combinations for which we could afford the exhaustive extraction of the local optima networks. In order to minimize the influence of the random generation of landscapes, 30 independent problem instances are considered for each combination of N and K, which accounts for a total of 270 instances in the problem set. The function-evaluations limit is set to 1/5 of S, i.e. F Emax 5.2 • 10 4 , success rate ps and running time of successful runs Ts are estimated on 500 random restarts per instance.

RESULTS

Descriptive Statistics

Table 1 summarises the LON metrics. Results are grouped according to the K value of the corresponding N K-landscapes, and present averages and standard deviations over 30 independent realizations per group. The number of local optima nv, is a metric familiar to any description of a rugged landscape. The other metrics are particular to the complex-network perspective provided by the LON model [START_REF] Newman | The structure and function of complex networks[END_REF].

From left to right in the table, lo represents the average length of the shortest paths that reach the global optimum starting from any other local optimum. The cost associated to an edge ij of the LON graph, is dij = 1/wij. This measure can be interpreted as the expected number of random perturbations for escaping LOi and entering exactly the basin of LOj. lv gives the average path length for the whole graph, which accounts for the weighted network characteristic length. lo is intuitively more directly related to the search difficulty. Indeed, lo increases steadily with the landscape ruggedness K, whilst the trend with lv is less clear. As a possible explanation, the network growth in terms of nodes might be counteracted by a growth in nodes connectivity: the number of weighted outgoing transitions from a given LO, i.e. its out-degree in the LON, increases with K (cf. column zout in Tab. 1).

Column f nn measures the correlation between the fitness of a node and the weighted average of the fitness of its nearest neighbors [START_REF] Barthélemy | Characterization and modeling of weighted networks[END_REF]. This is relevant as the ILS acceptance criterion takes the fitness values of LO into account. With respect to this metric, N Klandscapes behave at the LO level as they do at the solution level: they become the more and more uncorrelated with K approaching N -1 (see column f nn). In general, it is expected that a high and positive f nn correlation would help the search process.

Column wii reports the average number of perturbations that remain in the same basin of attraction, which is a proxy for the basin size. The larger the basins, the more difficult it is to escape them. Table 1 shows that average wii value decreases with the landscape ruggedness.

Column cc, reports the clustering coefficient [START_REF] Newman | The structure and function of complex networks[END_REF], which measures the ratio of connected triples in the LON graph. In a social network, this coefficient measures the probability for one's friends to be friends among each other. In the LO networks, it provides an index of topological locality for the transitions between local optima. Table 1 suggest that cc values decrease steadily with increasing K, but so does the LON density [START_REF] Verel | Local Optima Networks with Escape Edges[END_REF] and that might be the reason.

Column zout counts the number of transitions departing from a given LO. It is relevant to know whether all the transitions have the same rate, or if there is a preferred direction. To this aim, the disparity score y2 gauges the weight heterogeneity of out-going edges [START_REF] Barthélemy | Characterization and modeling of weighted networks[END_REF]. When all connections eij leaving a given LOi have the same probability wij, the disparity will be close to the inverse of the out-degree 1/zouti otherwise it will be higher than this value. Column y2 shows that the disparity monotonically decreases with increasing K. From the point of view of a metaheuristic dynamic, a low disparity means that transitions are almost equiprobable. The LON topology does not preferentially guide the search trajectory, which makes the search process harder.

Column knn, reports the nearest-neighbors degree correlation or assortativity, a classical description of the mixing pattern of nodes in a network. The assortativity measures the affinity to connect with high or low-degree neighbors. The LONs of N K-landscapes are strongly disassortative, i.e. LO with few connections tend to link to others with many, and conversely. The implications of this observation on the search difficulty are worth further investigation.

The last column in Table 1 gives the values of the success performance indicator E(T ) detailed in Sec 2.3, which we abbreviate to ets in the following. Clearly, the expected running time increases with the landscape ruggedness (problem non-linearity) K. In order to analyze the inter-correlations between all the studied metrics, Figure 1 displays a correlogram of the whole data set. The figure depicts any possible pairing among the observed variables with a scatter plot in the panel below the diagonal, and the corresponding correlation coefficient in the upper panel. Table 1 corresponds to the first column of the correlation matrix, i.e. the one showing how the different network metrics shape against the epistasis K.

We are mostly interested in the correlations with the performance metric ets. Thus, the most relevant scatter plots are those in the last row of Figure 1 together with the respective Pearson coefficients in the last column. By inspecting these values, the highest positive correlation is with the average length of the shortest path climbing to the global optimum (r(ets, lo) 0.52). and with the total number of LO (r(ets, nv) 0.50). This observation seems reasonable: the more rugged the landscape (i.e. the larger its LON), the higher the number of hops through the LON to reach the global optimum, and thus the longer the expected running time of a restarting local search. Conversely, the larger the LO basins, and the higher the nearest-neighbors fitness-fitness correlation, the shorter the running time ((r(ets, wii) r(ets, f nn) -0.40)). The catch here, is that those metrics are in turn correlated among themselves (e.g. r(f nn, wii) 0.92), which prevents one from drawing relationships of causality from simple pair-wise correlations.

The last column plots also suggest a strong non-linearity of the performance estimator w.r.t. the LON metrics. In order to further analyze these relationships, Figure 2 zooms on the relevant scatter plots, and displays the logarithm of ets as a function of the considered landscape measures. The log-transformation allows to approach linearity, highlighting and confirming the results in the last column of Figure 1. Namely, the relationships between ets and all the metrics but lv appear clearly. The picture suggests a positive exponential trend with lo and zout, and a negative exponential trend with f nn, wii, cc, and y2. In the case of nv, the relation could also be close to power-law.

Since data are far from being normal bivariate, a more robust measure of association would be the rank-based Spearman's ρ statistic, reported in Table 2. These results complement the visual inspection of scatter plots and confirm the previous observations. 2, address already some of the research questions asked in Section 1, but do not provide an explanatory model for the algorithm performance as a function of the landscape features. To this end, we perform a multiple linear Correlation matrix for all pairs of the observed variables, to be read from the diagonal. The lower panel displays scatter plots and smoothing splines for every possible pairing; the upper panel gives the corresponding Pearson correlation coefficient, with text size proportional to its absolute value. In the lower panel, the smoothing is done through locally-weighted polynomial regressions [START_REF] Cleveland | Lowess: A program for smoothing scatterplots by robust locally weighted regression[END_REF]. In the upper panel, the correlation coefficient is tested against the null hypothesis and the resulting p-value is symbolically encoded at the levels of 0.1 ('), 0.05 (*), 0.01 (**), and 0.001 (***). q q q q q q q q 2 4 6 k log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2000 4000 6000 8000 10000 12000 14000 6 8 10 12 14 16 nv log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 50 100 150 6 8 10 12 14 16 lo log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 100 150 200 250 6 8 10 12 14 16 lv log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.2 0.4 0.6 0.8 6 8 10 12 14 16 fnn log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 40 60 80 100 120 6 8 10 12 14 16 wii log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.1 0.2 0.3 0.4 0.5 0.6 0.7 6 8 10 12 14 16 cc log(ets) q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 40 60 80 6 8 10 12 14 16 zout log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.6 6 8 10 12 14 16 y2 log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 6 8 10 12 14 16 knn log(ets)

Statistical Modeling

Figure 2: Scatter plots of the logarithm of the estimated time to succeed against all measured variables (see labels on x-axis). The epistasis value K is treated as a category, thus results are in a box-and-whisker diagram per group (fist plot on the top-left corner).

Table 1: Group averages of all the observed variables, aggregated by the epistasis K of the corresponding N K-landscape. Standard deviations are given in subscripts. nv = number of vertices (Local Optima), lo = average shortest path to reach the global optimum, lv = average path length (dij = 1/wij), f nn = Spearman coefficient for the nearest-neighbors fitness-fitness correlation, wii = average non-normalized weight of self-loops, cc global clustering coefficient, zout = average out-going degree, y2 = average weight disparity for out-going edges, knn = degree assortativity, ets = estimated time to succeed.

yi = β0 + β1xi,1 + β2xi,2 + • • • + βpxi,p + i (1)
where the response variable y in our case would be ets and p different predictors xj are to be chosen among the LON metrics; is the usual random noise term. The least square regression produces estimates βj for the βj model coefficients; the difference between the predicted values and the actual, observed values, are the regression residuals:

ei = yi -ŷi = yi -β1xi1 -β2xi2 -• • • -βpxip.
The difficulty of this analysis is that there are several possible explanatory variables, which are in turn intercorrelated. In consequence, some of them could have a confounding effect on the regression. In general, when confounders are known, measurable, and measured, it is a good practice to include them in the model. We, therefore, start by fitting the following formula:

log(ets) = β0 +β1k +β2 log(nv)+β2lo+• • •+β10knn+ (2)
where, w.r.t the general expression 1, the response and one of the predictors have been log-transformed in order to better approach linearity, as seen in Figure 2. Moreover, the variable k is qualitative and enters the model as a fixed effect, which translates in one appropriate dummy variable for each class K = 2, K = 4, . . . , K = 17.

The summary statistics for this model are reported in Table 3. In the table caption, the multiple R 2 represents the proportion of variance explained by the linear regression. R 2 would be equal to 1 if all observed data points were lying on the regression plane. When comparing models with a different number of predictors, the adjusted R 2 should be used instead. The F statistic is the ratio of the variance explained by the parameters in the model, to the residual or unexplained variance. The p-value is the probability of achieving an F that large under the null hypothesis of no effect [START_REF] Chambers | Linear models[END_REF].

The estimated coefficient βj and their estimated standard error σj are given in the 2 nd and 3 rd columns, respectively. Their ratio is the t-statistic (4 th column) that is used to calculate a p-value for the significance of the estimation (last column). In this initial model, the average length of paths to the global optimum lo is the only predictor with a regression coefficient that is statistically-significant at the 0.05 threshold (β lo = 0.04882, p-value = 2.37e -19).

Therefore, we proceed to perform a step-wise model selection by backward elimination [START_REF] Venables | Modern Applied Statistics with S[END_REF]. From the initial formula, at each step we compute what change in the fit could be produced by dropping each predictor in turn, and then we eliminate the one that minimizes the AIC score of the resulting model [START_REF] Sakamoto | Akaike information criterion statistics[END_REF]. By iterating until all predictors become significant, we obtain the final model:

log(ets) = β0 + β lo lo + βzoutzout + βy2y2 + β knn knn + (3)
which is detailed in Table 4. This final model is able to explain 84.94% of the variance observed in log(ets) with a linear regression on four variables that are all LON network metrics. Among these metrics, the length of the paths to the global optimum, and the weight disparity have the highest relative importance [START_REF] Groemping | Relative importance for linear regression in R: The package relaimpo[END_REF].

Without a check on the model assumptions, this would remain an observational study and could not be used to make predictions. To this end, a combination of parametric tests (not reported for space reasons) provided a positive confirmation [START_REF] Pena | gvlma: Global Validation of Linear Models Assumptions[END_REF]. However, a visual diagnostics is more informative. In particular, Figure 3 helps to assess if the regression residuals follow a normal distribution with zero mean and homogeneous variance, whereas Figure 4 displays the contributions to the model of each predictor in turn, highlightq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 6 8 10 12 14 16 -2 -1 0 1 2 Fitted values Residuals q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 t-Quantiles Studentized Residuals q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 3: Top: residual plot, to asses the hypothesis of zero-mean and constant-variance of the regression residuals (circle dots) around the fitted values (dotted line); no visually-significant deviation appears (red smooth line). Bottom: quantile-quantile comparison of the studentized regression residuals (circle dots) against the theoretical quantiles (red thick line), to inspect the distribution of residuals; no significant deviation from normality appears (confidence intervals as dotted red lines). Comp+Res log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 40 60 80

-2 -1 0 1 2 3 zout
Comp+Res log(ets) q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

0.1 0.2 0.3 0.4 0.5 0.6 -4 -3 -2 -1 0 1 2 3 y2
Comp+Res log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 Comp+Res log(ets) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 4: Component+residual plots for the linear regression model. Circle dots are, for each observation, the corresponding residual error from the regression plus the value fitted by one explanatory variable alone, plotted against that variable (see labels on x-axis). No significant deviation from linearity appears (smooth green line against the red dotted line of the partial regression).

ing possible violations of the linearity hypothesis [START_REF] Fox | An R Companion to Applied Regression[END_REF]. All assumptions seem acceptable. Therefore, formula 3 could be used to make inferences. In other words, formula 3 coefficients can be interpreted as conditional expectations for the average change in the response when one predictor undergoes a unitary change, and all the others remain fix. Since the dependent variable is log-transformed, this effect would be multiplicative.

The only important limitation of the proposed model is multicollinearity: due to the complex intercorrelations among LON metrics (cf. Figure 1), predictors are not really independent. This does not invalidate the multiple linear regression analysis, but it inflates the variance of its coefficients and makes it harder to disentangle their respective contributions [START_REF] Fox | An R Companion to Applied Regression[END_REF].

CONCLUSIONS

This article explored correlations between local optima network features and the performance of a stochastic local search algorithm running on the underlying combinatorial optimization problem. The NK family of landscapes and the iterated local search meta-heuristic were considered. It has been shown in previous work e.g. [START_REF] Tomassini | Complex-network analysis of combinatorial spaces: The NK landscape case[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF] that some features of the LON networks are related to the N Klandscapes ruggedness, and thus to problem difficulty. However, no statistically testable model was presented. The contribution of this study was to investigate, with a statistically-sound approach, which features of the LONs have a strong influence on the search performance, expressed as expected running times to success. The results obtained through the use of a multiple linear regression model show that some LON metrics are more important than others. These are: the average length of the shortest paths to the optimum, the average out-degree, the average disparity, and the degree assortativity. This study confirms and provides significant evidence that LON modeling is a compressed-but-relevant view of the fitness landscape, and can be used to understand and predict search difficulty.

It is worth noticing that some network metrics can be estimated without knowing the global optimum beforehand, such as the average out-degree, the fitness-fitness correlation, the average disparity and assortativity. Using these metrics and an adequate statistical model, as we have done in this work, opens up exciting possibilities. With standard sampling methods, larger search spaces could be studied. Thereafter, using the performance model based on the estimated LON metrics, the search heuristic parameters or its operators can be off-line tuned, or even on-line controlled. We plan to continue working in this direction and to extend this analysis to other combinatorial problems such as QAP.

Figures 1 and 2 ,

 2 Figures 1 and 2, along with Table2, address already some of the research questions asked in Section 1, but do not provide an explanatory model for the algorithm performance as a function of the landscape features. To this end, we perform a multiple linear

Figure 1 :

 1 Figure1: Correlation matrix for all pairs of the observed variables, to be read from the diagonal. The lower panel displays scatter plots and smoothing splines for every possible pairing; the upper panel gives the corresponding Pearson correlation coefficient, with text size proportional to its absolute value. In the lower panel, the smoothing is done through locally-weighted polynomial regressions[START_REF] Cleveland | Lowess: A program for smoothing scatterplots by robust locally weighted regression[END_REF]. In the upper panel, the correlation coefficient is tested against the null hypothesis and the resulting p-value is symbolically encoded at the levels of 0.1 ('), 0.05 (*), 0.01 (**), and 0.001 (***).

Table 2 :

 2 Spearman's ρ statistic for the correlation between ets and the LON metrics (p-value < 2.2e -16 for all pairings).

	aggregate	nv	lo		lv	fnn		wii	cc	zout	y2	knn	ets (×10 4 )
	K = 2		4328	33.514	18751	0.7030.19	10511	0.4250.086	6.91.8	0.3920.075	0.1550.4	2.163.3
	K = 4		22139	53.712	21415	0.5870.07	83.93	0.2630.013	14.31	0.2190.016	-0.5360.13	8.397.74
	K = 6		74870	66.713	1884.8	0.5350.041	67.51.9	0.190.005	24.80.86	0.1240.0059	-0.7780.035	32.636.1
	K = 8		166973	76.69.1 1711.9	0.4310.025 53.30.88 0.1590.0012 35.70.57	0.07690.002	-0.8560.022	51.861.1
	K = 10	3148110	90.78.4 1661.2	0.3420.016 40.70.78 0.1430.00085 47.20.57	0.04910.0011	-0.9040.011	81.570.6
	K = 12	5270104	10812	1700.64 0.2550.015 30.80.35 0.1330.00054 57.80.39 0.03340.00046 -0.9280.0093	276544
	K = 14	8100121	1258.6	1810.6	0.190.011	23.50.25 0.1280.00032 66.90.33 0.02450.00022 -0.9440.0063	300288
	K = 16	11688101	14611	1970.42 0.1430.0073 18.20.11 0.1250.00023 74.60.17 0.01967.8e-05 -0.9480.0055	414632
	K = 17	1380174	15612	2050.42	0.1330.01	16.10.06 0.1250.00021 78.20.13 0.01796.5e-05 -0.9440.0063	793844
	nv	lo	lv	fnn	wii	cc	zout	y2	knn
	0.885 0.915 0.006 -0.830 -0.883 -0.875 0.885 -0.883 -0.850
	regression on the data, which has the general form:

Table 3 :

 3 Summary statistics of the linear regression model on all variables. Residual standard error: 0.8702 on 248 degrees of freedom (4 observations deleted due to missingness). Multiple R-squared: 0.8585, Adjusted R-squared: 0.8488. Fstatistic: 88.52 on 17 and 248 DF, p-value: < 2.2e -16.

	summary	Estimate	Std. Error	t value	Pr(>|t|)
	(Intercept)	16.06966	7.39640	2.1726 3.08e -02
	k4	0.01542	0.63987	0.0241 -01
	k6	-1.08926	1.39976 -0.7782 4.37e -01
	k8	-3.14529	2.46616 -1.2754 2.03e -01
	k10	-5.67316	3.76828 -1.5055 1.33e -01
	k12	-8.19327	5.04638 -1.6236 1.06e -01
	k14	-10.34765	6.17715 -1.6751 9.52e -02
	k16	-12.85523	7.14074 -1.8003 7.30e -02
	k17	-13.40456	7.59325 -1.7653 7.87e -02
	log(nv)	-1.91370	1.12656 -1.6987 9.06e -02
	lo	0.04882	0.00499	9.7919 2.37e -19
	lv	0.00198	0.00376	0.5265 5.99e -01
	fnn	0.54148	0.89574	0.6045 5.46e -01
	wii	-0.00302	0.02739 -0.1104 9.12e -01
	cc	-7.22853	5.00042 -1.4456 1.50e -01
	zout	0.29514	0.15838	1.8636 6.36e -02
	y2	-3.46837	5.00914 -0.6924 4.89e -01
	knn	-0.88961	0.49062 -1.8132 7.10e -02

Table 4 :

 4 Summary statistics of the final linear regression model. Residual standard error: 0.8751 on 261 degrees of freedom (4 observations deleted due to missingness). Multiple R-squared: 0.8494, Adjusted R-squared: 0.8471. F-statistic: 368.1 on 4 and 261 DF, p-value: < 2.2e -16.

	summary	Estimate Std. Error t value	Pr(>|t|)
	(Intercept) 10.3838	0.58512	17.75 9.24e -47
	lo	0.0439	0.00434	10.11 1.67e -20
	zout	-0.0306	0.00831 -3.68 2.81e -04
	y2	-7.2831	1.63038 -4.47 1.18e -05
	knn	-0.7457	0.40501 -1.84 6.67e -02

probability distribution of the number N -1 of failures before the first success (Bernoulli trials)

ACKNOWLEDGMENTS

Fabio Daolio and Marco Tomassini gratefully acknowledge the Swiss National Science Foundation for financial support under grant number 200021-124578.