A Fast Approximation Scheme for Low-Dimensional $k$-Means
Vincent Cohen-Addad












A Fast Approximation Scheme for Low-Dimensional k-Means
    Vincent Cohen-Addad
 University of Copenhagen
    
=======================================================


  We consider the popular k-means problem in d-dimensional
  Euclidean space. Recently Friggstad, Rezapour, Salavatipour [FOCS'16]
  and Cohen-Addad, Klein, Mathieu
  [FOCS'16] showed that the standard local search algorithm yields a 
  (1+)-approximation
  in time (n · k)^1/^O(d), giving the first 
  polynomial-time approximation
  scheme for the problem in low-dimensional Euclidean space.
  While local search achieves optimal approximation guarantees, it is not 
  competitive with the state-of-the-art heuristics such as the famous 
  k-means++ and D^2-sampling algorithms.


  In this paper, we aim at bridging the gap between theory and practice by
  giving a (1+)-approximation algorithm for low-dimensional 
  k-means 
  running in time n · k · (log n)^(d^-1)^O(d),
  and so  matching the running time of the k-means++ and 
  D^2-sampling heuristics 
  up to polylogarithmic factors.
  We speed-up
  the local search approach by making a non-standard use
  of randomized dissections that allows to 
  find the best local move efficiently using a quite simple
  dynamic
  program. We hope that our techniques could help
  design better
  local search heuristics for geometric problems.




§ INTRODUCTION

The k-means objective is arguably the most popular clustering 
objective among practitioners. While originally motivated by applications
in image compression, the k-means problem has proven to be a 
successful objective to optimize in order to pre-process and extract
information from datasets. 
Its most successful applications are now stemming from 
machine learning problems like for example learning 
mixture of Gaussians, Bregman clustering, or DP-means <cit.>. 
Thus, it has become a classic problem in both machine learning
and theoretical computer science.

Given a set of points in a metric space, the k-means
problem asks for a
 set of k points, called centers,
that minimizes the sum of the squares of the 
distances of the points to their closest center.

The most famous algorithm for k-means is arguably the 
Lloyd[Also referred to as Lloyd-Forgy] heuristic introduced
in the 80s <cit.> and sometimes 
referred to as “the k-means algorithm”.
While this algorithm is very competitive in
practice and yields empirically 
good approximate solutions to real-world inputs,
it is known that its running time can be exponential in the input 
size and that it can return arbitrarily bad solutions
in the worst-case (see <cit.>). This induces a gap between
theory and practice.

Thus, to fix this unsatisfactory situation,
Arthur and Vassilvitskii <cit.> have designed a variant
of the  Lloyd heuristic, called 
k-means++,


and proved that it achieves an O(log k) approximation.

The k-means++ algorithm has now become 
a standard routine that is part of several machine learning libraries
(<cit.>) and is widely-used in practice.
While this has been a major step for reducing the gap between
theory and practice, it has remained an important problem as to
conceive algorithms with nearly-optimal approximation guarantee.

Unfortunately, the k-means problem is known to be APX-Hard
even for (high dimensional) Euclidean inputs <cit.>. 
Hence, to design
competitive approximation schemes it is needed to restrict our
attention to classes of “more structured” inputs that are important 
in practice.
The low-dimensional Euclidean inputs form a class of 
inputs that naturally arise in image processing and
machine learning (see examples of <cit.> or
in <cit.>).
Thus, finding a polynomial
time approximation scheme (PTAS) for 
O(1)-dimensional
Euclidean inputs of k-means  has been
an important research problem for the last 20 years 
since the seminal work of <cit.>. 


Recently, Friggstad et
al. <cit.>
and Cohen-Addad et al. <cit.> both showed that the classic
local search heuristic with neighborhood of magnitude (d/)^O(d)
achieves a (1+)-approximation.
While this has been an important result for the theory community, it 
has a much weaker impact for practitioners since the running time
of the algorithm is n^(d/)^O(d). 
Therefore, to reduce the gap between
theory and practice, it is natural to ask for
near-optimal approximation algorithms with competitive
running time.
This is the goal of this paper.




  
Fast local search techniques are important.
The result of Friggstad et al. and Cohen-Addad et al. has been 
preceded by several result showing that local search achieves good
approximation bound or even exact algorithms for various problems 
(see <cit.>).
Furthermore, there is a close relationship
between local search and clustering showing that the standard local
search heuristics achieve very good approximation guarantees
in various settings 
(in addition to the two 
aforementioned papers, 
see <cit.>).
Moreover local search approaches are extremely popular in practice 
because they are easy to tune, easy to implement, and easy
to run in parallel. 

Thus, it has become part of the research agenda of the theory community
to develop fast local search approaches while preserving the 
guarantees on the quality of the output (see <cit.>). 





 §.§ Our Results

We show that our fast local search algorithm 
(Algorithm <ref>)
yields a PTAS for the slightly more
general variants of the k-means problem where centers can 
have an opening cost
(a.k.a. weight).


  
  There exists a randomized
  algorithm (Algorithm <ref>) 
  that returns a
  (1+) approximation to the center-weighted d-dimensional 
  Euclidean
  k-means problem in time n · k · 
  (log n)^(d ^-1)^O(d) with probability at least 1/2.






We would like to remark that the doubly exponential dependency 
in d is needed: 
Awasthi et al. <cit.> showed that the Euclidean k-means 
is APX-Hard when d = Ω(log n).
Note that it is possible to obtain an arbitrarily small
probability of failure p > 0 by repeating the algorithm log 1/p 
times.

As far as we know, this is the first occurrence of a local search 
algorithm whose neighborhood size only impacts the running time 
by polylogarithmic factors.




 §.§ Other Related Work

The k-means problem is known to be NP-Hard, even when restricted
to inputs lying in the Euclidean plane 
(Mahajan et al.  <cit.>, and Dasgupta and Freud <cit.>) and
was recently shown to be APX-hard in Euclidean space of dimension
Ω(log n) (<cit.>).

There has been a large body of work on approximation algorithms 
for the Euclidean k-means problem (see <cit.>),
very recently
Ahmadian et al. <cit.> gave a 6.357-approximation
improving over the 9-approximation of Kanungo et al. <cit.>.





Given the hardness results, researchers have focused on 
different scenarios. There have been various 
(1+ε)-approximation algorithms when k is 
considered a fixed-parameter
(see <cit.>). Another successful approach
has been through the definition of “stable instances” to characterize
the real-world instances stemming from machine learning and data 
analysis 
(see for example 
<cit.>) 
or in the context of smoothed analysis (see for 
example <cit.>). 
In the case of low-dimensional inputs, Bandyapadhyay and
Varadarajan showed that local search with neighborhood of size
^-O(d) achieves a (1+)-approximation <cit.> when
allowed to open O( k) extra centers. As mentioned before,
this results has been improved by Friggstad et al. <cit.>
and Cohen-Addad et al. <cit.> who showed that even when
constrained to open exactly k centers, local search achieves
a (1+)-approximation.




  
Related work on local search
Local Search heuristics belong to the toolbox of all practitioners,
(see Aarts and Lenstra <cit.> for a general introduction).
As mentioned before, there is a tight connection between local search
and clustering:
Arya et al. <cit.> proved that local search with a neighborhood size of 1/ε yields a 3+2ε approximation to k-median. For the k-means problem,
Kanungo et al. <cit.> showed a similar result by proving 
that the approximation guarantee of local search for Euclidean k-means
is 9+ε.
For more applied examples of local search an clustering 
see <cit.>.
For other theoretical example of local search for clustering, 
we refer to <cit.>.



  
Related work on k-median.
The k-median problem has been widely studied.
For the best known results in terms of approximation for general
metric spaces inputs we refer to Li and Svensson <cit.> and 
Byrka et al. <cit.>. More related to our results are the 
approximation schemes for k-median in low-dimensional Euclidean
space given by Arora et al. <cit.> who gave  a
(1+)-approximation algorithm running in time
n^^-O(d). This was later improved by Kolliopoulos and 
Rao <cit.> who obtained a running time of 2^^-O(d) n·polylog n. Quite surprisingly, it is pretty unclear whether 
the techniques used by Arora et al. and Kolliopoulos and Rao could be
used to obtain a (1+)-approximation for the k-means problem; 
this has induced a 20-year gap between the first PTAS for k-median
and the first PTAS for k-means in low-dimensional Euclidean space.
See Section <ref> for more details.



 §.§ Overview of the Algorithm and the Techniques


Our proof is rather simple.
Given a solution L, our goal is to identify – in near-linear
time – a minimum
cost solution L'
such that |L - L'| + |L' - L| ≤δ for some (constant) parameter 
δ. If (L) - (L') = O(()/k),
then we can immediately apply the result of  
Friggstad et al. <cit.>
or Cohen-Addad et al. <cit.>: the solution L is 
locally optimal and its cost is at most (1+O())().
Finding L has
to be done in near-linear time since we could repeat this 
process up to Θ(k) times
until reaching a local optimum.
Hence, the crux of the algorithm is to efficiently identify L',
it proceeds as follows (see Algorithm <ref> for a full description): 


  * Compute a random recursive decomposition of L
  (see Section <ref>); 

  * Apply dynamic programming on the recursive decomposition;
  we show that there exists a near-optimal solution
  whose interface between different regions
  has small complexity.


To obtain our recursive decomposition we make a quite non-standard
use of the classic quadtree dissection techniques 
(see Section <ref>).
Indeed, the k-means problem is famous for being “resilient” to 
random quadtree approaches – this is partly why a PTAS for the 
low-dimensional k-median problem was obtained 20 years ago but the 
first PTAS for the k-means was only found last year. More precisely,
the classic quadtree approach (which works well for k-median) 
defines portals on the boundary of 
the regions of the dissection and forces the clients of a given region
that are served by a center that is in a different region 
(in the optimal solution) to make a detour through the closest portal.
This is a key property as it allows to bound the complexity of 
 between regions. Unfortunately, when dealing with squared 
distances, making a detour could result in a dramatic cost increase
and it is not clear that it could be compensated by the fact that 
the event of separating a client form its center happens with small
probability (applying the analysis of 
Arora et al. <cit.>
or Kolliopoulos et al. <cit.>). This problem comes from
the fact that some facilities of  and L might be too
close from the boundary of the dissection (and so, their
clients might have to make too important detour (relative
to their cost in ) through
the portals). We call these facility the “moat” facilities
(as they fall in a bounded-size “moat” around the boundaries).

We overcome this barrier by defining a more structured
near-optimal solution as follows: (1) when a facility of  is too close from the boundary of our dissection we simply remove it and (2) if a facility of the current solution L is 
too close from a boundary of our dissection we add it to .















Of course, this induces two problems: first we have to bound
the cost of removing the facilities of  and second we
have to show that adding the facilities of L does not result
in a solution containing more than k centers.
This is done through some technical lemmas and using the concept
of isolated facilities inherited from
Cohen-Addad et al. <cit.>. Section <ref>
shows the existence of near-optimal solution S^* whose set
of  facilities (facilities that are too close from
the boundaries) is exactly the set of  facilities
of L (and so, we already know their location and so the
exact cost of assigning a given client to such a facility).

We now aim at using the result of Friggstad et al. <cit.>
and Cohen-Addad et al. <cit.>: if for any sets
Δ_1 ⊆ L and Δ_2 ⊆ S^*,
we have (L) - (L - Δ_1 ∪Δ_2) =
O(()/k), then we have (L) ≤
(1+O()(S^*).
Thus, we provide a dynamic program (in Section <ref>) 
that allows to find the best solution S that is such that
(1) its set of 
facilities is exactly the set of  facilities
of L and (2) |S-L|+|L-S| ≤δ for some fixed
constant δ. The dynamic program simply “guess” the
approximate location of the centers of (L-S) ∪ (S-L),
we show that since each such center is far from the boundary,
its location can be approximated.























































 §.§ Preliminaries

In this article, we consider the k-means problem in a d-dimensional Euclidean space:
Given a set A of points (also referred to as clients) and candidate centers C in ^d,
the goal is to output a set S⊆ C of size k that minimizes:
∑_a ∈ A(a,C)^2, where (a,C) = min_c ∈ C(a,c). We refer to S as a set of centers or facilities.
Our results naturally extends to any objective function of the form 
∑_a ∈ A(a,C)^p for constant p. For ease of exposition, we focus on the 
k-means problem. As Friggstad et al <cit.>, we also 
consider the more general version
called weighted k-means for which, in addition of the sets A and C, we are given a weight
function w : C ↦_+ and the goal is to minimize
∑_c ∈ C w(c) + ∑_a ∈ A(a,C)^2.

A classic result of Matousek <cit.> shows that, if C = ^d, it is possible to compute in linear
time a set C' of linear size (and polynomial dependency in d and ) such that the optimal solution using the 
centers in C' cost at most (1+) times the cost of the optimal solution using
C. Hence, we assume without loss of generality that 
|C| has size linear in |A| and we
let n = |A| + |C|.



  
Isolated Facilities. We make use of the notion of isolated facilities
introduced by Cohen-Addad et al. <cit.> defined as follows.
Let _0 <1/2 be a positive constant and  and  be two solutions for the 
Euclidean k-means problem.

 
  Let _0 <1/2 be a positive number and let  and  be two solutions for the k-clustering problem with parameter p.
  Given a facility f_0 ∈ and a facility ℓ∈, 
we say that the pair
  (f_0,ℓ) is  1-1 _0-isolated if 
  most of the clients served by ℓ in  are served by 
  f_0 in , and
  most of the clients served by f_0 in  are served by 
  ℓ in ; formally, if

  
    |(ℓ) ∩(f_0)  |   ≥max{[ (1-_0)   
                                                         |(ℓ)|,;                                                        (1-_0) |(f_0)| ]}


When _0 is clear from the context we refer to 1-1 _0-isolated 
pairs as 1-1 isolated pairs.


  Let _0 <1/2 be a positive number and let  and  be two solutions for the k-clustering problem with exponent p.
  Let k̅ denote the number of facilities f of   that
  are not in a 1-1 _0-isolated region. 
  There exists a constant c and a set S_0 of facilities of  of size at least _0^3 k̅/6 that can be removed 
  from  at low cost:  (∖ S_0) ≤ (1+c ·_0) () +  c ·_0 ().




Observe  that since _0 < 1/2, each facility of  belongs to at most one isolated region.
Let  denote the facilities of  that are not
in an isolated region.  
In the rest of the paper, we will use it with _0 = ^3.




 §.§ Fast Local Search

This section is dedicated to our the description of the fast local
search algorithm for the k-means problem.
It relies on a dynamic program called
FindImprovement that allows to find the best improvement of
the current solution in time n ·poly_,d(log n).
We then show that total number of iterations of the do-while loop is 
O(k).

SE[DOWHILE]DodoWhile[1] #1
MRepeatEndRepeatrepeat
EndRepeat





§ DISSECTION PROCEDURE


In this section, we recall the classic definition of 
quadtree dissection.
For simplicity we give the definition for _2, 
the definition directly generalizes to any fixed dimension d, 
see Arora <cit.> and Arora et al. <cit.> for a 
complete description.



Our definition of quadtree is standard and follows the 
definition of <cit.>, our contribution in the 
structural properties we extract from the dissection and 
is summarized by Lemma <ref>.

Let  be the length of the bounding box of the client
set A (the smallest square containing all the points in A).  
Applying standard preprocessing techniques, 
see in <cit.> and <cit.>, 
it is possible
to assume that the points lie on a unit grid of size polynomial in the 
number of input points. This incurs an addititive error of O(/n^c) 
for
some constant c and yields = n ·poly(^-d).

We define a quadtree dissection  of a set of 
points  as follows.
A dissection of (the bounding box of)  is a recursive 
partitioning into smaller squares. We view
it as a 4-ary tree whose root is the bounding
box of the input points. Each square in the tree is partitioned into
4 equal squares, which are its children. We
stop partitioning a square if it has size < 1 (and
therefore at most one input point). It follows that the
depth of the tree is log = O(log n). 




Standard techniques show that such a quadtree can be computed 
in n·log n ·poly(^-d), 
see <cit.> for more details. The total number
of nodes of the quadtree is 
n·log n ·poly(^-d).

Given two integers a,b ∈ [0,), the (a,b)-shifted
dissection consists in shifting the x- and y- coordinates
of all the vertical and horizontal lines by a and
b respectively. For a shifted dissection, we naturally
define the level of a bounding box to be its depth in the quadtree.
From this, we define the level of a line to be the level
of the square it bounds.

For a given square of the decomposition, each boundary of 
the square defines a subline of one of the 2 lines
of the grid. It follows that each line at level i consists
of 2^i sublines of length /2^i.

Given an (a,b)-shifted quadtree dissection of a set of n points, 
and given a set of points U,  
we say that a point p of U is an i,γ-
point
if it is at distance less than γ·/2^i of a line
of the dissection that is at level i.
We say that a point p of U is a γ-
point
if there exists an i such that p is a i,γ-
point. When γ is clear from the context, we simply
call such a point a  point.


We have:

  
  
  For any p ∈ U, the probability that 
  p is a γ- point
  is at most γlog = O(γlog n).


  
  Let i be an integer in [0,…,log] 
  and consider the horizontal lines at level i
  (an analogous 
  reasonning applies to the vertical lines).
  
  By definition, the number of dissection
  lines that are at distance 
  at most γ/2^i 
  from p is
  γ/2^i.
  
  
  We now bound the probability that one of them is 
  at level i (and so the probability of p being
  at distance less than γ/2^i of a horizontal
  line of length /2^i). 
  For any such line l,
  we have:
  Pr_a[l  is at level i] = 2^i/.
  Hence, Pr_a[p is a i,γ-  point]
  ≤∑_l: (l,p) ≤γ/2^i 
  Pr_a[l  is at level i] ≤γ.
  The lemma follows by taking a union bound over all i.























We now consider an optimal solution  and
any solution L. In the following, we will focus
on γ- centers of L and  for
γ = ^13/log n. In the rest of the
paper, γ is fixed to that value and so
γ-
centers are simply called  centers.

Let ι the facilities of  that are not
1-1 isolated.

We define a weigth function w̃ : 
L ∪∪↦_+ as follows.
For each facility s ∈ι 
we define (s) as the sum of w(s) and the cost
of serving all the clients served by s in 
by the closest facility ℓ in L plus w(ℓ).
For each facility s ∈ L 
we define (s) = w(s)
+ ∑_c served by s in L(c,s).
Similarly, for each s ∈ - ι,
we let (s) = w(s)
+ ∑_c served by s in (c,s).












We show:

  
  There exists a constant c_0 such that
  ∑_s ∈ι(s)
  ≤  c_0 (() + (L)).  


  Consider a facility s ∈ι
  and the closest facility l in L that (1)
  serves in L at least one client that is served by s in 
  and that (2) minimizes the
  following quantity:
  η = min_c (c, l)^2 + (c, s)^2. Let c^* be a client
  that minimizes the quantity (c,l)^2+(c,s)^2.
  Let N(s) be the set of all clients
  served by s in . We have:
  |N(s)| η≤∑_c ∈ N(s)(c,L)^2+(c,)^2.
  
  We have that the total cost of sending all the clients in N(s)
  is at most (by triangle inequality):
  ∑_c ∈ N(s) ((c, s) + (s,l))^2 ≤∑_c ∈ N(s) ((c, s) + (s,c) + (c,l))^2.
  Note that there exists a constant c_0 such that the above sum is at most
  c_0 ∑_c ∈ N(s)(c, s)^2 + (s,c^*)^2 + (c^*,l)^2 and so at
  most
  c_0 (|N(s)|η + ∑_c ∈ N(s)(c, )^2). The lemma follows
  by combining with the
  above bound on |N(s)|  η.


In the following we denote by ϕ : ι↦ L
the mapping from
each non-isolated facility of  to its
closest facility in L. We define ι_L to
be the set of non-isolated facilities of L.

We define Event (L ∪, w̃)
as follows:


  * The set of 
  centers S_1 of L ∪ι is such that
  (S_1) = ∑_c ∈ S_1(c) 
  ≤^9 ∑_c ∈ L ∪ι(c)
  = ^9 (L ∪ι), and

  * The set of 
  centers S_1 of ι_L ∪ι is such that 
  |S_0| ≤^9 |ι_L ∪ι| = ^9k̅
  (recall that k̅ is the number of non 1-1-isolated
  facilities of  (and so of L as well)).

The following lemma follows from Lemma <ref>,
applying Markov's
inequality and taking a union bound over the probability of failures of 
property (1) and (2).


  
  The probability that Event (L ∪, w̃)
  happens is at least 1/2.


  We apply Lemma <ref> and obtain that any
  element of L ∪
  is a  center
  with probability O(ρ^-1log n). Since 
  ρ^-1 =  (c ^-12log n)^-1,
  we obtain that this probability
  is at most ^12.
  Thus, taking linearity of expectation
  we have that the expected size 
  of S_0 is at most ^12 |S ∪ι| and
  that the expected
  value of (S_0) is at most ^12(S ∪ι).
  Applying Markov's inequality to 
  obtain concentration bounds on both quantities and
  then taking a union 
  bound over the probabilities of failure yields the lemma.



We finally conclude this section with some
additional definitions
that are used in the following sections.
We define a basic region of a decomposition
of a to be a region of the dissection that contains exactly
1 points of L. The other squares of the decomposition 
are simply called regions. 


































§ A STRUCTURED NEAR-OPTIMAL SOLUTION















This section is dedicated to the following proposition.

  
  Let L be any solution.
  Let _L be a random quadtree dissection 
  of L as per Sec. <ref>. 
  Suppose Event (L ∪,) happens. 
  Then there exists
  a constant c and a solution S^* 
  
  of cost at most
  (1+c ·) ()  + c·(L)
  and such that the set of  centers
  of S^* is equal to the set of  centers of L.


We prove Propositon <ref> by
explicitly constructing S^*.
We iteratively modify  in four main steps:


  * Modify  by replacing f_0 by ℓ_0 for 
  for each 1-1 isolated pair (f_0,ℓ_0) 
  where ℓ_0 or f_0 is a  center.
  This yields a near-optimal solution S_0 
  (Lemma <ref>).





  * Replace in  each  center s that is in
  ι by ϕ(s) (as per Section <ref>). This yields a near-optimal
  solution S_1 (Lemma <ref>).

  * Apply Theorem <ref> (Theorem III.7 in <cit.>) 
to obtain a near-optimal 
  solution S_2 that has at most k -c_2 ^9·k
  where k is the number of facilities of  that
  are not 1-1 isolated (Lemma <ref>).

  * Add the  centers of L that are non
  1-1 isolated to S_2. This yields a near-optimal 
  solution S_3 that has at most
  k centers.


See Section <ref> 
for a detailed proof.




 §.§ 1-1 Isolated Pairs


We start from  and for each 1-1 isolated
facility (f_0,ℓ_0),
f_0 ∈, ℓ_0 ∈ L,
where ℓ_0 or
f_0 is a  center, we replace 
f_0 by ℓ_0 in . 

This results in a solution S_0 whose structural properties
are captured by the following lemma. For any f_0 ∈
(resp. ℓ_0 ∈), let (f_0) (resp. 
(ℓ_0)) be 
the set of clients served by f_0 in  (resp. 
the set of clients served by ℓ_0 in L).















  
  Assuming Event (L∪,) happens, 
  there exists a constant c_0 such that 
  (S_0)≤ (1+c_0 ·) () + c_0 ·(L). 


  Since Event (L∪,) happens,
  we have by Lemma <ref> 
  that  the total opening cost of the  centers plus
  the total service cost induced by the clients served by the 
   centers is 
  
  
  bounded by c_1 ·· ((L) + ())
  for some constant c_1. More formally, we can write:
  
    ∑_(f_0,ℓ_0): 1-1 
          isolated pair and ℓ_0 or f_0 is a  center 
      (w(ℓ_0)  + ∑_a ∈(ℓ_0)(a, ℓ_0)) ≤ 
      c_1 ·· ((L) + ())


  Thus, we need to bound the cost for the clients that are
  in (f_0) - 
  (ℓ_0), for each 1-1 isolated pair f_0,ℓ_0 where 
  ℓ_0 or f_0 is a  center.
  
  Consider such a pair f_0,ℓ_0.
  We bound the cost of the clients served by f_0 in  
  by the cost of rerouting them toward ℓ_0. 

  We can thus write for each such client c:
  (c,ℓ_0)^2 ≤ ((c, f_0) + (ℓ_0,f_0))^2.
  Also,
  (c,ℓ_0)^2 ≤ (1+)^2(c, f_0)^2 + 
  (1+^-1)^2(ℓ_0,f_0))^2. Note that (c,f_0)^2 is
  the cost paid by c in .
  Thus we aim at bounding (ℓ_0,f_0).

  Applying the triangle inequality, we obtain
  (ℓ_0,f_0)^2 ≤ ((ℓ_0,c_1) + (c_1,f_0))^2, 
  for any c_1 in (f_0) ∩(ℓ_0).
  Hence,
  
    (ℓ_0,f_0)^2≤1/|(f_0) ∩(ℓ_0)|∑_c_1 ∈(f_0) ∩(ℓ_0) 
      ((ℓ_0,c_1) + (c_1,f_0))^2.

  Now,   
  
    ∑_c_1 ∈(f_0) ∩(ℓ_0) ((ℓ_0,c_1) + (c_1,f_0))^2
      
      
         ≤ 3∑_c_1 ∈(f_0) ∩(ℓ_0)(c_1,ℓ_0)^2 + 
        (c_1,f_0)^2.



  Combining, we obtain that the total cost for the clients
  in (f_0) - (ℓ_0) is at most
  
    ∑_c ∈(f_0) - (ℓ_0)(c,ℓ_0)^2    ≤
        (1+)^2 (∑_c ∈(f_0) - (ℓ_0)(c,f_0)^2 ) + 
         |(f_0) - (ℓ_0)|(1+^-1)^2(ℓ_0,f_0)^2
       ≤ (1+)^2 (∑_c ∈ N(f_0) - N(ℓ_0)(c,f_0)^2 )
          
       + (1+^-1)^2 
          |(f_0) - (ℓ_0)|/|(f_0) ∪(ℓ_0)|· 3∑_c_1 ∈(f_0) ∩(ℓ_0)(c_1,f_0)^2 + 
              (c_1,ℓ_0)^2

  The lemma follows from applying the definition of 1-1 isolation:
  |(f_0) - (ℓ_0)|/|(f_0) ∪(ℓ_0)| ≤^3, and
  summing over all 1-1 isolated pair.












 §.§ Replacing the Moat Centers of 

In this section, we consider the solution S_0 and define
a solution S_1 whose set of  centers is a
subset of the  centers of L. Namely:

  
  There exists a constant c_1 and a solution S_1 such that
  the  centers of S_1 are a subset of the  centers
  of L and (S_1) ≤ (1+c_1 ) (S_0) +
  c_1(L).


  Note that by Lemma <ref>,
  all the 1-1 isolated facilities
  of S_0 that are  centers are also in L.
  Thus, we focus on the  centers of S_0 that
  are not 1-1-isolated (and so, by definition, in ι).
  
  We replace each center s ∈ι by the center
  ϕ(s) ∈ L (as per Section <ref>):
  the bound on the cost follows
  immediately from the definition of  and
  by combining Lemma <ref> and
  Lemma <ref>.




 §.§ Making Room for Non-Isolated Moat Facilities

We now consider the solution S_1
described in the previous section
that satisfies the condition of Lemma <ref>.
The following lemma is a direct corollary of 
Theorem <ref> (from <cit.>).


  
  Given S_1 and L, there exists a solution S_2 ⊆ S_1 
  and constants c_2,c_3 
  such that 
  
  
  * |S_2| ≤ k - c_2 ·^9 k,
    where k is the number of facilities that 
    are not 1-1 isolated in  (or in L it is the same number), 
    and
  
  * (S_2) ≤ (1+c_3 ·) () + 
    c_3 ··(L).
  




 §.§ Adding the  Non-Isolated Moat Facilities
and Proof of Proposition <ref>


We now consider the set of centers S_3 consisting of 
the centers in S_2
and the non-1-1-isolated  centers of 
L.
























  Combining
  Lemmas <ref>, <ref>,
  and <ref>
  yields the bound on the cost of S_3.

  By definition of S_3 and applying
  Lemmas <ref> and <ref>,
  we have that the set of   centers of S_3
  is exactly the set of  centers of L.


  By Lemma <ref> (and because
  Event (L ∪,) happens), we have that
  the total number of centers in S_3 is at most
  k.
  

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  






§ PROOF OF THEOREM <REF>




We summarize:
By Proposition <ref>, we have that there exists
a near-optimal solution S^* whose set of  centers
is the set of  centers of L. 
By Proposition <ref>, we have that, 
FindImprovement identifies a solution S' 
that is δ-close w.r.t. L, whose set of
 centers is the set of  centers
of L,
and such that (L)-(S') ≥ (1-) ((L) - 
(_δ)), where _δ is the minimum cost 
solution S such that |S - L| + |L - S| ≤δ.
We now argue  that:

If FindImprovement outputs a solution S_4
such that (L) - (S_4) ≤()/k 
then there exists a constant c^* such that 
(L) ≤ (1+c^* ·)().



  
  
  






Assuming (L) - (S_4) ≤()/k implies 
by Proposition <ref> that 
(L) - (_δ) ≤ 2()/k, 
since < 1/2.
Now, consider the solution S^* defined
in Section <ref>. 


By Theorem 1
in <cit.> 


(see also <cit.> for 
a slightly better dependency in 
in the unweighted case), if for any pair of sets
Δ_1 ⊆ L, Δ_2 ⊆ S_2 such that 
|Δ_2| ≤ |Δ_1| = (d ^-1)^O(d), we have 
(L) - (L - Δ_1 ∪Δ_2) ≤/k,
then there exists a constant c_6 such that 
(L) ≤ (1+c_6 ·) (S_2).



To obtain such a bound we want to apply
Proposition <ref>, and so we need to show that
any such solution M = L - Δ_1 ∪Δ_2
is such that its 
 centers are the  centers of L. This follows
immediately from Proposition <ref>: the
 centers of S^* are the  centers of L.















Therefore, we can apply Proposition <ref> and we have
(L) ≤ (1+c^* ·) (),
for some constant c^*.

We now bound the running time of 
Algorithm <ref>.


  
  The running time of Algorithm <ref> is at most 
  n · k · (log n)^(d^-1)^O(d).


  
  
  
  By Proposition <ref>, 
  we only need to bound the number of iterations
  of the  loop (lines <ref> 
  to <ref>) of 
  Algorithm <ref>. 
  Let (S_0) denotes the cost of the initial solution.
  The number of iterations of the  loop is
  
    log((S_0) / ())/log(1/1-1/k).

  Assuming (S_0) ≤ O(), we have that 
  the total number of iterations is at most
  O(k). To obtain 
  (S_0) ≤  O() it is possible 
  to use the 
  algorithm of Guha et al. <cit.>
  which outputs an O(1)-approximation 
  in time n · k ·polylog(n), 
  as a preprocessing step (for the computations
  at line <ref>), without
  increasing the overall running time.


We conclude by bounding the probability of failure:
By lemma <ref>, Event  happens with 
probability at least 1/2. Since the random dissection is repeated 
independently c ·log(k) times, 
the probability of failure of Event 
for a given iteration of the while loop is at most 2^-c·log k
= k^-c. Now, by Lemma <ref>, the do-while loop 
is repeated
a total of at most O(k) times, thus the probability of failure
is at most O(k^-c+1).

















































  






  







  


  










  



  











§ A DYNAMIC PROGRAM TO FIND THE BEST IMPROVEMENT


For a given solution L, we define a solution L' to 
be δ-close from L if |L-L'| + |L'-L| ≤δ. Let 
_δ denote the cost of the best solution that is 
δ-close from L and whose set of  centers
is the set of  centers of L.
In the following we refer to this solution
by the best δ-close solution.

As a preprocessing step, we round the weights of the centers
to the 
closest (1+)^i /n, for some integer i.
It is easy to see
that this only modify the total value by a factor (1+).

For each region R, we define the center of the region c_R to
be the center of the square R.
For each point p that is outside of R and at distance at least
∂ R/log n of R, we define the coordinates
of p w.r.t. c_R as follows. Consider the coordinates of the vector 
c_R p,
rounded to the closest (1+/log n)^i ^14∂ R/log n,
for some integer i. Let c_R p be the
resulting list of coordinates. Let s
the point such that the coordinates of the vector c_R s
are equal to the coordinates of c_R p.
We define the coordinates of p rounded w.r.t. R to be the coordinates of s.

For each region R we also define the grid G_R of R as the d-dimensional grid of
size 2log n/^14× ...× 2log n/^14 on R. Note that the distance between
two consecutive points of the grid is ^14∂ R/2log n.
For each point p that is inside of R, we define the coordinates of 
p rounded w.r.t. R to be the coordinates of the closest grid point.

The following fact follows from the definition and
recalling that the region sizes
are in the interval [1,poly(n)].

  
  For any region R,
  the number of different coordinates rounded w.r.t. to R is at most 
  O((log n/^14)^2d.


We now describe the dynamic program. Each entry of the table is
defined
by the following parameters:


  * a region R,

  * a list of the
  rounded coordinates w.r.t. R of the centers
  of L - _δ and _δ - L,

  * a list of the (rounded) weights of the centers of L - _δ and _δ- L.

  * a boolean vector of length δ indicating
  whether ith center
  in the above lists is in L-_δ (value 0)or 
  in _δ-L (value 1).


The following fact follows from the definition and Fact <ref>.

  The total number of entries that are parameterized by region R is at most
  (log n/^14)^O(d ·δ).


We now explain how to fill-up the table. We maintain the
following constraint when we compute a (possibly partial)
solution L': there is no center of (L' - L) ∪ (L-L') that
is . Under this constraint, we proceeds
as follows, starting with the basic regions which define
the base-case of our DP.
The base-case regions contains only one single
candidate center. Hence, the algorithm proceeds as follows:
it fills up table entries that are parameterized by:


  * any boolean vector of length δ.

  * the rounded coordinates of the unique candidate center
  inside R and a set of
  δ-1 rounded coordinates for the centers of (L-_δ) ∪ (_δ-L) outside R, or
  
  
  
  

  * a set of
  δ rounded coordinates for the centers of (L-_δ) ∪ (_δ-L) outside R.

Additionally, we require that the boolean vector is consistent
with the rounded coordinates: the candidate center inside
R is already in L if and only if
its corresponding boolean entry is 0.

It iterates over all possible rounded coordinates for the at most δ centers
of (L-_δ) ∪ (_δ-L) outside R and for each of possibility, it computes
the cost. Note that this can be done in time n ·δ.

We now consider the general case which consists in merging table
entries of child regions.
Fix a table entry parameterized by a region R, and
the rounded coordinates of the centers of (L-_δ) ∪ (_δ-L). 

We define which 
tables entries of the child regions are compatible given the rounded
coordinates.



For a table entry of a child region R_1, with the rounded coordinates of the centers
of (L-_δ) ∪ (_δ-L), we require the following for all center c_0 ∈ (L-_δ) ∪ (_δ-L).
Denote by c̃_0^1 the coordinates of c_0 ∈ (L-_δ) ∪ (_δ-L)
rounded w.r.t. c_R_1, namely its values in the table entry
for R_1. Let c̃_0^R denote
its rounded coordinate w.r.t. R, namely its values
in the table entry for R.


We require:


  * If c̃_0^R is outside R, we say that the table entries are compatible for c_0
  if the coordinates of the vector c_R c̃_0^1 and
  c_R c̃_0^R are all within a (1±/log n) factor.

  * If c_0 is inside R, we say that the table entries are compatible for c_0
  if the point of the grid G_R that is the closest to c̃_0^1 is c̃_0^R.

  * the entries corresponding to c_0 in the boolean vectors
  are the same.


The following lemma follows immediately from the above facts and the definition.

   
  The running time of the dynamic program is n (log n/)^O(dδ).


We now turn to the proof of correctness. For a given region R, and a δ-close
solution S
we define the table entry of R induced by S to be the table entry parameterized
by R and the coordinates of the centers of (L-S) ∪ (S-L) rounded w.r.t. R.

  
  Consider the best δ-close 
  solution _δ. For any level i of the quadtree
  dissection, for any region R at level i,
  we have that the table entry induced by
  _δ has cost at most
  ∑_c ∈ R ((1+/log n)^i (c, _δ))^2.
  
  
  
  
  
  
  
  


  Observe that we consider a δ-close 
  solution that has the same set of  centers than L. Hence, if a client is served by a  center in _δ, we know exactly the position of this center (as it is also in L and so cannot be removed).
  Thus, for any region R, the set of clients in R is served
  by either a center in R
  or a center at distance at least ^13∂ R/log n from
  the boundary of R or a center of L that is a  center.
  
  We now proceed by induction. We consider the base case:
  let R be a region at the maximum level.
  Consider the table entry induced by _δ.
  We claim that for each client c in R, the cost induced by the solution for this table entry is
  at most (1+/log n) (c,_δ). Indeed, since _δ has the same set of  centers
  each client that is served by a center outside of R that
  is at distance at most ^13∂ R /log n from the boundary is served by a  center of L and so there is no approximation in its service cost. Each client that is served by a center of _δ - L is at distance at least ^13∂ R/log n
  and so, the error induced by the rounding is at most (c,_δ)/log n.
  Finally,
  the cost for the clients in R served by the unique center of R (if there is one) is exact.

  Thus, assume that this holds up to level i-1. Consider a region R at level i and the
  table entry induced by _δ. The inductive hypothesis implies that for each of the
  table entries of the child regions that are induced by _δ, the cost for the clients
  in each subregion is at most (1+/log n)^i-1∑_c ∈ R(c,_δ)^2.
  
  By definition, we have that each client of R that is served in _δ by a center
  that is outside R is at distance at least ^13∂ R/log n or a  center of L (and so the distance is
  known exactly).
  Thus the rounding error incurred for the cost of the clients
  of R served by a center outside R is at most
  (1+/log n)^i(c, _δ).

  We now turn to the rounding error introduced 
  for
  the centers that are inside R. Let c
  be such a center.
  We have that the error introduced is at most the distance
  between two consecutive grid points
  and so at most ^14∂ R/2log n. Now, observe that
  again because _δ shares the
  same  centers
  than L,  each client a pf a
  child region R_1 that suffers some rounding
  error and that is served by c,
  is at distance at least ^13∂ R_1/log n
  from c and so, combining with the inductive hypothesis
  the error incurred is at most (1+/log n)^i(a,c).


Proposition <ref> follows from combining Lemmas <ref> and <ref>.
 
  
  Let >0 be a small enough constant.       
  Let L be a solution and  be a decomposition
  and  be the  centers of L.
  The dynamic program FindImprovement output a solution 
  _δ such that 
  (L) - (_δ) ≥ 
  (1-)((L) - (_δ)),
  where _δ is a minimum-cost
  solution that is δ-close
  from L and
  whose set of  centers is the set of
   centers of L.
  Its running time is
  n · (log n/^14)^O(d δ).




























































































































































































































































































































































  












































































  




















  










  






































  













  






















































































































abbrv

