Bilal Piot
email: bilal.piot@inria.fr

Matthieu Geist
email: matthieu.geist@centralesupelec.fr

Olivier Pietquin
email: olivier.pietquin@univ-lille1.fr

Batch Policy Iteration Algorithms for Continuous Domains Batch Policy Iteration Algorithms for Continuous Domains

Keywords: Batch Reinforcement Learning, Policy gradient, 1, The notion

This paper establishes the link between an adaptation of the policy iteration method for Markov decision processes with continuous state and action spaces and the policy gradient method when the differentiation of the mean value is directly done over the policy without parameterization. This approach allows deriving sound and practical batch Reinforcement Learning algorithms for continuous state and action spaces.

Introduction

Policy Search (PS) (Ng and Jordan, 2000;[START_REF] Fix | Monte-carlo swarm policy search[END_REF] and more specifically Policy Gradient (PG) methods [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF][START_REF] Peters | Natural actor-critic[END_REF][START_REF] Degris | Model-free reinforcement learning with continuous action in practice[END_REF][START_REF] Silver | Deterministic policy gradient algorithms[END_REF] are well known to have a practical edge compared to value-based methods when it comes to Reinforcement Learning (RL) for Markov Decisions Processes (MDPs) with continuous state and action spaces. This principally comes from the fact that value-based methods, either inspired by the Policy Iteration (PI) [START_REF] Lagoudakis | Least-squares policy iteration[END_REF] or Value Iteration (VI) algorithms [START_REF] Ernst | Tree-based batch mode reinforcement learning[END_REF][START_REF] Riedmiller | Neural fitted q iteration-first experiences with a data efficient neural reinforcement learning method[END_REF], rely on the computation of the maximum of the action-value function over the set of actions: max a∈A Q(π, s, a) (greedy step). If this set is infinite, then those methods become intractable. However, it is possible to adapt Policy Iteration-based RL algorithms, such as Least Squares Policy Iteration (LSPI) [START_REF] Lagoudakis | Least-squares policy iteration[END_REF], by changing this global greedy step into a local improvement of the policy made by a gradient step. Indeed, let us take a closer look to the PI method for a finite action space. It consists in two steps, a first step of evaluation of the current policy π k in order to obtain the action-value function Q(π k , s, a) followed by a greedy step π k+1 (s) = argmax a∈A Q(π k , s, a) which improves the action-value function:

Discrete PI scheme

• Evaluation of Q(π k , ., .).

• Greedy step: π k+1 (s) = argmax a∈A Q(π k , s, a).

In practice, when faced with batch data, the evaluation step can be realized for instance by the Least Square Temporal Difference (LSTD) algorithm [START_REF] Bradtke | Linear least-squares algorithms for temporal difference learning[END_REF] if features for the action-value function are provided. When no features are provided, one could easily adapt for instance the Fitted-Q algorithm [START_REF] Ernst | Tree-based batch mode reinforcement learning[END_REF] that evaluates the optimal action value function in order to evaluate the current value function by replacing the optimality-operator by the evaluation operator.

Obviously, when the action space is continuous such a greedy step becomes cumbersome. However, a local improvement is still possible if we are able to compute the gradient of the actions-value function over the possible actions ∂Q(π k ,s,π k (s)) ∂a . 1 Then, the global greedy step can be replaced by a local gradient-like improvement step:

π k+1 (s) = π k (s)+α k ∂Q(π k ,s,π k (s)) ∂a
, where α k ∈ R + . Thus, a canonical adaptation of the PI method in the continuous context is:

Continuous PI scheme

• Evaluation of Q(π k , ., .). • Local improvement step: π k+1 (s) = π k (s) + α k ∂Q(π k ,s,π k (s)) ∂a .
Our main contribution (see Sec. 3) consists in showing, after introducing some notations relative to MDPs and differentiation in Sec. 2, that this PI method for continuous states and actions MDPs (called continuous PI) is in fact a sound algorithm. More precisely, it is a direct derivation of the PG method when the differentiation of the mean value J ν (π) = S V π (s)ν(ds) is directly done over the policy without parameterizing it. This establishes a strong link between PI and PG methods in the continuous setting. It should be noted that this link between PI and PG methods have already been studied in the finite action setting [START_REF] Scherrer | Local policy search in a convex space and conservative policy iteration as boosted policy search[END_REF] where the relation between Conservative Policy Iteration (CPI) [START_REF] Kakade | Approximately optimal approximate reinforcement learning[END_REF] and PG methods has been highlighted. Finally, in Sec. 4 as a second contribution, we provide several practical instantiations of this continuous PI method which can be either feature-based or no.

Background and notations

Before introducing notations specific to MDPs and differentiation of functions, we give some general notations. Let B be a Borel space, ν a measure on B and f ∈ R B a real function, the integral of f under ν, if it exists, is noted B f (b)ν(db). The set of probability measures of the Borel space B is noted ∆ B and we have E

ν [f (X)] = B f (b)ν(db), where E ν [f (X)]
is the expectation of the function f under the probability measure ν.

Let (H, ., . H ) (where ., . H is the dot product associated to H) and (G, ., . G ) be Hilbert spaces,

f ∈ G → R H and ν ∈ ∆ G , then the notation E ν [f (X)] = G f (g)ν(dg) means that G f (g)ν(dg) ∈ R H is a function such that ∀h ∈ H, ( G f (g)ν(dg)) (h) = G f (g)[h]ν(dg). In particular if f ∈ G H , then G f (g), . H ν(dg) ∈ R H is a function such that ∀h ∈ H, ( G f (g), . H ν(dg)) (h) = G f (g), h H ν(dg).

MDP

We give some brief definitions relative to MDPs with continuous state and action spaces. More precisely, we focus on MDPs where the state space S = R p and the action space A = R q with (p, q) ∈ N * 2 . In that particular case, S and A can be seen as both Hilbert spaces (where the dot product is chosen to be the canonical dot product) and Borel spaces with the canonical distance. A continuous MDP (a complete and more formal presentation is done by Hernández-Lerma and Lasserre (1996)) is a tuple (S, A, R, P, γ) where the state space S = R p represents the states of the environment, the action space A = R q represents the possible actions the agent can take, the reward R ∈ R S×A is a bounded function that represents the local benefit of doing action a in state s, the dynamics P (.|.) is a stochastic kernel on S given S × A and γ ∈]0, 1[ is a discount factor. More precisely, if B is a measurable set of S and (s, a) ∈ S × A, then P (B|s, a) is the probability that the agent is in the measurable set B after doing action a in state s. Here, we focus on deterministic policies (π ∈ A S ) which are mappings from states to actions. More precisely, we are interested in the set of all equivalence classes of policies π ∈ A S such that the Lebesgue integral S π(s) 2

A ds, where . A is the canonical norm in A, is finite. This space of policies will be noted L 2 (A, S) or in short L 2 and it is an Hilbert space with the following natural dot product

∀f, g ∈ L 2 , f, g L 2 = S f (s), g(s) A ds,
where ., . A is the canonical dot product relative to A. To evaluate a policy, we use the concepts of value function V ∈ R L 2 ×S and action-value function

Q ∈ R L 2 ×S×A : Q(π, s, a) = E π s,a [ +∞ t=0 γ t R(s t , a t )], V (π, s) = Q(π, s, π(s))
, where E π s,a is the expectation over the distribution of the admissible trajectories (s 0 , a 0 , s 1 , ., . . . ) obtained by executing the policy π starting from s 0 = s and a 0 = a. These functions are well-defined when the policies are deterministic and the reward bounded [START_REF] Hernández-Lerma | Discrete-time Markov control processes[END_REF]. Moreover, V (π, s) follows the Bellman equation:

V (π, s) = R(s, π(s)) + γ S V (π, s )P (ds , s, π(s)).
(1)

One can show that verifying the Bellman equation is equivalent to the existence of a distribution d ν,π ∈ ∆ S , called the γ-weighted occupancy distribution induced by policy π starting with the distribution ν ∈ ∆ S , such that:

∀ν ∈ ∆ S , S V (π, s)ν(ds) = 1 1 -γ E dν,π [R(X, π(X))] = 1 1 -γ S R(s, π(s))d ν,π (ds).
This result is important and will be used in one of our proofs. Explicit formulae (depending on ν, P, γ, π) of d ν,π exist [START_REF] Peters | Natural actor-critic[END_REF][START_REF] Scherrer | Local policy search in a convex space and conservative policy iteration as boosted policy search[END_REF] but require more definitions and here the existence is sufficient to prove our results.

Differentiability

In this section, we give some brief definitions relative to differentiable functions. More properties needed for some of our proofs are given in the appendix A. Let (X, . X ) and (Y, . Y ) be two normed vector spaces and f ∈ Y X a function. We say that f is Fréchetdifferentiable at the point x ∈ X if there exists a linear and continuous function Df (x) ∈ Y X and a function ∈ Y X such that:

∀h ∈ X, f (x + h) = f (x) + Df (x)[h] + h X (h),
where lim

h X →0 (h) = 0. A function f is differentiable if it is differentiable for all x ∈ X.
If f is linear and continuous then ∀x ∈ X, Df (x) = f . Let (H, ., . H ) be an Hilbert space and f ∈ R H a real and differentiable function. Let x ∈ H, Df (x) ∈ R H is a linear form and by the Riez theorem, there exists a vector noted

∂f (x) ∂h ∈ H such that Df (x) = ∂f (x)
∂h , .

H

. The function ∂f (x) ∂h , .

H ∈ R H is such that ∀y ∈ H, ∂f (x) ∂h , . H (y) = ∂f (x)
∂h , y

H

. The vector ∂f (x) ∂h is called the gradient of f at x (where h is a dummy variable canonically associated to H).

Formal link between policy-gradient and continuous LSPI

In this section, we show that continuous PI is in fact an uphill method that maximizes the mean value. Therefore, it is closely related to the PG method. For sake of simplicity, we assume that all operations of integration and differentiation are licit.

The PG method, introduced by [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF], consists in maximizing the mean value J ν ∈ R L 2 , where ν ∈ ∆ S and ∀π ∈ L 2 , J ν (π) = S V (π, s)ν(ds), by a gradient ascent. Thus, it consists in the following algorithm:

π k+1 = π k + α k ∂J ν (π k ) ∂π ,
where α k ∈ R * + is the gradient step. For now, we do not consider the problem of the evaluation of ∂Jν (π k ) ∂π which is key in all PG methods. A more general uphill algorithm which allows to maximize J ν [START_REF] Nocedal | Numerical optimization[END_REF] is:

π k+1 = π k + α k p k , where p k is such that ∂Jν (π k )
∂π , p k > 0. In practice, PG methods do not consider a general set of policies such as L 2 but a set of parameterized policies [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF][START_REF] Peters | Policy gradient methods for robotics[END_REF][START_REF] Silver | Deterministic policy gradient algorithms[END_REF]. Therefore, the differentiation is not done over a policy but over a vector of parameters. Here, we are going to keep this more general kind of differentiation as it will highlight the link with continuous PI. With the more formal notations described in Secs. 2.2 and A.3, the update of continuous PI can be written:

π k+1 = π k + α k ∂Q(π k , ., π k (.)) ∂a ,
where ∀π ∈ L 2 , ∂Q(π,.,π(.)) ∂a ∈ L 2 is a function such that:

∀s ∈ S, ∂Q(π, ., π(.)) ∂a (s) = ∂Q(π, s, π(s)) ∂a .

Therefore, if we manage to show that ∂Jν (π) ∂π , ∂Q(π,.,π(.)) ∂a L 2 > 0, it will imply that continuous PI is an uphill algorithm for the mean value J ν . This will establish the link between maximizing the mean value by gradient ascent (the PG method) and the continuous PI approach.

To do so, we suppose that the mean value J ν (π) = S V (π, s)ν(ds) has the following property:

∀(ν, π) ∈ ∆ S × L 2 , DJ ν (π) = ∂ S V (π, s)ν(ds) ∂π , . L 2 = S ∂V (π, s) ∂π , . L 2 ν(ds).
This property only implies that the permutation between the integral and differential operations is licit and it can be seen as a generalization of the Liebniz integral rule [START_REF] Flanders | Differentiation under the integral sign[END_REF]. Then, thanks to Th. 1, we find an expression of ∂Jν (π) ∂π and thanks to Th. 2 we calculate ∂Jν (π) ∂π , ∂Q(π,.,π(.)) ∂a L 2 . It is important to note that Th. 1 is just another variant of the classical calculus of the gradient of the mean value with respect to the policy originally done by [START_REF] Sutton | Policy gradient methods for reinforcement learning with function approximation[END_REF] for parameterized stationary policies in the discrete scenario and by [START_REF] Silver | Deterministic policy gradient algorithms[END_REF] for parameterized deterministic policies in the continuous scenario. However, here this calculus is done over a non-parameterized deterministic policy and a proof is given for sake of completeness.

Theorem 1 DJ ν (π) = ∂Jν (π) ∂π , . L 2 = 1 1-γ S ∂Q(π,s,π(s)) ∂a , . A • δ s d ν,π (ds)
,

where ∀s ∈ S, δ s ∈ A L 2 is a function such that ∀π ∈ L 2 , δ s (π) = π(s).
Proof The proof is provided in the appendix B.

Theorem 2 ∂Jν (π) ∂π , ∂Q(π,.,π(.)) ∂a L 2 = DJ ν (π) ∂Q(π,.,π(.)) ∂a ≥ 0,

where the equality happens when ∂Q(π,.,π(.)) ∂a = 0 almost everywhere (relatively to the measure d ν,π ).

Proof

The proof is provided in the appendix C.

We conclude that the continuous PI method is indeed an uphill method to maximize the mean value which makes it a sound algorithm. However, to make it a practical batch RL algorithm, we still need to show how we can efficiently compute ∂Q(π k ,.,π k (.)) ∂a from a set of batch data

D RL = (s i , a i , s i , r i ) N RL i=1
where s i ∼ P (.|s i , a i ) which is done in the next section.

Toward practical algorithms

To provide a practical algorithm for the proposed continuous PI scheme, one has to estimate state-action value functions and to differentiate them respectively to the actions. A natural approach consists in considering a linearly parameterized state-action value function Q θ (s, a) = θ φ(s, a), for some feature vector φ : R p+q → R d . For a given policy π k , the parameter vector θ k can be estimated using the LSTD algorithm [START_REF] Bradtke | Linear least-squares algorithms for temporal difference learning[END_REF]. Notice that considering continuous actions in LSTD causes no problem. Then, the policy is updated with ∂Q θ (s,a) ∂a = ∂φ(s,a) ∂a θ. If the components of the feature vector are radial basis functions, for example, this gradient can be easily analytically computed. The computed policy after K iterations is then π K (s) = ∂φ (s,a) ∂a K k=1 α k θ k , each parameter vector θ k being computed thanks to LSTD. This approach requires a linear parameterization. A nonlinearly parameterized stateaction value function could still be estimated by using an iterated projected fixed-point approach. For example, if the Q-function is represented as a neural network, the neural fitted-Q approach of [START_REF] Riedmiller | Neural fitted q iteration-first experiences with a data efficient neural reinforcement learning method[END_REF] could be easily adapted (this would amount to replace the Bellman optimality operator by the Bellman evaluation operator). Moreover, with such a neural network representation, the gradient respectively to the action can be efficiently computed using backpropagation.

In some case, it might be more natural to encode directly a policy, for example π Θ (s) = Θ φ(s), with Θ ∈ R d×q a parameter matrix. The proposed approach is critic-based, but such a parameterization can be envisioned. Consider the advantage function A(π, s, a) = Q(π, s, a) -V (π, s), it can be parameterized as A Θ (s, a) = a -π(s), π Θ (s) A . From this, we can parameterize a state-action value function by adding a linearly parameterized value function: Q Θ,w (s, a) = A Θ (s, a) + w ϕ(s). This Q-function can be estimated with LSTD and its gradient is simply ∂Q Θ,w (s,a) ∂a = π Θ (s). The policy computed after K iterations is

then π K (s) = K k=1 α k Θ k φ(s)
, a linear mixture of the policies π Θ k . This is reminiscent of compatible function approximation for (natural) policy gradient with a critic and of the related parameterization of the advantage function, see for example [START_REF] Peters | Natural actor-critic[END_REF].

Alternatively, this last approach can be motivated by a first-order Taylor expansion of the state-action value function when the only varying parameter is the action a:

Q(π, s, a) = V (π, s) + ∂Q(π, s, π(s)) ∂a , a -π(s) A + a -π(s) A (a -π(s)),
where lim a-π(s) A →0 (a-π(s)) = 0. Therefore, the policy π Θ can alternatively be seen as a parameterization of the partial derivative of the state-action value function ∂Q(π,s,π(s)) ∂a . We have considered a linear parameterization for π Θ , but as before, a nonlinear parameterization could be envisioned, by using an iterated projected fixed point approach.

Conclusion

We have shown that the natural adaptation of the Policy Iteration (PI) to the continuous scenario scheme, which consists in replacing the classical global greedy step by a local gradient improvement step, can be motivated by the close link it has to the policy gradient method. Indeed, continuous PI can be seen as an uphill method for the mean value as shown in Sec. 3. In addition, we sketch practical related algorithms, that rely mainly on LSTD or an iterated projected fixed point approach if no features are provides, which should make continuous-PI algorithms as easy to use as their discrete counterparts. Finally, as perspective, we would like to test those different algorithms on benchmark problems.

•δ S + ∂Q(π, s, π(s)) ∂a , . A •δ A .
Thus, we have:

DV s (π) = ∂Q(π, s, π(s)) ∂π , . L 2 + ∂Q(π, s, π(s)) ∂a , . A • δ s .
We recall that: Clearly, the equality case happens when ∂Q(π,.,π(.)) ∂a is almost everywhere null (relatively to d ν,π ).

DQ (s,π(s)) (π) = ∂Q(π, s, π(s)) ∂π , . L 2 , where ∀(s, a) ∈ S × A, Q s,a ∈ R L 2 is a function such that ∀π ∈ L 2 , Q (s,a) (π) = Q(π,

P

  s, a). Moreover:Q (s,π(s)) = R • (I s , I a ) + γJ P (.|s,π(s)) ,where ∀a ∈ A, I a ∈ A L 2 is such that ∀π ∈ L 2 , I a (π) = a. Using Eq. (2), Eq. (3) and the fact that I a and I s are constant, we have:DQ (s,π(s)) (π) = DR(s, π(s)) • (0(ds |s, π(s)), = ∂Q(π, s, π(s)) ∂a , . A • δ s + E P (.|s,π(s)) DV X (π) .

So,

  DV s (π) verifies a Bellman equation (1) which implies that:∀ν ∈ ∆ S , E ν [DV X (π

Appendix A. Differential Calculus

In this section, we recall some important properties of differential calculus that are used in our proofs.

A.1 Differential of a Composition

Let (X, . X ), (Y, . Y ), (Z, . Z ) be three normed vector spaces,

(2)

A.2 Differential of a multidimensional output

Let ((Y i , . Y i )) n i be a family of n ∈ N * normed vector spaces, (X, . X ) a normed vector space and (f i ) n i=1 a family of functions such that

, is differentiable at x and:

(3)

A.3 Differential of a multidimensional input form

Let ((H i , ., .

then (H, ., . H ) is an Hilbert space. Now, let consider f ∈ R H differentiable. For i ∈ (1, . . . , n), we use the notation x = (x 1 , . . . , x n ) = (x -i , x i ) where x -i = (x 1 , . . . , x i-1 , x i+1 , . . . , x n ). For a fixed x -i , we can define the function f x -i ∈ R H i such that:

H i .

To simplify the notations, we write ∂f (x) ∂h i in lieu of ∂f x -i (x i ) ∂h i

. ∂f (x) ∂h i is the partial derivative (or the partial gradient) of f evaluated in x with respect to the variable h i canonically associated to the Hilbert space H i . In addition, Df (x) is such that:

, .

where

Appendix B. Proof of Th. 1

Proof We recall that:

, where ∀s ∈ S, V s ∈ S L 2 such that ∀π ∈ L 2 , V s (π) = V (π, s). As, for all s ∈ S, we have V (π, s) = Q(π, s, π(s)) then:

where 2) and Eq. ( 3), we have:

As I s is a constant and I L 2 and δ s are linear and continuous, then:

DV s (π) = DQ(π, s, π(s)) • (I L 2 , 0, δ s ).

We also recall that: DQ(π, s, π(s)) = ∂Q(π, s, π(s)) ∂π , .

L 2

•δ L 2 + ∂Q(π, s, π(s)) ∂s , .


Guilherme Maeda
email: maeda@ias.tu-darmstadt.de

Marco Ewerton
email: ewerton@ias.tu-darmstadt.de

Takayuki Osa
email: osa@ms.k.u-tokyo.ac.jp

Baptiste Busch
email: baptiste.busch@inria.fr

Jan Peters
email: peters@ias.tu-darmstadt.de

Active Incremental Learning of Robot Movement Primitives

Keywords: Gaussian process, movement primitives, active learning

Robots that can learn over time by interacting with non-technical users must be capable of acquiring new motor skills, incrementally. The problem then is deciding when to teach the robot a new skill or when to rely on the robot generalizing its actions. This decision can be made by the robot if it is provided with means to quantify the suitability of its own skill given an unseen task. To this end, we present an algorithm that allows a robot to make active requests to incrementally learn movement primitives. A movement primitive is learned on a trajectory output by a Gaussian Process. The latter is used as a library of demonstrations that can be extrapolated with confidence margins. This combination not only allows the robot to generalize using as few as a single demonstration but more importantly, to indicate when such generalization can be executed with confidence or not. In experiments, a real robot arm indicates to the user which demonstrations should be provided to increase its repertoire of reaching skills. Experiments will also show that the robot becomes confident in reaching objects for whose demonstrations were never provided, by incrementally learning from the neighboring demonstrations.

Introduction

Robots will eventually make part of our daily lives, helping us at home, taking care of the elderly, and collaborating at work. Such robots cannot be pre-programmed prior to deployment, and will most likely start with an empty set of skills. One challenge that must be addressed is that of creating, augmenting and tailoring a robot's repertoire of skills as needed, over time. While motion planning offers a general approach to generate robot motions, its requirements are difficult to satisfy under unstructured and uncertain collaborative scenarios: the environment changes constantly, perception is prone to occlusions, and non-technical users are not expected to design cost functions required by trajectory optimizers. Imitation learning can avoid such drawbacks by relying on the presence of a human teacher. Imitation learning takes advantage of elementary, parameterized movements, referred here as movement primitives, to encode and generalize human demonstrations. However, imitation learning has primarily addressed how to endow and refine robots with motor skills. It usually does not address when the learning should take place. Reasoning when improvement is actually needed is, nevertheless, an essential and difficult problem to be solved. This paper focuses on an algorithm for such reasoning during the learning of movement primitives.

The usual practice in imitation learning is to rely on human judgment to decide when a skill must be added or corrected. Human judgment in robot skills, however, presents a number of issues. For example, it requires constant human attention: as the robot moves, the user must observe and decide if a correction or addition of a new skill is needed, which is not always obvious. Also, the human cannot judge if a generalized skill is suitable unless the robot actually executes it, which can lead to 1st Conference on Robot Learning (CoRL 2017), Mountain View, United States.

undesirable/dangerous situations. On the other hand, if the robot has means to quantify the quality of its own skills for a given task, it can then decide how confident it is to execute it. We frame this problem as an active learning problem where the human provides demonstrations rather than labels.

The principal contribution of this paper is an active learning algorithm that allows a robot to reason about the confidence of its movement primitives. The robot can decide when a demonstration is required, making active requests to a human depending on its confidence. This capability also sheds light onto the problem of deciding how many demonstrations are needed to construct a probabilistic model when learning from demonstrations. Under active learning, the number of demonstrations is indicated by the robot, on-demand. A secondary contribution is a methodology to train Dynamical Movement Primitives (DMPs) [START_REF] Ijspeert | Dynamical movement primitives: learning attractor models for motor behaviors[END_REF] with contextualized demonstrations encoded by Gaussian Processes (GPs).

Related Work

Active learning addresses the problem of which labels to ask in order to improve the sampling efficiency of a learning agent. This idea strongly connects with the way a physical robot should improve its skills by asking a human for demonstrations. Not surprisingly, active learning has been explored in different contexts within robotics, particularly when associated with imitation learning [START_REF] Shon | Active imitation learning[END_REF] and reinforcement learning. For example, an agent has been shown to learn navigation policies more efficiently by using both autonomous execution and demonstrations when compared to a pure reinforcement learning agent [START_REF] Chernova | Confidence-based policy learning from demonstration using gaussian mixture models[END_REF][START_REF] Chernova | Interactive policy learning through confidence-based autonomy[END_REF]. Active learning also plays an important role in developmental robotics for the lifelong acquisition of robot skills, either by self-exploration and self-improvement [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF], or by requesting human assistance via demonstration [START_REF] Grollman | Dogged learning for robots[END_REF].

The interplay between active learning and skill learning on physical robots presents many open problems. One of them regards the appropriate metric used to trigger a request for assistance, such as novelty and uncertainty reduction [START_REF] Silver | Active learning from demonstration for robust autonomous navigation[END_REF], or the confidence in executing an action [START_REF] Chernova | Interactive policy learning through confidence-based autonomy[END_REF]. Other problems relate to the frequency of requests [START_REF] Cakmak | Designing interactions for robot active learners[END_REF], the modality of the input [START_REF] Cakmak | Designing robot learners that ask good questions[END_REF], and the learning of the threshold that triggers active requests [START_REF] Chernova | Confidence-based policy learning from demonstration using gaussian mixture models[END_REF]. In lifelong learning, one problem is that of defining in which space (motor or goal) a learning agent learns more efficiently [START_REF] Moulin-Frier | Learning how to reach various goals by autonomous interaction with the environment: unification and comparison of exploration strategies[END_REF]. In this paper, uncertainty reduction will be used as the metric of choice, the threshold will be manually set, and skills will be learned in the robot's goal space. These practical decisions will allow us to focus on the algorithmic problem of programming movement primitives for a robot under active requests.

While incremental learning has been addressed by a number of authors, the active learning aspect has not been explored to the same level. Calinon and Billard [START_REF] Calinon | Active teaching in robot programming by demonstration[END_REF] addressed the problem of incrementally refining robot primitives with Gaussian Mixture Models (GMMs). In a similar fashion, Lee and Ott [START_REF] Lee | Incremental kinesthetic teaching of motion primitives using the motion refinement tube[END_REF] proposed a refinement tube defined by impedance control where trajectories are encoded with Hidden Markov Models. Ewerton et al. [START_REF] Ewerton | Incremental imitation learning of context-dependent motor skills[END_REF] introduce a method to incrementally improve policies by interleaving feedback control with corrective interventions in open-loop. Farraj et al. [START_REF] Farraj | A learningbased shared control architecture for interactive task execution[END_REF] uses a measure of information gain to decide when new data should be added to the robot's skill repertoire. Ahmadzadeh et al. [START_REF] Ahmadzadeh | Generalized cylinders for learning, reproduction, generalization, and refinement of robot skills[END_REF] use a generalized cylinder representation to encode multiple demonstrations and refine the parts that need correction. In common is the fact that in the cited works, the human is in charge of deciding when to teach and refine the robot skills. Here, our goal is to allow a robot to actively take such decisions. To achieve the desired confidence reasoning, the choice of the representation used to encode the training data is crucial. Probabilistic modeling seems, therefore, the natural choice for enabling such reasoning. Probabilistic Movement Primitives (ProMPs) [START_REF] Paraschos | Probabilistic movement primitives[END_REF] treat multiple demonstrations as multi-variate Gaussians, but present the issue that an unknown number of initial demonstrations are required to create an informative prior. The same problem is also found with GMMs [START_REF] Calinon | Active teaching in robot programming by demonstration[END_REF][START_REF] Calinon | Learning-based control strategy for safe humanrobot interaction exploiting task and robot redundancies[END_REF]. Additionally, as models based on mixtures of linear components, such methods rely on heuristics to define the number of components. Moreover, an issue of the linearity is that consistency of the generalization is only guaranteed at the vicinity of the training data. In this respect, GPs have been used to represent robot primitives with great generalization capabilities [START_REF] Schneider | Robot learning by demonstration with local gaussian process regression[END_REF][START_REF] Osa | Online trajectory planning and force control for automation of surgical tasks[END_REF] at the expense of a higher computational cost. Dynamical Movement Primitives (DMPs) have been known to generalize well around the vicinity of a demonstration and its capabilities have been extended in many ways with GPs. Ude et al. [START_REF] Ude | Task-specific generalization of discrete and periodic dynamic movement primitives[END_REF] use a combination of local weighted regression and GPs to predict the parameters of DMPs in a database. In [START_REF] Forte | On-line motion synthesis and adaptation using a trajectory database[END_REF], GPs are used to learn DMPs parameters in joint space while querying states in the input space of the GP as the final desired Cartesian position of the movement. Matsubara et al. [START_REF] Matsubara | Learning parametric dynamic movement primitives from multiple demonstrations[END_REF] use GPs to learn the parameters of DMPs conditioned by the style of the movement while Bitzer and Vijayakumar [START_REF] Bitzer | Latent spaces for dynamic movement primitives[END_REF] use GPs as a dimensionality reduction approach to encode DMPs in a latent space. These prior works show that the combination of GPs and DMPs are useful to contextualize parameters of the movement but do not exploit the use of uncertainty for active learning. In this paper, we focus on the uncertainty provided by the GPs. For this reason, we use GPs and DMPs in two separate steps. The idea is to obtain confidence bounds in the same space in which demonstrations are provided: as trajectories (as opposed to the parameters of DMPs) so that the value of the uncertainty is physically meaningful, allowing for empirical tunning.

A Movement Primitive with Uncertainty for Active Request

GPs have been known to extrapolate with an indication of uncertainty. However, extrapolations inevitably lead to errors, which for many fine motor skills (such as reaching a grasping position) will turn into execution failures. On the other hand, a DMP guarantees that any goal can be reached beyond the training set, but it cannot indicate how far a generalization to an unseen input can be made before a new demonstration is required. When combined in tandem, DMPs will learn trajectories extrapolated by the GPs with confidence bounds. Each component will be explained in the following subsections. To avoid confusion, hereinafter, a query made in the input space of a GP will be referred to as a context. A query that the robot makes to the human to obtain a demonstration will be referred to as a demonstration request.

Generalizing the Training Data with Gaussian Processes

Define the vector x 1:N as a trajectory comprised of N time steps. To provide a concrete, practical explanation, we focus on modeling trajectories in the Cartesian1 space,

x n = [x n , y n , z n ] .
For each time step, independently trained Gaussian Processes (GPs) are used to predict the triple

(f x,n (x * ), f y,n (x * ), f z,n (x * )), where f (•),n (x *
) is a scalar, and x * is a desired location on the input space of the model. We focus on reaching specific positions of the robot workspace, such that x * = [x N , y N , z N ] goal is the context. Since each Cartesian dimension will be treated independently, the subscripts x, y, z will be omitted and f n (x) will be used to represent a generic dimension.

The goal of the GPs is to predict a trajectory as a sequence f 1:N (x * ) for an unseen context x * by learning a map from context to trajectories M : x * → f 1:N (x * ). For active requests of demonstrations, the prediction is in the form of a distribution p(f 1:N |x * ). Assume that the positions of D demonstrated trajectories at each n-th step t n = [t 1 n , ..., t D n ] and the final positions (i.e. con-texts) achieved by their respective trajectories X = [x 1 , ..., x D ] are jointly Gaussian. The joint distribution is given as

t n f n (x * ) ∼ N 0, K(X, X) + σ 2 s I k(X, x * ) k(x * , X) k(x * , x * ) ,
where σ 2 s is a noise variance to be trained as a hyper-parameter. For a single context query, k(x * , x * ) is a scalar, and the covariance matrices K(X, X), k(X, x * ) have dimensions D × D, D × 1, respectively. Each element K i,j = k(x i , x j ) is computed with the squared exponential covariance function

k(x i , x j ) = σ 2 f exp - 1 2 (x i -x j ) W(x i -x j ) , (1)
where σ 2 f and W are hyper-parameters. The proposed GP encoding is illustrated in Figure 1. The GP prediction is

f n (x * ) ∼ N (µ n , σ 2 n ), (2)
where

µ n = k(x * , x)[K(x, x) + σ 2 s I] -1 t n , (3)
and

σ 2 n = k(x * , x * ) -k(x * , x)[K(x, x) + σ 2 s I] -1 k(x, x * ). (4)
The variance σ 2 n is the essential component for active requests. In this paper, the uncertainty is computed as the average magnitude over the whole trajectory

u = 1 N N n=1 (σ x,n + σ y,n + σ z,n ) 2 . (5)
The computation in the form of ( 5) was chosen solely due to its simplicity. More sophisticated ways to compute the uncertainty can potentially improve the performance of the active learning algorithm but such investigation is beyond the scope of this work. The predictions of the GPs do not guarantee that the resulting trajectory will accurately achieve a desired unseen context x * . This is inevitable when the underlying function is highly nonlinear and the GPs have to extrapolate significantly beyond its training set. Figure 2 illustrates two cases of generalization where a context x 1 * is closer to the training set than the context x 2 * . Although none of the extrapolations exactly satisfy the context, only a small correction is required for x 1 * while for x 2 * the accuracy is quite low. The error in prediction is caused by the function nonlinearity and a lack of observations at the vicinity of the desired context, which is also reflected in terms of the resulting uncertainty.

Accurate Executions with Dynamical Movement Primitives

We assume that the predicted GP trajectory is valid when its uncertainty ( 5) is below a specified threshold. This prediction is then encoded in the form of a DMP, guaranteeing that the corrected trajectory can exactly reach the goal. Regarding the input x 2 * , given the large uncertainty of the GPs, it seems prudent to rather request a human demonstration.

Taking y as one of the Cartesian coordinates, a DMP is defined as a dynamical system of the form

τ ÿ = α y (β y (y * -y) -ẏ) + f (z), (6)
where τ is a temporal scaling variable that controls the duration of the movement. The parameters α y , β y are set to emulate a critically damped, linear spring-damper dynamics with stiffness α y β y and damping α y . The context, in the form of the final position of the trajectory, is used as the goal attractor y * .

An arbitrary force f (z) acts on the linear spring-damper system and captures the nonlinearities of the demonstration. The phase variable z controls the evolution of the system and is by itself defined as a dynamical system of the form τ ż = -α z z (see e.g. [START_REF] Ijspeert | Dynamical movement primitives: learning attractor models for motor behaviors[END_REF] for details). go to line 1

Algorithm 2 GP-DMP({x d , t d } 1:D , x * , u trig ) 1: GP ← TRAIN_GP ({x d , t d } 1:D ) 2: (µ 1:N , σ 2 1:N ) ← PREDICT (GP (x * )) 3: if TRAJ_UNC(σ 2 1:N ) < u trig then 4: DMP ← TRAIN_DMP(µ 1:N ) 5: y ← REACH_GOAL(µ 1:N , x * ) 6: else 7: {x D+1 , t D+1 } ← ACTIVE_REQUEST Algorithm 3 InverseKinematics(X [k] 0:T , q i ) 1: q t ← q i 2: for t < T -1 do 3:
while ikErr < minErr do 4:

X t ← ForwarKinematics(q t ) 5: e ← computeFKError(X t , X t ) 6: J t ← Jacobian(q t ) 7:
∆q t ← αJ -1 t e.

8:

q t ← q t + α∆q t 9:

e t ← e
Algorithm 1 shows a pseudocode of the proposed method. As inputs, the algorithm requires a training set comprised of D demonstrations, an unseen context x * , and the trigger of the uncertainty u trig . In line 1, the GPs are initialized. In line 2, the mean and variance of the extrapolated trajectories are computed. In line 4, a DMP is trained on the GP mean and its goal attractor is adjusted (line 5) to accurately reach the desired context. In line 8, the robot requests a new demonstration if the GP variance is beyond the uncertainty trigger.

Experiments on Reaching Tasks with a Real Robot

In the following experiments, a 7 degree-of-freedom lightweight arm was used. The threshold on the uncertainty, computed by [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF], was empirically set to 9 cm. The GPs were trained on the (x,y,z) coordinates of the Cartesian trajectory of the robot hand. The orientation of the robot hand along the trajectory was given by spherical linear interpolation (slerp) between the initial and final orientations of the robot along the trajectory. The final orientation was defined by the reference frame of the object to be reached. Inverse kinematics was used to execute the joint-space tracking controllers. A video of the experiments reported in this section can be watched in https://youtu.be/s9kG_ IKzqO4.

Assembly Scenario

In this first experiment, we considered an industrial environment where a non-technical user has to program a robot collaborator to reach different objects in its workspace. As shown in Figure 3, there were 10 objects to be reached such as different parts of a product on a shelf and tools on a toolbox.

To simplify the realization of the experiment, we hand-coded the reaching locations (position and orientation) of each object. Obstacles and the geometry of the environment are unknown to the robot. The plot on the right shows the desired reference frames of the final positions of the hand, that is, the contexts. The robot indicated for which objects the human should provide demonstrations, with the goal of reaching all objects with confidence.

As the initial training set was empty, the robot randomly selected one object for the first demonstration. After each demonstration, the GPs were trained and predictions for each of the contexts were made. Their respective uncertainties {u 1 , ..., u 10 } were then ranked according to [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF]. The most uncertain reaching object was then used by the robot to request a human demonstration. This process was repeated until all predictions fell below the specified 9 cm threshold. Figure 4(a) shows the result of this procedure for a particular trial where the robot initial random request was to reach the lower plate on the shelf. The prediction for all contexts given this single demonstration is shown in Figure 4(a). The ellipses represent the projections of one standard deviation along the trajectory on the XY plane. In this particular case, the next requested demonstration was for the ratchet, followed by the screw driver. The predictions are shown in Figure 4(b). Finally, demonstrations were requested for the set of hex keys and the upper plate. The predictions are shown in Figure 4(c). The color of the trajectories for uncertainties above the threshold are shown in red, otherwise in blue.

Note from the Figure 4(c) that after five demonstrations all 10 trajectories were within the confidence limits. At this stage, the robot could execute all reaching primitives confidently. This result shows an important distinction of this work when compared to other probabilistic movement representations such as GMMs and ProMPs [START_REF] Calinon | Active teaching in robot programming by demonstration[END_REF][START_REF] Calinon | Learning-based control strategy for safe humanrobot interaction exploiting task and robot redundancies[END_REF][START_REF] Paraschos | Probabilistic movement primitives[END_REF]. Here, the robot dictates which demonstrations are required to generate a library that satisfies the given task while in prior work the user has to make such decision based on his/her own judgment.

A characteristic of the method is that the uncertainty decreases for all objects, and not only for the objects whose demonstrations were requested. This is particularly evident when plotting the decrease of uncertainty for the three screw positions (left, middle, right) as shown in Figure 5.

Number of requests

Left screw Middle screw Right screw During this particular trial, no demonstrations for any of the screws were requested. Nevertheless, the robot could increase its confidence on how to reach the screws because of the demonstrations provided for the neighboring objects. The runs differ by the first object whose trajectory demonstration was requested. We systematically used each of the 10 objects to initialize the first skill in each trial. The results show that, regardless of which object a demonstration was first requested-which led to different sequences of demonstrations-in all cases, a total of five requests showed sufficient to achieve high confidence for all trajectories. In Figure 6(b), we compare the uncertainty decrease criterion with uniform random requests. Note that the uncertainty after the third request decreases at a signifi-

Contextual Movement Primitives

Contextualization of movement primitives has been explored by many authors (e.g. [START_REF] Ude | Task-specific generalization of discrete and periodic dynamic movement primitives[END_REF][START_REF] Matsubara | Learning parametric dynamic movement primitives from multiple demonstrations[END_REF][START_REF] Osa | Online trajectory planning and force control for automation of surgical tasks[END_REF]) as a way to augment the versatility of robot skills. Our method brings the advantage that the uncertainty is used to collect demonstrations for new contexts in an active manner. In the particular case addressed here, the context is the final position of the robot hand. The assumption due to the GPs is that contexts that are close to each other must have trajectories with similar profiles. The goal is not to learn an underlying pattern that governs all demonstrations, but to learn how to generate different motions according to their contexts. Figure 7(a) shows a scenario where the robot has to reach different parts of the workspace in different ways. As the robot hand starts behind the table, to reach the blocks on the table, the hand must move vertically with and arc-like motion. Differently, to reach the blocks on the pole, the hand must contour the back of the table in order to avoid collision. To train the GPs, five demonstrations were provided, one for each red block. These demonstrations are shown in Figure 7(b). Note that the robot does not need to be aware of the table as an obstacle since the human is providing the collision-free demonstration via kinesthetic teaching. query is shown by the red trajectory which gives the robot with possibility to execute or not the solution. The robot could also execute the trajectory slowly while asking the human to supervise its execution, or to request the human for a new demonstration.

Increasing Coverage On-Demand

In a third set of experiments, the robot was trained on a collision-free workspace with the goal of expanding its skills on-demand. Six demonstrations were initially provided by moving the robot hand from its home position to an arbitrary location as shown in Figure 8(a). The user arbitrarily defined new goal positions using a tracked wand. The robot then tried to generate a trajectory to reach this goal given the current set of demonstrations. If the uncertainty (5) was larger than the specified threshold, the robot requested a new demonstration, which was then added to the training set and used to re-train the GPs.

In total, the user indicated 18 new queries whose generalized trajectories are shown in Figure 8(b). For 11 cases, the trajectories could be generalized with confidence. The remainder seven cases required additional demonstrations as the robot uncertainty was above the 9 cm threshold. These additional demonstrations are shown in Figure 8(c). To quantify the performance we consider the volume of the convex hull defined by the trajectories. The initial demonstration set used to initialize the GPs covered a volume of 0.036 m 3 (Figure 8(a)). The generalized trajectories covered a volume of 0.303 m 3 (Figure 8(b)). The additional set of demonstrations covered a volume of 0.197 m 3 (Figure 8(c)). The increased workspace volume achieved by the robot is eight times larger than the initial volume of the initial demonstration set. This experiment indicates that, under the active learning setting, a robot can increase its initial skill set while being efficient in the number of requests.

The previous experiments show that differently from methods based on mixture of models such as GMMs, here, there is no need for computing the correct number of components as the data grow, thus, avoiding heuristics or cross-validations. Also, the user does not need to know, in advance, how many demonstrations he/she need to provide to construct the probabilistic model.

Conclusions

This paper presented an algorithm that allows a robot to incrementally increase its skills by reasoning when to request demonstrations. The method uses GPs to provide a measure on the confidence in which the training set is being extrapolated, which is used to trigger a new demonstration request. DMPs are learned on the GP output to guarantee that the desired position can be achieved exactly.

The combination also provides a way to contextualize the demonstration, thus, offering DMPs the capability to address different types of demonstrations. In this paper, a few decision have been made by the user, such as the criterion for queries, the threshold on the uncertainty, and the space in which to encode demonstrations. We are currently working to automate many of such decisions. We focused on the reaching position as a context, but as a future work, we will consider other contexts such as the position of obstacles, the geometry of objects and the differences in interactions with humans.

Figure 1 :

 1 Figure1: A sketch of the GP approach. Given four trajectory demonstrations (left plot), Gaussian processes are trained for each time slice to predict positions and uncertainties given new queries (contexts). In this paper, queries are defined as the final trajectory position. The plot at the right shows the prediction of a trained GP for the time step 25. In this example, the input space (the trajectory goal or context) is one dimensional.

3 Figure 2 :

 32 Figure 2: The DMPs are learned on the GP output and can correct for inaccuracies. The large uncertainty for x 2* indicates that a demonstration should be requested instead.

Figure 3 :

 3 Figure 3: A robot starts with an empty skill set. It must reach 10 objects whose positions and orientations were given (right). The robot uses the GP-DMP predictions to predict reaching trajectories for all objects. Active requests for demonstrations are made to decrease the uncertainty of the reaching trajectories.

Figure 4 :

 4 Figure 4: Decrease of the trajectory uncertainties as demonstrations for objects are being requested. The ellipses show a single standard deviation. Uncertainties above the threshold are shown in red. The blue trajectories indicate that the robot is confident in executing the primitives. (a) The uncertainties for all contexts after one demonstration. (b) Predictions after two demonstrations. (c) Predictions after five demonstrations.

Figure 5 :

 5 Figure 5: Decrease of the robot's uncertainty in reaching the three screws.

Figure 6 (

 6 Figure 6(a) shows the average result over 10 different runs where the bars represent the averaged uncertainty with a single standard deviation.The runs differ by the first object whose trajectory demonstration was requested. We systematically used each of the 10 objects to initialize the first skill in each trial. The results show that, regardless of which object a demonstration was first requested-which led to different sequences of demonstrations-in all cases, a total of five requests showed sufficient to achieve high confidence for all trajectories. In Figure6(b), we compare the uncertainty decrease criterion with uniform random requests. Note that the uncertainty after the third request decreases at a signifi-

Figure 6 :

 6 Figure 6: Average result over 10 runs where the bars represent the averaged uncertainty of the robot to reach each object. A single standard deviation is shown on the top of each bar. (a) Requests for the most uncertain object. (b) Uniform random requests.

Figure 7 :

 7 Figure 7: Generating different movements according to the location of the context. (a) Demonstrations were provided for the red blocks. Yellow blocks were used for test. (b) Demonstrated trajectories for the red blocks. (c) Generalization towards the yellow blocks. Note that the robot is uncertain for the block away from the training data.

Figure 8 :

 8 Figure 8: Incrementally adding demonstrations on-demand. The numbers indicate the volume of the convex hull formed by the trajectories. (a) Initial training set. (b) Generalizations satisfying the confidence threshold. (c) Additional requested demonstrations.

  Algorithm 1 GP-DMP({x d , t d } 1:D , x * , u trig )

	1: GP.train({x d , t d } 1:D ) 2: (µ 1:N , σ 2 1:N ) ← GP.predict(x * ) 3: if TRAJ_UNC(σ 2 1:N ) < u trig then
	4: 5: 6: 7: else DM P.train(µ 1:N ) y ← DM P.reach_goal(x * ) ROBOT_EXECUTION(y) 8: {x D+1 , t D+1 } ← ACTIVE_REQUEST 9:

Orientations are not addressed in this paper, for simplicity and due to the following two premises. First, the initial and final orientations of the robot hand are given, the former by the robot's initial position, and the latter by the reference frame of the object to be reached. Second, experiments show that interpolation along the trajectory between these two orientations is sufficient for satisfying the particular tasks proposed in this paper.

Acknowledgments

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7-ICT-2013-10) under grant agreement 610878 (3rdHand), from the project BIMROB of the "Forum für interdisziplinäre Forschung" (FiF) of the TU Darmstadt, from the European Union's Horizon 2020 research and innovation programme under grant agreement 640554 (SKILLS4ROBOTS), and 645582 (RoMaNS). The authors would like to thank Riad Akrour for invaluable feedback.


Lucie Daubigney

Matthieu Geist

Olivier Pietquin

RANDOM PROJECTIONS: A REMEDY FOR OVERFITTING ISSUES IN TIME SERIES PREDICTION WITH ECHO STATE NETWORKS

Keywords: Time series, echo state network, random projection

Modelling time series is quite a difficult task. The last recent years, reservoir computing approaches have been proven very efficient for such problems. Indeed, thanks to recurrence in the connections between neurons, this approach is a powerful tool to catch and model time dependencies between samples. Yet, the prediction quality often depends on the trade-off between the number of neurons in the reservoir and the amount of training data. Supposedly, the larger the number of neurons, the richer the reservoir of dynamics. However, the risk of overfitting problem appears. Conversely, the lower the number of neurons is, the lower the risk of overfitting problem is but also the poorer the reservoir of dynamics is. We consider here the combination of an echo state network with a projection method to benefit from the advantages of the reservoir computing approach without needing to pay attention to overfitting problems due to a lack of training data.

INTRODUCTION

Physical processes (such as the evolution of the temperatures in a place or the concentration of a gas in the atmosphere), economical processes, etc., are phenomenons which can be modelled as time series. A reasonable question which may arise once the data are collected could be: would it be possible to predict future samples of a time series given those we already know? This is particulary true when important problems, like climate change, are at stake.

A first solution would be to start from a known model (e.g., a differential equation) and adapt its parameters from the collected samples. However, determining a model is not an easy task when little or even nothing is known about the underlying process. Moreover, this way of solving the prediction problem usually provides only with ad hoc solutions and the crafting task, which can be time consuming, must be done for each new problem.

Soft computing approaches have thus been proposed to find models able to predict time series. The idea is to propose general models for which no hypotheses are required either about the underlying process or about the time dependencies between samples. The reservoir computing approach seems to be particularly appropriate. Indeed, reservoir computing uses Recurrent Neural Networks (RNNs) [START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF][START_REF] Verstraeten | An experimental unification of reservoir computing methods[END_REF] which can catch dependencies over time because of the recurrent connections among its neurons.

Results have been computed with the InterCell cluster funded by the Région Lorraine. The authors want to thank the ANR project METHODEO for partial funding.

An RNN is made of an input layer, a hidden layer (or reservoir) and an output layer. The hidden layer is supposed to contain a large number of neurons. Each neuron is possibly connected to itself and to any of the other neurons on this layer. Its connection to itself induces a memory effect. This layer is thus responsible for catching the dynamics. Supposedly, the larger the number of neurons in the reservoir, the more complex the caught dynamics is. More details can be found in [START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF]. The problem while using RNNs is that classic training methods such as backpropagation do not work properly [START_REF] Bengio | Learning long-term dependencies with gradient descent is difficult[END_REF].

Because of that, Echo State Networks (ESNs) have been introduced [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF]. ESNs are RNNs where both the weights of the connections from the input to the reservoir and the weights of the connections inside the reservoir are set up randomly initially. Some aspects for defining good weights remain unclear but good heuristics have been found to build efficient ESNs [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]. The problem of learning the connections inside the reservoir is thus avoided and only the weights of the output layer are learnt. To do so, several algorithms have been proposed, such as a simple linear regression or most commonly least-squares [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]. Especially, an online version of this algorithm, the Recursive Least Squares (RLS) has been successfully experienced in [START_REF] Jaeger | Adaptive nonlinear system identification with echo state networks[END_REF] on a time series prediction task.

Related to the training of the ESN, the choice of the number of hidden neurons is of major importance when designing the network. One must remind that there are as many weights to learn as the number of neurons in the reservoir. As said previously, a large number increases the ability to represent a complex dynamics. However, if this number is too large (with respect to the amount of training data), the generalisation of the representation is not ensured (this is commonly known as the bias-variance trade-off). In other words, overfitting problems might appear.

The overfiting problem is underlined in [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]. In this tutorial, the authors propose to reduce the size of the reservoir. Yet, this may lead to poor results if the dynamics is complex. Other solutions like regularisation have also been proposed. In [START_REF] Wyffels | Stable output feedback in reservoir computing using ridge regression[END_REF], a ridge regression is used to solve a time series prediction task. Also, results with several schemes of pruning and regularisation have also been published in [START_REF] Dutoit | Pruning and regularization in reservoir computing[END_REF] and tested on two different problems. Notably, the Least Angle Regression (LAR) algorithm [START_REF] Efron | Least angle regression[END_REF], which performs an L1 regularisation, was tested (this approach performs somehow feature selection).

Yet, since the underlying dynamics is rarely known, a large number of neurons in the reservoir is preferable by default. Designing an approach where this number could be large is thus of interest. The work presented here proposes to increase the robustness to overfitting when adding neurons in the reservoir. To do so, a projection of the internal state of the network (the activity of the hidden neurons) on a space of a smaller dimension is added. The projection used is a random projection [START_REF] Achlioptas | Database-friendly random projections: Johnson-Lindenstrauss with binary coins[END_REF] and the data from the original space are projected on a random hyperplane. Under particular conditions, the distance between the points from the original space to the new ones are preserved with a given probability. The advantage of such a projection is that the dimension of the new space does not rely on the dimension of the initial space but only on the amount of data.

This new method has the advantage of decreasing the number of parameters to learn while keeping the number of hidden neurons large. A complex dynamics can thus be caught and the training can be made with a reasonnable amount of data. In what follows, a brief review of ESNs and random projections is proposed in Sec. 2. Then, in Sec. 3, the experimental settings are presented. Finally, in Sec. 4, the combination of ESNs and random projections is tested on a time series prediction task. The results are compared to the ones returned by an ESN trained with RLS and LAR algorithms.

ESN COMBINED WITH RANDOM PROJECTIONS

Principle of ESNs

Echo State Networks (ESNs) [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF] are Recurrent Neural Networks (RNNs) where only the weights to the output layer are learnt. The other connections are chosen randomly in order for the network to have the echo state property. This property states that if the internal state of the network is equal to another one after the network has been fed with two sequences of inputs, then the input sequences are identical. The reservoir of the ESN must be rich enough to match each history of inputs to a different internal state. The ESNs are very interesting networks since they avoid the difficult task of learning the recurrent connections.

Let u(k) ∈ R N i be a column vector representing the input at time k (Ni being the number of input neurons), x(k) ∈ R N h be a column vector representing the activity of the hidden neurons (N h being the number of hidden neurons), and y(k) ∈ R No be a column vector representing the activity of the output (No being the number of output neurons). Let W in be a (N h × Ni)-matrix containing the synaptic weights of the connections between the input neurons and the hidden layer, W hid a (N h × N h )-matrix containing the weights of the connections in the hidden layer and W out a (No ×N h )-matrix containing the weights of the connections from the hidden layer to the output. The activation of the hidden neurons and the output neurons is updated as follows:

x(k + 1) = f (W hid x(k) + W in u(k + 1)) y(k + 1) = f out (W out (u(k + 1), x(k + 1)))
with f the activation function of the hidden layer (for example f = tanh) and f out the activation function of the output layer (for example f out being the identity). The vector (u(k + 1), x(k + 1)) T is the concatenation of the vector u(k + 1) and x(k + 1). The input matrix W in is a possibly sparse random matrix and W hid is a sparse matrix. No precise conditions have been determined to choose the coefficients of the matrix W hid but for the choice of the spectral radius α of W hid . The spectral radius is defined as α = max i |λi|, λi being the eigenvalues of W hid . A spectral radius greater than one will just induce the echo state property not to be ensured [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF]. In practice, a network which W hid matrix has a spectral radius smaller than one seems to often exhibit the echo state property. In Fig. 1 the description of the ESN is summarised.

The main advantage of the ESN approach is that most of the synaptic weights are chosen randomly. Only the output weights are learnt. To do so, a linear regression can be used to learn the is provided to the network. The training error to be minimised is

W out matrix. A sequence of n ∈ N pairs (u teach (k), y teach (k)) Input u(k) [•] MN i ×1 Hidden x(k)     • • • . . . •     MN h ×1 Output y(k) [•] MN o ×1 W hid MN h ×N h W in MN h ×N i W out MN o ×N h Input u(k) [•] MN i ×1 Hidden x(k)     • • • . . . •     MN h ×1 Projection x Π (k) • • • Mm×1 Output y(k) [•] MN o ×1 W hid MN h ×N h W in MN h ×N i Π Mm×N h W out MN o ×m
train(k) = y teach (k) -W out (u teach (k), x(k)).
Here are the steps of the ESN training. First, the input matrix W in and the matrix of hidden connections W hid are randomly built. Then, for all of the examples of the training set, S = (u teach (k), y teach (k)), k ∈ 0..n, the network is fed with u teach (k) and the output y(k) = W out (n -1)x(k) is computed. Each time an example is provided, the weights are updated following a Widrow-Hoff rule:

W out (k) = W out (k -1) + K k (y teach (k) -y(k)).
The gain K k is computed so that the empirical cost k i=0 train(i) 2 is minimised. The algorithm used here is the Recursive Least Squares (RLS) algorithm.

Since finding a good compromise between the complexity of the model (which memory depth is rarely known) and the amount of training data is not so obvious, it can be of interest to find a new approach to get rid of this constraint. This approach would present the advantage of allowing a large number of neurons in the reservoir without paying attention to overfitting problems. The idea developed here is to project the x vector into a space of smaller dimension.

Random projections

As presented in the previous section, choosing the number of hidden neurons may be difficult in situations where neither the amount of training data is big nor a good intuition about the complexity of the dynamics of the system is known. The method proposed here allows the number of neurons to be increased while avoiding overfitting.

The method is a projection into a space of a smaller dimension with some guarantees about the loss of information. The method is called random projection since the coefficients of the projection matrix are chosen randomly (under certain conditions). The weights of the output matrix W out are not learnt direclty from the x(k) vector anymore but with the projection of the x(k) vector on this basis.

The results about the random projections are based on a theorem by Dimitri Achiloptas [START_REF] Achlioptas | Database-friendly random projections: Johnson-Lindenstrauss with binary coins[END_REF] based itself on a lemma by Johnson and Lindenstrauss [START_REF] Johnson | Extensions of Lipschitz maps into Banach spaces[END_REF]. The theorem states that it is possible to find a projection with a controlled loss of information.

Theorem Let M be an arbitrary set of n points in R m represented in a n × m matrix A. Given , β > 0, let K ,β be a constant and m0 = K ,β log(n).

For integer p > p0, let R be a m × p random matrix with R(i, j) = rij, where rij are independant random variables: rij = + √ 3 with probability (w. p.) 1/6, -√ 3 w. p. 1/6, 0 w. p. 2/3. Let Π = 1 √ p AR and let f : R m → R p map the i th row of A to the i th row of Π. W. p. at least 1 -p -β , for all u, v ∈ M : This result states that finding a projection from the initial space to a space with a smaller dimension while losing a bounded amount of information is possible. The information loss depends on the dimension of the projection which itself depends on the amount of data n (log(n) dependency). Such a projection is interesting since its dimension does not depend on the dimension of the initial space.

(1 -)||u -v|| 2 ≤ ||f (u) -f (v)|| 2 ≤ (1 + )||u -v|| 2 W hid MN h ×N h Input u(k) [•] MN i ×1 Hidden x(k)     • • • . . . •     MN h ×1 Projection x Π (k) • • • Mp×1 Output y(k) [•] MN o ×1 W hid MN h ×N h W in MN h ×N i Π Mp×N h W out MN o ×p
By means of this theorem, a large number of neurons can be used in the hidden layer. The x vector is then projected into a space of smaller dimension. Eventually, the coefficients of the model are learnt from the projected data. The computational complexity thus decreases and the risk of overfitting decreases too since the number of parameters, for a same amount of training data, is lower. Fig. 2 sums up how ESNs and random projections are combined.

ESNs combined with RPs are tested in the next section for a time series prediction task, a standard benchmark where ESNs have already been proven efficient [START_REF] Jaeger | Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication[END_REF][START_REF] Jaeger | Adaptive nonlinear system identification with echo state networks[END_REF].

EXPERIMENTAL SETTINGS

The time series used for the experiments is described first. ESNs parameters and the metrics for the evaluation of the results follow.

Description the time series

A 10 th order non-linear NARMA time series is used for the experiments. This series is the same as the one described in [START_REF] Jaeger | Adaptive nonlinear system identification with echo state networks[END_REF] and was the first time used in [START_REF] Atiya | New results on recurrent network training: Unifying the algorithms and accelerating convergence[END_REF]. The system is driven from an independent and identically distributed uniform input e from [0 , 0.5]:

d(k + 1) = tanh 0.3d(k) + 0.05d(k) 9 i=0 d(k -i) + 1.5e(k -9)e(k) + 0.1 .
To run the experiment, input sequences (e teach (k), d teach (k)) have been generated. The first 200 samples have been discarded to let the time series stabilises. An overview of the series and an example of prediction returned by the ESN is given in Fig. 3. Although the series is noisy, the prediction follows the actual series even if not perfectly.

Description of the ESN

The parameters of the ESN are the following: Ni = 1 and No = 1 (since it is a 1-dimension time series), N h varies. The input unit u is linked to the reservoir with a random vector W in sampled from a uniform distribution in [-0.1 , 0.1]. The reservoir matrix W hid has a 5% connectivity and α set to 0.8. The number of hidden neurons varies during the experiments to check the influence of the reservoir size on the quality of the prediction. The weights between ouput unit y and the reservoir are learnt. When the ESN is fed with new samples, the first 500 data are discarded to let the reservoir of the ESN reach the steady state.

Evaluation of the performance

The Mean Squared Error (MSE) between the testing data ytest and the prediction ŷ of the network over several trials serves as a performance metrics:

M SE = 1 T T i=1 (ytest(i) -ŷ(i)) 2 .
The amount of training data n is set to 100. The amount of testing data T is set to 100. Each of the curves presented is an average over the MSE computed for 50 different training sets and different ESNs (for each new training set, new weights for W in and W hid are generated). The associated variance also appears on the graphs. Both the abscissia and ordinates are in log scale.

RESULTS

The following experiments illustrate the advantages of adding a random projection to ESNs to solve a time series prediction task. The first advantage is a better prediction for a given number of parameters. Indeed, since the number of parameters to be learnt is determined by the projection dimension, the number of hidden neurons can be increased. In consequence, a more complex dynamics can be modelled while avoiding overfitting problems. Second, the proposed contribution is an online method which can be of interest while dealing with real-time systems. A comparison to RLS without random projections and to the Least Angle Regression (LAR) [START_REF] Efron | Least angle regression[END_REF] algorithm is proposed. LAR has been chosen because it particularly fits for problems where the number of training data is smaller than the number of parameters to be learnt. One can notice that an L2 regularisation is used in the RLS case while it is an L1 regularisation in the LAR case and that LAR is a batch algorithm. The first experiment presented in Fig. 4 shows how the MSE evolves when the number of parameters p increases. The performance of the RLS is compared to the one of the RLS combined with random projections. In the first case, the number of parameters to learn p is equal to the number of neurons in the reservoir. In the second case, the number of parameters p to learn is equal to the projection dimension. One has to pay attention that the projection dimension is never larger than the number of hidden neurons which explain why some points are missing in the figure. In Fig. 4, the error with the RLS algorithm decreases when the number of hidden neurons increases to reach a minimun (N h = 50). Until a reasonnable number of neurons in the reservoir, the representation of the series is richer and thus of better quality. Then, the error increases to stabilise. The regularisation performs quite well since the error does not become huge even when the number of parameters is far larger than the number of training data.

When an RP is added for the approximation, for a fixed dimension, say for example 10, the error decreases when the number of hidden neurons increases to reach a constant level. This is of interest to use RP since for a fixed number of parameters p, when not chosen too large (5 or 10), the random projections allow a better learning. The MSE is lower than when the RLS is direclty ran with the activities of the hidden neurons. One can notice that adding neurons in the hidden layer does not degrade performance.

Using a random projection after the ESN is thus interesting: for an equivalent small number of parameters to learn (p < 100), the testing error is lower.

The Fig. 5 presents this aspect. The projection dimension is fixed and the number of hidden neurons increases. Curves are plotted for different projection dimensions. The results got with the RLS algorithm are recalled. It appears that if the projection dimension is sufficient enough (p > 100), the quality of the prediction is nearly equivalent: the difference between the best and the worst trial is at most 0.001. Moreover, the number of hidden neurons does not influence the quality of the learning. It can be chosen very large with respect to the amount of training data (N h > 1000) without increasing the error. The overfitting problem is thus avoided.

The two presentations of the results show the interest of using random projections: for a fixed number of parameters to learn, the error is lower and the choice of the number of parameters is made simple. Indeed, this number can be large (to ensure the dynamics of the system to be caught) without overfitting problems.

The third aspect presented is the comparison between LAR and RLS combined with random projections. This comparison is interesting since each of these algorithms somehow performs feature selection. The first one because only the features which bring the most important information are selected. The second one because the data is projected into a space of smaller dimension. The results are shown in Fig. 6. The error made with LAR tends to deacrease when the number of neurons in the reservoir increases until a minimum is reached (N h = 50). Then, when the number goes on increasing, the error increases too until N h = 300 where a plateau is reached. When the random projection is used and the learning is made with the RLS algorithm, the results are of better quality when the number of hidden neurons tends to be large. Moreover, if we compare the LAR and the RLS results (Fig. 4), the L1 regularisation (LAR) is less efficient than the L2 regularisation (RLS).

CONCLUSION

Some encouraging results concerning time series prediction have been obtained by combining echo state networks and random projections. This method proposes a solution to limitate the overfitting problems. Indeed, the number of parameters to learn is decreased without needing the number of neurons in the reservoir to be decreased too (a large number of hidden neurons guarantees time dependencies to be caught). Consequently, the computational costs are also lowered. Moreover, the proposed approach has the advantage of being online.

Fig. 1 .

 1 Fig. 1. Structure of an ESN. For the example, Ni = 1 and No = 1.

Fig. 2 .

 2 Fig. 2. Combination of the ESN random projections. With xΠ(k) = Πx(k). For this example, Ni = 1, m = 3, and No = 1.

Fig. 3 .

 3 Fig. 3. Overview of the time series and an example of prediction returned by the ESN.

Fig. 4 .

 4 Fig. 4. Comparison of the ESNs performance with the ones of the ESN and RPs versus the number of parameters p to learn. The algorithm used for the learning is RLS.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Comparison of the ESNs performance with the ones of the ESN and RPs versus the number of neurons in the reservoir. Several projection dimensions are tested. The algorithm used for the learning is RLS.