BCMA-ES: A Bayesian approach to CMA-ES

TL;DR

This paper presents a Bayesian framework for CMA-ES, deriving optimal parameter updates using conjugate priors, leading to two new variants with demonstrated fast convergence in experiments.

Contribution

It introduces a theoretically grounded Bayesian approach to CMA-ES, resulting in two novel algorithm versions based on different prior assumptions.

Findings

Fast convergence of the proposed CMA-ES variants

Theoretical justification for CMA-ES updates

Introduction of normal-Wishart and normal-Inverse Wishart based CMA-ES variants

Abstract

This paper introduces a novel theoretically sound approach for the celebrated CMA-ES algorithm. Assuming the parameters of the multi variate normal distribution for the minimum follow a conjugate prior distribution, we derive their optimal update at each iteration step. Not only provides this Bayesian framework a justification for the update of the CMA-ES algorithm but it also gives two new versions of CMA-ES either assuming normal-Wishart or normal-Inverse Wishart priors, depending whether we parametrize the likelihood by its covariance or precision matrix. We support our theoretical findings by numerical experiments that show fast convergence of these modified versions of CMA-ES.

Equations

$$p(x_1,x_2,...,x_n)=p(x_{\pi 1},x_{\pi 2},...,x_{\pi n})$$

$$p(x_1,x_2,...,x_n)=\int \prod _{i=1}^{n}p(x_i|\theta )P(d\theta ),$$

$$\text{Posterior}\propto \text{Prior}\times \text{Likelihood}$$

$$\varphi (\theta ,x)=p(x|\theta )\pi (\theta )$$

$$m(x)=\int \varphi (\theta ,x)d\theta =\int p(x|\theta )\pi (\theta )d\theta$$

$$\pi (\theta |x)=\frac{p(x|\theta )\pi (\theta )}{\int p(x|\theta )\pi (\theta )d\theta }\propto p(x|\theta )\pi (\theta )$$

$$\pi (\theta |x)\propto p(x|\theta )\pi (\theta )$$

$$p(\mathbf{x}|\eta )=h(\mathbf{x})\exp (\eta \cdot T(\mathbf{x})-A(\eta )),$$

$$A(\eta )\triangleq \log {\int }_{\mathcal{X}}h(\mathbf{x})\exp (\eta \cdot T(\mathbf{x}))\mathrm{d}x.$$

$$\eta (\mu ,\Sigma )$$

$$h(\mathbf{x})$$

$$f(\boldsymbol{\mu},\boldsymbol{\Sigma}|\boldsymbol{\mu}_{0},\lambda,\boldsymbol{\Psi},\nu)=\mathcal{N}\left(\boldsymbol{\mu}\Big{|}\boldsymbol{\mu}_{0},\frac{1}{\lambda}\boldsymbol{\Sigma}\right)\mathcal{W}^{-1}(\boldsymbol{\Sigma}|\boldsymbol{\Psi},\nu)$$

$$\begin{array}{rl}{\mu }_0^{\star }& =\frac{{\kappa }_0{\mu }_0+n\overline{x}}{{\kappa }_0+n},\\ {\kappa }_0^{\star }& ={\kappa }_0+n,\\ v_0^{\star }& =v_0+n\\ {\psi }^{\star }& =\psi +\sum _{i=1}^n\left(x_i-\overline{x}\right){\left(x_i-\overline{x}\right)}^T+\frac{{\kappa }_{0}n}{{\kappa }_0+n}\left(\overline{x}-{\mu }_0\right){\left(\overline{x}-{\mu }_0\right)}^T\end{array}$$

$${\mathbb{E}}_t[\mu ]={\mu }_t$$

$${\mathbb{E}}_t[\Sigma ]=\frac{{\psi }_t}{v_t-n-1}$$

$${\overline{\mu }}_t=\underset{\mathrm{MC}\mathrm{mean}\mathrm{for}{X}_{f\uparrow }}{\underbrace{\sum _{i=1}^k{w}_{(i),w\downarrow }\cdot X_{(i),f\uparrow }}}-\underset{\mathrm{MC}\mathrm{bias}\mathrm{for}X}{\underbrace{\left(\sum _{i=1}^k{w}_i{X}_i-{\overline{\mu }}_t\right)}}$$

$${\overline{\mu }}_t=\underset{X\in \mathcal{X}}{\mathrm{argmin}}f(X)$$

$${\overline{\Sigma }}_t=$$

$$\begin{array}{rl}{\mu }_{t+1}& =\frac{{\kappa }_t{\mu }_t+n{\overline{\mu }}_t}{{\kappa }_t+n},\\ {\kappa }_{t+1}& ={\kappa }_t+n,\\ v_{t+1}& =v_t+n,\\ {\psi }_{t+1}& ={\psi }_t+{\overline{\Sigma }}_t+\frac{{\kappa }_{t}n}{{\kappa }_t+n}\left({\overline{\mu }}_t-{\mu }_t\right){\left({\overline{\mu }}_t-{\mu }_t\right)}^T\end{array}$$

$$\begin{array}{rl}{\mathbb{E}}_{t+1}[\mu ]& ={\mathbb{E}}_t[\mu ]+w_t^{\mu }\left({\overline{\mu }}_t-{\mathbb{E}}_t[\mu ]\right),\\ {\mathbb{E}}_{t+1}[\Sigma ]& =\underset{\text{discount factor}}{\underbrace{w_t^{\Sigma ,1}}}{\mathbb{E}}_t[\Sigma ]+w_t^{\Sigma ,2}\underset{\text{rank one matrix}}{\underbrace{\left({\overline{\mu }}_t-{\mathbb{E}}_t[\mu ]\right){\left({\overline{\mu }}_t-{\mathbb{E}}_t[\mu ]\right)}^T}}\\ & +w_t^{\Sigma ,3}\underset{\text{rank (n-1) matrix}}{\underbrace{{\overline{\Sigma }}_t}}\\ \mathrm{where}{w}_t^{\mu }& =\frac{n}{{\kappa }_t+n},\\ w_t^{\Sigma ,1}& =\frac{{\kappa }_{t}n}{({\kappa }_t+n)(v_t-1)}\\ w_t^{\Sigma ,2}& =\frac{v_t-n-1}{v_t-1},\\ w_t^{\Sigma ,3}& =\frac{1}{v_t-1}\end{array}$$

$$\begin{array}{rl}{m}_{t+1}& =m_t+\sum _{i=1}^{\mu }{w}_i(x_{i:\lambda }-m_t)\\ C_{t+1}& =\underset{\text{discount factor}}{\underbrace{(1-c_1-c_{\mu }+c_s)}}{C}_t+c_1\underset{\text{rank one matrix}}{\underbrace{p_c{p}_c^T}}\\ & +c_{\mu }\underset{\text{rank}\min (\mu ,n-1)\text{matrix}}{\underbrace{\sum _{i=1}^{\mu }{w}_i\frac{x_{i:\lambda }-m_k}{{\sigma }_k}{\left(\frac{x_{i:\lambda }-m_t}{{\sigma }_t}\right)}^T}}\end{array}$$

$$f(x)={\left(\sum _{i=1}^n{x}_i^2\right)}^{1/2}=\Vert x{\Vert }_2$$

$$f(x)=\sum _{i=1}^n\mid x_i\mid +\prod _{i=1}^n\mid x_i\mid$$

$$f(x)=10\times n+\sum _{i=1}^n\left[x_i^2-10\cos (2\pi x_i)\right]$$

$$f(x)=418.9829\times n$$

$$-\sum _{i=1}^n\left[x_i\sin (\sqrt{\mid x_i\mid })\text{mathbbm}{1}_{|x_i|<500}+500\sin (\sqrt{500})\text{mathbbm}{1}_{|x_i|\ge 500}\right]$$

$$\displaystyle p(\mathcal{X}|\theta,\lambda)\!=\!\Big{(}\prod_{j=1}^{n}h({\mathbf{x}}^{j})\Big{)}\exp\Big{(}\eta(\theta,\lambda)^{T}\sum_{j=1}^{n}T(x^{j})-nA(\eta(\theta,\lambda))\Big{)}.$$

$$\pi (\theta |\mathcal{X})$$

$$\propto \exp \left(\eta (\theta ,\lambda )\cdot \sum _{j=1}^{n}T(x^j)-nA(\eta (\theta ,\lambda ))+\mathcal{F}(\theta )\right).$$

$$\mathcal{F}(\theta )={\mu }_1\cdot \eta (\theta ,\lambda )-{\mu }_{0}A(\eta (\theta ,\lambda ))$$

$$\displaystyle\!\!\!\!p(\mathcal{X}|\theta,\lambda)\!\propto\!\exp\Big{(}\Big{(}\mu_{1}+\sum_{j=1}^{n}T(x^{j})\Big{)}^{T}\!\!\eta(\theta,\lambda)-(n+\mu_{0})A(\eta(\theta,\lambda))\Big{)}\!.$$

$$p(\mathcal{X}|\theta ,\lambda )=\frac{1}{Z}\exp \left({\mu }_1\cdot \eta (\theta ,\lambda )-{\mu }_{0}A(\eta (\theta ,\lambda ))\right),$$

Full text

CMA-ES, Bayesian, conjugate prior, normal-inverse-Wishart

††copyright: rightsretained††doi: 10.1145/nnnnnnn.nnnnnnn††isbn: 978-x-xxxx-xxxx-x/YY/MM††conference: A.I Square Working Paper; March 2019; France††journalyear: 2019††price: 15.00††ccs: Mathematics of computing Probability and statistics

1. Introduction

The covariance matrix adaptation evolution strategy (CMA-ES)  (Hansen and Ostermeier, 2001) is arguably one of the most powerful real-valued derivative-free optimization algorithms, finding many applications in machine learning. It is a state-of-the-art optimizer for continuous black-box functions as shown by the various benchmarks of the COmparing Continuous Optimisers INRIA platform for ill-posed functions. It has led to a large number of papers and articles and we refer the interested reader to (Hansen and Ostermeier, 2001; Auger et al., 2004; Igel et al., 2007; Auger and Hansen, 2009; Hansen and Auger, 2011; Auger and Hansen, 2012; Hansen and Auger, 2014; Akimoto et al., 2015, 2016; Ollivier et al., 2017) and (Varelas et al., 2018) to cite a few.

It has has been successfully applied in many unbiased performance comparisons and numerous real-world applications. In particular, in machine learning, it has been used for direct policy search in reinforcement learning and hyper-parameter tuning in supervised learning ((Igel et al., 2009; Heidrich-Meisner and Igel, 2009; Igel, 2010)), and references therein, as well as hyperparameter optimization of deep neural networks  (Loshchilov and Hutter, 2016)

In a nutshell, the ($\mu$ / $\lambda$) CMA-ES is an iterative black box optimization algorithm, that, in each of its iterations, samples $\lambda$ candidate solutions from a multivariate normal distribution, evaluates these solutions (sequentially or in parallel) retains $\mu$ candidates and adjusts the sampling distribution used for the next iteration to give higher probability to good samples. Each iteration can be individually seen as taking an initial guess or prior for the multi variate parameters, namely the mean and the covariance, and after making an experiment by evaluating these sample points with the fit function updating the initial parameters accordingly.

Historically, the CMA-ES has been developed heuristically, mainly by conducting experimental research and validating intuitions empirically. Research was done without much focus on theoretical foundations because of the apparent complexity of this algorithm. It was only recently that (Akimoto et al., 2010; Glasmachers et al., 2010) and (Ollivier et al., 2017) made a breakthrough and provided a theoretical justification of CMA-ES updates thanks to information geometry. They proved that CMA-ES was performing a natural gradient descent in the Fisher information metric. These works provided nice explanation for the reasons of the performance of the CMA-ES because of strong invariance properties, good search directions, etc

There is however another way of explanation that has been so far ignored and could also bring nice insights about CMA-ES. It is Bayesian statistics theory. At the light of Bayesian statistics, CMA-ES can be seen as an iterative prior posterior update. But there is some real complexity due to tricky updates that may explain why this has always been ignored. First of all, in a regular Bayesian approach, all sample points are taken. This is not the case in the ($\mu$/$\lambda$) CMA-ES as out of the $\lambda$ generated paths, only the $\mu$ best are selected. The updating weights are also constant which is not consistent with Bayesian updates. But more importantly, the covariance matrix update is the core of the problem. It appeals important remarks. The update is done according to a weighted combination of a rank one matrix referred to $p_{C}p_{C}^{T}$ with parameter $c_{1}$ and a rank $min(\mu,n)$ matrix with parameter $c_{\mu}$, whose details are given for instance in (Hansen, 2016). The two updates for the covariance matrix makes the Bayesian update interpretation challenging as these updates are done according to two paths: the isotropic and anisotropic evolution path. All this may explain why a Bayesian approach for interpreting and revisiting the CMA-ES algorithm have seemed a daunting task and not tackled before.

This is precisely the objective of this paper. Section 2 recalls various Bayesian concepts of updates for prior and posterior to highlight the analogy of an iterative Bayesian update. Section 3 presents in greater details the Bayesian approach of CMA-ES, with the corresponding family of derived algorithms, emphasizing the various design choices that can conduct to multiple algorithms. Section 4 provides numerical experiments and shows that Bayesian adapted CMA-ES algorithms perform well on convex and non convex functions. We finally conclude about some possible extensions and further experiments.

However, the analogy with a successive Bayesian prior posterior update has been so far missing in the landscape of CMA-ES for multiple reasons. First of all, from a cultural point of view, the evolutionary and Bayesian community have always been quite different and not overlapping. Secondly, the CMA-ES was never formulated in terms of a prior and posterior update making its connection with Bayesian world non obvious. Thirdly, when looking in details at the parameters updates, the weighted combination between the global and local search makes the interpretation of a Bayesian posterior update non trivial. We will explain in this paper that the global search needs to be addressed with a special dilatation techniques that is not common in Bayesian wold.

2. Framework

CMA-ES computes at each step an update of the mean and covariance of the distribution of the minimum. From a very general point of view this can be interpreted as a prior posterior update in Bayesian statistics.

2.1. Bayesian vs Frequentist probability theory

The justification of the Bayesian approach is discussed in (Robert, 2007). In Bayesian probability theory, we assume a distribution on unknown parameters of a statistical model that can be characterized as a probabilization of uncertainty. This procedure leads to an axiomatic reduction from the notion of unknown to the notion of randomness but with probability. We do not know the value of the parameters for sure but we know specific values that these parameters can take with higher probabilities. This creates a prior distribution that is updated as we make some experiments as shown in  (Gelman et al., 2004; Marin and Robert, 2007; Robert, 2007). In the Bayesian view, a probability is assigned to a hypothesis, whereas under frequentist inference, a hypothesis is typically tested without being assigned a probability. There are even some nice theoretical justification for it as presented in (Jordan, 2010).

Definition 2.0.

(Infinite exchangeability). We say that $(x_{1},x_{2},...)$ is an infinitely exchangeable sequence of random variables if, for any n, the joint probability $p(x_{1},x_{2},...,x_{n})$ is invariant to permutation of the indices. That is, for any permutation $\pi$,

[TABLE]

Equipped with this definition, the De Finetti’s theorem as provided below states that exchangeable observations are conditionally independent relative to some latent variable.

Theorem 2.2.

(De Finetti, 1930s). A sequence of random variables $(x_{1},x_{2},...)$ is infinitely exchangeable iff, for all n,

[TABLE]

for some measure P on $\theta$.

This representation theorem 2.2 justifies the use of priors on parameters since for exchangeable data, there must exist a parameter $\theta$, a likelihood $p(x|\theta)$ and a distribution $\pi$ on $\theta$. A proof of De Finetti theorem is for instance given in (Schervish, 1996) (section 1.5).

Remark 2.1.

The De Finetti is trivially satisfied in case of i.i.d. sampling as the sequence is clearly exchangeable and that the joint probability is clearly given by the product of all the marginal distributions. However, the De Finetti goes far beyond as it proves that the infinite exchangeability is enough to prove that the joint distribution is the product of some marginal distribution for a given parameter $\theta$. The sequence may not be independent neither identically distributed, which is a much stronger result!

2.2. Conjugate priors

In Bayesian statistical inference, the probability distribution that expresses one’s (subjective) beliefs about the distribution parameters before any evidence is taken into account is called the prior probability distribution, often simply called the prior. In CMA-ES, it is the distribution of the mean and covariance. We can then update our prior distribution with the data using Bayes’ theorem to obtain a posterior distribution. The posterior distribution is a probability distribution that represents your updated beliefs about the parameters after having seen the data. The Bayes’ theorem tells us the fundamental rule of Bayesian statistics, that is

[TABLE]

The proportional sign indicates that one should compute the distribution up to a renormalization constant that enforces the distribution sums to one. This rule is simply a direct consequence of Baye’s theorem. Mathematically, let us say that for a random variable $X$, its distribution $p$ depends on a parameter $\theta$ that can be multi-dimensional. To emphasize the dependency of the distribution on the parameters, let us write this distribution as $p(x|\theta)$ and let us assume we have access to a prior distribution $\pi(\theta)$. Then the joint distribution of $(\theta,x)$ writes simply as

[TABLE]

The marginal distribution of $x$ is trivially given by marginalizing the joint distribution by $\theta$ as follows:

[TABLE]

The posterior of $\theta$ is obtained by Bayes’s formula as

[TABLE]

Computing a posterior is tricky and does not bring much value in general. A key concept in Bayesian statistics is conjugate priors that makes the computation really easy and is described at length below.

Definition 2.0.

A prior distribution $\pi(\theta)$ is said to be a conjugate prior if the posterior distribution

[TABLE]

remains in the same distribution family as the prior.

At this stage, it is relevant to introduce exponential family distributions as this higher level of abstraction that encompasses the multi variate normal trivially solves the issue of founding conjugate priors. This will be very helpful for inferring conjugate priors for the multi variate Gaussian used in CMA-ES.

Definition 2.0.

A distribution is said to belong to the exponential family if it can be written (in its canonical form) as:

[TABLE]

where $\eta$ is the natural parameter, $T({\mathbf{x}})$ is the sufficient statistic, $A(\eta)$ is log-partition function and $h({\mathbf{x}})$ is the base measure. $\eta$ and $T({\mathbf{x}})$ may be vector-valued. Here $a\cdot b$ denotes the inner product of $a$ and $b$.

The log-partition function is defined by the integral

[TABLE]

Also, $\eta\in\Omega=\{\eta\in\mathbb{R}^{m}|A(\theta)<+\infty\}$ where $\Omega$ is the natural parameter space. Moreover, $\Omega$ is a convex set and $A(\cdot)$ is a convex function on $\Omega$.

Remark 2.2.

Not surprisingly, the normal distribution $\mathcal{N}({\mathbf{x}};\mu,\Sigma)$ with mean $\mu\in\mathbb{R}^{d}$ and covariance matrix $\Sigma$ belongs to the exponential family but with a different parametrisation. Its exponential family form is given by:

[TABLE]

where in equations (4a), the notation $\mathrm{vec}(\cdot)$ means we have vectorized the matrix, stacking each column on top of each other and hence can equivalently write for $a$ and $b$, two matrices, the trace result $\operatorname{Tr}(a^{\mathrm{T}}b)$ as the scalar product of their vectorization $\mathrm{vec}(a)\cdot\mathrm{vec}(b)$. We can remark the canonical parameters are very different from traditional (also called moment) parameters. We can notice that changing slightly the sufficient statistic $T(x)$ leads to change the corresponding canonical parameters $\eta$.

For an exponential family distribution, it is particularly easy to form conjugate prior.

Proposition 2.5.

If the observations have a density of the exponential family form p(x|\theta,\lambda)=h(x)\exp\Big{(}\eta(\theta,\lambda)^{T}T(x)-nA(\eta(\theta,\lambda))\Big{)}, with $\lambda$ a set of hyper-parameters, then the prior with likelihood defined by $\pi(\theta)\propto\exp\left(\mu_{1}\cdot\eta(\theta,\lambda)-\mu_{0}A(\eta(\theta,\lambda))\right)$ with $\mu\triangleq(\mu_{0},\mu_{1})$ is a conjugate prior.

The proof is given in appendix subsection 6.1. As we can vary the parameterisation of the likelihood, we can obtain multiple conjugate priors. Because of the conjugacy, if the initial parameters of the multi variate Gaussian follows the prior, the posterior is the true distribution given the information $\mathcal{X}$ and stay in the same family making the update of the parameters really easy. Said differently, with conjugate prior, we make the optimal update. And it is enlightening to see that as we get some information about the likelihood, our posterior distribution becomes more peak as shown in figure1.

2.3. Optimal updates for NIW

The two natural conjugate priors for the Multi variate normal that updates both the mean and the covariance are the normal-inverse-Wishart if we want to update the mean and covariance of the Multi variate normal or the normal-Wishart if we are interested in updating the mean and the precision matrix (which is the inverse of the covariance matrix). In this paper, we will stick to the normal-inverse-Wishart to keep things simple. The Normal-inverse-Wishart distribution is parametrized by $\boldsymbol{\mu}_{0},\lambda,\boldsymbol{\Psi},\nu$ and its distribution is given by

[TABLE]

where $\mathcal{W}^{-1}$ denotes the inverse Wishart distribution. The key theoretical guarantee of the BCMA-ES is to update the mean and covariance of our CMA-ES optimally as follows.

Proposition 2.6.

If our sampling density follows a $d$ dimensional multivariate normal distribution $\sim\mathcal{N}_{d}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ with unknown mean $\mu$ and covariance $\Sigma$ and if its parameters are distributed according to a Normal-Inverse-Wishart $(\mu,\Sigma)\sim\mathrm{NIW}(\mu_{0},\kappa_{0},v_{0},\psi)$ and if we observe $\mathcal{X}=(x_{1},..,x_{n})$ samples, then the posterior is also a Normal-Inverse-Wishart with different parameters $\mathrm{NIW}(\mu_{0}^{\star},\kappa_{0}^{\star},v_{0}^{\star},\psi^{\star})$ given by

[TABLE]

with $\overline{x}$ the sample mean.

Remark 2.3.

This proposition is the cornerstone of the BCMA-ES. It provides the theoretical guarantee that the updates of the parameters in the algorithm are accurate and optimal under the assumption of the prior. In particular, this implies that any other formula for the update of the mean and variance and in particular the ones used in the mainstream CMA-ES assumes a different prior.

Proof.

A complete proof is given in the appendix section 6.2. ∎

3. Bayesian CMA-ES

3.1. Main assumptions

Our main assumptions are the followings :

• •

the parameters of the multi-variate Gaussian follow a conjugate prior distribution.

• •

the minimum of our objective function $f$ follows a multi-variate normal law.

3.2. Simulating the minimum

One of the main challenge is to simulate the likelihood to infer the posterior. The key question is really to use the additional information of the function value $f$ for candidate points. At step $t$ in our algorithm, we suppose multi variate Gaussian parameters $\mu$ and $\Sigma$ follow a normal inverse Wishart denoted by $NIW(\mu_{t},\kappa_{t},v_{t},\psi_{t})$.

In full generality, we need to do a Monte Carlo of Monte Carlo as the parameters of our multi variate normal are themselves stochastic. However, we can simplify the problem and take their mean values. It is very effective in terms of computation and reduces Monte Carlo noise. For the normal inverse Wishart distribution, there exist closed form for these mean values given by:

[TABLE]

and

[TABLE]

We simulate potential candidates \mathcal{X}=\{X_{i}\}\sim\mathcal{N}\big{(}\mathbb{E}_{t}[\mu],\mathbb{E}_{t}[\Sigma]\big{)} and evaluate them $f(X_{i})$. If the distribution of the minimum was accurate, the minimum would concentrate around $\mathbb{E}_{t}[\mu]$ and be spread with a variance of $\mathbb{E}_{t}[\Sigma]$. When evaluating potential candidates, as our guess is not right, we do not get values centered around $\mathbb{E}_{t}[\mu]$ and spread with a variance of $\mathbb{E}_{t}[\Sigma]$. This comes from three things:

• •

Our assumed minimum is not right. We need to shift our normal to the right minimum!

• •

Our assumed variance is not right. We need to compute it on real data taken into additional information given by $f$.

• •

Last but not least, our Monte Carlo simulation adds some random noise.

For the last issue, we can correct any of our estimator by the Monte Carlo bias. This can be done using standard control variate as the simulated mean and variance are given: $\mathbb{E}_{t}[\mu]$ and $\mathbb{E}_{t}[\Sigma]$ respectively and we can compute for each of them the bias explicitly.

The first two issues are more complex. Let us tackle each issue one by one.

To recover the true minimum, we design two strategies.

• •

We design a strategy where we rebuild our normal distribution but using sorted information of our $X$’s weighted by their normal density to ensure this is a true normal corrected from the Monte Carlo bias. We need to explicitly compute the weights. For each simulated point $X_{i}$, we compute it assumed density denoted by $d_{i}=\mathcal{N}(\mathbb{E}_{t}[\mu],\mathbb{E}_{t}[\Sigma])(X_{i})$ where $\mathcal{N}(\mathbb{E}_{t}[\mu],\mathbb{E}_{t}[\Sigma])(.)$ denotes the p.d.f. of the multi-variate Gaussian.

We divide these density by their sum to get weights $(w_{i})_{i=1..k}$ that are positive and sum to one as follows. $w_{j}=d_{j}/\sum_{i=1}^{k}d_{i}$. Hence for $k$ simulated points, we get $\{X_{i},w_{i}\}_{i=1..k}$. We reorder jointly the uplets (points and density) in terms of their weights in decreasing order.

To insist we take sorted value in decreasing order with respect to the weights $(w_{i})_{i=1..k}$, we denote the order statistics $(i),w\downarrow$.

This first sorting leads to k new uplets $\{X_{(i),w\downarrow},w_{(i),w\downarrow}\}_{i=1..k}$. Using a stable sort (that keeps the order of the density), we sort jointly the uplets (points and weights) according to their objective function value (in increasing order this time) and get a k new uplets $\{X_{(i),f\uparrow},w_{(i),w\downarrow}\}_{i=1..k}$. We can now compute the empirical mean $\overline{\mu}_{t}$ as follows:

[TABLE]

The intuition of equation (8) is to compute in the left term the Monte Carlo mean using reordered points according to their objective value and correct our initial computation by the Monte Carlo bias computed as the right term, equal to the initial Monte Carlo mean minus the real mean. We call this strategy one.

• •

If we think for a minute about the strategy one, we get the intuition that when starting the minimization, it may not be optimal. This is because weights are proportional to $\exp\left\{\frac{1}{2}(X-\mathbb{E}_{t}[\mu])^{T}(\mathbb{E}_{t}[\Sigma])^{-1}(X-\mathbb{E}_{t}[\mu])\right\}$.

When we start the algorithm, we use a large search space, hence a large covariance matrix $\overline{\Sigma}_{t}$ which leads to have weights which are quite similar. Hence even if we sort candidates by their fit, ranking them according to the value of $f$ in increasing order, we will move our theoretical multi variate Gaussian little by little. A better solution is more to brutally move the center of our multi variate Gaussian to the best candidate seen so far, as follows:

[TABLE]

We call this strategy two. Intuitively, strategy two should be best when starting the algorithm while strategy one would be better once we are close to the solution.

To recover the true variance, we can adapt what we did in stratey one as follows:

• •

[TABLE]

where $\overline{X}_{(.),f\uparrow}=\sum_{i=1}^{k}w_{(i),w\downarrow}X_{(i),f\uparrow}$ and $\overline{X}=\sum_{i=1}^{k}w_{i}X_{i}$ are respectively the mean of the sorted and non sorted points.

• •

Again, we could design another strategy that takes part of the points but we leave this to further research.

Once we have the likelihood mean and variance using (9) and (10) or (8) and (10), we update the posterior law according to equation (5). This gives us the iterative conjugate prior parameters updates:

[TABLE]

The resulting algorithm is summarized in Algo 1.

Proposition 3.1.

Under the assumption of a NIW prior, the updates of the BCMA-ES parameters for the expected mean and variance write as a weighted combination of the prior expected mean and variance and the empirical mean and variance as follows

[TABLE]

Remark 3.1.

The proposition above is quite fundamental. It justifies that under the assumption of NIW prior, the update is a weighted sum of previous expected mean and covariance. It is striking that it provides very similar formulae to the standard CMA ES update. Recall that these updates given for the mean $m_{t}$ and covariance $C_{t}$ can be written as follows:

[TABLE]

where the notations $m_{t},w_{i},x_{i:\lambda},C_{t},c_{1},c_{\mu},c_{s},\mathrm{etc}...$ are given for instance in (Wikipedia, 2018).

Proof.

See 6.3 in the appendix section. ∎

3.3. Particularities of Bayesian CMA-ES

There are some subtleties that need to be emphasized.

• •

Although we assume a prior, we do not need to simulate the prior but can at each step use the expected value of the prior which means that we do not consume additional simulation compared to the standard CMA-ES.

• •

We need to tackle local minimum (we will give example of this in the numerical section) to avoid being trapped in a bowl! If we are in a local minimum, we need to inflate the variance to increase our search space. We do this whenever our algorithm does not manage to decrease. However, if after a while we do not get better result, we assume that this is indeed not a local minimum but rather a global minimum and start deflating the variance. This mechanism of inflation deflation ensures we can handle noisy functions like Rastrigin or Schwefel 1 or Schwefel 2 functions as defined in the section 4.

3.4. Differences with standard CMA-ES

Since we use a rigorous derivation of the posterior, we have the following features:

• •

the update of the covariance takes all points. This is different from $\lambda/\mu$ CMA-ES that uses only a subset of the point.

• •

by design, the update is optimal as we compute at each step the posterior.

• •

the contraction dilatation mechanism is an alternative to global local search path in standard CMA-ES.

• •

weights varies across iterations which is also a major difference between main CMA ES and Bayesian CMA ES. Weights are proportional to $\exp(\frac{1}{2}X^{T}\Sigma^{-1}X)$ sorted in decreasing order. Initially, when the variance is large,

3.5. Full algorithm

The complete Bayesian CMA ES algorithm is summarized in 2. It iterates until a stopping condition is met. We use multiple stopping conditions. We stop if we have not increase our best result for a given number of iterations. We stop if we have reached the maximum of our iterations. We stop if our variance norm is small. Additional stopping condition can be incorporated easily.

Last but not least, we have a dilatation contraction mechanism for the variance to handle local minima with multiple level of contractions and dilatation that is given in function 3. The overall idea is first to dilate variance if we do not make any progress to increase the search space so that we are not trapped in a local minimum. Should this not succeed, it means that we are reaching something that looks like the global minimum and we progressively contract the variance. In our implemented algorithm, we take $L_{1}=5,L_{2}=20,L_{3}=30,L_{4}=40,L_{5}=50$ and the dilatation, contraction parameters given by $k_{1}=1.5,k_{2}=0.9,k_{3}=0.7,k_{5}=0.5$ We have also a restart at previous minimum level $L_{*}=L_{2}$.

4. Numerical results

4.1. Functions examined

We have examined four functions to stress test our algorithm. They are listed in increasing order of complexity for our algorithm and correspond to different type of functions. They are all generalized function that can defined for any dimension $n$. For all, we present the corresponding equation for a variable $x=(x_{1},x_{2},..,x_{n})$ of $n$ dimension. Code is provided in supplementary materials. We have frozen seeds to have reproducible results.

4.1.1. Cone

The most simple function to optimize is the quadratic cone whose equation is given by (14) and represented in figure 2. It is also the standard Euclidean norm. It is obviously convex and is a good test of the performance of an optimization method.

[TABLE]

4.1.2. Schwefel 2 function

A slightly more complicated function is the Schwefel 2 function whose equation is given by (15) and represented in figure 3. It is a piecewise linear function and validates the algorithm can cope with non convex function.

[TABLE]

4.1.3. Rastrigin

The Rastrigin function, first proposed by (Rastrigin, 1974) and generalized by (Mühlenbein et al., 1991), is more difficult compared to the Cone and the Schwefel 2 function. Its equation is given by (16) and represented in figure 4. It is a non-convex function often used as a performance test problem for optimization algorithms. It is a typical example of non-linear multi modal function. Finding its minimum is considered a good stress test for an optimization algorithm, due to its large search space and its large number of local minima.

[TABLE]

4.1.4. Schwefel 1 function

The last function we tested is the Schwefel 1 function whose equation is given by (4.1.4) and represented in figure 5. It is sometimes only defined on $\left[-500,500\right]^{n}$. The Schwefel 1 function shares similarities with the Rastrigin function. It is continuous, not convex, multi-modal and with a large number of local minima. The extra difficulty compared to the Rastrigin function, the local minima are more pronounced local bowl making the optimization even harder.

[TABLE]

4.2. Convergence

For each of the functions, we compared our method using strategy one entitled B-CMA-ES S1: update $\overline{\mu}_{t}$ and $\overline{\Sigma}_{t}$ using (8) and (10) plotted in orange, or strategy two B-CMA-ES S2: same update but using (9) and (10), plotted in blue and standard CMA-ES as provided by the opensource python package pycma plotted in green. We clearly see that strategies one and two are quite similar to standard CMA-ES. The convergence graphics that show the error compared to the minimum are represented:

• •

for the cone function by figure 6 (case of a convex function), with initial point $(10,10)$

• •

for the Schwefel 2 function in figure 7 (case of piecewise linear function), with initial point $(10,10)$

• •

for the Rastrigin function in figure 8 (case of a non convex function with multiple local minima), with initial point $(10,10)$

• •

and for the Schwefel 1 function in figure 9 (case of a non convex function with multiple large bowl local minima), with initial point $(400,400)$

The results are for one test run. In a forthcoming paper, we will benchmark them with more runs to validate the interest of this new method.

For the four functions, BCMAES achieves convergence similar to standard CMA-ES. The intuition of this good convergence is that shifting the multi variate mean by the best candidate seen so far is a good guess to update it at the next run (standard CMA-ES or B-CMA-ES S1).

5. Conclusion

In this paper, we have revisited the CMA-ES algorithm and provided a Bayesian version of it. Taking conjugate priors, we can find optimal update for the mean and covariance of the multi variate Normal. We have provided the corresponding algorithm that is a new version of CMA-ES. First numerical experiments show this new version is competitive to standard CMA-ES on traditional functions such as cone, Schwefel 1, Rastrigin and Schwefel 2. This faster convergence can be explained on a theoretical side from an optimal update of the prior (thanks to Bayesian update) and the use of the best candidate seen at each simulation to shift the mean of the multi-variate Gaussian likelihood. We envisage further works to benchmark our algorithm to other standard evolutionary algorithms, in particular to use the COCO platform to provide more meaningful tests and confirm the theoretical intuition of good performance of this new version of CMA-ES, and to test the importance of the prior choice.

6. Appendix

6.1. Conjugate priors

Proof.

Consider $n$ independent and identically distributed (IID) measurements $\mathcal{X}\triangleq\{{\mathbf{x}}^{j}\in\mathbb{R}^{d}|1\leq j\leq n\}$ and assume that these variables have an exponential family density. The likelihood $p(\mathcal{X}|\theta,\lambda)$ writes simply as the product of each individual likelihood:

[TABLE]

If we start with a prior $\pi(\theta)$ of the form $\pi(\theta)\propto\exp(\mathcal{F}(\theta))$ for some function $\mathcal{F}(\cdot)$, its posterior writes:

[TABLE]

It is easy to check that the posterior (19) is in the same exponential family as the prior iff $\mathcal{F}(\cdot)$ is in the form:

[TABLE]

for some $\mu\triangleq(\mu_{0},\mu_{1})$, such that:

[TABLE]

Hence, the conjugate prior for the likelihood (18) is parametrized by $\mu$ and given by:

[TABLE]

where $Z={\int{\exp\left(\mu_{1}\cdot\eta(\theta,\lambda)-\mu_{0}A(\eta(\theta,\lambda))\right)\;\mathrm{d}x}}$. ∎

6.2. Exact computation of the posterior update for the Normal inverse Wishart

To make our proof simple, we first start by the one dimensional case and show that in one dimension it is a normal inverse gamma. We then generalize to the multi dimensional case.

Lemma 6.1.

The probability density function of a Normal inverse gamma (denoted by NIG) random variable can be expressed as the product of a Normal and an Inverse gamma probability density functions.

Proof.

we suppose that $x|\mu,\sigma^{2}\sim\mathcal{N}(\mu_{0},\sigma^{2}/v)$. We recall the following definition of conditional probability:

Definition 6.0.

*Suppose that events A,B and C are defined on the same probability space, and the event B is such that $\mathbb{P}(B)>0$. We have the following expression:

$\mathbb{P}(A\cap B|C)=\mathbb{P}(A|B,C)\mathbb{P}(B|C)$.*

Applying 6.2, we have:

[TABLE]

Using the definition of the Normal inverse gamma law, we end the proof. ∎

Remark 6.1.

If $\left(x,\sigma^{2}\right)\sim NIG(\mu,\lambda,\alpha,\beta)$, the probability density function is the following:

[TABLE]

Proposition 6.3.

The Normal Inverse Gamma NIG $\left(\mu_{0},v,\alpha,\beta\right)$ distribution is a conjugate prior of a normal distribution with unknown mean and variance.

Proof.

the posterior is proportional to the product of the prior and likelihood, then:

[TABLE]

Defining the empirical mean and variance as $\overline{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$ and $\overline{s}=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}$, we obtain that $\sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2}=n(\overline{s}+(\overline{x}-\mu)^{2})$.

So, the conditional density writes:

[TABLE]

Besides,

[TABLE]

Denoting $a=v+n$ and $b=v\mu_{0}+n\overline{x}$, we have :

[TABLE]

So we can express the proportional expression of the posterior :

[TABLE]

with

• •

$\alpha^{\star}=\alpha+\frac{n}{2}$

• •

$\beta^{\star}=\beta+\frac{1}{2}\left(\sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)^{2}+\frac{nv}{n+v}\frac{(\overline{x}-\mu_{0})^{2}}{2}\right)$

• •

$\mu^{\star}=\frac{v\mu_{0}+n\overline{x}}{v+n}$

• •

$\lambda^{\star}=v+n$

We can identify the terms with the expression of the probability density function given in 6.1 to conclude that the posterior follows a NIG$(\mu^{\star},\lambda^{\star},\alpha^{\star},\beta^{\star})$. ∎

We are now ready to prove the following proposition:

Proposition 6.4.

The Normal Inverse Wishart (denoted by NIW) $\left(\mu_{0},\kappa_{0},v_{0},\psi\right)$ distribution is a conjugate prior of a multivariate normal distribution with unknown mean and covariance.

Proof.

we use the fact that the probability density function of a Normal inverse Wishart random variable can be expressed as the product of a Normal and an Inverse Wishart probability density functions (we use the same reasoning that in 6.1). Besides, the posterior is proportional to the product of the prior and the likelihood.

We first express the probability density function of the multivariate Gaussian random variable in a proper way in order to use it when we write the posterior density function.

[TABLE]

We can inject the previous result and use the properties of the trace function to express the following probability density function of the multivariate Gaussian random variable of parameters $\mu$ and $\Sigma$. The density writes as:

[TABLE]

Hence, we can compute explicitly the posterior as follows:

[TABLE]

[TABLE]

We organize the terms and find the parameters of our Normal Inverse Wishart random variable NIW$(\mu_{0}^{\star},\kappa_{0}^{\star},v_{0}^{\star},\psi^{\star})$.

[TABLE]

which are exactly the equations provided in (5). ∎

6.3. Weighted combination for the BCMA ES update

Proof.

[TABLE]

[TABLE]

$\overline{\Sigma}_{t}$ is a covariance matrix of rank $n-1$ as we subtract the empirical mean (which removes one degree of freedom). The matrix ${\left(\hat{\mu}-\mathbb{E}_{t}[\mu]\right)\left(\hat{\mu}-\mathbb{E}_{t}[\mu]\right)^{T}}$ is of rank 1 as it is parametrized by the vector $\hat{\mu}$. ∎