Mohammad Gheshlaghi Azar
email: m.azar@science.ru.nl

Rémi Munos
email: remi.munos@inria.fr

Mohammad Ghavamzadeh
email: mohammad.ghavamzadeh@inria.fr

Hilbert J Kappen
email: b.kappen@science.ru.nl

Reinforcement Learning with a Near Optimal Rate of Convergence

We consider the problem of model-free reinforcement learning (RL) in the Markovian decision processes (MDP) under the PAC ("probably approximately correct") model. We introduce a new variant of Q-learning, called speedy Q-learning (SQL), to address the problem of the slow convergence in the standard Q-learning algorithm, and prove PAC bounds on the performance of SQL. The bounds show that for any MDP with n state-action pairs and the discount factor γ ∈ [0, 1) a total of O n log(n)/((1 -γ) 4 ǫ 2 ) step suffices for the SQL algorithm to converge to an ǫ-optimal action-value function with high probability. We also establish a lower-bound of Ω n log(n)/((1 -γ) 2 ǫ 2 ) for all RL algorithms, which matches the upper bound in terms of n and ǫ. Further, our results have better dependencies on ǫ, 1 -γ and the same dependency on n , and thus, are tighter than the best available results for Q-learning. SQL also improves on existing results for the batch Q-value iteration, so far considered to be more efficient than the incremental methods like Q-learning.

Introduction

The Markovian decision process (MDP) problem is a classical problem in the fields of operations research and decision theory. When an explicit model of the MDP, transition probabilities and reward function, is known, one can rely on dynamic programming (DP) [START_REF] Bellman | Dynamic Programming[END_REF] algorithms such as value iteration or policy iteration (see, e.g., Bertsekas, 2007a;[START_REF] Puterman | Markov Decision Processes, Discrete Stochastic Dynamic Programming[END_REF], to compute the optimal policy. For example, value iteration computes the optimal value function Q * by successive iterations of the Bellman operator T, which we define later in Section 2. One can show that, in the discounted infinite-horizon setting, the convergence of value iteration is exponentially fast since the Bellman operator T is a contraction mapping [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] on the action-value function Q. However, DP relies on an explicit knowledge of the MDP. In many real world problems the transition probabilities are not initially known, but one may observe transition samples using Monte-Carlo sampling, either as a single trajectory (rollout) obtained by following an exploration policy, or by simulating independent transition samples anywhere in the space by resorting to an oracle which draws samples anywhere in the state-action space. The field of reinforcement learning (RL) is concerned with the problem of approximating the optimal policy, or the optimal value function, from the observed reward and transition samples [START_REF] Cs | Algorithms for Reinforcement Learning[END_REF][START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF].

One may characterize RL methods as model-based or model-free. In model-based RL we first learn a model of the MDP and then we use this model for computing an approximation of the value functions by dynamic programming techniques. Model-free methods, in contrast, compute directly an approximation of the value function by making use of a sample-based estimate of the Bellman operator, without resorting to learning a model. Q-learning (QL) is a well-known model-free reinforcement learning (RL) algorithm that, incrementally, finds an estimate of the optimal action-value function (Watkins, 1989). The QL algorithm can be seen as a combination of the value iteration algorithm and stochastic approximation, where at each time step k a new estimate of the optimal action-value function is calculated by the following update rule:

Q k+1 = (1 -α k )Q k + α k (TQ k -ǫ k ) ,
where ǫ k and α k denote the estimation error and the learning step at time step k, respectively. One can show, using an induction argument, that for the choice of linear learning step, i.e., α k = 1 k+1 , QL is taking average over the estimates of Bellman operator throughout the learning process:

Q k+1 = 1 k + 1 k j=0 (TQ j -ǫ j ).
It is not then difficult to prove, using a law of large number argument, that the term 1/(k + 1) k j=0 ǫ j is asymptotically averaged out and, therefore, for k → ∞ the update rule of Q-learning becomes equivalent to Q k+1 = 1/(k + 1) k j=0 TQ j . The problem with this result is that the rate of convergence of the recursion Q k+1 = 1/(k + 1) k j=0 TQ j to Q * is significantly slower than the original Bellman recursion Q k+1 = TQ k . In fact, one can show that the asymptotic rate of convergence of Q-learning, with linear learning step, is of order Õ(1/k 1-γ ) (Szepesvári, 1997), which, in the case of γ close to 1, makes the convergence of Q-learning extremely slower than the standard value iteration, which enjoys fast convergence rate of Õ(γ k ). This slow rate of convergence, i.e., high sample complexity, may explain why the practitioners often prefer the batch RL methods such as approximate value iteration (AVI) [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] to Q-learning despite the fact that Q-learning has an equal, or in some cases even better, computational complexity (per sample) and memory requirements than the batch RL methods.

In this paper, we focus on RL problems that are formulated as finite state-action discounted infinite-horizon MDPs, and propose a new algorithm, called speedy Q-learning (SQL), that addresses the problem of slow convergence of Q-learning. At each time step k, SQL uses two successive estimates of the bellman operator TQ k and TQ k-1 to update the action value function Q k :

Q k+1 = α k Q k + (1 -α k ) [kTQ k -(k -1)TQ k-1 -ǫ k ] , (1)
that makes its space complexity twice as QL. However, this allows the SQL algorithm to achieve a significantly faster rate of convergence than QL, since it reduces the dependency on the previous Bellman operators from the average 1/(k + 1) k j=0 TQ j , in the case of QL, to only TQ k + O(1/(k + 1)):

Q k+1 = α k Q k + (1 -α k ) [kTQ k -(k -1)TQ k-1 -ǫ k ] = 1 k + 1 k j=0 (jTQ j -(j -1)TQ j-1 -ǫ j ) By induction = 1 k + 1 (kTQ k -(k -1)TQ k-1 + (k -1)TQ k-1 -(k -2)TQ k-2 + • • • + TQ -1 ) - 1 k + 1 k j=0 ǫ j = TQ k + 1 k + 1 (TQ -1 -TQ k ) - 1 k + 1 k j=0 ǫ j ,
with the choice of α k = 1/(k + 1). In words, the iterates of SQL, like QL, are expressed in terms of the average estimation error and, therefore, the SQL update rule asymptotically averages out the estimation errors (simulation noise) as well. However, speedy Q-learning has the advantage that at each time step k the iterate Q k+1 closely follows (up to a factor of O(1/(k + 1)) the latest Bellman iterate TQ k instead of the average 1/(k + 1) k j=0 TQ j in the case of QL. Thus, unlike QL, it does not suffer from the slow convergence due to slow down in the value iteration process (For a detailed comparison of the convergence rate of QL and SQL see Section 3.3).

The idea of using previous estimates of the action-values has already been used in order to improve the performance of QL. A popular algorithm of this kind is Q(λ) [START_REF] Peng | Incremental multi-step Q-learning[END_REF]Watkins, 1989), which incorporates the concept of eligibility traces in QL, and has been empirically shown to have a better performance than QL, i.e., Q(0), for suitable values of λ. Another recent work in this direction is Double Q-learning (van Hasselt, 2010), which uses two estimators for the action-value function in order to alleviate the over-estimation of action-values in Q-learning. This over-estimation is caused by a positive bias introduced by using the maximum action value as an approximation for the maximum expected action value (van Hasselt, 2010).

The rest of the paper is organized as follows. After introducing the notations used in the paper in Section 2, we present our Speedy Q-learning algorithm in Section 3. We first describe the synchronous and asynchronous version of the algorithm in Section 3.1, then state our main theoretical result, i.e., high-probability bounds on the performance of SQL as well as a new lower bound for the sample complexity of reinforcement learning, in Section 3.2, and finally compare our bound with the previous results on QL and Q-value iteration in Section 3.3. In Section 4, we evaluate the performance of SQL, numerically, on different problem domains. Section 5 contains the detailed proof of the performance bound of the SQL algorithm and a general new lower bound for RL. Finally, we conclude the paper and discuss some future directions in Section 6.

Preliminaries

In this section, we introduce some concepts and definitions from the theory of Markov decision processes (MDPs) and stochastic processes that are used throughout the paper. We start by the definition of supremum norm (ℓ ∞ norm). For a real-valued function g : Y → R, where Y is a finite set, the supremum norm of g is defined as g max y∈Y |g(y)|.

We consider the standard reinforcement learning (RL) framework [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF][START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF] in which a learning agent interacts with a stochastic environment and this interaction is modeled as a discrete-time discounted MDP. A discounted MDP is a quintuple (X, A, P, R, γ), where X and A are the set of states and actions, P is the state transition distribution, R is the reward function, and γ ∈ (0, 1) is a discount factor. We also denote the horizon of MDP by β defined by β = 1/(1 -γ). We denote by P (•|x, a) and r(x, a) the probability distribution over the next state and the immediate reward of taking action a at state x, respectively.1 To keep the representation succinct, we use Z for the joint state-action space X × A.

Assumption A1 (MDP regularity)

We assume Z is a finite set with the cardinality n. We also assume that the immediate rewards r(x, a) are in the interval [0, 1].2 A policy kernel π(•|•) determines the distribution of the control action given the past observations. The policy is called stationary if the distribution of the control action just depends on the last state x. It is deterministic if it assigns a unique action to all states x ∈ X. The value and the action-value functions of a policy π, denoted respectively by V π : X → R and Q π : Z → R, are defined as the expected sum of discounted rewards that are encountered when the policy π is executed. Given a MDP, the goal is to find a policy that attains the best possible values, V * (x) sup π V π (x), ∀x ∈ X. Function V * is called the optimal value function. Similarly the optimal action-value function is defined as where the max operator M over action-value functions is defined as (MQ)(y) = max a∈A Q(y, a) , ∀y ∈ X.3 Further, We define the cover time of MDP under the policy π as follows:

Q * (x, a) = sup π Q π (x,
Definition 1 (Cover time of MDP) Let t be an integer. Define τ π (x, t), the cover time of the MDP under the (non-stationary) policy π, as the minimum number of steps required to visit all state-action pairs z ∈ Z starting from state x ∈ X at time-step 0 ≤ t. Also, the state-action space Z is covered by the policy π if all the state-action pairs are performed at least once under the policy π.

The following assumption which bounds the expected cover time of the MDP guarantees that, asymptotically, all the state-action pairs are visited infinitely many times under the policy π.

Assumption A2 (Boundedness of the the expected cover time) Let 0 < L < ∞ and t be a integer. Then, under the policy π for all x ∈ X and t > 0 assume that:

E (τ π (x, t)) ≤ L.

Speedy Q-Learning

In this section, we introduce a new RL method, called speedy Q-Learning (SQL), derive performance bounds for the asynchronous and synchronous variant of SQL, and compare these bounds with similar results on standard Q-learning.

Algorithms

In this subsection, we introduce two variants of SQL algorithms, the synchronous SQL and asynchronous SQL. In the asynchronous version, at each time step, the action-value of only one state-action pair, the current observed state-action, is updated, while the action-values of rest of the state-action pairs remain unchanged. For the convergence of this instance of the algorithm, it is required that all the states and actions are visited infinitely many times, which makes the analysis slightly more complicated. On the other hand, having access to an oracle which can generate samples anywhere in the state-action space, the algorithm may be also formulated in a synchronous fashion, in which we first generate a next state y ∼ P (•|x, a) for each state-action pair (x, a), and then update the action-values of all the state-action pairs using these samples. The pseudo-code of the synchronous and asynchronous SQL are shown in Algorithm 1 and 2, respectively. One can show that asynchronous SQL is reduced to Algorithm 1 when the cover time τ π (x, t) = n for all x ∈ X and t ≥ 0, in which case the action-values of all state-action pairs are updated in a row. In other words, the Algorithm 1 can be seen as a special case of Algorithm 2. Therefore, in the sequel we only describe the more general asynchronous SQL algorithm.

As it can be seen from the update rule of Algorithm 2, at each time step, Algorithm 2 keeps track of the action-value functions of the two most recent iterations Q k and Q k+1 , and its main update rule is of the following form for all (x, a) ∈ Z at time step t and the iteration round k:

Q k+1 (X t , A t ) = (1 -α k )Q k (X t , A t ) + α k kT k Q k (X t , A t ) -(k -1)T k Q k-1 (X t , A t ) , (2) where T k Q(X t , A t ) = 1/|Y k | y∈Y k [r(X t , A t )+γMQ(y)
] is the empirical Bellman optimality operator using the set of next state samples Y k , where Y k is a short-hand notation for Y k,t (x, a) = {y|y ∼ P (•|x, a)}, the set of all samples generated up to time step t in round k conditioned on the state-action pair (x, a) ∈ Z. At each time step t, Algorithm 2 works as follows: (i) It simulates the MDP for one-step, i.e., it draws the action A t ∈ A from the distribution π(•|x) and then make a transition to a new state y k . (ii) it updates the two sample estimates T k Q k-1 (X t , A t ) and T k Q k (X t , A t ) of the Bellman optimality operator (for state-action pair (X t , A t ) using the next state y k ) applied to the estimates Q k-1 and Q k of the action-value function at the previous and current round k -1 and k, (iii) it updates the action-value function of (X t , A t ), generates Q k+1 (X t , A t ), using the update rule of Eq. 2, (iv) it checks for the condition that if all (x, a) ∈ Z have been visited at least one time at iteration k. If the condition is satisfied then we move to next round k + 1, otherwise k remains unchanged. Finally, (v) we replace X t+1 with y k and repeat the whole process until t ≥ T . Moreover, we let α k decays linearly with the number of iterations k, i.e., α k = 1/(k + 1), in Algorithm 2. Also, we notice that the update rule

T k Q k (X t , A t ) := (1 -η N )T k Q k (X t , A t ) + η N (r(X t , A t ) + γMQ k (y k )
) is being used to make an unbiased estimate of TQ k in an incremental fashion.

α k := 1 k + 1 ; foreach (x, a) ∈ Z do // Update the action-values for all (x, a) ∈ Z Generate the next state sample y k ∼ P (•|x, a); T k Q k-1 (x, a) := r(x, a) + γMQ k-1 (y k ); T k Q k (x, a) := r(x, a) + γMQ k (y k ); Q k+1 (x, a) := (1 -α k )Q k (x, a) + α k kT k Q k (x, a) -(k -1)T k Q k-1 (x, a) ; // SQL update rule t := t + 1; end k := k + 1; until t ≥ T ; return Q k 3.

Main Theoretical Results

The main theoretical results of this paper are expressed as high-probability bounds over the performance of the SQL algorithms for both synchronous and asynchronous cases. We also report a new lower bound on the number of transition samples required, for every reinforcement learning algorithm, to achieve an ǫ-optimal estimate of Q * with probability 4 The derived performance bound shows that SQL has a rate of convergence of T = O(nβ 4 log(n/δ)/ǫ 2 ), which matches the proposed lower bound of reinforcement learning in terms of n, ǫ and δ. However, the dependency of our bound on the horizon β is worse than than the lower-bound by a factor of O(β 2 ).

ηN := 1 N k (Xt, At) + 1 ; if N k (x, a) > 0 then T k Q k-1 (Xt, At) := (1 -ηN )T k Q k-1 (Xt, At) + ηN (r(Xt, At) + γMQ k-1 (y k )); T k Q k (Xt, At) := (1 -ηN )T k Q k (Xt, At) + ηN (r(Xt, At) + γMQ k (y k )); else T k Q k-1 (Xt, At) := r(Xt, At) + γMQ k-1 (y k ); T k Q k (Xt, At) := r(Xt, At) + γMQ k (y k ); end Q k+1 (Xt, At) := (1 -α k )Q k (Xt, At) + α k kT k Q k (Xt, At) -(k -1)T k Q k-1 (Xt, At) ; // SQL update rule N k (Xt, At) := N k (Xt, At) + 1; Xt+1 = y k ; if min (x,a)∈Z N k (x, a) > 0 then // Check if all (x, a) ∈ Z have been visited at round k k := k + 1; α k := α0 k + 1 ; foreach (x, a) ∈ Z do N k (x, a) := 0; end t := t + 1; until t ≥ T ; return Q k (w.p.) 1 -δ.
Theorem 2 Let A1 hold, T be a positive integer and Q T be the estimate of Q * by Algorithm 1 at time step T . Then,the uniform approximation error

Q * -Q T ≤ β 2   γn T + 2n log 2n δ T   , w.p. at least 1 -δ.
Theorem 3 Let A1 and A2 hold and T > 0 be the number of time steps. Then, at step T of Algorithm 2 w.p. at least 1 -δ

Q * -Q T ≤ β 2   γeL log 2 δ T + 2eL log 2 δ log 4n δ T   .
These results, combined with Borel-Cantelli lemma [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF], guarantee that Q T converges almost surely to Q * with the rate 1/T for both Algorithm 1 and 2. Further, the following PAC bounds which quantify the number of steps T required to reach the error ǫ > 0 in estimating the optimal action-value function, w.p. 1 -δ, are immediate consequences of Theorem 2 and 3, respectively.

Corollary 4 (Finite-time PAC bound of synchronous SQL) Under A1, after

T = 3.74nβ 4 log 2n δ ǫ 2 steps (transitions) the uniform approximation error of Algorithm 1 Q * -Q T ≤ ǫ, w.p. at least 1 -δ.
Corollary 5 (Finite-time PAC bound of asynchronous SQL) Under A1 and A2, after

T = 3.74eLβ 4 log 2 δ log 4n δ ǫ 2 steps (transitions) the uniform approximation error of Algorithm 2 Q * -Q T ≤ ǫ, w.p. at least 1 -δ.
The following general result provides a new lower bound on the number of transitions T for every reinforcement learning algorithm to achieve ǫ-optimal performance w.p. 1 -δ under the assumption that the reinforcement learning algorithm is (ǫ, δ)-correct on the class of MDPs M. 5 4. We report the proofs in Section 5. 5. This result improves on the state-of-the-art [START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF] in terms of dependency on β (by a factor of O(β 2 )). Also, our result is more general than those of [START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF] in the sense that It does not require a sequential update of the value functions or following a deterministic policy.

Definition 6 ((ǫ, δ)-correct reinforcement learning) Define A as the class of all reinforcement learning algorithms which rely on estimating the action-value function Q * . Let Q A T be the estimate of Q * using T ≥ 0 number of transition samples under the algorithm A ∈ A. We then call A ∈ A (ǫ, δ)-correct on the class of MDPs M, if there exists some

0 < T * < +∞ such that, for all M ∈ M, Q * -Q A T ≤ ǫ w.p. at least 1 -δ iff T > T * .
Theorem 7 (Lower bound on the sample complexity of reinforcement learning) There exists some ǫ 0 , γ 0 , n 0 , c 1 , c 2 and a class of MDPs M such that for all ǫ ∈ (0, ǫ 0 ), γ ∈ (γ 0 , 1), n > n 0 , δ ∈ (0, 1) and every (ǫ, δ)-correct reinforcement learning algorithm A ∈ A on the class of MDPs M the number of transitions

T > T * = nβ 2 c 1 ǫ 2 log n c 2 δ .

Relation to Existing Results

In this section, we first compare our results for SQL with the existing results on the convergence of standard QL. The comparison indicates that SQL accelerates the convergence of Q-learning, especially for large values of β and small α's. We then compare SQL with batch Q-value iteration (QVI) in terms of sample and computational complexities, i.e., the number of samples and the computational cost required to achieve an ǫ-optimal solution with high probability, as well as space complexity, i.e., the memory required at each step of the algorithm.

A Comparison with the Convergence Rate of Standard Q-Learning

There are not many studies in the literature concerning the convergence rate of incremental model-free RL algorithms such as QL. Szepesvári (1997) has provided the asymptotic convergence rate for QL under the assumption that all the states have the same next state distribution. This result shows that the asymptotic convergence rate of QL with a linearly decaying learning step has exponential dependency on β, i.e. T = Õ(1/ǫ β ).

Even-Dar and Mansour (2003) throughly investigated the finite-time behavior of synchronous QL for different time scales. Their main result indicates that by using the polynomial learning step α k = 1 (k + 1) ω , 0.5 < ω < 1, the synchronous variant of Q-learning achieves ǫ-optimal performance w.p. at least 1 -δ after

T = O   n   β 4 log nβ δǫ ǫ 2 1 w + β log β ǫ 1 1-ω     , (3)
steps, where the time-scale parameter ω may be tuned to achieve the best performance. When γ ≈ 1, β = 1/(1 -γ) becomes the dominant term in the bound of Eq. 3, and thus, the optimized bound is achieved by finding an ω which minimizes the dependency on β. This leads to the optimized bound of Õ β 5 /ǫ 2.5 with the choice of ω = 4/5. On the other hand, SQL is guaranteed to achieve the same precision in only O β 4 /ǫ 2 steps. The difference between these two bounds is substantial for large β 2 /ǫ. Even-Dar and Mansour (2003) also proved bounds for the asynchronous variant of Qlearning in the case that the cover time of MDP can be uniformly bounded from above by some finite constant. The extension of their results to the more realistic case that the expected value of the cover-time is bounded by some L > 0 (assumption A2) leads to the following PAC bound: Proposition 8 (Even-Dar and Mansour, 2003) Under A1 and A2, for all ω ∈ (0.5, 1), after

T = O   L log 1 δ 1+3ω β 4 log nβ δǫ ǫ 2 1 w + Lβ log 1 δ log β ǫ 1 1-ω   steps (transitions) the uniform approximation error of asynchronous QL Q * -Q T ≤ ǫ, w.p. at least 1 -δ. The dependence on L in this algorithm is of O(L 3+ 1 ω +L 1 1-ω
), which leads, with the choice of ω ≈ 0.77, to the optimized dependency of O(L 4.34 ), whereas asynchronous SQL achieves the same accuracy after just O(L) steps. This result shows that for MDPS with large expected cover-time, i.e., slow-mixing MDPs, asynchronous SQL may converge substantially faster to a near-optimal solution than its Q-learning counterpart.

SQL vs. Q-Value Iteration

Finite sample bounds for both model-based and model-free (Phased Q-learning) QVI have been derived in [START_REF] Kearns | Finite-sample convergence rates for Q-learning and indirect algorithms[END_REF] and [START_REF] Even-Dar | PAC bounds for multi-armed bandit and Markov decision processes[END_REF]. These algorithms can be considered as the batch version of Q-learning. They show that to quantify ǫ-optimal action-value functions with high probability, we need Õ nβ 5 /ǫ 2 and Õ nβ 4 /ǫ 2 ) samples in model-free and model-based QVI, respectively. 6 A comparison between their results and the main result of this paper suggests that the sample complexity of SQL, which is of order Õ nβ 4 /ǫ 2 , 7 is better than model-free QVI in terms of β. Although the sample complexities of SQL and model-based QVI are of the same order, SQL has a significantly better computational and space complexity than model-based QVI: SQL needs only 2n memory space, while the space complexity of model-based QVI is given by min( Õ(nβ 4 /ǫ 2 ), n(|X| + 1)) (see [START_REF] Kearns | Finite-sample convergence rates for Q-learning and indirect algorithms[END_REF]. Also, SQL improves the computational complexity by a factor of Õ(β) compared to both model-free and model-based QVI. 8 Table 1 summarizes the comparisons between SQL and the other RL methods discussed in this section.

Experiments

In this section, we analyze empirically the effectiveness of the proposed algorithms on different problem domains. We examine the convergence of synchronous SQL (Algorithm 1) as well as the asynchronous SQL (Algorithm 2) on several discrete state-action problems 6. For the sake of simplicity, in this subsection, we ignore the logarithmic dependencies of the bounds. 7. Note that at each round of SQL n new samples are generated. This combined with the result of Corollary 5 deduces the sample complexity of order Õ(nβ 4 /ǫ 2 ). 8. SQL has the sample and computational complexity of a same order since it performs only one Q-value update per sample, the same argument also applies to the standard Q-learning, whereas, in the case of model-based QVI, the algorithm needs to iterate the action-value function of all state-action pairs at least Õ(β) times using Bellman operator, which leads to a computational complexity of order Õ(nβ 5 /ǫ 2 ) given that only Õ(nβ 4 /ǫ 2 ) entries of the estimated transition matrix are non-zero. [START_REF] Kearns | Finite-sample convergence rates for Q-learning and indirect algorithms[END_REF]. The source code of all tested algorithms are freely available at http://www.mbfys.ru.nl/~mazar/Research Top.html.

ǫ 2 Õ nβ 5 ǫ 2.5 Õ nβ 4 ǫ 2 Õ nβ 5 ǫ 2 CC Õ nβ 4 ǫ 2 Õ nβ 5 ǫ 2.5 Õ nβ 5 ǫ 2 Õ nβ 5 ǫ 2 SPC Θ(n) Θ(n) min( Õ(nβ 4 /ǫ 2 ), n(|X| + 1)) Θ(n) x 1 x 2500 x k x k+1 x k-1 a = -1 a = +1
Linear MDP: this problem consists of states x k ∈ X, k = {1, 2, . . . , 2500} arranged in a one-dimensional chain (see Figure 1). There are two possible actions A = {-1, +1} (left/right) and every state is accessible from any other state except for the two ends of the chain, which are absorbing states. A state x k ∈ X is called absorbing if P (x k |x k , a) = 1 for all a ∈ A and P (x l |x k , a) = 0, ∀l = k. Any transition to one of these two states has associated reward 1.

The transition probability for an interior state x k to any other state x l is inversely proportional to their distance in the direction of the selected action, and zero for all states corresponding to the opposite direction. Formally, consider the following quantity n(x l , a, x k ) assigned to all non-absorbing states x k and to every (x l , a) ∈ Z: We can write the transition probabilities as:

n(x l , a, x k ) =    1 |l -k| for (l -k)a > 0 0 otherwise . x 1 x 2500 x k x k+1 x k-1 a = -1 a = +1
P (x l |x k , a) = n(x l , a, x k ) xm∈X n(x m , a, x k ) .
Any transition that ends up in one of the interior states has associated reward -1.

The optimal policy corresponding to this problem is to reach the closest absorbing state as soon as possible.

Combination lock: the combination lock problem considered here is a stochastic variant of the reset state space models introduced in [START_REF] Koenig | Complexity analysis of real-time reinforcement learning[END_REF], where more than one reset state is possible (see Figure 2).

In our case we consider, as before, a set of states x k ∈ X, k ∈ {1, 2, . . . , 2500} arranged in a one-dimensional chain and two possible actions A = {-1, +1}. In this problem, however, there is only one absorbing state (corresponding to the state lockopened) with associated reward of 1. This state is reached if the all-ones sequence {+1, +1, . . . , +1} is entered correctly. Otherwise, if at some state x k , k < 2500, action -1 is taken, the lock automatically resets to some previous state x l , l < k randomly (in the original combination lock problem, the reset state is always the initial state x 1 ).

For every intermediate state, the rewards of actions -1 and +1 are set to 0 and -0.01, respectively. The transition probability upon taking the wrong action -1 is, as before, inversely proportional to the distance of the states. That is

n(x k , x l ) =    1 k -l for l < k 0 otherwise , P (x l |x k , 0) = n(k, l) xm∈X n(k, m) .
Note that this problem is more difficult than the linear MDP since the goal state is only reachable from one state, x 2499 .

Grid world: this MDP consists of a grid of 50 × 50 states. A set of four actions {RIGHT, UP, DOWN, LEFT} is assigned to every state x ∈ X. The location of each state x of the grid is determined by the coordinates c x = (h x , v x ), where h x and v x are some integers between 1 and 50. There are 196 absorbing firewall states surrounding the grid and another one at the center of grid, for which a reward -1 is assigned. The reward for the firewalls is

r(x, a) = - 1 c x 2 , ∀a ∈ A.
Also, we assign reward 0 to all of the remaining (non-absorbing) states.

This means that both the top-left absorbing state and the central state have the least possible reward (-1), and that the remaining absorbing states have reward which increases proportionally to the distance to the state in the bottom-right corner (but are always negative).

The transition probabilities are defined in the following way: taking action a from any non-absorbing state x results in a one-step transition in the direction of action a with probability 0.6, and a random move to a state y = x with probability inversely proportional to their Euclidean distance 1/ c x -c y 2 .

The optimal policy then is to survive in the grid as long as possible by avoiding both the absorbing firewalls and the center of the grid. Note that because of the difference between the cost of firewalls, the optimal control prefers the states near the bottom-right corner of the grid, thus avoiding absorbing states with higher cost.

Experimental Setup and Results

We describe now our experimental setting. The convergence properties of SQL are compared with two other algorithms: a Q-learning (Even-Dar and Mansour, 2003) (QL) and the model-based Q-value iteration (QVI) of [START_REF] Kearns | Finite-sample convergence rates for Q-learning and indirect algorithms[END_REF]. QVI is a batch reinforcement learning algorithm that first estimates the model using the whole data set and then performs value iteration on the learned model. All algorithms are evaluated in terms of ℓ ∞ -norm performance loss of the action-value function Q * -Q T at time-step T . We choose this performance measure in order to be consistent with the performance measure used in Section 3.2. The optimal action-value function Q * is computed with high accuracy through value iteration.

We consider QL with polynomial learning step α k = 1/(k+1) ω where ω ∈ {0.51, 0.6, 0.8} and the linear learning step α k = 1/(k + 1). Note that ω needs to be larger than 0.5, otherwise QL can asymptotically diverge (see [START_REF] Even-Dar | Learning rates for Q-learning[END_REF], for the proof).

To have a fair comparison of the three algorithms, since each algorithm requires different number of computations per iteration, we fix the total computational budget of the algorithms to the same value for each benchmark. The computation time is constrained to 30 seconds in the case of linear MDP and the combination lock problems. For the grid world, which has twice as many actions as the other benchmarks, the maximum run time is fixed to 60 seconds. We also fix the total number of samples, per state-action, to 1 × 10 5 samples for all problems and algorithms. Significantly less number of samples leads to a dramatic decrease of the quality of the obtained solutions using all the approaches. Algorithms were implemented as MEX files (in C++) and ran on a Intel core i5 processor with 8 GB of memory. CPU time was acquired using the system function times() which provides process-specific CPU time. Randomization was implemented using gsl rng uniform() function of the GSL library, which is superior to the standard rand(). 9 Sampling time, which is the same for all algorithms, were not included in CPU time. At the beginning of every run (i) the action-value functions are randomly initialized in the interval [-V max , V max ], and (ii) a new set of samples is generated from P (•|x, a) for all (x, a) ∈ Z. The corresponding results are computed after a small fixed amount of iterations.

Figure 3 shows the performance-loss in terms of elapsed CPU time for the three problems and algorithms with the choice of β = 1000. First, we see that, in all cases, SQL improves the performance very fast, almost by an order of magnitude, in just a few seconds. SQL outperforms both QL and QVI in all the three benchmarks. The minimum and maximum errors are attained for the combination lock problem and the grid world, respectively. We also observe that the difference between the final outcome of SQL and Q-learning (second best method) is very significant, about 30 times in all domains.

Figure 4 and 5 show means and standard deviations of the final performance-loss, respectively, over 50 runs, as a function of the horizon β. We observe that for large values of β, i.e. β ≥ 100, SQL outperforms other methods by more than an order of magnitude in terms of both mean and standard deviation of performance loss. SQL performs slightly worse than QVI for β ≤ 10. However, the loss of QVI scales worse than SQL with β, e.g., for β = 1000 SQL has an advantage of two order of magnitude over QVI. QL performs for larger values of ω when the horizon β is small, whereas for large values of β smaller ω's are preferable.

9. http://www.gnu.org/s/gsl. These results show that, as proved before in Theorem 2, SQL manages to average out the simulation noise caused by sampling and converges, rapidly, to a near optimal solution, which is robust in comparison to other methods. In addition, we can conclude that SQL significantly improves the computational complexity of learning w.r.t. the standard QL and QVI in the three presented benchmarks for our choice of experimental setup.

Analysis

In this section, we give some intuition about the convergence of asynchronous variant of SQL and provide the full proof of the finite-time analysis reported in Theorem 2 and Theorem 3. Then, we prove Theorem 7 (RL lower bound) in Subsection 5.1. We start by introducing some notations.

Let Y k be the set of all samples drawn at round k of the SQL algorithms and F k be the filtration generated by the sequence {Y 0 , Y 1 , Y 2 , . . . , Y k }. We first notice that for all (x, a) ∈ Z, the update rule of Eq. 2 may be rewritten in the following more compact form:

Q k+1 (x, a) = (1 -α k )Q k (x, a) + α k D k [Q k , Q k-1 ](x, a), where D k [Q k , Q k-1 ](x, a) 1 α k [(1 -α k )T k Q k (x, a) -(1 -2α k )T k Q k-1 (x, a)] and α k = 1/(k + 1).
We then define the operator

D[Q k , Q k-1 ]
as the expected value of the empirical operator

D k conditioned on F k-1 : D[Q k , Q k-1 ](x, a) E( D k [Q k , Q k-1 ](x, a)| F k-1 ) = E 1 -α k α k T k Q k (x, a) - 1 -2α k α k T k Q k-1 (x, a) F k-1 = 1 -α k α k TQ k (x, a) - 1 -2α k α k TQ k-1 (x, a),
where the last line follows by the fact that, in both Algorithm 1 and 2, T k Q k (x, a) and T k Q k-1 (x, a) are unbiased empirical estimates of the Bellman optimality operators TQ k (x, a) and TQ k-1 (x, a), respectively. Thus, the update rule of SQL can be re-expressed as

Q k+1 (x, a) = (1 -α k )Q k (x, a) + α k (D[Q k , Q k-1 ](x, a) -ǫ k (x, a)) , (4)
where the estimation error ǫ k is defined as the difference between the operator

D[Q k , Q k-1 ] and its sample estimate D k [Q k , Q k-1 ] for all (x, a) ∈ Z: ǫ k (x, a) D[Q k , Q k-1 ](x, a) -D k [Q k , Q k-1 ](x, a).
We have the property that E[ǫ k (x, a)|F k-1 ] = 0 which means that for all (x, a) ∈ Z the sequence of estimation error {ǫ 1 (x, a), ǫ 2 (x, a), . . . , ǫ k (x, a)} is a martingale difference sequence w.r.t. the filtration F k . Finally, we define the martingale E k (x, a) to be the sum of the estimation errors:

E k (x, a) k j=0 ǫ j (x, a), ∀(x, a) ∈ Z. (5)
The following steps lead to the proof of Theorem 2 and Theorem 3 (i) Lemma 9 shows the stability of SQL (i.e., the sequence of Q k stays bounded). (ii) Lemma 10 states the key property that the iterate Q k+1 is close to the Bellman operator T applied to the previous iterate Q k plus an estimation error term of order E k /k. (iii) By induction, Lemma 11 provides a performance bound Q * -Q k in terms of a discounted sum of the cumulative estimation errors {E j } j=0:k-1 . The above-mentioned results hold for both Algorithm 1 and Algorithm 2. Subsequently, (iv) we concentrate on proving Theorem 2 by making use of a maximal Azuma's inequality, stated in Lemma 13. (v) We then extend this result for the case of asynchronous SQL by making use of the result of Lemma 14.

For simplicity of the notations, we remove the dependence on (x, a), e.g., writing Q for Q(x, a) and E k for E k (x, a), when there is no possible confusion. Also, we notice that, for all k ≥ 0, the following relations hold between α k and α k+1 in Algorithm 1 and Algorithm 2:

α k+1 = α k α k + 1 and α k = α k+1 1 -α k+1 .
Lemma 9 (Stability of SQL) Let A1 hold and assume that the initial action-value function Q 0 = Q -1 is uniformly bounded by V max = β, then we have

Q k ≤ V max , ǫ k ≤ 2V max , and D k [Q k , Q k-1 ] ≤ V max ∀k ≥ 0. Proof We first prove that D k [Q k , Q k-1 ] ≤ V max by induction.
For k = 0 we have:

D 0 [Q 0 , Q -1 ] = T 0 Q -1 ≤ r + γ MQ -1 ≤ R max + γV max = V max .
Now for any k ≥ 0, let us assume that the bound

D k [Q k , Q k-1 ] ≤ V max holds. Thus D k+1 [Q k+1 , Q k ] = 1-α k+1 α k+1 T k+1 Q k+1 - 1-2α k+1 α k+1 T k+1 Q k 1-α k+1 α k+1 - 1-2α k+1 α k+1 r + γ 1-α k+1 α k+1 MQ k+1 - 1-2α k+1 α k+1 MQ k ≤ r + γ 1-α k+1 α k+1 M((1 -α k )Q k + α k D k [Q k , Q k-1 ]) - 1-2α k+1 α k+1 MQ k = r + γ 1- α k α k +1 α k α k +1 M((1 -α k )Q k + α k D k [Q k , Q k-1 ]) - 1-2 α k α k +1 α k α k +1 MQ k = r + γ M 1-α k α k Q k + D k [Q k , Q k-1 ] -1-α k α k MQ k ≤ r + γ M 1-α k α k Q k + D k [Q k , Q k-1 ] -1-α k α k Q k ≤ r + γ D k [Q k , Q k-1 ] ≤ R max + γV max = V max ,
and by induction, we deduce that for all k ≥ 0,

D k [Q k , Q k-1 ] ≤ V max . Now the bound on ǫ k follows from ǫ k = E(D k [Q k , Q k-1 ]|F k-1 ) -D k [Q k , Q k-1 ] ≤ 2V max , and the bound Q k ≤ V max is deduced by noticing that Q k = 1 k k-1 j=0 D j [Q j , Q j-1 ]. The next lemma shows that Q k is close to TQ k-1 , up to a O( 1 k ) term minus the cumu- lative estimation error 1 k E k-1 .
Lemma 10 Under A1, for any k ≥ 1:

Q k = TQ k-1 + 1 k (TQ 0 -TQ k-1 -E k-1 ) . ( 6
)
Proof We prove this result by induction. The result holds for k = 1, where (6) reduces to (4). We now show that if the property (6) holds for k ≥ 1 then it also holds for k + 1. Assume that (6) holds for k. Then, from (4) we have:

Q k+1 = (1 -α k )Q k + α k 1 -α k α k TQ k - 1 -2α k α k TQ k-1 -ǫ k = (1 -α k ) (TQ k-1 + α k-1 (TQ 0 -TQ k-1 -E k-1 )) + α k 1 -α k α k TQ k - 1 -2α k α k TQ k-1 -ǫ k = (1 -α k ) TQ k-1 + α k 1 -α k (TQ 0 -TQ k-1 -E k-1 ) + (1 -α k )TQ k -(1 -2α k )TQ k-1 -α k ǫ k = (1 -2α k )TQ k-1 + α k (TQ 0 -E k-1 ) + (1 -α k )TQ k -(1 -2α k )TQ k-1 -α k ǫ k = (1 -α k )TQ k + α k (TQ 0 -E k-1 -ǫ k ) = TQ k + α k (TQ 0 -TQ k -E k ) = TQ k + 1 k + 1 (TQ 0 -TQ k -E k ) .
Thus (6) holds for k + 1, and is thus true for all k ≥ 1.

Now we bound the difference between Q * and Q k in terms of the discounted sum of the cumulative estimation errors {E 0 , E 1 , . . . , E k-1 }.

Lemma 11 (Error Propagation of SQL) Let A1 hold and assume that the initial actionvalue function Q 0 = Q -1 is uniformly bounded by V max = β, then for all k ≥ 1, we have:

Q * -Q k ≤ 1 k   2γβ 2 + k j=1 γ k-j E j-1   (7) ≤ β k 2γβ + max j=1:k E j-1 . (8)

Proof

We first notice that k j=1 γ k-j E j-1 ≤ β max j=1:k E j-1 for any sequence of cumulative errors errors {E 0 , E 1 , E 2 , . . . , E k-1 }. Therefore, we only need to prove ( 7) and ( 8) follows.

Again we prove this lemma by induction. The result holds for k = 1 as:

Q * -Q 1 = TQ * -TQ 0 -ǫ 0 ≤ Q * -Q 0 + ǫ 0 ≤ 2γV max + ǫ 0 ≤ 2γβ 2 + E 0 .
We now show that if the bound holds for k, then it also holds for k + 1. Thus, assume that (7) holds for k. By using Lemma 10:

Q * -Q k+1 = Q * -TQ k - 1 k + 1 (TQ 0 -TQ k -E k ) ≤ α k (TQ * -TQ 0 ) + (1 -α k )(TQ * -TQ k ) + α k E k ≤ α k TQ * -TQ 0 + (1 -α k ) TQ * -TQ k + α k E k ≤ 2γα k V max + γ(1 -α k ) Q * -Q k + α k E k .
By plugging ( 7) into (9) we then deduce:

Q * -Q k+1 ≤ 2α k γV max + γ(1 -α k )α k-1   2γβ 2 + k j=1 γ k-j E j-1   + α k E k ≤ 2α k γβ + γ(1 -α k ) α k 1 -α k   2γβ 2 + k j=1 γ k-j E j-1   + α k E k ≤ 2α k γβ + γα k 2γβ 2 + γα k k j=1 γ k-j E j-1 + α k E k = 1 k + 1   2γβ 2 + k+1 j=1 γ k+1-j E j-1   .
Thus, the result holds for k + 1 thus for all k ≥ 1 by induction.

In the next lemma, we prove a high probability bound over the estimation error term of Lemma 11 Lemma 12 Let A1 hold and assume that the initial action-value function Q 0 = Q -1 is uniformly bounded by V max = β then for all 0 < δ < 1 and k ≥ 1, we have:

max j=1:k E j-1 ≤ β 8k log 2n δ , w.p. 1 -δ. ( 9
)
Proof We start by providing a high probability bound for max 1≤j≤k |E j-1 (x, a)| for a given (x, a). First notice that:

P max 1≤j≤k |E j-1 (x, a)| > ǫ = P max max 1≤j≤k (E j-1 (x, a)), max 1≤j≤k (-E j-1 (x, a)) > ǫ = P max 1≤j≤k (E j-1 (x, a)) > ǫ max 1≤j≤k (-E j-1 (x, a)) > ǫ ≤ P max 1≤j≤k (E j-1 (x, a)) > ǫ + P max 1≤j≤k (-E j-1 (x, a)) > ǫ , (10
) and each term is now bounded by using a maximal Azuma-Hoeffding's inequality, reminded now (see e.g., Cesa-Bianchi and Lugosi ( 2006)).

Lemma 13 (Maximal Hoeffding-Azuma's Inequality) Let V = {V 1 , V 2 , . . . , V k } be a martingale difference sequence w.r.t. a sequence of random variables {X 1 , X 2 , . . . , X k } (i.e., E(V j+1 |X 1 , . . . X j ) = 0 for all 0 < j ≤ k) such that V is uniformly bounded by L > 0. If we define S k = k i=1 V i , then for any ǫ > 0, we have:

P max 1≤j≤k S j > ǫ ≤ exp -ǫ 2 2kL 2 .
As mentioned earlier, the sequence of random variables {ǫ 0 (x, a), ǫ 1 (x, a), • • • , ǫ k (x, a)} is a martingale difference sequence w.r.t. the filtration F k (generated by the random samples {y 0 , y 1 , . . . , y k }(x, a) for all (x, a)), i.e., E[ǫ k (x, a)|F k-1 ] = 0. It follows from Lemma 13 that for any ǫ > 0 we have:

P max 1≤j≤k (E j-1 (x, a)) > ǫ ≤ exp -ǫ 2 8kV 2 max P max 1≤j≤k (-E j-1 (x, a)) > ǫ ≤ exp -ǫ 2 8kV 2 max . (11)
By combining ( 11) with (10) we deduce that P (max

1≤j≤k |E j-1 (x, a)| > ǫ) ≤ 2 exp -ǫ 2 8kV 2 max ,
and by a union bound over the state-action space, we deduce that

P max 1≤j≤k E j-1 > ǫ ≤ 2n exp -ǫ 2 8kV 2 max . (12)
This bound can be rewritten, for any δ > 0, as:

P max 1≤k≤T E j-1 ≤ V max 8k log 2n δ ≥ 1 -δ,
which completes the proof.

The result of Theorem 2 then follows by plugging (9) into (8) and taking into account that after T steps of Algorithm 1 we have T ≤ n(k + 1). For the proof of Theorem 3 we rely on the following result which bounds the number of steps required to visit all states-actions k times with high probability :

Lemma 14 Under A2, from any initial state x 0 and for any integer k > 0, after a run of T = keL log 1 δ steps the state-action space Z is covered at least k times under the policy π w.p. at least 1 -δ.

Proof

We begin by defining the random event Q k as the number of steps required to cover the whole MDP for k times starting from any state x 0 ∈ X at any time t > 0. We then bound Q k for any x 0 in high probability using Markov inequality [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF]:

P(Q k > ekL) ≤ E (Q k ) ekL ≤ k sup t>0 max x∈X E(τ π (x, t)) ekL ≤ kL ekL = 1 e .
In words after a run of length ekL the probability that the entire state-action space is not covered for at least k times is less than 1 e . The fact that the bound holds for any initial state and time implies that after m > 0 intervals of length ekL the chance of not covering Z for k times is less than 1 e m : In the sequel, we restrict ourselves to the case that the only unknown factor is the statetransition probability p M . Therefore, the RL task is to learn only those value functions which depend on p M , i.e., Q * (z) for all z ∈ Z 0 . We also denote by E m ad P m the expectation and the probability under the model M m in the rest of this section.

P(Q k > mekL) ≤ 1 e m .
We follow the following steps in the proof: (i) we prove a lower-bound for the samplecomplexity of each state-action pair z ∈ Z 0 on the class of MDP M * . (ii) we then make use of the fact that the estimates of Q * (z) for different z ∈ Z 0 are independent of each others to combine these bounds and prove the result of Theorem 7.

We begin by proving that for all RL algorithms A ∈ A, there exists an MDP M m ∈ M * and constants c 0 > 0 and c > 0 such that for all ǫ < ǫ 0 = 1/8:

P m (|Q * (z) -Q A tz (z)|) > ǫ) > θ z /c 0 , (13)
where θ z exp -ct z ǫ 2 /β 2 and Q A tz (z) is an empirical estimate of the action-value function Q * (z) by the RL algorithm A using t z > 0 transition samples from the state-action pair z ∈ Z 0 . To prove this lower bound we make use of a contradiction argument, i.e., we assume that there exists an algorithm A ∈ A for which:

P m (|Q * (z)-Q A tz (z)|) > ǫ) ≤ θ z /c 0 , or equivalently P m (|Q * (z)-Q A tz (z)|) ≤ ǫ) ≥ 1-θ z /c 0 ,
for all M m ∈ M * and show that this assumption leads to a contradiction. To prove this we need to, first, introduce some notations: We define the event Ω 1 (z) {|Q * 0 -Q A t (z)| ≤ ǫ} for all z ∈ Z, where Q * 0 γ/(2(1 -γ)) is the optimal action-value function for all z ∈ Z 0 under the MDP M 0 . We then define K t r 1 + r 2 + • • • + r t as the sum of rewards of trying z ∈ Z 0 for t times. We also introduce the event Ω 2 (z), for all z ∈ Z 0 :

for some positive quantities T , ξ and δ ≤ 0.397ψ, then:

T ≤ ψ c 1 η log ψ c 2 δ ,
with c 1 = c and c 2 = c 0 /1.78.

Proof

In order to prove this lemma we make use of the inequality exp(x) ≥ 1 + x:

1 - 1 c 0 exp - cξT ψ ≤ (1 -δ) 1 ψ ≤ exp -δ ψ ≤ 1 - δ 1.78ψ
,

where in the last step we rely on Lemma 16, which imposes the condition δ ≤ 0.397ψ. We then deuce by collecting terms:

δ 1.78ψ ≤ 1 c 0 exp - cξT ψ . ( 17
)
The result then follows by solving (17) w.r.t. T .

We now make use of Eq. ( 13), Lemma 18 and Lemma 19, also the fact that for every z 1 , z 2 ∈ Z 0 the random events {|Q * (z 1 ) -Q A tz 1 (z 1 )| ≤ ǫ} and {|Q * (z 2 ) -Q A tz 2 (z 2 )| ≤ ǫ} are independent of each others,11 to prove that there exists an MDP in the class M * such that:

P( Q * -Q A T ≤ ǫ) ≤ P   z∈Z 0 |Q * (z) -Q A tz (z)| ≤ ǫ   = z∈Z 0 P(|Q * (z) -Q A tz (z)| ≤ ǫ) = z∈Z 0 1 - 1 c 0 exp -ct z ǫ 2 β 2 ≤ 1 - 1 c 0 exp -cT ǫ 2 β 2 n 1 n 1 ≤ 1 - 1 c 0 exp -cT ǫ 2 β 2 n/3 n/3 ≤ 1 -δ,
with T ≤ β 2 n/(c 1 ǫ 2 ) log(n/(c 2 δ)), δ ≤ 0.0662n, c 1 = 3c = 2100 and c 2 = c 0 /1.78 = 2.24.12 In words, if the total number of samples T is less than β 2 n/(c 1 ǫ 2 ) log(n/(c 2 δ)) then the probability of Q * -Q A T ≤ ǫ is at maximum 1 -δ on either M 0 or M 1 . This is equivalent to the statement that for every RL algorithm A ∈ A to be (ǫ, δ)-correct on the set M * , and subsequently on the class of MDPs M, the total number of transitions T needs to satisfy T > β 2 n/(c 1 ǫ 2 ) log(n/(c 2 δ)), which concludes the proof of Theorem 7.

Conclusions and Future Work

In this paper, we introduced a new reinforcement learning algorithm, called speedy Qlearning (SQL). We analyzed the finite-time behavior of this algorithm as well as its asymptotic convergence to the optimal action-value function. Our result is in the form of high probability bound on the performance loss of SQL, which suggests that the algorithm converges to the optimal action-value function in a faster rate than the standard Q-learning. The numerical experiments in Section 4 confirm our theoretical results showing that for large value of β SQL outperforms the other RL methods by a wide margin. Overall, SQL is a simple, efficient and theoretically well-founded reinforcement learning algorithm, which improves on existing RL algorithms such as Q-learning and the sample-based value iteration.

In this work, we are only interested in the estimation of the optimal action-value function and not the problem of exploration. Therefore, we did not compare our result to the PAC-MDP methods (Szita and Szepesvári, 2010;[START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF] and the upper-confidence bound based algorithms [START_REF] Jaksch | Near-optimal regret bounds for reinforcement learning[END_REF][START_REF] Bartlett | REGAL: A regularization based algorithm for reinforcement learning in weakly communicating MDPs[END_REF], in which the choice of the exploration policy impacts the behavior of the learning algorithms. However, we believe that it would be possible to gain w.r.t. the state of the art in PAC-MDPs, by combining the asynchronous version of SQL with a smart exploration strategy. This is mainly due to the fact that the bound for SQL has been proved to be tighter than the RL algorithms that have been used for estimating the value function in PAC-MDP methods, especially in the model-free case. Also, SQL has a better computational requirement in comparison to the standard RL methods. We consider this as a subject for future research.

Another possible direction for future work is to scale up SQL to large (possibly continuous) state and action spaces where function approximation is needed. We believe that it would be possible to extend our current SQL analysis to the continuous case along the same path as in the fitted value iteration analysis by [START_REF] Munos | Finite-time bounds for fitted value iteration[END_REF] and [START_REF] Antos | Fitted Q-iteration in continuous action-space MDPs[END_REF]. This would require extending the error propagation result of Lemma 11 to a ℓ 2 -norm analysis and combining it with the standard regression bounds.

  a), ∀(x, a) ∈ Z. The optimal action-value function Q * is the unique fixed-point of the Bellman optimality operator T defined as: (TQ)(x, a) r(x, a) + γ y∈X P (y|x, a)(MQ)(y), ∀(x, a) ∈ Z.

Algorithm 1 :

 1 Synchronous speedy Q-learning Input: Initial action-values Q0, discount factor γ and number of steps T Q-1 := Q0; // Initalization t := k := 0; repeat // Main loop

Algorithm 2 :

 2 Asynchronous speedy Q-learning Input: Initial action-values Q0, the policy π(•|•), discount factor γ, number of step T and the initial state X0. t := k := 0; // Initialization α0 = 1; foreach (x, a) ∈ Z do Q-1(x, a) := Q0(x, a); N0(x, a) := 0; end repeat // Main loop Draw the action At ∼ π(•|Xt); Generate the next state sample y k by simulating P (•|Xt, At);

Figure 1 :

 1 Figure 1: Linear MDP: Illustration of the linear MDP problem. Nodes indicate states. States x 1 and x 2500 are the two absorbing states and state x k is an example of interior state. Arrows indicate possible transitions of these three nodes only. From x k any other node is reachable with transition probability (arrow thickness) proportional to the inverse of the distance to x k (see the text for details).

Figure 2 :

 2 Figure 2: Combination lock: illustration of the combination lock MDP problem. Nodes indicate states. State x 2500 is the goal (absorbing) state and state x k is an example of interior state. Arrows indicate possible transitions of these two nodes only. From x k any previous state is reachable with transition probability (arrow thickness) proportional to the inverse of the distance to x k . Among the future states only x k+1 is reachable (arrow dashed).

  Figure 3: A comparison between SQL, QL and model-based QVI. Each plot compares the performance loss of the policy induced by the algorithms for a different MDP averaged over 50 different runs (see the text for details).

  Figure 4: A comparison between SQL, QL and the model-based Q-value iteration (QVI) given a fixed computational and sampling budget. The plot shows means at the end of the simulations for three different algorithm in terms of the horizon of MDP β.

Figure 5 :

 5 Figure 5: A comparison between SQL, QL and the model-based Q-value iteration (QVI) given a fixed computational and sampling budget. The plot shows variances at the end of the simulations for three different algorithms in terms of the horizon of MDP β.

Figure 6 :

 6 Figure 6: Illustration of the class of MDP considered for the lower-bounds. Nodes indicate state-(action) pairs. Arrows indicate possible transitions of these nodes. (see the text for details).

Table 1 :

 1 Comparison between SQL, Q-learning, model-based and model-free Q-value iteration in terms of sample complexity (SC), computational complexity (CC), and space complexity (SPC).

	Method	SQL	Q-learning	Model-based QVI	Model-free QVI
	SC	Õ nβ 4

In this work, for the sake of simplicity in notation, we assume that the reward r(x, a) is a deterministic function for all (x, a) ∈ Z. However, it is not difficult to extend our results to the case of stochastic reward under some mild assumptions, e.g., boundedness of the absolute value of reward.

Our results also hold if the rewards are taken from some interval [rmin, rmax] instead of [0, 1], in which case the bounds scale with the factor rmax -rmin.

It is important to note that T is a γ-contraction mapping w.r.t. to the ℓ∞-norm, i.e., for any pair of action-value functions Q and Q ′ , we have TQ -TQ ′ ≤ γ Q -Q ′(Bertsekas, 2007b, Chap. 1).

This is due to the disjoint structure of the class M.

The result holds for all δ ∈ (0, 1) with the choice of n ≥ n0 = 16.

By the choice of m = log 1 δ we deduce:

The bound can be then rewritten as follows:

which concludes the proof.

Plugging the results of Lemma 14 and Lemma 12 into (8) concludes the proof of theorem 3.

Proof of Theorem 7

In this subsection, we provide the full proof of Theorem 7. In our analysis, we rely on the likelihood-ratio method, which has been perviously used to prove a lower-bound for multi-armed bandits, (see [START_REF] Mannor | The sample complexity of exploration in the multi-armed bandit problem[END_REF], and extend this approach for reinforcement learning in Markovian decision problems. We also make use of some technical lemmas in [START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF]. We begin by defining a class of MDPs for which the proposed lower-bound holds. The class of MDPs M (see Figure 6) is defined as a set of all MDPs with n = 3n 1 ≥ 3 state-action pairs, in which the set of state-action pairs Z is made up of three smaller sets Z 0 , Z 1 and Z 2 each of size n 1 . We make the assumption that, for all M ∈ M, both Z 1 and Z 2 consist of only absorbing states, see Section 4 for the definition of absorbing states. We also assume that for all M ∈ M, every entry z l 0 , 1 ≤ l ≤ n 1 , is connected to only two other state-action pairs z l 1 ∈ Z 1 and z l 2 ∈ Z 2 , such that the probability of moving from z l 0 to z l 1 is given by p M (z l 0 ), for all 1 ≤ l ≤ n 1 . Also, we assume that the instant reward r(z) is set to 1 for all z ∈ Z 1 and 0 for the rest of the state-actions. One can then solve the Bellman equation, For all z ∈ Z, in close-form and compute the optimal action-value function Q * as follows:

Here, we only consider a subset of 2 MDPs in the class M denoted by M * = {M 0 , M 1 }. For all M m ∈ M * and z ∈ Z 1 , the state transition probability p M (z) is given by:

We also define S(z)

The following lemmas are also required for our analysis. 10 Lemma 15 [START_REF] Mannor | The sample complexity of exploration in the multi-armed bandit problem[END_REF]) For all ǫ < ǫ 0 = 1/8 and γ > γ 0 = 1/2 we have:

where d = 1.78. Now, we proceed with the proof of the lower-bound. By the assumption that P m (|Q * (z)-

with the choice of c 0 = 4. This combined with Lemma 15 implies that P 0 (S(z)) > 1/2, for all z ∈ Z. Based on this result, we prove the following lemma:

Lemma 17 For all z ∈ Z 0 : P 1 (Ω(z)) > θ z /c 0 .

Proof We define W as the history of all the outcomes of trying z for t z times and the likelihood function L m (w) for all M m ∈ M * as:

for every possible history w and M m ∈ M * . This function can be used to define a random variable L(W m ). The likelihood ratio of the event w between two MDPs M 1 and M 0 can be written as:

where k z is a short-hand notation for K tz .

We then determine the lower bounds of the terms in the RHS of ( 14), when the event S = Ω 1 ∩ Ω 2 occurs:

If S(z) occurs, then Ω 1 (z) occurs, and we have k ≤ t. Therefore, we deduce:

10. For the proofs see [START_REF] Mannor | The sample complexity of exploration in the multi-armed bandit problem[END_REF].

Likewise, if the event S(z) occurs, then Ω 2 (z) occurs which leads to:

which implies:

(16) By plugging ( 15) and ( 16) into ( 14), we deduce:

Then, by the choice of c = 700, we obtain:

where 1 S is the indicator function of the event S(z). Then by a change of measure we deduce:

where we make use of the fact that P 0 (S(z)) > 1 2 .

Lemma 17 states that, for all z ∈ Z 0 , P 1 (Ω 1 (z)) > θ z /c, which violates the assumption that, for all z ∈ Z 0 and M m ∈ M * , P m (|Q * (z) -Q A tz (z)|) ≤ ǫ) ≥ 1 -θ z /c 0 under Algorithm A ∈ A, due to the fact Ω 1 (z) and {|Q * (z) -Q A tz (z)| ≤ ǫ} are two separate events. The contradiction between the result of Lemma 17 and the assumption which leads to this result proves the lower bound of Eq. ( 13).

Based on the lower-bound of Eq. ( 13), we can prove the lower-bound on the total number of samples T = z∈Z 0 t z . To prove the lower bound, we rely on the following lemmas from [START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF]:

Lemma 18 [START_REF] Strehl | Reinforcement learning in finite MDPs: PAC analysis[END_REF] Let c and ∆ be constants in (0, 1), Y be a finite set with the cardinality N and t y > 0 for all y ∈ Y. Define T = y∈Y t y . We then have :


Cuong Cao Dang

Vincent Lefort

Vinh Sy Le

Quang Si Le

Olivier Gascuel
email: olivier.gascuel@lirmm.fr

ReplacementMatrix: a Web Server for Maximum Likelihood Estimation of Amino Acid Replacement Rate Matrices

 and provides a number of tools and facilities. Users upload a set of multiple protein alignments from their domain of interest and receive the resulting matrix by email, along with statistics and comparisons with other matrices. A nonparametric bootstrap is performed optionally to assess the variability of replacement rate estimates. Maximum-likelihood trees, inferred using the estimated rate matrix, are also computed optionally for each input alignment. Finely-tuned procedures and up-to-date ML software (PhyML 3.0, XRATE) are combined to perform all these heavy calculations on our clusters.

INTRODUCTION

Amino-acid replacement matrices contain estimates of the instantaneous substitution rates from any amino acid to another. These rates reflect the biological, chemical, and physical properties of amino acids. For example, we usually observe a high substitution rate between lysine (positively charged) and arginine (also positively charged) and a low substitution rate between lysine and aspartate (negatively charged). Amino-acid replacement matrices are an essential basis of protein phylogenetics. They are used to compute substitution probabilities along phylogeny branches, and thus the likelihood of the data. They are also closely related to score matrices, which are essential for aligning proteins and computing alignment scores.

*

To whom correspondence should be addressed.

Several general replacement matrices have been proposed, such as PAM [START_REF] Dayhoff | A model of evolutionary change in proteins[END_REF], JTT [START_REF] Jones | The rapid generation of mutation data matrices from protein sequences[END_REF], WAG [START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF], and LG [START_REF] Le | An improved general amino acid replacement matrix[END_REF]. These matrices were estimated from large and diverse sets of protein alignments. They tend to be robust and perform well in many cases. However, the performance of replacement matrices depends on life domains and protein groups [START_REF] Keane | Assessment of methods for amino acid matrix selection and their use on empirical data shows that ad hoc assumptions for choice of matrix are not justified[END_REF]). Replacement matrices have thus been estimated for specific domains (e.g. for HIV, [START_REF] Nickle | HIV-specific probabilistic models of protein evolution[END_REF]influenza, Dang et al. 2010), and protein groups (e.g. mitochondrial proteins, [START_REF] Adachi | Model of amino acid substitution in proteins encoded by mitochondrial DNA[END_REF]. It has been shown that specific replacement matrices often differ significantly from general matrices, and thus perform better when applied to the data to which they are dedicated (e.g. [START_REF] Adachi | Model of amino acid substitution in proteins encoded by mitochondrial DNA[END_REF]Hasegawa 1996, Dang et al. 2010).

Since the seminal work of [START_REF] Dayhoff | A model of evolutionary change in proteins[END_REF], a number of methods have been designed to estimate amino-acid replacement matrices from protein alignments. These methods belong to either counting (e.g. [START_REF] Jones | The rapid generation of mutation data matrices from protein sequences[END_REF] or maximum likelihood (ML) approaches (e.g. [START_REF] Adachi | Model of amino acid substitution in proteins encoded by mitochondrial DNA[END_REF][START_REF] Yang | Models of amino acid substitution and applications to Mitochondrial protein evolution[END_REF][START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF]. The former are limited to pairwise protein alignments, while the latter fully benefit from the information contained in multiple alignments and the corresponding phylogenies. Recently, we improved the ML method proposed by [START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF] by incorporating the variability of evolutionary rates across sites into the matrix estimation process [START_REF] Le | An improved general amino acid replacement matrix[END_REF]. This procedure was successfully applied to estimate the LG matrix from 3,912 alignments of the Pfam database, the FLU matrix from 992 influenza protein alignments, and a number of matrices corresponding to various structural configurations of the residues [START_REF] Le | Accounting for Solvent Accessibility and Secondary Structure in Protein Phylogenetics is Clearly Beneficial[END_REF].

The demand to estimate amino-acid replacement matrices for particular data is rising quickly because of the rapidly growing volume of sequence data and a desire to better understand the evolution and relationships of specific protein groups and species. However, up-to-date replacement matrix estimation procedures are complex and highly demanding in computational terms. Our method [START_REF] Le | An improved general amino acid replacement matrix[END_REF] involves complex data processing and alternates tree building using PhyML [START_REF] Guindon | New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0[END_REF]) and matrix estimation using XRATE [START_REF] Klosterman | XRate: a fast prototyping, training and annotation tool for phylo-grammars[END_REF]. It thus requires a huge amount of work to estimate a matrix from raw data sets. Here we describe an implementation of this method in a Web server. Users upload their alignments and receive the output matrix by email along with a number of additional statistics and comparisons. Optionally, the server performs a non-parametric bootstrap to assess the variability of rate estimations, and infers the phylogeny of every input alignment using the estimated replacement matrix.

MODEL AND METHODS

The amino-acid substitution process is assumed to be independent among sites and lineages, and homogeneous during the course of evolution. The standard model is Markovian, time-continuous, time-reversible and represented by a 20 × 20 rate matrix

  ij Q q 
, where is the number of substitutions from amino acid i to amino acid j per time unit. The diagonal elements ii are such that the row sums are all zero. Any time-reversible matrix can be decomposed into a symmetric exchangeability matrix

( ij q i j  ) q Q   ij
and an amino-acid equilibrium frequency vector

R r   i    , using equality ij j . Moreover, Q is normalized, that is (i j    ) ij q r 1 i ii q     .
Here we consider (as usual) the most general time-reversible model (GTR), which involves 189 (R) and 19 () free parameters to be estimated from the data (see textbooks for additional explanation, e.g. [START_REF] Felsenstein | Inferring phylogenies[END_REF].

Given a set of protein alignments

  a D D  , Q is estimated by maxi- mizing the likelihood     a a
a , where the product runs over all alignments a L D L T ,    ,Q;D D and the inner term is the likelihood of a D given the phylogenetic tree a , the rate across site model a , and the replacement matrix Q. Here we use the standard discrete gamma distribution with 4 rate categories, and is the gamma parameter associated with

T  a a  D .
Simultaneously optimizing T, Q and  parameters is computationally difficult. However, several authors have showed that substitution model parameters (Q and ) can be accurately estimated using nearly optimal trees T. [START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF] estimated their WAG matrix by: (1) inferring tree topologies using NJ, (2) estimating tree branch lengths by ML assuming a JTT replacement process, and (3) estimating Q from the data and thereby inferred trees using a standard optimization procedure.

We refined this approach by incorporating an across-site rate model in the matrix estimation, namely 4 gamma categories plus invariant sites (4+I). Our method [START_REF] Le | An improved general amino acid replacement matrix[END_REF] involves: (1) estimating tree topologies and branch-lengths using PhyML [START_REF] Guindon | New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0[END_REF], (2) processing alignment and trees to account for the rate model, (3) estimating Q from these processed data and trees using the expectation-maximization software XRATE [START_REF] Klosterman | XRate: a fast prototyping, training and annotation tool for phylo-grammars[END_REF], and (4) iterating this procedure until   L D reaches a plateau. This estimation procedure is started using an approximate matrix. WAG was used to learn LG, and a nearly identical matrix was obtained when starting from JTT. We observed that 3 iterations are enough in practice and that the invariant site category has little impact on Q estimation.

The above procedure is very heavy in computational terms. It is simplified here. The most time-consuming aspect is the ML estimation of tree topologies, which is performed only once here (instead of ~3 times in the original procedure). Moreover, the rate model is simplified by using 4 gamma rate categories but no invariant sites (4). The resulting matrix is nearly the same as that obtained using the full procedure (see results below) but the run time is 2-3 times faster. The simplified procedure has three main estimation steps (1, 2 and 3) and is as follows:

Step 0 Input a set of multiple alignments and a starting replacement matrix S; only exchangeabilities in S are used, frequencies are estimated from the data.

Step 1 (a) For each alignment, build a BioNJ tree and optimize the branch lengths and gamma rate parameter using PhyML with S and 4.

(b) Process the alignments and trees to account for the 4 model: every alignment is divided into 4 sub-alignments using the posterior probability of site rate categories, and the 4 corresponding trees are rescaled using the rates estimated for each category under the gamma model.

(c) Run XRATE with default options and S starting matrix to estimate a first matrix from the processed alignments and trees. (c) Same as (1c), but replace S by and output .

1

Q 2 Q Step 3 (a)
For each alignment, re-optimize the branch-lengths of the previously inferred ML tree and gamma rate parameter using PhyML with and 4. (c) Same as (1c), but replace S by and output final Q matrix.

2

Q

Step 4 For each alignment, re-optimize the branch-lengths of the previously inferred ML tree and the gamma rate parameter using PhyML with Q, with S, and with LG when ; output the corresponding log-likelihood and AIC values of every alignment and site for comparison purposes.

LG S  Only Step (2) in this procedure fully constructs an ML tree; Step (1) uses a distance-based tree topology (as with WAG estimation), while Step (3) reuses the ML topology inferred during Step (2) with a fairly accurate 1 matrix. Other parts are the same as in the original LG estimation procedure (except for the invariant site category, removed here).

Q

When the final matrix has been estimated, it is returned along with a number of results, statistics, and comparisons. Two additional options are available: (1) performing a bootstrap study to assess the variability of rate estimates; and (2) running PhyML 3.0 with Q and standard options to infer the phylogenies estimated with the new matrix for all input alignments. When the latter option is used, the pipeline simultaneously estimates the replacement matrix and the trees from the input alignments. These are expected to be significantly different from the phylogenies inferred with starting matrix S or LG. To save computing time, the starting trees and initial parameter values are taken from Step (4) in the above procedure.

The aim of the bootstrap procedure is to measure the variability of rate estimations. This should be useful, for example, when comparing the properties of amino acids in specific contexts [START_REF] Kosiol | A new criterion and method for amino-acid classification[END_REF], or when using replacement rate matrices in the search for non-standard genetic codes [START_REF] Abascal | MtArt: a new model of amino acid replacement for Arthropoda[END_REF]). The bootstrap is performed in a standard manner: for every alignment a D in D we draw with replacement a D sites and then run the estimation procedure to obtain a pseudo rate matrix; this is repeated several times and the pseudo matrices are used to compute several statistics (e.g. the standard deviation) for each of the frequency   i  and exchangeability   ij r parameters. This procedure is highly time-consuming, and we thus only perform 10 replicates. Moreover, the estimation scheme described above is still too heavy to be repeated 10 times. We therefore use the trees and site rate categories computed by PhyML with Q in Step (4), and run XRATE only once for each replicate, starting from the S matrix. Experimental studies show that these simplifications do not significantly affect the variability measures.

RESULTS WITH TWO SAMPLE DATA SETS

To illustrate the properties of the Web server, we re-estimated the LG matrix from the data used in original publication (3,912 alignments, ~6.5 millions residues) and the FLU matrix using 100 randomly selected alignments from the original data set (~1.8 million residues). We performed a bootstrap with 500 (LG) and 1,000 (FLU100) replicates to obtain accurate measures of the variability of parameter estimates, and 20 standard pipeline bootstrap runs with 10 replicates each. Detailed results are available as Supplementary Material from the Web site and summarized in Table 1.

We see that the new LG matrix estimated by the Web server is nearly identical to the published matrix, despite the simplifications in the estimation procedure. The new FLU100 matrix (estimated from 100 alignments) is also very close to the original one (estimated from 992 alignments). The relative deviations of equilibrium frequencies ( i i   ) are quite low, while those of exchangeabilities ( ij ij r 

) are higher, especially with FLU100 where their average is nearly equal to 20%. This finding shows that exchangeabilities are much more difficult to estimate than frequencies. Exchangeabilities measure instantaneous rates of change, which are not directly observable from the data (as frequencies are) and may correspond to hidden changes between ancestral states. With highly conserved alignments as FLU100's, some amino-acid pairs are rarely seen together in the same alignment site (e.g. Cys-Lys is present 4 times only among ~37,000 sites, see Sup. Mat.), and thus estimating their exchangeabilities is inevitably a difficult task. Lastly, we see that the bootstrap variability measures with 10 replicates are clearly correlated with those obtained using a large number of replicates (500 and 1000), and should thus be useful for analyzing the differences between rates or between matrices. ): average relative deviation of the frequencies (exchangeabilities) obtained with 500 (LG) and 1,000 (Flu100) bootstrap replicates.

R i i   ij ij r    2 i i R     2 ij ij R r  LG 0.

 

2 i i R   (   ij ij R r  2
): average and standard deviation (among 20 trials) of the Pearson correlations of the relative deviations of frequencies (exchangeabilities) computed with 10 replicates, and those computed with 500 (LG) and 1,000 (FLU100) replicates.

WEB SERVER, INPUT AND OUTPUT FILES

The main input is a set of multiple alignments in PHYLIP or Fasta format. This typically contains hundreds or even thousands of alignments. However, each alignment must contain less than 100 sequences to reduce the computational burden. Larger alignments must be divided into several sub-alignments that are given separately. A starting replacement matrix may also be provided, otherwise LG is the default. Two options allow for bootstrapping and running PhyML with the estimated matrix. The user receives a job status URL and the estimated matrix by email along with a number of files and statistics. These include (see user guide for details):

 The new rate matrix in PAML triangular format, where exchangeability   ij r and frequency   i  parameters are given separately. These parameter values are compared to those of the starting matrix S and of LG (when S ≠ LG), using Pearson correlation, histograms, and bubble graphs.

 A series of score matrices for various evolutionary distances (), derived from the rate matrix using standard log-odds:

    i i j j  
, where the probability of change from i to j given  is calculated by exponentiation of the rate matrix [START_REF] Felsenstein | Inferring phylogenies[END_REF]. As with PAM matrices, the distance ranges from 0.10 (corresponding to PAM10) to 2.5 (PAM250). These matrices can be used, for example, with MAFFT, CLUSTALW, or BLAST to search for homologs or compute multiple alignments of specific protein groups.

 The fit of the new rate matrix to the input data is compared to that of S and of LG (when S ≠ LG), using the log-likelihood difference on the whole data set, divided by the total number of sites. To account for the fact that the new matrix is estimated from these data and thus has to be penalized for its (189+19) additional parameters, we use the AIC difference divided by the number of sites. The AIC and log-likelihood differences are also provided for every alignment and every site, for example to detect atypical alignments or site classes. Z-tests are used to assess the significance of AIC differences.

When the bootstrap and/or PhyML options are selected, the user receives separate emails providing:

 The standard deviation, the relative deviation, the minimum and the maximum values (among 10 bootstrap estimates) for each of the frequency and exchangeability parameters.

 All trees inferred by PhyML 3.0 using the new matrix with SPR and standard options for each of the input alignments.

The current waiting time when all options are selected is ~10 days for the very large LG data set, and ~2 days with FLU100.

1 Q

 1 Step 2 (a) For each alignment, infer an ML tree using PhyML 3.0 with , 4 and the SPR tree search option. 1 Q (b) Same as (1b).

Table 1 :

 1 Results with LG and FLU100 data sets 2

ACKNOWLEDGMENTS

Thanks to Ian Holmes and Christophe Dessimoz for their help, and to the Vietnam National Foundation for Science and Technology Development and ANR SysBio MytoSys for financial support.