Solving a Special Type of Optimal Transport Problem by a Modified Hungarian Algorithm
Authors: Yiling Xie
Yiling Luo
Xiaoming Huo
Venue: Accepted by TMLR ()

Abstract: Computing the empirical Wasserstein distance in the Wasserstein-distance-based independence test is an optimal transport (OT) problem with a special structure. This observation inspires us to study a special type of OT problem and propose {\it a modified Hungarian algorithm} to solve it {\it exactly}. For the OT problem involving two marginals with $m$ and $n$ atoms ($m\geq n$), respectively, the computational complexity of the proposed algorithm is $\mathcal{O}(m^2n)$. Computing the empirical Wasserstein distance in the independence test requires solving this special type of OT problem, where $m=n^2$. The associated computational complexity of the proposed algorithm is $\mathcal{O}(n^5)$, while the order of applying the classic Hungarian algorithm is $\mathcal{O}(n^6)$. In addition to the aforementioned special type of OT problem, it is shown that the modified Hungarian algorithm could be adopted to solve a wider range of OT problems. Broader applications of the proposed algorithm are discussed---solving the one-to-many assignment problem and the many-to-many assignment problem. We conduct numerical experiments to validate our theoretical results. The experiment results demonstrate that the proposed modified Hungarian algorithm compares favorably with the Hungarian algorithm and the well-known Sinkhorn algorithm, and the network simplex algorithm.

Published in Transactions on Machine Learning Research (02/2023)
Solving Special Type of Optimal Transport Problem by Modified Hungarian Algorithm
Yiling Xie
ymie3508gotech. edu
School of Industrial and Systems Engineering Georgia Institute of Technology
Yiling Luo
yluo373 edu
School of Industrial and Systems Engineering Georgia Institute of Technology
Xiaoming Huo
huo @wwtecatte...... edu
School of Industrial and Systems Engineering Georgia Instätute of Technology
Reviewed on OpenReview: // net/
forum? k5m8aXTOrG
Abstract
Computing the empirical Wasserstein distance in the Wassersteinn-ddssaann--bbbee indepen- dence test an optimal transport (OT) problem with a special structure. This observation inspires to study a special type of OT problem and propose modified Hungarian algo- rilhm to solve it eraclly. For the OT problem involving two marginals with and n atoms (m n), respectively, the computational complexity of the proposed algorithm is O(m (m) Computing the empirical Wasserstein distance in the independence test requires solving this special type of OT problem, where m = n². The associated computational complexity of the proposed algorithm is O(n)), while the order of applying the classic Hungarian algorithm Onn)). In addition the aforementioned special type of OT problem, it is shown that the modified Hungarian algorithm could be adopted solve wider range of OT problems. Broader applications of the proposed algorithm discussed- solving the one-to-many assignment problem and the many-to-man assignment problem. We conduct numerical ex- periments to validate our theoretical results. The experiment results demonstrate that the proposed modified Hungarian algorithm compares favorably with the Hungarian algorithm, the well-known Sinkhorn algorithm, and the network simplex algorithm.
1 Introduction
One appealing application of optimal transport (OT) and Wasserstein distance Villani 2009 Peyré &
Cuturi 2019) is the independence test. The Wasserstein distance between two distributions 42 on Z is
defined as:
W(n, (2) = inf { d(2, 2'd,,,")) y is a distribution with marginals in and #22
where (Z,d) is metric space (1- (1-Wasserstein distance is considered in this paper). The Wasserstein distance is a metric on probability measures (Villani 2009). To test the independence between the variables Y ~ vi and 2 ~ U2, people utilize the Wasserstein distance between the joint distribution of Y, Z and the product distribution of Y, Z, i... ((T, , 111 @ 122,, where π denotes the joint distribution of Y, Z, and v 822 denotes the product distribution of Y, Z. While the statistical properties of this approach have been intensely investigated (Nies et al. 2021 Mordant & Segers 2022 Wiesel 2022), no existing literature focuses on the computational aspect. In this paper, discuss the following:

---

Published in Transactions on Machine Learning Research (02/2023)
How to compute the empirical Wasserstein distance in the independence test?
In practice, given n i..d. samples (1,,, .., (yn,, zn)} generated from (Y, 2), one can build the statistic W(", 11 @ 72),, where म, denote the corresponding empirical distributions of x and V, respectively to test the independence. Computing W 27 @ 122) is equivalent to solving the following optimization problem: (more details are presented in Section
min E
d((r)) (yk, 2A ifak IIp kj.==1
Σ
Xiju 2
Σ
a
Σ
Xisz
1 Vi,j, 1, 1, n.
n²
k=1
i,=1
where the metric d is usually chosen as d((vi, 2y), (yu. 2)) = 11u. - Ukllp + 1|ag - =2llp. and || llp denotes the lp norm.
1 is an OT problem involving two marginals. One marginal is uniform with n atoms (i.e., we 2 Aisk 1/n, Vk, 1 s & S n), and the other marginal is uniform with n² atoms (i.c., we have En= -==1- 2i = 1/n2, Vi,j, 1 S i S n). Motivated by this structure, we study the following special OT problem:
Problem
have
m n Ciy, u' X'ell' 1=
n1 Σ j=1
m
min
ΣΣ
X
X
X13
my Vi 1,... m:j 1,... , (2)
m
m
where 0 < <n m, my's are positive integers, and Σj=1 my m holds. One marginal of this OT problem is n-dimensional where the probability of each component is prescribed as my/m /m (i.e., we have E X mj/m, Vj, j S n), and the other marginal is uniform with m atoms (i.c., we have Σj= X 1/m, Vi, 1 S $ m). In essence, problem ( is special case of problem (2). where my n, m n Vj 1, 1,... n2 Throughout this paper, consider real-valued entries in the cost matrix and later propose strongly polynomial-time algorithm to solve problem
(2] precisely.
Per Birkhoff's s theorem (Birkhofi 1946). the solution to problem 2 is vertex (whose coordinates are zeros and ones). Then, we could prove the following proposition, and the proof is relegated to the Appendix. 2 is equivalent to the optimization problem
Proposition 1. The optimization problem
u { Xú X1y = 10.1)
m
n C.jo m
m
min X
Σ x,
= 1, Σ X,j = myj Vi I. m; 1, .... n (3)
One may recall the assignment problem, seeing the definition in Section 2 where the permutation matrix is the solution matrix. X E U is similar but different from the permutation matrix: X E 2 is an m x n matrix instead of a square matrix and has mulliple entries of 1 in each column instead of only one entry. In this case, we are not able to directly apply algorithms for the assignment problem, such as the Hungarian algorithm {Kuhn| 1955 Munkres 1957). An approach to obtain the precise solution to problem (3) is first to duplicate the columns of C' and X, then apply the Hungarian algorithm. The computational complexity of this approach is O(m²). Tn this paper, modified Hungarian algorithm is proposed. The algorithm specializes in solving the special type of OT problem 3 which is equivalent to problem 2 with provable lower order (m?n).
For the special type of OT problem 2 we require that one marginal of the OT problem should be uniform. We could further relax the uniform requirement and consider the following more general OT problems:
m7 x
m
min X=Eli
XiC ==1=1
ur
X
0
X
n4
X
mi
Vi 1,.... ,m;j 1,... n M
(4)
M
where 0 < n S m, ny n,'s, my's are positive integers, and Hungarian algorithm could be adapted solve problem and the associated computational complexity is OMP)).
j=1 1 mj
Σ ii n, = M holds. The modified
2

---

Published in Transactions on Machine Learning Research (02/2023)
Back to the Wasserstein-distaree based independence test problem 1 the resulting computational complex- ity of applying the proposed algorithm is O(n²) while the order of applying the classic Hungarian algorithm is (nn)). In this sense, the proposed algorithm is faster. In addition to the application in the Wasser- stein independence test, broader applications of the modified Hungarian algorithm, inchuding solving the one-to-many assignment problem and the many-to-many assignment. problem (Zhu et al. 2011 2016), are investigated. Two practical assignment problems involving the soccer game and agent-task assignment serve as examples to illustrate how to apply the proposed algorithm.
1.1
Related work
1.1.1
Semi-assignment problem
The problem is also called the semi-assignment problem in the literature (Barr et al.)
Kennington & Wang 1992.. Kennington & Wang (1922) proposes strongly polynomial-time to solve the semi-assignment
1977
problem exactly. The proposed modified Hungarian algorithm and the algorithm proposed in Kennington & Wang 1992) are fundamentally different. The algorithm in Kennington & Wang [9922 adjusts the
shortest path augmnenting algorithm (Jonker & Volgenanl, 1987). while our algorithm modilies the Hungarian algorithm Kuhn 1955).. As illustrated in Jonker & Volgenant (987); Kennington & Wang (1992), the primal-dual method based on maximum flow, while the shortest path augmenting algorithm considers the assignment problem minimum cost flow problem and is dual method based on the shortest path. More specifically, the modified Hungarian algorithm is based on modified Kuhn-Munkres theorem and iterates to identify a perfect pseudo-matching in the bipartite graph with some feasible dual variables, while the algorithm in Kennington & Wang (1992) involves constructing the shorlest augmenting path in the auxiliary graph, and the flow pushed along the path (Kennington & Wang 1992). Regarding computational complexity, both algorithms have the order of O(m²n). However, our algorithm is casier understand and implement because there are four phases (columm reduction, reduetion transfer, row reduction augmentation, and shortest path augmentation) in the algorithm in Kennington & Wang [1992 while our algorithm only involves two phases (updating the feasible labeling and improving the pseudo-matthhngg)
Hungarian algorithm
1.1.2 Approximation algorithms
To solve the OT problem, we could apply the exact algorithms or the approximation algorithus. While the proposed modified Hungarian algorithm is an exact OT solver, there are bunch of approximation algorithms. People usually consider the following OT problem:
N
No
No X
min D
C,y
U(a x X
Σ
X1 aj,
B
(5)
where ΣM, B, = and
No 2j=1 aj 1. The approximation algorithms
CCuturi 2013
Dvurechensky et al.
1i
2018 Lin et al. 2019a Xie et al. 2022 in the literature are to obtain an e-approximation X E U((,, B) to such that (X,C) < C² +e, where X" is the solution to (5)
The development. of efficient exact algorithms is meaningful. Notably, precise solutions are needed in some scenarios, and Dong et al. (2020) demonstrates the favorable numerical performanco of the exact solutions over the approximate solutions. Numerical experiments are conducted to compare the modified Hungarian algorithm with the most widely-used approximation algorithm the Sinkhorn algorithm, highlighting the efficiency of our exact algorithm.
1.1.3 Exact algorithms
We review the exact algorithms to solve the OT problem and compare our algorithm with them as follows.
The special type of OT problem 2 is a minimum-cost flow problem and could be solved by the network simplex algorithms. Note that Orlin (1997) proposes the first polynomial-time network simplex algorithm, and Tarjan (1977) further improves the result. The associated computational complexity of applying Tarjan's
3

---

Published in Transactions on Machine Learning Research (02/2023)
algorithm to problem (2) is O(m'n log(m) min minnlog((mmma)).. mn log(m))). where Cmar denotes the maximum
absolute value of the costs if all costs are integers and xo otherwise (Orlin 1997 Tarjan
1997 More
specifically, if the costs are integral, the resulting computational complexity is O log ( log(mO max))
which is comparable to the proposed algorithm; if the costs are not integral, the resulting computational complexity Om"*n² log (m)),, which is worse than our algorithm.
The interior point algorithms can also be customized to solve the minimun-cost flow problems, e.g., Yeh [1889]; Resende & Veiga (1993). . As discussed in Resende & Pardalos (1996), from the perspective of computational complexity, Vavasis & Ye (1994] proposes strongly polynomial-time interior point algorithm for solving the minimum-cost flow problem. The adaption of this method to problem (2) has the order of O 6,5 log(mm)), which is much worse than the proposed algorithm. As demonstrated in Resende &
Pardalos 1996), prior to Vavasis & Ye (1994), the fastest interior point method to solve the minimum- cost. flow problem comes from Vaidya (1989). The computational complexity of the adoption of Vaidya's algorithm to problem (2) is Om log(mGmaxx). The complexity is worse than our algorithm, and the algorithm requires that all costs are inlegers.
As shown in Peyré & Cuturi (019)), could apply other combinatorial algorithms, including the auction algorithm [Bertsekas & Eckstein 1988 p. and the dual ascent algorithm [Bertsimas & Tsitsiklis 1997). to solve the special type of OT problem 2 We compare the proposed modified Hungarian algorithm with them as
follows: Regarding the auction algorithm, the output suboptimal (Peyré & Cuturi 2019 p, while the output of our algorithum is optimal. Furthermore, [Betssias & Tsitsiklis (1997); Bertsekas (1988) illustrate that finite termination of the auction algorithm with an optimal solution could come from the nature of integer-valued valued cost coefficients. In terms of the dual ascent algorithm, [Bertsimas & Tsitsiklis (197)) demonstrates that the
input of integral costs is one of the conditions for finite termination.
In conclusion, finite termination with an optimal solution or fast computation of the aforementioned algo- rithms, including network simplex algorithm, interior point algorithm, auction algorithm, and dual ascent algorithm, require integer-valued costs, while the proposed modified Hungarian algorithm could handle real-valued costs. Although most of the problems in practice have rational costs, which can be scaled to integer- valued costs, is more convenient use the proposed algorithm since it allows direct input of any real-valued costs. Also, looking into OT problems with real-valued costs itself is of much theoretical interest.
Similar discussions and comparisons could be developed for the more general OT problem could be reformulated to a minimum-cost flow problem and can be solved by the exact solvers mentioned above. Among the algorithms, the computational complexity of the network simplex algorithm can be comparable to the proposed modified Hungarian algorithm O(M2)) under the assumption that all costs are integral.
Problem (4)
1.1.4 Independence criteria
There some other independence criteria based OT the Wasserstein distance. Shi et al. (2020) Deb & Sen (2021) design the independence criterion by combining the distance covariance and OT. Liu et al. (2022) applies the entropy-regularizee OT. Wiesel (2022) utilizes the nested Wasserstein distance. The criterion proposed in Mordant & Segers (2022) is based on the 2- 2-Wasserstein distance under the quasi-
Gaussian assumption. This paper considers the criterion based on the general Wasserstein distance as shown in the formulation
1.2 Our contributions:
We propose modified Hungarian algorithm to solve special type of OT problem 2 The modification enables us to deal with the scenario where two marginals have different sizes of atoms, and the aloms in one of the marginals have multiple assignments. Further, the proposed modified Hungarian algorithm could be extended to solve more general OT problems. Moreover, the applications of the proposed algorithm are explored: adopting the modified Hungarian algorithm to solve the Wasserstein independence test problem 1 the one-to-many assignment problem, and the many-to-many assignment problem. Finally, several
4

---

Published in Transactions on Machine Learning Research (02/2023)
numerical experiments are carried out using Python to show the favorability of our algorithm over the classic Hungarian algorithm, the Sinkhorn algorithm, and the network simplex algorithm.
1.3 Organization
The remainder of this paper is organized follows. In Section 2 we introduce some basics of graph theory. In Section g
we propose the modified Hungarian algorithm and compute its computational complexity. In independence test problem. In Section 5 we apply the modified Hungarian algorithm to the one-to assignment problem and the many-to-many assignment problem. In Section G we carry out various numorical experiments on both synthetic data and real data to validate our theoretical results and show the lavorability of our algorithm. We discuss some future work in Section
Section 4 apply the modified Hungarian algorithm to the Wasserstein-ddssaacee-lll
2 Preliminaries
Some definitions related to combinatorial optimization and graph theory (Ahuja et al. 1988 Suri 2006 Burkard et al. 2012) are introduced. They will be needed in the rest of this paper.
Definition (Assignment problem). Given an & x l cost matrix with components Cij 2 0, E [k] =
(1, ,), the assignment problem is to solve minş Σ Cip(i), where & is the permutation of set [k].
Definition 2 (Bipartite graph). A graph G = EV, is called a bipartite graph if its nodes can be portitioned into two subsets V1 and V2 so that for each edge (V1, 02) in E, 01 V1 and 02 E V2.
Definition 3 (Matching and perfect matching).
matching in the bipartite graph G (V, a subset M C E such that at most one edge in M is incident upon v, Vu E V. M is called a perfect matching if every node in G coincides with ezactly edge of M.
Desinition 4 (Weighted bipartite graph). A weighted bipartile graph is a bipartite graph where each edge has weight w(.) 2 0. The weight of a matching M is the sum of the weights of edges in M, i.e., EeEM w(e). Definition 5 (Tabeling and feasible labeling). For weighted bipartite graph G (V, F.,, where V V1UV, labeling is function 1 V -> RR A feasible labeling is one labeling such that [(01] + I(02) 2 w(v1, U2), Ve E V1, E V/2
Definition 6 (Equality graph and neighbor). The equality graph w.r.t. the labeling 1 is G' (V, Ex) where Ei {(11, 02) ((v1) + I(02) w(11, v2)}. The neighbor of U2 V2 and S C V2 is defined (1(22) = {u: (21,02) E E and (S) = UsEss (1)2), respectively.
Definition 7 (Alternating and augmenting path). Let M be matching of the bipartite graph G = E). A wath in G = (V, E) is sequence of distinct nodes and edges in, (iz, 12), ... 1,)-,,,,,,,, sulisfying (ix, i++1) E E for each k = 1, 1. A path in G = EV is alternating if its edges alternate between M and E - M. An alternating path is augmenting if both endpoints do not coincide with any edges in M.
3 Modified Hungarian algorithm
In this section, we propose modified Hungarian algorithm to solve the special type of OT problem which is equivalent to problem
2
3.1 Review of the Hungarian algorithm
We first review the Hungarian algorithm (Kuhn 1955 Munkres 1957). Recall that the assignment problem is to solve ming Em Coo())- If we negate the costs and add the maximum of the costs to each component, solving the assignment problem is equivalent to finding maximum weighted matching in the weighted bipartite graph with weights jii.j) = max, Ctj Cy. The Kuhn-Munkres theorem (Munkres) 1957) shows that finding a maximum weighted matching is equivalent to finding a perfect matching on the equality graph associated with some feasible labeling in the bipartite graph. In this regard, the Hungarian algorithm solves the assignment problem by identifying perfect matching on some equality graph in the weighted bipartite
5

---

Published in Transactions on Machine Learning Research (02/2023)
Figure 1: pseudo-malching (left) and perfect pseudo-malching (right), where = 3, m = 9, 9,m = 2,12 =
3, mg
graph. One first generates an initial feasible labeling and an associated equality graph; then proceeds to look for an augmenting path to augment the matching. (if an augmenting path does not exist, update the feasible labeling.) If the matching is perfect the equality graph concerning some feasible labeling, stop and output the matching.
3.2 Freudo-matching
In problem 3 X E U has one entry of 1 in each row, multiple entries of 1 in each column, and 0's elsewhere, Since a permutation matrix corresponds to a (perfect) matching in the bipartite graph, we define 'pseudo- matching' in the bipartite graph G = (Vi U V2, E) to describe X. Vi has m nodes representing the rows of X while V/2 has n nodes representing the columns of X. Notice that we usually have m > n. In this case, each node in Vi coincides with at most one edge, while multiple edges are allowed to connect with nodes in V2. See the formal definition in Definition
Definition 8 (pseudo-matchingg perfect n. PM C E is a pseudo-matchine if every node of V1 coincides with at most one edge of PM, and jth node of V22 coincides with al most my edges of PM, where Σド= my m. Purthermore, if every node of Vi coincides with eractly one edge of PM and jth node of V/2 coincides with exactly m, edges of PM, PM is called a perfect pseudo-matchingg
In the bipartite graph G, where 1 [V2l
Figure 1 is example of (perfect) pseudo-matching, where 3, 9, mi 2, m2 m3 4. Under this setting, each node in the left-hand side of the graph coincide with at most one edge, while each node in right-hand side can coincide with at most 2 edges, 3 edges, and edges, respectively.
3.3 Our algorithm
Solving problem (3) is equivalent to looking for a maximum weighted pseudo-matching in the bipartite graph. We develop a modified Kuhn-Munkres theorem based on the psoudo-matching. See Theorem 1 (The proof can be found in the Appendix.) demonstrates that we only need to find a perfect maaudd-matching on some equality graph to solve problem (3)
Theorem 1 (Modified Kuhn- Munkres theorem). If l is feasible labeling on the weighted bipartite graph G = (V, E), and PM C E₁ is a perfect pseudo-matching on the corresponding equality graph G' (V, PM is marimum weighted pseudo-matching matching.
Equipped with the modified Kuhn-Munkres thcorem, we design modified Hungarian algorithm (Algorithm D. The definitions used in the algorithm are specified in Definition g The modified Hungarian algorithm improves either the feasible labeling (adding edges to the associated equality graph) the pseudo-matching until the pseudo-matching is perfect some equality graph w.r.t. some feasible labeling, The algorithm improves the pseudo-matching by generating senndo-augmenting paths and then exchanging the edge status along the paths. This process is called the pseudo-aaugmenting process. Also, we force the pseudo-augmenting paths emanating from V2, which have a lower order of nodes.
Definition 9 (Free, matched, pseudo-matched psendo-alternating path, pscudo-augmenting path). Let PM
be pseudo-malching of G (V, F)).
6

---

Published in Transactions on Machine Learning Research (02/2023)
Figure 2: psendo-augnentingg process, where n 3, m 9, m. 2, m2 = 3, m.
If the node is in V1, it is pseudo-matched if it is endpoint of some edge in PM; if the node is the jth node in V2, il is pseudo-matchell if il is endpoint of m, edges in PM. Otherwise, the node is free.
If the node VEV, we say it is matched if it is on endpoint of some edge in PM.
A path is pseudo-ale alternating if its edge alternates between PM and E PP. A pseudo alternotitnnn path pseudo-augmenting if both endpoints are free.
An example of the pseudo-augmenting process given in Figure 2 The solid line means that the edge belongs to the pscudo-matching. The dashed line means that the edge belongs to the equality graph but does not belong to the pseudo-matching, Node B and node C pseudo-matched. Edge A-B and edge C-D are not in the pseudo-matchine matching. In this sense, A-B-C-D is pseudo-alternating path. Because node A and node D are free, A-B-C-D is a pscudo-augmenting status of the edges: delete B-C from the pseudo-matching and enter A-B, C-D into the pseudo-matchingg The pseudo-matching has been improved in this way.
The pseudo-augmenting process is to exchange the
3.4 Computational complexity
We now analyze the computational complexity of Algorithm [ Similar to the Hungarian algorithm (Suri 2006), we keep track slack,, = min,)es(((vr)) l(v2) w(1,, 02)}, Ye & T. The computational cost increases when computing o via slacks, updating the values of slacks, and calculating the labeling.
The number of edges of the pseudo-matching increases by after one loop, O(m) loops needed to form perfect pssuudomaaching. There two subroutines in each loop: the first is to update the feasible labeling (Step ) and the second is to improve the pscudo-matching (Step 3 In the procedure of updating the feasible labeling, since there are nodes in V/2 the improvement occurs Onn) times to build pseudo- alternating tree. Tn each time, computing on, updating the slacks, and caleulating the labeling cost ((m). the procedure of improving the pseudo-matee matching, when a new node has been added to S, it costs O(m) update slacks, and (n) nodes could be added. On the other hand, when node has been added to T, we just remove the corresponding slacke We conclude that each loop costs O(mn), so the total computational complexity of Algorithm 1 to solve problem (3) is Om22)). We summarize the analysis above in Theorem 2
Theorem 2. The compulational complexity of applying the modified Hungarian algoritlem to solve problem is O(m²n).
3.5 Solving more general OT problem
We could adapt the proposed modified Hungarian algorithm to solve the more general OT problem 4 We first rewrite problem { as the formulation of the special type of OT problem (2) by duplicating the rows of the cost matrix, seeing Proposition (the proof is relegated to the Appendix).
Proposition 2. Problem
4 is equivalent to the following optimization problem:
M
M
min xteut
cl. ut =
x1 2
Σ
x
M.
Vi = 1. m:j 1, .....,, (6)
M
M
1=
111

---

Published in Transactions on Machine Learning Research (02/2023)
Algorithm 1: Modified Hungarian Algorithm
Generate an initial feasible labeling I: Vuz E V2,((02) = V, E Vi,((v1) = max (10 U2)} and
initialize a pscudo-matching M in E;;
if M is perfect pseudo-malching then
Stop
else
Pick up free node Utree E V22 Set S [Utree]: 0;
for V1 E Vi is matched to Ufree do
= TUV1
end
end
if Ni(S) - = 0 then update labeling such that forcing Ni(S) T * 0:
((()) ay VES
a , sT. + [(02) w(n 02)), 1(v)
((v) + ar, ve T
((u),
otherwise
end
if (S(S) T * 0 then
pick U1 ((()) -;; if is free then
VLeee + 11 psendo-angmenting path. Pseudo-augment the matchingg M. Go to Step end
if 01 is matudo-matched to z then extend the pseudo-matching tree: = S (2) (ni}:
for V1 V1 is matched to z do
T = TU
end
end
Go to Step 2 end
dependent p
dependent
independent p
independ
(nt#operation)
In(sample size)
In(sample size)
Intsample size)
Intsample size)
mocMecd langedarr medt Harggarannee
modites langarion
Hurgatinn
Figure 3: Comparison with the Hungarian algorithm on synthetic data
where Cf is an M matriz generated by duplicating the ith row of C n, times:
Cis Ci ni
1515 n +1sts 2n,25
ct C1
Problem 6 belongs to the special type of OT problem 2 We could apply the proposed modified Hungarian algorithm problem (6) and get the exact solution to problem 4 The resulting computational complexity
is O(M2)..
8

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p
Independent case, p=
Independent case, =
Int#operation)
in(sample size)
In(sample size)
intsample size)
In(sample size)
modfice Nangetan vestage
modi od Nurgunan be...
modiieed Pangaatin ward
Plangeraan sherage
Murgat on Des
Her Hergunannnn
Figure 4: Comparison with the Hungarian algorithm on CIFAR10
dependent p
dependent case, p
independent p
independent
In(#operation))
10
In(sample size)
In(sapple size)
Intsample size)
In(sample size)
modi
Figure 5: Comparison with the Hungarian algorithm on Wisconsin breast cancer data
3.6 Compare with the Hungarian algorithm
We discuss the modification and the improvement of the proposed modified Hungarian algorithm
Our algorithm modifies the Hungarian algorithm to solve wider class of problems. Hungarian algorithm specializes in the assignment problem which corresponds to the matching in the bipartite graph. The matching is embedded with the 'one-to-one' structure. However, as shown the solution structure of problem (3), the problems considered in this paper have 'one-to-many' structure. To deal with the 'one- to-many' structure, pseudo-matching is defined. The proposed algorithm is developed to identify the perfect pseudo-matching while the Hungarian algorithm identifies the perfect traditional matching. In addition, instead of building the augmenting paths in the Hungarian algorithm, we establish the pseudo-augmentiigg path in the proposed algorithm as shown in Figure 2
Our algorithm has a lower order of computational complexity. The complexity converting the problems to the assignment problem and solving them by the Hungarian algorithm only depends on the larger size of the marginals- Om=)) To reduce the complexity, we force the beeudo-augmenting paths emanating from V22 which corresponds to the marginal with the smaller size. In this way, the complexity of the direct application of the proposed algorithm depends the sizes of both marginals O(m²n). Hence, the proposed modified Hungarian algorithm will outperform, especially when m > n.
4 Application to the independence test using the Wasserstein distance
In this section, apply the modified Hungarian algorithm to the Wassestein-distaane--aas independence test, which originally motivates to study the special type of OT problem 2
Suppose that there are n i.i.d. samples {311, ... (Uny 2n }, where (J), (Y,2),Y ~ v1, 2 ~ V2.
Recall æ denotes the joint distribution of Y, Z, and one could prove the following equivalence:
Y12
W(T, vn @ v2) = 0.
which follows from the fact that the Wasserstein distance valid metric between probability measures. Given the empirical data, we utilize the statistic W 2n @ 22) to test the independence between Y and Z,
9

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, =
dependent case, p p=
Independent case, p=
Independent case, p
Int#operation)
in(sample size)
In(sample size)
intsample size)
In(sample size)
modfice Rangetan vestage
modl Murguman Ena.
mede Pangartnn ward
Mangerian overage
Mangat on Des
Her Nurgenannnos
Figure 6: Comparison with the Hungarian algorithm on DOT-benchmark with 512 x 512 resolution
dependent case, p
dependent case,
independent case,
independent case, p=
in(time)
in(sample size)
Insample size)
In(sample size)
In (sample
nocffie Nangetaa everage
modi ooddled Hurgaran
modified Pangartin wose
Plangarian average
Hangar on best
Har garian woret
Figure 7: Comparison with the Hungarian algorithm synthetic data
where 7,0 G denote the empirical distributions and have the following expressions:
n
n²
11
02
0
Plug in the Wasserstein distance formula, the resulting optimization problem is:
min Xell
Σ
d((nr, 2i), (ykr 21) 21)) Xigitt,
ij..k1=1
where
1 k =1 x
II Xijua
Σ k.Tml
X
ijckl
X iiskl
vi.jk,,==, n.
n
1,j=1
It is worth noting that XijuM = 0, + 1. we let X = Ei= XijM, problem can be simplified as problem [ Problem (1] belongs to the special type of OT problem, where my n, Vj. S S n Adopting the Hungarian algorithm to problem (1) costs O(n²), while adopting the proposed Hungarian algorithm directly costs O(n²).
5 Application to the one-to-many assignment problem and the many-to-many assignment problem
In this section, we proceed explain how to apply the modified Hungarian algorithm to solve the one- to-many assignment problem (corresponding to problem ( and the many-to-many assignment problem (corresponding to problem
[4] The applications are shown based on two practical examples.
Example 1 (one-to-many assignment problem): An assignment problem involving the soccer ball game mentioned by Zhu et al. (2011) is considered here. Suppose a coach is tasked to choose players from a soccer team with m, players (a,, a.,). There are d roles (ri,-.. .ra). is assumed that m, > d. Suppose each player's performance evaluation of each role known. The overall performance evaluation of the team is the of each selected player's performance evaluation of its assigned role. The optimal strategy is to
10

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p=2
independent case, p
Independent case, p
Initimel
in(sample size)
in(sample size)
In(sample size)
Isamplee size)
modfice Nangetan vestage
-.. modi Hurgunan
modllee Pangatin ward
Plangerkan sherage
Murgar on Des
Her puraannnsss
Figure 8: Comparison with the Hungarian algorithm CIFAR1O
dependent case, p
dependent case, p
independent case, p
independent case, p
In(time)
In(sample size)
20
53
35
In(sapple size)
25 In(sample size)
In(ampeee size
Figure 9: Comparison with the Hungarian algorithm on Wisconsin breast cancer data
maximize the overall performance evaluation of the team. The coach should solve the following optimization problem:
m
m Σ i=1
max AEA
Pip A= Aij = 10,1)
Σ /=1
Aij 11
= 1, mmij = 1,...
where P.j 2 0 denotes player a,'s performance evaluation of role j, Asj = 1 means player aj is selected as role j while Aij 0 means the player is not selected as role j.
If mi = Σj 1 Tj, the optimization problem above belongs to the special type of OT problem (3 where n = d, m - プj= The modified Hungarian algorithm could be applicd to find the optimal strategy, and the resulting computational complexity is O(dn)).
mm Σラ= rj can't apply the modified Hungarian algorithm directly. To make the problem tractable, create one more role, and each player's performance evaluation of this role is 0. Players who are not selected are 'assigned' to this role by default. In this scenario, our goal is to solve the following optimization problem:
d+1 ΣΣ At= ATEA! ==1==1
A - (0.
d+1 Σ 1=1
1 =
1 E1
15j 1=1 mi En j=d+1
max
where we append P by adding one more column of zeros to get p It belongs to the special type of OT problem, where π d 1,m m.. Then, we could apply the modified Hungarian algorithm to solve the problem, and the resulting computational complexity is O(dd + 1)m3).
Note that the computational order of applying the algorithm developed by Zhu et al. 2011 is O(m)), which is worse than the proposed modified Hungarian algorithm.
The following example an agent-task assignment problem mentioned by Zhu et al. (2016). Assume there are my tasks (f,, (mm)) and n, agents (a,, ... an in total, It is assumed that ny < m2. Each task should be undertaken by many agents, and each agent can perform many tasks. To be more specific, task 4 must be assigned to agents, agent aj can perform at most
Example 2 (many-to-man tany assignment problem):
11

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p=2
Independent case, p
Independent case, p
Initimel
in(sample size)
in(sample size)
In(sample size)
Isampple size)
modfice Nangetan vestage
modi Hurguran Ens.
madh Pangarin ward
Pangerian overage
Mangar on Des
Her Nergenanentl
Figure 10: Comparison with the Hungarian algorithm on DOT-benchmark with 512 x 512 resolution
dependent case, p
dependent case, p
independent case, p
independent case, p
inttime)
22 inaample size)
35 In(sample size
In(sample size)
In(sample size)
wangarorr awarage
Hunjarion
3e0rage
Figure 11: Comparison with the Sinkhorn algorithm on synthetic data
8y tasks. Suppose the performance evaluation of each agent performing each task is known. The optimal assignment plan to maximize the overall performance. The resulting optimization problem is as follows:
mh AEASPh A = A' - (0.1) [A,
max
Σ 1=1
A 4,
A =
Σ
Aly S 8jz Vi - 1,,.. m2;j = 1,
1=
where Pj 2 0 denotes agent 's performance evaluation on task 4,, A,j means that agent aj is assigned to perform task t; while Ale- = means that agent aj is not assigned to perform task t,.
If Σ 4 = Ejn Sy, the optimization problem follows the formulation of the problem where n ni, M - 8j. We could apply the modified Hungarian algorithm to find the optimal assignment plan, and the resulting computational complexity (1)
th ビグ=1
If TE < - Sj. we create one more task which must be performed by - 4) agents, and each agent's 's performance of this new task equals This reformulation promises that each agent performs the maximum amount of tasks. Accordingly, we need solve the following optimization problem:
ma+1 ΣΣ 4 A - {AA,
10. -
4 1
1 1 S m2
ma+1
Σ
Σ A
max AFEAA
Sj.
89 Em
j= m m2
+1
1=1
where we append P' by adding one more row of zeros to get Pt. It follows the formulation of problem 4 where n nı, M = و sj. We could adopt the method introduced in Section 3.5 and the resulting computation complexity is O(n, s)),,)))
Note that the computational order of applying the algorithm developed by Zhu al. (016) which has higher computational burden than our proposed method.
O(E))
"1 8))3),
6 Numerical experiments
In this section, we carry out experiments on the Wasserstein independence test problem on a synthetic dataset, CIFARIC Krizhevsky et al. 2009), Wisconsin breast cancer dataser Dua & Graff 2017 and
12

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p=2
independent case, p
Independent case, p=
Initimel
In(sample size)
In(sapple size)
In(sample size)
Inlmaplee size)
modefed Nangeria unarage
modh ool Marguran
modlfee hungarin WO73
Sertron aearage
Figure 12: Comparison with the Sinkhorn algorithm CIFAR10
dependent case,
dependent case, p
independent case,
independent case, p
Inttime)
in(sample size)
inaample size)
In(sample size)
In(sample
mediiiee Hangeror averaee
modi Hurganan
madifes WOIS
Figure 13: Comparison with the Sinkhorn algorithm on Wisconsin breast cancer data.
DOT-benmark (Schrieber et al 2016). We compare the proposed modified Hungarian algorithm with the Hungarian algorithm, the Sinkhorn algorithm, and the network simplex algorithm. The numerical results show the favorability of applying the proposed algorithm over the Hungarian algorithm, the Sinkhorn algorithm, and the network simplex algorithm.
6.1 Experiment setting
Our algorithm is adaptive any metric. The foregoing experiments are based on the Ip norm- base-based metric: d((x;, vy). (Ik, 3²) ||= - =x lyy - 3 lpl More specifically, we examine how the modified Hungarian algorithm, the Hungarian algorithm and the Sinkhorn algorithm perform when create one dependent case and one independent case with different sample sizes for each dataset and run the algorithms on each case 10 times. We plot the worst, best and average number of numerical operations and/or running time for each case.
1 and p 2. We
Synthetic data: Suppose that there are independent variables X ~ 30, where 1.0 10- dimensional vector with all ones and To is the identity matrix; and Y = -Yi,,, Y25)",, where Y'' are We calculate the empirical Wasserstein distance in (1) independent
independent and follow Unif(10,,0;; case: between X and Y; (2) dependent case: between X and Z (where Z X1 Y1, X1 is the first 5 coordinates of X, Yi is the first 5 coordinates of Y).
Breast cancer data: There are 569 instances, and each instance possesses 30 features. Each instance is a 30-dimensiinal vector, and we rescale the components to [0, 1]. There are two classes of instances: benign and malignant. Let E R330 be the distribution generated uniformly from the benign class, and R30 be the distribution generated uniformly from the malignant class. We calculate empirical Wasserstein distance in (1) independent case: between X1 and Y2 (where X, is the first 5 coordinates of X, Y2 the last 25 coordinates of Y); (2) dependent case: between X and Z (where Z X1 Y1, X1 is the first 5 coordinates of X, Y is the first 5 coordinates of Y, mcans the coordinato-wise product).
CIFAR10: Each image in CIFARIO contains 32 x 32 pixels, and each pixel is composed of 3 color channels. Each image is essentially a 3072-dimensional vector. Then, we rescale the vector components to [o, 1]. Suppose X E 3072 is the distribution generated uniformly from the images of classcs: airplane, automobile, bird, cat, and deer; Y 3077 the distribution generated uniformly from the images of other five classes. We calculate the empirical Wasserstein distance in (1) independent case: between X and Vi (where Y1 is the
13

---

Published in Transactions on Machine Learning Research (02/2023)
de depddennt case. p=
dependent case, p
independent case.
Independent case, p
16.
Initimel
15.
in(sample size)
In(sample size)
In(sample size)
In(sample size)
modeied Nanpertar uearage
modh cnd Marguran
modlfed Pungatin Wasd
S4iton werage
Figure 14: Comparison with the Sinkhorn algorithm DOT-benchmark with 512 x 512 resolution
dependent case, p=
dependent case, p=2
independent case, p
Independent case, p
(n(cime)
20 in(sample size)
30
20 Inaamplee size)
24
Islample size)
In(sample size)
Figure 15: Comparison with the network simplex algorithm on synthetic data
first 1536 coordinates of Y); (2) dependent case: between X and Z (where Z = X2/2 11/2, X2 is the last 1536 coordinates of X, Y1 the first 1536 coordinates of Y).
DOT-benchmark contains images with different resolutions from 10 classes. Each image is essentially r rrx-- dimensional vector, where 32, 64, 128, 256, 512. Then, we rescale the vector components [0, 1]. Suppose X R is the distribution generated uniformly from the images of classes: GRFrough, RFmoderate, Cauchy Density, MicroscopyInages, Shapes; Y E gexr is the distribution generated uniformly from the images of other five classes. We calculate the empirical Wasserstein distance in (1) independent case: between X and Y (2) dependent case: between X and Z (where Z = X/2 + Y/2). To save space, we relegate the numerical results when r - 32, 64, 128, 256 to the Appendix.
6.2 Comparison with the Hungarian algorithm
We compare the modified Hungarian algorithm with the classic Hungarian algorithm. The results terms of numerical operations are presented Figure 3
The figures illustrate that the proposed algorithm gains factor in computational complexity when solving the proposed special type of OT problem. To be more specific, notice that the slope of In number of numerical operations) over In (sampll size) indicates the order of the associated algorithm, and the slope of our algorithm is around while the slope of the Hungarian algorithm around 6. This observation implies that the order of applying our algorithm is O(n)) while the order of applying the Hungarian algorithm is Onn)). Such observations are consistent with our theoretical results.
The results in terms of rumning time are presented in Figure 7 instances, especially for larger sample size the modified Hungarian algorithm is faster than the Hungarian algorithm. It indicates the practical improvement of the modified Hungarian algorithm over the classic Hungarian algorithm.
9 10
One may observe that for almost all
6.3 Comparison with the Sinkhorn algorithm
We compare the modified Hungarian algorithm with the Sinkhorn algorithm. Among the state-of the-art approximation solvers for OT problems, the first-order approximation algorithms (Dvurcchensky al. 2018 Lin et al., 2019bb Guo et al.) (2000) are mainly employed solve the balanced case (m = n). The Sinkhorn algorithm could deal with the unbalanced scenario (m # n) and is widely used in all kinds of OT-related
14

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case,
dependent case, p=2
Independent case, p
Independent case,
Initimel
in(sample size)
in(sample size)
In(sample size)
Isammplee size)
incal Nurgaran ererage
med hal hungar test
mediled Mangarian ward
Mitmo Sorglee uva ags
Rutwork arrpkko bes. Assmett
Figure 16: Comparison with the network simplex algorithm CIFARIO
dependent case, p=
dependent case, p=
independent case, p=
Independent case, p
in(time)
20 Insample size)
50
30
20 In(sample size)
25 20 24
Inlamplle size)
Intamplee size)
Figure 17: Comparison with the network simplex algorithm on Wisconsin breast cancer data
models. Therefore, we choose the Sinkhorn algorithm as the baseline and then investigate the performance of independence test problem. When we implement
the Sinkhorn algorithm in the Wasserstein-ddittance--ee the Sinkhorn algorithm, we set the regularization parameter as 0.1 and the accuracy as 0.0001.
The results in terms of running time are presented in Figure terms of numerical operations to the Appendix. For most of the scenarios, the average rumning time of the modified Hungarian algorithm is less than the Sinkhorn algorithm. Moreover, the performance of the modified Hungarian algorithm has lower variance than the Sinkhorn algorithm. The results demonstrate that our proposed modified Hungarian algorithm should be chosen if one is interested in obtaining solutions with a high accuracy.
We relegate the results in
6.4 Comparison with the network simplex algorithm
We compare the modified Hungarian algorithm with the network simplex algorithm. The results are pre- sented in Figure 15 the package requires integer-valued input, round the costs to the nearest integers from the below. Ac- cording to the experimental results, could conclude that the modified Hungarian algorithm is superior to the network simplex algorithm.
18 We run the network simplex from networkX library in Python. Considering
7 Discussion
A modified Hungarian algorithm is developed to efficiently solve a wide range of OT problems. Theoretical analysis and numerical experiments demonstrate that the proposed algorithm compares favorably with the Hungarian algorithm and the Sinkhorn algorithm. In addition the computational aspects, broad applica- tions explored, including the -inessrsteii-distnnceeeee based independence test, the one-to-mant assignment problem and the many-to- many assignment problem. The many- many-to-me taanyyttyynyyyyymm assignment problem closely relates to practical problems involving service assignment problems (Ng et al. 2008), sensor networks (Bhardwaj) & Chandrakasan 2002), and access control (Ahn & Hu 2007). Future work along this line is to apply the proposed algorithm to problems involving engineering and control. Also, there some possibility of applying the proposed algorithm to some unsupervised learning problems. For example, the clustering problem could formulated as an OT problem (Genevay al.) 2019). Assume that there are n clusters and samples in total, and each cluster has my samples. If we want to identify the cluster assignment. to minimize the
15

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case,
dependent case, p=2
Independent case, p
Independent case, p=
Initimel
in(sample size)
in(sample size)
In(sample size)
In(smplle size)
inca Nargaran meraage
red hal Fungar teal
modiled Mangarian Ward
Milmo pargles ava
Retwork arrpcoo Des
.... Ausnett strp
Figure 18: Comparison with the network simplex algorithm DOT-benmark with 512 512 resolution
'distance' between cluster 'centers' and the associated assigned samples, we are solving the special type of OT problem this paper. The future work along this line may be to find a scheme determine 'distance' and 'centers' to promise desirable model performances.
Acknowledgement
The authors would like thank the Action Editor and anonymous reviewers for their detailed and construc- tive comments, which enhanced the quality and presentation of the manuscript.
This project is partially supported by the Transdisciplinary Research Institute for Advancing Data Science (TRIAD), https: grrtcch.atcch...mm edu adu/data/triad which is a part of the TRIPODS program at NSF and locates at Georgia Tech, enabled by the NSF grant CCF-1740776. The authors are also partially sponsored by NSF grants 2015363.
References
Gail-Joon Ahn and Hongxin Hu. Towards realizing a formal RBAC model real systems. In Proceedings of the 12th ACM symposium on Access control models and technologies, pp. 215-224, 2007.
Ravindra K Ahuja, Thomas L. Magnanti, and James B Orlin. Network flows. 1988.
Richard S Barr, Fred Glover, and Darwin Klingman. A new allernating basis algorithm for semi-ossigrmmnt networks. Business Rescarch Division, Graduate School of Business Administration ...., 1977.
Dimitri P Bertsekas. The auction algorithm: A distributed relaxation method for the assignment, problem. Annals of operations research, 14(1):105 123, 1988.
Dimitri P Bertsekas and Jonathan Eckstein. Dual coordinate step methods for linear network flow problems. Mathematical Programming, 42(1):203-245 203223333332222 1988.
Dimitris Bertsimas and John Tsitsiklis. Introduction to linear optimization, volume 6. Athena Scientific Belmont, MA, 1997.
Manish Bhardwaj and Anantha P Chandrakasan. Bounding the lifetime of sensor networks via optimal role assignments. In Proceedings. Fuenly-First Arnual Joint Conference of the IEEE Computer and Communications Societies, volume 3, pp. 1587-1596. IEEE, 2002.
Garrelt BirkhofT. Three observations on linear algebra. Univ. Nuc. Tacuman, Rev. Ser. A, 5:147 151, 1946.
Rainer Burkard, Mauro Dell' Amico, and Silvano Martello. Assignment problems: revised reprint. SIAM, 2012.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural informalion processing systems, 26, 2013.
Nabarun Deb and Bodhisattva Sen. Multivariate
distribution-fre nonparametric testing using measure transportation. Journal of the American Statistical Association, pp. 1-16, 2021.
16

---

Published in Transactions on Machine Learning Research (02/2023)
Yihe Dong, Yu Gao, Richard Peng, Ilya Razenshteyn, and Saurabh Sawlani. A study of performance of optimal transport. Xiviv preprint arXim:2005.011 02.0011111 182, 2020.
Dheeru Dua and Casey Graff. UCT machine learning repository, 2017. URL http::/arccive..iess wwwpeeeeiis....eeeee.e.ee.ee.eee edu/
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Com- plexity by accelerated gradient descent is better than by Sinkhorn's algorithm. In International conference on machine learning, pp. 1367-1376. PMTR, 2018.
Aude Genevay, Gabriel Dulac-Aroodd, and Jean-Philipps Vert. Differentiable deep clustering with cluster size constraints. arXiv preprint arXiv:1910. 09036, 2019.
Wenshuo Guo, Nhat Ho, and Michael Jordan. Fast algorithms for computational optimal transport and Wasserstein barycenter. In International Conference on Artificial Intelligence and Statistics, pp. 2088- 2097. PMLR, 2020.
Roy Jonker and Anton Volgenant. A shortest augmenting path algorithm for dense and sparse linear assign- ment problems. Computing, 38(4):325-344 340, 1987.
Kennington and Zhiming Wang. A shortest augmenting path algorithm for the semi-assignment problem. Operations Research, 40(1):178 187, 1992.
Alex Krishevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
Harold W Kuhn. The hungarian method for the assignment problem. Naval research logistics quarterly, 2 (1-2):83-91 97, 1955.
Tianyi Lin, Nhat Ho, and Michael Jordan. On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms. In International Conference on Machine Learning, pp. 3982-3991. PMLR, 2019a.
Tianyi Lin, Nhat Ho, and Michael I Jordan. On the efficiency of Sinkhorn and Greenkhorn and their acceleration for optimal transport. arXiv preprint arXiv:1906.01449 2019b.
Lang Liu, Soumik Pal, and Zaid Harchaoui. Entropy regularized optimal transport independence criterion. In International Conference on Artificial Intelligence and Statistics, pp. 11247-11279. PMLR, 2022.
Gilles Mordant and Johan Segers. Mcasuring dependence between random vectors via optimal transport.
Journal of Multivariale Analysis, 189:104912, 2022.
James Munkres. Algorithms for the assignment and transportation problems. Journal of the sociely for industrial and applied mathematics, 5(1):32-33 38, 1957.
Vincent TY Ng, Boris Chan, Louis LY Shun, and Ringo Tsang, Quality service assignments for web services. In 2008 IEEE International Conference on Systems, Man and Cybernetics, pp. 2219-2224. IEEE, 2008.
Thomas Giacomo Nies, Thomas Staudt, and Axel Munk. Transport dependency: Optimal transport based dependency measures. arXio preprint arXiv:2105. 05.02073, 2021.
James B Orlin. A polynomial time primal network simplex algorithm for minimum cost flows. Mathematical Programminng 78(2):109 129, 1997.
Gabriel Peyré and Marco Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 11(5-6):355-600 607, 2019.
Mauricio GC Resende and Panos M Pardalos. Interior point algorithms for network flow problems. Advances in linear and inleger programming, pp. 147-187, 1996.
17

---

Published in Transactions on Machine Learning Research (02/2023)
Mauricio GC Resende and Geraldo Veiga. An implementation of the dual affine scaling algorithm for minimum-cost flow on bipartite uncapacitated networks. STAM Journal on Oplimizalion, 3(3):516-55) 537,, 1993.
Jörn Schrieber, Dominic Schuhmacher, and Carsten Gottschlich. Dotmark-a benchmark for discrete optimal transport. IEEE Access, 5:271-288 2016.
Hongjian Shi, Mathias Drton, and Fang Han. Distribution-fres consistent independence tests via center- outward ranks and signs. Journal of the American Statistical Association, pp. 1-16, 2020.
Subhash Suri. Bipartite matching & the Hungarian method. Notes, Department of Computer Science, University of California, Santa Barbara, 8, 2006.
Robert E. Tarjan. Dynamic trees search trees via euler tours, applied to the network simplex algorithm. Mathematical Programming, 78(2):169 177, 1997.
Pravin M Vaidya. Speeding-uy up linear programming using fast matrix multiplication. In 30th annual sympo- sium on foundations of computer science, pp. 832-337. IEEE Computer Society, 1989.
Stephen Vavasis and Yinyu Ye. An accelerated interior point method whose running time depends only In Proceedings of the twenty-sirth annual ACM symposium on Theory of Computing, pp. 512-521,
1994.
Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
Johannes CW Wiesel. Measuring association with Wasserstein distances. Bernoulli, 28(4):2816 22022228222222 2022.
Yiling Xie, Yiling Luo, and Xiaoming Huo. An accelerated stochastic algorithm for solving the optimal transport problem. arXiv preprint arXiv:2203, 008113, 2022.
Quey-Jen Yeh. A reduced dual affine scaling algorilhm for solving assignament and transportation problems.
Columbia University, 1989.
Haibin Zhu, MengChu Zhou, and Rob Alkins, Group role assignment via a Kuhn Munkres algorithm- based solution. IEEE Transactions on Systems, Man, and Cybernetics-Pare A: Systems and Humans, 42(3): 739 750, 2011.
Haibin Zhu, Dongning Liu, Siqin Zhang, Yu Zhu, Luyao Teng, and Shaohua Teng. Solving the many to many assignment problem by improving the Kuhn-Mann Munkres algorithm with backtracking. Theoretical Computer Science, 618:30 41, 2016.
A Appendix
A.1 Proof of Proposition
Proof. We first consider the following two optimization problems
g
1
min ΣExech u' { -
1
Vi,j
1,...
(8)
2=1
where C is an m x m matrix generated by duplicating the jth column of C my times:
Ca 1 mi,
c
Cj mi + ... + m5-1 1st m1 + ... + 2 5
min xtett
Σ>
ot u¹ = xxx xx/=100 10. =0.11
Σ
X 1, Vi,j = 1,
x =
(9)
=1
1=1
1=1
18

---

Published in Transactions on Machine Learning Research (02/2023)
Then, denote the objective functions of problems (x()), respectively.
2
33
and 9
by f'(X), f(X), f'(X', and
Firstly, we prove
On one hand, for any X EU', if we let.
X
Σ
X ibs
j == 1,
t=1
X/J
Σ
it
2 S n
+m--+41
then we have
X/3 0,
X
m1
X
m1 tu
Σ
X
X
1
m
1=1
1=1
(=m) m1
1=1
m1
1=1,
m
t=1
t=
mh +m
Σ
X
Σ
Σ
X
m4
2 2 S n.
1=1
1+1
Thus, X' EU'
For the objective functions, have the following:
M
m m
f(x'= = Σ X C.j = Σ
X +
Σ
X1OO
I=
+my-+++
On the other hand, for any X' E U, if we let
X/1/mi 1 1515 m
x1
Xly/m3 mi + m + 1sts my
25,
then we have
X 0,
X
n
1
m,
X1J
m,
m
1=1
8
n
x1
X1
my
my
m
Thus, X Eu!
For the objective functions, we have the following:
X Cum, fXX)..
-xx)= >
-x
my
1=1==1
Hence,
By Birkhoff's theorem, know
Therefore, we have
2
19

---

Published in Transactions on Machine Learning Research (02/2023)
Similarly, for any X ut if we let
Xij Σ
X ib
j j= = 1,
t=1
Xij
Σ
X 2 j S n
1= +m--1+1
then we have X E U and f(xt) = f(X).
For any X U, if we let
Xai/mi, 15 1 S mi (is/mm
X
m1 + mj-1 + 1 St m, + .... + my 2 2 S
then we have X u and (x(x)= f(X). Therefore,
In conclusion, we have 13
A.2 Proof of Theorem 1
Proof. Denote the edge e by e (ev exg)- Let PM' be any perfect pseudo-matching in G (not necessarily in the equality graph Ei). And vf, 1,.., ,m; vig.j 1,... n nodes from V1 and V2, respectively. Since of E Vi is covered exactly once by PM', and vi E V2 covered exactly my times by PM', have
w PM) = Σ w(c) Σ ((eg))
Σ ((t)) + Σ m
eEPM
cEPM'
where the first inequality comes from the definition of feasible labeling.
Thus, En ((vf) + +E=1 myl(e)) the upper bound of the weight of any perfect pseudo-matching. Then let PM be perfect gaeudo-matchin in the equality graph Ed, we have
w(PM)
Σ w(e) = Σ ((0)) + myllus).
dEPM
Hence wPPM') S wPPM), and PM is the maximum weighted pseudo-matching.
A.3 Proof of Proposition| 2
Proof. We denote the objective functions of problems and 6 On one hand, for any Xt E Ut, if we let
by g'("),, and g (xt), respectively.
xt
== 1,
1j
oj
1=1
X =
Σ
xt tj:
25is m
(=n1 to
1+1
then we have
X1J 0,
20

---

Published in Transactions on Machine Learning Research (02/2023)
M
Σ
X3 =
Σ
Σ
X
X
m2 M
1=1
t=1
[=
t=1
11
Σ
Xij
X
71 M
i=1,
j=1
/= 1=
70 +0
X
Σ
X
n
2 m.
m
Thus, X u
For the objective function, we have the following:
M
M =ΣΣ j=
((()) = Σ C
Σ
X1, +
E
=1j=1
1= +8
On the other hand, for any X Ut if we let
(i//mi, 1sts n1
x
Xi//mi ni + ... n4 +1lsssn
n, siis
then we have
X 2 0,
X
X
π4
M
j=1
X
m
T4
2
X
Thus, Xt ET.
For the objective function, we have the following:
M
X Cun; = (X²).
g (x
=
1=1 1=1
Hence,
A.4 Additional experiment results
dependent case,
dependent case, p
independent case, p 1
independent case, p
|n(#operation)
In(sample size)
Intsample
Intsample size)
Intsample size)
medffie Mangadar sanrage
modt Hurgerinn
midite Haivssrin ware
Hampgritn awrrage
- Honggrra herd
Hungaiin Aee!
Figure 19: Comparison with the Hungarian algorithm DOT-benchmark with 32 32 resolution w.r.t. numerical operations
21

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, =
dependent case, p=2
Independent case, p=
Independent case, p
In(#operation)
in(sample size)
In(sample size)
Intsample size)
In(sample size)
modfice Rangetan vestage
modl Murguman
madh Pangartn ward
Mangerian overage
Pangaa on Des
Her punannnsss
Figure 21: Comparison with the Hungarian algorithm DOT-benchmark with 128 x 128 resolution w.r.t. numerical operations
dependent case, p
dependent case, p=
independent case, p
Independent case, p
In@operatioo)
in(sample size)
In(sample size)
intsample size)
In(sample size)
modfice Hangetar aeatage
modi Hurguran
medi Pangartin ward
Mangerian overage
-.... Hangar on Dei
Hor guran nons
Figure 22: Comparison with the Hungarian algorithm DOT-benchmark with 256 x 256 resolution w.r.t. numerical operations
dependent case, p
dependent case, p=
independent case, p
independent case, p
31+
Ini ooperation)
10
55 In(sample size)
12 Intsample size)
in(sample size)
In(sample size)
Aangadbr
Hangarian
Figure 20: Comparison with the Hungarian algorithm on DOT-benchmark with 64 x 64 resolution w.r.t.
numerical operations
22

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p=2
independent case, p
independent case, p
intomme)
33
45
to
26
23 25
35
In(sample size)
In(sapple size)
In(sample size)
Inammpeee size)
langarbe
Figure 23: Comparison with the Hungarian algorithm on DOT-benchmark with 32 x 32 resolution w.r.t. running time
dependent case, p
dependent case, =
independent case, p
independent case, p
in(tlme)
5o
69
25
24 20 14
in(sample size)
In(sample size)
In(sample size)
In(sample size)
Figure 24: Comparison with the Hungarian algorithm on DOT-benchmark with 64 x 64 resolution w.r.t. running time
dependent case, p=
dependent case, p=2
Independent case, p
Independent case, p
In(tlme)
20
23
39
20 25
23
24
22
20
24
In(sample size)
In(sample size)
In(sample size)
In(sample size)
langartt
Figure 25: Comparison with the Hungarian algorithm on DOT-benchmark with 128 x 128 resolution w.r.t. running time
dependent case, p=
dependent case, p=
independent case, p
Independent case, p p=2
(n(time)
22 Figure 26: Comparison with the Hungarian algorithm on DOT-benchharrk with 256 x 256 resolution w.r.t. running time
50
22
20 In(sample size)
25
25 20 15
in(sample size)
In(sample size)
In(sample size)
23

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p
dependent case, p
Independent case, p
independent case, p
223
223
815-
223-
In(#operation
153
352
12.
5.0
In(sample size)
In(sample size)
intsample size)
In(sample size)
Figure 27: Comparison with the Sinkhorn algorithm on synthetic data w.r.t. numerical operations
dependent dependent case,
independent case,
independent p
Int#operationi
In(sample size)
23 In(sample size)
Intsample size)
In(sample size)
modified
medifted
Figure 28: Comparison with the Sinkhorn algorithm on CIFARIO w.r.t. numerical operations
dependent case,
dependent case,
independent case,
independent case,
Int#operation)
In(sample size)
Intsample
Intsample size)
In(sample size)
mocMec langararr avarage
medtt Hargeran zeer
modifled unngadin More
#itccom Sverage
disenom see
Figure 29: Comparison with the Sinkhorn algorithm Wisconsin cancer data w.r.t. numerical operations
dependent case,
dependent case, p
independent case, p
independent case,
In(#operation))
In(sample size)
intsample size)
Intsample size)
Intsample size)
mocflie Nangeriar average
modit od Hargaran
medifed trungarian wast
allkron average
strenorn 240t ... alntrommmmi
Figure 30: Comparison with the Sinkhorn algorithm DOT-benchmark with 32 x 32 resolution w.r.t. numerical operations
dependent case, p
dependent case, p
independent case, p
independent case, p
In(#operation)
15.2
Intsample size)
In(sample size)
intsample size)
intsample size)
Hanprror average
Hungarion
Seerage
Figure 31: Comparison with the Sinkhorn algorithm on DOT-benchmark with 64 x 64 resolution w.r.t. numerical operations
24

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case,
dependent case, p=2
independent case,
pendent case, p
222-
Int In(#operation
16.0
26 Intsample size)
12
20 In(sample size)
In(sample size)
in(sample size)
Figure 32: Comparison with the Sinkhorn algorithm on DOT-benchmark with 128 x 128 resolution w.r.t.
numerical operations
dependent case, p
dependent case, p=
independent case. pase.p=
independent case, p
Ini ooperratioon
16 in(sample size)
16
15 In(sample size)
In(sample size)
Intsample size)
Figure 33: Comparison with the Sinkhorn algorithm on DOT-benchmark with 256 x 256 resolution w.r.t. numerical operations
dependent case, p=
dependent case,
Independent p p=1
Independent case, =
Int #operatioo))
20
25
23
26 Intsample size)
24 24
in(sample size)
In(sample size)
In(sample size)
Figure 34: Comparison with the Sinkhorn algorithm on DOT-benchmark with 512 x 512 resolution w.r.t. numerical operations
dependent case, p=
dependent case, p=2
Independent case, p=
Independent case, p
(n(time)
50 30 in(sample size)
25 20 14
In(sample size)
In(sample size)
In(sample size)
Figure 35: Comparison with the Sinkhorn algorithm on DOT-benchmark with 32 x 32 resolution w.r.t.
running time
25

---

Published in Transactions on Machine Learning Research (02/2023)
dependent case, p=
dependent case, p=2
independent case, p
independent case, p=
intomm)
50
25 23 16 1
in(sample size)
Isaample size)
In(sample size)
Ini (samplle size)
Figure 36: Comparison with the Sinkhorn algorithm on DOT-benchmark with 64 x 64 resolution w.r.t. running time
dependent case, p
dependent case,
independent case, p
independent case, p
(n(tlme)
25 25 20 14 In(sample size)
in(sample size)
In(sample size)
In(sample size)
Figure 37: Comparison with the Sinkhorn algorithm on DOT-benchmark with 128 x 128 resolution w.r.t.
running time
dependent case, ==
dependent case,
Independent case, =
Independent case, p=2
In(tlme)
-200
+00
-10--
29
25
20
25
30
29
25
In(sample size)
In(sample size)
Intsample size)
In(sample size)
Figure 38: Comparison with the Sinkhorn algorithm on DOT-benchmark with 256 x 256 resolution w.r.t.
running time
dependent case, p=
dependent case, p=2
independent case, p=]
Independent case,
(n(time)
20 in(sample size)
20 30
22
20 In(sample size)
25
25
20
15
In In(sample size)
In(sample size)
Figure 39: Comparison with the network simplex algorithm on NOT-benchmark with 32 32 resolution w.r.t. running time
26

---

Published in Transactions on Machine Learning Research (02/2023)
dependent p
dependent case, p
independent p
independent p
Inttime)
In(sample size)
Inaaplle size)
In(sample size)
24
In(smplle size)
Figure 40: Comparison with the network simplex algorithm on DOT-benchmark with 64 64 resolution
w.r.t. running time
dependent case,
dependent case, 2
independent case,
independent case,
in(tlme)
In(sample size)
(sampll size)
In(sample size)
In(smmlle
modffec Nurgaran everage
med Fedd liungar an best
mocfled Hangarian word
Velwork timpies avarage
Retwork dinillor best
. ketwort dirp
Figure 41: Comparison with the network simplex algorithm w.r.t. running time
with 128 128 resolution
dependent case, p
dependent case, 2
independent case,
Independent case, p
Intomme]
15
in(sample size)
inaamlle size)
In(sample size)
In(sample size)
Hunggr
Figure 42: Comparison with the network simplex algorithm on DOT-benchmark with 256 x 256 resolution w.r.t running time
27


Reviews:
[{"review_id": "b6xT4auFmf", "reviewer": "Reviewer_F8ov", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n- clear presentation of a simple but effective idea\n- theoretical results along with numerical experiments\n- comparison with related work\n\nWeakness:\n- the generalization section might be a bit confusing", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper has a simple and clear contribution: a modification of Hungarian algorithm to solve a special class of optimal transport  (OT) problems. This is useful, for example, for the application of Wasserstein distance for independence test, reducing the computational complexity from $\\mathcal O(n^6)$ to $\\mathcal O(n^5)$. The special class of OT problems include the ones defined for empirical measures which can be equivalently expressed as a (pseudo) matching problem. The computational complexity is theoretically proven and illustrated on several numerical examples.  ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Overall, I like the paper and I think it is making a solid contribution. My only minor request for change is about the generalization section. The proposed idea is to increase the number of atoms for one marginal from $n$ to $M$ by creating $n_i$ atoms for each node. This seems in contrast to the original motivation of the paper, which was about exactly not performing this action on the other marginal. There should be a better way to do this. So I am not a fan of calling this a \"generalization\" of your proposed approach. ", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ofFHE7TyNt", "reviewer": "Reviewer_zbUw", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- despite the exposition (that needs to be re-organized in my opinion), the paper is globally clear and well written\n- the computational improvement is significant\n- the topic is of interest to the TMLR community\n\n**Weaknesses**\n- as mentioned above, I feel that the organization of this paper should be revised. If I can understand that the Wassertstein independence test is one of the main motivations to the present contribution, I would not center the exposition around it. In particular, by Section 4 it seems that the proposed algorithm can be applied to a class of OT problems that is way more general than independence tests. For instance, if we have two empirical distributions composed respectively of $m$ and $n$ atoms, cannot we choose $M = mn$, $n_i=m$ and $m_j=n$, such that computing the Wasserstein distance between the two distributions is actually a special case of Problem (4)?\n- I think it would benefit the paper to provide some reminders on the non-modified Hungarian algorithm, for the reader to be able to appreciate the modification and see where it plays a role in reducing computations. This point is far from being clear at the moment.\n- the one-to-many and many-to-many assignment problems are presented as two different applications although the latter seems to be a strict generalization of the first one, or am I missing something here?\n- in Example 1, the computation complexities $\\mathcal{O}(4m_1^2)$ and $\\mathcal{O}(5m_1^2)$ do not make sense, as 4 and 5 are constants already taken into account by the $\\mathcal{O}$ notation. Instead, it would make more sense to present a generalized setting with an arbitrary number $r$ of roles and present computational complexities that depend on $r$.\n- regarding the experiments, why do the authors present experiments with and without dependence? Is it supposed to influence the algorithm's speed? Same question for $p$. It seems to me that once the cost matrix $C$ is computed, the algorithm is oblivious to anything else. Could the authors explain the big variations, when any, between the best and worst run? When benchmarking Sinkhorn, I would like to see the precision achieved as a function of the number of operations. Indeed, the algorithm might attain a reasonable precision very quickly, which makes it an interesting tool if one is not interested in solving the problem exactly.\n- the related work section is very brief and do not provide any details about what is done in the cited articles", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a modified Hungarian algorithm to solve exactly Integer Linear Programs of the form described in Equation (3). These optimization problems are equivalent to find a pseudo-matching (Definition 8) on some equality graph (Theorem 1). The algorithm proceeds by iteratively improving the labeling and the pseudo-matching (Algorithm 1). The authors then detail 3 applications: an independence test using the Wasserstein distance, the one-to-many assignment problem, and the many-to-many assignment problem, for which the proposed algorithm improves over existing methods in terms of computational time. Some numerical experiments are presented, that confirm the computational superiority of the proposed algorithm (Section 7).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Address the weaknesses detailed above. In particular:\n- clarify the exposition and detail the broadest setting in which the proposed algorithm can be applied\n- provide reminders about the Hungarian algorithm to make the modification explicit\n- clarify the difference between one-to-many and many-to-many assignment problems\n- clarify the computational complexities\n- clarify the experiments\n- clarify the related work section\n\nOther points:\n- in the introduction $\\gamma(\\nu_1, \\nu_2)$ is never defined, nor $\\hat{\\nu}_1$ and $\\nu_2$. This should be clarified\n- Proposition 2 should appear explicitly in the text\n- the assignment problem should be defined when discussed right after eq.(3), or at least a pointer given", "broader_impact_concerns": "No concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fhCyd8xUiW", "reviewer": "Reviewer_ayht", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**:\n* The idea of using Wasserstein distances for computing independence test is interesting and relevant to the ML community (but it is not novel)\n* The proposed algorithm is better than a naïve application of the Hungarian algorithm\n* The computational results look reproducible by a knowledgeable and coding-skilled reader\n\n**Weakness**:\n* The paper has flaws in the presentation and writing. They never cite the semi-assignment problem, which was studied in several papers from the ’90, see for example:\n\n1. Jonker, R. and Volgenant, A., 1987. A shortest augmenting path algorithm for dense and sparse linear assignment problems. Computing, 38(4), pp.325-340.\n\n2. Kennington, J. and Wang, Z., 1992. A shortest augmenting path algorithm for the semi-assignment problem. Operations Research, 40(1), pp.178-187.\n\n3. Balas, E., Miller, D., Pekny, J. and Toth, P., 1991. A parallel shortest augmenting path algorithm for the assignment problem. Journal of the ACM (JACM), 38(4), pp.985-1004.\n\n* The presentation at page 2 in the introduction is imprecise. The problem the authors are presenting is a semi-assignment problem, which is a special case of minimum cost flow problem, and which can be solved by linear programming using simplex-based algorithms (worst case exponential) and/or interior point algorithms (polynomially solvable). In particular, the minimum cost flow problem is strongly polynomially solvable by primal Network Simplex algorithms, see for example:\n4. Orlin, J.B., 1997. A polynomial time primal network simplex algorithm for minimum cost flows. Mathematical Programming, 78(2), pp.109-129.\n\n* In addition, interior point algorithms were customized for bipartite min cost flow problems:\n\n5. Resende, M.G. and Veiga, G., 1993. An implementation of the dual affine scaling algorithm for minimum-cost flow on bipartite uncapacitated networks. SIAM Journal on Optimization, 3(3), pp.516-537.\n\n* About the references on Optimal Transport algorithms (i.e. “Approximation algorithms” in related work section), it would be better to distinguish among algorithms that compute the optimal solution value and the optimal decision variables values (i.e., the optimal transportation plan), and the algorithms that compute or approximate the optimal solution value without actually computing a feasible transportation plan. In the two cases, the (pseudo)polynomial worst-case complexity is different, with the second type of algorithm having a better worst-case complexity.\n\n* The list of definitions in section 2 can be improved for clarity. For notations and definitions of these types of networks flow problems, we suggest referring to: \n6. Ahuja, R.K., Magnanti, T.L. and Orlin, J.B., 1988. Network flows.\n\n* The instances used in the benchmark are not very large, and we would suggest also using the standard DOT-benchmark to different image resolutions to study the scalability of the proposed approach:\n7. Schrieber, J., Schuhmacher, D. and Gottschlich, C., 2016. Dotmark–a benchmark for discrete optimal transport. IEEE Access, 5, pp.271-282.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper shows how the Wasserstein distance-based independence test between two empirical discrete probability measures boils down to the solution of a semi-assignment problem. \n\nThe authors claim to have a novel modification of the Hungarian algorithm for computing minimum cost bipartite matching with worst-case complexity of $O(n^5)$, which is better than applying the Hungarian algorithm with a complexity of $O(n^6)$, where n is the number of i.i.d. samples being tested for independence.\n\nThe proposed algorithm is implemented in Python and compared with the Sinkhorn-based algorithm on the CIFAR10 and the Wisconsin breast cancer datasets. In the proposed instances, the runtime of their algorithm looks modestly interesting.\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The main request for changes is concerned with better positioning of the presented optimal transport problem in the cost of existing network problems.\n\nThe best sequential algorithm to solve the semi-assignment problem is presented in Kennington, J. and Wang, Z., 1992 cited before, which has the same worst-case complexity as the algorithm proposed in this paper. Moreover, that algorithm was further improved in parallel implementations in subsequent papers from other authors.\n\nGiven that the algorithm proposed in this paper does not have a better worst-case complexity, the authors should compare their approach with the algorithm in Kennington, J. and Wang, Z., 1992.\n\nIf they believe the effort to re-implement their algorithm is too significant, at least they should compare with the solution of the semi-assignment problem via the Network Simplex as implemented in the COIN-OR Lemon Graph library, which in the paper on arxiv by Dong et. all 2020 (that you cite) is clearly one of the best implementation to solve minimum cost network flow problems.\n\nMore importantly, you should clearly explain how your algorithm differs from the Kennington, J. and Wang, Z., 1992 algorithm because it does not look too different.\n\nSection 4 does not add too much value to the paper, and it could be merged with the introduction.\n\nOther requests for changes are implicit in the definition of the weakness of the paper (previous section).\n\nMinor comments:\n- In problem formulation (1), in the definition of PI^o is unclear the quantification over the constraints (for all ij in … ?)\n- In problem (2), are the m_j variable or parameters? The last constraint on the summation of m_j does not look like a constraint but an assumption on the data of the problem.\n- Problem (3) is known in the literature as a semi-assignment problem and in general, can be formulated and solved as a minimum cost network flow problem (see the related citation)\n- We suggest including the plots of the computational results with the running time in the paper, and postponing in the appendix the results with the number of operations. The number of operations can be a misleading measure of complexity. For instance, the Sinkhorn algorithm is easily parallelizable, while primal-dual algorithms, like the  Hungarian algorithm, are much harder to parallelize and do not get the same speed-up gain if parallelized.\n- The many-to-many problem is again a network flow problem, and it can be solved by the network simplex algorithm (or other specialized minimum cost network flow problems).\n", "broader_impact_concerns": "We have no major concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "PTW2kEYVPx", "reviewer": "Reviewer_F217", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strength of the work is that the proposed Hungarian algorithm indeed has a good complexity for this problem.\n\nThe weaknesses are that the paper is not well-written, and the comparison to the literature is inadequate. Indeed, in the (short) related work section, the proposed approach is only compared to low-accuracy OT solvers, but it seems like a more appropriate comparison is with the combinatorial algorithms described in Peyré et al. (2019) (by the way, the bibliography entry for this reference is incorrect, as it says Peyré, Cuturi, et al., but I am only aware of two authors total), sections 3.5-3.7. In particular, as mentioned at the end of section 3.5, it is already known that the network simplex algorithm only requires $\\widetilde O(m^2 n)$ complexity. Without a detailed comparison to prior work, the submission cannot be accepted for publication.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors prove computational complexity guarantees for a modified Hungarian algorithm to solve a special type of discrete optimal transport problem. They showcase applications to independence testing and other assignment problems.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As discussed above, please conduct a more thorough literature search and provide a detailed comparison with the network simplex algorithm, among other algorithms; adjust the claims accordingly.\n\nAlso, the phrase “w.l.o.g.” is used incorrectly throughout. To be clear, WLOG means we can focus on a special case, because the general case is reducible to this special case. In this paper, WLOG seems to be used to mean “we impose this assumption as this is the situation that we care about”.\n\nThe notation $\\gamma(\\nu_1, \\nu_2)$ is never defined; presumably it refers to the joint distribution of $(Y, Z)$. The notation does not make sense; the joint distribution is not a “function” of the marginal distributions. Moreover it is confusing as it is easily confused for the optimal transport plan with marginals $\\nu_1$ and $\\nu_2$.\n", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper makes simple yet clear and concrete algorithmic contributions to a broad class of problems. There were some concerns initially regarding the presentation and organization of the paper and comparison with the related work. The authors have addressed them adequately in the revision they posted. 
Decision: Accept (poster)