MULTI-OBJECTIVE optimization (MOO), also calledmultiple criteria optimization [1], is an extension of the standard single-objective optimization, where the objectives may be conflicting with each other [2], [3]. When a conflict exists, we are no more looking for a single optimal solution but for a set of solutions, each one providing a trade-off on the objectives and none being better than the others. This solution set is called the Pareto set and its counterpart in the objective space is denoted the Pareto frontier. The Pareto frontier is at the core of MOO algorithms, being the foundation of many methods devoted to evaluating the performance and comparing the solutions to each other [4]. However, the frontier is defined by the objectives, which can be expensive to compute [5], [6], [7]. This leads to a variety of surrogate methods that try to approximate the objectives, e.g. [8], [9], thus saving computational resources. Among the surrogates that directly or indirectly estimate the Pareto frontier, one introduced by Yun et al. [10] is the closest to the surrogate described in this paper. They used a one-class SVM to define a function over the objective space whose null

C. S. Miranda and F. J. Von Zuben are with the School of Electrical and Computer Engineering, University of Campinas (Unicamp), Brazil. E-mail: contact@conradomiranda.com, vonzuben@dca.fee.unicamp.br

space describes an approximation of the Pareto frontier. This function is used to select individuals, since its value increases as its argument becomes more distant from the frontier, which are then used for crossover in a genetic algorithm.

Loshchilov et al. [11] presented a similar SVM approach, but the function learnt is defined over the decision space, which allows direct comparison with the Pareto frontier approximation without requiring evaluation of the objectives. This direct comparison can also be achieved with estimates built over the objective space by integrating surrogates for the objectives. However, contrary to the one-class SVM that learns a model to fit all samples on one side of the approximate frontier, the proposed SVM is also able to consider points that dominate the frontier being approximated, allowing approximation of multiple Pareto frontiers, each defined by a class of points in non-dominated sorting [12].

In a different approach, Loshchilov et al. [13] approximated the Pareto dominance instead of the Pareto frontier by using a rank-based SVM. In this case, instead of providing only the data points, the algorithm is also informed about the preference for an arbitrary number of sample pairs and tries to find a function where higher evaluation represents higher preference. Using the Pareto dominance to establish the preference between points and learning directly from the decision space, candidate solutions can be compared in dominance using the learnt function. However, both [11] and [13] try to estimate the Pareto frontier using generic function approximation models, which do not take into account the particularities of the Pareto frontier.

It is possible to guarantee that the Pareto frontier’s estimate is valid by building conservative estimates. For instance, using a binary random field over the objective space to model the boundary between dominated and non-dominated regions, Da Fonseca and Fonseca [14] described a theory that can be used to assess the statistical performance of a stochastic optimization algorithm and compare different algorithms. The attainment function described in the paper defines the probability that a run of the stochastic algorithm will dominate the function’s arguments. Although the attainment function is hard to compute, it can be approximated by multiple runs of the underlying algorithm, which makes it a good candidate for analyzing the performance statistics of the optimization algorithm and for performing hypothesis testing between MOO algorithms.

If a single run is considered, then the approximate attainment function describes a valid estimate of the Pareto frontier and it is defined as the border of the region dominated by the points provided. Although valid, this estimate is very conservative and does not interpolate between the points provided, which means it cannot provide a good idea of the frontier’s shape and any evaluation of new points could be performed using only dominance comparison with the provided points.

In this paper, we develop a theory that defines necessary and sufficient conditions for a functional description of a Pareto frontier. Based on this theory, the search for approximations for the Pareto frontier using surrogate functions should be constrained to, or at least focused on, the ones that satisfy the results. If not, the resulting manifold obtained from the function may have any shape, possibly with many dominated points, which could result in reduced performance.

Moreover, the theory is developed on the objective space, allowing either accurate or approximate objective evaluations to be used, without restricting the format of the objectives’ surrogates. If parametric surrogate objectives are used, their association with the Pareto frontier surrogate can provide feedback on how to adjust their parameters so that the approximation is closer to the real objectives.

As an example of how to integrate the theoretical conditions in a surrogate design, we show how to introduce the theoretical conditions as soft constraints in Gaussian processes [15], which are nonparametric models, thus being able to adjust to variable number of samples, and whose hyper-parameters can be easily optimized.

To validate the hypothesis that surrogate methods that do not consider this theory may define invalid Pareto frontier approximations, the constrained Gaussian process is compared to a regular Gaussian process and to an existing SVM-based surrogate [10] and results show that the soft constrained Gaussian process finds good approximations maximally obeying the constraints according to the degree of flexibility of the model. On the other hand, the models that do not take into account the theory can violate greatly and arbitrarily the conditions for a valid Pareto frontier.

This paper is organized as follows. Section II introduces the notation and principles of multi-objective optimization used in this paper. Section III shows the conditions that a function must satisfy to define a Pareto frontier. These conditions are then used in Section IV to build a function to approximate a frontier given some points on it and the approximation is compared to an existing surrogate. Finally, Section V summarizes the findings and points out future directions for research.

A multi-objective optimization (MOO) problem is defined by a decision space X and a set of objective functions , where [16]. Since the framework is the same for maximization or minimization, we will consider that minimization is desired in all objectives. For a given point x in the decision space, the point defined by its evaluation using the objectives is its counterpart in the objective space . Although the objective space usually only makes sense when coupled with the decision space and objectives, which

allows for its infeasible region and Pareto frontier to be defined, we will work only with the objective space in this paper, which means that the results hold for any problem. We will also consider that , since any restriction for a specific problem is defined by means of the objectives and decision space constraints, and are handled transparently.

Furthermore, we assume that the optimal solutions describe a set of manifolds on , which correspond to curves in the 2D case and surfaces in the 3D case. Most multi-objective optimization problems have solutions with this property, with noticeable exceptions, such as: i) problems where some of the objectives do not conflict, so that only one of them should be used in the MOO problem with the other conflicting objectives, while the optimality of the ignored objectives is guaranteed because they were redundant; and ii) some problems with less decision variables D than objectives M, such as the Viennet function [17].

Since we are dealing with an optimization problem, we must define operators to compare solutions, like the operators < and are used in the mono-objective case. In MOO, this operator is the dominance.

Definition 1 (Dominance). Let y and be points in , the objective space. Then y dominates , denoted , if for all i.

The definition of dominance used in this paper is the same provided in [4], which allows a point to dominate itself. This relation is usually called weak dominance, but we call it “dominance” for simplicity, since it is the main dominance relation used in this paper. Another common definition is to require that for at least one i, and both definitions are consistent with the theory developed in this paper.

Definition 2 (Strong Dominance). Let y and be points in , the objective space. Then y strongly dominates , denoted , if for all i.

Once defined the comparison operator, we can divide the space Y in three sets: an estimated Pareto frontier, the set of points strongly dominated by the estimated frontier, and the set of points not strongly dominated by the estimated frontier.

Definition 3 (Estimated Pareto Frontier). A path-connected set of points is said to be an estimated Pareto frontier if no point in it strongly dominates another point also in F, that is, , and every point in the objective space except for F either strongly dominates or is strongly dominated by a point in F, that is, .

A set S is path-connected if there is a path joining any two points x and y in S and a path is defined by a continuous function with p(0) = x and p(1) = y. Therefore, if there is a continuous path of points in S that gets from any to , then S is path-connected. Based on this definition, an estimated Pareto frontier F divides the objective space in three disjoint sets: points strongly dominated by points in F, points that strongly dominates points in F, and F itself.

−2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0

Figure 1: Example of the definitions for a particular multi-objective problem. The estimated strict Pareto frontier is shown in a solid blue line, the estimated Pareto frontier F includes the solid and dashed blue lines, the dominated region D is shown on the top right red area, and the non-dominated regionD is shown on the bottom left green area.

Definition 4 (Estimated Strict Pareto Frontier). A set of points is said to be an estimated strict Pareto frontier if no point in it dominates another point also in , that is, , and every point in the objective space except for either dominates or is dominated by a point in , that is, .

Definition 5 (True Pareto Frontier). An estimated strict Pareto frontier is a true Pareto frontier if and only if, for all points in , there is no other feasible point in the objective space that dominates the point in the frontier, that is, . Moreover, for a given problem, the true Pareto frontier is unique.

The estimated Pareto frontier of Definition 3 is a generalization and an approximation of the true Pareto frontier in two ways: i) if the true Pareto frontier is discontinuous, then dominated points are added so that the estimated Pareto frontier F is path-connected while also guaranteeing that no point in it strongly dominates any other; and ii) the estimated Pareto frontier is simply a set of points that divide the space into dominated and non-dominated regions, without stating anything about the optimality of its points.

Consider, for instance, a problem where one of the objectives is given by

and the other is given by . Then the true Pareto frontier is given by

which clearly is not path-connected. However, if we add the set of points to , then the resulting path-connected set satisfies Definition 3, despite the fact that every point in is dominated by , but not strongly dominated by it.

Figure 1 shows an estimated strict Pareto frontier , which coincides with the true Pareto frontier in this example, and the path-connected estimated Pareto frontier F for this problem. This makes it clear that the estimated Pareto frontier F can contain the true Pareto frontier , i.e. , while providing a path-connected 1D manifold that splits the whole objective space . Of course, these properties of the estimated Pareto frontier are extensible to M > 2 objectives.

With the definition of an estimated Pareto frontier, the objective space is divided into two sets, named dominated and non-dominated sets, also shown in Fig. 1.

Definition 6 (Dominated Set). The dominated set D for an estimated Pareto frontier F is the set of all points in where, for each one of them, there is at least one point in F that strongly dominates it, that is, .

Definition 7 (Non-Dominated Set). The non-dominated setD for an estimated Pareto frontier F is the set of all points that are not in F or D. This implies that.

Note that, from the definition of strong dominance, both D andD are open and unbounded sets, with boundaries defined by the estimated Pareto frontier F. Furthermore, if F contains the true Pareto frontier, then the points inD are not achievable due to the objectives’ definitions.

From the partition of the objective space in three sets, one estimated Pareto frontier, one dominated and one non-dominated set, we can define a score function similarly to [11], [13].

Definition 8 (Score Function). A score function R for a given estimated Pareto frontier F is a function that satisfies

Therefore, a score function provides a single value that places its argument in relation to the estimated frontier. Moreover, for a given estimated Pareto frontier F, there are many possible choices of score functions f(y) that satisfy the definition and all of them uniquely define F based on their solution set f(y) = 0. This allows a score function to work as a surrogate for the estimated Pareto frontier.

In this section, we will show how a score function f(y) can induce an estimated Pareto frontier F and the conditions it must satisfy so that the set it defines is indeed an estimated Pareto frontier, that is, no point in it strongly dominates any other point in it.

The main theory developed is based on the most general notion of a function f, but the conditions may be hard to evaluate for a general case. Therefore, we will also provide corollaries that prove the results for functions with additional constraints, like continuous derivatives. Since some of these results depend on Taylor approximations and the first derivative at the required points may be zero, we must define a generalized gradient.

Definition 9 (Generalized Gradient). Let , where

is the class of functions where the first k derivatives exist and are continuous, with . Let be the first non-zero derivative of h evaluated at 0, that is,

where is not defined if h is a constant function or no i satisfies the inequality. Then

is the generalized gradient operator, which is undefined if there is no i that satisfies the inequality.

The role of the generalized gradient in the theory to be presented is to avoid issues with functions that may have null derivative at the points being evaluated but that are also increasing. Consider, for instance, the function , whose gradient is null at x = 0. This function is strictly increasing, but the first-order approximation using Taylor series is a constant. In order to consider small changes in the function’s argument, we must use first non-null derivative, which is the generalized gradient, as it will dominate the approximation.

The generalized gradient can be used in the Taylor approximation as , where , , and is the big-O notation. Since the result is based on being a small value, the exact power used to compute is not important for the approximation and the term is dominated by the other factors.

The extensions to continuous functions f rely on the generalized gradient of a single-parameter continuous function , derived from the original f, having different signs for opposite directions. However, it does not hold for functions where is even.

For example, consider , which has . The Taylor approximation is given by , which does not give different signs to different directions of x. Therefore, the two constraints on defined in the corollaries that follow can be viewed as a single constraint on plus the constraint that is odd.

A. Necessary Conditions

The necessary conditions derived are direct applications of the estimated Pareto frontier’s definition and establish the basic ground on how to define a function f from a given estimated frontier.

Lemma 1 (General Necessity). Let F be an estimated Pareto frontier. Let be a score function for F. Then and for all , and .

Proof. Assume there are y, u, and such that 0. Let , so that .

If , then from the definition of a score function there is some such that . From the transitivity of dominance, we have that , which is a contradiction, since the point in the frontier cannot strongly dominate the point y also in the frontier. Then we must have , which means and also creates a contradiction.

Assume that , and let . Then we can similarly prove that it also creates a contradiction.

Therefore, there are no such y, u, and with or

This result is intuitive, since moving in direction u from y we enter either D orD. If the function has the required derivatives, then the following result holds.

Corollary 1 (Differentiable Necessity). Let F be an estimated Pareto frontier. Let be a score function for F. Let and , with . Let and be defined for all and . Then and for all and .

Proof. Since f satisfies all conditions from Lemma 1, we have that and for all y, u, and .

In particular, let . Approximating using Taylor series, we have that and , where is the appropriate power of for the expansion. Since f(y) = 0 and , then and must hold.

Although this corollary may appear to provide weaker guarantees on f, its proof shows that the inequality constraints on the generalized gradient is equivalent to the direct inequalities on the function defined in the previous lemma.

B. Sufficient Conditions

Once defined how the estimated Pareto frontier relates to a given score function, we will show that a function that satisfies the results of the previous lemma and corollary in fact uniquely defines an estimated Pareto frontier F.

Lemma 2 (General Sufficiency). Let be a function. Let be a path-connected set. Let and for all , , and . Then F is an estimated Pareto frontier.

Proof. For F to be an estimated Pareto frontier, we have to prove that for any we have . Assume there are y and in F such that .

Let and . Then we have , which violates the first inequality on . Alternatively, we have , which violates the second inequality.

Therefore, there are no y and in F such that , and F is an estimated Pareto frontier.

The restrictions on may be hard to verify in general, since they must be valid for all . However, if the function has the appropriate derivatives, then it becomes easier to check if it satisfies the requirements.

Corollary 2 (Differentiable Sufficiency). Let be a function. Let be a path-connected set. Let and xu), with . Let and for all and . Then F is an estimated Pareto frontier.

Proof. To use Lemma 2, we must prove that and for all , and .

Suppose there is some y, u, and in the domain such that . Moreover, let be the smallest value for which this happens for a given y and u. Let and . Then and , where is the appropriate power of for the approximation. However, cannot go from positive to negative without passing through 0 due to its continuity. Then there must be some such that , which contradicts the definition of .

Therefore, the first inequality on Lemma 2 holds. We can use a similar method to prove the second inequality, and then use the lemma.

Again, this corollary shows the equivalence between the inequalities on the function and on the generalized gradient.

C. Necessary and Sufficient Conditions

Since the symmetry between Lemmas 1 and 2 is clear, we can build a theorem to merge those two and provide necessary and sufficient conditions for defining an estimated Pareto frontier F from a score function f(y).

Theorem 1 (General Score Function). Let be a function. Let be a path-connected set. Let and . Let , and f(y) < . Then F is an estimated Pareto frontier if and only if and for all , and .

Proof. Assume that the constraints on f are valid. Then, from Lemma 2, we have that F is an estimated Pareto frontier. Now assume that F is an estimated Pareto frontier. Then, from Lemma 1, we have that the constraints on f are valid.

Instead of requiring knowledge of the sign of f(y) over the sets, we can use a more strict definition, requiring continuity, to guarantee that the result holds.

Corollary 3 (Continuous Score Function). Let R be a continuous function where there are points and such that , and . Let be a path-connected set. Then F is an estimated Pareto frontier if and only if and for all , and .

Proof. Assume that F is an estimated Pareto frontier. Assume that there are such that f(y) > 0 and . From the continuity of f, we have that there is some such that f(z) = 0. However, since f(z) = 0, it is in F. From the definition of D, there is some such that , which violates the assumption that F is an estimated Pareto frontier. Therefore, all points in D have the same sign over f. The same can be shown forD.

Since , we have that and . Then f satisfies all conditions from Theorem 1.

Again, we can replace the constraints on by the constraint on the generalized gradient.

Corollary 4 (Differentiable Score Function). Let be a function where there are points and such that , and . Let be a path-connected set. Let and . Let and be defined for all and . Then F is an estimated Pareto frontier if and only if and for all and .

Proof. We can use Corollary 3 to show that the restrictions on must hold. From Corollaries 1 and 2, we know that the restrictions on are the same as the restrictions on , so this corollary is valid.

After showing what conditions the function f must satisfy, one could ask how to build such function for a given problem and specially how to learn one from a given set of non-dominated points. This can be a hard question to answer in general, but we can provide an additional lemma that can help in many cases.

Lemma 3 (Strictly Increasing Sufficiency). Let R be a strictly increasing function on each coordinate. Let . Then F is an estimated Pareto frontier.

Proof. For F to be an estimated Pareto frontier, we have to prove that for any we have . Assume there are y and in F such that .

Let be a path between y and that increments only one coordinate at a time. Since f is strictly increasing, we have that . Thus , which contradicts the premise that because they are both in the frontier.

Therefore, there are no y and in F where and F is an estimated Pareto frontier.

Note that, because f is strictly increasing, there is no point in F that even dominates another point in F, which was allowed in Definition 3. This restriction can be relaxed to be only monotonically non-decreasing if one can guarantee that f(y) = 0 is only a manifold, and not a subspace with volume. If f(y) = 0 is a subspace, then we can find two points in it where one dominates the other, which violates the basic definition of an estimated Pareto frontier. For instance, a function that is monotonically non-decreasing and is constant in at most one dimension at a time does not create a subspace on f(y) = 0.

Nonetheless, this lemma can be used as a guide on how to build a function for the general case. We will build a model that tries to approximate an estimated Pareto frontier from a few of its samples using an approximated monotonically increasing function based on Gaussian processes.

A. Gaussian Process As a Function Approximation Problem

Since the model should have enough flexibility to fit the given samples, an appropriate choice for a surrogate function is a Gaussian process, which always has enough capacity to fit the data. Before describing how a Gaussian process is used to approximate the Pareto frontier, we provide the reader with an overview of how they work. For a more detailed description, we refer the reader to [15].

A Gaussian process (GP) is a generalization of the multivariate normal distribution to infinite dimensions and can be used to solve a regression problem. A GP defines a probability distribution over functions, such that the outputs are jointly normally distributed.

To better understand this concept, consider an infinite column vector and an infinite matrix . Then a function can be described by associating the row indexes, such that . The GP relies on the fact that the relationship between x and y can be written as:

which states that all dimensions of y are distributed according to a multivariate normal distribution with mean and covariance K(x). Moreover, the mean for a given dimension is given by and the covariance is given by , where is a positive semi-definite kernel function.

Although continuous functions, and thus Gaussian processes, are defined for an infinite number of points, which caused the vectors x and y to have infinite dimensions, only a finite number of observations are actually made in practice. Let N be such number of observations. Then, by the marginalization property of the multivariate normal distribution, we only have to consider N observed dimensions of x and y. Furthermore, the finite-dimension version of y is still normally distributed according to Eq. (1) when considering only the observed dimensions.

Usual choices for the mean and covariance functions are the null mean [15], such that , and the squared exponential kernel, defined by:

where and are the scale parameters, which define a representative scale for the smoothness of the function.

The choice of the kernel function establishes the shape and smoothness of the functions defined by the GP, with the squared exponential kernel defining infinitely differentiable functions. Other choices of kernel are possible and provide different compromises regarding the shape of the function being approximated, such as faster changes and periodicity of values. However, in order to use the monotonicity constraints introduced in Section IV-B, the kernel has to be at least twice differentiable, which limits the kernels that can be used.

Figure 2a shows the prior distribution over functions using the squared exponential kernel with , and the zero mean. This highlights the fact that the GP defines a distribution over functions, not a unique function. Three sample functions from this GP are also shown in the same figure. Note that the functions are not shown as continuous, which would require an infinite number of points, but as finite approximations.

To use the GP to make predictions, the observed values of x are split into a training set X, whose output Y is known, and a test set , whose output we want to predict. Since all observations are jointly normally distributed, we have that the posterior distribution is given by:

where are matrices built by computing the kernel function for each combination of the arguments values.

The posterior distribution for the previous GP, after four observations marked as black dots, is shown in Fig. 2b. Note that the uncertainty around the observed points is reduced due to the observation themselves, and the mean function passes over the points, as expected. Again, three functions are sampled from the posterior, and all agree on the value the function must assume over the observations.

In order to avoid some numerical issues and to consider noisy observations, we can assume that the covariance has a noisy term. Assuming that , where is normally distributed with zero mean and variance , then the covariance of the observations is given by . The noiseless value can then be estimated by:

, X, XX, X(4c) X, X, (4d)

which is similar to Eq. (3), except for the added term in corresponding to the noise.

B. Gaussian Processes with Monotonicity Soft Constraint as Surrogates

Just like in the previous section, we consider the null mean function and the squared exponential kernel defined

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −4

Figure 2: Function distribution using a Gaussian process. Before the observations, the distribution is the same over all the space. After the observations, the distribution adapts to constraint the possible functions. The distribution mean is given by the black line and the 95% confidence interval is given by the shadowed region. Three function samples are also provided for each case.

in Eq. (2). Since we are mapping from the objective space to a value in R, according to Definition 8, the input values are the objectives y and the outputs the scores z. Let be a set of N input points and

their desired targets for training. We define the latent variable L between the two, such that

where . The latent variable then produces the observed values Z through

where I is the identity matrix.

This model is the same as the one described in Section IV-A. However, only the mean prediction will be used in this paper to describe the estimated Pareto frontier. Moreover, we will show how changing the allowed noise level affects the Pareto frontier approximation.

Besides the observations of f(y) at the desired points, the GP framework also accepts observations of its derivative, since differentiation is a linear operator [18], [19], that is, the derivative of a GP is also a Gaussian process. However, since we do not know the desired value of the gradient, only that it should be positive, from Corollary 4 and Lemma 3, forcing an arbitrary value may lead to reduced performance.

Another option is to introduce a probability distribution over the gradient in order to favor positive values, introducing monotonicity information [20]. This new distribution can be viewed as adding constraints to the Gaussian process, making it feasible to include the monotonicity information to the existing framework.

Ideally, the probability distribution over the gradient is the step function, which provides a probability of zero if the gradient is negative and the same probability for all positive gradients. However, the step function defines a hard threshold and does not allow small errors, which can cause some problems for the optimization. Therefore, a smooth function that approximates the step is used to define a soft constraint over the gradient.

Let be the indication that the i-th sample is monotonic in the direction . Then the following probability distribution can be used to approximate the step function:

where we assume the probit function as the derivative probability. Since the probit is a cumulative distribution function, its value ranges from 0 to 1 and it is monotonically increasing, which makes it a good approximation for the step function. The parameter allows us to define how strict the distribution should be, with approximating the step function or a hard constraint. In this paper, following the suggestion of [20], we use .

Since the monotonicity probability is not normal, it has to be approximated by a normal distribution to be used in the GP framework. To understand this, first consider the problem without the monotonicity constraints, which is given by Eq. (4). The probability distribution of the observation is given by:

where L is the latent variable for the training data, whose

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

Figure 3: Contours for the f(y) learned using a Gaussian process with derivative constraint. The black dots are the frontier points provided.

−0.5 0 0.5 1 1.5 −0.5

−0.5 0 0.5 1 1.5 −0.5

Figure 4: Contours for the f(y) learned using standard Gaussian process. The black dots are the frontier points provided.

probability distribution, computed by the Bayes’ rule, is

According to the model, the prior p(L|X) and the likelihoods p(Y |L) and are normal distributions, which makes all integrals tractable and all other distributions defined in the closed form presented in Eq. 4.

Now, considering the monotonicity constraints, let M be the monotonicity constraints and be the random variable associated with the derivative of the latent variables L. Then the probability distribution in Eq. (5) can be written as . Rewriting the posterior distribution over the latent variables, we get:

Because the distribution is not normal and every other distribution in Eq. (7) is normal, the integrals defined in Eqs. (6) and (7b) are intractable. Therefore, the distribution must be approximated by a normal distribution, which can be achieved using the expectation propagation algorithm [21], with the update equations described in [20]. The expectation propagation algorithm iteratively adjusts an unnormalized normal distribution to locally approximate the distribution defined by the soft constraints, such that , where is a normalization constant, is a mean vector with one value for each monotonicity constraint, and is a diagonal covariance matrix.

Besides this monotonicity constraint, we also would like that the errors between the provided values for the points z and their latent values l are small, so that the estimated shape of the Pareto frontier is closer to the true one. This can be achieved by placing a prior inverse-gamma distribution over , whose density is given by:

where is the gamma function. As , this prior is ignored, while indicates that there is no noise. In the results shown, we fix and vary .

We define f(y) as the final expected value , and the parameters are optimized to maximize the full likelihood, including gradient probability and prior, of the training data Y and Z. We also add the monotonicity constraint on all training data for all directions, but it should be noted that we can also add only monotonicity constraint at a point without defining its desired value. This allows us to find points that have f(y) = 0 but negative gradient and add the constraint on them, which in turn could improve the estimation.

To test the GP’s performance as a surrogate, we consider the two test frontiers whose samples are given by , which is a convex frontier, and , which is a concave frontier, both with . Note that the points were purposely selected to test the ability to model very sharp frontiers. However, using only the points defined by and leads to a solution where f(y) is almost 0 everywhere. To avoid this problem, we add a point (1, 1), with target value 1, to and a point (0, 0), with target value , to . The parameters for the Gaussian process are found using gradient ascent in the samples likelihood.

Figure 3 shows the resulting curves for different values of . The first thing we notice is that, although does not place any restriction on , which allows the observed points in the frontier to be far from their latent values that actually define the frontier, the resulting curve is still able to fit the general shape defined by the points provided.

As we reduce the value of , the observed variance is required to be smaller and the frontier shape gets better and better. Ideally, with , the latent points would be the same as the observed points, but this causes numeric problems due to the monotonicity information and can make it harder to satisfy the monotonicity constraint, due to the smoothness of the GP.

When we reduce the value to and beyond, the resulting frontier is not valid anymore, with noticeable points with negative derivative. However, the largest difference in the concave problem is between points (0.82, 1.055) and (0.2, 0.985), with a total reduction in of just 0.07, and a similar result is obtained for the convex case. Therefore, this approximation is still close to the correct frontier and could be used to evaluate proposed solutions because it was built with the theoretical developments of this paper in mind and tries to approximate them, which most likely provides better frontier estimates than methods that use traditional regression solutions, such as [10], [11], [13], where the manifold f(y) = 0 can have any shape.

To evaluate the effect of using the gradient constraint, Fig. 4 shows a similar GP but without any information on the gradients. Although the expected Pareto frontier is correctly identified, there are also many points that do not belong to the frontier and where f(y) = 0. Since the unconstrained GP had better frontier estimates for the extreme points than the constrained GP, as all points between them and the knee satisfy the conditions, it appears that not every point benefits from the gradient constraint.

Even though both GP models failed to fully satisfy the theoretical conditions, we consider that the GP with derivative restriction performed better, both because there are some parameter sets that are able to satisfy the frontier conditions and because it does not violate the restrictions as much. Moreover, if the variance, which is not shown but is higher for points far from the inputs provided, is taken into account, then the violations of the GP with derivatives occur in a region with higher uncertainty than the violations of the pure GP.

Therefore, despite the minor violations of the GP with derivative constraints, this approximation is still close to the correct frontier and could be used to evaluate the proposed solutions.

C. Comparison to Existing SVM Surrogate

The surrogate method introduced in [10], like the method proposed in this paper, is based on approximating the frontier directly from values in the objective space. This makes it a good candidate for comparison and validating the conjecture that existing methods may arbitrarily violate the conditions described in this paper.

The one-class SVM used in [10] is defined by the following optimization problem:

where and the feature-extraction function is defined implicitly by the kernel

which is similar to the kernel used for the GP.

One important difference between training an SVM and a GP is that the GP has a natural way to optimize its hyper-parameters by maximizing the data likelihood, which automatically defines a trade-off between fitting the data and model complexity. For the SVM, we must use cross-validation [22], which reduces the number of points available to fit the model, since the data must be divided in the training and validation sets.

To compare the surrogate methods, we use one test problem from [3], which is also used in [10] to show the behavior of the proposed SVM surrogate. The problem is given by:

We chose this problem because its true Pareto frontier is discontinuous, which creates sharp changes in its associated estimated Pareto frontier, just like in Fig. 1, and makes it harder to approximate.

We chose so that the samples provided should be almost perfectly classified and we constraint the scales in Eq. 2 to be equal, so that both methods can use the same features from the samples. The data set provided is composed of a grid with step 0.05 for both variables, which includes some points in the Pareto frontier. The full grid is used to fit the SVM because it provided better results than using just the non-dominated points, while only the non-dominated points are required for the GP.

Figure 5 shows the resulting approximations of the Pareto frontier using a GP with parameters learnt through gradient ascent in the data likelihood, like in Section IV-B, and an SVM with different values of . The GP learns an appropriate shape from the samples provided despite the discontinuity in the frontier, but also slightly violates the constraints during the gap in . Moreover, in the absence of any information about the shape in the interval , because no point was provided there, the GP extrapolates a valid shape for the Pareto frontier.

The SVM is highly dependent on the parameter . When it is small, the shape learnt is very conservative and does not follow the shape defined by the points in the frontier. On the other hand, when it is large, the surrogate fits the points in the frontier better but also may define a function that violates greatly the conditions to be a valid Pareto frontier. The best value for that does not violate the constraints in the interval is . However, for this value the GP provides a better approximation of the Pareto frontier, as shown in Fig. 5b. Increasing provides a better approximation, achieving a quality comparable to the GP, but also creates regions that violate the conditions to be a valid Pareto frontier more than the GP. Furthermore, defines a region that the SVM believes is part of the Pareto frontier but actually is very distant from it and inside the dominated region, as shown in Fig. 5c.

Besides these issues, the SVM also does not extrapolate well to the region . Close inspection shows that

Figure 5: Estimated frontiers using SVM with different values of and Gaussian process. The points in the data set that belong to the true Pareto frontier are shown as dots.

the dominated region defined by the SVM is finite, that is, it is described by a region in the objective space that is surrounded by an infinite region that the SVM believes is not dominated. This behavior shows that the learnt model carries no concept of the problem it is solving, which is to approximate a Pareto frontier, but describes a generic function approximation. The results in Fig. 5 provide evidence for the conjecture that existing methods proposed in the literature may arbitrarily violate the conditions described in this paper.

Furthermore, if only the points at the Pareto frontier were provided for learning, then the region defined by the SVM would enclose only these points and would ignore the dominated region. Thus the SVM method requires data in the dominated region while the GP method only requires the points at the frontier.

In this paper, we have introduced the necessary and suf- ficient conditions that functions must satisfy so that their solution space describes an estimated Pareto frontier. These conditions follow from the definition of an estimated Pareto frontier and are extended for differentiable functions, which allows easier verification of the conditions.

Based on these conditions, a Gaussian process (GP) was tested on toy problems with very sharp Pareto frontiers. The GP was extended to include the theoretical conditions as soft probabilistic constraints and a regularization term was added to avoid large deviations between the points and their latent values. The mean latent value is used as surrogate for the Pareto frontier, and some values of the regularization constant allowed a correct frontier estimate to be found.

However, when the regularization becomes too strong, the surrogate violates the constraints that define a valid estimated Pareto frontier on some points, but this occurs far from the given inputs and the deviation is small. This suggests that, even under these conditions, the proposed function could be used to provide insight on the shape of the true Pareto frontier, and possibly provide more realistic estimates than other methods that do not take the restrictions into consideration during their design.

To validate this hypothesis and the conjecture that existing surrogate methods may violate the conditions described in this paper, we compared the proposed GP with a one-class SVM used in [10] on one of the test problems described in the same paper. We showed that the GP again violates the constraints by small values and provide a good estimate for the Pareto frontier, while the SVM defined a worse estimate or violated the conditions more than the GP. Furthermore, the dominated region defined by the SVM is bounded by what it represents as the non-dominated region, while the GP correctly divides the space in two infinite areas.

Besides being a better surrogate for the Pareto frontier, the GP has the data likelihood as an innate measure that can be used to optimize its hyper-parameters and only requires data at the frontier. On the other hand, the SVM must use some method, like cross-validation [22], to optimize its hyper-parameters and it requires data in the dominated region to define better approximations.

We highlight that, although GP were used together with the theory on this paper to approximate the Pareto frontier, the theory is general and does not depend on the specific choice of the function descriptor. Therefore, other models that are able to deal with the constraints imposed by the theory, in either a soft or hard way, should be able to learn the desired shape of the Pareto frontier too. Nonetheless, we are not aware of any other method to create the score function in which the constraints are as easy to include as in the GP. Additionally, a GP provides robustness to changing the number of points used in the estimation.

Further investigations involve studying the behavior of the GP to approximate the Pareto frontier with real benchmarks and using some multi-objective optimization algorithm, such as NSGA-II [23], to provide the points. Since the objectives tend to be smoother than in the example frontier provided [24], we expect the estimated Pareto frontier described by a GP to fit the true Pareto frontier even better in these problems. If this is the case, we will investigate the possibility of integrating the frontier surrogate with other surrogate models for the objectives, so that all of them are learned directly and the number of function evaluations could be reduced.

Moreover, since the only requirement for the surrogate is that the Pareto frontier is approximated by the null space and the exact value on other parts of the objective space are not relevant, the GP could be used to fit a regression model on the individuals of a population where the target value is monotonically increasing in the objective space. Standard performance measures in multi-objective optimization, such as the class in non-dominated sorting [23] and the dominance count [25], satisfy this property and can be used as targets of the regression. In this case, the GP would not only define the Pareto frontier, but would also define a measure of the distance between a given point and the approximated Pareto frontier.

Another interesting line of research is to evaluate when the derivative constraints on the points provided is beneficial, since in some points it avoids incorrect association of other points with the frontier, like around the knee in the unconstrained GP shown in this paper, and in others it may make the estimated shape not satisfy the constraints, like the points in the constrained GP also shown in this paper. This could not only provide better fit, but may also increase the fitting speed, since less constraints need to be evaluated, which reduces the size of the GP and the number of expectation propagation steps required. Therefore an iterative algorithm that adds the constraints as needed should be pursued.

The authors would like to thank CNPq for the financial support.

[1] X. Gandibleux, Multiple Criteria Optimization: State of the Art Annotated Bibliographic Surveys, ser. International Series in Operations Research & Management Science. Springer US, 2006.

[2] K. Miettinen, Nonlinear Multiobjective Optimization. Springer US, 1999.

[3] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms. John Wiley & Sons, 2001.

[4] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G. Da Fonseca, “Performance Assessment of Multiobjective Optimizers: An Analysis and Review,” Evolutionary Computation, IEEE Transactions on, vol. 7, no. 2, pp. 117–132, 2003.

[5] Y. Jin, “A Comprehensive Survey of Fitness Approximation in Evolu- tionary Computation,” Soft computing, vol. 9, no. 1, pp. 3–12, 2005.

[6] J. Knowles and H. Nakayama, “Meta-Modeling in Multiobjective Opti- mization,” in Multiobjective Optimization. Springer, 2008, pp. 245–284.

[7] I. Voutchkov and A. Keane, “Multi-objective Optimization Using Surro- gates,” in Computational Intelligence in Optimization. Springer, 2010, pp. 155–175.

[8] B. Liu, Q. Zhang, and G. Gielen, “A Gaussian process surrogate model assisted evolutionary algorithm for medium scale expensive optimization problems,” Evolutionary Computation, IEEE Transactions on, vol. 18, no. 2, pp. 180–192, 2014.

[9] H. Karshenas, R. Santana, C. Bielza, and P. Larranaga, “Multiobjective estimation of distribution algorithm based on joint modeling of objectives and variables,” Evolutionary Computation, IEEE Transactions on, vol. 18, no. 4, pp. 519–542, 2014.

[10] Y. Yun, H. Nakayama, and M. Arakava, “Generation of Pareto frontiers using support vector machine,” in International Conference on Multiple Criteria Decision Making, 2004.

[11] I. Loshchilov, M. Schoenauer, and M. Sebag, “A Mono Surrogate for Multiobjective Optimization,” in Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation. ACM, 2010, pp. 471–478.

[12] X. Zhang, Y. Tian, R. Cheng, and Y. Jin, “An efficient approach to nondominated sorting for evolutionary multiobjective optimization,” Evolutionary Computation, IEEE Transactions on, vol. 19, no. 2, pp. 201–213, 2015.

[13] I. Loshchilov, M. Schoenauer, and M. Sebag, “Dominance-Based Pareto- Surrogate for Multi-Objective Optimization,” in Simulated Evolution and Learning. Springer, 2010, pp. 230–239.

[14] V. G. da Fonseca and C. M. Fonseca, “The attainment-function approach to stochastic multiobjective optimizer assessment and comparison,” in Experimental methods for the analysis of optimization algorithms. Springer, 2010, pp. 103–130.

[15] C. E. Rasmussen and C. K. I. Williams, Gaussian Process for Machine Learning. MIT Press, 2006.

[16] K. Deb, “Multi-objective optimization,” in Search methodologies. Springer, 2014, pp. 403–449.

[17] C. A. C. Coello, D. A. Van Veldhuizen, and G. B. Lamont, Evolutionary algorithms for solving multi-objective problems. Springer, 2002, vol. 242.

[18] A. O’Hagan, “Some Bayesian numerical analysis,” Bayesian statistics, vol. 4, pp. 345–363, 1992.

[19] C. E. Rasmussen, J. M. Bernardo, M. J. Bayarri, J. Berger, A. P. Dawid, D. Heckerman, A. F. M. Smith, and M. West, “Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals,” in Bayesian Statistics 7, 2003, pp. 651–659.

[20] J. Riihimäki and A. Vehtari, “Gaussian processes with monotonicity information,” in International Conference on Artificial Intelligence and Statistics, 2010, pp. 645–652.

[21] T. P. Minka, “Expectation propagation for approximate Bayesian infer- ence,” in Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence. Morgan Kaufmann Publishers Inc., 2001, pp. 362–369.

[22] C. M. Bishop, Pattern Recognition and Machine Learning. Springer, 2006.

[23] K. Deb, A. Pratap, S. Agarwal, and T. A. M. T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” Evolutionary Computation, IEEE Transactions on, vol. 6, no. 2, pp. 182–197, 2002.

[24] S. Huband, P. Hingston, L. Barone, and L. While, “A review of multiob- jective test problems and a scalable test problem toolkit,” Evolutionary Computation, IEEE Transactions on, vol. 10, no. 5, pp. 477–506, 2006.

[25] N. Beume, B. Naujoks, and M. Emmerich, “SMS-EMOA: Multiobjec- tive selection based on dominated hypervolume,” European Journal of Operational Research, vol. 181, no. 3, pp. 1653–1669, 2007.

Conrado S. Miranda received his M.S. degree in Mechanical Engineering and his B.S. in Control and Automation Engineering from the University of Campinas (Unicamp), Brazil, in 2014 and 2011, respectively. He is currently a Ph.D. student at the School of Electrical and Computer Engineering, Unicamp. His main research interests are machine learning, multi-objective optimization, neural networks, and statistical models.

Fernando J. Von Zuben received his Dr.E.E. degree from the University of Campinas (Unicamp), Campinas, SP, Brazil, in 1996. He is currently the header of the Laboratory of Bioinformatics and Bioinspired Computing (LBiC), and a Full Professor at the Department of Computer Engineering and Industrial Automation, School of Electrical and Computer Engineering, University of Campinas (Unicamp). The main topics of his research are computational intelligence, natural computing, multivariate data analysis, and machine learning. He coordinates open-ended research projects in these topics, tackling real-world problems in the areas of information technology, decision-making, pattern recognition, and discrete and continuous optimization.