Finding Wasserstein Ball Center: Efficient Algorithm and The Applications in Fairness
Authors: Yuntao Wang
Yuxuan Li
Qingyuan Yang
Hu Ding
Venue: ICML 2025 poster (2025)
Keywords: Wasserstein barycenter
fairness
Primary area: optimization
TL;DR: We introduce and solve a new type of probability center that retains more information than Wasserstein barycenter by considering fairness.

Abstract: Wasserstein Barycenter (WB) is a fundamental geometric optimization problem in machine learning, whose objective is to find a representative probability measure that minimizes the sum of Wasserstein distances to given distributions. WB has a number of applications in various areas.  However, 
WB may lead to unfair outcome towards underrepresented groups in some applications (e.g., a "minority'' distribution may be far away from the obtained WB under Wasserstein distance). 
To address this issue, we propose an alternative objective called  "Wasserstein Ball Center (WBC)''. Specifically, WBC is a distribution that encompasses all input distributions within the minimum Wasserstein distance, which can be formulated as a ``minmax'' optimization problem.
 We show that the WBC problem with fixed support is equivalent to solving a large-scale linear programming (LP) instance, which is quite different from the previously studied LP model for WB. By incorporating some novel observations on the induced normal equation, we propose an efficient algorithm that accelerates the interior point method by $O(\min(N^2m, Nm^2, m^4))$ times ("$N$'' is the number of distributions and 
 "$m$'' is the support size).  Finally,  we conduct a set of experiments on  both synthetic and real-world datasets, demonstrating the computational efficiency  of our algorithm, and showing its ability to provide more fairness for input distributions.

Finding Wasserstein Ball Center: Efficient Algorithm and The Applications in Fairness
Yuntao Wang" 1 Yuxuan 2 Qingyuan Yang! |
Hu Ding'
Abstract
1. Introduction
To find a representative of several given probability distribu- tions is natural problem in machine learning. One popular approach is compute the geometric center on probability space with induced distances between probabilities, such as the seminal Wasserstein distance in the field of optimal trans- port (Villani, 2021). Given a non-negative weight vector (wi (W,, W.,......... WN) with 2 2, the Wasserstein baryeenter (WB) of N probability measures {ux k= under norm p is defined the weighted Fréchet Mean mean under Wasser- stein distance. Namely, it is the solution of the following problem
Wasserstein Barycenter (WB) is fundamental geometric optimization problem in machine learn- ing, whose objective is to find representative probability measure that minimizes the sum of Wasserstein distances to given distributions. WB has a number of applications in various areas. However, WB may lead to unfair outcome to- wards underrepresented groups in some applica- tions (e.g., "minority" distribution may be far away from the obtained WB under Wasserstein distance). To address this issue, we propose an alternative objective called "Wasserstein Ball Cen- ter (WBC)". Specifically, WBC is a distribution that encompasses all input distributions within the minimum Wasserstein distance, which be formulated as a "minmax" optimization prob- lem. We show that the WBC problem with fixed support is equivalent to solving a large-scale lin- ear programming (L.P instance, which is quite different from the previously studied LP model for WB. By incorporating some novel observa- tions on the induced normal equation, we propose efficient algorithm that accelerates the inte- rior point method by O(min{ Nm²² m²)) times N*'" is the number of distributions and "m" is the support size). Finally, we conduct a set of experiments on both synthetic and world datasets, demonstrating the computational effi- ciency of our algorithm, and showing its ability to provide more fairness for input distributions.
N
min (EP,())
Σ we W (4,
(1)
where (2))) is the set of Borel probability measure on Ω with finite p-th moment, and W, is the Wasserstein distance under norm p. which will be formally defined in Section 2. WB has found various applications in many fields, such as economics (Carlier & Ekeland, 2010; Chiappori, 2017), physics (Benamou al., 2014; Koehl et al., 2019), statis- ties (Kroshnin et al., 2021; Backhoff- Veraguas al., 2022; Goldfeld et al., 2024), medical imaging (Gramfort al., 2015; Nadeem et al., 2020), and machine learning (Qin et al., 2021; Cheng et al., 2021; Zhuang et al., 2022).
As the Fréchet mean under Wasserstein distance, WB tends to assign more measure to the region where the input density functions "cluster". other words, minimize the average distance from the barycenter to the input probabilities, if the support of most distribution is concentrated with high prob- ability in a region, then the WB should also have measure concentrated in that region. But this property may behave "unfairly" to "minority", i.e., the distributions with support deviated from the majority of others could be too far away from the WB. Fig. gives an intuitive demonstration for this issue.
'Equal contribution 'School of Computer Science and Technol- ogy. University of Science and Technology of China, Hefei, China "Scoool of Artificial Intelligence and Data Science, University of Science and Technology of China, Hefei, China. Correspondence to: Hu Ding <huding wwustc.edu.cn...
The unfaimess could cause negative impact in some scenar- ios (Mehrabi et al., 2021; Caton & Haas, 2024). In medical applications, datasets usually reflect societal biases or histor- ical inequities in healthcare access. If deep learning model is trained on biased data, it can propagate or even amplify these biases, leading to unequal treatment outcomes. For
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Finding Wasserstein Ball Center
Cloud
Cloud
Wasserstein distance
Cloud
Outtier
WB
WBC
WBC Outlier Cloud1 Cloud2 Cloud3
Figure 1. Left: four probability measures (cloud 1-3, and an outlier cloud), with their WB (computed by (1)) enclosed in purple ellipse, WBC (computed by (2)) enclosed in brown ellipse. The red cloud (i.e. the outlier cloud) has measure distributed distinctly from the other three clouds. In the histogram on the right, the "WBC" bin denotes the maximal Wasserstein distance from WBC to four clouds, and each of the other four bins denotes the Wasserstein distance from WB the corresponding cloud. From this histogram, can that WBC treats the outlier more equally (for example, the "WBC" bin is significantly lower than the "Outlier" bin).
instance, many skin lesion datasets predominantly contain images from lighter-skinned individuals (Bhardwaj & Rege, 2021; Adegun & Viriri, 2021). As a result, AI systems trained on such data often perform poorly on darker-skinned patients, leading to underdiagnosis or misdiagnosis cer- tain populations (Montoya et al., 2025). In addition, the medical data may come from different sources, therefore fair alignment also enhance model's lss accuracy (Aayush- man et al., 2024; Lin et al., 2023).
et al., 2021). In this framework, the data set for clustering consists of a set of groups, where each group contains set of data items. Standard clustering algorithms could incur higher clustering costs for certain protected groups (e.g., the groups that are defined by a sensitive attribute such as gen- der or race). This lack of fairness has motivated the study of clustering that minimizes the maximum clustering cost across different demographic groups, such as gender, race, high socioeconomic status, etc. The sensitive groups thus can obtain more equitable attention in the algorithm. Com- paring with the traditional clustering models (e.g., k-means or k-median clustering s), the main difference (which also the main challenge) lies this "minmax" objective.
To address this unfair issue, we propose a different objec- tive function. Rather than minimizing the summation of Wasserstein distances, we try find a distribution that of minimal distance from the farthest input distribution:
min max Wp(,, Mc). (D) ISESN
1.1. Our Main Contributions
Similar with the aforementioned socially fair clustering problems, solving Problem (2) also not easy due to its inherent "minmax" nature. When all distributions are of discrete support, the problem can be formulated as linear programming (LP) problem, where the details are shown in Section 2. Partly inspired by the recent developments of interior point method (IPM) and their various applica- tions (Gondzio, 2012) (e.g., the IPM algorithm for solving WB (Ge al., 2019)), we also consider developing an effi- cient IPM based algorithm for the WBC problem, though its formulation is much more complicated than WB due to the minimax issue.
From geometric perspective, can think of it the center of the ball in Wasserstein space, who covers all in- put distributions with minimum radius. In this setting, the output distribution does not put extra measure to the region where input distributions cluster with high density. Please see Fig. for an illustrative comparison. We call the solu- tion for Problem (2) the Wasserstein Ball Center (WBC), and aim design efficient algorithm to solve it. It should be noted that "Wasserstein ball" is not new concept and actually has been studied by several works before (Yue et al., 2022; Pesenti & Jaimungal, 2023; Chen et al., 2024), yet these previous works usually assume the ball center is given and take the ball as feasible region for constraining some optimization objective. In this paper, we focus on how to compute optimal center that the induced radius (under Wasserstein distance) is minimized.
Technically, there is a major challenge for directly applying IPM to the WBC problem, viz. the computational cost and space complexity are both very large. The linear program- ming formulation of WBC has m ΣN m, m + + variables and Nm + Σw m, + N + constraints, where the integer N denotes number of distributions, m; and m denote the support sizes for the i-th distribution and WBC respectively. This brings the challenge for computing the in-
The WBC model (2) also shares the similar "minmax" prop- erty with the recently proposed clustering problem under the notion of group fairness, called social fairness or eq- uitable group representation (Abbasi et al., 2021; Ghadiri
2

---

Finding Wasserstein Ball Center
ner loop of IPM, which requires a time complexity as large m4 + m++ N+)N)) To
problem (Ahuja et al., 1994). In the past years, a number of elegant algorithms have been proposed in theoretical com- puter science, such as (Orlin, 1988; Lee & Sidford, 2014; Sherman, 2017), and the ones focusing on Euclidean space (Khesin et al., 2021; Fox & Lu, 2023) or parallelism (L.ahn et al., 2023).
as O(NNm + Ei= m, + m,+ Σ== * 2 tackle this difficulty, we grind the intrinsic information of constraint matrix to simplify the Newton normal equation, which is a linear system with a large positive definite con- straint matrix, and is also the most expensive part in each inner loop of IPM. Specifically, we simplify the matrix in- verse occurred in the solution of Newton path, based on an important observation:
Wassertein barycenter. The computation for Wasserstein barycenter is an NP-hard problem (Altschuler & Boix- Adsera, 2022). Cuturi & Doucet (2014) showed that the computation for WB can be improved by using entropic regularization. Further progress includes iterative Breg- man projection (IBP) algorithm (Benamou et al., 2015), the semi-dual gradient descent (Cuturi & Peyré, 2018), accel- erated primal-dual gradient descent (APDAGI) (Kroshnin et al., 2019), alternating direction method of multipliers (ADMM) (Ye al., 2017), deterministic IBP (Lin et al., 2020), and the IPM algorithm MAAIPM (Ge et al., 2019). In Euclidean space, WB can be solved in polynomial time with fixed dimension (Altschuler & Boix-Adsera, 2021; Agarwal et al., 2025).
After performing some carefully designed transforma- tions, the seemingly dense coefficient matrix of normal equa- tion (a core part in our IPM framework) can be converted into a block-diagonal matrix in a particularly efficient way, where the inverse of each block is also easy to obtain.
Through this observation,
we obtain O(NmS 3 + +N2 + N time complexity for each iteration, whereas the vanilla IPM requires O(N*m*) by straight matrix inversion (for simplicity we just assume m, = =m) here). The latter one is often beyond accept- able scope in real-world scenarios. For example, assume we are given an instance that N 1000 and m 100. In each loop of IPM, the complexity of our algorithm is O(101)), while the vanilla IPM requires 0(10"7, which is O(10²) times higher. The strict description on this result is presented in Theorem 3.2. As for the practical effectiveness, our algorithm is significantly faster than the popular commercial solver "Gurobi" (Gurobi Optimization, LLC, 2024). which provides a concurrent optimizer that run different state-of-the-art optimization methods simultaneously. In particular for LP models, Gurobi simultaneouslyy runs the "dual simplex", "parallel barrier", and "primal simplex" algorithms in independent threads. In other words, our proposed algorithm runs faster than all these commonly used LP algorithms. In Section 4, we provide a set of experiments on numerical performance and applications.
Interior Point Method. Dikin (1967) proposed the first inte- rior point method. Further, Karmarkar (1984) proposed the first polynomial time method for linear programming called "Karmarkar's s algorithm". Mehrotra's 's predictor-correct algorithm (Mehrotra, 1992) provides the basis for most im- plementations of this class of methods, which is also the type of IPM applied in this paper (the details of predictor- corrector IPM are presented in Section 3.2). The seminal Mizuno-Todde Ye predictor method has quadratic convergence (Mizuno et al., 1993; Ye et al., 1993). For more details on IPM, we refer the reader to the survey pa- per (Gondzio, 2012). Recently, there are also several new studies on reducing the IPM complexity (Jiang et al., 2021; Cohen et al., 2021), which rely on technique called "matrix maintenance" to reduce the update time of each iteration.
2. Preliminaries
1.2. Related works
Some notations. Throughout this paper, we always use low case letter to denote scalar value, and bold low case letter to denote vector; also, we use "[n]" to denote the set {1,2, .n.. Let "1m" " be the all-one vector in RTM For two discrete probability vectors u E R and e E T with my and m2 being their support sizes, define the set of matrices M(u,v) = {II E R xm2 IIIm, IT 1 mm 4) as the coupling matrices, which con- sists of all joint distributions of margin u and v. Let @ = {(a;, q:) : E [mi]] denote the discrete probabil- ity measure supported on mi points qi...., qm, in Id with weights ai ai,..., am, respectively. Another probability mea- sure P (,), j E [m²]] is defined similarly. The p-norm Wasserstein distance of these two discrete proba-
We present an overview on several important related works, and the detailed introduction is placed in Section A.
Wasserstein distance Wasserstein distance is classic topic in machine learning (Rüschendorf, 1985). Cuturi (2013) introduced the "Sinkhorn Distance", which incorporates an entropic regularization term to smooth the transportation problem, offering paralleled and significantly faster solu- tions than exact computation of the discrete Wasserstein distance. Following Cuturi's work, recent years have seen the development of several improved Sinkhorn algorithms (Lin et al., 2019; Altschuler et al., 2019; Benamou et al., 2015; Altschuler et al., 2017). The discrete Wasserstein distance is also closely related to the min-cost max flow
3

---

Finding Wasserstein Ball Center
bility measures Q and P is
IPM framework introduced in Section 3.2 can run properly, In Section 3.3 we illustrate how to eliminate unnecessary computations to speed up IPM, as the key part of our con- tribution in the current article. We take advantage of some critical observation on the symmetry behind the seemingly dense constraint matrix. With that observation, we can sig- nificantly reduce the complexity of the most time consuming part in each loop of IPM, which is the inversing for large matrix.
my m2
W2(2,P)= == min 2
To ull:: -P)1BB)
1=1Ø=1
11 = (Ta)) (m² (m)) xm2)
E EM(a, (OT)
where a = (a1,-..,, )my) ) and b b= bmy) The coupling matrix achieving the minimum value can be seen as an alignment between the two distributions, which provides the optimal transportation (OT) plan between their supports.
3.1. From Minmax to LP
Now, consider the scenario where there are multiple prob- ability measures in the same space, {P(()) , = 1,E [N]}
First, we apply a slack variable y E I, turning Problem (4) into the following LP problem:
(N 2 2), where each p(t) {(a(0) qi (0) : E [m,]] with probability vector (t) (a] ..... am,) am))T The optimal Wasserstein ball center (WBC) Popt {(us, 2x)
(t)
(1)
min 7 w.II bY s.t. (w, IT() IT() E ES, 21, Im E R"
[m]} is another probability measure such that the maximum Wasserstein distance to these given N probability measures is minimized, as defined in the objective function (2) when Q = {z11.... }m The probability w = (W1,-.. 30m)" of Pope and its coupling matrices with {alo t [N]]
(5)
(DI S 1S1SN.
In (5), the value y indicates the maximal Wasserstein dis- tance from the WBC to the given N probability measures.
must be in the solution set S
Then, we conduct the following transformation on Prob- lem (5). We use ""vec(U"" denote the vectorization of matrix U, which is a vector obtained by concatenating ev- ery column of U to its predecessor, Let diag(U1, denote the block diagonal matrix with the i-th block U2, and "Ø" denote tensor product. Let "Imm be the identity matrix of order m. To transform the problem to a standard LP form, we vectorize the constraint (ii) (3), "I[(") I.m. w" and "(IT()) 1m a11), to be:
{ww, IT(1) (..,,...I((A)))
ITM) E RT x R m x x R mxmm :
JT = 1, 20: 0:
Umm)
constraint
(T(0) 1 w. (TI(T)
1 a (0)
(7() 2 e [NT),
(3)
constraint
(ii)
where "z" means each entry of the vector/matrix is non- negative, the constraint (i) guarantees w to be a feasi- ble probability distribution, and the constraint (ii) fol- lows the coupling condition in (OT). For given sup- port 2 the distance matrices are defined as (12) (||2;- 9,"" [B)((.)) E mm for t = I,...,N. .N. Since in most practical applications, the measures in ( have finite support points, and we can assume the barycenter is supported on given finite set Q (in general, Ω does not need to be the same as the supports of p())s), the problem WBC in (2) can be re- re-written as:
(1 8l (mc(((( = w. (Jm, 11 II vec = a
for = 1...., N. We further define the following ma- trices: Fn = diag( Im @ m. .... Imx 8 1 m). Fg = diag(1 m; & lm...., 1 1 & 4m), E3 = 11 8 Im and D (Daa(v))((()))....)))))))))))))) ..... (..())))......((())))
Through some simple calculations, Problem (5) can be formulated standard LP model:
min c
T a s.t. Ax = b, x 2 0
(WBC-LP)
min max min (D( ) s.t. (w,ITO) ITO), ..... IT(N) E ES AE[N] 1199)
where x (vec IT (1) ...: vec [T (N))
w 71 YN. 7), % 7 (D(i) II (1)) b (a (1) ...aa 0 ... Om 1),
(4) where "(,,')" denotes the inner product in Euclidean space, and ((t) denotes D' (I2 for simplicity.
c = (0: o:...... ; 0: 1) and
E
3. Our Optimization Framework for WBC
E2
Eg
1
In this section, we introduce our optimization framework for solving the Problem (4). Specifically, in Section 3.1 we formalize Problem (4) to be the standard LP form, ensur- ing that the constraint matrix is full row-rank so that the
D
IN - -1
Let M := E m, and then the numbers of constraints
and variables are ne = Nm + M + N + and ny ==
4

---

Finding Wasserstein Ball Center
Mm + N + 1, respectively. For visual demonstration of these matrices, please see Eq. (12).
Because of the quadratic term "AXAs the right hand- side (RHS), Eq. (6) does not have explicit solution. To reduce this nonlinear system to be linear case, first we obtain a predictor step by removing the " "AXAs+ In." term on the RHS of Eq. (6), then compute the corrector step by assigning the solution of predictor step to the quadratic term in the RHS of Eq. (6). Through combining them to- gether, we obtain the final update step, hence the name "predditor-cooreetoor.. The value /4 is updated after each loop, descending as the primal-dual gap narrows. When p4 approaches to O, the current position convergences to the optimal solution of WBC-LP.
Preprocessing on (WBC-LP). To implement IPM, it is essential that the constraint matrix A should be of full row- rank, so that the problem's 's complexity is reduced and IPM is also guaranteed to work more efficiently. In particular, the "full row-rank" property can simplify the procedure for solving the "normal equation" in IPM (which is discussed in next section). The following lemma eliminates all re- dundant constraints, turning A into a full row-rank matrix A. Specifically,, A E (g(".-))x,, is the matrix obtained from A by removing the (M + 11)thh (M + m + 1)-th. ... (M + (N- - 1)m= + 1) rows of A, and E R= is the vector obtained from b by removing the (M + 1)-th, (M + m 1)-th ..... (M + (N 1)m++ 1) entries of b. Proposition 3.1. 1) A has full row-rank; 2) solving the equation Ax b in (WBC-LP) is equivalent to solving the equation Ax= = b.
N
Built upon the above idea, both the predictor and corrector steps can be obtained by sequentially computing Ay. As and Ax, where the complete details are shown in Algo- rithm 2 of Appendix C. Here, we only focus the most time consuming part, computing Ay. as the solution of normal equation (Wright, 1997). Namely, for both the predictor and corrector steps, we need to solve system of linear equations as the following form:
We leave the proof to Section B. Due to Theorem 3.1, we will focus on A instead of A in the following sections.
(ARA Ay
3.2. Predictor-Corrector IPM on WBC-LP
where R = diag(s)- -1 X, and the symbol "f" represents the vector obtained in either the predictor or corrector steps, as shown in Eqs. (13) and (14) of Algorithm 2 (here, we can temporarily ignore the detailed form of f since it does not affect designing our algorithm for solving Eq. (7)).
We choose the classic predictor cooereccoororrrooroooooooo scheme (Mehrotra, 1992), which a variant of primal-dual IPM, that applies the block coordinate descent on the Lagrangian of the linear program smoothed by logarithm regularization. As a second order method, it is proved to have quadratic convergence rate (Ye et al., 1993), which surpasses first-order methods in general (for a detailed analysis, please see (Wright, 1997)).
Other parts in one loop of Algorithm are just simple oper- ations like matrix-vector multiplications, and thus the whole time complexity is dominated by the cost for solving Eq. (7), We can solve Eq. (7) by inverting the matrix applied to Ay (this requires Á to be of full row-rank, and this is also the reason why we need Theorem 3.1). However, ordinary inverse to such a large matrix will take excessively high complexity, as shown in Remark 3.3 later. This motivates us to propose new method to further improve the complexity in next section.
The high level idea. When we deal with a primal-dual sys- tem of linear programming, from the Marssh-KunnTTucke theory (Wright, 1997), we can update the primal and dual variables step by step using Newton-like method to the Lagrangian of LP. In particular, the algorithm takes loga- rithm barrier function, which is determined by a coefficient 14 2 0, to measure how close the current solution is to that of the original problem described in WBC-LP.
3.3. Solving the Normal Equation (7) Efficiently
Let y E R and 8 E RR'V be the dual variables cor- responding the constraints "Ax 6" and "x 2 o", respectively. Writing in matrix form, the search direction at the current primal-dual tuple (a, y, 8) should be the solution of the following nonlinear system of equations:
We introduce efficient algorithm for solving the normal equation (ARAT) A (ARAT)Ay= f of Eq. (7), whose complexity is summarized in the following theorem.
Theorem 3.2. In each inner iteration of IPM, the time complexity in terms of flops is O(m² E my + ~Nm3 + N2m² + N3), and the memory usage in terms of doubles O(m Σn= m; + Nm= + N2).
To
AT
In.
Ax Ay
A y t A Ax - b
A
0
0
(6)
S
0
X
As
Xs+ AxAs 1
Remark 3.3. For comparison, a vanilla IPM takes O(NN+ Σ m²+N)² EN m,+++NN) time and O(NNm+ m4 + (m ΣN m, + + N)) storage.
where the notation "A" denotes the change when updating some variable or matrix. For the vector æ, in Eq. (6) we use the capitalized letter "X" to denote the diagonal matrix with X = I₁, where H; is the i-th term of æ. Similarly, we also have the diagonal matrix S for the vector 8 in (6).
Roadmap of the proof for Theorem 3.2. It is worth noting that we do not need to really compute the matrix (ARÄT)=
5

---

Finding Wasserstein Ball Center
for solving Eq. (7); instead, we only need find efficient way simplify the structure of Á ARAT, so that can easily compute the matrix-vector multiplication "( "(ARAT) 11 If". Proposition 3.4 illustrates the detailed structure of ARA", which paves the way for designing the following transforma- tions. Our method contains three steps. First, we transform ARA to be a nearly block-diagonal matrix. Then, we re- solve the most challenging part that hinders it from being block-diagonal, as shown in Theorem 3.6. Finally, ob- tain the desired block-diagonal form for ARAT, Namely, to compute the product (Á (-RAA-1A we only need to care about the inverses of several small matrices with sizes no bigger than max ma m, N). As consequence, we save the massive time and storage for inverting a matrix of size O + N max(m:)).
several transformations. Relying Proposition 3.4 we de-
I
sign the transformation matrix Q ==
1
IN
IM TM
1s By
-axt
BT
such that
In
B
G1 + G2
K₂
Q. DARA QQ=
(8)
K
W
where G₁ := B3 - By By B2, G2 == B₁ t T, W
W K7 B K1. See Section E for computation.
Step 2: dealing with (G1 + G2) the RHS of Eq. (8), we can see that the transformed matrix is still not block- diagonal due to the existence of "Kz". To further eliminate K2, we need to design another transformation matrix that involves the computation of (G₁ + G2) -1, the inverse of a large matrix of size M' x M'. To find a way to compute (G1 + G2) more efficiently, we have a useful observation on the matrices G1 and G2, which indicates the specificity of their structures (the proof of Proposition 3.5 is placed in Section F).
Let be the nn-diimnsioonll vector with its i-th entry Ra (recall that nh denotes the number of vari- ables as defined in Section 3.1). Also recall that M = E m,, and we further define M' = N(m - 1). Let = r(Mm- + 2: Mm + m). the sub-vector of r includ- ing from the ( (Mm- + 2)-th to (Mm + mh) entrie.. Proposition 3.4. ÄRÄT has the following form:
B
B2
0
Ki
ARA--
B2
B3
+
B1
K₂
Proposition 3.5. There exist positive definite matrices Y, Gu [Im-1)x(m-- i L...., N, such that G1 diag(Gn1.... ,NNN) and G2 = (1x15)@1 &Y
0
a T
1
KT
K.T
0
W
where the sub-man mabricasices are defined as follows (their illus- trations are shown in Section D):
With the above observation, we achieve the following key lemma.
Lemma 3.6. Let 7 diag(z) 122 (z is defined before
.Bi E RRMM is a diagonal matrix with positive di- agonal entries; B2 E Mxx' is a blagk-diagonal matrix with N blocks, where each i-th block is of size my x (mm1))) Bg E RRM xM' is diagonal matrix with positive diagonal entries; B4 = (1)1 & diag(z);
Theorem 3.4), and G = (1v1 (v) s (2-1 - G then we have
(9)
Further, the flops needed applying a vector to the RHS of Eq. (9) is O(Nm²).
a = -I & 2; 1 1al1-I][(mm [_ m + nn :N);
Ki E M is dlock-diagonal matrix with N
blocks, with the i-th block being of size my x 1;
Proof. Since Y is positive definite, there should exit ma- trix U E R(m--)x(m-1), such that Y = UTU. Then, from Proposition 3.5 we have G2 = (1v1w) @ Y (IN U (I] @ U). and it implies
K2 E MNN is a block-diagonal matrix with N blocks, with the i-th block being of size (m - 1) x 1;
W W1 1n,1 N. where Tn, is the n term of r, and W1 E NNN is a diagonal matrix with positive diagonal entries.
(G1 + G2) - 2 (OTGHH(NTT))
= = - GT (1N
(mm-1+(1) & UUG @ U (1) & U) GT!.
7
We left the proof for Proposition 3.4 to Appendix D. Then, we illustrate our idea for efficiently solving the normal equa- tion (7) via three steps.
where the second equality is obtained by Woodbury ma- trix identity (Hager, 1989), a commonly used equation in linear algebra (see Section H). Noticing that both G1 and I s U are divided into N N sub-matrices of size
Step 1: transforming ARA to be nearly block-diagonal. We simplify the coefficient matrix ARA by performing
6

---

Finding Wasserstein Ball Center
200
15000
WBC
WBC
150
Gurobi
object value
Gurobi
10000
time
time
100
5000
50
0
20
30
40
50
60
70
80
90
20
30
40
50
60
70
80
90
1000
2000
3000
4000
5000
6000
N
N
N
(
1200
WBC
WBC
10000
1000
Gurobi
object value
Gurobi
8000
800
time
time
6000
4000
20
2000
100
200 300 400 500
600
700
100
200
300
400
500
1000
2000
3000
4000
5000
6000
m
m
m
(b)
(d)
(f)
Figure The first two column figures are the computational time and obtained objective value Gurobi and our method. For (a) and (c), we fix m 100 and vary N. For (b) and (d), we fix N = 30 and vary m. The third column figures are the computational times of our method when the problem scale very large. (e), m 50. For (f), we N
(m-11) - 1) x (m - 1). together with Proposition 3.5, we can simplify the RHS of the above equality to
the first m- rows, then copy them N times down to obtain the multiplication. This part also takes O time.
- ((wID))
Step from nearly block block-diagonal to block-diagonal, and achieve the final solution Ay. We design another transformation matrix Q' to eliminate K₂ in the nearly block diagonal matrix obtained in (8). Specifically, we let Q' IM
(U-
LUG
N - ( Σ
1x
Through combining with the
G G
1
KJJ (G++G2) transformation Q introduced in (8), we have
IN
i=1
Therefore, obtain the equality (9).
Q' Q'QQRATO QT OT = diag B1, G1 + G2, W), (10)
As for the complexity, we notice the following two facts. (1) GT is a block diagonal matrix with only N(m - 1)2 non-zero terms, thus the operation multiplying vector to it takes no more than O time. (2) The matrix ă can be viewed being (Z-1+ - N - That is, only have to store (Z-I+EN G )-1 square matrix of size m instead of G with size N (m- 1). Immediately, if we multiply vector on the right of G. the i-th coordinate of the resulting vector would be the same as the (i+m- -1)-th coordinate for any < M' (m- - 1). Therefore we only have to compute
+ G2) -22 Based on Lemma 3.6, the construction of W can be done in m flops, since K₂ has N columns.
where W W W- K (G1 (
N N matrix with every term
Now, we can solve the normal equation (ÄRA A = of Eq. (7). Set ÿ as the solution to
diag Bi Gi + G2, W) = Q'QS,
(11)
where the left coefficient matrix from the obtained block- diagonal matrix (10). Similar with the method for comput- ing W, due the term "G1 + G2 in the coefficient matrix

---

Finding Wasserstein Ball Center
1200
WWCC - WWB
log(value)
Maximum
Max WD
800
630
000
200
WBC
Middle
Fastem
While
Southoul
Asian
Indian
Latino-
Hiquaiic
East
Avians
Black
0.2
0.4
0.6
0.8
Iteration (a)
Imbalance factor
Race
(b)
Figure 3. (a) N 90, m 2000. The objective value of our algorithm becomes stable after 80 iterations. For other settings and N. see Section J. (b) The maximum Wasserstein distance between the barycenter (WBC or WB) and input distributions, with fixing N 20, m = 500. (c) The dark bin marked "WBC" is maximum WD between WBC and all the faces, while the rest represent the maximum WD from WB to the faces in different races.
of y. we also need to apply Theorem 3.6 to obtain y. After that, from Eq. (10) it is easy to verify that Ay = T " y.
usage over the commercial solver Gurobi, a powerful opti- mization solver widely used across various fields such as operations research, finance, and data science. The baseline choose is Gurobi Optimizer version 11.0.0 (academic license). (2) The second experiment reflects the fairness of WBC over the standard WB. For these two experiments, we generate random datasets in Euclidean space, and the weights of (gi", .-- Ti (t) in each distribution p(() gener- ated uniformly at random.
Algorithm summarizes the above steps. Due to the space limit, we leave the detailed proof for the overall complexities Section G.
Algorithm Solver for (ARAT) A
Input: R X", ff E GR"-M as described in Eq. (7)
(3) The third experiment further illustrates the performance on world dataset FairFace (Karkkainen & Joo, 2021) with diverse racial representa- tion. We choose 700 (100 for each race) images includ- ing seven racial groups of "Black", "East Asian", "Indian", LLatinoHHispanie "Middle Eastern', "Southeast Asian" and "White" as 700 distributions, where each image actually can be regarded as a distribution of size 50 x 50 in R2. (4) The fourth experiment on 3D point-cloud gives a visual example to the fairness of WBC as a barycenter. All the experiments are implemented on a workstation, Intel(R) Core(TM) i5- 9400 CPU @ 2.90GHz and 8GB for RAM, equipped with win64 Windows 11+.0.
Output:
The solution Ay
Step 1: Compute blocks of ARÁ
1: Compute B1, B2, B3, K1, K2, W
Step 2: Simplify the constraint matrix into block diagonal
1: Compute Q
2: K2 f- K₂ B,J B Ki
3: W + W- K2 Bi K1
//Eliminate K1
4: Compute G₁, G2
5: Decompose (G1 + G2) according Lemma 3.6
6: Compute Q'
7: W + W K2 (G1 + -2))-KK,
//Eliminate K₂
Step 3: Solve Ay
1: æ1 + Qf. z2 t Q'211
(1) Comparison with Gurobi. We of all distributions to be equal for brevity. As shown in Fig. 2, our algorithm is always faster than Gurobi, and the gap between the two methods is expanding as the scale increases. Moreover, Gurobi cannot solve the instance with m > 500 due to memory limitation, which showcases the advantage of our space complexity in Theorem 3.2. We also illustrate the convergence speed of our algorithm. From Fig. 3(a) (and Figs. 5 to 8 in appendix), we see that algorithm converges in a super- linear rate, which is consistent with the theoretical result of (Ye et al., 1993).
2: 23(1 M) <- BT B,'1z(11 "2(11 M)
zn(n. N ::ee) t W 22 +1:100 ne)
4: 21(n, -N) -1 (n. N)
5: Solve (G1 + (z) (M +l:no - - =2(M++1 +1:
N 1)
6: +- Q"T zạ, Ay +- Q 23
Return: Ay
4. Experiments
(2) Fairness of WBC. To evaluate the performance on fair- ness, we first compare the maximum Wasserstein distance (Max WD) between the obtained centers and input distri-
We conduct four experiments to investigate the real perfor- mance of our algorithm. (1) The first experiment demon- strates our advantages on computational speed and memory
8

---

Finding Wasserstein Ball Center
Figure 4. WBC (left) and WB (middle) of 3 cars and one chair. WB looks just like car while our proposed WBC preserves characteristics chair, such as backrest legs, indicating that WBC indeed places greater emphasis on fairness the minority
butions. The experimental distributions are partitioned into two families, with each family's 's measure concentrated in distinet regions that are quite distant from each other. To quantify the disparity between these two families, we define the imbalanced factor ("imf"" as the ratio of the smaller fam- ily's size to the larger one. For instance, "imf I" indicates that both families contain equal number of distributions. As is shown in Fig. 3(b), the smaller the imf, the bigger the difference of the Max WD between WB and WBC (except for the case "imf=0", since the smaller family empty this case), indicating WBC's potential as the alternative to WB when fairness is concerned. Please Section J.2 for the details of this part.
which is far greater than the object value of our algorithm (45.37). Mcanwhile, our algorithm achieves the object value only slightly higher than the average of WD to WB (38.70).
(4) 3D Point-cloud Averaging. WBC provides means for computing the average shape of point-cloud data by casting point cloud into an empirical probability measure. We take the ShapeNet Core-55 dataset (Chang et al., 2015), with each cloud containing 2048 points in R3 As an example, in Fig. we show that compared with WB, the WBC preserves more characteristics of the "chair" (which can regarded as an outlier distribution from the other "car" clouds).
5. Conclusion
Another fairness measurement is the variance of Wasserstein distances from the obtained center to the input distributions: (W( pi)) = - Σ (F. W denoted as var(WDD, where W Σ (p.,),) Smaller var(WDD indicates reduced representative disparity among all inputs. Across all 2,700 experiments (detailed in Section J.2), WBC consistently achieves lower var(WD) compared to WB. For example, in L1 space of dimension 100, we sample the 30 supports uniformly in a cube with side length 3, then trans- late 3 supports by adding - 100. and translate 5 supports by adding 1 100. we have var(WD)=58.35 for WB, while var(WD)=411866 for WBC. Those observations suggest that WBC not only improves fairness for the minority distribu- tions, but also could enhance fairness across all distributions by limiting disparities between transport plans.
In this paper, we present an efficient algorithm to com- pute Wasserstein ball center, an alternative of Wasserstein barycenter with emphasis on minority distributions, which makes it more suitable for the tasks that are sensitive faimess. For future work, it would be interesting to explore further properties and applications of WBC. Moreover, it is also deserved to study other models for fair ensembling different distributions in the field of machine learning.
Acknowledgements
We want thank the anonymous reviewers for their con- structive comments. This work was supported in part by the National Key Research and Development Program of China FFA0000000,,0000,,, the National Natural Sci- ence Foundation of China (NO.62272432, NO.62432016) and the Natural Science Foundation of Anhui Province NN0.228085MFF66).
(3) Experiments on FairFace Dataset. From Fig. 3(c). we can observe that the standard WB has significant gap between the object values different races. The object value of WB corresponding to "Middle Eastern" is 78.30,

---

Finding Wasserstein Ball Center
Impact Statement
Altschuler, M. and Boix-Adsera, E. Wasserstein barycen- ters np-hard to compute. SIAM Journal on Mathemat- ics of Data Science, 4(1):179-205) 2022.
The proposed algorithm provides an effective alternative to the Wasserstein barycenter by offering a trade-off between transport efficiency and fairness. However, practitioners should be cautious of potential outliers arising from noise or data emros. To ensure robust performance, data cleansing or selection techniques should be applied prior to implement- ing the WBC.
Backhoff-Veraguaa J., Fontbona, J., Rios, G., and Tobar, F. Bayesian learning with wasserstein barycenters. ESAIM: Probability and Statistics, 26:436 472366472,440 2022.
Backurs, A., Indyk, P., Onak, K., Schieber, B., Vakilian, A., and Wagner, T. Scalable fair clustering. In International Conference Machine Learning, pp- 405-413. PMLR, 2019.
References
Aayushman, Gaddey, H., Mittal, V., Chawla, M., and Gupta, G. R. Fair and accurate skin disease image classifica- tion by alignment with clinical labels. International Conference Medical Image Computing and Computer- Assisted Intervention, pp. 394-404. Springer, 2024.
Benamou, J.-D., Froese, D., and Oberman, A. M. Numer- ical solution of the optimal transportation problem using the monge-ampère ampère equation. Journal of Computational Physics, 260:107-126, 2014.
Benamou, J.-D., Carlier, G., Cuturi, M., Nenna, L., and Peyré, G. Iterative bregman projections for regularized transportation problems. SIAM Journal Scientific Com- puting, 37())A1111-A1133 2015.
Abbasi, M., Bhaskara, A., and Venkatasubramannia, S. Fair clustering via equitable group representations. In Pro- ceedings of the 2021 ACM conference on fairness, ac- countability, and Iransparency, pp. 504-514, 2021.
Bhardwaj, A. and Rege, P. P. Skin lesion classification using deep learning. In Advances in Signal and Data Processing: Select Proceedings of ICSDP 2019, pp. 575- 589. Springer. 2021.
Adegun, A. and Virini, S. Deep learning techniques for skin lesion analysis and melanoma cancer detection: survey of sfate-of-the-art. Artificial Intelligence Review, 54(2): 811-841, 2021.
Borgwardt, S.
An Ip-based, strongly-polyooee 2- approximation algorithm for sparse wasserstein barycen- ters. Operational Research, 22(2):1511-155) 2022.
Agarwal, A., Dudík, M., and Wu, S. Fair regression: Quantitative definitions and reduction-based algorithms. In International Conference Machine Learning, pp- 120-129. PMLR, 2019.
Carlier, G. and Ekeland, T. Matching for teams. Economic theory, 42:397-418, 2010.
Agarwal, P. K., Raghvendra, S., Shirzadian, P., and Yao, K. Efficient approximation algorithm for computing wasser- stein barycenter under euclidean metric. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 4809 4826. SIAM, 2025.
Caton, S. and Haas, C. Fairness in machine learning: survey. ACM Computing Surveys, 56(7):1-38, 2024.
Chai, J. and Wang, X. Fairness with adaptive weights. In International Conference on Machine Learning, pp. 2853- 2866. PMLR, 2022.
Ahuja, R. K., Orlin, J. B., Stein, C., and Tarjan, R. E. Im- proved algorithms for bipartite network flow. SIAM Jour- nal on Computing, 23(5):906-9333 1994.
Chang, A. X., Funkhouser, T., Guibas, L., Hanrahan, P., Huang, Q., Li. Z., Savarese, S., Savva, M., Song, S., Su, H., al. Shapenet: information-rich 3d model repository. arXiv preprint a 2015.
Altschuler, J., Niles-Weed, J., and Rigollet, P. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. Advances in neural information pro- cessing systems, 30, 2017.
Chen, X., Fain, B., Lyu, L., and Munagala, K. Proportionally fair clustering. International Conference on Machine Learning, pp. 1032-1041. PMLR, 2019.
Altschuler, J., Bach, F., Rudi, A., and Niles-Weed. J. Mas- sively scalable sinkhorn distances via the nyström method. Advances in neural information processing systems, 32, 2019.
Chen, Z., Kuhn, D., and Wiesemann, W. Data-driven chance constrained programs over wasserstein balls. Operations Research, 72(1) ((1)):10-424, 2024.
Altschuler, J. M. and Aoix-Adsera, E. Wasserstein barycen- ters can be computed in polynomial time in fixed dimen- sion. Journal of Machine Learning Research, 22(44): 1-19, 2021.
Cheng, K., Aeron, S., Hughes, M. C., and Miller, E, L. Dy- namical wasserstein barycenters for time-series modeling. Advances in Neural Information Processing Systems, 34: 27991-28003, 2021.
10

---

Finding Wasserstein Ball Center
Chiappori, P.-.. Matching with transfers: The economics of love and marriage. Princeton University Press, 2017.
Gong, J. and Kim, H. Rhsboost: Improving classification performance in imbalance data. Computational Statistics & Data Analysis, 111:1-13, 2017.
Chierichetti, F,, Kumar, R., Lattanzi, S., and Vassilvitskii, S. Fair clustering through fairlets. Advances in neural information processing systems, 30, 2017.
Gramfort, A., Peyré, G., and Cuturi, M. Fast optimal trans- port averaging of neuroimaging data. In Information Processing in Medical Imaging: 24th International Con- ference, IPMI 2015, Sabhal Mor Ostaig. Isle of Skye, UK, June 28-July 3, 2015, Proceedings 24, pp. 261-272. Springer, 2015.
Cohen, M. B., Lee, Y. T., and Song. Z. Solving linear pro- grams the current matrix multiplication time. Journal of the ACM (JACM), 68(1):1-39, 2021.
Cuturi, M. Sinkhorn distances: Lightspeed computation of optimal transport. Advances neural information processing systems, 26, 2013.
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2024. URL ht ht /www://wwwww//w///////////wwo... gurobi......... com.
Hager, W. W. Updating the inverse of matrix. SIAM
Cuturi, M. and Doucet, A. Fast computation of wasserstein barycenters. International conference on machine learning, pp. 683-693. PMLR, 2014.
review, 31(2):221 1989.
Huang, M., Ma. S., and Lai, Projection robust wasserstein barycenters. In International Conference Machine Learning, pp. 4456-4465. 2021.
Cuturi, M. and Peyré, G. Semidual regularized optimal transport. SIAM Review, 60(4):941-96 2018
Jang, T., Zheng, F., and Wang, X. Constructing a fair clas- sifier with generated fair data. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pp. 1908-7916, 2021.
Dikin, I. Iterative solution of problems of linear and quadratic programming. In Doklady Akademii Nauk, vol- ume 174, pp. 747-748. Russian Academy of Sciences, 1967.
Jiang, H. and Nachum, O. Identifying and correcting label bias in machine learning. In International conference on artificial intelligence and statistics, pp. 702-712 PMLR, 2020.
Dognin, P., Melnyk, I., Mroueh, Y., Ross, J., Santos, C. D., and Sercu, T. Wasserstein barycenter model ensembling. In International Conference on Learning Representations, 2019.
Dwork, C., Hardt, M., Pitassi, T., Reingold, O., and Zemel, R. Fairness through awareness. In Proceedings of the 3rd innovations in theoretical computer science conference, pp. 216-226, 2012.
Jiang, S., Song. Z., Weinstein, O., and Zhang, H. A faster algorithm for solving general lps. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 923-832 2021.
Fox. E. and Lu, J. A deterministic near- near-linee linear time approx- imation scheme for geometric transportation. In 2023 IEEE 64th Annual Symposium on Foundations of Com- puter Science (FOCS), pp. 1301-1315. IEEE, 2023.
Joseph, M., Kearns, M., Morgenstern, J. H., and Roth, A. Fairnsss in learning: Classic and contextual bandits. Advances in neural information processing systems, 29, 2016.
Ge. D., Wang, H., Xiong, Z., and Ye, Y. Interior-point methods strike back: Solving the wasserstein barycenter problem. Advances in neural information processing systems, 32, 2019.
Karkkainen, K. and Joo, J. Fairface: Face attribute dataset for balanced race, gender, and age for bias measurement and mitigation. In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pp. 1548- 1558, 2021.
Ghadiri, M., Samadi, S., and Vempala, S. Socially fair k- means clustering. In Proceedings of the Conference on Fairness, Accountability, and Transparency, pp. 438 448, 2021.
Karmarkar, N. A new polynomial-tine algorithm for linear programming. In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pp. 302-311, 1984.
Goldfeld, Z., Kato. K., Rioux, G., and Sadhu, R. Statistical inference with regularized optimal transport. Informa- tion and Inference: A Journal of the IMA, 13(adiaa, 2024.
Khesin, A. B., Nikolov, A., and Paramonov, D. Precondi- tioning for the geometric transportation problem. Journal of Computational Geomelry, 1 11(2):234-22 2021.
Gondzio, J. Interior point methods 25 years later. European Journal of Operational Research, 21863),587-601, 2012.
Kim, K., Zhou, B., Zhu, C., and Chen, X. Sobolev gradient ascent for optimal transport: Barycenter optimization and

---

Finding Wasserstein Ball Center
convergence analysis. arXiv preprini arXiv:2505.13660, 2025.
Makarychev, Y. and Vakilian, A. Approximation algorithms for socially fair clustering. In Conference on Learning Theory, pp. 3246-3264. PMLR, 2021.
Kochl, P., Delarue, M., and Orland, H. Statistical physies approach to the optimal transport problem. Physical review letters, 123(4):040603, 2019.
Martinez, N., Bertran, M., and Sapiro, G. Minimax pareto
fairness: A multi objective perspective. In International conference on machine learning, pp. 6755-6764. PMLR, 2020.
Kolouri, S., Nadjahi, K., Simsekli, U., Badeau, R., and Rohde, G. Generalized sliced wasserstein distances, Ad- vances in neural information processing systems, 32, 2019.
Mehrabi, N., Morstatter, F., Saxena, N., Lerman, K., and Galstyan, A. A survey bias and fairness in machine learning. ACM computing surveys (CSUR), 54(6))1-33, 2021.
Krasanakis, E., prrromitros-Xioufis E., Papadopoulos, S., and Kompatsiaris, Y. Adaptive sensitive reweighting to mitigate bias in fairres--wware classification. In Proceed- ings of the 2018 world wide web conference, pp. 853-862, 2018.
Mehrotra, S. On the implementation of a primal-dual inte- rior point method. SIAM Journal on optimization, 2(4): 575-601, 1992.
Mizuno, S., Todd, M. J., and Ye, Y. On adaptive-step primal-dual interior-point algorithms for linear program- ming. Mathematics of Operations research, 18(4):964- 981, 1993.
Kroshnin, A., Tupitsa, N., Dvinskikh, D., Dvurechensky, P., Gasnikov, A., and Uribe, C. On the complexity of approximating wasserstein barycenters. In International conference on machine learning, pp. 3530-3540. PMLR, 2019.
Montoya, L. N., Roberts, S., and Hidalgo, B. S. To- wards fairness in ai for melanoma detection: Systemic review and recommendations. In Future of Information and Communication Conference, pp. 320-341. Springer, 2025.
Kroshnin, A., Spokoiny, V., and Suvorikova, A. Statistical inference for wasse-waser barycenters. The Annuls of Applied Probability, 3 31(3):1264-123 2021.
Lahn, N., Raghvendra, S., and Zhang, K. A combinato- rial algorithm for approximating the optimal transport in the parallel and mpc settings. Advances in Neural Information Processing Systems, 36:21675-21686, 2023.
Nadeem, S., Hollmann, T., and Tannenbaum, A. Multi- marginal wasserstein barycenter for stain normalization and augmentation. In Medical Image Computing and Computer Assisted Intervention-MIICAI 2020: 23rd In- ternational Conference, Lima, Peru, October 4-8, 2020, Proceedings, Part V 23, pp. 362-371. Springer, 2020.
Lee, Y. T. and Sidford, A. Path finding methods for lin- ear programming: Solving linear programs in o (vrank) iterations and faster algorithms for maximum flow. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pp. 424 433. IEEE, 2014.
Nietert, S., Goldfeld, Z., Sadhu, R., and Kato, K. Statisti- cal, robustness, and computational guarantees for sliced wasserstein distances. Advances in Neural Information 2022.
Processing Systems,
Lin, L., Shi, W., Ye, J., and Li, J. Multisource single-cell data integration by maw barycenter for gaussian mixture models. Biometrics, 79(2):866-877 2023.
Orlin, J. A faster strongly polynomial minimum cost flow algorithm. In Proceedings of the Twentieth annual ACM symposium Theory of Computing, pp. 377-387, 1988.
Lin, T., Ho, N., and Jordan, M. On efficient optimal trans- port: An analysis of greedy and accelerated mirror de- scent algorithms. In International Conference Ma- chine Learning, pp. 3982-3991, 3MM92-..1..MMRR
Pele, O. and Werman, M. Fast and robust earth mover's distances. In 2009 IEEE 12th international conference on computer vision, pp. 460-467. IEEE, 2009.
PMLR, 2019.
Lin, T., Ho, N., Chen, X., Cuturi, M., and Jordan, M. Fixed- support wasserstein barycenters: Computational hardness and fast algorithm. Advances in neural information pro- cessing systems, 33:5368-523 -333680-0.. 2020.
Pesenti, S. M. and Jaimungal, S. Portfolio optimization within wasserstein ball. SIAM Journal on Financial Mathematics, 14(4): 2023.
Qin, R., Li, M., and Ding, H. Solving soft clustering en- semble via k-sparse discrete wasserstein barycenter. Ad- vances in Neural Information Processing Systems, 34: 900-913, 2021.
Lohaus, M., Perrot, M., and Von Luxburg, U. Too relaxed to be fair. In International Conference on Machine Learning, pp. 1360-6366 PMLR, 2020.
12

---

Finding Wasserstein Ball Center
Rüschendorf, L. The wasserstein distance and approxima- tion theorems. Probability Theory and Related Fields, 70 (1):117-129, 1985.
Sherman, J. Generalized preconditioning and undirected minimum-cost flow. In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 772-780. SIAM, 2017.
Song, S., Mo, G., and Ding, H. Relax and merge: sim- ple yet effective framework for solving fair k-means and k-sparse wasserstein barycenter problems. In The Thir- teenth International Conference Learning Representa- tions, 2025.
Vaskevicius, T. and Chizat, L. Computational guarantees for doubly entropic wasserstein barycenters. Advances in Neural Information Processing Systems, 36, 2024.
Villani, C. Topics in optimal transportation, volume 58.
American Mathematical Soc., 2021.
Wang, F., Cheng, L., Guo, R., Liu, K., and Yu, P. Equal opportunity of coverage in fair regression. Advances in Neural Information Processing Systems, 36, 2024.
Wright, S. Primal-dual interior-point methods. SIAM, 1997.
Yang, Q. and Ding, H. Approximate algorithms for k-sparse wasserstein barycenter with outliers. In Proceedings of the Thirty-Third International Joint Conference Artifi- cial Intelligence, pp. 2016--335, 2024.
Ye, J., Wu, P., Wang, J. Z., and Li, Fast discrete distribu- tion clustering using wasserstein barycenter with sparse support. IEEE Transactions on Signal Processing, 65(9): 2317-2332, 2017.
Ye. Y., Güler, O., Tapia, A., and Zhang, Y. quadrati- cally convergent (I)-iteration algorithm for linear pro- gramming. Mathematical programming, 19(1:151-112, 1993.
Yue, M.-C., Kuhn, D., and Wiesemann, W. On linear opti- mization wasserstein balls. Mathematical Program- ming, 195 2022.
Zhuang, Y., Chen, X., and Yang. Y. Wasserstein k-means for clustering probability distributions. Advances in Neural Information Processing Systems, 35:1 1382-11395 2022.
13

---

Finding Wasserstein Ball Center
A. Related Works
Wasserstein distance The Wasserstein distance, also known as the Earth Mover's distance when p = 1, quantifies the disparity between two probability distributions, particularly when their supports are discrete sets. Computing the discrete Wasserstein metric reduces to solving a minimum-cost maximum-flow optimization problem (Ahuja al., 1994; Orlin, 1988). Consider bipartite graph with n vertices, m edges, and maximum edge weight U. Lee & Sidford (2014) designed an efficient general LP solver that can solve the minimum cost flow problem in On. polyy(og U)). Lahn et al. (2023) designed combinatorial algorithm that computes an c-approximate transport plan in O ((2) sequential time and O llog
parallel time. Recently, Fox& & Lu (2023) presented the first deterministic algorithm that computes optimal transportation in Euclidean space with -eaarl time.
Wasserstein distance is also a classic topic in machine learning (Rüschendorf, 1985; Pele & Werman, 2009; Dognin et al., 2019). By using matrix scaling technique, Cuturi (2013) introduced the "Sinkhorn Distance", which incorporates an entropic regularization term to smooth the transportation problem, offering paralleled and significantly faster solutions than exact computation of the discrete Wasserstein distance. Following Cuturi's work, recent years have seen the development of several improved Sinkhorn algorithms (Lin al., 2019; Altschuler et al., 2019; Benamou et al.. 2015; Altschuler et al., 2017). Another practical variant is "sliced Wasserstein distance", which saves computation time by taking expectation of Wasserstein distance of the distributions projected from high dimension to one-dimensional space (Kolouri et al., 2019; Nietert et al., 2022).
Wassertein barycenter. Different from Wasserstein distance, the computation for Wasserstein barycenter is an NP-hard problem (Altschuler & Boix-Adsera, 2022). Cuturi & Doucet (2014) showed that the computation for WB can be improved by using an entropic regularization, leading to a simple gradient-descent scheme that was later improved and generalized under the iterative Bregman projection (IBP) algorithm (Benamou et al., 2015). Further progress includes the semi-dual gradient descent (Cuturi& & Peyré, 2018), accelerated primal-dual gradient descent (APDAGD) (Kim et al., 2025; Kroshnin et al., 2019), alternating direction method of multipliers (ADMM) (Ye et al., 2017), deterministic IBP (Lin al., 2020), and the IPM algorithm MAAIPM (Ge et al., 2019). In Euclidean space, WB can be solved in polynomial time with fixed dimension (Altschuler& & Boix-Adsera, 2021). Vaskevicius & Chizat (2024) offers the first non-asymptotic convergence guarantees for approximating Wasserstein barycenter between discrete point clouds in the free-support/grid-fre free setting. Agarwal et al. (2025) gives the first polynomial algorithm for the case where the cost is I,-nomm. In applications, we sometimes need the support of WB to be sparse (Ye et al., 2017). Borgwardt (2022) proposed a geronggy-polynoone 2-approximation algorithm for sparse Wasserstein barycenters. Yang & Ding (2024) investigated the case with outliers, connecting this problem with k-means clustering.
Interior Point Method. Dikin (1967) proposed the first interior point method, which was reinvented in Karmarkar (1984) which also proposed the first poly-time method for linear programming called "Karmarkar's algorithm"". Since then IPM has attracted a great amount of attention, where one of the most successful IPM methods is the class of primal-dual approaches, which iteratively optimize the primal and the dual variable, until the dual gap converges to zero. In particular, Mehrotra's cedlictor-correttor algorithm (Mehrotra, 1992) provides the basis for most implementations of this class of methods, which is also the type of IPM applied in this paper (the details of predictor-correcta IPM are presented in Section 3.2). The seminal Mizuno- Toooo-Tood------ Ye predictor-correctt method has quadratic convergence (Mizuno et al., 1993; Ye et al., 1993). For more details on IPM, we refer the reader to the survey paper (Gondzio, 2012). Recently, there are also some new studies on reducing the exponent of IPM complexity in theoretical computer science (Jiang et al., 2021; Cohen et al., 2021), which relies on a technique called "matrix maintenance" to reduce the update time of each iteration.
Fairness. The fairness issue has attracted a great amount of attention in computer science (Joseph et al., 2016; Mehrabi et al., 2021: Caton & Haas, 2024). The proposed solutions include adjusting labels from sensitive groups to reconstruct unbiased mapping (Dwork et al., 2012; Jiang & Nachum, 2020), and removing sensitive attributions (Krasanakis al., 2018). As is mentioned before, the WBC problem is also related with socially fair clustering (Abbasi et al., 2021; Ghadiri et al., 2021). An (p)()) log log log ) >approximattion algorithm for socially fair clustering with the fp-objective is given by Makarychev & Vakilian (2021). Unfairness can result from the issue of representation bias, which arises due to insufficient amount of data in certain groups or subgroups (Lohaus et al., 2020; Chai & Wang, 2022). Existing methods include fair data generation (Jang et al., 2021), multi-objective optimization (Martinez et al., 2020) and boosting (Gong & Kim, 2017). There are other concepts on fairness studied in computer science, for example the problems of fair k-means (Chierichetti et al., 2017; Chen et al., 2019; Backurs et al., 2019; Song et al., 2025) and fair regression (Agarwal et al., 2019; Wang al., 2024).
14

---

Finding Wasserstein Ball Center
B. Proof of Theorem 3.1
reformulate as the following.
Fi
F2
FN
G
I
G
m
A
(12)
0)
I
1
D1
D2
62
DN
ex
where F; m, & I G - & m, D; vec(D(()), 64 R is the indicator vector of the i-th coordinate, for i N. N.
Noticing that the last N rows are linearly independent, and are linearly independent to the linear space spanned by other
rows due to the last N 1 rows, it suffices to prove that Ag == A(1 nh N, n, N - 1). the matrix consisting of non-zero columns first ni - rows in A, has full row rank. Now we only need to verify the following two fact: a) After series of row transformations, we can transform matrix Ag into a matrix Åg whose elements (M 1)-th, (M + m 1)-th.- (M + (N 1)m + 1)-th rows are zeros, and elements in other positions are the same as Ag. b) The matrix A has full row rank. The proof follows Lemma 3.1 of Ge al. (2019), which we include here for completeness.
a). Define
0
1
e1
T
S.
N
m.X.m
and
7m
T
ImN
Imx
In
-Ti
I
Ln
S1
C1
- Tx
I m
SN
15

---

Finding Wasserstein Ball Center
Then
Fi
F2
FN
12
G(
HOI
G2
1
H1)
G
H
I
m
where G. = - S,T.F, H = e11m m S, Note that elements in the first rows H()) and G( l,... N are zeros. We have proved the claims a).
b). Ag obtained by removing the (M 1)tth, (M m 1)-th ... .M + (N
1)m + 1)-th rows of Ao. That is,
Fi
P2
FN
Ao
G(2)
(2)
G2
(2)
(
2
H
1
where ((2)) ( m,
1T @ 10m-1. 1). H(2)
(2) m,:) [0,.
4- and F; Im @ 1 Let
m,
Nhrow M N (m- 1)
For = 1, N, let
U1
1
1
U2
U, = 1m, 8
UN+1=
R₁
UN+1
then
(F(3)
(33
F2
(3)
Ao
(3)
G
(3)
:
(3)
(3)
N
a
where p(O) F,U, Im. & |1,0T m 1. (3)
G U,
G ()
1 T
00-1,1m-11, i 1,.., N. H
H2)) = H2)) [0m-11 1, - - and
T U
N [1, m
1
16

---

Finding Wasserstein Ball Center
Let
mm
K
K
K1
m
K
K
R2
N
Im__
I I
mi mLXm m
mm, X.mm,
then
(4)
N
AoR R2
(
H(3)
G2
(4)
11)
G4
H
T
where F(1) = K1 F(2) = Im, @ [1.0 .1. G = G. K, = 10 Am-1- (m 1)x(mm;-m)l i 1,.., N,
Let Ao be the matrix comsisting of the first (mM + 1) columns of An R2. That is,
F(4t)
F(
Åo =
G (4)
(1)
G2
(4)
The matrix ÀD possesses the following two properties:
1. Every row of Åg contains exactly one nonzero entry equal to 1, while all remaining entries are zero;
2. Every column of Âg contains no more than one nonzero entry.
Given these properties, we can establish the existence of permutation matrices P1 E M+M and QI E RmM++ +1 such that the product P1A0Q: yields the form [((++M)) O,..., 0]. Consequently, we have rank(Ag) = rank(P; = - + +M, which implies rank(Ag) M + M¹. This demonstrates that Ag achieves full row-rank.
17

---

Finding Wasserstein Ball Center
C. Algorithm for (WBC-LP): Predictor-Corrector Inner Point Method
Algorithm 2 Predictor-Correctto Inner Point Method for Linear Programming (WBC-LP)
1: Input: The standard LP form of Problem (WBC-LP)
min cta x
s.t. Ar=1d b, 2 0
2: Initialization: Set initial feasible point (zo, No. so), where 20 > 0, 80 (dual variables). Choose tolerance 0 and
set iteration counter k = O.
3: while ||ra|| > and ||ra|l > do
4: Compute residuals:
Tb Ar - (primal residual)
A y - (dual residual)
Xs (complementarity residual)
where X = diag(1),, S diag(s), and =² is the duality measure.
5: Predictor Step: Solve the linear system for aftine scaling direction (Aaall Ay As
0
A
1
A
Tc
0
0
Ayft
Tb
0
X
A
Ta
We have
Ayfff (ARA ) (re - AR(s - rb)), solved by calling Algorithm
(13)
A = - A Aya Ar = R(s As).
6: Compute the step size aat by finding the maximum step length that maintains x x+ aat > and s s aattAs 0. Corrector Step: Compute the corrector directions using central path perturbation with updated r.:
7:
V. rs + AX As
pl
and solve the system again get oA AY As com)
Ay corr - aft - solved by calling Algorithm 1.
(14)
As cuIr - Ay
Cur corr 5 Ax A A -1s-11
8: Compute the total search direction:
Ar
Azall eorr
Ay= = oy As = Ag + cor
9:
Compute the step size by updating with both predictor and corrector directions. Update variables:
10:
#a+1 = 2x + AAr, yk+1 = yk + ady, Sk+1 = St + aAs.
Update the duality measure y and increment the iteration counter h: k end while
Output: Optimal solution (x*, ,*))))).
18

---

Finding Wasserstein Ball Center
D. Proof of Proposition 3.4 Proof. Let M == E my and M' = Nmm 1). Using the same notation as the preceding section,
F2
FN
G
(
(2)
H
A
G
(
1
D1
-1
D2
62
DN
e
where G((0 (2(2((2 m. 1 0m 4m H2 (1) m,;) 10m-1: -m-11 and F; Lm, 1m
Let
F1
F2
: A(1:
FN
G
(2)
G2
(2)
(
G2 A(M + M M M
G 2
(2)
Ag = A(M + M +1:)
[0, 0,
0,
D1
D2
A. A(M M + +2:n. N,;)=
DN
eN
Then
G
G.
G,RGT
G,
G,RA
G
G2R
G2RG,
G2RA
G,RA
Ā
and ARAT
A3
AgRG
AaRG.
AgRA
AgRA
A
RGR
A,RG,
A,RA
A,RA
Now we analyze the structure of each sub-matrix A,RAJ and rename them for conciseness. Let Rt diag(R), RN4a))
where R, E 1 mmm, Xmm,, = I,... N, RN+1 E _ RN42 E NxN RN43 = Tn.
Then
F.R.FT
G RG
B
Fy
19

---

Finding Wasserstein Ball Center
Each F. is diagonal matrix with positive diagonal entries.
G P.
G2
= B2
G R FN]
G R 11 T
(2)
H RN+1H (2)
H RN+1
G2 RG
(15)
G R 1T)
H RN+1
where 22 Ry- H()) and each G R.G. (2) a diagonal matrix with positive diagonal entries. We use B3 to denote the first matrix in the right hand side of (15) and B, to denote the second. addition, other blocks of ARA are
AsROT O,
AgRG.
RN41 (2)T
1 H = aT,
AgRAj - R =1
F1
AGROT
== K
D RN
( R
AGROT
KT
D RN FN
R. diag diag(D,R DI,.D , Ry DJ) RN42 Tn. 1y W,
here have W, diag Di Ri DJ D2RaD DNRNND + RN42. With all notations above, we have
B
B2
K1
ARÄT
B2
B3
B4
a
K₂
O
0
KT
KTT
0
W
E. Derivation of Eq. (8)
IM
T
B
K1
B.
B
ARA
B
Iast
G B4 a K2
IN
Ix
KT
K
IM
B
0 Ki
Ixt
B
Ki
Ixs
-afl
G Gi+BA
B4
a K₂
Ixe
ax/l
G G1+G2
Ka
0
IN
KT
KT
Is
KT
KT
1ss
Bi
B
Gr+G
K3
T
G1+Gg
K₂
RT
B
Ix
KT
K
K.
B
IN
KJ
W
20

---

Finding Wasserstein Ball Center
F. Proof for Theorem 3.5
The proof follows the same approach of Lemma 4.2 Ge et al. (2019), which to prove positive definite property by confirming strict diagonal dominance. This approach requires several foundational results about the underlying matrix structure.
Lemma F.1. The block matrices B1, B2, B3, B4 possess the fundammnnt properties:
(i) All entries are non-negative.
(ii The identity M,lM, B3 1M holds.
(iii) The strict inequality a B2lM, > 0 is satisfied.
Proof. Property (i) follows directly from the construction using non-negative matrices A and R.
For property (ii), we compute:
G R 1
2)
BylM,
2 RNG (2)
and
( GL R m1
B2 1M
G RNFN Imn
Using the definitions ((22) 1T, @ 0m 1 and F Im, X I along with the diagonal structure of R, we establish
that G.²) F.F. 1m, G.
(,ch)R,,,R R,G 1m-11 1. yielding the desired identity.
For property (iii), observe that:
Fi FiR FT Ims
G (2)T
1m-1
By 1M
B₂ B21 M2
RN L ma
FNRNG
1m-1
Since F.R.F.T 1m, F.R.G. 1m-1 for each block the component-wise strict inequality follows.
Now proceed establish Theorem 3.5.
Proof of Theorem 3.5. The biock-diagonal structure of G1 is immediate from its construction. We establish positive definiteness via contradiction.
Assume G₁ is not positive definite, Then there exists an eigenvalue -X < 0, making the matrix Ad + G1 singular. Consider the vector equation:
(A Gi) A1 + B3l (B) B--BBB B2) 1 Ma
(16)
Applying property (ii) from Lemma F.1:
A1Ms + B.J 1M - (B.J B, B2)
AlMs + B,J B B (B) B B2 1
A1M, + B3 (B B2l
(17)
21

---

Finding Wasserstein Ball Center
Since By BT' 2 0 with each row containing at least one strietly positive element, and property (iii) ensures B
M2 0, conclude:
(AlM; + G > OM2
Next, we analyze the structure of M/, + Gi = A/ + B3 B2 B 1 B2. Since B1, B2, B3 2 0 and B3 is diagonal, this matrix has positive diagonal entries and non-positive off-diagonal entries.
Define the projection operators:
DM2 = TM (diagonal mask)
(18)
Oz = Ms 1 1 M₂ JM (off-diagonal mask)
(19)
Then:
DM, o MM + Gil DM2 o (AMM, + Gi)
(20)
O2, o [M + G1l = -0,, o (Mass + G1)
(21)
where o denotes the Hadamard product. This yields:
[M, o |M + +Gl -M, [MMM, + +GI]JMM = (Mas, (lb,+Gi)2 GG)MMM, > 00
This establishes strict diagonal dominance of AMa, + G1, implying non-singularity and contradicting our assumption. Therefore, G₁ positive definite.
As for G2, direct computation shows that:
G2 = (1N13) @ diag(z) 7=27
Given that / - and by construction, set Y diag(z) f2zT, T, we know that Y is positive definite by similar proof that of G1-
G. Algorithm Complexity
Proof of Theorem 3.2. We analyze the time and storage complexities Algorithm separately.
Time complexity. We list the most time consuming part of each step corresponding to Algorithm 1.
In step 1, computing of B1, By, W, takes most flops of O m²m)), achieved when computing B, and By.
In step 2, computing of (Gi -1 takes most flops of O(Nm²), due to the need to invert Gi,, i = I, N; computing of W takes O(N2m²) flops, which may exceed other procedures on time complexity when N is sufficiently large.
In step 3, computing W- takes ((N3) flops, while other part takes most O(Nm²).
To sum up, the total time complexity is O(m² E m, + Nm + 2 + N
Space complexity, For the implementation of the whole interior-point methods, the major data that should be kept in the memory include the following four parts:
(a) Compute primal variables or dual variables. Note that the size of primal variable is O(m(E N m))). and a dual variable is of size O(M + M').
(b) Matrix Å. According to the structure of À shown in appendix D, the matrix contains several sparse blocks, and we only need to store the non-zero entries of each block. Overall, the space for storing A is bounded by 2m Σ m, m4++20 +N (m-1)-m. Since each column of E, and Eg has at most one non-zero term, the total number of non-zero entries in E, and Eg is bounded by 2m ΣM My. addition. Eg has m - 1 non-zero zero terms, D has m E m, non-zero terms as the i-th row of contains mmn; non-zero terms. So the total number of non-zero terms in Á is bounded by O(M).
N'
22

---

Finding Wasserstein Ball Center
(c) Diagonals of matrices B₃ and B3, and diagonal blocks of matrices K1, K2 Bg and G1. B1 and B3 takes ne N O(M M") doubles. The diagonal blocks of matrices B2 and G1 have Ei= m, elements and N(m 1)2 elements, respectively. Noticing that every row of D' only has one non-zero term, we have Ki of size smaller than nnn(Eii O(M), the number of nonzero terms of E and K₂ smaller than nna (E2) O(M).
(d) The matrix W, W and W. Takes no more than 2 to store.
(c) Other intermediate vectors or matrices, whose data scale is bounded by a constant time of the data scale (a), (b), (c) and (d).
With the analysis in (a) to (d), we know the memory usage of algorithm 1 is bounded by O(m E m, + Nm + N2).
H. Woodbury identity
Given the matrices A, T. C, and V, the standard formula for Woodbury identity (Hager, 1989) is
(A (A+T = -A- ---'-(--'-- + +VA-'T)-
In the proof of Lemma 3.6, the corresponding A, T, C, and V are G1, 1y & UT, ,ar, UT,,,, 13 & U, respectively.
I. The free support WBC
For the free support WB, many previous researches (Cuturi & Doucet, 2014; Ge et al., 2019; Huang et al., 2021) apply block coordinate descent, which requires minimizing the
min w.X. (n.))
Σ (D(X, Q19),
(22)
[T 1 m. w. (II)) 1m aa(1), (I() 2 O. vt 1,..,, N
s.t.
1 w w
where (Wi,'". Wm)) T RT. X = " Fm Rdxm. (T( R and D(X, Q(1)
[| q 10
mmx, for 1,... N. Problem (22) is nonconvex problem, where one needs to find the optimal support points X and the optimal weight vector of a barycenter simultaneously, But specified from the support points of (P(*) ).V
For the free support Wasserstein barycenter, many previous researches (Cuturi & Doucet, 2014; Ge et al., 2019; Huang et al., 2021) apply block coordinate descent. For our Wasserstein ball center, the objective becomes:
min max (D(X, (())) (IE
s.t. 11 m. w. (II0) T 1m = a0), 11 2 0, = LL...
(23)
m w = 1 20
where w = (W1, mm E RT, X = [2,, 2m E R.xms (1 E R and D(X, (()) = [|| 0 (2) |||
e mm for t = i.... N. Problem (22) is nonconvex problem. By block coordinate descent, optimise the support set X, and then the fixed support WBC to obtain the weight w of WBC and coupling matrices (I) alternately. The algorithm
will converges a local minima. For instance, l2 space, the minimization of min max Σm Σ || q, 1 m
tE[N]
This is quadratically constrained quadratic program (QCQP) in that it can be reformulated as the following problem:
min C X.5
m
m,
(24)
s.t.
a q
11 SC
(1
Since the quadratic forms of X in the constraints is positive semidefinite, the problem is convex [1], thus can be efficiently solved with convex programming, such as interior point method. Note that the size of variable X, Ç is m + 1, the number of constraints is N, the scale of solving X is much smaller than solving the fixed support WBC.
23

---

Finding Wasserstein Ball Center
J. Supplementary Experiment
Log(value)
12
log(value)
128
100
Iteration
Iteration
Figure 5. N=10, m=6000
Figure 6. N=2000, m=300
68-
16-
20-
log(value)
14-
log(value)
6-
12-
t0.
Iteration
Iteration
Figure 7. N=3, m=30000
Figure N=500, m=500
J.1. Super-linearity of our algorithm
A sequence { is said to converge super-linearly limit x² if:
sn++1 2:
lim 2-000
[nn - 1=1
This means that the error ja,, - 2=| decreases at a rate faster than linear convergence. With our logarithmic coordinate of objective value, we have shown that our algorithm converges in super-linearly speed,
In both case when m > N and m < N, the super-linear convergence speed maintains.
J.2. Fairness of WBC
To compare the imf, we set all distributions be on support [m]] C R.. For distributions in family we assign the points in 21, 10) with probability proportional independent samples from ( 4)), the uniform distribution on [3, 4). and other points with probability proportional to samples from U(0,, 1)). For distribution in family 2, assign the points in {in ........ m with probability proportional to independent samples U U[33, 4J). and other points with probability proportional to samples from U[0,, I]). See Fig. 9 for the instance tested in Fig. 3(b): We also conduct several other experiments with different settings, see Figs. 10 and 11.
To evaluate and compare the (aD( we consider a comprehensive set of experimental configurations. The experiments are conducted across the following settings:
24

---

Finding Wasserstein Ball Center
famiyl's probebiitty on 2011
Femily2 probability on (121:1500)
6.
Probabilityy
Probability
0.24
6662
080
600
010
D20
1025
2544
203
JAb8
.033
283
1055
2508
Support
Support
Figure 9. The left histogram shows prohability in family on [20], whose density is concentrated in [10]. The right is probability in family 2, concentrated in (m- m].
Ground space dimensions: 3, 10, and 100.
Ground metrics: (1, f2, and to..
Support selection regimes:
Uniform sampling within cube;
Gaussian sampling;
Each sampling regime is further augmented with cluster.diversiffinng transformations, which divide the support points
into 1, 2, clusters.
For every combination of the above three parameter categories, generate 50 independent instances, resulting in total of 2700 experimental runs. Every instance consists of 30 probability distributions, each of which has weights sampled uniformly and support size of 200.
WBC
WBC
1500
WB
1500
WB
Max WD
Max WD
1000
1000
500
500
0.2
0.4
0.6
Imbalance factor
Imbalance factor
Figure 10. N=40, m=150
Figure 11. N=15, m=600
J.3. The Point cloud data
This is free-support problem. We select the support with two-step strategy: firstly, choose the support from combination of dense grid and samples from normal distribution with large variance, computing the weight of WBC; then update the support of WBC by adding the support in last step with samples near the points who have particular large weights, computing the final WBC. Both WB and WBC are given grey scales proportional to the weight of the point, so more weight means paler color.
25

---

Finding Wasserstein Ball Center
Figure 12, WBC (left) and WB (middle) of 3 airplane, one which is a biplane (with a double-decker wing structure), while the other are both monoplanes. Compared to the WB, the WBC retains characteristies of the biplane.
26


Reviews:
[{"review_id": "JVm0DWwQ1y", "reviewer": "Reviewer_Xncv", "summary": "This paper considers the following problem: Given a set of N probability measures, find a probability distribution that, minimizes the maximum distance to any input distribution. Intuitively, we can think of this as the problem of finding the center and radius of the “smallest Wasserstein ball” that encloses all the input distributions. A closely related problem is the Wasserstein barycenter problem where the objective is to minimize the sum total of the Wasserstein distances from the center to the $N$ distributions. While Wasserstein barycenter has been extensively studied, there is (to my knowledge) nothing known about solving the smallest Wasserstein ball problem. \n\nSimilar to the Wasserstein barycenter, if we assume that the support for the center is fixed, then one can formulate the exact problem as an LP and use standard methods to solve the problem. The authors make some interesting observations that gives this LP a simpler structure resulting in a faster exact algorithm using IPM for the problem. The overall execution time of their algorithm is $O(Nm^3+ N^2m^2+N^3)$ as opposed to $O(N^3m^4)$ using vanilla IPM.   \n\nThe authors apply this novel optimization problem as a means to achieve fair centers to represent the distributions including those that may be underrepresented. The authors implement their solutions and show that they outperform state-of-the-art LP solvers (such as gurobi). They also provide some evidence on the impact of their approach on the problem of fairness.\n\nRebuttal Update: Thank you for your response. I will keep my positive score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1:** Comparison between our work and 1-center problem for arbitrary metric space.\n\nThanks for the question regarding the connection between our WBC problem and the 1-center problem in arbitrary metric spaces. We will add more references and explanations in our paper. In Euclidean space, 1-center problem is often called “minimum enclosing ball problem'' or ”minmax location problem''. Their numerical or approximate algorithms has been extensively studied, see [1-5].  However, the 1-center problem under Wasserstein metric remains unexplored—precisely the gap our work aims to address. Due to the inherent nature and complexity of  Wasserstein metric, those previous 1-center algorithms cannot be directly applied to handle our proposed WBC problem (to our best knowledge).\n\n\n**Q2:** How to deal with outlier distributions? \n\n Thank you for raising this interesting question. Actually, we think it is hard to precisely distinguish between \"outlier'' and \"minority'', if there is no any prior knowledge/assumption. For example,  a minority distribution could be far away from other distributions, and thus performs quite similarly as an outlier. So we believe an efficient way should take some prior knowledge on inliers and outliers into account. \n\nWhen no prior information is available, we can also utilize the following  iterative strategy to extract \"outlier'' distributions and construct the WBC for the remaining distributions (which is similar with the \"trim'' idea in statistics): initially, we compute the WBC of the input distributions; perform the following two steps in each iteration, until stable:\n\n1. Recognize the set of distributions farthest from the current WBC as outliers;\n2. Update the WBC as the WBC of the remaining distributions.\n\nNevertheless, since we do not have any prior knowledge, the above strategy could bring some error (e.g., mistaken a minority distribution as an outlier). So we believe it deserves more further study on how to distinguish between \"outlier'' and \"minority'', and which kinds of prior knowledge can be utilized. \n\n[1] Bâdoiu, Mihai, and Kenneth L. Clarkson. \"Smaller core-sets for balls.\" Proceedings of the fourteenth annual ACM-SIAM symposium on Discrete algorithms. 2003.\n\n[2] Tansel, B.Ç. (2011). Discrete Center Problems. In: Eiselt, H., Marianov, V. (eds) Foundations of Location Analysis. International Series in Operations Research & Management Science, vol 155. Springer, New York.\n\n[3] Abboud, Amir, et al. \"On Complexity of 1-Center in Various Metrics.\" *Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques* (2023).\n\n[4] Yildirim, E. Alper. \"Two algorithms for the minimum enclosing ball problem.\" *SIAM Journal on Optimization* 19.3 (2008): 1368-1391.\n\n[5] Kumar, Piyush, Joseph SB Mitchell, and E. Alper Yildirim. \"Approximate minimum enclosing balls in high dimensions using core-sets.\" *Journal of Experimental Algorithmics (JEA)* 8 (2003): 1-1.", "claims_and_evidence": "All claims are supported by proofs.", "methods_and_evaluation_criteria": "The authors show the benefits of proposed algorithm (in terms of time and memory usage) with a gurobi implementation. This is a solid comparison. The experiment to highlight fairness seems like a reasonable one.", "theoretical_claims": "I have skimmed through the proofs presented in the main text and I did not find any issues with them.", "experimental_designs_or_analyses": "Yes. The experimental set up seems sound to me.", "relation_to_broader_scientific_literature": "The problem studied is novel and I'm not aware of any paper that directly addresses the algorithmic question of computing the smallest Wasserstein ball for fixed support. However, the problem is similar to the classical 1-center problem (perhaps the discrete 1-center problem because of the fixed support requirement) in Wasserstein metric. Since 1-center for arbitrary metric spaces is extensively studied, perhaps adding a comparison of your work with that would be good.", "essential_references_not_discussed": "I did not find any missing references.", "other_strengths_and_weaknesses": "Overall, the optimization problem considered here is novel and a very natural question to study. The result provided is interesting both from a theoretical standpoint as well as from an empirical standpoint.\n\nMy main concern that the paper does not address is the sensitivity of the objective function to the presence of an outlier distribution. A single outlier distribution can disproportionately shift the center towards the outlier. Due to this, I worry that the actual practical value of the optimization problem in real-world may be limited. Could the authors comment on this?", "other_comments_or_suggestions": "NA", "questions_for_authors": "Please address the concern raised in strength/weakness section.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VeXqJVKxxE", "reviewer": "Reviewer_f8bo", "summary": "This paper introduces the concept of \"Wasserstein Ball Center\" (WBC) as an alternative to the traditional Wasserstein Barycenter (WB) for finding a representative probability distribution from multiple input distributions. While WB minimizes the sum of Wasserstein distances from the barycenter to all input distributions, WBC minimizes the maximum Wasserstein distance to any input distribution. This \"minmax\" approach makes WBC more fair to \"minority\" distributions that differ significantly from the majority.\nThe main contributions are:\n\nFormulating the WBC problem as finding the center of the minimum-radius ball that covers all input distributions in Wasserstein space\n\nShowing that the fixed-support WBC problem can be formulated as a large-scale linear programming (LP) instance\n\nDeveloping an efficient interior point method (IPM) to solve this LP problem by exploiting structure in the constraint matrix to achieve significant speed improvements over standard approaches\n\nDemonstrating through experiments that WBC provides more equitable treatment of outlier distributions compared to WB\n\nThe authors' algorithm accelerates the IPM by O(min{N²m, Nm², m⁴}) times compared to vanilla IPM, where N is the number of distributions and m is the support size. Experiments on both synthetic and real datasets show that their algorithm is more computationally efficient than commercial solvers and that WBC indeed provides more fair representation for minority distributions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1:** Missing references\n\n Thank you for your valuable suggestions on the references, and we will add them to the revised version. \n\n\n**Q2:** Algorithm for free-support Wasserstein ball center (WBC)\n\n For the free support Wasserstein barycenter, many previous researches [1-3] apply block coordinate descent. For our WBC, we can also utilize this approach, where the objective becomes:\n\\\\[\n\\begin{array}{cl}\n& \\min\\limits_{\\boldsymbol{w}, X, \\left\\\\{ \\Pi^{(t)} \\right\\\\}} \\max\\limits_{t\\in [N]} \\left\\langle D(X,Q^{(t)}), \\Pi^{(t)} \\right\\rangle \\\\\\\\\n\\text{s.t.} & \\Pi^{(t)} \\boldsymbol{1}_{m_t} = \\boldsymbol{w}, \\left( \\Pi^{(t)} \\right)^{\\top} \\boldsymbol{1}_m = \\boldsymbol{a}^{(t)}, \\Pi^{(t)} \\geq 0, \\forall t = 1, \\cdots, N \\\\\\\\\n& \\boldsymbol{1}_m^{\\top} \\boldsymbol{w} = 1, \\boldsymbol{w} \\geq 0\n\\end{array}\n\\\\]\n\nwhere $\\boldsymbol{w} := (w_1, \\cdots, w_m)^{\\top} \\in \\mathbb{R}_+^m$, $X:= [\\boldsymbol{x}\\_1, \\cdots, \\boldsymbol{x}\\_m] \\in \\mathbb{R}^{d \\times m\\_t}$, $\\Pi^{(t)} \\in \\mathbb{R}\\_+^{m\\times m\\_t}$ and $D(X, Q^{(t)}):= [ \\\\| \\boldsymbol{x}_i - \\boldsymbol{q}^{(t)}_j \\\\|^p ]\\in \\mathbb{R}^{m \\times m_t}$ for $t=1,\\cdots, N$.\n\nFree support WBC is nonconvex. By block coordinate descent, one optimizes the support set $X$, and then the fixed support WBC to obtain the weight $\\boldsymbol{w}$ of WBC and coupling matrices $\\Pi^{(t)}$ alternately. The algorithm will converge to a local minima. For instance, in $l_2$ space, the minimization of $X$ is $\\min\\limits\\_{X}\\max\\limits\\_{t\\in [N]}\\sum\\_{k=1}^m\\sum\\_{j=1}^{m\\_t}||\\boldsymbol{x}\\_k-\\boldsymbol{q}\\_j^{(t)}||^2\\pi\\_{kj}^{(t)}$. This is a quadratically constrained quadratic program (QCQP) in that it can be reformulated as the following problem:\n\n\\\\[\n\\begin{aligned}\n&\\min\\limits\\_{X,\\zeta} \\zeta \\\\\\\\\n\\text{s.t.} \\ \\ &\\sum_{k=1}^m\\sum\\_{j=1}^{m_t}||\\boldsymbol{x}\\_k-\\boldsymbol{q}\\_j^{(t)}||^2\\pi\\_{kj}^{(t)}\\leq \\zeta\n\\end{aligned}\n\\\\]\n\nSince the quadratic forms of $X$ in the constraints is positive semidefinite, the problem is convex [4], thus can be efficiently solved with convex programming, such as interior point method. Note that the size of variable $X, \\zeta$ is $m+1$, the number of constraints is $N$, the scale of solving $X$ is much smaller than solving the fixed support WBC.\n\nDue to space limitations, please refer to our response to **Q4** of Reviewer YCD8 for our approach in the point-cloud experiment (see also Appendix I.3).\n\nFrom the theoretical perspective, \"free-support WBC''  is indeed a much more challenging problem than the fixed-support case, especially when the dimensionality is high. We further note that even the ``fixed-support WB'' problem is also very hard and the polynomial-time algorithm with high precision is only for fixed-dimension [5-6]. We appreciate the reviewer for raising this question, and will spend more effort to investigate this problem as the future work.\n\n\n**Q3:** Provide more guidance on choosing between WB and WBC for different applications - when is the fairness benefit of WBC worth potential increases in average Wasserstein distance?\n\n Thanks for this question, and we will include a discussion on it. In general, the trade-off between **fairness** (WBC) and **average performance** (WB) depends on the application:\n\n- **WBC** is preferable when minimizing *worst-case deviation* is critical (e.g., equitable resource allocation).\n- **WB** is better for applications where *average accuracy* dominates (e.g., density estimation).\n\n**Q4:** How does the proposed WBC approach handle noisy distributions or outliers that might not represent meaningful minorities but rather data corruption? \n\n Thank you for raising this interesting question. Please see the answer to **Q2** in response to Reviewer Xncv.\n\n[1] Cuturi, Marco, and Arnaud Doucet. \"Fast computation of Wasserstein barycenters.\" *International conference on machine learning*. PMLR, 2014.  \n\n[2] Ge, Dongdong, et al. \"Interior-point methods strike back: Solving the wasserstein barycenter problem.\" *Advances in neural information processing systems* 32 (2019).\n\n[3] Huang, Minhui, Shiqian Ma, and Lifeng Lai. \"Projection robust Wasserstein barycenters.\" *International Conference on Machine Learning*. PMLR, 2021.\n\n[4] Floudas, Christodoulos A., and Viswanathan Visweswaran. \"Quadratic optimization.\" *Handbook of global optimization* (1995): 217-269.\n\n[5] Altschuler, Jason M., and Enric Boix-Adsera. \"Wasserstein barycenters can be computed in polynomial time in fixed dimension.\" *Journal of Machine Learning Research* 22.44 (2021): 1-19.\n\n[6] Lin, Tianyi, et al. \"Fixed-support Wasserstein barycenters: Computational hardness and fast algorithm.\" Advances in neural information processing systems 33 (2020): 5368-5380.", "claims_and_evidence": "The paper's claims are generally well-supported by evidence:\n\n1. The theoretical claims about computational complexity improvements are supported by detailed mathematical derivations and proofs that exploit the structure of the linear programming formulation.\n\n2. The fairness advantage of WBC over WB is demonstrated through:\n\n- Visual examples showing how WBC preserves characteristics of outlier distributions (Fig. 1)\n- Quantitative results showing reduced maximum Wasserstein distance for minority distributions (Fig. 3b)\n- Case studies using the FairFace dataset showing more equitable treatment across racial groups\n- 3D point cloud visualization showing how WBC preserves features from minority shapes\n\n\nThe computational efficiency claims are validated by comparing runtime performance against the Gurobi commercial solver, with clear advantages demonstrated across various problem scales.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria are appropriate for the problem at hand:\n\n1. The paper provides a sound mathematical formulation of the WBC problem as a linear program, which is a natural approach for solving the proposed minmax optimization.\n\n2. The interior point method (IPM) with specialized matrix manipulations is an appropriate choice for this problem, as it exploits the specific structure of the constraint matrix.\n\n3. The evaluation criteria are well-chosen: Computational efficiency is measured against a strong baseline (Gurobi solver)", "theoretical_claims": "The paper contains several theoretical claims and proofs, particularly around the computational complexity of their algorithm and the structure of the constraint matrix in the proposed LP formulation.\n\nThe main theoretical result is Theorem 3.2, which states that the time complexity of each inner iteration of the IPM is O(m²∑ᵢmᵢ + Nm³ + N²m² + N³), and the memory usage is O(m∑ᵢmᵢ + Nm² + N²).\n\nThe proof of this theorem is built on several intermediate results:\n\nProposition 3.1 (showing the constraint matrix can be made full row-rank)\n\nProposition 3.4 (describing the structure of the normal equation matrix)\n\nProposition 3.5 (further decomposing the structure of specific submatrices)\n\nLemma 3.6 (providing an efficient way to compute a key matrix inverse)", "experimental_designs_or_analyses": "The experimental designs and analyses in the paper are sound:\n\nCompare the proposed algorithm against Gurobi, a state-of-the-art commercial solver\n\nEvaluate performance with varying N (number of distributions) and m (support size)\n\nMeasure both computation time and objective value\n\nDemonstrate super-linear convergence rate (Fig. 3a)", "relation_to_broader_scientific_literature": "The paper builds upon and extends several areas in the scientific literature:\n\nWasserstein Barycenter\n\nFairness in Machine Learning\n\nOptimal Transport", "essential_references_not_discussed": "The paper is generally thorough in its literature review, but a few relevant references might strengthen the connections:\n\nRecent advances in fast approximation algorithms for Wasserstein distance computation. For exmaple (https://arxiv.org/abs/2312.01432)\n\nSome recent work on fair Wasserstein barycenters (e.g., weighted schemes that give more importance to minority groups) could provide interesting contrast to the proposed minmax approach.", "other_strengths_and_weaknesses": "Strengths:\n\nThe paper addresses an important problem - the inherent bias in Wasserstein barycenters against minority distributions - with a principled mathematical approach.\n\nThe experimental results convincingly demonstrate both computational efficiency and fairness benefits across differents scenarios.\n\nThe paper connects theoretical optimal transport concepts with practical fairness implications, especially in domains like medical imaging where fairness is crucial.\n\n\nWeaknesses:\n\nThe paper focuses primarily on the discrete, fixed-support case. Some discussion on how the approach might extend to continuous distributions or free-support problems would strengthen the work.\n\nThe paper could provide more guidance on choosing between WB and WBC for different applications - when is the fairness benefit of WBC worth potential increases in average Wasserstein distance?\n\nThe computational improvements, while significant, are still limited by the inherent complexity of the LP formulation. For very large-scale problems, approximation methods might be necessary.  (https://arxiv.org/abs/2312.01432)\n\nThe experimental evaluation could benefit from more downstream tasks to show how the improved representation of minority distributions impacts practical applications.", "other_comments_or_suggestions": "NA", "questions_for_authors": "How does the proposed WBC approach handle noisy distributions or outliers that might not represent meaningful minorities but rather data corruption?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PwTIAbnfYu", "reviewer": "Reviewer_YCD8", "summary": "This paper considers fairness in representing a set of distributions and proposes to use the Wasserstein Ball Centers (WBC) as a representative of a distribution instead of the Wasserstein Barycenter (WB). Given a set of distributions $\\mu_1, \\ldots, \\mu_N$, the Wasserstein Ball Center is defined as a distribution $\\mu$ that minimizes the maximum distance to all distributions $\\mu_1, \\ldots, \\mu_N$. Although the definition of WBC seems to be new, as the authors mentioned, the concept of Wasserstein balls is not new.\n\nThe main contribution of this paper is in designing an efficient algorithm for computing the WBC of a set of discrete distributions, given a support for the WBC distribution. In particular, they present an IPM-based algorithm that computes the WBC in $O(Nm^3 + N^2m^2 + N^3) time, where $N$ is the number of distributions and $m$ is the support size of each distribution.\n\n\n## Update after rebuttal\nThe authors resolved my concerns by adding additional experimental results comparing the disparity scores of their proposed barycenter distribution and the well-known Wasserstein barycenter. I increased my score from 3 to 4.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1** On the missing part of the vector $\\boldsymbol{b}$.\n\n We sincerely appreciate the reviewer’s meticulous reading. Yes, there is a missing $\\boldsymbol{0}_N$ as the last $N$ terms of vector $\\boldsymbol{b}$. We will correct this typo in the final version. Notebly, since the design of our algorithm focuses on the constraint matrix $A$, this typo does not affect any other part of the paper.\n\n**Q2** About the experiments on the running times (experiment 1): How did you form the distributions? What are the dimensions of the distributions? Are they all supported on the same support? How did you choose the support of the WBC?\n\n  For the first three question, see Line 403-405 in paper: \"we generate random datasets in $\\mathbb{R}$, and the weights of $(q_1^{(t)}, ..., q_m^{(t)})$ in each distribution $P^{(t)}$ are generated uniformly at random.\"  \n\nThe support points of the WBC is generated in the same way, and we will clarify this in the final draft. Thank you for pointing out this ambiguity.\n\n**Q3** About the fairness metric in paper\n\n We sincerely appreciate your valuable suggestion, which highlights the need for a more comprehensive discussion of fairness in the context of WBC.\n\nOne of the inspirations of this paper is the *social fairness* concept as mentioned in the last paragraph of Introduction and Appendix A. Apart from that, *Minimax fairness* is gaining increasing attention within the machine learning community [1-3], where fairness is achieved by minimizing the maximum error across groups—as opposed to minimizing the average error (i.e., replacing $\\min\\sum$ with $\\min\\max$ in the objective). This aligns precisely with our formulation.\n\nWe acknowledge that there are other fairness notions—such as demographic parity [4], equalized odds/opportunity [5-6], and individual fairness [7]. Extending these frameworks to probability-measure spaces (e.g., via Wasserstein metrics) presents a compelling direction for future research, and we will include a brief discussion of this in the revised manuscript.\n\n**Q4** How can one derive this support set, and it is not the same as the support of the distributions?\n\nDue to space limitations, please refer to our response to **Q2** of Reviewer f8bo to view block coordinate descent method. As for our point-cloud experiment (see Appendix I.3), we proceeds as follows:\n\n1. **Initialization**: Compute the fixed-support WBC using the vertices of a dense grid as the initial support set.\n2. **Support Update**:\n   - *Augmentation*: Sample new points from Gaussian distributions centered at high-weight support points.\n   - *Pruning*: Remove support points with negligible weights (below a threshold).\n3. **Iteration**: Recompute the fixed-support WBC on the updated support and repeat until convergence.\n\nThis design is motivated by two key observations:\n\n- The WB admits a sparse support representation, with a theoretical upper bound of $N\\max\\limits_{i\\in [N]}m_i$ [8]. \n- Adaptive refinement balances computational efficiency with solution accuracy by dynamically focusing on high-density regions, which avoid too many iterations.\n\nFrom the theoretical perspective, \"free-support WBC'' is indeed a much more challenging problem than the fixed-support case, especially when the dimensionality is high. We further note that even the ``fixed-support WB'' problem is also very hard and the polynomial-time algorithm with high precision is only for fixed-dimension [8]. We appreciate the reviewer for raising this question, and will spend more effort to investigate this problem as the future work.\n\n[1] Martinez, Natalia, Martin Bertran, and Guillermo Sapiro. \"Minimax pareto fairness: A multi objective perspective.\" *International conference on machine learning*. PMLR, 2020.\n\n[2] Abernethy, Jacob D., et al. \"Active Sampling for Min-Max Fairness.\" *International Conference on Machine Learning*. PMLR, 2022.\n\n[3] Singh, Harvineet, et al. \"When do minimax-fair learning and empirical risk minimization coincide?.\" *International Conference on Machine Learning*. PMLR, 2023.\n\n[4] Zemel, Rich, et al. \"Learning fair representations.\" *International conference on machine learning*. PMLR, 2013.\n\n[5] Hardt, Moritz, Eric Price, and Nati Srebro. \"Equality of opportunity in supervised learning.\" *Advances in neural information processing systems* 29 (2016).\n\n[6]  Woodworth, B., Gunasekar, S., Ohannessian, M., and Srebro, N. (2017). Learning non-discriminatory predictors. In Proceedings of the 2017 Conference on Learning Theory, pages 1920–1953.\n\n[7] Dwork, Cynthia, et al. \"Fairness through awareness.\" *Proceedings of the 3rd innovations in theoretical computer science conference*. 2012.\n\n[8] Altschuler, Jason M., and Enric Boix-Adsera. \"Wasserstein barycenters can be computed in polynomial time in fixed dimension.\" *Journal of Machine Learning Research* 22.44 (2021): 1-19.", "claims_and_evidence": "The Theorems and Lemmas are clear.", "methods_and_evaluation_criteria": "I have concerns with the fairness evaluation criteria, as the only criterion used for measuring the fairness of the proposed method is the maximum distance of each group to the representative distribution (which is exactly the objective function that the WBC tries to minimize). It is not surprising to have lower maximum Wasserstein distance of each group to WBC compared to WB, and although one of the metrics for computing fairness is the maximum distance, many other metrics need be used to show better fairness of the proposed method over WB. \n\nI understand that the authors claim the main contribution of their paper to be the design of fast algorithms for computing WBC; however, since the authors provide experiments on comparing the fairness of their method compared to WB, they would need to provide more evidence on that front to consider it an added value to the paper.", "theoretical_claims": "I checked some of the proofs, and they seem to be correct. I have a problem with the construction of b, and I would appreciate it if the authors would let me know whether I am missing something.\n\nIn the construction of the matrix A, the first M + Nm rows check the mass preservation of the transport plans. Then there is a 1_m, which is supposed to check that the total mass of w is 1. The next N rows ensure that gamma_i's are properly defined; *here, I think there is a missing 1_N in the vector b*.", "experimental_designs_or_analyses": "The experimental results show some benefits of the proposed method over the Gurobi solver. I have some questions about the experiments on the running times (experiment 1):\nHow did you form the distributions? What are the dimensions of the distributions? Are they all supported on the same support? How did you choose the support of the WBC?", "relation_to_broader_scientific_literature": "The Wasserstein barycenter is a widely used method for representing a set of distributions and has found applications in numerous scientific fields. However, it might not be the fairest representation as it might be far from some under-represented groups. This paper suggested the use of WBC instead to resolve the fairness and provided an algorithm for computing WBC.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "The main weakness of the presented work is that the paper assumes that the support of the WBC is already given. It is not clear at all how one can derive this support set, and it is not the same as the support of the distributions. For instance, suppose we have three distributions in 2D, each one supported on a single point, forming a triangle. Then, one needs to take points inside this triangle while forming the WBC.", "other_comments_or_suggestions": "No", "questions_for_authors": "I asked my questions in previous sections.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents the Wasserstein Ball Center  as a new method for identifying a representative probability distribution from multiple inputs. The proposed approach focuses on minimizing the maximum distance promoting fairness towards minority distributions compared to 
Wasserstein barycenter which minimizes weighted average distance to distributions.
All reviewers agree that the contributions are original and solid and are happy to accept the paper.
Decision: Accept (poster)