{"id":1233,"date":"2019-09-26T15:59:10","date_gmt":"2019-09-26T14:59:10","guid":{"rendered":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/?p=1233"},"modified":"2019-10-08T23:58:19","modified_gmt":"2019-10-08T22:58:19","slug":"connected-components","status":"publish","type":"post","link":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/connected-components\/","title":{"rendered":"Connected Components, Concurrently"},"content":{"rendered":"

Computing\u00a0connected components<\/a> in an undirected graph is one of the most basic graph problems. Given a graph with n<\/span> vertices and m<\/span> edges, you can find its components in linear time O(n + m)<\/span> using\u00a0depth-first<\/a> or breadth-first<\/a> search. But what if you need to go faster? In this blog post, I will describe a cool new\u00a0concurrent<\/em> algorithm for this problem, which I learned this week at\u00a0the Heidelberg Laureate Forum<\/a>\u00a0from Robert Tarjan himself. The algorithm distributes the work among\u00a0n + m<\/span>\u00a0tiny processors that work concurrently most of the time and requires\u00a0O(log n)<\/span>\u00a0global synchronisation rounds. The algorithm is remarkably simple but it’s far from obvious that it works correctly and efficiently. Happily, Tarjan and his co-author S. Cliff Liu have done all the hard proofs in their recent paper<\/a>, so we can simply take the algorithm and use it.<\/p>\n

<\/p>\n

First, let’s recap the classic solution based on the depth-first search. I’ll use my favourite graph library Alga<\/a>, so my examples will be in Haskell. Below I create an example<\/span> undirected graph and compute the number of connected components by counting trees in the depth-first search forest. The graph and the forest are shown in the figure; the edges that belong to the forest are directed to illustrate the order of graph traversal.<\/p>\n

\n
\u03bb><\/span> import<\/span> Algebra.Graph.Undirected<\/span>\r\n\u03bb><\/span> example =<\/span> edges [ (1<\/span>,6<\/span>), (2<\/span>,6<\/span>), (3<\/span>,7<\/span>), ...<\/span> ]\r\n\u03bb><\/span> length $<\/span> dfsForest $<\/span> fromUndirected example\r\n3<\/span>\r\n<\/pre>\n<\/div>\n
\"Example<\/a><\/p>\n
This approach is very simple and you should definitely use it, provided that the linear complexity O(n + m)<\/span> is fast enough for your application. But what if you need to go faster? Some time ago my colleagues and I wrote a paper<\/a> where we showed how to take advantage of concurrency by implementing the breadth-first search on an FPGA, resulting in better time complexity O(d)<\/span> where d<\/span> is the diameter of the graph. This was an excellent result for our application where graphs had a small diameter, but in general, the diameter can be as large as n<\/span>. So can we do better?<\/p>\n
Yes! Tarjan’s paper presents several faster algorithms. It also describes the algorithms in a very nice compositional manner:<\/p>\n
\"Concurrent<\/a><\/p>\n
To play with and understand these algorithms, let’s translate the above definitions to Haskell, trying to preserve their clarity and conciseness while also being explicit about details. Jumping ahead a little, here is how Algorithm P will look like in the end:<\/p>\n
\n
algorithmP<\/span> =<\/span> repeat (parentConnect >>><\/span> update >>><\/span> shortcut)\r\n<\/pre>\n<\/div>\n
Pretty close! Now, let’s define the primitives\u00a0parentConnect<\/span>\u00a0etc. in terms of the underlying computational model where a lot of tiny processing threads communicate via short messages, in rounds. In a round, threads concurrently receive some messages, then update their local states, and possibly send out new messages for the next round. We can capture the local view<\/em> of a computation round by the following type:<\/p>\n
\n
type<\/span> Local<\/span> s t i o =<\/span> s -><\/span> t -><\/span> [i] -><\/span> (s, [(t, o)])\r\n<\/pre>\n<\/div>\n
Here s<\/span> is a type of states<\/em>, t<\/span> is a type of threads<\/em>, and i<\/span> and o<\/span> are types of incoming and outgoing messages<\/em> in a round of computation. In words, given the current state of a thread and a list of incoming messages, a Local<\/span>\u00a0function returns an updated state and a list of outgoing messages, each tagged with a target thread identifier — these messages will be delivered in the next computation round.<\/p>\n
We can upload\u00a0Local<\/span> functions to our tiny processors, say on an FPGA, inject input messages into the communication network, and extract the outputs after a computation round, i.e. when all threads complete the computation. This requires specialised hardware\u00a0and non-trivial setup, so let’s find a way to simulate such computations on a big sequential machine that has enough memory to have the global view<\/em>\u00a0of a round:<\/p>\n
\n
type<\/span> Global<\/span> s t i o =<\/span> (Map<\/span> t s, [(t, i)]) -><\/span> (Map<\/span> t s, [(t, o)])\r\n<\/pre>\n<\/div>\n
The Global<\/span> function takes a map from threads to their current states and a list of all input messages in the network and returns the resulting global state: a map of new thread states and a list of newly generated output messages. Note that such global computations can be composed using function composition and form a category; the operator >>><\/span> in the above code snippet for\u00a0Algorithm P<\/span> is just the left-to-right composition defined in the standard module\u00a0Control.Category<\/span>.<\/p>\n
Converting from a Local<\/span> to a Global<\/span> view is relatively straightforward:<\/p>\n
\n
round<\/span> ::<\/span> Ord<\/span> t =><\/span> Local<\/span> s t i o -><\/span> Global<\/span> s t i o\r\nround<\/span> local (states, messages) =<\/span> collect (Map<\/span>.<\/span>mapWithKey update states)\r\n  where<\/span>\r\n    deliveries =<\/span> Map<\/span>.<\/span>fromAscList (groupSort messages)\r\n    update t s =<\/span> local s t (Map<\/span>.<\/span>findWithDefault []<\/span> t deliveries)\r\n    collect    =<\/span> runWriter .<\/span> traverse writer\r\n<\/pre>\n<\/div>\n
We first find all the\u00a0deliveries<\/span>, i.e. lists of incoming messages that should be delivered to each thread, then execute the local update<\/span> functions of each thread, sequentially (note that this is OK since threads interact only between\u00a0rounds), and finally,\u00a0collect<\/span>\u00a0outgoing messages of all threads.<\/p>\n
To express Tarjan’s algorithms, we will need the generic\u00a0repeat<\/span> function that executes a given Global<\/span> computation repeatedly until thread states stop changing.<\/p>\n
\n
repeat<\/span> ::<\/span> (Eq<\/span> s, Eq<\/span> t) =><\/span> Global<\/span> s t Void<\/span> Void<\/span> -><\/span> Global<\/span> s t Void<\/span> Void<\/span>\r\nrepeat<\/span> g (states, messages)\r\n    |<\/span> states ==<\/span> newStates =<\/span> (states, messages)\r\n    |<\/span> otherwise           =<\/span> repeat g (newStates, newMessages)\r\n  where<\/span>\r\n    (newStates, newMessages) =<\/span> g (states, messages)\r\n<\/pre>\n<\/div>\n
Note that we require that the network is quiescent<\/em> before and after a given computation, as indicated by the type of incoming and outgoing messages\u00a0i = o = Void<\/span>. Of course, the computation may be comprised of multiple rounds, which can exchange messages between each other!<\/p>\n
Now let’s use the above definitions to express Algorithms P and R from Tarjan’s paper. We’ll have n + m<\/span> threads corresponding to vertices and edges, whose states will be very minimalistic: edges will have no state at all, and every vertex will store its current parent<\/em>\u00a0i.e. the minimum vertex of the connected component to which the vertex is currently assigned.<\/p>\n
\n
type<\/span> Vertex<\/span> =<\/span> Int<\/span>\r\n\r\ndata<\/span> Thread<\/span> =<\/span> VertexThread<\/span> Vertex<\/span> |<\/span> EdgeThread<\/span> Vertex<\/span> Vertex<\/span>\r\ndata<\/span> State<\/span>  =<\/span> VertexState<\/span>  Vertex<\/span> |<\/span> EdgeState<\/span>\r\n<\/pre>\n<\/div>\n
A vertex that is its own parent is called\u00a0root<\/em>; all vertices are initially roots. If this reminds you of the disjoint-set data structure<\/a> you are on the right track!\u00a0All algorithms follow the same general idea: we start by assigning each vertex to a separate component and then use edges to inform their neighbouring vertices about other reachable vertices in the component, maintaining the invariant that the root of a component is its smallest vertex. As rounds progress, we grow a parent forest, not dissimilar to the forest produced by the depth-first search for our earlier example graph, but now taking advantage of concurrency.<\/p>\n
Let’s also create a convenient type synonym for denoting global views of computation rounds involving s = State<\/span> and t = Thread<\/span>:<\/p>\n
\n
type<\/span> (~><\/span>) i o =<\/span> Global<\/span> State<\/span> Thread<\/span> i o\r\n<\/pre>\n<\/div>\n
We can now express computation primitives, such as\u00a0connect<\/span>, where edges inform their neighbours about each other:<\/p>\n
 <\/p>\n
\n
connect<\/span> ::<\/span> Void<\/span> ~><\/span> Vertex<\/span>\r\nconnect<\/span> =<\/span> round $<\/span> \\<\/span>s t _<\/span> -><\/span> case<\/span> t of<\/span>\r\n    VertexThread<\/span> _<\/span> -><\/span> (s, []<\/span>)\r\n    EdgeThread<\/span> x y -><\/span> (s, [(VertexThread<\/span> (max x y), min x y)])\r\n<\/pre>\n<\/div>\n
The type Void\u00a0~> Vertex<\/span>\u00a0says that the round starts with no messages in the network and ends with messages carrying vertices. Vertex threads are dormant in the round, whereas edge threads generate one message each, sending a smaller vertex to the thread corresponding to the larger vertex so that the latter could update its parent. Note that we express the behaviour locally, and then use the function\u00a0round<\/span>\u00a0defined above to obtain its Global<\/span> semantics.<\/p>\n
This connect<\/span> round can be followed by the update<\/span> round, where vertex threads process incoming messages, updating their parents accordingly:<\/p>\n
\n
update<\/span> ::<\/span> Vertex<\/span> ~><\/span> Void<\/span>\r\nupdate<\/span> =<\/span> round $<\/span> \\<\/span>s _<\/span> i -><\/span> case<\/span> s of<\/span>\r\n    VertexState<\/span> p -><\/span> (VertexState<\/span> $<\/span> minimum (p :<\/span> i), []<\/span>)\r\n    EdgeState<\/span>     -><\/span> (s, []<\/span>)\r\n<\/pre>\n<\/div>\n
The primitives\u00a0connect<\/span> and\u00a0update<\/span>\u00a0are not sufficient on their own; we need another crucial ingredient — the shortcut<\/span> primitive that halves the depths of trees in the parent forest, similarly to the path compression technique<\/a> used in the disjoint-set data structure.<\/p>\n
\n
shortcut<\/span> ::<\/span> Void<\/span> ~><\/span> Void<\/span>\r\nshortcut<\/span> =<\/span> request >>><\/span> respondParent >>><\/span> update\r\n  where<\/span>\r\n    request ::<\/span> Void<\/span> ~><\/span> Thread<\/span>\r\n    request =<\/span> round $<\/span> \\<\/span>s t _<\/span> -><\/span> case<\/span> s of<\/span>\r\n        VertexState<\/span> p -><\/span> (s, [(VertexThread<\/span> p, t)])\r\n        EdgeState<\/span>     -><\/span> (s, []<\/span>)\r\n\r\nrespondParent<\/span> ::<\/span> Thread<\/span> ~><\/span> Vertex<\/span>\r\nrespondParent<\/span> =<\/span> round $<\/span> \\<\/span>s _<\/span> i -><\/span> case<\/span> s of<\/span>\r\n    VertexState<\/span> p -><\/span> (s, map (,p) i)\r\n    EdgeState<\/span>     -><\/span> (s, []<\/span>)\r\n<\/pre>\n<\/div>\n
In shortcut<\/span>, every vertex request<\/span>s the parent of its parent sending itself as a “respond-to” address. In the subsequent round, each vertex thread respond<\/span>s by sending its parent. The process completes by running the update<\/span> primitive defined above. Note how we can compose simple primitives together in a type-safe way, obtaining a non-trivial shortcut<\/span>.<\/p>\n
The last primitive that I’ll cover is parentConnect<\/span>, a variation of connect<\/span> that informs the parents of two edge vertices about each other:<\/p>\n
\n
parentConnect<\/span> ::<\/span> Void<\/span> ~><\/span> Vertex<\/span>\r\nparentConnect<\/span> =<\/span> request >>><\/span> respondParent >>><\/span> receive\r\n  where<\/span>\r\n    request ::<\/span> Void<\/span> ~><\/span> Thread<\/span>\r\n    request =<\/span> round $<\/span> \\<\/span>s t _<\/span> -><\/span> case<\/span> t of<\/span>\r\n        VertexThread<\/span> _<\/span> -><\/span> (s, []<\/span>)\r\n        EdgeThread<\/span> x y -><\/span> (s, [(VertexThread<\/span> x, t), (VertexThread<\/span> y, t)])\r\n\r\n    receive ::<\/span> Vertex<\/span> ~><\/span> Vertex<\/span>\r\n    receive =<\/span> round $<\/span> \\<\/span>s _<\/span> i -><\/span> case<\/span> s of<\/span>\r\n        VertexState<\/span> _<\/span> -><\/span> (s, []<\/span>)\r\n        EdgeState<\/span>     -><\/span> case<\/span> i of<\/span>\r\n            [x, y] -><\/span> (s, [(VertexThread<\/span> (max x y), min x y)])\r\n            _<\/span>      -><\/span> error<\/span> \"Unexpected number of responses\"<\/span>\r\n<\/pre>\n<\/div>\n
We start by request<\/span>ing parents of edge vertices, continue by sending the responses using the\u00a0respondParent<\/span>\u00a0primitive defined above, and finally receive<\/span> and handle the responses in a way similar to connect<\/span>.<\/p>\n
All the ingredients required for expressing Algorithm P are now in place:<\/p>\n
\n
algorithmP<\/span> ::<\/span> Void<\/span> ~><\/span> Void<\/span>\r\nalgorithmP<\/span> =<\/span> repeat (parentConnect >>><\/span> update >>><\/span> shortcut)\r\n<\/pre>\n<\/div>\n
Great! My favourite algorithm from the paper is Algorithm R, which is a slight variation of Algorithm P:<\/p>\n
\n
algorithmR<\/span> ::<\/span> Void<\/span> ~><\/span> Void<\/span>\r\nalgorithmR<\/span> =<\/span> repeat (parentConnect >>><\/span> rootUpdate >>><\/span> shortcut)\r\n<\/pre>\n<\/div>\n
Here we use rootUpdate<\/span> instead of update<\/span>, whose only difference is that updates to non-root vertices are ignored. This makes the growth of the parent forest monotonic<\/em>: we never exchange subtrees between trees and only graft whole trees to other trees. This greatly simplifies the\u00a0analysis of the performance of Algorithm R compared to Algorithm P. In fact, according to Tarjan, analysis of Algorithm P is still an open problem! I encourage you to read the paper<\/a>, which presents five (!) algorithms for computing connected components concurrently, and also to implement the two remaining primitives, extendedConnect<\/span> and alter<\/span>, required for expressing Algorithms E, A and RA using our little modelling framework.<\/p>\n
Let’s check if our implementation of Algorithm P works correctly on the example graph:<\/p>\n
\n
initialise<\/span> ::<\/span> Graph<\/span> Int<\/span> -><\/span> Map<\/span> Thread<\/span> State<\/span>\r\ninitialise<\/span> g =<\/span> Map<\/span>.<\/span>fromList (vs ++<\/span> es)\r\n  where<\/span>\r\n    vs =<\/span> [ (VertexThread<\/span> x, VertexState<\/span> x) |<\/span> x      <-<\/span> vertexList g ]\r\n    es =<\/span> [ (EdgeThread<\/span> x y, EdgeState<\/span>    ) |<\/span> (x, y) <-<\/span> edgeList   g ]\r\n\r\nrun<\/span> ::<\/span> Global<\/span> s t Void<\/span> Void<\/span> -><\/span> Map<\/span> t s -><\/span> Map<\/span> t s\r\nrun<\/span> g m =<\/span> fst $<\/span> g (m, []<\/span>)\r\n\r\ncomponents<\/span> ::<\/span> Map<\/span> Thread<\/span> State<\/span> -><\/span> [(Int<\/span>, [Int<\/span>])]\r\ncomponents<\/span> m =<\/span> groupSort\r\n    [ (p, x) |<\/span> (VertexThread<\/span> x, VertexState<\/span> p) <-<\/span> Map<\/span>.<\/span>toList m ]\r\n\r\n\u03bb><\/span> mapM_ print $<\/span> components $<\/span> run algorithmP $<\/span> initialise example\r\n(1<\/span>,[1<\/span>,2<\/span>,3<\/span>,4<\/span>,5<\/span>,6<\/span>,7<\/span>,8<\/span>,9<\/span>,10<\/span>,11<\/span>,13<\/span>,15<\/span>,18<\/span>,22<\/span>,23<\/span>,24<\/span>,25<\/span>,27<\/span>,29<\/span>,31<\/span>,32<\/span>,35<\/span>,38<\/span>])\r\n(12<\/span>,[12<\/span>,14<\/span>,16<\/span>,17<\/span>])\r\n(19<\/span>,[19<\/span>,20<\/span>,21<\/span>,26<\/span>,28<\/span>,30<\/span>,33<\/span>,34<\/span>,36<\/span>,37<\/span>])\r\n<\/pre>\n<\/div>\n
As expected, there are three components with roots 1<\/span>, 12<\/span> and 19<\/span>, which matches the figure at the top of the blog post.<\/p>\n
And here is an animation of how Algorithm R builds up the parent forest in a monotonic manner. Note that the edges in the parent forest are now shown as undirected (unlike in the earlier depth-first search figure) since they may change their logical direction during the execution of the algorithm (e.g. see the edge 4-9<\/span>).<\/p>\n
\"Parent<\/a><\/p>\n
That’s all for now! I hope to find time soon to try these algorithms on Tinsel<\/a>, a multi-FPGA hardware platform with thousands of processors connected by a fast network. Tinsel is being developed at the University of Cambridge as part of the POETS research project<\/a>, which I worked on while at Newcastle University. Tinsel is very cool! Have a look at it if you are into FPGAs and unconventional computing architectures.<\/p>\n
P.S.: If you haven’t heard about the Heidelberg Laureate Forum before, you can read about it here<\/a>. If you are a young researcher in computer science or mathematics, you should consider applying — it’s an amazing event where you get a chance to talk to the leaders of the two fields, as well as to young researchers from all over the world. This year I’m here as part of the HLF blogging team, writing about people I met and things I learned here, and it’s been a lot of fun too — if you’d like to help cover and promote this event, get in touch with the HLF media team<\/a>.<\/p>\n
P.P.S.: And here are a couple of links: (i) a video recording<\/a> of Tarjan’s HLF lecture and (ii) the complete source code<\/a> from this blog post.<\/p>\n","protected":false},"excerpt":{"rendered":"
Computing\u00a0connected components in an undirected graph is one of the most basic graph problems. Given a graph with n vertices and m edges, you can find its components in linear time O(n + m) using\u00a0depth-first or breadth-first search. But what if you need to go faster? In this blog post, I will describe a cool … Continue reading Connected Components, Concurrently<\/span> →<\/span><\/a><\/p>\n","protected":false},"author":1174,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[11,6,12],"tags":[16,4],"class_list":["post-1233","post","type-post","status-publish","format-standard","hentry","category-algorithms","category-coding","category-research","tag-graphs","tag-haskell"],"_links":{"self":[{"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/posts\/1233","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/users\/1174"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/comments?post=1233"}],"version-history":[{"count":39,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/posts\/1233\/revisions"}],"predecessor-version":[{"id":1278,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/posts\/1233\/revisions\/1278"}],"wp:attachment":[{"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/media?parent=1233"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/categories?post=1233"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.ncl.ac.uk\/andreymokhov\/wp-json\/wp\/v2\/tags?post=1233"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}