lazy me, BronKerbosh algorithm revisited.

A while ago I wrote about the Bron - kerbosch algorithm giving a brain dead implementation based on simple lists. Of course since the amount of maximal cliques can be exponential in the size of graph, my previous implementation is useless on big graphs (say more then 80 nodes and 150 edges). It just takes forever.

Now I give a more efficient version of the algorithm that uses a simple lazy data structure to hold solution. Each element of the data structure is a set of nodes. The first part is the definition of a lazy list and a couple of ocamlgraph data structures. Pretty standard :

module PkgV = struct
    type t = int
    let compare = Pervasives.compare
    let hash i = i
    let equal = (=)
end

module UG = Graph.Imperative.Graph.Concrete(PkgV)
module N = Graph.Oper.Neighbourhood(UG)
module O = Graph.Oper.Make(Graph.Builder.I(UG))
module S = N.Vertex_Set
module RG = Graph.Rand.I(UG)

type 'a llist = 'a cell Lazy.t
and 'a cell = LList of 'a * 'a llist | Empty

exception LListEmpty

let empty = lazy(Empty)
let push e l = lazy(LList(e,l))

let hd s =
    match Lazy.force s with
    | LList (hd, _) -> hd
    | Empty -> raise LListEmpty

let rec append s1 s2 =
    lazy begin
        match Lazy.force s1 with
        | LList (hd, tl) -> LList (hd, append tl s2)
        | Empty -> Lazy.force s2
    end

let rec iter f s =
    begin
        match Lazy.force s with
        | LList (hd, tl) -> f hd ; iter f tl
        | Empty -> ()
    end

let rec map f s =
    lazy begin
        match Lazy.force s with
        | LList (hd, tl) -> LList (f hd, map f tl)
        | Empty -> Empty
    end

let rec flatten ss =
    lazy begin
        match Lazy.force ss with
        | Empty -> Empty
        | LList (hd, tl) ->
            match Lazy.force hd with
            | LList (hd2, tl2) -> LList (hd2, flatten (lazy (LList (tl2, tl))))
            | Empty -> Lazy.force (flatten tl)
    end

This second part is the algorithm itself. The fold function creates a lazy list of sets of vertex to be check by the algorithm. The body of the function each time choose a new pivot and invoke the bronKerbosch2 function on each of them in turn. When we finally manage to build a clique we return it (push r empty) embedded in a list. The flatten (map f m) takes care of unwinding the results in a flat lazy list.

let rec bronKerbosch2 gr r p x =
  let n v = N.set_from_vertex gr v in
  let rec fold acc (p,x,s) =
    if S.is_empty s then acc
    else
      let v = S.choose s in
      let rest = S.remove v s in
      let r' = S.union r (S.singleton v) in
      let p' = S.inter p (n v) in
      let x' = S.inter x (n v) in
      let acc' = (push (r',p',x') acc) in
      fold acc' (S.remove v p, S.add v x, rest)
  in
  if (S.is_empty p) && (S.is_empty x) then (push r empty)
  else
    let s = S.union p x in
    let u = S.choose s in
    let l = fold empty (p,x,S.diff p (n u)) in
    flatten (map (fun (r',p',x') -> bronKerbosch2 gr r' p' x') l)
;;

let max_independent_sets gr =
  let cgr = O.complement gr in
  let r = S.empty in
  let p = UG.fold_vertex S.add cgr S.empty in
  let x = S.empty in
  bronKerbosch2 cgr r p x
;;

To test it we just create a random graph and we iterates thought the results. Noticed that if the algorithm was not lazy, for big graphs, the iteration will get stuck waiting for the computation to be over.

let main () =
  let v = int_of_string Sys.argv.(1) in
  let e = int_of_string Sys.argv.(2) in
  let gr = RG.graph ~loops:true ~v:v ~e:e () in
  let mis = max_independent_sets gr in
  let i = ref 0 in
  iter (fun s ->
    Printf.printf "%d -> %s\n" !i (String.concat " , " (List.map string_of_int (S.elements s)));
    incr i
  ) mis
;;
main ();;

UPDATE

Following the lead from the comment below, I modified my lazy algorithm to use a smarter pivot selection strategy. On a random graph of 50 vertex and 80 edges I get almost 50% speed up ! This is great !

$./mis.native 50 80
--------------------------
Timer Mis classic. Total time: 15.035211. 
Timer Mis pivot selection. Total time: 8.708919. 
--------------------------

This is the code :

let choose gr cands s =
  let n v = N.set_from_vertex gr v in
  fst (
    S.fold (fun u (pivot, pivot_score) ->
      let score = S.cardinal (S.inter (n u) cands) in
      if score > pivot_score then (u, score)
      else (pivot, pivot_score)
    ) s (-1, -1)
  )

let rec bronKerbosch_pivot gr clique cands nots =
  let n v = N.set_from_vertex gr v in
  let rec fold acc (cands,nots,s) =
    if S.is_empty s then acc
    else
      let pivot = choose gr cands s in
      let rest = S.remove pivot s in
      let clique' = S.add pivot clique in
      let cands' = S.inter cands (n pivot) in
      let nots' = S.inter nots (n pivot) in
      let acc' = (push (clique',cands',nots') acc) in
      fold acc' (S.remove pivot cands, S.add pivot nots, rest)
  in
  if (S.is_empty cands) && (S.is_empty nots) then (push clique empty)
  else
    let s = S.union cands nots in
    let pivot = choose gr cands s in
    let rest = S.diff cands (n pivot) in
    let l = fold empty (S.remove pivot cands,S.add pivot nots,rest) in
    flatten (map (fun (clique',cands',nots') -> bronKerbosch_pivot gr clique' cands' nots') l)

let max_independent_sets_pivot gr =
  let cgr = O.complement gr in
  let clique = S.empty in
  let cands = UG.fold_vertex S.add cgr S.empty in
  let nots = S.empty in
  bronKerbosch_pivot cgr clique cands nots