Efficient computation of expected hypervolume improvement using box decomposition algorithms

In the field of multi-objective optimization algorithms, multi-objective Bayesian Global Optimization (MOBGO) is an important branch, in addition to evolutionary multi-objective optimization algorithms. MOBGO utilizes Gaussian Process models learned from previous objective function evaluations to decide the next evaluation site by maximizing or minimizing an infill criterion. A commonly used criterion in MOBGO is the Expected Hypervolume Improvement (EHVI), which shows a good performance on a wide range of problems, with respect to exploration and exploitation. However, so far, it has been a challenge to calculate exact EHVI values efficiently. This paper proposes an efficient algorithm for the exact calculation of the EHVI for in a generic case. This efficient algorithm is based on partitioning the integration volume into a set of axis-parallel slices. Theoretically, the upper bound time complexities can be improved from previously O(n2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O (n^2)$$\end{document} and O(n3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(n^3)$$\end{document}, for two- and three-objective problems respectively, to Θ(nlogn)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varTheta (n\log n)$$\end{document}, which is asymptotically optimal. This article generalizes the scheme in higher dimensional cases by utilizing a new hyperbox decomposition technique, which is proposed by Dächert et al. (Eur J Oper Res 260(3):841–855, 2017). It also utilizes a generalization of the multilayered integration scheme that scales linearly in the number of hyperboxes of the decomposition. The speed comparison shows that the proposed algorithm in this paper significantly reduces computation time. Finally, this decomposition technique is applied in the calculation of the Probability of Improvement (PoI).


Introduction
In multi-objective design optimization, the objective function evaluations are generally computationally costly, mainly due to the long convergence times of simulation models.A simple and common remedy to this problem is to use a statistical model learned from previous evaluations as the fitness function, instead of the 'true' objective function.This method is also known as Bayesian Global Optimization (BGO) [30].In BGO, a Gaussian Process (GP) model is used as a statistical model.In each iteration, the algorithm evaluates a new solution and updates the Gaussian Process model.A new solution is chosen by the score of an infill criterion, given a statistical model.For multi-objective problems, the family of these algorithms is called Multi-objective Bayesian global optimization (MOBGO).Compared to evolutionary multi-objective optimization algorithms (EMOAs), MOBGO requires only a small budget of function evaluations to achieve a similar result with respect to hypervolume indicator, and it has already been used in real-world applications to solve expensive evaluation problems [40].According to the authors' knowledge, BGO was used for the first time in the context of airfoil optimization in [27], and then applied in the field of biogas plant controllers [16], detection in water quality management [41], machine learning algorithm configuration [23], and structural design optimization [33].
In the context of Bayesian global optimization, an infill or pre-selection criterion is used to evaluate how promising a new point is.In single-objective optimization, the Expected Improvement (EI) is widely used as the infill criterion, and it was first introduced by Mockus et al. [30] in 1978.The EI exploits both the Kriging prediction and the variance in order to give a quantitative measure of a solution's improvement.Later, the EI became more popular due to the work of Jones et al [21].In MOBGO, a commonly used criterion is Expected Hypervolume Improvement (EHVI), which is a straightforward generalization of the EI and was proposed by Emmerich [10] in 2005.Compared with other criteria, EHVI leads to an excellent convergence to -and coverage of -the true Pareto front [5,37].Nevertheless, the calculation of EHVI itself so far has been time-consuming [23,32,36,41], even in the 2-D case 1 .Moreover, EHVI has to be computed many times by an optimizer to search for a promising solution in every iteration.For these reasons, a fast algorithm for computing EHVI is needed.
The first method suggested for EHVI calculation was Monte Carlo integration and was proposed by Emmerich [10,13].This method is simple and straightforward.However, the accuracy of EHVI highly depends on the number of iterations.The first exact EHVI calculation algorithm in the 2-D case was derived in [12], with time complexity of O(n3 log n).Here, n is the number of non-dominated points in the archive.The EHVI calculation algorithm in [12] partitions an objective space into boxes and then calculates the EHVI by summing all the EHVI values of each box.Couckuyt et al. [5] introduced an exact EHVI calculation algorithm (CDD13) for d > 2 by representing a non-dominated space with three types of boxes, where d represents the number of objective functions.The method in [5] was also practically much faster than those discussed in [12], though a detailed complexity analysis was missing.Hupkens et al. [19] reduced the time complexity to O(n2 ) and O(n 3 ) in the 2-D and 3-D cases, respectively.The algorithms in [19] improve the algorithms in [12] by two ways: 1) only summing the EHVI values of each box in a non-dominated space; 2) reusing some intermediate integrations during the EHVI calculation.The algorithms in [19] further improve the practical efficiency of EHVI on test data in comparison to [5].Recently, Emmerich et al. [14] proposed an asymptotically optimal algorithm with time complexity of Θ(n log n) in the bi-objective case.More recently, Yang et al. proposed an asymptotically optimal algorithm with time complexity Θ(n log n) in the 3-D case [38].The algorithm, KMAC2 in [38], partitions a non-dominated space by slices linearly and re-derives the EHVI calculation formulas.However, a generalization of this technique to more than three dimensions/objectives and the empirical testing of MOBGO algorithms on benchmark optimization problems, are still missing so far.
This paper mainly contributes to extending the state-of-the-art EHVI calculation methods into higher dimensional cases.The paper is structured as follows: Section 2 introduces the nomenclature, Kriging, and the framework of MOBGO; Section 3 provides some fundamental definitions used in this paper; Section 4 describes how to partition an integration space into (hyper)boxes efficiently, and how to calculate EHVI based on this partitioning method; Section 5 shows experimental results of speed comparison and MOBGO based algorithms' performance on 10 well-known scientific benchmarks in 6-and 18-dimensional search spaces; Section 6 draws the main conclusions and discusses some potential topics for further research.

Multi-objective Bayesian Global Optimization
A multi-objective optimization (MOO) problem is an optimization problem that involves multiple objective functions.A MOO problem can be formulated as: where d is the number of objective functions, y i (i = 1, • • • , d) are the objective functions, and a decision vector x is in an m-dimensional space.

Notations
The following table summarizes the notations used in this paper.The vectors in P, where P = (y (1) , . . ., Upper bounds of integration slices The multivariate independent normal distribution The integration slices in an i-dimensional case y Objective functions

Kriging
As a statistical interpolation method, Kriging is a Gaussian process based multivariate regression method.Compared with simulator-based evaluations in design optimization, one prediction/evaluation of the Kriging model is typically cheap [28].Therefore, Kriging is widely used as a popular surrogate model to approximate noise-free data in computer experiments.Kriging models are fitted from previously evaluated points.Given a set of n decision vectors X = {x (1) , x (2) , and associated function values Y(X) = y(x (1) ), y(x (2) ), • • • , y(x (n) ) , Kriging assumes Y(X) to be a realization of a random process Y of the following form [3,21]: where µ(x) is the estimated mean value over all given sampled points, and (x) is a realization of a Gaussian process with zero mean and variance σ 2 .The regression part µ(x) approximates the function Y globally and the Gaussian process (x) takes local variations into account.Opposed to other regression methods (such as support vector machine), Kriging/GP also provides an uncertainty qualification of a prediction.The correlation between the deviations at two decision vectors (x and x ) is defined as: Here R(., .) is the correlation function, which decreases with the distance between two points.It is common practice to use a Gaussian correlation function (also known as a squared exponential kernel): where θ i are parameters of the correlation model.They can be interpreted as a measurement of the variables' importance.The optimal θ = (θ opt 1 , • • • , θ opt m ) in the Kriging models are usually optimized by a continuous optimization algorithm.In this paper, the optimal θ is optimized by the simplex search method of Lagarias et al. (f minsearch) [26], with the parameter of max function evaluations equal to 1000.The covariance matrix can then be expressed through the correlation function: When µ(x) is assumed to be an unknown constant, the unbiased prediction is called ordinary Kriging (OK).In OK, the Kriging model determines the hyperparameters θ = [θ 1 , θ 2 , • • • , θ n ] by maximizing the likelihood over the observed dataset.The expression of the likelihood function is: The maximum likelihood estimates of the mean μ and the variance σ2 are given by: Then the predictor of the mean and the variance at a target point x t can be derived.They are shown in [21]:

Structure of MOBGO
In MOBGO, it is assumed that d objective functions are mutually independent in an objective space.Each objective function is approximated by a Kriging model individually, based on the η existing evaluated data D = (x (1) , y (1) = y(x (1) )), . . ., (x (η) , y (η) = y(x (η) )) .Each Kriging model is a one-dimensional normal distribution, with a mean µ and a standard deviation σ.Given a target solution x t , the Kriging models can predict the multivariate outputs by means of an independent joint normal distribution with means µ 1 , . . ., µ d and standard deviations σ 1 , . . ., σ d .These predictive means and standard deviations are used to calculate the score of an infill criterion, which can quantitatively measure how promising the target point x t is when compared with the current Pareto-front approximation set.A promising solution x * can be found by maximizing/minimizing3 the score of the infill criterion.Then, this promising solution x * is evaluated by the 'true' objective functions, and both the dataset D and the Pareto-front approximation set P are updated.The basic structure of the MOBGO algorithm is shown in Algorithm 1.It mainly contains three parts: initialization of a sampling dataset, searching for an optimal solution and updating the Kriging models, and returning the Pareto-front approximation set P.
The second part of MOBGO is the main loop, as shown in Algorithm 1 from Step 6 to Step 12.This main loop starts with training Kriging models M i based on dataset D (Step 7).Note that M contains d independent models for each objective function, and these models will be used as temporary objective functions instead of 'true' objective functions in Step 8.Then, an optimizer finds a promising solution x * by maximizing or minimizing an infill criterion C (Step 8).Here, an infill criterion is calculated by its corresponding calculation formula, whose inputs include Kriging models M , the current Pareto-front approximation set P, a target decision vector x t , etc. Theoretically, any single-objective optimization algorithm can be utilized as an optimizer to search for a promising solution x * .In this paper, the BI-population CMA-ES is chosen for its favorable performance on BBOB function testbed [18].Step 9 and Step 10 will update the dataset D by adding (x * , y(x * )) into D and update the Pareto-front approximation set P. The main loop from Step 6 to Step 12 will continue until g meets the termination criterion T c .The last part of MOBGO returns Pareto-front approximation set P.
The choice of infill criterion C at Step 8 distinguishes different types of MOBGO based algorithms.In this paper, EHVI-MOBGO and PoI-MOBGO, which set EHVI and PoI [22,24,35] as the infill criterion C respectively, are compared in Section 5.2.
From the perspectives of searching and optimization, non-dominated points are of greater interest.The concept of non-dominance is defined as: Definition 2 (Non-dominance [14]) Given a decision vector set X ∈ R m , and the image of the vector set Y = {y(x)|x ∈ X}, the non-dominated subset of Y is defined as: A vector y ∈ nd(Y) is called a non-dominated point of Y.
Definition 3 (Dominated Subspace of a Set) Let P be a subset of R d .The dominated subspace of P in R d , notation dom (P), is then defined as: Definition 4 (Non-Dominated Space of a Set) Let P be a subset of R d and let r ∈ R d be such that ∀p ∈ P : p ≺ r.The non-dominated space of P with respect to r, denoted as ndom(P), is then defined as: Note that the notion of dominated space as well as the notion of non-dominated space of a set can also be defined for (countably and non-countably) infinite sets P.
The Hypervolume Indicator (HV), introduced in [42], is one of the essential unary indicators for evaluating the quality of a Pareto-front approximation set.Its theoretical properties are discussed in [43].Notably, HV does not require the knowledge of the Pareto front in advance.The maximization of HV leads to a high-qualified and diverse Pareto-front approximation set.The Hypervolume Indicator is defined as follows: Definition 5 (Hypervolume Indicator) Given a finite approximation to a Pareto front, say P = {y (1) , . . ., y (n) } ⊂ R d , the Hypervolume Indicator of P is defined as the ddimensional Lebesgue measure of the subspace dominated by P and bounded below by a reference point r: with λ d being the Lebesgue measure on R d .
The hypervolume indicator measures the size of the dominated subspace bounded below by a reference point r.This reference point needs to be provided by users.Theoretically, in order to get the extreme non-dominated points, this reference point should be chosen in a way that it is dominated by all elements of a Pareto-front approximation set P during the optimization process.However, there is no requirement of setting the reference point in practice if the user is not interested in extreme non-dominated points.
Another important infill criterion is Hypervolume Improvement, which is also called the Improvement of Hypervolume in [11].The definition of Hypervolume Improvement is: Definition 6 (Hypervolume Improvement) Given a finite collection of vectors P ⊂ R d , the Hypervolume Improvement of a vector y ∈ R d is defined as: HVI(y, P) = HV(P ∪ {y}) − HV(P) (15) When we want to emphasize the reference point r, the notation HVI(y, P, r) will be used to denote Hypervolume Improvement.
Example 1 Figure 1 illustrates the concept of Hypervolume Improvement using two examples.The first example, on the left, is a 2-D example: Suppose a Pareto-front approximation set is P, which is composed by y (1) = (1, 2.5) , y (2) = (2, 1.5) and y (3) = (3, 1) .When a new point y (+) = (2.8,2.3) is added, the Hypervolume Improvement HV I(y (+) , P, r) is the area of the yellow polygon.The second example (on the right in Figure 1) illustrates the Hypervolume Improvement by means of a 3-D example.Assume a Pareto-front approximation set is P =( y (1) = (4, 4, 1) , The Hypervolume Improvement of y (+) = (3, 3, 2) relative to P is given by the joint volume covered by the yellow slices.
y (1)   y (2)   y Probability of Improvement (PoI) is an important criterion in MOBGO.It was first introduced by Stuckman in [34].Later, Emmerich et al. [13] generalized it to multi-objective optimization.PoI is defined as: Definition 7 (Probability of Improvement) Given parameters of the multivariate predictive distribution µ, σ and the Pareto-front approximation set P, the Probability of Improvement is defined as: where ξ µ,σ is a multivariate independent normal distribution with the mean values µ ∈ R d and the standard deviations σ ∈ R d + .Here, (y impr P) represents y ∈ R d as an improvement with respect to P, if and only if the following holds: y ≺ r and ∀p ∈ P : ¬(p ≺ y).
In Equation ( 16), I(y impr P) = 1 means that y is an element of the non-dominated space of P. In other words, y ∈ [r, ∞ d ] \ dom (P) if I(y impr P) = 1.A reference point r is not indicated in Equation ( 16) because r must be chosen as [−∞] d in PoI.Therefore, PoI is a reference-free infill criterion.Definition 8 (Expected Hypervolume Improvement) Given parameters of the multivariate predictive distribution µ, σ and the Pareto-front approximation set P, the expected hypervolume improvement is defined as: Example 2 An illustration of the 2-D EHVI is shown in Figure 2. The light gray area is the dominated subspace of P = {y (1)  y (2)   y (3)   y (+) HV I(y (+) ,P,r) For computing integrals of EHVI in Section 4, it is useful to define ∆ and Ψ ∞ functions.
Definition 9 (∆ function (see also [14,39])) For a given vector of objective function values y ∈ R d and y ∈ P, ∆(y, P, r) is the subset of the vectors in R d which are exclusively dominated by the vector y but not by elements in P, and which dominate the reference point r, that is: Definition 10 (Ψ ∞ function (see also [19]) denote its cumulative probability distribution function (CDF), and erf denote the Gaussian error function.The general normal distribution with mean µ and standard deviation σ has PDF ξ µ,σ (s One can easily show that

Efficient EHVI Calculation
This section mainly discusses an efficient partitioning method for a non-dominated space and how to employ this partitioning method to calculate EHVI and PoI.

Partitioning a non-dominated space
The efficiency of an infill criterion calculation is determined by a non-dominated search algorithm and the number of integration slices.The main idea of the partitioning method is to separate the integration volume (a non-dominated space) into as few integration slices as possible.Then, the integral of the criterion is calculated within each integration slice.The value of the criterion is the sum of its contribution in every integration slice.

The 2-D case
In the 2-D case, the partitioning method is simple and has already been published by Emmerich et al. [14].Given a Pareto-front approximation set P with n elements, the algorithm in [14] adopts a new way to derive the EHVI calculation formulas and only partitions a nondominated space into n + 1 integration slices, instead of (n + 1) 2 grids in [5,19].For the sake of completeness, we will introduce this integration technique here briefly.Suppose Y = y (1) , . . ., y (n) and d = 2, an integration space (a non-dominated space) of Y can be divided into n + 1 disjoint integration slices (S (i) 2 , i = 1, . . ., n + 1) by drawing lines parallel to y 2 -axis at each element in y, as indicated in Figure 3.Then, each integration slice can be expressed by its lower bound (l 2 ) and upper bound (u 2 ).In order to define the slices formally, we argue a Pareto-front approximation set P with two sentinels: . Then, the integration slices for the 2-D case are defined by: y ( 2) (3) 2
In the 2-D case, the number of integration slices is straightforward, namely, N 2 = n+1.

The 3-D case
Similar to the 2-D partitioning method, in the 3-D case, each integration slice can be defined by its lower bound (l 3 ) and upper bound (u 3 ).Since the upper bound of each integration slice is always ∞ in the y 3 axis, we can describe each integration slice as follows: Example 3 An illustration of integration slices is shown in with the corresponding lower and upper bounds (l 3 ).The partitioning algorithm is similar to the sweep line algorithm described in [15].The basic idea of this algorithm is to use an AVL tree to process points in descending order of the y 3 coordinate.For each such point, say Lower center: The projection of 3-D integration slices into the y 1 y 2 -plane, each slice can be described by lower bound and upper bound. ), the algorithm finds all the points (y (d [1]) , . . ., y (d[s]) ) which are dominated by y (i) in the y 1 y 2 -plane and inserts y (i) into the tree.Moreover, because of y (i) , the algorithm will also discard all the points (y (d [1]) , . . ., y (d[s]) ) from the AVL tree.See Figure 5 for describing one such iteration.In each iteration, s + 1 slices are created by coordinates of the points y (t) , y (d [1]) , . . ., y (d[s]) , y (r) , and y (i) as illustrated in Figure 5.
The number of the integration slices in the 3-D case N 3 is 2n + 1 where all points are in general position (for each i, i = 1, . . ., d: the i-th coordinate is different for each pair of points in Y).Otherwise, 2n + 1 provides an upper bound for the obtained number of slices.
Proof: In the algorithm, each point y (i) | i=1,...,n creates two slices.The first one, say slice A (i) , is created when the point y (i) is added to the AVL tree.Another slice, say slice 3 , is created when the point y (i) is discarded from the AVL tree due to domination by another point, say y (s) , in the y 1 y 2 -plane.These two slices are defined as follows 2 , ∞) ) and y (u) denote either the right neighbor among the newly dominated points in the y 1 y 2 -plane, or y (s) if y (i) is the rightmost point among all newly dominated points.In this way, each slice can be attributed to exactly one point in P, except the slice that is created in the final iteration.In the final iteration, one additional point y (n+1) = Algorithm 2: Integration slice acquisition in 3-D case Input: (y (1) , • • • , y (n) ): mutually non-dominated R 3 -points sorted by third coordinate (y 3 ) in descending order Output: .l 3 = y 13 .l 1 = y ) from tree T; 24 Insert y (i) in tree T. 25 Return S (1) (∞, ∞, ∞) is added to the AVL tree.This point will create a new slice when it is added, but because it is never discarded, it adds only a single slice.Therefore, 2n + 1 slices are created in total.

Higher dimensional cases
In higher dimensional cases, the non-dominated space can be partitioned into axis aligned hyperboxes, similar to the 3-D case.In the d-dimensional case (d ≥ 4), the hyperboxes can be denoted by S with their lower bounds (l (1) , . . ., l (N d ) ) and upper bounds (u (1) , . . ., u (N d ) ).Here, N d is the number of hyperboxes and has the same definition as N 2 and N 3 .The hyper-integral box S (i) d is defined as: An efficient algorithm for partitioning a higher dimensional, non-dominated space is proposed in this section, which is based on two state-of-the-art algorithms DKLV17 [6] by In other words, LKF17 is also efficient in partitioning the dominated space and provides the boundary information for each hyperbox in the dominated space.
The idea behind the proposed algorithm is transforming the problem of partitioning a non-dominated space into the problem of partitioning the dominated space, by means of introducing an intermediate Pareto-front approximation set P .This transformation is done by the following steps.Suppose that we have a current Pareto-front approximation set P for a maximization problem and we want to partition the non-dominated space of P. Firstly, DKLV17 is applied to locate the local lower bound points (L d ) of P in the dominated space.Secondly, regard L d as a new Pareto-front approximation P for a minimization problem with a reference point {∞} d .The dominated space of P is actually the non-dominated space of P.Then, LKF17 can be applied to partition the dominated space of P by locating the lower bound points l d and the upper bound points u d .These bound points (l d ,u d ) of P in the dominated space for a minimization problem are exact the lower/upper bound points of the partitioned, non-dominated hyperboxes of P for a maximization problem.The pseudo code of partitioning non-dominated space in higher dimensional cases is shown in Algorithm 3.
Set a new reference point r = (∞, ∞) and utilize LKF17 to partition the dominated space of P , by considering minimization, see Figure 6 (below left).The partitioned nondominated space of P is then the partitioned dominated space of P , see Figure 6 (below right).

EHVI calculation 6
This section discusses the problem of exact EHVI calculation.Moreover, a new and efficient algorithm is also derived.Section 4.2.1 and 4.2.2 introduce the proposed method in the 2-D and 3-D cases, respectively.Section 4.2.3 illustrates the general calculation formulas in higher dimensional cases, based on the proposed method.
In order to simplify the notation, ∆(y) is used whenever P, r are given by the context.Based on ∆(y), the expected hypervolume improvement function can be re-defined as: For the convenience of expressing the EHVI formula in the remaining parts of this paper, two functions ( and ϑ) are defined as follows: Definition 11 ( function) Given the parameters of an integration slice S (i) 6 Both C++ and MATLAB source code for computing the EHVI are available on http://liacs.leidenuniv.nl/˜csmoda/index.php?page=code or on request from the authors.
Fig. 6: The illustration of partitioning a space in higher dimensional cases.Above left: Pareto-front approximation set P. Above right: Locating L 2 points using DKLV17.Below left: Partitioning the dominated space of P using LKF17.Below right: The partitioned non-dominated space of P.
Definition 12 (ϑ function) Given the parameters of an integration slice S (i) d in a d-dimensional space and multivariate predictive distribution µ, σ, the function ϑ(l

2-D EHVI calculation
According to the definition of the 2-D integration slice in Equation 4.1.1,the Hypervolume Improvement y ∈ R 2 in the 2-D case is: HVI 2 gives rise to the compact integral for the original EHVI: Here y = (y 1 , y 2 ), the intersection of 2 with ∆(y 1 , y 2 ) is non-empty if and only if (y) dominates the lower left corner of S (i) 2 .Therefore: In Equation ( 28), the summation is done after integration.This operation is allowed, because integration is a linear mapping.Moreover, the integration interval ), because the HVI in one dimension λ 1 [S (i) 2 ∩ ∆(y 1 )] differs in these two integration intervals.Here λ 1 [B i ∩ ∆(y k )] is the HVI in dimension k, i.e., a 1-D HVI.Equation ( 28) can then be expressed as: According to the definition of HVI, (u ) is constant and is (u . Therefore, the Expected Improvement in dimension y 1 is also a constant and it is: ϑ(l Recall the Ψ ∞ function, by which the terms ( 29) and ( 30) can be expressed as follows: According to Equation ( 28), the exact EHVI calculation needs to compute the terms ( 29) and ( 30) n + 1 times, and each calculation requests O(1) computation.To keep P sorted in the first coordinate requires an effort of amortized time complexity O(log n) per iteration.Hence, the time complexity of the expected hypervolume improvement in the 2-D case is in O(n log n).In the case when P is sorted, we can show that the time complexity is in Θ(n).

3-D EHVI calculation
Given a partitioning of the non-dominated space into integration slices S , the EHVI integrations over each slice can be computed separately.To see how this calculation can be done, the Hypervolume Improvement of a point y ∈ R 3 is rewritten as: where ∆(y) is the part of the objective space that is dominated by y.The HVI expression in the definition of EHVI in Equation ( 17) can be replaced by HVI 3 in Equation ( 33): ), respectively.Also, again we can swap integration and summation based on the fact that integration is a linear mapping.Based on this subdivision, Equation ( 34) can be expressed as: Recalling the definition of the ϑ function and calculation of λ 1 [B i ∩ ∆(y k )], the term (35) can be rewritten as follows: ϑ(l ϑ(l ) ) ) The EHVI value for a hyperbox S (i) d , denoted as EHV I An integration for a divided interval.
Here, N d is the number of integration slices, and N 2 = n + 1, N 3 = 2n + 1 in the 2-D and 3-D case respectively.Since PoI is a reference-free indicator, a reference point r = {−∞} d is only used in order to obtain the correct boundary information (l d , u d ).

Speed Comparison
The test benchmarks from Emmerich and Fonseca [15] were used to generate Pareto-front sets.The Pareto-front sets and evaluated points were randomly generated based on CON-VEXSPHERICAL and CONCAVESPHERICAL functions.Two EHVI calculation algorithms, CDD13 [5] and KMAC, were compared using the same benchmarks in this experiment. 7 Note that KMAC algorithm in this paper includes the KMAC 2D for 2-D EHVI calculation in [14], the KMAC 3D algorithm for 3-D EHVI calculation in [38], and the extended KMAC algorithm for higher dimensional cases (d ≥ 4) in this paper.The experimental results in Figure 7 show that KMAC is much faster than CDD13, especially when |P| is increased.Moreover, CDD13 can not calculate the exact EHVI value within 30 minutes when |P| is bigger than 30 in the 5-D case.

Benchmark Performance
Five state-of-the-art algorithms are compared in this section, namely: EHVI-MOBGO, PoI-MOBGO 8 , NSGA-II [8], NSGA-III [7,17] and SMS-EMOA [1].The benchmarks are DTLZ1, DTLZ2, DTLZ3, DTLZ4, DTLZ5, DTLZ7 [9], MaF1, MaF5, MaF12 and MaF13 [2].The dimension m of all the benchmarks is 6 or 18.The parameter settings for all of these algorithms are shown in Table 2.The number of function evaluations (T c in Algorithm 1) for MOBGO based algorithms is 300.The Reference points for each benchmark are shown in Table 3.For DTLZ1 and DTLZ2, the reference points are chosen from the article [5].For the other test problems, the reference points are revised from the articles [2,5].All experiments were repeated 10 times.
The final Pareto fronts are evaluated by means of the Hypervolume indicator.Table 4 and Table 5 show the experimental results (statistical means and standard deviations) in 6-and 18-dimensional search spaces, respectively.Compared with EMOAs in 6-and 18dimensional search spaces, either EHVI-MOBGO or PoI-MOBGO yields the best result on the 10 benchmark function with 300 function evaluations.On the test problems of DTLZ7 7 Another EHVI calculation algorithm, IRS fast, is not compared in this paper, as it only works when d = 2, 3.For the detailed speed comparisons between IRS fast and KMAC in the 3-D case, see [38]. 8Using EHVI or PoI as the infill criterion in Algorithm 1.    generalizes them to higher dimensional cases with d ≥ 4. By using the fast box decomposition techniques, which were recently developed by Dächert et al. [6] and Lacour et al. [25], a non-dominated space can be partitioned with only Θ(n d/2 ) hyperboxes.The time complexity of our new EHVI computation algorithm scales linearly with the number of hyperboxes of arbitrary box decompositions.Unlike previous EHVI computation algorithms, the new algorithm does not require full grid partitionings with Θ(n d ) boxes.The new algorithm is, therefore, a significant improvement in terms of asymptotic time complexity.In addition, our benchmarks on random non-dominated sets show that this new algorithm is also many orders of magnitude faster for computations of typically sized problems where n ≤ 1000.
This paper also compares the performance of MOBGO based algorithms with three other state-of-the-art EMOAs on 10 benchmark test problems.For budgets of function evaluations up to 300, MOBGO based algorithms can yield the Pareto-front approximation sets with higher HV values than that of EMOAs in both 6-and 18-dimensional search spaces.In most cases, EHVI-MOBGO yields better performance than PoI-MOBGO.However, PoI-MOBGO performs better than EHVI-MOBGO on up to 3 test problems in 6-and 18-dimensional search spaces.The reason is that the PoI can imply an improvement of an evaluated point in the whole non-dominated space, but the EHVI can only indicate the improvement in the non-dominated space bounded by a reference point.A remedy to the EHVI's disadvantage can be achieved by setting a large reference point or using the dynamic reference point strategy.However, the selected reference point must not be too large.Otherwise, EHVI at any evaluated points would be similar, even identical, because of the insufficient numerical stability of the computations involved.
For future research, it is recommended to further investigate on reference-free computation of EHVI.Moreover, it is still an open question of how to obtain fast EHVI calculations for a larger number of objective functions.Although it is conjectured that in the worst case time complexity will increase superpolynomially with the number of objectives, a better average case time complexity could be obtained by adopting concepts from recently proposed divide-and-conquer algorithms for computing the hypervolume indicator [20,31].

Fig. 1 :
Fig. 1: The left and right figures illustrate Hypervolume Improvement in a 2-D and a 3-D example, respectively.
5) } bounded by the reference point r = (0, 0).The bivariate Gaussian distribution has the parameters µ 1 = 2.5, µ 2 = 2, σ 1 = 0.7, σ 2 = 0.8.The probability density function (ξ) of the bivariate Gaussian distribution is indicated as a 3-D plot.Here y (+) is a sample from this distribution and the area of improvement relative to P is indicated by the dark shaded area.Variables y 1 and y 2 stand for the first and the second objective values, respectively.

Figure 4 . 4 ) 3 = ( 1 , 2 , 2 )
A Pareto front set P is composed by 4 points (y (1) = (1, 3, 4) , y (2) = (4, 2, 3) , y (3) = (2, 4, 2) and y (4) = (3, 5, 1) ), and this Pareto front is shown in the upper left figure.The upper right figure represents the partitioned integration slices of P. The lower center figure illustrates the projection of the upper right figure onto the y 1 y 2 -plane with rectangle slices and l, u.The rectangular slices, which share a similar color but of different opacity, represent integration slices with the same value of y 3 in their lower bound.The lower bound of the 3-D integration slice B 4 is l (, and the upper bound of the slice is u

4 B 5 B 6 B 7 B 8 B 9 Fig. 4 :
Fig. 4: Upper left: 3-D Pareto-front approximation.Upper right: Integration slices in 3-D.Lower center: The projection of 3-D integration slices into the y 1 y 2 -plane, each slice can be described by lower bound and upper bound.

2 = 5
For the definition of the local lower bound points L d , see[6].

Algorithm 3 : 1 ) 2 ,
Partitioning a non-dominated space in higher dimensional cases Input: Pareto-front approximation set P (maximization problem), a reference point r Output: Hyperboxes S d 1 Locate local lower bound points L d : L d = DKLV 17(P, r); 2 Set new Pareto front P using L d : P = L d ; 3 Set a new reference point r : r = {∞} d ; 4 Get lower bound points l d and upper bound points u d : (l d , u d ) = LKF 17(P , r ) ; 5 S d = (l d , u d ) ; 6 Return S d (2, 1) and L (4) 2 = (3, 0) , see Figure 6 (above right).Regard all of the local lower bound points L 2 as the elements of a new Pareto-front approximation set P = (L (• • • , L The parameters: σ d = 2.5, µ d = 10, d = 2, • • • , 5 are used in the experiments.Pareto front sizes are |P| ∈ {10, 20, • • • , 200} and Batch Size is 1, which represents the number of the evaluated points under the same Pareto-front approximation set.Ten P sets are randomly generated by the same parameters.Average runtimes over 100 repetitions (10 repetitions for 10 P sets) were computed.All the experiments were performed on the same computer, and the hardware is: Intel(R) Xeon(R) CPU I7 3770 3.40GHz, RAM 16GB.The operating system is Ubuntu 16.04 LTS (64 bit), the compiler of KMAC is g++ 4.9.2 with compiler flag -Ofast, and CDD13 is based on MATLAB 8.4.0.150421 (R2014b), 64 bit.The experiments are set to halt when the algorithms cannot finish the EHVI computation within 30 minutes.

Table 1
+Number of the non-dominated points in P y(1), . . ., y (n) R d • , d The Kriging model for the i-th objective function

Table 2 :
Algorithm Parameter Settings

Table 3 :
Reference Points ∞) d \ dom (P)9, which is bigger than that of the EHVI y ∈ (r, ∞ d )\ dom (P).Therefore, PoI performs better when searching for extreme non-dominated points.In other words, EHVI is a reference based infill criterion and it cannot indicate any improvement of an evaluated point in a discarded part of the non-dominated space, namely, y ∈ (−∞ d , r) \ dom (P).

Table 4 :
Empirical Comparisons w.r.t HV in the case of m = 6This paper describes an efficient algorithm for EHVI calculation.It reviews and benchmarks recently proposed asymptotically optimal algorithms with Θ(n log n) time complexity and

Table 5 :
Empirical Comparisons w.r.t HV in the case of m = 18