On Robust Optimal Transport: Computational Complexity and Barycenter Computation

Abstract:
We consider robust variants of the standard optimal transport, named robust optimal transport, where marginal constraints are relaxed via Kullback-Leibler divergence. We show that Sinkhorn-based algorithms can approximate the optimal cost of robust optimal transport in Õ( 2 ε ) time, in which n is the number of supports of the probability distributions and ε is the desired error. Furthermore, we investigate a fixed-support robust barycenter problem between m discrete probability distributions with at most n number of supports and develop an approximating algorithm based on iterative Bregman projections (IBP). For the specific case m = 2, we show that this algorithm can approximate the optimal barycenter value in Õ( 2 ε ) time, thus being better than the previous complexity Õ( 2 ε2 ) of the IBP algorithm for approximating the Wasserstein barycenter.

text: 2
ε ) time, in which n is the number of supports of the probability distributions and ε is the desired error. Furthermore, we investigate a fixed-support robust barycenter problem between m discrete probability distributions with at most n number of supports and develop an approximating algorithm based on iterative Bregman projections (IBP). For the specific case m = 2, we show that this algorithm can approximate the optimal barycenter value in Õ(mn 2
ε )
time, thus being better than the previous complexity Õ(mn 2
ε2 ) of the IBP algorithm for approximating the Wasserstein barycenter.

heading: 1 Introduction
text: The recent advance in computation with optimal transport (OT) problem [12, 3, 13, 7, 22, 26, 20] has led to a surge of interest in using that tool in various domains of machine learning and statistics. The range of its applications is broad, including deep generative models [4, 16, 36], scalable Bayes [33, 34], mixture and hierarchical models [24], and other applications [32, 29, 10, 17, 37, 35, 8].
The goal of optimal transport is to find a minimal cost of moving masses between (supports of) probability distributions. It is known that the estimation of transport cost is not robust when there are outliers. To deal with this issue, [38] proposed a trimmed version of optimal transport. In particular, they search for truncated probability distributions such that the transport cost between them is minimized. However, their trimmed optimal transport is non-trivial to compute, which hinders its usage in practical applications. Another line of works proposed using unbalanced optimal transport (UOT) to solve the sensitivity of optimal transport to outliers [5, 31]. More specifically, their idea is to assign as small as possible masses to outliers by relaxing the marginal constraints of OT through a penalty function such as the Kullback-Leibler (KL) divergence. This direction of robust optimal transport has been shown to have good performance in generative models and domain adaptation [5]. Although this approach achieved considerable success, the full picture of its computational complexity has remained missing.
? Khang Le and Huy Nguyen contributed equally to this work.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Our Contribution: In the paper, we provide a comprehensive study of the computational complexity of robust optimal transport and its corresponding barycenter problem when the probability distributions are discrete and have at most n components. Our contribution is twofold and can be summarized as follows:
(1) On robust optimal transport, we consider two versions corresponding to two ways of relaxing marginal constraints in the standard optimal transport problem via the KL divergence. We show that two scaling algorithms computing these robust formulations have the complexities Õ(n2/ε), where ε denotes the desired error for the computed cost. These complexities are lower than the complexity of the Sinkhorn algorithm for solving the optimal transport problem, which is Õ(n2/ε2) [13], and match the complexity of the Sinkhorn algorithm that solves the UOT problem [27]. Furthermore, we show how the above complexity can be improved by utilizing the low-rank approximation method to speed up the matrix-vector computations in the loop similar to [2], and obtain the improved computing time of Õ(nr2 + nrε ), where r is the approximated rank.
(2) On robust barycenter problem, where the goal is to determine a probability measure that minimizes its robust optimal cost to a given set of m ≥ 2 probability measures, we propose ROBUSTIBP algorithm for solving the robust barycenter problem, which is inspired by the iterative Bregman projection (IBP) algorithm for solving the traditional barycenter problem [6]. We show that whenm = 2, the complexity of ROBUSTIBP algorithm is at the order of Õ(mn2/ε), better than that of the IBP algorithm for solving the traditional barycenter problem [19], which is Õ(mn2/ε2). To the best of our knowledge, the ROBUSTIBP is also the first practical algorithm obtaining the near-optimal complexity Õ(mn2/ε) for solving the barycenter problem even under only the setting m = 2.
Organization: The paper is organized as follows. In Section 2, we provide the background on the optimal transport problem and some of its variants that have robust effects. In Section 3, we discuss in-depth the variant where only one marginal constraint is relaxed, study the computational complexity of a Sinkhorn-based algorithm that solves it, and then briefly introduce the fully-relaxed formulation. We also establish the complexities of these algorithms after applying Nyström method. Subsequently, we present our study of the robust barycenter problem in Section 4. In Section 5, we carry out empirical studies to illustrate the theories before concluding with a few discussions in Section 6. The proofs of our theoretical results are in the supplementary material.
Notation: We let [n] stand for the set {1, 2, . . . , n} while Rn+ indicates the set of all vectors with nonnegative entries. For a vector x ∈ Rn and p ∈ [1,∞), we denote ‖x‖p as its `p-norm and diag(x) as the diagonal matrix with x on the diagonal. The natural logarithm of a vector a = (a1, ..., an) ∈ Rn+ is denoted by log a = (log a1, ..., log an), 1n stands for a vector of length n that all of its entries equal to 1, and ∂xf refers to the partial differentiation of function f with respect to x. For any given space X ⊂ Rd, we denote by P(X ) the space of all probability measures on X . Given an integer n > 0 and a real number ε > 0, the notation a = O (b(n, ε)) means that a ≤ C · b(n, ε) where C is independent of n and ε. Meanwhile, the notation a = Õ(b(n, ε)) indicates the previous inequality may depend on a logarithmic function of n and ε. For any two probability measures x = (x1, . . . , xn) and y = (y1, . . . , yn) with the same supports, the generalized Kullback-Leibler divergence is defined as KL(x‖y) = ∑n i=1 [ xi log ( xi yi ) − xi + yi ] . Finally, the entropy of a matrix
X is given by H(X) = ∑n i,j=1−Xij(logXij − 1).

heading: 2 Background on Optimal Transport
text: In this section, we review optimal transport and its unbalanced formulation, then from that deriving formulations for robust optimal transport. For any P and Q in P(X ) for a space X , the OT distance between P and Q takes the following form
OT(P,Q) := min π∈Π(P,Q)
∫ c(x, y)dπ(x, y), (1)
where Π(P,Q) is the set of joint probability distributions in X × X such that their marginal distributions are P and Q, and c : X × X → [0,∞) is a cost function.
Unbalanced Optimal Transport: When P or Q is not a probability distribution, the OT formulation between P and Q in equation (1) is no longer valid. One solution to this issue is using the unbalanced optimal transport (UOT) [9], which is given by:
UOT(P,Q) := min π∈M+(X×X )
∫ c(x, y)dπ(x, y) + τ1KL(π1‖P ) + τ2KL(π2‖Q), (2)
whereM+(X × X ) denotes the set of joint non-negative measures on the space X × X ; π1, π2 are the marginal distributions of π and respectively correspond to P and Q; τ1, τ2 are regularized positive parameters. Note that, we can replace the KL divergence in equation (2) by any Csiszár-divergence [11]. However, we only consider the case of KL divergence in this work.
Robust Optimal Transport: Optimal transport is well-known for not being robust in the present of outliers. A way to deal with this issue is using the approach of unbalanced optimal transport (UOT), which has demonstrated favorable practical performance in generative models and domain adaptation [5]. More specifically, when P and Q are probability distributions in X , the Robust Unconstrained Optimal Transport (ROT) admits the following form
ROT(P,Q) := inf P1,Q1∈P(X ) min π∈Π(P1,Q1)
∫ c(x, y)dπ(x, y) + τ1KL(P1‖P ) + τ2KL(Q1‖Q), (3)
where τ1, τ2 > 0 are some given regularized parameters. The reason to name it robust unconstrained optimal transport is that instead of looking for an optimal transport plan moving masses from P to Q, we seek another plan that optimally transports masses between their approximations, which are probability measures P1 and Q1, under the KL divergence. This formulation is closely related to the ones studied in [5] and [23]: the former used χ2-divergence for the relaxation and the latter used total variation distance (note that those three divergences all together belong to the family of f -divergence).
By relaxing only one marginal constraint regarding (presumably) on P , we have another version of ROT, named Robust Semi-constrained Optimal Transport (RSOT), which is given by
RSOT(P,Q) := inf P1∈P(X ) min π∈Π(P1,Q)
∫ c(x, y)dπ(x, y) + τKL(P1‖P ), (4)
where τ > 0 is a regularized parameter. We could also define RSOT(Q,P ) similarly with a remark that although RSOT(P,Q) can be different from RSOT(Q,P ), the techniques for obtaining the computational complexity of both are similar.
UOT vs ROT/RSOT: Though the formulations ROT/RSOT and UOT seem to be similar, they serve different purposes. The goal of UOT is to deal with unbalanced measures, thus there is no condition on the “transport plan”. Hence, the meaning of the optimal plan π of UOT problem is dependent on the interpreter. For example, in applications such as [30], the UOT is used to figure out the developmental trajectory of cells. Meanwhile, ROT/RSOT aim to seek an accurate transport plan between two possibly corrupted probability distributions. The toy example in Figure 1 illustrates this difference. In particular, the marginals of the “transport plan" obtained by the latter (see plots (b), (d)) are very different from the two original probability measures a,b. On the other hand, the solution of the former leads to good approximations of a and b (see plots (a), (c)) while removing some bumps in both tails which are presumably outliers.

heading: 3 Discrete Robust Optimal Transport and its Computational Complexity
text: When P and Q are discrete measures, the KL penalties in equations (3) and (4) suggest that the probability distributions P1 and Q1 need to share the same set of supports as that of P and Q, respectively. Therefore, throughout this section, we implicitly require this condition in our formulations of RSOT and ROT and we denote the masses of P and Q by a and b, respectively.

heading: 3.1 Robust Semi-constrained Optimal Transport
text: Assume that the marginal constraint associating with Q is kept and that of P is relaxed and P1 and P share the same set of supports, the formulation of RSOT in equation (4) can be rewritten as follows
min X∈Rn×n+ ,X>1n=b
frsot(X) := 〈C,X〉+ τKL(X1n||a), (5)
Algorithm 1: ROBUST-SEMISINKHORN Input: C,a,b, η, τ, niter Initialization: u0 = v0 = 0, k = 0 while k < niter do ak ← B(uk, vk)1n, bk ← ( B(uk, vk) )> 1n
if k is even then uk+1 ← ητη+τ
[ uk η + log(a)− log(a k) ]
vk+1 ← vk else uk+1 ← uk vk+1 ← η [ vk η + log(b)− log(b k) ] end if k ← k + 1
end while return B(uk, vk)
where a,b are the masses of P and Q respectively, and C is the cost matrix whose entries are distances between the supports of these distributions. Solving directly problem (5) by traditional linear programming solvers can be expensive and not scalable in terms of n. Therefore, we utilize the entropic regularization approach proposed by [12] to the objective function of RSOT, leading to
min X∈Rn×n+ ,X >1n=b grsot(X) := frsot(X)− ηH(X). (6)
Here, η > 0 is a given regularization parameter, and we refer the problem (6) to as entropic RSOT. The dual problem of entropic RSOT is
min u,v∈Rn
hrsot(u, v) := η‖B(u, v)‖1 + τ 〈 e−u/τ ,a 〉 − 〈 v,b 〉 , (7)
where B(u, v) is defined as a matrix of size n × n with entries [B(u, v)]ij := e(ui+vj−Cij)/η. Since equation (7) is an unconstrained convex optimization problem, we can perform alternating minimization for u and v by setting ∂h(u, v)/∂u = 0 and ∂h(u, v)/∂v = 0, resulting in closed-form updates of a Sinkhorn-like procedure (see [12]) in Algorithm 1. This procedure is known to converge to the optimal solution (u∗, v∗) := arg minhrsot(u, v). As strong duality holds for the convex optimization problem (6), the optimal transport plan of the entropic RSOT is exactly B(u∗, v∗).
Since no assumptions are made on the cost matrix, except its entries are non-negative, closed-form solutions of OT and UOT generally do not exist. Therefore, we introduce the definition of an ε-approximation solution of an optimization problem, which will be used for all the subsequent complexity analyses.
Definition 1 (ε-approximation). For any ε > 0, a transportation plan X is called an εapproximation of the minimizer X̂ of some objective function f if f(X) ≤ f(X̂) + ε.
Based on this concept, we then state our main theorem on the runtime complexity of Algorithm 1 in solving the RSOT problem (5). Theorem 1. For Ursot := max{3 log(n), ε/τ} and η = ε/Ursot, Algorithm 1 returns an εapproximation of the optimal solution X̂rsot of the problem (5) in time
O ( τn2
ε log(n)
[ log ( τ‖C‖∞ ε ) + log(log(n)) ]) .
Proof Sketch. The full proof of Theorem 1 is in Appendix B. Note that, this result is not achieved by directly applying Theorem 2 in [27] with τ2 →∞ as the nature of the dual function changes in that limit, invalidating many previous results. Let Xkrsot be the output of Algorithm 1 at the k-th step while X̂rsot and X∗rsot denotes the minimizers of equations (5) and (6), respectively. The goal is to find k that guarantees frsot(Xkrsot)− frsot(X̂rsot) ≤ ε = ηUrsot. We start by decomposing
frsot(X k rsot)︸ ︷︷ ︸
grsot(Xkrsot)+ηH(X k rsot) − frsot(X̂rsot)︸ ︷︷ ︸ grsot(X̂rsot)+ηH(X̂rsot)
≤ [ grsot(X k rsot)− grsot(X∗rsot) ] + η [ H(Xkrsot)−H(X̂rsot) ] ,
and try to bound each term by a linear function of η. Dealing with the entropy term is simple as the η factor is already presented, and the entropy difference can be bounded by a constant due to the fact that 1 ≤ H(X) ≤ 2 log(n) + 1 for all X ∈ Rn×n+ , ‖X‖1 = 1. The non-trivial part is bounding the difference between grsot values, which hinges upon two results. The first one is the value of grsot at optimality:
grsot(X ∗ rsot) = −η − τ(1− α) + 〈v∗rsot, b∗rsot〉. (8)
The second result is the geometric convergence rate of the updates on u and v (Lemma 6 in Appendix B):
max { ‖uk+1 − u∗‖∞, ‖vk+1 − v∗‖∞ } ≤ (const) ( τ τ + η )k/2 =: ∆k.
The final step is using equation (8) to tailor the grsot difference to be bounded by a linear function of ∆k, which is an exponential function of k, then solving for the minimum k at which this exponential function is small enough compared to η. The main technical difficulty here is to deal with the unknown term 〈v∗rsot, b∗rsot〉 in equation (8), which causes the deviation from the previous techniques.
Remark 1. The result of Theorem 1 indicates that the complexity of ROBUST-SEMISINKHORN algorithm for computing RSOT is at the order of Õ(n 2
ε ). This complexity is near-optimal and faster than the complexity of the standard Sinkhorn algorithm for computing the optimal transport problem [13, 22], which is at the order of Õ(n 2
ε2 ).

heading: 3.2 Robust Unconstrained Optimal Transport
text: In this section, we briefly present another version of robust optimal transport, abbreviated by ROT, when two distributions are contaminated. We first show that the approach of using the duality of the objective function of ROT problem with entropic regularizer does not produce a Sinkhorn algorithm as in the cases of RSOT and UOT. However, a second thought of the problem finds an interesting link between the optimal solutions of ROT and UOT, which results in a nice algorithm for the ROT. We also discuss some technical difficulties when analysing the complexity for the ROT problem. At the end of this section, we show that the result could be extended to the case of low-rank cost matrix, which will significantly reduce the computation.
Recall that the masses of P and Q are a and b, respectively, the ROT problem (3) becomes
min X∈Rn×n+ ,‖X‖1=1
frot(X) := 〈C,X〉+ τKL(X1n||a) + τKL(X>1n||b). (9)
Here we set τ1 = τ2 = τ for the sake of simplicity, since there are no more technical difficulties to work with finite τ1 6= τ2. As noted in Section 2, the formulation (9) bears some resemblance to the unbalanced optimal transport problem studied in [27], except the additional norm condition forcing X to be a transportation plan (i.e., a joint probability distribution), which shows the different nature
of two problems. Following the approach of using the Sinkhorn algorithm of UOT, the duality of formulation (9) has the form
η log ‖B(u, v)‖1 + τ { 〈eu/τ ,a〉+ 〈ev/τ ,b〉 } .
By taking derivatives of the above function with respect to u and v and set the derivatives to be zero, we obtain
B(u, v)1n ‖B(u, v)‖1
= e−u/τ a, B(u, v) >1n
‖B(u, v)‖1 = e−v/τ b,
where denotes element-wise multiplication. Unfortunately, the above equations do not have closedform solutions to produce update as the Sinkhorn algorithms do because of the term ‖B(u, v)‖1 in the denominator. However, the objective function of UOT is not homogeneous with respect to X , but could be written as a linear function of ROT and another function of ‖X‖1 due to some special properties of the KL divergence. This observation leads to the interesting result summarized in the below lemma.
Lemma 1 (Connections with UOT). The optimal solution of problem (9), denoted X∗rot, is the normalized version of X∗uot which is the minimizer of UOT in entropic formulation. More specifically, we have X∗rot = X∗uot ‖X∗uot‖1 .
The proof of Lemma 1 is in Appendix D. Based on this result, we can utilize the Sinkhorn algorithm that solves UOT (see [27]) with a normalizing step at the end to produce a solution for the ROT. Although the normalizing step is convenient in finding ROT’s solution, it introduces new challenge in the proof compared to that of UOT since the normalizing constant does not have a lower bound. Even so, we are still able to obtain an ε-approximation solution for the ROT in Õ(n2/ε) time without any additional constraints on the setting. For more technical details, please refer to Appendix D.
Further Improving Complexities by Low-Rank Approximation: As a consequence of our complexity analysis, we can show that by using low-rank approximation method studied in [2] to the kernel matrix K := exp(−C/η), we could further reduce the complexities of both robust semi/unconstrained optimal transport problem to Õ(nr2 + nr/ε) time, given the same ε-approximation and the approximated-rank r. This result is essentially different from the complexity studied in [2], where the ε-approximation is considered regarding the optimal value of the entropic-regularized problem, not the original one in our analysis. For a more detailed discussion, please refer to Appendix E.

heading: 4 The Robust Barycenter Problem
text: In this section, we consider the problem of computing the barycenter of a set of possibly corrupted probability measures. The semi-constrained formulation arises as a natural candidate for this goal, when potential outliers only appear in the given probability measures and the desired barycenter is the barycenter of the uncontaminated probability measures. In particular, assume that we have m ≥ 2 discrete probability measures P1, . . . , Pm: each has at most n fixed support points and the associated positive weights are given by ω1, . . . , ωm ( ∑m i=1 ωi = 1). The barycenter problem then aims to
find the probability measure that minimizes ∑m i=1 ωiRSOT(Pi, P ), which is a linear combination of RSOT divergence from the barycenter to all given probability measures. We refer it as Robust Semi-constrained Barycenter Problem (RSBP). In this work, we consider the fixed-support settings where all the probability measures Pi share the same set of support points. This setting had been widely used in the previous works to study the computational complexity of Wasserstein barycenter problem [19, 21]. Let pi be the mass of probability measure Pi for i ∈ [m], the discrete RSBP reads
min p∈Rn+,‖p‖1=1 m∑ i=1 ωi [ min Xi∈Rn×n+ ,X>i 1n=p 〈Ci, Xi〉+ τKL(Xi1n‖pi) ] ,
which is equivalent to
min X∈D1(X) frsbp(X) := m∑ i=1 ωi [ 〈Ci, Xi〉+ τKL(Xi1n‖pi) ] , (10)
Algorithm 2: ROBUSTIBP Input: {Ci}mi=1, {pi}mi=1, τ, η, niter Initialization: u0i = v0i = 0n for i ∈ [m], k = 0 while k < niter do aki ← B(uki , vki ;Ci)1n; bki ← ( B(uki , v k i ;Ci) )> 1n ∀i ∈ [m]
if k is even then uk+1i ← ητ η+τ [uki η + log(pi)− log(a k i ) ] ∀i ∈ [m]
vk+1i ← vki ∀i ∈ [m] else uk+1i ← uki ∀i ∈ [m] vk+1i ← η [ vki η − log(b k i )− ∑m t=1 ωt( vkt η − log(b k t )) ] ∀i ∈ [m] end if k ← k + 1
end while Xki ← B(uki , vki ;Ci) ∀i ∈ [m] return (Xk1 , . . . , Xkm) for equation (14) or ( Xk1 ‖Xk1 ‖1 , . . . , Xkm ‖Xkm‖1 ) for equation (11).
where D1(X) := {
(X1, . . . , Xm) : Xi ∈ Rn×n+ and ‖Xi‖1 = 1 ∀i ∈ [m]; X>i 1n = X>i+11n ∀i ∈ [m− 1] } . Note that the objective function of RSBP is different from that of Wasserstein barycenter [19]: here we relax the marginal constraints Xi1n = pi by using the KL divergence to deal with the contaminated Pi. Finally, the constraints X>i 1n = X > i+11n = p are to guarantee that the transportation plans Xi have one common marginal which turns out to be a feasible barycenter p. Similar to RSOT, we consider an entropic-regularized formulation of (10), named entropic RSBP:
min X∈D1(X) grsbp(X) := m∑ i=1 ωigrsot(Xi; pi, Ci). (11)
Since some functions like grsot(X), depends on some parameters like Ci and pi, we sometimes abuse the notation by including these parameters next to variables, e.g., grsot(Xi;Ci,pi). A general approach to deal with (11) is to consider its dual function, which admits the following form:
min u=(u1,...,um),v=(v1,...,vm)∑m
i=1 ωivi=0
hrsbp(u,v) := m∑ i=1 ωi [ η log ‖B(ui, vi;Ci)‖1 + τ 〈 e−ui/τ ,pi 〉] . (12)
We could use the alternating minimization method to find the minimizer of (12). In particular, starting at an initialization u0 and v0, we update them alternatively as follows:
uk+1 = arg min u hrsbp(u,v k), vk+1 = arg min v: ∑m
i=1 ωivi=0
hrsbp(u k+1,v). (13)
In some problems (e.g., RSOT), closed-form updates can be acquired if the system of equations ∂hrsbp(u,v
k)/∂u = 0 and ∂hrsbp(uk,v)/∂v = 0 could be solved exactly by some simple formulas. However, this is not the case with the formulation of hrsbp in equation (12) because the logarithmic term leads to an intractable system of equations of the partial derivative of hrsbp. Instead, we propose to solve the optimization problem (11) via another objective function, whose dual form can be solved effectively by alternating minimization.

heading: 4.1 ROBUSTIBP Algorithm
text: We consider a similar problem to the entropic RSBP in (11), with its feasible set D(X) := {(X1, . . . , Xm) : Xi ∈ Rn×n+ ,∀i ∈ [m];X>i 1n = X>i+11n∀i ∈ [m− 1]} which does not have the
norm constraint. The primal objective function and its dual are as follows:
Primal: min X∈D(X) grsbp(X) := m∑ i=1 ωigrsot(Xi; pi, Ci), (14)
Dual: min u,v: ∑m i=1 ωivi=0 h̄rsbp(u,v) := m∑ i=1 ωi [ η‖B(ui, vi;Ci)‖1 + τ 〈 e−ui/τ ,pi 〉] . (15)
The dual formulation (15) has a closed form updates for u and v. Based on these, we develop Algorithm 2, namely ROBUSTIBP, since this procedure resembles the iterative Bregman projections studied in [6] and [19]. The updates of u and v are known to converge to the optimal solution (u∗,v∗) of the problem (15), and strong duality suggests that X∗ = (B(u∗i , v ∗ i ;Ci)) m i=1 is the optimal solution of the problem (14). Furthermore, there is an intriguing relation between the optimal solution of the problem (14) to that of the problem (11), presented in the following lemma. Lemma 2. Let X̄∗ = (X̄∗1 , . . . , X̄∗m) and X∗ = (X∗1 , . . . , X∗n) be the optimizers of grsbp with the feasible set D(X) and with the feasible set D1(X), respectively. Then, X∗i = X̄∗i ‖X̄∗i ‖1 for all i ∈ [m].
The proof of Lemma 2 is in Appendix C. This result indicates that we can approximate the solution of equation (11) by the solution of equation (14), using the same Algorithm 2 with an additional normalizing step at the end.

heading: 4.2 Complexity Analysis
text: In this section, we provide the analysis of ROBUSTIBP algorithm for obtaining an ε-approximation of the robust semi-constrained barycenter problem (11) when m = 2. We also discuss the challenges of extending the current proof technique to m ≥ 3 at the end of this section. First, we present the complexity of the ROBUSTIBP algorithm in the following theorem.
Theorem 2. For m = 2 and η = εU−1rsbp where Ursbp := max{2 + 2 log(n), 2ε, 3ε log(n)/τ}, the ROBUSTIBP algorithm returns an ε-approximation of the optimal solution (X̂1, . . . , X̂m) of the
RSBP (10) in time O (τn2 ε log(n) [ log ( τ m∑ i=1 ‖Ci‖∞ ) + log ( log(n) ε )]) .
Remark 2. The complexity Õ(n2/ε) of ROBUSTIBP algorithm is near-optimal and better than that of IBP algorithm for solving the Wasserstein barycenter problem, which is Õ(n2/ε2) when m = 2 in [19]. It is also better than the complexity of FASTIBP algorithm in [21], which is Õ(n7/3/ε4/3). To the best of our knowledge, the ROBUSTIBP is also the first practical algorithm obtaining the near-optimal complexity Õ(n2/ε) for solving the barycenter problem under the setting m = 2.
The main ingredient in the proof of Theorem 2 is the convergence rate of vectors u and v of the problem (15), which is captured as follows:
max { m∑ i=1 ‖∆uk+1i ‖∞, m∑ i=1 ‖∆vk+1i ‖∞ } ≤ (constant) ( τ τ + η )k/2 , (16)
where ∆uki := u k+1 i − u∗i and ∆vki := v k+1 i − v∗i . The result can be achieved by alternatively applying two following inequalities.
For the first inequality, with even k, from the update of uk+1 in the Algorithm 2, we obtain ‖∆uk+1i ‖∞ ≤ ττ+η‖∆v k i ‖∞.
The second inequality is obtained from the update of vk in Algorithm 2 as follows: m∑ i=1 ‖∆vki ‖∞ ≤ m∑ i=1 ((m− 2)ωi + 1)‖∆uk−1i ‖∞.
Thus, when m = 2, we can achieve inequality (16), though this approach is inapplicable for the case m > 2. For a formal statement regarding the above convergence rate, please refer to Lemma 12 in
Appendix C. Note that for m ≥ 3, the result of Theorem 2 still holds if uk and vk converge at the rate of the order ( ττ+η )
k/2. So next we will take a closer look at this case to see whether the rate remains geometric. On m ≥ 3: In Figure 2, we plot the values of two ratios: Ruv := ∑m i=1 ‖∆u k+1 i ‖∞∑m
i=1 ‖∆vki ‖∞ and Ruu :=∑m
i=1 ‖∆u k+1 i ‖∞∑m
i=1 ‖∆u k−1 i ‖∞ . When k is even, we have that Ruu ≤ ττ+η for all m, while the inequality Ruv ≤ τ τ+η was only proved for the case m = 2. From this figure, both these bounds are true in all considered cases. However, while the bound on Ruv (which is theoretically true for all m) is only tight when m = 2 and seems to be loose in several trials with larger values of m, the bound Ruu (which is only showed for the case m = 2) appears to be tight in all reported scenarios. Thus, we conjecture that the geometric convergence rate at equation (16) may still hold for m greater than 2. We leave the case m ≥ 3 for the future work.

heading: 5 Experiments
text: In this section, we provide numerical evidences regarding our presented complexities for ROBUSTSEMISINKHORN and ROBUST-IBP algorithms. We put additional experiments (including the runtime comparison of ROT/RSOT on synthetic and real datasets, as well as some applications for the studied robust formulations) in Appendix F. All the optimal solutions for convex problems in the following part are computed using the cvxpy library [1]. All the experiments are conducted on a server with 32 GB RAM, 8 cores Intel(R) Core(TM) i7-9700K and 1 GeForce RTX 2080 GPU.
Runtime Demonstration: For each algorithm, we investigate the number of iterations required to obtain an ε-approximation. We compare the theoretical values in Theorems 1 and 2 with the empirical values computed by running the corresponding algorithms to obtain the first iterations from where the algorithm always returns an ε-approximation.
For RSOT, we let n = 100, τ = 1, generate entries of C uniformly from the interval [1, 50] and draw entries a, b uniformly from [0.1, 1] then normalizing them to form probability vectors. η is set according to Theorem 1. For each ε varying from 5× 10−2 to 5× 10−5, we calculate the number of theoretical and empirical iterations described above, as well as their ratio. This experiment is run 10 times and we report their mean and standard deviation values in Figure 3 (a). We also carry out a similar experiment on MNIST data, which is reported in the Appendix F.
For RSBP, we run the ROBUSTIBP algorithm with the following setup: n = 10; τ = 1; p1, . . . , pm, [ω1, . . . , ωm] are randomly-initialized probability vectors; {Ci}mi=1 is a set of n× n matrices whose entries drawn uniformly in [0.01, 0.1]; five chosen values of ε vary from 10−3 to 10−5 (which are relatively small compared to the optimal cost frsbp(X∗) is about 0.019± 0.001 when m = 2 and is about 0.021± 0.001 when m = 3); and the corresponding values of η are set according to Theorem 2. The results are shown in Figure 3 (b) and (c). Note that the complexity for the case m ≥ 3 is still an open problem, and we use the formula in Theorem 2 to compute the (hypothetical) theoretical number of iterations in that case.
In all three experiments, it is noticeable that the ratios between theoretical and empirical values decrease as ε→ 0, indicating the our complexity bounds get tighter.

heading: 6 Conclusion
text: In the paper, we study the complexity of Sinkhorn-based algorithms for approximately solving robust versions of optimal transport between two discrete probability measures with at most n components, and show that they return ε-approximated solutions in Õ(n2/ε) time. Low-rank approximation technique is also analysed to further reduce the dependency of these complexities on n, resulting in Õ(nr2 + nr/ε) complexities. Finally, we investigate a robust barycenter problem between m probability measures and develop the IBP-based algorithm for solving it. When m = 2, the complexity of the ROBUSTIBP algorithm is proved to be at the order of Õ(mn2/ε), while in the case m ≥ 3 we believe that a novel proof technique needs to be developed to establish the geometric convergence of the updates from the algorithm. We leave this direction for the future work.

Review summaries:
This paper proposes two robust formulation of the optimal transport problem (namely RSOT and ROT) based on the relaxation of the marginal constraints of the optimal transport plan using Kullback-Leibler divergence. The authors introduces algorithms based on Sinkhorn iterations in order to solve RSOT and ROT problems, and analyse their computational complexities. The fixed-support Wasserstein barycenter problem is also tackled using the proposed RSOT as a divergence, and the algorithm shows a better complexity than the one based on entropy regularised OT. Numerical experiments are also provided to support the theoretical results and applied task are performed (color transfer and generative modeling).

This paper presents the complexity analysis of two problems: robust optimal transport and robust barycenter problem. For both these problems, the paper analyzes a sinkhorn-based algorithm and obtains a complexity of O(n^2 / \epsilon). Finally, the authors show some numerical analysis to support the theoretical claims.  

The paper proposes the Robust Optimal Transport(ROT), Robust Semi-constrained Optimal Transport (RSOT), Robust Unconstrained Optimal Transport (RUOT) and Robust Semi-constrained Barycenter Problem (RSBP). In these problems, the equality constraints in the original OT problem is replaced by the soft KL divergence. Then the Sinkhorn-like algorithms are proposed to solve the corresponding dual problems. Finally, the detailed convergence analysis is given. 

This paper focuses on the computational complexity of robust optimal transport between discrete probability measures. The authors propose two Sinkhorn-based algorithms with near-optimal complexity to approximate the optimal cost. Furthermore, for solving a robust barycenter problem, they propose an algorithm based on iterative Bregman projections (IBP) and show that it also has near-optimal complexity between two discrete probability distributions. Experiments verify the presented complexities for the above algorithms. The main limitations and future work directions are also discussed.