# On the Convergence of Reinforcement Learning in Nonlinear Continuous State
Space Problems

Raman Goyal, Suman Chakravorty, Ran Wang, Mohamed Naveed Gul Mohamed The
authors are with the Department of Aerospace Engineering, Texas A&M
University, College Station, TX 77843 USA. {ramaniitrgoyal92, schakrav,
rwang0417, naveed} @tamu.edu

###### Abstract

We consider the problem of Reinforcement Learning for nonlinear stochastic
dynamical systems. We show that in the RL setting, there is an inherent “Curse
of Variance” in addition to Bellman’s infamous “Curse of Dimensionality”, in
particular, we show that the variance in the solution grows factorial-
exponentially in the order of the approximation. A fundamental consequence is
that this precludes the search for anything other than “local” feedback
solutions in RL, in order to control the explosive variance growth, and thus,
ensure accuracy. We further show that the deterministic optimal control has a
perturbation structure, in that the higher order terms do not affect the
calculation of lower order terms, which can be utilized in RL to get accurate
local solutions.

###### Index Terms:

RL, Optimal control, Nonlinear systems

##  I Introduction

A large class of decision making problems under uncertainty can be posed as a
nonlinear stochastic optimal control problem that requires the solution of an
associated Dynamic Programming (DP) problem, however, as the state dimension
increases, the computational complexity goes up exponentially in the state
dimension [1]: the manifestation of Bellman’s infamous “curse of
dimensionality (CoD)” [2]. To understand the CoD better, consider the simpler
problem of estimating the cost-to-go function of a feedback policy
$\mu_{t}(\cdot)$. Let us further assume that the cost-to-go function can be
“linearly parametrized” as:
$J_{t}^{\mu}(x)=\sum_{i=1}^{M}\alpha^{i}_{t}\phi_{i}(x)$, where the
$\phi_{i}(x)$’s are some _a priori_ basis functions. Then the problem of
estimating $J_{t}^{\mu}(x)$ becomes that of estimating the parameters
$\bar{\alpha}_{t}=\\{\alpha^{1}_{t},\cdots,\alpha^{M}_{t}\\}$. This can be
shown to be the recursive solution of the linear equations
$\bar{\alpha}_{t}=\bar{c}_{t}+L_{t}\bar{\alpha}_{t+1},\text{ where
}\bar{c}_{t}=[c^{i}_{t}],$ with $c^{i}_{t}=\int
c(x,\mu_{t}(x))\phi_{i}(x)\,dx,$ and $L_{t}^{ij}=\int\int
p^{\mu_{t}}(x^{\prime}|x)\phi_{i}(x^{\prime})\phi_{j}(x)\,dx^{\prime}\,dx,\;i,j\in\\{1,\dots,M\\},$
where $p^{\mu_{t}}(./.)$ is the transition density of the Markov chain under
policy $\mu_{t}$. This can be done using numerical quadratures given knowledge
of the model $p^{\mu}(x^{\prime}|x)$, termed Approximate DP (ADP), or
alternatively, in Reinforcement Learning (RL), simulations of the process
under the policy $\mu_{t}$,
$x_{t}\xrightarrow{\mu_{t}(x_{t})}x_{t+1}\rightarrow\cdots$, is used to get an
approximation of the $L^{ij}_{t}$ by sampling, and solve the equation above
either batchwise or recursively [3, 1]. But, as the dimension $d~{}$
increases, the number of basis functions and the number of evaluations
required to evaluate the integrals go up exponentially. There has been recent
success using the Deep RL paradigm where deep neural networks are used as
nonlinear function approximators to keep the parametrization tractable [4, 5,
6, 7, 8], however, the training times required for these approaches is still
prohibitive. Hence, the primary problem with ADP/ RL techniques is the CoD
inherent in the complex representation of the cost-to-go function, and the
exponentially large number of evaluations required for its estimation. In this
paper, we show that there is an additional “Curse of Variance” that afflicts
the RL solution, the fact that the variance grows at a factorial-exponential
rate in the order of the approximation, that precludes us from solving for
higher order approximations of the feedback law. Prior research has focused in
some detail on the sample complexity of RL for finite control problems [9, 10,
11], and the case of optimal Linear Quadratic Control (LQR) in continuous
state and control spaces [12, 13]. We study the general nonlinear problem, and
show the scale of the variance inherent in an RL estimate. Albeit anecdotal
and empirical evidence of the variance phenomenon has always existed in the RL
literature [14], we believe we are the first to exactly enumerate the
factorial-exponential growth and its consequences: it is necessary that we
look for local solutions in order to find accurate solutions, that stochastic
control problems are fundamentally intractable, and the best we can hope for
is a suitably accurate deterministic approximation (see points 1-5 of
contributions below). However, this does not mean we need to give up on global
optimality. In [15], we established the local optimality of the deterministic
feedback law, in that the nominal (zero noise) action, and the linear feedback
action, of the optimal stochastic and deterministic policies are close to
fourth order in a small noise parameter, starting at any given state, which,
when allied with replanning, recovers a near-optimal solution. Thus, a local
solution allied with replanning is an efficient and near-optimal way to solve
nonlinear stochastic control problem rather than solve for a global (higher
order) solution. This gets us to the context of Model Predictive Control
(MPC).

In the case of continuous state, control and observation space problems, the
MPC [16, 17] approach has been used with a lot of success in the control
system and robotics community. However, stochastic control problems, and the
control of uncertain systems in general, is still an unresolved problem in
MPC. As noted in [16], the problem arises due to the fact that in stochastic
control problems, the MPC optimization at every time step cannot be over
deterministic control sequences, but rather has to be over feedback policies,
which is, in general, intractable. Thus, the tube-based MPC approach, and its
stochastic counterparts, typically consider linear systems [18, 19, 20] for
which a linear parametrization of the feedback policy suffices but the methods
become intractable when dealing with nonlinear systems [21]. In this paper, we
show that the stochastic problem is fundamentally intractable, and if solved
in an RL fashion, leads to a very high variance in the solution. In
particular, it becomes necessary to look for local deterministic solutions to
ensure accuracy (which are locally optimal due to the results of [15]), and
re-plan when necessary as in MPC, to recover global optimality.

We summarize our contributions as follows.

1\. It is fundamentally intractable to solve for a high order approximation of
a feedback law for optimal control via RL (global/ nonlocal), since the
variance of the solution grows factorial-exponentially in the order of the
approximation.
2\. The deterministic problem has a perturbation structure, in that higher
order terms do not affect the calculation of lower order terms, and thus, when
a model is known, the calculations can be closed at any order without
affecting the accuracy of the lower order terms.
3\. If the deterministic problem is solved in an RL fashion, then an accurate
solution can be found, if and only if we concentrate on a suitably local
solution, enforced via constraining the random exploration around a nominal
trajectory.
4\. The stochastic problem is intractable in the sense that it lacks a
perturbation structure, and if solved via RL, this necessarily implies very
high bias as well as variance, and hence, inaccuracy in the result.
5\. The perturbation structure and locality of the solution are key to an
accurate RL implementation.

Outline of Paper. The rest of the document is organized as follows: Section II
outlines the Problem Formulation, Section III studies the convergence of
Policy evaluation in a finite time RL setting and the resulting variance in
the solution. Section IV derives a perturbation structure inherent to the
deterministic policy evaluation problem, and shows how to leverage this for
accurate local RL solutions. Section V concentrates on the stochastic policy
evaluation problem and shows its fundamental intractability in terms of the
lack of a perturbation structure, and the resulting high variance for RL
solutions. Section VI gives empirical results in a simple example to validate
the theoretical development.

##  II Problem Formulation

The problem of control under uncertainty can be formulated as a stochastic
optimal control problem in the space of feedback policies. We assume here that
the uncertainty in the problem lies in the system’s process model.
System Model: For a dynamic system, we denote the state and control vectors by
$x_{t}\in\ \mathbb{X}\subset\ \mathbb{R}^{n_{x}}$ and $u_{t}\in\
\mathbb{U}\subset\ \mathbb{R}^{n_{u}}$ respectively at time $t$. The motion
model
$h:\mathbb{X}\times\mathbb{U}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{X}$ is
given by the equation

$x_{t+1}=h(x_{t},u_{t},w_{t});\ w_{t}\sim\mathcal{N}(0,{\Sigma}_{w_{t}}),$ (1)

where {$w_{t}$} are zero mean independent, identically distributed (i.i.d)
random sequences with variance ${\Sigma}_{w_{t}}$.

Stochastic optimal control problem: The stochastic optimal control problem for
a dynamic system with initial state $x_{0}$ is defined as:

$J_{\pi^{*}}(x_{0})=\min_{\pi}\
E\left[{\sum^{T-1}_{t=0}c(x_{t},\pi_{t}(x_{t}))+g(x_{T})}\right],$ (2)

$s.t.\ x_{t+1}=h(x_{t},\pi_{t}(x_{t}),w_{t}),$ where: the optimization is over
feedback policies $\pi:=\\{\pi_{0},\pi_{1},\ldots,\pi_{T-1}\\}$ and
$\pi_{t}(\cdot)$: $\mathbb{X}\rightarrow\mathbb{U}$ specifies an action given
the state, $u_{t}=\pi_{t}(x_{t})$;
$J_{\pi^{*}}(\cdot):\mathbb{X}\rightarrow\mathbb{R}$ is the cost function on
executing the optimal policy $\pi^{*}$;
$c_{t}(\cdot,\cdot):\mathbb{X}\times\mathbb{U}\rightarrow\mathbb{R}$ is the
one-step cost function; $g(\cdot):\mathbb{X}\rightarrow\mathbb{R}$ is the
terminal cost function; $T$ is the horizon of the problem. The solution to the
above problem is given by the Dynamic Programming equation:

$J_{t}(x)=\min_{u}[c(x,u)+E[J_{t+1}(x^{\prime})]],$ (3)

where $x^{\prime}\sim p(./x,u)$, and $p(./x,u)$ denotes the transition density
of the state at the next time step arising form the system dynamics, given the
control $u$ is taken at state $x$, solved with the terminal condition
$J_{T}(x)=g(x)$. The DP equation can be solved using the so-called Policy
Iteration method, where given a time varying feedback policy
$\pi^{(k)}_{t}(\cdot)$, one first solves for the cost function corresponding
to it:

$J^{(k)}_{t}(x)=c(x,\pi^{(k)}_{t}(x))+E[J_{t+1}^{(k)}(x^{\prime})],$ (4)

where $x^{\prime}\sim p(./x,\pi^{(k)}_{t}(x))$, and the above equation is
solved with the terminal condition $J_{T}^{(k)}(x)=g(x)$, which is followed by
a policy improvement step:

$\pi^{(k+1)}_{t}(x)=\operatorname*{arg\,min}_{u}[c(x,u)+E[J_{t+1}^{(k)}(x^{\prime})]],$
(5)

where $x^{\prime}\sim p(./x,u)$. This process is followed till convergence,
starting with some initial time varying policy $\pi^{(0)}_{t}(x)$ [1].

##  III Convergence of Policy Evaluation

In the following, we shall concentrate on the Policy Evaluation (PE) part of
Policy iteration, in particular, a single Policy Evaluation step, to show the
convergence issues inherent, and at the end of this section, outline the
issues arising from the dynamic recursion. We shall consider a synchronous
model of computing, i.e., where all the experiments are done first, and the
cost functions at any step updated using all the experiments. We consider the
deterministic scalar state case for simplicity, the generalization to the
vector state case is straightforward and the stochastic case is treated in
Section 5. Let us rewrite the policy evaluation equation from above as:

$J_{t}(x_{t})=c_{t}(x_{t})+J_{t+1}(f(x_{t})),$ (6)

where the deterministic dynamics are $x_{t+1}=f(x_{t})$, and given that the
terminal cost $J_{T}(x_{T})=g(x_{T})$. In the context of Policy Iteration, the
dynamics corresponds to the closed loop under some feedback policy $\pi(x)$,
i.e., $f(x)=h(x,\pi(x),0)$. Note that the policy is, in general, time varying,
and thus, the dynamics should also be time varying. But we consider time
invariant dynamics for simplicity, and all the results obtained below
generalize to the time varying case in a straightforward fashion.

#### Computing Model

Suppose that we have basis functions $\\{\phi^{1}(x),\cdots\phi^{N}(x)\\}$
such that any $J_{t}(x)=\sum_{i}\alpha_{t}^{i}\phi^{i}(x)$ for suitably chosen
coefficients $\alpha_{t}^{i}$, and such that
$c_{t}(x)=\sum_{i}c_{t}^{i}\phi^{i}(x)$. Suppose now that we are given $R$
samples from the dynamical system, say $\\{x_{1}^{(k)},\cdots x_{t}^{(k)}\\}$,
for $k=1,\cdots R,$ that are sampled from some, in general, time varying
density $p_{t}(\cdot)$. We leave the question of what this density ought to be
to later on in our development. Given the samples, we write:

$J_{t}(x^{(k)}_{t})=c(x_{t}^{(k)})+J_{t+1}(f(x_{t}^{(k)}))+v^{(k)}_{t},$ (7)

where $v^{(k)}_{t}$ is an independent identically distributed (i.i.d.) noise
sequence for all time steps $t$. Representing the cost functions in terms of
the basis functions, we obtain:

$\bar{\alpha}_{t}\phi_{t}^{(k)}=\bar{c}_{t}\phi_{t}^{(k)}+\bar{\alpha}_{t+1}\phi_{t+1}^{(k)}+v^{(k)}_{t},$
(8)

where
$\bar{\alpha}_{t}^{(k)}=\begin{bmatrix}\alpha_{t}^{1}\cdots\alpha_{T}^{N}\end{bmatrix}$,
$\bar{c}_{t}=\begin{bmatrix}c_{t}^{1}\cdots c_{t}^{N}\end{bmatrix}$,
$\phi_{t}^{(k)}=\begin{bmatrix}\phi^{1}(x_{t}^{(k)})\\\ \vdots\\\
\phi^{N}(x_{t}^{(k)})\end{bmatrix}$, and
$\phi_{t+1}^{(k)}=\begin{bmatrix}\phi^{1}(x_{t+1}^{(k)})\\\ \vdots\\\
\phi^{N}(x_{t+1}^{(k)})\end{bmatrix}$, where note that
$x_{t+1}^{(k)}=f(x_{t}^{(k)})$.

#### RL as Least Squares

Then, we may view the above as the following least squares problem:

$\bar{\alpha}_{t}^{R}=\operatorname*{arg\,min}_{\bar{\alpha}_{t}}||\bar{\alpha_{t}}\Phi_{t}^{R}-\bar{c}_{t}\Phi_{t}^{R}-\bar{\alpha}_{t+1}^{R}\Phi_{t+1}^{R}||^{2},$
(9)

where
$\Phi_{t}^{R}=\begin{bmatrix}\phi_{t}^{(1)},\cdots,\phi_{t}^{(R)}\end{bmatrix}$,
and the supercase $R$ is used to denote the solution after $R$ samples. In
this case, we are sweeping back in time starting at the final time $T$, and
thus, it is assumed above that we have solved for $\alpha_{t+1}^{R}$ already.
The solution to this problem is standard and given by:

$\bar{\alpha}_{t}^{R}=\bar{c}_{t}+\bar{\alpha}_{t+1}^{R}\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}(\Phi_{t}^{R}\Phi_{t}^{R^{\prime}})^{-1},$
(10)

where $A^{\prime}$ denotes the transpose of a matrix $A$.
Now, we shall establish some properties of the least squares (LS) solution
above in the context of RL. First, let us find the “true” solution of the
Policy Evaluation equation (6). First, we make the following assumption.

###### Assumption 1

There exist a set of constants $\beta_{ij}$, $i=1,2,\cdots N$, and $j=1,\cdots
N^{\prime}$, where $N^{\prime}>N$, such that for any
$\phi^{i}(f(x))=\sum_{j=1}^{N^{\prime}}\beta^{ij}\phi^{j}(x)$.

The reason $N^{\prime}>N$ is that, in general, unless $f(\cdot)$ is linear, it
will require more basis functions to represent $J_{t}(x)$ than $J_{t+1}(x)$.
The reason is that if $\phi^{i}$ is an $i$ degree polynomial, then
$\phi^{i}(f(x))$ will be a $ki$ degree polynomial if $f(\cdot)$ is a $k$
degree polynomial. Thus, in general, we will need an expanding basis to
represent the functions $J_{t}(\cdot)$ as we sweep back in time from $T$ to
$0$. Then, we can characterize the “true” solution to the Policy Evaluation
equation (6) as follows.

###### Proposition 1

The true solution to the policy evaluation equation (6) is given by:
$\bar{\alpha}_{t}^{*}=\bar{c}_{t}+\bar{\alpha}_{t+1}^{*}B_{t}$, where
$B_{t}=\begin{bmatrix}\beta^{11},\cdots,\beta^{1N_{t}}\\\ \ddots\\\
\beta^{N_{t+1}1},\cdots,\beta^{N_{t+1}N_{t}}\end{bmatrix}$, where $N_{t+1}$ is
the number of basis functions required to represent $J_{t+1}(\cdot)$ and
$N_{t}$ is the number of basis functions required to represent $J_{t}(\cdot)$.

_Proof:_ The proof follows by substituting for function $\phi(f(x))$ using
Assumption 1 to equation (6) and then solving for the coefficients
$\bar{\alpha}^{*}_{t}$. (Detailed proof: See Appendix).

###### Remark 1

The solution above is the discrete time analog of the classical Galerkin
procedure for solving Partial Differential Equations (PDEs) [22]. This should
not come as a surprise since the Policy Evaluation equation (6) is really a
discrete time analog of the continuous time PDE: $\frac{\partial J}{\partial
t}+\bar{c}+\bar{f}\frac{\partial J}{\partial x}=0$, solved with the terminal
condition $J_{T}(x)=g(x)$, and where $\dot{x}=\bar{f}(x)$ represents the
continuous dynamics and the continuous time cost is given by the integral
$\int_{0}^{T}\bar{c}(x)dt$.

Next, we show that the RL least squares solution (10) converges to the above
true solution in the mean square sense as the number of samples $R$ becomes
large.

###### Proposition 2

PE convergence. Let Assumption 1 hold. Further, let the number of basis
functions at time $t$ required be $N_{t}$. Given that all the required basis
functions at time $t$ are considered, the RL least square estimate (10)
converges to the true solution in the mean square sense.

_Proof:_ See Appendix.

The above result shows that the RL least squares procedure is a randomized
approximation of the PE equation. However, the above result follows under an
idealized situation when all necessary basis functions are considered, and $R$
becomes very large. Thus, in the following, we characterize the bias and the
variance of the estimate, which in turn will allow us to find the sample
complexity of the estimates.

###### Corollary 1

Bias in RL estimate (10). Suppose that the number of basis functions required
at time $t$ is $N_{t}$ and only $N<N_{t}$ are used. Then, the RL least squares
estimate (10) is biased for all $\tau<t$.

_Proof:_ See Appendix.

Next, we consider the variance of the estimate. The key role here is played by
the Gram matrix $\mathcal{G}_{t}=[<\phi_{i},\phi_{j}>_{t}]$, $i,j=1\cdots
N_{t}$, where $<.,.>_{t}$ denotes the inner product with respect to the
sampling distribution $p_{t}(.)$. In general, the variance of the solution is
determined by the variance in the R-sample empirical Gram matrix estimate
$\mathcal{G}^{R}_{t}=\frac{1}{R}\Phi_{t}^{R}\Phi_{t}^{R^{\prime}}$. In the
following, we characterize the variance of the empirical Gram matrix
$\mathcal{G}^{R}_{t}$ for suitable choice of basis functions and sampling
distribution.

###### Assumption 2

Let the basis functions used at time $t$ be
$\\{\phi^{1},\cdots\phi^{N_{t}}\\}$. We assume that there exists a constant
matrix $H_{t}$ such that $\begin{bmatrix}\phi^{1}\\\ \vdots\\\
\phi^{N_{t}}\end{bmatrix}=H_{t}\begin{bmatrix}1\\\ \vdots\\\
x^{M_{t}}\end{bmatrix},$ i.e., the basis functions at time $t$ can be
represented as $M_{t}$ degree polynomials.

Let us define

$H_{t}=[H_{1},\cdots H_{M_{t}}]=\begin{bmatrix}H^{1}\\\ \vdots\\\
H^{N_{t}}\end{bmatrix},$ (11)

i.e., we define the rows and columns of the matrix $H_{t}$. The covariance of
the error in the LS estimate is given by
$P_{t}^{R}=\sigma_{v}^{2}(G^{R}_{t})^{-1}$, where
$G^{R}_{t}=R\mathcal{G}^{R}_{t}$, and
$\mathcal{G}^{R}_{t}=\frac{1}{R}\Phi^{R}_{t}\Phi^{R^{\prime}}_{t}$ (see Proof
of Proposition 2). Now, we can characterize the size of the error in the LS
estimate as a function of number of samples $R$ required.

###### Theorem 1

Variance of RL least squares estimate. Let
$x_{t}\sim\mathcal{N}(0,\sigma_{X}^{2}),$ i.e. the sampling distribution is
zero mean Gaussian with variance $\sigma_{X}^{2}$. Let $\beta<1,\delta>0$, and
$n<\infty$ be given. Then, to probabilistically bound the norm of the error
covariance:

$Prob(||P_{t}^{R}||\leq\delta)>1-2e^{-n^{2}/2},$ (12)

the number of samples required are:

$R>\max[(\frac{n}{\beta}CC^{\prime})^{2}\sigma_{2M_{t}}^{2},\frac{\sigma_{v}^{2}C}{\delta(1-\beta)}],$
(13)

where

$\sigma_{2M_{t}}^{2}=[(4M_{t}-1)!!-(2M_{t}-1)!!^{2}]\sigma_{X}^{4M_{t}}$ (14)

and $C$ and $C^{\prime}$ are constants such that $||\mathcal{G}_{t}^{-1}||\leq
C$, and $||H_{M_{t}}||||H^{N_{t}}||\leq C^{\prime}$.

_Proof:_ See Appendix.

The above result establishes the number of samples required to get an accurate
RL least squares estimate. Next, we shall see the implications of the above
result for particular choices of basis functions.

###  III-A Sample Complexity

#### Monomial Basis

For a monomial basis, $C^{\prime}=1$ since $H$ is the identity matrix and
$N_{t}=M_{t}$. Thus, the number of samples has to satisfy:

$R\sim O([(4N_{t}-1)!!-(2N_{t}-1)!!^{2}]\sigma_{X}^{4N_{t}}),$ (15)

for the LS error to be small enough.

#### Hermite (Orthonormal) Basis

It is well known that the Hermite polynomials form the orthonormal basis for
Gaussian sampling distributions [22]. Noting that $M_{t}=N_{t}$, in the case
of Hermite polynomials, owing to their orthonormality, one can show that
$C^{\prime}=\frac{1}{N_{t}!^{2}\sigma_{X}^{4N_{t}}}$, and this, in turn,
implies that for the Hermite basis, the number of samples $R$ need to satisfy:

$R\sim O\left(\frac{(4N_{t}-1)!!-(2N_{t}-1)!!^{2}}{N_{t}!^{2}}\right),$ (16)

for the LS error to be small enough.

#### Nonlinear Basis

In this case, we mean that the cost function $J_{t}(x)=h(x,\theta_{t})$, where
$h(x,\theta)$ is a suitable nonlinear approximation architecture, such as a
(deep) neural net, parametrized nonlinearly by the (vector) parameter
$\theta$. The PE equation in this case becomes:

$h(x,\theta_{t})=c(x)+h(f(x),\theta_{t+1}),$ (17)

where $\theta_{t}$ parametrizes the cost function at time $t$, and the same
holds for $\theta_{t+1}$. It is reasonable to assume, given the change in the
parameter $\theta$ between consecutive steps is small enough, that:

$h(x,\theta_{t})\approx
h(x,\theta_{t+1})+\sum_{i=1}^{N_{t}}H_{t+1}^{i}(x)\delta\theta_{i},$ (18)

where $H_{t+1}^{i}(x)=\frac{\partial
h(x,\theta)}{\partial\theta_{i}}|_{\theta_{t+1}}$, where $\theta_{i}$ are the
components of the vector parameter $\theta$. Rewriting the policy evaluation
equation, one obtains:

$H_{t+1}^{1}(x)\delta\theta_{1}+\cdots+H_{t+1}^{N_{t}}(x)\delta\theta_{N_{t}}\\\
=c(x)+\underbrace{[h(f(x),\theta_{t+1}]-h(x,\theta_{t+1})]}_{\delta
h_{t+1}(x)}+v,$ (19)

where $v$ is a noise term, which may be written in matrix form for $R$ samples
as:

$[\delta\theta_{1}\cdots\delta\theta_{N_{t}}]\underbrace{\begin{bmatrix}H_{t+1}^{1}(x^{(1)}_{t})&\cdots&H_{t+1}^{1}(x^{(R)_{t}})\\\
\vdots&\vdots&\vdots\\\
H_{t+1}^{N_{t}}(x^{(1)}_{t})&\cdots&H_{t+1}^{N_{t}}(x^{(R)}_{t})\end{bmatrix}}_{\mathcal{H}^{R}_{t}}\\\
=\underbrace{\begin{bmatrix}c(x^{(1)})~{}\cdots~{}c(x^{(R)})\end{bmatrix}}_{\mathcal{C}^{R}_{t}}+\underbrace{\begin{bmatrix}\delta
h_{t+1}(x^{(1)})~{}\cdots~{}\delta
h_{t+1}(x^{(R)})\end{bmatrix}}_{\delta\mathcal{H}^{R}_{t}}\\\
+\underbrace{\begin{bmatrix}v^{(1)}_{t}&\cdots&v^{(R)}_{t}\end{bmatrix}}_{V^{R}_{t}},$

where we assume that the noise terms $v^{(i)}_{t}$ are i.i.d with variance
$\sigma_{v}^{2}$ as before. The least squares solution to the above equation
is, as usual, given by:

$\delta\theta^{R}_{t}=\mathcal{C}^{R}_{t}\mathcal{H}^{R^{\prime}}_{t}(\mathcal{H}^{R}_{t}\mathcal{H}^{R^{\prime}}_{t})^{-1}+\delta
H^{R}_{t}\mathcal{H}^{R^{\prime}}_{t}(\mathcal{H}^{R}_{t}\mathcal{H}^{R^{\prime}}_{t})^{-1}.$
(20)

Denote the empirical Gram matrix for the instantaneous basis functions,
$H_{t+1}^{i}(x)$, at time $t+1$ as:

$\mathcal{G}_{t}^{R}=\frac{1}{R}\mathcal{H}^{R}_{t}\mathcal{H}^{R^{\prime}},$
(21)

and thus, the covariance of the error in the LS solution is given by
$\frac{\sigma_{v}^{2}}{R}(\mathcal{G}_{t}^{R})^{-1}$. Note that this situation
is no different from the linear case considered previously, in that we are
approximating the change in the cost function via the change in the parameter
$\theta$ at time $t+1$, and the instantaneous basis functions
$H_{t+1}^{i}(x)$. The primary difference with the linear case is that our
basis functions $H_{t+1}^{i}(x)$ change with time unlike the fixed basis in
the linear case. Suppose that we have: $\begin{bmatrix}H_{t+1}^{1}(x)\\\
\vdots\\\ H_{t+1}^{N_{t}}(x)\end{bmatrix}=\mathcal{D}_{t+1}\begin{bmatrix}1\\\
x\\\ \vdots\\\ x^{M_{t}}\end{bmatrix},$ then we are back to the conditions
enumerated in Theorem 1. Therefore, it follows that the number of samples $R$
required such that we get a small enough error in the solution should be:

$R\sim O([(4M_{t}-1)!!-(2M_{t}-1)!!^{2}]\sigma_{X}^{4M_{t}}).$ (22)

Discussion. It can be seen clearly from above that choosing orthonormal (o.n.)
polynomials to the sampling distribution greatly reduces the variance of the
RL least squares estimate (10). Further, the variance of the estimate if we
use unnormalized polynomials is very high. The situation is no different even
when using a nonlinear basis. Typically in RL, one uses rollouts and the basis
functions are almost never chosen to satisfy orthonormality. Further, in
general, since the rollout sampling distributions $p_{t}(\cdot)$ need not be
Gaussian, finding such o.n. basis functions is a challenge in itself since one
does not have an idea about these sampling distributions in advance, or they
are never known explicitly. However, in our opinion, the sampling
distributions can, and should, be chosen at our convenience and need not arise
from rollouts if viewed in the context of the Galerkin interpretation of
Proposition 1. In fact, the above analysis suggests that if we are to reduce
the variance of the estimates, the sampling distributions should not be chosen
from rollouts.

A significantly more intractable problem arises due to the exponential growth
of the basis functions required at every time step. Suppose that the dynamics
could be represented by a $k$ degree polynomial, and the terminal cost
function was degree $N_{T}$, then the number of basis functions required at
time $t$ is $k^{T-t}N_{T}$. When coupled with the variance estimate above, one
can see that, even with o.n. basis, this leads to an explosive variance
growth: a factorial-exponential rate of growth.

##  IV A Perturbation Structure for PE

In this section, we present a construct that ideally allows a perturbation
structure to the Policy Evaluation equations, i.e., a structure where the
higher order terms do not affect the evolution of the lower order terms which
implies that we can close our computations at any desired order without
incurring an error.

Figure 1: Time-varying versus time-invariant perturbation expansion.

Recall the Policy evaluation equation (6). Now, consider that we are given a
nominal trajectory under the dynamics, say $\bar{x}_{t}=f(\bar{x}_{t-1})$,
$t=0,1,\cdots T$ given some initial condition $\bar{x}_{0}=x_{0}$. Next,
expand the dynamics about this nominal trajectory as:

$f(x_{t})=f(\bar{x}_{t}+\delta x_{t})=f(\bar{x}_{t})+F^{1}_{t}\delta
x_{t}+F^{2}_{t}\delta x_{t}^{2}+\cdots,$ (23)

where $F_{t}^{i}$ denotes the $i^{th}$ term in the expansion of the dynamics
around the nominal trajectory. Similarly:

$J_{t+1}(f(x_{t}))=J_{t+1}(f(\bar{x}_{t}+\delta x_{t}))\\\
=J_{t+1}(f(\bar{x}_{t}))+K_{t+1}^{1}(F_{t}^{1}\delta
x_{t}+\frac{1}{2}F_{t}^{2}\delta x_{t}^{2}+\cdots)\\\
+\frac{1}{2}K_{t+1}^{2}(F_{t}^{1}\delta x_{t}+F_{t}^{2}\delta
x_{t}^{2}+\cdots)^{2}+\cdots,$ (24)

where $K_{t+1}^{i}$ denotes the $i^{th}$ term in the Taylor expansion of
$J_{t+1}(f(x_{t}))$ around the nominal. Similarly the incremental cost
function $c_{t}(x_{t})=c_{t}(\bar{x}_{t})+C_{t}^{1}\delta
x_{t}+\frac{1}{2}C_{t}^{2}\delta x_{t}^{2}+\cdots$, and the optimal cost
function at time $t$, $J_{t}(x_{t})=J_{t}(\bar{x}_{t})+K_{t}^{1}\delta
x_{t}+\frac{1}{2}K_{t}^{2}\delta x_{t}^{2}+\cdots$. Substituting the above
expressions into the policy evaluation equation (6), we obtain:

$J_{t}(x_{t})=\bar{J}_{t}+K_{t}^{1}\delta x_{t}+\frac{1}{2}K_{t}^{2}\delta
x_{t}^{2}+\cdots\\\ =(\bar{c}_{t}+C_{t}^{1}\delta
x_{t}+\frac{1}{2}C_{t}^{2}\delta x_{t}^{2}+\cdots)+\bar{J}_{t+1}\\\
+K_{t+1}^{1}(F_{t}^{1}\delta x_{t}+\frac{1}{2}F_{t}^{2}\delta
x_{t}^{2}+\cdots)\\\ +\frac{1}{2}K_{t+1}^{2}(F_{t}^{1}\delta
x_{t}+F_{t}^{2}\delta x_{t}^{2}+\cdots)^{2}+\cdots,$ (25)

where $\bar{J}_{t}=J(\bar{x}_{t})$, $c(\bar{x}_{t})=\bar{c}_{t}$ and
$\bar{J}_{t+1}=J_{t+1}(\bar{x}_{t+1})=J_{t+1}(f(\bar{x}_{t}))$. Now, grouping
the different terms of $\delta x_{t}^{i}$ on both sides of the equation allows
for writing the following vector-matrix form equation as:

$[K_{t}^{1}\,K_{t}^{2}\cdots]=[C_{t}^{1}\,C_{t}^{2}\cdots]\\\
+[K_{t+1}^{1}\,K_{t+1}^{2},\cdots]\underbrace{\begin{bmatrix}F_{t}^{1}&F_{t}^{2}&\cdots\\\
0&(F_{t}^{1})^{2}&\cdots\\\ 0&0&\cdots\\\ \vdots&\vdots&\vdots\\\
0&0&\cdots\end{bmatrix}}_{\mathcal{B}_{t}},$ (26)

note the equations have a beautiful perturbation/ upper triangular structure.
Practically, this means that we can close our computations at any order we
desire without worrying about the effect of the higher order terms on the
lower order terms, given we have knowledge of the dynamics, and hence, the
Taylor coefficients $F_{t}^{i}$.

A special case of the above equation is when $f(0)=0$, and the nominal
trajectory is simply $\bar{x}_{t}=0$. In such a case the expansion is about a
nominal trajectory that stays at the origin. In such a case, the Taylor
coefficients, rather than being time varying, will be time invariant, i.e.,
$F^{1},F^{2},\cdots$ etc. Typically, we are given problems where the initial
state $x_{0}\neq 0$, and can be far from the origin. In such a case, we may
see that the number of terms required for a static expansion, i.e., about
$\bar{x}_{t}=0$, will require far more terms than would an expansion that was
centered on a nominal trajectory starting at $x_{0}$. The situation is
illustrated in Fig. 1. Thus, it is much more efficient to seek the time
varying expansion above. In particular, we shall explore the implication
further when we solve the perturbed policy evaluation (PPE) equation (26).

###### Remark 2

Note that given the Taylor coefficients $F_{t}^{i}$, the second and higher
rows of the “dynamics” matrix encoding the structure of the PPE equation are
perfectly known. This knowledge can be used to solve the PPE equation in a
“model based” fashion, as opposed to a model-free approach, where this
structure is actually teased out of the data from the system.

###  IV-A RL type solution to the PPE

First, we show a model-based solution to the PPE equation (26), i.e., one
where we explicitly estimate the first $M$ Taylor coefficients of the dynamics
$\mathcal{F}_{t}\equiv[F_{t}^{1},F_{t}^{2},\cdots F_{t}^{M}]$ which are then
substituted into (26) to solve the PE equation. Next, we show how this can be
extended to the model-free case: one where we do not solve for
$\mathcal{F}_{t}$, and instead directly solve (26) and infer the matrix
$\mathcal{B}_{t}$ from the system data.

We can write the following approximation, after neglecting the higher order
terms beyond $\delta x_{t}^{M}$:
$\delta x_{t+1}^{(i)}\approx\mathcal{F}_{t}\begin{bmatrix}\delta
x_{t}^{(i)}\\\ (\delta x_{t}^{(i)})^{2}\\\ \vdots\\\ (\delta
x_{t}^{(i)})^{M}\end{bmatrix}+v_{t}^{(i)}$, where as before $v_{t}^{(i)}$ is
an i.i.d. noise sequence and $i=1,2\cdots R$. A least squares estimate of
$\mathcal{F}_{t}$ is quite straightforward and may be written as:

$\mathcal{F}_{t}^{R}=\delta
X_{t+1}^{R}\delta\chi_{t}^{R^{\prime}}(\delta\chi_{t}^{R}\delta\chi_{t}^{R^{\prime}})^{-1},$
(27)

where $\delta\chi_{t}^{R}=\begin{bmatrix}\delta x_{t}^{(1)}&\cdots&\delta
x_{t}^{(R)}\\\ \vdots&\vdots&\vdots\\\ (\delta x_{t}^{(1)})^{M}&\cdots&(\delta
x_{t}^{(R)})^{M}\end{bmatrix}$, and $\delta X_{t+1}^{R}=\begin{bmatrix}\delta
x_{t+1}^{(1)}&\cdots&\delta x_{t+1}^{(R)}\end{bmatrix}$, where $\delta
x_{t}^{(i)}=x_{t}^{(i)}-\bar{x}_{t}$, and $\delta
x_{t+1}^{(i)}=f(x_{t}^{(i)})-\bar{x}_{t+1}$, where
$\bar{x}_{t+1}=f(\bar{x}_{t})$. Further, we assume that $\delta
x_{t}^{(i)}\sim\mathcal{N}(0,\sigma_{X}^{2})$. Thus, the data is obtained by
perturbing the system from the nominal trajectory.

The following development characterizes the error in the LS solution (27)
incurred from neglecting the higher order terms of the dynamics (beyond
$\delta x_{t}^{M}$) and shows that it can be made arbitrarily small by
choosing the perturbation $\delta x_{t}$ to be suitably small. In the
following, unlike in Section III, the Gram matrix $\mathcal{G}$ size will not
change, since we are looking at order $M$ approximation throughout time.

###### Lemma 1

Let
$\Delta_{t}^{R}=[\Delta_{t}^{R,1},\Delta_{t}^{R,2}\cdots\Delta_{t}^{R,M}]$,
where $\Delta_{t}^{R,l}=\frac{1}{R}\sum_{k>M}\sum_{i=1}^{R}F_{t}^{k}(\delta
x_{t}^{(i)})^{k}(\delta x_{t}^{(i)})^{l}$. Let the empirical Gram matrix
$\mathcal{G}^{R}=[\mathcal{G}^{R}_{ij}]$, where $i,j=1,2,\cdots M$, and
$\mathcal{G}^{R}_{ij}=\frac{1}{R}\sum_{k=1}^{R}(\delta x_{t}^{(k)})^{i}(\delta
x_{t}^{(k)})^{j}$. Then,

$\mathcal{F}_{t}^{R}=\mathcal{F}_{t}+\Delta_{t}^{R}{(\mathcal{G}^{R})^{-1}}+\frac{1}{R}V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}.$
(28)

_Proof:_ See Appendix.

The above result makes it clear that our estimates are biased for any finite
$R$ (in fact, even for the limit), but if $\Delta^{R}_{t}$ is small enough,
then this bias can be made small. In the following, we show precisely such a
result.

###### Proposition 3

Let $\Delta_{t}=\lim_{R}\Delta_{t}^{R}$ and let
$\mathcal{G}=\lim_{R}\mathcal{G}^{R}$. Given any $M$, any $\epsilon>0$, there
exists a variance $\sigma_{X}^{2}<\infty$ such that
$|\Delta_{t}^{j}|\leq\epsilon|\tilde{f}_{t}^{j}|$, where
$\mathcal{F}_{t}\mathcal{G}=[\tilde{f}_{t}^{1},\cdots,\tilde{f}_{t}^{j},\cdots\tilde{f}_{t}^{M}]$.

_Proof:_ See Appendix.

Next, we have the following consequence.

###### Corollary 2

The least squares estimate (27),
$\mathcal{F}_{t}^{R}\rightarrow\mathcal{F}_{t}+\Delta_{t}\mathcal{G}^{-1}$ as
$R\rightarrow\infty$ in mean square sense.

_Proof:_ See Appendix.

Accuracy of the Solution. To understand the result above, let us rewrite the
limiting solution as
$\mathcal{F}_{t}^{R}=(\mathcal{F}_{t}\mathcal{G}+\Delta_{t})\mathcal{G}^{-1}$.
Thus, the limiting solution can be understood as the solution to the linear
equation $\mathcal{F}_{t}\mathcal{G}=\tilde{\mathcal{F}}_{t}+\Delta_{t}$ where
$\tilde{\mathcal{F}}_{t}=\mathcal{F}_{t}\mathcal{G}$. We have shown in
Proposition 3 that for small enough variance,
$||\Delta_{t}||\leq\epsilon||\tilde{\mathcal{F}}_{t}||$, and thus we can
expect a small error on the right hand side of the equation above. However,
albeit the error $\Delta_{t}$ may be small compared to the “signal
$\tilde{\mathcal{F}}_{t}$”, the actual error in the solution is affected by
the conditioning of the Gram matrix $\mathcal{G}$. In fact, one can show that:

$\displaystyle{\frac{||\mathcal{F}_{t}^{\infty}-\mathcal{F}_{t}||}{||\mathcal{F}_{t}||}\leq\kappa(\mathcal{G})\epsilon},$
(29)

where $\mathcal{F}_{t}^{\infty}=\lim_{R}\mathcal{F}_{t}^{R}$, and
$\kappa(\mathcal{G})$ denotes the condition number of the Gram matrix
$\mathcal{G}$. Thus, the true pacing item in the accuracy of the solution is
the conditioning of the matrix $\mathcal{G}$. In fact, if the input is
Gaussian, then the conditioning of the matrix rapidly deteriorates as $M$
increases since the higher moments of a Gaussian increase as
$(2M-1)!!\sigma_{X}^{2M}$, i.e, in factorial-exponential fashion. Thus, in
practice, one cannot make $M$ large as the solution becomes highly sensitive
due to the ill conditioning of the Gram matrix.

Variance of Solution. As shown previously, the variance of the solution is
directly proportional to the variance in the empirical Gram matrix
$\mathcal{G}^{R}$, and thus, the number of samples required is
$R\sim\mathcal{O}([(4M-1)!!-(2M-1)!!^{2}]\sigma_{X}^{4M})$ as in (15). Hence,
the central role in the accuracy and in the variance of the least squares
estimate is played by the empirical Gram matrix $\mathcal{G}^{R}$.

Model-based or Model-free? The number of computations in the model-free method
is $\approx 2\times$ the computation in the model-based method, since $R$ is
typically large. Thus, the model-free approach will result in higher
computational efforts or $\approx 2\times$ the variance as compared to model-
based methods. _(Detailed Discussion: See Appendix. )_

Orthonormal and Nonlinear Basis Functions: As we saw previously, the variance
of the RL-LS solution (10) is highly significantly reduced when using o.n.
basis functions such as the Hermite polynomials, see (15) vs (16). Thus, it is
of interest to see if any added advantage can be gained by using such o.n.
functions in the PPE case.
We may write the Hermite polynomials in terms of the monomials as:
$\begin{bmatrix}h_{1}(\delta x)\\\ \vdots\\\ h_{N}(\delta
x)\end{bmatrix}=H\begin{bmatrix}\delta x\\\ \vdots\\\ \delta
x^{N}\end{bmatrix}$, where $h_{i}(\cdot)$ represent the Hermite polynomials,
and $H$ is a lower triangular matrix that encodes the linear transformation
from the monomials to the Hermite polynomials. Similarly, let $H^{-1}$
represent the inverse transformation from the Hermite polynomials to the
monomials. Then, the PPE equation 26 may be written in the Hermite basis as:

$\mathcal{K}_{t}^{H}=\mathcal{C}^{H}_{t}+\mathcal{K}_{t+1}^{H}\mathcal{B}_{t}^{H},$
(30)

where $\mathcal{K}_{t}^{H}$, and $\mathcal{C}_{t}^{H}$ consist of the Hermite
coefficients (rather than the Taylor coefficients) and
$\mathcal{B}_{t}^{H}=H\mathcal{F}_{t}H^{-1}$. Note that $H$ and $H^{-1}$ are
lower triangular, while $\mathcal{F}_{t}$ is upper triangular, and therefore,
$\mathcal{B}_{t}^{H}$ is fully populated, i.e., the PPE equations in the
Hermite basis lose their perturbation/ upper triangular structure. This is the
primary shortcoming of the o.n. representation, in particular, we shall show
the necessity of the perturbation structure to an accurate solution in the
following section after we outline the stochastic case.

Further, in the nonlinear basis case, nothing changes from the LS problem in
(20), since the formulation cannot distinguish deviations from a trajectory.
Thus, the equations for the parameters, in general, will be fully coupled, and
not have a perturbation structure.

##  V The Stochastic Case

Now, we consider the stochastic case. The policy evaluation equation in the
stochastic case becomes:

$J_{t}(x_{t})=c_{t}(x_{t})+E[J_{t+1}(x_{t+1})],$ (31)

where $x_{t+1}=f(x_{t})+\omega_{t}$, and $\omega_{t}$ is a Gaussian white
noise sequence.
Let us define the nominal trajectory as the noise free case, i.e.,
$\bar{x}_{t+1}=f(\bar{x}_{t})$, given some nominal initial condition
$\bar{x_{0}}$. Next, let us consider perturbations $\delta x_{t}$ about this
nominal. Thus, we get

$J_{t+1}(x_{t+1})=J_{t+1}(f(x_{t})+\omega_{t})=J_{t+1}(f(\bar{x}_{t}+\delta
x_{t})+\omega_{t})$. Using:

$f(\bar{x}_{t}+\delta x_{t})=f(\bar{x}_{t})+F_{t}^{1}\delta
x_{t}+F_{t}^{2}\delta x_{t}^{2}+\cdots$,

$J_{t}(x_{t})=J_{t}(\bar{x}_{t})+K_{t}^{1}\delta x_{t}+K_{t}^{2}\delta
x_{t}^{2}+\cdots$,

and thus,

$J_{t+1}(x_{t+1})=J_{t+1}(\bar{x}_{t+1})+K_{t+1}^{1}(\delta
f_{t}+\omega_{t})+K_{t+1}^{2}(\delta f_{t}+\omega_{t})^{2}+\cdots$,

where $\delta f_{t}=F_{t}^{1}\delta x_{t}+F_{t}^{2}\delta x_{t}^{2}+\cdots$.
Thus, using the fact that the noise sequence is white and Gaussian:

$E[J_{t+1}(x_{t+1})]=(\bar{J}_{t+1}+K_{t+1}^{2}+\cdots)+\\\
(K_{t+1}^{1}F_{t}^{1}+K_{t+1}^{3}(3F_{t}^{1})+\cdots)\delta x_{t}+\cdots.$
(32)

In particular, we may write:

$E[J_{t+1}(x_{t+1})]=\sum_{k}G^{k}(\mathcal{F}_{t},\mathcal{K}_{t+1})\delta
x_{t}^{k},$ (33)

where recall that $\mathcal{F}_{t}$ and $\mathcal{K}_{t+1}$ are the Taylor
coefficients of the dynamics at time $t$ and the cost function at time $t+1$
respectively, and $G^{k}(\mathcal{F}_{t},\mathcal{K}_{t+1})$ are suitably
defined functions of the Taylor coefficients
$\mathcal{F}_{t},\mathcal{K}_{t+1}$. Therefore, the policy evaluation
equations in the stochastic case become (equating the coefficients of the
different powers of $\delta x_{t}$ on both sides of the equation):
$\bar{J}_{t}=\bar{c}_{t}+G^{0}(\mathcal{F}_{t},\mathcal{K}_{t+1}),~{}K_{t}^{1}=C_{t}^{1}+G^{1}(\mathcal{F}_{t},\mathcal{K}_{t+1}),~{}K_{t}^{2}=C_{t}^{2}+G^{2}(\mathcal{F}_{t},\mathcal{K}_{t+1}),~{}\cdots.$
The import of the above equations is that, unlike in the deterministic case,
there is no perturbation structure in the stochastic case, and as a result all
terms affect all other terms. This is the fundamental issue with the
stochastic case, and makes the computation of a solution, in general,
intractable.

###  V-A The Necessity of the Perturbation Structure.

Now, let us consider why the perturbation structure is critical to a solution.
Write the $G^{k}(\cdot)$ above as
$G^{k}(\mathcal{F}_{t},\mathcal{K}_{t+1})=G^{k,M}(\mathcal{F}_{t}^{M},\mathcal{K}_{t+1}^{M})+\delta
G^{k,M}$, i.e., we keep only the first $M$ terms in the Taylor series
$\mathcal{F}_{t}$ and $\mathcal{K}_{t+1}$, denoted by $\mathcal{F}_{t}^{M}$
and $\mathcal{K}_{t+1}^{M}$ respectively, and $\delta G^{k,M}$ denotes the
error resulting from neglecting the higher order terms. Then, for an accurate
solution, we require that $M$ is large enough such that $||\delta
G^{k,M}||\leq\epsilon||G^{k,M}||$, for suitably small $\epsilon$. In the case
of the PPE eq. (26), we can choose any $M$, no matter how small, and still be
assured of zero error in computing those coefficients, given the model
knowledge. However, even with model knowledge, the computations are
intractable in the stochastic case since the number of terms $M$ can be
expected to be large for most problems (except linear problems).

Next, consider the RL scenario: the stochastic case is essentially similar to
the deterministic case in that the LS solution is still given by a slight
modification of eq. (10) where now:
$\Phi^{R}_{t+1}=\begin{bmatrix}\phi^{1}(x_{t+1}^{(1)})&\cdots&\phi^{1}(x_{t+1}^{(R)})\\\
\vdots&\vdots&\vdots\\\
\phi^{N}(x_{t+1}^{(1)})&\cdots&\phi^{N}(x_{t+1}^{(R)}),\end{bmatrix}$ and the
random next states $x_{t+1}^{(i)}$ are sampled from the stochastic dynamics:
$x_{t+1}^{(i)}=f(x_{t}^{(i)})+\omega_{t}^{(i)}$. Thus, there will be a higher
variance in the elements of the $\Phi_{t+1}^{R}$ matrix owing to the input
noise $\omega_{t}$. In this case, due to the additional variance,
$\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}$ will take longer to converge than
$\Phi_{t}^{R}\Phi_{t}^{R^{\prime}}$. Nonetheless, the error covariance
estimates do not change from the deterministic case, and thus, the sample
complexity depends on the empirical Gram matrix $\mathcal{G}^{R}$, and remains
the same (at least in the $\mathcal{O}(\cdot)$ sense).

Bias and Variance. It is inevitable, even with the perturbation structure of
the deterministic case, that the higher order terms corrupt the computations
of the lower order terms in policy evaluation when done via RL. However, as
shown in Proposition 3, this error can be made arbitrarily small by
controlling the input variance, and we can close our computations at any
suitable $M$ such that the Gram matrix is well conditioned. However, in the
stochastic case, we do not have this luxury since for a good solution, we will
have to choose a suitably large $M$. This, in turn, implies that the bias is
adversely affected due to the ill-conditioning of the Gram matrix,
$\mathcal{G}$, for large $M$. Furthermore, owing to the same reason, the
variance of the empirical Gram matrix $\mathcal{G}^{R}$ is also bound to be
very high. Thus, in the stochastic case, we see that the RL solution will be
overwhelmed by both the bias as well as the variance, thereby making the
solution highly unreliable, i.e., one can always obtain a solution to the
stochastic case by solving a least squares problem, but it is bound to be
inaccurate owing to the fundamental structure of the problem. Finally, the
situations does not change if we use an o.n. basis like the Hermite basis or a
nonlinear basis, due to the lack of a perturbation structure.

##  VI Empirical Results

Thus far in this paper, we have established theoretical results that show that
finding a global (higher order) solution for the policy evaluation problem is
subject to very high error and variance. In this section, we provide empirical
evidence of this explosive growth of errors. In order to accomplish this, we
need a system for which we know the true solution to the PE problem. To this
end, we assume a discretized system of the form:

$\displaystyle x_{t+1}=x_{t}+\delta(-x_{t}+\epsilon x_{t}^{3}),$ (34)

where $\delta$ is the time discretization. Let us define the terminal cost to
be of the quadratic form: $J_{T}(x_{T})=\alpha x_{T}^{2}$ and the incremental
cost to be: $c_{t}(x_{t})=cx_{t}^{2}$, so that the DP equation can be written
as:

$\displaystyle J_{t}(x_{t})=cx_{t}^{2}+J_{t+1}(x_{t+1}).$ (35)

Since this system has a polynomial nonlinearity, it is straightforward to find
the cost functions backward in time as polynomials. Given that we have this
true answer, we can now find the error in our solution as we sweep back in
time, as a function of the number of samples, as well as the exploration
parameter $\sigma_{X}$ which determines whether we explore locally or
globally. We show results for polynomial approximations of order
$N=M=\\{6,12,18\\}$ as we back propagate 3 steps in time. The required number
of basis functions at time $t=\\{T,T-1,T-2,T-3\\}$ are
$M_{t}=N_{t}=\\{2,6,18,54\\}$, respectively. The parameter values assumed are:
$\epsilon=1,\delta=0.1,c=10,$ and $\alpha=10$.

Figure 2: Results for varying number of samples for exploration parameter
$\sigma_{X}=0.1$ and different orders of approximation, $M=6,12,18$, for the
backward sweep in time till $T-3$.

For a “small” exploration noise of $\sigma_{X}=0.1$, as we increase the number
of basis functions $M_{t}$, better results in the mean error are obtained as
we propagate back in time due to the increase in number of required basis
functions (1st column of Fig. 2). Although the mean error improves till M= 12,
there is a marked increase in the variance of the error with increased basis
function (2nd column of Fig. 2) from 12 to 18. Also, notice the large number
of samples required to reduce the variance of the error even for a one-
dimensional problem.

Figure 3: Results for varying number of samples for $\sigma_{X}=1$.

For the large exploration case with $\sigma_{X}=1$, the mean error and its
variance become very high (Fig. 3). Also, notice that even with a large number
of samples, the error values do not decrease as in the small noise case of
$\sigma_{X}=0.1$. Another important observation here is that albeit the
performance at $T-1$ is very good, as we sweep back in time, the errors and
their variance show explosive growth and become unacceptable.
Thus, albeit preliminary, this empirical evidence suggests that solution
accuracy, measured via the mean error and its variance, is indeed adversely
affected by higher order approximations and/ or large exploration. Conversely,
the only way to ensure accuracy is to have a suitably low order approximation
which needs to be enforced by a suitably local exploration.

##  VII Conclusions and Implications

In this paper, we have studied the inherent structure of the Reinforcement
Learning problem. Concentrating on the policy evaluation problem, we have
shown that unless we seek local solutions, the answers found are bound to
suffer from high variance, and thus, be inaccurate. In particular, the
deterministic problem has a perturbation structure that can be exploited to
obtain arbitrary accurate local solutions. It is also shown that the
stochastic problem lacks the perturbation structure, and thus, is bound to
suffer from high variance when solved in an RL fashion.
The primary issue one has to worry about now: “what now for stochastic
control?” It is intractable so the best seems to be the deterministic
approximation. However, the deterministic solution is optimal locally, and
thus, when allied with replanning of the nominal trajectory, we can recover at
least near-optimal solution for the stochastic case. We also conjecture that
this local replanning based approach is “fundamentally” the best that one can
hope to achieve via computation/ RL. Our future research will concentrate on
doing an extensive version of the experiments that we started in this paper,
and provide further evidence that local RL methods are the ones to pursue.

## References

  * [1] D. P. Bertsekas, _Dynamic Programming and Optimal Control, vols I and II_. Cambridge, MA: Athena Scientific, 2012\.
  * [2] R. E. Bellman, _Dynamic Programming_. Princeton, NJ: Princeton University Press, 1957.
  * [3] M. Lagoudakis and R. Parr, “Least squares policy iteration,” _Journal of Machine Learning Research_ , vol. 4, pp. 1107–1149, 2003.
  * [4] R. Akrour, A. Abdolmaleki, H. Abdulsamad, and G. Neumann, “Model free trajectory optimization for reinforcement learning,” in _Proc. of the ICML_ , 2016.
  * [5] E. Todorov and Y. Tassa, “Iterative local dynamic programming,” in _Proc. of the IEEE Int. Symposium on ADP and RL._ , 2009.
  * [6] E. Theodorou, Y. Tassa, and E. Todorov, “Stochastic differential dynamic programming,” in _Proc. of the ACC_ , 2010.
  * [7] S. Levine and P. Abbeel, “Learning neural network policies with guided search under unknown dynamics,” in _Advances in NIPS_ , 2014.
  * [8] S. Levine and K. Vladlen, “Learning complex neural network policies with trajectory optimization,” in _Proc. of the ICML_ , 2014.
  * [9] M. G. Azar, R. Munos, and B. Kappen, “On the sample complexity of reinforcement learning with a generative model,” _arXiv preprint arXiv:1206.6461_ , 2012.
  * [10] S. M. Kakade, “On the sample complexity of reinforcement learning,” Ph.D. dissertation, UCL (University College London), 2003.
  * [11] R. Munos and C. Szepesvári, “Finite-time bounds for fitted value iteration.” _Journal of Machine Learning Research_ , vol. 9, no. 5, 2008\.
  * [12] S. Dean, H. Mania, N. Matni, B. Recht, and S. Tu, “On the sample complexity of the linear quadratic regulator,” _Foundations of Computational Mathematics_ , pp. 1–47, 2019.
  * [13] B. Recht, “A tour of reinforcement learning: The view from continuous control,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 2, pp. 253–279, 2019.
  * [14] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, “Deep reinforcement learning that matters,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [15] M. N. G. Mohamed, S. Chakravorty, and R. Wang, “Optimality and tractability in stochastic nonlinear control,” _arXiv preprint arXiv:2004.01041_ , 2020.
  * [16] D. Q. Mayne, “Model predictive control: Recent developments and future promise,” _Automatica_ , vol. 50, pp. 2967–2986, 2014.
  * [17] J. B. Rawlings and D. Q. Mayne, _Model Predictive Control: Theory and Design_. Madison, WI: Nob Hill, 2015.
  * [18] L. Chisci, J. A. Rossiter, and G. Zappa, “Systems with persistent disturbances: predictive control with restricted contraints,” _Automatica_ , vol. 37, pp. 1019–1028, 2001.
  * [19] J. A. Rossiter, B. Kouvaritakis, and M. J. Rice, “A numerically stable state space approach to stable predictive control strategies,” _Automatica_ , vol. 34, pp. 65–73, 1998.
  * [20] D. Q. Mayne, E. C. Kerrigan, E. J. van Wyk, and P. Falugi, “Tube based robust nonlinear model predictive control,” _International journal of robust and nonlinear control_ , vol. 21, pp. 1341–1353, 2011.
  * [21] D. Mayne, “Robust and stochastic mpc: Are we going in the right direction?” _IFAC-PapersOnLine_ , vol. 48, no. 23, pp. 1 – 8, 2015, 5th IFAC Conference on NMPC 2015.
  * [22] R. Courant and D. Hilbert, _Methods of Mathematical Physics, vol. II_. New York: Interscience publishers, 1953, vol. 336.
  * [23] B. Gnedenko, _Theory of Probability_. New York, NY: Chelsea, 1968.
  * [24] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.

## APPENDIX : Detailed Proofs and Supplementary Material

The supplementary materials contain detailed proofs of the results that are
missing in the main paper.

###  VII-A Proof of Proposition 1

###### Proof:

Due to Assumption 1, we may write:

$\sum_{i=1}^{N_{t}}\alpha_{t}^{i}\phi^{i}(x)=\sum_{i=1}^{N}c_{t}^{i}\phi^{i}(x)+\sum_{i=1}^{N_{t+1}}\alpha_{t+1}^{i}\sum_{j=1}^{N_{t}}\beta^{ij}\phi^{j}(x).$
(36)

Taking inner products on both sides of the above equation, we obtain:

$\sum_{i=1}^{N_{t}}\alpha_{t}^{i}<\phi^{i},\phi^{k}>=\sum_{i=1}^{N_{t}}c_{t}^{i}<\phi^{i},\phi^{k}>+\sum_{i=1}^{N_{t+1}}\alpha_{t+1}^{i}\sum_{j}\beta^{ij}<\phi^{j},\phi^{k}>,\,\forall\,k,$
(37)

where $<\phi^{i},\phi^{k}>=\int\phi^{i}(x)\phi^{k}(x)p_{t}(x)dx$, i.e., the
weighted inner product using the density $p_{t}(\cdot)$. Denote the Gram
matrix $\mathcal{G}_{t}=[<\phi^{i},\phi^{k}>]$, $i,k=1\cdots N_{t}$, and
$\alpha_{t+1}^{i}=0$, for all $i>N_{t+1}$. The above equations can then be
written as:

$\bar{\alpha}_{t}\mathcal{G}_{t}=\bar{c}_{t}\mathcal{G}_{t}+\bar{\alpha}_{t+1}B_{t}\mathcal{G}_{t}.$
(38)

Taking the inverse of $\mathcal{G}_{t}$ on both sides, the answer follows. ∎

###  VII-B Proof of Proposition 2

###### Proof:

Owing to Assumption 1, $\Phi_{t+1}^{R}=B_{t}\Phi_{t}^{R}$, and thus,
$\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}=B_{t}\Phi_{t}^{R}\Phi_{t}^{R^{\prime}}.$
Let $\Phi_{t}^{R}\Phi_{t}^{R^{\prime}}\equiv G^{R}_{t}$. Then, the RL estimate
can be written as:

$\displaystyle\bar{\alpha}_{t}^{R}=\bar{c}_{t}G_{t}^{R}(G_{t}^{R})^{-1}+\bar{\alpha}_{t+1}^{R}B_{t}G_{t}^{R}(G_{t}^{R})^{-1}+V_{t}^{R}\Phi_{t}^{R^{\prime}}(G_{t}^{R})^{-1},$
(39)

where $V_{t}^{R}=[v_{t}^{(1)},\cdots v_{t}^{(R)}]$ is the noise sequence.
Noting that the noise sequence is i.i.d., the covariance of $V^{R}_{t}$ is
some $\sigma_{v}^{2}I_{R}$, where $I_{R}$ is the $R$ dimensional identity
matrix and $\sigma_{v}^{2}$ is the variance of the i.i.d. noise sequence.
Thus, the covariance of the least squares estimate above is given by

$\displaystyle P_{t}^{R}$
$\displaystyle=\sigma_{v}^{2}(G^{R}_{t})^{-1},~{}G^{R}_{t}=R\mathcal{G}^{R}_{t},$
$\displaystyle\mathcal{G}^{R}_{t}$
$\displaystyle\equiv\begin{bmatrix}\frac{1}{R}\sum_{i}\phi^{1}(x_{t}^{(i)})\phi^{1}(x_{t}^{(i)}),\cdots,\frac{1}{R}\sum_{i}\phi^{1}(x_{t}^{(i)})\phi^{N}(x_{t}^{(i)})\\\
\ddots\\\
\frac{1}{R}\sum_{i}\phi^{N}(x_{t}^{(i)})\phi^{1}(x_{t}^{(i)}),\cdots,\frac{1}{R}\sum_{i}\phi^{N}(x_{t}^{(i)})\phi^{N}(x_{t}^{(i)})\end{bmatrix}.$

We need to show that $\bar{\alpha}_{t}^{R}\rightarrow\bar{\alpha}_{t}$ in m.s.
sense. We proceed by induction. Let us assume that
$\bar{\alpha}_{t+1}^{R}\rightarrow\bar{\alpha}_{t+1}$ in m.s.s. Define
$S=\frac{1}{R}\bar{\alpha}_{t+1}^{R}\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}(\mathcal{G}_{t}^{R})^{-1}-\bar{\alpha}_{t+1}B_{t}$.
Then:

$E||\bar{\alpha}_{t}^{R}-\bar{\alpha}_{t}||^{2}=E||S||^{2}+E[S(V_{t}^{R}\Phi_{t}^{R^{\prime}}(\Phi_{t}^{R}\Phi_{t}^{R^{\prime}})^{-1})^{\prime}]+\text{Tr}\frac{1}{R^{2}}E[(\mathcal{G}_{t}^{R})^{-1}\Phi_{t}^{R}V_{t}^{R^{\prime}}V_{t}^{R}\Phi_{t}^{R^{\prime}}(\mathcal{G}_{t}^{R})^{-1}],$
(40)

where Tr denotes the trace of a matrix. Since
$\bar{\alpha}_{t+1}^{R}\rightarrow\bar{\alpha}_{t+1}^{R}$ in m.s.s by
assumption, and
$\frac{1}{R}\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}=B_{t}\mathcal{G}_{t}^{R}$,
$E||S||^{2}\rightarrow 0$ as $R\rightarrow\infty$. The second term is zero
since $V_{t}^{R}$ is independent of $\Phi_{t}^{R}$ and zero mean. For the last
term, note that:

$\displaystyle\text{Tr}\frac{1}{R^{2}}E[(\mathcal{G}_{t}^{R})^{-1}\Phi_{t}^{R}V_{t}^{R^{\prime}}V_{t}^{R}\Phi_{t}^{R^{\prime}}(\mathcal{G}_{t}^{R})^{-1}]=\frac{1}{R^{2}}E_{\mathcal{G}_{t}^{R}}[(\mathcal{G}_{t}^{R})^{-1}\Phi_{t}^{R}E[V_{t}^{R^{\prime}}V_{t}^{R}]\Phi_{t}^{R^{\prime}}(\mathcal{G}_{t}^{R})^{-1}|\mathcal{G}_{t}^{R}]=\frac{\sigma_{v}^{2}}{R}E_{\mathcal{G}_{t}^{R}}[(\mathcal{G}_{t}^{R})^{-1}]$
(41)

Noting that $\mathcal{G}_{t}^{R}\rightarrow\mathcal{G}_{t}$ almost surely, and
$\mathcal{G}_{t}$ is a finite matrix, it follows that
$\frac{\sigma_{v}^{2}}{R}E_{\mathcal{G}_{t}^{R}}[(\mathcal{G}_{t}^{R})^{-1}]\rightarrow
0$ as $R\rightarrow\infty$. Therefore, the above means that
$\bar{\alpha}_{t}^{R}\rightarrow\bar{\alpha}_{t}$ in m.s.s. if
$\bar{\alpha}_{t+1}^{R}\rightarrow\bar{\alpha}_{t+1}$ in m.s.s. A similar
argument as above can be used to show that
$\bar{\alpha}_{T}^{R}\rightarrow\bar{\alpha}_{T}$ in m.s.s. for the final time
time $T$, and thus, by backwards induction the result follows for all $t$.
∎

###  VII-C Proof of Corollary 1

###### Proof:

In general, since $N<N_{t}$,
$\phi^{i}(f(x))=\sum_{j=1}^{N_{t}}\beta^{ij}\phi^{j}(x)=\sum_{j=1}^{N}\beta^{ij}\phi^{j}(x)+\sum_{j=N+1}^{N_{t}}\beta^{ij}\phi^{j}(x)$.
Thus, $\Phi_{t+1}^{R}=B_{t}\Phi_{t}^{R}+\tilde{B}_{t}\tilde{\Phi}_{t}^{R}$,
where $B_{t}=[\beta^{ij}],i,j=1,\cdots N,$
$\tilde{B}_{t}=[\beta^{ij}],i=1,\cdots N;j=N+1,\cdots N_{t}$,

$\displaystyle\Phi_{t}^{R}=\begin{bmatrix}\phi^{1}(x_{t}^{(1)}),\cdots,\phi^{1}(x_{t}^{(R)})\\\
\ddots\\\
\phi^{N}(x_{t}^{(1)}),\cdots,\phi^{N}(x_{t}^{(R)})\end{bmatrix},~{}\text{and},~{}\tilde{\Phi}_{t}^{R}=\begin{bmatrix}\phi^{N+1}(x_{t}^{(1)}),\cdots,\phi^{N+1}(x_{t}^{(R)})\\\
\ddots\\\
\phi^{N_{t}}(x_{t}^{(1)}),\cdots,\phi^{N_{t}}(x_{t}^{(R)})\end{bmatrix}.$ (42)

Therefore,
$\Phi_{t+1}^{R}\Phi_{t}^{R^{\prime}}=B_{t}\Phi_{t}^{R}\Phi_{t}^{R^{\prime}}+\tilde{B}_{t}\tilde{\Phi}_{t}^{R}\Phi_{t}^{R^{\prime}}$,
and the second term
$\Delta\equiv\tilde{B}_{t}\tilde{\Phi}_{t}^{R}\Phi_{t}^{R^{\prime}}$ on the
right, in general, biases the coefficients $\alpha_{t}^{i}$, for $i=1\cdots
N$. Further, even if $\Delta=0$, which is feasible if the basis functions are
orthonormal, i.e., $<\phi^{i},\phi^{j}>=\delta_{ij}$, where $\delta_{ij}$ is
the Kronecker delta function, the coefficients $\alpha_{t}^{i}$, $i>N$, are
missed, which then affects the calculation of $J_{t-1}(\cdot)$, in particular,
the coefficients $\alpha_{t-1}^{i}$, $i=1,\cdots N$. Thus, the solution, even
in such a case, gets biased at the next time step $t-1$.
∎

###  VII-D Proof of Theorem 1

###### Proof:

Recall that
$\Phi^{R}_{t}=\begin{bmatrix}\phi^{1}(x^{(1)}_{t})\cdots\phi^{1}(x^{(R)}_{t})\\\
\ddots\\\
\phi^{N_{t}}(x^{(1)}_{t})\cdots\phi^{N_{t}}(x^{(R)}_{t})\end{bmatrix}.$

We have $\mathcal{G}^{R}_{t}=\mathcal{G}_{t}+\mathcal{E}^{R}_{t}$, where

$\displaystyle\mathcal{G}_{t}=\begin{bmatrix}E[\phi^{1}]^{2}\cdots
E[\phi^{1}\phi^{N_{t}}]\\\ \ddots\\\ E[\phi^{N_{t}}\phi^{1}]\cdots
E[\phi^{N_{t}}]^{2}\end{bmatrix}=\lim_{R}\frac{1}{R}\Phi^{R}_{t}\Phi_{t}^{R^{\prime}},$
(43)

is the Gram matrix for the basis functions $\phi^{i}$ at time $t$. Assuming
that $R$ is large enough, we may write:

$\displaystyle(G^{R}_{t})^{-1}=\frac{1}{R}(\mathcal{G}_{t}+\mathcal{E}^{R}_{t})^{-1}=\frac{1}{R}\sum_{k=0}^{\infty}(-1)^{k}\mathcal{G}_{t}^{-1}(\mathcal{E}^{R}_{t}\mathcal{G}_{t}^{-1})^{k}.$
(44)

Thus, taking (spectral) norms on both sides, we obtain:
$||(G^{R}_{t})^{-1}||\leq\frac{1}{R}\sum_{k}||\mathcal{G}_{t}^{-1}||^{k+1}||\mathcal{E}^{R}_{t}||^{k}.$
Using the assumption that $||\mathcal{G}_{t}^{-1}||<C$, we obtain:

$||(G^{R}_{t})^{-1}||\leq\frac{C}{R}\sum_{k=0}^{\infty}C^{k}||\mathcal{E}^{R}_{t}||^{k}.$
(45)

Next, we find a bound on $||\mathcal{E}^{R}_{t}||$.
Note that $\mathcal{E}^{R}_{t}=H_{t}\bar{\mathcal{E}}^{R}_{t}H^{\prime}_{t},$
where
$\bar{\mathcal{E}}^{R}_{t}=\begin{bmatrix}\bar{\epsilon}_{2}^{R}\cdots\bar{\epsilon}_{M_{t}+1}^{R}\\\
\ddots\\\
\bar{\epsilon}_{M_{t}+1}^{R}\cdots\bar{\epsilon}_{2M_{t}}^{R},\end{bmatrix}$
where
$\bar{\epsilon}_{p}^{R}=E[x^{p}]-\frac{1}{R}\sum_{i=1}^{R}(x^{(i)})^{p},$ i.e,
the error between the moment $E[x^{p}]$ and its R-sample empirical estimate.
Then, it is clear that $E[\bar{\epsilon}_{p}^{R}]=0$, and
$\text{var}[\bar{\epsilon}_{p}^{R}]=\frac{\text{var}(x^{p})}{R}.$ Also note
that $\text{var}[x^{p}]=E[x^{2p}]-(E[x^{p}])^{2}=(2p-1)!!\sigma_{X}^{2p}$, if
$p$ is odd, and $[(2p-1)!!-(p-1)!!^{2}]\sigma_{X}^{2p}$, if $p$ is even.
If $R$ is large enough, then owing to the Central Limit Theorem [23],
$\bar{\epsilon}_{p}^{R}\sim\mathcal{N}(0,\frac{var[x^{p}]}{R})$. Further,
given that $R$ is large enough, it is reasonable to assume that all
$\bar{\epsilon}_{p}^{R}$, $p<2M_{t}$, have converged, and thus, the error
matrix $\bar{\mathcal{E}}^{R}_{t}\approx\begin{bmatrix}0&\cdots&0\\\
\vdots&\vdots&\vdots\\\ 0&\cdots&\bar{\epsilon}_{2M_{t}}^{R}\end{bmatrix}.$
And thus, it follows that
$\mathcal{E}^{R}_{t}=\bar{\epsilon}_{2N}^{R}H_{M_{t}}H^{N_{t}}$, which implies
that:

$||\mathcal{E}^{R}_{t}||\leq|\bar{\epsilon}_{2M_{t}}^{R}|||H_{M_{t}}||||H^{N_{t}}||\leq
C^{\prime}|\bar{\epsilon}_{2M_{t}}^{R}|,$ (46)

using the assumption that $||H_{M_{t}}||||H^{N_{t}}||\leq C^{\prime}$. Hence,
using (45) and (46), it follows that if

$CC^{\prime}|\bar{\epsilon}_{2M_{t}}^{R}|\leq\beta<1,$ (47)

then $||(G^{R}_{t})^{-1}||\leq\frac{C}{R(1-\beta)}$, i.e.
$||P^{R}_{t}||\leq\frac{\sigma_{v}^{2}C}{R(1-\beta)}$ Thus, it suffices to
have:

$R\geq\frac{\sigma_{v}^{2}C}{\delta(1-\beta)},$ (48)

for $||P^{R}_{t}||\leq\delta.$ However, note that condition (47) still needs
to be satisfied.
Next, using the tail bounds for a Gaussian random variable [23], we have that
$Prob(|\bar{\epsilon}_{2M_{t}}^{R}|>\frac{n\sigma_{2M_{t}}}{\sqrt{R}})\leq
2e^{-n^{2}/2}$. Thus, if we require that $R$ is such that:

$\frac{n\sigma_{2M_{t}}}{\sqrt{R}}\leq\frac{\beta}{CC^{\prime}},$ (49)

then owing to the tail bounds, we have that:
$Prob(||P^{R}_{t}||\leq\delta)>1-2e^{-n^{2}/2}$. Thus, using the conditions
(47) and (49), if choose
$R=\max[\frac{\sigma_{v}^{2}C}{\delta(1-\beta)},(\frac{nCC^{\prime}}{\beta})^{2}\sigma_{2M_{t}}^{2}]$,
the result follows. ∎

###  VII-E Proof of Lemma 1

###### Proof:

It is clear that
$\delta\chi_{t}^{R}\delta\chi_{t}^{R^{\prime}}=R\mathcal{G}^{R}$. Next,
$\delta x_{t+1}^{(i)}=\sum_{k}F_{t}^{k}(\delta x_{t}^{(i)})^{k}+v_{t}^{(i)}$,
and $\delta X_{t+1}^{R}\delta\chi_{t}^{R^{\prime}}=[d^{1},\cdots d^{M}]$,
where $d^{l}\equiv\sum_{k=1}^{M}\sum_{i=1}^{R}F_{t}^{k}(\delta
x_{t}^{(i)})^{k}(\delta
x_{t}^{(i)})^{l}+\sum_{k>M}\sum_{i=1}^{R}F_{t}^{k}(\delta
x_{t}^{(i)})^{k}(\delta x_{t}^{(i)})^{l}$. However,
$\sum_{k>M}\sum_{i=1}^{R}F_{t}^{k}(\delta x_{t}^{(i)})^{k}(\delta
x_{t}^{(i)})^{l}=R\Delta_{t}^{R,l}$. Thus, using the above identities, it
follows that: $\delta
X_{t+1}^{R}\delta\chi_{t}^{R^{\prime}}=\mathcal{F}_{t}(R\mathcal{G}^{R})+R\Delta_{t}^{R}+V_{t}^{R}\delta\chi_{t}^{R^{\prime}}$,
and substituting this into (27), the result follows. ∎

###  VII-F Proof of Proposition 3

First, we need the following basic result.

###### Lemma 2

Let $h(z)=\sum_{k}H^{k}z^{k}$ be the Taylor series expansion of the function
$h(\cdot)$, and let it be convergent for any finite $z$. Let the lowest power
in the expansion be $L$. Then, given any $M\geq L$, and $\epsilon>0$, there
exists a $Z$ such that if $z<Z$,
$|\sum_{k>M}H^{k}z^{k}|\leq\epsilon|\sum_{k=1}^{M}H^{k}z^{k}|$.

###### Proof:

Let $S_{M}(z)=\sum_{k=1}^{M}H^{k}z^{k}$, and $\delta
S_{M}(z)=\sum_{k>M}H^{h}z^{k}$. Since the series is convergent, given any
$\epsilon>0$, there exists an $M^{\epsilon}(z)<\infty$, such that for all
$M>M^{\epsilon}(z)$, $\delta S_{M}(z)\leq\epsilon S_{M}(z)$. Let
$\bar{M}^{\epsilon}(z)=\max_{z^{\prime}\leq z}M^{\epsilon}(z^{\prime})$. That
this value is finite is trivial. Then, it follows that if $M^{\prime}$ is
chosen such that $M^{\prime}>\bar{M}^{\epsilon}(z)$, then for any
$z^{\prime}<z$, $\delta S_{M^{\prime}}(z^{\prime})\leq\epsilon
S_{M^{\prime}}(z^{\prime})$.
Next consider a $Z$ such that $M^{\epsilon}(Z)>M$ for the first time. If this
never happens, then $Z=\infty$, and the result trivially holds for the entire
domain. Then, if $z<Z$, due to the above analysis, it follows that $\delta
S_{M}(z)\leq\epsilon S_{M}(z)$, for all $z<Z$. The only case remaining is if
$z=0$: this is never possible unless the lowest degree in the expansion is
greater than $M$ which is precluded by the assumption of there being lower
order terms than $M$. This concludes the proof of the result.
∎

Using the result above, we can now prove Proposition 3.

###### Proof:

Recall that $\mathcal{F}_{t}\mathcal{G}+\Delta_{t}=[f_{t}^{j}]$ where
$f_{t}^{j}=E[\delta f\delta x_{t}^{j}]$, and $\delta
f=\sum_{k>0}F_{t}^{k}\delta x_{t}^{k}$. Thus, it follows that
$\tilde{f}_{t}^{j}=\sum_{k=1}^{M}F_{t}^{k}E[\delta x_{t}^{j+k}]$, and
$\Delta_{t}^{j}=\sum_{k>M}F_{t}^{k}E[\delta x_{t}^{j+k}]$.
First, suppose that $j$ is odd, then:
$\tilde{f}_{t}^{j}=\sum_{k=1}^{M}F_{t}^{k}C_{j+k}\sigma_{X}^{j+k}$ where
$E[\delta x_{t}^{p}]=C_{p}\sigma_{x}^{p}$, for some constant $C_{p}$. And
similarly $\Delta_{t}^{j}=\sum_{k>M}F_{t}^{k}C_{j+k}\sigma_{X}^{j+k}$. Thus,
$f_{t}^{j}(\sigma_{X})=E[\delta f\delta
x_{t}^{j}]=\sigma_{X}^{j}\sum_{k}F_{t}^{k}C_{j+k}\sigma_{X}^{k}$. Note that if
the input distribution is zero mean Gaussian, then we only have even powers of
$\sigma_{X}$. The above is simply a Taylor series expansion of the function
$f_{t}^{j}(\sigma_{X}^{2})$. Thus, using Lemma 2, it follows that there exists
$(\sigma_{X}^{j})^{2}<\infty$ such that
$|\Delta_{t}^{j}(\sigma^{2})|\leq\epsilon|\tilde{f}_{t}^{j}(\sigma^{2})|$, for
all $\sigma^{2}\leq(\sigma_{X}^{j})^{2}$. A similar argument holds for even
$j$.
Next, choose $\sigma_{X}^{2}=\min_{j}(\sigma_{X}^{j})^{2}$. Then, using the
above result, it follows that:
$|\Delta_{t}^{j}|\leq\epsilon|\tilde{f}_{t}^{j}|$, for all $j=1,..M$, thereby
proving the result.
∎

###  VII-G Proof of Corollary 2

###### Proof:

Note that $\delta\chi_{t}^{R}\delta\chi_{t}^{R^{\prime}}=R\mathcal{G}^{R}$,
and
$\mathcal{F}_{t}^{R}=\mathcal{F}_{t}+R\Delta_{t}^{R}(\frac{(\mathcal{G}^{R})^{-1}}{R})+\frac{1}{R}V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}$.
Defining
$S\equiv\Delta_{t}^{R}(\mathcal{G}^{R})^{-1}-\Delta_{t}\mathcal{G}^{-1}$, and
noting that $\mathcal{F}_{t}$, $S$ etc. are all row vectors, we can show that:

$\displaystyle
E||\mathcal{F}_{t}^{R}-(\mathcal{F}_{t}+\Delta_{t}\mathcal{G}^{-1})||^{2}=E||S||^{2}+\frac{2}{R}E[S(V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1})^{\prime}]+\frac{1}{R^{2}}E[V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}(\mathcal{G}^{R})^{-1}\delta\chi_{t}^{R}V_{t}^{R^{\prime}}].$
(50)

Since $\mathcal{G}^{R}\rightarrow\mathcal{G}$,
$\Delta_{t}^{R}\rightarrow\Delta_{t}$ in m.s.s. as $R\rightarrow\infty$, and
$\mathcal{G},\Delta_{t}$ are finite matrices, $E||S||^{2}\rightarrow 0$ as
$R\rightarrow\infty$, while
$E[S^{\prime}V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}]=0$,
since $V_{t}^{R}$ is a zero mean sequence, and independent from
$\delta\chi_{t}^{R}$, and hence from $S$ and $\mathcal{G}^{R}$. The last term
can be written as:

$\displaystyle\frac{1}{R^{2}}E[V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}(\mathcal{G}^{R})^{-1}\delta\chi_{t}^{R}V_{t}^{R^{\prime}}]$
$\displaystyle=\frac{1}{R^{2}}\text{Tr}(E[(\mathcal{G}^{R})^{-1}\delta\chi_{t}^{R}V_{t}^{R^{\prime}}V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}])$
(51)
$\displaystyle=\frac{1}{R^{2}}\text{Tr}(E_{\mathcal{G}^{R}}E[(\mathcal{G}^{R})^{-1}\delta\chi_{t}^{R}V_{t}^{R^{\prime}}V_{t}^{R}\delta\chi_{t}^{R^{\prime}}(\mathcal{G}^{R})^{-1}|\mathcal{G}^{R}]])$
(52)
$\displaystyle=\text{Tr}([E_{\mathcal{G}^{R}}[\frac{\sigma_{v}^{2}}{R}(\mathcal{G}^{R})^{-1}]),$
(53)

where $\text{Tr}(A)$ denotes the trace of a matrix $A$, and
$E_{\mathcal{G}^{R}}[.]$ represents the expectation with respect to
$\mathcal{G}^{R}$. Next, we know that $\mathcal{G}^{R}\rightarrow\mathcal{G}$
almost surely as $R\rightarrow\infty$, and since $\mathcal{G}$ is a finite
matrix, it follows that
$\text{Tr}([E_{\mathcal{G}^{R}}[\frac{\sigma_{v}^{2}}{R}(\mathcal{G}^{R})^{-1}])\rightarrow
0$ as $R\rightarrow\infty$. Therefore,
$E||\mathcal{F}_{t}^{R}-(\mathcal{F}_{t}+\Delta_{t}\mathcal{G}^{-1})||^{2}\rightarrow
0$ as $R\rightarrow\infty$, proving the result. ∎

###  VII-H Model Free Solution of the PPE

In the model-free method, rather than estimating $\mathcal{F}_{t}$ and
substituting into $\mathcal{B}_{t}$ in (26), one can directly try to estimate
the Taylor coefficients $\mathcal{K}_{t}$ given those at the next time
$\mathcal{K}_{t+1}$. A Least squares estimate of the same is easily seen to
be:

$\mathcal{K}_{t}^{R}=\mathcal{C}_{t}+\mathcal{K}^{R}_{t+1}(\delta\chi_{t+1}\delta\chi_{t}^{\prime})(\delta\chi_{t}\delta\chi_{t}^{\prime})^{-1},$
(54)

where note that unlike in the model-based case where we use $\delta X_{t+1}$,
we use
$\delta\chi_{t+1}^{R}=\begin{bmatrix}\delta x_{t+1}^{(1)}&\cdots&\delta
x_{t+1}^{(R)}\\\ \vdots&\vdots&\vdots\\\ (\delta
x_{t+1}^{(1)})^{M}&\cdots&(\delta x_{t+1}^{(R)})^{M}\end{bmatrix}$. In
essence, the correlation $\delta\chi_{t+1}\delta\chi_{t}^{\prime}$ allows one
to deduce the matrix $\mathcal{B}_{t}$ in the PPE equation (26) rather than
using the known structure of $\mathcal{B}_{t}$.

The accuracy and variance of the solution are amenable to an identical
analysis as we did for the model-based case, and it should not come as a
surprise that even in this case, it is determined by the convergence and
variance of the empirical Gram matrix $\mathcal{G}^{R}$.

Model-based or Model-free?
The key to computing the solution is doing the product
$\delta\chi_{t+1}\delta\chi_{t}$, the number of computations are
$\mathcal{O}(2R^{2}M^{2})$, while that in the model-based case for the product
$\delta X_{t+1}\delta\chi_{t}^{\prime}$ is $\mathcal{O}((R^{2}+R)M^{2})$.
Since $R$ is typically going to be large, this shows that we need to do
approximately $2\times$ the computation to get the answer in the model-free
method, since it is trying to infer the matrix $\mathcal{B}_{t}$ rather than
using the inherent structure. Further, given the same number of computations,
we would be able to do $\approx 2R$ simulations to the $R$ of the model-free
method, and thus, the variance in the empirical Gram matrix $\mathcal{G}^{R}$
is $\frac{1}{2}$ times that of the model-free empirical Gram matrix. Again,
this shows that we should prefer the model based approach, either to save
computation or reduce the variance of the solution.

##  VIII The Discounted Infinite Horizon Case

Thus far, we have considered the finite horizon optimal control case. Now, we
consider the more prevalent discounted infinite horizon case. Again, we
consider the deterministic case first, and then the stochastic case. It is
always possible to turn any discounted infinite horizon case into a finite
horizon problem with a sufficiently long horizon, and a zero terminal cost. In
this section, we concentrate on the infinite horizon policy evaluation
problem, given the dynamics $x_{t+1}=f(x_{t})$. In essence, we want to
evaluate the cost of the above dynamics in the discounted case. To that end,
let $c(x)$ denote the incremental cost and let $\beta<1$ denote a discount
factor. the discounted cost-to-go given an initial condition $x_{0}$ is given
by: $J(x_{0})=\sum_{t=0}^{\infty}\beta^{t}c(x_{t})$, where $x_{t}$ evolves
according to the dynamics above. The cost function $J$ obeys the following
stationary equation, given any state $x$:

$J(x)=c(x)+\beta J(f(x)).$ (55)

Next, consider a nominal trajectory, given some initial condition $x_{0}$, say
$\bar{x}_{t+1}=f(\bar{x}_{t})$. Next, let us do expansions of the cost about
this nominal trajectory, as we did previously, $J(x_{t})=J(\bar{x}_{t}+\delta
x_{t})=J(\bar{x}_{t})+K_{t}^{1}\delta x_{t}+K_{t}^{2}\delta x_{t}^{2}+\cdots$,
$c(x_{t})=c(\bar{x}_{t}+\delta x_{t})=c(\bar{x}_{t})+C_{t}^{1}\delta
x_{t}+C_{t}^{2}\delta x_{t}^{2}+\cdots$, and

$J_{t+1}(f(x_{t}))=J_{t+1}(f(\bar{x}_{t})+\delta
f_{t})=J_{t+1}(f(\bar{x}_{t}))+K_{t+1}^{1}\delta f_{t}+K_{t}^{2}\delta
f_{t}^{2}+\cdots,$ (56)

where $\delta f_{t}=F_{t}^{1}\delta x_{t}+F_{t}^{2}\delta x_{t}^{2}+\cdots$.
Substituting into (55), and equating the different powers of $\delta x_{t}$ on
both sides, we obtain:

$\displaystyle\bar{J}_{t}=\bar{c}_{t}$ $\displaystyle+\beta\bar{J}_{t+1},$
$\displaystyle K_{t}^{1}=C_{t}^{1}$ $\displaystyle+\beta
K_{t+1}^{1}F_{t}^{1},$ $\displaystyle K_{t}^{2}=C_{t}^{2}$
$\displaystyle+\beta(K_{t+1}^{1}F_{t}^{2}+K_{t+1}^{2}(F_{t}^{1})^{2}),$ (57)
$\displaystyle\vdots$

where $\bar{J}_{t}=J(\bar{x}_{t})$, $\bar{c}_{t}=c(\bar{x}_{t})$, and
$\bar{J}_{t+1}=J_{t+1}(f(\bar{x}_{t}))$, with terminal condition $J_{T}(x)=0$
for any $x$. Again, note the beautiful perturbation structure of the above
equations in that the higher order terms do not affect the calculations of the
lower order terms, and thus, we may close our computations at any desired
order.
However, one has to be careful about using the above perturbation equations.
An implicit assumption here is that $T$ is large enough such that $J_{t}(x)$
becomes stationary, i.e., converges to some $J_{\infty}(x)$, $J_{t}(x)\approx
J_{\infty}(x)$, for all $t<\bar{t}$ where $T-\bar{t}\ll T$. Moreover, we
should use the solution to the above equations (VIII) only after the initial
transient time $\bar{t}$ has passed, i.e., we use it only $J_{t}$ for
$t<\bar{t}$. The situation is illustrated in Fig. 4. Further, we should note
that albeit $J_{\infty}(x)$ is time invariant, when we expand it locally
around a nominal trajectory, the result expansions are time-varying (again,
see Fig. 4).

Figure 4: The Discounted Infinite Horizon Case

Moreover, just like in the finite horizon case, One may also choose not to do
a time varying decomposition about a trajectory, and instead do an expansion
about the time invariant trajectory, $\bar{x}(t)=0$, assuming that $f(0)=0$.
However, identical to the situation in the finite horizon case, most initial
conditions $x_{0}$ would be far from the origin, and thus, we would need far
more terms to make an accurate approximation when compared to the time varying
case, thereby drastically increasing the computational burden of the method.
The policy evaluation can also be done in an RL fashion as was done previously
for the finite horizon case in an almost identical fashion by estimating
$\mathcal{F}_{t}$, or directly by using $\delta\chi_{t+1}$, the only caveat
again being that we use the solution only after it has equilibriated, i.e.,
for $t<\bar{t}$. The stochastic case remains intractable owing again to the
lack of a perturbation structure, and thus, any RL solution is bound to have
high variance, and thus, be inaccurate.

###### Remark 3

The Method of Characteristics. The equation 6, when written in continuous time
becomes the PDE: $\frac{\partial J}{\partial t}+c+f\frac{\partial J}{\partial
x}=0$, where the dynamics are now in continuous time, i.e., $\dot{x}=f(x)$ and
the cost function is given by: $\int_{0}^{T}c(x_{t})dt$. The classical method
of characteristics reduces the above PDE into a family of ordinary
differential equations (ODE) called the characteristic ODEs/ Lagrange-Charpit
equations, in terms of the state $x$ and the co-state $q=\frac{\partial
J}{\partial x}$, given terminal conditions $\bar{x}_{T}$ and
$\bar{q}_{t}=\frac{\partial g}{\partial x}|_{\bar{x}_{T}}$ [22] (see Section
VII B). The PPE equations (26) are the discrete time analogs of the
perturbation expansion of the characteristic ODEs about the nominal
characteristic curve $(\bar{x}_{t},\bar{q}_{t})$. Next, consider the infinite
horizon equation (55): this is the discrete time analog of the PDE
$\frac{\partial J}{\partial t}=c+\beta f\frac{\partial J}{\partial x}$, where
the cost function now is given by $\int_{0}^{\infty}c(x_{t})e^{-\beta t}dt$.
In this case, the perturbation equations (VIII) are the discrete time analogs
of the perturbation expansions of the characteristic ODEs around the nominal
path $\dot{\bar{x}}=f(\bar{x})$. In this case, an additional condition is that
after a suitable amount of time, the solution equilibriates, i.e.,
$\frac{\partial J_{\infty}}{\partial t}=c+\beta f\frac{\partial
J_{\infty}}{\partial x}=0$. The stationary PDE is difficult to solve and
typically is solved by evolving the time varying PDE till the solution becomes
stationary, which is exactly what is done when using (VIII) in discrete time.

###### Remark 4

Koopman Operator. In recent years, there has been an increasing interest in
the Koopman operator approach to the study of nonlinear systems [24]. The idea
is to consider the recursive dynamic map $\mathcal{K}_{f}\cdot g\equiv
g(f(x))$, where $g(\cdot)$ is termed an observable and to solve for this
recursive map in a data based fashion. We note that the Koopman recursive map
is equivalent to the policy evaluation equation (6), with the cost $c(x)=0$,
and the terminal condition given by $g(x)$. Thus, any data based solution to
the Koopman problem also suffers from the same variance and convergence issues
as those afflicting policy evaluation, when it is solved in an RL/ data based
fashion.