Impact of Computation in Integral Reinforcement Learning for Continuous-Time Control
Authors: Wenhan Cao
Wei Pan
Venue: ICLR 2024 spotlight (2024)
Keywords: Integral Reinforcement Learning
Bayesian Quadrature
Newton's Method
Primary area: reinforcement learning

Abstract: Integral reinforcement learning (IntRL) demands the precise computation of the utility function's integral at its policy evaluation (PEV) stage. This is achieved through quadrature rules, which are weighted sums of utility functions evaluated from state samples obtained in discrete time. Our research reveals a critical yet underexplored phenomenon: the choice of the computational method -- in this case, the quadrature rule -- can significantly impact control performance. This impact is traced back to the fact that computational errors introduced in the PEV stage can affect the policy iteration's convergence behavior, which in turn affects the learned controller. To elucidate how computation impacts control, we draw a parallel between IntRL's policy iteration and Newton's method applied to the Hamilton-Jacobi-Bellman equation. In this light, computational error in PEV manifests as an extra error term in each iteration of Newton's method, with its upper bound proportional to the computational error. Further, we demonstrate that when the utility function resides in a reproducing kernel Hilbert space (RKHS), the optimal quadrature is achievable by employing Bayesian quadrature with the RKHS-inducing kernel function. We prove that the local convergence rates for IntRL using the trapezoidal rule and Bayesian quadrature with a Matérn kernel to be $O(N^{-2})$ and $O(N^{-b})$, where $N$ is the number of evenly-spaced samples and $b$ is the Matérn kernel's smoothness parameter. These theoretical findings are finally validated by two canonical control tasks.

Published as a conference paper at ICLR 2024
IMPACT OF COMPUTATION IN INTEGRAL REINFORCE- MENT LEARNING FOR CONTINUOUS- TIME CONTROL
Wenhan Cao
Wei Pan"
I DDpaarmmett of Computer Science, University of Manchester Schooo of Vehicle and Mobility, Tsinghua University
owh1 191mmmaass.....m.mm ts tsinghua. edu. en wei.ppanemanchesterraa.. uk
ABSTRACT
Integral reinforcement learning (IntRL) demands the precise computation of the utility function's 's integral at its policy evaluation (PEV) stage. This is achieved through quadrature rules, which are weighted sums of utility functions evaluated from state samples obtained in discrete time. Our research reveals a critical yet underexplored phenomenon: the choice of the computational method in this case, the quadrature rule can significantly impact control performance. This im- pact traced back to the fact that computational errors introduced in the PEV stage can affect the policy iteration's convergence behavior, which in turn af- fects the learned controller. To elucidate how computation impacts control, we draw a parallel between IntRL's policy iteration and Newton's method applied to the Hamilton- Jacoon--acoi---a Belbi-Bellman equation. In this light, computational error in PEV manifests an error term each iteration of Newton's method, with its up- per bound proportional to the computational error. Further, we demonstrate that when the utility function resides in reproducing kernel Hilbert space (RKHS), the optimal quadrature achievable by employing Bayesian quadrature with the RKHSSinducing kernel function. We prove that the local convergence rates for In- tRL using the trapezoidal rule and Bayesian quadrature with Matérn kernel to be O(N-2) and O(N~ where N is the number of evenly-spaced samples and b is the Matérn kernel's smoothness parameter. These theoretical findings finally validated by two canonical control tasks.
INTRODUCTION
Recent advancements in reinforcement learning (RL) have prominently focused on discrete-time (DT) systems. Notable applications range from Atari games [Schrittwisser et al. (2020) the game
modees of Go Silver et al. (2016) 2017) as well as large language models Bubeck et al. (2023). However,
most physical and biological systems are inherently continuous in time and driven by differential equation dynamics. This inherent discrepancy underscores the need for the evolution of continuous- time RL (CTRL) algorithms/ Baird (1994); [ewis et al. (1998); Abu-Khalaf & Lewis (2005); Vra-
bie & Lewis 2009); | Vrabie et al. 2009); Lewis & Vrabie (2009); Vamvoudakis & Lewis (010)):
Modares al. ( L et al. (2014); 4 Modares & Lewis (2014); Vamvoudakis al. (2014 Yildiz
ct al. (2021) Holt et al. (2023): y Wallace & (2023).
Unfortunately, adopting CTRL presents both conceptual and algorithmic challenges. First, Q- functions are known to vanish in CT systems Baird (1994); A
Lewis (2005); Lewis & Vrabie (2009), which makes even simple RI. algorithms such as Q-learning infeasible for CT systems. Second, the one-step transition model in DT system needs to be replaced with time deriva- tives, which leads to the CT Bellman equation (also called nonlinear Lyapunov equation or general- ized Hamilton-Jcobi Bellman equation) governed by a complex partial differential equation (PDE) rather than simple algebraic equation (AE) as in DT systems Lewis & Vrabic (2009):
Vrabie &
Lewis (2009): Vamvoudakis & Lewis 2010):
CT Bellman Equation (PDE): (=((()) = (xx(t), u(x(t)), DT Bellman Equation (AE): (V"((()) = -[(()((),,(((())))
(la)
(Ib)
In Гap and (Ib), the symbols x, u, 1, V represent the state, control policy, utility function and value function, respectively. While traditional DTRL. focuses maximizing "reward function".

---

Published as a conference paper at ICLR 2024
work, in line with most CTRL literature, employs "utility function" to represent costs or penalties, aiming to minimize associated value. Solving the Bellman equation refers to the policy evaluation (PEV), which is a vital step in policy iteration (PI) of RL [Sutton et al. (1998); Lewis & Vrabie
2009 However, CT Bellman equation cannot be solved directly because the explicit form of
(F(() u(x(())) hinges on the explicit form of the state trajectory a æ((), which is generally unknown.
The CT Bellman equation can be formulated as the interval reinforcement form Vrabie & Lewis
2009); : Lewis & Vrabie (2009); Modares et al.] (2014); MMdarrs & Lewis (2014); Vanvoudakis ct al. 2014). which computes the value function through the integration of the utility function:
()((()) = E(1) + (=(( + AT)), E(1) =
[x((s), u(8))ds.
(2)
Here, E(1) represents the integral of the utility function. When the system's 's internal dynamics is unknown, E(U) can only be computed by quadrature rule solely utilizing state samples in discrete time, i.e..
N {(1) ~ E(1) Σ (.:(a(tt),urx(tt))
(3)
Here, the time instants for collecting state samples are 11 = < l2 < <... < iN = t AT, with N being the sample size. The quadrature rule characterized by a set of weights ] which can be chosen from classical methods like the trapezoidal rule or advanced probabilistic methods such as Bayesian quadrature (BQ) O'Hagan (1991): Karvonen & Särkkä 20 Cockayne et al. (2019); Briol et al. (2019): Hennig et al. (2022). For simplicity, we denote the computational error for the quadrature rule as Err(E(t)) == [())- - É()). In real-world applications, sensors of autonomous systems are the primary source of state samples. For example, state samples are gathered from var- ious sensors when the CTRL algorithm trains a drone navigation controller. These sensors include the gyroscope for orientation, the accelerometer for motion detection, the optical flow sensor for positional awareness, the barometer for altitude measurement, and GPS for global positioning, If these sensors operate at a sampling frequency of 10 Hz and the time interval AT in (2) is set second, we would obtain N = 11 state samples within that duration.
The impact of computational methods on the solution of the CT Bellman equation is emphasized by Yildiz jct al. (202I). Echoing this finding, our study further reveals that the choice of computational method can also impact the performance of the learned controller. As illustrated in Figure 1 we investigate the perfor- mance of the Integral Reinforcement Learning (InRL) algorithm Vrabie & Lewis (200)). We utilize both the trapezoidal rule and the BQ with a Matérn kernel to compute the PEV step, applying different sample sizes N. Our results for canonical control task (Example 1 of Vrabie & Lewis 2009)) indicate that larger sample sizes diminish accumulated costs. In addition, we ob- serve notable differences in accumulated costs between these two quadrature methods. This trend highlights a crucial insight: the computational method itself can be a determining factor in the control performance. This phenomenon is not exclusive to the IntRL algo- rithm but applies to the CTRL algorithm when internal dynamics is known, as elaborated in Appendix
10
N=2 (Matèm, Trapz)
N=3 (Matëm, Trapz)
N=4 (Matêm, Trapz)
100²
10"
0
2
Iteration
Figure 1: Evaluation of accumu- lated cost J for controllers, computed through the trapezoidal rule and BQ with a Matérn kernel, compared to the optimal controller cost J" This sim- ulation is performed on different sam- ple sizes N and relies on evenly-spaced samples within a canonical control task in Vrabie & Lewis 2009).
The impact of computational methods - specifically, the quadrature rule as described in for unknown internal dynamics scenarios - on control performance is a compelling yet underexplored topic. This gap prompts crucial questions: How does the computation impact control performance? How can we quantitatively describe this impact? Despite the importance of these questions, existing literature has yet to address them. Our paper aims to bridge this gap by investigating the intricate relationship between computational methods and controller performance. We choose to focus our analysis on the IntRI. algorithm Vrabie & Lewis (2009). This choice is motivated by IntRL's role as a foundational CTRL algorithm in unknown internal dynamics scenarios Modares et al. (2014); Lee et al. (2014); Modares & Lewis (014)) Vamvoudakis et al. (2014); Wallace & Si (023)) As
2

---

Published as a conference paper at ICLR 2024
will illustrate in Section 2 understanding the impact of computation control in scenarios with unknown internal dynamics is more crucial and merits deeper examination.
The contributions of this paper can be summarized as follows: (1) We show that the PEV step in IntRL essentially computes the integral of the utiliiyy function. This computation relies solely on state samples obtained at discrete time intervals, introducing an inherent computational error. This error is bounded by the product of the integrand's norm in the reproducing kernel Hilbert space (RKHS), and the worst-case error. When employed as the quadrature rule, BQ minimizes the worst- case error, which coincides precisely with BQ's posterior covariance. Additionally, we analyze the computational error's s convergence rate concerning the sample size, focusing on both Wiener and Matérn kernels for evenly spaced samples. (2) We demonstrate that PI of IntRL can be interpreted as Newton's 's method applied to the Hamilton-Jacobr Jacobi-Belllaa (HJB) equation. Viewed from this perspective, the computational error in PEV acts as extra error term in each Newton's method iteration, with its upper bound proportional convergence rate of PI for IntRI., in terms of sample size, for both the trapezoidal rule and BQ using the Matérn kernel. Specifically, the convergence rate for PI is O(N when using the trapezoidal rule, and O DN- when using BQ with a Matérn kernel, where b is the smoothness parameter. Our paper the first to explore the impact of computational errors on the controller learned by IntRL.
the computational error. (3) We present the local
The remainder of this paper is organized as follows: Section II provides a detailed problem formu- lation. Our main results are presented in Section III, followed by simulations in Section IV. Section V draws conclusions and discussions on this new topic.
2 PROBLEM FORMULATION
We assume that the governing system control affine and the system's internal dynamics f is not necessarily known throughout the paper:
a = f(x) + g(z)u,
(4)
where x E Ω C R' is the system state, u E I is the control input, f R + R== and g IR + IR Xn. are locally Lispchitz mappings with f(0) O. The objective is to find a control policy u that minimizes the performance index:
J (To, u) = 1/' (((ss),(s)) ds, (0) = Fo.
(5)
where 1(x, = ((x) + u R(x)u is the utility function, with g(z) and R(x) being the positive
definite function and positive definite matrix, respectively. Furthermore, the system (4) is assumed
to be zero-state observable through q Beard et al.] (1997). The optimal control policy u" and its
associated optimal value function V are defined as:
(zo) == min JIH, u), u =x(x)) AC)(
arg min JIo, u).
(6)
(EA(
Here, u E A(2)) represents that u is an admissible policy on Ω. For more details of the definition of an admissible policy. please refer to Appendix B The optimal value function can be determined by the HJB equation, which is illustrated in the following theorem:
Theorem 1 (HJB Equation Properties Vrabie & Lewis (2009)). unique, positive definite, and continuous function V* serves as the solution to the HJB equation:
G(V) (V.V) - gR JV.V +
HJB Equation: G(V*)=
(7)
In this case, V is the optimal value function defined in (6). Consequently, the optimal control policy can be represented as (x) = 1 R- g(x)) V.
-
While solving the HJB equation poses a challenge as it lacks an analytical solution, IntRI. tackles this issue by utilizing PI to itcratively solve the infinite horizon optimal control problem. This approach is taken without necessitating any knowledge of the system's internal dynamics / Vrabie
& Lewis (2009). Consider u (x((t) A(2)) be an admissible policy, PI of IntRL performs PEV and policy improvement (PIM) iteratively:
PEV: V (x(1)=
(xx(s), ( (x(s)) ds + V (x(1+ + AT)), (0(x)= 00, (8a)
=
PIM: n (++)(()
= -- R /1x)
1
V. (
(8b)
3

---

Published as a conference paper at ICLR 2024
assesses the performance of the currently derived control pol- icy, while the PIM step (8b) refines this policy based on the evaluated value function. According to Vrabie & Lewis (209), these iterative steps converge to an optimal control policy and its corre- sponding optimal value function when computations, Le.,, the integration of the utility function {()),
In the Pl process, the PEV step (8a)
are performed flawlessly. However, any computational error Err(E(1) of the quadrature rule (3) can adversely affect the accuracy of PEV. This inaccuracy can subsequently influence the updated control policies during the PIM step. Such computational errors may accumulate over successive iterations of the Pl process and finally impact the learned controller, leading to the phenomenon of "computation impacts control" in IntRL.
The impact of computation requires more consideration for systems with unknown internal dynam- ics. When the system's 's internal dynamics is known, adaptive numerical solvers such as the Runge- Kutta series, offer an accurate solution to the CT Bellman equation (refer to Appendix C These solvers adjust their step sizes according to the local behavior of the equation they solve, promising high computational accuracy. However, this step-size adjustment mechanism becomes a limitation in unknown internal dynamics cases. For instance, the Runge-Kutta (4,5) method produces both fourth- and fifth-order approximations, and the difference between the two approximations refines the step size. In real-world autonomous systems, state samples are typically acquired from sensors at evenly spaced intervals, challenging alignment with adaptive solver steps. Hence, adaptive numer- ical solvers are not suitable for unknown internal dynamics. In such cases, approximation methods such as quadrature rules become essential, but might introduce significant computational errors in the PEV step. This fact underscores the importance of considering the computational impact control performance when facing unknown internal dynamics.
Given the preceding discussion, it is evident that computational errors significantly influence control performance by affecting the convergence behavior of PI. With this understanding the backdrop, this paper aims to examine the convergence rate of PI for IntRL in the presence of computational errors. Specifically, we seek to address two key questions: (1) How does the computational error Er)(((1) in the PEV step of PI impact the IntRL algorithm's 's convergence? (2) How to quantify the computational error Err(E ((()) in the PEV step and what is the optimal quadrature rule?
THEORETICAL ANALYSIS
3.1
CONVERGENCE ANALYSIS OF PI
Newton's Method
((++)) = pio (or00--((vv)) d'(vin)"" ve (10)
PIM
PEV
y000
y(2)
v" satisfies: G(V') =0
Newton's Method with Bounded Emor
p(-+)) pco C(Ppon" pw g00 (11)
|| (m)) V
2
Mr.
PIM
PEV
0000
g00)
(00)
F S E
Figure 2: Relationship between PI and Newton's 's method. The standard PI can be regarded as performing Newton's 's method to solve the HJB equation, while Pl incorporated by the computational error can be seen as Newton's 's method with bounded error.
In this subsection, we answer the first question, i.c., examine the impact of computational error the PEV step on the convergence behavior of IntRL. Initially, we utilize the Fréchet differential to demonstrate the correspondence between the standard PI process in IntRL (without computational error) and Newton's method within the Banach space to which the value function belongs.
4

---

Published as a conference paper at ICLR 2024
Lemma 1 (PI as Newton's Method in Banach Space). Consider Banach space V (V(z)VV(>) (()) 0), equipped with the norm lla. For G(V) as defined its Gâteaux and Fréchet derivatives at V can be expressed as
Ow 1 Jw
OV
G'(V) dx
2 dx Tx T9 1og9R111
9
Consider the value function VLi at the thh iteration of Pl, executing PIM and PEV (aa) itera- rively corresponds to Newton's 's method applied to solve the HJB equation This Newton method iteration is formulated as:
v(+1)= =V() - [[][v]] G(V()).
(10)
The proof of the lemma is presented in Appendix D In this lemma, we demonstrate that a single iteration consisting of both the PIM and PEV in IntRL is equivalent method for solving the HJB equation in a Banach space. This parallel allows us to view the com-
iteration of Newton's
putational error introduced in the PEV step of the 7th iteration an extra error term, denoted E
in Newton's method iteration. We define the value function with the incorporated error term E (highlighted in blue) as V(O) Its iterative update is given by:
((+1) p( - [V())))] G((())+B) E v(0) v(())
(1I)
Assume that the error term ( has a bounded norm, denoted as | I C, where the bound ē is proportional to the computational error introduced during the PEV step, which will be shown later in Section 3.3 The subsequent theorem describes the convergence behavior of Newton's 's method when an extra bounded error term is included.
Theorem (Convergence of Newton's Method with Bounded Error). Define Bo = {VV V IIV y(t) S ro). B (V V IVV V lle S ro d} with d > o. To = [VV l(v and Lo = supver, [ GVVVilla If we assume
the following conditions: (i) G is twice Fréchet continuous differentiable B: (ii) roLo };; (iii) GVV is nonsingular for VV E B; (iv) There exists number & M > O such that
+ 2 suPvEB [IG'(V)) M 2 supvEB [V)(())ll and + Mro Verr>V2MM00 2M (v) > T Mrooo
20
Then we have (0() E B and
224
2 22(2roLo)
||((0 V S
+
(12)
① Mro
Lo
where 5 is optimal value function.
The proof can be found in Appendix ton's method, when augmented with an extra bounded error term, hinges on two key elements: the cumulative effect of the bounded eтгоr term over the iterations and the inherent convergence prop- erties of the standard PI algorithm. It should be noted that the convergence property described in Theorem 2 is local and applies when the underlying operator G(V) is twice continuously differen- tiable. This iteration commences with an estimate of the constant d to ascertain the value of B. If condition (iv) of the theorem is not met, we are likely encountering pathological situation, which necessitates choosing a smaller d until the condition (iv) is satisfied or (iv) and (v) are found to be incompatible in which case our analysis breaks down. This theorem extends existing insights from earlier work on the error analysis of Newton's 's method mrabod (1956), which assumes striet Lips- chitz condition for the operator G. However, this Lipschitz condition can be difficult to verify for utility functions that do not include time discount factor. A comprehensive representation of the convergence analysis conducted in this subsection is provided in Figure 2
B Theorem 2 clarifies that the convergence behavior of New-
3.2
COMPUTATIONAL ERROR QUANTIFICATION
After analyzing the convergence property, our subsequent focus is quantifying the computational error inherent to each PEV step of IntRL. As discussed in Section [I] when internal dynamics is unknown, the integral of the utility function E(l) (2) in the PEV step must be approximated using state samples in discrete time. This approximation will inevitably introduce computational integra- tion error, especially when dealing with sparse samples, such as those collected from sensors on
5

---

Published as a conference paper at ICLR 2024
True Integral
50)
Computatioal Error
for e 9i
Error Bound
Br ((())) -(0) t0)
Mibe. =06))
Approximate
Quadrature Rule
( = Peet)
Bayesian Quadrature
Wiener Kernehtz O(V-3)
(( My ((ttt))(xtt)))
Optimal
Posterian:g() - -MM (Ptto) Fsvo)
Prior: ~ (P(0.))
Matárr Kemnke (N--P)
(smoothness parameter b)
Figure 3: Flowchart illustrating the quantification of computational error. The computational error defined as the absolute difference between the true integral and its approximation derived from the quadrature rule. The computational eтror is bounded by the product of the integrand's norm in the RKHS, and the worst-case error. When employed as the quadrature rule, BQ minimizes the worst-case error, which coincides precisely with BQ's posterior covariance.
real-world autonomous systems. A tight upper bound of the computational error, given a specific quadrature rule É()), can be determined by the Bauchy-Schwarz inequality, provided that / belongs the RKHS Hx Kanagawa al. (2016) 2020); Briol et al. (2019):
Err (É(I) I Hx <(Hc).
e(Hxx)
sup LEHx I S1
Err (É(I):
(13)
Here, HK represents the RKHS induced by the kernel function K with ||-1 being its correspond- ing norm, and e(Hxx) is the worst-case error. From [33. it is clear that if quadrature rule satisfies
e(Hxx ON with b > then the computational error satisfies Err (É(I) - O(N-t). Next,
we will discuss which quadrature rule minimizes the worst-case error.
Minimum of the worst error BQ: The optimal quadrature rule that minimizes the worst-case
error e(Hx) ) is the maximum a posteriori (MAP) estimate of BQ Briol et al. (2019). Specifically, the minimization can be achieved when the kernel function of BQ is set as the kernel function K that induces RKHS Hk. BQ aims to approximate the integral {()) by specifying a Gaussian process prior to the integrand (in this case, the integrand is the utility function 1) denoted as GP(0, K)
and conditioning this prior on the integrand evaluated at the discrete time (l(r(44), u(x((,)) 1= to
obtain the Gaussian posterior E(T) ~ N (ue(n) Σgmn O'Hagan (1991) Karvonen & Särkkä 201 Cockayne al. (2019); Briol et (019 Hennig al. 2022
++ AT
Ks, ds KTT((T))
(E(L)
H+T
I+AT
Σg(D
[K(,,6') K(s,T)K s')] ds ds
Here, T' ) represents the set of time instants, and ((T') corresponds to the evaluated values of the utility function at these sample points, with (I(T)) ((()));=1l =0 ((((((t)),((((())))))) being the evaluation at the 1th sample point. addition, (K(,T)); == K(,,,, (4) and (KTT); == K ((,,,))) 1,). Note that the posterior mean (e(!) takes the form of a quadrature rule, with weights wi being the elements of the vector ft ++TT ((s, ds KTJ. The BQ estimate of {()) be obtained using the MAP estimate, i.e., the posterior mean (l) (g(()- In this case, the worst-case error equals the square root of the posterior variance e(Hx) = Eg(t) Kanagawa et al. 2016 2020); Briol et al. (011)).
BQ with Wiener Kernel (Trapezoidal Rule) When the kernel function of the BQ is chosen as the Wiener process (we call Wiener kernel for simplicity), i.e., Kwiener (s s') : max {s, s"
with " t, the posterior mean and covariance of BQ are calculated as Hennig et al. (2022):
N ((((t)), (((t)),((((())))) + (x((4+1), u(x(te+1))) 8,, Eg(y
1 N
83.
(4E(1)
2
12
i=1
Here, 8, := t
t. The BQ estimate is equal to the integral approximated by the trapezoidal rule
Sul din 1959 Hennig et al. (2022). For data points evenly spaced at the interval, i.e., ò; N-I- the
AT

---

Published as a conference paper at ICLR 2024
posterior variance of BQ becomes Ee(()
12(N-1j2: implying linear convergence rate O(N-1) of the computational error bound when / belongs to the RKHS induced by Wiener kernel. However, this might not represent a tight upper bound, since the convergence rate of the trapezoidal rule can achieve O(N when the integrand function is twice differentiable Atkinson 1191]).
AT
BQ with Matérn kernel (In Sobolev Space) The Sobolev space W2 of order b E N is defined by W == {F E 1/2 : DaF E La exists V/ai S b). where o == (01, Qo,,... Odd) with og 2 0 is multi-index and |a| Di= O'k Adams & Fournier (2003). Besides, DaF is the ah weak
Σ jalso [ The Sobolev space W is the RKHS with the reproducing kernel K being the Matérn kernel with smoothness parameter b SSeeerr (2004); Matérn (2013); Kanagawa et al. (2020); [Pagliana et al. (2020). More details about the Matern kernel can be found in Appendix [F] In addition, the computational error for the optimal quadrature rule, i.e., BQ with Matérn kemel, achieves ONN
derivative of FF and its norm is defined by ||Fllwg
>>>)).. if E W Novak (2066) Kanagawa
et al. (2020). The illustration of BQ can be found in Appendix G and the technical flowchart of
computational error quantification can be summarized in Figure 3
3.3 CONVERGENCE RATE OF INTRL FOR DIFFERENT QUADRATURE RULES
After illustrating how to quantify the computational error, will delve deeper into understanding its implications on the solution of PEV (8a). And, more importantly, we aim to uncover its impact on the convergence of PI. Consider the case where the value function is approximated by the linear combination of basis functions Vrabie & Lewis (2009); Vamvoudakis et al. (2014):
v (x) = ( ((x), (xx) = (0) (z), (x)), Ono (x))
(14)
Here, d(x) represent the basis functions with e = 0, Ak = 1,2, .... ng. and w E Re denotes the parameters determined at the 2th iteration. A typical assumption of value function and the basis function is expressed as Vrabie & Lewis (2009): Vamvoudakis et al. (2014):
Assumption 1 (Continuity, Differentiaaiility and Linear Independence). The value function obtained at each PEV iteration (8a). is continuous and differentiable across D, denoted by vii C (2)). Besides, (x) E CC (2)), Vk- 1,2,..., Vag and the set {ox(x)}²' are linearly independent.
Such an assumption commonly encountered in optimal control Lewis et al. (1998); Abu-Khalaf & Lewis (2005); Vrabie & Lewis (2009); Vrabie et al.] (200)) suggesting that the value function can be uniformly approximated by when its dimension nis approaches infinity. In our analysis, we focus primarily on the impact of computational error. To isolate this effect, we make the assumption that the learning errors introduced by the approximation in (14) can be neglected. This assumption based on the premise that such errors can be effectively minimized or rendered negligible through adequate training duration and optimal hyperparameter tuning, among other techniques. With this approximation
[4], the parameter (+)) is found via the linear matrix equation:
E(Ti, (()) E(T2, T3 20)
(x((T)) 0 (x((T))) (T(((T)) o (x(())
(00 = =(0 =0
O
E(TT T 1- (0))
Φ (x (T(Tm)))) - Tt (x((T)))
(15)
Here, O has full column rank given approximate selections of Tí, T2. .... Im+1- This condi- tion can be seen as the persistence condition Wallace & Si (023). ensuring that equation (15)
yields a unique solution. For more discussions, see Appendix H Besides, E(Te, T&+1, (t)
TT++2 ((x(s), (z(s)) ds represents the integral over the time interval from Tr Tk+1 of the
utility function l. As previously discussed, achieving the exact value of this integral is not feasible. is instead approximated using quadrature rule, with its computational error bounded by (33. We symbolize the computational error's upper bound for E(TK,, T ()) as 5T Th+1, ()). Therefore, the integral value for E(T). T2+1,
lies within the interval [É(Tk. Te41, ())
og,T,, Tk+1, (8) E(T,,, Ta+1, (() + (E Th T (.)) Here, É(T,, Tk+1, (() signifies the in-
tegral's estimate using the quadrature rule. For simplicity, this can be formulated in matrix form
20)) E [50 5=() (0 + 8=(0). Consequently, due the computational error, the actual linear
7

---

Published as a conference paper at ICLR 2024
matrix equation solve in the PEV step becomes 166 instead of (15)
c - =
(16)
Here, 6 is the data matrix constructed by the states generated by the control policy a() == - g(x) T V V ()(i))) Incorporating the convergence analysis of Pl in Section 3.1 with the com- putational error quantification for each PI step in Section 3.2, we can present our mai results.
Theorem 3 (Convergence Rate of IntRL Concerning Computational Error), Assume that Q is a bounded set, and the conditions outlined in Theorem hold. Under Assumption]] we have:
2
2 (2rg
- (x)) Φ - -ro
+
Vx E E..
(17)
Lo
is the upper bound of the extra error term E(d) in (IT), satisfying E S.. where |10||00 == max sasis supren lox(x)), Besides, the definitions of M, ro, Lo consistent with Theorem
Here, = ||ollo sup; [1(((TTT)) -6 11a100(0)12 (1
83() is a vector with its htt row defined as 8E(T,, (TTT,TTTI, T a()), serving as upper bound for the computational error of the integral E(T), Tk+,, ((0)= = Sr. [[ M(r(s), al) (x(s))) By invoking Theorem 3 we can relate the com- putational error 6E(Tk. Tk+1, (()) to the upper bound of the extra error term È in the modified Newton's 's method (TT). This relationship establishes the convergence rate of IntRI., Combining with prior discussions of the computational error for the trapezoidal rule and the BQ with Matérn kernel, we can end our analysis with the subsequent corollary:
The proof is given in appendix In 17).
Corollary (Convergence Rate of IntRL for Trapezoidal Rule and BQ with Matern Kernel). As- suming that the sample points are uniformly distributed within each time interval [Tk. Tg+11 for k = 1,2. m, and that the IntRI. algorithm has run for sufficient iterations (i.e., i + 0o). Under the conditions in Theorem if belongs to the RKHS induced by the Wiener kernel and the integrals in the PEV step are calculated by the trapezoidal rule, we have (T)(00)()) V*(x) Here, (1x)) (x) is defined as (x) := lim (x))). If belongs to the RKHS W induced by the Matérn kernel and the integrals in the PEV step are calculated by the BQ with Matérn kernel
O(N-2).
(having smoothness parameter b), we have (1(20) (x) - V*(x)] (=]]===000)) O(N-))
4 EXPERIMENTAL RESULTS
In this section ! we validate the convergence rate of the IntRL for the trapzoidal rule and BQ with Matérn kernel as demonstrated in Corollary [
Example 1 (Linear System): First, we will consider the canonical linear-quadratic regulator prob- lem defined by system dynamics - = Ax + Bu and the utility function 1(x, u) x Qxx u Ru. The associated Pl can be obtained by setting <xx) = x & x. The parameter of value function (0) equals vec(P.), where P, is the solution to the Lyapunov equation Vrabie et al. (2009). In this case, the optimal parameter of the value function wt vec( can be determined by solving the Ricatti equation Vrabie et al. (2009). Consider the third-order system Jiang & Jiang (2017):
0
0
A =
0
0
B
=
0
Q = Iaxs, R = 1.
00.1
-0.5
-0.7
1
2.36
2.24
0.90
The solution of the Ricatti equation p
2.24
4.24
1.89
The initial policy of the PI is
0.90
1.89
1.60
chosen admissible policy u -ox with Kg
0,. We the trapezoidal rule and
the BQ with Matérn kernel (smoothness parameter b
4) for evenly spaced samples with size
5 < N < 15 to compute the PEV step of IntRL. The m2-norm of the difference of the parameters
( w"Il concerning the sample size can be shown in Figure 4 where ((()) == lim a00)
the learned parameter after sufficient iterations and w is the optimal parameter.
'The code is available at http://scceeeeeeeooooooieeeeeeee
8

---

Published as a conference paper at ICLR 2024
x10
x10-5
12
2
0.01
0.02
0.03
0.04
0.0000
0.0005
0.0010
0.0015
1/N²
1/V
(a) Trapezoidal Rule
(b) BQ with Matémn Kernel
Figure 4: Simulations for Example 1. The convergence rates of the learned parameters ((o) solved by the trapezoidal rule and BQ with Matérn Kernel (b 4) are shown to be O(N and ONN-+)
We demonstrate that the convergence rate of the parameter (=)) aligns with the derived upper bound for the value function outlined in Corollary [
This observation is sufficient to validate Corollary
since the convergence rate of (()) inherently serves as upper bound for the convergence rate of the value function (=c), Specifically, we have:
(x)0 - (x)) [|(()(x=) - "-"=""")""" (0())l2 S [|,( - w" (121)(()))12 Var E Ω.
Example (Nonlinear System): Then we consider a canonical nonlinear system in the original paper of IntLL (Example in Vrabie & Lewis (2009)):
1 = -1 + 32, d 0.5 + x2) + 0.5zz sin (x1) + sin(x1)u,
with the utility function defined as l(x, u) = 2 + 23 + u². The optimal value function for this system V (x) 0.5.f + x² and the optimal controller is (x))
sin(x1 )a-2. Due to the simple formulation of the value function, we set the basis function to be <(z) = x @ æ as in Vrabie & Lewis 2009). The initial policy is chosen as (0) - 1.5x, sin(x) (x) + x2) and we use the trapezoidal rule and the BQ with the Matérn kernel (smoothness parameter b 4) for evenly spaced samples with size < < 15 to compute the PEV step. The convergence rate of the learned parameter for Example 2 can be shown in Figure 5 and more simulation results regarding the convergence rate of the controller and the average accumulated cost be found in Appendix
x -10-5
x10-6
12
3
0.01
0.02
0.03
0.04
0.0000
0.0005
0.0010
0.0015
1/N2
1/M
(a) Trapezoidal Rule
(b) BQ with Matém Kernel
Figure 5: Simulations for Example 2. The convergence rates of the learned parameters (oo) solved by the trapezoidal rule and BQ with Matérn Kernel (b 4) are shown to be O(N- and O(N -.).
5 CONCLUSION AND DISCUSSION
This paper primarily identifies the phenomenon termed "computation impact control". Yet. is pos- sible to harness computational uncertainty, specifically the posterior variance in the BQ, to devise "computation-aww awaree control policy. One approach is to view computational uncertainty with pes- simism. In situations where computational uncertainty poses significant challenges, robust control methodologies can be employed to derive conservative policy. Alternatively. a more optimistic perspective would embrace computational uncertainty as opportunity for exploration.

---

Published as a conference paper at ICLR 2024
ETHICS STATEMENT
Our research strictly adheres to ethical guidelines and does not present any of the issues related to human subjects, data set releases, harmful insights, conflicts of interest, discrimination, privacy, legal compliance, or research integrity. All experiments and methodologies were conducted in simulation environment, ensuring no ethical or fairness concerns arise.
REPRODUCIBILITY STATEMENT
The source code associated with this paper is available for download through anonymous GitHub repository at wutttpss//gituub.com https///gihuu..ooo/aaaamao wttpa/aoymiyy888oompaaaa
Impacts-Contree All simulation re- sults have been uploaded to the repository, ensuring that readers and researchers can reproduce the findings presented in this paper with ease. Besides, all the theoretical proof's for the theorems and lemmas are provided the appendix.
REFERENCES
Murad Abu-Khalaf and Frank I. Lewis. Nearly optimal control laws for nonlinear systems with saturating actuators using a neural network hjb approach. Automatica, 41(5))77--791, 2005.
Robert A Adams and John JF Fournier. Sobolev Spuces (Pure and applied mathemutics; 140).
Elsevier, 2003.
Kendall Atkinson. An introduction to numerical analysis. John wiley & sons, 1991.
Leemon C Baird. Reinforcement learning in continuous time: Advantage updating. In Proceedings of 1994 IEEE International Conference Neural Networks 94), volume 4, pp. 2448- 2453. IEEE, 1994.
Randal W Beard, George Saridis, and John T Wen. Galerkin approximations the generalized hamilton-jacobi-bellaaa equation. Automatica, 33(12):2159-2177 1997.
François-Xavier Briol, Chris J Oates, Mark Girolami, Michael A Osborne, and Dino Sejdinovic.
Probabilistic integration: a role in statistical computation? 2019.
Sébastien Bubeck, Varun Chandrasekaran, Ronen Eldan, Johannes Gehrke, Eric Horvitz, Ece Ka- mar, Peter Lee, Yin Tat Lee, Yuanzhi Li, Scott Lundberg, et al. Sparks of artificial general intelligence: Early experiments with gpt-4. arXiv preprint arXiv:2303.12712. 2023.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Jon Cockayne, Chris J Oates, Timothy John Sullivan, and Mark Girolami. Bayesian probabilistic numerical methods. SIAM review, 61(4):756-789 2019.
Kenji Doya. Reinforcement learning in continuous time and space. Neural computation, 12(1): 219-245, 2000.
Philipp Hennig, Michael A Osborne, and Hans P Kersting, Probabilistic Numerics: as Machine Learning. Cambridge University Press, 2022.
Computation
Samuel Holt, Alihan Hüyük, Zhaozhi Qian, Hao Sun, and Mihaela van der Schaar. Neural laplace control for continuous-time delayed systems. In International Conference Artificial Intelli- gence and Statistics, pp. 1747-1778. PMLR, 2023.
Yu Jiang and Zhong-Ping Jiang. Robust adaptive dynamic programming. John Wiley Sons, 2017.
Motonobu Kanagawa, Bharath K Sriperumbudur, and Kenji Fukumizu, Convergence guarantees for kernel-based quadrature rules in misspecified settings. Advances in Neural Information Process- ing Systems, 29, 2016.
Motonobu Kanagawa, Bharath Sriperumbudur, and Kenji Fukumizu. Convergence analysis of de- terministic kernel-based quadrature rules in misspecified settings. Foundations of Computational Mathematics, 20: 155-194, 2020.
10

---

Published as a conference paper at ICLR 2024
Leonid Vital 'viihhvichhh Kantorovich and Gleb Pavlovich Akilov. Functional analysis. Elsevier, 2016.
Toni Karvonen and Simo Särkkä. Classical quadrature rules via gaussian processes. In 2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP), pp. 1-6, IEEE, 2017.
Peter Lancaster. Error analysis for the newton-raphson method. Numerische Mathematik, 9(1): 55-68, 1966.
Jae Young Lee, Jin Bae Park, and Yoon Ho Choi. Integral reinforcement learning for continuous- time input-affine nonlinear systems with simultancous invariant explorations. IEEE Transactions on Neural Networks and Learning Systems, 26(5):916-953) 2014.
Frank L. Lewis and Draguna Vrabie. Reinforcement learning and adaptive dynamic programming for feedback control. IEEE circuits and systems magazine, 9())322-50, 2009.
FW Lewis, Suresh Jagannathan, and Aydin Yesildirak. Neural network control of robot manipulators and non-linear systems. CRC press, 1998.
Timothy P Lillicrap, Jonathan Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra. Continuous control with deep reinforcement learning. arXiv preprint arXiv:1509.029) 2015.
Michael Lutter, Christian Ritter, and Jan Peters, Deep lagrangian networks: Using physics as model prior for deep learning. In International Conference on Learning Representations, 2018.
Bertil Matérn. Spatial variation, volume 36. Springer Science & Business Media, 2013.
Hamidreza Modares and Frank L Lewis. Optimal tracking control of nonlinear partially-unknown constrained-inpel -nnut systems using integral reinforcement learning. Automatica, 50(7):180--179, 2014.
Hamidreza Modares, Frank L Lewis, and Mohammad Baaaheead-Bagher Naghibi -istani. Integral rein- forcement learning and experience replay for adaptive optimal control of partially-unknon constrained-input continuous-time systems. Automatica, 50(1)2193-202, 2014.
Erich Novak. Deterministic and stochastic emoor bounds in numerical analysis, volume 1349. Springer, 2006.
Anthony O' D'Hagan. Bayes- hermite quadrature. Journal of statistical planning and inference, 29(3): 245-260, 1991.
Alexander M Ostrowski. Solution of equations and systems of equations: Pure and applied mathe- matics: A series of monographs and textbooks, vol. 9, volume 9, Elsevier, 2016.
Nicolò Pagliana, Alessandro Rudi, Ernesto De Vito, and Lorenzo Rosasco. Interpolation and learn- ing with scale dependent kernels. arXiv preprint wxXi:22006.09984 2020.
Julian Schrittwieser, Ioannis Antonoglou, Thomas Hubert, Karen Simonyan, Laurent Sifre, Simon Schmitt, Arthur Guez, Edward Lockhart, Demis Hassabis, Thore Graepel, et al. Mastering atari, go, chess and shogi by planning with a learned model. Nature, 588(7839):604 609, 2020.
Matthias Seeger. Gaussian processes for machine learning. International journal of neural systems, 14(02):69-106, 2004.
David Silver, Aja Huang, Chris Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schritwieserr Joannis Antonoglou, Veda Panneershelvam, Marc Lanctot, al. Mastering the game go with deep neural networks and tree search. nature, 529(7587):484 489, 2016.
David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Mare Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, et al. Mastering chess and shogi by self- self-play with general reinforcement learning algorithm. arXiv preprint 172 2017.
Al Al'bert Valentinovich Sul'din Wiener measure and its applications to approximation methods. Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika, (6):145-158, 1959.

---

Published as a conference paper at ICLR 2024
Richard S Sutton, Andrew G Barto, et al. Introduction to reinforcement learning, volume 135. MIT press Cambridge, 1998.
Minoru Urabe. Convergence of numerical iteration in solution of equations. Journal of Science of the Hiroshima University, Series A (Mathematics, Physics, Chemistry), 19(3):4794489 1956.
Kyriakos G Vamvoudakis and Frank L. Lewis. Online actor-critic algorithm to solve the continuous- time infinite horizon optimal control problem. Automalica, 46(5):878-888, 2010.
Kyriakos G Vamvoudakis, Draguna Vrabie, and Frank L. Lewis. Online adaptive algorithm for opti- mal control with integral reinforcement learning. International Journal of Robust and Nonlinear Control, 24(17):2686-2) 2710,, 2014.
Draguna Vrabie and Frank Lewis. Neural network approach to continuous-time direct adaptive optimal control for partially unknown nonlinear systems. Neural Networks, 22(3):237-246, 2009.
Draguna Vrabie, O Pastravanu, Murad Abu-Khalaf, and Frank L Lewis. Adaptive optimal control for continuous-time linear systems based on policy iteration. Automatica, 45(2):477-484, 2009.
Brent A Wallace and Jennie Si. Continuous-time reinforcement learning control: A review of theo- retical results, insights on performance, and needs for new designs. IEEE Transactions on Neural Networks and Learning Systems, 2023.
Haoran Wang and Xun Yu Zhou, Continuous- timee mean- variance portfolio selection: A reinforce- ment learning framework. Mathematical Finance, 12(3)):273-1308 2020.
Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
Cagatay Yildiz, Markus Heinonen, and Harri Lähdesmäki. Continuous-time model-based reinforce- ment learning. International Conference on Machine Learning, pp. 12009-12018. PMLR, 2021.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. In International Conference Learning Representations, 2019.
12

---

Published as a conference paper at ICLR 2024
APPENDIX
A MOTIVATING EXAMPLE FOR KNOWN INTERNAL DYNAMICS
The phenomenon "computation impacts control performance" also exists for the CTRL algorithm with known internal dynamics Yildiz al. (2021). Figure/6 we evaluate the impact of different numerical ordinary differential equation (ODE) solvers on the accumulated costs when applying the CTRL algorithm with known internal dynamics. Specifically,
consider the fixed-step Euler
method with varying step sizes and the adaptive Kunge-Kutta (4.5) method. Our results from the Cartpole task show that the adaptive ODE solver, in this case the Runge- (4.5) method, consis- tently outperforms the Euler method, yielding lower accumulated costs. Additionally, smaller step size is associated with improved performance within the Euler method.
Euler Method (0.01)
Accumulated Cast
Euler Method (0.02)
10
Euler Method (0.04)
RK45
10²
10"
10"
0.0
20.0
25.0
Thousand lteration
Figure 6: The illustration of the phenomenon "computation impacts control performance" when internal dynamics is known. This figure shows the accumulated costs for the Euler method with different step sizes (indicated in parentheses) and the Runge-Kutte (4,5) method in the Cartpole task with known internal dynamics.
B DEFINITION OF ADMISSIBLE POLICIES
Definition 1 (Admissible Policies Vrabie & Lewis (2009). A control policy u is admissible with respect to the cost, denoted as u E A(2)), if u is continuous on Ω u(0) = 0, stabilizes the system Qt and J (Io, u) in (5) is finite for all 2o Ω.
C CT BELLMAN EQUATION FOR KNOWN INTERNAL DYNAMICS
When the internal dynamics is considered to be known, which can be derived from fundamental principles such Lagrangian or Hamiltonian mechanics, inferred using deep learning models Chen al. (2018), Lutter al. (2018), and Zhong al. (2019). In this case, the CT Bellman equation can be addressed by tormulating an augmented ODE:
[
=(t) [((()) +(0)))())))0)))))) (x((t)
((((),((((())
Here, both the state and value function V act independent variables of the augmented ODE. This can be solved by numerical ODE solvers ranging from fixed-step solvers like the Fuler method to adaptive solvers such as the Runge-Kutta series Yildiz et al. (2021). The computational error for solving the CT bellman equation can be sufficiently small by choosing the adaptive ODE solvers.
13

---

Published as a conference paper at ICLR 2024
D PROOF OF LEMMA I
Proof. For AV V and W Vn C CV,, where v, denotes the neighborhhod V. By the definition of G(V) in 7. we have
G(V + sW) G(V)
a(V + sW)
1 a(V + +W)
dr g 0(V +
&W) 2v
f + 1 OV
g g OV
dr
4
dr
dr
4 dx
dr
OW
1 OW
9 T oW
OW
0V
f - dr
T
g
dr 2 du 9 T9R dr
0r
Thus, the Gâteaux differential Kantorovich & Akilov 2016) V is
L(W) GV(V)W
G(V + +sW G(V)
lim
s-30
8
OW
1 oW
0V
2 dr gl "g" Ōr
Ox
Then we will prove that L(W) is continuous on Vn. For WWo E Vn, we have I(W)- I(Wo)
[C' GVV
8(W - WWo)
O@( - WW)
OV
gR- 9
Or
Ox
Ör
[I.fll +
2 2
1
UV
8(W Wo)
dx
Oz
1
0V
IIfla +
79 g
m1
|| - Wolla.
dx
where m, 0. Then, for Ve >> 0, there exists ō mi a+ 5aR 1 2V lln)
satisfying |[W
L(Wo)lle when I Wo < Due to the continuity, L. is also the Fréchet derivative at V
Kantorovich & Akilov (0116).
Then, prove PI defined in 8a 8b
equals 10
From 10 obtain
G'(V (r
(1(())))) G(((O).
(18)
According to and 8b
we have
G((()))( @v((-1
f (1((+))
2 Ox gR- -V(())
OV((
dx (19a)
dr
@v(-11)
(i+ (+ggit gu((t)))
dr
aV()
2 d. gR OV()
(V(()
G'(V = dx
(19b)
dr 2u (+1) Ruti+1)
dx
f + q 4 Ox. gl T Ov(1)
10v(1
G(V()
(19c)
Ox
Or
q - (i+1) T (1)
Ox
Pluging 19 into
18
we obtain
(v
(f + ((((t())) q- ((++)) (+(++))
dr
Which means
f(i+1) = = - (i+1) Rit
(20)
Finally, (aa) can be obtained by integrating 20
n
14

---

Published as a conference paper at ICLR 2024
E PROOF OF THEOREM 2
First, we introduce useful lemma for the convergence of the standard Newton's method (without computational error).
Lemma 2 (Convergence of the Standard Newton's Method Ostrowski (0116)
Lancaster (966)).
Suppose these exists an v V for which [G'(v["))]
-1 exists and G is twice Fréchet differen-
tiable on Bo. where Bo = {V E V |V - lle rro.. Define ro := =[[()) V( lla = Il [G' (V[(V)))) lla and Lo s, [V'(V()))] (V) we have ro S 2 then the sequence v E Bo generated by Newton's method converges to the unique solution of G(V) = 0, i.e., V* in Bo. Moreover, we have [|[(0- V* < 22ro20))))"
Le
Then we will analyze the convergence property of Newton's method for solving G(V) = 0 with an extra error term E(" < in (II) We first study the iteration error between the value function for standard Newton's 's method V(a) and the value function for Newton's 's method incorporated with an extra error term (0) Specifically, we have the subsequent proposition:
Proposition 1 (Iteration Error for Newton's 's Method). If (0 B. G" exists and continuous in B, where = (V E V IIV V I|a I ro + d} with d 2 0. If GV((i) has an inverse then
((41) 1100--111 (((11) [c((11)]]
{1re
pC0
-"(vv--V,,) dV + f (b( EX0,,
(21)
where hi - [V'( G(V()). -(((i)))
Proof. Subtract (10 with
we have
v(i+1) - ((+11) = v(i - pr()) [cr G(V())+ [c](p(]]]]] G(D(0)- - E(()
= [c((] -1 [C000))+(000000000_0000000000000 ((0) + - -E()
(22)
Besides, we notice that
G(V()) G(VV)) v(0)+ Lr -v.)),,
G'(DO)) G((V()) + f ("(V)dv.
pw
(23)
D(
G((( -G((()) po
G"(h,.)
By substituting 23 into 22
obtain
21
Then, we will use Proposition to achieve the bound of the iteration error [[V() ( ||. as summarized in the subsequent lemma.
Lemma 3 (Iteration Error Bound for Newton's Method), Under the Assumption in Lemma 2 and
additionally assume that G" exists and is continuous in B and G' is nonsingular for VV E B. If
there exists number 40 M >> O, such that 1 supveb [VG lla M 2 supves G (V)lka
and +- Mr > 22M &D and d >
D-Min We have V(0) E B and
20
2
|V0 - v:0
(24)
②b Mro
15

---

Published as a conference paper at ICLR 2024
Proof. The Lemma is proved by induction. Recall that V(O) (((0),
obtain
|| - 100 1|
||e
Mro
27
Mro
Thus 24
holds for i 1. Now suppose
24 is proved for i. Because d > Mro and V Bo
220
according to Lemma we have f7(1) E B. Then by setting V = v (V()) E B, we get
Live ("(V()-V,,)
p()
- | - 1 I GG + + - dr
a
1 p00 -
Denoting 5 = --Mmm
25t
by
21 we have
[ - ((+ lla
pio
V,-) + for p
("(ha-)aV } E(()
vin
12
+ +1/000 G"(h,,))
pro
[c(
c"(vVo-V,
+ E
12
17
M8 + M1thall0d + x
12
According to Lemmal
Ihilla ro. 1, 1,2,... so we have
MR 2M + 2
p(++1) - -((+2)
24
8 M8 + 2Mro + ( Mro)
2
8 M8 + Mro + +
2
The condition - Mro > 22M implies that
M8 + Mro 4.
Hence
[ ((+1)) n 5
2
क G-Mro
This completes the proof of Lemma
Using a direct combination of Lemma and Lemma 3. and perform triangle incquality in Banach space, we can obtain Theorem 2
F DETAILS FOR MATÉRN KERNEL
We focus on scale-dependent Matérn kernels, which are radial basis kernels with Fourier decay characterized by smoothness parameter b > 1/2 for samples belonging to Wendland (2004); Pagliana al. (2020):
||x - 20||
6-3
| - x
K Matörn (x,
Bo-
p
16

---

Published as a conference paper at ICLR 2024
where || zx|] is the Euclidean distance between and x', is parameter that controls the smoothness of the function, p is the length scale parameter and 3. is the modified Bessel function of the second kind with parameter a |Seeger] (2004). A specific and important property of the Matérn kernel is that for 5. the RKHS corresponding to different scales are all equivalent to the Sobolev
space W Seeger 2004);
Matérn )013): Kanagawa et al. 2020
G ILLUSTRATION OF BQ
Here, we will take the integral
[ 10
3tm t) 5
sin
dt
(25)
10
example to illustrate BQ using Wiener kernel and Matérn kemne,, see Figure 7 and respec- tively. can be seen from figures that BQ can be understood as first fitting the integrand using Gaussian process regression, and then performing integration. The accuracy of BQ is related to the choice of the Gaussian process kernel and the sample size. For evenly spaced samples, as the sam- ple size increases, the accuracy of BQ improves. Besides, we also plot the integral values and their computational errors for further illustration; see figure
GP Regression Actual Integral
GP Covariance
GP Regression
GP Covariance Sample Points
Sample Points
Actual Integral
((x()) ((x()))
1lxtO, ((x(t)))
2
2
4
6
8
10
2
4
8
10
N
6
(b) N
4
GP Regression
GP Covariance
Regression
GP Covariance
((x()) u(x(t)))
1.
Actual Integral
Sample Points
HxM.. u(x(()))
Actual Integral
Sample Points
2
2
4
6
8
10
2
10
t
(c) N 10
(d) N = 12
Figure 7: Illustration of BQ with Wiener kernel (equal to the trapezoidal rule) for the integral in 25
H DISCUSSION ABOUT FULL COLUMN RANK CONDITION
Lemma 4 (Full Column Rank of ( Vrabie & Lewis (009)). Suppose (x) A(Q), and the set { no is linearly independent. Then, 3T >> O such that for Vx(t) generated by (x)), (ob( + T)) - (x(x(())) is also linear independent.
Based on this lemma, there exist values of Tı. T2, Tm+11 such that ( is invertible and thus 15 has unique solution.
17

---

Published as a conference paper at ICLR 2024
GP Regression
GP Covariance
GP Regression
GP Covariance Sample Points
Actual Integral
Sample Points
Actual Integral
((x()), u(x(e)))
((x(()) [0x(t), u(x((t))
2
0
2
4
6
8
10
2
4
8
10
(a) N
(b)
4
GP Regression
GP Covariance
GP Regression
GP Covariance Sample Points
Actual Integral
Sample Points
Actual Integral
((x(t)) ((xttt))
2
1
2
4
6 t
8
10
2
4
8
10
(c) N 10
(d) N = 12
Figure Illustration of BQ with Matérn kernel the integral
16.0
BQ Integral
Actual Integral
0.5
Value
15.5
Integral Emor
Integral PregralValle
0.4
0.3
14.5
0.2
14.0
6
10
12
0.010
0.015
0.020
0.025
1/N²
(a) Integral value Wiener kernel
(b) Integral error Wiener kernel
16.0
Integral
Actual Integral
0.6
Integral Value
15.5
Integral Error
15.0
0.2
14.5
0.0
14.0
6
10
12
0.0002
0.0004
0.0006
0.0008
1/V
(c) Integral value Matérn kernel
(d) Integral enorr Matérn kernel
Figure 9: Integral values and errors of BQ for the integral in (5
I PROOF OF THEOREM 3
First, we will give proposition which presents the upper bound of the error between w(() and (0)) the same control policy (1) is utilized.
18

---

Published as a conference paper at ICLR 2024
Proposition 2 (Solution Error of PEV). For the same control policy ( the difference of the solu- tion of the (5) and (16) can be bounded by:
||
1: 1l2 = ( |2
Proof. For the same control policy (1)), we have ) ( Consider the linear matrix equation
15 and (6), we have
||0( - 10 -1) (E - -
s (((6(66(0)) | = ( - 1|2
< [1( 12
In this equation, (5(() is vector with its hth row defined as ö5(Tk. Tk+1, al0). Note
that 65 T4+1; a()) serves as an upper bound for the computational enror of the integral E(T,, Tk+1; (()) Sr. ((r(s), al0) (x(s)) ds.
Then we present the proof Theorem
Proof. We choose the norm in the Banach space V the infinity norm, i.e., [ ||
supzen [r)).. Because Ω is bounded and 6x C!, vkk 1.2.... ng. thus || llae = supzen (z)] and || := maxisksn, { exists. From Proposition and Cauchy- Schwarz inequality, we have
E (x) = (x)ll
< |1000 - 12))())))))
s |"6
s sup (16@0T 1 1= 12] l1(2 Yx E f.
Thus, we have
|| |xx = sup (x))
< sup {
[b |1o1
(26)
Combining 266 with 12 in Theorem
we can obtain
[|((] 1 < 2ΦE
2 2-((2roo))
Mro
Lo
which directly leads
2Φ
2 (20o2r)
(=)))) -Mr
+
Vo E Ω.
Lo
J ADDITIONAL SIMULATION RESULTS
In this subsection, we show additional simulations results for the convergence rate of the controller and the accumulated costs for Example and Example 2 in Section
Example 1:
Control Gain Matrix: For linear systems, the learned control policy adheres to (<s) (x:)
(o0) where (00) is the learned control gain matrix which is computed as Jiang & Jiang 201 7)
ke(o)) R- BBB p(a))
19

---

Published as a conference paper at ICLR 2024
x10- 8
x10
[1( K ||F
2
0.01
0.02
0.03
0.04
0.0000
0.0005
0.0010
0.0015
1/N2
1/N4
(a) Trapezoidal Rule
BQ with Matém Kernll
Figure 10: Simulations for Example showing convergence rates of ((0)) computed via the trape- zoidal rule and BQ with Matérn Kernl (b 4) as O(N-2) and O(N respectively.
Here, p(x) ensures @((x)) = vec(P(oe)). The Frobenius norm difference between the learned and optimal control gain matrix K" = [0.90 1.89 1.60], denoted || K |r is depicted in Figure I0
Average Accumulated Cost:
The average accumulated costs of the learned and optimal policies, represented as J and J respectively, with initial state Igg ~ NO, 1002. T3x3), are defined as:
Eeg (0.1000-aa0)
(1((0) (xo)) J
V00 (V (ro))
Their difference concerning the sample size N is shown in Figure 11
0.8
12
0.6
1
0.4
3
0.2
0.01
0.02
0.03
0.04
0.0000
0.0005
0,0010
0.0015
1/N2
1/N4
(a) Trapezoidal Rule
(b) BQ with Matém Kerel
Figure 11: Simulations for Example 1 illustrating convergence rates of the average accumulated cost computed via the trapezoidal rule and BQ with Matérn Kernel (b = 4) as O(N-2) and O respectively.
Example
Control Policy: For nonlinear systems, the control gain matrix does not exist. Instead, the learned control policy is given by
(0)) R- )(=) (T (o0(
Therefore, present the average control difference across state æ (O0 1002. 12x2):
Eg=~A {lat>>
uz (x))
where u² (x) represents the optimal controller. The average difference between the learned and optimal control policies is presented in Figure 12
Average Accumulated Cost:
The average accumulated costs of the learned and optimal policies,
represented and respectively, with initial state 2i ~ (O, 1002.. T2x2), are defined as:
= Egg (10.00000000)
(Zo))
J
Ee, {V' (xo))
20

---

Published as a conference paper at ICLR 2024
x10-3
x10
3
2
0.01
0.02
0.03
0.04
0.0000
0.0005
0.0010
0.0015
1/N2
1/N4
(a) Trapezoidal Rule
BQ with Matém Kemel
Figure 12: Simulations for Example displaying convergence rates ((c)) via the trapezoidal rule and BQ with Matérn Kernel (b = 4) as O(N-²) and O respectively.
Their difference concerning the sample size N is shown in Figure
x10-2
1.5
12
1.0
0.5
0.01
0.02
0.03
0.04
0.0000
0.0005
0.0010
0.0015
1/N2
1/N
(a) Trapezoidal Rule
(b) BQ with Matérn Kernel
Figure 13: Simulations for Example 2 illustrating convergence rates of the average accumulated cost via the trapezoidal rule and BQ with Matérn Kernel (h = 4) as O( and O(N-1) respectively.
K DISCUSSION ABOUT THE LIMITATIONS OF INTRL
In this appendix, want discuss the limitations of IntRL algorithm. IntRL is a canonical control task in CTRL, just like the Q learning in DTRL. However, IntRL is still underdeveloped for high- dimensional systems compared to Q learning. As shown in Vrabie & Lewis (2009); Modares et al. (2014); Wallace & (023). IntRL is typically applied systems with no more than dimensions. This is because IntRL is hard to converge in dine-siimniii systems. The difficulties in achieving convergence in high-dimensional systems be attributed several factors:
Limited Approximation Capability: IntRL uses linear combination of basis functions approximate the value function, which is less powerful than the neural network-based ap- proaches prevalent in DTRL. Besides, in Wallace Si (2023), the author observes that the condition number of the pseudo-inverted matrices for IntRI. will degrade significantly due an additional basis function. This conditioning issue limits IntRL leverage complex basis functions to obtain larger approximation capability. is emphasised in Wallace Si (2023) that severe numerical breakdowns to even small increments in problem dimension.
Lack of Effective Exploration Mechanisms: IntRL struggles with poor data distribution due to inadequate exploration strategies. Common DTRL methods improve exploration by adding noise during data collection, but in continuous systems, this leads to complex stochastic differential equations (SDEs) that are challenging to manage. In Wallace & (2023), the author emphasises that IRL's lack of exploration noise causes data quality degradation especially when the state is regulated to the origin. Besides, it is observed that lack of exploration noise will result in the phenomenon of "hyperparameter deadlock".
21

---

Published as a conference paper at ICLR 2024
Absence of Advanced Training Techniques: Techniques like replay buffers, parallel ex- ploration, delayed policy updates, etc., which enhance sample efficiency and training in DTRL, are lacking in IntRL.
No Discount Factor:
The absence of discount factor in IntRL makes it difficult to en- sure that policy iteration acts as a contraction mapping, which can affect the stability of algorithm's 's training process.
Non-existence of Q function: In CT systems, the concept of Q function is not directly ap- plicable. The absence of Q function impedes the direct translation of many DTRL method- ologies and insights to IntRL.
Considering these challenges, CTRL algorithms should be investigated and developed for complex and high-dimensional scenarios in future research.
L FURTHER DISCUSSION ABOUT THE MOTIVATION OF CTRL
Many physical and biological systems inherently continuous in time, governed by ODEs. Com- pared to discretising time and then applying DTRL algorithms, directly employing CTRL offers these ad vantages Jiang & Jiang (2017) Wallace & Si (202)):
Smoother Control Law: Direct application of CTRL typically results smoother control outputs. This contrasts with the outcomes of coarse discretization, where control is less smooth and can lead to suboptimal performance Doya (000).
Time Partitioning:
When addressing CTRL. directly, there's no need to predefine time partitioning. Instead, it is efficiently managed by numerical integration algorithms, which find the appropriate granularity Yildiz et al. (202T). Thus, CTRL often a better choice when the time interval is uneven.
Enhanced Precision with CT Transition Models: Employing CT transition models in modelling CT systems offers a higher degree of precision compared to DT transition mod- els. As evidenced in Wallace & Si (2023), comparisons between CT and DT trajectories, such as in the CartPole and Acrobat tasks, reveal that CT models align more closely with true solutions, especially in scenarios involving irregularly sampled data.
While direct comparisons of performance between DTRL CTRL in literature are rare, CTRL's prominence in the financial sector, especially in portfolio management, is notable. Many portfolio management models inherently rely Stochastic Differential Equations (SDEs). A notable exam- ple includes the comparison in Wang & Zhou (2020) between the EMV algorithm (a CTRL-based approach) and DDPG (a DTRL approach) Lillicrap et al. (2015). This comparison showed the EMV algorithm's 's superior performance, highlighting CTRL's advantages in scenarios that utilize SDE- based models.
The distinct benefits of CTRL. across various domains, from its enhanced precision in modeling to its adaptability handling irregular time intervals, establish a vital and influential methodology. Its applicability diverse sectors, notably in complex and dynamic fields like engineering and finance, make CTRL increasingly relevant and powerful tool.
22


Reviews:
[{"review_id": "gx41VRsL4T", "reviewer": "Reviewer_8W8E", "summary": "This paper studies the impact of computation methods (quadrature rule for solving integrals) when applying reinforcement learning in continuous control tasks. Building upon the connections between HJB equation and Newton's method, the authors show that the computation error is an extra error term in each iteration of Newton's method. With the bounded error assumption, they provide a convergence results for Newton's methods with an extra error term (Theorem 1). Furthermore, the computation error bounds are also discussed by minimizing the worst case error under different quadrature rules and kernels. Finally, an end to end convergence result is provided (Theorem 3, Corollary 1).", "questions": "1. Are the affine nonlinear systems necessary? Can you consider a more general class of systems?  \n2. The assumptions of Theorem 2 need more explanations. To guarantee those assumptions, what properties do you need for the systems functions $f$, $g$, and cost function $J$?\n3. The experimental examples are basically toy examples (3d linear system and 2d nonlinear system). These are far from high-dimensional real-time applications. How the computation error and convergence will behave in higher dimensional cases need to be examined.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper is well written and easy to follow. The problem of studying the impact of computational errors on continuous control in integral RL is well motivated and interesting. As the authors claim, this problem is widespread but understudied.", "weaknesses": "Although the problem studied in this work is interesting. I feel like the paper is mainly a combination of existing results (convergence of Newton's method with an extra error term and error bound on the computation step). What are the novel techniques appled in the analysis such that this work is not simply A+B?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Y2mBo85Ffb", "reviewer": "Reviewer_3VMw", "summary": "This paper studies the continuous time RL and provides a detailed convergence rate discussion on the impact of policy iteration's computation errors when approximating the integration in the policy evaluation step. By showing that PI can be viewed as Newton updates and PI with computation errors can be viewed as Newton updates with errors, this paper established local convergence rates and uses simulation to demonstrate the tightness of the order of the convergence rate.", "questions": "Q1: In Theorem 3, the computation error is treated as a constant value. How does this constant decay with the number of samples? \n\nQ2: Corollary 1 assumes that $i \\to +\\infty$. By using the decay rate of the computation error in Theorem 1, can the authors comment on a more realistic convergence rate based on different number of iterations $i$? Further, can the theoretical results provide some guidelines on how to choose the number of iterations to terminate at?\n\nQ3: How does this result compare with LSPI, which is also based on a linear combination of basis functions?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is very well written. Even though I am not an expert in this area, the detailed motivation and the clear illustration diagrams help me understand the importance of this problem and the key ideas behind the proofs. The connection between Newton updates and PI, and using approximate Newton to analyze PI with computation errors are also fascinating. Further, the proofs are quite involved too. Lastly, the numerical results demonstrate the tightness of the order of the convergence rate with respect to N.", "weaknesses": "See below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8ygSJZ27JE", "reviewer": "Reviewer_jx6Y", "summary": "The paper addresses the impact of computational methods on control performance in Integral Reinforcement Learning. The authors focus on the policy evaluation stage of IntRL, where the integral of the utility function needs to be computed using quadrature rules. They demonstrate that computational errors introduced during PEV can significantly influence the convergence behavior of policy iteration and the performance of the learned controller.\n\nThe authors show that computational errors in PEV manifest as an extra error term in each iteration of Newton's method. They provide a theoretical analysis, proving that the upper bound of this error term is proportional to the computational error. The paper further explores the case where the utility function resides in a reproducing kernel Hilbert space (RKHS), presenting local convergence rates for IntRL using both the trapezoidal rule and Bayesian quadrature with a Matern kernel.", "questions": "Could you provide more details on the choice of the canonical control tasks used for experimental validation? Were any real-world scenarios or more complex tasks considered?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The authors demonstrated how computational errors in the PEV stage of IntRL affect the convergence behavior of policy iteration and the performance of the learned controller, which is previously unexplored.\n- They also provided a solid theoretical analysis of the impact of computational errors, providing bounds and convergence rates that relate the computational method to control performance.\n- Validation of the theoretical findings is also offered through simulations on canonical control tasks, showing the practical implications of the choice of computational method in IntRL.\n\nThe paper sheds light on the impact of computational methods on control performance in IntRL, providing both theoretical insights and practical guidelines for improving controller learning in continuous-time reinforcement learning scenarios.", "weaknesses": "- The paper provides theoretical claims about the impact of computational methods on control performance in IntRL. However, the experimental validation seems to be limited in scope. The authors only consider canonical control tasks to validate their findings. The authors could consider a broader set of experiments, including more complex and real-world scenarios, to showcase the practical implications of their findings.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xYDPDKzzFA", "reviewer": "Reviewer_EcSW", "summary": "This work investigates the effects of choice of quadrature rules for integral RL; especially when the true dynamics is unknown.\nIt discusses how the computational error in policy evaluation stage affects each iteration of PI, which is shown to be corresponding to the Newton’s method: theoretically, the work proves the local convergence rates for IntRL; and the findings are validated by some control tasks.\nFurthermore, the work shows that the case where the utility function lives in an RKHS as corollary.", "questions": "1. For Appendix G, not only for the figure for the utility itself, a plot for the integral of the utility and the worst case error would be informative too.\n2. Are there any answer to the weakness points?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "1. Conceptually this work may bring up a new research direction of studying the effects of “approximation error” for ODE or problems including integrals in ML field (e.g. for neural ODE).  This is the point I particularly find value in this work.\n2. The claims are validated by simple yet informative experiments.", "weaknesses": "1. While CT formulation helps in some analysis or for some applications, it was a bit unclear what the motivations behind studying CTRL if the task can be done with DT formulations.\n(Especially when the time interval is even; in which case DT can well manage.)\nIn particular, for DT and CT systems, there should be different conditions for the solutions to exist.\nFor contact-rich dynamics for example, this kind of analysis becomes harder for example.\n(Does PIM ensures the existence of solutions throughout the whole process?)\nAlso, for stochastic systems, CT requires more conditions for certain analysis.\n2. About approximation error of value functions; if we know the utility lives in certain RKHS, can we say anything about the value function which may validate the assumptions?  At least, there should be a trivial case for this assumption: if you know the value function exactly, that becomes a single basis function; when is it continuously differentiable?  A bit more discussions needed.\n(also there is a type (additional “]”) for the interval for the integral value.)\n3. In Appendix H; no approach to find a suitable T?  There may require some discussions on how rare the independence property fails for a random T.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)