smallequation

$\BODY$

(1)

smallalign

$\BODY$

(2)

#  $k$-Median Clustering via Metric Embedding: Towards Better Initialization
with Differential Privacy

Chenglin Fan, Ping Li, Xiaoyun Li

Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004, USA

{chenglinfan2020, pingli98<EMAIL_ADDRESS>
###### Abstract

111This paper was initially made public in October 2021 at
https://openreview.net/pdf?id=beUek8ku1Q.

When designing clustering algorithms, the choice of initial centers is crucial
for the quality of the learned clusters. In this paper, we develop a new
initialization scheme, called HST initialization, 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
proposed HST initialization can produce initial centers achieving lower errors
than those from another popular initialization method, $k$-median++, with
comparable efficiency. The HST initialization can also be extended to the
setting of differential privacy (DP) to generate private initial centers. We
show that the error from applying DP local search followed by our private HST
initialization improves previous results on the approximation error, and
approaches the lower bound within a small factor. Experiments justify the
theory and demonstrate the effectiveness of our proposed method. Our approach
can also be extended to the $k$-means problem.

##  1 Introduction

Clustering is an important problem in unsupervised learning that has been
widely studied in statistics, data mining, network analysis, etc. (Punj and
Stewart, 1983; Dhillon and Modha, 2001; Banerjee et al., 2005; Berkhin, 2006;
Abbasi and 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)\leq\gamma OPT_{k}(D)+\xi,$

where $C$ is the centers output by the algorithm and $cost(C,D)$ is a cost
function defined for $C$ on dataset $D$. $OPT_{k}(D)$ is the cost of optimal
(oracle) clustering solution on $D$. When everything is clear from context, we
will use $OPT$ for short. Here, $\gamma$ is called multiplicative error and
$\xi$ 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. Particularly, $k$-median is
usually the one used for clustering on graph/network data. In general, there
are two popular frameworks for clustering. 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 (continuous) 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 ($\gamma=O(1)$) to the optimal solution for $k$-median
(Arya et al., 2004). In this paper, we will focus on clustering under this the
general metric space setting.

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 and Vassilvitskii, 2007) iteratively selects $k$ data points
as initial centers, favoring distant points in a probabilistic way.
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. Follow-up works of $k$-means++ further
improved its efficiency and scalability, e.g., Bahmani et al. (2012); Bachem
et al. (2016); Lattanzi and Sohler (2019). In this work, we propose a new
initialization framework, called Hierarchically Well-Separated Tree (HST)
initialization, based on metric embedding techniques. Our method is built upon
a novel search algorithm on metric embedding trees, with comparable
approximation error and running time as $k$-median++. Moreover, importantly,
our initialization scheme can be conveniently combined with the notion of
differential privacy (DP).

Clustering with Differential Privacy. The concept of differential privacy
(Dwork, 2006; McSherry and 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 continuous 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 yields $O(\log^{3}n)$
multiplicative error and $O((k^{2}+d)\log^{5}n)$ additive error. Later on, the
Euclidean $k$-means errors are further improved to $\gamma=O(1)$ and
$\xi=O(k^{1.01}\cdot d^{0.51}+k^{1.5})$ by Stemmer and Kaplan (2018), with
more advanced candidate set selection. Huang and 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. 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(\log^{3/2}n)$ multiplicative error and
$O((k^{2}+d)\log^{3}n)$ additive error in high dimensional Euclidean space.
However, one main limitation of these methods is that they cannot be applied
to general metric spaces (e.g., on graphs). In discrete space, Gupta et al.
(2010) proposed a private method for the classical local search heuristic,
which applies to both $k$-medians and $k$-means. To cast privacy on each
swapping step, the authors applied the exponential mechanism of McSherry and
Talwar (2007). Their method produced an $\epsilon$-differentially private
solution with cost $6OPT+O(\triangle k^{2}\log^{2}n/\epsilon)$, where
$\triangle$ is the diameter of the point set. In this work, we will show that
our HST initialization can improve DP local search for $k$-median (Gupta et
al., 2010) in terms of both approximation error and efficiency.

The main contributions of this work include :

  * •

We introduce the HST (Fakcharoenphol et al., 2004) to the $k$-median
clustering problem for initialization. We design an efficient sampling
strategy to select the initial center set from the tree, with an approximation
factor $O(\log\min\\{k,\triangle\\})$ in the non-private setting, which is
$O(\log\min\\{k,d\\})$ when $\triangle=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.

  * •

We propose a differentially private version of HST initialization under the
setting of Gupta et al. (2010) in discrete metric space. The so-called DP-HST
algorithm finds initial centers with $O(\log n)$ multiplicative error and
$O(\epsilon^{-1}\triangle k^{2}\log^{2}n)$ additive error. Moreover, running
DP-local search starting from this initialization gives $O(1)$ multiplicative
error and $O(\epsilon^{-1}\triangle k^{2}(\log\log n)\log n)$ additive error,
which improves previous results towards the well-known lower bound
$O(\epsilon^{-1}\triangle k\log(n/k))$ on the additive error of DP $k$-median
(Gupta et al., 2010) within a small $O(k\log\log n)$ factor. To our knowledge,
this is the first initialization method with differential privacy guarantee
and improved error rate in general metric space.

  * •

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 initialization approach achieves smaller initial cost than
$k$-median++ (i.e., finds better initial centers), which may also lead to
improvements in the final clustering quality.

##  2 Background and Setup

###  2.1 Differential Privacy (DP)

######  Definition 2.1 (Differential Privacy (DP) (Dwork, 2006)).

If for any two adjacent data sets $D$ and $D^{\prime}$ with symmetric
difference of size one, for any $O\subset Range(\mathbbm{A})$, an algorithm
$\mathbbm{A}$ satisfies

$Pr[\mathbbm{A}(D)\in O]\leq e^{\epsilon}Pr[\mathbbm{A}(D^{\prime})\in O],$

then algorithm $\mathbbm{A}$ is said to be $\epsilon$-differentially private.

Intuitively, differential privacy requires that after removing any data point
(graph node in our case), the output of $D^{\prime}$ should not be too
different from that of the original dataset $D$. Smaller $\epsilon$ 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, one approach is to add noise to the algorithm output. The
Laplace mechanism adds Laplace($\eta(f)/\epsilon$) noise to the output, which
is known to achieve $\epsilon$-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\times O\rightarrow\mathbb{R}$ is what we aim to maximize. The
exponential mechanism outputs an element $o\in O$ with probability
$P[\mathbbm{A}(D)=o]\propto\exp(\frac{\epsilon q(D,o)}{2\eta(q)})$, where $D$
is the input dataset and $\eta(f)=\sup_{|D-D^{\prime}|=1}|f(D)-f(D^{\prime})|$
is the sensitivity of $f$. Both mechanisms will be used in our paper.

###  2.2 Metric $k$-Median Clustering

Following Arya et al. (2004); Gupta et al. (2010), the problem of metric
$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 $\rho:U\times
U\rightarrow\mathbb{R}$, the goal of $k$-median to pick $F\subseteq U$ with
$|F|=k$ to minimize

$\displaystyle\text{{$k$-median:}}\hskip 14.45377ptcost_{k}(F,U)=\sum_{v\in
U}\min_{f\in F}\rho(v,f).$ (3)

Let $D\subseteq U$ be a set of demand points. The goal of DP $k$-median is to
minimize

$\displaystyle\text{{DP $k$-median:}}\hskip 14.45377ptcost_{k}(F,D)=\sum_{v\in
D}\min_{f\in F}\rho(v,f).$ (4)

At the same time, the output $F$ is required to be $\epsilon$-differentially
private to $D$. We may drop “$F$” and use “$cost_{k}(U)$” or “$cost_{k}(D)$”
if there is no risk of ambiguity.

To better understand the motivation of the DP clustering, we provide a 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
feature 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\in F$ and $y\in U$,
and do the swap $F\leftarrow F-\\{x\\}+\\{y\\}$ such that $F-\\{x\\}+\\{y\\}$
improves the cost of $F$ the most (if more than factor $(1-\alpha/k)$ where
$\alpha>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)\leq 5OPT$.

Input: Data points $U$, parameter $k$, constant $\alpha$

Initialization: Randomly select $k$ points from $U$ as initial center set $F$

while  _$\exists\ x\in F,y\in U$ s.t.
$cost(F-\\{x\\}+\\{y\\})\leq(1-\alpha/k)cost(F)$_ do

Select $(x,y)\in F_{i}\times(D\setminus F_{i})$ with
$\arg\min_{x,y}\\{cost(F-\\{x\\}+\\{y\\})\\}$

Swap operation: $F\leftarrow F-\\{x\\}+\\{y\\}$

Output: Center set $F$

Algorithm 1 Local search for $k$-median clustering (Arya et al., 2004)

###  2.3 $k$-median++ Initialization

Although local search is able to find a solution with constant error, it takes
$O(n^{2})$ per iteration Resende and Werneck (2007) in expected $O(k\log n)$
steps (in total $O(kn^{2}\log n)$) when started from random center set, which
would be slow for large datasets. Indeed, we do not need such complicated
algorithm to reduce the cost at the beginning, i.e., when the cost is large.
To accelerate the process, efficient initialization methods find a “roughly”
good center set as the starting point for local search. In the paper, we
compare our new initialization scheme mainly with a popular (and perhaps most
well-known) initialization method, the $k$-median++ (Arthur and Vassilvitskii,
2007) (see Algorithm 2).

Input: Data points $U$, number of centers $k$

Randomly pick a point $c_{1}\in U$ and set $F=\\{c_{1}\\}$

for  _$i=2$ to $k$_ do

Select $c_{i}=u\in U$ with probability $\frac{\rho(u,F)}{\sum_{u^{\prime}\in
U}\rho(u^{\prime},F)}$

$F=F\cup\\{c_{i}\\}$

Output: $k$-median++ initial center set $F$

Algorithm 2 $k$-median++ initialization (Arthur and Vassilvitskii, 2007)

Here, the function $D(u,C)$ is the shortest distance from a data point $u$ to
the closest (center) point in set $C$. Arthur and Vassilvitskii (2007) showed
that the output centers $C$ by $k$-median++ achieves $O(\log k)$ approximation
error with time complexity $O(nk)$. Starting from the initialization, we only
need to run $O(k\log\log k)$ steps of the computationally heavy local search
to reach a constant error solution. Thus, initialization may greatly improve
the clustering efficiency.

##  3 Initialization via Hierarchically Well-Separated Tree (HST)

In this section, we propose our novel initialization scheme for $k$-median
clustering, and provide our analysis in the non-private case solving (3). The
idea is based on the metric embedding theory. We will start with an
introduction to the main tool used in our approach.

###  3.1 Hierarchically Well-Separated Tree (HST)

In this paper, for an $L$-level tree, we will count levels in descending order
down the tree. We use $h_{v}$ to denote the level of $v$, and $n_{i}$ 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,\rho)$ be a metric space with
$|U|=n$, and we will refer to this metric space without specific
clarification. A $\beta$–padded decomposition of $U$ is a probabilistic
distribution of partitions of $U$ such that the diameter of each cluster
$U_{i}\in U$ is at most $\beta$, i.e., $\rho(u,v)\leq\beta$, $\forall u,v\in
U_{i}$, $i=1,...,k$. The formal definition of HST is given as below.

######  Definition 3.1.

Assume $\min_{u,v\in U}\rho(u,v)=1$ and denote $\triangle=\max_{u,v\in
U}\rho(u,v)$. An $\alpha$-Hierarchically Well-Separated Tree ($\alpha$-HST)
with depth $L$ 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 most
$\triangle/\alpha^{L-i}$.

In this paper, we consider $\alpha=2$-HST for simplicity, as $\alpha$ only
affects the constants in our theoretical analysis. As presented in Algorithm
3, the construction starts by applying a permutation $\pi$ on $U$, such that
in following steps the points are picked in a random sequence. We first find a
padded decomposition $P_{L}=\\{P_{L,1},...,P_{L,n_{L}}\\}$ of $U$ with
parameter $\beta=\triangle/2$. The center of each partition in $P_{L,j}$
serves as a root node in level $L$. Then, we re-do a padded decomposition for
each partition $P_{L,j}$, to find sub-partitions with diameter
$\beta=\triangle/4$, and set the corresponding centers as the nodes in level
$L-1$, and so on. Each partition at level $i$ is obtained with
$\beta=\triangle/2^{L-i}$. This process proceeds until a node has a single
point, or a pre-specified tree depth is reached.

Figure 1: An illustrative example of a 3-level padded decomposition and its
corresponding 2-HST. Left: The thickness of the ball represents the level. The
color corresponds to the levels in the HST in the right panel. “$\triangle$”’s
are the center nodes of partitions (balls), and “$\times$”’s are non-center
data points. Right: The resulting 2-HST generated from the padded
decomposition.

In Figure 1, we provide an example of $L=3$-level 2-HST (left panel), along
with its underlying padded decompositions (right panel). Besides this basic
implementation for better illustration, Blelloch et al. (2017) proposed an
efficient HST construction in $O(m\log n)$ time, where $n$ and $m$ are the
number of nodes and the number of edges in a graph, respectively.

Input: Data points $U$ with diameter $\triangle$, $L$

Randomly pick a point in $U$ as the root node of $T$

Let $r=\triangle/2$

Apply a permutation $\pi$ on $U$

// so points will be chosen in a random sequence

for  _each $v\in U$_ do

Set $C_{v}=[v]$

for  _each $u\in U$_ do

Add $u\in U$ to $C_{v}$ if $d(v,u)\leq r$ and $u\notin\bigcup_{v^{\prime}\neq
v}C_{v^{\prime}}$

Set the non-empty clusters $C_{v}$ as the children nodes of $T$

for  _each non-empty cluster $C_{v}$_ do

Run 2-HST$(C_{v},L-1)$ to extend the tree $T$; stop until $L$ levels or
reaching a leaf node

Output: 2-HST $T$

Algorithm 3 Build 2-HST($U,L$)

Input: $U$, $\triangle$, $k$

Initialization: $L=\log\triangle$, $C_{0}=\emptyset,C_{1}=\emptyset$

Call Algorithm 3 to build a level-$L$ 2-HST $T$ using $U$

for  _each node $v$ in $T$_ do

$N_{v}\leftarrow|U\cap T(v)|$

$score(v)\leftarrow N_{v}\cdot 2^{h_{v}}$

while  _$|C_{1}| <k$ _ do

Add top $(k-|C_{1}|)$ nodes with highest score to $C_{1}$

for  _each $v\in C_{1}$_ do

$C_{1}=C_{1}\setminus\\{v\\}$, if $\exists\ v^{\prime}\in C_{1}$ such that
$v^{\prime}$ is a descendant of $v$

$C_{0}=\textrm{FIND-LEAF}(T,C_{1})$

Output: Initial center set $C_{0}\subseteq U$

Algorithm 4 NDP-HST initialization

The first step of our method is to embed the data points into an HST (see
Algorithm 4). Next, we will describe our proposed new strategy to search for
the initial centers on the tree (w.r.t. the tree metric). Before moving on, 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 (e.g.,
$O(kn^{2})$), making them unsuitable for the purpose of fast initialization.
Moreover, it is unknown how to apply them effectively to the private case. As
will be shown, our new algorithm 1) is very efficient, 2) gives $O(1)$
approximation error in the tree metric, and 3) can be effectively extended to
DP easily (Section 4).

###  3.2 HST Initialization Algorithm

Let $L=\log\Delta$ and suppose $T$ is a level-$L$ 2-HST in $(U,\rho)$, where
we assume $L$ is an integer. For a node $v$ at level $i$, we use $T(v)$ to
denote the subtree rooted at $v$. Let $N_{v}=|T(v)|$ be the number of data
points in $T(v)$. The search strategy for the initial centers, NDP-HST
initialization (“NDP” stands for “Non-Differentially Private”), is presented
in Algorithm 4 with two phases.

Subtree search. The first step is to identify the subtrees that contain the
$k$ centers. To begin with, $k$ initial centers $C_{1}$ are picked from $T$
who have the largest $score(v)=N(v)\cdot 2^{h_{v}}$. 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 $C_{1}$: if there is
any ancestor-descendant pair of nodes, we remove the ancestor from $C_{1}$. If
the current size of $C_{1}$ is smaller than $k$, we repeat the process until
$k$ centers are chosen (we do not re-select nodes in $C_{1}$ and their
ancestors). This way, $C_{1}$ contains $k$ root nodes of $k$ disjoint
subtrees.

Input: $T$, $C_{1}$

Initialization: $C_{0}=\emptyset$

for  _each node $v$ in $C_{1}$_ do

while  _$v$ is not a leaf node_ do

$v\leftarrow\arg_{w}\max\\{N_{w},w\in ch(v)\\}$, where $ch(v)$ denotes the
children nodes of $v$

Add $v$ to $C_{0}$

Output: Initial center set $C_{0}\subseteq U$

Algorithm 5 FIND-LEAF ($T,C_{1}$)

Leaf search. After we find $C_{1}$ the set of $k$ subtrees, the next step is
to find the center in each subtree using Algorithm 5 (“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 complexity of our search algorithm is given as follows.

######  Proposition 3.1 (Complexity).

Algorithm 4 takes $O(dn\log n)$ time in the Euclidean space.

######  Remark 3.1.

The complexity of HST initialization is in general comparable to $O(dnk)$ of
$k$-median++. Our algorithm would be faster if $k>\log n$, i.e., the number of
centers is large. Similar comparison also holds for general metrics.

###  3.3 Approximation Error of HST Initialization

Firstly, we show that the initial center set produced by NDP-HST is already a
good approximation to the optimal $k$-median solution. Let
$\rho^{T}(x,y)=d_{T}(x,y)$ denote the “2-HST metric” between $x$ and $y$ in
the 2-HST $T$, where $d_{T}(x,y)$ is the tree distance between nodes $x$ and
$y$ in $T$. By Definition 3.1 and since $\triangle=2^{L}$, in the analysis we
assume equivalently that the edge weight of the $i$-th level $2^{i-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 adapted to the general metrics
by the following Lemma (Bartal, 1996).

######  Lemma 3.2.

In a metric space $(U,\rho)$ with $|U|=n$ and diameter $\triangle$, it holds
that $E[\rho^{T}(x,y)]=O(\min\\{\log n,\log\triangle\\})\rho(x,y)$. In the
Euclidean space $\mathbb{R}^{d}$, $E[\rho^{T}(x,y)]=O(d)\rho(x,y)$.

Recall $C_{0},C_{1}$ from Algorithm 4. We define

$\displaystyle cost_{k}^{T}(U)$ $\displaystyle=\sum_{y\in U}\min_{x\in
C_{0}}\rho^{T}(x,y),$ (5) $\displaystyle{cost_{k}^{T}}^{\prime}(U,C_{1})$
$\displaystyle=\min_{\begin{subarray}{c}|F\cap T(v)|=1,\\\ \forall v\in
C_{1}\end{subarray}}\sum_{y\in U}\min_{x\in F}\rho^{T}(x,y),$ (6)
$\displaystyle OPT_{k}^{T}(U)$ $\displaystyle=\min_{F\subset
U,|F|=k}\sum_{y\in U}\min_{x\in F}\rho^{T}(x,y)\equiv\min_{C_{1}^{\prime}}\
{cost_{k}^{T}}^{\prime}(U,C_{1}^{\prime}).$ (7)

For simplicity, we will use ${cost_{k}^{T}}^{\prime}(U)$ to denote
${cost_{k}^{T}}^{\prime}(U,C_{1})$. Here, $OPT_{k}^{T}$ (7) is the cost of the
global optimal solution with 2-HST metric. The last equivalence in (7) holds
because the optimal centers set can always located in $k$ disjoint subtrees,
as each leaf only contain one point. (5) is the $k$-median cost with 2-HST
metric of the output $C_{0}$ of Algorithm 4. (6) is the oracle cost after the
subtrees are chosen. That is, it represents the optimal cost to pick one
center from each subtree in $C_{1}$. Firstly, we bound the approximation error
of subtree search and leaf search, respectively.

######  Lemma 3.3 (Subtree search).

${cost_{k}^{T}}^{\prime}(U)\leq 5OPT^{T}_{k}(U)$.

######  Lemma 3.4 (Leaf search).

$cost^{T}_{k}(U)\leq 2{cost_{k}^{T}}^{\prime}(U)$.

Combining Lemma 3.3 and Lemma 3.4, we obtain

######  Theorem 3.5 (2-HST error).

Running Algorithm 4, we have $cost^{T}_{k}(U)\leq 10OPT^{T}_{k}(U)$.

Thus, HST-initialization produces an $O(1)$ approximation to $OPT$ in the
2-HST metric. Define $cost_{k}(U)$ as (3) for our HST centers, and the optimal
cost w.r.t. $\rho$ as

$\displaystyle OPT_{k}(U)=\min_{|F|=k}\sum_{y\in U}\min_{x\in F}\rho(x,y).$
(8)

We have the following result based on Lemma 3.2.

######  Theorem 3.6.

In general metric space, the expected $k$-median cost of Algorithm 4 is
$E[cost_{k}(U)]=O(\min\\{\log n,\log\triangle\\})OPT_{k}(U)$.

######  Remark 3.2.

In the Euclidean space, Makarychev et al. (2019) proved $O(\log k)$ random
projections suffice for $k$-median to achieve $O(1)$ error. Thus, if
$\triangle=O(d)$ (e.g., bounded data), by Lemma 3.2, HST initialization is
able to achieve $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 are interested in the approximation quality of
standard local search (Algorithm 1), when initialized by our NDP-HST.

######  Theorem 3.7.

NDP-HST local search achieves $O(1)$ approximation error in expected
$O(k\log\log\min\\{n,\triangle\\})$ number of iterations for input in general
metric space.

Before ending this section, we remark that the initial centers found by NDP-
HST can be used for $k$-means clustering analogously. For general metrics,
$E[cost_{km}(U)]=O(\min\\{\log n,\log\triangle\\})^{2}OPT_{km}(U)$ where
$cost_{km}(U)$ is the optimal $k$-means cost. See Appendix B for the detailed
(and similar) analysis.

##  4 HST Initialization with Differential Privacy

In this section, we consider initialization method with differential privacy
(DP). Recall (4) that $U$ is the universe of data points, and $D\subset U$ is
a demand set that needs to be clustered with privacy. Since $U$ is public,
simply running initialization algorithms on $U$ would preserve the privacy of
$D$. However, 1) this might be too expensive; 2) 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 proposed 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). Again, suppose $T$ is
an $L=\log\triangle$-level 2-HST of universe $U$ in a general metric space.
Denote $N_{v}=|T(v)\cap D|$ for a node point $v$. Our private HST
initialization (DP-HST) is similar to the non-private Algorithm 4. To gain
privacy, we perturb $N_{v}$ by adding i.i.d. Laplace noise:

$\hat{N_{v}}=N_{v}+Lap(2^{(L-h_{v})}/\epsilon),$

where $Lap(2^{(L-h_{v})}/\epsilon)$ is a Laplace random number with rate
$2^{(L-h_{v})}/\epsilon$. We will use the perturbed $\hat{N}_{v}$ for node
sampling instead of the true value $N_{v}$, as described in Algorithm 6. The
DP guarantee of this initialization scheme is straightforward by the
composition theory (Dwork, 2006).

Input: $U,D$, $\triangle$, $k$, $\epsilon$

Build a level-$L$ 2-HST $T$ based on input $U$

for  _each node $v$ in $T$_ do

$N_{v}\leftarrow|D\cap T(v)|$

$\hat{N_{v}}\leftarrow N_{v}+Lap(2^{(L-h_{v})}/\epsilon)$

$score(v)\leftarrow\hat{N}(v)\cdot 2^{h_{v}}$

Based on $\hat{N}_{v}$, apply the same strategy as Algorithm 4: find $C_{1}$;
$C_{0}=\text{FIND-LEAF(}C_{1}\text{)}$

Output: Private initial center set $C_{0}\subseteq U$

Algorithm 6 DP-HST initialization

######  Theorem 4.1.

Algorithm 6 is $\epsilon$-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 $\epsilon/2^{(L-i)}$, and the total privacy is
$\sum_{i}\epsilon/2^{(L-i)}<\epsilon$. ∎

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 A.

######  Theorem 4.2.

Algorithm 6 finds initial centers such that

$E[cost_{k}(D)]=O(\log n)(OPT_{k}(D)+k\epsilon^{-1}\triangle\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 7. After initialization, the DP local search procedure follows
Gupta et al. (2010) using the exponential mechanism.

Input: $U$, demand points $D\subseteq U$, parameter $k,\epsilon$, $T$

Initialization: $F_{1}$ the private initial centers generated by Algorithm 6
with privacy $\epsilon/2$

Set parameter $\epsilon^{\prime}=\frac{\epsilon}{4\triangle(T+1)}$

for  _$i=1$ to $T$_ do

Select $(x,y)\in F_{i}\times(V\setminus F_{i})$ with prob. proportional to
$\exp(-\epsilon^{\prime}\times(cost(F_{i}-\\{x\\}+\\{y\\}))$

Let $F_{i+1}\leftarrow F_{i}-\\{x\\}+\\{y\\}$

Select $j$ from $\\{1,2,...,T+1\\}$ with probability proportional to
$\exp(-\epsilon^{\prime}\times~cost(F_{j}))$

Output: $F=F_{j}$ the private center set

Algorithm 7 DP-HST local search

######  Theorem 4.3.

Algorithm 7 achieves $\epsilon$-differential privacy. With probability
$(1-\frac{1}{poly(n)})$, the output centers admit

$\displaystyle cost_{k}(D)\leq
6OPT_{k}(D)+O(\epsilon^{-1}k^{2}\triangle(\log\log n)\log n)$

in $T=O(k\log\log n)$ iterations.

The DP local search with random initialization (Gupta et al., 2010) has 6
multiplicative error and $O(\epsilon^{-1}\triangle k^{2}\log^{2}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)$ to $O(k\log\log n)$ (see Section 5.3 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 $\epsilon$-DP
algorithm are $O(1)$ and $O(\epsilon^{-1}\triangle k\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 would be small in many cases. To our knowledge, Theorem 4.3 is the first
result in literature to improve the error of DP local search in general metric
space.

##  5 Experiments

###  5.1 Datasets and Algorithms

Figure 2: Example of synthetic graphs: subgraph of 50 nodes. Left: $r=1$.
Right: $r=100$. Darker and thicker edged have smaller distance. When $r=100$,
the graph is more separable.

Discrete Euclidean space. Following previous work ., 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 $l_{1}$ and $l_{2}$ 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 $\rho$
for this experiment is the shortest (weighted) path distance. 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. In Figure 2, we plot two example graphs (subgraphs of
50 nodes) with $r=100$ and $r=1$. We present two cases: $r=1$ and $r=100$. For
this task, $U$ has 3000 nodes, and the private set $D$ (500 nodes) 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: (1) NDP-rand: Local search with random initialization; (2) NDP-
kmedian++: Local search with $k$-median++ initialization (Algorithm 2); (3)
NDP-HST: Local search with NDP-HST initialization (Algorithm 4), as described
in Section 3; (4) DP-rand: Standard DP local search algorithm (Gupta et al.,
2010), which is Algorithm 7 with initial centers randomly chosen from $U$; (5)
DP-kmedian++: DP local search with $k$-median++ initialization run on $U$; (6)
DP-HST: DP local search with HST-initialization (Algorithm 7). For non-DP
tasks, we set $L=6$. For DP clustering, we use $L=8$.

For non-DP methods, we set $\alpha=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 $\epsilon=1$. We test $k\in\\{2,5,10,15,20\\}$.
The average cost over $T$ iterations is reported for more robustness. All
results are averaged over 10 independent repetitions.

###  5.2 Results

The results on MNIST dataset are given in Figure 3. The comparisons are
similar for both $l_{1}$ and $l_{2}$:

  * •

From the left column, the initial centers found by HST has lower 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 right column, we also observe lower final cost of HST followed by
local search in DP clustering. In the non-DP case, the final cost curves
overlap, which means that despite HST offers better initial centers, local
search can always find a good solution eventually.

  * •

The advantage of DP-HST, in terms of both the initial and the final cost, is
more significant when $D$ is an imbalanced subset of $U$. As mentioned before,
this is because our DP-HST initialization approach also privately incorporates
the information of $D$.

The results on graphs are reported in Figure 4, which give similar
conclusions. In all cases, our proposed HST scheme finds better initial
centers with smaller cost than $k$-median++. Moreover, HST again considerably
outperforms $k$-median++ in the private and imbalanced $D$ setting, for both
$r=100$ (highly separable) and $r=1$ (less separable). The advantages of HST
over $k$-median++ are especially significant in the harder tasks when $r=1$,
i.e., the clusters are nearly mixed up.

Figure 3: Initial and final $k$-median cost on MNIST dataset. 1st column:
$l_{1}$ distance. 2nd column: $l_{2}$ distance.

Figure 4: Initial and final $k$-median cost on graph dataset. 1st column:
$l_{1}$ distance. 2nd column: $l_{2}$ distance.

###  5.3 Improved Iteration Cost of DP-HST

In Theorem 4.3, we show that under differential privacy constraints, the
proposed DP-HST (Algorithm 7) 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 5. We see that on all the tasks (MNIST with $l_{1}$ and
$l_{2}$ distance, and graph dataset with $r=1$ and $r=100$), DP-HST has
significantly smaller iterations cost. In Figure 6, 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.

Figure 5: Iteration cost to reach a locally optimal solution, on MNIST and
graph datasets with different $k$. The demand set is an imbalanced subset of
the universe.

Figure 6: The $k$-median cost of the best solution found by each
differentially private algorithm. The demand set is an imbalanced subset of
the universe. Same comparison holds on graph data.

###  5.4 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 2 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=\Omega(\log n)$.

Figure 7: Empirical time comparison of HST initialization v.s. $k$-median++,
on MNIST dataset with $l_{2}$ distance. Left: The running time against $k$, on
a subset of $n=2000$ data points. Right: The running time against $n$, with
$k=20$ centers.

In Figure 7, we plot the empirical running time of HST initialization against
$k$-median++, on MNIST dataset with $l_{2}$ distance (similar comparison holds
for $l_{1}$). 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 results suggest that
in general, the proposed HST initialization would have similar efficiency as
$k$-median++ in common practical scenarios.

##  6 Conclusion

In this paper, we propose a new initialization framework for the metric
$k$-median problem in general (discrete) metric space. Our approach is called
HST initialization, which leverages tools from metric embedding theory. Our
novel tree search approach has comparable efficiency and approximation error
to the popular $k$-median++ initialization. Moreover, we propose the
differentially private (DP) HST 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 additive error of DP local search and our result is close to the
theoretical lower bound within a small factor. Experiments with Euclidean
metrics and graph metrics verify the effectiveness of our method, which
improves the cost of both the initial centers and the final $k$-median output.

## References

  * Abbasi and Younis (2007) Ameer Ahmed Abbasi and Mohamed F. Younis.  A survey on clustering algorithms for wireless sensor networks.  _Comput. Commun._ , 30(14-15):2826–2841, 2007\.
  * Arthur and Vassilvitskii (2007) David Arthur and Sergei Vassilvitskii.  k-means++: the advantages of careful seeding.  In _Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , pages 1027–1035, New Orleans, LA, 2007.
  * Arya et al. (2004) Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and Vinayaka Pandit.  Local search heuristics for k-median and facility location problems.  _SIAM J. Comput._ , 33(3):544–562, 2004.
  * Bachem et al. (2016) Olivier Bachem, Mario Lucic, S. Hamed Hassani, and Andreas Krause.  Approximate k-means++ in sublinear time.  In Dale Schuurmans and Michael P. Wellman, editors, _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 1459–1467, Phoenix, AZ, 2016.
  * Bahmani et al. (2012) Bahman Bahmani, Benjamin Moseley, Andrea Vattani, Ravi Kumar, and Sergei Vassilvitskii.  Scalable k-means++.  _Proc. VLDB Endow._ , 5(7):622–633, 2012.
  * Balcan et al. (2013) Maria-Florina Balcan, Steven Ehrlich, and Yingyu Liang.  Distributed k-means and k-median clustering on general communication topologies.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 1995–2003, Lake Tahoe, NV, 2013.
  * Balcan et al. (2017) Maria-Florina Balcan, Travis Dick, Yingyu Liang, Wenlong Mou, and Hongyang Zhang.  Differentially private clustering in high-dimensional euclidean spaces.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pages 322–331, Sydney, Australia, 2017.
  * Banerjee et al. (2005) Arindam Banerjee, Srujana Merugu, Inderjit S. Dhillon, and Joydeep Ghosh.  Clustering with bregman divergences.  _J. Mach. Learn. Res._ , 6:1705–1749, 2005.
  * Bartal (1996) Yair Bartal.  Probabilistic approximations of metric spaces and its algorithmic applications.  In _Proceedings of the 37th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 184–193, Burlington, VT, 1996.
  * Beasley (1990) John E Beasley.  OR-Library: distributing test problems by electronic mail.  _Journal of the Operational Research Society_ , 41(11):1069–1072, 1990.
  * Berkhin (2006) Pavel Berkhin.  A survey of clustering data mining techniques.  In _Grouping Multidimensional Data_ , pages 25–71. Springer, 2006\.
  * Blelloch et al. (2017) Guy E. Blelloch, Yan Gu, and Yihan Sun.  Efficient construction of probabilistic tree embeddings.  In _Proceedings of the 44th International Colloquium on Automata, Languages, and Programming (ICALP)_ , pages 26:1–26:14, Warsaw, Poland, 2017.
  * Dhillon and Modha (2001) Inderjit S. Dhillon and Dharmendra S. Modha.  Concept decompositions for large sparse text data using clustering.  _Mach. Learn._ , 42(1/2):143–175, 2001.
  * Dwork (2006) Cynthia Dwork.  Differential privacy.  In _Proceedings of the 33rd International Colloquium on Automata, Languages and Programming (ICALP),Part II_ , pages 1–12, Venice, Italy, 2006.
  * Fakcharoenphol et al. (2004) Jittat Fakcharoenphol, Satish Rao, and Kunal Talwar.  A tight bound on approximating arbitrary metrics by tree metrics.  _J. Comput. Syst. Sci._ , 69(3):485–497, 2004\.
  * Feldman et al. (2009) Dan Feldman, Amos Fiat, Haim Kaplan, and Kobbi Nissim.  Private coresets.  In _Proceedings of the 41st Annual ACM Symposium on Theory of Computing (STOC)_ , pages 361–370, Bethesda, MD, 2009.
  * Feldman et al. (2017) Dan Feldman, Chongyuan Xiang, Ruihao Zhu, and Daniela Rus.  Coresets for differentially private k-means clustering and applications to privacy in mobile sensor networks.  In _Proceedings of the 16th ACM/IEEE International Conference on Information Processing in Sensor Networks (IPSN)_ , pages 3–15, Pittsburgh, PA, 2017.
  * Gupta et al. (2010) Anupam Gupta, Katrina Ligett, Frank McSherry, Aaron Roth, and Kunal Talwar.  Differentially private combinatorial optimization.  In _Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , pages 1106–1125, Austin, TX, 2010.
  * Huang and Liu (2018) Zhiyi Huang and Jinyan Liu.  Optimal differentially private algorithms for k-means clustering.  In _Proceedings of the 37th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems (PODS)_ , pages 395–408, Houston, TX, 2018.
  * Kanungo et al. (2002) Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman, and Angela Y. Wu.  A local search approximation algorithm for k-means clustering.  In _Proceedings of the 18th Annual Symposium on Computational Geometry (CG)_ , pages 10–18, Barcelona, Spain, 2002.
  * Kaufman et al. (1977) Leon Kaufman, Marc Vanden Eede, and Pierre Hansen.  A plant and warehouse location problem.  _Journal of the Operational Research Society_ , 28(3):547–554, 1977.
  * Lattanzi and Sohler (2019) Silvio Lattanzi and Christian Sohler.  A better k-means++ algorithm via local search.  In _Proceedings of the 36th International Conference on Machine Learning (ICML)_ , pages 3662–3671, Long Beach, CA, 2019.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Lloyd (1982) Stuart P. Lloyd.  Least squares quantization in PCM.  _IEEE Trans. Inf. Theory_ , 28(2):129–136, 1982\.
  * Makarychev et al. (2019) Konstantin Makarychev, Yury Makarychev, and Ilya P. Razenshteyn.  Performance of johnson-lindenstrauss transform for _k_ -means and _k_ -medians clustering.  In _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC)_ , pages 1027–1038, Phoenix, AZ, 2019.
  * McSherry and Talwar (2007) Frank McSherry and Kunal Talwar.  Mechanism design via differential privacy.  In _Proceedings of the 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS)_ , pages 94–103, Providence, RI, 2007.
  * Nock et al. (2016) Richard Nock, Raphaël Canyasse, Roksana Boreli, and Frank Nielsen.  k-variates++: more pluses in the k-means++.  In _Proceedings of the 33nd International Conference on Machine Learning (ICML)_ , pages 145–154, New York City, NY, 2016.
  * Punj and Stewart (1983) Girish Punj and David W Stewart.  Cluster analysis in marketing research: Review and suggestions for application.  _Journal of Marketing Research_ , 20(2):134–148, 1983.
  * Resende and Werneck (2007) Mauricio G. C. Resende and Renato Fonseca F. Werneck.  A fast swap-based local search procedure for location problems.  _Ann. Oper. Res._ , 150(1):205–230, 2007.
  * Shah (2003) Rahul Shah.  Faster algorithms for k-median problem on trees with smaller heights.  _Technical report_ , 2003.
  * Stemmer and Kaplan (2018) Uri Stemmer and Haim Kaplan.  Differentially private k-means with constant multiplicative error.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 5436–5446, Montréal, Canada, 2018.
  * Tamir (1996) Arie Tamir.  An $o(pn^{2})$ algorithm for the p-median and related problems on tree graphs.  _Oper. Res. Lett._ , 19(2):59–64, 1996.
  * Todo et al. (2019) Keisuke Todo, Atsuyoshi Nakamura, and Mineichi Kudo.  A fast approximate algorithm for k-median problem on a graph.  In _Proceedings of the 15th International Workshop on Mining and Learning with Graphs (MLG)_ , Anchorage, AK, 2019.

##  Appendix A Proofs

The following composition result of differential privacy will be used in our
proof.

######  Theorem A.1 (Composition Theorem [Dwork, 2006]).

If Algorithms $\mathbbm{A}_{1},\mathbbm{A}_{2},...,\mathbbm{A}_{m}$ are
$\epsilon_{1},\epsilon_{2},...,\epsilon_{m}$ differentially private
respectively, then the union
$(\mathbbm{A}_{1}(D),\mathbbm{A}_{2}(D),...,\mathbbm{A}_{m}(D))$ is
$\sum_{i=1}^{m}\epsilon_{i}$-DP.

###  A.1 Proof of Lemma 3.3

###### Proof.

Consider the intermediate output of Algorithm 4,
$C_{1}=\\{v_{1},v_{2},...,v_{k}\\}$, which is the set of roots of the minimal
subtrees each containing exactly one output center $C_{0}$. Suppose one of the
optimal “root set” that minimizes (6) is
$C^{*}_{1}=\\{v^{\prime}_{1},v^{\prime}_{2},...,v^{\prime}_{k}\\}$. If
$C_{1}=C^{*}_{1}$, the proof is done. Thus, we prove the case for $C_{1}\neq
C_{1}^{*}$. Note that $T(v),v\in C_{1}$ are disjoint subtrees. We have the
following reasoning.

  * •

Case 1: for some $i,j^{\prime}$, $v_{i}$ is a descendant node of
$v_{j}^{\prime}$. Since the optimal center point $f^{*}$ is a leaf node by the
definition of (6), we know that there must exist one child node of
$v_{j}^{\prime}$ that expands a subtree which contains $f^{*}$. Therefore, we
can always replace $v_{j}^{\prime}$ by one of its child nodes. Hence, we can
assume that $v_{i}$ is not a descendant of $v_{j}^{\prime}$.

Note that, we have $score(v^{\prime}_{j})\leq score(v_{i})$ if
$v^{\prime}_{j}\notin C^{*}_{1}\cap C_{1}$. Algorithm 4 sorts all the nodes
based on cost value, and it would have more priority to pick $v^{\prime}_{j}$
than $v_{i}$ if $score(v^{\prime}_{j})>score(v_{i})$ and $v_{i}$ is not a
child node of $v^{\prime}_{j}$.

  * •

Case 2: for some $i,j^{\prime}$, $v_{j}^{\prime}$ is a descendant of $v_{i}$.
In this case, optimal center point $f^{*}$, which is a leaf of $T(v_{i})$,
must also be a leaf node of $T(v^{\prime}_{j})$. We can simply replace $C_{1}$
with the swap $C_{1}\setminus\\{v_{i}\\}+\\{v_{j}^{\prime}\\}$ which does not
change ${cost_{k}^{T}}^{\prime}(U)$. Hence, we can assume that
$v_{j}^{\prime}$ is not a descendant of $v_{i}$.

  * •

Case 3: Otherwise. By the construction of $C_{1}$, we know that
$score(v^{\prime}_{j})\leq\min\\{score(v_{i}),i=1,...,k\\}$ when
$v^{\prime}_{j}\in C^{*}_{1}\setminus C_{1}$. Consider the swap between
$C_{1}$ and $C^{*}_{1}$. By the definition of tree distance, we have
$OPT^{T}_{k}(U)\geq\sum_{v_{i}\in C_{1}\setminus
C^{*}_{1}}N_{v_{i}}2^{h_{v_{i}}}$, since $\\{T(v_{i}),v_{i}\in C_{1}\setminus
C^{*}_{1}\\}$ does not contain any center of the optimal solution determined
by $C^{*}_{1}$ (which is also the optimal “root set” for $OPT_{k}^{T}(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^{\prime}_{j}$), and
$S^{\prime}_{j}$ be the leaves assigned to $c^{*}_{j}$. Let $S_{j}$ denote the
set of leaves in $S^{\prime}_{j}$ whose distance to $c^{*}_{j}$ is strictly
smaller than its distance to any centers in $C_{1}$. Let $P_{j}$ denote the
union of paths between leaves of $S_{j}$ to its closest center in $C_{1}$. Let
$v^{\prime\prime}_{j}$ be the nodes in $P_{j}$ with highest level satisfying
$T(v^{\prime\prime}_{j})\cap C_{1}=\emptyset$. The score of
$v^{\prime\prime}_{j}$ is
$2^{h_{v^{\prime\prime}_{j}}}N(v^{\prime\prime}_{j})$. That means the swap
with a center $v^{\prime}_{j}$ into $C_{1}$ can only reduce $4\cdot
2^{h_{v^{\prime\prime}_{j}}}N(v^{\prime\prime}_{j})$ to
${cost_{k}^{T}}^{\prime}(U)$ (the tree distance between any leaf in $S_{j}$
and its closest center in $C_{1}$ is at most $4\cdot
2^{h_{v^{\prime\prime}_{j}}}$). We just use $v^{\prime}_{j}$ to represent
$v^{\prime\prime}_{j}$ for later part of this proof for simplicity. By our
reasoning, summing all the swaps over $C^{*}_{1}\setminus C_{1}$ gives

${cost_{k}^{T}}^{\prime}(U)-OPT^{T}_{k}(U)\leq 4\sum_{v^{\prime}_{j}\in
C^{*}_{1}\setminus C_{1}}N_{v^{\prime}_{j}}2^{h_{v^{\prime}_{j}}},$
$OPT^{T}_{k}(U)\geq\sum_{v_{i}\in C_{1}\setminus
C^{*}_{1}}N_{v_{i}}2^{h_{v_{i}}}.$

Also, based on our discussion on Case 1, it holds that

$N_{v^{\prime}_{j}}2^{h_{v^{\prime}_{j}}}-N_{v_{i}}2^{h_{v_{i}}}\leq 0.$

Summing them together, we have ${cost_{k}^{T}}^{\prime}(U)\leq
5OPT^{T}_{k}(U)$. ∎

###  A.2 Proof of Lemma 3.4

###### Proof.

Since the subtrees in $C_{1}$ are disjoint, it suffices to consider one
subtree with root $v$. With a little abuse of notation, let
${cost_{1}^{T}}^{\prime}(v,U)$ denote the optimal $k$-median cost within the
point set $T(v)$ with one center in 2-HST:

$\displaystyle{cost_{1}^{T}}^{\prime}(v,U)=\min_{x\in T(v)}\sum_{y\in
T(v)}\rho^{T}(x,y),$ (9)

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 ${cost_{1}^{T}}^{\prime}(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\neq w$. We have the following two observations:

  * •

Since one needs to pick a leaf of $T(u)$ to minimize
${cost_{1}^{T}}^{\prime}(v,U)$, we have
${cost_{1}^{T}}^{\prime}(v,U)\geq\sum_{x\in ch(v),x\neq u}N_{x}\cdot
2^{h_{x}}$ where $ch(u)$ denotes the children nodes of $u$.

  * •

By our greedy strategy, $cost^{T}_{1}(v,U)\leq\sum_{x\in ch(u)}N_{x}\cdot
2^{h_{x}}\leq{cost_{1}^{T}}^{\prime}(v,U)+N_{u}\cdot 2^{h_{u}}$.

Since $h_{u}=h_{w}$, we have

$2^{h_{u}}\cdot(N_{u}-N_{w})\leq 0,$

since our algorithm picks subtree roots with highest scores. Then we have
$cost^{T}_{1}(v,U)\leq{cost_{1}^{T}}^{\prime}(v,U)+N_{w}\cdot 2^{h_{w}}\leq
2{cost_{1}^{T}}^{\prime}(v,U)$. Since the subtrees in $C_{1}$ are disjoint,
the union of centers for $OPT_{1}^{T}(v,U)$, $v\in C_{1}$ forms the optimal
centers with size $k$. Note that, for any data point $p\in U\setminus C_{1}$,
the tree distance $\rho^{T}(p,f)$ for $\forall f$ that is a leaf node of
$T(v)$, $v\in C_{1}$ 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. ∎

###  A.3 Proof of Proposition 3.1

###### Proof.

It is known that the 2-HST can be constructed in $O(dn\log n)$ [Bartal, 1996].
The subtree search in Algorithm 4 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 $C_{1}$. 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 $C_{0}$ 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.
∎

###  A.4 Proof of Theorem 4.2

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
$cost_{k}^{T}(D)$, ${cost_{k}^{T}}^{\prime}(D)$ and $OPT_{k}^{T}(D)$ be
defined analogously to (5), (6) and (7), where “$y\in U$” in the summation is
changed into “$y\in D$” since $D$ is the demand set. That is,

$\displaystyle cost_{k}^{T}(D)$ $\displaystyle=\sum_{y\in D}\min_{x\in
C_{0}}\rho^{T}(x,y),$ (10) $\displaystyle{cost_{k}^{T}}^{\prime}(D,C_{1})$
$\displaystyle=\min_{|F\cap T(v)|=1,\forall v\in C_{1}}\sum_{y\in D}\min_{x\in
F}\rho^{T}(x,y),$ (11) $\displaystyle OPT_{k}^{T}(D)$
$\displaystyle=\min_{F\subset D,|F|=k}\sum_{y\in D}\min_{x\in
F}\rho^{T}(x,y)\equiv\min_{C_{1}^{\prime}}\
{cost_{k}^{T}}^{\prime}(D,C_{1}^{\prime}).$ (12)

We have the following.

######  Lemma A.2.

$cost^{T}_{k}(D)\leq 10OPT^{T}_{k}(D)+10ck\epsilon^{-1}\triangle\log n$ with
probability $1-4k/n^{c}$.

###### Proof.

The result follows by combining the following Lemma A.4, Lemma A.5, and
applying union bound. ∎

######  Lemma A.3.

For any node $v$ in $T$, with probability $1-1/n^{c}$, $|\hat{N}_{v}\cdot
2^{h_{v}}-N_{v}\cdot 2^{h_{v}}|\leq c\epsilon^{-1}\triangle\log n$.

###### Proof.

Since $\hat{N_{v}}=N_{v}+Lap(2^{(L-h_{v})/2}/\epsilon)$, we have

$Pr[|\hat{N}_{v}-N_{v}|\geq x/\epsilon]=exp(-x/2^{(L-h_{v})}).$

As $L=\log\triangle$, we have

$Pr[|\hat{N}_{v}-N_{v}|\geq x\triangle/(2^{h_{v}}\epsilon)]\leq exp(-x).$

Hence, for some constant $c>0$,

$\displaystyle Pr[|\hat{N}_{v}\cdot 2^{h_{v}}-N_{v}\cdot 2^{h_{v}}|\leq
c\epsilon^{-1}\triangle\log n]\geq 1-exp(-c\log n)=1-1/n^{c}.$

∎

######  Lemma A.4 (DP Subtree Search).

With probability $1-2k/n^{c}$, ${cost_{k}^{T}}^{\prime}(D)\leq
5OPT^{T}_{k}(D)+4ck\epsilon^{-1}\triangle\log n$.

###### Proof.

The proof is similar to that of Lemma 3.3. Consider the intermediate output of
Algorithm 4, $C_{1}=\\{v_{1},v_{2},...,v_{k}\\}$, which is the set of roots of
the minimal disjoint subtrees each containing exactly one output center
$C_{0}$. Suppose one of the optimal “root set” that minimizes (6) is
$C^{*}_{1}=\\{v^{\prime}_{1},v^{\prime}_{2},...,v^{\prime}_{k}\\}$. Assume
$C_{1}\neq C_{1}^{*}$. By the same argument as the proof of Lemma 3.3, we
consider for some $i,j$ such that $v_{i}\neq v_{j}^{\prime}$, where $v_{i}$ is
not a descendent of $v_{j}^{\prime}$ and $v_{j}^{\prime}$ is either a
descendent of $v_{i}$. By the construction of $C_{1}$, we know that
$score(v^{\prime}_{j})\leq\min\\{score(v_{i}),i=1,...,k\\}$ when
$v^{\prime}_{j}\in C^{*}_{1}\setminus C_{1}$. Consider the swap between
$C_{1}$ and $C^{*}_{1}$. By the definition of tree distance, we have
$OPT^{T}_{k}(U)\geq\sum_{v_{i}\in C_{1}\setminus
C^{*}_{1}}N_{v_{i}}2^{h_{v_{i}}}$, since $\\{T(v_{i}),v_{i}\in C_{1}\setminus
C^{*}_{1}\\}$ does not contain any center of the optimal solution determined
by $C^{*}_{1}$ (which is also the optimal “root set” for $OPT_{k}^{T}$). 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^{\prime}_{j}$), and $S^{\prime}_{j}$ be the
leaves assigned to $c^{*}_{j}$. Let $S_{j}$ denote the set of leaves in
$S^{\prime}_{j}$ whose distance to $c^{*}_{j}$ is strictly smaller than its
distance to any centers in $C_{1}$. Let $P_{j}$ denote the union of paths
between leaves of $S_{j}$ to its closest center in $C_{1}$. Let
$v^{\prime\prime}_{j}$ be the nodes in $P_{j}$ with highest level satisfying
$T(v^{\prime\prime}_{j})\cap C_{1}=\emptyset$. The score of
$v^{\prime\prime}_{j}$ is
$2^{h_{v^{\prime\prime}_{j}}}N(v^{\prime\prime}_{j})$. That means the swap
with a center $v^{\prime}_{j}$ into $C_{1}$ can only reduce $4\cdot
2^{h_{v^{\prime\prime}_{j}}}N(v^{\prime\prime}_{j})$ to
${cost_{k}^{T}}^{\prime}(U)$ (the tree distance between any leaf in $S_{j}$
and its closest center in $C_{1}$ is at most $4\cdot
2^{h_{v^{\prime\prime}_{j}}}$). We just use $v^{\prime}_{j}$ to represent
$v^{\prime\prime}_{j}$ for later part of this proof for simplicity. Summing
all the swaps over $C^{*}_{1}\setminus C_{1}$, we obtain

${cost_{k}^{T}}^{\prime}(U)-OPT^{T}_{k}(U)\leq 4\sum_{v^{\prime}_{j}\in
C^{*}_{1}\setminus C_{1}}N_{v^{\prime}_{j}}2^{h_{v^{\prime}_{j}}},$
$OPT^{T}_{k}(U)\geq\sum_{v_{i}\in C_{1}\setminus
C^{*}_{1}}N_{v_{i}}2^{h_{v_{i}}}.$

Applying union bound with Lemma A.3, with probability $1-2/n^{c}$, we have

$N_{v^{\prime}_{j}}2^{h_{v^{\prime}_{j}}}-N_{v_{i}}2^{h_{v_{i}}}\leq
2c\epsilon^{-1}\triangle\log n.$

Consequently, we have with probability, $1-2k/n^{c}$,

$\displaystyle{cost_{k}^{T}}^{\prime}(D)$ $\displaystyle\leq
5OPT^{T}_{k}(D)+4c|C_{1}\setminus C_{1}^{*}|\epsilon^{-1}\triangle\log n$
$\displaystyle\leq 5OPT^{T}_{k}(D)+4ck\epsilon^{-1}\triangle\log n.$

∎

######  Lemma A.5 (DP Leaf Search).

With probability $1-2k/n^{c}$, Algorithm 6 produces initial centers with
$cost^{T}_{k}(D)\leq 2{cost_{k}^{T}}^{\prime}(D)+2ck\epsilon^{-1}\triangle\log
n$.

###### Proof.

The proof strategy follows Lemma 3.4. We first consider one subtree with root
$v$. Let ${cost_{1}^{T}}^{\prime}(v,U)$ denote the optimal $k$-median cost
within the point set $T(v)$ with one center in 2-HST:

$\displaystyle{cost_{1}^{T}}^{\prime}(v,D)=\min_{x\in T(v)}\sum_{y\in T(v)\cap
D}\rho^{T}(x,y).$ (13)

Suppose $v$ has more than one children $u,w,...$, and the optimal solution of
${cost_{1}^{T}}^{\prime}(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\neq
w$. We have the following two observations:

  * •

Since one needs to pick a leaf of $T(u)$ to minimize
${cost_{1}^{T}}^{\prime}(v,U)$, we have
${cost_{1}^{T}}^{\prime}(v,U)\geq\sum_{x\in ch(v),x\neq u}N_{x}\cdot
2^{h_{x}}$ where $ch(u)$ denotes the children nodes of $u$.

  * •

By our greedy strategy, $cost^{T}_{1}(v,U)\leq\sum_{x\in ch(u)}N_{x}\cdot
2^{h_{x}}\leq{cost_{1}^{T}}^{\prime}(v,U)+N_{u}\cdot 2^{h_{u}}$.

As $h_{u}=h_{w}$, leveraging Lemma A.3, with probability $1-2/n^{c}$,

$\displaystyle 2^{h_{u}}\cdot(N_{u}-N_{w})$ $\displaystyle\leq
2^{h_{u}}(\hat{N}_{u}-\hat{N}_{w})+2c\epsilon^{-1}\triangle\log n$
$\displaystyle\leq 2c\epsilon^{-1}\triangle\log n.$

since our algorithm picks subtree roots with highest scores. Then we have
$cost^{T}_{1}(v,D)\leq{cost_{k}^{T}}^{\prime}(v,D)+N_{w}\cdot
2^{h_{u}}+2c\epsilon^{-1}\triangle\log n\leq
2{cost_{k}^{T}}^{\prime}(v,D)+2c\epsilon^{-1}\triangle\log n$ with high
probability. Lastly, applying union bound over the disjoint $k$ subtrees gives
the desired result. ∎

###  A.5 Proof of Theorem 4.3

###### Proof.

The privacy analysis is straightforward, by using the composition theorem
(Theorem A.1). Since the sensitivity of $cost(\cdot)$ is $\triangle$, in each
swap iteration the privacy budget is $\epsilon/2(T+1)$. Also, we spend another
$\epsilon/2(T+1)$ privacy for picking a output. Hence, the total privacy is
$\epsilon/2$ for local search. Algorithm 6 takes $\epsilon/2$ DP budget for
initialization, so the total privacy is $\epsilon$.

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 A.6 (Gupta et al. [2010]).

Assume the solution to the optimal utility is unique. For any output $o\in O$
of $2\triangle\epsilon$-DP exponential mechanism on dataset $D$, it holds for
$\forall t>0$ that

$Pr[q(D,o)\leq\max_{o\in O}q(D,o)-(\ln|O|+t)/\epsilon]\leq e^{-t},$

where $|O|$ is the size of the output set.

######  Lemma A.7 (Arya et al. [2004]).

For any set $F\subseteq D$ with $|F|=k$, there exists some swap $(x,y)$ such
that the local search method admits

$cost_{k}(F,D)-cost_{k}(F-\\{x\\}+\\{y\\},D)\geq\frac{cost_{k}(F,D)-5OPT(D)}{k}.$

From Lemma A.7, we know that when $cost_{k}(F_{i},D)>6OPT(D)$, there exists a
swap $(x,y)$ s.t.

$cost_{k}(F_{i}-\\{x\\}+\\{y\\},D)\leq(1-\frac{1}{6k})cost_{k}(F_{i},D).$

At each iteration, there are at most $n^{2}$ possible outputs (i.e., possible
swaps), i.e., $|O|=n^{2}$. Using Lemma A.6 with $t=2\log n$, for $\forall i$,

$Pr[cost_{k}(F_{i+1},D)\geq cost_{k}(F^{*}_{i+1},D)+4\frac{\log
n}{\epsilon^{\prime}}]\geq 1-1/n^{2},$

where $cost_{k}(F^{*}_{i+1},D)$ is the minimum cost among iteration
$1,2,...,t+1$. Hence, we have that as long as
$cost(F_{i},D)>6OPT(D)+\frac{24k\log n}{\epsilon^{\prime}}$, the improvement
in cost is at least by a factor of $(1-\frac{1}{6k})$. By Theorem 4.2, we have
$cost_{k}(F_{1},D)\leq C(\log n)(6OPT(D)+6k\triangle\log n/\epsilon)$ for some
constant $C>0$. Let $T=6Ck\log\log n$. We have that

$\displaystyle E[cost(F_{i},D)]$
$\displaystyle\leq(6OPT(D)+6k\epsilon^{-1}\triangle\log n)C(\log
n)(1-1/6k)^{6Ck\log\log n}$ $\displaystyle\leq
6OPT(D)+6k\epsilon^{-1}\triangle\log n\leq 6OPT(D)+\frac{24k\log
n}{\epsilon^{\prime}}.$

Therefore, with probability at least $(1-T/n^{2})$, there exists an $i\leq T$
s.t. $cost(F_{i},D)\leq 6OPT(D)+\frac{24k\log n}{\epsilon^{\prime}}$. Then by
using the Lemma A.7, one will pick an $F_{j}$ with additional additive error
$4\ln n/\epsilon^{\prime}$ to the $\min\\{cost(F_{j},D),j=1,2,...,T\\}$ with
probability $1-1/n^{2}$. Consequently, we know that the expected additive
error is

$\displaystyle 24k\triangle\log n/\epsilon^{\prime}+4\log
n/\epsilon^{\prime}=O(\epsilon^{-1}k^{2}\triangle(\log\log n)\log n),$

with probability $1-1/poly(n)$.

∎

##  Appendix B Extending 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.

$\displaystyle cost_{km}^{T}(U)$ $\displaystyle=\sum_{y\in U}\min_{x\in
C_{0}}\rho^{T}(x,y)^{2},$ (14)
$\displaystyle{cost_{km}^{T}}^{\prime}(U,C_{1})$ $\displaystyle=\min_{|F\cap
T(v)|=1,\forall v\in C_{1}}\sum_{y\in U}\min_{x\in F}\rho^{T}(x,y)^{2},$ (15)
$\displaystyle OPT_{km}^{T}(U)$ $\displaystyle=\min_{F\subset
U,|F|=k}\sum_{y\in U}\min_{x\in
F}\rho^{T}(x,y)^{2}\equiv\min_{C_{1}^{\prime}}\
{cost_{km}^{T}}^{\prime}(U,C_{1}^{\prime}).$ (16)

For simplicity, we will use ${cost_{km}^{T}}^{\prime}(U)$ to denote
${cost_{km}^{T}}^{\prime}(U,C_{1})$ if everything is clear from context. Here,
$OPT_{km}^{T}$ (16) is the cost of the global optimal solution with 2-HST
metric.

######  Lemma B.1 (Subtree search).

${cost_{km}^{T}}^{\prime}(U)\leq 17OPT^{T}_{km}(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^{\prime}_{j}$), and $S^{\prime}_{j}$ be the
leaves assigned to $c^{*}_{j}$ in optimal k-means clustering in tree metric.
Let $S_{j}$ denote the set of leaves in $S^{\prime}_{j}$ whose distance to
$c^{*}_{j}$ is strictly smaller than its distance to any centers in $C_{1}$.
Let $P_{j}$ denote the union of paths between leaves of $S_{j}$ to its closest
center in $C_{1}$. Let $v^{\prime\prime}_{j}$ be the nodes in $P_{j}$ with
highest level satisfying $T(v^{\prime\prime}_{j})\cap C_{1}=\emptyset$. The
score of $v^{\prime\prime}_{j}$ is
$2^{h_{v^{\prime\prime}_{j}}}N(v^{\prime\prime}_{j})$. That means the swap
with a center $v^{\prime}_{j}$ into $C_{1}$ can only reduce $(4\cdot
2^{h_{v^{\prime\prime}_{j}}})^{2}N(v^{\prime\prime}_{j})$ to
${cost_{km}^{T}}^{\prime}(U)$. We just use $v^{\prime}_{j}$ to represent
$v^{\prime\prime}_{j}$ for later part of this proof for simplicity. By our
reasoning, summing all the swaps over $C^{*}_{1}\setminus C_{1}$ gives

${cost_{km}^{T}}^{\prime}(U)-OPT^{T}_{km}(U)\leq\sum_{v^{\prime}_{j}\in
C^{*}_{1}\setminus C_{1}}N_{v^{\prime}_{j}}\cdot(4\cdot
2^{h_{v^{\prime}_{j}}})^{2},$ $OPT^{T}_{km}(U)\geq\sum_{v_{i}\in
C_{1}\setminus C^{*}_{1}}N_{v_{i}}(2^{h_{v_{i}}})^{2}.$

Also, based on our discussion on Case 1, it holds that

$N_{v^{\prime}_{j}}2^{h_{v^{\prime}_{j}}}-N_{v_{i}}2^{h_{v_{i}}}\leq 0.$

Summing them together, we have ${cost_{km}^{T}}^{\prime}(U)\leq
17OPT^{T}_{km}(U)$. ∎

Next, we show that the greedy leaf search strategy (Algorithm 5) only leads to
an extra multiplicative error of 2.

######  Lemma B.2 (Leaf search).

$cost^{T}_{km}(U)\leq 2{cost_{km}^{T}}^{\prime}(U)$.

###### Proof.

Since the subtrees in $C_{1}$ are disjoint, it suffices to consider one
subtree with root $v$. With a little abuse of notation, let
${cost_{1}^{T}}^{\prime}(v,U)$ denote the optimal $k$-means cost within the
point set $T(v)$ with one center in 2-HST:

$\displaystyle{cost_{1}^{T}}^{\prime}(v,U)=\min_{x\in T(v)}\sum_{y\in
T(v)}\rho^{T}(x,y)^{2},$ (17)

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 ${cost_{1}^{T}}^{\prime}(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\neq w$. We have the following two observations:

  * •

Since one needs to pick a leaf of $T(u)$ to minimize
${cost_{1}^{T}}^{\prime}(v,U)$, we have
${cost_{1}^{T}}^{\prime}(v,U)\geq\sum_{x\in ch(v),x\neq
u}N_{x}\cdot(2^{h_{x}})^{2}$ where $ch(u)$ denotes the children nodes of $u$.

  * •

By our greedy strategy, $cost^{T}_{1}(v,U)\leq\sum_{x\in
ch(u)}N_{x}\cdot(2^{h_{x}})^{2}\leq{cost_{1}^{T}}^{\prime}(v,U)+N_{u}\cdot(2^{h_{u}})^{2}$.

Since $h_{u}=h_{w}$, we have

$2^{h_{u}}\cdot(N_{u}-N_{w})\leq 0,$

since our algorithm picks subtree roots with highest scores. Then we have
$cost^{T}_{1}(v,U)\leq{cost_{1}^{T}}^{\prime}(v,U)+N_{w}\cdot(2^{h_{w}})^{2}\leq
2{cost_{1}^{T}}^{\prime}(v,U)$. Since the subtrees in $C_{1}$ are disjoint,
the union of centers for $OPT_{1}^{T}(v,U)$, $v\in C_{1}$ forms the optimal
centers with size $k$. Note that, for any data point $p\in U\setminus C_{1}$,
the tree distance $\rho^{T}(p,f)$ for $\forall f$ that is a leaf node of
$T(v)$, $v\in C_{1}$ 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 4), which is a natural combination of Lemma B.1 and Lemma B.2.

######  Theorem B.3 (HST initialization).

$cost^{T}_{km}(U)\leq 34OPT^{T}_{km}(U)$.

We have the following result based on Lemma 3.2.

######  Theorem B.4.

In a general metric space,

$E[cost_{km}(U)]=O(\min\\{\log n,\log\triangle\\})^{2}OPT_{km}(U).$