# An illustrated tutorial on global optimization in nanophotonics

Pauline Bennet  Denis Langevin  Chaymae Essoual  Université Clermont Auvergne,
Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France
Abdourahman Khaireh-Walieh  LAAS, Université de Toulouse, CNRS, Toulouse,
France  Olivier Teytaud  Meta AI Research Paris, France  Peter Wiecha  LAAS,
Université de Toulouse, CNRS, Toulouse, France  Antoine Moreau  Université
Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000
Clermont-Ferrand, France<EMAIL_ADDRESS>
###### Abstract

Numerical optimization for the inverse design of photonic structures is a tool
which is providing increasingly convincing results – even though the wave
nature of problems in photonics makes them particularly complex. In the
meantime, the field of global optimization is rapidly evolving but is prone to
reproducibility problems, making it harder to identify the right algorithms to
use. This paper is thought as a tutorial on global optimization for photonic
problems. We provide a general background on global optimization algorithms
and a rigorous methodology for a physicist interested in using these tools –
especially in the context of inverse design. We suggest algorithms and provide
explanations for their efficiency. We provide codes and examples as an
illustration than can be run online, integrating quick simulation code and
Nevergrad, a state-of-the-art benchmarking library. Finally, we show how
physical intuition can be used to discuss the results of the optimization and
to determine whether the solutions are satisfying or not.

##  I Introduction

While efficient automated design methods for multilayered structures have
emerged in the 1970s, typically, numerical optimization has been used only
more recently, thanks to the increase in the available computational power and
the progress in simulation techniques. These developments lead to methods
providing original and efficient designs for three-dimensional structures, for
which no design rules exist[1]. In photonics, the most promising approaches so
far are inspired by successful methods from mechanics and are based on local
optimization algorithms[2]. However, in photonics, the wave nature of the
problem typically generates a very large number of local minima, making global
optimization algorithms valuable tools, while they are in many cases
unreasonably expensive in mechanics[3].

Numerical optimization is a domain in which significant progress has been made
in the last two decades, with enormous practical implications. Recent results
suggest that modern global optimization algorithms are able to provide us with
original and efficient photonic designs[4] that are particularly convincing as
they can be understood from a physical point of view[5]. However,
reproducibility problems have made the important results in the optimization
field harder to identify and to trust[6] – especially for researchers in other
communities.

The aim of this paper is to serve as a tutorial for photonics specialists who
are interested in using modern numerical optimization tools. We provide
insights, practical tips, and guidance to help researchers navigate the
challenges and pitfalls associated with optimization. We demonstrate how
simulation tools and state-of-the-art optimization libraries can be easily
integrated to effectively tackle inverse design problems.

Specifically, we provide examples of multi-layer photonics problems, simulated
with the PyMoosh toolkit [7] and optimized using the Nevergrad python library
[8]. We present a comprehensive methodology that includes defining relevant
observables, choosing optimization strategies, and computing specific criteria
to assess the reliability of the obtained solutions. We offer practical
examples inspired by real-world problems involving multilayered structures.
These examples effectively illustrate our methodology and make it easy to
transpose to other situations. For easy reproducibility, the codes are
provided as an online Jupyter Notebook, hosted on google’s Colab plateform[9].

In the first part, we provide essential background information on automated
design of photonic structures and optimization. We also introduce the test
cases that we use to demonstrate the concepts, discussed in this paper.

The second part focuses on the seamless integration of the PyMoosh simulation
library [7] and the Nevergrad optimization library [8] using Python. We
explain how these libraries work together to define the cost function and to
enable efficient benchmarking of optimization algorithms. In the final part,
we present general guidelines for optimizing photonic structures. These
guidelines serve as a methodology to avoid common pitfalls and maximize the
potential of the optimization process.

##  II General background and description of the test cases

In the fields of photonics and optimization, numerous concepts and
terminologies have been introduced over the years. As we now require tools and
concepts from both domains to optimize complex photonic structures, it is
important to present an overview of the vocabulary that has been developed and
provide clear definitions for key terms.

###  II.1 Fundamentals of optimization in photonics

Figure 1: Parametrization. The process of selecting the parameters to describe
the structure plays a crucial role in determining the types of solutions that
can be obtained. (a) Less than 10 parameters are optimized. The geometry of
the structure is fully constrained. We refer to problems of such low degrees
of freedom as optimization, but not inverse design. Image from [10]. (b)
Around a hundred parameters are optimized. The geometry is constrained. We
consider problems of this or larger degrees of freedom to be inverse design
because very different patterns can emerge, not necessarily a checkerboard.
Image from [4]. (c) Around a thousand parameters are optimized on a pixel-
based optimization. Each pixel is filled with one material or another to
create the final design. Image from [11] (d) Tens of thousands of parameters
are optimized. We call this type of image parametrization topology
optimization. Image from [12].

Defining the optimization problem. To apply optimization techniques for
improving a photonic structure, the structure must be represented by a finite
set of parameters. This process of selecting the parameters to describe the
structure is known as parametrization, and it plays a crucial role in
determining the types of solutions that can be obtained, potentially
introducing biases. Figure. 1 presents typical examples of problems with
increasing complexity in their parametrization.

In photonics, many problems are continuous, meaning that valid structures can
be achieved by making the parameters change continuously. If this is not the
case, the problem is said to be discrete, requiring specialized algorithms
that will not be discussed here [13, 14]. In the present paper we focus on
continuous optimization.

It is then necessary to define a cost function that quantifies the discrepancy
between the actual optical response of the structure and the desired response.
This cost function serves as a measure of how close the performance of the
structure is to the target, even if this is not a measure in a strict
mathematical sense. Finally, an optimization domain must be defined, typically
by setting bounds on the parameter values. Together, the cost function and the
optimization domain form the optimization problem. The global optimum is the
point of the optimization domain where the value of the cost function is the
lowest. While it is simple to determine whether a solution corresponds to a
minimum of the cost function, called a local minimum, it is generally
impossible to know whether such a minimum is the global optimum, i.e., the
solution with the lowest possible value of the cost function, and hence the
closest to the desired response. The surface representing the cost function as
a function of the parameters is called the cost function landscape.

From optimization to inverse design. Optimization techniques have found
applications in the improvement of photonic structures. They can also be
employed in retrieving experimental parameters, as for instance commonly done
in ellipsometry.

More recently, optimization has been used in cases when the geometric model is
sufficiently unconstrained such that totally new solutions can emerge. Such an
approach is commonly called inverse design. It is difficult to determine
exactly when an optimization problem becomes an inverse design problem. For
instance, optimizing a multilayered structure consisting of 20 layers with
alternating refractive index can already be considered an inverse design
problem, even though the geometry is relatively constrained: multilayers
present a wide diversity of optical response of a 20 layer structure, allowing
to propose very diverse and unexpected solutions to different problems [15].

On the contrary a complex 3D structure characterized by less than 10
parameters is typically not an inverse design problem[16] because the
structure geometry is already fully constrained and no new pattern can emerge
from the optimization process. A rule of thumb would be to consider that an
optimization problem presents the characteristics of an inverse design problem
when more than 20 parameters are used to describe the structure. As explained
below, in order to constrain the structures as little as possible, several
thousand parameters are often considered. Inverse design problems are thus
particularly complex, requiring advanced methods and increased computational
power to solve them effectively.

Local optimization. A local optimization algorithm starts with one structure
(either generated randomly within the optimization domain or given by the
user), and tries to improve it at each step by changing its parameters by a
small amount – often guided by the gradient of the cost function landscape.
Typically the new structure is kept, if the associated cost function is better
than for the previous structure. This makes for a simple and fast optimization
technique, but has the problem of easily getting stuck in local minima of the
cost function. Several strategies are possible to avoid the problem of local
minima[17, 18, 19], but they often make the algorithms more computationally
expensive. In the end, performing multiple iterations of the optimization
process with different initial structures can provide insight into the
presence of local minima, which indicates the level of difficulty of the
problem. Put simply, the greater the number of local minima, the less likely
it is for a solution found by local optimization to be the global optimum.

An attraction basin is the region of the optimization domain for which a
solution will be found whenever a local algorithm starts in this region. An
optimization domain can typically be divided into attraction basins.

Global optimization and Black Box Optimization. The search for a global
optimum is called global optimization (GO). Algorithms that optimize without
using the gradient or other side information, are called Black Box
Optimization (BBO). BBO and GO are not synonymous: the former refers to the
absence of additional information (gradient or other), whereas the latter
refers to caring about global minima rather than local minima. However, there
is a link, because local methods in continuous domains frequently use
gradients. BBO algorithms are generalist in nature and can be applied to a
large variety of problems. A wide range of algorithms exist, and new ones are
continually proposed[20]. This includes genetic algorithms, mathematical
programming, Bayesian optimization, machine learning, reinforcement learning
or other BBO methods. Most of these algorithms are heuristic in nature, making
them non-deterministic. This means that two different runs, even with the same
starting points, can yield different solutions. Moreover, the performance of
an algorithm often vastly depends on the specific optimization problem at
hand[21], making it challenging to compare different algorithms. Consequently,
BBO suffers from poor reproducibility, hindering the identification of the
most efficient algorithms[22]. Combining rigorous benchmarking[23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33] and BBO libraries is the unique approach, able to
address the reproducibility crisis in optimization, ensuring transparency and
reliability of results. A lot of efforts have thus been made to design
benchmarks for BBO[34, 35, 6, 8].

Topology optimization. With the increase in available computational power, it
has become recently possible to divide a given region of space into
pixels/voxels that can typically be either filled with a material or left
empty and to look for a structure presenting a desired optical response. Then,
the number of parameters is intentionally very large, to offer the algorithms
a very large number of degrees of freedom. Such an approach is called Topology
Optimization (TO). TO draws inspiration from its successful application in the
field of mechanics, where it has been widely used. Given the large number of
parameters used for representing a structure, TO usually employs gradient-
based methods. First, the problem is made continuous instead of discrete
(continuous properties are allowed for a given pixel instead of a binary
choice between filled and void) and since the gradient can be computed for a
low computational cost (through the use of the adjoint method [36]), steepest
descent algorithms seem a natural choice. This approach has been extremely
successful in mechanics, to the point that global optimization has been
considered obsolete [3].

When applied to photonics problems, this approach has at first shown
remarkable success, demonstrating that photonic structures can be miniaturized
to an unprecedented degree while maintaining high optical performance[12].
However, the physics fundamentally differ in some aspects between mechanics
and photonics. Mechanical problems can be regarded as comparatively simpler
since a continuous deformation of an object results in a continuous and nearly
monotonic deformation of its properties. The optimization process in photonics
poses greater challenges compared to mechanical cases, primarily due to the
wave nature of the problem. Photonic structures can exhibit resonant behavior,
with each resonance corresponding to a local maximum or minimum of the cost
function. This characteristic poses challenges for gradient descent methods,
which are better suited for problems with smoother landscapes of the cost
function. In photonics, two different starting points in the optimization
process most often lead to two different outcomes[14]. This is in contrast to
mechanics, where different starting points do not typically result in
significant variations in the outcomes[37, 3]. Moreover, the structures
produced by these algorithms often exhibit excessive intricacies, which can
pose challenges for fabrication and hinder their commercialization
potential[1, 38].

Global optimization for photonics. It is now evident that early attempts at
using genetic algorithms for optimizing simple photonic structures were
unsuccessful[39]. Specifically designed heuristics have in general failed to
produce structures that are convincing enough to persuade the community to
embrace global numerical optimization as a tool for years. Moreover, for
continuous problems, even modern global optimization algorithms often yield
unsatisfactory solutions when the parameter space has a high dimensionality,
as in TO. The first successes of gradient-based TO may thus suggest that
global optimization, as in mechanics, cannot compete with TO on inverse design
problems.

However, in our experience with continuous problems and modern optimization
algorithms, we have found that keeping the number of parameters relatively
low, typically below 100, usually suffices to give BBO algorithms enough
degrees of freedom to yield innovative and convincing results[4]. This
approach requires limiting the number of degrees of freedom, a strategy also
referred to as Parametric Optimization (PO). It becomes however challenging to
make a fair comparison between BBO and TO, as they operate in very different
regimes.

In addition, global optimization algorithms may be valuable in problems that
are typically considered in TO, due to the discretization requirement. Many
topology optimization problems are inherently discrete, and the use of
continuous parameters is primarily for leveraging gradient-based methods. The
variety of the results obtained using continuous parameters[14] suggests that
making the problem continuous is also making it more difficult, for example
because local optimization algorithms get stuck in minima with intermediate
values of the refractive index. In such cases, discrete global optimization
algorithms may actually offer advantages and recent results seem to indicate
that such methods are able to yield efficient solutions in a relatively
consistent way, even spontaneously generating welcome features like symmetry
[13] and the dissapearance of intermediate refractive index values. Once
again, the challenge is to identify the most suitable algorithm for a given
problem, which can only be done with comprehensive algorithm libraries. In the
present paper, we focus on PO, due to its advantages: suggested designs are
frequently smooth and possible to manufacture, and the dimensionality
typically remains moderate (dozens to hundreds).

In our opinion, global optimization algorithms are too often overlooked due to
the challenges posed by high-dimensional parameter spaces. We underline that
most often, finding the right algorithm makes a tremendous difference and
that, given the inherent complexity of optimization, a rigorous methodology
must be applied to achieve satisfactory solutions. While these methods may be
computationally demanding, the advent of parallelism has made many global
algorithms (which are usually intrinsically parallel) significantly more
efficient, making them increasingly relevant.

###  II.2 Typical photonics test cases

We have chosen three test cases that are typical of problems that can be
encountered in photonics, on which we applied the methods we recommend. A
Jupyter Notebook[9] hosted on the Colab platform is provided in which we show
how cost functions can be easily defined in the context of multilayered
structures using the PyMoosh module[40] and how all the algorithms implemented
in the comprehensive Nevergrad library[8] can be subsequently tested. These
three test cases are: a dielectric mirror optimized for high reflectivity, an
ellipsometric problem, and a photovoltaic cell optimized for absorption inside
the active material.

High reflectivity problem. The cost function is defined as $1-R$ where $R$
represents the reflectance at the working wavelength of 600 nm of a
multilayered structure with alternating refractive indexes ($n_{1}=1.4$ and
$n_{2}=1.8$, starting with $n_{2}$). Therefore, the optimization algorithms
will aim to find highly reflective structures. As represented in Fig. 2(b) we
considered three sub-cases (i) a structure with only 10 layers (minibragg),
which is representative of low dimensionality problems in photonics (ii) a
structure with 20 layers (bragg), a number we feel marks the beginning of the
domain of inverse design problems and (iii) a structure with 40 layers
(bigbragg), where some algorithms may face challenges. The starting points of
each optimization are chosen randomly.

Figure 2: Illustration of the Bragg mirror benchmark problem. The objective is
to maximize the reflectance of a multilayered structure composed of two
alternating materials at a wavelength of 600 nm, as illustrated in a). The
figures presented to illustrate the different problem size in b) are typical
of an initial geometry for the optimization.

The thickness of each layer is taken between 0 and the thickness of a half-
wave layer for the lower index. When considering only a single working
wavelength, adding or removing a half-wave layer to a given layer has no
impact on the optical response (this is called an absent layer). Therefore,
considering larger thicknesses would only introduce degenerate minima. We
underline that letting the refractive index vary would lead to the same
solutions, as the optimization would converge to the highest index contrast
between two layers[4]. This can be considered as physically expected, as this
is known to be the most efficient choice to modulate the optical response of a
multilayer[41]. The Bragg mirror, a periodic solution, has been identified as
the best solution to this problem so far[4]. It is a local optimum, and it
outperforms any disordered solution, suggesting that it might be a regular,
principled solution that might be the optimum – even though it is not, as of
now, proven. This is why we have selected it as a test case, for which we know
the (likely) global optimum, and called it “Bragg”.

Ellipsometry problem. The objective here is to find the material and the
thickness of a reference layer, knowing its reflectance properties, obtained
using spectroscopic ellipsometry. This step is required to extract the desired
information from ellipsometric measurements. For a tested layer with a given
material and thickness, we compute the relation $e=\frac{r_{p}}{r_{s}}$ where
$r_{p}$ is the reflectance of the whole structure in TE polarization and
$r_{s}$ is the reflectance in TM polarization, for a wavelength range of 400 -
800 nm and for a given incidence angle of 40°. The quantity we minimize here
is the difference between the $e$ computed from the tested layer and the $e$
computed from the reference layer.

The thickness of the tested layer is taken between 30 and 250 nm, and we
assume the refractive index is comprised between 1.1 and 3. This simple
problem is illustrated in Fig.3. Problems in ellipsometry are generally more
complicated but highly dependent on the response model assumed for the
material. Our simple dispersion-less and absorption-less test case can however
be adapted easily thanks to the architecture of PyMoosh[7].

We underline that in many ways, this case mirrors the practical challenges
faced by photonics researchers. It illustrates the common situation where
researchers design structures, characterized by a small number of parameters,
and then engage in optimization to determine the upper limits of their
performance[42, 43, 44].

Figure 3: Illustration of the Ellipsometry benchmark problem. The objective is
to retrieve the thickness and the material of an unknown layer, knowing its
reflectance in TE and TM polarization for wavelengths comprised between 400 nm
and 800 nm.

Photovoltaics problem. The objective here is to find the best possible
antireflective coating in normal incidence with a multilayer structure made of
two alternating materials (permittivities of 2 and 3 respectively) on an
amorphous hydrogenated silicon substrate (30000 nm). The quantity to be
maximized here is the short circuit current in the 375 - 750 nm range,
assuming a quantum yield equal to 1, as described in [4, 5].

The cost function $f$ is one minus the efficiency defined as the ratio between
the short circuit current $j_{sc}$ and the maximum achievable current
$j_{max}$ if all photons are converted into electron-hole pairs :

$f=1-\frac{j_{sc}}{j_{max}}$ (1)

with

$j_{sc}=\int
A(\lambda)\frac{dI}{d\lambda}.\frac{\mbox{e}\lambda}{\mbox{h}\mbox{c}}\,\mbox{d}\lambda,$
(2)

where $A(\lambda)$ is the absorbance of the active layer, e, h and c are
respectively the elementary charge, the Planck constant and the speed of light
in vacuum and the spectral density of the illumination $\frac{dI}{d\lambda}$
is given by the solar spectrum[45].

The problem is illustrated in Fig. 4. As for the high reflectivity problem we
considered three sub-cases (i) a structure with only 10 layers
(photovoltaics), (ii) a structure with 20 layers (bigphotovoltaics), and (iii)
a structure with 32 layers (hugephotovoltaics).

Figure 4: Illustration of the Photovoltaic benchmark problem. The objective is
to maximise the absorption in the silicon layer, using an antireflective
coating composed of an alternating two materials multilayered structure.

##  III Basic tools for optimization: algorithms and observables

As discussed before, many algorithms and methods are available to perform
optimization and eventually inverse design. It is therefore necessary to
define some observables and basic tools to compare the performances of
algorithms in a reliable way. This section presents first some well known
categories of algorithms. It is shown how these algorithms can all be run
through the Nevergrad platform[8]. Then we defined relevant observables, i.e.
quantities allowing discussions about the performances of algorithms and their
results. All the discussions in this section are illustrated by results of the
Jupyter Notebook provided in the supplementary material[9].

###  III.1 Algorithms categories

Many BBO platforms exist, including CMA (Covariance Adaptation Algorithm)[46],
AX (Adaptive eXperiments)[47], HyperOpt [48], SMAC (Sequential Model-based
Algorithm Configuration [49]), NLOPT[50] and Nevergrad[8]: Nevergrad is a free
Python library which imports all those ones and others. We present in this
section the algorithms from Nevergrad used in our Jupyter Notebook
experiments[9], and organize them based on their respective categories.

Mathematical programming. Methods available for gradient-based optimization
have been adapted to the black-box setting. For example the limited-memory
Broyden-Fletcher-Goldfarb-Shanno LBFGS [51] method can be applied, with
derivatives computed by finite differences. Note that BFGS and LBFGS based on
finite differences are parallel: the $d+1$ evaluations needed for evaluating
the gradient by finite differences can be done in parallel. We also include
GradBFGS, which is using the gradient without having to “pay” the cost of
$d+1$ evaluations per iteration: GradBFGS is therefore not directly
comparable, in the sense that it uses more information than other methods and
is not as parallel as others.

Genetic and Evolutionary Algorithms. One of the most well known families of
algorithms are evolutionary and/or genetic algorithms. In an evolutionary
algorithm, each step of the process creates and evaluates the cost function of
a “population” of structures, preserving the “individuals” (i.e., the
structures) that are better than those belonging to the population of the
previous step (also called the previous “generation”). The different genetic
algorithms have different ways of creating the new generation of structures.
Most include steps such as mutations and crossovers between structure
information.

We underline that the historical algorithms, in which for instance the binary
coding of the parameters is considered as its genetic code, are considered
obsolete in the optimization community because of their well documented lack
of efficiency[52, 53] and their use should be avoided. A lot of more efficient
algorithms, inspired by these early ideas, have emerged over the years, with
an improved efficiency. One of the most well known and efficient overall seems
to be Differential Evolution (DE [54]), an algorithm that has many variants.
One of the most classical variants is presented in Table 1 and the different
formulas that can be used for the mutation to define new variants are
presented in Table 2.

In our experiments, we often use the quasi-opposite variant QODE [55], which
randomly draws half the initial population, and then, for each point $x$, adds
$-r\times x$ to the population with $r$ randomly uniformly drawn in $[0,1]$
(assuming that the center is $0$: otherwise $c-r\times(x-c)$ with $c$ the
center of the optimization domain). We also include QNDE (Quasi-Newton DE),
which runs QODE for half the budget and then finishes with BFGS with finite
differences: this is a memetic algorithm, which runs first a global algorithm
and then a local one. Both seem particularly efficient on our test cases.

Randomly draw the population $x_{1}$,…,$x_{30}$ in the search space.

while There is time do

{$//$ This is one iteration}

for each $x$ in the population do

Randomly draw $a$, $b$, $c$ and $d$ individuals (i.e. for us designs) in the
population.

Compute $y$ using one of the equations in Table 2.{$//$Mutation operator}

{$//$Crossover operator: for loop below.}

Randomly draw $CR$ in $\\{1,\dots,d\\}$. {$//$ $d$ is the dimension}

for each $i$ in $\\{1,\dots,d\\}$ do

if $i=CR$ or with probability $\frac{1}{2}$ then

$z_{i}\leftarrow y_{i}$.

else

$z_{i}\leftarrow x_{i}$.

end if

if $z$ has better fitness than $x$ then

$x$ is replaced by $z$ in the population.{End of crossover operator}

end if

end for

end for

end while

Table 1: Pseudo-code of DE. $a,b,c,d$ are distinct, randomly drawn individuals
in the population and $best$ refers to the best design so far. The $y$
individual is obtained by a mutation formula (see Tab. 2) and $CR$ is the
crossover rate, indicating the percentage of the mutant $y$ used in the
creation of a new individual $z$. DE/rand/1: $\displaystyle y(x)=$
$\displaystyle a+F_{1}(b-c)$ DE/best/1: $\displaystyle y(x)=$ $\displaystyle
best+F_{1}(a-b)$ DE/randToBest/1: $\displaystyle y(x)=$ $\displaystyle
c+F_{1}(a-b)+F_{2}(best-c)$ DE/currToBest/1: $\displaystyle y(x)=$
$\displaystyle x+F_{1}(a-b)+F_{2}(best-x)$ DE/rand/2: $\displaystyle y(x)=$
$\displaystyle a+F_{1}(a-b+c-d)$ DE/best/2: $\displaystyle y(x)=$
$\displaystyle best+F1(a-b+c-d)$ Table 2:  Various DE formulas. $a,b,c,d$ are
distinct, randomly drawn individuals in the population. We see that the
mutated variant of $x$, namely $y(x)$, is in some cases, by design,
independent of $x$. $best$ refers to the best design so far. Typically,
$F_{1}=F_{2}=0.8$ and $CR=\frac{1}{2}$.

All the variants of DE typically perform well for a large parameter space
dimension (but smaller than topological optimization), including quite
irregular cost functions such as the ones which appear in photonics. Most
winning algorithms for Large-Scale Global Optimization (LSGO) competitions use
a variant of DE[27].

There are many reasons why DE can be a sensible default choice - among them
its simplicity, its low computational overhead, and its good performance
overall. We also point out that the rise of parallel computing resources makes
DE and other global methods like PSO (see below) way faster: whereas
parallelizing BFGS or other gradient-based methods is tricky, DE (as
implemented and parametrized in the present paper and in most implementations)
is just 30 times faster with 30 cores than on a sequential machine, and can be
way more parallel than that with a simple adaptation of parameters (typically
increasing the population size). Besides this natural parallelism advantage,
many balck-box optimization methods can include arbitrary parameters
(including discrete parameters, e.g. choosing a material in a list), addingin
non-differentiable stochastic manufacturing errors, adding multiple objective
functions and handling worst-case analysis over input angle or over a range of
wavelength.

Moreover, DE is built on a differential formula to create new structures based
on existing ones. This means that if multiple structures in the current
population share the same characteristics (e.g., in the present framework, the
same block of layers with identical thicknesses and materials), those
characteristics will probably be preserved in the creation of the new
candidate solution. In photonics, this leads to the emergence of modular
structures with distinct blocks of layers, each having well-defined optical
functions. We believe that this specific property of the optimization
algorithm contributes to its remarkable efficiency addressing photonic
problems [4, 56].

We have thus natively included DE in PyMoosh[7]. We invite the interested
reader to consult our Notebook to run simple optimizations with it. The Bragg
and Photovoltaics cases show how to use PyMoosh to perform an optimization
(one run), while the Ellipsometry case illustrates how to perform several
runs. We underline also that parameters can be grouped and evolve together
during the optimization, which is especially relevant if they are physically
linked (like the thickness and refractive index of the same layer). This can
be considered as a variant of DE and is called structured optimization[56]. It
is not shown here because in the cases we study, grouping parameters does not
make sense.

Other evolutionary methods include Particle Swarm Optimization (PSO [57]): a
population of particles with random initial positions and velocities is
attracted towards good search points found so far. PSO is known as a robust
method, dedicated to rugged cost functions. SQOPSO (Special Quasi-Opposite
PSO) is a variant in Nevergrad, based on quasi-opposite sampling. Many other
search methods are based on updating a Gaussian probability distribution for
matching the best points so far[58].

Typically, Evolution Strategies (ES) iteratively update a probability
distribution model, that is used to optimized the likelihood for generation of
good candidates. A well-known example is Covariance Matrix Adaptation (CMA
[59]) which updates the entire covariance matrix of a multivariate normal
distribution (MND). The CMA algorithm samples candidates using a MND,
evaluates them, and then modifies the MND using the evaluation results.

While CMA typically performs best on the Black Box Optimization Benchmark
(BBOB[26]), simpler methods such as the $(1+1)$ evolution strategy with one-
fifth rule[60] are still good in many cases due to its ability to quickly
adapt the search scale.

Bayesian Optimization. In a Bayesian optimization process[61, 62], the
algorithm uses a probability distribution (a Gaussian process) for modeling
the cost function and the uncertainty, and updates it with each new structure
tested. The model provides, for each point in the search space, an estimated
probability distribution for the cost function value of that search point.

Therefore, it is possible, at each iteration, to define the Expected
Improvement (EI) for a search point $x$:
$EI(x)=\mathbb{E}_{w}\max(0,m-CostModel(x,w))$, where $m$ is the best cost so
far and $CostModel(x,w)$ is the current statistical model of the cost
function. $CostModel(x,w)$ depends on a random $w$, because this model is
statistical. The value of $w\mapsto CostModel(x,w)$ is updated each time a new
cost value is available: it is frequently a Gaussian process.

This gives an estimation of what structure to try next: we search for $x$ in
the search space such that $EI(x)$ is approximately maximal – this requires a
side optimization for each iteration, hopefully orders of magnitude faster
than the original problem.

Many variants of this approach have been proposed in the literature, with many
parameters for the stochastic model. By design, Bayesian optimization (BO)
works best with problems with a small parameter space dimension and for
relatively smooth functions. Also, the probabilistic interpolation process,
which is necessary to know which structure to try next, can be expensive,
possibly becoming computationally more expensive than the cost function
itself. The authors in [63] tested Bayesian Optimization in a limited budget
scenario but other, computationally cheaper algorithms were still performing
well. In our code, we include two Bayesian Optimization methods for the low-
dimensional and low budget problems, namely BO[64] and PCABO[65].

Other black-box optimization methods. Other methods include pattern search
methods, such as e.g. the Nelder-Mead method[66], which iteratively updates a
simplex of search points. Also, optimization wizards (also known as hyper-
Heuristics) are defined by automatically selecting and combining existing
methods[67]: NGOpt and NGOptRW are examples of wizards included in Nevergrad.
Methods based on reinforcement learning (RL) have also been defined and
applied on machine learning problems[68], though simple classical evolutionary
methods were actually competitive [69]. Please note that[70] mentions
difficulties for reproducing results in some RL papers.

###  III.2 Observables

For the vast majority of continuous optimization problems, proving that a
solution provided by an algorithm is the best possible solution in the
optimization domain is essentially impossible. In addition, many optimization
algorithms are non-deterministic or sensitive to the initialization, which
means each run of the algorithm will lead to a different outcome.

As a consequence, it is possible to perform many optimization runs and still
not be able to firmly determine whether the best of these solutions is good
enough to stop looking for other, potentially better, outcomes. Yet, observing
how exactly each different run progresses towards a solution and considering
the solutions that are produced statistically, yields crucial information that
may help the user gain confidence in the best solutions produced or, on the
contrary, indicate that these solutions cannot be trusted.

Convergence curves. The first observable, which is widely used, is the
convergence curve that can be produced for each run of a given algorithm by
plotting the value of the cost function of the solution that the algorithm
recommends (typically the best solution found so far) as a function of the
number of iterations. When multiple runs are launched, convergence curves can
be either drawn for each run or an averaged convergence can be obtained. Both
essentially provide the same kind of information: whether most of the runs
have settled on a given solution (if they have almost all reached a plateau
for many iterations), or if further iterations have a chance to improve the
outcome. Fig. 5 (resp. Fig. 6) presents such an example of individual curves
for each run (resp. aggregated curves with average convergence).

Budget. Not all iterations of all algorithms have the same computational cost
or evaluate the cost function the same number of times. An iteration for a
genetic or evolutionary algorithm typically corresponds to a new generation
and thus the evaluation of a few dozen of solutions. For some other
algorithms, an iteration requires the evaluation of a single new solution. A
way to compare two algorithms fairly is to discuss their performances for a
given number of evaluations of the cost function. The maximal number of
evaluations allowed is called the budget and it has to be chosen before
running the algorithms. Each run of an algorithm is then stopped when the
budget has been reached. Of course, this does not take into account the
computational cost of the optimization algorithm itself, which should be
discussed when not negligible compared to the cost function values.

Consistency curve. Convergence curves for different runs allow to determine
whether the chosen budget is large enough to reach a satisfactory convergence
for each run. However, since the different runs can produce completely
different results, the variability of the different results has to be
discussed. This can be done by plotting what we call the consistency curve
(Fig. 7). This curve is obtained by plotting the different values of the cost
function reached at the end of each run, sorted from the lowest to the
highest. This is generally done for the highest budget allowed. When such a
curve displays a plateau, then the same solution, or at least solutions with
the same cost, has been found several times by the algorithm. A large plateau
for the lowest value of the cost function thus means the best solution has
been found a large number of times. This reinforces the confidence that can be
placed in that solution.

Box plots. It is not always relevant to draw the consistency curve, especially
when comparing a large number of algorithms. This curve can be summarized by a
box plot as e.g. in Fig. 8. While a box plot will not allow to observe any
plateau and thus bears less information than a consistency curve, it gives an
immediate and easily readable access to the statistical properties of the
values of the cost function.

Illustration of local optima with local optimization. In addition, local
optimization algorithms can be used to estimate whether there is a limited or
a large number of local minima and whether the best solution found by any
other algorithm has a large attraction basin or not. A simple way to do so is
to launch multiple runs of the algorithm with starting points chosen randomly
in the optimization domain, then make sure all the runs have converged (e.g.
BFGS has a stop criterion essentially guaranteeing a local minimum has been
reached) and finally plot the resulting consistency curve. If we were to run a
very large number of such optimizations, the consistency curve should present
plateaus corresponding different local minima – eventually several minima can
present identical values of the cost function, because of symmetries of the
problem and not be distinguishable. The width of a plateau would allow to
estimate the volume of the corresponding attraction basin. Running BFGS with
randomly chosen starting point may appear here, beyond a sound optimization
approach, as a way of estimating the difficulty of a problem. A lot of
different results suggest a lot of different minima and if the best result is
rarely found, this also means it has a relatively small attraction basin
compared to the size of the optimization domain. These characteristics are
indicative of a difficult problem.

##  IV Running and discussing experiments

With the idea that anyone should be able to build easily on our work, we
provide a Jupyter Notebook hosted on the Colab platform[9] which can be
downloaded and run, or run online allowing to reproduce all the results below.
The objective of this Notebook is to show how to optimize the typical
photonics test cases, either using a DE algorithm natively implemented in
PyMoosh or, for comparing a wide range of methods and/or creating new ones,
relying on Nevergrad. Launching a sufficient number of runs to be able to
compare algorithms afterwards is called running an experiment in the
optimization community. The first part of the online code can be modified to
choose the experimental setup. We have chosen a limited set of representative
algorithms to be included in our experiment for readability, but modifying the
set is straightforward. For instance, Bayesian Optimization results are shown
only in the lowest dimensional setting, the only one in which it was
performing well.

###  IV.1 Algorithm comparisons with Nevergrad

Using the results produced by our code, we first illustrate here how the
observables defined above can be used to discuss the results of the
experiments.

It is common to represent the convergence curves for all the different runs of
a single algorithm, as shown in Fig. 5. While this does indeed provide useful
information about the robustness of an algorithm’s convergence, as discussed
above, it is convenient to represent the average of the convergence curves
(rather than all curves) when benchmarking a lot of different algorithms – and
adding the standard deviation can also be done as shown in Fig. 6 to reduce
the information loss.

For Photovoltaics we see that convergence seems reached only for variants of
DE. For the other algorithms, convergence is not completely clear and
increasing the budget might thus help. BFGS and LBFGS use finite differences,
and therefore the $d+1$ first points are essentially horizontal because they
are the finite differences computations. GradBFGS uses the real gradient, so
it does not suffer from this, and any flat convergence curve can be attributed
to the fact that, at a random starting point and at early iterations,
following the gradient improves the efficiency only slightly.

Figure 5:  Convergence curves. Convergence curves for different runs for
different algorithms showing the value of the cost function (the lower, the
better) as a function of the budget consumed by the optimization. The
information brought by such individual curves is useful but difficult to read
when comparing even a small number of algorithms, that is why the average
convergence is often plotted, as shown in Fig. 6.

For small budgets, local algorithms, that tend to converge faster, can be a
good choice. This holds especially in cases with lower dimensionality like the
Ellipsometry example, because local optima are often less dangerous in a low
dimensional setting.

However, Fig. 7 shows that for more difficult problems and larger budgets, DE
variants are efficient despite the large number of local minima with high cost
function values shown on the consistency curves of BFGS and the like.

Figure 6: Averaged convergence curves. Averaged convergence curves for
different algorithms (26 runs each), with standard deviation visible in
ellipses, as a function of the number of iterations. The algorithms are
ordered by average performance for the maximum budget (the lower the better).
Using an average convergence curves allows to compare more algorithms[16].
Figure 7:  Consistency curves. Each algorithm has been run several times (from
15 to 25 times) with the highest budget (from $10^{3}$ to 32000 depending on
the case, as shown on the convergence curves Fig. 5). The fitness values for
each run are sorted from left to right in ascending order. The lower the curve
and the flatter, the more efficient and reliable the algorithm can be
considered (good reproducibility the convergence value). A plateau means, that
a certain fitness value and probably the same solution has been found multiple
times. The results for BFGS and its variants correspond to different local
minima if the algorithm has converged (this is often the case when the budget
is large).

Showing on the same graph the consistency curves for a large variety of
algorithms is not convenient, according figures can be difficult to interpret
(see Fig. 7). Each consistency curve can however be summarized using a box
plot, allowing for a fair benchmarking between various algorithms, as shown in
Fig. 8 (for clarity, the results are not shown for all the algorithms that we
have tested). As can be seen, for a case like Ellipsometry with few
parameters, local methods perform excellently. When the landscape is simple
enough, mathematical programming tools (ranging from BFGS and the like to
Nelder-Meade) seem the most adapted. Differential Evolution variants are
frequently good, though slower than other methods on these problems.

For higher dimensions, and cases considered as inverse problems, DE variants
perform well overall, in full agreement with previous results[4]. We underline
that DE variants are tested with population size 30, which means that they
would run, without change, 30 times faster on a machine with 30 cores. In
addition, they can easily be adapted to parallel machines above this 30
threshold just by increasing the population size parameter, whereas local
algorithms are way harder to parallelize. Combining DE and local methods
sometimes brings improvement (e.g. QNDE vs QODE in the low dimensional case)
and is thus recommended[67]. QNDE is almost as good as QODE in most cases, and
competitive in the low-dimensional case.

Judging from the consistency curves and the boxplots, the CMA algorithm seems
as efficient as DE on the Bragg problems, finding the Bragg mirror most of the
time and outperforming local algorithms even if it presents a slower
convergence. However, increasing the dimension of the problem for
Photovoltaics makes the problem significantly harder for this algorithm, to
the point where it is not better than local optimization algorithms for the
largest problem. We see in that case, how comparing a global optimization
algorithm to local algorithms with a random start, allows to assess whether an
algorithm has an interest for a given problem.

Bayesian Optimization (BO) is computationally expensive, and this
computational overhead does not provide significant improvements, especially
in the more complex test cases. Furthermore, different BO implementations lead
to different results because there are many free parameters in Bayesian
optimization. We underline that the cases studied here are limited to
relatively fast cost functions and representations by a moderate number of
parameters. This confirms that BO is well suited for problems with relatively
low degree of freedom and for which the cost function is expensive to compute,
which imposes a hard limit on the budget when the landscape is still
complex[16, 71, 62].

Figure 8:  Performance box-plots for comparing algorithms. Box-plots
representing the distribution of the minimum values of the cost function
reached at the end of each run for a given algorithm for different test cases.
The results are shown for the highest budget, ranging from 1000 for
Ellipsommetry to 32000 for HugePhotovoltaics, as shown Fig. 5. Each box
presents the first quartile, the third quartile, a line at the median, and
dashed lines are added at 1.5 times the inter-quartile range.

###  IV.2 Physical analysis

Analyzing physically both the solutions produced and how they have been
produced (by which algorithm, in which conditions and how fast) sheds a new
light both on the solution and on the optimization process itself. When a
problem is grounded in physics, its distinct characteristics can and have to
be leveraged to gain deeper insights.

Low dimension. In the Ellipsometry case, local algorithms perform perfectly
and find the solution reliably and quickly. This is typical for a simple
landscape with few local minima. Physically, this means the geometry is simple
enough and the considered layers sufficiently thin, to have a very low number
of resonances.

Multiple resonances and local minima. The Bragg test case is more interesting.
Even for ten layers, local algorithms (versions of BFGS) most often fail to
find the Bragg mirror. Since the algorithms have converged and since local
algorithms tend to converge to local minima, this means local minima are
numerous, even for such a relatively simple problem with an intermediate
number of parameters. The fact that BFGS most often misses the Bragg mirror
means that the optimization domain is large compared to the attraction basin
of the Bragg mirror. This can be expected because each resonance in the
structure is likely to produce a local minimum or maximum in any considered
physical quantity. With a global thickness for the structure which is beyond 2
µm, at least as many resonances as for a homogeneous layer of the lower index
material with a thickness ranging from 0 to 2 µm can be expected. This means
at the very least a dozen resonances and thus a dozen local minima exist. The
algorithms CMA and DE perform noticeably better than local algorithms, which
means they are adapted to the problem.

Regularity and modularity. As explained above, for DE, if the same values for
some parameters can be found in all the individuals they will be conserved
throughout the optimization. If a subset of parameters then describes a
subpart of the structure, DE will keep this part unchanged. This makes DE
uniquely equipped to produce modular structures, in which specific parts play
a well defined role.

In the Bragg case for instance, we suspect that once a few layers begin to
present the characteristics of a Bragg mirror at any place in the structure,
the algorithm keeps them intact and explores the sub-space of the optimization
domain that this defines. As this sub-space has a smaller volume, it is then
easier for the algorithm to find the rest of the Bragg mirror. This allows to
understand why DE seems to generate easily structures that are periodic also.
Somehow, periodic structures can be seen as modular also.

This is illustrated Fig. 9, where the best results (whatever the algorithm)
are shown for different budgets. For budgets ranging from 100 to 1000, the
best solutions are produced by GradBFGS, as shown Fig 5 c). For the lower
budget, the structure appears disordered (no algorithm has converged at that
point) and the spectrum indicates that a relatively narrow anti-resonance
producing a large reflectance has been found. When the budget is around 1000,
GradBFGS has converged and proposes typically the solutions shown Fig. 9 b),
that show some regularity and decent performances. The spectrum (presenting a
broad reflectance peak) and the structure seem to indicate that an anti-
resonance is not enough to explain such performances. However for larger
budget, DE has converged towards the Bragg mirror shown in c), which is
perfectly periodic and has the highest efficiency as well as the broadest
reflectance spectrally.

Figure 9: Best structures (right) and their associated reflectance spectrum
(left) obtained for the Bragg case with 20 layers for a budget of a) 100 b)
1000 c) 20000\. The light, resp. dark color represents the high (1.8), resp.
low (1.4) refractive index material. The grey color represents substrate and
superstrate.

The Photovoltaics case is a perfect example of a modular structure. Fig. 10
shows the best results for all the Photovoltaics cases, all produced by QODE.
It is important to notice that, whatever the number of layers and even though
the structures have been obtained for independent runs, the three upper and up
to 5 lower layers are common to all the structures. For the largest number of
layers, periodical patterns (alternating 120 nm resp. 150 nm thick layers for
permittivy 3.0 resp. 2.0) are appearing. A previous study has shown that the
upper and lower layers allow light to be coupled in and out of the central
photonic crystal more easily[5]. The fact that they appear consistently
indicates their physical importance. They allow to diminish the oscillations
characteristics of Bragg mirrors outside the bandgap that can be seen Fig 9 c)
on the reflectance and are detrimental to the absorption.

Figure 10: Best structures (right) and their absorptance spectrum (left)
obtained for the Photolvoltaics case with a) 10 layers b) 20 layers and c) 32
layers. The light, resp. dark color represents the high (3), resp. low (2)
permittivity material. The grey color represents substrate and superstrat.

Sometimes DE is able to generate structures like chirped dielectric
mirrors[4], which locally resembles a Bragg mirror but whose parameters
changes gradually and smoothly. Such a structure cannot be considered
periodic. It is however modular, as the different parts of the structure are
obviously able to reflect different parts of the spectrum. We call this kind
of structures regular. A periodic structure is obviously regular, but a
regular structure is not periodic in that sense.

From the cases presented here and our experience, it can generally be expected
that periodic or regular structures will also be the most efficient to have an
impact on light because it is a wave – so that it is particularly sensitive to
periodicity. The Bragg mirror, for instance, has a lower value of the cost
function used here than any disordered structure that we have generated. The
thickness values for the AR coating point towards a structure containing a
photonic crystal, even if they are not as precise as with the Bragg case
because of the irregularities of the solar spectrum. We underline that we have
run these optimizations a very large number of times on these two cases and
never found any better solution that the ones presented here. This strongly
suggests that regularity should be expected and even sought after[11] and we
underline that in many results from the literature, periodic pattern often
seem to emerge spontaneously (see Fig. 11 below) .

Robustness. In photonics, robustness to fabrication defaults is always
desirable. A robust structure presents an efficiency which will not change
much when the parameters are slightly modified. From an optimization point of
view, this simply means that the attraction basin of the desired solution is
large. Evaluating the robustness of a solution is thus computationally costly
because it involves modifying a large number of parameters and computing the
change in the cost function. It can be tempting to include the robustness in
the cost function, however in our experience this is not especially necessary
to produce regular or periodic structures.

In the Bragg case, a robust solution is a solution for which the peak of
reflectance at the working wavelength will not be displaced too much when the
structure is slightly modified. As a consequence, the larger such a peak, the
more robust the solution can be considered. The largest reflectance on the
broadest spectral range is usually produced by regular structures like
photonic crystals. Which means that the most robust solutions will also be, at
least in that case, the most periodic. And in general, periodic solutions are
indeed associated with robustness. We underline thus, that finding regular
solutions should be the priority, and that this will probably also ensure the
robustness of the design.

##  V Good practices for optimizing photonic structures

The optimization of photonic structures thus presents quite a few specific
characteristics, which influence the strategy to be followed in this
particular context. While the computational cost may put a limit to the
problems that can be tackled, we give below a list of strategies that, applied
together, constitute a methodology for the optimization of photonic
structures.

One of the most important questions is, when we should stop the optimization.
It is impossible to prove that a solution is optimal. However, there are ways
do determine the quality of a solution. We give below criteria that can help
to determine whether a solution is satisfying – meaning the optimization can
be stopped – or not.

###  V.1 Optimization methodology

Our methodology consists in maximizing the information extracted from the
observables we have defined above, and to make use of specific characteristics
of photonic problems to gradually establish confidence in the generated
solutions. We leave other strategies that could also make interesting
solutions emerge for further work (e.g. multi-objective optimization, or the
use of manufacturing noise in the cost function).

Convergence curves for the determination of the budget. The presence of
plateaus in a convergence curve indicates that the algorithm has converged,
suggesting that the budget is adequate. On the contrary, the absence of
plateaus on the consistency curve suggests that the budget should be
increased.

Systematic use of consistency curves for checking local minima. Since global
algorithms are usually non deterministic, it is useful to launch multiple
runs, despite the computational cost, in order to conduct a statistical
analysis with a consistency curve. Ideally, the consistency curve exhibits a
plateau at its minimum value (see Fig. 7a). However, this is not always the
case (see Fig. 7f), and then this indicates that the problem might not be
fully solved.

Adding degrees of freedom. In photonics problems, it is generally
straightforward to gradually increase the number of elements that can be
optimized without changing the nature of the problem. In the cases presented
above, this can be done by increasing the number of layers of the structure,
but it could for example also be through a decrease in the discretization
stepsize, e.g. in topology optimization. Structures with different numbers of
layers can then be compared in terms of performances. It can be generally
expected that increasing the number of degrees of freedom leads to improved
optimized performances. Plotting the minimum value of the cost function as a
function of complexity, represented by the number of layers or elements in the
structure, can provide valuable insights: e.g. if increasing the number of
layers does not improve performance, which indicates that the difficulty of
the problem has also increased, continuing that path is likely pointless[5].

Parametrization bias awareness: meaningful representations make the
optimization problem easier. In parametric optimization, when a handful of
parameters are needed to describe a relatively complex device, the choice of
these parameters and of their limits (which defines the optimization domain)
is crucial. More precisely, these initial choices may introduce biases, favor
certain types of algorithms or make the convergence more difficult. When, for
instance, the parametrization is chosen such, that subsets of parameters
correspond to components of the structure, algorithms like DE are particularly
efficient. In DE, when a component of a structure is widely spread in the
whole population, it might be exactly preserved through the iterations,
whereas many other algorithms keep perturbing all variables. DE has the
distinctive property to produce, spontaneously, modular structures which are
exactly preserved through iterations. When the different parameters do not
describe a part of the structure but a more global property, other kinds of
algorithms might be more relevant, as has been underlined in previous
works[16].

Sensitivity to the optimization domain: choosing bounds. In many cases, the
imposed constraints strongly control the emerging solutions. For example,
using a medium with an extremely high refractive index (typically infinite) is
a simple but not realistic way to reflect light completely. The constraints on
the refractive index values are therefore the fundamental reason for the
production of Bragg mirrors as a solution. However, there are instances where
the constraints become too demanding, making it difficult to find satisfactory
solutions. It is important in that case to verify whether some parameters are
stuck at the optimization domain boundary (i.e. if the boundary constraints
are active). On the other hand, when a satisfactory solution is produced, it
can be informative to add or remove constraints or expand the optimization
domain. Bragg mirrors tend to emerge, whether the refractive indices are
allowed to vary within certain limits or are imposed, with the latter case
being straightforward. A clear understanding of the conditions under which a
solution is generated also contributes to building confidence.

Leverage your physical intuition. We underline that in optimization, there are
no rules a priori. If, for instance, it makes sense physically to modify a
solution by hand, this should not be considered forbidden or ”cheating”,
especially when it seems that the algorithm is stuck in a local optimum that
can be criticized based on a physical reasoning[72]. The limits of the
optimization range can also be set to encourage the algorithm to explore areas
where promising solutions are likely based on physics, or to stay within
specific functional ranges[73]. This approach usually makes the task easier
for the algorithm and provides solutions that are easier to understand and
more satisfying.

The physical intuition is what often determines the conditions in which the
optimization takes place and what allows to detect parametrization biases, or
even that a problem is not well posed enough for any algorithm to find a
satisfying solution. It should never be overlooked.

###  V.2 Assessing the quality of a solution

Usually, no solution can be proven optimal, due to the impossibility to
explore the entire space of possible solutions or to locate all the local
minima. Therefore, it becomes necessary to establish criteria that enhance the
confidence in a solution. Besides the optimization criteria above based on
optimization observables, a physical analysis is possible, as developped in
the present section. When enough confidence in a solution has been built, it
can be deemed satisfying.

Consistency. A solution that has been obtained at least more than once
inspires greater confidence. If it is obtained repeatedly, it might correspond
to a plateau of similarly good solutions on the consistency curve of the most
efficient algorithm. In that case the solution can be deemed truly consistent
and particularly trustworthy. This is most often not the case. When the best
solution is obtained only for a single run, this should be considered
indicative of a local minimum. We regret that, except in a few cases[14],
elements allowing to assess the consistency of a solution, even if this is not
a consistency curve, are generally not given.

Spontaneous emergence of regularity. In photonics, periodical or regular
structures are ubiquitous. This can be directly attributed to the wave nature
of the underlying phenomena. As underlined in a pioneering optimization
study[11], “the emergence of periodicity suggests that periodicity is a
principal condition for strong light manipulation.” Many studies have shown
the emergence of partially periodic structures, even when the solution lacks
complete periodicity or regularity (like for chirped dielectric mirrors
typically[4]). However, when an algorithm proposes completely periodic
structures as a solution, they naturally inspire more satisfaction. Based on
our experience, we have yet to encounter a simple problem where fully
disordered structures outperform regular ones in terms of efficiency.

A symmetric structure can also be considered more regular. We believe that the
spontaneous emergence of symmetry also reinforce the confidence that can be
placed in a structure. We underline that this is rare, as symmetry is almost
always imposed (this tend to simplify problems noticeably).

Overall, we disagree with the notion that prioritizing performance over
aesthetics is necessary, as we believe that both aspects are inherently
intertwined in photonics. We underline again that in the Bragg mirror
benchmark, no disordered structure has ever presented a better performance
than a Bragg mirror with the same number of layers. In the Photovoltaics case,
the irregularities can be linked to the noise in the solar spectrum and, as a
consequence, in the cost function. Irregularities in that case improve the
performances, but the periodic pattern is still distinguishable and is central
for the overall efficiency. In the case of multilayered structures, regularity
or periodicity may be more likely to emerge, due to the relative simplicity of
the geometry. However, we underline that in the literature, relatively regular
patterns (in the sense defined above) emerge all the time, as shown Fig. 11.
Sometimes the patterns look unfinished, perhaps indicating the solution can be
further improved – which is likely if a local algorithm has been used. In our
experience, very regular or periodic pattern, including spontaneously
symmetrical ones, can be generated also in more complex setups[13].

Physical interpretability. The solutions that are most satisfying are those
that can be readily understood from a physical point of view. We underline
that, generally, only periodic, regular or modular structures can be truly
understood. This is more difficult for completely disordered structures,
except if the disorder itself is tailored, which cannot be ruled out. The
absence of physical interpretability is likely what hinders the widespread
adoption of optimization as a daily research tool within the community.
Sometimes, the solutions can be studied and fully understood a posteriori[5].
Although this does not guarantee optimality, this is at least a good reason to
stop looking for alternative solutions: any solution that is comprehensible
and understandable can serve as a valuable source of inspiration for manually
generated structures and can offer valuable design rules. In rare situations,
algorithms can produce solutions that resemble patterns found in nature, on
insect wings for instance, which have evolved for optical purposes. Although
these occurrences are uncommon, they can be highly satisfying, as they align
with the concept of evolution as a form of optimization. However, due to their
infrequency, they cannot be included in the above criteria. In the case of the
reflection problem which we have considered in this work, this criterion is
obviously fulfilled too as Bragg mirror are commonly found in nature.

Figure 11: Spontaneous emergence of regularity or periodicity. When
periodicity or regularity emerge, solutions appear much more satisfying and
are more likely to be analyzed physically. (a) When the position of
nanopillars of silicon are optimized to enhanced the magnetic Purcell factor,
using Differential Evolution algorithm, 2 rings cavities spontaneously
emerge[74]. (b) Two-dimensional dielectric metalenses designs obtained by
topology optimization[75] with a relatively low number of parameters. (c)
Silicon carbide optical cavity obtained using gradient-based optimization [76]
showing a very regular pattern (holes of increasing dimensions in the
waveguide). (d) When gratings of rectangular blocks of chitin are optimized to
maximizing the scattered reflection of a 450 nm wavelength in the first
diffraction orders, while minimizing the specular reflection of a range of
visible wavelengths quasi perfect checkboards with interdigitated blocks are
produced[4].

We underline that these criteria to determine whether a solution is satisfying
or not can be applied to any inverse design problem, whether it is solved by
topology optimization, shape optimization, parametric optimization or any
other technique. We advise all authors, as far as possible, to publish all the
necessary information so that other researchers can reproduce and verify the
quality of optimization results. When this is done, particularly interesting
and thorough discussions become possible[14].

##  VI Conclusion

In this paper, we present different types of popular optimization algorithms
and compare them using the Nevergrad platform on three typical photonic
problems. For problems with low dimensionality, like ellipsometry or more
generally parameters retrieving, a lot of algorithms seem adapted, including
local optimization algorithms or Bayesian optimization. For problems closer to
inverse design, not all approaches are effective, because of the large number
of local minima, frequent in photonics problems. We have shown how algorithms
can be rigorously and thoroughly compared. The illustrated examples show that
variants of the Differential Evolution algorithm are highly efficient on a
large range of photonic problems. Additionally, we developed a methodology and
offer some advice for conducting high-performance design optimizations and
evaluating the quality of a solution. Finally, in supplement to the
manuscript, we provide python jupyter notebooks that illustrate our workflow
and can be used to reproduce the presented benchmarks.

We meant this work as a tutorial, but also as a warning. Optimization is
difficult because of its unique curse: it is never possible to guarantee that
a solution is optimum, or even close to it, making science much more
difficult. This has not been too much of a problem in mechanics because the
problem landscape is generally simple enough for understandable physical
reasons. However, photonics problems are comparatively more complicated as our
results here show: even the simplest problems present numerous local minima.
This must lead to extra caution. The field of optimization is awash in dubious
claims of novelty and efficiency. In order to avoid a similar reproducibility
crisis in photonics, adopting an open science approach is imperative: data
regarding the different runs of optimization should be published, codes should
be shared and both should be discussed[14, 77]. We are convinced that the
optimization of photonics structures is a work intensive domain, and that
neither a single method nor a single team will be enough to uncover what
optimization can bring in terms of innovation.

There also is a danger to deem a solution satisfactory when it should not: it
is to miss innovative and more efficient structures. We underline that the
question of how reflectance can be maximized using numerical optimization,
central in the present work, has been considered 25 years ago[39], exactly
when DE was published, but because of the lack of efficient algorithms at the
time, the resulting structures were disordered and not convincing enough to
foster further work in that direction. Perhaps because the problem had already
been studied, more than 20 years have passed until it has been examined
again[4], producing the convincing results this work is based on. Given the
potential of inverse design but the difficulty to find structures that would
be commercially viable[1], there is a danger in giving up too soon and miss
particularly efficient structures. Fortunately, physical analysis of
structures seems to be a powerful tool to discuss both the solutions and the
optimization process itself.

We have shown in the present work how modern numerical tools have made the use
of optimization much simpler and efficient in photonics. Even for well defined
functioning regimes suggested by physics itself, with a relatively low number
of parameters, numerically optimizing a photonic structure often yields
unexpectedly high performances. We underline that numerical optimization is
now able to produce photonics structures that can be understood. This
constitutes a complete change compared to the times when inefficient
algorithms (as the first genetic algorithms) were producing disordered and
impossible to understand results. Open science approaches now allow any
researcher in the field to use such tools easily. We hope that our work will
encourage fellow researchers within the community to seamlessly integrate
optimization tools into their routine practice and to join the effort in
discovering novel and more efficient structures to address the challenges of
the future.

## Acknowledgements

A.M. is an Academy CAP 20-25 chair holder. He acknowledges the support
received from the Agence Nationale de la Recherche of the French government
through the program Investissements d’Avenir (16-IDEX-0001 CAP 20-25). This
work was supported by the International Research Center ”Innovation
Transportation and Production Systems” of the Clermont-Ferrand I-SITE CAP
20-25. P.R.W. acknowledges the support of the French Agence Nationale de la
Recherche (ANR) under grant ANR-22-CE24-0002 (project NAINOS), and from the
Toulouse high performance computing facility CALMIP (grant p20010).

## References

  * Molesky  _et al._ [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, Nature Photonics 12, 659 (2018).
  * Bendsoe and Sigmund [2003] M. P. Bendsoe and O. Sigmund,  _Topology optimization: theory, methods, and applications_ (Springer Science & Business Media, 2003).
  * Sigmund [2011] O. Sigmund, Structural and Multidisciplinary Optimization 43, 589 (2011).
  * Barry  _et al._ [2020] M. A. Barry, V. Berthier, B. D. Wilts, M.-C. Cambourieux, P. Bennet, R. Pollès, O. Teytaud, E. Centeno, N. Biais, and A. Moreau, Scientific reports 10, 1 (2020).
  * Bennet  _et al._ [2021] P. Bennet, P. Juillet, S. Ibrahim, V. Berthier, M. A. Barry, F. Réveret, A. Bousquet, O. Teytaud, E. Centeno, and A. Moreau, Physical Review B 103, 125135 (2021).
  * Sörensen [2015] K. Sörensen, International Transactions in Operational Research 22, 3 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1111/itor.12001 .
  * pym [2023] Pymoosh, https://pypi.org/project/PyMoosh/ (2023), verified: 2022-06-10.
  * Rapin and Teytaud [2018] J. Rapin and O. Teytaud, Nevergrad - A gradient-free optimization platform, https://GitHub.com/FacebookResearch/Nevergrad (2018).
  * Wiecha [2023] P. Wiecha, Tuto Multilayer Optimization, https://doi.org/10.5281/zenodo.8347023 (2023).
  * Kim  _et al._ [1996] J. Kim, H.-B. Lee, H. K. Jung, S.-Y. Hahn, C. Cheon, and H. Kim, IEEE Transactions on Magnetics 32, 1250 (1996).
  * Gondarenko  _et al._ [2006] A. Gondarenko, S. Preble, J. Robinson, L. Chen, H. Lipson, and M. Lipson, Physical review letters 96, 143904 (2006).
  * Piggott  _et al._ [2015] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, Nature Photonics 9, 374 (2015).
  * Teytaud  _et al._ [2022] O. Teytaud, P. Bennet, and A. Moreau, Photonics and Nanostructures - Fundamentals and Applications 52, 101072 (2022).
  * Su  _et al._ [2020] L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A. Petykiewicz, and J. Vučković, Applied Physics Reviews 7 (2020).
  * Tikhonravov  _et al._ [1996] A. V. Tikhonravov, M. K. Trubetskov, and G. W. DeBell, Applied optics 35, 5493 (1996).
  * Schneider  _et al._ [2019] P.-I. Schneider, X. Garcia Santiago, V. Soltwisch, M. Hammerschmidt, S. Burger, and C. Rockstuhl, ACS Photonics 6, 2726 (2019).
  * Mahfoud [1995] S. W. Mahfoud,  _Niching Methods for Genetic Algorithms_ , Ph.D. thesis, University of Illinois at Urbana-Champaign (1995).
  * Shir [2012] O. M. Shir, Niching in evolutionary algorithms, in  _Handbook of Natural Computing_ (Springer Berlin Heidelberg, 2012) pp. 1035–1069.
  * Petrowski [1996] A. Petrowski, in  _Proceedings of IEEE International Conference on Evolutionary Computation_ (1996) pp. 798–803.
  * Jakšić  _et al._ [2023] Z. Jakšić, S. Devi, O. Jakšić, and K. Guha, Biomimetics 8, 10.3390/biomimetics8030278 (2023).
  * Wolpert and Macready [1997] D. Wolpert and W. Macready, IEEE Transactions on Evolutionary Computation 1, 67 (1997).
  * Markov [2023] I. L. Markov, The false dawn: Reevaluating google’s reinforcement learning for chip macro placement (2023), arXiv:2306.09633 [cs.LG] .
  * Gould  _et al._ [2003] N. I. Gould, D. Orban, and P. L. Toint, ACM Transactions on Mathematical Software (TOMS) 29, 373 (2003).
  * Gould  _et al._ [2015] N. I. Gould, D. Orban, and P. L. Toint, Computational optimization and applications 60, 545 (2015).
  * DIMACS [2021] DIMACS, Dimacs implementation challenge, http://dimacs.rutgers.edu/programs/challenge (2021).
  * Hansen  _et al._ [2009] N. Hansen, A. Auger, S. Finck, and R. Ros,  _Real-Parameter Black-Box Optimization Benchmarking 2009: Experimental setup_ , Tech. Rep. RR-6828 (INRIA, France, 2009).
  * Li  _et al._ [2013] X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, in  _CEC 2013 proceedings_ (2013).
  * Gallagher and Saleem [2018] M. Gallagher and S. Saleem, in  _PPSN’18 workshop_ (2018).
  * Häse  _et al._ [2021] F. Häse, M. Aldeghi, R. J. Hickman, L. M. Roch, M. Christensen, E. Liles, J. E. Hein, and A. Aspuru-Guzik, Machine Learning: Science and Technology 2, 035021 (2021).
  * Lee  _et al._ [2018] J. Lee, M. X. Grey, S. Ha, T. Kunz, S. Jain, Y. Ye, S. S. Srinivasa, M. Stilman, and C. K. Liu, Journal of Open Source Software 3, 500 (2018).
  * Coumans and Bai [2017] E. Coumans and Y. Bai, Pybullet, a python module for physics simulation in robotics, games and machine learning (2017).
  * Todorov  _et al._ [2012] E. Todorov, T. Erez, and Y. Tassa, in  _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems._ (2012) pp. 5026–5033.
  * Brockman  _et al._ [2016] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, arXiv preprint arXiv:1606.01540 (2016).
  * Musgrave  _et al._ [2020] K. Musgrave, S. J. Belongie, and S. Lim, CoRR abs/2003.08505 (2020), arXiv:2003.08505 .
  * Kapoor and Narayanan [2022] S. Kapoor and A. Narayanan, Leakage and the reproducibility crisis in ml-based science (2022).
  * Cea [1986] J. Cea, ESAIM: Mathematical Modelling and Numerical Analysis 20, 371 (1986).
  * Wang  _et al._ [2003] M. Y. Wang, X. Wang, and D. Guo, Computer methods in applied mechanics and engineering 192, 227 (2003).
  * Wang  _et al._ [2018] F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, arXiv preprint arXiv:1810.02417 (2018).
  * Martin  _et al._ [1995] S. Martin, J. Rivory, and M. Schoenauer, applied Optics 34, 2247 (1995).
  * Moreau [2023] A. Moreau, PyMoosh (2023).
  * Tikhonravov [1993] A. V. Tikhonravov, Applied Optics 32, 5417 (1993).
  * Smaali  _et al._ [2021] R. Smaali, T. Taliercio, A. Moreau, and E. Centeno, Applied Physics Letters 119 (2021).
  * Centeno  _et al._ [2021] E. Centeno, E. Alvear-Cabezón, R. Smaali, A. Moreau, and T. Taliercio, Semiconductor Science and Technology 36, 085014 (2021).
  * Moreau  _et al._ [2007] A. Moreau, C. Lafarge, N. Laurent, K. Edee, and G. Granet, Journal of Optics A: Pure and Applied Optics 9, 165 (2007).
  * Santbergen  _et al._ [2010] R. Santbergen, J. Goud, M. Zeman, J. van Roosmalen, and R. C. van Zolingen, Solar energy materials and solar cells 94, 715 (2010).
  * Hansen  _et al._ [2019] N. Hansen, Y. Akimoto, and P. Baudis, CMA-ES/pycma on Github, Zenodo, DOI:10.5281/zenodo.2559634 (2019).
  * FacebookResearch [2020] FacebookResearch, Ax - adaptive experimentation, ax.dev (2020).
  * Bergstra  _et al._ [2015] J. Bergstra, B. Komer, C. Eliasmith, D. Yamins, and D. D. Cox, Computational Science and Discovery 8, 014008 (2015).
  * Hutter  _et al._ [2011] F. Hutter, H. H. Hoos, and K. Leyton-Brown, in  _LION_, Lecture Notes in Computer Science, Vol. 6683, edited by C. A. C. Coello (Springer, 2011) pp. 507–523.
  * Johnson [1994] S. G. Johnson, The nlopt nonlinear-optimization package (1994).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, Math. Program. 45, 503 (1989).
  * Janikow  _et al._ [1991] C. Z. Janikow, Z. Michalewicz, _et al._ , in  _ICGA_ , Vol. 1991 (1991) pp. 31–36.
  * Herrera  _et al._ [1998] F. Herrera, M. Lozano, and J. L. Verdegay, Artificial intelligence review 12, 265 (1998).
  * Storn and Price [1997] R. Storn and K. Price, J. of Global Optimization 11, 341 (1997).
  * Rahnamayan  _et al._ [2007] S. Rahnamayan, H. R. Tizhoosh, and M. M. A. Salama, in  _2007 IEEE Congress on Evolutionary Computation_ (2007) pp. 2229–2236.
  * Rapin  _et al._ [2020] J. Rapin, P. Bennet, E. Centeno, D. Haziza, A. Moreau, and O. Teytaud, in  _Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion_ (2020) pp. 1599–1607.
  * Kennedy and Eberhart [1995] J. Kennedy and R. C. Eberhart, in  _Proceedings of the IEEE International Conference on Neural Networks_ (1995) pp. 1942–1948.
  * Beyer [2001] H.-G. Beyer,  _The Theory of Evolution Strategies_ , Natural Computing Series (Springer, Heideberg, 2001).
  * Hansen and Ostermeier [2003] N. Hansen and A. Ostermeier, Evolutionary Computation 11 (2003).
  * Rechenberg [1973] I. Rechenberg,  _Evolutionsstrategie Optimierung technischer Systeme nach Prinzipien der biologischen Evolution_ (Friedrich Frommann Verlag, Stuttgart-Bad Cannstatt, 1973).
  * Jones  _et al._ [1998] D. R. Jones, M. Schonlau, and W. J. Welch, Journal of Global Optimization 13, 455 (1998).
  * Elsawy  _et al._ [2021] M. M. Elsawy, A. Gourdin, M. Binois, R. Duvigneau, D. Felbacq, S. Khadir, P. Genevet, and S. Lanteri, ACS photonics 8, 2498 (2021).
  * Hutter  _et al._ [2013] F. Hutter, H. Hoos, and K. Leyton-Brown, in  _Proceedings of the 15th Annual Conference Companion on Genetic and Evolutionary Computation_, GECCO ’13 Companion (Association for Computing Machinery, New York, NY, USA, 2013) p. 1209–1216.
  * Nogueira [14 ] F. Nogueira, Bayesian Optimization: Open source constrained global optimization tool for Python (2014–).
  * Raponi  _et al._ [2020] E. Raponi, H. Wang, M. Bujny, S. Boria, and C. Doerr, CoRR abs/2007.00925 (2020), 2007.00925 .
  * Nelder and Mead [1965] J. A. Nelder and R. Mead, Computer Journal 7, 308 (1965).
  * Liu  _et al._ [2020] J. Liu, A. Moreau, M. Preuss, J. Rapin, B. Roziere, F. Teytaud, and O. Teytaud, in  _Proceedings of the 2020 Genetic and Evolutionary Computation Conference_ , GECCO ’20 (2020) p. 620–628.
  * Zoph and Le [2017] B. Zoph and Q. V. Le, in  _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_ (OpenReview.net, 2017).
  * Real  _et al._ [2017] E. Real, S. Moore, A. Selle, S. Saxena, Y. L. Suematsu, Q. V. Le, and A. Kurakin, CoRR abs/1703.01041 (2017), arXiv:1703.01041 .
  * Cheng  _et al._ [2023] C.-K. Cheng, A. B. Kahng, S. Kundu, Y. Wang, and Z. Wang, in  _Proceedings of the 2023 International Symposium on Physical Design_, ISPD ’23 (Association for Computing Machinery, New York, NY, USA, 2023) p. 158–166.
  * Elsawy  _et al._ [2020] M. M. Elsawy, S. Lanteri, R. Duvigneau, J. A. Fan, and P. Genevet, Laser & Photonics Reviews 14, 1900445 (2020).
  * Frellsen  _et al._ [2016] L. F. Frellsen, Y. Ding, O. Sigmund, and L. H. Frandsen, Optics Express 24, 16866 (2016), publisher: Optical Society of America.
  * Moreau  _et al._ [2012] A. Moreau, R. Smaali, E. Centeno, and C. Seassal, Journal of Applied Physics 111 (2012).
  * Brûlé  _et al._ [2022] Y. Brûlé, P. Wiecha, A. Cuche, V. Paillard, and G. C. Des Francs, Optics Express 30, 20360 (2022).
  * Christiansen and Sigmund [2021] R. E. Christiansen and O. Sigmund, JOSA B 38, 510 (2021).
  * Yang  _et al._ [2023] J. Yang, M. A. Guidry, D. M. Lukin, K. Yang, and J. Vučković, arXiv preprint arXiv:2303.17079 (2023).
  * Jiang  _et al._ [2020] J. Jiang, R. Lupoiu, E. W. Wang, D. Sell, J. P. Hugonin, P. Lalanne, and J. A. Fan, Optics express 28, 13670 (2020).