{'fraction_non_alphanumeric': 0.048251688950448034, 'fraction_numerical': 0.020239138423269357, 'mean_word_length': 4.100634371395617, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Temporal-Differential Learning in Continuous Environments

In this paper, a new reinforcement learning (RL) method known as the method of temporal differential is introduced. Compared to the traditional temporal-difference learning method, it plays a crucial role in developing novel RL techniques for continuous environments. In particular, the continuous-time least squares policy evaluation (CT-LSPE) and the continuous-time temporal-differential (CT-TD) learning methods are developed. Both theoretical and empirical evidences are provided to demonstrate the effectiveness of the proposed temporal-differential learning methodology.


Introduction
Over the past decade, reinforcement learning (RL) has quickly become one of the most prominent technologies in the field of artificial intelligence (AI). The popularity of RL is mainly driven by its recent success in various application fields (Barto et al., 2017;Silver et al., 2017;Wang et al., 2018;Kolm and Ritter, 2019). Despite the rapid growth of RL in recent years, a common feature in most existing RL techniques, or broadly speaking most AI technologies, is that the learning task is performed in discrete learning environments (Russell and Norvig, 2010, Section 2.3.2), that is, either the time space or the state-action space or both are associated with the discrete topology. Although most practical physical environments are continuous in nature, in the sense that both the time space and the state-action spaces are continuous, RL and AI problems in continuous learning environments have rarely been touched upun in the past literature. In fact, the ability of learning in continuous learning environments has long been visioned as an important component of a futuristic AI system (Russell, 1997, Section 7). Admittedly, building AI and RL agents in continuous environments is a challenging task. One challenge is that most processes evolving in continuous environments are described by ordinary differential equations (ODEs), stochastic differential equations (SDEs), or other complex mathematical models. These models are still not fully-understood from the AI and RL perspective. This raises the question on how to develop a generic learning framework applicable to all these models. In addition, the state-action spaces encountered in practical RL problems often exhibit complex shapes and structures. For example, the presence of singularity points in robot workspace (Merlet, 2006) is a notorious issue in robot control and learning. These irregularities create a serious challenge for the convergence and robustness analysis of learning algorithms. Besides its potential value in practical applications, understanding the continuous limits of RL and AI algorithms also provides us a unique insight into the dynamics and theoretical nature of these algorithms. In fact, dynamical systems and differential equation theories have gained great popularity recently in understanding ML, especially neural network (NN) models (E, 2017;Chaudhari et al., 2018;Chen et al., 2018;Ruthotto and Haber, 2020). In addition, the ODE method has already been used in stochastic approximation (Ljung, 1977) to study the convergence of RL algorithms (Kushner and Yin, 2003). Hence, in order to expand the application of RL to a broader range of problems, and to enrich our understanding of RL theory, it is necessary to investigate RL methods in general continuous learning environments.
The purpose of this paper is to extend the classical temporal difference (TD) error (Sutton, 1988) to the setting of continuous environments. To overcome the shortcomings of traditional RL in dealing with problems arising from continuous environments, we will introduce a new concept under the name of temporal-differential error. Similar to the TD in discrete environments, the temporal-differential error plays a key role in developing new RL methods in continuous environments. In particular, we introduce two continuous RL algorithms, namely, the continuous-time least squares policy evaluation (CT-LSPE) (Algorithm 1) and the continuous-time temporal-differential (CT-TD) learning (Algorithm 2). These algorithms are extensions of the LSPE (Nedić and Bertsekas, 2003) and the TD learning (Sutton and Barto, 2018, Chapter 6), respectively, to continuous environments. As a continuous RL methodology, our algorithms also possess several unique features that distinguish them from most traditional RL methods. First, instead of using NN approximation, kernel approximation is used in our design. The linkage between the kernel function and the Hilbert space theory allows our algorithms to achieve convergence over the entire state space. Second, systems and control theory will be applied to show the stability and robustness of the proposed learning algorithms. Finite-time error bounds are also obtained, without assuming the stationarity of the training data as in the past literature (Antos et al., 2008;.

Related Works
The study on RL in continuous environments can be traced back to the concept of advantage updating (Baird, III, 1993;Harmon et al., 1995). Shortly after Baird's advantage updating, the continuous-time TD learning and actor-critic algorithms were developed by Doya (1996) and Doya (2000). In addition, Munos (2000) developed a continuous RL technique to solve viscosity solutions in order to tackle the case where the classical optimal solution does not exist. Continuous RL has also been applied in different fields including game theory (Börgers and Sarin, 1997) and neurobiology (Frémaux et al., 2013). A common feature of these early results is that the discrete-time sampling is performed to convert the continuous-time problem into the discrete-time setting. More recently, the research on deep RL opens possibilities for developing practical learning algorithms in continuous state-action spaces (van Hasselt and Wiering, 2007;van Hasselt, 2012;Lillicrap et al., 2016;Duan et al., 2016;Recht, 2019). Indeed, computer experiments from balancing inverted pendulums to learning the locomotion have shown very promising results.
Besides the above collective efforts, a learning-based control design methodology, usually under the name of adaptive (or approximate) dynamic programming (ADP) (Jiang and Jiang, 2012;Vrabie et al., 2013;Bian and Jiang, 2016;Kiumarsi et al., 2018), has been developed to solve optimal control problems for continuous-time dynamical systems. Instead of discretizing a continuous-time system into a discrete-time system, ADP algorithms solve directly the optimal control problem by utilizing the continuous-time data flow generated by dynamical systems. Exploiting explicitly systems and control techniques, the stability and optimality of the controlled system can be guaranteed. Finally, to perform the learning task in more complex environments, robust and decentralized extensions of ADP have also been proposed by Jiang and Jiang (2017); Bian and Jiang (2019).

Problem Formulation
Throughout this paper, R denotes the set of reals. I denotes the identity matrix with appropriate dimension. | · | denotes the Euclidean norm for vectors, or the induced matrix norm for matrices. id denotes the identity mapping. Let (Ω, F, P) be a probability space and {x t } t≥0 be a càdlàg and time-homogeneous strong Markov process with state space X. Here X is a locally compact Hausdorff space with a countable base, equipped with the Borel σ-algebra B. Denote by E x the expectation conditionally on the initial state x and by E µ the expectation with respect to µ. The Hilbert space L 2 µ and the norm · µ are induced from the inner product ·, · µ defined via f, µ . For more details on the notations used in this paper, see Appendix A and references therein.
Our goal in this paper is to estimate the following value function V * that depends on the future path of {x t } t≥0 : where r ∈ L 2 µ is the running-reward/cost function, and γ > 0 is a scalar representing the discounting effect. To be mathematically concrete, we assume E x [r(x t )] is measurable and locally essentially bounded.
To see why the classical discrete-time TD learning is no longer a good choice to solve V * here, we write the discrete-time TD error for the discretized process (with sampling period ∆t) as Obviously, when ∆t is small, the discounting term e −γ∆t is close to 1, resulting in a poor convergence performance (see Tsitsiklis and Van Roy (1997) for convergence analysis). If ∆t is too large, the updating frequency is reduced, which could also compromise the on-line learning performance. In light of the above difficuities, it is desirable to solve V * without discretizing the performance index and the underlying Markov process in the first place. To proceed, note from Ito (1960, Chapter 1) that V * satisfies the following linear functional equation: where T = A − γ. The converse statement is also true by Dynkin's formula (Kushner, 1967a, Chapter 1). In this paper, we will introduce two on-line learning methods to solve (2).

Value Approximation in Continuous Time
In this section, we consider the following approximated equation: where Π is defined by the following integral equation: Here K is a square-integrable continuous symmetric positive-definite kernel. Corresponding to each K, we can uniquely define a reproducing kernel Hilbert space (RKHS) H ⊆ L 2 µ (Aronszajn, 1950, Section 2). In addition, Π is an orthogonal projection from L 2 µ to H (Aronszajn, 1950, Section 2). To fit into our learning framework, the reproducing kernel K is designed to ensure H is a subspace of D(A). This requirement can be satisfied for many practical processes by selecting a sufficiently smooth K; see Section 5 for details. Finally, we assume the initial value function V 0 ∈ H.
From the perspective of systems and control theory, (3) can be viewed as a reduced-order model (Scarciotti and Astolfi, 2017) As a result, we can rewrite (3) in the form of system (4) with an approximation error (Π − id)(T V t + r) added to its right-hand side.
The following theorem shows that V t converges into a neighbourhood of V * exponentially.
Theorem 1 V t is well-defined and is exponentially stable (at V ∞ ) in H. In addition, See Appendix B.1 for the proof of Theorem 1. The finite-time error bounds in Theorem 1 are composed with two parts. The first part represents the error induced from the initial guess V 0 . The second part is due to the approximation error introduced by the projection Π. If we can choose H properly, such that for any f ∈ H, Af ∈ H, then one can show that The above inequality implies that in this case the error bound is purely controlled by the initial guess V 0 and the "observation error" in the running-reward/cost r. Note that the first error bound in Theorem 1 is tighter for large t, and the second error bound is tighter for small t. From the perspective of robust control theory, the second error bound quantifies the impact of (Π − id)(T V t + r) on the approximation error between V t and V * . In addition, the discounting factor γ is linked to the robustness of (3), in the sense that a larger γ results in a smaller error bound and a faster convergence rate.
Finally, it is worth pointing out that the results in Theorem 1 can also be extended to evaluate the ergodic cost (Arapostathis et al., 2012). In this case, there is no discounting term in (1), and we have T = A in (2) and (3). Note from Bhattacharya (1982, Proposition 2.2) that as long as {x t } t≥0 is ergodic, A has a simple zero eigenvalue corresponding to the constant eigenfunction. Hence, by excluding the constant function from H, we can ensure there still exists a γ > 0, such that V, ΠT V µ ≤ −γ V µ . As a result, the two inequalities in Theorem 1 still hold.

Model-free On-line Learning
The results presented in the last section provide an efficient way to approximate V * through the linear functional equation (3). Unfortunately, in order to solve (3) directly, we must have access to the knowledge of µ, A, and r. This is not an easy task in practice, since µ, A, and r are associated with the mathematical model of the learning environment, which is usually not known.
In this section, we propose two on-line learning algorithms to estimate V * using the data observed directly from the environment. The knowledge of µ, A, and r is not required in our learning algorithm design. Throughout this section, we consider the following realization of K: where φ = [φ 1 φ 2 · · · φ N ] T , and {φ j } N j=1 are N independent functions in D(A). The independence of {φ j } N j=1 guarantees the matrix inverse in K is well-defined. Obviously, K satisfies our definition of reproducing kernel in Section 4. In fact, the RKHS H associated with K is the N -dimensional space spanned by {φ j } N j=1 , and V t in (3) can be parameterized as V t (x) = φ T (x)c t for some c t ∈ R N . Directly plugging the definitions of K and V t in (3), we have the following ODE in c t : In the above equation, T is applied element-wise to φ. By Theorem 1, (5) admits an equilibrium point c * satisfying E µ φT φ T c * + E µ [φr] = 0. In addition, by Dynkin's formula, we can replace Aφ in (5) with where M t is a vector-valued martingale, and dφ t φ t+dt − φ t represents the infinitesimal change in the value of φ t . By abuse of notation we denote φ t and r t for φ(x t ) and r(x t ), respectively, when there is no confusion.
Algorithm 1: Continuous-time least squares policy evaluation (CT-LSPE) Update c t through the following ODEs:

Continuous-time least squares policy evaluation
We first replace the integrations with respect to µ in (5) by their empirical estimations. This, together with (6), leads to the CT-LSPE algorithm in Algorithm 1. The convergence of Algorithm 1 relies on the existence of three functions ψ i , i = 1, 2, 3, defined as These functions quantify the differences between the expectations in (5) and their corresponding empirical estimations. Now, we impose the following assumptions on x t and φ.
Assumption 1 ensures that x t is an ergodic process. Assumption 2 requires that ψ i , i = 1, 2, 3, exist and are solutions to certain Poisson equations; see Glynn and Meyn (1996) for conditions under which ψ i exists. Note that the discrete-time versions of Assumptions 1 and 2 for MDPs have been widely used in ADP and RL literature (Tsitsiklis, 1994;Tsitsiklis and Van Roy, 1997). Assumption 3 essentially requires the second moment of the stochastic integral t 0 φ j (x s )dM s , j = 1, 2, · · · , N , does not grow too fast. Since we can always choose a sufficiently smooth φ (such as the softmax function, the Gaussian function, or the wavelets) to saturate the noise in the environment, it is possible to satisfy this assumption even for a very noisy process (such as a fat-tailed process).
The convergence of Algorithm 1 is guaranteed in the following theorem.
Theorem 2 Consider Algorithm 1. Under Assumptions 1, 2, and 3, c t converges to c * with probability one. In addition, there exist constants C i > 0, i = 0, 1, 2, such that for a sufficiently large , Algorithm 2: Continuous-time temporal-differential (CT-TD) learning Update c t through the following SDE: with probability at least 1 − C 0 / 2 , where λ M and λ m represent the largest and smallest eigenvalues of E µ [φφ T ], respectively, and ε < γ/λ M is an arbitrary positive constant. In particular, (7) implies that |c t − c * | = O t −0.5 with probability one.
See Appendix B.2 for the proof of Theorem 2 and the detailed formulation of the error bound. The first term on the right-hand-side of (7) represents the error induced from the prior guess c 0 . The second term is driven by the estimation error induced from the empirical estimation of the expectations in (5). By the law of large numbers, this estimation error converges to 0 with rate O t −0.5 . Since t −0.5 decreases slower than the exponential rate, it dominates the error bound of |c t − c * |. Note that the obtained convergence rate is consistent with the result in Yu and Bertsekas (2009) for discrete-time LSPE.

Continuous-time temporal-differential learning
An alternative way to approximate (5) via on-line data is to use stochastic approximation (Benveniste et al., 1990). In this case, E µ φT φ T c t + E µ [φr] is directly replaced by the on-line data φ t c T t dφ t + φ t r t at time t. This, together with (6), leads to the CT-TD learning algorithm in Algorithm 2. To accommodate the noise induced from stochastic approximation, the matrix inverse in (5) is replaced by a slowly decreasing step size α t . In particular, we name as the temporal-differential error.
Assumption 4 (Finite moments) There exist C > 0 and real-valued function ϕ, such that for the j-th element in φ and ψ i , and i = 1, 2, 3.
A discrete-time version of the boundedness condition on the high-order moments of ψ i and Aψ i is also required in Tsitsiklis and Van Roy (1997) and Benveniste et al. (1990). The conditions on (dφ j (x t )) 2 and (dψ i,j (x t )) 2 are new in this paper and essentially require ψ i (x t ) and φ t do not vary too fast.
Theorem 3 Consider Algorithm 2. Under Assumptions 1, 2, and 4, c t converges to c * with probability one. In addition, there exist constants C 0 and C 1 , such that for sufficiently large and integer k, In particular, (8) implies that |c t − c * | = O e −γ t 0 αsds with probability one.
See Appendix B.3 for the proof of Theorem 3 and the detailed formulation of the error bound. Note that as k goes to the infinity, A k converges to 0 by the conditions on α t in Algorithm 2. In addition, the convergence rate of A k is slower than the exponential rate, due to the slow decreasing rate of α t . As a result, lim k→∞ √ A k e 2C 1 k = ∞, and the error bound in (8) grows with k.
A key difference between Algorithms 1 and 2 is that the time integration of the whole time series is incorporated in Algorithm 1, while only the on-line data at current time is used in Algorithm 2. As a result, Algorithm 2 may produce a much noisier estimation on the value function, especially at the beginning of the learning process. Indeed, the experimental results in Section 6.1 have shown that the point estimations produced by Algorithm 1 have smaller standard errors. In fact, to ensure the convergence of Algorithm 2, a more restrictive condition (Assumption 4) is required. In contrast, the presence of decreasing step size α t in Algorithm 2 provides more freedom to adjust the convergence rate of the CT-TD learning.
We can see from the proofs of Algorithms 1 and 2 that the constants in the error bounds in Theorem 2 and Theorem 3 are polynomials of N , that is, the number of bases. As a result, if N increases, the error bounds in Theorems 2 and 3, which represent the estimation errors, increase as well. On the other hand, a small N also leads to a large divergence from H to L 2 µ , which is quantified by (Π − id)(T V s + r) µ in Theorem 1. As a result, a small N causes a large approximation error. This trade-off between approximation error and estimation error, or better known as the bias-variance trade-off, is a well-known phenomenon in various machine learning methods.
Finally, note from Algorithms 1 and 2 that the state x t of the environment does not appear explicitly in the updating equations. Instead, the updating equations only depend on r(x t ) and φ(x t ). From the perspective of control theory, the pair (φ(x), r(x)) can be viewed as the output of the environment at state x. Therefore, Algorithms 1 and 2 can also be applied to output-feedback control problems and partially observable Markov processes.

Computer-based Experiments
In this section, we present two computer-based experiments to illustrate the effectiveness of the two learning algorithms presented in the previous section.
6.1 Continuous-time ARMA(2,1) process Consider {y t } t≥0 as a real-valued continuous-time ARMA(2,1) process governed bÿ where θ 1 , θ 2 , σ 1 , σ 2 , and σ 3 are positive reals, and ξ i,t , i = 1, 2, 3, are three independent white noises (Arnold, 1974, Section 3.2). Denote x t = (y t ,ẏ t − ξ 1,t ). Then, we can derive the following two-dimensional Ornstein-Uhlenbeck (OU) process: We know from Kushner (1967b) that the discounted value function V * defined by (1) solves the following partial differential equation: In this example, we aim at using Algorithms 1 and 2 to estimate the discounted value functions with respect to three types of running-costs. First, a quadratic function r(x) = 10 2 i=1 |x i | 2 is used, under which (9) reduces to a Lyapunov equation. Then, we select the L 1 vector norm , which is essentially the 0-1 loss. See Figure 1 for an illustration of different choices of r. In all three cases, the discounting factor is fixed at γ = 1, and the step size α t in Algorithm 2 is chosen as α t = 0.7/(1 + t). Model parameters of the OU process are chosen as θ 1 = θ 2 = σ 1 = σ 2 = σ 3 = 1, and the initial state is chosen as x 0 = (1, 0). Six basis functions are selected: Figures 2 and 3 show the experimental results of our algorithms. To evaluate the performance of the algorithms, we conduct learning processes over 35 sample paths. Figure 2 illustrates the 35 paths of the weights {c i } 6 i=1 corresponding to the six basis functions. Figures 3a, 3b, and 3c show the box plots of these weights at t = 1000. Overall, the weights obtained from both the CT-LSPE algorithm (Algorithm 1) and the CT-TD learning algorithm (Algorithm 2) converge to same values in all three cases. However, compared with the CT-LSPE, the CT-TD learning generates more dispersed samples, and the outliers are scattered further from the center of the clusters of samples. This is not surprising, as the time integration in the CT-LSPE tends to smooth out the noise in the learning process, especially at the beginning of the learning phase. Indeed, we can see from Figure 2 that CT-TD learning generates noisier paths of basis-function weights. To further illustrate the effectiveness of the proposed algorithms, we compare the value functions obtained from the two learning algorithms with that from Monte-Carlo (MC) simulations in Figures 3d, 3e, and 3f. The point estimations of value functions at t = 1000 with 95% confidence intervals are given in Figures 3g, 3h, and 3i. Note that in all three cases, these point estimations are approximately at the same level. In addition, compared with the CT-TD learning, the CT-LSPE produces

Control of the benchmark double inverted pendulum on a cart
In this example, we apply the CT-LSPE together with the actor-critic algorithm (Konda and Tsitsiklis, 2000) to improve the control performance of a double inverted pendulum on a cart (DIPC) (Figure 4). It is well-known that the DIPC is a nonlinear under-actuated plant, and designing a controller for DIPC is a challenging task in the field of nonlinear control (Khalil, 2015). Figure 4: Double inverted-pendulum on a cart.
The mathematical model of the DIPC is given below (Moysis, 2016), which is derived from the Lagrange equations of the DIPC: where θ = [θ 0 θ 1 θ 2 ] T , θ 0 denotes the horizontal displacement of the cart, θ 1 and θ 2 denote the angles of lower and upper pendulum links with respect to the vertical position, respectively, u is an input force applied to the cart, and The model parameters are listed in Table 1. Denote the state of the DIPC as x = [θ 0 θ 1 θ 2θ0θ1θ2 ]. Now, we can rewrite the above model in the standard state-space form: Note that D, C, and G are also nonlinear functions of x. Hence, the above system is an affine nonlinear system. In the learning task, at each state-action pair (x, u), the following r(x, u) = x T Qx + u 2 , Q = diag{100, 1000, 1000, 50, 10, 10}.
In another word, we can always write u t = Kx t for some real control gain matrix K.
The entire learning process is composed with 50 learning trials, indexed by k = 1, 2, · · · , 50. Each learning trial is performed over a fixed time interval [0, t f ]. The initial controller parameter K 1 is adopted from Moysis (2016): The CT-LSPE algorithm is employed in each learning trial to estimate the value function. To facilitate the learning process, we add exploration noises in the system input in each learning trial. As a result, the control action applied to the DIPC in the k-th learning trial is u t = K k x t + noise t , where noise t is a stationary Gaussian process with fixed distribution N (0, σ) for all t ∈ [0, t f ]. Then, u t has a distribution N (K k x t , σ) π k (a t ) in the k-th trial. Extending the actor-critic algorithm with eligibility trace in Konda (2002) to the continuous-time setting, we update the control gain matrix after finishing the k-th learning trial via where h > 0 is the step size, V t f is the value function learned in the k-th trial, and is the eligibility trace. Note that we only update the control gain matrix K k at the end of each learning trial in order to increase computational efficiency. The experimental results are given in Figures 5 and 6. We fix h = 0.0005 and t f = 800 during the learning process. After 50 learning trials, the control gain matrix becomes    Trial 1 Trial 50 Figure 6: The weights of basis functions ({c i,t } t≥0 for i = 1, 2, · · · , 10) for the DIPC, in the first and the 50-th learning trials.
The cost values V t f (x 0 ) estimated in each learning trial are plotted in Figure 5a. Clearly the cost values show a decreasing trend as the trial number increases. We also plot the system trajectories under controllers u 1 = K 1 x t and u 50 = K 50 x t in Figure 5b, respectively. Obviously, u 50 leads to a much better transient performance, in the sense that the DIPC shows less oscillations. Finally, the weights of basis functions ({c i,t } t≥0 for i = 1, 2, · · · , 10) before and after conducting the learning process are plotted in Figure 6.

Conclusions
This paper is motivated to provide a new RL framework for learning problems in continuous environments. A novel concept, known as the temporal-differential error, and a new class of temporal-differential learning methods, are proposed. In particular, two RL algorithms with detailed convergence analysis are designed for RL in continuous environments. Next, we point out several future research directions that deserve further investigations. First, the proposed algorithms are purely based on the temporal-differential error. It is interesting to see how to incorporate eligibility traces in RL algorithms driven by the temporal-differential error. Second, it is worth checking how to extend other existing discrete-time TD-learning-based RL algorithms to the continuous setting, by means of the techniques developed in this paper. Finally, besides the prediction problem which is the main focus of this paper, it is also important to investigate the off-policy learning problem in continuous environments.

Acknowledgments
The work of Z. P. Jiang has been supported partially by the National Science Foundation under Grants ECCS-1501044 and EPCN-1903781. The opinions expressed in this paper are those of the authors and do not necessarily reflect the views and policies of Bank of America.

Appendix A. Review of Continuous-time Markov Processes
Let (Ω, F, P) be a probability space. Consider {x t } t≥0 as a homogeneous càdlàg (right continuous with left limits) strong Markov process on state space X, which is a locally compact Hausdorff space with a countable base and is equipped with the Borel σ-algebra B. The definition of X covers most spaces that will be encountered in RL tasks, including Euclidean spaces, countable discrete spaces, and topological manifolds. Starting with x 0 = x, x t admits a stationary probability measure µ on B. Denote by E x the expectation conditionally on the initial state x and by E µ the expectation with respect to µ.
[M ] t denotes the quadratic variation of a real-valued stochastic process M t defined on (Ω, F, P).
Given the stationary distribution µ, we can define a Hilbert space L 2 µ equipped with the inner product ·, · µ and the norm · µ : In addition, we can define an operator P t : L 2 µ → L 2 µ as One can easily check that P t is a contraction semigroup on L 2 µ (see Lemma 4 in Section C). In addition, the infinitesimal generator A associated with P t is defined as provided that the above limit exists. Denote the domain and the range of A as D(A) ⊆ L 2 µ and R(A) ⊆ L 2 µ , respectively. By Hille-Yosida theorem (Pazy, 1983), we know A is closed, and D(A) is a dense subset of L 2 µ . Since x t is right continuous, its transition probability is stochastically continuous, and is uniquely determined by its infinitesimal generator. For additional properties on P t and A, see Nelson (1958); Pazy (1983); Hansen and Scheinkman (1995).

Appendix B. Proofs of Theorems
In this section, we present proofs of our three main theorems.

B.1 Proof of Theorem 1
We first show that V t exists on H. For any V ∈ H, we have from the definition of K that Since A generates a contraction semigroup P t in L 2 µ (Lemma 4), we know by Lumer-Phillips theorem (Lumer and Phillips, 1961) that A is closed and dissipative. Hence, V, AV µ ≤ 0. Thus, ΠA is dissipative on H, and −ΠA is maximal monotone (Borwein, 2010, Proposition 1). Then, we have by Minty surjectivity theorem (Minty, 1962) that ΠA − id is surjective. In addition, since A is closed, we know by the definition of Π and the dominated convergence theorem that ΠA is closed on H. Hence, again by Lumer-Phillips theorem, ΠA generates a contraction semigroup on H. Then we have from Phillips (1954) that the abstract Cauchy problem (3) admits a unique solution V t . In addition, by Hille-Yosida theorem, (ΠA − γ) −1 is bounded, and hence 0 = Π(TV * + r) = (ΠA − γ)V * + Πr admits a unique solutionV * on H.
DenoteṼ t = V t −V * . We next show that V t converges toV * exponentially. Note from (3) and the above analysis that This implies that V t is globally exponentially stable atV * . Now, to derive the first finite-time error bound, we have Denote 0 = T (V * ) + r + ε, where ε is the approximation error due to Π. Using Dynkin's formula and the tower property, we havê Since µ is the stationary distribution, we have by Jenson's inequality and Fubini's theorem that To derive the second finite-time error bound, we rewrite (3) as where ε t is the approximation error at time t. Using Dynkin's formula, we have Then, where (10) and (11) above are due to Jenson's inequality. This completes the proof.

B.2 Proof of Theorem 2
We first show solves the ODE of R t in Algorithm 1. Since the right-hand-side of the ODE of R t is locally Lipschitz, and t 0 φ s φ T s ds + ρI −1 is bounded, we can select a sufficiently large subset S of the space of all symmetric positive definite matrices, such that is a unique solution of the ODE of R t on S, and it also remains in S.
We next denote Then, we can rewrite the updating equation of c t in Algorithm 1 aṡ where By Assumptions 2, we have where i = 2, 3, 4, ψ 4 = ψ 1 − γψ 3 , and Then, we have by Lemma 6 (element-wise) that Now consider the boundedness of η t . The quadratic variation of the j-th element in M t is given below: By Assumption 3 and Ito's isometry, we have sup t≥0 (t + 1)E x |η t | 2 < ∞.
Applying Lemma 5 with α t = 1/ √ t + 1, we know with probability one that A t , b t , and Σ t converge to A, b, and Σ, respectively. In addition, we know with a probability at least of As a result, by the boundedness of A t , b t , Σ t , and η t , we have with a probability at least Applying Lemma 7, we know c t converges to c * with probability one. In addition, with a probability at least 1 − C 0 / 2 , where λ M and λ m represent the largest and the smallest eigenvalues of Σ, respectively, ε < 2γ/λ M is an arbitrary positive constant, and by Lemma 9, In particular, we know by Lemma 10 that |c t − c * | 2 = O t −1 . This completes the proof.

B.3 Proof of Theorem 3
First, we introduce the following auxillary system: where A and b follow the definitions in Appendix B.2. We first show thatĉ t converges to c * . Indeed, a direct calculation shows that By the comparison lemma (Sontag, 1998, Lemma C.3.1), we have Since ∞ 0 α t dt = ∞,ĉ t converges to c * as t goes to the infinity. Denotec = c −ĉ. We can rewrite the updating equation in Algorithm 2 as where ∆ 3,t and ∆ 4,t follow the definitions in Appendix B.2.
Note that c T dM = j c j dM j , where dM j denotes the j-th element in dM . Now the proof is completed by applying Lemma 8.

Appendix C. Technical Lemmas
In this section, we present several supporting lemmas that have been used in the proofs of our main results.
Lemma 4 P t is a contraction semigroup in L 2 µ .
Proof We have by Jenson's inequality that This completes the proof.
Lemma 5 Suppose S t is a time series satisfying for some C < ∞, where α t > 0 is deterministic and decreases monotonically to 0. Then, for any > √ C, In particular, α 2 t S t = O(α t ) with probability one.
Proof Denote events By Markov's inequality, P (E T ) ≤ C/ 2 . Since E T is increasing, by monotone convergence theorem, Letting go to the infinity, we have This completes the proof.
Then m is a martingale, and we have by the martingale property that where 0 = t 1 < t 2 < · · · = t is a set of partition points on the interval [0, t], ∆t i = t i+1 − t i , and Aψ(x s )ds.
We next quantify E x [m 2 t ] by inspecting E x [∆m 2 i,t ]. By Hölder's inequality, where the last equality holds since the integrand is bounded. In addition, Putting the above items back into the summation in E x [m 2 t ] and letting ∆t i go to 0, we have Since both ψ and f are square integrable, lim sup t→∞ t −1 E x [m 2 t ] < ∞. This completes the proof.
Lemma 7 Consider the following system: where Σ is a symmetric positive definite matrix, A is a negative definite real matrix, and ∆A t and ∆ t are random variables bounded with probability one. Assume for all ω ∈ E ∈ F. Then c is input-to-state stable (ISS) (Sontag, 2008) where γ > 0 satisfies A T + A ≤ −γI, and ε < γ|Σ| −1 is a positive constant. Reformulating the above inequality, we have Hence, where λ M and λ m represent the largest and the smallest eigenvalues of Σ, respectively. One can see from the above inequality and our assumption on ∆A that c t is ISS, from ∆ t to c t . This completes the proof.
Lemma 8 Consider the following system: where A is a real matrix satisfying A + A T < −2γI for some γ > 0, α t > 0 is continuously differentiable withα t ≤ 0, lim t→∞ α t = 0 and ∞ 0 α t = ∞, g j is a vector-valued function, M j is a real-valued martingale, and ∆ 1,t and ∆ 2,t are random variables that are bounded with probability one and satisfy for some martingales m i , i = 1, 2. In addition, M j and m i are adapted to a filtration {F t } t≥0 , and there exist a function ϕ and a constant C > 0, such that d[M j ] t ≤ ϕ t dt and d[m i,k ] t ≤ ϕ t dt for the k-th elements in m i , and for all c, where E 0 denotes the expectation conditional on F 0 . Given the above conditions, we have C 0 is a positive constant, and n and are sufficiently large integers and reals, respectively.
This completes the proof.
Lemma 9 For any 0 ≤ s < t, a > 0 and b > 0, we have If a > (2 √ b) −1 , thenẏ t > 0 for all t. As a result, y t > y s for all t > s. If a ≤ (2 √ b) −1 , then y t − y s ≥ min t y t − y 0 , where min t y t = −ab − (4a) −1 . This completes the proof.
Lemma 10 Suppose a is a positive real. Then for any t 0 > 0,