Anne Auger
email: anne.auger@inria.fr

Johannes Bader

Dimo Brockhoff

Eckart Zitzler

Articulating User Preferences in Many-Objective Problems by Sampling the Weighted Hypervolume

Keywords: I.2.8 [Artificial Intelligence]: Problem Solving, Control Methods, and Search-Heuristic Hypervolume indicator, preference articulation, Monte Carlo sampling

The hypervolume indicator has become popular in recent years both for performance assessment and to guide the search of evolutionary multiobjective optimizers. Two critical research topics can be emphasized with respect to hypervolume-based search: (i) the hypervolume indicator inherently introduces a specific preference and the question is how arbitrary user preferences can be incorporated; (ii) the exact calculation of the hypervolume indicator is expensive and efficient approaches to tackle many-objective problems are needed. In two previous studies, we addressed both issues independently: a study proposed the weighted hypervolume indicator with which user-defined preferences can be articulated; other studies exist that propose to estimate the hypervolume indicator by Monte-Carlo sampling. Here, we combine these two approaches for the first time and extend them, i.e., we present an approach of sampling the weighted hypervolume to incorporate user-defined preferences into the search for problems with many objectives. In particular, we propose weight distribution functions to stress extreme solutions and to define preferred regions of the objective space in terms of so-called preference points; sampling them allows to tackle problems with many objectives. Experiments on several test functions with up to 25 objectives show the usefulness of the approach in terms of decision making and search.

INTRODUCTION

The hypervolume indicator has been proposed about a decade ago to compare the performance of multiobjective optimization algorithms [START_REF] Zitzler | Multiobjective Optimization Using Evolutionary Algorithms -A Comparative Case Study[END_REF]. It measures the quality of a set of solutions quantitatively as the "size of the space covered" [START_REF] Zitzler | Multiobjective Optimization Using Evolutionary Algorithms -A Comparative Case Study[END_REF]. Equivalent definitions based on the Lebesgue measure [START_REF] Laumanns | Approximating the Pareto Set: Concepts, Diversity Issues, and Performance Assessment[END_REF], or the attainment function [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] have been proposed later. The hypervolume indicator has also gained interest as a selection criterion within multiobjective evolutionary algorithms in recent years [START_REF] Huband | An Evolution Strategy with Probabilistic Mutation for Multi-Objective Optimisation[END_REF][START_REF] Nicolini | A Two-Level Evolutionary Approach to Multi-criterion Optimization of Water Supply Systems[END_REF][START_REF] Beume | SMS-EMOA: Multiobjective Selection Based on Dominated Hypervolume[END_REF][START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF][START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF], especially if the number of objectives is large [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF]. The main reason for the hypervolume's popularity is its property of refining the Pareto dominance relation: whenever a solution set is dominating another, the hypervolume indicator value of the former is strictly larger than the one of the latter [START_REF] Zitzler | SPAM: Set Preference Algorithm for Multiobjective Optimization[END_REF]-resulting in the fact that a set that optimizes the hypervolume indicator covers the entire Pareto front [START_REF] Fleischer | The measure of Pareto optima. Applications to multi-objective metaheuristics[END_REF]. Up to now, no unary indicator is known that is not based on the hypervolume and that has the refinement property.

Clearly, using the hypervolume indicator for guiding search introduces a certain preference. We investigated this recently by characterizing the distribution of µ points that maximize the hypervolume indicator [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF]. One question that is of special interest in practice is whether the inherent preference of the hypervolume indicator can be changed to arbitrary user preferences, e.g., towards extreme solutions or towards so-called preference points 1 . There already exist various studies that focus on the issue of preference articulation in EMO, in particular integrating preferences such as priorities, goals, and preference points [START_REF] Branke | Guidance in Evolutionary Multi-Objective Optimization[END_REF][START_REF] Deb | Reference Point Based Multi-Objective Optimization Using Evolutionary Algorithms[END_REF][START_REF] Rachmawati | Preference Incorporation in Multi-objective Evolutionary Algorithms: A Survey[END_REF]. However, none of these methods leads to a Pareto-compliant relation-to obtain such relations we proposed a generalized version of the hypervolume indicator where arbitrary user preferences can be incorporated by stressing certain regions of the objective space in terms of a weight distribution function [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF].

In the same paper [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF], we also propose three exemplary weight distribution functions to incorporate user preferences towards the extremes of the front and towards predefined preference points for biobjective problems. The study is only preliminary and inhibits two drawbacks: (i) the proposed weight distribution function for articulating preference points is not easily extendable to more than two objectives and (ii) the exact computation of the hypervolume indicator is expensive if the number of objectives is high 2 .

In this paper, we tackle these two drawbacks by estimating the weighted hypervolume indicator by Monte Carlo sampling and propose according weight distribution functions that can be used for incorporating user preferences for an arbitrary number of objectives. In particular, we

• introduce the approach of sampling the weighted hypervolume indicator based on the work for the standard hypervolume indicator in [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF] to avoid the exponential running time of the hypervolume indicator, • propose two weight distributions to articulate preferences towards extremes and towards predefined preference points which can be arbitrarily combined and applied to problems with any number of objectives and • show the potentials of the new approach experimentally for several test problems with up to 25 objectives by means of both visual inspection and statistical tests. The paper is organized as follows. We start with a brief review of the weighted hypervolume indicator of [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] (Sec. 2). Section 3 presents the general ideas of sampling the weighted hypervolume after which two examples of weight distributions and their combination are introduced in Sec. 4. An extensive experimental validation is presented in Sec. 5 before Sec. 6 concludes the paper.

THE WEIGHTED HYPERVOLUME

Throughout this study, we assume without loss of generality, that k objective functions fi :

X → R (1 ≤ i ≤ k)
need to be minimized that map solutions x from the decision space X to their objective vector f (x) = (f1(x), . . . , f k (x)). Furthermore, we assume the underlying dominance structure is given by the weak Pareto dominance relation where a solution a ∈ X weakly dominates another solution b ∈ X, written a b, if and only if fi(a) ≤ fi(b) holds for all 1 ≤ i ≤ k. We call a solution x * ∈ X Pareto optimal if there is no other x ∈ X that dominates x * , i.e., ∄x ∈ X : x x * ∧ x * x. The set of all Pareto optimal solutions is called Pareto set and its image Pareto front.

In this paper, we consider the following optimization goal: find a set A of µ compromise solutions such that the weighted hypervolume indicator value I w H (A) of [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] is maximized. To this end, we recall the definition of the weighted hypervolume indicator I w H (A) for a solution set A ⊆ X from [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF]:

I w H (A) = Z r l w( z) • αA( z)d z (1)
where

αA( z) =  1 if f (A) { z} 0 else is the attainment func- tion 3 of a solution set A ⊆ X, z ∈ R k and w : R k → R>0
is a weight distribution function where without loss of generality 4 R r l w( z)d z = 1; the limits l ∈ R k and r ∈ R k of the integral are parameters of the hypervolume, where usually l is set to (-∞, . . . , -∞) and r is called the reference point of the indicator. We will see in the following section, that these 2 Meanwhile, the #P-hardness proof in [START_REF] Bringmann | Approximating the Volume of Unions and Intersections of High-Dimensional Geometric Objects[END_REF] has theoretically shown that the hypervolume computation is indeed exponential in the number of objectives unless P = N P. 3 Where is generalized to sets as follows: A B iff ∀b ∈ B : ∃a ∈ A : a b 4 If R r l w( z)d z = 1 we can normalize w by R r l w( z)d z without changing the ordering of solutions with respect to the hypervolume indicator. points do not need to be fixed when the weighted hypervolume is sampled according to Gaussian distributions where the integral over the entire objective space is estimated.

SAMPLING THE HYPERVOLUME INDICATOR

Monte Carlo sampling is a well-known and easy-to-use approach to solve problems numerically by using random numbers. Monte Carlo sampling is used within several application areas such as atomic physics or finance. However, its most popular field of application is the computation of integrals [START_REF] Caflisch | Monte Carlo and Quasi-Monte Carlo Methods[END_REF].

Sampling the Non-Weighted Hypervolume

Using Monte Carlo methods to evaluate the hypervolume indicator is not new. Everson et al. [START_REF] Everson | Full Elite-Sets for Multiobjective Optimisation[END_REF] sampled the standard hypervolume for performance assessment whereas Bader and Zitzler [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF] proposed to sample the hypervolume within multiobjective search. The basic idea behind sampling the hypervolume is to draw N independent random samples X1, . . . , XN from the objective space uniformly at random within a defined (hyper-)box and use the fraction of dominated objective vectors as an estimate for the hypervolume indicator value IH(A). The attainment function αA of a solution set A ⊆ X formalizes the domination of objective vectors by assigning 1 to an objective vector if it is dominated by A and 0 otherwise. The unweighted hypervolume indicator can, thus, be estimated as5

IH(A) ≈ 1 N N X k=1 αA(X k ) . (2)
As indicated in [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF], the space from which the objective vectors are sampled should be as small as possible to avoid unnecessary samples, e.g., the (hyper-)box could have the hypervolume's reference point as one corner and the ideal point as the opposite one. When used during the selection step of a hypervolume-based multiobjective evolutionary algorithms, e.g., in SMS-EMOA [START_REF] Beume | SMS-EMOA: Multiobjective Selection Based on Dominated Hypervolume[END_REF] or MO-CMA-ES [START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF], the boxes can be much smaller if only the hypervolume loss of a single solution has to be estimated, see [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF] for details.

Sampling the Weighted Hypervolume

The straightforward way of approximating the weighted hypervolume indicator of Eq. 1 is to sample X1, . . . , XN uniformly at random as before and use

I w H (A) ≈ 1 N N X k=1 w(X k )αA(X k )
as an estimate of the weighted hypervolume. In this approach, however, the precision of the estimation heavily depends on the weight distribution w: if the support of w is small, the number of samples N needs to be large to have a reliable estimation. Using Hoeffding's inequality [START_REF] Hoeffding | Probability Inequalities for Sums of Bounded Random Variables[END_REF], one can show that the length of a confidence interval for a given confidence level is proportional to the supremum of w. In the extreme case of a dirac "function" as suggested in [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] this would result in an infinite length for the confidence interval.

Therefore, we propose a different approach to sample the weighted hypervolume indicator here. Since w is positive and R r l w( z)d z = 1, the weight function w is a density function. In principle, any density function can be used as w. For an efficient way of sampling, however, we choose w in the following such that we are able to sample efficiently random numbers distributed according to w. For this reason, multivariate normal distributions and exponential distributions will be used for sampling non-uniformly.

To give the explicit expression of the Monte Carlo estimator let X w denote a random variable admitting w as probability density function 6 . Let X w 1 , . . . , X w N be N independent samples of random variables distributed as X w . The weighted hypervolume I w H (A) corresponds to the expected value of the random variable αA(X w ) and can be approximated by

I w H (A) ≈ 1 N N X k=1 αA(X w k ) . (3)
This technique of sampling according to the weight distribution function instead of uniformly has the advantage that the accuracy of the estimate, i.e., the confidence interval, is independent of the weight distribution 7 .

In order to sample weight distributions that are defined as a mixture of several independent distributions wi (1 ≤ i ≤ m) as proposed in the following section, we propose to distribute the number of samples among the different distributions in the following way: a weight distribution w( z) = P m i=1 pi • wi( z) with P 1≤i≤m pi = 1 is estimated by sampling each of the distributions wi independently with N • pi samples and summing up the resulting estimates.

INTEGRATING USER PREFERENCES

In this section, we present two different weight distribution functions to express user preferences. Both distributions are continuous probability densities that enable to draw samples according to the procedure presented above. The first distribution allows to attribute importance to one objective and the second to emphasize a preference point in the objective space. We complete the section by demonstrating how any number of the two distributions can be combined, e.g., to use more than one preference point.

Stressing the Extremes

One potential preference a user might have is to optimize preferably one objective, say fs. The corresponding weight distribution should therefore increase for decreasing values of fs. In terms of the rest of the objectives, the weight distribution should stay constant for changing values in order 6 We refer to [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF] for an extensive overview how random samples from those distributions can be generated. 7 Denoting the estimator of I w H (A) (Eq. 3) as c

I w H (A, N ) = 1 N P N k=1 αA(X w k ), Hoeffding's inequality implies that with probability larger than 1 -α, I w H (A) belongs to [ c I w H (A, N ) - tN,α, c I w H (A, N ) + tN,α]
where tN,α = `8 N log(2/α) ´1/2 which is independent of w and the same confidence interval than for the non-weighted case. In other words, it is not more expensive to do a Monte Carlo integration of the weighted hypervolume than for the standard hypervolume indicator.

not to introduce additional preference.

In [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] the authors proposed to use an exponential function as the weight distribution. Here, the same distribution is represented by the probability density function whose marginal distribution for objective fs is an exponential distribution with rate parameter λ and whose marginal distributions of the remaining objectives is a uniform distribution:

w(z1,• • •, z k ) = 8 < : " Q i =s (b u i -b l i ) " -1 λe -λ(zs-b l s ) z ∈ B 0 z / ∈ B where B = [b l 1 , b u 1 ] × . . . × [b l k , b u k ]
denotes the space with non-zero probability density.

Figure 1(a) shows the weight distribution for a biobjective problem when stressing f1 with an exponential distribution in f1 (λ = 5) together with a uniform distribution in the interval [0, 0.95] in the second objective

(B = [b l s , b u s ] × [b l 2 , b u 2 ] = [0, ∞] × [0, 0.95]
). The spread of the distribution is inversely proportional to the parameter λ. Hence, the smaller λ the steeper the weight distribution increases at the border of the objective space and the smaller the weight farther away (see Figure 2(a) for contour plots of the exponential weight distribution for distinct values of λ).

Preference Points

Another user preference is the preference point [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. This point, as well as to a lesser extent the adjacent region, represents the most important part of the objective space for the user. Together with the location of the preference point, denoted by µ = (µ1,

• • • , µ k ) T ∈ R k , the user has to define a direction t = (t1, • • • , t k ) T ∈ R k .
The solutions should preferably lie along this direction if the preference point cannot be reached or, on the contrary, even better solutions are found. The corresponding weight distribution function reflects this preference by having the largest values at the preference point and along t while decreasing fast perpendicular to t. To this end, [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] proposes a bivariate ridge-like function that cannot be easily extended to an arbitrary number of objectives. Therefore, we propose to use the following multivariate Gaussian distribution here, which allows an efficient sampling according to Eq. 3 and which can be used for any number of objectives. Besides µ and t, let σε, σt ∈ R denote standard deviations of the distribution. Then the following probability density function describes a multivariate normal distribution centered at µ

w( z) = 1 (2π) k/2 |C| 1/2 e -1 2 ( z-µ) T C -1 ( z-µ))
where the covariance matrix C := σ 2 ε I + σ 2 t t t T / t 2 is nonsingular with orthogonal eigenvectors t, t2, • • • , t k where the vectors t2, . . . , t k can be taken from an arbitrary orthogonal basis of the hyperplane orthogonal to t. The eigenvalues associated to t, t2,

• • • , t k are σ 2 ε +σ 2 t , σ 2 ε , • • • , σ 2 ε ; |C| denotes the determinant of C.
The equidensity contours of the distributions are ellipsoids whose principal axis are t, t2, • • • , t k , see Fig. 1(b). The lengths of the axes are given by the two standard deviations (i) σt for the axis spanned by t and (ii) σε for the remaining k -1 axes perpendicular to t. The larger σt is chosen the farther the objective vectors can lie from the preference point in direction of ± t while they are still affected by the The second variance, σε, influences the extension of points close to the preference point. The smaller σε, the less widespread the solutions are (see Figures 2(b) for contour plots of three different choices of σε).

Combinations

A mixture of m weight distributions admitting the probability density functions w1( z), . . . , wm( z) yields the distribution with density

w( z) = p1w1( z) + . . . + pmwm( z)
where the pi are positive real numbers with p1 +. . .+pm = 1. Though it is not possible to translate any user preference directly to a weight distribution function as in [START_REF] Auger | Investigating and Exploiting the Bias of the Weighted Hypervolume to Articulate User Preferences[END_REF], a wide range of different user preferences can be represented by combining weight distributions. These are-in contrast to the weight distributions in [START_REF] Auger | Investigating and Exploiting the Bias of the Weighted Hypervolume to Articulate User Preferences[END_REF]-also applicable to problems with more than two objectives. In the next section we will examine mixtures of the two distributions presented above.

EXPERIMENTAL VALIDATION

In order to test the approach of articulating user preferences presented in Sec. 4, the sampling strategy is employed by an existing hypervolume-based algorithm called "Hypervolume Estimation Algorithm for Multiobjective Optimization (HypE)" [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF]. The application to different multiobjective test problems investigates three important aspects of the approach.

First, we investigate visually the influence of the different parameters on the distribution of the resulting Pareto front approximations for both approaches preferring preference points and extremes. In particular, we investigate for a preference point its location µ, the direction t and the influence of the standard deviations σε and σt. When stressing extremes, we show the effects of changing the parameter λ.

Secondly, the weighted hypervolume approach is visually compared to existing reference algorithms that do not optimize any user preference explicitly for problems with more than two objectives. This demonstrates that our approach is-in contrast to [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF]-also applicable to problems involving more than two objectives.

Finally, a short statistical comparison on problems with up to 25 objectives is carried out to investigate whether the generated Pareto front approximations obtained by HypE, as a matter of fact, better fulfill the underlying user preference than Pareto fronts resulting from reference algorithms.

Experimental Setup

To maximize the weighted hypervolume, the evolutionary algorithm HypE [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF] is used. While mating selection is performed uniformly, the next population is chosen in environmental selection as follows: first, the offspring and parent population are merged and allocated to fronts by nondominated sorting. Thereafter, the worst front is removed entirely as long as the resulting population size does not fall below the desired population size. If excluding entire fronts is no longer possible and there are still too many individuals in the population, the algorithm then starts removing individuals one-by-one from the worst front. Instead of other hypervolume-based selection schemes such as the ones in the SMS-EMOA [START_REF] Beume | SMS-EMOA: Multiobjective Selection Based on Dominated Hypervolume[END_REF] or the MO-CMA-ES [START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF], HypE not only considers the hypervolume solely dominated by a solution as its fitness but estimates the expected loss if more than one solution is deleted from the population by means of Monte Carlo sampling. To this end, 10, 000 samples are generated according to the probability density functions presented in Sec. 4 using the corresponding built-in functions of MATLAB R version 2008a. These samples are then used to calculate a fitness value for each individual of the worst front. Iteratively, the individual with the smallest value is deleted until the desired population size is reached where the fitness values are recalculated after each removal using the same samples as in the previous step. A detailed description of the fitness calculation of HypE can be found in [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF].

The evolutionary multiobjective algorithms NSGA-II [START_REF] Deb | A Fast Elitist Non-Dominated Sorting Genetic Algorithm for Multi-Objective Optimization: NSGA-II[END_REF] and IBEA [START_REF] Zitzler | Indicator-Based Selection in Multiobjective Search[END_REF] serve as reference algorithms. For the latter, the ε-indicator has been used since preliminary experiments showed this variant to be superior to the one using the hypervolume indicator. The parameters of IBEA are set as κ = 0.05 and ρ = 1.1. All algorithms are run for 100 generations. New individuals are generated by the SBX crossover operator with ηc = 15 and the variable-wise polynomial mutation operator with ηm = 20 [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF]. The crossover and mutation probabilities are set to 1 and 1/20 respectively. For the biobjective test problems both the population size and the number of offspring are set to 25 while for more objectives these numbers are doubled. For the biobjective investigations, the test problems ZDT1 (convex Pareto front), ZDT3 (discontinuous Pareto front) [START_REF] Zitzler | Comparison of Multiobjective Evolutionary Algorithms: Empirical Results[END_REF] and DTLZ2 [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF] (concave Pareto front) are utilized with 20 decision variables. For higher dimensions only DTLZ2 is employed.

(a) (b) (c) (d)

Visual Inspection of Parameter Choices

In this section, the influence of different parameters on the weight distribution functions and the resulting Pareto front approximations are investigated. Unless noted otherwise, we use σt = 0.5, σε = 0.05 and t = (1, 1) when stressing a preference point and B = [0, ∞] × [0, 3] when stressing the first objective (fs = f1). The weight distributions are indicated by contour lines at the intervals of 10% of the maximum value that arises. The contour lines do not reflect the actual weight but only the relative distribution thereof. We tested multiple runs for each testcase that led to similar results such that we mostly display only one run to illustrate the influence of the weight on the distribution of points.

Spread of the Distribution

Both proposed weight distribution functions have parameters that cause the weight to be more or less spread.

Direction of the Distribution

By t, the user can define the desired trade-off between the objectives for the case that either the preference point cannot be reached or if solutions dominating the preference point are obtained. In Fig. 2(c) the preference point is positioned at µ = (0, 0.4) which lies below the Pareto front and can therefore not be reached. In this case, the direction t predetermines where the resulting points lie. In the topmost example, a choice of t = (cos(80 • ), sin(80 • )) reflects a higher preference of the first objective at the expense of the second. On the other hand, the bottom figure is obtained for t = (1, 0), i.e., values of 0.4 are preferred for the second objective and only increases of the first objective are desired. The figure in the middle presents a compromise where the focus lies close to the diagonal, t = (cos(40 • ), sin(40 • )).

Location of the Preference Point

Since the preference point can be placed both too optimistically (as in the previous section) or too pessimistically, the parameter σt allows to tune how far away the individuals can be from the preference point and still be influenced by it. For a fixed σt however, the location of the preference point has a high impact on the resulting distribution of solutions, see Fig. 2(d). If none to only a few samples are dominated by the individuals (top, µ = (-1.2, -1.4)), no pressure towards the preference point is active-in fact only non-dominated sorting operates. In this case, the preference point should be combined with a uniform distribution, e.g., as in Fig. 3(a) 90% of the samples are used for the preference point and 10% to sample uniformly in the objective space within [0, 3] × [0, 3]. This causes the solutions to be distributed according to the unweighted hypervolume indicator as long as the preference point has no influence.

As soon as a couple of samples are dominated, the corresponding individuals are promoted which leads to an accumulation in that area (middle, µ = (-0.3, -0.5)). If the preference point is chosen very pessimistically (bottom, µ = (1.5, 1.3)) individuals are able to dominate all or most of the samples even if they are not located where the direction t dictates. This leads to a much ampler arrangement of solutions than expected considering the chosen σε.

Combinations of Distributions

As demonstrated in Sec. 4, any number of weight distribution functions can be combined as a weighted sum, even assigning them different weights or focus. For example, the user might define different preference points he or she is interested in as depicted in Fig. 3(b): three preference points are positioned at µ = (0.2, 0.8), at µ = (0.5, 0.5) and at µ = (0.8, 0.2). The one in the middle is declared to be the most important one by assigning the largest weight p2 = 0.5, the preference points to the left and right use p1 = 0.2 and p3 = 0.3 respectively. As expected, in this case the points are partitioned into disjoint regions around the three preference points. 10 individuals cluster around the center where the most samples emerge, 7 are associated with the preference point on the left and 8 with the one on the right.

To promote individuals at the border of the objective space, two exponential weight distributions can be added up as in Fig. 3(c) where λ = 10 with p1 = 0.3 for the first objective and p2 = 0.7 for the second.

Comparison Between Different Problems

In addition to ZDT1, the tests of the previous sections were also carried out for other test problems, namely ZDT3 which has a discontinuous Pareto front shape, DTLZ2 and ZDT6 (both non-convex). These three test problems are much harder to optimize and neither HypE nor the reference algorithms used were able to find Pareto optimal solutions. The points are nevertheless clustered at regions with the largest weight, see Figures 3(d), 3(e) and 3(f) where one preference point with µ = (0.7, 0.3) and σε = 0.1 is used.

High-Dimensional Spaces

For illustrative reasons in the previous section we applied the sampling procedure to biobjective problems. The advantage of the method, however, is that an arbitrary number of objectives can be tackled. Figure 4 DTLZ2 problem with 3 objectives. While NSGA-II and IBEA do not optimize any user defined preference, HypE uses two preference points at µ1 = (0.8, 0.2, 0.6) (p1 = 0.2) and µ2 = (0.2, 0.9, 0.5) (p2 = 0.8) with σε = 0.1 (shown as ellipsoids). This leads to a cluster of points at each preference point. The Pareto front approximation on DTLZ2 with 7 objectives is depicted in Fig. 5 by means of parallel coordinates plots for IBEA and HypE with σε = 0.05. The plot for NSGA-II is omitted due to space limitations; it can be noted that it looks similar to the one of IBEA except that NSGA-II does not come as close to the Pareto front as IBEA, i.e., the objective values are spread between 0 and 4.5. Both IBEA and NSGA-II generate solutions at the boundary of the objective space while only the former finds solutions near the Pareto front. To get solutions near the center of the Pareto front, HypE is applied with a preference point at 0.3780 • (1, . . . , 1). A second preference point is set at a random location near the Pareto front8 . The spread σε is set for both preference points to 0.05 and the probabilities of the mixture are set as 0.8 and 0.2 respectively leading to a population of solutions grouped around the two preference points (Fig. 5).

To investigate further whether HypE actually optimizes the weight distribution used during search, two versions of HypE are run in another experiment: the first version (denoted by HypEuni) uses a uniform weight distribution wuni with reference point r = (3, . . . , 3), the second (HypEp) optimizes a weight distribution wp consisting of a combination of two preference points at µ1 = (0.9, 0.1, 0.9, 0.1, . . . ) with p1 = 0.3, σε = 0.05 and at µ2 = (0.1, 0.9, 0.1, 0.9, . . . ) with p2 = 0.7, σε = 0.2. Both versions together with IBEA and NSGA-II then optimized the DTLZ2 test problem with different numbers of objectives ranging from 2 to 25. The Pareto front approximations for 50 independent runs of both versions of HypE as well as of NSGA-II and IBEA are then compared in terms of the weighted hypervolume indicator with the weight distribution function wp, see Table 1 where the hypervolume values are normalized to the minimum and maximum value observed for the respective number of objectives. In each case, HypEp outperforms statistically significant the other algorithms in terms of the hypervolume indicator with wp-assessed by Kruskal-Wallis and the Conover-Inman post hoc tests with a significance level of 5% according to [START_REF] Conover | Practical Nonparametric Statistics[END_REF]. This indicates that applying the weighted integration technique during search will generate Pareto front approximations that score better on the corresponding hypervolume indicator than using general purpose algorithms with no user defined search direction.

CONCLUSIONS

This paper has described two procedures to approximate the weighted hypervolume integration developed in [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] by means of Monte Carlo sampling. This offers the possibility to incorporate user preferences in a Pareto-compliant way, where problems with an arbitrary number of objectives can be considered. Two types of user preferences have been expressed by probability density functions that ease the fast generation of samples-one stressing certain objectives and the second emphasizing a preference point. Additionally, any combination of the two is possible.

The combination of the sampling approach with HypE, an existing hypervolume based algorithm, has been applied to various test problems. It has turned out by both visual inspection and statistical tests that the generated Pareto front approximations reflect the underlying weight distribution better than methods with no user defined preference.

Figure 1 :

 1 Figure 1: Illustrates the exponential distribution that corresponds to stressing one objective (a) and the Gaussian distribution representing a preference point (b). The left parts of the two subplots indicate the notation used along with a contour plot at intervals of 10% of the maximum value observed (which occurs on the second axis and at µ respectively). The right parts of the subplots show the actual value of the distribution as a third component.

Figure 2 :

 2 Figure 2: Shows the Pareto front approximations (dots) found by HypE using different weight distribution functions, shown as contour lines at intervals of 10% of the maximum weight value. For each column of three figures one parameter of the sampled distribution was modified, i.e., (a) the rate parameter of the exponential distribution λ, (b) the spread σε, (c) the direction t and (d) the location µ of the preference point (see text for details and the values used). The test problem is ZDT1 where the Pareto front is shown as a solid line.

  Figure 2(a) shows the weight distribution and the resulting Pareto front approximation using the exponential distribution proposed in Sec. 4.1 for λ = 100 (top), λ = 20 (center) and λ = 5 (bottom).

  Figure 2(b) shows the distribution of points for a preference point located at µ = (0.7, 0.3) where σε is set to 0.2 (top), 0.05 (middle) and 0.01 (bottom).

Figure 3 :

 3 Figure 3: These figures use the same visual elements as Fig. 2 which explains them in its caption. In (a) the same preference point is used as in the upper plot of Fig. 2(d) but spending 10% of the samples on a uniform distribution. Figure (b) shows the combination of three preference points, and Fig. (c) stressing both the first and second objective. The figures on the right show the distribution of the objective vectors when applying the same preference point to different test problems, i.e., (d) ZDT3 (only the positive part shown) (e) DTLZ2 and (f ) ZDT6.

Figure 4 :

 4 Figure 4: Pareto front approximations of five runs (depicted by different symbols) for the 3-objective DTLZ2 test problem using NSGA-II (left), IBEA (middle), and HypE with two preference points displayed as ellipsoids (right).

Figure 5 :

 5 Figure 5: Distribution of points, plotted in parallel coordinates, for the 7-objective DTLZ2 test problem for IBEA (top) and HypE with two preference points (solid black lines) and σε = 0.05 (bottom).

Table 1 :

 1 Normalized hypervolume values using a weight distribution wp that reflects two preference points. The mean value and twice the standard deviation of 50 runs are listed. HypEp uses the weight distribution function wp while HypEuni uses a uniform weight distribution.

		NSGA-II	IBEA	HypEuni	HypEp
	2d	0.443±0.321	0.658±0.156	0.487±0.202	0.910±0.065
	3d	0.480±0.353	0.956±0.015	0.917±0.035	0.996±0.005
	5d	0.021±0.071	0.807±0.033	0.959±0.121	0.995±0.004
	10d	0.215±0.311	0.688±0.121	0.964±0.049	0.999±0.002
	25d	0.171±0.280	0.613±0.322	0.826±0.152	0.999±0.003

Instead of the standard term reference point, see for example[START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF], we use the term preference point throughout this paper to reduce the likelihood of confusion with the hypervolume's reference point.

The hypervolume indicator satisfies IH(A) = E(αA(X)), where X is a random variable uniformly distributed on the objective space. Therefore from the strong law of large numbers the right-hand side of Eq. 2 converges to IH (A) for the number of samples N to infinity.

µ = (0.1377, 0.4131, 0.0688, 0.6196, 0.2065, 0.2754, 0.5507)


Anne Auger
email: anne.auger@inria.fr

Johannes Bader

Dimo Brockhoff

Eckart Zitzler

Investigating and Exploiting the Bias of the Weighted Hypervolume to Articulate User Preferences

Keywords: I.2.8 [Artificial Intelligence]: Problem Solving, Control Methods, and Search Algorithms

Optimizing the hypervolume indicator within evolutionary multiobjective optimizers has become popular in the last years. Recently, the indicator has been generalized to the weighted case to incorporate various user preferences into hypervolume-based search algorithms. There are two main open questions in this context: (i) how does the specified weight influence the distribution of a fixed number of points that maximize the weighted hypervolume indicator? (ii) how can the user articulate her preferences easily without specifying a certain weight distribution function?

In this paper, we tackle both questions. First, we theoretically investigate optimal distributions of µ points that maximize the weighted hypervolume indicator. Second, based on the obtained theoretical results, we propose a new approach to articulate user preferences within biobjective hypervolumebased optimization in terms of specifying a desired density of points on a predefined (imaginary) Pareto front. Within this approach, a new exact algorithm based on dynamic programming is proposed which selects the set of µ points that maximizes the (weighted) hypervolume indicator. Experiments on various test functions show the usefulness of this new preference articulation approach and the agreement between theory and practice.

INTRODUCTION

In the context of evolutionary multiobjective optimization, the hypervolume indicator is a set quality measure being a refinement of the Pareto dominance, i.e., whenever a Pareto set approximation entirely dominates another one, then the measure of the former will be larger. Hence, maxc ACM, 2009. This is the authors' version of the work. It is posted here by permission of ACM for your personal use. Not for redistribution. The definitive version was published at GECCO'09, July 8-12, 2009, Montréal Québec, Canada. http://doi.acm.org/10.1145/1569901.1569980 imizing the hypervolume indicator is, in general, in accordance to the Pareto dominance relation [START_REF] Knowles | On Metrics for Comparing Non-Dominated Sets[END_REF][START_REF] Fleischer | The measure of Pareto optima. Applications to multi-objective metaheuristics[END_REF][START_REF] Zitzler | On Set-Based Multiobjective Optimization[END_REF]. In recent years, this significant property has made the indicator popular as a selection criterion in evolutionary multiobjective optimizers [START_REF] Beume | SMS-EMOA: Multiobjective Selection Based on Dominated Hypervolume[END_REF][START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF][START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF].

An important question in this context concerns the bias of the hypervolume indicator when introduced into search. In a recent study [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF], this bias has been investigated theoretically by characterizing the distribution of µ points that maximize the hypervolume indicator, so-called optimal µdistributions. When hypervolume-based algorithms are employed, the inherent bias of the hypervolume indicator is often unwanted but rather the incorporation of the user's own preference into search is desired. Several approaches for articulating these user preferences are known from the literature, e.g., by defining preference points [START_REF] Deb | Reference Point Based Multi-Objective Optimization Using Evolutionary Algorithms[END_REF], specifying preferred search directions [START_REF] Deb | Interactive Evolutionary Multi-Objective Optimization and Decision-Making using Reference Direction Method[END_REF] or defining linear minimum and maximum tradeoffs [START_REF] Branke | Guidance in Evolutionary Multi-Objective Optimization[END_REF]. For a general overview of articulating user preferences, we refer to [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF][START_REF] Coello Coello | Handling Preferences in Evolutionary Multiobjective Optimization: A Survey[END_REF][START_REF] Rachmawati | Preference Incorporation in Multi-objective Evolutionary Algorithms: A Survey[END_REF].

Recently, a general framework to incorporate user preferences into hypervolume-based search algorithms has been proposed [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF]. The so-called weighted hypervolume indicator allows stressing certain regions of the objective space by applying a weight distribution. It has been shown for three different weight distribution functions that optimizing the weighted hypervolume indicator results in solutions clustered in regions with higher weight whereas regions with low weight contain only a few solutions. Further work improved the applicability of this weighted hypervolume approach to problems with many objectives by means of Monte Carlo sampling [START_REF] Auger | Articulating User Preferences in Many-Objective Problems by Sampling the Weighted Hypervolume[END_REF].

Two main questions arise in this context: [START_REF] Auger | Articulating User Preferences in Many-Objective Problems by Sampling the Weighted Hypervolume[END_REF] what is the bias of the weighted hypervolume indicator, i.e., how can the distribution of µ points maximizing the weighted hypervolume be characterized ? [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] how can the user articulate and integrate her preferences easily into the search without defining the weight distribution functions by hand?

Both questions are tackled in this paper. First, we investigate the optimal distribution of µ points optimizing the weighted hypervolume indicator; generalizing the results for the unweighted case presented in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] to the weighted case. Second, we propose a new approach to articulate user preferences by defining a reference front and a desired density of points on that front. This user preference is automatically transferred into a weight distribution function for the weighted hypervolume of [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF], the optimization of which guarantees that the resulting distribution of solutions on the actual front follows the desired density. This paper is organized as follows: In Sec. 3 we generalize the theoretical results on optimal µ-distributions from [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] to the weighted hypervolume. Based on the theoretical results on optimal µ-distributions, we propose in Sec. 4 a new approach for articulating user preferences in terms of defining a desired density of points on a predefined front for biobjective problems. In Sec. 4.2 we incorporate the weighted hypervolume into search and propose a new exact algorithm for selecting the best q individuals from a set of p solutions that maximize the overall (weighted) hypervolume by means of dynamic programming for biobjective problems. In Sec. [START_REF] Beume | Effects of 1-Greedy S-Metric-Selection on Innumerably Large Pareto Fronts[END_REF], we apply the new approach to several test problems for scenarios with known and unknown Pareto fronts.

PRELIMINARIES

In this study we consider biobjective optimization problems F : X → Z. The two objective functions (F1(x), F2(x)) = F(x) map the decision space X to the objective space Z ⊆ R 2 and are, without loss of generality, to be minimized. Let denote the weak Pareto dominance relation defined as x y if and only if F1(x) ≤ F1(y) and F2(x) ≤ F2(y). The Pareto-optimal set Ps then consists of all solutions x * ∈ X, such that ∄ x ∈ X with x x * and x *

x. The image of Ps under F is called Pareto front or front for short. For the sake of simplicity, we describe the Pareto front in terms of a function f mapping the image of the Pareto set under F1which is presumed to be a closed interval [xmin, xmax]-onto the image of the Pareto set under F2, see Fig. 1(a).

Throughout this study, we assume the optimization goal of maximizing the weighted hypervolume indicator of [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF], i.e., finding a set A of µ solutions that maximizes the weighted hypervolume. A set of µ points maximizing the unweighted hypervolume indicator has been denoted in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] as an optimal µ-distribution and we use the same term here for the weighted case as well. For simplicity, we define the weighted hypervolume indicator only for sets of objective vectors instead for solution sets A ⊆ X as it was already done in [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] and [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF]. For a set of µ biobjective objective vectors S = {(x µ 1 , y µ 1 ), . . . , (x µ µ , y µ µ )}, we define the weighted hypervolume indicator I µ H,w according to [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] and as a generalization of [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] as

I µ H,w (S) := µ X i=1 Z x µ i+1 x µ i Z y µ 0 y µ i w(x, y)dy ! dx (1)
where x µ µ+1 = r1, y µ 0 = r2, and r = (r1, r2) is the reference point of the hypervolume indicator, see Fig. 1(a). Note that for points (xi, f (xi)) that all lie on the front described by f , Eq. 1 results in a µ-dimensional function I µ H,w (x µ 1 , . . . , x µ µ ) that is equivalent to Eq. 2 in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] if we choose w(x, y) = 1.

THEORETICAL BACKGROUND

The theoretical results presented in this section build the basis of the new preference articulation approach proposed later on in this paper. To this end, we investigate the influence of the weighted hypervolume indicator on optimal µ-distributions by deriving the limit density when µ goes to infinity. We then build on this result and provide for a desired density the corresponding weight that will allow a user to articulate his preferences.

Limit Density of Optimal µ-Distributions for the Weighted Hypervolume

In [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF], we argued that, besides for linear fronts, it is very difficult to determine the optimal µ-distributions for arbitrary fronts exactly and therefore investigated optimal µdistributions in the limit case of µ to ∞. This limit is expressed as a density-approximating the percentage of points in a particular segment of the front. Here, we generalize this approach so as to obtain the density of points associated to the weighted hypervolume. Before we do so, note that optimal µ-distributions exist in the weighted case under the same assumption as for the unweighted case in [2, Theorem 1]: if the function f describing the Pareto front is continuous, there exists (at least) one set of µ points maximizing the weighted hypervolume 1 .

We assume without loss of generality that xmin = 0 and that f : 5]. We also assume that f is continuous within [0, xmax], is differentiable and that its derivative is a continuous function f ′ defined in the interval ]0, xmax[. An optimal µ distribution is defined as a set of µ points maximizing the weighted hypervolume indicator I µ H,w ((x µ 1 , . . . , x µ µ )). However, instead of maximizing the weighted hypervolume indicator I µ H,w , it is easy to see that, since r1r2 is constant, one can equivalently minimize

x ∈ [0, xmax] → f (x) with f (xmax) = 0 [2, Fig.
r1r2 -I µ H,w ((x µ 1 , . . . , x µ µ )) = µ X i=0 Z x µ i+1 x µ i Z f (x µ i ) 0 w(x, y) dy dx with x µ 0 = 0, f (x µ 0 ) = r2, and x µ µ+1 = r1 (see Fig. 1(b)
). If we subtract the area below the front curve, i.e., the integral

R xmax 0 ( R f (x) 0 w(x, y)dy)dx of constant value (Fig. 1(c)), we see that minimizing µ X i=0 x µ i+1 Z x µ i f (x µ i ) Z 0 w(x, y) dy dx - xmax Z 0 f (x) Z 0 w(x, y) dy dx (2)
is equivalent to maximizing the weighted hypervolume indicator (Fig. 1(d)).

For a fixed integer µ, we consider a sequence of µ ordered points in [0, xmax], x µ 1 , . . . , x µ µ that lie on the Pareto front. We assume that the sequence converges-when µ goes to ∞-to a density δ(x) that is regular enough. Formally, the density in x ∈ [0, xmax] is defined as the limit of the number of points contained in a small interval [x, x + h[ normalized by the total number of points µ when both µ goes to ∞ and h to 0, i.e., δ(x) = limµ→∞ h→0 "

1 µh P µ i=1 1 [x,x+h[ (x µ i )
" . As explained above, maximizing the weighted hypervolume is equivalent to minimizing Eq. 2, which is also equivalent to minimizing

Eµ = µ " µ X i=0 Z x µ i+1 x µ i Z f (x µ i ) 0 w(x, y)dy ! dx- Z xmax 0 Z f (x) 0 w(x, y)dy ! dx # . (3)
1 The proof is similar to the one in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF]: the µdimensional function in Eq. 1 is continuous and upper bounded by the hypervolume of the entire front, i.e. R r 1

x min R r 2 f (xmax) 1 y>f (x) w(x, y)dydx. (a) (b) (c) (d)
Figure 1: Illustration of the idea behind deriving the optimal density: Instead of maximizing the weighted hypervolume indicator I µ H,w ((x µ 1 , . . . , x µ µ )) (a), one can minimize (b) which is equivalent to minimize the integral between the attainment surface of the solution set and the front itself (c) which can be expressed with the help of the integral of f (d).

We now conjecture that the equivalence between minimizing Eµ and maximizing the hypervolume also holds for µ going to infinity. Therefore, our proof consists of two steps: (1) compute the limit of Eµ when µ goes to ∞. This limit is going to be a function of a density δ. (2) Find the density δ that minimizes E(δ) := limµ→∞ Eµ. The first step therefore consists in computing the limit of Eµ.

Lemma 1. If f is continuous, differentiable with the deri- vative f ′ continuous, if x → w(x, f (x)) is continuous, if x µ 1 , . . . , x µ µ converge to a continuous density δ, with 1 δ ∈ L 2 (0, xmax) 2 , and ∃ c ∈ R + such that µ sup "" sup 0≤i≤µ-1 |x µ i+1 -x µ i | « , |xmax -x µ µ | « → c then Eµ converges for µ → ∞ to E(δ) := - 1 2 Z xmax 0 f ′ (x)w(x, f (x)) δ(x) dx . (4)
Proof. For the sake of readability, this proof has been sent to the appendix.

The limit density of the µ-distribution, as explained before, is minimizing E(δ). It remains therefore to find the density minimizing E(δ). This optimization problem is posed in a functional space, the Banach space L2 (0, xmax) and is also a constraint problem since the density δ has to satisfy the constraint J(δ

) := R xmax 0 δ(x)dx = 1.
The constraint optimization problem (P) that needs to be solved is summarized in:

minimize E(δ), δ ∈ L 2 (0, xmax) subject to J(δ) = 1 (P)
Theorem 1. The density solution of the constraint optimization problem (P) equals

δ(x) = p -f ′ (x)w(x, f (x)) R xmax 0 p -f ′ (x)w(x, f (x))dx . (5)
Proof. This proof is similar to the proof of Theorem 7 in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF] where -f ′ needs to be replaced everywhere by -f ′ w.

Deriving the Weight from Given Density

In the previous section, we characterized the density δ(x) of points that maximize the weighted hypervolume indicator for a given weight distribution w(x, y) and front shape f (x). Here, we use the result from the opposite direction: given user-defined preferences, expressed by a density3 , we are interested in the corresponding weight distribution. This allows to model user preferences in a concise manner by optimizing the weighted hypervolume indicator as we will see in the following.

Instead of specifying the density δ(x), which is the density projected onto the x-axis, it is more convenient for the user to define the density δF of points on the front itself, which results from δ(x) by dividing it by the norm of the tangent at the corresponding point on the front, i.e., δF (x) = δ(x)/ p 1 + f ′ (x) 2 . Transforming the density this way, we see that the density on the front δF (x) is propor-

tional to q -f ′ (x)w(x,f (x)) 1+f ′ (x) 2
which implies that

w `x, f (x) ´∝ 1 + f ′ (x) 2 -f ′ (x) • δF (x) 2 . ( 6
)
Note that the weight is a strictly positive finite function if -f ′ (x) is and that it peaks to infinity if the derivative of f either goes to zero or to -∞.

Since it is more convenient to define the desired density as a function δ φ F (φ) of the angle φ rather than as a function of x, in the following we use δF (x) = δ φ F (arctan (f (x)/x)). The weight distribution is only defined for points (x, f (x)) on the front. In order to make the weight useful for search, it is extended to the rest of the objective space in the following section.

ARTICULATING USER PREFERENCES

This section builds on the theoretical results from the previous section and proposes an approach to articulate user preferences in hypervolume-based evolutionary algorithms. First we present a way to generalize the weight distribution of Eq. 6 to the entire objective space. Second, we propose a method to perform environmental selection based on the weight distribution derived. Thereby, the weighted hypervolume-and at the same time the user's preferenceis optimized.

Extending the Weight Distribution to the Entire Objective Space

In principle, any non-negative function can be chosen as the extension of w to points not on the front. For practical applications, however, a few things need to be considered:

1. The number of solutions is finite which leads to contributions that extend beyond the front. If the weight varies greatly in this region, the resulting distribution of points might be a bad approximation of δ φ F . 2. During search, the solutions are not mapped to the Pareto front where the weight w(x, f (x)) is defined. In those cases, the resulting distribution should possibly meet the user's expectations and not hinder convergence to the front. 3. The corresponding hypervolume indicator should be calculable in reasonable time. To address all three issues, we propose to use a two-dimensional piecewise constant weight distribution function, where straight lines separate the objective space into stripes Si of constant weight denoted by wi, see Fig. 2 4 . A stripe Si usually covers more than one Pareto-optimal point such that the weight wi is not defined uniquely. Instead, we set the weight wi to the mean

wi = R γ i w(s) ds R γ i ds , γi : {x ∈ R|(x, f (x)) ∈ Si} → R 2 , t → (t, f (t)) (7)
where the curve γi denotes the part of the Pareto front Si. Since the weight usually varies along γi, the weight according to [START_REF] Coello Coello | Handling Preferences in Evolutionary Multiobjective Optimization: A Survey[END_REF] does not reflect (6) precisely. However, the approximation by wi can be made arbitrarily accurate by increasing the number of stripes.

Besides facilitating the calculation of the hypervolume indicator, the stripes have the advantage of the weight remaining constant if a solution is moving along the direction of the stripe. Many different arrangements of stripes are possible, see Fig. 2. Here, we propose using rays starting from the ideal point z I = (z I 1 , z I 2 ) at equally separated angles ranging from 0 • to 90 • to separate the stripes

S rays i =  (x, y)|φi ≤ arctan " x -z I 1 y -z I 2 « < φi+1 ff (8)
where φi = (i -1)/n • 90 • and 1 ≤ i ≤ n. For example, on the left-hand side of Fig. 2 nine rays and z I = (0, 0) are used. 4 In this setting, the corresponding hypervolume calculation can be done easily by determining intersections of straight lines and calculating the area of triangles and trapezoids.

Optimizing the Hypervolume in 2D

The weighted hypervolume indicator represented by the stripes proposed before, shall be used as the optimization criterion of an evolutionary algorithm in order to find solutions that meet the corresponding user preference as good as possible. Several evolutionary algorithms aiming at maximizing the hypervolume indicator in their environmental selection step which can be formulated as solving the hypervolume subset selection problem (HSSP) [START_REF] Bader | HypE: An Algorithm for Fast Hypervolume-Based Many-Objective Optimization[END_REF]: given a set of solutions A and 0 ≤ q ≤ |A|, find a subset A * ⊆ A with |A * | = q, such that the weighted hypervolume indicator of A * is maximal.

While for more than two objectives the HSSP problem is expected to be difficult and for this reason greedy heuristics are used to tackle the HSSP, e.g., in [START_REF] Beume | SMS-EMOA: Multiobjective Selection Based on Dominated Hypervolume[END_REF][START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF][START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF], we propose an efficient exact algorithm for the case of 2 objectives here-using the fact that the hypervolume contribution of an objective vector only depends on its two adjacent neighbors 5 . Exploiting this property, dynamic programming [START_REF] Cormen | Introduction to Algorithms[END_REF] can be used to solve the problem exactly in time O(|A| 3 ) as opposed to O(|A| 2 ) for the greedy approach by combining solutions of smaller subproblems P t-1 c in a bottom-up fashion to solutions for larger subproblems P t c : for a fixed solution ac ∈ A and a t ∈ {0, . . . , |A|}, we define the subproblem P t c as finding the set A t c ⊆ A of t solutions maximizing the hypervolume such that A t c contains ac and in addition, only elements a k lying to the right of ac, i.e., f1(ac) ≤ f1(a k ).

Obviously, {ac} is the solution for P 1 c . According to the above made statement, the solution for P t c with t > 1 can now be easily found when considering the unions of {ac} with the solutions of all P t-1 k with f1(ac) ≤ f1(a k ) and taking the resulting solution set with the highest hypervolume. Once the solutions for t = q are determined, the subset which then has the largest hypervolume corresponds to the solution to the overall problem. Algorithm 1 shows the pseudo code of the procedure where sets S t c of indices instead of the sets A t i are considered for clarity. We illustrate the algorithm by means of an example:

Example 1. Consider six objective vectors o1 to o6 of which one likes to choose those q = 3 that maximize the hypervolume, see Fig. 3. In the first stage (a), the optimal subsets of size 1 and their hypervolume value are calculated (Lines 1 and 2 in Alg. 1). Please note that some subsets do not exist or will not be used to build the overall solution and can therefore be neglected (dashes).

In the next stages, the subsets of size t = 2 to q (Lines 3-13) are determined for all individuals oc (Lines 4-12). To this end, the hypervolume of combining oc with any subset to its right of size t -1 (Lines 6-8) are calculated. For example, in Fig. 3(b) o3 is combined with the subset S 1 5 to form S 2 3 = {3, 5} with hypervolume h 1 3 = 127. In this way, all subsets of size 2 (c) and then of size 3 (d) are determined.

Reaching t = q, the optimal solution to the overall problem corresponds to the set with the largest hypervolume, in this example S 3 1 = {1, 3, 5} with value h 3 1 = 142 (Line 14). Note that the advantage of the exact algorithm over often used greedy approaches for HSSP is that it overcomes the non-convergence of greedy algorithms, see [START_REF] Beume | Effects of 1-Greedy S-Metric-Selection on Innumerably Large Pareto Fronts[END_REF] for details. 5 For example, in Figure 3(a) the hypervolume contribution of o2 is bounded by o1 and o3 but not by o4, o5 or o6. This in turn means, that the increase in hypervolume, when adding o2 to any subset whose left-most element is o3, is equal.

1 3 S 3 { } = 1 4 S 4 { } = 1 5 S 5 { } = 1 6 S 6 { } = 112 = 1 3 h = 1 4 h = 1 5 h 1 5 h = 1 5 h 1 6 h (a) 1 3 S 3 { } = 1 5 S 5 { } = 2 3 h = (b) 2 2 S 2 3 { , } = 2 3 S 3 5 { , } = = 2 2 h = 2 3 h = 2 4 h = 2 5 h + + (c) 3 1 S = 3 1 S = 3 1 S = = 3 3 h = 3 4 h = 3 2 h = 3 1 h (d)

EXPERIMENTS

In this section, the theoretical concept to model (Sec. 3) and to optimize (Sec. 4) user preferences by a weighted hypervolume indicator are examined by simulations. To this end, the algorithm and the test problems used are presented first. After that, the Pareto front approximations obtained by the algorithm are visually compared to the expected density that corresponds to the user's preference.

Experimental Setup

The approach derived in Sec. 4 is applied to different test problems and user preferences. Unless noted otherwise, the latter is modeled by 90 stripes starting from the origin which corresponds to a resolution of 1 • , see Eq. 8. The reference point of the hypervolume is set to (1000, 1000) to ensure the extremes in the optimal µ-distribution if possible [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF].

The evolutionary algorithm in use first generates an initial population of α individuals that code potential solutions as real vectors. Thereafter, µ parent individuals are selected uniformly at random to generate λ offspring using the SBX crossover operator and the variable-wise polynomial mutation operator [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF]. The recombination and mutation probabilities were set according to [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF]. After that, non-dominated sorting and front-wise removal is carried out as in [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF]. From the set F of equally ranked solutions, the omission of which would result in a population P ′ that is smaller than α, the q = α -|P ′ | solutions computed by Alg. 1 that solve the HSSP on F are taken to fill up the population to |P ′ | = α. The above steps are repeated for N generations. In this study N = 1, 000 and α= λ= µ = 50. The main difference between the described algorithm and standard hypervolumebased algorithms such as SIBEA [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF] is the environmental selection step where Alg. 1 is used instead of a greedy strategy to solve the HSSP.

As test problems, ZDT1 (front shape f

(x) = 1 - √ x), ZDT2 (f (x) = 1 -x 2 )
, and ZDT3 (discontinuous front) [START_REF] Zitzler | Comparison of Multiobjective Evolutionary Algorithms: Empirical Results[END_REF] as well as DTLZ2 (f (x) = √ 1 -x 2 ) [START_REF] Deb | Scalable Test Problems for Evolutionary Multi-Objective Optimization[END_REF] are used. For some simulations, the objective functions are offset or scaled.

Algorithm 1 HSSP-Solver(O, q, r,C w H ) Require: A matrix O := (oi,j)p×2, where the rows represent the objective vectors sorted in ascending order according to the first objective; subset size q ≥ 1; reference point r = (r1, r2). C w H ( l, u) returns the weighted hypervolume of the rectangle from l to u. Ensure: Returns a set S that references rows of O that maximize the weighted hypervolume. 1:

h 1 i ← hi := C w H `(oi,1, oi,2), r ´∀ 1 ≤ c ≤ p 2: S 1 c ← {c} ∀ q ≤ c ≤ p
⊲ optimal subsets 3: for t = 2 to q do ⊲ bottom-up approach 4:

for c = q -t + 1 to p -t + 1 do ⊲ subproblem 5:

l ← (0, . . . , 0) 6:

for d = c + 1 to p -t + 1 do 7: l d = h t-1 d + C w H `(oc,
δ φ F (φ) ∝ 1 -(φ/90 • ) 3 4(b) δ φ F (φ) ∝  1, φ∈{[30 • ,50 • ],[70 • ,90 • ]} 0.1, otherwise 4(c) δ φ F (φ) ∝ 3 • `5 • 10 -3 + (φ -45 • ) 2 ´-1 + `5 • 10 -3 + (φ -10 • ) 2 ´-1 4(d) as in 4(b) 4(e) δ φ F (φ) ∝ 1 4(f) δ φ F (φ) ∝ 8 > < > :

Results

In the following, the distribution of points obtained when using different weight functions and test problems are presented. To this end, a graphical representation is used that shows multiple elements (see Fig. 4(a)): (i) the solutions found ( ) together with the Pareto front ( ); (ii) the front, specified by the user, that is used to calculate the weight rays ( ) if differing from the Pareto front; (iii) the weight per ray proportional to the length of gray bars ( ) and at intervals of 10 • the direction of the rays ( ) (iv) the density desired by the user ( ) as a function of the angle in polar coordinates with 0 lying at the radius where the weight bars start (e.g., at r = 1.2 for Fig. 4(a)); and (v) the obtained density represented by the solutions in terms of a histogram ( ) in the same coordinate system than the desired density.

Test Problems with Known Pareto Fronts

If the Pareto front is known, e.g., for test problems, the weight rays directly follow from Eq. 6 and 7. Figures 4(a Difficulties arise mainly where the front is very steep or flat and therefore a high weight is necessary to compensate this, e.g., at the extremes of DTLZ2 in Fig. 4(e).

The Typical Case of Unknown Fronts

If the front is unknown like in practice, one has to assume a front shape or find an approximate front shape by running an algorithm with respect to the unweighted hypervolume. For simplicity, we assume the front to be linear here. Then, the obtained density does not correspond to the desired density if the real front differs from the assumption or if the front is not reached. In Fig. 5(a) to 5(c) the weight distribution is derived for an assumed linear front, but the real front was (a) convex, (b) concave, and (c) concave while being scaled. In the latter example, the density is not affected by the scaling due to the radial nature of the weight rays.

If the assumed front lies at a different position than the real front, the obtained density differs from the desired one, see Fig. 5(d) and (e), even if the shape is correct as in (d).

Number of Weight Rays

Setting the number of rays to 90 does not lead to any visually noticeable impreciseness. However, if the resolution is decreased by using only a small number of rays, e.g., 15 in Fig. 5(f), the density obtained clearly reflects the step like character of the weight distribution.

CONCLUSIONS

The hypervolume indicator has an inherent characteristic in terms of how a set of µ solutions, optimizing the hypervolume indicator, is distributed on the Pareto front, which has been expressed in terms of density in [START_REF] Auger | Theory of the Hypervolume Indicator: Optimal µ-Distributions and the Choice of the Reference Point[END_REF]. However, the bias can be modified by applying weight distribution functions as shown in [START_REF] Zitzler | The Hypervolume Indicator Revisited: On the Design of Pareto-compliant Indicators Via Weighted Integration[END_REF]. Here, the two precursor studies have been combined to derive the density for the weighted case. On the one hand, this allows to determine the bias of the weighted hypervolume in a concise manner. On the other hand, the dual viewpoint has been taken to determine the weight distribution for a given density and front shape. This enables to express user preferences in a very precise way by a unary indicator.

To utilize the indicator for search, an approximation of the weight distribution function by means of piecewise constant rays has been suggested. It allows simplifying the hypervolume calculation while maintaining desirable properties. Additionally, an algorithm based on dynamic programming that selects solutions yielding maximal hypervolume value has been proposed. A visual comparison of the obtained Pareto front approximations with the expected densityprovided by the user-has shown the practical usefulness of the theoretical concept as well as the efficiency of the proposed algorithm. The approach is restricted to biobjective problems, though. Therefore, we should aim at extending it to more objectives in future studies. The approach has also to show its usefulness in real-world applications.

(a) (b) (c) (d) (e) (f)

APPENDIX Proof of Lemma 1

Let us first note that the Cauchy-Schwarz inequality implies that

Z xmax 0 |f ′ (x)w `x, f (x) ´| |δ(x)| dx ≤ s Z xmax 0 `f ′ (x)w(x, f (x)) ´2dx Z xmax 0 (1/δ(x)) 2 dx (9)
and since x → f ′ (x)w(x, f (x)) ∈ L 2 (0, xmax) and 1 δ ∈ L 2 (0, xmax), the right-hand side of Eq. 9 is finite and Eq. 4 is well-defined.

Step 1. In a first step we are going to prove that Eµ defined in Eq. 3 satisfies

Eµ = µ µ P i=0 " -1 2 f ′ (x µ i )w(x µ i , f (x µ i ))(x µ i+1 -x µ i ) 2 + O((x µ i+1 -x µ i ) 3 ) " (10)
To this end, we elongate the front to the right such that f equals f w(x, y) dy dx [START_REF] Deb | Reference Point Based Multi-Objective Optimization Using Evolutionary Algorithms[END_REF] At the first order, we have that

Z f (x µ i ) f (x)
w(x, y)dy = w(x

µ i , f (x µ i ))(f (x µ i )-f (x))+O((x -x µ i )) ( 12
)
Since f is differentiable, we can use a Taylor approximation of f in each interval [x µ i , x µ i+1 ] and write f (x) = f (x µ i ) + f ′ (x µ i )(x -x µ i ) + O((x -x µ i ) 2 ), which thus implies that f (x µ i ) -f (x) = -f ′ (x µ i )(x -x µ i ) + O((x -x µ i ) 2 ) and thus the left hand side of Eq. 12 becomes -w(x µ i , f (x µ i ))f ′ (x µ i )(x -x µ i ) +O((x -x µ i ) 2 ) By integrating the previous equation between x µ i and x µ i+1 we obtain

Z x µ i+1 x µ i Z f (x µ i ) f (x)
w(x, y) dy dx = -0.5w(x µ i , f (x µ i ))f ′ (x µ i )(x µ i+1 -x µ i ) 2 + O((x µ i+1 -x µ i ) 3 )

Figure 2 :

 2 Figure 2: Different arrangements of stripes to represent a weight distribution function: radial (left), parallel (middle) or freely (right).

Figure 3 :

 3 Figure 3: Three out of six objective vectors need to be selected. Algorithm 1 starts by calculating subsets of size 1 (a). Then the results (b) are combined to sets of size of 2 (c) and finally of size 3 (d).

5

  (a) to 5(f) as in 4(c)

  ) to 4(c) show the results on ZDT1 for (a) stressing the second objective (b) desiring a piecewise constant density and, (c) preferring two preference points. Analytic descriptions of (a) (b)

Figure 4 :

 4 Figure 4: Shows the solutions found together with the Pareto front for different test problems, i.e., ZDT1 (a)-(c), ZDT2 (d), DTLZ2 (e) and ZDT3 (f ). The underlying densities (listed in Tab. 1) are shown as a function of the angle, as well as the corresponding weight and the obtained density (see text for details).

Figure 5 :

 5 Figure 5: See Fig. 4(b) for a legend to the visual elements. Additionally, the front used to calculate the weight is depicted. The test problems are (modified versions) of ZDT1 for (a), (d) and (f ); ZDT2 (b); and DTLZ2 (e). The underlying density is the same in all six cases (see Tab. 1), consult the text for details.

  (xmax) = 0 for x ∈ [xmax, x µ µ+1 ]. Like that, we can decompose , y)dydx = 0. Using the right hand side of the previous equation in Eq. 3, we find that Eµ = µ

Table 1 :

 1 1, oc,2), (o d,1 , r2) Densities used for the experiments.

	8:

L 2 (0, xmax) is a functional space (Banach space) defined as the set of all functions whose square is integrable in the sense of the Lebesgue measure.

Practically speaking, the user can define the desired density for example with the help of a graphical user interface in terms of piecewise linear or polynomial functions. For the user, this might be more convenient than specifying a weight distribution function over the entire objective space.

Summing up for i = 0 to i = µ, multiplying by µ and using Eq. 11 we obtain Eq. 10, which concludes Step 1.

Step 2. We now decompose

For the sake of convenience in the notations, for the remainder of the proof, we redefine x µ µ+1 as xmax such that the previous equation becomes 1 2

For µ to ∞, the assumption µ sup((sup 0≤i≤µ-1 |x µ i+1 -x µ i |), |xmaxx µ µ |) → c implies that the distance between two consecutive points |x µ i+1 -x µ i | as well as |x µ µ -xmax| converges to zero. Let x ∈ [0, xmax] and let us define for a given µ, ϕ(µ) as the index of the points such that x µ ϕ(µ) and

) and therefore by integrating between x µ ϕ(µ) and x µ ϕ(µ)+1 we obtain

Moreover by definition of the density δ,

Using Eq. 14 and Eq. 15, we thus have

Therefore for every i we have that

Since x → f ′ (x)w(x, f (x))/δ(x) is continuous, we also obtain

2 ) Injecting Eq. 16 in the previous equation, we obtain

Multiplying by 1/2 and summing up for i from 0 to µ and using Eq. 13 we obtain 1 2

Let us define ∆µ as sup((sup 0≤i≤µ-

By assumption, we know that µ∆µ converges to a positive constant c. The last term of Eq. 17 satisfies

where K > 0. Since µ∆µ converges to c, (µ∆µ) 2 converges to c 2 . With ∆µ converges to 0, we therefore have that µ 2 ∆ 3 µ converges to 0. Taking the limit in Eq. 17 we therefore obtain