An efficient algorithm for the single facility location problem with polyhedral norms and disk-shaped demand regions

The single facility location problem with demand regions seeks for a facility location minimizing the sum of the distances from n demand regions to the facility. The demand regions represent sales markets where the transportation costs are negligible. In this paper, we assume that all demand regions are disks of the same radius, and the distances are measured by a rectilinear norm, e.g. ℓ1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _1$$\end{document} or ℓ∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _\infty $$\end{document}. We develop an exact combinatorial algorithm running in time O(nlogcn)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(n\log ^c n)$$\end{document} for some c dependent only on the space dimension. The algorithm is generalizable to the other polyhedral norms.


Introduction
In the conventional facility location problem each customer is associated with a single point in the space. Often in applications, the number of points is too large to consider them individually and the clients are combined in demand regions, see Dinler et al. [6]. In the basic problem considered in the present paper, we assume the transportation costs to a demand region are proportional to the shortest distance between any point in the region and the facility, also referred to as a center. The transportation costs within the demand region are negligible. While Brimberg and Wesolowsky [3] have used the same cost function with the regions being polytopes, we consider a special case where the demand regions are disks (balls) of radius R ≥ 0 and the distances are measured by a rectilinear norm, e.g., 1 or ∞ . Notice, under 1 -and ∞ -norm a disk shape is not a conventional disk but a diamond and a square, respectively.
A special case of the problem with R = 0 and distances defined by the Euclidean norm 2 is the well-known geometric 1-median problem. It has been shown that for the geometric 1-median problem, neither a closed form formula, nor an exact algorithm involving only ruler and compass can be constructed, see Bajaj [1]. Kumar et al. [9] presented a randomized sampling algorithm, which for any ε > 0 with probability at least 1/2 finds a (1 + ε)-approximation to the geometric k-median problem in where d is the space dimension and n is the number of points. Very recently, Cohen et al. [5] presented a deterministic algorithm computing the geometric median to arbitrary precision in nearly linear time. For R = 0 and under 1 -or ∞ -norm, Kalcsics et al. [7] constructed an algorithm solving the 1-median problem in O(n log d n) time. In that work the authors utilized a general framework of Cohen and Megiddo [4] developed for minimization of a broad class of convex functions over polytopes.
In this paper we develop an exact algorithm for the single facility location problem with disk-shaped demand regions in R d under rectilinear distance measures. The running time of our exact algorithm matches the running time of the algorithm by Kalcsics et al. [7], but solves the problem for arbitrary R ≥ 0. Then, we generalize the algorithm to generic polyhedral norms.

Problem definition
First, we derive the results for the 1 norm. Let ρ(x) = x 1 be the 1 -norm. Since a disk in R d under 1 -norm has a diamond shape, we refer to the disk-shaped demand regions as the diamond-shaped regions.
Given an input set X = {x 1 , x 2 , . . . , x n } ⊆ R d and R ≥ 0, the single facility location problem (SFLP) with diamond-shaped regions is to find a center c ∈ R d , which minimizes where (x) + = max{0, x}, x ∈ R. This problem can be reformulated as a linear program as follows. Let index j ∈ {1, 2, . . . , d} specify the jth coordinate of a point. Then, the linear program reads min c,y,z n i=1 z i where for each data point i ∈ {1, 2, . . . , n}, variable y i, j is the distance between c and x i projected on the jth coordinate, and z i is the distance between c and x i minus R if this value is non-negative, otherwise it equals 0. This linear program has nd + n + d variables and 2nd + n constraints. Vaidya [14] have shown that the time needed to solve such a linear program is O((nd) 2.5 L), where L is a parameter characterizing the constraint matrix. In this paper we introduce a simple combinatorial algorithm as an alternative to the solution of the above linear program. Moreover, our algorithm outperforms the classical linear programming algorithm.
The idea of the algorithm is purely geometrical. First, we observe, that the objective function f (c) is convex and piecewise linear. Second, we can define polynomially many polytopal regions, such that the objective function is linear in each region; see Fig. 1 for an illustration of a two-dimensional case with four data points.Third, the algorithm runs a binary search in each dimension over these polytopal regions and finds an optimal solution.

A two-dimensional 1 case
To simplify the presentation, we first describe the algorithm for the planar case, i.e. for d = 2. For one data point x, the basic division of the space into regions, where f is linear, is depicted in Fig. 2a. This basic drawing is centered at x. For the case of more than one data point, say for x 1 , x 2 , . . . , x n , we embed such basic drawings in the plane centering them at points x 1 , x 2 , . . . , x n ( Fig. 1), respectively. Similar to the graph drawing terminology, we call the union of the basic drawings an embedding. For Notice, division of the space into the regions where the objective function is linear is very well-known approach in Location Theory, see e.g. [10][11][12][13]. These regions are often referred to as linearity regions, domains of linearity, or cells, see the above references. For simplicity, in this paper we follow the latter notion of cells introduced in Nickel et al. [12]. The most important property of the problem is that within any cell of the embedding the function f is linear. Therefore the evaluation of the function in any point of a cell takes only O(n) time. In

Notation and definitions
Notice that under the 1 -norm for any x ∈ R d we have ρ(x) = max a∈{−1,1} d a T x, where {−1, 1} d is the set of d-dimensional vectors consisting entries 1 and −1. Therefore, we rewrite the objective function f as follows: Similar to the two-dimensional case, we define the notions of grid, cells and blocks in R d . We first define the set of directions A, which is the set of orthogonal vectors to the norm defining hyperplanes. Let being the set of positions/levels of the hyperplanes orthogonal to a k ∈ A. The intuition for the set B k is that the hyperplanes a T k x = b k j , x ∈ R d , b k j ∈ B k , are the delimiters of the linearity regions (cells) in the space, and these hyperplanes will be used in the construction of the high dimensional grid, blocks and cells, similarly to the twodimensional case in Sect. 3. Again, without loss of generality, assume b k

Definition 1
The hyperplanes x ∈ R d | a T k x = b k j , a k ∈ A, b k j ∈ B k form the grid. A cell is an inclusion maximal polyhedron in R d not crossed by a grid hyperplane. A set of indices { j 1 , j 2 , . . . , j | A| }, where j k ∈ {0, 1, . . . , |B k | + 1} for any k ∈ {1, 2, . . . , |A|}, completely specifies the cell: Leaving out l constraints in the definition of a cell, we arrive at the definition of a block of cells parameterized by l:

Observations
The following observations are necessary ingredients of the algorithm and its further analysis. First of all, notice that for any cell C and for any x ∈ X , either the cell is within the range R from x and then or there exists a vector a 0 ∈ {−1, 1} d such that In other words, every individual contribution of a data point x ∈ X to the objective function f is linear in any cell C. As the sum of linear functions is a linear function, f is linear in any cell C.
The second observation is that every individual contribution of a data point x ∈ X to the objective function f is convex as it is a maximum of two linear functions. Since the sum of convex functions is convex, the function f is convex.
The third observation is that for any cell C the minimum of f over C can be found in O(dn +

High dimensional recursive binary search algorithm
Now, we are ready to present the algorithm and to analyze its running time. Algorithm 1 works also for l ∞ and for more general polyhedral norms. We use the definition of a polyhedral norm as in Barvinok et al. [2]. Let P be a centrally symmetric bounded polyhedron in R d with the origin in its center. By symmetry, the number of facets is even. Let the number of facets be 2m. Then, there is a set H P = {h 1 , h 2 , . . . , h m } of points in R d such that P is the intersection of a collection of half-spaces defined by H P : Then, the polyhedral norm is defined by ρ(x) = max 1≤i≤m |h T i x|. For instance, for 1 in R 2 we can take H P = {(1, 1), (−1, 1)} and for ∞ in R 2 we can take For a polyhedral norm, the only adjustment we have to make in the Algorithm 1 is to generalize sets A and B in the definition of the grid. For this generalization, let Again, we assume all vectors in A are indexed: A = {a 1 , a 2 , . . . , a |A| }. For every vector a k ∈ A, we define Now, to adjust Algorithm 1 to a generic polyhedral norm we define the grid by the sets A and B k , a k ∈ A, according to the definitions (10) and (11). The rest of the algorithm remains intact.
Finally, let us estimate the run-time of Algorithm 1 adjusted to a polyhedral norm. The new set A determines the types of the grid hyperplanes. For each a k ∈ A the set B k determines all hyperplanes of type a k . Notice, |A| = O(m 2 ) as the size of A 2 is quadratic in the number of facets. By definition, for every a k ∈ A, the size of the set B k is at most 2n. Thus, the number of grid cells created in this construction is n O(m 2 ) . Using the adjusted number of grid cells and following literally the same arguments as in Theorem 1, we generalize the result as follows.
Theorem 2 Under a polyhedral norm determined by a polyhedron with 2m facets, the adjusted Algorithm 1 returns an optimal solution c * ∈ R d to the single facility location problem with disk-shaped demand regions in O (dn + d 4.5 m 2 ) log O(m 2 ) n time.
Since 2m ≥ d + 1 for any bounded polyhedron in R d , we have the following straightforward corollary: Corollary 1 Under a polyhedral norm determined by a polyhedron with 2m facets, the adjusted Algorithm 1 returns an optimal solution to the single facility location problem with disk-shaped demand regions in O (mn + m 6.5 ) log O(m 2 ) n time.