October 2011 (Revised Feb. 2012) Report LIDS - 2874

Lambda-Policy Iteration: A Review and a

New Implementation††† To appear in Reinforcement Learning and Approximate
Dynamic Programming for Feedback Control, by F. Lewis and D. Liu (eds.), IEEE
Press Computational Intelligence Series.

Dimitri P. Bertsekas ††‡ The author is with the Dept. of Electr. Engineering
and Comp. Science, M.I.T., Cambridge, Mass., 02139. His research was supported
by NSF Grant ECCS-0801549, and by the Air Force Grant FA9550-10-1-0412. Thanks
are due to Bruno Scherrer for helpful comments, to Huizhen Yu for related
collaboration, and to Mengdi Wang for assistance with computational
experimentation.

Abstract

In this paper we discuss $\lambda$-policy iteration, a method for exact and
approximate dynamic programming. It is intermediate between the classical
value iteration (VI) and policy iteration (PI) methods, and it is closely
related to optimistic (also known as modified) PI, whereby each policy
evaluation is done approximately, using a finite number of VI. We review the
theory of the method and associated questions of bias and exploration arising
in simulation-based cost function approximation. We then discuss various
implementations, which offer advantages over well-established PI methods that
use LSPE($\lambda$), LSTD($\lambda$), or TD($\lambda$) for policy evaluation
with cost function approximation. One of these implementations is based on a
new simulation scheme, called geometric sampling, which uses multiple short
trajectories rather than a single infinitely long trajectory.

1\. INTRODUCTION

Approximate dynamic programming (DP for short) has attracted substantial
research interest, and has a wide range of applications, because of its
potential to address large and complex problems that may not be treatable in
other ways. The literature on the subject is very extensive, and includes
several textbooks, research monographs, and surveys that relate to the
computational context of this paper. For a nonexhaustive list, we mention the
books by Bertsekas and Tsitsiklis [BeT96], Sutton and Barto [SuB98], Gosavi
[Gos03], Cao [Cao07], Chang, Fu, Hu, and Marcus [CFH07], Meyn [Mey07], Powell
[Pow07], Borkar [Bor08], Haykin [Hay08], Busoniu, Babuska, De Schutter, and
Ernst [BBD10], and the author’s text in preparation [Ber11a]; the edited
volumes and special issues by White and Sofge [WhS92], Si, Barto, Powell, and
Wunsch [SBP04], Lewis, Lendaris, and Liu [LLL08], and the 2007-2009
Proceedings of the IEEE Symposium on Approximate Dynamic Programming and
Reinforcement Learning; and the recent surveys by Borkar [Bor09], Lewis and
Vrabie [LeV09], Werbos [Wer09], Szepesvari [Sze10], and Bertsekas [Ber11b].

The purpose of this paper is to critically review and extend a class of
methods for exact and approximate DP, which are based on the $\lambda$-policy
iteration ($\lambda$-PI) method, proposed by Bertsekas and Ioffe [BeI96]. This
method is intermediate between the classical value iteration (VI) and policy
iteration (PI) methods, and it is closely related to optimistic (also known as
modified) PI, whereby each policy evaluation is done approximately, using a
finite number of VI. It was originally used as the starting point for the
development of approximate simulation-based DP methods of the temporal
difference (TD) type, such as LSPE($\lambda$) (see [BeI96], and also [BeT96],
Sections 2.3.1 and 8.3). The emphasis in this paper is on implementations of
$\lambda$-PI, which provide alternatives to approximate PI methods that use
other more established methods for policy evaluation.

We will focus on the $\alpha$-discounted $n$-state Markovian Decision Problem
(MDP), although the main ideas are more broadly applicable. The problem
involves states $1,\ldots,n$, controls $u\in U(i)$ at state $i$, transition
probabilities $p_{ij}(u)$, and cost $g(i,u,j)$ for transition from $i$ to $j$
under control $u$. A (stationary) policy $\mu$ is a function from states $i$
to admissible controls $u\in U(i)$, and $J_{\mu}(i)$ is the cost starting from
state $i$ and using policy $\mu$. It is well-known (see e.g., Puterman [Put94]
or Bertsekas [Ber07]) that the vector $J_{\mu}\in\Re^{n}$, which has
components $J_{\mu}(i)$,††† In our notation, $\Re^{n}$ is the $n$-dimensional
Euclidean space, all vectors in $\Re^{n}$ are viewed as column vectors, and a
prime denotes transposition. The identity matrix is denoted by $I$. is the
unique fixed point of the mapping $T_{\mu}:\Re^{n}\mapsto\Re^{n}$, which maps
$J\in\Re^{n}$ to the vector $T_{\mu}J\in\Re^{n}$ that has components

$(T_{\mu}J)(i)=\sum_{j=1}^{n}p_{ij}\big{(}\mu(i)\big{)}\big{(}g(i,\mu(i),j)+\alpha
J(j)\big{)},\qquad i=1,\ldots,n.$ (1.1)

Similarly, the optimal costs starting from $i=1,\ldots,n$, are denoted
$J^{*}(i)$, and the optimal cost vector $J^{*}\in\Re^{n}$, which has
components $J^{*}(i)$, is the unique fixed point of the mapping
$T:\Re^{n}\mapsto\Re^{n}$ defined by

$(TJ)(i)=\min_{u\in U(i)}\sum_{j=1}^{n}p_{ij}(u)\big{(}g(i,u,j)+\alpha
J(j)\big{)},\qquad i=1,\ldots,n.$ (1.2)

An important property is that $T_{\mu}$ and $T$ are sup-norm contractions. In
particular, the iterations $J_{k+1}=T_{\mu}J_{k}$ and $J_{k+1}=TJ_{k}$
converge to $J_{\mu}$ and $J^{*}$, respectively, from any starting point
$J_{0}$ \- this is the VI method.

A major alternative to VI is PI. It produces a sequence of policies and
associated cost functions through iterations that have two phases: policy
evaluation (where the cost function of a policy is evaluated), and policy
improvement (where a new policy is generated). In the exact form of the
algorithm, the current policy $\mu$ is improved by finding $\bar{\mu}$ that
satisfies $T_{\bar{\mu}}J_{\mu}=TJ_{\mu}$ [i.e., by minimizing in the right-
hand side of Eq. (1.2) with $J_{\mu}$ in place of $J$]. The improved policy
$\bar{\mu}$ is evaluated by solving the linear system of equations
$J_{\bar{\mu}}=T_{\bar{\mu}}J_{\bar{\mu}}$, and $(J_{\bar{\mu}},\bar{\mu})$
becomes the new cost vector-policy pair, which is used to start a new
iteration. Thus, the exact form of PI can be succinctly defined as

$T_{\mu_{k+1}}J_{k}=TJ_{k},\qquad J_{k+1}=T_{\mu_{k+1}}J_{k+1},$ (1.3)

with the equation on the left describing the policy improvement and the
equation on the right describing the evaluation of ${\mu_{k+1}}$.

In a variant of the method, a policy $\mu_{k+1}$ is evaluated by a finite
number of applications of $T_{\mu_{k+1}}$ to an approximate evaluation of the
preceding policy. This is known as “optimistic” or “modified” PI, and its
motivation is that in problems with a large number of states, the linear
system $J_{k+1}=T_{\mu_{k+1}}J_{k+1}$ cannot be practically solved directly by
matrix inversion, so it is best solved iteratively by VI. The method can be
succinctly defined as

$T_{{\mu_{k+1}}}J_{k}=TJ_{k},\qquad J_{k+1}=T_{\mu_{k+1}}^{m_{k}}J_{k}.$ (1.4)

If the number $m_{k}$ of applications of $T_{{\mu_{k+1}}}$ is very large, the
exact form of PI is essentially obtained, but practice has shown that it is
most efficient to use a moderate value of $m_{k}$. In this case, the algorithm
looks like a hybrid of VI and PI, involving a sequence of alternate
applications of $T$ and $T_{\mu_{k}}$, with $\mu_{k}$ changing over time.
Optimistic PI is generally believed to be more computationally efficient that
either VI or PI. This is particularly so for problems where $n$ is very large
and implementation of exact PI is difficult due to the associated $n\times n$
matrix inversion, and also for problems with a large number of controls, where
the overhead due to minimization over all controls $u\in U(i)$ in the mapping
$T$ [cf. Eq. (1.2)] is substantial.

We note that the convergence properties of the optimistic PI method (1.4) are
quite complicated and have been the subject of continuing research. The
convergence $J_{k}\to J^{*}$ has been established by Rothblum [Rot79] (see
also the more recent work by Canbolat and Rothblum [CaR11], which extends some
of the results of [Rot79]). On the other hand, when optimistic PI is
implemented asynchronously (as it normally would be when simulation is used),
it may oscillate as shown by the convergence counterexamples of Williams and
Baird [WiB93]. Recent work of Bertsekas and Yu [BeY10a], [BeY10b], [YuB11] has
developed convergent variants of synchronous and asynchronous optimistic PI
and Q-learning, based on a new way to perform policy evaluation: by solving
approximately an optimal stopping problem rather than a system of linear
equations.

The $\lambda$-PI method is a form of optimistic PI, given by

$T_{\mu_{k+1}}J_{k}=TJ_{k},\qquad J_{k+1}=T_{\mu_{k+1}}^{(\lambda)}J_{k},$
(1.5)

where for any $\mu$ and $\lambda\in[0,1)$, $T_{\mu}^{(\lambda)}$ is the linear
mapping given by

$T_{\mu}^{(\lambda)}=(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T_{\mu}^{\ell+1}.$
(1.6)

Note that the mapping $T_{\mu}^{(\lambda)}$ is central in much recent research
on approximate DP, simulation-based PI, and TD methods, as will be discussed
in the sequel.

To compare the optimistic PI method (1.4) and the $\lambda$-PI method (1.5),
note that both mappings $T_{\mu_{k+1}}^{m_{k}}$ and
$T_{\mu_{k+1}}^{(\lambda)}$ appearing in Eqs. (1.4) and (1.5), involve
multiple applications of the VI mapping $T_{\mu_{k+1}}$: a fixed number
$m_{k}$ in the former case (with $m_{k}=1$ corresponding to VI and
$m_{k}\to\infty$ corresponding to PI), and a geometrically weighted number in
the latter case (with $\lambda=0$ corresponding to VI and $\lambda\to 1$
corresponding to PI). Thus optimistic PI and $\lambda$-PI are similar: they
just control the accuracy of the approximation $J_{k+1}\approx J_{\mu_{k+1}}$
by applying VI in different ways. In a classical DP/non-simulation-based
setting, $\lambda$-PI is far more complicated relative to optimistic PI, since
exact computations using the mapping $T_{\mu}^{(\lambda)}$ are unwieldy.
However, this advantage of optimistic PI is dissipated in a simulation
context, where computations involving $T_{\mu}^{(\lambda)}$ can be performed
conveniently, as extensive analytical and experimental work with TD methods
has demonstrated.

Recent research on DP has focused on the use of simulation, in order to deal
with model-free situations where the transition probabilities and/or the cost
per stage are not known explicitly, and also to deal with the associated high-
dimensional linear algebra operations. For problems with very large number of
states, the evaluation of various fixed points of mappings, such as $T_{\mu}$
or $T_{\mu}^{(\lambda)}$, is typically done by approximation with a vector
$\Phi r$ from the subspace $S=\\{\Phi r\mid r\in\Re^{s}\\}$ that is spanned by
the columns of an $n\times s$ matrix $\Phi$. In this paper we will focus on
the projected equation approach, whereby given a generic mapping
$L:\Re^{n}\mapsto\Re^{n}$ (such as for example $T_{\mu}$) we approximate its
fixed point by solving the equation

$\Phi r=\Pi L(\Phi r),$

where $\Pi$ denotes projection onto the subspace $S$. The projection is with
respect to a Euclidean norm $\|\cdot\|_{\xi}$, weighted by a suitable vector
$\xi$ of positive weights. An alternative possibility is to solve instead the
equation

$\Phi r=\Pi L^{(\nu)}(\Phi r),$ (1.7)

where, similar to Eq. (1.6),

$L^{(\nu)}=(1-\nu)\sum_{\ell=0}^{\infty}\nu^{\ell}L^{\ell+1},$

and $\nu\in[0,1)$ is a parameter [not necessarily the same as the $\lambda$
parameter in Eqs. (1.5)-(1.6)]. In our context we will encounter several
different types of mappings $L$, and in all cases $L$ is a contraction with
respect to the projection norm $\|\cdot\|_{\xi}$, with fixed point $\hat{J}$,
while $\Pi L^{(\nu)}$ are contractions with respect to $\|\cdot\|_{\xi}$ for
all $\nu\in[0,1)$. It is well-known that the fixed point of $\Pi L^{(\nu)}$,
denoted $\Phi r(\nu)$, converges to $\Pi\hat{J}$ as $\nu\to 1$. The norm of
the difference $\Phi r(\nu)-\Pi\hat{J}$ is known as the bias. Its size/norm
depends on $\nu$ and is generally smaller as $\nu$ gets closer to 1 (see
[BeT96], [TsV97], [YuB10] for error bound analyses).

A common example of fixed point approximation in PI is when $L=T_{\mu}$ for a
policy $\mu$, in which case the fixed point of $\Pi L$ or $\Pi L^{(\nu)}$ is
an approximation to the fixed point of $T_{\mu}$, i.e., the cost vector
$J_{\mu}$. If the Markov chain corresponding to $\mu$ is irreducible and $\xi$
is the corresponding steady-state distribution vector, the mapping $\Pi
T^{(\lambda)}_{\mu}$ is a contraction with respect to $\|\cdot\|_{\xi}$ for
all $\lambda\in[0,1)$, and is unique fixed point, denoted $\Phi
r_{\mu}(\lambda)$, converges to $\Pi J_{\mu}$ as $\lambda\to 1$. Generally,
the projected equation $\Phi r=\Pi T^{(\lambda)}_{\mu}(\Phi r)$ is solved by a
simulation process that generates a sequence of states according to a sampling
scheme to be discussed later, and then by matrix inversion [this is the Least
Squares Temporal Differences [LSTD($\lambda$)] method, proposed by Bradtke and
Barto [BrB96]], or by iteration, using the TD($\lambda$) method, proposed by
Sutton [Sut87] and analyzed by Tsitsiklis and VanRoy [TsV97] among others, or
the Least Squares Policy Evaluation [LSPE($\lambda$)] method, proposed by
Bertsekas and Ioffe [BeI96].††† The paper [BeI96] as well as the book [BeT96]
used the name “$\lambda$-policy iteration” for both the lookup table and the
compact representation versions of the method described here, and tested a
compact representation version on the game of tetris, a challenging SSP
problem. The name “LSPE” was first used in the subsequent paper by Nedić and
Bertsekas [NeB03] to describe a specific iterative implementation of the
$\lambda$-PI method with cost function approximation for discounted MDP
(essentially the discounted version of the implementation used in [BeI96] and
[BeT96] for the aforementioned tetris case study). Reference [NeB03] proved
convergence of the LSPE($\lambda$) method, as described in Section 3.1, for
the case of a diminishing stepsize. Convergence for a stepsize equal to 1 was
proved shortly afterwards by Bertsekas, Borkar, and Nedić [BBN04]. The use of
two different names for essentially the same method has been a source of some
confusion. While in practical implementations these two names refer to
algorithms that are closely related, we reserve the name “$\lambda$-policy
iteration” for the more abstract form (1.5)-(1.6), and we will view
LSPE($\lambda$) as an implementation of $\lambda$-PI (see Section 4.1). These
methods are extensively discussed in the literature, and exhibit complex and
sometimes pathological behavior, particularly when embedded within PI (see
[Ber95], [SzL06], [ThS09] for some notable failures, and [Ber10] for a recent
assessment). Moreover matrix inversion and iterative methods, like
TD($\lambda$), LSTD($\lambda$), and LSPE($\lambda$), can be used for solving
not only the projected equation $\Phi r=\Pi T^{(\lambda)}_{\mu}(\Phi r)$, but
also the more general equation $\Phi r=\Pi L^{(\nu)}(\Phi r)$ of Eq. (1.7), as
long as $L$ is a linear mapping that is convenient for the use of simulation
[and in the case of TD($\lambda$) and LSPE($\lambda$), $\Pi L^{(\nu)}$ is a
contraction; see [BeY09] or [Ber11c]].

In this paper we will review some of the basic issues in approximate PI using
the projected equation approach, thereby setting the stage for assessing the
relative strengths and weaknesses of the $\lambda$-PI methodology. We will
then focus on three alternative implementations of $\lambda$-PI, which involve
simulation and cost function approximation. The first is basically the
LSPE($\lambda$) method as implemented in [BeI96]. The second is an interesting
recent proposal by Thiery and Scherrer [ThS10a], who gave extensive and quite
successful computational results, as well as error bounds [ThS10b]. The third
implementation is new and may have some advantages over the first two. We will
argue that it deals better with the combined issues of bias and exploration.
This implementation embodies a new idea for $\lambda$-methods: a simulation
scheme, called geometric sampling, that uses multiple short trajectories with
random geometrically distributed length, and exploration-enhanced restart,
rather than a single infinitely long trajectory.

The three implementations are described in Section 4, following a discussion
of the generic properties of exact $\lambda$-PI in Section 2, and the
LSTD($\lambda$) and LSPE($\lambda$) methods in Section 3. In our description,
these implementations are model-based and use cost function approximation, but
there are versions that are model-free and use Q-factor approximation; these
can be straightforwardly constructed by the reader.

2\. LAMBDA-POLICY ITERATION WITHOUT COST FUNCTION APPROXIMATION

We first recall a central result from [BeI96]. It provides a helpful
characterization of the $\lambda$-PI method (1.5), which will later become the
basis for cost function approximations.

Proposition 2.1: Given $\lambda\in[0,1)$, $J_{k}$, and $\mu_{k+1}$, consider
the mapping $W_{k}$ defined by W_kJ=(1-λ)T_μ_k+1J_k+λT_μ_k+1J. (a) $W_{k}$ is
a sup-norm contraction of modulus $\lambda\alpha$. (b) The vector
$J_{k+1}=T_{\mu_{k+1}}^{(\lambda)}J_{k}$ generated next by the $\lambda$-PI
method (1.5) is the unique fixed point of $W_{k}$. (2.1)

Proof: (a) For any two vectors $J$ and $\bar{J}$, using the definition (2.1)
of $W_{k}$, we have

$\|W_{k}J-W_{k}\bar{J}\|=\bigl{\|}\lambda(T_{\mu_{k+1}}J-T_{\mu_{k+1}}\bar{J})\bigr{\|}=\lambda\|T_{\mu_{k+1}}J-T_{\mu_{k+1}}\bar{J}\|\leq\lambda\alpha\|J-\bar{J}\|,$

where $\|\cdot\|$ denotes the sup-norm, so $W_{k}$ is a sup-norm contraction
with modulus $\lambda\alpha$.

(b) We have

$J_{k+1}=T_{\mu_{k+1}}^{(\lambda)}J_{k}=(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T_{\mu_{k+1}}^{\ell+1}J_{k},$

so the fixed point property to be shown, $J_{k+1}=W_{k}J_{k+1}$, is written as

$(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T_{\mu_{k+1}}^{\ell+1}J_{k}=(1-\lambda)T_{\mu_{k+1}}J_{k}+\lambda
T_{\mu_{k+1}}(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T_{\mu_{k+1}}^{\ell+1}J_{k},$

and evidently holds. Q.E.D.

From part (b) of the preceding proposition, we see that
$J_{k+1}=W_{k}J_{k+1}$, or equivalently

$J_{k+1}(i)=\sum_{j=1}^{n}p_{ij}\big{(}\mu_{k+1}(i)\big{)}\Big{(}g\big{(}i,\mu_{k+1}(i),j\big{)}+(1-\lambda)\alpha
J_{k}(j)+\lambda\alpha J_{k+1}(j)\Big{)},\qquad i=1,\ldots,n.$ (2.2)

The solution of this fixed point equation can be obtained by viewing it as
Bellman’s equation for two equivalent MDP.

(a) As Bellman’s equation for an infinite-horizon $\lambda\alpha$-discounted
MDP where $\mu_{k+1}$ is the only policy, and the cost per stage is

$g\big{(}i,\mu_{k+1}(i),j\big{)}+(1-\lambda)\alpha J_{k}(j).$

(b) As Bellman’s equation for an infinite-horizon stopping problem where
$\mu_{k+1}$ is the only policy. In particular, $J_{k+1}$ is the cost vector of
policy $\mu_{k+1}$ in a stopping problem that is derived from the given
$\alpha$-discounted problem by introducing transitions from each state $j$ to
an artificial termination state as follows: at state $i$ we first make a
transition to $j$ with probability $p_{ij}\big{(}\mu_{k+1}(i)\big{)}$ and
transition cost $g\big{(}i,\mu_{k+1}(i),j\big{)}$; then we either stay at $j$
and wait for the next transition (this occurs with probability $\lambda$), or
else we move from $j$ to the termination state with an additional termination
cost $\alpha J_{k}(j)$ (this occurs with probability $1-\lambda$). All
transition costs as well as the termination cost are discounted by an
additional factor $\alpha$ with each transition.

The convergence and rate of convergence of the $\lambda$-PI method (1.5) was
given in [BeI96] and also in [BeT96], Prop. 2.8. We will simply quote the
results for completeness.

Proposition 2.2: Assume that $\lambda\in[0,1)$, and let $\\{J_{k},\mu_{k}\\}$
be the sequence generated by the $\lambda$-PI method (1.5). Then $J_{k}$
converges to $J^{*}$. Furthermore, for all $k$ greater than some index
$\bar{k}$, $\mu_{k}$ is optimal.

Proposition 2.3: Let the assumptions of Prop. 2.2 hold and let $\bar{k}$ be
the index such that for all $k\geq\bar{k}$, $\mu_{k}$ is optimal. The sequence
$\\{J_{k}\\}$ generated by the $\lambda$-PI method (1.5) satisfies for all
$k>\bar{k}$ ∥J_k+1-J^*∥≤α(1-λ)1-λα∥J_k-J^*∥, where $\|\cdot\|$ denotes the
sup-norm. (2.3)

Note that the convergence rate estimate (2.3) holds only for $k\geq\bar{k}$,
essentially after an optimal policy has been identified, as per Prop. 2.2.
Nonetheless, this rate estimate is qualitatively correct, and supports the
empirical observation that the iterates $(J_{k},\mu_{k})$ generated by
$\lambda$-PI converge faster as $\lambda$ increases. Indeed in the limit, as
$\lambda\to 1$, $\lambda$-PI becomes exact PI, and converges to the optimum in
a finite number of iterations. On the other hand, the computation of
$J_{k+1}=T_{\mu_{k+1}}^{(\lambda)}J_{k}$ [cf. Eq. (1.5)] becomes more time-
consuming as $\lambda$ increases, particularly when simulation is used,
because the simulation-based calculation of $T_{\mu_{k+1}}^{(\lambda)}J_{k}$
involves more simulation noise as $\lambda$ gets larger.

We finally note that Props. 2.2 and 2.3 apply to synchronous implementations
of $\lambda$-PI. When implemented asynchronously, $\lambda$-PI has similar
convergence difficulties to optimistic PI. To see this, note that asynchronous
implementations of these two methods essentially coincide when $m_{k}=1$ in
Eq. (1.5) and $\lambda=0$ in Eq. (1.4), and the counterexamples of Williams
and Baird [WiB93] apply. Thus the development of convergent versions of
asynchronous $\lambda$-PI is an open research question.

3\. APPROXIMATE POLICY EVALUATION USING PROJECTED EQUATIONS

In PI methods with cost function approximation, we evaluate $\mu$ by
approximating $J_{\mu}$ with a vector $\Phi r_{\mu}$ from the subspace
$S=\\{\Phi r\mid r\in\Re^{s}\\}$, spanned by the columns of an $n\times s$
matrix $\Phi$, which may be viewed as basis functions. We generate an
“improved” policy $\bar{\mu}$ using the formula $T_{\bar{\mu}}(\Phi
r_{\mu})=T(\Phi r_{\mu})$, i.e.,

$\bar{\mu}(i)\in\arg\min_{u\in
U(i)}\sum_{j=1}^{n}p_{ij}(u)\bigl{(}g(i,u,j)+\alpha\phi(j)^{\prime}r_{\mu}\bigr{)},\qquad
i=1,\ldots,n,$

where $\phi(j)^{\prime}$ is the row of $\Phi$ that corresponds to state $j$
[the method terminates with $\mu$ if $T_{\mu}(\Phi r_{\mu})=T(\Phi r_{\mu})$].
We then repeat with $\mu$ replaced by $\bar{\mu}$. For the purposes of this
paper, we assume that $\Phi$ has rank $s$, and that the Markov chain
corresponding to $\mu$ is irreducible.

As noted earlier, in the projected equation approach to approximate PI, we
approximate $J_{\mu}$ with a vector of the form $\Phi r_{\mu}(\lambda)$ that
solves the fixed point problem

$\Phi r=\Pi T_{\mu}^{(\lambda)}(\Phi r).$ (3.1)

Here $\Pi$ denotes projection onto the subspace $S$ with respect to a weighted
Euclidean norm $\|\cdot\|_{\xi}$, where $\xi=(\xi_{1},\ldots,\xi_{n})$ is a
probability distribution with positive components (i.e.,
$\|J\|_{\xi}^{2}=\sum_{i=1}^{n}\xi_{i}x_{i}^{2}$, where $\xi_{i}>0$ for all
$i$). In nonoptimistic PI methods, the projected equation (3.1) is solved
exactly, while in optimistic PI methods it is solved approximately. We note
that this approach has a long history in the context of Galerkin methods for
the approximate solution of high-dimensional or infinite-dimensional linear
equations (partial differential, integral, inverse problems, etc; see e.g.,
[Kra72], [Fle84]). In fact some of the policy evaluation theory referred to in
this paper applies to general projected equations arising in contexts beyond
DP (see [BeY09], [Ber09], [Yu10a,b], [Ber11c]). However, Monte Carlo
simulation is not part of the Galerkin methodology, as currently practiced in
the numerical analysis field. For this reason much of the extensive available
knowledge about Galerkin methods does not apply to the approximate DP context,
which is primarily simulation-oriented.

We now discuss some of the issues relating to projected equations. While we
focus on Eq. (3.1), much of our discussion also applies to the more general
projected equations.

Exploration-Contraction Tradeoff

An important choice in the projected equation approach is the distribution
$\xi$ that defines the projection norm $\|\cdot\|_{\xi}$. This distribution is
sometimes chosen to be the steady-state probability vector $\xi_{\mu}$ of the
Markov chain corresponding to $\mu$, in which case the mapping $\Pi
T_{\mu}^{(\lambda)}$ can be shown to be a contraction with respect to
$\|\cdot\|_{\xi_{\mu}}$ with modulus

$\alpha_{\lambda}={\alpha(1-\lambda)\over 1-\lambda\alpha}$ (3.2)

(see [BeT96], Lemma 6.6, or [Ber07], Prop. 6.3.3).

On the other hand the choice of $\xi$ is related to exploration, i.e., the
need to collect an adequately rich set of samples from a broad and
representative set of states. This is a critical issue in simulation-based PI,
and results in a well-known tradeoff: to evaluate a policy $\mu$, we may need
to generate cost samples using $\mu$, but this may affect the simulation
results by underrepresenting states that are unlikely to occur under $\mu$
(more weight is placed on states that are visited more frequently under
$\mu$). As a result, the cost-to-go estimates of the underrepresented states
may be highly inaccurate, causing potentially serious errors in the
calculation of the improved control policy.

A well-known approach for exploration is to choose $\xi$ to be a mixture of
the form

$\xi=(1-\beta)\xi_{\mu}+\beta\tilde{\xi},$ (3.3)

where $\beta\in(0,1)$ and $\tilde{\xi}$ is another distribution (often
referred to as the off-policy distribution), which is added to enhance
exploration (see the discussion of Section 1). Unfortunately, with such a
choice the contraction property of $\Pi T_{\mu}^{(\lambda)}$ comes into doubt:
it depends on the size of the parameters $\lambda$ and $\beta$ [it can be
shown that $\Pi T_{\mu}^{(\lambda)}$ is a contraction for any $\beta\in[0,1)$
provided $\lambda$ is close enough to 1, and it is a contraction for any
$\lambda\in[0,1)$ provided $\beta$ is close enough to 0]. This is important
because for convergence of iterative methods such as TD($\lambda$) and some
forms of LSPE($\lambda$), it is critical that $\Pi T_{\mu}^{(\lambda)}$ be a
contraction. Thus there is a tradeoff between exploration enhancement using
the mixture distribution (3.3) and ability to use a broader range of methods
for solution of the projected equation.

Bias

While the Bellman equation $J=T_{\mu}^{(\lambda)}J$ has the same fixed point
$J_{\mu}$ for all $\lambda\in[0,1)$, the fixed point $\Phi r_{\mu}(\lambda)$
of the projected version (3.1) depends on $\lambda$. The difference of $\Phi
r_{\mu}(\lambda)$ and the closest point of $S$ to $J_{\mu}$, $\Phi
r_{\mu}(\lambda)-\Pi J_{\mu}$, is generally nonzero. Its norm, the bias, tends
to decrease to 0 as $\lambda\uparrow 1$ and tends to increase as
$\lambda\downarrow 0$. It is known that the bias can be very large and may
seriously degrade the practical value of the approximate policy evaluation for
small values of $\lambda$; see [Ber95] for some examples.

The following is a well-known error bound for the case $\xi={\xi_{\mu}}$:

$\|J_{\mu}-\Phi
r_{\mu}(\lambda)\|_{\xi_{\mu}}\leq{1\over\sqrt{1-\alpha_{\lambda}^{2}}}\,\|J_{\mu}-\Pi
J_{\mu}\|_{\xi_{\mu}},$ (3.4)

where $\alpha_{\lambda}$ is given by Eq. (3.2), and $\|\cdot\|_{\xi_{\mu}}$ is
the weighted Euclidean norm corresponding to $\xi={\xi_{\mu}}$, the steady-
state probability vector of the Markov chain corresponding to $\mu$. Thus the
error bound becomes worse as $\lambda$ decreases (and $\alpha_{\lambda}$
increases), suggesting a larger size of bias. While the bound is rather
conservative, the paper by Yu and Bertsekas [YuB10] (see also Scherrer
[Sch10]) derives sharper error bounds, which also apply to cases where
$\xi\neq{\xi_{\mu}}$ and $\Pi T_{\mu}^{(\lambda)}$ is not a contraction. These
error bounds and the bound (3.4) are consistent in suggesting that the bias
increases as $\lambda$ decreases, and they are also largely consistent with
the results of computational experimentation.

Bias-Variance Tradeoff

In simulation-based methods for solving the projected equation (3.1), one must
deal with the effects of simulation error. Generally as $\lambda$ increases,
the methods become more vulnerable to simulation noise, and hence require more
sampling for good performance. Indeed, the noise in a simulation sample of an
$\ell$-stages cost vector $T^{\ell}_{\mu}J$ tends to be larger as $\ell$
increases, and from the formula

$T^{(\lambda)}_{\mu}=(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T^{\ell+1}_{\mu}$

it can be seen that simulation samples of $T^{(\lambda)}_{\mu}(\Phi r_{k})$
tend to contain more noise as $\lambda$ increases. This is consistent with
practical experience, and gives rise to the so called bias-variance tradeoff:
a large value of $\lambda$ to reduce bias results in slower and less reliable
computation because of higher simulation noise (and consequently, a larger
number of samples to achieve the same accuracy of various simulation-based
estimates). Generally, there is no rule of thumb for selecting $\lambda$,
which is usually chosen with some trial and error.

In summary, the preceding discussion suggests that if simulation noise is not
an issue (i.e., one can afford many simulation samples) one should choose
large values of $\lambda$, since then the bias is reduced and one may afford
greater exploration without losing the contraction property of $\Pi
T_{\mu}^{(\lambda)}$. In the contrary case, however, the degradation of the
estimate of $J_{\mu}$ due to simulation noise may offset whatever
bias/contraction benefits a large value of $\lambda$ may bring.

3.1  TD Methods

Most of the simulation-based methods for solving the projected equation use
explicitly or implicitly the notion of temporal difference (TD), which
originated in reinforcement learning with the works of Samuel [Sam59], [Sam67]
on a checkers-playing program. The first TD method is TD($\lambda$), which can
be viewed as an iterative stochastic approximation-type algorithm. The
LSTD($\lambda$) method is based on batch simulation: it first generates a
batch of state and cost samples, it approximates the projected equation $\Phi
r=\Pi T_{\mu}^{(\lambda)}(\Phi r)$ using these samples, and then solves the
equation directly by matrix inversion. Another TD method is LSPE($\lambda$),
which while being more iterative, shares much of the simulation philosophy of
LSTD($\lambda$).

To describe more specifically the LSTD($\lambda$) and LSPE($\lambda$) methods,
we first note that the orthogonality condition that characterizes the
projection in the projected equation $\Phi r=\Pi T_{\mu}^{(\lambda)}(\Phi r)$
is

$\Phi^{\prime}\Xi\big{(}\Phi r-T_{\mu}^{(\lambda)}(\Phi r)\big{)}=0,$ (3.5)

where $\Xi$ is the diagonal matrix with the vector $\xi$ along the diagonal
(see e.g., [Ber07]). Thus the projected equation (3.1) is equivalent to the
lower-dimensional equation (3.5), which can in turn be written in matrix form
as

$C^{(\lambda)}r=d^{(\lambda)},$ (3.6)

with

$C^{(\lambda)}=\Phi^{\prime}\Xi\big{(}I-P_{\mu}^{(\lambda)}\big{)}\Phi,\qquad
d^{(\lambda)}=\Phi^{\prime}\Xi g_{\mu}^{(\lambda)},$ (3.7)

and

$P_{\mu}^{(\lambda)}=(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}\alpha^{\ell+1}P_{\mu}^{\ell+1},\qquad
g_{\mu}^{(\lambda)}=\sum_{\ell=0}^{\infty}\lambda^{\ell}\alpha^{\ell}P_{\mu}^{\ell}g_{\mu},$
(3.8)

where $P_{\mu}$ and $g_{\mu}$ are the transition probability matrix and
expected single-stage cost vector corresponding to $\mu$. The LSTD($\lambda$)
and LSPE($\lambda$) methods use simulation-based approximations of
$C^{(\lambda)}$ and $d^{(\lambda)}$. This is done by simulating a state
sequence $(i_{0},\ldots,i_{t})$ and corresponding transition cost sequence,
using the current policy $\mu$ (perhaps with exploration enhancement, as
discussed earlier). Then after each simulated state $i_{\ell}$,
$\ell=0,\ldots,t$, is generated, estimates $C_{\ell}^{(\lambda)}$ and
$d_{\ell}^{(\lambda)}$ are obtained using the simulation samples up to time
$\ell$, using formulas that we will not give here, as they are not important
for our purposes. Such formulas, in various alternative forms, can be found in
several sources, including the textbooks cited earlier. The papers [NeB03],
[BeY09], [Yu10a], [Yu10b] discuss the conditions for the convergence
$\lim_{\ell\to\infty}C_{\ell}^{(\lambda)}=C^{(\lambda)}$,
$\lim_{\ell\to\infty}d_{\ell}^{(\lambda)}=d^{(\lambda)}$ to hold with
probability 1.

The LSTD($\lambda$) method is based on simple matrix inversion: after the last
state $i_{t}$ of the simulation trajectory is generated, it computes the
solution

$\hat{r}=\big{(}C_{t}^{(\lambda)}\big{)}^{-1}d_{t}^{(\lambda)}$ (3.9)

of the corresponding simulation-based approximation to Eq. (3.6),

$C_{t}^{(\lambda)}r=d_{t}^{(\lambda)},$ (3.10)

and approximates the cost vector $J_{\mu}$ by $\Phi\hat{r}$. An important
point is that $\hat{r}$ can be obtained regardless of whether $\Pi
T^{(\lambda)}_{\mu}$ is a contraction. It is only required that
$C_{t}^{(\lambda)}$ is invertible, a much less restrictive condition.

One version of the LSPE($\lambda$) method consists of iterative solution of
the system (3.10). It approximates the cost vector $J_{\mu}$ by $\Phi
r_{t+1}$, where $r_{t+1}$ is obtained at the last step of the iteration

$r_{\ell+1}=r_{\ell}-\gamma
G_{\ell}\bigl{(}C_{\ell}^{(\lambda)}r_{\ell}-d_{\ell}^{(\lambda)}\bigr{)},\qquad\ell=0,\ldots,t,$
(3.11)

where $r_{0}$ is some initial vector, likely the vector obtained from the
preceding policy evaluation, $\gamma$ is a positive stepsize, $G_{\ell}$ is
the matrix

$G_{\ell}=\left({1\over\ell}\sum_{m=0}^{\ell-1}\phi(i_{m})\phi(i_{m})^{\prime}\right)^{-1},\qquad\ell=0,\ldots,t,$
(3.12)

and as earlier, $\phi(i)^{\prime}$ denotes the $i$th row of the matrix $\Phi$.
In the original proposal of [BeI96] the stepsize is $\gamma=1$; convergence of
$\Phi r_{t}$ to the fixed point of $\Pi T^{(\lambda)}_{\mu}$ for this stepsize
was shown in [BBN04]. The matrix $G_{\ell}$ is a simulation-based
approximation of $(\Phi^{\prime}\Xi\Phi)^{-1}$ (alternative choices of
$G_{\ell}$ have been discussed recently in [Ber11b], [Ber11c]). There is also
an equivalent implementation of this iteration, which is based on solution of
a least squares problem (see Section 4.1).

The choice (3.12) for $G_{\ell}$ and the use of $\gamma=1$ are based on a view
of the method as an approximation to the projected value iteration method

$\Phi r_{\ell+1}=\Pi T_{\mu}^{(\lambda)}(\Phi r_{\ell}),$

which after some calculation can be written as

$\Phi
r_{\ell+1}=\Phi\big{(}r_{\ell}-(\Phi^{\prime}\Xi\Phi)^{-1}\big{(}C^{(\lambda)}r_{\ell}-d^{(\lambda)})\big{)},$

or equivalently, since $\Phi$ has full rank, as

$r_{\ell+1}=r_{\ell}-(\Phi^{\prime}\Xi\Phi)^{-1}\bigl{(}C^{(\lambda)}r_{\ell}-d^{(\lambda)}\bigr{)};$

cf. Eq. (3.11)-(3.12) with $\gamma=1$.

Note that the matrix inversion in Eq. (3.12) is not so onerous, because it can
be formed incrementally, with a rank-one correction as each sample becomes
available. On the other hand, contrary to LSTD($\lambda$) [and similar to
TD($\lambda$)], the LSPE($\lambda$) method (3.11)-(3.12) requires that $\Pi
T_{\mu}^{(\lambda)}$ be a contraction for convergence. Indeed if the
simulation is performed using the steady-state distribution $\xi_{\mu}$, it
can be shown that $\Pi T_{\mu}^{(\lambda)}$ is a contraction, but if the
simulation is performed using a mixture/off-policy distribution (3.3) for the
purpose of exploration-enhancement, the contraction property may be lost and
repeated iterations of the form (3.11) may diverge.

We finally note that in iteration (3.11) the underlying assumption is that we
update $r$ as simulation samples are collected and used to form ever improving
approximations to $C$ and $d$. An alternative is to use batch simulation, like
in LSTD: first simulate to obtain $C_{t}^{(\lambda)}$, $d_{t}^{(\lambda)}$,
and $G_{t}$, and then solve the system $C_{t}^{(\lambda)}r=d_{t}^{(\lambda)}$
iteratively rather than through the direct matrix inversion (3.9), by using
any number of iterations of the type (3.11). In fact, we may use only one
iteration, in which case the method takes the form

$r_{1}=r_{0}-\gamma
G_{t}\bigl{(}C_{t}^{(\lambda)}r_{0}-d_{t}^{(\lambda)}\bigr{)}.$ (3.13)

A single (or very few) iterations may be sufficient if $\lambda$ is close to
1, since then the contraction modulus of $\Pi T_{\mu}^{(\lambda)}$ is close to
0 (see e.g., [BeT96], Lemma 6.6, or [Ber07], Prop. 6.3.3), so a single
iteration with $\Pi T_{\mu}^{(\lambda)}$ is very effective, yielding a vector
that is close to its fixed point. We will return to this variant of the method
later.

3.2  Comparison of LSTD($\lambda$) and LSPE($\lambda$)

There has been speculation about the relative merits of LSTD($\lambda$) and
LSPE($\lambda$). Generally speaking, it is difficult to reach definitive
conclusions, as there are several complex factors to consider, such as the
length of the simulation sequence $(i_{0},\ldots,i_{t})$, and the potential
near-singularity of $C^{(\lambda)}$, which affects the error in the matrix
inversion in the LSTD($\lambda$) formula (3.9). As an illustration, consider a
few different situations:

(a) Assume, as an idealization, that an infinite number of samples is
collected. Then both methods yield in the limit the same result, the fixed
point of the projected equation $J=\Pi T_{\mu}^{(\lambda)}J.$ However, in
contrast to LSTD($\lambda$), in order to guarantee convergence,
LSPE($\lambda$) requires that $\Pi T_{\mu}^{(\lambda)}$ is a contraction,
which interferes with the freedom to do exploration, as discussed earlier.

(b) Assume that $C^{(\lambda)}$ is invertible, but is nearly singular. Then
the matrix inversion in the LSTD($\lambda$) formula (3.9) may require a very
large number of samples to yield a reasonably accurate solution of
$C^{(\lambda)}r=d^{(\lambda)}$.††† It is well-known from fundamental error
analyses of linear equation solvers that small errors in a nearly singular
matrix $C^{(\lambda)}$ will cause large errors in the solution of
$C^{(\lambda)}r=d^{(\lambda)}$. Near-singularity of $C^{(\lambda)}$ may be due
either to the columns of $\Phi$ being nearly linearly dependent or to the
matrix $\Xi(I-\alpha P^{(\lambda)})$ being nearly singular [cf. Eq. (3.7)].
Near-linear dependence of the columns of $\Phi$ will not affect the error in
the solution of the high-dimensional projected equation, which can be written
as $\Phi C^{(\lambda)}r=\Phi d^{(\lambda)}$. The reason is that this error
depends only on the subspace $S$ and not its representation in terms of the
matrix $\Phi$. In particular, if we replace $\Phi$ with a matrix $\Phi B$
where $B$ is an $s\times s$ invertible scaling matrix, the subspace $S$ will
be unaffected and the error in the solution of the projected equation will
also be unaffected. On the other hand, near singularity of the matrix
$I-\alpha P^{(\lambda)}$ may affect significantly the error. Note that
$I-\alpha P^{(\lambda)}$ is nearly singular in the case where $\alpha$ is very
close to 1, or in the corresponding undiscounted case where $\alpha=1$ and $P$
is substochastic with some eigenvalues very close to 1. Large variations in
the size of the diagonal components of $\Xi$ may also affect significantly the
error, although this dependence is complicated by the fact that $\Xi$ appears
not only in the formula $C^{(\lambda)}=\Phi^{\prime}\Xi(I-\alpha
P^{(\lambda)})\Phi$ but also in the formula $d^{(\lambda)}=\Phi^{\prime}\Xi
g^{(\lambda)}$. To correct the sensitivity of LSTD($\lambda$) to simulation
noise, it may be necessary to turn it into an iterative method through some
form of regularization, which then brings it close to a form of
LSPE($\lambda$) (see [Ber09], [WPB09], [Ber11a], [Ber11b], [Ber11c] for such
regularization methods and their connection to LSPE). Of course, the situation
becomes even more complex if $C^{(\lambda)}$ is singular, perhaps due to
inadvertent rank deficiency of $\Phi$ (see [WaB11a], [WaB11b] for a discussion
of this possibility).

(c) When LSTD($\lambda$) and LSPE($\lambda$) are embedded within a PI
framework, the number of samples collected using any one policy is often
relatively small. Then the behavior of the two methods becomes very
complicated, and it is hard to reach any kind of reliable conclusion [Ber10].
Computational studies indicate that LSPE($\lambda$) being an iterative method,
is less sensitive to the matrix inversion errors that afflict LSTD($\lambda$)
in the presence of high simulation noise.

The preceding discussion is also relevant to the implementations of
$\lambda$-PI to be discussed in the next section, since these implementations
bear strong relations to both LSTD($\lambda$) and LSPE($\lambda$).

4\. LAMBDA-POLICY ITERATION WITH COST FUNCTION APPROXIMATION

We saw in Section 2 that the policy evaluation portion of $\lambda$-PI,

$J_{k+1}=T_{\mu_{k+1}}^{(\lambda)}J_{k},$ (4.1)

[cf. Eq. (1.5)] can be implemented in two ways:

(1) By computing $T_{\mu_{k+1}}^{(\lambda)}J_{k}$.

(2) By finding the fixed point of the mapping $W_{k}$ [cf. Eq. (2.1)] through
solution of the equation

$J=W_{k}J,$ (4.2)

which can be viewed as Bellman’s equation associated with the current policy
for the two equivalent DP problems discussed in Section 2 [cf. Eq. (2.2)]: a
$\lambda\alpha$-discounted problem and a stopping problem.

Let us now consider approximation of $\lambda$-PI on the subspace $S=\\{\Phi
r\mid r\in\Re^{s}\\}$. A natural possibility is to introduce projection in the
preceding approaches. In particular, we may approximate the $\lambda$-PI
iterate $J_{k+1}$ of Eq. (4.1) by $\Phi r_{k+1}$ in three ways:

(a) By using a single projected value iteration for the original
$\alpha$-discounted problem,

$\Phi r_{k+1}=\Pi T_{\mu_{k+1}}^{(\lambda)}(\Phi r_{k}).$ (4.3)

This is the original proposal of [BeI96]. It is the variant of the
LSPE($\lambda$) method (3.11)-(3.12), which involves just the last iteration.

(b) By solving a projected version of Eq. (4.2), viewing it as Bellman’s
equation for the $\lambda\alpha$-discounted problem of Section 2, and setting
$r_{k+1}$ equal to its solution. This is the proposal of [ThS10a], and
implements by simulation the solution of this projected equation, essentially
by applying LSTD(0) to Bellman’s equation for the $\lambda\alpha$-discounted
problem formulated in Section 2.

(c) By solving a projected version of Eq. (4.2), viewing it as Bellman’s
equation for the stopping problem formulated in Section 2, and setting
$r_{k+1}$ equal to its solution.

In the following three subsections, we will describe three alternative
implementations of $\lambda$-PI corresponding to the possibilities (a)-(c)
above. Of course when linear cost function approximation of the form $\Phi
r_{k}$ is used to represent $J_{k}$, the $\lambda$-PI method need not
converge, and the cost vectors $J_{\mu_{k}}$ of the generated policies
typically oscillate within some suboptimality threshold from $J^{*}$. We do
not address this issue, but we note that related error bounds, which also
apply to other forms of optimistic PI are given by Bertsekas and Yu [BeY10a],
Thiery and Scherrer [ThS10b], and Scherrer [Sch11].

4.1  The LSPE($\lambda$) Implementation

A variant of the LSPE($\lambda$) method (3.11)-(3.12) is to form batches of
simulation samples and perform iteration (3.11) at the end of each batch. In
an extreme case, we treat the entire simulation trajectory
$(i_{0},\ldots,i_{t})$ as a single simulation batch, and we perform a single
iteration (3.11), for $\ell=t$, yielding the method

$r_{k+1}=r_{k}-G_{t}\bigl{(}C_{t}^{(\lambda)}r_{k}-d_{t}^{(\lambda)}\bigr{)},$
(4.4)

where $\Phi r_{k}$ is the approximate evaluation of the cost vector of the
preceding policy $\mu_{k}$ [cf. Eq. (3.13)]. As $t\to\infty$ and the
simulation becomes exact in the limit, i.e.,

$\lim_{t\to\infty}C_{t}^{(\lambda)}=C^{(\lambda)},\qquad\lim_{t\to\infty}d_{t}^{(\lambda)}=d^{(\lambda)},$

and if $G_{t}$ is given by the formula (3.12), it can be verified that

$\Phi r_{k+1}\to\Pi T_{\mu_{k+1}}^{(\lambda)}(\Phi r_{k}).$ (4.5)

Thus the method (4.4) with $G_{t}$ given by Eq. (3.12) can be viewed as a
simulation-based implementation of Eq. (4.3), the projected version of
$\lambda$-PI, which becomes exact in the limit as $t\to\infty$. In practice of
course $t$ is finite, and one may consider variants of the method, whereby
multiple iterations of the form (4.4) are performed, with each iteration using
additional simulation samples.

We note a mathematically equivalent description of this method, which is given
in terms of a least-squares optimization (see [Ber07], Section 6.3.3 for a
more detailed textbook account): we set

$r_{k+1}=\arg\min_{r\in\Re^{s}}\sum_{\ell=0}^{t}\left(\phi(i_{\ell})^{\prime}r-\phi(i_{\ell})^{\prime}r_{k}-\sum_{m=\ell}^{t}(\lambda\alpha)^{m-\ell}q(i_{m},i_{m+1})\right)^{2},$
(4.6)

where $q(i_{m},i_{m+1})$ is the temporal difference

$q(i_{m},i_{m+1})=g\big{(}i_{m},\mu_{k+1}(i_{m}),i_{m+1}\big{)}+\alpha\phi(i_{m+1})^{\prime}r_{k}-\phi(i_{m})^{\prime}r_{k},\qquad
m=0,\ldots,t.$ (4.7)

In fact this is how the method was originally described in [BeI96] and
[BeT96].

A positive aspect of this method is that it approximates directly $\Pi
T_{\mu_{k+1}}^{(\lambda)}(\Phi r_{k})$, so it is not subject to bias in the
evaluation of the fixed point of $W_{k}$; cf. Eq. (4.5). However, in the form
given here, the method does not address the issue of exploration. Despite this
fact, this implementation [in the form (4.6)] has been successful in several
challenging computational studies, including the one involving the game of
tetris in the original paper [BeI96] and some followup works, and a recent one
by Foderaro et. al. [FRF11] involving the game of pac-man, a benchmark problem
of pursuit-evasion.

4.2  $\lambda$-PI(0) - An Implementation Based on a Discounted MDP

This implementation, suggested and tested by Thiery and Scherrer [ThS10a],
[ThS10b], is based on the fixed point property of $J_{k+1}$ [cf. Prop.
2.1(b)]. It produces an approximation $\Phi r_{k+1}$ to $J_{k+1}$ within the
subspace $S$, by solving the projected equation

$\Phi r=\Pi W_{k}(\Phi r),$

with $W_{k}$ given by

$W_{k}J=(1-\lambda)T_{\mu_{k+1}}(\Phi r_{k})+\lambda T_{\mu_{k+1}}J,\qquad
J\in\Re^{n}.$ (4.8)

We may find the solution $r_{k+1}$ of this equation by using an LSTD(0)-like
simulation approach. In particular, $r_{k+1}$ satisfies the orthogonality
condition

$Cr=d(k),$

where

$C=\Phi^{\prime}\Xi(I-\lambda\alpha P_{\mu_{k+1}})\Phi,\qquad
d(k)=\Phi^{\prime}\Xi\big{(}g_{\mu_{k+1}}+(1-\lambda)\alpha P_{\mu_{k+1}}\Phi
r_{k}\big{)},$

so that

$r_{k+1}=C^{-1}d(k).$ (4.9)

We refer to this method as $\lambda$-PI(0) to distinguish it notationally from
the method of the next subsection (the name LS$\lambda$PI was introduced for
this method in [ThS10a]).

In a simulation-based implementation, the matrix $C$ and the vector $d(k)$ are
approximated by estimates $C_{t}$ and $d_{t}(k)$. Thus this method does not
require that $\Pi T_{\mu_{k+1}}^{(\lambda)}$ is a contraction, and like LSTD,
it can deal well with the issue of exploration. The simulation samples need
not depend on the policy $\mu_{k+1}$ being evaluated, so they can be generated
only once within a PI process. On the other hand the objective of the
implementation is to approximate the next iterate of $\lambda$-PI, i.e.,
$T^{(\lambda)}_{\mu_{k+1}}(\Phi r_{k})$, and it is not clear that it is doing
this well. To see this, suppose that the iteration (4.9), or equivalently
$\Phi r_{k+1}=\Pi W_{k}(\Phi r_{k})$, is repeated an infinite number of times
so it converges to a limit $\bar{r}$, which must satisfy $\Phi\bar{r}=\Pi
W_{k}(\Phi\bar{r})$. Then using Eq. (4.8), we have

$\Phi\bar{r}=(1-\lambda)\Pi T_{\mu_{k+1}}(\Phi\bar{r})+\lambda\Pi
T_{\mu_{k+1}}(\Phi\bar{r}),$

which shows that $\Phi\bar{r}=\Pi T_{\mu_{k+1}}(\Phi\bar{r})$. Thus
$\lambda$-PI(0) aims at $\bar{r}$, which is the limit of TD(0) independent of
the value of $\lambda$. Indeed as $\lambda\to 1$, $\Pi W_{k}$ tends to $\Pi
T_{\mu_{k+1}}$ [cf. Eq. (4.8)], so its fixed point $\Phi r_{k+1}$ tends to the
fixed point of $\Pi T_{\mu_{k+1}}$, i.e., the limit of TD(0). It follows that
while this implementation deals well with the issue of exploration, it may be
subject to significant bias-related error.

4.3  $\lambda$-PI(1) - An Implementation Based on a Stopping Problem

The third implementation is based on the property mentioned in Section 2: the
fixed point equation $J=W_{k}J$ [or equivalently, Eq. (2.2)] is Bellman’s
equation for the policy $\mu_{k+1}$ in the context of a stopping problem. Here
there is an artificial termination state 0, and for all states $j$, there is
probability $1-\lambda$ that a transition to $j$ will be followed by an
immediate transition to state 0, with cost $\alpha J_{k}(j)$, cf. Eq. (2.2).
Note that if $\lambda$ is not too close to 1, the trajectories of this problem
tend to be short, and in fact if $\lambda=0$ all trajectories consist of a
single transition.

To compute an approximation $\Phi r_{k+1}$ to the fixed point of $W_{k}$ by
using the stopping problem, we may use any policy evaluation algorithm with
cost function approximation over the subspace $S=\\{\Phi r\mid
r\in\Re^{s}\\}$. An interesting choice is to use the LSPE(1) method, which
consists of a least squares fit of $\Phi r$ to the simulated costs of the
trajectories of the stopping problem whose Bellman equation mapping is
$W_{k}$. The use of LSPE(1) not only involves minimum bias relative to all
LSPE($\nu$) methods with $\nu\in[0,1]$, but also leads to a simple least
squares implementation.

To this end, we introduce a simulation procedure, called geometric sampling,
which departs from the single infinitely long simulation trajectory format of
the implementation of Section 4.1, and has the following characteristics:

(a) It uses multiple relatively short simulation trajectories.

(b) The initial state of each trajectory is chosen essentially as desired,
thereby allowing flexibility to generate a richer mixture of state visits.

(c) The length of each trajectory is random and is determined by a
$\lambda$-dependent geometric distribution [a probability
$(1-\lambda)\lambda^{\ell}$ that the number of transitions is $\ell+1$].

In particular, given the current representation $\Phi r_{k}$ of $J_{k}$ and
the current policy $\mu_{k+1}$, we update the parameter vector from $r_{k}$ to
$r_{k+1}$ after generating $t$ simulated trajectories. The states of a
trajectory are generated according to the transition probabilities
$p_{ij}\big{(}\mu_{k+1}(i)\big{)}$, the transition cost is discounted by an
additional factor $\alpha$ with each transition, and following each transition
to a state $j$, the trajectory is terminated with probability $1-\lambda$ and
with an extra cost $\alpha\phi(i)^{\prime}r_{k}$. Once a trajectory is
terminated, an initial state for the next trajectory is chosen according to a
fixed probability distribution
$\zeta_{0}=\bigl{(}\zeta_{0}(1),\ldots,\zeta_{0}(n)\bigr{)}$, and the process
is repeated. Note that the sequence of restart states need not depend on the
policy being evaluated, so that it can be simulated only once within a PI
process. Of course, the simulated trajectories have to be recalculated for
each new policy. The details are as follows.

Let the $m$th trajectory, $m=1,\ldots,t$, have the form
$(i_{0,m},i_{1,m},\ldots,i_{N_{m},m})$, where $i_{0,m}$ is the initial state,
and $i_{N_{m},m}$ is the state at which the trajectory is completed (the last
state prior to termination). For each state $i_{\ell,m}$,
$\ell=0,\ldots,N_{m}-1$, of the $m$th trajectory, the simulated cost is

$c_{\ell,m}(r_{k})=\alpha^{N_{m}-\ell}\phi(i_{N_{m},m})^{\prime}r_{k}+\sum_{q=\ell}^{N_{m}-1}\alpha^{q-\ell}g(i_{q,m},u_{q,m},i_{q+1,m}),$
(4.10)

where

$u_{q,m}=\mu_{k+1}(i_{q,m}),\qquad m=1,\ldots,t,\ q=0,\ldots,N_{m}-1.$

Once the costs $c_{\ell,m}(r_{k})$ are computed for all states $i_{\ell,m}$ of
the $m$th trajectory and all trajectories $m=1,\ldots,t$, the vector $r_{k+1}$
is obtained by a least squares fit of these costs:

$r_{k+1}=\arg\min_{r\in\Re^{s}}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\big{(}\phi(i_{\ell,m})^{\prime}r-c_{\ell,m}(r_{k})\big{)}^{2},$
(4.11)

cf. Eqs. (4.6)-(4.7). Equivalently, we can write the solution of the least
squares problem explicitly as

$r_{k+1}=\left(\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})\phi(i_{\ell,m})^{\prime}\right)^{-1}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})c_{\ell,m}(r_{k}).$
(4.12)

We refer to the resulting implementation as $\lambda$-PI(1).

Note the extreme special case when $\lambda=0$. Then all the simulated
trajectories consist of a single transition, and there is a restart at every
transition. This means that the simulation samples are from states that are
generated independently according to the restart distribution $\zeta_{0}$.

Convergence of the Simulation Process

We will now show that in the limit, as $t\to\infty$, the vector $r_{k+1}$ of
Eq. (4.12) satisfies

$\Phi r_{k+1}=\hat{\Pi}T_{\mu_{k+1}}^{(\lambda)}(\Phi r_{k}),$ (4.13)

where $\hat{\Pi}$ denotes projection with respect to the weighted sup-norm
$\|\cdot\|_{\zeta}$ with weight vector
$\zeta=\big{(}\zeta(1),\ldots,\zeta(n)\big{)}$, where

$\zeta(i)={\hat{\zeta}(i)\over\sum_{j=1}^{n}\hat{\zeta}(j)},\qquad
i=1,\ldots,n,$ (4.14)

and

$\hat{\zeta}(i)=\sum_{\ell=0}^{\infty}\zeta_{\ell}(i),$

with $\zeta_{\ell}(i)$ being the probability of the state being $i$ after
$\ell$ transitions of a randomly chosen simulation trajectory. This is the
underlying norm in TD methods such as LSTD, LSPE, and TD, as applied to SSP
problems (see [BeT96], Section 6.3.4). Note that $\zeta(i)$ is the long-term
occupancy probability of state $i$ during the simulation process. We assume
that the restart distribution $\zeta_{0}$ is chosen so that $\zeta(i)>0$ for
all $i=1,\ldots,n$, implying that $\|\cdot\|_{\zeta}$ is a legitimate norm
[this is guaranteed if we require that $\zeta_{0}(i)>0$ for all $i$].

Indeed, let us view $T^{\ell+1}_{\mu_{k+1}}J$ as the vector of total
discounted costs over a horizon of $(\ell+1)$ stages with the terminal cost
function being $J$, and write

$T^{\ell+1}_{\mu_{k+1}}J=\alpha^{\ell+1}P^{\ell+1}_{\mu_{k+1}}J+\sum_{q=0}^{\ell}\alpha^{q}P_{\mu_{k+1}}^{q}g_{\mu_{k+1}},$

where $P_{\mu_{k+1}}$ and $g_{\mu_{k+1}}$ are the transition probability
matrix and cost vector, respectively, under $\mu_{k+1}$. As a result the
vector
$T_{\mu_{k+1}}^{(\lambda)}J=(1-\lambda)\sum_{\ell=0}^{\infty}\lambda^{\ell}T_{\mu_{k+1}}^{\ell+1}J$
can be expressed as

$\big{(}T_{\mu_{k+1}}^{(\lambda)}J\big{)}(i)=\sum_{\ell=0}^{\infty}(1-\lambda)\lambda^{\ell}E\left\\{\alpha^{\ell+1}J(i_{\ell+1})+\sum_{q=0}^{\ell}\alpha^{q}g\big{(}i_{q},\mu_{k+1}(i_{q}),i_{q+1}\big{)}\
\Big{|}\ i_{0}=i\right\\},\qquad i=1,\ldots,n.$

Thus $\big{(}T_{\mu_{k+1}}^{(\lambda)}J\big{)}(i)$ may be viewed as the
expected value of the $(\ell+1)$-stages cost of policy $\mu_{k+1}$ starting at
state $i$, with the number of stages being random and geometrically
distributed with parameter $\lambda$ [probability of $\kappa+1$ transitions is
$(1-\lambda)\lambda^{\kappa}$, $\kappa=0,1,\ldots$]. It follows that the cost
samples $c_{\ell,m}(r_{k})$ of Eq. (4.10), produced by the simulation process
described earlier, can be used to estimate
$\big{(}T_{\mu_{k+1}}^{(\lambda)}(\Phi r_{k})\big{)}(i)$ for all $i$ by Monte
Carlo averaging. The estimation formula is

$D_{t}(i)={1\over\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)}\cdot\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)c_{\ell,m}(r_{k}),$
(4.15)

where $\delta(i_{\ell,m}=i)=1$ if $i_{\ell,m}=i$ and $\delta(i_{\ell,m}=i)=0$
otherwise, and we have

$\big{(}T_{\mu_{k+1}}^{(\lambda)}(\Phi
r_{k})\big{)}(i)=\lim_{t\to\infty}D_{t}(i),\qquad i=1,\ldots,n,$

(see also the discussion on the consistency of Monte Carlo simulation for
policy evaluation in [BeT96], Section 5.2).

Let us now compare the $\lambda$-PI iteration (4.13) with the simulation-based
implementation (4.12). Using the definition of projection, Eq. (4.13) can be
written as

$r_{k+1}=\arg\min_{r\in\Re^{s}}\sum_{i=1}^{n}\zeta(i)\Big{(}\phi(i)^{\prime}r-\big{(}T_{\mu_{k+1}}^{(\lambda)}(\Phi
r_{k})\big{)}(i)\Big{)}^{2},$

or equivalently

$r_{k+1}=\left(\sum_{i=1}^{n}\zeta(i)\phi(i)\phi(i)^{\prime}\right)^{-1}\sum_{i=1}^{n}\zeta(i)\phi(i)\big{(}T_{\mu_{k+1}}^{(\lambda)}(\Phi
r_{k})\big{)}(i).$ (4.16)

Let $\tilde{\zeta}(i)$ be the empirical relative frequency of state $i$ during
the simulation, given by

$\tilde{\zeta}(i)={1\over
N_{1}+\cdots+N_{t}}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i).$
(4.17)

Then the simulation-based estimate (4.12) can be written as

$\eqalign{r_{k+1}&=\left(\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})\phi(i_{\ell,m})^{\prime}\right)^{-1}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})c_{\ell,m}(r_{k})\cr&=\left(\sum_{i=1}^{n}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)\phi(i)\phi(i)^{\prime}\right)^{-1}\sum_{i=1}^{n}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)\phi(i)c_{\ell,m}(r_{k})\cr&=\left(\sum_{i=1}^{n}\tilde{\zeta}(i)\phi(i)\phi(i)^{\prime}\right)^{-1}\sum_{i=1}^{n}{1\over
N_{1}+\cdots+N_{t}}\cdot\phi(i)\cdot\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)c_{\ell,m}(r_{k})\cr&=\left(\sum_{i=1}^{n}\tilde{\zeta}(i)\phi(i)\phi(i)^{\prime}\right)^{-1}\sum_{i=1}^{n}{\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)\over
N_{1}+\cdots+N_{t}}\cdot\phi(i)\cdot\cr&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\cdot{1\over\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)}\cdot\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\delta(i_{\ell,m}=i)c_{\ell,m}(r_{k})\cr}$

and finally, using Eqs. (4.15) and (4.17),

$r_{k+1}=\left(\sum_{i=1}^{n}\tilde{\zeta}(i)\phi(i)\phi(i)^{\prime}\right)^{-1}\sum_{i=1}^{n}\tilde{\zeta}(i)\phi(i)D_{t}(i).$
(4.18)

We can now compare the $\lambda$-PI iteration (4.16) and the simulation-based
implementation (4.18). Since $\big{(}T_{\mu_{k+1}}^{(\lambda)}(\Phi
r_{k})\big{)}(i)=\lim_{t\to\infty}D_{t}(i)$ and
$\zeta(i)=\lim_{t\to\infty}\tilde{\zeta}(i)$, we see that these two iterations
asymptotically coincide.

The expression (4.18) provides some insight on how $\lambda$-PI(1)
approximates the $\lambda$-PI iteration (4.16) [or equivalently $\Phi
r_{k+1}=\hat{\Pi}T^{(\lambda)}_{\mu_{k+1}}(\Phi r_{k})$; cf. Eq. (4.13)].
Generally the simulation process of $\lambda$-PI(1) (many short trajectories)
involves more noise than the simulation process of the other implementations
(a single long trajectory), because the length of each simulation trajectory
is random (exponentially distributed). This can be seen from iteration (4.18),
which involves considerable simulation noise due to the presence of
$\tilde{\zeta}(i)$ and $D_{t}(i)$. However, we will argue that from a
practical point of view much of this noise does not play a significant role.
To see this, first note that the deviation of $\tilde{\zeta}(i)$ from
$\zeta(i)$, is not important since $\tilde{\zeta}(i)$ simply redefines the
projection norm. Next note that $D_{t}(i)$ can be written as

$D_{t}(i)=\sum_{\ell=0}^{\infty}\tilde{f}_{\ell}(i)\tilde{E}_{\ell}(i),$
(4.19)

where $\tilde{f}_{\ell}(i)$ and $\tilde{E}_{\ell}(i)$ are the following
empirical averages over the entire simulation process:

(a) $\tilde{f}_{\ell}(i)$ is the empirical relative frequency of cost samples
that start at state $i$, and correspond to trajectories consisting of $\ell+1$
transitions. As $t\to\infty$ it converges to $(1-\lambda)\lambda^{\ell}$ based
on the way the simulation is structured.

(b) $\tilde{E}_{\ell}(i)$ is the Monte Carlo estimate of the cost of
trajectories that start at state $i$, consist of $\ell+1$ transitions, and
have terminal cost vector $\Phi r_{k}$. As $t\to\infty$ it converges to
$T^{\ell+1}_{\mu_{k+1}}(\Phi r_{k})(i)$.

While both $\tilde{f}_{\ell}(i)$ and $\tilde{E}_{\ell}(i)$ contribute to the
variance of $D_{t}(i)$, only $\tilde{E}_{\ell}(i)$ has practical significance.
To see this note that based on Eq. (4.19), $D_{t}(i)$ can also be viewed as an
estimate of

$\tilde{T}_{\mu_{k+1}}(\Phi
r_{k})(i)=\sum_{\ell=0}^{\infty}\tilde{f}_{\ell}(i)T^{\ell+1}_{\mu_{k+1}}(\Phi
r_{k})(i).$ (4.20)

Thus iteration (4.18) may also be viewed as a simulation-based implementation
of the optimistic PI method

$\Phi r_{k+1}=\tilde{\Pi}\tilde{T}_{\mu_{k+1}}(\Phi r_{k}),$

where $\tilde{\Pi}$ is projection with respect to the weighted sup-norm
defined by $\tilde{\zeta}$. From a practical point of view, this iteration and
the $\lambda$-PI iteration $\Phi
r_{k+1}=\hat{\Pi}T^{(\lambda)}_{\mu_{k+1}}(\Phi r_{k})$ perform similarly:
there is only a difference in the projection norm ($\tilde{\Pi}$ rather than
$\hat{\Pi}$), and a difference in the weights of the terms
$T^{\ell+1}_{\mu_{k+1}}$ [$\tilde{f}_{\ell}(i)$ rather than
$(1-\lambda)\lambda^{\ell}$]; compare $\tilde{T}_{\mu_{k+1}}(\Phi r_{k})(i)$
as given by Eq. (4.20) with

$T^{(\lambda)}_{\mu_{k+1}}(\Phi
r_{k})(i)=\sum_{\ell=0}^{\infty}(1-\lambda)\lambda^{\ell}T^{\ell+1}_{\mu_{k+1}}(\Phi
r_{k})(i),$

the definition of $T^{(\lambda)}_{\mu_{k+1}}$. Neither difference should
affect significantly the quality of the obtained approximation $\Phi r_{k+1}$.

In conclusion, with the $\lambda$-PI(1) implementation (4.10)-(4.12), as
$t\to\infty$, we obtain in the limit the $\lambda$-PI iteration Eq. (4.13),
with comparable performance degradation due to simulation noise as for the
LSPE($\lambda$) implementation of Section 4.1. A key characteristic of the
implementation is that it deals with the issue of exploration flexibly and
effectively. Since a trajectory of the stopping problem is completed at each
transition with the potentially large probability $1-\lambda$, a restart with
a new initial state $i_{0}$ is frequent and the length of each of the
simulated trajectories is relatively small. The restart mechanism can be used
as a “natural” form of exploration, by choosing appropriately the restart
distribution $\zeta_{0}$ so that $\zeta(i)$ reflects a “substantial” weight
for all states $i$. Thus $\lambda$-PI(1) is like LSPE($\lambda$) (Section
4.1), but with built-in exploration enhancement. Compared to $\lambda$-PI(0)
(Section 4.2) it involves reduced bias since it aims to find the limit point
of TD($\lambda$), not TD(0). In particular, as $\lambda\to 1$, it produces an
evaluation $\Phi r_{k+1}$ that tends to the fixed point of TD(1), i.e., the
projection $\hat{\Pi}J_{\mu_{k+1}}$.

4.4  Comparison with Alternative Approximate PI Methods

The preceding $\lambda$-PI implementations are in direct competition with
approximate PI methods that use LSTD($\lambda$) for policy evaluation. A
popular method, often referred to as LSPI (Lagoudakis and Parr [LaP03]), can
be simply described as approximate PI combined with LSTD(0) for policy
evaluation. The LSPI and $\lambda$-PI(0) methods have been compared in
[ThS10a] in terms of four characteristics.

(a) Bias: Both methods are subject to qualitatively similar bias [they aim to
find the limit point of TD(0)].

(b) Sample efficiency: Both methods can reuse the same set of sample state
trajectories over all policies. (In the model-free case where Q-factors are
approximated, again the set of sample state-control trajectories is reusable.)

(c) Exploration: Both methods provide the same options for exploration, since
the validity of these methods does not depend on whether the simulation
trajectories are obtained by using the current policy [in fact these
trajectories are reusable as per (b) above].

(d) Optimistic operation: Since $\lambda$-PI(0) has an iterative character
($r_{k+1}$ depends on $r_{k}$), it is less susceptible to simulation noise and
has an advantage over LSPI in the case where the number of samples per policy
is low. Indeed this assertion is made by Thiery and Scherrer [ThS10a] based on
experimentation, who also found that the effect of the choice of $\lambda$ is
more pronounced in this case.

Note that (b) and (c) above are the advantages of LSPI and $\lambda$-PI(0)
over the LSPE($\lambda$) implementation of Section 4.1 (which in turn involves
less bias because of the use of $\lambda>0$, and also has an optimistic
character).

Let us now compare $\lambda$-PI(1) with LSPI and $\lambda$-PI(0) in terms of
the characteristics (a)-(d) above. It has better bias characteristics as noted
earlier. It has worse sample efficiency as it cannot reuse simulation
trajectories (it can only reuse the restart state sequence). It deals with
exploration about as well, thanks to the restart mechanism of the SSP
formulation. Finally, like $\lambda$-PI(0), $\lambda$-PI(1) has an optimistic
character, and has a similar advantage over LSPI in this regard, cf. (d)
above.

4.5  Exploration-Enhanced LSTD($\lambda$) with Geometric Sampling

The geometric sampling idea underlying the $\lambda$-PI(1) implementation of
Eqs. (4.10)-(4.12) may also be modified to obtain an exploration-enhanced
version of LSTD($\lambda$). In particular, we use the same simulation
procedure, and in analogy to Eq. (4.10) we define

$c_{\ell,m}(r)=\alpha^{N_{m}-\ell}\phi(i_{N_{m},m})^{\prime}r+\sum_{q=\ell}^{N_{m}-1}\alpha^{q-\ell}g(i_{q,m},u_{q,m},i_{q+1,m}).$

We then obtain an approximation $\Phi\hat{r}$ to the solution of the projected
equation

$\Phi r=\hat{\Pi}T_{\mu_{k+1}}^{(\lambda)}(\Phi r),$

[cf. Eq. (4.13)] by finding $\hat{r}$ such that

$\hat{r}=\arg\min_{r\in\Re^{s}}\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\big{(}\phi(i_{\ell,m})^{\prime}r-c_{\ell,m}(\hat{r})\big{)}^{2}.$
(4.21)

By writing the optimality condition

$\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})\big{(}\phi(i_{\ell,m})^{\prime}\hat{r}-c_{\ell,m}(\hat{r})\big{)}=0$

for the least squares minimization in Eq. (4.21) and solving for $\hat{r}$, we
obtain the following implementation of LSTD($\lambda$):

$\hat{r}=\hat{C}^{-1}\hat{d},$ (4.22)

where

$\hat{C}=\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})\big{(}\phi(i_{\ell,m})-\alpha^{N_{m}-\ell}\phi(i_{N_{m},m})\big{)}^{\prime},$
(4.23)

and

$\hat{d}=\sum_{m=1}^{t}\sum_{\ell=0}^{N_{m}-1}\phi(i_{\ell,m})\sum_{q=\ell}^{N_{m}-1}\alpha^{q-\ell}g(i_{q,m},u_{q,m},i_{q+1,m}).$
(4.24)

For a large number of trajectories $t$, the exploration-enhanced
LSTD($\lambda$) method (4.21) [or equivalently (4.22)-(4.24)] and
$\lambda$-PI(1) [cf. Eq. (4.12)] yield similar results, particularly when
$\lambda\approx 1$. However, $\lambda$-PI(1) has an iterative character
($r_{k+1}$ depends on $r_{k}$), so it is reasonable to expect that it is less
susceptible to simulation noise in an optimistic PI setting where the number
of samples per policy is low.

As an example, when $\lambda=0$, all the simulation trajectories consist of a
single transition, so $N_{m}=1$ for all $m=1,\ldots,t$. Then, using Eqs.
(4.23) and (4.24), the equation $\hat{C}r=\hat{d}$ becomes

$\sum_{m=1}^{t}\phi(i_{0,m})\big{(}\phi(i_{0,m})-\alpha\phi(i_{1,m})\big{)}^{\prime}r=\sum_{m=1}^{t}\phi(i_{0,m})g(i_{0,m},u_{0,m},i_{1,m}).$

It yields the same vector $\hat{r}=\hat{C}^{-1}\hat{d}$ as the LSTD(0) method
that simulates $t$ independent transitions according to the restart
distribution $\zeta_{0}$, rather than simulating a single long trajectory. In
fact this is the policy evaluation process in the LSPI method mentioned in
Section 4.4. The geometric sampling procedure described here allows
exploration-enhancement for any $\lambda$.

5\. CONCLUSIONS

We discussed a few implementations of $\lambda$-PI with linear cost function
approximation, which have different strengths and weaknesses with respect to
dealing with the critical issues of bias and exploration. Out of the three
implementations, the one of Section 4.3, $\lambda$-PI(1), is new and seems
capable of dealing well with both issues, although it has worse sample
complexity than the $\lambda$-PI(0) implementation of Section 4.2.

On the other hand, our discussion has been somewhat speculative, and our
assessments, while relying on past computational experience, still require
supportive experimentation. Moreover, the $\lambda$-PI implementations should
be compared to other approximate PI methods based on projected equations, such
as the exploration-enhanced LSTD($\lambda$) method for policy evaluation,
discussed in Section 3, and the LSPI method discussed in Section 4.4. A
computational comparison of $\lambda$-PI(0) with this latter method is given
in [ThS10a], and a similar comparison with $\lambda$-PI(1) would be desirable.

Fundamentally, $\lambda$-PI(1) is based on geometric sampling, a new
simulation idea for $\lambda$-methods that uses multiple short trajectories
with exploration-enhanced restart, rather than a single infinitely long
trajectory. This idea can also be applied to LSTD($\lambda$), thereby
obtaining a new exploration-enhanced version of this method, which has been
described in Section 4.5.

6\. REFERENCES

[BBD10] Busoniu, L., Babuska, R., De Schutter, B., and Ernst, D., 2010.
Reinforcement Learning and Dynamic Programming Using Function Approximators,
CRC Press, N. Y.

[BBN04] Bertsekas, D. P., Borkar, V. S., and Nedić, A., 2004. “Improved
Temporal Difference Methods with Linear Function Approximation,” in Learning
and Approximate Dynamic Programming, by J. Si, A. Barto, W. Powell, and D.
Wunsch (Eds.), IEEE Press, N. Y.

[BSA83] Barto, A. G., Sutton, R. S., and Anderson, C. W., 1983. “Neuronlike
Elements that Can Solve Difficult Learning Control Problems,” IEEE Trans. on
Systems, Man, and Cybernetics, Vol. 13, pp. 835-846.

[BeI96] Bertsekas, D. P., and Ioffe, S., 1996. “Temporal Differences-Based
Policy Iteration and Applications in Neuro-Dynamic Programming,” Lab. for
Info. and Decision Systems Report LIDS-P-2349, MIT.

[BeT96] Bertsekas, D. P., and Tsitsiklis, J. N., 1996\. Neuro-Dynamic
Programming, Athena Scientific, Belmont, MA.

[BeY09] Bertsekas, D. P., and Yu, H., 2009. “Projected Equation Methods for
Approximate Solution of Large Linear Systems,” Journal of Computational and
Applied Mathematics, Vol. 227, pp. 27-50.

[BeY10a] Bertsekas, D. P., and Yu, H., 2010. “Q-Learning and Enhanced Policy
Iteration in Discounted Dynamic Programming,” Lab. for Information and
Decision Systems Report LIDS-P-2831, MIT.

[BeY10b] Bertsekas, D. P., and Yu, H., 2010. “Asynchronous Distributed Policy
Iteration in Dynamic Programming,” Proc. of Allerton Conf. on Information
Sciences and Systems.

[Ber95] Bertsekas, D. P., 1995. “A Counterexample to Temporal Differences
Learning,” Neural Computation, Vol. 7, pp. 270-279.

[Ber07] Bertsekas, D. P., 2007. Dynamic Programming and Optimal Control, 3rd
Edition, Vol. II, Athena Scientific, Belmont, MA.

[Ber09] Bertsekas, D. P., 2009. “Projected Equations, Variational
Inequalities, and Temporal Difference Methods,” Lab. for Information and
Decision Systems Report LIDS-P-2808, MIT.

[Ber10] Bertsekas, D. P., 2010. “Pathologies of Temporal Difference Methods in
Approximate Dynamic Programming,” Proc. 2010 IEEE Conference on Decision and
Control.

[Ber11a] Bertsekas, D. P., 2011. Approximate Dynamic Programming, on-line at
http://web.mit.edu/dimitrib/www/dpchapter.html.

[Ber11b] Bertsekas, D. P., 2011. “Approximate Policy Iteration: A Survey and
Some New Methods,” J. of Control Theory and Applications, Vol. 9, pp. 310-335.

[Ber11c] Bertsekas, D. P., 2011. “Temporal Difference Methods for General
Projected Equations,” IEEE Trans. on Aut. Control, Vol. 56, pp. 2128-2139.

[Bor08] Borkar, V. S., 2008. Stochastic Approximation: A Dynamical Systems
Viewpoint, Cambridge Univ. Press.

[Bor09] Borkar, V. S., 2009. “Reinforcement Learning: A Bridge Between
Numerical Methods and Monte Carlo,” in World Scientific Review, Vol. 9,
Chapter 4.

[BrB96] Bradtke, S. J., and Barto, A. G., 1996. “Linear Least-Squares
Algorithms for Temporal Difference Learning,” Machine Learning, Vol. 22, pp.
33-57.

[CFH07] Chang, H. S., Fu, M. C., Hu, J., Marcus, S. I., 2007. Simulation-Based
Algorithms for Markov Decision Processes, Springer, N. Y.

[CaR11] Canbolat, P. G., and Rothblum, U. G., 2011. “(Approximate) Iterated
Successive Approximations Algorithm for Sequential Decision Processes,”
Technical Report, The Technion - Israel Institute of Technology, May 2011.

[Cao07] Cao, X. R., 2007. Stochastic Learning and Optimization: A Sensitivity-
Based Approach, Springer, N. Y.

[FRF11] Foderaro, G., Raju, V., and Ferrari, S., 2011. “A Model-Based
Approximate $\lambda$-Policy Iteration Approach to Online Evasive Path
Planning and the Video Game Ms. Pac-Man,” J. of Control Theory and
Applications, Vol. 9, pp. 391-399.

[Fle84] Fletcher, C. A. J., 1984. Computational Galerkin Methods, Springer-
Verlag, N. Y.

[Gos03] Gosavi, A., 2003. Simulation-Based Optimization Parametric
Optimization Techniques and Reinforcement Learning, Springer-Verlag, N. Y.

[Hay08] Haykin, S., 2008. Neural Networks and Learning Machines (3rd Edition),
Prentice-Hall, Englewood-Cliffs, N. J.

[Kra72] Krasnoselskii, M. A., et. al, 1972. Approximate Solution of Operator
Equations, Translated by D. Louvish, Wolters-Noordhoff Pub., Groningen.

[LLL08] Lewis, F. L., Lendaris, G. G., and Liu, D., 2008. Special Issue on
Adaptive Dynamic Programming and Reinforcement Learning in Feedback Control,
IEEE Transactions on Systems, Man, and Cybernetics, Vol. 38.

[LaP03] Lagoudakis, M. G., and Parr, R., 2003. “Least-Squares Policy
Iteration,” J. of Machine Learning Research, Vol. 4, pp. 1107-1149

[LeV09] Lewis, F. L., and Vrabie, D., 2009. “Reinforcement Learning and
Adaptive Dynamic Programming for Feedback Control,” IEEE Circuits and Systems
Magazine, 3rd Q. Issue.

[Mey07] Meyn, S., 2007. Control Techniques for Complex Networks, Cambridge
University Press, N. Y.

[NeB03] Nedić, A., and Bertsekas, D. P., 2003. “Least Squares Policy
Evaluation Algorithms with Linear Function Approximation,” Discrete Event
Dynamic Systems: Theory and Applications, Vol. 13, pp. 79-110.

[Pow07] Powell, W. B., 2007. Approximate Dynamic Programming: Solving the
Curses of Dimensionality, Wiley, N. Y.

[Put94] Puterman, M. L., 1994. Markov Decision Processes: Discrete Stochastic
Dynamic Programming, J. Wiley, N. Y.

[Rot79] Rothblum, U. G., 1979. “Iterated Successive Approximation for
Sequential Decision Processes,” in Stochastic Control and Optimization, by J.
W. B. van Overhagen and H. C. Tijms (eds), Vrije University, Amsterdam.

[SBP04] Si, J., Barto, A., Powell, W., and Wunsch, D., (Eds.) 2004. Learning
and Approximate Dynamic Programming, IEEE Press, N. Y.

[Sam59] Samuel, A. L., 1959. “Some Studies in Machine Learning Using the Game
of Checkers,” IBM Journal of Research and Development, pp. 210-229.

[Sam67] Samuel, A. L., 1967. “Some Studies in Machine Learning Using the Game
of Checkers. II – Recent Progress,” IBM Journal of Research and Development,
pp. 601-617.

[Sch10] Scherrer, B., 2010. “Should One Compute the Temporal Difference Fix
Point or Minimize the Bellman Residual? The Unified Oblique Projection View,”
in ICML’10: Proc. of the 27th Annual International Conf. on Machine Learning.

[Sch11] Scherrer, B., 2011. “Performance Bounds for Lambda Policy Iteration
and Application to the Game of Tetris,” Report RR-6348, INRIA.

[SuB98] Sutton, R. S., and Barto, A. G., 1998. Reinforcement Learning, MIT
Press, Cambridge, MA.

[Sut88] Sutton, R. S., 1988. “Learning to Predict by the Methods of Temporal
Differences,” Machine Learning, Vol. 3, pp. 9-44.

[SzL06] Szita, I., and Lorinz, A., 2006. “Learning Tetris Using the Noisy
Cross-Entropy Method,” Neural Computation, Vol. 18, pp. 2936-2941.

[Sze10] Szepesvari, C., 2010. “Reinforcement Learning Algorithms for MDPs,”
Morgan and Claypool Publishers.

[ThS09] Thiery, C., and Scherrer, B., 2009. “Improvements on Learning Tetris
with Cross-Entropy,” International Computer Games Association Journal, Vol.
32, pp. 23-33.

[ThS10a] Thiery, C., and Scherrer, B., 2010. “Least-Squares Policy Iteration:
Bias-Variance Trade-off in Control Problems,” Proc. of 2010 ICML, Haifa,
Israel.

[ThS10b] Thiery, C., and Scherrer, B., 2010. “Performance Bound for
Approximate Optimistic Policy Iteration,” Technical Report, INRIA.

[TsV97] Tsitsiklis, J. N., and Van Roy, B., 1997. “An Analysis of Temporal-
Difference Learning with Function Approximation,” IEEE Transactions on
Automatic Control, Vol. 42, pp. 674–690.

[WPB09] Wang, M., Polydorides, N., and Bertsekas, D. P., 2009. “Approximate
Simulation-Based Solution of Large-Scale Least Squares Problems,” Lab. for
Information and Decision Systems Report LIDS-P-2819, MIT.

[WaB11a] Wang, M., and Bertsekas, D. P., 2011. “Stabilization of Simulation-
Based Iterative Methods for Singular and Nearly Singular Linear Systems,” Lab.
for Information and Decision Systems Report LIDS-P-2878, MIT.

[WaB11b] Wang, M., and Bertsekas, D. P., 2011. “On the Convergence of
Iterative Simulation-Based Methods for Singular Linear Systems,” Lab. for
Information and Decision Systems Report LIDS-P-2879, MIT.

[Wer09] Werbos, P. J., 2009. “Intelligence in the Brain: A Theory of how it
Works and how to Build it,” Neural Networks, Vol. 22, pp. 200-212.

[WhS92] White, D., and Sofge, D., 1992. Handbook of Intelligent Control, Van
Nostrand Reinhold, N.Y.

[WiB93] Williams, R. J., and Baird, L. C., 1993. “Analysis of Some Incremental
Variants of Policy Iteration: First Steps Toward Understanding Actor-Critic
Learning Systems,” Report NU-CCS-93-11, College of Computer Science,
Northeastern University, Boston, MA.

[YuB09] Yu, H., and Bertsekas, D. P., 2009. “Convergence Results for Some
Temporal Difference Methods Based on Least Squares,” IEEE Trans. on Aut.
Control, Vol. 54, 2009, pp. 1515-153.

[YuB10] Yu, H., and Bertsekas, D. P., 2010. “Error Bounds for Approximations
from Projected Linear Equations,” Mathematics of Operations Research, Vol. 35,
pp. 306-329.

[YuB11] Yu, H., and Bertsekas, D. P., 2011. “Q-Learning and Policy Iteration
Algorithms for Stochastic Shortest Path Problems,” Lab. for Information and
Decision Systems Report LIDS-P-2871, MIT.

[Yu10a] Yu, H., 2010. “Least Squares Temporal Difference Methods: An Analysis
Under General Conditions,” Technical report C-2010-39, Dept. Computer Science,
Univ. of Helsinki.

[Yu10b] Yu, H., 2010. “Convergence of Least Squares Temporal Difference
Methods Under General Conditions,” Proc. of the 27th ICML, Haifa, Israel.