Freek Stulp

Pierre-Yves Oudeyer
email: pierre-yves.oudeyer@inria.fr

Proximodistal Exploration in Motor Learning as an Emergent Property of Optimization

Keywords: human motor learning, proximo-distal exploration, stochastic optimization, modelling, evolution strategies, cross-entropy methods, policy search, morphology. Research highlights

To harness the complexity of their high-dimensional bodies during sensorimotor development, infants are guided by patterns of freezing and freeing of degrees of freedom. For instance, when learning to reach, infants free the degrees of freedom in their arm proximodistally, i.e. from joints that are closer to the body to those that are more distant. Here, we formulate and study computationally the hypothesis that such patterns can emerge spontaneously as the result of a family of stochastic optimization processes (evolution strategies with covariance-matrix adaptation), without an innate encoding of a maturational schedule. In particular, we present simulated experiments with an arm where a computational learner progressively acquires reaching skills through adaptive exploration, and we show that a proximodistal organization appears spontaneously, which we denote PDFF (ProximoDistal Freezing and Freeing of degrees of freedom). We also compare this emergent organization between different arm morphologies -from human-like to quite unnatural ones -to study the effect of different kinematic structures on the emergence of PDFF.

Introduction

As Bernstein emphasized [START_REF] Bernstein | The coordination and regulation of movements[END_REF], a great mystery in infant motor development is to understand how they can learn motor skills efficiently given a a complex non-linear body with many degrees of freedoms. Robots face the same problems, and this issue has similarly been the object of many studies in the recent years [START_REF] Vijayakumar | Incremental online learning in high dimensions[END_REF][START_REF] Baranes | The interaction of maturational constraints and intrinsic motivations in active motor development[END_REF][START_REF] Kober | Policy search for motor primitives in robotics[END_REF]Baranes & Oudeyer, 2013;[START_REF] Stulp | Robot skill learning: From reinforcement learning to evolution strategies[END_REF]. Learning motor skills involves experimenting with one's own body under limited time resources, and thus only a small fraction of physically possible movements can be sampled within the first years of life. Thus, as argued for example in [START_REF] Berthier | Proximodistal structure of early reaching in human infants[END_REF] and theoretically analyzed from a machine learning perspective [START_REF] Oudeyer | Intrinsically Motivated Learning of Real World Sensorimotor Skills with Developmental Constraints[END_REF], learning strategies based on simple forms of trial and error cannot lead to efficient learning in such contexts.

Several strands of research have studied families of mechanisms that could constrain and guide motor learning processes. In particular, Bernstein established a motor development perspective based on staged learning processes where some degrees of freedoms were first frozen, transforming a complex learning problem in a simpler one, and then progressively freed, allowing the learner to take advantage of the full potential of its body [START_REF] Bernstein | The coordination and regulation of movements[END_REF]. A number of experimental studies allowed to confirm this perspective. For example, Berthier et al. [START_REF] Berthier | Proximodistal structure of early reaching in human infants[END_REF] showed that the development of early reaching in infants [START_REF] Bertenthal | Eye, head and trunk control: The foundation for manual development[END_REF]) followed a proximodistal structure, where infants first learnt to reach by freezing the elbow and the hand, while varying shoulder and trunk movements, and then progressively used more distal joints of the elbow and hand. Studies in adult motor skill acquisition showed similar patterning of freezing and freeing of degrees of freedom, applied to the acquisition of racket skills [START_REF] Southard | Changing movement patterns: Effects of demonstration and practice[END_REF], soccer [START_REF] Hodges | Changes in coordination, control and outcome as a result of extended practice on a novel motor skill[END_REF] or skiing [START_REF] Vereijken | Free (z) ing degrees of freedom in skill acquisition[END_REF]. Other experimental observations have shown the complexity and contextdependance of this form of patterning, where for example infant reaching with different postural constraints could lead to higher use of elbow with respect to the shoulder [START_REF] Thelen | The transition to reaching: Mapping intention and intrinsic dynamics[END_REF].

Several hypotheses explaining the underlying mechanisms leading to such staged motor learning schedules were formulated so far. For example, Berthier et al. [START_REF] Berthier | Proximodistal structure of early reaching in human infants[END_REF] suggested that these learning schedules could be innate and due to the progressive neuromuscular development, where physiological maturation of motor neurons along the corticospinal tract [START_REF] Kuypers | Anatomy of the descending pathways[END_REF][START_REF] Jansen | The perinatal reorganization of the innervation of skeletal muscle in mammals[END_REF] could potentially lead to an initial limitation in the control of distal degrees of freedom. Yet, the extent to which physiological maturation can constrain motor exploration is still unclear in the infant [START_REF] Adolph | Physical and motor development[END_REF], and does not provide an explanation of the underlying mechanisms which drive freezing and freeing of degrees of freedom in adult motor learning.

In this article, we formulate, explore and analyze another (possibly complementary) hypothesis from a computational modelling perspective. This hypothesis is formulated within the optimal control framework of motor learning, where the learner uses exploration to find a motor program which minimizes a given cost (or maximizes an objective function) [START_REF] Todorov | Optimality principles in sensorimotor control[END_REF][START_REF] Berthier | Approximate optimal control as a model for motor learning[END_REF]. The hypothesis we study states that staged learning schedules with freezing and progressive freeing of degrees of freedom can self-organize spontaneously as a result of the interaction between certain families of stochastic optimization methods (which drive exploration of the learner) with physical properties of the body, and without involving physiological maturation. In particular, we present simulated experiments with a 6-DOF arm where a computational learner progressively acquires reaching skills (i.e. minimizing a cost to reach), and we show that a proximodistal organization appears spontaneously, which we denote PDFF (Proximo Distal Freezing and Freeing of degrees of freedom). We also compare the emergent structuration as different arm structures are used -from human-like to quite unnatural ones -to study the effect of different kinematic structures on the emergence of PDFF.

In these experiments1 , the reaching task is learned by applying stochastic optimization to optimize the parameters of a movement policy. The algorithm we use -PI BB , a special case of PI 2 -is based on covariance matrix adaptation through weighted averaging, which is a concept present in a wide range of optimization frameworks [START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF][START_REF] Rubinstein | The cross-entropy method: A unified approach to combinatorial optimization, monte-carlo simulation, and machine learning[END_REF][START_REF] Stulp | Robot skill learning: From reinforcement learning to evolution strategies[END_REF]. Covariance matrix adaptation allows the algorithm to determine dynamically the appropriate exploration magnitude and direction for each joint in order to progress fastest towards the goal at any given point in the development. In the context of PDFF, increasing and decreasing the exploration corresponds to freeing and freezing joints respectively.

However, to our knowledge, methods for exploration through covariance matrix adap-tation were so far analyzed only from an engineering perspective and in terms of speed to find optimal controllers. Here, on the contrary, we use these general methods as tools to modeling processes of exploration during motor learning in infants and study the patterns of freeing and freezing of DOFs that they generate. Preliminary work in this direction was presented in [START_REF] Stulp | Emergent proximo-distal maturation through adaptive exploration[END_REF], 2012a), but was based on more complex and specific optimization algorithms, did not include detailed analysis of results, and did not study how different morphologies of the body impacted the resulting patterns of exploration.

Here, we use a simple and generic form of covariance matrix adaptation -PI BB -and study how it spontaneously generates PDFF exploration patterns in the context of several arm morphologies. In the two analysis sections, we further study these results by considering the effect of joints on the cost in a static context. We first perform a sensitivity analysis by quantifying the effect of perturbing individual joint on the main cost component. We then analyze the interactions between joints by determining the effects of perturbing distal joints in the context of perturbations to proximal joints. The results of this analysis provide a deeper understanding of why PDFF arises during the stochastic optimization.

Limitations of Prior Research

Many computational models have studied how prewired stages or patterns of freezing and freeing of degrees of freedom could contribute or hinder learning of motor skills in highdimensions. Some studies considered the impact of alternation of freezing and freeing phases upon robot learning of swinging skills [START_REF] Berthouze | Motor skill acquisition under environmental perturbations: On the necessity of alternate freezing and freeing degrees of freedom[END_REF], studied how the pace of the sequencing of discrete stages [START_REF] Bongard | Morphological change in machines accelerates the evolution of robust behavior[END_REF][START_REF] Grupen | A developmental organization for robot behavior[END_REF][START_REF] Lee | Developmental learning for autonomous robots[END_REF] or of the continuous increase of explored values of DOFs along a proximodistal scheme [START_REF] Baranes | The interaction of maturational constraints and intrinsic motivations in active motor development[END_REF]) could be adaptively and non-linearly controlled by learning progress and lead to efficient motor learning in high-dimensional robots. Other related models have explored how the progressive freeing of degrees of freedom in the perceptual space [START_REF] Nagai | Learning for joint attention helped by functional development[END_REF][START_REF] French | The importance of starting blurry: Simulating improved basic-level category learning in infants due to weak visual acuity[END_REF], in the environment [START_REF] Uchibe | Environmental complexity control for vision-based learning mobile robot[END_REF], or in the structure of neural networks for learning abstractions [START_REF] Elman | Learning and development in neural networks: The importance of starting small[END_REF][START_REF] Westermann | Neuroconstructivism[END_REF] could guide the acquisition of sensorimotor and cognitive skills.

In all these models, the global scheduling of freezing and freeing degrees of freedom is encoded by the engineer (but the rhythm of progression from stage to stage can be adaptive as in [START_REF] Baranes | The interaction of maturational constraints and intrinsic motivations in active motor development[END_REF][START_REF] Lee | Developmental learning for autonomous robots[END_REF]). Some models have explored explicitly the evolutionary mechanisms that could generate and select such innate maturational schedules [START_REF] Cangelosi | Heterochrony and adaptation in developing neural networks[END_REF][START_REF] Matos | Heterochrony and evolvability in neural network development[END_REF].

A related model is presented in Schlesinger et al. [START_REF] Schlesinger | Learning to reach by constraining the movement search space[END_REF]. It is most similar to ours in that it also uses a kinematically simulated arm, and explores how evolutionary-like stochastic optimization methods can lead "several constraints to appear to fall out as a consequence of a relatively simple trial-and-error learning algorithm" [START_REF] Schlesinger | Learning to reach by constraining the movement search space[END_REF], one of them being the locking of joints. Movement policies are represented as four-layer feedforward neural networks, which are trained through evolutionary learning. In terms of the experimental setup, one main difference to our work is that we consider higher-dimensional systems -10 DOF instead of 3DOF -and use one learning agent instead of a population of 100. A second difference is that we employ a different family of stochastic optimization techniques, which has more flexibility in that it can dynamically update ranges of exploration during single agent learning. While the model in Schlesinger et al. [START_REF] Schlesinger | Learning to reach by constraining the movement search space[END_REF] only accounted for the freezing of some degrees of freedom as a result of optimization, the flexibility of our learning model allows us to find the entire developmental pattern outlined by Bernstein [START_REF] Bernstein | The coordination and regulation of movements[END_REF]: freezing of degrees of freedom followed by progressive and ordered freeing of degrees of freedom. We also consider various arm morphologies to show how the emergent scheduling adapts to the peculiarities of a given kinematic structure.

Learning to Reach with Stochastic Optimization

The methods, results and discussions of the experiments are distributed over three sections, corresponding to the three experiments conducted. In this first section, we describe an experiment in which a parameterized policy generates reaching movements, where the parameters of the policy are optimized through stochastic optimization.

Methods

Arm Model

The evaluation task in this paper consists of a kinematically simulated arm with M = 6 degrees of freedom, and a normalized length of 1. To study the effect of different kinematic structures on maturation, we use three sets of relative link lengths, depicted at the bottom of Figure 1: 1) typical relative link lengths of a human arm; 2) equidistant link lengths; 3) 'inverted' human arm, i.e. with short link lengths first.

Figure 1: Visualization of the task, which is to reach for a specific target location, in this case [0, 0.85]. The arm starts horizontally, and the movement to the target is visualized with a stroboscopic snapshots.

Task Specification

The main aim of the task is to reach for a specific target with a 0.5s movement, visualized in Figure 1. The angles and angular velocities of all joints are initially 0, which corresponds to a completely stretched 'horizontal' arm. The cost function in equations ( 1)-(2) consists of three parts, expressing different criteria to be optimized during the learning process:

Terminal cost ||x t N -x g || 2 .
The distance between the 2-D Cartesian coordinates of the end-effector (x t N ) at the end of the movement at t N , and the goal x g . This expresses that we want to reach to the target x g as closely as possible. Terminal cost max(q t N ). A cost that corresponds to the largest angle over all the joints at the end of the movement. This expresses an end-state comfort effect [START_REF] Cohen | Where grasps are made reveals how grasps are planned: generation and recall of motor plans[END_REF]. Immediate cost r t The immediate costs at each time step r t in (2) penalize joint accelerations. The weighting term (M + 1 -m) penalizes DOFs closer to the origin, the underlying motivation being that wrist movements are less costly than shoulder movements for humans, cf. [START_REF] Theodorou | A generalized path integral control approach to reinforcement learning[END_REF] 2 .

φt N = 10 2 ||xt N -x g || 2 + max(qt N ) Terminal cost (1) rt = 10 -5 M m=1 (M + 1 -m)(qt,m) 2 M m=1 (M + 1 -m) Immediate cost (2)
The factors 10 2 and 10 -5 have two purposes: 1) a scaling factor to compensate for different range of values the different cost components have. 2) a weighting factor enabling the prioritization of tasks. The order of priorities is: reach close to the target, achieve end-state comfort, minimize accelerations.

Figure 2: Explanation and visualization of the PI BB algorithm, using a 2D search space. For illustratory purposes, the cost of a sample θ is simply the distance to the origin: J(θ) = ||θ||. Left: PI BB pseudocode. Center: Visualization of one parameter update with PI BB . Right: Evolution of the parameters over several updates, demonstrating how the distribution converges towards the minimum θ * = [0, 0]. The algorithm is initialized by setting the mean and covariance parameters θ, Σ to θ init and λ init I respectively, visualized as a dark blue dot and circle. After this initialization, PI BB then iteratively updates these parameters with the following steps: 1) Explore. Sample K parameter vectors θ k from N (θ, Σ), and determine the cost J k of each sample. In the visualization of our illustratory example task K = 15, and the cost J(θ) is the distance to the origin ||θ||, which lies approximately between 8 and 19 in this example. 2) Evaluate. Determine the weight (probability) P k of each sample, given its cost. Essentially, low-cost samples have higher weights, and vice versa. The normalized exponentiation function that maps costs to weights is taken directly from the PI 2 algorithm, and is visualized in the center graph. Larger green circles correspond to higher weights. 3) Update. Updating the parameters θ, Σ with weighted averaging. In the visualization, the updated parameters are depicted in red. Because low-cost samples (e.g. a cost of 8-10) have higher weights, they contribute more to the update, and θ therefore moves in the direction of the optimum θ * = [0, 0].

Policy Representation

The policy representation encodes how movements are generated, by specifying the acceleration profiles of each joint. It is represented as a linear combination of normalized Gaussian basis functions. The acceleration qm,t of the m th joint at time t is determined as in (3), where the parameter vector θ m represents the weights of joint m. Intuitively, different basis functions are active at different times during the movement. The first basis function Ψ b=0 (t) is most active at the beginning of the movement, and the last Ψ b=B (t) at the end of the movement, with a cascade of basis functions in between. Setting different weights in the parameter vector θ thus leads to different acceleration profiles during the movement. If θ = 0, then there is no acceleration, and thus no movement. qm,t = g t θ m

Acceleration of joint m (3)

[g t ] b = Ψ b (t) B b=1 Ψ b (t) Time-dependent basis functions (4) Ψ b (t) = exp -(t -c b ) 2 /w 2 Kernel (5)
The centers c b=1...B of the kernels Ψ are spaced equidistantly in the 0.5s duration of the movement, and all have a width of w = 0.05s. The number of kernels per joint is B = 5. Since we do not simulate arm dynamics, the joint velocities and angles are acquired by integrating the accelerations. The end-effector "hand" position x is computed with the forward kinematics of the arm.

Policy Improvement through Stochastic Optimization

Stochastic optimization is based on iteratively exploring and updating parameters in a search space θ (θ is the vector of parameters of a movement policy). At each iteration, stochastic optimization algorithms generate K perturbations of the parameter vector {θ k = θ + k } K k=1 , compute the cost J k for each perturbation, and update the parameters θ → θ new based on these costs. This process continues until the costs have converged, or some termination condition is reached. In this article, the parameter space θ corresponds to the parameters of the policy. Optimizing policy parameters is known as direct policy search.

The specific stochastic optimization algorithm we use is PI BB , short for "Policy Improvement with Black-Box optimization" [START_REF] Stulp | Policy improvement methods: Between black-box optimization and episodic reinforcement learning[END_REF]. The PI BB algorithm is explained and visualized in Figure 2. We recommend readers to consider Figure 2 in detail, as it is important to understanding the rest of this paper. The main equations from Figure 2 are repeated in ( 6)-( 9).

{θ k = θ + k } K k=1
Explore by sampling K exploratory parameter vectors N (θ, Σ) (6)

J k = J(θ), P k = f (J k )
Evaluate by computing K costs and map costs to weights ( 7)

Σ new = K k=1 [P k (θ k -θ)(θ k -θ) ]
Update Σ for future exploration with weighted averaging.

(8)

θ new = K k=1 [P k θ k ]
Update the mean for future sampling with weighted averaging.

(9)

The core underlying principle in PI BB relevant to our experiments is using weighted averaging to update the mean θ (9) and covariance matrix Σ (8) of the sampling distribution. It shares this principle with many other evolution strategies algorithms such as "Cross-Entropy Methods" [START_REF] Rubinstein | The cross-entropy method: A unified approach to combinatorial optimization, monte-carlo simulation, and machine learning[END_REF], "Covariance Matrix Adaptation -Evolutionary Strategies" [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]) (CMA-ES) and "Policy Improvement with Path Integrals" [START_REF] Theodorou | A generalized path integral control approach to reinforcement learning[END_REF] (PI 2 ). Hoffman et al. have shown that weighted averaging better explains human motion learning than gradient-based algorithms [START_REF] Hoffmann | Optimization strategies in human reinforcement learning[END_REF].

The intuitive meaning is that the perturbations that are sampled at each iteration of learning are themselves adaptive, fostering exploration in directions where the cost decreases fastest. Indeed, PI BB is a special case of both PI 2 (without temporal averaging and with covariance matrix adaptation) and CMA-ES (without evolution paths), which are state-ofthe-art in direct policy search and black-box optimization respectively.

Here we select PI BB over PI 2 because it is the simplest such algorithm that implements the principle of weighted averaging to update policy parameters; we are interested in studying the formation of staged patterns of freezing and freeing of degrees of freedom (PDFF), not in achieving for instance the fastest possible convergence. For a full explanation of PI BB , and its relationship to PI 2 and CMA-ES, we refer to [START_REF] Stulp | Policy improvement methods: Between black-box optimization and episodic reinforcement learning[END_REF].

Adaptive Exploration through Covariance Matrix Adaptation

Covariance matrix adaptation allows PI BB to automatically adapt the exploration so as to generate more samples in the direction of the minimum. Because this property is important for PDFF, we highlight and illustrate it with the two simple examples in Figure 3.

In the left graph in Figure 3 (repeated from Figure 2), the current parameters θ = [10, 10] are far from the optimum θ * = [0, 0]. The samples that are closer to θ * (to the lower left in the graph) have larger weights (visualized as green circles) than those that are further away from θ * . Because the new parameters are weighted with these weights when averaging, the mean θ moves in the direction of these high-weight, low-cost samples when updating, bringing θ closer to θ * . The same principle applies to the covariance matrix, which becomes elongated -its largest eigenvalue λ increases, and corresponds to the eigenvector pointing more in the direction of θ * (visualized as an arrow). More formally, it follows the natural gradient of the cost with respect to the parameters [START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]. The effect is that exploration increases, and in the direction of low-cost regions in the parameter space.

Figure 3: Left: when the current parameters are far from the optimum (θ * = [0, 0] in this example), the covariance matrix is updated (red ellipse) so that it tends to elongate in the direction of steepest descent, leading to increased exploration along this direction (see arrow). Right: when the current parameters are exactly at the optimum, the covariance matrix tends to shrink, leading to decreased exploration.

In the right graph, the initial covariance matrix is the same, but the current parameters θ = [0, 0] are now perfectly at the optimum θ * = [0, 0]. For this reason, samples closer to θ will have a lower cost J(θ) = ||θ||, and thus a higher weight. Note the larger green circles are all close to the center. Therefore, θ hardly moves after updating; this is desirable, because θ is already at the optimum. However, the covariance matrix shrinks (see smaller eigenvalue as arrow) because closer samples have higher weights. Therefore, exploration decreases in all directions.

This adaptive exploration behavior may also be observed in the right graph in Figure 2; in the first few updates the covariance matrix elongates towards the optimum θ * = [0, 0], but once it is reached, the covariance matrix shrinks and exploration decreases. This behavior is especially apparent in this idealized example, which uses only a 2-D search space and a very benign cost function; however, the general principle also applies to high-dimensional spaces and discontinuous cost functions, as demonstrated in [START_REF] Stulp | Adaptive exploration for continual reinforcement learning[END_REF].

Application of PI BB to the Reaching Task

In this article, we apply PI BB to the parameters of the policy representation previously described. Each joint has its own parameters θ m and covariance matrix Σ m , which are iteratively updated with PI BB . The input parameters of PI BB are set as follows. The initial parameter vector is θ init = 0, which means the arm is completely stretched, and not moving at all over time. The initial and minimum exploration magnitude of each joint m is set to λ init = λ min = 0.05. The number of trials per update is K = 20, and the eliteness parameter is h = 10, the default value suggested by [START_REF] Theodorou | A generalized path integral control approach to reinforcement learning[END_REF] 3 A separate stochastic optimization session was run 10 times for each of the 20 target points, i.e. 200 sessions per arm structure.

The exploration magnitude λ m of a particular joint m at some point during the learning process is defined as the maximum eigenvalue of the covariance matrix Σ m . Initially, λ m is λ init , because the all eigenvalues of the initial diagonal Σ = λ init I are λ init . The length of the two arrows in Figure 3 visualize λ m for non-diagonal covariance matrices, which may arise as Σ is updated.

In the context of this paper, we consider joints that have low exploration magnitudes λ m to be 'frozen', whereas those with high λ m are 'free'. The thick vertical bars indicates when a joint reached its maximum relative exploration magnitude (position on x-axis), as well as the magnitude itself (the height of the bar). For example, for the human arm structure (top graph), the first joint achieves a maximum relative exploration magnitude of 0.56 at update 7. This means that 56% of the exploration is accounted for by only the most proximal (first) joint.

Results

To investigate the robustness of the method against initial conditions, Figure 5 plots the variability around the mean of the exploration in the first joint in the human arm for the 200 learning sessions for this arm type. In generating this figure, we noticed that the exact onset of increasing exploration (freeing degrees of freedom) is influenced by the stochastic nature of the optimization algorithm. To factor this out, we have applied dynamic time warping [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] to the exploration curves of the 200 individual learning sessions before computing the mean and variance.

Discussion

For all arm structures, we see that the total exploration (thick black/yellow line) initially increases, indicating that DOFs are globally freed. After achieving a maximum total exploration at around update 20-25, exploration then decreases again once the task has been learned. This behavior is a direct consequence of the adaptive exploration described in Figure 3, and is also observed in [START_REF] Stulp | Emergent proximo-distal maturation through adaptive exploration[END_REF][START_REF] Stulp | Adaptive exploration for continual reinforcement learning[END_REF].

The relative exploration magnitude between the joints, however, shows quite a different development for the different arm structures. For the human arm (top graph), the most proximal joint is already responsible for more than 50% of the exploration (0.56) after 7 updates. This joint has been freed, whereas the others are frozen; the three distal joints account for less than 20% of exploration at update 7. For the human arm, we see that joints 1, 2, 3 achieve their maximum relative exploration of 0.56, 0.42 and 0.26 at updates 7, 18 and 27 respectively. In conclusion, we clearly see that the first three joints are freed in a proximodistal order4 .

For the equidistant arm, the order in which the joints achieve their maxima is 2,3,4,1,5,6. Thus, apart from the most proximal joint, we again see a proximodistal freeing of joints. When considering the maximum relative magnitudes of the exploration (vertical bars), we see however that the freeing/freezing of joints is much less less pronounced. None of the maxima exceeds 0.3, which is in contrast with the human arm, where the most proximal joint is responsible for 0.56 if the exploration. For the equidistant arm, the exploration is thus spread out over the joints much more, rather than being focused in only one or several joints.

Finally, for the inverted human arm, the order in which joints are freed is 3,4,5,6,1,2. This time, the bulk of the exploration is being done by joint 6, which accounts for more than 42% of exploration at update 50. The other joints again never exceed 0.3. Thus, for this rather unnatural arm, we see more exploration in the distal joints; only later on do the proximal joints 1 and 2 achieve their maximum values.

In summary, the results show that there is a consistent emergent organization of exploration over time in all arm structures. The PDFF organization is quite pronounced in the human arm, where the order of freeing is 1,2,3,5,4, and the exploration switches most clearly (i.e. high relative magnitudes of exploration) from one joint to the other. It is important to realize that this effect emerges solely from adaptive exploration through covariance matrix adaptation, and the order and/or stages in which degrees-of-freedom are freed is not pre-defined, as in for instance [START_REF] Berthouze | Motor skill acquisition under environmental perturbations: On the necessity of alternate freezing and freeing degrees of freedom[END_REF][START_REF] Bongard | Morphological change in machines accelerates the evolution of robust behavior[END_REF][START_REF] Baranes | The interaction of maturational constraints and intrinsic motivations in active motor development[END_REF].

Analysis: Individual Joints

In this section, we analyze how and why PDFF arises when applying stochastic optimization, and how and why this depends on the arm structure. We perform the analysis in a static context -static because we do not perturb the parameters of the policy that determines the joint angles over time, but rather perturb the joint angles directly without a temporal component. This analysis helps to understand why PDFF arises within an optimization context. This figure illustrates that the variance is higher when exploration is highest, as is to be expected. The variances for the other joints, not plotted here, are lower.

Here, we consider the effect that perturbations of individual joints have on the cost through sensitivity analysis. Sensitivity analysis aims at "providing an understanding of how the model response variables respond to changes in the input." [START_REF] Saltelli | Sensitivity analysis[END_REF]. We use sensitivity analysis to investigate how the variation in individual joint angles -the input -influences the variation in the cost -the response variables. This provides a first indication of why PDFF arises.

Methods

In the default posture, all joint angles are zero. This posture is perturbed by setting one of the 6 joint angles to π 10 . The 6 possible perturbations, one for each joint, are visualized in the top row of Figure 6. For the default and perturbed configuration, we then compute the distance of the end-effector to the target ||x t N -x g ||. Because this is a static context there are no joint accelerations, and the immediate costs (2) are not included. The lower row plots the difference in cost between the outstretched arm (where all joints are 0), and the slightly bent arm (where one joint angle is π 10 ). To acquire a value that is representative for the whole workspace, the differences in the lower row of Figure 6 is the average over the 20 target positions x g depicted in the top row of Figure 6.

Results

For all arm configurations, we see that proximal joints lead to a higher average difference in the distance to the target than more distal ones. This should not come as a surprise, as rotating more proximal joint leads to smaller movement in the end-effector space, and it is the end-effector space that determines the distance to the target. As a consequence, the same magnitude of perturbation will lead to a larger difference in cost for more proximal joints.

Figure 6: Results of the sensitivity analysis. For all arm configurations, we see that proximal joints lead to a higher average difference in the distance to the target than more distal ones.

Discussion

The goal of stochastic optimization is to minimize costs through exploring and updating in parameter space. The results in Figure 6 demonstrate that perturbing proximal joints leads to larger differences in costs than distal joints. Therefore, an optimizer can be expected to minimize costs more quickly if it initially focuses exploration on proximal joints, rather than distal ones. This may be an explanation why exploration is larger in more proximal joints, but does not explain why distal joints are not also freed. This is the aim of our second analysis, which now follows.

Analysis: Interactions Between Joints

Whereas the previous section on sensitivity analysis considered joints individually, we now turn to the interaction between pairs of joints. We especially focus on how perturbations in proximal joints affect the influence of perturbations in more distal joints on the cost.

Methods

Pairs of joints are considered. For the more proximal joint, two perturbations {P 1, P 2} are sampled from N (0, π 10 ). For each perturbation of the proximal joint, the distal joint is perturbed twice, also by sampling from N (0, π 10 ). This leads to the four arm configurations in Figure 7. The question we ask for these four configurations is: does the perturbation of the distal joint change the cost ranking? This is the same question underlying uncertainty handling in ranked-based evolutionary direct policy search [START_REF] Heidrich-Meisner | Uncertainty handling in evolutionary direct policy search[END_REF]. Figure 7 depicts examples where the answer is no (left) and yes (right). The question is asked for 100 different samples, and for each of the 20 target points. The average of these 2000 values represents the ratio that the answer was 'no', i.e. a ratio of 1 implies that, no matter how much the distal joint is perturbed, it does not affect the ranking. A value of 0.5 implies that the perturbation of the distal joint affects the ranking half the time.

Results

Figure 8 depicts the ratio for pairs of joints for all three arm configurations. For proximal joint 1 and distal joint 3 -the case depicted in Figure 7 -this value is 0.89, labeled (A).

Figure 7: Examples of perturbing pairs of joints for the human arm, and the effect of these perturbations on the cost ranking (numbers at the end of the arm).

Thus, when joint 1 is perturbed, the perturbation of joint 3 affect the cost ranking in only 11% of the samples and target points.

Discussion

The main implication of these results for stochastic optimization is that if we are performing exploration with for instance joint 1, exploration in more distal joints is less relevant to the cost. For the human arm for instance, it is not sensible to explore with joint 6 when exploring with joint 1, because joint 6 only affects the result 2% of the time. Thus, from the point of view of stochastic optimization, joint 6 may well be frozen when searching in joint 1.

It is interesting to see that for the human arm configuration, proximal joints dominate distal joint much more (median=0.89) than an equidistant arm (median=0.79), and even more so for the 'inverted' human arm (median=0.70). Thus, exploring with more distal joints has a lower impact on the cost ranking for the human arm, and thus the effect of PDFF may be expected to be stronger for this configuration. This is confirmed by the empirical results in Figure 4, where PDFF is most pronounced in the human arm configuration.

General Discussion and Conclusion

Staged motor learning, and in particular the progressive freezing and freeing of degrees of freedom observed in infant and adult motor exploration, has been argued to facilitate human acquisition of high-dimensional motor skills, and was also shown to be efficient for robot motor learning. Several hypotheses explaining the underlying mechanisms leading to such staged motor learning schedules were formulated so far, but have mostly relied on forms of innate scheduling of patterns of freeing and freezing.

Here, in the framework of approximate optimal control, we have studied the hypothesis that staged learning schedules with freezing and progressive freeing of degrees of freedom can self-organize spontaneously as a result of the interaction between certain families of stochastic optimization methods with the physical properties of the body, and without involving physiological maturation. In particular, we have presented simulated experiments with a 6-DOF arm where a computational learner progressively acquired reaching skills, and we showed that a proximodistal organization appeared spontaneously. We also compared the emergent structuration as different arm structures are used -from human-like to quite unnatural ones -to study the effect of different kinematic structures on the emergence of PDFF. We analyzed further these results through a sensitivity analysis, providing a deeper understanding of why PDFF arises during the stochastic optimization.

Parsimony and biological plausibility Overall, this model does not invalidate the hypothesis that an an innate maturational scheduling for freezing and freezing DOFs can be involved in infant motor learning. However, it shows that relatively simple stochastic optimization processes with adaptive exploration, which Bernstein already suggested were at play in infants, can already account for the formation of patterns of staged motor exploration. Yet, while the form of adaptive stochastic optimization we have considered is simple and general, one can wonder whether such mechanisms could be actually implemented in biologically plausible neural networks. These mechanisms are a form of evolutionary optimization algorithms that are based on two complementary principles: the capacity to make variations/mutations of current good solutions (and to select the most useful ones), and the capacity to identify which directions of variation/mutation are currently most improving the current good solutions, involving a form of memory of past explorations. These two principles, and other more complex forms of Darwinian search processes in the brain, have been shown to be neurally plausible by several lines of research (C. T. [START_REF] Fernando | Selectionist and evolutionary approaches to brain function: a critical appraisal[END_REF], building on Edelman theory of neuronal group selection [START_REF] Edelman | Neural darwinism: The theory of neuronal group selection[END_REF]), Changeux's theory of synaptic selection and selective stabilization [START_REF] Changeux | A theory of the epigenesis of neuronal networks by selective stabilization of synapses[END_REF]), Calvin's replicating activity patterns [START_REF] Calvin | The brain as a darwin machine[END_REF]. In particular, the recently developped Neuronal Replicator Hypothesis (C. [START_REF] Fernando | The neuronal replicator hypothesis[END_REF] has shown how various known neuronal physiological mechanisms could implement such general genetic algorithms, including the mechanisms of adaptive exploration that we have been using in the model presented in this paper. An implementation of these mechanisms was shown to work with a neural network using realistic Izhikevich spiking neurons (C. [START_REF] Fernando | Evolvable neuronal paths: a novel basis for information and search in the brain[END_REF]. To summarize, the mechanisms used to allow incremental exploration and learning in the model presented here are not only simple and parsimonious, but they are also implementable in realistic neural networks.

Open questions and research perspectives. Our general hypothesis therefore forms a baseline against which more complex, domain-specific hypotheses should be compared. Also, the spontaneous formation of PDFF patterns through stochastic optimization appears to be compatible with observation of the patterns of motor exploration in adult motor learning, where maturational mechanisms have little probability to be at play. This new hypothesis also points to several open questions and new experimental investigations. In adults, several experiments have shown that a structuration of exploration through freezing and freezing of degrees of freedom happened [START_REF] Southard | Changing movement patterns: Effects of demonstration and practice[END_REF][START_REF] Hodges | Changes in coordination, control and outcome as a result of extended practice on a novel motor skill[END_REF][START_REF] Vereijken | Free (z) ing degrees of freedom in skill acquisition[END_REF], however to our knowledge these experiments did not systematically study to what extent the freeing or freezing of particular degrees of freedom was correlated to their current usefulness in progressing towards a goal. Also, they considered tasks such as skiing, soccer or racket skills that were culturally known by subjects and thus could involve other mechanisms such as imitation learning, complicating the analysis of the exploration strategies. In infants, it is also an open question to know whether the adaptive exploration mechanisms shown by adults are already at play, or whether they are not yet in place and exploration is rather controlled by maturational mechanisms such as myelination.

To address these open questions, one could imagine using experimental setups in which human subjects have to learn and explore new sensorimotor mapping that are highly different from what they already know, and where the relation between degrees of freedom and their usefulness for progressing towards the goal can be controlled systematically. [START_REF] Miard | A new experimental setup to study the structure of curiositydriven exploration in humans[END_REF] have proposed an experimental setup that could be used with adults in this context (but was not used for these specific questions so far): subjects have to learn how to control an abstract visual shape on a screen by using movements of their body as measured by a 3D camera (and the mathematical form of the relation between body movements and the abstract visual shape can be changed systematically). Similar setups could be imagined for infants, taking inspiration from the famous Rovee-Collier task (Rovee-Collier, 1999) (developped for studying other questions): infants' movements (arms, legs) could be tracked with sensors and used to control the intensity/frequency of a sound or color of a light, through a mapping where the relation between degrees of freedom and their impact on the sound/light could be systematically changed to cancel out possible effects of myelination in the structuration of exploration. For example, it could be possible to program an inverse proximodistal relationship between arm movements and the sound/light: movements of the tip of the arm could be made to have more impact on the sound/light than movements of the shoulder. In such a case, observing an infant exploration with a corresponding inverse proximodistal law would reinforce our hypothesis that adaptive exploration plays an important role, while observing a standard proximodistal exploration would invalidate the hypothesis that such a mechanism is present or can have a leading role in exploration.

Complementarity with other mechanisms It is also an open question to undertand how such adaptive exploration based on stochastic optimization could interact with other mechanisms guiding exploration such as maturation through myelination, imitation/social guidance, and intrinsic motivation. Several models in robotics have began to explore these links. [START_REF] Baranes | The interaction of maturational constraints and intrinsic motivations in active motor development[END_REF] have studied the efficiency of combining stochastic optimization to reach goals with maturational mechanism which progressively grow the limits within which stochastic optimization can physically explore, showing an increase in efficiency from a machine learning point of view. Several works have shown how human demonstration of movements could bootstrap this optimization process (e.g. [START_REF] Stulp | Simultaneous on-line discovery and improvement of robotic skill options[END_REF], or how humans can progressively shape subparts of the movements to complement autonomous exploration [START_REF] Chernova | Robot learning from human teachers[END_REF]. Finally, exploration in infants is also highly driven by mechanisms of intrinsic motivation (also called curiosity), where instead of trying to reach a goal imposed by social peers or the experimenter (as in the model presented in this paper), they use intrinsic criteria such as information gain or surprise to set their own goals and choose how to practice these self-selected goals [START_REF] Gottlieb | Information-seeking, curiosity, and attention: computational and neural mechanisms[END_REF][START_REF] Moulin-Frier | Self-organization of early vocal development in infants and machines: the role of intrinsic motivation[END_REF]. An important side effect of exploring multiple self-selected goals is that transfer learning across goals happens and in turn can shape the selection of future goals and ways to use degrees of freedom to explore them. Several computational architectures of intrinsically motivated learning have used stochastic optimization as the lower-level mechanism to learn how to reach self-selected goals in sensorimotor learning and achieve transfer learning across goals (Baranes & Oudeyer, 2013). These integrated architecture have also shown the self-organization of developmental structure, such as the transition from nonarticulated speech sounds to articulated vowels to proto-syllables in models if infant vocal development [START_REF] Moulin-Frier | Self-organization of early vocal development in infants and machines: the role of intrinsic motivation[END_REF][START_REF] Oudeyer | How evolution may work through curiositydriven developmental process[END_REF].

Figure 4

 4 Figure 4 presents the results of applying PI BB to the three different arm structures. Each graph plots the total exploration magnitude over all joints M m=1 λ m (thick yellow/black line) and the relative exploration magnitude λ m / M m=1 λ m per joint (colored patches) as learning progresses over 100 updates, which corresponds to 2000 = 100 • K roll-outs. All values are averaged over the 20 target points and 10 optimization sessions per target point.The thick vertical bars indicates when a joint reached its maximum relative exploration magnitude (position on x-axis), as well as the magnitude itself (the height of the bar). For example, for the human arm structure (top graph), the first joint achieves a maximum relative exploration magnitude of 0.56 at update 7. This means that 56% of the exploration is accounted for by only the most proximal (first) joint.To investigate the robustness of the method against initial conditions, Figure5plots the variability around the mean of the exploration in the first joint in the human arm for the 200 learning sessions for this arm type. In generating this figure, we noticed that the exact onset of increasing exploration (freeing degrees of freedom) is influenced by the stochastic nature of the optimization algorithm. To factor this out, we have applied dynamic time warping[START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] to the exploration curves of the 200 individual learning sessions before computing the mean and variance.

Figure 4 :

 4 Figure 4: Results of stochastic optimization for the three arm structures.The thick black/yellow line represents the total exploration magnitude (right y-axis) and relative exploration magnitudes of the 6 individual joints as learning progresses (left y-axis).

Figure 5 :

 5 Figure5: Mean and standard deviation in the exploration magnitude of the first joint in the human arm at each update, after aligning the 200 sessions for this arm with dynamic time warping. This figure illustrates that the variance is higher when exploration is highest, as is to be expected. The variances for the other joints, not plotted here, are lower.

Figure 8 :

 8 Figure 8: For pairs of joints, the height of the bars represents the ratio that the perturbation of the distal joint does not affect the cost ranking of the resulting postures.

The Matlab code used to generate and visualize the results in this article is available as open source, and can be downloaded here: https://github.com/stulp/dmp bbo/archive/proximodistalmaturation.zip

This cost term was taken from[START_REF] Theodorou | A generalized path integral control approach to reinforcement learning[END_REF]. In the context of this paper, it cannot be the reason for the PDFF we shall see in later sections. Rather than favoring PDFF, this cost term actually works against it, as proximal joints are penalized more for the accelerations that arise due to exploration.

High values of the eliteness h lead only a few samples to contribute to the weighted averaging. h = 0 would give all samples equal weight independent of the cost, and no learning would occur. As[START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] note: "In general, the selection related parameters [such as h] are comparatively uncritical and can be chosen in a wide range without disturbing the adaptation procedure." In fact, the same parameter settings have been used in entirely different domains, for instance to optimize robot control policies[START_REF] Theodorou | A generalized path integral control approach to reinforcement learning[END_REF].

Note that the 6th joint achieves its maximum at update 1. This is not because it is freed very early, but rather because it is almost frozen throughout the learning process, and thus achieves its maximum when the exploration is initially the same for all joints