Lucas Ammann

Pascal Barla

Xavier Granier

Gael Guennebaud

Patrick Reuter

Surface Relief Analysis for Illustrative Shading

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Features are extracted and combined at multiple scales to depict relevant relief details (c). Finally, radiance scaling is added to enhance the relief based on the curvature at each feature (d).

Introduction

During the last decade, the bulk of computer graphics applications has made use of 3D objects exhibiting an everincreasing amount of shape details. Such geometric complexity may either come from densely scanned real-world objects, or it may be created by artists in modern digital modeling and sculpting software. In both cases, there is a growing need for shading tools that help visualize, enhance, or exaggerate complex surface relief patterns automatically.

This need firstly arises in scientific illustration domains. In Cultural Heritage, example applications include the study of detailed inscriptions found in engraved stones or of shallow traces over human bones. In Cartography, mountain relief is acquired with an ever-increasing resolution, † {ammann,barla,guenneba,granier,reuter}@labri.fr which then requires specific coloring and shading techniques for a proper depiction. The same requirement has recently emerged with manually-created geometry. Artists nowadays spend a lot of time using modeling and sculpting software to work out fine shape details such as skin pores, wrinkles or grain. Surface shading techniques that properly depict this geometric complexity without requiring additional manual effort are then necessary.

In this paper, we consider relief in a standard way as a height field relative to a smooth (and not necessarily planar) base surface. This is a reasonable assumption both in scientific applications where surface details are often very shallow (as in Archeology) or naturally defined as height data (as in Cartography), and in artistic applications where most surface details are stored in displacement maps.

A relief feature is then meant to be a region on the base surface where the relief is prominent. In particular, we are in-terested in three types of features -convexities, inflections and concavities -and the ways these are nested at multiple scales. Take for instance a relief ridge pitted with small cavities, as shown on the shoulder of the statue in Figure 1(a): a point anywhere inside one of the small pits belongs to a concave feature at a small scale, and to a convex feature at a larger scale. Conveying surface relief through illustrative shading then amounts to assigning a proper shading value to each surface point, in a way that depends on the multi-scale features it belongs to. This raises two major issues: 1) how to locate relief features in the neighborhood of a surface point? 2) how to combine features found at various scales (i.e., at increasing neighborhood sizes)?

We address these issues by introducing an efficient bottom-up, multi-scale relief analysis technique based on cubic polynomial fitting. It locates relief features (concavities, convexities and inflections) in a neighborhood around each surface point, and provides a straightforward solution to the multi-scale combination problem. We illustrate the benefits of our approach by adapting shading techniques to take advantage of non-local relief information available at each surface point. An example is shown in Figure 1(b-d), where shading color is assigned based on the distance to the closest feature center, first at a single (large) scale, then at multiple scales. The result is a far more legible depiction of relief features occuring at various scales (including the small pits on the shoulder). Our system only requires a short pre-process for the completion of the bottom-up analysis step. It then allows users to explore surface relief in real-time, by varying the contribution of each scale, or by modifying shading and filtering parameters.

Previous work

One of the earliest shading methods explicitly designed to convey shape features is accessibility shading [Mil94]. Together with its most popular variant ambient occlusion (e.g., [PG04]), these methods produce images that depict cavities by darkening hardly-accessible surface points. However, they have two main drawbacks. First, they demand long pre-processing times when done accurately in object space. But most importantly, they provide limited control in terms of surface feature depiction: shading and analysis are not separated and mainly convey deep concavities with very smooth shading variations.

More accurate surface measurements are obtained with differential geometry operators [dC76]. In particular, the second-and third-order tensors provide information about curvature and curvature variations at a surface point. Such measurements only apply to small surface neighborhoods though (infinitesimal neighborhoods for smooth surfaces). They have thus traditionally been used for identifying surface points that exactly lie on extremal curves such as ridges and valleys [START_REF] Ohtake | Ridge-valley lines on meshes via implicit surface fitting[END_REF] or demarcating curves [START_REF] Kolomenkin | Demarcating Curves for Shape Illustration[END_REF]. Some methods have used similar operators for larger surface neighborhoods (e.g., [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF]CPG09]) by first fitting a local plane to the surface, and then fitting a quadratic function on relative height values. When the surface "folds-over" itself, complex issues arise that can be overcome using total leastsquare fits [GGG08]. Such analysis methods are limited to purely local surface feature measurements, while we are interested in locating salient relief features, which we recall are to be found in the neighborhood of each surface point. In particular, these methods provide no simple way to know the distance to the nearest feature center of a given type, to combine features found at multiple scales, or to filter out less prominent features (e.g., shallow relief).

Another approach for analysing relief at multiple scales consists in decomposing surface normals into different layers, which also avoids fold-over issues. Such a decomposition has first been applied to normal enhancement [CST05, ZCF * 10], where a single relief layer is manipulated to exaggerate or attenuate surface details through shading. Similar enhancement techniques have been applied to polynomial texture maps [MWGA06]. In both cases, only two scales are considered. These methods have been extended to multiple scales in the exaggerated shading technique [RBD06]. Different levels of smoothed normals are used to align a single light source at grazing angles. Shading values are then computed at each scale through half-Lambertian shading, and combined with a weighted sum to exaggerate multiscale surface relief. The drawback with this approach is that scale combination depends on the choice of light direction and shading parameters, which not only produces artifacts when the light is moved around, but also makes it difficult to control which type of feature is depicted.

Normal variations in image space have also been used to convey surface details by letting them drive variations of incoming radiance (e.g., [VPB * 10b]). This is an improvement compared to exaggerated shading, since it works for arbitrary materials and illumination and is devoid of temporal artifacts. However, the method is confined to a single scale per point, and provides limited control over the type of depicted feature as before.

When an accurate depiction of surface relief is targeted, it is preferable to decompose surface geometry into base and relief layers. Various methods have been proposed to perform this decomposition (e.g., [ZTS09]). The goal of this paper is not to present a novel decomposition though, but rather to provide analysis solutions specifically targeted to the relief layer. Only a few existing techniques have tackled this problem. The prominent field technique [KST09b] combines output of second-and third-order tensors to identify a direction field along which prominent relief features are likely to be located. This direction field has proven to be useful for feature-aware smoothing, curvature-based shading, and more recently line-based rendering [KST09a]. A similar solution has been proposed to identify feature lines in image-space, using fitting techniques [VVC * 11]. The advantage of using a direction field is that it locates more accurately surface features around a point of interest. Unfortu-nately, as depicted in Figure 2 using a single direction also introduces artifacts around field singularities, and since the field is likely to change at different scales, it is not clear how to perform scale combination.

Overview

Our approach focuses on surface relief and assumes input geometry to be decomposed into base and relief layers. We first locate relief features in the neighborhood of a surface point, at a number of scales and directions. This is done thanks to cubic polynomial fitting as explained in Section 4.1. We then combine features together and filter them in a way that preserves only the most pertinent data. This is permited by a weighting scheme based on the information conveyed by individual fits, as explained in Section 4.2. Shading techniques are finally adapted to depict properties of the relief features found around each surface point, as described in Section 5.

Details of our GPU implementation are given in Section 6: it requires a short pre-process for analysis, after which it runs in real-time with interactive scale combination, filtering and shading abilities. Our solution does not rely on any parametrization of the base surface, and we demonstrate our system on both height-fields and detailed 2D manifolds. Moreover, for dynamic scenes, there is no need to update the analysis during animation since in most cases, only the base surface is deformed. The scalar function h corresponds to the relief of the surface. For instance, sculpting tools are able to directly provide a base surface with a displacement map. With acquired objects though, the base and relief layers must be extracted in pre-process. The presentation of this step is out of the scope of our paper, and we refer to, e.g., Zatzarinni et al. [ZTS09] for such a base/relief decomposition technique. For the sake of simplicity, and without loss of generality, we assume that B is provided as a dense (regular or irregular) triangular mesh, while h is given as a scalar value per vertex.

Relief analysis

1D Analysis

The key idea of our approach is to analyse relief along 1D neighborhoods of the base surface, in multiple directions around each point b b b. The main reasons for this choice are that 1) the localisation of relief feature centers (e.g., curvature extrema) around b b b is made considerably simpler in one dimension, and 2) no 2D parametrization is required. This strategy is also significantly faster for large scales.

In practice, we perform the analysis at mesh vertices, in multiple directions θ i uniformly distributed in the [0, π[ range. We define t t t i = (cos θ i , sin θ i ) as the corresponding direction vector that lies in the tangent plane of the base surface. A 1D neighborhood is obtained as the intersection of the base surface with the plane spanned by n n n B and t t t i , and going through b b b (see Figure 3(a)). The reference direction θ 0 = 0 can be chosen arbitrarily for each point b b b without any impact on quality because, in our approach, we combine multipe directions to compute a single feature property. A polynomial is given by h i j (t) = u 0 + u 1 t + u 2 t 2 + u 3 t 3 . As illustrated in Figure 3(b), h i j is computed by a standard least square minimization applied to a finite set of samples t k uniformly spread in the range [-s j /2, s j /2]:

i (t)) = h •b b b i (t).
h i j = arg min P ∑ k (P(t k ) -h(b b b i (t k ))) 2 .
(1)

Details on how to efficiently solve these minimization problems are given in Section 6, where we show how to take advantage of sampling regularity and incremental calculus.

Feature extraction One of the main reasons for making use of an analytic representation is that it permits to identify entire relief features. The fitted cubic polynomial provides us with analytic expressions of the first three order derivatives h ′ i j , h ′′ i j , and h ′′′ i j , as well as the curvature κ i j :

κ i j (t) = h ′′ i j (t) (1 + h ′ i j (t) 2 ) 3 2 . ( 2
)
We use these differential quantities to identify the positions of the convexity (t ⌢ i j ), concavity (t ⌣ i j ), and inflexion (t ∼ i j ) points. We define the first two as the locations of the curvature extrema which are obtained as the zero-crossings of the curvature variation. More precisely, t ⌢ i j and t ⌣ i j are the roots of the rational polynomial κ ′ i j (t) for which explicit formulas are given in the Appendix. We observe that these roots do not coincide with the local extrema of the cubic defined as the zerocrossings of the first-order derivative h ′ i j , even though they are often very close (see Figure 3(b)). This choice is particularly important when the cubic is monotone but not linear: curvature extrema can still be found whereas the value extrema do not exist. Curvature extrema thus constitute a more robust choice in general. The position of the inflexion point t ∼ i j is given by the zero-crossing of the second-order derivative (P ′′ i j (t ∼ i j ) = 0).

Normalized distances and error measure Features extracted from the polynomial approximation may not always be pertinent. First, specific feature locations may be identified outside the support used to fit the relief data. We thus compute a truncated and normalized distance for each type of relief feature (⋆ ∈ {⌢, ⌣, ∼}):

d ⋆ i j = 2|t ⋆ i j | s j , ( 3
)
where ⌊x⌉ is a function that clamps x to the [0, 1] range.

Second, the detailed surface normals may substantially deviate from the plane used to intersect the base surface geometry. Indeed, with a large normal deviation, a small change in direction implies a large change in the fitted heights. We take into account this source of error by measuring for each fit the average normal deviation η i j :

η i j = 1 s j ∑ k |n n n(t k ) • (t t t i ×n n n B )| , ( 4
)
where n n n is the normal of the detailed surface S, and t t t i × n n n B is the normal of the plane supporting the cubic. Both types of error measures are visualized in Figure 3(c) as a function of direction θ i at a single scale j.

Combination and Filtering

Our multi-scale and multi-direction analysis provides a dense relief description around a surface point, which ensures that most nearby relevant relief features are properly identified by at least one fitted polynomial. We now present combination and filtering mechanisms that make such information exploitable by subsequent shading techniques.

Feature combination

The key idea of our approach is to combine scalar feature values f i j instead of polynomials themselves, to yield a single scalar F. For instance f i j could be taken to be the convexity curvature κ i j (t ⌢ i, j ), yielding an average curvature K. Combination is done in two steps. We first average feature values at each scale over all directions:

F j = ∑ i w i j f i j , w i j = w i j ∑ i ′ w i ′ j , (5)
using w i j = (1 -η i j )(1 -d ⋆ i j
) for the weights. The factor 1 -η i j favors directions of low normal deviation for which extracted feature values are most pertinent. This may be seen as a generalization of principal curvature directions, since for a quadratic surface, they exactly correspond to the directions of zero normal deviation. For choices of f i j involving either t ⌢ i j , t ⌣ i j or t ∼ i j , the weights 1d ⋆ i j also decrease as the feature location approaches the limit of the support. In practice, any decreasing functions of η i j and d ⋆ i j could be used. The F j from Equation 5 are then averaged over scales:

F = ∑ j α j W j F j , W j = W j ∑ j ′ W j ′ . ( 6
)
with W j = ∑ i w i j . This weighting provides for a natural balance between features at each scale: the influence of a given scale depends on the average pertinence of the angular analysis at this scale. The parameter α j is introduced for usercontrolled filtering. By default it is set to 1, and its adjustment will be discussed below.

Figure 4 illustrates our approach at individual scales (Equation 5). If we take f i j := κ i j (0), it boils down to a local estimation of mean curvature. Like total least square fits [GGG08], it permits to avoid fold-over issues at large scales (we use 8 directions for each scale). However, the main benefit of our approach is in estimating properties of features in the neighborhood of each surface point. To illustrate this, we choose f i j := κ i j (t • i, j ), which corresponds to the curvature of the closest relevant feature, i.e., the closest convexity (t • i, j = t ⌢ i j ) or concavity (t • i, j = t ⌣ i j ). It results in a much more segmented measurement of curvature, each point displaying the curvature of the convexity it belongs to.

Figure 5 illustrates our approach at multiple scales (Equation 6). Here we use another feature value: f i j := 1d ∼ i j , the distance to the closest inflection point. Contrary to exaggerated shading [RBD06], our scale combination is independent of the choice of lighting or shading model. More importantly, our blending weights (the W j in Equation 6) are locally adjusted over the object surface based on fitting pertinence, while their are based solely on scale, thus leading to a more legible end-result.

User-controlled filtering Users may desire to select specific features during the combination process; this is made possible through the α j parameter introduced in Equation 6. In practice, α j decreases the contribution of the scale s j relative to the other scales and to a fall-back value F B . This is done thanks to a simple linear interpolation:

Fα = lerp(F B , F, α) = F B + (F -F B )α, ( 7
)
where α is the average of all the α j . In other words, this procedure ensures that in regions where none, or little, of the user-selected features are found, we fall back to a "background" value F B .

The choice of α j depends on the kind of features the user wants to preserve. We suggest the following feature filter:

α j = ψ std i (η i j ) ∑ i w i j κ i j (t ⌢ i j ) , ( 8
)
where ψ is a user controlled sigmoid function allowing for a smooth transition between the selected and rejected parts.

In this paper we used a cubic smooth step function. The first factor retains anisotropic relief features defined as the standard deviation of the normal deviation η i j over the set of directions. The second factor selects regions for which the closest relevant feature has a high mean curvature. Because we make use of non-local curvature measurements, entire

(in red, green and blue for concavities, inflections and convexities) serve to identify features for direct visualization. The curvature function c i j (in blue) is used instead of κ i j (in orange) as it offers better user control. The thin blue lines represent variations of the exponent γ.

features are selected even if filtering is applied per surface point. As illustrated in Figure 6, it preserves much better relief structure compared to filtering based on local mean curvature, which tends to cut-off relief features. In our experience, this filter works remarkably well to select engraved symbols.

Shading techniques

In the previous section, extracted relief features have been displayed with simple colors for clarity of exposition. We now explain how more complex illustrative shading techniques are adapted to take advantage of our analysis. Thus, all subsequent figures present variants of our approach. In all results, we use 8 directions and between 3 and 5 scales.

Normal enhancement

The availability of relief normals is necessary for all shading techniques regardless of the target application. However, using the original relief normal raises coherence issues whenever filtering is activated (i.e., α < 1) . Our solution consists in estimating an averaged normal N N N alongside the computation of F. We compute one local normal n n n i j for each fitted polynomial h i j (see Figure 7 world space. The coordinates of n n n i j in the plane spanned by h i j are directly given by:

n n n i j = (-u 1 , 1) (-u 1 , 1) , ( 9
)
where u 1 is the linear coefficient of h i j .

The averaged normal N N N is then obtained by applying Equation 6 to each coordinate of n n n i j in world-space and then renormalizing. The final filtered normal is computed via spherical interpolation: N N Nα = slerp(n n n B ,N N N, βα), where β ∈ R + is a user-controlled parameter reproducing the normal enhancement technique [CST05]. With β = 1, we obtain the estimated normal. For β ∈ [0, 1], relief is attenuated, while for β > 1, relief is exaggerated. This is illustrated in Figure 8 on an acquired pipe model. The sole benefits of our relief analysis in this case is that it provides a precise control over which surface features are enhanced through scale combination and filtering. Note that the normals of the original detailed surface are easily reconstructed by using the finest scale only and α = 1.

Accessibility shading

The entire fitted polynomial may also be used for local visibility evaluation. We propose a variant of accessibility shading inspired by ambient occlusion. It computes a local visibility value v i j ∈ [0, 1] for each 1D polynomial as the integra- tion of the cosine term over the range of visible directions:

v i j = φ + i j φ - i j cos(φ)| sin(φ)|dφ = 1 - 1 2 cos 2 (φ - i j ) + cos 2 (φ + i j ) . (10)
The integration bounds φ ± i j correspond to the minimum and maximum visibility angles around the normal n n n i j (Figure 7(a)). Their derivation is given in the Appendix. The local visibility terms v i j are combined using Equation 6 to yield an average visibility V which can be arbitrarily scaled by a factor λ ∈]0, 1] to exaggerate the effect. The final filtered occlusion is computed via linear interpolation:

Vα = lerp(1,V, α).
Figure 9 illustrates this technique on a spherical box. Our approach has two advantages over classical accessibility techniques when applied to surface relief: it is fast; and it provides an accurate control over which cavities are conveyed thanks to both multi-scale combination and usercontrolled filtering. This effect emphasizes the non purelylocal aspect of our analysis, since it exploits the entire polynomial information included in the support size.

Feature visualization

In scientific illustration, shading is used to help visualize relief features unambiguously. To this end, we assign one shading function per feature type (ρ ⌢ , ρ ⌣ and ρ ∼ ), and blend them based on proximity. To do so, we define three membership functions m ⌢ i j , m ⌣ i j and m ∼ i j that tell us how close we are relative to each feature center inside each 1D fitted polynomial. Let's suppose that t ⌣ i j < t ⌢ i j ; then we have:

m ⌢ i j (t) = t ∼ i j -t ∆ i j ; m ⌣ i j (t) = t-t ∼ i j ∆ i j ; m ∼ i j (t) = 1- |t ∼ i j -t| ∆ i j
where ∆ i j = |t ⌣ i j -t ⌢ i j |/2, and recall that ⌊x⌉ is a function that clamps x to the [0, 1] range. In the case t ⌣ i j > t ⌢ i j , memberships are obtained by exchanging values for m ⌣ i j and m ⌢ i j . As shown in Figure 7(b), the membership functions form a partition of unity: ∀t, m ⌢ i j (t) + m ⌣ i j (t) + m ∼ i j (t) = 1. In practice, we evaluate these memberships only at t = 0 and plug the resulting values into Equation 6. Since feature combination is a linear process, the resulting averaged memberships M ⌢ , M ⌣ and M ∼ also form a partition of unity. This is a handy property, since it permits to unambiguously assign a reflectance to each membership triplet:

ρ = M ⌢ ρ ⌢ + M ⌣ ρ ⌣ + M ∼ ρ ∼ . ( 11
)
The final filtered reflectance is obtained by linear interpolation: ρ α = lerp(ρ B , ρ, α), with ρ B the base reflectance.

Figure 10 demonstrates feature visualization on a statue model using lit-sphere shading [SMGG01]. We illustrate another important advantage of our approach here: it provides distance information to nearby concavities, convexities and inflections within a single analysis. Again, this is made possible thanks to our ability to detect these features in the neighborhood of each surface point.

Radiance scaling

Recent shading techniques have demonstrated enhancement capabilities with arbitrary material and illumination. For instance, radiance scaling [VPB * 10a] correlates shading intensity variations produced by each incoming light direction with surface normal variations. Thanks to our approach, curvature extrema are easily located around a surface point. It allows us to not only convey the curvature of nearby features, but also provides a more precise control. Indeed, as can be seen in Figure 7(c), the highest values of κ i j (in orange) are mostly located in a narrow neighborhood around each feature location. Thus, instead of directly using κ i j , we remap it to a curvature function c i j given by

c i j (t) = κ i j (t) 1 - t-t • i j ∆ i j 2 γ , (12)
where ∆ i j is the same as in Equation 11, γ ∈ R + is a userspecified parameter that controls the sharpness of concave and convex regions, and, as before, t • i j corresponds to the location of the closest relevant feature (t ⌣ i j or t ⌢ i j ). As illustrated in Figure 7(c), c i j (t) = κ i j (t) for t ∈ {t ⌣ i j ,t ⌢ i j ,t ∼ i j }, but exhibits a smoother behavior, in a way controlled by γ.

The curvatures c i j are combined using Equation 6 to obtain the average curvature C. The final filtered curvature is computed via simple interpolation (Cα = lerp(0,C, α) as in Equation 7), and it is used in the radiance scaling technique to modulate incoming light intensities at each surface point.

Figure 11 illustrates this technique applied to a detailed face model sculpted in Autodesk Mudbox c . The main advantage of our approach here is to provide a direct control over the magnitude (via filtering) and the thickness (through γ) of conveyed features. Each of them cannot be achieved using a purely local curvature estimator.

Figures 1 and12 show more sophisticated shading examples, where our illustrative shading techniques are combined in different ways. The accompanying video and image archive present all our results with more parameter variations and in high resolution.

Implementation

In our prototype implementation, the analysis (Section 4.1) is performed on mesh vertices in a preprocess and updated every time a new scale is added. The cubic coefficients for each h i j along with normal deviations η i j are stored in vertex buffers. By trading memory for speed, this strategy enables real-time feature combination (Section 4.2) in a vertex shader followed by per-pixel illustrative shading (Section 5).

The fitting of the polynomials is a critical step of our algorithm. With a naive implementation, complexity would depend on the number of scales. We propose a more efficient strategy that makes complexity depend on the largest scale only (see Table 6 for performance statistics). Given a vertex v v v and a direction θ i , the key idea is to perform fitting for all scales incrementally, from the smallest to the largest scale. Relief samples t k are then taken at equal intervals ∆t , which we set to half the average edge-length in our system. In this way, the samples t k are the same at all scales and the fitting of a cubic (Equation 1) boils down to a normal equation:

A A A j u u u i j (v v v) = r r r i j (v v v) , A A A j = ∑ k L(t k ) T L(t k ) , r r r i j (v v v) = ∑ k h(b b b i (t k ))L(t k ) T ,
where u u u i j (v v v) is the vector of unknown polynomial coefficients, and L(t k ) = (1, t k , t 2 k , t 3 k ) is the cubic basis vector. From these equations, it is clear that the covariance matrices A A A j only depend on their corresponding scale s j . They are thus computed and pre-inverted once and for all. Likewise, the right hand sides r r r i j (v v v) and deviation measures η i j (Equation 4) are easily computed incrementally from r r r i( j-1) (v v v) and η i( j-1) . Fitting is performed on the GPU using a minimalistic half-edge data structure to walk on the surface. Table 6 reports some performance statistics of our OpenCL implementation running on a Nvidia GTX 480 graphics card.

Discussion

Depicting surface relief through illustrative shading requires to gain knowledge and control over relief features at various scales. Our approach permits to identify, filter and combine three types of relief features at the small expense of a short, automatic pre-process; after which all subsequent manipulations and rendering occur in real-time. Our current implementation works at the vertex level, but it would be easily extended to take 3D models with displacement maps as input, performing the analysis and storing its outcome per texel.

Another limitation is that we interpolate averaged feature values F per pixel, which lets appear tesselation in closeup views. One solution could be to interpolate h i j per pixel; however, the interpolation of such directional cubic polynomials is not straightforward. Our method could also be adapted to other kinds of surface representations, provided we can intersect the surface with a plane, and march on the resulting curve. This includes for instance raw point clouds reconstructed with implicit fitting methods.

Apart from these technical issues, the main limitation of our approach is that it assumes that relief information can be represented by a single height-field. Although this is true for many kinds of surface reliefs, this is not always the case. One solution would be to interleave our multi-scale analysis with 3D multi-scale decompositions (e.g., [START_REF] Kobbelt | Point-based multiscale surface representation[END_REF]), though it might require longer preprocessing times. Our approach also neglects the curvature of the base surface. While it is assumed to be flat most of the time, one could imagine to cut the neighborhoods at regions of sharp transitions [START_REF] Schmidt R | Interactive decal compositing with discrete exponential maps[END_REF].

Conclusion

We have presented a novel analysis technique of surface relief that identifies three types of features in the neighborhood of each surface point, and provides multi-scale combination mechanisms that favor most pertinent features. The direct benefit of our approach is to gain a direct global control over the type of relief features depicted through illustrative shading techniques: the influence of each scale is easily tuned, features are efficiently filtered out based on their importance and/or based on our anisotropic metric, and a variety of feature properties can be selected for illustration. This is made possible by our cubic polynomial fitting approach that identifies features at multiple scales and orientations, without the necessity of a complex parametrization.

Our system thus provides a convenient solution for the depiction of surface relief, either for its exploration as in scientific illustration, or for its exaggeration as in visual effects. It relies on a base/relief decomposition though, and we believe that different solutions to this problem may be chosen depending on the target application, considering how relief has been generated. Has it been engraved or eroded? Does it include cracks? Is it composed of different materials? Although we believe concavities, convexities and inflections are very common features, they might not always be adapted. An exciting avenue of future work would be to consider alternative feature types for expert applications. As long as the fitting functions used to identify them are one-dimensional, our method will be applicable with little changes. Finally, we believe such relief analysis technique could be used in many other applications, for instance, for geometric-aware texture synthesis [MKC * 06].

Eurographics

  Symposium on Rendering 2012 Fredo Durand and Diego Gutierrez (Guest Editors) Volume 31 (2012), Number 4 Surface Relief Analysis for Illustrative Shading Lucas Ammann † , Pascal Barla, Gaël Guennebaud, Xavier Granier, Patrick Reuter Inria -Univ. Bordeaux -IOGS -CNRS (a) Diffuse shading (b) Coarse scale only (c) All scales (d) + Radiance scaling

Figure 1 :

 1 Figure 1: Feature-based shading: Given a detailed surface (a), we analyze its relief to locate relief features in the neighborhood of each surface point (b). We focus on three types of features: convexities, concavities, and inflexions, shown on the right half with blue, red and white colors respectivelly. Extracted information is used to assign them different shading functions: here we use three different lit-spheres, shown on the left half. An additional accessibility shading effect helps convey relief cavities. Features are extracted and combined at multiple scales to depict relevant relief details (c). Finally, radiance scaling is added to enhance the relief based on the curvature at each feature (d).

Figure 2 :

 2 Figure 2: Artifacts due to singularities. As opposed to isotropic measurements (a), the use of a directional field to analyse relief produces artifacts at it singularities (b).

  Our analysis takes as input a surface S defined as a smooth base surface B displaced along its normal field n n n B by a height function h. Mathematically, any point b b b ∈ B yields a point p p p on the surface S with: p p p = b b b + h(b b b)n n n B (b b b) .

For

  each direction, this intersection yields a base polygonal parametric curve b b b i which is parametrized in term of arc-length on the base surface. Every b b b i (t) corresponds to a point on the base surface at a distance t along the polygonal curve, with associated height value h(b b b

  Polynomial fitting The 1D relief function h • b b b i may contain a variety of features at multiple scales. In order to identify them, we approximate h • b b b i by a set of cubic polynomials fitted for multiple support sizes s j . Higher values of s j correspond to coarser scales. Cubic-polynomials are particularly well suited to our purpose since they exactly exhibit a pair of convex and concave regions separated by an inflection. It permits to locate and characterize such neighboring features analytically.

  Figure 3: 1D analysis. (a) The curves b b b i on the detailed surface S are extracted in multiple directions by intersecting planes defined relative to the base surface B. (b) The height along a curve h •b b b i on a support of size s j is fitted with a cubic polynomial h i j from which we compute the curvature function κ i j and feature center locations t ⋆ i j (⋆ ∈ {⌢, ⌣, ∼}) analytically. (c) For each direction, we compute the distance d ⋆ i j to a particular feature center and the normal deviation η i j of the detailed surface.

Figure 4 :

 4 Figure 4: Curvature comparisons. The mean curvature obtained from the least square fit of an algebraic sphere (a) yields similar results than averaging our curvature κ i, j (0) at the current point (b). On the other hand, averaging the curvature of the closest convexity or concavity (κ i, j (t • i, j )) exhibits a much more segmented curvature measurement (c).

Figure 5 :

 5 Figure 5: Scale combinations: Diffuse shading is darkened proportional to the distance to curvature inflections, in (a) with only a coarse scale, in (b) using only a fine scale, and in (c) and (d) combining all 3 scales together using our adaptive weights and uniform weights respectively. The use of a single scale fails to convey all the relief features in a legible way.

Figure 6 :Figure 7 :

 67 Figure6: Feature filtering. In (a), surface relief height is displayed with a red (lowest) to blue (highest) gradient (i.e., we use f i j := h i j (0) in Equation5). It is filtered based on (b) the average curvature at the evaluation point (α j = ψκ(∑ i w i j κ i j (0)), and (c) our feature filter. The zoom insets (d) shows how our non-local criteria better preserves feature integrity.

Figure 8 :

 8 Figure 8: Normal enhancement: The original normals of an object (a) are enhanced at a coarse scale (b), then at all scales (c-d) using different enhancement parameter values.

L

  . Ammann & P.Barla & G. Guennebaud & X. Granier & P.Reuter / Surface Relief Analysis for Illustrative Shading (a) Diffuse shading (b) Coarse scale only (c) All scales (d) + feature filter

Figure 9 :

 9 Figure 9: Accessibility shading: Starting from (a) standard diffuse shading, we add accessibility shading (using λ = 8) (b) at a coarse scale first, then (c) at all scales, and finally (d) we apply our feature filter.

  (a) Original shading (b) Reflectance functions (c) Fine scale only (d) All scales

Figure 10 :

 10 Figure 10: Feature visualization: When using a single litsphere (a) to render a statue, we miss many of its features. Instead, we assign one lit-sphere per feature, showing results either at a coarse scale (b), or a fine scale (c), or all scales combined (d), which provides a more balanced depiction of surface details.

LFigure 11 :

 11 Figure 11: Radiance Scaling: We enhance the original shading (a) of a detailed face model by applying radiance scaling to all scales (b-c) using different values for sharpness, and then adding our anisotropic feature filter (d) to select the most prominent features. The grey images show the remapped curvature function c i j .

LFigure 12 :

 12 Figure 12: The original diffuse shading (a) of the scanned portion of a cave wall hardly reveals the engraved shape. After the analysis, and assignment of different lit-spheres to the extracted features (b) we can provide a more legible rendering using radiance scaling (RS) and accessibility shading (AS) (c). Filtering removes most of the unwanted features (d).

Table 1 :

 1 Analysis timings for different scales (given in terms of number of samples) and 8 directions. The last column reports the time in second to decompose the input mesh into a base/relief layers with our implementation of [ZTS09].

	model	scales		base/relief
	(#vertices)	{35} {75} {35,55,75} decomposition
	Box (180k)	2.4 5.2	5.8	12.5
	Head (667k) 11.8 25.2	26.9	66.3

c 2012 The Author(s) c 2012 The Eurographics Association and Blackwell Publishing Ltd.

Acknowledgments

This work has been supported by the ANR SeARCH project (ANR-09-CORD-019), and the European Community's Seventh Framework Program [FP7 -2007/2013] under the Grant Agreement 270404. We would also like to thanks Robert Vergnieux, Romain Vergne and Nicolas Mellado for their insightful discussions. The models of figures 4, 8, 9, 10 are provided courtesy of AIMSHAPE Shape Repository, and figure 5 courtesy of M. Kolomenkin.

Appendix

Relief analysis -feature locations As explained in Section 4.1, the convexity and concavity locations t ⌢ , t ⌣ are the roots of the rational polynomial κ ′ (t), where κ (Equation 2) corresponds to the curvature of a cubic polynomial h(t) = u 0 + u 1 t + u 2 t 2 + u 3 t 3 . If u 3 = 0, then there exist only two real roots, t ⌢ and t ⌣ , given by:

The inflexion position t ∼ is then obtained as the unique root of the second-order derivative polynomial (h ′′ (t ∼ ) = 0):

On the other hand, if u 3 = 0 but u 2 = 0, then there is only a unique concavity or convexity point at t = -u 1 /2u 2 , and no inflexion point can be detected. If u 2 is also zero, then no feature points can be extracted. In these cases, the respective undefined d ⋆ i j values are set to 1.

Accessibility shading -integration bounds

The computation of the accessibility term of Equation 10relies on the determination of the cosine of the two integration bounding angles φ ± . In the case of a cubic polynomial h(t), they are given by the clamped dot product between the normal direction n n n (as defined in Equation 9), and the two normalized vectors (x x x ±o o o), where x x x ± = (t ± , h(t ± )) are the two extremity points of the curve visible from the observer location o o o = (0, -u 0 ) (see Figure 7). Since the curve is clamped to the range [-s j /2, s j /2], they either coincide with the boundaries of the curves, or with the point y y y = (ty, h(ty)) for which the vector (y y yo o o) is tangent to the curve at y y y. It is obtained by solving for the cubic equation (-h ′ (ty), 1) • (y y y -o o o) = 0, which admits a single non null root: ty = -u 2 /(2u 3 ). The two positions t -, t + of the visible extremities x x x ± = (t ± , h(t ± )) are finally given by:


Alexandre Chotard

Anne Auger

Nikolaus Hansen

Cumulative Step-size Adaptation on Linear Functions: Technical Report

Keywords: CSA, cumulative path, evolution path, evolution strategies, step-size adaptation

The CSA-ES is an Evolution Strategy with Cumulative Step size Adaptation, where the step size is adapted measuring the length of a so-called cumulative path. The cumulative path is a combination of the previous steps realized by the algorithm, where the importance of each step decreases with time. This article studies the CSA-ES on composites of strictly increasing functions with affine linear functions through the investigation of its underlying Markov chains. Rigorous results on the change and the variation of the step size are derived with and without cumulation. The step-size diverges geometrically fast in most cases. Furthermore, the influence of the cumulation parameter is studied.

Introduction

Evolution strategies (ESs) are continuous stochastic optimization algorithms searching for the minimum of a real valued function f : R n → R. In the (1, λ)-ES, in each iteration, λ new children are generated from a single parent point X ∈ R n by adding a random Gaussian vector to the parent, X ∈ R n → X + σN (0, C) .

Here, σ ∈ R * + is called step-size and C is a covariance matrix. The best of the λ children, i.e. the one with the lowest f -value, becomes the parent of the next iteration. To achieve reasonably fast convergence, step size and covariance matrix have to be adapted throughout the iterations of the algorithm. In this paper, C is the identity and we investigate the so-called Cumulative Step-size Adaptation (CSA), which is used to adapt the step-size in the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Ostermeier | Step-size adaptation based on non-local use of selection information[END_REF][START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF]. In CSA, a cumulative path is introduced, which is a combination of all steps the algorithm has made, where the importance of a step decreases exponentially with time. Arnold and Beyer studied the behavior of CSA on sphere, cigar and ridge functions [START_REF] Arnold | Performance analysis of evolutionary optimization with cumulative step length adaptation[END_REF][START_REF] Arnold | On the behaviour of evolution strategies optimising cigar functions[END_REF][START_REF] Arnold | Cumulative step length adaptation on ridge functions[END_REF][START_REF] Arnold | Evolution strategies with cumulative step length adaptation on the noisy parabolic ridge[END_REF] and on dynamical optimization problems where the optimum moves randomly [START_REF] Arnold | Random dynamics optimum tracking with evolution strategies[END_REF] or linearly [START_REF] Arnold | Optimum tracking with evolution strategies[END_REF]. Arnold also studied the behaviour of a (1, λ)-ES on linear functions with linear constraint [START_REF] Arnold | On the behaviour of the (1,λ)-es for a simple constrained problem[END_REF].

In this paper, we study the behaviour of the (1, λ)-CSA-ES on composites of strictly increasing functions with affine linear functions, e.g. f : x → exp(x 2 -2). Because the CSA-ES is invariant under translation, under change of an orthonormal basis (rotation and reflection), and under strictly increasing transformations of the f -value, we investigate, w.l.o.g., f : x → x 1 . Linear functions model the situation when the current parent is far (here infinitely far) from the optimum of a smooth function. To be far from the optimum means that the distance to the optimum is large, relative to the step-size σ. This situation is undesirable and threatens premature convergence. The situation should be handled well, by increasing step widths, by any search algorithm (and is not handled well by the (1, 2)-σSA-ES [START_REF] Hansen | An analysis of mutative σ-self-adaptation on linear fitness functions[END_REF]). Solving linear functions is also very useful to prove convergence independently of the initial state on more general function classes.

In Section 2 we introduce the (1, λ)-CSA-ES, and some of its characteristics on linear functions. In Sections 3 and 4 we study ln(σ t ) without and with cumulation, respectively. Section 5 presents an analysis of the variance of the logarithm of the stepsize and in Section 6 we summarize our results.

Notations In this paper, we denote t the iteration or time index, n the search space dimension, N (0, 1) a standard normal distribution, i.e. a normal distribution with mean zero and standard deviation 1. The multivariate normal distribution with mean vector zero and covariance matrix identity will be denoted N (0, I n ), the i th order statistic of λ standard normal distributions N i:λ , and Ψ i:λ its distribution. If x = (x 1 , • • • , x n ) ∈ R n is a vector, then [x] i will be its value on the i th dimension, that is [x] i = x i . A random variable X distributed according to a law L will be denoted X ∼ L. If A is a subset of X , we will denote A c its complement in X .

The (1, λ)-CSA-ES

We denote with X t the parent at the t th iteration. From the parent point X t , λ children are generated: Y t,i = X t +σ t ξ t,i with i ∈ [ [1, λ]], and ξ t,i ∼ N (0, I n ), (ξ t,i ) i∈[ [1,λ]] i.i.d. Due to the (1, λ) selection scheme, from these children, the one minimizing the function f is selected: X t+1 = argmin{f (Y ), Y ∈ {Y t,1 , ..., Y t,λ }}. This latter equation implicitly defines the random variable ξ t as

X t+1 = X t + σ t ξ t . (1)
In order to adapt the step-size, the cumulative path is defined as

p t+1 = (1 -c)p t + c(2 -c) ξ t (2)
with 0 < c ≤ 1. The constant 1/c represents the life span of the information contained in p t , as after 1/c generations p t is multiplied by a factor that approaches 1/e ≈ 0.37 for c → 0 from below (indeed (1-c) 1/c ≤ exp(-1)). The typical value for c is between 1/ √ n and 1/n. We will consider that p 0 ∼ N (0, I n ) as it makes the algorithm easier to analyze.

The normalization constant c(2 -c) in front of ξ t in Eq. ( 2) is chosen so that under random selection and if p t is distributed according to N (0, I n ) then also p t+1 follows N (0, I n ). Hence the length of the path can be compared to the expected length of N (0, I n ) representing the expected length under random selection.

The step-size update rule increases the step-size if the length of the path is larger than the length under random selection and decreases it if the length is shorter than under random selection:

σ t+1 = σ t exp c d σ p t+1 E( N (0, I n ) ) -1
where the damping parameter d σ determines how much the step-size can change and is set to d σ = 1. A simplification of the update considers the squared length of the path [START_REF] Arnold | Random dynamics optimum tracking with evolution strategies[END_REF]:

σ t+1 = σ t exp c 2d σ p t+1 2 n -1 . ( 3
)
This rule is easier to analyse and we will use it throughout the paper. We will denote η t the random variable for the step-size change, i.e.

η t = exp(c/(2d σ )( p t+1 2 /n -1)), and for u ∈ R n , η (u) = exp(c/(2d σ )( u 2 /n -1)).
Preliminary results on linear functions. Selection on the linear function,

f (x) = [x] 1 , is determined by [X t ] 1 + σ t [ξ t ] 1 ≤ [X t ] 1 + σ t ξ t,i 1 for all i which is equivalent to [ξ t ] 1 ≤ ξ t,i 1
for all i where by definition ξ t,i 1 is distributed according to N (0, 1). Therefore the first coordinate of the selected step is distributed according to N 1:λ and all others coordinates are distributed according to N (0, 1), i.e. selection does not bias the distribution along the coordinates 2, . . . , n. Overall we have the following result.

Lemma 1. On the linear function f (x) = x 1 , the selected steps (ξ t ) t∈N of the (1, λ)-ES are i.i.d. and distributed according to the vector ξ := (N 1:λ , N 2 , . . . , N n ) where N i ∼ N (0, 1) for i ≥ 2.

Because the selected steps ξ t are i.i.d. the path defined in Eq. 2 is an autonomous Markov chain, that we will denote P = (p t ) t∈N . Note that if the distribution of the selected step depended on (X t , σ t ) as it is generally the case on non-linear functions, then the path alone would not be a Markov Chain, however (X t , σ t , p t ) would be an autonomous Markov Chain. In order to study whether the (1, λ)-CSA-ES diverges geometrically, we investigate the log of the step-size change, whose formula can be immediately deduced from Eq. 3:

ln σ t+1 σ t = c 2d σ p t+1 2 n -1 (4)
By summing up this equation from 0 to t -1 we obtain

1 t ln σ t σ 0 = c 2d σ 1 t t k=1 p k 2 n -1 . (5)
We are interested to know whether 1 t ln(σ t /σ 0 ) converges to a constant. In case this constant is positive this will prove that the (1, λ)-CSA-ES diverges geometrically. We recognize thanks to (5) that this quantity is equal to the sum of t terms divided by t that suggests the use of the law of large numbers to prove convergence of (5). We will start by investigating the case without cumulation c = 1 (Section 3) and then the case with cumulation (Section 4).

Divergence rate of (1, λ)-CSA-ES without cumulation

In this section we study the (1, λ)-CSA-ES without cumulation, i.e. c = 1. In this case, the path always equals to the selected step, i.e. for all t, we have p t+1 = ξ t . We have proven in Lemma 1 that ξ t are i.i.d. according to ξ. This allows us to use the standard law of large numbers to find the limit of 1 t ln(σ t /σ 0 ) as well as compute the expected log-step-size change.

Proposition 1. Let ∆ σ := 1 2dσn E N 2 1:λ -1 .
On linear functions, the (1, λ)-CSA-ES without cumulation satisfies (i) almost surely lim t→∞ 1 t ln (σ t /σ 0 ) = ∆ σ , and (ii) for all t ∈ N, E(ln(σ t+1 /σ t )) = ∆ σ .

Proof. We have identified in Lemma 1 that the first coordinate of ξ t is distributed according to N 1:λ and the other coordinates according to N (0, 1), hence

E ξ t 2 = E [ξ t ] 1 2 + n i=2 E [ξ t ] 2 i = E N 2 1:λ + n -1. Therefore E ξ t 2 /n -1 = (E N 2 1:λ -1)/n.
By applying this to Eq. ( 4), we deduce that E(ln

(σ t+1 /σ t ) = 1/(2d σ n)(E(N 2 1:λ ) -1). Furthermore, as E(N 2 1:λ ) ≤ E((λN (0, 1)) 2 ) = λ 2 < ∞, we have E( ξ t
2 ) < ∞. The sequence ( ξ t 2 ) t∈N being i.i.d according to Lemma 1, and being integrable as we just showed, we can apply the strong law of large numbers on Eq. ( 5). We obtain

1 t ln σ t σ 0 = 1 2d σ 1 t t-1 k=0 ξ k 2 n -1 a.s. -→ t→∞ 1 2d σ E ξ • 2 n -1 = 1 2d σ n E N 2 1:λ -1
The proposition reveals that the sign of E N 2 1:λ -1 determines whether the step-size diverges to infinity. In the following, we show that E N 2 1:λ increases in λ for λ ≥ 2 and that the (1, λ)-ES diverges for λ ≥ 3. For λ = 1 and λ = 2, the step-size follows a random walk on the log-scale. To prove this we need the following lemma: Lemma 2 ( [START_REF] Ostermeier | Sizing the population with respect to the local progress in (1,lambda)-evolution strategies -a theoretical analysis[END_REF]). Let g be a real valued function on R. For λ ≥ 2,

(λ + 1) E (g (N 1:λ )) = E (g (N 2:λ+1 )) + λE (g (N 1:λ+1 )) . (6)
Proof. of Lemma 2 This method can be found with more details in [START_REF] Ostermeier | Sizing the population with respect to the local progress in (1,lambda)-evolution strategies -a theoretical analysis[END_REF]. Let χ i = g (ξ i ), and χ i:λ = g (ξ i:λ ). Note that in general

χ 1:λ = min i∈[[1,λ]] χ i .
The sorting is made on (ξ i ), not on (χ i ).

We will also note χ {j} i:λ the i t h order statistic after that the variable χ j has been taken away. χ

[j] i:λ will be i t h order statistic after χ j:λ has been taken away : if i = 1 then we have χ From the first equation we deduce that λE (χ 1:λ-1 ) = λE(χ

{i} 1:λ ) = λ i=1 E(χ {i} 1:λ ) = E( λ i=1 χ {i} 1:λ ).
With the second equation, we get that E(

λ i=1 χ {i} 1:λ ) = E( λ i=1 χ [i] 1λ ) = E(χ 2:λ ) + (λ -1)E(χ 1:λ ).
By combining both, we get the final equation:

(λ -1)(E(χ 1:λ ) -E(χ 1:λ-1 )) = E(χ 1:λ-1 ) -E(χ 2:λ )
We are now ready to prove the following result. [1,λ]] be independent random variables, distributed according to N (0, 1), and N i:λ the i th order statistic of

Lemma 3. Let (N i ) i∈[
(N i ) i∈[[1,λ]] . Then E N 2 1:1 = E N 2 1:2 = 1. In addition, for all λ ≥ 2, E N 2 1:λ+1 > E N 2 1:λ .
Proof. of Lemma 3 The strict monotony of E(N 2 1:λ ) in λ from the previous proposition is equivalent to show that

E(N 2 1:λ ) > E N 2 2:λ for λ ≥ 3. Indeed E(N 2 1:λ ) - E(N 2 1:λ-1 ) = E(N 2 1:λ -N 2
2:λ )/λ which follows from Lemma 2 taking g as the square function.

Let

E 1 = {ω ∈ Ω|N 2 1:λ (ω) < N 2 2:λ (ω)}, where Ω = R λ and P (ω) = exp(-ω 2 /2)/ √ 2π.
For ω ∈ Ω, let us note ω i:λ the i th order statistic of ([ω] j ) j∈[ [1,λ]] . Let g be a function that maps ω ∈ Ω to ω ∈ Ω, where ω1:λ = -ω 2:λ , ω2:λ = ω 1:λ and for i ≥ 3, ωi:λ = ω i:λ . The function g is bijective between E 1 and its image by g, E 2 . Let us note that for ω ∈ E 1 , N 2 2:λ (ω) -N 2 1:λ (ω) = N 2 1:λ (g(ω)) -N 2 2:λ (g(ω)), and P (ω) = P (g(ω)) since the standard normal distribution is symmetric. That is E1 (N 2 2:λ (ω) -N 2 1:λ (ω))P (ω)dω = E2 (N 2 1:λ (ω)-N 2 2:λ (ω))P (ω)dω by a change of variables õmega = g(ω). As according to the definition of E 1 , for all ω ∈ Ω\E 1 N 2 1:λ (ω) ≥ N 2 2:λ (ω), and that E 1 is properly counterweighted by E 2 in the expected value of N 2 1:λ -N 2 2:λ , we do have

E(N 2 1:λ ) ≥ E(N 2 2:λ ) for all λ ≥ 2. For λ ≥ 3, let E 3 = {ω ∈ Ω|ω 3:λ ∈] -|ω 1:λ |, |ω 1:λ |[ and ω 1:λ < ω 2:λ }. Then, for ω ∈ E 3 we also have ω 2:λ ∈] -|ω 1:λ |, |ω 1:λ |[, so N 2 1:λ (ω) > N 2 2:λ (ω) which means ω / ∈ E 1 , or E 1 ∩ E 3 = ∅. For ω ∈ E 1 , as ω 2 1:λ < ω 2 2:λ and ω 1:λ < ω2 : λ, ω 2:λ > 0, so ω 3:λ / ∈ [-ω2 : λ, ω2 : λ]. Hence, as g(ω) 3:λ = ω3 : λ and [-ω2 : λ, ω2 : λ] = [-|g(ω) 1:λ )|, |g(ω) 1:λ )|], g(ω) / ∈ E 3 . That is E 2 ∩ E 3 = ∅. So E 3 is disjoint with E 1 and E 2 .
Furthermore, for every ω ∈ Ω, except when ω 1:λ = 0 which is a negligible subset of events, there exists a non negligible set of (ω

i:λ ) i∈[[1,λ]] such that ω 3:λ ∈ ] -|ω 1:λ |, |ω 1:λ |[ and ω 1:λ < ω 2:λ }. So E 3 is a non negligible subset of Ω, where N 2 1:λ (ω) > N 2 2:λ (ω). Hence E(N 2 1:λ (ω)) > E(N 2 2:λ (ω)), which is the monotony of Lemma 3. For λ = 1, N 1:1 ∼ N (0, 1) so E(N 2 1:1 ) = 1. For λ = 2 we have E(N 2 1:2 + N 2 2:2 ) = 2E(N (0, 1) 2 ) = 2, and since the normal distribution is symmetric E(N 2 1:2 ) = E(N 2 2:2 ), hence E(N 2 1:2 ) = 1.
We can now link Proposition 1 and Lemma 3 into the following theorem:

Theorem 1. On linear functions, for λ ≥ 3, the step-size of the (1, λ)-CSA-ES without cumulation (c = 1) diverges geometrically almost surely and in expectation at the rate

1/(2d σ n)(E(N 2 1:λ ) -1), i.e. 1 t ln σ t σ 0 a.s. -→ t→∞ E ln σ t+1 σ t = 1 2d σ n E N 2 1:λ -1 . (7)
For λ = 1 and λ = 2, without cumulation, the logarithm of the step-size does an additive unbiased random walk i.e. ln σ t+1 = ln

σ t + W t where E[W t ] = 0. More precisely W t ∼ 1/(2d σ )(χ 2 n /n -1) for λ = 1, and W t ∼ 1/(2d σ )((N 2 1:2 + χ 2 n-1 )/n - 1) for λ = 2, where χ 2
k stands for the chi-squared distribution with k degree of freedom.

Proof. For λ > 2, from Lemma 3 we know that 7) is strictly positive, and with Proposition 1 we get that the step-size diverges geometrically almost surely at the rate 1/(2d σ )(E(N 2 1:λ ) -1). With Eq. 4 we have ln(σ t+1 ) = ln(σ t ) + W t , with W t = 1/(2d σ )( ξ t 2 /n -1). For λ = 1 and λ = 2, according to Lemma 3, E(W t ) = 0. Hence ln(σ t ) does an additive unbiased random walk. Furthermore

E(N 2 1:λ ) > E(N 2 1:2 ) = 1. Therefore E(N 2 1:λ ) -1 > 0, hence Eq. (
ξ 2 = N 2 1:λ + χ 2 n-1 , so for λ = 1, since N 1:1 = N (0, 1), ξ 2 = χ 2 n . 3.1 Geometric divergence of ([X t ] 1 ) t∈N
As the selection occurs only on the first dimension, if there is geometric divergence for

X t , it is on [X t ] 1 . From Eq (1) ln [X t+1 ] 1 [X t ] 1 = ln 1 + σ t [X t ] 1 [ξ t ] 1 .
Summing previous equation from 0 till t -1 and dividing by t gives us that

1 t ln [X t ] 1 [X 0 ] 1 = 1 t t-1 k=0 ln 1 + σ k [X k ] 1 [ξ t ] 1 . (8)
Although it is not obvious at first sight, it is important to take the logarithm, as we intuitively know that the speed of σ t and the speed of X t are connected. The divergence rate of σ t being log-linear, so should be the one of X t . Let Z -1 = 0, and

Z t = [Xt+1] 1 -[X0] 1 σt for t ≥ 0. Z t+1 = [X t+2 ] 1 -[X 0 ] 1 σ t+1 = [X t+1 ] 1 -[X 0 ] 1 + σ t+1 ξ t+1 1 σ t+1 Z t+1 = Z t η t + ξ t+1 1 using that σ t+1 = σ t η t . According to Lemma 1 (ξ t ) t∈N is independent over time. As η t = exp(( ξ t 2 /n -1)/(2d σ )), (η t ) t∈N is also independent over time. Therefore, Z = (Z t ) t∈N , is a Markov chain.
By introducing Z in Eq (8), we obtain:

1 t ln [X t ] 1 [X 0 ] 1 = 1 t t-1 k=0 ln 1 + σ k-1 η k-1 [X k ] 1 [ξ k ] 1 = 1 t t-1 k=0 ln 1 + η k-1 Z k-1 [ξ k ] 1 = 1 t t-1 k=0 ln Z k-1 η k-1 + ξ k Z k-1 η k-1 = 1 t t-1 k=0 ln |Z k | -ln |Z k-1 | + ln η k-1 (9)
The right hand side of this equation reminds us again of the law of large numbers. There is no independence over time, but Z being a Markov chain, if it follows some specific stability properties of Markov chains, then a law of large numbers may apply.

Study of the Markov chain Z

To apply a law of large numbers to a Markov chain, it has to satisfies some stability properties: in particular, the Markov chain P has to be ϕ-irreducible, that is, there exists a measure ϕ such that every Borel set A of R n with ϕ(A) > 0 has a positive probability to be reached in a finite number of steps by P starting from any p 0 ∈ R n . In addition, the chain P needs to be (i) positive, that is the chain admits an invariant probability measure π, i.e., for any borelian A, π(A) = R n P (x, A)π(dx) with P (x, A) being the probability to transition in one time step from x into A, and (ii) Harris recurrent which means for any borelian A such that ϕ(A) > 0, the chain P visits A an infinite number of times with probability one. Under those conditions, P satisfies a law of large numbers, more precisely: Lemma 4. [11, 17.0.1] Suppose that Φ is a positive Harris chain with stationary measure π, and let g be a π-integrable function that is such that π

(|g|) = R n |g(x)|π(dx) < ∞. Then 1/t t k=1 g(Φ k ) a.s -→ t→∞ π(g) . (10)
To show that a Markov defined in a space X is positive Harris recurrent, we generally show that the chain follows a so-called drift condition over a small set, that is for a function V , an inequality over the drift operator ∆V : x → X V (y)P (x, dy) -V (x). A small set is a borel set such that there exists a m ∈ N * and a non-trivial measure ν m on β(X) such that for all x ∈ C, B ∈ β(X), P m (x, B) ≥ ν m (B). The set C is then called a ν m -small set. The chain also needs to be aperiodic, that is there is no

d-cycle, that is disjoint Borel sets (D i ) i∈[[1,d]] , such that for x ∈ D i , P (x, D i+1 ) = 1 for i = 0 • • • d -1(modd), and [∪ d i=1 ] c is ϕ-negligible.
If there exists a ν 1 -small-set A such that ν 1 (A) > 0, then the chain is strongly aperiodic (and therefore aperiodic). We then have the following lemma.

Lemma 5. [11, 14.0.1] Suppose that the chain Φ is ϕ-irreductible and aperiodic, and f ≥ 1 a function on X. Let us assume that there exists V some extended-valued nonnegative function finite for some x 0 ∈ X, a small set C and b ∈ R such that

∆V (x) ≤ -f (x) + b1 C (x) , x ∈ X. ( 11
)
Then the chain Φ is positive Harris recurrent with invariant probability measure π and

π(f ) = X π(dx)f (x) < ∞ . ( 12
)
To prove the irreducibility, aperiodicity and to exhibit the small sets of the Markov chain Z through its transition kernel would be difficult. Instead, it can be done by showing some properties of its underlying control model. In our case, the model associated to Z is called a non-linear state space model. We will, in the following, define this non-linear state space model and some of its properties.

Suppose We define its associated control model CM(F ) the deterministic system

X = {X k }, X k ∈ X . If there is a smooth function (C ∞ ) F such that X k+1 = F (X k , W k+1 ) with (W i ) i∈N
x k = F k (x 0 , u 1 , • • • , u k ), where F k is given by F k (x 0 , u 1 , • • • , u k ) = F (F k-1 (x 0 , u 1 , • • • , u k-1 ), u k ), and F 0 (x 0 ) = x 0 , provided that (u i ) i∈N lies in the control set O w .
For a point x ∈ X , and k ∈ N we define

A k + (x) = {F k (x, u 1 , • • • , u k )|u i ∈ O w ∀i ∈ N},
the set of points reachable from x after k steps of time. And A + (x) = i∈N A i + (x). The associated control model CM(F ) is called forward accessible if for each x 0 ∈ X , the set A + (x 0 ) has non empty-interior.

Let E be a subset of X . We note A + (E) = x∈E A + (x), and we say that E is invariant if A + (E) ⊂ E. We call a set minimal if it is closed, invariant, and does not strictly contain any closed and invariant subset. Restricted to a minimal set, a Markov chain has strong properties, as stated in the following lemma. Furthermore, if the control set O w and M are connected, and that M is the unique minimal set of the CM(F ), then the NSS(F ) model is a ψ-irreducible aperiodic T-chain for which every compact set is a small set.

We can now prove the following lemma: Proof. This is exactly the result of Theorem 6 when all conditions are fulfilled. We then have to show the right properties of the underlying control model.

If we note

F (X k , W k+1 ) = X k exp -1/2d σ W k+1 2 /n -1 + [W k+1 ] 1 , then we do have Z t+1 = F (Z t , ξ t ).
The function F is smooth (it is not smooth along the instances ξ t,i , but along the chosen step). Furthermore, the distribution of ξ t admits a continuous density, whose support is R n . Therefore the process Z is a NSS(F ) model of control set R n .

We now have to show that the associated control model is forward accessible. Let z ∈ R. When [ξ t ] 1 → ±∞, F (z, ξ t ) → ±∞. As F is continuous, for the right value of [ξ t ] 1 any point of R can be reach. Therefore for any z ∈ R, A + (z) = R. The set R has a non-empty interior, so the CM(F ) is forward accessible.

As from any point of R, all of R can be reached, the only invariant set is R itself. It is therefore the only minimal set. Finally, the control set O w = R n is connected, and so is the only minimal set, so all the conditions of Lemma 6 are met. So the Markov chain Z is ψ-irreducible, aperiodic, and compacts of R are small-sets.

We may now show Foster-Lyapunov drift conditions to ensure the Harris positive recurrence on the chain Z. In order to do so, we will need the following Lemma:

Lemma 8. Let exp(-1 2dσ ( ξ 2 n - 1
)) be denoted η • . For all λ > 2 there exists α > 0 such that

E η • -α -1 < 0 . (13)
Proof. Using the Taylor series of the exponential function we have

E η • -α = E exp - α 2d σ ξ 2 n -1 = E    ∞ i=0 -α 2dσ ξ 2 n -1 i i!    = 1 -α 1 2d σ n E N 2 1:λ -1 -o α 2 .
According to Lemma 3 E N 2 1:λ > 1 for λ > 2, so when α goes to 0 we have

E η • -α < 1.
We are now ready to prove the following lemma:

Lemma 9. The Markov chain Z is Harris recurrent positive, and admits a unique invariant measure µ such that for f :

x → |x| α ∈ R, µ(f ) = R µ(dx)f (x) < ∞
, with α such that Eq. (13) holds true.

Proof. By using Lemma 7 and Lemma 5, we just need the drift condition [START_REF] Meyn | Markov chains and stochastic stability[END_REF] to prove Lemma 9. Let V be such that for x ∈ R, V (x) = |x| α + 1.

∆V (x) = R P (x, dy)V (y) -V (x) = R P x η • + [ξ • ] 1 ∈ dy (1 + |y| α ) -(1 + |x| α ) = E x η • + [ξ • ] 1 α -|x| α ≤ |x| α E η • -α -1 + E [ξ • ] α 1 ∆V (x) V (x) = |x| α 1 + |x| α E η • -α -1 + 1 1 + |x| α E [ξ t ] α 1 lim |x|-→∞ ∆V (x) V (x) = E η • -α -1
We take α such that Eq. ( 13) holds true (as according to Lemma 8, there exists such a α). As E(η • -α -1) < 0, there exists > 0 and M > 0 such that for

all |x| ≥ M , ∆V /V (x) ≤ -. Let b be equal to E([ξ t ] 2 1 ) + V (M ). Then for all |x| ≤ M , ∆V (x) ≤ -V (x) + b. Therefore, if we note C = [-M, M ],
which is according to Lemma 7 a small-set, we do have ∆V (x) ≤ -V (x) + b1 C (x) which is Eq. ( 11) with f = V . Therefore from Lemma 5 the chain Z is positive Harris recurrent with invariant probability measure µ, and

V is µ-integrable. As R µ(dx)|x| α = 1/ R µ(dx) V (x) -1 < ∞, the function x → |x| α is also µ-integrable.
In order to use Lemma 4 on Z with the function g : x → E (ln |x/x -[ξ • ] 1 |), we must prove that this function is µ-integrable, that is R g(u)µ(du) < ∞. To do so we will need the following lemma on the existence of moments for stationary Markov chains: Lemma 10. Let Z be a Harris-recurrent Markov chain with stationary measure µ, on a state space (S, F), with F is σ-field of subsets of S. Let f be a positive measurable function on S.

In order that S f (z)µ(dz) < ∞, it suffices that for some set A ∈ F such that 0 < µ(A) and A f (z)µ(dz) < ∞, and some measurable function g with g(z) ≥ f (z) for z ∈ A c , 1. We may now prove the following theorem:

A c P (z, dy)g(y) ≤ g(z) -f (z) , ∀x ∈ A c

sup

Theorem 2. On linear functions, for λ ≥ 3, the absolute value of the first dimension of the parent point in the (1, λ)-CSA-ES without cumulation (c = 1) diverges geometrically almost surely at the rate of 1/(2d σ n)E(N 2 1:λ -1), i.e.

1 t ln [X t ] 1 [X 0 ] 1 a.s -→ t→∞ 1 2d σ n E N 2 1:λ -1 . (14)
Proof. We will first prove here that the function g : x → ln |x| is µ-integrable. From Lemma 9 we know that the function f : x → |x| α is µ-integrable, and as for any M > 0, and any x ∈ [-M, M ] c there exists

K > 0 such that K|x| α > | ln |x||, then g1 A c is µ-integrable, with A = [-M, M ].
So what is left is to prove that g1 A is also µ-integrable. We will now check the conditions to use Lemma 10.

According to Lemma 7 the chain Z is open-set irreducible, so µ(A c ) > 0. For C > 0, if we take h : z → C/ |z|, with M small enough we do have for all z ∈ A c , h(z) ≥ |g(z)|. Furthermore, if we study the inequality

A c P (z, dy)h(y) ≤ h(z) -|g(z)| S P z η • + [ξ • ] 1 ∈ dy 1 A c (y) C |y| ≤ C |z| -|ln|z|| E     1 z η • + [ξ • ] 1 1 A c z η • + [ξ • ] 1     ≤ 1 |z| - | ln |z|| C
We can increase C up until | ln |z||/C is negligible compared to 1/ |z|, and we can decrease M to make

E 1/ |z/η • + [ξ • ] 1 |1 A c (z/η • + [ξ • ] 1 )
as small as we would like it to be, as it decreases the size of A c , so the inequality holds if we choose M and 1/C small enough. The second inequality for Lemma 10 holds as well:

A c P (u, dv)h(v) ≤ A c C |v| dv = 4C √ M < ∞
Finally, according to Lemma 9, the chain Z is Harris recurrent. So Lemma 10 shows that g is µ-integrable. This allows us to apply Lemma 4 to the function g:

1/t t k=1 g(z k ) a.s -→ t→∞ µ(g).
With Lemma 1 we can apply a strong law of large numbers upon 1/t

t-1 k=0 ln |η k-1 | = 1/t t-1 k=0 1/(2d σ )(ξ k-1 /n-1), to get as in the proof of Proposition 1 1/(2d σ n)(E(N 2 1:λ )- 1).
By inserting these results into Eq. ( 9), we get that

1/t ln | [X t ] 1 / [X 0 ] 1 | a.s -→ t→∞ µ(g) -µ(g) + 1/(2d σ n)(E(N 2 1:λ ) -1)
, which with Lemma 3 is strictly positive for λ ≥ 3.

Divergence rate of CSA-ES with cumulation

We are now investigating the (1, λ)-CSA-ES with cumulation, i.e. 0 < c < 1.

According to Lemma 1, the random variables (ξ t ) t∈N are i.i.d., hence the path P = (p t ) t∈N is a Markov chain. By a recurrence on Eq. ( 2) we see that the path follows the following equation

p t = (1 -c) t p 0 + c(2 -c) t-1 k=0 (1 -c) k ξ t-1-k i.i.d. . (15)
For i = 1, [ξ t ] i ∼ N (0, 1) and, as also [p 0 ] i ∼ N (0, 1), by recurrence [p t ] i ∼ N (0, 1) for all t ∈ N. For i = 1 with cumulation (c < 1), the influence of [p 0 ] 1 vanishes with (1 -c) t . Furthermore, as from Lemma 1 the sequence ([ξ t ] 1 ]) t∈N is independent, we get by applying the Kolgomorov's three series theorem that the series

t-1 k=0 (1 -c) k ξ t-1-k 1 converges almost surely.
Therefore, the first component of the path becomes distributed as the random variable

[p ∞ ] 1 = c(2 -c) ∞ k=0 (1 - c) k [ξ k ] 1 (by re-indexing the variable ξ t-1-k in ξ k , as the sequence (ξ t ) t∈N is i.i.d.).
As in Subsection 3.1 we will show that P has the right stability properties to apply a law of large numbers to it. First we will extract from P the part of interest as stated in the following lemma.

Lemma 11. On linear functions, for any λ the step-size of the (1, λ)-CSA-ES follows almost surely

1 t ln σ t σ 0 - c 2d σ n 1 t t i=1 ([p i ] 2 1 -1 a.s -→ t→∞ 0 , (16)
and in expectancy

E ln σ t+1 σ t = c 2d σ n E [p t+1 ] 2 1 -1 (17)
Proof. We separate Eq. ( 5) over the dimensions, which gives us that 1/t ln( 16).

σ t /σ 0 ) = c/(2d σ n)( n i=1 1/t t j=1 [p j ] 2 i -n), so 1/t ln(σ t /σ 0 ) -c/(2d σ n)(1/t t j=1 [p j ] 2 1 - 1) = n i=2 c/(2d σ n)(1/t t j=1 [p j ] 2 i -1). As for i = 1, [p 0 ] i N (0, I n ) and [ξ 0 ] i has no selection pressure, then [p 1 ] i N (0, I n ),
The same reasoning over Eq. ( 4) gives Eq. ( 17).

The part of P left to analyse is its first dimension

[P] 1 = ([p i ] 1 ) i∈N .
We start the study of [P] 1 with the following lemma.

Lemma 12. The Markov chain [P] 1 is ϕ-irreducible, aperiodic, and compacts of R are small-sets.

Proof. We have the following transition kernel:

P (p, A) = R 1 A (1 -c)p + c(2 -c)u P (N 1:λ = u)du .
With a change of variables ũ = (1 -c)p c(2 -c)u, we get that

P (p, A) = 1 (2 -c)c R 1 A (ũ) P N 1:λ = ũ -(1 -c)p
(2 -c)c dũ .

As P (N 1:λ = x) > 0 for all x ∈ R, for all A non-µ Leb -negligible we have P (p, A) > 0, thus the chain [P] 1 is µ Leb -irreducible. Furthermore, if we take C a non-µ Leb -negligible compact of R, and ν C a measure such that for A a borel set of R,

ν C (A) = 1 √ (2-c)c R 1 A (ũ) min p∈C P ξ t = (ũ -(1 -c)p) / (2 -c)c dũ, we see that P (p, A) ≥ ν C (A)
for all p ∈ R, while ν C is not a trivial measure (indeed, P (N 1:λ = x) > k > 0 for all x ∈ C); Therefore compact sets of R are small sets for [P] 1 . Finally, ν C (C) > 0, so the chain [P] 1 is strongly aperiodic.

We use this new lemma with Lemma 5 to prove what is needed to apply the law of large numbers on [P] 1 .

Lemma 13. The chain [P] 1 is Harris recurrent positive with invariant measure µ path , and the function x → x 2 is µ path -integrable.

Proof. We now have to get the right drift condition for the chain. Let V : x → x 2 + 1.

∆V (x) = R V (y)P (x, dy) -V (x) ∆V (x) = R y 2 + 1 P (1 -c)x + c(2 -c) [ξ • ] 1 ∈ dy -x 2 + 1 ∆V (x) = E 1 -c)x + c(2 -c) [ξ • ] 1 2 + 1 -x 2 -1 ∆V (x) ≤ ((1 -c) 2 -1)x 2 + 2|x| c(2 -c)E ([ξ • ] 1 ) + c(2 -c)E [ξ • ] 1 2 ∆V (x) V (x) ≤ -c(2 -c) x 2 1 + x 2 + 2|x| c(2 -c) 1 + x 2 E (|[ξ • ] 1 |) + c(2 -c) 1 + x 2 E [ξ • ] 1 2 lim |x|→∞ ∆V (x) V (x) ≤ -c(2 -c) As 0 < c ≤ 1, c(2 -c
) is strictly positive and therefore, for > 0 there exists

C = [-M, M ] with M > 0 such that for all x ∈ C c , ∆V (x)/V (x) ≤ -. If we take b = V (M ) + 2M c(2 -c)E(| [ξ • ] 1 |) + c(2 -c)E([ξ t ] 1 2 ), then for all x ∈ C we have ∆V (x) ≤ b. Hence the drift condition ∆V (x) ≤ -V (x) + b1 C is satisfied for all x ∈ R.
According to Lemma 12 the chain [P] 1 is ϕ-irreducible and aperiodic, so with Lemma 5 it is positive Harris recurrent, with invariant measure µ path , and V is µ pathintegrable. Therefore the function x → x 2 is also µ path -integrable.

To obtain an equality between the rate we get through almost sure divergence, and the rate in expectation, we need to define the f -norm, which for a signed measure ν and a function f ≥ 1 is equal to ν f = sup g:|g|≤f |ν(g)|, and we need the following lemma.

Lemma 14. [11, 14.3.5] Suppose Φ is an aperiodic positive Harris chain on a space X with stationary measure π, and that there exists some non-negative function V , a function f ≥ 1, a small-set C and b ∈ R such that for all x ∈ X , ∆V (x) ≤ -f (x) + b1 C (x). Then for all initial probability distribution ν, νP n -π f -→ t→∞ 0.

We now obtain geometric divergence of the step-size and get an explicit estimate of the expression of the divergence rate.

Theorem 3. The step-size of the (1, λ)-CSA-ES with λ ≥ 2 diverges geometrically fast if c < 1 or λ ≥ 3. Almost surely and in expectation we have for 0 < c ≤ 1,

1 t ln σ t σ 0 -→ t→∞ 1 2d σ n 2(1 -c) E (N 1:λ ) 2 + c E N 2 1:λ -1
>0 for λ≥3 and for λ=2 and c<1

.

Proof. We will start by the convergence in expectation. From Eq. ( 17) we see that the part to develop is

E( p t+1 2 1
). By recurrence

p t+1 1 = (1-c) t+1 [p 0 ] 1 + c(2 -c) t i=0 (1- c) i ξ t-i 1 .
When t goes to infinity, the influence of [p 0 ] 1 in this equation goes to 0 with (1 -c) t+1 , so we can remove it when taking the limit:

lim t→∞ E p t+1 2 1 = lim t→∞ E c(2 -c) t i=0 (1 -c) i ξ t-i 1 2 (19)
We will now develop the sum with the square, such that we have either a product ξ t-i 1 ξ t-j 1 with i = j, or ξ t-j

. This way, we can separate the variables by using Lemma 1 with the independence of ξ i over time. To do so, we use the development formula

( n i=1 a n ) 2 = 2 n i=1 n j=i+1 a i a j + n i=1 a 2 i .
We take the limit of

E( p t+1 2
1
) and find that it is equal to

lim t→∞ c(2-c)      2 t i=0 t j=i+1 (1-c) i+j E ξ t-i 1 ξ t-j 1 =E[ξ t-i ] 1 E[ξ t-j ] 1 =E[N 1:λ ] 2 + t i=0 (1-c) 2i E ξ t-i 2 1 =E[N 2 1:λ ]     
(20) Now the expected value does not depend on i or j, so what is left is to calculate

t i=0 t j=i+1 (1 -c) i+j and t i=0 (1 -c) 2i . We have t i=0 t j=i+1 (1 -c) i+j = t i=0 (1 -c) 2i+1 1-(1-c) t-i 1-(1-c
) and when we separates this sum in two, the right hand side goes to 0 for t → ∞. Therefore, the left hand side converges to

lim t→∞ t i=0 (1 - c) 2i+1 /c, which is equal to lim t→∞ (1 -c)/c t i=0 (1 -c) 2i . And t i=0 (1 -c) 2i is equal to (1 -(1 -c) 2t+2 )/(1 -(1 -c) 2 ), which converges to 1/(c(2 -c)). So, by inserting this in Eq. (20) we get that E p t+1 2 1 -→ t→∞ 2 1-c c E (N 1:λ ) 2 + E N 2 1:
λ , which gives us the right hand side of Eq. (18).

By summing E(ln(σ i+1 /σ i )) for i = 0, . . . , t -1 and dividing by t we have the Cesaro mean 1/tE(ln(σ t /σ 0 )) that converges to the same value that E(ln(σ t+1 /σ t )) converges to when t goes to infinity. Therefore we have in expectation Eq. (18).

We will now focus on the almost sure convergence. From Lemma 13, we see that we have the right conditions to apply Lemma 4 to the chain [P] 1 with the µ path -integrable function g 16) we obtain that

: x → x 2 . So 1/t t k=1 [p k ] 2 1 a.s -→ t→∞ µ path (g). With Eq. (
1/t ln(σ t /σ 0 ) a.s -→ t→∞ c/(2d σ n)(µ path (g) -1).
We will now prove that µ path (g) = lim t→∞ E( p t+1 2 1 ). Let ν be the initial distribution of [p 0 ]1, so we have

|E( p t+1 2 1 ) -µ path (g)| ≤ νP t+1 -µ path h , with h : x → 1 + x 2 .
From the proof of Lemma 13 and from Lemma 12 we have all conditions for Lemma 14. Therefore νP t+1 -µ path h -→ t→∞ 0, which shows that

µ path (g) = lim t→∞ E( p t+1 2 1 ) = (2 -2c)/cE(N 1:λ ) 2 + E(N 2 1:λ ). According to Lemma 3, for λ = 2, E(N 2 1:2 ) = 1
, so the RHS of Eq. ( 18) is equal to

(1 -c)/(d σ n)E(N 1:2 ) 2 .
The expected value of N 1:2 is strictly negative, so the previous expression is strictly positive. Furthermore, according to Lemma 3, E(N 2 1:λ ) increases strictly with λ, as does E(N 1:2 ) 2 . Therefore we have geometric divergence for λ ≥ 2 if c < 1, and for λ ≥ 3.

From Eq. ( 1) we see that the behaviour of the step-size and of (X t ) t∈N are directly related. Geometric divergence of the step-size, as shown in Theorem 3, means that also the movements in search space and the improvements on affine linear functions f increase geometrically fast. Analyzing (X t ) t∈N with cumulation would require to study a double Markov chain, which is left to possible future research.

5 Study of the variations of ln (σ t+1 /σ t )

The proof of Theorem 3 shows that the step size increase converges to the right hand side of Eq. ( 18), for t → ∞. When the dimension increases this increment goes to zero, which also suggests that it becomes more likely that σ t+1 is smaller than σ t . To analyze this behavior, we study the variance of ln (σ t+1 /σ t ) as a function of c and the dimension.

Theorem 4. The variance of ln (σ t+1 /σ t ) equals to Var ln

σ t+1 σ t = c 2 4d 2 σ n 2 E p t+1 4 1 -E p t+1 2 1 2 + 2(n -1) . (21) Furthermore, E p t+1 2 1 -→ t→∞ E N 2 1:λ + 2-2c c E (N 1:λ ) 2 and with a = 1 -c lim t→∞ E p t+1 4 1 = (1 -a 2 ) 2 1 -a 4 (k 4 + k 31 + k 22 + k 211 + k 1111 ) , (22)
where

k 4 = E N 4 1:λ , k 31 = 4 a(1+a+2a 2 ) 1-a 3 E N 3 1:λ E (N 1:λ ), k 22 = 6 a 2 1-a 2 E N 2 1:λ 2 , k 211 = 12 a 3 (1+2a+3a 2 ) (1-a 2 )(1-a 3 ) E N 2 1:λ E(N 1:λ ) 2 and k 1111 = 24 a 6 (1-a)(1-a 2 )(1-a 3 ) E (N 1:λ ) 4 .
Proof.

Var ln

σ t+1 σ t = Var c 2d σ p t+1 2 n -1 = c 2 4d 2 σ n 2 Var p t+1 2 E( p t+1 4 )-E( p t+1 2 ) 2 (23) The first part of Var( p t+1 2 ), E( p t+1 4
), is equal to E((

n i=1 p t+1 2 i ) 2
). We develop it along the dimensions such that we can use the independence of [p t+1 ] i with [p t+1 ] j for i = j, to get E(2

n i=1 n j=i+1 p t+1 2 i p t+1 2 j + n i=1 p t+1 4 i ). For i = 1 p t+1 i is distributed according to a standard normal distribution, so E p t+1 2 i = 1 and E p t+1 4 i = 3. E p t+1 4 = 2 n i=1 n j=i+1 E p t+1 2 i E p t+1 2 j + n i=1 E p t+1 4 i =   2 n i=2 n j=i+1 1   + 2 n j=2 E p t+1 2 1 + n i=2 3 + E p t+1 4 1 = 2 n i=2 (n -i) + 2(n -1)E p t+1 2 1 + 3(n -1) + E p t+1 4 1 = E p t+1 4 1 + 2(n -1)E p t+1 2 1 + (n -1)(n + 1)
The other part left is E( p t+1 2 ) 2 , which we develop along the dimensions to get

E( n i=1 p t+1 2 i ) 2 = (E( p t+1 2 1
) + (n -1)) ) is the same than the one done in the proof of Theorem 3, that is 4 ). We neglect in the limit when t goes to ∞ the part with (1 -c) t [p 0 ] 1 , as it converges fast to 0. So

E( p t+1 2 1 ) = (2 -2c)/cE(N 1:λ ) 2 + E(N 2 1:λ ). We now develop E( p t+1 4 1 ). We have E( p t+1 4 1 ) = E(((1-c) t [p 0 ] 1 + c(2 -c) t i=0 (1-c) i ξ t-i 1 )
lim t→∞ E p t+1 4 1 = lim t→∞ E   c 2 (2 -c) 2 t i=0 (1 -c) i ξ t-i 1 4   . ( 24
)
To develop the RHS of Eq.( 24) we use the following formula: for

(a i ) i∈[[1,m]]
This formula will allow us to use the independence over time of

[ξ t ] 1 from Lemma 1, so that E([ξ i ] 3 1 ξ j 1 ) = E([ξ i ] 3
1 )E( ξ j 1 ) = E(N 3 1:λ )E(N 1:λ ) for i = j, and so on. We apply Eq (25) on Eq (22), with a = 1 -c.

lim t→∞ E p t+1 4 1 c 2 (2 -c) 2 = lim t→∞ t i=0 a 4i E N 4 1:λ + 4 t i=0 t j=0 j =i a 3i+j E N 3 1:λ E (N 1:λ ) + 6 t i=0 t j=i+1 a 2i+2j E N 2 1:λ 2 + 12 t i=0 t j=0 j =i t k=j+1 k =i a 2i+j+k E N 2 1:λ E (N 1:λ ) 2 + 24 t i=0 t j=i+1 t k=j+1 t l=k+1 a i+j+k+l E (N 1:λ ) 4 (26)
We now have to develop each term of Eq. (26).

t i=0 a 4i = 1 -a 4(t+1) 1 -a 4 lim t→∞ t i=0 a 4i = 1 1 -a 4 (27) t i=0 t j=0 j =i a 3i+j = t-1 i=0 t j=i+1 a 3i+j + t i=1 i-1 j=0 a 3i+j (28) t-1 i=0 t j=i+1 a 3i+j = t-1 i=0 a 4i+1 1 -a t-i 1 -a lim t→∞ t-1 i=0 t j=i+1 a 3i+j = lim t→∞ a 1 -a t-1 i=0 a 4i = a (1 -a)(1 -a 4 ) (29) t i=1 i-1 j=0 a 3i+j = t i=1 a 3i 1 -a i 1 -a = 1 1 -a a 3 1 -a 3t 1 -a 3 -a 4 1 -a 4t 1 -a 4 lim t→∞ t i=1 i-1 j=0 a 3i+j = 1 1 -a a 3 1 -a 3 - a 4 1 -a 4 = a 3 (1 -a 4 ) -a 4 (1 -a 3 ) (1 -a)(1 -a 3 )(1 -a 4 ) = a 3 -a 4 (1 -a)(1 -a 3 )(1 -a 4 ) (30)
By combining Eq (29) with Eq (30) to Eq (28) we get

lim t→∞ t i=0 t j=0 j =i a 3i+j = a(1 -a 3 ) + a 3 -a 4 (1 -a)(1 -a 3 )(1 -a 4 ) = a(1 + a 2 -2a 3 ) (1 -a)(1 -a 3 )(1 -a 4 ) = a(1 -a)(1 + a + 2a 2 )) (1 -a)(1 -a 3 )(1 -a 4 ) = a(1 + a + 2a 2 )) (1 -a 3 )(1 -a 4 ) (31) t-1 i=0 t j=i+1 a 2i+2j = t-1 i=0 a 4i+2 1 -a 2(t-i) 1 -a 2 lim t→∞ t-1 i=0 t j=i+1 a 2i+2j = a 2 1 -a 2 t-1 i=0 a 4i = a 2 (1 -a 2 )(1 -a 4 ) (32) t i=0 t-1 j=0 j =i t k=j+1 k =i a 2i+j+k = t i=2 i-2 j=0 i-1 k=j+1 a 2i+j+k + t-1 i=1 i-1 j=0 t k=i+1 a 2i+j+k + t-2 i=0 t-1 j=i+1 t k=j+1 a 2i+j+k (33) t i=2 i-2 j=0 i-1 k=j+1 a 2i+j+k = t i=2 i-2 j=0 a 2i+2j+1 1 -a i-j-1 1 -a = 1 1 -a t i=2 a 2i+1 1 -a 2(i-1) 1 -a 2 -a 3i 1 -a i-1 1 -a = 1 1 -a a 5 1 -a 2 1 -a 2(t-1) 1 -a 2 - a 7 (1 -a 2 ) 1 -a 4(t-1) 1 -a 4 - a 6 1 -a 1 -a 3(t+1) 1 -a 3 + a 7 1 -a 1 -a 4(t+1) 1 -a 4 -→ t→∞ a 5 1 -a 1 (1 -a 2 ) 2 - a 2 (1 -a 2 )(1 -a 4 ) - a (1 -a)(1 -a 3 ) + a 2 (1 + a) (1 + a)(1 -a)(1 -a 4 ) -→ t→∞ a 5 1 -a (1 + a 2 ) (1 -a 2 ) 2 (1 + a 2 ) + a 3 (1 -a 2 )(1 -a 4 ) - a (1 -a)(1 -a 3 ) -→ t→∞ a 5 1 -a 1 + a 2 + a 3 (1 + a)(1 -a)(1 -a 4 ) - a (1 -a)(1 -a 3 ) -→ t→∞ a 5 (1 -a) 2 (1 + a 2 + a 3 )(1 -a 3 ) -a(1 + a)(1 -a 4 ))) (1 + a)(1 -a 3 )(1 -a 4 ) -→ t→∞ a 5 1 + a 2 -a 5 -a 6 -(a + a 2 -a 5 -a 6 ) (1 -a)(1 -a 2 )(1 -a 3 )(1 -a 4 ) -→ t→∞ a 5 (1 -a 2 )(1 -a 3 )(1 -a 4 ) (34) t-1 i=1 i-1 j=0 t k=i+1 a 2i+j+k = t-1 i=1 i-1 j=0 a 3i+j+1 1 -a t-i 1 -a -→ t→∞ lim t→∞ a 1 -a t-1 i=1 a 3i 1 -a i 1 -a -→ t→∞ lim t→∞ a (1 -a) 2 a 3 1 -a 3t 1 -a 3 -a 4 1 -a 4t 1 -a 4 -→ t→∞ a (1 -a) 2 a 3 (1 -a 4 ) -a 4 (1 -a 3 ) (1 -a 3 )(1 -a 4 ) -→ t→∞ a 4 -a 5 (1 -a) 2 (1 -a 3 )(1 -a 4 ) = a 4 (1 -a)(1 -a 3 )(1 -a 4 ) (35) t-2 i=0 t-1 j=i+1 t k=j+1 a 2i+j+k = t-2 i=0 t-1 j=i+1 a 2i+2j+1 1 -a t-j 1 -a -→ t→∞ lim t→∞ a 1 -a t-2 i=0 a 4i+2 1 -a 2(t-i-1) 1 -a 2 -→ t→∞ lim t→∞ a 3 (1 -a)(1 -a 2 ) 1 -a 4(t-1) 1 -a 4 -→ t→∞ a 3 (1 -a)(1 -a 2 )(1 -a 4 ) (36)
We now combine Eq (34), Eq. ( 35) and Eq. (34) in Eq. (33). Figure 1 shows the time evolution of ln(σ t /σ 0 ) for 5001 runs and c = 1 (left) and c = 1/ √ n (right). By comparing Figure 1a and Figure 1b we observe smaller variations of ln(σ t /σ 0 ) with the smaller value of c. Fig. 1: ln(σ t /σ 0 ) against t. The different curves represent the quantiles of a set of 5.10 3 + 1 samples, more precisely the 10 i -quantile and the 1 -10 -i -quantile for i from 1 to 4; and the median. We have n = 20 and λ = 8.

t i=0 t-1 j=0 j =i t k=j+1 k =i a 2i+j+k -→ t→∞ a 5 (1 -a) + a 4 (1 -a 2 ) + a 3 (1 -a 3 ) (1 -a)(1 -a 2 )(1 -a 3 )(1 -a 4 ) -→ t→∞ a 3 + a 4 + a 5 -3a 6 (1 -a)(1 -a 2 )(1 -a 3 )(1 -a 4 ) -→ t→∞ a 3 (1 + 2a + 3a 2 ) ((1 -a 2 )(1 -a 3 )(1 -a 4 ) (37)
the relative standard deviation. To get a value below one, c must be smaller for larger dimension. In agreement with Theorem 4, In Figure 2, right, the relative standard deviation increases like √ n with the dimension for constant c (three increasing curves). A careful study [START_REF] Chotard | Cumulative step-size adaptation on linear functions[END_REF] of the variance equation of Theorem 4 shows that for the choice of c = 1/(1 + n α ), if α > 1/3 the relative standard deviation converges to 0 with (n 2α + n)/n 3α . Taking α = 1/3 is a critical value where the relative standard deviation converges to 1/( √ 2E(N 1:λ ) 2 ). On the other hand, lower values of α makes the relative standard deviation diverge with n (1-3α)/2 .

Summary

We investigate throughout this paper the (1, λ)-CSA-ES on affine linear functions composed with strictly increasing transformations. We find, in Theorem 3, the limit distribution for ln(σ t /σ 0 )/t and rigorously prove the desired behaviour of σ with λ ≥ 3 for any c, and with λ = 2 and cumulation (0 < c < 1): the step-size diverges geometrically fast. In contrast, without cumulation (c = 1) and with λ = 2, a random walk on ln(σ) occurs, like for the (1, 2)-σSA-ES [START_REF] Hansen | An analysis of mutative σ-self-adaptation on linear fitness functions[END_REF] (and also for the same symmetry reason). We derive an expression for the variance of the step-size increment. On linear functions when c = 1/n α , for α ≥ 0 (α = 0 meaning c constant) and for n → ∞ the standard deviation is about (n 2α + n)/n 3α times larger than the step-size increment. From this follows that keeping c < 1/n 1/3 ensures that the standard deviation of ln(σ t+1 /σ t ) becomes negligible compared to ln(σ t+1 /σ t ) when the dimensions goes to infinity. That means, the signal to noise ratio goes to zero, giving the algorithm strong stability. The result confirms that even the largest default cumulation parameter c = 1/ √ n is a stable choice. The curves were plotted using Eq. ( 21) and Eq. ( 22). On the left, curves for (right to left) n = 2, 20, 200 and 2000. On the right, different curves for (top to bottom) c = 1, 0.5, 0.2, 1/(1 + n 1/4 ), 1/(1 + n 1/3 ), 1/(1 + n 1/2 ) and 1/(1 + n).

  = χ 1:λ , and for i = 1 χ

  being a sequence of i.i.d. random variables, whose marginal distribution Γ possesses a semi lower-continuous density γ w which is supported on an open set O w ; then X is called a non-linear state space model driven by F or NSS(F ) model, with control set O w .

Lemma 6 .

 6 [11, 7.2.4, 7.3.5] Let M ⊂ X be a minimal set for CM(F ). If CM(F ) is forward accessible then the NSS(F ) model restricted to M is an open set irreducible T-chain.

Lemma 7 .

 7 The Markov chain Z is open set and ψ-irreducible, aperiodic, and compacts of R are small-sets.

P

  (z, dy)g(y) < ∞ .

  and per recurrence [p k ] i N (0, I n ) for any k ∈ N. Therefore, we can apply the strong law of large numbers and 1/t t j=1 [p j ]

a

  i a j a k a l .

a-a 3 -→ t→∞ lim t→∞ a 6 ( 1 -a 4 -→ t→∞ a 6 ( 1

 361461 4i+3 1 -a 3(t-2-i) 1 -a)(1 -a 2 )(1 -a 3 ) 1 -a 4(t-2) 1 -a)(1 -a 2 )(1 -a 3 )(1 -a 4 )(38)By factorising Eq. (27), Eq. (31), Eq. (32), Eq. (37) and Eq. (38) by 1 1-a 4 we get the coefficients of Theorem 4.

Figure 2

 2 Figure1shows the time evolution of ln(σ t /σ 0 ) for 5001 runs and c = 1 (left) and c = 1/ √ n (right). By comparing Figure1aand Figure1bwe observe smaller variations of ln(σ t /σ 0 ) with the smaller value of c.Figure2shows the relative standard deviation of ln (σ t+1 /σ t ) (i.e. the standard deviation divided by its expected value). Lowering c, as shown in the left, decreases

Fig. 2 :

 2 Fig.2: Standard deviation of ln (σ t+1 /σ t ) relatively to its expectation. Here λ = 8. The curves were plotted using Eq. (21) and Eq. (22). On the left, curves for (right to left) n = 2, 20, 200 and 2000. On the right, different curves for (top to bottom) c = 1, 0.5, 0.2, 1/(1 + n 1/4 ), 1/(1 + n 1/3 ), 1/(1 + n 1/2 ) and 1/(1 + n).

Acknowledgments

This work was partially supported by the ANR-2010-COSI-002 grant (SIMINOLE) of the French National Research Agency and the ANR COSINUS project ANR-08-COSI-007-12.