{'fraction_non_alphanumeric': 0.050507305698196524, 'fraction_numerical': 0.01719133735208752, 'mean_word_length': 3.899963534702808, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characterizing the Exact Behaviors of Temporal Difference Learning Algorithms Using Markov Jump Linear System Theory

In this paper, we provide a unified analysis of temporal difference learning algorithms with linear function approximators by exploiting their connections to Markov jump linear systems (MJLS). We tailor the MJLS theory developed in the control community to characterize the exact behaviors of the first and second order moments of a large family of temporal difference learning algorithms. For both the IID and Markov noise cases, we show that the evolution of some augmented versions of the mean and covariance matrix of the TD estimation error exactly follows the trajectory of a deterministic linear time-invariant (LTI) dynamical system. Applying the well-known LTI system theory, we obtain closed-form expressions for the mean and covariance matrix of the TD estimation error at any time step. We provide a tight matrix spectral radius condition to guarantee the convergence of the covariance matrix of the TD estimation error, and perform a perturbation analysis to characterize the dependence of the TD behaviors on learning rate. For the IID case, we provide an exact formula characterizing how the mean and covariance matrix of the TD estimation error converge to the steady state values. For the Markov case, we use our formulas to explain how the behaviors of TD learning algorithms are affected by learning rate and the underlying Markov chain. For both cases, upper and lower bounds for the mean square TD error are provided. The mean square TD error is shown to converge linearly to an exact limit.


Introduction
Reinforcement learning (RL) has shown great promise in solving sequential decision making tasks [5,48]. One important topic for RL is policy evaluation whose objective is to evaluate the value function of a given policy. A large family of temporal difference (TD) learning methods including standard TD, GTD, TDC, GTD2, DTD, and ATD [47,50,49,38] have been developed to solve the policy evaluation problem. These TD learning algorithms have become important building blocks for RL algorithms. See [17] for a comprehensive survey. Despite the popularity of TD learning, the behaviors of these algorithms have not been fully understood from a theoretical viewpoint. The standard ODE technique [51,9,7,36,8] can only be used to prove asymptotic convergence. Finite sample bounds are challenging to obtain and typically developed in a case-by-case manner. Recently, there have been intensive research activities focusing on establishing finite sample bounds for TD learning methods with linear function approximations under various assumptions. The IID noise case is covered in [16,37,41]. In [6], the analysis is extended for a Markov noise model but an extra projection step in the algorithm is required. Very recently, finite sample bounds for the TD method (without the projection step) under the Markov assumption have been obtained in [45]. The bounds in [45] actually work for any TD learning algorithm that can be modeled by a linear stochastic approximation scheme. It remains unclear how tight these bounds are (especially for the large learning rate region). To complement the existing analysis results and techniques, we propose a general unified analysis framework for TD learning algorithms by borrowing the Markov jump linear system (MJLS) theory [14] from the controls literature. Our approach is inspired by a recent research trend in applying control theory for analysis of optimization algorithms [39,30,31,29,21,52,15,46,28,32,22,3,40,26,4,18,43], and extends the jump system perspective for finite sum optimization methods in [31] to TD learning.
Our key insight is that TD learning algorithms with linear function approximations are essentially just Markov jump linear systems. Notice that a MJLS is described by a linear state space model whose state/input matrices are functions of a jump parameter sampled from a finite state Markov chain. Since the behaviors of MJLS have been well established in the controls field [14,23,1,12,13,33,34,19,20,44], we can borrow the analysis tools there to analyze TD learning algorithms in a more unified manner. Our main contributions are summarized as follows.
1. We present a unified Markov jump linear system perspective on a large family of TD learning algorithms including TD, TDC, GTD, GTD2, ATD, and DTD. Specifically, we make the key observation that these methods are just MJLS subject to some prescribed input. 2. By tailoring the existing MJLS theory, we show that the evolution of some augmented versions of the mean and covariance matrix of the estimation error in all above TD learning methods exactly follows the trajectory of a deterministic linear time-invariant (LTI) dynamical system for both the IID and Markov noise cases. As a result, we obtain unified closed-form formulas for the mean and covariance matrix of the TD estimation error at any time step. 3. We provide a tight matrix spectral radius condition to guarantee the convergence of the covariance matrix of the TD estimation error under the general Markov assumption. By using the matrix perturbation theory [42,35,2,24], we perform a perturbation analysis to show the dependence of the behaviors of TD learning on learning rate in a more explicit manner. For the IID case, we provide an exact formula characterizing how the mean and covariance matrix of the TD estimation error converge to the steady state values at a linear rate. For the Markov case, we use our formulas to explain how the behaviors of TD learning algorithms are affected by learning rate and the underlying Markov chain. For both cases, we have shown that the mean square error of TD learning converges linearly to an exact limit. In addition, upper and lower bounds for the mean square error of TD learning are simultaneously derived.
We view our proposed analysis as a complement rather than a replacement for existing analysis techniques. The existing analysis focuses on upper bounds for the TD estimation error. Our closedform formulas provide both upper and lower bounds for the mean square error of TD learning. Our analysis also sheds light on the exact limit of the final TD error and the related convergence rate.

Notation
The set of m-dimensional real vectors is denoted as R m . The Kronecker product of two matrices A and B is denoted by when the matrices have compatible dimensions. Let vec denote the standard vectorization operation that stacks the columns of a matrix into a vector. We have vec(AXB) = (B T ⊗ A) vec(X). Let sym denote the symmetrization operation, i.e. sym(A) = A T +A 2 . Let diag(H i ) denote a matrix whose (i, i)-th block is H i and all other blocks are zero. Given H i for i = 1, . . . , n, we have A square matrix is Schur stable if all its eigenvalues have magnitude strictly less than 1. A square matrix is Hurwitz if all its eigenvalues have strictly negative real parts. The spectral radius of a matrix H is denoted as σ(H). The eigenvalue with the largest magnitude of H is denoted as λ max (H) and the eigenvalue with the largest real part of H is denoted as λ max real (H).

Linear time-invariant systems
A linear time-invariant (LTI) system is typically governed by the following state-space model where x k ∈ R nx , u k ∈ R nu , H ∈ R nx×nx , and G ∈ R nx×nu . The LTI system theory has been well documented in standard control textbooks [27,10]. Here we briefly review several useful results.
• Closed-form formulas for x k : Given an initial condition x 0 and an input sequence {u k }, the sequence {x k } can be determined using the following closed-form expression where (H) k stands for the k-th power of the matrix H. • Necessary and sufficient stability condition: When H is Schur stable, we know (H) k x 0 → 0 for any arbitrary x 0 . When σ(H) ≥ 1, there always exists x 0 such that (H) k x 0 does not converge to 0. When σ(H) > 1, there even exists x 0 such that (H) k x 0 → ∞. See Section 7.2 in [27] for a detailed discussion. A well-known result in the controls literature is that the LTI system (1) is stable if and only if H is Schur stable. • Exact limit for x k : If H is Schur stable and u k converges to a limit u ∞ , then x k will converge to an exact limit. This is formalized as follows. Proposition 1. Consider the LTI system (1). If σ(H) < 1 and lim k→∞ u k = u ∞ , then lim k→∞ x k exists and we have • Response for constant input: If u k = u ∀k and σ(H) < 1, then the closed-form expression for x k can be further simplified to give the following tight convergence rate result. Proposition 2. Suppose σ(H) < 1, and x k is determined by (1). If u k = u ∀k, then x k converges to a limit point In addition, x k − x ∞ ≤ C 0 (σ(H) + ε) k for some C 0 and any arbitrarily small ε > 0.
From the above proposition, we can clearly see that now x k is a sum of a constant steady state term x ∞ and a matrix power term that decays at a linear rate specified by σ(H) (see Section 2.2 in [39] for more explanations). The convergence rate characterized by (σ(H) + ε) is tight. More discussions on the tightness of this convergence rate are provided in the supplementary material. • Response for exponentially shrinking input: When u k itself converges at a linear ratẽ ρ and H is Schur stable, x k will converge to its limit point at a linear rate specified by max{σ(H) + ε,ρ}. A formal statement is provided as follows. Proposition 3. Suppose σ(H) < 1, and x k is determined by (1). If u k converges to u ∞ as u k − u ∞ ≤ Cρ k , then we have x ∞ = lim k→∞ x k = (I − H) −1 Gu ∞ and x k − x ∞ ≤ C 0 (max{σ(H) + ε,ρ}) k for some C 0 and any arbitrarily small ε > 0.
The results in Propositions 1, 2, and 3 are well known in the control community. For completeness, we will include their proofs in the supplementary material.

Markov jump linear systems
Another important class of dynamic systems that have been extensively studied in the controls literature is the so-called Markov jump linear system (MJLS) [14]. Let {z k } be a Markov chain sampled from a finite state space S. A MJLS is governed by the following state-space model: where H(z k ) and G(z k ) are matrix functions of z k . Here, ξ k is the state, and y k is the input. There is a one-to-one mapping from S to the set N := {1, 2, . . . , n} where n = |S|. We can assume H(z k ) is sampled from a set of matrices {H 1 , H 2 , . . . , H n } and G(z k ) is sampled from {G 1 , G 2 , . . . , G n }.
We have H(z k ) = H i and G(z k ) = G i when z k = i. The MJLS theory has been well developed in the controls community [14]. We will apply the MJLS theory to analyze TD learning algorithms.

A general Markov jump system perspective for TD learning
In this section, we provide a general jump system perspective for TD learning with linear function approximations. Notice that many TD learning algorithms including TD, TDC, GTD, GTD2, A-TD, and D-TD can be modeled by the following linear stochastic recursion: where {z k } forms a finite state Markov chain and b(z k ) satisfies lim k→∞ Eb(z k ) = 0. 1 We have A(z k ) = A i and b(z k ) = b i when z k = i. For simplicity, we mainly focus on analyzing (5). Other models including two time-scale schemes [25,54] will be discussed in the supplementary material.
The above model is a special case of (4) if we set H(z k ) = I + αA(z k ), G(z k ) = αb(z k ), and y k = 1 ∀k. Consequently, many TD learning algorithms can be analyzed using the MJLS theory.
We briefly highlight the analysis idea here. Our analysis is built upon the fact that some augmented versions of the mean and the covariance matrix of {ξ k } for the MJLS model (4) actually follow the dynamics of a deterministic LTI model in the form of (1) [14,Chapter 3]. To see this, we denote the transition probabilities for the Markov chain {z k } as p ij := P(z k+1 = j|z k = i) and specify the transition matrix P by setting its (i, j)-th entry to be p ij . Obviously, we have p ij ≥ 0 and n j=1 p ij = 1 for all i. Next, the indicator function 1 {z k =i} is defined as 1 {z k =i} = 1 if z k = i and 1 {z k =i} = 0 otherwise. Now we define the following key quantities: [14,Proposition 3.35], q k and Q k can be iteratively calculated as where p k i := P(z k = i). If we further augment q k i and Q k i as then it is straightforward to rewrite (7) (8) as the following LTI system q k+1 vec(Q k+1 ) = where H 11 , H 21 , H 22 , u k q , and u k Q are given by A detailed derivation for the above result is presented in the supplementary material. A key implication here is that q k and vec(Q k ) follow the LTI dynamics (9) and can be analyzed using the standard LTI theory reviewed in Section 2.2. Obviously, we have Hence the mean, covariance, and mean square norm of ξ k can all be calculated using closed-form expressions. We will present a detailed analysis for (6) and provide related implications for TD learning in the next two sections.
For illustrative purposes, we explain the jump system perspective for the standard TD method.
Example 1: TD method. The standard TD method (or TD(0)) uses the following update rule: where {s k } is the underlying Markov chain, φ is the feature vector, r is the reward, γ is the discounting factor, and θ k is the weight vector to be estimated. Suppose θ * is the vector that solves the projected Bellman equation. We can set z k = (s k+1 ) T (s k ) T T and then rewrite the TD update as where Since the projected Bellman equation and the equation n i=1 p ∞ i b i = 0 are actually equivalent, we have naturally enforced lim k→∞ Eb(z k ) = 0. Therefore, the TD update can be modeled as (6) with b(z k ) satisfying lim k→∞ Eb(z k ) = 0. See Section 3.1 in [45] for a similar formulation. Now we can apply the MJLS theory and the LTI model (9) to analyze the covariance E (θ k − θ * )(θ k − θ * ) T and the mean square error E θ k − θ * 2 . In this case, we have Then we can easily analyze q k and Q k by applying the LTI model (9). In general, the covariance matrix E (θ k − θ * )(θ k − θ * ) T and the mean value E(θ k −θ * ) do not directly follow an LTI system. However, when working with the augmented covariance matrix Q k and the augmented mean value vector q k , we do obtain an LTI model in the form of (1). Once the closed-form expression for Q k is obtained, the mean square estimation error for the TD update can be immediately calculated as Here we omit the detailed formulations for other TD learning methods. The key message is that {z k } can be viewed as a jump parameter and TD learning methods are essentially just MJLS. Notice that all the TD learning algorithms that can be analyzed using the ODE method are in the form of (6). Jump system perspectives for other TD learning algorithms are discussed in the supplementary material.
In this paper, we will assumeĀ is Hurwitz. This assumption is standard and even required by the basic ODE approach. For the standard TD method,Ā is Hurwitz when the discount factor is smaller than 1, p ∞ i is positive for all i, and the feature matrix is full column rank. It is worth emphasizing that the assumption This case can still be handled using a shifting argument sinceĀ is Hurwitz. Notice the iteration ξ k+1 = (I +αA(z k ))ξ k +αb(z k ) can be (4)). Notice that (4) provides a general jump system model for linear stochastic schemes that can be in more complicated forms than (5). However, (4) can not be directly used to cover nonlinear stochastic approximation schemes. See [53,11] for recent finite sample results on nonlinear stochastic approximation over non-IID data.

Analysis under the IID assumption
For illustrative purposes, we first present the analysis for (6) under the IID assumption (P(z k = i) = p i ∀i). In this case, the analysis is significantly simpler, since {Eξ k } and {E ξ k (ξ k ) T } directly form LTI systems with much smaller dimensions. We denote µ k := Eξ k and Q k := E ξ k (ξ k ) T . Then the following equations hold for the general jump system model (4) There are many ways to derive the above formulas. One way is to first show q k i = p i µ k and Q k i = p i Q k in this case and then apply (7) and (8). Another way is to directly modify the proof of Theorem 1. Now consider the jump system model (6) We can directly obtain the following result. Theorem 1. Consider the jump system model (6) with Then µ k and vec(Q k ) are governed by the following LTI system: where H 11 , H 21 and H 22 are determined as In addition, if σ(H 22 ) < 1, we have where µ ∞ = lim k→∞ µ k = 0, and vec(Q ∞ ) is given as Proof. For completeness, a detailed proof is presented in the supplementary material.
Now we discuss the implications of the above theorem for TD learning. For simplicity, we denote Stability condition for TD learning. From the LTI theory, the system (14) is stable if and only if H is Schur stable. We can apply Proposition 3.6 in [14] to show that H is Schur stable if and only if H 22 is Schur stable. Hence, a necessary and sufficient stability condition for the LTI system (14) is that H 22 is Schur stable. Under this condition, the first term on the right side of (16) converges to 0 at a linear rate specified by σ(H), and the second term on the right side of (16) is a constant matrix quantifying the steady state covariance. An important question for TD learning is how to choose α such that σ(H 22 ) < 1 for some given {A i }, {b i }, and {p i }. We provide some clue to this question by applying an eigenvalue perturbation analysis to the matrix H 22 . We assume α is small. Then under mild technical condition 2 , we can ignore the quadratic term Then we immediately obtain σ(H 22 ) ≈ 1 + 2 real(λ max real (Ā))α + O(α 2 ). Therefore, as long asĀ is Hurwitz, there exists sufficiently small α such that σ(H 22 ) < 1. More details of the perturbation analysis are provided in the supplementary material.
Exact limit for the mean square error of TD learning. Obviously, µ k converges to 0 at the rate specified by σ(I + αĀ) due to the relation µ k = (I + αĀ) k µ 0 . Applying Proposition 3 and making use of the block structure in H, one can show vec( , which leads to the result in (17). A key message here is that the covariance matrix converges linearly to an exact limit under the stability condition σ(H 22 ) < 1. We can clearly see lim k→0 vec(Q k ) = O(α) and can be controlled by decreasing α. When α is large, we need to keep the quadratic term . Therefore, our theory captures the steady-state behavior of TD learning with both small and large α, and complement the existing finite sample bounds in literatures. To further compare our results with existing finite sample bounds, we obtain the following result for the mean square error of TD learning.
where Q ∞ is given by (17). In addition, the following mean square TD error bounds hold for some constant C 0 and any arbitrary small positive ε: Finally, for sufficiently small α, one has is a semisimple eigenvalue, then σ(H) = σ(H 11 ) = 1 + real(λ max real (Ā))α.
Proof. Recall that we have E θ k − θ * 2 = trace(Q k ). Taking limits on both sides leads to the expression for δ ∞ . Then we can apply Proposition 2 to obtain a linear convergence bound for Q k which eventually leads to (19). Finally, we can apply standard matrix perturbation theory to show δ ∞ = O(α) and σ(H) = σ(H 11 ) = 1 + real(λ max real (Ā))α.
The above corollary gives both upper and lower bounds for the mean square error of TD learning. From the above result, the final TD estimation error is actually exactly on the order of O(α). This justifies the tightness of the existing upper bounds for the final TD error up to a constant factor. From the above corollary, we can also see that one can increase the convergence rate at the price of increasing the steady state error. This is consistent with the finite sample bound in the literature [6,45]. It is also worth mentioning that (σ(H) + ε) is a tight estimation for the convergence rate of the matrix power (H) k . See the supplementary material for more explanations.

Analysis under the Markov assumption
Now we can analyze the behaviors of TD learning under the general assumption that {z k } is a Markov chain. Recall that the augmented mean vector q k and the augmented covariance matrix Q k have been defined in Section 3. We can directly modify (9) to obtain the following result.
Then q k and vec(Q k ) are governed by the following state-space model: where H 11 , H 21 and H 22 are given by H 11 = (P T ⊗ I n ξ ) diag(I n ξ + αA i ), In addition, the following closed-form solution holds for any k where H 11 , H 21 and H 22 are determined by (21).

Proof.
A detailed proof is presented in the supplementary material. We present a proof sketch here. Notice (20) is a direct consequence of (7) and (8) (which are special cases of Proposition 3.35 in [14]). Specifically, one can verify the following equations using the Markov assumption Then we can apply the basic property of the vectorization operation vec to obtain (20). Applying (2) to iterate (20) directly leads to (22).
Therefore, the evolutions of q k and Q k can be fully understood via the well-established LTI system theory. Now we discuss the implications of Theorem 2 for TD learning.
Stability condition for TD learning. Similar to the IID case, the necessary and sufficient stability condition is σ(H 22 ) < 1. Now H 22 becomes a much larger matrix depending on the transition matrix P . An important question is how to choose α such that σ(H 22 ) < 1 for some given {A i }, {b i }, P , and {p 0 }. Again, we perform an eigenvalue perturbation analysis for the matrix H 22 . This case is quite subtle due to the fact that we are no longer perturbing an identity matrix. We are perturbing the matrix (P T ⊗ I n 2 ξ ) and the eigenvalues here are not simple. Under the ergodicity assumption, the largest eigenvalue for (P T ⊗ I n 2 ξ ) (which is 1) is semisimple. Hence we can directly apply the results in Section II of [35] or Theorem 2.1 in [42] and p ∞ is the unique stationary distribution of the Markov chain under the ergodicity assumption. Then we still have σ(H 22 ) ≈ 1 + 2 real(λ max real (Ā))α + o(α). Therefore, as long asĀ is Hurwitz, there exists sufficiently small α such that σ(H 22 ) < 1. This is consistent with Assumption 3 in [45]. To understand the details of our perturbation argument, we refer the readers to the remark placed after Theorem 2.1 in [42]. Notice we have . The largest eigenvalue of P T ⊗ I n 2 ξ is semisimple due to the ergodicity assumption. Then the perturbation result directly follows as a consequence of Theorem 2.1 in [42]. More explanations are also provided in the supplementary material.
Exact limit for the mean square TD error and related convergence rate. Assume the Markov chain is geometrically ergodic, and then p t → p ∞ at some linear rate where p ∞ is the stationary distribution. In this case, we can apply Proposition 3 to show that the mean square error of TD learning converges linearly to an exact limit. Corollary 2. Consider the TD update (12) withĀ being Hurwitz. Let {z k } be a Markov chain sampled from N using the transition matrix P . Suppose σ(H 22 ) < 1. We set N = nn 2 ξ . If we assume p k → p ∞ where p ∞ is the stationary distribution for {z k }, then we have If we further assume the geometric ergodicity, i.e. p k − p ∞ ≤ Cρ k , then we have where C 0 is some constant and ε is an arbitrary small positive number. For sufficiently small α, we have δ ∞ = O(α) and σ(H) = 1 + real(λ max real (Ā))α + o(α).
Proof. Notice E θ k − θ * 2 = (1 T n ⊗vec(I n ξ ) T ) vec(Q k ). We can directly apply Theorem 2, Proposition 1, and Proposition 3 to prove (26) and (27). When α is small, we can apply the Laurent series trick in [2,24] to show that lim k→0 vec(Q k ) = O(α) and δ ∞ = O(α). The difficulty here is that ξ is a singular matrix and hence (I N − H 22 ) −1 does not have a Taylor series around α = 0. Therefore, we need to apply some recent matrix inverse perturbation result to perform a Laurent expansion of (I N − H 22 ) −1 . From the ergodicity assumption, we know the singularity order of I N − P T ⊗ I n 2 ξ is just 1. Applying Theorem 1 in [2], we can obtain the Laurent expansion of (I Due to the assumption n i=1 p ∞ i b i = 0, we have lim k→∞ q k = 0 in general but µ ∞ = 0. Again, we have obtained both upper and lower bounds for the mean square TD error. From the above result, the final TD error is actually exactly on the order of O(α). This justifies the tightness of the existing upper bounds for the final TD error [6,45] up to a constant factor. From the above corollary, we can also see the trade-off between the convergence rate and the steady state error. Clearly, the convergence rate in (27) also depends on the initial distribution p 0 and the mixing rate of the underlying Markov jump parameter {z k } (which is denoted asρ). If the initial distribution is the stationary distribution, i.e. p 0 = p ∞ , the input to the LTI dynamical system (20) is just a constant for all k and then we will be able to obtain an exact formula similar to (16). However, for a general initial distribution p 0 , the mixing rateρ matters more and may affect the overall convergence rate. One resultant guideline for algorithm design is that increasing α may not increase the convergence rate when the mixing rateρ dominates the convergence process. When α becomes smaller and smaller, eventually σ(H) is going to become the dominating term and the mixing rate does not affect the convergence rate any more.
Algorithm design. Here we make a remark on how our proposed MJLS framework can be further extended to provide clues for designing fast TD learning. When α (or even other hyperparameters including momentum term) is changing with time, we can still obtain expressions of vec(Q k ) and q k in an iterative form. However, both H and G depend on k now. Then given a fixed time budget T , in theory it is possible to minimize the mean square estimation error at T subject to some optimization constraints in the form of a time-varying iteration q k+1 vec(Q k+1 ) = H(k) q k vec(Q k ) + G(k)u k . One may use this control-oriented optimization formulation to gain some theoretical insights on how to choose hyperparameters adaptively for fast TD learning. Clearly, solving such an optimization problem requires knowing the underlying Markov model. However, this type of theoretical study may lead to new hyperparameter tuning heuristics that do not require the model information.

Supplementary Material
A More facts about LTI systems A.1 Tightness of the spectral radius stability condition Technically speaking, the condition σ(H) < 1 is necessary and sufficient for the asymptotic stability of the LTI system (1). Since exponential stability and asymptotic stability are equivalent notions of stability for LTI systems, the condition σ(H) < 1 is also necessary and sufficient for the exponential stability of (1). See Theorem 8.3 and Theorem 8.4 in [27] for formal statements of these facts. We will give a more intuitive explanation here. Specifically, we look at the behaviors of the matrix power term (H) k x 0 . This is the homogeneous state response of (1). There are three possible behaviors for this term.
1. When H is Schur stable (or equivalently σ(H) < 1), the term (H) k converges to a zero matrix and (H) k x 0 → 0 for any arbitrary x 0 . The convergence rate is linear and is completely specified by the spectral radius σ(H). See Section 7.2 in [27] for a detailed explanation of the above fact. Consequently, we can only guarantee (H) k x 0 to converge for all x 0 when σ(H) < 1. Therefore, the condition σ(H) < 1 is a tight condition for the stability of the LTI system (1). As mentioned above, when σ(H) < 1, (H) k x 0 converges at a linear rate completely determined by σ(H). Technically speaking, the convergence rate is either equal to σ(H) + ε for some arbitrary small positive ε or just equal to σ(H) itself. Now we provide a detailed discussion on this convergence rate.

A.2 Convergence rate of the matrix power
Notice we have (H) k x 0 = ρ k ρ −1 H k x 0 . As long as ρ −1 H k x 0 stays bounded for any x 0 , the term (H) k x 0 will converge at the linear rate ρ.
From Section 7.2 in [27], we can almost directly see how to determine the convergence rate of (H) k x 0 .  In this paper, for simplicity we do not want to further look at the Jordan decomposition structure of H and hence we always set the rate as ρ = σ(H) + ε. One can also use the relationship between spectral radius and other matrix norms to obtain the above convergence rate. See Section 2.2 in [39] for a detailed discussion. The arbitrarily small number ε also appears in the argument there.

B More discussions about Markov jump linear systems
First, we verify that (7) and (8) are equivalent to the LTI model (9). Rewriting (7) as an LTI model is quite trivial. Rewriting (8) as an LTI system requires applying the vectorization operation to obtain the following formula, Then we can augment the update rules for q k and vec(Q k ) to obtain the desired LTI model for (q k , vec(Q k )).
Now we briefly review how to analyze q k and Q k using the LTI model (9). We denote H = H11 0 H21 H22 .
First, we will have the following closed-form formula for computing (q k , vec(Q k )): (B. 2) The first term on the right side of the above equation will be guaranteed to converge to 0 if we have the stability condition σ(H) < Finally, if p k − p ∞ ≤ Cpρ k , then there exists a constant C such that the following inequality holds.
If we know σ(H22) < 1, then we can directly apply Proposition 3 to obtain the following linear convergence result: where C0 is some constant and ε is an arbitrarily small positive number. We can see that the convergence rates of vec(Q k ) and q k depend on both σ(H) and the mixing rate of the underlying Markov jump parameter {z k } (which is denoted asρ).
Therefore, when the underlying Markov chain {z k } is aperiodic and irreducible, the mean and covariance of the MJLS just converges to the steady state values at a linear rate specified by max{σ(H) + ε,ρ}. This is a powerful result that can be potentially applied to more general stochastic approximation schemes other than (5).
We also want to mention that there are other proofs for the convergence of (q k , Q k ). See Proposition 3.36 in [14] for an alternative proof. Here, our result is a little bit stronger than Proposition 3.36 in [14] since we also specify the convergence rate of (q k , Q k ).
Finally, it is worth mentioning that under the assumption σ(H22) < 1, one can further prove {Q k } converges to a stationary process in some sense. This is a stronger result. Specifically, Proposition 3.37 In [14] shows that the MJLS is "asymptotically wide sense stationary" under the assumption σ(H22) < 1. We are not that interested in the correlation between the updates at different steps since our main purpose is analyzing TD learning. Hence we will skip a detailed discussion on this topic. See Chapter 3.4 in [14] for a thorough treatment.
C Detailed proofs C.1 Detailed proofs of Propositions 1, 2, and 3 We believe that all the statements in Propositions 1, 2, and 3 are known in the controls field. Since we are not able to find a reference to exactly match the statements, we provide a proof here for completeness.
We will need the following lemma.
Lemma C.1. Consider the LTI model (1). Suppose σ(H) < 1. We set ρ = σ(H) + ε where ε is an arbitrary small positive number. Then there exists a positive definite matrix V and a positive constant γ s.t. the following inequality holds for all k, Proof. We know ρ −1 H is Schur stable. Based on Theorem 8.4 in [27], there exists a positive definite matrix V such that where the matrix inequality holds in the negative definite sense. The above condition is actually equivalent to By Schur complement lemma, this is equivalent to Now we left and right multiply the right side of the above matrix inequality with (x k ) T (u k ) T and (x k ) T (u k ) T T . This leads to One can verify that the first term on the left side of the above inequality is just equal to ( Now we are ready to prove Proposition 1. Since σ(H) < 1, x ∞ can still be well defined as x ∞ = (I −H) −1 Gu. Notice we have not shown the existence of lim k→∞ x k at this point. We will show lim k→∞ x k exists and is equal to x ∞ . Applying the relation (I − H)x ∞ = Gu ∞ , we still have By Lemma C.1, there exists a positive definite matrix V and a positive number γ such that is bounded for all k and then we apply lim sup to the above inequality. Clearly there exists a constant U such that u k − u ∞ 2 ≤ U for all k. Then we have is finite. Now we take lim sup on both sides of (C.2) and will immediately be able to show lim sup k→∞ ( which directly leads to the desired conclusion Finally, we will still use (C.2) to prove Proposition 3. It is assumed that the arbitrary small ε is chosen in a way that ε + σ(H) =ρ since one can always decrease ε by a tiny bit. Then iterating (C.2) leads to Obviously the right side of the above inequality is on the order of O (max{ρ,ρ}) 2k . Hence we have This completes the proof of this proposition.
An interesting thing is that when ρ =ρ, the convergence rate is actually on the order of O(kρ k ). Specifically, we have Of course this rate is always bounded above by O((ε + ρ) k ). In addition, if it happens ε + σ(H) =ρ, one can always decrease ε by a tiny bit and the convergence rate becomes linear again.

C.2 A detailed proof for Theorem 1
The underlying probability space is denoted by (Ω, F, P). We denote by F k the σ-algebra generated by (z 0 , z 1 , . . . , z k ). Clearly, z k is F k -adapted and we obtain a filtered probability space (Ω, F, {F k }, P) on which the stochastic optimization method is defined.
Next, we prove We can use a similar argument. We have Taking full expectation and applying the fact n i=1 pibi = 0 leads to This proves the recursive formula for Q k . Now we can apply the vectorization operation to this formula. For any matrices A, X, and B, we have vec(AXB) = (B T ⊗ A) vec(X). Hence we can directly show vec(Q k+1 ) = vec(Q k ) + α(vec(ĀQ k ) + vec(Q kĀT )) + α 2 Therefore, we have vec(Q k+1 ) = H22 vec(Q k ) + H21µ k + α 2 n i=1 pi(bi ⊗ bi) where H21 and H22 are determined by (15). Putting this together with µ k+1 = (I + αĀ)µ k gives us the LTI model in (14). Then notice we have .
Recall that we have H = H11 0 H21 H22 . We can apply Proposition 3.6 in [14] to show that H is Schur stable if and only if H22 is Schur stable. Therefore, a direct application of Proposition 3 will lead to (16). This completes the proof for Theorem 1.

C.3 A detailed proof for Theorem 2
One may prove this theorem as a corollary of Proposition 3.35 in [14]. For completeness, we add more detailed calculations and present the proof in a self-contained manner. One can update q k+1 j as q k+1 which can be compactly written as q k+1 = (P T ⊗ I) diag(I + αAi)q k + α( This proves q k+1 = H11q k + α((P T diag(p k i )) ⊗ In ξ )b, where H11 is given by (21).
Next, we perform similar steps to obtain the iterative formula for Q k . One can update Q k+1 We can compactly rewrite the above equation as vec(Q k+1 ) = H22 vec(Q k ) + H21q k + α 2 diag(p k i )(P T ⊗ I n 2 ξ )B, where H22 and H21 are given by (21), andB = (b1 ⊗ b1) T . . . (bn ⊗ bn) T T . Putting the recursion formulas for q k and vec(Q k ) together leads to the desired state-space model (20). The rest of the theorem statement follows from direct applications of Equation (2).

D Details for perturbation analysis under the Markov assumption
The perturbation analysis in Section 5 relies on a few technical lemmas from matrix perturbation theory. We provide more details here. We will use the following fact. Proposition D.1. Suppose λ is a simple eigenvalue of A with left eigenvector y and right eigenvector x. Suppose B and A ⊗ Im have the same dimension. Let c be an eigenvalue of the m × m matrix (y ⊗ Im)B(x ⊗ Im). Then the matrix (A ⊗ Im) + αB has an eigenvalue yielding the first-order expansion λ + cα + O(α 2 ) for small α.
The above proposition is a special case of Theorem 2.1 in [42]. See the remark placed after Theorem 2.1 in [42] for explanations. Now we can directly apply the above proposition to analyze the spectral radius of H11 and H22. First recall that H11 = (P T ⊗ In ξ ) diag(In ξ + αAi) = P T ⊗ In ξ + α(P T ⊗ In ξ ) diag(Ai). Based on the ergodicity assumption on {z k }, 1 is a simple eigenvalue of P T with left eigenvector y = 1 1 . . . 1 and right eigenvector p ∞ which is the unique stationary distribution of {z k }. Since we have (y ⊗ In ξ )(P T ⊗ In ξ ) diag(Ai)(p ∞ ⊗ In ξ ) = n i=1 p ∞ i Ai =Ā, we can directly apply the above proposition to show λmax(H11) = 1 + λ max real (Ā)α + O(α 2 ).

E Connections to existing finite sample bounds on mean square errors
Most existing finite sample bounds for TD learning with a constant learning rate have the following form: where C0 is a constant, ρ 2 is the convergence rate, and C1 quantifies the final error level. Typically one proves ρ 2 = 1 − cα + O(α 2 ) for some c, and C1 = O(α). One most relevant result of this nature for the Markov noise model was presented as Theorem 7 in [45]. Our result justifies the tightness of the result in [45] for the following reasons.
• In [45], the constant c in the rate ρ 2 is a constant determined byĀ. Using our perturbation analysis, we can see eventually c is going to be determined by the real part of λ max real (Ā). Actually one could modify the argument in [45] to match the constant c with our perturbation analysis result by choosing a slightly better Lyapunov function based on the solution of an linear matrix inequalitȳ A T V + VĀ + 2ρV 0.
• In [45], the constant C1 is at the order of O(α) which matches the perturbation analysis result obtained in our paper. It is possible that the constant C1 may be improved to match the steady state mean square error lim k→∞ trace(Q k ) obtained by our perturbation analysis, although we have not pursued such an analysis. • In [45], the rate ρ does not depend on the mixing time property. This is consistent with our theory.
Based on our theory, as α gets smaller, the rate ρ becomes independent of the mixing rateρ, although the constant C0 still has some dependence onρ.
• Our results further show that E ξ k 2 converges to an exact limit at a linear rate. Hence both upper and lower bounds for the mean square TD error are simultaneously provided.
It is worth mentioning that the bounds in the form of (E.1) capture the behaviors of TD learning quite well for small α, but can be conservative for large α. Our formulas are exact for all α. The generality comes at the price of loosing some interpretability for the large learning rate region. How to interpret σ(H22) for larger α in a better way remains unclear at this moment.
F More discussions on jump system formulations for variants of TD(0) Finally, we present some extra details for the jump system formulations of several TD learning algorithms other than TD(0). Specifically, all the methods that can be analyzed using the ODE method has the form ξ k+1 = (I + αA(z k ))ξ k + αb(z k ). Then taking expectation of A(z k ) and b(z k ) under the stationary distribution and making α arbitrarily small leads to the ODEξ =Āξ. As commented in Section 3, the linear stochastic approximation scheme ξ k+1 = (I + αA(z k ))ξ k + αb(z k ) is just a MJLS. Now we give detailed references for this type of formulations for various TD learning algorithms. The detailed linear stochastic approximation form for GTD is given in Section 4 of [50]. The detailed linear stochastic approximation form for GTD2 is given in Section 5 of [49]. TDC yields a similar formulation. The double temporal difference (DTD) learning method and the average temporal difference (ATD) learning method are proposed in [38]. The ODE formulations for both DTD and ATD are presented in the supplementary material of [38], yielding straightforward jump system formulations.
It is also possible to model two time-scale methods or off-policy TD learning using the general jump system model (4). One needs to slightly modify the definitions of {Hi} and {Gi}. Then one can immediately apply the LTI model (9) to obtain closed-form formulas for the mean square error of these methods.