
\section{Introduction}
\label{sec:intro}

In this paper we propose a new quality diversity algorithm that bridges the gap between single-objective optimization and quality diversity optimization~\citep{pugh2015confronting, chatzilygeroudis2021quality}.

Consider an example problem of searching for specific celebrity faces in the latent space of a generative model. As a single-objective optimization problem, we specify an objective $f$ that targets a celebrity such as Beyonc\'e. A single-objective optimizer, such as CMA-ES~\citep{hansen:cma16}, will converge to a single solution of high objective value, an image that looks like Beyonc\'e as much as possible. 

However, this objective has ambiguity. At what age did we want a photo of Beyonc\'e? Did we want the person in the image to have short or long hair? By instead framing the problem as a quality diversity optimization problem, we additionally specify a measure function $m_1$ that quantifies age and a measure function $m_2$ that quantifies hair length. A quality diversity algorithm, such as \mbox{CMA-ME}~\citep{fontaine2020covariance}, can then optimize for a collection of images that are diverse with respect to age and hair length, but all look like Beyonc\'e. While previous work~\citep{fontaine2020covariance, fontaine2020illuminating, earle2021illuminating, pontes2022unified, fontaine2021importance} has shown that CMA-ME solves such QD problems efficiently, the algorithm prioritizes exploration over exploitation, sometimes returning solutions of lower quality in order to explore new solutions that are diverse in measure space, the vector space defined by the measure function outputs.


We propose a new algorithm, CMA-MAE, that smoothly blends between CMA-ES and CMA-ME through a single parameter $\alpha$. To show how to smoothly blend between the two algorithms, first consider how CMA-ES and CMA-ME differ. At each step CMA-ES's objective ranking maximizes the objective function $f$ by approximating the natural gradient of $f$ at the current solution point~\citep{akimoto2010bidirectional}. In contrast, CMA-ME's improvement ranking moves in the direction of the natural gradient of $f-f_A$ at the current solution point, where $f_A$ is a discount function, equal to the objective of the best solution so far that has the same measure values as the current solution point. The function $f-f_A$ quantifies the gap between a candidate solution and the best solution so far at the candidate solution's position in measure space.



Our key insight is to \textit{anneal the function $f_A$ by a learning rate $\alpha$}. We observe that when $\alpha = 0$, then our discount function $f_A$ never increases and our algorithm behaves like CMA-ES. However, when $\alpha = 1$, then our discount function always maintains the best solution for each region in measure space and our algorithm behaves like CMA-ME. For $0 < \alpha < 1$, CMA-MAE smoothly blends between the two algorithms' behavior, allowing for an algorithm that spends more time on the optimization of $f$ before transitioning to exploration. Figure~\ref{fig:dual} is an illustrative example of varying the learning rate $\alpha$.

Our proposed annealing method has an additional benefit. Both CMA-ES and CMA-ME struggle on flat objective functions $f$ as the natural gradient becomes $\bm{0}$ in this case and each algorithm will restart. However, we show that, when CMA-MAE optimizes $f-f_A$ for $0 < \alpha < 1$, the algorithm becomes a density descent method in measure space, allowing for continued exploration in this case.


Our work makes three main contributions. 1) We propose a new QD algorithm Covariance Matrix Adaptation MAP-Annealing (CMA-MAE) that smoothly blends between CMA-ES and CMA-ME.\\ 2) We prove that CMA-MAE blends between CMA-ES and CMA-ME, and that \mbox{CMA-MAE} becomes a density descent method in measure space for any $0 < \alpha < 1$ when optimizing flat objectives.\\ 3) We show in three QD benchmark domains (Linear Projection, Arm Repertoire, and Latent Space Illumination) that CMA-MAE outperforms state-of-the-art QD algorithms, and CMA-MAE maintains robust performance across different archive resolutions and archive learning rates.



\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/dual.pdf}
\caption{An example of how different $\alpha$ values affect the exploration properties of CMA-MAE. In this example $f$ is a bimodal objective, where each modality is separated by measure $m_1$ and mode $A$ is harder to optimize than $B$, requiring more optimization steps. For $\alpha = 0$, the objective $f$ is equivalent to $f-f_A$, as $f_A$ remains constant. As $\alpha$ increases, CMA-MAE places heavier discounts on optimizing region $B$ in favor of more optimization steps being taken in region $A$.}
\label{fig:dual}
\end{figure}



\section{Problem Definition } \label{sec:problem}
\noindent\textbf{Quality Diversity.} We adopt the quality diversity (QD) problem definition from previous work~\citep{fontaine2021differentiable}. A quality diversity problem consists of an objective $f: \mathbb{R}^n \rightarrow \mathbb{R}$ that maps \mbox{$n$-dimensional} solution parameters to a scalar value denoting the quality of the solution and $k$ measures $m_i: \mathbb{R}^n \rightarrow \mathbb{R}$ or, as a vector function, $\bm{m}: \mathbb{R}^n \rightarrow \mathbb{R}^k$ that quantify behavior or attributes of each solution\footnote{In agent-based settings, such as reinforcement learning, the measure functions are sometimes called behavior functions and the outputs of each measure function are called behavioral characteristics or behavior descriptors.}.  The range of $\bm{m}$ forms a measure space $S = \bm{m}(\mathbb{R}^n)$. The QD objective is to find a set of solutions $\bm{\theta} \in \mathbb{R}^n$, such that  $\bm{m}(\bm{\theta}) = \bm{s}$ for each $\bm{s}$ in $S$ and $f(\bm{\theta})$ is maximized.

The measure space $S$ is continuous, but solving algorithms need to produce a finite collection of solutions. Therefore, QD algorithms in the MAP-Elites~\citep{mouret2015illuminating,cully:nature15} family relax the QD objective by discretizing the space $S$. Given $T$ as the tessellation of $S$ into $M$ cells, the QD objective becomes to find a solution $\bm{\theta_i}$ for each of the $i\in \{1,\ldots,M\}$ cells, such that each $\bm{\theta_i}$ maps to the cell corresponding to $\bm{m}(\bm{\theta_i})$ in the tesselation $T$. The QD objective then becomes maximizing the objective value $f(\bm{\theta_i})$ of all cells:

   \vspace{-2em}
\begin{equation}
  \max \sum_{i=1}^M f(\bm{\theta_i})
\label{eq:objective}
\end{equation}

The \textit{differentiable quality diversity} (DQD) problem~\citep{fontaine2021differentiable} is a special case of the QD problem where both the objective $f$ and measures $m_i$ are \textit{first-order differentiable}.











\section{Preliminaries} \label{sec:preliminaries}

We present several QD algorithms that solve derivative-free QD problems to provide context for our proposed CMA-MAE algorithm. Appendix~\ref{sec:cma-maega} contains information about the DQD algorithm CMA-MEGA, which solves problems where exact gradient information is available. 

 

\noindent\textbf{MAP-Elites and MAP-Elites (line)}. The MAP-Elites QD algorithm produces an archive of solutions, where each cell in the archive corresponds to the provided tesselation $T$ in the QD problem definition. The algorithm initializes the archive by sampling solutions from the solution space $\mathbb{R}^n$ from a fixed distribution. After initialization, MAP-Elites produces new solutions by selecting occupied cells uniformly at random and perturbing them with isotropic Gaussian noise, $\bm{\theta}' = \bm{\theta_i} + \sigma \mathcal{N}(\bm{0},I)$. For each new candidate solution $\bm{\theta}'$, the algorithm computes an objective $f(\bm{\theta}')$ and measures $\bm{m}(\bm{\theta}')$. MAP-Elites places $\bm{\theta}'$ into the archive if the cell corresponding to $\bm{m}(\bm{\theta}')$ is empty or $\bm{\theta}'$ obtains a better objective value $f(\bm{\theta}')$ than the current occupant. The MAP-Elites algorithm results in an archive of solutions that are diverse with respect to the measure function $\bm{m}$, but also high quality with respect to the objective $f$. Prior work~\citep{vassiliades:gecco18} proposes the MAP-Elites (line) algorithm by augmenting the isotropic Gaussian perturbation with a linear interpolation between two solutions $ \bm{\theta_i}$ and $ \bm{\theta_j}$: \mbox{$\bm{\theta'} = \bm{\theta_i} + \sigma_1 \mathcal{N}(\bm{0},I) +  \sigma_2 \mathcal{N}(\bm{0},1) (\bm{\theta_i} - \bm{\theta_j})$}.











\noindent\textbf{CMA-ME}. \label{subsec:CMA-ME} 
Covariance Matrix Adaptation MAP-Elites (\mbox{CMA-ME})~\citep{fontaine2020covariance} combines the archiving mechanisms of MAP-Elites with the adaptation mechanisms of CMA-ES~\cite{hansen:cma16}. Instead of perturbing archive solutions with Gaussian noise, CMA-ME maintains a multivariate Gaussian of search directions $\mathcal{N}(\bm{0},\Sigma)$ and a search point $\bm{\theta} \in \mathbb{R}^n$. The algorithm updates the archive by sampling $\lambda$ solutions around the current search point $\bm{\theta_i} \sim \mathcal{N}(\bm{\theta},\Sigma)$. After updating the archive, \mbox{CMA-ME} ranks solutions via a two stage ranking. Solutions that discover a new cell are ranked by the objective $\Delta_i = f(\bm{\theta_i})$, and solutions that map to an occupied cell $e$ are ranked by the improvement over the incumbent solution $\bm{\theta_e}$ in that cell: $\Delta_i = f(\bm{\theta_i}) - f(\bm{\theta_e})$. CMA-ME prioritizes exploration by ranking all solutions that discover a new cell before all solutions that improve upon an existing cell. Finally, CMA-ME moves $\bm{\theta}$ towards the largest improvement in the archive, according to the CMA-ES update rules. Previous work~\citep{fontaine2021differentiable} has shown that the improvement ranking of CMA-ME approximates a natural gradient of a modified QD objective (see Eq.~\ref{eq:objective}). 








\begin{figure}[!t]
\centering
\includegraphics[width=0.9\linewidth]{figures/front.pdf}
\caption{Our proposed Covariance Matrix Adaptation MAP-Annealing (CMA-MAE) algorithm smoothly blends between the behavior of CMA-ES and CMA-ME via an archive learning rate $\alpha$.}
\label{fig:blend}
\end{figure}










\section{Proposed Algorithms} 

We present the CMA-MAE algorithm, which smoothly blends between the single-objective optimizer CMA-ES and the QD algorithm CMA-ME. While we focus on CMA-MAE, the same augmentations apply to CMA-MEGA to form the novel CMA-MAEGA algorithm (see Appendix~\ref{sec:cma-maega}).

\noindent\textbf{CMA-MAE.} CMA-MAE is an algorithm that smoothly blends between CMA-ES and CMA-ME. First, consider at a high level how CMA-ME explores the measure space and discovers high quality solutions. The CMA-ME algorithm maintains a solution point $\bm{\theta}$ and an archive $A$ with previously discovered solutions. When CMA-ME samples a new solution $\bm{\theta}'$, the algorithm computes the solution's objective value $f(\bm{\theta'})$ and maps the solution to a cell $e$ in the archive based on the measure $\bm{m}(\bm{\theta'})$. CMA-ME then computes the improvement of the objective value $f(\bm{\theta'})$ of the new solution, over a discount function $f_A : \mathbb{R}^n \rightarrow \mathbb{R}$. In CMA-ME, we define $f_A(\bm{\theta'})$ by computing the cell $e$ in the archive corresponding to $\bm{m}(\bm{\theta}')$  and letting $f_A(\bm{\theta'}) = f(\bm{\theta_e})$, where $\bm{\theta_e}$ is the incumbent solution of cell $e$. The algorithm ranks candidate solutions by improvement $f(\bm{\theta'})-f_A(\bm{\theta'})=f(\bm{\theta'})-f(\bm{\theta_e})$ and moves the search in the direction of higher ranked solutions.


Assume that CMA-ME samples a new solution $\bm{\theta'}$ with a high objective value of $f(\bm{\theta'}) = 99$. If the current occupant $\bm{\theta_e}$ of the corresponding cell has a low objective value of $f(\bm{\theta_e}) = 0.3$, then the improvement in the archive $\Delta = f(\bm{\theta'}) - f(\bm{\theta_e}) = 98.7$ is high and the algorithm will move the search point $\bm{\theta}$ towards $\bm{\theta'}$. Now, assume that in the next iteration the algorithm discovers a new solution $\bm{\theta''}$ with objective value $f(\bm{\theta''}) = 100$ that maps to the same cell as $\bm{\theta'}$. The improvement then is $\Delta = f(\bm{\theta''}) - f(\bm{\theta'}) = 1$ as $\bm{\theta'}$ replaced $\bm{\theta_e}$ in the archive in the previous iteration. CMA-ME would likely move $\bm{\theta}$ away from $\bm{\theta''}$ as the solution resulted in low improvement. In contrast, CMA-ES would move towards $\bm{\theta''}$ as it ranks only by the objective $f$, ignoring previously discovered solutions with similar measure values.







In the above example, CMA-ME moves away from high performing solutions in order to maximize how the archive changes. However, in domains with hard-to-optimize objective functions, it is beneficial to perform more optimization steps in high-performing regions~\citep{tjanaka2022approximating}. 

Like CMA-ME, CMA-MAE maintains a discount function $f_A(\bm{\theta'})$ and ranks solutions by improvement $f(\bm{\theta'})-f_A(\bm{\theta'})$. However, instead of setting $f_A(\bm{\theta'})$ equal to $f(\bm{\theta_e})$, we set $f_A(\bm{\theta'})=t_e$, where $t_e$ is an acceptance threshold maintained for each cell in the archive $A$.  When we generate a new candidate solution, we control the rate that $t_e$ changes when the archive is updated by the archive learning rate $\alpha$ as follows: $t_e \leftarrow (1-\alpha) t_e + \alpha f(\bm{\theta'})$.




The archive learning rate $\alpha$ in CMA-MAE allows us to control how quickly we leave a high-performing region of measure space. For example, consider discovering solutions in the same cell with objective value 100 in 5 consecutive iterations. The improvement values computed by CMA-ME would be $100,0,0,0,0$, thus CMA-ME would move rapidly away from this cell. The improvement values computed by CMA-MAE with $\alpha=0.5$ would diminish smoothly as follows:  $100,50,25,12.5,6.25$, enabling further exploitation of the high-performing region. 


Next, we walk through the CMA-MAE algorithm step-by-step. Algorithm~\ref{alg:cma-mae} shows the pseudo-code for CMA-MAE with the differences from CMA-ME highlighted in yellow. First, on line~\ref{mae:initialize_t} we initialize the acceptance threshold to $min_f$. In each iteration we sample $\lambda$ solutions around the current search point $\bm{\theta}$ (line~\ref{mae:sample}). For each candidate solution $\bm{\theta_i}$, we evaluate the solution and compute the objective value $f(\bm{\theta_i})$ and measure values $\bm{m}(\bm{\theta_i})$ (line~\ref{mae:evaluate}). Next, we compute the cell $e$ in the archive that corresponds to the measure values and the improvement $\Delta_i$ over the current threshold $t_e$ (lines~\ref{mae:calc_cell}-\ref{mae:calc_imp}). If the objective crosses the acceptance threshold $t_e$, we replace the incumbent $\bm{\theta_e}$ in the archive and increase the acceptance threshold $t_e$ (lines~\ref{mae:start_arc_update}-\ref{mae:inc_thresh}). Next, we rank all candidate solutions $\bm{\theta_i}$ by their improvement $\Delta_i$. Finally, we step our search point $\bm{\theta}$ and adapt our covariance matrix $\Sigma$ towards the direction of largest improvement (lines~\ref{mae:ranking}-\ref{mae:adaptation}) according to CMA-ES's update rules~\citep{hansen:cma16}.




\begin{algorithm}[t!]
\SetAlgoLined
\caption{Covariance Matrix Adaptation MAP-Annealing (CMA-MAE)}
\label{alg:cma-mae}
\SetKwInOut{Input}{input}
\SetKwInOut{Result}{result}
\SetKwProg{CMAMAE}{CMA-MAE}{}{}
\DontPrintSemicolon
\CMAMAE{$(evaluate, \bm{\theta_0}, N, \lambda, \sigma,$ \hl{$min_f, \alpha$})}
{
\Input{An evaluation function $evaluate$ that computes the objective and measures, an initial solution $\bm{\theta_0}$, a desired number of iterations $N$, a branching population size $\lambda$, an initial step size $\sigma$, \hl{ a minimal acceptable solution quality $min_f$, and an archive learning rate $\alpha$.}}
\Result{Generate $N \lambda$ solutions storing elites in an archive $A$.}

\BlankLine
Initialize solution parameters $\bm{\theta}$ to $\bm{\theta_0}$, CMA-ES parameters $\Sigma=\sigma I$ and $\bm{p}$, where we let $\bm{p}$ be the CMA-ES internal parameters. 

\hl{Initialize the archive A and the acceptance threshold $t_e$ with $min_f$ for each cell $e$.} \label{mae:initialize_t}

\For{$iter\leftarrow 1$ \KwTo $N$}{


\For{$i\leftarrow 1$ \KwTo $\lambda$}{ \label{mae:startforloop} 
  $\bm{\theta_i} \sim \mathcal{N}(\bm{\theta},\Sigma)$\; \label{mae:sample} 
  

 $f, \bm{m} \gets \mbox{evaluate}(\bm{\theta_i})$\; \label{mae:evaluate}

 $e \gets \mbox{calculate\_cell}(A, \bm{m})$\; \label{mae:calc_cell}

 \hl{$\Delta_i \gets f - t_e$} \label{mae:calc_imp}

    \If{ \hl{$f > t_e$}}{
    \label{mae:start_arc_update}
    Replace the current occupant in cell $e$ of the archive $A$ with $\bm{\theta_i}$\;
    
    \hl{$t_e \gets (1 - \alpha) t_e + \alpha f$\;} \label{mae:inc_thresh}
 }


}

  \label{mae:endforloop}
\hl{rank $\bm{\theta}_i$ by $\Delta_i$ \;} \label{mae:ranking}

Adapt CMA-ES parameters $\bm{\theta},\Sigma,\bm{p}$ based on improvement ranking $\Delta_i$\label{mae:adaptation}

\If{\mbox{CMA-ES converges}}{
Restart CMA-ES with $\Sigma = \sigma I$.

Set $\bm{\theta}$ to a randomly selected existing cell $\bm{\theta_i}$  from the archive

}
}

}
\end{algorithm}

\noindent\textbf{CMA-MAEGA.} We note that our augmentations to the CMA-ME algorithm only affects how we replace solutions in the archive and how we calculate $\Delta_i$. CMA-ME and CMA-MAEGA replace solutions and calculate $\Delta_i$ identically, so we apply the same augmentations to CMA-MEGA to form a new DQD algorithm, CMA-MAEGA, in Appendix~\ref{sec:cma-maega}.



\section{Theoretical Properties of CMA-MAE}

We provide insights about the behavior of CMA-MAE for different $\alpha$ values. We include all proofs in Appendix~\ref{sec:cmame-theory}. CMA-MAEGA has similar theoretical properties discussed in Appendix~\ref{sec:cmamaega-theory}.


\begin{theorem}
The CMA-ES algorithm is equivalent to CMA-MAE when $\alpha=0$, if CMA-ES restarts from an archive solution.
\end{theorem}


The next theorem states that CMA-ME is equivalent to CMA-MAE when $\alpha = 1$ with the following caveats: First, we assume that CMA-ME restarts only by the CMA-ES restart rules, rather than the additional ``no improvement'' restart rule in prior work~\citep{fontaine2020covariance}. Second, we assume that both CMA-ME and CMA-MAE leverage $\mu$ selection~\citep{hansen:cma16} rather than filtering selection~\citep{fontaine2020covariance}.




\begin{theorem}
\label{thm:CMA-ME}
The CMA-ME algorithm is equivalent to CMA-MAE when $\alpha = 1$ and $min_f$ is an arbitrarily large negative number.
\end{theorem}

We next provide theoretical insights on how the discount function $f_A$ smoothly increases from a constant function $min_f$ to the discount function used by CMA-ME, as $\alpha$ increases from 0 to 1. We focus on the special case of a fixed sequence of candidate solutions.


\begin{theorem}
\label{thm:monot}
Let $\alpha_i$ and $\alpha_j$ be two archive learning rates for archives $A_i$ and $A_j$ such that $0 \leq \alpha_i < \alpha_j \leq 1$. For two runs of CMA-MAE that generate the same sequence of $m$ candidate solutions $\{S\} = \bm{\theta_1}, \bm{\theta_2}, ..., \bm{\theta_m}$, it follows that $f_{A_i}(\bm{\theta}) \leq f_{A_j}(\bm{\theta})$ for all $\bm{\theta}\in \mathbb{R}^n$.
\end{theorem}


Finally, we wish to provide insights about the exploration properties of CMA-MAE for an archive learning rate $\alpha$ between 0 and 1, when the objective $f$ is constant. Consider an approximate density descent algorithm that is identical to CMA-ME, but differs by how solutions are ranked. Specifically, we assume that this  algorithm maintains  an occupancy count $o_e$ for each cell $e$, which is the number of times a solution was generated in that cell. This algorithm ranks solutions based on the occupancy count of the cell that the solution maps to, where solutions that discover less frequently visited cells are ranked higher. 



\begin{theorem}
\label{thm:density_descent}
The CMA-MAE algorithm optimizing a constant objective function $f(\bm{\theta}) = C$ for all $\bm{\theta} \in \mathbb{R}^n$ is equivalent to the approximate density descent algorithm, when $0 < \alpha < 1$ and $min_f < C$.
\end{theorem}

While Theorem~\ref{thm:density_descent} assumes a constant objective $f$, we conjecture that the theorem holds true generally when threshold $t_e$ in each cell $e$ approaches the local optimum within the cell boundaries. 




\section{Experiments}
\label{sec:experiments}

We compare the performance of CMA-MAE with the state-of-the-art QD algorithms MAP-Elites, MAP-Elites (line), and CMA-ME, using existing Pyribs~\citep{pyribs} QD library implementations. We set $\alpha=0.01$ for CMA-MAE and include additional experiments for varying $\alpha$ in section~\ref{sec:robustness}.
Because annealing methods replace solutions based on the threshold, we retain the best solution in each cell for comparison purposes. We include additional comparisons between CMA-MEGA and CMA-MAEGA -- the gradient-based counterpart of CMA-MAE -- in Appendix~\ref{sec:additional-results}.




We select the benchmark domains from prior work~\citep{fontaine2021differentiable}: linear projection~\citep{fontaine2020covariance}, arm repertoire~\citep{cully2017quality}, and latent space illumination~\citep{fontaine2020illuminating}. To evaluate the good exploration properties of CMA-MAE on flat objectives, we introduce a variant of the linear projection domain to include a ``plateau'' objective function that is constant everywhere for solutions within a fixed range and has a quadratic penalty for solutions outside the range. We describe the domains in detail in Appendix~\ref{sec:domain}. 







\subsection{Experiment Design}

\label{subsec:design}
\noindent\textbf{Independent Variables.} We follow a between-groups design with two independent variables: the algorithm and the domain. 

\noindent\textbf{Dependent Variables.} We use the sum of $f$ values of all cells in the archive, defined as the \mbox{QD-score}~\cite{pugh2015confronting}, as a metric for the the quality and diversity of solutions. Following prior work~\citep{fontaine2021differentiable}, we normalize the QD-score metric by the archive size (the total number of cells from the tesselation of measure space) to make the metric invariant to archive resolution. We additionally compute the coverage, defined as the number of occupied cells in the archive divided by the total number of cells. 


\begin{table*}[t]
\centering
\resizebox{1.0\linewidth}{!}{
\begin{tabular}{l|rr|rr|rr|rr|rr}
\hline
             & \multicolumn{2}{l|}{LP (sphere)}     & \multicolumn{2}{l|}{LP (Rastrigin)}   &
             \multicolumn{2}{l|}{LP (plateau)}        & \multicolumn{2}{l|}{Arm Repertoire} &\multicolumn{2}{l}{LSI} \ \\ 
    \toprule
Algorithm     & QD-score & Coverage & QD-score & Coverage & QD-score & Coverage  & QD-score & Coverage & QD-score & Coverage \\
    \midrule
MAP-Elites  &   41.64 & 50.80\% & 31.43 & 47.88\% & 47.07 &   47.07\% & 71.40 & 74.09\% & 12.85 & 19.42\%\\
MAP-Elites (line)  & 49.07 & 60.42\% & 38.29 & 56.51\% & 52.20 & 52.20\% & 74.55  & 75.61\% & 14.40 & 21.11\% \\
CMA-ME    & 36.50 & 42.82\% & 38.02 & 53.09\% & 34.54  & 34.54\% & 75.82 & 75.89\% &  14.00 & 19.57\% \\
 CMA-MAE   & \textbf{64.86} & \textbf{83.31\%} & \textbf{52.65} & \textbf{80.46\%} & \textbf{79.27} & \textbf{79.29\%} & \textbf{79.03} & \textbf{79.24\%} & \textbf{17.67} & \textbf{25.08\%} \\

  \bottomrule
\end{tabular}
 }
\caption{Mean QD-score and coverage values after 10,000 iterations for each algorithm per domain.}
\label{tab:results}
\end{table*}

\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/sphere.pdf}
\includegraphics[width=1.0\columnwidth]{figures/plateau.pdf}
\includegraphics[width=1.0\columnwidth]{figures/arm.pdf}

\caption{QD-score plot with 95$\%$ confidence intervals and heatmaps of generated archives by CMA-MAE and CMA-ME for the linear projection sphere (top), plateau (middle), and arm repertoire (bottom) domains.}
\label{fig:summary}
\end{figure*}


\subsection{Analysis}
Table~\ref{tab:results} shows the QD-score and coverage values for each algorithm and domain, averaged over 20 trials for the linear projection (LP) and arm repertoire domains and over 5 trials for the LSI domain. Fig.~\ref{fig:summary} shows the QD-score values  for increasing number of iterations and example archives for CMA-MAE and CMA-ME, with 95\% confidence intervals.

We conducted a two-way ANOVA to examine the effect of the algorithm and domain (LP (sphere), LP (Rastrigin), LP (plateau), arm repertoire, and LSI) on the QD-score. There was a significant interaction between the search algorithm and the domain ($F(12,320)=1958.34,p<0.001$). Simple main effects analysis with Bonferroni corrections showed that CMA-MAE outperformed all baselines in all benchmark domains. 

We note that for the arm repertoire domain we can compute the optimal archive coverage for a given resolution, by testing whether each cell overlaps with a circle of radius equal to the maximum arm length (see Appendix~\ref{sec:domain}). We observe that CMA-MAE approaches the computed optimal coverage $80.24\%$ for a resolution of $100 \times 100$ and  outperforms CMA-MEGA~\citep{fontaine2021differentiable} (see Appendix~\ref{sec:additional-results}).







These results show that the archive learning rate $\alpha$ is particularly beneficial for CMA-MAE. We observe that CMA-MAE initially explores regions of the measure space that have high-objective values. Once the archive becomes saturated, CMA-MAE reduces to approximate density descent, as we prove in Theorem~\ref{thm:density_descent} for flat objectives. On the other hand, CMA-ME does not receive any exploration signal when the objective landscape becomes flat, resulting in poor performance.
Fig.~\ref{fig:collage_Beyonce} shows an example collage generated by a single run of CMA-MAE for LSI.

\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/collage_beyonce_mae.pdf}
\caption{A latent space illumination collage generated by CMA-MAE for the objective ``A photo of Beyonce.'' and for measures ``A small child.'' and ``A woman with long blonde hair.'' The axis values indicate the score returned by the CLIP model, where lower score indicates a better match.}

\label{fig:collage_Beyonce}
\end{figure}










\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/alpha_results.pdf}
\caption{Final QD-score of each algorithm for 25 different archive resolutions.}
\label{fig:results_resolutions}
\end{figure*}

\section{On the Robustness of CMA-MAE} \label{sec:robustness}



Next, we present two studies that evaluate the robustness of CMA-MAE across two hyperparameters that may affect algorithm performance: the archive learning rate $\alpha$ and the archive resolution. 

\begin{table*}[t]
\centering
\resizebox{1.0\linewidth}{!}{
\begin{tabular}{r|rr|rr|rr|rr}
\hline
             & \multicolumn{2}{l|}{LP (sphere)}     & \multicolumn{2}{l|}{LP (Rastrigin)}   &
             \multicolumn{2}{l|}{LP (plateau)}        & \multicolumn{2}{l}{Arm Repertoire} \\ \toprule
$\alpha$ (CMA-MAE)  & QD-score & Coverage & QD-score & Coverage  & QD-score & Coverage   & QD-score & Coverage \\
    \midrule
0.000  & 5.82 & 6.06\%  & 5.33 & 6.24\% &  19.49 & 19.49\% & 65.91 & 66.25\% \\
0.001  & 62.65 & 79.36\% & 47.87 & 68.10\% & 77.60  &  77.68\% & 78.63 & 79.07\% \\
0.010  & \textbf{64.86} & \textbf{83.31\%} & \textbf{52.65} & \textbf{80.56\%} & 79.27   & 79.29\% & \textbf{79.03} & \textbf{79.24\%} \\
0.100  & 60.42 & 76.19\% & 48.74  & 72.50\% & \textbf{83.21}  & \textbf{83.21\%} & 78.74 & 78.85\% \\
1.000  & 37.01 & 43.50\% &  37.86 & 52.82\% &  34.00 & 34.00\% & 75.94 & 76.01\%\\
  \bottomrule
\end{tabular}
 }
\caption{Mean QD metrics after 10,000 iterations for CMA-MAE at different learning rates.}
\label{tab:lr_results}
\end{table*}


\noindent\textbf{Archive Learning Rate.}
We examine the effect of different archive learning rates on the performance of CMA-MAE in the linear projection and arm repertoire domains. We vary the learning rate from 0 to 1 on an exponential scale, while keeping the resolution constant in each domain. 

Table~\ref{tab:lr_results} shows that running CMA-ME with the different $0<\alpha<1$ results in relatively similar performance, showing that CMA-MAE is fairly robust to $\alpha$ values. On the other hand, if $\alpha=0$ or $\alpha=1$ the performance drops drastically. Setting $\alpha=1$ results in very similar performance with CMA-ME, which supports our insight from Theorem~\ref{thm:CMA-ME}. 

\noindent\textbf{Archive Resolution.} As noted by previous work~\citep{cully2021multi, fontaine2021differentiable}, quality diversity algorithms in the \mbox{MAP-Elites} family sometimes perform differently when run with different archive resolutions. For example, in the linear projection domain presented in previous work~\citep{fontaine2020covariance}, CMA-ME outperformed MAP-Elites and \mbox{MAP-Elites (line)} for archives of resolution $500 \times 500$, while in this paper we observe that it performs worse for resolution $100 \times 100$. In this study, we investigate how CMA-MAE performs at different archive resolutions, compared to our baselines.


First, we note that the optimal archive learning rate $\alpha$ is dependent on the resolution of the archive. Consider as an example a sequence of solution additions to two archives $A_1$ and $A_2$ of resolution $100 \times 100$ and $200 \times 200$, respectively. $A_2$ subdivides each cell in $A_1$ into four cells, thus archive $A_2$'s thresholds $t_e$ should increase at a four times faster rate than $A_1$. To account for this difference, we compute $\alpha_2$ for $A_2$ via a conversion formula $\alpha_2 = 1 - (1 - \alpha_1) ^ r$ (see derivation in Appendix~\ref{sec:conversion}), where $r$ is the ratio of cell counts between archives $A_1$ and $A_2$. We initialize $\alpha_1=0.01$ for $A_1$. In the above example, $\alpha_2 = 1 - (1 - 0.01) ^ 4 = 0.0394$.


Fig.~\ref{fig:results_resolutions} shows the QD-score of CMA-MAE with the resolution-dependent archive learning rate and the baselines for each benchmark domain. CMA-ME performs worse as the resolution decreases because the archive changes quickly at small resolutions, affecting CMA-ME's adaptation mechanism. On the contrary, MAP-Elites and MAP-Elites (line) perform worse as the resolution increases due to having more elites to perturb.  CMA-MAE's performance is invariant to the resolution of the archive.  










\section{Related Work}

\noindent\textbf{Quality Diversity Optimization.}  The predecessor to quality diversity optimization, simply called diversity optimization, originated with the Novelty Search algorithm~\citep{lehman2011abandoning}, which searches for a collection of solutions that are diverse in measure space. Later work introduced the Novelty Search with Local Competition (NSLC)~\citep{lehman2011evolving} and \mbox{MAP-Elites}~\citep{cully:nature15,mouret2015illuminating} algorithms, which combined single-objective optimization with diversity optimization and were the first QD algorithms. Since then, several QD algorithms have been proposed, based on a variety of single-objective optimization methods, such as Bayesian optimization~\citep{kent2020bop}, evolution strategies~\citep{conti2018improving, colas2020scaling, fontaine2020covariance}, differential evolution~\citep{choi:gecco21}, and gradient ascent~\citep{fontaine2021differentiable}. Several works have improved selection mechanisms~\citep{sfikas2021monte, cully2017quality}, archives~\citep{fontaine2019mapping,vassiliades2016using,smith2016rapid}, perturbation operators~\citep{vassiliades:gecco18, nordmoen2018dynamic}, and emitter scheduling~\citep{cully2021multi, paolo2021sparse}. 

\noindent\textbf{QD with Gradient Information.} Several works combine gradient information with quality diversity optimization in ways that do not leverage the objective and measure gradients directly. For example, in model-based quality diversity optimization~\citep{gaier2018data, hagg2020designing, cazenille2019exploring, Keller2020ModelBasedQS, Lim2021DynamicsAwareQF, zhang2021deep, gaier2020blackbox}, prior work~\citep{rakicevic2021policy} trains an autoencoder on the archive of solutions and leverages the Jacobian of the decoder network to compute the covariance of the Gaussian perturbation. In quality diversity reinforcement learning (QD-RL), several works~\citep{parkerholder2020dvd, pierrot2020diversity, nilsson2021policy, tjanaka2022approximating} approximate a reward gradient or diversity gradient via a critic network, action space noise, or evolution strategies and incorporate those gradients into a QD-RL algorithm.

\noindent\textbf{Acceptance Thresholds.} Our proposed archive learning rate $\alpha$ was loosely inspired by simulated annealing methods~\citep{bertsimas1993simulated} that maintain an acceptance threshold that gradually becomes more selective as the algorithm progresses. The notion of an acceptance threshold is also closely related to minimal criterion methods in evolutionary computation~\citep{lehman2010revising, brant2017minimal, brant2020diversity, stanley2016strictness}. Our work differs by both 1) maintaining an acceptance threshold per archive cell rather than a global threshold and 2) annealing the threshold.

\section{Societal Impacts} 
\label{sec:societal_impacts}
By controlling the trade-off between exploration and exploitation in quality diversity algorithms, we aim towards improving their performance and robustness, thus making these algorithms easier to apply in a wide range of domains and applications. 


One promising application is synthetically extracting datasets from generative 
models to train machine learning algorithms~\cite{jahanian2021generative,besnier2020dataset}. This can raise ethical considerations because generative models can reproduce and exacerbate existing biases in the datasets that they were trained on
~\citep{jain2020imperfect,menon2020pulse}. On the other hand, quality diversity algorithms with carefully selected measure functions can target diversity with desired attributes, thus we hypothesize that they can be effective in generating balanced datasets. 

By attempting to find diverse solutions, QD algorithms are a step towards open-endedness in AI~\cite{stanley2017open} and will often result in unexpected and often surprising emergent behaviors~\citep{lehman2020surprising}. We recognize that this presents several challenges in predictability and monitoring of AI systems~\citep{hendrycks2021unsolved}, and we highlight the importance of future work on balancing the tradeoff between open-endedness and control~\citep{ecoffet2020open}.


\section{Limitations and Future Work}  \label{sec:limitations}

Quality diversity optimization is a rapidly growing branch of stochastic optimization with applications in generative design~\citep{hagg:gecco21, gaier2020blackbox, gaier2018data}, automatic scenario generation in robotics~\citep{fontaine:thri21, fontaine2021importance, fontaine:sa_rss2021}, reinforcement learning~\citep{parkerholder2020dvd, pierrot2020diversity, nilsson2021policy, tjanaka2022approximating}, damage recovery in robotics~\citep{cully:nature15}, and procedural content generation~\citep{gravina2019procedural, fontaine2020illuminating, zhang2021deep, earle2021illuminating, khalifa2018talakat, steckel2021illuminating, schrum2020cppn2gan, sarkar2021generating}. Our paper introduces a new quality diversity algorithm, CMA-MAE, that bridges the gap between single-objective optimization and quality diversity optimization.

Our approach introduced two hyperparameters, $\alpha$ and $min_f$, to allow a smooth blending between CMA-ES's and CMA-ME's algorithmic behavior. We observed that an $\alpha$
 set strictly between $0$ and $1$ yields theoretical exploration improvements and that CMA-MAE is robust with respect to the exact choice of $\alpha$. We additionally derived a conversion formula that converts an $\alpha_1$ for a specific archive resolution to an equivalent $\alpha_2$ for a different resolution. However, the conversion formula still requires practitioners to specify a good initial value of $\alpha_1$. Future work will explore ways to automatically initialize $\alpha$, similar to how CMA-ES automatically assigns internal parameters~\citep{hansen:cma16}.



In addition to CMA-MAE, these improvements can be applied to the DQD algorithm CMA-MEGA to form the CMA-MAEGA algorithm (see Appendix~\ref{sec:cma-maega}). While DQD greatly accelerates QD optimization when exact gradient information is available, this is not the case  in all application settings. Moreover, more DQD application domains are needed to accurately evaluate the benefits of annealing $\alpha$ in CMA-MAEGA, since currently the only DQD application domains are the linear projection, arm repertoire, and latent space illumination domains. Given the prevalence of gradients in deep learning~\citep{goodfellow2016deep}, we anticipate more application domains for DQD in the near future.




{
\small
\bibliographystyle{plainnat}

\section{Hyperparameter Selection}
\label{sec:hyperparams}
For all domains we mirror the hyperparameter selection in prior work~\citep{fontaine2021differentiable}. For CMA-MAE and \mbox{CMA-MAEGA}, we duplicate the hyperparameter selections of CMA-ME and CMA-MEGA, respectively. Following prior work~\citep{fontaine2020covariance}, we run all algorithms with 15 emitters on the linear projection and arm repertoire domains. In the latent space illumination domain, we run experiments with only one emitter, due to the computational expense of the domain. Emitters are independent CMA-ES instances that run in parallel with a shared archive. For each algorithm, we select a batch size $\lambda = 36$ following prior work~\citep{fontaine2021differentiable}. For MAP-Elites and MAP-Elites (line), we initialize the archive with $100$ random solutions, sampled from the distribution $\mathcal{N}(\bm{0},I)$. These initial solutions do not count in the evaluation budget for MAP-Elites and MAP-Elites (line). For algorithms in the CMA-ME family (CMA-ME, CMA-MAE, CMA-MEGA, and CMA-MEGA), we initialize $\bm{\theta_0} = \bm{0}$ for every domain.

In our experiments we want to directly compare the ranking mechanisms of \mbox{CMA-ME} and \mbox{CMA-MAE}. However, \mbox{CMA-ME} is typically run with a ``no improvement'' restart rule, where the algorithm will restart if no solution changes the archive. Due to \mbox{CMA-MAE}'s annealed acceptance threshold $t_e$, a ``no improvement'' restart rule would cause \mbox{CMA-ME} and \mbox{CMA-MAE} to restart at different rates, confounding the effects of restarts and rankings. Filter selection also has a similar confounding effect as solutions are selected if they change the archive. For these reasons, in the main paper we run CMA-ME with a basic restart rule (CMA-ES style restarts only~\citep{hansen:cma16}) and $\mu$ selection~\citep{hansen:cma16} (selecting the top half of the ranking). In Appendix Section~\ref{sec:additional-results}, we run an extra CMA-ME with filter selection and the ``no improvement'' restart rule, which we denote CMA-ME*.

In the latent space illumination domain, due to the computational expense of the domain, we compare directly against the results from prior work~\citep{fontaine2021differentiable}, where we obtained the data (MIT license) with consent from the authors. For CMA-MAE and CMA-MAEGA we include the ``no improvement'' restart rule to match CMA-ME and CMA-MEGA as closely as possible. For this domain, we take gradient steps with the Adam optimizer~\citep{kingma2014adam}, following the recommendation of prior work~\citep{fontaine2021differentiable}. However, we run CMA-MAE with $\mu$ selection, since we found that  small values of the archive learning rate $\alpha$ makes filter selection worse. %



\noindent\textbf{Linear Projection (sphere, Rastrigin, plateau).} 
\begin{itemize}
    \item MAP-Elites: $\sigma = 0.5$
    \item MAP-Elites (line): $\sigma_1 = 0.5$,  $\sigma_2 = 0.2$
    \item CMA-ME: $\sigma = 0.5$, $\mu$ selection, basic restart rule
    \item CMA-ME*: $\sigma = 0.5$, filter selection, no improvement restart rule
    \item CMA-MAE: $\sigma = 0.5$, $\alpha = 0.01$, $min_f = 0$, $\mu$ selection, basic restart rule
    \item CMA-MEGA:  $\sigma_g = 10.0$, $\eta = 1.0$, basic restart rule, gradient ascent optimizer
    \item CMA-MAEGA:  $\sigma_g = 10.0$, $\eta = 1.0$,  $\alpha = 0.01$, $min_f = 0$, basic restart rule, gradient ascent optimizer
\end{itemize}

\noindent\textbf{Arm Repertoire.}
\begin{itemize}
    \item MAP-Elites:  $\sigma = 0.1$
    \item MAP-Elites (line):  $\sigma_1 = 0.1$, $\sigma_2 = 0.2$
    \item CMA-ME: $\sigma = 0.2$, $\mu$ selection, basic restart rule
    \item CMA-ME*: $\sigma = 0.2$, filter selection, no improvement restart rule
    \item CMA-MAE: $\sigma = 0.2$, $\alpha = 0.01$, $min_f = 0$, $\mu$ selection, basic restart rule
    \item CMA-MEGA: $\sigma_g= 0.05$, $\eta = 1.0$, basic restart rule, gradient ascent optimizer
    \item CMA-MAEGA: $\sigma_g= 0.05$, $\eta = 1.0$,  $\alpha = 0.01$, $\min_f = 0$, basic restart rule, gradient ascent optimizer
\end{itemize}

\newpage

\noindent\textbf{Latent Space Illumination.}
\begin{itemize}
    \item MAP-Elites:  $\sigma = 0.2$
    \item MAP-Elites (line):  $\sigma_1 = 0.1$, $\sigma_2 = 0.2$ 
    \item CMA-ME: $\sigma = 0.02$, filter selection, no improvement restart rule
    \item CMA-MAE: $\sigma = 0.02$, $\alpha = 0.1$, $min_f = 55$, $\mu$ selection, no improvement restart rule
    \item CMA-MEGA:  $\sigma_g = 0.002$, $\eta = 0.002$, Adam optimizer, no improvement restart rule
    \item CMA-MAEGA:  $\sigma_g = 0.002$, $\eta = 0.002$, $\alpha = 0.1$, $min_f = 55$, Adam optimizer, no improvement restart rule
\end{itemize}

\noindent\textbf{Adam Hyperparameters.}
We use the same hyperparameters as previous work~\cite{styleclip, fontaine2021differentiable}.
\begin{itemize}
    \item $\beta_1=0.9$
    \item $\beta_2=0.999$
\end{itemize}

\noindent\textbf{Archives.} For the linear projection and arm repertoire domains, we initialize an archive of $100 \times 100$ cells for all algorithms. For latent space illumination we initialize an archive of $200 \times 200$ cells for all algorithms, following previous work~\citep{fontaine2021differentiable}.


\section{Domain Details}
\label{sec:domain}

To experimentally evaluate both CMA-MAE and CMA-MAEGA, we select domains from prior work~\citep{fontaine2021differentiable}: linear projection~\citep{fontaine2020covariance}, arm repertoire~\citep{cully2017quality}, and latent space illumination~\citep{fontaine2020illuminating}. While many quality diversity optimization domains exist, we select these because gradients of $f$ and $\bm{m}$ are easy to compute analytically and allow us to evaluate DQD algorithms in addition to derivative-free QD algorithms. To evaluate the good exploration properties of CMA-MAE on flat objectives, we introduce a variant of the linear projection domain to include a ``plateau'' objective function. 


\noindent\textbf{Linear Projection.} The linear projection domain~\citep{fontaine2020covariance} was introduced to benchmark distortions caused by mapping a high-dimensional search space to a low-dimensional measure space. The domain includes a 2D measure space formed by a linear projection that bounds the contribution of each component $\theta_i$ of the projection to the range $[-5.12, 5.12]$. QD algorithms must adapt the step size of each component $\theta_i$ to slowly approach the extremes of the measure space, with a harsh penalty for components outside $[-5.12, 5.12]$. As QD domains must also provide an objective, the linear projection domain included two objectives from the black-box optimization benchmarks~\citep{hansen:arxiv16,hansen:gecco10}: sphere and Rastrigin. Following prior work~\citep{fontaine2020covariance}, we run all experiments for $n=100$.

\begin{figure}[t!]
    \centering
    \begin{tabular}{cc}
    \begin{subfigure}{0.45\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/bates.pdf}
      \caption{}
      \end{subfigure}& 
    \begin{subfigure}{0.45\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/clip.pdf}
      \caption{}
      \end{subfigure}
     \end{tabular}
   \caption{Measure function figures reproduced from prior work~\citep{fontaine2020covariance, fontaine2021differentiable} with the authors' permission. (a) a Bates distribution. (b) $clip$ function for defining the measures in the linear projection domain.}
    \label{fig:lp_examples}
\end{figure}


Formally, the measure functions are defined as a linear projection, a weighted sum of the components $\theta_i \in \mathbb{R}$ of a solution $\bm{\theta} \in \mathbb{R} ^ n$. The first measure function $m_1$ is a weighted sum of the first half of the solution $\bm{\theta}$, and the second measure function $m_2$ is a weighted sum of the second half of the solution $\bm{\theta}$ (see Eq.~\ref{eq:projection}). To ensure that all solutions mapped to measure space occupy a finite volume, the contribution in measure space of each component $\theta_i$ is bounded to the range $[-5.12, 5.12]$ via a clip function (see Eq.~\ref{eq:clip}) that applies a harsh penalty for solution components $\theta_i$ stepping outside the range $[-5.12, 5.12]$.

\begin{equation} \label{eq:clip}
  clip(\theta_i) =
  \begin{cases}
    \theta_i & \text{if $-5.12 \leq \theta_i \leq 5.12$} \\
    5.12 / \theta_i & \text{otherwise}
  \end{cases}
\end{equation}

\begin{equation} \label{eq:projection}
m(\bm{\theta}) = \left( \sum_{i=1}^{\lfloor{\frac{n}{2}}\rfloor}  {clip(\theta_i)}, \sum_{i=\lfloor{\frac{n}{2}}\rfloor+1}^{n}  {clip(\theta_i)} \right)
\end{equation}


Fig.~\ref{fig:lp_examples} visualizes \emph{why} the linear projection domain is challenging. First, we note that the density of solutions in search space mapped to measure space mostly occupies the region close to $\bm{0}$. To justify why, consider sampling uniformly in the hypercube $[-5.12, 5.12]^n$ in search space. We note that each of these points maps to the linear region of the measure functions and each of our measures becomes a sum of random variables. If we divide by $n$, we normalize by the dimensions of the search space, then the measure functions become an \emph{average} of random variables. The average of $n$ uniform random variables is the Bates distribution~\citep{johnson1995continuous}, a distribution that narrows in variance as $n$ grows larger. Without the $clip$ function, a QD algorithm could simply increase a single $\bm{\theta_i}$ to reach any point in the measure space. However, the clip function prevents this by bounding the contribution of each component of $\bm{\theta}$ to the range $[-5.12, 5.12]$. To reach the extremes of measure space all components $\theta_i$ must converge to the extremums $\pm 5.12$. The linear projection domain is challenging to explore due to both the clustering of solutions in a small region of measure space and the heavy measure space penalties applied by the clip function when a component $\theta_i$ leaves the region $[-5.12, 5.12]$.

\begin{figure}[t!]
    \centering
    \begin{tabular}{ccc}
    \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/sphere_func.pdf}
      \caption{}
      \end{subfigure}& 
    \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/rastrigin_func.pdf}
      \caption{}
      \end{subfigure}&
      \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/plateau_func.pdf}
      \caption{}
      \end{subfigure}
     \end{tabular}
   \caption{Objective functions for the linear projection domain in their minimization form for $n=1$. (a) a sphere function. (b) the Rastrigin function. (c) a plateau function.}
    \label{fig:obj_examples}
\end{figure}

Next, we describe the linear projection domain's objective functions visualized in Fig.~\ref{fig:obj_examples}.

The objectives of the linear projection domain satisfy the requirements that a QD domain needs to have an objective and are of lesser importance than the measure function definitions, since the benchmark primarily evaluates exploration capabilities. Previous work~\citep{fontaine2020covariance} introduces two objectives from the black-box optimization benchmarks ~\citep{hansen:arxiv16,hansen:gecco10}: sphere and Rastrigin. The sphere function (Eq.~\ref{eq:sphere}) is a quadratic function\footnote{In derivative-free optimization many of the benchmark functions are named after the shape of the contour lines. In the case of quadratic functions with an identity Hessian matrix, the contour lines form hyperspheres.}, while the Rastrigin function (Eq.~\ref{eq:rastrigin}) is a multi-modal function that when smoothed is quadratic. The domain shifts the global optimum to the position $\theta_i = 5.12 \cdot 0.4 = 2.048$.

\begin{equation} \label{eq:sphere}
f_{sphere}(\bm{\theta}) = \sum_{i=1}^{n}  {\theta_i ^ 2}
\end{equation}

\begin{equation} \label{eq:rastrigin}
f_{Rastrigin}(\bm{\theta}) =  10n + \sum_{i=1}^{n} [{\theta_i ^ 2} - 10 \cos(2\pi \theta_i ^ 2)]
\end{equation}

We introduce an additional objective to evaluate the good exploration properties of CMA-MAE on flat objectives. Our ``plateau'' objective function (Eq.~\ref{eq:plateau}) is constant everywhere, but with a quadratic penalty for each component outside the range $[-5.12, 5.12]$. The penalty acts as a regularizer to encourage algorithms to search in the linear region of measure space.

\begin{equation} \label{eq:plateau_single}
  f_{plateau}(\theta_i) =
  \begin{cases}
    0 & \text{if $-5.12 \leq \theta_i \leq 5.12$} \\
    (\abs{x} - 5.12) ^ 2 & \text{otherwise}
  \end{cases}
\end{equation}





\begin{equation} \label{eq:plateau}
f_{plateau}(\bm{\theta}) = \frac{1}{n} \sum_{i=1}^{n}  {f_{plateau}(\theta_i)}
\end{equation}

\begin{figure}[t!]
    \centering
    \begin{tabular}{ccc}
    \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/arm_opt_grasp.pdf}
      \caption{}
      \end{subfigure}& 
    \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/arm_sub_grasp.pdf}
      \caption{}
      \end{subfigure}&
      \begin{subfigure}{0.3\linewidth}
      \includegraphics[width = 1.0\linewidth]{figures/arm_fan.pdf}
      \caption{}
      \end{subfigure}
     \end{tabular}
   \caption{Examples of the Arm Repertoire domain for $n=5$. The figures are reproduced from previous work~\citep{fontaine2021differentiable} with the authors' permission.  (a) An optimal grasp with 0 variance between joint angles. (b) A sub-optimal grasp. (c) An ensemble of 0 variance optimal grasps.}
    \label{fig:arm_examples}
\end{figure}


\noindent\textbf{Arm Repertoire.} The arm repertoire domain~\cite{cully2017quality,vassiliades:gecco18} tasks QD algorithms to find a diverse collection of arm positions for an $n$-dimensional planar robotic arm with revolute joints. The measures in this domain are the 2D coordinates of the robot's end-effector and the objective is to minimize the variance of the joint angles.

In Fig.~\ref{fig:arm_examples}, we visualize example arms for $n=5$ (5-DOF). The optimal solutions in this domain have 0 variance between all joint angles. The measure functions are bounded to the range $[-n, n]$ as each arm segment has a unit length. The reachable cells form a circle of radius $n$. Therefore, the optimal archive coverage is approximately $\frac{\pi n^2}{4n^2} \approx 78.5\%$. An archive can achieve an upper-bound of this ratio that becomes tighter at higher resolutions.  We select $n=100$ (100-DOF) arms for the experiments.

\noindent\textbf{Latent Space Illumination.} Prior work introduced the latent space illumination problem~\citep{fontaine2020illuminating}, the problem of searching the latent space of a generative model with a quality diversity algorithm. We evaluate on the StyleGAN+CLIP version of this problem~\citep{fontaine2021differentiable}, by searching the latent space of StyleGAN~\citep{karras2019style} with a QD algorithm. We form the differentiable objective and measures in this domain by specifying text prompts to the CLIP model~\citep{radford2021learning} that can determine the similarity of an image and text. We specify an objective prompt of ``A photo of Beyonce''. For measures, we would like to have CLIP quantify abstract concepts like the hair length or age of the person in the photo. However, CLIP can only determine similarity of an image and a text prompt. As surrogates for age and hair length, we specify the measure prompts of ``A small child'' and ``A woman with long blonde hair''. The objective and measure functions guide the QD algorithms towards discovering a collection of photos of Beyonc\'e with varying age and hair length.

\noindent\textbf{Transformations of the Objective Function.} We highlight two issues that must be addressed by transforming the objective in each domain. First, we note that the problem definition in each of our domains contains an objective $f$ that must be minimized. In contrast, the QD problem definition specifies an objective $f$ that must be maximized. Second, the QD-score metric, which measures the performance of QD algorithms, requires a non-negative objective function. Following prior work~\citep{fontaine2020covariance, fontaine2021differentiable}, we transform the objective $f$ via a linear transformation: $f' = af+b$. The linear transformation maps function outputs to the range $[0, 100]$.

In the linear projection domain, we estimate the largest objective value for the sphere and Rastrigin function within the region $[-5.12, 5.12]$ for each solution component $\theta_i$. We compute $f(-5.12, -5.12, ..., -5.12)$ for each objective as the maximum. The minimum of each function is $0$. We calculate the linear transformation as:

\begin{equation}
f'(\bm{\theta}) = 100 \cdot \frac{f(\bm{\theta}) - f_{max}}{f_{min} - f_{max}}
\label{eq:obj_transform}
\end{equation}

For our new plateau objective, all solution points within the region $[-5.12, 5.12]^n$ have objective value of $0$. For this objective we set $f_{min} = 0$ and $f_{max} = 100$ and apply the transformation in Eq.~\ref{eq:obj_transform}.

For the arm domain we select $f_{min}=0$ and $f_{max}=1$, and in the LSI domain we select $f_{min} = 0$ and $f_{max} = 10$. We select these values to match prior work~\citep{fontaine2021differentiable}.

\section{Implementation}
\label{sec:implementation}

We replicate the implementation details of prior work~\citep{fontaine2021differentiable}.

\noindent\textbf{Archives.} For the linear projection and arm repertoire domains, we initialize an archive of $100 \times 100$ cells for all algorithms. For latent space illumination we initialize an archive of $200 \times 200$ cells for all algorithms, following previous work~\citep{fontaine2021differentiable}.

\noindent\textbf{Metrics.}  We use the sum of $f$ values of all cells in the archive, defined as the \mbox{QD-score}~\cite{pugh2015confronting}, as a metric for the the quality and diversity of solutions.  Following previous work~\citep{fontaine2021differentiable}, we normalize the QD-score by the total number of cells, both occupied and unoccupied, to make QD-score invariant to the resolution of the archive. We additionally compute the coverage, defined as the number of occupied cells in the archive divided by the total number of cells. 

\noindent\textbf{Computational Resources.} We ran all trials of the linear projection and arm repertoire domains on an AMD Ryzen Threadripper 32-core (64 threads) processor. A run of 20 trials in parallel takes about 20 minutes for the linear projection domain and 25 minutes for the arm repertoire domain. For the latent space illumination domain, we accelerate the StyleGAN+CLIP pipeline on a GeForce RTX 3090 Nvidia GPU. One trial for latent space illumination takes approximately 2 hours and 30 minutes. In all domains, runtime increases when an algorithm obtains better coverage, because we iterate over the archive when QD statistics are calculated.

\noindent\textbf{Software Implementation.} We use the open source Pyribs~\citep{pyribs} library for all algorithms. We implemented the CMA-MAE and CMA-MAEGA algorithms using the same library.


\section{Covariance Matrix Adaptation MAP-Elites via a Gradient Arborescence (CMA-MAEGA)}
\label{sec:cma-maega}

In this section, we provide information of the CMA-MEGA differentiable quality diversity (DQD) algorithm, and we derive CMA-MAE's DQD counterpart: CMA-MAEGA.

\noindent\textbf{CMA-MEGA.} Covariance Matrix Adaptation MAP-Elites via Gradient Arborescence (CMA-MEGA) solves the DQD problem, where the objective $f$ and measures $\bm{m}$ are first-order differentiable. Like \mbox{CMA-ME}, the algorithm maintains a solution point $\bm{\theta} \in \mathbb{R}^n$ and a MAP-Elites archive. \mbox{CMA-MEGA} samples new solutions by perturbing the search point $\bm{\theta}$ via the objective and measure gradients. However, the contribution of each gradient is balanced by gradient coefficients $\bm{c}$: \mbox{$\bm{\theta_i} = \bm{\theta} + c_0 \bm{\nabla} f(\bm{\theta}) + \sum_{j=1}^{k} c_j \bm{\nabla} m_j(\bm{\theta})$}. These coefficients are sampled from a multivariate Gaussian distribution $N(\bm{\mu},\Sigma)$ maintained by the algorithm. After sampling new candidate solutions $\bm{\theta_i}$, the solutions are ranked via the improvement ranking from \mbox{CMA-ME}. \mbox{CMA-MEGA} updates $N(\bm{\mu},\Sigma)$ via the CMA-ES update rules and the algorithm steps $\bm{\theta}$ also in the direction of largest archive improvement. The authors showed that CMA-MEGA approximates a natural gradient step of the QD objective (Eq.~\ref{eq:objective}), but with respect to the gradient coefficients.


\begin{algorithm}[t!]
\SetAlgoLined
\setcounter{AlgoLine}{0}
\caption{Covariance Matrix Adaptation MAP-Annealing via a Gradient Arborescence (\mbox{CMA-MAEGA})}
\label{alg:cma-maega}
\SetKwInOut{Input}{input}
\SetKwInOut{Result}{result}
\SetKwProg{CMAMAEGA}{CMA-MAEGA}{}{}
\DontPrintSemicolon
\CMAMAEGA{$(evaluate, \bm{\theta_0}, N, \lambda, \eta, \sigma_g,$ \hl{$min_f,\alpha$})}
{
\Input{An evaluation function $evaluate$ that computes the objective, the measures, and the gradients of the objective and measures, an initial solution $\bm{\theta_0}$, a desired number of iterations $N$, a branching population size $\lambda$, a learning rate $\eta$, an initial step size for CMA-ES $\sigma_g$, \hl{a minimal acceptable solution quality $min_f$, and an archive learning rate $\alpha$}.}
\Result{Generate $N (\lambda + 1)$ solutions storing elites in an archive $A$.}

\BlankLine
Initialize solution parameters $\bm{\theta}$ to $\bm{\theta_0}$, CMA-ES parameters $\bm{\mu}=\bm{0}$, $\Sigma=\sigma_g I$, and $\bm{p}$, where we let $\bm{p}$ be the CMA-ES internal parameters.

\hl{Initialize the archive $A$ and the acceptance threshold $t_e$ with $min_f$ for each cell $e$.}

\For{$iter\leftarrow 1$ \KwTo $N$}{
$f, \bm{\nabla}_{f}, \bm{m}, \bm{\nabla_m} \gets  \mbox{evaluate}(\bm{\theta})$\; \label{maega:initeval}

$\bm{\nabla}_{f} \gets \mbox{normalize}(\bm{\nabla}_{f}),\bm{\nabla}_{m} \gets \mbox{normalize}(\bm{\nabla}_{m})$\;\label{maega:norm}


\If{\hl{$f > t_e$}}{
    Replace the current elite in cell $e$ of the archive $A$ with $\bm{\theta_i}$\;
    
    \hl{$t_e \gets (1 - \alpha) t_e + \alpha f$}\; 
 }

\For{$i\leftarrow 1$ \KwTo $\lambda$}{ \label{maega:startforloop} 
  $\bm{c} \sim \mathcal{N}(\bm{\mu},\Sigma)$\;  

  
  $\bm{\nabla}_i \gets c_0 \bm{\nabla}_f + \sum_{j=1}^{k} c_j \bm{\nabla}_{m_j}$\;\label{lst:line2}

  $\bm{\theta_i}' \gets \bm{\theta} + \bm{\nabla}_i$\;

 $f', *, \bm{m}',* \gets \mbox{evaluate}(\bm{\theta_i}')$\;

 \hl{$\Delta_i \gets f' - t_e$} \label{maega:calc_imp}

 \If{\hl{$f' > t_e$}}{
    \label{maega:start_arc_update}
    Replace the current occupant in cell $e$ of the archive $A$ with $\bm{\theta_i}$\;
    
    \hl{$t_e \gets (1 - \alpha) t_e + \alpha f'$}\; \label{maega:inc_thresh}
 }  %

}  \label{maega:endforloop}
\hl{rank $\bm{\nabla}_i$ by $\Delta_i$ \;} \label{maega:ranking}

$\bm{\nabla}_{\textrm{step}} 
\gets \sum_{i=1}^{\lambda}w_i  \bm{\nabla}_{\textrm{rank[i]}}$\; \label{maega:gradient_linear_comb}

$\bm{\theta} \gets \bm{\theta} +  \eta \bm{\nabla}_{\textrm{step}}$\; \label{maega:gradient_step}

Adapt CMA-ES parameters $\bm{\mu},\Sigma,\bm{p}$ based on improvement ranking $\Delta_i$\label{maega:adaptation}

\If{\mbox{there is no change in the archive}}{
Restart CMA-ES with $\bm{\mu}=0, \Sigma = \sigma_g I$.

Set $\bm{\theta}$ to a randomly selected existing cell $\bm{\theta_i}$  from the archive

}
}

}

\end{algorithm}


\noindent\textbf{CMA-MAEGA.} We note that our augmentations to the CMA-ME algorithm only affects how we replace solutions in the archive and how we calculate $\Delta_i$. Both CMA-ME and CMA-MAEGA replace solutions and calculate $\Delta_i$ identically, so we apply the same augmentations from CMA-ME to CMA-MEGA to form a new DQD algorithm, CMA-MAEGA. Algorithm~\ref{alg:cma-maega} shows the pseudo-code for CMA-MAEGA with the differences from CMA-MEGA highlighted in yellow.

\begin{table*}[t]
\centering
\resizebox{1.0\linewidth}{!}{
\begin{tabular}{l|rr|rr|rr|rr|rr}
\hline
             & \multicolumn{2}{l|}{LP (sphere)}     & \multicolumn{2}{l|}{LP (Rastrigin)}   &
             \multicolumn{2}{l|}{LP (plateau)}        & \multicolumn{2}{l|}{Arm Repertoire} &\multicolumn{2}{l}{LSI} \ \\ 
    \toprule
Algorithm     & QD-score & Coverage & QD-score & Coverage & QD-score & Coverage  & QD-score & Coverage & QD-score & Coverage \\
    \midrule

 CMA-MEGA  & 75.32 & \textbf{100.00\%} & \textbf{63.07} & \textbf{100.00\%} & \textbf{100.00}  &\textbf{100.00\%} & 75.21 & 75.25\% & 16.08 & 22.58\% \\
 CMA-MAEGA  & \textbf{75.39} & \textbf{100.00\%} & 63.06 & \textbf{100.00\%} & \textbf{100.00} & \textbf{100.00\%} & \textbf{79.27} & \textbf{79.35\%} & \textbf{16.20} & \textbf{23.83\%}\\

  \bottomrule
\end{tabular}
 }
\caption{Mean QD-score and coverage values after 10,000 iterations for each algorithm per domain.}
\label{tab:dqd_results}
\end{table*}



\noindent\textbf{Experiments.} We compare CMA-MEGA and CMA-MAEGA in the five benchmark domains.  Table~\ref{tab:dqd_results} shows the QD-score and coverage values for each algorithm and domain, averaged over 20 trials for the linear projection (LP) and arm repertoire domains and over 5 trials for the LSI domain. We conducted a two-way ANOVA to examine the effect of the algorithm and domain (LP (sphere), LP  (Rastrigin), LP (plateau), arm repertoire and LSI) on the QD-score. There was a significant interaction between the search algorithm and the domain ($F(4,160)=403.30,p<0.001$). Simple main effects analysis with Bonferroni corrections showed that CMA-MAEGA outperformed CMA-MEGA in the LP (sphere) and arm repertoire domains. There was no statistically significance difference between the two algorithms in the LP (Rastrigin), plateau, and LSI domains.  We attribute this to the large variance of the runs and to the fact that exact gradient information enables both algorithms to efficiently explore the search space. Overall, these preliminary results show that CMA-MAEGA can outperform CMA-MEGA in some domains and warrant future research on a larger number of DQD application domains.




\section{Theoretical Properties of CMA-MAE}
\label{sec:cmame-theory}

\begin{theorem}
\label{thm_apx:alpha0}
The CMA-ES algorithm is equivalent to CMA-MAE when $\alpha=0$, if CMA-ES restarts from an archive solution.
\end{theorem}

\begin{proof}
CMA-ES and CMA-MAE differ only on how they rank solutions.
CMA-ES ranks solutions purely based on the objective $f$, while CMA-MAE ranks solutions by $f-t_e$, where $t_e$ is the acceptance threshold initialized by $min_f$. Thus, to show that CMA-ES is equivalent to CMA-MAE for $\alpha=0$, we only need to show that they result in identical rankings. 

In CMA-MAE, $t_e$ is updated as follows: $t_e \gets (1-\alpha) t_e + \alpha f$. For $\alpha=0$, $t_e = min_f$ is invariant for the whole algorithm: $t_e \leftarrow 1 t_e + 0 f = t_e$. Therefore, CMA-MAE ranks solutions based on $f-min_f$. However, comparison-based sorting is invariant to order-preserving transformations of the values being sorted~\cite{hansen:cma16}. Thus, CMA-ES and CMA-MAE rank solutions identically.
\end{proof}

Next, we prove that CMA-ME is equivalent to CMA-MAE with the following caveats. First, we assume that CMA-ME restarts only with the CMA-ES restart rules, rather than the additional ``no improvement'' restart condition in prior work~\citep{fontaine2020covariance}. Second, we assume that both CMA-ME and CMA-MAE leverage $\mu$ selection rather than filtering selection.

\begin{lemma}
During execution of the CMA-MAE algorithm with $\alpha=1$, the threshold $t_e$ is equal to $f(\bm{\theta_e})$ for cells that are occupied by a solution $\bm{\theta_e}$ and to $min_f$ for all empty cells. 
\label{Le_1}
\end{lemma}

\begin{proof}
We will prove the lemma by induction. All empty cells are initialized with $t_e=min_f$, satisfying the basis step. Then, we will show that if the statement holds after $k$ archive updates, it will hold after a subsequent update $k+1$. 

Assume that at step $k$ we generate a new solution $\bm{\theta_i}$ mapped to a cell $e$.
We consider two cases: 

\noindent\textbf{Case 1}: The archive cell $e$ is empty. Then, $f(\bm{\theta_i})>min_f$ and both CMA-ME and CMA-MAE will place $\bm{\theta_i}$ in the archive as the new cell occupant $\bm{\theta_e}$.  The threshold $t_e$ is updated as $t_e = (1-\alpha) t_e + \alpha f(\bm{\theta_e}) = 0 min_f + 1 f(\bm{\theta_e}) = f(\bm{\theta_e})$. 

\noindent\textbf{Case 2}: The archive cell $e$ contains an incumbent solution $\bm{\theta_e}$. Then, either $f(\bm{\theta_i})\leq f(\bm{\theta_e})$ or $f(\bm{\theta_i})>f(\bm{\theta_e})$. If $f(\bm{\theta_i})\leq f(\bm{\theta_e})$, then the archive does not change and the inductive step holds via the inductive hypothesis. If $f(\bm{\theta_i})>f(\bm{\theta_e})$, then $\bm{\theta_i}$ becomes the new cell occupant $\bm{\theta_e}$ and $t_e$ is updated as $t_e = (1-\alpha) t_e + \alpha f(\bm{\theta_e}) = 0 t_e + 1 f(\bm{\theta_e}) = f(\bm{\theta_e})$. 
\end{proof}

\begin{theorem}
\label{thm_apx:alpha1}
The CMA-ME algorithm is equivalent to CMA-MAE when $\alpha = 1$ and $min_f$ is an arbitrarily large negative number.
\end{theorem}

\begin{proof}
Both CMA-ME and CMA-MAE rank candidate solutions $\bm{\theta_i}$ based on improvement values  $\Delta_i$. While CMA-ME and CMA-MAE compute $\Delta_i$ differently, we will show that for $\alpha=1$, the rankings are identical for the two algorithms.

We assume a new candidate solution  mapped to a cell $e$.  We describe first the computation of $\Delta_i$ for CMA-ME. CMA-ME ranks solutions that discover an empty cell based on their objective value. Thus, if $\bm{\theta_i}$ discovers an empty cell, $\Delta_i = f(\bm{\theta_i})$. On the other hand, if $\bm{\theta_i}$ is mapped to a cell occupied by another solution $\bm{\theta_e}$, it will rank $\bm{\theta_i}$ based on the improvement $\Delta_i = f(\bm{\theta_i})-f(\bm{\theta_e})$. CMA-ME performs a two-stage ranking, where it ranks all solutions that discover empty cells before solutions that improve occupied cells.  

We now show the computation of $\Delta_i$ for CMA-MAE with $\alpha = 1$. If $\bm{\theta_i}$ discovers an empty cell
$\Delta_i = f(\bm{\theta_i}) - t_e$ and by Lemma~\ref{Le_1} $\Delta_i = f(\bm{\theta_i}) - min_{f}$.
If $\bm{\theta_i}$ is mapped to a cell occupied by another solution $\bm{\theta_e}$, $\Delta_i = f(\bm{\theta_i}) - t_e$ and by Lemma~\ref{Le_1} $\Delta_i = f(\bm{\theta_i}) - f(\bm{\theta_e})$.

Comparing the values $\Delta_i$ between the two algorithms we observe the following: (1) If $\bm{\theta_i}$ discovers an empty cell, $\Delta_i = f(\bm{\theta_i}) - min_{f}$ for CMA-MAE. However, $min_{f}$ is a constant and comparison-based sorting is invariant to order preserving transformations~\citep{hansen:cma16}, thus ranking by $\Delta_i = f(\bm{\theta_i}) - min_{f}$ is identical to ranking by $\Delta_i = f(\bm{\theta_i})$ performed by CMA-ME. (2) If $\bm{\theta_i}$ is mapped to a cell occupied by another solution $\bm{\theta_e}$, $\Delta_i = f(\bm{\theta_i}) - f(\bm{\theta_e})$ for both algorithms. (3) Because $min_f$ is an arbitrarily large negative number $f(\bm{\theta}_i)-min_f > f(\bm{\theta_i})-f(\bm{\theta_e})$. Thus, CMA-MAE will always rank solutions that discover empty cells  before solutions that are mapped to occupied cells, identically to CMA-ME. 
\end{proof}


We next provide theoretical insights on how the discount function $f_A$ smoothly increases from a constant function $min_f$ to CMA-ME's discount function as $\alpha$ increases from 0 to 1. We show this for the special case of a fixed sequence of candidate solutions. 


\begin{theorem}
\label{thm_apx:monot}
Let $\alpha_i$ and $\alpha_j$ be two archive learning rates for archives $A_i$ and $A_j$ such that $0 \leq \alpha_i < \alpha_j \leq 1$. For two runs of CMA-MAE that generate the same sequence of $m$ candidate solutions $\{S\} = \bm{\theta_1}, \bm{\theta_2}, ..., \bm{\theta_m}$, it follows that $f_{A_i}(\bm{\theta}) \leq f_{A_j}(\bm{\theta})$ for all $\bm{\theta}\in \mathbb{R}^n$.
\end{theorem}

\begin{proof}
We prove the theorem via induction over the sequence of solution additions. $f_A$ is the histogram formed by the thresholds $t_e$ over all archive cells $e$ in the archive. Thus, we prove $f_{A_i} \leq f_{A_j}$ by showing that $t_e(A_i) \leq t_e(A_j)$ for all archive cells $e$ after $m$ archive additions.


As a basis step, we note that $A_i$ equals $A_j$ as both archives are initialized with $min_f$.

Our inductive hypothesis states that after $k$ archive additions we have $t_e(A_i) \leq t_e(A_j)$, and we need to show that $t_e(A_i) \leq t_e(A_j)$ after solution $\bm{\theta_{k+1}}$ is added to each archive.

Our solution $\bm{\theta_{k+1}}$ has three cases with respect to the acceptance thresholds:

\textbf{Case 1:} $f(\bm{\theta_{k+1}}) \leq t_e(A_i) \leq t_e(A_j)$. The solution is not added to either archive and our property holds from the inductive hypothesis.

\textbf{Case 2:} $t_e(A_i) \leq f(\bm{\theta_{k+1}}) \leq t_e(A_j)$. The solution is added to $A_i$, but not $A_j$, thus $t'_e(A_j) = t_e(A_j)$. We follow the threshold update: $t'_e(A_i) = (1-\alpha_i) t_e(A_i) + \alpha_i f(\bm{\theta_{k+1}})$. Next, we need to show that $t'_e(A_i) \leq t'_e(A_j)$ to complete the inductive step: 

\begin{align*}
(1-\alpha_i) t_e(A_i) + \alpha_i f(\bm{\theta_{k+1}}) &\leq f(\bm{\theta_{k+1}})&\iff&\\
(1-\alpha_i) t_e(A_i) &\leq (1-\alpha_i) f(\bm{\theta_{k+1}})&\iff&\\
t_e(A_i) &\leq f(\bm{\theta_{k+1}})~~~\text{as}~ 1-\alpha_i \geq 0 &   &
\end{align*}

The last inequality holds true per our initial assumption for Case 2. From the inductive hypothesis, we have $f(\bm{\theta_{k+1}}) \leq t_e(A_j) =t'_e(A_j)$.


\textbf{Case 3:} $t_e(A_i) \leq t_e(A_j) \leq f(\bm{\theta_{k+1}})$. The solution $\bm{\theta_{k+1}}$ is added to both archives. We need to show that $t'_e(A_i) \leq t'_e(A_j)$:

\begin{align}
t'_e(A_i) &\leq t'_e(A_j) &\iff \nonumber\\
(1-\alpha_i) t_e(A_i) + \alpha_i f(\bm{\theta_{k+1}}) &\leq (1-\alpha_j) t_e(A_j) + \alpha_j f(\bm{\theta_{k+1}}) 
 \label{eq:case_3a}
\end{align}

We can rewrite Eq.~\ref{eq:case_3a} as:
\begin{align}
(1-\alpha_j) t_e(A_j) - (1-\alpha_i) t_e(A_i) + \alpha_j f(\bm{\theta_{k+1}}) - \alpha_i f(\bm{\theta_{k+1}}) \geq 0
\label{eq:case_3b}
\end{align}

 First, note that:
 \begin{align*}
 (1-\alpha_j) t_e(A_j) - (1-\alpha_i) t_e(A_i) &\geq (1-\alpha_j) t_e(A_i) - (1-\alpha_i) t_e(A_i)\\
 &= (1-\alpha_j - 1 + \alpha_i)t_e(A_i)
 \\  &= (\alpha_i - \alpha_j) t_e(A_i).
 \end{align*}
 Thus:
 \begin{equation}
 (1-\alpha_j) t_e(A_j) - (1-\alpha_i) t_e(A_i) \geq (\alpha_i - \alpha_j) t_e(A_i)
 \label{eq:case_3c}
 \end{equation}
 

From Eq.~\ref{eq:case_3b} and~\ref{eq:case_3c} we have: 

\begin{align*}
  (1-\alpha_j) t_e(A_j)& + \alpha_j f(\bm{\theta_{k+1}}) - (1-\alpha_i) t_e(A_i) - \alpha_i f(\bm{\theta_{k+1}})\\ &\geq (\alpha_i - \alpha_j) t_e(A_i) + (\alpha_j - \alpha_i) f(\bm{\theta_{k+1}})\\
    &=(\alpha_j - \alpha_i)(f(\bm{\theta_{k+1}}) - t_e(A_i))
 \end{align*}
 
  As $\alpha_j > \alpha_i$ and $f(\bm{\theta_{k+1}}) \geq t_e(A_i)$, we have $(\alpha_j - \alpha_i)(f(\bm{\theta_{k+1}}) - t_e(A_i)) \geq 0$. This completes the proof that Eq.~\ref{eq:case_3b} holds. 

 










As all cases in our inductive step hold, our proof by induction is complete.
\end{proof}


Next, we wish to provide insights about the exploration properties of CMA-MAE for an archive learning rate $\alpha$ between 0 and 1, when the objective $f$ is constant. Consider an approximate density descent algorithm that is identical to CMA-ME, but differs by how solutions are ranked. Specifically, the algorithm maintains  an occupancy count $o_e$ for each cell $e$, which is the number of times a solution was generated in that cell. We let the algorithm rank solutions based on the occupancy count of the cell that the solution is mapped to, where solutions that discover less frequently visited cells are ranked higher. An approximate density descent algorithm with this ranking mechanism will search for solutions that cover the least explored regions of the measure space.

\begin{lemma}
The threshold $t_e$ after $k$ additions to cell $e$ forms a strictly increasing  sequence for a constant objective function $f(\bm{\theta})=C$ for all $\bm{\theta} \in \mathbb{R}^n$, when $0<\alpha<1$ and $min_f < C$.
\label{Le_2}
\end{lemma}

\begin{proof}
To show that $t_e$ after $k$ additions to cell $e$ forms a strictly increasing sequence, we write a recurrence relation for $t_e$ after $k$ solutions have been added to cell $e$. Let $t_e(k) = (1-\alpha) t_e(k-1) + \alpha f(\theta_i)$ and $t_e(0) = min_f$ be that recurrence relation. To show the recurrence is an increasing function, we need to show that $t_e(k) > t_e(k-1)$ for all $k \geq 0$.

We prove the inequality via induction over cell additions $k$. As a basis step, we show $t_e(1) > t_e(0)$: $(1-\alpha) min_f + \alpha C > min_f \iff min_f - min_f - \alpha \cdot min_f + \alpha C \iff \alpha C > \alpha \cdot min_f$. As $C > min_f$ and $\alpha > 0$, the basis step holds.

For the inductive step, we assume that $t_e(k) > t_e(k-1)$ and need to show that $t_e(k+1) > t_e(k)$: $t_e(k+1) > t_e(k) \iff (1-\alpha) t_e(k) + \alpha C > (1-\alpha) t_e(k-1) + \alpha C \iff (1-\alpha) t_e(k) > (1-\alpha) t_e(k-1) \iff t_e(k) > t_e(k-1)$.
\end{proof}


\begin{theorem}
\label{thm_apx:density_descent}
The CMA-MAE algorithm optimizing a constant objective function $f(\bm{\theta}) = C$ for all $\bm{\theta} \in \mathbb{R}^n$ is equivalent to the approximate density descent algorithm, when $0 < \alpha < 1$ and $min_f < C$.
\end{theorem}

\begin{proof}
We will prove that for an arbitrary archive $A$ with both the occupancy count for each cell $o_e$ and the threshold value $t_e$ computed with arbitrary learning rate $0 < \alpha < 1$, CMA-MAE results in the same ranking for an arbitrary batch of solutions $\{ \bm{\theta_i} \}$ as the approximate density descent algorithm.


We let $\bm{\theta_i}$ and $\bm{\theta_j}$ be two arbitrary solutions in the batch mapped to cells $e_i$ and $e_j$. Without of loss of generality, we let $o_{e_i} \leq o_{e_j}$. The approximate density descent algorithm will thus rank  $\bm{\theta_i}$ before $\bm{\theta_j}$. We will show that CMA-MAE results in the same ranking. 

If $o_{e_i} \leq o_{e_j}$, and since $t_e$ is a strictly increasing function from Lemma~\ref{Le_2}: $t_{e_i}(o_{e_i})\leq t_{e_j}(o_{e_j})$. We have
$t_{e_i}(o_{e_i})\leq t_{e_j}(o_{e_j}) \iff C-t_{e_i}(o_{e_i})\geq C- t_{e_j}(o_{e_j})$. Thus, the archive improvement by adding $\bm{\theta_i}$ to the archive is larger than the improvement by adding $\bm{\theta_j}$ and CMA-MAE will rank $\bm{\theta_i}$ higher than $\bm{\theta_j}$, identically with density descent. 
\end{proof}


While Theorem~\ref{thm_apx:density_descent} assumes a constant objective $f$, we conjecture that the theorem holds true generally when threshold $t_e$ in each cell $e$ approaches the local optimum within the cell boundaries. 

\begin{conjecture}
\label{conj:density_descent_convex}
The CMA-MAE algorithm becomes equivalent to the density descent algorithm for a subset of archive cells for an arbitrary convex objective $f$, where the cardinality of the subset of cells increases as the number of iterations increases. 
\end{conjecture}

We provide intuition for our conjecture through the lense of the elite hypervolume hypothesis~\citep{vassiliades:gecco18}. The elite hypervolume hypothesis states that optimal solutions for the MAP-Elites archive form a connected region in search space. Later work~\citep{rakicevic2021policy}, connected the elite hypervolume hypothesis to the manifold hypothesis~\citep{fefferman2016testing} in machine learning, stating that the elite hypervolume can be represented by a low dimensional manifold in search space.

For our conjecture, we assume that the elite hypervolume hypothesis holds and there exists a smooth manifold that represents the hypervolume. Next, we assume in the conjecture that $f$ is an arbitrary convex function. As $f$ is convex, early in the CMA-MAE search the discount function $f_A$ will be flat and the search point $\bm{\theta}$ will approach the global optimum following CMA-ES's convergence properties~\citep{hansen1997convergence, hansen2003reducing}, where the precision of convergence is controlled by archive learning rate $\alpha$. By definition, the global optimum $\bm{\theta^*}$ is within the elite hypervolume as no other solution of higher quality exists within its archive cell. Assuming the elite hypervolume hypothesis holds, a subset of adjacent solutions in search space will also be in the hypervolume due to the connectedness of the hypervolume. As $f_A$ increases around the global optimum, we conjecture that the function $f(\bm{\theta^*})-f_A(\bm{\theta^*})$ will form a plateau around the optimum, since it will approach the value $f(\bm{\theta_i})-f_A(\bm{\theta_i})$ of adjacent solutions $\bm{\theta_i}$. By Theorem~\ref{thm_apx:density_descent} we have a density descent algorithm within the plateau, pushing CMA-MAE to discover solutions on the frontier of the known hypervolume.

Finally, we remark that our conjecture implies that $f-f_A$ tends towards a constant function in the limit, resulting in a density descent algorithm across the elite hypervolume manifold as the number of generated solutions approaches infinity. We leave a formal proof of this conjecture for future work.

\section{Theoretical Properties of CMA-MAEGA}
\label{sec:cmamaega-theory}

In this section, we investigate how the theoretical properties of CMA-MAE apply to CMA-MAEGA. While many of the properties are nearly a direct mapping, we note that, while CMA-MAE is equivalent to the single-objective optimization algorithm CMA-ES for $\alpha=0$, there is no single-objective counterpart to CMA-MAEGA. To make the direct mapping easier, we introduce a counterpart: the gradient arborescence ascent algorithm.

The gradient arborescence ascent algorithm is similar to CMA-MEGA, but without an archive. Like CMA-MEGA, the algorithm assumes a differentiable objective $f$ and differentiable measures $\bm{m}$. However, the algorithm leverages the objective and measure function gradients only to improve the optimization of the objective $f$, rather than to find solutions that are diverse with respect to measures $\bm{m}$. As with CMA-MEGA, the gradient arborescence algorithm branches in objective-measure space. However, the algorithm ranks solutions purely by the objective function $f$ and adapts the coefficient distribution $N(\bm{\mu},\Sigma)$ towards the natural gradient of the objective $f$.

Next, we prove properties of CMA-MAEGA that directly follow from the properties of CMA-MAE.

\begin{theorem}
The gradient arborescence ascent algorithm is equivalent to CMA-MAEGA when $\alpha=0$, if gradient arborescence ascent restarts from an archive elite.
\end{theorem}

\begin{proof}
We note that CMA-MAEGA and the gradient arborescence ascent algorithm differ only in how they rank solutions, and we note that the differences between CMA-MAE and CMA-ES mirror the differences between CMA-MAEGA and gradient arborescence ascent algorithm. So by directly adapting the proof of Theorem~\ref{thm_apx:alpha0}, we complete our proof.
\end{proof}

\begin{theorem}
The CMA-MEGA algorithm is equivalent to CMA-MAEGA when $\alpha = 1$ and $min_f$ is an arbitrarily large negative number.
\end{theorem}

\begin{proof}
We note that CMA-MAEGA and the CMA-MEGA algorithm differ only in how they rank solutions and how they update the archive $A$, and we note that the differences between CMA-MAE and CMA-ME mirror the differences between CMA-MAEGA and CMA-MEGA. So by directly adapting the proof of Theorem~\ref{thm_apx:alpha1}, we complete our proof.
\end{proof}

\begin{theorem}
Let $\alpha_i$ and $\alpha_j$ be two archive learning rates for archives $A_i$ and $A_j$ such that $0 \leq \alpha_i < \alpha_j \leq 1$. For two runs of CMA-MAEGA that generate the same sequence of $m$ candidate solutions $\{S\} = \bm{\theta_1}, \bm{\theta_2}, ..., \bm{\theta_m}$, it follows that $f_{A_i}(\bm{\theta}) \leq f_{A_j}(\bm{\theta})$ for all $\bm{\theta}\in \mathbb{R}^n$.
\end{theorem}

\begin{proof}
We note that CMA-MAE and CMA-MAEGA update the archive $A$ in exactly the same way. Therefore, the proof follows directly by adapting the proof of Theorem~\ref{thm_apx:monot} to CMA-MAEGA.
\end{proof}


Next, we wish to show that CMA-MAEGA results in density descent in measure space. However, we need a counterpart to the approximate density descent algorithm we defined in Theorem~\ref{thm_apx:density_descent}. 

Consider an approximate density descending arborescence algorithm that is identical to CMA-MEGA, but differs by how solutions are ranked. Specifically, we assume that this  algorithm maintains an occupancy count $o_e$ for each cell $e$, which is the number of times a solution was generated in that cell. The density descent algorithm ranks solutions based on the occupancy count of the cell that the solution maps to, where solutions that discover less frequently visited cells are ranked higher. The algorithm takes steps in search space $\mathbb{R}^n$ that minimize the approximate density function defined by the archive and adapts the coefficient distribution $N(\bm{\mu},\Sigma)$ towards coefficients that minimize the density function.

\begin{theorem}
The CMA-MAEGA algorithm optimizing a constant objective function $f(\bm{\theta}) = C$ for all $\bm{\theta} \in \mathbb{R}^n$ is equivalent to the approximate density descending arborescence algorithm, when $0 < \alpha < 1$ and $min_f < C$.
\end{theorem}

\begin{proof}
The proof of Theorem~\ref{thm_apx:density_descent} relies only on how CMA-MAE updates the archive $A$ and acceptance threshold $t_e$. The proof of this theorem follows directly by adapting the proof of Theorem~\ref{thm_apx:density_descent} to CMA-MAEGA.
\end{proof}

\section{Derivation of the Conversion Formula for the Archive Learning Rate}
\label{sec:conversion}

In this section, we derive the archive learning rate conversion formula $\alpha_2 = 1 - (1 - \alpha_1)^r$ mentioned in Section~\ref{sec:robustness} of the main paper, where $r$ is the ratio between archive cell counts, and $\alpha_1$ and $\alpha_2$ are archive learning rates for two archives $A_1$ and $A_2$.

Given an archive learning rate $\alpha_1$ for $A_1$, we want to derive an equivalent archive learning rate $\alpha_2$ for $A_2$ that results in robust performance when CMA-MAE is run with either $A_1$ or $A_2$. A principled way to derive a conversion formula for $\alpha_2$ is to look for an invariance property that affects the performance of CMA-MAE and that holds when CMA-MAE generates solutions in archives $A_1$ and $A_2$.



Since CMA-MAE ranks solutions by $f-f_A$, we wish for $f_A$ to increase at the same rate in the two archives. Since $f_A(\bm{\theta})=t_e$, where $t_e$ is the cell that a solution $\bm{\theta}$ maps to, we select the \textit{average value} of the acceptance thresholds $t_e$ over all cells in each archive as our invariant property. 



We assume an arbitrary sequence of $N$ solution additions $\bm{\theta_1}, \bm{\theta_2}, ..., \bm{\theta_N}$, evenly dispersed across the archive cells. We then specify $t_e$ as a function that maps $k$ cell additions to a value $t_e$ in archive cell $e$.\footnote{Here we abuse notation and view $t_e$ as a function instead of threshold for simplicity and to highlight the connection to the threshold value $t_e$.}  Equation~\ref{eq:avg_te} then defines the average value of $t_e$ across the archive after $N$ additions to an archive $A$ with $M$ cells. 


\begin{equation} \label{eq:avg_te}
 \frac{1}{M} \sum_{i=1}^{M}  {t_e \left( \frac{N}{M} \right)}
\end{equation}

Then, equation~\ref{eq:alr_invariance} defines the invariance we want to guarantee between archives $A_1$ and $A_2$.

\begin{equation} \label{eq:alr_invariance}
 \frac{1}{M_1} \sum_{i=1}^{M_1}  {t_e \left( \frac{N}{M_1} \right)} =  \frac{1}{M_2} \sum_{i=1}^{M_2}  {t_e \left( \frac{N}{M_2} \right)}
\end{equation}

In Eq.~\ref{eq:alr_invariance}, we let $M_1$ and $M_2$ the number of cells in archives $A_1$ and $A_2$, and we assume that $M_1$ and $M_2$ divide $N$.
To solve for a closed form of $\alpha_2$ subject to our invariance, we need a formula for the function $t_e$. Similar to Lemma~\ref{Le_2}, we can represent the function $t_e$ as a recurrence relation after adding $k$ solutions to cell $e$ of an archive $A$. 

\begin{align}
t_e(0) &= min_f \nonumber \\
t_e(k) &= (1-\alpha) t_e(k-1) + \alpha f(\bm{\theta_k})
\end{align}

Next, we look to derive a closed form for $t_e(k)$ for an archive $A$ as a way to manipulate Equation~\ref{eq:alr_invariance}. However, solving for $t_e(k)$ when $f$ is an arbitrary function is difficult, because different regions of the archive will change at different rates. Instead, we solve for the special case when $f(\bm{\theta}) = C$ and $min_f < C$, where $C\in \mathbb{R}$ is a constant scalar. To solve for a closed form of the recurrence $t_e(k)$, we leverage the recurrence unrolling method~\citep{graham1989concrete}, allowing us to guess the closed form in Equation~\ref{eq:closed_te}.

\begin{align}
 \label{eq:closed_te}
t_e(1) &= (1-\alpha) t_e(0) + \alpha C = (1-\alpha) min_f + \alpha C \nonumber \\
t_e(2) &= (1-\alpha) t_e(1) + \alpha C = (1-\alpha)[(1-\alpha) min_f + \alpha C] + \alpha C  \nonumber \\
&= (1-\alpha)^2 min_f + (1-\alpha) \alpha C + \alpha C \nonumber \\
t_e(3) &= (1-\alpha) t_e(2) + \alpha C \nonumber \\ 
&= (1-\alpha)[(1-\alpha)^2 min_f + (1-\alpha) \alpha C+ \alpha C] + \alpha C \nonumber \\
&= (1-\alpha)^3 min_f + (1-\alpha)^2 \alpha C + (1-\alpha) \alpha C + \alpha C \nonumber\\
&~~\vdots\nonumber\\
t_e(k) &= (1-\alpha)^k min_f + \sum_{i=0}^{k-1} (1-\alpha)^i \alpha C
\end{align}


We recognize the summation in Equation~\ref{eq:closed_te} as a geometric series. As $0 < \alpha < 1$, we rewrite the summation as follows.

\begin{align}
t_e(k) &= (1-\alpha)^k min_f + \sum_{i=0}^{k-1} (1-\alpha)^i \alpha C \nonumber \\
&= (1-\alpha)^k min_f + \alpha C \left( \frac{1-(1-\alpha)^k}{1-(1-\alpha)} \right) \nonumber \\
&= (1-\alpha)^k min_f + \alpha C \left( \frac{1-(1-\alpha)^k}{\alpha} \right) \nonumber \\
&= (1-\alpha)^k min_f + C-C(1-\alpha)^k \nonumber \\
&= (min_f - C) (1 - \alpha)^k + C \nonumber \\
&= C - (C - min_f)(1 - \alpha)^k
\end{align}

Next, we prove that the closed form we guessed is the closed form of the recurrence relation.

\begin{theorem}
\label{thm:closed_form_te}
The recurrence relation $t_e(0) = min_f$ and $t_e(k) = (1-\alpha) t_e(k-1) + \alpha C$ has the closed form $t_e(k) = C - (C - min_f)(1 - \alpha)^k$, where $0 < \alpha < 1$ and $min_f < C$.
\end{theorem}

\begin{proof}
We show the closed form holds via induction over cell additions $k$.

As a basis step we show that $t_e(0) = C - (C - min_f)(1 - \alpha)^0 = C - (C - min_f) = min_f$. 

For the inductive step, suppose after $j$ insertions into the archive $A$ in cell $e$ our closed form holds. We show that the closed form holds for $j+1$ insertions.

\begin{align}
t_e(j+1) &= (1-\alpha) t_e(j) + \alpha C \nonumber \\
&= (1-\alpha) [C - (C-min_f)(1-\alpha)^j] + \alpha C \nonumber \\
&= C(1-\alpha) - (C-min_f)(1-\alpha)^{j+1} + \alpha C \nonumber \\
&= C - \alpha C + \alpha C - (C-min_f)(1-\alpha)^{j+1} \nonumber \\
&= C - (C-min_f)(1-\alpha)^{j+1} 
\end{align}

As our basis and inductive steps hold, our proof is complete.

\end{proof}

The closed form from Theorem~\ref{thm:closed_form_te} allows us to derive a conversion formula for $\alpha_2$ via our invariance formula in Equation~\ref{eq:alr_invariance}.

\begin{align}
\frac{1}{M_1} \sum_{i=1}^{M_1}  {t_e \left( \frac{N}{M_1} \right)} &=  \frac{1}{M_2} \sum_{i=1}^{M_2}  {t_e \left( \frac{N}{M_2} \right)} \nonumber \\
\frac{M_1}{M_1} \left( C - (C - min_f)(1 - \alpha_1)^\frac{N}{M_1} \right) &= \frac{M_2}{M_2} \left( C - (C - min_f)(1 - \alpha_2)^\frac{N}{M_2} \right) \nonumber \\
(C - min_f)(1 - \alpha_1)^\frac{N}{M_1} &= (C - min_f)(1 - \alpha_2)^\frac{N}{M_2} \nonumber \\
(1 - \alpha_1)^\frac{N}{M_1} &= (1 - \alpha_2)^\frac{N}{M_2} \nonumber \\
(1 - \alpha_1)^\frac{M_2}{M_1} &= (1 - \alpha_2) \nonumber \\
\alpha_2 &= 1 - (1 - \alpha_1)^\frac{M_2}{M_1}
\end{align}

We remark that our conversion formula is not dependent on the number of archive additions $N$. 

Although our conversion formula assumes $f$ to be a constant objective, we conjecture that the formula holds generally for a convex objective $f$. %

\begin{conjecture}
The archive learning rate conversion formula results in invariant behavior of CMA-MAE for two arbitrary archives $A_1$ and $A_2$ with archive resolutions $M_1$ and $M_2$, for a convex objective $f$.
\end{conjecture}

Our intuition is similar to the intuition behind Conjecture~\ref{conj:density_descent_convex}, where we assume the elite hypervolume hypothesis holds~\citep{vassiliades:gecco18}. At the beginning of the CMA-MAE search, $f_A$ is a constant function and \mbox{CMA-MAE} optimizes for the global optimum, following the convergence properties of \mbox{CMA-ES}~\citep{hansen1997convergence, hansen2003reducing}. Eventually, the cells around the global optimum become saturated and the function $f-f_A$ forms a plateau around the global optimum. The invariance described in Eq.~\ref{eq:alr_invariance} implies that the $f_{A_1}$ and $f_{A_2}$ will increase at the same rate within the flat region of the plateau. Let $\bm{\theta_p}$ be an arbitrary solution in the plateau and $\bm{\theta'}$ be a solution on the frontier of the known hypervolume. The plateau of each archive $A_i$ expands when the solutions on the frontier of the elite hypervolume achieve a larger $f(\bm{\theta'})-f_{A_i}(\bm{\theta'})$ than the plateau $f(\bm{\theta_p})-f_{A_i}(\bm{\theta_p})$. We conjecture that the plateau will expand at the same rate in the two archives as $f_{A_1}$ and $f_{A_2}$ increase at the same rate for the plateau region, due to our invariance in Eq.~\ref{eq:alr_invariance}. 

We speculate that our conjecture explains why we observe invariant behavior across archive resolutions in the experiments of Section~\ref{sec:robustness}, even though $f$ is not a constant function in the linear projection and arm repertoire domains.






\section{On Improving the Quality of Latent Space Illumination}
\label{sec:improving-lsi}

In the main latent space illumination (LSI) experiments in section~\ref{sec:experiments}, we showed that CMA-MAE outperformed the other QD algorithms according to standard QD metrics following the exact hyperparameters of prior work~\citep{fontaine2021differentiable}. In this section, we show how to improve the quality of images generated by CMA-MAE in the LSI domain.

In the main experiments, we used latent space illumination as purely an optimization benchmark. However, obtaining high performance on LSI as a benchmark can be a competing objective with producing high quality images.

First, finding solutions that result in a high objective value does not always result in high quality images that match the text prompt. For example, a QD algorithm can find images that result in CLIP reporting a high similarity score by leaving the training distribution of StyleGAN. 

Furthermore, we use the CLIP loss as a measure function, thus a QD algorithm attempts to both decrease and increase the loss function to cover the measure space. Increasing the loss function results in minimizing similarity with the text prompt, which can be attained by unrealistic images.  


In the main LSI experiments all derivative-free QD algorithms would drift out of the latent distribution and produce archives of low image quality. We found that CMA-MAE would stay in the latent distribution longer before drifting out of distribution during exploration, due to the low archive learning rate $\alpha$ prioritizing the objective.

To address drifting out of distribution, we adopt a ``timeout'' restart rule proposed by other evolution strategy-based quality diversity algorithms~\citep{colas2020scaling, paolo2021sparse}. A timeout restart rule runs for a fixed number of iterations before restarting. We add to the basic restart rule of CMA-MAE (Appendix~\ref{sec:hyperparams}) an additional criterion for restarting based on the timeout restart rule. To generate all collages in the paper, we use a timeout of 50 iterations for both CMA-MAE and CMA-MAEGA.

Beyond the specifics of the QD optimization algorithm, many aspects of latent space illumination can be improved. For example, prior work~\citep{frans2021clipdraw} on CLIP guiding single-objective optimization notes that the gradients that CLIP provides can be noisy and recommends data augmentations of the generated images, such as tiling or translating each image, before being passed to CLIP. This change can help smooth the gradients for gradient descent optimizers like Adam and make the generated images retained by the archive match their text prompts more accurately.

Prior work on optimizing the latent space of VQ-GAN~\citep{crowson2022vqgan} also notes that CLIP will not always provide smooth optimization gradients, nor accurate objective values. The authors recommend creating a batch of random cutouts of the generated image and passing those images to CLIP, which produces smoother gradients and objective values. The paper also recommends regularizing the latent codes so that they become attracted to the Gaussian ball that captures the training disribution of the GAN. Both these techniques could improve the qualitative performance of latent space illumination.

Finally, we used the first version of StyleGAN~\citep{karras2019style} that was used in previous work~\citep{fontaine2021differentiable}. We anticipate recent versions of StyleGAN~\citep{karras:stylegan2,karras:stylegan2ada,karras:stylegan3} to further improve the quality of the generated images.


\section{Additional Results}
\label{sec:additional-results} 

\begin{table*}[t]
\centering
\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{Linear Projection (sphere)}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 41.64 $\pm$ 0.06& 50.80  $\pm$ 0.09\% &  98.63 $\pm$ 0.01 \\
MAP-Elites (line)  & 49.07 $\pm$ 0.03& 60.42 $\pm$ 0.05\% &  99.43 $\pm$ 0.01\\
CMA-ME    &  36.50 $\pm$ 0.31  & 42.82 $\pm$ 0.40\%  & 100.00 $\pm$ 0.00\\
CMA-ME$^*$  & 45.07 $\pm$ 0.07 & 55.11 $\pm$ 0.10\%  & 99.23  $\pm$ 0.02 \\
CMA-MAE    &  64.86 $\pm$ 0.04  & 83.31 $\pm$ 0.07\% & 99.59 $\pm$ 0.01\\
CMA-MEGA &  75.32 $\pm$ 0.00& 100.00 $\pm$ 0.00\% & 100.00 $\pm$ 0.00 \\
CMA-MAEGA & 75.39 $\pm$ 0.00 & 100.00 $\pm$ 0.00\% & 100.00 $\pm$ 0.00\\
 \end{tabular}
}
\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
\hline
             & \multicolumn{3}{c}{Linear Projection (Rastrigin)}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 31.43 $\pm$ 0.07 & 47.88 $\pm$ 0.12\% &  82.16 $\pm$ 9,11 \\
MAP-Elites (line)  & 38.29 $\pm$ 0.05&  56.51 $\pm$ 0.09\% & 85.71 $\pm$ 0.07 \\
CMA-ME    & 38.02 $\pm$ 0.11 & 53.09 $\pm$ 0.16\%  & 97.59 $\pm$ 0.06\\
CMA-ME$^*$  & 35.06 $\pm$ 0.05 & 53.01 $\pm$ 0.12\%  & 83.47 $\pm$ 0.08 \\
CMA-MAE    & 52.65 $\pm$  0.06 &  80.46 $\pm$ 0.11\% & 95.90 $\pm$ 0.25\\
CMA-MEGA & 63.07 $\pm$ 0.00 & 100.00 $\pm$ 0.00\% & 100.00 $\pm$ 0.00\\
CMA-MAEGA & 63.06 $\pm$ 0.00 & 100.00 $\pm$ 0.00\% & 100.00 $\pm$ 0.00\\
  \bottomrule
\end{tabular}
}

\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{Linear Projection (plateau)}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 47.07 $\pm$ 0.17 & 47.07 $\pm$ 0.17\% &   100.00 $\pm$ 0.00\\
MAP-Elites (line)  & 52.20 $\pm$ 0.19 & 52.20 $\pm$ 0.19\% & 100.00  $\pm$ 0.00 \\
CMA-ME    & 34.54 $\pm$  0.35 & 34.54 $\pm$ 0.35\% & 100.00 $\pm$ 0.00 \\
CMA-ME$^*$ & 51.11 $\pm$ 0.25 & 51.11 $\pm$ 0.25\%  &  100.00 $\pm$ 0.00 \\
CMA-MAE    & 79.27 $\pm$  0.21 & 79.29 $\pm$ 0.21\% & 100.00 $\pm$ 0.00 \\
CMA-MEGA & 100.00 $\pm$ 0.00 & 100.00 $\pm$ 0.00\% & 100.00 $\pm$ 0.00 \\
CMA-MAEGA & 100.00 $\pm$ 0.00  & 100.00  $\pm$ 0.00\% & 100.00 $\pm$ 0.00 \\
  \bottomrule
\end{tabular}
}

\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{Arm Repertoire}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 71.40 $\pm$ 0.03 & 74.09 $\pm$ 0.04\% &  97.38 $\pm$ 0.03 \\
MAP-Elites (line)  & 74.55 $\pm$ 0.02 & 75.61 $\pm$ 0.02\%  & 99.16 $\pm$ 0.01 \\
CMA-ME    &  75.82 $\pm$ 0.11 &  75.89 $\pm$ 0.11\% & 100.00 $\pm$ 0.00\\
CMA-ME$^*$    & 75.68 $\pm$ 0.04 & 76.13 $\pm$ 0.03\%  & 99.78 $\pm$ 0.01 \\
CMA-MAE    & 79.03 $\pm$ 0.02 & 79.24 $\pm$ 0.02\%  & 99.93 $\pm$ 0.00\\
CMA-MEGA &  75.21 $\pm$ 0.13 & 75.25 $\pm$ 0.13\% & 100.00 $\pm$ 0.00\\
CMA-MAEGA & 79.27 $\pm$ 0.02 & 79.35 $\pm$ 0.02\% & 100.00 $\pm$ 0.00\\
  \bottomrule
\end{tabular}
}
\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{Latent Space Illumination}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 12.85 $\pm$ 0.10 &  19.42 $\pm$ 0.16\% &  71.42  $\pm$ 0.14 \\
MAP-Elites (line)  & 14.40 $\pm$ 0.09 & 21.11 $\pm$ 0.11\% & 73.04 $\pm$ 0.05\\
CMA-ME    & 14.00  $\pm$  0.62 & 19.57 $\pm$ 0.90\% & 74.11 $\pm$ 0.08\\
CMA-MAE    & 17.67 $\pm$ 0.27  & 25.08  $\pm$ 0.40\% & 73.48 $\pm$ 0.18\\
CMA-MEGA &  16.08 $\pm$ 0.37 & 22.58  $\pm$ 0.57\% & 74.95  $\pm$ 0.27 \\
CMA-MAEGA & 16.20 $\pm$ 0.41 &  23.83 $\pm$ 0.46\% & 75.52 $\pm$ 0.22\\
  \bottomrule
\end{tabular}
}

\caption{Results: The QD-score, coverage, and best solution after 10,000 iterations for each algorithm and domain with standard errors.
 }
\label{tab:results_detailed}
\end{table*}

\subsection{Generated Archives and Additional Metrics}

Table~\ref{tab:results_detailed} presents the values of the QD-score, coverage, and best solution for each algorithm and domain. We used $\alpha=0.01$ for CMA-MAE, identically to the main experiments. 
Similarly to previous work~\citep{fontaine2021differentiable}, we disambiguate the quality of solutions found and coverage by showing for MAP-Elites, MAP-Elites (line), CMA-ME and CMA-MAE the percentage of cells (y-axis) that have objective value greater than the threshold specified in the x-axis (Fig.~\ref{fig:cumulative_plots}).


\subsection{Example Archives}

Fig.~\ref{fig:results_all_sphere}  -\ref{fig:results_all_lsi} show example archives for each algorithm and domain. 

 \begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/cumulative_plots.pdf}
\caption{The percentage of cells in the y-axis with objective values larger than or equal to a threshold specified in the x-axis, with 95$\%$ confidence intervals. The percentage is the number of filled cells (filtered by the threshold) over the archive size. A larger area under each curve indicates better performance.}
\label{fig:cumulative_plots}
\end{figure}


\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/sphere_heatmaps.pdf}
\caption{Example archives for each algorithm for the linear projection (sphere) domain.}
\label{fig:results_all_sphere}
\end{figure*}

\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/rast_heatmaps.pdf}
\caption{Example archives for each algorithm for the linear projection (Rastrigin) domain.}
\label{fig:results_all_rast}
\end{figure*}

\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/plateau_heatmaps.pdf}
\caption{Example archives for each algorithm for the linear projection (plateau) domain.}
\label{fig:results_all_plateau}
\end{figure*}

\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/arm_heatmaps.pdf}
\caption{Example archives for each algorithm for the arm repertoire domain.}
\label{fig:results_all_arm}
\end{figure*}


\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/lsi_heatmaps.pdf}
\caption{Example archives for each algorithm for the latent space illumination domain.}
\label{fig:results_all_lsi}
\end{figure*}


\subsection{Additional Experiments in the LSI Domain}


We include the same additional experiments in the LSI domain as with previous work~\citep{fontaine2021differentiable}. The first additional experiment has  objective prompt ``A photo of Jennifer Lopez'' and measure prompts ``A small child.'' and ``A woman with long blonde hair.'' The second has objective prompt ``A photo of Elon Musk'' and measure prompts ``A person with red hair.'' and ``A man with blue eyes.'' Table~\ref{tab:lsi_additional_runs} shows the results of the additional runs, as well as the Beyonc\'e run from the main paper.

\subsection{Additional Results for Varying Resolutions}
Fig.~\ref{fig:alpha_results_all} shows the QD-score and coverage of CMA-MAE with resolution-dependent archive learning rate and the baselines, for each benchmark domain. For CMA-MAE, we set the resolution dependent archive learning rate $\alpha$ using the conversion formula from Appendix~\ref{sec:conversion}, with $\alpha_1 = 0.01$ for resolution $100 \times 100$.



\begin{figure*}[t!]
\includegraphics[width=1.0\columnwidth]{figures/alpha_results_all.pdf}
\caption{Final QD-score and coverage of each algorithm for 25 different archive resolutions.}
\label{fig:alpha_results_all}
\end{figure*}



\begin{table}
\centering


\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{LSI: Beyonc\'e}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  & 12.85 $\pm$ 0.10 &  19.42 $\pm$ 0.16\% &  71.42  $\pm$ 0.14 \\
MAP-Elites (line)  & 14.40 $\pm$ 0.09 & 21.11 $\pm$ 0.11\% & 73.04 $\pm$ 0.05\\
CMA-ME    & 14.00  $\pm$  0.62 & 19.57 $\pm$ 0.90\% & 74.11 $\pm$ 0.08\\
CMA-MAE    & 17.67 $\pm$ 0.27  & 25.08  $\pm$ 0.40\% & 73.48 $\pm$ 0.18\\
CMA-MEGA &  16.08 $\pm$ 0.37 & 22.58  $\pm$ 0.57\% & 74.95  $\pm$ 0.27 \\
CMA-MAEGA & 16.20 $\pm$ 0.41 &  23.83 $\pm$ 0.46\% & 75.52 $\pm$ 0.22\\
  \bottomrule
\end{tabular}
r}


\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{LSI: Jennifer Lopez}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  &  12.51 $\pm$ 0.28 & 19.18 $\pm$ 0.48\%    &  70.87 $\pm$ 0.27 \\
MAP-Elites (line)  &  14.73 $\pm$ 0.06 & 21.60 $\pm$ 0.08\%  &  73.50 $\pm$ 0.13\\
CMA-ME    &  15.24 $\pm$ 0.37 & 20.86 $\pm$ 0.50\% & 75.39 $\pm$ 0.09\\
CMA-MAE   & 18.33 $\pm$ 0.16 & 25.42 $\pm$ 0.24\% & 75.10 $\pm$ 0.17\\
CMA-MEGA  &  17.06 $\pm$  0.10 & 23.40 $\pm$ 0.14\% & 76.02 $\pm$ 0.08\\
CMA-MAEGA &  16.45 $\pm$ 0.27 & 23.60 $\pm$ 0.49\% & 76.42 $\pm$ 0.13\\
  \bottomrule
\end{tabular}
}~\\

\resizebox{0.7\linewidth}{!}{
\begin{tabular}{l|rrr}
             & \multicolumn{3}{c}{LSI: Elon Musk}   \ \\ 
    \toprule
Algorithm     & QD-score & Coverage  & Best \\
    \midrule
MAP-Elites  &  13.88 $\pm$ 0.11 & 23.15 $\pm$ 0.14\%    &  69.76 $\pm$ 0.07\\
MAP-Elites (line)  & 16.54 $\pm$ 0.28& 25.73 $\pm$ 0.31\%  &  72.63 $\pm$ 0.28\\
CMA-ME    & 18.96 $\pm$ 0.17 & 26.18 $\pm$ 0.24\% & 75.84 $\pm$ 0.10\\
CMA-MAE & 22.10 $\pm$ 0.31 & 30.89 $\pm$ 0.44\% & 75.25 $\pm$ 0.20\\
CMA-MEGA   & 21.82 $\pm$ 0.18 & 30.73 $\pm$ 0.15\% & 76.89 $\pm$ 0.15\\
CMA-MAEGA & 19.99 $\pm$ 0.21 & 30.12 $\pm$ 0.42\% & 77.25 $\pm$ 0.18\\
  \bottomrule
\end{tabular}
}

\caption{Results from additional runs for Beyonc\'e, Jennifer Lopez,  and Elon Musk} 
\label{tab:lsi_additional_runs}
\end{table}



\subsection{Qualitative Results in the LSI Domain} 
We present two additional example collages for CMA-MAE (Fig.~\ref{fig:mae_lopez},~\ref{fig:mae_musk}) and for CMA-MAEGA (Fig.~\ref{fig:maega_lopez},~\ref{fig:maega_musk}). We also include collages of each run of CMA-MAE and CMA-MEGA for all three LSI domains in the link: \url{https://www.dropbox.com/sh/7e22190k3p4zh69/AACcAKV7_Xgi4IMrhzxkCz5ca?dl=0}. 

\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/collage_lopez_mae.pdf}
\caption{A latent space illumination collage generated by CMA-MAE for the objective ``A photo of Jennifer Lopez.'' and for measures ``A small child.'' and ``A woman with long blonde hair.'' The axes values indicate the score returned by the CLIP model, where lower score indicates a better match.}
\label{fig:mae_lopez}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/collage_musk_mae.pdf}
\caption{A latent space illumination collage generated by CMA-MAE for the objective ``Elon Musk with short hair.'' and for measures ``A man with blue eyes.'' and ``A person with red hair.'' The axes values indicate the score returned by the CLIP model, where lower score indicates a better match.}
\label{fig:mae_musk}
\end{figure}


\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/collage_lopez_maega.pdf}
\caption{A latent space illumination collage generated by CMA-MAEGA for the objective ``A photo of Jennifer Lopez.'' and for measures ``A small child.'' and ``A woman with long blonde hair.'' The axes values indicate the score returned by the CLIP model, where lower score indicates a better match.}
\label{fig:maega_lopez}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=1.0\linewidth]{figures/collage_musk_maega.pdf}
\caption{A latent space illumination collage generated by CMA-MAEGA for the objective ``Elon Musk with short hair.'' and for measures ``A man with blue eyes.'' and ``A person with red hair.'' The axes values indicate the score returned by the CLIP model, where lower score indicates a better match.}
\label{fig:maega_musk}
\end{figure}
