Approximate Dynamic Programming And The Krylov Subspace Methods

1. Introduction

Dynamic programming (DP) is both an optimization method and a computer programming method which was developed by Richard Bellman in the 1950s. The main idea is to break a complicated problem down into simpler sub-problems recursively. If there exists a recursive relation between the original problem and its sub-problems, then DP methods are applicable. In the optimization literature this relationship is called the Bellman equation.

In machine learning, one of the problems that we are interested in is the sequential decision making problem regarding the way in which the software agents ought to take actions to maximize the cumulative rewards from the stochastic environment. We formulate the interaction bewteen agents and the environment by Markov Decision Process (MDP) whose structure can be modeled by the Bellman equation. But when the state space of the environment is too large or an exact model of the dynamics is absent, MDP can only be solved approximately with simulation and function approximation and this is what we call approximate dynamic programming (neuro-dynamic programming or reinforcement learning). A large body of methods have been developed, e.g. temporal difference learning (TD) [4] , least-squares temporal difference learning (LSTD) [3,4], least-squares policy evaluation (LSPE) [9] and least-squares policy iteration (LSPI) [6].

Computing the value function of a given policy (policy evaluation) by the Bellman equation plays an important role in solving approximate DP and it is equivalent to solving a linear system [13]. In matrix computations, Krylov subspace methods [5,8] are the most powerful tools to solve large scale linear systems of equations and least-squares problems. So there is a connection between these two areas but this connection lacks deep exploration.

Instead of using reinforcement learning, we use approximate DP in this report to emphasize its strong relationship with DP, because DP provides theoretical foundations and important guidelines for designing approximate DP algorithms.

My report will review some fundamental concepts and ideas of approximate DP and Krylov subspace methods and their connections.

2. Approximate dynamic programming methods

2.1. Problem description

Markov decision process (MDP) is a framework to model the learning process that the agent learns from the interaction with the environment. Given the current state $S_t\in \mathcal{S}$ at time t, the agent interacts with the environment by selecting actions $A_t\in \mathcal{A}\left(S_t\right)$ according to the policy $\pi \left(A_t\right|S_t)$ and the environment would take the agent to another state $S_{t+1}\in \mathcal{S}$, and then provide a numerical feedback (immediate reward) $R_{t+1}=R(S_t,A_t,S_{t+1})\in \mathbb{R}$ according to the dynamics (transition probability) $p(S_{t+1},R_{t+1}|S_t,A_t)$. The policy $\pi$ is a probability distribution from which we select the action. Our goal is to find an optimal policy (often deterministic) to maximize the long-term cumulative rewards. We assume that the policy is stationary. We only consider the inifite horizon discounted problem because the ideas can also be applied to the simpler absorbing Markov chains [17]. So the objective function is $G_t={\displaystyle \sum _{k=0}^{\infty }}{\gamma }^k{R}_{t+k+1}$. It’s easy to see the recursive structure

$$G_t=R_{t+1}+\gamma G_{t+1}.$$

(1)

The value function of the given policy $\pi$ at state s is defined as

$$v_{\pi }\left(s\right)=E_{\pi }\left[G_t\right|S_t=s]=E_{\pi }\left[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s\right],\forall s\in \mathcal{S}$$

(2)

It can be written as

$$\begin{array}{cc}\hfill v_{\pi }\left(s\right)& =E_{\pi }[R_{t+1}+\gamma G_{t+1}|S_t=s]\hfill \\ & =\sum _a\pi (A_t=a|S_t=s)\sum _{s^{\prime }}\sum _{r}p(S_{t+1}=s^{\prime },R_{t+1}=r|S_t=s,A_t=a)\left[r+\gamma E_{\pi }\left[G_{t+1}\right|S_{t+1}=s^{\prime }]\right]\hfill \\ & =\sum _a\pi \left(a\right|s)\sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+\gamma v_{\pi }\left(s^{\prime }\right)],\forall s\in \mathcal{S}\hfill \end{array}$$

(3)

The Equation (3) is the Bellman equation for MDP, which depicts the recursive relationships between the value of s and the value of its possible successor states $s^{\prime }$ and is consistent to (1). It can be written in a more compact form.

$$\begin{array}{cc}\hfill v_{\pi }\left(s\right)& =\sum _a\pi \left(a\right|s)\sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+\gamma v_{\pi }\left(s^{\prime }\right)]\hfill \\ & =\sum _a\pi \left(a\right|s)\sum _{s^{\prime },r}r·p(s^{\prime },r|s,a)+\sum _a\pi \left(a\right|s)\sum _{s^{\prime },r}p(s^{\prime },r|s,a)\gamma v_{\pi }\left(s^{\prime }\right)\hfill \\ & =\sum _a\pi \left(a\right|s)r(s,a)+\gamma \sum _{s^{\prime }}\{\sum _a\pi \left(a\right|s)·p\left(s^{\prime }\right|s,a)\}{v}_{\pi }\left(s^{\prime }\right)\hfill \\ & =r_{\pi }\left(s\right)+\gamma \sum _{s^{\prime }}{p}_{\pi }\left(s^{\prime }\right|s)v_{\pi }\left(s^{\prime }\right),\forall s\in \mathcal{S}\hfill \end{array}$$

(4)

The Equation (4) can be written in the matrix form

$$v_{\pi }=r_{\pi }+\gamma P_{\pi }{v}_{\pi }$$

(5)

where $v_{\pi }\in {\mathbb{R}}^{\left|\mathcal{S}\right|},r_{\pi }\in {\mathbb{R}}^{\left|\mathcal{S}\right|}$, $P_{\pi }\in {\mathbb{R}}^{\left|\mathcal{S}\right|\times \left|\mathcal{S}\right|}$ and $P_{\pi }(i,j)=p_{\pi }\left(j\right|i)$. Here $\left|\mathcal{S}\right|$ can be infinite, but for simplicity, we assume it is finite. Since we only consider a fixed policy, we can write (5) as $v=r+\gamma Pv$ or $Lv=r\mathrm{where}L=I-\gamma P$.

We define a short hand notation for DP mapping:

$$Tv=r+\gamma Pv,\forall v\in \mathcal{V}$$

(6)

where $\mathcal{V}$ is the value function space and T is called the Bellman operator. It can be proved that T is a contraction mapping [1] and its fixed point $v^{*}$ is the solution for (5). Finding the fixed point is exactly what we do in the policy evaluation step in the policy iteration method [13].

Solving (5) by a direct linear solver is not feasible when the state space is very large or when the transition matrix is unavailable. But we can find approximate solutions by value function approximation (VFA). Linear approximation is a popular VFA method [11] and it has the form of $\widehat{v}=\Phi \theta ,\Phi =({\Phi }_1,{\Phi }_2,\cdots ,{\Phi }_K)\in {\mathbb{R}}^{\left|\mathcal{S}\right|\times K},$ where ${\Phi }_j$ is the j-th feature function or basis function over the state space and K is the number of features, see [14]. For each state i, $\widehat{v}\left(i\right)={\displaystyle \sum _{j=1}^K}{\Phi }_j\left(i\right){\theta }_j=\varphi {\left(i\right)}^T\theta$, where $\varphi {\left(i\right)}^T$ is the $i^{th}$ row of $\Phi$ and represent the feature vector for state i. Linear function approximations have better convergence property [17] than nonlinear function approximation and is easier to analyze, interpret and implement. We can introduce nonlinearity when doing feature extraction, e.g. deep neural networks. Suppose we already have well-chosen features, then by the linear approximation, (5) becomes

$$\Phi \theta =r+\gamma P\Phi \theta .$$

(7)

We define a $\xi -$weighted quadratic norm on $\mathcal{V}$,

$${||v||}_{\xi }=\sqrt{\sum _{s\in \mathcal{S}}\xi \left(s\right){\left(v\left(s\right)\right)}^2}=\sqrt{v^T\Xi v},\forall v\in \mathcal{V}$$

(8)

where $\xi ={\left({\xi }_i\right)}_{i=1}^{\left|\mathcal{S}\right|}$ with ${\xi }_i>0$ is the weight vector and $\Xi$ is the diagonal matrix with $\xi$ on the diagonal. Then we can define the operator $\Pi$ that projects any value function to its closest representable function in range$(\Phi )$ (i.e. column space of $\Phi$) with respect to $\Xi$. That means $\Pi$ is an orthogonal projector onto $range(\Phi )$ with respect to $\Xi$, i.e. ${\Phi }^T\Xi (v-\Pi v)=0$, which is equivalent to

$$\Pi v\doteq \Phi {\theta }^{\prime },\mathrm{where}{\theta }^{\prime }=\underset{\theta \in {\mathbb{R}}^K}{argmin}{||v-\Phi \theta ||}_{\xi }$$

(9)

We assume that the Markov chain we are studying is ergodic, then we have the stationary distribution,

$$\mu \left(s\right)=\underset{N\leftarrow \infty }{lim}\frac{1}{N}\sum _{k=1}^{N}P(S_k=s|S_0=s_0)>0,\forall s,s_0\in \mathcal{S}$$

(10)

We usually use the vector $\mu$ as our weight vector $\xi$, otherwise it may cause a norm mismatch problem which will make the iterative methods described in the following sections diverge [1,14].

Bertsekas [1] gives a short summary of the main approximate DP algorithms:

• Direct policy evaluation: Find $\theta$ that minimize ${||v-\Phi \theta ||}_{\xi }^2$. We compute the approximation of the projection by Monte Carlo simulation. We can introduce nonlinear architectures in $\widehat{v}$, e.g. deep neural networks.
• Indirect policy evaluation:

  -

  Solve the projected equation (Galerkin approximation [16,18])

  $$\Phi \theta =\Pi T(\Phi \theta )$$

  (11)

  by simulation. Equation (11) is equivalent to minimizing the expected temporal difference error (TDE) ||${\widehat{v}-\Pi T\widehat{v}}_{\xi }$||_ξ down to 0 for $\widehat{v}\in range(\Phi )$ [14].

  *

  TD learning: Stochastic iterative algorithm for solving (11).

  *

  Least-squares temporal difference learning (LSTD): TD fixed point method [14].

  *

  Least-squares policy evaluation (LSPE): A simulation-based form of projected iteration; essentially $\Phi {\theta }_{k+1}=\Pi T(\Phi {\theta }_k)+$simulation noise.

  -

  Bellman equation error (Bellman residual [14]) method: minimize ${||\Phi \theta -T(\Phi \theta )||}_{\xi }$

2.2. Direct policy evaluation

We are interested in finding the ${\theta }^{*}=\underset{\theta \in {\mathbb{R}}^K}{argmin}{||v-\widehat{v}||}_{\xi }$ with $\widehat{v}=\Phi \theta$. Suppose we can sample the exact values or give an accurate estimation of v, then we can learn the parameters in supervised learning paradigm. The samples or estimation of v could be considered as the true label, and $\widehat{v}$ is the approximator. Then, it becomes a standard weighted least-squares problem. The solution is the projection of v onto $range(\Phi )$ with respect to $\Xi$ [1]. Then,

$${\theta }^{*}=\underset{\theta \in {\mathbb{R}}^K}{argmin}{||v-\Phi \theta ||}_{\xi }^2=\underset{\theta \in {\mathbb{R}}^K}{argmin}\sum _{i=1}^{\left|\mathcal{S}\right|}{\xi }_i{(v\left(i\right)-\varphi {\left(i\right)}^T\theta )}^2$$

Set the gradient of $\theta$ to 0, we have

$${\theta }^{*}={\left(\sum _{i=1}^{\left|\mathcal{S}\right|}{\xi }_i\varphi \left(i\right)\varphi {\left(i\right)}^T\right)}^{-1}\sum _{i=1}^{\left|\mathcal{S}\right|}{\xi }_i\varphi \left(i\right)v\left(i\right)$$

Suppose we have samples $\{(s_1,v\left(s_1\right)),\cdots ,(s_m,v\left(s_m\right))\}$ generated by the environment or a simulator of the environment according to $\xi$, where $s_j\in \{1,\cdots ,|\mathcal{S}\left|\right\}\mathrm{for}j=1:m,$ then we can get the simulation-based approximation:

$${\widehat{\theta }}^{*}=\underset{\theta \in {\mathbb{R}}^K}{argmin}\sum _{j=1}^m{(v\left(s_j\right)-\varphi {\left(s_j\right)}^T\theta )}^2={\left(\sum _{j=1}^m\varphi \left(s_j\right)\varphi {\left(s_j\right)}^T\right)}^{-1}\sum _{j=1}^m\varphi \left(s_j\right)v\left(s_j\right)$$

Note that $\varphi \left(s_j\right)\varphi {\left(s_j\right)}^T$ is a $K\times K$ matrix. Thus, we are using low-dimensional linear algebra (of order K, the number of basis functions) compared with (5).

We can obtain the unbiased estimate of the true value function v by Monte Carlo sampling, i.e. perform long-term lookahead through the Markov chain and collect all the future rewards. However, the Monte-Carlo estimate of v requires a large number of samples as it suffers from a high variance.

2.3. Indirect policy evaluation

2.3.1. Stochastic algorithms: a general view

Consider the approximate solution ${\widehat{x}}_m$ of a linear system of equations $x=b+Ax,A\in {\mathbb{R}}^{n\times n}$ by using m simulation samples $b+w_k$ and $A+W_k,k=1,\cdots ,m,$ where $w_k,W_k$ are simulation noise. Think of $x=b+Ax$ as approximate policy evaluation which is the same as (5) [1]:

• Stochastic approximation (SA) approach: starting from an initial point $x_0$, for $k=1,\cdots ,m$

  $$x_k=(1-{\alpha }_k)x_{k-1}+{\alpha }_k((b+w_k)+(A+W_k)x_{k-1})$$

  (12)
  ${\alpha }_k,k=1,\cdots ,m$ are the weights.
• Monte Carlo estimation (MCE) approach: Form Monte Carlo estimation of b and A,

  $$b_m=\frac{1}{m}\sum _{k=1}^m(b+w_k),A_m=\frac{1}{m}\sum _{k=1}^m(A+W_k)$$

  (13)
  Then solve $x=b_m+A_{m}x$ by a direct method or an iterative method. $A_m,b_m$ can be built incrementally.

TD is SA method, LSTD and LSPE are MCE methods. We will introduce them in section 2.3.2-2.3.4.

2.3.2. Temporal difference (TD) learning

In 1988, Sutton developed temporal difference learning as a class of incremental learning methods for general prediction problems. It assigns credit through the difference of the predictions between current state and its successor, which makes it different with the conventional supervised learning methods. It has the form [4]

$$v\left(S_t\right)=v\left(S_t\right)+\alpha [R_{t+1}+\gamma v\left(S_{t+1}\right)-v\left(S_t\right)]$$

(14)

where $\alpha$ is a constant step-size parameter. This form is equivalent to (12). TD learning is a combination of Monte Carlo ideas and DP ideas. It can learn directly from the experience $R_{t+1}$ sampled from environment which is similar to Monte Carlo methods. It also updates the estimates partially based on other learned estimates $v\left(S_{t+1}\right)$ without collecting all future rewards (bootstrapping) which is similar to DP.

When we use linear function approximation $\widehat{v}\left(s_t\right)=\varphi {\left(s_t\right)}^T{\theta }_t$, the parameter update can be written as [4],

$${\theta }_{t+1}={\theta }_t+\alpha [R_{t+1}+\gamma \varphi {\left(s_{t+1}\right)}^T{\theta }_t-\varphi {\left(s_t\right)}^T{\theta }_t]\varphi \left(s_t\right)$$

(15)

Equation (15) is like we are doing ${min}_{\theta }J\left(\theta \right)={(v\left(s_t\right)-\widehat{v}\left(s_t\right))}^2$. We apply gradient descent ${\theta }_{t+1}={\theta }_t-\frac{1}{2}\alpha \nabla J\left(\theta \right)={\theta }_t+\alpha (v\left(s_t\right)-\widehat{v}\left(s_t\right))\nabla \widehat{v}\left(s_t\right)$ and substitute $v\left(s_t\right)$ by its approximation $R_{t+1}+\gamma \varphi {\left(s_{t+1}\right)}^T{\theta }_t$. We call them semi-gradient methods because we only take derivative of $\widehat{v}\left(s_t\right)$ but not of $R_{t+1}+\gamma \varphi {\left(s_{t+1}\right)}^T{\theta }_t$. TD is proved to be convergent with probability one [17].

TD requires less memory, less peak computation and produce more accurate predictions than conventional supervised learning methods.

2.3.3. Least-squares temporal difference (LSTD) learning

As $\Pi T$ is a contraction mapping on $range(\Phi )$ with respect to $\Xi$ [1,14], there exists an unique fixed point $\Phi {\theta }^{*}$ that satisfies (11). That means ${\theta }^{*}$ is the unique solution to

$$\Phi \theta =\Pi T(\Phi \theta )$$

(16)

which means $\Phi {\theta }^{*}$ is the unique point in $range(\Phi )$ that satisfies the condition that when we use the projector $\Pi T$ on it, the projection is itself. According to the definition of $\Pi$, ${\theta }^{*}$ should satisfies the condition

$${\Phi }^T\Xi (\Phi {\theta }^{*}-(r+\gamma P\Phi {\theta }^{*}))=0.$$

(17)

The orthogonal condition (17) is equivalent to the Galerkin method [16,18]. Equation (17) can be written as a linear systems of equations $C{\theta }^{*}=d$, where

$$C={\Phi }^T\Xi (I-\gamma P)\Phi ,d={\Phi }^T\Xi r$$

(18)

So we can have ${\theta }^{*}=C^{-1}d$. The simulation machanism for C and d is

$$d_k=\frac{1}{k+1}\sum _{t=0}^k\varphi \left(s_t\right)r_t\approx \sum _{i,j}{\xi }_{i}P(i,j)\varphi \left(i\right)R(i,j)={\Phi }^T\Xi r=d$$

(19)

$$C_k=\frac{1}{k+1}\sum _{t=0}^k\varphi \left(s_t\right){(\varphi \left(s_t\right)-\gamma \varphi \left(s_{t+1}\right))}^T\approx {\Phi }^T\Xi (I-\gamma P)\Phi =C$$

(20)

where $C_k\in {\mathbb{R}}^{K\times K},d_k\in {\mathbb{R}}^K$ can be computed with low-dimensional calculation. Then ${\theta }^{*}$ can be approximated by ${\widehat{\theta }}_k=C_k^{-1}{d}_k$.

LSTD makes more efficient use of data because it extract more information from each additional observation, and thus would be expected to converge with fewer training samples. In addition, it does not require the user to manually tune a stepsize schedule for good performance. LSTD converges to the same point as TD [4].

TD’s convergence can be slowed dramatically by a poor choice of the stepsize parameters $\alpha$. Furthermore, TD’s performance is sensitive to the initial estimate for $\theta$. LSTD does not rely on an arbitrary initial estimate. TD is also sensitive to the ranges of the individual features. LSTD is not.

2.3.4. Least-squares policy evaluation

An alternative way to solve the projected equation is projected value iteration (PVI), that is

$$\Phi {\theta }_{k+1}=\Pi T(\Phi {\theta }_k)=\Pi (r+\gamma P\Phi {\theta }_k)$$

(21)

It converges to ${\theta }^{*}$ because $\Pi T$ is a contraction mapping with respect to $\Xi$. Rather than directly solve the projected equation, it uses an iterative way to get close to the fixed point starting from an arbitrary initial point. It is equivalent to the minimization problem

$${\theta }_{k+1}=\underset{\theta \in {\mathbb{R}}^K}{argmin}{||\Phi \theta -(r+\gamma P\Phi {\theta }^k)||}_{\xi }^2$$

(22)

Set the gradient with respect to $\theta$ to be 0, we could have

$${\Phi }^T\Xi (\Phi {\theta }_{k+1}-(r+\gamma P\Phi {\theta }_k))=0$$

(23)

Then we could have

$${\theta }_{k+1}={\theta }_k-{({\Phi }^T\Xi \Phi )}^{-1}(C{\theta }_k-d)$$

(24)

where $C,d$ are defined as (18) and could be estimated by (19) and (20). $G={\Phi }^T\Xi \Phi$ could be estimated by

$$G_k=\frac{1}{k+1}\sum _{t=0}^k\varphi \left(s_t\right)\varphi {\left(s_t\right)}^T\approx {\Phi }^T\Xi \Phi =G$$

(25)

$C_k,d_k$ and $G_k$ can be formed incrementally or in temporal difference style [1].

2.3.5. Bellman equation error methods

The main idea of the Bellman equation error (BE) method is to find a $\widehat{v}\in range(\Phi )$ so that ${||\widehat{v}-T\widehat{v}||}_{\xi }$ can be minimized. Recall that $L=I-\gamma P$. With the notation $\Psi =L\Phi$, we have ${||\widehat{v}-T\widehat{v}||}_{\xi }={||\Phi \theta -r-\gamma P\Phi \theta ||}_{\xi }={||\Psi \theta -r||}_{\xi }$. Then it becomes a standard least squares problem with regard to ${||·||}_{\xi }$. It is easy to see that the solution to the weighted least squares problem is

$${\theta }_{BE}={({\Psi }^T\Xi \Psi )}^{-1}{\Psi }^T\Xi r$$

(26)

One reason that we are interested in BE method is the well-known result: $\forall \widehat{v},{||v-\widehat{v}||}_{\infty }\le \frac{1}{1-\gamma }{||T\widehat{v}-\widehat{v}||}_{\infty }$ and for weighted norm we have $\forall \widehat{v},{||v-\widehat{v}||}_{\xi }\le \frac{\sqrt{C\left(\xi \right)}}{1-\gamma }{||T\widehat{v}-\widehat{v}||}_{\xi }$, where $C\left(\xi \right):=ma{x}_{i,j}\frac{P(i,j)}{{\xi }_i}$ is a “concentration coefficient”, which can be seen as a measure of the stochasticity of the MDP [14]. This means when BE is small, we are close to the true value function v. But for TD, there is no similar result. Actually, it is proved [14] that BE is an upper bound of the TDE. More exactly, since $\widehat{v}=\Phi \theta =\Pi \Phi \theta$ and from Pythagorean Theorem,

$$B{E}^2={||\widehat{v}-T\widehat{v}||}_{\xi }^2={||\Pi (\widehat{v}-T\widehat{v})||}_{\xi }^2+{||(I-\Pi )(\widehat{v}-T\widehat{v})||}_{\xi }^2=TD{E}^2+{||T\widehat{v}-\Pi T\widehat{v}||}_{\xi }^2$$

(27)

which means minimizing BE is not only minimizing TDE, but also minimizing ${||T\widehat{v}-\Pi T\widehat{v}||}_{\xi }^2$, which is part that the TD method misses. And also we can see that $TD{E}^2={||\Pi (\widehat{v}-T\widehat{v})||}_{\xi }^2$ and this is why people also call it the projected Bellman error (PBE).

BE has the additional component which has to do with the features. Hence, if the feature space is fixed, the term we really want to minimize is TDE. Only when features can be modified we start to look at BE.

2.4. Least-squares policy iteration (LSPI)

Motivated by LSTD, Lagoudakis [6] extended it to control problems. Instead of evaluating the state value function v, he did evaluation on the state-action pairs ($s,a$), i.e. the Q function. The resulting method is referred to as LSTDQ. The Q function is defined as

$$Q(s,a)=E\left[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a\right],\forall s\in \mathcal{S},a\in \mathcal{A}\left(s\right)$$

(28)

We can directly choose action and improve our policy from the Q function without any knowledge of the underlying model of the environment. Everything about LSTDQ looks similar to LSTD after we change $\Phi$ to be the feature matrix for all state-action pairs. Then we put LSTD into the general policy iteration (GPI) frameworks to update our policy and get LSPI algorithm.

2.5. A unified oblique projection view

In the LSTD, we attempt to find the fixed point of the composite projector $\Pi T$, where $\Pi$ is an orthogonal projection onto $range(\Phi )$ with respect to $\Xi$. Now we consider using a projector $\Pi$ that is not necessarily orthogonal, which is a general operator that satisfies ${\Pi }^2=\Pi$ and whose range is $range(\Phi )$. In the most general case, such operators are called oblique projections and can be written ${\Pi }_X=\Phi {\pi }_X$ with ${\pi }_X={(X^T\Phi )}^{-1}{X}^T$. The parameter X specifies the projection direction: precisely, ${\Pi }_X$ is the projection onto $range(\Phi )$ orthogonal to $range\left(X\right)$ or we could write it as $\forall v,(v-{\Pi }_{X}v)\perp X$ (the $\Pi$ defined in (9) can be written as ${\Pi }_{\Xi \Phi }$). The vector ${\pi }_{X}v$ is the coordinate of the projection of v under the basis $\Phi$.

Scherrer [14] provided a unified oblique projection view of TD fixed point and BE minimization methods. Recall that ${\theta }_{TD}={({\Phi }^T\Xi L\Phi )}^{-1}{\Phi }^T\Xi r,{\theta }_{BE}={({\Psi }^T\Xi \Psi )}^{-1}{\Psi }^T\Xi r$. If we write $X_{TD}=\Xi \Phi$ and $X_{BE}=\Xi L\Phi$, we can have (1) The TD fixed point and the BE minimizer are solutions to the projected equation ${\widehat{v}}_X={\Pi }_{X}T{\widehat{v}}_X$ ($X=X_{TD}$ and $X=X_{BR}$ respectively). (2) When it exists, the solution of this projected equation is the projection of v onto $range(\Phi )$ orthogonally to $range\left(L^{T}X\right)$, i.e. formally ${\widehat{v}}_X={\Pi }_{L^{T}X}v.$ From this perspective, we put the two methods under the same framework of projected fixed point equation.

From the results above, it’s easy to have

$${||v-{\widehat{v}}_X||}_{\xi }\le {||{\Pi }_{L^{T}X}||}_{\xi }{||v-\Pi v||}_{\xi }=\sqrt{\sigma \left(ABC{B}^T\right)}{||v-\Pi v||}_{\xi },\forall X.$$

(29)

where $A={\Phi }^T\Xi \Phi ,B={(X^{T}L\Phi )}^{-1}$, $C=X^{T}L{\Xi }^{-1}{L}^{T}X$ are matrices of size $K\times K$ and $\sigma$ is the spectral radius of the given matrix. $A,B,C$ can be estimated by simulations. This inequality could also give us some inspiration on how to choose X and extract $\Phi$, e.g. if we choose $X={\left(L^T\right)}^{-1}\Xi \Phi$, we could have $\widehat{v}={\Pi }_{\Xi \Phi }v$, which is the orthogonal projection of true v onto $range(\Phi )$ with respect to $\Xi$ .

Finally, Scherrer gives a comparison of TD and BE methods: BE is more robust than TD because it has guarantee of its performance. Although TD often outperforms BE, it could sometimes fail dramatically. So TD is more risky. In practice, TD is easier to implement by sampling and has lower variance. BE needs more samples to eliminate the effect of noise.

3. Krylov subspace methods

The earliest Krylov subspace methods were developed by Lanczos, Hestenes and Stiefel in the late 1940s and early 1950s for solving large sparse linear systems of equations $Ax=b$ (with nonsingular $A\in {\mathbb{R}}^{n\times n}$). Krylov subspace methods form a the class of iterative method. An iterative method that solves the linear systems by producing a sequence of iterates $\left\{x^{\left(k\right)}\right\}$ starting from an initial point $x^{\left(0\right)}$, which is expected to converge to the true solution. Iterative methods have advantages over the direct methods (e.g. matrix factorization methods) when A is a large sparse matrix whose structure is hard to exploit. Moreover, we can stop an iterative method early if we just need an approximate solution, but we cannot do this with a direct method. Krylov subspace methods are different from classical iterative methods (e.g. matrix splitting methods) because they have memory, i.e. they construct $x^{(k+1)}$ not only from $x^{\left(k\right)}$ but also from the information of the previous iterates. So they are faster. Furthermore, Krylov subspace methods require only the ability to determine $Ax$ for any x (some also require the computation of $A^{T}x$).

3.1. Projection Process

First we will introduce a general projection process to find the solution $\widehat{x}$. Let $x_0$ denote an initial approximation to $\widehat{x}$ and then the initial residual $r_0\equiv b-A{x}_0$. At step k, we impose some conditions on $x_k$:

$$x_k\in x_0+{\mathcal{S}}_k,r_k\equiv b-A{x}_k\perp {\mathcal{C}}_k,k=1,2,\cdots .$$

(30)

where ${\mathcal{S}}_k\subseteq {\mathbb{R}}^n$ is a k-dimensional subspace, referred to as the search subspace, and ${\mathcal{C}}_k\subseteq {\mathbb{R}}^n$ is also a $k-$dimensional subspace, referred to as the constraint subspace. Since $S_k$ has dimension k, we have k degrees of freedom to construct $x_k$. In order to determine $x_k$, we need ${\mathcal{C}}_k$ to give us k constraints which are imposed on the residual $r_k$ in (30).

Let $V_k\in {\mathbb{R}}^{n\times k}$ and $U_k\in {\mathbb{R}}^{n\times k}$ be two full-column rank matrices such that

$${\mathcal{S}}_k=range\left(V_k\right),{\mathcal{C}}_k=range\left(U_k\right)$$

(31)

Then from (1) there exists $z_k\in {\mathbb{R}}^k$ such that

$$x_k=x_0+V_k{z}_k$$

(32)

and

$$U_k^T(b-A{x}_0-A{V}_k{z}_k)=0$$

(33)

which can be rewritten as

$$U_k^{T}A{V}_k{z}_k=U_k^T{r}_0$$

(34)

Theorem 1. 

1(see Theorem 2.1.2 in [8]) The matrix $U_k^{T}A{V}_k$ in Equation (34) is nonsingular if and only if

$${\mathbb{R}}^n=A{\mathcal{S}}_k\oplus {\mathcal{C}}_k^{\perp }$$

(35)

Now we assume that (35) holds. The linear system (34) can be regarded as a projected one of the original linear system $Ax=b$. As the dimensions of the matrix $U_k^{T}A{V}_k$ is $k\times k$, (34) can also be regarded as a reduced linear system.

From (32) and (34)

$$x_k=x_0+V_k{\left(U_k^{T}A{V}_k\right)}^{-1}{U}_k^T{r}_0$$

(36)

Then

$$\begin{array}{cc}\hfill r_k& =b-A{x}_k\hfill \\ \hfill & =b-A{x}_0-A{V}_k{\left(U_k^{T}A{V}_k\right)}^{-1}{U}_k^T{r}_0=(I-P_k)r_0\hfill \end{array}$$

where

$$P_k\equiv A{V}_k{\left(U_k^{T}A{V}_k\right)}^{-1}{U}_k^T$$

(37)

The matrix $P_k$ is an oblique projector because $P_k^2=P_k$. So $r_k$ is the projection of $r_0$ onto ${\mathcal{C}}_k^{\perp }$. And it is easy to verify that in (35) $A{\mathcal{S}}_k=\mathcal{R}\left(P_k\right)$ and ${\mathcal{C}}_k^{\perp }=\mathcal{N}\left(P_k\right)=\mathcal{R}(I-P_k).$

Theorem 2. 

2 (see Lemma 2.1.7 in [8]) Suppose that (34) holds. If $r_0\in {\mathcal{S}}_k$ and $A{\mathcal{S}}_k={\mathcal{S}}_k$, then $x_k$ is the solution of $Ax=b$

Theorem 2 gives a sufficent condition for $x_k$ to be the solution of $Ax=b$ at step k, that is in order to to make $r_k=0$ we need to have $r_0=P_k{r}_0$[8]. It seems like a good idea to start the projection process with the search space ${\mathcal{S}}_1=span\left\{r_0\right\}$, and to build up a nested sequence of search spaces.

$$span\left(r_0\right)={\mathcal{S}}_1\subseteq {\mathcal{S}}_2\subseteq {\mathcal{S}}_3\subseteq \cdots$$

(38)

and these spaces eventually satisfy $A{\mathcal{S}}_k={\mathcal{S}}_k$ for some k. Then we can obtain the solution in finite number of steps and we hope $k\ll n$.

3.2. Krylov Subspaces

Given $A\in {\mathbb{R}}^{n\times n}$ and a nonzero vector $v\in {\mathbb{R}}^n$, consider the following Krylov sequence,

$$v,Av,A^{2}v,\cdots .$$

(39)

There exists an integer $d=d(A,v)$ such that $\{v,\cdots A^{d-1}v\}$ are linearly independent, but $\{v,\cdots A^{d}v\}$ are linearly dependent. It’s easy to see that $d<n$. The matrix ${\mathcal{K}}_k(A,v)=span[v,Av,\cdots A^{k-1}v]$ is called a Krylov subspace. We will show that $d\le m$, where m is the degree of the minimal polynomial of A to be introduced below.

The characteristic polynomial of matrix A is defined by

$$P_A\left(t\right)=det(tI-A)=\prod _{i=1}^k{(t-{\lambda }_i)}^{s_i}$$

(40)

where ${\lambda }_1,\cdots ,{\lambda }_k$ are distinct eigenvalues of A and $s_i$ is the algebraic multiplicity for the corresponding ${\lambda }_i$. The minimal polynomial of A is defined as

$$m_A\left(t\right)=\prod _{i=1}^k{(t-{\lambda }_i)}^{m_i}$$

(41)

where $m_i$ is the largest size of Jordan block associated with ${\lambda }_i$. We all know that for a Jordan block $J_i$ with eigenvalue ${\lambda }_i$ whose size is $m_i$, ${(J-{\lambda }_{i}I)}^{m_i}=0$ and ${(J-{\lambda }_{i}I)}^{m_i-1}\ne 0$. And for a polynomial $q\left(t\right)$ and a matrix B whose Jordan blocks are $B_1,\cdots ,B_k$, $q\left(B\right)$ has Jordan blocks $q\left(B_1\right),\cdots ,q\left(B_k\right)$. Then for matrix A whose minimal polynomial is $m_A\left(t\right)$, we have $m_A\left(A\right)=0$. Then for all v, we have $m_A\left(A\right)v=0$. $m_A\left(t\right)$ could be written as

$$m_A\left(t\right)=\prod _{i=1}^k{(t-{\lambda }_i)}^{m_i}=\sum _{j=0}^m{\alpha }_j{t}^j$$

(42)

where $m=\sum _{i=1}^k{m}_i$. It it easy to see that ${\alpha }_m=1$ from (41), so $m_A\left(A\right)v=0$ means $A^{m}v$ can be represented by $v,Av,\cdots ,A^{m-1}v$ which means $A^{m}v\in {\mathcal{K}}_m(A,v)$ and we must have $d\le m.$ Combining this conclusion with theorem 2 and the fact that $r_0\in {\mathcal{K}}_i(A,r_0)$, it tells us that ${\mathcal{K}}_i(A,r_0)$ will eventually stop growing and become invariant under A and ${\mathcal{K}}_i(A,r_0),i=1,\cdots ,m$ is a choice of nested sequence of spaces satisfy (38). So we have

Theorem 3. 

3 If the minimal polynomial of nonsingular matrix A has degree m, then the solution to $Ax=b$ lies in the space ${\mathcal{K}}_m(A,b)$.

Actually, with different choice of ${\mathcal{S}}_k$ and ${\mathcal{C}}_k$, we have different Krylov subspace methods as Table 1 shows.

Table 1. Different Krylov subspace methods.

Name	Problem	${\mathcal{S}}_{\mathit{k}}$	${\mathcal{C}}_{\mathit{k}}$	Property
CG	SPD${}^1$ A	${\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{x-x_{*}}_A$
SYMMLQ	$A^T=A$	$A{\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{x-x_{*}}_2$${}^2$
MINRES	$A^T=A$	${\mathcal{K}}_{k(A,r_0)}$	$A{\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$
GMRES	general A	${\mathcal{K}}_{k(A,r_0)}$	$A{\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$
BiCG	general A	${\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A^T,{\tilde{r}}_0)}$	Does not solve an optimization problem
QMR	general A	${\mathcal{K}}_{k(A,r_0)}$	no ${\mathcal{C}}_k$	$x_k\approx \arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$

${}^1$ Symmetric positive definite A. $\forall x\ne 0,x^{T}Ax>0$. ${}^2$ Vector 2-norm. $x_2=\sqrt{x^{T}x}$.

Table 1. Different Krylov subspace methods.

Name	Problem	${\mathcal{S}}_{\mathit{k}}$	${\mathcal{C}}_{\mathit{k}}$	Property
CG	SPD${}^1$ A	${\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{x-x_{*}}_A$
SYMMLQ	$A^T=A$	$A{\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{x-x_{*}}_2$${}^2$
MINRES	$A^T=A$	${\mathcal{K}}_{k(A,r_0)}$	$A{\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$
GMRES	general A	${\mathcal{K}}_{k(A,r_0)}$	$A{\mathcal{K}}_{k(A,r_0)}$	$x_k=\arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$
BiCG	general A	${\mathcal{K}}_{k(A,r_0)}$	${\mathcal{K}}_{k(A^T,{\tilde{r}}_0)}$	Does not solve an optimization problem
QMR	general A	${\mathcal{K}}_{k(A,r_0)}$	no ${\mathcal{C}}_k$	$x_k\approx \arg {\min }_{x_0+{\mathcal{K}}_{k(A,r_0)}}{b-Ax}_2$

${}^1$ Symmetric positive definite A. $\forall x\ne 0,x^{T}Ax>0$. ${}^2$ Vector 2-norm. $x_2=\sqrt{x^{T}x}$.

3.3. The Arnoldi algorithm for orthonormal basis

If we form the matrix $K_j(A,r_0)$,

$$K_j(A,r_0)=[r_0,A{r}_0,A^2{r}_0,\cdots ,A^{j-1}{r}_0]$$

(43)

because $A^j{r}_0$ converges to an eigenvector corresponding to the largest eigenvalue in magnitude, it is likely that $K_j(A,r_0)$ is very ill-conditioned for large j. In order to avoid numerical troubles, we want to find orthonormal vectors ${\left\{v_i\right\}}_{i=1}^j$ such that

$${\mathcal{K}}_j(A,r_0)=range\left(K_j(A,r_0)\right)=span\{v_1,v_2,\cdots ,v_j\}$$

(44)

The idea is to use the Gram–Schmidt orthogonalization process and we show one step of it in the following.

Suppose we want to find $v_1$ and $v_2$ such that

$${\mathcal{K}}_2(A,r_0)=span\{v_1,v_2\},v_i^T{v}_j={\delta }_{ij},i,j=1,2$$

(45)

Then this can be done by:

1. Normalize $v_1$: $v_1:=r_0/{r_0}_2$, then ${v_1}_2=1.$

2. Let $w_1:=A{v}_1$, then

$${\mathcal{K}}_2(A,r_0)=span\{v_1,w_1\}$$

(46)

3. Orthogonalize $w_1$ against $v_1$: $v_2:=w_1-\left(w_1^T{v}_1\right)v_1$, then

$$v_2^T{v}_1=w_1^T{v}_1-\left(v_1^T{v}_1\right)\left(w_1^T{v}_1\right)=0$$

(47)

and we still have

$${\mathcal{K}}_2(A,r_0)=span\{v_1,v_2\}$$

(48)

4. Normalize $v_2$: $v_2:=v_2/{v_2}_2$, then ${v_2}_2=1$ and we still have ${\mathcal{K}}_2(A,r_0)=span\{v_1,v_2\}$.

We can write Arnoldi algorithm in terms of matrix. Let $V_k=[v_1,\cdots ,v_k]$ whose columns are composed by the orthonormal vectors $v_1,\cdots ,v_k$. After k steps, the algorithm has [8]

$$A{V}_k=V_k{H}_k+h_{k+1,k}{v}_{k+1}{e}_k^T,k=1,2,\cdots$$

(49)

$$H_k=\left[\begin{array}{ccccc}{h}_{11}& h_{12}& \cdots & h_{1,k-1}& h_{1,k}\\ h_{21}& h_{22}& \cdots & h_{2,k-1}& h_{2,k}\\ & h_{32}& \cdots & h_{3,k-1}& h_{3,k}\\ & & \ddots & & ⋮\\ & & & h_{k,k-1}& h_{k,k}\end{array}\right]$$

where $H_k$ is the $k\times k$ unreduced upper Hessenberg matrix with $h_{ij}=w_j^T{v}_i,i=1,\cdots ,j$ and $h_{j+1,j}={w_j}_2$ and $e_k$ is the k-th column of $I_k$. At the m-th iteration, we get $h_{m+1,m}=0$. Then Arnoldi algorithm stops and we have

$$A{V}_m=V_m{H}_m,V_m^T{V}_m=I_m$$

(50)

4. Connections

One direct and most common application of the Krylov subspace methods onto VFA is to construct basis functions [11]. Recall as (5) shows, we need to solve linear system

$$(I-\gamma P)v=r,0\le \gamma \le 1$$

(51)

Suppose $I-\gamma P\in {\mathbb{R}}^{n\times n}$ is nonsingular, since P is a stochastic matrix and $0\le \gamma <1$, we have $\sigma \left(\gamma P\right)<1$. If we want to approximate ${(I-\gamma P)}^{-1}$, one way is to use the Neumann series expansion [13],

$${(I-\gamma P)}^{-1}=\sum _{i=0}^{\infty }{\left(\gamma P\right)}^i$$

(52)

Then the solution of (5) lies in ${\mathcal{K}}_m(P,r)$, where m is the degree of the minimal polynomial of P, as what shows below,

$$v={(I-\gamma P)}^{-1}r=\sum _{i=0}^{\infty }{\left(\gamma P\right)}^{i}r\in span\{r,Pr,P^{2}r,\cdots \}=span\{r,Pr,P^{2}r,\cdots ,P^{m-1}r\}$$

(53)

Therefore, we can set our feature matrix as $\Phi =K_m(P,r)$.

In practice, we could first start with a random initial guess ${\theta }_0$ and $K_1(P,r)$ and progressively add $P^{i}r$ into the basis function matrix.

Parr [10] shows that the basis that constructed by Krylov subspace is equivalent the Bellman Error Basis Functions (BEBFs) and Model Error Basis Functions (MEBFs) which respectively add the Bellman error vectors and model error vectors progressively into the current basis.

Senda [15] applied the Conjugate Gradient (CG) algorithm, the Bi-Conjugate Gradient (BiCG) algorithm and the Block Product-type Krylov Subspace Method (BPKSM) for policy evaluation and did a comparison with the existing stationary iterative methods. What he concluded was that the algorithms based on the Krylov subspace methods were tens to hundreds times more efficient than conventional algorithms. He found the convergence rate of Krylov subspace algorithms improved as the iteration number increased, whereas the other algorithms remained the same. Algorithms based on Krylov subspace methods are far more efficient than the classical methods. The Krylov subspace methods can also be used in other machine learning area [19,20,21,22,23,24].

5. Conclusion and Discussion

From the last section we can see that, constructing the basis functions is only a simple and direct usage of Krylov subspace methods. It works better on model-based learning, in which we have enough information of the transition function P. But whether or not it can be extended to model-free learning has not been well studied. Bertsekas [2] provided a good example on how to connect proximal iterations for solving $x=Ax+b$ with TD($\lambda$) methods and got some interesting results. We can get some inspiration from it and try to use Krylov subspace methods for linear systems to solve approximate DP problems as well.