Yoshiyuki Asai
email: asai@bpe.es.osaka-u.ac.jp

Alessandro E P Villa
email: alessandro.villa@ujf-grenoble.fr

Transmission of distributed deterministic temporal information through a diverging/converging three-layers neural network

Keywords: spiking neural networks, synfire chains, adaptive threshold neurons, computational neuroscience, preferred firing sequences

This study investigates the ability of a diverging/converging neural network to transmit and integrate a complex temporally organized activity embedded in afferent spike trains. The temporal information is originally generated by a deterministic nonlinear dynamical system whose parameters determine a chaotic attractor. We present the simulations obtained with a network formed by simple spiking neurons (SSN) and a network formed by a multiple-timescale adaptive threshold neurons (MAT). The assessment of the temporal structure embedded in the spike trains is carried out by sorting the preferred firing sequences detected by the pattern grouping algorithm (PGA). The results suggest that adaptive threshold neurons are much more efficient in maintaining a specific temporal structure distributed across multiple spike trains throughout the layers of a feed-forward network.

Introduction

A neuronal network can be considered as a highly complex nonlinear dynamical system able to exhibit deterministic chaotic behavior, as suggested by the experimental observations of single unit spike trains, which are sequences of the exact timing of the occurrences of action potentials [START_REF] Mpitsos | Chaos in brain function and the problem of nonstationarity: a commentary[END_REF][START_REF] Celletti | Determination of chaotic attractors in the rat brain[END_REF]. Previous studies [START_REF] Tetko | A comparative study of pattern detection algorithm and dynamical system approach using simulated spike trains[END_REF][START_REF] Asai | Detection of a dynamical system attractor from spike train analysis[END_REF] showed that deterministic nonlinear dynamics in noisy time series could be detected by applying algorithms aimed at finding preferred firing sequences with millisecond order time precision from simultaneously recorded neural activities. A neural network is also characterized by the presence of background activity of unspecified or unknown origin that is often represented by stochastic inputs to each cell of the network. Then, a neuron belonging to a cell assembly, somehow associated to a deterministic nonlinear system, within the network is expected to receive inputs characterized by an embedded temporal structure as well as inputs corresponding to the stochastic background activity. It has been shown that the characteristic transfer function of a neuronal model and the statistical feature of the the background activity may affect the transmission of temporal information through synaptic links [START_REF] Asai | Deterministic neural dynamics transmitted through neural networks[END_REF].

In the current paper we extend our previous analysis to diverging/converging feed-forward neuronal networks-synfire chains-which are supposed to represent the most appropriate circuits able to transmit information with the best temporal accuracy [START_REF] Abeles | Local Cortical Circuits[END_REF]. Moreover the temporally organized activity was fed to the network in a distributed way across the input spike trains [START_REF] Asai | Spatio temporal filtering of the distributed spike train with deterministic structure by ensemble of spiking neurons[END_REF]. We suggest that adaptive threshold neurons are much more efficient in maintaining a specific temporal structure throughout the layers of a synfire chain.

Methods

Spiking neuron model

We investigated two neuron models aimed to reproduce the dynamics of regular spiking neurons. The first is a simple spiking neuron (SSN) [START_REF] Izhikevich | Which model to use for cortical spiking neurons?[END_REF] described as:

dv dt = 0.04v 2 + 5v + 140 -u + I ext (t) , ( 1
)
du dt = a(bv -u) ,
with the auxiliary after-spike resetting, v ← c and u ← u + d when v ≥ +30 mV . v represents the membrane potential [mV ], u is a membrane recovery variable, a and b control the time scale of the recovery variable and its sensitivity to the subthreshold fluctuation of the membrane potential. This model generates an action potential with a continuous dynamics followed by a hyperpolarization modeled as a discontinuous resetting. Parameters were set as a = 0.02, b = 0.2, c = -65, d = 8 so to mimic the behavior of a regular spiking neuron [START_REF] Izhikevich | Which model to use for cortical spiking neurons?[END_REF].

The second model is a multiple-timescale adaptive threshold (MAT) model [START_REF] Kobayashi | Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold[END_REF] derived from [START_REF] Brette | Adaptive exponential integrate-and-fire model as an effective description of neuronal activity[END_REF]. In this model, the dynamics of the membrane potential is described as a non-resetting leaky integrator,

τ m dV dt = -V + R A I ext (t) , (2)
where τ m , V, R and A are the membrane time constant, membrane potential, membrane resistance, and scaling factor, respectively. A spike is generated when the membrane potential V reaches the adaptive spike threshold θ(t),

θ(t) = ω + H 1 (t) + H 2 (t) , dH 1 dt = -H 1 /τ 1 , ( 3
)
dH 2 dt = -H 2 /τ 2 ,
where ω is the resting value, H 1 and H 2 are components of the fast and slow threshold dynamics (characterized by decaying time constants τ 1 and τ 2 , respectively) which has a discrete jump when V (t) ≥ θ(t),

H 1 = H 1 + α 1 , H 2 = H 2 + α 2 . ( 4
)
Parameters were set to values τ m = 5 ms, R = 50 MΩ, A = 0.106, ω = 19 mV, τ 1 = 10 ms, τ 2 = 200 ms, α 1 = 37 mV, α 2 = 2 mV. The model with the above parameter values reproduces the activity of a regular spiking neuron [START_REF] Kobayashi | Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold[END_REF]. Let us denote I ext the input synaptic current, defined as

I ext = -A ext ∑ k g syn (t -t k ) , ( 5
)
where A ext is an intensity of the synaptic transmission of the spike received as an external input (A ext = 1 was used here for all simulations), t k represents time when the k-th spike arrives to the neuron model, and g syn is the post synaptic conductance represented by

g syn (t) = C 0 e -t/τ1 -e -t/τ2 τ1 -τ2 , ( 6
)
where τ1 and τ2 are rise and decay time constants given by 0.17 and 4 ms, respectively, and C 0 is a coefficient used to normalize the maximum amplitude of g syn (t) to 1. Notice that a single synaptic current given to a neuron is not strong enough to evoke post-synaptic neuronal discharges. Hence, it is necessary for a post-synaptic neuron to integrate several arriving synaptic currents for a spike generation.

Input spike train

We consider the deterministic dynamical system described by Zaslavskii [START_REF] Zaslavskii | The simplest case of a strange attractor[END_REF]:

{ x n+1 = x n + v(1 + µy n ) + εvµ cos x n , (mod. 2π) y n+1 = e -γ (y n + ε cos x n ) , (7)
where x, y, µ, v ∈ R, the parameters are µ = 1-e -γ γ , v = 400 3 and initial conditions set to x 0 = y 0 = 0.3. With this parameter set the system exhibits a chaotic behavior. Time series {x n } are generated by iterative calculation. A new time series {w n } corresponding to the sequence of the inter-spike-intervals is derived by

w n = x n+1 -x n + C, where C = min{(x n+1 -x n )} + 0.1 is a constant to make sure w n > 0.
The dynamics was rescaled in milliseconds time units with an average rate of 5 events/s (i.e., 5 spikes/s) in order to let the mean rate of the Zaslavskii spike train be comparable to neurophysiological experimental data. We calculated N = 10000 points of time series {w n } which corresponds to a spike train lasting L = 2000 seconds.

Given a dynamical information ratio D, where 0 ≤ D ≤ 1, a percentage of spikes corresponding to (1 -D) × 100 % are selected at random (uniformly

Neuronal network

We consider a diverging/converging neural network composed of three layers (Fig. 2). Each layer includes 20 neurons characterized by the same neuronal model with identical parameter values. Each neuron belonging to the first layer receives fifteen input spike trains randomly selected out of the twenty that were generated for a given dynamical information ratio D. Hence, a neuron in the first layer receives afferences from 15 input spike trains (each one firing on average at = 5 spikes/s) and an independent Poissonian spike train with a mean firing rate of 425 spikes/s as background activity. This means a neuron of the first layer integrates about 500 spikes in 1000 millisecond by the fourth order Runge-Kutta numerical integration method with 0.01 ms time steps. Each neuron of the next layer receives afferences from 15 neurons randomly selected in the previous layer.

In addition, each neuron receives an independent Poissonian spike train with a mean firing rate of 425 spikes/s as background activity. We observed that those neurons integrated between 490 and 540 spikes in 1000 ms The explicit synaptic transmission delay is not considered here. All connections were hardwired, and no synaptic plasticity was taken into account. Each simulation run lasted 2000 s.

Pattern detection and reconstruction of time series

Subsets of spike trains were obtained by using the Pattern Grouping Algorithm (PGA) [START_REF] Villa | Spatiotemporal activity patterns detected from single cell measurements from behaving animals[END_REF][START_REF] Tetko | A pattern grouping algorithm for analysis of spatiotemporal patterns in neuronal spike trains. 1. detection of repeated patterns[END_REF][START_REF] Abeles | Detecting precise firing sequences in experimental data[END_REF] as follows. Firing sequences repeating at least 5 times and above the chance level (p = 0.05) are detected by PGA. The interval between the first and the last spike of the firing sequence defines the duration of the pattern that was set to ≤ 600 ms. Given a maximum allowed jitter in spike timing accuracy (±3 ms) clusters of firing sequences are represented by a template pattern. For example, if there are 9 triplets (i.e., firing sequences formed by 3 spikes) belonging to the same cluster, a subset of the original spike train that includes 27 spikes (= 9 × 3) can be determined by a template pattern. Then, the subset of the original spike train referred to as "reconstructed spike train" is obtained by pooling all spikes belonging to all template pattern clusters [START_REF] Asai | Detection of a dynamical system attractor from spike train analysis[END_REF]. The reconstructed spike train from the original Zaslavskii series included 92.3% of the original spikes and its return map is shown in Fig. 3a. In a case of a Poissonian spike train with an average rate of 5 spikes/s the reconstructed spike train included only 0.4% spikes of the original series (Fig. 3b). Moreover, we have measured the dispersion of spike distribution by the Fano factor [START_REF] Sacerdote | On the classification of experimental data modeled via a stochastic leaky integrate and fire model through boundary values[END_REF], which is F = 1 for a Poissonian spike train, and the similarity ratio S between two spike trains defined as follows. Suppose that spike trains A and B contain N A and N B spikes and M spikes occur in A and B at the same time. The similarity ratio is defined by S = 2M/(N A + N B ), which is S = 1 for two identical spike trains. If we allow the coincidence to occur within a given jitter (∆ = 5 ms here), then the condition t n B -∆ ≤ t k A ≤ t n B + ∆ satisfies the coincidence of the n-th spike in train B with the k-th spike in train A.

Results

We investigated the continuous dynamics of the membrane potential for neurons characterized by the models SSN and MAT and analyzed their output spike trains at all layers. Table 1 summarizes the mean firing rates as a function of the layer and of dynamical information ratio D. The rates increased with an increase of D and for the same D they increased with the order of the layer.

In the 1st layer we analyzed the effect of the model by comparing cells that received the same inputs. Figure 4 shows the example of two different neurons (cells no. 114 and 115) located in the 1st layer. In the bottom panel the input spike trains with dynamical information ratio D = 0.5 and the Poissonian background are sorted in order to emphasize the spikes belonging to Zaslavskii. Zaslavskii spikes increase the chance to overlap and to produce a stronger postsynaptic current by temporal summation with an increase in D. In this example, eight spikes belonging to the original Zaslavskii spike train arrive simultaneously at t = 2150 ms (see the upward arrow in the last panel of Fig. 4) and evoke a suprathreshold current that generates a spike.

The return maps of the raw output spike trains of one representative neuron of each layer and for each neuronal model are shown in Fig. 5a as a function of D. As D decreased, the attractor contour become blurred. Notice that for exclusive Poissonian input spike trains (D = 0) the return maps of the SSN model (Fig. 5a (rightmost column) show a bias in the distribution of points, with empty bands Table 1. Mean firing rate (spikes/s) of a neuron of SSN and MAT models as a function of the order of the layer (1st-2nd-3rd) and of the dynamical information ratio D. SD ranged between 0.02 and 0.03 spikes/s.

SSN model

MAT model D 1 0.7 0.5 0.4 0.3 0.2 0 1 0.7 0.5 0.4 0.3 0.2 0 1st 6.4 6.1 5.9 5.7 5.4 5.1 4.8 6.6 6.4 6.1 5.7 5.3 4.9 4.5 2nd 6.8 6.6 6.4 6.2 5.8 5.4 5.0 7.7 7.1 6.7 6.4 5.8 5.2 4.5 3rd 7.0 6.9 6.7 6.5 6.1 5.7 5.2 8.9 7.9 7.3 6.9 6. near the axis, due to an internal temporal structure embedded within the model dynamics. In the MAT model it is interesting to observe that with an increase in the order of the layer the attractor contour become clearer even for D as low as D = 0.3. The "reconstructed spike trains" statistics are summarized in Table 2 and the return maps illustrated by Fig. 5b clearly show the noise filtering effect obtained by applying PGA, thus revealing the underlying attractor contour. With a decrease of D, the number of spikes detected by PGA decreased (i.e. the firing rate of the reconstructed spike trains decreased). In the case of SSN significant amount of spikes were detected by PGA even for D ≤ 0.3, but the return maps don't show the contour of the Zaslavskii attractor and the preferred firing sequences detected by PGA can be attributed to the intrinsic dynamics of the model. On the opposite, the MAT model seldom introduced a temporal structure in the output spike train due to intrinsic model dynamics. With the MAT model notice that the similarity ratio and the firing rate of the reconstructed spike train increased from the 1st to the higher order layers with D = 0.4. In both models, the Fano factor was larger for small values of D and became less than 1 at the third layer for both models with D ≥ 0.4. Looking at the similarity ratio the two models behaved very differently. Furthermore, for the MAT model only the similarity ratio tended to be preserved across the layers for D ≥ 0.5 and was even near 0.5 in the 3rd layer with d = 0.4.

Discussion

The deterministic sequence of spikes generated by a chaotic attractor was distributed and embedded in the input spike trains fed to a partially convergent/divergent feed-forward layered network. We have provided evidence that a multiple-timescale adaptive threshold (MAT) neuronal model [START_REF] Kobayashi | Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold[END_REF] was able to retain and transmit a sizable amount of the initial temporal information up to the 3rd layer with dynamical information ratio as low as D = 0.4. Conversely, a simple spiking neuron (SSN) model [START_REF] Izhikevich | Which model to use for cortical spiking neurons?[END_REF] introduced a bias in the temporal pattern of the output spike train associated to its model dynamics which interfered with the input temporally organized information. It is interesting to notice that by passing through the successive layers, the similarity ratio of the SSN neurons decreased drastically despite the fact the reconstructed spike train and the Fano factor were kept rather high.

The current study does not pretend to exclude SSN models from being able to preserve and transmit temporal information through complex neural network circuits because we did not carry out a parameter search of that class of models in order to optimize the performance. The MAT model is interesting because in presence of a pure stochastic input very few spikes were detected by the PGA filtering procedure, thus indicating that this model did not introduce a bias. We consider that this work may be viewed as seminal addressing the novel problem because it suggests that MAT class of models might represent a good candidate for integrating a distributed deterministic temporal information and preserve its dynamics through networks of cell assemblies. Our further work is aimed to determine the limits of this performance by increasing the number of layers, designing inhomogeneous and diverging/converging networks with recurrent connections and with the introduction of explicit synaptic delays and spike timing dependent plasticity.

Fig. 1 .

 1 Fig. 1. Return maps of input spike trains with an average rate of 5 spikes/s as a function of the dynamical information ratio (D). The (n + 1)-st inter-spike-interval are plotted against the n-th inter-spike-interval. The axes are scaled in ms time units. (a) D = 1, (b) D = 0.7, (c) D = 0.5, (d) D = 0.3, (e) D = 0.

Fig. 2 .

 2 Fig. 2. Convergent/divergent feed-forward circuit formed by three neuron layers. Each cell receives 15 afferent spike trains randomly selected out of 20 and a PST (independent Poissonian spike train).

Fig. 3 .

 3 Fig. 3. Return maps of reconstructed spike trains with mean firing rate at 5 spikes/s. (a) from the original Zaslavskii spike train; (b) from a Poissonian spike train

Fig. 4 .

 4 Fig. 4. Left and right panels shows data from two neurons belonging to the 1st layer. Dynamics of the membrane potentials for model SSN (first row) and for model MAT (second row), and the total post-synaptic input current (third row) are shown as a function of the input spike trains (bottom panels) where Zaslavskii and Poissonian spike trains are sorted out. The dynamical information ratio was set to D = 0.5.

Fig. 5 .

 5 Fig.5. Return maps of neuronal output spike trains and spike trains reconstructed from them. One neuron from each layer was selected as an example for several dynamical information ratio D and for both of the SSN and MAT models.

Table 2 .

 2 Firing rate statistics of the reconstructed spike trains of SSN and MAT neurons shown in

Acknowledgments

This study was partially funded by the bi-national JSPS/INSERM grant SYR-NAN and Japan-France Research Cooperative Program.


Jérémie Cabessa
email: jcabessa@nhrg.org

Alessandro E P Villa
email: avilla@nhrg.org

A Hierarchical Classification of First-Order Recurrent Neural Networks

We provide a refined hierarchical classification of first-order recurrent neural networks made up of McCulloch and Pitts cells. The classification is achieved by first proving the equivalence between the expressive powers of such neural networks and Muller automata, and then translating the Wadge classification theory from the automata-theoretic to the neural network context. The obtained hierarchical classification of neural networks consists of a decidable pre-well ordering of width 2 and height ω ω , and a decidability procedure of this hierarchy is provided. Notably, this classification is shown to be intimately related to the attractive properties of the networks, and hence provides a new refined measurement of the computational power of these networks in terms of their attractive behaviours.

Introduction

In neural computability, the issue of the computational power of neural networks has often been approached from the automata-theoretic perspective. In this context, McCulloch and Pitts, Kleene, and Minsky already early proved that the class of first-order recurrent neural networks discloses equivalent computational capabilities as classical finite state automata [START_REF] Kleene | Representation of events in nerve nets and finite automata[END_REF][START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF][START_REF] Minsky | Computation: finite and infinite machines[END_REF]. Later, Kremer extended this result to the class of Elman-style recurrent neural nets, and Sperduti discussed the computational power of different other architecturally constrained classes of networks [START_REF] Kremer | On the computational power of elman-style recurrent networks[END_REF][START_REF] Sperduti | On the computational power of recurrent neural networks for structures[END_REF].

Besides, the computational power of first-order recurrent neural networks was also proved to intimately depend on both the choice of the activation function of the neurons as well as the nature of the synaptic weights under consideration. Indeed, Siegelmann and Sontag showed that, assuming rational synaptic weights, but considering a saturated-linear sigmoidal instead of a hard-threshold activation function drastically increases the computational power of the networks from finite state automata up to Turing capabilities [START_REF] Siegelmann | Turing computability with neural nets[END_REF][START_REF] Siegelmann | On the computational power of neural nets[END_REF]. In addition, Siegelmann and Sontag also nicely proved that real-weighted networks provided with a saturated-linear sigmoidal activation function reveal computational capabilities beyond the Turing limits [START_REF] Siegelmann | Computation beyond the turing limit[END_REF][START_REF] Siegelmann | Neural and super-turing computing[END_REF][START_REF] Siegelmann | Analog computation via neural networks[END_REF].

This paper concerns a more refined characterization of the computational power of neural nets. More precisely, we restrict our attention to the simple class of rational-weighted first-order recurrent neural networks made up of Mc-Culloch and Pitts cells, and provide a refined classification of the networks of this class. The classification is achieved by first proving the equivalence between the expressive powers of such neural networks and Muller automata, and then translating the Wadge classification theory from the automata-theoretic to the neural network context [START_REF] Duparc | A hierarchy of deterministic context-free ω-languages[END_REF][START_REF] Finkel | An effective extension of the wagner hierarchy to blind counter automata[END_REF][START_REF] Selivanov | Wadge degrees of ω-languages of deterministic Turing machines[END_REF][START_REF] Wagner | On ω-regular sets[END_REF]. The obtained hierarchical classification of neural networks consists of a decidable pre-well ordering of width 2 and height ω ω , and a decidability procedure of this hierarchy is provided. Notably, this classification is shown to be intimately related to the attractive properties of the considered networks, and hence provides a new refined measurement of the computational capabilities of these networks in terms of their attractive behaviours.

The Model

In this work, we focus on synchronous discrete-time first-order recurrent neural networks made up of classical McCulloch and Pitts cells.

Definition 1. A first-order recurrent neural network consists of a tuple N = (X, U, a, b, c), where X = {x i : 1 ≤ i ≤ N } is a finite set of N activation cells, U = {u i : 1 ≤ i ≤ M } is a finite set of M external input cells, and a ∈ Q N ×N , b ∈ Q N ×M ,
and c ∈ Q N ×1 are rational matrices describing the weights of the synaptic connections between cells as well as the incoming background activity.

The activation value of cells x j and u j at time t, respectively denoted by x j (t) and u j (t), is a boolean value equal to 1 if the corresponding cell is firing at time t and to 0 otherwise. Given the activation values x j (t) and u j (t), the value x i (t + 1) is then updated by the following equation

x i (t + 1) = σ   N j=1 a i,j • x j (t) + M j=1 b i,j • u j (t) + c i   , i = 1, . . . , N (1)
where σ is the classical hard-threshold activation function defined by σ(α) = 1 if α ≥ 1 and σ(α) = 0 otherwise. Note that Equation (1) ensures that the whole dynamics of network N is described by the following governing equation

x(t + 1) = σ (a • x(t) + b • u(t) + c) , (2)
where x(t) = (x 1 (t), . . . , x N (t)) and u(t) = (u 1 (t), . . . , u M (t)) are boolean vectors describing the spiking configuration of the activation and input cells, and σ denotes the classical hard threshold activation function applied component by component. An example of such a network is given below.

Example 1. Consider the network N depicted in Figure 1. The dynamics of this network is then governed by the following equation:

x1(t+1) x2(t+1) x3(t+1) = σ 0 0 0 1 2 0 0 1 2 0 0 • x1(t) x2(t) x3(t) + 1 0 0 0 0 1 2 • u1(t) u2(t) + 0 1 2 0 1/2 1/2 x 3 x 2 x 1 u 1 u 2 1 1/2 1/2
Fig. 1. A simple neural network

Attractors

The dynamics of recurrent neural networks made of neurons with two states of activity can implement an associative memory that is rather biological in its details [START_REF] Hopfield | unlearning' has a stabilizing effect in collective memories[END_REF]. In the Hopfield framework, stable equilibrium reached by the network that do not represent any valid configuration of the optimization problem are referred to as spurious attractors. According to Hopfield et al., spurious modes can disappear by "unlearning" [START_REF] Hopfield | unlearning' has a stabilizing effect in collective memories[END_REF], but Tsuda et al. have shown that rational successive memory recall can actually be implemented by triggering spurious modes [START_REF] Tsuda | Memory dynamics in asynchronous neural networks[END_REF]. Here, the notions of attractors, meaningful attractors, and spurious attractors are reformulated in our precise context. Networks will then be classified according to their ability to switch between different types of attractive behaviours. For this purpose, the following definitions need to be introduced. As preliminary notations, for any k > 0, we let the space of k-dimensional boolean vectors be denoted by B k , and we let the space of all infinite sequences of k-dimensional boolean vectors be denoted by [B k ] ω . Moreover, for any finite sequence of boolean vectors v, we let the expression v ω = vvvv • • • denote the infinite sequence obtained by infinitely many concatenations of v.

Now, let N be some network with N activation cells and M input cells. For each time step t ≥ 0, the boolean vectors x(t) = (x 1 (t), . . . , x N (t)) ∈ B N and u(t) = (u 1 (t), . . . , u M (t)) ∈ B M describing the spiking configurations of both the activation and input cells of N at time t are respectively called the state of N at time t and the input submitted to N at time t. An input stream of N is then defined as an infinite sequence of consecutive inputs s = (u(i)

) i∈N = u(0)u(1)u(2) • • • ∈ [B M ] ω .
Moreover, assuming the initial state of the network to be x(0) = 0, any input stream s = (u(i)

) i∈N = u(0)u(1)u(2) • • • ∈ [B M ] ω induces via Equation (2) an infinite sequence of consecutive states e s = (x(i)) i∈N = x(0)x(1)x(2) • • • ∈ [B N ] ω that is called the evolution of N induced by the input stream s.
Along some evolution e s = x(0)x(1)x(2) • • • , irrespective of the fact that this sequence is periodic or not, some state will repeat finitely often whereas other will repeat infinitely often. The (finite) set of states occurring infinitely often in the sequence e s is denoted by inf(e s ). It can be observed that, for any evolution e s , there exists a time step k after which the evolution e s will necessarily remain confined in the set of states inf(e s ), or in other words, there exists an index k such that x(i) ∈ inf(e s ) for all i ≥ k. However, along evolution e s , the recurrent visiting of states in inf(e s ) after time step k does not necessarily occur in a periodic manner. Now, given some network N with N activation cells, a set A = {y 0 , . . . , y k } ⊆ B N is called an attractor for N if there exists an input stream s such that the corresponding evolution e s satisfies inf(e s ) = A. Intuitively, an attractor can be seen a trap of states into which some network's evolution could become forever confined. We further assume that attractors can be of two distinct types, namely meaningful or optimal vs. spurious or non-optimal. In this study we do not extend the discussion about the attribution of the attractors to either type. From this point onwards, we assume any given network to be provided with the corresponding classification of its attractors into meaningful and spurious types. Now, let N be some network provided with an additional type specification of each of its attractors. The complementary network N is then defined to be the same network as N but with an opposite type specification of its attractors. 3In addition, an input stream s of N is called meaningful if inf(e s ) is a meaningful attractor, and it is called spurious if inf(e s ) is a spurious attractor. The set of all meaningful input streams of N is called the neural language of N and is denoted by L(N ). Note that the definition of the complementary network implies that L(N ) = L(N ) . Finally, an arbitrary set of input streams L ⊆ [B M ] ω is defined as recognizable by some neural network if there exists a network N such that L(N ) = L. All preceding definitions are now illustrated in the next example.

Example 2. Consider again the network N described in Example 1, and suppose that an attractor is meaningful for N if and only if it contains the state (1, 1, 1) T (i.e. where the three activation cells simultaneously fire). The periodic input stream s = [( 11 ) ( 1 1 ) ( 1 1 ) ( 0 0 )] ω induces the corresponding periodic evolution

e s = 0 0 0 1 0 0 1 1 1 1 1 1 0 1 0 1 0 0 ω .
Hence, inf(e s ) = {(1, 1, 1) T , (0, 1, 0) T , (1, 0, 0) T }, and the evolution e s of N remains confined in a cyclic visiting of the states of inf(e s ) already from time step t = 2. Thence, the set {(1, 1, 1) T , (0, 1, 0) T , (1, 0, 0) T } is an attractor of N . Moreover, this attractor is meaningful since it contains the state (1, 1, 1) T .

Recurrent Neural Networks and Muller Automata

In this section, we provide an extension of the classical result stating the equivalence of the computational capabilities of first-order recurrent neural networks and finite state machines [START_REF] Kleene | Representation of events in nerve nets and finite automata[END_REF][START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF][START_REF] Minsky | Computation: finite and infinite machines[END_REF]. More precisely, here, the issue of the expressive power of neural networks is approached from the point of view of the theory of automata on infinite words, and it is proved that first-order recurrent neural networks actually disclose the very same expressive power as finite Muller automata. Towards this purpose, the following definitions first need to be recalled.

A finite Muller automaton is a 5-tuple A = (Q, A, i, δ, T ), where Q is a finite set called the set of states, A is a finite alphabet, i is an element of Q called the initial state, δ is a partial function from Q × A into Q called the transition function, and T ⊆ P(Q) is a set of set of states called the table of the automaton. A finite Muller automaton is generally represented by a directed labelled graph whose nodes and labelled edges respectively represent the states and transitions of the automaton.

Given a finite Muller automaton A = (Q, A, i, δ, T ), every triple (q, a, q ) such that δ(q, a) = q is called a transition of A. A path in A is then a sequence of consecutive transitions ρ = ((q 0 , a 1 , q 1 ), (q 1 , a 2 , q 2 ), (q 2 , a 3 , q 3 ), . . .), also denoted by ρ :

q 0 a1 -→ q 1 a2 -→ q 2 a3 -→ q 3 • • • .
The path ρ is said to successively visit the states q 0 , q 1 , . . .. The state q 0 is called the origin of ρ, the word a 1 a 2 a 3 • • • is the label of ρ, and the path ρ is said to be initial if q 0 = i. If ρ is an infinite path, the set of states visited infinitely often by ρ is denoted by inf(ρ). Besides, a cycle in A consists of a finite set of states c such that there exists a finite path in A with same origin and ending state that visits precisely all the sates of c. A cycle is called successful if it belongs to T , and non-succesful otherwise. Moreover, an infinite initial path ρ of A is called successful if inf(ρ) ∈ T . An infinite word is then said to be recognized by A if it is the label of a successful infinite path in A, and the ω-language recognized by A, denoted by L(A), is defined as the set of all infinite words recognized by A. The class of all ω-languages recognizable by some Muller automata is precisely the class of ω-rational languages. Now, for each ordinal α < ω ω , we introduce the concept of an α-alternating tree in a Muller automaton A, which consists of a tree-like disposition of the successful and non-successful cycles of A induced by the ordinal α (see Figure 2). We first recall that any ordinal 0 < α < ω ω can uniquely be written of the form α = ω np •m p +ω np-1 •m p-1 +. . .+ω n0 •m 0 , for some p ≥ 0, n p > n p-1 > . . . > n 0 ≥ 0, and m i > 0. Then, given some Muller automata A and some ordinal

α = ω np • m p + ω np-1 • m p-1 + . . . + ω n0 • m 0 < ω ω ,
an α-alternating tree (resp. α-co-alternating tree) is a sequence of cycles of A (C i,j k,l ) i≤p,j<2 i ,k<mi,l≤ni such that: firstly, C 0,0 0,0 is successful (resp. not successful); secondly, C i,j k,l C i,j k,l+1 , and

C i,j k,l+1 is successful iff C i,j k,l is not successful; thirdly, C i,j k+1,0 is strictly accessible from C i,j
k,0 , and C i,j k+1,0 is successful iff C i,j k,0 is not successful; fourthly, C i+1,2j 0,0 and C i+1,2j+1 0,0 are both strictly accessible from C i,j mi-1,0 , and each C i+1,2j 0,0 is successful whereas each C i+1,2j+1 0,0 is not successful. An α-alternating tree is said to be maximal in A if there is no β-alternating tree in A such that β > α.

We now come up to the equivalence of the expressive power of recurrent neural networks and Muller automaton. First of all, we prove that any firstorder recurrent neural network can be simulated by some Muller automaton.

Proposition 1. Let N be a network provided with a type specification of its attractors. Then there exists a Muller automaton

A N such that L(N ) = L(A N ). C 0,0 0,n 0 C 0,0 1,n 0 C 0,0 m 0 -1,n 0 . . . . . . . . . C 0,0 0,1 C 0,0 1,1 C 0,0 m 0 -1,1 C 0,0 0,0 -→ C 0,0 1,0 -→ • • • -→ C 0,0 m 0 -1,0 -→ -→ C 1,0 0,n 1 C 1,0 1,n 1 C 1,0 m 1 -1,n 1 . . . . . . . . . C 1,0 0,1 C 1,0 1,1 C 1,0 m 1 -1,1 C 1,0 0,0 -→ C 1,0 1,0 -→ • • • -→ C 1,0 m 1 -1,0 • • • -→ -→ • • • C 1,1 0,n 1 C 1,1 1,n 1 C 1,1 m 1 -1,n 1 . . . . . . . . . C 1,1 0,1 C 1,1 1,1 C 1,1 m 1 -1,1 C 1,1 0,0 -→ C 1,1 1,0 -→ • • • -→ C 1,1 m 1 -1,0 • • • -→ -→ • • • Fig. 2.
The inclusion and accessibility relations between cycles in an α-alternating tree

Proof. Let N be given by the tuple (X, U, a, b, c), with card(X) = N , card(U ) = M , and let the meaningful attractors of N be given by A 1 , . . . , A K . Now, consider the Muller automaton According to the construction given in the proof of Proposition 1, any evolution of the network N naturally induces a corresponding infinite initial path in the Muller automaton A N , and conversely, any infinite initial path in A N corresponds to some possible evolution of N . This observation ensures the existence of a biunivocal correspondence between the attractors of the network N and the cycles in the graph of the corresponding Muller automaton A N . Consequently, a procedure to compute all possible attractors of a given network N is simply obtained by first constructing the corresponding Muller automaton A N and then listing all cycles in the graph of A N .

A N = (Q, A, i, δ, T ), where Q = B N , A = B M , i is the N -dimensional zero vector, δ : Q × A → Q is defined by δ(x, u) = x if and only if x = σ (a • x + b • u + c),
Conversely, we now prove that any Muller automaton can be simulated by some first-order recurrent neural network. For the sake of convenience, we choose to restrict our attention to Muller automata over the binary alphabet B 1 . Proposition 2. Let A be some Muller automaton over the alphabet B 1 . Then there exists a network N A such that L(A) = L(N A ).

Proof. Let A be given by the tuple (Q, A, q 1 , δ, T ), with Q = {q 1 , . . . , q N } and T ⊆ P(Q). Now, consider the network N A = (X, U, a, b, c) defined as follows: First of all, X = {x i : 1 ≤ i ≤ 2N } ∪ {x 1 , x 2 , x 3 , x 4 }, U = {u 1 }, and each state q i in the automaton A gives rise to a two cell layer {x i , x N +i } in the network N A as illustrated in Figure 3. Moreover, the synaptic weights between u 1 and all activation cells, between all cells in {x 1 , x 2 , x 3 , x 4 }, as well as the background activity are precisely as depicted in Figure 3. Furthermore, for each 1 ≤ i ≤ N , both cells x i and x N +i receive a weighted connection of intensity 1 2 from cell x 4 (resp. x 2 ) if and only if δ(q 1 , (0)) = q i (resp. δ(q 1 , (1)) = q i ), as also shown in Figure 3. Farther, for each 1 ≤ i, j ≤ N , there exist two weighted connection of intensity 1 2 from cell x i (resp. from cell x N +i ) to both cell x j and x N +j if and only if δ(q i , (1)) = q j (resp. δ(q i , (0)) = q j ), as partially illustrated in Figure 3 only for the k-th layer. This description of the network N A ensures that, for any possible evolution of N A , the two cells x 1 and x 3 are firing at each time step t ≥ 1, and furthermore, one and only one cell of {x i : 1 ≤ i ≤ 2N } are firing at each time step t ≥ 2. According to this observation, for any 1 ≤ j ≤ N , let 1 j ∈ B 2N +4 (resp. 1 N +j ∈ B 2N +4 ) denote the boolean vector describing the spiking configuration where only the cells x 1 , x 3 , and x j (resp. x 1 , x 3 , and x N +j ) are firing. Hence, any evolution x(0)x( 1

)x(2) • • • of N A satisfies x(t) ∈ {1 k : 1 ≤ k ≤ N } ∪ {1 N +l : 1 ≤ l ≤ N } for all t ≥ 2,

and thus any attractor A of N can necessarily be written of the form

A = {1 k : k ∈ K} ∪ {1 N +l : l ∈ L}, for some K, L ⊆ {1, 2, . . . , N }. Now, any infinite sequence s = u(0)u(1)u(2) • • • ∈ [B 1 ] ω induces both a corresponding infinite path ρ s : q 1 u(0) ---→ q j1 u(1) ---→ q j2 u(2) ---→ q j3 • • • in A as well as a corresponding evolution e s = x(0)x(1)x(2) • • • in N A .
The network N A is then related to the automaton A via the following important property: for each time step t ≥ 1, if u(t) = (1), then x(t + 1) = 1 jt , and if u(t) = (0), then x(t + 1) = 1 N +jt . In other words, the infinite path ρ s and the evolution e s evolve in parallel and satisfy the property that the cell x j is spiking in N A if and only if the automaton A is in state q j and reads letter (1), and the cell x N +j is spiking in N A if and only if the automaton A is in state q j and reads letter (0). Finally, an attractor

u 1 x 1 x 2 x 3 x 4 -1/2 -1 -1 -1 1/2 1/2 +1 +1 +1 1/2 1/2 1/2 1/2 x 1 x i x N 1/2 1/2 x 2N x j x k x N+1 +1 1/2 1/
L(A) = L(N A ).
Finally, the following example provides an illustration of the two translating procedures described in the proofs of propositions 1 and 2.

Example 3. The translation from the network N described in Example 2 to its corresponding Muller automaton A N is illustrated in Figure 4. Proposition 1 ensures that L(N ) = L(A N ). Conversely, the translation from some given Muller automaton A over the alphabet B 1 to its corresponding network N A is illustrated in Figure 5. Proposition 2 ensures that L(A) = L(N A ).

1/2 1/2 x 3 x 2 x 1 u 1 u 2 1 1/2 1/2 0 0 0 0 1 0 1 1 0 1 1 1 0 1 1 1 0 0 ( 1 0 ) ( 0 0 ) ( 1 1 )
( 0 0 ) ( 0 1 ) ( 1 0 )( 1 1 ) ( 0 0 ) ( 0 1 ) ( 1 0 ) ( 1 1 ) ( 0 0 ) ( 1 0 )
(

A ⊆ B 3 is meaningful for N Table T = {A ∈ B 3 : A is meaningful for N } if and only if (1, 1, 1) T ∈ A Fig. 4. Translation from a given network N provided with a type specification of its attractors to a corresponding Muller automaton AN

q 1 q 2 q 3 (1) (0) (1)
(1)

(0) (0) u 1 x 3 x 4 x 5 x 1 x 2 x 3 x 4 -1/2 -1 -1 -1 1 1/2 1/2 1/2 1/2 1/2 1/2 +1 +1 +1 1/2 1/2 1/2 1/2 x 6 x 1 x 2 Table T = {{q2}, {q3}}
Meaningful attractors: A1 = {15} and A2 = {13}. In the theory of automata on infinite words, abstract machines are commonly classified according the topological complexity of their underlying ω-language, as for instance in [START_REF] Duparc | A hierarchy of deterministic context-free ω-languages[END_REF][START_REF] Finkel | An effective extension of the wagner hierarchy to blind counter automata[END_REF][START_REF] Selivanov | Wadge degrees of ω-languages of deterministic Turing machines[END_REF][START_REF] Wagner | On ω-regular sets[END_REF]. Here, this approach is translated from the automata to the neural networks context, in order to obtain a refined classification of firstorder recurrent neural networks. Notably, the obtained classification actually refers to the ability of the networks to switch between meaningful and spurious attractive behaviours. For this purpose, the following facts and definitions need to be introduced. To begin with, for any k > 0, the space [B k ] ω can naturally be equipped with the product topology of the discrete topology over B k . Thence, a function f :

[B k ] ω → [B l ] ω is
[B M1 ] ω → [B M2 ] ω such that any input stream s of N 1 satisfies s ∈ L(N 1 ) ⇔ f (s) ∈ L(N 2 ).
The corresponding strict reduction, equivalence relation, and incomparability relation are then naturally defined by

N 1 < W N 2 iff N 1 ≤ W N 2 ≤ W N 1 , as well as N 1 ≡ W N 2 iff N 1 ≤ W N 2 ≤ W N 1 , and N 1 ⊥ W N 2 iff N 1 ≤ W N 2 ≤ W N 1 . Moreover, a network N is called self-dual if N ≡ W N ; it is non-self-dual if N ≡ W N ,
which can be proved to be equivalent to saying that N ⊥ W N . By extension, an ≡ W -equivalence class of networks is called self-dual if all its elements are self-dual, and non-self-dual if all its elements are non-self-dual. Now, the Wadge reduction over the class of neural networks naturally induces a hierarchical classification of networks. Formally, the collection of all first-order recurrent neural networks ordered by the Wadge reduction "≤ W " is called the RNN hierarchy.

Propositions 1 and 2 ensure that the RNN hierarchy and the Wagner hierarchy -the collection of all ω-rational languages ordered by the Wadge reduction [START_REF] Wagner | On ω-regular sets[END_REF] -coincide up to Wadge equivalence. Accordingly, a precise description of the RNN hierarchy can therefore be given as follows. First of all, the RNN hierarchy is well founded, i.e. there is no infinite strictly descending sequence of networks N 0 > W N 1 > W N 2 > W . . .. Moreover, the maximal strict chains in the RNN hierarchy have length ω ω , meaning that the RNN hierarchy has a height of ω ω . Furthermore, the maximal antichains of the RNN hierarchy have length 2, meaning that the RNN hierarchy has a width of 2. 4 More precisely, any two networks N 1 and N 2 satisfy the incomparability relation N 1 ⊥ W N 2 if and only if N 1 and N 2 are non-self-dual networks such that N 1 ≡ W N 2 . These properties imply that, up to Wadge equivalence and complementation, the RNN hierarchy is actually a well-ordering. In fact, the RNN hierarchy consists of an alternating succession of non-self-dual and self-dual classes with pairs of non-self-dual classes at each limit level, as illustrated in Figure 6, where circle represent the Wadge equivalence classes of networks and arrows between circles represent the strict Wadge reduction between all elements of the corresponding classes. For convenience reasons, the degree of a network N in the RNN hierarchy is now defined in order to make the non-self-dual (n.s.d.) networks and the self-dual ones located just one level above share the same degree, as illustrated in Figure 6:

d(N ) =      1 if L(N ) = ∅ or ∅ , sup {d(M) + 1 : M n.s.d. and M < W N } if N is non-self-dual, sup {d(M) : M n.s.d. and M < W N } if N is self-dual.
Also, the equivalence between the Wagner and RNN hierarchies ensure that the RNN hierarchy is actually decidable, in the sense that there exists a algorithmic procedure computing the degree of any network in the RNN hierarchy. All the above properties of the RNN hierarchy are summarized in the following result.

Theorem 1. The RNN hierarchy is a decidable pre-well-ordering of width 2 and height ω ω .

Proof. The Wagner hierarchy consists of a decidable pre-well-ordering of width 2 and height ω ω [START_REF] Wagner | On ω-regular sets[END_REF]. Propositions 1 and 2 ensure that the RNN hierarchy and Wagner hierarchy coincide up to Wadge equivalence. The following result provides a detailed description of the decidability procedure of the RNN hierarchy. More precisely, it is shown that the degree of a network N in the RNN hierarchy corresponds precisely to the largest ordinal α such that there exists an α-alternating tree or an α-co-alternating tree in the Muller automaton A N . Theorem 2. Let N be a network provided with a type specification of its attractors, A N be the corresponding Muller automaton of N , and α be an ordinal such that 0 < α < ω ω .

• If there exists in A N a maximal α-alternating tree and no maximal α-coalternating tree, then d(N ) = α and N is non-self-dual.

• If there exists in A N a maximal α-co-alternating tree and no maximal αalternating tree, then d(N ) = α and N is non-self-dual. • If there exist in A N both a maximal α-alternating tree as well as a maximal α-co-alternating tree, then d(N ) = α and N is self-dual.

Proof. For any ω-rational language L, let d W (L) denote the degree of L in the Wagner hierarchy. On the one hand, propositions 1 and 2 ensure that d(N ) = d W (L(A N )). On the other hand, the decidability procedure of the Wagner hierarchy states that d W (L(A N )) corresponds precisely to the largest ordinal α such that there exists a maximal α-(co)-alternating tree in A N [START_REF] Wagner | On ω-regular sets[END_REF].

The decidability procedure of the degree of a network N in the the RNN hierarchy thus consists in first translating the network N into its corresponding Muller automaton A N (as described in Proposition 1), and then returning the ordinal α associated to the maximal α-(co)-alternating tree(s) in contained in A N (which can be achieved by some graph analysis of the automaton A N ). In other words, the complexity of a network N is directly related to the relative disposition of the successful and non-successful cycles in the Muller automaton A N , or in other words, to how some infinite path in A N could maximally alternate between successful and non-successful cycles along its evolution. Therefore, according to the biunivocal correspondence between cycles in A N and attractors of N , as well as between infinite paths in A N and evolutions of the network N , it follows that the complexity of a network N in the RNN hierarchy actually refers to the capacity of this network to maximally alternate between punctual visitings of meaningful and spurious attractors along some possible evolutiona concept close to chaotic itinerancy [START_REF] Tsuda | Chaotic itinerancy as a dynamical basis of hermeneutics of brain and mind[END_REF][START_REF] Kaneko | Chaotic itinerancy[END_REF].

Example 4. Let N be the network of Example 2. Then d(N ) = ω and N is non-self-dual. Indeed, {(0, 0, 0) T } {(0, 0, 0) T , (1, 0, 0) T , (1, 1, 1) T , (0, 1, 1) T } is a maximal ω 1 -co-alternating tree in the Muller automaton A N of Figure 4.

Conclusion

The present work proposes a new approach of neural computability from the point of view infinite word reading automata theory. More precisely, the Wadge classification of infinite word languages is translated from the automata-theoretic to the neural network context, and a transfinite decidable hierarchical classification of first-order recurrent neural network is obtained. This classification provides a better understanding of this simple class of neural networks that could be relevant for implementation issues. Moreover, the Wadge hierarchies of deterministic pushdown automata or deterministic Turing Machines both with Muller conditions [START_REF] Duparc | A hierarchy of deterministic context-free ω-languages[END_REF][START_REF] Selivanov | Wadge degrees of ω-languages of deterministic Turing machines[END_REF] ensure that such Wadge-like classifications of strictly more powerful models of neural networks could also be described; however, in these cases, the decidability procedures of the obtained hierarchies remain hard open problems.

Besides, this work is envisioned to be extended in several directions. First of all, it could be of interest to study the same kind of hierarchical classification applied to more biologically oriented models, like neural networks provided with some additional simple STDP rule. In addition, neural networks' computational capabilities should also rather be approached from the point of view of finite word instead of infinite word reading automata, as for instance in [START_REF] Kremer | On the computational power of elman-style recurrent networks[END_REF][START_REF] Siegelmann | Computation beyond the turing limit[END_REF][START_REF] Siegelmann | Neural and super-turing computing[END_REF][START_REF] Siegelmann | Turing computability with neural nets[END_REF][START_REF] Siegelmann | Analog computation via neural networks[END_REF][START_REF] Siegelmann | On the computational power of neural nets[END_REF][START_REF] Sperduti | On the computational power of recurrent neural networks for structures[END_REF]. Unfortunately, as opposed to the case of infinite words, the classification theory of finite words reading machines is still a widely undeveloped, yet promising issue. Finally, the study of hierarchical classifications of neural networks induced by more biologically oriented reduction relations than the Wadge reduction would be of specific interest.

  and T = {A 1 , . . . , A K }. According to this construction, any input stream s of N is meaningful for N if and only if s is recognized by A N . In other words, s ∈ L(N ) if and only if s ∈ L(A N ), and therefore L(N ) = L(A N ).

Fig. 5 .

 5 Fig. 5. Translation from a given Muller automaton A to a corresponding network NA provided with a type specification of its attractors

1 degree ω • 2 Fig. 6 .

 126 Fig. 6. The RNN hierarchy

  2 K} ∪ {1 N +l : l ∈ L} with K, L ⊆ {1, 2, . . . , N } is set to be meaningful if and only if {q k : k ∈ K} ∪{ q l : l ∈ L} ∈ T . Consequently, for any infinite infinite sequence s ∈ [B 1 ] ω , the infinite path ρ s in A satisfies inf(ρ s ) ∈ T if and only if the evolution e s in N A is such that inf(e s ) is a meaningful attractor. Therefore,

	...		...			...
	...	x N+i	...	x N+j	...	x N+k

Fig. 3. The network NA A = {1 k : k ∈

  said to be continuous if and only if the inverse image by f of every open set of [B l ] ω is an open set of [B k ] ω . Now, given two first-order recurrent neural networks N 1 and N 2 with M 1 and M 2 input cells respectively, we say that N 1 Wadge reduces [17] (or continuously reduces or simply reduces) to N 2 , denoted by N 1 ≤ W N 2 , if any only if there exists a continuous function f :

More precisely, A is a meaningful attractor for N if and only if A is a spurious attractor for N .

A strict chain (resp. an antichain) in the RNN hierarchy is a sequence of neural networks (N k ) k∈α such that Ni <W Nj iff i < j (resp. such that Ni ⊥W Nj for all i, j ∈ α with i = j). A strict chain (resp. an antichain) is said to be maximal if its length is at least as large as the length of every other strict chain (resp. antichain).


Camille Salaün
email: camille.salaun@upmc.fr

Vincent Padois
email: vincent.padois@upmc.fr

Olivier Sigaud
email: olivier.sigaud@upmc.fr

A Two-Level Model of Anticipation-Based Motor Learning for Whole Body Motion

We present a model of motor learning based on a combination of Operational Space Control and Optimal Control. Anticipatory processes are used both in the learning of the dynamics model of the system and in the coordination between both types of control. In order to illustrate the proposed model and associated control method, we apply these principles to the control of a simplified virtual humanoid performing a stand-up task starting from a crouching posture.

Introduction

Early in the history of motor control studies, Bernstein raised the following paradox: thanks to the redundancy of their motor system, human beings can realize a task with an infinity of ways. However, for a given task, they always reproduce the same kind of stereotypic motion. The problem consists in explaining where these regularities come from, given the diversity of possible solutions [START_REF] Bernstein | The Co-ordination and Regulation of Movements[END_REF].

An attractive solution to this problem consist in stating that motor control optimises some criterion, thus the performed motion among the many possible ones is the optimal one with respect to that particular criterion. Several potential criteria have been proposed such as the Minimum jerk [START_REF] Flash | The Coordination of Arm Movements: An Experimentally Confirmed Mathematical Model[END_REF], the Minimum torque-change [START_REF] Uno | Formation and control of optimal trajectory in human multijoint arm movement[END_REF] and their variants (Minimum motor command change [START_REF] Kawato | Optimization and learning in neural networks for formation and control of coordinated movement[END_REF] and Minimum commanded torque change [START_REF] Nakano | Quantitative examinations of internal representations for arm trajectory planning: Minimum commanded torque change model[END_REF]). These criteria are all based on the idea that human motion must optimise smoothness, but a principled explanation of why motion should optimise smoothness is missing.

By contrast, the Minimum end-point variance criterion comes with a convincing explanation: natural movements must reach their target accurately, in spite of a neural noise that is proportional to the motor signal. Thus motor control must minimise the error on the terminal position of the controlled effector. Harris and Wolpert [START_REF] Harris | Signal-dependent noise determines motor planning[END_REF] showed that an optimal controller based on this principle can reproduce all motor invariants on which previous criteria were based, but also additional invariants such as Fitts' law [START_REF] Fitts | The information capacity of the human motor system in controlling the amplitude of movement[END_REF] which relates accuracy to the speed of execution.

The direct consequence of this criterion is the Minimum intervention principle [START_REF] Todorov | A minimal intervention principle for coordinated movement[END_REF] that stipulates that the motor system activates muscles as few as possible to achieve its task, so as to minimise the motor signal related noise. As a result, the generated control only acts in the dimensions that are relevant with respect to the task and rejects the noise to the other dimensions. This fact is validated experimentally by [START_REF] Scholz | The uncontrolled manifold concept: identifying control variables for a functional task[END_REF] who showed that fluctuations are larger on the Uncontrolled manifold consisting of the irrelevant dimensions than on the relevant ones.

Fig. 1. A global view of SOFC: the controller must estimate the state of the system and maximise movement accuracy in the presence of sensory and motor noise Based on these findings, the Stochastic Optimal Feedback Control (SOFC) method [START_REF] Todorov | Optimal feedback control as a theory of motor coordination[END_REF][START_REF] Todorov | Optimality principles in sensorimotor control[END_REF], illustrated in Fig. 1, and some variants (Terminal Optimal Feedback Controller [START_REF] Guigon | Computational motor control: Redudancy and invariance[END_REF][START_REF] Guigon | Optimality, stochasticity and variability in motor behavior[END_REF][START_REF] Guigon | Computational motor control: Feedback and accuracy[END_REF] and Task Optimisation in the Presence of Signal dependent noise [START_REF] Miyamoto | Computing the optimal trajectory of arm movement: the TOPS (task optimization in the presence of signal-dependent noise) model[END_REF]) have received a wide agreement as good models of human motor control for elementary tasks for a single arm. However the computational cost of Optimal Control (OC) methods such as SOFC make them unsuitable to solve larger problems such as the synthesis of whole body motion. Thus an important issue in this domain consists in finding computationally less expensive approaches endowed with the same properties.

Another issue is concerned about learning and adapting to changing dynamic conditions. Several authors [START_REF] Wolpert | Computational principles of movement neuroscience[END_REF][START_REF] Wolpert | Multiple paired forward and inverse models for motor control[END_REF] propose a general framework in which human motor control is model-based, combines feedforward and feedback control processes [START_REF] Davidson | Widespread access to predictive models in the motor system: a short review[END_REF] and calls upon optimisation processes as explained above. In this framework, motor adaptation results from learning the dynamics model of the system.

These principles were first implemented in Modular Selection And Identification for Control (MOSAIC) [START_REF] Haruno | MOSAIC model for sensorimotor learning and control[END_REF], a modular architecture where each module gets specialised to deal with particular dynamic circumstances. In MOSAIC, control learning is based on an handcrafted corrector that requires a priori knowledge on the dynamics of the system. Multiple Model-Based Reinforcement Learning (MMRL) [START_REF] Doya | Multiple model-based reinforcement learning[END_REF] solves this limitation by introducing a reinforcement learning process. Each module in MMRL calls upon an optimal control module, thus it still suffers from the computational cost of OC approaches.

One way to overcome this computational cost problem is suggested by Shadmehr's global view of the motor system [START_REF] Shadmehr | The Computational Neurobiology of Reaching and Pointing[END_REF]. From studies based on reaching and pointing movements, Shadmehr considers that motion generation can be decomposed into two parts: at the higher level, the central nervous system computes, in a visual frame of reference centred on the fixation point, a vector from the end-point effector to the target. This vector defines a trajectory in the operational visual space. Then this trajectory is converted into a muscular, low-level control of the joints angles by calling upon learned visuo-articular velocity and position mappings.

The most adequate robotics control tool to formalise Shadmehr's intuition consists in projection methods, such as Operational Space Control (OSC) [START_REF] Khatib | A unified approach for motion and force control of robot manipulators: The operational space formulation[END_REF]. These methods perform their computation in the so-called operational space, a space relative to the task, which is usually smaller than the joints space. They can thus be applied to large robotic systems [START_REF] Sentis | Control of free-floating humanoid robots through task prioritization[END_REF]. They also give rise to a mathematically straightforward way of decoupling a set of tasks ranked by priority [START_REF] Chiaverini | Singularity-robust task-priority redundancy resolution for real-time kinematic control of robot manipulators[END_REF], but are sensitive to these priorities and cannot directly perform a global optimisation of a cost function associated to the motion of the system over time.

By contrast, optimisation methods, such as OC perform their computation in the joints space, but benefit from the easy definition of the constraints and performance criteria [START_REF] Barthlemy | Stability measure of postural dynamic equilibrium based on residual radius[END_REF] of the tasks. However, the only way to express priorities among tasks within this framework consists in tuning the relative weights of the corresponding performance criteria.

In this contribution, we present an hierarchical model combining the assets of both control methods. Our model is consistent with Shadmehr's view in that a high level control loop computes a global trajectory in the operational space and then a low level control loop dynamically realizes it at the joints level. Besides, our model calls upon a velocity kinematics model and a dynamics model of the system, respectively. We show how these models can be learned from experience within an anticipatory process. Finally, the coordination of both control levels also results from an anticipatory process, the control set point being chosen as a point that is forward in time with respect to the current time. Our model is functional rather than neuromimetic: It focuses on the level of computational principles of human motor control and does not claim anything about the underlying neuroanatomy.

The contribution is organised as follows. Section 2 gives the technical background of each component of our model. The global structure of the controller itself is presented in section 3. In section 4, we describe an empirical study about the control of a simplified virtual humanoid performing a stand-up task starting from a crouching posture. We discuss the results in section 5.

Background

We present a model of human motor learning that combines the definition of tasks in the operational space, the use of OC at the joints level and experiencebased, incremental model learning processes. An inverse velocity kinematics model is learned with the Locally Weighted Projection Regression (lwpr) nonlinear function approximation algorithm [START_REF] Vijayakumar | LWPR: A scalable method for incremental online learning in high dimensions[END_REF] whereas a direct linearised model of the dynamics is learned with a classical Recursive Least Squares algorithm. Learning is performed on-line, during the control of the system. The different components of the resulting adaptive control architecture are described below.

Forward and Inverse Velocity Kinematics

The operational space or task space is the most natural space to describe the motion of the end effector of a system, whereas the joints space is the most natural space to control the system. The OSC approach consists in specifying the goal in the operational space and then using a linear transformation from the operational space to the joints space to control the system. Let ξ be an m-dimensional vector of operational coordinates used to describe the task (e.g. position and orientation of the end-point effector) and let q be an n-dimensional vector of generalised coordinates -i.e. the vector of parameters necessary to describe the configuration of the system without ambiguity. For systems with a fixed base, the generalised coordinates are often chosen as the joints angles. The relation ξ = f (q) is called forward kinematics model (FKM). It describes the operational coordinates as a function of the configuration of the system. Differentiating with respect to time, we obtain

ξ = J(q) q, (1)
where J(q) = ∂f (q) ∂q is the Jacobian matrix of f . To control the system, we need q as a function of ξ. For a redundant system (rank(J) < n) for which the inversion of (1) has an infinity of solutions, one particular solution is:

q = J(q) ξ (2)
where J(q) is a weighted pseudoinverse of J(q) [START_REF] Golub | Matrix Computations[END_REF].

Weighted pseudo-inverse of J (q) are written:

J(q) = W -1 J (q) T J (q) W -1 J (q) T -1
where W is a symmetric and positive matrix of dimension n. This particular set of solutions to (1) minimises the Euclidean W-weighted norm qT W q of the solution. In the case where W = I n , this solution is called the Moore-Penrose or pseudoinverse of the Jacobian matrix and is noted J + .

The general form of the solution can be written as:

q = J(q) ξ + I n -J (q) J (q) q0 , ( 3
)
where I n -J (q) J (q) is a projector onto the nullspace of J (q) (i.e. the space of internal mobility with respect to the main task) and q0 any vector of dimension n. The second term of the right-hand term of equation ( 3) gives access to the redundancy of the system in a hierarchical manner: any secondary task or constraint projected onto the nullspace of the Jacobian matrix is given less priority than the main task. It is achieved as long as it does not disturb this main task.

Operational Space Control

OSC consists in reaching a target ξ ref from an initial position. This control method uses the inverse velocity kinematics model described by equation ( 2) as well as a proportional controller to compute the desired operational velocity:

ξdes = K P (ξ des -ξ) , ( 4
)
where K p is a positive definite matrix used as a proportional gain. This method is called resolved rate motion control. Using the inverse velocity kinematics model described by equation ( 2), we have: qdes = J(q) ξdes .

(

Thus, given a desired task velocity ξ ref and knowing J(q) , one can derive the desired joints velocities ξdes . If one wishes to specify the desired joints accelerations, one may use a progressive differentiation:

qdes = qdes - q Δt (6)
Independently of the formalism used (Lagrange or Newton-Euler), the equations of motion of a fully-actuated, holonomic system in the contact-free case can be written: Γ = M (q)q + b(q, q) + g(q) [START_REF] Fitts | The information capacity of the human motor system in controlling the amplitude of movement[END_REF] where M is the symmetric positive definite inertia matrix of the system, b is the vector of nonlinear effects modelling centrifugal and Coriolis forces, g is the gravity vector and Γ represents the torques applied to the system. This inverse dynamics model is used to compute, for a given state, the torque that must be applied to the system to get the desired joint accelerations:

Γ = M (q)q des + b(q, q) + g(q) ( 8
)
where qdes is computed from ξdes using equations ( 4), ( 5) and (6).

Linearisation of the Dynamics Model

It is possible to linearise the dynamics model around a state o:

F 1 o q + F 2 o q + F 3 o q = Γ ,
We can reformulate this as:

q q = 0 Id -F -1 1o F 3 o -F -1 1o F 2 o q q + 0 F -1 1 o Γ .
This formula is called state representation and can be rewritten ẋ = Ax+Bu. The state and action matrices A and B generally depend on the state. This will not be the case in the study hereafter, the dynamics model being linearised. In its discrete form, this equation is written:

x k+1 = A x k + B u k (9)
where

x k = q k qk and u k = Γ k .
x k+1 is the future state whereas x k and u k are the current state and control input. A and B are the current state and action matrices.

In this study, we learn a model of the dynamics using the state space form of equation ( 9).

Optimal Control

Optimal control is a family of control methods that optimize a cost function over the achievement of a task. In the Linear Quadratic Regulator (LQR) framework in discrete time, the cost function (or optimality criterion) is written:

J = ∞ k=0 (x k T Qx k + u k T Ru k ),
where Q and R are semi-definite positive and definite positive matrices associated to the cost on the state error and on the control input respectively. The minimization of the cost function is obtained by solving the Riccati equation whose unknown is P :

A T P A -A T P B R + B T P B -1 B T P A + Q = P,
where A and B are the state representation matrices defined above.

Solving this equation provides L = R + B T P B -1 B T P A , a constant state feedback matrix used to generate the following feedback controller:

u = L (x des -x) ,
where x = q q and where x des = q des qdes . qdes is computed using the inverse velocity kinematics model described by equation ( 2) whereas q des is extrapolated from the current configuration:

q des = q + qdes Δt.
A schematic view of LQR in shown in Fig. 2.

Fig. 2. LQR Feedback control:

The feedback matrix L is computed from A' and B' using a discrete time Riccati equation

Model Learning

We want to learn models of our system incrementally and from experience. The nonlinear function we want to approximate can be written y = f (x). We use two techniques, one for obtaining a unique linear model β with y ≈ βx and the second for obtaining a nonlinear one f with y ≈ f (x) where f is a combination of linear models.

Least squares and recursive least squares. In the linear case, given some input-output data

X k = [x 1 x 2 ... x k ] and Y k = [
Y k = β k X k by X T k (X k X T k ) -1 on both sides, one gets Y k X + k = β k where X + k = X T k (X k X T k ) -1
and β k represents the average for each term over β 1 ... β k in the least squares sense.

This algorithm can be made incremental through a recursive formulation:

β k = β k-1 + (y k -β k-1 x k ) λ
Thus, from this so-called Recursive Least Squares (rls) method, we can get β k for any k given the previous matrix and a forgetting factor λ that must be tuned.

LWPR.

For most systems, linear models will not be accurate enough and must be replaced by nonlinear or piecewise-linear approximators such as neural networks, decision trees [START_REF] Potts | Incremental learning of linear model trees[END_REF] or radial basis function networks [START_REF] Sun | A fast and efficient model of learning to reach[END_REF].

The lwpr algorithm [START_REF] Vijayakumar | Local dimensionality reduction for locally weighted learning[END_REF] is an incremental radial basis function approximator which provides accurate approximation in very large spaces in a O(k) complexity, where k is the number of sample data. We use it here to learn the FKM of our robot. The algorithm uses a combination of linear models, which are valid on an elliptic zone of the input spaces and whose relative strengths are weighted by a gaussian. This space may evolve during training to match the training data. The domain of validity of each model is called a Receptive Field (RF). The prediction of an entire lwpr model on an input vector is the weighed sum of the results of all the active surrounding RFs. The RFs of a model are created when new input data are not part of any existing RF. Conversely, when a field overlaps another, some criterion can be used to delete it.

Fig. 3. Example of lwpr learning on a sample function

Each RF first projects the input vector on the most relevant dimensions for estimating the output vector by using Partial Least Squares [START_REF] Tenenhaus | La régression PLS: théorie et pratique[END_REF]. During each update, the projector is updated and the algorithm checks whether increasing the complexity by adding another dimension to the input projection significantly reduces the estimation error. If it is so, it modifies the projector accordingly. The projected vector is then used in the n dimension linear model (n being the output dimension) to give the output of the RF.

During prediction with an input vector, the distance between the vector and the RF area is tested for activation and only the significant RFs are activated (see Fig. 3

for an example of RF repartition for function approximation).

The latest version of the algorithm also computes the first and second differential of the learned function with respect to each input dimension (Jacobian and Hessian matrices) when learning the model. lwpr can estimate them for any input value. This calculation is made easier by the fact that the model is a simple sum of multiple linear functions which are easily differentiated. We use this method to extract the Jacobian matrix from the learned model.

Our Motor Learning Model

Our global control scheme, illustrated in Fig. 4, consists of:

a definition of the task in the operational space, resulting on the easy definition of tasks; -learning the inverse velocity kinematics model with lwpr and a high level control law based on a proportional controller given by equation ( 4); -learning the dynamics model with rls and a low level optimal control law optimising a quadratic cost.

Our presentation below distinguishes three layers: high level control, low level control and the coupling between the two.

High Level Control

The operational space control part of our controller consists in the proportional controller given by equation ( 4) based on the difference between ξ des and ξ. ξ is the output of the FKM learned with lwpr, taking the joints positions as input and the operational space position as output.

After equation ( 4), we apply equation ( 5) to transform the error in the operational space into a desired velocity in the joints space. The inverse velocity kinematics model is learned with lwpr taking ξ as input and q as output, as shown in Fig. 5. This approach is similar to the one used in the work of D 'Souza et al. in [32] who also choose to directly learn a specific inverse of the Jacobian matrix using the lwpr algorithm. Given the fact that there is an an infinity of inverses for the Jacobian matrix, the choice of a specific inverse leads to a loss of information regarding the redundant nature of the system. However, the learned Fig. 5. Learning a generalised inverse with lwpr inverse which we can expect to be a weighted pseudoinverse J (q) , i.e. to have minimum norm properties, is the one that best fits the trajectories generated so far.

Low Level Control

At the low level, we use an LQR controller to track the desired joints velocities q computed by the high level control. In the linearised dynamics model x k+1 = Ax k + Bu k , the state x k is [q q] T , the action u k corresponds to the vector of joints torques, and the matrices A and B are learned incrementally with the rls algorithm as shown in Fig. 6, taking (x k ,u k ) as input and x k+1 as output. Note that, in order to bootstrap the learning process of rls, rather than initialising A and B randomly, we rather store samples (x k ,u k ,x k+1 ) and apply a standard LS method until A and B are full rank.

Coupling Both Levels with Anticipation

If we take the current q and q coming from the high level controller as reference, the tracking delay combined with the estimation errors results in an oscillatory behaviour. Thus, we couple the high level and low level control loops by anticipating the values of q and q after some horizon H in the future and using these vectors in the low level control objective function to track them. Algorithm 1. Anticipatory coupling algorithm.

ξ 0 = ξ ; 1 for i = 0 to H do 2 ξi+1 = Kp (ξ des -ξ i ) ; 3 qi+1 = J ξi+1 ; 4 q i+1 = q i + qi+1 Δt ; 5 ξi+1 = F KM q i+1 ; 6 u0H = LQR A, B, [q H qH ] T ; 7
To perform such an anticipation, taking the current operational point ξ as initial value ξ 0 , we perform a loop on ξ i . First, we infer its derivative ξi with a proportional controller (Algol. 1, line 3). This operational velocity can be translated into joints velocities through the weighted pseudoinverse of the Jacobian matrix (line 4). From these joints velocities, we estimate the next configuration (line 5) which, from the FKM, can be translated into an estimation ξ i+1 of the next operational position (line 6). We iterate this anticipation H times and this results in values for q and q after H time steps. Finally, u 0H is the control input applied at moment 0 with horizon H.

Note that the way we estimate the next configuration does not take dynamics effects into account. We assume that our controller will reach the desired state at each step, which may not be true in practice if dynamics effects are not fully balanced by the low level controller. By using the state representation x i = q i qi , one can improve the estimation by using the x i+1 = A x i + B u i into the anticipation loop, but preliminary experiments have shown that this does not result in a significant difference in our context.

Empirical Study

A Simplified Virtual Humanoid System

To evaluate our control architecture, we work with a simplified 3 degrees of freedom (DOFs) mannequin model whose feet are fixed to the ground so as to alleviate the need to care about equilibrium. Thechoice of this simple model is intented to facilitate the presentation of our approach which we expect to be easily scalable to much higher dimension systems. The mannequin is simulated with Arboris [START_REF] Barthlemy | Stability measure of postural dynamic equilibrium based on residual radius[END_REF], a Matlab simulator dedicated to the study of the dynamics of poly-articulated tree-like systems. Its parameters of our simplified mannequin model are extracted from anthropometric tables [START_REF] Wieber | The HuMAnS toolbox, a homogeneous framework for motion capture, analysis and simulation[END_REF]. It is 1.74 meters tall and weighs 73 kg.

Task and Parameters

The task consists in making the mannequin stand up from a crouching posture, so that the vertex (upper extremity of the head) reaches 1.70 meters and stabilises over 1.65 meters, as shown in Fig. 7. Thus our operational task consists in reaching ξ des = ξ x ξ y T = 1.70 0.05 T .

As for parameters, in (4), after testing diverse values, we take K P = 10. Furthermore, a preliminary study of the low level control loop resulted in tuning the Q and R matrices of the LQR controller as follows:

Q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
0.4 0 0 0 0 0 0 0.6 0 0 0 0 0 0 0.3 0 0 0 0 0 0 0.05 0 0 0 0 0 0 0.2 0 0 0 0 0 0 0.3

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ R = 10 -4 ⎡ ⎢ ⎢ ⎣ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ⎤ ⎥ ⎥ ⎦

Empirical Results

Fig. 8 shows that, for a controller without an anticipation loop, the system is able to stand-up, but a vertical oscillation may be observed. With a short anticipation (h = 2 to 6), the system still reaches the operational target but is able to stay there for a longer time period and without oscillation. For a larger anticipation, the system is not able to maintain the desired position any more. Fig. 9 shows an example of the input/output data used by lwpr to learn. After building the corresponding model, lwpr can predict the output q from the input ξ. Fig. 10 shows that, if we perform only one learning epoch along the trajectory, the predicted output cannot be used. We must perform 20 learning epochs where each set of input/output data is presented ten times so as to get a satisfactory model of the weighted pseudoinverse of the Jacobian matrix. Note that, once the system is stabilised, we must stop learning with lwpr. Indeed, if the system is stable, its output is constant and the learned models degenerate. The persistent activation principle used in adaptive control [START_REF] Sastry | Adaptive control: Stability, convergence, and robustness[END_REF] can be used here to stop the learning process.

Discussion and Perspectives

Learning the Dynamics Model

Learning a linear model of the dynamics with rls is a limited approach that can tackle easy problems such as the control of our 3 DOFs system, but as soon as the dynamics gets more complicated, the approximation error will get too large. Our most immediate future work will consist in using lwpr to learn a piecewise linear model of the dynamics. As explained in section 2.5, one can learn a nonlinear evolution of the A and B matrices with the state [q k qk ] T and the control input [u k ] as input at any instant and the next state q k+1 , qk+1 T as output.

However, it is not trivial to exploit the explicit dynamics model learned by lwpr under the form of matrices as it was done with rls. The output given by lwpr is a weighted sum of the outputs of linear models multiplied by a gaussian. The global dynamics model matrices cannot be reconstructed easily from this collection of linear models. Nevertheless, we are working on directly learning an inverse dynamics model in order to control our system with a computed torque [START_REF] Siciliano | Springer Handbook of Robotics[END_REF] approach. Furthermore, we need to control the internal mobility (see equation 3) of the system when it is redundant, which boils down to learning a specific weighted pseudoinverse. This is a matter for future work.

Related Models

An important related model is the work of Mitrovic [START_REF] Mitrovic | Adaptive optimal control for redundantly actuated arms[END_REF], who combines learning a dynamics model of a 6 DOFs arm with lwpr and performing Iterative Linear Quadratic Gaussian (iLQG) optimal control [START_REF] Todorov | A generalized iterative LQG method for locally-optimal feedback control of constrained nonlinear stochastic systems[END_REF] based on that learned model. To our knowledge, Mitrovic et al. are the only authors who succeeded in using a dynamics model learned with lwpr so far. By contrast, N'Guyen-Tuong [START_REF] Nguyen-Tuong | Learning inverse dynamics: a comparison[END_REF] performed such a learning process on more general systems in the context of a comparison with other learning approaches, but without using the model to control a system. Note that iLQG does not use the operational versus joints space distinction, thus it does not benefit from the dimensionality reduction properties of OSC and is limited to addressing systems with no more than 10 DOFs.

Another important related model is MMRL. As in MMRL, the low level controller of our model is based on LQR. So far, our model does not benefit from the modular specialisation mechanism that MMRL shares with MOSAIC. However, such a modularity property will come into play as soon as we will use lwpr at the dynamics level, since lwpr calls upon a collection of local linear models. But here again, the main originality of our approach results from the way we use OSC to specify the task in the operational space. To our knowledge, the only other model that is endowed with the same property is presented in [START_REF] Butz | Exploiting redundancy for flexible behavior: Unsupervised learning in a modular sensorimotor control architecture[END_REF] and is based on neural networks.

State Estimation

In our model, we considered the state of the system exactly known at any moment. In robots, if the FKM of the system is known and if fast enough sensors in servo-motors can give an accurate estimation of the configuration, one can get a good estimation of the current operational state of the system. But when the FKM is learned and inaccurate, state estimation must be used. To perform such an estimation, one must combine informations coming from several exteroceptive and proprioceptive sensors whose accuracy is varying under the environmental conditions and whose delay is generally not compatible with the constraints of feedback control-based motion. Even motor control itself is noisy, as we highlighted in the introduction. Thus maintaining an accurate estimation of the joints state of the musculoskeletal system is a central issue that we must address in the future.

Conclusion

We have presented a model of human motor learning that combines two levels. At the higher level, we specify a task in the operational space and use the OSC approach to derive a target trajectory that will be tracked in the joints space by the lower level. The transformation from the operational space to the joints space is performed with an inverse Jacobian matrix that learned with the lwpr algorithm. At the lower level, we use LQR control to compute the optimal joints torques that realizes the trajectory specified in the operational space, based on a simple dynamics model learned with rls. In our approach, optimal control does not lead to a globally optimal motion but rather to a piecewise optimal solution with respect to the learned system dynamics. In that respect, we can expect our method to be more easily scalable to large dimension systems and it would be of interest to study the best compromise between the computational complexity and the chosen optimal control horizon.

We illustrated this model on a stand-up task. One advantage of the combination of OSC with learning is its flexibility in the definition of the task. Indeed, if we change the end-point effector with respect to which the goal is expressed, or if we change the actuators with which the low-level control is performed, we should just need to learn again the corresponding dynamics to obtain a different controller. Validating this property is in our future work agenda.

Furthermore, we have shown that, due to the approximation error of the linear model learned with rls, we need to stabilise the control architecture by tuning the anticipation horizon of the low level control loop with respect to the high level one. Thus anticipation appears in our control architecture both in model learning processes and in the coupling between the high and the low level control loops.

Fig. 4 .

 4 Fig. 4. Global control scheme

Fig. 6 .

 6 Fig. 6. Learning the dynamics model with rls

Fig. 7 .

 7 Fig. 7. Our 3 DOFs mannequin in its initial and final configurations

Fig. 8 .Fig. 9 .

 89 Fig. 8. Trajectory as a function of the anticipation horizon H

Fig. 10 .

 10 Fig. 10. Error in prediction of joints velocities after one (a) and twenty (b) learning epochs where each input/output data was presented ten times