Efficient Computation of Expected Hypervolume Improvement Using Box Decomposition Algorithms

TL;DR

This paper introduces an efficient algorithm for exactly computing the Expected Hypervolume Improvement in multi-objective Bayesian optimization, significantly reducing computational complexity and enabling faster optimization processes.

Contribution

The paper presents a novel hypervolume decomposition algorithm that improves the computational complexity of EHVI calculation from quadratic and cubic to near-linear, applicable to higher dimensions.

Findings

Algorithm reduces EHVI computation time significantly.

Complexity improved to Θ(n log n) for two-objective problems.

Applicable to higher dimensions with hyperbox decomposition.

Abstract

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 (EMOAs). 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 common 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. In this paper, an efficient algorithm for the computation of the exact EHVI for a generic case is proposed. This efficient algorithm is based on partitioning the integration volume into a set of axis-parallel slices. Theoretically, the upper bound time complexities are improved from previously $O (n^2)$ and $O(n^3)$, for two- and three-objective problems respectively, to $\Theta(n\log n)$, which is asymptotically optimal. This article generalizes the scheme in higher dimensional case by utilizing a new hyperbox decomposition technique, which was proposed by D{\"a}chert et al, EJOR, 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).

Tables

Table 1. Table 1 : Notations

Symbol	Type	Description
$m$	${\mathbb{N}}^{+}$	Dimension of a search space
$d$	${\mathbb{N}}^{+}$	Dimension of an objective space
$𝝁$	${ℝ}^d$	Mean values of predictive distribution
$𝝈$	${\left({ℝ}_0^{+}\right)}^d$	Standard deviations of predictive distribution
$𝒫$	${\left({ℝ}^d\right)}^n$	A Pareto-front approximation set
$n$, $\left|𝒫\right|$	${\mathbb{N}}^{+}$	Number of the non-dominated points in $𝒫$
${𝐲}^{\left(1\right)},\dots ,{𝐲}^{\left(n\right)}$	${ℝ}^d$	The vectors in $𝒫$, where $𝒫=({𝐲}^{\left(\mathrm{𝟏}\right)},\dots ,{𝐲}^{\left(𝐧\right)})$
$𝐫$	${ℝ}^d$	Reference point
$S_d$	${\left({ℝ}^d\right)}^2$	Integration slices in a $d$-dimensional space
$N_i$	${\mathbb{N}}^{+}$	Number of integration slice in $i$-dimensional case
${𝐥}_d^{\left(1\right)},\dots ,{𝐥}_d^{\left(N_d\right)}$	${ℝ}^d$	Lower bounds of integration slices
${𝐮}_d^{\left(1\right)},\dots ,{𝐮}_d^{\left(N_d\right)}$	${ℝ}^d$	Upper bounds of integration slices
${𝐱}^t$	${ℝ}^d$	A target solution in a search space
${𝐱}^{\ast }$	${ℝ}^d$	A promising solution in a search space
$M_i$	$i=1,\cdots ,d$	The Kriging model for the $i$-th objective function
$g$	${\mathbb{N}}^{+}$	Counter/Number of function evaluations
$T_c$	${\mathbb{N}}^{+}$	Termination criterion
$\eta$	${\mathbb{N}}^{+}$	The number of initialized solutions
$D$	${\left({ℝ}^{\left(d+m\right)}\right)}^g$	Training dataset for the Kriging models $M$
${\lambda }_i$	$i=1,\cdots ,d$	The Lebesgue measure on ${ℝ}^i$
$𝝃$		The multivariate independent normal distribution
$S_i$	$i=1,\cdots ,d$	The integration slices in an i-dimensional case
$y$		Objective functions

Table 2. Table 2 : Algorithm Parameter Settings

	EHVI-MOBGO	PoI-MOBGO	NSGA-II	NSGA-III	SMS-EMOA
$\mu$	30	30	30	/	30
$\lambda$	1	1	30		/
Evaluation	300	300	300/2000	300/2000	300/2000
divisions_outer	/	/	/	12	/
$p_c$	/	/	0.9		0.9
$p_m$	/	/	1/6		1/6
Platform	MATLAB	MATLAB	MATLAB	Python	MATLAB

Equations

$$\displaystyle``\max"\big{(}{y}_{1}(\mathbf{x}),{y}_{2}(\mathbf{x}),\cdots,{y}_{d}(\mathbf{x})\big{)}\qquad\qquad\mathbf{x}\in\mathbb{R}^{m}$$

$$Y(\mathbf{x})=\mu (\mathbf{x})+ϵ(\mathbf{x})$$

$$Corr[ϵ(\mathbf{x}),ϵ({\mathbf{x}}^{\prime })]=R(\mathbf{x},{\mathbf{x}}^{\prime })=\prod _{i=1}^m{R}_i(x_i,x_i^{\prime })$$

$$R(\mathbf{x},{\mathbf{x}}^{\prime })=\prod _{i=1}^m\text{exp}(-{\theta }_i(x_i-x_i^{\prime }{)}^2)({\theta }_i\ge 0)$$

$$Cov(ϵ(\mathbf{x}))={\sigma }^2\mathbf{\Sigma },\text{where}{\mathbf{\Sigma }}_{i,j}=R({\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{j}})$$

$$L=-\frac{n}{2}\ln ({\sigma }^2)-\frac{1}{2}\ln (|\mathbf{\Sigma }|)$$

$$\hat{\mu }=\frac{1_n^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{y}}{1_n^{\top }{\mathbf{\Sigma }}^{-1}{1}_n}$$

$${\hat{\sigma }}^2=\frac{1}{n}(\mathbf{y}-1_n\hat{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{y}-1_n\hat{\mu })$$

$$\mu ({\mathbf{x}}^t)=\hat{\mu }+{\mathbf{c}}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\hat{\mu }{1}_n)$$

$${\sigma }^2({\mathbf{x}}^t)={\hat{\sigma }}^2[1-{\mathbf{c}}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{c}+\frac{1-{\mathbf{c}}^T{\Sigma }^{\top }\mathbf{c}}{1_{\mathbf{n}}^{\top }{\mathbf{\Sigma }}^{-1}{1}_{\mathbf{n}}}]$$

$$nd(\mathbf{Y}):=\{\mathbf{y}\in \mathbf{Y}|\nexists \mathbf{z}\in \mathbf{Y}:\mathbf{z}\prec \mathbf{y}\}$$

$$\text{mbox}dom(\mathcal{P}):=\{\mathbf{y}\in {\mathbb{R}}^d|\exists \mathbf{p}\in \mathcal{P}\text{mbox}with\mathbf{p}\prec \mathbf{y}\}$$

$$\text{mbox}ndom(\mathcal{P}):=\{\mathbf{y}\in {\mathbb{R}}^d|\mathbf{y}\prec \mathbf{r}\text{mbox}and\nexists \mathbf{p}\in \mathcal{P}\text{mbox}suchthat\mathbf{p}\prec \mathbf{y}\}$$

$$\text{mbox}HV(\mathcal{P})={\lambda }_d({\cup }_{\mathbf{y}\in \mathcal{P}}[\mathbf{r},\mathbf{y}])$$

$$\text{mbox}HVI(\mathbf{y},\mathcal{P})=\text{mbox}HV(\mathcal{P}\cup \{\mathbf{y}\})-\text{mbox}HV(\mathcal{P})$$

$$\text{mbox}PoI(\mathbf{\mu },\mathbf{\sigma },\mathcal{P}):$$

$$EHVI(\mathbf{\mu },\mathbf{\sigma },\mathcal{P},\mathbf{r}):={\int }_{{\mathbb{R}}^d}\text{mbox}HVI(\mathcal{P},\mathbf{y},\mathbf{r})\cdot {\mathbf{\xi }}_{\mathbf{\sigma },\mathbf{\mu }}(\mathbf{y})d\mathbf{y}$$

$$\Delta (\mathbf{y},\mathcal{P},\mathbf{r})=\{\mathbf{z}\in {\mathbb{R}}^d|\mathbf{y}\prec \mathbf{z}\text{mbox}and\mathbf{z}\prec \mathbf{r}\text{mbox}and\nexists \mathbf{q}\in \mathcal{P}:\mathbf{q}\prec \mathbf{z}\}$$

$${\Psi }_{\infty }(a,b,\mu ,\sigma ):$$

$$\displaystyle S_{2}^{(i)}=(\mathbf{l}_{2}^{(i)},\mathbf{u}_{2}^{(i)})=\left(\left(\begin{array}[]{c}l_{1}^{(i)}\\ l_{2}^{(i)}\end{array}\right),\left(\begin{array}[]{c}u_{1}^{(i)}\\ u_{2}^{(i)}\end{array}\right)\right)=\left(\left(\begin{array}[]{c}y_{1}^{(i-1)}\\ y_{2}^{(i)}\end{array}\right),\left(\begin{array}[]{c}y_{1}^{(i)}\\ \infty\end{array}\right)\right),i=1,\dots,N_{2}$$

$$\displaystyle S_{3}^{(i)}=(\mathbf{l}_{3}^{(i)},\mathbf{u}_{3}^{(i)})=\left(\left(\begin{array}[]{c}l_{1}^{(i)}\\ l_{2}^{(i)}\\ l_{3}^{(i)}\end{array}\right),\left(\begin{array}[]{c}u_{1}^{(i)}\\ u_{2}^{(i)}\\ \infty\end{array}\right)\right),\quad\quad i=1,\dots,N_{3}$$

$$\displaystyle S_{d}^{(i)}=(\mathbf{l}_{d}^{(i)},\mathbf{u}_{d}^{(i)})=\big{(}(l_{1}^{(i)},\cdots,l_{d}^{(i)})^{\top},(u_{1}^{(i)},\cdots,\infty)^{\top}\big{)}\quad\quad i=1,\dots,N_{d}$$

$$\text{mbox}EHVI(\mathbf{\mu },\mathbf{\sigma },\mathcal{P},\mathbf{r})$$

$$={\int }_{y_1=-\infty }^{\infty }\cdots {\int }_{y_d=-\infty }^{\infty }{\lambda }_d[S_d\cap \Delta (\mathbf{y})]{\mathbf{\xi }}_{\mathbf{\mu },\mathbf{\sigma }}(\mathbf{y})d\mathbf{y}$$

$$\ell (u_k^{(i)},y_k,l_k^{(i)}):={\lambda }_1[S_d^{(i)}\cap \Delta (y_k)]=|[l_k^{(i)},u_k^{(i)}]\cap [l_k^{(i)},y_k]|=\min \{u_k^{(i)},y_k\}-l_k^{(i)}$$

$$\vartheta (l_k^{(i)},u_k^{(i)},{\sigma }_k,{\mu }_k):$$

$$={\int }_{y_k=u_k^{(i)}}^{\infty }(u_k^{(i)}-l_k^{(i)})\cdot {\mathbf{\xi }}_{{\mu }_k,{\sigma }_k}(y_k)d{y}_k$$

$$\displaystyle=(u_{k}^{(i)}-l_{k}^{(i)})\cdot\left(1-\Phi\Big{(}\frac{u_{k}^{(i)}-\mu_{k}}{\sigma_{k}}\Big{)}\right)\qquad k=1,\cdots,d-1$$

$$\text{mbox}{HVI}_2(\mathbf{y},\mathcal{P},\mathbf{r})=\sum _{i=1}^{N_2}{\lambda }_2[S_2^{(i)}\cap \Delta (\mathbf{y})]$$

$$\text{mbox}EHVI(\mathbf{\mu },\mathbf{\sigma },\mathcal{P},\mathbf{r})=$$

$$\text{mbox}EHVI(\mathbf{\mu },\mathbf{\sigma },\mathcal{P},\mathbf{r})$$

$$\sum _{i=1}^{N_2}{\int }_{y_1=u_1^{(i)}}^{\infty }{\int }_{y_2=l_2^{(i)}}^{\infty }{\lambda }_2[S_2^{(i)}\cap \Delta (\mathbf{y})]\cdot {\mathbf{\xi }}_{\mathbf{\mu },\mathbf{\sigma }}(\mathbf{y})d\mathbf{y}$$

Full text

1 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) Mockus1978 . 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 kyang2015cec . According to the authors’ knowledge, BGO was used for the first time in the context of airfoil optimization in laniewski2010development , and then applied in the field of biogas plant controllers gaida2014dynamic , detection in water quality management zaefferer2013case , machine learning algorithm configuration koch2015efficientnoise , and structural design optimization shimoyama2013updating .

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. Mockus1978 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 jones1998efficient . In MOBGO, a commonly used criterion is Expected Hypervolume Improvement (EHVI), which is a straightforward generalization of the EI and was proposed by Emmerich emmerich2005single in 2005. Compared with other criteria, EHVI leads to an excellent convergence to – and coverage of – the true Pareto front couckuyt2014fast ; kyang2016cec . Nevertheless, the calculation of EHVI itself so far has been time-consuming shimoyama2013kriging ; zaefferer2013case ; wagner2010expected ; koch2015efficientnoise , even in the 2-D case222In this paper, 2-D and 3-D represent two and three objective functions, respectively.. 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 emmerich2005single ; emmerich2006single . 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 emmerich2011hypervolume , with time complexity of $O(n^{3}\log n)$. Here, $n$ is the number of non-dominated points in the archive. The EHVI calculation algorithm in emmerich2011hypervolume partitions an objective space into boxes and then calculates the EHVI by summing all the EHVI values of each box. Couckuyt et al. couckuyt2014fast 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 couckuyt2014fast was also practically much faster than those discussed in emmerich2011hypervolume , though a detailed complexity analysis was missing. Hupkens et al. hupkens2015faster reduced the time complexity to $O(n^{2})$ and $O(n^{3})$ in the 2-D and 3-D cases, respectively. The algorithms in hupkens2015faster improve the algorithms in emmerich2011hypervolume 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 hupkens2015faster further improve the practical efficiency of EHVI on test data in comparison to couckuyt2014fast . Recently, Emmerich et al. Michael2016book proposed an asymptotically optimal algorithm with time complexity of $\Theta(n\log n)$ in the bi-objective case. More recently, Yang et al. proposed an asymptotically optimal algorithm with time complexity $\Theta(n\log n)$ in the 3-D case yang2017computing . The algorithm, KMAC333KMAC stands for the authors’ given names of the EHVI exact calculation algorithm. in yang2017computing , 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.

2 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:

[TABLE]

where $d$ is the number of objective functions, ${y}_{i}(i=1,\cdots,d)$ are the objective functions, and a decision vector $\mathbf{x}$ is in an $m$-dimensional space.

2.1 Notations

The following table summarizes the notations used in this paper.

2.2 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 li2008metamodel . 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 $\mathbf{X}=\{\mathbf{x}^{(1)},\mathbf{x}^{(2)},\cdots,\mathbf{x}^{(n)}\},\mathbf{x}^{(i)}\in\mathbb{R}^{m}|_{i,=1,\cdots,n}$ in an $m$-dimensional search space, and associated function values \mathbf{Y}(\mathbf{X})=\big{(}{y}(\mathbf{x}^{(1)}),{y}(\mathbf{x}^{(2)}),\cdots,{y}(\mathbf{x}^{(n)})\big{)}^{\top}, Kriging assumes $\mathbf{Y}(\mathbf{X})$ to be a realization of a random process $Y$ of the following form chugh2017handling ; jones1998efficient :

[TABLE]

where $\mu(\mathbf{x})$ is the estimated mean value over all given sampled points, and $\epsilon(\mathbf{x})$ is a realization of a Gaussian process with zero mean and variance $\sigma^{2}$. The regression part $\mu(\mathbf{x})$ approximates the function $Y$ globally and the Gaussian process $\epsilon(\mathbf{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 ($\mathbf{x}$ and $\mathbf{x^{\prime}}$) is defined as:

[TABLE]

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

[TABLE]

where $\theta_{i}$ are parameters of the correlation model. They can be interpreted as a measurement of the variables’ importance. The optimal $\boldsymbol{\theta}=(\theta_{1}^{opt},\cdots,\theta_{m}^{opt})$ in the Kriging models are usually optimized by a continuous optimization algorithm. In this paper, the optimal $\boldsymbol{\theta}$ is optimized by the simplex search method of Lagarias et al. ($fminsearch$) lagarias1998fminsearch , with the parameter of max function evaluations equal to 1000. The covariance matrix can then be expressed through the correlation function:

[TABLE]

When $\mu(\mathbf{x})$ is assumed to be an unknown constant, the unbiased prediction is called ordinary Kriging (OK). In OK, the Kriging model determines the hyperparameters $\mathbf{\theta}=[\theta_{1},\theta_{2},\cdots,\theta_{n}]$ by maximizing the likelihood over the observed dataset. The expression of the likelihood function is:

[TABLE]

The maximum likelihood estimates of the mean $\hat{\mu}$ and the variance $\hat{\sigma}^{2}$ are given by:

[TABLE]

Then the predictor of the mean and the variance at a target point $\mathbf{x}^{t}$ can be derived. They are shown in jones1998efficient :

[TABLE]

where $\mathbf{c}=(Corr[y(x^{t}),y(x_{1})],\cdots,Corr[y(x^{t}),y(x_{n})])^{\top}$.

2.3 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 $\eta$ existing evaluated data D=\big{(}(\mathbf{x}^{(1)},\mathbf{y}^{(1)}=y(\mathbf{x}^{(1)})),\dots,(\mathbf{x}^{(\eta)},\mathbf{y}^{(\eta)}=y(\mathbf{x}^{(\eta)}))\big{)}. Each Kriging model is a one-dimensional normal distribution, with a mean $\mu$ and a standard deviation $\sigma$. Given a target solution $\mathbf{x}^{t}$, the Kriging models can predict the multivariate outputs by means of an independent joint normal distribution with means $\mu_{1}$, $\dots$, $\mu_{d}$ and standard deviations $\sigma_{1}$, $\dots$, $\sigma_{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 $\mathbf{x}^{t}$ is when compared with the current Pareto-front approximation set. A promising solution $\mathbf{x}^{*}$ can be found by maximizing/minimizing444It depends on which criterion is chosen. the score of the infill criterion. Then, this promising solution $\mathbf{x}^{*}$ is evaluated by the ’true’ objective functions, and both the dataset $D$ and the Pareto-front approximation set $\mathcal{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 $\mathcal{P}$.

First, a dataset $D$ is initialized and a Pareto-front approximation set $\mathcal{P}$ is computed, as shown in Algorithm 1 from Step 1 to Step 5. The initialization of $D$ contains the generation of the decision vectors using Latin Hypercube Sampling method (LHS) lhs1979 (Step 1), calculation of the corresponding objective values (Step 2) and storage of this information in dataset $D$ (Step 3). This dataset $D$ will be utilized to build the Kriging models in the second part.

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 $\mathbf{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 $\mathcal{P}$, a target decision vector $\mathbf{x}^{t}$, etc. Theoretically, any single-objective optimization algorithm can be utilized as an optimizer to search for a promising solution $\mathbf{x}^{*}$. In this paper, the BI-population CMA-ES is chosen for its favorable performance on BBOB function testbed hansen2009benchmarking . Step 9 and Step 10 will update the dataset $D$ by adding $(\mathbf{x}^{*},y(\mathbf{x}^{*}))$ into $D$ and update the Pareto-front approximation set $\mathcal{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 $\mathcal{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 kushner1964poi ; ulmer2003evolution ; keane2006statistical as the infill criterion C respectively, are compared in Section 5.2.

3 Definitions555For the convenience of visualization and consistency, this paper only considers maximization problems. Minimization problems can always be re-written as maximization problems by multiplying the corresponding objective functions by $(-1)$.

Pareto dominance, or briefly dominance, is a fundamental concept in MOO and provides an ordering relation on the set of potential solutions. Dominance is defined as follows:

Definition 1 (Dominance CoelloCoello2011 )

Given two decision vectors $\mathbf{x}^{(1)},\mathbf{x}^{(2)}\in\mathbb{R}^{m}$ and their corresponding objective values $\mathbf{y}^{(1)}=y(\mathbf{x}^{(1)})$, $\mathbf{y}^{(2)}=y(\mathbf{x}^{(2)})$ in a maximization problem, it is said that $\mathbf{y}^{(1)}$ dominates $\mathbf{y}^{(2)}$, being represented by $\mathbf{y}^{(1)}\prec\mathbf{y}^{(2)}$, iff  $\forall i\in\{1,2,\cdots,d\}:{y}_{i}(\mathbf{x}^{(1)})\geq{y}_{i}(\mathbf{x}^{(2)})$ and $\exists j\in\{1,2,\cdots,d\}:{y}_{j}(\mathbf{x}^{(1)})>{y}_{j}(\mathbf{x}^{(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 Michael2016book )

Given a decision vector set $\mathbf{X}\in\mathbb{R}^{m}$, and the image of the vector set $\mathbf{Y}=\{y(\mathbf{x})|\mathbf{x}\in\mathbf{X}\}$, the non-dominated subset of $\mathbf{Y}$ is defined as:

[TABLE]

A vector $\mathbf{y}\in nd(\mathbf{Y})$ is called a non-dominated point of $\mathbf{Y}$.

Definition 3 (Dominated Subspace of a Set)

Let $\mathcal{P}$ be a subset of $\mathbb{R}^{d}$. The dominated subspace of $\mathcal{P}$ in $\mathbb{R}^{d}$, notation $\mbox{ dom }(\mathcal{P})$, is then defined as:

[TABLE]

Definition 4 (Non-Dominated Space of a Set)

Let $\mathcal{P}$ be a subset of $\mathbb{R}^{d}$ and let $\mathbf{r}\in\mathbb{R}^{d}$ be such that $\forall\mathbf{p}\in\mathcal{P}:\mathbf{p}\prec\mathbf{r}$. The non-dominated space of $\mathcal{P}$ with respect to $\mathbf{r}$, denoted as $\mbox{ndom}(\mathcal{P})$, is then defined as:

[TABLE]

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 $\mathcal{P}$.

The Hypervolume Indicator (HV), introduced in zitzler1999multiobjective , is one of the essential unary indicators for evaluating the quality of a Pareto-front approximation set. Its theoretical properties are discussed in zitzler2003performance . 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 $\mathcal{P}=\{\mathbf{y}^{(1)},\dots,\mathbf{y}^{(n)}\}\subset\mathbb{R}^{d}$, the Hypervolume Indicator of $\mathcal{P}$ is defined as the $d$-dimensional Lebesgue measure of the subspace dominated by $\mathcal{P}$ and bounded below by a reference point $\mathbf{r}$:

[TABLE]

with $\lambda_{d}$ being the Lebesgue measure on $\mathbb{R}^{d}$.

The hypervolume indicator measures the size of the dominated subspace bounded below by a reference point $\mathbf{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 $\mathcal{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 emmerich2008computation . The definition of Hypervolume Improvement is:

Definition 6 (Hypervolume Improvement)

Given a finite collection of vectors $\mathcal{P}\subset\mathbb{R}^{d}$, the Hypervolume Improvement of a vector $\mathbf{y}\in\mathbb{R}^{d}$ is defined as:

[TABLE]

When we want to emphasize the reference point $\mathbf{r}$, the notation $\mbox{HVI}(\mathbf{y},\mathcal{P},\mathbf{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 $\mathcal{P}$, which is composed by $\mathbf{y}^{(1)}=(1,2.5)^{\top}$, $\mathbf{y}^{(2)}=(2,1.5)^{\top}$ and $\mathbf{y}^{(3)}=(3,1)^{\top}$. When a new point $\mathbf{y^{(+)}}=(2.8,2.3)^{\top}$ is added, the Hypervolume Improvement $HVI(\mathbf{y^{(+)}},\mathcal{P},\mathbf{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 $\mathcal{P}=$( $\mathbf{y}^{(1)}=(4,4,1)^{\top}$, $\mathbf{y}^{(2)}=(1,2,4)^{\top}$, $\mathbf{y}^{(3)}=(2,1,3)^{\top}$). The Hypervolume Improvement of $\mathbf{y}^{(+)}=(3,3,2)^{\top}$ relative to $\mathcal{P}$ is given by the joint volume covered by the yellow slices.

Probability of Improvement (PoI) is an important criterion in MOBGO. It was first introduced by Stuckman in stuckman1988global . Later, Emmerich et al. emmerich2006single generalized it to multi-objective optimization. PoI is defined as:

Definition 7 (Probability of Improvement)

Given parameters of the multivariate predictive distribution $\boldsymbol{\mu}$, $\boldsymbol{\sigma}$ and the Pareto-front approximation set $\mathcal{P}$, the Probability of Improvement is defined as:

[TABLE]

where $\boldsymbol{\xi}_{\boldsymbol{\mu},\boldsymbol{\sigma}}$ is a multivariate independent normal distribution with the mean values $\boldsymbol{\mu}\in\mathbb{R}^{d}$ and the standard deviations $\boldsymbol{\sigma}\in\mathbb{R}^{d}_{+}$. Here, ($\mathbf{y}\mbox{ impr }\mathcal{P}$) represents $\mathbf{y}\in\mathbb{R}^{d}$ as an improvement with respect to $\mathcal{P}$, if and only if the following holds: $\mathbf{y}\prec\mathbf{r}$ and $\forall\mathbf{p}\in\mathcal{P}:\neg(\mathbf{p}\prec\mathbf{y})$.

In Equation (16), $\mbox{I}(\mathbf{y}\mbox{ impr }\mathcal{P})=1$ means that $\mathbf{y}$ is an element of the non-dominated space of $\mathcal{P}$. In other words, $\mathbf{y}\in[\mathbf{r},\infty^{d}]\setminus\mbox{ dom }(\mathcal{P})$ if $\mbox{I}(\mathbf{y}\mbox{ impr }\mathcal{P})=1$. A reference point $\mathbf{r}$ is not indicated in Equation (16) because $\mathbf{r}$ must be chosen as $[-\infty]^{d}$ in PoI. Therefore, PoI is a reference-free infill criterion.

Definition 8 (Expected Hypervolume Improvement)

Given parameters of the multivariate predictive distribution $\boldsymbol{\mu}$, $\boldsymbol{\sigma}$ and the Pareto-front approximation set $\mathcal{P}$, the expected hypervolume improvement is defined as:

[TABLE]

Example 2

An illustration of the 2-D EHVI is shown in Figure 2. The light gray area is the dominated subspace of $\mathcal{P}=\{{\mathbf{y}^{(1)}}=(3,1)^{\top},$ $\mathbf{y}^{(2)}=(2,1.5)^{\top},$ ${\mathbf{y}^{(3)}}=(1,2.5)^{\top}\}$ bounded by the reference point $\mathbf{r}=(0,0)$. The bivariate Gaussian distribution has the parameters $\mu_{1}=2.5,$ $\mu_{2}=2,$ $\sigma_{1}=0.7,$ $\sigma_{2}=0.8$. The probability density function ($\boldsymbol{\xi}$) of the bivariate Gaussian distribution is indicated as a 3-D plot. Here $\mathbf{y}^{(+)}$ is a sample from this distribution and the area of improvement relative to $\mathcal{P}$ is indicated by the dark shaded area. Variables $y_{1}$ and $y_{2}$ stand for the first and the second objective values, respectively.

For computing integrals of EHVI in Section 4, it is useful to define $\Delta$ and $\Psi_{\infty}$ functions.

Definition 9 ($\Delta$ function (see also Michael2016book ; yang2017phdthesis ))

For a given vector of objective function values $\mathbf{y}\in\mathbb{R}^{d}$ and $\mathbf{y}\not\in\mathcal{P}$, $\Delta(\mathbf{y},\mathcal{P},\mathrm{r})$ is the subset of the vectors in $\mathbb{R}^{d}$ which are exclusively dominated by the vector $\mathbf{y}$ but not by elements in $\mathcal{P}$, and which dominate the reference point $\mathbf{r}$, that is:

[TABLE]

Definition 10 ($\Psi_{\infty}$ function (see also hupkens2015faster ))

Let $\phi(s)=1/\sqrt{2\pi}e^{-\frac{1}{2}s^{2}}(s\in\mathbb{R})$ denote the PDF ($\xi$) of the standard normal distribution. Moreover, let $\Phi(s)=\frac{1}{2}\left(1+\operatorname{erf}\left(\frac{s}{\sqrt{2}}\right)\right)$ denote its cumulative probability distribution function (CDF), and $\operatorname{erf}$ denote the Gaussian error function. The general normal distribution with mean $\mu$ and standard deviation $\sigma$ has PDF $\xi_{\mu,\sigma}(s)=\phi_{\mu,\sigma}(s)=\frac{1}{\sigma}\phi(\frac{s-\mu}{\sigma})$ and its CDF is $\Phi_{\mu,\sigma}(s)=\Phi(\frac{s-\mu}{\sigma})$. Then the function $\Psi_{\infty}(a,b,\mu,\sigma)$ is defined as:

[TABLE]

One can easily show that $\psi_{\infty}(a,b,\mu,\sigma)=\sigma\phi\left(\dfrac{b-\mu}{\sigma}\right)+(\mu-a)\left[1-\Phi\left(\dfrac{b-\mu}{\sigma}\right)\right]$.

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

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

4.1.1 The 2-D case

In the 2-D case, the partitioning method is simple and has already been published by Emmerich et al. Michael2016book . Given a Pareto-front approximation set $\mathcal{P}$ with $n$ elements, the algorithm in Michael2016book adopts a new way to derive the EHVI calculation formulas and only partitions a non-dominated space into $n+1$ integration slices, instead of $(n+1)^{2}$ grids in hupkens2015faster ; couckuyt2014fast . For the sake of completeness, we will introduce this integration technique here briefly.

Suppose $\mathbf{Y}={\mathbf{y}^{(1)},\dots,\mathbf{y}^{(n)}}$ and $d=2$, an integration space (a non-dominated space) of $\mathbf{Y}$ can be divided into $n+1$ disjoint integration slices ($S_{2}^{(i)},i=1,\dots,n+1$) by drawing lines parallel to $y_{2}$-axis at each element in $\mathbf{y}$, as indicated in Figure 3. Then, each integration slice can be expressed by its lower bound ($\mathbf{l}_{2}^{(i)}$) and upper bound ($\mathbf{u}_{2}^{(i)}$). In order to define the slices formally, we argue a Pareto-front approximation set $\mathcal{P}$ with two sentinels: $\mathbf{y}^{(0)}=(r_{1},\infty)^{\top}$ and $\mathbf{y}^{(n+1)}=(\infty,r_{2})^{\top}$. Then, the integration slices for the 2-D case are defined by:

[TABLE]

In the 2-D case, the number of integration slices is straightforward, namely, $N_{2}=n+1$.

4.1.2 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 ($\mathbf{l}_{3}$) and upper bound ($\mathbf{u}_{3}$). Since the upper bound of each integration slice is always $\infty$ in the $y_{3}$ axis, we can describe each integration slice as follows:

[TABLE]

Example 3

An illustration of integration slices is shown in Figure 4. A Pareto front set $\mathcal{P}$ is composed by 4 points ($\mathbf{y}^{(1)}=(1,3,4)^{\top},\mathbf{y}^{(2)}=(4,2,3)^{\top},\mathbf{y}^{(3)}=(2,4,2)^{\top}$ and $\mathbf{y}^{(4)}=(3,5,1)^{\top}$), and this Pareto front is shown in the upper left figure. The upper right figure represents the partitioned integration slices of $\mathcal{P}$. The lower center figure illustrates the projection of the upper right figure onto the $y_{1}y_{2}$-plane with rectangle slices and $\mathbf{l},\mathbf{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 $\mathbf{l}_{3}^{(4)}=(1,2,2)^{\top}$, and the upper bound of the slice is $\mathbf{u}_{3}^{(4)}=(2,4,\infty)^{\top}$.

Algorithm 2 describes how to obtain the slices $S_{3}^{(1)}$, $\dots$, $S_{3}^{(i)}$, $\dots$, $S_{3}^{(N_{3})}$ with the corresponding lower and upper bounds ($\mathbf{l}_{3}^{(i)}$ and $\mathbf{u}_{3}^{(i)}$). The partitioning algorithm is similar to the sweep line algorithm described in emmerich2011computing . 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 $\mathbf{y}^{(i)}$, the algorithm finds all the points $(\mathbf{y}^{(d[1])},\dots,\mathbf{y}^{(d[s])})$ which are dominated by $\mathbf{y}^{(i)}$ in the $y_{1}y_{2}$-plane and inserts $\mathbf{y}^{(i)}$ into the tree. Moreover, because of $\mathbf{y}^{(i)}$, the algorithm will also discard all the points ($\mathbf{y}^{\text{(d[1])}}$, $\dots$, $\mathbf{y}^{\text{(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 $\mathbf{y}^{\text{(t)}}$, $\mathbf{y}^{\text{(d[1])}}$, $\dots$, $\mathbf{y}^{\text{(d[s])}}$, $\mathbf{y}^{(r)}$, and $\mathbf{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,\dots,d$: the $i$-th coordinate is different for each pair of points in $\mathbf{Y}$). Otherwise, $2n+1$ provides an upper bound for the obtained number of slices.

Proof: In the algorithm, each point $\mathbf{y}^{(i)}|_{i=1,\dots,n}$ creates two slices. The first one, say slice $A^{(i)}$, is created when the point $\mathbf{y}^{(i)}$ is added to the AVL tree. Another slice, say slice $S_{3}^{(i)}$, is created when the point $\mathbf{y}^{(i)}$ is discarded from the AVL tree due to domination by another point, say $\mathbf{y}^{(s)}$, in the $y_{1}y_{2}$-plane. These two slices are defined as follows $A^{(i)}=((y^{(t)},y_{2}^{(l2)},y_{3}^{(i)})^{\top},(y_{1}^{(u1)},y_{2}^{(i)},\infty)^{\top})$ whereas $y_{2}^{(l2)}$ is either $y_{2}^{(r)}$ if no point is dominated by $\mathbf{y}^{(i)}$ in the $y_{1}y_{2}$-plane, or $y_{2}^{(d[1])}$, otherwise. Moreover, $S_{3}^{(i)}=((y^{(i)}_{1},y^{(r)}_{2},y^{(s)}_{3})^{\top},$ $(y^{(u)}_{1},y^{(s)}_{2},\infty)^{\top})$ and $\mathbf{y}^{(u)}$ denote either the right neighbor among the newly dominated points in the $y_{1}y_{2}$-plane, or $\mathbf{y}^{(s)}$ if $\mathbf{y}^{(i)}$ is the rightmost point among all newly dominated points. In this way, each slice can be attributed to exactly one point in $\mathcal{P}$, except the slice that is created in the final iteration. In the final iteration, one additional point $\mathbf{y}^{(n+1)}=(\infty,\infty,\infty)^{\top}$ 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.

4.1.3 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\geq 4$), the hyperboxes can be denoted by $S_{d}^{(1)},\dots,S_{d}^{(i)},\dots,S_{d}^{N_{d}}$ with their lower bounds ($\mathbf{l}^{(1)},\dots,\mathbf{l}^{(N_{d})}$) and upper bounds ($\mathbf{u}^{(1)},\dots,$ $\mathbf{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_{d}^{(i)}$ is defined as:

[TABLE]

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 dachert2017efficient by Dächert et al. and LKF17 lacour2017box by Lacour et al. Here, algorithm DKLV17 is an efficient algorithm to locate the local lower bound points666For the definition of the local lower bound points $\mathbf{L}_{d}$, see dachert2017efficient . ($\mathbf{L}_{d}$) in a dominated space for maximization problems, based on a specific neighborhood structure among local lower bounds. Moreover, LKF17 is an efficient algorithm to calculate the HVI by partitioning the dominated space. 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 $\mathcal{P}^{{}^{\prime}}$. This transformation is done by the following steps. Suppose that we have a current Pareto-front approximation set $\mathcal{P}$ for a maximization problem and we want to partition the non-dominated space of $\mathcal{P}$. Firstly, DKLV17 is applied to locate the local lower bound points ($\mathbf{L}_{d}$) of $\mathcal{P}$ in the dominated space. Secondly, regard $\mathbf{L}_{d}$ as a new Pareto-front approximation $\mathcal{P}^{{}^{\prime}}$ for a minimization problem with a reference point $\{\infty\}^{d}$. The dominated space of $\mathcal{P}^{{}^{\prime}}$ is actually the non-dominated space of $\mathcal{P}$. Then, LKF17 can be applied to partition the dominated space of $\mathcal{P}^{{}^{\prime}}$ by locating the lower bound points $\mathbf{l}_{d}$ and the upper bound points $\mathbf{u}_{d}$. These bound points ($\mathbf{l}_{d}$,$\mathbf{u}_{d}$) of $\mathcal{P}^{{}^{\prime}}$ in the dominated space for a minimization problem are exact the lower/upper bound points of the partitioned, non-dominated hyperboxes of $\mathcal{P}$ for a maximization problem. The pseudo code of partitioning non-dominated space in higher dimensional cases is shown in Algorithm 3.

Example 4

Figure 6 illustrates Algorithm 3. In the 2-D case, suppose a Pareto-front approximation set is $\mathcal{P}$, which consists of $\mathbf{y}^{(1)}=(1,2.5)^{\top}$, $\mathbf{y}^{(2)}=(2,1.5)^{\top}$ and $\mathbf{y}^{(3)}=(3,1)^{\top}$. The reference point is $\mathbf{r}=(0,0)$, see Figure 6 (above left). Use DKLV17 to locate the local lower bound points $\mathbf{L}_{2}$, which consist of $\mathbf{L}_{2}^{(1)}=(0,2.5)^{\top}$, $\mathbf{L}_{2}^{(2)}=(1,1.5)^{\top}$, $\mathbf{L}_{2}^{(3)}=(2,1)^{\top}$ and $\mathbf{L}_{2}^{(4)}=(3,0)^{\top}$, see Figure 6 (above right). Regard all of the local lower bound points $\mathbf{L}_{2}$ as the elements of a new Pareto-front approximation set $\mathcal{P}^{{}^{\prime}}=(\mathbf{L}_{2}^{(1)},\cdots,\mathbf{L}_{2}^{(4)})$. Set a new reference point $\mathbf{r}^{{}^{\prime}}=(\infty,\infty)$ and utilize LKF17 to partition the dominated space of $\mathcal{P}^{{}^{\prime}}$, by considering minimization, see Figure 6 (below left). The partitioned non-dominated space of $\mathcal{P}$ is then the partitioned dominated space of $\mathcal{P}^{{}^{\prime}}$, see Figure 6 (below right).

4.2 EHVI calculation

777Both 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.

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, $\Delta(\mathbf{y})$ is used whenever $\mathcal{P},\mathbf{r}$ are given by the context. Based on $\Delta(\mathbf{y})$, the expected hypervolume improvement function can be re-defined as:

[TABLE]

For the convenience of expressing the EHVI formula in the remaining parts of this paper, two functions ($\ell$ and $\vartheta$) are defined as follows:

Definition 11 ($\ell$ function)

Given the parameters of an integration slice $S_{d}^{(i)}$ in a $d$-dimensional space, the Hypervolume Improvement of slice $S_{d}^{(i)}$ in dimension $k\leq d$ is defined as:

[TABLE]

Definition 12 ($\vartheta$ function)

Given the parameters of an integration slice $S_{d}^{(i)}$ in a $d$-dimensional space and multivariate predictive distribution $\boldsymbol{\mu},\boldsymbol{\sigma}$, the function $\vartheta(l_{k}^{(i)},u_{k}^{(i)},\sigma_{k},\mu_{k})$ is then defined as:

[TABLE]

4.2.1 2-D EHVI calculation

According to the definition of the 2-D integration slice in Equation 4.1.1, the Hypervolume Improvement $\mathbf{y}\in\mathbb{R}^{2}$ in the 2-D case is:

[TABLE]

$\mbox{HVI}_{2}$ gives rise to the compact integral for the original EHVI:

[TABLE]

Here $\mathbf{y}=(y_{1},y_{2})$, the intersection of $S_{2}^{(i)}$ with $\Delta(y_{1},y_{2})$ is non-empty if and only if $(\mathbf{y})$ dominates the lower left corner of $S_{2}^{(i)}$. Therefore:

[TABLE]

In Equation (28), the summation is done after integration. This operation is allowed, because integration is a linear mapping. Moreover, the integration interval $\int_{y_{1}=l_{1}^{(i)}}^{\infty}$ can be divided into $(\int_{y_{1}=l_{1}^{(i)}}^{u_{1}^{(i)}}+\int_{y_{1}=u_{1}^{(i)}}^{\infty})$, because the HVI in one dimension $\lambda_{1}[S_{2}^{(i)}\cap\Delta(y_{1})]$ differs in these two integration intervals. Here $\lambda_{1}[B_{i}\cap\Delta(y_{k})]$ is the HVI in dimension $k$, i.e., a 1-D HVI. Equation (28) can then be expressed as:

[TABLE]

According to the definition of HVI, $\ell(u_{1}^{(i)},y_{1},l_{1}^{(i)})$ is constant and is $(u_{1}^{(i)}-l_{1}^{(i)})$ in $\int_{y_{1}=u_{1}^{(i)}}^{\infty}$. Therefore, the Expected Improvement in dimension $y_{1}$ is also a constant and it is: $\vartheta(l_{1}^{(i)},u_{1}^{(i)},\sigma_{1},\mu_{1})$. Recall the $\Psi_{\infty}$ function, by which the terms (29) and (30) can be expressed as follows:

[TABLE]

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 $\mathcal{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 $\mathcal{P}$ is sorted, we can show that the time complexity is in $\Theta(n)$.

4.2.2 3-D EHVI calculation

Given a partitioning of the non-dominated space into integration slices $S_{3}^{(1)}$, $\dots$, $S_{3}^{(i)}$, $\dots$, $S_{3}^{(2n+1)}$, the EHVI integrations over each slice can be computed separately. To see how this calculation can be done, the Hypervolume Improvement of a point $\mathbf{y}\in\mathbb{R}^{3}$ is rewritten as:

[TABLE]

where $\Delta\mathbf{(y)}$ is the part of the objective space that is dominated by $\mathbf{y}$. The HVI expression in the definition of EHVI in Equation (17) can be replaced by $\mbox{HVI}_{3}$ in Equation (33):

[TABLE]

Similar to the 2-D case, we can divide the integration interval $\int_{y_{1}=l_{1}^{(i)}}^{\infty}$ and $\int_{y_{2}=l_{2}^{(i)}}^{\infty}$ into $(\int_{y_{1}=l_{1}^{(i)}}^{u_{1}^{(i)}}+\int_{y_{1}=u_{1}^{(i)}}^{\infty})$ and $(\int_{y_{2}=l_{2}^{(i)}}^{u_{2}^{(i)}}+\int_{y_{2}=u_{2}^{(i)}}^{\infty})$, 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:

[TABLE]

Recalling the definition of the $\vartheta$ function and calculation of $\lambda_{1}[B_{i}\cap\Delta(y_{k})]$, the term (35) can be rewritten as follows:

[TABLE]

Similar to the derivation of the term (35), the terms (36), (37) and (38) can be written as follows:

[TABLE]

The final EHVI formula is the sum of the terms (39), (40), (41) and (42).

During the EHVI calculation, $y_{1}y_{2}$-projections are mutually non-dominated. Moreover, the points are sorted by the $y_{2}$ coordinate in the AVL tree. Therefore, identifying a neighboring/discard point takes time $O(\log n)$. Then the EHVI for these integration slices is calculated by the summation of the term (39), (40), (41) and (42), with the parameters of $\boldsymbol{\mu},\boldsymbol{\sigma}$ and $S_{3}^{(N_{3})}$. The EHVI time complexity for each slice is $O(1)$. Moreover, the dominated points ($\mathbf{y}^{\text{(d[s])}}$) are removed from the AVL tree, and the new points ($\mathbf{y}^{(j)}$) are inserted into the AVL tree. Since the points dominated by the new point $\mathbf{y}^{(j)}$ are deleted at the end of the current loop, they will not occur again in later computations. Hence, the total number of open slices does not exceed $N_{3}$, as mentioned before, and the total computational cost is $O(n\log n)$.

4.2.3 Higher dimensional EHVI

The interval of integration in each coordinate (except the last) can be divided into 2 parts: $[l,u]$ and $[u,\infty]$. Therefore, the equation for EHVI for each hyperbox can be decomposed into $2^{d-1}$ parts. For the interval of $[u,\infty]$, the improvements ($\lambda_{k}[S_{d}^{(i)}\cap\Delta(y_{k})]$) are constant values, and the $\Psi_{\infty}$ function can be simplified by calculating function $\Phi$ and the improvement in these coordinates. For the last coordinate, there is no need to separate the interval, because the improvement in this coordinate ($\lambda_{d}[S_{d}^{(i)}\cap\Delta(y_{d})]$) is a variable in $[l,\infty]$.

According to the definition of higher dimensional integral boxes in Section 4.1.3, the EHVI ($d\geqslant 4$) can be calculated by the following equation:

[TABLE]

[TABLE]

[TABLE]

In the term (43), the integral of each dimension $\int_{y_{k}=l_{k}^{(i)}}^{u_{k}^{(i)}}\lambda_{k}[S_{k}^{(i)}\cap\Delta(y_{1},\cdots,y_{k})]\cdot\boldsymbol{\xi}$ ${}_{\boldsymbol{\mu},\boldsymbol{\sigma}}(\mathbf{y})d\mathbf{y}|_{1\leq k\leq d-1}$ has two and only two different expressions ($\Psi_{\infty}$ or $\vartheta$), except that in the last dimension ($k=d$), the expression of $\int_{y_{d}=l_{d}^{(i)}}^{u_{d}^{(i)}}\lambda_{d}[S_{d}^{(i)}\cap\Delta(y_{1},\cdots,y_{d})]\cdot\boldsymbol{\xi}_{\boldsymbol{\mu},\boldsymbol{\sigma}}(\mathbf{y})d\mathbf{y}|_{k=d}$ is always $\Psi_{\infty}$. The final expression of the EHVI is the sum of the $EHVI_{S_{d}^{(i)}}$ in Equation (44) for all the partitioned, non-dominated hyperboxes $S_{d}^{(i)}|_{i=1,\cdots,N_{d}}$. In Equation (44), $EHVI_{S_{d}^{(i)}}$ has $2^{d-1}$ terms because the integration has two different expressions in dimension $k$ (${1\leq k\leq d-1}$) and has only one expression in $d$-th dimension ($k=d$).

In Equation (44), $(j)_{2}$ stands for the binary string representation of the integer $j$. The length of $(j)_{2}$ is $d-1$. $C^{(j)_{2}}_{k}$ is a bit and represents the $k$-th bit of $(j)_{2}$ in the binary string. For example, if $d=5$, $j=8$, then $(j)_{2}=(1\enskip 0\enskip 0\enskip 0)$, $C^{(j)_{2}}_{k=4}=1$ and $C^{(j)_{2}}_{k=1,2,3}=0$. In Equation (44), $\omega(i,k,C_{k}^{(j)_{2}})$ is defined as:

[TABLE]

Equation (44) shows how to calculate EHVI in the case of $d$ objectives. According to Equation (44), the runtime complexity of the proposed algorithm can be calculated. The exact EHVI is given by $\sum_{j=0}^{2^{d-1}-1}\left(\prod_{k=1}^{d-1}\omega(i,k,C_{k}^{(j)_{2}})\cdot\Psi_{\infty}(l_{d}^{(i)},l_{d}^{(i)},\mu_{d},\sigma_{d})\right)$, which requires $O(1)$ computation steps for each hyperbox calculation. Currently, the exact number of hyperboxes $N_{d}|_{d\geq 4}$ for a non-dominated space is still unknown. It is hypothesized by the authors that $N_{d}$ is the exact number of the local lower bound points, which can be calculated by the DKLV17 algorithm. The LKF17 algorithm partitions the non-dominated space into $\Theta(n^{\lfloor d/2\rfloor})$ hyperboxes and the time complexity for computing these hyperboxes grows linearly with the number of boxes (see, e.g., Lacour et al. lacour2017box ). Given a fixed dimension $d$, the time complexity of our EHVI computation algorithm is therefore also in $\Theta(n^{\lfloor d/2\rfloor})$.

Note that every hyperbox requires $O(2^{d-1})$ time complexity. Therefore, the complexity in terms of $d$ and $n$ is given by $O(2^{d-1}\cdot n^{\lfloor d/2\rfloor})$. Due to the exponential dependence on $d$, the EHVI computation algorithm is only useful in moderate dimensional cases. Note, however, that a much faster computation cannot be expected, as the time complexity of the hypervolume indicator itself scales superpolynomially with the number of objectives $d$, under the assumption P $\neq$ NP. It is easy to show, that the EHVI computation has at least the same time complexity than the hypervolume indicator computation Michael2016book , and it is therefore also an NP hard problem in $d$ – but polynomial in $n$ for any fixed value of $d$.

4.3 Probability of Improvement (PoI)

According to the partitioning method in Section 4.1, PoI can be calculated as follows:

[TABLE]

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 $\mathbf{r}=\{-\infty\}^{d}$ is only used in order to obtain the correct boundary information ($\mathbf{l}_{d},\mathbf{u}_{d}$).

5 Experiments

5.1 Speed Comparison

The test benchmarks from Emmerich and Fonseca emmerich2011computing were used to generate Pareto-front sets. The Pareto-front sets and evaluated points were randomly generated based on convexSpherical and concaveSpherical functions. Two EHVI calculation algorithms, CDD13 couckuyt2014fast and KMAC, were compared using the same benchmarks in this experiment.888Another 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 yang2017computing . Note that KMAC algorithm in this paper includes the KMAC_2D for 2-D EHVI calculation in Michael2016book , the KMAC_3D algorithm for 3-D EHVI calculation in yang2017computing , and the extended KMAC algorithm for higher dimensional cases ($d\geq 4$) in this paper.

The parameters: $\sigma_{d}=2.5$, $\mu_{d}=10,d=2,\cdots,5$ are used in the experiments. Pareto front sizes are $|\mathcal{P}|\in{\{10,20,\cdots,200\}}$ and Batch Size is 1, which represents the number of the evaluated points under the same Pareto-front approximation set. Ten $\mathcal{P}$ sets are randomly generated by the same parameters. Average runtimes over 100 repetitions (10 repetitions for 10 $\mathcal{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.

The experimental results in Figure 7 show that KMAC is much faster than CDD13, especially when $|\mathcal{P}|$ is increased. Moreover, CDD13 can not calculate the exact EHVI value within 30 minutes when $|\mathcal{P}|$ is bigger than 30 in the 5-D case.

5.2 Benchmark Performance

Five state-of-the-art algorithms are compared in this section, namely: EHVI-MOBGO, PoI-MOBGO999Using EHVI or PoI as the infill criterion in Algorithm 1., NSGA-II Deb2002NSGA2 , NSGA-III Deb2014NSGA3 ; Platypus and SMS-EMOA beume2007sms . The benchmarks are DTLZ1, DTLZ2, DTLZ3, DTLZ4, DTLZ5, DTLZ7 dtlz2002a , MaF1, MaF5, MaF12 and MaF13 Cheng2017 . 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 couckuyt2014fast . For the other test problems, the reference points are revised from the articles couckuyt2014fast ; Cheng2017 . 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 18-dimensional 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 and MaF13, even if the function evaluation budget of the EMOAs is increased to 2000, EHVI-MOBGO can still outperform the EMOAs in 6- and 18-dimensional search spaces. Among EHVI-MOBGO and PoI-MOBGO, EHVI-MOBGO outperforms PoI-MOBGO in most cases, but PoI-MOBGO yields better results on two and three (out of ten) test problems when $m=6$ and $m=18$, respectively. The reason is that the dominated space of the PoI is $\mathbf{y}\in(-\infty,\infty)^{d}\setminus\mbox{ dom }(\mathcal{P})$101010A reference point in Equation (16) is defined as $\mathbf{r}=(-\infty,\infty,\infty)$ in this paper. For PoI-MOBGO, a reference point in Table 3 is only used to evaluate the experimental results., which is bigger than that of the EHVI $\mathbf{y}\in(\mathbf{r},\infty^{d})\setminus\mbox{ dom }(\mathcal{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, $\mathbf{y}\in(-\infty^{d},\mathbf{r})\setminus\mbox{ dom }(\mathcal{P})$.

6 Conclusions and Outlook

This paper describes an efficient algorithm for EHVI calculation. It reviews and benchmarks recently proposed asymptotically optimal algorithms with $\Theta(n\log n)$ time complexity and generalizes them to higher dimensional cases with $d\geq 4$. By using the fast box decomposition techniques, which were recently developed by Dächert et al. dachert2017efficient and Lacour et al. lacour2017box , a non-dominated space can be partitioned with only $\Theta(n^{\lfloor d/2\rfloor})$ 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 $\Theta(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\leq 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 russo2014quick ; JASZKIEWICZ201872 .

Acknowledgements.

The authors gratefully acknowledge the contribution of Carlos M. Fonseca (University of Coimbra, Portugal) on the discussion of the integration techniques.