text:
\begin{document}

\title{A Covariance Matrix Self-Adaptation Evolution
  Strategy for Optimization under Linear Constraints}

\author{Patrick~Spettel,~Hans-Georg~Beyer,~and~Michael~Hellwig
  \thanks{Manuscript received Month xx, xxxx; revised Month xx, xxxx;
    accepted Month xx, xxxx. Date of publication Month xx, xxxx; date
    of current version Month xx, xxxx. The authors thank Asma Atamna
    for providing the simulation code for the ES with
    augmented Lagrangian constraint handling. This work was supported by the
    Austrian Science Fund FWF under grant P29651-N32.}
  \thanks{The authors are with the
    Research Center Process and Product Engineering
    at the Vorarlberg University of Applied Sciences,
    Dornbirn, Austria (e-mail: Patrick.Spettel@fhv.at; Hans-Georg.Beyer@fhv.at;
    Michael.Hellwig@fhv.at).}
  \thanks{Digital Object Identifier xx.xxxx/TEVC.xxxx.xxxxxxx}}

\markboth{IEEE Transactions on Evolutionary Computation,~Vol.~xx,
  No.~x, Month~xxxx}
{Spettel \MakeLowercase{\textit{et al.}}:
  A CMSA-ES for Linear Constrained Optimization}

\IEEEpubid{\begin{minipage}{\textwidth}\centering
    xxxx--xxxX~\copyright~xxxx IEEE. Personal use is permitted, but
    republication/redistribution requires IEEE permission.\\
    See \url{http://www.ieee.org/publications_standards/publications/
      rights/index.html} for more information.\end{minipage}}

{\maketitle}
\IEEEpeerreviewmaketitle

\begin{abstract}
  This paper addresses the development of a covariance
  matrix self-adaptation evolution strategy (CMSA-ES) for solving
  optimization problems with linear constraints.
  The proposed algorithm is referred to as Linear Constraint CMSA-ES
  (lcCMSA-ES).
  It uses a specially built
  mutation operator together with repair by projection
  to satisfy the constraints. The lcCMSA-ES evolves
  itself on a linear manifold defined by the constraints.
  The objective function is only evaluated
  at feasible search points (interior point method). This is
  a property often required in application domains such as simulation
  optimization
  and finite element methods. The algorithm is tested on a
  variety of different test problems revealing considerable results.
\end{abstract}

\begin{IEEEkeywords}
  Constrained Optimization,
  Covariance Matrix Self-Adaptation Evolution Strategy,
  Black-Box Optimization Benchmarking,
  Interior Point Optimization Method
\end{IEEEkeywords}

\section{Introduction}
\IEEEPARstart{T}{he} Covariance Matrix Self-Adaptation Evolution Strategy
(CMSA-ES)~\cite{Beyer2008} variant called Constraint CMSA-ES (cCMSA-ES)
was proposed in~\cite{Beyer2012}. It showed promising
results in portfolio optimization applications.
Therefore, further research on ES design principles
for constrained optimization problems is of interest.
The CMSA-ES and linear constraints are chosen
as a first step. This is because the
CMSA-ES is arguably one of the most simple variants
of Covariance Matrix Adaptation (CMA)
ESs~\cite{Hansen2001,Hansen2003}.
In addition, linear constraints are the most simple
constraints after box constraints. But this does not mean that it is only
of theoretical interest. Such problems occur in practical
applications. Examples include risk management
in finance~\cite{McCarl1997AppliedMathProgramming,
Pflug2007ModelingMeasuringManagingRisk,Yan2007Multi},
agriculture~\cite{Kaylen1987RiskModeling},
hybrid dynamic systems~\cite{MINLP27},
model predictive control~\cite{Buijs2002},
controlled perturbation for tabular data~\cite{Castro2006}, and
optimization of heat exchanger
networks~\cite{Gorji2011OptimizationHeatExchanger}. Further,
the CEC 2011 real world optimization problem competition contains
an electrical transmission pricing
problem based on the IEEE 30 bus
system~\cite[Prob. 9 in Sec. 8]{CEC2011RealWorldProblems}.
Another example is a problem from the area of chemistry.
The chemical composition of a
complex mixture under chemical equilibrium conditions has to be
determined. This problem is described in detail
in~\cite[pp. 47-49]{Bracken1968ApplicationsNonlinearProgramming}.

Evolutionary Algorithms (EAs) in general are well-suited for scenarios
in which objective function and/or constraint functions
cannot be expressed in terms of (exact) mathematical expressions.
Moreover, if that information is incomplete or if that information
is hidden in a black-box, EAs are a good choice as well.
Such methods are commonly referred
to as direct search, derivative-free, or zeroth-order
methods~\cite{Wright1996, Trosset1997, Kolda2003, Conn2009}.
In fact, the unconstrained case has been studied well.
In addition, there is a wealth of proposals in the field of Evolutionary
Computation dealing with
constraints in real-parameter optimization, see e.g.~\cite{MezuraCoello2011}.
This field is mainly dominated by Particle Swarm Optimization (PSO)
algorithms and Differential
Evolution~(DE)~\cite{Huang2006ConSADE,MalliPeddi2010DE,Takahama2010EpsilonDEGA}.
For the case of constrained discrete optimization, it has been shown
that turning constrained optimization problems into multi-objective optimization
problems can achieve better performance than the single-objective
variant with a penalty approach
for some constrained combinatorial optimization problems,
e.g.,~\cite{Neumann2006,Friedrich2010,Qian2015}.

ESs for constrained optimization have not yet been studied
extensively. Early work includes the $(1+1)$-ES
for the axis-aligned corridor model~\cite{Rechenberg1973},
the $(1,\lambda)$-ES for the same environment~\cite{Schwefel1981},
and the $(1+1)$-ES for a constrained, discus-like
function~\cite{Beyer1989Evolutionsverfahren}.
Moreover, a stochastic ranking approach was proposed
in~\cite{RunarssonStochasticRanking2000}. An ES
for constrained optimization was proposed
in~\cite{MezuraMontesConstrainedES2005}.
For the CMA-ESs, in addition to the cCMSA-ES~\cite{Beyer2012},
a $(1+1)$-CMA-ES based on active covariance matrix adaptation
is presented in~\cite{Arnold2012ActiveCovMatrixUpdate}.
There exists an extension of this idea to a $(\mu,\lambda)$-CMA-ES
motivated by an application in the area of rocket
design~\cite{Chocat2015ModifiedCMAESRocketDesign}.
In~\cite{Arnold2016ActiveSetES} an Active-Set ES that is able to handle
constraints is described.
Further, an ES with augmented Lagrangian constraint
handling is presented in~\cite{Atamna2017LagrangianES}.

\IEEEpubidadjcol
The cCMSA-ES~\cite{Beyer2012} uses
two mechanisms to ensure the feasibility
with respect to box and equality constraints. First, mutations are generated in
such a way that they lie on the hypersurface defined by the equality constraint.
Second, a repair mechanism is applied to offspring violating the
box-constraints.

Being based on these ideas, the contribution of this work
is a theoretically motivated and principled algorithm design.
The proposed algorithm is an interior point ES. It is able to optimize
a black-box objective function subject to general linear constraints.
The peculiarity of this design is that the algorithm treats
the function $f$ to be optimized as a black-box.
However,
only \emph{feasible} candidate solutions are used to query
the black-box $f$. This is in contrast to most of the evolutionary
algorithms proposed for constrained black-box optimization.
However, it is a property often required in the field of
simulation optimization, e.g. in Computational Fluid Dynamics (CFD)
optimizations.
In CFD optimizations, constraint violations on simulator
input parameters
may result in simulator crashes.\footnote{Note that although interior point
methods only run the simulator with feasible input parameters,
problems can still occur.
Feasible input parameters can lead to crashes because of
possible parameter combinations that were never thought of.
And constraints are often defined on simulator outputs.
Such cases need a different treatment not subject of this paper.}
In~\cite{LeDigabelWildTaxonomy2015}, concrete real-world examples are
provided for different constraint types.
Among those examples, a ground water optimization
problem~\cite{KannanWild2012Groundwater}
is provided for which (some) constraints are not allowed to be violated.
Because the simulator only supports extraction but not injection, the
lower bounds on the pumping rate values must hold for the simulation.
Physical requirements like that usually prohibit the violation of
(some) constraints. Further, it is an important
property for problems that cannot tolerate even small infeasibility
rates. Such problems can be a topic in finance or business applications,
i.e., the optimization of a function subject to a constant amount of
total money in the system.
Moreover, the mutation operator and
the repair method are specially designed. The ES moves completely
on a linear manifold defined by the constraints. For this design, the
theory is an essential part.

The rest of the paper is organized as follows.
In \Cref{sec:optproblem} the optimization problem
is presented. Then, the proposed algorithm
is described in \Cref{sec:algo} and simulation results are
presented in \Cref{sec:expeval}. Finally, \Cref{sec:conclusion}
summarizes the main results and provides an outlook.

\vspace{0.2cm}\noindent\emph{Notations} Boldface $\mathbf{x} \in \mathbb{R}^D$
is a
column vector with $D$ real-valued components. $\mathbf{x}^T$ is its
transpose. $x_d$ and equivalently ($\mathbf{x})_d$ denote the $d$-th element
of a vector $\mathbf{x}$. $x_{(k:D)}$ and equivalently $(\mathbf{x})_{(k:D)}$
are the order statistic notations, i.e., they denote the $k$-th smallest
of the $D$ elements of the vector $\mathbf{x}$.
$||\mathbf{x}|| = \sqrt{\sum_{d = 1}^D {x_d}^2}$
denotes the euclidean norm ($\ell_2$ norm) and
$||\mathbf{x}||_1 = \sum_{d = 1}^D |x_d|$ the $\ell_1$ norm.
$\mathbf{X}$ is a matrix, $\mathbf{X}^T$ its transpose, and
$\mathbf{X}^+$ its pseudoinverse. $\mathbf{0}$ is the vector
or matrix (depending on the context) with all elements equal to zero.
$\mathbf{I}$ is the identity matrix. $\mathcal{N}(\bm{\mu}, \mathbf{C})$
denotes the multivariate normal distribution with mean $\bm{\mu}$
and covariance matrix
$\mathbf{C}$. $\mathcal{N}(\mu, \sigma^2)$
is written for the normal distribution with mean $\mu$ and
variance $\sigma^2$.
$\mathcal{U}[a,b]$ represents the continuous uniform distribution
with lower bound $a$
and upper bound $b$. The symbol $\sim$ means ``distributed according to'',
$\gg$ ``much greater than'', and $\simeq$ ``asymptotically equal''.
A superscript $\mathbf{x}^{(g)}$ stands for the element
in the $g$-th generation. $\langle \mathbf{x} \rangle$ denotes the
mean (also centroid) of a parameter
of the $\mu$ best individuals of a population,
e.g. $\langle \tilde{\mathbf{z}} \rangle
= \frac{1}{\mu}\sum_{m = 1}^\mu \tilde{\mathbf{z}}_{(m:\lambda)}$.

\section{Optimization Problem}
\label{sec:optproblem}
We consider the (non-linear) optimization problem
\begin{subequations}
  \label{sec:optproblem:eqn:problem}
  \begin{align}
    &f(\mathbf{x}) \rightarrow \text{min!}
    \label{sec:optproblem:eqn:objective}\\
    \text{s.t. }&\mathbf{A}\mathbf{x} = \mathbf{b}
    \label{sec:optproblem:eqn:linearsystem}\\
    &\mathbf{x} \ge \mathbf{0}
    \label{sec:optproblem:eqn:nonnegative}
  \end{align}
\end{subequations}
where $f: \mathbb{R}^D \rightarrow \mathbb{R}$,
$\mathbf{A} \in \mathbb{R}^{K \times D}$, $\mathbf{x} \in \mathbb{R}^D$,
$\mathbf{b} \in \mathbb{R}^K$.
Note that \Cref{sec:optproblem:eqn:linearsystem,sec:optproblem:eqn:nonnegative}
form a linear constraint system in standard form. Any linear
inequality constraints and bounds can be transformed into an equivalent
problem of this form.
In particular, a vector satisfying the constraints in the
transformed system, also
satisfies the constraints in the original system.
A method for this transformation is presented in the supplementary
material (\Cref{sec:appendix:standardoptproblem2}).

\section{Algorithm}
\label{sec:algo}
Based on the Covariance Matrix Self-Adaptation
Evolution Strategy (CMSA-ES)~\cite{Beyer2008},
we propose an algorithm for dealing with problem
\eqref{sec:optproblem:eqn:problem}.
In particular, we describe the design of the linear constraint
$(\mu/\mu_I, \lambda)$-CMSA-ES (lcCMSA-ES).

The $(\mu/\mu_I, \lambda)$-CMSA-ES makes use of
$\sigma$-self-adaptation and a simplified covariance update.
It starts from a parental individual $\mathbf{x}$
(\cref{sec:algo:alg:lccmsaes:findinhsol} in
\Cref{sec:algo:alg:lccmsaes} corresponds to this step)
and an initial mutation strength $\sigma$
(part of \cref{sec:algo:alg:lccmsaes:initparams} in
\Cref{sec:algo:alg:lccmsaes}).
The generational loop consists of two main steps.
First, $\lambda$ offspring are generated. For every offspring $l$, its
mutation strength $\tilde{\sigma}_l$ is sampled from a log-normal distribution.
Then, the offspring's parameter vector is sampled from a multi-variate
normal distribution $\mathcal{N}(\mathbf{x}, \tilde{\sigma}_l \mathbf{C})$
(\crefrange{sec:algo:alg:lccmsaes:lambdaloopstart}
{sec:algo:alg:lccmsaes:lambdaloopend} in \Cref{sec:algo:alg:lccmsaes} is the
offspring generation (extended for the constrained case)).
Second, the parental individual is updated for the next
generation. For this, the parameter vectors and the mutation strengths
of the best $\mu$ offspring are averaged. Then, the covariance is
updated (\crefrange{sec:algo:alg:lccmsaes:rankoffspring}
{sec:algo:alg:lccmsaes:updateC} in \Cref{sec:algo:alg:lccmsaes}).
This completes the steps for one iteration of the inner loop.

The method proposed in this paper represents an interior point method,
i.e., the individuals
evaluated in \Cref{sec:optproblem:eqn:objective}
are always feasible.
In other words, the ES moves inside the feasible region
while searching for the optimum.
Concretely, this is realized by starting off with an initial
centroid that is feasible (\Cref{sec:algo:subsec:initcentroid}).
Mutation (\Cref{sec:algo:subsec:mutation})
is performed in the null space of $\mathbf{A}$ in order
to keep the mutated individuals feasible with respect to
\Cref{sec:optproblem:eqn:linearsystem}.
It is possible that after mutation
\Cref{sec:optproblem:eqn:nonnegative} is violated.
Repair (\Cref{sec:algo:subsec:repair})
by projection to the positive orthant is performed in such cases.

\subsection{Initialization of the Initial Centroid}
\label{sec:algo:subsec:initcentroid}
The problem for the initialization of the initial parental centroid
is to find an $\mathbf{x}$ such that $\mathbf{A}\mathbf{x} = \mathbf{b}$
and $\mathbf{x} \ge \mathbf{0}$. For the first part,
the system of linear equations can be solved. This solution
$\mathbf{x}$ possibly violates the second part. In that case the repair
approach of the ES (\Cref{sec:algo:subsec:repair})
is applied to this initial solution.
Under the assumption that the linear system solver and the repair operator
are deterministic, this whole initialization is deterministic. This means
that for the same $\mathbf{A}$ the same $\mathbf{x}$ is computed every time.
In order to use the algorithm in a restarted fashion, random
initialization is important. For this, an initial random movement
of $\mathbf{x}$ can be obtained in a similar way as for the
mutation (\Cref{sec:algo:subsec:mutation}).

\subsection{Mutation}
\label{sec:algo:subsec:mutation}
The goal of mutation is to introduce variation to the population
of candidate solutions. In the unconstrained
case one could simply add a random vector to the parental centroid.
But the constraints make this more complicated. Since the proposed
method is an interior point method, mutated individuals that violate
the constraints must be repaired. One option would be to design a mutation
operator that does not violate any constraints. Here, the approach is a
mixture of both. The mutation operator does not violate the linear equality
constraints but it does possibly violate the non-negativity constraint.
The latter case is handled through repair by projection. For the former
note that
$\mathbf{A}(\mathbf{x}_{\text{inh}} + \mathbf{x}_{\text{h}}) = \mathbf{b}$
where $\mathbf{x}_{\text{inh}}$ is an inhomogeneous
and $\mathbf{x}_{\text{h}}$ is a homogeneous solution.
Thus, $\mathbf{x}_{\text{h}} \in \text{null}(\mathbf{A})$ and
therefore $\mathbf{A}\mathbf{x}_{\text{h}} = \mathbf{0}$.
Let $N$ be the dimension and $\mathbf{B} \in \mathbb{R}^{D \times N}$
an orthonormal basis of the null space
$\text{null}(\mathbf{A})$,\,i.e.,
$\mathbf{B}^{T}\mathbf{B} = \mathbf{I}$ and
$\mathbf{A}\mathbf{B} = \mathbf{0}$ hold.
Mutations are performed in $\text{null}(\mathbf{A})$ and
therefore do not violate
\Cref{sec:optproblem:eqn:linearsystem}.
This means that a mutation vector $\mathbf{s}$ in the null space is sampled
from a normal distribution with zero mean and the covariance matrix
$\mathbf{C}$,
i.e., $\mathbf{s} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}^{N \times N})$.
Transforming this $\mathbf{s}$ into the parameter space
and scaling it with the mutation strength $\sigma$
yields a mutation vector $\mathbf{z} = \sigma \mathbf{B} \mathbf{s}$
in the parameter space.
This $\mathbf{z}$ can be added to the parental centroid (or the initial
centroid for initial value randomization)
\begin{equation}
  \mathbf{x}^{(g + 1)} = \mathbf{x}^{(g)} + \mathbf{z} =
                        \mathbf{x}^{(g)} + \sigma \mathbf{B} \mathbf{s}.
\end{equation}
Assuming the parental centroid satisfies the linear constraints
$\mathbf{A}\mathbf{x}^{(g)} = \mathbf{b}$, a short calculation shows
that the mutated offspring fulfills them as well:
\begin{align}
  \begin{split}
    \mathbf{A}\mathbf{x}^{(g + 1)}
       &= \mathbf{A}(\mathbf{x}^{(g)} + \sigma \mathbf{B} \mathbf{s})
       = \mathbf{A}\mathbf{x}^{(g)} + \mathbf{A}(\sigma\mathbf{B}\mathbf{s})\\
       &= \mathbf{A}\mathbf{x}^{(g)} +
          \sigma
          \underbrace{(\mathbf{A}\mathbf{B})}_{\mathbf{0}^{K \times N}}
                                                          \mathbf{s}
       = \mathbf{b} + \mathbf{0}^{K \times 1}
       = \mathbf{b}.
  \end{split}
\end{align}

\subsection{Repair}
\label{sec:algo:subsec:repair}
\Cref{sec:optproblem:eqn:linearsystem}
is not violated through mutation. But violation of
\Cref{sec:optproblem:eqn:nonnegative}
must be dealt with. The approach followed here is repair by projection
onto the positive orthant.

Although repair by minimal change is intuitively the most plausible
approach, it is worth noting that
the repair in the ES does not have to be optimal.
It is enough to find a point on the positive orthant
that is approximately at minimal
distance to the infeasible point and the evolution strategy is
still able to move. Regarding the minimal distance, there come
different distance definitions into mind, e.g. the $\ell_2$
and the $\ell_1$ norm, respectively.
We use the latter instead of the squared euclidean distance
(squared $\ell_2$ norm).
In addition, we propose another projection method based on random
reference points.

The projection formulated as the
minimization of the squared $\ell_2$ norm leads to a Quadratic Program (QP).
This is, however, computationally expensive. In particular, if the evolution
strategy moves near the boundary of the feasible region
the probability of repair is high.
For this reason minimizing the $\ell_1$ distance for repair
and a new projection approach based on random reference points are investigated
with the goal of having a projection method with a
faster asymptotic runtime.
The runtime scaling behavior of the different projection approaches
has been experimentally compared. The comparison plot is provided in
the supplementary
material (\Cref{chapter:conssatrepairlinear:projmethodruntimecomp}).
Additionally, the projection quality has been experimentally compared
for the different projection approaches
(see \Cref{chapter:conssatrepairlinear:ecdfcomparison}).

\subsubsection{Projection by minimizing the
  \texorpdfstring{$\ell_1$}{\$ell\_1\$} norm}
Projection by minimizing the $\ell_1$ norm
results in a Linear Program (LP) and thus an LP solver can be used.
The optimization problem is
\begin{align}
  \label{sec:algo:subsec:repair:eqn:projection}
  \begin{split}
  \mathbf{\hat{x}} &= \argmin_{\mathbf{x'}}
                         \lVert\mathbf{x'} - \mathbf{x}\rVert_1
                   = \argmin_{\mathbf{x'}}
                   \left(
                       \sum_k{\left|x'_k - x_k\right|}
                   \right)
  \end{split}
\end{align}
where $\mathbf{x}$ is the individual to be repaired.
We introduce the convenience function
\begin{equation}
  \mathbf{\hat{x}} =
      \text{projectToPositiveOrthant}(\mathbf{x}, \mathbf{A}, \mathbf{b})
\end{equation}
returning the solution $\hat{\mathbf{x}}$ of the
problem~\eqref{sec:algo:subsec:repair:eqn:projection}.
Technically,
problem~\eqref{sec:algo:subsec:repair:eqn:projection}
can be turned into an LP

\begin{align}
\begin{split}
  &\mathbf{1}^T\mathbf{z} \rightarrow \text{min!}\\
  \text{s.t. }
  &\mathbf{z} - \mathbf{x'} \ge -\mathbf{x}\\
  &\mathbf{z} + \mathbf{x'} \ge \mathbf{x}\\
  &\mathbf{A}\mathbf{x'} = \mathbf{b}\\
  &\mathbf{x'} \ge \mathbf{0}.
\end{split}
\end{align}
The introduced vector
$\mathbf{z}$ is used to deal with the cases
of the absolute value operator in~\eqref{sec:algo:subsec:repair:eqn:projection}.
If $x'_k - x_k > 0$, then
$-x'_k + x_k < 0$, if $x'_k - x_k < 0$, then
$-x'_k + x_k > 0$ and if $x'_k - x_k = 0$, then
$-x'_k + x_k = 0$. Consequently, the absolute value operator,
the linear constraints, and the non-negativity
constraint are handled.
Depending on the format the LP solver expects as input,
additional slack variable vectors can be introduced to turn
it into an LP in standard form.

\subsubsection{Projection based on random reference points}
We propose a further alternative projection
idea, the ``Iterative Projection''.
The main idea is to create a set of
one or more points
\begin{equation}
  P = \{\mathbf{p}_k | k \in \{1, \ldots, \#\text{points}\},
  \mathbf{A}\mathbf{p}_k = \mathbf{b}, \mathbf{p}_k \ge \mathbf{0}\}
\end{equation}
inside the feasible region once in the beginning.
These points are computed as follows. For each one a random point
in the null space is chosen.
It is then projected by the
$\ell_1$ minimization approach to get $\mathbf{p}_k$
fulfilling the constraints (\Cref{sec:algo:alg:itprojinitp}
shows the pseudo-code).
\begin{algorithm}[t]
  \begin{algorithmic}[1]
    \Function{initReferencePointsForIterativeProjection}
             {$\mathbf{x} \in \mathbb{R}^D$,
               numberOfPoints, $\mathbf{A}, \mathbf{b}$}
    \label{sec:algo:alg:itprojinitp:functiondef}
    \State{$\mathbf{B}^{D \times N} \gets
      \text{orthonormalize}(\text{null}(\mathbf{A}))$}
    \label{sec:algo:alg:itprojinitp:computenullspace}
    \For{$k \gets 1$ \textbf{to} numberOfPoints}
    \label{sec:algo:alg:itprojinitp:fornumpointsbegin}

      \State{$\mathbf{p}_k \gets$}
        \Statex{\qquad\qquad\Call{projectToPositiveOrthant}
        {$\mathcal{U}[-||\mathbf{x}||, ||\mathbf{x}||],
          \mathbf{A}, \mathbf{b}$}}
    \EndFor
    \label{sec:algo:alg:itprojinitp:fornumpointsend}
    \Return{($\{\mathbf{p}_k | k \in \{1, \ldots, \text{numberOfPoints}\}\}$)}
    \label{sec:algo:alg:itprojinitp:return}
    \EndFunction
  \end{algorithmic}
  \caption{Initialization of the set $P$ of reference points
    for the Iterative Projection.}
  \label{sec:algo:alg:itprojinitp}
\end{algorithm}
With this pre-processing in mind, we now consider
a point $\mathbf{x}$ that needs to be repaired. For this point, movement in the
null space towards a randomly chosen $\mathbf{p} \in P$
is possible without violating the linear constraints.
As soon as the positive orthant is reached, the point is considered repaired.
Intuitively, the points in $P$ should be ``far'' inside the feasible
region. This results in a movement that
yields \emph{different points on the boundary} for
\emph{different points outside the positive orthant}.
In other words, the positive orthant should be reached
before getting ``too close'' to $\mathbf{p}$.
More formally, let $\mathbf{d} = \mathbf{p} - \mathbf{x}$ be the direction
of the movement towards $\mathbf{p}$. If the movement's starting point
$\mathbf{x}$
fulfills $\mathbf{A}\mathbf{x} = \mathbf{b}$ movement in the null space
to the positive orthant is possible without violating
the linear constraints.
Since all the negative elements should be zero,
a factor $\alpha$ is necessary to compute
\begin{equation}
  \label{sec:algo:eqn:projectionalpha}
  \mathbf{x}_{\text{projected}} = \mathbf{x} + \alpha\mathbf{d}
\end{equation}
such that
\begin{equation}
  \label{sec:algo:eqn:projectionnonnegative}
  \mathbf{x}_{\text{projected}} \ge \mathbf{0}.
\end{equation}
The projection $\mathbf{x}_{\text{projected}}$ fulfills the
linear equality constraints
\begin{align}
\begin{split}
  \mathbf{A}\mathbf{x}_{\text{projected}} &=
      \mathbf{A}(\mathbf{x} + \alpha\mathbf{d})
        = \mathbf{A}\mathbf{x} + \alpha\mathbf{A}(\mathbf{p} - \mathbf{x})\\

  &= \mathbf{b} + \alpha(\mathbf{b} - \mathbf{b})
        = \mathbf{b}.
\end{split}
\end{align}
Note that the target point $\mathbf{p} \in P$ is located in the
positive orthant and satisfies the linear equality constraints
$\mathbf{A}\mathbf{p}=\mathbf{b}$.
Consequently, there exists an $\alpha$ that
fulfills\footnote{With $\alpha = 1$ we get
  $\mathbf{x}_{\text{projected}} =
  \mathbf{x} + (\mathbf{p} - \mathbf{x}) = \mathbf{p}$.
  And by construction we know that
  $\mathbf{A}\mathbf{p} = \mathbf{b}$ and $\mathbf{p} \ge 0$.}
\Cref{sec:algo:eqn:projectionalpha,sec:algo:eqn:projectionnonnegative}.
One way would be to approach the positive orthant iteratively in the direction
of $\mathbf{d}$ with a small $\alpha$.
But note that the $\alpha$ can also be computed such that all the
negative elements are non-negative after the projection.
The idea is to move towards the chosen reference point with
an $\alpha$
that yields $0$ for the component with the largest deviation from $0$.
This leads to an algorithm with a running time
that is linear in the dimension of the vector.
\Cref{sec:algo:alg:itproj} shows the pseudo-code.
The input is an $\mathbf{x}$ that needs to
be projected and the linear constraint system $\mathbf{A}$ and $\mathbf{b}$
(\Cref{sec:algo:alg:itproj:functiondef}).
The precondition is checked by the assertion in
\Cref{sec:algo:alg:itproj:assert}.
The result $\mathbf{x}_{\text{projected}}$ is initialized with $\mathbf{x}$
(\Cref{sec:algo:alg:itproj:init}) and the direction
vector $\mathbf{d}$ is computed
(\Cref{sec:algo:alg:itproj:d})
using a $\mathbf{p} \in P$ that is chosen randomly according to a uniform
distribution (\Cref{sec:algo:alg:itproj:choosep}).
Then, the worst alpha is computed in
\Crefrange{sec:algo:alg:itproj:alpha}
{sec:algo:alg:itproj:loopend}.
After the loop, the final projected vector is computed in
\Cref{sec:algo:alg:itproj:finalmovement}
using the calculated $\alpha$.
The loop requires $D$ steps. Every statement inside
the loop can be implemented in constant time. This
leads to a running time of $O(D)$.
\begin{algorithm}
  \begin{algorithmic}[1]
    \Function{projectToPositiveOrthantIter}
             {$\mathbf{x} \in \mathbb{R}^D, \mathbf{A}, \mathbf{b}, P$}
    \label{sec:algo:alg:itproj:functiondef}
    \State{\Call{\textbf{assert}}
               {$D > 0 \land \mathbf{A}\mathbf{x} = \mathbf{b}$}}
    \label{sec:algo:alg:itproj:assert}
    \State{$\mathbf{x}_{\text{projected}} \gets \mathbf{x}$}
    \label{sec:algo:alg:itproj:init}
    \State{Choose a $\mathbf{p}$ uniformly at random from $P$}
    \label{sec:algo:alg:itproj:choosep}
    \State{$\mathbf{d} \gets \mathbf{p} - \mathbf{x}$}
    \label{sec:algo:alg:itproj:d}
    \State{$\alpha \gets 0$}
    \label{sec:algo:alg:itproj:alpha}
    \For{$k \gets 1 \textbf{ to } D$}
      \label{sec:algo:alg:itproj:loopstart}
      \If{$(\mathbf{x}_{\text{projected}})_k < 0 \wedge |(\mathbf{d})_k| > 0$}
        \label{sec:algo:alg:itproj:movement1}
        \State{$\alpha \gets \max(\alpha,
          -\frac{(\mathbf{x})_k}{(\mathbf{d})_k})$}
        \label{sec:algo:alg:itproj:movement2}
      \EndIf
      \label{sec:algo:alg:itproj:movement3}
    \EndFor
    \label{sec:algo:alg:itproj:loopend}
    \State{$\mathbf{x}_{\text{projected}} \gets
            \mathbf{x}_{\text{projected}} +
                \alpha\mathbf{d}$}
    \label{sec:algo:alg:itproj:finalmovement}
    \Return{($\mathbf{x}_{\text{projected}}$)}
    \EndFunction
  \end{algorithmic}
  \caption{Iterative Projection (runtime $O(D)$).}
  \label{sec:algo:alg:itproj}
\end{algorithm}

Experimental results for the different projection methods regarding
runtime and projection quality are provided in the supplementary material
(\Cref{sec:expeval:additionalresults:projectionmethodcomparison}).
\Cref{chapter:conssatrepairlinear:projmethodruntimecomp} shows the
scaling behavior of the different projection methods.
\Cref{chapter:conssatrepairlinear:ecdfcomparison} shows the quality of the
different projection methods. For this, \Cref{sec:algo:alg:lccmsaes} was
configured with the different projection methods and run on different
test problems.

\subsection{lcCMSA-ES Pseudo-Code}
\label{sec:algo:subsec:pseudocode}
\Cref{sec:algo:alg:lccmsaes} shows the lcCMSA-ES
in pseudo-code for the optimization problem described in
\Cref{sec:optproblem}.
It makes use of the ideas described in
\Cref{sec:algo:subsec:initcentroid,
  sec:algo:subsec:mutation,sec:algo:subsec:repair}.

An individual is represented as a tuple $\mathfrak{a}$.
It consists of the objective function value $f(\mathbf{x})$,
the parameter vector $\mathbf{x}$ for achieving this
function value $f(\mathbf{x})$, the null space mutation vector $\mathbf{s}$,
the mutation vector $\mathbf{z}$ and the mutation strength $\sigma$.
The best-so-far (bsf) individual is tracked in $\mathfrak{a}_{\text{bsf}}$
and the corresponding generation in $g_\text{bsf}$
(\cref{sec:algo:alg:lccmsaes:bsfinit,
  sec:algo:alg:lccmsaes:bsfupdate,
  sec:algo:alg:lccmsaes:bsfupdate2} of
\Cref{sec:algo:alg:lccmsaes}).

In \cref{sec:algo:alg:lccmsaes:initparams}
of \Cref{sec:algo:alg:lccmsaes}
all the necessary parameters are initialized.
The covariance matrix $\mathbf{C}$ is initialized
to the identity matrix with dimension
of the null space $N$ (\cref{sec:algo:alg:lccmsaes:initC}).
The initial solution is found
by solving the linear system of equations
$\mathbf{A}\mathbf{x} = \mathbf{b}$ for an
$\mathbf{x}_{\text{inh}}$
in \cref{sec:algo:alg:lccmsaes:findinhsol}.
It is then
randomized in \cref{sec:algo:alg:lccmsaes:initcentroid}
as described in \Cref{sec:algo:subsec:initcentroid}.
In case the non-negativity constraint
(\Cref{sec:optproblem:eqn:nonnegative}) is violated,
the initial solution
is repaired (\cref{sec:algo:alg:lccmsaes:repairinhsol1,
  sec:algo:alg:lccmsaes:repairinhsol2,
  sec:algo:alg:lccmsaes:repairinhsol3}).
Next, the generation loop is entered in
\cref{sec:algo:alg:lccmsaes:generationloopstart}.

In every generation $\lambda$ offspring are created
(\cref{sec:algo:alg:lccmsaes:lambdaloopstart,
sec:algo:alg:lccmsaes:offspringsigma,
sec:algo:alg:lccmsaes:offsprings,
sec:algo:alg:lccmsaes:offspringz,
sec:algo:alg:lccmsaes:offspringx,
sec:algo:alg:lccmsaes:offspringrepair1,
sec:algo:alg:lccmsaes:offspringrepair2,
sec:algo:alg:lccmsaes:offspringrepair3,
sec:algo:alg:lccmsaes:offspringrepair4,
sec:algo:alg:lccmsaes:offspringrepair5,
sec:algo:alg:lccmsaes:offspringf,
sec:algo:alg:lccmsaes:offspringindividual,
sec:algo:alg:lccmsaes:bsfupdate,
sec:algo:alg:lccmsaes:lambdaloopend}).
The offspring's mutation strength is sampled from a log-normal distribution
(\cref{sec:algo:alg:lccmsaes:offspringsigma}).
The mutation direction in the null space is sampled from a normal distribution
with the learned covariance and zero
mean (\cref{sec:algo:alg:lccmsaes:offsprings}).
Transformation of this mutation direction in the null space into the
problem space yields the mutation direction in the problem space
(\cref{sec:algo:alg:lccmsaes:offspringz}).
Using this, the new offspring solution is calculated in
\cref{sec:algo:alg:lccmsaes:offspringx}.
It is repaired if
it violates the non-negativity constraint. This is done by projection.
The projection yields a new
solution (\cref{sec:algo:alg:lccmsaes:offspringrepair2}).
From this, the mutation vector
and the mutation vector in the null space are calculated back
(\cref{sec:algo:alg:lccmsaes:offspringrepair3,
  sec:algo:alg:lccmsaes:offspringrepair4}).

The offspring are ranked
according to the order relation ``$\succ$''
(\cref{sec:algo:alg:lccmsaes:rankoffspring})
to update the values $\mathbf{x}$, $\sigma$, and $\mathbf{C}$.
They are updated with the mean values (denoted by $\langle \cdot \rangle$)
of the corresponding quantities of the $\mu$ best individuals in
\cref{sec:algo:alg:lccmsaes:updatex,
  sec:algo:alg:lccmsaes:updatesigma,
  sec:algo:alg:lccmsaes:updateC}.

Since the goal is to minimize $f$ and $f(\mathbf{x})$ is stored
in the individual, the order relation is defined as
\begin{equation}
  \label{sec:algo:eqn:lccmsaes:order}
  \mathfrak{a}_l \succ \mathfrak{a}_m \Leftrightarrow
      f(\tilde{\mathbf{x}}_l) < f(\tilde{\mathbf{x}}_m).
\end{equation}

There are multiple termination
criteria (\cref{sec:algo:alg:lccmsaes:generationloopend}).
The generation loop is terminated if a maximum number of generations
is reached or the $\sigma$ value falls below a threshold.
In addition, the loop is stopped if
the absolute or relative difference of $\mathbf{x}^{(g)}$
and $\mathbf{x}^{(g - G)}$ is below the threshold
$\varepsilon_{\text{abs}}$ or $\varepsilon_{\text{rel}}$, respectively.
Further, if the best-so-far individual has not been updated for the last
$G_\text{lag}$ generations, the generational loop is quit.

The runtime of the generational loop of \Cref{sec:algo:alg:lccmsaes}
is mainly dominated by two computation steps.
The first expensive step is the eigendecomposition in the computation of
$(\sqrt{\mathbf{C}})_{\text{normalized}}$. Second, the offspring generation step
can be bounded as $O(\lambda \cdot t_{\text{proj}})$.
This represents the worst case
assuming every generated offspring has to be repaired. The runtime cost
of the repair step is denoted as $t_{\text{proj}}$. Consequently, assuming
$\lambda=O(D)$, the projection starts to matter if $t_{\text{proj}}$ gets about
asymptotically quadratic in $D$.

Details concerning the covariance matrix $\mathbf{C}$ are explained
in the following two subsections.

\subsubsection{Computation of \texorpdfstring{$\sqrt{\mathbf{C}}$}
  {\$\textbackslash sqrt\{\textbackslash mathbf\{C\}\}\$}}
The covariance matrix $\mathbf{C}$ is symmetric and positive definite.
Therefore, it holds that $\mathbf{C} = \mathbf{M}\mathbf{M}^T$ where
$\mathbf{M} = \sqrt{\mathbf{C}}$.
The computation of $(\sqrt{\mathbf{C}})_{\text{normalized}}$
for \cref{sec:algo:alg:lccmsaes:sqrtCnormalized}
can be done every $\lfloor{\tau_c\rfloor}$-th generation to save time.
This is possible because the changes to the covariance matrix are small
in between these generations.
\Cref{sec:algo:alg:sqrtC} outlines the $\sqrt{\mathbf{C}}$ calculation steps.
Note that
$\text{det}(\mathbf{M}_r)^{-\frac{1}{N}}$ is a normalization factor such that
the determinant of $(\sqrt{\mathbf{C}})_{\text{normalized}}$ is one.
The idea behind this is that the resulting
transformation is volume-preserving.

\subsubsection{Regularization of \texorpdfstring{$\mathbf{C}$}
  {\$\textbackslash mathbf\{C\}\$}
  for Computing
  \texorpdfstring{$\sqrt{\mathbf{C}}$}
                 {\$\textbackslash sqrt\{\textbackslash mathbf\{C\}\}\$}}
\label{sec:algo:subsec:regularization}
\mbox{}\\
When the strategy approaches the boundary, the selected (and repaired)
mutation steps toward
the boundary decrease rapidly. But the other directions
are not affected. Consequently, the condition number of $\mathbf{C}$
increases rapidly.

Therefore, regularization of $\mathbf{C}$ to delimit the condition
number is a way to overcome this.
This prevents the ES from evolving in a degenerated
subspace of the null space when approaching
the boundary. The regularization is done by
adding a small positive value to the diagonal
elements if the condition number exceeds a threshold $t$, i.e.,
\begin{equation}
    \mathbf{M}_r = \sqrt{\mathbf{C}} + r\mathbf{I} \text{ with }
    r = 0 \text{ if } \text{cond}(\mathbf{C}) \le t.
\end{equation}

The regularized covariance matrix $\tilde{\mathbf{C}}$ is then
$\mathbf{M}_r\mathbf{M}_r^T$.
Let $\lambda_i$ denote the $i$-th eigenvalue\footnote{Note
that we use $\lambda_i$ here to denote an eigenvalue. In
\Cref{sec:algo:alg:lccmsaes} we use $\lambda$ to denote
the number of offspring.}
of $\mathbf{C}$
such that $\lambda_1 \le \lambda_i \le \lambda_N$.
Accordingly, the $i$-th eigenvalue of
$\sqrt{\mathbf{C}}$ is $\sqrt{\lambda_i}$.
The eigenvalues of $\mathbf{M}_r$ and
$\mathbf{M}_r\mathbf{M}_r^T$ are
$\sqrt{\lambda_i} + r$ and $(\sqrt{\lambda_i} + r)^2$,
respectively.

In case the condition number exceeds the threshold $t$,
i.e.,
$\text{cond}(\mathbf{C}) = \text{cond}(\mathbf{M}\mathbf{M}^T)
=\lambda_N/\lambda_1 > t$, the
factor $r$ is chosen to limit the condition number to $t$.
That is, the corresponding $r$ value is determined by
\begin{align}
  \label{sec:algo:eqn:targetcondition}
  \begin{split}
    \text{cond}(\tilde{\mathbf{C}})
    & = \text{cond}(\mathbf{M}_r\mathbf{M}_r^T)
      = \frac{(\sqrt{\lambda_N} + r)^2}{(\sqrt{\lambda_1} + r)^2}
    \stackrel{!}{=} t
  \end{split}
\end{align}
The detailed steps solving \Cref{sec:algo:eqn:targetcondition}
for $r$ are provided in the supplementary material
(\Cref{sec:appendix:regularizationofsqrtC}). They result in
\begin{equation}
    r = \frac{\sqrt{\lambda_N}}{t} - \sqrt{\lambda_1}
      + \sqrt{\frac{\lambda_N}{t^2} + \frac{\lambda_N}{t}
        - \frac{2\sqrt{\lambda_1\lambda_N}}{t}}.
\end{equation}

\begin{algorithm}[t]
  \begin{algorithmic}[1]
    \State{Input $\mathbf{A}, \mathbf{b}, f$}
    \State{Initialize parameters
      $\mu$, $\lambda$, $\sigma$, $\tau$, $\tau_c$,
      $G$, $G_{\text{lag}}$, $g_{\text{stop}}$, $\sigma_{\text{stop}}$,
      $\epsilon_{\text{abs}}$, $\epsilon_{\text{rel}}$, $t$}
    \label{sec:algo:alg:lccmsaes:initparams}
    \State{$\mathbf{B}^{D \times N} \gets
      \text{orthonormalize}(\text{null}(\mathbf{A}))$}
    \label{sec:algo:alg:lccmsaes:computenullspace}
    \State{$\mathbf{C} \gets \mathbf{I}^{N \times N}$}
    \label{sec:algo:alg:lccmsaes:initC}
    \State{$\mathbf{x}^{(0)} \gets$
      \Call{findInhomogeneousSolution}{$\mathbf{A}, \mathbf{b}$}}
    \label{sec:algo:alg:lccmsaes:findinhsol}
    \State{$P \gets$ \Call{initReferencePointsForIterativeProjection}
      {$\mathbf{x}^{(0)}$, $10N$, $\mathbf{A}$, $\mathbf{b}$}}
    \label{sec:algo:alg:lccmsaes:initp}
    \State{Randomize $\mathbf{x}^{(0)}$,}
    \Statex{\qquad e.g.: $\mathbf{x}^{(0)} \gets
        \mathbf{x}^{(0)} + ||\mathbf{x}^{(0)}||
        \mathbf{B}\mathcal{N}(\mathbf{0}, \mathbf{I}^{N \times N})$}
    \label{sec:algo:alg:lccmsaes:initcentroid}
    \If{$\left(\mathbf{x}^{(0)}\right)_{1:D} < 0$}
    \label{sec:algo:alg:lccmsaes:repairinhsol1}
    \State{$\mathbf{x}^{(0)} \gets$ \Call{projectToPositiveOrthantIter}
      {$\mathbf{x}^{(0)}$, $\mathbf{A}$, $\mathbf{b}$, $P$}}
    \label{sec:algo:alg:lccmsaes:repairinhsol2}
    \EndIf
    \label{sec:algo:alg:lccmsaes:repairinhsol3}
    \State{$\left(\mathfrak{a}_{\text{bsf}}, g_{\text{bsf}}\right)
      \gets
      \left((f(\mathbf{x}^{(0)}), \mathbf{x}^{(0)},
      \mathbf{0}, \mathbf{0},
      \sigma),
      0
      \right)$}
    \label{sec:algo:alg:lccmsaes:bsfinit}
    \State{$g \gets 0$}
    \label{sec:algo:alg:lccmsaes:generationcounterinit}
    \Repeat
    \label{sec:algo:alg:lccmsaes:generationloopstart}
    \State{$(\sqrt{\mathbf{C}})_{\text{normalized}} \gets$
      \Call{computeSqrtCNormalized}{$\mathbf{C}$, $t$}}
    \label{sec:algo:alg:lccmsaes:sqrtCnormalized}
    \For{$l \gets 1 \textbf{ to } \lambda$}
    \label{sec:algo:alg:lccmsaes:lambdaloopstart}
    \State{$\tilde{\sigma}_l \gets
      \sigma e^{\tau \mathcal{N}_l(0, 1)}$}
    \label{sec:algo:alg:lccmsaes:offspringsigma}
    \State{$\tilde{\mathbf{s}}_l \gets
      (\sqrt{\mathbf{C}})_{\text{normalized}}
      \mathcal{N}_l(\mathbf{0},
      \mathbf{I}^{N \times N})$}
    \label{sec:algo:alg:lccmsaes:offsprings}
    \State{$\tilde{\mathbf{z}}_l \gets
      \tilde{\sigma}_l \mathbf{B} \tilde{\mathbf{s}}_l$}
    \label{sec:algo:alg:lccmsaes:offspringz}
    \State{$\tilde{\mathbf{x}}_l \gets
      \mathbf{x}^{(g)} + \tilde{\mathbf{z}}_l$}
    \label{sec:algo:alg:lccmsaes:offspringx}
    \If{$(\tilde{\mathbf{x}}_l)_{1:D} < 0$}
    \label{sec:algo:alg:lccmsaes:offspringrepair1}
    \State{$\tilde{\mathbf{x}}_l \gets$}
      \Statex{\qquad\qquad\qquad\Call{projectToPositiveOrthantIter}
      {$\tilde{\mathbf{x}}_l$, $\mathbf{A}$, $\mathbf{b}$, $P$}}
    \label{sec:algo:alg:lccmsaes:offspringrepair2}
    \State{$\tilde{\mathbf{z}}_l \gets \tilde{\mathbf{x}}_l -
      \mathbf{x}^{(g)}$}
    \label{sec:algo:alg:lccmsaes:offspringrepair3}
    \State{$\tilde{\mathbf{s}}_l \gets \mathbf{B}^T
      \tilde{\mathbf{z}}_l / \tilde{\sigma}_l$}
    \label{sec:algo:alg:lccmsaes:offspringrepair4}
    \EndIf
    \label{sec:algo:alg:lccmsaes:offspringrepair5}
    \State{$\tilde{f}_l \gets f(\tilde{\mathbf{x}}_l)$}
    \label{sec:algo:alg:lccmsaes:offspringf}
    \State{$\tilde{\mathfrak{a}}_l \gets (\tilde{f}_l, \tilde{\mathbf{x}}_l,
      \tilde{\mathbf{z}}_l, \tilde{\mathbf{s}}_l, \tilde{\sigma}_l)$}
    \label{sec:algo:alg:lccmsaes:offspringindividual}
    \State{$\left(\mathfrak{a}_{\text{bsf}}, g_{\text{bsf}}\right)
      \gets
      \begin{cases}
        \left(\tilde{\mathfrak{a}}_l, g + 1\right) &
        \text{if } \tilde{\mathfrak{a}}_l \succ
        \mathfrak{a}_{\text{bsf}} \\
        \left(\mathfrak{a}_{\text{bsf}}, g_\text{bsf}\right) &
        \text{otherwise}
      \end{cases}$}
    \label{sec:algo:alg:lccmsaes:bsfupdate}
    \EndFor
    \label{sec:algo:alg:lccmsaes:lambdaloopend}
    \State{\Call{rankOffspringPopulation}{$\tilde{\mathfrak{a}}_1,
      \ldots, \tilde{\mathfrak{a}}_\lambda$}}
    \Statex{\hspace{4cm}acc. to ``$\succ$''(\Cref{sec:algo:eqn:lccmsaes:order})}
    \label{sec:algo:alg:lccmsaes:rankoffspring}
    \State{$\mathbf{x}^{(g + 1)} \gets \mathbf{x}^{(g)} +
      \langle \tilde{\mathbf{z}} \rangle$}
    \label{sec:algo:alg:lccmsaes:updatex}
    \State{$\mathfrak{a} \gets
      \left(f(\mathbf{x}^{(g + 1)}), \mathbf{x}^{(g + 1)},
      \langle\tilde{\mathbf{z}}\rangle,
      \langle\tilde{\mathbf{s}}\rangle,
      \langle\tilde{\sigma}\rangle
      \right)$}
    \label{sec:algo:alg:lccmsaes:centroidindividual}
    \State{$\left(\mathfrak{a}_{\text{bsf}}, g_{\text{bsf}}\right)
      \gets
      \begin{cases}
        \left(\mathfrak{a}, g + 1\right) &
        \text{if } \mathfrak{a} \succ
        \mathfrak{a}_{\text{bsf}} \\
        \left(\mathfrak{a}_{\text{bsf}}, g_\text{bsf}\right) &
        \text{otherwise}
      \end{cases}$}
    \label{sec:algo:alg:lccmsaes:bsfupdate2}
    \State{$\sigma \gets \langle \tilde{\sigma} \rangle$}
    \label{sec:algo:alg:lccmsaes:updatesigma}
    \State{$\mathbf{C} \gets \left(1 - \frac{1}{\tau_c}\right)\mathbf{C} +
      \frac{1}{\tau_c}
      \langle\tilde{\mathbf{s}}\tilde{\mathbf{s}}^T\rangle$}
    \label{sec:algo:alg:lccmsaes:updateC}
    \State{$g \gets g + 1$}
    \label{sec:algo:alg:lccmsaes:generationcounterupdate}
    \Until{$g > g_{\text{stop}} \vee \sigma < \sigma_{\text{stop}}
      \vee ||\mathbf{x}^{(g)} - \mathbf{x}^{(g - G)}||< \epsilon_{\text{abs}}
      \vee \left|\frac{||\mathbf{x}^{(g)}||}{||\mathbf{x}^{(g - G)}||}
      - 1\right| < \epsilon_{\text{rel}}
      \vee g - g_\text{bsf} \ge G_\text{lag}$}
    \label{sec:algo:alg:lccmsaes:generationloopend}
  \end{algorithmic}
  \caption{The $(\mu/\mu_I, \lambda)$-lcCMSA-ES.}
  \label{sec:algo:alg:lccmsaes}
\end{algorithm}

\begin{algorithm}
  \begin{algorithmic}[1]
    \Function{computeSqrtCNormalized}{$\mathbf{C}$, $t$}
    \State{$\mathbf{C} \gets
      \frac{1}{2} \left(\mathbf{C} + \mathbf{C}^T\right)$}
    \State{\parbox[t]{\dimexpr\linewidth-\algorithmicindent}
      {Perform eigendecomposition to get
        $\mathbf{U}$ and $\mathbf{D}$ such that
        $\mathbf{C} = \mathbf{U}\mathbf{D}\mathbf{U}^T$
        with $\mathbf{D}$ being the diagonal matrix of eigenvalues
        and the columns of $\mathbf{U}$ being the corresponding
        eigenvectors}}
    \State{$(\lambda_1, \ldots, \lambda_N)^T \gets \text{diag}(\mathbf{D})$}
    \State{$r \gets 0$}
    \If{$\text{cond}(\mathbf{C}) > t$}
      \State{$r =
        \frac{\sqrt{\lambda_N}}{t} - \sqrt{\lambda_1}
        + \sqrt{\frac{\lambda_N}{t} + \frac{\lambda_N}{t^2}
        - \frac{2\sqrt{\lambda_1\lambda_N}}{t}}$}
    \EndIf
    \State{$\mathbf{M}_r \gets \mathbf{U}\sqrt{\mathbf{D}} + r\mathbf{I}$}
    \State{$(\sqrt{\mathbf{C}})_{\text{normalized}}
      = \text{det}(\mathbf{M}_r)^{-\frac{1}{N}}\mathbf{M}_r$}
    \Return{($(\sqrt{\mathbf{C}})_{\text{normalized}}$)}
    \EndFunction
  \end{algorithmic}
  \caption{Computation of
    $(\sqrt{\mathbf{C}})_{\text{normalized}}$.}
  \label{sec:algo:alg:sqrtC}
\end{algorithm}

\section{Experimental Evaluation}
\label{sec:expeval}
The lcCMSA-ES is tested on a variety of different test functions.
Linear objective functions are considered as a first step and
non-linear objective functions as a second step.
For the linear objective function tests,
the Klee-Minty cube~\cite{KleeMinty1972} is used.
For the non-linear objective function experiments, the BBOB
COCO framework~\cite{cocodoc} with adaptions is used.
In addition, the performance of the lcCMSA-ES is compared
with other methods that are able to deal with problem
\eqref{sec:optproblem:eqn:problem}.

The algorithms are implemented in
Octave with mex-extensions\footnote{We provide the
code in a GitHub repository (\url{https://github.com/patsp/lcCMSA-ES}).}
and the experiments were run on a cluster with 5 nodes.
Every node has an Intel 8-core Xeon E5420 2.50GHz processor
with 8GiB of RAM running a GNU/Linux system.
For the BBOB COCO tests, the post-processing tool with slight adjustments
was used to generate the figures. This post-processing tool is
part of the BBOB COCO framework.

In the experiments, the parameters for the lcCMSA-ES are set as shown in
\Cref{chapter:conssatrepairlinear:parameters}.
The six parameters $G$, $G_{\text{lag}}$,
$g_{\text{stop}}$, $\sigma_{\text{stop}}$
$\epsilon_{\text{abs}}$, and $\epsilon_{\text{rel}}$ are used
for the stopping criteria. The chosen values turned out to be good
choices in initial experiments. The initial $\sigma$, $\tau$, and $\tau_c$
were set according to the suggestions in~\cite{Beyer2008}.
The parameters $\mu$ and $\lambda$ were chosen to have a truncation ratio
$\mu/\lambda=1/4$ (similar as in~\cite{Beyer2008}). The value of $t$ was
set as a trade-off between numerical accuracy and the toleration of approaching
the boundary in the ES.

The sum of objective and constraint function evaluations are
considered for the performance measure.
In the BBOB COCO framework one call to the constraint
evaluation function yields the values of all the constraints for a given
query point.

\begin{table}
  \caption{Parameter settings for the lcCMSA-ES experiments.}
  \label{chapter:conssatrepairlinear:parameters}
  \centering
  \ra{1.3}
  \begin{tabular}{lr|lr}
    \hline
    \multicolumn{2}{c|}{Core ES param.} & \multicolumn{2}{c}{Stopping criteria param.}\\
    \hline
    $\lambda$ & $4D$ & $G$ & $10$\\
    $\mu$ & $\lfloor \frac{\lambda}{4} \rfloor$ & $G_{\text{lag}}$ & $50N$\\
    $\sigma$ (initial value) & $\frac{1}{\sqrt{D}}$ & $g_{\text{stop}}$ & $10000$\\
    $\tau$ & $\frac{1}{\sqrt{2N}}$ & $\sigma_{\text{stop}}$ & $10^{-6}$\\
    $\tau_c$ & $1 + \frac{N(N-1)}{2\mu}$ & $\epsilon_{\text{abs}}$ & $10^{-9}$\\
    $t$ & $10^{12}$ & $\epsilon_{\text{rel}}$ & $10^{-9}$\\
    \hline
  \end{tabular}
\end{table}

\subsection{Performance on the Klee-Minty cube}
Klee and Minty formulated a special LP~\cite{KleeMinty1972}
to show that the Simplex
algorithm~\cite{Dantzig1951Simplex}, although working well in practice,
has an exponential runtime in the worst case. To this end,
they invented the so-called Klee-Minty cube.
The $n$-dimensional Klee-Minty cube is a distorted hypercube with $2^n$
corners. The inside of
the cube represents the feasible region. The objective function is constructed
in such a way
that the Simplex algorithm visits all the corners in the worst case and thus
its worst case runtime is exponential.
Formally, the inequalities of the feasible region write
\begin{equation}
  \begin{array}{r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r@{\hspace{3pt}}r}
    x_1    &   &                         & &&&&&&       \le & 5 \\
    4x_1   & + & x_2                       &&&&&&&          \le & 25 \\
    \vdots &   & \vdots             & &&&&&& \vdots & \vdots  \\
    2^nx_1 & +  & 2^{(n-1)}x_2& + & \cdots & + &
                         4x_{n-1} & + & x_n & \le& 5^n\\
  \end{array}
\end{equation}
where $x_1 \ge 0, \ldots, x_n \ge 0$.
The objective function is
$2^{(n-1)}x_1 + 2^{(n-2)}x_2 + \cdots + 2x_{n-1} + x_n \rightarrow \text{max!}$.
The maximum is reached for the vector
$\mathbf{x}_{\text{opt}} = \left(0, 0, \ldots, 0, 5^n\right)^T$
yielding $f(\mathbf{x}_{\text{opt}}) = 5^n$.
The Klee-Minty problem has been chosen because it is known to be
also a hard problem for interior point
methods~\cite{Megiddo1989Boundary,Deza2008HowGoodInterior}.

\Cref{sec:expeval:resultstable} shows the results of
single runs of the lcCMSA-ES with the Iterative Projection
on the Klee-Minty problem with different dimensions.
The optimal value is reached
up to a small error for all the dimensions from $1$ to $15$.
For dimensions larger than $15$ we have observed numerical
instabilities.

We also tested interior point LP solvers on
the Klee-Minty problem. We applied glpk's interior
point LP algorithm using Octave and Mathematica's interior point LP algorithm
to the Klee-Minty problem. We have observed
that the absolute error to the optimum increases with increasing dimension
for both LP solvers. The supplementary material contains
detailed results (\Cref{sec:expeval:additionalresults:kleeminty}).
\Cref{sec:expeval:resultstableglpkoctave,
sec:expeval:resultstablemathematica} display the results of
single runs of the glpk LP solver and the Mathematica LP solver, respectively.
Both solvers were run with default parameters and interior point
methods. The number of generations and the number of function evaluations
are not comparable. For example, according to the documentation,
the default number of maximum iterations for the interior point algorithm
glpk in Octave is 200. This is independent of the number of variables
and constraints.

\begin{table*}
  \caption{Results of single runs of the lcCMSA-ES with the
    Iterative Projection (linear runtime version)
    on the Klee-Minty cube.}
  \label{sec:expeval:resultstable}
  \centering
\begin{tabular}{l|l|l|l|l|l|l}
  Name & $f_{\text{opt}}$ & $\text{ES}_{f_{\text{best}}}$ & $|f_{\text{opt}} - \text{ES}_{f_{\text{best}}}|$ & $|f_{\text{opt}} - \text{ES}_{f_{\text{best}}}| / |f_{\text{opt}}|$ & \#generations & \#f-evals \\
  \hline
  Klee-Minty $D = 1$ & -5.000000 & -5.000000 & 2.910383e-11 & 5.820766e-12 & 97 & 874 \\
  Klee-Minty $D = 2$ & -25.000000 & -25.000000 & 2.693810e-10 & 1.077524e-11 & 104 & 1769 \\
  Klee-Minty $D = 3$ & -125.000000 & -125.000000 & 1.987161e-09 & 1.589729e-11 & 153 & 3826 \\
  Klee-Minty $D = 4$ & -625.000000 & -625.000000 & 2.280285e-08 & 3.648456e-11 & 201 & 6634 \\
  Klee-Minty $D = 5$ & -3125.000000 & -3125.000000 & 2.121087e-07 & 6.787479e-11 & 251 & 10292 \\
  Klee-Minty $D = 6$ & -15625.000000 & -15625.000003 & 2.568122e-06 & 1.643598e-10 & 301 & 14750 \\
  Klee-Minty $D = 7$ & -78125.000000 & -78125.000030 & 3.049150e-05 & 3.902912e-10 & 351 & 20008 \\
  Klee-Minty $D = 8$ & -390625.000000 & -390625.000303 & 3.030710e-04 & 7.758617e-10 & 403 & 26196 \\
  Klee-Minty $D = 9$ & -1953125.000000 & -1953125.001656 & 1.656145e-03 & 8.479462e-10 & 451 & 32924 \\
  Klee-Minty $D = 10$ & -9765625.000000 & -9765625.000914 & 9.139776e-04 & 9.359131e-11 & 501 & 40582 \\
  Klee-Minty $D = 11$ & -48828125.000000 & -48828125.000000 & 0.000000e+00 & 0.000000e+00 & 551 & 49040 \\
  Klee-Minty $D = 12$ & -244140625.000000 & -244140625.000000 & 2.980232e-08 & 1.220703e-16 & 602 & 58395 \\
  Klee-Minty $D = 13$ & -1220703125.000000 & -1220703125.000000 & 0.000000e+00 & 0.000000e+00 & 650 & 68251 \\
  Klee-Minty $D = 14$ & -6103515625.000000 & -6103515625.000001 & 9.536743e-07 & 1.562500e-16 & 735 & 83056 \\
  Klee-Minty $D = 15$ & -30517578125.000000 & -30517578125.000004 & 3.814697e-06 & 1.250000e-16 & 755 & 91356 \\
\end{tabular}
\end{table*}

\subsection{Performance on the BBOB COCO constrained suite}
For the non-linear objective function experiments the BBOB
COCO
framework~\cite{cocodoc} with adaptions is used.
The adapted version\footnote{We provide the adapted code in a GitHub fork
of the BBOB COCO framework, \url{https://github.com/patsp/coco}.
The changes are in the new branch \texttt{development-sppa-2}.
This branch is based on the \texttt{development} branch of
\url{https://github.com/numbbo/coco} with changes up to and including
Dec 10, 2017.
A list of the changes is also provided in the supplementary material
(\Cref{sec:appendix:coco_changes}).}
is based on the code in the
branch
\texttt{development}\footnote{Because the \texttt{bbob-constrained} suite
is still under development, we provide a fork. This makes our results
reproducible. Even though it is still under development,
this suite gives a good indication of the algorithm's performance
in comparison to other methods. We use this suite
instead of defining our own test problems with linear constraints for this
work.}
in~\cite{cococoderepo}.
A documentation can be found in~\cite{cocodocrepo} under
\path{docs/bbob-constrained/functions/build}
after building it according to the instructions.

The BBOB COCO framework
provides a test suite, bbob-constrained, for constrained
black-box optimization
benchmarking. It contains 48 constrained functions with
dimensions $D \in \{2, 3, 5, 10, 20, 40\}$. For every problem, random instances
can be generated. The 48 problems are constructed by combining
8 functions of the standard BBOB COCO
suite for single-objective optimization with 6 different numbers of
constraints, namely $1$, $2$, $6$, $6 + D / 2$, $6 + D$, and $6 + 3D$
constraints.
The 8 functions are Sphere, Separable Ellipsoid, Linear Slope,
Rotated Ellipsoid, Discus, Bent Cigar, Sum of Different Powers, and the
Separable Rastrigin.
The constraints are linear with nonlinear perturbations
and defined by their gradient. These constraints are generated
by sampling their gradient vectors from a normal distribution
and ensuring that the feasible region is not empty.
The generic algorithm of generating a constrained problem
is outlined in~\cite{cocodocrepo},
\path{docs/bbob-constrained/functions/build}.

The optimization problem in the BBOB COCO
framework is stated as
\begin{subequations}
  \label{sec:expeval:eqn:bbobprob}
  \begin{align}
    &f'(\mathbf{x}) \rightarrow \text{min!}\\
    \text{s.t. }& \mathbf{g}(\mathbf{x}) \le \mathbf{0}\\
    & \check{\mathbf{x}} \le \mathbf{x} \le \hat{\mathbf{x}}
  \end{align}
\end{subequations}
where $f' : \mathbb{R}^{D'} \rightarrow \mathbb{R}$
and $\mathbf{g} : \mathbb{R}^{D'} \rightarrow \mathbb{R}^{K'}$.
In order for the lcCMSA-ES to be applicable
this must be transformed into
\begin{subequations}
  \label{sec:expeval:eqn:bbobprobtransformed}
  \begin{align}
    &f(\mathbf{x}) \rightarrow \text{min!}\\
    \text{s.t. }& \mathbf{A}\mathbf{x} = \mathbf{b}\\
    & \mathbf{x} \ge \mathbf{0}
  \end{align}
\end{subequations}
where $f: \mathbb{R}^D \rightarrow \mathbb{R}$,
$\mathbf{A} \in \mathbb{R}^{K \times D}$, $\mathbf{x} \in \mathbb{R}^D$,
$\mathbf{b} \in \mathbb{R}^K$.
It is known that the constraints in the bbob-constrained
suite of the BBOB COCO framework are linear
with non-linear perturbations.
Using this fact in addition with the enhanced ability to disable the
non-linear perturbations, a pre-processing step is used.
It transforms \Cref{sec:expeval:eqn:bbobprob} into
\Cref{sec:expeval:eqn:bbobprobtransformed}. This pre-processing
step is based on the idea
of querying the constraint function at enough positions in the parameter
space.
This allows constructing a system of equations
that can be solved for the underlying
coefficients of the linear constraints. The resulting coefficients and
the bounds can be put into matrix form. Slack variables are added for
transforming the inequalities into equalities to arrive at the form
in \Cref{sec:expeval:eqn:bbobprobtransformed}.
Due to space limitations,
\Cref{sec:preprocessbbobprob,
chapter:conssatrepairnonlinear:sec:approxunknownconstraintsglobally}
in the supplementary material
describe how this can be done and show pseudo-code.

In the following, the performance of different algorithms
is visualized by use of bootstrapped Empirical Cumulative Distribution
Functions (ECDF). These plots show
the percentages of function target values reached for a given
budget of function and constraint evaluations
per search space dimensionality. The x-axis shows the sum of
objective function and constraint evaluations normalized by dimension
(log-scaled).
The y-axis shows the percentage of so-called targets that were reached
for the given sum of objective function and constraint evaluations.
Every target is defined as a particular distance from the optimum.
In the plots shown, the
standard BBOB ones are used:
$f_{\text{target}} = f_{\text{opt}} + 10^k$ for 51 different values of
k between $-8$ and $2$.
The crosses indicate the medians of the sum of objective function and
constraint evaluations of instances that did not reach the most
difficult target.
Note that the steps at the beginning of the lines of
some variants are due to the pre-processing step that requires
an initial amount of constraint evaluations.
Furthermore, the top-left corner in every plot shows information
about the experiments. The first line indicates the functions
of the BBOB COCO framework that have been used in the experiment.
The second line specifies the targets. The number of runs (instances)
are indicated in the third line. A line (with a marker) to every entry in the
legend is drawn. This shows which line in the plot belongs to which entry
in the legend.

The ECDF plots shown in
\Cref{sec:expeval:lccmsaes:ecdf,
sec:expeval:compwithother:ecdfbbobcococomparisonall} show results
aggregated over multiple problems. For
this, the results of an algorithm over all the problems are considered for a
specific dimension. It is referred to~\cite{Hansen2016CocoPerformance}
for all the details. In the supplementary material, we provide single function
plots (\Cref{chapter:conssatrepairlinear:ecdf}).

\Cref{sec:expeval:lccmsaes:ecdf} presents the ECDF
of runs of the lcCMSA-ES with
the Iterative Projection. All the constrained problems of the BBOB COCO
bbob-constrained test suite are shown aggregated. The considered dimensions
are 2, 3, 5, 10, 20, 40. For these, the performance of the algorithm
is evaluated on 15 independent randomly generated
instances of each constrained test problem. Based
on the observed run lengths, ECDF
graphs are generated. Every line corresponds to the aggregated
ECDF over all problems of a dimension
(2, 3, 5, 10, 20, 40 from top to bottom).

Additional simulation results of the lcCMSA-ES
are presented in the supplementary material
(\Cref{sec:expeval:additionalresults:lccmsaes}).
They include ECDF and average runtime graphs of the lcCMSA-ES
for all the single functions of the BBOB COCO constrained suite.
Graphs showing the evolution dynamics of the lcCMSA-ES are presented as well.

We see that for the dimensions 2, 3, and 5 the most difficult target
is reached with about $10^5D$
function and constraint evaluations.
For the higher dimensions the most difficult target is not reached. But
about 90\% of the targets are reached with about $10^6D$
function and constraint evaluations. One can see that the
performance in the higher dimensions
is low in particular for the Rastrigin functions (functions 43-48 shown
in the last six subplots of \Cref{chapter:conssatrepairlinear:ecdf}).
The Rastrigin function is multimodal. In order to deal with such a function,
an extension of the lcCMSA-ES is possible. An example could be an integration
of the lcCMSA-ES into a restart meta ES.

\providecommand{\bbobdatapath}{}
\renewcommand{\bbobdatapath}{figures/bbob_constrained_coco_2/linear_constraints/itprojlccmsaes/cocoBenchmark/ppdata/}
\subimport*{\bbobdatapath}{cocopp_commands.tex}
\providecommand{\pprldmanydatapath}{}
\renewcommand{\pprldmanydatapath}{\bbobdatapath itprojlccmsaes_on_bbob-constrained_f1_48/pprldmany-single-functions/}
\providecommand{\ppfigdimdatapath}{}
\renewcommand{\ppfigdimdatapath}{\bbobdatapath itprojlccmsaes_on_bbob-constrained_f1_48/}
\providecommand{\DIM}{}
\renewcommand{\DIM}{\ensuremath{\mathrm{DIM}}}
\providecommand{\aRT}{}
\renewcommand{\aRT}{\ensuremath{\mathrm{aRT}}}
\providecommand{\FEvals}{}
\renewcommand{\FEvals}{\ensuremath{\mathrm{FEvals}}}
\providecommand{\nruns}{}
\renewcommand{\nruns}{\ensuremath{\mathrm{Nruns}}}
\providecommand{\Dfb}{}
\renewcommand{\Dfb}{\ensuremath{\Delta f_{\mathrm{best}}}}
\providecommand{\Df}{}
\renewcommand{\Df}{\ensuremath{\Delta f}}
\providecommand{\nbFEs}{}
\renewcommand{\nbFEs}{\ensuremath{\mathrm{\#FEs}}}
\providecommand{\fopt}{}
\renewcommand{\fopt}{\ensuremath{f_\mathrm{opt}}}
\providecommand{\ftarget}{}
\renewcommand{\ftarget}{\ensuremath{f_\mathrm{t}}}
\providecommand{\CrE}{}
\renewcommand{\CrE}{\ensuremath{\mathrm{CrE}}}
\providecommand{\change}[1]{}
\renewcommand{\change}[1]{{\color{red} #1}}

\begin{figure}
  \centering
  \includegraphics[width=0.35\textwidth]{\pprldmanydatapath pprldmany}
  \caption{Bootstrapped empirical cumulative distribution
    function of the number
    of objective function and constraint evaluations
    divided by dimension for the lcCMSA-ES with the
    Iterative Projection.
    \sppabbobecdfaxisexplanation}
  \label{sec:expeval:lccmsaes:ecdf}
\end{figure}

\subsection{Comparison with other approaches}
To compare the lcCMSA-ES proposed in this work
a selection of other algorithms is benchmarked
on the same adapted BBOB COCO suite.
Three variants of DE that showed promising results in benchmarks
are tested, namely
``Self-adaptive Differential Evolution Algorithm for Constrained
Real-Parameter Optimization'' (conSaDE)~\cite{Huang2006ConSADE},
``Differential Evolution with Ensemble of Constraint Handling
Techniques'' (ECHT-DE)~\cite{MalliPeddi2010DE},
and
``Constrained Optimization by the $\varepsilon$ Constrained
Differential Evolution with an Archive and Gradient-Based
Mutation'' ($\varepsilon$DEag)~\cite{Takahama2010EpsilonDEGA}.
Further, an Active-Set ES~\cite{Arnold2016ActiveSetES},
an ES with augmented Lagrangian constraint
handling~\cite{Atamna2017LagrangianES} and a method
based on surrogate
modeling with adaptive parameter control~\cite{Bagheri2015SACOBRA}
(SACOBRA) are benchmarked and compared
to the approach presented in this work.

For the conSaDE, the ECHT-DE, the $\varepsilon$DEag, the Active-Set ES,
the ES with augmented Lagrangian constraint
handling, and the SACOBRA algorithm, the implementations provided by
the respective authors were
used\footnote{conSaDE: \url{http://web.mysites.ntu.edu.sg/epnsugan/PublicSite/Shared\%20Documents/Codes/2006-CEC-Const-SaDE.rar},
ECHT-DE: \url{http://web.mysites.ntu.edu.sg/epnsugan/PublicSite/Shared\%20Documents/Codes/2010-TEC-Ens-Con-EP-DE.zip},
$\varepsilon$DEag: \url{http://www.ints.info.hiroshima-cu.ac.jp/~takahama/download/eDEa-2010.0430.tar.gz},
Active-Set ES: \url{https://web.cs.dal.ca/~dirk/AS-ES.tar},
SACOBRA: \url{https://cran.r-project.org/web/packages/SACOBRA/index.html},
ES with augmented Lagrangian constraint handling: Code provided by
Asma Atamna.} (adapted for the BBOB COCO framework).
All the algorithms for the comparison were run with default parameters.

The goal is to have a direct comparison to the method proposed in this work.
The non-linear
transformations are turned off in order to be able to compare
the approaches to the lcCMSA-ES.
For the lcCMSA-ES the
inequality constraints are first pre-processed
to transform the problem into one that the lcCMSA-ES is able to handle.
Hence, the lcCMSA-ES
solves~\eqref{sec:expeval:eqn:bbobprobtransformed}.
Similarly, for the active-set ES, the linear constraints are
determined but no slack variables are added because inequalities
can be handled by the algorithm. The optimization problem
\begin{subequations}
  \label{sec:expeval:eqn:bbobprobtransformedforactivesetes}
  \begin{align}
    &f(\mathbf{x}) \rightarrow \text{min!}\\
    \text{s.t. }& \mathbf{A}\mathbf{x} \le \mathbf{b}\\
    &\check{\mathbf{x}} \le \mathbf{x} \le \hat{\mathbf{x}}
  \end{align}
\end{subequations}
is passed to the active-set ES, i.e., the active-set ES
solves~\eqref{sec:expeval:eqn:bbobprobtransformedforactivesetes}.
$\mathbf{A}$ and $\mathbf{b}$ represent the BBOB COCO
constraint system of the current problem and $f$ is the
current problem's objective function.
As the DE variants, the CMA with augmented Lagrangian handling
and the SACOBRA are able to handle the form of the
BBOB COCO problem directly,
no problem transformation is necessary for them. Therefore, they
solve~\eqref{sec:expeval:eqn:bbobprob}.

\Cref{sec:expeval:compwithother:ecdfbbobcococomparisonall} shows ECDF graphs
for all the different algorithms. For every algorithm the ECDF aggregated over
all functions and dimensions in the BBOB COCO constrained suite is displayed.
Our approach (named \path{itprojlccmsaes} for the variant with
Iterative Projection and \path{l1lpsolvelccmsaes} for the variant
with the $\ell_1$ projection in the plots)
is among the best for all the dimensions. The CMA-ES
with augmented Lagrangian constraint handling
(named \path{cma_es_augmented_lagrangian} in the plots)
is able to reach about 50-60\% of the targets.
The Active-Set ES (named \path{activesetES} in the plots)
performs similarly to our approach for all the dimensions.
The DE variants, the Active-Set ES, and the SACOBRA
approach perform similarly as our approach in the smaller dimensions
but not as well in the larger dimensions. The exception to this
is the conSaDE that performs better for dimension $20$ and similarly
as our algorithm for dimension $40$. Similar to
\Cref{sec:expeval:lccmsaes:ecdf}, a closer look at the single
function ECDF plots for all the algorithms (not shown here)
reveals more insight. The lower performance in the higher dimensions
20 and 40 is mainly due to the Rastrigin problem for all the algorithms.
The CMA-ES with augmented Lagrangian constraint handling is only able
to solve a subset of the problems for all dimensions. In particular,
it is able to solve the Sphere, the Linear Slope, and the Different
Powers problems with $1$, $2$, and $6$ constraints. The SACOBRA
algorithm has problems with the higher number of constraints as well.

\Cref{sec:expeval:compwithother:ecdfkleemintycomparisonall} shows
ECDF graphs for all the different algorithms for the Klee-Minty problem with
dimensions $9$, $12$, and $15$ (plots for all the dimensions are presented
in the supplementary material (\Cref{sec:expeval:additionalresults:kleeminty})).
The SACOBRA approach, the $\varepsilon$DEag, and the
CMA with augmented Lagrangian constraint handling are only able to reach
20\% of the targets. All the other approaches perform well with some
difficulties in higher dimensions (13-15). In the large dimensions it is
worth noting that the numbers get large quickly and therefore numerical
stability can become an issue.

\providecommand{\bbobdatapath}{}
\renewcommand{\bbobdatapath}{figures/bbob_constrained_coco_2/linear_constraints/comparison_lccmsaes_and_others/ppdata/}
\subimport*{\bbobdatapath}{cocopp_commands.tex}
\providecommand{\pprldmanyalldatapath}{}
\renewcommand{\pprldmanyalldatapath}{\bbobdatapath activ_cma_e_conSa_ECHT-_epsDE_itpro_l1lps_et_al/}
\providecommand{\DIM}{}
\renewcommand{\DIM}{\ensuremath{\mathrm{DIM}}}
\providecommand{\aRT}{}
\renewcommand{\aRT}{\ensuremath{\mathrm{aRT}}}
\providecommand{\FEvals}{}
\renewcommand{\FEvals}{\ensuremath{\mathrm{FEvals}}}
\providecommand{\nruns}{}
\renewcommand{\nruns}{\ensuremath{\mathrm{Nruns}}}
\providecommand{\Dfb}{}
\renewcommand{\Dfb}{\ensuremath{\Delta f_{\mathrm{best}}}}
\providecommand{\Df}{}
\renewcommand{\Df}{\ensuremath{\Delta f}}
\providecommand{\nbFEs}{}
\renewcommand{\nbFEs}{\ensuremath{\mathrm{\#FEs}}}
\providecommand{\fopt}{}
\renewcommand{\fopt}{\ensuremath{f_\mathrm{opt}}}
\providecommand{\ftarget}{}
\renewcommand{\ftarget}{\ensuremath{f_\mathrm{t}}}
\providecommand{\CrE}{}
\renewcommand{\CrE}{\ensuremath{\mathrm{CrE}}}
\providecommand{\change}[1]{}
\renewcommand{\change}[1]{{\color{red} #1}}

\begin{figure*}
  \centering
  \begin{tabular}{@{\hspace*{-0.000\textwidth}}c@{\hspace*{-0.000\textwidth}}c@{\hspace*{-0.000\textwidth}}c}
    $D=2$ & $D=3$ & $D=5$\\
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_02D_noiselessall}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_03D_noiselessall}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_05D_noiselessall}\\
    $D=10$ & $D=20$ & $D=40$\\
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_10D_noiselessall}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_20D_noiselessall}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_40D_noiselessall}
  \end{tabular}
  \caption
      {Bootstrapped empirical cumulative distribution function
          of the number
        of objective function and constraint evaluations
        divided by dimension: comparison of all the approaches.
        \sppabbobecdfaxisexplanation}
  \label{sec:expeval:compwithother:ecdfbbobcococomparisonall}
\end{figure*}

\providecommand{\bbobdatapath}{}
\renewcommand{\bbobdatapath}{figures/kleeminty_coco/comparison/ppdata/}
\subimport*{\bbobdatapath}{cocopp_commands.tex}
\providecommand{\pprldmanyalldatapath}{}
\renewcommand{\pprldmanyalldatapath}{\bbobdatapath pprldmany-single-functions/}
\providecommand{\DIM}{}
\renewcommand{\DIM}{\ensuremath{\mathrm{DIM}}}
\providecommand{\aRT}{}
\renewcommand{\aRT}{\ensuremath{\mathrm{aRT}}}
\providecommand{\FEvals}{}
\renewcommand{\FEvals}{\ensuremath{\mathrm{FEvals}}}
\providecommand{\nruns}{}
\renewcommand{\nruns}{\ensuremath{\mathrm{Nruns}}}
\providecommand{\Dfb}{}
\renewcommand{\Dfb}{\ensuremath{\Delta f_{\mathrm{best}}}}
\providecommand{\Df}{}
\renewcommand{\Df}{\ensuremath{\Delta f}}
\providecommand{\nbFEs}{}
\renewcommand{\nbFEs}{\ensuremath{\mathrm{\#FEs}}}
\providecommand{\fopt}{}
\renewcommand{\fopt}{\ensuremath{f_\mathrm{opt}}}
\providecommand{\ftarget}{}
\renewcommand{\ftarget}{\ensuremath{f_\mathrm{t}}}
\providecommand{\CrE}{}
\renewcommand{\CrE}{\ensuremath{\mathrm{CrE}}}
\providecommand{\change}[1]{}
\renewcommand{\change}[1]{{\color{red} #1}}

\begin{figure*}
  \centering
  \begin{tabular}{@{\hspace*{-0.000\textwidth}}c@{\hspace*{-0.000\textwidth}}c@{\hspace*{-0.000\textwidth}}c}
    $D=9$ & $D=12$ & $D=15$\\
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_f001_09D}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_f001_12D}&
    \includegraphics[width=0.26\textwidth]{\pprldmanyalldatapath pprldmany_f001_15D}
  \end{tabular}
  \caption
      {Bootstrapped empirical cumulative distribution function
          of the number
        of objective function and constraint evaluations
        divided by dimension for the Klee-Minty problem with
        dimensions $D = 9$, $D = 12$ and $D = 15$:
        comparison of all the approaches.
        \sppabbobecdfaxisexplanation}
  \label{sec:expeval:compwithother:ecdfkleemintycomparisonall}
\end{figure*}

\section{Conclusion}
\label{sec:conclusion}
We have proposed the lcCMSA-ES; a CMSA-ES for solving optimization
problems with linear constraints. The algorithm is based on the CMSA-ES.
It is an interior point approach that repairs infeasible candidate solutions
if necessary. The mutation operator and repair
method are specially designed. They allow the ES to evolve itself on a
linear manifold. This distinguishes the proposed algorithm from the other
methods considered for the comparison.
It has been experimentally shown that the method works
well on the Klee-Minty optimization problem as well as the bbob-constrained
suite of the BBOB COCO framework (with disabled non-linear perturbations).
Additionally, the proposed lcCMSA-ES has been compared to other
evolutionary approaches for constrained optimization.
Experiments have shown that the lcCMSA-ES is among the best
for the BBOB COCO constrained suite and the Klee-Minty problem.
It is worth noting that not all algorithms that have been compared
are interior point methods. Consequently, they have the advantage
of evaluating the objective function outside the feasible region.
In particular, they do not move on the linear manifold defined
by the constraints.
All the three DE variants
considered (conSaDE, ECHT-DE, and $\varepsilon$DEag) allow infeasible
candidates. The ES with augmented Lagrangian constrained handling
works with a penalty. Thus, infeasible candidates are involved
during the evolution as well. The surrogate modeling with adaptive
parameter control does most of the optimization work on the surrogate
models. But already the computation of the initial surrogate models
involves evaluating objective and constraint functions.
For the initialization this is done at random points in the search space.
Those random points are not necessarily feasible.
The optimization on the surrogate models involves infeasible solutions
with respect to the real functions. But the results of the optimization
on the surrogate models are tried to be repaired if necessary with
respect to the real functions.
The Active-Set ES only considers feasible candidates. It starts
with an initial feasible candidate solution
and uses repair by projection for dealing with infeasible solutions.
This repair approach is not designed to move on the linear space
defined by the constraints. This is in contrast to our proposed ES.
Consequently, the objective function is only evaluated for feasible
candidate solutions.

\bibliographystyle{IEEEtran}


TITLE: How to go about solving $((1+iz)/(1-iz))^4 = 1/2 + i\sqrt3/2$?
QUESTION [4 upvotes]: I have problem solving this equation:
$$
 \left(\frac{1+iz}{1-iz}\right)^4 = \frac12 + i {\sqrt{3}\over 2}
$$
I know how to solve equations that are like:
$$
 w^4 = \frac12 + i {\sqrt{3}\over 2}
$$
And I have solved it to:
$$
    w = \cos(-\frac{\pi}{12} + \frac{\pi k}{2})) + i\sin(-\frac{\pi}{12} + \frac{\pi k}{2}))
$$
But now is:
$$
    w = \frac{1+iz}{1-iz}
$$
How does one get the complex z? Or am I solving it wrong?

REPLY [0 votes]: $$
w=\frac{1+iz}{1-iz}
$$
First, multiply both sides by $1-iz$:
$$
w(1-iz) = 1+iz
$$
Expand the left side:
$$
w-wiz = 1+iz
$$
Put all terms involving $z$ on one side and those not involving $z$ on the other side:
$$
w-1=iz+wiz
$$
Factor
$$
w-1 = iz(1+w)
$$
Divide both sides by $i(1+w)$:
$$
z= \frac{w-1}{i(w+1)}
$$
Multiply the numerator and denominator by the conjugate, $-i$:
$$
z = -i\frac{w-1}{w+1} = i\frac{1-w}{1+w}.
$$


TITLE: How to prove the equality of the floors of two expressions?
QUESTION [2 upvotes]: I was doing this problem from Olimpiada de Mayo (Argentina)
problem
when I wonder if
$$
\bigg\lfloor\frac{-1+\sqrt{1+8n}}{2}\bigg\rfloor = \bigg\lfloor\frac{-1+\sqrt{9+8n}}{2}\bigg\rfloor
$$
if $n\geq 3$ and $9+8n$ isn't a perfect square.

REPLY [1 votes]: There is no odd square strictly between $1+8n$ and $9+8n$ in any case. When there is also no even square between the two numbers, (and the larger one not a square) the two square roots have the same integer part.
When $$ 1+8n < 4 m^2 < 9 + 8n   $$
we find
$$  1+8n = 2m - \delta \; , \; \; \;  9+8n = 2m + \epsilon \; ,$$
with
$$ 0 < \delta , \epsilon < 1 $$
This is not as bad as it looks, because of the $-1$ additions. $$  \left\lfloor \frac{2m-1 \pm \delta}{2}\right\rfloor =  \left\lfloor \frac{2m-2}{2}\right\rfloor = m-1 $$
Example, with $n=4$  we get $33 < 36 < 41,$ $\sqrt{33} \approx 5.744 \approx 6 - 0.255 \; , \; \; \; $ $\sqrt{41} \approx 6.403 \approx 6 + 0.403 \; . \; \; \; $ $$ $$
Then $\sqrt{33} -1\approx 4.744 \approx 5 - 0.255 \; , \; \; \; $ $\sqrt{41}-1 \approx 5.403 \approx 5 + 0.403 \; . \; \; \; $ $$ $$
Then $ \frac{\sqrt{33} -1}{2} \approx \frac{5 - 0.255}{2} \approx 2.5 - 0.127 \; , \; \; \; $ $\frac{\sqrt{41}-1}{2}  \approx \frac{ 5 + 0.403}{2} \approx 2.5 + 0.2015 \; . \; \; \; $
and both have $2$ as integer part.


TITLE: For the minimum can we find a value of $n>0$ where $f(n)=0$?
QUESTION [0 upvotes]: Let $r_1, r_2, . . . , r_m$ be a given set of positive rational numbers whose sum is $1$. Define the function $f$ by $f(n) = n −\sum_{k=1}^{m} \lfloor r_kn \rfloor$ for each
positive integer $n$. Determine the minimum and maximum values of
$f(n)$.
Solution: Of course $\lfloor r_kn \rfloor$ ≤ $r_kn$, so $f(n) ≥ 0$, with equality for
$n = 0$, so $0$ is the minimum value. On the other hand, we have
$r_kn − \lfloor r_kn \rfloor$ < 1, so $f(n) ≤ m − 1$. Here equality holds for $n = t − 1$
if t is the least common denominator of the $r_k$.
For the minimum can we find a positive integer $n>0$ such that for fixed rational numbers $q_is$ satisfying $q_1 +q_2 +..=1$ we obtain $f(n)=0$ ?

REPLY [0 votes]: We have $f(n)=0$ for any $n$ such that $nr_k$ is an integer for each $k$. Clearly infinitely many such $n$ exist (if $r_k=\frac{p_k}{q_k}$, then take $n=q_1q_2\ldots q_m$ or any multiple thereof). In fact, it is for the maximum value that you have a problem: for any $x$, $x-\lfloor x\rfloor<1$ is a strict inequality, so we have always have $f(n)<m+1$.


TITLE: Sum of a rearranged alternating harmonic series, with three positive terms followed by one negative term
QUESTION [10 upvotes]: The series is this:
$$ 1 + 1/3 + 1/5 - 1/2 + 1/7 + 1/9 + 1/11 - 1/4 + 1/13 + 1/15 + 1/17 - 1/6 ...$$
The hint is to consider partial sums to $4n$ terms. I did that, and got the summation
$$ \sum 1/(6x-5) + 1/(6x-3) + 1/(6x-1) - 1/(2x) $$
But then I got stuck. The terms don't seem to simplify to anything... Wolfram alpha gives the result of this summation as $\log(12)/2$, but I have no idea how to get there.
Any help greatly appreciated.

REPLY [6 votes]: Hint:
Try writing the sum as
$$
\sum_{k=0}^\infty \frac{1}{6k+1}+\frac{1}{6k+3}+\frac{1}{6k+5}-\frac{3}{6k+6}\tag{1}
$$
Then use the fact that
$$
\begin{align}
\sum_{k=0}^n\frac{1}{2k+1}
&=\log(2n+1)+\gamma-\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)\\
&=\log(2)+\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)\tag{2}
\end{align}
$$
and that
$$
\begin{align}
\sum_{k=0}^n\frac{1}{2k+2}
&=\tfrac{1}{2}(\log(n+1)+\gamma)+O\left(\tfrac{1}{n}\right)\\
&=\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)\tag{3}
\end{align}
$$
Addendum:
Both formulas $(2)$ and $(3)$ follow from the asymptotic expansion for the partial harmonic series:
$$
\sum_{k=1}^n\frac{1}{k}=\log(n)+\gamma+O\left(\tfrac{1}{n}\right)
$$
where $\gamma$ is the Euler-Mascheroni Constant. Thus,
$$
\begin{align}
\sum_{k=0}^n\frac{1}{2k+1}
&=\sum_{k=1}^{2n+1}\frac{1}{k}-\sum_{k=1}^{n}\frac{1}{2k}\\
&=\log(2n+1)+\gamma-\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)\\
&=\log(2)+\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)
\end{align}
$$
and
$$
\begin{align}
\sum_{k=0}^n\frac{1}{2k+2}
&=\sum_{k=1}^{n+1}\frac{1}{2k}\\
&=\tfrac{1}{2}(\log(n+1)+\gamma)+O\left(\tfrac{1}{n}\right)\\
&=\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)
\end{align}
$$
Using $(2)$, we get
$$
\begin{align}
\sum_{k=0}^n \frac{1}{6k+1}+\frac{1}{6k+3}+\frac{1}{6k+5}
&=\log(2)+\tfrac{1}{2}(\log(3n+2)+\gamma)+O\left(\tfrac{1}{n}\right)\\
&=\tfrac{1}{2}(\log(n)+\log(12)+\gamma)+O\left(\tfrac{1}{n}\right)\tag{4}
\end{align}
$$
and using $(3)$, we get
$$
\begin{align}
\sum_{k=0}^n \frac{3}{6k+6}
&=\tfrac{1}{2}(\log(n)+\gamma)+O\left(\tfrac{1}{n}\right)\tag{5}
\end{align}
$$
Combining $(4)$ and $(5)$, we get
$$
\sum_{k=0}^n \frac{1}{6k+1}+\frac{1}{6k+3}+\frac{1}{6k+5}-\frac{3}{6k+6}=\tfrac{1}{2}\log(12)+O\left(\tfrac{1}{n}\right)
$$


TITLE: Powers of a Positive Matrix in the Limit
QUESTION [8 upvotes]: I'm trying to prove a standard result: for a positive $n \times n$ matrix $A$, the powers of $A$ scaled by its leading eigenvalue $\lambda$ converge to a matrix whose columns are just scalar multiples of $A$'s leading eigenvector $\mathbf{v}$. More precisely,
$$ \lim_{k \rightarrow \infty} \left(\frac{A}{\lambda}\right)^k = \mathbf{v}\mathbf{u},$$
where $\mathbf{v}$ and $\mathbf{u}$ are the leading right and left eigenvectors of $A$, respectively (scaled so that $\mathbf{u}\mathbf{v} = 1$).
These notes give a nice compact proof, pictured below. (It covers the more general case where $A$ is primitive, but I'm happy to assume it's positive for my purposes.)
I'm stuck on the highlighted step. Having established the existence of a matrix $M$ that (a) fixes $\mathbf{u}$ and $\mathbf{v}$, and (b) annihilates all other generalized eigenvalues of $A$, how do we know $M$ is unique? Why couldn't there be other matrices satisfying (a) and (b)?
I don't know much about generalized eigenvectors. I gather they're linearly independent, hence form a basis for $\mathbb{R}^n$. So each column of $M$ must be a unique linear combination of $A$'s generalized right eigenvectors. Is there some path I'm not seeing from there to the conclusion that only one matrix can satisfy both (a) and (b)?

REPLY [1 votes]: You may dispense with left eigenvectors when proving the result,
but calling for the one of them is a comfortable path to the explicit expression for
$$M=\lim_{k\to\infty}J^k$$
where $\,J=\lambda^{-1}\!A\,\in M_n\big(\mathbb R^{>0}\big)\,$ for notational convenience.
All the good properties from $A$, due to Perron & Frobenius, carry over to $J$:
Its leading eigenvalue is $1$, the eigenspace is spanned by one positive vector $v$, and all its other eigenvalues live in the interior of the complex unit disc.
This is the reason why $M\,$ has an $n\!-\!1$-dimensional kernel, as worked out in M. Boyle's notes, cf the question-embedded screen-shot, and only the fixed points of $J$ are not annihilated in the limit, i. e., the $1$-dimensional eigenspace to the eigenvalue $1$ which equals $\operatorname{span}(v)$.
Hence $M\,$ has rank $1$, and it is uniquely determined because $M$ is known on all of $\,\mathbb R^n$.
Notice that left eigenvectors of a matrix are in 1-to-1 correspondence to the right eigenvectors of its transpose. The transpose $J^T$ satisfies the assumptions for Perron & Frobenius as well as $J$, whence there exists a positive left eigenvector $u$ of $J$ to the eigenvalue $1$. Analogous to $v$, it survives in the limit.
Since $M\,$ has rank $1$, let's make the ansatz $M=v\,\langle ?|\cdot\rangle\,$. Then $M^T\!=\, ?\,\langle v|\cdot\rangle\,$, and
$$Mv =v\,\langle ?|v\rangle = v\quad\text{and}\quad
M^T\!u=\, ?\,\langle v|u\rangle = u$$
yield $M=v\,\langle u|\cdot\rangle \equiv vu\,$ with $\,uv\equiv\langle u|v\rangle =1$.
Note that $M$ is idempotent. It is selfadjoint (hence an orthogonal projector) if and only if $v$ and $u$ are scalar multiples of each other (thus linear dependent).
Furthermore, one has $v=\lim\limits_{k\to\infty}\,(J^kw)\,$ where $w>0\,$ is an arbitrary positive vector.
PS & IMHO
$(1)\:$ If you seek to master this result and its proof, you shouldn't consider
projectors, kernels, images as "additional machinery" but integrate them into your picture.
$(2)\:$ In this same spirit I'd like to add that loup blanc's answer follows M. Boyle's path more than it's different to it.


TITLE: $\sum_1^{\infty} \frac{(p+1)(p+2)(p+3)...(p+n)}{(q+1)(q+2)(q+3)...(q+n)}$ convergence
QUESTION [7 upvotes]: I need to determine for which values of $p$ and $q$, both greater than $0$, the following series converges: $$\sum_1^{\infty} \frac{(p+1)(p+2)(p+3)...(p+n)}{(q+1)(q+2)(q+3)...(q+n)}$$
I've tried using the ratio test, comparison test, and I've also tried partial fraction decomposition but I can't get to anything.
Could you give me a hint on how to solve this?
Any help will be appreciated

REPLY [5 votes]: As you have probably noticed, the ratio test is useless here. The next thing to try when the ratio test fails is the Raabe-Duhame test. In our case this amounts to computing the following limit: $\lim \limits _{n \to \infty} n(\frac {q+n+1} {p+n+1} -1) = \lim \limits _{n \to \infty} \frac {(q-p)n} {p+n+1} = q-p$. Now, if $q-p < 1$ the series diverges and if $q-p>1$ the series converges. If $q-p=1$ then the Raabe-Duhamel test fails, but in this case the series becomes $\sum \frac {p+1} {p+1+n}$. Since $p \geq 0$, the series has the same behaviour as the series $\sum \frac 1 n$ which is divergent.
Therefore the series converges only if $q-p>1$.