Unified Polynomial Dynamic Programming Algorithms for P-Center Variants in a 2D Pareto Front

With many efficient solutions for a multi-objective optimization problem, this paper aims to cluster the Pareto Front in a given number of clusters K and to detect isolated points. K-center problems and variants are investigated with a unified formulation considering the discrete and continuous versions, partial K-center problems, and their min-sum-K-radii variants. In dimension three (or upper), this induces NP-hard complexities. In the planar case, common optimality property is proven: non-nested optimal solutions exist. This induces a common dynamic programming algorithm running in polynomial time. Specific improvements hold for some variants, such as K-center problems and min-sum K-radii on a line. When applied to N points and allowing to uncover M<N points, K-center and min-sum-K-radii variants are, respectively, solvable in O(K(M+1)NlogN) and O(K(M+1)N2) time. Such complexity of results allows an efficient straightforward implementation. Parallel implementations can also be designed for a practical speed-up. Their application inside multi-objective heuristics is discussed to archive partial Pareto fronts, with a special interest in partial clustering variants.


Introduction
This paper is motivated by real-world applications of multi-objective optimization (MOO). Some optimization problems are driven by more than one objective function, with conflicting optimization directions. For example, one may minimize financial costs, while maximizing the robustness to uncertainties or minimizing the environmental impact [1,2]. In such cases, higher levels of robustness or sustainability are likely to induce financial over-costs. Pareto dominance, preferring one solution to another if it is better for all the objectives, is a weak dominance rule. With conflicting objectives, several non-dominated points in the objective space can be generated, defining efficient solutions, which are the best compromises. A Pareto front (PF) is the projection in the objective space of the efficient solutions [3]. MOO approaches may generate large sets of efficient solutions using Pareto dominance [3]. Summarizing the shape of a PF may be required for presentation to decision makers. In such a context, clustering problems are useful to support decision making to present a view of a PF in clusters, the density of points in the cluster, or to select the most central cluster points as representative points. Note than similar problems are of interest for population MOO heuristics such as evolutionary algorithms to archive representative points of a partial Pareto fronts, or in selecting diversified efficient solutions to process mutation or cross-over operators [4,5].
With N points in a PF, one wishes to define K N clusters while minimizing the measure of dissimilarity. The K-center problems, both in the discrete and continuous versions, define the cluster costs in this paper, covering the PF with K identical balls while minimizing the radius of the balls used. By definition, the ball centers belong to the PF for the discrete K-center version, whereas the continuous version is similar to geometric covering problems, without any constraint for the localization of centers. Furthermore, sum-radii or sum-diameter are min-sum clustering variants, where the covering balls are not necessarily identical. For each variant, one can also consider partial clustering variants, where a given percentage (or number) of points can be ignored in the covering constraints, which is useful when modelling outliers in the data.
The K-center problems are NP-hard in the general case, [6] but also for the specific case in R 2 using the Euclidean distance [7]. This implies that K-center problems in threedimensional (3D) PF are also NP-hard, with the planar case being equivalent to an affine 3D PF. We consider the case of two-dimensional (2D) PF in this paper, focusing on the polynomial complexity results. It as an application to bi-objective optimization, the 3D PF and upper dimensions are shown as perspectives after this work. Note that 2D PF are a generalization of one-dimensional (1D) cases, where polynomial complexity results are known [8,9]. A preliminary work proved that K-center clustering variants in a 2D PF are solvable in polynomial time using a Dynamic Programming (DP) algorithm [10]. This paper improves these algorithms for these variants, with an extension to min-sum clustering variants, partial clustering, and Chebyshev and Minkowski distances. The properties of the DP algorithms are discussed for efficient implementation, including parallelization. This paper is organized as follows. The considered problems are defined formally with unified notations in Section 2. In Section 3, related state-of-the-art elements are discussed. In Sections 4 and 5, intermediate results and specific complexity results for sub-cases are presented. In Section 6, a unified DP algorithm with a proven polynomial complexity is designed. In Section 7, specific improvements are presented. In Section 8, the implications and applications of the results of Sections 5-7 are discussed. In Section 9, our contributions are summarized, with a discussion of future research directions.

Problem Statement and Notation
In this paper, integer intervals are denoted as [[a, b] [1,N]] a set of N elements of R 2 , such that for all i = j, x i I x j defining the binary relations I, ≺ for all y = (y 1 , y 2 ), z = (z 1 , z 2 ) ∈ R 2 with y ≺ z ⇐⇒ y 1 < z 1 and y 2 > z 2 (1) y z ⇐⇒ y ≺ z or y = z (2) These hypotheses on E define 2D PF considering the minimization of two objectives [3,11]. Such configuration is illustrated by Figure 1. Without loss of generality, transforming objectives to maximize f into − f allows for the consideration of the minimization of two objectives. This assumption impacts the sense of the inequalities of I, ≺. A PF can also be seen as a Skyline operator [12]. A 2D PF can be extracted from any subset of R 2 using an output-sensitive algorithm [13], or using any MOO approach [3,14].
The results of this paper will be given using the Chebyshev and Minkowski distances, generically denoting d(y, z) the l 8 and l m norm-induced distance, respectively. For a given m > 0, the Minkowski distance is denoted d m , and given by the formula ∀y = (y 1 , y 2 ), z = (z 1 , z 2 ) ∈ R 2 , d m (y, z) = m |y 1 The case m = 2 corresponds to the Euclidean distance; it is a usual case for our application. The limit with m → 8 defines the Chebyshev distance, denoted d 8 and given by the formula ∀y = (y 1 , y 2 ), z = (z 1 , z 2 ) ∈ R 2 , d 8 (y, z) = max y 1 − z 1 , y 2 − z 2 (5) Once a distance d is defined, a dissimilarity among a subset of points E ⊂ E is defined using the radius of the minimal enclosing ball containing E . Numerically, this dissimilarity function, denoted as f C , can be written as A discrete variant considers enclosing balls with one of the given points as the center. Numerically, this dissimilarity function, denoted f D , can be written as For the sake of having unified notations for common results and proofs, we define γ ∈ {0, 1} to indicate which version of the dissimilarity function is considered. γ = 0 (respectively, 1) indicates that the continuous (respectively, discrete) version is used, f γ , thus denoting f 1 = f C (respectively, f 0 = f D ). Note that γ ∈ {0, 1} will be related to complexity results which motivated such a notation choice.
For each a subset of points E ⊂ E and integer K 1, we define Π K (E), as the set of all the possible partitions of E in K subsets. Continuous and discrete K-center are optimization problems with Π K (E) as a set of feasible solutions, covering E with K identical balls while minimizing the radius of the balls used min π∈Π K (E) max P∈π f γ (P) (8) The continuous and discrete K-center problems in the 2D PF are denoted K-γ-CP2DPF. Another covering variant, denoted min-sum-K-radii problems, covers the points with non-identical balls, while minimizing the sum of the radius of the balls. We consider the following extension of min-sum-K-radii problems, with α > 0 being a real number min π∈Π K (E) ∑ P∈π f γ (P) α (9) α = 1 corresponds to the standard min-sum-K-radii problem. α = 2 with the standard Euclidean distance is equivalent to the minimization of the area defined by the covering disks. For the sake of unifying notations for results and proofs, we define a generic operator ⊕ ∈ {+, max} to denote, respectively, sum-clustering and max-clustering. This defines the generic optimization problems Lastly, we consider a partial clustering extension of problems (10), similarly to the partial p-center [15]. The covering with balls mainly concerns the extreme points, which make the results highly dependent on outliers. One may consider that a certain number M < N of the points may be considered outliers, and that M points can be removed in the evaluation. This can be written as Problem (11) is denoted K-M-⊕-(α, γ)-BC2DPF. Sometimes, the partial covering is defined by a maximal percentage of outliers. In this case, if M is much smaller than N, we have M = Θ(N), which we have to keep in mind for the complexity results. K-center problems, K-γ-CP2DPF, are K-M-max-(α, γ)-BC2DPF problems for all α > 0; the value of α does not matter for max-clustering, defining the same optimal solutions as α = 1. The standard min-sum-k-radii problem, equivalent to the min-sum diameter problem, corresponds to k-0-+-(1, γ)-BC2DPF problems for discrete and continuous versions, k-M-+-(1, γ)-BC2DPF problems consider partial covering for min-sum-k-radii problems.

Related Works
This section describes works related to our contributions, presenting the state-of-theart for p-center problems and clustering points in a PF. For more detailed survey on the results for the p-center problems, we refer to [16].

Solving P-Center Problems and Complexity Results
Generally, the p-center problem consists of locating p facilities among possible locations and assigning n clients, called c 1 , c 2 , . . . , c n , to the facilities in order to minimize the maximum distance between a client and the facility to which it is allocated. The continuous p-center problem assumes that any place of location can be chosen, whereas the discrete p-center problem considers a subset o m potential sites, denoted f 1 , f 2 , . . . , f m , and distances d i,j for all i ∈ [ [1, n]] and j ∈ [ [1, m]]. Discrete p-center problems can be formulated with bipartite graphs, modeling that si unfeasibile for some assignments. In the discrete p-center problem defined in Section 2, points f 1 , f 2 , . . . , f m are exactly c 1 , c 2 , . . . , c n , and the distances are defined using a norm, so that triangle inequality holds for such variants.
P-center problems are NP-hard [6,17]. Furthermore, for all α < 2, any α-approximation for the discrete p-center problem with triangle inequality is NP-hard [18]. Two approximations were provided for the discrete p-center problem running in O(pn log n) time and in O(np) time, respecitvely [19,20]. The discrete p-center problem in R 2 with a Euclidean distance is also NP-hard [17]. Defining binary variables x i,j ∈ {0, 1} and y j ∈ {0, 1} with x i,j = 1 if and only if the customer i is assigned to the depot j, and y j = 1 if and only if location f j is chosen as a depot, the following Integer Linear Programming (ILP) formulation models the discrete p-center problem [21] min x,y,z z (12a) Constraints (12b) are implied by a standard linearization of the min-max original objective function. Constraint (12c) fixes the number of open facilities to p. Constraints (12d) assign each client to exactly one facility. Constraints (12e) are necessary to induce that any considered assignment x i,j = 1 implies that facility j is open with y j = 1. Tighter ILP formulations than (12) were proposed, with efficient exact algorithms relying on the IP models [22,23]. Exponential exact algorithms were also designed for the continuous p-center problem [24,25]. An n O( √ p) -time algorithm was provided for the continuous Euclidean p-center problem in the plane [26]. An n O(p 1−1/d ) -time algorithm is available for the continuous p-center problem in R d under Euclidean and L 8 -metric [27].
Specific cases of p-center problems are solvable in polynomial time. The continuous 1-center problem is exactly the minimum covering ball problem, which has a linear complexity in R 2 . Indeed, a "prune and search" algorithm finds the optimum bounding sphere and runs in linear time if the dimension is fixed as a constant [28]. In dimension d, its complexity is in O((d + 1)(d + 1)!n) time, which is impractical for high-dimensional applications [28]. The discrete 1-center problem is solved in time O(n log n), using furthestneighbor Voronoi diagrams [29]. The continuous and planar 2-center problem is solved in randomized expected O(n log 2 n) time [30,31]. The discrete and planar 2-center problem is solvable in O(n 4/3 log 5 n) time [32].
1D p-center problems, or those with equivalent points that are located in a line, have specific complexity results with polynomial DP algorithms. The discrete 1D k-center problem is solvable in O(n) time [33]. The continuous and planar k-centers on a line, finding k disks with centers on a given line l, are solvable in polynomial time, in O(n 2 log 2 n) time in the first algorithm by [29], and in O(nk log n) time and O(n) space in the improved version provided by [34]. An intensively studied extension of the 1D sub-cases is the p-center in a tree structure. The continuous p-center problem is solvable in O(n log 3 n) time in a tree structure [7]. The discrete p-center problem is solvable in O(n log n) time in a tree structure [35].
Rectilinear p-center problems, using the Chebyshev distances, were less studied. Such distance is useful for complexity results; however, it has fewer applications than Euclidean or Minkowski norms. For the planar and rectangular 1-center and 2-center problems, O(n) algorithms are available for the 1-center problem, and such 3-center problems can be solved in O(n log n) time [36]. In a general dimension d, continuous and discrete versions of rectangular p-center problems are solvable in O(n) and O(n log d−2 n log log n + n log n) running time, respectively. Specific complexity results for rectangular 2-center problems are also available [37].

Solving Variants of P-Center Problems and Complexity Results
Variants of p-center problems were studied less intensively than the standard p-center problems. The partial variants were introduced in 1999 by [15], whereas a preliminary work in 1981 considered a partial weighted one-center variant and a DP algorithm to solve it running in O(n 2 log n) time [38]. The partial discrete p-center can formulated as an ILP starting from the formulation provided by [21] as written in (12). Indeed, considering that n 0 points can be uncovered, constraints (12.4) become inequalities ∑ m j=1 x i,j 1 for all i, j and the maximal number of unassigned points is set to n 0 , adding one constraint ∑ n j=i ∑ m j=1 x i,j n − n 0 . Similarly, the sum-clustering variants K-M-+-(α, γ)-BC2DPF can be written as the following ILP min z,r 0 In this ILP formulation, binary variables z n,n ∈ {0, 1} are defined such that z n,n = 1 if and only if the points x n and x n are assigned in the same cluster, with x n being the discrete center. Continuous variables r n 0 denote the powered radius of the ball centered in x n , if x n is chosen as a center, and r n = 0 otherwise. Constraint (13b) is a standard linearization of the non-linear objective function. z n,n indicates that if point x n is chosen as the center, then this implies with (13c) that K such variables are nonzero, and with (13f) that a nonzero variable z n,n implies that the corresponding z n ,n is not null. (13d) and (13e) allow the extension with partial variants, as discussed before.
Min-sum radii or diameter problems were rarely studied. However, such objective functions are useful for meta-heuristics to break some "plateau" effects [39]. Min-sum diameter clustering is NP-hard in the general case and polynomial within a tree structure [40]. The NP-hardness is also proven, even in metrics induced by weighted planar graphs [41]. Approximation algorithms were studied for min-sum diameter clustering. A logarithmic approximation with a constant factor blowup in the number of clusters was provided by [42]. In the planar case with Euclidean distances, a polynomial time approximation scheme was designed [43].

Clustering/Selecting Points in Pareto Frontiers
Here, we summarize the results related to the selection or the clustering of points in PF, with applications for MOO algorithms. Polynomial complexity resulting in the use of 2D PF structures is an interesting property; clustering problems have a NP-hard complexity in general [17,44,45].
To the best of our knowledge, no specific work focused on PF sub-cases of k-center problems and variants before our preliminary work [10]. A Distance-Based Representative Skyline with similarities to the discrete p-center problem in a 2D PF may not be fully available in the Skyline application, which makes a significant difference [46,47]. The preliminary results proved that K-γ-CP2DPF is solvable in O(KN log γ N) time using O(N) additional memory space [10]. Partial extensions and min-sum-k-radii variants were not considered for 2D PF. We note that the 2D PF case is an extension of the 1D case, with 1D cases being equivalent to the cases of an affine 2D PF. In the study of complexity results, a tree structure is usually a more studied extension of 1D cases. The discrete k-center problem on a tree structure, and thus the 1D sub-case, is solvable in O(N) time [33]. 3F PF cases are NP-complete, as already mentioned in the introduction, this being a consequence of the NP-hardness of the general planar case.
Maximization of the quality of discrete representations of Pareto sets was studied with the hypervolume measure in the Hypervolume Subset Selection (HSS) problem [48,49]. The HSS problem is known to be NP-hard in dimension 3 (and greater dimensions) [50]. HSS is solvable with an exact algorithm in N O( √ K) and a polynomial-time approximation scheme for any constant dimension d [50]. The 2D case is solvable in polynomial time with a DP algorithm with a complexity in O(KN 2 ) time and O(KN) space [49]. The time complexity of the DP algorithm was improved in O(KN + N log N) by [51], and in O(K(N − K) + N log N) by [52].
The selection of points in a 2D PF, maximizing the diversity, can also be formulated using p-dispersion problems. Max-Min and Max-Sum p-dispersion problems are NPhard problems [53,54]. Max-Min and Max-Sum p-dispersion problems are still NP-hard problems when distances fulfill the triangle inequality [53,54]. The planar (2D) Max-Min p-dispersion problem is also NP-hard [9]. The one-dimensional (1D) cases of Max-Min and Max-Sum p-dispersion problems are solvable in polynomial time, with a similar DP algorithm running in O(max{pN, N log N}) time [8,9]. Max-Min p-dispersion was proven to be solvable in polynomial time, with a DP algorithm running in O(pN log N) time and O(N) space [55]. Other variants of p-dispersion problems were also proven to be solvable in polynomial time using DP algorithms [55].
Similar results exist for k-means, k-medoid and k-median clustering. K-means is NP-hard for 2D cases, and thus for 3D PF [44]. K-median and K-medoid problems are known to be NP hard in dimension 2, since [17], where the specific case of 2D PF was proven to be solvable in O(N 3 ) time with DP algorithms [11,56]. The restriction of k-means to 2D PF would be also solvable in O(N 3 ) time with a DP algorithm if a conjecture was proven [57]. We note that an affine 2D PF is a line in R 2 , where clustering is equivalent to 1D cases. 1D k-means were proven to be solvable in polynomial time with a DP algorithm in O(KN 2 ) time and O(KN) space. This complexity was improved for a DP algorithm in O(KN) time and O(N) space [58]. This is thus the complexity of K-means in an affine 2D PF.

Lemma 1.
is an order relation, and ≺ is a transitive relation ∀x, y, z ∈ R 2 , x ≺ y and y ≺ z =⇒ x ≺ z Proposition 1 implies an order among the points of E, for a re-indexation in O(N log N) time Proof.
We index E such that the first coordinate is increasing. This sorting procedure runs The re-indexation also implies monotonic relations among distances of the 2D PF Lemma 2. We suppose that E is re-indexed as in Proposition 1. Letting d be a Minkowski, Euclidean or Chebyshev distance, we obtain the following monotonicity relations Proof. We first note that the equality cases are trivial, so we can suppose that i 1 < i 2 < i 3 in the following proof. We prove the propriety (17); the proof of (18) is analogous.
. Hence, the result is proven for Euclidean, Minkowski and Chebyshev distances.

Lemmas Related to Cluster Costs
This section provides the relations needed to compute or compare cluster costs. Firstly, one notes that the computation of cluster costs is easy in a 2D PF in the continuous clustering case.

Lemma 3.
Let P ⊂ E, such that card(P) 1. Let i (resp i ) be the minimal (respective maximal) index of points of P with the indexation of Proposition 1. Then, f C (P) can be computed with To prove the Lemma 3, we use the Lemmas 4 and 5.

Lemma 4.
Let P ⊂ E, such that card(P) 1. Let i (resp i ) the minimal (resp maximal) index of points of P with the indexation of Proposition 1. We denote with O =

Proof of Lemma 4:
with the equality being trivial as points O, x i , x i are on a line and d is a distance. Let x ∈ P. We calculate the distances using a new system of coordinates, translating the original coordinates such that O, is a new origin (which is compatible with the definition of Pareto optimality). x i and x i have coordinates (−a, b) and (a, −b) in the new coordinate system, with a, b > 0 and a m + b m = r m if a Minkowski distance is used, otherwise it is max(a, b) = r for the Chebyshev distance. We use (a , b ) to denote the coordinates of x. x i ≺ x ≺ x i implies that −a a a and −b b b, i.e., |a | a and |b | b, which implies d(x, O) r, using Minkowski or Chebyshev distances.

Lemma 5.
Let P ⊂ E such that card(P) 1. Let i (respective i ) be the minimal (respective maximal) index of points of P with the indexation of Proposition 1. We denote, using O = x i +x i 2 , the midpoint of x i , x i . Then, using a Minkowski or Chebyshev distance d, we have for all y ∈ R 2 :

Proof of Lemma 3:
We first note that f C (P) = min y∈R 2 max x∈P d(x, y) max x∈P d(x, O), Reciprocally, for all y ∈ R 2 , r max(d(x i , y), d(x i , y)) using Lemma 5, and thus r max x∈P d(x, y). This implies that r min y∈R 2 max x∈P d(x, y) = f C (P). Lemma 6. Let P ⊂ E such that card(P) 3. Let i (respective i ) the minimal (respective maximal) index of points of P.
Lastly, we notice that extreme points are not optimal centers. Indeed, Proof. Using the order of Proposition 1, let i (respectively, i ) the minimal index of points of P (respectively, P ) and let j (respectively, j ) the maximal indexes of points of P (respectively, P ). f C (P) f C (P ) is trivial using Lemmas 2 and 3. To prove f D (P) f D (P ), we use i i j j, and Lemmas 2 and 6 Let P ⊂ E, such that card(P) 1. Let i (respectively, i ) the minimal (respectively, maximal) index of points of P. For all P ⊂ P, such that

Optimality of Non-Nested Clustering
In this section, we prove that non-nested clustering property, the extension of interval clustering from 1D to 2D PF, allows the computation of optimal solutions, which will be a key element for a DP algorithm. For (partial) p-center problems, i.e., K-M-max-(α, γ)-BC2DPF, optimal solutions may exist without fulfilling the non-nested property, whereas for K-M-+-(α, 0)-BC2DPF problems, the nested property is a necessary condition to obtaining an optimal solution.
Proof. Let π ∈ Π K (E) represent an optimal solution of 1-M-⊕-(α, γ)-BC2DPF, let OPT be the optimal cost, and C ⊂ E with |C| N − M and f γ (C) = OPT. Let i (respectively, i ) be the minimal (respectively maximal) index of C using order of Proposition 1. C ⊂ C i,i , so Lemma 8 applies and Proof. We prove the results by the induction on K ∈ N * . For K = 1, Lemma 9 gives the initialization.
Let us suppose that K > 1 and the Induction Hypothesis (IH) that Proposition 2 is true for K-M-⊕-(α, γ)-BC2DPF. Let π ∈ Π K (E) be an optimal solution of K-M-⊕-(α, γ)-BC2DPF; let OPT be the optimal cost. Let X ⊂ E be the subset of the non-selected points, |X| M, and C 1 , . . . , C K be the K subsets defining the costs, so that X, We consider the subsets . . , C K−1 , C K is an optimal solution of K-M-⊕-(α, γ)-BC2DPF in E using only clusters C i,i . Hence, the result is proven by induction.
Proof. Starting with an optimal solution of K-M-+-(α, 0)-BC2DPF, let OPT be the optimal cost, and let X ⊂ E be the subset of the non-selected points, |X| M, and C 1 , . . . , C K , the K subsets defining the costs, so that X, C 1 , . . . , C K is a partition of E. Removing random M − |X| points in C 1 , . . . , C K , we have clusters C 1 , . . . , and thus the clusters C 1 , . . . , C K and outliers X = E \ ∪ k C k define and provide the optimal solution of K-M-⊕-(α, 0)-BC2DPF with exactly M outliers.
Reciprocally, one may investigate if the conditions of optimality in Propositions 2 and 3 are necessary. The conditions are not necessary in general. For instance, with E = {(3, 1); (2, 2); (1, 3)}, K = M = 1 and the discrete function F D , ie γ = 1, the selection of each pair of points defines an optimal solution, with the same cost as the selection of the three points, which do not fulfill the property of Proposition 3. Having an optimal solution with the two extreme points also does not fulfill the property of Proposition 2. The optimality conditions are necessary in the case of sum-clustering, using the continuous measure of the enclosing disk.

Proposition 4.
Let an optimal solution of K-M-+-(α, 0)-BC2DPF be defined with X ⊂ E as the subset of outliers, with |X| M, and C 1 , . . . , C K as the K subsets defining the optimal cost. We therefore have Proof. Starting with an optimal solution of K-M-+-(α, 0)-BC2DPF, let OPT be the optimal cost, and let X ⊂ E be the subset of the non-selected points, |X| M, and C 1 , . . . , C K be the K subsets defining the costs, so that X, C 1 , . . . , C K is a partition of E. We prove (i) and (ii) ad absurdum.
If |X| < M, one may remove one extreme point of the cluster C 1 , defining C 1 . With Lemmas 2 and 3, we have This is in contraction with the optimality of C 1 , . . . , C K , C 1 , C 2 . . . , C K , defining a strictly better solution for we have nested clusters C l and C k . We suppose that i k < i l (otherwise, reasoning is symmetrical). We define a better solution than the optimal one with C

Computation of Cluster Costs
Using Proposition 2, only cluster costs C i,i are computed. This section allows the efficient computation of such cluster costs. Once points are sorted using Proposition 1, cluster costs f C (C i,i ) can be computed in O(1) using Lemma 3. This makes a time complexity in O(N 2 ) to compute all the cluster costs f C (C i,i ) for 1 i i N. Equation (19) ensures that cluster costs f D (C i,i ) can be computed in O(i − i) for all i < i . Actually, Algorithm 1 and Proposition 5 allow for computations in O(log(i − i)) once points are sorted following Proposition 1, with a dichotomic and logarithmic search.
Let i < i . Proposition 2, applied to i and any j, j + 1 with j i and j < i , ensures that g is decreasing. Similarly, Proposition 2, applied to i and any j, j + 1, ensures that h is increasing.
A is a non-empty and bounded subset of N, so that A has a maximum. We note that l = max A.h i,i (i ) = 0 and g i, l + 1 / ∈ A and g i,i (l + 1) > h i,i (l + 1) have to be coherent with the fact that l = max A.
Lastly, the minimum of f can be reached in l or in l + 1, depending on the sign of , there are two minimums l, l + 1. Otherwise, there is a unique minimum l 0 ∈ {l, l + 1}, f i,i , which decreases before increasing.
, and thus i j * , using Lemma 10. Similarly, we always obtain Remark 1. Algorithm 1 improves the previously proposed binary search algorithm [10]. If it has the same logarithmic complexity, this leads to two times fewer calls of the distance function. Indeed, in the previous version, the dichotomic algorithm is computed at each iteration f i,i (i ) and f i,i (i + 1) to determine if i is in the increasing or decreasing phase of f i,i . In Algorithm 1, the computations that are provided for each iteration are equivalent to the evaluation of only f i,i (i ), Proof. We prove (i); we suppose that i < N and we prove that, for all c < c f i,i +1 (c) f i,i +1 (c ), so that either c is an argmin of the minimization, and the superior minimum to Proof. The validity of Algorithm 2 is based on Lemmas 10 and 11: once a discrete center c is known for a f D (C j,j ) α , we can find a center c of f D (C j,j +1 ) α with c c, and Lemma 10 gives the stopping criterion to prove a discrete center. Let us prove the time complexity; the space complexity is obviously within O(N) memory space. In Algorithm 2, each computation f j ,j (curCtr) is in O(1) time; we have to count the number of calls for this function. In each loop in j , one computation is used for the initialization; the total number of calls for this initialization is N − j N. Then, denoting, with c N N, the center found for C j,N , we note that the number of loops is c N − j N. Lastly, there are less that 2N computations calls f j ,j (curCtr); Algorithm 2 runs in O(N) time. Proof. The proof is analogous with Proposition 6, applied to Algorithm 2.

Particular Sub-Cases
Some particular sub-cases have specific complexity results, which are presented in this section.  (N log N). In the continuous case, i.e., γ = 0, one requires only the M minimal and maximal points with the total order ≺ to compute the cluster costs using. If M < log N, one may use one traversal of E, storing the current m minimal and extreme points, which has a complexity in O(MN). Choosing among the two possible algorithms, the time complexity is in O (N min(M, log N)).

Sub-Cases with K = 2
Specific cases with K = 2 define two clusters, and one separation as defined in Proposition 2. For these cases, specific complexity results are provided, enumerating all the possible separations. Proof. After the re-indexation phase running in O(log N) time), Proposition 2 ensures that there is an optimal solution for 2-M-⊕-(α, γ)-BC2DP, removing the m 1 0 first indexes, the m 3 0 last indexes, and m 2 0 points between the two selected clusters, with m 1 + m 2 + m 3 M. Using Proposition 3, there is an optimal solution, exactly defining the M outliers, so that we can consider that m 1 + m 2 + m 3 = M. Denoting i as the last index of the first cluster, the first selected cluster is C 1+m 1 ,i ; the second one is C i+m 2 +1,N−M+m 1 +m 2 . We have i m 1 + 1 and i + m 2 + 1 N − M + m 1 + m 2 i.e., i N − M + m 1 . We denote, with X, the following feasible i, m 1 , m 2 Computing an optimal solution for 2-M-⊕-(α, γ)-BC2DP brings the following enumeration OPT = min In In the discrete case (i.e., γ = 1), we use O((M + 1) 2 ) computations to enumerate the possible m 1 , m 2 , and O(N) computations to enumerate the possible i once m 1 , m 2 are fixed. This uses O(1) additional memory space, and the total time complexity is O(N log N((M + 1) 2 ). To decrease the time complexity, one can use two vectors of size N to store a given m 1 , m 2 , for which the cluster costs f γ (C 1+m 1 ,i ) α and f γ (C i+m 2 +1,N−M+m 1 +m 2 ) α are given by Algorithms 2 and 3, so that the total time complexity remains in O(N((M + 1) 2 + log N)) with an O(N) additional memory space. These two variants, using O(1) or O(N) additional memory space, induce the time complexity announced in Proposition 11.

Continuous Min-Sum K-Radii on A Line
To the best of our knowledge, the 1D continuous min-sum k-radii and the min-sum diameter problems were not previously studied. The specific properties hold, as proven in Lemma 12. This allows a time complexity of O(N log N).

Lemma 12.
Let E = {x 1 , . . . , x N } be N points in a line of R 2 , indexed such that for all i < j, x i ≺ x j . The min-sum k-radii in a line, K-0-+-(1, 0)-BC2DPF, is equivalent to selecting the K − 1 highest values of the distance among consecutive points, with the extremity of such segments defining the extreme points of the disks.

Proof. Let a feasible and non-nested solution of
Using the alignment property, we can obtain Reciprocally, this is equivalent to considering K-0-+-(1, 0)-BC2DPF or the maximization of the sum of K − 1 sa a different distance among consecutive points. The latter problem is just computing the K − 1 highest distances among consecutive points.

Algorithm 4: Continuous min-sum K-radii on a line
return OPT the optimal cost and the partition of selected clusters P

Unified DP Algorithm and Complexity Results
Proposition 2 allows the design of a common DP algorithm for p-center problems and variants, and to prove polynomial complexities. The key element is to design Bellman equations. Proof.

Specific Improvements
This section investigates how the complexity results of Theorem 2 may be improved, and how to speed up Algorithm 5, from a theoretical and a practical viewpoint.

Improving Time Complexity for Standard and Partial P-Center Problems
In Algorithm 5, the bottleneck for complexity are the computations Proof. We yfirst note that the case k = 1 is implied by the Lemma 7, so that we can suppose in the following, that . Let π ∈ Π K (E) be an optimal solution of k-m-⊕-(α, γ)-BC2DPF among the points indexed in [ [1, j]]; its cost is O j,k,m . Let X ⊂ E, the subset of the non-selected points, |X| M, and C 1 , . . . , C k with the k subsets defining the costs, so that X, C 1 , . . . , C k is a partition of E and k k =1 f γ (C k ) α = O j,k,m . If x j ∈ X, then O j,k,m = O j−1,k,m−1 O j−1,k,m using Lemma 13, which is the result. We suppose to end the proof that x j / ∈ X and re-index the clusters such that x j ∈ C k . We consider the clusters C 1 , . . . , C k = C 1 , . . . , C k−1 , C k − x k . With X, a partition of (x l ) l∈[[1,j−1]] is defined, with, at most, M outliers, so that it defines a feasible solution of the optimization problem, defining O j−1,k,m as a cost OPT O j−1,k,m . Using Lemma 7, OPT O j,k,m , so that O j−1,k,m−1 O j−1,k,m .

Proof.
Similarly to the proof of Lemma 10, the following applications are monotone:

Computing min
in the proof of Theorem 1 for p-center problem and variants, the complexity results are updated for these sub-problems.

Remark 2.
For the standard discrete p-center, Theorem 2 improves the time complexity given in the preliminary paper [10], from O(pN log 2 N) to O(pN log N). Another improvement was given by Algorithm 1; the former computation of cluster costs has the same asymptotic complexity but requires two times more computations. tTis proportional factor is non negligible in practice.

Improving Space Complexity for Standard P-Center Problems
For standard p-center problems, Algorithm 5 has a complexity in memory space in O(KN), the size of the DP matrix. This section proves it is possible to reduce the space complexity into an O(N) memory space.
One can compute the DP matrix for k-centers "line-by-line", with k increasing. This does not change the validity of the algorithm, with each computation using values that were previously computed to the optimal values. Two main differences occur compared to Algorithm 5. On one hand, the k + 1-center values use only k-center computations, and the computations with k < k can be deleted once all the required k-center values are computed when having to compute only the K-center values, especially the optimal cost. On the other hand, the computations of cluster costs are not factorized, as in Algorithm 5; this does not make any difference in the continuous version, where Lemma 3 can to recompute cluster costs in O(1) time when needed, whereas recomputing each cost induces the computations running in O(log N) for the discrete version with Algorithm 1.
The search order of operations slightly degrades the time complexity for the discrete variant, without inducing a change in the continuous variant. This allows only for computations of the optimal value; another difficulty is that the backtracking operations, as written in Algorithm 5, require storage of the whole stored values of the whole matrix. The issue is obtaining alternative backtracking algorithms that allow the computation of an optimal solution of the standard p-center problems using only the optimal value provided by the DP iterations, and with a complexity of, at most, O(KN log γ N) time and O(N) memory space. Algorithms 7 and 8 have such properties. -N points of a 2D PF, E = {z 1 , . . . , z N }, sorted such that for all i < j, z i ≺ z j ; -K ∈ N the number of clusters; -OPT, the optimal cost of K-γ-CP2DPF; output: P an optimal partition of K-γ-CP2DPF. -N points of a 2D PF, E = {z 1 , . . . , z N }, sorted such that for all i < j, z i ≺ z j ; -K ∈ N the number of clusters; -OPT, the optimal cost of K-γ-CP2DPF; output: P an optimal partition of K-γ-CP2DPF. Proof. This lemma is proven a decreasing induction on k, starting from k = K − 1. The case k = K − 1 is furnished by the first step of Algorithm 4, and j ∈ [ [1, N]] −→ f γ (C j,N ) decreaswa with Lemma 7. WIth a given k, i k i k , i k−1 i k−1 is implied by Lemma 2 and d(z i k , z i k−1 −1 ) > OPT.
Algorithm 8 is similar to Algorithm 7, with iterations increasing the indexes of the points of E. The validity is similarly proven, and this provides the upper bounds for the indexes of any optimal solution of K-center problems.
Analyzing the complexity, Algorithm 7 calls for a maximum of (K + N) 2N times the clustering cost function, without requiring stored elements; the complexity is in O(N log γ N) time.  (KN log 2 N). Depending on the value of K, it may be preferable, with stronger constraints in memory space, to have this second version.

Improving Space Complexity for Partial P-Center Problems?
This section tries to generalize the previous results for the partial K-center problems, i.e., K-M-max-(α, γ)-BC2DPF with M > 0. The key element is to obtain a backtracking algorithm that does not use the DP matrix. Algorithm 10 extends Algorithm 7 by considering all the possible cardinals of outliers between clusters k and k + 1 for k ∈ [[0, K − 1]] and the outliers after the last cluster. A feasible solution of the optimal cost should be feasible by iterating Algorithm 7 for at least one of these sub-cases.  We suppose there exist j 0 ∈ [ [1, i]], such that f γ (C j 0 ,i ) α β. Then, each optimal index j * , such Proof. With f γ (C j 0 ,i ) α β, Lemma 7 implies that for all j < j 0 , f γ (C j,i ) α > f γ (C j 0 ,i ) α β. Using O i ,k ,m 0 for all i , k implies that for all j < j 0 , f γ (C j 0 ,i ) α + O j 0 −1,k−1,m > β, and the optimal index gives O i,k,m = min j∈[[k+m,i]] O j−1,k−1,m + f γ (C j,i ) α , which is superior to j 0 .
Proposition 16 can be applied to compute each value of the DP matrix using fewer computations than the naive enumeration. In the enumeration, β is updated to the best current value of O j−1,k−1,m + f γ (C j,i ) α . The index would be enumerated in a decreasing way, starting from j = i until an index is found, such that f γ (C j 0 ,i ) α β, and no more enumeration is required with Proposition 16, ensuring that the partial enumeration is sufficient to find the wished-for minimal value. This is a practical improvement, but we do not furnish proof of complexity improvements, as it is likely that this would not change the worst case complexity.

Importance of the 2D PF Hypothesis, Summarizing Complexity Results
Planar p-center problems were not studied previously in the PF case. The 2D PF hypothesis is crucial for the complexity results and the efficiency of the solving algorithms. Table 1 compares the available complexity results for 1D and 2D cases of some k-center variants.
The complexity for 2D PF cases is very similar to the 1D cases; the 2D PF extension does not induce major difficulties in terms of complexity results. 2D PF cases may induce significant differences compared to the general 2D cases. The p-center problems are NP-hard in a planar Euclidean space [17], since adding the PF hypothesis leads to the polynomial complexity of Theorem 1, which allows for an efficient, straightforward implementation of the algorithm. Two properties of 2D PF were crucial for these results: The 1D structure implied by Proposition 1, which allows an extension of DP algorithms [58,59], and Lemmas 3 and 6, which allow quick computations of cluster costs. Note that rectangular p-center problems have a better complexity using general planar results than using our Theorems 2 and 3. Our algorithms only use common properties for Chebyshev and Minkowski distances, whereas significant improvements are provided using specificities of Chebyshev distance. Note that our complexity results are given considering the complexity of the initial re-indexation with Proposition 1. This O(N log N) phase may be the bottleneck for the final complexity. Some papers mention results which consider that the data are already in the memory (avoiding an O(N) traversal for input data) and already sorted. In our applications, MOO methods such as epsilon-constraint provide already sorted points [3]. Using this means of calculating the complexity, our algorithms for continuous and discrete 2-center problems in a 2D PF would have, respectively, a complexity in O(log N) and O(log 2 N) time. A notable advantage of the specialized algorithm in a 2D PF instead of the general cases in 2D is the simple and easy to implement algorithms.

Equivalent Optimal Solutions for P-Center Problems
Lemmas 16 and 17 emphasize that many optimal solutions may exist; the lower and upper bounds may define a very large funnel. We also note that many optimal solutions can be nested, i.e., non-verifying the Proposition 2. For real-world applicationa, having well-balanced clusters is more natural, and often wished for. Algorithms 7 and 8 provide the most unbalanced solutions. One may balance the sizes of covering balls, or the number of points in the clusters. Both types of solutions may be given using simple and fast post-processing. For example, one may proceed with a steepest descent local search using two-center problem types for consecutive clusters in the current solution. For balancing the size of clusters, iterating two-center computations induces marginal computations in O(log 1+γ N) time for each iteration with Algorithm 6. Such complexity occurs once the points are re-indexed using Proposition 1; one such computation in O(N log N) allows for many neighborhood computations running in O(log 1+γ N) time, and the sorting time is amortized.

Towards a Parallel Implementation
Complexity issues are raised to speed-up the convergence of the algorithms in practice. An additional way to speed up the algorithms in practice is to consider implementation issues, especially parallel implementation properties in multi-or many-core environments. In Algorithm 5, the values of the DP matrix O i,k,m for a given i ∈ [[1; N]] requires only to compute the values O j,k,m for all j < i . Independent computations can thus be operated at the iteration i of Algorithm 5, once the cluster costs f γ (C i ,i ) α for all i ∈ [[1; i]] have been computed, which is not the most time-consuming part when using Algorithms 2 and 3. This is a very useful property for a parallel implementation, requiring only N − 1 synchronizations to process O(KN 2 (1 + M)) operations. Hence, a parallel implementation of Algorithm 5 is straightforward in a shared memory parallelization, using OpenMP for instance in C/C++, or higher-level programming languages such as Python, Julia or Chapel [60]. One may also consider an intensive parallelization in a many-core environment, such as General Purpose Graphical Processing Units (GPGPU). A difficulty when using this may be the large memory size that is required in Algorithm 5.
Section 7 variants, which construct the DP matrix faster, both for k-center and min-sum k-radii problems, are not compatible with an efficient GPGPU parallelization, and one would prefer the naive and fixed-size enumeration of Algorithm 5, even with its worse time complexity for the sequential algorithm. Comparing the sequential algorithm to the GPGPU parallelization, having many independent parallelized computations allows a huge proportional factor with GPGPU, which can compensate the worst asymptotic complexity for reasonable sized instances. Shared memory parallelization, such as OpenMP, is compatible with the improvements provided in Section 7. Contrary to Algorithm 5, Algorithms 9 and 11 compute the DP matrix with index k increasing, with O(N) independent computation induced at each iteration. With such algorithms, there are only K − 2 synchronizations required, instead of N − 1 for Algorithm 5, which is a better property for parallelization. The O(N) memory versions are also useful for GPGPU parallelization, where memory space is more constrained than when storing a DP matrix on the RAM.
Previously, the parallelization of the DP matrix construction was discussed, as this is the bottleneck in time complexity. The initial sorting algorithm can also be parallelized on GPGPU if needed; the sorting time is negligible in most cases. The backtracking algorithm is sequential to obtain clusters, but with a low complexity in general, so that a parallelization of this phase is not crucial. We note that there is only one case where the backtracking Algorithm has the same complexity as the construction of the DP matrix: the DP variant in O(N) memory space proposed in Algorithm 11 with Algorithm 10 as a specific backtrack. In this specific case, the O(K) tests with different positions of the chosen outlier are independent, which allows a specific parallelization for Algorithm 10.

Applications to Bi-Objective Meta-Heuristics
The initial motivation of this work was to support decision makers when an MOO approach without preference furnishes a large set of non-dominated solutions. In this application, the value of K is small, allowing for human analyses to offer some preferences. In this paper, the optimality is not required in the further developments. Our work can also be applied to a partial PF furnished by population meta-heuristics [5]. A posteriori, the complexity allows for the use of Algorithms 5, 9 and 11 inside MOO meta-heuristics. Archiving PF is a common issue of population meta-heuristics, facing multi-objective optimization problems [4,5]. A key issue is obtaining diversified points of the PF in the archive, to compute diversified solutions along the current PF. Algorithms 5, 9 and 11 can be used to address this issue, embedded in MOO approaches, similarly to [49]. Archiving diversified solutions of Pareto sets has application for the diversification of genetic algorithms, to select diversified solutions for cross-over and mutation phases [61,62], but also for swarm particle optimization heuristics [63]. In these applications, clustering has to run quickly. The complexity results and the parallelization properties are useful in such applicationas.
For application to MOO meta-heuristics like evolutionary algorithms, the partial versions are particularly useful. Indeed, partial versions may detect outliers that are isolated from the other points. For such points, it is natural to process intensification operators to look for efficient solutions in the neighborhood, which will make the former outlier less isolated. Such a process is interesting for obtaining a better balanced distribution of the points along the PF, which is a crucial point when dealing with MOO meta-heuristics.

How to Choose K, M?
A crucial point in clustering applications the selection of an appropriate value of K. A too-small value of K can miss that instances which are well-captured with K + 1 representative clusters. Real-world applications seek the best compromise between the minimization of K, and the minimization of the dissimilarity among the clusters. Similarly, with [11], the properties of DP can be used to achieve this goal. With the DP Algorithm 9, many couples {(k, O N,k )} k are computed, using the optimal k-center values with k clusters. Having defined a maximal value K , the complexity for computing these points is seen in O(NK log 1+γ N). When searching for good values of k, the elbow technique, may be applied. Backtracking operations may be used for many solutions without changing the complexity. Rhe same ideas are applicable along the M index. In the previoulsy described context of MOO meta-heuristics, the sensitivity with the M parameter is more important than the sensitivity for the parameter K, where the number of archived points is known and fixed regarding other considerations, such as the allowed size of the population.

Conclusions and Perspectives
This paper examined the properties of p-center problems and variants in the special case of a discrete set of non-dominated points in a 2D space, using Euclidean, Minkowski or Chebyshev distances. A common characterization of optimal clusters is proven for the discrete and continuous variants of the p-center problems and variants. Thie can solve these problems to optimality with a unified DP algorithm of a polynomial complexity. Some complexity results for the 2D PF case improve the general ones in 2D. The presented algorithms are useful for MOO approaches. The complexity results, in O(KN log N) time for the standard K-center problems, and in O(KN 2 ) time for the standard min-sum kradii problems, are useful for application with a large PF. When applied to N points and able to ncover M < N points, partial K-center and min-sum-K-radii variants are, respectively, solvable in O(K(M + 1)N log N) and O(K(M + 1)N 2 ) time. Furthermore, the DP algorithms have interesting properties for efficient parallel implementation in a shared memory environment, such as OpenMP or using GPGPU. This allows their application for a very large PF with short solving times. For an application for MOO meta-heuristics such as evolutionary algorithms, the partial versions are useful for the detection of outliers where intensification phases around these isolated solutions may be processed in order to obtain a better distribution of the points along the PF.
Future perspectives include the extension of these results to other clustering algorithms. The weighted versions of p-center variants were not studied in this paper, which was motivated by MOO perspectives, and future perspectives shall consider extending our algorithms to weighted variants. Regarding MOO applications, extending the results to dimension 3 is a subject of interest for MOO problems with three objectives. However, clustering a 3D PF will be an NP-hard problem as soon as the general 2D cases are proven to be NP-hard. The perspective in such cases is to design specific approximation algorithms for a 3D PF.

Conflicts of Interest:
The authors declare no conflict of interest.