The hypervolume subset selection problem consists of finding a subset, with a given cardinality k, of a set of nondominated points that maximizes the hypervolume indicator. This problem arises in selection procedures of evolutionary algorithms for multiobjective optimization, for which practically efficient algorithms are required. In this article, two new formulations are provided for the two-dimensional variant of this problem. The first is a (linear) integer programming formulation that can be solved by solving its linear programming relaxation. The second formulation is a k-link shortest path formulation on a special digraph with the Monge property that can be solved by dynamic programming in time. This improves upon the result of in Bader (2009), and slightly improves upon the result of in Bringmann et al. (2014b), which was developed independently from this work using different techniques. Numerical results are shown for several values of n and k.
Given a set of nondominated points in objective space, the hypervolume indicator measures the dominated region of this space bounded by some reference point (Zitzler and Thiele, 1998). Because of its properties (Knowles and Corne, 2002; Zitzler et al., 2003; 2008), this indicator has been instrumental both in the assessment of the performance of multiobjective evolutionary algorithms (Knowles et al., 2000; Zitzler et al., 2000) and in the development of multiobjective selection and archiving procedures for such algorithms (Fleischer, 2003; Knowles et al., 2003; Knowles and Corne, 2003; Zitzler and Künzli, 2004; Beume et al., 2007; Bader and Zitzler, 2011).
The hypervolume subset selection problem (HSSP) consists of finding a subset of k elements from a set of n nondominated points that maximizes the hypervolume indicator. This has been considered computationally expensive for an arbitrary number of objectives (Bader and Zitzler, 2011). Therefore, the main focus has been on finding the n−k elements that contribute the least in terms of hypervolume in a greedy manner, e.g., in Huband et al. (2003); Beume et al. (2007); Igel et al. (2007); see also approximation results by Bringmann and Friedrich (2009) and an asymptotically optimal algorithm for finding all contributions of every element in the given set in for two and three dimensions (Emmerich and Fonseca, 2011).
So far, the tightest time complexity bound for solving the HSSP to optimality for dimensions and arbitrary k is (Bringmann and Friedrich, 2010). For the particular case of , to the best of our knowledge, only two approaches have been proposed. Bader (2009) introduced a dynamic programming algorithm with time complexity. This algorithm is based on the fact that the contribution to the hypervolume indicator of the leftmost point of a given nondominated subset depends only on its immediate neighbor. Independently from the present work, Bringmann et al. (2014b) proposed another approach to the same problem with a time complexity of . This algorithm computes, in the th iteration, all maximal hypervolume indicator values using at most points with respect to n appropriately chosen reference points. This can be done for each reference point by computing the maximum of different linear function evaluations. The running time is achieved by using a linear time algorithm to compute the upper envelope of lines.
In this article, we propose two different formulations for the two-dimensional case of the HSSP: an integer programming formulation and a k-link shortest path formulation. Both formulations are based on a preprocessing step, which makes a partition of the dominated region into different areas induced by the set of nondominated points. We show that the polyhedron of the linear programming relaxation of the first formulation is integral, which allows linear programming methods, such as simplex and interior-point methods, to solve the HSSP. In the k-link shortest path formulation, the arc costs have a special property, called the Monge property. This property allows us to solve the HSSP with a simple dynamic programming approach in time, which slightly improves upon the result of Bringmann et al. (2014b). This time complexity suggests that our approach should perform better for large values of k.
The remainder of this article is organized as follows. In Section 2 we introduce concepts, definitions, and notation used throughout this article. In Section 3 we explain the crucial preprocessing step for calculating the weights used in the integer programming formulation and for the introduction of the arc costs in the k-link shortest path formulation. In Section 4 we introduce an integer programming formulation for the HSSP and prove the integrality of the polyhedron of its linear programming relaxation. In Section 5 we present the k-link shortest path formulation, which is used to achieve an improved complexity bound for the HSSP. An experimental analysis is provided in Section 6. Finally, in Section 7 we provide some conclusions and avenues for future research.
A point dominates if . Let denote a set of nondominated points, in which no point in N is dominating another point in N.
The hypervolume indicator maps a set of nondominated points to the size of the region in the corresponding space dominated by this set and bounded below by a reference point.
In the following sections, we assume a set of nondominated points , some reference point , and a desired cardinality to be given. We further assume that the points in set N are sorted in increasing order of the first component, i.e., for , which can be achieved in time.
3 Preprocessing: Decomposition of the Dominated Region
4 An Integer Programming Formulation
For proving that this IP can be efficiently solved with the help of the corresponding linear programming (LP) relaxation, we need the notion of a totally unimodular matrix, which will lead to a famous result from integer optimization.
A matrix is called totally unimodular if the determinant of each square submatrix of A belongs to .
Let b, , . If is totally unimodular and , then P is an integral polyhedron, i.e., each of its nonempty faces contains an integral point.
The constraint matrix of (LPk) is totally unimodular.
Let B denote an arbitrary squared submatrix of the constraint matrix of (LPk).
Case 1:B is completely contained in or completely contained in D and therefore, since both matrices are totally unimodular.
Case 2:B possesses and columns from matrix and matrix D, respectively, w.l.o.g. no duplicate column from D.
Choose some column from B belonging to D, and expand the determinant of B with respect to the jth column (Laplace expansion). Since this column has only one nonzero entry, say bij, we get , where Mij is the minor of matrix B formed by eliminating row i and column j from B. The minor Mij also corresponds to a squared submatrix of the constraint matrix. If we follow the above Laplace expansion, after t steps, we get a submatrix of B matching Case 1, i.e., . Then, by construction, we get .
Since B was an arbitrarily chosen squared submatrix, we have shown the totally unimodular property of the constraint matrix.
The polyhedron corresponding to (LPk) is integral. In particular, (IPk) can be solved by solving its linear relaxation (LPk) with an appropriate LP solver.
5 A k-link Shortest Path Formulation with the Monge Property
In the following, we show that the HSSP can be modeled using a k-link shortest path formulation in a directed graph (digraph). The corresponding shortest path problem with a cardinality constraint can then be solved using a dynamic programming (DP) approach. This digraph has a special structure, the Monge property, that allows us to solve the HSSP in time. Without loss of generality, we assume , since for the case , only the preprocessing step is needed to obtain the optimal objective value (see also Algorithm 2).
By definition, we get that the cost cuv describes the contribution to the hypervolume indicator of the whole set , which is the exclusive volume of the set and denoted by . An example for the graph construction is depicted in Figure 2.
From the construction of the graph, we immediately get the following property.
Each choice in the HSSP of a subset of N with cardinality corresponds to a path in our constructed digraph with exactly k arcs that starts in node 0, visits the nodes s1 to , and ends in the node . Since the cost of a used arc euv corresponds to the exclusive volume of the jumped-over nodes to , the hypervolume contribution minus the total cost of the path corresponds then to the hypervolume contribution of the corresponding subset of N. Hence, the k-link shortest path problem on our constructed digraph models the HSSP with desired cardinality .
Since in our special k-link shortest path problem the Bellman principle of optimality is valid, we can use a straightforward DP approach to solve this problem (see Algorithm 1). In each iteration, the length of the optimal path for the problem of finding the -link shortest path from 0 to v is calculated. However, this would not lead directly to a much better running time than Bader’s DP algorithm, since finding the minimum in line 4 is in a naïve way done in , resulting in an overall running time in .
In the following, we show that the time complexity can be improved by proving some special structure for this digraph, the so-called (concave) Monge property (Aggarwal et al., 1994).
For , we define as the area of the rectangle induced by the two corner points zh and the special reference point with and (compare Figure 4).
Adapting a proof of Aggarwal and Park (1989), we can state the following equivalent property.
Going back to Algorithm 1, to find for a fixed the new entries , , we have to find in a matrix , where only the entries , , are relevant, for each column v the minimal value, which is then assigned to . Ignoring all irrelevant columns and rows, matrix can be considered a square matrix in . For an efficient searching algorithm that finds these minimal values for each column, we first need to define the notion of a totally monotone matrix.
Consider a matrix . For a column , let denote the index of the greatest row containing the minimum value of column j. Matrix A is called monotone if implies . Moreover, matrix A is called totally monotone if each submatrix is monotone.
Let , , denote a totally monotone matrix. Then the matrix-searching algorithm of Aggarwal et al. (1987) finds the minimal entries in all columns in time.
The matrix-searching algorithm considers recursively submatrices while reducing in each iteration, with the help of clever comparisons of selected entries, the currently considered submatrix to a square matrix where irrelevant rows are deleted.
is totally monotone for a fixed .
Choose some arbitrary submatrix with rows , w.l.o.g. without empty columns.
Choose two columns from the submatrix . Suppose now that and let and .
Using the matrix-searching algorithm from Theorem 13 in the DP approach (see Algorithm 1), the k-link shortest path problem in the constructed digraph can be solved in time.
The two-dimensional HSSP can be solved in time.
Note that a simple modification of the algorithm induced from Theorem 16 can obtain the solution of the HSSP for all in time.
6 Experimental Analysis
The algorithm from Section 5 has been implemented and is available at https://eden.dei.uc.pt/∼paquete/HSSP. We ran our approach on randomly generated instances in order to study the influence of n and k on its performance. For comparison purposes, we also ran the implementation described in Bringmann et al. (2014b).1 Both programs were written in Java and used the same data structures as much as possible. The instances were generated according to the description given in Bringmann et al. (2014b). Each instance was run three times before measuring the running time in order to eliminate fluctuations. We considered and . For each pair , we ran each implementation on 10 different random instances. All experiments were conducted on a computer with Intel Pentium processor, running at 2.0 GHz with 3 GB RAM and 64-bit Linux (ArchLinux, kernel 3.18.5). Both implementations were compiled with OpenJDK version 1.7.0_75.
The plots of Figure 6 show the computational results obtained with both approaches for different values of n and fixed ratio , in a log-log scale. The plots indicate that the running time of both approaches increases in a very similar way as the instance size grows, which is expected because of their similar time complexity. Noteworthy, our approach performs better as k becomes closer to n, as suggested by the time complexity of our approach; see plots (c) and (d). This observation becomes clearer in Figure 7, for different values of k and fixed n, which shows that our approach performs better for .
In this article, we considered the two-dimensional hypervolume subset selection problem. Based on the investigation of structural properties of the problem, we proposed a new integer programming formulation, and showed that the polyhedron of its linear programming relaxation is integral. Moreover, we developed a k-link shortest path formulation on a simple, directed acyclic graph. Exploiting the special structure of the arc costs, we stated a dynamic programming approach which solves the HSSP in .
Finally, we remark that the hypervolume indicator can also be used as a measure of the quality of discrete representations of an optimal set of nondominated points (see Sayn, 2000; Ruzika and Wiecek, 2005). Given a large nondominated set in the biobjective case, the algorithm described in this article can be used to find an optimal cardinality-constrained representation with respect to the hypervolume indicator very efficiently (see also Bringmann et al., 2014a).
The authors gratefully acknowledge support under the bilateral cooperation research project “Tractability in multiobjective combinatorial optimization,” funded by the Deutscher Akademischer Austausch Dienst and Fundação para a Ciência e a Tecnologia. Authors Kuhn and Ruzika acknowledge the support of the Federal Ministry of Education and Research Germany, grant DSS_Evac_Logistic, FKZ 13N12229. Fonseca and Paquete acknowledge the iCIS project (CENTRO-07-ST24-FEDER-002003), which is co-financed by QREN, in the scope of the Mais Centro Program and European Union’s FEDER. Fonseca, Kuhn, and Paquete also acknowledge the Lorentz Center Workshop on Set-Oriented and Indicator-Based Multi-Criteria Optimization (SIMCO 2013), where some of the results presented in this work were developed. Finally, the authors acknowledge Karl Bringmann for pointing out the result in Bringmann et al. (2014b) as well as the technique mentioned in Remark 11.