Alexandre Chotard

Anne Auger

Nikolaus Hansen

Markov Chain Analysis of Evolution Strategies on a Linear Constraint Optimization Problem

This paper analyses a (1, λ)-Evolution Strategy, a randomised comparison-based adaptive search algorithm, on a simple constraint optimization problem. The algorithm uses resampling to handle the constraint and optimizes a linear function with a linear constraint. Two cases are investigated: first the case where the step-size is constant, and second the case where the step-size is adapted using path length control. We exhibit for each case a Markov chain whose stability analysis would allow us to deduce the divergence of the algorithm depending on its internal parameters. We show divergence at a constant rate when the step-size is constant. We sketch that with step-size adaptation geometric divergence takes place. Our results complement previous studies where stability was assumed.

I. INTRODUCTION

Derivative Free Optimization (DFO) methods are tailored for the optimization of numerical problems in a black-box context, where the algorithms can only query the objective function to optimize f : R n → R, and no properties on f , such as convexity or differentiability, is exploited.

Evolution Strategies (ES) are comparison-based randomised DFO algorithms. At iteration t, solutions are sampled from a multivariate normal distribution centered in a vector X t . The candidate solutions are ranked according to f , and update of X t and other parameters of the distribution (usually a step-size σ t and a covariance matrix) is performed using the ranking information given by the candidate solutions. Since ES do not directly use the function values of the new points, but only how f ranks the different samples, they are invariant to the composition of the objective function by a strictly increasing function h : R → R.

This property and the black-box scenario make Evolution Strategies suited for a wide class of real-world problems, where constraints on the variables are often given. Different techniques for handling constraints in randomised algorithms have been proposed, see [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF] for a survey. For ES, common techniques are resampling, i.e. resample a solution till it lies in the feasible domain, repair of solutions that project unfeasible points onto the feasible domain (e.g. [START_REF] Arnold | Analysis of a repair mechanism for the (1, λ)-ES applied to a simple constrained problem[END_REF]), penalty methods where unfeasible solutions are penalised either by a quantity that depends on the distance to the constraint (e.g. [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF] with adaptive penalty weights) (if this latter one can be computed) or by the constraint value itself (e.g. stochastic ranking [START_REF] Runarsson | Stochastic ranking for constrained evolutionary optimization[END_REF]) or methods inspired from multiobjective optimization (e.g. [START_REF] Efrén | A simple multimembered evolution strategy to solve constrained optimization problems[END_REF]).

Alexandre Chotard, Anne Auger and Nikolaus Hansen work in TAO, at INRIA-Saclay and LRI in University Paris-Sud, France (mail at name@lri.fr).

Our thanks to Dirk Arnold for suggesting this work during PPSN is Sicily.

In this paper we focus on the resampling method and study it on a simple constraint problem. More precisely, we study a (1, λ)-ES optimizing a linear function with a linear constraint and resampling any unfeasible solution until a feasible solution is sampled. The linear function models the situation where the current point is, relatively to the step-size, far from the optimum and "solving" this function means diverging. The linear constraint models being close to the constraint relatively to the step-size and far from other constraints. Due to the invariance of the algorithm to the composition of the objective function by a strictly increasing map, the linear function could be composed by a function without derivative and with many discontinuities without any impact on our analysis.

The problem we address was studied previously for different step-size adaptation mechanisms: with constant step-size, self-adaptation and cumulative step-size adaptation [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF], [START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF]. The drawn conclusion is that when adapting the step-size the (1, λ)-ES fails to diverge unless some requirements on internal parameters of the algorithm are met. However, the approach followed in the aforementioned studies relies on finding simplified theoretical models to explain the behaviour of the algorithm: typically those models arise by doing some approximations (considering some random variables equal to their expected value, ...) and assuming some mathematical properties like the existence of stationary distributions of underlying Markov chains.

In contrast, our motivation is to study the real-in the sense not simplified-algorithm and prove rigorously different mathematical properties of the algorithm allowing to deduce the exact behaviour of the algorithm, as well as to provide tools and methodology for such studies. Our theoretical studies need to be complemented by simulations of the convergence/divergence rates. The mathematical properties that we derive show that these numerical simulations converge fast.

As for the step-size adaptation mechanism, our aim is to study the cumulative step-size adaptation (CSA), default step-size mechanism for the CMA-ES algorithm [START_REF] Hansen | Completely derandomized selfadaptation in evolution strategies[END_REF]. The mathematical object to study for this purpose is a discrete time, continuous state space Markov chain that is defined as the couple: evolution path and normalized distance to the constraint. More precisely stability properties like irreducibility, existence of a stationary distribution of this Markov chain need to be studied to deduce the geometric divergence of the CSA and have a rigorous mathematical framework to perform Monte Carlo simulations allowing to study the influence of different parameters of the algorithm. We start however by illustrating in details the methodology on the simpler case where the step-size is constant. We deduce in this case the divergence at a constant speed. We keep-due to some space limitation-the details of the generalization to the CSA study for a future publication and give only a sketch of the results.

This paper is organized as follows. In Section II we define the (1, λ)-ES using resampling and the problem. In Section III we provide some preliminary derivations on the distributions that come into play for the analysis. In Section IV we analyze the constant step-size case: exhibit the Markov chain, prove its stability and deduce the divergence of the (1, λ)-ES on the constraint problem. In Section V we sketch out our results when the step-size is adapted using cumulative step-size adaptation. Finally we discuss our results and our methodology in Section VI.

Notations

Throughout this article, we denote by ϕ the density function of the standard multivariate normal distribution, and Φ the cumulative distribution function of a standard univariate normal distribution. The standard (unidimensional) normal distribution is denoted N (0, 1), the (n-dimensional) multivariate normal distribution with covariance matrix identity is denoted N (0, Id n ) and the i th order statistic of λ i.i.d. standard normal random variables is denoted N i:λ . The uniform distribution on an interval I is denoted U I . We denote µ Leb the Lebesgue measure. The set of natural numbers (including 0) is denoted N, and the set of real numbers R. We denote R + the set {x ∈ R|x ≥ 0}, and for A ⊂ R n , the set A * denotes A\{0} and 1 A denotes the indicator function of A. For two vectors x ∈ R n and y ∈ R n , we denote [x] i the i th -coordinate of x, and x.y the scalar product of x and y. Take (a, b) ∈ N 2 with a ≥ b, we denote [a..b] the interval of integers between a and b. For a topological set X , B(X ) denotes the Borel algebra of X . For X and Y two random vectors, we denote X d = Y if X and Y are equal in distribution. For (X t ) t∈N a sequence of random variables and X a random variable we denote X t a.s. → X if X t converges almost surely to X and X t P → X if X t converges in probability to X.

II. PROBLEM STATEMENT AND ALGORITHM DEFINITION

A. (1, λ)-ES with resampling

In this paper, we study the behaviour of a (1, λ)-Evolution Strategy maximizing a function f : R n → R, λ ≥ 2, n ≥ 2, with a constraint defined by a function g : R n → R restricting the feasible space to X feasible = {x ∈ R n |g(x) ≥ 0}. To handle the constraint, the algorithm resamples any unfeasible solution until a feasible solution is found.

From iteration t ∈ N, given the vector X t ∈ R n and stepsize σ t ∈ R * + , the algorithm generates λ new candidates: until it lies within the feasible domain. The first feasible i th candidate solution is denoted Ỹi t and the realization of the multivariate normal distribution giving Ỹi t is Ñi t , which is called a feasible step. Note that Ñi t is not distributed as a multivariate normal distribution, further details on its distribution are given later on.

Y i t = X t + σ t N i t , (1)
We define

⋆ = argmax i∈[1..λ]
f ( Ỹi t ) as the index realizing the maximum objective function, and call Ñ⋆ t the selected step. The vector X t is then updated as the solution realizing the maximum value of the objective function, i.e.

X t+1 = Ỹ⋆ t = X t + σ t Ñ⋆ t . (2)
The step-size and other internal parameters are then adapted. We denote for the moment in a non specific manner the adaptation as

σ t+1 = σ t ξ t (3)
where ξ t is a random variable whose distribution is a function of the selected steps ( Ñ⋆ i ) i≤t . We will define later on specific rules for this adaptation.

B. Linear fitness function with linear constraint

In this paper, we consider the case where f , the function that we optimize, and g, the constraint, are linear functions. W.l.o.g., we assume that ∇f = ∇g = 1. We denote n := -∇g a vector normal to the constraint hyperplane. We choose an orthonormal Euclidean coordinate system with basis (e i ) i∈ [1..n] with its origin located on the constraint hyperplane where e 1 is equal to the gradient ∇f , hence

f (x) = [x] 1 (4)
and the vector e 2 lives in the plane generated by ∇f and n and is such that the angle between e 2 and n is positive. We define θ the angle between ∇f and n, and restrict our study to θ ∈ (0, π/2). The function g can be seen as a signed distance to the linear constraint as

g(x) = x.∇g = -x.n = -[x] 1 cos θ -[x] 2 sin θ . ( 5
)
A point is feasible if and only if g(x) ≥ 0 (see Figure 1).

Overall the problem reads

maximize f (x) = [x] 1 subject to g(x) = -[x] 1 cos θ -[x] 2 sin θ ≥ 0 . (6)
Although Ñi t and Ñ⋆ t are in R n , due to the choice of the coordinate system and the independence of the sequence

([N i t ] k ) k∈[1.
.n] , only the two first coordinates of these vectors are affected by the resampling implied by g and the selection according to f . Therefore

[ Ñ⋆ t ] k ∼ N (0, 1) for k ∈ [3..n].
With an abuse of notations, the vector Ñi t will denote the 2-dimensional vector

([ Ñi t ] 1 , [ Ñi t ] 2 ), likewise Ñ⋆ t will also denote the 2-dimensional vector ([ Ñ⋆ t ] 1 , [ Ñ⋆ t ] 2 )
, and n will denote the 2-dimensional vector (cos θ, sin θ). The coordinate system will also be used as (e 1 , e 2 ) only.

Following [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF], [START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF], [START_REF] Arnold | On the behaviour of the (1 + 1)-ES for a simple constrained problem[END_REF], we denote the normalized signed distance to the constraint as δ t , that is

δ t = g(X t ) σ t . (7)
We initialize the algorithm by choosing X 0 = -n and σ 0 = 1, which implies that δ 0 = 1.

III. PRELIMINARY RESULTS AND DEFINITIONS

Throughout this section we derive the probability density functions of the random vectors Ñi t and Ñ⋆ t and give a definition of Ñi t and of Ñ⋆ t as a function of δ t and of an i.i.d. sequence of random vectors.

A. Feasible steps

The random vector Ñi t , the i th feasible step, is distributed as the standard multivariate normal distribution truncated by the constraint, as stated in the following lemma.

Lemma 1: Let a (1, λ)-ES with resampling optimize a function f under a constraint function g. If g is a linear form determined by a vector n as in [START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF], then the distribution of the feasible step Ñi t only depends on the normalized distance to the constraint δ t and its density given that δ t equals δ reads

p δ (x) = ϕ(x)1 R+ (δ -x.n) Φ(δ) . ( 8
)
Proof: A solution Y i t is feasible if and only if g(Y i t ) ≥ 0, which is equivalent to -(X t + σ t N i t )
.n ≥ 0. Hence dividing by σ t , a solution is feasible if and only if δ t = -X t .n/σ t ≥ N i t .n. Since a standard multivariate normal distribution is rotational invariant, N i t .n follows a standard (unidimensional) normal distribution. Hence the probability that a solution Y i t or a step N i t is feasible is given by Pr(N (0, 1)

≤ δ t ) = Φ (δ t ) .

Therefore the density probability function of the random variable

Ñi t .n for δ t = δ is x → ϕ(x)1 R+ (δ -x)/Φ(δ).
For any vector n ⊥ orthogonal to n the random variable Ñi t .n ⊥ was not affected by the resampling and is therefore still distributed as a standard (unidimensional) normal distribution. With a change of variables using the fact that the standard multivariate normal distribution is rotational invariant we obtain the joint distribution of Eq. [START_REF] Hansen | Completely derandomized selfadaptation in evolution strategies[END_REF].

Then the marginal density function p 1,δ of [ Ñi t ] 1 can be computed by integrating Eq. ( 8) over [x] 2 and reads

p 1,δ (x) = ϕ (x) Φ δ-x cos θ sin θ Φ (δ) , (9)
(see [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF]Eq. 4] for details) and we denote F 1,δ its cumulative distribution function.

It will be important in the sequel to be able to express the vector Ñi t as a function of δ t and of a finite number of random samples. Hence we give an alternative way to sample Ñi t rather than the resampling technique that involves an unbounded number of samples.

Lemma 2: Let a (1, λ)-ES with resampling optimize a function f under a constraint function g, where g is a linear form determined by a vector n as in [START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF]. Let the feasible step Ñi t be the random vector described in Lemma 1 and Q be the 2-dimensional rotation matrix of angle θ. Then

Ñi t d = F -1 δt (U i t )n + N i t n ⊥ = Q -1 F -1 δt (U i t ) N i t ( 10
)
where

F -1 δt denotes the generalized inverse of the cumu- lative distribution of Ñi t .n 1 , U i t ∼ U [0,1] , N i t ∼ N (0, 1) with (U i t ) i∈[1..λ],t∈N i.i.d. and (N i t ) i∈[1..λ],t∈N i.i.d. random variables.
Proof: We define a new coordinate system (n, n ⊥ ) (see Figure 1). It is the image of (e 1 , e 2 ) by Q. In the new basis (n, n ⊥ ), only the coordinate along n is affected by the resampling. Hence the random variable Ñi t .n follows a truncated normal distribution with cumulative distribution function Fδt equal to min(1, Φ(x)/Φ(δ t )), while the random variable Ñi t .n ⊥ follows an independent standard normal distribution, hence Ñi

t d = ( Ñi t .n)n + N i t n ⊥ .
Using the fact that if a random variable has a cumulative distribution F , then for F -1 the generalized inverse of F , F -1 (U ) with U ∼ U [0,1] has the same distribution as this random variable, we get that

F -1 δt (U i t ) d = Ñi t .
n, so we obtain Eq. [START_REF] Efrén | A simple multimembered evolution strategy to solve constrained optimization problems[END_REF]. We now extend our study to the selected step Ñ⋆ t .

B. Selected step

The selected step Ñ⋆ t is chosen among the different feasible steps ( Ñi t ) i∈ [1..λ] to maximize the function f , and has the density described in the following lemma.

Lemma 3: Let a (1, λ)-ES with resampling optimize the problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF]. Then the distribution of the selected step Ñ⋆ t only depends on the normalized distance to the constraint δ t and its density given that δ t equals δ reads

p ⋆ δ (x) = λp δ (x) F 1,δ ([x] 1 ) λ-1 , (11)
= λ ϕ(x)1 R+ (δ -x.n) Φ(δ) [x]1 -∞ ϕ(u) Φ( δ-u cos θ sin θ ) Φ(δ) du λ-1
where p δ is the density of Ñi t given that δ t = δ given in Eq. ( 8) and F 1,δ the cumulative distribution function of [ Ñi t ] 1 whose density is given in Eq. ( 9) and n the vector (cos θ, sin θ).

Proof:

The function f being linear, the rankings on ( Ñi t ) i∈ [1..λ] corresponds to the order statistic on

([ Ñi t ] 1 ) i∈[1..λ] . If we look at the joint cumulative distribution F ⋆ δ of Ñ⋆ t F ⋆ δ (x, y) = Pr [ Ñ⋆ t ] 1 ≤ x, [ Ñ⋆ t ] 2 ≤ y = λ i=1 Pr Ñi t ≤ x y , [ Ñj t ] 1 < [ Ñi t ] 1 for j = i
by summing disjoints events. The vectors ( Ñi t ) i∈[1.

.λ] being independent and identically distributed

F ⋆ δ (x, y) = λ Pr Ñ1 t ≤ x y , [ Ñj t ] 1 < [ Ñ1 t ] 1 for j = 1 = λ x -∞ y -∞ p δ (u, v) λ j=2 Pr([ Ñj t ] 1 < u)dvdu = λ x -∞ y -∞ p δ (u, v)F 1,δ (u) λ-1 dvdu .
Deriving F ⋆ δ on x and y yields the density of Ñ⋆ t of Eq. ( 11).

We may now obtain the marginal of [ Ñ⋆ t ] 1 and [ Ñ⋆ t ] 2 . Corollary 1: Let a (1, λ)-ES with resampling optimize the problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF]. Then the marginal distribution of [ Ñ⋆ t ] 1 only depends of δ t and its density given that δ t equals δ reads

p ⋆ 1,δ (x) = λp 1,δ (x)F 1,δ (x) λ-1 , (12)
= λϕ(x) Φ δ-x cos θ sin θ Φ(δ) F 1,δ (x) λ-1 ,
and the same holds for [ Ñ⋆ t ] 2 whose marginal density reads

p ⋆ 2,δ (y) = λ ϕ(y) Φ(δ) δ-y sin θ cos θ -∞ ϕ(u)F 1,δ (u) λ-1 du . ( 13
)
Proof: Integrating Eq. ( 11) directly yields Eq. ( 12). The conditional density function of

[ Ñ⋆ t ] 2 is p ⋆ 2,δ (y|[ Ñ⋆ t ] 1 = x) = p ⋆ δ ((x, y)) p ⋆ 1,δ (x)
.

As p ⋆ 2,δ (y) = R p ⋆ 2,δ (y|[ Ñ⋆ t ] 1 = x)p ⋆ 1,δ (x
)dx, using the previous equation with Eq. [START_REF] Runarsson | Stochastic ranking for constrained evolutionary optimization[END_REF] gives that p ⋆ 2,δ (y) = R λp δ ((x, y))F 1,δ (x) λ-1 dx, which with Eq. ( 8) gives

p ⋆ 2,δ (y) = λ ϕ(y) Φ(δ) R ϕ(x)1 R+ δ - x y .n F 1,δ (x) λ-1 dx.
The condition δx cos θy sin θ ≥ 0 is equivalent to x ≤ (δy sin θ)/ cos θ, hence Eq. ( 13) holds. We will need in the next sections an expression of the random vector Ñ⋆ t as a function of δ t and a random vector composed of a finite number of i.i.d. random variables. To do so, using notations of Lemma 2, we define the function

G : R + × ([0, 1] × R) → R 2 as G(δ, w) = Q -1 F -1 δ ([w] 1 ) [w] 2 . ( 14
)
According to Lemma 2, given that U ∼ U [0,1] and N ∼ N (0, 1), ( F -1 δ (U ), N ) (resp. G(δ, (U, N ))) is distributed as the resampled step Ñi t in the coordinate system (n, n ⊥ )

(resp. (e 1 , e 2 )). Finally, let (w i ) i∈[1..λ] ∈ ([0, 1] × R) λ and let G : R + × ([0, 1] × R) λ → R 2 be the function defined as G(δ, (w i ) i∈[1..λ] ) = argmax N∈{ G(δ,wi)|i∈[1..λ]} f (N) . (15)
As shown in the following proposition, given that

W i t ∼ (U [0,1] , N (0, 1)) and W t = (W i t ) i∈[1.
.λ] , the function G(δ, W t ) is distributed as the selected step Ñ⋆ t . Proposition 1: Let a (1, λ)-ES with resampling optimize the problem defined in Eq. ( 6), and let (W i t ) i∈[1..λ],t∈N be an i.i.d. sequence of random vectors with W i t ∼ (U [0,1] , N (0, 1)), and

W t = (W i t ) i∈[1..λ] . Then Ñ⋆ t d = G(δ t , W t ) , (16)
where the function G is defined in Eq. ( 15). Proof:

Since f is a linear function f ( Ỹi t ) = f (X t ) + σ t f ( Ñi t ), so f ( Ỹi t ) ≤ f ( Ỹj t ) is equivalent to f ( Ñi t ) ≤ f ( Ñj t ). Hence ⋆ = argmax i∈[1..λ] f ( Ñi t ) and therefore Ñ⋆ t = argmax N∈{ Ñi t |i∈[1..λ]} f (N).
From Lemma 2 and Eq. ( 14), Ñi

t d = G(δ t , W i t ), so Ñ⋆ t d = argmax N∈{ G(δt,W i t )|i∈[1..λ]} f (N), which from (15) is G(δ t , W t ).

IV. CONSTANT STEP-SIZE CASE

We illustrate in this section our methodology analysis on the simple case where the step-size is constantly equal to σ and prove that then (X t ) t∈N diverges almost surely at constant speed (Theorem 1). The analysis of the CSA will then be a generalisation of the results presented here, with a few more technical results to derive.

As suggested in [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF], the sequence (δ t ) t∈N plays a central role for the analysis, and we will show that it admits a stationary measure. We first prove that this sequence is an homogeneous Markov chain.

Proposition 2: Consider the (1, λ)-ES with resampling and with constant step-size σ optimizing the constraint problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF]. Then the sequence δ t = g(X t )/σ is an homogeneous Markov chain on R + and

δ t+1 = δ t -Ñ⋆ t .n d = δ t -G(δ t , W t ).n , ( 17
)
where G is the function defined in (15) and The Markov Chain (δ t ) t∈N comes into play for investigating the divergence of f (X t ) = [X t ] 1 . Indeed, we can express

(W t ) t∈N = (W i t ) i∈[1..λ],t∈N is an i.i.d. sequence with W i t ∼ (U [0,1] , N (0, 1)) for all (i, t) ∈ [1..λ] × N
[Xt-X0]1 t
in the following manner:

[X t -X 0 ] 1 t = 1 t t-1 k=0 [X k+1 ] 1 -[X k ] 1 = σ t t-1 k=0 [ Ñ⋆ k ] 1 d = σ t t-1 k=0 [G(δ k , W k )] 1 . (18)
The latter term suggests the use of a Law of Large Numbers (LLN) to prove the convergence of [Xt-X0]1 t which will in turn imply-if the limit is positive-the divergence of f (X t ) at a constant rate. Sufficient conditions on a Markov chain to be able to apply the LLN include the existence of an invariant probability measure π. The limit term is then expressed as an expectation over the stationary distribution. More precisely, assume the LLN can be applied, the following limit will hold

lim t→∞ [X t -X 0 ] 1 t = σ R + E ([G(δ, W)] 1 ) π(dδ) (19) = lim t→∞ E δ0∼µ ([X t+1 ] 1 -[X t ] 1 ) , (20)
with µ any initial distribution. The latter term corresponds to the limit of the progress rate (see [2, Eq. 2]). The invariant measure π is also underlying the study carried out in [2, Section 4] where more precisely it is stated: "Assuming for now that the mutation strength σ is held constant, when the algorithm is iterated, the distribution of δ-values tends to a stationary limit distribution.". We will now provide a formal proof that indeed (δ t ) t∈N admits a stationary limit distribution π, as well as prove some other useful properties that will allow us in the end to conclude to the divergence of (f (X t )) t∈N .

A. Study of the stability of (δ t ) t∈N

We study in this section the stability of (δ t ) t∈N . We first derive its transition kernel P (δ, A) := Pr(δ t+1 ∈ A|δ t = δ) for all δ ∈ R + and A ∈ B(R + ). Since Pr(δ t+1 ∈ A|δ t = δ) = Pr(δ t -Ñ⋆ t .n ∈ A|δ t = δ) ,

P (δ, A) = R 2 1 A (δ -u.n) p ⋆ δ (u) du (21)
where p ⋆ δ is the density of Ñ⋆ t given in [START_REF] Runarsson | Stochastic ranking for constrained evolutionary optimization[END_REF]. For t ∈ N * , the t-step transition kernel P t is defined by P t (δ, A) := Pr(δ t ∈ A|δ 0 = δ).

From the transition kernel, we will now derive the first properties on the Markov chain (δ t ) t∈N . First of all we investigate the so-called ψ-irreducible property.

A Markov chain (δ t ) t∈N on a state space R + is ψirreducible if there exists a non-trivial measure ψ such that for all set A ∈ B(R + ) with ψ(A) > 0 and for all δ ∈ R + , there exists t ∈ N * such that P t (δ, A) > 0. We denote B + (R + ) the set of Borel sets of R + with strictly positive ψ-measure.

We also need the notion of small sets: a set C ∈ B(R + ) is called a small set if there exists m ∈ N * and a non trivial measure ν m such that for all set A ∈ B(R + ) and all δ ∈ C P m (δ, A) ≥ ν m (A) .

(

) 22
If there exists C a ν 1 -small set such that ν 1 (C) > 0 then the Markov chain is said strongly aperiodic. Proposition 3: Consider a (1, λ)-ES with resampling and with constant step-size optimizing the constraint problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF] and let (δ t ) t∈N be the Markov chain exhibited in (17). Then (δ t ) t∈N is µ Leb -irreducible, strongly aperiodic, and compact sets are small sets.

Proof: Using Eq. ( 21) and Eq. ( 11) the transition kernel can be written

P (δ, A) = λ R 2 1 A (δ- x y .n) ϕ(x)ϕ(y) Φ(δ) F 1,δ (x) λ-1 dydx .
We remove δ from the indicator function by a substitution of variables u = δx cos θy sin θ, and v = x sin θy cos θ. As this substitution is the composition of a rotation and a translation the determinant of its Jacobian matrix is 1. We denote

h δ : (u, v) → (δ -u) cos θ + v sin θ, h ⊥ δ : (u, v) → (δ -u) sin θ -v cos θ and g(δ, u, v) → λϕ(h δ (u, v))ϕ(h ⊥ δ (u, v))/Φ(δ)F 1,δ (h δ (u, v)) λ-1 . Then x = h δ (u, v), y = h ⊥ δ (u, v) and P (δ, A) = R R 1 A (u)g(δ, u, v)dvdu . ( 23
)
For all δ, u, v the function g(δ, u, v) is strictly positive hence for all A with µ Leb (A) > 0, P (δ, A) > 0. Hence (δ t ) t∈N is irreducible with respect to the Lebesgue measure.

In addition, the function (δ, u, v) → g(δ, u, v) is continuous as the composition of continuous functions (the continuity of δ → F 1,δ (x) for all x coming from the dominated convergence theorem). Given a compact C we hence know that there exists g C > 0 such that for all

(δ, u, v) ∈ C × [0, 1] 2 , g(δ, u, v) ≥ g C > 0. Hence for all δ ∈ C, P (δ, A) ≥ g C µ Leb (A ∩ [0, 1]) :=νC (A)
.

The measure ν C being non-trivial, the previous equation shows that compact sets are small and that for C a compact such that µ Leb (C ∩ [0, 1]) > 0, we have ν C (C) > 0 hence the chain is strongly aperiodic. The application of the LLN for a ψ-irreducible Markov chain (δ t ) t∈N on a state space R + requires the existence of an invariant measure π, that is satisfying for all A ∈ B(R + )

π(A) = R+ P (δ, A)π(dδ) . (24)
If a Markov chain admits an invariant probability measure then the Markov chain is called positive.

A typical assumption to apply the LLN is positivity and Harris-recurrence. A ψ-irreducible chain (δ t ) t∈N on a state space R + is Harris-recurrent if for all set A ∈ B + (R + ) and for all δ ∈ R + , Pr(η A = ∞|δ 0 = δ) = 1 where η A is the occupation time of A, i.e. η A = ∞ t=1 1 A (δ t ). We will show that the Markov chain (δ t ) t∈N is positive and Harris-recurrent by using so-called Foster-Lyapunov drift conditions: define the drift operator for a positive function V as

∆V (δ) = E[V (δ t+1 )|δ t = δ] -V (δ) .
Drift conditions translate that outside a small set, the drift operator is negative. We will show a drift condition for Vgeometric ergodicity where given a function f ≥ 1, a positive and Harris-recurrent chain (δ t ) t∈N with invariant measure π is called f -geometrically ergodic if π(f ) < ∞ and there exists r f > 1 such that

t∈N r t f P t (δ, •) -π f < ∞ , ∀δ ∈ R + , (25)
where for ν a signed measure ν f denotes sup g:|g|≤f | R+ g(x)ν(dx)|.

To prove V -geometric ergodicity, we will prove that there exists a small set C, constants b ∈ R, ǫ ∈ R * + and a function V ≥ 1 finite for at least some δ 0 ∈ R + such that for all

δ ∈ R + ∆V (δ) ≤ -ǫV (δ) + b1 C (δ) . (26)
If the Markov chain (δ t ) t∈N is ψ-irreducible and aperiodic, this drift condition implies that the chain is V -geometrically ergodic [9, Theorem 15.0.1]2 as well as positive and Harrisrecurrent 3 .

Because compacts are small sets and drift conditions investigate the negativity outside a small set, we need to study the chain for δ large. The following lemma is a technical lemma studying the limit of E(exp(G(δ, W).n)) for δ to infinity.

Lemma 4: Consider the (1, λ)-ES with resampling optimizing the constraint problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF], and let G be the function defined in (15). We denote K and K the ran-

dom variables exp(G(δ, W).(a, b)) and exp(a|[G(δ, W)] 1 | + b|[G(δ, W)] 2 |). For W ∼ (U [0,1]
, N (0, 1)) λ and any (a, b) ∈ R 2 lim δ→+∞ E(K) = E(exp(aN λ:λ ))E(exp(bN (0, 1))) < ∞ and lim δ→+∞ E( K) < ∞ For the proof see the appendix. We are now ready to prove a drift condition for geometric ergodicity.

Proposition 4: Consider a (1, λ)-ES with resampling and with constant step-size optimizing the constraint problem (6) and let (δ t ) t∈N be the Markov chain exhibited in (17). The Markov chain (δ t ) t∈N is V -geometrically ergodic with V : δ → exp(αδ) for α > 0 small enough, and is Harrisrecurrent and positive with invariant probability measure π.

Proof: Take the function

V : δ → exp(αδ) then ∆V (δ) = E (exp (α (δ -G(δ, W).n))) -exp (αδ), ∆V V (δ) = E (exp (-αG(δ, W).n)) -1. With Lemma 4 we obtain lim δ→+∞ E (exp (-αG(δ, W).n)) = E (exp(-αN λ:λ cos θ)) E(exp(-αN (0, 1) sin θ)) < ∞ .
As the right hand side of the previous equation is finite we can invert integral with series with Fubini's theorem, so with Taylor series the limit equals to

i∈N (-α cos θ) i E N i λ:λ i! i∈N (-α sin θ) i E N (0, 1) i i! ,
which in turns yields

lim δ→+∞ ∆V V (δ) = (1 -αE(N λ:λ ) cos θ + o(α)) (1 +o(α))-1 = -αE(N λ:λ ) cos θ + o(α) .
Since for λ ≥ 2, E(N λ:λ ) > 0, for α > 0 and small enough we get lim δ→+∞ ∆V V (δ) < -ǫ < 0. Hence there exists ǫ > 0, M > 0 and b ∈ R such that ∆V (δ) ≤ -ǫV (δ) + b1 [0,M] (δ) .

According to Proposition 3, [0, M ] is a small set, hence it is petite [START_REF] Meyn | Markov chains and stochastic stability[END_REF]Proposition 5.5.3]. Furthermore (δ t ) t∈N is a ψ-irreducible aperiodic Markov chain so (δ t ) t∈N satisfies the conditions of Theorem 15.0.1 from [START_REF] Meyn | Markov chains and stochastic stability[END_REF], which with Lemma 15.2.2, Theorem 9.1.8 and Theorem 14.0.1 of [START_REF] Meyn | Markov chains and stochastic stability[END_REF] proves the proposition.

We now proved rigorously the existence (and unicity) of an invariant measure π for the Markov chain (δ t ) t∈N , which provides the so-called steady state behaviour in [2, Section 4]. As the Markov chain (δ t ) t∈N is positive and Harrisrecurrent we may now apply a Law of Large Numbers [START_REF] Meyn | Markov chains and stochastic stability[END_REF]Theorem 17.1.7] in Eq (18) to obtain the divergence of f (X t ) and an exact expression of the divergence rate.

Theorem 1: Consider a (1, λ)-ES with resampling and with constant step-size optimizing the constraint problem [START_REF] Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF] and let (δ t ) t∈N be the Markov chain exhibited in (17). The sequence ([X t ] 1 ) t∈N diverges in probability to +∞ at constant speed, that is

[X t -X 0 ] 1 t P -→ t→+∞ σE π×µW ([G (δ, W)] 1 ) > 0 , (27)
with G defined in (15) and W = (W i ) i∈ [1..λ] where

(W i ) i∈[1.
.λ] is an i.i.d. sequence such that W i ∼ (U [0,1] , N (0, 1)) and µ W is the probability measure of W.

Proof: From Proposition 4 the Markov chain (δ t ) t∈N is Harris-recurrent and positive, and since (W t ) t∈N is i.i.d., the chain (δ t , W t ) is also Harris-recurrent and positive with invariant probability measure π × µ W , so to apply the Law of Large Numbers [START_REF] Meyn | Markov chains and stochastic stability[END_REF]Theorem 17

.0.1] to [G] 1 we only need [G] 1 to be π × µ W -integrable. With Fubini-Tonelli's theorem E π×µW (|[G(δ, W)] 1 |) equals to E π (E µW (|[G(δ, W)] 1 |)). As δ ≥ 0, we have Φ(δ) ≥ Φ(0) = 1/2,
and for all x ∈ R as Φ(x) ≤ 1, F 1,δ (x) ≤ 1 and ϕ(x) ≤ exp(-x 2 /2) with Eq. ( 12) we obtain that

|x|p ⋆ 1,δ (x) ≤ 2λ|x| exp(-x 2 /2) so the function x → |x|p ⋆ 1,δ (x) is integrable. Hence for all δ ∈ R + , E µW (|[G(δ, W)] 1 |) is finite. Using the dominated convergence theorem, the function δ → F 1,δ (x) is continuous, hence so is δ → p ⋆ 1,δ (x). From (12) |x|p ⋆ 1,δ (x) ≤ 2λ|x|ϕ(x)
, which is integrable, so the dominated convergence theorem implies that the function δ → E µW (|[G(δ, W] 1 |) is continuous. Finally, using Lemma 4 with Jensen's inequality shows that

lim δ→+∞ E µW (|[G(δ, W)] 1 |) is finite. Therefore the function δ → E µW (|[G(δ, W] 1 |) is bounded by a constant M ∈ R + . As π is a probability measure E π (E µW (|[G(δ, W)] 1 |)) ≤ M < ∞, meaning [G] 1 is π × µ W -integrable.
Hence we may apply the LLN on Eq. ( 18)

σ t t-1 k=0 [G(δ k , W k )] 1 a.s. -→ t→+∞ σE π×µW ([G(δ, W)] 1 ) < ∞ .
The equality in distribution in (18) allows us to deduce the convergence in probability of the left hand side of (18) to the right hand side of the previous equation.

As the measure π is an invariant measure for the Markov chain (δ t ) t∈N , using (17), E π×µW (δ) = E π×µW (δ -G(δ, W).n), hence E π×µW (G(δ, W).n) = 0 and thus

E π×µW ([G(δ, W)] 1 ) = -tan θE π×µW ([G(δ, W)] 2 ) .
We see from Eq. ( 13) that for y > 0, p ⋆ 2,δ (y) < p ⋆ 2,δ (-y) hence the expected value E π×µW ([G(δ, W)] 2 ) is strictly negative. With the previous equation it implies that E π×µW ([G(δ, W)] 1 ) is strictly positive.

We showed rigorously the divergence of [X t ] 1 and gave an exact expression of the divergence rate, which is the limit of the progress rate defined in [2, Eq. ( 2)]. The fact that the chain (δ t ) t∈N is V -geometrically ergodic gives that t r t V P t (δ, •)π V < ∞. This implies that the distribution π can be simulated efficiently by a Monte Carlo simulation allowing to have precise estimations of the divergence rate of [X t ] 1 . Assuming a CLT could be applied, confidence intervals on the Monte Carlo simulations could also be obtained.

A Monte Carlo simulation of the right hand side of Eq. (27) for 10 6 time steps gives the progress rate ϕ ⋆ = E([X t+1 -X t ] 1 ), which once normalized by σ and λ yields 2. We normalize per λ as in evolution strategies the cost of the algorithm is assumed to be the number of f -calls. We see that for small values of θ, the normalized serial progress rate assumes roughly ϕ ⋆ /λ ≈ θ 2 . Only for larger constraint angles the serial progress rate depends on λ where smaller λ are preferable. Fig. 3 is obtained through simulations of the Markov chain (δ t ) t∈N defined in Eq. (17) for 10 6 time steps where the values of (δ t ) t∈N are averaged over time. We see that when θ → π/2 then E π (δ t ) → +∞ since the selection does not attract X t towards the constraint anymore, while the resampling still repels X t from the constraint. With a larger population size the algorithm is closer to the constraint, as better samples are more likely to be found close to the constraint.

V. CUMULATIVE STEP-SIZE ADAPTATION CASE

We generalise the previous results to the cumulative stepsize adaptation mechanism. However due to space limitation we only sketch the results that we plan to present in details in an extended version of the paper. CSA introduces a new variable, p t , called the evolution path. It is a weighted recombination of the previous selected steps, where the weight of Ñ⋆ k is proportional to (1c) t-1-k with c ∈ (0, 1]

E π ([N ⋆ t ] 1 )/λ λ =5 λ =10 λ =20 Fig. 2.
Normalized progress rate ϕ ⋆ = E([ Ñ⋆ t ] 1 ) divided by λ for the (1, λ)-ES with constant step-size and resampling, plotted against the constraint angle θ, for λ ∈ {5, 10, 20}. being the cumulation parameter. For c = 1 the algorithm has "no memory" and the evolution path p t is Ñ⋆ t-1 . The stepsize is adapted depending on the norm of p t [START_REF] Hansen | Completely derandomized selfadaptation in evolution strategies[END_REF]. The Markov chain to study in this case is (δ t , p t ) t∈N , except when c = 1 where it is (δ t ) t∈N .

As in Section IV if the Markov chain is ψ-irreducible, aperiodic, and compact sets are small, then for c = 1 the Markov chain (δ t ) t∈N is positive, Harris recurrent and V -geometrically ergodic, and a LLN can be applied on ln(σ t /σ 0 ) to obtain that

1 t ln σ t σ 0 a.s. -→ t→∞ E πc×µW G (δ, W) 2 -2 2d σ n , (28)
with π c the stationary measure of (δ t ) t∈N , G defined in (15), W = (W i ) i∈ [1..λ] where (W i ) i∈[1..λ] is an i.i.d. sequence such that W i ∼ (U [0,1] , N (0, 1)) and µ W the probability measure of W. So the step-size converges (resp. diverges) exponentially fast when the right hand side of Eq. ( 28) is strictly negative (resp. strictly positive).

VI. DISCUSSION

We investigated the (1, λ)-ES with constant step-size optimizing a linear function under a linear constraint handled by resampling unfeasible solutions. We prove the stability (formally V-geometric ergodicity) of the Markov chain (δ t ) t∈N defined as the normalised distance to the constraint, which was pressumed in [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF]. This property implies the divergence of the algorithm at a constant speed (see Theorem 1). In addition, it ensures (fast) convergence of Monte Carlo simulations of the divergence rate, justifying their use.

We believe that with the same approach, the CSA can be analysed. Simulations suggest that geometric divergence occurs for a small enough cumulation parameter, c, or large enough population size, λ. However, smaller values of the constraint angle seem to increase the difficulty of the problem arbitrarily, i.e. no given values for c and λ solve the problem for every θ ∈ (0, π/2).

Using a different covariance matrix to generate new samples can be interpreted as a change of the constraint angle. Therefore a correct adaptation of the covariance matrix will render the problem arbitrarily close to the one with θ = π/2. The unconstrained linear function case has been shown to be solved by a (1, λ)-ES with cumulative step-size adaptation for a population size larger than 3, regardless of other internal parameters [START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF]. We believe this is a strong argument for using covariance matrix adaptation with ES when dealing with constraints, as pure step-size adaptation has been shown to be liable to fail on even a very basic problem.

This work provides a methodology that can be applied to many ES variants. It demonstrates that a rigorous analysis of the constrained problem can be achieved. It relies on the theory of Markov chains for a continuous state space that once again proves to be a natural theoretical tool for analysing ESs, complementing particularly well previous studies [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF], [START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF], [START_REF] Arnold | On the behaviour of the (1 + 1)-ES for a simple constrained problem[END_REF].

Fig. 1 .

 1 Fig.1. Linear function with a linear constraint, in the plane generated by ∇f and n, a normal vector to the constraint hyperplane with angle θ ∈ (0, π/2) with ∇f . The point x is at distance g(x) from the constraint.

Fig. 3 .

 3 Fig. 3. Average normalized distance δ from the constraint for the (1, λ)-ES with constant step-size and resampling plotted against the constraint angle θ for λ ∈ {5, 10, 20}.

.

  Proof: It follows from the definition of δ t that δ t+1 = , W t ). Since δ t+1 has the same distribution as a time independent function of δ t and of W t where (W t ) t∈N are i.i.d., it is an homogeneous Markov chain.

	g(Xt+1) σt+1	=	-(Xt+σ Ñ⋆ t ).n σ	= δ t -Ñ⋆ t .n, and in Proposition 1
	we state that Ñ⋆ t	d = G(δ t

The generalized inverse of Fδ is F -1 δt (y) := inf x∈R { Fδt (x) ≥ y}.

The condition π(V ) < ∞ is given by [9, Theorem 14.0.1].

The function V of (26) is unbounded off petite sets [9, Lemma 15.2.2], hence with [9, Theorem 9.1.8] the Markov chain is Harris-recurrent.

ACKNOWLEDGMENTS

This work was supported by the grants ANR-2010-COSI-002 (SIMINOLE) and ANR-2012-MONU-0009 (NumBBO) of the French National Research Agency.

APPENDIX

Proof of Lemma 4.

Proof: From Proposition 1 the density probability function of G(δ, W) is p ⋆ δ , and from Eq. ( 11)

From Eq. ( 9) p 1,δ (x) = ϕ(x)Φ((δx cos θ)/ sin θ)/Φ(δ), so as δ ≥ 0 we have 1 ≥ Φ(δ) ≥ Φ(0) = 1/2, hence p 1,δ (x) ≤ 2ϕ(x). So p 1,δ (x) converges when δ → +∞ to ϕ(x) while being bounded by 2ϕ(x) which is integrable. Therefore we can apply Lebesgue's dominated convergence theorem: F 1,δ converges to Φ when δ → +∞ and is finite.

being dominated by 2λ exp(ax)ϕ(x)ϕ(y), which is integrable. Therefore by the dominated convergence theorem and as the density of N λ:λ is x → λϕ(x)Φ(x) λ-1 , when δ → +∞, R h δ,y (x)dx converges to ϕ(y)E(exp(aN λ:λ )) < ∞.

So the function y → exp(by) R h δ,y (x)dx converges to y → exp(by)ϕ(y)E(exp(aN λ:λ )) while being dominated by y → 2λϕ(y) exp(by) R exp(ax)ϕ(x)dx which is integrable. Therefore we may apply the dominated convergence theorem: E(exp(G(δ, W).(a, b))) converges to R exp(by)ϕ(y)E(exp(aN λ:λ ))dy which equals to E(exp(aN λ:λ ))E(exp(bN (0, 1))); and this quantity is finite.

The same reasoning gives that lim δ→∞ E( K) < ∞.


Alexandre Chotard
email: alexandre.chotard@lri.fr

Martin Holeňa
email: martin@cs.cas.cz

A Generalized Markov-Chain Modelling Approach to (1, λ)-ES Linear Optimization

Keywords: evolution strategies, continuous optimization, linear optimization, linear constraint, linear function, Markov chain models, Archimedean copulas

Several recent publications investigated Markov-chain modelling of linear optimization by a (1, λ)-ES, considering both unconstrained and linearly constrained optimization, and both constant and varying step size. All of them assume normality of the involved random steps, and while this is consistent with a black-box scenario, information on the function to be optimized (e.g. separability) may be exploited by the use of another distribution. The objective of our contribution is to complement previous studies realized with normal steps, and to give sufficient conditions on the distribution of the random steps for the success of a constant step-size (1, λ)-ES on the simple problem of a linear function with a linear constraint. The decomposition of a multidimensional distribution into its marginals and the copula combining them is applied to the new distributional assumptions, particular attention being paid to distributions with Archimedean copulas.

Introduction

Evolution Strategies (ES) are Derivative Free Optimization (DFO) methods, and as such are suited for the optimization of numerical problems in a black-box context, where the algorithm has no information on the function f it optimizes (e.g. existence of gradient) and can only query the function's values. In such a context, it is natural to assume normality of the random steps, as the normal distribution has maximum entropy for given mean and variance, meaning that it is the most general assumption one can make without the use of additional information on f . However such additional information may be available, and then using normal steps may not be optimal. Cases where different distributions have been studied include so-called Fast Evolution Strategies [START_REF] Yao | Fast evolution strategies[END_REF] or SNES [START_REF] Schaul | Benchmarking Separable Natural Evolution Strategies on the Noiseless and Noisy Black-box Optimization Testbeds[END_REF][START_REF] Schaul | High dimensions and heavy tails for natural evolution strategies[END_REF] which exploits the separability of f , or heavy-tail distributions on multimodal problems [START_REF] Hansen | When do heavy-tail distributions help?[END_REF][START_REF] Schaul | High dimensions and heavy tails for natural evolution strategies[END_REF].

In several recent publications [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF][START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF][START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF][START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF], attention has been paid to Markov-chain modelling of linear optimization by a (1, λ)-ES, i.e. by an evolution strategy in which λ children are generated from a single parent X ∈ R n by adding normally distributed n-dimensional random steps M , X ← X + σC 1 2 M , where M ∼ N (0, In).

(1)

Here, σ is called step size, C is a covariance matrix, and N (0, I n ) denotes the ndimensional standard normal distribution with zero mean and covariance matrix identity. The best among the λ children, i.e. the one with the highest fitness, becomes the parent of the next generation, and the step-size σ and the covariance matrix C may then be adapted to increase the probability of sampling better children. In this paper we relax the normality assumption of the movement M to a more general distribution H.

The linear function models a situation where the step-size is relatively small compared to the distance towards a local optimum. This is a simple problem that must be solved by any effective evolution strategy by diverging with positive increments of ∇f.M . This unconstrained case was studied in [START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF] for normal steps with cumulative step-size adaptation (the step-size adaptation mechanism in CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]).

Linear constraints naturally arise in real-world problems (e.g. need for positive values, box constraints) and also model a step-size relatively small compared to the curvature of the constraint. Many techniques to handle constraints in randomised algorithms have been proposed (see [START_REF] Coello Coello | Constraint-handling techniques used with evolutionary algorithms[END_REF]). In this paper we focus on the resampling method, which consists in resampling any unfeasible candidate until a feasible one is sampled. We chose this method as it makes the algorithm easier to study, and is consistent with the previous studies assuming normal steps [START_REF] Arnold | On the behaviour of the (1 + 1)-ES for a simple constrained problem[END_REF][START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF][START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF][START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF], studying constant step-size, self adaptation and cumulative step-size adaptation mechanisms (with fixed covariance matrix).

Our aim is to study the (1, λ)-ES with constant step-size, constant covariance matrix and random steps with a general absolutely continuous distribution H optimizing a linear function under a linear constraint handled through resampling. We want to extend the results obtained in [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF][START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF] using the theory of Markov chains. It is our hope that such results will help in designing new algorithms using information on the objective function to make non-normal steps. We pay a special attention to distributions with Archimedean copulas, which are a particularly well transparent alternative to the normal distribution. Such distributions have been recently considered in the Estimation of Distribution Algorithms [START_REF] Cuesta-Infante | Bivariate empirical and n-variate archimedean copulas in estimation of distribution algorithms[END_REF][START_REF] Wang | Copula estimation of distribution algorithms based on exchangeable archimedean copula[END_REF], continuing the trend of using copulas in that kind of evolutionary optimization algorithms [START_REF] Salinas-Gutierrez | Using copulas in estimation of distribution algorithms[END_REF].

In the next section, the basic setting for modelling the considered evolutionary optimization task is formally defined. In Section 3, the distributions of the feasible steps and of the selected steps are linked to the distribution of the random steps, and another way to sample them is provided. In Section 4, it is shown that, under some conditions on the distribution of the random steps, the normalized distance to the constraint defined in ( 5) is a ergodic Markov chain, and a law of large numbers for Markov chains is applied. Finally, Section 5 gives properties on the distribution of the random steps under which some of the aforementioned conditions are verified.

Due to a lack of space proofs were not included in this paper, and can instead be found at http://hal.inria.fr/docs/01/00/30/15/PDF/ppsn2014TRlinear constraintgeneraldistributions.pdf.

Notations

For (a, b) ∈ N 2 with a < b, [a..b] denotes the set of integers i such that a ≤ i ≤ b.

For X and Y two random vectors, X

= Y denotes that these variables are equal in distribution, X a.s.

→ Y and X P → Y denote, respectively, almost sure convergence and convergence in probability. For (x, y) ∈ R n , x.y denotes the scalar product between the vectors x and y, and for i ∈ [1..n], [x] i denotes the i th coordinate of x. For A a subset of R n , ½ A denotes the indicator function of A. For X a topological set, B(X ) denotes the Borel algebra on X .

Problem setting and algorithm definition

Throughout this paper, we study a (1, λ)-ES optimizing a linear function f : R n → R where λ ≥ 2 and n ≥ 2, with a linear constraint g : R n → R, handling the constraint by resampling unfeasible solutions until a feasible solution is sampled.

Take (e k ) k∈[1.

.n] a orthonormal basis of R n . We may assume ∇f to be normalized as the behaviour of an ES is invariant to the composition of the objective function by a strictly increasing function (e.g. h : x → x/ ∇f ), and the same holds for ∇g since our constraint handling method depends only on the inequality g(x) ≤ 0 which is invariant to the composition of g by a homothetic transformation. Hence w.l.o.g. we assume that ∇f = e 1 and ∇g = cos θe 1 + sin θe 2 with the set of feasible solutions X feasible := {x ∈ R n |g(x) ≤ 0}. We restrict our study to θ ∈ (0, π/2). Overall the problem reads

maximize f (x) = [x] 1 subject to g(x) = [x] 1 cos θ + [x] 2 sin θ ≤ 0 . (2)
At iteration t ∈ N, from a so-called parent point X t ∈ X feasible and with step-size σ t ∈ R * + we sample new candidate solutions by adding to X t a random vector σ t M i,j t where M i,j t is called a random step and (M i,j t ) i∈[1..λ],j∈N,t∈N is a i.i.d. sequence of random vectors with distribution H. The i index stands for the λ new samples to be generated, and the j index stands for the unbounded number of samples used by the resampling. We denote M i t a feasible step, that is the first element of (M i,j t ) j∈N such that X t + σ t M i t ∈ X feasible (random steps are sampled until a suitable candidate is found). The i th feasible solution

Y i t is then Y i t := X t + σ t M i t . (3)
Then we denote

⋆ := argmax i∈[1..λ] f (Y i t )
the index of the feasible solution maximizing the function f , and update the parent point where M ⋆ t is called the selected step. Then the step-size σ t , the distribution of the random steps H or other internal parameters may be adapted.

X t+1 := Y ⋆ t = X t + σ t M ⋆ t , (4)
Following [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF][START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF][START_REF] Arnold | On the behaviour of the (1 + 1)-ES for a simple constrained problem[END_REF][START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF] we define δ t as

δ t := - g(X t ) σ t . ( 5
)
3 Distribution of the feasible and selected steps

In this section we link the distributions of the random vectors M i t and M ⋆ t to the distribution of the random steps M i,j t , and give another way to sample M i t and M ⋆ t not requiring an unbounded number of samples.

Lemma 1. Let a (1, λ)-ES optimize the problem defined in (2) handling constraint through resampling. Take H the distribution of the random step M i,j t , and for δ ∈ R * + denote L δ := {x ∈ R n |g(x) ≤ δ}. Providing that H is absolutely continuous and that H(L δ ) > 0 for all δ ∈ R + , the distribution Hδ of the feasible step and H⋆ δ the distribution of the selected step when δ t = δ are absolutely continuous, and denoting h, hδ and h⋆ δ the probability density functions of, respectively, the random step, the feasible step M i t and the selected step

M ⋆ t when δ t = δ hδ (x) = h(x)½ L δ (x) H(L δ ) , (6)
and

h⋆ δ (x) = λ hδ (x) Hδ ((-∞, [x] 1 ) × R n-1 ) λ-1 = λ h(x)½ L δ (x)H((-∞, [x] 1 ) × R n-1 ∩ L δ ) λ-1 H(L δ ) λ . ( 7
)
The vectors (M i t ) i∈[1.

.λ] and M ⋆ t are functions of the vectors (M i,j t ) i∈[1..λ],j∈N and of δ t . In the following Lemma an equivalent way to sample M i t and M ⋆ t is given which uses a finite number of samples. This method is useful if one wants to avoid dealing with the infinite dimension space implied by the sequence (M i,j t ) i∈[1..λ,j∈N .

Lemma 2. Let a (1, λ)-ES optimize problem (2), handling the constraint through resampling, and take δ t as defined in [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF]. Let H denote the distribution of M i,j t that we assume absolutely continuous, ∇g ⊥ :=sin θe 1 + cos θe 2 , Q the rotation matrix of angle θ changing (e 1 , e 2 , . . . , e n ) into (∇g, ∇g ⊥ , . . . , e n ). Take

F 1,δ (x) := Pr(M i t .∇g ≤ x|δ t = δ), F 2,δ (x) := Pr(M i t .∇g ⊥ ≤ x|δ t = δ) and F k,δ (x) := Pr([M i t ] k ≤ x|δ t = δ) for k ∈ [3.
.n], the marginal cumulative distribution functions when δ t = δ, and C δ the copula of (M i t .∇g, M i t .∇g ⊥ , . . . , M i t .e n ). We define

G : (δ, (u i ) i∈[1..n] ) ∈ R + × [0, 1] n → Q    F -1
1,δ (u 1 ) . . .

F -1 n,δ (u n )    , (8)
G ⋆ : (δ, (v i ) i∈[1..λ] ) ∈ R + × [0, 1] nλ → argmax G∈{G(δ,vi)|i∈[1..λ]} f (G) . (9)
Then, if the copula C δ is constant in regard to δ, for

W t = (V i,t ) i∈[1..λ] a i.i.d. sequence with V i,t ∼ C δ G(δ t , V i,t ) (d) = M i t , (10)
G ⋆ (δ t , W t ) (d) = M ⋆ t . (11)
We may now use these results to show the divergence of the algorithm when the step-size is constant, using the theory of Markov chains [START_REF] Meyn | Markov chains and stochastic stability[END_REF].

Divergence of the (1, λ)-ES with constant step-size

Following the first part of [START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF], we restrict our attention to the constant step size in the remainder of the paper, that is for all t ∈ N we take σ t = σ ∈ R * + . From Eq. ( 4), by recurrence and dividing by t, we see that

[X t -X 0 ] 1 t = σ t t-1 i=0 M ⋆ i . (12)
The latter term suggests the use of a Law of Large Numbers to show the convergence of the left hand side to a constant that we call the divergence rate. The random vectors (M ⋆ t ) t∈N are not i.i.d. so in order to apply a Law of Large Numbers on the right hand side of the previous equation we use Markov chain theory, more precisely the fact that (M ⋆ t ) t∈N is a function of a (δ t , (M i,j t ) i∈[1..λ],j∈N ) t∈N which is a geometrically ergodic Markov chain. As (M i,j t ) i∈[1..λ],j∈N,t∈N is a i.i.d. sequence, it is a Markov chain, and the sequence (δ t ) t∈N is also a Markov chain as stated in the following proposition.

Proposition 1. Let a (1, λ)-ES with constant step-size optimize problem (2), handling the constraint through resampling, and take δ t as defined in [START_REF] Arnold | On the behaviour of the (1,λ)-ES for a simple constrained problem[END_REF]. Then no matter what distribution the i.i.d. sequence (M i,j t ) i∈[1..λ],(j,t)∈N 2 have, (δ t ) t∈N is a homogeneous Markov chain and

δ t+1 = δ t -g(M ⋆ t ) = δ t -cos θ[M ⋆ t ] 1 -sin θ[M ⋆ t ] 2 . ( 13
)
We now show ergodicity of the Markov chain (δ t ) t∈N , which implies that the t-steps transition kernel (the function A → Pr(δ t ∈ A|δ 0 = δ) for A ∈ B(R + )) converges towards a stationary measure π, generalizing Propositions 3 and 4 of [START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF].

Proposition 2. Let a (1, λ)-ES with constant step-size optimize problem (2), handling the constraint through resampling. We assume that the distribution of M i,j t is absolutely continuous with probability density function h, and that h is continuous and strictly positive on R n . Denote µ + the Lebesgue measure on (R + , B(R + )), and for α > 0 take the functions V : δ → δ, V α : δ → exp(αδ) and r 1 : δ → 1. Then (δ t ) t∈N is µ + -irreducible, aperiodic and compact sets are small sets for the Markov chain.

If the following two additional conditions are fulfilled

E(|g(M i,j t )| | δ t = δ) < ∞ for all δ ∈ R + , and (14)
lim δ→+∞ E(g(M ⋆ t )|δ t = δ) ∈ R * + , (15)
then (δ t ) t∈N is r 1 -ergodic and positive Harris recurrent with some invariant measure π. Furthermore, if

E(exp(g(M i,j t ))|δ t = δ) < ∞ for all δ ∈ R + , (16)
then for α > 0 small enough, (δ t ) t∈N is also V α -geometrically ergodic.

We now use a law of large numbers ( [START_REF] Meyn | Markov chains and stochastic stability[END_REF] Theorem 17.0.1) on the Markov chain (δ t , (M i,j t ) i∈[1..λ],j∈N ) t∈N to obtain an almost sure divergence of the algorithm.

Proposition 3. Let a (1, λ)-ES optimize problem (2), handling the constraint through resampling. Assume that the distribution H of the random step M i,j t is absolutely continuous with continuous and strictly positive density h, that conditions ( 16) and (15) of Proposition 2 hold, and denote π and µ M the stationary distribution of respectively (δ t ) t∈N and (M i,j t ) i∈[1..λ],(j,t)∈N 2 . Then

[X t -X 0 ] 1 t a.s. -→ t→+∞ σE π×µM ([M ⋆ t ] 1 ) . ( 17
)
Furthermore if E([M ⋆ t ] 2 ) < 0, then the right hand side of Eq. (17) is strictly positive.

Throughout this section we give cases where the assumptions on the distribution of the random steps H used in Proposition 2 or Proposition 3 are verified.

The following lemma shows an equivalence between a non-identity covariance matrix for H and a different norm and constraint angle θ. Lemma 3. Let a (1, λ)-ES optimize problem (2), handling the constraint with resampling. Assume that the distribution H of the random step M i,j t has positive definite covariance matrix C with eigenvalues (α 2 i ) i∈[1.

.n] and take B = (b i,j ) (i,j)∈ [1..n] 2 such that BCB -1 is diagonal. Denote A H,g,X0 the sequence of parent points (X t ) t∈N of the algorithm with distribution H for the random steps M i,j t , constraint angle θ and initial parent X 0 . Then for all k ∈ [1.

.n]

β k [A H,θ,X0 ] k (d) = A C -1/2 H,θ ′ ,X ′ 0 k , (18)
where

β k = n j=1 b 2 j,i α 2 i , θ ′ = arccos( β1cosθ βg ) with β g = β 2 1 cos 2 θ + β 2 2 sin 2 θ, and [X ′ 0 ] k = β k [X 0 ] k for all k ∈ [1..n].
Although Eq. ( 17) shows divergence of the algorithm, it is important that it diverges in the right direction, i.e. that the right hand side of Eq. ( 17) has a positive sign. This is achieved when the distribution of the random steps is isotropic, as stated in the following proposition. Proposition 4. Let a (1, λ)-ES optimize problem (2) with constant step-size, handling the constraint with resampling. Suppose that the Markov chain (δ t ) t∈N is positive Harris, that the distribution H of the random step M i,j t is absolutely continuous with strictly positive density h, and take C its covariance matrix. If the distribution C -1/2 H is isotropic then E π×µM ([M ⋆ t ] 1 ) > 0.

Lemma 3 and Proposition 4 imply the following result to hold for multivariate normal distributions.

Proposition 5. Let a (1, λ)-ES optimize problem (2) with constant step-size, handling the constraint with resampling. If H is a multivariate normal distribution with mean 0, then (δ t ) t∈N is a geometrically ergodic positive Harris Markov chain, Eq. (17) holds and its right hand side is strictly positive.

To obtain sufficient conditions for the density of the random steps to be strictly positive, it is advantageous to decompose that distribution into its marginals and the copula combining them. We pay a particular attention to Archimedean copulas, i.e., copulas defined

(∀u ∈ [0, 1] n ) C ψ (u) = ψ(ψ -1 ([u] 1 ) + • • • + ψ -1 ([u] n )), (19)
where ψ : [0, +∞] → [0, 1] is an Archimedean generator, i.e., ψ(0) = 1, ψ(+∞) = lim t→+∞ ψ(t) = 0, ψ is continuous and strictly decreasing on [0, inf{t : ψ(t) = 0}), and ψ -1 denotes the generalized inverse of ψ,

(∀u ∈ [0, 1]) ψ -1 (u) = inf{t ∈ [0, +∞] : ψ(t) = u}. (20)
The reason for our interest is that Archimedean copulas are invariant with respect to permutations of variables, i.e.,

(∀u ∈ [0, 1] n ) C ψ (Qu) = C ψ (u). ( 21
)
holds for any permutation matrix Q ∈ R n,n . This can be seen as a weak form of isotropy because in the case of isotropy, (19) holds for any rotation matrix, and a permutation matrix is a specific rotation matrix.

Proposition 6. Let H be the distribution of the two first dimensions of the random step M i,j t , H 1 and H 2 be its marginals, and C be the copula relating H to H 1 and H 2 . Then the following holds:

1. Sufficient for H to have a continuous strictly positive density is the simultaneous validity of the following two conditions. (i) H 1 and H 2 have continuous strictly positive densities h 1 and h 2 , respectively. (ii) C has a continuous strictly positive density c. Moreover, if (i) and (ii) are valid, then

(∀x ∈ R 2 ) h(x) = c(H 1 ([x] 1 ), H 2 ([x] 2 ))h 1 ([x] 1 )h 2 ([x] 2 ). (22)
2. If C is Archimedean with generator ψ, then it is sufficient to replace (ii) with (ii') ψ is at least 4-monotone, i.e., ψ is continuous on [0, +∞], ψ ′′ is decreasing and convex on R + , and (∀t ∈ R + ) (-1) k ψ (k) (t) ≥ 0, k = 0, 1, 2. In this case, if (i) and (ii') are valid, then

(∀x ∈ R 2 ) h(x) = ψ ′′ (ψ -1 (H 1 ([x] 1 )) + ψ -1 (H 2 ([x] 2 ))) ψ ′ (ψ -1 (H 1 ([x] 1 )) + ψ -1 (H 2 ([x] 2 ))) h 1 ([x] 1 )h 2 ([x] 2 ).
(23)

Discussion

The paper presents a generalization of recent results of the first author [START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF] concerning linear optimization by a (1, λ)-ES in the constant step size case. The generalization consists in replacing the assumption of normality of random steps involved in the evolution strategy by substantially more general distributional assumptions. This generalization shows that isotropic distributions solve the linear problem. Also, although the conditions for the ergodicity of the studied Markov chain accept some heavy-tail distributions, an exponentially vanishing tail allow for geometric ergodicity, which imply a faster convergence to its stationary distribution, and faster convergence of Monte Carlo simulations. In our opinion, these conditions increase the insight into the role that different kinds of distributions play in evolutionary computation, and enlarges the spectrum of possibilities for designing evolutionary algorithms with solid theoretical fundamentals. At the same time, applying the decomposition of a multidimensional distribution into its marginals and the copula combining them, the paper attempts to bring a small contribution to the research into applicability of copulas in evolutionary computation, complementing the more common application of copulas to the Estimation of Distribution Algorithms [START_REF] Cuesta-Infante | Bivariate empirical and n-variate archimedean copulas in estimation of distribution algorithms[END_REF][START_REF] Salinas-Gutierrez | Using copulas in estimation of distribution algorithms[END_REF][START_REF] Wang | Copula estimation of distribution algorithms based on exchangeable archimedean copula[END_REF].

Needless to say, more realistic than the constant step size case, but also more difficult to investigate, is the varying step size case. The most important results in [START_REF] Chotard | Markov chain analysis of evolution strategies on a linear constraint optimization problem[END_REF] actually concern that case. A generalization of those results for non-Gaussian distributions of random steps for cumulative step-size adaptation ( [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]) is especially difficult as the evolution path is tailored for Gaussian steps, and some careful tweaking would have to be applied. The σ self-adaptation evolution strategy ( [START_REF] Beyer | Toward a theory of evolution strategies: Self-adaptation[END_REF]), studied in [START_REF] Arnold | On the behaviour of the (1, λ)-σSA-ES for a constrained linear problem[END_REF] for the same problem, appears easier, and would be our direction for future research.

Fig. 1 .

 1 Fig.1. Linear function with a linear constraint, in the plane spanned by ∇f and ∇g, with the angle from ∇f to ∇g equal to θ ∈ (0, π/2). The point x is at distance g(x) from the constraint hyperplan g(x) = 0.

Acknowledgment

The research reported in this paper has been supported by grant ANR-2010-COSI-002 (SIMINOLE) of the French National Research Agency, and Czech Science Foundation (GA ČR) grant 13-17187S.


Quentin Grimonprez

Alain Celisse

Guillemette Marot

Analysis of genomic markers: Make it easy with the R package MPAgenomics

Data

Affymetrix genome-wide SNP6 arrays.

About 200 biological samples with two types of profiles :

• copy-number: ∼ 1.8 million probes (SNPs + CN) • allele B fraction: proportion of total signal from allele B (∼ 930.000 SNPs).

Goal

• Create an R package : pipeline for beginners in R to easily perform data analysis from genome-wide SNP arrays. • Calibration method for the segmentation parameter.

Data Normalization

Packages aroma

Allele B fraction

Heterozygous SNPs are kept and the signal is symmetrized. Then, the signal is segmented the same way as the copy-number signal.

Calibration of λ parameter in PELT

• PELT depends on a parameter to calibrate.

• MPAgenomics: automatic calibration of λ.

Calling method

• Assign labels (loss, normal or gain) to segments (copy-number).

• CGHcall package (van de Wiel et al., 2007).

Markers selection Strategy

• Select genomic markers (e.g. SNPs or CNV) associated with a response y.

• Lasso method for sparse selection (few markers) with ρ > 0 :

I i=1 (y i -(Xβ) i ) 2 + ρ P p=1 |β p |

Implementation in MPAgenomics

• Linear regression: HDPenReg for huge amount of variables (HDPenReg: R package, C++ implementation of LARS [START_REF] Efron | Least angle regression[END_REF]).

• Logistic regression: wrapper of glmnet R package [START_REF] Friedman | Regularization Paths for Generalized Linear Models via Coordinate Descent[END_REF].

• Choice of ρ by k-fold cross validation.

Calibration of λ (segmentation)

• PELT default parameter is misleading.

• MPAgenomics: automatic data-driven choice of λ Strategy 1. Grid of λ: 0 < λ 1 < λ 2 < • • • < λ max .
2. Run PELT for each λ i (see Figure 2 left). 3. Choose λ corresponding to the widest range such that the number of segments is constant (> 1).

Sample-specific parameter versus common λ

1. Common λ:

• Compute the signal-to-noise ratio (SNR) for each profile.

• Cluster profiles according to SNR (Gaussian mixture).

• For each cluster, choose λ.

Sample-specific λ:

MPAgenomics provides an automatic choice of λ for each profile.

Sample-specific parameter versus common λ

Common λ within each cluster is misleading (Figure 2

Figure 1 :

 1 Figure 1: Top : Copy-number profile. Bottom : Allele B fraction profile : homozygous SNPs (red), heterozygous SNPs (black).

  right).

Figure 2 :

 2 Figure 2: Left: Number of segments versus λ (chromosome 1). Right: Ranges of optimal λ for each profile. Color indicates a given signal-to-noise ratio (red¡green¡blue).