A Distributional Analysis of Sampling-Based Reinforcement Learning Algorithms
Philip Amortila Doina Precup Prakash Panangaden Marc G. Bellemare
cs.LG cs.AI stat.ML














 Philip Amortila, Doina Precup, Prakash Panangaden, Marc G. Bellemare

[

A Distributional Analysis of Sampling-Based Reinforcement Learning Algorithms

  Philip Amortila^α Doina Precup^α,β Prakash Panangaden^α Marc G. Bellemare^α,β,γ 
 ^αMcGill University;  ^βCIFAR Fellow;  ^γGoogle Research] 



We present a distributional approach to theoretical analyses of reinforcement learning algorithms for constant step-sizes. We demonstrate its effectiveness by presenting simple and unified proofs of convergence for a variety of commonly-used methods. 
We show that value-based methods such as TD(λ) and Q-Learning have update rules which are contractive in the space of distributions of functions, thus establishing 
their exponentially fast convergence to a stationary distribution. 
We demonstrate that the stationary distribution obtained by any algorithm whose target is an expected Bellman update has a mean which is equal to the true value function.
Furthermore, we establish that the distributions concentrate around their mean as the step-size shrinks.
We further analyse the optimistic policy iteration algorithm, for which the contraction property does not hold, and formulate a probabilistic policy improvement property which entails the convergence of the algorithm. 





§ INTRODUCTION


Basic results in the theory of Markov decision processes (MDPs) and dynamic programming (DP) rely on the two fundamental properties of the Bellman operator: contraction and monotonicity. For instance, proofs of convergence for value iteration and policy iteration follow immediately from the contractive properties of the Bellman operators and the Banach fixed point theorem <cit.>.

However, proving the convergence of sample-based algorithms such as TD-learning <cit.> or optimistic policy iteration <cit.> requires substantially more effort. The typical stochastic approximation approach relies on hitting-time or martingale arguments
to bound the sequence of value function iterates within progressively smaller regions (see, e.g., <cit.>).

In this work we present a distributional framework for analyzing sample-based reinforcement learning algorithms. Rather than consider the evolution of the random point estimate produced by the learning process, we study the dynamics of the distribution of these point estimates. As a concrete example, we view the TD(0) algorithm as defining a sequence of random iterates (V_n)_n ∈ whose distributions are recursively defined by the distributional equation

    V_n+1(s) D= (1-α) V_n(s) + α(R(s,A) + γ V_n(S')) ,

where s is the initial state and (A, R, S') is the random action-reward-next-state transition sampled from the underlying Markov Decision Process. The equation recursively describes the distribution of a random variable undergoing the TD learning process.

We study the constant step-size case. Our main contribution is to show that, for a variety of algorithms, the random iterates converge in distribution to a fixed point of the corresponding distributional equation, even though the random point estimate may not converge. We further characterize this fixed point, showing that it depends on both the step-size and the specific Markov Decision Process under consideration. 
Our framework views the learning process as defining a time-homogeneous Markov chain over the space of value functions. We prove convergence by establishing the existence of a stationary distribution for this Markov chain and demonstrating that the sequence of random iterates generated by a sample-based learning rule must converge to this stationary distribution, using tools from optimal transport <cit.>.

We first analyze sample-based algorithms whose corresponding distributional operator is a contraction mapping in the infinity norm, including TD(λ), Q-learning, and double Q-learning. Following a proof technique of <cit.>, we lift these stochastic algorithms to the distributional setting. We show that this lifting recovers contractive guarantees, now in the Wasserstein metric using the infinity norm as a cost function. The contraction coefficient depends on the discount factor, as usual, but also on the step-size: updates with smaller step-sizes converge more slowly to their distributional fixed point. TD(0), for example, is a contraction mapping with coefficient 1 - α + αγ.

We also analyze the sample-based equivalent of policy iteration, called optimistic policy iteration  tsitsiklis2002convergence or Monte Carlo control <cit.>. The convergence of policy iteration is not driven by a contraction mapping, but rather by the monotonicity of the policy iteration operator <cit.>. We derive a similar, weaker property for the sample-based setting which we call probabilistic policy improvement. We use this property to show that optimistic policy improvement also converges to a distributional fixed point.

By recovering the contraction mapping that underlies many dynamic programming algorithms, our distributional analysis significantly simplifies existing proofs of convergence for stochastic algorithms, at least for constant step-sizes. Our approach easily allows us to quantify the limiting behaviour of these algorithms; the same tool even provides us with confidence bounds over the true value function. We believe this type of analysis should prove useful going forward, including for the study of reinforcement learning with function approximation.




§ BACKGROUND



We write () for the set of probability distributions on a space . We consider an agent interacting with an environment modelled as a finite Markov decision process (,,,,γ). As usual,  is a finite state space,  is a finite set of actions, : ×→([0,Rmax]) is a bounded reward distribution function, : ×→() is a transition distribution function, and γ∈ [0,1) is a discount factor. The strategy of the agent is captured by a policy π : →(). The value function v^π : → of a policy π is the expected discounted sum of rewards observed when starting at state s and following policy π. The value function is the fixed point of the Bellman operator ^π defined by 

    ^π v(s) _   a∼π(·|s)
       r∼(·|s,a)[r + γ_s' ∼(·|s,a)[v(s')]].

The value function of the optimal policy π^⋆ is also the fixed point of the Bellman optimality operator ^⋆, defined by 

    ^⋆ v(s) max_a{_   r∼(·|s,a)
       s'∼(·|s,a)[r + γ v(s')]}.

A closely-related object is the action-value function q^π, the expected discounted return of first taking action a and thereafter following policy π. The action-value function satisfies the Bellman equations q^π(s,a) = ^π q^π(s,a) and q^⋆(s,a) = ^⋆ q^⋆(s,a), where ^π and ^⋆ are defined analogously to Equations (<ref>) and (<ref>) sutton1998introduction.
The Bellman operators for value functions (resp. action-value functions) are contractions on ^|| (resp. ^||×||) with respect to the infinity norm v := v_∞ = max_i |v_i|  puterman1994markov. A policy π is called greedy with respect to an action-value function Q ∈^||×|| if π(s) ∈_aQ(s,a) for each s ∈.



 §.§ Couplings and the Wasserstein Metric


To establish convergence in distribution, we will use the Wasserstein metric  between distributions Villani08. As a cost function, we use the infinity norm. For two distributions μ,ν∈(^)̣, a pair of random vectors (X,Y) is a coupling of (μ,ν) if X ∼μ and Y ∼ν. We write Ξ(μ,ν) for the set of such couplings. The Wasserstein metric on (^)̣ with the infinity norm as a cost function is defined as: 

    (μ,ν) = inf_(X,Y) ∈Ξ(μ,ν)[X-Y_∞].

The metric is defined over the set (^)̣ = {μ∈(^)̣: ∫x_∞μ(x) < + ∞} of measures with finite first moment. Assuming a bounded reward function, we will always be dealing with finite-moment measures. 
The Wasserstein metric characterizes the weak convergence of measures [Theorem 6.9]Villani08, or equivalently the convergence in distribution of the associated random variables.




§ MARKOV CHAINS ON THE SPACE OF FUNCTIONS

With many value-based RL algorithms, the stochasticity of the algorithm depends only on the sampled transition and the random current estimate. For example, recalling the update rule for TD(0) (Equation (<ref>)), the value of V_n+1(s) for a particular state is fully determined by knowledge of V_n and the action, reward, and successor state which was sampled from s:

    {V_n+1| V_n, V_n-1,...,V_1,V_0} = {V_n+1| V_n}.

We therefore view these methods as inducing Markov chains on the space of value functions. We note that their state space is continuous rather than discrete – we take it to be ^|| when modelling value functions or ^||×|| when modelling action-value functions. When results hold for both cases, we will write the discussion in terms of ^, ∈̣. Whenever needed, we may also restrict ourselves to the subset of realizable functions [0, Rmax/1 - γ]^⊂^.

For a given update rule  and step-size α, the transition function for an induced Markov chain is as follows. Given f_k ∈^$̣, letf_k+1be the random function obtained by applyingwith step-sizeα. For a Borel set∈(^), we define the Markov kernelK_,αas:K_,α(f_k,ℬ) = {f_k+1∈ℬ | f_k }.This Markov kernel describes the probability of transitioning fromf_kto some function in the setunder the update rule. In the sequel, we omit the subscripts on the kernel when the update rule is clear from context. 
 For a given probability measureμ∈(^), the distribution of functions after one transition of the Markov chain is given byμ K () = ∫_^ K(θ,)μ(θ).The distribution of functions afterntransitions is given byK^n, which is defined inductively as:K^n(θ, ) = ∫_^ K(θ',)K^n-1(θ,θ').A probability measureψis a stationary distribution for a Markov chain with kernelKifψ= ψ K.An algorithm updates synchronously when all states or state-actions pairs are updated at every iteration. In the regime of constant step-sizes and synchronous updates, the Markov kernels are time-homogeneous (or time-independent). Thus, the lawμ_n() = f_n ∈of the random variablef_nis given by:μ_n←μ_0 (K)^n.§.§ Stochastic operators

In this section, we introduce the notion of a stochastic operator and provide a general formalism for the analysis of stochastic update rules. We will distinguish two classes of stochastic operators which will require different analyses.

We model the sampling space as a probability space(Ω, ,η). A stochastic operator is a map between functions which depends on a randomly sampled eventω∈Ω. 

A stochastic operator is a function : ^×̣Ω→^$̣.

When operating on functions, a stochastic operatoroutputs a random function. We will write a number of stochastic value-based algorithms as

    f_n+1 = (1-α) f_n + α(f_n,ω ),

wheref_n,f_n+1∈^$̣ are functions, α is a step-size, and  is some algorithm-dependent stochastic operator. In this notation, the operator  is the target of the algorithm. We say that  is an empirical Bellman operator if it behaves like a Bellman operator in expectation.

The stochastic operator  is an empirical Bellman operator for a policy π if 
_ω∼η[(f,ω)] = ^π f   ∀ f ∈^.̣
Similarly,   is an empirical Bellman optimality operator if _ω∼η[(f,ω)] = ^⋆ f.


In general, the sampling distribution of the stochastic operator may depend on the function which it is acting on. 
Two examples of methods for which the sampling distribution is independent of the current function estimates are TD(0), which applies an empirical Bellman operator, and Q-Learning, which applies an empirical Bellman optimality operator. TD(0) is defined by the stochastic operator 

    (V,(a_s,r_s,s'_s)_s ∈)(s) = r_s + γ V(s'_s),
 where (a_s,r_s,s'_s) is a transition sampled for every state, and Q-Learning is defined by the operator

    (Q,(r_s,a,s'_s,a)_(s,a))(s,a) = r_s,a + γmax_a'Q(s'_s,a,a'),

where (r_s,a,s'_s,a) is a transition sampled for every state-action pair. Convergence of these methods is covered in Section <ref>. On the other hand, methods for which the sampling of the update rule depends on the function being updated are more akin to policy iteration, which applies Bellman operators that depend on the current greedy policy. We will see an example of such a method in Section <ref>. 



§ CONVERGENCE VIA CONTRACTION TO A STATIONARY DISTRIBUTION

In this section we demonstrate that common value-based algorithms converge to a stationary distribution when updated synchronously and with constant step-sizes. The convergence follows by showing that their Markov kernels are contractive with respect to the Wasserstein metric. To illustrate our approach, we provide a proof of convergence for TD(0). With the same proof method, we also establish convergence and give convergence rates for Monte Carlo evaluation, Q-Learning, TD(λ), SARSA, Expected SARSA van2009theoretical, and Double Q-Learning hasselt2010double. The proofs for these other algorithms are given in Appendix <ref>.

Recall the update rule of the synchronous TD(0) algorithm given by Equation (<ref>).

For any step size 0<α≤ 1, the TD(0) algorithm has a contractive Markov kernel K_α:

    (μ K_α, ν K_α) ≤ (1-α+αγ)(μ,ν),

for all μ,ν∈(^||).

Let μ^(1), μ^(2)∈(^||) be two distributions of function estimates. Let V^(1)_0 ∼μ^(1), V^(2)_0 ∼μ^(2) be the coupling which minimizes the Wasserstein metric, i.e.:
 (μ^(1),μ^(2)) = inf_(X,Y)[ X-Y ] = [V^(1)_0-V^(2)_0]. 
 Such an optimal coupling always exists [Theorem 4.1]Villani08. We couple the updates (V^(1)_1,V^(2)_1) to sample identical transitions at each state, i.e.:

    V^(1)_1(s)    = (1-α)V^(1)_0(s) + α(r_s + γ V^(1)_0(s'_s)) 
    
    V^(2)_1(s)    = (1-α)V^(2)_0(s) + α(r_s + γ V^(2)_0(s'_s)),

for the samea ∼π(·|s),r_s ∼(·|s,a), and s'_s ∼(·|s,a). Note that this is a valid coupling of (μ^(1)K_α,μ^(2)K_α) since V_1^(1) and V_1^(2) sample transitions from the same distributions. 
We upper-bound ( μ^(1)K_α,  μ^(2)K_α) with the coupling above. 

    (μ^(1) K_α,μ^(2) K_α)    ≤𝔼[V^(1)_1-V^(2)_1] 
       ≤ (1-α)[ V^(1)_0-V^(2)_0] 
       +α[max_s |(r_s-r_s) + γ( V^(1)_0(s'_s) - V^(2)_0(s'_s)) |]

We note that the expectation is over the pair (V^(1)_0,V^(2)_0) as well as the random samples a_s,r_s,s'_s. By our coupling construction,

    [max_s |(r_s-r_s) + γ( V^(1)_0(s'_s) - V^(2)_0(s'_s)) |] 
       =γ[ max_s | V^(1)_0(s'_s) - V^(2)_0(s'_s) |] 
       ≤γ[ max_s  | V^(1)_0(s) - V^(2)_0(s) |] = γ(μ^(1),μ^(2))

The inequality follows since V_0^(1) and V_0^(2) sample the same set of successor states – the maximum is the same if each s samples a different s'_s and is lesser otherwise. Using Equation (<ref>) in Equation (<ref>) gives:

    (μ^(1) K_α,μ^(2) K_α) ≤ (1-α + αγ)(μ^(1),μ^(2)).
 
Since 1-α + αγ< 1, the kernel K_α  is a contraction mapping.


The contraction property readily entails the convergence to a stationary distribution. We initialize with any V_0 drawn from an arbitrary distribution of finite first moment. 

For any constant step size 0<α≤ 1 and initialization V_0 ∼μ_0 ∈(^||), the sequence of random variables (V_n)_n ≥ 0 defined by the recursion (<ref>) converges in the Wasserstein metric to a unique stationary distribution ψ^TD(0)_α∈(^||).

 The space of probability measures (^||) metrized with  is a complete metric space (<cit.>), and therefore it follows from Banach's fixed point theorem  that (μ_0 K_α^n)_n ≥ 0 converges exponentially quickly to a unique fixed point ψ^TD(0)_α. The distribution ψ^TD(0)_α is a stationary distribution by the fixed point property: ψ^TD(0)_α K_α = ψ^TD(0)_α. 

As evidenced by the above, lifting the analysis to distributions over value functions greatly simplifies the proof. The key is in the choice of a proper coupling. The same technique extends to a broad class of algorithms, with relatively few modifications. This avoids, for example, the additional hurdles caused by the greedy probability kernel in Q-learning <cit.>. We further note some surprising connections with distributional reinforcement learning <cit.>. For α = 1, the fixed point of TD(0) is in fact 's return distribution. The same coupling, which forces two processes to sample the same transitions, has also been implicitly used to study the behaviour of distributional algorithms <cit.>.



To demonstrate the power of the approach, we summarize in Table <ref> a series of results regarding common sampling-based RL algorithms. Under similar conditions to Theorem <ref>, each algorithm listed in Table <ref> converges to a stationary distribution (which is in general different for different algorithms, as we show in the next section). Each proof only requires small adjustments to the basic proof template, for example an extended state space (Double Q-Learning). Full details, along with the proof template, are given in the appendix.




§ THE STATIONARY DISTRIBUTIONS

In this section, we characterize the stationary distributions which are attained by any algorithm whose target is a Bellman operator or Bellman optimality operator in expectation. In our notation, these algorithms are defined in terms of empirical Bellman (optimality) operators. As before, we write the discussion in terms of ^$̣ since results will hold for both value functions and action-value functions.

What do these distributions look like? We first consider the case of policy evaluation algorithms, which have as expected operator^π. In that case, their mean corresponds to the fixed point of^π, i.e. the value functionsv^πorq^π. Second, they concentrate around this mean in inverse proportion to the step-sizeα. Hence, as expected, small step sizes lead to a more accurate distribution at the cost of a larger contraction factor. The full distributions are not symmetric or easily described, however; as a simple example, takeα = 1in TD(0), corresponding to the return distribution <cit.>. In the case of optimality operators, we show that the mean of the stationary distributions is in fact greater than the fixed pointsv^⋆orq^⋆. 



 §.§ Sample-based Evaluation Algorithms

Suppose ^π is an empirical Bellman operator for some policy π and that the updates (<ref>) with step-size α converge to a stationary distribution ψ_α. Let f_α∼ψ_α and f^π be the fixed point of ^π. 
Then [f_α] = f^π. 

Let f_0 be distributed according to ψ_α. By stationarity,

    f_1  = (1-α)f_0 + α^π (f_0,ω)

is also distributed according to ψ_α. We write f_α[f_0]. Taking expectations on both sides,
and using that _ω[^π(f,ω))] = ^π(f) for any f: 

    f_α   = (1-α) f_α  + α_ω,ψ_α[^π(f_0,ω)] 
    f_α   = _ψ_α[^π f_0] 
    f_α   = ^π_ψ_α[f_0] = ^πf_α

And therefore f_α = f^π since it is the unique fixed point of  ^π.








We remark again that this characterization will hold for any algorithm which converges and performs Bellman updates in expectation. Although they have the same mean, the stationary distributions will depend on the update rule. These differences will be reflected in their higher moments. To this effect, we next derive a closed-form expression for the covariance of the stationary distribution. We writeA^for the transpose of a matrixA. The outer product of two vectorsx,y ∈^is the matrixxy^𝖳∈^×̣defined by(xy^𝖳)_i,j = x_iy_j.
Thus,[(X⃗-μ⃗)(X⃗-μ⃗)^]is the covariance of a random vectorX⃗with meanμ⃗. The proof of the following result is provided in Appendix <ref>. 

Let ^π be an empirical Bellman operator for some policy π. Suppose ^π is such that the updates (<ref>) with step-size α converge to a stationary distribution ψ_α. Define ξ_ω(f) = ^π(f,ω) - ^π f, and  
(f) _ω[ξ_ω(f)ξ_ω(f)^]
to be the covariance of the zero-mean noise term ξ_ω(f) for a given function f. Define C = (1-(1-α))^2. The covariance of f_α∼ψ_α is given by

    C[(f_α-f^π)(f_α-f^π)^]     = 
       α^2 (γ^π)[(f_α-f^π)(f_α-f^π)^](γ^π)^
        + α(1-α)(γ^π) [(f_0-f^π)(f_0-f^π)^] 
        + α(1-α) [(f_α-f^π)(f_α-f^π)^](γ^π)^
        + α^2 ∫(f)ψ_α(f).

Theorem <ref> provides a recursive definition for the covariance of the stationary distributionψ_α. The integral in the final line corresponds to the expected covariance of the empirical Bellman operator when sampling from the distribution. Under the assumption that states are updated independently, this “one-step” covariance is diagonal. More generally, the covariance matrix is scaled byα, suggesting that the distribution concentrates around its mean whenαis close to 0. The following makes this precise. We writeA_op=sup{Av : v≤1, v ∈^}for the operator norm of a matrixA. 


Assume that the state space of the Markov chain is bounded. Let C  (2Rmax/1-γ)^2. Then, we have that [(f_α-f^π)(f_α-f^π)^]_op is monotonically decreasing with respect to α.  In particular, lim_α→ 0[(f_α-f^π)(f_α-f^π)^]_op = 0, and we have that: 

    {min_i | f_α(i) - f^π(i)|≥ε}   ≤C/ε̣^2α^2/1-(1-α+αγ)^2
       α→ 0⟶ 0.


We remark that the boundedness of the state space (e.g. by[0,Rmax/1-γ]^⊂̣^$̣) is easily satisfied in the presence of bounded rewards in the MDP. Furthermore, the above results about the mean and covariance can easily be extended beyond Bellman operators to any operator which has a unique fixed point and commutes with expectation. 



 §.§ Sample-Based Control Algorithms


Above we saw that the mean of the stationary distribution of a sample-based method using a fixed policy is the value function for that policy. This no longer holds in the presence of optimality operators, for example in what is called the control setting<cit.>. To conclude this section, we use our distributional approach to highlight behavioural characteristics of control algorithms. 


Suppose ^⋆ is an empirical Bellman optimality operator such that the updates (<ref>) with step-size α converge to a stationary distribution ψ^⋆_α. Let f_α∼ψ^⋆_α  and f^⋆ is the fixed point of ^⋆. Then [f_α] ≥ f^⋆.  Equality holds if and only if the expectation and the maximum commute, i.e. f =  f
As before, let f_0 be distributed according to ψ^⋆_α. Taking expectations on both sides of f_1  = (1-α)f_0 + α^⋆(f_0,ω) and writing f_α[f_α] gives:

    f_α   = (1-α) f_α  + α_ω,f_0[^⋆(f_0,ω)] 
    f_α   = _f_0[max_π^π f_0] 
    f_α   ≥max_π_f_0[^π f_0 ] 
    f_α   ≥max_π^πf_α = ^⋆f_α

By the linear programming formulation of MDPs <cit.>, we conclude that f̅_̅α̅≥ f^⋆ = min_f {f ≥^⋆ f }.


The theorem shows that in general, sample-based control methods such as Q-learning produces a biased (in an expected sense) estimate of the optimal Q-value. This brings fresh evidence about the algorithm's well-known overestimation problem, which double Q-learning seeks to correct. 





















§ CONVERGENCE VIA MONOTONICITY: OPTIMISTIC POLICY ITERATION

In a previous section, we showed that a number of sampling-based algorithms induce a contraction mapping in the Wasserstein metric between distributions over value functions. In this section we analyze a non-contractive example, namely the optimistic policy iteration (OPI) algorithm. The OPI algorithm is a sampling-based analogue of the classic policy iteration (PI) algorithm. The latter is driven to convergence by the monotonicity of the greedy policy updates. We show in this section that our Markov chain approach can regain a distributional analogue of the monotonicity property, which we call probabilistic policy improvement, and that this property can be used to analyze the algorithm in a restricted setting. 

The convergence of optimistic policy iteration is more difficult to prove than that of most sample-based algorithms, and has been previously been established for Robbins-Monro decreasing stepsizes by using monotonicity arguments and assumptions on the sampling distribution <cit.>. To the best of our knowledge, the convergence of OPI for more general conditions (including constant step-sizes) remained an open problem.

Optimistic policy iteration proceeds by constructing a greedy policy from its current value function, sampling one trajectory per state-action pair from this policy, then updating its value function towards the return of these trajectories. 
We will write 

    ^π(s_0,a_0)=∑_t=0^∞γ^tr_t(s_t,a_t)

for a sampled discounted return starting at state s_0, taking first action a_0, and thereafter following policy π. For any Q ∈^||×||, we write π_Q for the greedy policy corresponding to Q (assuming a consistent tie-breaking so that this is well-defined). Let Q_0 be some initial estimate and π_0=π_Q_0. The update rule of OPI is as follows (α∈ (0, 1]):

    Q_n+1(s,a)    = (1-α)Q_n(s,a) + α^π_n(s,a) 
    π_n+1   = π_Q_n+1.

Analyzing optimistic policy iteration in the distributional setting poses a few challenges. First, the distribution of sampled trajectories depends on the exact value function. Informally, the greedy mapping from value functions to policies induces a greedy partition<cit.>, with a different empirical Bellman operator corresponding to each region of this partition. This rules out a simple coupling argument, as functions with different greedy policies may have arbitrarily different return distributions.  point out that optimistic policy iteration can lead to chattering, where the greedy policy fails to converge even the value function converges. For our analysis, we consider the simpler case α = 1; we discuss the extension to α < 1 at the end of the section.

For α = 1 and initialization Q_0 ∼μ_0 ∈(^||×||),
the sequence of random variables (Q_n)_n ≥ 0 defined by the recursion (<ref>)
converges to a unique stationary distribution φ_1 ∈(^||×||).

The key lemma is to extend the monotonicity property of the policy iteration operator to the distributional case. In policy iteration, the greedy policy π'=π_Q^π with respect to Q^π leads to an improved value function:

    Q^π'≥ Q^π .

The role of Q^π is therefore to provide us with the improved policy π'. We will show that the same holds true for the sampled returns: there is some probability that the greedy policy with respect to ^π is π'. This allows us to argue that there is a chance that optimistic policy iteration follows the correct “greedy path” to π^*. We will call this property probabilistic policy improvement. 

We analyze the case α=1 by considering a Markov chain over policies. Formally, Π={π: →} will be our state space, with the Markov kernel:

    K(π, π')    π'= π_^π
       = π' ^π .

This Markov chain reflects the OPI process since, at every step, the greedy policy π_n corresponding to Q_n is sufficient to determine the distribution of Q_n+1. Since the set Π of deterministic policies is finite, K is a discrete Markov chain.

Suppose 
π' = π_Q^π.
Then K(π,π') > 0.

The proof of Lemma <ref> is given in Appendix <ref>.
This shows that there is a nonzero probability that the chain improves on the current policy. 
This implies that there is some probability that OPI applied from π^* produces π^* as a greedy policy.

The optimal policy π^⋆ is aperiodic. In particular: K(π^⋆,π^⋆)>0. 

Since the optimal policy π^⋆ is greedy with respect to Q^⋆, from Lemma <ref> we conclude that K(π^⋆,π^⋆) > 0. 

All that remains to show is that the optimal policy π^⋆ is reachable from any other policy with positive probability.

For every π_0 ∈Π, there exists an n(π_0) ∈ such that K^n(π_0)(π_0, π^⋆) > 0.

Let π_0 be an initial policy. Let Q^π_0, Q^π_1,...,Q^π^⋆ be the sequence of action-value functions obtained from classical PI. Since PI converges in a finite number of steps (say n_π_0) and is a deterministic process, this sequence is well-defined. For every i ∈{1,..,n(π_0)}, we have that K(π_i,π_i+1) > 0 by Lemma <ref> (since π_i+1 is greedy with respect to Q^π_i by construction). Thus we have that K(π_0,π_1)K(π_1,π_2)⋯ K(π_n(π_0)-1,π^⋆) > 0 and in particular K^n(π_0)(π_0 , π^⋆ ) > 0.

Finally, the reachability and aperiodicity of π^⋆ allow us to apply the ergodic theorem for finite Markov chains. 

The policy π^⋆ must be contained in a communicating class 𝒞^⋆ of policies (perhaps consisting of only π^⋆) which is aperiodic since π^⋆ is. There may be other communicating classes in the Markov chain, but by Lemma <ref> they must all be transient since they can reach π^⋆. By the Markov chain convergence theorem [Theorem 4.3]levin2017markov, any initial distribution converges to a stationary distribution φ_1 ∈(Π) with support over 𝒞^⋆. 

Our result shows that optimistic policy iteration, applied with a step-size of α = 1, converges to a stationary distribution φ_1 over aperiodic policies (and thus to a stationary distribution over value functions through the possible returns of these policies).
Since K(π^*, π^*) < 1 in general, we know that this distribution has support on suboptimal policies; in fact, we know that

    φ_1(π^*) = 1/1 - K(π^*, π^*)∑_ππ^*φ_1(π) K(π, π^*) .

By “continuity”, this suggests that the algorithm should also converge for the general case α∈ [0, 1). Unfortunately, our proof technique does not immediately carry over. The issue is that, for α < 1, we no longer have a Markov chain over policies: the greedy policy depends on the history of past policies, through the value function. One path forward may be to study the Markov chain over value functions, but the known brittleness of optimistic policy iteration suggests that its distributional behaviour may be quite complex. In particular, the transition kernel fails to satisfy many basic properties (such as the weak Feller property) which are typically used to establish convergence in Markov chains over continuous spaces <cit.>. We leave as open questions whether the algorithm does converge, and to which distribution.





§ RELATED WORK

Some of our methods are inspired from the work of <cit.>, which develops the theory of constant step-size stochastic gradient descent (in the context of supervised learning). In particular, the proof method we present in Section <ref> is inspired from the proof of their Proposition 2, although simplified and adapted to the RL setting, and the results in Section <ref> follow the methods of their Proposition 3.

In RL, convergence in distribution results for constant step-sizes are typically derived using tools common to stochastic approximation theory such as the mean ODE method and Lyapunov functions (see, e.g., <cit.> and <cit.>). Examples of works which feature these methods include <cit.>. The results and methods of these works are different, as they neither exploit the Markov chain perspective nor establish the convergence of the iterates to a stationary distribution. 

Some works do make explicit use of the Markov chain perspective, most related are <cit.>. The first of these establishes the convergence of the Markov chains with respect to the Total Variation metric using tools from [Chapters 13-16]meyn2012markov. In applications to the analysis of RL algorithms, this type of convergence does not hold without restrictive assumptions such their Assumption (2.6) – see Appendix <ref> for a simple counterexample featuring a bandit with a single deterministic arm.
On the other hand, results about weak convergence of RL algorithms yu2016weak have established the convergence of the averaged iterates rather than the full sequence of distributions. The methods are also different, and rely on the weak Feller property meyn2012markov amongst other stochastic approximation techniques kushner2003stochastic.
As far as we are aware, the use of the Wasserstein metric and the result that RL algorithms are contractive with respect to this metric are novel.  


 



§ CONCLUSION AND FUTURE WORK

We studied the convergence properties of sample-based reinforcement learning algorithms by considering how they induce distributions over value functions.
Many of these algorithms are in fact contractive not in the space of functions but in the lifted space of distributions of functions. The proof methods relies on coupling the events sampled by two executions of the algorithm, and can be re-used for many algorithms. Using the same Markov chain approach, we also analyzed a restricted version of optimistic policy iteration, which is not amenable to a contraction mapping-type analysis. One of the key results is to make explicit that constant step-size reinforcement learning algorithms do converge, albeit in the weaker distributional sense. As an upside of using a constant step size, we obtain exponentially fast convergence (as indicated by the presence of a contraction factor). By controlling the step-sizes, the stationary distributions thus obtained can be tailored to yield values close to the true value function with high confidence. In the control setting, this should enable us to better explain the performance of practical reinforcement learning schemes.

Our work opens a number of interesting avenues for future research. First, it would be valuable to fully characterize the stationary distribution of sample-based methods, for example by deriving a closed-form expression for their characteristic functions. A deeper understanding of the distributions obtained by control algorithms is also of interest. Second, we did not analyze the case of decaying step-sizes or online updates, which would correspond to time-inhomogenenous Markov processes. More broadly, the coupling method has historically been invaluable for many applications in probability theory. It would be interesting to see if our approach can be applied to policy-based methods, for example policy gradient or actor critic, which are closer in spirit to optimistic policy iteration. Finally, the simplicity of our analysis suggests that it may be carried to the function approximation setting, perhaps eventually shedding light on the behaviour of reinforcement learning with nonlinear approximation methods such as deep networks.





§ ACKNOWLEDGEMENTS
 
We gratefully acknowledge funding from the CIFAR Learning in Machines and Brains program. We thank Adam Oberman for early discussions on this idea, and Nan Jiang for helpful conversations. We also thank the anonymous reviewers, Robert Dadashi, and Pablo Samuel Castro for feedback on earlier drafts.





§ LAUNDRY LIST OF CONVERGENT ALGORITHMS


We outline the general proof recipe, which will be re-using for the following examples.



  
Proof strategy
  (P1) Let μ^(1),μ^(2) be initial distributions and (f^(1)_0, f^(2)_0) be the optimal coupling which minimizes (μ^(1),μ^(2));

  (P2) Define an appropriate coupling f^(1)_1 ∼μ^(1) K, f^(2)_1 ∼μ^(2) K– e.g. by defining them to follow the same trajectories if the updates sample from the same distributions;

  (P3) Use the upper bound (μ^(1) K ,μ^(2) K) ≤[f^(1)_1-f^(2)_2] and bound [f^(1)_1-f^(2)_1] ≤ρ[f^(1)_0-f^(2)_0] for some ρ which depends on γ, α, and other parameters of the algorithm. Pick the step-size α such that ρ < 1 to get that μ↦μ K is a contraction.


 §.§ Convergence of synchronous Monte Carlo Evaluation with constant step-sizes


We prove that Monte Carlo Evaluation with synchronous updates & constant step-size converges to a stationary distribution. The algorithm aims to evaluate the value function of a given policy π using Monte Carlo returns. The update rule is given by:

    ∀ s ∈:    V_n+1(s) = (1-α)V_n(s) + α^π_n(s) MCE

where ^π_n(s) = ∑_n ≥ 0γ^n r_n(s_n,a_n) is the return of a random trajectory (s_n,a_n,r_n)_n ≥ 0 starting from s, following a_n ∼π(·|s_n), r_n ∼(·|s_n,a_n), and s_n+1∼(·|s_n,a_n). 

For any constant step size 0<α≤ 1 and initialization V_0 ∼μ_0 ∈(^||), the sequence of random variables (V_n)_n ≥ 0 defined by the recursion (<ref>) converges in distribution to a unique stationary distribution φ_α∈(^||).

Following the proof strategy outlined above, we skip to step (P2) of the proof. We define the coupling of the updates
(V^(1)_1,V^(2)_1) to sample the same trajectories:

    .
    V^(1)_1(s)    = (1-α)V^(1)_0(s) + α^π_k(s) 
    
    V^(2)_1(s)    = (1-α)V^(2)_0(s) + α^π_k(s).
    }for the same^π_k(s)

Note that this is a valid coupling of (μ^(1)K_α,μ^(2)K_α), since V^(1)_1(s) and V^(2)_1(s) have access to the same sampling distributions. 
We upper bound ( μ^(1)K_α,  μ^(2)K_α) by the coupling defined in Equation (<ref>). This gives:

    (μ^(1) K_α,μ^(2) K_α)    ≤𝔼[V^(1)_1-V^(2)_1]  
       = 𝔼[(1-α)V^(1)_0 + α^π_1 - ((1-α)V^(2)_0 + α^π_1 )] 
       = [(1-α)(V^(1)_0-V^(2)_0)] 
       = (1-α)[V^(1)_0-V^(2)_0] = (1-α) (μ^(1),μ^(2))

Since 1-α <1, K_α is a contraction mapping and we are done.


 §.§ Convergence of synchronous Q-Learning with constant step-sizes

We prove that Q-Learning with synchronous updates & constant step-sizes converges to a stationary distribution. The algorithm aims to learn the optimal action-value function Q^⋆. The updates are given by:

    ∀  (s,a) ∈×:     Q_n+1(s,a) = (1-α)Q_n(s,a) + α(r + γmax_a' Q_n(s',a')), QL

where r ∼(·|s,a), s' ∼(·|s,a), and α > 0. 

For any constant step size 0<α≤ 1 and initialization Q_0 ∼μ_0 ∈(^||×||), the sequence of random variables (Q_n)_n ≥ 0 defined by the recursion (<ref>) converges in distribution to a unique stationary distribution ξ_α∈(^||).

We use the proof outline given above, and jump straight to step (P2). We witness the same-sampling coupling again:

    .
    Q^(1)_1(s,a)    = (1-α)Q^(1)_0(s,a) + α(r + γmax_a'Q^(1)_0(s',a')) 
    
    Q^(2)_1(s.a)    = (1-α)Q^(2)_0(s,a) + α(r + γmax_a'Q^(2)_0(s',a'))
    }for the same6.3emr ∼(s,a),
     s' ∼(·|s,a)

The bound follows similarly, but with one additional step. Again we write (Q)(s,a) = r + γmax_a' Q(s'_(s,a),a') for the empirical Bellman (optimality) operator. 

    [(Q^(1)) - (Q^(2))]    =  [max_s,a| r-r + γ(max_a'Q^(1)(s'_(s,a),a') - max_a'Q^(2)(s'_(s,a),a'))|] 
       = γ[max_s,a|max_a'Q^(1)(s'_(s,a),a') - max_a'Q^(2)(s'_(s,a),a') |]
       ≤γ[max_s,amax_a'| Q^(1)(s'_(s,a),a') - Q^(2)(s'_(s,a),a') |]
       ≤γ[ max_s,a| Q^(1)(s,a)-Q^(2)(s,a)|] = γ[ Q^(1)-Q^(2)]

The first inequality follows from |max_a' Q_1(s,a') - max_a'Q_2(s,a')|≤max_a'| Q_1(s,a')-Q_2(s,a')|, and the second inequality follows since Q^(1) and Q^(2) sampled the same s'. Concluding the proof as before we see that the kernel is contractive with Lipschitz constant 1+α-αγ < 1, and we are done.


 §.§ TD(λ)


We prove that TD(λ) with synchronous updates & constant step-size converges to a stationary distribution. The algorithm aims to evaluate the value function of a given policy π using a convex combination of n-step returns. The update rule is given by:

    ∀ s:   V_n+1(s) = (1-α)V_n(s,a) + α(1-λ)∑_k=1^∞λ^k-1(∑_i=0^k γ^i r(s_i,a_i) + γ^k V_n(s_k)) TD(λ)

where each n-step trajectory is sampled starting from s and following policy π.

For any constant step size 0<α≤ 1 and initialization V_0 ∼μ_0 ∈(^||), the sequence of random variables (V_n)_n ≥ 0 defined by the recursion (<ref>) converges in distribution to a unique stationary distribution ζ_α∈(^||).

Again, we jump straight to step (P2) of the template given above. We couple every n-step trajectory to sample the same n rewards, actions, and successors states.

    .
    V^(1)_k+1(s)    = (1-α)V^(1)_k(s) + α(1-λ)∑_n=1^∞λ^n-1(∑_i=0^n-1γ^i r_i(s_i,a_i) + γ^n V^(1)_k(s_n)) 
    
    V^(2)_k+1(s)    = (1-α)V^(2)_k(s) + α(1-λ)∑_n=1^∞λ^n-1(∑_i=0^n-1γ^i r_i(s_i,a_i) + γ^n V^(2)_k(s_n))
    }14.5emsame
    (s_i,a_i,r_i)_i=0^n
    ∀ n


By the coupling, the reward terms will cancel in every n-step trajectory. We write R^(i)_n=∑_i=0^n-1γ^i r_i(s_i,a_i) + γ^n V^(i)_k(s_n) for the n-step return and (V)(s) = ∑_k=1^∞λ^k-1(∑_i=0^k γ^i r(s_i,a_i) + γ^k V_n(s_k)) for the empirical Bellman operator of TD(λ).

    [ (V^(1))-(V^(2))]    =  [max_s | ∑_n=1^∞λ^n-1R^(1)_n- ∑_n=1^∞λ^n-1R^(2)_n| ] 
       = [max_s | ∑_n=1^∞λ^n-1(R^(1)_n - R^(2)_n) | ] 
       =  [max_s | ∑_n=1^∞λ^n-1γ^n(V^(1)(s_n) - V^(2)(s_n)) | ]  reward terms cancel
       ≤[ ∑_n=1^∞λ^n-1γ^n max_s | (V^(1)(s_n) - V^(2)(s_n)) | ] triangle inequality
       ≤∑_n=1^∞λ^n-1γ^n [max_s | V^(1)(s) - V^(2)(s) | ] by the coupling
       =∑_n=1^∞λ^n-1γ^n [V^(1)-V^(2)] = γ1/1-λγ[ V^(1)-V^(2)]
  
Concluding the proof as before, we have (μ^(1)K,μ^(2)K) ≤ (1-α + αγ1-λ/1-λγ )(μ^(1),μ^(2)). Since 1-α + αγ1-λ/1-λγ < 1 we are done. 


 §.§ SARSA with ε-greedy policies


In this example we will example the use of ε-greedy policies for control. In particular, we examine SARSA updates with ε-greedy policies. Let π(·|s) be some base policy. The updates are as follow:


    Q_k+1(s,a) = 
    	(1-α)Q_k(s,a) + α(r(s,a) + γ Q_k(s',a'))     w.p. ε
    
    	(1-α)Q_k(s,a)+α(r(s,a)+γmax_a' Q_k(s',a'))       w.p.  1-ε
    SARSA

where r ∼(·|s,a) and s' ∼(·|s,a) in both cases and a' ∼π(·|s') in the first case.


For any constant step size 0<α≤ 1 and initialization Q_0 ∼μ_0 ∈(^||× ||), the sequence of random variables (Q_n)_n ≥ 0 defined by the recursion (<ref>) converges in distribution to a unique stationary distribution θ_α∈(^||×||).

We jump straight to step (P2) of the proof template. We use the same-sampling coupling, where Q_1^(1) takes the greedy action if and only if Q_1^(2) does. In the non-greedy case, they sample the same a' ∼π(·|s'). In all cases, both functions  sample the same r(s,a) and s'.
We write (Q)(s,a) =  r+ γ Q(s',a')  w.p. ε

r+γmax_a'Q(s',a')  w.p.  1-ε

The bound follows similarly to the examples of Q-learning and TD(0). We omit the subscripts on the Q-functions.

    [(Q^(1)) - (Q^(2))]    = {greedy action chosen}[ max_s,aγ| (max_a'Q^(1)(s',a') - max_a'Q^(2)(s',a')  |] 
          + {non-greedy action chosen}[max_s,a|γ (Q^(1)(s',a')-Q^(2)(s',a'))|] 
       ≤εγ[Q^(1)-Q^(2)] + (1-ε)γ[Q^(1)-Q^(2)] 
       = γ[Q^(1)-Q^(2)]

The bound [ max_s,aγ| (max_a'Q^(1)(s',a') - max_a'Q^(2)(s',a')  |] ≤γ[Q^(1)-Q^(2)] follows from |max_a' Q_1(s,a') - max_a'Q_2(s,a')|≤max_a'| Q_1(s,a')-Q_2(s,a')|, and since Q^(1) and Q^(2) sampled the same s' in the greedy case. The bound [max_s,a|γ (Q^(1)(s',a')-Q^(2)(s',a'))|] ≤[Q^(1)-Q^(2)] follows since Q^(1) and Q^(2) sampled the same state-action pair in the non-greedy case. Concluding the proof as before, we have that [Q_1^(1)-Q_1^(2)]≤ (1-α+αγ) [Q_0^(1)-Q_0^(2)], and thus the kernel is a contraction.


 §.§ Expected SARSA with ε-greedy policies


In this example we examine the Expected SARSA updates with ε-greedy policies. Let π(·|s) be some base policy. Define π_ε(·|s) as the ε-greedy policy which takes the greedy action with probability 1-ε and π otherwise. The updates are as follow:


    Q_k+1(s,a) = (1-α)Q_k(s,a) + α(r(s,a) + γ∑_a'π_ε(a'|s) Q_k(s',a'))
    	Expected-SARSA

where r ∼(·|s,a) and s' ∼(·|s,a) in both cases and a' ∼π(·|s') in the first case.


For any constant step size 0<α≤ 1 and initialization Q_0 ∼μ_0 ∈(^||× ||), the sequence of random variables (Q_n)_n ≥ 0 defined by the recursion (<ref>) converges in distribution to a unique stationary distribution β_α∈(^||×||).

We jump straight to step (P2) of the proof template. We use the same-sampling coupling. 


We write (Q)(s,a) = r+ γ∑_a'π(a'|s) Q(s',a').
The bound follows similarly to the examples of Q-learning and TD(0). We omit the subscripts on the Q-functions.

    [(Q^(1)) - (Q^(2))]    = [ max_s,aγ|∑_a'π_ε(a')Q^(1)(s',a') - ∑_a'π_ε(a')Q^(2)(s',a')  |] 
       ≤[ max_s,aγ∑_a'π_ε(a') | Q^(1)(s',a') - Q^(2)(s',a') |] 
       ≤[ max_s,aγ∑_a'π_ε(a') Q^(1)(s',a') - Q^(2)(s',a')] 
       ≤γ[Q^(1)-Q^(2)]

Concluding the proof as before, we have that [Q_1^(1)-Q_1^(2)]≤ (1-α+αγ) [Q_0^(1)-Q_0^(2)], and thus the kernel is a contraction.


 §.§ Double Q-Learning


In this example we will have to modify our state-space and introduce a new metric on pairs of Q-functions.  The Double Q-Learning algorithm hasselt2010double[This is the original algorithm, not the deep reinforcement learning version given in van2016deep.] maintains two random estimates (Q^A,Q^B) and updates  Q^A with probability p and Q^B with probability 1-p. Should Q^A be chosen to be updated, the update is:

    Q^A_n+1(s,a) = (1-α)Q^A_n(s,a) + α(r(s,a) + γ Q^B_n (s,_a'Q^A_n(s',a'))).

Analogously, the update for Q^B is:

    Q^B_n+1(s,a)=(1-α)Q^B_n(s,a) + α(r(s,a) + γ Q^A_n (s,_a'Q^B_n(s',a'))).

In both cases, we have s' ∼(·|s,a). For this algorithm, the updates are Markovian on pairs of action-value functions. Thus we set the state space to be ^||×||×^||×||. We choose the product metric defined by d_1((Q^A,Q^B),(R^A,R^B)) = Q^A-R^A+Q^B-R^B.


For any constant step size 0<α≤ 1 and initialization (Q^A_0,Q^B_0) ∼μ_0 ∈(^||×||×^||×||), the sequence of random variables (Q^A_n,Q^B_n)_n ≥ 0 defined by the Double Q-Learning recursion converges in distribution to a unique stationary distribution χ_α∈(^||×||×^||×||).

As before, let μ^(1),μ^(2)(^||×||×^||×||) be arbitrary initializations and (Q_0^A,Q_0^B) and (R_0^A,R_0^B) be the optimal coupling of (μ^(1),μ^(2)). We couple (Q_1^A,Q_1^B) and (R_1^A,R_1^B) to sample the same function to be updated and the same s'. Assume for a moment that Q^A and R^A are chosen to be updated. Proceeding as in the proof of Q-Learning (cf. Theorem <ref>), we find that

    [Q_1^A-R_1^A] ≤ (1-α)[Q_0^A-R_0^A] + αγ[Q_0^B-R_0^B].

Analogously, if Q^B and R^B are chosen to updated, we have:

    [Q_1^B-R_1^B] ≤ (1-α)[Q_0^B-R_0^B] + αγ[Q_0^A-R_0^A].

Putting everything together, the full expectation is:

    [d((Q_1^A,Q_1^B),(R_1^A,R_1^B))]   = [Q_1^A-R_1^A+Q_1^B-R_1^B] 
       ={A is updated}[ Q_1^A-R_1^A+Q_1^B-R_1^B] 
          +{B is updated}[ Q_1^A-R_1^A+Q_1^B-R_1^B] 
       =p[ Q_1^A-R_1^A+Q_0^B-R_0^B] 
          +(1-p)[ Q_0^A-R_0^A+Q_1^B-R_1^B] 
       ≤ p((1-α)[Q_0^A-R_0^A] + (1+αγ)[Q_0^B-R_0^B]) 
          +(1-p)((1+αγ)[ Q_0^A-R_0^A] + (1-α)[Q_0^B-R_0^B] ) 
       ≤1/2(2+αγ - α)([ Q_0^A-R_0^A] + [Q_0^B-R_0^B] )  
       = 1/2(2+αγ - α) [ d((Q_0^A,Q_0^B),(R_0^A,R_0^B)) ]

Since 0 ≤ 1/2(2+αγ - α) < 1, so we are done. We note that the first equality only follows since, under the coupling, either A or B is updated for both functions.


§ PROOFS OF SECTION <REF>

Suppose ^π is such that the updates (<ref>) with step-size α converge to a stationary distribution ψ_α. If  is an empirical Bellman operator for some policy π, then [f_α] = f^π where f_α∼ψ_α and f^π is the fixed point of ^π. 

Let f_0 be distributed according to ψ_α. Rewriting equation (<ref>):

    f_1  = (1-α)f_0 + α^π f_0 + αξ (f_0),

where ξ(f_0) = ^π(f_0,ω) - ^π f_0 is a zero-mean noise term. Taking expectations on both sides, and using that f_1 is also distributed according to ψ_α by stationarity and that [ξ(f)] = 0 for any f: 

    f_α   = (1-α) f_α  + α[^π f_0] 
    αf_α   = α[^π + γ^π f_0] 
    f_α   = ^π + γ^π[f_0] 
    f_α   = ^πf_α

And therefore f_α = f^π since it is the unique fixed point of  ^π.

Suppose ^π is such that the updates (<ref>) with step-size α converge to a stationary distribution ψ_α, and that ^π is an empirical Bellman operator for some policy π. Define (f) _ω[(^π(f,ω)-^π f)(^π(f,ω)-^π f)^] to be the covariance of the zero-mean noise term ^π(f,ω)-^π f for a given function f.  Then, the covariance of f_α∼ψ_α is given by 

    (1-(1-α)^2)[(f_α-f^π)(f_α-f^π)^]    = α^2 (γ^π)[(f_α-f^π)(f_α-f^π)^](γ^π)^
          + α(1-α)(γ^π) [(f_0-f^π)(f_0-f^π)^] 
          + α(1-α) [(f_α-f^π)(f_α-f^π)^](γ^π)^
          + α^2 ∫(f)ψ_α(f)


Furthermore, we have that [(f_α-f^π)(f_α-f^π)^]_op is monotonically decreasing with respect to α, where ·_op denotes the operator norm of a matrix. In particular, lim_α→ 0[(f_α-f^π)(f_α-f^π)^]_op = 0, and we have that:  
{min_i | f_α(i) - f^π(i)|≥ε}α→ 0⟶ 0   ∀ ε>0

We preface the proof with some useful identities. We will write the covariance in terms of the tensor product for ease of manipulations

Write ξ(f)  (^π(f,ω)-^π f). In the same setup as Theorem <ref>: 
[(f_α-f^π)(^π f_α - f^π + ξ(f_0))^] = [(f_α-f^π)(f_α-f^π)^](γ^π)^
and 

    [((^π f_α-f^π)+ξ(f_α))((^π f_α-f^π)+ξ(f_α))^]     =  (γ^π)[(f_α-f^π)(f_α-f^π)^](γ^π)^
          + ∫ C(v)ψ_α(v)

Let f_0 ∼ψ_α, by (<ref>) we have f_1 = (1-α)f_0 + α(^π f_0 + ξ(f_0)) and f_1 ∼ψ_α. Furthermore, the distribution of f_0  is independent of the distribution of ω. By independence,  

    [(f_0-f^π)ξ(f_0)^]    = _f_0_ω[(f_0-f^π) ξ(f_0)^] by independence of f_0 and ξ(·)
       = _f_0[(f_0-f^π) (_ωξ(f_0))^] = 0 _ω[ξ(f)] = 0 for every f
 
For the first identity, note that 

    [(f_0-f^π)(^π f_0-f^π))^]    = [(f_0-f^π)(^π + γ^π(f_0)- ^π - γ^π(f^π))^] 
       = [(f_0-f^π)(γ^π(f_0-f^π))^] 
       = [(f_0-f^π)(f_0-f^π)^(γ^π)^] 
       = [(f_0-f^π)(f_0-f^π)^](γ^π)^

The first identity then follows by using [(f_0-f^π)ξ(f_0)^] = 0 and linearity of expectations.

 For the second identity, expanding the outer product gives:

    [((^π f_0-f^π)+ξ(f_0))((^π f_0-f^π)+ξ(f_0))^]    = [(^π f_0-f^π)(^π f_0-f^π)^] 
          + [(ξ(f_0))(ξ(f_0)))^] 
           + [(^π f_0-f^π)(ξ(f_0))^]
          + [ξ(f_0)(^π f_0-f^π)^]
       = [(γ^π(f_0-f^π))(γ^π(f_0-f^π))^]  
          + ∫(v)ψ_α(v) 
       = (γ P^π)[(f_0-f^π)(f_0-f^π)^](γ P^π)^
          + ∫(v)ψ_α(v)

where we used [(^π f_0-f^π)(ξ(f_0))^] = 0. 

Again let f_0 be distributed according to ψ_α. Subtracting f^π  from equation (<ref>), 

    f_1 - f^π = (1-α)(f_0-f^π) + α(^π f_0-f^π+ξ (f_0)).

and taking outer products:

    (f_1 - f^π)(f_1-f^π)^ =    (1-α)^2(f_0-f^π)(f_0-f^π)^
       + α^2(^π f_0-f^π+ξ(f_0))(^π f_0-f^π+ξ(f_0))^
       + α(1-α)(f_0-f^π)(^π f_0-f^π+ξ(f_0))^
       + α(1-α)(^π f_0 - f^π+ξ(f_0))(f_0-f^π)^.

Taking expectations on both sides, and using Lemma <ref>:


    [(f_1 - f^π)(f_1 - f^π)^] =    (1-α)^2[(f_0-f^π)(f_0-f^π)^] + α^2(γ^π)[(f_0-f^π)](γ^π)^
       + α^2∫(v)ψ_a(v) 
       + α(1-α)(γ^π)[(f_0-f^π)(f_0-f^π)^]  
       + α(1-α)[(f_0-f^π)(f_0-f^π)^](γ^π)^

Since [(f_1 - f^π)(f_1 - f^π)^] = [(f_0 - f^π)(f_0 - f^π)^] by stationarity, re-arranging to the LHS and factoring gives: 

    (1-(1-α)^2)[(f_α-f^π)(f_α-f^π)^]    = α^2 (γ^π)[(f_α-f^π)(f_α-f^π)^](γ^π)^
          + α(1-α)(γ^π) [(f_0-f^π)(f_0-f^π)^] 
          + α(1-α) [(f_α-f^π)(f_α-f^π)^](γ^π)^
          + α^2 ∫(f)ψ_α(f)


For the remainder of the proof we re-write the above expression in terms of tensor products. The tensor product of two vectors x,y is the matrix defined by x ⊗ y = xy^. By extension, the tensor product of two matrices A,B is the operator defined by (A ⊗ B)X = AXB^. Then, the above expression can be re-written as:

    (1-(1-α)^2)[(f_α-f^π)(f_α-f^π)^]    = α^2 (γ^π)^⊗ 2[(f_α-f^π)(f_α-f^π)^] 
          + α(1-α)(γ^π⊗) [(f_0-f^π)(f_0-f^π)^] 
          + α(1-α)(⊗γ^π) [(f_α-f^π)(f_α-f^π)^] 
          + α^2 ∫(f)ψ_α(f).

Factoring the tensor products further gives:
[I-((1-α)I + αγ P^π)^⊗ 2][(f_α - f^π)^⊗ 2] = α^2 ∫(f)ψ_α(f)
We show that the matrix on the LHS is invertible. By [Corollary C.4]puterman2014markov it will follow from showing that ρ(((1-α)I + αγ P^π)^⊗ 2) < 1, where ρ(A) is the spectral radius of matrix A. Writing A_op = max_i ∑_j |A(i,j)| for the operator norm of a matrix A, and using that ρ(A) ≤A_op, A ⊗ B_op = A_opB_op, and P^π_op = I_op = 1:

    ((1-α)I + αγ P^π)^⊗ 2_op = (1-α)I + αγ P^π^2_op≤((1-α) + αγ)^2 < 1,

where the last inequality followed since γ < 1.
Finally, for the limit α→ 0,  we use the following identity: if A is such that I-A≤ 1 then A^-1≤1/1-I-A. We let A = I - ((1-α)I + αγ^π)^⊗ 2, by the calculation in (<ref>) we have I-A < 1. So we calculate the operator norm of the covariance matrix:

    [(f_0-f^π)(f_0-f^π)^]   = α^2[I-((1-α)I + αγ P^π)^⊗ 2]^-1∫(v)ψ_α(v)
       ≤α^2 [I-((1-α)I + αγ P^π)^⊗ 2]^-1∫(v)ψ_α(v)
       ≤α^2 1/1-I - I + ((1-α)I + αγ P^π)^⊗ 2∫(v)ψ_α(v)
       = α^2 1/1-((1-α)I + αγ P^π)^⊗ 2∫(v)ψ_α(v)
       = α^2 1/1-((1-α)I + αγ P^π)^2∫(v)ψ_α(v)
       ≤α^2 1/1-(1-α+αγ)^2∫(v)ψ_α(v)


Finally, since the state space is bounded in [0,Rmax/(1-γ)]^n, we have (f)_i ≤Rmax/(1-γ) and ( f)_i ≤Rmax/(1-γ) for each i. Then, we have |ξ_ω(f)_i ξ_ω(f)_j|= |(f)_i( f)_j - ( f)_i(f)_j - ( f)_j(f)_j + ( f)_j ( f)_i| ≤ 4Rmax^2/(1-γ)^2  
Thus we have (f)≤ 4 Rmax^2/(1-γ)^2 M and thus

    [(f_0-f^π)(f_0-f^π)^]   ≤ M α^2/1-(1-α+αγ)^2α→ 0⟶ 0
 


For the concentration inequality, we will use a multivariate Chebyshev inequality [Theorem 3.1]marshall1960multivariate, whos statement is as follows:

Let X = (X_1,...,X_n) be a random vector with X = 0 and [ X^T X] =. Let T = T_+ ∪{x: -x ∈ T_+}, where T_+ ⊆^n is a closed, convex set. If A= { a ∈^n : ⟨ a, x ⟩≥ 1  ∀ x ∈ T_+ }, then 
{ X ∈ T }≤inf_a ∈ A a^ a  


Let ε > 0. We first bound a^ a with the operator norm of . Note that 

    a^ a    = ∑_i a_i ( a)_i 
       ≤∑_i a_i  a≤ n _opa^2

We define T_+ to be the intersection of half-planes the {x | x_i ≥ε}, so that T_+ = { x | x_i ≥ε ∀ i}. Since the half-planes are closed and convex, T_+ is also closed and convex since it is an intersection of closed and convex sets.Then, T = T_+ ∪{x: -x ∈ T_+} = { x | x_i ≥ε ∀ i  or  x_i ≤ -ε ∀ i }. Note that x ∈ T min_i | x_i |≥ε. We define X = f_α - f^π which has zero-mean. Finally, Theorem <ref> states that 
{X ∈ T } = { f_α - f^π∈ T }≤inf_a ∈ A a^T  a ≤ n _opinf_a ∈ Aa^2.
Note that inf_a a^2 is bounded since a = (1/nε,1/nε,....,1/nε) is in A and a^2 = 1/(nε)^2. So n inf_a ∈ Aa^2 ≤ C for some constant C independent of α. From the previous result, we can take the limit of α→ 0 of _op = [ (f_α-f^π)(f_α-f^π)^]_op and obtain:
{ f_α - f^π∈ T } = {min_i | f_α (i) - f^π (i) |≥ε}≤ C ·[ (f_α-f^π)(f_α-f^π)^]_op→ 0 


§ PROOFS OF SECTION <REF>

Suppose π'(s) = _aQ^π(s,a) for each s. Then K(π,π') = {π'  is greedy with respect to ^π} > 0.


We will prove an intermediate probability lemma. Let X_1,...,X_n be mutually independent random variables bounded in [a,b], and F_i(x) = X_i ≤ x denote the cumulative density functions of X_i for i=2,..,n. Note that 

    X_1 ≥ X_2, X_1 ≥ X_3, ..., X_1 ≥ X_n    = ∫_a^b ∫_a^x_1⋯∫_a^x_1(x_1,...,x_n) 
       = ∫_a^b ∫_a^x_1⋯∫_a^x_1_1(x_1)_2(x_2)_n(x_n)   by mutual independence
       = ∫_a^b F_2(x_1) ⋯ F_n(x_1) _1(x_1) 
       = [F_2(X_1)F_3(X_1)⋯ F_n(X_1) ] .
 
Then, we have:

Suppose that [F_i(X_1)] > 0  ∀ i=2,...,n. Then also
[F_2(X_1) ⋯ F_n(X_1) ] >0

It is easy to see that H(x_1) = Π_i = 2^n F_i(x_1) is also a CDF. In particular, H starts at 0, ends at 1, and it monotone and right-continuous. In fact, by Equation (<ref>) it corresponds to the CDF of max(X_2,...,X_n). Assume for a contradiction that [F_2(X_1) ⋯ F_n(X_1) ] =0. By positivity, monotonicity, and right-continuity, we have that H(x_1) = 0  ∀ x_1 ∈ [a,b). Then, for every x we have 
 
H(x) = 0  F_i(x) = 0  for some i.

Since we have H(b)=1 and H(x)=0 otherwise, note that there must exist one i' such that F_i'(b)=1 and F_i'(x)=0 otherwise. If not, then for all i there exists a ε_i > 0 such that F_i(b-ε_i) > 0. By monotonicity, F_i(b-min_i ε_i) > 0  ∀ i, and thus H(b-min_iε_i)>0. Thus we have [F_i'(x)]=0, a contradiction.

Note that

    K(π,π') = π'  is greedy with respect to ^π = for each  s, ^π(s,π'(s)) ≥^π(s,a)  ∀ a .

Fix a state s, write X_i(s) := G^π(s, a_i), and without loss of generality assume that π'(s) = a_1. We first show that [ F_i(X_1)] > 0, i.e. G^π(s, a_1) ≥ G^π(s, a) > 0 for all a. Suppose that it is not so, and pick a such that G^π(s, a_1) ≥ G^π(s, a) = 0. Then

    Q^π(s,a_1)    = [^π(s,a_1)] 
       = {^π(s,a_1) ≥^π(s,a)}[^π(s,a_1) |{^π(s,a_1) ≥^π(s,a)}] 
       + {^π(s,a_1) < ^π(s,a)}[^π(s,a_1) |{^π(s,a_1) < ^π(s,a)}] 
       = 0 + [^π(s,a_1) |{^π(s,a_1) < ^π(s,a)}] 
       < [^π(s,a)] = Q^π(s,a),

which contradicts the fact that π' is greedy wrt Q^π. Hence [ F_i(X_1)] > 0, and we apply Lemma <ref> to this set to conclude that for each s,

    G^π(s, a_1) ≥ G^π(s, a), ∀ a > 0 .

Because the returns are mutually independent, we further know that

    G^π(s, a_1) ≥ G^π(s, a), ∀ s, a = ∏_s ∈𝒮G^π(s, a_1) ≥ G^π(s, a), ∀ a > 0,

completing the proof.


§ ON WASSERSTEIN CONVERGENCE VS. TOTAL VARIATION CONVERGENCE
 
Recall the definition of the Total Variation metric:

The total variation metric between probability measures is defined by:

d_(μ,ν) = sup_∈(^)̣|μ()-ν()|,

for μ,ν∈(^)̣.


Consider a bandit with a single arm that has a deterministic reward of 0. Consider any of the classic algorithms covered in this paper, which will sample a target of 0 at every iteration. It is easy to see that the unique stationary distribution of the algorithm in this instance is a Dirac distribution at 0 (denoted δ_0). 

Suppose a step-size of α<1. If we initialize with some f_0 ≠ 0 then the algorithm will never converge to the true stationary distribution in Total Variation distance. This is because a Dirac distribution at any x≠ 0 is always a constant distance of 1 away from a Dirac at 0. In other words, 

d_(δ_0,δ_f_n) = 1   ∀ n

despite the fact that f_n → 0. On the other hand, we have
(δ_0,δ_f_n) → 0,

since the Wasserstein metric takes into consideration the underlying metric structure of the space.[In particular, the Wasserstein metric isometrically embeds the original metric space into the space of probability measures mardare2018free.]