Year: 2021

Abstract:
This paper studies the equitable and optimal transport (EOT) problem, which has many applications such as fair division problems and optimal transport with multiple agents etc. In the discrete distributions case, the EOT problem can be formulated as a linear program (LP). Since this LP is prohibitively large for general LP solvers, Scetbon et al. Scetbon et al. (2021) suggests to perturb the problem by adding an entropy regularization. They proposed a projected alternating maximization algorithm (PAM) to solve the dual of the entropy regularized EOT. In this paper, we provide the first convergence analysis of PAM. A novel rounding procedure is proposed to help construct the primal solution for the original EOT problem. We also propose a variant of PAM by incorporating the extrapolation technique that can numerically improve the performance of PAM. Results in this paper may shed lights on block coordinate (gradient) descent methods for general optimization problems.

text: This paper studies the equitable and optimal transport (EOT) problem, which has many applications such as fair division problems and optimal transport with multiple agents etc. In the discrete distributions case, the EOT problem can be formulated as a linear program (LP). Since this LP is prohibitively large for general LP solvers, Scetbon et al. Scetbon et al. (2021) suggests to perturb the problem by adding an entropy regularization. They proposed a projected alternating maximization algorithm (PAM) to solve the dual of the entropy regularized EOT. In this paper, we provide the first convergence analysis of PAM. A novel rounding procedure is proposed to help construct the primal solution for the original EOT problem. We also propose a variant of PAM by incorporating the extrapolation technique that can numerically improve the performance of PAM. Results in this paper may shed lights on block coordinate (gradient) descent methods for general optimization problems.

heading: 1 INTRODUCTION
text: Optimal transport (OT) is a classical problem that recently finds many emerging applications in machine learning and artificial intelligence, including generative models Arjovsky et al. (2017), representation learning Ozair et al. (2019), reinforcement learning Bellemare et al. (2017) and word embeddings Alvarez-Melis et al. (2019) etc. More recently, Scetbon et al. Scetbon et al. (2021) proposed an equitable and optimal transport (EOT) problem that targets to fairly distribute the workload of OT when there are multiple agents. In this problem formulation, there are multiple agents working together to move mass from measures µ to ν and each agent has its unique cost function. A very important issue that needs to be considered here is the fairness, which aims at finding transportation plans such that the workloads among all the agents are equal to each other. This can be achieved by minimizing the largest transportation cost among all agents, which leads to a convex-concave saddle point problem. The EOT problem has wide applications in economics and machine learning, such as fair division or the cake-cutting problem Moulin (2003); Brandt et al. (2016), multi-type resource allocation Mackin & Xia (2015), internet minimal transportation time and sequential optimal transport Scetbon et al. (2021). We now describe the EOT problem formally. Given two discrete probability measures µn =∑n i=1 aiδxi and νn = ∑n i=1 biδyi , the EOT studies the problem of transporting mass from µ to ν by N agents. Here, {x1, x2, ..., xn} ⊂ Rd and {y1, y2, ..., yn} ⊂ Rd are the support points of each measure and a = [a1, a2, ..., an]⊤ ∈ ∆n, b = [b1, b2, ..., bn]⊤ ∈ ∆n are corresponding weights for each measure, where ∆n denotes the probability simplex in Rn. Moreover, throughout this paper, we assume vector b > 0. For each agent k, we denote its unique cost function as ck(x, y), k ∈ [N ] = {1, . . . , N} and its cost matrix as Ck, where Cki,j = ck(xi, yj). Moreover, we define the following coupling decomposition set
ΠNa,b := { π = (πk)k∈[N ] ∣∣∣∣∣ r (∑
k
πk ) = a, c (∑ k πk ) = b, πkij ≥ 0,∀i, j ∈ [n] } ,
where r(π) = π1, c(π) = π⊤1 are the row sum and column sum of matrix π respectively. Mathematically, the EOT problem can be formulated as
min π∈ΠNa,b max 1≤k≤N
⟨πk, Ck⟩. (1)
When N = 1, (1) reduces to the standard OT problem. Note that (1) minimizes the point-wise maximum of a finite collection of functions. It is easy to see that (1) is equivalent to the following constrained problem:
min π∈ΠNa,b max λ∈∆N+ ℓ(π, λ) := N∑ k=1 λk⟨πk, Ck⟩. (2)
The following proposition shows an important property of EOT: at the optimum of the minimax EOT formulation (2), the transportation costs of the agents are equal to each other.
Proposition 1 (Scetbon et al., 2021, Proposition 1) Assume that all cost matrices Ck, k ∈ [N ] have the same sign. Let π∗ ∈ ΠNa,b be the optimal solution of (2). It holds that
⟨(π∗)i, Ci⟩ = ⟨(π∗)j , Cj⟩, ∀i, j ∈ [N ]. (3)
Note that Proposition 1 requires all cost matrices to have the same sign. When the cost matrices are all non-negative, (2) solves the transportation problem with multiple agents. When the cost matrices are all non-positive, the cost matrices are interpreted as the utility functions and (2) solves the fair division problem Moulin (2003).
The discrete OT is a linear programming (LP) problem (in fact, an assignment problem) with a complexity of O(n3 log n) Tarjan (1997). Due to this cubic dependence on the dimension n, it is challenging to solve large-scale OT in practice. A widely adopted compromise is to add an entropy regularizer to the OT problem Cuturi (2013). The resulting problem is strongly convex and smooth, and its dual problem can be efficiently solved by the celebrated Sinkhorn’s algorithm Sinkhorn & Knopp (1967); Cuturi (2013). This strategy is now widely used in the OT community due to its computational advantages as well as improved sample complexity Genevay et al. (2019). Similar ideas were also used for computing the Wasserstein barycenter Benamou et al. (2015), projection robust Wasserstein distance Paty & Cuturi (2019); Lin et al. (2020); Huang et al. (2021a), projection robust Wasserstein barycenter Huang et al. (2021b). Motivated by these previous works, Scetbon et al. Scetbon et al. (2021) proposed to add an entropy regularizer to (2), and designed a projected alternating maximization algorithm (PAM) to solve its dual problem. However, the convergence of PAM has not been studied. Scetbon et al. Scetbon et al. (2021) also proposed an accelerated projected gradient ascent algorithm (APGA) for solving a different form of the dual problem of the entropy regularized EOT. Since the objective function of this new dual form has Lipschitz continuous gradient, APGA is essentially the Nesterov’s accelerated gradient method and thus its convergence rate is known. However, numerical experiments conducted in Scetbon et al. (2021) indicate that APGA performs worse than PAM. We will discuss the reasons in details later.
Our Contributions. There are mainly three issues with the PAM and APGA algorithms in Scetbon et al. (2021), and we will address all of them in this paper. Our results may shed lights on designing new block coordinate descent algorithms. Our main contributions are given below.
(i) The PAM algorithm in Scetbon et al. (2021) only returns the dual variables. How to find the primal solution of (2), i.e., the optimal transport plans π , was not discussed in Scetbon et al. (2021). In this paper, we propose a novel rounding procedure to find the primal solution. Our rounding procedure is different from the one widely used in the literature Altschuler et al. (2017).
(ii) We provide the first convergence analysis of the PAM algorithm, and analyze its iteration complexity for finding an ϵ-optimal solution to the EOT problem (2). In particular, we show that it takes at most O(Nn2ϵ−2) arithmetic operations to find an ϵ-optimal solution to (2). This matches the rate of the Sinkhorn’s algorithm for computing the Wasserstein distance Dvurechensky et al. (2018).
(iii) We propose a variant of PAM that incorporates the extrapolation technique as used in Nesterov’s accelerated gradient method. We name this variant as PAM with Extrapolation (PAME). The iteration complexity of PAME is also analyzed. Though we are not able to prove a better complexity over PAM at this moment, we find that PAME performs much better than PAM numerically.
Notation. For vectors a and b with the same dimension, a./b denotes their entry-wise division. We denote c∞ := maxk ∥Ck∥∞. Throughout this paper, we assume vector b > 0, and we denote ι := minj log(bj). We use 1n to denote the n-dimensional vector whose entries are all equal to one. We use IX (x) to denote the indicator function of set X , i.e., IX (x) = 0 if x ∈ X , and IX (x) = ∞ otherwise. We denote ct = c( ∑N k=1 π
k(f t+1, gt, λt)). For integer N > 0, we denote [N ] := {1, . . . , N}. We also denote π(f, g, λ) = [πk(f, g, λ)]k∈[N ].

heading: 2 PROJECTED ALTERNATING MAXIMIZATION ALGORITHM
text: The PAM algorithm proposed in Scetbon et al. (2021) aims to solve the entropy regularized EOT problem, which is given by
min π∈ΠNa,b max λ∈∆N+ ℓη(π, λ) := N∑ k=1 pkη(π k, λ) (4)
where η > 0 is a regularization parameter, pkη(π k, λ) := λk⟨πk, Ck⟩ − ηH(πk), and the entropy function H is defined as H(π) = − ∑
i,j πi,j(log πi,j − 1). The entropy regularization was first introduced into the OT problem by Cuturi (2013) and is now widely used in the OT community. By adding an entropy regularizer, the primal problem becomes strongly convex and the dual problem is unconstrained and is suitable for alternating maximization. This leads to the Sinkhorn’s algorithm which has low per-iteration complexity and thus is scalable. The PAM algorithm proposed by Scetbon et al. (2021) used the same idea for the EOT problem. Note that (4) is a strongly-convexconcave minimax problem whose constraint sets are convex and bounded, and thus the Sion’s minimax theorem Sion (1958) guarantees that
min π∈ΠNa,b max λ∈∆N+ ℓη(π, λ) = max λ∈∆N+ min π∈ΠNa,b ℓη(π, λ). (5)
Now we consider the dual problem of minπ∈ΠNa,b ℓη(π, λ). First, we add a redundant constraint∑ k,i,j π k i,j = 1 and consider the dual of
min π∈ΠNa,b, ∑ k,i,j π k i,j=1 ℓη(π, λ). (6)
The reason for adding this redundant constraint is to guarantee that the dual objective function is Lipschitz smooth. It is easy to verify that the dual problem of (6) is given by
max f,g min∑ k,i,j π k i,j=1,
π∈(Rn×n+ ) N
N∑ k=1
λk⟨πk, Ck⟩−ηH(π)+f⊤ ( a− r (∑ k πk )) +g⊤ ( b− c (∑ k (πk)⊤ )) ,
(7) where f and g are the dual variables and H(π) = ∑ k H(π
k). It is noted that problem (7) admits the following solution:
πk(f, g, λ) = ζk(f, g, λ)∑
k ∥ζk(f, g, λ)∥1 , ∀k ∈ [N ], (8)
where
ζk(f, g, λ) = exp
( f1⊤n + 1ng ⊤ − λkCk
η
) , ∀k ∈ [N ]. (9)
By plugging (8) into (7), we obtain the following dual problem of (6):
max f∈Rn, g∈Rn
⟨f, a⟩+ ⟨g, b⟩ − η log
( N∑
k=1
∥ζk(f, g, λ)∥1 ) − η. (10)
Plugging (10) into (5), we know that the entropy regularized EOT problem (4) is equavalent to a pure maximization problem:
max f∈Rn, g∈Rn, λ∈∆N
F (f, g, λ) := ⟨f, a⟩+ ⟨g, b⟩ − η log
( N∑
k=1
∥ζk(f, g, λ)∥1 ) − η. (11)
Function F (f, g, λ) is a smooth concave function with three block variables (f, g, λ). We use (f∗, g∗, λ∗) to denote an optimal solution of (11), and we denote F ∗ = F (f∗, g∗, λ∗). The PAM algorithm proposed in Scetbon et al. (2021) is essentially a block coordinate descent (BCD) algorithm for solving (11). More specifically, the PAM updates the three block variables by the following scheme:
f t+1 ∈ argmax f F (f, gt, λt), (12a)
gt+1 ∈ argmax g F (f t+1, g, λt), (12b)
λt+1 := Proj∆N ( λt + τ∇λF (f t+1, gt+1, λt) ) . (12c)
Each iteration of PAM consists of two exact maximization steps followed by one projected gradient step. Importantly, the two exact maximization problems (12a)-(12b) have numerous optimal solutions, and we choose to use the following ones:
f t+1 = f t + η log  a r (∑N
k=1 ζ k(f t, gt, λt)
)  , (13)
gt+1 = gt + η log  b c (∑N
k=1 ζ k(f t+1, gt, λt)
)  . (14)
Furthermore, the optimiality conditions of (12a)-(12b) imply that a− r (∑N k=1 ζ k(f t+1, gt, λt) ) ∑
k ∥ζk(f t+1, gt, λt)∥1 = 0, b−
c (∑N
k=1 ζ k(f t+1, gt+1, λt) ) ∑
k ∥ζk(f t+1, gt+1, λt)∥1 = 0, ∀t. (15)
However, we need to point out that the PAM (12) only returns the dual variables (f t, gt, λt). One can compute the primal variable π using (8), but it is not necessarily a feasible solution. That is, π computed from (8) does not satisfy π ∈ ΠNa,b. How to obtain an optimal primal solution from the dual variables was not discussed in Scetbon et al. (2021). For the OT problem, i.e., N = 1, a rounding procedure for returning a feasible primal solution has been proposed in Altschuler et al. (2017). However, this rounding procedure cannot be applied to the EOT problem directly. In the next section, we propose a new rounding procedure for returning a primal solution based on the dual solution (f t, gt, λt). This new rounding procedure involves a dedicated way to compute the margins.

heading: 2.1 THE ROUNDING PROCEDURE AND THE MARGINS
text: Given a ∈ ∆n, b ∈ ∆n, and π = {πk}k∈[N ] satisfying r( ∑ k π k) = a, we construct vectors ak, bk ∈ Rn, k ∈ [N ] from the procedure
(ak, bk)k∈[N ] = Margins(π, a, b). (16)
The details of this procedure is given below. First, we set ak = r(πk), which immediately implies∑N k=1 a
k = a. We then construct bk such that the following properties hold (these properties are required in our convergence analysis later):
(i) bk ≥ 0; (ii) ∑N
k=1 b k = b; (iii) ∑n
i=1 a k i = ∑n j=1 b k j , ∀k ∈ [N ];
(iv) For any fixed j ∈ [n], the quantities bkj − [c(πk)]j have the same sign for all k ∈ [N ]. That is, for any k and k′, we have
(bkj − [c(πk)]j) · (bk ′ j − [c(πk ′ )]j) ≥ 0, (17)
which provides the following identity that is useful in our convergence analysis later:
N∑ k=1 ∥bk − c(πk)∥1 = N∑ k=1 n∑ j=1 |bkj − [c(πk)]j | = n∑ j=1 ∣∣∣∣∣ N∑ k=1 (bkj − [c(πk)]j) ∣∣∣∣∣ =
n∑ j=1 ∣∣∣∣∣∣bj − [ c ( N∑ k=1 πk )] j ∣∣∣∣∣∣ = ∥∥∥∥∥b− c ( N∑ k=1 πk )∥∥∥∥∥ 1 .
(18)
The procedure on constructing (bk)k∈[N ] satisfying these four properties is provided in Appendix A.
After (ak, bk)k∈[N ] are constructed from (16) with π = π(fT , gT−1, λT−1), we adopt the rounding procedure proposed in Altschuler et al. (2017) to output a primal feasible solution (π̂k)k∈[N ]. The rounding procedure is described in Algorithm 2.
With this new procedure for rounding and computing the margins ak, bk, we now formally describe our PAM algorithm in Algorithm 1. Note that the algorithm is terminated when the following criteria are met:
∥ct−1 − b∥1 ≤ ϵ/(6(6c∞ − ηι)), (19a)∥∥λt − λt−1∥∥ 2 ≤ ηϵ/(18c2∞), (19b) F̃ (f t, gt−1, λt−1) ≤ ϵ/6, (19c)
where F̃ (f, g, λ) is the suboptimality defined as: F̃ (f, g, λ) = F (f∗, g∗, λ∗)− F (f, g, λ).
Algorithm 1 Projected Alternating Maximization Algorithm
1: Input: Cost matrices {Ck}1≤k≤N , vectors a, b ∈ ∆n+ with b > 0, accuracy ϵ. f0 = g0 = [1, ..., 1]⊤, λ0 = [1/N, ..., 1/N ]⊤ ∈ ∆N+ . t = 0 2: Choose parameters as
η = min
{ ϵ
3(log(n2N) + 1) , c∞
} , τ = η/c2∞. (20)
3: while (19) is not met do 4: Compute f t+1 by (13) 5: Compute gt+1 by (14) 6: Compute λt+1 by (12c) 7: t← t+ 1 8: end while 9: Assume stopping condition (19) is satisfied at the T -th iteration. Compute (ak, bk)k∈[N ] =
Margins(π(fT , gT−1, λT−1), a, b) as in Section 2.1. 10: Output: (π̂, λ̂) where π̂k = Round(πk(fT , gT−1, λT−1), ak, bk), ∀k ∈ [N ], λ̂ = λT−1.
Algorithm 2 Round(π, a, b)
1: Input: π ∈ Rn×n, a ∈ Rn+, b ∈ Rn+. 2: X = Diag (x) with xi = air(π)i ∧ 1 3: π′ = Xπ 4: Y = Diag (y) with yj =
bj c(π′)j ∧ 1 5: π′′ = π′Y 6: erra = a− r(π′′), errb = b− c(π′′) 7: Output: π′′ + erraerr⊤b /∥erra∥1.

heading: 2.2 CONNECTIONS WITH BCD AND BCGD METHODS
text: We now discsuss the connections between PAM and the block coordinate descent (BCD) method and the block coordinate gradient descent (BCGD) method. For the ease of presentation, we now assume that we are dealing with the following general convex optimization problem with m block variables:
min xi∈Xi,i=1,...,m J(x1, x2, . . . , xm), (21)
where Xi ⊂ Rdi and J is convex and differentiable. The BCD method for solving (21) iterates as follows:
xt+1i = argmin xi∈Xi J(xt+11 , x t+1 2 , . . . , x t+1 i−1, xi, x t i+1, . . . , x t m), (22)
and it assumes that these subproblems are easy to solve. The BCGD method for solving (21) iterates as follows:
xt+1i = argmin xi∈Xi ⟨∇xiJ(xt+11 , x t+1 2 , . . . , x t+1 i−1, xi, x t i+1, . . . , x t m), xi − xti⟩+
1
2τ ∥xi − xti∥22, (23)
where τ > 0 is the step size. The PAM (12) is a hybrid of BCD (22) and BCGD (23), in the sense that some block variables are updated by exactly solving a maximization problem (the f and g steps), and some other block variables are updated by taking a gradient step (the λ step). Though this hybrid idea has been studied in the literature Hong et al. (2017); Xu & Yin (2013), their convergence analysis requires the blocks corresponding to exact minimization to be strongly convex. However, in our problem (11), the negative of the objective function is merely convex. Hence we need to develop new convergence proofs to analyze the convergence of PAM (Algorithm 1). How to extend our convergence results of PAM (Algorithm 1) to more general settings is a very interesting topic for future study.

heading: 3 CONVERGENCE ANALYSIS OF PAM
text: In this section, we analyze the iteration complexity of Algorithm 1 for obtaining an ϵ-optimal solution to the original EOT problem (2). The ϵ-optimal solution to (2) is defined as follows.
Definition 2 (see, e.g., Nemirovski (2005)) We call (π̂, λ̂) ∈ ΠNa,b × ∆N an ϵ-optimal solution to the EOT problem (2) if the following inequality holds:
max λ∈∆N ℓ(π̂, λ)− min π∈ΠNa,b ℓ(π, λ̂) ≤ ϵ.
Note that the left hand side of the inequality is the duality gap of (2).

heading: 3.1 MAIN RESULT
text: We now present our main theorem, which gives the iteration complexity of PAM such that (55) is satisfied, and as a result of Lemma 15, an ϵ-optimal solution to the original EOT problem (2) is obtained.
Theorem 3 Define ϵ′ = ϵ/(6c∞ − ηι), and set T as
T = 5 + 36
η √ γ0ϵ′ + 648c2∞ ηϵ + 28 ηγ0ϵ = O
( c2∞ϵ −2) , (24) where γ0 = min { 1
(2c∞−ηι)2 , 19c2∞
} is a constant and we know γ0 = O(c−2∞ ). The output pair of
Algorithm 1 is an ϵ-optimal solution of the EOT problem (2).
Proof. See Appendix B. □
Remark 4 Though our complexity result matches the rate of the Sinkhorn’s algorithm in terms of the dependence on ϵ, we argue that EOT is a more difficult problem than the entropic regularized OT, and thus our results are promising. First, EOT is a saddle-point problem while entropic regularized
OT is a minimization problem. Second, the extra variable λ in EOT requires a gradient projection step in the PAM algorithm, which introduces significant difficulty to the analysis of the convergence behavior. While for Sinkhorn’s algorithm it is much easier to analyze, because the dual is unconstrained. Third, since there are multiple agents in EOT, it is more difficult to design the rounding procedure to obtain the primal solution. We also note that the dependence of c∞ in our result and in the result of Sinkhorn’s algorithm Dvurechensky et al. (2018) are both c2∞.

heading: 4 PROJECTED ALTERNATING MAXIMIZATION WITH EXTRAPOLATION
text: In this section, we discuss how to accelerate the PAM algorithm (Algorithm 1). It can be shown that the gradient of F in (11) is Lipschitz continuous1. Therefore, Scetbon et al. Scetbon et al. (2021) proposed to adopt Nesterov’s accelerated gradient method Nesterov (2004) to solve (11). Their algorithm, named APGA (Accelerated Projected Gradient Ascent algorithm), iterates as follows:
(v, w, z)⊤ ← (f t−1, gt−1, λt−1)⊤ + t− 2 t+ 1
( (f t−1, gt−1, λt−1)⊤ − (f t−2, gt−2, λt−2)⊤ ) (25a)
(f t, gt)⊤ ← (v, w)⊤ + 1 L ∇(f,g)F (v, w, z) (25b)
(λt)⊤ ← Proj∆N ( z + 1
L ∇λF (v, w, z)
) , (25c)
where L is the Lipschitz constant of∇F . Note that APGA treats the problem (11) as a generic convex and smooth problem, and does not take advantage of the special structures of (11). In particular, f and g are updated using gradient ascent steps. This is in contrast to PAM in which f and g are obtained by exact maximizations, which is expected to improve the function value of F more significantly. In the following, we will design an accelerated algorithm that utilizes this property. Our method is called PAME (PAM with Extrapolation) and it incorporates the extrapolation technique to the gradient step for updating λ, and f and g are still updated using exact maximizations. We note that currently we are not able to prove a better complexity for PAME. Our iteration complexity result in Theorem 22 is in the same order as that of PAM, but numerically we have observed great improvement of PAME over PAM. It is an interesting future topic to study other accelerations to PAM that can provably achieve improved complexity.
A typical iteration of our PAME algorithm is given below:
f t+1 = f t + η log
( a
r ( ∑
k ζ k(f t, gt, λt))
) , (26a)
gt+1 = gt + η log
( b
c ( ∑
k ζ k(f t+1, gt, λt))
) , (26b)
yt+1 = Proj∆N ( λt + (1− θ)(λt − λt−1) ) , (26c)
λt+1 = Proj∆N ( yt+1 + τ∇λF (f t+1, gt+1, yt+1) ) . (26d)
Here θ ∈ (0, 1) is a given parameter for the extrapolation step. We see that steps (26a)-(26b) are the same as (13)-(14) and they are solutions to the exact maximizations (12a)-(12b). Steps (26c)-(26c) give extrapolation to the gradient step for λ, similar to Nesterov’s accelerated gradient method. Note that PAME (26) solves the dual entropy-regularized EOT problem (11). We use the same rounding procedure in Section 2.1 to generate a primal solution to the original EOT problem (1). The complete PAME algorithm is described in Algorithm 3. Note that the algorithm is terminated when the following criteria are met:
∥ct−1 − b∥1 ≤ ϵ/(6(6c∞ − ηι), (27a) ∥λt−1 − λt−2∥2 ≤ ηϵ/(60(1− θ)c2∞), (27b)
∥λt − yt∥2 ≤ ηϵ/(42c2∞), (27c) F̃ (f t, gt−1, λt−1) ≤ ϵ/6. (27d)
1In Lemma 8 we proved that ∇λF is Lipschitz continuous. The Lipschitz continuity of ∇fF and ∇gF can be proved similarly.
Algorithm 3 Projected Alternating Maximization with Extrapolation Algorithm
1: Input: Cost matrices {Ck}1≤k≤N , accuracy ϵ, θ ∈ (0, 1). 2: Initialization: f0 = g0 = [1, ..., 1]⊤, λ0 = [1/N, ..., 1/N ]⊤ ∈ ∆N . 3: Choose parameters as
η = min
{ ϵ
3(log(n2N) + 1) , c∞
} , τ = η
2c2∞ . (28)
4: while (27) is not met do 5: Compute f t+1 by (26a) 6: Compute gt+1 by (26b) 7: Compute yt+1 by (26c) 8: Compute λt+1 by (26d) 9: t← t+ 1
10: end while 11: Assume the stop condition (27) is satisfied at the T -th iteration. Compute (ak, bk)k∈[N ] =
Margins(π(fT , gT−1, λT−1), a, b) as in Section 2.1. 12: Output: (π̂, λ̂) where π̂k = Round(πk(fT , gT−1, λT−1), ak, bk), ∀k ∈ [N ], λ̂ = λT−1.

heading: 4.1 CONVERGENCE ANALYSIS OF PAME ALGORITHM
text: In this section, we analyze the iteration complexity of PAME (Algorithm 3) for obtaining an ϵoptimal solution to the original EOT problem (1). The proof for PAME is different from that of PAM, and here we need to analyze the behavior of the following Hamiltonian, inspired by Jin et al. Jin et al. (2018).
E(f, g, λ1, λ2) = F (f, g, λ1)− 1 2τ ∥λ1 − λ2∥22. (29)
Theorem 5 Define ϵ′ = ϵ/(6c∞ − ηι), and set T to be
T = 8 + 48
η √ γ1ϵ′ +
( 3600(1− θ)2 + 882 ) c2∞s
ηϵ2 +
48
ηγ1ϵ = O
( c2∞ϵ −2) , (30) where γ1 = min { 1
(2c∞−ηι)2 , 2(2θ−θ 2) (7−5θ)2c2∞ , 149c2∞ } and we know γ1 = O(c−2∞ ). At least one of the
iterations in Algorithm 3, after rounding, is an ϵ-saddle point of the EOT problem (2).
Proof. See Appendix C. □
Remark 6 We are not able to analytically prove that PAME has an improved complexity bound at this moment yet. The APGA proposed in Scetbon et al. (2021) in fact has better complexity than PAM and PAME. However, as demonstrated in Scetbon et al. (2021) and in our numerical experiments (Sections 5 and D), APGA performs worse than PAM. We believe the reason is that APGA takes gradient step for the variables f and g, while PAM exactly minimizes the subproblems corresponding to these two variables. It is the exact minimization step that led to the improvement. Developing a provably better algorithm is definitely important and interesting, and we will work on it in the future.

heading: 5 NUMERICAL EXPERIMENTS
text: We compare the performance of PAME with PAM and APGA (25) Scetbon et al. (2021) on a synthetic dataset: the Gaussian distributions. We also conduct numerical comparison on another synthetic dataset: the fragmented hypercube dataset, and the results will be given in the Appendix.
Gaussian Distribution: Consider the case when two sets of discrete support {xi}i∈[n], {yj}j∈[n] are independently sampled from Gaussian distributions
N ((
1 1
) , ( 10 1 1 10 )) and N (( 2 2 ) , ( 1 −0.2 −0.2 1 )) (31)
respectively. The base cost matrix Cbase is computed by Cbasei,j = ∥xi − yj∥22. Assume we have N agents. The cost matrix of each agent can be obtained by adding Gaussian noise sampled from N (0, 10) to each element of the base cost. For instance, for the k-th agent with a cost matrix Ck, we have Cki,j = |Cbasei,j +N (0, 10)|.
We then set a = b = [1/n, ..., 1/n] for all experiments. For all algorithms, we set τ = 5ηc2∞ and we set θ = 0.1 for the PAME algorithm. We consider the EOT error as a measure of optimality. The EOT error at iteration t is defined by
Error = |ℓ(π(f t, gt, λt), λt)− ℓ∗|, (32)
where ℓ∗ is the approximated optimal value of EOT (2) obtained by running the PAM algorithm for 20000 iterations. Figures 1 plots the EOT error against the execution time for Gaussian distributions. We run each algorithm for 2000 iterations for different parameter settings. In all cases, the PAME and PAM perform significantly better than APGA, and PAME also shows significant improvement over PAM.

heading: 6 CONCLUSION
text: In this paper, we provided the first convergence analysis of the PAM algorithm for solving the EOT problem. Specifically, we have shown that it takes at most O(ϵ−2) iterations for the PAM algorithm to find an ϵ-saddle point. We proposed a PAME algorithm which incorporates the extrapolation technique to PAM. The PAME shows significant numerical improvement over PAM. Results in this paper might shed lights on designing new BCD type algorithms.

heading: B PROOF OF THEOREM 16

heading: B.1 TECHNICAL PREPARATIONS
text: We first give the partial gradients of F . [∇fF (f, g, λ)]i = ai − ∑
k,j exp((fi + gj − λkCkij)/η)∑ k ∥ζk(f, g, λ))∥1 = ai −
[ r (∑ k πk(f, g, λ) )] i ,
(33a) [∇gF (f, g, λ)]j = bj − ∑
k,i exp((fi + gj − λkCkij)/η)∑ k ∥ζk(f, g, λ))∥1 = bj −
[ c (∑ k πk(f, g, λ) )] j ,
(33b) [∇λF (f, g, λ)]k = ∑ i,j C k ij exp((fi + gj − λkCkij)/η)∑
k ∥ζk(f, g, λ))∥1 = ⟨πk(f, g, λ), Ck⟩. (33c)
Since (13) and (14) renormalize the row sum and column sum of ∑
k ζ k(f, g, λ) to be a and b, we
immediately have N∑
k=1
∥ζk(f t+1, gt, λt)∥1 = 1, N∑
k=1
∥ζk(f t+1, gt+1, λt)∥1 = 1,∀t, (34)
which, combined with (8), yields
πk(f t+1, gt, λt) = ζk(f t+1, gt, λt), πk(f t+1, gt+1, λt) = ζk(f t+1, gt+1, λt),∀t. (35)
The following lemma gives an error bound for Algorithm 2 (see Altschuler et al. (2017)).
Lemma 7 (Rounding Error) Let a, b ∈ Rn+ with ∑n i=1 ai = ∑n j=1 bj = q, π ∈ R n×n + , and π̂ = Round(π, a, b). The following inequality holds:
∥π̂ − π∥1 ≤ 2(∥r(π)− a∥1 + ∥c(π)− b∥1).
Proof. The proof is a slight modification from (Altschuler et al., 2017, Lemma 7). Note that Lines 2-5 in Algorithm 2 renormalize the row sum and column sum that are larger than the corresponding ai and bj . It is easy to verify that π̂, π′′, erra and errb are nonnegative with ∥erra∥1 = ∥errb∥1 = q − ∥π′′∥1 and
r(π̂) = r(π′′) + r(erraerr ⊤ b /∥erra∥1) = r(π′′) + erra = a ,
and likewise c(π̂) = b. Denote ∆ = ∥π∥1 − ∥π′′∥1. Since we remove mass from a row of π when ri(π) ≥ ai, and from a column when cj(π′) ≥ bj , we have
∆ = n∑ i=1 (ri(π)− ai)+ + n∑ j=1 (cj(π ′)− bj)+ .
Firstly, a simple calculation shows n∑
i=1
(ri(π)− ai)+ = 1
2 [∥r(π)− a∥1 + ∥π∥1 − q] .
Secondly, the fact that the vector c(π) is entrywise larger than c(π′) leads to n∑
j=1
(cj(π ′)− bj)+ ≤ n∑ j=1 (cj(π)− bj)+ ≤ ∥c(π)− b∥1.
Therefore we conclude
∥π̂ − π∥1 ≤ ∆+ ∥erraerr⊤b ∥1/∥erra∥1 = ∆+ q − ∥π′′∥1 = 2∆+ q − ∥π∥1 ≤ ∥r(π)− a∥1 + 2∥c(π)− b∥1 ≤ 2 [ ∥r(π)− a∥1 + ∥c(π)− b∥1 ] .
□
The following lemma shows that∇λF is Lipschitz continuous.
Lemma 8 For any f, g ∈ Rn and λ1, λ2 ∈ ∆N , the following inequality holds
∥∇λF (f, g, λ1)−∇λF (f, g, λ2)∥2 ≤ c2∞∥λ1 − λ2∥2/η, (36)
which immediately implies
F (f, g, λ1) ≥ F (f, g, λ2) + ⟨∇λF (f, g, λ2), λ1 − λ2⟩ − c2∞ 2η ∥λ1 − λ2∥22. (37)
Proof. The proof essentially follows Scetbon et al. (2021). It is easy to verify that the (q, k)-th entry of the Hessian of F (f, g, λ) with respect to λ is
∂2F
∂λq∂λk =
1
ην2 [σq,1(λ)σk,1(λ)− ν(σk,2(λ)11k=q)]
where 11k=q = 1 iff k = q and 0 otherwise, for all k ∈ {1, ..., N} and p ≥ 1
σk,p(λ) = ∑ i,j (Cki,j) p exp
( fi + gj − λkCki,j
η
) ,
ν = N∑ k=1 ∑ i,j exp
( fi + gj − λkCki,j
η
) .
Let v ∈ RN satisfying ∥v∥2 = 1, and by denoting∇2λF the Hessian of F with respect to λ for fixed f, g, we obtain
v⊤∇2λFv = 1
ην2 ( N∑ k=1 vkσk,1(λ) )2 − ν N∑ k=1 v2kσk,2  ≤ 1 ην2 ( N∑
k=1
vkσk,1(λ)
)2 − 1
ην2  N∑ k=1 |vk| √√√√∑ i,j exp ( fi + gj − λkCki,j η )√√√√∑ i,j (Cki,j) 2 exp ( fi + gj − λkCki,j η )2
≤ 1 ην2 ( N∑ k=1 vkσk,1(λ) )2 −  N∑ k=1 |vk| ∑ i,j |Cki,j | exp ( fi + gj − λkCki,j η )2  ≤ 0,
where the last three inequalities come from Cauchy Schwartz inequality. Moreover we have
vT∇2λFv = 1
ην2 ( N∑ k=1 vkσk,1(λ) )2 − ν N∑ k=1 v2kσk,2  ≥ −∑Nk=1 v2kσk,2 ην ≥ −c 2 ∞ η ,
which completes the proof. □
The next lemma gives a bound for g.
Lemma 9 Let (f t, gt, λt) be the sequence generated by Algorithm 1. For any t ≥ 0, it holds that max
j gtj −min j gtj ≤ c∞ − ηι, (38a)
max j g∗j −min j g∗j ≤ c∞ − ηι. (38b)
Proof. We prove (38a) first. When t = 0, (38a) holds because of the initialization g0. When t ≥ 1, from (15) we have
N∑ k=1 e−λ t−1 k C k ij/η ≥ N∑ k=1 e−λ t−1 k ∥C k∥∞/η ≥ N∑ k=1 e−∥C k∥∞/η ≥ Ne−c∞/η, (39)
where the second inequality is due to 0 ≤ λt−1k ≤ 1. Combining (39) and (34) we get
eg t j/η ·Ne−c∞/η⟨1, ef t/η⟩ ≤ ∑ i eg t j/η
( N∑
k=1
e−λ t−1 k C k ij/η ) ef t i /η = bj ≤ 1,
which leads to max
j gtj ≤ c∞ − η log(N⟨1, ef t/η⟩). (40)
Moreover, note that e−λ t−1 k C k ij/η ≤ 1, therefore 1N ∑N k=1 e −λt−1k C k ij/η ≤ 1. This fact leads to:
eg t j/η · ⟨1, ef t/η⟩ ≥ ∑ i eg t j/η
( 1
N N∑ k=1 e−λ t−1 k C k ij/η
) ef t i /η = 1
N bj ,
which gives min j gtj ≥ ηι− η log(N⟨1, ef t/η⟩). (41) Combining (40) with (41) yields (38a). The bound for g∗ (38b) can be obtained similarly, by noting that π∗ ∈ ΠNa,b. We omit the details for brevity. □
Lemma 10 Let {f t, gt, λt} be generated by PAM (Algorithm 1). The following equality holds. N∑ k ∥πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt)∥1 = ∥∥ct − b∥∥ 1 ,∀t.
Proof. By (35), we have N∑ k ∥πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt)∥1
= N∑ k ∑ i,j |e(f t+1 i +g t+1 j −λ t kC k i,j)/η − e(f t+1 i +g t j−λ t kC k i,j)/η|
= N∑ k ∑ i,j [πk(f t+1, gt, λt)]i,j ∣∣bj/ctj − 1∣∣ =∑ j ctj |bj/ctj − 1| = ∥∥ct − b∥∥ 1 .
□

heading: B.2 KEY LEMMAS
text: In this subsection, we provide a few useful lemmas that will lead to our main theorem on the iteration complexity of PAM (Algorithm 1). These lemmas yield the following results: the function F is monotonically increasing (Lemmas 11), the suboptimality of the dual problem can be upper bounded (Lemma 12-14), and the PAM returns an ϵ-optimal solution under conditions (55) (Lemma 15). In Theorem 16 we will show that these conditions can indeed be satisfied.
Lemma 11 [Increase of F ] Let {f t, gt, λt} be generated by PAM (Algorithm 1). The following inequalities hold:
F (f t+1, gt, λt)− F (f t, gt, λt) ≥ 0 (42a)
F (f t+1, gt+1, λt)− F (f t+1, gt, λt) ≥ η 2 ∥∥ct − b∥∥2 1
(42b)
F (f t+1, gt+1, λt+1)− F (f t+1, gt+1, λt) ≥ c2∞∥λt+1 − λt∥2/(2η). (42c)

heading: Proof.
text: First, (42a) is a direct consequence of (12a).
Next, we prove (42b). We have F (f t+1, gt+1, λt)− F (f t+1, gt, λt)
= ⟨gt+1 − gt, b⟩ − η log
( N∑
k=1
∥ζk(f t+1, gt+1, λt)∥1 ) + η log ( N∑
k=1
∥ζk(f t+1, gt, λt)∥1
)
= ⟨gt+1 − gt, b⟩ = η n∑
j=1
bj log(bj/c t j) = ηK(b||ct) ≥
η 2 ∥ct − b∥21,
where K(x||y) denotes the KL divergence of x and y, the second equality is due to (34), the third equality is due to (14), and the last inequality follows the Pinsker’s inequality.
Finally, we prove (42c). From the optimality condition of (12c), we know that there exists h(λt+1) ∈ ∂I∆N (λt+1) (43)
such that ∇λF (f t+1, gt+1, λt)− 1
τ (λt+1 − λt)− h(λt+1) = 0. (44)
From (37) we have
F (f t+1, gt+1, λt+1)− F (f t+1, gt+1, λt) ≥ ⟨∇λF (f t+1, gt+1, λt), λt+1 − λt⟩ − c2∞ 2η ∥λt+1 − λt∥2
= ⟨1 τ (λt+1 − λt) + h(λt+1), λt+1 − λt⟩ − c 2 ∞ 2η ∥λt+1 − λt∥2
≥ ⟨1 τ (λt+1 − λt), λt+1 − λt⟩ − c 2 ∞ 2η ∥λt+1 − λt∥2
= c2∞∥λt+1 − λt∥2/(2η),
where the first equality is due to (44), the second inequality is due to (43), and the last equality is due to the definition of τ in (20).
□
Before we bound the suboptimality gap, we need the following lemma.
Lemma 12 Let {f t, gt, λt} be generated by PAM (Algorithm 1). For any λ ∈ ∆N , the following inequality holds:〈
λ− λt,∇λF (f t+1, gt, λt) 〉 ≤ 3c2∞∥λt+1 − λt∥2/η + c∞ ∥∥ct − b∥∥ 1 . (45)
Proof. The optimality condition of (12c) is given by:
⟨λ− λt+1, 1 τ (λt+1 − λt)−∇λF (f t+1, gt+1, λt)⟩ ≥ 0, ∀λ ∈ ∆N , (46)
which implies that
⟨λt+1 − λ,−∇λF (f t+1, gt+1, λt)⟩
≤ ⟨λ− λt+1, 1 τ (λt+1 − λt)⟩ ≤ 1 τ ∥λ− λt+1∥2∥λt+1 − λt∥2 ≤ 2c2∞∥λt+1 − λt∥2/η,
(47)
where the last inequality is due to the fact that the diameter of ∆N is bounded by √ 2 ≤ 2. Moreover, we have ⟨λt − λ,∇λF (f t+1, gt+1, λt)−∇λF (f t+1, gt, λt)⟩
= N∑ k (λtk − λk) · ⟨πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt), Ck⟩
≤ N∑ k ∥πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt)∥1∥Ck∥∞
≤ c∞∥ct − b∥1,
(48)
where the equality is due to (33c), and the last inequality is due to Lemma 10. Finally, we have〈 λt − λ,−∇λF (f t+1, gt, λt) 〉 = 〈 λt − λt+1,−∇λF (f t+1, gt+1, λt) 〉 + 〈 λt+1 − λ,−∇λF (f t+1, gt+1, λt) 〉 +〈
λt − λ,∇λF (f t+1, gt+1, λt)−∇λF (f t+1, gt, λt) 〉
≤∥λt − λt+1∥2 · ∥∇λF (f t+1, gt+1, λt)∥2 + 2c2∞∥λt+1 − λt∥2/η + c∞∥ct − b∥1,
(49)
where the first inequality is due to (47) and (48). From (33c) we have ∥∇λF (f t+1, gt+1, λt)∥2 ≤ c∞, which, combined with (49) and the fact that η ≤ c∞, yields the desired result. □
The suboptimality of (11) is defined as: F̃ (f, g, λ) = F (f∗, g∗, λ∗) − F (f, g, λ). Note that F̃ (f, g, λ) ≥ 0,∀f, g, λ ∈ ∆N .
Lemma 13 Let (f t, gt, λt) be generated by PAM (Algorithm 1). The following inequality holds:
F̃ (f t+1, gt, λt) ≤ (2c∞ − ηι)∥ct − b∥1 + 3c2∞∥λt+1 − λt∥2/η.
Proof. Denote ut = (maxj gtj +minj g t j)/2, u ∗ = (maxj g ∗ j +minj g ∗ j )/2. From (35) we get
⟨1, ct − b⟩ = n∑
i=1
ai − n∑
j=1
bj = 0,
which further implies〈 gt − g∗, ct − b 〉 = 〈 (gt − ut1)− (g∗ − u∗1), ct − b 〉 ≤ ( ∥gt − ut1∥∞ + ∥g∗ − u∗1∥∞ ) ∥∥ct − b∥∥ 1 ≤ (c∞ − ηι) ∥∥ct − b∥∥ 1 , (50)
where the last inequality is due to Lemma 9. Now we set λ = λ∗ in (45), and we obtain〈 λt − λ∗,−∇λF (f t+1, gt, λt) 〉 ≤ 3c2∞∥λt+1 − λt∥2/η + c∞ ∥∥ct − b∥∥ 1 . (51)
Since F (f, g, λ) is a concave function, we have
F (f∗, g∗, λ∗) ≤ F (f t+1, gt, λt) + ⟨∇F (f t+1, gt, λt), (f∗, g∗, λ∗)− (f t+1, gt, λt)⟩,
which, combining with (33) yields
F̃ (f t+1, gt, λt) = F (f∗, g∗, λ∗)− F (f t+1, gt, λt) ≤ ⟨f t+1 − f∗, r( ∑N
k=1 π k(f t+1, gt, λt))− a⟩+ ⟨gt − g∗, ct − b⟩ + 〈 λt − λ∗,−∇λF (f t+1, gt, λt) 〉 ≤ (2c∞ − ηι)∥ct − b∥1 + 3c2∞∥λt+1 − λt∥2/η,
where the last inequality follows from (15), (35), (50) and (51). □
The next lemma shows that the suboptimality gap F̃ (f, g, λ) can be bounded by O(1/t).
Lemma 14 Let (f t, gt, λt) be generated by PAM (Algorithm 1). The following inequality holds:
F̃ (f t+1, gt+1, λt+1) ≤ 4/(ηγ0) t+ 1 + 4/(ηγ0F̃ (f0, g0, λ0)) ,
where γ0 = min {
1 (2c∞−ηι)2 , 19c2∞
} is a constant.
Proof. Combining (42b) and (42c), we have
F (f t+1, gt+1, λt+1)− F (f t+1, gt, λt) ≥ η 2 ∥∥ct − b∥∥2 1 + c2∞ ∥∥λt+1 − λt∥∥2 2 /(2η). (52)
Therefore, we have
F̃ (f t+1, gt+1, λt+1)− F̃ (f t+1, gt, λt)
≤− η 2 ∥∥ct − b∥∥2 1 − c2∞ ∥∥λt+1 − λt∥∥2 2 /(2η) ≤− η 2 γ0 · ( ((2c∞ − ηι)∥ct − b∥1)2 + (3c2∞∥λt+1 − λt∥2/η)2
) ≤− η
4 γ0 ( (2c∞ − ηι)∥ct − b∥1 + 3c2∞∥λt+1 − λt∥2/η )2 ≤− η
4 γ0F̃ (f
t+1, gt, λt)2,
(53)
where the last inequality is from Lemma 13. Dividing both sides of (53) by F̃ (f t+1, gt+1, λt+1) · F̃ (f t+1, gt, λt), we have
1 F̃ (f t+1, gt+1, λt+1) ≥ 1 F̃ (f t+1, gt, λt) + η 4 γ0 ·
F̃ (f t+1, gt, λt)
F̃ (f t+1, gt+1, λt+1)
≥ 1 F̃ (f t+1, gt, λt) + η 4 γ0 ≥
1
F̃ (f t, gt, λt) +
η 4 γ0,
(54)
where the second inequality is due to (53) and the last inequality is from (42a). Summing (54) from 0 to t leads to
1 F̃ (f t+1, gt+1, λt+1) ≥ 1 F̃ (f0, g0, λ0) + η(t+ 1) 4 γ0,
which implies the desired result. □
The next lemma gives sufficient conditions for the PAM algorithm to return an ϵ-optimal solution to the original EOT problem (2).
Lemma 15 Assume PAM terminates at the T -iteration, i.e.,
∥cT−1 − b∥1 ≤ ϵ/(6(6c∞ − ηι)), (55a)∥∥λT − λT−1∥∥ 2 ≤ ηϵ/(18c2∞), (55b) F̃ (fT , gT−1, λT−1) ≤ ϵ/6. (55c)
Then the output (π̂, λ̂) of PAM (Algorithm 1), i.e., π̂k = Round(πk(fT , gT−1, λT−1), ak, bk), ∀k ∈ [N ], λ̂ = λT−1, is an ϵ-optimal solution of the original EOT problem (2).
Proof. According to Definition 2, it is sufficient to show that the output (π̂, λ̂) ∈ ΠNa,b×∆N satisfies the following two inequalities:
max λ∈∆N ℓ (π̂, λ)− ℓ(π̂, λ̂) ≤ ϵ 2 , (56a)
ℓ(π̂, λ̂)− min π∈ΠNa,b ℓ(π, λ̂) ≤ ϵ 2 . (56b)
We prove (56a) first. For ease of presentation, we denote π̃ = π(fT , gT−1, λT−1), π∗ = π(f∗, g∗, λ∗). Note that π̂k = Round(π̃k, ak, bk),∀k ∈ [N ]. We also denote
λ̄(π) := argmax λ∈∆N
{ ℓ(π, λ) =
N∑ k=1 λk⟨πk, Ck⟩
} . (57)
Note that the term on the left hand side of (56a) can be rewritten as
ℓ ( π̂, λ̄(π̂) ) − ℓ ( π̂, λ̂ ) =(ℓ(π̂, λ̄(π̂))− ℓ(π̃, λ̄(π̃)))︸ ︷︷ ︸
(I)
+([ℓ(π̃, λ̄(π̃))− ηH(π̃)]− [ℓ(π∗, λ∗)− ηH(π∗)])︸ ︷︷ ︸ (II)
+ ([ℓ(π∗, λ∗)− ηH(π∗)]− [ℓ(π̃, λ̂)− ηH(π̃)])︸ ︷︷ ︸ (III) +(ℓ(π̃, λ̂)− ℓ(π̂, λ̂))︸ ︷︷ ︸ (IV ) .
(58)
We now provide upper bounds for these four terms. Denote
k̂∗ = argmax k∈[N ] ⟨π̂k, Ck⟩, k̃∗ = argmax k∈[N ] ⟨π̃k, Ck⟩. (59)
Since (1) and (2) are equivalent, we have the following for the term (I):
(I) = ∑ k [λ̄(π̂)]k⟨π̂k, Ck⟩ − ∑ k [λ̄(π̃)]k⟨π̃k, Ck⟩ = ⟨π̂k̂ ∗ , C k̂ ∗ ⟩ − ⟨π̃k̃ ∗ , C k̃ ∗ ⟩
≤⟨π̂k̂ ∗ , C k̂ ∗ ⟩ − ⟨π̃k̂ ∗ , C k̂ ∗ ⟩ ≤ ∥π̂k̂ ∗ − π̃k̂ ∗ ∥1∥Ck∥∞ ≤ c∞ ∑ k ∥π̂k − π̃k∥1
≤2c∞ ∑ k (∥r(π̃k)− ak∥1 + ∥c(π̃k)− bk∥1) = 2c∞∥cT−1 − b∥1,
(60)
where the first inequality follows from the definition of k̃∗ in (59), the fourth inequality is from Lemma 7, and the last equality follows from (15) and (17).
For the term (II), recall that H(π) = − ∑
k,i,j π k i,j(log π k i,j − 1) and π̃k = exp ( fT 1⊤+1(gT−1)⊤−λT−1k C k
η
) due to (35), and define uT−1 = maxj g T−1 j +minj g T−1 j
2 . We
have (II) = ∑ k λ̄(π̃)k⟨π̃k, Ck⟩+ η ∑ k,i,j π̃ki,j ( fTi + g T−1 j − λ T−1 k C k ij η − 1 ) − F ∗
= ∑ k (λ̄(π̃)k − λ̂k)⟨π̃k, Ck⟩+ ∑ k,i,j π̃ki,j ( fTi + g T−1 j − η ) − F ∗
= 〈 λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1) 〉 + 〈 fT , a 〉 + 〈 gT−1, cT−1 〉 − η ∑ i,j,k π̃ k i,j − F ∗
= 〈 λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1) 〉 + 〈 fT , a 〉 + 〈 gT−1, cT−1 〉 − log (∑ k ∥π̃k∥1 ) − η − F ∗
= 〈 λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1) 〉 + ⟨gT−1, cT−1 − b⟩+ F (fT , gT−1, λT−1)− F ∗
≤ 〈 λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1) 〉 + ⟨gT−1, cT−1 − b⟩
≤c∞∥cT−1 − b∥1 + 3c2∞∥λT − λT−1∥2/η + ⟨gT−1 − uT−11, cT−1 − b⟩ ≤c∞∥cT−1 − b∥1 + 3c2∞∥λT − λT−1∥2/η + ∥gT−1 − uT−11∥∞∥cT−1 − b∥1 ≤(3c∞/2− ηι/2)∥cT−1 − b∥1 + 3c2∞∥λT − λT−1∥2/η,
(61) where the third equality uses (35), (33c) and (15), the second inequality follows from Lemma 12 by setting λ = λ̄(π̃) and t = T − 1, and the last inequality uses Lemma 9. For the term (III), we have
(III) ≤ ∣∣∣∣∣∣ ∑ k λ̂k⟨π̃k, Ck⟩+ η ∑ i,j π̃ki,j ( fTi + g T−1 j − λ T−1 k C k η − 1 ) − F ∗ ∣∣∣∣∣∣ =|⟨gT−1, cT−1 − b⟩+ F (fT , gT−1, λT−1)− F ∗| ≤(c∞/2− ηι/2)∥cT−1 − b∥1 + |F (fT , gT−1, λT−1)− F ∗|,
(62)
where the last inequality follows from Lemma 9.
Finally, for the term (IV), we have (IV ) = ∑ k ⟨π̃k − π̂k, λ̂kCk⟩ ≤ ∑ k ∥π̃k − π̂k∥1∥Ck∥∞
≤2c∞ ∑ k (∥r(π̃k)− ak∥1 + ∥c(π̃k)− bk∥1) = 2c∞∥cT−1 − b∥1, (63)
where the first inequality uses |λ̂k| ≤ 1, the second inequality uses Lemma 7 and (17). Plugging (60) - (63) into (58), and using (55), we obtain (56a).
Now we prove (56b). For ease of presentation, we denote
π̄(λ) := argmin π∈ΠNa,b ℓ(π, λ). (64)
We also denote b̃ = cT−1 = ∑ k c ( π̃k ) and π′k = Round(π̄(λ̂)k, ãk, b̃k), where
(ãk, b̃k)k∈[N ] := Margins(π̄(λ̂), a, b̃),
as defined in (16). From (18) we know that∑ k ∥∥∥c((π̄(λ̂))k)− b̃k∥∥∥ 1 = ∥∥∥∥∥∑ k c ( (π̄(λ̂))k ) − ∑ k b̃k ∥∥∥∥∥ 1 = ∥b− b̃∥1 = ∥b− cT−1∥1, (65)
where the second equality is due to π̄(λ̂) ∈ ΠNa,b and thus c( ∑
k(π̄(λ̂)) k) = b, and the fact that∑
k b̃ k = b̃ due to Property (ii) of the Margins procedure in Section 2.1. By the Sinkhorn’s theorem
Sinkhorn (1967), π̃ is the unique optimal solution of minπ∈ΠN a,b̃ ℓη(π, λ̂). Therefore∑ k λ̂k⟨π̃k, Ck⟩ − ηH(π̃) ≤ ∑ k λ̂k⟨π′k, Ck⟩ − ηH(π′). (66)
Now, note that the left hand side of (56b) can be arranged into three parts:
ℓ(π̂, λ̂)− ℓ(π̄(λ̂), λ̂)
= (∑ k λ̂k⟨π̂k, Ck⟩ − ∑ k λ̂k⟨π̃k, Ck⟩ ) ︸ ︷︷ ︸
(V )
+ (∑ k λ̂k⟨π̃k, Ck⟩ − ∑ k λ̂k⟨π′k, Ck⟩ ) ︸ ︷︷ ︸
(V I)
+ (∑ k λ̂k⟨π′k, Ck⟩ − ∑ k λ̂k⟨(π̄(λ̂))k, Ck⟩ ) ︸ ︷︷ ︸
(V II)
.
(67)
We now upper bound these three terms. First note that the term (V) is the same as the term (IV) and thus has the same upper bound in (63). Since 0 ≤ H(π) ≤ log(n2N) + 1, from (66) we have that
(V I) = ∑ k λ̂k⟨π̃k, Ck⟩ − ∑ k λ̂k⟨π′k, Ck⟩ ≤ η ∣∣H(π̃)−H(π′)∣∣ ≤ 1 3 ϵ, (68)
where the last step uses the definition of η in (20).
For the term (VII), we have (V II) = ∑ k λ̂k⟨π′k, Ck⟩ − ∑ k λ̂k⟨(π̄(λ̂))k, Ck⟩ ≤ ∑ k ∥π′k − (π̄(λ̂))k∥1∥Ck∥∞
≤2c∞ ∑ k (∥r((π̄(λ̂))k)− ãk∥1 + ∥c((π̄(λ̂))k)− b̃k∥1) = 2c∞∥cT−1 − b∥1, (69)
where the second inequality follows from Lemma 7, the second equality uses (65) and the fact that r((π̄(λ̂))k) = ãk due to the property of the Margins procedure in (16).
Finally, plugging (63) (note (V)=(IV)), (68) and (69) into (67), and using (55a) and noting ι < 0, we obtain (56b). This completes the proof.
□

heading: B.3 MAIN RESULT
text: Theorem 16 Define ϵ′ = ϵ/(6c∞ − ηι), and set T as
T = 5 + 36
η √ γ0ϵ′ + 648c2∞ ηϵ + 28 ηγ0ϵ = O
( c2∞ϵ −2) , (70) where γ0 is defined in Lemma 14 and we know γ0 = O(c−2∞ ). The output pair of Algorithm 1 is an ϵ-optimal solution of the EOT problem (2).
Proof. According to Lemma 15, we only need to show that (55) holds after T iterations as defined in (70). To guarantee (55a) and (55b), we follow the ideas of Dvurechensky et al. Dvurechensky et al. (2018) and construct a switching process. We first reduce F̃ from F̃ (f0, g0, λ0) to a constant s by running t1 steps. In this process, Lemma 14 indicates
t1 ≤ 1 + 4 ηγ0s − 4 ηγ0F̃ (f0, g0, λ0) . (71)
Secondly, starting from s, we continue running the algorithm, and assume that there are t2 iterations in which (55a) fails. By (42b) we have
t2 ≤ 1 + 72s
ηϵ′2 .
Therefore, we know that the total iteration number that (55a) fails is upper bounded by
T1 = t1 + t2 ≤ 2 + 72s
ηϵ′2 +
4 ηγ0s − 4 ηγ0F̃ (f0, g0, λ0)
iterations. By choosing s = ϵ′/(6 √ γ0), we know that
T1 ≤
{ 2 + 12η√γ0ϵ′ + 24 η √ γ0ϵ′ − 4 ηγ0F̃ (f0,g0,λ0) ≤ 2 + 36η√γ0ϵ′ if F̃ (f 0, g0, λ0) ≥ ϵ ′ 6 √ γ0
2 + 12η√γ0ϵ′ + 24 η √ γ0ϵ′ − 4 ηγ0F̃ (f0,g0,λ0) ≤ 2 + 12η√γ0ϵ′ otherwise.
Therefore, we have T1 ≤ 2 + 36η√γ0ϵ′ . Similarly, starting from s, the number of iterations that (55b) fails can be bounded by
t3 ≤ 1 + 648sc2∞ ηϵ2 ,
where we apply (42b). By choosing s = ϵ, we know that the total iteration number that (55b) fails is upper bounded by
T2 = t1 + t3 ≤ 2 + 648c2∞ ηϵ + 4 ηγ0ϵ − 4 ηγ0F̃ (f0, g0, λ0)
iterations. Finally, by letting s = ϵ/6 in (71), we know that
F̃ (fT3−1, gT3−1, λT3−1) ≤ ϵ/6
after T3 = 1 + 24
ηγ0ϵ
iterations. From (42a) we know that after T3 iterations, we have
F̃ (fT3 , gT3−1, λT3−1) ≤ F̃ (fT3−1, gT3−1, λT3−1) ≤ ϵ/6,
i.e., (55c) holds. Combining the above discussions, we know that after T = T1 + T2 + T3 + 1 iterations, there must exist at least one iteration such that (55) holds, and thus the output of PAM is an ϵ-optimal solution to the original EOT problem (2).
□

heading: C PROOF OF THEOREM 22
text: The following simple fact is useful for our analysis later. ∥yt+1 − λt∥2 = ∥Proj∆N ( λt + (1− θ)(λt − λt−1) ) − Proj∆N ( λt ) ∥2 ≤ (1− θ)∥λt − λt−1∥2, (72) where the equality follows from the definition of yt+1 in (26c), and the inequality is due to the non-expansiveness of the projection operator.
The following lemma shows that the Hamiltonian E(f t, gt, λt, λt−1) is monotonically increasing when updating λ in Algorithm (3).
Lemma 17 [Sufficient increase in λ] Let {f t, gt, yt, λt} be generated by PAME (Algorithm 3). The following inequality holds:
E(f t+1, gt+1, λt+1, λt)− E(f t+1, gt+1, λt, λt−1) ≥ 2θ − θ 2
2τ ∥λt − λt−1∥22 +
1
4τ ∥λt+1 − yt+1∥22.
(73) Note that since θ ∈ (0, 1), the right hand side of (73) is always nonnegative.
Proof. From the optimality condition of (26d) we know that, there exists h(λt+1) ∈ ∂I∆N (λt+1) such that
∇λF (f t+1, gt+1, yt+1)− 1
τ (λt+1 − yt+1)− h(λt+1) = 0. (74)
By the convexity of the indicator function I∆N (λt+1), we have ⟨yt+1 − λt+1, h(λt+1)⟩ ≤ 0, ⟨λt − λt+1, h(λt+1)⟩ ≤ 0. (75)
Moreover, we have the following inequality:
∥λt+1 − λt∥22 = ∥λt+1 − yt+1 + yt+1 − λt∥22 = ∥yt+1 − λt∥22 + 2⟨λt+1 − yt+1, yt+1 − λt⟩+ ∥λt+1 − yt+1∥22 ≤ (1− θ)2∥λt − λt−1∥22 + 2⟨λt+1 − yt+1, yt+1 − λt⟩+ ∥λt+1 − yt+1∥22,
(76)
where the inequality is from (72).
We then have the following inequality:
F (f t+1, gt+1, λt)− F (f t+1, gt+1, λt+1) ≤ ( F (f t+1, gt+1, yt+1) + ⟨∇λF (f t+1, gt+1, yt+1), λt − yt+1⟩ ) −(
F (f t+1, gt+1, yt+1) + ⟨∇λF (f t+1, gt+1, yt+1), λt+1 − yt+1⟩ − c2∞∥λt+1 − yt+1∥22/(2η) )
= ⟨∇λF (f t+1, gt+1, yt+1), λt − λt+1⟩+ c2∞∥λt+1 − yt+1∥22/(2η) ≤ ⟨∇λF (f t+1, gt+1, yt+1)− h(λt+1), λt − λt+1⟩+ c2∞∥λt+1 − yt+1∥22/(2η)
= 1 τ ⟨λt+1 − yt+1, λt − λt+1⟩+ 1 4τ ∥λt+1 − yt+1∥2 = 1
τ ⟨λt+1 − yt+1, λt − yt+1 + yt+1 − λt+1⟩+ 1 4τ ∥λt+1 − yt+1∥2
= − 1 τ ⟨λt+1 − yt+1, yt+1 − λt⟩ − 3 4τ ∥λt+1 − yt+1∥2,
(77) where the first inequality is from the concavity of F with respect to λ and (37), the second inequality is due to (75), the second equality is due to (74). Combining (76) and (77) leads to
E(f t+1, gt+1, λt+1, λt) = F (f t+1, gt+1, λt+1)− 1 2τ ∥λt+1 − λt∥22
≥ F (f t+1, gt+1, λt) + 1 τ ⟨λt+1 − yt+1, yt+1 − λt⟩+ 3 4τ ∥λt+1 − yt+1∥2
− (1− θ) 2
2τ ∥λt − λt−1∥22 −
1 τ ⟨λt+1 − yt+1, yt+1 − λt⟩ − 1 2τ ∥λt+1 − yt+1∥22
= F (f t+1, gt+1, λt)− 1 2τ ∥λt − λt−1∥22 +
2θ − θ2
2τ ∥λt − λt−1∥22 +
1
4τ ∥λt+1 − yt+1∥2
= E(f t+1, gt+1, λt, λt−1) + 2θ − θ2
2τ ∥λt − λt−1∥22 +
1
4τ ∥λt+1 − yt+1∥2,
which completes the proof. □
Now we define the following function Ẽ, and later we will prove that Ẽ(f t, gt, λt, λt−1) can be upper bounded by O(1/t).
Ẽ(f, g, λ1, λ2) = F (f∗, g∗, λ∗)− E(f, g, λ1, λ2).
The next lemma is useful for obtaining the upper bound for Ẽ(f t, gt, λt, λt−1). Moreover, it is noted that Ẽ(f, g, λ1, λ2) ≥ 0,∀f, g, λ1, λ2, and Ẽ(f, g, λ, λ) = F̃ (f, g, λ),∀f, g, λ.
Lemma 18 Let {f t, gt, yt, λt} be generated by PAME (Algorithm 3). For any λ ∈ ∆N , the following inequality holds〈 λ− λt,∇λF (f t+1, gt, λt) 〉 ≤ c∞∥ct − b∥1 + 7c2∞∥λt+1 − yt+1∥2/η + 5(1− θ)c2∞∥λt − λt−1∥2/η.
(78)
Proof. From the optimality condition of (26d), we have the following inequality:〈 λ− λt+1, 1
τ (λt+1 − yt+1)−∇λF (f t+1, gt+1, yt+1)
〉 ≥ 0, ∀λ ∈ ∆N . (79)
The left hand side of (78) can be rearranged to three terms.
⟨λ− λt,∇λF (f t+1, gt, λt)⟩ = ⟨λt − λ,−∇λF (f t+1, gt+1, yt+1)⟩︸ ︷︷ ︸
(I)
+ ⟨λt − λ,∇λF (f t+1, gt+1, yt+1)−∇λF (f t+1, gt+1, λt)⟩︸ ︷︷ ︸ (II)
+ ⟨λt − λ,∇λF (f t+1, gt+1, λt)−∇λF (f t+1, gt, λt)⟩︸ ︷︷ ︸ (III) .
(80) We now bound these three terms one by one. To bound the term (I), we first note that from (33c) and (15), we have ∥∇λF (f t+1, gt+1, λt)∥2 ≤ c∞ ≤ c2∞/η, (81) where the second inequality is due to the definition of η (28). Now we can bound the term (I) as follows:
(I) = 〈 λt − λt+1,−∇λF (f t+1, gt+1, λt) 〉 + 〈 λt − λt+1,∇λF (f t+1, gt+1, λt)−∇λF (f t+1, gt+1, yt+1) 〉 + 〈 λt+1 − λ,−∇λF (f t+1, gt+1, yt+1)
〉 ≤ ∥∥λt − λt+1∥∥
2 · ∥∥∇λF (f t+1, gt+1, λt)∥∥2 + c2∞η ∥∥λt − λt+1∥∥2 · ∥∥λt − yt+1∥∥2
+ 1
τ
∥∥λt+1 − λ∥∥ 2 · ∥∥λt+1 − yt+1∥∥ 2
≤3c2∞∥λt − λt+1∥2/η + 4c2∞∥λt+1 − yt+1∥2/η, (82)
where the first inequality uses Lemma 8 and (79), the second inequality uses (81) and the facts that∥∥λt − yt+1∥∥ 2 ≤ 2 and ∥λt − λ∥2 ≤ 2.
For the term (II), Lemma 8 yields: (II) ≤2 ∥∥∇λF (f t+1, gt+1, yt+1)−∇λF (f t+1, gt+1, λt)∥∥2 ≤ 2c2∞ ∥∥yt+1 − λt∥∥2 /η. (83)
For the term (III), it can be bounded as:
(III) = N∑ k=1 (λtk − λk) · ⟨πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt), Ck⟩
≤ N∑
k=1
∥πk(f t+1, gt+1, λt)− πk(f t+1, gt, λt)∥1∥Ck∥∞ ≤ c∞∥ct − b∥1,
(84)
where the last inequality is due to Lemma (10). Plugging (82) - (84) into (80) and applying the triangle inequality, we obtain〈
λ− λt,∇λF (f t+1, gt, λt) 〉 ≤ c∞∥ct − b∥1 + 7c2∞∥λt+1 − yt+1∥2/η + 5c2∞∥yt+1 − λt∥2/η,
which immediately implies (78) by noting (72).
□
Lemma 19 Let (f t, gt, yt, λt) be generated by PAME (Algorithm 3). The following inequality holds:
Ẽ(f t+1, gt, λt, λt−1) ≤ (2c∞−ηι)∥ct−b∥1+7c2∞∥λt+1−yt+1∥2/η+(7−5θ)c2∞∥λt−λt−1∥2/η.
Proof. Since F (f, g, λ) is a concave function, we have F (f∗, g∗, λ∗) ≤ F (f t+1, gt, λt) + 〈 ∇F (f t+1, gt, λt), (f∗, g∗, λ∗)− (f t+1, gt, λt) 〉 ,
which implies that
F̃ (f t+1, gt, λt) ≤⟨gt − g∗, ct − b⟩+ ⟨λt − λ∗,−∇λF (f t+1, gt, λt)⟩ ≤(c∞ − ηι)∥ct − b∥1 + c∞∥ct − b∥1 + 7c2∞∥λt+1 − yt+1∥2/η + 5(1− θ)c2∞∥λt − λt−1∥2/η.
(85)
where in the first inequality we have used (35), and the second inequality follows from (50) and setting λ = λ∗ in (78). From (85) we immediately get
Ẽ(f t+1, gt, λt, λt−1) = F̃ (f t+1, gt, λt) + 1
2τ ∥λt − λt−1∥22
≤ c2∞∥λt − λt−1∥22/η + (2c∞ − ηι)∥ct − b∥1 + 7c2∞∥λt+1 − yt+1∥2/η + 5(1− θ)c2∞∥λt − λt−1∥2/η ≤ 2c2∞∥λt − λt−1∥2/η + (2c∞ − ηι)∥ct − b∥1 + 7c2∞∥λt+1 − yt+1∥2/η + 5(1− θ)c2∞∥λt − λt−1∥2/η,
where the second inequality is due to ∥λt − λt−1∥2 ≤ 2. This completes the proof. □
The following lemma bounds Ẽ(f t+1, gt+1, λt+1, λt) by O(1/t).
Lemma 20 Let {f t, gt, ytλt} be generated by PAME (Algorithm 3). The following inequality holds:
Ẽ(f t+1, gt+1, λt+1, λt) ≤ 6/(ηγ1) t+ 1 + 6/(ηγ1F̃ (f0, g0, λ0)) ,∀t ≥ 0,
where we assume λ−1 = λ0, and
γ1 = min
{ 1
(2c∞ − ηι)2 , 2(2θ − θ2) (7− 5θ)2c2∞ , 1 49c2∞
} (86)
is a constant.
Proof. Combining (42b) and Lemma 17, we have
E(f t+1, gt+1, λt+1, λt)− E(f t+1, gt, λt, λt−1) = ( E(f t+1, gt+1, λt+1, λt)− E(f t+1, gt+1, λt, λt−1) ) + ( F (f t+1, gt+1, λt)− F (f t+1, gt, λt) ) ≥ η
2 ∥ct − b∥21 +
2θ − θ2
2τ
∥∥λt − λt−1∥∥2 2 + 1
4τ
∥∥λt+1 − yt+1∥∥2 2 ,
which implies that
Ẽ(f t+1, gt+1, λt+1, λt)− Ẽ(f t+1, gt, λt, λt−1)
≤ − η 2 ∥ct − b∥21 − (2θ − θ2) c2∞ η ∥λt − λt−1∥22 − c2∞ 2η ∥λt+1 − yt+1∥22
≤ − η 2 γ1
[( (2c∞ − ηι) ∥∥ct − b∥∥ 1 )2 + ( (7− 5θ)c2∞∥λt − λt−1∥2/η )2 + ( 7c2∞∥λt+1 − yt+1∥2/η )2] ≤ − η
6 γ1 [ (2c∞ − ηι) ∥∥ct − b∥∥ 1 + (7− 5θ)c2∞∥λt − λt−1∥2/η + 7c2∞∥λt+1 − yt+1∥2/η ]2 ≤ − η
6 γ1Ẽ(f
t+1, gt, λt, λt−1)2,
(87) where the last inequality applies Lemma 19. We then divide both sides of (87) by Ẽ(f t+1, gt+1, λt+1, λt) · Ẽ(f t+1, gt, λt, λt−1), and we obtain
1 Ẽ(f t+1, gt+1, λt+1, λt) ≥ 1 Ẽ(f t+1, gt, λt, λt−1) + η 6 γ1 ·
Ẽ(f t+1, gt, λt, λt−1)
Ẽ(f t+1, gt+1, λt+1, λt)
≥ 1 Ẽ(f t+1, gt, λt, λt−1) + η 6 γ1 ≥
1
Ẽ(f t, gt, λt, λt−1) +
η 6 γ1,
(88)
where the second inequality holds because (87) implies that Ẽ(f t+1, gt, λt, λt−1) ≥ Ẽ(f t+1, gt+1, λt+1, λt), and the last inequality follows from (42a). Summing (88) from 0 to t leads to
1 Ẽ(f t+1, gt+1, λt+1, λt) ≥ 1 Ẽ(f0, g0, λ0, λ−1) + η(t+ 1) 6 γ1 =
1
F̃ (f0, g0, λ0) +
η(t+ 1)
6 γ1,
which immediately leads to the desired result. □
Similar to Lemma 15, the following lemma provides some sufficient conditions for the PAME algorithm to return an ϵ-optimal solution to the original EOT problem (2).
Lemma 21 Assume PAME terminates at the T -iteration, i.e.,
∥cT−1 − b∥1 ≤ ϵ/(6(6c∞ − ηι), (89a) ∥λT−1 − λT−2∥2 ≤ ηϵ/(60(1− θ)c2∞), (89b)
∥λT − yT ∥2 ≤ ηϵ/(42c2∞), (89c) F̃ (fT , gT−1, λT−1) ≤ ϵ/6. (89d)
Then the output (π̂, λ̂) of PAME (Algorithm 3), i.e., π̂k = Round(πk(fT , gT−1, λT−1), ak, bk), ∀k ∈ [N ], λ̂ = λT−1, is an ϵ-optimal solution of the original EOT problem (2).
Proof. The proof is essentially the same as that of Lemma 15. More specifically, we again need to show that the output of PAME (π̂, λ̂) satisfies (56). The proof of (56b) is exactly the same as the proof of Lemma 15. The proof of (56a) only requires to develop a new bound for〈
λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1) 〉
(90)
that is used in (61). Other parts are again exactly the same as the ones in Lemma 15. The new bound of (90) can be obtained by applying Lemma 18 with λ = λ̄(π̃) and t = T − 1, which yields
⟨λ̄(π̃)− λ̂,∇λF (fT , gT−1, λT−1)⟩ ≤ c∞∥cT−1 − b∥1 + 5(1− θ)c2∞∥λT−1 − λT−2∥2/η + 7c2∞∥λT − yT ∥2/η.
(91)
By combining (91) with (60)-(63), we can bound the left hand side of (56a) by ℓ ( π̂, λ̄(π̂) ) − ℓ(π̂, λ̂)
≤ (6c∞ − ηι) ∥∥cT−1 − b∥∥
1 + 5(1− θ)c2∞ ∥∥λT−1 − λT−2∥∥ 2 /η + 7c2∞∥λT − yT ∥2/η
+ ∣∣F (fT , gT−1, λT−1)− F ∗∣∣
≤ ( 1
6 +
1
12 +
1
12 +
1
6
) ϵ = 1
2 ϵ,
(92)
where in the last inequality we have used all the sufficient conditions (89a)-(89d). □
Theorem 22 Define ϵ′ = ϵ/(6c∞ − ηι), and set T to be
T = 8 + 48
η √ γ1ϵ′ +
( 3600(1− θ)2 + 882 ) c2∞s
ηϵ2 +
48
ηγ1ϵ = O
( c2∞ϵ −2) , (93) where γ1 is defined in (86) and we know γ1 = O(c−2∞ ). At least one of the iterations in Algorithm 3, after rounding, is an ϵ-saddle point of the EOT problem (2).
Proof. According to Lemma 21, we only need to show that (89) holds after T iterations as defined in (93).
We follow the same idea as the proof of Theorem 16. First we reduce Ẽ(f t+1, gt+1, λt+1, λt) from Ẽ(f0, g0, λ0, λ−1) = F̃ (f0, g0, λ0) to a constant s by running t1 steps. By Lemma 20, we have
t1 ≤ 1 + 6 ηγ1s − 6 ηγ1F̃ (f0, g0, λ0) . (94)
Secondly, starting from s, we continue running the algorithm, and assume that there are t2 iteration in which (89a) fails. By (42b) we have
t2 ≤ 1 + 72s
ηϵ′2 .
Therefore, we know that the total iteration number that (89a) fails can be upper bounded by
T1 = t1 + t2 ≤ 2 + 72s
ηϵ′2 +
6 ηγ1s − 6 ηγ1F̃ (f0, g0, λ0)
iterations. By choosing s = ϵ ′
6 √ γ1 , we know that
T1 ≤
{ 2 + 12η√γ1ϵ′ + 36 η √ γ1ϵ′ − 6 ηγ1F̃ (f0,g0,λ0) ≤ 2 + 48η√γ1ϵ′ if F̃ (f 0, g0, λ0) ≥ ϵ ′ 6 √ γ1 ,
2 + 12η√γ1ϵ′ + 36 η √ γ1ϵ′ − 6 ηγ1F̃ (f0,g0,λ0) ≤ 2 + 12η√γ1ϵ′ otherwise.
Therefore, we have T1 ≤ 2 + 48η√γ1ϵ′ . Similarly, from Lemma 17 we know that, starting from s, the number of iterations that (89b) and (89c) fail can be respectively bounded by
t3 ≤ 1 + 3600(1− θ)2c2∞s ηϵ2(2θ − θ2) , and t4 ≤ 1 + 3528c2∞s ηϵ2 .
By choosing s = ϵ, we have the total iteration numbers that (89b) and (89c) fail can be respectively bounded by
T2 = t1 + t3 ≤ 2+ 3600(1− θ)2c2∞ ηϵ(2θ − θ2) + 6 ηγ1ϵ − 6 ηγ1F̃ (f0, g0, λ0) ≤ 2+ 3600(1− θ) 2c2∞ ηϵ(2θ − θ2) + 6 ηγ1ϵ and
T3 = t1 + t4 ≤ 2 + 3528c2∞
ηϵ +
6 ηγ1ϵ − 6 ηγ1F̃ (f0, g0, λ0) ≤ 2 + 3528c
2 ∞
ηϵ +
6
ηγ1ϵ .
Finally, by letting s = ϵ/6 in (94), we know that
Ẽ(fT4−1, gT4−1, λT4−1, λT4−2) ≤ ϵ/6 (95) after
T4 = 1 + 36
ηγ1ϵ iterations. From (95) we know that F̃ (fT4−1, gT4−1, λT4−1) ≤ ϵ/6, which implies that (89d) holds with T = T4 by noting (42a).
Combining the above discussions, we know that after T = T1 + T2 + T3 + T4 + 1 iterations, there must exist at least one iteration such that the sufficient condition (89) holds, and thus the output of PAME is an ϵ-optimal solution to the original EOT problem (2). □

heading: D ADDITIONAL NUMERICAL RESULTS

heading: D.1 GAUSSIAN DISTRIBUTION
text: Figure 2 shows the optimal couplings obtained from the standard OT and EOT of two Gaussian distributions under three different metrics: the Euclidean cost (∥ · ∥2), the square Euclidean cost (∥ · ∥22) and the L1.51 norm (∥ · ∥1.51 ) respectively. We set n = 4, η = 0.05 and generate samples independently according to (31). For the EOT problem, we consider three agents with cost matrices computed by the three metrics mentioned above. Note that the entropy regularized models lead to a dense transportation plan and Figure 2 only plots the couplings with a probability larger than 10−3. We see that all the agents have the same total cost in the EOT model, and as expected, the cost is smaller than the other three OT costs obtained by using the same metric. The sub figures in the first row imply that if we split the workload to three parts evenly, then the three agents will each have costs 5.935/3, 2.158/3 and 5.030/3, which is not fair because they have different costs. But the EOT model can indeed guarantee the fairness.
We further compare the computational time for Gaussian distributions. We generate the data as in Section 5 and set the parameters as η = 0.5, τ = 5η/c2∞. We stop all the algorithms when the EOT error (32) is less than 10−4. Tables 1 and 2 show the CPU time (in seconds) for different (n,N) pairs. The reported computational time is averaged over 5 runs. In Table 1, the APGA algorithm fails to reach an error of 10−4 in 500000 iterations when n = 100 and n = 500. We conclude that the APGA algorithm converges much slower than PAM and PAME algorithms. The PAME algorithm performs the best among all three algorithms.

heading: D.2 NUMERICAL RESULTS ON FRAGMENTED HYPERCUBE DATASET
text: In this section, we compare the performance of PAME with PAM and APGA (25) Scetbon et al. (2021) on the fragmented hypercube dataset.
Fragmented Hypercube: We now consider transferring mass between a uniform distribution over a hypercube µ = U([−1, 1]d) and a distribution ν obtained by a pushforward ν = T♯µ defined by T (x) = x + 2sign(x) ⊙ (∑m∗ m=1 em ) . Here sign(·) is taken elementwisely, m∗ ∈ [d] and ei, i ∈ [d] is the canonical basis of Rd. In our experiments, we set d = 10,m∗ = 2 and sample two base support sets {xbasei }i∈[n], {ybasej }j∈[n] independently from µ, ν. To obtain the cost matrix for one agent, we first add Gaussian noise sampled from N (0, 1) to the base support sets to get {xnoisyi }i∈[n], {y noisy j }j∈[n] and compute the cost using the noisy support sets. For instance, for the k-th agent, we have (xnoisyi ) k = xbasei + N (0, 1) , (y noisy j )
k = ybasej + N (0, 1) and Cki,j = ∥(xnoisyi )k − (y noisy j )
k∥22. Figures 3 plots the EOT error versus the CPU time for Fragmented Hypercube dataset. We run PAM for 20000 iterations to get an approximate optimal ℓ∗ and run all algorithms for 2000 iterations for different parameter settings. In all cases, the PAME and PAM perform significantly better than APGA, and PAME also shows significant improvement over PAM.
We then compare the computational time for Fragmented Hypercube dataset. We set the parameters as η = 0.2, τ = 5η/c2∞. We stop all the algorithms when the EOT error (32) is less than 10
−4. Tables 3 and 4 show the CPU time (averaged over 5 runs) for different (n,N) pairs. We see that the PAME algorithm still performs the best among all three algorithms. Note that in Table 3 the APGA algorithm fails to reach an error of 10−4 in 500000 iterations when n = 50, n = 100 and n = 500.
Table 2: CPU time (in seconds) comparison for Gaussian Distributions. Fixed n = 50.
Algorithms N = 2 N = 3 N = 5 N = 10 N = 20 PAM 0.180343 0.177209 1.021598 0.719909 0.903429 PAME 0.105775 0.091593 0.560785 0.385495 0.618381 APGA 70.959300 106.768840 169.213121 178.637889 212.697866
Figure 3: CPU time comparison between PAM, PAME and APGA algorithms on the Fragmented Hypercube dataset. Upper Left: N = 5, n = 100, η = 0.2, Upper Right: N = 5, n = 500, η = 0.2, Bottom Left: N = 5, n = 100, η = 0.1, Bottom Right: N = 10, n = 100, η = 0.2.

Review summaries:
This paper develops new algorithms, convergence rate analysis and additional tools for the equitable and optimal transport problem. Following the same history of OT, where we moved from the OT formulation, to the entropic regularization, to primal-dual analysis to explicit convergence analysis, this manuscript extends the idea of EOT and provides useful tools with provable performance for the practitioner. The contributions of this paper are on the complexity analysis, as well as new algorithms for the EOT problem.

The paper addresses the equitable and optimal transport (EOT) problem. One approach to this problem involves perturbation via the addition of an entropy regularization. For this formulation, researchers have proposed a projected alternating maximization algorithm (PAM). The authors of this paper provide a convergence analysis of the PAM algorithm. They also introduce a novel rounding procedure to construct the primal solution of the EOT problem. The paper also contributes a variant of PAM that leads to numerical performance improvements. Numerical experiments are provided for a synthetic dataset.

The authors focus on the equitable and optimal transport (EOT) problem. In particular, EOT is a linear programming (LP) problem, which is expensive to solve in practice. The existing work suggests adding an entropy term to the objective and solving the resulting problem in the dual domain via an algorithm named Projected Alternating Maximization algorithm (PAM). However, there are two main issues with PAM: that is its convergence is unknown and it only provides a solution in the dual domain. The authors precisely address these issues by providing a convergence rate for PAM and proposing a novel method to compute a feasible primal solution from the dual solution. Furthermore, they propose a variant of PAM, i.e., PAM with Extrapolation (PAME), which appears to perform better in practice. 

This paper considers the problem of approximate optimal transport, where population probability measures are replaced by discrete PMFs at collections of samples, and there is an additional "equitable" constraint that seeks to ``fairly" distribute the probability mass. This admits a formulation as a linear program and may be solved efficiently through its dual. However, recovering the primal solution is nontrivial, and there exists no rate analysis to-date. This work develops a stochastic dual ascent method to solve an entropy-regularized version of this LP, where entropy regularization is introduced to ensure strong concavity, as well as avoid a notion of over-fitting. The crux of this method is an alternating block coordinate ascent algorithm, which is then combined with an extrapolation/Nesterov acceleration scheme to accelerate convergence. Numerical experiments demonstrate the practical utility of the proposed approach.