A Box Decomposition Algorithm
to Compute the Hypervolume Indicator
Renaud Lacour∗1 , Kathrin Klamroth1 , and Carlos M. Fonseca2

arXiv:1510.01963v2 [cs.DM] 8 Oct 2015

1

Department of Mathematics and Computer Science,
University of Wuppertal, Germany
{lacour,klamroth}@math.uni-wuppertal.de
2
CISUC, Department of Informatics Engineering, University of Coimbra, Portugal,
cmfonsec@dei.uc.pt
October 9, 2015

Abstract
We propose a new approach to the computation of the hypervolume indicator, based
on partitioning the dominated region into a set of axis-parallel hyperrectangles or boxes.
We present a nonincremental algorithm and an incremental algorithm, which allows inp−1
p
sertions of points, whose time complexities are O(nb 2 c+1 ) and O(nb 2 c+1 ), respectively.
While the theoretical complexity of such a method is lower bounded by the complexity
of the partition, which is, in the worst-case, larger than the best upper bound on the
complexity of the hypervolume computation, we show that it is practically efficient. In
particular, the nonincremental algorithm competes with the currently most practically
efficient algorithms. Finally, we prove an enhanced upper bound of O(np−1 ) and a lower
p
bound of Ω(nb 2 c log n) for p ≥ 4 on the worst-case complexity of the WFG algorithm.

1

Introduction

In multi-objective optimization (MOO), since objective functions are often conflicting in practice, there is typically no single solution that simultaneously optimizes all objectives. Instead
there are several efficient solutions, i.e. that cannot be improved on one objective without degrading at least another objective. Due to the possible large number of efficient solutions or
even nondominated points – their images in the objective space – approximation algorithms
are often favored in practice. These algorithms are able to generate several discrete approximations or representations of the nondominated set and quality indicators are required to
compare such approximations. Among them is the hypervolume indicator or Lebesgue measure (see e.g. Zitzler and Thiele, 1998), which measures the volume of the part of the objective
space dominated by the points of the approximation and bounded by some reference point.
Practically efficient algorithms to approximate the nondominated set include in particular
evolutionary multi-objective (EMO) algorithms. Most often the hypervolume indicator is used
∗

Corresponding author.

1

within EMO algorithms to compare the quality of the computed approximations (Beume et al.,
2007; Wagner et al., 2007). Because these algorithms repeatedly generate approximations of
large cardinality there is a clear need for efficient algorithms to compute the hypervolume
indicator. It has been shown (Bringmann and Friedrich, 2010) that, under “P 6= NP”, there
is no algorithm to compute the hypervolume indicator in time polynomial in the number of
objectives. Therefore, we will assume throughout the paper that the number of objectives,
although arbitrary, is fixed.
In this paper, we are interested in the efficient computation of the hypervolume associated
to a given set of points. We make a difference between the nonincremental case and the
incremental case. In the nonincremental case, the whole set of points has, in general, to be
known in advance because the computation approach imposes an order on the points to be
processed. In the incremental case, the hypervolume is updated each time a new point is
considered without restriction on the order new points come up.

1.1

Terminology and notations

We consider a minimization problem formulated as follows:
min f (x) = (f1 (x), . . . , fp (x))
s.t. x ∈ X

(1)

where X 6= ∅ is the feasible set and f1 , . . . , fp are p ≥ 2 objective functions mapping from X
to R. We assume that all feasible points f (x) : x ∈ X are located in some open hyperrectangle
of Rp , namely Z = (0, z r ). In order to compare points of the objective space, we define the
following binary relations. For z 1 , z 2 ∈ Rp :
z 1 5 z 2 (z 1 weakly dominates z 2 ) ⇔ zj1 ≤ zj2 ,
j = 1, . . . , p,
z1 ≤ z2
(z 1 dominates z 2 )
⇔ z 1 5 z 2 and z 1 6= z 2 ,
z 1 < z 2 (z 1 strictly dominates z 2 ) ⇔ zj1 < zj2 ,
j = 1, . . . , p.
For any set N of points of Rp , we define the dominated region as
D(N ) = {z 0 ∈ Rp : z 5 z 0 5 z r , for some z ∈ N }
where z r is used as a reference point. Considering the set
Nnd = {z ∈ N : z 0 6≤ z, for all z 0 ∈ N }
of all nondominated points of N , we have D(N ) = D(Nnd ). Therefore, we assume in the
remainder that N is a stable set of points for the dominance relation, i.e. for all z 1 , z 2 ∈ N ,
z 1 6≤ z 2 .
We denote by V (N ) the volume of the polytope D(N ) which is also referred to as the
hypervolume associated to N .
For any z ∈ Rp , we let z−j be the (p − 1)-dimensional vector of all components of z
excluding component j, for a given j ∈ {1, . . . , p}. Finally, for any z, a ∈ Rp and any
j ∈ {1, . . . , p}, (zj , a−j ) denotes the vector (a1 , . . . , aj−1 , zj , aj+1 , . . . , ap ).

2

1.2

Literature review on the computation of the hypervolume indicator

The currently most efficient algorithm for the computation of the hypervolume indicator
in the case p ≥ 4 in terms of theoretical worst-case complexity is by Chan (2013). For
the computation of the dominated hypervolume of n p-dimensional points, his algorithm
p
runs in O(n 3 polylog(n)) time. To our knowlege, there is currently, however, no available
implementation of this approach and no evidence of its practical efficiency. The complexity
of computing the hypervolume indicator is Θ(n log n) (see Beume et al., 2009).
On the practical point of view, the “HV4D” algorithm of Guerreiro et al. (2012), specialized for the case p = 4, is the most efficient in this case (see e.g. the computational results
of Russo and Francisco, 2014; Nowak et al., 2014). Their algorithm achieves O(n2 ) time
complexity. Above p = 4, the Walking Fish Group (WFG) (While et al., 2012) and Quick
Hypervolume (QHV) (Russo and Francisco, 2014) algorithms are currently two of the most
efficient according to the computational experiments conducted by the authors. The best
upper bounds on their complexity are, however, in both cases exponential.

1.3

Goals and outline

In this paper, we propose to compute the hypervolume indicator by partitioning the dominated
region into hyperrectangles. Given that this partition is determined by a set of corner points or
local upper bounds of the dominated region, we present approaches to compute these points.
Specifically, we propose a nonincremental approach which requires that the points for which
the hypervolume is computed are sorted with respect to one of the objective functions, and
an incremental approach which relaxes this assumption at some extra cost. Both approaches
have a good worst-case complexity and perform very well in practice. We also provide new
insight into the theoretical complexity of the WFG algorithm.
The remainder of this paper is organized as follows. In Section 2, we present the proposed
approach. Section 3 provides an enhanced analysis of the complexity of the WFG algorithm.
Section 4 describes computational experiments conducted to compare the proposed algorithms
to the state-of-the-art and presents the results. Section 5 concludes the paper.

2

The Hypervolume Box Decomposition Algorithm

This section describes the proposed new approach to compute the hypervolume indicator: the
Hypervolume Box Decomposition Algorithm (HBDA). Section 2.1 defines the decomposition
scheme. Section 2.2 presents two algorithms to compute an auxiliary set of points, referred
to as an upper bound set, that is necessary to obtain the decomposition. In Section 2.3, the
worst-case time complexity of HBDA is discussed. The general position assumption that is
made in the first two sections is relaxed in Section 2.4. Section 2.5 is concerned with the
implementation of HBDA.

2.1

Partitioning the dominated region into disjoint hyperrectangles

The dominated region D(N ) is a union of hyperrectangles of the type [z, z r ] where z ∈ N .
The idea of the proposed approach is to rely on another description of D(N ) as a union
of pairwise disjoint hyperrectangles so that the hypervolume can be computed as the sum of

3

z2

zr

ẑ 1

Points in N
All local upper bounds of U (N )
D(N )

z1
z2
z3
z4
z5

ẑ 2
z1

0

Figure 1: Decomposition of the dominated region in the bi-objective case
their volumes. We illustrate the concepts of this section on a bi-objective instance represented
in Figure 1.
In Kaplan et al. (2008, Section 3), a decomposition of D(N ) into a set of pairwise disjoint
hyperrectangles is established. Their approach is based on the computation of a set U (N ) of
auxiliary points associated to N , which are identical to local nadir points in the bi-objective
case, or to local upper bounds in the general multi-objective case (see Klamroth et al., 2015).
The set U (N ) is defined as the set of all the points u ∈ [0, z r ] that satisfy the two following
properties:
(P1 ) [0, u) does not contain any point of N .
(P2 ) u is maximal for property (P1 ), i.e. for any u0 ∈ [0, z r ] such that u ≤ u0 , u0 does not
satisfy (P1 ).
Following Klamroth et al. (2015), we denote by local upper bounds the elements of U (N ) and
by upper bound set the set U (N ) itself.
We first make a simplifying “general position” assumption that will be relaxed later.
Under this assumption, no two distinct points, among the points of N ⊂ Z to be considered,
share the same value in any dimension. We also define additional dummy points ẑ 1 , . . . , ẑ p
where ẑ j = (zjr , 0−j ) for all j ∈ {1, . . . , p}. These points are all the points of [0, z r ] that (1)
are not dominated by and do not dominate any point of Z and (2) are minimal with respect
to the dominance relation ≤ for (1).
Let N̂ = N ∪ {ẑ 1 , . . . , ẑ p }. We have U (N̂ ) = U (N ) and each vector u ∈ U (N ) is defined
j
by p distinct points {z 1 (u), . . . , z p (u)} of N̂ in the sense that zjj (u) = uj and z−j
(u) < u−j
for all j ∈ {1, . . . , p}. We refer to Klamroth et al. (2015) for more details on this aspect.
Then, according to Kaplan et al. (2008), the set {B(u) : u ∈ U (N )} where
B(u) =

[z11 (u), z1r ]

×

p 
Y
j=2

4

max{zjk (u)}, uj
k<j


(2)

is a partition of D(N ). We refer to Figure 1 for an illustration of this partition in the
bi-objective case. Therefore, the hypervolume of D(N ) is obtained as the sum of the hypervolumes of the B(u)’s, for all u ∈ U (N ), which requires a computation time linear in
|U (N )|.
In the next section, we discuss the computation of the set U (N ).

2.2

Computing upper bound sets

Several approaches exist for the computation of the set U (N ). Kaplan et al. (2008) propose
a nonincremental algorithm which requires that the points of N are sorted in increasing
order according to any fixed component. Their algorithm is output-sensitive and achieves
O(|U (N )| logp−1 |N |) time complexity using dynamic range trees. Przybylski et al. (2010) also
provide an algorithm to compute U (N ) which is used to define search zones for an algorithm
to compute the nondominated set of MOCO problems. In Klamroth et al. (2015) we propose
another algorithm which uses, as in Dächert and Klamroth (2015) for the tri-objective case,
the relation between points of U (N ) and their defining points to avoid a filtering step with
respect to Pareto dominance found in Przybylski et al. (2010). Both Przybylski et al. (2010)
and Klamroth et al. (2015) propose incremental algorithms. We note that since defining points
are tracked in Kaplan et al. (2008) and Klamroth et al. (2015), the corresponding algorithms
make it directly possible to compute the hypervolume indicator using the partition described
by (2).
We first present the algorithm of Klamroth et al. (2015), which is used when arbitrary
insertions into N are required (Section 2.2.1). Then we present the nonincremental algorithm
based on both Kaplan et al. (2008) and a theorem of Klamroth et al. (2015), which is expected
to be more efficient when the whole set N is known in advance (Section 2.2.2).
2.2.1

Incremental algorithm

The incremental algorithm to compute an upper bound set is presented in Algorithm 1. Given
a stable set N and a new point z̄ such that N ∪ {z̄} is also a stable set, it identifies from the
upper bound set for N the set A of all local upper bounds that no longer satisfy property (P1 )
with respect to z̄ (Step 4). From the set A, Step 2 generates the valid new local upper bounds
using the result provided in Theorem 2.1.
Theorem 2.1 (Klamroth et al., 2015). Let z̄ be a point of (0, z r ) such that N ∪ {z̄} is a stable
set of points in general position. Consider a local upper bound u ∈ U (N ) such that z̄ < u.
Then, for any j ∈ {1, . . . , p}, (z̄j , u−j ) is a local upper bound of U (N ∪ {z̄}) if, and only
if, z̄j ≥ maxk6=j {zjk (u)}.
To be able to use this result, it is required to keep track of the associated defining points
for each local upper bound. This is done by setting z k (z r ) ← ẑ k at Step 1 and
(
z̄
if k = j
z k (z̄j , u−j ) ←
k
z (u) otherwise
at Step 2, for all k ∈ {1, . . . , p}.

5

Algorithm 1: Incremental algorithm to compute an upper bound set
1
2
3
4
5

ubsI(∅)

= {z r }

ubsI(N ∪ {z̄}) = {(z̄j , u−j ) : z̄j ≥ maxk6=j {zjk (u)}, u ∈ A, j = 1, . . . , p}
∪A
where A = {u ∈ ubsI(N ) : z̄ < u}
A = ubsI(N ) \ A

2.2.2

Nonincremental algorithm

In the case where points of N are given in increasing order of some component, say component
p, then the computation of U (N ) can take advantage of this. The approach is described in
Algorithm 2 and corresponds to the algorithm of (Kaplan et al., 2008, Section 3.1) with the
addition of the condition from Theorem 2.1 at Step 3. All local upper bounds u of A satisfy
up = zpr , therefore all (z̄p , u−p ) with u ∈ A are valid new local upper bounds. The other new
local upper bounds are obtained as in Algorithm 1 considering only the first p − 1 components
of z̄. In addition to the update of defining points described in Section 2.2.1 that also has to
be performed for the nonincremental algorithm, we have to set z p (z̄p , u−p ) ← z̄ at Step 2.
Algorithm 2: Nonincremental algorithm to compute an upper bound set – assumes
that z̄p > zp , for all z ∈ N
= {z r }

1

ubsNI(∅)

2

ubsNI(N ∪ {z̄}) = {(z̄p , u−p ) : u ∈ A}
∪{(z̄j , u−j ) : z̄j ≥ maxk6=j {zjk (u)}, u ∈ A, j = 1, . . . , p − 1}
∪A
where A = {u ∈ ubsNI(N ) : z̄ < u}
A = ubsNI(N ) \ A

3
4
5
6

For the computation of the hypervolume indicator, the local upper bounds of A need not
be kept. Indeed they are not modified later, which implies that the associated hyperrectangles
according to (2) will not change.
Algorithms 1 and 2 and the decomposition of the dominated region (Section 2.1) yield
two Hypervolume Box Decomposition Algorithms: an incremental version (HBDA-I) and a
nonincremental version (HBDA-NI), respectively.

2.3

Time complexity of the algorithms

The time complexity of HBDA-I and HBDA-NI mainly depends on the size of the current
upper bound set U (N ) and of the set A in Algorithms 1 and 2, respectively. Indeed for both
algorithms, a constant time is spent on each element of A. Let tI (n, p) and tNI (n, p) be upper
bounds on the time complexity of HBDA-I and HBDA-NI, respectively, applied to n points
of dimension p. We obtain the following relations:
tI (n, p) ≤ tI (n − 1, p) +
s(n − 1, p)
tNI (n, p) ≤ tNI (n − 1, p) + s(n − 1, p − 1)
6

where s(n, p) is the worst-case size of an upper bound set on n points of dimension p, which
p
is equal to Θ(nb 2 c ) (Kaplan et al., 2008). (We recall that, in the case of HBDA-NI, only an
upper bound set for at most n (p − 1)-dimensional points needs to be maintained.) Therefore,
we have:
p

tI (n, p) = O(nb 2 c+1 )
tNI (n, p) = O(nb

p−1
c+1
2

).

Note that here we assume the worst-case for |U (N )| but Algorithm 2 is an output-sensitive
algorithm (Kaplan et al., 2008).

2.4

Relaxing the simplifying “general position” assumption

Real or generated instances may contain points that are not in general position. Therefore, it is
important to allow, in algorithms to compute the hypervolume indicator, points having equal
component values in the same dimension. The partition of the dominated region presented in
Section 2.1 is based on the existence and uniqueness, for each local upper bound u, of a p-uple
of points that define the p components of u. The uniqueness in particular is guaranteed by the
general position assumption and Theorem 2.1 assumes general position. We show, however,
that Algorithms 1 and 2 can be applied without any modification to non-general position
instances.
Proposition 2.2. Algorithms 1 and 2 are still valid when the input points are in non-general
position.
Proof. Comparison between component values of points appear at three different places in
the algorithms, namely (a) in the strict dominance tests at Steps 4 and 5 of Algorithms 1
and 2, respectively, (b) in the application of the condition of Theorem 2.1 at Steps 2 and 3
of Algorithms 1 and 2, respectively, and (c) in the computation of the hypervolume of a box
according to (2).
For (a), note that the strict dominance tests at Steps 4 and 5 of Algorithms 1 and 2,
respectively, should remain the same, since they are related to property (P1 ), which does not
assume general position.
For (b) and (c), we follow the idea of symbolically perturbating the component values of the
input points as suggested in Kaplan et al. (2008). More precisely, given a set N = {z 1 , . . . , z n }
of points in non-general position, one can define for each component j ∈ {1, . . . , p} a total
order <j on the values of the points of N on component j:
zji1 <j zji2 if, and only if, zji1 < zji2 or (zji1 = zji2 and i1 > i2 )
for all i1 , i2 ∈ {1, . . . , n}, i1 6= i2 . This relation is obviously compatible with the natural strict
ordering on real numbers, in the sense that if zji1 < zji2 then zji1 <j zji2 . Thus we can use <j ,
or more precisely the symmetric >j in place of ≥ to apply the condition of Theorem 2.1 in
the non-general position case. In fact, if we label the points of N so that the current point
always gets the largest index among the points considered so far, the relation ≥ can be used
equivalently to >j , since the left-hand side of the comparison is always a component value of
the current point.
Finally, it is equivalent to use <j or ≤ to compute the maxima in (2).
7

2.5

Implementation and data structures

Computing the set A in ubsI and ubsNI requires to determine the subset of a set of local
upper bounds that are strictly dominated by a given point. We expect that the size of the
output subset is much smaller than the cardinality of the input set, therefore a specialized
data structure could be used instead of a simple linked list.
Options for this are range trees, kd-trees and generalized quadtrees (de Berg et al., 2008).
Generalized quadtrees become inefficient for large dimensional point sets, since the number
of children of an internal node is equal to 2p . Besides, the chosen data structure needs to
handle both insertions and deletions, which is costly to achieve with range trees and kd-trees.
Therefore, for the incremental algorithm ubsI, we still suggest to store local upper bounds
in a linked list. The list is sorted in nondecreasing order of the sum of the component values
of each local upper bound to avoid some of the dominance tests.
For the case of the nonincremental algorithm ubsNI, we propose to take advantage of
the information provided by the point set N , which is assumed to be known in advance.
We suggest to use a combination of a kd-tree and a set of linked lists. Namely, we build a
balanced kd-tree from the points of N . Then we consider the partition of the objective space
induced by this kd-tree. For each new local upper bound, we identify the cell of the partition
it belongs to by traversing the tree. Instead of creating a new node for this local upper bound,
we insert it in a linked list located in place of this potential new node. To perform Step 5
of ubsNI, local upper bounds strictly dominated by a given point z̄ are identified by first
searching the tree with the query interval (z̄, z r ) and then the linked lists containing local
upper bounds in some cell intersecting (z̄, z r ). The corresponding local upper bounds can
then be removed in constant time without altering the tree structure.

3

A better bound on the worst-case time complexity of the
WFG algorithm

In this section, we propose an improved analysis of the worst-case time complexity of the
WFG algorithm of While et al. (2012) While the upper bound provided by the authors is
O(2n ), we show that it can be lowered at least to O(np−1 ), matching the upper bound of the
Hypervolume by Slicing Objectives (HSO) of While et al. (2006). Moreover, we show that its
p
complexity is at least Ω(nb 2 c log n) for p ≥ 4.
In Section 3.1, we briefly describe the WFG algorithm. Then in Section 3.2 we prove the
new upper and lower bounds.

3.1

Brief description of the WFG algorithm

The WFG algorithm computes the hypervolume in a recursive way. Given a stable set N ∪{z̄}
of points and a reference point z r , the volume of D(N ∪ {z̄}) is computed as the sum of the
volume of D(N ) and the exclusive hypervolume associated to z̄, i.e. V (N ∪ {z̄}) − V (N ).
This last quantity is equal to V ({z̄}) − V (N 0 ), where N 0 is the stable subset of all orthogonal
projections of the points of N onto the dominance cone {z ∈ Rp : z̄ 5 z}. This idea is
illustrated with a 3-dimensional example represented in Figure 2, where N = {z 1 , . . . , z 7 } and
it is assumed that z3 ≤ z̄3 , for all z ∈ N . The basic algorithm is summarized in Algorithm 3,
where pmax(z, z 0 ) is the parallel maximum of z and z 0 , i.e. pmax(z, z 0 ) = (max{zj , zj0 })j=1,...,p .
8

z2

zr
Points in N ∪ {z̄}
Points in {pmax(z̄, z) : z ∈ N }
D(N )
D(N 0 )

z1
z2

z6
z3
z7
z4
z̄
z5
z1

0

Figure 2: Orthogonal projection on the hyperplane z3 = z̄3 of a 3-dimensional instance

Algorithm 3: The basic WFG algorithm
wfg(∅)
wfg({z̄})

=0
Q
= pj=1 zjr − z̄j

wfg(N ∪ {z̄}) = wfg(N ) + wfg({z̄}) − wfg(N 0 )
where N 0 = {pmax(z̄, z) : z ∈ N }nd

9

If the points are sorted in non-decreasing order of some component, say component p, then,
as remarked by While et al. (2006), wfg(N 0 ) is a (p − 1)-dimensional subproblem. Indeed,
the dominated region of N 0 , D(N 0 ), can be described by sliding its orthogonal projection on
the hyperplane of equation fp = z̄p along the fp axis, from level z̄p to level zpr . Therefore,
the hypervolume of D(N 0 ) can be obtained by multiplying by (z r − z̄p ) the hypervolume of
D(N 0 ) projected on the hyperplane defined by fp = z̄p . Besides, a simple algorithm is used
for subproblems with p = 2. The basic algorithm together with this important enhancement
is what we refer to as WFG algorithm hereafter.

3.2

New upper and lower bounds

We first show an improved upper bound in Proposition 3.1.
Proposition 3.1. The worst-case time complexity of the WFG algorithm in the case where
the points are sorted in non-decreasing order of some arbitrary component is bounded by
O(np−1 ).
Proof. Let t(n, p) be the worst-case time complexity of the WFG algorithm as described in
Algorithm 3, with n = |N | and n0 = |N 0 |. We have:
t(n + 1, p) ≤ t(n, p) + c1 + c2 n2 + t(n0 , p − 1)
where t(n, p), c1 , and t(n0 , p − 1) are the worst-case time complexities of wfg(N ), wfg({z̄}),
and wfg(N 0 ), respectively, and c2 n2 is an upper bound on the complexity of computing N 0 ,
c1 and c2 being constants with respect to n and n0 . We have t(n, 2) = O(n) if the points are
sorted in non-decreasing order of some component, therefore we obtain t(n, 3) = O(n3 ). It
follows, since n0 ≤ n that
t(n, p) = O(np )
for any integers n and p. If an O(n log n)-algorithm (Beume et al., 2009) or even just an
O(n2 )-algorithm is used for the case p = 3, then the complexity becomes
t(n, p) = O(np−1 )
for any integers n and p ≥ 3.
Now we prove a non-trivial lower bound in Proposition 3.2
p

Proposition 3.2. The worst-case complexity of the WFG algorithm is Ω(nb 2 c log n) for p ≥ 4.

|
k k − 1 ··· 1
Proof. Let Ak =
and 0k be a null matrix of the same order as Ak . For
1
2
··· k

any k ≥ 1 and any even p ≥ 4 we define the following p2 k × p -dimensional block-matrix:


Ak 0k · · · 0k

. 
 0k . . . . . . .. 

(3)
Mk,p = 

 . .
.. ... 0 
 ..
k
0k · · · 0k Ak
and consider the p-dimensional instance where the points are the n = p2 k rows of Mk,p .
Note that the rows are already sorted in non-decreasing order of component p. Consider any
10

point z̄ taken among the last k rows of Mk,p , denoted by z̄ 1 , . . . , z̄ k , and the computation
of wfg(N ∪ {z̄}), where N is the set of all rows above z̄ in Mk,p . The set N 0 involved in
this computation contains all the first ( p2 − 1)k rows of Mk,p with the last two components
replaced by z̄p−1 and z̄p , respectively, and, in the case z̄ 6= z̄ 1 , the vector (0, . . . , 0, z̄p−1 +1, z̄p )
(and possibly further points, depending on the value of j ∈ {1, . . . , k} satisfying z̄ = z̄ j , that
0 are sorted in
are however dominated by (0, . . . , 0, z̄p−1 + 1, z̄p ) in N 0 ). If the points of Nnd
non-decreasing order of component p − 1, all (case z̄ = z̄ 1 ) or all but the last point (case
0 have identical values on components p − 1 and p, the other components being
z̄ 6= z̄ 1 ) of Nnd
those of Mk,p−2 . Therefore, for each of these k computations, there will be a call of wfg on
an instance of the same structure with ( p2 − 1)k points of dimension p − 2. Thus we have:
 p


t(n, p) ≥ k · t
− 1 k, p − 2
2
which, given that n = p2 k, is equivalent to
2
t(n, p) ≥ n · t
p


p
2


n
,p − 2 .
−1
p

Given that the recursions with p = 2 are solved in Θ(n log n), we obtain t(n, p) = Ω(n 2 log n)
when p is even.
For an odd p ≥ 5, we apply the above analysis for p − 1, adding a zero-column to Mk,p−1 .
p−1
We then obtain t(n, p) = Ω(n 2 log n) in this case.
p
Overall we have shown that t(n, p) = Ω(nb 2 c log n) for any p ≥ 4.
Note that we assumed in the proof of Proposition 3.2 that, for each recursive call of wfg,
the sorted component is imposed to the algorithm. This corresponds to the description given
by the authors as well as their implementation. Other choices for the sorted component than
the one given in the proof can make the solution to the instance type given in the proof a lot
more efficient.

4

Computational experiments

In this section, we provide the setup (Sections 4.1 and 4.2) and the results (Section 4.3) of
computational experiments conducted to compare the efficiencies of Algorithms 1 and 2, as
well as the HV4D (Guerreiro et al., 2012), QHV (Russo and Francisco, 2014), and WFG
(While et al., 2012) algorithms to compute the hypervolume indicator.

4.1

Implementations and general setup

We implemented HBDA-I and HBDA-NI in C. We used for the HV4D, QHV, and WFG
algorithms the implementations provided by the authors, namely Fonseca et al. (2010, version
2.0 RC 2), Russo and Francisco (2013, retrieved on April, 2015), and While et al. (2014, version
1.10), respectively. We note that, in the implementation of the QHV algorithm, the number
of objectives is fixed at compile time, thus the final executable may take advantage of this
information, which the other implementations do not.
From the implementation of the WFG algorithm, we derived an incremental version. In
this version, the initial set of points is not assumed to be sorted, thus it may not be known
11

entirely in advance. The auxiliary sets N 0 (see Algorithm 3) are, however, sorted since they
are built from points for which the exclusive hypervolume has already been computed. In
other words, this corresponds to considering that the initial set of points is sorted according
to an extra (p + 1)-th component (compatible with the order in which the points are actually
processed) that is, however, not taken into account in the computation of the hypervolume.
We denote by ”WFG incremental“ this approach.
All implementations were compiled using GCC 4.3.4 with the same options -O3 -DNDEBUG
-march=native. Compilation and tests were both performed under SUSE Linux Enterprise
Server 11 on identical workstations equipped with four Intel Xeon E7540 CPU at 2.00GHz
and with 128GB of RAM.

4.2

Instances

All algorithms were tested on stable sets of points generated according to several schemes.
We considered the following four instance types:
(C) Concave or so-called spherical (Deb et al., 2002) instances
(X) Convex instances
(L) Linear instances
(H) Hard instances
Instances of types (C), (X), and (L) are obtained by drawing uniformly points from the open
hypercube (0, 1)p . Then each point z is modified as follows:
z
(C) zj ← √Ppj

2
k=1 zk

their `2 -norm,

z
(X) zj ← 1 − √Ppj

, for each j ∈ {1, . . . , p}, i.e. the component values of z are divided by

2
k=1 zk

, for each j ∈ {1, . . . , p},

z

(L) zj ← Pp j z , for each j ∈ {1, . . . , p}, i.e. the component values of z are divided by their
k=1 k
`1 -norm.
Note, in particular, that for type (C) and (X) instances, we followed the suggestion of Russo
and Francisco (2014) to project uniformly distributed points on a hypersphere, instead of
using the instances of Deb et al. (2002).
The points of the instances of types (C), (X), and (L) are randomly distributed on some
subset of P
(0, 1)p , namely S(1, 1) ∩ (0, 1)p for type (C), S(0, 0) ∩ (0, 1)p for type (X), and
{z ∈ Rp : pj=1 zj = 1} ∩ (0, 1)p for type (L), where S(z, r) denotes a hypersphere of Rp with
center z ∈ Rp and radius r ∈ R+ .
Instances of type (H) are motivated by the instance type built to derive a lower bound
on the worst-case complexity of the WFG algorithm in Section 3.2 (Equation 3). We slightly
modified these instances for the computational experiments since the presence of many identical component values could perturb the comparison depending on the way algorithms handle this case. Therefore, we defined the following instance type in general position which
yields the same lower bound on the worst-case complexity of the WFG algorithm. Let

12

Algorithm
Optimization direction Coordinates range Ref. point # points
HBDA-I, HBDA-NI
any
any
any
any
HV4D
minimize
any
any
any
QHV
maximize
[0, 1]
1
≤ 1 000
WFG
maximize
any
any
any
Table 1: Restrictions of the tested implementations of the hypervolume algorithms considered
for computational experiments

|
k + lk k − 1 + lk · · · 1 + lk
Ak,l =
. For any even p and k, a hard instance is de1 + lk
2 + lk
· · · k + lk
fined by the set of all rows of the following block matrix:

Ak, p −1 Ak, p −2 · · ·
Ak,1
Ak,0
2
2

..
.. 
.
 Ak,0
Ak, p −1 . .
.
. 


2
 .
.. 
.

.
.
= .
. Ak, p −2
. 
Ak,0

2


..
 ..

..
. Ak, p −1 Ak, p −2 
 .
.
2
2
Ak, p −2 Ak, p −3 · · ·
Ak,0
Ak, p −1


0
Mk,p

2

2

2

and consists of p2 k points of dimension p. Matrices Ak,l have the same property of matrices
Ak , which is to have decreasing coefficients on the first column and increasing in the second
0 is also built by repeating along the diagonal the same matrix, which
column. Moreover, Mk,p
has strictly larger coefficients than the other matrices of the block matrix.
Some implementations of the algorithms we consider in this section have certain restrictions on the instances that can be solved. We summarize these restrictions in Table 1.
Because of the restrictions of the implementation of the QHV algorithm, we normalized
the hard instances so that the points are all in the open hypercube (0, 1)p . Moreover, for all
instances, we chose as reference point the all-ones vector and the null vector in the minimization and maximization cases, respectively. To cope with the restrictions on the optimization
direction, we generated instances primarily for the minimization case, and for each instance, a
symmetric instance with points in (0, 1)p for the maximization case, by taking for each point
z and component j the complement 1 − zj .
For types (C), (X), and (L) we generated instances for each p ∈ {4, . . . , 10} and n ∈
{100, 200, . . . , 1 000}. For type (H) we considered p ∈ {4, 6, 8, 10} and for each value of p
we chose 10 values for the parameter k so as to obtain 10 sizes of instances. Since for any
algorithm, type (H) instances are significantly harder to solve than the other types considered
in this paper, we limited the maximal number of points to 1 000, 900, 300, and 150 for
p = 4, 6, 8, 10, respectively.
For a fixed type, number of objectives, and number of points, we generated 10 instances.
The implementations were run up to 100 times on small instances to obtain significant computation times. The results we report are thus averaged.

13

4.3

Results

Figures 3, 4, 5, and 6 show computation times for nonincremental algorithms on instances of
type (C), (X), (L) and (H), respectively.
Our approach HBDA-NI performs better than all other algorithms on types (C), (X), and
(L), for p ∈ {5, 6, 7}, and on type (H) for all the tested dimensions above 4. For 4-dimensional
instances of any type, it is confirmed that the HV4D algorithm performs the best while the
computation times obtained with HBDA-NI are very close. HBDA-NI is still faster than the
QHV algorithm for p ∈ {8, 9, 10} on concave instances.
We also show computation times for the incremental algorithms (HBDA-I and the incremental implementation of the WFG algorithm) in Figures 7, 8, 9, and 10.
According to these results, the incremental WFG algorithm performs significantly better
than Algorithm 1 for almost all instances except on 6 and 8 objectives hard instances, where
both algorithms behave similarly. The relative poorer efficiency of our incremental approach
can be explained by the fact that its implementation lacks an efficient data structure to identify
local upper bounds strictly dominated by the point that is currently processed. Also the
incremental WFG algorithm is still able to reduce the dimension of the points in subproblems.

5

Conclusion

In this paper, we investigated a new way of computing the hypervolume indicator by calculating a partition of the dominated region into hyperrectangles. This decomposition is based on
the computation of local upper bounds. We proposed an incremental and a nonincremental
approach. These approaches provide a good worst-case complexity and the nonincremental version, through an efficient implementation is very competitive in practice. In fact, we
demonstrate that computing explicitly the dominated region, in the sense of computing all its
vertices, i.e. all local upper bounds additionally to feasible points, is an interesting approach
with respect to the computation time.
Future work includes improving the incremental approach. This could be done by implementing an efficient dynamic data structure to identify the local upper bounds that have to
be updated or removed when a new point is considered, which is subject to current work.
Alternatively, a special property of the dominated region such as the neighborhood relation
between local upper bounds elaborated in Dächert et al. (2015) could be exploited.
It would also be interesting to refine the analysis of the worst-case time complexity of
the WFG algorithm, because the gap between the lower and upper bounds we showed is still
large. Finally, one could think of using the concept of local upper bounds and the associated
decomposition of the dominated region to compute hypervolume contributions.

References
N. Beume, B. Naujoks, and M. Emmerich. SMS-EMOA: Multiobjective selection based on
dominated hypervolume. European Journal of Operational Research, 181(3):1653–1669,
2007. doi: 10.1016/j.ejor.2006.08.008 .
N. Beume, C. M. Fonseca, M. López-Ibáñez, L. Paquete, and J. Vahrenhold. On the complexity of computing the hypervolume indicator. Trans. Evol. Comp, 13(5):1075–1082, 2009.
doi: 10.1109/TEVC.2009.2015575 .
14

p=4

p=5
0.1

0.01
0.01
0.001
0.001
100

200

100

1 000

200

Computation time (seconds)

p=6

1 000
p=7

1

10

0.1

1
0.1

0.01
100

200

100

1 000

200

p=8

1 000
p=9

100

1 000

10

100
10

1

1
0.1

100

200

1 000

100

200

1 000

p = 10
10 000

100

1

100

200

1 000

Number of points
HBDA-NI

WFG

HV4D

QHV

Figure 3: Nonincremental algorithms on type (C) instances

15

p=4

p=5
0.1

0.01
0.01
0.001
0.001
0.0001

100

200

1 000

100

200

p=6

p=7

1

0.1

Computation time (seconds)

1 000

0.1

0.01

0.01
100

200

100

1 000

200

1 000
p=9

p=8
100

10

10
1
1
0.1
0.1
100

200

100

1 000

200

1 000

p = 10

100

1
100

200

1 000

Number of points
HBDA-NI

WFG

HV4D

QHV

Figure 4: Nonincremental algorithms on type (X) instances

16

p=5

p=4
0.1
0.01
0.01
0.001
0.001
100

200

100

1 000

200

Computation time (seconds)

p=6

1 000
p=7

1

10

0.1

1
0.1

0.01
100

200

1 000

100

200

1 000
p=9

p=8

100

10

10

1

1
0.1
100

200

1 000

100

200

1 000

p = 10

100

1
100

200

1 000

Number of points
HBDA-NI

WFG

HV4D

QHV

Figure 5: Nonincremental algorithms on type (L) instances

17

p=4

p=6

Computation time (seconds)

1

100

1

0.01

0.01
100

200

1 000

90

180

p=8

900

p = 10
1 000

100
1

1

0.01

0.001
28

60

300

15

30

150

Number of points
HBDA-NI

WFG

HV4D

QHV

Figure 6: Nonincremental algorithms on type (H) instances
p=4

p=5

0.1

Computation time (seconds)

1
0.01

0.1
0.01

0.001
100

200

1 000

100

200

1 000
p=7

p=6
100

10

10

1

1

0.1

0.1
100

200

1 000

100

200

1 000

Number of points
HBDA-I

WFG incremental

Figure 7: Incremental algorithms on type (C) instances

18

p=4

p=5

Computation time (seconds)

0.1

1
0.1

0.01

0.01
0.001
100

200

1 000

100

200

p=6

1 000
p=7

100
1
1
0.01
100

200

1 000

100

200

1 000

Number of points
HBDA-I

WFG incremental

Figure 8: Incremental algorithms on type (X) instances
p=4

p=5

Computation time (seconds)

0.1

1
0.1

0.01

0.01
0.001
100

200

1 000

100

200

p=6

1 000
p=7

100

10
1

1
0.1
100

200

1 000

100

200

1 000

Number of points
HBDA-I

WFG incremental

Figure 9: Incremental algorithms on type (L) instances

19

p=6

p=4
1 000

Computation time (seconds)

1
10
0.01

0.1
100

200

1 000

90

180

900

p=8

100
1
0.01
28

60

300

Number of points
HBDA-I

WFG incremental

Figure 10: Incremental algorithms on type (H) instances
K. Bringmann and T. Friedrich. Approximating the volume of unions and intersections of
high-dimensional geometric objects. Computational Geometry, 43(6–7):601–610, 2010.
doi: 10.1016/j.comgeo.2010.03.004 .
T. M. Chan. Klee’s Measure Problem Made Easy. In Foundations of Computer Science
(FOCS), 2013 IEEE 54th Annual Symposium on, pages 410–419, 2013.
doi: 10.1109/
FOCS.2013.51 .
K. Dächert and K. Klamroth. A linear bound on the number of scalarizations needed to solve
discrete tricriteria optimization problems. Journal of Global Optimization, 61(4):643–676,
2015. doi: 10.1007/s10898-014-0205-z .
K. Dächert, K. Klamroth, R. Lacour, and D. Vanderpooten. Efficient computation of the
search region in multi-objective optimization. Technical report, University of Wuppertal,
2015.
M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars. Computational Geometry: Algorithms and Applications. Springer, Santa Clara, CA, USA, 3rd edition, 2008.
K. Deb, L. Thiele, M. Laumanns, and E. Zitzler. Scalable multi-objective optimization test
problems. In Evolutionary Computation, 2002. CEC ’02. Proceedings of the 2002 Congress
on, volume 1, pages 825–830, 2002. doi: 10.1109/CEC.2002.1007032 .
C. M. Fonseca, L. Paquete, M. Lopez-Ibanez, and A. P. Guerreiro. Computation of the
Hypervolume Indicator. http://iridia.ulb.ac.be/~manuel/hypervolume, 2010.
20

A. P. Guerreiro, C. M. Fonseca, and M. T. Emmerich. A fast dimension-sweep algorithm
for the hypervolume indicator in four dimensions. In Proceedings of the 24th Canadian
Conference on Computational Geometry (CCCG 2012), pages 77–82, 2012.
H. Kaplan, N. Rubin, M. Sharir, and E. Verbin. Efficient Colored Orthogonal Range Counting.
SIAM Journal on Computing, 38(3):982–1011, 2008. doi: 10.1137/070684483 .
K. Klamroth, R. Lacour, and D. Vanderpooten. On the representation of the search region in
multi-objective optimization. European Journal of Operational Research, 245(3):767–778,
2015. doi: 10.1016/j.ejor.2015.03.031 .
K. Nowak, M. Märtens, and D. Izzo. Empirical Performance of the Approximation of the
Least Hypervolume Contributor. In T. Bartz-Beielstein, J. Branke, B. Filipič, and J. Smith,
editors, Parallel Problem Solving from Nature – PPSN XIII, volume 8672, chapter Lecture
Notes in Computer Science, pages 662–671. Springer International Publishing, 2014. doi:
10.1007/978-3-319-10762-2 65 .
A. Przybylski, X. Gandibleux, and M. Ehrgott. A two phase method for multi-objective
integer programming and its application to the assignment problem with three objectives.
Discrete Optimization, 7(3):149–165, 2010. doi: 10.1016/j.disopt.2010.03.005 .
L. M. S. Russo and A. P. Francisco. Quick Hypervolume. http://web.tecnico.ulisboa.
pt/luis.russo/QHV/#down, 2013.
L. M. S. Russo and A. P. Francisco. Quick Hypervolume. Evolutionary Computation, IEEE
Transactions on, 18(4):481–502, 2014. doi: 10.1109/TEVC.2013.2281525 .
T. Wagner, N. Beume, and B. Naujoks. Pareto-, Aggregation-, and Indicator-Based Methods
in Many-Objective Optimization. In S. Obayashi, K. Deb, C. Poloni, T. Hiroyasu, and
T. Murata, editors, Evolutionary Multi-Criterion Optimization, volume 4403, chapter Lecture Notes in Computer Science, pages 742–756. Springer Berlin Heidelberg, 2007. doi:
10.1007/978-3-540-70928-2 56 .
L. While, P. Hingston, L. Barone, and S. Huband. A faster algorithm for calculating hypervolume. Evolutionary Computation, IEEE Transactions on, 10(1):29–38, 2006. doi:
10.1109/TEVC.2005.851275 .
L. While, L. Bradstreet, and L. Barone. A Fast Way of Calculating Exact Hypervolumes.
Evolutionary Computation, IEEE Transactions on, 16(1):86–95, 2012.
doi: 10.1109/
TEVC.2010.2077298 .
L. While, L. Bradstreet, and L. Barone. Walking Fish Group: hypervolume project. http:
//www.wfg.csse.uwa.edu.au/hypervolume/, 2014.
E. Zitzler and L. Thiele. Multiobjective optimization using evolutionary algorithms — A
comparative case study. In A. Eiben, T. Bäck, M. Schoenauer, and H.-P. Schwefel, editors,
Parallel Problem Solving from Nature — PPSN V, volume 1498, chapter Lecture Notes
in Computer Science, pages 292–301. Springer Berlin Heidelberg, 1998.
doi: 10.1007/
BFb0056872 .

21

