## Abstract

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*.

## 1 Introduction

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.

## 2 Terminology

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

*k*-link shortest path formulation. The dominated region can be partitioned into certain rectangles. Let

*A*, , be the rectangle defined by the subregion of , which is exclusively dominated by all points in and no other point in

_{ij}*N*. An example of this partition is given in Figure 1. For every such rectangle, we define

*w*as the area of rectangle

_{ij}*A*. If we define and , the rectangle

_{ij}*A*can be written as Hence, we get and we can calculate all the weights

_{ij}*w*, , in time.

_{ij}## 4 An Integer Programming Formulation

*w*the area of

_{ij}*A*, , where

_{ij}*A*is the rec-tangle defined by the subregion of , which is exclusively dominated by and no other point in

_{ij}*N*. The following IP formulation models the corresponding HSSP: Thereby, variable is equal to 1 if and only if is selected, and variable

*x*determines whether the subregion

_{ij}*A*is covered by some point in , which is guaranteed by the constraints (3). Constraint (2) ensures the compliance of the selection of exactly

_{ij}*k*points, and the objective function (1) calculates the value of the current hypervolume indicator, which has to be maximized.

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.

*s*are

_{ij}*surplus variables*(Bazaraa et al., 2011).

*x*and

_{ij}*s*, , , according to the ordering of the constraints (6), the structure of the constraint matrix corresponding to the constraints (5) and (6) of (

_{ij}*LP*) is given by where

_{k}*C*is a -matrix and is the negative of the identity matrix. Let us denote by the submatrix , where is the vector of all 1s, and by

*D*the submatrix . Observe that

*D*is totally unimodular and has the

*consecutive 1s property*(Nemhauser and Wolsey, 1999) and thus is also totally unimodular.

The constraint matrix of (*LP _{k}*) is totally unimodular.

Let *B* denote an arbitrary squared submatrix of the constraint matrix of (*LP _{k}*).

*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 *j*th column (*Laplace expansion*). Since this column has only one nonzero entry, say *b _{ij}*, we get , where

*M*is the minor of matrix

_{ij}*B*formed by eliminating row

*i*and column

*j*from

*B*. The minor

*M*also corresponds to a squared submatrix of the constraint matrix. If we follow the above Laplace expansion, after

_{ij}*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 (*LP _{k}*) is integral. In particular, (

*IP*) can be solved by solving its linear relaxation (

_{k}*LP*) with an appropriate LP solver.

_{k}## 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).

*N*. The graph construction is based on the observation that for each choice of a subset , for , the contribution to the hypervolume indicator of the consecutive points for two indices with only depends on the coordinates of the points and . For each element , we create a node . In addition, we also add two other nodes 0 and to

*V*, as source and target nodes, respectively. We add the arcs , for all with to

*E*. According to the notation in the preprocessing step (see Section 3), the cost

*c*of an arc

_{uv}*e*is defined as follows: where for all .

_{uv}By definition, we get that the cost *c _{uv}* 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 *s*_{1} to , and ends in the node . Since the cost of a used arc *e _{uv}* 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 *z ^{h}* 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.

*c*on the fly. By setting and , we get for Hence, we only need to precompute all costs and for all , which can be done in time (see Algorithm 2).

_{uv}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 .

^{10}, we immediately get This leads together with expression (9) to a contradiction and we get .

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.

Hence, incorporating the idea from Remark ^{11} and the complexity for the sorting mentioned at the end of Section 2, we can state the following result.

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 .

## 7 Conclusion

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 .

*IP*) from Section 4 can be extended to an arbitrary number of dimensions . In this case, any appropriate IP solver can be used to obtain an exact solution, even if the corresponding LP formulation will not, in general, define an integral polyhedron, as in the following example: Here, the solution of the linear programming relaxation (w.r.t. ) has objective value , whereas the solution of the integer program has objective value . Nevertheless, the linear programming relaxation can still be used to obtain an upper bound on the optimal hypervolume indicator.

_{k}*k*-link shortest path problem. This can be observed in the following example. Let us consider three points: Looking at the corresponding dominated regions/boxes (w.r.t. ), one can infer that each pair out of the three induced boxes possesses a nonempty intersection belonging only to the two boxes considered. Hence, there is no unique sorting of the points as in the two-dimensional case. Nevertheless, suppose that we assign an arbitrary order to the three nodes in the corresponding digraph with five nodes (including source and target nodes). Clearly, to model the subset selection corresponding to all points except one, the arc jumping over this node must have cost equal to the exclusive volume of the corresponding point. Then, the path corresponding to the subset including only the midpoint (w.r.t. the digraph nodes) will have the wrong weight, since it contains only the two arcs that jump over the second and second to last nodes, respectively. We would miss subtracting the part of the exclusive volume of the two points not selected that is dominated by both. On the other hand, with two objectives, the application of the obtained formulations and algorithms to other indicators such as, for instance, the generalized hypervolume indicator introduced by Grunert da Fonseca and Fonseca (2012), should follow naturally.

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).

## Acknowledgments

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}.

## References

*k*-link path in graphs with the concave Monge property and applications

## Note

^{1}

Code available at http://docs.theinf.uni-jena.de/code/ssp.zip