Polynomial algorithms for p-dispersion problems in a planar Pareto Front

In this paper, p-dispersion problems are studied to select $p\geqslant 2$ representative points from a large 2D Pareto Front (PF), solution of bi-objective optimization. Four standard p-dispersion variants are considered. A novel variant, Max-Sum-Neighbor p-dispersion, is introduced for the specific case of a 2D PF. Firstly, $2$-dispersion and $3$-dispersion problems are proven solvable in $O(n)$ time in a 2D PF. Secondly, dynamic programming algorithms are designed for three p-dispersion variants, proving polynomial complexities in a 2D PF. Max-min p-dispersion is solvable in $O(pn\log n)$ time and $O(n)$ memory space. Max-Sum-Neighbor p-dispersion is proven solvable in $O(pn^2)$ time and{$O(n)$} space. Max-Sum-min p-dispersion is solvable in $O(pn^3)$ time and $O(pn^2)$ space, this complexity holds also in 1D, proving for the first time that Max-Sum-min p-dispersion is polynomial in 1D. Furthermore, properties of these algorithms are discussed for an efficient implementation {and for a practical application inside bi-objective meta-heuristics.


Introduction
In real-world applications, optimization problems may be driven by several conflicting objectives. Designing network or (system of) system architectures, financial costs must be traded off with quality of service or robustness [12,43]. Dealing with complex maintenance planning problems, stability and robustness of the planning matter as well as financial costs [19]. Multi-objective optimization (MOO) supports such decision making. Many efficient (i.e. best compromise) solutions of MOO problems may exist with Pareto dominance [20]. A Pareto Front (PF) denotes the projection of efficient solutions in the objective space. This work aims to select p solutions from n p non dominated solutions, while maximizing the diversity of these p solutions in the objective space. Firstly, such problem occurs when selecting alternatives for decision makers. Secondly, MOO approaches use such operators to represent large PF [4,47], and MOO meta-heuristics archive diversified solutions during the heuristic search [49,52,57]. Thirdly, a similar problem occurs in a Skyline for databases [5,8,36]. When selecting alternatives for decision makers, p is small, p 5 is realistic. Otherwise, p is larger, having p = 100 or p = 1000 is realistic.
The hypervolume measure is often used in such a context [3,24,29]. Covering and clustering algorithms [14,16,57] are also used to select points in a PF. In this paper, we consider (discrete) p-dispersion problems, to select p points and while maximizing dispersion measures among selected points [22]. Although p-dispersion is mentioned to have relevant applications for MOO [23,45], no specific studies concerned the p-dispersion in a PF to the best of our knowledge. Four variants of discrete p-dispersion problems are defined [23]. Max-min and Max-Sum p-dispersion problems, the most studied variants, are N P-hard [22,30]. Max-min p-dispersion maximizes the minimal distance between each pair of selected points. Max-Sum p-dispersion maximizes the total sum of the distances between selected points. Max-Sum-min p-dispersion variant maximizes the sum of the distances between each selected points to the closest different selected points. Max-min-Sum p-dispersion variant maximizes the minimal sum of distances between each selected point to the other selected points. This paper studies p-dispersion variants in the case of a two-dimensional (2D) PF, which is an extension of one-dimensional (1D) cases. A novel variant of Max-Sum-min p-dispersion, denoted Max-Sum-Neighbor p-dispersion, is specifically introduced for 2D PFs. For these five p-dispersion problems, the cases p = 2 and p = 3 are proven to be solvable in O(n) time and the cases p = 4 and p = 5 are solvable respectively in O(n 2 ) and O(n 3 ) time. Generally, the Max-min, Max-Sum-min and Max-Sum-Neighbor p-dispersion problems are proven to be solvable in polynomial time in a 2D PF with dynamic programming (DP) algorithms. For Max-Sum-min p-dispersion, it is the first time that this problem is proven to be polynomially solvable in 1D.
The remainder of this paper is organized as follows. In Section 2, we introduce the notation and formally describe the problems. In Section 3, we discuss related state-of-the-art elements to situate our contributions. In Section 4, intermediate results are presented. In Sections 5, 6 and 7, DP algorithms with a polynomial complexity are respectively presented for the Max-Sum-Neighbor, Max-min and Max-min-Sum variants. In Section 8, our results are discussed from both a theoretical and a practical point of view. In Section 9, our contributions and the open perspectives are summarized.

Problem statement
Let E = {x 1 , . . . , x n } be a set of n points in a 2D PF, considering the minimization of two objectives (this is not a loss of generality). We denote discrete intervals [[a, b]] = [a, b]∩Z, so that we can use the notation of discrete index sets and write E = {x i } i∈[ [1,n]] . As in [14,16], PFs are formally defined using binary relations: relation I expresses Pareto incomparability, whereas relation ≺ defines an order from left to right, as illustrated in Figure 1.
Binary relations I, ≺ are defined 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) y I z ⇐⇒ y ≺ z or z ≺ y. A 2D PF can be the projected costs of efficient solutions using exact approaches in discrete MOO problems [20], or using population meta-heuristics like an evolutionary algorithm (EA) [52]. In the context of databases, Skyline operators are also PFs [5] A 2D PF can be extracted from any subset of R 2 using an output-sensitive algorithm [42]. MOO problems with continuous variables may have as solution a continuous PF. We will discuss in Section 8.4 how to use the results of this paper in this last case. Finally, an affine 2D PF is similar to a 1D instance, we will formalize it.
A strong assumption in this paper is that the 2D PF E is known a priori: p-dispersion problems are computed knowing E. This version of the problem is denoted explicit, unlike the so-called implicit versions where points of the PF are selected by simultaneously calculating the PF. In the context of MOO optimization, implicit p-dispersion would be combined with MOO approaches as in [4]. Context of databases are less structured: an implicit version of p-dispersion would consider the total number of rows in the database h n that can be enumerated where n is the size of the Skyline, which is unkown at the beginning of the search, as in [8].
where α > 0 and d(y, z) denotes a Chebyshev or a Minkowski distance, induced by the 8 and m norms. For a given m > 0, Minkowski distance d m is defined by the following formula for y = (y 1 , y 2 ), z = (z 1 , z 2 ) ∈ R 2 : The case m = 2 corresponds to the Euclidean distance. The limit with m → 8 defines the Chebyshev distance, denoted d 8 : In p-dispersion problems, the task is to select p 2 out of n given candidate points, while maximizing a dispersion function f : where D p denotes the set of all the p-tuples with distinct points of E: The most standard p-dispersion problem is also denoted Max-min pdispersion problem or p-dispersion-Mm. The dispersion function is in this case the minimum of distances d ij between pairs of the selected points. The Max-min p-dispersion problem for p 2 is written as: Max-Sum(-Sum) dispersion problem, denoted p-dispersion-MS, considers as dispersion function the total sum of distances among selected points: We consider also the Max-Sum-min p-dispersion variant, denoted p-dispersion-MSm [23], where dispersion is measured as the sum of the distance of each selected points to its closest (and different) selected point: The Max-min-Sum p-dispersion problem [23], denoted p-dispersion-MmS, is defined with a dispersion calculated as the sum of distances of each selected point to the other selected points: Note that 1D instances are special cases of 2D PF, equivalent to aligned points in a 2D PF. Considering any variant of p-dispersion problems, only the relative distance matters. Hence, p-dispersion problems for aligned points in the plane are equivalent to 1D instances and in these cases, Minkowski and Chebyshev distances are the same.
The previous definitions are generic and do not use specificities of 2D PFs. Lemmas 1 and 2, mentioned and proven in [16], allows to reformulate the definitions of some p-dispersion variants in a 2D PF, as well as defining a new variant that will be specific for 2D PFs: [1,n]] be a 2D PF. The relation is an order and ≺ is transitive. E can be re-indexed in O(n log n) time such that: [1,n]] be a re-indexed 2D PF usinf Lemma 1. The following monotony relations are valid considering a Minkowski or Chebyshev distance d, and any real number α > 0: Lemma 1 defines a 1D structure and a total order in a 2D PF. The reindexing in Lemma 1 is equivalent to a lexicographic sort minimizing hierarchically the two objectives, thus running inO(n log n) time. We will not use Lemma 1 when the O(n log n) time complexity degrade the overall complexity of some algorithms. When the indexing of a 2D PF fulfill equations (12), it will be mentioned as a re-indexed 2D PF. Without additional precision, a 2D PF will not be considered as re-indexed with Lemma 1.
Illustration of a 2D PF with 15 points and the indexing implied by Lemma 1 Using Lemma 2, Max-min and Max-Sum-min p-dispersion problems can be reformulated in a 2D PF considering only consecutive distances: [1,n]] be a re-indexed 2D PF. The Max-min and Max-Sum-min p-dispersion problems in E are also defined as: Proof : In the inner minimization of (8) and (10), distances d i j ,i j are considered. Such distances are higher than d i j ,i j+1 and d i j −1 ,i j using Lemma 2, and it remains only distances among consecutive points. For the two extreme points x 1 and x p , it remains only d i 1 ,i 2 and d i p−1 ,ip .
Furthermore, Lemma 1 allows to define a new variant of Max-Sum-min p-dispersion, the Max-Sum-Neighbor p-dispersion problem: [1,n]] be a re-indexed 2D PF. Max-Sum-Neighbor p-dispersion, denoted p-dispersion-MSN, is defined in E summing only the distances between neighbor points: Dispersion variants are illustrated in Figure 3. Lemma 3 shows that pdispersion-MSm is not a symmetric expression of distances, inducing more importance to extreme distances AB and CD. On the contrary, p-dispersion-MSN induces symmetrical expressions. :

Related works
3.1. Complexity results for p-dispersion problems. Max-min and Max-Sum p-dispersion problems are N P-hard for general metric spaces. This is proven in both cases by a polynomial reduction from the maximum independent set problem [22,30]. Max-min and Max-Sum p-dispersion problems are still N P-hard problems when distances fulfill the triangle inequality [22,30]. The planar Max-min p-dispersion problem is also N P-hard [56], the N P-hardness of the planar Max-Sum p-dispersion problem is still an open question to our knowledge.
Approximability and non approximability with constant factors were studied for p-dispersion problems. Unless P = N P, Max-min p-dispersion cannot be approximated with a constant factor [45]. For general metric spaces, Max-min and Max-Sum p-dispersion problems can be approximated with a 1/2 factor [53,54]. For 2D spaces, Max-min and Max-Sum p-dispersion problems can be approximated with a constant factor [45]. The 1/2 factor is the tightest approximation factor for Max-min p-dispersion with triangle inequality, unless P = N P [45].
Some sub-cases of p-dispersion problems are solvable in polynomial time. The 1D cases of Max-min and Max-Sum p-dispersion problems are solvable in polynomial time, with DP algorithms running in O(max{pn, n log n}) time [45,56]. Within a tree structure, Max-min p-dispersion is solvable in O(n 2 log n) time [10]. Although, p-dispersion is mentioned to have relevant applications for MOO [23,45], no specific studies concerned p-dispersion in a PF besides 1D cases to the best of our knowledge. We note an interest for implicit versions of dispersion problems in of Skyline Operators [36,37,55].

3.2.
Exact methods for p-dispersion problems. Max-Sum p-dispersion can be formulated as a quadratic optimization problem, defining binary variables z j ∈ {0, 1} with z j = 1 if and only if the point x j is selected: Linearizing (19.1) leads to the Integer Linear Programming (ILP) formulation provided in [33]. Exact Branch&Bound (B&B) algorithms were also provided computing iteratively higher and lower bounds, with a Lagrangian relaxation [2] or with tailored higher bounds computable in O(n 3 ) time [44].
Max-min p-dispersion is also a non linear optimization problem [44]: The standard linearization of constraints (20.3) leads to the Mixed Integer Linear Programming (MILP) formulation [33]. An alternative MILP formulation and specific cuts for a Branch&Cut algorithm is provided by [46]. Decomposition schemes from [2] and [44] have been also extended for the Max-min p-dispersion problem. Similarly, MILP formulations were designed for the Max-Sum-min and Max-min-Sum p-dispersion variants [23]. Such variants were less studied. A recent work proposed a unified MILP formulation and B&B algorithm for the four variants of p-dispersion problems [35].
3.3. Clustering/selecting points in PFs. We summarize here results related to the explicit versions of selection or clustering of points in a PF.
Maximizing the quality of discrete representations of Pareto sets was studied with hypervolume measure in the Hypervolume Subset Selection (HSS) problem [3,47]. The HSS problem, maximizing representativeness of k solutions among a PF of size n, is N P-hard in 3D (and higher dimensions) [6]. An exact algorithm in n O( √ k) in 3D and a polynomial-time approximation scheme for any constant dimension d were also provided [6]. 2D PF instances are solvable in polynomial time, a DP algorithm running in O(kn 2 ) time and using O(kn) space was firstly provided in [3]. The time complexity was improved in O(kn + n log n) by [7] and in O(k(n − k) + n log n) by [34]. Some similar results exist also for clustering problems. The k-median and k-medoid problems are N P hard in 2D since [40]. The 2D PF cases are solvable in O(n 3 ) time with DP algorithms [17,14]. The 1D cases are solvable in O(nk) time [31]. Min Sum of Square Clustering (MSSC), also denoted k-means problem, is also N P-hard for 2D cases [38]. The 1D cases of k-means are also solvable by a DP algorithm, with a complexity in O(kn) using memory space in O(n) [28].
Similar results are available for variants of p-center problems. The pcenter problems minimize the radius of a ball, to cover all the points with p such identical balls. Contrary to the continuous version, discrete p-center variant consider the additional constraint to have a discrete set of points as candidates for ball centers. It makes sense to consider the original points as such candidates for centers, as in [16]. The dual of a Max-min p-dispersion is similar to a min-max optimization as in p-center problems. Duality relations hold between p-dispersion-Mm and p-center problems in 1D and in tree structures: p-dispersion-Mm is the dual of the continuous (p − 1)-center for such cases [50]. Max-min p-dispersion and p-center problems have similar complexity results. The discrete and continuous p-center problems are N P-hard in general, the discrete p-center problem in R 2 with a Euclidean distance is also N P-hard [40]. The 1D and 2D PF sub-cases of p-center problems are polynomially solvable with DP algorithms. The 1D continuous p-center is solvable in O(n log 3 n) time [41], whereas the time complexity in a 2D PF is in O(pn log n) [15,16]. The discrete p-center problem in a 2D PF is solvable in O(n log n) time and O(n) space [8], whereas it is solvable in O(n) time in 1D [25]. The Min-sum p-radii problems, denoted also min-sum-diameter, are p-center variants, where the covering balls may not be identical. It is a min-Sum-Max optimization, that we can compare to the Max-Sum-min optimization with p-dispersion-MSm. Min-sum p-radii is N P-hard in the general case and polynomial within a tree structure [13]. The N P-hardness is also proven even in metrics induced by weighted planar graphs [26]. In a 2D PF, Min-sum p-radii is solvable in O(pn 2 ) time with a DP algorithm [16]. In 1D, a specific algorithm runs in O(n log n) time [16].
Lastly, partial variants were studied for p-center problems and variants in a 2D PF, allowing that m points are not covered. These points can be outliers to remove or isolated points that are wished to be detected in the context of EAs [16]. DP algorithms for 2D PF can be extended allowing to uncover m > 0 points. The time and space complexity of DP algorithms are then multiplied by a factor m [16]. In particular, partial p-center problems in a 2D PF are solvable in O(mpn log n) time and O(mn) space [16].

Summary of contributions and relations to state of the art.
From this literature review, some common results appear. General 2D instances of the clustering and selection problems are often N P-hard. DP algorithms induce a polynomial complexity in 1D and some 2D PF cases. Many DP algorithms use a p × n matrix to store results of k p selection/clustering among the n n first points. This induces p × n computations of a partial solution using O(n) previous ones. A linear enumeration induces time complexities in O(pn 2 ) whereas some O(pn log n) time complexities can be obtained using a logarithmic search. Our DP algorithms for p-dispersion-MSN and p-dispersion-Mm have this form. Other techniques make it possible to divide the complexity by a factor of p, as in [8,28]. It is an open perspective for the DP algorithms provided in this paper. Table 1 summarizes the complexity results for selection and clustering problems for the 1D, 2D PF and 2D sub-cases. A first comparison between these problems situates the selection problems p-center, p-dispersion-Mm and HSS as the fastest to solve in 2D PFs. Clustering with p-median or p-medoids is more accurate for measuring cluster similarity, however the time complexity in O(n 3 ) is a burden for application. Time complexity in O(pn 3 ) for p-dispersion-MSm is a more a theoretical than a practical result. Surprisingly, p-dispersion-MSN, which seems to be a slight variant of p-dispersion-MSm, has a much better complexity. Min-sum p-radii and p-dispersion-MSN have the same time complexity in O(pn 2 ) with similar DP algorithms. Both variants improve weaknesses of max-min and minmax optimization with p-center and p-dispersion-Mm, such as the possible large number of optimal solutions, including solutions potentially very unbalanced, to obtain fewer and better balanced solutions. Such weaknesses of p-dispersion-Mm are highlighted in Proposition 10, similar algorithms and results for p-center problems are described in [16].
Complexity results can be significantly worse for 2D PFs than the 1D sub-cases, as for the Min-sum p-radii and p-medoids. In such cases, triangle inequality instead of additivity of distances is crucial. For p-center and p-dispersion problems, the time complexity for 2D PFs is not significantly worse than for 1D sub-cases. Within a tree structure, which is another extension of 1D cases, p-dispersion-Mm is solvable in O(n 2 log n) time [10] instead of O(pn log n) in 2D PF, the time complexity in a 2D PF is significantly better. Note that Max-Sum-min p-dispersion was not studied before in 1D to the best of our knowledge, so that Theorem 3 and Corollary 1 proves for the first time that Max-Sum-min p-dispersion in 1D is solvable in polynomial time. Perspectives may be to improve this complexity result using distance additivity in 1D.

Intermediate results
This section presents intermediate results that will be a basis for future developments. A key element is that the extreme points of a re-indexed 2D PF are natural candidates for p-dispersion problems: Table 1. Comparison of the time complexity obtained after our study for dispersion problems (in bold) with related problems on 1D, 2D PF and 2D cases and their reference. Some cases are still open questions. BF denotes brute force naive enumeration. Note that p-dispersion-MSN is defined only for 1D and 2D PF instances, not for general 2D instances Proof : Any 2-dispersion variant consider the same problem:  Proof : Let 1 i 1 < i 2 < · · · < i p p be the indexes defining an optimal solution of a p-dispersion variant. Considering new indexes . Points x i 1 , . . . , x i p−1 , x i p have at least the same dispersion than the original points x i 1 , . . . , x i p−1 , x ip . Hence, it defines an optimal solution of the considered p-dispersion variant.
. , x i p−1 , x i p would induce a strictly better solution than an optimal solution of Max-Sum or Max-Sum-Neighbor p-dispersion. By contradiction, the extreme points are in any optimal solution for both variants.
Remark: For Max-min p-dispersion, optimal solutions exist without containing x 1 and x n . We consider 3-dispersion-Mm and 4 points x 1 = (0, 10), Remark: Clustering measure like k-means, k-medoids or k-center variants do not return the extreme points in general.
Proof : Considering variants (8), (9), (10), (11) or (18)  A specific result holds for p-dispersion-MSN in 1D. Proposition 5 and its proof show that it makes sense to consider this problem only for α = 1: Proof : Using Proposition 2, one selects the two extreme points. Let a = min E and b = max E, a and b are computed in O(n) time. With α = 1, adding any subset of size p − 2 of distinct elements of E \ {a, b}, we will have the same MSN dispersion of b − a > 0 because of the distance additivity in 1D, which is trivially optimal.
. Because of Proposition 2, we can assume that j 1 = 1 and j k = i. Necessarily, j 1 , j 2 , . . . , j k−1 defines an optimal solution of (k−1)-dispersion-MSN among points indexed in [[1, j k−1 ]]. On the contrary, a strictly better solution for C M SN k,i would be constructed adding the index i. We have thus: Proof : Let p 3 and let us prove the validity of Algorithm 1. Induction formula (23) uses only values C i,j with j < k. Hence, C k,n is at the end of each loop in k the optimal value of k-dispersion-MSN among the n points of E, and the optimal cost is given by C p,n . The remaining operations consist in a standard backtrack algorithm to return an optimal solution. This proves the validity of Algorithm 1 to solve optimally p-dispersion-MSN.
Let us analyze the complexity. Re-indexing E following Lemma 1 has a time complexity in O(n log n). Computing the line k = 2 of the DP matrix has also a time complexity in O(n). for i = k to n − 1: add j in J i := j end for return OP T the optimal cost and the set of selected indexes J .
to delete the line k−1 once the line k is completed. Such implementation has a spatial complexity in O(n) with at most 2n elements in memory. One may have only one line in memory, with an "in-place" implementation, computing values C k,m for a given k with index m decreasing: in-place m -th values for m < m are still C k−1,m that are needed to compute C k,m . Algorithm 1 has a spatial complexity in O(pn) because the backtracking operations use the full DP matrix. Algorithm 2 has a O(n) memory space algorithm by adapting techniques that were used in [11,21]. Algorithm 2 stores an intermediate value in the middle of the path of an optimal solution, to recover an optimal solution with a recursive divide and conquer strategy which will not be penalizing for the asymptotic time complexity. Such recursion applied to index j = argmax j ∈[[p−1,N −1]] (C p−1,j + d j ,N ) is valid, but it would lead to a O(p 2 n 2 ) time complexity to have a space complexity on O(n).  2 ). By induction, the optimal solution is concatenated using that x b is an optimal selected point, having an optimal solution of p -dispersion-MSN calling DivideConquer(E, a, h, p ) and it remains to compute a (p−p −1)dispersion-MSN between x h and x j . This proves the validity by induction.
Space complexity is in O(n) with C and H vectors in the first cost computation, and thereafter the space usage decreases (C and H are deleted before the recursive calls). Key point is the time complexity. Let T (n, p) the computation time to compute p-dispersion-MSN among n points. Using Proposition 7, there exists β such that βpn 2 is an upper bound for the computation of OPT. Hence, We have following induction relation: (27) T (n, p) βpn 2 + T (H p,n , p ) + T (j − H p,n , p − p − 1) By induction, we can prove that it exists γ 2 × β such that for all n, p, T (n, p) γpn 2 . It is true for terminal conditions and by induction: Hence, Algorithm 2 runs in O(pn 2 ) time using O(n) space.
Remark: Using Algorithm 2 instead of Algorithm 1, if we have the same asymptotic time complexity, the number of operations (and thus the CPU time) is approximatively doubled. As mentioned by [21], this should be used only if memory is missing to use Algorithm 1. ] min C M m k−1,j , d j,i . Let j 1 < j 2 < · · · < j k−1 < j k be indexes defining an optimal solution of k-dispersion-Mm, its cost is C M m k,i . Using Proposition 2, we can assume that j 1 = 1 and j k = i. Let c be the (k − 1)-dispersion-Mm of points indexed by j 1 < j 2 < · · · < j k−1 . We have c C M m k−1,j k−1 . If c d j k−1 ,j k , the bottleneck distance is given by points indexed by j k−1 , j k and we have C M m k,i = d j k−1 ,j k = min(C M m k−1,j k−1 , d j k−1 ,i ). Otherwise, j 1 , j 2 , . . . , j k−1 define an optimal solution of (k

p-dispersion-Mm is polynomially solvable in a 2D PF
k−1,j is increasing: any feasible solution of (k − 1)-dispersion-Mm among the j first points, is a feasible solution for (k − 1)-dispersion-Mm considering the j + 1 first points, and the optimal value C M m k−1,j is increasing. Hence, For j α, ψ i,k (j) = d j,i , and ψ i,k is strictly decreasing for j α. For j < α, ψ i,k (j) = C M m k−1,j , and ψ is increasing for j < α. Hence, ψ i,k reaches a maximum for j = α or j = α − 1. The computation of α, as the minimal value such that ϕ i,k (j) 0, can be solved with a dichotomic search presented in Algorithm 3, for a time complexity in O(log(i + 1 − k)).
To have a linear space complexity, one can design a recursive DP algorithm as in Algorithm 2. For Max-min p-dispersion, simple greedy algorithms in Algorithms 4 and 4' are valid as backtracking procedures: Once the optimal cost of Max-min p-dispersion problem is computed, Algorithms 4 and 4' compute an optimal solution in O(p log n) time using O(p) additional memory space. Furthermore, let j 1 = 1, j 2 , . . . , j p−1 , j p = n be the indexes of an optimal solution, let 1, i 2 , . . . , i p−1 , n (resp 1, i 2 , . . . , i p−1 , n) be the indexes given by Algorithm 4 and 4'. We have: In other words, the indexes given by Algorithm 4 and 4' are lower and upper bounds of the indexes of any optimal solution of Max-min p-dispersion considering the extreme points x 1 and x n . Proof : We prove the result for Algorithm 4, proof for Algorithm 4' is similar. Let j 1 = 1, j 2 , . . . , j p−1 , j p = n be the indexes of an optimal solution, let i 1 = 1, i 2 , . . . , i p−1 , i p = n be the indexes given by Algorithm 4. Firstly, we prove by induction on k that for all k ∈ [[1, p − 1]], i k j k . The case k = 1 is given by j 1 = i 1 = 1. We suppose k > 1 and that the induction hypothesis is true for k − 1, i.e. i k−1 j k−1 . The index i k is the smallest index such that d(x i k , x i k−1 ) OP T . Using Lemma 2 and i k−1 j k−1 , d(x j k , x j k−1 ) d(x j k , x i k−1 ). Having i k > j k would be in contradiction with d(x j k , x j k−1 ) OP T and the definition of i k as the smallest index such that d(x i k , x i k−1 ) OP T . We have also i k j k , which terminates the induction proof, indexes i k are lower bounds of indexes j k .
Let us prove that indexes i 1 = 1, i 2 , . . . , i p−1 , i p = n define an optimal solution. By construction d( Let us analyze the complexity. Algorithm 4 calls at most p − 2 times the computation of smallest index i k such that d(x i k , x i k−1 ) OP T , which can be proceeded with a dichotomic search, it runs in O(p log n) time.
Using Proposition 10, Algorithm 5 is a valid DP algorithm for Max-min p-dispersion running in O(n) memory space. Using Algorithms 4 and 4' instead of a divide-and-conquer strategy avoids to double the computation times , as noticed for p-dispersion-MSN. Theorem 2 summarizes the complexity results for Max-min p-dispersion: [1,n]] be a 2D PF. Max-min p-dispersion is polynomially solvable to optimality in the 2D PF E. Cases p = 2, 3 are re-index E following the order of Lemma 1 initialize vector C with

p-dispersion-MSm is polynomially solvable in a 2D PF
To have a DP algorithm for Max-Sum-min p-dispersion, more adaptation is needed as the cost computation of adding a new point i depends on the choice of the next selected point i > i. To design a DP algorithm, we do not store partial optimal solution of Max-Sum-min k-dispersion in a subset of points. We define C M Sm k,i,i as the best partial cost of k-dispersion-MSm in with i < i knowing that i is the last selected point before i and without counting a cost for point i . Proof : Using Proposition 2, C M Sm

3,i,i
is defined selecting 1, i, i , this makes the partial dispersion removing the d i,i term as given in (32).
an optimal partial solution of (k − 1)-dispersion-MSm among points indexed in [ [1, i]] with j as last selected point before i, and adding point i, it makes a feasible solution for (k − 1)-dispersion-MSm among points indexed in [ [1, i]] with a partial cost C M Sm k−1,j,i + min(d j,i , d i ,i ). This last cost is lower than the optimum C M Sm Let j 1 , . . . , j k be indexes such that 1 j 1 < j 2 < · · · < j k−1 = i < j k = i defining an optimal partial solution of k-dispersion-MSm in [ [1, i ]] with i as last selected point before i , its cost is C M Sm k,i,i . Necessarily, j 1 , j 2 , . . . , j k−1 defines an optimal partial solution of (k − 1)-dispersion-MSm among points indexed in [ [1, i]] with j k−2 as last selected point before i. On the contrary, a strictly better solution for C M Sm k,i,i would be constructed adding the index i . We have thus: Bellman equations of Proposition 11 allow to solve p-dispersion-MSm in a 2D PF with a DP algorithm detailed in Algorithm 6. Once DP matrix C M Sm k,i,i is computed, the optimal value of the complete Max-Sum-min p-dispersion is the best value C p,j,n + d j,n for j < n. A corollary is that these algorithms also apply in 1D, proving for the first time the polynomial complexity of Max-Sum-min p-dispersion in 1D: re-index E following the order of Lemma 1   Proof : By applying Algorithm 6 and Theorem 3 for x 1 ), . . . , (0, x n )}, it would give the result with a valid algorithm. However, this degenerate case of PF is not considered by our assumptions. Therefore, we use an alternative definition of E to conform to the assumptions of this article. Using the Euclidean distance and denoting M = max j∈[[1,n]] x j , we consider the affine 2D PF: .
With this definition, we still have: so that p-dispersion-MSm in E is the same problem than considering pdispersion-MSm in E .
Algorithm 6 uses a memory in O(pn 2 ). To compute only the optimal cost, a space complexity in O(n 2 ) is obtained deleting elements C k−1,j ,i when all the elements C k,i,i are computed. As for Algorithm 2, it is possible to decrease this space complexity in O(n 2 ), using recursion and alternative backtrack operations. Such algorithm is presented in Appendix.

Discussions
In this section, we discuss some theoretical insights and practical applications of Theorems 1, 2, 3 and Algorithms 1, 2, 5 and 6. Another way to have well-balanced solutions is to design a polishing heuristic starting from the solutions given by Algorithms 4 or 4'. One may use a local polishing procedure for a better balance considering 3dispersion-Mm optimizations for consecutive points. Each 3-dispersion-Mm computation is running in O(log n) using Proposition 9 as points are already sorted. Computing the optimal consecutive Max-min 3-dispersions runs in O(p log n), so that even with O(n) such iterations, the total complexity remains in O(pn log n) time and O(n) space.

Implementation and parallelization issues. Time complexity in
O(pn log n) or O(pn 2 ) can be satisfactory for large scale p-dispersion computations. For a practical speed-up, Algorithms 1, 2, 5 and 6 have useful properties for efficient implementation and parallelization. The bottleneck in time complexity is the construction of the DP matrices. The backtracking algorithm is essentially sequential, but with a low complexity; this phase is not crucial for the global efficiency. The initial sorting algorithm is significant in the computation times only for Max-min p-dispersion with small values of p. Standard parallelization of sorting algorithms applies, even using General Purpose Graphical Processing Units (GPGPU) [51].
The construction of the DP matrix requires independent computations to compute (k + 1)-dispersion costs from k-dispersion values. In Algorithm In all cases, the parallelization is straightforward in a shared memory environment like OpenMP or in a distributed environment using Message Passign Interface (MPI). Parallel implementation requires only p−3 synchronizations; this is a good property for the practical efficiency. Applying LPT (Lowest Processing Times) rules from [27] for load balancing among the operations that can be computed in parallel, it is better to calculate the most time-consuming operations first, starting from the highest values of i down to the lowest. DP Algorithms 1, 2 and 5 are cache-friendly for an efficient implementation. Indeed, it requires only k-dispersion values to compute (k + 1)dispersion costs. Since these k-dispersion previous values are called several times, it is crucial to access them quickly. Having the k-dispersion values in cache allows such quick access. As written in Algorithm 1, one have to cache two vectors of size n for an implementation keeping the previous and current lines of the DP matrix in cache. With in-place implementations as in Algorithms 2 and 4, one may cache only one vector of size n, which is useful when cache size becomes a limiting factor. In backtracking operations of Algorithm 1, the elements are likely no longer cached, inducing more access time. With the lower complexity of backtracking, this is not a problem.
Finally, we discuss the possibility of massive parallelization under GPGPU. Exhaustive enumerations as in p-dispersion-MSN (and the lexicographic optimization) is compatible with GPGPU parallelization. However, this is not the case for the dichotomic search in Algorithm 3. For Max-min p-dispersion, one may parallelize the O(pn log n) time version with OpenMP or the O(pn 2 ) time under GPGPU. Finally, note that a line by line computation of the DP matrix as in Algorithm 2 is useful for GPGPU parallelization: memory in the GPU hardware may be a limiting factor. 8.3. From 2D PF to 3D PF. In this paper, we considered 2D PF generated by bi-objective optimization. We analyze here the possible extension of the results for larger dimensions, for an application to MOO problems with three or more objectives. 3D PF are generated in many real-world optimization problems, for instance maximizing robustness and stability and minimizing financial costs of maintenance planning [19]. General 2D pdispersion problems is a sub-case of 3D PF: these are affine 3D PF. As in the proof of Corollary 1, a similar transformation applies to consider general 2D instances as 3D PFs. N P-hard complexity for such cases with Max-min p-dispersion implies that Max-min p-dispersion is also N P-hard for 3D PF. Unless P = N P, there is no hope to generalize DP algorithms with a polynomial complexity to PF in dimension three and more. The 1/2 approximation factor is valid in 3D (and higher dimensions) PF for Max-min and Max-Sum p-dispersion problems, as it holds for any metric space [53,54].
Lemmas 1 and 2 are fundamental results to design DP algorithms and are highly specific for 2D PF cases, no such total order exists in 3D and larger dimensions. Generally, this explains why clustering and selection problems are polynomial in 2D PF and N P-hard in larger dimensions. Proposition 1 eases the calculation of p-dispersion problems with the selection of extreme points, which is crucial for Propositions 2, 3, 4 and 5. An open question is whether Proposition 1 could be extended to 3 dimensions (and higher dimensions). Extreme points are generically defined by different lexicographic optimization with permutations of the considered objectives. An extension of Proposition 1 would be to analyze if optimal solutions of p-dispersion problems in a PF should necessarily contain such extreme points. Actually, the answer is negative, as shown by the following counter-example. For n = 5 and with 3-dispersion problems, we consider points (10, 0, 0) ; (5, 5, 0); (0, 10, 0); (9.99, 0.01, 1) and (0.01, 9.99, 1). There are here only two extreme points (out of the 3! = 6 lexicographic minimizations): (10, 0, 0) and (0, 10, 0), whereas each 3-dispersion variant has the same and unique optimal solution (5, 5, 0); (9.99, 0.01, 1) and (0.01, 9.99, 1).
Lastly, the possible use of 2D PF DP algorithms as heuristics is discussed for 3D PF and higher dimensions. In general, one can project any PF structure into 1D structures, for example with weighted-sum scalarization, Principal Component Analysis (PCA) or linear regression. It allows to initialize a local search approach with solutions obtained after a projection in 1D, as in [32]. The projection of a 3D PF to a representative 2D PF is difficult in general, affine 3D PF represent any planar instance without any regularity. This perspective would hold only for some specific 3D PFs.

8.4.
Having continuous 2D PFs. The first hypothesis of our paper was to have a 2D PF of size n. To address complexity results for p-dispersion problems (and also for k-means, k-medoids, k-center variants) in the general and in the 2D PF case, a finite number of points shall be considered. This hypothesis can be reformulated as "let n points from a 2D PF" to define the problem. In the context of PFs, this finite hypothesis may be in contradiction with the possibility of having infinite PFs in MOO problems. In continuous PFs, p-dispersion problems make sense as continuous optimization problems. If bounded discrete MOO problems induce finite PFss, this is not the case with MO Linear Programs (MOLP) or MO MILPs [20]. For MOLPs, PFs are given as a connected subset of the frontier of a polyhedron, which is described by a finite number of extreme points of this polyhedron. In 2D, extreme points define the extremity of consecutive segments. For MO MILPs, PFs may be composed of isolated points and PFs of sub MOLPs. Continuous MOO problems may also give an analytic formula of PFs. To address p-dispersion problems in such a context, one may sample regularly the continuous parts of PFs to have a good discrete approximation of the dispersion problem. By analyzing which maximum value of n induces reasonable computation times depending on the context, this guides the choice of the granularity of such a discretization. 8.5. Application to k-medoids/k-means clustering in a 2D PF. Exact DP algorithms for k-medoids run in O(n 3 ) time [14]. Such complexity is a bottleneck for the practical applications with large values of n. Classical heuristics for k-means/k-medoids problems may be used in such cases [9,48].
Such heuristics have no guarantee of optimality, and depend heavily on initialization strategies [9]. One may initialize local search with p-dispersion solutions to have a quick initialization procedure. Several such initialization strategies are possible. Firstly, one can initialize the k centroids using k-dispersion. Secondly, one can solve 2k + 1-dispersion, and select the k intermediate points following order of Lemma 1. Thirdly, one can solve 3kdispersion, and select the 3k + 1 intermediate points for k ∈ [[0; k − 1]]. For these strategies, one may use optimal DP algorithms or local search iterations of 3-dispersion as in Section 8.1. Based on the preprint version of this paper, numerical results were provided for randomly generated 2D PFs with n 5000. Good results and very quick computation times are obtained with such heuristics, whereas computing optimal solution with exact DP is time consuming for n = 5000 [32]. Having many different initialization procedures is useful for the final quality of solutions, implementing several local searches with different initialization strategies in parallel environments, as in [18]. This is also an interest of the numerical results presented in [32]. 8.6. Applications to MOO meta-heuristics. In the case of population MOO meta-heuristics like EAs, selection operators may be called iteratively among the current n non-dominated points to operate cross-overs or mutations (or applying a trajectory local search) among a restricted number of p n solutions [52]. Randomized selection operators may be used in that goal [52]. Deterministic strategies can also provide representative and diverse solutions, in addition to or as a replacement for random selection operators for some iterations. In such context, the number of points to select p may be large, regarding the maximal size of the archive, the maximal number of points that would be in memory. Selecting such p points algorithms running in O(pn log n) time are of major interest in such context, it is the case with Max-min p-dispersion, continuous p-center problems [16], and hypervolume subset selection [34]. Discrete p-center is even faster to compute in O(n log n) time [8]. Note that it is not required to have optimal solutions of such problems inside EAs. Even k-medoids heuristics may be used as in [32] as long as the calculation time of the heuristic remains fast, complexity calculations guide for such an heuristic design. Polishing procedures are of interest to have better balanced solutions as presented in Section 8.1. The fastest DP algorithms may be used to initialize local search heuristics for other problems, as in [32].
Having several types of deterministic selection operators is of interest for EAs to diversify the points where mutations or intensification are operated among consecutive iterations. On the contrary, it is a source of inefficiency to call the operators of mutations and selections always on the same points. Varying the value of p gives a first solution of this problem for deterministic operators. Selecting points in a 2D PF, many optimization measures and algorithms are available, and have different properties; another solution is to solve iteratively different selection problems. When selecting points in a PF, the rule of thumb is using hypervolume measures and variants [24,29]. HSS has useful properties to explain such popularity [29]. Clustering with kmeans or k-medoids define dense zones in the PF where little intensification is needed [14]. Partial p-center variants detect outliers (isolated points) simultaneously with a rough definition of clusters [16]. In a 2D PF, isolated points should be preferred to try intensification strategies in such zones, to have better-balanced points along the 2D PF [16]. Without the partial extension, p-center problems are faster to compute, but the clustering and selection are less relevant; fast DP algorithms for p-center problems are useful to design quick heuristics fr other selection or clustering optimization. HSS and p-dispersion measures find diverse points in the 2D PF, but can lead to very different solutions: HSS would avoid points that are close to local Nadir points (of knee-points) contrary to p-dispersion problems.
An open perspective is to design EAs combining randomized selection operators with different deterministic and complementary strategies.

Conclusion and perspectives
The properties of the four standard p dispersion problems have been examined in a 2D PF using Euclidean, Minkowski or Chebyshev distances. A novel variant, namely Max-Sum-Neighbor p-dispersion, is defined specifically for 2D PF. Cases p = 2 and p = 3 induce a complexity in O(n) time. Cases p = 4 and p = 5 have respectively time complexities in O(n 2 ) and O(n 3 ) time. Such results are useful to select a small number of representative solutions for decision makers.
Three variants are proven solvable in polynomial time in a 2D PF, with the design of DP algorithms. Standard Max-min p-dispersion problem is solvable in a 2D PF in O(pn log n) time and using O(n) memory space. Max-Sum-Neighbor p-dispersion is solvable in a 2D PF in O(pn 2 ) time and O(n) space. A lexicographic optimization considering Max-min and Max-Sum-Neighbor p-dispersion is also solvable in O(pn 2 ) time using O(pn) space. Max-Sum-min p-dispersion is solvable in a 2D PF in O(pn 3 ) time. This last result and the DP algorithm proves also that Max-Sum-min p-dispersion is polynomially solvable in 1D, which was never studied before. Perspectives may be to improve this complexity result using specificity of 1D instances. Considering Max-Sum p-dispersion, the N P-hardness of 2D PF and also 2D sub-cases are still open questions.
These results are not only of a theoretical interest, but raise also practical perspectives.Complexity for Max-min and Max-Sum-Neighbor p-dispersion allows a straightforward application for large 2D PF. Furthermore, DP algorithms have useful properties for an efficient implementation, including efficient parallelization in a multi or many-core environment. It allows an application inside MOO population meta-heuristics to archive partial PF at each iteration. In this context, p-dispersion DP algorithms may be used also to initialize k-medoids clustering in a 2D PF. The extension of results to higher dimensions was discussed. 3D PF cases are N P-hard and approximation algorithms with factor 1/2 are available. Perspectives are only to design quick heuristics for PF in such dimensions. Lastly, the results of this paper may be extended to implicit versions of the problem related to Skyline Operator and MOO applications. thanks the anonymous reviewers for their helpful comments that allowed to improve the paper very significantly.