Density Descent for Diversity Optimization

David H. Lee Anishalakshmi V. Palaparthi Matthew C. Fontaine

University of Southern California Los Angeles, CA<EMAIL_ADDRESS>University of
Southern California Los Angeles, CA<EMAIL_ADDRESS>University of Southern
California Los Angeles, CA<EMAIL_ADDRESS>
Bryon Tjanaka Stefanos Nikolaidis

University of Southern California Los Angeles, CA<EMAIL_ADDRESS>University
of Southern California Los Angeles, CA<EMAIL_ADDRESS>
###### Abstract

Diversity optimization seeks to discover a set of solutions that elicit
diverse features. Prior work has proposed Novelty Search (NS), which, given a
current set of solutions, seeks to expand the set by finding points in areas
of low density in the feature space. However, to estimate density, NS relies
on a heuristic that considers the $k$-nearest neighbors of the search point in
the feature space, which yields a weaker stability guarantee. We propose
Density Descent Search (DDS), an algorithm that explores the feature space via
gradient descent on a continuous density estimate of the feature space that
also provides stronger stability guarantee. We experiment with DDS and two
density estimation methods: kernel density estimation (KDE) and continuous
normalizing flow (CNF). On several standard diversity optimization benchmarks,
DDS outperforms NS, the recently proposed MAP-Annealing algorithm, and other
state-of-the-art baselines. Additionally, we prove that DDS with KDE provides
stronger stability guarantees than NS, making it more suitable for adaptive
optimizers. Furthermore, we prove that NS is a special case of DDS that
descends a KDE of the feature space.

##  1 Introduction

We study how stable, continuous approximations of density can accelerate the
search for diverse solutions to optimization problems.

There is a range of applications where searching directly for behavioral
diversity results in finding suboptimal solutions that act as “stepping
stones,” mitigating convergence to local optima (Gaier et al.,, 2019). A
classic example is the problem of training an agent to reach a target position
in a _deceptive maze_ (Lehman and Stanley,, 2011). There, directly minimizing
the distance between the agent’s final position and a target goal causes the
agent to get stuck. On the other hand, the problem can be solved by ignoring
the objective of directly reaching the target and instead attempting to find a
diverse range of agents, each of which reaches a different region of the maze.

We refer to the problem of finding a range of solutions that are diverse with
respect to prespecified features as diversity optimization (DO). We
characterize DO as a special instance of quality diversity optimization (QD)
(Pugh et al.,, 2016). Like DO, QD seeks a diverse set of solutions, but QD
also considers an objective function that the solutions must optimize. For
instance, in deceptive maze, one could add an energy consumption objective so
that the goal is to find a set of agents that minimize energy consumption
while moving to diverse regions in the maze. Thus, DO is an instance of QD
where the objective value of all solutions is constant.

Figure 1: We propose density descent search (DDS) for solving diversity
optimization (DO) problems. DDS first draws solutions from a Gaussian
$\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$. After computing the solution features
(in this case, the final position of the robot in a maze), DDS ranks solutions
by density. This density ranking is passed to CMA-ES, which updates the search
distribution to sample solutions with lower density on the next iteration.
Concurrently, solutions are stored in a buffer that forms the basis for
density estimates, and in a passive archive that tracks all discovered
solutions.

In DO and QD, discovering solutions with new features is challenging because
the mapping from solutions to features is complex and non-invertible (Chenu et
al.,, 2021). For example, in maze exploration, it is not known _a priori_ how
to produce an agent that travels to a given $(x,y)$ position.

A classic method applied to DO is Novelty Search (NS) (Lehman and Stanley,,
2011). NS retains an _archive_ of previously found individual solutions. It
aims to expand the archive by finding solutions that are far in the feature
space from existing solutions in the archive. Specifically, NS optimizes for
solutions with a high _novelty score_ , where novelty score is the average
distance in feature space from a solution to its $k$-nearest neighbors in the
archive.

While the novelty score is framed in terms of distance, we study its
interpretation as an approximation of _density_ (Chenu et al.,, 2021), where
density is directly proportional to the number of solutions in a region of
feature space. A high novelty score indicates that a solution’s features are
far from those of its $k$-nearest neighbors in the archive, i.e., it is
located in an area of the feature space with low density.

An alternative approach to DO is to apply general-purpose QD algorithms and
assume that all solutions in the search space have identical quality, i.e.,
equal to some constant. Covariance Matrix Adaptation MAP-Annealing (CMA-MAE)
(Fontaine and Nikolaidis,, 2023) is a state-of-the-art QD algorithm that
optimizes for archive improvement with the CMA-ES (Hansen,, 2016) optimizer.
CMA-MAE retains a discrete archive of solutions in feature space. When applied
to solve a DO problem by ignoring the solution quality, this archive naturally
becomes a histogram that represents the distribution of solutions in feature
space, with lower values of the histogram indicating lower density. CMA-MAE
then performs _density descent_ on this histogram, where it continually seeks
to fill the areas of low density.

We emphasize that, to efficiently explore feature space, both NS and CMA-MAE
leverage density estimates of solutions in feature space — NS is guided by
novelty score, while CMA-MAE reads from its histogram. However, both of these
density estimates have drawbacks. On one hand, the novelty score in NS is
_continuous_ but provides a _weaker_ stability guarantee, meaning that its
value can change arbitrarily when new features are discovered. On the other
hand, the histogram in CMA-MAE provides a _stronger_ stability guarantee but
utilizes a _discrete_ approximation that gives flat gradient signals on its
bins.

_Our key insight is to overcome the drawbacks of current density estimation
methods in DO by introducing continuous, stable approximations of the solution
density in feature space._ This insight results in the following
contributions: (1) We propose two new algorithms that query continuous density
estimates to guide an optimizer, in our case CMA-ES, to search for solutions
that are diverse in the feature space (Sec. 4). The first algorithm, DDS-KDE,
leverages non-parametric Kernel Density Estimation (KDE) (Chen,, 2017). The
second algorithm, DDS-CNF, leverages continuous normalizing flow (Lipman et
al.,, 2023) to learn the underlying density function. (2) We show that, when
combined with a ranking-based optimizer like CMA-ES, NS reduces to a special
case of DDS-KDE (Sec. 5). (3) We prove that KDE provides stronger algorithmic
stability guarantees than novelty score (Sec. 5). (4) We demonstrate that our
DDS algorithms outperform prior work on 3 out of 4 domains (Sec. 6). (5) We
show that our algorithms perform well on multi-dimensional feature spaces,
which currently present significant challenges to QD and DO algorithms.

##  2 Problem Definition

Quality Diversity Optimization (QD). For solution
${\bm{\theta}}\in\mathbb{R}^{n}$, QD assumes an _objective_ function
$f:\mathbb{R}^{n}\to\mathbb{R}$ and $m$ _feature_ functions, jointly
represented as a vector-valued function
${\bm{\phi}}:\mathbb{R}^{n}\to\mathbb{R}^{m}$. We refer to the image of
${\bm{\phi}}$ as the _feature space_ $S$. The goal of QD is to find, for each
${\bm{s}}\in S$, a solution ${\bm{\theta}}$ where $f({\bm{\theta}})$ is
maximized and ${\bm{\phi}}({\bm{\theta}})={\bm{s}}$.

Diversity Optimization (DO). DO is a special case of QD where the objective is
constant, i.e., $f({\bm{\theta}})=C$. The goal of DO is to find, for each
${\bm{s}}\in S$, a solution ${\bm{\theta}}$ where
${\bm{\phi}}({\bm{\theta}})={\bm{s}}$.

##  3 Background

###  3.1 QD and DO Algorithms

Various methods address the QD and DO problems by relaxing them to the problem
of finding an _archive_ (i.e., a set) $\mathcal{A}$ of representative
solutions. The structure of the archive defines two major families of
algorithms: those based on Novelty Search and those based on MAP-Elites.

Novelty Search (NS). The key insight of NS (Lehman and Stanley,, 2011) is to
discover diverse solutions in feature space by optimizing for solutions that
are “novel” with respect to a current set of solutions. Given a set of
features $\mathcal{X}\subseteq\mathbb{R}^{m}$, the novelty of a solution
${\bm{\theta}}\in\mathbb{R}^{n}$ and its features
${\bm{y}}={\bm{\phi}}({\bm{\theta}})$ relative to $\mathcal{X}$ is
encapsulated in its _novelty score_ , denoted $\rho({\bm{y}};\mathcal{X})$:

$\rho({\bm{y}};\mathcal{X})=\frac{1}{k}\sum_{{\bm{y}}^{\prime}\in
N_{k}({\bm{y}};\mathcal{X})}\mathrm{dist}({\bm{y}},{\bm{y}}^{\prime})$ (1)

where $N_{k}({\bm{y}};\mathcal{X})$ is the set of $k$-nearest neighbors to
${\bm{y}}$ in $\mathcal{X}$, and $\mathrm{dist}$ is a distance function —
henceforth, we consider $\mathrm{dist}$ to be some norm $\lVert\cdot\rVert$.

NS maintains an archive $\mathcal{A}$ of unbounded size and generates
solutions with an underlying optimizer, traditionally a genetic algorithm. For
each solution ${\bm{\theta}}$ produced by the optimizer, NS computes the
novelty score with respect to the current archive, i.e.,
$\rho({\bm{y}};\mathcal{A})$. If $\rho({\bm{y}};\mathcal{A})$ exceeds an
_acceptance threshold_ , then $({\bm{\theta}},{\bm{y}})$ is added to the
archive. In this manner, NS gradually adds novel solutions to the archive and
explores the feature space.

Importantly, the novelty score is non-stationary in that the novelty of a
solution ${\bm{\theta}}$ changes as the archive $\mathcal{A}$ is updated. We
prove in Sec. 5 that the degree of non-stationarity in novelty score is
unbounded for general feature spaces, resulting in significant optimization
challenges for adaptive optimizers such as CMA-ES.

Multi-dimensional Archive of Phenotypic Elites (MAP-Elites). While NS was
developed for DO, MAP-Elites (Mouret and Clune,, 2015) was developed for the
general QD setting. Compared to the unstructured archive of NS, MAP-Elites
divides the feature space into a predefined tessellation
$T:\mathbb{R}^{m}\to\\{1,\dots,l\\}$, where $e\in\\{1,\dots,l\\}$ is a _cell_
in the tessellation and $l$ is the total number of cells. Given a cell $e$,
the MAP-Elites archive associates a solution ${\bm{\theta}}$ with cell $e$ if
and only if the solution’s features are contained in $e$, i.e.,
$T({\bm{\phi}}({\bm{\theta}}))=e$. Moreover, MAP-Elites stores at most one
solution for every cell in the tessellation.

During a QD search, MAP-Elites first draws solutions from a predefined
distribution, e.g., a Gaussian, and inserts the solutions into the archive.
Subsequently, MAP-Elites generates and inserts new solutions by mutating
existing archive solutions with a genetic operator, such as the Iso+LineDD
operator (Vassiliades and Mouret,, 2018). Importantly, when solutions inserted
into the archive land in the same tessellation cell, MAP-Elites only retains
the solution with the greatest objective value. Thus, MAP-Elites gradually
collects high-performing solutions that have diverse features. We apply MAP-
Elites to DO by setting a constant objective $f(\cdot)=C$ over the solution
space.

The choice of tessellation in MAP-Elites can significantly impact scalability.
The most common tessellation is a grid tessellation (Mouret and Clune,, 2015),
which divides the feature space into equally-sized hyperrectangles. In high-
dimensional feature spaces, grid tessellations require exponentially more
memory due to the curse of dimensionality. Thus, a common alternative is the
centroidal Voronoi tessellation (CVT) (Vassiliades et al.,, 2018), which
divides the feature space into $l$ evenly-sized polytopes.

Covariance Matrix Adaptation Evolution Strategy (CMA-ES). One recent line of
QD algorithms has combined CMA-ES (Hansen,, 2016) with MAP-Elites. CMA-ES is a
state-of-the-art single-objective optimizer that represents a population of
solutions with a multivariate Gaussian
$\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$. Each iteration, CMA-ES draws
$\lambda$ solutions from the Gaussian. Based on the rankings (rather than the
raw objectives) of the solutions, CMA-ES adapts the covariance matrix
${\bm{\Sigma}}$ to regions of higher-performing solutions. While CMA-ES is a
derivative-free optimizer, it has been shown to approximate a natural gradient
(Akimoto et al.,, 2010).

Covariance Matrix Adaptation MAP-Elites (CMA-ME). The first work to integrate
CMA-ES with MAP-Elites was CMA-ME (Fontaine et al.,, 2020). The key idea of
CMA-ME is to optimize for archive improvement with CMA-ES. Namely, in addition
to a MAP-Elites-style archive, CMA-ME maintains an instance of CMA-ES. The
CMA-ES instance samples solutions from a Gaussian, and the solutions are
ranked based on how much they improve the archive, e.g., solutions that found
new cells in the archive are ranked high, while those that were not added at
all are ranked low. With this _improvement ranking_ , CMA-ES adapts the
Gaussian to sample solutions that will further improve the archive. Note that
in DO settings, since the objective is a constant value, the improvement
ranking places solutions that filled a new cell ahead of solutions that did
not add a new cell.

Covariance Matrix Adaptation MAP-Annealing (CMA-MAE). One limitation of CMA-ME
is that it focuses too much on exploring for new cells in feature space rather
than optimizing the objective (Tjanaka et al.,, 2022). CMA-MAE (Fontaine and
Nikolaidis,, 2023) addresses this problem by introducing an _archive learning
rate_ $\alpha$. When $\alpha=1$, CMA-MAE maintains the same exploration
behavior as CMA-ME, but when $\alpha=0$, CMA-MAE focuses solely on optimizing
the objective, like CMA-ES. For $0<\alpha<1$, CMA-MAE blends between these two
extremes, enabling it to both explore feature space and optimize the
objective.

In DO settings, where the objective is constant ($f(\cdot)=C$), CMA-MAE
naturally performs _density descent_ (Fontaine and Nikolaidis,, 2023).
Intuitively, the archive becomes a histogram that represents how many
solutions have been found in each region of feature space. Lower values of the
histogram indicate lower _density_ of solutions, i.e., if a histogram bin has
a low value, few solutions have been discovered in that region of feature
space. CMA-MAE then seeks to _descend_ the histogram by searching for
solutions that fill the low-valued histogram bins.

###  3.2 Density Estimation

We consider two density estimation methods: Kernel Density Estimation and
Continuous Normalizing Flows.

Kernel Density Estimation (KDE). KDE is a non-parametric density estimation
method that does not make any assumptions on the underlying probability
density distribution from which samples are drawn (Chen,, 2017). Given a set
of features $\mathcal{X}\subseteq\mathbb{R}^{m}$, bandwidth parameter $h$, and
kernel function $K(\cdot)$, KDE computes the density function
$\hat{D}:\mathbb{R}^{m}\to\mathbb{R}$ for a given feature
${\bm{y}}\in\mathbb{R}^{m}$ as:

$\displaystyle\hat{D}_{h}({\bm{y}};\mathcal{X})=\frac{1}{|\mathcal{X}|h}\sum_{{\bm{y}}^{\prime}\in\mathcal{X}}K\left(\frac{\lVert{\bm{y}}-{\bm{y}}^{\prime}\rVert}{h}\right)$
(2)

Prior work has thoroughly studied the problem of selecting bandwidth that
accurately estimate the underlying density function (Silverman,, 1986;
Rudemo,, 1982; Bowman,, 1984). While accurate density estimation is important
for KDE, we show in Theorem 5.2 that larger bandwidth results in more stable
density estimates (albeit at the cost of accuracy), which is beneficial to the
proposed algorithms in Sec. 4.

When optimizing a density estimate (e.g., with gradient descent), KDEs have
several advantages over histograms. First, the shape of a histogram depends on
its bin size (Weglarczyk,, 2018). Second, the binning procedure leads to a
discontinuous optimization landscape and flat gradient signals inside each bin
(Weglarczyk,, 2018). In contrast, KDEs produce a smooth, continuous density
function based on the location of the samples in the underlying distribution
(Chen,, 2017).

Continuous Normalizing Flow (CNF). CNF (Lipman et al.,, 2023) is a generative
modeling method that constructs a diffusion path between a simple distribution
(e.g., a Gaussian distribution) and an unknown distribution. The diffusion
path describes a mapping from a point on the simple distribution to a
corresponding point on the unknown distribution such that their probability
densities are roughly equal. CNFs enable sampling from probability
distributions where direct sampling is difficult (e.g., distributions of
images). This is accomplished by sampling from the simple distribution and
transforming the sample via the learned diffusion path of the CNF. However, we
utilize CNF to estimate the density of feature space for our proposed
algorithms in Sec. 4.

##  4 Density Descent Search

We present Density Descent Search (DDS), a DO algorithm that efficiently
explores the feature space by leveraging continuous density estimates. DDS
maintains a buffer of features ${\mathcal{B}}$ that represents the
distribution of features discovered so far. Based on this buffer, DDS models a
density estimation of the feature distribution, and by querying this density
estimate, DDS guides an adaptive optimizer to discover solutions in less-dense
areas of feature space. As it searches for such solutions, DDS expands the
discovered set of features. Below, we describe the components of DDS.

Density Estimation. The core of DDS is its representation of feature space
density. We propose two variants of DDS that differ in their density
representation. The first variant, DDS-KDE, represents density with KDE. With
the KDE, we compute density via Eq. 2, with $\mathcal{X}$ set to the buffer of
features ${\mathcal{B}}$.

The second variant, DDS-CNF, estimates the density in feature space by
learning a CNF between the standard normal distribution
$\mathcal{N}(\bm{0},{\bm{I}})$ and the observed feature space distribution.
Similar to DDS-KDE, the feature space distribution is represented by the
buffer ${\mathcal{B}}$. We compute a density estimate at any location in the
feature space by integrating an ordinary differential equation (Lipman et
al.,, 2023). Applying techniques from prior work (Lipman et al.,, 2023), we
represent the CNF with a neural network, and on each iteration, we finetune
the network on the new distribution of features contained in the buffer.

KDE and CNF differ in how easily their smoothness can be controlled. On one
hand, the shape of the KDE can be controlled with the bandwidth
hyperparameter. Higher bandwidth leads to a smoother density estimate but
conceals modes of the true density function, as illustrated in 2(a). On the
other hand, while CNF does not require selecting a bandwidth hyperparameter,
the smoothness of the density estimation cannot be easily controlled, which
undermines the performance of the algorithm (Sec. 6.3). We discuss the impacts
that the smoothness of the density estimate has on DDS in Sec. 5.

Feature Buffer. To provide the basis for density estimates in feature space,
DDS maintains a buffer ${\mathcal{B}}$ that stores the features of sampled
solutions (Algorithm 1). This buffer represents the _observed_ distribution of
the feature space and is updated every time new solutions and features are
discovered. In theory, the buffer can have infinite capacity, storing every
feature ever encountered. In practice, the buffer can only retain a finite
number of features due to computation and memory limitations. To decide which
features to retain in the buffer, we manage the buffer with an optimal
reservoir sampling algorithm (Li,, 1994). This algorithm updates the buffer
with online samples that accurately represent the distribution of features
discovered so far by DDS (Algorithm 1).

Optimizer. DDS guides an adaptive optimizer to discover solutions in low-
density regions of the feature space. Examples of such optimizers include xNES
(Glasmachers et al.,, 2010) and Adam (Kingma and Ba,, 2015). We select CMA-ES
as the underlying optimizer due to its reputation as a state-of-the-art
optimizer (Hansen,, 2016) and its high performance in existing QD algorithms
(Fontaine et al.,, 2020; Fontaine and Nikolaidis,, 2021, 2023).

Passive Archive. Following prior work (Pugh et al.,, 2015), to track how much
of the feature space has been explored, DDS inserts all discovered solutions
and features into a passive MAP-Elites-style (Sec. 3) archive $\mathcal{A}$.
While DDS itself only uses the archive to store solutions and features, the
archive is useful when computing metrics in experiments (Sec. 6.1).

Summary. Algorithm 1 shows how the components of DDS come together. DDS begins
by initializing its various components (Algorithm 1-1). During the main loop
(Algorithm 1), DDS samples solutions with the underlying optimizer (Algorithm
1). Since our optimizer is CMA-ES, sampling consists of drawing from a
multivariate Gaussian $\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$. Subsequently,
DDS computes the features of the solutions (Algorithm 1) and estimates their
feature space density (Algorithm 1). After sampling solutions, DDS updates its
various components. For instance, in DDS-KDE, the density update (Algorithm 1)
consists of replacing the feature buffer, and in DDS-CNF, the update involves
fine-tuning the neural network on the feature buffer. Furthermore, solutions
with lower density are ranked first on Algorithm 1, causing CMA-ES to adapt
(Algorithm 1) towards less-dense regions of the feature space.

1 DDS  _( $\bm{\phi}(\cdot)$, $density$, $b,N,\lambda,{\bm{\mu}}_{0},\sigma$)_

input : Feature function $\bm{\phi}(\cdot)$, density estimator object
$density$, buffer size $b$, number of iterations $N$, batch size $\lambda$,
initial solution ${\bm{\mu}}_{0}$, and initial step size $\sigma$

result : Generates $N\cdot\lambda$ solutions, storing a representative subset
of them in a passive archive $\mathcal{A}$.

2

3 Initialize CMA-ES search point ${\bm{\mu}}:={\bm{\mu}}_{0}$, search
direction ${\bm{\Sigma}}=\sigma{\bm{I}}$, and internal parameters ${\bm{p}}$

4 Initialize empty archive $\mathcal{A}$

5 Initialize empty buffer ${\mathcal{B}}$ of size $b$

6 for  _$iter\leftarrow 1$ to $N$_ do

7

8 for  _$i\leftarrow 1$ to $\lambda$_ do

9 ${\bm{\theta}}_{i}\sim\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$

10 ${\bm{\phi}}_{i}\leftarrow{\bm{\phi}}({\bm{\theta}}_{i})$

11 $D_{i}\leftarrow density({\bm{\phi}}_{i})$

12 Add $({\bm{\theta}}_{i},{\bm{\phi}}_{i})$ to archive $\mathcal{A}$

13 end for

14 Update buffer ${\mathcal{B}}$ with ${\bm{\phi}}_{1..\lambda}$ via reservoir
sampling

15 Update $density$ with the new buffer ${\mathcal{B}}$

16 Rank ${\bm{\theta}}_{i}$ in _ascending_ order by $D_{i}$

17 Adapt CMA-ES parameters ${\bm{\mu}},{\bm{\Sigma}},{\bm{p}}$ based on
density ranking $D_{i}$

18 end for

19

20

Algorithm 1 Density Descent Search (DDS)

##  5 Connection between Novelty Search and Kernel Density Estimates

We provide theoretical insight into the connection between NS and DDS-KDE and
delineate the advantage of KDE over novelty score. KDE and novelty score are
non-stationary since they change as more solutions are discovered by their
respective optimizers. However, the magnitude of non-stationarity is
potentially different in KDE and novelty score (Theorem 5.2, 5.3).
Furthermore, when $k\geq|{\mathcal{B}}|$ in the $k$-nearest neighbor
calculation for novelty score, novelty search becomes a special case of DDS-
KDE (Theorem 5.4). Finally, we conjecture that any meaningful density
estimator utilizing a point’s distance to its $k$-nearest neighbors will incur
similar weak stability guarantees as novelty score (5.5). We include all
proofs in Appendix E.

Stability Under Non-stationarity. DO algorithms such as NS, DDS, and CMA-MAE
gradually learn a non-stationary density representation (e.g., KDE or
histogram) as they explore the feature space. However, drastic changes in the
representation present significant challenges for adaptive optimizers, such as
Adam (Kingma and Ba,, 2015) and CMA-ES. To characterize the extent of change
in the density estimate, we appeal to the notion of _uniform stability_ (Hardt
et al.,, 2016), defined as follows:

######  Definition 5.1.

Let a function $F({\bm{x}};{\mathcal{B}}):\mathbb{R}^{d}\to\mathbb{R}$ be
parameterized by a set ${\mathcal{B}}\in\mathbb{R}^{d}$. We say that $F$ is
$\epsilon$-uniformly stable if for all
${\mathcal{B}},{\mathcal{B}}^{\prime}\subseteq\mathbb{R}^{d}$, where
${\mathcal{B}}$ and ${\mathcal{B}}^{\prime}$ differ by at most one element, we
have

$\sup_{{\bm{x}}\in\mathbb{R}^{d}}\left|F({\bm{x}};{\mathcal{B}})-F({\bm{x}};{\mathcal{B}}^{\prime})\right|\leq\epsilon$
(3)

We note that uniform stability has close connections to influence functions in
statistics, which measures how much an estimator deviates from the ground
truth when given a subset of the data (Cook and Weisberg,, 1980). In this
interpretation, uniform stability is an upper-bound of the (empirical)
influence function.

We prove that KDE has uniform stability directly proportional to the size of
its feature buffer and bandwidth:

######  Theorem 5.2.

A kernel density estimate $\hat{D}_{h}({\bm{x}};{\mathcal{B}})$ managed with
reservoir sampling, such that features in the buffer are exchanged one at a
time, is $\frac{1}{|{\mathcal{B}}|h}$-uniformly stable, where $h$ is the
bandwidth.

Higher bandwidth makes the function more stationary and thus more suited for
adaptive optimizers. However, higher bandwidth also leads to _over-smoothing_
of the density estimate, which conceals modes of the true density function
(Chen,, 2017) (2(a)). Thus, selecting an optimal bandwidth for DDS requires a
fine balance between the accuracy and uniform stability of the KDE.

In contrast, novelty score becomes less uniformly stable as the diameter of
the feature space $W$ increases:

######  Theorem 5.3.

Novelty score $\rho({\bm{x}};{\mathcal{B}})$ is $\frac{W}{k}$-uniformly
stable, where $k$ is the nearest neighbors parameter in novelty score and
$W=\max_{{\bm{s}}_{1},{\bm{s}}_{2}\in S}\lVert{\bm{s}}_{1}-{\bm{s}}_{2}\rVert$
is the diameter of the feature space $S$.

Therefore, for unbounded feature spaces, the uniform stability of novelty
score is also unbounded, and for bounded feature spaces, novelty score has
uniform stability directly proportional to $k$.

When the feature space is bounded, as in our experiments, KDE has a stronger
stability guarantee than novelty score for bandwidth
$h\geq\frac{k}{|{\mathcal{B}}|}$. Following this insight, we select a
bandwidth satisfying this inequality with the bandwidth selection mechanism in
Appendix B. Our theoretical results are corroborated by our experiments in
Sec. 6, which demonstrate that DDS-KDE outperforms NS in all domains on all
metrics.

Equivalence of NS and DDS-KDE. We observe that, when $h=1$ and as
$k\to|{\mathcal{B}}|$, the uniform stability upper-bound of novelty score
approaches that of KDE. We prove in Theorem 5.4 that when all points are
considered in the novelty score computation (i.e., let $k=\infty$), NS is a
special case of DDS-KDE under ranking-based optimizers such as CMA-ES.
Specifically, we demonstrate that under these conditions, the ranking of
solutions based on their novelty score is identical to the ranking based on
their kernel density estimate.

######  Theorem 5.4.

Let ${\mathcal{B}}\subseteq\mathbb{R}^{m}$ be a set of features. Consider the
rankings $\pi_{\mathrm{NS}}$ and $\pi_{\mathrm{KDE}}$ on another set of
features $\\{{\bm{\phi}}_{1},\dots,{\bm{\phi}}_{m}\\}\subseteq\mathbb{R}^{m}$
where ${\bm{\phi}}_{i}={\bm{\phi}}({\bm{\theta}}_{i})$. We define the rankings
as follows:
$\pi_{\mathrm{NS}}(i)\geq\pi_{\mathrm{NS}}(j)\iff\rho({\bm{\phi}}_{i};{\mathcal{B}})\geq\rho({\bm{\phi}}_{j};\mathcal{X})$
and
$\pi_{\mathrm{KDE}}(i)\geq\pi_{\mathrm{KDE}}(j)\iff\hat{D}({\bm{\phi}}_{i};{\mathcal{B}})\leq\hat{D}({\bm{\phi}}_{j};{\mathcal{B}})$.
We show that $\pi_{NS}=\pi_{KDE}$, when NS has $k=\infty$ (or, equivalently,
$k=|{\mathcal{B}}|$) and KDE has triangular kernel $K(u)=1-|u|$ with support
over the entire feature space $S$.

Stability of ${\bm{k}}$-NN. We conjecture that only considering $k$-nearest
neighbors in the density representation, for $k<|{\mathcal{B}}|$, will
inevitably result in poor uniform stability bounds. Notice that the maximum
distance between any two features is the diameter of the feature space, $W$.
Thus, when a new feature replaces a former $k$-nearest neighbors of
${\bm{x}}$, the total distance between a feature ${\bm{x}}$ to its $k$-nearest
neighbors can change by $W$ in the worst case.

######  Conjecture 5.5.

Let $\mathcal{X}\subseteq\mathbb{R}^{m}$ and $D:\mathbb{R}^{m}\to\mathbb{R}$
be a density estimator. If $D({\bm{x}};\mathcal{X})$ for
${\bm{x}}\in\mathbb{R}^{m}$ is computed based on the distance to its
$k$-nearest neighbors of ${\bm{x}}$, then $D({\bm{x}};\mathcal{X})$ is
$O\left(W\right)$-uniformly stable for $k<|{\mathcal{B}}|$.

A proof of this conjecture would imply that density estimators based on
$k$-nearest neighbor metric have unbounded uniform stability in unbounded
feature spaces, for $k<|{\mathcal{B}}|$. Thus, the potential for rapid changes
in the optimization landscape makes such estimators incompatible with adaptive
optimizers.

##  6 Experiments

We compare the performance of DDS-KDE and DDS-CNF with the QD algorithms MAP-
Elites (line), CMA-ME, and CMA-MAE, and with NS using CMA-ES as the underlying
optimizer. MAP-Elites (line) describes the implementation of MAP-Elites with
the Iso+LineDD operator (Vassiliades and Mouret,, 2018). All algorithms are
implemented with pyribs (Tjanaka et al.,, 2023).

Our experiments include canonical benchmark domains from QD and DO: Linear
Projection (Fontaine et al.,, 2020), Arm Repertoire (Vassiliades and Mouret,,
2018), and Deceptive Maze (Lehman and Stanley,, 2011). As discussed in Sec. 2,
we convert QD domains into DO domains by setting the objective to be constant,
effectively removing the importance of solution quality. For Deceptive Maze,
we use the implementation in Kheperax (Grillotti and Cully,, 2023).
Furthermore, to address the challenges posed by high-dimensional feature
spaces, we introduce a new domain, Multi-feature Linear Projection, that
generalizes Linear Projection to high-dimensional feature spaces.

All experiments run on a 128-core workstation with an NVIDIA RTX A6000 GPU.
While all algorithms are single-threaded, we use the GPU for training the CNF
in DDS-CNF and for evaluating the Deceptive Maze domain.

###  6.1 Experimental Design

Independent Variable. In each domain, we conduct a between-groups study with
the algorithm as the independent variable.

Dependent Variables. We compute the archive coverage as the number of occupied
cells in the archive divided by the total number of cells. To compare the
coverage of tessellation-based algorithms (CMA-MAE, CMA-ME, and MAP-Elites)
with that of non-tessellation-based algorithms (DDS-KDE, DDS-CNF, and NS), we
track the coverage of all algorithms on a passive archive $\mathcal{A}$
tessellated by a $100\times 100$ grid. For the high-dimensional feature space
experiment, we track coverage with a centroidal Voronoi tessellation with
10,000 cells (Vassiliades et al.,, 2018).

Following prior work (Gomes et al.,, 2015), we also assess the ability of each
algorithm to uniformly explore the feature space. Thus, we measure the cross-
entropy between a uniform distribution and the distribution of sampled
features. Let $N_{e}$ be the total times throughout the entire search that
solutions were discovered in cell $e\in\\{1,\dots,l\\}$ in the passive
archive, and $N_{\mathrm{total}}=\sum_{e=1}^{l}N_{e}$. The cross-entropy score
is defined as:

$\mathrm{CE}=-\sum_{e=1}^{l}\frac{1}{l}\log\left(\frac{N_{e}}{N_{\mathrm{total}}}\right)$
(4)

$CE$ achieves its minimum value when $N(\cdot)$ is uniformly distributed. For
all passive archives in our experiments, this minimum is 9.21 for $l=$10,000$$
cells.

###  6.2 Domain Details

Linear Projection (LP). LP is a QD domain where ${\bm{\phi}}$ induces high
distortion by mapping an $n$-dimensional solution space to a 2D feature space
(Fontaine et al.,, 2020). A harsh penalty is applied outside the bounds of the
feature space to hinder exploration near the bounds. The feature function
${\bm{\phi}}:\mathbb{R}^{n}\to\left[-5.12\cdot\frac{n}{2},5.12\cdot\frac{n}{2}\right]^{2}$
is defined as:

$\displaystyle{\bm{\phi}}({\bm{\theta}})$
$\displaystyle=\left(\sum_{i=1}^{\frac{n}{2}}\mathrm{clip}(\theta_{i}),\sum_{i=\frac{n}{2}+1}^{n}\mathrm{clip}(\theta_{i})\right)$
(5) $\displaystyle\mathrm{clip}(\theta_{i})$
$\displaystyle=\begin{cases}\theta_{i}&\text{if}\;|\theta_{i}|\leq 5.12\\\
5.12/\theta_{i}&\text{otherwise}\\\ \end{cases}$ (6)

where $\theta_{i}$ is the $i$th component of ${\bm{\theta}}$, and we assume
that $n$ is divisible by $2$. ${\bm{\phi}}({\bm{\theta}})$ applies
$\mathrm{clip}(\cdot)$ to each $\theta_{i}$ and sums the two halves of
${\bm{\theta}}$. Since $\mathrm{clip}(\theta_{i})$ restricts $\theta_{i}$ to
the interval $[-5.12,5.12]$, ${\bm{\phi}}({\bm{\theta}})$ is bounded by the
closed interval $\left[-5.12\cdot\frac{n}{2},5.12\cdot\frac{n}{2}\right]^{2}$.

For DO experiments, we set the objective function $f(\bm{\theta})=1$.
Following prior work (Fontaine and Nikolaidis,, 2021, 2023), we let the
solution dimension be $n=100$.

Arm Repertoire. The goal of Arm Repertoire (Cully and Demiris,, 2017;
Vassiliades and Mouret,, 2018) is to search for a diverse collection of arm
positions for a planar robotic arm with $n$ revolving joints. In this domain,
${\bm{\theta}}\in[-\pi,\pi]^{n}$ represents the angles of the $n$ joints, and
${\bm{\phi}}({\bm{\theta}})$ computes the $(x,y)$ position of the arm’s end-
effector using forward kinematics. While all other domains in this paper have
a maximum of $100\%$ coverage, Arm Repertoire has a maximum coverage of
$80.24\%$ when using a $100\times 100$ grid archive (Fontaine and Nikolaidis,,
2023), since the arm can only move in a circle of radius $n$. Similar to LP,
we set $f({\bm{\theta}})=1$ and $n=100$.

(a) The effect of the bandwidth $h$ on a one-dimensional KDE. Red dots
represent the data, and black lines depict the KDE. When $h$ is too small, the
KDE reveals many misleading local maxima. When $h$ is too large, the KDE
conceals modes from the underlying distribution.

(b) The layout in the deceptive maze domain. The red dot indicates the start,
and the green dot indicates the goal. While there is a goal, DO seeks to find
agents that reach every position in the maze. Adapted from prior work (Lehman
and Stanley,, 2011).

Figure 2:

Deceptive Maze. Deceptive Maze is a DO domain that challenges the algorithm to
discover a diverse set of final positions for robots navigating a maze (2(b))
(Lehman and Stanley,, 2011). In this domain, ${\bm{\theta}}$ parameterizes the
robot’s neural network controller. ${\bm{\phi}}({\bm{\theta}})$ is the final
position of the robot after evaluating its path in the maze. As a DO domain,
this has no objective function. In our experiments, the neural network is a
MLP with $n=66$ parameters.

Multi-feature Linear Projection (Multi-feature LP). We introduce a generalized
version of LP that scales to $m$-dimensional feature spaces. The feature
function
${\bm{\phi}}:\mathbb{R}^{n}\to\left[-5.12\cdot\frac{n}{m},5.12\cdot\frac{n}{m}\right]^{m}$
is defined as

${\bm{\phi}}({\bm{\theta}})=\left(\sum_{i=\frac{jn}{m}+1}^{\frac{(j+1)n}{m}}\mathrm{clip}(\theta_{i}):j\in\\{0,\dots,m\\}\right)$
(7)

where we assume that $n$ is divisible by $m$. When $m=2$, this is equivalent
to LP definition in Eq. 5. Our experiments use $n=100$ and $m=10$.

###  6.3 Analysis

| LP | Arm Repertoire | Deceptive Maze | Multi-feature LP
---|---|---|---|---
| Coverage | Cross-Entropy | Coverage | Cross-Entropy | Coverage | Cross-Entropy | Coverage | Cross-Entropy
DDS-KDE | 67.67 ±2.13% | 17.57 ±0.41 | 80.22 ±0.01% | 14.14 ±0.02 | 94.93 ±0.84% | 11.84 ±0.21 | 50.22 ±0.45% | 22.20 ±0.09
DDS-CNF | 63.65 ±1.38% | 18.50 ±0.29 | 79.82 ±0.19% | 13.48 ±0.12 | 85.17 ±2.28% | 14.01 ±0.38 | 16.57 ±0.22% | 28.07 ±0.03
NS | 22.96 ±0.59% | 24.78 ±0.06 | 78.67 ±0.23% | 15.32 ±0.07 | 77.70 ±2.47% | 15.78 ±0.50 | 10.36 ±0.16% | 28.39 ±0.03
CMA-MAE | 56.47 ±1.06% | 19.09 ±0.23 | 77.77 ±0.10% | 13.14 ±0.01 | 98.83 ±0.18% | 10.55 ±0.06 | 11.68 ±0.37% | 29.06 ±0.06
CMA-ME | 42.90 ±0.25% | 22.45 ±0.05 | 39.29 ±0.24% | 20.13 ±0.09 | 79.97 ±2.17% | 15.33 ±0.39 | 0.56 ±0.14% | 28.10 ±0.28
MAP-Elites (line) | 42.30 ±0.31% | 22.43 ±0.07 | 74.62 ±0.07% | 13.60 ±0.01 | 90.32 ±0.87% | 13.28 ±0.16 | 6.28 ±0.34% | 29.61 ±0.02

Figure 3: Coverage and cross-entropy (CE) after 5,000 iterations of each
algorithm in all domains. We report the mean over 10 trials, with error bars
showing the standard error of the mean. Higher coverage and lower cross-
entropy are better.

Fig. 3 shows the mean coverage and cross-entropy over 10 trials. In each
domain, we conducted a one-way ANOVA to check if the algorithms differed in
their coverage and cross-entropy. Since all ANOVAs were significant
($p<0.001$), we followed up with pairwise comparisons via Tukey’s HSD test. We
include full statistical analysis in Appendix F.

Coverage. DDS-KDE and DDS-CNF outperform all baselines on LP, Arm Repertoire,
and Multi-feature LP. They exhibit no statistical difference in performance
with the best-performing algorithm on Deceptive Maze (CMA-MAE). Notably, DDS-
KDE solves Arm Repertoire nearly perfectly, as the maximum coverage in the
domain is $80.24\%$.

We attribute the high coverage of DDS-KDE and DDS-CNF on these domains to the
continuity of their density estimate, which prevents DDS-KDE from converging
prematurely. For example, on LP, our algorithms discover more solutions near
the edges of the feature space than CMA-MAE (Appendix G). The passive archive
maintained by CMA-MAE converges as all the solutions fall in previously
discovered cells (Fontaine et al.,, 2020). In contrast, the continuity of our
density estimate and the online buffer updates facilitate DDS algorithms to
achieve slow but continual progress in exploring the feature space (Fig. 3, LP
and Arm Repertoire). This is because the shape of the continuous density
estimate always changes slightly after each iteration of the algorithm.

For multi-feature LP, we observe that algorithms leveraging continuous
representations of the feature space, i.e., DDS and NS, explore the feature
space much faster than other algorithms driven by discrete feature space
representations, e.g., CMA-MAE (Fig. 3). This is because CMA-MAE is optimizing
on a centroidal Voronoi tessellation with 10,000 cells, where each cell has
significantly more volume compared to that of a $100\times 100$ grid due to
the increased dimensionality of the feature space. Hence, more solutions map
to the same cells, making it more difficult to find solutions outside of
previously explored cells. These results show that DO algorithms based on
continuous density estimates scale better with the dimensionality of the
feature space.

Cross-Entropy. DDS-KDE and DDS-CNF outperform111Recall that lower CE is better
as it indicates a more uniform exploration of the feature space (Sec. 6.1) all
baselines in LP, with the exception of DDS-CNF, whose performance is not
significantly different from CMA-MAE. On Arm Repertoire and Deceptive Maze,
DDS-KDE and DDS-CNF outperform NS. However, while DDS-KDE is on par with CMA-
MAE on Arm Repertoire , CMA-MAE outperforms DDS-KDE on Deceptive Maze and DDS-
CNF on both Arm Repertoire and Deceptive Maze. Finally, DDS-KDE outperforms
all algorithms on Multi-feature LP, while DDS-CNF outperforms CMA-MAE and MAP-
Elites (line).

We attribute the improved performance of CMA-MAE to the nature of the cross-
entropy metric. Cross-entropy is intended to approximate the distributional
similarity between the exploration of the feature space and the uniform
distribution. CMA-MAE directly optimizes a passive archive with uniform
tessellation, unlike DDS and NS, which are agnostic to the passive archive.
This experimental setup naturally favors CMA-MAE with respect to the cross-
entropy metric.

DDS-KDE performs as well as DDS-CNF in terms of coverage and cross-entropy
across most domains. However, on Multi-feature LP, DDS-KDE outperforms DDS-CNF
on both metrics; on LP and Deceptive Maze, DDS-KDE outperforms DDS-CNF in
terms of cross-entropy.

We attribute the difference between the performance of DDS-KDE and DDS-CNF to
the bandwidth parameter in KDE, which allows us to adjust the smoothness of
the KDE (2(a)). On the other hand, CNF lacks explicit control over its
smoothness. Thus, while we adjust the smoothness of the KDE to improve
performance for DDS-KDE (Sec. 5), we can not apply the same techniques to
boost the performance of DDS-CNF.

##  7 Limitations and Future Work

We introduce two new diversity optimization algorithms, DDS-KDE and DDS-CNF.
We provide theoretical insight into the connection between DDS-KDE and novelty
search and show that both DDS algorithms excel at discovering diverse
solutions.

We envision several extensions to DDS. While we proposed two DDS variants, DDS
is a general method that can integrate with a wide variety of distribution
fitting techniques. DDS can be combined with parametric estimators like
mixture models (McLachlan et al.,, 2019) and moment matching (Han et al.,,
2018), as well as with non-parametric estimators like vector quantization
(Kohonen,, 1990) and generative models (Goodfellow et al.,, 2020). In general,
we believe DDS will interface with any method that learns a continuous density
representation of the feature space.

For DDS-KDE, while we assumed a constant bandwidth, prior work has adapted the
bandwidth online (Kristan et al.,, 2011) and varied the bandwidth based on
query location (Terrell and Scott,, 1992). These methods could be applied to
secure tighter uniform stability bounds for the KDE, making feature space
exploration more efficient.

A key limitation of our algorithms is their significant computational cost.
Querying the KDE incurs a runtime of $O(|{\mathcal{B}}|)$, where
$|{\mathcal{B}}|$ is the buffer size, and fine-tuning the CNF is
computationally expensive due to backpropagation, even when accelerated on a
GPU. These issues could be mitigated by accelerating KDE computation (O’Brien
et al.,, 2016; Yang et al.,, 2003) and CNF training efficiency (Onken et al.,,
2021; Huang and Yeh,, 2021).

Finally, we are excited to see how the underlying insight of DDS — leveraging
continuous density representations of the feature space — can improve the
exploration power of QD algorithms, especially in high-dimensional feature
spaces. Applying DDS to QD will improve its performance on applications such
as scenario generation for complex systems (Fontaine et al.,, 2021),
generative design (Gaier et al.,, 2018), damage recovery in robotics (Cully et
al.,, 2015), reinforcement learning (Batra et al.,, 2023), and procedural
content generation (Gravina et al.,, 2019).

## References

  * Akimoto et al., (2010) Akimoto, Y., Nagata, Y., Ono, I., and Kobayashi, S. (2010).  Bidirectional relation between cma evolution strategies and natural evolution strategies.  In International Conference on Parallel Problem Solving from Nature, pages 154–163. Springer.
  * Batra et al., (2023) Batra, S., Tjanaka, B., Fontaine, M. C., Petrenko, A., Nikolaidis, S., and Sukhatme, G. (2023).  Proximal policy gradient arborescence for quality diversity reinforcement learning.  arXiv preprint arXiv:2305.13795.
  * Bowman, (1984) Bowman, A. W. (1984).  An alternative method of cross-validation for the smoothing of density estimates.  Biometrika, 71(2):353–360.
  * Chen, (2017) Chen, Y.-C. (2017).  A tutorial on kernel density estimation and recent advances.  Biostatistics & Epidemiology, 1(1):161–187.
  * Chenu et al., (2021) Chenu, A., Perrin-Gilbert, N., Doncieux, S., and Sigaud, O. (2021).  Selection-expansion: A unifying framework for motion-planning and diversity search algorithms.  In Artificial Neural Networks and Machine Learning–ICANN 2021: 30th International Conference on Artificial Neural Networks, Bratislava, Slovakia, September 14–17, 2021, Proceedings, Part IV 30, pages 568–579. Springer.
  * Cook and Weisberg, (1980) Cook, R. D. and Weisberg, S. (1980).  Characterizations of an empirical influence function for detecting influential cases in regression.  Technometrics, 22(4):495–508.
  * Cully et al., (2015) Cully, A., Clune, J., Tarapore, D., and Mouret, J.-B. (2015).  Robots that can adapt like animals.  Nature, 521(7553):503–507.
  * Cully and Demiris, (2017) Cully, A. and Demiris, Y. (2017).  Quality and diversity optimization: A unifying modular framework.  IEEE Transactions on Evolutionary Computation, 22(2):245–259.
  * Fontaine and Nikolaidis, (2021) Fontaine, M. and Nikolaidis, S. (2021).  Differentiable quality diversity.  Advances in Neural Information Processing Systems, 34.
  * Fontaine and Nikolaidis, (2023) Fontaine, M. and Nikolaidis, S. (2023).  Covariance matrix adaptation map-annealing.  In Proceedings of the Genetic and Evolutionary Computation Conference, pages 456–465.
  * Fontaine et al., (2021) Fontaine, M. C., Hsu, Y.-C., Zhang, Y., Tjanaka, B., and Nikolaidis, S. (2021).  On the importance of environments in human-robot coordination.  Proceedings of Robotics: Science and Systems 17.
  * Fontaine et al., (2020) Fontaine, M. C., Togelius, J., Nikolaidis, S., and Hoover, A. K. (2020).  Covariance matrix adaptation for the rapid illumination of behavior space.  In Proceedings of the 2020 genetic and evolutionary computation conference, pages 94–102.
  * Gaier et al., (2018) Gaier, A., Asteroth, A., and Mouret, J.-B. (2018).  Data-efficient design exploration through surrogate-assisted illumination.  Evolutionary computation, 26(3):381–410.
  * Gaier et al., (2019) Gaier, A., Asteroth, A., and Mouret, J.-B. (2019).  Are quality diversity algorithms better at generating stepping stones than objective-based search?  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’19, page 115–116, New York, NY, USA. Association for Computing Machinery.
  * Glasmachers et al., (2010) Glasmachers, T., Schaul, T., Yi, S., Wierstra, D., and Schmidhuber, J. (2010).  Exponential natural evolution strategies.  In Proceedings of the 12th annual conference on Genetic and evolutionary computation, pages 393–400.
  * Gomes et al., (2015) Gomes, J., Mariano, P., and Christensen, A. L. (2015).  Devising effective novelty search algorithms: A comprehensive empirical study.  In Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, GECCO ’15, page 943–950, New York, NY, USA. Association for Computing Machinery.
  * Goodfellow et al., (2020) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. (2020).  Generative adversarial networks.  Communications of the ACM, 63(11):139–144.
  * Gravina et al., (2019) Gravina, D., Khalifa, A., Liapis, A., Togelius, J., and Yannakakis, G. N. (2019).  Procedural content generation through quality diversity.  In 2019 IEEE Conference on Games (CoG), pages 1–8. IEEE.
  * Grillotti and Cully, (2023) Grillotti, L. and Cully, A. (2023).  Kheperax: A lightweight jax-based robot control environment for benchmarking quality-diversity algorithms.  In Proceedings of the Companion Conference on Genetic and Evolutionary Computation, GECCO ’23 Companion, page 2163–2165, New York, NY, USA. Association for Computing Machinery.
  * Han et al., (2018) Han, Y., Jiao, J., and Weissman, T. (2018).  Local moment matching: A unified methodology for symmetric functional estimation and distribution estimation under wasserstein distance.  In Conference On Learning Theory, pages 3189–3221. PMLR.
  * Hansen, (2016) Hansen, N. (2016).  The cma evolution strategy: A tutorial.  arXiv preprint arXiv:1604.00772.
  * Hardt et al., (2016) Hardt, M., Recht, B., and Singer, Y. (2016).  Train faster, generalize better: Stability of stochastic gradient descent.  In International conference on machine learning, pages 1225–1234. PMLR.
  * Huang and Yeh, (2021) Huang, H.-H. and Yeh, M.-Y. (2021).  Accelerating continuous normalizing flow with trajectory polynomial regularization.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 7832–7839.
  * Kingma and Ba, (2015) Kingma, D. P. and Ba, J. (2015).  Adam: A method for stochastic optimization.  In Bengio, Y. and LeCun, Y., editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015.
  * Kohonen, (1990) Kohonen, T. (1990).  Improved versions of learning vector quantization.  In 1990 ijcnn international joint conference on Neural networks, pages 545–550. IEEE.
  * Kristan et al., (2011) Kristan, M., Leonardis, A., and Skočaj, D. (2011).  Multivariate online kernel density estimation with gaussian kernels.  Pattern Recognition, 44(10):2630–2642.  Semi-Supervised Learning for Visual Content Analysis and Understanding.
  * Lehman and Stanley, (2011) Lehman, J. and Stanley, K. O. (2011).  Abandoning objectives: Evolution through the search for novelty alone.  Evolutionary computation, 19(2):189–223.
  * Li, (1994) Li, K.-H. (1994).  Reservoir-sampling algorithms of time complexity o(n(1 + log(n/n))).  ACM Trans. Math. Softw., 20(4):481–493.
  * Lipman et al., (2023) Lipman, Y., Chen, R. T. Q., Ben-Hamu, H., Nickel, M., and Le, M. (2023).  Flow matching for generative modeling.  In The Eleventh International Conference on Learning Representations.
  * McLachlan et al., (2019) McLachlan, G. J., Lee, S. X., and Rathnayake, S. I. (2019).  Finite mixture models.  Annual review of statistics and its application, 6:355–378.
  * Mouret and Clune, (2015) Mouret, J.-B. and Clune, J. (2015).  Illuminating search spaces by mapping elites.  arXiv preprint arXiv:1504.04909.
  * O’Brien et al., (2016) O’Brien, T. A., Kashinath, K., Cavanaugh, N. R., Collins, W. D., and O’Brien, J. P. (2016).  A fast and objective multidimensional kernel density estimation method: fastKDE.  Computational Statistics & Data Analysis, 101:148–160.
  * Onken et al., (2021) Onken, D., Fung, S. W., Li, X., and Ruthotto, L. (2021).  Ot-flow: Fast and accurate continuous normalizing flows via optimal transport.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 9223–9232.
  * Pugh et al., (2016) Pugh, J. K., Soros, L. B., and Stanley, K. O. (2016).  Quality diversity: A new frontier for evolutionary computation.  Frontiers in Robotics and AI, 3:40.
  * Pugh et al., (2015) Pugh, J. K., Soros, L. B., Szerlip, P. A., and Stanley, K. O. (2015).  Confronting the challenge of quality diversity.  In Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, GECCO ’15, page 967–974, New York, NY, USA. Association for Computing Machinery.
  * Rudemo, (1982) Rudemo, M. (1982).  Empirical choice of histograms and kernel density estimators.  Scandinavian Journal of Statistics, 9(2):65–78.
  * Silverman, (1986) Silverman, B. W. (1986).  Density estimation for statistics and data analysis, volume 26.  CRC press.
  * Terrell and Scott, (1992) Terrell, G. R. and Scott, D. W. (1992).  Variable Kernel Density Estimation.  The Annals of Statistics, 20(3):1236 – 1265.
  * Tjanaka et al., (2023) Tjanaka, B., Fontaine, M. C., Lee, D. H., Zhang, Y., Balam, N. R., Dennler, N., Garlanka, S. S., Klapsis, N. D., and Nikolaidis, S. (2023).  Pyribs: A bare-bones python library for quality diversity optimization.  In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’23, page 220–229, New York, NY, USA. Association for Computing Machinery.
  * Tjanaka et al., (2022) Tjanaka, B., Fontaine, M. C., Togelius, J., and Nikolaidis, S. (2022).  Approximating gradients for differentiable quality diversity in reinforcement learning.  In Proceedings of the Genetic and Evolutionary Computation Conference, pages 1102–1111.
  * Vassiliades et al., (2018) Vassiliades, V., Chatzilygeroudis, K., and Mouret, J.-B. (2018).  Using centroidal voronoi tessellations to scale up the multidimensional archive of phenotypic elites algorithm.  IEEE Transactions on Evolutionary Computation, 22(4):623–630.
  * Vassiliades and Mouret, (2018) Vassiliades, V. and Mouret, J.-B. (2018).  Discovering the elite hypervolume by leveraging interspecies correlation.  In Proceedings of the Genetic and Evolutionary Computation Conference, pages 149–156.
  * Weglarczyk, (2018) Weglarczyk, S. (2018).  Kernel density estimation and its application.  ITM Web Conf., 23:00037.
  * Yang et al., (2003) Yang, Duraiswami, and Gumerov (2003).  Improved fast gauss transform and efficient kernel density estimation.  In Proceedings ninth IEEE international conference on computer vision, pages 664–671. IEEE.

##  Appendix A Hyperparameters

We provide an overview of the hyperparameters of our experiments. Table 1
lists hyperparameters that varies across algorithms and domains.
Hyperparameters that are fixed for the algorithms are described below.

Following prior work (Fontaine et al.,, 2020), all algorithms run with 15
emitters — emitters (Tjanaka et al.,, 2023) are objects that generate
solutions for QD algorithms. For MAP-Elites (line), the emitter applies the
Iso+LineDD operator to random solutions from the archive (Vassiliades and
Mouret,, 2018). For other algorithms, the emitters are individual instances of
CMA-ES. In each experiment, we run the algorithm for 5000 iterations, where
each emitter generates $\lambda=36$ solutions per iteration. Following prior
work (Fontaine and Nikolaidis,, 2023), CMA-ES is configured with the $\mu$
selection rule, basic restart rule, and an initial solution
${\bm{x}}_{0}={\bm{0}}$.

For LP, Arm Repertoire, and Deceptive Maze, MAP-Elites-based algorithms (CMA-
MAE, CMA-ME, and MAP-Elites (line)) run on a grid archive of $100\times 100$
following (Fontaine and Nikolaidis,, 2023). For Multi-feature LP, since
maintaining a grid tessellation would require an exponential number of cells
in the number of dimensions, we use a CVT archive with 10,000 cells. Due to
randomness in the generation of the CVT cells, we pre-generate the centroids
and use the same centroids for all algorithms.

Our DDS algorithms (DDS-KDE and DDS-CNF) maintain a buffer of capacity 10,000.
DDS-KDE uses a Gaussian kernel with a bandwidth parameter of
$h=25.6,10,0.01,2.56$ for Linear Projection, Arm Repertoire, Deceptive Maze,
and Multi-feature LP, respectively. The bandwidth is found through empirical
fine-tuning (Appendix B).

Table 1: Detailed hyperparameter specifications for each algorithm and domain. Algorithm | Parameter | LP | Arm Repertoire | Deceptive Maze | Multi-feature LP
---|---|---|---|---|---
DDS-KDE | Initial step size $\sigma_{0}$ | 1.5 | 0.5 | 1.5 | 1.5
Bandwidth $h$ | 25.6 | 10 | 0.01 | 2.56
DDS-CNF | Initial step size $\sigma_{0}$ | 0.5 | 1.5 | 1.5 | 0.5
NS | Initial step size $\sigma_{0}$ | 0.5 | 0.5 | 0.5 | 0.5
Nearest neighbors $k$ | 100 | 100 | 100 | 100
Acceptance threshold | $0.042\cdot 512$ | $0.042\cdot 200$ | $0.042$ | $0.042\cdot 51.2$
CMA-MAE | Initial step size $\sigma_{0}$ | 0.5 | 0.2 | 0.5 | 0.5
Threshold min | 0 | 0 | 0 | 0
Learning rate $\alpha$ | 0.01 | 0.01 | 0.01 | 0.01
CMA-ME | Initial step size $\sigma_{0}$ | 0.5 | 0.2 | 0.5 | 0.5
MAP-Elites (line) | Isotropic variance $\sigma_{\mathrm{iso}}$ | 0.5 | 0.1 | 0.1 | 0.5
Line variance $\sigma_{\mathrm{line}}$ | 0.2 | 0.2 | 0.2 | 0.2

##  Appendix B Bandwidth Selection for DDS-KDE

The bandwidth parameter $h$ in KDE controls the smoothness of the estimation
(2(a)). While prior work (Silverman,, 1986; Rudemo,, 1982; Bowman,, 1984) has
developed bandwidth selection methods for accurate density estimations, the
optimal bandwidth for DDS-KDE may require an over-smoothing of the KDE (Sec.
5).

To examine the effect of the bandwidth parameter on the performance of DDS-
KDE, we perform ablation experiments on LP, Arm Repertoire, and Deceptive
Maze. While applying the same hyperparameters as described in Appendix A, we
vary the bandwidth $h\in(0,W]$, where $W$ is the bound of the feature space
defined for each domain (Sec. 6.2). For instance, in the LP domain,
$W_{\mathrm{LP}}=5.12\cdot\frac{n}{2}$ where $n$ is the dimensionality of the
solution space. For conciseness and consistency, we report the results with
respect to the normalized bandwidth $h_{0}=\frac{h}{W_{\mathrm{domain}}}$ so
that $h_{0}\in(0,1]$, where $W_{\mathrm{domain}}$ is the feature space bound
for the particular domain.

The experimental results shown in Fig. 4 highlight a peak across all domains
around $h_{0}=0.05$. For LP, Deceptive Maze, and Multi-feature LP, the
performance decays as the bandwidth increases. However, for Arm Repertoire,
DDS-KDE achieves near-optimal coverage for all $h_{0}\geq 0.01$.

| LP | Arm Repertoire | Deceptive Maze | Multi-feature LP
---|---|---|---|---
| Coverage | Cross-Entropy | Coverage | Cross-Entropy | Coverage | Cross-Entropy | Coverage | Cross-Entropy
$h_{0}=0.01$ | 44.57 ±0.75% | 21.99 ±0.14 | 75.75 ±1.06% | 15.07 ±0.22 | 94.93 ±0.84% | 11.84 ±0.21 | 9.50 ±0.15% | 29.33 ±0.02
$h_{0}=0.035$ | 64.00 ±1.82% | 18.10 ±0.37 | 79.84 ±0.20% | 14.39 ±0.09 | 92.22 ±1.44% | 12.56 ±0.37 | 33.17 ±0.26% | 25.57 ±0.05
$h_{0}=0.05$ | 67.67 ±2.13% | 17.57 ±0.41 | 80.22 ±0.01% | 14.14 ±0.02 | 87.70 ±1.60% | 13.98 ±0.35 | 50.22 ±0.45% | 22.20 ±0.09
$h_{0}=0.1$ | 52.41 ±1.18% | 20.75 ±0.22 | 80.20 ±0.01% | 14.14 ±0.04 | 69.63 ±8.59% | 17.60 ±1.69 | 43.61 ±0.81% | 23.13 ±0.15
$h_{0}=0.2$ | 41.90 ±1.69% | 22.50 ±0.33 | 80.15 ±0.00% | 13.99 ±0.02 | 57.95 ±5.43% | 20.41 ±1.03 | 14.88 ±0.18% | 28.08 ±0.02
$h_{0}=0.3$ | 43.48 ±1.79% | 22.22 ±0.32 | 80.10 ±0.01% | 14.01 ±0.02 | 41.98 ±5.51% | 23.54 ±0.98 | 9.53 ±0.12% | 27.97 ±0.03
$h_{0}=0.4$ | 39.37 ±2.00% | 22.92 ±0.40 | 80.10 ±0.01% | 14.00 ±0.01 | 29.77 ±2.80% | 25.81 ±0.47 | 9.24 ±0.09% | 27.78 ±0.03
$h_{0}=0.5$ | 35.74 ±1.49% | 23.68 ±0.28 | 80.09 ±0.00% | 13.99 ±0.02 | 27.90 ±2.12% | 26.14 ±0.37 | 9.39 ±0.19% | 27.74 ±0.04
$h_{0}=0.6$ | 33.19 ±2.39% | 24.17 ±0.44 | 80.10 ±0.01% | 13.99 ±0.01 | 25.04 ±0.24% | 26.63 ±0.03 | 9.42 ±0.11% | 27.78 ±0.02
$h_{0}=0.7$ | 36.20 ±1.31% | 23.52 ±0.28 | 80.09 ±0.01% | 13.99 ±0.01 | 30.23 ±4.19% | 25.73 ±0.77 | 9.02 ±0.16% | 27.78 ±0.03
$h_{0}=0.8$ | 31.35 ±2.65% | 24.46 ±0.54 | 80.09 ±0.00% | 14.00 ±0.01 | 25.74 ±0.32% | 26.53 ±0.06 | 9.29 ±0.17% | 27.77 ±0.04
$h_{0}=0.9$ | 33.66 ±1.61% | 24.00 ±0.30 | 80.10 ±0.00% | 13.98 ±0.01 | 31.08 ±3.55% | 25.56 ±0.64 | 9.13 ±0.18% | 27.82 ±0.04
$h_{0}=1$ | 34.69 ±2.13% | 23.75 ±0.41 | 80.10 ±0.00% | 14.02 ±0.01 | 25.34 ±0.27% | 26.61 ±0.04 | 9.23 ±0.13% | 27.77 ±0.03

Figure 4: Coverage and cross-entropy (CE) after 5,000 iterations of DDS-KDE in
all domains for each normalized bandwidth $h_{0}$. We report the mean over 10
trials (3 trials for Deceptive Maze) with error bars showing the standard
error of the mean. Higher coverage and lower cross-entropy are better.
Highlighted cells are results from the main experiments in Fig. 3.

##  Appendix C Conversion Formula for Affine Transformations of the Feature
Space

We derive a formula that allows one to empirically find an optimal bandwidth
on a domain, then apply the formula to compute the optimal bandwidths for any
affine transformations of that domain.

Consider feature mappings ${\bm{\phi}}_{1}(\cdot)$ and
${\bm{\phi}}_{2}(\cdot)$. Let ${\bm{\phi}}_{2}(\cdot)$ be an affine
transformation of ${\bm{\phi}}_{1}(\cdot)$, i.e.,
${\bm{\phi}}_{2}({\bm{x}})={\bm{A}}{\bm{\phi}}_{1}({\bm{x}})+{\bm{b}}$.
Furthermore, assume there is a KDE $\hat{D}_{h_{1}}(\cdot)$ over
${\bm{\phi}}_{1}(\cdot)$. We prove that when ${\bm{A}}$ is some constant $c$,
$\hat{D}_{h_{2}}$ on ${\bm{\phi}}_{2}(\cdot)$ where $h_{2}=ch_{1}$ is an
identical density function as $\hat{D}_{h_{1}}(\cdot)$, in expectation.

######  Theorem C.1.

Consider $\hat{D}_{h_{0},{\bm{X}}}({\bm{x}};\mathcal{X})$ and
$\hat{D}_{h_{1},{\bm{Y}}}({\bm{y}};\mathcal{Y})$, where
$\hat{D}_{h_{0},{\bm{X}}}$ is a KDE with bandwidth $h_{0}$ modeling the
density of a random vector ${\bm{X}}$ and $\hat{D}_{h_{1},{\bm{Y}}}$ is a KDE
with bandwidth $h_{1}$ modeling the density of a random vector ${\bm{Y}}$. If
${\bm{Y}}=a{\bm{X}}+{\bm{b}}$, then
$\hat{D}_{h_{0},{\bm{X}}}({\bm{x}})=\hat{D}_{h_{1},{\bm{Y}}}({\bm{y}})$ when
$h_{0}=ah_{1}$.

###### Proof.

The density function $p_{X}$ estimated by a KDE for a random vector ${\bm{X}}$
with a given bandwidth $h$ is represented as $\hat{D}_{h,{\bm{X}}}$. Recall
the definition of $\hat{D}_{h,{\bm{X}}}$ according to the definition of a KDE
from Eq. 2

$\hat{D}_{h,{\bm{X}}}({\bm{x}};\mathcal{X})=\frac{1}{|\mathcal{X}|h}\sum_{{\bm{x}}^{\prime}\in\mathcal{X}}K\left(\frac{\lVert{\bm{x}}-{\bm{x}}^{\prime}\rVert}{h}\right)$
(8)

for kernel $K(\cdot)$ and $\mathcal{X}\subseteq\mathbb{R}^{m}$. Consider the
norm as a random variable $R^{\prime}=\lVert{\bm{X}}-{\bm{x}}^{\prime}\rVert$.
Since the kernel function $K(\cdot)$ maintains the property
$\int_{-\infty}^{\infty}K_{h}(u)\,du=1$ (Chen,, 2017), we can define the
probability density function of $R^{\prime}$ to be:

$p_{R^{\prime}}(r^{\prime})=\frac{1}{h}K\left(\frac{r^{\prime}}{h}\right)=\frac{1}{h}K\left(\frac{||{\bm{x}}-{\bm{x}}^{\prime}||}{h}\right)$
(9)

Hence, the density estimate for $X$ can be represented in terms of
$r_{i}=||{\bm{x}}-{\bm{x}}_{i}||$ as shown below:

$\hat{D}_{h,{\bm{X}}}({\bm{x}};\mathcal{X})=\frac{1}{|\mathcal{X}|h}\sum_{{\bm{x}}^{\prime}\in\mathcal{X}}K\left(\frac{\lVert{\bm{x}}-{\bm{x}}^{\prime}\rVert}{h}\right)=\frac{1}{|\mathcal{X}|}\sum_{{\bm{x}}^{\prime}\in\mathcal{X}}p_{R^{\prime}}(r^{\prime})$
(10)

Now, when we transform our random vector ${\bm{X}}$ into
${\bm{Y}}=a{\bm{X}}+{\bm{b}}$, where $a\in\mathbb{R}^{+}$ and
$b\in\mathbb{R}^{m}$, we transform the norm $R^{\prime}$ into $Q^{\prime}$ as
follows:

$\displaystyle R^{\prime}=||{\bm{X}}-{\bm{x}}^{\prime}||,Q^{\prime}$
$\displaystyle=||{\bm{Y}}-{\bm{y}}^{\prime}||$
$\displaystyle=||a{\bm{X}}+{\bm{b}}-\left(a{\bm{x}}^{\prime}+{\bm{b}}\right)||$
$\displaystyle=||a{\bm{X}}-a{\bm{x}}^{\prime}||$
$\displaystyle=a||{\bm{X}}-{\bm{x}}^{\prime}||$ $\displaystyle=aR^{\prime}$

So, we can represent $Q^{\prime}$ as $Q^{\prime}=g(R^{\prime})=aR^{\prime}$.
Since $g$ is a monotonically increasing function, we can represent the density
function of $Q^{\prime}$ as follows:

$\displaystyle p_{Q^{\prime}}(q^{\prime})$
$\displaystyle=p_{R^{\prime}}\left(g^{-1}(q^{\prime})\right)\left|\frac{d}{dq^{\prime}}g^{-1}(q^{\prime})\right|$
$\displaystyle=\left(\frac{1}{a}\right)\left(\frac{1}{h}\right)K\left(\frac{\frac{q^{\prime}}{a}}{h}\right)$
$\displaystyle=\frac{1}{ah}K\left(\frac{q^{\prime}}{ah}\right)$
$\displaystyle=\frac{1}{ah}K\left(\frac{||{\bm{y}}-{\bm{y}}^{\prime}||}{ah}\right)$

We can apply this relationship now to the density function estimated by the
KDE for ${\bm{Y}}$, given that $q^{\prime}=||{\bm{Y}}-{\bm{y}}^{\prime}||$:

$\displaystyle\hat{D}_{h,{\bm{Y}}}({\bm{y}};\mathcal{Y})$
$\displaystyle=\frac{1}{|\mathcal{Y}|}\sum_{{\bm{y}}^{\prime}\in\mathcal{Y}}p_{q^{\prime}}(q^{\prime})$
$\displaystyle=\frac{1}{|\mathcal{Y}|}\sum_{{\bm{y}}^{\prime}\in\mathcal{Y}}\frac{1}{ah}K\left(\frac{||{\bm{y}}-{\bm{y}}^{\prime}||}{ah}\right)$
$\displaystyle=\frac{1}{|\mathcal{Y}|(ah)}\sum_{{\bm{y}}^{\prime}\in\mathcal{Y}}K\left(\frac{||{\bm{y}}-{\bm{y}}^{\prime}||}{ah}\right)$
$\displaystyle=\hat{D}_{ah,{\bm{Y}}}({\bm{y}};\mathcal{Y})$

Hence, if a KDE with a bandwidth of $h$ models the density of a given random
vector ${\bm{X}}$, then a random vector ${\bm{Y}}=a{\bm{X}}+{\bm{b}}$ should
select a bandwidth of $ah$ to maintain a comparable density function. ∎

##  Appendix D Reservoir Sampling for Maintaining a Representative Buffer

As outlined in Sec. 4, DDS maintains a buffer ${\mathcal{B}}$ that stores the
features of sampled solutions. While it would be ideal to maintain an infinite
buffer that stores all features ever encountered, due to runtime and memory
constraints of KDE and CNF, we maintain a finite-sized buffer ${\mathcal{B}}$
of size $|{\mathcal{B}}|=b$. Moreover, the subset of features retained by
${\mathcal{B}}$ must accurately represent all the features ever encountered.
This can be described as the _reservoir sampling_ problem.

######  Problem D.1 (Reservoir Sampling (Li,, 1994)).

Suppose there is a finite population of unknown size $N$. Derive an algorithm
that draws a simple random sample of size $b$ from this population by scanning
it exactly once, with the memory constraint such that the algorithm can
remember at most $b$ samples at a time.

The population consists of the features that are discovered by DDS throughout
its entire execution. Since the buffer ${\mathcal{B}}$ is a simple random
sample of all features, the density function estimated on ${\mathcal{B}}$ will
be close to the density function of all features. Thus, with reservoir
sampling, we are only required to maintain a subset of all discovered features
without severely altering the shape of the density estimation of the feature
space. In DDS, this is achieved by leveraging algorithm L, the optimal
reservoir sampling algorithm (Li,, 1994).

##  Appendix E Proofs

######  Theorem E.1.

A kernel density estimate $\hat{D}_{h}({\bm{x}};{\mathcal{B}})$ managed with
reservoir sampling, such that features in the buffer are exchanged one at a
time, is $\frac{1}{|{\mathcal{B}}|h}$-uniformly stable, where $h$ is the
bandwidth.

###### Proof.

Let buffer ${\mathcal{B}}\subseteq\mathbb{R}^{m}$ and
${\mathcal{B}}^{\prime}={\mathcal{B}}\setminus\\{{\bm{x}}_{i}\\}\cup\\{{\bm{x}}_{j}\\}$
for ${\bm{x}}_{i}\in{\mathcal{B}}$ and ${\bm{x}}_{j}\in\mathbb{R}^{m}$.
Consider the maximum difference between $\hat{D}_{h}({\bm{x}};{\mathcal{B}})$
and $\hat{D}_{h}({\bm{x}};{\mathcal{B}}^{\prime})$:

$\displaystyle\sup_{{\bm{x}}\in\mathbb{R}^{d}}\left|\hat{D}_{h}({\bm{x}};{\mathcal{B}})-\hat{D}_{h}({\bm{x}};{\mathcal{B}}^{\prime})\right|$
$\displaystyle=\sup_{{\bm{x}}\in\mathbb{R}^{m}}\left|\frac{1}{|\mathcal{{\mathcal{B}}}|h}\sum_{{\bm{a}}\in{\mathcal{B}}}K\left(\frac{\lVert{\bm{x}}-{\bm{a}}\rVert}{h}\right)-\frac{1}{|\mathcal{{\mathcal{B}}}|h}\sum_{{\bm{a}}^{\prime}\in{\mathcal{B}}^{\prime}}K\left(\frac{\lVert{\bm{x}}-{\bm{a}}^{\prime}\rVert}{h}\right)\right|$
(11)
$\displaystyle=\frac{1}{|\mathcal{{\mathcal{B}}}|h}\sup_{{\bm{x}}\in\mathbb{R}^{m}}\left|\sum_{{\bm{a}}\in{\mathcal{B}}}K\left(\frac{\lVert{\bm{x}}-{\bm{a}}\rVert}{h}\right)-\sum_{{\bm{a}}^{\prime}\in{\mathcal{B}}^{\prime}}K\left(\frac{\lVert{\bm{x}}-{\bm{a}}^{\prime}\rVert}{h}\right)\right|$
(12)
$\displaystyle=\frac{1}{|\mathcal{{\mathcal{B}}}|h}\sup_{{\bm{x}}\in\mathbb{R}^{m}}\left|K\left(\frac{\lVert{\bm{x}}-{\bm{x}}_{i}\rVert}{h}\right)-K\left(\frac{\lVert{\bm{x}}-{\bm{x}}_{j}\rVert}{h}\right)\right|$
(13) $\displaystyle\leq\frac{1}{|\mathcal{{\mathcal{B}}}|h}$ (14)

Since ${\mathcal{B}}$ and ${\mathcal{B}}^{\prime}$ differ by exactly one
element, each term in the summation for Eq. 12 cancels except for
${\bm{x}}_{i}$ and ${\bm{x}}_{j}$. Eq. 14 is due to the fact that any kernel
$K(\cdot)$ is bounded between $[0,1]$ by definition. ∎

######  Theorem E.2.

Novelty score $\rho({\bm{x}};{\mathcal{B}})$ is $\frac{W}{k}$-uniformly
stable, where $k$ is the nearest neighbors parameter in novelty score and
$W=\max_{{\bm{s}}_{1},{\bm{s}}_{2}\in S}\lVert{\bm{s}}_{1}-{\bm{s}}_{2}\rVert$
is the diameter of the feature space $S$.

###### Proof.

Let buffer ${\mathcal{B}}\subseteq\mathbb{R}^{m}$ and
${\mathcal{B}}^{\prime}={\mathcal{B}}\cup\\{{\bm{x}}_{j}\\}$ for
${\bm{x}}_{j}\in\mathbb{R}^{m}$. Since we constructed ${\mathcal{B}}^{\prime}$
by adding a single element to ${\mathcal{B}}$, $N_{k}({\bm{x}},{\mathcal{B}})$
and $N_{k}({\bm{x}},{\mathcal{B}}^{\prime})$ differ by at most one element
${\bm{x}}_{i}\in{\mathcal{B}}$ and ${\bm{x}}_{j}\in{\mathcal{B}}^{\prime}$.
Hence,

$\displaystyle\sup_{{\bm{x}}\in\mathbb{R}^{d}}\left|\rho({\bm{x}};{\mathcal{B}})-\rho({\bm{x}};{\mathcal{B}}^{\prime})\right|$
$\displaystyle=\sup_{{\bm{x}}\in\mathbb{R}^{d}}\left|\frac{1}{k}\sum_{{\bm{a}}\in
N_{k}({\bm{x}};{\mathcal{B}})}\lVert{\bm{x}}-{\bm{a}}\rVert-\frac{1}{k}\sum_{{\bm{a}}\in
N_{k}({\bm{x}};{\mathcal{B}}^{\prime})}\lVert{\bm{x}}-{\bm{a}}\rVert\right|$
$\displaystyle=\frac{1}{k}\sup_{{\bm{x}}\in\mathbb{R}^{d}}\left|\lVert{\bm{x}}-{\bm{x}}_{i}\rVert-\lVert{\bm{x}}-{\bm{x}}_{j}\rVert\right|$
$\displaystyle\leq\frac{1}{k}\sup_{{\bm{x}}\in\mathbb{R}^{d}}\lVert{\bm{x}}_{i}-{\bm{x}}_{j}\rVert$
$\displaystyle\leq\frac{W}{k}$

∎

######  Theorem E.3.

Let ${\mathcal{B}}\subseteq\mathbb{R}^{m}$ be a set of features. Consider the
rankings $\pi_{\mathrm{NS}}$ and $\pi_{\mathrm{KDE}}$ on another set of
features $\\{{\bm{\phi}}_{1},\dots,{\bm{\phi}}_{m}\\}\subseteq\mathbb{R}^{m}$
where ${\bm{\phi}}_{i}={\bm{\phi}}({\bm{\theta}}_{i})$. We define the rankings
as follows:
$\pi_{\mathrm{NS}}(i)\geq\pi_{\mathrm{NS}}(j)\iff\rho({\bm{\phi}}_{i};{\mathcal{B}})\geq\rho({\bm{\phi}}_{j};\mathcal{X})$
and
$\pi_{\mathrm{KDE}}(i)\geq\pi_{\mathrm{KDE}}(j)\iff\hat{D}({\bm{\phi}}_{i};{\mathcal{B}})\leq\hat{D}({\bm{\phi}}_{j};{\mathcal{B}})$.
We show that $\pi_{NS}=\pi_{KDE}$, when NS has $k=\infty$ (or, equivalently,
$k=|{\mathcal{B}}|$) and KDE has triangular kernel $K(u)=1-|u|$ with support
over the entire feature space $S$.

###### Proof.

Recall the novelty score function,

$\rho({\bm{x}};{\mathcal{B}})=\frac{1}{|{\mathcal{B}}|}\sum_{{\bm{b}}\in{\mathcal{B}}}\lVert{\bm{x}}-{\bm{b}}\rVert$

Consider the kernel density estimator with $h=1$ and triangular kernel,

$\displaystyle\hat{D}({\bm{x}};{\mathcal{B}})$
$\displaystyle=\frac{1}{|{\mathcal{B}}|}\sum_{{\bm{b}}\in{\mathcal{B}}}K\left(\lVert{\bm{x}}-{\bm{b}}\rVert\right)$
$\displaystyle=\frac{1}{|{\mathcal{B}}|}\sum_{{\bm{b}}\in{\mathcal{B}}}(1-\lVert{\bm{x}}-{\bm{b}}\rVert)$
$\displaystyle=1-\frac{1}{|{\mathcal{B}}|}\sum_{{\bm{b}}\in{\mathcal{B}}}\lVert{\bm{x}}-{\bm{b}}\rVert$
$\displaystyle=1-\rho({\bm{x}};{\mathcal{B}})$

Thus,
$\rho({\bm{\phi}}_{i};{\mathcal{B}})\leq\rho({\bm{\phi}}_{j};{\mathcal{B}})\iff\hat{D}({\bm{\phi}}_{i};{\mathcal{B}})\geq\hat{D}({\bm{\phi}}_{j};{\mathcal{B}})$.
∎

##  Appendix F Statistical Analysis

We include results from our statistical analysis in Sec. 6. We performed one-
way ANOVAs in each domain, shown in Table 2. We also performed pairwise
comparisons with Tukey’s HSD, shown in Table 3 and Table 4.

Table 2:  One-way ANOVA results in each domain. All $p$-values are less than 0.001.  | Coverage | Cross-Entropy
---|---|---
LP | $F(5,54)=206.61$ | $F(5,54)=156.04$
Arm Repertoire | $F(5,54)=9691.05$ | $F(5,54)=1523.35$
Deceptive Maze | $F(5,54)=23.89$ | $F(5,54)=39.33$
Multi-feature LP | $F(5,54)=3350$ | $F(5,54)=466.49$
Table 3: Pairwise comparisons for coverage in each domain. Each entry compares the method in the row to the method in the column. For instance, we can see that DDS-KDE had significantly higher coverage than NS in LP. Comparisons were performed with Tukey’s HSD test. $>$ indicates significantly greater, $<$ indicates significantly less, and $-$ indicates no significant difference. $\varnothing$ indicates an invalid comparison. | LP | Arm Repertoire | Deceptive Maze | Multi-dim LP
---|---|---|---|---
|  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line)
DDS-KDE | $\varnothing$ | $-$ | $>$ | $>$ | $>$ | $>$ | $\varnothing$ | $-$ | $>$ | $>$ | $>$ | $>$ | $\varnothing$ | $>$ | $>$ | $-$ | $>$ | $-$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$ | $>$
DDS-CNF | $-$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$ | $-$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$ | $<$ | $\varnothing$ | $>$ | $<$ | $-$ | $-$ | $<$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$
NS | $<$ | $<$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $<$ | $\varnothing$ | $>$ | $>$ | $>$ | $<$ | $<$ | $\varnothing$ | $<$ | $-$ | $<$ | $<$ | $<$ | $\varnothing$ | $<$ | $>$ | $>$
CMA-MAE | $<$ | $<$ | $>$ | $\varnothing$ | $>$ | $>$ | $<$ | $<$ | $<$ | $\varnothing$ | $>$ | $>$ | $-$ | $>$ | $>$ | $\varnothing$ | $>$ | $>$ | $<$ | $<$ | $>$ | $\varnothing$ | $>$ | $>$
CMA-ME | $<$ | $<$ | $>$ | $<$ | $\varnothing$ | $-$ | $<$ | $<$ | $<$ | $<$ | $\varnothing$ | $<$ | $<$ | $-$ | $-$ | $<$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $<$ | $\varnothing$ | $<$
MAP-Elites (line) | $<$ | $<$ | $>$ | $<$ | $-$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $>$ | $\varnothing$ | $-$ | $-$ | $>$ | $<$ | $>$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $>$ | $\varnothing$
Table 4: Pairwise comparisons for cross-entropy in each domain. Note that lower cross-entropy is better, so significantly less ($<$) indicates that a method is significantly better. | LP | Arm Repertoire | Deceptive Maze | Multi-dim LP
---|---|---|---|---
|  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line) |  DDS-KDE |  DDS-CNF |  NS |  CMA-MAE |  CMA-ME |  MAP-Elites (line)
DDS-KDE | $\varnothing$ | $-$ | $<$ | $<$ | $<$ | $<$ | $\varnothing$ | $>$ | $<$ | $>$ | $<$ | $>$ | $\varnothing$ | $<$ | $<$ | $-$ | $<$ | $<$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $<$
DDS-CNF | $-$ | $\varnothing$ | $<$ | $-$ | $<$ | $<$ | $<$ | $\varnothing$ | $<$ | $>$ | $<$ | $-$ | $>$ | $\varnothing$ | $<$ | $>$ | $-$ | $-$ | $>$ | $\varnothing$ | $-$ | $<$ | $-$ | $<$
NS | $>$ | $>$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$ | $>$ | $\varnothing$ | $>$ | $<$ | $>$ | $>$ | $>$ | $\varnothing$ | $>$ | $-$ | $>$ | $>$ | $-$ | $\varnothing$ | $<$ | $-$ | $<$
CMA-MAE | $>$ | $-$ | $<$ | $\varnothing$ | $<$ | $<$ | $<$ | $<$ | $<$ | $\varnothing$ | $<$ | $<$ | $-$ | $<$ | $<$ | $\varnothing$ | $<$ | $<$ | $>$ | $>$ | $>$ | $\varnothing$ | $>$ | $<$
CMA-ME | $>$ | $>$ | $<$ | $>$ | $\varnothing$ | $-$ | $>$ | $>$ | $>$ | $>$ | $\varnothing$ | $>$ | $>$ | $-$ | $-$ | $>$ | $\varnothing$ | $>$ | $>$ | $-$ | $-$ | $<$ | $\varnothing$ | $<$
MAP-Elites (line) | $>$ | $>$ | $<$ | $>$ | $-$ | $\varnothing$ | $<$ | $-$ | $<$ | $>$ | $<$ | $\varnothing$ | $>$ | $-$ | $<$ | $>$ | $<$ | $\varnothing$ | $>$ | $>$ | $>$ | $>$ | $>$ | $\varnothing$

##  Appendix G Heatmaps

Figure 5: Heatmaps of DDS-KDE, DDS-CNF, NS, and CMA-MAE in the three domains
that have 2D feature spaces (LP, Arm Repertoire, Deceptive Maze).