Coresets for Clustering with Missing Values
Year: 2021

Abstract:
We provide the first coreset for clustering points in R that have multiple missing values (coordinates). Previous coreset constructions only allow one missing coordinate. The challenge in this setting is that objective functions, like k-MEANS, are evaluated only on the set of available (non-missing) coordinates, which varies across points. Recall that an -coreset of a large dataset is a small proxy, usually a reweighted subset of points, that (1 + )-approximates the clustering objective for every possible center set. Our coresets for k-MEANS and k-MEDIAN clustering have size (jk)( −1d log n), where n is the number of data points, d is the dimension and j is the maximum number of missing coordinates for each data point. We further design an algorithm to construct these coresets in near-linear time, and consequently improve a recent quadratic-time PTAS for k-MEANS with missing values [Eiben et al., SODA 2021] to near-linear time. We validate our coreset construction, which is based on importance sampling and is easy to implement, on various real data sets. Our coreset exhibits a flexible tradeoff between coreset size and accuracy, and generally outperforms the uniformsampling baseline. Furthermore, it significantly speeds up a Lloyd’s-style heuristic for k-MEANS with missing values.

heading: 1 Introduction
text: We consider coresets and approximation algorithms for k-clustering problems, particularly k-MEANS1 and more generally (k, z)-CLUSTERING (see Definition 2.1), for points in Rd with missing values (coordinates). The presence of missing values in data sets is a common phenomenon, and dealing with it is a fundamental challenge in data science. While data imputation is a very popular method for handling missing values, it often requires prior knowledge which might not be available, or statistical assumptions on the missing values that might be difficult to verify [All01, LR19]. In contrast, our worst-case approach does not requires any prior knowledge. Specifically, in our context of clustering, the distance dist(x, c) between a clustering center point c and a data point x is evaluated only on the available (i.e., non-missing) coordinates. Similar models that aim to minimize clustering costs using only the available coordinates have been proposed in previous work [HB01, Wag04, CCB16, WLH+19], and some other relevant works were discussed in a survey [HC10].
Clustering under this distance function, which is evaluated only on the available coordinates, is a formidable computational challenge, because distances do not satisfy the triangle inequality, and
1In the usual k-MEANS problem (without missing coordinates), the input is a data set X ⊂ Rd and the goal is to find a center set C ⊂ Rd, |C| = k that minimizes the sum of squared distances from every x ∈ X to C.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
therefore many classical and effective clustering algorithms, such as k-MEANS++ [AV07], cannot be readily applied or even be defined properly. Despite the algorithmic interest in clustering with missing values, the problem is still not well understood and only a few results are known. In a pioneering work, Gao, Langberg and Schulman [GLS08] initiated the algorithmic study of the k-CENTER problem with missing values. They took a geometric perspective and interpreted the k-CENTER with missing values problem as an affine-subspace clustering problem, and followup work [GLS10, LS13] has subsequently improved and generalized their algorithm. Only very recently, approximation algorithms for objectives other than k-CENTER, particularly k-MEANS, were obtained for the limited case of at most one missing coordinate in each input point [MF19] or for constant number of missing coordinates [EFG+21].
We focus on designing coresets for clustering with missing values. Roughly speaking, an -coreset is a small proxy of the data set, such that the clustering objective is preserved within (1± ) factor for all center sets (see Definition 2.2 for formal definition). Efficient constructions of small -coresets usually lead to efficient approximations schemes, since the input size is reduced to that of the coreset, see e.g. [HJLW18, FRS19, MF19]. Moreover, apart from speeding up approximation algorithms in the classical setting (offline computation), coresets can also be applied to design streaming [HM04, FS05, BFL+17], distributed [BEL13, RPS15, BLK18], and dynamic algorithms [Cha09, HK20], which are effective methods/models for dealing with big data, and recently coresets were used even in neural networks [MOB+20].

heading: 1.1 Our Results
text: Coresets. Our main result, stated in Theorem 1.1, is a near-linear time construction of coresets for k-MEANS with missing values. Here, an -coreset for k-MEANS for a data set X in Rd with missing coordinates is a weighted subset S ⊆ X with weights w : S → R+, such that
∀C ⊂ Rd, |C| = k, ∑ x∈S w(x) · dist2(x,C) ∈ (1± ) ∑ x∈X dist2(x,C),
where dist(x, c) := √∑
i:xi not missing (xi − ci) 2, and dist(x,C) := minc∈C dist(x, c); note that the
center set C does not contain missing values. More generally, our coreset also works for (k, z)CLUSTERING, which includes k-MEDIAN (see Definition 2.1 and Definition 2.2). Throughout, we use Õ(f) to denote O(f poly log f).
Theorem 1.1 (Informal version of Theorem 3.1). There is an algorithm that, given 0 < < 1/2, integers d, j, k ≥ 1, and a set X ⊂ Rd of n points each having at most j missing values, it constructs with constant probability an -coreset for k-MEANS onX of sizem = (jk)O(min{j,k}) ·( −1d log n)2, and runs in time Õ ( (jk)O(min{j,k}) · nd+m ) .
Our coreset size is only a low-degree polynomial of d, and log n, and can thus deal with moderatelyhigh dimension or large data set. The dependence on k (number of clusters) and j (maximum number of missing values per point) is also a low-degree polynomial as long as at least one of k and j is small. Actually, we justify in Theorem 1.2 that this exponential dependence in min{j, k} cannot be further improved, as long as the coreset size is in a similar parameter regime, i.e., the coreset size is of the form f(j, k) · poly( −1d log n). We provide the proof of Theorem 1.2 in the full version. Theorem 1.2. Consider the k-MEANS with missing values problem in Rd? where each point can have at most j missing coordinates. Assume there is an algorithm that constructs an -coreset of size f(j, k) · poly( −1d log n), then f(j, k) can not be as small as 2o(min(j,k)).
Furthermore, the space complexity of our construction algorithm is near-linear, and since our coreset is clearly mergeable, it is possible to apply the merge-and-reduce method [HM04] to convert our construction into a streaming algorithm of space poly log n. Prior to our result, the only known coreset construction for clustering with missing values is for the special case j = 1 [MF19]2 and has size kO(k) · ( −2d log n). Since our coreset has size poly(k −1d log n) when j = 1, it improves the dependence on k over that of [MF19] by a factor of kO(k).
2In fact, [MF19] considers a slightly more general setting where the input are arbitrary lines that are not necessarily axis-parallel.
Near-linear time PTAS for k-MEANS with missing values. Very recently, a PTAS for kMEANS with missing values, was obtained by Eiben, Fomin, Golovach, Lochet, Panolan, and Simonov [EFG+21]. Its time bound is quadratic, namely O(2poly(jk/ ) · n2d), and since our coreset can be constructed in near-linear time, we can speedup this PTAS to near-linear time by first constructing our coreset and then running this PTAS on the coreset.
Corollary 1.3 (Near-linear time PTAS for k-MEANS with missing values). There is an algorithm that, given 0 < < 1/2, integers d, j, k ≥ 1, and a set X ⊂ Rd of n points each having at most j missing values, it finds with constant probability a (1 + )-approximation for k-MEANS on X , and runs in time Õ ( (jk)O(min{j,k}) · nd+ 2poly(jk/ ) · dO(1) ) .
Experiments. We implement our algorithm and validate its performance on various real and synthetic data sets in Section 4. Our coreset exhibits flexible tradeoffs between coreset size and accuracy, and generally outperforms a uniform-sampling baseline and a baseline that is based on imputation, in both error rate and stability, especially when the coreset size is relatively small. In particular, on each data set, a coreset of moderate size 2000 (which is 0.5%-5% of the data sets) achieves low empirical error (5%-20%). We further demonstrate an application and use our coresets to accelerate a Lloyd’s-style heuristic adapted to the missing-values setting. The experiments suggest that running the heuristic on top of our coresets gives equally good solutions (error < 1% relative to running on the original data set) but is much faster (speedup > 5x).

heading: 1.2 Technical Overview
text: Our coreset construction is based on the importance sampling framework introduced by Feldman and Langberg [FL11] and subsequently improved and generalized by [FSS20, BJKW21]. In the framework, one first computes an importance score σx for every data point x ∈ X , and then draws independent samples with probabilities proportional to these scores. When no values are missing, the importance scores can be computed easily, even for general metric spaces [VX12b, FSS20, BJKW21]. However, a significant challenge with missing values is that distances do not satisfy the triangle inequality, hence importance scores cannot be easily computed.
We overcome this hurdle using a method introduced by Varadarajan and Xiao [VX12a] for projective clustering (where the triangle inequality similarly does not hold). They reduce the importancescore computation to the construction of a coreset for k-CENTER objective; this method is quite different from earlier approaches, e.g. [FL11, VX12b, FSS20, BJKW21], and yields a coreset for k-MEANS whose size depends linearly on log n and of course on the size of the k-CENTER coreset. (Mathematically, this arises from the sum of all importance scores.) We make use of this reduction, and thus focus on constructing (efficiently) a small coreset for k-CENTER with missing values.
An immediate difficulty is how to deal with the missing values. We show that it is possible to find a collection of subsets of coordinates I (so each I ∈ I is a subset of [d]), such that if we construct k-CENTER coresets SI on the data set “restricted” to each I ∈ I, then the union of these SI ’s is a k-CENTER coreset for the original data set with missing values. Crucially, we ensure that each “restricted” data set does not contain any missing value, so that it is possible to use a classical coreset construction for k-CENTER. Finally, we show in a technical lemma how to find a collection as necessary of size |I| ≤ (jk)O(min{j,k}). Since a “restricted” data set does not contain any missing values, we can use a classical k-CENTER coreset construction, and a standard construction has size O(k −d) [AP02], which is known to be tight. We bypass this −d limitation by observing that actually Õ(1)-coreset for k-CENTER suffices, even though the final coreset error is . We observe that an Õ(1)-coreset can be constructed using a variant of Gonzalez’s algorithm [Gon85].
To implement Gonzalez’s algorithm, a key step is to find the furthest neighbor of a given subset of at most O(k) points, and a naive implementation of this runs in linear time, which overall yields a quadratic-time coreset construction, because the aforementioned reduction of [VX12a] actually requires Θ(n/k) successive runs of Gonzalez’s algorithm. To resolve this issue, we propose a fully-dynamic implementation of Gonzalez’s algorithm so that a furthest-point query is answered in time poly(k log n), and the point-set is updated between successive runs instead of constructed from scratch. Our dynamic algorithm is based on a random-projection method that was proposed for furthest-point queries in the streaming setting [Ind03]. Specifically, we project the (restricted) data
set onto several random directions, and on each projected (one-dimensional) data set we apply a data structure for intervals.

heading: 1.3 Additional Related Work
text: Coresets for k-MEANS and k-MEDIAN clustering have been studied extensively for two decades, and we only list a few notable results. The first strong coresets for Euclidean k-MEANS and kMEDIAN were given in [HM04]. In the last decade, most work on coresets for clustering follows the importance sampling framework initiated in [LS10, FL11]. In Euclidean space, recent work showed that coresets for k-MEANS and k-MEDIAN clustering can have size that is independent of the Euclidean dimension [FSS20, SW18, HV20]. Beyond Euclidean space, coresets of size independent of the data-set size were constructed also for many important metric spaces [HJLW18, BJKW21, CASS21]. A more comprehensive overview can be found in recent surveys [Phi17, Fel20].
Recently, attention was given also to non-traditional settings of coresets for clustering, including coresets for Gaussian mixture models (GMM) [LFKF17, FKW19]; simultaneous coresets for a large family of cost functions that include both k-MEDIAN and k-CENTER [BJKW19]; and coresets for clustering under fairness constraints [HJV19]. Also considered were settings that capture uncertainty, for example when each point is only known to lie in a line (i.e., clustering lines) [MF19], and when each point comes from a finite set (i.e., clustering point sets) [JTMF20].

heading: 2 Preliminaries
text: We represent a data point as a vector in (R ∪ {?})d, and a coordinate takes “?” if and only if it is missing. Let Rd? be a shorthand for (R ∪ {?})d. Throughout, we consider a data set X ⊂ Rd? . The distance is evaluated only on the coordinates that are present in both x, y, i.e.,
∀x, y ∈ Rd? , dist(x, y) := √ ∑ i:xi,yi 6=? (xi − yi)2.
For x ∈ Rd? , we denote the set of coordinates that are not missing by Ix := {i : xi 6=?}. For integer m ≥ 1, let [m] := {1, . . . ,m}. For two points p, q ∈ Rd? and an index set I ⊆ Ip ∩ Iq , we define the I-induced distance to be distI(p, q) := √∑ i∈I(pi − qi)2. A point x ∈ Rd? is called a j-point if it has at most j missing coordinates, i.e., |Ix| ≥ d− j. We consider a general k-clustering problem called (k, z)-clustering, which asks to minimize the following objective function. This objective function (and problem) is also called k-MEDIAN when z = 1 and k-MEANS when z = 2.
Definition 2.1 ((k, z)-CLUSTERING). For data set X ⊂ Rd? and a center set C ⊂ Rd containing k (usual) points, let
costz(X,C) := ∑ x∈X distz(x,C).
Definition 2.2 ( -Coreset for (k, z)-CLUSTERING). For data set X ⊂ Rd? , we say a weighted set S ⊆ X with weight function w : S → R+ is an -coreset for (k, z)-CLUSTERING, if
∀C ⊂ Rd, |C| = k, ∑ x∈S w(x) · distz(x,C) ∈ (1± ) · costz(X,C).

heading: 3 Coresets
text: Theorem 3.1. There is an algorithm that, given as input a data set X ⊂ Rd? of size n = |X| consisting of j-points and parameters k, z ≥ 1 and 0 < < 1/2, constructs with constant probability an -coreset of size m = Õ ( zz · (j+k) j+k+1
jjkk−z−2 · −2(d log n) z+22
) for (k, z)-CLUSTERING of X , and
runs in time Õ ( (j+k)j+k+1
jjkk−2 · nd+m
) .
Theorem 3.1 is the main theorem of this paper, and we only present a sketch of the proof in this section due to the space limitation. Please see the full version for a more detailed and self-contained proof, as well as a complete description of our algorithm. We remark that (j+k) j+k
jjkk ≤ (jk)O(min(j,k))
which is used in the statement of Theorem 1.1.
As mentioned in Section 1, we use importance sampling method which is a well-known technique for constructing coresets [FL11, FSS20]. A key step is to compute for every data point x ∈ X an importance score σx ≥ 0 that estimates its maximum relative contribution to any solution. The computation of {σx} is standard in metric spaces, see e.g. [FSS20, BJKW21], but this is not applicable for us because distances with missing values do not satisfy the triangle inequality. Hence, we employ an alternative approach proposed by Varadarajan and Xiao [VX12a, Lemma 3.1], which reduces the computation of {σx} to finding coresets for k-CENTER. This coreset concept, adapted to our setting, is defined as follows.
Definition 3.1. An α-coreset for k-CENTER of a data set X ⊂ Rd? is a subset Y ⊆ X such that
∀C ⊂ Rd, |C| = k, max x∈X dist(x,C) ≤ α ·max y∈Y dist(y, C).
We focus on an efficient construction of an Õ(1)-coreset for k-CENTER. The main concern is that the reduction in [VX12a, Lemma 3.1] requires constructing a k-CENTER coreset for multiple data sets. Fortunately, these data sets are related — each data set is a subset of the previous one — and thus to execute the reduction in near-linear time, we need a k-CENTER coreset construction that supports efficient point deletions. Such a dynamic coreset for k-CENTER with missing values is our main technical contribution. We stated it next, and outline its proof in Section 3.1.
Lemma 3.2. There is a randomized dynamic algorithm with the following guarantees. The input is a dynamic set X ⊂ Rd? of j-points, such that X undergoes q adaptive updates (point insertions and deletions) and the points ever added are fixed in advance (non-adaptively). The algorithm maintains in time Õ ( (j+k)j+k+1
jjkk · (j + k log q)(d+ k2 log q)
) per update, a subset Y ⊆ X of size
|Y | ≤ O ( (j+k)j+k+1
jjkk−1 · log d
) such that with constant probability, Y is an O(k √ d log q)-coreset for
k-CENTER on X after every update.
3.1 Proof of Lemma 3.2: Dynamic Õ(1)-Coresets for k-Center Clustering
As mentioned, the high level idea is to identify a collection I of subsets of coordinates (so each I ∈ I satisfies I ⊆ [d]), construct for each Ii ∈ I an α-coreset Yi (for α determined later) for k-CENTER on the data set X with coordinates restricted on Ii (as defined below), and then their union ⋃ i Yi
would be the overall α √ d-coreset for k-CENTER on X .
Definition 3.2. For a point p ∈ Rd? and a subset I ⊆ Ip, define p|I ∈ RI in the obvious way, by selecting the coordinates {pi}i∈I . Define the I-restricted data set to beX|I := {p|I : p ∈ X, I ⊆ Ip}. Since each vector in X|I arises from a specific vector in X , a subset Y ⊆ X|I corresponds to a specific subset of X , and we shall denote this subset by Y −1.
We observe that X|I ⊂ R|I|, namely, has no missing values (because of the condition I ⊆ Ip). Thus, the metric space on the restricted data set is an ordinary metric space that satisfies the triangle inequality, and so our goal is reduced to constructing k-CENTER coresets for this ordinary setting. However, another key step is to identify a small collection I such that the union of the coresets restricted on I yields a coreset. To this end, we consider the so-called (j, k, d)-family of coordinates as in Definition 3.3. We show in Lemma 3.3 that such a family guarantees the correctness of the coreset, and in Lemma 3.4 that a small family exists and moreover can be constructed efficiently.
Definition 3.3. A collection of sets I ⊂ 2[d] is called a (j, k, d)-family if for every two disjoint subsets J,K ⊂ [d], |J | = j, |K| = k, the family includes I ∈ I that misses J and contains K, i.e., I ∩ J = ∅ and K ⊂ I . Lemma 3.3. Suppose I is a (j, k, d)-family Let X ⊆ Rd? be a set of j-points, and for every I ∈ I, let YI be an α-coreset for k-CENTER on X|I . Then ∪I∈IY −1I is an α √ d-coreset for k-Center on X .
Proof. It suffices to show that for any center set C = {c1, . . . , ck} ⊆ Rd with k points and x ∈ X , if dist(x,C) ≥ r for some r ≥ 0, then we can find a coreset point y ∈ ∪I∈IY −1I such that dist(y, C) ≥ r
α √ d .
For i ∈ [k], let ti ∈ arg maxt∈Ix |xt − cit|, i.e., ti is the index of coordinate that contributes the most in distance dist(x, ci), so |xti − citi | ≥ r√ d
. Let K be any k-subset such that K ⊆ Ix and {t1, . . . , tk} ⊆ K. Since I is a (j, k, d)-family and |Ix| ≥ d− j, by definition, there exists an I ⊆ I such that K ⊆ I ⊆ Ix. We note that
dist(x|I , C|I) = distI(x,C) = min i∈[k]
distI(x, c i) ≥ min
i∈[k] distK(x, c i) ≥ min i∈[k] |xti − citi | ≥ r√ d .
Since I ⊆ Ix, we know that x|I ∈ X|I . As YI is an α-coreset for X|I , we know that there exists y ∈ Y −1I such that
dist(y, C) ≥ distI(y, C) = dist(y|I , C|I) ≥ dist(x|I , C|I) α ≥ r α √ d .
Lemma 3.4 asserts the existence of a small (j, k, d)-family. We remark that this combinatorial structure has been employed in designing fault-tolerant data structures and algorithms (cf. [DK11, DGR21, KP21]), and they obtained similar bounds although their context and language is different.
Lemma 3.4. There is a (j, k, d)-family I of size |I| = O ( (j+k)j+k+1
jjkk log d
) . Moreover, there is a
randomized algorithm that constructs such I in time O(d · |I|) with probability at least 1− 1 dj+k .
k-CENTER coreset for restricted data set via Gonzalez’s algorithm. Finally, the k-CENTER coreset for the restricted data set on each I ∈ I is constructed using an approximate version of Gonzalez’s algorithm [Gon85], where we first pick an arbitrary data point as the initial coreset, and in every iteration an Õ(1)-approximate furthest neighbor in the dataset is picked into the coreset, and we do this for k times. We show that the resulting coreset, consisting of k + 1 points, is an Õ(1)-coreset for k-CENTER. The assumption that the input forms a metric space is crucial, and this is guaranteed since we always run on a restricted data set that satisfies the triangle inequality.
Dynamic implementation of Gonzalez’s algorithm. To make this k-CENTER coreset construction dynamic, we adapt the random projection technique to Gonzalez’s algorithm. We project the (restricted) point set onto several random lines and use a one-dimensional data structure to construct k-CENTER coreset for each of these (projected) one-dimensional data set. Note that the key step in Gonzalez’s algorithm is finding a furthest neighbor, and we can show that our projection method yields an O(k √ log n)-approximation of the furthest neighbor with high probability.
Lemma 3.5. Let A ⊂ Rd, |A| = n, δ > 0 and integer k ≥ 1. Let V be a collection of O(k log n+ log 1δ ) random vectors, each drawn independently fromN (0, Id). Then with probability at least 1− δ, for every P ⊆ A and Q ⊆ A, |Q| ≤ k, if p? is a furthest point in P from Q, and for every v ∈ V we let 〈pv, v〉 be a furthest point in 〈P, v〉 from 〈Q, v〉, then
dist(p?, Q) ≤ O(k √
log n) ·max v∈V dist(pv, Q),
where we denote 〈X, v〉 := {〈x, v〉 : x ∈ X}.
For each one-dimensional line, we use a balanced search tree structure to support the point update and the furthest neighbor query. All operations can be done in O(k logm) time where m is the number of currently inserted elements. This combining with the above lemmas implies Lemma 3.2.

heading: 4 Experiments
text: We implement our proposed coreset construction algorithm, and we evaluate its performance on real and synthetic datasets. We focus on k-MEANS with missing values, and we examine the speedup
for a Lloyd’s-style heuristic. In addition to measuring the absolute performance of our coreset, we also compare it with a) uniform sampling baseline, which is a naive way to construct coresets, and b) an imputation-based baseline where missing values are filled in by random values and then a standard importance-sampling coreset construction (cf. [FL11]) is run on top of it. We implement the algorithms using C++ 11, on a laptop with Intel i5-8350U CPU and 8GB RAM.
Datasets. We run our experiments on three real datasets and one synthetic dataset. Below, we briefly describe how we process and choose the attributes of the dataset, and the parameters of the datasets after processing are summarized in Table 1.
1. Russian housing [Rus17] is a dataset on Russian house market. We pick four main numerical attributes of the houses which are the full area, the live area, the kitchen area and the price, and the price attribute is divided by 105 so as it lies in the similar range of other attributes. Three columns regarding area contain missing values, and the price column doesn’t contain any missing value. 2. KDDCup 2009 [KDD09] is a dataset on customer relationship prediction. We pick 31 numerical attributes that have similar magnitudes. Each column contains missing values. 3. Vertical farming [Sam21] is a dataset about cubes which are used for advanced vertical farming. We include all of four numerical attributes of the dataset. Each column contains missing values. 4. Synthetic dataset. We generate a large synthetic dataset to validate our algorithm’s scalability. Data points are randomly generated so that 97% of them are in a square and 3% of them are far away from the square. After that, we delete 25% of attributes at random. We remark that the 3% far away points is to make the dataset less uniform which prevents it from being trivial for clustering.
Implementation notes. In our experiments, we follow a standard practice of fixing coreset size in each experiment (cf. [BBH+20, JTMF20]). Recall that when computing the importance score, our algorithm chooses a family I of subsets of coordinates and work on each restricted data set X|I for I ∈ I (see Section 3.1). For a fixed size coreset, the family size |I| is a parameter that needs to be optimized. In Figure 1, we plot the empirical error (defined in (1), Section 4.1) for the Russian housing dataset with respect to the family size |I|. Although Lemma 3.4 gives a theoretical upper bound on |I| but our experiments suggest that a much smaller size |I| = 20 is optimal in this case.

heading: 4.1 Accuracy of Coresets
text: We evaluate the accuracy versus size tradeoff of our coresets. Since the coreset should preserve the clustering cost for all centers, we evaluate the accuracy by testing the empirical error on a selected set of centers C. Namely, for a data set X , a coreset D ⊆ X and a collection of center sets C, we define the empirical error of D as
err(D) = max C∈C |cost(D,C)− cost(X,C)| cost(X,C) . (1)
We use a randomly selected collection of centers C that consists of 100 randomly generated k-subset C ⊂ Rd. Since both the evaluation method and the algorithm has randomness, we run the experiment for T = 103 times with independent random bits and report the average empirical error to make it stable. We choose 20 different coreset sizes from 200 to 9700 in a step size of 500, and report the corresponding average empirical error.
Results. We report the size versus accuracy tradeoff of our coreset for all four datasets in Figure 2, and record the standard deviation in Figure 3. We compare these results against the abovementioned uniform sampling and imputation baseline. As can be seen from the figures, the accuracy of our coreset improves when the size increases, and we achieve 5%-20% error using only 2000 coreset points (which is within 0.5%− 5% of the datasets). This 5%-20% error is likely to be enough for practical use, since practical algorithms for k-MEANS are approximation algorithms anyway. Our coresets generally outperform both the uniform sampling and imputation baselines on almost every coreset sample size, and the advantage is more significant when the coreset size is relatively small. Moreover, our coresets have a much lower variance.

heading: 4.2 Speedup of Lloyd’s-style Heuristic
text: Coresets often help to speed up existing approximation algorithms. Before our work, the only algorithm for k-MEANS with provable guarantees for multiple missing values was [EFG+21]. Unfortunately, [EFG+21] is not practical even when combined with coresets, since it contains several
enumeration procedures that require Θ(exp(poly( −1jk))) time. We consider a variant of Lloyd’s heuristic [Llo82] that is adapted to the missing-value setting, and we evaluate its speedup with coresets. The algorithm is essentially the same as the original Lloyd’s algorithm, except that the distance as well as the optimal 1-mean for a cluster (which can be computed optimally in O(d|P |) for a cluster P [EFG+21]), is computed differently. We show that our coreset can significantly accelerate this algorithm. In particular, we run the modified Lloyd’s heuristic directly on the original dataset, and take its running time and objective value as the comparison reference. Then we run this modified Lloyd’s heuristic again, but on top of our coreset and the uniform sampling baseline respectively, and we compare both the speedup and the relative error3 against the reference. The experiments are run on the Russian housing data set where the number of iterations of the modified Lloyd’s is set to T = 5000 and the number of clusters is set to a small value k = 3 so as the heuristic is likely to find a local minimum faster. Again, to obtain a stable result, we run the experiments for 40 times with independent random bits and report the average relative errors and running time.
Results. The relative error with respect to varying coreset sizes can be found in Figure 4a. We can see that the relative error of Lloyd’s algorithm running on our coreset is consistently low, while the uniform sampling baseline has several times higher error and the error does not seem to improve even when improving the size. We note that relative errors for both our coreset and uniform sampling are significantly lower than that we observe from the empirical error in Figure 2a. In fact, they are not necessarily comparable since the empirical error in Figure 2a is always evaluated on a same center, while what we compare in Figure 4a is the center sets found by the modified Lloyd’s running on different data sets. This also helps to explain why improving the size of uniform sampling may not result in a better solution, since as shown in Figure 2a, uniform sampling has a large empirical error (around 50%), so a good solution for the uniform sample may not be a good solution for the original data set.
The running time of the modified Lloyd’s on top of our coresets can be found in Figure 4b, and the running time of Lloyd’s on the original dataset is 22.9s (which is not drawn on the figure). To make a fair comparison, we also take the coreset construction time into account. Note that coreset size is not a dominating factor in the running time of coreset construction, since the majority of time is
3For x ∈ R+, the relative error of x against a reference x? > 0 is defined as |x−x ?|
x? .
spent on computing the importance scores and the coreset size only affects the number of samples. A coreset of size only 1000 can achieve < 1% error, and the running time of constructing the coreset and applying Lloyd’s on top of it are 3s and 0.8s, respectively, which offers more than 5 times of speedup. We remark that our experiments only demonstrate the speedup in a single-machine scenario, and the speedup will increase in the parallel or distributed setting.

heading: 5 Conclusion
text: Our coreset construction builds upon the sensitivity-sampling method (cf. [FL11]). However, a central technical challenge is that the standard method to compute the sensitivity scores breaks, because distances between points with missing values do not satisfy the triangle inequality. We overcome this using another known method, of [VX12a], that requires a coreset for k-CENTER. Our main innovation is a near-linear time algorithm that computes an O(1)-approximate k-CENTER coreset for points with missing values. To this end, we need the following key steps, which constitute our main technical contribution.
• We reduce the k-CENTER coreset construction with missing values, to the construction of traditional k-CENTER coresets (i.e., without missing values) on a series of instances. These instances are built by restricting data points with missing values to a carefully-chosen collection of subspaces. The guarantee needed from this collection is a certain combinatorial structure, and we indeed prove it exists.
• The method of Varadarajan and Xiao executes the k-CENTER coreset algorithm many times, and overall takes quadratic time. To improve the running time, we design an efficient dynamic algorithm for the well-known Gonzales’ algorithm (which computes an O(1)approximate k-CENTER coreset). The main idea in this dynamic algorithm is to project the data points onto (data-oblivious) random 1D lines, and build on each line a dynamic data structure that supports furthest-neighbor queries (in 1D).
Finally, we implemented our algorithm and the experiments indicate that our algorithm is efficient and accurate enough to be potentially applicable in practice.
Future directions. As an immediate follow-up, one could try to improve our coreset size, e.g., removing the dependence in log n. Our input can be viewed as axis-parallel affine-subspaces. Hence, another an interesting direction is to obtain coresets for the more general setting where the input consists of general affine-subspaces.
Potential negative societal impacts. Our paper focuses on computational issues (improving time and space) of known clustering tasks. Clustering methods in general have potential issues with fairness and privacy, which applies also to our work, but our research is not expected to introduce new negative societal impact beyond what is already known.
Acknowledgments and Disclosure of Funding
The majority of this work was done when Shaofeng Jiang was at Aalto University. This work is partially supported by ONR Award N00014-18-1-2364, by the Israel Science Foundation grant #1086/18, by a Weizmann-UK Making Connections Grant, and by a Minerva Foundation grant.

Review summaries:
The paper provides the first coreset for clustering tasks, where each of the input d-dimensional points might contain up to j missing values. The coreset construction time is near-linear. The main technical result of the paper is a reduction from the above problem, to the problem of computing a dynamic coreset for the k-center problem with missing values. The authors utilize the sensitivity sampling framework for this task. Using this coreset, the authors improve the running time of an existing PTAS for the problem, from quadratic, to near-linear time. The authors provide experimental results which demonstrate the effectiveness of the coreset in practice, as compared to a uniform random sample.

The paper describes the first algorithms for coresets for k-means clustering with missing data coordinates for some points in d-dimensions, with provably guarantees.  The algorithm runs in linear time O(nd) times factor depending only on j (most possible missing values per coordinate) and k (number of allowed clusters).  A SODA 2021 paper provides a PTAS for this problem, but it has worse dependence on j, k, epsilon and has a base cost of O(n^2 d), so quadratic in n.  This paper demonstrates their approach is implementable, and practical.