Exact Algorithms for Maximum Clique: a computational study

We investigate a number of recently reported exact algorithms for the maximum clique problem (MCQ, MCR, MCS, BBMC). The program code used is presented and critiqued showing how small changes in implementation can have a drastic effect on performance. The computational study demonstrates how problem features and hardware platforms influence algorithm behaviour. The minimum width order (smallest-last) is investigated, and MCS is broken into its consituent parts and we discover that one of these parts degrades performance. It is shown that the standard procedure used for rescaling published results is unsafe.


Introduction
The purpose of this paper is to investigate a number of recently reported exact algorithms for the maximum clique problem. The actual program code used is presented and critiqued. The computational study aims to show how implementation details, problem features and hardware platforms influence algorithmic behaviour in those algorithms.
The Maximum Clique Problem (MCP): A simple undirected graph G is a pair (V, E) where V is a set of vertices and E a set of edges. An edge {u, v} is in E if and only if {u, v} ⊆ V and vertex u is adjacent to vertex v. A clique is a set of vertices C ⊆ V such that every pair of vertices in C is adjacent in G. Clique is one of the six basic NP-complete problems given in [10]. It is posed as a decision problem [GT19]: given a simple undirected graph G = (V, E) and a positive integer k ≤ |V | does G contain a clique of size k or more? The optimization problems is then to find the maximum clique, where ω(G) is the size of a maximum clique.
A graph can be coloured, by that we mean that any pair of adjacent vertices must be given different colours. We do not use colours, we use integers to label the vertices. The minimum number of different colours required is then the chromatic number of the graph χ(G), and ω(G) ≤ χ(G). Finding the chromatic number is NP-complete.
Exact Algorithms for MCP: We can address the decision and optimization problems with an exact algorithm, such as a backtracking search [8,20,31,6,18,17,22,21,12,26,27,14,5]. Backtracking search incrementally constructs the set C (initially empty) by choosing a candidate vertex from the candidate set P (initially all of the vertices in V ) and then adding it to C. Having chosen a vertex the candidate set is then updated, removing vertices that cannot participate in the evolving clique. If the candidate set is empty then C is maximal (if it is a maximum we save it) and we then backtrack. Otherwise P is not empty and we continue our search, selecting from P and adding to C.
There are other scenarios where we can cut off search, i.e. if what is in P is insufficient to unseat the champion (the largest clique found so far) search can be abandoned. That is, an upper bound can be computed. Graph colouring can be used to compute an upper bound during search, i.e. if the candidate set can be coloured with k colours then it can contain a clique no larger than k [31,8,22,12,26,27]. There are also heuristics that can be used when selecting the candidate vertex, different styles of search, different algorithms to colour the graph and different orders in which to do this.
Structure of the Paper: In the next section we present in Java the following algorithms: Fahle's Algorithm 1 [8], Tomita's MCQ [26], MCR [25] and MCS [27] and San Segundo's BBMC [22]. By using Java and its inheritance mechanism algorithms are presented as modifications of previous algorithms. Three vertex orderings are presented, one being new to these algorithms. Starting with the basic algorithm MC we show how minor coding details can significantly impact on performance. Section 3 presents a chronological review of exact algorithms, starting at 1990. Section 4 is the arXiv:1207.4616v1 [cs.DS] 19 Jul 2012 computational study. The study investigates MCS and determines where its speed advantage comes from, measures the benefits resulting from the bit encoding of BBMC, the effectiveness of three vertex orderings and the potential to be had from tie-breaking within an ordering. New benchmark problems are then investigated. Finally an established technique for calibrating and scaling results is put to the test and is shown to be unsafe. Finally we conclude.

The Algorithms: MC, MCQ, MCS and BBMC
We start by presenting the simplest algorithm [8] which I will call MC. This sets the scene. It is presented as a Java class, as are all the algorithms, with instance variables and methods. And we might say that I am abusing Java by using the class merely as a place holder for global variables and methods for the algorithms, and inheritance merely to over-ride methods so that I can show a program-delta, i.e. to present the differences between algorithms.
Each algorithm is first described textually and then the actual implementation is given in Java. Sometimes a program trace is given to better expose the workings of the algorithm. It is possible to read this section skipping the Java descriptions, however the Java code makes it explicit how one algorithm differs from another and shows the details that can severely affect the performance of the algorithm.
We start with MC. MC is essentially a straw man: it is elegant but too simple to be of any practical worth. Nevertheless, it has some interesting features. MCQ is then presented as an extension to MC, our first algorithm that uses a tight integration of search algorithm, search order and upper bound cut off. Our implementation of MCQ allows three different vertex orderings to be used, and one of these corresponds to MCR. The presentation of MCQ is somewhat laborious but this pays off when we present two variants of MCS as minor changes to MCQ. BBMC is presented as an extension of MCQ, but is essentially MCSa with sets implemented using bit strings. Figure  1 shows the hierarchical structure for the algorithms presented. Appendix 1 gives a description of the execution environment (the MaxClique class) that allows us to read clique problems and solve them with a variety of algorithms and colour orderings. It also shows how to actually run our programs from the command line.

MC
We start with a simple algorithm, similar to Algorithm 1 in [8]. Fahle's Algorithm 1 uses two sets: C the growing clique (initially empty) and P the candidate set (initially all vertices in the graph). C is maximal when P is empty and if |C| is a maxima it is saved, i.e. C becomes the champion. If |C| + |P | is too small to unseat the champion search can be terminated. Otherwise search iterates over the vertices in P in turn selecting a vertex v, creating a new growing clique C where C = C ∪ {v} and a new candidate set P as the set of vertices in P that are adjacent to v (i.e. P = P ∩ neighbours(v)), and recursing. We will call this MC.
MC in Java Listing 1.1 can be compared to Algorithm 1 in [8]. The constructor, lines 14 to 22, takes three arguments: n the number of vertices in the graph, A the adjacency matrix where A[i][j] equals 1 if and only if vertex i is adjacent to vertex j, and degree where degree[i] is the number of vertices adjacent to vertex i (and is the sum of A[i]). The variables nodes and cpuT ime are used as measures of search performance, timeLimit is a bound on the run-time, maxSize is the size of the largest clique found so far, style is used as a flag to customise the algorithm with respect to ordering of vertices and the array solution is the largest clique found such that solution[i] is equal to 1 if and only if vertex i is in the largest clique found.
The method search() finds a largest clique 1 or terminates having exceeding the allocated timeLimit. Two sets are produced: the candidate set P and the current clique C 2 . Vertices from P may be selected and added to the growing clique C. Initially all vertices are added to P and C is empty (lines 27 and 28). The sets P and C are represented using Java's ArrayList, a re-sizablearray implementation of the List interface. Adding an item is an O(1) operation but removing an arbitrary item is of O(n) cost. This might appear to be a damning indictment of this simple data structure but as we will see it is the cost we pay if we want to maintain order in P and in many cases we can work around this to enjoy O(1) performance.
The search is performed in method expand 3 . In line 34 a test is performed to determine if the cpu time limit has been exceeded, and if so search terminates. Otherwise we increment the number of nodes, i.e. a count of the size of the backtrack search tree explored. The method then iterates over the vertices in P (line 36), starting with the last vertex in P down to the first vertex in P . This form of iteration over the ArrayList, getting entries with a specific index, is necessary when entries are deleted (line 45) as part of that iteration. A vertex v is selected from P (line 38), added to C (line 39), and a new candidate set is then created (line 40) newP where newP is the set of vertices in P that are adjacent to vertex v (line 41). Consequently all vertices in newP are adjacent to all vertices in C and all pairs of vertices in C are adjacent (i.e. C is a clique). If newP is empty C is maximal and if it is the largest clique found it is saved (line 42). If newP is not empty then C is not maximal and search can proceed via a recursive call to expand (line 43). On returning from the recursive call v is removed from P and from C (lines 44 and 45).
There is one "trick" in expand and that is at line 37: if the combined size of the current clique and the candidate set cannot unseat the best clique found so far this branch of the backtrack tree can be abandoned. This is the simplest upper bound cut-off and corresponds to line 3 from Algorithm 1 in [8]. The method saveSolution does as it says: it saves off the current maximal clique and records its size.
Observations on MC There are several points of interest. This first is the search process itself. If we commented out lines 37 and changed line 41 to add to newP all vertices in P other than v, method expand would produce the power set of P and at each depth k in the backtrack tree we would have n k calls to expand. That is, expand produces a binomial backtrack search tree of size O(2 n ) (see page 6 and 7 of [11]). This can be compared to a bifurcating search process, where on one side we take an element and make a recursive call, and on the other side reject it and make a recursive call, terminating when P is empty. This generates the power set on the leaf nodes of the backtrack tree and explores 2 n+1 − 1 nodes. This is also O(2 n ) but in practice is often twice as slow as the binomial search. In Figure 2 we see a binomial search produced by a simplification of MC, generating the power set of {0, 1, 2, 3}. Each node in the tree contains two sets: the set that will be added to the power set and the set that can be selected from at the next level. We see 16 nodes and at each depth k we have n k nodes. The corresponding tree for the bifurcating search (not shown) has 31 nodes with the power set appearing on the 16 leaf nodes at depth 4.
The second point of interest is the actual Java implementation. Java gives us an elegant construct for iterating over collections, the for-each loop, used in line 41 of Listing 1.1. This is rewritten § 1 import j a v a . u t i l . A r r a y L i s t <I n t e g e r > C = new A r r a y L i s t <I n t e g e r >() ; 28 A r r a y L i s t <I n t e g e r > P = new A r r a y L i s t <I n t e g e r >(n ) ; 29 f o r ( i n t i =0; i <n ; i ++) P . add ( i ) ; 30 expand (C, P) ; 31 } 32 33 void expand ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > P)  and if it is adjacent to v it is added to newP (line 17 of Listing 1.2). In MC (line 41 of Listing 1.1) the for-each statement implicitly creates an iterator object and uses that for selecting elements. This typically results in a 10% reduction in runtime for MC0.
Our third point is how we create our sets. In MC0 line 15 the new candidate set is created with a capacity of i. Why do that when we can just create newP with no size and let Java work it out dynamically? And why size i? In the loop of line 10 i counts down from the size of the candidate set, less one, to zero. Therefore at line 14 P is of size i + 1 and we can set the maximum size of newP accordingly. If we do not set the size Java will give newP an initial size of 10 and when additions exceed this newP will be re-sized. By grabbing this space we avoid that. This results in yet another measurable reduction in run-time.
Our fourth point is how we remove elements from our sets. In MC we remove the current vertex v from C and P (lines 44 and 45) whereas in MC0 we remove the last element in C and P (lines 21 and 22). Clearly v will always be the last element in C and P . The code in MC results in a sequential scan to find and then delete the last element, i.e. O(n), whereas in MC0 it is a simple O(1) task. And this raises another question: P and C are really stacks so why not use a Java Stack? The Stack class is represented using an ArrayList and cannot be initialised with a size, but has a default initial size of 10. When the stack grows and exceeds its current capacity the capacity is doubled and the contents are copied across. Experiments showed that using a Stack increased run time by a few percentage points.
Typically MC0 is 50% faster than MC. In many cases a 50% improvement in run time would be considered a significant gain, usually brought about by changes in the algorithm. Here, such a gain can be achieved by moderately careful coding. And this is our first lesson: when comparing published results we need to be cautious as we may be comparing programmer ability as much as differences in algorithms.
The fifth point is that MC makes more recursive calls that it needs to. At line 37 |C| + |P | is sufficient to proceed but at line 43 it is possible that |C| + |newP | is actually too small and will generate a failure at line 37 in the next recursive call. We should have a richer condition at line 43 but as we will soon see the algorithms that follow do not need this.
The sixth point is a question of space: why is the adjacency matrix an array of integers when we could have used booleans, surely that would have been more space efficient? In Java a boolean is represented as an integer with 1 being true, everything else false. Therefore there is no saving in space and only a minuscule saving in time (more code is generated to test if A[i][j] equals 1 than to test if a boolean is true). Furthermore by representing the adjacency matrix as integers we can sum a row to get the degree of a vertex.
And finally, Listing 1.1 shows exactly what is measured. Our run time starts at line 25, at the start of search. This will include the times to set up the data structures peculiar to an algorithm, and any reordering of vertices. It does not include the time to read in the problem or the time to write out a solution. There is also no doubt about what we mean by a node: a call to expand counts as one more node. void expand ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > P) i n t v = P . g e t ( i ) ; 13 C . add ( v ) ; 14 A r r a y L i s t <I n t e g e r > newP = new A r r a y L i s t <I n t e g e r >( i ) ; 15 f o r ( i n t j =0; j<=i ; j ++){ 16 i n t w = P . g e t ( j ) ; 17 i The problem is referred to as g10-50, and is a randomly generated graph with 10 vertices with edge probability 0.5. This is shown in Figure 3 and has at top a cartoon of the search process and immediately below a trace of our program. The program prints out the arguments C and P on each call to expand (between lines 33 and 34), in line 37 if a FAIL occurs, and between lines 38 and 39 when a vertex is selected. The indentation corresponds to the depth of recursion. The Line dd boxes in the cartoon of Figure 3 corresponds to the line numbers in the trace of Figure  3, each of those a call to expand. Green coloured vertices are in P , blue vertices are those in C and red vertices are those removed from P and C in lines 44 and 45 of Listing 1.1. Also shown is the backtrack tree. The boxes correspond to calls to expand and contain C and P . On arcs we have numbers with a down arrow ↓ if that vertex is added to C and an up arrow ↑ if that vertex is removed from C and P . A clear white box is a call to expand that is an interior node of the backtrack tree leading to further recursive calls or the creation of a new champion. The green shriek! is a champion clique and a red shriek! a fail because |C| + |P | was too small to unseat the champion. The blue boxes correspond to calls to expand that fail first time on entering the loop at line 36 of Listing 1.1. By looking at the backtrack tree we get a feel for the nature of binomial search.
Method expand (line 23 to 43) corresponds to Figure 2 in [26]. The array colour is local to the method and holds the colour of the i th vertex in P . The candidate set P is then sorted in non-decreasing colour order by the call to numberSort in line 28, and colour[i] is then the colour of integer vertex P.get(i). The search then begins in the loop at line 29. We first test to see if the combined size of the candidate set plus the colour of vertex v is sufficient to unseat the champion (the largest clique found so far) 4 (line 30). If it is insufficient the search terminates. Note that the loop starts at m − 1, the position of the last element in P , and counts down to zero. The i th element of P is selected and assigned to v. As in MC we create a new candidate set newP , the set of vertices (integers) in P that are adjacent to v (lines 33 to 37). We then test to see if C is maximal (line 38) and if it unseats the champion. If the new candidate set is not empty we recurse (line 39). Regardless, v is removed from P and from C (lines 40 and 41) as in (lines 21 and 22 of MC0).
Method numberSort can be compared to Figure 3 in [26]. numberSort takes as arguments an ordered ArrayList of integers ColOrd corresponding to vertices to be coloured in that order, an ArrayList of integers P that will correspond to the coloured vertices in non-decreasing colour order, and an array of integers colour such that if v = P.get(i) (v is the i th vertex in P ) then the colour of v is colour[i]. Lines 45 to 64 differs from Tomita's NUMBER-SORT method because we use the additional arguments ColOrd and the growing clique C as this allows us to easily implement our next algorithm MCS.
Rather than assign colours to vertices explicitly numberSort places vertices into colour classes, i.e. if a vertex is not adjacent to any of the vertices in colourClass[i] then that vertex can be placed into that class and given colour number i + 1 (i + 1 so that colours range from 1 upwards). The vertices can then be sorted into colour order via a pigeonhole sort, where colour classes are the pigeonholes.
numberSort starts by clearing out the colour classes that might be used (line 48). In lines 49 to 55 vertices are selected from ColOrd and placed into the first colour class in which there are no conflicts, i.e. a class in which the vertex is not adjacent to any other vertex in that class (lines 51 to 53, and method conf licts). Variable colours records the number of colours used. Lines 56 to 63 is in fact a pigeonhole sort, starting by clearing P and then iterating over the colour classes (loop start at line 58) and in each colour class adding those vertices into P (lines 59 to 63). The boolean method conf licts, lines 66 to 72, takes a vertex v and an ArrayList of vertices colourClass where vertices in colourClass are not pair-wise adjacent and have the same colour i.e. the vertices are an independent set. If vertex v is adjacent to any vertex in colourClass the method returns true (lines 67 to 70) otherwise false. Note that if vertex v needs to be added into a new colour class in numberSort the size of that colourClass will be zero, the for loop of lines 67 to 70 will not be performed and conf licts returns true. The complexity of numberSort is quadratic in the size of P .
Vertices need to be sorted at the top of search, line 19. To do this we use the class Vertex in Listing 1.4 and the comparator MCRComparator in Listing 1.5. If i is a vertex in P then the corresponding Vertex v has an integer index equal to i. The Vertex also has attributes degree and nebDeg. degree is the degree of the vertex index and nebDegree is the sum of the degrees of the vertices in the neighbourhood of vertex index. Given an array V of class Vertex this can be sorted § A r r a y L i s t <I n t e g e r > C = new A r r a y L i s t <I n t e g e r >(n ) ; 17 A r r a y L i s t <I n t e g e r > P = new A r r a y L i s t <I n t e g e r >(n ) ; void expand ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > P) void numberSort ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > ColOrd , A r r a y L i s t <I n t e g e r > P , using Java's Arrays.sort(V ) method in O(n. log(n)) time, and is ordered by default using the compareT o method in class Vertex. Our method forces a strict ordering of V by non-increasing degree, tie-breaking on index. And this is one step we take to ensure reproducibility of results. If we allowed the compareT oM ethod to deliver 0 when two vertices have the same degree then Arrays.sort would then break ties. If the sort method was unstable, i.e. did not maintain the relative order of objects with equal keys [3], results may be unpredictable. § The class MCRComparator (Listing 1.5) allows us to sort vertices by non-increasing degree, tie breaking on the sum of the neighbourhood degree nebDeg and then on index, giving again a strict order. This is the MCR order given in [25], where MCQ uses the simple degree ordering and MCR is MCQ with tie-breaking on neighbourhood degree. § Vertices can also be sorted into a minimum-width order via method minW idthOrder of lines 86 to 98 of Listing 1.6. The minimum width order (mwo) was proposed by Freuder [9] and also by Matula and Beck [16] where it was called "smallest last", and more recently in [7] as a degeneracy ordering. The method minW idthOrder, lines 86 to 98 of Listing 1. is computed in lines 77 to 79. ColOrd is then sorted in one of three orders: style == 1 in nonincreasing degree order, style == 2 in minimum width order, style == 3 in non-increasing degree tie-breaking on sum of neighbourhood degree. MCQ then uses the ordered candidate set P for colouring, initially in one of the initial orders, thereafter in the order resulting from numberSort and that is non-decreasing colour order. In [26] it is claimed that this is an improving order (however, no evidence was presented for this claim). In [25] Tomita proposes a new algorithm, MCR, where MCR is MCQ with a different initial ordering of vertices. Here MCR is MCQ with style = 3. Figure 4 shows a cartoon and trace of MCQ over graph g10-50. Print statements were placed immediately after the call to expand (Listing 1.3 line 24), after the selection of a vertex v (line 31) and just before v is rejected from P and C (line 40). Each picture in the cartoon gives the corresponding line numbers in the trace immediately below. Line 0 of the trace is a print out of the ordered array V just after line 83 in method orderV ertices in Listing 1.6. This shows for each vertex the pair < index, degree >: the first call to expand has P = {3, 0, 4, 6, 1, 2, 5, 8, 9, 7}, i.e. non-decreasing degree order. MCQ makes 3 calls to expand whereas MC makes 9 calls, and the MCQ colour bound cut off in line 30 of Listing 1.3 is satisfied twice (listing lines 9 and 11).

A trace of MCQ
Observations on MCQ We noted above that MC can make recursive calls that immediately fail. Can this happen in MCQ? Looking at lines 39 of Listing 1.3, |C| + |newP | must be greater than maxSize. Since the colour of the vertex selected colour[i] was sufficient to satisfy the condition of line 30 it must be that integer vertex v (line 31) is adjacent to at least colour[i] vertices in P and thus in newP , therefore the next recursive call will not immediately fail. Consequently each call to expand corresponds to an internal node in the backtrack tree.
We also see again exactly what is measured as cpu time: it includes the creation of our data structures, the reordering of vertices at the top of search and all recursive calls to expand (lines 12 to 20).
Why is colourClass an ArrayList[] rather than an ArrayList<ArrayList<Integer>> That would have done away with the explicit cast in line 60. When using an ArrayList<ArrayList<Integer>> we can do away with the cast but Java generates an implicit cast, so nothing is gained. At the top of MCQ's search Tomita sorts vertices in non-increasing degree order and the first ∆ vertices are given colours 1 to ∆ respectively, where ∆ is the maximum degree in the graph, thereafter vertices are given colour ∆ + 1. This is done to prime the candidate set for the initial call to EXPAND. Thereafter Tomita calls NUMBER-SORT immediately before the recursive call to EXPAND. A simpler option is taken here: colouring and sorting of the candidate set is done at the start of expand (call to numberSort at line 28). This strategy is adopted in all the algorithms here.

MCS
Tomita's MCR [25] is MCQ with a richer initial ordering, essentially non-decreasing degree with tie-breaking on the sum of the degrees of adjacent vertices. This ordering is then modified during search via the colouring routine numberSort as in MCQ. MCR is compared to MCQ in [25] over 8 of the 66 instances of the DIMACS benchmarks [1] showing an improvement in MCR over MCQ. As previously stated, MCR is our MCQ with style = 3.
MCS [27] is MCR with two further modifications. The first modification is that MCS uses "... an adjunct ordered set of vertices for approximate coloring". This is an ordered list of vertices to be used in the sequential colouring, and was called V a . This order is static, set at the top of search. Therefore, rather than use the order in the candidate set P for colouring the vertices in P the vertices in P are coloured in the order of vertices in V a .
The second modification is to use a repair mechanism when colouring vertices (this is called a Re-NUMBER in Figure 1 of [27]). When colouring vertices an attempt is made to reduce the colours used by performing exchanges between vertices in different colour classes. This is similar to the Sequential-X algorithm of Maffray and Preissmann [15], i.e. a bi-chromatic exchange is performed between a subset of vertices in a pair of colour classes indifferent to a given vertex. In [27] a recolouring of a vertex v occurs when a new colour class is about to be opened for v and that colour class exceeds the search bound, i.e. if the number of colours can be reduced this could result in search being cut off. In the context of colouring I will say that vertex u and v conflict if they are adjacent, and that v conflicts with a colour class C if there exists a vertex u ∈ C that is in conflict with v. Assume vertex v is in colour class C k . If there exists a lower colour class C i (i < k − 1) and v conflicts only with a single vertex w ∈ C i and there also exists a colour class C j , where i < j < k, and w does not conflict with any vertex in C j then we can place v in C i and w in C j , freeing up colour class C k . This is given in Figure 1 of [27] and the procedure is named Re-NUMBER. Figure 5 illustrates this procedure. The boxes correspond to colour classes i, j and k where i < j < k. The circles correspond to vertices in that colour class and the red arrowed lines as conflicts between pairs of vertices. Vertex v has just been added to colour class k, v conflicts only with w in colour class i and w has no conflicts in colour class j. We can then move w to colour class j and v to colour class i. Experiments were then presented in [27] comparing MCR against MCS in which MCS is always the champion. But it is not clear where the advantage of MCS comes from: does it come from the static colour order (the "adjunct ordered set") or does it come from the colour repair mechanism?
I now present two versions of MCS. The first I call MCSa and this uses the static colouring order. The second, MCSb, uses the static colouring ordering and the colour repair mechanism (so MCSb is Tomita's MCS). Consequently we will be able to determine where the improvement in MCS comes from: static colour ordering or colour repair.
MCSa in Java In Listing 1.7 we present MCSa as an extension to MCQ. Method search creates an explicit colour ordering ColOrd and the expand method is called with this in line 18 (compare this to line 20 of M CQ). Method expand now takes three arguments: the growing clique C, the candidate set P and the colouring order ColOrd. In line 26 numberSort is called using ColOrd (compare to line 28 in MCQ) and lines 27 to 45 are essentially the same as lines 29 to 42 in MCQ with the exception that ColOrd must also be copied and updated (line 32, 36 and 37) prior to the recursive call to expand (line 40) and then down-dated after the recursive call (line 43). Therefore MCSa is a simple extension of MCQ and like MCQ has three styles of ordering.
MCSb in Java In Listing 1.8 we present MCSb as an extension to MCSa: the difference between MCSb and MCSa is in numberSort, with the addition of lines 10 and 20. At line 10 we compute delta as the minimum number of colour classes required to match the search bound. At line 20, if we have exceeded the number of colour classes required to exceed the search bound and a new colour class k has been opened for vertex v and we can repair the colouring such that one less colour class is used we can decrement the number of colours used. This repair is done in the boolean method repair of lines 43 to 57. The repair method returns true if vertex v in colour class k can be recoloured into a lower colour class, false otherwise, and can be compared to Tomita's Re-NUMBER procedure. We search for a colour class i, where i < k − 1, in which there exists only one vertex in conflict with v and we call this w (line 45). The method getSingleConf lictV ariable, lines 32 to 41, takes as arguments a vertex v and a colour class and counts and records the conflict vertex. If this exceeds 1 we return a negative integer otherwise we deliver the single conflicting vertex. The repair method then proceeds at line 46 if a single conflicting vertex w was found, searching for a colour class j above i (for loop of line 47) in which there are no conflicts with w. If that was found (line 48) vertex v is removed from colour class k, w is removed from colour class i, v is added to colour class i and w to colour class j (lines 49 to 52) and repair delivers true (line 53). Otherwise, no repair occurred (line 56). A r r a y L i s t <I n t e g e r > C = new A r r a y L i s t <I n t e g e r >(n ) ; 14 A r r a y L i s t <I n t e g e r > P = new A r r a y L i s t <I n t e g e r >(n ) ; 15 A r r a y L i s t <I n t e g e r > ColOrd = new A r r a y L i s t <I n t e g e r >(n ) ; void expand ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > P , A r r a y L i s t <I n t e g e r > ColOrd ) A r r a y L i s t <I n t e g e r > newP = new A r r a y L i s t <I n t e g e r >( i ) ; 32 A r r a y L i s t <I n t e g e r > newColOrd = new A r r a y L i s t <I n t e g e r >( i ) ; 33 f o r ( i n t j =0; j<=i ; j ++){ 34 i n t u = P . g e t ( j ) ; super ( n , A, d e g r e e , s t y l e ) ; 7 } 8 9 void numberSort ( A r r a y L i s t <I n t e g e r > C, A r r a y L i s t <I n t e g e r > ColOrd , A r r a y L i s t <I n t e g e r > P , Observations on MCS Tomita did not investigate where MCS's improvement comes from and neither did [5], coding up MCS in Python in one piece. We can also tune MCS. In MCSb we repair colourings when we open a new colour class that exceeds the search bound. We could instead repair unconditionally every time we open a new colour class, attempting to maintain a compact colouring. We do not investigate this here.

BBMC
San Segundo's BB-MaxClique algorithm [22] (BBMC) is similar to the earlier algorithms in that vertices are selected from the candidate set to add to the current clique in non-increasing colour order, with a colour cut-off within a binomial search. BBMC is at heart a bit-set encoding of MCSa with the following features.
1. The "BB" in "BB-MaxClique" is for "Bit Board". Sets are represented using bit strings. 2. BBMC colours the candidate set using a static sequential ordering, the ordering set at the top of search, the same as MCSa. 3. BBMC represents the neighbourhood of a vertex and its inverse neighbourhood as bit strings, rather than use an adjacency matrix and its compliment 4. When colouring takes place a colour class perspective is taken, determining what vertices can be placed in a colour class together, before moving on to the next colour class. Other algorithms (e.g [26,27]) takes a vertex perspective, deciding on the colour of a vertex.
The candidate set P is encoded as a bit string, as is the currently growing clique C. For a given vertex v we also have its neighbourhood N [v], again encoded as a bit-string, and its inverse ) and this is used in colouring. When a vertex v is selected from the candidate set P to be added to C a new candidate set is produced newP , where newP is the set of vertices in the candidate set P that are adjacent to v i.e. newP = P ∧ N [v]. For the set elements that reside in word boundaries this operation is fast. BBMC takes a "colour-class" perspective. BBColour starts off building the first colour class, i.e. the first set of vertices in P that form an independent set. It then finds the next colour class, and so on until P is exhausted. Given a set of vertices Q, when a vertex v is selected and removed from Q and added to a colour class, Q becomes the set of vertices that are in Q but not adjacent to v. That is Q = Q ∧ invN [v]. Colour classes are then combined using a pigeonhole sort delivering a list of vertices in non-decreasing colour order and this is then used in the BBM axClique method (the BBMC equivalent of expand) to cut off search as in MCQ and MCSa.
MCSa colours vertices in a static order. This is achieved in BBMC by a renaming of the vertices, and this is done by re-ordering the adjacency matrix at the top of search.
BBMC in Java We implement sets using Java's BitSet class, therefore from now on we refer to P as the candidate BitSet and an ordered array of integers U as the ordered candidate set. In Listing 1.9, lines 5 to 7, we have the an array of BitSet N for representing neighbourhoods, invN as the inverse neighbourhoods (the compliment of N ) and V an array of Vertex. N [i] is then a BitSet representing the neigbourhood of the i th vertex in the array V , and invN as its compliment. The array V is used at the top of search for renaming vertices (and we discuss this later).
The search method (lines 16 to 30) creates the candidate BitSet P , current clique (as a BitSet) C, and Vertex array V . The orderV ertices method renames the vertices and will be discussed later. The method BBM axClique corresponds to the procedure in Figure 3 of [22] and can be compared to the expand method in Listing 1.7. In a BitSet we use cardinality rather than size (line 35, 40 and 44). The integer array U (same name as in [22] is essentially the colour ordered candidate set such that if v = U If the current clique is both maximal and a maximum it is saved via BBMC's specialised save method (described later) otherwise if C is not maximal (i.e. newP is not empty) a recursive call is made to BBM axClique. Regardless, v is removed from the current candidate BitSet and the current clique (line 46) and the for loop continues.
Method BBColour corresponds to the procedure of the same name in Figure 2 of [22] but differs in that it does not explicitly represent colour classes and therefore does not require a pigeonhole sort as in San Segundo's description. Our method takes the candidate BitSet P (see line 38), ordered candidate set U and array of colour as parameters. Due to the nature of Java's BitSet the and operation is not functional but actually modifies bits, consequently cloning is required (line 51) and we take a copy of P . In line 52 colourClass records the current colour class, initially zero, and i is used as a counter for adding coloured vertices into the array U . BBMC requires its own saveSolution method (lines 86 to 90 of Listing 1.10) due to C being a BitSet. Again the solution is saved into the integer array solution and again we need to use the Vertex array V to map bits to vertices. This is done in line 88: if the i th bit of C is true then integer vertex V [i].index is in the solution. This explains why V is global to the BBM C class.
Observations on BBMC In our Java implementation we might expect a speed up if we did away with the in-built BitSet and did our own bit-sting manipulations explicitly. It is also worth noting that in [21] comparisons are made with Tomita's results in [27] by re-scaling tabulated results, i.e. Tomita's code was not actually run. This is not unusual.
In [21] there is the bit-board version, BB ReCol in Fig 1, of Tomita's Re-NUMBER. I believe that version reported is flawed and can result in colour classes not being pair-wise disjoint, and that Fig. 1 should have a return statement between lines 7 and 8, similar to the return statement in Re-NUMBER. As it stands BB ReCol can result in the candidate set becoming a multi-set resulting in redundant re-exploration of the search space with subsequent poor performance.

Summary of MCQ, MCR, MCS and BBMC
Putting aside the chronology [26,25,27,22], MCSa is the most general algorithm presented here. BBMC is in essence MCSa with a BitSet encoding of sets. MCQ is MCSa except that we do away § 1 import j a v a . u t i l . with the static colour ordering and allow MCQ to colour and sort the candidate set using the candidate set, somewhat in the manner of Uroborus the serpent that eats itself. And MCSb is MCSa with an additional colour repair step. Therefore we might have an alternative hierarchy of the algorithms, such as the hierarchy of Figure 6 within the image of Uroborus.

Exact Algorithms for Maximum Clique: a brief history
We now present a brief history of complete algorithms for the maximum clique problems, starting from 1990. The algorithms are presented in chronological order.
1990: In 1990 [6] Carraghan and Pardalos present a branch and bound algorithm. Vertices are ordered in non-decreasing degree order at each depth in the binomial search with a cut-off based on the size of the largest clique found so far. Their algorithm is presented in Fortran 77 along with code to generate random graphs, consequently their empirical results are entirely reproducible. Their algorithm is similar to MC (Listing 1.1) but sorts the candidate set P using current degree in each call to expand.

1992:
In [18] Pardalos and Rodgers present a zero-one encoding of the problem where a vertex v is represented by a variable x v that takes the value 1 if search decides that v is in the clique and 0 if it is rejected. Pruning takes place via the constraint ¬adjacent(u, v) → x u + x v ≤ 1 (Rule 4). In addition, a candidate vertex adjacent to all vertices in the current clique is forced into the clique (Rule 5) and a vertices of degree too low to contribute to the growing clique is rejected (Rule 7). The branch and bound search selects variables dynamically based on current degree in the candidate set: a non-greedy selection chooses a vertex of lowest degree and greedy selects highest degree. The computational results showed that greedy was good for (easy) sparse graphs and non-greedy good for (hard) dense graphs.

1994:
In [19] Pardalos and Xue reviewed algorithms for the enumeration problem (counting maximal cliques) and exact algorithms for the maximum clique problem. Although dated, it continues to be an excellent review.

1997:
In [31] graph colouring and fractional colouring is used to bound search. Comparing again to MC (Listing 1.1) the candidate set is coloured greedily and if the size of the current clique plus the number of colours used is less than or equal to the size of the largest clique found so far that branch of search is cut off. In [31] vertices are selected in non-increasing degree order, the opposite of that proposed by [18]. We can get a similar effect to [31] in MC if we allow free selection of vertices, colour newP between lines 42 and 43 and make the recursive call to expand in line 43 conditional on the colour bound.
2002: Patric R. J.Östergård proposed an algorithm that has a dynamic programming flavour [17]. The search process starts by finding the largest clique containing vertices drawn from the set S n = {v n } and records it size in c[n]. Search then proceeds to find the largest clique in the set S i = {v i , v i+1 , ..., v n } using the value in c[i + 1] as a bound. The vertices are ordered at the top of search in colour order, i.e. the vertices are coloured greedily and then ordered in non-decreasing colour order, similar to that in numberSort Listing 1.3.Östergård's algorithm is available as Cliquer 5 . In the same year Torsten Fahle [8] presented a simple algorithm (Algorithm 1) that is essentially MC but with a free selection of vertices rather than the fixed iteration in line 36 of Listing 1.1 and dynamic maintenance of vertex degree in the candidate set. This is then enhanced (Algorithm 2) with forced accept and forced reject steps similar to Rules 4, 5 and 7 of [18] and the algorithm is named DF (Domain F iltering). DF is then enhanced to incorporate a colouring bound, similar to that in Wood [31].
2003: Jean-Charles Régin proposed a constraint programming model for the maximum clique problem [20]. His model uses a matching in a duplicated graph to deliver a bound within search, a Not Set as used in the Bron Kerbosch enumeration Algorithm 457 [4] and vertex selection using the pivoting strategy similar to that in [4,2,28,7]. That same year Tomita reported MCQ [26].
2007: Tomita proposed MCR [25] and in the same year Janez Konc and Dusanka Janezic proposed the MaxCliqueDyn algorithm [12] 6 . The algorithm is essentially MCQ [26] with dynamic reordering of vertices in the candidate set, using current degree, prior to colouring. This reordering is expensive and takes place high up in the backtrack tree and is controlled by a parameter T limit . Varying this parameter influences the cost of the search process and T limit must be tuned on an instance-by-instance basis. [23] and Tomita presented MCS [27]. In the same year Li and Quan proposed new max-SAT encodings for maximum clique [14,13].

2010: Pablo San Segundo and Cristóbal Tapia presented an early version of BBMC (BB-MCP)
2011: Pablo San Segundo proposed BBMC [22] and in [21] a version of BBMC with the colour repair steps from Tomita's MCS.
2012: Renato Carmo and Alexandre P. Züge [5] reported an empirical study of 8 algorithms including those from [6] and [8] along with MCQ, MCR (equivalent to MCQ3), MCS (equivalent to MCSb1) and MaxCliqueDyn. The claim is made that the Bron Kerbosch algorithm provides a unified framework for all the algorithms studied, although a Not Set is not used, neither do they use pivoting as described in [4,2,28,7]. Their algorithms can be viewed as an iterative (non-recursive) version of MC. Referring to algorithm MaxClique(G) in [5], the stack S holds pairs (Q,K) where Q is the currently growing clique and K is the candidate set. In line 4 a pair (Q,K) is popped from S. The while loop lines 5 to 8 selects and removes a vertex v from K (line 6), pushes the pair (Q,K) onto S (line 7), and then updates Q and K (line 8) such that Q has vertex v added to it and K becomes all vertices in K adjacent to v, i.e. K becomes the updated candidate set for updated Q. Therefore line 6 and 7 give a deferred iteration over the candidate set with v removed from K and v not added to Q. Therefore a pop in line 4 is equivalent to going once round the for loop in method expand in MC. All algorithms are coded in Python, therefore the study is objective (the authors include none of their own algorithms) and fair (all algorithms are coded by the authors and run in the same environment) and in that regard is both exceptional and laudable.

The Computational Study
The computational study attempts to answer the following questions.
1. Where does the improvement in MCS come from? By comparing MCQ with MCSa we can measure the contribution due to static colouring and by comparing MCSa with MCSb we can measure the contribution due to colour repair. 2. How much benefit can be had from the BitSet encoding? We compare MCSa with BBMC over a variety of problems. 3. We have three possible initial orderings (styles). Is any one of them better than the others and is this algorithm independent? 4. The candidate set is ordered in non-decreasing colour order. Could a tie-breaking rule influence performance? 5. Most papers use only random problems and the DIMACS benchmarks. What other problems might we use in our investigation? 6. Is it safe to recalibrate published results? Throughout our study we use a reference machine (named Cyprus), a machine with two Intel E5620 2.4GHz quad-core processors with 48 GB of memory, running linux centos 5.3 and Java version 1.6.0 07.

MCQ versus MCS: static ordering and colour repair
Is MCS faster than MCQ, and if so why? Just to recap, MCSa is MCQ with a static colour ordering set at the top of search and MCSb is MCSa with the colour repair mechanism. By comparing these algorithms we can determine if indeed MCSb is faster than MCQ and where that gain comes from: the static colouring order or the colour repair. We start our investigation with Erdós-Rënyi random graphs G(n, p) where n is the number of vertices and each edge is included in the graph with probability p independent from every other edge. The code for generating these graphs is given in Appendix 2.
The  Figure 7 shows on the left average number of nodes against edge probability and on the right average run time in milliseconds against edge probability for MCQ1, MCSa1 and MCSb1. The top row has n = 100, middle row n = 150 and bottom row n = 200 7 . As we apply the modifications to MCQ we see a reduction in nodes with MCSb1 exploring less states than MCSa1 and MCSa1 less than MCQ1. But on the right we see that reduction in search space does not always result in a reduction in run time. MCSb1 is always slower than MCSa1, i.e. the colour repair is too expensive and when n = 100 MCSb1 is often more expensive to run than MCQ! Therefore it appears that MCS gets its advantage just from the static colour ordering and that the colour repair slows it down.
We also see a region where problems are hard for all our algorithms, at n = 100 and n = 150, both in terms of nodes and run time. However at n = 200 there is a different picture. We see a hard region in terms of nodes but an ever increasing run time. That is, even though nodes are falling cpu time is climbing. This agrees with the tabulated results in [22] (Tables 4 and 5 on page 580) for BB-MaxClique. It is a conjecture that run time increases because the cost of each node (call to expand) incurs more cost in the colouring of the relatively larger candidate set. In going from G(200, 0.90) to G(200, 0.95) maximum clique size increased on average from 41 to 62, a 50% increase, and for MCSa1 the average number of nodes fell by 20% (30% for MCSb1): search space has fallen, clique size has increased, this increases the cost of colouring and this results in an overall increase in run time. Figure 7 shows erratic behaviour on G(100, p) with 0.85 ≤ p ≤ 0.90. Why is this? Figure  8 shows on the left a plot of the number of edges in each of the 100 instances generated for G(100, p) with p ∈ {0.86, 0.87, 0.88}, one contour for each. G(100, 0.87) instances are very much a mix between G(100, 0.86) and G(100, 0.88) instances. On the right is a scatter plot of search effort against edge probability for MCSa1 on G(100, p), sample size 100. We see a large variation and it is this that gives us the erratic behaviour with ≤ 0.85p ≤ 0.90. Most likely, a larger sample size would smooth out the average.
We now report on the 66 DIMACS instances [1] in Table 1. For each algorithm we have 3 entries: the number of nodes, cpu time in seconds, and in brackets the size of the largest clique found. Each algorithm was allowed 14,400 cpu seconds, i.e. 4 hours, and if that was exceeded we have a table entry of "-". The best cpu time in a row is in bold font, and when cpu time limit is exceeded the largest maximum clique size is emboldened. A time entry of 0 corresponds to a run time of less than a second and we then consider the problem as being too easy to be of interest. Overall, we see that MCQ1 is rarely the best choice with MCSa1 or MCSb1 performing better. There are 11 problems where MCSb1 beats MCSa1 and 8 problems where MCSa1 beats MCSb1. Therefore the DIMACS benchmarks don't significantly separate the behaviour of these two algorithms.

BBMC versus MCSa: a change of representation
What advantage is to be had from the change of representation between MCSa and BBMC, i.e. representing sets as ArrayList in MCSa and as a BitSet in BBMC? MCSa and BBMC are at heart the same algorithm. They both produce the same colourings, order the candidate set in the same way and explore the same backtrack tree. The only difference is in the implementation of sets and how these are exploited. Figure 9 shows on the left run time of MCSa1 (x-axis) against run time of BBMC1 (y-axis) in milliseconds on each of the G(100, p) random instances. On the right we have the same but for G (200, p). The dotted line is the reference x = y. If points are below the line then BBMC1 is faster than MCSa1. BBMC1 is typically twice as fast as MCSa1.   In Table 2 we tabulate Goldilocks instances from the DIMACS benchmark suite: we remove the instances that are too easy (take less than a second) and those that are too hard (take more than 4 hours) leaving those that are "just right" for both algorithms. Under each algorithm we have: nodes visited (and this is the same for both algorithms), run time in seconds and in brackets size of the maximum clique. The column on the far right is the ratio of MCSa1's run time over BBMC1's run time, and a value greater than 1 shows that BBMC1 was faster by that amount. Again we see similar behaviour to that observed over the random problems: BBMC1 is typically twice as fast as MCSa1.

MCQ and MCS: the effect of style
What effect does the initial ordering of vertices have on performance? First, we investigate MCQ, MCSa and MCSb with our three orderings: style 1 being non-decreasing degree, style 2 a minimum width ordering, style 3 non-decreasing degree tie-breaking on the accumulated degree of neighbours. At this stage we do not consider BBMC as it is just a BitSet encoding of MCSa.
We use random problems G(n, p) with n equal to 100 and 150 with sample size of 100. This is shown graphically in Figure 10: on the left G(100, p) and on the right G(150, p) with average nodes visited plotted against edge probability. Plots on the first row are for MCQ, middle row MCSa and bottom MCSb. For MCQ style 3 is the winner and style 2 is worst, whereas in MCSa and MCSb style 2 is always best. Why is this? In MCQ the candidate set is ordered as the result of colouring and this order is then used in the next colouring. Therefore MCQ gradually disrupts the initial minimum width ordering but MCSa and MCSb do not (and neither does BBMC). The minimum width ordering (style 2) is best for MCSa, MCSb and BBMC. Note that MCQ3 is Tomita's MCR [25] and our experiments on G(n, p) show that MCR (MCQ3) beats MCQ (MCQ1). We now report on the 66 DIMACS instances [1], Tables 3 and 4. Table 3 gives run times in seconds. An entry of "-" corresponds to the cpu time limit of 14,400 seconds being exceeded and search terminating early and an entry of 0 (zero) when run time is less than a second. For each algorithm we have three columns, one for each style: first column s 1 is style 1 with vertices in non-increasing degree order, s 2 is style 2 with vertices in minimum width order, s 3 is style 3 with vertices in non-increasing degree order tie-breaking on sum of neighbouring degrees. Table 4 is the number of nodes, in thousands, for the experiments in Table 3. In Table 3 a bold entry is the best run time for that algorithm against the problem instance, and this is done only when run  MCSb  BBMC  instance  s1  s2  s3  s1  s2  s3  s1  s2  s3  s1  s2  s3  brock200-1  7  5  4  4  3  3  3  3  3  2  1  1 Table 3. DIMACS instances: the effect of style on run time in seconds times differ significantly. For MCQ there is no particular style that is a consistent winner. This is a surprise as MCQ3 is Tomita's MCR and in [25] it is claimed that MCR was faster than MCQ. The evidence that supports this claim is Table 2 of [25], 8 of the 66 DIMACS instances. For MCSa and BBMC style 2 is best more often than not, and in MCSb style 1 is best more often than not.
Overall we see that the BBMC2 is our best algorithm, i.e. BBMC with a minimum width ordering.

MCSa: tie breaking in colour classes
It has often been reported that the order that vertices are chosen for expansion can have a profound effect on search effort. This is demonstrated in Figure 11 where the algorithm MC0 was applied to G(70, p) with a sample size of 100 using all three styles (MC1, MC2, MC3), using index order (MCi) and the reverse orderings (MC-1, MC-2, MC-3). Plotted is edge probability against logarithm of average run time. MC (and MC0) expand vertices in the candidate set from last to first therefore with a style of 1 vertices are expanded in non-decreasing degree order (smallest first). Figure  11 shows that the order can have an enormous effect. Styles -1 and -3 (largest degree first) are orders of magnitude worse than styles 1 and 3 (smallest degree first). In fact the MC-1 and MC-3 experiments were abandoned after 72 hours with G(70, 0.96) incomplete. This suggests that ordering vertices in colour classes may be worthwhile.  Fig. 11. MC applied to G(70,p) with different initial orderings.
The numberSort method used by MCS delivers a colour-ordered candidate set. Vertices are picked out of colour classes and added to the candidate set in the order they were added to the colour class. If vertices are coloured in non-increasing degree order, vertices in a colour class will also be in non-increasing degree order. Consequently the candidate set will be in non-decreasing colour order tie-breaking on non-increasing degree. The expand method iterates over the candidate set from last to first and expands vertices of the same colour class in non-increasing degree order. What might happen if this was reversed so that we visit vertices from the same colour class in non-decreasing degree order?

More Benchmarks (not DIMACS)
In [7] experiments are performed on counting maximal cliques in exceptionally large sparse graphs, such as the Pajek data sets (graphs with hundreds of thousands of vertices 8 ) and SNAP data sets (graphs with vertices in the millions 9 ). Those graphs are out of the reach of the exact algorithms reported here. The initial reason for this is space consumption. To tackle such large sparse problems we require a change of representation, away from the adjacency matrix and towards the adjacency lists as used in [7]. Therefore we explore large random instances as in [22,27] to further investigate ordering and the effect of the BitSet representation, the hard solvable instances in BHOSLIB to see how far we can go, and structured graphs produced via the SNAP (Stanford Network Analysis Project) graph generator. But first, we start with BHOSLIB.
In Table 5 we have the only instances from the BHOSLIB suite (Benchmarks with Hidden Optimum Solutions 10 ) that could be solved in 4 hours. Each instance has a maximum clique of size 30. A bold entry is the best run time. For this suite we see that with respect to style there is no clear winner. Table 6 shows results on large random problems. Similar experiments are reported in Table 4 and 5 of [22] and Table 2 in [27]. The first three columns are the nodes visited, and this is the same for MCSa and BBMC. Run times are then given in seconds for MCSa and BBMC using each of the tree styles. Highlighted in bold is the search of fewest nodes and this is style 2 (minimum width ordering) in all but one case. Comparing the run times we see that as problems get larger, involving more vertices, the relative speed difference between BBMC and MCSa diminishes and at n = 15, 000 MCSa and BBMC's performances are substantially the same.   Table 6. Large random graphs, sample size 10 The graphgen program was downloaded from the SNAP web site and modified to use a random seed so that generated graphs with the same parameters were actually different. This allows us to generate a variety of graphs, such as complete graphs, star graphs, 2D grid graphs, Erdós-Rënyi random graphs with an exact number of edges, k-regular graphs (each vertex with degree k), Albert-Barbasi graphs, power law graphs, Klienberg copying model graphs and small-world graphs. Finding maximum cliques in a complete graph, star graph and 2D grid graph is trivial. Similarly, and surprisingly, small scale experiments suggested that Albert-Barbasi and Klienberg's graphs are also easy with respect to maximum clique. However k-regular and small world are a challenge.
The SNAP graphgen program was used to generated k-regular graphs KR(n, k), i.e. random graphs with n vertices each with degree k. Graphs were generated with n = 200 and 50 ≤ k ≤ 160, with k varying in steps of 5, 20 instances at each point. BBMC1 and BBMC2 were then applied to each instance. Obviously, with style equal to 1 or 3, there is no heuristic information to be exploited at the top of search. But would a minimum width ordering, style 2, have an advantage? Figure  13 shows average search effort in nodes plotted against uniform degree k. We see that minimum width ordering does indeed have an advantage. What is also of interest is that KR(n, k) instances tend to be harder than their G(n, p) equivalents. For example, we can compare KR(200, 160) with G(200, 0.8) in Figure 7: MCSa1 took on average 1.9 million nodes for G(200, 0.8) and BBMC1 took on average 4.7 million nodes on the twenty KR(200, 160) instances.
Small-world graphs SW (n, k, p) were then generated using graphgen. This takes three parameters: n the number of vertices, k where each vertex is connected to k nearest neighbours to the right in a ring topology (i.e. vertices start with uniform degree 2k), and p is a rewiring probability. This corresponds to the graphs in Figure 1 of [29]. Small-world graphs were generated with n = 1, 000, 50 ≤ k ≤ 100 in steps of 5, 0.0 ≤ p ≤ 0.25 in steps of 0.01, 10 graphs at each point. BBMC1 was then applied to each instance to investigate how difficulty of finding a maximum clique varies with respect to k and p and also how size of maximum clique varies, i.e. this is an investigation of the problem. The results are shown as three dimensional plots in Figure 14: on the left average search effort and on the right average maximum clique size. Looking at the graph on the left: when p = 0.0 problems are easy, as p increases and randomness is introduced problems quickly get hard, but as p continues to increase the graphs tend to become predominantly random and behave more like large sparse random graphs and get easier. We also see that as neigbourhood size k increases problems get harder. We can compare the SW (1000, 100, p) to the graphs G(1000, 0.2) in Table 6 small-world instances are relatively hard. Looking at the graph on the right (average maximum clique size) we see that as rewiring probability p increases maximum cliques size decreases and as k increases so too does maximum clique size.

Calibration of results
To compare computational results across publications authors compile and run a standard C program, dfmax, against a set of benchmarks. These run times are then used as a conversion factor, and results are then taken from one publication, scaled accordingly, and then included in another publication. Recent examples of this are [17] including rescaled results from [24]; [20] including rescaled results from [17], [31] and [8]; [25] including rescaled results from [17] and [24]; [22] including rescaled results from [12]; [21] including rescaled results from [22]; [14] including rescaled results from [25] and [20]. Is this procedure safe? To test this we take two additional machines, Fais and Daleview, and calibrate them with respect to our reference machine Cyprus. We then run experiments on each machine using the Java implementations of the algorithms implemented here against some of the DIMACS benchmarks. These results are then rescaled. If the rescaling gives substantially different results from those on the reference machine this would suggest that this technique is not safe. Table 7 gives a "Rosetta Stone" for the three machines used in this study. The standard program dfmax 11 was compiled using gcc and the -O2 compiler option on each machine and then run on the benchmarks r* on each machine. Run times in seconds are tabulated for the five benchmark instances, each machine's /proc/cpuinfo is given and a conversion factor relative to the reference machine Cyprus is then computed in the same manner as that reported in [22] "... the first two graphs from the benchmark were removed (user time was considered too small) and the rest of the times averaged ...". Therefore when rescaling the run times from Fais we multiply actual run time by 0.41 and for Daleview by 0.50. machine r100.5 r200.5 r300.5 r400.5 r500.    Table 8. Calibration experiments using 3 machines, 2 algorithms and a subset of DIMACS benchmarks Table 8 shows the results of the calibration experiments. Tabulated are DIMACS benchmark instances that took more than 1 second and less than 2 minutes to solve using MCSa1 on our second slowest machine (Fais). Run times are tabulated in milliseconds (in brackets) and the actual ratio of Cyprus-time over Fais-time (expected to be 0.41) is given as well as Cyprus-time over Daleviewtime (expected to be 0.50) for each data point. Two algorithms are used, MCSa1 and BBMC1. The last row of Table 8 gives the relative performance ratios computed using the sum of the run times in the table. Referring back to Table 7 we expect a Cyprus/Fais ratio of 0.41 but empirically get 0.21 when using MCSa1 and 0.13 when using BBMC1, and expect a Cyprus/Daleview ratio of 0.50 but empirically get an average 0.23 with MCSa1 and 0.11 with BBMC1. The conversion factors in Table 7 consistently over-estimate the speed of Fais and Daleview. For example, we would expect MCSa1 applied to brock200-1 on Fais to have a run time of 19, 343 × 0.41 = 7, 930 milliseconds on Cyprus. In fact it takes 4,777 milliseconds. If we use the derived ratio in the last row of Table 8 we get 19, 343 × 0.21 = 4, 062 milliseconds, closer to actual performance on Cyprus. As another example consider san1000 using BBMC1 on Daleview. We would expect this to take 54, 816 × 0.50 = 27, 408 milliseconds on Cyprus. In fact it takes 5,927 milliseconds! If we use the conversion ratio from the last row of Table 8 we get a more accurate estimate 54, 816×0.11 = 6, 030 milliseconds.  Table 9. Calibration experiments for Cliquer and dfmax using 3 machines.
But maybe this is because we have used a C program (dfmax) to calibrate a Java program. Would we get a reliable calibration if a C program was used?Östergård's Cliquer program was downloaded and compiled on our three machines and run against DIMACS benchmarks, i.e. the experiments in Table 8 were repeated using Cliquer and dfmax with a different, and easier, set of problems. The results are shown in Table 9 12 . What we see is an actual scaling factor of 0.62 for Cliquer on Fais when dfmax predicts 0.41 and for Cliquer on Daleview 0.26 when we expect 0.50; again we see that the rescaling procedure fails. The last three columns show a dfmax calibration using problems other than the r* benchmarks and here we see an error of about 5% on Fais (expected 0.41, actual 0.39) and about 16% on Daleview (expected 0.50, actual 0.43). Therefore it appears that rescaling results using dfmax and the five r* benchmarks is not a safe procedure and can result in wrong conclusions being drawn regarding the relative performance of algorithms.

Relative algorithmic performance on different machines
But is it even safe to draw conclusions on our algorithms when we base those conclusions on experiments performed on a single machine? Previously, in Table 2 we compared MCSa against BBMC on our reference machine Cyprus and concluded that BBMC was typically twice as fast as MCSa. Will that hold on Fais and on Daleview? Table 10 takes the data from Table 8 and divides the run time of MCSa by BBMC for each instance on our three machines. On Fais BBMC is rarely more than 50% faster than MCSa and on Daleview BBMC is slower than MCSa more often than not! If experiments were performed only on Daleview using only the DIMACS instances we might draw entirely different conclusions and claim that BBMC is slower than MCSa. This change in relative algorithmic ordering has been observed on five different machines (four using the Java 1.6.0) using all of the algorithms. 13 instance  Table 10. Calibration experiment part 2, does hardware affect relative algorithmic performance? Values greater than 1 imply BBMC is faster than MCSa, less than 1 MCSa is faster.

Conclusion
We have seen that small implementation details (in MC) can result in large changes in performance. Modern programming languages with rich constructs and large libraries of utilities makes it easier for the programmer to do this. We have also drifted away from the days when algorithms were presented along with their implementation code (examples here are [4] and [18]) to presenting algorithms only in pseudo-code. Fortunately we are moving into a new era where code is being made publicly available (examples here areÖstergård's Cliquer and Konc and Janezic's MaxCliqueDyn). Hopefully this will grow and allow Computer Scientist to be better able to perform reproducible empirical studies. Tomita [27] presented MCS as an improvement on MCR brought about via two modifications: (1) a static colour ordering and (2) a colour repair step. Our study has shown that modification (1) improves performance and (2) degrades performance, i.e. MCSa is better than MCSb.
BBMC is algorithm MCSa with sets represented as bit strings, i.e. BitSet is used rather than ArrayList. Experiments on the reference machine showed a speed up typically of a factor of 2. The three styles of ordering were investigated. The orderings were quickly disrupted by MCQ, but in the other algorithms minimum width ordering was best in random problems but in the DIMACS instances there was no clear winner.
It was demonstrated in our basic algorithm MC (which does not use a colouring bound) that the order vertices were selected can have an enormous effect on search effort (and this is well known). The best order was to select vertices of low degree first, i.e. the worst-out heuristic in [19]. Incorporating this into MCSa as a tie breaker had negligible effect and the reason for this was because of the small size of colour classes in hard (dense) instances left little scope for tie-breaking.
New benchmark problems (i.e. problems rarely investigated by the maximum clique community) were investigated such as BHOSLIB, k-regular and small-world graphs. Motivation for this study was partly to compare algorithms but also to explore these problems to determine if and when they are hard.
Finally we demonstrated that the standard procedure for calibrating machines and rescaling results is unsafe, and that running our own code on different machines can lead to different relative algorithmic performance. This is disturbing. First, it suggests that to perform a fair and reliable empirical study we should not rescale other's results: we must either code up the algorithms ourselves, as done here and also by Carmo and Züge [5], or download and run code on our machines. And secondly, we should run our experiments on different machines.
All the code used in this study is available at http://www.dcs.gla.ac.uk/∼pat/maxClique