3.2.1 Greedy and Path Growing algorithms

In the following we give an overview of various efforts at designing efficient linear or close to linear time algorithms for the edge-weighted matching problem.

Perhaps the simplest matching heuristic is to visit each vertex $u$ once and unless $u$ is already incident on a matched edge, add an edge $(u,v)$ of maximum weight to the matching where $v$ is not already incident on a matched edge. The heuristic is shown in Algorithm 3.

As an example consider the graph in Figure 3.2 consisting of six vertices and nine edges with integer weights. If the graph is traversed in lexicographical order then Local Greedy will compute the matching $M=\{(a,e),(b,d)\}$ of weight 11. Note that the optimal matching is given by $M^{\ast }=\{(a,b),(c,d),(e,f)\}$ of weight 20.

Figure 3.2.

Example graph for approximate weighted matching.

While this strategy has linear run time in the size of $G$ , it does not offer any bound on the quality of the obtained solution. To see this consider an even length (number of edges) path $v_{1},v_{2},\ldots ,v_{k}$ where $w(v_{i},v_{i+1})=\unicode[STIX]{x1D716}$ if $i$ is odd, $w(v_{i},v_{i+1})=\unicode[STIX]{x1D6FF}$ if $i$ is even, and $\unicode[STIX]{x1D716}<\unicode[STIX]{x1D6FF}$ . If the Local Greedy algorithm processes the vertices in the given order, it will compute a matching with weight $k/2\cdot \unicode[STIX]{x1D716}$ , while the optimal solution has weight $k/2\cdot \unicode[STIX]{x1D6FF}$ .

One way to obtain a quality guarantee on the computed matching is to visit the edges in a predefined order. This gives rise to the classical greedy matching algorithm as shown in Algorithm 4. Here a matching is computed by considering each edge $e$ by non-increasing weights. Then if $e$ is not adjacent to an edge already in the matching it is added to the matching; otherwise it is discarded.

In the graph in Figure 3.2 Greedy will compute $M=\{(c,e),(b,d)\}$ of weight 16. The following result shows that Greedy always gives a solution of weight at least half of the optimal one.

It is clear that the Greedy algorithm runs in time $O(m\log n)$ due to the sorting of the edges. However, it is possible to compute a $1/2$ -approximation in linear time. One such is the Path Growing algorithm (Drake and Hougardy Reference Drake and Hougardy2003b ), shown in Algorithm 5. In line 9 of the algorithm, we concatenate the edge $(u,v)$ to the path  $P$ .

In its original form the algorithm starts at an arbitrary vertex $u$ and then grows a path $P$ from $u$ by traversing a heaviest edge incident on $u$ that leads to an unvisited vertex  $v$ . The process is then continued from $u$ and terminates when a vertex is reached such that there is no edge leading to an unvisited vertex. The result is a path $P=e_{0},e_{1},\ldots ,e_{k}$ . This is then repeated starting from every unvisited vertex. For each such path one selects the heavier set of the odd or even numbered edges to add to  $M$ .

In the graph in Figure 3.2 starting from vertex $a$ Path Growing would follow the path $a,e,c,d,b,f$ . The edges of this path would then be divided into sets $\{(a,e),(c,d),(b,f)\}$ (odd) of weight 15 and $\{(c,e),(b,d)\}$ (even) of weight 16. The final solution consists of the even set.

The run time of Path Growing is clearly $O(n+m)$ as the search for a path is only started once from each vertex and each edge is only considered twice (once from each of its endpoints). Note that the resulting matching might contain 1-augmentations. This is the case for a path of length four, $e_{0},e_{1},e_{2},e_{3}$ with edge weights 3, 1, 2 and 3. Then the resulting matching will consist of the edges $e_{0}$ and $e_{2}$ ; however, $e_{2}$ and $e_{3}$ form a 1-augmentation. Thus it is not possible to use Lemma 3.1 to show the performance guarantee. Still, it is not hard to prove that the resulting matching is a $1/2$ -approximation.

Proof. Let $K=P_{0},P_{1},\ldots ,P_{k}$ be the non-empty paths generated by Path Growing. We show that $\sum _{i}w(P_{i})\geq w(M^{\ast })$ . The result then follows since

$$\begin{eqnarray}w(M)\geq \mathop{\sum }_{i}\displaystyle \frac{1}{2}\cdot w(P_{i})\geq \displaystyle \frac{1}{2}\cdot w(M^{\ast }).\end{eqnarray}$$

Consider an edge $(u,v)\in M^{\ast }$ . Then if $(u,v)$ is not part of some $P_{i}$ at least one of $u$ and $v$ must be part of some path  $P_{j}$ , otherwise Path Growing would have started a new path from either $u$ or  $v$ . Assume therefore that $u$ is part of $P_{i}$ and that if $v$ is part of $P_{j}$ then $i<j$ . This ensures that $u$ is visited before $v$ irrespective of if $v$ is part of some $P_{j}$ or not. Thus when Path Growing visits $u$ it is possible to choose the edge $(u,v)$ . Since this is not done it must choose an edge $(u,x)$ to add to $P_{i}$ such that $w(u,x)\geq w(u,v)$ . Note that $u$ can only be incident on at most one edge in $M^{\ast }$ . It therefore follows that for every edge $e\in M^{\ast }$ either $e$ is in some $P_{i}$ or there is a unique edge $e^{\prime }\in P_{i}$ such that $w(e)\leq w(e^{\prime })$ .◻

The Path Growing algorithm can be improved further without increasing its asymptotic run time, by noting that one can compute an optimal matching for each $P_{i}$ instead of just choosing between the even and odd numbered edges. This is done using dynamic programming as follows. Let $P=e_{0},e_{1},\ldots ,e_{k}$ be the edges on a selected path and let $\operatorname{Opt}(i)$ denote the weight of an optimal matching on $e_{0},e_{1},\ldots ,e_{i-1}$ . Then the weight of an optimal matching on $P$ can be computed in linear time using the recursion

$$\begin{eqnarray}\operatorname{Opt}(i)=\max \{\operatorname{Opt}(i-1),\operatorname{Opt}(i-2)+w(e_{i-1})\},\end{eqnarray}$$

with $\operatorname{Opt}(0)=0$ and $\operatorname{Opt}(1)=w(e_{0})$ . The actual matching can be calculated by storing a Boolean value for each $i$ to indicate if $e_{i-1}$ is used in the computation of $\operatorname{Opt}(i)$ . We denote this algorithm as Path Growing-DP.

The first linear time $1/2$ -approximation algorithm for weighted matching was given by Preis (Reference Preis1999), and it is based on the notion of a dominating edge. An edge $e$ is dominating if it is heavier than all of its incident edges. This Locally Dominant Edge algorithm repeatedly locates and adds dominant edges to the matching. When an edge is added to the matching, all adjoining edges are removed from the graph before the next dominant edge is found. The outline of this strategy is shown in Algorithm 6.

To be able to break ties we rank edges of equal weight first by the ID of the higher-numbered endpoint and then by the lower-numbered one. In practice we will only be comparing edges of equal weight incident on the same vertex and thus we always break ties by comparing the ID of the opposing vertex.

To achieve the linear time run time Algorithm 6 performs a carefully designed depth-first search (DFS) of $G$ , locating and removing dominating edges as it progresses. Since the actual algorithm is quite involved, we do not list the details. Moreover, a DFS search is inherently linear, thus prohibiting parallel execution. Hence we present variants of Algorithm 6 that have worse theoretical run times, but which are faster in practice and that also lend themselves better to parallel execution.

The simplest algorithm based on discovering dominating edges is the Local Max algorithm of Birn et al. (Reference Birn, Osipov, Sanders, Schulz and Sitchinava2013). The algorithm operates in stages. At each stage a maximal set of dominating edges are discovered and added to the matching. These edges, along with their vertices, are then removed from the graph before the process is repeated. The algorithm terminates when the graph is empty. It is not hard to see that the run time of each stage is linear in the size of the remaining graph. It is possible that only a constant number of vertices and edges are removed at each stage, thus leading to a worst-case run time of $O(nm)$ . However, it can be shown that if the edge weights are drawn from a uniform random distribution, then half of the remaining edges will be removed at each iteration (Birn et al. Reference Birn, Osipov, Sanders, Schulz and Sitchinava2013). Thus in this case the expected number of rounds is $O(\log n)$ and the expected run time is $O(m)$ .

3.2.2 The Suitor algorithm

The Suitor algorithm, designed by Manne and Halappanavar (Reference Manne and Halappanavar2014), is another variant of Algorithm 6 that also computes the same matching as Greedy. The algorithm is based on considering the vertices as players that are making bids to match with one of their neighbours. The value of a bid from a vertex $u$ to a vertex $v$ is $w(u,v)$ . We denote the highest bid a vertex $u$ has received as its suitor value, denoted by  $s(u)$ .

The algorithm works by each vertex $u$ proposing to its neighbour $v$ such that $w(u,v)$ is maximized under the constraint that $v$ has not already received a higher offer than $w(u,v)$ . (We say that a vertex $v$ that satisfies this property is an eligible neighbour of  $u$ .) If at a later stage in the algorithm $v$ receives a better offer, then the bid from $u$ to $v$ is annulled, and $u$ has to make a new bid to its next heaviest eligible neighbour. The algorithm terminates when every vertex $u$ is either a suitor of one of its neighbours or when $u$ has no more vertices to propose to. At this point the edges where $s(u)=v$ and $s(v)=u$ constitute the matching.

There is no restriction on the order in which the vertices in $V$ are processed nor on the order in which rejected suitors are processed. This leaves considerable freedom in how the algorithm is implemented. Two of the most natural ways to handle this is to use either a queue or a stack to organize the vertices that still need to be processed. If a queue is used then all vertices are initially put in the queue and the next considered vertex is always taken as the head of the queue. Any vertex that gets dislodged as a suitor because its chosen neighbour receives a better offer will be added at the end of the queue. The algorithm terminates when the queue is empty.

When using a stack all vertices are initially put on the stack and the next vertex to process is taken from the top of the stack. Any vertex that gets dislodged as a suitor is put on the top of the stack. Note that this is equivalent to processing each vertex exactly once and then immediately starting to process any vertex that gets dislodged. This variant of the Suitor algorithm is shown in Algorithm 7. In the algorithm, if a vertex $x$ does not have a suitor, i.e. $s(x)=\operatorname{NULL}$ , then $w(u,s(x))=0$ .

As to correctness consider any initial dominating edge $(u,v)$ , for instance the heaviest edge in $E$ . Then the first time $u$ is processed, the value of $s(v)$ will be set to $u$ . This follows since $u$ has no better alternative and it cannot be prevented from doing so by any other neighbour $x\in N(v)$ as $w(u,v)>w(v,x)$ . Using the same argument it follows that $s(u)$ will be set to $v$ when $v$ is initially processed. The suitor values of $u$ and $v$ will remain unchanged throughout the execution for the rest of the algorithm, again because $(u,v)$ is dominant. It follows that for the final result of the vertices in $V-\{u,v\}$ , the net effect is the same as if $u$ and $v$ were removed initially along with their adjacent edges. The only difference is that some vertices might become temporary suitors of $u$ or $v$ , but when they are dislodged, as they will be, they will behave as if $u$ and $v$ had been removed initially. Thus Algorithm 7 follows the same pattern as Algorithm 6 in computing the same matching as Greedy.

For the time complexity note first that a vertex can only propose once to each of its neighbours. Thus it follows that there can at most be $2m$ proposals. The time complexity therefore depends on how the next vertex to propose to is determined. If a linear search is performed each time, the time complexity is $O(n+m\unicode[STIX]{x1D6E5})$ . If the weights of the edges incident on each vertex are sorted initially, then the rest of the algorithm runs in $O(n+m)$ time. The run time will then be dominated by the sorting, which takes time $O(m\log \unicode[STIX]{x1D6E5})$ . Other strategies include using a partial Quicksort partitioning to find the heaviest neighbours, and then only when these do not contain any more potential candidates to match with are more candidates found.

In order to say something about the expected run time of Suitor we note that there is a strong relationship between computing a greedy matching and the stable marriage problem (SM).

In the SM problem we are given two equal-sized sets $L=\{l_{1},l_{2},\ldots ,l_{n}\}$ and $R=\{r_{1},r_{2},\ldots ,r_{n}\}$ , typically referred to as ‘men’ and ‘women’ respectively. Every man and woman has a total ranking of all the members of the opposite sex. These give the ‘desirability’ for each participant to match with a member of the other set. The object is to find a complete matching $M$ (i.e. a pairing) between the entries in $L$ and $R$ such that no two $l_{i}\in L$ and $r_{j}\in R$ would both obtain a higher-ranked partner if they were to abandon their current partner in $M$ and rematch with each other. Any such solution is stable.

The main differences between computing a matching in a bipartite graph and an SM instance is that for SM one only considers how attractive two partners are relative to each other and also that attractiveness might not be symmetric. Thus one can have cycles such as $l_{1},r_{1},l_{2},r_{2},l_{1}$ where each one prefers to match with the next one in the list. This does not happen in the matching problem. Moreover, the overall objective is to find any stable solution, whereas in a matching instance possible solutions can be ranked.

Gale and Shapley (Reference Gale and Shapley1962) formulated the stable marriage problem and also proposed the first algorithm for solving it. The algorithm operates in rounds as follows. In the first round each man in $L$ proposes to his most preferred woman in  $R$ . Each woman will then reject all proposals she has received in this round except the one that is the highest in her ranking. In subsequent rounds each man that was rejected in the previous round will again propose to the woman that he has ranked highest, but now disregarding any woman that he has already proposed to in previous rounds. Gale and Shapley showed that this process will terminate with each man in $L$ being matched to a woman in $R$ and that this solution is stable. The algorithm also converges to a stable solution even if each participant has only ranked a subset of the opposing participants. For a recent discussion of stable matching and generalizations such as stable room-mates, stable fixtures, etc., see the book by Manlove (Reference Manlove2013).

Manne, Naim, Lerring and Halappanavar (Reference Manne, Naim, Lerring and Halappanavar2016) show that a greedy matching on a graph $G$ can be computed by reformulating the problem as an appropriately constructed instance of SM and then solving this using the Gale–Shapley algorithm. In fact, the Gale–Shapley algorithm has a strong resemblance to the Suitor algorithm when using a queue. Similarly, there is a variation of the Gale–Shapley algorithm due to McVitie and Wilson (Reference McVitie and Wilson1971) corresponding to the Suitor algorithm when using a stack.

Wilson (Reference Wilson1972) showed that for any profile of women’s preferences, if the men’s preferences are random, then the expected sum of men’s rankings of their mates as assigned by the Gale–Shapley algorithm is bounded above by $nH_{n}$ , where $H_{n}$ is the $n$ th harmonic number. Knoblauch (Reference Knoblauch2007) showed that this is also an approximate lower bound in the sense that the ratio of the expected sum of men’s rankings of their assigned mates and $(n+1)H_{n}-n$ has limit 1 as $n\rightarrow \infty$ . Thus if the men’s preferences are random then this sum is $\unicode[STIX]{x1D6E9}(n\ln n)$ for large  $n$ . However, it is not hard to design instances where this sum is $\unicode[STIX]{x1D6E9}(n^{2})$ . One such case is when the men have identical preferences. These results carry over to the Suitor algorithm in that for a graph with random weights the expected number of proposals made is also bounded by $O(n\ln n)$ , and if the graph is sufficiently dense by $\unicode[STIX]{x1D6E9}(n\ln n)$ .

The depth of a parallel algorithm is the length of the critical path in the algorithm, and the work in the algorithm is the total number of operations performed by the algorithm. Blelloch, Fineman and Shun (Reference Blelloch, Fineman and Shun2012) have shown that a maximal matching in a graph can be computed with $O(m)$ work and $O(\log m\log \unicode[STIX]{x1D6E5})$ depth. Their proof technique was adapted by Khan, Pothen and Ferdous (Reference Khan, Pothen and Ferdous2018b ) to show that the parallel Suitor algorithm has $O(m)$ work and $O(\log m\log \unicode[STIX]{x1D6E5})$ depth when the edge weights are chosen uniformly at random.

Among the algorithms presented thus far, experimental evaluations have shown that the highest-weight matching is obtained by Path Growing-DP. This is due to the use of dynamic programming to select an optimal solution on the discovered paths. The Global Paths Algorithm by Maue and Sanders (Reference Maue and Sanders2007) further expands this idea, by trying to get a heavier set of edges on which to perform the dynamic programming. The algorithm starts by sorting the edges and then considers each edge by order of non-increasing weight, similar to the Greedy algorithm. An edge is kept for further consideration in a set $S$ if it either connects at most two paths already contained in  $S$ , or if it completes a cycle of even length contained in  $S$ . The algorithm then uses dynamic programming on the paths and cycles in $S$ to obtain the final matching. Even though the Global Paths Algorithm in most instances gives higher weight matchings compared to Path Growing-DP, it does not improve the approximation ratio beyond half.

The idea of computing optimal matchings on restricted subgraphs was also used by Manne and Halappanavar (Reference Manne and Halappanavar2014) in the M1M2 algorithm. This algorithm is based on the observation that the union of the edges from two matchings always consists of paths and even length cycles. The M1M2 algorithm first computes a greedy matching $M_{1}$ on a graph $G(V,E,w)$ . It then computes a second greedy matching $M_{2}$ on $G(V,E\setminus M_{1},w)$ before performing dynamic programming on $M_{1}\cup M_{2}$ in the same way as was done in Global Paths Algorithm. This gave results of a similar quality to Global Paths Algorithm. Further use of this idea was shown in Idelberger and Manne (Reference Idelberger and Manne2014) where the dynamic programming was expanded to maximum weight spanning trees, while also performing mergers of multiple matchings following a tree-like structure.

While the algorithms presented thus far tend to perform quite well in practice, often producing optimal or close to optimal solutions, they all have in common that they cannot guarantee a performance ratio higher than $1/2$ . To do so requires that one considers augmentations containing more than one unmatched edge.

Drake and Hougardy (Reference Drake and Hougardy2005) presented the first such algorithm which computes a $(2/3-\unicode[STIX]{x1D716})$ -approximation in time $O(m\unicode[STIX]{x1D716}^{-1})$ . The algorithm is based on increasing the weight of a maximal matching by repeatedly performing 2-augmentations. This algorithm was subsequently simplified by Pettie and Sanders (Reference Pettie and Sanders2004) who also improved the run time to $m\log \unicode[STIX]{x1D716}^{-1}$ . Their algorithm is based on finding the best 2-augmentation centred in a node  $u$ . For a matching $M$ and 2-augmenting path $P$ , a centre node $u\in P$ has the property that each edge $(x,y)\in P\setminus M$ is either incident on $u$ or on the matching partner $v$ of  $u$ . For a particular vertex $u$ one can then find the best 2-augmenting path centred in $u$ in time $O(\deg (u)+\deg (v))$ .

Pettie and Sanders presented both a simple randomized algorithm RAMA and a slightly more involved deterministic one. In RAMA one picks a random vertex $u$ and then finds the best 2-augmentation centred at $v$ and augments if this gives a positive gain. By repeating this $(1/3)\log \unicode[STIX]{x1D716}^{-1}$ times the algorithm has an expected run time of $O(m\log \unicode[STIX]{x1D716}^{-1})$ and an expected performance ratio of $2/3-\unicode[STIX]{x1D716}$ (Pettie and Sanders Reference Pettie and Sanders2004). The algorithm can either start with an empty matching or from an existing one. A variant of RAMA, named ROMA was later determined to be more effective in practice (Maue and Sanders Reference Maue and Sanders2007). In ROMA the algorithm iterates multiple times through all vertices in a random order and performs the same augmentation step as in RAMA. This algorithm computed a higher weight matching when initialized with a matching given by Global Paths Algorithm, albeit at the cost of higher run times. Subsequent developments have given $(3/4-\unicode[STIX]{x1D716})$ -approximation algorithms running in time $O(m\log n\log \unicode[STIX]{x1D716}^{-1})$ (Duan and Pettie Reference Duan and Pettie2010, Hanke and Hougardy Reference Hanke and Hougardy2010). However, these have not been tested in practice.

Finally, Duan and Pettie (Reference Duan and Pettie2014) presented a $(1-\unicode[STIX]{x1D716})$ -approximation algorithm that runs in time $O(m\unicode[STIX]{x1D716}^{-1}\log \unicode[STIX]{x1D716}^{-1})$ . This algorithm employs the scaling approach to computing weighted matchings, with the subproblem at each scale solved by a primal–dual linear programming formulation of matching. The feasibility and complementary slackness conditions are enforced approximately, with the violation of these conditions dynamically dependent on the scale of the computation.

Table 3.1.

Structural properties of graphs used for weighted matching, sorted in increasing order of edges.

3.3.1 Comparison of exact and approximation algorithms

We begin by comparing the performance of exact algorithms for maximum edge-weighted matching with several approximation algorithms. We report on the exact algorithm implemented in LEDA (Mehlhorn and Näher Reference Mehlhorn and Näher1999). We report on three variant algorithms: LEDA1 uses no initialization, LEDA2 employs a greedy $1/2$ -approximation matching, and LEDA3 uses a fractional matching initialization. The latter initialization computes a fractional $\{0,1/2,1\}$ solution to a linear programming formulation of the edge-weighted matching problem; the solution is computed by a combinatorial algorithm, not an LP solver. The fractional solution is then rounded to obtain an initial matching. The $(2/3-\unicode[STIX]{x1D716})$ -approximation algorithm ROMA is initialized with the Global Paths Algorithm and the Suitor algorithms. We report results for the RAMA algorithm without any initialization since the initialized versions computed lower weights than the corresponding variants of the ROMA algorithm; however, it is faster than the latter. The Duan and Pettie $(1-\unicode[STIX]{x1D716})$ -approximation algorithm for weighted matching was implemented by Al-Herz and Pothen (Reference Al-Herz and Pothen2019), and we report results for it as well.

For the experiments we used an Intel Xeon E5-2660 processor-based system (part of the Purdue University Community Cluster), called Rice.Footnote ¹ The machine consists of two processors, each with ten cores running at 2.6 GHz (20 cores in total) with 25 MB unified L3 cache and 64 GB of memory. The operating system is Red Hat Enterprise Linux release 6.9. All code was developed using C++ and compiled using the g++ compiler (version: 4.4.7) using the –O3 flag.

Table 3.2.

Comparing the run times of exact and approximation algorithms for maximum edge-weighted matching. Edge weights are from a uniform random distribution in the range $[1\,1000]$ . Run times of LEDA without any initialization are reported in LEDA1, and for all other algorithms the ratios of the run time of LEDA1 to the other algorithm are presented. For the last problem in the set, LEDA1 did not complete in four hours, and LEDA with greedy initialization (LEDA2) is the baseline. LEDA3 is the LEDA matching algorithm with fractional matching initialization.

Table 3.3.

Comparison of weights obtained by the exact and approximation algorithms for maximum edge-weighted matching. For each algorithm, we report the difference between one and the ratio of the weight computed by the approximation algorithm and the maximum weight obtained by the exact algorithm, expressed as a percentage.

Table 3.4.

Comparison of cardinalities obtained by the maximum edge-weighted matching algorithm and the approximation algorithms. Note that this is not the cardinality obtained by a maximum cardinality matching. For each algorithm, we report the difference between one and the ratio of the cardinality of the matching computed by the approximation algorithm and the cardinality of the maximum edge-weighted matching, expressed as a percentage.

Run times and relative performances are reported in Table 3.2. Note that LEDA1 does not terminate on the last problem in the set, whereas LEDA2 and LEDA3 do. For the other nine problems, there is not much difference between the uninitialized version and the greedy initialization, whereas the fractional matching initialization is more effective, about four times faster than the uninitialized version. The $(2/3-\unicode[STIX]{x1D716})$ -algorithms and the $(1-\unicode[STIX]{x1D716})$ -approximation algorithms are 20–30 times faster than LEDA1, while the Suitor algorithm is about 1700 times faster, in geometric mean. For the last problem, we use LEDA2 as the baseline algorithm, the Suitor algorithm is 2000 times faster than LEDA2, while for the other approximation algorithms the factor is in the range 35–50.

Now we compare the weights computed by the algorithms. All the exact algorithmic variants compute the same weight, and this is reported in Table 3.3. For the approximation algorithms we report the gap to optimality expressed as a percentage. This value is computed as the ratio of the difference in weight between the maximum and approximate weights and the maximum weight. The Suitor algorithm computes more than 94% of the maximum weight in geometric mean, with the lowest value being 78%. ROMA obtains a lower gap than RAMA, and the ROMA algorithm initialized with the Suitor algorithm computes the highest weight among the $(2/3-\unicode[STIX]{x1D716})$ -approximation algorithms, but the $(1-\unicode[STIX]{x1D716})$ -algorithm with $\unicode[STIX]{x1D716}=1/4$ lowers the gap further to about 1.6%. Note that all these algorithms obtain much better weights than their worst-case approximation ratios.

Finally we compare the cardinalities of the maximum weight matching and approximate matchings in Table 3.4. Note that the cardinality of the exact algorithm is not necessarily equal to that of the maximum cardinality matching. The gap in cardinality is computed analogously to the gap in weights. The cardinality of the approximation algorithms is lower than that of the exact algorithm. The gap in cardinalities is about 3% for Suitor, while it is less than 0.5% for the other approximation algorithms.

The relative performance of these algorithms has a strong dependence on the weights. When we use these algorithms to solve vertex-weighted matchings (with vertex weights chosen uniformly at random, and edge weights computed by adding the vertex weights on the endpoints of an edge), there are problems for which the exact algorithms do not terminate in hundreds of hours. For these problems the $(2/3-\unicode[STIX]{x1D716})$ -algorithms tend to obtain better weights than the $(1-\unicode[STIX]{x1D716})$ -approximation algorithm. We discuss these in Section 5 on vertex-weighted matchings.

3.3.2 Comprehensive comparison of approximation algorithms

Now we compare the various approximation algorithms using original edge weights from the graphs rather than the random integer weights used in the previous experiments. We first compare the $1/2$ -approximation algorithms. Figure 3.3 shows the run time of the Global Paths Algorithm (GPA), and the Greedy, Suitor, Path Growing-DP and M1M2 algorithms for the ten graphs. For Greedy and Path Growing-DP the time to sort the edges is shown separately. This is done using the standard $qsort$ routine in C. Thus for an application where the edges are not already sorted, then this time must be added to get the total run time.

For these experiments we used a Dell computer running GNU/Linux equipped with four 10 core 2.00 GHz Intel Xeon E7-4850 processors with a total of 128 GB memory. Algorithms were implemented in C using OpenMP for parallelization and compiled with gcc using the –O3 flag.

Figure 3.3.

Run times for five $1/2$ -approximation algorithms.

As can be seen from the figure, for all but one graph the time to sort dominates the time taken by all other algorithms, in many cases by as much as an order of magnitude. When one excludes the time to sort, Greedy is the fastest algorithm followed by Suitor and Path Growing-DP. These again outperform M1M2 and the Global Paths Algorithm. If one includes the time to sort the edges then Suitor outperforms Greedy while M1M2 outperforms the Global Paths Algorithm.

To measure the quality of the solutions given by these algorithms we compute the average performance ratio for each algorithm compared to the best algorithm for each graph. This shows that Greedy, Path Growing-DP, M1M2 and Path Growing algorithms have an average performance ratio of 2.44%, 1.67%, 0.21% and 0.04%, respectively, higher than the best algorithm for each graph. Thus it follows that the more time-consuming algorithms also gives better results.

Next, we compare the use of either RAMA or ROMA for post-processing a solution. We first note that in terms of run time RAMA is about 10% slower than ROMA. Each application of ROMA is on average 4.9 times slower than the corresponding edge sorting. Thus there is a fairly large cost for using this post-processing. As there is also a slight advantage of using ROMA over RAMA in terms of quality of solution, we only present results from using ROMA.

How much improvement ROMA gives depends on the starting configuration. The average improvement for using Suitor followed by ROMA was 4.9% while both M1M2 and Global Paths Algorithm were improved by 2.6%. In terms of absolute quality of the final solution Suitor followed by ROMA always gave the best solution. On average M1M2 was 0.19% worse while Global Paths Algorithm was 2.11%. Thus, in terms of both run time and quality, the preferred strategy is to use Suitor followed by ROMA.

In Figure 3.4 we show the speedup curves for Suitor, M1M2 and Local Max algorithms when run on the graph nlpkkt200. For both Suitor and M1M2 the speedups increase almost linearly as long as there are available threads. Although the CPUs can perform hyper-threading, this does not give any additional speedup when using 50 threads. The speedup for Local Max drops off earlier than Suitor. The sequential run time for Local Max is about 70% higher than that of Suitor. Combined with the lower speedup it follows that the parallel run time of Local Max is even worse.

Figure 3.4.

Speedup for Suitor, M1M2 and Local Max algorithms on the nlpkkt200 graph.

Approximation algorithms for edge-weighted matchings have been implemented on GPUs by several authors. We list some of these papers: Cohen and Castonguay (Reference Cohen and Castonguay2012), Fagginger Auer and Bisseling (Reference Fagginger Auer, Bisseling, Keller, Kramer and Weiss2012), Halappanavar et al. (Reference Halappanavar, Feo, Villa, Dobrian and Pothen2012) and Naim et al. (Reference Naim, Manne, Halappanavar, Tumeo and Langguth2015).

3.4.1 Maximum cardinality matching

Maximum cardinality matchings in bipartite graphs have been employed in sparse matrix algorithms for several purposes. A maximum matching can be used to permute a sparse matrix to place the largest number of non-zeros on the diagonal. This number is the structural rank of the matrix, which is an upper bound on the numerical rank.

In the following we say that a matrix has a row-perfect matching to mean that its bipartite graph consisting of row vertices and column vertices has a matching in which all rows are matched. The concept of a column perfect matching is similar, and a perfect matching is both row-perfect and column-perfect.

A maximum matching can be used to compute the block (lower) triangular form (BTF) of a sparse matrix, which has the block matrix structure

$$\begin{eqnarray}A=\left(\begin{array}{@{}ccc@{}}A_{11} & A_{12} & A_{13}\\ 0 & A_{22} & A_{23}\\ 0 & 0 & A_{33}\end{array}\right)\!,\end{eqnarray}$$

where $A_{11}$ has fewer rows than columns and has a row-perfect matching, $A_{22}$ is a square matrix with a perfect matching, and $A_{33}$ has fewer columns than rows and has a column-perfect matching. The submatrices $A_{12}$ , $A_{13}$ and $A_{23}$ could be zero. The submatrix $A_{11}$ has the strong Hall property with respect to its rows (Coleman, Edenbrandt and Gilbert Reference Coleman, Edenbrandt and Gilbert1986, Pothen and Fan Reference Pothen and Fan1990). This property states that every set of $k$ rows has non-zeros in at least $(k+1)$ columns. The submatrix $A_{33}$ has the strong Hall property with respect to its columns, and $A_{22}$ has the strong Hall property with respect to both rows and columns (in this last case $k<n$ , where $n$ is the dimension of  $A_{22}$ ).

The BTF of a sparse matrix is derived from the Dulmage–Mendelsohn decomposition of the bipartite graph of the matrix, which could be computed from a maximum matching in the graph. The subgraph corresponding to $A_{11}$ is obtained by following all alternating paths from unmatched columns, and including all rows and columns reached in the subgraph. The subgraph corresponding to $A_{33}$ is obtained by following all alternating paths from unmatched rows, and including all rows and columns reached. The submatrix corresponding to $A_{22}$ consists of all rows and columns not reached thus far, and they must be matched to each other since all unmatched columns and rows have been accounted for. (There is a finer decomposition for $A_{22}$ which we do not describe here.) The decomposition is unique for a bipartite graph and is independent of the specific maximum matching used to induce it. Further details may be found in Pothen and Fan (Reference Pothen and Fan1990).

The BTF can be used to reduce the work in solving sparse linear systems of equations by a factorization-based algorithm, since only the diagonal blocks in the BTF need to be factored. The BTF has been used to solve Kirchoff’s equations in circuit design within the Xyce circuit simulation code (Keiter et al. Reference Keiter, Thornquist, Hoekstra, Russo, Schiek, Rankin, Li, Silveira and Feldmann2011), and it is also used to solve nonlinear systems of equations in systems modelling software such as Modelica (Fritzson Reference Fritzson2014). The strong Hall matrices in the BTF could be used to correctly predict the non-zero structures of orthogonal-triangular ( $QR$ ) factors of sparse matrices (Coleman et al. Reference Coleman, Edenbrandt and Gilbert1986, Pothen Reference Pothen1993). Matchings have also been used to compute sparse bases for the null space of sparse, under-determined matrices, as also sparse bases for the column space of such matrices (Coleman and Pothen Reference Coleman and Pothen1987, Pinar, Chow and Pothen Reference Pinar, Chow and Pothen2006).

In applications such as the BTF, one needs to compute a maximum cardinality matching in bipartite graphs, and an approximation will not suffice. Duff, Kaya and Uçar (Reference Duff, Kaya and Uçar2011) have addressed the several choices one needs to obtain an efficient practical algorithm. The Karp–Sipser algorithm is used as an initialization, and then there are choices to be made on if the augmenting path searches should use breadth-first search (BFS) or depth-first search (DFS) or some combination of them. It is also clear from their work and other authors, that the theoretically less efficient $O(nm)$ time algorithms are faster than the Hopcroft–Karp algorithm with $O(n^{1/2}m)$ time complexity. Azad, Buluç and Pothen (Reference Azad, Buluç and Pothen2017) have designed efficient shared-memory parallel algorithms for this problem. A distributed-memory parallel algorithm scaling to hundred-thousand cores has been implemented using sparse matrix vector products using the Graph-BLAS operations by Azad and Buluç (Reference Azad and Buluç2016).

3.4.2 Maximum edge-weighted matchings

An important application we consider is sparse Gaussian elimination ( $LU$ factorization). Here, permuting the sparse matrix to have large elements on the diagonal before the numerical factorization makes it less likely that numerical pivoting (row permutations) will be required during the computation to prevent large element growth in the factors. Olschowka and Neumaier (Reference Olschowka and Neumaier1996) described a pivoting strategy using a primal–dual algorithm for the assignment problem, which was then adapted and implemented for sparse bipartite graphs by Duff and Koster (Reference Duff and Koster2001). A perfect matching that has the maximum product of matched elements is computed in this context, and the resulting code, MC64, is widely used. Hogg and Scott (Reference Hogg and Scott2013) discuss several pivoting strategies, including matching-based orderings, to solve symmetric, indefinite systems of equations.

As parallel computers are able to solve systems of linear equations with millions of rows and columns, one challenge that has remained here is that the primal–dual algorithm used to compute perfect matchings does not have much concurrency. Hogg and Scott (Reference Hogg and Scott2015) have employed the auction algorithm on shared-memory machines to compute the matching in parallel, but scaling to the large numbers of cores on distributed-memory machines remained an open problem. This raised the possibility that approximation algorithms for matching could be employed in this context. In recent work Azad et al. (Reference Azad, Buluç, Li, Wang and Langguth2018) have described an approach that scales to 17 000 cores, and we consider this next.

Azad et al. (Reference Azad, Buluç, Li, Wang and Langguth2018) first compute a perfect matching in the bipartite graph and then seek to increase its weight. (A perfect matching must exist if the matrix is non-singular.) The perfect matching is computed by a distributed-memory parallel algorithm that recasts the problem in terms of matrix–vector products in a suitable semi-ring using the Graph-BLAS operations (Azad and Buluç Reference Azad and Buluç2016). The increase in weight is accomplished by means of the $(2/3-\unicode[STIX]{x1D716})$ -approximation algorithm of Pettie and Sanders (Reference Pettie and Sanders2004), which can be simplified here because the graph is bipartite and the matching is perfect. Hence they search for weight-increasing cycles with four edges (see Figure 3.1), finding for each vertex a cycle that leads to the largest increase in matching weight. From the set of cycles obtained, a maximal set of vertex-disjoint cycles is chosen using a greedy algorithm, and the matching weight is increased without affecting the cardinality of the perfect matching. The approximation guarantee of this algorithm is $(2/3-\unicode[STIX]{x1D716})$ and it is computed in $O(m\log \unicode[STIX]{x1D716}^{-1})$ time.

In a distributed-memory parallel implementation, it is not easy to find a maximal set of vertex-disjoint cycles without extensive communication, and hence these authors instead use only local comparisons to obtain a heavy set of vertex-disjoint cycles. Thus the parallel algorithm does not have the approximation guarantee, but it computes a heavy-weight matching that performs well in practice; it runs in $O((m/p)(1+\unicode[STIX]{x1D6FD})+\unicode[STIX]{x1D6FC}p)$ time, where $p$ is the number of processors, $\unicode[STIX]{x1D6FC}$ is the latency and $\unicode[STIX]{x1D6FD}$ is the inverse bandwidth. The authors report results for the algorithm on distributed-memory machines with 17 000 cores. Their code has been incorporated into the distributed-memory SuperLU solver for sparse, unsymmetric, linear systems of equations (Li and Demmel Reference Li and Demmel2003). Additional details on the applications of matchings in sparse matrix algorithms is provided in Duff and Uçar (Reference Duff, Uçar, Naumann and Schenk2012).

Exact algorithms for maximum cardinality matchings and maximum edge-weighted matchings in bipartite and non-bipartite graphs are available in the LEDA library (Mehlhorn and Näher Reference Mehlhorn and Näher1999). We compare the performance of approximation algorithms for vertex- and edge-weighted matchings with the LEDA codes in Section 5.

4.2.1 The Greedy algorithm

Algorithm 8 describes the Greedy algorithm for computing a $b$ -matching in a graph  $G$ .

The Greedy algorithm is a $1/2$ -approximation algorithm for the maximum weight $b$ -matching problem. We prove this by an argument that employs matroid theory due to Mestre (Reference Mestre2006).

Let $E$ be a set of elements, and ${\mathcal{I}}$ a collection of subsets of  $E$ . The tuple $(E,{\mathcal{I}})$ is a matroid if

1. (i) for all $A\subseteq B$ , if $B\in {\mathcal{I}}$ , then $A\in {\mathcal{I}}$ , and

2. (ii) for all $A,B\in {\mathcal{I}}$ with $|A|<|B|$ , there exists an element $x\in B\setminus A$ such that $A+x\in {\mathcal{I}}$ .

The sets in ${\mathcal{I}}$ are called the independent sets of the matroid, and by the first property, the empty set is independent. Maximal independent sets are called the bases of the matroid, and by the second property, all bases have the same cardinality. The reader unfamiliar with matroids may find two examples helpful. If $E$ is the set of columns of a matrix, and independence corresponds to linear independence in a vector space, then we have a matric matroid. If $E$ is the set of edges of a connected undirected graph without loops, and a subset of edges is defined to be independent if the edges do not induce a cycle, then we have a graphic matroid. A basis here corresponds to a spanning tree of the graph.

Assign every element in $E$ a non-negative weight, define the weight of a set as the sum of the weights of its elements, and consider the problem of finding an independent set of largest weight. The Greedy algorithm begins with the empty set, and at each step adds an element of largest weight if its addition will preserve independence, and rejects it otherwise. The Greedy algorithm finds an independent set (a basis) of maximum weight if and only if the tuple $(E,{\mathcal{I}})$ is a matroid.

It is well known that a matching does not correspond to a matroid, but a matching in a bipartite graph is the intersection of two matroids. We now consider a more general system called a $k$ -extendible system related to a matroid. A $k$ -extendible system is a tuple $(E,{\mathcal{I}})$ that satisfies property (i) of a matroid, but property (ii) is replaced by

1. (ii′) for all $A\subseteq B\in {\mathcal{I}}$ and $x\in E$ , if $A+x\in {\mathcal{I}}$ but $B+x\notin {\mathcal{I}}$ , then there exists $Y\subseteq B\setminus A$ , such that $B\setminus Y+x\in {\mathcal{I}}$ and $|Y|\leq k$ .

In words, this means that if we can augment a smaller independent set $A$ with a new element $x$ and preserve its independence, then it should be possible to remove a subset of size at most $k$ from any superset of $A$ that is independent, add $x$ to the new set and still preserve independence. Maximal independent sets do not necessarily have unique maximum cardinality for $k$ -extendible systems.

As a concrete example, let us consider an undirected graph $G=(V,E,w)$ with vertex set  $V$ , edge set $E$ and non-negative weights on its edges $W$ . Further, we are given a set of natural numbers $b(v)\leq \deg (v)$ for each $v\in V$ . Let the collection of independent sets ${\mathcal{I}}$ be defined as subsets of edges that consist of a $b$ -matching in  $G$ , i.e. subsets of edges $M$ such that $\deg _{M}(v)\leq b(v)$ for all $v\in V$ . Let $A$ be a $b$ -matching in $G$ and let $B\supset A$ be another $b$ -matching. Let $x=(u,v)$ be an edge that could be added to $A$ but not to $B$ to obtain a larger $b$ -matching. The reason that this edge cannot be added to $B$ is that the values of one or both of $b(u)$ and $b(v)$ would be exceeded. By removing one edge incident on $u$ and another edge incident on $v$ from $B$ , we would be able to add the edge $(u,v)$ to $B$ and preserve a $b$ -matching. Thus we have shown that a $b$ -matching is a 2-extendible system.

We now provide some concepts and prove a preliminary lemma in order to prove that the Greedy algorithm computes a $1/2$ -approximation algorithm. We prove a more general result for $k$ -extendible systems. An extension of an independent set $A\in {\mathcal{I}}$ is a superset $B$ of $A$ with $B\in {\mathcal{I}}$ . We denote an extension of maximum weight of a set $A$ by $\text{OPT}(A)$ .

Let the Greedy algorithm choose elements $x_{1},x_{2},\ldots ,x_{l}$ on a maximization problem on a $k$ -independence system $(E,{\mathcal{I}})$ . We denote the set of chosen elements at the end of the $i$ th iteration of the Greedy algorithm by  $S_{i}$ , with $S_{0}=\emptyset$ and $S_{i}=S_{i-1}+x_{i}$ . Note that $\text{OPT}(\emptyset )$ is the optimal solution to the maximization problem, and $\text{OPT}(S_{l})=S_{l}$ , since the set $S_{l}$ is maximal. Furthermore, we have

$$\begin{eqnarray}w(\text{OPT}(S_{0}))\geq w(\text{OPT}(S_{1}))\geq \cdots \geq w(\text{OPT}(S_{l})),\end{eqnarray}$$

since with increasing index $i$ , $S_{i}$ has fewer extensions available than $S_{i-1}$ .

Lemma 4.1. If $(E,{\mathcal{I}})$ is a $k$ -extendible system, then the element $x_{i}$ chosen by the Greedy algorithm at the $i$ th step satisfies

$$\begin{eqnarray}w(\text{OPT}(S_{i-1}))\leq w(\text{OPT}(S_{i}))+(k-1)\,w(x_{i}).\end{eqnarray}$$

Proof. By the definition of $k$ -extendibility, $S_{i-1}\in {\mathcal{I}}$ and $\text{OPT}(S_{i-1})\in {\mathcal{I}}$ . Now since the Greedy algorithm chooses the element $x_{i}$ at the $i$ th step, it does not belong to $S_{i-1}$ , but it could belong to $\text{OPT}(S_{i-1})$ . If it does, then $\text{OPT}(S_{i-1})=\text{OPT}(S_{i})$ , and the lemma holds trivially. Hence consider the situation when the element $x_{i}\not \in \text{OPT}(S_{i-1})$ . By $k$ -extendibility, there exists $Y\in \text{OPT}(S_{i-1})\setminus S_{i-1}$ such that $\text{OPT}(S_{i-1})\setminus Y+x_{i}\in {\mathcal{I}}$ , where $|Y|\leq k$ . We have the relationships

$$\begin{eqnarray}\displaystyle w(\text{OPT}(S_{i-1}) & = & \displaystyle w(\text{OPT}(S_{i-1})\setminus Y+x_{i})+w(Y)-w(x_{i}),\nonumber\\ \displaystyle & \leq & \displaystyle w(\text{OPT}(S_{i}))+w(Y)-w(x_{i}).\nonumber\end{eqnarray}$$

The first line of the equation above follows since $Y$ belongs to $\text{OPT}(S_{i-1})$ but $x_{i}$ does not. The second follows because $\text{OPT}(S_{i-1})\setminus Y+x_{i}$ is an extension of the set $S_{i}=S_{i-1}+x_{i}$ by the choice of $Y$ , and $\text{OPT}(S_{i})$ is an extension of maximum weight of the set  $S_{i}$ .

Now we show that any element $y\in Y$ is not heavier than the element $x_{i}$ , i.e. $w(y)\leq w(x_{i})$ . Assume for the sake of contradiction that $w(y)>w(x_{i})$ . If $y$ is heavier than  $x_{i}$ , then since $y\notin S_{i-1}$ by the choice of  $Y$ , the element $y$ must have been considered by the Greedy algorithm before $x_{i}$ was included in  $S_{i}$ , but was not included in the greedy solution. Thus there exists $j\leq i$ such that $S_{j}+y\notin {\mathcal{I}}$ , but $S_{j}+y\subseteq \text{OPT}(S_{i-1})\in {\mathcal{I}}$ . But this contradicts property (i) of a $k$ -independence system. Thus $w(y)\leq w(x_{i})$ , and since all weights are positive, $w(Y)\leq kw(x_{i})$ .◻

Proof. Let $\{x_{i}:i=1,2,\ldots ,l\}$ be the elements chosen by the Greedy algorithm, and the corresponding sets of chosen elements be denoted by $\{S_{i}:i=1,2,\ldots ,l\}$ . We apply Lemma 4.1  $l$ times, beginning with the empty set $S_{0}$ , to get

$$\begin{eqnarray}\displaystyle w(\text{OPT}(S_{0})) & \leq & \displaystyle w(\text{OPT}(S_{1}))+(k-1)\,w(x_{1})\nonumber\\ \displaystyle & \leq & \displaystyle w(\text{OPT}(S_{l}))+(k-1)\,\mathop{\sum }_{i=1}^{l}w(x_{i})\nonumber\\ \displaystyle & = & \displaystyle w(S_{l})+(k-1)\,w(S_{l})=k\,w(S_{l}).\nonumber\end{eqnarray}$$

The first term in last line of the equation follows from the fact that the set $S_{l}$ is maximal, and the second term is obtained by the summing the weight of the elements in  $S_{l}$ . Since $\text{OPT}(S_{0})$ is an independent set of maximum weight, we obtain the $1/k$ -approximation property of the Greedy algorithm.◻

4.2.2 The $b$ -Suitor algorithm

The $b$ -Suitor algorithm, designed by Khan et al. (Reference Khan, Pothen, Patwary, Satish, Sundaram, Manne, Halappanavar and Dubey2016b ), computes a $1/2$ -approximate $b$ -matching; indeed it computes a matching identical to the one obtained by the Greedy algorithm, provided ties in weights are broken consistently in both algorithms. This algorithm is based on proposals, much like the algorithms for the stable marriage problem and its variants (here the stable fixtures problem), and is a generalization of the Suitor algorithm. Vertices can propose to their heaviest neighbours, and these proposals may be reciprocated or annulled by other vertices. Two vertices are matched when both propose to each other. Unlike the Greedy algorithm, the $b$ -Suitor algorithm does not need to process edges in non-increasing order of weights; instead, it can process vertices in any order, although each vertex has to make proposals to its available neighbours in non-increasing order of weights.

The pseudo-code for the $b$ -Suitor algorithm is described in Algorithm 9. It maintains a priority queue $S$ to track the proposals made in the algorithm. The queue $S(u)$ consists of the suitors of a vertex $u$ , i.e. those neighbours of $u$ that currently have an active proposal to $u$ . The operation $S(u).insert(v)$ adds a vertex $v$ to the priority queue $S(u)$ , and $S(u).remove(v)$ removes it. The value $r(u)$ is the number of proposals that a vertex $u$ currently has made to its neighbours. We keep track of the lowest weight of a proposal received by a vertex $u$ (made by a suitor of $u$ ) in $S(u).\text{last}$ . If $u$ has received fewer than $b(u)$ proposals, this value is $\operatorname{NULL}$ .

In each iteration of the outer while loop, the algorithm processes vertices that have their current $b(.)$ values unsatisfied, and makes the requisite number of proposals to satisfy the $b(.)$ values; these vertices are stored in a set  $Q$ . During the iteration it collects vertices whose $b(.)$ values may be decremented in a set $Q^{\prime }$ so that they could be processed in the next iteration. For each vertex $u\in Q$ , while its $b(u)$ value is not satisfied and $\text{adj}(u)$ has not been exhaustively searched, the algorithm finds eligible neighbours to propose to. A neighbour $v$ is an eligible neighbour of $u$ if it holds fewer than $b(v)$ proposals, or $u$ can beat the lowest offer that $v$ currently has, which is $S(v).\text{last}$ . If $u$ cannot beat this offer, then it considers its next heaviest neighbour, and this prevents proposals being extended which at the current stage of the algorithm have no chance of success, saving work. But if $v$ has fewer than $b(v)$ proposals, or if $u$ can beat the lowest offer that $v$ has, then $u$ proposes to  $v$ , and becomes a suitor of  $v$ . In the latter case, $u$ annuls the lowest-weight proposal of  $v$ , say from a vertex  $y$ , and $y$ has to make another proposal in the next iteration. The value of $S(v).\text{last}$ is also updated.

Algorithm 9 shows that when a proposal made by a vertex $y$ is annulled, it is placed into a queue to be processed in the next iteration. This corresponds to the Gale–Shapley algorithm for stable marriage which processes proposals in rounds. Instead, one could consider the vertex $y$ making a proposal immediately, and this would correspond to the McVitie–Wilson algorithm. It is the possibility of proposals being annulled that permits vertices to be processed in any order, thus increasing the concurrency in the algorithm. We have shown that if the edge weights are chosen uniformly at random, then the total expected work in the algorithm is linear in the number of edges in the graph, so that the algorithm does not create a lot of unnecessary work (Khan et al. Reference Khan, Pothen and Ferdous2018b ). This paper also shows that the depth (the length of the critical path) is $O(\log m\log \unicode[STIX]{x1D6E5})$ .