Next Article in Journal
An Efficient Algorithm for Automatic Peak Detection in Noisy Periodic and Quasi-Periodic Signals
Next Article in Special Issue
Extracting Co-Occurrence Relations from ZDDs
Previous Article in Journal
Finite Element Quadrature of Regularized Discontinuous and Singular Level Set Functions in 3D Problems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Exact Algorithms for Maximum Clique: A Computational Study

Computing Science, University of Glasgow, Glasgow G12 8QQ, UK
Algorithms 2012, 5(4), 545-587; https://doi.org/10.3390/a5040545
Submission received: 11 September 2012 / Revised: 29 October 2012 / Accepted: 29 October 2012 / Published: 19 November 2012
(This article belongs to the Special Issue Graph Algorithms)

Abstract

:
We investigate a number of recently reported exact algorithms for the maximum clique problem. The program code is presented and analyzed to show 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 effect of vertex ordering is investigated. One of the algorithms (MCS) is broken into its constituent parts and we discover that one of these parts frequently degrades performance. It is shown that the standard procedure used for rescaling published results (i.e., adjusting run times based on the calibration of a standard program over a set of benchmarks) is unsafe and can lead to incorrect conclusions being drawn from empirical data.

1. 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.

1.1. 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, where vertex u is adjacent to vertex v if and only if { u , v } is in E. 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 [1]. 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 colouring of the graph is an upper bound on the size of the maximum clique. When colouring a graph any pair of adjacent vertices are given different colours. We do not use colours but 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.

1.2. Exact Algorithms for MCP

We can address the decision and optimization problems with an exact algorithm, such as a backtracking search [2,3,4,5,6,7,8,9,10,11,12,13,14]. 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 [4,5,11,12,13,14]. 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.

1.3. Structure of the Paper

In the next section, we present in Java the following algorithms: Fahle’s Algorithm 1 [4], Tomita’s MCQ [12], MCR [15] and MCS [13] and San Segundo’s BBMC [11]. By using Java and its inheritance mechanism, algorithms are presented as modifications of previous algorithms. Three vertex orderings are then presented. 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 computational study. The study investigates MCS, determines where its speed advantage comes from, and measures the benefits resulting from the bit encoding of BBMC and the effectiveness of three vertex orderings. 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. We then conclude.

2. The Algorithms: MC, MCQ, MCR, MCS and BBMC

We start by presenting the simplest algorithm [4] 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. 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.
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 [12] 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 [15]. The presentation of MCQ is somewhat laborious but this pays off when we present two variants of MCS [13] (MCSa and MCSb) as minor changes to MCQ. BBMC [11] 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. In the code presented, we endeavour to use the same procedure names as in the original publications.
Figure 1. The hierarchy of algorithms.
Figure 1. The hierarchy of algorithms.
Algorithms 05 00545 g001

2.1. MC

MC is similar to Algorithm 1 in [4]. 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 maximum it is saved, i.e., C becomes the champion. If | C | + | P | is too small to unseat the champion search can be terminated. Otherwise the 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 n e i g h b o u r s ( v ) ), and recursing. We will call this MC.
Listing 1. The basic clique solver.
Algorithms 05 00545 i001

2.1.1. MC in Java

Listing 1 can be compared with Algorithm 1 in [4]. 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 d e g r e e where d e g r e e [ i ] is the number of vertices adjacent to vertex i (and is the sum of A [ i ] ). The variables n o d e s and c p u T i m e are used as measures of search performance, t i m e L i m i t is a bound on the run-time, m a x S i z e is the size of the largest clique found so far, s t y l e is used as a flag to customise the algorithm with respect to ordering of vertices (and is not used till we get to MCQ), and the array s o l u t i o n is the largest clique found such that s o l u t i o n [ i ] is equal to 1 if and only if vertex i is in the largest clique found.
The method s e a r c h ( ) finds a largest clique or terminates when having exceeded the allocated t i m e L i m i t . Two sets are produced: The candidate set P and the current clique C. 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 to 29). The sets P and C are represented using Java’s ArrayList, a resizable-array 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 e x p a n d . 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 n e w P is then created (line 40) where n e w P is the set of vertices in P that are adjacent to vertex v (line 41). Consequently all vertices in n e w P are adjacent to all vertices in C and all pairs of vertices in C are adjacent (i.e., C is a clique). If n e w P is empty C is maximal and if it is the largest clique found it is saved (line 42). If n e w P is not empty then C is not maximal and search can proceed via a recursive call to e x p a n d (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 e x p a n d and that is at line 37: If the combined size of the current clique and the candidate set cannot unseat the champion, 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 [4]. The method s a v e S o l u t i o n saves off the current maximal clique and records its size.

2.1.2. Observations on MC

There are several points of interest. First, there is the search process itself. If we commented out lines 37 and changed line 41 to add to n e w P all vertices in P other than v, method e x p a n d would produce the power set of P and at each depth k in the backtrack tree we would have n k calls to e x p a n d . That is, e x p a n d produces a binomial backtrack search tree of size O ( 2 n ) (see pages 6 and 7 of [17]). 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 (i.e., generating a binary backtrack tree such as in [6,8,18]). 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.
Figure 2. A binomial search tree producing the power set of { 0 , 1 , 2 , 3 } .
Figure 2. A binomial search tree producing the power set of { 0 , 1 , 2 , 3 } .
Algorithms 05 00545 g002
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. This is rewritten in class MC0 (extending MC, overwriting the e x p a n d method) Listing 2 lines 15 to 18: One line of code is replaced with 4 lines. MC0 gets the j t h element of P, calls it w (line 16 of Listing 2) and if it is adjacent to v it is added to n e w P (line 17 of Listing 2). In MC (line 41 of Listing 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 14 the new candidate set is created with a capacity of i. Why do that when we can just create n e w P with no size and let Java work it out dynamically? And why size i?
Listing 2. Inelegant but 50% faster, MC0 extends MC.
Algorithms 05 00545 i002
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 n e w P accordingly. If we do not set the size Java will give n e w P an initial size of 10 and when additions exceed this n e w P 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. 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 than it needs to. At line 37 | C | + | P | is sufficient to proceed but at line 43 it is possible that | C | + | n e w P | 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.
Finally, Listing 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 e x p a n d counts as one more node.

2.1.3. A Trace of MC

We now present two views of the MC search process over a simple problem. 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, to be read from left to right and top to bottom. 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. Also shown is the backtrack tree. The boxes correspond to calls to e x p a n d 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, therefore C represent the path from the root to the current node. A clear white box is a call to e x p a n d 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!” is a failure because | C | + | P | was too small to unseat the champion. The blue boxes correspond to calls to e x p a n d that fail first time on entering the loop at line 36 of Listing 1. By looking at the backtrack tree we get a feel for the nature of binomial search.

2.2. MCQ and MCR

We now present Tomita’s algorithm MCQ [12] as Listings 3–7. MCQ is at heart an extension of MC, performing a binomial search, with two significant advances. First, the graph induced by the candidate set is coloured using a greedy sequential colouring algorithm. This gives an upper bound on the size of the clique in P that can be added to C. Vertices in P are then selected in decreasing colour order, that is, P is ordered in non-decreasing colour order (highest colour last). And this is the second advance. Assume we select the i t h entry in P and call it v. We then know that we can colour all the vertices in P corresponding to the 0th entry up to and including the i t h entry using no more than the colour number of v. Consequently that sub-graph can contain a clique no bigger than the colour number of v, and if this is too small to unseat the largest clique, the search can be abandoned.
Figure 3. Cartoon, trace and backtrack-tree for MC on graph g10–50.
Figure 3. Cartoon, trace and backtrack-tree for MC on graph g10–50.
Algorithms 05 00545 g003

2.2.1. MCQ in Java

MCQ extends MC, Listing 3 line 3, and has an additional instance variable c o l o u r C l a s s (line 5) such that c o l o u r C l a s s [ i ] is an ArrayList of integers (line 15) and will contain all the vertices of colour i + 1 and is used when sorting vertices by their colour (lines 45 to 64, Listing 4). At the top of search (method s e a r c h , lines 12 to 21, Listing 3) vertices are sorted (call to o r d e r V e r t i c e s ( P ) at line 19) into some order, and this is described later.
Method e x p a n d (line 23 to 43 Listing 3) corresponds to the method of the same name in Figure 2 of [12]. The array c o l o u r is local to the method and holds the colour of the i t h vertex in P. The candidate set P is then sorted in non-decreasing colour order by the call to n u m b e r S o r t in line 28, and c o l o u r [ i ] is then the colour of integer vertex P . g e t ( 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 (line 30). If it is insufficient, the search terminates. Note that the loop starts at m 1 (line 29), the position of the last element in P, and counts down to zero. The i t h element of P is selected and assigned to v. As in MC we create a new candidate set n e w P , 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).
Listing 3. MCQ (part 1), Tomita 2003.
Algorithms 05 00545 i003
Listing 4. MCQ (part 1 continued), Tomita 2003.
Algorithms 05 00545 i004
Listing 5. Vertex.
Algorithms 05 00545 i005
Listing 6. MCRComparator.
Algorithms 05 00545 i006
Listing 7. MCQ (part 2), Tomita 2003.
Algorithms 05 00545 i007
Method n u m b e r S o r t (Listing 4) can be compared to the method of the same name in Figure 3 of [12]. n u m b e r S o r t takes as arguments the current clique C, an ordered A r r a y L i s t of integers C o l O r d corresponding to vertices to be coloured in that order, an A r r a y L i s t of integers P that will correspond to the coloured vertices in non-decreasing colour order, and an array of integers c o l o u r such that if v = P . g e t ( i ) (v is the i t h vertex in P) then the colour of v is c o l o u r [ i ] . Lines 45 to 64 (Listing 4) differs from Tomita’s NUMBER-SORT method because we use the additional arguments C o l O r d and the growing clique C as this allows us to easily implement our next algorithm MCS (clique C is not used in n u m b e r S o r t until we get to MCSb, therefore we carry it for convenience only.)
Rather than assigning colours to vertices explicitly, n u m b e r S o r t places vertices into colour classes, i.e., if a vertex is not adjacent to any of the vertices in c o l o u r C l a s s [ 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.
n u m b e r S o r t starts by clearing out the colour classes that might be used (line 48). In lines 49 to 55 vertices are selected from C o l O r d 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 c o n f l i c t s ). The variable c o l o u r s records the number of colours used. Lines 56 to 63 is 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 c o n f l i c t s , lines 66 to 72, takes a vertex v and an ArrayList of vertices c o l o u r C l a s s where vertices in c o l o u r C l a s s 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 c o l o u r C l a s s , 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 n u m b e r S o r t , the size of that c o l o u r C l a s s will be zero, and the for loop of lines 67 to 70 will not be performed and c o n f l i c t s returns true. The complexity of n u m b e r S o r t 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 5 and the comparator MCRComparator in Listing 6. If i is a vertex in P then the corresponding Vertex v has an integer i n d e x equal to i. The Vertex also has attributes d e g r e e and n e b D e g . d e g r e e is the degree of the vertex i n d e x and n e b D e g is the sum of the degrees of the vertices in the neighbourhood of vertex i n d e x . Given an array V of class Vertex, this can be sorted using Java’s A r r a y s . s o r t ( V ) method in O ( n . log ( n ) ) time, and is ordered by default using the c o m p a r e T o method in class Vertex. Our method forces a strict ordering of V by non-increasing degree, tie-breaking on i n d e x . This ensures reproducibility of results. If we allowed the c o m p a r e T o M e t h o d to deliver 0 when two vertices have the same degree, then A r r a y s . s o r t would break ties. If the sort method was unstable, i.e., it did not maintain the relative order of objects with equal keys [19], results may be unpredictable.
The class MCRComparator (Listing 6) allows us to sort vertices by non-increasing degree, tie breaking on the sum of the neighbourhood degree n e b D e g and then on i n d e x , giving again a strict order. This is the MCR order given in [15], 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. Given an ordered set of vertices, the width of a vertex is the number of edges that lead from that vertex to previous vertices in the order, and the width of the ordering is the maximum width of its vertices. The minimum width order (mwo) was proposed by Freuder [20] and also by Matula and Beck [21] where it was called “smallest last”, and more recently in [16] as a degeneracy ordering. The method m i n W i d t h O r d e r , lines 86 to 98 of Listing 7, sorts the array V of Vertex into an “mwo”. The vertices of V are copied into an ArrayList L (lines 87 and 89). The while loop starting at line 90 selects the vertex in L with smallest degree (lines 91 and 92) and calls it v. Vertex v is pushed onto the stack S and removed from L (line 93) and all vertices in L that are adjacent to v have their degree reduced (line 94). On termination of the while loop, vertices are popped off the stack and placed back into V, giving a minimum width (smallest last) ordering.
Method o r d e r V e r t i c e s (Listing 7 lines 74 to 84) is then called once, at the top of search. The array of Vertex V is created for sorting in lines 75 and 76, and the sum of the neighbourhood degrees is computed in lines 77 to 79. C o l O r d is then sorted in one of three orders: s t y l e = = 1 in non-increasing degree order, s t y l e = = 2 in minimum width order, s t y l e = = 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 n u m b e r S o r t and that is non-decreasing colour order. In [12] it is claimed that this is an improving order (however, no evidence was presented for this claim). In [15] Tomita proposes a new algorithm, MCR, where MCR is MCQ with a different initial ordering of vertices, i.e., MCQ with s t y l e = = 3 .

2.2.2. A Trace of MCQ

Figure 4 shows a cartoon and trace of MCQ over graph g10-50. Print statements were placed immediately after the call to e x p a n d (Listing 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 o r d e r V e r t i c e s in Listing 7. This shows for each vertex the pair < i n d e x , d e g r e e > : the first call to e x p a n d has P = { 3 , 0 , 4 , 6 , 1 , 2 , 5 , 8 , 9 , 7 } , i.e., non-decreasing degree order. MCQ makes 3 calls to e x p a n d whereas MC makes 9 calls, and the MCQ colour bound cut-off in line 30 of Listing 3 is satisfied twice (Figure 4 lines 9 and 11).
Figure 4. Trace of MCQ1 on graph g10–50.
Figure 4. Trace of MCQ1 on graph g10–50.
Algorithms 05 00545 g004

2.2.3. 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 3, | C | + | n e w P | must be greater than m a x S i z e . Since the colour of the vertex selected c o l o u r [ i ] was sufficient to satisfy the condition of line 30, it must be that integer vertex v (line 31) is adjacent to at least c o l o u r [ i ] vertices in P and thus in n e w P , therefore the next recursive call will not immediately fail. Consequently each call to e x p a n d 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 e x p a n d (lines 12 to 20).
Why is c o l o u r C l a s s 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>> Java generates an implicit cast, so nothing is gained—it is merely syntactic sugar.
Tomita’s presentation of MCQ [12] differs from Listing 3 in that it initially colours the sorted vertices prior to calling e x p a n d and thereafter colour-sorts the new candidate set immediately before making a recursive call to e x p a n d . Appendix 1 explains this in detail and investigates the effect on performance.

2.3. MCS

In [13] MCS is presented as two modifications to MCQ. The first modification is to use “... 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 using 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 [13]). When colouring vertices, an attempt is made to reduce the colours used by performing exchanges between vertices in different colour classes. In [13] 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 with only 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 [13] 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 up to colour class j and v down to colour class i.
Figure 5. A repair scenario with colour classes i, j and k.
Figure 5. A repair scenario with colour classes i, j and k.
Algorithms 05 00545 g005
Experiments were then presented in [13] 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, which I call MCSa, 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.

2.3.1. MCSa in Java

In Listing 8 we present MCSa as an extension to MCQ. Method s e a r c h creates an explicit colour ordering C o l O r d and the e x p a n d method is called with this in line 18 (compare this to line 20 of MCQ). Method e x p a n d now takes three arguments: The growing clique C, the candidate set P and the colouring order C o l O r d . In line 26 n u m b e r S o r t is called using C o l O r d (compare with line 28 in MCQ) and lines 27 to 45 are essentially the same as lines 29 to 42 in MCQ with the exception that C o l O r d must also be copied and updated (lines 32, 36 and 37) prior to the recursive call to e x p a n d (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.

2.3.2. MCSb in Java

In Listing 9 we present MCSb as an extension to MCSa: The difference between MCSb and MCSa is in n u m b e r S o r t , with the addition of lines 10 and 20. At line 10 we compute d e l t a 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 r e p a i r of lines 43 to 57. The r e p a i r 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 g e t S i n g l e C o n f l i c t V a r i a b l e , lines 32 to 41, searches for such a vertex. It takes as arguments a vertex v and a colour class c o l o u r C l a s s . If v is adjacent to only one vertex in c o l o u r C l a s s the index of that vertex is returned (line 40), where 0 c o n f l i c V a r < n , otherwise a negative number is returned (line 39). The r e p a i r 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 r e p a i r delivers t r u e (line 53). Otherwise, no repair occurred (line 56).
Listing 8. MCSa, Tomita 2010.
Algorithms 05 00545 i008
Listing 9. MCSb, Tomita 2010.
Algorithms 05 00545 i009

2.3.3. Observations on MCS

Tomita did not investigate where MCS’s improvement comes from and neither did [2], coding up MCS in Python in one piece. However San Segundo did [10], incrementally adding colour repair to BBMC. 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.

2.4. BBMC

San Segundo’s BB-MaxClique algorithm [11] (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.
  • The “BB” in “BB-MaxClique” is for “Bit Board”. Sets are represented using bit strings.
  • BBMC colours the candidate set using a static sequential ordering, the ordering set at the top of search, the same as MCSa.
  • BBMC represents the neighbourhood of a vertex and its inverse neighbourhood as bit strings, rather than using a row of an adjacency matrix and its complement.
  • 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., [12,13]) takes a vertex perspective, deciding on the colour of a vertex.

2.4.1. BBMC in Java

We implement sets using Java’s BitSet class (a vector of bits with associated methods) and 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 10, lines 5 to 7, we have an array of BitSet N for representing neighbourhoods, i n v N as the inverse neighbourhoods (the complement of N) and V an array of Vertex. N [ i ] is then a BitSet representing the neigbourhood of the i t h vertex in the array V, and i n v N [ i ] as its complement. The array V is used at the top of search for renaming vertices (and we discuss this later).
The s e a r c h method (lines 16 to 30) creates the candidate BitSet P, current clique (as a BitSet) C, and Vertex array V. The o r d e r V e r t i c e s method renames the vertices and will be discussed later. The method B B M a x C l i q u e corresponds to the procedure in Figure 3 of [11] and can be compared to the e x p a n d method in Listing 8. In a BitSet we use c a r d i n a l i t y rather than s i z e (line 35, 40 and 44). The integer array U (same name as in [11]) is essentially the colour ordered candidate set such that if v = U [ i ] then c o l o u r [ i ] corresponds to the colour given to v and c o l o u r [ i ] c o l o u r [ i + 1 ] . The method call of line 38 colours the vertices and delivers those colours in the array c o l o u r and the sorted candidate set in U. The for loop, lines 39 to 47 (again, counting down from m 1 to zero), first tests to see if the colour cut-off occurs (line 40) and if it does the method returns.
Listing 10. San Segundo’s BB-MaxClique in Java (part 1).
Algorithms 05 00545 i010
Otherwise a new candidate BitSet is created, n e w P on line 41, as a clone of P. The current vertex v is then selected (line 42) and in line 43 v is added to the growing clique C and n e w P becomes the BitSet corresponding to the vertices in the candidate BitSet that are in the neighbourhood of v. The operation n e w P . a n d ( N [ v ] ) (line 43) is equivalent to the for loop in lines 34 to 37 of Listing 3 of MCQ. 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., n e w P is not empty) a recursive call is made to B B M a x C l i q u e . Regardless, v is removed from the current candidate BitSet and the current clique (line 46) and the for loop continues.
Method B B C o l o u r (Listing 11) corresponds to the procedure of the same name in Figure 2 of [11] 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 c o l o u r as parameters. Due to the nature of Java’s BitSet the a n d operation is not functional but actually modifies bits, consequently cloning is required (line 51 Listing 11) and we take a copy of P. In line 52 c o l o u r C l a s s records the current colour class, initially zero, and i is used as a counter for adding coloured vertices into the array U. The while loop, lines 54 to 64, builds up colour classes whilst consuming vertices in c o p y P . The BitSet Q (line 56) is the candidate BitSet as we are about to start a new colour class. The while loop of lines 57 to 64 builds a colour class: The first vertex in Q is selected (line 58) and is removed from the candidate BitSet c o p y P (line 59) and BitSet Q (line 60), Q then becomes the set of vertices that are in the current candidate BitSet (Q) and in the inverse neighborhood of v (line 61), i.e., Q becomes the BitSet of vertices that can join the same colour class with v. We then add v to the ordered candidate set U (line 62), record its colour and increment our counter (line 63). When Q is exhausted (line 57) the outer while loop (line 54) starts a new colour class (lines 55 to 64).
Listing 11. San Segundo’s BB-MaxClique in Java (part 1 continued).
Algorithms 05 00545 i011
Listing 12 shows how the candidate BitSet is renamed/reordered. In fact it is not the candidate BitSet that is reordered, rather it is the description of the neighbourhood N and its inverse i n v N that is reordered. Again, as in MCQ and MCSa, a Veretx array is created (lines 69 to 73) and is sorted into one of three possible orders (lines 74 to 76). Once sorted, a bit in position i of the candidate BitSet P corresponds to the integer vertex v = V [ i ] . i n d e x . The neighbourhood and its inverse are then reordered in the loop of lines 77 to 83. For all pairs ( i , j ) , we select the corresponding vertices u and v from V (lines 79 and 80) and if they are adjacent then the j t h bit of N [ i ] is set true, otherwise false (line 81). Similarly, the inverse neighbourhood is updated in line 82. The loop could be made twice as fast by exploiting symmetries in the adjacency matrix A. In any event, this method is called once at the top of search and is generally an insignificant contribution to run time.
Listing 12. San Segundo’s BB-MaxClique in Java (part 2).
Algorithms 05 00545 i012
BBMC requires its own s a v e S o l u t i o n method (lines 86 to 90 of Listing 12) due to C being a BitSet. Again the solution is saved into the integer array s o l u t i o n and again we need to use the Vertex array V to map bits to vertices. This is done in line 88: If the i t h bit of C is true then integer vertex V [ i ] . i n d e x is in the solution. This explains why V is global to the B B M C class.

2.4.2. Observations on BBMC

In our Java implementation, we might expect a speedup if we did away with the in-built BitSet and did our own bit-string manipulations explicitly. It is also worth noting that in [10] comparisons are made with Tomita’s results in [13] by rescaling tabulated results, i.e., Tomita’s code was not actually run, but this is not unusual.

2.5. Summary of MCQ, MCR, MCS and BBMC

Putting aside the chronology [11,12,13,15], 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 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.

3. 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 [3] 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) but sorts the candidate set P using current degree in each call to e x p a n d .
1992: In [8] 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 ¬ a d j a c e n t ( 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 was good for (hard) dense graphs.
1994: In [22] 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 [14] graph colouring and fractional colouring is used to bound search. Comparing again to MC (Listing 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 [14] vertices are selected in non-increasing degree order, the opposite of that proposed by [8]. We can get a similar effect to [14] in MC if we allow free selection of vertices, colour n e w P between lines 42 and 43 and make the recursive call to e x p a n d in line 43 conditional on the colour bound.
2002: Patric R. J. Östergård proposed an algorithm that has a dynamic programming flavour [7]. 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 n u m b e r S o r t Listing 4. Östergård’s algorithm is available as Cliquer [7]. In the same year, Torsten Fahle [4] 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 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 [8] and the algorithm is named DF (Domain Filtering). DF is then enhanced to incorporate a colouring bound, similar to that in Wood [14].
2003: Jean-Charles Régin proposed a constraint programming model for the maximum clique problem [9]. 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 [23] and vertex selection using the pivoting strategy similar to that in [16,23,24,25]. That same year Tomita reported MCQ [12].
2004: Faisal N. Abu-Khzam et al. [26] presented a number of kernelization steps to reduce a graph before and during search in the vertex cover problem, where a minimum vertex cover of the complement graph is a maximal clique in the original graph. Some of the kernelization steps are similar to the pruning rules in [4,8] although Crown Reduction appears to be novel and effective.
2007: Tomita proposed MCR [15] and in the same year Janez Konc and Dus̆anka Janez̆ic̆ proposed the MaxCliqueDyn algorithm [5]. The algorithm is essentially MCQ [12] 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 l i m i t . Varying this parameter influences the cost of the search process and T l i m i t must be tuned on an instance-by-instance basis.
2010: Pablo San Segundo and Cristóbal Tapia presented an early version of BBMC (BB-MCP) [27] and Tomita presented MCS [13]. In the same year Li and Quan proposed new max-SAT encodings for maximum clique [6,18].
2011: Pablo San Segundo proposed BBMC [11] and BBMCR [10], where BBMCR includes a colour repair step. In [10] it is noted that in [13] “... the concrete contribution of recolouring is unfortunately not made explicit.” San Segundo’s colour repair, B B _ R e C o l differs from that in Listing 9 in that a single swap can occur after a double swap (as in lines 49 to 52 of Listing 9). This cannot occur in Listing 9 because r e p a i r (line 43) is called only when a new colour class k is opened for vertex v (line 20); consequently v must have been adjacent to at least one vertex in each colour class less than k and therefore c o u n t (line 34) cannot be equal to zero at line 40.
2012: Renato Carmo and Alexandre P. Züge [2] reported an empirical study of 8 algorithms including those from [3] and [4] along with MCQ, MCR, MCS 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 [16,23,24,25]. 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). BBMC is not include in the study, MCS is not broken into its constituent parts (MCSa and MCSb), and the study uses only the DIMACS benchmarks.

4. The Computational Study

The computational study attempts to answer the following questions.
  • 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.
  • How much benefit can be obtained from the BitSet encoding? We compare MCSa with BBMC over a variety of problems.
  • We have three possible initial orderings (styles). Is any one of them better than the others and is this algorithm independent?
  • Most papers use only random problems and the DIMACS benchmarks. What other problems might we use in our investigation?
  • Is it safe to recalibrate published results?
Throughout our study we use a reference machine (named Cyprus), a machine with two Intel E5620 2.4 GHz quad-core processors with 48 GB memory, running Linux CentOS 5.3 and Java version 1.6.0_07.

4.1. MCQ vs. MCS: Static Ordering and Colour Repair

Is MCS faster than MCQ, and if so, why? 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 first experiments are on random G ( n , p ) , first with n = 100 , 0 . 40 p 0 . 99 , p varying in steps of 0.01, sample size of 100, then with n = 150 , 0 . 50 p 0 . 95 , p varying in steps of 0.05, sample size of 100, and n = 200 , 0 . 55 p 0 . 95 , p varying in steps of 0.05, sample size of 100. Unless otherwise stated, all experiments are carried out on our reference machine. The algorithms MCQ, MCSa and MCSb all use s t y l e = 1 (i.e., MCQ1, MCSa1, MCSb1). Figure 6 shows on the left the average number of nodes against the edge probability and on the right the average run time in milliseconds against the edge probability, for MCQ1, MCSa1 and MCSb1. The top row has n = 100 , middle row n = 150 and bottom row n = 200 . For MCQ1 the sample size at G ( 200 , 0 . 95 ) was reduced to 28, i.e., the MCQ1-200 job was terminated after 60 hours. 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. However, on the right we see that reduction in search space does not always result in 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.
Figure 6. G ( n , p ) , sample size 100. MCQ vs. MCS, where is the win? (left) Search effort in nodes visited (i.e., decisions made by the search process); (right) run time in milliseconds.
Figure 6. G ( n , p ) , sample size 100. MCQ vs. MCS, where is the win? (left) Search effort in nodes visited (i.e., decisions made by the search process); (right) run time in milliseconds.
Algorithms 05 00545 g006
Table 1. DIMACS instances: MCQ vs. MCS, nodes, run time in seconds and (clique size).
Table 1. DIMACS instances: MCQ vs. MCS, nodes, run time in seconds and (clique size).
instanceMCQ1MCSa1MCSb1
brock200-1868,2137(21)524,7234(21)245,1463(21)
brock400-1342,473,9504,471(27)198,359,8292,888(27)142,253,3192,551(27)
brock400-2224,839,0702,923(29)145,597,9942,089(29)61,327,0561,199(29)
brock400-3194,403,0552,322(31)120,230,5131,616(31)70,263,8461,234(31)
brock400-482,056,0861,117(33)54,440,888802(33)68,252,3521,209(33)
brock800-11,247,519,247(23)1,055,945,239(23)911,465,283(21)
brock800-21,387,973,191(21)1,171,057,646(24)914,638,570(21)
brock800-31,332,309,827(21)1,159,165,900(21)914,235,793(21)
brock800-4804,901,115(26)640,444,53612,568(26)659,145,64213,924(26)
hamming10-4636,203,658(40)950,939,457(37)858,347,653(37)
johnson32-2-410,447,210,976(16)8,269,639,389(16)7,345,343,221(16)
keller5603,233,453(27)596,150,386(27)523,346,613(27)
keller6285,704,599(48)226,330,037(52)240,958,450(54)
MANN-a2738,0199(126)38,0196(126)38,5978(126)
MANN-a452,851,5724,989(345)2,851,5723,766(345)2,545,1314,118(345)
MANN-a81550,869(1100)631,141(1100)551,612(1100)
p-hat1000-1237,4372(10)176,5762(10)151,0332(10)
p-hat1000-2466,616,845(45)34,473,9781,401(46)166,655,5437,565(46)
p-hat1000-3440,569,803(52)345,925,712(55)298,537,771(56)
p-hat1500-11,642,98116(12)1,184,52614(12)990,24614(12)
p-hat1500-2414,514,960(52)231,498,292(60)259,771,137(57)
p-hat1500-3570,637,417(56)220,823,126(69)176,987,047(69)
p-hat300-33,829,00574(36)624,94713(36)713,10721(36)
p-hat500-21,022,19023(36)114,0093(36)137,5685(36)
p-hat500-3515,071,375(47)39,260,4581,381(50)104,684,0544,945(50)
p-hat700-218,968,155508(44)750,90327(44)149,052274(44)
p-hat700-3570,423,439(48)255,745,746(62)243,836,191(62)
san1000302,89520(15)150,72510(15)53,2153(15)
san200-0.9-21,149,56420(60)229,5675(60)62,7761(60)
san200-0.9-38,260,345154(44)6,815,145111(44)1,218,31732(44)
san400-0.7-155,0101(40)119,3562(40)134,7723(40)
san400-0.7-2606,15914(30)889,12519(30)754,14616(30)
san400-0.7-3582,64611(22)521,41010(22)215,7855(22)
san400-0.9-1523,531,417(56)4,536,723422(100)582,44554(100)
sanr200-0.7206,2621(18)152,8821(18)100,9771(18)
sanr200-0.944,472,276892(42)14,921,850283(42)9,730,778245(42)
sanr400-0.5380,1512(13)320,1102(13)190,7062(13)
sanr400-0.7101,213,527979(21)64,412,015711(21)46,125,168650(21)
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, and in [10] it is suggested that this behaviour is a “... phase transition to triviality ...”. 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 [11] (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 ) , the 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). The search space has fallen and the clique size has increased, which increases the cost of colouring and results in an overall increase in run time. Therefore it does not appear to be a phase transition in the sense of [28,29,30], i.e., a feature of the problem that is algorithm independent.
We now report on the 66 DIMACS instances [31] 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 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. Easy instances are not tabulated, i.e., those that took less than a second. 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 9 problems where MCSa1 beats MCSb1. Therefore, the DIMACS benchmarks do not significantly separate the behaviour of these two algorithms.
These results conflict somewhat with those in [10]. There it is claimed that colour repair, when added to BBMC, results in a performance gain in dense graphs ( p 0 . 8 ). Results are presented for a subset of the DIMACS instances with some of the difficult instances absent (brock800-*, hamming10-4, keller5, keller6, johnson32-2-4, MANN-a81, p-hat1000-3, p-hat1500-2, p-hat1500-3) and for random graphs with a sample size of 10.

4.2. BBMC vs. MCSa: A Change of Representation

What advantage is to be obtained from the change of representation between MCSa and BBMC, i.e., representing sets as A r r a y L i s t 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.
Figure 7. Run time of MCSa1 against BBMC1, on the left ( G 100 , p ) and on the right G ( 200 , p ) .
Figure 7. Run time of MCSa1 against BBMC1, on the left ( G 100 , p ) and on the right G ( 200 , p ) .
Algorithms 05 00545 g007
Figure 7 shows on the left the run time of MCSa1 (x-axis) against the run time of BBMC1 (y-axis) in milliseconds on each of the G ( 100 , p ) random instances and on the right 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 h), 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 the 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.
Table 2. DIMACS Goldilocks instances: MCSa1 vs. BBMC1, nodes, run time in seconds and clique size.
Table 2. DIMACS Goldilocks instances: MCSa1 vs. BBMC1, nodes, run time in seconds and clique size.
instanceMCSa1BBMC1MCSa1/BBMC1
brock200-1524,7234(21)524,7232(21)2.03
brock400-1198,359,8292,888(27)198,359,8291,421(27)2.03
brock400-2145,597,9942,089(29)145,597,9941,031(29)2.03
brock400-3120,230,5131,616(31)120,230,513808(31)2.00
brock400-454,440,888802(33)54,440,888394(33)2.03
brock800-4640,444,53612,568(26)640,444,5366,908(26)1.82
MANN-a2738,0196(126)38,0191(126)4.12
MANN-a452,851,5723,766(345)2,851,572542(345)6.94
p-hat1000-1176,5762(10)176,5761(10)1.80
p-hat1000-234,473,9781,401(46)34,473,978720(46)1.95
p-hat1500-11,184,52614(12)1,184,5269(12)1.52
p-hat300-3624,94713(36)624,9475(36)2.36
p-hat500-2114,0093(36)114,0091(36)2.56
p-hat500-339,260,4581,381(50)39,260,458606(50)2.28
p-hat700-2750,90327(44)750,90312(44)2.20
san1000150,72510(15)150,7255(15)1.76
san200-0.9-2229,5675(60)229,5672(60)2.36
san200-0.9-36,815,145111(44)6,815,14550(44)2.20
san400-0.7-1119,3562(40)119,3561(40)2.04
san400-0.7-2889,12519(30)889,1259(30)2.12
san400-0.7-3521,41010(22)521,4105(22)2.10
san400-0.9-14,536,723422(100)4,536,723125(100)3.37
sanr200-0.914,921,850283(42)14,921,850123(42)2.30
sanr400-0.5320,1102(13)320,1101(13)1.85
sanr400-0.764,412,015711(21)64,412,015365(21)1.95
Figure 8. The effect of style on MCQ, MCSa and MCSb. On the left G ( 100 , p ) and on the right G ( 150 , p ) . Plotted is search effort in nodes against edge probability. The top two plots are for MCQ, middle plots MCSa and bottom MCSb.
Figure 8. The effect of style on MCQ, MCSa and MCSb. On the left G ( 100 , p ) and on the right G ( 150 , p ) . Plotted is search effort in nodes against edge probability. The top two plots are for MCQ, middle plots MCSa and bottom MCSb.
Algorithms 05 00545 g008

4.3. MCQ and MCS: The Effect of Initial Ordering

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 a sample size of 100. This is shown graphically in Figure 8: 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 [15] and our experiments on G ( n , p ) show that MCR (MCQ3) beats MCQ (MCQ1).
Table 3. DIMACS instances: The effect of style on run time in seconds.
Table 3. DIMACS instances: The effect of style on run time in seconds.
MCQMCSaMCSbBBMC
instance s 1 s 2 s 3 s 1 s 2 s 3 s 1 s 2 s 3 s 1 s 2 s 3
brock200-1754433333211
brock400-14,4713,6405,6102,8881,9993,7522,5513,7482,1521,4219831,952
brock400-22,9234,5731,8242,0892,4151,2041,1992,6952,6471,0311,230616
brock400-32,3222,6961,4911,6161,4041,0271,2342,8172,117808711534
brock400-41,1175741,8728023381,2831,2091,154607394158651
brock800-1
brock800-2
brock800-39,47912,815
brock800-412,56813,50213,9246,9087,75012,992
hamming10-4
johnson32-2-4
keller5
keller6
MANN-a27999676878111
MANN-a454,9895,3694,9993,7663,5393,7334,1183,9524,242542580554
MANN-a81
p-hat1000-1221222222111
p-hat1000-21,4018611,4817,5658,4596,606720431763
p-hat1000-3
p-hat1500-11616151415151414169910
p-hat1500-2
p-hat300-37412769131012212418545
p-hat500-31,3816601,1224,9456,9825,167606282500
p-hat700-25085513532725247493108121111
p-hat700-312,2446,7545,6937,000
san1000201918101010333555
san200-0.9-2207335515111202
san200-0.9-3154459111065323850027
san400-0.7-11522174301181
san400-0.7-2144716192623169491110
san400-0.7-3113841102239513195918
san400-0.9-14228,8545401253,799
sanr200-0.7121111111000
sanr200-0.98921,7821,083283229364245227444123104164
sanr400-0.5222222222111
sanr400-0.79791,075975711608719650660674365326369
We now report on the 66 DIMACS instances [31], Table 3 and Table 4. Table 3 gives run times in seconds. An entry of “—” corresponds to the CPU time limit of 14,400 s being exceeded and the search terminating early. Problems that took less than a second have been excluded from the tables. 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. DIMACS instances: The effect of style on search nodes in 1,000,000’s.
Table 4. DIMACS instances: The effect of style on search nodes in 1,000,000’s.
MCQMCSaMCSbBBMC
instance s 1 s 2 s 3 s 1 s 2 s 3 s 1 s 2 s 3 s 1 s 2 s 3
brock200-10.860.590.5150.520.300.320.240.260.270.520.300.32
brock400-1342.5266.2455.3198.4132.8278.9142.3208.6114.8198.4132.8278.9
brock400-2224.8381.9125.2145.6178.576.461.3151.8154.3145.6178.576.4
brock400-3194.4214.0114.7120.2101.672.870.3163.5125.5120.2101.672.8
brock400-482.136.5148.354.419.390.968.362.731.954.419.390.9
brock800-1
brock800-2
brock800-3949.41,369.1
brock800-4640.4773.3659.1640.4773.31,440.8
hamming10-4
johnson32-2-4
keller5
keller6
MANN-a270.0380.0380.0380.0380.0380.0380.0380.0340.0380.0380.0380.038
MANN-a452.92.92.92.92.92.82.52.42.52.92.92.8
MANN-a81
p-hat1000-12.42.52.41.81.71.81.51.481.481.81.71.8
p-hat1000-234.519.236.9166.7177.9142.034.519.236.9
p-hat1000-3
p-hat1500-11.61.81.91.21.21.41.00.91.21.21.21.4
p-hat1500-2
p-hat300-33.87.14.00.620.490.640.710.820.640.620.490.64
p-hat500-339.316.930.9104.7152.8111.639.316.930.9
p-hat700-218.919.012.80.750.630.591.51.92.20.750.630.59
p-hat700-3216.5282.4216.5297.1
san10000.300.310.290.150.150.150.050.050.050.150.150.15
san200-0.9-21.14.32.10.230.060.240.060.050.030.230.060.23
san200-0.9-38.30.233.26.80.013.61.20.120.246.80.013.6
san400-0.7-10.060.120.090.120.660.150.130.010.050.120.660.15
san400-0.7-20.611.70.670.890.880.930.750.310.160.890.880.93
san400-0.7-30.581.92.30.520.921.90.220.550.990.520.921.9
san400-0.9-14.5220.20.580.024.5220.2
sanr200-0.70.210.290.220.150.180.160.100,120.110.150.180.16
sanr200-0.944.5101.062.214.912.520.69.78.119.014.912.520.6
sanr400-0.50.380.420.350.320.320.300.190.180.200.320.320.30
sanr400-0.7101.2106.7101.564.454.464.146.144.948.764.454.464.1
Table 4 gives the number of nodes, in millions, 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 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 [15] it is claimed that MCR was faster than MCQ. The evidence that supports this claim is Table 2 of [15], 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.

4.4. More Benchmarks (not DIMACS)

In [16] 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) and SNAP data sets (graphs with vertices in the millions) [35]. 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 [16]. Therefore we explore large random instances as in [11,13] 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. We start with BHOSLIB.
In Table 5 we have the only instances from the BHOSLIB suite (Benchmarks with Hidden Optimum Solutions [36]) 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 5. BHOSLIB using BBMC: 1,000’s of nodes and run time in seconds. Problems have 450 vertices and graph density 0.82.
Table 5. BHOSLIB using BBMC: 1,000’s of nodes and run time in seconds. Problems have 450 vertices and graph density 0.82.
instancenedgesBBMC1BBMC2BBMC3
frb30-15-145083,198292,0953,099626,8336,503361,9493,951
frb30-15-245083,151557,2525,404599,5436,136436,1104,490
frb30-15-345083,126167,1161,707265,1572,700118,4951,309
frb30-15-445083,194991,4609,663861,3918,5131,028,1299,781
frb30-15-545083,231282,7632,845674,9877,033281,1522,802
Table 6 shows results on large random problems. Similar experiments are reported in Tables 4 and 5 of [11] and Table 2 in [13]. 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 three 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 the performances of MCSa and BBMC are essentially the same. This was also observed in [10] and is expected: As problems get larger the BitSet requires more words to represent the set, and with more words the number of iterations within the BitSet increases.
Table 6. Large random graphs, sample size 10.
Table 6. Large random graphs, sample size 10.
instancenodesMCSaBBMC
np s 1 s 2 s 3 s 1 s 2 s 3 s 1 s 2 s 3
1,0000.14,5364,4724,563000000
0.239,47838,25038,838000000
0.3400,018371,360404,948444222
0.43,936,7613,780,7374,052,677403938262526
0.579,603,71275,555,47880,018,645860910859570574604
3,0000.1144,375142,719145,487333222
0.22,802,0112,723,4432,804,830383838323232
0.373,086,97871,653,88973,354,584964960978926930931
10,0000.15,351,5915,303,6155,432,812236252245212216214
15,0000.122,077,21221,751,10021,694,0361,1791,1171,0811,2491,2351,208
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.
Figure 9. k-Regular SNAP instances K R ( 200 , k ) , 130 k 160 , sample size 20.
Figure 9. k-Regular SNAP instances K R ( 200 , k ) , 130 k 160 , sample size 20.
Algorithms 05 00545 g009
The SNAP graphgen program was used to generated k-regular graphs K R ( 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 9 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 K R ( n , k ) instances tend to be harder than their G ( n , p ) equivalents. For example, we can compare K R ( 200 , 160 ) with G ( 200 , 0 . 8 ) in Figure 6: 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 K R ( 200 , 160 ) instances.
Figure 10. Small world graphs S W ( 200 , k , p ) : (upper) search effort, (lower) maximum clique size.
Figure 10. Small world graphs S W ( 200 , k , p ) : (upper) search effort, (lower) maximum clique size.
Algorithms 05 00545 g010
Small-World graphs S W ( 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 2 k ), and p a rewiring probability. This corresponds to the graphs in Figure 1 of [32]. 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 10: The graph above is average search effort and below average maximum clique size. Looking at the graph above: When p = 0 . 0 problems are easy; as p increases and randomness is introduced, the 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 neighbourhood size k increases, the problems get harder. We can compare the S W ( 1000 , 100 , p ) to the graphs G ( 1000 , 0 . 2 ) in Table 6: G ( 1000 , 0 . 2 ) took on average 39,478 nodes whereas S W ( 1000 , 100 , 0 . 01 ) took 709,347 nodes, S W ( 1000 , 100 , 0 . 08 ) took 2,702,199 nodes and S W ( 1000 , 100 , 0 . 25 ) 354,430 nodes. Clearly small-world instances are relatively hard. Looking at the graph below (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.

4.5. Calibration of Results

To compare computational results across publications a standard C program, dfmax, is compiled and run against a set of benchmarks. These run times are then used as a conversion factor, and the results are then taken from one publication, scaled accordingly, and then included in another publication. Recent examples of this are [7] including rescaled results from [33]; [9] including rescaled results from [7], [14] and [4]; [15] including rescaled results from [7] and [33]; [11] including rescaled results from [5]; [10] including rescaled results from [11]; [6] including rescaled results from [9,15]. 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. Conversion factors using dfmax on three machines: Cyprus, Fais and Daleview.
Table 7. Conversion factors using dfmax on three machines: Cyprus, Fais and Daleview.
machiner100.5r200.5r300.5r400.5r500.5Intel(R)GHzcacheJavascaling factor
Cyprus0.00.020.241.495.58Xeon(R) E56202.4012,288KB1.6.0_071
Fais0.00.080.583.5613.56XEON(TM) CPU2.40512KB1.5.0_060.41
Daleview0.00.090.533.0010.95Atom(TM) N2801.66512KB1.6.0_180.50
Table 7 gives a “Rosetta Stone” for the three machines used in this study. The standard program dfmax [37] 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 [11] (“... 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 the actual run time by 0.41 and for Daleview by 0.50.
Table 8 shows the results of the calibration experiments. Tabulated are a subset of DIMACS instances that took more than 1 s and less than 2 h 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 Daleview-time (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.12 when using MCSa1 and 0.14 when using BBMC1. We expect a Cyprus/Daleview ratio of 0.50 but empirically get an average 0.26 with MCSa1 and 0.10 with BBMC1. The conversion factors in Table 7 consistently overestimate 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 . 12 = 2 , 321 milliseconds. 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 . 10 = 5 , 481 milliseconds.
Table 8. Calibration experiments using 3 machines, 2 algorithms and a subset of DIMACS.
Table 8. Calibration experiments using 3 machines, 2 algorithms and a subset of DIMACS.
MCSa1BBMC1
instanceFaisDaleviewCyprusFaisDaleviewCyprus
brock200-10.25(19,343)0.27(17,486)1.00(4,777)0.15(15,365)0.09(25,048)1.00(2,358)
brock200-40.40(1,870)0.43(1,765)1.00(755)0.20(1,592)0.13(2,464)1.00(321)
hamming10-20.18(1,885)0.14(2,299)1.00(333)0.25(608)0.21(710)1.00(151)
hamming8-40.24(1,885)0.28(1,647)1.00(455)0.23(1,625)0.19(1,925)1.00(367)
johnson16-2-40.35(2,327)0.38(2,173)1.00(823)0.26(1,896)0.14(3,560)1.00(495)
MANN-a270.21(32,281)0.22(31,874)1.00(6,912)0.14(12,335)0.10(16,491)1.00(1,676)
p-hat1000-10.25(8,431)0.28(7,413)1.00(2,108)0.14(8,359)0.12(9,389)1.00(1,169)
p-hat1500-10.19(77,759)0.22(66,113)1.00(14,421)0.11(90,417)0.10(92,210)1.00(9,516)
p-hat300-30.25(53,408)0.26(51,019)1.00(13,486)0.14(41,669)0.09(60,118)1.00(5,711)
p-hat500-20.27(13,400)0.30(12,091)1.00(3,659)0.14(10,177)0.11(13,410)1.00(1,428)
p-hat700-10.40(1,615)0.51(1,251)1.00(641)0.29(1,169)0.24(1,422)1.00(344)
san10000.11(94,107)0.12(89,330)1.00(10,460)0.10(57,868)0.11(54,816)1.00(5,927)
san200-0.9-10.29(4,918)0.31(4,705)1.00(1,444)0.18(4,201)0.11(6,588)1.00(748)
san200-0.9-20.22(23,510)0.25(20,867)1.00(5,240)0.15(14,572)0.09(23,592)1.00(2,218)
san400-0.7-10.25(10,230)0.27(9,607)1.00(2,573)0.15(8,314)0.12(10,206)1.00(1,260)
san400-0.7-20.23(84,247)0.27(72,926)1.00(19,565)0.13(71,360)0.11(87,325)1.00(9,219)
san400-0.7-30.24(45,552)0.27(40,792)1.00(10,839)0.13(39,840)0.11(46,818)1.00(5,162)
sanr200-0.70.31(5,043)0.33(4,676)1.00(1,548)0.19(4,079)0.12(6,652)1.00(795)
sanr200-0.90.23(1,249,144)0.23(1,211,762)1.00(283,681)0.15(844,487)0.09(1,409,428)1.00(123,461)
sanr400-0.50.28(9,898)0.31(8,754)1.00(2,745)0.16(9,177)0.12(12,658)1.00(1,484)
sanr400-0.70.10(7,292,771)0.28(2,544,196)1.00(711,861)0.14(2,698,444)0.10(3,737,833)1.00(365,629)
ratio (total)0.12(9,033,624)0.26(4,202,746)1.00(1,098,326)0.14(3,937,554)0.10(5,622,663)1.00(539,439)
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 were an entry “—” was a run of dfmax that was terminated after 2 minutes. 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.

4.6. 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. The -server and -client options were also tried. The -server option sometimes gave speedups of a factor of 2, sometimes a factor of 0.5, and this can also affect relative algorithmic performance.
Table 9. Calibration experiments for Cliquer and dfmax using 3 machines.
Table 9. Calibration experiments for Cliquer and dfmax using 3 machines.
Cliquerdfmax
instanceFaisDaleviewCyprusFaisDaleviewCyprus
brock200-10.66(9,760)0.43(18,710)1.00(6,490)0.39(25,150)0.42(23,020)1.00(9,730)
brock200-40.64(690)0.47(1,190)1.00(440)0.41(1,510)0.46(1,360)1.00(620)
p-hat1000-10.62(1,750)0.36(3,020)1.00(1,090)0.41(1,680)0.45(1,540)1.00(690)
p-hat700-10.67(150)0.37(270)1.00(100)
san10000.75(120)0.30(300)1.00(90)
san200-0.7-10.48(1,750)0.20(4,220)1.00(840)
san200-0.9-20.61(18,850)0.21(53,970)1.00(11,530)
san400-0.7-30.62(6,800)0.26(16,100)1.00(4,230)
sanr200-0.70.65(2,940)0.36(5,270)1.00(1,900)0.40(5,240)0.44(4,770)1.00(2,080)
sanr400-0.50.62(1,490)0.38(2,420)1.00(930)0.41(3,550)0.47(3,080)1.00(1,460)
ratio (total)0.62(44,300)0.26(105,470)1.00(27,640)0.39(37,130)0.43(33,770)1.00(14,580)

5. Conclusions

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 make 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 [8,23]) 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 Janežič’s MaxCliqueDyn) via the web. Hopefully this will grow and allow Computer Scientist to be better able to perform reproducible empirical studies.
Tomita [13] 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.
Table 10. Calibration experiment part 2, does hardware affect relative algorithmic performance? Values greater than 1 imply BBMC is faster than MCSa, and values less than 1 imply MCSa is faster.
Table 10. Calibration experiment part 2, does hardware affect relative algorithmic performance? Values greater than 1 imply BBMC is faster than MCSa, and values less than 1 imply MCSa is faster.
instance   Fais   DaleviewCyprus
brock200-11.260.702.03
brock200-41.170.722.35
hamming10-23.103.242.21
hamming8-41.160.861.24
johnson16-2-41.230.611.66
MANN-a272.621.934.12
p-hat1000-11.010.791.80
p-hat1500-10.860.721.52
p-hat300-31.280.852.36
p-hat500-21.320.902.56
p-hat700-11.380.881.86
san10001.631.631.76
san200-0.9-11.170.711.93
san200-0.9-21.610.882.36
san400-0.7-11.230.942.04
san400-0.7-21.180.842.12
san400-0.7-31.140.872.10
sanr200-0.71.240.701.95
sanr200-0.91.480.862.30
sanr400-0.51.080.691.85
sanr400-0.72.700.681.95
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 speedup 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 the best in random problems, whereas in the DIMACS instances there was no clear winner.
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 others’ results: We must either code up the algorithms ourselves, as done here and also by Carmo and Züge [2], or download and run code on our machines. Secondly, we should run our experiments on different machines.
All the codes used in this study are available online [34] along with instructions on how to run the code, the DIMACS instances, random problem generator and runtime results.

Appendix

At the top of MCQ’s search, Tomita [12] 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 only at the start of e x p a n d . Using graph g10–50 as an example, in Figure 2 of [12] vertices would initially be selected in the order [7,9,8,5,2,1,6,4,0,3] with colours respectively [6,6,6,6,6,5,4,3,2,1], i.e., using 6 colours. Here vertices are selected in order [9,8,5,7,2,1,6,4,0,3] with colours [4,3,3,2,2,2,2,1,1,1], i.e., using 4 colours, a tighter upper bound and vertices no longer in degree order.
Figure 11. The effect of Tomita’s initial colour ordering.
Figure 11. The effect of Tomita’s initial colour ordering.
Algorithms 05 00545 g011
Listing 13 presents a Java implementation of MCQ as described in [12] but in our framework. Lines 18 to 20 give an initial colour to the sorted vertices. Method n u m b e r S o r t is now called after the selection of a vertex (line 40). A similar change was made to MCSa. Figure 11 shows the effect of the initial colour ordering, using G ( 100 , p ) on calls to e x p a n d (nodes). We see on the left that Tomita’s MCQ is marginally better than MCQ1 and on the right that MCSa1 is better than Tomita’s equivalent (and we see a similar improvement in BBMC1). In conclusion, the approach adopted here is simpler, using a single colour-ordering procedure. In MCQ1 the effect on performance is detrimental but small, and in MCSa1 (and BBMC1a) it is beneficial.

Acknowledgements

I would like to thank Pablo San Segundo, Jeremy Singer, Ciaran McCreesh and my reviewers.

Listing 13. MCQTomita.
Algorithms 05 00545 i013

References

  1. Garey, M.R.; Johnson, D.S. Computers and Intractability; W.H. Freeman and Co.: New York, NY, USA, 1979. [Google Scholar]
  2. Renato, C.; Alexandre, P. Z. Branch and bound algorithms for the maximum clique problem under a unified framework. J. Braz. Comp. Soc. 2012, 18, 137–151. [Google Scholar]
  3. Randy, C.; Panos, M.P. An exact algorithm for the maximum clique problem. Oper. Res. Lett. 1990, 9, 375–382. [Google Scholar]
  4. Torsten, F. Simple and Fast: Improving a Branch-and-Bound Algorithm for Maximum Clique. In Proceedings of the ESA 2002, LNCS 2461, Rome, Italy, 17–21 September 2002; pp. 485–498.
  5. Janez, K.; Dušanka, J. An improved branch and bound algorithm for the maximum clique problem. MATCH Commun. Math. Comput. Chem. 2007, 58, 569–590. Available online: http://www.sicmm.org/ konc/ (accessed on 12 November 2012). [Google Scholar]
  6. Chu, M.; Li, Z.Q. An Efficient Branch-and-Bound Algorithm Based on Maxsat for the Maximum Clique Problem. In Proceedings of the AAAI’10, Atlanta, GA, USA, 11–15 July 2010; pp. 128–133.
  7. Östergård, P.R.J. A fast algorithm for the maximum clique problem. Discret. Appl. Math. 2002, 120, 197–207. Available online: http://users.tkk.fi/pat/cliquer.html/ (accessed on 12 November 2012). [Google Scholar] [CrossRef]
  8. Pardalos, P.M.; Rodgers, G.P. A branch and bound algorithm for the maximum clique problem. Comput. Oper. Res. 1992, 19, 363–375. [Google Scholar] [CrossRef]
  9. Régin, J.-C. Using Constraint Programming to Solve the Maximum Clique Problem. In Proceedings CP 2003, LNCS 2833, Kinsale, Ireland, 29 September–3 October 2003; pp. 634–648.
  10. Segundo, P.S.; Matia, F.; Diego, R.-L.; Miguel, H. An improved bit parallel exact maximum clique algorithm. Optim. Lett. 2011. [Google Scholar] [CrossRef]
  11. Segundo, P.S.; Diego, R.-L.; Augustín, J. An exact bit-parallel algorithm for the maximum clique problem. Comput. Oper. Res 2011, 38, 571–581. [Google Scholar] [CrossRef]
  12. Tomita, E.; Sutani, Y.; Higashi, T.; Takahashi, S.; Wakatsuki, M. An Efficient Branch-and-Bound Algorithm for Finding a Maximum Clique. In Proceedings of the DMTCS 2003, LNCS 2731, Dijon, France, 7–12 July 2003; pp. 278–289.
  13. Tomita, E.; Sutani, Y.; Higashi, T.; Takahashi, S.; Wakatsuki, M. A Simple and Faster Branch-and-Bound Algorithm for Finding Maximum Clique. In Proceedings of the WALCOM 2010, LNCS 5942, Dhaka, Bangladesh, 10–12 February 2010; pp. 191–203.
  14. Wood, D.R. An algorithm for finding a maximum clique in a graph. Oper. Res. Lett. 1997, 21, 211–217. [Google Scholar] [CrossRef]
  15. Tomita, E.; Toshikatsu, K. An efficient branch-and-bound algorithm for finding a maximum clique and computational experiments. J. Glob. Optim. 2007, 37, 95–111. [Google Scholar] [CrossRef]
  16. David, E.; Darren, S. Listing all maximal cliques in large sparse real-world graphs. Experimental Algorithms, LNCS 6630. Comput. Sci. 2011, 6630, 364–375. [Google Scholar]
  17. Knuth, D.E. Generating all Combinations and Permutations. In The Art of Computer Programming; Pearson Education Inc.: Stoughton, MA, USA, January 2006; Volume 4, pp. 1–3. [Google Scholar]
  18. Li, C.M.; Quan, Z. Combining Graph Structure Exploitation and Propositional Reasoning for the Maximum Clique Problem. In Proceedings of the ICTAI’10, Arras, France, 27–29 October 2010; Volume 1, pp. 344–351.
  19. Bentley, J.L.; McIlroy, M.D. Engineering a sort function. Softw.-Pract. Exp. 1993, 23, 1249–1265. [Google Scholar] [CrossRef]
  20. Eugene, C.F. A sufficient condition for backtrack-free search. J. Assoc. Comput. Mach. 1982, 29, 24–32. [Google Scholar]
  21. David, W.M.; Beck, L.L. Smallest-Last ordering and clustering and graph coloring algorithms. J. Assoc. Comput. Mach. 1983, 30, 417–427. [Google Scholar]
  22. Pardalos, P.M.; Xue, J. The maximum clique problem. J. Glob. Optim. 1994, 4, 301–324. [Google Scholar] [CrossRef]
  23. Bron, C.; Kerbosch, J. Algorithm 457: Finding all cliques of an undirected graph [h]. Commun. ACM 1973, 16, 575–579. [Google Scholar] [CrossRef]
  24. Akkoyunlu, E.A. The enumeration of maximal cliques of large graphs. SIAM J. Comput. 1973, 2, 1–6. [Google Scholar] [CrossRef]
  25. Tomita, E.; Tanaka, A.; Takahashi, H. The worst-case time complexity for generating all maximal cliques and computational experiments. Theor. Comput. Sci. 2006, 363, 28–42. [Google Scholar] [CrossRef]
  26. Abu-Khzam, F.N.; Collins, R.L.; Fellows, M.R.; Langston, M.A.; Suters, W.H.; Symons, C.T. Kernelization algorithms for the vertex cover problem: Theory and experiments. In ALENEX/ANALC, New Orleans, LA, USA, 10–13 January 2004; pp. 62–69.
  27. Segundo, P.S.; Tapia, C. A New Implicit Branching Strategy for Exact Maximum Clique. In Proceedings of ICTAI’10, Arras, France, 27–29 October 2010; Volume 1, pp. 352–357.
  28. Cheeseman, P.; Kanefsky, B.; Taylor, W.M. Where the Really Hard Problems are. In Proceedings of the IJCAI’91, Sidney, Australia, 24–30 August 1991; pp. 331–337.
  29. Gent, I.P.; MacIntyre, E.; Prosser, P.; Walsh, T. The Constrainednss of Search. In Proceedings of the AAAI’96, Portland, OR, USA, 4–8 August 1996; pp. 246–252.
  30. Zweig, K.A.; Palla, G.; Vicsek, T. What makes a phase transition? Analysis of the random satisfiability problem. Physica A 2010, 389, 1501–1511. [Google Scholar]
  31. DIMACS instances. Available online: ftp://dimacs.rutgers.edu/pub/challenge/graph/benchmarks/clique (accessed on 12 November 2012).
  32. Watts, D.J.; Strogatz, S.H. Collective dynamics of small world networks. Nature 1998, 394, 440–442. [Google Scholar] [CrossRef] [PubMed]
  33. Sewell, E.C. A branch and bound algorithm for the stability number of a sparse graph. INFORMS J. Comput. 1998, 10, 438–447. [Google Scholar] [CrossRef]
  34. Prosser, P. Maximum Clique Algorithms in Java. Available online: http://www.dcs.gla.ac.uk/ pat/maxClique (accessed on 12 November 2012).
  35. Stanford Large Network Dataset Collection. Available online: http://snap.stanford.edu/data/index.html (accessed on 12 November 2012).
  36. Benchmarks with Hidden Optimum Solutions. Available online: http://www.nlsde.buaa.edu.cn/kexu/benchmarks/graph-benchmarks.htm (accessed on 12 November 2012).
  37. Dfmax. Available online: ftp://dimacs.rutgers.edu/pub/dsj/clique (accessed on 12 November 2012).

Share and Cite

MDPI and ACS Style

Prosser, P. Exact Algorithms for Maximum Clique: A Computational Study . Algorithms 2012, 5, 545-587. https://doi.org/10.3390/a5040545

AMA Style

Prosser P. Exact Algorithms for Maximum Clique: A Computational Study . Algorithms. 2012; 5(4):545-587. https://doi.org/10.3390/a5040545

Chicago/Turabian Style

Prosser, Patrick. 2012. "Exact Algorithms for Maximum Clique: A Computational Study " Algorithms 5, no. 4: 545-587. https://doi.org/10.3390/a5040545

Article Metrics

Back to TopTop