Layla El Asri
email: layla.elasri@maluuba.com

Bilal Piot
email: bilal.piot@univ-lille1.fr

Matthieu Geist
email: matthieu.geist@centralesupelec.fr

Romain Laroche
email: romain.laroche@orange.com

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

Score-based Inverse Reinforcement Learning

Keywords: Reinforcement Learning, Inverse Reinforcement Learning, Markov Decision Processes, Learning from Demonstration, Spoken Dialogue Systems

This paper reports theoretical and empirical results obtained for the score-based Inverse Reinforcement Learning (IRL) algorithm. It relies on a non-standard setting for IRL consisting of learning a reward from a set of globally scored trajectories. This allows using any type of policy (optimal or not) to generate trajectories without prior knowledge during data collection. This way, any existing database (like logs of systems in use) can be scored a posteriori by an expert and used to learn a reward function. Thanks to this reward function, it is shown that a near-optimal policy can be computed. Being related to least-square regression, the algorithm (called SBIRL) comes with theoretical guarantees that are proven in this paper. SBIRL is compared to standard IRL algorithms on synthetic data showing that annotations do help under conditions on the quality of the trajectories. It is also shown to be suitable for real-world applications such as the optimisation of a spoken dialogue system.

INTRODUCTION

Reinforcement Learning (RL) has seen a growing interest this last decade both from a theoretical and a practical point of view, especially in the domain of Human-Machine Interaction (HMI) where it is now part of the state of the art [START_REF] Singh | Reinforcement learning for spoken dialogue systems[END_REF][START_REF] Lemon | Data-Driven Methods for Adaptive Spoken Dialogue Systems: Computational Learning for Conversational Interfaces[END_REF][START_REF] Young | POMDP-based statistical spoken dialog systems: A review[END_REF]. In this paradigm, an agent steps from states to states in an environment of unknown dynamics by selecting actions. The action selection process is encoded in a so-called policy, which is a mapping between states and ac-tions. Learning occurs via a reward function providing a numerical feedback to the agent after each transition from one state to another. An optimal policy, which should be the result of the learning process, leads to the highest cumulative reward in average. This paradigm is very attractive to learn optimal behaviors in a model-free manner, which allows dealing with complex systems (such as HMI where the human is in the loop and hard to model). Yet, correctly defining the reward function so that it leads to the desired behavior remains a tricky task. For this reason, Inverse Reinforcement Learning (IRL) [START_REF] Russell | Learning agents for uncertain environments[END_REF] received a lot of attention recently, including in HMI [START_REF] Paek | Reinforcement learning for spoken dialogue systems: Comparing strengths and weaknesses for practical deployment[END_REF][START_REF] Boularias | Learning the reward model of dialogue pomdps from data[END_REF][START_REF] Pietquin | Inverse reinforcement learning for interactive systems[END_REF], and many algorithms can now be found in the literature [START_REF] Abbeel | Apprenticeship learning via inverse reinforcement learning[END_REF][START_REF] Boularias | Relative entropy inverse reinforcement learning[END_REF][START_REF] Klein | Inverse reinforcement learning through structured classification[END_REF]. IRL consists of finding the reward optimised by an expert that provides optimal demonstrations. Yet, these algorithms require numerous expert demonstrations, some of them have to solve the direct RL problem many times [START_REF] Abbeel | Apprenticeship learning via inverse reinforcement learning[END_REF] or also require many non-expert (and if possible random) trajectories to estimate the dynamics [START_REF] Boularias | Relative entropy inverse reinforcement learning[END_REF]. It is often hard and potentially impossible to comply with these requirements, especially when dealing with human experts. Collecting enough expert demonstrations is time consuming and expensive. Besides, it is also hard for a person to consciously act sub-optimally on demand and, above all, randomly so as to cover the widest range of possible behaviors. This is even more difficult when facing other persons (in the case of HMI problems).

In this paper, we study a different setting for IRL based on scored trajectories and especially an associated algorithm based on the minimisation of the distance between actual and predicted scores [START_REF] Asri | Reward shaping for statistical optimisation of dialogue management[END_REF]. In this paradigm, a database of trajectories is first collected. Each trajectory can be obtained by applying any policy (optimal or not). Each trajectory is then annotated a posteriori in terms of global performance by an external expert, who might be error-prone. This is a very practical setting where no expert is needed during the collection of trajectories and there is no requirement on their quality. The expert's job is reduced since s/he only has to rate trajectories instead of providing optimal decisions in each state. It also applies in tasks where only subjective ratings can be obtained, which is often the case for HMI [START_REF] Möller | Predicting the quality and usability of spoken dialogue services[END_REF],

or where a set of unreliable data already exists (such as Mechanical Turk data collection).

Compared to preference-based reinforcement learning [START_REF] Akrour | Preference-based policy learning[END_REF], this method doesn't require a user to explicitly compare each pair of trajectories in the database and therefore results in a lighter effort for the expert. Yet, each trajectory is implicitly compared to all the others through the scores which results in additional information. Moreover, preference-based RL learns a policy, not a reward whereas learning a reward is the purpose of this paper. Indeed, we argue that IRL is a technique which has many advantages such as the transferability to other environments [START_REF] Piot | Learning from demonstrations: is it worth estimating a reward function?[END_REF] or the ability for the agent to learn online by reinforcement after the reward has been computed, making the policy improve with time. In addition, in the domain of HMI, collecting expert ratings is common and well-researched [START_REF] Möller | Predicting the quality and usability of spoken dialogue services[END_REF][START_REF] Asri | DINASTI: Dialogues with a Negotiating Appointment Setting Interface[END_REF] It is important to notice that global scores cannot be used directly as rewards (for instance on the final state). Indeed, similar scores can be associated to very different trajectories, finishing on the same state, but having different lengths for instance. Scores must thus be considered as a value associated to the global trajectory. In practice, the IRL problem is cast into a special case of least-square regression and thus results in a very efficient method (renamed SBIRL for Scorebased IRL instead of Distance-Minimisation IRL), robust to noisy scoring and coming with theoretical guarantees. Especially, the error propagation analysis is an original contribution of this paper which is accompanied by a finite-sample analysis relying on very recent results on linear least-square estimation. The practical properties of the algorithm are studied on synthetic problems as well as on a real-world application: the optimisation of a spoken dialogue system strategy.

FORMAL DESCRIPTION

A Markov Decision Process (MDP) is a tuple {S, A, P, γ, r} where S is the state space, A is the finite action space, P (ds |s, a) the Markovian transition kernel on S, γ ∈ (0, 1) the discount factor and r : S → R the (bounded) reward function. A policy is a mapping π : S → A. For a policy π, we define the transition kernel Pπ(ds |s) = P (ds |s, π(s)). The quality of a policy is quantified by the related value function, that associates to each state the expected discounted sum of rewards received from starting in state s and following the policy π:

v r π (s) = E ∞ t=0 γ t r(St)|S0 = s, St+1 ∼ Pπ(.|St) .
An optimal policy π r * (with respect to the reward r) is such that the related value function v r π r * = v r * satisfies componentwise v r * ≥ v r π , for any policy π. In this work, the true reward r is unknown and has to be estimated. To do so, we adopt a linear parameterisation of the reward and we assume the availability of a set of trajectories scored by a human expert. More formally, the reward is parameterised by a feature vector φ :

S → R d , r θ (s) = θ φ(s).
The available data is a set

D = {(hi, vi) 1≤i≤n }, where hi = (s i 0 , . . . , s i T i ) = (s i j ) T i j=0
is a trajectory of length Ti +1 and vi is the score (seen as the discounted sum of rewards of hi) given by a human expert for this trajectory. Notice that we do not make any specific assumption about how the trajectory is generated (except that it is a trajectory obtained by applying successively Ti actions on encountered states, next states being drawn according to the dynamics). Trajectories in the dataset are assumed independently and identically distributed (which does not mean that states of a given trajectory are independent, obviously). For a given trajectory h = (st) T t=0 and a given reward r θ , the discounted sum of rewards can be written as

T t=0 γ t r θ (st) = θ µ(h) with µ(h) = T t=0 γ t φ(st).
We interpret the scores given by the human expert as noisy estimates of the discounted sum of rewards of the associated trajectories. The expert is thus not assumed to give perfect scores but is error-prone. The problem thus consists of regressing the scores vi on the mappings of the histories µ(hi). Formally, we require that the underlying procedure asymptotically minimizes the risk based on the classic 2-loss:

R(θ) = E[(V -θ µ(H)) 2 ]. (1)
The joint distribution on V (scores) and µ(H) (mappings of histories) is imposed: it is the distribution used to sample the dataset D. Given any procedure that asymptotically minimises risk (1), and denoting by θn the estimate computed from the dataset D, we obtain an estimate r θn of the reward. For example, the SBIRL algorithm estimates the reward by solving the following linear least-squares problem:

θn = argmin θ∈R d 1 n n i=1 vi -θ µ(hi) 2 = 1 n n i=1 µ(hi)µ(hi) -1 1 n n i=1 µ(hi)vi, (2)
assuming that the matrix 1 n n i=1 µ(hi)µ(hi) is invertible (ordinary least-squares). If this is not the case, one can use for example 2-regularisation [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF] or 1-regularisation [START_REF] Tibshirani | Regression shrinkage and selection via the LASSO[END_REF], among others. Hence, given a set of trajectories scored by a human expert, we have a reward estimate r θn . In the next section, we analyze the quality of this reward.

ANALYSIS

To study the estimate r θn , we assume that the scores provided by the expert correspond to a discounted cumulative sum of rewards, up to some noise η, for a reward function that lies in the hypothesis space H = {θ φ(s), θ ∈ R d }.

Assumption 1

There exists a vector parameter θ * ∈ R d and a centered noise η such that for any trajectory h = (st) T t=0 and any associated score v, one has

v = T t=0 γ t r θ * (st) + η(h).
Under this assumption, the minimiser of risk (1) is

θ * = argmin θ∈R d R(θ) = A -1 b, with A = E[µ(H)µ(H) ] and b = E[µ(H)V ],
assuming that the matrix A is invertible. First, we assume that the estimator θn satisfies (with high probability) that R(θn) ≤ R(θ * ) + , for some error . Then, we instantiate the related results in the case of ordinary least-squares.

Propagation of errors

First, we want to control the risk Eν [(r θ * (S) -r θn (S)) 2 ] for some distribution ν over states.

Theorem 1 (Propagation of errors: rewards) Write λm the minimum eigenvalue of E[µ(H)µ(H) ] and λM the maximum eigenvalue of Eν

[φ(S)φ(S) ]. Assume that θn satisfies R(θn) ≤ R(θ * ) + , then the associated reward satisfies Eν [(r θ * (S) -r θn (S)) 2 ] ≤ λM λm .
Proof. The proof is provided in the appendix.

Theorem 1 shows that if is small, the estimated reward r θn will be close to the true reward r θ * . However, we are ultimately interested in computing an optimal policy π r θn * (optimal policy relatively to the reward r θn ) for the control problem at hand. The question we answer now is: how close to optimal is the policy π r θn * respectively to the one of the true reward π * r θ * ? To do so, we bound the quantity

E ν [(v r θ * * (S) -v r θ * π r θn * (S)) 2 ],
states being sampled according to an arbitrary distribution ν of interest (possibly different from the distribution ν).

Before stating the result, we provide some notations. For a distribution ν and a function f ∈ R S , we write

νf = Eν [f (S)].
For a stochastic kernel Q, we have

[Qf ](s) = E Q(.|s) [f (S )].
Therefore, νQ is a distribution such that

EνQ[f ] = Eν [Qf ] = ES∼ν [E S ∼Q(.|S) ][f (S )].

Theorem 2 (Propagation of errors: values)

For a policy π, define Cπ as the smallest coefficient such that

(1 -γ)ν (I -γPπ) -1 ≤ Cπν.
Assume that θn satisfies R(θn) ≤ R(θ * ) + , then the associated value satisfies

E ν [(v r θ * * (S) -v r θ * π r θn * (S)) 2 ] ≤ 2 Cπ * + C π r θn * (1 -γ) 2 λM λm .
Proof. The proof is provided in the appendix.

Theorem 2 shows that, if the error is small, the optimal policy with respect to the learned reward function will be close to the optimal policy with respect to the unknown reward (closeness being measured in terms of value functions). There is an additional multiplicative term compared to Th. 1. The concentrability coefficient Cπ measures the dissimilarity between the distribution of data ν and the distribution (1 -γ)ν (I -γPπ) -1 , the γ-weighted occupancy measure induced by policy π when the initial state is sampled from ν (the distribution of interest for controlling the value function) and the coefficient (1 -γ) -1 is the average optimisation horizon. Both terms are standard when bounding value functions.

A finite sample analysis

In this section, we provide a finite sample analysis in the case where θn ∈ R d is the ordinary linear least-squares estimate of Eq. ( 2), based on the bound of [START_REF] Hsu | Random design analysis of ridge regression[END_REF]. For this, we need some technical assumptions.

Assumption 2

The noise is subgaussian: there exists σ ≥ 0 such that, almost surely, for all λ ∈ R,

E[e λη(H) |H] ≤ e λ 2 σ 2 2
.

Moreover, there exists a finite ρ ≥ 1 such that, almost surely,

A -1 2 µ(H) 2 √ d ≤ ρ.
The subgaussian assumption is easily satisfied if the scores provided by the expert and the basis functions are bounded (in this case, the noise is bounded, and a bounded random variable is subgaussian [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF]). The other assumption corresponds to condition 1 of [START_REF] Hsu | Random design analysis of ridge regression[END_REF]. For example, if for any s ∈ S we have

φ(s) ∞ ≤ φmax,
then almost surely we have

µ(H) 2 ≤ √ d 1 -γ φmax
and the assumption holds for

ρ ≥ φmax (1 -γ) √ λm .

Corollary 1

Let δ be such that

ln 9 δ > max(0, 2.6 -ln d).
If assumptions 1 and-2 hold and if

n ≥ 6ρ 2 d ln 9d δ ,
then with probability at least 1 -δ we have

Eν [(r θ * (S) -r θn (S)) 2 ] ≤ λM λm σ 2 (d + 2 d ln 9 δ + 2 ln 9 δ ) n + o 1 n .
Under the same assumptions, w.p. at least 1 -δ:

E ν [(v r θ * * (S) -v r θ * π r θn * (S)) 2 ] ≤ 2 Cπ * + C π r θn * (1 -γ) 2 λM λm σ 2 (d + 2 d ln 9 δ + 2 ln 9 δ ) n + o 1 n .
Proof. The proof is provided in the appendix.

This result mainly tells that in the case of ordinary leastsquares, we have = O( d n ), and thus the rates for the errors of the reward and of the value (of the optimal policy) are the same, up to the additional constants arising because of error propagation. Since d is the number of features and n the number of trajectories, this gives an idea of the amount of data required to attain a given accuracy (however, notice that the concentrability coefficients can be hardly estimated, a standard problem in reinforcement learning). A wiser (but much more difficult) analysis would take into account the number of transitions instead of the number of trajectories, yet this is beyond the scope of this paper.

EXPERIMENTS

Experiments are organised in two sets. The first set is executed on synthetic data, i.e. randomly-constructed finite MDPs called Garnets [START_REF] Archibald | On the generation of Markov decision processes[END_REF][START_REF] Piot | Learning from demonstrations: is it worth estimating a reward function?[END_REF][START_REF] Piot | Boosted and reward-regularized classification for apprenticeship learning[END_REF]. We lead three main canonical experiments aiming at supporting the theoretical results. First, we verify the improvement in performance of SBIRL when the number of trajectories grows. Second, we evaluate the sensitivity to the noise η. Third, we determine if sampling the trajectories with an expert or a random policy influences the performance of SBIRL and compare with standard IRL methods provided with the same data. In a second set of experiments, we investigate the suitability of SBIRL for real-world problems by applying it to the optimisation of spoken dialogue system management.

Experiments on Garnets

In this section, experiments are done on a set of artificially built MDPs so as to support the theoretical findings. Testing SBIRL requires three key elements: first, to build a finite MDP; second, to generate trajectories of random lengths in this MDP; and finally, to score these trajectories with a given noise η. Garnets [START_REF] Archibald | On the generation of Markov decision processes[END_REF] are an abstract class of finite MDPs, easy to build. Here, we consider a special case of Garnets specified by three parameters: (NS , NA, NB). Parameters NS and NA are respectively the number of states and of actions. Thus, S = (si) N S i=1 and A = (ai) N A i=1 are, respectively, the state and action spaces. The parameter NB (NB ≤ NS ), called the branching factor, defines for each pair (s, a) the number of next states. NB states are drawn uniformly and without replacement from S and form the set of next states of (s, a) noted Ss,a = (s i ) N B i=1 . Then, to define the Markovian kernel P , for each state-action (s,a), we draw randomly and uniformly in [0, 1] NB -1 cutting points. Let us note (pi) N B -1 i=1 this set of cutting points sorted in increasing order, p0 = 0 and pN B = 1. To completely define the dynamics, one assigns P (s i |s, a) according to the following rule: ∀i ∈ {1, . . . , NB}, P (s i |s, a) = pi -pi-1. Finally, for each state s ∈ S, the reward r(s) is drawn randomly and uniformly in [0, 1], and γ = 0.9. This choice imposes that the non-perturbed score of a trajectory h = (st) T t=0 is bounded as follows

0 ≤ T t=0 γ t r(st) ≤ 1 1 -γ = 10.
As we choose finite MDPs, a canonical choice of features φ is the tabular basis φ : S → R N S where φ(s) ∈ R N S is a vector which is null excepted in s where it is equal to 1.

The strategy chosen to generate a trajectory h of random length in a Garnet following a policy π consists of first choosing randomly and uniformly a starting state s0 ∈ S. Then, starting from s0 we apply, with probability (1 -p) (p ∈ (0, 1)), the policy π to the current state si in order to go to the next state si+1 and with probability p we stop the trajectory. Doing so, we obtain a trajectory h = (st) T t=0 where the length of the trajectory T + 1 is a geometrically distributed random variable: for all k ∈ N * ,

P r(T = k) = (1 -p) k-1 p, E[T ] = 1 p and V ar[T ] = 1 -p p 2 .
We choose this strategy because it is easy to implement and allows us to control the mean of the length of the trajectories.

Finally, to give a score to a trajectory h = (st) T t=0 , we proceed as follows. Let η be a Gaussian distributed and real valued random variable with mean ν = 0 and variance σ 2 (η ∼ N (0, σ 2 )), the score of the trajectory h is

v(h) = T t=0 γ t r(st) + η = T t=0 γ t r(st) + η + ηq,
where x is the nearest integer from x and ηq is a quantification noise which is supposed to be a uniformly distributed random variable of mean µ = 0 and variance σ 2 = 1 12 . We deliberately add a quantification noise in order to model the fact that experts are often asked to use an integer scale. More precisely, we have seen that the choice of the reward r and the discount factor γ imposes, for a trajectory h = (st) T t=0 , that the non-perturbed score is bounded: 0 ≤ T t=0 γ t r(st) ≤ 1 1-γ = 10. Thus, we have constructed the score v(h) such that it models the choice of a score in the set of integers {0, 1, . . . , 10}.

Performance evaluation of SBIRL

So as to evaluate the SBIRL performance, we want first to measure the mean performance of this algorithm over several Garnets when the number of trajectories, NT , grows. To do so, we create (Gq) N G q=1 Garnets of size (NS = 100, NA = 5, NB = 10) where we compute an optimal policy π e q w.r.t the real reward rq of Gq via the policy iteration algorithm. For each Gq and each iteration j ∈ {1, . . . , NIt} (NIt = 10), we generate NT trajectories (hi) N T i=1 according to a random policy πr (at each state s, the probability to choose action a is 1 N A ) with p = 0.01 such that the mean length is around 100. We compute, for each trajectory the score v(hi) with η ∼ N (0, σ 2 = 1). Here, the standard deviation of the noise represents 10% of the maximum non perturbed-score achievable by a trajectory. Thus, for each Gq, each iteration j and a given number of trajectories NT , we obtain a set (hi, vi = v(hi)) N T i=1 that we use as an input for SBIRL which outputs the reward rq,j. If π r q,j * is the optimal policy w.r.t rq,j, then the performance of SBIRL is measured via the normalised error between the value functions (w.r.t rq) of π e q and π r q,j * :

Tq,j(NT ) = v rq π e q -v rq π r q,j * 2 v rq π e q -1 2 ,
where . 2 is the Euclidean norm. The lower the error, the better the performance. Finally, the mean error T (NT ) is the mean of the (Tq,j(NT )) 1≤j≤10 1≤q≤50 . In Fig. 1, we plot the SBIRL error, T (NT ), where NT varies from 100 to 1000 and the error corresponding to the random policy. We observe, as predicted by the analysis, that the error converges to zero as NT grows.

Another interesting aspect to evaluate is the tolerance of our method to the noise η ∼ N (0, σ) over a number of Garnets. In order to evaluate this aspect, we realise exactly the same experiment as before except that NT is fixed to 500 and σ (standard deviation of the noise) is now the varying parameter. For each Gq, each iteration j and for a given σ, SBIRL outputs the reward rq,j and the error Tq,j(σ) of SBIRL is defined as before. The mean error T (σ) is the mean of the (Tq,j(σ)) 1≤j≤10 1≤q≤50 . In Fig. 2, we plot the SBIRL error, T (σ), where σ varies from 0 to 5 which represents 50% of the maximum nonperturbed score given to a trajectory. We observe that the performance is good when the noise is low and deteriorates as the noise gets bigger which is also observed in the bound of Corollary 1. We also remark that the standard deviation of the error, which appears as a shade in Fig. 2, gets bigger as the noise increases.

Comparison to IRL

This experiment is designed to provide insights about how unreliable (w.r.t an optimal policy) the data should be to require annotation, preventing the use of standard IRL. To compare IRL to SBIRL, we need to generate near-optimal trajectories to simulate a wide range of expertise levels. To do so, we consider a noisy expert where the noise is controlled by a parameter β such that the generated policy could range from optimal to random policies. More precisely, the trajectories (hi) N T i=1 are sampled according to a noisy-expert policy π n (β) where β ∈ [0, 1] and given as an input to SBIRL and two IRL algorithms which are Projection [START_REF] Abbeel | Apprenticeship learning via inverse reinforcement learning[END_REF] and SCIRL [START_REF] Klein | Inverse reinforcement learning through structured classification[END_REF]. To do so, βNS states are drawn randomly without replacement from S. On those states, π n (β) is random and on the others, π n (β) is optimal. We perform the same experiment as before except that NT is fixed to 300 and 500, NIt = 3 and the sampling policy, for each Gq, is π n q (β). For each Gq, each iteration j and for a given β, each algorithm A outputs the reward r A q,j with an error T A q,j (β). The mean error T A (β) is the mean of the (T A q,j (β)) 1≤j≤3 1≤q≤50 . In Fig. 3(a) and Fig. 3(b), errors are plot w.r.t β (proportion of random vs optimal policy) varying from 0 to 1 for NT = 300 and NT = 500 respectively. When NT = 300, the collected trajectories must be optimal on approximately 60% (β = 0.4) of states for Projection and 40% (β = 0.6) for SCIRL to achieve the same result as SBIRL. And, when NT = 500, the collected trajectories must be optimal on approximately 80% (β = 0.2) of states for Projection and 55% (β = 0.45) for SCIRL to achieve the same result as SBIRL. This shows that when NT grows, the data should be of high quality to reach the level of performance of SBIRL. So, according to the amount of available data, their reliability has to be assessed so as to decide whether directly using IRL or spending time for annotations and use SBIRL.

Application to Spoken Dialogue Systems

Here, SBIRL is applied to the optimization of a spoken dialogue system (SDS) implementing an appointment scheduling task. This experiment is mainly proposed to show the suitability of SBIRL for use in a real-world problem. Of course we do not have information about the quality of the scores or the noise and we do not study the sensitivity to these variables here.

An SDS is an HMI which converses in spoken natural language with users. Speech is transcribed to the machine by Automatic Speech Recognition (ASR). A dialogue is mod- To do so, it keeps track of the dialogue context. This process is modeled as an MDP where a state s encodes the current dialogue context and an action a represents a System Dialogue Act (SDA). An SDA describes the meaning of an utterance at the level of illocutionary force [START_REF] Austin | How to do Things with Words[END_REF]. Unlike [START_REF] Asri | Reward shaping for statistical optimisation of dialogue management[END_REF] which simulated dialogues, this experiment relies on interactions collected using an actual SDS. Scores are computed thanks to annotations given by real users. The collected corpus, DINASTI [START_REF] Asri | DINASTI: Dialogues with a Negotiating Appointment Setting Interface[END_REF], contains 1734 dialogues. For each dialogue, the user had a two-week calendar and was asked to set an appointment at one of the available slots on the calendar. There was always only one slot available for both the system and the user so each dialogue needed several turns of negotiation to find the convenient appointment. After each dialogue, users were asked to fill in a questionnaire and rate the dialogue on a scale of 1 to 10. We trained an ordinal logistic regression model on these ratings and computed the score of new dialogues based on the numbers of dialogue turns and ASR rejections, the average ASR confidence score and task completion.

The system Sys and the user Usr have 30 time slots in common. Each slot is composed of a day of the week, the day number and an hour. Usr and Sys have different preferences concerning these slots. This is modeled by uniformly drawing in [0, 1] two values, v Usr (sli) for Usr and v Sys (sli) for Sys, for each slot sli. These values are drawn at the beginning of each dialogue and Usr is set to only accept the 5 slots with the highest values. If Sys proposes one of these slots, Usr accepts it, otherwise s/he suggests one. Sys might ask Usr to confirm the last utterance it understood. For an implicit confirmation, Sys repeats the utterance and then moves on with the dialogue. If the utterance is correct, Usr does not say anything. Otherwise, s/he barges in before Sys has moved on and repeats the correct utterance. For an explicit confirmation, Sys asks Usr whether the understood utterance is correct or not. Usr then either confirms or repeats the correct utterance. Usr also repeats the last utterance if Sys informs her/him that the utterance was not understood (ASR rejection) or if one of the three elements of the slot was not understood. Each dialogue begins by a proposition by Sys. Then, during the dialogue, when there is no ASR problem, Sys can choose between 4 actions: accept a slot, refuse and propose another slot, ask for an implicit confirmation and ask for an explicit confirmation. A dialogue ends when one participant has accepted the other's proposition. Sys proposes slots in decreasing order of v Sys (sli). ASR confidence scores were simulated with a word error rate of 10%. The score was computed by drawing from a normal distribution N (-1, 1) in the case of a misrecognition (N (1, 1) otherwise) and by then applying the sigmoid function to this value. The threshold for ASR rejection was set to 0.5.

For evaluation purposes, SBIRL was applied to simulated dialogues with a uniformly random policy. Given sli the selected slot, task completion was defined as v Usr (sli) + 0.5 × v Sys (sli). The state space was composed of the latest User Dialogue Act (UDA, for instance ACCEPT) as well as the numbers of dialogue turns and ASR rejections, the average ASR confidence score and the value v Sys (sli) of the slot sli being currently discussed. The state space was built by applying entropy-based discretisation [START_REF] Fayyad | Multi-interval discretization of continuous-valued attributes for classification learning[END_REF] with these parameters on the random corpus and then by aggregating the re-sulting intervals. We performed 100 runs on 500, 1000 and 1500 dialogues. For each run, we computed the state space, ran SBIRL and then learnt a policy with the new reward function. On average, the state space had respectively 48.5, 71.8 and 79.5 states. The learning algorithm was SARSA with -greedy exploration, where was set to 0.01. The learning factor α was set to 0.1. We learnt a policy on 5000 dialogues and tested the learnt policy on 2000 dialogues. We compare learning with a handcrafted policy tested on 100000 dialogues which always asks the user for an implicit confirmation and then accepts a slot sli if v Sys (sli) ≥ 0.5 or if the number of dialogue turns is superior or equal to 15. Table 1 contains 95% confidence intervals for the estimated user satisfaction score and dialogue statistics with the handcrafted policy and with the policies learnt with SBIRL. The policy learnt with SBIRL improves the performance score by finding a better compromise between the number of dialogue turns and task completion. The policy learnt with SBIRL produces dialogues with, on average, 2 fewer dialogue turns than the ones resulting of the handcrafted policy and with sacrificing less than 10% of the task completion value. This is an important result since the average dialogue duration is well-known to highly influence user's satisfaction [START_REF] Walker | PARADISE: A framework for evaluating spoken dialogue agents[END_REF].

CONCLUSION

This paper studied the SBIRL algorithm which applies to a non-standard setting for Inverse Reinforcement Learning (IRL) relying on annotated trajectories. Trajectories can be collected using policies that are not assumed to be generated by an expert but are globally scored a posteriori by such an expert. This setting is much more suitable for many applications where datasets are already available but with unreliable quality. Expert's annotations are supposed to be noisy and robustness to noise was assessed both theoretically and practically. The theoretical analysis shows that the algorithm produces a policy that approximates the optimal one at a rate that increases linearly with the number of scored trajectories. Results of synthetic experiments were in accordance with theory and also showed that the noise in the scores is not too influential. These experiments also provided insights about how unreliable the data should be to require annotation, preventing the use of standard IRL methods. An application to spoken dialogue systems, using real data collected with an actual system, showed that this method is applicable in practice to learn more efficient interaction policies in a realistic domain (appointment scheduling). In future work, this methodology will be applied to large scale systems such as human-robot interaction.

APPENDIX

This appendix recalls the stated theoretical results and proves them. First, we recall the required assumptions.

Assumption 1

There exists a vector parameter θ * ∈ R d and a centered noise η such that for any trajectory h = (st) T t=0 and any associated score v, one has

v = T t=0 γ t r θ * (st) + η(h).

Assumption 2

The noise is subgaussian: there exists σ ≥ 0 such that, almost surely, for all λ ∈ R, E[e λη (H) |H] ≤ e λ 2 σ 2 2

.

Moreover, there exists a finite ρ ≥ 1 such that, almost surely,

A -1 2 µ(H) 2 √ d ≤ ρ.
We also recall that for a distribution ν and a function f ∈ R S , we write

νf = Eν [f (S)].
For a stochastic kernel Q, we have

[Qf ](s) = E Q(.|s) [f (S )].
Therefore, νQ is a distribution such that Proof. First, we use the fact that the risk is a quadratic form (and that Aθ * = b): Proof. To simplify the notations, write r = r θ * the unknown reward function (of associated optimal policy π * = π r θn *

EνQ[f ] = Eν [Qf ] = ES∼ν [E S ∼Q(.|S) ][f (S )].

A. PROPAGATION OF ERRORS: REWARDS

R(θ) = E[(v -θ µ(H)) 2 ] = θ Aθ -2θ b + E[V 2 ] = Aθ -b 2 A -1 + E[V 2 ] -b A -1 b = Aθ -b 2 A -1 + E[V 2 ] -2b θ * + θ * Aθ * = Aθ -b 2 A -1 + R(θ * ) = θ -θ * 2 A + R(θ * ).
) and r = r θn the estimated reward function (of associated optimal policy π = π r θn *

). With these notations, the quantity to be bounded is E ν [(v r π * (S) -v r π (S)) 2 ]. Componentwise, we have

v r π * -v r π = v r π * -v r π * + v r π * -v r π + v r π -v r π .
As π is optimal for r, we have v r π * -v r π ≤ 0, and a value function satisfies

C. FINITE SAMPLE ANALYSIS

Corollary 1

Let δ be such that ln Proof. The assumptions made allow using Theorem 1 of [START_REF] Hsu | Random design analysis of ridge regression[END_REF], which states that with probability at least 1 -δ we have Plugging this result in Th. 1 and 2 allows concluding.

Figure 1 :

 1 Figure 1: SBIRL error as NT grows

Figure 2 :

 2 Figure 2: SBIRL error as σ grows

Figure 3 :

 3 Figure 3: SBIRL vs IRL

Theorem 1 (

 1 Propagation of errors: rewards) Write λm the minimum eigenvalue of E[µ(H)µ(H) ] and λM the maximum eigenvalue of Eν [φ(S)φ(S) ]. Assume that θn satisfies R(θn) ≤ R(θ * ) + , then the associated reward satisfies Eν [(r θ * (S) -r θn (S)) 2 ] ≤ λM λm .

2 = 1 A 2 Theorem 2 (

 2122 Similarly, writing C = Eν [φ(S)φ(S) ], one can easily show thatEν [(r θ * (S) -r θn (S)) 2 ] = θ -θ * 2 C . Now, let x ∈ R d .For a symmetric and definite positive matrix M , write x M = √ x M x and SM its Cholesky decomposition. Write also M 2 the spectral norm of M . We havex C = x S C SC x = SC x SC S -1 A SAx 2 ≤ SC 2 S -SAx Eν [(r θ * (S) -r θn (S)) 2 ] = θ -θ * 2 C ≤ λM λm θ -θ * 2 A = λM λm (R(θn) -R(θ * )).Using the fact that R(θn) -R(θ * ) ≤ proves the result.B. PROPAGATION OF ERRORS: VALUES Propagation of errors: values)For a policy π, define Cπ as the smallest coefficient such that (1 -γ)ν (I -γPπ) -1 ≤ Cπν. Assume that θn satisfies R(θn) ≤ R(θ * ) + , then the associated value satisfiesE ν [(v r θ * * (S) -v

4 ( 1 2 ≤ 4 (= 2 (

 41242 t π r = (I -γPπ) -1 r, therefore:v r π * -v r π ≤ v r π * -v r π * + v r π -v r π = (I -γPπ * ) -1 (rr) + (I -γP π ) -1 (r -r) -γPπ * ) -1 + (I -γP π ) -1 ).Notice that Q is a stochastic kernel. We haveE ν [(v r π * (S) v r π (S)) 2 ] ≤ -γ) 2 E S∼ν [ E S ∼Q(.|S) [|r(S ) -r(S )|] 2 ] ≤ 4 (1 -γ) 2 E S∼ν [E S ∼Q(.|S) [(r(S ) -r(S )) 2 ]] = 4E S∼ν Q [(r(S) -r(S)) 2 ] (1 -γ) Cπ * + C π ) 2(1 -γ) 2 Eν [(r(S) -r(S)) 2 ] Cπ * + C π ) (1 -γ) 2 Eν [(r(S) -r(S)) 2 ],the second inequality being due to the Jensen's inequality and the third one to the definition of Cπ. Using Th. 1 allows concluding.

Table 1 : Performance of the handcrafted and SBIRL policies in the appointment scheduling simulations.

 1 ASR rejections Task Completion Handcrafted 8.266 ± 0.008 6.547 ± 0.036 0.73 ± 0.008 0.878 ± 0.006 SBIRL 500 8.373 ± 0.006 5.568 ± 0.035 0.541 ± 0.006 0.796 ± 0.001 SBIRL 1000 8.424 ± 0.006 4.938 ± 0.027 0.469 ± 0.005 0.791 ± 0.001 SBIRL 1500 8.446 ± 0.005 4.776 ± 0.025 0.44 ± 0.005 0.79 ± 0.001

	Policy	Score	#turns	#

  9 δ > max(0, 2.6-ln d). If assumptions 1 and-2 hold and if n ≥ 6ρ 2 d ln 9d δ , then with probability at least 1 -δ we have Eν [(r θ * (S) -r θn (S)) 2 ] ≤ Under the same assumptions, with probability at least 1 -δ: E ν [(v r θ * * (S) -v

						r θ * * π r θn	(S)) 2 ] ≤
					2 Cπ * + C π (1 -γ) 2 r θn *	λM λm	σ 2 (d + 2 d ln 9 δ + 2 ln 9 δ ) n	+ o	1 n	.
	λM λm	σ 2 (d + 2 d ln 9 δ + 2 ln 9 δ ) n	+ o	1 n	.


Hadrien Glaude

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

PAC learning of Probabilistic Automaton based on the Method of Moments

Probabilitic Finite Automata (PFA) are generative graphical models that define distributions with latent variables over finite sequences of symbols, a.k.a. stochastic languages. Traditionally, unsupervised learning of PFA is performed through algorithms that iteratively improves the likelihood like the Expectation-Maximization (EM) algorithm. Recently, learning algorithms based on the so-called Method of Moments (MoM) have been proposed as a much faster alternative that comes with PAC-style guarantees. However, these algorithms do not ensure the learnt automata to model a proper distribution, limiting their applicability and preventing them to serve as an initialization to iterative algorithms. In this paper, we propose a new MoM-based algorithm with PAC-style guarantees that learns automata defining proper distributions. We assess its performances on synthetic problems from the PAutomaC challenge and real datasets extracted from Wikipedia against previous MoM-based algorithms and EM algorithm.

, if we allow the number of samples to depend polynomially on some parameters measuring the complexity of the target 1 now with Google DeepMind

Introduction

In this paper, we address the problem of learning a distribution with latent variables over sequences of symbols from independent samples drawn from it. In particular, we are interested in distributions realized by generative graphical models called Probabilistic Finite Automata (PFA). Traditionally, algorithms to learn PFA rely on iterative procedures that maximize the joint likelihood like the gradient ascent or EM and its variants. However, these algorithms do not scale well with the number of samples and latent variables to the point where obtaining good solutions for Proceedings of the 33 rd International Conference on Machine Learning, New York, NY, USA, 2016. JMLR: W&CP volume 48. Copyright 2016 by the author(s). large models becomes intractable. In addition, they are prone to get stuck in local optima. There exist also full Bayesian methods like variational Bayes or Gibbs sampling, but these methods are computationally expensive and strongly rely on assumptions made on the prior.

A recent alternative line of work consists in modeling the distribution at sight by a Multiplicity Automaton (MA), also called weighted finite automaton [START_REF] Balle | Learning finite-state machines: statistical and algorithmic aspects[END_REF]. MA are graphical models that realize functions over finite sequences of symbols. Hence, they encompass a large variety of linear sequential systems in addition to stochastic languages [START_REF] Thon | Links between multiplicity automata, observable operator models and predictive state representations-a unified learning framework[END_REF]. For example, when MA model stochastic processes, they are equivalent to Observable Operator Models [START_REF] Thon | Links between multiplicity automata, observable operator models and predictive state representations-a unified learning framework[END_REF]. When considering action-observation pairs as symbols, MA can model controlled processes and are equivalent to Predictive State Representation [START_REF] Glaude | Subspace identification for predictive state representation by nuclear norm minimization[END_REF]. In fact, MA are strictly more general and infinitely more compact than PFA, Hidden Markov Models (HMMs) and Partially Observable Markov Decision Processes (POMDPs). Casting the learning problem into the one of learning a more general class of models allows using the MoM. This method leverages the fact that low order moments of distributions contain most of the distribution information and are typically easy to estimate. MoM-based algorithms have several pros over iterative methods. First, they are extremely fast. Their complexity is linear in the number of samples as estimated moments can be computed in one pass on the training set. The time complexity is also polynomial in the learnt model size as these algorithms rely only on few linear algebra operations to recover the parameters. In addition, MoM-based algorithms are often consistent with theoretical guarantees in the form of finite-sample bounds on the distribution. These parameters are typically small singular values of matrices defined by the target distribution.

However, current MoM-based algorithms have a major con. Although errors in the estimated parameters are bounded, they may correspond to automatons that lie outside the class of models defining proper distributions. Hence, the learnt models can output negative values or values that does not sum to one. In the PAC terminology, the learning is said to be improper. As mentioned in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF][START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF], this is a longstanding issue called the Negative Probability Problem (NPP). For some applications, the NPP is a major issue. For example, in reinforcement learning [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], the value iteration algorithm can diverge when planning with an unbounded measure. Although some heuristics, that perform a local normalization, exist to recover probabilities on a finite set of events, they prevent the theoretical guarantees to hold. In addition, NPP prevents the use of MoM-based algorithms to initialize a local search with an iterative algorithm like EM. Some attempts try to perform a global normalization by projecting the learnt model onto the space of valid model parameters [START_REF] Mossel | Learning nonsingular phylogenies and hidden markov models[END_REF][START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF][START_REF] Hsu | A spectral algorithm for learning hidden Markov models[END_REF][START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF]. While the resulting model is usable, these steps create an additional error and perform poorly in the experiments.

In this paper, we adopt the opposite approach. Instead of considering a more general class of automaton, we identify a subclass of models called Probabilistic Residual Finite Automaton (PRFA). We show that PRFA are PAClearnable using a MoM-based algorithm that returns PFA and thus avoids the NPP. Although PRFA are strictly less general than PFA, their expressiveness is large enough to closely approximate many distributions used in practice. In addition, learnt models can serve as a good initialization to iterative algorithms that perform a local search in the more general class of PFA. The paper is organized as follows: in Section 2, we recall the definition of a PFA and the basic Spectral learning algorithm; in Section 3 we define PRFA and a provable learning algorithm, CH-PRFA, that runs in polynomial time; finally, we assess the performance of CH-PRFA on synthetic problems and a real large dataset that cannot be handled by traditional methods like EM.

Background

Probabilistic Finite Automaton

PFA are graphical models constrained to represent distributions over sequences of symbols. Let Σ be a set of symbols, also called an alphabet. We denote by Σ , the set of all finite words made of symbols of Σ, including the empty word ε. Words of length k form the set Σ k . Let u and v ∈ Σ , uv is the concatenation of the two words and uΣ is the set of finite words starting by u. We are interested in capturing a distribution over Σ . Let p be such a distribution, for a set of words S, we define p(S) = u∈S p(u), in particular we have that p(Σ ) = 1. In addition, for any word u we define p such that p(u) = p(uΣ ). Thus, p defines distributions over prefixes of fix length : ∀n, u∈Σ n p(u) = 1. Some of these distributions can be modeled by graphical models called Probabilistic Finite Automaton (PFA).

Definition. A PFA is a tuple Σ, Q, {A o } o∈Σ , α 0 , α ∞ ,
where Σ is an alphabet and Q is a finite set of states. Matrices A o ∈ IR + |Q|×|Q| contain the transition weights. The vectors α ∞ ∈ IR + |Q| and α 0 ∈ IR + |Q| contain respectively the terminal and initial weights. These weights should verify,

1 α 0 = 1 α ∞ + o∈Σ A o 1 = 1 (1)
A PFA realizes a distribution over Σ , [START_REF] Denis | On rational stochastic languages[END_REF], defined by

p(u) = p(o 1 . . . o k ) = α 0 A u α 0 = α 0 A o1 . . . A o k α ∞ .
(2) Because of the constraints defined in Equation ( 1), the weights belong to [0, 1] and can be viewed as probabilities over initial states, terminal states and transitions with symbol emission. For a word, we define a path as a sequence of states starting in an initial state, transiting from state to state, emitting symbols of the word in each state and exiting in a final state. The probability of a path is the product of the weights along the path including initial and final weights. Hence, the probability of a word is given by the sum of all paths probabilities, as written in Equation (2). A path (resp. word) with a positive probability is called an accepting path (resp. word).

Spectral Learning

Actually, PFA define a particular kind of Multiplicity Automaton (MA). A MA is also a tuple Σ, Q, {A o } o∈Σ , α 0 , α ∞ but without any constraints on the weights, which can be negative and thus lose their probabilistic meaning. Thus, the function realized by a MA is not constrained to be a distribution. In the sequel, we call a Stochastic MA (SMA) a MA that realizes a distribution. In [START_REF] Denis | On rational stochastic languages[END_REF], the authors showed that SMA are strictly more general than PFA and can be infinitely more compact [START_REF] Esposito | Contribution à l'inférence d'automates probabilistes[END_REF].

The Spectral algorithm presented in this section relies on the Hankel matrix representation of a function to learn a MA. Let f : Σ → IR be a function, we define H ∈ IR Σ ×Σ the bi-infinite Hankel matrix whose rows and columns are indexed by Σ such that H[u, v] = f (uv).

Hence, when f is a distribution, H contains occurrence probabilities that can be estimated from samples by counting occurrences of sequences. Let for all o ∈ Σ, H o ∈ IR Σ ×Σ and h ∈ IR Σ be such that H o (u, v) = f (uov), h(u) = f (u). These vectors and matrices can be extracted from H. The Hankel representation lies at the heart of all MoM-based learning algorithms, because of the following fundamental theorem.

Theorem 1 (See [START_REF] Carlyle | Realizations by stochastic finite automata[END_REF]). Let f be a function realized by a MA with n states, then rank(H) ≤ n. Conversely, if the Hankel matrix H of a function f : Σ → IR has a finite rank n, then f can be realized by a MA with exactly n states but not less.

For a MA with n states, observe that H[u, v] = (α 0 A u )(A v α ∞ ). Let P ∈ K Σ ×n and S ∈ K n×Σ be matrices defined as follows,

P = ((α 0 A u ) ) u∈Σ , S = (A v α ∞ ) v∈Σ ,
then H = P S. Moreover, we have that,

H o = P A o S, h = α 0 S, h = P α ∞ . (3)
So, the MA parameters can be recovered by solving Equation (3). Hopefully, we do not need to consider the biinfinite Hankel matrix to recover the underlying MA. Given a basis B = (P, S) of prefixes and suffixes, we denote by H B the sub-block of H. Similarly, H o B is a sub-block of H o and h P and h S are sub-blocks of h. We say, that a basis B is complete if H B has the same rank than H. In [START_REF] Balle | Learning finite-state machines: statistical and algorithmic aspects[END_REF], the author shows that if B = (P, S) is a complete basis, by also defining P over P, S over S, we can recover a MA using Equation (3). Several methods are proposed in the literature to build a complete basis from data. In the experiments, we used the most frequent prefixes and suffixes.

Once a basis is chosen, the Spectral algorithm first estimates ĤB . Then, it approximates ĤB with a low dimensional factorized form ĤB ≈ Û D V through a truncated Singular Value Decomposition (SVD). Finally, setting P = Û D and Ŝ = V , the algorithm solves Equation (3), through linear regression. Because of the properties of the SVD, this leads to the following equations :

Âo = D-1 Û Ĥo B V , α 0 = ĥ S V , α∞ = D-1 Û ĥP .
Although the Spectral algorithm can return a MA arbitrary close to a SMA that realizes the target distribution, it does not ensure that the returned MA will be a SMA (and so a PFA). This causes the NPP explained in introduction.

Recalling that our goal is to learn proper distributions, one would like to add constraints to ensure that the MA learned by regression realizes a proper distribution. Unfortunately, this would require two things : 1) the non-negativity of series for any word (∀ ∈ Σ , p(u) ≥ 0), 2) the convergence of the series to one ( u∈Σ p(u) = 1). Although 2) can be checked in polynomial time, 1) requires adding an infinite set of constraints during the linear regression step. That is why, in general, verifying if a MA is a SMA is undecidable [START_REF] Denis | On rational stochastic languages[END_REF]. Actually, it has been shown in [START_REF] Esposito | Contribution à l'inférence d'automates probabilistes[END_REF] that no algorithm can learn a SMA in the limit of an infinite number of samples with probability 1.

So, in a first attempt we could restrict ourselves to the learning of PFA that are identifiable in the limit with probability 1 [START_REF] Denis | Learning classes of probabilistic automata[END_REF]. However, in [START_REF] Kearns | On the learnability of discrete distributions[END_REF], the authors showed that PFA are not PAClearnable by reduction to the learning of noisy parity functions, which is supposed to be difficult. Note that from [START_REF] Abe | On the computational complexity of approximating distributions by probabilistic automata[END_REF], we know that this negative result comes from the computational complexity, as only a polynomial number of samples could suffice. Thus, in this work we will focus on a smaller, but still rich, set of automata, called Probabilistic Residual Finite Automata (PRFA) that have been introduced in [START_REF] Denis | On rational stochastic languages[END_REF].

Probabilistic Residual Finite Automata

This section defines a particular kind of MA called Probabilistic Residual Finite Automata that realizes distributions. First, for any word u, we define the linear operator u on functions of IR Σ such that ∀v ∈ Σ , up(v) = p(uv). Then, for any distribution p, we denote, for each word u such that p(u) > 0, p u the conditional distribution defined by p u = up p(u) . In addition, for a PFA Σ, Q, {A o } o∈Σ , α 0 , α ∞ realizing a distribution p, we denote, for all q ∈ Q, p q the distribution defined by p q (u) = 1 q A u α ∞ . Thus, p q (u) is the probability of observing u starting from the state q. Similarly, p u (v) is the probability to observe v after u.

Definition. A PRFA is a PFA (Σ, Q, α 0 , A, α ∞ ) such that for all state q ∈ Q, there exists a word u ∈ Σ such that p(u) > 0 and p q = p u .

In particular, a PFA such that, for all state q there exists at least one prefix of an accepted word that ends only in state q, is a PRFA. In addition, if the PFA is reduced (there is no PFA with strictly less state realizing the same language), the converse is true. Note that as a PRFA is a PFA, it realizes a distribution, which we denote by p. In addition, for all q ∈ Q, p q is also a distribution (see [START_REF] Denis | On rational stochastic languages[END_REF]).

From that definition, we see that PRFA are more general than PFA with deterministic transition (PDFA). In fact, PRFA are strictly more general than PDFA but strictly less general than PFA [START_REF] Esposito | Contribution à l'inférence d'automates probabilistes[END_REF]. Similarly, the equiva-lent of PRFA for stochastic processes lies between HMMs and any finite order Markov chains. Now, we show how the existence of a set of words verifying some properties characterizes a PRFA. This equivalence allows us to design a learning algorithm. Proposition 2. Let p be a distribution, if there exists a set of words R and two associated sets of non-negative reals {a v u,o } u,v∈R,o∈Σ and {a v ε } v∈R , such that ∀u ∈ R, p(u) > 0 and,

∀u ∈ R, o ∈ Σ, ȯp u = v∈R a v u,o p v and p = v∈R a v ε p v , (4) then, Σ, Q, {A o } o∈Σ , α 0 , α ∞ defines a PRFA realizing p, where Q = R, α 0 = (a u ε ) u∈R , (5) α ∞ = (p u (ε)) u∈R , (6) ∀u, v ∈ R, A o [u, v] = a v u,o . (7)
Proof. First, we show by induction on the length of u that ∀u ∈ Σ , (p v (u)) v∈R = A u α ∞ . By the definition of α ∞ , the property is verified for u = ε. Assume the property is true for u ∈ Σ ≤n , then for any o ∈ Σ, let v = ou and we have,

A v α ∞ = A o A u α ∞ = A o (p w (u)) w∈R = w∈R a w w ,o p w (u)
w ∈R

(by Equation ( 5))

= ( ȯp w (u)) w ∈R (by Equation ( 4)) = (p w (v)) w∈R . Now, for all u ∈ Σ , with Equations ( 4) and ( 5), we have

α ∞ A u α ∞ = w∈R a w ε p w (u) = p(u).
It remains to show that Σ, Q, {A o } o∈Σ , α 0 , α ∞ defines a PRFA, meaning that it satisfies Equation (1) and

∀q ∈ Q, ∃u ∈ Σ , (p(u) > 0) ∧ (p q = p u ). (8)
From Equation (4), we have p = u∈R a u ε p u . In particular p(Σ ) = u∈R a u ε p u (Σ ) and p(Σ ) = p u (Σ ) = 1 for all u ∈ R implies that u∈R a u ε = 1. Finally, from Equation (4), we have for all u ∈ R that

ȯp u (Σ ) = v∈R a v u,o p v (Σ ), o∈Σ p u (oΣ ) = o∈Σ v∈R a v u,o , p u (ΣΣ ) = o∈Σ v∈R a v u,o .
As p u is a distribution, p u (Σ ) = p u (ε)+p u (ΣΣ ). So, we obtain that α ∞ + A Σ 1 = 1. Hence, (α 0 , A, α ∞ ) satisfies Equation (1). We showed that for any word v ∈ Σ ,

(p q (v)) q∈Q = A v α ∞ = (p u (v)) u∈R .
So by taking Q = R, we have that Equation ( 8) is satisfied and

Σ, Q, {A o } o∈Σ , α 0 , α ∞ is a PRFA.
In Proposition 3, we show the converse, i.e. the existence of the R, {a v u,o } u,v∈R,o∈Σ and {a v ε } v∈R for a PRFA. In fact, the existence of coefficients satisfying Equation ( 4) can be reformulated using the notion of conical hull. We denote by coni(E) the set defined by,

coni(E) = { e∈E α e e|α e ∈ IR + }. Proposition 3. Let Σ, Q, {A o } o∈Σ , α 0 , α ∞ be a PRFA
and p the distribution it realizes, then there exists a set of words R such that

∀u ∈ R, p(u) > 0, p ∈ coni{p u |u ∈ R}, ∀u ∈ R, o ∈ Σ, ȯp u ∈ coni{p v |v ∈ R}.
Proof. From the definition of a PRFA given Equation ( 8), there exists a finite set of words R of size bounded by |Q|, the number of state, such that ∀q ∈ Q, ∃u ∈ R, (p(u) > 0) ∧ (p q = p u ).

To prove the existence of {a v u,o } u,v∈R,o∈Σ and {a v ε } v∈R , we write that for any w ∈ Σ ,

p(w) = α 0 A w α ∞ = α 0 (p q (w)) q∈Q = α 0 (p u (w)) u∈R .
As for a PRFA, α 0 is a vector of non-negative coefficients, we have that p ∈ coni{p u |u ∈ R}. Similarly, for all u ∈ R, o ∈ Σ, we have for any

w ∈ Σ ȯp(w) = α 0 A ow α ∞ = α 0 A o (p q (w)) q∈Q = α 0 A o (p u (w)) u∈R .
As α 0 and A o are non-negative, α 0 A o is a vector with non-negative coefficients and we have that ∀u ∈ R, o ∈ Σ, ȯp u ∈ coni{p v |v ∈ R}.

Learning PRFA

As in the Spectral algorithm, our method assumes that a complete basis B = (P, S) is provided. To simplify the discussion, we will also assume that the empty word ε is included as a prefix and a suffix in the basis. For convenience, we denote 1 ε a vector on P or S, depending on the context, filled with zeros but a one at the index of ε. In addition, this basis must contain the prefixes and the suffixes allowing the identification of a set R generating a conical hull containing p and, for all u ∈ R, ȯp u . In the sequel, we denote p (resp. p u , ȯp u ), the vector representation of p (resp. p u , ȯp u ) on the basis of suffixes S. Thus, in addition to be complete, the basis B must be residual.

Definition (Residual basis).

A basis B = (P, S) is residual if the conical hull coni{p u |u ∈ Σ , p(u) > 0} projected on S coincide with coni{p u |u ∈ P, p u > 0}.

In addition, we assume the Algorithm 1 is provided with the minimal dimension d of PRFA realizing p. Hence, as the basis is complete and residual, we have that the hypothesis in Proposition 2 are satisfied, i.e. there exists

R ⊂ P such that ∀u ∈ R, p(u) > 0, p ∈ coni{p u |u ∈ R} and ∀u ∈ R, o ∈ Σ, ȯp u ∈ coni{p v |v ∈ R}.
The CH-PRFA algorithm works by first estimating the p u with u ∈ P. Then, it finds the set R. We will see this step can be solved using near-separable Non-negative Matrix Factorization (NMF). To identify R, instead of using {p u |u ∈ P}, we used { du |u ∈ P}, where du is defined in Algorithm 1 because it improves the robustness and helps in the proof of Theorem 4. Finally, the parameters of a PFA are retrieved through linear regressions. Because of estimation errors, we need to add non-negativity and linear constraints to ensure the parameters defines a PFA. Hence, CH-PRFA returns a PFA but not necessarily a PRFA. In the PAC terminology, our algorithm is improper like the Spectral one. However, it is inconsequential as a PFA realizes a proper distribution. In contrast, we recall that the Spectral learning algorithm returns a MA that does not necessarily realizes a proper distribution. In addition, our algorithm is proper in limit whereas the Spectral algorithm is not.

In the literature, many algorithms for NMF have been proposed. Although in its general form NMF is NP-Hard and ill-posed [START_REF] Gillis | The why and how of nonnegative matrix factorization[END_REF], in the near-separable case the solution is unique and can be found in O(k |S| |P|) steps. Stateof-the-art algorithms for near-separable NMF comes with convergence guarantees and robustness analysis. In our experiment, the Successive Projection Algorithm (SPA), analyzed in [START_REF] Gillis | Fast and robust recursive algorithmsfor separable nonnegative matrix factorization[END_REF], gave good results. In addition, SPA is very efficient and can be easily distributed.

In Algorithm 1, the minimization problems can be cast as quadratic optimization problems under linear constraints with convex costs. This kind of problem can be solved in polynomial time using a solver like [START_REF] Mosek | The MOSEK Python optimizer API manual Version 7.1[END_REF]. As the cost is convex, solvers converge to a stationary point. In addition, all stationary points have the same cost and are optimal [START_REF] Lőtstedt | Perturbation bounds for the linear least squares problem subject to linear inequality constraints[END_REF]. To get a unique solution, it is possible to look for the minimal norm solution, as a Moore-Pseudo inverse does.

u ∈ R do {â v u,o } ← argmin {a v u,o } o∈Σ ȯp u - v∈ R a v u,o pv 2 s.t. v∈ R,o∈Σ a v u,o ȯ = 1 -pu 1 and a v u,o ≥ 0 .
end for

{â u ε } ← argmin {a u ε } p - u∈ R a u ε pu 2 s.t. u∈ R a u ε = 1 and a u ε ≥ 0. α 0 ← (â u ε ) u∈ R, α∞ ← (p u 1 ) u∈ R. for o ∈ Σ do Âo ← âv u,o u,v∈ R.

end for

Using the perturbation analysis of SPA [START_REF] Gillis | Fast and robust recursive algorithmsfor separable nonnegative matrix factorization[END_REF] and the one of quadratic optimization [START_REF] Lőtstedt | Perturbation bounds for the linear least squares problem subject to linear inequality constraints[END_REF], we show the following non-asymptotic bound. . For all 0 < δ < 1, there exists a constant K such that, for all t > 0, > 0, with probability 1 -δ, if

n ≥ K t 4 d 4 |Σ| 2 σ 10 d log |P| δ ,
CH-PRFA returns a PFA realizing a proper distribution p such that

u∈Σ ≤t |p(u) -p(u)| ≤ .
Proof. In the supplementary material. Now, we compare this result to previous bounds on the Spectral algorithm. First our bound depends on n, instead of N = |D|. Using Hoeffding inequality, we could obtain a bound on N depending on the inverse of the probability of the least frequent prefix of P. This dependence comes from the use of conditional distributions (the p u ) instead of joint distributions. Indeed, the Hankel matrix used in Spectral contains joint distributions. Removing this dependency on the prefix set seems possible by changing how the NMF is computed. This direction will be explored in future researches. Then, we address the convergence speed. First, the dependency on log(|P|) could be removed using recent dimension-free concentration bounds on Hankel matrices [START_REF] Denis | Dimension-free concentration bounds on Hankel matrices for spectral learning[END_REF]. Secondly, although the term |Σ| is better than |Σ| 2 in the classical error bounds for Spectral [START_REF] Hsu | A spectral algorithm for learning hidden Markov models[END_REF][START_REF] Balle | Learning finite-state machines: statistical and algorithmic aspects[END_REF], the results from [START_REF] Foster | Spectral dimensionality reduction for HMMs[END_REF] on HMMs suggest that the error bound could be independent of |Σ|. Third, the convergence speed in O( -2 log(δ -1 )) comes directly from concentration results and seems optimal. Finally, the required number of sample depends on d 4 which is worst than in the bounds for the Spectral algorithm. This strong dependency on d comes from the constraints in the optimization problems. Finally, the high polynomial order of σ -1 d in the bound is a direct consequence of SPA.

Experimental results

Near-Separable NMF

In this Section, we give more details on the identification of a conical hull containing a set of vectors. This problem has often been addressed as a NMF problem, where a nonnegative matrix as to be decomposed in two non-negative matrices of reduced dimensions. One of these contains the vectors supporting the conical hull and the other the conical combinations to recover all the vectors of the original matrix. This problem, in its general form, is ill-posed and NP-Hard [START_REF] Gillis | The why and how of nonnegative matrix factorization[END_REF] and algorithms often rely on alternated optimization that converge only to a local optimum. A decade ago [START_REF] Donoho | When does nonnegative matrix factorization give a correct decomposition into parts?[END_REF], a sufficient condition, called separability, have been identified to ensure the uniqueness of the solution. Geometrically, separability implies that the vectors supporting the conical hull are contained in the original matrix. Since, many algorithms relying on different additional assumptions have been proposed to solve NMF for separable matrices.

In particular, the Successive Projection Algorithm used both in the finite sample analysis and in the experiments, identifies recursively the supporting vectors among the column of the original matrix. It assumes that the matrix formed with the supporting vectors, (p u (v)) u∈R in our case, has full rank. Although the full rank assumption is not generally true, as we work with empirical estimate it is satisfied with probability 1. Moreover, experimental results does not seem to suffer from that assumption. Additionally, SPA assumes the supporting vectors form a convex hull instead of a conical hull. This assumption can be made without loss of generality as columns from the original matrix can be normalized to belong to the simplex. In our case, this assumption is already satisfied because of the constraints given in Equation ( 1). However, we will explain why a direct application of the SPA algorithm constrains which prefixes can be included in the basis.

As suggested by the finite sample analysis and experiences, incorporating prefixes with low occurring probabilities in the basis degrades strongly the performances. In fact, the robustness of SPA depends on the maximum error in norm made on p u . However, probabilities conditioned on rare prefixes are not well estimated. Therefore, taking only the most frequent prefixes is needed to achieve good results. In order to fix this issue, we propose another way to use the SPA algorithm. By rescaling the vectors, we can obtain the same amount of uncertainty for all prefixes. So, instead of working with estimates of p u = p(uv) p(u) v∈S , we use the vectors q u = (p(uv)) v∈S (where u ∈ P). In other words, the SPA algorithm is executed on the same Hankel matrix than the Spectral algorithm. The Figure 1 shows the differences between the conical hull and the convex hull as well as the effect of normalizing the vectors to the simplex.

One can wonder if the algorithm is still valid because, even with the true vectors, the result of SPA will differ depending on whether we use p u or q u . Let R be the set of prefixes identified by SPA using p u . Let Q be the set of prefixes identified by SPA using q u . The algorithm is still valid if only a finite number of vectors q u lies outside the convex hull generated by {q u |u ∈ R}. Thus, as long as the model dimension is large enough Q will contain R. Including more prefixes than needed is inconsequential for the remaining part of the algorithm. So, using q u instead of p u leave the validity CH-PRFA unchanged. In the experiments, we used this variation of SPA as it increased the performances without increasing models sizes.

PAutomaC Challenge

The Probabilistic Automata learning Competition (PAu-tomaC) deals with the problem of learning probabilistic distributions from strings drawn from finite-state automata.

From the 48 problems available, we have selected the same 12 problems than in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF], to provide a fair comparison with other algorithms. The generating model can be of three kinds: PFA, HMMs or PDFA. A detailed description of each problem can be found in [START_REF] Verwer | Results of the pautomac probabilistic automaton learning competition[END_REF]. We compared CH-PRFA and CH-PRFA+BW (BW initialized with CH-PRFA) to Baum-Welch (BW) with 3 random restarts of 3 iterations (then the best run is continued for a maximum of 30 iterations) and other MoMbased algorithms : CO [START_REF] Balle | Local loss optimization in operator models: A new insight into spectral learning[END_REF], Tensor [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF], NNSpectral [START_REF] Glaude | Non-negative spectral learning for linear sequential systems[END_REF] and Spectral with variance normalization [START_REF] Cohen | Experiments with spectral learning of latent-variable PCfgs[END_REF].

In the experiments, negative values outputted by the algorithms are set to zero, then we normalize to obtain probabilities. These MoM-based algorithms have been trained using statistics on sequences and subsequences as proposed in [START_REF] Balle | Learning finite-state machines: statistical and algorithmic aspects[END_REF]. The best result is then selected. For all MoM-based algorithm, we used basis sizes varying between 50 and 10000 (except for CO where the computation time limits the basis size to 200 and NNSpectral to 500).

For BW, we stopped the iterations if after 4 iteration the score was not improving. To obtain probability-like values from MoM-based algorithms we used several tricks. For Spectral and CO, we zeroed negative values and normalized. For NNSpectral, we just normalized. For Tensor, we projected the transition and observation matrix onto the simplex, as described in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF]. Finally, we assessed the quality of the learned distribution p by the perplexity. It corresponds to the average number of bits needed to represent a word using the optimal code p .

Perplexity(M) = 2 -u∈T p (u) log(p M (u)) .

We also measured the quality of p by computing the Word Error Rate (WER) which is the average number of incorrect predictions of the next symbols given the past ones. A grid search was performed to find the optimal dimension and basis size for each of the performance metric. On Tables 1 and2, we ranked the algorithm according to their average performances on the twelve problems. For the WER, the average corresponds to the mean. For the perplexity, the average corresponds to 2 -12 i=1 u∈T p i (u) log(p M i (u)) .

Wikipedia

We also evaluated CH-PRFA and CH-PRFA+BW on raw text extracted from English Wikipedia pages. The training set is made of chunks of sequences of 250 characters randomly extracted from the 2GB corpus used in [START_REF] Sutskever | Generating text with recurrent neural networks[END_REF]. Each character stands for a symbol. We restricted ourself to 85 different symbols. For the training phase, we used a small set of 500 chunks and a medium one of 50000 characters. For testing, an independent set of 5000 chunks has been used. The algorithms are the same than the ones evaluated for PAutomaC but slightly modified to fit a stochastic process. For example, we used the suffix-history algorithm [START_REF] Wolfe | Learning predictive state representations in dynamical systems without reset[END_REF] to estimate the occurrence probabilities in the Hankel matrices. For CH-PRFA, we changed the constraints defined in Equation ( 1) and used during the optimization problems to be the ones of a stochastic process. As the results of CO and Tensor were very poor (between 0.05 and 0.10 for the likelihood), we did not report them for clarity. The time taken by BW were excessive (more than a day in comparison to the tens of minutes needed by the others) for model sizes above 60 on the small training set and for all dimensions on the medium one. Performance is measured by the average likelihood of the next observation given the past. More precisely, on chunk of 250 characters, the 50 firsts served to initialize the belief, then the next observations are predicted based on the characters observed so far. The likelihood of all these predictions is then averaged over observations and chunks. We also computed the number of bits per character (BPC) by averaging log(P (o t+1 |o 1:t )) over sequences o 1:t of symbols.

Discussion

On PAutomaC, the top performer is NNSpectral both for the perplexity and the WER. CH-PRFA performed almost as good as NNSpectral considering the perplexity.

The mean WER of CH-PRFA is slightly less good than NNSpectral, Spectral and BW. Although CH-PRFA is not the top performer, it has PAC-style guarantees in contrast to NNSpectral and BW and is much faster as shown on Figure 4. Finally, in contrast to Spectral, CH-PRFA does not suffer from the NPP.

On Wikipedia, for the likelihood, the top performer is Spectral. We believe that its good scores come from the fact that MA are more strictly more general and more compact than PFA and so PRFA. On a training set like Wikipedia, it can be a big asset as a natural language model is likely to be very complex. The mixed performances of CH-PRFA can be explained by the lack of expressiveness of PRFA. The lack of expressiveness of PRFA is then filled using BW to improve the model quality as BW learns a PFA. Hence, when trained on the small set, CH-PRFA+BW achieves the best performances and beats Spectral. However, on the medium set the BW algorithm takes too much time to be a decent alternative. For the BPC, results are quite surprising as Spectral is the least performer. In fact, the BPC gives more importance to rare events than the conditional likelihood. So, Spectral predicts better frequent events than rare events. Indeed, small probabilities associated to rare events are likely to be a sum of products of small parameters of the MA. As parameters are not constrained to be non-negative, a small error can flip the sign of small parameters which in turn leads to large errors. That is why, NNSpectral and CH-PRFA performed better than Spectral for the BPC.

Finally, using two sets of different sizes shows that BW do not scale well. The running time of BW between the two sets has changed from few minutes to at least a day, whereas the one of CH-PRFA has only increased by few seconds.

Conclusions

In this paper, we proposed a new algorithm based on a nearseparable NMF and constraint quadratic optimization that can learn in polynomial time a PRFA from the distribution p it realizes. In addition, even if p is not realized by a PRFA or is empirically estimated, our algorithm returns a PFA that realizes a proper distribution closed to the true distribution. We established PAC-style bounds that allowed us to tweak the NMF to achieve better results. Then, we empirically demonstrated its good performances in comparison to other MoM-based algorithms, and, its scalability in comparison to BW. Finally, experiments have shown that initializing BW with CH-PRFA can substantially improve the performances of BW. For future works, extending Algorithm 1 to handle controlled processes would allow designing consistent reinforcement learning algorithms for non-Markovian environments.

Theorem 4 .

 4 Let p be a distribution realized by a minimal PRFA of size d, B = (P, S) be a complete and residual basis, we denote by σ d the d-th largest singular values of (p u (v)) u∈R . Let D be a training set of words generated by p, we denote by n the number of time the least occurring prefix of P appears in D (n = min u∈P |{∃v ∈ Σ |uv ∈ D}|)

Figure 1 .

 1 Figure 1. On the left, a conical hull of a set of points and their projection onto the simplex. On the right, a convex hull of a set of points.

Figure 2 .Figure 3 .

 23 Figure2. Likelihood on Wikipedia (the higher is the better). Asterisks denote scores on the medium training set.

Figure 4 .

 4 Figure 4. Mean learning time on twelve problems of PAutomaC.

  Algorithm 1 CH-PRFA Input: A targeted dimension d, a separable complete basis B = (P, S) and a training set. Estimate pu and ȯp u from the training set. du ← p u ȯ1 p u . . . ȯ|Σ| p u end for Find R a subset of d prefixes of P such that ∀u ∈ R, du > 0 and d ∈ coni{ du |u ∈ R}. for

	Output: A PFA Σ, R, Âo

o∈Σ

, α0 , α∞ .

for u ∈ P do

Table 1 .

 1 Average WER on twelve problems of PAutomaC.

	RANK ALGORITHM	WER
	1. NNSPECTRAL	64.618
	2. CH-PRFA+BW 65.160
	3. SPECTRAL	65.529
	4. BW	66.482
	5. CH-PRFA	67.141
	6. CO	71.467
	7. TENSOR+BW	73.544
	8. TENSOR	77.433

Table 2 .

 2 Average perplexity on twelve problems of PAutomaC.

	RANK ALGORITHM	PERPLEXITY
	1. NNSPECTRAL	30.364
	2. CH-PRFA+BW	30.468
	3. CH-PRFA	30.603
	4. CO	35.641
	5. BW	35.886
	6. SPECTRAL	40.210
	7. TENSOR+BW	47.655
	8. TENSOR	54.000