Approximate Algorithms for k-Sparse Wasserstein Barycenter with Outliers

1 Introduction

Let P = {p1, p2, · · · , pn1 } and Q = {q1, q2, · · · , qn2 } be two sets of weighted points in the Euclidean space Rd , where we use wP (·) (resp., wQ(·)) to denote the non-negative weight function for each point of P (resp., Q); we also assume that P and Q have the same total weight, i.e., ∑n1 i=1 wP (pi) = ∑n2 j=1 wQ(qj) = n > 0. For any l ≥ 1, the seminal Wasserstein Distance [Rubner et al., 2000] is to measure their minimum transportation cost:

where || · || is the Euclidean distance, and the flow set F = {fij | 1 ≤ i ≤ n1, 1 ≤ j ≤ n2} from P to Q should satisfy: ∑n1 i=1 fij = wQ(qj) for any j, and ∑n2 j=1 fij = wP (pi)

∗Corresponding author.

for any i. We usually set l = 1 or 2 for most practical applications. Also note that "||pi - qj||l " in (1) can be replaced by other distance function if the input P and Q are in some other metric space rather than the Euclidean space. The Wasserstein distance is one of the most fundamental topics in mathematics in past decades, and recently finds a broad range of applications in machine learning, such as image retrieval [Rubner et al., 2000], generative model [Arjovsky et al., 2017], and robust optimization [Kuhn et al., 2019].

In this paper, we focus on an important optimization problem from Wasserstein distance, k-sparse Wasserstein Barycenter (WB) [Borgwardt and Patterson, 2021]. Suppose the input contains m ≥ 1 weighted point sets P1, P2, · · · , Pm, the problem of k-sparse WB for a given k ∈ Z+ is to construct a new weighted point set S with the support size |supp(S)| = k, such that

is minimized. The Wasserstein barycenter [Agueh and Carlier, 2011] is a natural representation for the average of a set of given distributions. Recently, it has been applied to a number of real-world problems, such as medical imaging [Gramfort et al., 2015], Bayesian learning [Srivastava et al., 2018], clustering [Ho et al., 2017], and natural language processing [Singh et al., 2020]. In practice, we often prefer a simple representation of the barycenter with low complexity, and thus we have the sparsity requirement "|supp(S)| = k".

Note that it is quite challenging to achieve a high-quality solution with theoretical guarantee for the objective (2); this is mainly due to two aspects, where one is from the inherent hardness for computing the Wasserstein distance, and the other reason is due to the restriction of "|supp(S)| = k" which is a troublesome combinatorial constraint for the optimization. For example, [Borgwardt and Patterson, 2021] proved that the k-sparse WB problem is NP-hard, even if m = 3 and the dimensionality d = 2; on the other hand, if the k-sparse restriction is removed, one can compute Wasserstein barycenter in polynomial time in low dimensional space [Altschuler and Boix-Adser`a, 2021].

We consider an even harder but also more practical variant of the k-sparse WB problem called "k-sparse WB with outliers". Roughly speaking, we allow a certain fraction of outliers for the mapping from each Pj to S in (2) (the formal definition is shown in Section 2). The motivation of allowing outliers is very natural in practical scenarios. Suppose we want to compute the WB for a set of images [Cuturi and Doucet, 2014]; it is common that the input images may contain noises and/or some irrelevant objects. To achieve a more robust solution, it is better to compute a barycenter that matches each input image partially (the un-matched part can be viewed as outliers).

Actually, the study on Wasserstein distance with outliers has already attracted attentions recently in machine learning [Chapel et al., 2020; Mukherjee et al., 2021; Le et al., 2021; Nietert et al., 2022]. However, the existing algorithms for computing WB with outliers are still quite limited, to the best of our knowledge. Though several robust models and algorithms for WB have been studied [Cazelles et al., 2021; Le et al., 2021], their methods do not explicitly handle outliers. In the current article, we consider developing approximate algorithms for k-sparse WB with outliers. For the sake of simplicity, we fix l = 2 for the Wasserstein distance in (1); actually our results can be easily extended for any l ≥ 1. Our main ideas rely on the vanilla clustering algorithms. Suppose we have an α-approximation algorithm A for k-means clustering (e.g., the algorithms from [Kanungo et al., 2002; Arthur and Vassilvitskii, 2007]) and a β-approximation algorithm B for k-means clustering with outliers (e.g., the algorithms from [Friggstad et al., 2019; Gupta et al., 2017]), where α, β ≥ 1. Our contributions are as follows.

• Our first contribution is to illustrate the relation between k-sparse WB with outliers and the problems of k-means clustering and k-means clustering with outliers. We are aware that the relation between the vanilla WB (without outliers) problem and k-means clustering has been discussed before [Cuturi and Doucet, 2014]. However, it is much more challenging to analyze the case with outliers; for example, even one single outlier can seriously destroy the clustering performance. In particular, we need to develop some significantly new insight to investigate the influence of outliers in the context of the complicated Wasserstein flows. We show that one can achieve an O(α)-approximate solution by utilizing the algorithm A, but the returned barycenter has the support size larger than k. If we want to keep the support size being equal to k, we can take advantage of B instead and achieve an O(β)-approximate solution.
• We further study the problem in low-dimensional space (e.g., computing the WB for a set of 2D images). Our idea follows the aforementioned clustering method, but in a more sophisticated manner. We utilize the lowdimensional coreset technique [Har-Peled and Mazumdar, 2004] to generate a set of "anchor" points, and then build a set of non-uniform grids surrounding them. For any ϵ > 0, if we relax the "k-sparse" requirement, we can compute a WB that achieves a (1 + ϵ) approximation factor based on those grids. Moreover, our result can be generalized to any metric space with constant doubling dimension (e.g., the input distributions could have small intrinsic dimension even in high-dimensional

space [Roweis and Saul, 2000]).

Remark 1. The algorithms A and B can be bi-criteria approximation algorithms, that is, they can return more than k cluster centers (if we relax the k-sparse requirement). The benefit of using bi-criteria approximation algorithms is that they often achieve lower α and β. For example, the popular k-means++ algorithm yields an O(log k) approximation factor [Arthur and Vassilvitskii, 2007]; but if we run the kmeans++ seeding procedure more than k steps (say λk with some constant integer λ > 1), the approximation factor can be reduced to be O(1) [Aggarwal et al., 2009].

1.1 RelatedWorks

Wasserstein distance. The research on computing the Wasserstein distance (1) has gained a great amount of attentions in theory and various practical applications. It is easy to see that (1) can be viewed as a min-cost flow problem that one can apply the network simplex algorithm to solve it [Ahuja et al., 1988]. In machine learning community, [Cuturi, 2013] proposed a new variant called "Sinkhorn distance" that can be computed much faster than the original Wasserstein distance. Following Cuturi's work, [Altschuler et al., 2017] proposed a nearly-linear time Wasserstein distance algorithm. [Genevay et al., 2016] proposed a stochastic algorithm for solving large-scale optimal transportation. Some recent improvements also include [Dvurechensky et al., 2018; Lin et al., 2019]. As mentioned before, several works on the Wasserstein distance with outliers problem and its applications (e.g., outlier detection and shape matching) were also proposed recently [Chapel et al., 2020; Mukherjee et al., 2021; Le et al., 2021; Nietert et al., 2022].

Wasserstein barycenter. The study of Wasserstein barycenter mainly focuses on two different types. One is called "fixed-support WB", where the barycenter S in (2) has a fixed support (the support size can be very large) and the task is to determine the weight distribution over the support such that the average Wasserstein distance to the given m input distributions is minimized. The other one is called "free-support WB" where the barycenter S can have the support that locates anywhere in the space. The former problem is relatively easier to solve, since the weight distribution can be obtained by computing a linear programming (LP) [Auricchio et al., 2019]. A number of efficient algorithms for fixed-support WB have been proposed. For example, [Claici et al., 2018] presented a stochastic algorithm for WB; [Ge et al., 2019] developed a novel interior-point method by removing redundant constraints for the LP; [Lin et al., 2020] provided a fast iterative Bregman projection algorithm.

On the other hand, the free-support WB problem is more challenging. Recently, [Altschuler and Boix-Adsera, 2022] showed that it is NP-hard to compute even a WB with ϵ additive error in Euclidean space (if d is not constant); only for low-dimensional space, one can obtain the optimal WB in polynomial time [Altschuler and Boix-Adser`a, 2021]. [Borgwardt, 2022] provided a 2-approximate WB in Euclidean space. But those algorithms cannot guarantee small support for the obtained barycenter (e.g., the support can be as large as O(∑mj=1 |supp(Pj)|)).

If we further require |supp(S)| = k, the problem can be much more challenging; as mentioned before the k-sparse WB problem is NP hard even for m = 3 and d = 2 [Borgwardt and Patterson, 2021]. Most existing algorithms for k-sparse WB rely on the idea of alternating minimization, that is, they iteratively update the location of the k sparse support of S and the weight distribution, until the solution converges to some local optimum [Cuturi and Doucet, 2014; Ye et al., 2017; Claici et al., 2018; Ge et al., 2019].

2 Preliminaries

To formally define the k-sparse WB with outliers problem, we should provide the definition for Wasserstein distance with outliers first (the vanilla Wasserstein distance (1) is the special case with zero outlier). For any weighted point set P ⊆ Rd , we use wP (p) and wP (S) to denote the non-negative weight of a point p ∈ P and the total weight of a subset S ⊂ P, respectively. We also define the relation "⪯" between two point sets P and P ′ : P ′ ⪯ P if supp(P ′ ) = supp(P) and wP ′(p) ≤ wP (p) for any p ∈ supp(P). Further, we define the set PP-z = {P ′ ⊂ Rd | P ′ ⪯ P and wP ′ (P ′ ) = wP (P )- z} for any given non-negative value z ≤ wP (P ).

Definition 1 (Wasserstein distance with z outliers). Let n > z ≥ 0. Suppose P = {p1, p2, · · · , pn1 } and Q = {q1, q2, · · · , qn2 } are two sets of weighted points in Rd ; also ∑n1 i=1 wP (pi) = n and ∑n2 j=1 wQ(qj ) = n - z. The Wasserstein distance with z outliers from P to Q is

The set P Q = argminP ′∈PP-z W(P ′ ,Q) is called "the inliers of P induced by Q".

Remark 2. The objective function (3) was similarly defined as "unbalanced optimal transport" before [Chapel et al., 2020; Pham et al., 2020]. This definition is similar with the "trim" idea widely used in robust statistics [Rousseeuw and Leroy, 1987]. Also our definition for the function W-z(P, Q) is one-side, i.e., only P contains outliers; actually, a more general definition can be two-side that both P and Q may contain outliers. We refer the reader to the recent articles [Mukherjee et al., 2021; Nietert et al., 2022] for more detailed discussion. Due to the space limit, we only present our results of one-side here, and leave the results for two-side (which can be easily extended from the one-side results) to our full paper [Yang and Ding, 2024].

Computing W-z(P,Q). Obviously, the total flow ∑n1 i=1 ∑n2 j=1 fij = n - z in Definition 1. So the missing flows with the total weight z can be viewed as the outliers. Actually the problem of Wasserstein distance with z outliers can be easily reduced to the vanilla Wasserstein distance problem (1) via a "dummy point" idea that was studied in [Chapel et al., 2020; Ding et al., 2023] before. We add a dummy point q∗ to Q with the weight equal to z; also we force the "distance" between q∗ and each pi to be 0. Note that we cannot find such a real point q∗ in the space, where in reality we just need to set all the entries corresponding to the line of q∗ to be 0 in the n1 × (n2 + 1) distance matrix. Then we can run any off-the-shelf Wasserstein distance algorithm, e.g., the simplex network algorithm [Ahuja et al., 1988] or the sinkhorn distance algorithm [Cuturi, 2013], to compute the solution. Intuitively, the dummy point q∗ absorbs the furthest z outliers from P .

Claim 1. ComputingW-z(P,Q) is equivalent to computing W(P, Q ∪ {q∗}).

Definition 2 (k-sparse WB with z outliers). Let n > z ≥ 0. Suppose the input P contains m ≥ 1 weighted point sets P1, P2, · · · , Pm where each Pj has the total weight n. Then the problem of k-sparse Wasserstein Barycenter with z outliers for a given k ∈ Z+ is to construct a new weighted point set S with total weight n - z and the size |supp(S)| = k, such that

is minimized. Throughout this paper, we always use "Sopt" to denote the optimal solution of (4).

Remark 3. (Fixed-support WBwith outliers) As mentioned in Section 1.1, if we remove the "k-sparse" requirement and let the support of S be fixed to a given set G (|G| can be larger than k), the problem of WB can be solved by a linear programming [Ge et al., 2019; Lin et al., 2020]. Through the idea of Claim 1, the fixed-support WB with outliers problem can be also solved by using LP. Namely, we add a dummy point g∗ to G and set its weight to be z; then the problem is exactly equivalent to the vanilla fixed-support WB problem on G ∪ {g∗}. For completeness, we provide the detailed formulation in our full paper [Yang and Ding, 2024].

Relation to k-means clustering with z outliers. To better illustrate our algorithms for k-sparse WB with outliers, we need to elaborate on its relation to k-means clustering with outliers first. Given a set P of weighted points in Rd with the total weight n, the goal of the vanilla k-means is to find k cluster centers C = {c1, c2, · · · , ck}, such that each point of P is assigned to its nearest cluster center and the total weighted squared distances S(P, C) = ∑p∈P wP (p) · min1≤s≤k ||p - cs||2 is minimized. If we allow to discard z outliers, the goal becomes to find not only the k cluster centers, but also a set P ′ ∈ PP-z , such that S(P ′ ,C) is minimized. We denote this optimal cost as Meank-z (P), and let Ck-z (P) denote the set of optimal cluster centers {c1, c2, · · · , ck} with the weight wCk-z (cs) = the total weight of the s-th cluster, 1 ≤ s ≤ k. If z = 0, we use Meank (P) and Ck (P) for simplicity.

1. First, we consider the basic case m = 1 for k-sparse WBwith z outliers. It is easy to see that it is equivalent to the k-means clustering with z outliers on P1.

Claim 2. Suppose m = 1. The set Ck-z (P1) forms the optimal solution for k-sparse WB with z outliers. Namely, for any |supp(Q)| = k, Meank-z (P1) = W2-z (P1, Ck-z (P1)) ≤ W2-z (P1, Q).

(2) Then we consider the general case m ≥ 2. From Claim 2, we know that the barycenter actually induces k clusters and z outliers on P1 for the case m = 1. When m ≥ 2, the k points of the barycenter also induce k clusters on ∪mj=1Pj but with additional constraint.

Claim 3. The optimal solution of k-sparse WB with z outliers is equivalent to solving the k-means clustering with mz outliers on ∪mj=1Pj with the following constraint: for each obtained cluster Us, 1 ≤ s ≤ k, wUs (Us ∩ P1) = wUs (Us ∩ P2) = · · · = wUs (Us ∩ Pm).

Actually the above constrained k-means clustering with outliers problem is a special fairness clustering with outliers problem [Bera et al., 2019]: suppose each Pj has a unique color, and we require that each color only takes 1m of the total weight in each cluster Cs. Though a number of algorithms have been proposed for fairness clustering, the study on the case with outliers is still quite limited, to the best of our knowledge.

3 Our Clustering Based LP Algorithm

In this section we propose a clustering based LP algorithm for solving the k-sparse WB with outliers problem, where our main idea is inspired by the observations of Claim 2 and Claim 3. Though the algorithm is simple, the analysis on the quality is the major challenge since the inliers and outliers are mixed without any prior knowledge.

The clustering based LP algorithm. Let P = {P1, P2, · · · , Pm} be an instance of the k-sparse WB with outliers problem as Definition 2. Assume wmin = min1≤j≤m,p∈Pj wPj (p), and denote by ˆ z = ⌈z/wmin⌉ for convenience. Our algorithm has the following three steps.

1. We first run an α-approximate (k + ˆ z)-means clustering algorithm A on each Pj and obtain the set Tj of its λ(k + ˆ z) cluster centers with some integer λ ≥ 1 (as discussed in Remark 1, we can run a bi-criteria approximation algorithm).
2. Then, for each Tj we consider the following fixedsupport WB with outliers problem (as described in Remark 3): the support of the barycenter is fixed to be the O(k + ˆ z) points of Tj , and compute the optimal weight distribution over Tj via LP.
3. Let ˜ T1, · · · , ˜ Tm be the obtained m candidate WBs from Step (2). We return the best one, say ˜ Tj0 , which has the smallest cost over the m candidates with respect to the cost (4).

The theoretical quality guarantee of ˜ Tj0 is given in Theorem 1. Since |supp( ˜ Tj0 )| = O(k + ˆ z) that violates the k-sparse requirement, we can replace the algorithm A by a β-approximate k-means clustering with z outliers algorithm B in the above method. Then each Tj should have the support size exactly equal to k for j = 1, 2, · · · ,m. We still select the best candidate WB ˜ Tj0 by the same manner, and return it as the solution for k-sparse WB with outliers. The improved result is shown in Theorem 2.

Theorem 1. Our clustering based LP Algorithm returns a solution ˜ Tj0 for k-sparse WB with outliers and achieves the following quality guarantee:

We have multiple choices for the algorithm A. For example, we can run the (9 + ϵ)-approximate local search algorithm [Kanungo et al., 2002] (but its running time is super linear since there are too many swap combinations that should be tested in the local search procedure). We also can run the k-means++ based algorithms [Aggarwal et al., 2009; Lattanzi and Sohler, 2019] to achieve an O(1)-approximation as discussed in Remark 1. So the approximation factor in Theorem 1 can be O(1) as well.

To prove Theorem 1, we need several key lemmas. Lemma 1 shows that each obtained Tj can approximately represent the corresponding Pj, even in the presence of outliers. Lemma 2 further shows that the set T Sopt j , i.e., the inliers of Tj induced by Sopt (see Definition 1), should yield an upper bound for the total cost where the bound is determined by the distance between Pj and Sopt. Also note that we cannot obtain T Sopt j in reality since the optimal solution Sopt is always unknown to us; in fact we only use T Sopt j in our analysis for bridging the gap between ˜ Tj0 and Sopt. Through Lemma 2 we can prove that the selected best candidate ˜ Tj0 yields the desired quality guarantee.

Lemma 1. For each 1 ≤ j ≤ m, suppose the obtained cluster centers from A is Tj = {t1, t2, · · · , tλ(k+ˆ z)}; also each weight wTj (ts) = the total weight of the s-th cluster. Then W-z (Tj , Sopt) ≤ (1 + √α)W-z (Pj , Sopt).

Proof. First, we consider the relationship between k-means clustering with z outliers and (k + ˆ z)-means clustering. Intuitively, we can regard the result of k-means clustering with z outliers as a special solution for the (k + ˆ z)-means clustering, where each outlier actually is a cluster of single point. We then have the following claim (due to the space limit, the proof is shown in the full paper [Yang and Ding, 2024]).

From Claim 2 we know that Meank-z (Pj) ≤ W2-z (Pj, Sopt). Also, because Tj is obtained from the α-approximate algorithm A, we have

According to Definition 1, we know that W-z(Tj, Sopt) = W(T Sopt j , Sopt), where T Sopt j is the inliers of Tj induced by Sopt. Note P Sopt j is a set with total weight = n - z, so we have W(T Sopt j , Sopt) ≤W(T P Sopt j j , Sopt). Thus,

We also have the following bound

where the first inequality follows from the triangle inequality of Wasserstein Distance, and the last inequality follows from (6). Finally, we complete the proof by combining (7) and (8).

Lemma 2. For any 1 ≤ j ≤ m, Cost-z(P, T Sopt j ) ≤ (2 + √α)Cost-z (P, Sopt) + (2 + 3√α + α)W2-z (Pj , Sopt).

Proof. For any 1 ≤ j1 ≤ m, we have

where the second inequality follows from the triangle inequality of Wasserstein distance and the third inequality follows from Lemma 1. Then we obtain the following bound by using (9):

where the second inequality follows from the fact that (a + δb)2 ≤ (1 + δ)a2 + (δ2 + δ)b2 for any numbers a, b, and δ.

Proof. [of Theorem 1] Because ˜ Tj is the optimal weight distribution over Tj, we have Cost-z(P, ˜ Tj) ≤ Cost-z (P, T Sopt j ). For the best candidate ˜ Tj0 , we have

based on Lemma 2. Also it is easy to know min1≤j≤mW2-z (Pj, Sopt) ≤ 1m ∑mj=1 W2-z (Pj, Sopt) = Cost-z(P, Sopt), so (11) implies Cost-z(P, ˜ Tj0 ) ≤ (2 + √α)2Cost-z(P, Sopt) which completes the proof.

Theorem2. If we run the β-approximate k-means clustering with z outliers algorithm B instead of A, and compute the optimal weight distribution with 2z outliers over Tj in the clustering based LP algorithm, we have

Comparing with Theorem 1, it is guaranteed that the output ˜ Tj0 in Theorem 2 is exactly k-sparse, with only a violation on the size of outliers. The total weight of discarded outliers is increased to 2z; but usually z is a value much smaller than n and thus we believe this influence is acceptable in practice. For the algorithm B, we can use the O(1)-approximate algorithms [Friggstad et al., 2019; Krishnaswamy et al., 2018] (the algorithm of [Friggstad et al., 2019] slightly violates the number of returned cluster centers to be (1 + ϵ)k with an arbitrarily small value ϵ > 0). In practice, we can also use some faster algorithms like [Chawla and Gionis, 2013; Gupta et al., 2017] (though their theoretical guarantees are weaker than [Friggstad et al., 2019; Krishnaswamy et al., 2018]).

The proof of Theorem 2 is similar with that of Theorem 1, but the only major challenge is that we have to provide a more complicated version for Lemma 1 (which is Lemma 3 below). For each 1 ≤ j ≤ m, let the obtained cluster centers from B be Tj = {t1, t2, · · · , tk}; the key difficult problem is that Tj and Sopt may induce different inliers on Pj. To resolve this issue, we should make a deep analysis on the distribution of Tj and prove the existence of a set T ′ j ⪯ Tj who can play the same role as Tj in Lemma 1. Due to the space limit, we leave the proofs of Lemma 3 and Lemma 4 to our full paper.

Lemma 3. There exists a weighted point set T ′ j satisfying W-z(Sopt, T ′ j ) ≤ (1 + √β)W-z(Pj, Sopt), where supp(T ′ j ) = supp(Tj ) and wT ′ j (T ′ j ) = n - 2z.

Lemma 4. For any 1 ≤ j ≤ m, the weighted point set T ′ j in Lemma 3 also satisfies Cost-2z(P, T ′ j ) ≤ (2 + √β)Cost-z (P, Sopt) + (2 + 3√β + β)W2-z (Pj , Sopt).

Proof. [of Theorem 2] Because of ˜ Tj is the optimal weight distribution over Tj, so we have Cost-2z(P, ˜ Tj) ≤ Cost-2z (P, T ′ j ). For the best candidate ˜ Tj0 , we have

where the second inequality follows from Lemme 4 and the last inequality follows from min1≤j≤mW2-z (Pj, Sopt) ≤ Cost-z (P, Sopt).

Analysis on running time. Let Γ1 be the time complexity of A or B on each input set Pj, and let Γ2 be the time complexity for solving the fixed-support WB for each ˜ Tj as described in Remark 3. Then the total time complexity of our algorithm is m(Γ1 + Γ2). Usually Γ1 can be linear in n if using the previous clustering algorithms (e.g., [Aggarwal et al., 2009]) ; the complexity Γ2 can be ˜O(mn7/3 ϵ-4/3 ) if using the recent algorithm [Lin et al., 2020] (ϵ > 0 is a prespecified small error).

4 Improvement in Low-Dimensional Space

We further consider a common case that the dimensionality d is small (e.g., the input P is a set of 2D images). In the previous section, we show that the quality of our solution heavily depends on the clustering performance of the algorithm A or B. So a natural question is

Can we remove the dependence on the factors α and β in our quality guarantee?

We answer this question in the affirmative. Our main idea is to generate the support through a more sophisticated approach. Our algorithm contains the following two steps.

(1) Generate the anchor points. We still run the αapproximate (k+ˆ z)-means clustering algorithm A on each Pj and obtain the set Tj = {t1, · · · , tλ(k+ˆ z)} of its cluster centers. To have a higher-quality fixed support for replacing Tj, we utilize the low-dimensional coreset 1 technique [Har-Peled and Mazumdar, 2004] to generate a set of "anchor" points ˆ Tj as follows. Denote by B(c, r) the ball centered at point c with radius r ≥ 0. We fix a j, and let ¯ rj = √W(Pj, Tj)/n and ϵ1 > 0. Then for s = 1, 2, · · · , λ(k + ˆ z), we partition the ball B(ts, n¯ rj) into ⌈logn⌉ + 1 layers: Tj,s,0 = B(ts, ¯ rj) and Tj,s,h = B(ts, ¯ rj2h )\B(ts, ¯ rj2h-1 ) for h = 1, · · · , ⌈log(n)⌉. For each h-th layer, we can build a grid with the side length ¯ rjϵ12h-1 /√αd; each point of Pj ∩Tj,s,h is assigned to its nearest grid point, and each grid point has the weight equal to the total weight of the points assigned to it. Finally, we have the set ˆ Tj which contains all the weighted grid points of the ⌈logn⌉ + 1 layers. The size | ˆ Tj| = O((k + ˆ z) log(n)αd/2 /ϵd1 ). We call the union set ∪mj=1 ˆ Tj as the "anchor points" (which is actually the coreset in [Har-Peled and Mazumdar, 2004]).

(2) Construct the support. Without loss of generality, we assume that the minimum and maximum pairwise distances of ∪mj=1Pj are 1 and ∆, respectively. For each anchor point q ∈ ∪mj=1 ˆ Tj , we draw ⌈log ∆⌉ + 1 concentric balls B(q, 2h ), h = 0, 1, · · · , ⌈log∆⌉. Let ϵ2 > 0. Inside each ball B(q, 2h ), we build a grid of side length ϵ22h-1 /√d. We denote the union set of the ⌈log∆⌉ + 1 grids as Gq, and denote ¯G = ⋃q∈⋃mj=1 ˆ Tj Gq.

Finally, we solve the fixed-support WB with outliers problem by using the LP method on ¯G instead of the Tjs. The obtained solution is denoted by ˜G.

Theorem 3. Our Algorithm returns a solution ˜G that has the following quality guarantee by setting ϵ1 = ϵ2 = ϵ/16:

1 Coreset is an algorithmic technique for representing large-scale data, which has been widely used for the optimization problems like clustering and regression [Feldman, 2020].

Before proving Theorem 3, we provide the following two key lemmas first. Lemma 5 shows that each obtained ˆ Tj can efficiently preserve the Wasserstein distance error for Pj within any arbitrarily small bound, even in the presence of outliers. Lemma 6 further shows that the instance T = { ˆ T1, ˆ T2, · · · , ˆ Tm} yields a barycenter on the fix support ¯G which can approximately represent the barycenter on the input instance P. The detailed proofs are placed to our full paper [Yang and Ding, 2024].

Lemma 5. For each 1 ≤ j ≤ m, we have

Lemma 6. Let T = { ˆ T1, ˆ T2, · · · , ˆ Tm} be a new instance of k-sparse WB with outliers, then we have

Proof. [of Theorem 3] We can combine Lemma 5 and Lemma 6 to complete the proof. Let δ > 0 be a parameter that will be determined later. First, we have

where the second inequality follows from the generalized triangle inequality for any real numbers δ > 0, a, and b: (a + b)2 ≤ (1 + δ)a2 + (1 + 1δ )b2 . Then from (17) we have Cost-z(P, ˜G) ≤ (1 + δ) 1m ∑mj=1 W2-z (Pj , Sopt) + (1 + 1δ )Cost-z(T, ˜G) ≤ ((1 + δ)1.25ϵ21 + (1 + 1δ )(1 + ϵ2)2 (1 + √1.25ϵ1)2 )Cost-z(P, Sopt), where the second inequality follows from Lemma 5 and Lemma 6. Finally we choose δ = 1/(√1.25ϵ1(1 + ϵ2)(1 + √1.25ϵ1)), then we have Cost-z(P, ˜G) ≤ (√1.25ϵ1 + (1 + ϵ2)(1 + √1.25ϵ1))2Cost-z(P, Sopt). By setting ϵ1 = ϵ2 = ϵ/16 we obtain Theorem 3.

Running time analysis. The running time is similar with the complexity of Section 3, where the only difference is adding the time complexity for building the anchor points and constructing the support in our algorithm. Note that the complexity of [Har-Peled and Mazumdar, 2004] is linear in the input size n for each Pj. The total complexity for this extra part is ˜O( log(∆)(k + ˆ z)αd/2 ϵ-2d + n).

Extension in doublingmetric. Our result can be easily extended to the more general case in doubling metric. Informally speaking, the "doubling dimension" measures the intrinsic dimension of data (Euclidean dimension is one kind of special doubling dimension) [Gupta et al., 2003]. We show the extension with details in our full paper.

5 Experiments

In this section, we illustrate the practical performance of our algorithms and study the significance of considering outliers for WB. Our experiments contain three parts. Firstly, we conduct the experiments on synthetic datasets, where the positions of the barycenter supports are predefined, allowing us to compute the exact optimal objective value for measuring the approximation ratio of our algorithm. Secondly, we compare our algorithms with several baselines on real-world datasets. Finally, we provide the visualized results on the MNIST dataset [LeCun et al., 2010]. Some omitted experimental results are placed to our full paper [Yang and Ding, 2024].

Datasets. In our synthetic datasets, we set the supports size k ∈ [10, 40] and the dimensionality d ∈ [10, 40]; each instance comprises m ∈ [2, 10] different distributions, where each distribution consists of n = 20, 000 points. The true barycenter supports are uniformly sampled within a hypercube with a side length of 10, and random weights are assigned to each center point. The points are randomly generated within Gaussian balls around the centers based on the assigned weights. We introduce outliers by uniformly sampling z points for each distribution within the cube, with z ranging from 0 to 0.15 × n.

We also select three widely-used datasets from the UCI repository [Dua and Graff, 2017]: Bank [Moro et al., 2012] (4, 521 points in R3 ) represents the individual telephone calls during a marketing campaign, which contains the information of the customers. We have m = 3 distributions categorized based on marital status. Credit card [Yeh, 2016] (30, 000 points in R14 ) includes the information about the credit card holders. We partitioned the data into m = 9 distributions based on marriage and education. Adult [Becker and Kohavi, 1996] (32, 561 points in R5 ) represents the individual information from the 1994 U.S. Census. We partitioned it into m = 10 distributions based on sex and race. Finally, 5% random noise are added to each dataset as outliers.

Baselines and our implementation. It is worth noting that there is no method that explicitly addresses k-sparse WB with outliers or fair clustering with outliers, to the best of our knowledge. We employ three baselines. First, following Remark 3, we consider the fixed-support WB with outliers algorithm, utilizing k random centers as support, and compute the optimal weight distribution via LP (denoted as "Random O"). The other two baselines include a fair clustering algorithm that does not consider outliers (denoted as "FC O") [Bera et al., 2019], and a non-fair clustering method considering outliers "k-means- - O" [Chawla and Gionis, 2013]. For the FC algorithm, we identify the farthest points in each class as outliers; for k-means- -, after obtaining the support positions, a new fair clustering solution can be obtained through LP. Additionally, we also test their three "plain" versions that do not discard outliers, aiming to study the significance of considering outliers (denoted as "Random", "FC", "k-means- -", respectively).

In our implementation, we use the k-means++[Arthur and Vassilvitskii, 2007] as Algorithm A and k-means- -[Chawla and Gionis, 2013] as Algorithm B; we also employ the LP

solver [Gurobi Optimization, LLC, 2023] as the subroutine for solving fixed-support WB with outliers. To ensure a fair comparison, although Theorem 2 suggests removing 2z outliers, we remove only z outliers in reality. Also, to keep ksparsity for the result returned by A, we only retain the top k centers with the largest cluster sizes. We use "Our A" and "Our B" to denote them.

Results on synthetic datasets. We compute the optimal cost Costopt for the WB problem by using the pre-specified barycenter support. Subsequently, we execute our algorithm under various parameters to obtain the Cost and calculate the approximation ratio, defined as Cost Costopt . Part of the results obtained by Algorithm A is presented in Table 1. We can see that our algorithm consistently achieves favorable approximation ratios across different dimensions and outlier proportions, where more than 70% of them are less than 1.5.

Results on real datasets. The results are illustrated in Figure 1. As can be seen, even with only 5% outliers, the plain versions of the three baselines take almost double costs than their counterparts who consider outliers. Moreover, our algorithms demonstrate even lower costs across all the datasets with different values of k.

Visualized results. In Figure 2 and Figure 3, we show the 40-sparse barycenters obtained by Our A and Our B for digit 0-9 in the MNIST dataset, with 2% of outliers removed from each digit. It is evident that the obtained set of 40 points effectively captures the distinctive features for each digit.

6 Conclusions

In this paper, we study the problem of k-sparse WB with outliers and present several efficient approximate algorithms with theoretical quality guarantees. Some omitted proofs are placed to our full paper [Yang and Ding, 2024]. Following this work, there are several interesting problems deserved to study in future. For example, inspired by the local search method for designing the PTAS algorithm for ordinary k-means clustering with outliers [Friggstad et al., 2019], an interesting theoretical question is that whether we can also apply it to achieve a PTAS for k-sparse WB with outliers in low-dimensional space. Our current paper focuses more on the theoretical quality for computing WB with outliers, so an ignored question is that how to improve the time complexity (e.g., using some randomization techniques). From the perspective of applications, it is deserved to consider that using WB with outlier to handle some complex data fusion problems [Cheng et al., 2021; Mao et al., 2018; Buchin et al., 2019].

Acknowledgments

The research of this work was supported in part by the National Key Research and Development Program of China 2021YFA1000900, the National Natural Science Foundation of China 62272432, and the Natural Science Foundation of Anhui Province 2208085MF163.