\section[Using HyperNOMAD]{Using {\sf HyperNOMAD}\xspace}
\label{sec-appendix}

{\sf HyperNOMAD}\xspace is a {\sf C++} and {\sf Python} package dedicated to the hyperparameter optimization of deep neural networks. The package contains a blackbox specifically designed for this problematic and provides a link with the {\sf NOMAD}\xspace software~\cite{Le09b} used for the optimization. The blackbox takes as inputs the hyperparameters discussed in Section~\ref{sec:hpos}, builds a corresponding deep neural network in order to train, validate and test it on a specific data set before returning the test accuracy as a mesure of performance. {\sf NOMAD}\xspace is then used to minimize this error. The following appendix provides an overview of how to use the {\sf HyperNOMAD}\xspace package.

\subsection*{Prerequisites}

{\sf HyperNOMAD}\xspace relies on:
\begin{itemize}
	\item A compiled version of the {\sf NOMAD}\xspace software available at \url{https://www.gerad.ca/nomad/} for the optimization;
	\item The {\sf PyTorch}\xspace library available at \url{https://pytorch.org/} for modeling the neural network within the blackbox;
	\item A version of {\sf Python} superior to~3.6;
	\item A version of {\sf gcc} superior to~3.8.
\end{itemize}

\subsection*{Installation of {\sf HyperNOMAD}\xspace}

{\sf HyperNOMAD}\xspace is available at
\href{https://github.com/DouniaLakhmiri/HyperNOMAD}{
{\tt https://github.com/DouniaLakhmiri/Hyper} {\tt NOMAD}
}. The user must produce the executable \texttt{hypernomad.exe} using the provided makefile as follows:

\begin{lstlisting}
> make
	building HYPERNOMAD ...

    	To be able to run the example
        the HYPERNOMAD_HOME environment variable
        must be set to the HyperNOMAD home directory
\end{lstlisting}

\vspace*{5mm}
When the compilation is successful, a message appears asking to set an environment variable {\tt HYPERNOMAD\_HOME} which can be done by adding a line in the
{\tt .profile} or {\tt .bashrc} files:

\begin{lstlisting}
export HYPERNOMAD_HOME=hypernomad_directory
\end{lstlisting}

\vspace*{5mm}
The user can check that the installation is successful by trying to run the command:\\

\begin{lstlisting}
> $HYPERNOMAD_HOME/bin/hypernomad.exe -i

--------------------------------------------------
  HYPERNOMAD - version 1.0
--------------------------------------------------
  Using Nomad version 3.9.0 - www.gerad.ca/nomad
--------------------------------------------------

Run           : hypernomad.exe parameters_file
Info          : hypernomad.exe -i
Help          : hypernomad.exe -h
Version       : hypernomad.exe -v
Usage         : hypernomad.exe -u
Neighboors    : hypernomad.exe -n parameters_file

\end{lstlisting}

\subsection*{Using {\sf HyperNOMAD}\xspace}

The next phase is to create a parameter file that contains the necessary information to specify the classification problem, the search space and the initial starting point. {\sf HyperNOMAD}\xspace allows for a good flexibility of tuning a convolutional network by considering multiple aspects of a network at once such as the architecture, the dropout rate, the choice of the optimizer and the hyperparameters related to the optimization aspect (learning rate, weight decay, momentum, etc.), the batch size, etc. The user can choose to optimize all these aspects or select a few and fix the others to certain values. The user can also change the default range of each hyperparameter. This information is passed through the parameter file by using a specific syntax where ``LB" represents the lower bound and ``UB" the upper bound.

\begin{lstlisting}
KEYWORD  INITIAL_VALUE  LB  UB  FIXED/VAR
\end{lstlisting}

\vspace*{5mm}While the hyperparameters have default values in {\sf HyperNOMAD}\xspace, the data set must be explicitly provided by the user in a separate file in order to specify the considered optimization problem. The following section explains how to specify the necessary parameter file before running an optimization.

\subsubsection*{Choosing a data set}

The library can be used on different data sets whether they are already incorporated in {\sf HyperNOMAD}\xspace, such as the ones listed in Table~\ref{tab:datasets}, or are provided by the user. In the latter case, please refer to the user guide in \url{https://hypernomad.readthedocs.io/en/latest/} for details on how to link a personal data set to the library. The rest of the section describes how to run an optimization on a data set provided with {\sf HyperNOMAD}\xspace.

Because of the nature of the applications considered by {\sf HyperNOMAD}\xspace, the computing time can become constraining, especially during the training phase of each configuration, which is why ``{\sf TOYMNIST}" is created as a subset of MNIST containing 300 training images, 100 for the validation and another 100 for testing. It is added to the package for experimenting with {\sf HyperNOMAD}\xspace without having to wait several hours for each blackbox evaluation.

\subsubsection*{Specifying the search space}

In order to specify the problem to optimize and its parameters, the user must provide a parameter file that contains all the necessary informations to run an optimization.  As shown below, the parameter file consists of a list of keywords, each corresponding to a hyperparameter, and the values that the user wishes to attribute them. Some of these key words are mandatory such as the data set, in order to specify the problem, and the number of blackbox evaluations. Other keywords are optional and have default values if they do not appear on the parameter file. Table~\ref{tab:keywords} summarizes all the possible keywords with their default values and ranges. The user can change the lower and upper bounds of a hyperparameter and decide to maintain a hyperparameter to a fixed value during the entire optimization.

Below is a first example of a parameter file that corresponds to the one provided
in {\tt \$HYPERNOMAD\_HOME/examples/mnist\_first\_example.txt}.
First, the MNIST data set is chosen and {\sf HyperNOMAD}\xspace is allowed to try a maximum of 100 configurations. Then, the number of convolutional layers is fixed throughout the optimization to five. The two ``{\tt -}" appearing after the ``{\tt 5}" mean that the default lower and upper bounds are maintained. The kernels, number of fully connected layers, and activation function, are respectively initialized to 3, 6, and 2 (Sigmoid) and the dropout rate is initialized to 0.6 with a new lower bound of 0.3 and upper bound of 0.8  instead of the default range of [0;1]. Finally, all the remaining hyperparameters from Table~\ref{tab:keywords} that are not explicitly mentioned in this file are fixed to their default values.

\begin{lstlisting}
# Mandatory information
DATASET                 MNIST
MAX_BB_EVAL             100

# Optional information
NUM_CON_LAYERS          5  -  -  FIXED # The initial value is fixed
                                       # lower and upper bounds have
                                       # no influence when parameter
                                       # is fixed.

KERNELS                 3 # Only the initial value is set (not fixed)
                          # the lower bound and upper bound
                          # have default values.

NUM_FC_LAYERS           6
ACTIVATION_FUNCTION     2
DROPOUT_RATE            0.6  0.3 0.8  # The lower and upper bounds
                                      # are set to values that are not
                                      # the default ones
REMAINING_HPS           FIXED
\end{lstlisting}

\vspace*{5mm}
Below is a second example of a parameter file where the user is only interested in optimizing the fully connected block of a CNN on a the MNIST data set. All the remaining aspects of the network are fixed to their default values throughout the execution of {\sf HyperNOMAD}\xspace. The optimization starts from a point with 10 fully connected layers of the same size of 500 neurones. This parameter file is provided with the package in \texttt{\$HYPERNOMAD\_HOME/examples/mnist\_fc\_optim.txt}.

\begin{lstlisting}
# Mandatory information
DATASET                 MNIST
MAX_BB_EVAL             150

# Optional information
NUM_FC_LAYERS           10  # Initial value is set to 10
                            # the lower and upper bounds are
                            # the default ones

SIZE_FC_LAYER           500  -  2000 # Initial value is set to 500
                                     # the lower bound is the default one
                                     # the upper bound in now 2000

REMAINING_HPS           FIXED
\end{lstlisting}

\vspace*{5mm}Finally, below is a minimal parameter file where only the mandatory information is specified. The execution of {\sf HyperNOMAD}\xspace starts from the default starting point and all the hyperparameters of Table~\ref{tab:keywords} can be changed. The last line of this file can actually be removed without changing the behavior of {\sf HyperNOMAD}\xspace since the default value for \texttt{REMAINING\_HPS} is set to \texttt{VAR}. Executing {\sf HyperNOMAD}\xspace with this file should return the same values obtained in Figure~\ref{fig:comp}. This file is provided in \texttt{\$HYPERNOMAD\_HOME/examples/cifar10\_default.txt}.

\begin{lstlisting}
# Mandatory information
DATASET                 CIFAR10
MAX_BB_EVAL             100

REMAINING_HPS           VAR
\end{lstlisting}

\subsubsection*{Running an execution}

The user can run the previous example  by executing the following command from the {\tt examples} directory:
\begin{lstlisting}[breaklines=true]
> $HYPERNOMAD_HOME/bin/hypernomad.exe ./mnist_fc_optim.txt
\end{lstlisting}

During the optimization, a window appears to plot the training and validation accuracies of the network corresponding to the current point at each epoch as shown in Figure~\ref{fig:training}. When the optimization is done, {\sf HyperNOMAD}\xspace produces the two files \texttt{history.txt} and \texttt{stats.txt}. The first contains each evaluated point and the corresponding testing accuracy, and the second contains the list of successful points.

\begin{figure}[ht]
   \centering
       \includegraphics[page=1,width=.5\textwidth]{training.png} \\
 \caption{Example of a window that appears during one evaluation of the blackbox in {\sf HyperNOMAD}\xspace. This figure shows in real time the training and validation accuracies of the current evaluated set of hyperparameters, at each epoch.}
 \label{fig:training}
\end{figure}

\begin{table}[ht]
\caption{Keywords for the {\sf HyperNOMAD}\xspace parameters file.}
  \centering
  \resizebox{0.89\textwidth}{!}{
   \begin{tabular}{lccr} \toprule
   {Name } & {Description} & \makecell{Default \\ value} & \makecell{Scope} \\ \midrule
   {\tt DATASET} 	   	&   \makecell{Name of the data set \\used for the optimization}	     & No default	 & \makecell{ A data set from Table~\ref{tab:datasets} or\\ {\tt CUSTOM} for a custom data set}\\ \midrule
   {\tt NUMBER\_OF\_CLASSES}&   \makecell{Number of classes of \\ the classification problem}	     & \makecell{No default. Only use if \\ {\tt DATASET} $=$ {\tt CUSTOM}}	          & $\mathbb{N}$ \\ \midrule
   {\tt MAX\_BB\_EVAL}  	   	&  \makecell{Maximum number of calls\\ to the blackbox}	     &  No default  & [1; $\infty$]\\ \midrule
   {\tt NUM\_CON\_LAYERS} 	   		&  Number of convolutional layers 	     & {\tt 2}	          & [0;100]\\ \midrule
   {\tt OUTPUT\_CHANNELS}	& \makecell{Number of output channels\\ for each convolutional layer}  & {\tt 6} & [1;100] \\ \midrule
   {\tt KERNELS}	   		&  \makecell{Size of the kernel applied \\to each convolutional layer}	     & {\tt 5}	          & [1;20]\\ \midrule
   {\tt STRIDES} 	   		&  \makecell{Step of the kernel for\\ each convolutional layer}     & {\tt 1}	          & [1;3]\\ \midrule
   {\tt PADDINGS} 	   		&  \makecell{Size of the padding for \\  each convolutional layer}     & {\tt 0}	          & [0;2]\\ \midrule
   {\tt DO\_POOLS}			&  \makecell{Apply a pooling after \\  each convolutional layer}	     & {\tt 0}	          & $\{$0,1$\}$\\  \midrule
   {\tt NUM\_FC\_LAYERS} 	   		&  Number of fully connected layers 	     & {\tt 2}	          & [0;500]\\ \midrule
   {\tt SIZE\_FC\_LAYER} 	& Size of each fully connected layer & {\tt 128} & [1;1,000]\\ \midrule
   {\tt BATCH\_SIZE} 	& \makecell{Size of batch for the \\mini-batch gradient}	 & {\tt 128} & [1;400] \\ \midrule
   {\tt OPTIMIZER\_CHOICE} & \makecell{Optimizer to use \\ from Table~\ref{tab:optimizer}} & {\tt 3} & $\{$1,2,3,4$\}$ \\ \midrule
   {\tt OPT\_PARAM\_1} & \makecell{Learning rate} & {\tt 0.1}  & [0;1]\\ \midrule
   {\tt OPT\_PARAM\_2} & \makecell{Second hyperparameter \\ related to the optimizer.} &  {\tt 0.9}  & [0;1]\\ \midrule
   {\tt OPT\_PARAM\_3} & \makecell{Third hyperparameter \\related to the optimizer.} &  {\tt 0.005} &  [0;1]\\ \midrule
   {\tt OPT\_PARAM\_4} & \makecell{Fourth hyperparameter \\related to the optimizer.} &  {\tt 0}  &  [0;1]\\ \midrule
   {\tt DROPOUT\_RATE} & \makecell{Probability that a node \\ will be dropped out} & {\tt 0.5} & [0;0.95] \\ \midrule
   {\tt ACTIVATION\_FUNCTION} &  \makecell{Choice of the activation function\\ from ReLU, Sigmoid and Tanh} & {\tt 1} & $\{$1,2,3$\}$ \\ \midrule
   {\tt REMAINING\_HPS} & \makecell{ Allows to fix or to vary all the \\  hyperparameters not explicitly \\ mentioned in the parameter file } & {\tt VAR} &
   $\{$ {\tt FIXED} , {\tt VAR}$\}$ \\
   \bottomrule
   \end{tabular}}
   \label{tab:keywords}
\end{table}

\end{document}

\section{Introduction}

Neural networks are mathematical structures used to solve supervised classification problems such as images, sounds and speech, to name a few. In the recent years, neural networks gained in popularity and were declined in different versions: deep, convolutional, recurrent, etc. in order to adapt to specific problematics. This popularity is due to the emergence of large size databases and the development of computational power of contemporary machines, through the use of GPUs in particular. These favorable conditions have allowed neural networks to learn complex structures and achieve a level of precision that can surpass human performance across multiple instances such as robotics~\cite{levine2018learning}, medical diagnostics~\cite{litjens2017survey}, and more.

However, the performance of a neural network is strongly linked to its structure and to the values of the parameters of
the optimization algorithm used to minimize the error between the predictions of the network and the data during its training.
The choices of the neural network hyperparameters can greatly affect its ability to learn from the training data and to generalize with new data. The algorithmic hyperparameters of the optimizer must be chosen a priori and cannot be modified during optimization. Hence, to obtain a neural network, it is necessary to fix several hyperparameters of various types: real, integer and categorical. A variable is categorical when it describes a class, or category, without a relation of order between these categories. The search for an optimal configuration is a very slow process that, along with the training, takes up the majority of the time when developing a network for a new application. It is a relatively new problem that is often solved randomly or empirically.

Derivative free optimization (DFO)~\cite{AuHa2017, CoScVibook} is the field that aims to solve optimization problems where the derivatives are unavailable, although they might exist. This is the case for example when the objective and/or constraints functions are non differentiable, noisy or expensive to evaluate. In addition, the evaluation in some points may fail especially if the values of the objective and/or contraints are the outputs of a simulation or an experience. Blackbox optimization (BBO) is a subfield of DFO where the derivatives do no exist and the problem is modeled as a blackbox. This term refers to the fact that the computing process behind the output values is unknown. The general DFO problem is described as follows:
\begin{align*}
	\min_{x  \in \Omega} f(x)
\end{align*}
where $f$ is the objective function to minimize over the domain $\Omega$.

There are two main classes of DFO methods: model-based and direct search methods.
The first uses the value of the objective and/or the constraints at some already evaluated points to build a model able to guide the optimization by relying on the predictions of the model. For example, this class includes methods based on trust regions~\cite[Chapter~10]{CoScVibook} or interpolations models~\cite{bobyqa}. This differentiates them from direct search methods~\cite{HoJe61a} that adopt a more straightforward strategy to optimize the blackbox. At each iteration, direct search methods generate a set of trial points that are compared to the ``best solution" available. For example, the GPS algorithm~\cite{Torc97a} defines a mesh on the search space and determines the next point to evaluate by choosing a search direction. DFO algorithms usually include a proof of convergence that ensures a good quality solution under certain hypotheses on the objective function. BBO algorithms extend beyond this scope by including heuristics such as evolutionary algorithms, sampling methods and so on.

In~\cite{AuDaOr2014, AuOr06a}, the authors explain how a hyperparameter optimization (HPO) problem can be seen as a blackbox one. Indeed, the HPO problem is equivalent to a blackbox that takes the hyperparameters of a given algorithm and returns some measure of performance defined in advance such as the time of resolution, the value of the best point found or the number of solved problems. In the case of neural networks, the blackbox can return the accuracy on the test data set as a mesure of performance. With this formulation, DFO techniques can be applied to solve the original HPO problem.

This work presents {\sf HyperNOMAD}\xspace, a package that applies {MADS}\xspace, a direct search method behind the {\sf NOMAD}\xspace software, to tune the hyperparameters that affect the architecture and the learning process of a deep neural network.
Figure~\ref{fig-chart} illustrates the workflow when solving HPO problems with {\sf HyperNOMAD}\xspace.
For a given set of hyperparameters, the construction of the network, the network training, validation and testing, are all wrapped as a single blackbox evaluation. One specificity of {\sf HyperNOMAD}\xspace is its ability to explore a large search space by exploiting categorical variables.

\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{hypernomad_acc.pdf}
\caption{The {\sf HyperNOMAD}\xspace workflow.}
\label{fig-chart}
\end{figure}

The manuscript is structured as follows. Section~\ref{sec-lit} presents and discusses some of the main approaches used to solve the HPO problem of neural networks.
In Section~\ref{sec:exp}, the experimental setup is explicitly defined, and the instances used to test the proposed approach are presented.
Section~\ref{sec:hypernomad} introduces the {\sf HyperNOMAD}\xspace package and gives an overview of  {MADS}\xspace, the algorithm selected to carry out the optimization task including the handling of categorical variables. Computational results are provided and discussed in Section~\ref{sec-results}.
Finally, Appendix~\ref{sec-appendix} describes the basic usage of {\sf HyperNOMAD}\xspace.

\section{Literature review}
\label{sec-lit}

Tuning the hyperparameters of a deep neural network is a critical and time consuming process that was mainly done manually relying on the knowledge of the experts. However, the rising popularity of deep neural networks and their usage for diverse applications called for the automatization of this process in order to adapt to each problematic.

The hyperparameters that define a deep neural network can be separated into two categories: The ones that define the architecture of the network and the ones that affect the optimization process of the training phase. Tuning the hyperparameters of the first category alone has led to a separate field of research called Neural Architecture Search (NAS)~\cite{elsken2018neural} that allowed to achieve state of the art performances~\cite{real2018regularized, zoph2016neural} on some benchmark problems, although at a massive computational cost of 800 GPUs for a few weeks. Typically, one would perform a NAS first and then start tuning the other hyperparameters with the optimized architecture.
 However, Zela et al.~\cite{zela2018towards} argue that this separation is not optimal since the two aspects are not entirely independent from one another. Therefore, the proposed research considers both aspects at once.

One of the first scientific approach used to tackle the HPO problem of neural networks is the grid search. This method consists of discretizing the hypercube defined by the range of each hyperparameter and then evaluating each points on the grid. This technique is still used today and is implemented in several HPO libraries such as {\sf scikit-learn}  and {\sf Spearmint}~\cite {scikit-learn, Bayesian_HPO}. It has the advantage of being easy to understand, implement and parallelize. However, it becomes very expensive when training large networks, which is the case of deep neural networks, or when one seeks to optimize several hyperparameters at once. In addition, the grid search ignores the impact of each hyperparameter on the overall performance of the network.

To avoid the drawbacks of the grid search, an alternative is to use random search~\cite{bergstra2012random}. Indeed, a random exploration of the space allows to evaluate more different values for each of the hyperparameters. This has the advantage of increasing the chances of finding a better configuration, but also to highlight the importance of some hyperparameters compared to the others. In addition, the random search makes it possible to highlight these properties with fewer evaluations than an exhaustive grid search. More recently, the {\sf Hyperband} algorithm~\cite{li2017hyperband} was introduced, which is a variant of a random search that uses an early stopping criteria to detect a non promising point early on in order to save computational resources and time. Thus achieving an important speedup compared to other methods. However, despite its advantages over the grid search, a random approach is limited because it is not adaptive and it does not exploit the performance scores of each configuration to direct the search. This can also waste resources that could have been better exploited by another optimization approach.

Genetic algorithms are evolutionary heuristics that are also used for the HPO problem. Inspired by biology, a genetic algorithm generates an initial population, i.e. a set of configurations, then, it combines the best parents to create a new generation of children. It also introduces random mutations to ensure a certain diversity in the population. These heuristics are therefore adaptive, thus allowing to explore the space more wisely even if they remain impregnated with a random character. These algorithms are often used to optimize hyperparameters~\cite{elsken2018efficient, suganuma2017genetic, young2015optimizing}. In~\cite {lorenzo2017particle}, a method based on particle swarm optimization is able to provide networks with higher performance than those defined by experts in less time than what would have required a grid search or a completely random search. Another approach using the evolutionary algorithm CMA-ES~\cite{loshchilov2016cma} was proposed with satisfactory results.

Other approaches based on machine learning can be found in the literature. For example, the HPO problem can be seen as reinforcement learning~\cite{baker2016designing, zoph2016neural, zoph2017learning} where the main difference between each method relies on how the agents are defined and dealt with. In~\cite{smithson2016neural}, a neural network is able to design other neural networks by learning to explore the possible configurations. Here, the HPO of neural networks is seen as a multiobjective problem where one seeks to improve the performance of the network while minimizing the computing power required. This approach, although successful, solves a different problem from the one considered in the context of this study. Also,~\cite{bosc2016} uses a network of long-term memory neurones to learn the parameters of another multi-layer network that is tested on a binary classification problem.

Derivative-Free Optimization (DFO) is naturally adapted to the HPO problem since it aims at solving problems typically given in the form of blackboxes that can be computationally costly to evaluate, with nonexistent or inoperable derivatives. In~\cite{AuOr06a}, the authors propose a general way of modeling hyperparameter optimization problems as a blackbox optimization problem. This formulation is used in~\cite {Liu2018}  to optimize 11 hyperparameters (3 real and 8  integer) of the BARON solver. This study compared the solutions found by  27 DFO algorithms on a total of 126  problems. The results show that the DFO methods have reduced the average resolution time, sometimes by more than 50\%. Another formulation inspired by robust optimization is used in~\cite{PoTo2017}, in addition to that of~\cite{AuOr06a}, to optimize the hyperparameters of the BFO algorithm~\cite{PoTo2017}. BBO methods are also at the heart of {\sf Google Vizier}\xspace~\cite{golovin2017google} which is a tool that can be used for the HPO problem of machine learning algorithms, and especially for deep neural networks.

Bayesian optimization (BO) can be seen as a subclass of DFO methods and as such, can be used to solve the HPO problem. The BO methods use informations collected during previous assessments to diagnose the search space and predict which areas to explore first. Among them, Gaussian processes (GP) are models that seek to explain the collected observations that supposedly come from a stochastic function. GPs are a generalization of multi-variate Gaussian distributions, defined by a mean and a covariance function. GPs are popular models for optimizing the hyperparameters of neural networks~\cite {Bayesian_HPO, wistuba2018scalable}. However, the disadvantage of GPs is that they do not fit well to the categorical features, and its performance depends on the choice of the kernel function that defines it.  Tree-structured Parzen Estimator (TPE) is also a Bayesian method that can be used as a model instead of a GP. After a certain number of evaluations, this method separates the evaluated points into two sections: a portion  ($<$25\%) of the points with the best performances, and the remaining. This method seeks to find a distribution of the best observations to determine the next candidates.TPEs are also used for the HPO of neural networks~\cite{bergstra2011algorithms}, even if it has the disadvantage of ignoring the interactions between the hyperparameters.

Other model-based DFO methods were also applied to the HPO problem. In~\cite{Nannicini_HPO}, the authors applied radial basis functions to model the blackbox as previously defined. This article presents the results obtained on the MNIST data set~\cite{mnist}, then on a problem of interactions between drugs. These tests show that this model provides comparable or better results than popular configurations. In~\cite{ghsc2017}, a trust-region DFO algorithm is applied to optimize the hyperparameters of a SVM model. Here again, this approach obtained a more efficient model than those defined by the experts or by a Bayesian algorithm.

Thus, the positive results of these methods suggest that the DFO approach is well suited to solve the HPO problem of deep neural networks. This motivated the idea of using the {MADS}\xspace algorithm~\cite{AuDe2006}, which is implemented into the {\sf NOMAD}\xspace package~\cite{Le09b}, especially as it can handle integer and categorical variables~\cite{AACW09a, AuLeDTr2018}.
Using the {MADS}\xspace algorithm for hyperparameter tuning
has been  validated in~\cite{mello2019novel}
where a SVM model is calibrated using MADS combined with
the Nelder-Mead and VNS search strategies~\cite{AuBeLe08,AuTr2018}.

A non exhaustive list of open source librairies for HPO is given in Table~\ref{tab:lib_hpo} along with the optimization algorithms implemented in each library and the types of variables handled.

\begin{table}[ht]
  \caption{Selection of open source libraries for the hyperparameter optimization problem.}
    \centering
   \begin{tabular}{l c c c c c c c r} \toprule
   	{  } & \multicolumn{5}{c }{Optimization method} & \multicolumn{3}{c }{Type of variables} \\ \midrule
	 Package & \makecell{Grid \\ search}  &  \makecell{Random\\ search} &  \makecell{Bayesian \\ optimization} &  \makecell{Model- \\ based} &  \makecell{Direct- \\ search} & Real & Int. & Cat. \\ \midrule
	 {\sf scikit-learn}~\cite{scikit-learn} &  \cmark  &  \cmark  &   - & - &  -  & \cmark & \cmark &  \cmark  \\
	{\sf hyperopt}\xspace~\cite{bergstra2013making} &   -   & \cmark & \cmark &  -  & -   & \cmark & \cmark & \cmark \\
	{\sf Spearmint}~\cite{Bayesian_HPO} & \cmark & \cmark & \cmark & -   &   -  & \cmark & \cmark & \cmark \\
	{\sf SMAC}~\cite{hutter2011sequential} &  -  &  -  &  -  & \cmark &  -  & \cmark & \cmark & \cmark \\
	{\sf MOE}~\cite{moe}&  -  &   - & \cmark &  -  &  -   & \cmark &  -  &  -  \\
	{\sf RBFOpt }~\cite{Nannicini_HPO} &  -  &  -  &   - & \cmark &  -  & \cmark & \cmark &  -  \\
	{\sf DeepHyper}~\cite{deephyper} &  -  &  \cmark  &   - & \cmark &  -  & \cmark & \cmark &  \cmark  \\
	{\sf Or\'ion}~\cite{orion} &  -  &  \cmark  &   - & - &  -  & \cmark & \cmark &  \cmark  \\
	{\sf Google Vizier}\xspace~\cite{golovin2017google} & \cmark & \cmark & \cmark & \cmark & - & \cmark & \cmark & \cmark \\
	{\sf HyperNOMAD}\xspace  & -   &  -  &   - &  -  & \cmark & \cmark & \cmark & \cmark \\  \bottomrule
   \end{tabular}
   \label{tab:lib_hpo}
\end{table}

\section{Experimental setup} \label{sec:exp}

This section first defines the blackbox approach used for modeling the HPO problem. This is done by listing the different hyperparameters considered to construct, train and validate a deep neural network (DNN) in order to obtain its test accuracy. The second part of the section gives an overview of the data sets provided with {\sf HyperNOMAD}\xspace.

\subsection{Hyperparameters of the framework}
\label{sec:hpos}

A variety of hyperparameters must be chosen to tune a DNN for a given application. These hyperparameters affect different aspects of the network: the architecture, the optimization process and the handling of the data. The following section lists the hyperparameters considered in this study along with their respective types and scopes.

\subsubsection{The network architecture}

A convolutional neural network (CNN) is a deep neural network consisting of a succession of convolutional layers followed by fully connected layers as illustrated in Figure~\ref{fig:convnet}.

\begin{figure}[ht]
   \centering
   \begin{tabular}{@{}c@{\hspace{.5cm}}c@{}}
      \includegraphics[scale=0.65]{convnet_fig.png} \\
   \end{tabular}
 \caption{Example of a convolutional neural network. Image taken from~\cite{schema_cnn}.}
 \label{fig:convnet}
\end{figure}

To define a new CNN,  one must first decide on the number of convolutional layers. These layers can be seen as matrices in a two dimensional convolution. The size of the first convolutional layer is determined by the size of the images the network is fed. The size of the remaining layers is computed by taking into account the different operations applied from layer to layer. These operations can be divided into two categories: a convolution or a pooling. Figure~\ref{fig:1a} represents the steps of a convolution operation. The initial image is a 5$\times$5 matrix whose borders are padded with zeros. The convolution consists of choosing a kernel that is passed over the image in order to compute the coefficients of the feature map. Each coefficient is equal to the sum of the products between the coefficients of the image and the ones of the kernel situated in the same position. In Figure~\ref{fig:1a}, the coefficient (2,1) of the feature map is obtained by the following operation: $(0\times0) + (0\times0) + (0\times1) + (0\times0) + (21\times1) + (0\times0) + (85\times1) + (71\times0) + (0\times0) = 106$.
In general, a convolution can be determined with few factors such as the number of feature maps - or output channels - generated, the size of the kernel which in turn will affect the size of the feature map, the stride which corresponds to the step by which the kernel is moved over the image and the padding. In Figure~\ref{fig:1a}, the image is padded with one layer of zeros.

When the feature map is obtained, one can decide to apply a pooling operation to decrease the size of the output by keeping only the biggest coefficients in a certain area. Figure~\ref{fig:1b} illustrates a 2$\times$2 pooling that results in an output of half the size of the feature map.

\begin{figure}[ht]
   \centering
   \begin{subfigure}{0.52\textwidth}
	\includegraphics[width=\linewidth]{convolution.png}
	\caption{Convolution} \label{fig:1a}
   \end{subfigure}
   \hspace*{\fill}
   \begin{subfigure}{0.45\textwidth}
	\includegraphics[width=\linewidth]{maxpooling.png}
	\caption{Pooling} \label{fig:1b}
   \end{subfigure}
 \caption{Illustration of a convolution operation in (a) and a pooling operation in (b). Images taken from~\cite{schema_cnn2}.}
 \label{fig:convolution}
\end{figure}

Each fully connected layer that follows the convolutional ones is determined by the number of neurones it contains.  The neurones of a layer are connected to all of the ones in the next layer through weighted arcs. Let $x_1, x_2, \ldots, x_{n_l}$ be the values of the neurones of the layer $l$ and $a_j$ be the value of the neurone $j$ in the layer $l+1$, then $ a_j = \phi(\sum\limits_{i =1}^{s_l} w_{ij} x_i)$, where $w_{1j}, w_{2j}, \ldots, w_{{n_l}j}$ are the weights of the arcs from the neurones of the layer $l$ to the $j-th$ neurone of the following layer and $\phi$ is an activation function used to introduce a non linearity in the outputs. Figure~\ref{fig:activation} presents some examples of activation functions.

\begin{figure}[ht]
   \centering
    \begin{subfigure}{0.32\textwidth}
	\includegraphics[width=\linewidth]{activation_relu.pdf}

   \end{subfigure}
   \hspace*{\fill}
   \begin{subfigure}{0.32\textwidth}
	\includegraphics[width=\linewidth]{activation_sigmoid.pdf}
   \end{subfigure}
    \begin{subfigure}{0.32\textwidth}
	\includegraphics[width=\linewidth]{activation_tanh.pdf}
   \end{subfigure}
 \caption{Examples of activation functions.}
 \label{fig:activation}
 \end{figure}

Table~\ref{tab:HPNN} summarizes the hyperparameters responsible for defining the structure of the network.
Hyperparameters 2 to 6 must be defined for each convolutional layer and the hyperparameter 8 must also be defined for each fully connected layer. Therefore, if $n_1$ is the number of convolutional layers and $n_2$ the number of fully connected layers, the total number of hyperparameters responsible for defining the structure of the neural network is $5n_1 + n_2 + 4$.

\begin{table}[ht]
  \caption{Hyperparameters that define the architecture of a neural network.}
    \centering

   \begin{tabular}{p{0.08 \textwidth} p{0.45 \textwidth} p{0.2 \textwidth} p{0.22 \textwidth}  }  \toprule

   	 {\#} & \textbf{Hyperparameter} & \textbf{Type} &  \hspace*{1cm}\textbf{Scope} \\
	\hline
	1 & Number of convolutional layers $(n_1)$& Categorical &  \hspace*{1cm}$\{$0,1,\ldots,20$\}$\\
	2 & \hspace*{5mm}Number of output channels & Integer &  \hspace*{1cm}$\{$0,1,\ldots,50$\}$\\
	3 & \hspace*{5mm}Kernel size & Integer &  \hspace*{1cm}$\{$0,1,\ldots,10$\}$\\
	4 & \hspace*{5mm}Stride & Integer &  \hspace*{1cm}$\{$1,2,3$\}$ \\
	5 & \hspace*{5mm}Padding & Integer &  \hspace*{1cm}$\{$0,1,2$\}$\\
	6 & \hspace*{5mm}Do a pooling & Boolean &  \hspace*{.9cm} $\{$0,1$\}$ \\

	7 & Number of full layers $(n_2)$& Categorical &  \hspace*{1cm}$\{$0,1,\ldots,30$\}$\\
	8 & \hspace*{5mm}Size of the full layer & Integer &  \hspace*{1cm}$\{$0,1,\ldots,500$\}$\\

	9 & Dropout rate & Real &  \hspace*{1cm}[0;1] \\

	10 & Activation function & Categorical/Integer &  \makecell{$\{$ReLU (1), \\ Sigmoid (2), \\ Tanh (3)$\}$}  \\
	\bottomrule
   \end{tabular}
   \label{tab:HPNN}
\end{table}

\subsubsection{The optimizer}

For a given network architecture, the training phase is conducted to minimize
the error between the predictions of the network and the correct values of the labels attached to the validation data. Let $\Theta$ be a multi-dimensional matrix that stores the weights of the arcs that link each layer of the network $l$ with the next one $l+1$, and let $J(\Theta)$ the sum of the errors between the predictions and the labels for all the data. The optimizer must then solve $\min_{\Theta } J(\Theta).$  Before starting the training phase, the optimizer that carries out this task must be selected along with its specific algorithmic hyperparameters.

A stochastic gradient approach is more suitable in this case because of the high dimension of this problem which is usually in the millions. At each iteration, the weights of the network are updated by following a stochastic direction with a particular step size which is called a learning rate in the machine learning context. Similarly to any gradient descent method, the learning rate must be chosen and updated accordingly to avoid oscillations or divergence. Substantial research and tricks of the trade are developed to solve this problematic~\cite{bengio2012practical, stochastic-gradient-tricks, LeCun2012}. The optimizers Adam~\cite{Adam}, Adagrad~\cite{duchi2011adaptive} and RMSProp~\cite{rmsprop} have embedded strategies to adapt the learning rate at each iteration and for each weight. SGD however does not require external management. In {\sf HyperNOMAD}\xspace, the learning rate of SGD is divided by 10 every 100 epochs as long as its value is greater than $10^{-6}$.

Table~\ref{tab:optimizer} presents the list of selectable optimizers considered in the blackbox along with their corresponding hyperparameters. There is one categorical hyperparameter that determines which optimizer is chosen and always four real hyperparameters related to it. This aspect of the network relies on defining five hyperparameters in total.

\begin{table}[ht]
\caption{Choices of the optimizer and the corresponding hyperparameters.}
  \centering
   \begin{tabular}{p{0.5 \textwidth}  p{0.25 \textwidth}  p{0.1 \textwidth}  p{0.1 \textwidth}} \toprule
   	\textbf{Optimizer } & \textbf{Hyperparameter} & \textbf{Type} & \textbf{Scope} \\ \midrule
	Stochastic Gradient Descent (SGD)	& Initial learning rate 		& Real 	& [0;1]\\
	 						       	& Momentum 			& Real 	& [0;1]\\
							       	& Dampening 			& Real 	& [0;1]\\
							       	& Weight decay 		& Real 	& [0;1]\\ \midrule
	Adam						& Initial learning rate 		& Real 	& [0;1]\\
								& $\beta_1$			& Real 	& [0;1]\\
								& $\beta_2$			& Real 	& [0;1]\\
								& Weight decay 		& Real 	& [0;1]\\ \midrule
	Adagrad						& Initial learning rate 		& Real 	& [0;1]\\
								& Learning rate decay	& Real 	& [0;1]\\
								& Initial accumulator		& Real 	& [0;1]\\
								& Weight decay 		& Real 	& [0;1]\\ \midrule
	RMSProp						& Initial learning rate 		& Real 	& [0;1]\\
								& Momentum 			& Real 	& [0;1]\\
								& $\alpha$			& Real 	& [0;1]\\
								& Weight decay 		& Real 	& [0;1]\\
	 \bottomrule
   \end{tabular}
   \label{tab:optimizer}
\end{table}

\subsubsection{The training phase}
Before training a network, the data must be separated into three groups, each one is responsible for the training, the validation and the testing of the network. During the training phase, the network is fed with the training data, performs a forward pass, computes the prediction error and do a back-propagation in order to update the weights using the optimizer. The way the network is fed is also of great importance. One can choose to input the training data one by one, all at once, or by sending subsets or mini-batches of the data. The size of the mini-batches is an integer hyperparameter that varies between $[1, n_{train}]$, where $n_{train}$ is the size of the training data.

When the network has been fed all of the training data, it is said to have performed an epoch. Usually, the training data has to be passed more than once in order to obtain good weights and a good testing accuracy. Therefore, the number of epochs must be chosen as well. This hyperparameter is dealt with as follows. The validation accuracy is evaluated after each epoch and the weights of the network responsible for the best validation accuracy are stored. This process is repeated as long as the number of epochs is lower than a certain maximum number of epochs (usually 500) and as long as an early stopping condition has not been satisfied. These early stopping criteria depend on the evolution of the training and validation of the network. When the validation accuracy staggers or when it stays lower than 20\% after 50 epochs then the training can be interrupted in order to save time and computational ressources. Once the training is done, the test accuracy is evaluated using the weights that gave the best validation accuracy.

Finally, the blackbox optimization problem is obtained following the model in~\cite{AuOr06a}. This blackbox takes $5n_1 + n_2 + 10$ mixed variable inputs, where $n_1$ is the number of convolutional layers and $n_2$ the number of fully connected layers of the network, and returns the value of the accuracy on the test data set. This blackbox problem is solved using the {\sf NOMAD}\xspace software~\cite{Le09b} described in Section~\ref{sec:nomad}.

\subsection{Data sets}

The {\sf HyperNOMAD}\xspace package comes with a selection of data sets all meant for classification problems. Table~\ref{tab:datasets} lists the data sets embedded so far through {\sf PyTorch}\xspace~\cite{paszke2017automatic}, a relatively complete tool to model and manipulate deep neural networks. {\sf HyperNOMAD}\xspace also allows the usage of a personal data set by following the instructions given in Appendix~\ref{sec-appendix}.
When loading a data set from Table~\ref{tab:datasets}, {\sf HyperNOMAD}\xspace applies a normalization and a random horizontal flip to regulate and augment the data.

\begin{table}[ht]
\caption{Data sets embedded in {\sf HyperNOMAD}\xspace.}
  \centering
   \begin{tabular}{l r r r r} \toprule
   \textbf{Data set } & \textbf{\makecell{Training \\ data}} & \textbf{ \makecell{Validation\\ data}} &\textbf{\makecell{Testing \\data}} & \textbf{\makecell{Number \\ of classes}} \\ \midrule
   MNIST 	   			&  \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000   & \hspace*{1.25cm}10\\
   Fashion-MNIST		& \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000	   & \hspace*{1.25cm}10\\
   EMNIST 	   		&  \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000	   & \hspace*{1.25cm}10\\
   KMNIST 	   		&  \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000	   & \hspace*{1.25cm}10\\
   CIFAR-10	   		&  \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000	   & \hspace*{1.25cm}10\\
   CIFAR-100 	   		&  \hspace*{1cm}40,000 	     & \hspace*{1cm}10,000	          & \hspace*{1cm}10,000	   &\hspace*{1.25cm}100\\
   STL-10 	   			&  \hspace*{1cm}4,000 	     & \hspace*{1cm}1,000	          & \hspace*{1cm}8,000	   &\hspace*{1.25cm}10\\
   \bottomrule
   \end{tabular}
   \label{tab:datasets}
\end{table}

The rest of the section describes the data sets used for benchmarking {\sf HyperNOMAD}\xspace. First, a validation is done using MNIST~\cite{mnist} and once positive results are obtained, the second and more complex data set, CIFAR-10~\cite{krizhevsky2009learning}, is considered.

\subsubsection{MNIST}\label{sec:dataset_mnist}

MNIST~\cite{mnist} is a data set containing 60,000 images of hand written digits that is usually divided into three categories: 40,000 for training, 10,000 for validation and the remaining 10,000 for testing. The set is used for developing a convolutional neural network capable of recognizing the digits in each image and assigning it to the correct class.
The relative simplicity of this task does not require complex neural networks to obtain a good accuracy. Therefore, this data set is usually considered as a first validation of a concept and not a sufficient proof of the quality of a method among the machine learning community.

\subsubsection{CIFAR-10}

The second set of tests are performed with CIFAR-10~\cite{krizhevsky2009learning}. This data set contains 60,000 colored images of objects that belong to ten different and independent categories. The data is once again divided into three sets: 40,000 for training, 10,000 for validation and 10,000 for testing.

For this test, the blackbox within {\sf HyperNOMAD}\xspace is used to construct the convolutional neural network corresponding to the values of the hyperparameters described in Section~\ref{sec:hpos}. This network is trained, validated and tested on CIFAR-10 according to the mode of operation of {\sf HyperNOMAD}\xspace detailed in Section~\ref{sec:hypernomad}.

\section[HyperNOMAD]{{\sf HyperNOMAD}\xspace}\label{sec:hypernomad}

The  {\sf HyperNOMAD}\xspace package is available on GitHub\footnote{\url{https://github.com/DouniaLakhmiri/HyperNOMAD}}.
It contains a series of {\sf Python} programs wrapped into a blackbox responsible for constructing, training and evaluating the test accuracy of a neural network depending on the values of the hyperparameters described in Section~\ref{sec:exp}. This blackbox uses the {\sf PyTorch}\xspace package~\cite{paszke2017automatic}  for its simplicity. {\sf HyperNOMAD}\xspace also contains an interface that runs the optimization of the blackbox using the {\sf NOMAD}\xspace software~\cite{Le09b} described in the rest of this section. The basic usage of  {\sf HyperNOMAD}\xspace is described in Appendix~\ref{sec-appendix}.

\subsection{Overview of {\sf NOMAD}\xspace}\label{sec:nomad}

The {\sf NOMAD}\xspace software~\cite{Le09b} is a {\sf C++} implementation of the {MADS}\xspace algorithm~\cite{AuDe2006,AuLeDTr2018} which is a direct search method that generates, at each iteration $k$, a set of points on the {\em mesh}
$M^k = \left\{ x +  \mathop{\mathrm{diag}}({\boldsymbol \delta}^k) z : x \in V^k, z \in {\mathbb{Z}}^n \right\}$
where $V^k$ contains the points that were previously evaluated (including the current iterate $x^k$)
and  ${\boldsymbol \delta}^k \in {\mathbb{R}}^n$ is the {\em mesh size vector}.

Each iteration of {MADS}\xspace is divided into two steps: The \textit{search} and the \textit{poll}. The \textit{search} phase is optional and can contain different strategies to explore a wider space in order to generate a finite number of possible mesh candidates. This step can be based on surrogate functions, latin hyper-cube sampling, etc. The \textit{poll}, on the other hand, is strictly defined since the convergence theory of {MADS}\xspace relies solely on this phase. Here, the algorithm generates directions around the current iterate $x^k$ to search for candidates locally in a region centered around  $x^k$ and of radius, in each dimension, of ${\boldsymbol \Delta}^k \in {\mathbb{R}}^n$, which is called the {\em poll size vector}. The set of candidates in this phase defines the {\em poll set} $P_k$.

If {MADS}\xspace finds a better point then the iteration is declared a success and the mesh and poll sizes are increased, however, if the iteration fails then both parameters are reduced so that
${\boldsymbol \delta}^k \leq {\boldsymbol \Delta}^k$ is maintained.  This relation insures that the set of search directions becomes dense in the unit sphere asymptotically.
The {MADS}\xspace algorithm is summarized in Algorithm~\ref{algo-mads}.

\begin{algorithm}
\caption{Mesh adaptive direct search ({MADS}\xspace)}
\label{algo:mads}
\begin{flushleft}
$k = 0, {\boldsymbol \delta}^0$, $x^0$

\textbf{[1] Search (optional)}\\
\hspace*{10mm} Construct a set of mesh points and evaluate them \\
\hspace*{10mm} If there is a success, go to \textbf{[3]}

\textbf{[2] Poll}\\
\hspace*{10mm} Evaluate the points in the poll set $P_k$

\textbf{[3] Updates}\\
\hspace*{10mm} Update ${\boldsymbol \delta}^k, x^k, M^k, V^k$ depending on the success of the previous phases\\
\hspace*{10mm} If no stopping condition is satisfied: $k \leftarrow k+1$ and go to \textbf{[1]}
\end{flushleft}
\label{algo-mads}
\end{algorithm}

In addition, {\sf NOMAD}\xspace can handle categorical variables by adding a step in the basic {MADS}\xspace algorithm.
A variable is categorical when it can take a finite number of nominal or numerical values that express a qualitative property that assign the variable to a class (or category). The algorithm relies on an ad~hoc neighborhood structure, provided in practice by the user as a list of neighbors  for any given point. The poll step of {MADS}\xspace is augmented with the so-called \textit{extended poll} that links the current iterate $x^k$ with the independent search spaces where the neighbors can be found. The first neighbor that improves the objective function is chosen and the optimization carries on in the corresponding search space.
For more detail on how {MADS}\xspace handles categorical variables, the reader is referred to the following list of articles~\cite{Abra04,AACW09a,AbAuDe2007a,AuDe01a,KoAuDe01a}.

\subsection[Hyperparameters in HyperNOMAD]{Hyperparameters in {\sf HyperNOMAD}\xspace}

The selected neighborhood structure in {\sf HyperNOMAD}\xspace relies on blocks of categorical variables with their associated variables.
The following subsections describe this structure.

\subsubsection{Blocks of hyperparameters}

 {\sf HyperNOMAD}\xspace splits the hyperparameters (HPs) defined in Section~\ref{sec:hpos} into different blocks: one for the convolution layers, the fully connected layers, the optimizer and one for each of the other HPs. A block is an implemented structure that stores a list of values, each one starting with a header and followed by the associated variables, when applicable, that are gathered into groups.
 %
For example, consider a CNN with two convolutional layers, each one defined with the number of output channels, the kernel size, the stride, the padding and whether a pooling is applied or not as stated in Table~\ref{tab:HPNN}. Then consider the values $(16, 5, 1, 1, 0)$ and $(7, 3, 1, 1, 1)$. Each set of values corresponds to a group of variables that describes one convolutional layer and both groups constitute the convolution block. The header of the convolution block is the categorical variable that represents the number of convolutional layers ($n_1$) that the CNN contains as shown in Figure~\ref{fig:conv_block} (top).
The convolution block is followed by the fully connected block. The header of this block also corresponds to the categorical variable that describes the number of fully connected layers. Here, each layer is defined with the number of neurones it contains. Therefore, if $n_2$ is the value in the header, then  there are $n_2$ groups of a single variable as illustrated in Figure~\ref{fig:full_block} (top).
The optimizer block always possesses a fixed size since there is always five HPs that describe the optimizer: The choice of the algorithm and four related HPs as summarized in Table~\ref{tab:optimizer}. The header of this block is the categorical variable corresponding to the choice of the optimizer and the four associated variables are gathered into one group as shown in Figure~\ref{fig:opt_block}.
The other HPs are put as the headers of their individual block with no associated variable.

\begin{figure}[!htb]
  \flushleft
\hspace*{3.25em}\begin{tabular}{| p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} |}
   	\hline
   	\cellcolor{green!25}{2 } & \cellcolor{blue!25}{16} & \cellcolor{blue!25}{5} & \cellcolor{blue!25} {1} & \cellcolor{blue!25}{1} & \cellcolor{blue!25}{0} & \cellcolor{yellow!25}{7} & \cellcolor{yellow!25}{3} & \cellcolor{yellow!25} {1} & \cellcolor{yellow!25}{1} &\cellcolor{yellow!25} {1} \\
	\hline
   \end{tabular}
   \hspace*{3em}
   \begin{tabular}{| p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} |p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} |}
   	\hline
   	\cellcolor{green!25}{3 } & \cellcolor{blue!25}{16} & \cellcolor{blue!25}{5} & \cellcolor{blue!25} {1} & \cellcolor{blue!25}{1} & \cellcolor{blue!25}{0} & \cellcolor{yellow!25}{7} & \cellcolor{yellow!25}{3} & \cellcolor{yellow!25} {1} & \cellcolor{yellow!25}{1} &\cellcolor{yellow!25} {1} & \cellcolor{red!25}{7} & \cellcolor{red!25}{3} & \cellcolor{red!25} {1} & \cellcolor{red!25}{1} &\cellcolor{red!25} {1} \\
	\hline
   \end{tabular}
   \hspace*{3em}
   \begin{tabular}{| p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} | p{5mm} |}
   	\hline
   	\cellcolor{green!25}{1} & \cellcolor{blue!25}{16} & \cellcolor{blue!25}{5} & \cellcolor{blue!25} {1} & \cellcolor{blue!25}{1} & \cellcolor{blue!25}{0}  \\
	\hline
\end{tabular}
      \caption{Example of a convolution block (top). Its first neighbor is obtained by adding a convolutional layer (middle) and the second neighbor is obtained by subtracting a convolutional layer (bottom).}
\label{fig:conv_block}
   \end{figure}

\begin{figure}[!htb]
       \flushleft
\hspace{8em}\begin{tabular}{| p{5mm} | p{10mm} | p{10mm} | p{10mm} | }
   		\hline
   		\cellcolor{green!25}{3 } & \cellcolor{blue!25}{1200} & \cellcolor{yellow!25}{512} & \cellcolor{red!25} {20}\\
		\hline
   	\end{tabular}

\hspace{8em}\begin{tabular}{| p{5mm} | p{10mm} | p{10mm} | p{10mm} | p{10mm} | }
   		\hline
   		\cellcolor{green!25}{4 } & \cellcolor{gray!25}{1200}  & \cellcolor{blue!25}{1200} & \cellcolor{yellow!25}{512} & \cellcolor{red!25} {20}\\
		\hline
   	\end{tabular}

\hspace{8em}\begin{tabular}{| p{5mm} | p{10mm} | p{10mm} | }
   		\hline
   		\cellcolor{green!25}{2 }& \cellcolor{yellow!25}{512} & \cellcolor{red!25} {20}\\
		\hline
   	\end{tabular}
		\caption{Example of a fully connected block (top). Its first neighbor is obtained by adding a fully connected layer (middle) and the second neighbor is obtained by subtracting a fully connected layer (bottom).}
		\label{fig:full_block}
\end{figure}

\begin{figure}[!htb]
      \centering
        \begin{tabular}{| p{5mm} | p{8mm} | p{8mm} | p{8mm} |  p{8mm} | }
   		\hline
   		\cellcolor{green!25}{1 } &  \cellcolor{blue!25}{0.2} & \cellcolor{blue!25}{0.95} & \cellcolor{blue!25} {$1e^{-4}$
		} & \cellcolor{blue!25} {0.03}\\
		\hline
 	  \end{tabular}

	  \begin{tabular}{| p{5mm} | p{8mm} | p{8mm} | p{8mm} |  p{8mm} | }
   		\hline
   		\cellcolor{green!25}{2 } & \cellcolor{blue!25}{0.1} & \cellcolor{blue!25}{0.9} & \cellcolor{blue!25} {$5e^{-4}$
		} & \cellcolor{blue!25} {0} \\
		\hline
 	  \end{tabular}
	  	   \caption{Example of an optimizer block. Its neighbor is obtained by selecting the next optimizer and by initializing the associated variables to their default values.}
		   \label{fig:opt_block}
\end{figure}

\subsubsection{Neighborhood structure}

The extended poll of {MADS}\xspace with categorical variables constructs one or more neighbor points from any given point, and evaluates them.
There are up to three categorical variables that are exploited using an ad~hoc generation of neighbor points. A neighborhood structure considering coupled effect between the categorical variables may find promising search spaces but it certainly increases the resources needed to carry out the optimization. To limit the number of neighbor points, the neighborhood structure considers each categorical variable independently. Hence, to create a neighbor point related to a given block, all the remaining values are fixed at the current iterate values.
The neighbor structure of the convolution block is obtained by adding and subtracting a group of associated variables at the right of the block.  These operations can only be performed if the resulting size is within the bounds for the variable $n_1$. When adding a group of associated variables, the values of the associated variables are copied from the most right group. Adding or subtracting a group to the convolution block is illustrated in~Figure~\ref{fig:conv_block}.
The neighbor structure of the fully connected block is obtained by adding and subtracting one associated variable at the left of the block. These operations can only be performed if the resulting size is within the bounds for the variable $n_2$. When adding a group, the value of the associated variable is copied and inserted from the most left value (see example in Figure~\ref{fig:full_block}).
The structure of the network varies when adding or subtracting a convolutional or a full layer, and so does the dimension of the HPO problem. Varying the remaining categorical variables has not such effect.
The categorical variable controlling the choice of optimizer has four possible values: SGD, Adam, Adagrad or RMSprop. The choice of optimizer does not change the dimension of the optimization problem but it affects the interpretation of the four associated variables related to the optimizer as illustrated in Table~\ref{tab:optimizer}. A single neighbor point is obtained by looping between the optimizers listed in Table~\ref{tab:optimizer} from top to bottom. For each possible optimizer, there are four associated variables controlling the algorithm with different interpretations. When the optimizer is changed, these variables are reset to their initial values.
In some cases, a variable controlling a category can be well handled as an integer variable by ordering the categories with a predefined order.  This is the case for the variable selecting among the three possible activation functions (see Section~\ref{sec:hpos} and Table~\ref{tab:HPNN}): ReLU, Sigmoid and Tanh, with corresponding values between 1 and 3.
The choice of activation function and the remaining variables are not treated as categorical variables and no neighborhood structure is required.

\section{Computational results}
\label{sec-results}

This section summarizes the results obtained by  {\sf HyperNOMAD}\xspace and compares them to other methods when applied to the MNIST~\cite{mnist} and CIFAR-10~\cite{krizhevsky2009learning} data sets. For both series of tests, all the hyperparameters discussed in Section~\ref{sec:exp} are allowed to vary. However,  the user of the framework can fix some hyperparameters and choose to focus on others as described in Appendix~\ref{sec-appendix}.
All the following tests are allowed a maximum of 100 blackbox evaluations due to time limitations since one call to the blackbox takes, in average, three to four hours.

\subsection{MNIST}

The first tests are performed on the same blackbox provided by the authors of~\cite{Nannicini_HPO} which considers the MNIST data set with the {\sf Caffe} library~\cite{jia2014caffe}.
The {\sf NOMAD}\xspace software is directly used instead of {\sf HyperNOMAD}\xspace,
 in order to compare the different methods of Table~\ref{tab:res_mnist}. The blackbox takes a simplified set of hyperparameters as described in Table~\ref{tab:simpl_hpo}, constructs a convolutional neural network that is trained on the MNIST data set~\cite{mnist},
and finally returns the validation accuracy as a measure of performance.

\begin{table}[ht]
\caption{Hyperparameters considered for the  tests on the MNIST data set with the simplified {\sf Caffe} blackbox.}
  \centering
    \begin{tabular}{p{0.08 \textwidth} p{0.45 \textwidth} p{0.15 \textwidth} p{0.25 \textwidth}  }  \toprule
   	 {\# } & \textbf{Hyperparameter} & \textbf{Type} &  \hspace*{1cm}\textbf{Scope} \\ \midrule
	1 & Number of convolutional layers & Categorical & \hspace*{1cm}$\{$0, 1, 2$\}$\\
	2 & Number of output channels & Integer &  \hspace*{1cm}$\{$1, 2, \ldots, 50$\}$\\
	3 & Number of full layers & Categorical  &  \hspace*{1cm}$\{$0, 1, 2$\}$\\
	4 & Size of the full layer & Integer &  \hspace*{1cm}$\{$1, 2, \ldots, 50$\}$\\
	5 & Learning rate & Real &  \hspace*{1cm}[0;1] \\
	6 & Momentum    & Real 	&  \hspace*{1cm}[0;1] \\
	7 & Weight decay & Real 	&  \hspace*{1cm}[0;1] \\
	8 & Learning decay & Real  &  \hspace*{1cm}[0;1] \\\bottomrule
   \end{tabular}
   \label{tab:simpl_hpo}
\end{table}

The results are obtained by choosing five random seeds and executing the optimization five times for each seed. Table~\ref{tab:res_mnist} presents the results obtained by a random sampling (RS), RBFOpt, and SMAC, that are taken from~\cite{Nannicini_HPO}, and {\sf NOMAD}\xspace. These results show that using {\sf NOMAD}\xspace surpasses all of the other methods in terms of both the validation and the test accuracies.

\begin{table}[ht]
\caption{Results on MNIST with the simplified {\sf Caffe} blackbox.}
  \centering
   \begin{tabular}{p{0.3 \textwidth} p{0.3 \textwidth} p{0.3 \textwidth}}\toprule
   	{Algorithm } & \makecell{Average accuracy on \\  the validation set} & \makecell{Average accuracy \\ on the test set} \\ \midrule
	RS & \hspace*{2cm}94.02 & \hspace*{2cm}89.07 \\
	SMAC & \hspace*{2cm}95.48 & \hspace*{2cm}97.54 \\
	RBFOpt & \hspace*{2cm}95.66 & \hspace*{2cm}97.93\\
	{\sf NOMAD}\xspace  & \hspace*{2cm}\textbf{96.81} & \hspace*{2cm}\textbf{97.98} \\

	\bottomrule
   \end{tabular}
   \label{tab:res_mnist}
\end{table}

The next phase consists to test {\sf HyperNOMAD}\xspace,
this time with {\sf PyTorch}\xspace and its embedded MNIST data set,
and to compare it with other methods such as a random search and a Bayesian method. The {\sf hyperopt}\xspace library~\cite{bergstra2013making} is the one used in this comparison since it contains a random search in addition to TPE, a Bayesian method that relies on Parzen trees~\cite{bergstra2011algorithms}.  The blackbox used for this comparison is the one embedded in {\sf HyperNOMAD}\xspace which allows for a greater flexibility than the previous {\sf Caffe} blackbox since it takes into account all of the hyperparameters described in Section~\ref{sec:hpos}.
The optimization is launched from the same point that corresponds to the default values for the hyperparameters in {\sf HyperNOMAD}\xspace. This initial point contains 22 hyperparameters and obtains a test accuracy of 93.36\%. Figure~\ref{fig:mnist} shows the evolution of {\sf HyperNOMAD}\xspace versus the two variants of {\sf hyperopt}\xspace (RS and TPE) for this test where, after 100 blackbox evaluations,  {\sf HyperNOMAD}\xspace finds the best configuration with a final test accuracy of 99.61\%. The best solution found by TPE obtains a test accuracy of 99.17\% and the random search fails to improve the initial point.

\begin{figure}[H]
   \centering
       \includegraphics[page=1,width=.5\textwidth]{comparison_mnist_default.pdf}
 \caption{Comparison between {\sf HyperNOMAD}\xspace, TPE and RS when launched from the default starting point of {\sf HyperNOMAD}\xspace,
 on the MNIST data set.}
 \label{fig:mnist}
\end{figure}

\subsection{CIFAR-10}

Similarly to the previous test, {\sf HyperNOMAD}\xspace is compared to TPE and the random search. These tests are launched using different starting points, the first being the default values of the hyperparameters in {\sf HyperNOMAD}\xspace with 22 hyperparameters and the second being a network with the VGG-13 architecture. The VGG networks~\cite{Simonyan14c} are very deep convolutional neural networks with small kernels. Figure~\ref{fig:vgg16} illustrates the architecture of the VGG 16 network.

\begin{figure}[ht]
   \centering
       \includegraphics[scale=0.35]{vgg16-neural-network}
 \caption{Architecture of the VGG-16 network. Image taken from~\cite{vgg_architecture}.}
 \label{fig:vgg16}
\end{figure}

Figure~\ref{fig:comp} compares {\sf HyperNOMAD}\xspace, TPE and the random search starting from the default settings of {\sf HyperNOMAD}\xspace which achieve a test accuracy of 28.3\%. Once again, the random search could not bring any improvements to the initial point whereas the best solution of TPE obtains a test accuracy of 64.12\% and {\sf HyperNOMAD}\xspace finds a solution that achieves 77.6\%.

Figure~\ref{fig:HN_vgg13} shows the results of a second test performed using a starting point with a VGG-13 architecture, which corresponds to 62 hyperparameters, that achieves a test accuracy of 90.8\%. In this example neither the random search nor TPE are able to improve on the initial point given. Moreover, they spend all their evaluation budget sampling non feasible architectures. An architecture is infeasible when the size of the image passed through the convolutional layers becomes nil. This behavior can be explained by the sampling strategy of both methods since they tend to change multiple hyperparameters at once thus increasing the probability of obtaining a non feasible architecture.  {\sf HyperNOMAD}\xspace is much more conservative when choosing a new point to evaluate which is why 50 evaluated points are feasible.
The best configuration found by {\sf HyperNOMAD}\xspace achieves a final test accuracy of 92.54\%.

\begin{figure}[ht]
   \centering
      \begin{subfigure}{0.49\textwidth}
	\includegraphics[width=\linewidth]{comparison_cifar10_default_2.pdf}
	\caption{Default starting point.} \label{fig:comp}
   \end{subfigure}
    \begin{subfigure}{0.49\textwidth}
	\includegraphics[width=\linewidth]{comparison_cifar10_vgg16_2.pdf}
	\caption{Starting from a VGG architecture.} \label{fig:HN_vgg13}
   \end{subfigure}
   \hspace*{\fill}

\caption{Comparison between {\sf HyperNOMAD}\xspace, TPE and RS, on the CIFAR-10 data set.}
 \label{fig:num_res_cifar10}
\end{figure}

\section{Discussion}

This work introduces {\sf HyperNOMAD}\xspace, a framework package for hyperparameter optimization of DNNs using the {\sf NOMAD}\xspace software~\cite{Le09b}. The key aspects of this framework is its ability to optimize both the architecture and the optimization phase of a deep neural network simultaneously on the one hand, and to explore different search spaces during a single execution by taking advantage of categorical variables. The framework obtains good results for both the MNIST and CIFAR-10 data sets and finds better solutions than TPE and a random search as is illustrated in Figure~\ref{fig:comp}.
Future work aims at considering different techniques of data augmentation as additional hyperparameters of the blackbox, adding more flexibility in the way the learning rate is updated and expanding the framework to other types of problems than classification and provide interfaces compatible with other tools such as {\sf Tensorflow} or {\sf Caffe2}.

\section*{Acknowledgments}

The authors would like to thank the Nvidia GPU Grant Program for providing a GPU used in this research and Dr. Giacomo Nannicini for providing the initial blackbox used for the preliminary testings of {\sf HyperNOMAD}\xspace.


\section{Introduction}

Kerr frequency combs have revolutionized several fields in optical sciences since they were first proposed in monolithic microresonators \cite{kippenberg0}. To date, they have been successfully applied to a vast variety of state-of-the-art technologies and highly sophisticated measurement techniques as ultra-high resolution spectroscopy \cite{swann,hansch}, massively parallel coherent optical communications \cite{koos1}, light detection and ranging (LIDAR) \cite{koos2}, optical-frequency synthesizers \cite{scott} and microphotonic astrocombs \cite{benedick,herr}. Frequency combs in optical microresonators are generated by the interaction of {either, a continuous-wave (cw) pump laser or pumped optical pulses \cite{herr2},} with the resonator modes via parametric four-wave mixing (FWM), assisted by the Kerr non-linearity of the material. This parametric process fulfills the energy conservation and is enhanced when the side bands created by FWM coincide with the resonator frequencies. Since the intensity-dependent refractive index induces a frequency shift on the modes, an increasing free spectral range (FSR) with frequency (anomalous dispersion) is required to compensate this effect and effectively produce equidistant spectral lines suitable to support the cascaded FWM \cite{kippenberg5}. When a high number of mode-locked modes are involved in this non-linear interaction and dissipation cannot be neglected, the complex non-linear dynamics may give rise to disipative Kerr solitons (DKS), which arise as a double balance between non-linearity and dispersion (preserving their shape), as well as dissipation and parametric gain (preserving their amplitude) \cite{kippenberg4}. {Moreover, multiple-soliton formations have become very interesting in the context of localized-state interactions \cite{tlidi,rivas2,torner}, breather solitons \cite{gaeta1,kippenberg6,kippenberg7} and soliton crystals \cite{cole,karpov}.} DKS are notably relevant because their corresponding spectra exhibit highly coherent frequency combs with perfect single-FSR spacing between side bands \cite{kippenberg1,gaeta5}, and dispersive waves in presence of high-order dispersion \cite{erkintalo,skryabin,rivas,kippenberg2}. The latter is particularly important because it allows a further spanning of the frequency combs, highlighting the important role of dispersion engineering on parametric gain, and consequently, on comb generation \cite{gaeta2}. Such task has been successfully achieved in non-linear microresonators by varying the ring geometrical parameters in order to tailor the dispersion for specific comb functionalities \cite{michel,gaeta3,yu}. Nevertheless, this geometry has a limited parameter space, thus restricting the choice of materials, operation wavelength and size of the final device. Recently, Fabry-Perot resonators have been proposed as appealing candidates for frequency comb generation since different methods to reshape the cavity dispersion can be applied \cite{lugiato,papp}, however, advanced dispersion engineering still remains challenging.

In this {work}, we propose a coupled-cavity waveguide (CCW, also known as coupled-resonator optical waveguide or CROW) as a new candidate to efficiently generate low-threshold Kerr frequency combs and DKS. Similar to ring resonators, CCWs define a discrete spectrum of propagating modes which may trigger cascaded parametric FWM in the anomalous dispersion regime, and lead to a highly coherent comb if the waveguide dispersion is conveniently optimized. In fact, since their appearance two decades ago \cite{yariv}, CCWs have shown to be extremely flexible for advanced dispersion engineering and enhancement of non-linear phenomena via slow-light \cite{fan,blair,krauss,notomi1,notomi2,momchil1}. {We show in Fig.~\ref{fig3}(a) a schematic representation of the CCW studied in this work. The periodic waveguide is formed by a staggered distribution of cavities where first (dashed red) and second (dashed green) neighbor cavity-cavity coupling is assumed. The system is pumped with a cw source and a DKS is propagating along the waveguide direction $y$.}

\begin{figure}[t!]
\centering
\includegraphics[width=0.47\textwidth]{fig3v2.eps}
\caption{{(a) Schematic representation of a CCW system with first and second neighbor cavity coupling. The system is pumped with a cw source and a DKS propagates along the waveguide direction $y$.} (b) Photonic crystal CCW formed by coupled L3 cavities where the red and green holes are allowed to vary in size to tune the first and second neighbor coupling, respectively, while the out-of-plane losses of the waveguide Bloch modes is optimized by varying the position and size of the white holes. The intensity profile of the displacement field [$\mathbf{D}(\mathbf{r})=\epsilon(\mathbf{r})\vec{{\cal E}}(\mathbf{r})$] is shown at the boundary of the Brillouin zone.}\label{fig3}
\end{figure}

\begin{figure*}[t!]
\centering\includegraphics[width=0.97\textwidth]{fig4v3.eps}
\caption{(a) Dispersion relation, (b) second order dispersion, and (c) group index in the Brillouin zone of a photonic crystal CCW of 400 L3 cavities. Red and blue curves are where the dispersion is normal and anomalous, respectively. The waveguide is pumped at $a\alpha_0=1.2742$. (d) Threshold (red curve) and bistability boundaries (blue curves) determined by the external pump intensity $|{\cal F}_{\alpha_0}|^2$ as a function of the laser detuning $\sigma_{\alpha_0}$. Dashed curves are for $\beta_{\rm TPA}=0$. (e) Averaged intra-waveguide power of the CCW in the steady state as a function of $\sigma_{\alpha_0}$ for $|{\cal F}_{\alpha_0}|^2=6|{\cal A}_{\alpha_0}|^2_{\rm th}$. Hysteresis arises in the yellow region {$\sigma_{\alpha_0}<-\gamma_{\alpha_0}\sqrt{3}\rho(\kappa)/2$} with $\rho(\kappa=0.2236)=1.84$, and the inset corresponds to a zoom of the rectangular gray area where the discrete steps, signature of switching between soliton states, appear. The corresponding frequency combs at the marked points I-IV in (e) are respectively shown in (f) supercritical Turing pattern of 40-FSR (355.2~GHz) repetition rate, (g) supercritical Turing pattern of 41-FSR (363.3~GHz) repetition rate, (h) soliton molecule of two pulses with single FSR (9.1~GHz) repetition rate, and (i) soliton pulse with single FSR repetition rate. All power quantities are given in {dB units relatively to the} threshold {$|{\cal A}_{\alpha_0}|_{\rm th}^2$}, while detunings are in $\gamma_{\alpha_0}$ units.}\label{fig4}
\end{figure*}

\section{Non-linear photonic crystal CCW}

{The system of Fig.~\ref{fig3}(a) is realized with} the waveguide of coupled L3 photonic crystal slab cavities shown in Fig.~\ref{fig3}(b), where we have plotted the intensity profile of the displacement field [$\mathbf{D}(\mathbf{r})=\epsilon(\mathbf{r})\vec{{\cal E}}(\mathbf{r})$] at the boundary of the Brillouin zone. The photonic crystal is formed by a hexagonal lattice of holes, pitch $a=400$~nm and hole radii $r=0.25a$, etched in a silicon membrane of thickness $d=0.55a$, { and the L3 cavity is introduced by removing three holes along the $\Gamma K$ direction of the lattice}. This CCW configuration allows us to tune first and second neighbor coupling between the L3 cavities by varying the size of the red and green holes, respectively, and optimize the out-of-plane losses by varying the position and size of the white ones. The waveguide of Fig.~\ref{fig3}(b) has been previously proposed as a compact CCW with outstanding figures of merit obtained via automated global optimization \cite{momchil1} and successfully measured experimentally \cite{mohamed,mohamed2}. Here, we adopt the design with largest averaged group index and small second order dispersion reported in Ref.~\citenum{momchil1} with parameters $(\Delta r_1,\Delta r_2, \Delta r_3, \Delta x)=(-0.0049,-0.0340,-0.1016,0.2204)a$, where $r_1=r+\Delta r_1$, $r_2=r+\Delta r_2$ and $r_3=r+\Delta r_3$ are the radii of the red, green and white holes, respectively, and $\Delta x$ is the outward displacement of the latter ones. The dispersion relation $\nu_\alpha$ of the waveguide and decay rates $\gamma_\alpha$ of the Bloch modes of momentum $\alpha$ are calculated with the guided mode expansion method (GME) \cite{gme} for a system length of $M=400$ cavities (400 normal modes within the Brillouin zone of the CCW) with period $l=\sqrt{3}a$ (waveguide length $L=Ml$), while the non-linear mode volume is computed with a commercial FDTD solver \cite{lumerical} and found to be $V_c=0.345~\mu$m$^3$. {This sample length represents a good compromise between total number of normal modes [avoiding finite-size effects in the coupled mode equations Eq.~(\ref{cme7})] and robustness to fabrication disorder \cite{song,mohamed2}.} Silicon parameters at telecom frequencies are considered for the material, namely, dielectric constant $\epsilon=12.04$, Kerr coeficient $n_2=5.52\times10^{-18}$~m$^2/$W and two-photon absorption (TPA) coefficient $\beta_{\rm TPA}=1\times10^{-11}$~m$/$W \cite{jalali}. The GME dispersion relation, second order dispersion and group index are shown in Fig.~\ref{fig4}(a)-\ref{fig4}(c), respectively, for which red and blue segments correspondingly highlight the regions of normal and anomalous dispersion. {Notice that these dispersion curves are valid for a straight waveguide with periodic boundary conditions, nevertheless, we expect them to describe very well a close loop of coupled cavities as bending losses is usually very small, with respect to out-of-plane losses, in photonic crystal geometries \cite{zhang}}. The system is pumped at $a\alpha_0=1.2742$, where $\nu_{\alpha_0}=193.39$~THz, $(ac)^{-1}(d^2\nu/\alpha^2)_{\alpha_0}=8.63\times10^{-4}$, $n_{g,\alpha_0}=119.34$ and the normal mode quality factor is $Q_{\alpha_0}=5.72\times10^4$. This sets an internal mode threshold given by
\begin{equation}\label{Thm}
\left|{\cal A}_{\alpha_0}\right|_{\rm th}^2=\frac{\epsilon V_c}{2ln_{g,\alpha_0}n_2Q_{\alpha_0}}f(\kappa)=138~\mbox{mW},
\end{equation}

\noindent where ${\cal A}_{\alpha_0}$ is the slowly-varying amplitude of the Bloch mode $\alpha_0$, and $f(\kappa)=(\sqrt{1+\kappa^2}+2\kappa)/(1-3\kappa^2)$, with $\kappa=c\beta_{\rm TPA}/(2n_2\omega_{\alpha_0})$, is a function of the material parameters only which gives $f(\kappa=0.2236)=1.73$ for silicon at $\omega_{\alpha_0}/2\pi=\nu_{\alpha_0}=193.39$~THz, and reduces to $f(\kappa=0)=1$ for $\beta_{\rm TPA}=0$. This result is particularly remarkable because, even by considering that TPA increases the threshold by a factor of 1.73, $|{\cal A}_{\alpha_0}|_{\rm th}^2=138~$mW is still one order of magnitude smaller than in typical mm-size crystalline non-linear ring resonators \cite{chembo2}. Equation~(\ref{Thm}) undoubtedly strengths the potential capabilities of CCW structures for low-threshold comb generation. While the factor $V/Q$ also enters in the expression for the internal mode threshold of ring resonators \cite{chembo0,chembo1,chembo3}, CCWs can drastically decrease this minimal power in the slow-light regime, i.e., for $n_{g,\alpha_0} \gg \sqrt{\epsilon}$. Such structural effect allows to access the rich physics of DKS in CCWs at much lower powers than in their monolithic counterparts even in presence of strong TPA, which is the main source of non-linear losses in silicon structures at telecom wavelengths \cite{gaeta4,agrawal}. It is worth emphasizing that, since we are not specifying any particular coupling architecture with the external pump, we are referring to the internal mode threshold which scales as $1/Q$, and not to the external power threshold which is found to scale as $1/Q^2$ for ring resonators at critical coupling with the external source \cite{kippenberg3}. {While the latter is an important parameter to set the excitation power in an experimental setup, the former quantity, also called the minimum intracavity power for comb generation \cite{chembo2}, gives the corresponding effective power into the specific normal mode of the system.} We plot in Fig.~\ref{fig4}(d) the driven intensity $|{\cal F}_{\alpha_0}|^2$, in $|{\cal A}_{\alpha_0}|_{\rm th}^2$ units, evaluated at threshold (continuous red lines) and at the bistability boundaries (continuous blue curves), as a function of the laser detuning $\sigma_{\alpha_0}=\Omega_0-\omega_{\alpha_0}$ in units of $\gamma_{\alpha_0}$, where the dashed lines are the corresponding curves for $\beta_{\rm TPA}=0$. The effects of TPA on the threshold are clearly seen in this plot, moreover, the hysteresis region (where the system has bistable states) is red-shifted and the area between the bistability boundaries is significantly reduced. Hence, the system needs to be pumped stronger and the laser frequency has to be further decreased to see Kerr frequency combs and DKS in presence of TPA.

\begin{figure}[t!]
\centering\includegraphics[width=0.47\textwidth]{fig5v4.eps}
\caption{Steady state envelope functions at the marked points I-IV in Fig.~\ref{fig4}(e). (a) Supercritical Turing pattern with 40 rolls. (b) Supercritical Turing pattern with 41 rolls. (c) Soliton molecule of two pulses. (d) Single soliton pulse. {$|\psi^s|^2$ is given in units of the power threshold $|{\cal A}_{\alpha_0}|_{\rm th}^2$.}}\label{fig5}.
\end{figure}

\section{Frequency comb and DKS solutions}

In order to find the possible steady state solutions determined by the non-linear CCW dynamics, we carry out a frequency scan at $|{\cal F}_{\alpha_0}|^2=6|{\cal A}_{\alpha_0}|_{\rm th}^2$ in Fig.~\ref{fig4}(d) within the $\sigma_{\alpha_0}/\gamma_{\alpha_0}$ interval $[-3,1]$. The dynamical non-linear coupled-mode equations Eq.~(\ref{cme7}) are thus propagated in time for each $\sigma_{\alpha_0}$ value along this trajectory using an explicit Runge-Kutta integrator and fast Fourier transform (see Ref.~\citenum{hansson}) until the steady state is reached. Our simulation takes into account the momentum-dependent decay rate $\gamma_{\alpha}$ of the Bloch modes, and the initial condition of the integrator is a sharp Gaussian pulse $\psi(y,0)=\exp[-0.5(y/l)^2]$, where $\psi(y,t)=\sum_\alpha{\cal A}_\alpha(t)e^{i\sigma_\alpha t}e^{-i(\alpha-\alpha_0)y}$ is the envelope function along the CCW direction. Results of this analysis are presented in Fig.~\ref{fig4}(e) where we show the averaged waveguide power in the steady state $\langle|\psi^s|^2\rangle=\int_L|\psi^s(y)|^2dy$, in units of threshold, as a function of the laser detuning. When the laser frequency is red-shifted the intra-waveguide power start to increase until a clear series of discrete steps (see inset), decreasing the average power, are seen in the yellow area, which corresponds to the region where there is hysteresis in the system and the soliton solutions are expected to appear \cite{kippenberg4,chembo4}. In fact, these steps have been previously measured experimentally in the transmission of non-linear ring resonators, and they are associated to the formation of different soliton states within the system, in excellent agreement with the predictions of the non-linear coupled mode equations and the Lugiato-Lefever formalism \cite{kippenberg1}. We show in Figs.~\ref{fig4}(f)-\ref{fig4}(i) the respective frequency combs at the four representative points marked in Fig.~\ref{fig4}(e). Figures~\ref{fig4}(f) and \ref{fig4}(g) displays a sequence of primary combs separated by 40-FSR and 41-FSR, respectively, which correspond to supercritical Turing patterns as they are excited above threshold. In Figs.~\ref{fig4}(h)-\ref{fig4}(i), both combs are subcritical (excited right below threshold) with a single FSR spacing, and they are the signature of soliton complexes. Notably, because of the presence of non-trivial high order terms in the dispersion relation of the photonic crystal CCW, these frequency combs {have} signatures of dispersive wave formation or soliton Cherenkov radiation {at the linear phase matching condition $\omega_\alpha-[\omega_{\alpha_0}+\zeta_1(\alpha-\alpha_0)]=0$, where $\zeta_1$ represents the group velocity \cite{kippenberg2}. Particularly, the dispersive peak clossest to the zero comb line is predicted around $a\alpha\simeq1.63$, which is in perfect agreement with our numerical simulations}. The corresponding steady state intensity of the envelope functions are shown in Fig.~\ref{fig5}. For the I (40 FSR repetition rate) and II (41 FSR repetition rate) states, we see 40 and 41 Turing rolls, respectively, in agreement with the Lugiato-Lefever theory of Kerr frequency combs in the anomalous dispersion regime \cite{chembo4}. Moreover, for the III and IV states we clearly identify, respectively, a soliton molecule composed of two pulses and a single soliton pulse, propagating along the waveguide direction while keeping their shape and amplitude. It must be said that because these soliton structures are subcritical (pumped below threshold), they are extremely sensitive to the initial conditions, and therefore, the system may follow different trajectories in Fig~\ref{fig4}(e) when slightly changing the Bloch mode amplitudes at $t=0$ \cite{chembo4,kippenberg1}. {Furthermore, the characteristic width of the soliton pulse is much smaller than the waveguide length thus ruling out finite-size effects in the basis expansion.} {While results presented in Fig.~\ref{fig4} are strictly
valid for a CCW system with periodic boundary conditions, we still expect highly coherent combs in finite waveguides where non-negligible corrections to the non-linear interaction are predicted by the theory of frequency combs in Fabry-Perot resonators \cite{lugiato}. Nevertheless, a rigorous quantitative analysis is out of the scope of this paper, and will be focus of a future work.}

\begin{figure}[t!]
\centering\includegraphics[width=0.47\textwidth]{fig6.eps}
\caption{Frequency spectrum of the DKS in Fig.~\ref{fig5}(d), where the dashed vertical lines represent the normal modes of the photonic crystal CCW. A small repetition rate of 9.1~GHz is obtained.}\label{fig6}
\end{figure}

Figure~\ref{fig6} shows the frequency spectrum of the DKS in Fig.~\ref{fig5}(d), where the unloaded modes (normal modes of the CCW) are represented by the vertical dashed lines. As expected from energy conservation, the spectral lines are equally spaced with a repetition rate of 9.1~GHz. Such rates have been previously achieved in non-linear ring resonators of 7 mm diameter \cite{vahala}, in contrast to our system length of $\sim277~\mu$m. { Notice that the minimal external power for comb generation may significantly vary depending on the material, geometry and more importantly, the coupling architecture with a bus waveguide \cite{chembo2}, however,} this important result notably highlights the potential capabilities of $\mu$m-scale low-threshold CCW systems to generate low-rate frequency combs, which are desirable in high-precision comb applications \cite{kippenberg4}.

{It is important to stress that the specific choice of photonic crystal CCWs is mainly motivated by their characteristic diffraction-limited mode volumes \cite{hughes1} which play an important role on decreasing the threshold power and final device size. Nevertheless, the results presented in this section only rely on the dispersion relation of the CCW and are therefore expected for any CCW system, e.g. coupled microrings, displaying anomalous dispersion with, more generally, $n$-th neighbor coupling.}

\section{Conclusions}

In conclusion, we have presented results on Kerr frequency combs and DKS in CCWs where the main source of non-linear loses is given by TPA. We found that the internal mode threshold for comb generation depends on the main CCW figures of merit, namely, cavity volume, normal mode quality factor and group index. While TPA losses has the main effect of increasing the threshold power and inducing a red shift in the optimal laser detuning, structural slow-light plays an important role on reducing the minimal power to trigger FWM phenomena between the CCW normal modes, evidencing the capabilities of CCW systems for low-threshold frequency combs. Specifically, we have demonstrated the possibility of DKS in a realistic dispersion-engineered silicon photonic crystal CCW at telecom wavelengths, where highly-complex combs were seen with signatures of soliton Cherenkov radiation. Repetition rates of few gigahertz in a $\sim277~\mu$m-length photonic crystal CCW have been obtained into the soliton regime, which are commonly achieved in $\sim7$~mm size ring resonators, thus demonstrating the potential of CCW systems for high-precision comb applications in ultra-compact devices. Finally, although the spectral spanning of the combs generated in CCW is limited to the waveguide bandwidth, our results open the way to Kerr frequency comb and DKS generation via advanced dispersion engineering and slow-light non-linear enhancement.

\onecolumngrid