# Engineering fast algorithms for the bottleneck matching problem

Ioannis Panagiotas, Grégoire Pichon, Somesh Singh, Bora Uçar

► To cite this version:

Ioannis Panagiotas, Grégoire Pichon, Somesh Singh, Bora Uçar. Engineering fast algorithms for the bottleneck matching problem. ESA 2023 - The 31st Annual European Symposium on Algorithms, Sep 2023, Amsterdam (Hollande), Netherlands. hal-04146298v2

HAL Id: hal-04146298

https://inria.hal.science/hal-04146298v2

Submitted on 14 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Engineering fast algorithms for the bottleneck matching problem

Ioannis Panagiotas

Neo4j, Sweden

Grégoire Pichon

Université Claude Bernard Lyon 1 and LIP, France

UMR5668 (CNRS, ENS de Lyon, Inria, UCBL1) France

Somesh Singh

CNRS and LIP, ENS de Lyon, France

UMR5668 (CNRS, ENS de Lyon, Inria, UCBL1) France

Bora Uçar

CNRS and LIP, ENS de Lyon, France

UMR5668 (CNRS, ENS de Lyon, Inria, UCBL1) France

## Abstract

We investigate the maximum bottleneck matching problem in bipartite graphs. Given a bipartite graph with nonnegative edge weights, the problem is to determine a maximum cardinality matching in which the minimum weight of an edge is the maximum. To the best of our knowledge, there are two widely used solvers for this problem based on two different approaches. There exists a third known approach in the literature, which seems inferior to those two which is presumably why there is no implementation of it. We take this third approach, make theoretical observations to improve its behavior, and implement the improved method. Experiments with the existing two solvers show that their run time can be too high to be useful in many interesting cases. Furthermore, their performance is not predictable, and slight perturbations of the input graph lead to considerable changes in the run time. On the other hand, the proposed solver's performance is much more stable; it is almost always faster than or comparable to the two existing solvers, and its run time always remains low.

**2012 ACM Subject Classification** Theory of computation

**Keywords and phrases** bipartite graphs, assignment problem, matching

**Related Version** Published version: https://doi.org/10.4230/LIPIcs.ESA.2023.87

**Supplementary Material** Source code

**Software:** https://gitlab.inria.fr/bora-ucar/bottled

## 1 Introduction

A matching in a graph is a set of edges without any common vertices. A maximum cardinality matching in a graph has the largest number of edges among all matchings. We investigate algorithms for finding a maximum cardinality matching whose minimum edge weight is the maximum on bipartite graphs with edge weights. This is called the bottleneck matching problem or linear bottleneck assignment problem when all vertices can be matched [3, Section 6.2]. Formally, the bottleneck matching problem is to find a maximum cardinality matching $M$ which maximizes

$$ \min_{(r_i, c_j) \in M} w_{i,j}. $$

This problem can be solved in polynomial time.
## 2 Algorithms for bottleneck matching

The bottleneck matching problem arises in different contexts [3, Section 6.2.8]. We are motivated by the Birkhoff-von Neumann (BvN) decomposition of doubly stochastic matrices [2], which arises in practical applications [1, 4, 16, 21]. In this case, the bipartite graphs associated with matrices have equal number of vertices on both sides and contain perfect matchings. A known heuristic for the BvN decomposition [10] repeatedly calls a bottleneck matching algorithm on the bipartite graph of a dynamically changing matrix.

<table><caption>Table 1 Run time, in seconds, of MC64J2 and MC64J3 on 12 problems. **A** is the original matrix, **AP** is the same matrix with a random column permutation; **S** = **DP(A)E**, where **P(A)** is the 0-1 matrix with 1s at the nonzero positions of **A**; and **D** and **E** are positive diagonal matrices scaling **P(A)** to be doubly stochastic.</caption><thead><tr><th rowspan="2">matrix</th><th colspan="2">A</th><th colspan="2">AP</th><th colspan="2">S</th></tr><tr><th>J2</th><th>J3</th><th>J2</th><th>J3</th><th>J2</th><th>J3</th></tr></thead><tbody><tr><td>atmosmodm</td><td>0.05</td><td>0.26</td><td>0.13</td><td>0.40</td><td>684.87</td><td>6487.07</td></tr><tr><td>CurlCurl_3</td><td>0.05</td><td>0.31</td><td>0.15</td><td>0.48</td><td>578.83</td><td>2958.07</td></tr><tr><td>ss</td><td>2.28</td><td>0.90</td><td>2317.46</td><td>1.73</td><td>810.10</td><td>16501.00</td></tr><tr><td>vas_stokes_2M</td><td>0.25</td><td>1.87</td><td>6310.44</td><td>727.57</td><td>1420.01</td><td>4571.14</td></tr></tbody></table>

The software MC64 [7, 8] is the state-of-the-art and implements two algorithms, denoted MC64J2 and MC64J3, for the bottleneck matching problem. To the best of our knowledge these are the only available codes that can handle bipartite graphs corresponding to large sparse matrices. Their worst-case run time are $O(n(m + n) \log_2 n)$ and $O(nm \log_2 n)$, on bipartite graphs with $n$ vertices on each side and $m$ edges [8]. MC64, especially the newer MC64J2 [8], is well engineered. It works very well for graphs corresponding to matrices from numerical applications. However, it does not have stable run time behavior in two senses. First, when run on the same bipartite graph twice with different edge weights the difference in run time can be in the order of hours. Second, on two equivalent problem instances, where one is obtained from the other by just reordering the vertices, the run time can change dramatically. We report the run time of MC64 on four matrices, from the SuiteSparse Matrix Collection [5], in Table 1 to explain this—more experiments of similar nature are in Section 4. Here, the set of rows and the set of columns of a matrix correspond to the two parts of the bipartite graph with an edge between two vertices if the corresponding entry in the matrix is nonzero, and the nonzero values are the edge weights. The table contains results for **A**, for **AP** where **P** is a random permutation, and for the doubly stochastic matrix **S** which is obtained by scaling the pattern of **A** with Sinkhorn-Knopp algorithm [20].

The bipartite graphs of **A** and **AP** are the same apart from renumbering of the vertices in one part. While on **A** both MC64J2 and MC64J3 are fast, both methods suffer on **AP**; the run time of MC64J2 is not acceptable for the last two instances, and that of MC64J3 is high for the last one in Table 1. The bottleneck matching problems on **A** and **AP** are essentially the same, as permuting the columns does not change the values, nor the bottleneck matching and its value. The bipartite graph of $S = DP(A)E$ is the same as that of **A** with different edge weights, hence the problems are not equivalent. Now, the run time of both methods is too much for all instances. The wildly varying run time of both MC64 routines on **S** in comparison to those on **A** further highlight the instability in their performance.

Our aim in this paper is to develop an algorithm for the bottleneck matching problem which is better than the state-of-the-art codes in MC64. For this purpose, we study an overlooked alternative from the literature. We make observations that pave the way for an efficient algorithm, implement and compare it against the codes from MC64. We conduct a large set of experiments to show that our approach is usually much faster than MC64 and in
addition exhibits stable and robust performance.

Section 2 gives a brief background and a summary of the known algorithms. Section 3 contains the proposed algorithm. Section 4 presents the experimental results, and Section 5 concludes the paper. The appendix contains detailed description of the experiments and further experiments.

## 2 Background and related work

A matrix $\mathbf{A}$ with $nnz(\mathbf{A})$ nonzero entries can be represented with a bipartite graph $G = (R \cup C, E)$ where each row of $\mathbf{A}$ corresponds to a unique vertex in $R$, each column of $\mathbf{A}$ corresponds to a unique vertex in $C$, and there is an edge $(r_i, c_j)$ whenever $a_{ij} \neq 0$. The number $m$ of edges in $G$ is thus equivalent to $nnz(\mathbf{A})$. When the edges are weighted, the weight of the edge $(r_i, c_j)$ is $|a_{ij}|$, that is the magnitude of the nonzeros of $\mathbf{A}$. For a vertex $v$, we use $\mathrm{adj}(v)$ to denote the set of its neighbors.

A matching is a set of edges with no common vertices. A matching is of maximum cardinality if it has the largest number of edges. Given a matching $\mathcal{M}$, a vertex is matched if an edge from $\mathcal{M}$ is incident on it and free otherwise. A matching is perfect if it matches all vertices. The deficiency of a matching $\mathcal{M}$ is the difference between the maximum cardinality of a matching and $|\mathcal{M}|$. Given a matching $\mathcal{M}$ in the graph $G$, a path in $G$ is $\mathcal{M}$-alternating if its edges are alternately in $\mathcal{M}$. An $\mathcal{M}$-alternating path $\mathcal{P}$ is $\mathcal{M}$-augmenting if the start and end vertices of $\mathcal{P}$ are both free. A vertex cover is a set of vertices that includes at least one vertex from each edge. In a bipartite graph the maximum cardinality of a matching is equal to the minimum cardinality of a vertex cover [3, Th. 2.7].

Given a bipartite graph, any of its maximum cardinality matchings can be used to obtain a canonical decomposition called Dulmage-Mendelsohn (DM) decomposition [11]. Based on the DM decomposition, Pothen and Fan [18] describe algorithms to permute sparse matrices in a block upper triangular form (BTF):

$$ \mathbf{A} = \begin{array}{c} H_R \begin{pmatrix} \mathbf{A}_H & * & * \\ O & \mathbf{A}_S & * \\ O & O & \mathbf{A}_V \end{pmatrix} \\ S_R \\ V_R \end{array}. \qquad (1) $$

In a BTF, the submatrix $\mathbf{A}_H$ has more columns than rows, and all rows in $H_R$ are matched to a column in $H_C$ in any maximum cardinality matching; the submatrix $\mathbf{A}_S$ is square with at least one perfect matching; the submatrix $\mathbf{A}_V$ has more rows than columns, and all columns in $V_C$ are matched to a row in $V_R$. The rows/columns in each block are defined as follows

$H_R = \{\text{row vertices reachable from free column vertices via alternating paths}\}$,

$H_C = \{\text{free column vertices or column vertices reachable from free column vertices via alternating paths}\}$,

$V_R = \{\text{free row vertices or row vertices reachable from free row vertices via alternating paths}\}$,

$V_C = \{\text{column vertices reachable from free row vertices via alternating paths}\}$,

$S_R = R \setminus (H_R \cup V_R)$, and

$S_C = C \setminus (H_C \cup V_C)$.

A standard BFS/DFS-based graph traversal algorithm will find these sets in linear time.
4 Algorithms for bottleneck matching

We make some observation on the BTF form of a matrix. First, the DM decomposition
reveals a minimum cover [3, Alg. 3.1]. As all nonzeros in the BTF (1) are confined in the
rows $H_R \cup S_R$ and the columns in $V_C$, the vertex set $C = H_R \cup S_R \cup V_C$ is a cover. Since
the cardinality of $C$ is equal to the maximum cardinality of a matching, it is a minimum
cover. Second, if one adds new nonzeros to the diagonal blocks, upper diagonal blocks, or to
the blocks $(S_R, H_C)$ and $(V_R, S_C)$, the maximum cardinality of a matching does not change.
This is so, as the new nonzeros cannot create augmenting paths.

Hall's theorem [13] states that for a bipartite graph $G$ to have a column-perfect matching, the relation $|S| \le |\bigcup_{c \in S} \text{adj}(c)|$ must hold for any subset $S$ of columns. If a similar relation holds for all subsets of rows, then $G$ will have perfect matchings.

An $n \times n$ matrix $\mathbf{A} \neq 0$ is called *doubly stochastic* if every entry is nonnegative and the sum of entries in each row and each column is equal to 1. Any nonnegative square matrix, whose bipartite graph has perfect matchings, can be scaled with two diagonal matrices to be doubly stochastic [20]. A *permutation matrix* is a square matrix where each row/column contains exactly one nonzero value equal to 1. A perfect matching in the bipartite graph representation of $\mathbf{A}$ corresponds to a permutation matrix. The bipartite graphs of doubly stochastic matrices have perfect matchings.

Henceforth, we assume that the given bipartite graph contains perfect matchings. We comment on rectangular matrices and matrices without perfect matchings in Section 3.3.

**2.1 Related work**

We review three algorithms from the literature [3, 7, 8]. The first two are implemented in MC64, and to the best of our knowledge are currently the best practical algorithms. Burkard et al. [3, Section 6.2.4] describe two other algorithms [12, 19], which are more theoretical.

**2.1.1 Shortest-augmenting path based algorithms**

Algorithms based on shortest augmenting paths start with a matching which has the maximum
bottleneck value for the currently matched vertices $C'$ in one part, say $C$. In order to augment
the matching, a shortest augmenting path from a free vertex $c$ of $C$ is found with a variant
of Dijkstra's shortest path algorithm. Augmenting along the shortest paths maintains the
invariant that the current matching has the maximum bottleneck value for any matching that
matches the vertices $C' \cup \{c\}$. The process continues until a perfect matching is obtained.

The state-of-the-art implementation in MC64 [8], MC64J2, starts by computing an upper bound $\omega$ on the bottleneck value, which is the minimum of maximum in each column and row. It then computes a maximal matching on the graph containing only edge weights no smaller than $\omega$, which is then improved by length-three augmenting paths in a preprocess step. Then, a shortest-augmenting path is sought from each free column vertex to solve the problem. MC64J2 implements an efficient adaptation of Dijkstra's algorithm to find these paths. Depending on the edge weights, the structure of the bipartite graph, or the visit order many edges and vertices may be visited while finding an augmenting path. As seen in Table 1, this can accumulate and result in very long run time.

**2.1.2 Threshold-based algorithms**

Let $G$ be a weighted bipartite graph. For a value $\omega$, let $G[\omega]$ contain only the edges of $G$
with weight at least $\omega$. Threshold-based algorithms find the largest $\omega$ for which $G[\omega]$ has a
perfect matching. They do so by considering different values for $\omega$, testing whether $G[\omega]$
has a perfect matching or not and tuning the next $\omega$ to a higher or lower value accordingly. MC64's implementation, MC64J3, discusses initialization and algorithmic choices to reduce the number of tests [7]. However, MC64 uses a depth-first-search algorithm to find the augmenting paths for each test which can slow it down, as seen in Table 1.

### 2.1.3 An algorithm based on duality

This algorithm is also threshold-based and uses the duality of matchings and coverings to find the next threshold. We explain this algorithm in more detail as ours improves upon it. Let $\omega$ be a value where $G[\omega]$ does not contain a perfect matching. Let $M$ be a maximum cardinality matching in $G[\omega]$. From $M$ one can define a minimum vertex cover $C = H_R \cup S_R \cup V_C$ of $G[\omega]$ with vertex sets $I = H_R \cup S_R$ and $J = V_C$ [3, Section 6.2.3]. Since $M$ is not perfect, there must be an edge in $G$ from a vertex in $\bar{I} = R \setminus I$ to another vertex in $\bar{J} = C \setminus J$. The value $\max_{i \in \bar{I}, j \in \bar{J}} a_{ij}$ thus cannot be smaller than the bottleneck matching value of $G$ and can be used as the next threshold. This approach is shown in Algorithm 1 [3, Section 6.2.3].

**Algorithm 1** Duality-based algorithm

**Input** : $G$, an edge weighted bipartite graph having perfect matchings
**Output**: $M$, a bottleneck perfect matching
Let $\omega$ be an upper bound on the bottleneck matching value
$M \leftarrow \emptyset$
**while** $|M| < n$ **do**
$M \leftarrow$ a maximum cardinality matching in $G[\omega]$ 
**if** $|M| < n$ **then**
Let $I \subseteq R$ and $J \subseteq C$ be the vertex sets of the associated cover of $G[\omega]$ 
$\omega \leftarrow \max_{i \in R \setminus I, j \in C \setminus J} w_{i,j}$
/* in $G$, not $G[\omega]$ */
**1**
**2**

The maximum cardinality matching in $G[\omega]$ can be found in $O(\sqrt{n} \text{nnz}(\mathbf{A}[\omega]))$ time in the worst case [14]. Once such a matching is found, the associated minimum cover and the maximum uncovered value at Line 2 can be obtained in linear time. Therefore, the worst case time complexity of Algorithm 1 is $O(\sqrt{n} \text{nnz}(\mathbf{A}))$ times the number of iterations. The worst case run time for Algorithm 1 can hence be too high. This is so as a new edge, due to the reduced $\omega$, does not mean one more edge in the maximum matching, and hence the while loop can even run for more than $n$ iterations.

## 3 The proposed algorithm

Our algorithm is based on Algorithm 1 and integrates threshold techniques. Let $G = (R \cup C, E)$ be an edge weighted bipartite graph, with $w_{i,j}$ being the weight of the edge $(r_i, c_j)$. Let $\omega$ be a nonnegative value, $G[\omega]$ be as before, and the bottleneck matching value $b^*$ be the largest $\omega$ for which $G[\omega]$ has a perfect matching. We call a value $\omega$ *safe*, when $\omega \ge b^*$. Algorithm 1 produces decreasing safe values that converge to $b^*$. We make a series of observations to find $b^*$ faster. The first observation is that there are several minimum vertex covers associated with a given maximum cardinality matching. Using the BTF (1), let $C' = H_R \cup S_C \cup V_C$. As all nonzeros are covered by $C'$ and $|C'| = |H_R \cup S_R \cup V_C|$, $C'$ is also a minimum cover. If we choose this cover, the sets $I$ and $J$ at Line 1 of Algorithm 1 become $H_R$ and $S_C \cup V_C$, in which case, the maximum value can be different. This leads to the following proposition.

**Proposition 1.** Let $C_1 = H_R \cup S_R \cup V_C$ and $C_2 = H_R \cup S_C \cup V_C$ be two minimum covers of $G[\omega]$ revealed by the BTF (1). Let $\omega_1$ and $\omega_2$ be the maximum values defined in Line 2 of Algorithm 1 for $C_1$ and $C_2$, respectively. Then, $\min(\omega_1, \omega_2)$ is safe.
6 Algorithms for bottleneck matching

**Proof.** Let $\omega$ be the current value, and $\omega_1 < \omega_2$ without loss of generality. This means that $\omega_2$ is in $\mathbf{A}(V_R, S_C)$ and all entries in $\mathbf{A}(V_R, H_C)$ are smaller than $\omega_2$. The maximum cardinality of a matching in $G[\omega_2]$ cannot be larger than that in $G[\omega]$, as the set $C_1 = H_R \cup S_R \cup V_C$ still covers all edges of $G[\omega_2]$, including those that arise in $\mathbf{A}(S_R, H_C)$. Hence the cover $C_1$ can be used to get the next $\omega$ in Line 2, which concludes the proof. ◀

Based on Proposition 1, one can use the smaller of the largest uncovered element in $(S_R \cup V_R, H_C)$ and that in $(V_R, H_C \cup S_C)$. We propose to exploit the two identified covers as much as possible for a faster convergence of $\omega$ to $b^*$. As we reason in Lemma 2 below, one can find a tighter bound on $b^*$, depending on the deficiency of the current matching as well as the (original) adjacencies of the vertices in the identified covers.

**Lemma 2.** Let $\mathcal{M}$ be a maximum cardinality matching in $G[\omega]$ with a deficiency of $k$ in $G$; $\mathbf{A}[\omega]$ and $\mathbf{A}$ be, respectively, the matrices associated with $G[\omega]$ and $G$; $C = H_R \cup S_C \cup V_C$ be a minimum vertex cover of $G[\omega]$ associated with $\mathcal{M}$ when $\mathbf{A}[\omega]$ is permuted in a BTF (1); Let $\omega_k^r$ be the $k$th largest element in $\bigcup_{i \in S_R \cup V_R} \max\{w_{i,j} : j \in H_C\}$, and $\omega_k^c$ the $k$th largest element in $\bigcup_{j \in H_C} \max\{w_{i,j} : i \in S_R \cup V_R\}$. Then, $\omega_k = \min(\omega_k^r, \omega_k^c)$ is safe.

**Proof.** Consider first the set $H_C$ of columns, and note that $|H_C| - |H_R| = k$ as all other columns are matched. By Hall's theorem, $|H_C| \le |\bigcup_{c \in H_C} \mathrm{adj}(c)|$ must hold in $\mathbf{A}$ as there is a perfect matching. Among all rows in $\bigcup_{c \in H_C} \mathrm{adj}(c)$, we have $|H_R|$ in the set $H_R$. Therefore there must be at least $k$ other nonzero rows in $\mathbf{A}(S_R \cup V_R, H_C)$. The element $\omega_k^r$ from $\bigcup_{i \in S_R \cup V_R} \max\{w_{i,j} : j \in H_C\}$ is safe as any value greater than that will cover less than $k$ rows and Hall's conditions cannot be satisfied. A similar argument applies to $\omega_k^c$ by considering the set $S_R \cup V_R$ of rows. In $\mathbf{A}[\omega]$ we have $\mathrm{adj}(S_R \cup V_R) = S_C \cup V_C$, and $|S_R \cup V_R| - |S_C \cup V_C| = k$. The element $\omega_k^c$ must be safe since we need at least $k$ nonzero columns in $\mathbf{A}(S_R \cup V_R, H_C)$. The value $\omega_k$ is thus safe as the minimum of two safe values. ◀

One can identify several minimum covers, collect the $k$th largest uncovered element with respect to each, and use the minimum of the collected elements as the next $\omega$. As finding these minimum covers can be expensive, we propose using the two which are readily revealed by the BTF. That is, we use $\omega = \min(\omega_1, \omega_2)$ where

$$ \omega_1^r, \omega_1^c = k\text{-th largest row and column maximum entries in } \mathbf{A}(S_R \cup V_R, H_C), \\ \omega_1 = \min(\omega_1^r, \omega_1^c), \qquad (2) $$

$$ \omega_2^r, \omega_2^c = k\text{-th largest row and column maximum entries in } \mathbf{A}(V_R, H_C \cup S_C), \\ \omega_2 = \min(\omega_2^r, \omega_2^c). \qquad (3) $$

This corresponds to applying Hall's theorem to the sets $H_C$ and $H_C \cup S_C$ of columns and to the sets $V_R$ and $S_R \cup V_R$ of rows.

## 3.1 Putting it all together

The proposed algorithm BOTTLED is shown in Algorithm 2. The input is a sparse matrix represented in the compressed storage by columns (CSC) format. BOTTLED creates a compressed storage by rows (CSR) representation of the input matrix. It then sorts the nonzeros in each row and each column in non-increasing order of their values. Then the threshold $\omega$ is initialized as the minimum of the $2n$ nonzero values consisting of the maximum in each row and maximum in each column. The algorithm then updates the threshold $\omega$ in a while-loop as in Algorithm 1. In the while-loop there are three subroutines:
MaximumCardinalityMatching, SAP, and DM-dec. These subroutines correspond to a maximum cardinality matching algorithm, a shortest augmenting path-based method to match a column, and an algorithm obtaining the row and column blocks of the BTF (1).

<caption>**Algorithm 2 BOTTLED:** The proposed bottleneck matching algorithm</caption>

**Input** : **A**, stored by columns
**Output**: $\mathcal{M}$, a bottleneck perfect matching
Create a CSR representation of **A**
Sort the nonzeros in each row and in each column in non-increasing order of values
1. $\omega \leftarrow \min$ of the maximum in each row, maximum in each column
   $G[\omega] \leftarrow (R \cup C, \emptyset)$ and $\mathcal{M} \leftarrow \emptyset$
**while** $|\mathcal{M}| < n$ **do**
2. for each $i$, release new $a_{ij} \ge \omega$ and for each $j$ release new $a_{ij} \ge \omega$ into $G[\omega]$
   if $|\mathcal{M}| = n-1$ **then**
      $\mathcal{M} \leftarrow \text{SAP}(G, \mathcal{M}, c)$ with the last free vertex $c$
   else
      $\mathcal{M}' \leftarrow \text{MaximumCardinalityMatching}(G[\omega], \mathcal{M})$
      if $|\mathcal{M}'| = n$ **then** $\mathcal{M} \leftarrow \mathcal{M}'$; break
      if $|\mathcal{M}'| = |\mathcal{M}|$ **then**
         select a vertex $c$
         $\mathcal{M} \leftarrow \text{SAP}(G, \mathcal{M}', c)$
      else
         $\mathcal{M} \leftarrow \mathcal{M}'$
         $\langle H_R, S_R, V_R, H_C, S_C, V_C \rangle \leftarrow \text{DM-dec}(G[\omega], \mathcal{M})$
         $\omega \leftarrow \min(\omega_1, \omega_2)$ with $\omega_1$ as in (2) and $\omega_2$ as in (3)

Algorithm 2 stores the edges of $G[\omega]$ over the storage of **A** in the CSC and CSR formats, without explicitly building adjacency lists. The start address of each row and column are the same as those of **A**. For each row/column of $G[\omega]$, we keep an end-pointer which points to the smallest nonzero of **A** in that row/column that is no smaller than $\omega$. These end-pointers are initialized before the while-loop in $O(n)$ time, and incremented at Line 2 at each iteration of the while-loop. Therefore the total run time cost of building $G[\omega]$ is $O(nnz(\mathbf{A}))$.

In Algorithm 2, $\text{SAP}(G, \mathcal{M}, c)$ refers to the algorithm summarized in Section 2.1.1. As stated before, SAP needs the current matching to have the bottleneck value among all matchings covering the same set of column vertices. The approach outlined in Algorithm 1 produces such matchings, that is why, at any point, one can resort to SAP. We invoke SAP in two cases: (i) when the deficiency is one; (ii) when an update of $\omega$ did not yield an increase in the cardinality of the current matching. The first case is straightforward. For the second case, we apply a simple heuristic to help the algorithm converge faster. As any free column vertex can be the start of an augmenting path, we choose $c \in C$ whose largest edge weight not included in the current $G[\omega]$ is minimum. Matching $c$ will lead to a reduced $\omega$, and the reduction will hopefully be large with this choice of $c$ (some empirical results are in the appendix of the related version). The most common algorithms for the maximum cardinality matching problem take an initial matching as input and augment it. This is very suitable at Line 3 of Algorithm 2, as we have a maximum cardinality matching on a graph, we add new edges, and then ask for a maximum cardinality matching in the new graph. We have used the code-base of MatchMaker [6, 15] to implement this step in the implicit representation of $G[\omega]$.

At Line 5 of Algorithm 2, we use a binary heap with a limit $k$ on its size. For $\omega_1^r$, the nonzeros of each row in $S_R \cup V_R$ with value smaller than $\omega$ are visited, and the largest element is used as a key in the heap. When the heap is full, keys are added only if they are larger than the current minimum, which is then removed. The minimum of the heap is returned as
8 Algorithms for bottleneck matching

$\omega_1^r$. Other quantities of (2) and (3) are computed similarly. All nonzeros of $\mathbf{A} - \mathbf{A}[\omega]$ in the rows $S_R \cup V_R$ and columns $H_C$ can be visited and the heap operations can be performed in $O(\text{nnz}(\mathbf{A}(S_R \cup V_R, H_C)) + (|H_C| + |S_R \cup V_R|) \log k)$ time. A similar analysis holds for $\mathbf{A}(S_R \cup V_R, H_C)$. The run time of one iteration of the while loop is thus dominated by that of the maximum cardinality matching algorithm. We do not have an estimate on the number of iterations of the while-loop; in the empirical results in Section 4.3 the number is less than what a binary search approach would yield.

3.2 In the context of a BvN decomposition method

The Birkhoff–von Neumann (BvN) theorem [2] states that a doubly stochastic matrix **A** can be written as **A** = $\sum_{i=1}^{l} \alpha_i \mathbf{P}_i$ where each **P**_i is a permutation matrix, and $\alpha_i$s are positive coefficients summing up to one. Such a decomposition is not unique, and the problem of finding a decomposition with the smallest $\ell$ value is NP-Complete [10].

One heuristic [10] for obtaining a BvN decomposition of a doubly stochastic matrix **A** works as follows. It finds the value **b** of a bottleneck matching whose pattern is the permutation matrix **P**, replaces **A** with **A** – **b****P**, and continues until a zero matrix is obtained. A property of this heuristic is that the successive bottleneck values are in a non-increasing order [9]. Our bottleneck matching algorithm is very fitting in this case. One can create the CSR representation and sort each row and column once. Then, executing the while loop of Algorithm 2 will obtain a bottleneck matching for the current matrix. Once **b** and **P** are obtained, replacing **A** with **A** – **b****P** can be done by subtracting **b** from each matched entry and updating that entry's position in the sorted list of both rows and columns in overall $O(n + \text{nnz}(\mathbf{A}))$ time. While doing so, one can update the end-pointers used for $G[\omega]$, and avoid the preprocessing in Algorithm 2 at subsequent invocations.

3.3 Rectangular matrices or matrices without perfect matchings

The case in which there are perfect matchings in the given bipartite graph is common in
applications where the bipartite graphs correspond to sparse matrices. This is especially so
in the BvN decomposition, which was our motivation. Nonetheless MC64's J2 and J3 work
for cases in which one part of the bipartite graph has more vertices than the other, where the
smaller side can be perfectly matched. This corresponds to $n_R \times n$ matrices for $n_R > n$ that
have column-perfect matchings. Our algorithm can handle this case either by initializing $\omega$
at Line 1 using only the column values, or by using those and only the $n$th maximum of the
set of $n_R$ maximum entries, one from in each row.

Consider now the most general case corresponding to $n_R \times n$ matrices, with $n_R \ge n$ and without column perfect matchings. In this case, MC64J2 returns a maximum cardinality matching, without necessarily finding the correct bottleneck value. That is so because not all vertices from which the shortest-augmenting paths are sought can be matched in a bottleneck maximum cardinality matching. We do not know of a suitable fix for this. MC64J3 on the other hand works correctly. The presented algorithm BOTTLED needs four minor modifications to handle this general case: (i) the maximum size of a matching $n'$ should be computed at the beginning; (ii) the initialization should use the smallest of the $n'$th maximum of $n$ maximum entries, one from each column, and the $n'$th maximum of the $n_R$ maximum entries, one from each row; (iii) at Line 5, the deficiency is $k = n' - |M|$; and (iv) the shortest-augmenting path method should not be used.
# 4 Experiments

We observed in a preliminary set of experiments that MC64J2 is generally faster than MC64J3; the geometric mean of the ratios of the run time of the latter to that of the former was 4.6. As already highlighted in Table 1, the run time of MC64J3 can be too large to experiment (this happens more frequently than with MC64J2). That is why we have implemented a threshold-based approach, referred here as THRESH, using the same code base of BOTTLED, and use it in our experiments instead of MC64J3. THRESH uses the initialization procedure common to MC64J3 and BOTTLED to find a large initial matching and an upper bound $\omega$ on the bottleneck value. If this matching is not maximum, then the edge weights are sorted to find the bottleneck value by a binary search between the smallest edge weight and the initial value $\omega$. The binary search uses only the available edge weights and each time half of the edges are discarded. Our codes, written in C, are available at https://gitlab.inria.fr/bora-ucar/bottled; the codes used in the experiments are elsewhere [17].

The next subsection describes the data set, and Section 4.2 conducts a performance analysis of BOTTLED. Section 4.3 compares BOTTLED with MC64J2 and THRESH. Last, some experiments using BOTTLED within a BvN decomposition heuristic are discussed in Section 4.4. Appendix of the related version contains detailed information about experiments.

## 4.1 Data set and measurements

We have experimented with all square matrices with perfect matchings, at least 100,000 rows, less than 250,000,000 nonzeros, and with no explicit zeros from SuiteSparse [5]. There were 113 such matrices at the time of experimentation. From each matrix, we created six types of problem instances, which are denoted as **A**, **DAE**, **DP(A)E**, **AP**, **DAPE**, and **DP(A)PE**. The **A**-type instance corresponds to the bipartite graph of the original matrix with the magnitudes of the entries as edge weights. The **DAE**-type and **DP(A)E**-type instances correspond, respectively, to the scaled version of the matrices and their patterns with 20 iterations of the Sinkhorn-Knopp [20] algorithm, and the other three types of instances are obtained from the first three by random column permutations. We discarded the instances **A** and **AP** for 0-1 matrices; for the same set of matrices **DAE**-type and **DP(A)E**-type instances are identical, and hence we kept only one of them. We discarded the instances in which the initialization algorithm found the bottleneck value, in which case all three methods are equivalent. We report the experiments with 14 **A**-type, 18 **DAE**-type, 58 **DP(A)E**-type instances, and the same number of instances with column permutations. For each instance, each algorithm is run five times and the geometric mean of the run time is reported as that algorithm's performance on that instance; when column permutations are applied, these correspond to five different permutations. When comparing two algorithms' run time, we do not include cases where both algorithms run in less than one second (as both are very small and the difference between the algorithms is insignificant). Appendix of the related version contains all the results.

We carry out the experiments on a machine having Intel(R) Xeon(R) CPU E7-8890 v4 with a clock-speed of 2.20 GHz, and 1.5TB memory. The machine runs Debian GNU/Linux 11 (64 bit). All the codes are compiled with GCC version 10.2.1, with option -O3. We ran MatchMaker [6, 15] to verify that there were perfect matchings. We used MatchMaker with options “no cheap matching”, “Push-Relabel + fairness” as the core algorithm. For the sake of completeness, we report that the maximum run time of MatchMaker for any matrix was 2.65 seconds for `vas_stokes_4M` and 28.19 for the column permuted version of the same
matrix.

## 4.2 Performance analysis of BOTTLED

We first analyze the percentage of the total run time of BOTTLED spent in the preprocessing step—creating the CSR representation or sorting the entries. Table 2 summarizes the results for six different problem instances, where BOTTLED took more than one second for A-, DAE- and/or DP(A)E-type instances. This table presents the geometric mean of the percentage of the time spent creating the CSR representation and sorting (with respect to the total time of BOTTLED). The row “tTime” contains the geometric mean of the run time of BOTTLED on the instances A, DAE and DP(A)E in seconds, and for the instances AP, DAPE and DP(A)PE it contains the geometric mean of the ratio of the run time of BOTTLED on the permuted instances to the original ones. For example, the run time under the column DP(A)PE is obtained by taking the geometric mean of ratio of the run time of BOTTLED on the 33 DP(A)PE instances to that on the 33 DP(A)E instances.

As we can see from Table 2, the average run time of BOTTLED on any of the six problem instance types is below 10 seconds. A first observation is that for all instances, the two preprocessing steps account for a non-negligible part of the total run time. At two extremes, they take 7% and 37% of the total time in DP(A)E and DAPE instances, respectively. We further observe that the absolute run time of CSR and sorting increases for the permuted instances. The percentage of the total time spent in CSR also increases for the permuted instances while that of sorting either increases or remains the same. If the CSR representation is available, its creation can be skipped and one can reduce the run time considerably. As discussed in Section 3.2, for the targeted BvN application BOTTLED can skip not only the CSR creation, but also the sorting phase across different runs of the algorithm.

<table><caption>■ Table 2 Percentage of the total time spent in creating a CSR matrix and sorting the nonzeros in BOTTLED; “tTime”: the geometric mean of the run time of BOTTLED on the instances A, DAE and DP(A)E in seconds, and the ratios of the others to their counterparts.</caption><thead><tr><th></th><th>A</th><th>AP</th><th>DAE</th><th>DAPE</th><th>DP(A)E</th><th>DP(A)PE</th></tr></thead><tbody><tr><th>CSR</th><td>5%</td><td>13%</td><td>7%</td><td>25%</td><td>4%</td><td>17%</td></tr><tr><th>sort</th><td>5%</td><td>5%</td><td>10%</td><td>12%</td><td>3%</td><td>3%</td></tr><tr><th>tTime</th><td>6.54(s)</td><td>1.28</td><td>3.43(s)</td><td>1.02</td><td>3.40(s)</td><td>1.28</td></tr></tbody></table>

<table><caption>■ Table 3 Breakdown of the total time (tTime) of BOTTLED in seconds.</caption><thead><tr><th>matrix</th><th>instance-type</th><th>tTime</th><th>CSR</th><th>sort</th></tr></thead><tbody><tr><td rowspan="4">vas_stokes_4M</td><td>A</td><td>14.81</td><td>6.70</td><td>3.17</td></tr><tr><td>AP</td><td>32.36</td><td>20.13</td><td>5.51</td></tr><tr><td>DAE</td><td>23.44</td><td>6.93</td><td>3.66</td></tr><tr><td>DAPE</td><td>52.99</td><td>19.93</td><td>5.57</td></tr><tr><td rowspan="2">vas_stokes_2M</td><td>DP(A)E</td><td>11.67</td><td>2.96</td><td>2.01</td></tr><tr><td>DP(A)PE</td><td>19.30</td><td>8.81</td><td>2.57</td></tr></tbody></table>

Another observation from Table 2 is that the run time of BOTTLED consistently increases for the permuted instances. To put this into perspective, we present the run time of BOTTLED on a few instances in Table 3. We see that the increase in run time for the permuted matrices can be attributed, in part, to the creation of the CSR. For example, for vas_stokes_4M, excluding CSR and sort times from the total time yields 4.94 and 6.72 seconds for the while loop for A- and AP-type respectively. A similar calculation for the pairs DP(A)E and DP(A)PE for vas_stokes_2M shows that the while loop takes 6.7 and 7.92 seconds, respectively. In the instance pairs DAE and DAPE for vas_stokes_4M, the increase in the CSR time is still a contributing factor. From these two tables, we conclude that BOTTLED's preprocessing takes up a significant portion of the total run time.

We have also investigated the number of iterations of the while-loop of BOTTLED in different types of instances to see how stable and robust it is with respect to random column permutations. The number of iterations of the while loop were almost always the same for
all 90 instances; in a few cases there was a difference of one in the number of iterations. The small changes in the number of iterations confirm the robustness of BOTTLED.

## 4.3 Comparison of different algorithms

We compared BOTTLED against MC64J2 and THRESH. We first provide a broad overview of the comparisons in Table 4, which lists the geometric mean of the run time of the three methods as well as their maximum run time in different problem instance types (all instances). As can be seen, BOTTLED exhibits the best performance overall, despite being slightly slower on average for A-, DAE-, DAPE- instance types than MC64J2—the margin is not large enough for MC64J2 to make up for the performance loss on the other instance types. Furthermore, BOTTLED's maximum run time is consistently smaller than those of both MC64J2 and THRESH. In contrast, the maximum run time of MC64J2 is prohibitive in all but the A and DAE-type instances. While THRESH avoids the prohibitive run time of its equivalent MC64J3 (see Table 1) by using a better cardinality matching algorithm [15], its maximum run time is still noticeably larger than that of BOTTLED in all but two cases.

We now look more into the performance of BOTTLED versus that of MC64J2. Figure 1a shows the natural logarithm (in the y-axis) of the ratio of the run time of MC64J2 to that of BOTTLED for the instances on which either method has a run time greater than 1 second. Here, MC64J2 is faster than BOTTLED on 22 instances, and BOTTLED is faster on the remaining 117 instances. In all 139 instances, the geometric mean of the ratio of the run time of MC64J2 to that of BOTTLED is 8.5, confirming that BOTTLED is faster than MC64J2 in average.

<figure><img src="image_2.png" /><figcaption>(a) The natural logarithm of the ratio of the run time of MC64J2 to that of BOTTLED in the y-axis in sorted order on 139 problem instances in the x-axis.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) The run time behavior of MC64J2 on four DP(A)E-type instances measured at 20 uniform steps (x-axis), normalized to the total time (y-axis).</figcaption></figure>

<caption>**Figure 1** Performance comparison between MC64J2 and BOTTLED, and investigation on MC64J2.</caption>

In order to put these numbers into a perspective with the run time, we present Table 5. Table 5 lists six matrices in which the ratio of the run time of MC64J2 to that of BOTTLED was the smallest for certain instance types with or without permutation. It next lists six matrices with different instance types in which the ratio of the run time of MC64J2 to that of BOTTLED was the largest for the original matrices and three for the permuted ones (three others were already in the list). Here we see that MC64J2 can have equivalently good performance on both the A- and AP-type instances (see the first block in Table 5). Still, there were some AP instances where MC64J2 ran prohibitively long (see the last two rows with ss and vas_stokes_4M). MC64J2 also struggled on several DP(A)E and DP(A)PE instances as seen in the second block. The prohibitively high run time of MC64J2 highlight
12 Algorithms for bottleneck matching

<table><caption>Table 4 The geometric mean and the maximum of the run time of the three methods on six different instance types.</caption><thead><tr><th></th><th></th><th>A</th><th>AP</th><th>DAE</th><th>DAPE</th><th>DP(A)E</th><th>DP(A)PE</th></tr></thead><tbody><tr><td rowspan="3">geomean</td><td>MC64J2</td><td>0.49</td><td>6.23</td><td>0.16</td><td>0.74</td><td>12.47</td><td>20.16</td></tr><tr><td>THRESH</td><td>2.58</td><td>3.80</td><td>1.09</td><td>1.54</td><td>2.24</td><td>3.16</td></tr><tr><td>BOTTLED</td><td>1.23</td><td>2.38</td><td>0.44</td><td>0.81</td><td>1.23</td><td>2.07</td></tr><tr><td rowspan="3">maximum</td><td>MC64J2</td><td>123.04</td><td>18054.00</td><td>31.42</td><td>9139.63</td><td>1419.99</td><td>6813.20</td></tr><tr><td>THRESH</td><td>43.55</td><td>64.90</td><td>59.42</td><td>91.74</td><td>63.51</td><td>77.57</td></tr><tr><td>BOTTLED</td><td>14.52</td><td>31.81</td><td>23.22</td><td>51.55</td><td>58.64</td><td>76.20</td></tr></tbody></table>

the issue in the approaches based on the shortest augmenting paths: one may visit many edges and vertices to find augmenting paths at different stages during the execution. As such a behavior is also seen for DFS-based-cardinality matching algorithms [6], it cannot be attributed on the particular implementation of the shortest path algorithm in MC64J2.

Figure 1b shows the run time behavior of MC64J2 on the four DP(A)E-type instances of Table 1. The total number of augmentations is divided by 20 to obtain a step size, the run time is measured after each step and normalized by the total time. As this figure shows, in some instances the augmentations took about the same time all throughout, whereas in others later augmentations took more time. As the greedy matching approach of MC64J2 does not always find a maximum matching, it can needlessly result in many additional augmentations. Those augmentations may lead to large run time, not solely because of their number. For example, in the DP(A)E-type instance of CurlCurl_3, a square matrix with $n = 1219574$ rows and $nnz = 13544618$ nonzeros, there are 20040 augmentations with a total run time of about 600 seconds, even though the maximum cardinality matching on $G[\omega_0]$ has a deficiency of one. Obviously, detecting this would lead to much better run time. On the other hand, for the DP(A)E-type instance of rajat31, a square matrix with $n = 4690002$ rows and $nnz = 20316253$ nonzeros, MC64J2 needs 1562500 augmentations (its greedy approach finds a maximum cardinality matching initially), and the whole run time is 48.39 seconds. As rajat31 is much bigger than CurlCurl_3 and needs more augmentations, its shorter run time attests that the process of augmenting one-by-one can take large time due to instance specific properties. It hence cannot constitute a reliable method for the bottleneck matching problem in general. In passing we note that MC64J2 is well-engineered and is faster than using our own SAP implementation for the augmentations.

We compare now BOTTLED and THRESH. In total, BOTTLED beat THRESH in 165 instances out of 180. The largest differences (in seconds) in favor of THRESH were in the A- and AP-type instances of c-73b and c-73. THRESH obtained a run time of 2.88 and 3.48 seconds on c-73b, and 6.35 and 5.11 seconds on c-73. On these instances BOTTLED was relatively close to THRESH (see Table 5)—the largest difference is 10.67 seconds on the A-type instance of c-73b. Since both methods utilize the same core matching algorithm, the superiority of BOTTLED over THRESH should come from doing fewer iterations. Figure 2a supports this reasoning by plotting the difference between the number of iterations of THRESH and that of BOTTLED in nondecreasing order. As seen in this figure, THRESH's number of iterations is always larger than BOTTLED's; the theoretical observations of Section 3 translate to practical gains. We compute the ratio of the run time of THRESH to that of BOTTLED for the instances on which either THRESH or BOTTLED has a run time greater than 1 second (123 instances), and present the natural logarithm of this ratio in the y-axis in Figure 2b. As seen in this figure, BOTTLED fares better than THRESH in the majority of cases.
<table><caption>Table 5 The run time of MC64J2 and BOTTLED on selected instances, where their ratios were the lowest or the highest in different problem instance types or the permuted versions.</caption><thead><tr><th>matrix</th><th>type</th><th>MC64J2</th><th>BOTTLED</th><th>type</th><th>MC64J2</th><th>BOTTLED</th></tr></thead><tbody><tr><td>c-73b</td><td>A</td><td>0.09</td><td>13.55</td><td>AP</td><td>0.10</td><td>9.42</td></tr><tr><td>c-73</td><td>A</td><td>0.19</td><td>14.52</td><td>AP</td><td>0.34</td><td>14.68</td></tr><tr><td>dielFilterV3real</td><td>A</td><td>0.33</td><td>11.72</td><td>AP</td><td>0.81</td><td>21.48</td></tr><tr><td>boyd2</td><td>DP(A)E</td><td>0.03</td><td>3.31</td><td>DP(A)PE</td><td>0.07</td><td>1.84</td></tr><tr><td>boyd2</td><td>DAE</td><td>0.03</td><td>2.43</td><td>DAPE</td><td>0.10</td><td>1.50</td></tr><tr><td>dielFilterV3clx</td><td>A</td><td>0.12</td><td>3.97</td><td>AP</td><td>0.49</td><td>7.19</td></tr><tr><td>atmosmodd</td><td>DP(A)E</td><td>555.76</td><td>7.01</td><td>DP(A)PE</td><td>1757.81</td><td>12.11</td></tr><tr><td>vas_stokes_1M</td><td>DP(A)E</td><td>452.36</td><td>4.20</td><td>DP(A)PE</td><td>2427.71</td><td>7.42</td></tr><tr><td>vas_stokes_2M</td><td>DP(A)E</td><td>1419.99</td><td>11.93</td><td>DP(A)PE</td><td>6813.20</td><td>18.37</td></tr><tr><td>ss</td><td>DP(A)E</td><td>824.91</td><td>4.66</td><td>DP(A)PE</td><td>2621.10</td><td>10.79</td></tr><tr><td>CurlCurl_3</td><td>DP(A)E</td><td>610.16</td><td>1.54</td><td>DP(A)PE</td><td>211.56</td><td>3.80</td></tr><tr><td>atmosmodj</td><td>DP(A)E</td><td>565.60</td><td>7.27</td><td>DP(A)PE</td><td>1683.32</td><td>11.97</td></tr><tr><td>vas_stokes_4M</td><td>DAE</td><td>0.83</td><td>23.22</td><td>DAPE</td><td>9139.63</td><td>51.55</td></tr><tr><td>ss</td><td>A</td><td>2.25</td><td>2.65</td><td>AP</td><td>2532.09</td><td>8.65</td></tr><tr><td>vas_stokes_4M</td><td>A</td><td>0.53</td><td>14.23</td><td>AP</td><td>18054.00</td><td>31.81</td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>

(a) The number of iterations of THRESH minus that of BOTTLED in the y-axis in sorted order on all 180 instances.

<figure><img src="image_5.png" /></figure>

(b) The natural logarithm of the ratio of the run time of THRESH to that of BOTTLED in the y-axis in sorted order on 123 instances in the x-axis.

<caption>Figure 2 Performance comparison between THRESH and BOTTLED.</caption>

## 4.4 Inside a BvN decomposition method

Table 6 presents the run time of the BvN decomposition method on the DAE- and DP(A)E- type instances of four matrices (obtained with 2500 scaling iterations). We run the BvN decomposition method until 50 permutation matrices or a coefficient of 0.92 are obtained. As each permutation matrix is obtained by a call to BOTTLED, we show their number in the column “num.perm”. Further, the table also presents the time for the initial preprocessing of BOTTLED, and the maximum time taken in a call to BOTTLED subsequently for obtaining a permutation matrix. We observe that the maximum run time of BOTTLED, in an iteration, multiplied by the number of calls is at least 2.92 (for DAE of ss) and up to 10.67 (for DP(A)E of ss) times the total BvN time. This suggests that the run time of the subsequent bottleneck matching calls reduces appreciably, and BOTTLED works well inside the decomposition method by avoiding the preprocessing.
14 Algorithms for bottleneck matching

<table><caption>Table 6 Run time, in seconds, of BOTTLED and the BvN decomposition heuristic, along with the number of permutation matrices in the decomposition.</caption><thead><tr><th rowspan="2">matrix</th><th rowspan="2">instance</th><th rowspan="2">BOTTLED<br>initialization</th><th rowspan="2">longest run time of<br>BOTTLED in an iteration</th><th colspan="2">BvN</th></tr><tr><th>num.perm</th><th>time</th></tr></thead><tbody><tr><td rowspan="2">atmosmodm</td><td>DAE</td><td>0.22</td><td>58.73</td><td>46</td><td>452.14</td></tr><tr><td>DP(A)E</td><td>0.21</td><td>14.43</td><td>50</td><td>200.80</td></tr><tr><td rowspan="2">CurlCurl_3</td><td>DAE</td><td>0.35</td><td>17.30</td><td>50</td><td>183.43</td></tr><tr><td>DP(A)E</td><td>0.34</td><td>4.33</td><td>50</td><td>68.45</td></tr><tr><td rowspan="2">ss</td><td>DAE</td><td>1.15</td><td>13.37</td><td>50</td><td>228.92</td></tr><tr><td>DP(A)E</td><td>1.46</td><td>39.77</td><td>50</td><td>186.31</td></tr><tr><td rowspan="2">vas_stokes_2M</td><td>DAE</td><td>5.04</td><td>4.74</td><td>50</td><td>78.32</td></tr><tr><td>DP(A)E</td><td>5.10</td><td>4.87</td><td>50</td><td>74.36</td></tr></tbody></table>

5 Conclusion

We have investigated the problem of finding a maximum bottleneck matching in bipartite graphs. Existing implementations for the problem suffer from unpredictable run time that can get prohibitively large, i.e., requiring thousands or even tens of thousands of seconds to complete. We have proposed a new algorithm called BOTTLED that converts an inefficient, duality-based approach into an efficient one through theoretical findings. Experimental results show that BOTTLED is almost always faster than the state-of-the-art methods. Furthermore, its run time is reliable and always remains within reasonable time limits. We have also explored its use inside a heuristic for the Birkhoff-von Neumann decomposition of doubly stochastic matrices and experimentally established the suitability of the proposed algorithm for this purpose.

Currently the proposed approach resorts to an augmenting-path-based method in few corner cases and only when there are perfect matchings. We plan to explore the possibility to use them more effectively, along with potential data reduction rules.

References

1 M. Benzi and B. Uçar. Preconditioning techniques based on the Birkhoff-von Neumann decomposition. *Computational Methods in Applied Mathematics*, 17:201-215, 2017.

2 G. Birkhoff. Tres observaciones sobre el algebra lineal. *Univ. Nac. Tucumán Rev. Ser. A*, (5):147-150, 1946.

3 R. Burkard, M. Dell'Amico, and S. Martello. *Assignment Problems*. SIAM, Philadelphia, PA, USA, 2009.

4 C. Chang, W. Chen, and H. Huang. On service guarantees for input-buffered crossbar switches: A capacity decomposition approach by Birkhoff and von Neumann. In *Quality of Service, 1999. IWQoS '99. 1999 Seventh International Workshop on*, pages 79-86, 1999.

5 T. A. Davis and Y. Hu. The University of Florida sparse matrix collection. *ACM Trans. Math. Softw.*, 38(1):1:1-1:25, 2011.

6 I. S. Duff, K. Kaya, and B. Uçar. Design, implementation, and analysis of maximum transversal algorithms. *ACM Transactions on Mathematical Software*, 38:13:1-13:31, 2011.

7 I. S. Duff and J. Koster. The design and use of algorithms for permuting large entries to the diagonal of sparse matrices. *SIAM Journal on Matrix Analysis and Applications*, 20(4):889-901, 1999.

8 I. S. Duff and J. Koster. On algorithms for permuting large entries to the diagonal of a sparse matrix. *SIAM Journal on Matrix Analysis and Applications*, 22:973-996, 2001. doi:http://dx.doi.org/10.1137/S0895479899358443.
9 F. Dufossé, K. Kaya, I. Panagiotas, and B. Uçar. Further notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices. *Linear Algebra and its Applications*, 554:68-78, 2018. URL: http://www.sciencedirect.com/science/article/pii/S0024379518302568, doi:https://doi.org/10.1016/j.laa.2018.05.017.

10 F. Dufossé and B. Uçar. Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices. *Linear Algebra and its Applications*, 497:108-115, 2016. URL: http://www.sciencedirect.com/science/article/pii/S0024379516001257, doi:http://dx.doi.org/10.1016/j.laa.2016.02.023.

11 A. L. Dulmage and N. S. Mendelsohn. Coverings of bipartite graphs. *Canadian Journal of Mathematics*, 10:517-534, 1958.

12 H. N. Gabow and R. E. Tarjan. Algorithms for two bottleneck optimization problems. *J. Algorithms*, 9(3):411-417, 1988. Hopcroft-Karp algorithm can be used as an approximation algorithm p415.

13 P. Hall. On representatives of subsets. *Journal of the London Mathematical Society*, s1-10(37):26-30, 1935. arXiv: http://jlms.oxfordjournals.org/cgi/reprint/s1-10/37/26.pdf.

14 J. E. Hopcroft and R. M. Karp. An $n^{5/2}$ algorithm for maximum matchings in bipartite graphs. *SIAM Journal on Computing*, 2(4):225-231, 1973.

15 K. Kaya, J. Langguth, F. Manne, and B. Uçar. Push-relabel based algorithms for the maximum transversal problem. *Computers & Operations Research*, 40(5):1266-1275, 2013.

16 J. Kulkarni, E. Lee, and M. Singh. Minimum Birkhoff-von Neumann decomposition. In *International Conference on Integer Programming and Combinatorial Optimization*, pages 343-354. Springer, 2017.

17 I. Panagiotas, G. Pichon, S. Singh, and B. Uçar. Bottled: Fast algorithms for the bottleneck matching problem, April 2023. doi:10.5281/zenodo.7871464.

18 A. Pothen and C.-J. Fan. Computing the block triangular form of a sparse matrix. *ACM Transactions on Mathematical Software*, 16(4):303-324, December 1990. URL: http://doi.acm.org/10.1145/98267.98287, doi:10.1145/98267.98287.

19 A. P. Punnen and K. Nair. Improved complexity bound for the maximum cardinality bottleneck bipartite matching problem. *Discret. Appl. Math.*, 55(1):91-93, 1994.

20 R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic matrices. *Pacific J. Math.*, 21:343-348, 1967.

21 V. Valls, G. Iosifidis, and L. Tassiulas. Birkhoff's decomposition revisited: Sparse scheduling for high-speed circuit switches. *IEEE/ACM Transactions on Networking*, 29(6):2399-2412, 2021. doi:10.1109/TNET.2021.3088327.
A Description of the appendix

A detailed description of the data set is given in Appendix B. This same section also contains the whole set of experiments with MC64J2, THRESH, and BOTTLED on all the data set. Appendix C contains an investigation of the parts of the proposed BOTTLED algorithm, where we present results to assess the effects of using the shortest-augmenting path methods.

B Detailed tables

We give comprehensive results for the presented experiments. We start with the properties of matrices used in experiments in Table 7 for the A-type problem instances, in Table 8 for the DAE-type problem instances, and in Table 10 for the DP(A)E-type problem instances. Recall that there are 14 A-type, 18 DAE-type, 58 DP(A)E-type instances, and the same number of instances with column permutations, as we discarded the instances in which the initialization algorithm found the bottleneck value. In these three tables, we give the number of rows/columns and the number of nonzeros of the matrices; some are repeated in different tables for the sake of convenience. We then give information about the graph G[ω₀], where ω₀ is the value obtained in the initialization (the minimum of maximum in each row and in each column): (i) the number of edges in G[ω₀]; (ii) the deficiency obtained by the fast greedy algorithm of MC64J2 under the header “def MC64J2 ”; and (iii) the difference between n and the maximum cardinality of a matching in G[ω₀]. The following three tables Table 11, Table 12 and Table 13 give the run time of MC64J2, THRESH, and BOTTLED on, respectively, (i) A- and AP-type instances; (ii) DAE- and DAPE-type instances; and (iii) DP(A)E- and DP(A)PE-type instances.

<table><caption>Table 7 n × n matrices with nnz nonzeros, the number of edges in the graph G[ω₀], where ω₀ is the first bottleneck value of the A-type instances, the deficiency “def MC64J2 ” of the initialization of MC64J2, and the deficiency “def.” of the maximum cardinality of a matching in G[ω₀] with respect to that of G.</caption><thead><tr><th rowspan="2">matrix</th><th rowspan="2">n</th><th rowspan="2">nnz</th><th colspan="3">G[ω<sub>0</sub>]</th></tr><tr><th>#Edges</th><th>def<sub>MC64J2</sub></th><th>def.</th></tr></thead><tbody><tr><td>boyd2</td><td>466316</td><td>1500397</td><td>1220614</td><td>93252</td><td>93252</td></tr><tr><td>c-73</td><td>169422</td><td>1279274</td><td>701684</td><td>3440</td><td>3410</td></tr><tr><td>c-73b</td><td>169422</td><td>1279274</td><td>409370</td><td>3004</td><td>3004</td></tr><tr><td>cont-300</td><td>180895</td><td>988195</td><td>896402</td><td>453</td><td>299</td></tr><tr><td>d_pretok</td><td>182730</td><td>1641672</td><td>638067</td><td>52425</td><td>52419</td></tr><tr><td>dielFilterV3clx</td><td>420408</td><td>32886208</td><td>18823657</td><td>7</td><td>5</td></tr><tr><td>dielFilterV3real</td><td>1102824</td><td>89306020</td><td>54707197</td><td>35</td><td>18</td></tr><tr><td>Goodwin_095</td><td>100037</td><td>3226066</td><td>2705970</td><td>94</td><td>94</td></tr><tr><td>Goodwin_127</td><td>178437</td><td>5778545</td><td>4848818</td><td>126</td><td>126</td></tr><tr><td>Hamrle3</td><td>1447360</td><td>5514242</td><td>3664040</td><td>214288</td><td>202281</td></tr><tr><td>PR02R</td><td>161070</td><td>8185136</td><td>3410413</td><td>481</td><td>11</td></tr><tr><td>ss</td><td>1652680</td><td>34753577</td><td>9288800</td><td>24964</td><td>24964</td></tr><tr><td>turon_m</td><td>189924</td><td>1690876</td><td>546760</td><td>72826</td><td>72786</td></tr><tr><td>vas_stokes_4M</td><td>4382246</td><td>131577616</td><td>94445792</td><td>4</td><td>4</td></tr></tbody></table>

C Investigating the effect of resorting to SAP in BOTTLED

We briefly investigate the effect of resorting to the SAP subroutine, and the effects of the vertex selection heuristic (Line 4 of Algorithm 2). By design, BOTTLED resorts to SAP on rare cases; indeed it takes place only in 2 of the A-type, 8 of the DAE-type, and 16 of the DP(A)E-type instances. We present results with six instances in Table 9, where we give the
<table><thead><tr><th rowspan="2">matrix</th><th rowspan="2">n</th><th rowspan="2">nnz</th><th colspan="3">G[ω<sub>0</sub>]</th></tr><tr><th>#Edges</th><th>def.<sub>MC64J2</sub></th><th>def.</th></tr></thead><tbody><tr><td>analytics</td><td>303813</td><td>2006126</td><td>824483</td><td>20</td><td>3</td></tr><tr><td>boyd2</td><td>466316</td><td>1500397</td><td>1278564</td><td>10</td><td>2</td></tr><tr><td>c-73</td><td>169422</td><td>1279274</td><td>882004</td><td>383</td><td>1</td></tr><tr><td>crashbasis</td><td>160000</td><td>1750416</td><td>163658</td><td>2</td><td>2</td></tr><tr><td>d_pretok</td><td>182730</td><td>1641672</td><td>299061</td><td>4103</td><td>3917</td></tr><tr><td>Hamrle3</td><td>1447360</td><td>5514242</td><td>2569988</td><td>18499</td><td>17456</td></tr><tr><td>kim2</td><td>456976</td><td>11330020</td><td>1687449</td><td>2643</td><td>2628</td></tr><tr><td>lung2</td><td>109460</td><td>492564</td><td>110229</td><td>2</td><td>2</td></tr><tr><td>mac_econ_fwd500</td><td>206500</td><td>1273389</td><td>469290</td><td>11094</td><td>4413</td></tr><tr><td>mc2depi</td><td>525825</td><td>2100225</td><td>656454</td><td>1</td><td>1</td></tr><tr><td>PR02R</td><td>161070</td><td>8185136</td><td>514760</td><td>966</td><td>59</td></tr><tr><td>RM07R</td><td>381689</td><td>37464962</td><td>1706225</td><td>4231</td><td>7</td></tr><tr><td>stomach</td><td>213360</td><td>3021648</td><td>220920</td><td>35</td><td>15</td></tr><tr><td>torso1</td><td>116158</td><td>8516500</td><td>1023697</td><td>1423</td><td>1024</td></tr><tr><td>torso3</td><td>259156</td><td>4429042</td><td>263564</td><td>1</td><td>1</td></tr><tr><td>TSOPF_FS_b39_c30</td><td>120216</td><td>3121160</td><td>390053</td><td>825</td><td>9</td></tr><tr><td>turon_m</td><td>189924</td><td>1690876</td><td>399669</td><td>781</td><td>474</td></tr><tr><td>vas_stokes_4M</td><td>4382246</td><td>131577616</td><td>18666723</td><td>172</td><td>170</td></tr></tbody></table>

■ Table 8 n × n matrices with nnz nonzeros, the number of edges in the graph G[ω₀], where ω₀ is the first bottleneck value of the DAE-type instances, the deficiency “def MC64J2” of the initialization of MC64J2, and the deficiency “def.” of the maximum cardinality of a matching in G[ω₀] with respect to that of G.

number of iterations of the while-loop when SAP is not called (no aug.), when the vertex selection is not applied within SAP (no. vs, a vertex with the largest value is selected instead), and with the presented version of BOTTLED. As seen in Table 9, resorting to SAP in general helps, and the vertex selection heuristic makes BOTTLED more robust; for example in c-73 not using SAP results in 63 additional iterations, and in lung2 not selecting the vertex with the minimum value results in 46 additional iterations.

<table><thead><tr><th rowspan="2">matrix</th><th rowspan="2">instance</th><th rowspan="2">no aug.</th><th colspan="2">with aug.</th></tr><tr><th>no. vs</th><th>BOTTLED</th></tr></thead><tbody><tr><td>c-73b</td><td><b>A</b></td><td>38</td><td>10</td><td>10</td></tr><tr><td>c-73</td><td><b>DAE</b></td><td>66</td><td>3</td><td>3</td></tr><tr><td>Hamrle3</td><td><b>DP(A)E</b></td><td>22</td><td>11</td><td>11</td></tr><tr><td>lung2</td><td><b>DP(A)E</b></td><td>15</td><td>55</td><td>9</td></tr><tr><td>PR02R</td><td><b>DAE</b></td><td>19</td><td>14</td><td>16</td></tr><tr><td>cont-300</td><td><b>DP(A)E</b></td><td>8</td><td>10</td><td>10</td></tr></tbody></table>

■ Table 9 The number of iterations of the while loop, without augmenting paths, with augmenting paths but without vertex selection heuristic, and the presented version of BOTTLED.
18 Algorithms for bottleneck matching

<table><caption>Table 10 $n \times n$ matrices with nnz nonzeros, the number of edges in the graph $G[\omega_0]$, where $\omega_0$ is the first bottleneck value of the DP(A)E-type instances, the deficiency “def.$_{MC64J2}$” of the initialization of MC64J2, and the deficiency “def.” of the maximum cardinality of a matching in $G[\omega_0]$ with respect to that of G.</caption><thead><tr><th rowspan="2">matrix</th><th rowspan="2">n</th><th rowspan="2">nnz</th><th colspan="3">G[ω<sub>0</sub>]</th></tr><tr><th>#Edges</th><th>def<sub>MC64J2</sub></th><th>def.</th></tr></thead><tbody><tr><td>af_0_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>af_1_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>af_2_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>af_3_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>af_4_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>af_5_k101</td><td>503625</td><td>17550675</td><td>12175825</td><td>7939</td><td>6520</td></tr><tr><td>analytics</td><td>303813</td><td>2006126</td><td>1575964</td><td>665</td><td>1</td></tr><tr><td>atmosmodd</td><td>1270432</td><td>8814880</td><td>3469408</td><td>60277</td><td>59496</td></tr><tr><td>atmosmodj</td><td>1270432</td><td>8814880</td><td>3469408</td><td>60277</td><td>59496</td></tr><tr><td>atmosmodl</td><td>1489752</td><td>10319760</td><td>2836248</td><td>86407</td><td>85376</td></tr><tr><td>atmosmodm</td><td>1489752</td><td>10319760</td><td>2836248</td><td>86407</td><td>85376</td></tr><tr><td>BenElechi1</td><td>245874</td><td>13150496</td><td>9301880</td><td>7097</td><td>4470</td></tr><tr><td>boyd2</td><td>466316</td><td>1500397</td><td>1313873</td><td>10</td><td>2</td></tr><tr><td>Chevron3</td><td>381381</td><td>3413113</td><td>2849503</td><td>2695</td><td>2672</td></tr><tr><td>Chevron4</td><td>711450</td><td>6376412</td><td>5586092</td><td>3741</td><td>3718</td></tr><tr><td>circuit5M</td><td>5558326</td><td>59524291</td><td>44231380</td><td>1448</td><td>3</td></tr><tr><td>CO</td><td>221119</td><td>7666057</td><td>7620031</td><td>4566</td><td>1</td></tr><tr><td>cont-300</td><td>180895</td><td>988195</td><td>813457</td><td>2239</td><td>1134</td></tr><tr><td>crashbasis</td><td>160000</td><td>1750416</td><td>1292128</td><td>1490</td><td>1456</td></tr><tr><td>CurlCurl_3</td><td>1219574</td><td>13544618</td><td>12157436</td><td>20040</td><td>1</td></tr><tr><td>d_pretok</td><td>182730</td><td>1641672</td><td>516597</td><td>75581</td><td>75581</td></tr><tr><td>dc1</td><td>116835</td><td>766396</td><td>485343</td><td>1</td><td>1</td></tr><tr><td>dc2</td><td>116835</td><td>766396</td><td>485343</td><td>1</td><td>1</td></tr><tr><td>dc3</td><td>116835</td><td>766396</td><td>485343</td><td>1</td><td>1</td></tr><tr><td>ecology1</td><td>1000000</td><td>4996000</td><td>4623256</td><td>3947</td><td>3944</td></tr><tr><td>ecology2</td><td>999999</td><td>4995991</td><td>4616887</td><td>4067</td><td>4064</td></tr><tr><td>FEM_3D_thermal2</td><td>147900</td><td>3489300</td><td>1890900</td><td>1205</td><td>1200</td></tr><tr><td>Ga10As10H30</td><td>113081</td><td>6115633</td><td>3692549</td><td>2026</td><td>2</td></tr><tr><td>Ga19As19H42</td><td>133123</td><td>8884839</td><td>4261426</td><td>1792</td><td>5</td></tr><tr><td>Ge99H100</td><td>112985</td><td>8451395</td><td>3395058</td><td>4606</td><td>67</td></tr><tr><td>Goodwin_095</td><td>100037</td><td>3226066</td><td>1001631</td><td>8502</td><td>8210</td></tr><tr><td>Goodwin_127</td><td>178437</td><td>5778545</td><td>1841347</td><td>15260</td><td>14988</td></tr><tr><td>Hamrle3</td><td>1447360</td><td>5514242</td><td>2581382</td><td>123188</td><td>122750</td></tr><tr><td>Hardesty1</td><td>938905</td><td>12143314</td><td>10661294</td><td>6394</td><td>4922</td></tr><tr><td>hcircuit</td><td>105676</td><td>513072</td><td>506361</td><td>1</td><td>1</td></tr><tr><td>iChem_Jacobian</td><td>274087</td><td>4137369</td><td>1716099</td><td>96864</td><td>88302</td></tr><tr><td>kim2</td><td>456976</td><td>11330020</td><td>8963433</td><td>3974</td><td>2604</td></tr><tr><td>Lin</td><td>256000</td><td>1766400</td><td>541648</td><td>15933</td><td>15576</td></tr><tr><td>lung2</td><td>109460</td><td>492564</td><td>267598</td><td>887</td><td>176</td></tr><tr><td>mac_econ_fwd500</td><td>206500</td><td>1273389</td><td>702581</td><td>12009</td><td>1005</td></tr><tr><td>majorbasis</td><td>160000</td><td>1750416</td><td>1292128</td><td>1490</td><td>1456</td></tr><tr><td>mc2depi</td><td>525825</td><td>2100225</td><td>1939125</td><td>3023</td><td>3016</td></tr><tr><td>parabolic_fem</td><td>525825</td><td>3674625</td><td>3288209</td><td>8798</td><td>2996</td></tr><tr><td>PR02R</td><td>161070</td><td>8185136</td><td>3596918</td><td>38349</td><td>37621</td></tr><tr><td>rajat31</td><td>4690002</td><td>20316253</td><td>7815002</td><td>1562500</td><td>1562500</td></tr><tr><td>RM07R</td><td>381689</td><td>37464962</td><td>35794471</td><td>5406</td><td>1555</td></tr><tr><td>Si41Ge41H72</td><td>185639</td><td>15011265</td><td>5791343</td><td>4724</td><td>41</td></tr><tr><td>ss</td><td>1652680</td><td>34753577</td><td>21786247</td><td>32086</td><td>1</td></tr><tr><td>ss1</td><td>205282</td><td>845089</td><td>590921</td><td>7862</td><td>3175</td></tr><tr><td>stomach</td><td>213360</td><td>3021648</td><td>1712993</td><td>42033</td><td>41776</td></tr><tr><td>t2em</td><td>921632</td><td>4590832</td><td>4149720</td><td>4512</td><td>4512</td></tr><tr><td>tmt_unsym</td><td>917825</td><td>4584801</td><td>4278035</td><td>2986</td><td>2599</td></tr><tr><td>torso1</td><td>116158</td><td>8516500</td><td>1840272</td><td>1799</td><td>1791</td></tr><tr><td>torso2</td><td>115967</td><td>1033473</td><td>781379</td><td>3861</td><td>2618</td></tr><tr><td>turon_m</td><td>189924</td><td>1690876</td><td>534458</td><td>77704</td><td>77704</td></tr><tr><td>vas_stokes_1M</td><td>1090664</td><td>34767207</td><td>23969241</td><td>3175</td><td>6</td></tr><tr><td>vas_stokes_2M</td><td>2146677</td><td>65129037</td><td>47869836</td><td>5686</td><td>71</td></tr><tr><td>xenon2</td><td>157464</td><td>3866688</td><td>1945788</td><td>28313</td><td>25849</td></tr></tbody></table>
<table><caption>Table 11 Run time (in seconds) of MC64J2, THRESH, and BOTTLED on all A- and AP-type instances. The geometric mean of the run time of MC64J2, THRESH, and BOTTLED over all A- and AP-type instances are 1.74 3.13 and 1.71, respectively.</caption><thead><tr><th rowspan="2">matrix</th><th colspan="3">A-type</th><th colspan="3">AP-type</th></tr><tr><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th></tr></thead><tbody><tr><td>boyd2</td><td>0.02</td><td>0.25</td><td>0.08</td><td>0.11</td><td>0.59</td><td>0.28</td></tr><tr><td>c-73</td><td>0.19</td><td>6.35</td><td>14.52</td><td>0.34</td><td>5.11</td><td>14.68</td></tr><tr><td>c-73b</td><td>0.09</td><td>2.88</td><td>13.55</td><td>0.10</td><td>3.48</td><td>9.42</td></tr><tr><td>cont-300</td><td>1.51</td><td>0.31</td><td>0.10</td><td>3.89</td><td>0.45</td><td>0.16</td></tr><tr><td>d_pretok</td><td>7.34</td><td>0.46</td><td>0.38</td><td>3.91</td><td>0.70</td><td>0.58</td></tr><tr><td>dielFilterV3clx</td><td>0.12</td><td>9.66</td><td>3.97</td><td>0.49</td><td>13.03</td><td>7.19</td></tr><tr><td>dielFilterV3real</td><td>0.33</td><td>28.04</td><td>11.72</td><td>0.81</td><td>38.78</td><td>21.48</td></tr><tr><td>Goodwin_095</td><td>0.02</td><td>0.61</td><td>0.12</td><td>4.76</td><td>0.90</td><td>0.29</td></tr><tr><td>Goodwin_127</td><td>0.03</td><td>1.09</td><td>0.25</td><td>21.65</td><td>2.01</td><td>0.75</td></tr><tr><td>Hamrle3</td><td>123.04</td><td>7.85</td><td>1.62</td><td>13.86</td><td>14.34</td><td>5.31</td></tr><tr><td>PR02R</td><td>3.59</td><td>2.85</td><td>0.84</td><td>42.93</td><td>3.63</td><td>1.73</td></tr><tr><td>ss</td><td>2.25</td><td>10.23</td><td>2.65</td><td>2532.09</td><td>17.17</td><td>8.65</td></tr><tr><td>turon_m</td><td>0.89</td><td>0.49</td><td>0.41</td><td>2.10</td><td>0.75</td><td>0.61</td></tr><tr><td>vas_stokes_4M</td><td>0.53</td><td>43.55</td><td>14.23</td><td>18054.00</td><td>64.90</td><td>31.81</td></tr><tr><td>geomean</td><td>0.49</td><td>2.58</td><td>1.23</td><td>6.23</td><td>3.80</td><td>2.38</td></tr></tbody></table>

<table><caption>Table 12 Run time (in seconds) of MC64J2, THRESH, and BOTTLED on all DAE- and DAPE-type instances. The geometric mean of the run time of MC64J2, THRESH, and BOTTLED over all DAE- and DAPE-type instances are 0.35, 1.30 and 0.59 respectively.</caption><thead><tr><th rowspan="2">matrix</th><th colspan="3">DAE-type</th><th colspan="3">DAPE-type</th></tr><tr><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th></tr></thead><tbody><tr><td>analytics</td><td>0.01</td><td>0.42</td><td>0.16</td><td>0.06</td><td>0.52</td><td>0.21</td></tr><tr><td>boyd2</td><td>0.03</td><td>2.37</td><td>2.43</td><td>0.10</td><td>1.63</td><td>1.50</td></tr><tr><td>c-73</td><td>0.08</td><td>0.30</td><td>0.12</td><td>0.15</td><td>0.38</td><td>0.15</td></tr><tr><td>crashbasis</td><td>0.01</td><td>0.27</td><td>0.06</td><td>0.02</td><td>0.54</td><td>0.24</td></tr><tr><td>d_pretok</td><td>0.09</td><td>0.43</td><td>0.35</td><td>0.40</td><td>0.60</td><td>0.60</td></tr><tr><td>Hamrle3</td><td>0.10</td><td>2.64</td><td>0.95</td><td>0.46</td><td>4.26</td><td>2.08</td></tr><tr><td>kim2</td><td>31.42</td><td>4.79</td><td>1.81</td><td>37.35</td><td>7.18</td><td>3.76</td></tr><tr><td>lung2</td><td>0.01</td><td>0.08</td><td>0.02</td><td>0.01</td><td>0.12</td><td>0.04</td></tr><tr><td>mac_econfwd500</td><td>1.52</td><td>0.75</td><td>0.64</td><td>3.94</td><td>1.15</td><td>0.97</td></tr><tr><td>mc2depi</td><td>0.02</td><td>0.35</td><td>0.10</td><td>0.06</td><td>0.74</td><td>0.31</td></tr><tr><td>PR02R</td><td>1.80</td><td>2.49</td><td>1.43</td><td>6.70</td><td>3.13</td><td>2.15</td></tr><tr><td>RM07R</td><td>16.45</td><td>11.25</td><td>3.98</td><td>43.19</td><td>14.65</td><td>6.89</td></tr><tr><td>stomach</td><td>0.01</td><td>0.66</td><td>0.13</td><td>0.03</td><td>0.92</td><td>0.37</td></tr><tr><td>torso1</td><td>10.76</td><td>1.94</td><td>0.85</td><td>25.14</td><td>2.37</td><td>1.15</td></tr><tr><td>torso3</td><td>0.02</td><td>1.03</td><td>0.20</td><td>0.05</td><td>1.61</td><td>0.70</td></tr><tr><td>TSOPF_FS_b39_c30</td><td>0.66</td><td>0.59</td><td>0.24</td><td>0.21</td><td>0.84</td><td>0.43</td></tr><tr><td>turon_m</td><td>0.01</td><td>0.39</td><td>0.21</td><td>0.73</td><td>0.59</td><td>0.38</td></tr><tr><td>vas_stokes_4M</td><td>0.83</td><td>59.42</td><td>23.22</td><td>9139.63</td><td>91.74</td><td>51.55</td></tr><tr><td>geomean</td><td>0.16</td><td>1.09</td><td>0.44</td><td>0.74</td><td>1.54</td><td>0.81</td></tr></tbody></table>
20 Algorithms for bottleneck matching

<table><caption>Table 13 Run time (in seconds) of MC64J2, THRESH, and BOTTLED on all DP(A)E- and DP(A)PE-type instances. The geometric mean of the run time of MC64J2, THRESH, and BOTTLED over all DP(A)E- and DP(A)PE-type instances are 15.86, 2.66 and 1.59 respectively.</caption><thead><tr><th rowspan="2">matrix</th><th colspan="3">DP(A)E-type</th><th colspan="3">DP(A)PE-type</th></tr><tr><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th><th>MC64J2</th><th>THRESH</th><th>BOTTLED</th></tr></thead><tbody><tr><td>af_0_k101</td><td>79.91</td><td>8.87</td><td>3.56</td><td>418.20</td><td>13.52</td><td>7.68</td></tr><tr><td>af_1_k101</td><td>80.76</td><td>8.92</td><td>3.59</td><td>426.10</td><td>14.00</td><td>7.72</td></tr><tr><td>af_2_k101</td><td>84.49</td><td>9.38</td><td>3.81</td><td>447.46</td><td>13.94</td><td>8.02</td></tr><tr><td>af_3_k101</td><td>81.85</td><td>9.03</td><td>3.67</td><td>425.69</td><td>13.43</td><td>7.92</td></tr><tr><td>af_4_k101</td><td>80.37</td><td>8.90</td><td>3.59</td><td>428.32</td><td>13.66</td><td>7.73</td></tr><tr><td>af_5_k101</td><td>80.95</td><td>8.95</td><td>3.63</td><td>447.20</td><td>14.13</td><td>7.88</td></tr><tr><td>analytics</td><td>0.11</td><td>0.43</td><td>0.18</td><td>0.15</td><td>0.53</td><td>0.23</td></tr><tr><td>atmosmodd</td><td>555.76</td><td>9.18</td><td>7.01</td><td>1757.81</td><td>13.67</td><td>12.11</td></tr><tr><td>atmosmodj</td><td>565.60</td><td>9.38</td><td>7.27</td><td>1683.32</td><td>14.20</td><td>11.97</td></tr><tr><td>atmosmodl</td><td>621.46</td><td>13.25</td><td>8.70</td><td>1580.69</td><td>20.32</td><td>16.21</td></tr><tr><td>atmosmodm</td><td>699.66</td><td>14.14</td><td>9.38</td><td>1768.57</td><td>22.60</td><td>18.84</td></tr><tr><td>BenElechi1</td><td>31.45</td><td>2.92</td><td>2.53</td><td>20.36</td><td>4.68</td><td>4.56</td></tr><tr><td>boyd2</td><td>0.03</td><td>3.37</td><td>3.31</td><td>0.07</td><td>2.09</td><td>1.84</td></tr><tr><td>Chevron3</td><td>8.38</td><td>2.60</td><td>1.12</td><td>45.28</td><td>3.49</td><td>1.65</td></tr><tr><td>Chevron4</td><td>22.24</td><td>5.11</td><td>1.99</td><td>150.66</td><td>6.84</td><td>4.12</td></tr><tr><td>circuit5M</td><td>4.99</td><td>63.51</td><td>58.64</td><td>57.68</td><td>77.57</td><td>76.20</td></tr><tr><td>CO</td><td>45.63</td><td>1.88</td><td>0.62</td><td>0.36</td><td>2.49</td><td>1.19</td></tr><tr><td>cont-300</td><td>3.17</td><td>0.54</td><td>0.35</td><td>3.93</td><td>0.74</td><td>0.53</td></tr><tr><td>crashbasis</td><td>4.82</td><td>0.73</td><td>0.26</td><td>14.12</td><td>0.90</td><td>0.40</td></tr><tr><td>CurlCurl_3</td><td>610.16</td><td>4.34</td><td>1.54</td><td>211.56</td><td>7.08</td><td>3.80</td></tr><tr><td>d_pretok</td><td>8.36</td><td>1.11</td><td>1.20</td><td>6.33</td><td>1.34</td><td>1.50</td></tr><tr><td>dc1</td><td>0.01</td><td>0.14</td><td>0.06</td><td>0.04</td><td>0.19</td><td>0.10</td></tr><tr><td>dc2</td><td>0.01</td><td>0.14</td><td>0.06</td><td>0.04</td><td>0.19</td><td>0.10</td></tr><tr><td>dc3</td><td>0.01</td><td>0.14</td><td>0.06</td><td>0.03</td><td>0.19</td><td>0.10</td></tr><tr><td>ecology1</td><td>43.10</td><td>3.65</td><td>1.92</td><td>14.36</td><td>5.55</td><td>3.44</td></tr><tr><td>ecology2</td><td>66.06</td><td>2.54</td><td>1.93</td><td>24.35</td><td>4.13</td><td>3.57</td></tr><tr><td>FEM_3D_thermal2</td><td>0.39</td><td>0.83</td><td>0.29</td><td>15.53</td><td>1.26</td><td>0.57</td></tr><tr><td>Ga10As10H30</td><td>18.84</td><td>1.46</td><td>0.55</td><td>0.36</td><td>1.82</td><td>0.75</td></tr><tr><td>Ga19As19H42</td><td>30.42</td><td>2.39</td><td>0.92</td><td>0.98</td><td>2.73</td><td>1.22</td></tr><tr><td>Ge99H100</td><td>25.74</td><td>2.29</td><td>1.00</td><td>15.39</td><td>2.55</td><td>1.21</td></tr><tr><td>Goodwin_095</td><td>3.72</td><td>0.95</td><td>0.83</td><td>8.53</td><td>1.33</td><td>1.18</td></tr><tr><td>Goodwin_127</td><td>6.28</td><td>2.16</td><td>2.08</td><td>50.78</td><td>2.95</td><td>2.81</td></tr><tr><td>Hamrle3</td><td>0.29</td><td>1.74</td><td>1.24</td><td>0.81</td><td>4.42</td><td>3.24</td></tr><tr><td>Hardesty1</td><td>198.57</td><td>6.98</td><td>4.78</td><td>14.75</td><td>11.30</td><td>9.23</td></tr><tr><td>hcircuit</td><td>0.01</td><td>0.07</td><td>0.02</td><td>0.02</td><td>0.11</td><td>0.04</td></tr><tr><td>iChem_Jacobian</td><td>75.27</td><td>1.70</td><td>1.51</td><td>99.14</td><td>2.48</td><td>2.32</td></tr><tr><td>kim2</td><td>112.82</td><td>5.64</td><td>1.61</td><td>345.17</td><td>8.74</td><td>3.75</td></tr><tr><td>Lin</td><td>34.97</td><td>0.92</td><td>0.93</td><td>64.27</td><td>1.21</td><td>1.09</td></tr><tr><td>lung2</td><td>0.01</td><td>0.12</td><td>0.07</td><td>0.02</td><td>0.18</td><td>0.11</td></tr><tr><td>mac_econ_fwd500</td><td>3.89</td><td>0.65</td><td>0.30</td><td>4.58</td><td>0.89</td><td>0.48</td></tr><tr><td>majorbasis</td><td>4.69</td><td>0.74</td><td>0.26</td><td>15.00</td><td>0.97</td><td>0.42</td></tr><tr><td>mc2depi</td><td>5.70</td><td>1.32</td><td>1.08</td><td>27.13</td><td>1.57</td><td>1.57</td></tr><tr><td>parabolic_fem</td><td>57.21</td><td>1.66</td><td>1.12</td><td>39.10</td><td>2.09</td><td>1.60</td></tr><tr><td>PR02R</td><td>37.99</td><td>2.48</td><td>1.92</td><td>61.01</td><td>3.22</td><td>2.93</td></tr><tr><td>rajat31</td><td>48.39</td><td>11.00</td><td>7.83</td><td>55.91</td><td>26.56</td><td>28.76</td></tr><tr><td>RM07R</td><td>259.28</td><td>9.13</td><td>5.56</td><td>468.04</td><td>12.83</td><td>8.77</td></tr><tr><td>Si41Ge41H72</td><td>65.22</td><td>4.32</td><td>1.72</td><td>6.51</td><td>5.07</td><td>2.34</td></tr><tr><td>ss</td><td>824.91</td><td>10.77</td><td>4.66</td><td>2621.10</td><td>18.19</td><td>10.79</td></tr><tr><td>ss1</td><td>5.74</td><td>0.45</td><td>0.36</td><td>5.33</td><td>0.62</td><td>0.44</td></tr><tr><td>stomach</td><td>1.45</td><td>1.16</td><td>0.24</td><td>11.23</td><td>1.69</td><td>0.49</td></tr><tr><td>t2em</td><td>53.66</td><td>3.52</td><td>2.11</td><td>89.95</td><td>5.57</td><td>4.81</td></tr><tr><td>tmt_unsym</td><td>47.59</td><td>3.28</td><td>4.83</td><td>106.41</td><td>5.57</td><td>8.59</td></tr><tr><td>torso1</td><td>10.62</td><td>1.45</td><td>0.64</td><td>17.61</td><td>1.77</td><td>0.87</td></tr><tr><td>torso2</td><td>0.45</td><td>0.24</td><td>0.19</td><td>1.19</td><td>0.38</td><td>0.32</td></tr><tr><td>turon_m</td><td>6.22</td><td>0.60</td><td>0.71</td><td>3.75</td><td>0.80</td><td>1.02</td></tr><tr><td>vas_stokes_1M</td><td>452.36</td><td>12.00</td><td>4.20</td><td>2427.71</td><td>15.56</td><td>7.42</td></tr><tr><td>vas_stokes_2M</td><td>1419.99</td><td>26.46</td><td>11.93</td><td>6813.20</td><td>34.32</td><td>18.37</td></tr><tr><td>xenon2</td><td>6.66</td><td>1.14</td><td>0.81</td><td>4.01</td><td>1.54</td><td>1.21</td></tr><tr><td>geomean</td><td>12.47</td><td>2.24</td><td>1.23</td><td>20.16</td><td>3.16</td><td>2.07</td></tr></tbody></table>