Minimal Radius Enclosing Polyellipsoids
aa r X i v : . [ m a t h . O C ] O c t Minimal Radius Enclosing Polyellipsoids
V´ıctor Blanco † and Justo Puerto ‡ † IEMath-GR, Universidad de Granada ‡ IMUS, Universidad de Sevilla
Abstract.
In this paper we analyze the extension of the classical smallest enclosing diskproblem to the case of the location of a polyellipsoid to fully cover a set of demand pointsin R d . We prove that the problem is polynomially solvable in fixed dimension and analyzemathematical programming formulations for it. We also consider some geometric approachesfor the problem in case the foci of the polyellipsoids are known. Extensions of the classicalalgorithm by Elzinga-Hearn are also derived for this new problem. Moreover, we addressseveral extensions of the problem, as the case where the foci of the enclosing polyellipsoid arenot given and have to be determined among a potential set of points or the induced coveringproblems when instead of polyellipsoids, one uses ordered median polyellipsoids. For theseproblems we also present Mixed Integer (Non) Linear Programming strategies that lead toefficient ways to solve it. This article is dedicated to Prof. Marco A. L´opez on the occasion of his 70th birthday. Introduction
Given a sets of demand points in a given finite dimensional normed space R d , ContinuousFacility Location Problems (CFLP) deal with the determination of the optimal placement ofsome new points in order to minimize certain measure of the distances from the given to thenew points. The most popular CFLP is the Weber Problem [35] in which a single facility is tobe located minimizing the overall (weighted) sum of the Euclidean distances from the demandpoints to the facility. A common approach to solve the Weber Problem is via the Weiszfeldalgorithm [36]. The extension to multiple facilities may either consider that the overall sumfrom the demand points to all facilities (multiple-allocation case) or to its closest facility (single-allocation case) has to be minimized. Different aggregation measures, apart from the overallsum, have been also proposed to find optimal facilities for CFLP (see e.g., [3]).The level curves of the sum of Euclidean distances to given points on the plane (the Eu-clidean planar Weber level curve), in case the ( n ) demand points (also known as foci withinthis framework) are not collinear, are called in the literature polyellipses , n -ellipses or multifocalellipses . The regions bounded by those curves are known to be convex bodies and can be alsodefined in higher dimensional spaces and with different norm measures, inducing the notion of polyellipsoids . Several authors have been interested in the analysis of these convex bodies, fromMaxwell [22] to Er¨dos and Vincze [17], mainly from an algebraic or geometric point of view(see [27]). In particular, the problem whether all the convex plane curves can be arbitrarilyapproximated by polyellipses under a sufficiently large number of points (also known as foci )was posed by Weiszfeld. In this sense, the Erd¨os-Vincze’s theorem states that regular triangles Date : October 30, 2019.This research has been partially supported by Spanish Ministry of Econom´ıa and Competitividad/FEDERgrants number MTM2016-74983-C02-01. cannot be presented as the Hausdorff limit of polyellipses even if the number of the foci can bearbitrary large. Recently, in [32], the authors extend the analysis for regular polygons in theplane. However, these sets (polyellipsoids) have only been partially analyzed under optimizationlens beyond their connections with the Weber problem (see, e.g. [24]).On the other hand, Covering Location Problems consist of locating facilities to satisfy thedemand of a set of customers whenever they are only allowed to be served within certain coverageareas. In the continuous case there are two main families of problems in this field: Full CoveringLocation problems (FCP) and Maximal Covering Location problems (MCP). In the former allthe demands points have to be covered and the goal is to minimize the set up costs, whichmay consists of the cost of opening a facility and/or the cost of enlarging the coverage areaof the facilities to reach the demand points (see [2, 28]). In the latter (MCP), the goal is tolocate the facilities such that the (weighted) number of points minus the facilities set-up costsis maximized. Both problems have been analyzed in the literature, mainly on the plane. If thenumber of facilities to open p is given and the regions are defined as distance-based balls, the FCPis also called the continuous p -center problem or the minimum volume enclosing sphere, since theproblem can be interpreted as the one of locating p facilities in the space, such that the maximumdistance from the demand points to its closest facility (center) is minimized. The particular casein which a single facility is to be located and the distance measure is the Euclidean, Elzinga andHearn provided one of the most popular algorithms in continuous location for this problem [16].More recently, linear time algorithms have been proposed for the weighted minimum sphereenclosing problems ([13] and [23]). Also, the problem of finding the minimum volume enclosingellipsoid in fixed dimension is known to be solved in linear time [14]. On the other hand, for theFCP in case of locating a single circle-shaped facility on the plane there is always an optimalsolution of the problem whose position lies either in the intersection of certain circles centered atthe demand points, or at the demand points, the so-called circle intersection points (CIP) [8, 11].This property allows one to transform the original continuous problem into a discrete coveringlocation problem. The multifacility case is addressed in [10], proving the same CIP property, forboth the Euclidean and the ℓ -norm distance measures. In [1] and [6] the case of locating ellipticalregions in the plane is studied, by means of Mixed Integer Non Linear Programming formulationsand heuristic approaches. Apart from that, as far as we know, there are no further advances onthe topic. In particular, the definition and use of weighted d -dimensional polyellipsoids to fullycover a set of demand points have not been previously considered.In this paper, we extend the classical smallest enclosing disk problem (also known as minimaxfacility location problem ) to the case in which instead of spheres or ellipsoids, polyellipsoids(higher dimensional polyellipses) are to be located. Since polyellipses with a single focus coincidewith disks, our approach naturally extends the classical problems.Turning to the applications of these models, it is natural to think about locating several alarmsirens or facilities (which play the role of the foci) from where some emergency services shouldbe dispatched. In these situations the service quality decreases with distance, so the closer thebetter. Therefore, the problem can be seen as to adequately locate a set of facilities whoserelative positions are given and determine the minimum coverage distance for which the demandof all customers is satisfied. Note that the coverage area in this problem is defined as the regionfor which the sum of the distances from a customer to all the foci is less than or equal to thegiven threshold. In Figure 1-right we draw a situation in which the demand of a set of customers( ◦ ) has to be covered by nine facilities ( N ). The minimum radius coverage region is limited bythe drawn polyellipse. As opposed to the case where the coverage has to be done with respectto a single focus (circle) which is shown in Figure 1-left. Another field of applications, differentfrom Location Analysis, comes from Data Analysis and outliers detection. Minimum enclosingballs allow one to determine possible outliers in a dataset, looking at the points that do not inimal Radius Enclosing Polyellipsoids 3 Figure 1.
Minimum radius circle (left) and 9-ellipse (right) enclosing a set ofdemand points.fall inside the ball. However, many clouds of points do not fit well to a sphere but to a morecomplex shape. Thus, polyellipsoids with an adequate set of foci allow one a more appropriatefitting and a better outlier detection tool. A more sophisticated tool to detect outliers (andalso for one-class classification) that is becoming popular in Data Science is Support VectorData Description (SVDD) [31]. In SVDD, outlier observations are detected by enclosing part ofthe points inside a Euclidean ball, after introducing misclassifying errors and transforming thepoints into a higher dimensional space using kernels. Our methodology could be adapted to thisframework determining outliers observations for non-spherical data, avoiding the use of kerneltransformations.The paper is organized in seven sections. In Section 2 we introduce the framework of thepaper, recall the definitions and some preliminary results of polyellipsoids and formulate theproblem of determining the optimal position of a given set of foci such that a set of demandspoints is fully covered by the polyellipse with minimum radius. We also derive a polynomial-timecomplexity result for the problem. In Section 3 we describe several solution approaches for theproblems, based either on dualization of the problem or in decomposition approaches inspiredon the Elzinga-Hearn algorithm for the Euclidean 1-center problem. Section 4 is devoted toreport the results of an extensive computational experience which shows the efficiency of ourdecomposition approach. We analyze in Section 5 the very particular case in which the demandpoints are one-dimensional and derive a linear time algorithm (in the number of foci) for theminimal enclosing one-dimensional polyellipsoid. In Section 6 we show two possible extensions forthe problem. First we study the problem of selecting the foci from a potential set of candidateswhen they are not “a priori” given. Second, we present some results generalizing the notion ofpolyellipsoid, and the problem under analysis, to the ordered median framework. Finally, wedraw some conclusions and point out some further research topics in Section 7.2.
Minimal enclosing polyellipsoids with given foci
In this section we introduce the notation and concepts to be used for the rest of the paper.We recall the definition of polyellipsoid and derive some of its useful properties. We also analyzethe problem of locating a polyellipsoid with given relative positions of the foci (but unknownradius and position) to minimally cover a given set of demand points.We are given a finite set of points U = { u , . . . , u k } in R d and weights ω ∈ R k + , and a distancemeasure induced by a norm k · k . Without loss of generality, we assume that P u ∈U ω u = 1. The minisum or Weber location problem consists of finding the placement of a facility x ∈ R d thatminimizes the ω -weighted distances from x to the points in U , i.e., the function:Φ U ,ω ( x ) = X u ∈U ω u k u − x k , x ∈ R d . inimal Radius Enclosing Polyellipsoids 4 If the points in U are not collinear and the norm is strictly convex, Φ U ,ω is strictly convex, andtherefore has a unique minimum. The levels curves of Φ U ,ω , for r ≥
0, are given by the followingregions: E U ,ω ( r ) = n x ∈ R d : X u ∈U ω u k u − x k = r o . The set E U ,ω ( r ) is called a (weighted) polyellipsoid with foci U , weights ω , and radius r . Equal-weights Euclidean polyellipsoids on the plane are called polyellipses . Polyellipses were introducedfor the first time by Tschirnhaus in 1686. In case the number of foci is exactly one, polyellipsescoincide with circles, and for two foci, one obtains standard ellipses. In Figure 2 we draw someplanar polyellipses with five foci ( N ) and different radii. Figure 2.
Polyellipses with five foci and different radii.Let us denote by P U ,ω ( r ) the region bounded by the polyellipsoid E U ,ω ( r ), i.e., P U ,ω ( r ) = n x ∈ R d : X u ∈U ω u k x − u k ≤ r o . Abusing of notation, we will also call the sets P U ,ω ( r ) polyellipsoids.In the following we recall some known properties of the regions P U ,ω ( r ): Proposition 1.
Let U = { u , . . . , u k } ⊆ R d , ω ∈ R k + with P u ∈U ω u = 1 and r ≥ . Then: (1) P U ,ω ( r ) is a closed, bounded and convex set. (2) If x ∗ ∈ arg min x ∈ R d X u ∈U ω u k u − x k and r ∗ = P u ∈U ω u k u − x ∗ k : (a) x ∗ ∈ P U ,ω ( r ∗ ) . If k·k is strictly convex and U are not collinear then P U ,ω ( r ∗ ) = { x ∗ } . (b) P U ,ω ( r ) = ∅ for all r < r ∗ . (c) If P U ,ω ( r ) = { ¯ x } for some r > and ¯ x ∈ R d , then r = r ∗ and ¯ x = x ∗ . (d) x ∗ ∈ P U ,ω ( r ) for all r ≥ r ∗ . (3) If d = 2 and k · k is the Euclidean norm. Then, P U ,ω ( r ) can be represented as theset of solutions of a polynomial equation with degree k (if k is odd) or k - (cid:0) k k (cid:1) (if k is even) which can be expressed as the determinant of a symmetric matrix of linearpolynomials [27] . (4) P U ,ω ( r ) is contained in the ring with center at x ∗ and radii r − r ∗ P u ∈U ω u and r + r ∗ P u ∈U ω u The point x ∗ of the above result is known as the center of the polyellipsoid P U ,ω ( r ∗ ) and it canbe efficiently computed with classical location analysis algorithms as the hyperbolic modificationof Weiszfeld algorithm [36].Although the above result gives us information about the complexity of determining the pointswhich belong to a given polyellipsoid in terms of the number of foci, if the considered norm is inimal Radius Enclosing Polyellipsoids 5 a ℓ p -norm with p ∈ Q , p ≥
1, these convex bodies are not only convex sets, but Second OrderCone (SOC) representable sets (see e.g.[4]). The main implications of such a representation isthat optimization problems with linear constraints and SOC-representable feasible sets can beefficiently solved via interior point methods [25] available in most common optimization solvers.2.1.
The Minimum Radius Enclosing Pollyellipsoid Problem with Given Foci.
Givena finite set of demand points, A = { a , . . . , a n } ⊆ R d , the goal of the Minimum Radius EnclosingPollyellipsoid Problem with Given Foci (MEP, for short) is to determine the placement andthe minimum radius of the polyellipsoid with foci U = { u , . . . , u k } and weights ω ∈ R k + with P kj =1 ω j = 1, such that the points in A are fully covered. In other word, the optimal location ofthe polyellipsoid such that the largest sum of the weighted distances from a the demand point toall the translated foci is minimized. This problem naturally extends the widely studied 1-centerproblem in case a single focus is considered.Observe that the demand point a ∈ A belongs to the x -translation, for some x ∈ R d , ofa polyellipsoid P U ,ω ( r ) if a ∈ P x + U ,ω ( r ), where x + U = { x + u : u ∈ U} , i.e., if ϕ U a ( x ) := P u ∈U ω u k a − u − x k ≤ r .Thus, the problem can be formulated as:min x ∈ R d max a ∈A ϕ U a ( x ) . (MEP)Note that if x ∗ is an optimal solution of the above problem, it induces a polyellipsoid with foci U ∗ = { u + x ∗ : u ∈ U} , weights ω , and radius r ∗ = max a ∈A ϕ U a ( x ). Thus, x ∗ represents the optimal translation of the foci in U to minimally cover the points in A , with respect to theweighted ω -sum of the functions ϕ U a ( x ), ∀ a ∈ A . Observe also, that, as usual in minimaxproblems, one can reformulate the problem by introducing the auxiliary variable r (representingthe minimum radius) as follows:min x ∈ R d ,r ∈ R + r s.t. X u ∈U ω u k a − u − x k ≤ r, ∀ a ∈ A . It easily follows that the problem above can be formulated as the minimization of a linearobjective function subject only to linear or norm-type constraints:min r s.t. X u ∈U ω u d ua ≤ r, ∀ a ∈ A , k a − u − x k ≤ d ua , ∀ u ∈ U , a ∈ A , (MEP ) x ∈ R d , r ∈ R ,d ua ∈ R |U|×|A| , ∀ u ∈ U , a ∈ A . By the comments above, for ℓ p -norms, the nonlinear constraints in (MEP ) can be rewritten asa set of SOC constraints. The case in which the norm is a block norm is even simpler since theconstraints can be written as linear constraints. Thus, in the two cases the problem is solvableby interior point algorithms using any of the available optimization solvers.Note also that finding a feasible solution, ( x, r ), to (MEP ) is equivalent to find a point x inthe intersection of the n polyellipsoids P a + U ,ω ( r ) , . . . , P a n + U ,ω ( r ), since such a solution, ( x, r )must verify all the polyellipsoid constraints in (MEP ). In Figure 3 (left) we show a toy examplewith three demand points (stars) to be covered by a (Euclidean) polyellipse with three foci. Ifwe fix a radius and center the polyellipse at each of the demand points we obtain the three inimal Radius Enclosing Polyellipsoids 6 gray polyellipses in the picture. The point x ∗ belongs to the intersection of those polyellipsesand centering the polyellipse at that point, the polyellipse (dashed line) covers all the demandpoints (in this case being clearly a non optimal solution). In the right picture we draw the samesituation but with the optimal radius of the polyellipse. a a a x ∗ x ∗ + u x ∗ + u x ∗ + u a a a x ∗ x ∗ + u x ∗ + u x ∗ + u Figure 3.
Equivalence on the intersection of the polyellipses centered at thedemand points and a feasible covering polyellipse.The following result allows us to decompose (MEP) into smaller problems:
Theorem 1.
Let r ∗ be the solution of (MEP) . Then, r ∗ = min S ⊂A : | S | = d +1 r ∗ S , where r ∗ S is the optimalsolution of the following problem: r ∗ S := min r (MEP S ) s.t. X u ∈U ω u d ua ≤ r, ∀ a ∈ S, k a − u − x k ≤ d ua , ∀ u ∈ U , a ∈ S,x ∈ R d , r ∈ R ,d ua ∈ R |U|×|A| , ∀ u ∈ U , a ∈ S. for all S ⊆ A with | S | = d + 1 .Proof. The proof follows by a standard application of Helly’s theorem on the intersection ofconvex bodies [21]. The optimal solution of (MEP) reduces to finding the smallest r among thesolutions of the above collection of subproblems for all the subsets S ⊂ A with | S | = d + 1 (cid:3) This reformulation of the problem will be used both to derive a polynomial-time complexityfor the problem and also to develop an Elzinga-Hearn based solution approach.2.2.
On the complexity of (MEP) . In what follows we prove a polynomiality result for theminimum enclosing polyellisoid with a fixed number k of given foci problem in fixed dimension d .As usual in computational geometry, the model of computation is that of algebraic computationsand comparisons over the reals. The reader is referred to the papers by Renegar [30] and Dyer[14]for further details on the underlying results and the algorithms on which we base our construction.The main tool to be applied is the linear time algorithm of Dyer for a special class of convexprogram with few nonlinear constraints. We will prove that an iterative use of that algorithmimplies the polynomial-time complexity of the (MEP). Observe that the reformulation of (MEP)as (MEP ) has a number of nonlinear constraints that are essentially equivalent to polynomials.However, even assuming fixed dimension d and constant number k of foci, the number of nonlinear inimal Radius Enclosing Polyellipsoids 7 constraints is high as compared with the number of linear ones, and therefore one cannot applythe result of Dyer[14]. In spite of that, we will prove that a convenient decomposition of thatproblem would allow one to polynomially solve (MEP) in fixed dimension. Theorem 2.
Let the number of foci, k , and the dimension, d , be fixed. Then, (MEP) withEuclidean norm is solvable in polynomial time in |A| .Proof. Recall that by Theorem 1, solving (MEP ) (and then (MEP)) is equivalent to solve allproblems (MEP S ) for S ⊆ A with | S | = d + 1. Observe that in (MEP S ) the number of nonlinearconstraints is constant (( d + 1) × k ). Let us denote by K ( S ) = { ( x, d ) ∈ R d × R k × ( d +1) : k a − u − x k ≤ d ua , ∀ u ∈ U , a ∈ S } , the feasible domain induced by the nonlinear constraints in(MEP S ). Then, we have:
1. :
For any δ > a ∈ S and u ∈ U : g ua ( x, d ) := k a − u − x k − d ua ≤ ε, ∀ ε ∈ [0 , δ ] , and ∀ ( x, d ) ∈ K ( S ) . (1)Moreover, since k · k is assumed to be the Euclidean norm, the inequality g ua ( x, d ) ≤ ε can be easily transformed into the following polynomial inequality: P dj =1 ( a j − u j − x j ) ≤ ( ε − d ua ) .
2. :
The gradient ∇ g ua ( x, d ua ) = ( − − a − u + x k a − u − x k , − t . Therefore, the inequality ∇ g tau ( x, d au ) y ≤ y ∈ R d +1 is equivalent to the following polynomial inequality [ P dj =1 ( a j − u j − x j ) y j ] ≤ P dj =1 ( a j − u j − x j ) y d +1 .
3. :
The largest degree of all the polynomials involved in (MEP S ) and the transformationsmade explicitly in items [1.] and [2.] above is 4 (constant) and much smaller than |A| .Now, we are in position to apply the algorithm of Dyer to problem (MEP S ) because the numberof nonlinear constraints is constant. Thus, that problem is solvable in linear time of the numberof linear constrains, O ( d ), which in this applications turns out to be constant. Hence, because forsolving (MEP), one needs to solve O ( |A| d +1 ) problems in the collection and taking the maximum r value, the original problem can be solved in O ( |A| d +1 )-time. ✷ . (cid:3) For the sake of readability, the above result has been presented for the Euclidean norm.However, it extends with minimal changes to any ℓ p -norm with p ∈ Z and 1 < p < + ∞ . Observethat instead of squaring the norms expressions in g ua , one can reformulate them as the followingset of inequalities: v aj ≥ a j − u j − x j , ∀ j = 1 , . . . , d, a ∈ S,v aj ≥ − a j + u j + x j , ∀ j = 1 , . . . , d, a ∈ S, d X j =1 v paj ≤ d pua , ∀ a ∈ A . where the v -variables represent the absolute values of each of the coordinates involved in thenorm expressions. Thus, similar arguments to those in the proof of Theorem 2 can be derived incase p is small (compared to |A| ).Note also that if k · k is polyhedral the constraints in (MEP ) can be written as the followingset of linear inequalities: e t ( a − u − x ) ≤ d ua , ∀ u ∈ U , a ∈ A , e ∈ Ext k·k o where Ext k·k o = { e o , . . . , e og } are the extreme points of the polar ball of the unit ball of k · k (seee.g., [26, 34]). Thus, (MEP) can be solved using linear programming tools. inimal Radius Enclosing Polyellipsoids 8 Solution approaches for solving (MEP)In this section we describe different solution approaches for solving (MEP) beyond the solutionof the SOCP model provided in Section 2. We describe here two types of approaches. First, wedescribe those which exploit the strong duality of the problem under two different formulations :conic and Lagrangean dual. Then, we provide a geometric construction based on the classicalElzinga-Hearn algorithm for solving the problem.While the dual approaches allow us to provide alternative mathematical programming formu-lation for the problem and relate it with other classical problems (as the Weber problem), theElzinga-Hearn based approach allows, in practice, to efficiently solve geometrically the problem.3.1.
Conic Dual.
Again, for the sake of simplicity we consider the (MEP) with Euclidean norm.In such a case, the problem can be reformulated as the SOCP (MEP ). Therefore, one can deriveits conic dual, which results in:max X u ∈U X a ∈ A ( a − u ) (cid:0) λ ua + λ ua (cid:1) s.t. X a ∈ A µ a ≤ , X u ∈U X a ∈ A λ ua ≥ , X u ∈U X a ∈ A λ ua ≤ , (MEP ConicDual ) p ( λ ua ) + ( λ ua ) ≤ ω u µ a , ∀ u ∈ U , a ∈ A,µ a ∈ R + , λ ua , λ ua ∈ R , ∀ a ∈ A , u ∈ U . Slater condition ensures strong duality between both problems (see e.g. [5]) giving rise to simpleapproaches to solve the (MEP). Again, the case in which the norm is not Euclidean but ℓ p orpolyhedral permits a similar representation which results in a SOCP problem (as the above) ora linear problem, respectively.The main advantage, again, of the primal and the dual SOCP formulations of the problem isthat off-the-shelf software packages are capable to solve this type of problems efficiently usinginterior-point based algorithms. These algorithms applied to SOCP problems are known to havea polynomial time complexity for a given tolerance factor ε , assuring convergence in at most agiven number of iterations.Based on our computational experience, as we will see in Section 4, the SOCP formulationsare not only interesting because its theoretical complexity, but also because current solvers areable to handle medium size instances in a reasonable CPU time with reduced implementationeffort.3.2. Lagrangean Dual.
Now, we analyze a primal-dual approach based on the Lagrangeandual to solve (MEP). For the sake of simplicity, we assume that the norm is strictly convex,although as we will point out at the end of the section, the results extend to non strictly convexnorms. Let us denote by f a ( x ) = P u ∈U ω u Φ ua ( x ) = P u ∈U ω u k a − u − x k , for all a ∈ A . Clearly,each Φ ua is a convex function, and then f a is convex. inimal Radius Enclosing Polyellipsoids 9 Thus, the Lagrangean dual is: max α ∈ R |A| min x ∈ R d X a ∈A α a f a ( x )s.t. X a ∈A α a = 1 ,α a ≥ , ∀ a ∈ A . where α a for a ∈ A are the dual multipliers.Note that each of the functions f a attains its minimum in a solution of the Weber problemwith demand points A a = { a − u : u ∈ U} and weights ω , which is unique provided that thenorm is strictly convex and the foci in U are not collinear. Furthermore, if we denote by F ( α ) = min x ∈ R d X a ∈A α a f a ( x ) = min x ∈ R d X a ∈A α a X u ∈U ω u k a − u − x k , i.e., the optimal value of the Weber Problem with the nk demand points S a ∈ A A a and weights { α a ω u : a ∈ A , u ∈ U} , the dual problem becomes:max α ∈ R |A| F ( α )s.t. X a ∈A α a = 1 , (MEP LagDual ) α a ≥ , ∀ a ∈ A . Note that the feasible region is nothing but the |A| -dimensional probabilistic simplex, ∆ A = { α ∈ R |A| + : X a ∈A α a = 1 } .By strong duality: min x ∈ R d max a ∈A X u ∈U ω u k a − u − x k = max α ∈ ∆ A F ( α )It is straighforward that F is concave on ∆ A . Observe also that if k · k is a strictly convexnorm and the points in A (and in U ) are not collinear, for a fixed α ∈ ∆ A , the function h α ( x ) = P a ∈A α a f a ( x ), for all x ∈ R d , is strictly convex. Denoting by x ∗ α = arg min x ∈ R d h α ( x )(which is well-defined by the strict convexity of h α ) we have that F ( α ) = P a ∈A α a f a ( x ∗ α ) for all α ∈ ∆ A . Furthermore, the gradient of F can be easily derived: ∂F∂α a ( α ) = f a ( x ∗ α ) = X u ∈U ω u k a − u − x ∗ α k , ∀ a ∈ A . With such an information, since (MEP
LagDual ) is an optimization problem with a concave ob-jective function and a linear feasible set, a projected gradient descent method can be applied tosolve it. Let start with an initial solution α ∈ ∆ and the iterations are in the form: α k +1 = Π ∆ ( α k − η k ∇ α F ( α k )) . where Π ∆ ( β ) = arg min γ ∈ ∆ k β − γ k is the orthogonal projection onto the unit simplex which canbe efficiently computed in O ( n log n ) computation time (see [9, 12, 33]).Observe that at each step of this approach one has to compute the gradient ∇ α F ( α ) whichin turns implies solving the following Weber problem:min x ∈ R d X a ∈A α a X u ∈U ω u k a − u − x k . These problems can be arbitrarily approximated using Weiszfeld algorithm and its relatives.Thus, the optimal translation of the covering polyellipsoid, x , coincides with an optimal solution inimal Radius Enclosing Polyellipsoids 10 of a weighted Weber problem with demand points { a − u : a ∈ A , u ∈ U} and weights { α a ω u : a ∈ A , u ∈ U} , and then, as stated in [19], x belongs to the convex hull of those demand points.In case the norm k ·k is not strictly convex, uniqueness of optimal solutions cannot be ensured,but still a similar approach can be followed replacing gradients by subgradients. Actually, ∂F ( α ) = conv [ ¯ x ∈ arg min h α ( x ) ( f a (¯ x )) a ∈A . A decomposition approach for (MEP) . One of the most popular approaches for solvingthe planar 1-center problem is due to Elzinga and Hearn [16]. The Elzinga-Hearn approach(EH, for short) is based on constructing the minimum enclosing disks of three demand pointsuntil the whole sets of points is assured to be fully covered. At each iteration, full coverageis checked, and the triplet of points is changed by incorporating the demand point which isfurthest from the actual covering disk replacing one of the three points. The overall worst casecomplexity of the algorithm is O ( n ), but its performance in practice outperforms other existingstrategies. Although the EH algorithm was initially designed to solve the planar unweighted1-center problem (taking advantage of the geometry of disks), the approach can be adequatelyextended to the higher-dimensional weighted case [20]. In what follows we describe how the EHdecomposition paradigm can be applied to solve (MEP).First, let us assume that r ∗ ≥ X ∗ , is convex. Let us denote by: P a U ,ω ( r ∗ ) = { z ∈ R d : X u ∈U ω u k a − u − z k ≤ r ∗ } , ∀ a ∈ A , the polyellipse with foci { a − u : u ∈ U} and radius r ∗ . Then, any optimal solution x ∗ ∈ T a ∈A P a U ,ω ( r ∗ ). Actually, since r ∗ is minimum: \ a ∈A P a U ,ω ( r ∗ ) = X ∗ since otherwise, it would contradict the optimality of r ∗ .The following result is the main component to apply the decomposition approach to (MEP). Theorem 3.
There exists S ⊆ A with | S | = d + 1 such that \ a ∈ S P a U ,ω ( r ∗ ) = X ∗ .Proof. Let A d +1 = { S ⊆ A : | S | = d + 1 } . We define Q r ∗ ( S ) = T a ∈ S P a U ,ω ( r ∗ ) for all S ∈ A d +1 .Since T a ∈A P a U ,ω ( r ∗ ) = X ∗ , in particular we have that X ∗ ⊆ Q r ∗ ( S ) = ∅ for all S ∈ A d +1 .On the other hand, since T a ∈A P a U ,ω ( r ∗ ) = ∅ , then, for all sets S ⊆ A with | S | = d + 1, T a ∈ S P a U ,ω ( r ∗ ) = ∅ . Suppose that for all S ∈ A d +1 , int ( Q r ∗ ( S )) includes X ∗ . Since int ( Q r ∗ ( S ))is convex, there exists ε > int ( Q r ∗ − ε ( S )) = ∅ . Thus, T a ∈ S P a U ,ω ( r ∗ − ε ) = ∅ , andby Helly’s theorem also T a ∈A P a U ,ω ( r ∗ − ε ) = ∅ , contradicting the optimality of r ∗ .Alternatively, if there exists S ∈ A d +1 such that X ∗ ⊂ Q r ∗ ( S ) and int ( Q r ∗ ( S )) = ∅ , it impliesthat X ∗ is included in one face of dimension at most d − P a U ,ω ( r ∗ ), for all a ∈ S . Then forany ε > T a ∈ S P a U ,ω ( r ∗ − ε ) = ∅ , and this implies that r ∗ is the minimum value for a fullcoverage of A and also that X ∗ is a set of optimal solutions. ✷ (cid:3) Based on the above result, a decomposition approach can be derived to solve the problem, asdescribed in Algorithm 1.
Theorem 4.
Algorithm 1 computes the minimum-radius polyellipsoid in O ( | A | d +2 ) time forstrictly convex norms. inimal Radius Enclosing Polyellipsoids 11 Algorithm 1
A decomposition algorithm for solving (MEP).Compute the minimum radius polyellipsoid covering S : r , x . Input: A and S ⊆ A with | S | = d + 1. k = 0. Let ρ k = max a ∈A X u ∈U ω u k a − u − x k k and a k ∈ A reaching such a maximum. ρ k = r k then STOP. else S k +1 = S k ∪ { a k }\{ b k } with b k reaching r k +1 = max b k ∈ S k max a ∈ S k ∪{ a k }\{ b k } X u ∈U ω u k a − u − x k k end ifOutput: r k and x k . Go to Proof.
By Theorem 3, an optimal polyellipsoid is defined by d + 1 demand points. Such apolyellipsoid will be determined by the minimum-radius polyellipsoid covering these points. Atany iteration k of the algorithm, one computes in the minimum-radius polyellipsoid for d + 1given points (which takes constant time provided that d is fixed, see Theorem 2). If the optimalpolyellipsoid cover all the points, then, the furthest points (at distance ρ k ) are those at sum of thedistances exactly r k . Then, we have constructed a polyellipsoid, that covers all the points in A and which is the one which minimally covers a subset of A . Otherwise, a new d + 1-points subsetof A is constructed by replacing one (the most convenient) of its elements by the furthest pointto the previous polyellipsoid. Observe that at each iteration the solution of the subproblem isunique because of the strict convexity of the norm. Then, since the number of subsets with d + 1elements is O ( |A| d +1 ) and the polyellipsoids generated at each iteration are monotone increasingin radius (the assumption r k +1 > r k is considered because otherwise, the algorithm stops) thealgorithm stops in at most O ( |A| d +1 ) iterations and each iteration requires O ( |A| ) (the numberof comprobations for susbtitution in Step 3). ✷ (cid:3) Observe that the main advantage of this approach is that one only solves minimum enclosingpolyellipsoid problems for d + 1 points, which in practice is very convenient, in particular whenthe number of demand points in A is large. Moreover, in the planar case, no matter the numberof demands points in the problem this approach only needs to solve MEP problems coveringthree of the demand points. Also, we will see that the number of iterations of the approach, asfor the 1-center problem, is small, reducing significatively the computation times as comparedwith a plain SOCP formulation.Note also, that in Algorithm 1, one has to solve (MEP) for subsets of A with d + 1 points. Inthe planar Euclidean case with a single focus, one can explicilty derive the the minimum enclosingdisk of three points [16], being this oracle doable in constant time. In the general case, one hasto solve a SOC-programming problem, increasing the theoretical complexity of the procedure.In case the norm is not strictly convex, the above procedure is not assured to reach the optimalsolution of the problem since the solution of the subproblems with d + 1 demand points may notbe unique. In the example of Figure 4, one has four demand points that want to be covered bya square with minimum edge length (i.e., the minimum enclosing polyellipse for ℓ norm and asingle foci). If at the first iteration of Algorithm 1 one chooses S = { a , a , a } , the optimalcenters for the minimal enclosing squares for those three points are drawn in the left picture(solid line). If the solver chooses x ∗ ∈ X ∗ as the optimal center center, clearly, it does not cover a , so in the next iteration, one element in S is replaced by a . In this case, the three possible inimal Radius Enclosing Polyellipsoids 12 a a a a X ∗ x ∗ a a a a x ∗ a a a a Figure 4.
Example of non-monotonicity in the radius for non strictly convexnorms on Algorithm 1.replacements result in the same square length. For instance, if a is replaced by a , one gets thatone of the minimum enclosing squares for S = { a , a , a } is the one drawn in the right picture,which has smaller radius than the one obtained in the previous iteration.To overcome the above undesired fact, we adapt Algorithm 1. In this modification, at the k thiteration, if r k +1 = max b ∈ S k max a ∈ S k ∪{ a k }\{ b k } X u ∈U ω u k a − u − x k k < r k , then, instead of removing b k from S k , we keep both a k and b k in S k +1 , increasing by one the cardinality of the set S k .Observe that this modification assures that the sequence of radii { r k } is monotone and thealgorithm converges in a finite number of steps. Nevertheless, the complexity of the algorithmis higher now since one may need to perform O ( |A| d +1 ) iterations but in each iteration a linearproblem with up to O ( |A| ) constraints could have to be solved.4. Computational Experiments
We have run a battery of experiments in order to show the performance of the Second Or-der Cone programming formulation and the decomposition approach (Algorithm 1) for solving(MEP). We used the classical planar dataset in [15] and also the datasets for location prob-lems recently proposed in [7]. The first instance, consists of the classical 50-points (
EWC ) ofEilon, Watson-Gandy and Christofides, while the rest of the instances consist of geographicalcoordinates of different population areas in Slovakia. The sizes of these instances are 4873(
Partiz´anske ), 9562 (
Koˇsice ), 79612 ( ˇZilina ) and 663203 (
Slovakia ). We compute the mini-mal radius enclosing polyellipse with number of foci ranging in { , , , } and several norms.In particular we use four strictly convex norms: ℓ , ℓ , ℓ and ℓ , and three block norms: ℓ , ℓ ∞ and hex , where hex is a block norm whose unit ball is a hexagon with extreme points {± (2 , , ± (1 , , ± ( − , } (see [26]). The foci were randomly chosen from the set of demandpoints. We analyze both the unweighted and weighted problems. For the latter, we use eitherrandom weights (for the EWC instance) or the weights provided in [7] for the rest of the instances.We implemented both the SOCP model and the decomposition approach in Python 3, usingGurobi 8 as the optimization solver. The experiments were run in a Mac OSX Mojave with anIntel Core i7 processor at 3.3 GHz and 16GB of RAM.In Tables 1-3 we report the results of our computational experience. Table 1 and Table 2show the CPU times required for solving the problems for ℓ p -norm and block norms polyellipses,respectively. We provide, for each instance its size ( |A| ), the norm used ( k · k ) and the numberof foci ( |U| ), the CPU running times (in seconds) for solving (MEP) both with the Second Or-der Cone programming formulation ( Time
SOC ) and with the decomposition approach (
Time
DEC ).As one can observe, except for the smallest instance ( |A| = 50), the decomposition approachoutperforms the results obtained with the SOC formulation. Observe that both approaches opti-mally find minimum enclosing polyellipsoids but, while the SOC formulation solves the problems inimal Radius Enclosing Polyellipsoids 13 handling simultaneously all the points in A (adding all the constraints enforcing covering all itselements), Algorithm 1 solves (iteratively) the problems considering only d + 1 points at a time(in the strictly convex norm case). This fact is hard to be observed while the number of demandpoints is small to medium but it clearly pops up as the number of points increases. The readercan also observe in Table 3 that the number of iterations needed to solve the problems with thedecomposition approach is rather small. This table shows the average number of iterations ofthe decomposition approach for solving the instances and also the average cardinality of the sets S k required for solving the block-norm cases. The number of iterations of the decompositionapproach for all the instances range in [2 , ℓ p -norm polyellipses, it means that at most6 (MEP) problems with 3 demand points had to be solved, to obtain the solution of the originalproblem. In the case block-norm case, it may happen that one has to increase the number ofpoints for which the (MEP) problem has to be solved (the sets S k ). However, the maximumnumber of points that we obtained in our experiments was 6, that combined with the maximumnumber of iterations, 6, gives rise to the highly competitive CPU times needed for solving theproblems.We have marked with OoM those instances for which the SOC approach is not able to load/solvingthe problem because of an
Out of Memory flag.We also show the number of iterations performed by our decomposition procedure ( it dec ),and for the block norms, the maximum cardinality of the sets S k used in the modification ofAlgorithm 1 for non strictly convex norms.5. Case study: One dimensional covering polyellipsoids
In this section we analyze problem (MEP) in case d = 1, and describe closed formulas for theminimal covering one-dimensional polyellipsoid covering a finite set of demand points. Observethat while the solution of the 1-center problem on the line can be trivially obtained ( x ∗ = (max a ∈A a + min a ∈A a ) and r = (max a ∈A a − min a ∈A a )), a further analysis is needed whenthe number of foci is greater than 1.Let us first analyze the shape of a polyellipsoid on the real line. Let U ⊂ R , ω ∈ R |U| + and r ≥
0. The one-dimensional polyellipse with foci U , weights ω and radius r is: P U ,ω ( r ) = { z ∈ R : X u ∈U ω u | z − u | ≤ r } . By Proposition 1, since polyellipsoids are closed and bounded sets, P U ,ω ( r ) is, either empty, a sin-gle point or an interval. On the other hand, the Weber problem on the line, i.e., min x ∈ R P u ∈U ω u | x − u | , is solved at the median interval, X ∗ = [ x ∗ , x ∗ f ]. Let us also denote by r ∗ = P u ∈U ω u | x ∗ − u | for x ∗ ∈ X ∗ and ¯ u = P u ∈U ω u u .The following result whose proof is straightforward provides the explicit shape of P U ,ω ( r ). Lemma 1.
There exist u , u f ∈ U such that P U ,ω ( r ) = ∅ r < r ∗ , (cid:2) u , u f (cid:3) if r = r ∗ , " r − ¯ u + 2 P u ∈U ω u u P u ∈U ω u − , r − ¯ u + 2 P u ∈U f ω u u P u ∈U f ω u − if r > r ∗ . where U = { u ∈ U : u ≤ u } and U f = { u ∈ U : u ≤ u f } .Proof. The first statement is straightforward. Let us analyze the case r ≥ r ∗ .Assume, w.l.o.g., that the elements in U are sorted in increasing order, i.e. U = { u < · · · < u k } , i n i m a l R a d i u s E n c l o s i n g P o l y e lli p s o i d s Unweighted Weighted ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ |A| |U| t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC
50 1 0.01 0.14 0.00 0.02 0.01 0.01 0.01 0.02 0.01 0.02 0.00 0.02 0.01 0.01 0.01 0.025 0.07 0.10 0.04 0.05 0.08 0.10 0.05 0.10 0.07 0.10 0.03 0.08 0.06 0.13 0.05 0.0710 0.14 0.14 0.04 0.14 0.14 0.24 0.20 0.20 0.18 0.14 0.03 0.09 0.16 0.15 0.19 0.0925 0.56 0.36 0.12 0.20 0.57 0.49 0.48 0.47 0.53 0.55 0.14 0.43 0.44 0.35 0.57 0.644873 1 1.88 0.12 0.62 0.15 2.80 0.16 3.09 0.13 1.74 0.14 0.58 0.17 2.46 0.13 2.51 0.175 11.99 0.59 5.29 1.08 13.59 0.90 11.80 0.77 11.78 0.55 3.79 0.51 11.57 0.76 12.72 0.7710 31.61 1.79 7.28 1.05 31.53 1.60 37.75 1.64 31.87 1.50 15.20 1.07 45.39 1.49 46.34 2.0325 117.96 2.77 23.30 3.34 112.68 3.88 129.25 3.53 148.39 2.60 39.00 2.45 167.49 3.48 147.35 3.839562 1 4.38 0.30 1.56 0.29 4.36 0.30 4.25 0.35 4.12 0.33 1.38 0.30 4.19 0.32 4.16 0.375 30.58 1.27 7.11 1.38 27.46 1.70 35.00 1.36 33.42 1.28 9.90 1.30 50.72 1.33 42.58 1.3810 97.25 2.82 16.65 1.80 71.98 2.04 84.53 2.67 205.41 3.19 33.87 2.47 245.38 2.67 160.70 3.3325 213.34 6.12 47.93 5.88 317.86 6.27 365.62 6.81 503.54 6.35 368.24 6.12 531.46 6.41 463.49 6.5579612 1 52.86 1.72 16.51 1.67 56.13 1.80 67.14 1.91 52.48 1.73 16.56 1.72 55.40 1.80 66.80 1.805 379.13 10.04 84.89 9.75 432.17 10.71 413.29 10.58 7200 12.77 377.48 7.72 1490.33 7.69 1376.67 7.8510 1427.21 29.41 186.11 19.63 1262.71 15.21 1637.17 16.29 4448.48 29.64 1040.49 15.45 2496.65 17.05 3479.00 17.4725 7200 71.20 389.17 58.66 7200 38.77 7200 40.07 7200 38.43 2734.22 37.40 7200 37.55 7200 37.64663203 1 1125.80 14.39 271.17 18.82 788.36 14.70 951.78 15.22 1149.75 13.57 281.24 17.97 806.99 14.09 949.27 14.415 7200.01 84.43 1785.31 64.11
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM T a b l e . C P U T i m e s f o r M EP w i t h s t r i c t l y c o n v e x n o r m s . i n i m a l R a d i u s E n c l o s i n g P o l y e lli p s o i d s Unweighted Weighted ℓ ℓ ∞ hex ℓ ℓ ∞ hex |A| |U| t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC t SOC t DEC
50 1 0.00 0.02 0.00 0.02 0.00 0.05 0.00 0.08 0.00 0.02 0.00 0.055 0.01 0.03 0.01 0.05 0.02 0.17 0.01 0.03 0.01 0.03 0.03 0.0910 0.02 0.09 0.01 0.09 0.04 0.18 0.02 0.08 0.01 0.06 0.05 0.1825 0.10 0.09 0.06 0.16 0.06 0.46 0.10 0.14 0.07 0.26 0.06 0.454873 1 0.17 0.16 0.08 0.19 0.21 0.52 0.18 0.14 0.10 0.18 0.22 0.515 0.92 0.61 0.86 0.44 1.44 3.13 0.93 0.73 0.95 0.71 1.33 2.4410 2.05 1.17 2.08 1.13 2.87 5.04 2.15 1.46 1.81 1.10 2.43 4.8325 7.06 3.61 5.64 2.79 8.63 9.10 6.50 2.11 5.59 3.08 8.78 9.069562 1 0.34 0.44 0.19 0.41 0.44 0.98 0.35 0.41 0.20 0.41 0.44 1.035 2.31 1.14 1.86 1.09 2.85 5.88 2.58 0.83 2.40 1.11 2.97 4.6910 5.62 3.34 4.17 2.10 6.45 7.11 3.01 2.78 2.83 2.76 4.21 11.6925 15.08 6.63 11.54 5.08 20.35 35.70 13.16 6.63 10.68 6.64 17.86 29.1979612 1 5.25 2.65 1.55 2.64 5.66 6.06 6.75 2.59 1.74 3.21 5.56 6.005 26.92 11.13 22.10 8.65 34.63 47.97 29.45 8.97 18.69 6.43 32.93 28.7510 68.77 13.15 43.87 12.78 94.67 57.86 71.63 17.19 40.54 12.35 88.22 56.7125 189.67 31.70 143.21 30.64 308.12 243.43 144.43 32.56 110.38 32.37 224.39 175.40663203 1 55.84 26.05 16.00 17.37 56.48 47.34 53.39 26.07 23.01 20.67 57.84 44.195 390.16 53.92 364.48 54.29 798.43 219.50 1232.96 69.97 533.16 52.13 1234.53 213.3310 3674.59 141.34
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM
OoM T a b l e . C P U T i m e s f o r M EP w i t hp o l y h e d r a l n o r m s . a ndd e n o t e b y u = − ∞ a nd u k + = ∞ , a nd ω j = ω u j , f o r j = ,..., k . F o r a n y z ∈ R ,l e t inimal Radius Enclosing Polyellipsoids 16 ℓ p norms block norms |A| |U| it DEC it DEC | S k |
50 1 2.8 4.7 4.75 3.8 3.5 3.010 3.4 3.5 3.225 3.9 3.2 3.04873 1 3.5 4.3 4.75 3.8 4.3 3.710 4 4.2 3.325 3.6 3.7 39562 1 4.0 5.3 4.75 4.1 4.0 3.010 4.0 4.7 3.725 4.0 5.0 3.079612 1 3 4.3 4.35 3.8 4 3.710 3.9 3.2 325 3.6 3.5 3.2663203 1 3.3 4.5 4.55 3.3 3.2 3.210 4.0 3.5 3.025 3.3 3.5 3.0
Table 3.
Average number of iterations in the decomposition approach for allthe instances. s = s ( z ) ∈ { , , . . . , k } such that u s < z ≤ u s +1 , then: X u ∈U ω u | z − u | = s X j =1 ω j ( z − u j ) − k X j = s +1 ω j ( z − u j )= ( s X j =1 ω j − k X j = s +1 ω j ) z − s X j =1 ω j u j + k X j = s +1 ω j u j = (2 s X j =1 ω j − z + ¯ u − s X j =1 ω j u j . The extremes of the nonempty closed interval P U ,ω ( r ) are those z ∈ R verifying that X u ∈U ω u | z − u | = r . Thus, z is an extreme of P U ,ω ( r ) if (2) is equal to r .(1) If s X j =1 ω j = 12 (= k X j = s +1 ω j ), then r = X u ∈U ω u | z − u | = ¯ u − s X j =1 ω j u j . In case |U| iseven, all the points in [ u s , u s +1 ] have the same sum of weighted distances to the foci.Otherwise, P U ,ω ( r ) = { u s } . Thus, u = u s and u f = u s +1 or u f = u s in the statementof the Lemma. Note that in this case r coincides with r ∗ and P U ,ω ( r ) is the polyellipsoidinduced by the ω -weighted median of the points. inimal Radius Enclosing Polyellipsoids 17 (2) If s X j =1 ω j = 12 , then z = r − ¯ u + 2 s X j =1 ω j u j (2 s X j =1 ω j −
1) . In this case, P U ,ω ( r ) is a proper closedinterval with nonempty interior, and then there exist two points defining the extremes ofthe interval, thus, there exist s and s f in with s < s f such that, u = u s and u f = u s f and for which the equation in the Lemma holds. ✷ (cid:3) Le us consider that a set of demand points
A ⊆ R is given. Finding the minimum enclosingone-dimensional polyellipsoid with foci U ⊆ R , consists of finding the minimum value of r suchthat all the elements in A belong to a x -translation of P U ,ω ( r ), for some x ∈ R . Since r determinesthe length of the interval, it is clear that P x + U ,ω ( r ) = [min a ∈A a, max a ∈A a ], since otherwise itwould not cover the points or it would not be minimal.Applying Lemma 1 to the x -translated polyellipse, we get that, after some algebra, the x -translated polyellipsoid P x + U ,ω ( r ) is in the form: x + r − ¯ u + 2 X u ∈U : x + u ≤ a ω u u X u ∈U : x + u ≤ a ω u − , x + r − ¯ u + 2 X u ∈U : x + u ≤ af ω u X u ∈U : x + u ≤ af ω u − Equaling the extremes of that interval to [min a ∈A a, max a ∈A a ] =: [ a , a f ] we get that: r ∗ = (1 − X u ∈U : x + u ≤ a ω u )( a f − a − X u ∈U : a There exists an unique pair ( s , s f ) ∈ { , . . . , q } × { s + 1 , . . . , k } such that ( x ∗ , r ∗ ) is the optimal center and radius of the minimal enclosing polyellipsoid.Proof. The proof follows by noting that the problem is equivalent to find x ∗ ∈ P a −U ,ω ( r ) ∩ P a f −U ,ω ( r ). Thus, if two different solutions exist, it contradicts the minimallity of r . Thus, x ∗ isunique, so there exists an unique combination of indices ( s , s f ) fullfiling the conditions. ✷ (cid:3) The above result allows us to terminate the search as soon as a solution x ∗ in the form aboveis found verifying that max { j ∈ { , . . . , k } : u j ≤ a − x } = s . Corollary 1. Provided that the sets A and U are sorted, the one-dimensional minimum enclosingpolyellipse can be found in O ( |U| ) time for any number of points. Also, if all the demand points fall inside the interval determined by some translation of thefoci, the optimal polyellipsoid can be explicitly stated. Corollary 2. If max u ∈U u − min u ∈U u ≤ max a ∈A a − min a ∈A a . The minimum enclosing polyel-lipse is determined by: r ∗ = max a ∈A a − min a ∈A a and x ∗ = max a ∈A a + min a ∈A a − X u ∈U ω u u Extensions This section addresses two extensions of (MEP) where we can still exploit the methods andtools previously developed. In particular, we analyze the problem that, apart from finding theposition and the radius of the covering polyellipsoid, incorporates the foci selection among agiven set of potential candidates. Also, we extend the notion of polyellipsoid to ordered medianpolyellipsoid, following the analogy between the Weber and Ordered Median Location problems.6.1. Foci Selection. In many situations, the foci, which are assumed to be fixed in (MEP),may be unknown. In the following we address the question of how to optimally select a givennumber of foci from a finite set of candidate points.Let A ⊆ R d be the set of demand points and B ⊆ R d a set of potential foci. The goal is toselect from B a subset of k foci U ⊆ B with |U| = k to adequately cover the set of demand pointsby the translated polyellipsoid with foci in U and minimal radius. The problem can be stated,in a mathematical programming manner, as follows:min x ∈ R d, U⊆B : |U| = k max a ∈A ϕ U a ( x ) (MEP-FS)Recall that ϕ U a = X u ∈U ω u k x − a − u k for all U ⊆ B and a ∈ A .Observe that once the foci, U , are selected, the problem reduces to (MEP). However, theenumeration of the (cid:0) |B| k (cid:1) combinations of possible foci becomes intractable in practice. Thus, wepropose, first, a mixed integer non linear programming formulation for the problem, and then,we develop a decomposition approach for solving the problem geometrically. inimal Radius Enclosing Polyellipsoids 19 In the mathematical programming formulation, apart from the variables used to reformulate(MEP), we use the following set of binary decision variables: y u = (cid:26) u is selected as a focus0 otherwise. , ∀ u ∈ B . The following model, (MEPFS), allows us to determine the optimal foci and the polyellipsoid.min r (2)s.t. X u ∈B y u = k, (3) d au ≥ k a − u − x k y u , ∀ a ∈ A , u ∈ B , (4) r ≥ X u ∈B ω ua d ua , ∀ a ∈ A , (5) r, d ua ≥ , x ∈ R d ,y u ∈ { , } . The objective function (2), as in (MEP), minimizes the radius, r , by choosing an adequate subsetof k foci from B (3). The way the distances are accounted in the problem depends on the selectedfoci: in constraint (4), the distance between a demand point a ∈ A and a translated foci u + x with u ∈ B is either k a − u − x k if u is chosen from B (i.e., if y u = 1) or 0 otherwise. Finally,as in (MEP), constraint (5) ensures that r is defined as the maximum among all the sum of thedistances from each demand point to the translated selected foci.The formulation above corresponds to a mixed integer non-linear program. The discretecharacter comes from the y -variables, whereas the non-linearity appears by the constraints (4),which can be rewritten using big- M constants as follows: d au + M au (1 − y u ) ≥ k a − u − x k , ∀ a ∈ A , u ∈ B , where M au is an upper bound on the value of the norm k a − u − x k . As already mentioned inSection 3.2, the optimal x -value coincides with the solution of a weighted Weber problem for theset of the demand points { a − u : u ∈ U} for some a ∈ A , thus, these big M -constants can bederived explicitly.The above formulation reduces to a Mixed Integer Second Order Cone Optimization (MIS-OCO) problem provided that the considered norm is polyhedral or in the ℓ p ( p ≥ 1) family.Then, medium size instances can be solved with nowadays available off-the-shell software.In what follows, we describe an adaptation of the Elzinga-Hearn based approach described inAlgorithm 1 to this new case where the foci have to be selected from B . Recall that the successof the EH-based approach comes from decomposing the original problem into smaller ones, bysolving the minimal covering polyellipsoid problem on a reduced subset of demand points, inparticular on subsets of demand points S ⊆ A with cardinality d + 1. Following a similarscheme, for each subset S ⊂ A with | S | = d + 1, one has to compute not only the polyellipsoidbut also the optimal k foci from B . Let us denote by U R ( S ) the optimal set of foci from B thatminimally cover the points in S , with the condition that the points in R ⊆ B are not in U R ( S ),i.e., the solution of the following problem: r S,R := min U⊆B\{ R } : |U| = k max a ∈ S ϕ U a ( x ) (6)Observe that the problem above is a reduced version of (MEP-FS) in which the set of potentialfoci is B\{ R } and the set of demand points is S . Our decomposition approach for solving(MEP-FS) consists of solving problems in the form of (6) , sequentially for different sets S , until inimal Radius Enclosing Polyellipsoids 20 a termination covering criterion is met. The pseudocode of this method is described in Algorithm2. First, we initialize the set S to a set of d + 1 demands points from A and R = ∅ . At eachiteration, it , we compute the set of optimal foci (solving (6)) with demand points in the actualset, S it , excluding those in the set R . Apart from the set of optimal foci for those points, theproblem gives us a lower bound of the optimal radius for the covering polyellipsoid, r S it ,R , whichis updated. Since the original problem is always feasible, the restricted problem (6) is feasible ifand only if one can choose k elements from B\ R , i.e., if |B\ R | ≥ k . If this were not possible, theprocess is finished. Once a set of optimal foci, U it +1 , is computed, one has to solve (MEP) forthose foci using Algorithm 1 which also gives the d + 1 points defining the covering polyellipsoid(the set S it +1 in the pseudocode). The solution of that problem provides an upper bound onthe optimal value of (MEP-FS) ( U B it ) which is updated at each iteration. The procedure isrepeated until the upper and lower bounds coincide or the number of reduced potential foci isless than k . Observe that Algorithm 2 terminates in a finite number of iterations since B and Algorithm 2 Decomposition approach for solving (MEP-FS). Input: A , S ⊆ A with | S | = d + 1, it = 0 R = ∅ , U B = ∞ , LB = 0. while U B > LB and |B\ R | ≥ k do U it +1 = U R it ( S it ), LB = max { LB, r S it ,R } .Solve (MEP) for U = U it +1 and set S it +1 and U B = min { U B, U B it } ). R = R ∪ U it +1 .it → it + 1. end whileOutput: U = U it and r = U B . A are finite sets. Note also that two situations may lead to the termination of the algorithm.On the one hand, if at some iteration U B ≤ LB , then, it would imply that a feasible solution of(MEP) has been found (with radius U B ) and that for a subset S it +1 ⊆ A , any choice among thenon explored potential foci, gets a covering polyellipsoid with a larger radius than one alreadycomputed. Thus, the best solution found at that iteration must be optimal. On the otherhand, if at some iteration |B\ R | < k , then all the possible foci have already been explored, andthe solution must be among those already computed. Thus, the algorithm outputs the optimalsolution of (MEP-FS), as desired.To conclude, we would like to add some comments on the complexity of Algorithm 2. At eachiteration, for solving (MEP-FS), one has to solve a MISOCO which in general could be NP-hard.Thus, although (MEP) can be solved in O ( |A| d +1 )-time, the overall complexity of Algorithm 2is in general non polynomial.We have run some experiments in order to analyze the performance of Algorithm 2 comparedto the MISOCO formulation. We have used two of the datasets used in Section 4, the one with50 points from [15] and the one with 4873 points from [7]. We use unweighted instances andnorms in { ℓ , ℓ , ℓ , ℓ } . For each dataset we have randomly chosen the set of potential foci, B ,from the demand points with sizes 10, 15 and 20. We consider the foci to choose, k , ranging in { , } , and we set a time limit of 30 minutes in Gurobi for solving the MISOCO problems.We report in Table 4 the results of our experiments. We show, for each instance the CPU timesdevoted to the MISOCO formulation (t MISOCO ) and for Algorithm 2 (t DEC ). For the latter, wealso provide the number of iterations needed to solve the problems, i.e. the number of MISOCOproblems with three demand points ( It ) as well as the average number of iterations of Algorithm inimal Radius Enclosing Polyellipsoids 21 it DEC ). W also report the MINLP gap obtained withthe MISOCO formulation at the time limit (% Gap MISOCO ).As one can observe from the results, Algorithm 2 outperforms the MISOCO formulation in allthe instances. Moreover, in most of the instances for the 4873-dataset, the solver was not ableto find the optimal solution within the time limit with the MISOCO formulation. We observethat loading the MISOCO problem in Gurobi is highly time consuming, particularly when non-Euclidean or polyhedral norms are used because the large number of auxiliary variables, andSOC and linear constraints that have to be introduced to represent the norms. This can beseen in columns LinCtrs , SOCCtrs and BinVars , that indicate the number of linear,SOC constraints, and binary variables of the entire MISOCO problem, respectively. We havehighlighted, with TL ∗ , the instances for which the MISOCO problem was not able to be loadedin Gurobi within 4 hours (in those instances the gap is not available).6.2. Ordered Median Polyellipsoids. As mentioned in the introduction, polyellipsoids areidentified with the levels curves of the classical Weber (median) problem. Several extensionsof this problem has been analyzed in the Location Analysis literature. One of them is the so-called ordered median continuous location problem [4, 18, 26]. In this problem, apart from thedemand points, U = { u , . . . , u k } ∈ R d , a distance measure induced by a norm k · k and weights ω , . . . , ω k ∈ R + , one is also given a set of modelling weights λ , . . . , λ k ∈ R and the goal is tofind the optimal placement x minimizing the following aggregating function of the distances: F λ U ,ω ( x ) = k X i =1 λ i ω σ ( i ) k x − u σ ( i ) k , where σ is a permutation of the indices { , . . . , k } such that ω σ ( i ) k x − u σ ( i ) k ≥ ω σ ( i +1) k x − u σ ( i +1) k for all i = 1 , . . . , k − 1. The ordered median function generalizes most of the objective functionsconsidered in facility location. For instance, choosing λ = (1 , . . . , λ = (1 , , . . . , 0) one gets themax criterion which is used in the center problem. It is known that F λ U ,ω is convex if and onlyif the λ -weights verify λ ≥ · · · ≥ λ k . Thus, we will assume such a condition on the weights.In what follows we draw some comments on the extension of the notion of polyellipsoid to thecase in which they are defined as the level curves of functions in the form of F λ U ,ω .Let U = { u , . . . , u k } ⊆ R d be a set of foci and ω , . . . , ω k ∈ R + a set of weights for the foci.Let us also consider a set of weights λ , . . . , λ k ∈ R + such that λ ≥ · · · ≥ λ k , the λ -orderedmedian polyellipsoid with radius r is defined as: E λ U ,ω ( r ) = { z ∈ R d : k X i =1 λ i ω σ ( i ) k z − u σ ( i ) k = r } We denote by P λ U ,ω ( r ) the region bounded by E λ U ,ω ( r ). Since F λ U ,ω is convex, then P λ U ,ω ( r ) is anonempty bounded and convex set for r ≥ min x ∈ R d F λ U ,ω ( x ). In Figure 5 we show four differentchoices for the lambda weights and the resulting ordered median polyellipses for the same set ofthree foci on the plane.Analogously to (MEP), the ordered median minimal enclosing polyellipsoid problem consists offinding the minimum radius, r , and the translation of the ordered median polyellipsoid P λ U ,ω ( r )such that all the demand points are covered by the obtained polyellipsoid. A mathematicalprogramming formulation for the problem is the following:min x ∈ R d ,r ∈ R + max a ∈A k X i =1 λ i ω σ a ( i ) k a − u σ a ( i ) − x k (MEP λ OM ) inimal Radius Enclosing Polyellipsoids 22 |A| |B| k Norm t MISOCO t DEC it it DEC % Gap LinCtrs SOCCtrs BinVars 50 10 5 ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ ℓ TL ℓ ℓ TL ℓ ℓ TL ℓ TL ℓ TL ℓ ℓ TL ℓ ℓ TL ℓ TL ℓ TL ℓ TL ℓ TL ℓ TL ℓ TL ∗ ℓ TL ∗ ℓ TL ∗ ℓ TL ℓ TL ∗ ℓ TL ℓ TL ∗ ℓ TL ℓ TL ∗ ℓ TL ℓ TL ∗ Table 4. Results of the Experiments for the Foci Selection Problem. inimal Radius Enclosing Polyellipsoids 23 λ = (1 , , λ = (1 , , λ = (1 , , λ = (1 , . , . Figure 5. Shapes of ordered median polyellipses for the same foci U = { ( − , , (2 , , ( − , } and different λ -weights.where for all a ∈ A , σ a is a permutation of the indices { , . . . , k } such that ω σ ( i ) k a − u aσ ( i ) − x k ≥ ω σ ( i +1) k a − u aσ ( i +1) − x k for all i = 1 , . . . , k − 1. The problem can be reformulated, using theresults proved in [4] as:min r s.t. r ≥ k X i =1 ( v ai + t ai ) , ∀ a ∈ A ,v ai + t aj ≥ λ j ω i k a − u i − x k , ∀ a ∈ A , i, j = 1 , . . . , k,r ≥ ,v ai , t ai ∈ R , ∀ a ∈ A , i = 1 , . . . , k. where the auxiliary variables v and t allow one to represent the sorting factor in the problemwithout using binary variables. As in (MEP), in case the norm is polyhedral or in the ℓ p familywith p ≥ 1, the problem above can be reformulated as a second order cone programming problem.A further geometrical analysis for this problem is in order. For the sake of presentation,let us assume that the considered norm is strictly convex. The optimal solution of the problem(MEP λ OM ) is given by the smallest r such that ∩ a ∈A P λ U− a,ω ( r ) = ∅ . Nevertheless, the geometricalstructure of the sets P λ U− a,ω ( r ) is intricate since the evaluation depends on the relative sortingof the points with respect to the translated foci. In order to apply an “a la” Elzinga-Hearnalgorithm to solve this problem we have to introduce the concept of ordered region (see [29]).An ordered region O σ is a closed, connected region of R d such that for all z ∈ O σ a permutationthat sorts the vector ( ω k a − u − z k , . . . , ω k k a − u k − z k ) in non-increasing sequence is σ . Thatis, the chain of inequalities ω σ k a − u σ − z k ≥ ω σ k a − u σ − z k ≥ . . . ≥ ω σ k k a − u σ k − z k isvalid for all z ∈ O σ . It is straightforward to observe that a subdivision, S a , of R d into orderedregions with respect to the point a ∈ A can be obtained with the following arrangement offunctions: ω i k a − u i − z k = ω j k a − u j − z k , ∀ i < j ∈ { , . . . , k } . Denote by S the intersectionof all the subdivisions S a , a ∈ A . In each cell of the subdivision S the order of the vectors inimal Radius Enclosing Polyellipsoids 24 ( ω k a − u − z k , . . . , ω k k a − u k − z k ) , ∀ a ∈ A is constant and thus, P λ U− a,ω ( r ) for all a ∈ A are standard polyellipsoids. In order to have a valid application of a decomposotion algorithmfor problem (MEP λ OM ), similar to Algorithm 1, one would have to replicate it on each cell ofthe subdivision S . This would imply to solve constrained (to the cell) (MEP) until an optimalrestricted solution on that cell is found. Then, choosing the minimum value among all for thedifferent cells will result in the optimal solution for (MEP λ OM ).To have an idea of the complexity of these subdivisions we present the analysis for the intuitive,important case of the unweighted Euclidean norm. In this case, for each a ∈ A the elements thatinduce the arrangement S a are k a − u i − z k = k a − u j − z k , ∀ i < j ∈ { , . . . , k } . It is well-known that for each i = j these equations define Euclidean bisectors (hyperplanes) and thereforethe number of ordered regions is the number of cells in an arrangement of O ( |A| ) hyperplanesin dimension d . This number is known to be O ( |A| d ). Intersecting these subdivision for all a ∈ A gives rise to the following upper bound, O (( |A| ) d ), on the number of elements in thesubdivision S . The discussion above allows us to state the following results, whose proof followsfrom the application of Theorem 4, with complexity O ( |A| d +2 ), on each one of the O ( |A| d ) fulldimensional cells of the subdivision S . Theorem 5. Let us assume that the dimension d is fixed. The unweighted version of Problem (MEP λ OM ) with Euclidean norm can be solved in polynomial time O ( |A| d +2 ) . Conclusions In this paper we analyze minimum radius covering problems with shapes based on polyellip-soids. We study the problem of covering a finite set of d -dimensional demand points by optimallyfinding the dilation and translation of a polyellipsoid with given foci. We provide several for-mulations for the problem in general dimension d : a primal SOC formulation, a conic dual SOCmodel and a Lagrangean dual program whose resolution needs to solve different Weber problems.We analyze the worst case complexity of these problems and develop an Elzinga-Hearn decom-position approach to solve them geometrically. We report the results of an extensive batteryof experiments to show the performance of the proposed approaches. We further analyze theone-dimensional problem and derive a linear time algorithm (in the number of foci) for solvingthe problem as well as some closed formulas for the covering polyellipse in terms of the foci andthe demand points. Finally we draw some comments on two extensions of the problem: 1) theproblem of simultaneously finding the minimum radius and selecting the foci from a potentialfinite set of candidates; and 2) the use of more general shapes, which we call ordered medianpolyellipsoids, to cover the demand points. As far as we know, both extensions, have been firstintroduced and analyzed here.Further research on the topic includes, among others, the analysis of the multiple facility case,in which several polyellipsoids have to be located by minimizing the maximum sum of the closestdistances from the demand points to the foci, extending the multifacility center problems; ormaximal covering location problems using polyellipsoidal shapes. References [1] Andretta, M., Birgin, E.G.: Deterministic and stochastic global optimization techniques for planar coveringwith ellipses problems. European Journal of Operational Research (1), 23–40 (2013)[2] Berman, O., Drezner, Z., Krass, D.: Generalized coverage: New developments in covering location models.Computers & Operations Research (10), 1675–1687 (2010)[3] Blanco, V., ElHaj-BenAli, S., Puerto, J.: Minimizing ordered weighted averaging of rational functions withapplications to continuous location. Computers & Operations Research (5), 1448–1460 (2013)[4] Blanco, V., Puerto, J., ElHaj-BenAli, S.: Revisiting several problems and algorithms in continuous locationwith ℓ τ norms. Computational Optimization and Applications (3), 563–595 (2014)[5] Boyd, S., Vandenberghe, L.: Convex optimization. Cambridge university press (2004) inimal Radius Enclosing Polyellipsoids 25 [6] Canbolat, M.S., von Massow, M.: Planar maximal covering with ellipses. Computers & Industrial Engineering (1), 201–208 (2009)[7] Cebecauer, M., Buzna, L.: Large-scale test data set for location problems. Data in brief , 267–274 (2018)[8] Chazelle, B.M., Lee, D.T.: On a circle placement problem. Computing (1-2), 1–16 (1986)[9] Chen, Y., Ye, X.: Projection onto a simplex. arXiv preprint arXiv:1101.6081 (2011)[10] Church, R.L.: The planar maximal covering location problem. Journal of Regional Science (2), 185–201(1984)[11] Drezner, Z.: Note—on a modified one-center model. Management Science (7), 848–851 (1981)[12] Duchi, J., Shalev-Shwartz, S., Singer, Y., Chandra, T.: Efficient projections onto the l 1-ball for learningin high dimensions. In: Proceedings of the 25th international conference on Machine learning, pp. 272–279.ACM (2008)[13] Dyer, M.: On a multidimensional search technique and its application to the euclidean one-centre problem.SIAM Journal on Computing (3), 725–738 (1986)[14] Dyer, M.: A class of convex programs with applications to computational geometry. In: Proceedings of theEighth Annual Symposium on Computational Geometry, SCG ’92, pp. 9–15. ACM, New York, NY, USA(1992)[15] Eilon, S., Watson-Gandy, C., Christofides, N.: Distribution management: mathematical modelling andpractical analysis. NEw York, Hafner (1971)[16] Elzinga, J., Hearn, D.W.: Geometrical solutions for some minimax location problems. Transportation Science (4), 379–394 (1972)[17] Erd¨os, P., Vincze, I.: On the approximation of convex, closed plane curves by multifocal ellipses. Journal ofApplied Probability , 89–96 (1982)[18] Espejo, I., Rodr´ıguez-Ch´ıa, A.M., Valero, C.: Convex ordered median problem with ℓ p -norms. Computers &Operations Research (7), 2250–2262 (2009)[19] Hansen, P., Perreur, J., Thisse, J.F.: Location theory, dominance, and convexity: Some further results.Operations Research (5), 1241–1250 (1980)[20] Hearn, D.W., Vijay, J.: Efficient algorithms for the (weighted) minimum circle problem. Operations Research (4), 777–795 (1982)[21] Helly, E.: ¨Uber systeme von abgeschlossenen mengen mit gemeinschaftlichen punkten. Monatshefte f¨ur Math-ematik (1), 281–302 (1930)[22] Maxwell, C., et al.: 1. on the description of oval curves, and those having a plurality of foci. Proceedings ofthe Royal Society of Edinburgh , 89–91 (1851)[23] Megiddo, N.: The weighted euclidean 1-center problem. Mathematics of Operations Research (4), 498–504(1983)[24] Melzak, Z., Forsyth, J.: Polyconics. i. polyellipses and optimization. Quarterly of Applied Mathematics (2),239–255 (1977)[25] Nesterov, Y., Nemirovskii, A.: Interior-Point Polynomial Algorithms in Convex Programming. Society forIndustrial and Applied Mathematics (1994)[26] Nickel, S., Puerto, J.: Location theory: a unified approach. Springer Science & Business Media (2006)[27] Nie, J., Parrilo, P.A., Sturmfels, B.: Semidefinite representation of the k-ellipse. In: Algorithms in algebraicgeometry, pp. 117–132. Springer (2008)[28] Plastria, F.: Continuous covering location problems. Facility location: applications and theory , 37–79(2002)[29] Puerto, J., Fern´andez, F.: Geometrical properties of the symmetric single facility location problem. Journalof Nonlinear and Convez Analysis , 321–242 (2000)[30] Renegar, J.: A faster pspace algorithm for deciding the existential theory of the reals. Tech. rep., CornellUniversity Operations Research and Industrial Engineering (1988)[31] Tax, D.M., Duin, R.P.: Support vector data description. Machine learning (1), 45–66 (2004)[32] Vincze, C., Kov´acs, Z., Csorv´assy, Z.F.: On the generalization of erd¨os-vincze’s theorem about the approxi-mation of closed convex plane curves by polyellipses. arXiv preprint arXiv:1705.04318 (2017)[33] Wang, W., Carreira-Perpin´an, M.A.: Projection onto the probability simplex: An efficient algorithm with asimple proof, and an application. arXiv preprint arXiv:1309.1541 (2013)[34] Ward, J.E., Wendell, R.E.: Using block norms for location modeling. Operations Research (5), 1074–1090(1985)[35] Weber, A.: Theory of the Location of Industries. University of Chicago Press (1929)[36] Weiszfeld, E.: Sur le point pour lequel la somme des distances de n points donn´es est minimum. TohokuMathematical Journal , 355–386 (1937) inimal Radius Enclosing Polyellipsoids 26 IEMath-GR, Universidad de Granada, SPAIN. E-mail address : [email protected] IMUS, Universidad de Sevilla, SPAIN. E-mail address ::