Eric Benhamou
email: eric.benhamou@aisquareconnect.com

David Saltiel
email: david.saltiel@aisquareconnect.com

Sébastien Verel
email: verel@univ-littoral.fr

Fabien Teytaud
email: fabien.teytaud@univ-littoral.fr

Eric Benhamou Lamsade

Sebastien Verel

Bayesian CMA-ES: a new approach

Keywords: Mathematics of computing → Probability and statistics; CMA-ES, Bayesian, conjugate prior, normal-inverse-Wishart

come

INTRODUCTION

The covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] 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 [1, 2, 4-6, 10-12, 15, 21] and [START_REF] Varelas | A Comparative Study of Large-Scale Variants of CMA-ES[END_REF] 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 ( [START_REF] Heidrich | Neuroevolution strategies for episodic reinforcement learning[END_REF][START_REF] Igel | Evolutionary Kernel Learning[END_REF][START_REF] Igel | Shark[END_REF]), and references therein, as well as hyperparameter optimization of deep neural networks [START_REF] Loshchilov | CMA-ES for Hyperparameter Optimization of Deep Neural Networks[END_REF] In a nutshell, the (𝜇 / 𝜆) CMA-ES is an iterative black box optimization algorithm, that, in each of its iterations, samples 𝜆 candidate solutions from a multivariate normal distribution, evaluates these solutions (sequentially or in parallel) retains 𝜇 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 [START_REF] Akimoto | Bidirectional Relation between CMA Evolution Strategies and Natural Evolution Strategies[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF] and [START_REF] Ollivier | Information-geometric Optimization Algorithms: A Unifying Picture via Invariance Principles[END_REF] 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 (𝜇/𝜆) CMA-ES as out of the 𝜆 generated paths, only the 𝜇 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 𝑝 𝐶 𝑝 𝑇 𝐶 with parameter 𝑐 1 and a rank 𝑚𝑖𝑛(𝜇, 𝑛) matrix with parameter 𝑐 𝜇 , whose details are given for instance in [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF]. 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.

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.

Bayesian vs Frequentist probability theory

The justification of the Bayesian approach is discussed in [START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. 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 [START_REF] Gelman | Bayesian Data Analysis[END_REF][START_REF] Marin | Bayesian Core: A Practical Approach to Computational Bayesian Statistics[END_REF][START_REF] Robert | The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation[END_REF]. 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 [START_REF] Jordan | Lecture notes: Justification for Bayes[END_REF]. Definition 2.1. (Infinite exchangeability). We say that (𝑥 1 , 𝑥 2 , ...) is an infinitely exchangeable sequence of random variables if, for any n, the joint probability 𝑝 (𝑥 1 , 𝑥 2 , ..., 𝑥 𝑛 ) is invariant to permutation of the indices. That is, for any permutation 𝜋,

𝑝 (𝑥 1 , 𝑥 2 , ..., 𝑥 𝑛 ) = 𝑝 (𝑥 𝜋 1 , 𝑥 𝜋 2 , ..., 𝑥 𝜋𝑛 )
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.1. (De Finetti, 1930s). A sequence of random variables (𝑥 1 , 𝑥 2 , ...) is infinitely exchangeable iff, for all n,

𝑝 (𝑥 1 , 𝑥 2 , ..., 𝑥 𝑛 ) = ∫ 𝑛 𝑖=1 𝑝 (𝑥 𝑖 |𝜃 )𝑃 (𝑑𝜃 ),
for some measure P on 𝜃 . This representation theorem 2.1 justifies the use of priors on parameters since for exchangeable data, there must exist a parameter 𝜃 , a likelihood 𝑝 (𝑥 |𝜃 ) and a distribution 𝜋 on 𝜃 . A proof of De Finetti theorem is for instance given in [START_REF] Schervish | Theory of Statistics[END_REF] (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 𝜃 . The sequence may not be independent neither identically distributed, which is a much stronger result!

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

Posterior ∝ Prior × Likelihood
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 𝑋 , its distribution 𝑝 depends on a parameter 𝜃 that can be multi-dimensional. To emphasize the dependency of the distribution on the parameters, let us write this distribution as 𝑝 (𝑥 |𝜃 ) and let us assume we have access to a prior distribution 𝜋 (𝜃 ). Then the joint distribution of (𝜃, 𝑥) writes simply as

𝜙 (𝜃, 𝑥) = 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 )
The marginal distribution of 𝑥 is trivially given by marginalizing the joint distribution by 𝜃 as follows:

𝑚(𝑥) = ∫ 𝜙 (𝜃, 𝑥)𝑑𝜃 = ∫ 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 )𝑑𝜃
The posterior of 𝜃 is obtained by Bayes's formula as

𝜋 (𝜃 |𝑥) = 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 ) ∫ 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 )𝑑𝜃 ∝ 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 )
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.2. A prior distribution 𝜋 (𝜃 ) is said to be a conjugate prior if the posterior distribution 𝜋 (𝜃 |𝑥) ∝ 𝑝 (𝑥 |𝜃 )𝜋 (𝜃 )

(1) 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.3. A distribution is said to belong to the exponential family if it can be written (in its canonical form) as:

𝑝 (x|𝜂) = ℎ(x) exp(𝜂 • 𝑇 (x) -𝐴(𝜂)), ( 2
)
where 𝜂 is the natural parameter, 𝑇 (x) is the sufficient statistic, 𝐴(𝜂) is log-partition function and ℎ(x) is the base measure. 𝜂 and 𝑇 (x) may be vector-valued. Here 𝑎 • 𝑏 denotes the inner product of 𝑎 and 𝑏. The log-partition function is defined by the integral

𝐴(𝜂) ≜ log ∫ X ℎ(x) exp(𝜂 • 𝑇 (x)) d𝑥 . ( 3
)
2023-04-18 13:30. Page 2 of 1-11.

Also, 𝜂 ∈ Ω = {𝜂 ∈ R 𝑚 |𝐴(𝜃 ) < +∞} where Ω is the natural parameter space. Moreover, Ω is a convex set and 𝐴(•) is a convex function on Ω.

Remark 2.2. Not surprisingly, the normal distribution N (x; 𝜇, Σ) with mean 𝜇 ∈ R 𝑑 and covariance matrix Σ belongs to the exponential family but with a different parametrisation. Its exponential family form is given by:

𝜂 (𝜇, Σ) = Σ -1 𝜇 vec(Σ -1 ) , 𝑇 (x) = x vec(-1 2 xx T ) , (4a)
ℎ(x) = (2𝜋) -𝑑 2 , 𝐴(𝜂 (𝜇, Σ)) = 1 2 𝜇 T Σ -1 𝜇 + 1 2 log |Σ|. ( 4b
)
where in equations (4a), the notation vec(•) means we have vectorized the matrix, stacking each column on top of each other and hence can equivalently write for 𝑎 and 𝑏, two matrices, the trace result Tr(𝑎 T 𝑏) as the scalar product of their vectorization vec(𝑎) • vec(𝑏). We can remark the canonical parameters are very different from traditional (also called moment) parameters. We can notice that changing slightly the sufficient statistic 𝑇 (𝑥) leads to change the corresponding canonical parameters 𝜂.

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

Proposition 2.2. If the observations have a density of the exponential family form 𝑝 (𝑥 |𝜃, 𝜆) = ℎ(𝑥) exp 𝜂 (𝜃, 𝜆) 𝑇 𝑇 (𝑥) -𝑛𝐴(𝜂 (𝜃, 𝜆)) , with 𝜆 a set of hyper-parameters, then the prior with likelihood defined by

𝜋 (𝜃 ) ∝ exp (𝜇 1 • 𝜂 (𝜃, 𝜆) -𝜇 0 𝐴(𝜂 (𝜃, 𝜆))) with 𝜇 ≜ (𝜇 0 , 𝜇 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 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.

Figure 1: As we get more and more information using the likelihood, the posterior becomes more peak.

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 normalinverse-Wishart to keep things simple. The Normal-inverse-Wishart distribution is parametrized by 𝝁 0 , 𝜆, 𝚿, 𝜈 and its distribution is given by

𝑓 (𝝁, 𝚺|𝝁 0 , 𝜆, 𝚿, 𝜈) = N 𝝁 𝝁 0 , 1 𝜆 𝚺 W -1 (𝚺|𝚿, 𝜈)
where 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.3. If our sampling density follows a 𝑑 dimensional multivariate normal distribution ∼ N 𝑑 (𝝁, 𝚺) with unknown mean 𝜇 and covariance Σ and if its parameters are distributed according to a Normal-Inverse-Wishart (𝜇, Σ) ∼ NIW(𝜇 0 , 𝜅 0 , 𝑣 0 ,𝜓 ) and if we observe X = (𝑥 1 , .., 𝑥 𝑛 ) samples, then the posterior is also a Normal-Inverse-Wishart with different parameters NIW(𝜇 ★ 0 , 𝜅 ★ 0 , 𝑣 ★ 0 ,𝜓 ★ ) given by

𝜇 ★ 0 = 𝜅 0 𝜇 0 + 𝑛𝑥 𝜅 0 + 𝑛 , 𝜅 ★ 0 = 𝜅 0 + 𝑛, 𝑣 ★ 0 = 𝑣 0 + 𝑛 𝜓 ★ = 𝜓 + 𝑛 𝑖=1 (𝑥 𝑖 -𝑥) (𝑥 𝑖 -𝑥) 𝑇 + 𝜅 0 𝑛 𝜅 0 + 𝑛 (𝑥 -𝜇 0 ) (𝑥 -𝜇 0 ) 𝑇 (5)
with 𝑥 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. □

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 𝑓 follows a multivariate normal law.

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 𝑓 for candidate points. At step 𝑡 in our algorithm, we suppose multi variate Gaussian parameters 𝜇 and Σ follow a normal inverse Wishart denoted by 𝑁 𝐼𝑊 (𝜇 𝑡 , 𝜅 𝑡 , 𝑣 𝑡 ,𝜓 𝑡 ).

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:

E 𝑡 [𝜇] = 𝜇 𝑡 (6
) and . This comes from three things:

E 𝑡 [Σ] = 𝜓 𝑡 𝑣 𝑡 -𝑛 -1 (7)
• 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 𝑓 . • 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: E 𝑡 [𝜇] and E 𝑡 [Σ] 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 𝑋 '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 𝑋 𝑖 , we compute it assumed density denoted by

𝑑 𝑖 = N (E 𝑡 [𝜇], E 𝑡 [Σ])(𝑋 𝑖 ) where N (E 𝑡 [𝜇], E 𝑡 [Σ])(.
) denotes the p.d.f. of the multi-variate Gaussian.

We divide these density by their sum to get weights (𝑤 𝑖 ) 𝑖=1..𝑘 that are positive and sum to one as follows. 𝑤 𝑗 = 𝑑 𝑗 / 𝑘 𝑖=1 𝑑 𝑖 . Hence for 𝑘 simulated points, we get {𝑋 𝑖 , 𝑤 𝑖 } 𝑖=1..𝑘 . 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 (𝑤 𝑖 ) 𝑖=1..𝑘 , we denote the order statistics (𝑖), 𝑤 ↓. This first sorting leads to k new uplets {𝑋 (𝑖),𝑤↓ , 𝑤 (𝑖),𝑤↓ } 𝑖=1..𝑘 . 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 {𝑋 (𝑖),𝑓 ↑ , 𝑤 (𝑖),𝑤↓ } 𝑖=1..𝑘 . We can now compute the empirical mean 𝜇 𝑡 as follows:

𝜇 𝑡 = 𝑘 𝑖=1 𝑤 (𝑖),𝑤↓ • 𝑋 (𝑖),𝑓 ↑ MC mean for 𝑋 𝑓 ↑ - 𝑘 𝑖=1 𝑤 𝑖 𝑋 𝑖 -𝜇 𝑡 MC bias for 𝑋 (8)
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 1 2 (𝑋 -

E 𝑡 [𝜇]) 𝑇 (E 𝑡 [Σ]) -1 (𝑋 -E 𝑡 [𝜇]
) . When we start the algorithm, we use a large search space, hence a large covariance matrix Σ 𝑡 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 𝑓 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:

𝜇 𝑡 = arg min 𝑋 ∈X 𝑓 (𝑋 ) (9)
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:

•

Σ 𝑡 = 𝑘 𝑖=1 𝑤 (𝑖),𝑤↓ • 𝑋 (𝑖),𝑓 ↑ -𝑋 (.),𝑓 ↑ 𝑋 (𝑖),𝑓 ↑ -𝑋 (.),𝑓 ↑ 𝑇 MC covariance for 𝑋 𝑓 ↑ - 𝑘 𝑖=1 𝑤 𝑖 • 𝑋 𝑖 -𝑋 𝑋 𝑖 -𝑋 𝑇 -Σ 𝑡 MC covariance for simulated 𝑋 (10)
where 𝑋 (.),𝑓 ↑ = 𝑘 𝑖=1 𝑤 (𝑖),𝑤↓ 𝑋 (𝑖),𝑓 ↑ and 𝑋 = 𝑘 𝑖=1 𝑤 𝑖 𝑋 𝑖 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 [START_REF] Hansen | CMA-ES: evolution strategies and covariance matrix adaptation[END_REF], we update the posterior law according to equation [START_REF] Auger | Tutorial CMA-ES: evolution strategies and covariance matrix adaptation[END_REF]. This gives us the iterative conjugate prior parameters updates:

𝜇 𝑡 +1 = 𝜅 𝑡 𝜇 𝑡 + 𝑛𝜇 𝑡 𝜅 𝑡 + 𝑛 , 𝜅 𝑡 +1 = 𝜅 𝑡 + 𝑛, 𝑣 𝑡 +1 = 𝑣 𝑡 + 𝑛, 𝜓 𝑡 +1 = 𝜓 𝑡 +Σ 𝑡 + 𝜅 𝑡 𝑛 𝜅 𝑡 + 𝑛 𝜇 𝑡 -𝜇 𝑡 𝜇 𝑡 -𝜇 𝑡 𝑇 (11)
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

E 𝑡 +1 [𝜇] = E 𝑡 [𝜇] + 𝑤 𝜇 𝑡 𝜇 𝑡 -E 𝑡 [𝜇] , E 𝑡 +1 [Σ] = 𝑤 Σ,1 𝑡 discount factor E 𝑡 [Σ] + 𝑤 Σ,2 𝑡 𝜇 𝑡 -E 𝑡 [𝜇] 𝜇 𝑡 -E 𝑡 [𝜇] 𝑇 rank one matrix + 𝑤 Σ,3 𝑡 Σ 𝑡 rank (n-1) matrix
where 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 𝑚 𝑡 and covariance 𝐶 𝑡 can be written as follows:

𝑤 𝜇 𝑡 = 𝑛 𝜅 𝑡 + 𝑛 , 𝑤 Σ,1 𝑡 = 𝜅 𝑡 𝑛 (𝜅 𝑡 + 𝑛)(𝑣 𝑡 -1) 𝑤 Σ,2 𝑡 = 𝑣 𝑡 -𝑛 -1 𝑣 𝑡 -1 , 𝑤 Σ,3 𝑡 = 1 𝑣 𝑡 -1 (12
𝑚 𝑡 +1 = 𝑚 𝑡 + 𝜇 𝑖=1 𝑤 𝑖 (𝑥 𝑖:𝜆 -𝑚 𝑡 ) 𝐶 𝑡 +1 = (1 -𝑐 1 -𝑐 𝜇 + 𝑐 𝑠 ) discount factor 𝐶 𝑡 +

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.

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 𝜆/𝜇 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( 1 2 𝑋 𝑇 Σ -1 𝑋 ) sorted in decreasing order. Initially, when the variance is large,

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. end if 12: end while 13: return best solution 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 𝐿 1 = 5, 𝐿 2 = 20, 𝐿 3 = 30, 𝐿 4 = 40, 𝐿 5 = 50 and the dilatation, contraction parameters given by 𝑘 1 = 1.5, 𝑘 2 = 0.9, 𝑘 3 = 0.7, 𝑘 5 = 0.5 We have also a restart at previous minimum level 𝐿 * = 𝐿 2 .

NUMERICAL RESULTS

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 𝑛. For all, we present the corresponding equation for a variable 𝑥 = (𝑥 1 , 𝑥 2 , .., 𝑥 𝑛 ) of 𝑛 dimension. Code is provided in supplementary materials. We have frozen seeds to have reproducible results. Set 𝑓 𝑚𝑖𝑛 = 𝑓 (𝑋 )

4:
Memorize current point and its variance:

5:

• 𝑋 𝑚𝑖𝑛 = 𝑋 6:

• Σ 𝑡,𝑚𝑖𝑛 = Σ 𝑡 7:

Set retrial = 0 8: else Restart at previous best solution:

12:

• 𝑋 = 𝑋 𝑚𝑖𝑛 13:

• Σ 𝑡 = Σ 𝑡,𝑚𝑖𝑛 [START_REF] Rastrigin | Systems of extremal control[END_REF] and generalized by [START_REF] Mühlenbein | The parallel genetic algorithm as function optimizer[END_REF], 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.

Convergence

For each of the functions, we compared our method using strategy one entitled B-CMA-ES S1: update 𝜇 𝑡 and Σ 𝑡 using ( 8) and ( 10) plotted in orange, or strategy two B-CMA-ES S2: same update but using ( 9) and ( 10 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).

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

If we start with a prior 𝜋 (𝜃 ) of the form 𝜋 (𝜃 ) ∝ exp(F (𝜃 )) for some function F (•), its posterior writes:

𝜋 (𝜃 |X) ∝ 𝑝 (X|𝜃 ) exp(F (𝜃 )) ∝ exp 𝜂 (𝜃, 𝜆) • 𝑛 𝑗=1 𝑇 (𝑥 𝑗 ) -𝑛𝐴(𝜂 (𝜃, 𝜆)) + F (𝜃 ) . ( 19
)
It is easy to check that the posterior ( 19) is in the same exponential family as the prior iff F (•) is in the form:

F (𝜃 ) = 𝜇 1 • 𝜂 (𝜃, 𝜆) -𝜇 0 𝐴(𝜂 (𝜃, 𝜆)) (20)
for some 𝜇 ≜ (𝜇 0 , 𝜇 1 ), such that:

𝑝 (X|𝜃, 𝜆) ∝ exp 𝜇 1 + 𝑛 𝑗=1 𝑇 (𝑥 𝑗 ) 𝑇 𝜂 (𝜃, 𝜆) -(𝑛 + 𝜇 0 )𝐴(𝜂 (𝜃, 𝜆)) . (21)
Hence, the conjugate prior for the likelihood ( 18) is parametrized by 𝜇 and given by:

𝑝 (X|𝜃, 𝜆) = 1 𝑍 exp (𝜇 1 • 𝜂 (𝜃, 𝜆) -𝜇 0 𝐴(𝜂 (𝜃, 𝜆))) , (22)
where 𝑍 = ∫ exp (𝜇 1 • 𝜂 (𝜃, 𝜆) -𝜇 0 𝐴(𝜂 (𝜃, 𝜆))) d𝑥. □

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 𝑥 |𝜇, 𝜎 2 ∼ N (𝜇 0 , 𝜎 2 /𝑣). We recall the following definition of conditional probability: Definition 6.1. Suppose that events A,B and C are defined on the same probability space, and the event B is such that P(𝐵) > 0. We have the following expression: P(𝐴 ∩ 𝐵|𝐶) = P(𝐴|𝐵, 𝐶)P(𝐵|𝐶).

Using the definition of the Normal inverse gamma law, we end the proof. □ Remark 6.1. If 𝑥, 𝜎 2 ∼ 𝑁 𝐼𝐺 (𝜇, 𝜆, 𝛼, 𝛽), the probability density function is the following:

𝑓 (𝑥, 𝜎 2 |𝜇, 𝜆, 𝛼, 𝛽) = √ 𝜆 𝜎 √ 2𝜋 𝛽 𝛼 Γ(𝛼) 1 𝜎 2 𝛼+1 exp - 2𝛽 + 𝜆(𝑥 -𝜇) 2 2𝜎 2 . ( 24
)
Proposition 6.2. The Normal Inverse Gamma NIG (𝜇 0 , 𝑣, 𝛼, 𝛽) 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:

𝑝 𝜇, 𝜎 2 |𝑋 ∝ √ 𝑣 √ 2𝜋 1 𝜎 2 1/2 exp -𝑣 (𝜇 -𝜇 0 ) 2 2𝜎 2 × 𝛽 𝛼 Γ(𝛼) 1 𝜎 2 𝛼+1 exp -𝛽 𝜎 2 × 1 2𝜋𝜎 2 𝑛/2 exp - 𝑛 𝑖=1 (𝑥 𝑖 -𝜇) 2 2𝜎 2 . ( 25
)
Defining the empirical mean and variance as 𝑥

= 1 𝑛 𝑛 𝑖=1 𝑥 𝑖 and 𝑠 = 1 𝑛 𝑛 𝑖=1 (𝑥 𝑖 -𝑥) 2 , we obtain that 𝑛 𝑖=1 (𝑥 𝑖 -𝜇) 2 = 𝑛(𝑠 + (𝑥 -𝜇) 2 ).
So, the conditional density writes:

𝑝 𝜇, 𝜎 2 |𝑋 ∝ √ 𝑣 1 𝜎 2 𝛼+𝑛/2+3/2 × exp - 1 𝜎 2 𝛽 + 1 2 𝑣 (𝜇 -𝜇 0 ) 2 + 𝑛 𝑠 + (𝑥 -𝜇) 2 . (26)
Besides, • 𝜇 ★ = 𝑣𝜇 0 +𝑛𝑥 𝑣+𝑛

𝑣 (𝜇 -𝜇 0 ) 2 + 𝑛 𝑠 + (𝑥 -𝜇) 2 = 𝑣 𝜇 2 -2𝜇𝜇 0 + 𝜇 2 𝑜 + 𝑛𝑠 + 𝑛 𝑥 2 -2𝑥 𝜇 + 𝜇 2 = 𝜇 2 (𝑣 + 𝑛) -2𝜇 (𝑣𝜇 0 + 𝑛𝑥) + 𝑣𝜇 2 𝑜 + 𝑛𝑠 + 𝑛𝑥 2 . ( 27
• 𝜆 ★ = 𝑣 + 𝑛
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(𝜇 ★ , 𝜆 ★ , 𝛼 ★ , 𝛽 ★ ). □

We are now ready to prove the following proposition: Proposition 6.3. The Normal Inverse Wishart (denoted by NIW) (𝜇 0 , 𝜅 0 , 𝑣 0 ,𝜓 ) 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.

We organize the terms and find the parameters of our Normal Inverse Wishart random variable NIW(𝜇 ★ 0 , 𝜅 ★ 0 , 𝑣 ★ 0 ,𝜓 ★ ). .

  ) 2023-04-18 13:30. Page 4 of 1-11.

  : A Bayesian approach to CMA-ES Conference'17, July 2017, Washington, DC, USA Remark 3.1.

Algorithm 2 8 : 9 :

 289 Bayesian update of CMA-ES parameters: 1: Initialization 2: Start with a prior distribution Π on 𝜇 and Σ 3: Set retrial to 0 4: Set 𝑓 𝑚𝑖𝑛 to max float 5: while stop criteria not satisfied do 6: 𝑋 ∼ N (𝜇, Σ) 7: update the parameters of the Gaussian thanks to the posterior law Π(𝜇, Σ|𝑋 ) following details given in algorithm 1 Handle dilatation contraction variance for local minima as explained in algorithm 3 if DilateContractFunc(𝑋, Σ 𝑡 , 𝑋 𝑚𝑖𝑛 , 𝑓 𝑚𝑖𝑛 , Σ 𝑡,𝑚𝑖𝑛 ) == 1 then 10:

  2023-04-18 13:30. Page 5 of 1-11.

Conference' 17 , 1 :

 171 July 2017, Washington, DC, USA Benhamou, Saltiel, Verel, Teytaud Algorithm 3 Dilatation contraction variance for local minima: Function DilateContractFunc(𝑋, Σ 𝑡 , 𝑋 𝑚𝑖𝑛 , 𝑓 𝑚𝑖𝑛 , Σ 𝑡,𝑚𝑖𝑛 ) 2: if 𝑓 (𝑋 ) ≤ 𝑓 𝑚𝑖𝑛 then 3:

  if retrial == 𝐿 * then 11:

20 :Figure 2 :Figure 3 :

 2023 Figure 2: A simple convex function: the quadratic norm. Minimum in 0

Figure 4 :

 4 Figure 4: Rastrigin function: a non convex function multimodal and with a large number of local minima

.

Figure 5 :

 5 Figure 5: Schwefel 1 function: a non convex function multimodal and with a large number of local pronounced bowls

  ), 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 figure6(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)

Figure 6 :

 6 Figure 6: Convergence for the Cone function

Figure 7 :

 7 Figure 7: Convergence for the Schwefel 2 function

Figure 8 :

 8 Figure 8: Convergence for the Rastrigin function

6 APPENDIX 6 . 1

 661 Conjugate priorsProof. Consider 𝑛 independent and identically distributed (IID) measurements X ≜ {x 𝑗 ∈ R 𝑑 |1 ≤ 𝑗 ≤ 𝑛} and assume that these variables have an exponential family density. The likelihood 2023-04-18 13:30. Page 7 of 1-11.

.

Conference' 17 ,Figure 9 :

 179 Figure 9: Convergence for the Schwefel 1 function

Applying 6 . 1 ,

 61 we have: 𝑝 𝜇, 𝜎 2 |𝜇 0 , 𝑣, 𝛼, 𝛽 = 𝑝 𝜇|𝜎 2 , 𝜇 0 , 𝑣, 𝛼, 𝛽 𝑝 𝜎 2 |𝜇 0 , 𝑣, 𝛼, 𝛽 = 𝑝 𝜇|𝜎 2 , 𝜇 0 , 𝑣 𝑝 𝜎 2 |𝛼, 𝛽 .

  ) 2023-04-18 13:30. Page 8 of 1-11.

𝑛𝑖=1(𝜇× |𝜓 | 𝑣/ 2 2 2 𝑒𝑥𝑝 - 1 2 𝑡𝑟 𝜓 Σ - 1 ×( 2 𝜇

 221212 𝑥 𝑖 -𝜇) 𝑇 Σ -1 (𝑥 𝑖 -𝜇) = 𝑛 (𝑥 -𝜇) 𝑇 Σ -1 (𝑥 -𝜇) + 𝑛 𝑖=1 (𝑥 𝑖 -𝑥) 𝑇 Σ -1 (𝑥 𝑖 -𝑥) . (29)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 𝜇 and Σ.The density writes as:|Σ | -𝑛/2 √ (2𝜋 ) 𝑝𝑛 exp -𝑛 2 (𝑥 -𝜇) 𝑇 Σ -1 (𝑥 -𝜇) -1 2 𝑡𝑟 Σ -1 𝑛 𝑖=1 (𝑥 𝑖 -𝑥) (𝑥 𝑖 -𝑥) 𝑇 . (30)Hence, we can compute explicitly the posterior as follows:-𝜇 0 𝑇 Σ -1 𝜇 -𝜇 0 𝑣𝑝/2 Γ 𝑝 (𝑣 0 /2) |Σ| -𝑣 0 +𝑝+1 |Σ| -𝑛/2 exp -𝑛 2 (𝑥 -𝜇) 𝑇 Σ -1 (𝑥 -𝜇) 𝑥 𝑖 -𝑥) (𝑥 𝑖 -𝑥) 𝑇 -𝜇 0 𝑇 Σ -1 𝜇 -𝜇 0 -𝑛 2 (𝑥 -𝜇) 𝑇 Σ -1 (𝑥 -𝜇) -1 2 𝑡𝑟 Σ -1 𝜓 + 𝑛 𝑖=1(𝑥 𝑖 -𝑥) (𝑥 𝑖 -𝑥) 𝑇 .

𝜇 ★ 0 = 6 . 3 1 = 1 𝑣 𝑡 - 1 𝜓 𝑡 + 1 𝑣𝑣 𝑡 -𝑛 - 1 𝑣 𝑡 - 1 ,

 063111111 𝜅 0 𝜇 0 + 𝑛𝑥 𝜅 0 + 𝑛 , 𝜅 ★ 0 = 𝜅 0 + 𝑛, 𝑣 ★ 0 = 𝑣 0 + 𝑛 𝜓 ★ = 𝜓 + 𝑛 𝑖=1 (𝑥 𝑖 -𝑥) (𝑥 𝑖 -𝑥) 𝑇 + 𝜅 0 𝑛 𝜅 0 + 𝑛 (𝑥 -𝜇 0 ) (𝑥 -𝜇 0 ) 𝑇 (33)which are exactly the equations provided in[START_REF] Auger | Tutorial CMA-ES: evolution strategies and covariance matrix adaptation[END_REF].□ Weighted combination for the BCMA ES updateProof.E 𝑡 +1 [𝜇] = 𝜇 𝑡 +1 = 𝜅 𝑡 𝜇 𝑡 + 𝑛𝜇 𝑡 𝜅 𝑡 + 𝑛 = E 𝑡 [𝜇] + 𝑤 𝜇 𝑡 ( μ -E 𝑡 [𝜇]) (34) E 𝑡 +1 [Σ] = 𝜓 𝑡 +1 𝑣 𝑡 +1 -𝑛 -𝑡 -1 Σ 𝑡 + 𝜅 𝑡 𝑛 (𝜅 𝑡 + 𝑛)(𝑣 𝑡 -1) 𝜇 𝑡 -𝜇 𝑡 𝜇 𝑡 -𝜇 𝑡 𝑇 = 𝑤 Σ,1 𝑡 discount factor E 𝑡 [Σ] + 𝑤 Σ,2 𝑡 ( μ -E 𝑡 [𝜇]) ( μ -E 𝑡 [𝜇]) 𝑇 𝑤 Σ,2 𝑡 = 𝜅 𝑡 𝑛 (𝜅 𝑡 + 𝑛)(𝑣 𝑡 -1)(35)2023-04-18 13:30. Page 9 of 1-11.

  We simulate potential candidates X = {𝑋 𝑖 } ∼ N E 𝑡 [𝜇], E 𝑡 [Σ] and evaluate them 𝑓 (𝑋 𝑖 ). If the distribution of the minimum was accurate, the minimum would concentrate around E 𝑡 [𝜇] and be spread with a variance of E 𝑡[Σ]. When evaluating potential candidates, as our guess is not right, we do not get values centered around E 𝑡 [𝜇] and spread with a variance of E 𝑡[Σ]

		Conference'17, July 2017, Washington, DC, USA	Benhamou, Saltiel, Verel, Teytaud
	349
	350
	351
	352
	353
	354
	355
	356
	357
	358
	359
	360
	361
	362
	363
	364
	365
	366
	367
	368
	369
	370
	371
	372
	373
	374
	375
	376
		.
	377
	378
	379
	380
	381
	382
	383
	384
	385
	386
	387
	388
	389
	390
	391
	392
	393
	394
	395
	396
	397
	398
	399
	400
	401
	402
	403
	404
	405
	406	2023-04-18 13:30. Page 3 of 1-11.

  𝑐 1 𝑝 𝑐 𝑝 𝑇 𝑡 [𝜇] and Σ 𝑡 using (8) and (10) (strategy one) 11: Update 𝜇 𝑛+1 , 𝜅 𝑛+1 , 𝑣 𝑛+1 ,𝜓 𝑛+1 using[START_REF] Hansen | Evolution strategies and CMA-ES (covariance matrix adaptation)[END_REF]

					𝑐
				rank one matrix	(13)
	+ 𝑐 𝜇	𝜇	𝑤 𝑖	𝑥 𝑖:𝜆 -𝑚 𝑘	𝑥 𝑖:𝜆 -𝑚 𝑡	𝑇
		𝑖=1		𝜎 𝑘	𝜎 𝑡
			rank min(𝜇,𝑛-1) matrix

where the notations 𝑚 𝑡 , 𝑤 𝑖 , 𝑥 𝑖:𝜆 , 𝐶 𝑡 , 𝑐 1 , 𝑐 𝜇 , 𝑐 𝑠 , etc... are given for instance in [26]. Proof. See 6.3 in the appendix section. □ Algorithm 1 Predict and Correct parameters at step t 1: Simulate candidate 2: Use mean values E 𝑡 [𝜇] = 𝜇 𝑡 and Σ 𝑡 = E[Σ] = 𝜓 𝑡 /(𝑣 𝑡 -𝑛 -1) 3: Simulate k points X = {𝑋 𝑖 } = 1..𝑘 ∼ N (E 𝑡 [𝜇], Σ 𝑡 ) 4: Compute densities (𝑑 𝑖 ) 𝑖..𝑘 = (N (E 𝑡 [𝜇], Σ 𝑡 )(𝑋 𝑖 )) 𝑖..𝑘 = 5: Sort in decreasing order with respect to 𝑑 to get {𝑋 (𝑖),𝑑 ↓ , 𝑑 (𝑖),𝑑 ↓ } 𝑖=1..𝑘 6: Stable Sort in increasing order order with respect to 𝑓 (𝑋 𝑖 ) to get {𝑋 (𝑖),𝑓 ↑ , 𝑑 (𝑖),𝑑 ↓ } 𝑖=1..𝑘 7: 8: Correct E 𝑡 [𝜇] and Σ 𝑡 9: Either Update E 𝑡 [𝜇] and Σ 𝑡 using (9) and (10) (strategy two) 10: Or Update E

  if 𝐿 2 > retrial and retrial > 𝐿 1 then else if 𝐿 3 > retrial and retrial ≥ 𝐿 2 then

	14:	end if
	15:
	16:	Dilate variance by 𝑘 1
	17:
	18:	Contract variance by 𝑘 2

19:

else if 𝐿 4 > retrial and retrial ≥ 𝐿 3 then

  BCMA-ES: A Bayesian approach to CMA-ES Conference'17, July 2017, Washington, DC, USA Denoting 𝑎 = 𝑣 + 𝑛 and 𝑏 = 𝑣𝜇 0 + 𝑛𝑥, we have : 𝑣 (𝜇 -𝜇 0 ) 2 + 𝑛 𝑠 + (𝑥 -𝜇) 2

	𝛽 +	1
			2
	= 𝛽 +	1 2	𝑎𝜇 2 -2𝑏𝜇 + 𝑣𝜇 2 𝑜 + 𝑛𝑠 + 𝑛𝑥 2
	= 𝛽 +	1 2	𝑎 𝜇 2 -	2𝑏 𝑎	𝜇 + 𝑣𝜇 2 𝑜 + 𝑛𝑠 + 𝑛𝑥 2
	= 𝛽 +	1 2	𝑎 𝜇 -	𝑏 𝑎	2	-	𝑏 2 𝑎	+ 𝑣𝜇 2 𝑜 + 𝑛𝑠 + 𝑛𝑥 2 .	(28)
	So we can express the proportional expression of the posterior :
	𝑝 𝜇, 𝜎 2 |𝑋 ∝		1 𝜎 2	𝛼 ★ +3/2	× exp -	2𝛽 ★ + 𝜆 ★ 𝜇 -𝜇 ★ 2 2𝜎 2	,
	with
	• 𝛼 ★ = 𝛼 + 𝑛 2
	• 𝛽 ★ = 𝛽 + 1 2		𝑛 𝑖=1 (𝑥 𝑖 -𝑥) 2 + 𝑛𝑣 𝑛+𝑣	(𝑥-𝜇 0 ) 2 2