summaries: This paper considers the problem of finding good initial centers for the fundamental problem of $k$-median clustering using a randomized embedding of the original metric into a tree metric.  After setting the initial centers, a standard local search algorithm is applied to produce an improved solution.  This is explored in both the standard context of $k$-median clustering, as well as in the relevant context of differentially private clustering.  In the latter setting, the goal is to minimize the amount of additive error introduced by the algorithm subject to being $\epsilon$-differentially private.  An extension to $k$-means is given in the appendix.  In the standard setting of $k$-median clustering, the main theoretical result is an initialization algorithm which is an $O(\log(\min(\Delta,k)))$-approximation to the optimal k-median clustering.  This is an improvement over k-median++ (which gives $O(\log k)$) when $\Delta$ is small, e.g. for $\Delta = O(d)$ and $d$ is small.  Using this as a seed results for a local search method results in an $O(1)$-approximation overall.  At a high level, their algorithm first constructs an embedding of the original metric into a hierarchically well-separated tree (HST).  From there, the initialization can be seen as finding an $O(1)$-approximate solution on the HST efficiently.  The overall guarantee follows from standard results about HST's  In  the differentially private setting, the main result is a similar guarantee on the quality of the initial solution and also a bound on the quality of the final solution when using a known private local search algorithm.  The quality of the final solution has $O(1)$-multiplicative error and $O(\epsilon^{-1}k^2\Delta\log(n)\log\log(n))$ additive error.  This is an improvement over the additive error of $O(\epsilon^{-1}k^2\Delta\log^2(n)$ due to Gupta et al. 2010.  The number of local search iterations is also improved from $O(k\log n)$ to $O(k\log\log n)$.  The main idea for the initialization is similar to the standard setting, but here they use the structure of the HST to ensure the initial solution is private by injecting a different amount of noise at each level of the tree.  An empirical study is done on a class of synthetic graphs as well as the MNIST dataset.  For the synthetic graphs, the metric space is given by the weighted shortest path distance in each graph, while for the MNIST dataset the metric is given by either $\ell_1$ or $\ell_2$.  The authors compare both the initial costs and the final costs (after running a local search method) for several initialization methods in both the standard and differentially private settings.  The main observation is that the proposed initialization methods tend to have better initial cost and the proposed differentially private method often outperforms the other methods in both initial cost and final cost.     The paper proposes a new initialization scheme for the k-median problem on graph input (or general metric spaces) using metric embedding tree structure. The paper proposes an algorithm that finds initialization of good centers using HST that gets an approximation factor of O(log min{k,d}) if the data is in Euclidean space  where d is the number of dimensions. Then, the paper studies clustering with differential privacy guarantee and hows that the initialization method could be adapted to give a slightly stronger muliplicative and additive errors. The work complemented these theoretical findings with experiments and show that the proposed initialization imporves the performance of k-median++ initialization. This paper introduces a new initialization scheme for the k-medians clustering problem in the general metric space setting. This is based on the construction of metric embeddings via 2-HST’s (Hierarchically well-separated trees). The authors also extend this to the differential privacy (DP) setting. They prove approximation guarantees in both the non-DP and DP settings, improving upon the literature. Finally, they empirically validate algorithms against a number of baselines with both real world and synthetic datasets for multiple metrics. The paper suggests an algorithm for the metric k-median problem using ideas from Metric embedding theory. The suggested use of the algorithm is as an initialization routine for the local search based algorithm for k-median. The differentially private version of the algorithm is also given along with bounds on k-median approximation factor. Experiments are conducted over datasets such as MNIST and results compared against the k-means++ algorithm (a popular initialisation algorithm).
abstractText: In clustering algorithms, the choice of initial centers is crucial for the quality of the learned clusters. We propose a new initialization scheme for the k-median problem in the general metric space (e.g., discrete space induced by graphs), based on the construction of metric embedding tree structure of the data. From the tree, we propose a novel and efficient search algorithm for good initial centers that can be used subsequently for the local search algorithm. Our method, named the HST initialization, can also be easily extended to the setting of differential privacy (DP) to generate private initial centers. Theoretically, the initial centers from HST initialization can achieve lower error than those from another popular initialization method, k-median++, in the non-DP setting. Moreover, with privacy constraint, we show that the error of applying DP local search followed by our private HST initialization improves previous results, and approaches the known lower bound within a small factor. Empirically, experiments are conducted to demonstrate the effectiveness of our methods.
references.authors: Ameer Ahmed Abbasi Mohamed F. Younis
references.title: A survey on clustering algorithms for wireless sensor networks
references.venue: Comput. Commun.,
references.authors: David Arthur Sergei Vassilvitskii
references.title: k-means++: the advantages of careful seeding
references.venue: In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms,
references.authors: Vijay Arya Naveen Garg Rohit Khandekar Adam Meyerson Kamesh Munagala Vinayaka Pandit
references.title: Local search heuristics for k-median and facility location problems
references.venue: SIAM J. Comput.,
references.authors: Olivier Bachem Mario Lucic S. Hamed Hassani Andreas Krause
references.title: Approximate k-means++ in sublinear time
references.venue: Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, February 12-17,
references.authors: Bahman Bahmani Benjamin Moseley Andrea Vattani Ravi Kumar Sergei Vassilvitskii
references.title: Scalable k-means++
references.venue: Proc. VLDB Endow.,
references.authors: Maria-Florina Balcan Steven Ehrlich Yingyu Liang
references.title: Distributed k-means and k-median clustering on general communication topologies. In Advances in Neural Information Processing Systems
references.venue: Annual Conference on Neural Information Processing Systems
references.authors: Maria-Florina Balcan Travis Dick Yingyu Liang Wenlong Mou Hongyang Zhang
references.title: Differentially private clustering in high-dimensional euclidean spaces
references.venue: In Proceedings of the 34th International Conference on Machine Learning,
references.authors: Arindam Banerjee Srujana Merugu Inderjit S. Dhillon Joydeep Ghosh
references.title: Clustering with bregman divergences
references.venue: J. Mach. Learn. Res.,
references.authors: Yair Bartal
references.title: Probabilistic approximations of metric spaces and its algorithmic applications
references.venue: Annual Symposium on Foundations of Computer Science, FOCS ’96,
references.authors: John E Beasley
references.title: Or-library: distributing test problems by electronic mail
references.venue: Journal of the operational research society,
references.authors: Inderjit S. Dhillon Dharmendra S. Modha
references.title: Concept decompositions for large sparse text data using clustering
references.venue: Mach. Learn.,
references.authors: Cynthia Dwork
references.title: Differential privacy. In Automata, Languages and Programming, 33rd International Colloquium, ICALP 2006, Venice, Italy
references.venue: July 10-14,
references.authors: Jittat Fakcharoenphol Satish Rao Kunal Talwar
references.title: A tight bound on approximating arbitrary metrics by tree metrics
references.venue: J. Comput. Syst. Sci.,
references.authors: D. Feldman C. Xiang R. Zhu D. Rus
references.title: Coresets for differentially private k-means clustering and applications to privacy in mobile sensor networks
references.venue: 16th ACM/IEEE International Conference on Information Processing in Sensor Networks (IPSN),
references.authors: Dan Feldman Amos Fiat Haim Kaplan Kobbi Nissim
references.title: Private coresets
references.venue: In Proceedings of the 41st Annual ACM Symposium on Theory of Computing,
references.authors: Anupam Gupta Katrina Ligett Frank McSherry Aaron Roth Kunal Talwar
references.title: Differentially private combinatorial optimization
references.venue: In Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms,
references.authors: Zhiyi Huang Jinyan Liu
references.title: Optimal differentially private algorithms for k-means clustering
references.venue: In Proceedings of the 37th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems,
references.authors: Tapas Kanungo David M. Mount Nathan S. Netanyahu Christine D. Piatko Ruth Silverman Angela Y. Wu
references.title: A local search approximation algorithm for k-means clustering
references.venue: In Proceedings of the 18th Annual Symposium on Computational Geometry,
references.authors: Leon Kaufman Marc Vanden Eede Pierre Hansen
references.title: A plant and warehouse location problem
references.venue: Journal of the Operational Research Society,
references.authors: Silvio Lattanzi Christian Sohler
references.title: A better k-means++ algorithm via local search
references.venue: In Proceedings of the 36th International Conference on Machine Learning,
references.authors: Yann LeCun Léon Bottou Yoshua Bengio Patrick Haffner
references.title: Gradient-based learning applied to document recognition
references.venue: Proceedings of the IEEE,
references.authors: S. Lloyd
references.title: Least squares quantization in pcm
references.venue: IEEE Transactions on Information Theory,
references.authors: Konstantin Makarychev Yury Makarychev Ilya P. Razenshteyn
references.title: Performance of johnsonlindenstrauss transform for k-means and k-medians clustering
references.venue: In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing,
references.authors: Frank McSherry Kunal Talwar
references.title: Mechanism design via differential privacy
references.venue: Annual IEEE Symposium on Foundations of Computer Science (FOCS 2007), October 20-23,
references.authors: Richard Nock Raphaël Canyasse Roksana Boreli Frank Nielsen
references.title: k-variates++: more pluses in the k-means++
references.venue: In Proceedings of the 33nd International Conference on Machine Learning,
references.authors: Girish Punj David W Stewart
references.title: Cluster analysis in marketing research: Review and suggestions for application
references.venue: Journal of marketing research,
references.authors: Rahul Shah
references.title: Faster algorithms for k-median problem on trees with smaller heights
references.venue: Technical report,
references.authors: Uri Stemmer Haim Kaplan
references.title: Differentially private k-means with constant multiplicative error
references.venue: Advances in Neural Information Processing Systems
references.authors: Keisuke Todo Atsuyoshi Nakamura Mineichi Kudo
references.title: A fast approximate algorithm for kmedian problem on a graph
references.venue: 15th International Workshop on Mining and Learning with Graphs,
references.authors: Gupta
references.title: Assume the solution to the optimal utility is unique. For any output o ∈ O of 24 -DP exponential mechanism on dataset D, it holds for ∀t
sections.heading: 1 INTRODUCTION
sections.text: Clustering is an important problem in unsupervised learning that has been widely used in social network analysis, sensor networks, etc. (Punj & Stewart, 1983; Dhillon & Modha, 2001; Banerjee et al., 2005; Abbasi & Younis, 2007). The goal of clustering is to partition a set of data points into clusters such that items in the same cluster are expected to be similar, while items in different clusters should be different. This is concretely measured by the sum of distances (or squared distances) between each point to its nearest cluster center. One conventional notion to evaluate a clustering algorithms is: with high probability,
  cost(C,D) ≤ γOPTk(D) + ξ, where C is the centers output by the algorithm and cost(C,D) is some cost function defined for C on the dataset D. OPTk(D) is the cost of optimal (oracle) k-median solution on D. When everything is clear from context, we will use OPT for short. Here, γ is called multiplicative error and ξ is called additive error. Alternatively, we may also use the notion of expected cost.
  Two popularly studied clustering problems are 1) the k-median problem, and 2) the k-means problem. The origin of k-median dates back to the 1970’s (e.g., Kaufman et al. (1977)), where one tries to find the best location of facilities that minimizes the cost measured by the distance between clients and facilities. Formally, given a set of points D and a distance measure, the goal is to find k center points minimizing the sum of absolute distances of each sample point to its nearest center. In k-means, the objective is to minimize the sum of squared distances instead. There are two popular clustering algorithms. One heuristic is the Lloyd’s algorithm (Lloyd, 1982), which is built upon an iterative distortion minimization approach. In most cases, this method can only be applied to numerical data, typically in the Euclidean space. Clustering in general metric spaces (discrete spaces) is also important and useful when dealing with, for example, the graph data, where Lloyd’s method is no longer applicable. A more broadly applicable approach, the local search method (Kanungo et al., 2002; Arya et al., 2004), has also been widely studied. It iteratively finds the optimal swap between the center set and non-center data points to keep lowering the cost. Local search can achieve a constant approximation ratio (γ = O(1)) to the optimal solution for k-median (Arya et al., 2004).
  Initialization of cluster centers. It is well-known that the performance of clustering can be highly sensitive to initialization. If clustering starts with good initial centers (i.e., with small approximation error), the algorithm may use fewer iterations to find a better solution. The k-median++
  algorithm (Arthur & Vassilvitskii, 2007) iteratively selects k data points as initial centers, favoring distant points in a probabilistic way (see Appendix A for more details). Intuitively, the initial centers tend to be well spread over the data points (i.e., over different clusters). The produced initial center is proved to have O(log k) multiplicative error. Followup works of k-means++ further improved its efficiency and scalability, e.g., Bahmani et al. (2012); Bachem et al. (2016); Lattanzi & Sohler (2019). In this work, we propose a new initialization framework based on metric embedding methods, with comparable approximation error and running time as k-median++. Importantly, our initialization scheme can be conveniently combined with the notion of differential privacy.
  Clustering with Differential Privacy. The concept of differential privacy (Dwork, 2006; McSherry & Talwar, 2007) has been popular to rigorously define and resolve the problem of keeping useful information for model learning, while protecting privacy for each individual. Private k-means problem has been widely studied, e.g., Feldman et al. (2009); Nock et al. (2016); Feldman et al. (2017), mostly in the Euclidean space. The paper (Balcan et al., 2017) considered identifying a good candidate set (in a private manner) of centers before applying private local search, which yieldsO(log3 n) multiplicative error and O((k2 + d) log5 n) additive error. Later on, the Euclidean k-means errors are further improved to γ = O(1) and ξ = O(k1.01 ·d0.51+k1.5) by Stemmer & Kaplan (2018), with more advanced candidate set selection. Huang & Liu (2018) gave an optimal algorithm in terms of minimizing Wasserstein distance under some data separability condition.
  For private k-median clustering, Feldman et al. (2009) considered the problem in high dimensional Euclidean space. However, it is rather difficult to extend their analysis to more general metrics in discrete spaces (e.g., on graphs). The strategy of Balcan et al. (2017) to form a candidate center set could as well be adopted to k-median, which leads to O(log3/2 n) multiplicative error and O((k2+d) log3 n) additive error in high dimensional Euclidean space. Gupta et al. (2010) proposed a private method for the classical local search heuristic in discrete space, which applies to both kmedians and k-means. In order to cast privacy on each swapping step, the authors applied the exponential mechanism of McSherry & Talwar (2007). Their method produced an -differentially private solution with cost 6OPT +O(4k2 log2 n/ ), where4 is the diameter of point set. The main contributions of this work include :
  • We introduce Hierarchically Well-Separated Tree (HST), a metric embedding tree structure, to the k-median clustering problem. Once the HST is constructed using the data samples, we provide an efficient sampling strategy to select the initial center set, with an approximation factor O(log min{k,4}) in the non-private setting, which is O(log min{k, d}) when 4 = O(d) (e.g., bounded data). This improves the O(log k) error of k-means/median++ in e.g., the lower dimensional Euclidean space.
  • The main strength of our HST initialization method is that it could be effectively adapted to differential privacy (DP). The so-called DP-HST algorithm is -DP and outputs initial centers withO(log n) multiplicative error andO( −14k2 log2 n) additive error. Moreover, running DP-local search starting from this initialization givesO(1) multiplicative error and O( −14k2(log log n) log n) additive error, which improves previous results towards the well-known lower bound O( −14k log(n/k)) on the additive error (Gupta et al., 2010) within a small O(k log log n) factor. To our knowledge, this is the first center initialization method with differential privacy guarantee and improved error rate in general metrics.
  • We conduct experiments on simulated and real-world datasets to demonstrate the effectiveness of our methods. In both non-private and private settings, our proposed HST-based approach achieves smaller cost at initialization than k-median++, which may also lead to improvements in the final clustering quality. Our private algorithm can also save computational costs by reaching a good solution with fewer iterations.
sections.heading: 2 PRELIMINARIES
sections.text: Definition 2.1 (Differential Privacy (DP) (Dwork, 2006)). If for any two adjacent data sets D and D′ with symmetric difference of size one, for any O ⊂ Range(A), an algorithmA satisfies
  Pr[A(D) ∈ O] ≤ e Pr[A(D′) ∈ O],
  then algorithmA is said to be -differentially private.
  Intuitively, differential privacy requires that after removing any observation, the output of D′ should not be too different from that of the original dataset D. Smaller indicates stronger privacy, which, however, usually sacrifices utility. Thus, one of the central topics in differential privacy literature is to balance the utility-privacy trade-off.
  To achieve DP, it is common to introduce noise to the data or the algorithm output. The Laplace mechanism adds Laplace(η(f)/ ) noise to the output, which is known to achieve -DP. The exponential mechanism is also a tool for many DP algorithms. Let O be the set of feasible outputs. The utility function q : D × O → R is what we aim to maximize. The exponential mechanism outputs an element o ∈ O with probability P [A(D) = o] ∝ exp( q(D,o)2η(q) ), where D is the input dataset and η(f) = sup|D−D′|=1 |f(D) − f(D′)| is the sensitivity of f . Both mechanisms will be used in the design of our proposed DP approach.
  k-Median Clustering. Following the framework of Arya et al. (2004); Gupta et al. (2010), the problem of k-median clustering (DP and non-DP) studied in our paper is stated as below.
  Definition 2.2 (k-median). Given a universe point set U and a metric ρ : U × U → R, let D ⊆ U be a set of demand points. The goal of DP k-median is to pick F ⊆ U with |F | = k to minimize
  Private k-median: cost(F,D) = ∑ v∈D min f∈F ρ(v, f). (1)
  At the same time, the output F is required to be -differentially private to D. For standard nonprivate k-median, we assume U = D and minimize
  k-median: cost(F,U) = ∑ v∈U min f∈F ρ(v, f). (2)
  We may drop “F” and use “cost(D)” or “cost(U)” if there is no risk of ambiguity.
  To better understand the DP clustering problem, we provide an real-world example as follows.
  Example 2.1. Consider U to be the universe of all users in a social network (e.g., Twitter). Each user (account) is public, but also has some private information that can only be seen by the data holder. Let D be users grouped by some information that might be set as private. Suppose a third party plans to collaborate with the most influential users in D for e.g., commercial purposes, thus requesting the cluster centers of D. In this case, we need a strategy to safely release the centers, while protecting the individuals in D from being identified (since the membership of D is private).
  The local search procedure for k-median proposed by Arya et al. (2004) is summarized in Algorithm 1. First we randomly pick k points in U as the initial centers. In each iteration, we search over all x ∈ F and y ∈ U , and do the swap F ← F − {x}+ {y} such that F − {x}+ {y} improves the cost of F the most (if more than factor (1−α/k) where α > 0 is a hyper-parameter). We repeat the procedure until no such swap exists. Arya et al. (2004) showed that the output centers F achieves 5 approximation error to the optimal solution, i.e., cost(F ) ≤ 5OPT .
  Algorithm 1: Local search for k-median clustering (Arya et al., 2004) Input: Data points U , parameter k, constant α Initialization: Randomly select k points from U as initial center set F while ∃ x ∈ F, y ∈ U s.t. cost(F − {x}+ {y}) ≤ (1− α/k)cost(F ) do
  Select (x, y) ∈ Fi × (D \ Fi) with arg minx,y{cost(F − {x}+ {y})} Swap operation: F ← F − {x}+ {y}
  Output: Center set F
sections.heading: 3 INITIALIZATION VIA HIERARCHICALLY WELL-SEPARATED TREE (HST)
sections.text: In this section, we propose our new initialization scheme for k-median clustering, and provide our analysis in the non-private case solving (2). The idea is based on the metric embedding theory. We will start with an introduction to the main tool used in our approach.
sections.heading: 3.1 HIERARCHICALLY WELL-SEPARATED TREE (HST)
sections.text: In this paper, for an L-level tree, we will count levels in descending order down the tree. That is, the root is level L, and its children form level L − 1, and etc. We use hv to denote the level of v, and ni be the number of nodes at level i. The Hierarchically Well-Separated Tree (HST) is based on the padded decompositions of a general metric space in a hierarchical manner (Fakcharoenphol et al., 2004). Let (U, ρ) be a metric space with |U | = n, and we will refer to this metric space without specific clarification. A β–padded decomposition of U is a probabilistic distribution of partitions of U such that the diameter of each cluster Ui ∈ U is at most β, i.e., ρ(u, v) ≤ β, ∀u, v ∈ Ui, i = 1, ..., k. Based on the this, we will adopt the standard 2-HST for our algorithm design.
  Definition 3.1. Assume the smallest distance between u, v ∈ U is 1, and4 = maxu,v∈U ρ(u, v). A 2-Hierarchically Well-Separated Tree (2-HST) is an edge-weighted rooted tree T , such that an edge between any pair of two nodes of level i− 1 and level i has length at most4/2L−i.
  We briefly describe the construction of 2-HST, and place the detailed algorithm in Algorithm 7 in the Appendix A. We first find a padded decomposition PL = {PL,1, ..., PL,nL} of U with parameter β = 4/2. The center of each partition in PL,j serves as a root node in level L. Then, we re-do a padded decomposition for each partition PL,j , to find sub-partitions with diameter β = 4/4, and set the corresponding centers as the nodes in level L − 1, and so on. Each partition at level i is obtained with β = 4/2L−i. This process proceeds until a node has a single point, or a pre-specified tree depth is reached. In Figure 1, we provide an example of L = 3-level 2-HST (left panel), along with its underlying padded decompositions (right panel).
  HST induces a tree metric, which will be related to the general metric in our analysis. It is worth mentioning that, there are polynomial time algorithms for computing an exact k-median solution in the tree metric (Tamir (1996); Shah (2003)). However, the dynamic programming algorithms have high complexity (O(kn2) and O(k2nh) where h is at least O(log n), respectively), making them unsuitable to serve the purpose of fast initialization. Moreover, it is unknown how to apply them effectively to the private case. Hence, next, we propose a novel tree search algorithm that is very efficient, produces a O(1) approximation error in the tree metric, and can be extended to DP easily.
sections.heading: 3.2 HST INITIALIZATION ALGORITHM
sections.text: Let L = log ∆ and suppose T is a level-L 2-HST in (U, ρ), where for simplicity we assume L is an integer. For a node v at level i, we use T (v) to denote the subtree rooted at v. LetNv = |T (v)| be the number of data points in T (v). The sampling strategy for initial centers is presented in Algorithm 2, which has two phases: subtree search and leaf search.
  Subtree search. The first step is to identify the subtrees that contain the k centers. To begin with, k initial centers C1 are picked from T who have the largest score(v) = N(v) · 2hv . This is intuitive, since to get a good clustering, we typically want the ball surrounding each center to include more data points. Next, we do a screening over C1: if there is any ancestor-descendant pair of nodes, we remove the ancestor from C1. If the current size of C1 is smaller than k, we repeat the process until k centers are chosen (we do not re-select nodes in C1 and their ancestors). This way, C1 contains k
  Algorithm 2: HST initialization - NDP (Non-Differentially Private) Input: U ,4, k Initialization: L = log4, C0 = ∅, C1 = ∅ Call Algorithm 7 to build a level-L 2-HST T based on input U for each node v in T do
  Nv ← |U ∩ T (v)| score(v)← Nv · 2hv
  while |C1| < k do Add top (k − |C1|) nodes with highest score in T to C1 for each v ∈ C1 do C1 = C1 \ {v}, if ∃ v′ ∈ C1 that v′ is a descendant of v C0 = FIND-LEAF(T,C1) Output: Initial center set C0 ⊆ U
  Algorithm 3: FIND-LEAF (T,C1) Input: T , C1 Initialization: C0 = ∅ for each node v in C1 do
  while v is not a leaf node do v ← argw max{Nw, w ∈ ch(v)}, where ch(v) denotes the children nodes of v
  Add v to C0 Output: Initial center set C0 ⊆ U
  root nodes of k disjoint subtrees. For any node w with T (w) ∩ C1 = ∅ and a node v ∈ C1, if w is not an ancestor of v, then score(w) ≤ score(v). Leaf search. After we find C1 the set of k subtrees, the next step is to find the center in each subtree using Algorithm 3 (“FIND-LEAF”). We employ a greedy search strategy, by finding the child node with largest score level by level, until a leaf is found. This approach is intuitive since the diameter of the partition ball exponentially decays with the level. Therefore, we are in a sense focusing more and more on the region with higher density (i.e., with more data points).
  The time complexity is given as below, considering the Euclidean space as an example.
  Proposition 3.1 (Complexity). Algorithm 2 takes O(dn log n) time in the Euclidean space. Remark 3.1. The complexity of HST initialization is in general comparable to the O(dnk) of kmedian++. Our algorithm would be faster if k, the number of centers, is larger. See Appendix B for a numerical comparison of the running time.
sections.heading: 3.3 APPROXIMATION ERROR OF HST INITIALIZATION
sections.text: Firstly, we show that the initial center set produced by Algorithm 2 is already a good approximation to the optimal k-median solution. Let ρT (x, y) = dT (x, y) denote the “2-HST metric” between x and y in the 2-HST T , where dT (x, y) is the tree distance between nodes x and y in T . In 2-HST, by Definition 3.1 and 4 = 2L, in the analysis we assume equivalently that the edge weight of the i-th level 2i−1. The crucial step of our analysis is to examine the approximation error in terms of the 2-HST metric, after which the error can be easily adapted to the general metrics by Lemma 3.2.
  Lemma 3.2. (Bartal, 1996). In a metric space (U, ρ) with |U | = n and diameter 4, it holds that E[ρT (x, y)] = O(min{log n, log4})ρ(x, y). In Rd, E[ρT (x, y)] = O(d)ρ(x, y).
  Recall C0, C1 from Algorithm 2. We define costTk (U) = ∑ y∈U min x∈C0 ρT (x, y), (3)
  costTk ′(U,C1) = min
  |F∩T (v)|=1,∀v∈C1 ∑ y∈U min x∈F ρT (x, y), (4)
  OPTTk (U) = min F⊂U,|F |=k ∑ y∈U min x∈F ρT (x, y) ≡ min C′1 costTk ′(U,C ′1). (5)
  For simplicity, we will use costTk ′(U) to denote costTk ′(U,C1). Here, OPTTk (5) is the cost of the global optimal solution with 2-HST metric. The last equality in (5) holds because the optimal centers set can always located in k disjoint subtrees, as each leaf only contain one point. (3) is the k-median cost with 2-HST metric of the output C0 of Algorithm 2. (4) is the oracle cost after the subtrees are chosen. That is, it represents the optimal cost to pick one center from each subtree in C1. Firstly, we bound the approximation error of the subtree search process. Lemma 3.3 (Subtree search). costTk ′(U) ≤ 5OPTTk (U).
  Next, we show that the greedy leaf search (Algorithm 3) only has constant extra multiplicative error. Lemma 3.4 (Leaf search). costTk (U) ≤ 2costTk ′(U).
  Combining Lemma 3.3 and Lemma 3.4, we have the next Theorem. Theorem 3.5 (2-HST error). Running Algorithm 2, we have costTk (U) ≤ 10OPTTk (U).
  Thus, HST-initialization produces an O(1) approximation to OPT in the 2-HST metric. Define costk(U) as (2) for our HST centers, and the optimal cost w.r.t. ρ as
  OPTk(U) = min |F |=k ∑ y∈U min x∈F ρ(x, y). (6)
  We have the following result in the general metric space based on Lemma 3.2. Theorem 3.6. Running Algorithm 2 gives E[costk(U)] = O(min{log n, log4})OPTk(U). Remark 3.2. In the Euclidean space, Makarychev et al. (2019) provedO(log k) random projections suffice for k-median to achieve O(1) error. Thus, if4 = O(d) (e.g., bounded data), by Lemma 3.2, there exists an algorithm with HST initialization that achieves O(log(min{d, k})) error, which is better than O(log k) of k-median++ when d is small.
  NDP-HST Local Search. We can apply standard local search starting from the HST initialization in Algorithm 1. We call it the NDP-HST (“NDP” stands for “Non-Differentially Private”) method, with the following guarantee. Theorem 3.7. NDP-HST achieves O(1) approximation in expected O(k log min{log n, log4}) number of iterations for input in general metric space.
  Before ending this section, we remark that HST initialization and the analysis can be extended to k-means clustering analogously (see Appendix D). In a general metric space, E[costkm(U)] = O(min{log n, log4})2OPTkm(U) where costkm(U) is the optimal k-means cost.
sections.heading: 4 HST INITIALIZATION WITH DIFFERENTIAL PRIVACY
sections.text: In this section, we consider private HST initialization method. Recall in this setting,U is the universe of data points, andD ⊂ U is a demand set that needs to be clustered with privacy. Since U is public, running initialization algorithms on U would preserve the privacy of D. Yet, this might be too expensive, and in many cases one would probably want to incorporate some information about D in the initialization, since D could be a very imbalanced subset of U . For example, D may only contain data points from one cluster, out of tens of clusters in U . In this case, initialization on U is likely to pick initial centers in multiple clusters, which would not be helpful for clustering on D. Next, we show how our HST initialization can be easily combined with differential privacy that at the same time contains information about the demand set D, leading to improved approximation error (Theorem 4.3) and empirical clustering performance (Section 5).
  Again, suppose T is an L = log4-level 2-HST of universe U in a general metric space. Denote Nv = |T (v) ∩ D| for a node point v. Our private HST initialization is similar to the non-private Algorithm 2. To gain privacy, we perturb Nv by adding i.i.d. Laplace noise:
  N̂v = Nv + Lap(2 (L−hv)/ ),
  Algorithm 4: HST initialization with differential privacy Input: U,D,4, k, Build a level-L 2-HST T based on input U for each node v in T do
  Nv ← |D ∩ T (v)| N̂v ← Nv + Lap(2(L−hv)/ ) score(v)← N̂(v) · 2hv
  Based on N̂v , apply the same strategy as Algorithm 2: find C1; C0 = FIND-LEAF(C1) Output: Private initial center set C0 ⊆ U
  where Lap(2(L−hv)/ ) is a Laplace random number with rate 2(L−hv)/ . We will use the perturbed N̂v for node sampling instead of the true value Nv , as described in Algorithm 4. The DP guarantee of this initialization scheme is straightforward by the composition of the Laplace mechanisms.
  Theorem 4.1. Algorithm 4 is -differentially private.
  Proof. For each level i, the subtrees T (v, i) are disjoint to each other. The privacy used in i-th level is /2(L−i), and the total privacy is ∑ i /2 (L−i) < .
  We now consider the approximation error. As the structure of the analysis is similar to the non-DP case, we present the main result here and defer the detailed proofs to Appendix C.
  Theorem 4.2. Algorithm 4 outputs an initial center set such that
  E[costk(D)] = O(log n)(OPTk(D) + k −14 log n).
  DP-HST Local Search. Similarly, we can use private HST initialization to improve the performance of private k-median local search, which is presented in Algorithm 5. After initialization, the DP local search procedure follows Gupta et al. (2010) using the exponential mechanism.
  Algorithm 5: DP-HST local search Input: U , demand points D ⊆ U , parameter k, , T Initialization: F1 the private initial centers generated by Algorithm 4 with privacy /2 Set parameter ′ = 24(T+1) for i = 1 to T do
  Select (x, y) ∈ Fi× (V \Fi) with prob. proportional to exp(− ′× (cost(Fi−{x}+{y})) Let Fi+1 ← Fi − {x}+ {y}
  Select j from {1, 2, ..., T + 1} with probability proportional to exp(− ′ × cost(Fj)) Output: F = Fj the private center set
  Theorem 4.3. Algorithm 5 achieves -differential privacy. The output centers admit
  costk(D) ≤ 6OPTk(D) +O( −1k24(log log n) log n)
  with probability (1− 1/poly(n)), with T = O(k log log n) iterations.
  The DP local search with random initialization (Gupta et al., 2010) has 6 multiplicative error and O( −14k2 log2 n) additive error. Our result improves the log n term to log log n in the additive error. Meanwhile, the number of iterations needed is improved from T = O(k log n) toO(k log log n) (see Appendix B for an empirical justification). Notably, it has been shown in Gupta et al. (2010) that for k-median problem, the lower bounds on the multiplicative and additive error of any -DP algorithm areO(1) andO( −14k log(n/k)), respectively. Our result matches the lower bound on the multiplicative error, and the additive error is only worse than the bound by a factor of O(k log log n) which is typically small in many cases. Thus, our private HST initialization method pushes the approximation error of private local search closer to the lower bound. To our knowledge, Theorem 4.3 is the first result in literature to improve the error of DP local search in general metric space.
sections.heading: 5 EXPERIMENTS
sections.text: We numerically test the proposed methods on two problems—clustering in an Euclidean space and on a graph. Our results show that the proposed HST initialization can improve the performance of using k-median++ initialization in both non-private and private clustering tasks.
sections.heading: 5.1 DATASETS AND ALGORITHMS
sections.text: Discrete Euclidean space. Following Balcan et al. (2017), we test k-median clustering on the MNIST hand-written digit dataset (LeCun et al., 1998) with 10 natural clusters (digit 0 to 9). We set U as 10000 randomly chosen data points. We choose the demand set D using two strategies: 1) “balance”, where we randomly choose 500 samples from U ; 2) “imbalance”, where D contains 500 random samples from U only from digit “0” and “8” (two clusters). We note that, the imbalanced D is a very practical setting in real-world scenarios, where data are typically not uniformly distributed. On this dataset, we test clustering with both l1 and l2 distance as the underlying metric.
  Metric space induced by graph. Random graphs have been widely considered in testing k-median methods (Balcan et al., 2013; Todo et al., 2019). The construction of graphs follows a similar approach as the synthetic pmedinfo graphs provided by the popular OR-Library (Beasley, 1990). The metric ρ for this experiment is the shortest (weighted) path distance on graph. To generate a size n graph, we first randomly split the nodes into 10 clusters. Within each cluster, each pair of nodes is connected with probability 0.2 and weight drawn from standard uniform distribution. For each pair of clusters, we randomly connect some nodes from each cluster, with weights following uniform [0.5, r]. A larger r makes the graph more separable, i.e., clusters are farther from each other (see Appendix B for example graphs). We present two cases: r = 1 and r = 100. For this task, U has 3000 nodes, and the private set D is chosen using similar “balanced” and “imbalanced” scheme as described above. In the imbalanced case, we choose D randomly from only two clusters.
  Algorithms. We compare the following clustering algorithms in both non-DP and DP setting.
  • NDP-rand: Local search on D with random initialization (Algorithm 1). • NDP-kmedian++: Algorithm 1 with k-median++ initialization. • NDP-HST: Algorithm 1 with HST initialization, as described in Section 3. • DP-rand: Standard private local search algorithm (Gupta et al., 2010), which is Algo-
  rithm 5 with initial centers randomly chosen from U .
  • DP-kmedian++: Algorithm 5 with k-median++ initialization run on U . • DP-HST: Private local search with HST-initialization (Algorithm 5). For non-private tasks,
  we set L = 6. For private clustering, we use L = 8.
  For non-DP methods, we set α = 10−3 in Algorithm 1 and the maximum number of iterations as 20. To examine the quality of initialization as well as the final centers, We report both the cost at initialization and the cost of the final output. For DP methods, we run the algorithms for T = 20 steps and report the results with = 1. We test k ∈ {2, 5, 10, 15, 20}. The average cost over T iterations is reported for more robustness. All results are averaged over 10 independent repetitions.
sections.heading: 5.2 RESULTS
sections.text: The results on MNIST dataset are given in Figure 2 for l1 (left two columns) and l2 (right two columns) metric. The comparisons are similar in both cases.
  • We see that the initial centers generated by HST has lower k-median cost than k-median++ and random initialization, for both non-DP and DP setting, and for both balanced and imbalanced demand set D. This confirms that the proposed HST initialization is more powerful than k-median++ in finding good initial centers.
  • From the final k-median cost plots, we also observe lower cost of HST approaches in DP clustering. In the non-DP case, the curves overlap, which means despite that HST offers better initial centers, local search can always find a good solution eventually.
  2 5 10 15 20 k
  4
  5 6 in iti al c os t
  104 Balanced D
  MNIST - l 1
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3.5
  4
  4.5
  5
  5.5
  km
  ed ia
  n co
  st
  104 Balanced D
  MNIST - l 1
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3500
  4000
  4500
  5000
  in iti
  al c
  os t
  Balanced D
  MNIST - l 2
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3500
  4000
  4500
  km
  ed ia
  n co
  st
  Balanced D
  MNIST - l 2
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  4
  5
  6
  7
  in iti
  al c
  os t
  104 Imbalanced D
  MNIST - l 1
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3.5
  4
  4.5
  5
  5.5
  6
  km
  ed ia
  n co
  st
  104 Imbalanced D
  MNIST - l 1
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3500
  4000
  4500
  5000
  5500
  in iti
  al c
  os t
  Imbalanced D
  MNIST - l 2
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  3500
  4000
  4500
  km
  ed ia
  n co
  st
  Imbalanced D
  MNIST - l 2
  NDP-HST NDP-kmedian++ NDP-rand DP-HST DP-kmedian++ DP-rand
  2 5 10 15 20 k
  2 5 10 15 20 k
  2 5 10 15 20 k
  2 5 10 15 20 k
  The set of results on the graph dataset is reported in Figure 3, which gives us similar conclusion. In all cases, our proposed HST method finds better initial centers with smaller cost than k-median++. In terms of the final clustering output, HST again considerably outperforms k-median++ in the private and imbalanced D setting, for both r = 100 (highly separable graph) and r = 1 (less separable graph). The improvement of HST over k-median++ is especially significant in the harder tasks when r = 1, i.e., the clusters are nearly mixed up.
sections.heading: 6 CONCLUSION
sections.text: In this paper, we propose a new initialization framework for the k-median problem in general discrete metric space. Our approach is called HST, which leverages tools from metric embedding theory. Our novel tree search approach has comparable efficiency and cost to k-median++ initialization. Moreover, we propose differentially private (DP) HST center initialization algorithm, which adapts to the private demand point set, leading to better clustering performance. When combined with subsequent DP local search heuristic, our algorithm is able to improve the addition error of prior works, which is close to the theoretical lower bound within a small factor. Experiments with Euclidean metrics and graph metrics verify the effectiveness of our methods.
sections.heading: A POSTPONED ALGORITHMS
sections.text: A.1 k-MEDIAN++
  In the paper, we compared our HST initialization mainly with another (perhaps most well-known) initialization algorithm for clustering, the k-median++ (Arthur & Vassilvitskii, 2007). For reference, we present the concrete procedures in Algorithm 6. Here, the function D(u,C) is the shortest distance from a data point u to the closest (center) point in set C. Arthur & Vassilvitskii (2007) showed that the output centers C by k-median++ achievesO(log k) approximation error, inO(dnk) time.
  Algorithm 6: k-median++ (Arthur & Vassilvitskii, 2007) Input: Data points U , number of centers k Set C = [ ] Randomly pick a point c1 in U and set C = C ∪ {c1} for i = 2 to k do
  Select ci = u ∈ U with probability D(u,C)∑ u′∈U D(u ′,C)
  C = C ∪ {ci} Output: k-median++ initial center set C
sections.heading: A.2 CONSTRUCTING A 2-HST
sections.text: As presented in Algorithm 7, the construction starts by applying a permutation π on U , such that in following steps the points are picked in a random sequence. We first find a padded decomposition PL = {PL,1, ..., PL,nL} of U with parameter β = 4/2. The center of each partition in PL,j serves as a root node in level L. Then, we re-do a padded decomposition for each partition PL,j , to find sub-partitions with diameter β = 4/4, and set the corresponding centers as the nodes in level L−1, and so on. Each partition at level i is obtained with β = 4/2L−i. This process proceeds until a node has a single point, or a pre-specified tree depth is reached. In Figure 1, we provide an example of L = 3-level 2-HST (left panel), along with its underlying padded decompositions (right panel).
  Algorithm 7: Build 2-HST(U,L) Input: Data points U with diameter4, L Randomly pick a point in U as the root node of T Let r = 4/2 Apply a permutation π on U // so points will be chosen in a random sequence for each v ∈ U do
  Set Cv = [v] for each u ∈ U do Add u ∈ U to Cv if d(v, u) ≤ r and u /∈ ⋃ v′ 6=v Cv′
  Set the non-empty clusters Cv as the children nodes of T for each non-empty cluster Cv do
  Run 2-HST(Cv, L− 1) to extend the tree T ; stop until L levels or reaching a leaf node Output: 2-HST T
sections.heading: B MORE EXPERIMENTS
sections.heading: B.1 EXAMPLES OF GRAPH DATA
sections.text: In Figure 4, we plot two example graphs (subgraphs of 50 nodes) with r = 100 and r = 1. When r = 100, the graph is highly separable (i.e., clusters are far from each other). When r = 1, the clusters are harder to be distinguished from each other.
  B.2 RUNNING TIME COMPARISON WITH k-MEDIAN++
  In Proposition 3.1, we show that our HST initialization algorithm admits O(dn log n) complexity when considering the Euclidean space. With a smart implementation of Algorithm 6 where each data point tracks its distance to the current closest candidate center in C, k-median++ has O(dnk) running time. Therefore, the running time of our algorithm is in general comparable to k-median++. Our method would run faster if k = Ω(log n). In Figure 5, we plot the empirical running time of HST initialization against k-median++, on MNIST dataset with l2 distance (similar comparison holds for l1). From the left subfigure, we see that k-median++ becomes slower with increasing k, and our method is more efficient when k > 20. In the right panel, we observe that the running time of both methods increases with larger sample size n. Our HST algorithm has a slightly faster increasing rate, which is predicted by the complexity comparison (n log n v.s. n). However, this difference in log n factor would not be too significant unless the sample size is extremely large. Overall, our numerical results suggest that in general, the proposed HST initialization would have similar efficiency as k-median++ in common practical scenarios.
  B.3 IMPROVED ITERATION COST OF DP-HST
  In Theorem 4.3, we show that under differential privacy constraints, the proposed DP-HST (Algorithm 5) improves both the approximation error and the number of iterations required to find a good solution of classical DP local search (Gupta et al., 2010). In this section, we provide some numerical results to justify the theory.
  First, we need to properly measure the iteration cost of DP local search. This is because, unlike the non-private clustering, the k-median cost after each iteration in DP local search is not decreasing monotonically, due to the probabilistic exponential mechanism. To this end, for the cost sequence with length T = 20, we compute its moving average sequence with window size 5. Attaining the minimal value of the moving average indicates that the algorithm has found a “local optimum”, i.e., it has reached a “neighborhood” of solutions with small clustering cost. Thus, we use the number of iterations to reach such local optimum as the measure of iteration cost. The results are provided in Figure 6. We see that on all the tasks (MNIST with l1 and l2 distance, and graph dataset with r = 1 and r = 100), DP-HST has significantly smaller iterations cost. In Figure 7, we further report the k-median cost of the best solution in T iterations found by each DP algorithm. We see that DP-HST again provide the smallest cost. This additional set of experiments again validates the claims of Theorem 4.3, that DP-HST is able to found better initial centers in fewer iterations.
sections.heading: C PROOFS
sections.text: The following composition result of differential privacy will be used in our proof.
  Theorem C.1 (Composition Theorem (Dwork, 2006)). If Algorithms A1,A2, ...,Am are 1, 2, ..., m differentially private respectively, then the union (A1(D),A2(D), ...,Am(D)) is∑m i=1 i-DP.
sections.heading: C.1 PROOF OF LEMMA 3.3
sections.text: Proof. Consider the intermediate output of Algorithm 2, C1 = {v1, v2, ..., vk}, which is the set of roots of the minimal subtrees each containing exactly one output center C0. Suppose one of the optimal “root set” that minimizes (4) is C∗1 = {v′1, v′2, ..., v′k}. If C1 = C∗1 , the proof is done. Thus, we prove the case for C1 6= C∗1 . Note that T (v), v ∈ C1 are disjoint subtrees. We have the following reasoning.
  • Case 1: for some i, j′, vi is a descendant node of v′j . Since the optimal center point f ∗ is a
  leaf node by the definition of (4), we know that there must exist one child node of v′j that expands a subtree which contains f∗. Therefore, we can always replace v′j by one of its child nodes. Hence, we can assume that vi is not a descendant of v′j .
  Note that, we have score(v′j) ≤ score(vi) if v′j /∈ C∗1 ∩C1. Algorithm 2 sorts all the nodes based on cost value, and it would have more priority to pick v′j than vi if score(v ′ j) > score(vi) and vi is not a child node of v′j .
  • Case 2: for some i, j′, v′j is a descendant of vi. In this case, optimal center point f ∗, which
  is a leaf of T (vi), must also be a leaf node of T (v′j). We can simply replace C1 with the swap C1 \ {vi} + {v′j} which does not change costTk ′(U). Hence, we can assume that v′j is not a descendant of vi.
  • Case 3: Otherwise. By the construction of C1, we know that score(v′j) ≤ min{score(vi), i = 1, ..., k} when v′j ∈ C∗1 \ C1. Consider the swap between C1 and C∗1 . By the definition of tree distance, we have OPT T k (U) ≥ ∑ vi∈C1\C∗1 Nvi2 hvi , since
  {T (vi), vi ∈ C1 \ C∗1} does not contain any center of the optimal solution determined by C∗1 (which is also the optimal “root set” for OPT T k (U)).
  Thus, we only need to consider Case 3. Let us consider the optimal clustering with center set be C∗ = {c∗1, c∗2, ..., c∗k} (each center c∗j is a leaf of subtree whose root be c′j), and S′j be the leaves assigned to c∗j . Let Sj denote the set of leaves in S ′ j whose distance to c ∗ j is strictly smaller than its distance to any centers in C1. Let Pj denote the union of paths between leaves of Sj to its closest center in C1. Let v′′j be the nodes in Pj with highest level satisfying T (v ′′ j ) ∩ C1 = ∅. The score of v′′j is 2 hv′′ j N(v′′j ). That means the swap with a center v ′ j into C1 can only reduce 4 · 2 hv′′ j N(v′′j ) to
  costTk ′(U) (the tree distance between any leaf in Sj and its closest center in C1 is at most 4 · 2 hv′′ j ). We just use v′j to represent v ′′ j for later part of this proof for simplicity. By our reasoning, summing all the swaps over C∗1 \ C1 gives
  costTk ′(U)−OPTTk (U) ≤ 4 ∑ v′j∈C∗1 \C1 Nv′j2 hv′ j ,
  OPTTk (U) ≥ ∑
  vi∈C1\C∗1
  Nvi2 hvi .
  Also, based on our discussion on Case 1, it holds that
  Nv′j2 hv′ j −Nvi2hvi ≤ 0.
  Summing them together, we have costTk ′(U) ≤ 5OPTTk (U).
sections.heading: C.2 PROOF OF LEMMA 3.4
sections.text: Proof. Since the subtrees in C1 are disjoint, it suffices to consider one subtree with root v. With a little abuse of notation, let costT1
  ′(v, U) denote the optimal k-median cost within the point set T (v) with one center in 2-HST:
  costT1 ′(v, U) = min
  x∈T (v) ∑ y∈T (v) ρT (x, y), (7)
  which is the optimal cost within the subtree. Suppose v has more than one children u,w, ..., otherwise the optimal center is clear. Suppose the optimal solution of costT1
  ′(v, U) chooses a leaf node in T (u), and our HST initialization algorithm picks a leaf of T (w). If u = w, then HST chooses the optimal one where the argument holds trivially. Thus, we consider u 6= w. We have the following two observations:
  • Since one needs to pick a leaf of T (u) to minimize costT1 ′(v, U), we have costT1 ′(v, U) ≥∑ x∈ch(v),x 6=uNx · 2hx where ch(u) denotes the children nodes of u.
  • By our greedy strategy, costT1 (v, U) ≤ ∑ x∈ch(u)Nx · 2hx ≤ costT1 ′(v, U) +Nu · 2hu .
  Since hu = hw, we have 2hu · (Nu −Nw) ≤ 0,
  since our algorithm picks subtree roots with highest scores. Then we have costT1 (v, U) ≤ costT1
  ′(v, U) + Nw · 2hw ≤ 2costT1 ′(v, U). Since the subtrees in C1 are disjoint, the union of centers for OPTT1 (v, U), v ∈ C1 forms the optimal centers with size k. Note that, for any data point p ∈ U \ C1, the tree distance ρT (p, f) for ∀f that is a leaf node of T (v), v ∈ C1 is the same. That is, the choice of leaf in T (v) as the center does not affect the k-median cost under 2-HST metric. Therefore, union bound over k subtree costs completes the proof.
sections.heading: C.3 PROOF OF PROPOSITION 3.1
sections.text: Proof. It is known that the 2-HST can be constructed in O(dn log n) (Bartal, 1996). The subtree search in Algorithm 2 involves at most sorting all the nodes in the HST based on the score, which takes O(nlogn). We use a priority queue to store the nodes in C1. When we insert a new node v into queue, its parent node (if existing in the queue) would be removed from the queue. The number of nodes is O(n) and each operation (insertion, deletion) in a priority queue based on score has O(log n) complexity. Lastly, the total time to obtain C0 is O(n), as the FIND-LEAF only requires a top down scan in k disjoint subtrees of T . Summing parts together proves the claim.
sections.heading: C.4 PROOF OF THEOREM 4.2
sections.text: Similarly, we prove the error in general metric by first analyzing the error in 2-HST metric. Then the result follows from Lemma 3.2. Let costTk (D), cost T k ′(D) and OPTTk (D) be defined analogously to (3), (4) and (5), where “y ∈ U” in the summation is changed into “y ∈ D” since D is the demand set. That is,
  costTk (D) = ∑ y∈D min x∈C0 ρT (x, y), (8)
  costTk ′(D,C1) = min
  |F∩T (v)|=1,∀v∈C1 ∑ y∈D min x∈F ρT (x, y), (9)
  OPTTk (D) = min F⊂D,|F |=k ∑ y∈D min x∈F ρT (x, y) ≡ min C′1 costTk ′(D,C ′1). (10)
  We have the following.
  Lemma C.2. costTk (D) ≤ 10OPTTk (D) + 10ck −14 log n with probability 1− 4k/nc.
  Proof. The result follows by combining the following Lemma C.4, Lemma C.5, and applying union bound.
  Lemma C.3. For any node v in T , with probability 1− 1/nc, |N̂v · 2hv −Nv · 2hv | ≤ c −14 log n.
  Proof. Since N̂v = Nv + Lap(2(L−hv)/2/ ), we have
  Pr[|N̂v −Nv| ≥ x/ ] = exp(−x/2(L−hv)).
  As L = log4, we have
  Pr[|N̂v −Nv| ≥ x4/(2hv )] ≤ exp(−x).
  Hence, for some constant c > 0,
  Pr[|N̂v · 2hv −Nv · 2hv | ≤ c −14 log n] ≥ 1− exp(−c log n) = 1− 1/nc.
  Lemma C.4 (DP Subtree Search). With probability 1 − 2k/nc, costTk ′(D) ≤ 5OPTTk (D) + 4ck −14 log n.
  Proof. The proof is similar to that of Lemma 3.3. Consider the intermediate output of Algorithm 2, C1 = {v1, v2, ..., vk}, which is the set of roots of the minimal disjoint subtrees each containing exactly one output center C0. Suppose one of the optimal “root set” that minimizes (4) is C∗1 = {v′1, v′2, ..., v′k}. Assume C1 6= C∗1 . By the same argument as the proof of Lemma 3.3, we consider for some i, j such that vi 6= v′j , where vi is not a descendent of v′j and v′j is either a descendent of vi. By the construction of C1, we know that score(v′j) ≤ min{score(vi), i = 1, ..., k} when v′j ∈ C∗1 \ C1. Consider the swap between C1 and C∗1 . By the definition of tree distance, we have OPTTk (U) ≥ ∑ vi∈C1\C∗1 Nvi2 hvi , since {T (vi), vi ∈ C1 \ C∗1} does not contain any center of the optimal solution determined byC∗1 (which is also the optimal “root set” forOPT T k ). Let us consider the optimal clustering with center set be C∗ = {c∗1, c∗2, ..., c∗k} (each center c∗j is a leaf of subtree whose root be c′j), and S ′ j be the leaves assigned to c ∗ j . Let Sj denote the set of leaves in S ′ j whose distance to c∗j is strictly smaller than its distance to any centers in C1. Let Pj denote the union of paths between leaves of Sj to its closest center in C1. Let v′′j be the nodes in Pj with highest level satisfying T (v′′j ) ∩ C1 = ∅. The score of v′′j is 2 hv′′ j N(v′′j ). That means the swap with a center v ′ j into C1 can only reduce 4 · 2 hv′′ j N(v′′j ) to cost T k ′(U) (the tree distance between any leaf in Sj and its closest center in C1 is at most 4 · 2 hv′′ j ). We just use v′j to represent v ′′ j for later part of this proof for simplicity. Summing all the swaps over C∗1 \ C1, we obtain
  costTk ′(U)−OPTTk (U) ≤ 4 ∑ v′j∈C∗1 \C1 Nv′j2 hv′ j ,
  OPTTk (U) ≥ ∑
  vi∈C1\C∗1
  Nvi2 hvi .
  Applying union bound with Lemma C.3, with probability 1− 2/nc, we have
  Nv′j2 hv′ j −Nvi2hvi ≤ 2c −14 log n.
  Consequently, we have with probability, 1− 2k/nc,
  costTk ′(D) ≤ 5OPTTk (D) + 4c|C1 \ C∗1 | −14 log n
  ≤ 5OPTTk (D) + 4ck −14 log n.
  Lemma C.5 (DP Leaf Search). With probability 1 − 2k/nc, Algorithm 4 produces initial centers with costTk (D) ≤ 2costTk ′(D) + 2ck −14 log n.
  Proof. The proof strategy follows Lemma 3.4. We first consider one subtree with root v. Let costT1 ′(v, U) denote the optimal k-median cost within the point set T (v) with one center in 2-HST:
  costT1 ′(v,D) = min
  x∈T (v) ∑ y∈T (v)∩D ρT (x, y). (11)
  Suppose v has more than one children u,w, ..., and the optimal solution of costT1 ′(v, U) chooses a leaf node in T (u), and our HST initialization algorithm picks a leaf of T (w). If u = w, then HST chooses the optimal one where the argument holds trivially. Thus, we consider u 6= w. We have the following two observations:
  • Since one needs to pick a leaf of T (u) to minimize costT1 ′(v, U), we have costT1 ′(v, U) ≥∑ x∈ch(v),x 6=uNx · 2hx where ch(u) denotes the children nodes of u.
  • By our greedy strategy, costT1 (v, U) ≤ ∑ x∈ch(u)Nx · 2hx ≤ costT1 ′(v, U) +Nu · 2hu .
  As hu = hw, leveraging Lemma C.3, with probability 1− 2/nc,
  2hu · (Nu −Nw) ≤ 2hu(N̂u − N̂w) + 2c −14 log n ≤ 2c −14 log n.
  since our algorithm picks subtree roots with highest scores. Then we have costT1 (v,D) ≤ costTk
  ′(v,D) + Nw · 2hu + 2c −14 log n ≤ 2costTk ′(v,D) + 2c −14 log n with high probability. Lastly, applying union bound over the disjoint k subtrees gives the desired result.
sections.heading: C.5 PROOF OF THEOREM 4.3
sections.text: Proof. The privacy analysis is straightforward, by using the composition theorem (Theorem C.1). Since the sensitivity of cost(·) is4, in each swap iteration the privacy budget is /2(T + 1). Also, we spend another /2(T + 1) privacy for picking a output. Hence, the total privacy is /2 for local search. Algorithm 4 takes /2 DP budget for initialization, so the total privacy is .
  The analysis of the approximation error follows from Gupta et al. (2010), where the initial cost is reduced by our private HST method. We need the following two lemmas.
  Lemma C.6 (Gupta et al. (2010)). Assume the solution to the optimal utility is unique. For any output o ∈ O of 24 -DP exponential mechanism on dataset D, it holds for ∀t > 0 that
  Pr[q(D, o) ≤ max o∈O
  q(D, o)− (ln |O|+ t)/ ] ≤ e−t,
  where |O| is the size of the output set.
  Lemma C.7 (Arya et al. (2004)). For any set F ⊆ D with |F | = k, there exists some swap (x, y) such that the local search method admits
  costk(F,D)− costk(F − {x}+ {y}, D) ≥ costk(F,D)− 5OPT (D)
  k .
  From Lemma C.7, we know that when costk(Fi, D) > 6OPT (D), there exists a swap (x, y) s.t.
  costk(Fi − {x}+ {y}, D) ≤ (1− 1
  6k )costk(Fi, D).
  At each iteration, there are at most n2 possible outputs (i.e., possible swaps), i.e., |O| = n2. Using Lemma C.6 with t = 2 log n, for ∀i,
  Pr[costk(Fi+1, D) ≥ costk(F ∗i+1, D) + 4 log n ′ ] ≥ 1− 1/n2,
  where costk(F ∗i+1, D) is the minimum cost among iteration 1, 2, ..., t + 1. Hence, we have that as long as cost(Fi, D) > 6OPT (D) + 24k logn ′ , the improvement in cost is at least by a factor of
  (1− 16k ). By Theorem 4.2, we have costk(F1, D) ≤ C(log n)(6OPT (D)+6k4 log n/ ) for some constant C > 0. Let T = 6Ck log log n. We have that
  E[cost(Fi, D)] ≤ (6OPT (D) + 6k −14 log n)C(log n)(1− 1/6k)6Ck log logn
  ≤ 6OPT (D) + 6k −14 log n ≤ 6OPT (D) + 24k log n ′ .
  Therefore, with probability at least (1 − T/n2), there exists an i <= T s.t. cost(Fi, D) ≤ 6OPT (D) + 24k logn ′ . Then by using the Lemma C.7, one will pick an Fj with additional additive error 4 lnn/ ′ to the min{cost(Fj , D), j = 1, 2, ..., T} with probability 1 − 1/n2. Consequently, we know that the expected additive error is
  24k4 log n/ ′ + 4 log n/ ′ = O( −1k24(log log n) log n),
  with probability 1− 1/poly(n).
  D EXTEND HST INITIALIZATION TO k-MEANS
  Naturally, our HST method can also be applied to k-means clustering problem. In this section, we extend the HST to k-means and provide some brief analysis similar to k-median. We present the analysis in the non-private case, which can then be easily adapted to the private case. Define the following costs for k-means.
  costTkm(U) = ∑ y∈U min x∈C0 ρT (x, y)2, (12)
  costTkm ′(U,C1) = min
  |F∩T (v)|=1,∀v∈C1 ∑ y∈U min x∈F ρT (x, y)2, (13)
  OPTTkm(U) = min F⊂U,|F |=k ∑ y∈U min x∈F ρT (x, y)2 ≡ min C′1 costTkm ′(U,C ′1). (14)
  For simplicity, we will use costTkm ′(U) to denote costTkm ′(U,C1) if everything is clear from context. Here, OPTTkm (14) is the cost of the global optimal solution with 2-HST metric. Lemma D.1 (Subtree search). costTkm′(U) ≤ 17OPTTkm(U).
  Proof. The analysis is similar with the proof of Lemma 3.3. Thus, we mainly highlight the difference. Let us just use some notations the same as in Lemma 3.3 here. Let us consider the clustering with center set be C∗ = {c∗1, c∗2, ..., c∗k} (each center c∗j is a leaf of subtree whose root be c′j), and S′j be the leaves assigned to c ∗ j in optimal k-means clustering in tree metric. Let Sj denote the set of leaves in S′j whose distance to c ∗ j is strictly smaller than its distance to any centers in C1. Let Pj denote the union of paths between leaves of Sj to its closest center in C1. Let v′′j be the nodes in Pj with highest level satisfying T (v′′j ) ∩ C1 = ∅. The score of v′′j is 2 hv′′ j N(v′′j ). That means the swap with a center v′j into C1 can only reduce (4 · 2 hv′′ j )2N(v′′j ) to cost T km ′(U). We just use v′j to represent v′′j for later part of this proof for simplicity. By our reasoning, summing all the swaps over C∗1 \ C1 gives
  costTkm ′(U)−OPTTkm(U) ≤ ∑ v′j∈C∗1 \C1 Nv′j · (4 · 2 hv′ j )2,
  OPTTkm(U) ≥ ∑
  vi∈C1\C∗1
  Nvi(2 hvi )2.
  Also, based on our discussion on Case 1, it holds that
  Nv′j2 hv′ j −Nvi2hvi ≤ 0.
  Summing them together, we have costTkm ′(U) ≤ 17OPTTkm(U).
  Next, we show that the greedy leaf search strategy (Algorithm 3) only leads to an extra multiplicative error of 2. Lemma D.2 (Leaf search). costTkm(U) ≤ 2costTkm′(U).
  Proof. Since the subtrees in C1 are disjoint, it suffices to consider one subtree with root v. With a little abuse of notation, let costT1
  ′(v, U) denote the optimal k-means cost within the point set T (v) with one center in 2-HST:
  costT1 ′(v, U) = min
  x∈T (v) ∑ y∈T (v) ρT (x, y)2, (15)
  which is the optimal cost within the subtree. Suppose v has more than one children u,w, ..., otherwise the optimal center is clear. Suppose the optimal solution of costT1
  ′(v, U) chooses a leaf node in T (u), and our HST initialization algorithm picks a leaf of T (w). If u = w, then HST chooses the optimal one where the argument holds trivially. Thus, we consider u 6= w. We have the following two observations:
  • Since one needs to pick a leaf of T (u) to minimize costT1 ′(v, U), we have costT1 ′(v, U) ≥∑ x∈ch(v),x 6=uNx · (2hx)2 where ch(u) denotes the children nodes of u.
  • By our greedy strategy, costT1 (v, U) ≤ ∑ x∈ch(u)Nx ·(2hx)2 ≤ costT1 ′(v, U)+Nu ·(2hu)2.
  Since hu = hw, we have 2hu · (Nu −Nw) ≤ 0,
  since our algorithm picks subtree roots with highest scores. Then we have costT1 (v, U) ≤ costT1
  ′(v, U) + Nw · (2hw)2 ≤ 2costT1 ′(v, U). Since the subtrees in C1 are disjoint, the union of centers for OPTT1 (v, U), v ∈ C1 forms the optimal centers with size k. Note that, for any data point p ∈ U \ C1, the tree distance ρT (p, f) for ∀f that is a leaf node of T (v), v ∈ C1 is the same. That is, the choice of leaf in T (v) as the center does not affect the k-median cost under 2-HST metric. Therefore, union bound over k subtree costs completes the proof.
  We are ready to state the error bound for our proposed HST initialization (Algorithm 2), which is a natural combination of Lemma D.1 and Lemma D.2. Theorem D.3 (HST initialization). costTkm(U) ≤ 34OPTTkm(U).
  We have the following result based on Lemma 3.2. Theorem D.4. In a general metric space,
  E[costkm(U)] = O(min{log n, log4})2OPTkm(U).