2023

[1]Matilde Gargiani

[1]Automatic Control Laboratory, ETH, Physikstrasse 3, 8092 Zurich,
Switzerland

2]Control and Automation Group, Inspire AG, Zurich, Switzerland

3]Department of Mechanical Engineering, UBC, 2054-6250 Applied Science Lane,
Vancouver, British Columbia, Canada

# Inexact Policy Iteration Methods for Large-Scale Markov Decision Processes

<EMAIL_ADDRESS>Robin Sieber<EMAIL_ADDRESS>Efe Balta
<EMAIL_ADDRESS>Dominic Liao-McPherson<EMAIL_ADDRESS>John Lygeros
<EMAIL_ADDRESS>*  [  [

###### Abstract

We consider inexact policy iteration (iPI) methods for large-scale infinite-
horizon discounted Markov decision processes (MDPs) with finite state and
action spaces, a variant of policy iteration where the policy evaluation step
is implemented inexactly using an iterative solver for linear systems. In the
classical dynamic programming literature, a similar principle is deployed in
optimistic policy iteration, where an a-priori fixed-number of iterations of
value iteration is used to inexactly solve the policy evaluation step.
Inspired by the connection between policy iteration and semismooth Newton’s
method, we investigate a class of iPI methods that mimic the inexact variants
of semismooth Newton’s method by adopting a parametric stopping condition to
regulate the level of inexactness of the policy evaluation step. For this
class of methods we discuss local and global convergence properties and derive
a practical range of values for the stopping-condition parameter that provide
contraction guarantees. We also show that, when the iterative solver for
policy evaluation enjoys linear contraction, the stopping criterion is
verified in a finite-number of iterations. Our analysis is general and
therefore encompasses a variety of iterative solvers for policy evaluation,
including the standard value iteration as well as more sophisticated ones such
as GMRES. As underlined by our analysis, the selection of the inner solver is
of fundamental importance for the performance of the overall method. We
therefore consider different iterative methods to solve the policy evaluation
step and analyse their applicability and contraction properties when used for
policy evaluation. In particular, we study Richardson’s method, steepest
descent, the minimal residual method and GMRES for policy evaluation. We show
that the contraction properties of these methods tend to be enhanced by the
specific structure of policy evaluation and that there is margin for
substantial improvement over value iteration in terms of convergence rate.
Finally, we study the numerical performance of different instances of inexact
policy iteration on large-scale MDPs for the design of health policies to
control the spread of infectious diseases in epidemiology.

###### keywords:

Dynamic Programming, Inexact Semismooth Newton-Type Methods, Policy Iteration,
Regular Markov Decision Processes, GMRES

##  1 Introduction

Markov decision processes (MDPs) are a mathematical model for dynamic
decision-making problems under uncertainty puterman_mdp  and have been used
over the years in a variety of different applications, from agriculture
agriculture_MDP  to epidemiology sis_model  and finance finance_MDP , to name
a few. We refer to White1993ASO  for a comprehensive survey on applications of
MDPs.

The size of these problems is usually very large. For systems modeled with a
set of propositional state variables, the state space grows exponentially with
the number of variables. This problem becomes even more relevant for MDPs with
continuous state-spaces. If the continuous space is discretized to find a
solution, the discretization causes another level of exponential blow up. This
phenomenon is known as the curse of dimensionality Bellman:1957  and limits
the exact solution methods to problems with a compact representation which
fits into the memory RAM of the hardware at hand. For these scenarios, it is
fundamental to design exact algorithms that scale well with the dimensionality
of the problem. This is the main motivation behind optimistic policy iteration
(OPI), a variant of policy iteration (PI) which economizes in computation by
approximating the policy evaluation step using a fixed finite number of value
iteration updates (DB_book, , Chapter 2). Indeed, even if it locally enjoys a
quadratic rate of convergence gargiani_2022 , PI has expensive iterations. In
the large-scale setting, the computational costs per-iteration are mostly due
to the policy evaluation step, which requires the exact solution of a linear
system with size equal to the state space of the MDP. Consequently, even
though PI converges in very few iterations, its time performance is not
scalable as it degrades rapidly with the size of the state space of the MDP.

A similar problem arises in Newton’s method and its semismooth variants, where
inexactness is introduced in the solution of the Newtonian linear systems in
order to reduce the computational costs and, in general, the level of
inexactness is regulated via a parametric stopping condition that depends on
locally available information inexact_Newton ; MARTINEZ1995127 . Different
iterative methods for linear systems can be deployed to produce the
approximate solution and generally the choice is dictated by the specific
structure of the linear systems inexactnewtonMR ; inexactgmresnewton ;
inexactnewtoncg . Local convergence results which are independent from the
selected iterative method are available in the literature inexact_Newton ;
MARTINEZ1995127 . These inexact variants are collectively known as inexact
Newton methods and are often deployed in large-scale scenarios as their
performance is more scalable than that of Newton’s method DEUFLHARD1991366 .

In light of the connection between PI and semismooth Newton’s method
gargiani_2022 , we analyse the class of inexact policy iteration methods,
which are characterized by an approximation of the solution of the policy
evaluation step and where the level of inexactness is regulated by a
parametric stopping condition that is borrowed from the inexact variants of
Newton’s method. This class of methods has been originally proposed in
gargiani_2023 , where the authors focus on the deployment of GMRES as inner
solver, showcasing the promising performance of the resulting inexact scheme
on large-scale artificial MDPs. Building upon these results, we extend the
analysis of inexact policy iteration methods. In particular, the main
contributions of this paper can be summarized as follow:

  * •

a convergence analysis of inexact policy iteration methods which encompasses
the choice of the inner solver,

  * •

a characterization of the convergence properties of different iterative
methods for policy evaluation, which proves that value iteration for policy
evaluation is not optimal in the sense of contraction rate for a certain class
of MDPs and underlines the advantages of deploying different iterative methods
than value iteration as inner solvers,

  * •

an empirical analysis of the performance of inexact policy iteration methods
on MDPs arising from a dynamic extension of SIS models for epidemiology.

The paper is organized as follows. In Section 2 we describe the problem
setting and propose a classification of MDPs based on the properties of the
transition probability matrices of their induced Markov chains. We also
briefly discuss the classic dynamic programming methods and their convergence
and scalability properties. In Section 3 we discuss the connection between
semismooth Newton’s method and policy iteration, summarizing the main results
of gargiani_2022  and gargiani_2023 . In addition, we provide an upper-bound
on the global semismooth constant at the root of the Bellman residual function
(Lemma 10) and derive an estimate for the region of attraction of semismooth
Newton-type methods for DP (Theorem 11). Section 4 is dedicated to inexact
policy iteration methods. After a description of the methods in this class, we
provide local and global convergence guarantees. We also prove that, under
some general assumptions on the inner solver contraction properties, the
overall scheme is well-posed as the stopping condition is always met in a
finite number of inner iterations when the parameter in the stopping condition
is constant. This result stresses the importance of selecting a well-
performing inner solver for the approximate policy evaluation step. We
therefore proceed by characterizing the convergence properties of different
iterative methods for linear systems when deployed for policy evaluation. In
particular, we analyse Richardson’s method 4.2.1, the steepest descent method
4.2.2, the minimal residual method 4.2.3 and GMRES 4.2.4. Finally, in Section
5 we propose a dynamic extension of SIS models for epidemiology and illustrate
how to reformulate them as MDPs. We then study the empirical performance of
different instances of inexact policy iteration on MDPs constructed from
dynamic SIS models. Our benchmarks showcase the advantages of inexact policy
iteration methods on large-scale MDPs over policy iteration as well as the
importance of a careful selection of the inner solver.

###  1.1 Notation

We denote with $\mathbb{S}^{n\times n}\subset\mathbb{R}^{n\times n}$ the set
of $n\times n$ row-stochastic matrices. Given an $n\times n$ real matrix $A$
with $W\leq n$ distinct eigenvalues, we use $A\succ 0$ to indicate that the
matrix is symmetric positive definite, $\rho(A)$ to denote its spectral
radius, $\kappa(A)$ for its condition number, $\mu_{A}(x)$ for its minimal
polynomial, $A_{s}$ for its symmetric part, i.e.,
$A_{s}=\frac{1}{2}\left(A+A^{\top}\right)$ , $\Lambda(A)$ for its spectrum and
$\lambda_{1},\lambda_{2},\dots,\lambda_{W}$ for its distinct eigenvalues,
where

$|\lambda_{1}|\geq|\lambda_{2}|\geq\dots\geq|\lambda_{W}|\,.$

We use $a(\lambda)$ to denote the algebraic multiplicity of
$\lambda\in\Lambda(A)$, and $b(\lambda)$ for its associated degree in the
minimal polynomial, i.e.,

$\mu_{A}(x)=\Pi_{\lambda\in\Lambda(A)}(x-\lambda)^{b(\lambda)}\,.$

If $A$ is symmetric, then we use
$\lambda_{\max}(A)=\max\left\\{\lambda_{1},\,\dots,\,\lambda_{W}\right\\}$ and
$\lambda_{\min}(A)=\min\left\\{\lambda_{1},\,\dots,\,\lambda_{W}\right\\}$. In
addition, we use $a_{ij}$ and $a^{(q)}_{ij}$ for the $(i,j)$-element of $A$
and $A^{q}$, respectively. Analogously, given a vector $v\in\mathbb{R}^{n}$,
we denote with $v_{i}$ its $i$-th entry. We call non-negative (positive) a
square matrix whose elements are non-negative (positive). Given $\delta>0$ and
$x^{*}\in\mathbb{R}^{n}$, we use $\mathcal{B}(x^{*},\delta)$ to denote the set
of points $x\in\mathbb{R}^{n}$ such that $\|x-x^{*}\|_{\infty}<\delta$.

##  2 Infinite-Horizon Discounted Markov Decision Processes with Finite Spaces

A discounted MDP with finite spaces is a 5-tuple
$\left\\{\mathcal{S},\,\mathcal{A},\,P,\,g,\,\gamma\right\\}$, where
$\mathcal{S}$ and $\mathcal{A}$ are the finite sets of states and actions,
respectively;
$P:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$ is the
transition probability function, where $P(i,a,j)$ defines the probability of
transitioning to state $j$ when applying action $a$ in state $i$;
$g:\mathcal{S}\times\mathcal{A}\rightarrow[-R,\,R]$ is the stage-cost
function, which associates to each state-action pair a bounded cost; and
$\gamma\in(0,1)$ is a discount factor. Without loss of generality, we consider
$\mathcal{S}=\left\\{1,\dots,n\right\\}$ and
$\mathcal{A}=\left\\{1,\dots,m\right\\}$. Because of the presence of actions,
MDPs can be regarded as an extension of Markov chains, where the actions are
generally selected according to a given criterion called policy. There exists
different types of policies, e.g., stochastic, non-stationary etc…. In this
work we consider the subclass of deterministic stationary control policies. In
particular, a deterministic stationary control policy
$\pi:\mathcal{S}\rightarrow\mathcal{A}$ is a function that maps states to
actions. We denote with $\Pi$ the finite set comprising all the deterministic
stationary control policies, from now on simply called policies.

Consider now a generic discounted MDP with finite spaces. At time step $t$ of
the decision process under policy $\pi\in\Pi$, the system is in some state
$s_{t}$ and the action $a_{t}=\pi(s_{t})$ is applied. The discounted cost
$\gamma^{t}g(s_{t},a_{t})$ is accrued and the system transitions to a state
$s_{t+1}$ according to the probability distribution $P(s_{t},a_{t},\cdot)$.
Starting from state $i$, this process is repeated over an infinite-horizon,
leading to the following cumulative discounted cost

$V^{\pi}(i)=\lim_{T\rightarrow\infty}\mathbb{E}\left[\sum_{t=0}^{T-1}\gamma^{t}g(s_{t},\pi(s_{t}))\,\Big{|}\,s_{0}=i\right]\,,$
(1)

where $\left\\{s_{0},\pi(s_{0}),s_{1},\dots\right\\}$ is the state-action
sequence generated by the MDP under policy $\pi$ and with initial state
$s_{0}$. The expected value is taken with respect to the corresponding
transition probability measure over the space of sequences. The optimal cost
is defined as

$V^{*}(i)=\min_{\pi\in\Pi}V^{\pi}(i)\quad\forall\,i\in\mathcal{S}\,.$ (2)

Any policy $\pi^{*}\in\Pi$ that attains the optimal cost is called optimal
policy. Notice that we restrict our attention to deterministic stationary
control policies as in the considered setting there provably exists a policy
in this class that attains (2) (DB_book, , Chapter 1).

Given the cost function $V:\mathcal{S}\rightarrow\mathbb{R}$, we call greedy
with respect to the cost $V$ any policy that satisfies the following equation

$\pi(i)\in\arg\min\left\\{g(i,\pi(i))+\gamma\,\mathbb{E}_{s^{\prime}\sim
P(i,\pi(i),\cdot)}\left[V(s^{\prime})\right]\right\\}\quad\forall\,i\in\mathcal{S}.$
(3)

We use $\text{GreedyPolicy}(V)$ to denote the operator which extracts a greedy
policy associated with $V$.

Equations (1) and (2) admit recursive formulations which are also known as the
Bellman equations. In particular, given $\pi\in\Pi$

$V^{\pi}(i)=g(i,\pi(i))+\gamma\,\mathbb{E}_{s^{\prime}\sim
P(i,\pi(i),\cdot)}\left[V^{\pi}(s^{\prime})\right]\quad\forall\,i\in\mathcal{S}$
(4)

is the Bellman equation associated with policy $\pi$, and

$V^{*}(i)=\min_{\pi\in\Pi}\left\\{g(i,\pi(i))+\gamma\,\mathbb{E}_{s^{\prime}\sim
P(i,\pi(i),\cdot)}\left[V^{*}(s^{\prime})\right]\right\\}\quad\forall\,i\in\mathcal{S}$
(5)

is the Bellman equation associated with the optimal cost. Let
$\mathcal{F}(\mathcal{S},\,\mathbb{R})$ denote the space of functions mapping
from $\mathcal{S}$ to $\mathbb{R}$. Starting from the Bellman equations, we
can define the mappings
$T_{\pi}:\mathcal{F}(\mathcal{S},\,\mathbb{R})\rightarrow\mathcal{F}(\mathcal{S},\,\mathbb{R})$
with

$(T_{\pi}V)(i)=g(i,\pi(i))+\gamma\,\mathbb{E}_{s^{\prime}\sim
P(i,\pi(i),\cdot)}\left[V(s^{\prime})\right]\quad i\in\mathcal{S}\,,$ (6)

and
$T:\mathcal{F}(\mathcal{S},\,\mathbb{R})\rightarrow\mathcal{F}(\mathcal{S},\,\mathbb{R})$
with

$(TV)(i)=\min_{\pi\in\Pi}\left\\{g(i,\pi(i))+\gamma\,\mathbb{E}_{s^{\prime}\sim
P(i,\pi(i),\cdot)}\left[V(s^{\prime})\right]\right\\}\quad i\in\mathcal{S}\,.$
(7)

These mappings are known as the Bellman operators and are very practical for
shorthand notation and algorithmic description (DB_book, , Chapter 1). The
Bellman operators enjoy three fundamental properties: monotonicity,
$\gamma$-contractivity in the infinity-norm and shift-invariance. For
completeness, we report them for the $T$-operator in the following
propositions and we refer to (DB_book, , Chapter 1) for the analogous
properties for the $T_{\pi}$-operator and for the technical proofs.

######  Proposition 1 (Monotonicity).

For any $V:\mathcal{S}\rightarrow\mathbb{R}$ and
$V^{\prime}:\mathcal{S}\rightarrow\mathbb{R}$, and for $k=0,1,\dots$, then

$\left((T^{k}V)(i)-(T^{k}V^{\prime})(i)\right)\,\cdot\,\left(V(i)-V^{\prime}(i)\right)\geq
0\quad i\in\mathcal{S}\,.$ (8)

######  Proposition 2 ($\gamma$-Contractivity).

For any $V:\mathcal{S}\rightarrow\mathbb{R}$ and
$V^{\prime}:\mathcal{S}\rightarrow\mathbb{R}$, and for $k=0,1,\dots$, then

$\max_{i\in\mathcal{S}}\,|(T^{k}V)(i)-(T^{k}V^{\prime})(i)|\leq\gamma^{k}\max_{i\in\mathcal{S}}\,|V(i)-V^{\prime}(i)|\,.$
(9)

######  Proposition 3 (Shift-Invariance).

For any $V:\mathcal{S}\rightarrow\mathbb{R}$, $r\in\mathbb{R}$ and
$k=0,1,\dots$, then

$(T^{k}(V+re))(i)=(T^{k}V)(i)+\gamma^{k}r\quad i\in\mathcal{S}\,,$ (10)

where $e:\mathcal{S}\rightarrow\mathbb{R}$ is the unit function that takes
value 1 identically on $\mathcal{S}$.

From their definitions we evince that the Bellman operators (6) and (7) admit
(1) and (2) as fixed-points, respectively. Finally, uniqueness of their fixed-
points follows directly from the Banach’s fixed-point theorem rockafellar_1976
.

Because we are dealing with finite state and action sets, any function $V$ on
$\mathcal{S}$, including the mappings $T_{\pi}V$ and $TV$, can be represented
as an $n$-dimensional vector. With a slight abuse of notation, this leads to
the following vector-notation

$V=\begin{bmatrix}V(1)\\\ \vdots\\\ V(n)\end{bmatrix}\,,\quad
T_{\pi}V=\begin{bmatrix}(T_{\pi}V)(1)\\\ \vdots\\\
(T_{\pi}V)(n)\end{bmatrix}\,,\quad TV=\begin{bmatrix}(TV)(1)\\\ \vdots\\\
(TV)(n)\end{bmatrix}\,,$ (11)

which allows for a compact formulation of the Bellman equations, i.e.,
$V^{\pi}=T_{\pi}V^{\pi}$ and $V^{*}=TV^{*}$. Finally, given a policy
$\pi\in\Pi$, we can represent its associated transition probability
distributions and stage-costs in matrix and vector-form, respectively. In
particular, we denote with $P^{\pi}\in\mathbb{R}^{n\times n}$ its associated
transition probability matrix, where $p^{\pi}_{ij}=P^{\pi}(i,\pi(i),j)$; and
with $g^{\pi}\in\mathbb{R}^{n}$ its associated cost-vector, where
$g^{\pi}_{i}=g^{\pi}(i,\pi(i))$. Notice that, for any $\pi\in\Pi$,
$P^{\pi}\in\mathbb{S}^{n\times n}$.

Row-stochastic matrices are an important subset of non-negative matrices
(nonnegative_matrices, , Chapter 2), i.e., square matrices all of whose
elements are non-negative. We now recall some fundamental definitions and
properties of non-negative matrices that will be later used to produce a
classification of discounted MDPs with finite spaces. We refer to
(nonnegative_matrices, , Chapter 2) for a more thorough review on non-negative
matrices as well as for the technical proofs of the following propositions.

######  Proposition 4.

If $A$ is non-negative, then $(\rho(A),0)$ is an eigenvalue.

We refer to the proof of statement (a) in Theorem 1.1 in
(nonnegative_matrices, , Chapter 2) for a proof of Proposition 4.

Since $\rho(A)=1$ when $A$ is stochastic, from the previous proposition we can
deduce that stochastic matrices always have $(1,0)$ as an eigenvalue. In
addition, by the Gershgorin circle theorem gerschgorin_31 , the spectrum of
stochastic matrices is contained in the circle of radius $1$ and center
$(0,0)$ in the complex plane. Combining these two properties, we obtain that
$\lambda_{1}=(1,0)$ (see Theorem 5.3 in (nonnegative_matrices, , Chapter 2)
for a formal proof).

The following definition characterizes the subset of irreducible matrices.

######  Definition 1 (Irreducible Non-Negative Matrices).

A non-negative matrix $A$ is irreducible if and only if for every $(i,j)$
there exists a natural number $q$ such that $a^{(q)}_{ij}>0$.

As we can infer from the following proposition, irreducibility of a non-
negative matrix has important implications on the distribution of its
eigenvalues.

######  Proposition 5.

If $A$ is irreducible, then $\rho(A)$ is a simple eigenvalue and any other
eigenvalue of $A$ with modulus $\rho(A)$ is also simple.

We refer to the proof of statement (b) in Theorem 1.4 in
(nonnegative_matrices, , Chapter 2) for a proof of Proposition 6.

Within the set of irreducible matrices it is possible to define the subset of
primitive matrices, which verify a more stringent condition than the one in
Definition 1. We formalize this subset with the following definition.

######  Definition 2 (Primitive Non-Negative Matrices).

An irreducible matrix $A$ is primitive if and only if there exists a natural
number $n_{p}$ such that $A^{n_{p}}$ is positive.

From the previous definition, it is clear that every primitive matrix is also
irreducible, but not viceversa. For instance, it is easy to verify that the
following matrix is irreducible but not primitive

$A=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}\,.$ (12)

The following proposition characterizes primitiveness of an irreducible matrix
via a necessary and sufficient condition on its spectrum.

######  Proposition 6.

$A$ is primitive if and only if $\rho(A)$ is greater in magnitude than any
other eigenvalue.

We refer to the proof of statement (a) in Theorem 1.7 in
(nonnegative_matrices, , Chapter 2) for a proof of Proposition 6.

Irreducibility and primitiveness have an intuitive graph theoretical
interpretation. In particular, denoting with $G(A)$ the graph associated with
the non-negative matrix $A$, we can conclude that $A$ is irreducible if and
only if $G(A)$ is strongly connected, i.e., there exists a finite sequence of
edges which connects any two nodes in the graph. Primitiveness instead
translates into the stronger requirement of reachability of any node from any
other node after a sequence of exactly $n_{p}$ edges.

We conclude this brief recap on non-negative matrices and some of their
fundamental properties with the following definition, which formalizes the
index of cyclicity of an irreducible matrix.

######  Definition 3 (Index of Cyclicity).

Let $A$ be an irreducible matrix. We call index of cyclicity of $A$, and we
denote it with $h(A)$, the number of eigenvalues of $A$ with modulus equal to
$\rho(A)$.

Going back to (12) as example, $h(A)=2$ since
$\Lambda(A)=\left\\{+1,\,-1\right\\}$. In addition, notice that, because of
Proposition 6, $h(A)=1$ if and only if $A$ is primitive.

As mentioned, row-stochastic matrices are an important subset of non-negative
matrices. We can therefore deploy the structural properties of the latter to
further characterize discounted MDPs with finite spaces. In particular, we
refer to any discounted MDP with finite spaces as general and distinguish two
sub-classes, ergodic and regular discounted MDPs with finite spaces, based on
the structural properties of $P^{\pi}$ for all $\pi\in\Pi$. We formalize this
classification with the following definition.

######  Definition 4 (Ergodic & Regular MDPs).

A discounted MDP with finite spaces is called ergodic if and only if the
transition probability matrix corresponding to every policy is irreducible.
If, in addition, the transition probability matrix corresponding to every
policy is primitive, then it is called regular.

A similar principle for the classification of finite MDPs is adopted by
Puterman in (puterman_mdp, , Chapter 8).

By relying on the discussed properties of irreducible matrices, we can
conclude that, for ergodic MDPs, $P^{\pi}$ has $h(P^{\pi})\geq 1$ simple
eigenvalues with modulus 1 for all $\pi\in\Pi$, while all the other
eigenvalues are contained strictly inside the unitary circle in the complex
plane. In addition, $h(P^{\pi})=1$ for all $\pi\in\Pi$ if and only if the MDP
is regular. To keep the notation compact, in the following we will use
$h_{\pi}$ in place of $h(P^{\pi})$.

###  2.1 Dynamic Programming

Dynamic programming (DP) comprises the methods to solve the Bellman equations.
In particular, from now on we refer to the solution of (4) for a given policy
$\pi\in\Pi$ as policy evaluation. In general terms, we can recognize three
fundamental DP methods: value iteration (VI), policy iteration (PI) and the
linear programming approach (DB_book, , Chapter 2). In this work we only focus
on the first two methods and some of their variants.

VI, also known as the successive approximation method, is an iterative method
to solve the Bellman equations. In particular, VI is based on repeated
applications of the Bellman operators starting from an arbitrary cost-vector
$V_{0}\in\mathbb{R}^{n}$, where the operator in (6) is used to solve (4) and
the operator in (7) to solve (5). In the first case, we refer to it as VI for
policy evaluation. Algorithm 1 provides an algorithmic description of VI. The
following proposition characterizes the convergence properties of VI and VI
for policy evaluation when the infinity-norm is considered to measure the
distance of the generated iterates from the solution.

Algorithm 1 Value Iteration

1:$V_{0}\in\mathbb{R}^{n},\,\textit{tol}>0$

2:$k\leftarrow 0$

3:$V_{1}\leftarrow TV_{0}$

4:while $\|V_{k+1}-V_{k}\|_{\infty}>tol$ do

5: $k\leftarrow k+1$

6: $V_{k+1}\leftarrow TV_{k}$

7:end while

######  Proposition 7.

Consider a generic infinite-horizon discounted MDP with finite spaces and any
policy $\pi\in\Pi$ and let $\left\\{V_{k}\right\\}$ and
$\left\\{V^{\pi}_{k}\right\\}$ denote the sequence generated by VI and VI for
policy evaluation, respectively. Then for any $k\geq 0$

$\|V_{k+1}-V^{*}\|_{\infty}\leq\gamma\|V_{k}-V^{*}\|_{\infty}$ (13)
$\|V^{\pi}_{k+1}-V^{\pi}\|_{\infty}\leq\gamma\|V^{\pi}_{k}-V^{\pi}\|_{\infty}$
(14)

Proof: We report the proof for VI. Analogous steps can be done to prove the
result for VI for policy evaluation.

We consider an arbitrary starting point $V_{0}\in\mathbb{R}^{n}$ and we prove
the final result by upper-bounding the infinity-norm of the distance between
the $(k+1)$-iterate and the solution as follow

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}$
$\displaystyle=\|TV_{k}-TV^{*}\|_{\infty}$ (15)
$\displaystyle\leq\gamma\,\|V_{k}-V^{*}\|_{\infty}\quad k=0,1,\dots\,,$

where the equality is obtained by using the definition of VI iterate and the
fact that the $T$-operator has $V^{*}$ as unique fixed-point, and the
inequality follows from Proposition 2. Extensive numerical evidence shows that
the results of Proposition 7 are tight, which translates in an exacerbating
slow convergence for values of $\gamma$ close to 1. By using monotonic error
bounds in some cases it is possible to improve the convergence rate of VI to
$\gamma\,|\lambda_{2}|$, where $\lambda_{2}$ denotes the subdominant
eigenvalue of $P^{\pi^{*}}$ (DB_book, , Chapter 2). Unfortunately, there are
often situations where $|\lambda_{2}|=1$, e.g., ergodic but not regular MDPs
with $h(P^{\pi^{*}})>1$, which results in VI with monotonic error bounds
having the same contraction as the original VI.

On a final note, VI only produces a sequence of cost-iterates, but it is
possible to produce a sequence of policies too by extracting at every
iteration a greedy policy associated with the current cost-iterate.

Equation (5) can also be solved with PI, an iterative method that alternates
two steps: policy evaluation and policy improvement. In particular, we start
by extracting a greedy policy associated with an arbitrary cost
$\tilde{V}\in\mathbb{R}^{n}$ and, until convergence, we first evaluate the
cost associated with the current policy (policy evaluation step) and use this
value to update the cost-iterate, and then we update the current policy with a
greedy policy associated with the current cost-iterate (policy improvement
step). See Algorithm 2 for an algorithmic description of PI. Even if PI with
exact arithmetic is guaranteed to converge in a finite-number of iterations,
the upper bound on the number of iterations is exponential in the number of
states and it is therefore important to characterize its improvement-per-
iteration. The following proposition characterizes the convergence properties
of PI.

######  Proposition 8.

Consider a generic infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ denote the sequence of iterates generated by PI. Then
for any $k\geq 0$

$\|V_{k+1}-V^{*}\|_{\infty}\leq\gamma\|V_{k}-V^{*}\|_{\infty}\,.$ (16)

In addition, PI with exact arithmetic converges to the optimal policy in at
most $m^{n}$ iterations.

Proof: We start by proving the second statement. By exploiting the definition
of greedy policy and the properties of the Bellman operators, it is possible
to show that $V_{k+1}\leq V_{k}$ for $k=0,1,\dots$. For a step-by-step proof
we refer to (DB_book, , Proposition 2.3.1, Chapter 2). The final result
follows directly by considering that there only exists $m^{n}$ policies.

Since $V_{k+1}\leq V_{k}$ for $k=0,1,\dots$, by exploiting the monotonicity of
$T_{\pi_{k+1}}$ and the fact that $V_{k+1}$ is its unique fixed-point, we
obtain the following upper-bound on $V_{k+1}$

$\displaystyle V_{k+1}$ $\displaystyle=T_{\pi_{k+1}}V_{k+1}$ (17)
$\displaystyle\leq T_{\pi_{k+1}}V_{k}$ $\displaystyle=TV_{k}\,,$

where the last equality follows from the fact that $\pi_{k+1}$ is a greedy
policy associated to $V_{k}$, therefore $T_{\pi_{k+1}}V_{k}=TV_{k}$.
Considering the previous derivations and since, by definition, $V^{*}\leq
V^{\pi}$ for any $\pi\in\Pi$, we obtain that

$0\leq V_{k+1}-V^{*}\leq TV_{k}-TV^{*}\,.$ (18)

Finally, by taking the infinity-norm on both sides of (18) and considering
that $V^{*}=TV^{*}$, we obtain the final result

$\|V_{k+1}-V^{*}\|_{\infty}\leq\|TV_{k}-TV^{*}\|_{\infty}\leq\gamma\,\|V_{k}-V^{*}\|_{\infty}\,,$
(19)

where the last inequality follows from Proposition 2.

Algorithm 2 Policy Iteration

1:$\tilde{V}\in\mathbb{R}^{n},\,\textit{tol}>0$

2:$k\leftarrow 0$

3:$\pi_{0}\leftarrow\text{GreedyPolicy}(\tilde{V})$

4:$V_{0}\leftarrow V^{\pi_{0}}$

5:$\pi_{1}\leftarrow\text{GreedyPolicy}(V_{0})$

6:$V_{1}\leftarrow V^{\pi_{1}}$

7:while $\|V_{k+1}-V_{k}\|_{\infty}>tol$ do

8: $k\leftarrow k+1$

9: $\pi_{k+1}\leftarrow\text{GreedyPolicy}(V_{k})$

10: $V_{k+1}\leftarrow V^{\pi_{k+1}}$

11:end while

Despite the results in Proposition 8, Figure 1 and extensive empirical
observations suggest that PI enjoys a qualitatively better rate of convergence
than VI. These observations are also confirmed by the Newton-based theoretical
analysis of the local contraction properties of DP methods, which is discussed
in Section 3. This leads to conclude that, at least locally, the results of
Proposition 8 are loose and not able to describe properly the practical
convergence behaviour of PI-iterates.

Regarding the computational complextiy of PI, the policy evaluation (step 3, 5
and 9 in Algorithm 2) comprises the exact solution of an $n$-dimensional
linear system. In particular, given a policy $\pi\in\Pi$, $V^{\pi}$ is
obtained by solving exactly the following linear system

$\left(I-\gamma P^{\pi}\right)V^{\pi}=g^{\pi}\,.$ (20)

In general, when the number of states is large, the exact solution of (20) is
computationally expensive, e.g., $\frac{2}{3}n^{3}+\mathcal{O}(n^{2})$ with
LU-factorization. Despite its fast rate of convergence, this results in an
overall poor computational scalability of PI, as also confirmed by the results
reported in Table 1.

Figure 1: We consider an infinite-horizon discounted MDP with $\gamma=0.9$,
$n=10000$ and $m=40$ and we plot the distance in infinity-norm of the iterates
from the solution versus number of iterations for policy iteration, value
iteration and optimistic policy iteration with different values of inner
iterations.

A better trade-off between computational complexity and convergence rate is
achieved by optimistic policy iteration (OPI) (DB_book, , Chapter 2), also
known in the literature as modified policy iteration (puterman_mdp, , Chapter
7). This variant of PI consists in substituting the policy evaluation step
with the inexact solution of (20) obtained with a finite number
$w\in\mathbb{N}$ of iterations of VI for policy evaluation starting from the
current iterate. This results in the following cost-iterates

$V_{k+1}=T_{\pi_{k+1}}^{w}V_{k}\quad k=0,1,\dots\,,$ (21)

where $\pi_{k+1}=\text{GreedyPolicy}(V_{k})$ and $V_{0}\in\mathbb{R}^{n}$ is
the initial iterate. Clearly, when $w=1$ we recover VI and when $w=\infty$ we
recover PI. In practice, this value is adjusted heuristically in order to
produce the best trade-off between computational complexity and convergence
rate for the specific problem at hand. Indeed, the more inner iterations are
run and the faster is the resulting convergence in terms of iterations (see
Figure 1), but the more computationally expensive is the overall method. The
convergence properties of OPI are discussed in the following proposition.

######  Proposition 9.

Consider a general infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ denote the sequence of iterates generated by OPI with
$V_{0}\in\mathbb{R}^{n}$. Then for any $k\geq 0$

$\|V_{k}-V^{*}\|_{\infty}\leq\gamma^{k}\|V_{0}-V^{*}\|_{\infty}+\gamma^{k}b\,.$
(22)

In addition, if $TV_{0}\leq V_{0}$, then for any $k\geq 0$

$\|V_{k+1}-V^{*}\|_{\infty}\leq\gamma\|V_{k}-V^{*}\|_{\infty}\,.$ (23)

Proof: We start by defining the auxiliary bar-iterates, where
$\bar{V}_{0}=V_{0}+be$ and $b\in\mathbb{R}_{+}$. In particular, we set $b=0$
if $TV_{0}\leq V_{0}$ and
$b=\frac{1}{(1-\gamma)}\,\max_{i\in\mathcal{S}}\,(TV_{0}-V_{0})(i)$ otherwise.
From the definition of $b$ it follows directly that
$T\bar{V}_{0}\leq\bar{V}_{0}$. Indeed, if $b=0$ then the bar-quantities
coincide with the original iterates and $TV_{0}\leq V_{0}$ by assumption. If
$b>0$, then $TV_{0}-V_{0}\leq(1-\gamma)be$. By adding $\gamma be$ on both
sides of the previous inequality and by exploiting Proposition 3 and the
definition of the bar-quantities, we obtain that $TV_{0}+\gamma
be=T\bar{V}_{0}\leq V_{0}+be=\bar{V}_{0}$. Since for any $k\geq 0$
$\bar{V}_{k}=V_{k}+\gamma^{kw}be$, the bar-quantities produce the same
sequence of greedy-policies $\left\\{\pi_{k}\right\\}$ as the original
iterates. We can therefore proceed with studying the convergence in terms of
the bar-iterates.

By following the same steps as in the proof of Proposition 2.3.2 in (DB_book,
, Chapter 2), we are able to conclude that the following inequalities hold for
any $k\geq 0$

$\bar{V}_{k+1}\leq T\bar{V}_{k}$ (24) $V^{*}\leq\bar{V}_{k}\leq
T^{k}\bar{V}_{0}\,.$ (25)

By using the fact that $\bar{V}_{k}=V_{k}+\gamma^{kw}be$ and re-arranging the
terms in (25), we obtain that for any $k\geq 0$ the following inequalities
hold

$\displaystyle-\gamma^{kw}be$ $\displaystyle\leq V_{k}-V^{*}\leq
T^{k}V_{0}-V^{*}+\gamma^{k}be-\gamma^{kw}be\,.$ (26)

Since $\gamma^{k}be\geq\gamma^{kw}be\geq 0$ for any $w\in\mathbb{N}$ and by
exploiting the fact that $V^{*}$ is the unique fixed-point of the
$T$-operator, we obtain that the following inequalities hold for any $k\geq 0$

$\displaystyle-\gamma^{k}\|V_{0}-V^{*}\|_{\infty}e-\gamma^{k}be$
$\displaystyle\leq V_{k}-V^{*}\leq T^{k}V_{0}-T^{k}V^{*}+\gamma^{k}be\,.$ (27)

Starting from (27) and exploiting the fact that
$T^{k}V_{0}-T^{k}V^{*}\leq\|T^{k}V_{0}-T^{k}V^{*}\|_{\infty}e\leq\gamma^{k}\|V_{k}-V^{*}\|_{\infty}e$,
where the last inequality follows from Proposition 2, and since
$\|V_{0}-V^{*}\|_{\infty}\geq 0$, we obtain that the following inequalities
hold for any $k\geq 0$

$\displaystyle-\gamma^{k}\|V_{0}-V^{*}\|_{\infty}e-\gamma^{k}be$
$\displaystyle\leq
V_{k}-V^{*}\leq\gamma^{k}\|V_{0}-V^{*}\|_{\infty}e+\gamma^{k}be\,.$ (28)

From (28), we can conclude the final result

$\|V_{k}-V^{*}\|_{\infty}\leq\gamma^{k}\|V_{0}-V^{*}\|_{\infty}+\gamma^{k}b\,.$
(29)

In addition, if $b=0$, then, from (24) and (25), we can conclude that for any
$k\geq 0$ the following holds

$0\leq V_{k+1}-V^{*}\leq TV_{k}-V^{*}=TV_{k}-TV^{*}\,,$ (30)

where the equality follows from the fact that $V^{*}$ is the unique fixed-
point of the $T$-operator. By switching to the infinity-norm and exploiting
Proposition 2, we obtain the final result

$\|V_{k+1}-V^{*}\|_{\infty}\leq\gamma\|V_{k}-V^{*}\|_{\infty}\,.$ (31)

Results on the convergence of OPI can also be found in (DB_book, , Proposition
2.3.2). Our results give also an upper-bound on the improvement-per-iteration.
Unfortunately, Proposition 9 as well as Proposition 2.3.2 in DB_book  fail in
capturing the improvement in terms of convergence rate that is observed when
increasing the number of inner iterations $w$, and therefore the level of
accuracy that is used to solve the policy evaluation tasks.

Table 1: We consider an infinite-horizon discounted MDPs with $\gamma=0.9$, $m=40$ and different sizes of the state space and we report the total number of iterations to reach an accuracy of $10^{-9}$ and the cost-per-iteration in terms of CPU time for policy iteration and value iteration. The benchmarks are run on Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz enforcing single core execution. | Policy Iteration | Value Iteration
---|---|---
$\gamma=0.9$ | Iterations | CPU Time | Iterations | CPU Time
$n=100$ | 6 | $3.67\times 10^{-4}$ | 270 | $1.96\times 10^{-4}$
$n=1000$ | 6 | $4.55\times 10^{-2}$ | 271 | $1.81\times 10^{-2}$
$n=10000$ | 7 | $17.07$ | 273 | $1.84$
|  |  |  |

##  3 The Bellman Residual Function

Solving the Bellman equation (5) can be recast as a nonlinear non-smooth root-
finding problem, where the residual function is the so-called Bellman residual
function $r:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, with $r(V)=V-TV$
gargiani_2022 . Because of the relation between $T$ and $T_{\pi}$, we can also
rewrite the Bellman residual function as follows

$r(V)=V-\min_{\pi\in\Pi}\left\\{T_{\pi}V\right\\}=V-\min_{\pi\in\Pi}\left\\{g^{\pi}+\gamma
P^{\pi}V\right\\}\,.$ (32)

The function in (32) is piecewise affine with at most $|\Pi|$ selection
functions. In particular, we denote with $\tilde{\Pi}\subseteq\Pi$ the non-
empty subset of policies that are associated to the selection functions of
$r$, i.e. for any $\pi\in\tilde{\Pi}$ there exists $V\in\mathbb{R}^{n}$ such
that $r(V)=V-T_{\pi}V$, and we use $r^{\pi}(V)=V-T_{\pi}V$ for any
$\pi\in\Pi$. For a general piecewise affine function, it is always possible to
add extra spurious selection functions to the minimal set of effective ones
(facchinei_bookI, , Chapter 4). The following assumption excludes the presence
of spurious selection functions by requiring all the policies in $\tilde{\Pi}$
to be associated to effective selection functions.

######  Assumption 1.

Consider a discounted MDP with finite spaces. Let $\tilde{\Pi}\subseteq\Pi$
denote the non-empty subset of policies associated to the selection functions
of the Bellman residual function. We assume that for any $\pi\in\tilde{\Pi}$
then

$\text{int}\left(\\{V\in\mathbb{R}^{n}\,\,\,\text{s.t.}\,\,\,r(V)=V-T_{\pi}V\\}\right)\neq\emptyset\,.$
(33)

Under Assumption 1, the generalized Jacobian of (32) at $V\in\mathbb{R}^{n}$
is given by the following set

$\partial r(V)=\text{Conv}\left(\\{I-\gamma P^{\pi},\,\text{where
}\pi\in\text{GreedyPolicy}(V)\\}\right)\,.$ (34)

Since for any $V\in\mathbb{R}^{n}$ every $J\in\partial r(V)$ is non-singular,
(32) is globally CD-regular (gargiani_2022, , Proposition 2). In addition, as
shown in (gargiani_2023, , Lemma 3), for any $\pi\in\Pi$

$\|(I-\gamma P^{\pi})^{-1}\|_{\infty}\leq 1/(1-\gamma)\,.$ (35)

Because of its piecewise affine structure, the Bellman residual function is
globally $(1+\gamma)$-Lipschitz continuous in the infinity-norm
(gargiani_2023, , Lemma 2) and strongly semismooth everywhere
(facchinei_bookII, , Proposition 7.4.7). In particular, we are interested in
strong semismoothness at the root $V^{*}$.

######  Lemma 10.

Consider a general infinite-horizon discounted MDP with finite spaces and let
$r:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ denote its associated Bellman
residual function as defined in (32). Then there exists $\delta>0$ such that,
for any $C\geq\frac{2R}{(1-\gamma)\delta^{2}}$, the following inequality holds
for all $V\in\mathbb{R}^{n}$

$\|r(V)-r(V^{*})-\left(I-\gamma
P^{\pi}\right)\left(V-V^{*}\right)\|_{\infty}\leq
C\,\|V-V^{*}\|^{2}_{\infty}\,,$ (36)

where $\pi\in\text{GreedyPolicy}(V)$. We call
$C_{\min}=\frac{2R}{(1-\gamma)\delta^{2}}$ the global strong semismoothness
constant of $r$ at $V^{*}$.

Proof: We start the proof by upper-bounding the LHS in (36)

$\displaystyle\|r(V)-r(V^{*})-(I-\gamma P^{\pi})(V-V^{*})\|_{\infty}$
$\displaystyle=\|(I-\gamma P^{\pi})V^{*}-g^{\pi}\|_{\infty}$ (37)
$\displaystyle\leq(1+\gamma)\|V^{*}\|_{\infty}+R$
$\displaystyle\leq\frac{2}{1-\gamma}R\,,$

where the last inequality follows from the definition of $V^{*}$ and the limit
of the geometric series. We start by considering the case where $V$ is “close-
enough” to $V^{*}$. In particular, there exists $\delta>0$ such that, for all
$V\in\mathcal{B}(V^{*},\delta)$, then $\pi$ is an optimal policy, since by
continuity the selection functions that are active at points in this
neighborhood are a subset of the selection functions active at the root (a
similar argument is also used in the proof of Theorem 7.2.15 in
facchinei_bookII ). Consequently, whenever $V\in\mathcal{B}(V^{*},\delta)$,
then the LHS in (37) is equal to zero and (36) trivially holds for any $C\geq
0$. Otherwise, we can always select $C\geq 2R/((1-\gamma)\delta^{2})$ and (36)
is verified since $\|V-V^{*}\|_{\infty}\geq\delta$.

###  3.1 Semismooth Newton-Type Methods for DP

The semismooth adaptation of Newton’s method and its variants constitute the
main class of algorithms to solve general nonlinear non-smooth root finding
problems. If the residual function enjoys certain structural properties, i.e.,
Lipschitz continuity and strong semismoothness in a neighborhood of the root,
the semismooth variant of Newton’s method maintains the fast local quadratic
convergence typical of Newton’s method (izmailov_book, , Chapter 2).

Instead of solving the original problem directly, semismooth Newton method
starts from an initial guess $V_{0}\in\mathbb{R}^{n}$ and, for every $k\geq
0$, solves the following linear approximation of (32) at the current iterate
$V_{k}$

$r(V_{k})+J_{k}(V_{k+1}-V_{k})=0\,,$ (38)

where $J_{k}\in\partial r(V_{k})$. The iterate is then updated with the
derived solution of (38) as follows

$V_{k+1}=V_{k}-J_{k}^{-1}r(V_{k})\,.$ (39)

As proved in gargiani_2022 , under Assumption 1 policy iteration is an
instance of semismooth Newton method. By selecting $J_{k}=I-\gamma
P^{\pi_{k+1}}$, with $\pi_{k+1}=\text{GreedyPolicy}(V_{k})$, and plugging
these quantities back into (39), we obtain the following expression for the
update of semismooth Newton method

$V_{k+1}=(I-\gamma P^{\pi_{k+1}})^{-1}g^{\pi_{k+1}}\,,$ (40)

which is equivalent to the policy iteration cost update. In light of this
connection and because of the structural properties of the Bellman residual
function, it is possible to conclude local quadratic contraction for policy
iteration (gargiani_2022, , Proposition 3.5).

The semismooth Newton-type variants are obtained by replacing
$J_{k}\in\partial r(V_{k})$ in (38) with a non-singular matrix
$B_{k}\in\mathbb{R}^{n\times n}$. In our setting, the main reason behind this
approximation of the original method are the computational costs. Indeed,
whenever $B_{k}$ is selected with a favorable sparsity structure,
factorization and inversion operations are generally less expensive, leading
to lower computational costs for the iteration update. In addition, whenever
$B_{k}$ approximates “well-enough” an element in $\partial r(V_{k})$ for all
$k\geq 0$, then the method preserves local linear contraction guarantees
gargiani_2022 . As shown in gargiani_2022 , VI can be interpreted as a
semismooth Newton-type method, where $B_{k}=I$ for all $k\geq 0$. A simple
variant of VI, called $\alpha$-VI, is obtained by approximating
$J_{k}\in\partial r(V_{k})$ with a scaled identity matrix for all $k\geq 0$,
where $\alpha>0$ is the scaling parameter and plays a fundamental role in
determining the convergence properties of the method gargiani_2022 .

The following proposition summarizes the local convergence properties of
semismooth Newton-type methods when used to find the root of the Bellman
residual function.

######  Theorem 11.

Consider a general infinite-horizon discounted MDP with finite spaces and let
$r:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ denote its associated Bellman
residual function as defined in (32). In addition, let
$V_{0}\in\mathbb{R}^{n}$, $L>0$ and $\kappa\in[0,1)$ be a constant. For any
sequence of nonsingular matrices
$\left\\{B_{k}\right\\}\subseteq\mathbb{R}^{n\times n}$, such that, for all
$k\geq 0$, $\|B_{k}^{-1}\|_{\infty}\leq L$ and $\exists\,\,J_{k}\in\partial
r(V_{k})$ for which the kappa-condition

$\|B_{k}^{-1}(B_{k}-J_{k})\|_{\infty}\leq\kappa_{k}\leq\kappa$ (41)

is verified, the sequence $\left\\{V_{k}\right\\}\subseteq\mathbb{R}^{n}$,
with

$V_{k+1}=V_{k}-B_{k}^{-1}r(V_{k})\,,$ (42)

verifies the following inequality

$\|V_{k+1}-V^{*}\|_{\infty}\leq\kappa_{k}\|V_{k}-V^{*}\|_{\infty}+LC_{\min}\,\|V_{k}-V^{*}\|^{2}_{\infty},$
(43)

where $C_{\min}$ is the global strong semismooth constant of $r$ at $V^{*}$.
Therefore, if
$V_{0}\in\mathcal{B}\left(V^{*},\frac{1-\kappa}{LC_{\min}}\right)$, the
sequence of iterates generated by (42) is Q-linearly converging to $V^{*}$ in
the infinity-norm with rate at least
$\kappa+LC_{\min}\,\|V_{0}-V^{*}\|_{\infty}$.

Proof: A proof of (43) can be derived by following similar steps as in
(gargiani_extended_2022, , Theorem 2.6).

We derive the first result by upper-bounding the infinity-norm of the
difference between the $(k+1)$-iterate and the root $V^{*}$ as follow

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}=\|V_{k}-B_{k}^{-1}r(V_{k})-V^{*}\|_{\infty}$
$\displaystyle=\|B_{k}^{-1}\left(B_{k}-J_{k}\right)\left(V_{k}-V^{*}\right)-B_{k}^{-1}\left(r(V_{k})-r(V^{*})-J_{k}\left(V_{k}-V^{*}\right)\right)\|_{\infty}$
$\displaystyle\overset{(a)}{\leq}\|B_{k}^{-1}\left(B_{k}-J_{k}\right)\|_{\infty}\|V_{k}-V^{*}\|_{\infty}+\|B_{k}^{-1}\|_{\infty}\|r(V_{k})-r(V^{*})-J_{k}\left(V_{k}-V^{*}\right)\|_{\infty}$
$\displaystyle\overset{(b)}{\leq}\kappa_{k}\|V_{k}-V^{*}\|_{\infty}+LC_{\min}\,\|V_{k}-V^{*}\|^{2}_{\infty}\,,$

where inequality $(a)$ follows from the fundamental properties of the
infinity-norm and inequality $(b)$ from the assumptions on the sequence
$\left\\{B_{k}\right\\}$ and Lemma 10. From the previous derivations we can
conclude that

$\|V_{k+1}-V^{*}\|_{\infty}\leq\left(\kappa_{k}+LC_{\min}\|V_{k}-V^{*}\|_{\infty}\right)\|V_{k}-V^{*}\|_{\infty}\,.$
(44)

Finally, we can use (44) to obtain a sufficient condition on the initial
iterate $V_{0}$ such that the sequence $\left\\{V_{k}\right\\}$ enjoys local
Q-linear convergence to $V^{*}$. In particular, starting from (44), we obtain
that if

$\|V_{0}-V^{*}\|_{\infty}<\frac{1-\kappa}{LC_{\min}}\,,$ (45)

then the iterates are Q-linearly contracting to $V^{*}$ in the infinity-norm
with rate at least $\kappa+LC_{\min}\,\|V_{0}-V^{*}\|_{\infty}$.

Notice that for exact semismooth Newton method for all $k\geq 0$
$\kappa_{k}=0$ since $B_{k}\in\partial r(V_{k})$. Consequently, the generated
sequence of iterates enjoys local Q-quadratic convergence to $V^{*}$ in the
infinity-norm as (43) reduces to

$\|V_{k+1}-V^{*}\|_{\infty}\leq LC_{\min}\,\|V_{k}-V^{*}\|^{2}_{\infty}\,.$
(46)

This consideration, together with the fact that, under Assumption 1, PI is an
instance of semismooth Newton method, leads to the results discussed in the
following corollary.

######  Corollary 12.

Consider a general infinite-horizon discounted MDP with finite spaces. Under
Assumption 1, the sequence $\left\\{V_{k}\right\\}\subseteq\mathbb{R}^{n}$
generated by PI enjoys local Q-quadratic convergence to $V^{*}$.

##  4 Inexact Policy Iteration Methods

Even if, as proved in Theorem 11, the semismooth Newton method enjoys a fast
rate of convergence, computing the exact solution of (38) using a direct
method can be expensive when the number of unknowns is large. A more
computationally efficient solution in the large-scale case consists in solving
(38) only approximately with some iterative solver for linear systems and
using a certain stopping rule. These are the principles behind the inexact
variants of semismooth Newton method MARTINEZ1995127 . In particular,
$V_{k+1}$ is no longer required to exactly solve (38), but only to satisfy

$\|r(V_{k})+J_{k}(V_{k+1}-V_{k})\|\leq\alpha_{k}\|r(V_{k})\|\,,$ (47)

for some $\alpha_{k}\in[0,1)$. The sequence $\left\\{\alpha_{k}\right\\}$ is
called forcing sequence and it greatly affects both local convergence
properties and robustness of the method MARTINEZ1995127 . Different iterative
solvers for linear systems can be used to approximately solve (38). Often
Krylov subspace methods, such as the generalized minimal residual method
(GMRES) gmres , are deployed in large-scale scenarios. Based on these
observations, we define a novel variant of PI for large-scale scenarios, which
we call inexact policy iteration (iPI) methods. This class of methods is based
on approximately solving the policy evaluation step with an iterative solver.
The methods in this class start with an initial guess of the optimal cost
$V_{0}\in\mathbb{R}^{n}$ and then at every iteration extract a greedy policy
associated with the current iterate $V_{k}\in\mathbb{R}^{n}$, which is used to
compute an element in Clarke’s generalized Jacobian. The next iterate
$V_{k+1}\in\mathbb{R}^{n}$ is selected as an approximate solution of the
Newtonian linear system

$(I-\gamma P^{\pi_{k+1}})V=g^{\pi_{k+1}}$ (48)

which verifies the stopping condition in (47) with the infinity-norm. Because
of the specific structure of the Bellman residual function, (47) simplifies to

$\|g^{\pi_{k+1}}-(I-\gamma
P^{\pi_{k+1}})V_{k+1}\|\leq\alpha_{k}\|g^{\pi_{k+1}}-(I-\gamma
P^{\pi_{k+1}})V_{k}\|\,.$ (49)

In principle, any iterative solver for linear systems with non-singular
coefficient matrices can be used to generate an approximate solution of (48),
such as VI, its mini-batch version gargiani_minibatch  and GMRES. Notice that,
when VI is deployed as inner solver, we obtain a variant of OPI where the
number of inner iterations is not selected a priori, but dictated by the
stopping condition. The parametric stopping condition allows for a dynamic
selection of the number of inner iterations based on locally available
information. Finally, as we discuss in Section 4.2 and 5.2, deploying a
different inner solver than VI for policy evaluation may lead to a superior
convergence performance.

See Algorithm 3 for a pseudocode description of a general iPI method with a
constant $\alpha$-parameter and where $\text{IterativeSolver}(A,b,\tilde{x})$
is used to denote the operator that applies one iteration with a well-defined
iterative solver for linear systems with non-singular coefficient matrices to
the linear system $Ax=b$ with $A$ being non-singular and starting from
$\tilde{x}$.

Algorithm 3 Inexact Policy Iteration

1:$V_{0}\in\mathbb{R}^{n},\,\textit{tol}>0,\alpha\in(0,1)$

2:$k\leftarrow 0$

3:while $\|r(V_{k})\|_{\infty}>tol$ do

4: $\pi_{k+1}\leftarrow\text{GreedyPolicy}(V_{k})$

5: $J_{k}\leftarrow I-\gamma P^{\pi_{k+1}}$

6: $\theta_{0}\leftarrow V_{k}$

7: $i\leftarrow 0$

8: while
$\|g^{\pi_{k+1}}-J_{k}{\theta}_{i}\|_{\infty}\leq\alpha\|g^{\pi_{k+1}}-J_{k}V_{k}\|_{\infty}$
do

9:
${\theta}_{i+1}\leftarrow\text{IterativeSolver}(J_{k},g^{\pi_{k+1}},{\theta}_{i})$

10: $i\leftarrow i+1$

11: end while

12: $V_{k+1}\leftarrow\theta_{i}$

13: $k\leftarrow k+1$

14:end while

###  4.1 Local and Global Convergence

This section is dedicated to the analysis of the local and global convergence
properties of inexact policy iteration methods. Local results were originally
studied in gargiani_2023 , but we will report them for completeness. On that
side, we complement the analysis by adding an estimate for the region of
attraction. The latter is later exploited in conjunction with additional
results to show global convergence of a general iPI method.

######  Theorem 13 (local convergence).

Consider a general infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ be the sequence of iterates generated by a general
inexact policy iteration method as described in Algorithm 3 with
$V_{0}\in\mathbb{R}^{n}$. In addition, let
$\omega=\frac{(1-\gamma)-(1+\gamma)\alpha}{C_{\min}}$. Under Assumption 1, if
$\left\\{\alpha_{k}\right\\}\subseteq[0,\,\alpha]$ with
$\alpha\in(0,\frac{1-\gamma}{1+\gamma})$ and
$V_{0}\in\mathcal{B}(V^{*},\omega)$, then
$\left\\{\|V_{k}-V^{*}\|_{\infty}\right\\}$ is Q-linearly converging to zero
with rate $\frac{1+\gamma}{1-\gamma}\alpha$. If, in addition,
$\lim_{k\rightarrow\infty}\alpha_{k}=0$, then
$\left\\{\|V_{k}-V^{*}\|_{\infty}\right\\}$ is converging to zero
Q-superlinearly.

Proof: We refer to the proof of Theorem 4 in gargiani_2023  for the derivation
of the following inequality

$\|V_{k+1}-V^{*}\|_{\infty}\leq\frac{1}{1-\gamma}\left[\|r(V_{k})-r(V^{*})-J_{k}(V_{k}-V^{*})\|_{\infty}\\!+(1+\gamma)\alpha_{k}\|V_{k}-V^{*}\|_{\infty}\\!\right]\\!.$
(50)

Making use of Lemma 10, we can upper-bound the first term in the RHS of (50)
and obtain the following inequality

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}$
$\displaystyle\leq\frac{1}{1-\gamma}\left[C_{\min}\|V_{k}-V^{*}\|^{2}_{\infty}+(1+\gamma)\alpha_{k}\|V_{k}-V^{*}\|_{\infty}\right]$
(51)
$\displaystyle=\left[\frac{C_{\min}}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}+\frac{(1+\gamma)}{(1-\gamma)}\alpha_{k}\right]\|V_{k}-V^{*}\|_{\infty}\,.$

Finally, since $\alpha_{k}\leq\alpha$, we obtain the following upper-bound

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}$
$\displaystyle\leq\underbrace{\left[\frac{C_{\min}}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}+\frac{(1+\gamma)}{(1-\gamma)}\alpha\right]}_{:=\delta_{k}}\|V_{k}-V^{*}\|_{\infty}\,.$
(52)

Convergence of the sequence $\left\\{\|V_{k}-V^{*}\|_{\infty}\right\\}$ to
zero follows from the fact that
$\|V_{1}-V^{*}\|_{\infty}\,\leq\,\delta_{0}\,\|V_{0}-V^{*}\|_{\infty}$ and
$\delta_{0}<1$, which is an implication of the assumption that
$V_{0}\in\mathcal{B}(V^{*},\omega)$. Starting from (52) and using these
inequalities recursively lead to conclude that

$\|V_{k}-V^{*}\|_{\infty}\leq\,\delta_{0}^{k}\,\|V_{0}-V^{*}\|_{\infty}\,.$
(53)

This means that, under the considered assumptions, the sequence
$\left\\{\|V_{k}-V^{*}\|_{\infty}\right\\}$ converges to zero at least
linearly with rate $\delta_{0}$. To obtain the local contraction rate we take
the limit for $k\rightarrow\infty$ of (52) as follow

$\displaystyle\lim_{k\rightarrow\infty}\frac{\|V_{k+1}-V^{*}\|_{\infty}}{\|V_{k}-V^{*}\|_{\infty}}=\lim_{k\rightarrow\infty}\frac{C_{\min}}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}+\frac{(1+\gamma)}{(1-\gamma)}\alpha=\frac{(1+\gamma)}{(1-\gamma)}\alpha<1\,,$
(54)

where the second equality follows from (53) and the strict inequality follows
from the assumption that $\alpha\in(0,\frac{1-\gamma}{1+\gamma})$. If, in
addition, $\lim_{k\rightarrow\infty}\alpha_{k}=0$, then we can repeat similar
steps starting from (51) and conclude local Q-superlinear convergence. While
the convergence results of OPI (see for instance Proposition 9) fails in
capturing the improvement in convergence rate that results from running more
inner iterations, and, consequently, from solving with higher accuracy the
policy evaluation step, the results of Theorem 13 show that smaller values of
the $\alpha$-parameter, which translates into solving with a higher accuracy
the policy evaluation step, lead to a better local convergence rate of the
method.

While it is possible to show local contraction with relatively simple steps,
proving global convergence of a general iPI method is not equally simple.
Intuitively, from the convergence results of OPI in Proposition 9 (or also
Proposition in DB_book ), one may expect that, with no additional assumptions
on the structural properties of the underlying MDP, for any $\gamma\in(0,1)$
there exists a non-empty set of positive $\alpha$ for which we can get global
convergence guarantees for a general iPI method. As we discuss in the
following, despite its intuitiveness, proving that is quite challenging.
Differently from the global convergence proof of OPI (see Proposition 9), we
can not exploit the fundamental properties of the Bellman operators tout court
as, depending on the specific selection, the inner solver may not be a
contraction in the infinity-norm.

We start the discussion on global convergence of iPI methods by illustrating
possible technical complications. Then, by introducing additional assumptions
on the structure of the underlying problem, we prove global convergence of a
general iPI method.

######  Proposition 14.

Consider a general infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ and $\left\\{\pi_{k+1}\right\\}$ be the sequences
generated by a general iPI method as described in Algorithm 3, where
$V_{0}\in\mathbb{R}^{n}$. Then for any $k\geq 0$

$\|V_{k+1}-V^{*}\|_{\infty}\leq\frac{(1+\gamma)}{(1-\gamma)}\alpha\|V_{k}-V^{*}\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}\,.$
(55)

Proof: In the following, we use $J_{k}=I-\gamma P^{\pi_{k+1}}$ and $\Delta
V_{k}=V_{k+1}-V_{k}$. We start by algebraically manipulating $V_{k+1}-V^{*}$
and then we use the fundamental properties of the infinity-norm, the
structural properties of the Bellman residual function and the stopping
condition to obtain the final result as follow

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}$ $\displaystyle=\|V_{k}+\Delta
V_{k}-V^{*}\|_{\infty}$ $\displaystyle=\|V_{k}+\Delta
V_{k}-V^{\pi_{k+1}}+V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\|V_{k}+\Delta
V_{k}-V^{\pi_{k+1}}\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle=\|J_{k}^{-1}g^{\pi_{k+1}}\\!-J_{k}^{-1}g^{\pi_{k+1}}+J_{k}^{-1}J_{k}V_{k}+J_{k}^{-1}J_{k}\Delta
V_{k}-J_{k}^{-1}J_{k}V^{\pi_{k+1}}\|_{\infty}$
$\displaystyle\phantom{=\,}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle=\|J_{k}^{-1}\left(g^{\pi_{k+1}}\\!-J_{k}V^{\pi_{k+1}}\right)+J_{k}^{-1}\left(J_{k}V_{k}-g^{\pi_{k+1}}+J_{k}\Delta
V_{k}\right)\|_{\infty}$
$\displaystyle\phantom{=\,}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle=\|J_{k}^{-1}\left(J_{k}V_{k}-g^{\pi_{k+1}}+J_{k}\Delta
V_{k}\right)\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\|J_{k}^{-1}\|_{\infty}\|J_{k}V_{k}-g^{\pi_{k+1}}+J_{k}\Delta
V_{k}\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\overset{(a)}{\leq}\frac{1}{1-\gamma}\|J_{k}V_{k}-g^{\pi_{k+1}}+J_{k}\Delta
V_{k}\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\frac{1}{1-\gamma}\alpha\|r(V_{k})-r(V^{*})\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\overset{(b)}{\leq}\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\|V^{\pi_{k+1}}-V^{*}\|_{\infty}\,,$

where $(a)$ and $(b)$ are obtained by deploying Lemma 3 and Lemma 2 in
gargiani_2023 , respectively.

Starting from (55) it is very tempting to proceed by deploying the convergence
results from policy iteration and, in particular, Inequality (16) in
Proposition 8, as $V^{\pi_{k+1}}$ is obtained by applying a full policy
iteration step (policy improvement and policy evaluation) starting from
$V_{k}$. That would allow one to conclude a sufficient condition on $\alpha$,
similar to the one in Proposition 13, and that guarantees global convergence
of iPI methods. Unfortunately, without additional assumptions, for the
following upper-bound to hold

$\|V^{\pi_{k+1}}-V^{*}\|_{\infty}\leq\gamma\,\|V_{k}-V^{*}\|_{\infty}\,,$ (56)

$V_{k}$ would need to be the cost associated to a policy $\pi\in\Pi$. Because
of the inexactness in the policy evaluation step, this condition is not
necessarily verified, therefore we can not deploy (56).

Without introducing additional assumptions, a tight upper-bound on
$\|V^{\pi_{k+1}}-V^{*}\|_{\infty}$ is given by Proposition 2.3.3 in DB_book .
In particular, by deploying

$\|V^{\pi_{k+1}}-V^{*}\|_{\infty}\leq\frac{2\gamma}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}\,$
(57)

in (55), we obtain the following upper-bound

$\|V_{k+1}-V^{*}\|_{\infty}\leq\frac{(1+\gamma)\alpha+2\gamma}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}$
(58)

which allows one to conclude global convergence if $\gamma<1/3$ and
$\alpha\in(0,(1-3\gamma)/(1+\gamma))$. These results are unsatisfactory as
they are only valid for a limited range of discount factor values. By
exploiting also the contractivity of the Bellman operator, in the next
proposition we provide a similar sufficient condition to guarantee global
convergence of a general iPI method, but for a slightly wider range of
discount factor values.

######  Proposition 15.

Consider a general infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ be the sequence of iterates generated by a general
iPI method as described in Algorithm 3. If $\gamma\in(0,\,-1+\sqrt{2})$ and
$\alpha\in(0,\,(1-2\gamma-\gamma^{2})/(1+\gamma))$, then
$(\alpha+\gamma)(1+\gamma)+\gamma<1$ and for any $k\geq 0$

$\|V_{k+1}-V^{*}\|_{\infty}\leq((\alpha+\gamma)(1+\gamma)+\gamma)\|V_{k}-V^{*}\|_{\infty}\,.$
(59)

Proof: We derive an upper-bound on $\|V_{k+1}-V^{*}\|_{\infty}$ by deploying
the properties of the infinity-norm, the stopping condition and the structural
properties of the Bellman residual function and the Bellman operators. This
leads to the following upper-bound

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}=\|V_{k+1}-V_{k}+V_{k}-TV_{k}+TV_{k}-V^{*}\|_{\infty}$
(60) $\displaystyle=\|\Delta V_{k}+r(V_{k})+TV_{k}-V^{*}\|_{\infty}$
$\displaystyle=\|\Delta V_{k}+\left(I-\gamma
P^{\pi_{k+1}}\right)r(V_{k})+\gamma
P^{\pi_{k+1}}r(V_{k})+TV_{k}-V^{*}\|_{\infty}$ $\displaystyle\leq\|\Delta
V_{k}+\left(I-\gamma P^{\pi_{k+1}}\right)r(V_{k})\|_{\infty}+\|\gamma
P^{\pi_{k+1}}r(V_{k})\|_{\infty}+\|TV_{k}-V^{*}\|_{\infty}$
$\displaystyle\overset{(a)}{\leq}\alpha\|r(V_{k})\|_{\infty}+\|\gamma
P^{\pi_{k+1}}\left(r(V_{k})-r(V^{*})\right)\|_{\infty}+\|TV_{k}-TV^{*}\|_{\infty}$
$\displaystyle\overset{\phantom{(c)}}{\leq}\alpha\|r(V_{k})\|_{\infty}+\gamma\|P^{\pi_{k+1}}\|_{\infty}\|\left(r(V_{k})-r(V^{*})\right)\|_{\infty}+\|TV_{k}-TV^{*}\|_{\infty}$
$\displaystyle\overset{(b)}{\leq}\alpha(1+\gamma)\|V_{k}-V^{*}\|_{\infty}+\gamma(1+\gamma)\|V_{k}-V^{*}\|_{\infty}+\gamma\|V_{k}-V^{*}\|_{\infty}$
$\displaystyle=((\alpha+\gamma)(1+\gamma)+\gamma)\|V_{k}-V^{*}\|_{\infty}\,,$

where $(a)$ follows from the stopping condition and $(b)$ from Lemma 2 in
gargiani_2023  and Proposition 2. The final result follows from observing that
$(\alpha+\gamma)(1+\gamma)+\gamma\in(0,1)$ if and only if
$\gamma\in(0,-1+\sqrt{2})$ and
$\alpha\in(0,(1-2\gamma-\gamma^{2})/(1+\gamma))$.

Going back to Inequality (55), we now assume that the sequence of iterates
$\left\\{V_{k}\right\\}$ generated by a general iPI method is such that
$V_{k}\geq TV_{k}$ for $k\geq 0$. This requirement has a strong relation with
the concept of region of decreasing as introduced in (yuchao_phd, , Section
5.2). In particular, it corresponds to assuming that $\left\\{V_{k}\right\\}$
is contained in this subset. The assumption is also comparable to some sort of
policy improvement requirement as it implies that $V_{k}\geq V^{\pi_{k+1}}$.
Indeed, by definition of greedy-policy, we have that
$TV_{k}=T_{\pi_{k+1}}V_{k}$. Because of the monotonicity of the
$T_{\pi}$-operator for any $\pi\in\Pi$, we have that $V_{k}\geq
T_{\pi_{k+1}}V_{k}$ implies that $T_{\pi_{k+1}}V_{k}\geq V^{\pi_{k+1}}$. In
addition, by definition of optimal cost, $V^{\pi}\geq V^{*}$ for any
$\pi\in\Pi$. Putting these pieces together we have that

$TV_{k}\geq V^{\pi_{k+1}}\geq V^{*}\,.$ (61)

By subtracting the term $V^{*}$ using the fact that that $V^{*}=TV^{*}$, we
obtain that

$TV_{k}-TV^{*}\geq V^{\pi_{k+1}}-V^{*}\geq 0\,,$ (62)

which implies the following

$\|TV_{k}-TV^{*}\|_{\infty}\geq\|V^{\pi_{k+1}}-V^{*}\|_{\infty}\,.$ (63)

By using (63) in (55) and exploiting the fact that the $T$-operator is
$\gamma$-contractive in the infinity-norm, we obtain the following upper-bound

$\|V_{k+1}-V^{*}\|_{\infty}\leq\frac{(1+\gamma)\alpha}{1-\gamma}\|V_{k}-V^{*}\|_{\infty}+\gamma\|V_{k}-V^{*}\|_{\infty}\,.$
(64)

From (64) we can then conclude global Q-linear convergence of a general iPI
method with rate $\frac{(1+\gamma)\alpha}{1-\gamma}+\gamma$ if
$\alpha\in(0,\,(1-\gamma)^{2}/(1+\gamma))$. Having an assumption on the
generated sequence of iterates is not ideal as it is not possible to verify it
a priori. It would therefore be more practical to modify the algorithm in
order to enforce this condition, but this would obviously radically change the
general iPI scheme.

In the following, we discuss global convergence of iPI methods. In particular,
first we show global convergence to a neighborhood of the solution for the
same interval of $\alpha$-values for which local contraction is guaranteed.
Then, by introducing an extra requirement on the stage-cost bound, we show
that there exists a non-empty interval of $\alpha$-values for which global
convergence to the solution is guaranteed. Differently from Proposition 15,
these results hold for any $\gamma\in(0,1)$.

######  Theorem 16 (global convergence to a neighborhood of $V^{*}$).

Consider a general infinite-horizon discounted MDP with finite spaces and let
$\left\\{V_{k}\right\\}$ be the sequence of iterates generated by a general
iPI method as described in Algorithm 3. If
$\alpha\in(0,\,(1-\gamma)/(1+\gamma))$, then, for any
$V_{0}\in\mathbb{R}^{n}$, $\left\\{V_{k}\right\\}$ asymptotically converges to
$\mathcal{B}(V^{*},2R/((1-\gamma)(1-\zeta)))$, where
$\zeta=\frac{1+\gamma}{1-\gamma}\alpha$.

Proof: In the following, we denote with $\left\\{\pi_{k+1}\right\\}_{k\geq 0}$
and $\left\\{\pi^{\text{\tiny{PI}}}_{k}\right\\}_{k\geq 0}$ the sequences of
greedy policies that are generated by iPI and exact PI starting from $V_{0}$
and $\tilde{V}$, respectively. Let $\bar{V}\in\mathbb{R}^{n}$, we assume that
$V_{0}=\tilde{V}=\bar{V}$ and that the first policy for both sequences is
initialized with the same greedy policy with respect to $\bar{V}$. Therefore,
$\pi_{1}=\pi^{\text{PI}}_{0}$ while for $k>0$ the equality does not
necessarily hold since in general $V_{k}\neq V^{\pi_{k}}$. We also use
$J_{\pi_{k+1}}$ to denote $I-\gamma P^{\pi_{k+1}}$ and
$J_{\pi^{\text{PI}}_{k+1}}$ to denote $I-\gamma P^{\pi^{\text{PI}}_{k+1}}$.
Following similar steps as in the proofs of Proposition 14 and 15, we upper-
bound $\|V_{k+1}-V^{*}\|_{\infty}$ by using the fundamental properties of the
infinity-norm, the structural properties of the Bellman residual function and
the stopping-condition. In particular, we obtain the following upper-bound

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}=\|V_{k}-V^{*}+\Delta
V_{k}\|_{\infty}$ (65)
$\displaystyle=\|J_{\pi_{k+1}}^{-1}\left(J_{\pi_{k+1}}V_{k}-g^{\pi_{k+1}}+J_{\pi_{k+1}}\Delta
V_{k}\right)+J_{\pi_{k+1}}^{-1}g^{\pi_{k+1}}+V^{\pi^{\text{PI}}_{k+1}}-V^{*}-V^{\pi^{\text{PI}}_{k+1}}\|_{\infty}$
$\displaystyle\leq\|J_{\pi_{k+1}}^{-1}\left(J_{\pi_{k+1}}V_{k}-g^{\pi_{k+1}}+J_{\pi_{k+1}}\Delta
V_{k}\right)\|_{\infty}+\|J_{\pi_{k+1}}^{-1}g^{\pi_{k+1}}-V^{\pi^{\text{PI}}_{k+1}}\|_{\infty}$
$\displaystyle\phantom{\leq}\,\,+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle=\|J_{\pi_{k+1}}^{-1}\left(J_{\pi_{k+1}}V_{k}-g^{\pi_{k+1}}+J_{\pi_{k+1}}\Delta
V_{k}\right)\|_{\infty}+\|J_{\pi_{k+1}}^{-1}g^{\pi_{k+1}}-J_{\pi^{\text{PI}}_{k+1}}^{-1}g^{\pi^{\text{PI}}_{k+1}}\|_{\infty}$
$\displaystyle\phantom{\leq}\,\,+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\max_{\pi\,\,\in\,\,\left\\{\pi_{k+1},\,\pi^{\text{PI}}_{k+1}\right\\}}\|J_{\pi}^{-1}(g^{\pi_{k+1}}-g^{\pi^{\text{PI}}_{k+1}})\|_{\infty}$
$\displaystyle\phantom{\leq}\,\,+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\max_{\pi\,\,\in\,\,\left\\{\pi_{k+1},\,\pi^{\text{PI}}_{k+1}\right\\}}\|J_{\pi}^{-1}\|_{\infty}\|g^{\pi_{k+1}}-g^{\pi^{\text{PI}}_{k+1}}\|_{\infty}$
$\displaystyle\phantom{\leq}\,\,+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\leq\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\frac{1}{1-\gamma}\|g^{\pi_{k+1}}-g^{\pi^{\text{PI}}_{k+1}}\|_{\infty}+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}$
$\displaystyle\overset{(a)}{\leq}\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\|V^{\pi^{\text{PI}}_{k+1}}-V^{*}\|_{\infty}+\frac{2R}{1-\gamma}$
$\displaystyle\overset{(b)}{\leq}\frac{1+\gamma}{1-\gamma}\alpha\|V_{k}-V^{*}\|_{\infty}+\gamma^{k+1}\|V^{\pi_{0}}-V^{*}\|_{\infty}+\frac{2R}{1-\gamma}\,,$

where Inequality $(a)$ follows from the fact that, for any $\pi\in\Pi$,
$g^{\pi}:\mathcal{S}\times\mathcal{A}\rightarrow[-R,\,R]$, and Inequality
$(b)$ follows from Proposition 8.

By applying Inequality (65) recursively, we obtain the following upper-bound

$\displaystyle\|V_{k+1}-V^{*}\|_{\infty}\leq$
$\displaystyle\zeta^{k+1}\|V_{0}-V^{*}\|_{\infty}+\|V^{\pi_{0}}-V^{*}\|_{\infty}\gamma^{k+1}\sum_{j=0}^{k}\left(\frac{\zeta}{\gamma}\right)^{j}$
(66) $\displaystyle\quad\,+\frac{2R}{1-\gamma}\sum_{j=0}^{k}\zeta^{j}\,,$

where $\zeta=\frac{1+\gamma}{1-\gamma}\alpha$. We now take the limit for
$k\rightarrow\infty$ of (66). Since $\alpha\in(0,(1-\gamma)/(1+\gamma))$, then
$\zeta<1$ and therefore
$\lim_{k\rightarrow\infty}\zeta^{k+1}\|{V}_{0}-V^{*}\|_{\infty}=0$ and
$\lim_{k\rightarrow\infty}\frac{2R}{1-\gamma}\sum_{j=0}^{k}\zeta^{j}=\frac{2R}{(1-\gamma)(1-\zeta)}$.
To study the limit of the second term in the RHS of (66), we distinguish
between three scenarios: $\zeta<\gamma$, $\zeta>\gamma$ and $\zeta=\gamma$. In
the first scenario, by using the properties of the geometric series with
common ratio strictly less than $1$, we can directly conclude that
$\lim_{k\rightarrow\infty}\gamma^{k+1}\sum_{j=0}^{k}\left(\frac{\zeta}{\gamma}\right)^{j}=0\cdot\frac{1}{1-\frac{\zeta}{\gamma}}=0$.
For the second scenario, we can apply a similar argument starting from the
observation that
$\gamma^{k+1}\sum_{j=0}^{k}\left(\frac{\zeta}{\gamma}\right)^{j}$ can be
equivalently rewritten as
$\gamma\zeta^{k}\sum_{j=0}^{k}\left(\frac{\gamma}{\zeta}\right)^{j}$. Finally,
when $\zeta=\gamma$,
$\gamma^{k+1}\sum_{j=0}^{k}\left(\frac{\zeta}{\gamma}\right)^{j}=(k+1)\gamma^{k+1}$.
By using l’H$\hat{\text{o}}$pital’s rule, we can conclude that
$\lim_{k\rightarrow\infty}(k+1)\,\gamma^{k+1}=\lim_{k\rightarrow\infty}\frac{\gamma^{k-1}}{k}=0$.
Therefore, by combining these derivations, we obtain the following upper-bound

$\lim_{k\rightarrow\infty}\|V_{k+1}-V^{*}\|_{\infty}\leq\frac{2R}{(1-\gamma)(1-\zeta)}\,,$
(67)

from which we can conclude that the sequence of iterates generated by a
general iPI method asymptotically converges to
$\mathcal{B}(V^{*},\,2R/((1-\gamma)(1-\zeta)))$.

######  Corollary 17.

Consider the same conditions of Theorem 16. If
$R<\frac{(1-\gamma)^{3/2}\delta}{2}$, then
$\left\\{\|V_{k}-V^{*}\|_{\infty}\right\\}$ asymptotically converges to zero
for any $V_{0}\in\mathbb{R}^{n}$.

Proof:

We deploy the results of Theorem 13 and 16 and derive a condition on $\alpha$
such that the estimated radius of the global region of convergence is smaller
than the radius of the local region of attraction. This boils down to solving
the following inequality in $\alpha$

$\frac{2R}{(1-\gamma)(1-\zeta)}\leq\frac{(1-\gamma)-(1+\gamma)\alpha}{C_{\min}}\,.$
(68)

With some algebraic manipulation, we can rewrite (68) as follow

$(1+\gamma)^{2}\alpha^{2}-2(1-\gamma)(1+\gamma)\alpha+(1-\gamma)^{2}-2RC_{\min}\geq
0\,.$ (69)

Inequality (69) is verified for $\alpha\leq\alpha_{1}$ and
$\alpha\geq\alpha_{2}$, where

$\alpha_{1,2}=\frac{(1-\gamma)\pm\sqrt{2RC_{\min}}}{(1+\gamma)}\,.$ (70)

The second interval can not be accepted since to deploy the local contraction
results we need $\alpha\in(0,(1-\gamma)/(1+\gamma))$ and
$\alpha_{1}\geq(1-\gamma)/(1+\gamma)$ by construction since $R,C_{\min}\geq
0$. The first interval is valid as it is always contained in
$(0,(1-\gamma)/(1+\gamma))$ if $R<(1-\gamma)^{3/2}\delta/2$, where the latter
requirement indeed ensures that $\sqrt{2RC_{\min}}<(1-\gamma)$.

Algorithm 3 is well-posed if the stopping condition for the inexact policy
evaluation step can be verified by the iterative linear solver in a finite
number of iterations. The following proposition provides guarantees that the
stopping condition is verified in a finite number of iterations whenever the
iterative linear solver generates a sequence of residuals (iterates) that is
linearly contractive to zero (the solution) in some norm.

######  Proposition 18.

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. Let $z,\,\alpha\in(0,1)$ and $C_{*}\geq 1$ be constants,
$\|\cdot\|_{*}$ an arbitrary norm and
$\Phi^{\pi}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ the function with
$\Phi^{\pi}(V)=g^{\pi}-\left(I-\gamma P^{\pi}\right)V$. Consider an iterative
solver that generates a sequence of iterates $\left\\{\theta_{i}\right\\}$
such that $\|\Phi^{\pi}(\theta_{i})\|_{*}\leq
z^{i}\|\Phi^{\pi}(\theta_{0})\|_{*}$ or $\|\theta_{i}-V^{\pi}\|_{*}\leq
z^{i}\|\theta_{0}-V^{\pi}\|_{*}$ for all $i\geq 0$ and
$\theta_{0}\in\mathbb{R}^{n}$. Then

$\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq\alpha\,\|\Phi^{\pi}(\theta_{0})\|_{\infty}$
(71)

for all
$i\geq\lceil\log_{z}\frac{\alpha}{C_{*}}\frac{1-\gamma}{1+\gamma}\rceil$.

Proof: Because of the equivalence of norms on finite-dimensional spaces, we
know that there exists $C_{*}\geq 1$ such that, if
$\|\Phi^{\pi}(\theta_{i})\|_{*}\leq z^{i}\,\|\Phi^{\pi}(\theta_{0})\|_{*}$,
then $\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq
z^{i}\,C_{*}\,\|\Phi^{\pi}(\theta_{0})\|_{\infty}$. Analogously, if
$\|\theta_{i}-V^{\pi}\|_{*}\leq z^{i}\,\|\theta_{0}-V^{\pi}\|_{*}$, then
$\|\theta_{i}-V^{\pi}\|_{\infty}\leq
z^{i}\,C_{*}\,\|\theta_{0}-V^{\pi}\|_{\infty}$. From these relations we can
easily see that, in the first scenario, the condition
$\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq\alpha\|\Phi^{\pi}(\theta_{0})\|_{\infty}$
holds if $z^{i}C_{*}\leq\alpha$, which leads to
$i\geq\lceil\log_{z}\left({\alpha}/{C_{*}}\right)\rceil$. To find a lower-
bound on the number of iterations such that
$\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq\alpha\|\Phi^{\pi}(\theta_{0})\|_{\infty}$
for the second scenario, we start by upper-bounding the infinity-norm of
$\Phi^{\pi}(\theta_{i})$

$\displaystyle\|\Phi^{\pi}(\theta_{i})\|_{\infty}$
$\displaystyle=\|g^{\pi}-\left(I-\gamma P^{\pi}\right)\theta_{i}\|_{\infty}$
(72) $\displaystyle=\|\left(I-\gamma P^{\pi}\right)\left(\left(I-\gamma
P^{\pi}\right)^{-1}g^{\pi}-\theta_{i}\right)\|_{\infty}$
$\displaystyle\overset{(a)}{\leq}\|I-\gamma P^{\pi}\|_{\infty}\|\left(I-\gamma
P^{\pi}\right)^{-1}g^{\pi}-\theta_{i}\|_{\infty}$ $\displaystyle=\|I-\gamma
P^{\pi}\|_{\infty}\|V^{\pi}-\theta_{i}\|_{\infty}$
$\displaystyle\overset{(b)}{\leq}z^{i}\,C_{*}\,\|I-\gamma
P^{\pi}\|_{\infty}\|\theta_{0}-V^{\pi}\|_{\infty}$
$\displaystyle=z^{i}\,C_{*}\,\|I-\gamma P^{\pi}\|_{\infty}\|\left(I-\gamma
P^{\pi}\right)^{-1}\left(\left(I-\gamma
P^{\pi}\right)\theta_{0}-g^{\pi}\right)\|_{\infty}$
$\displaystyle\overset{(c)}{\leq}z^{i}\,C_{*}\,\|I-\gamma
P^{\pi}\|_{\infty}\|\left(I-\gamma
P^{\pi}\right)^{-1}\|_{\infty}\|\left(I-\gamma
P^{\pi}\right)\theta_{0}-g^{\pi}\|_{\infty}$
$\displaystyle=z^{i}\,C_{*}\,\|I-\gamma P^{\pi}\|_{\infty}\|\left(I-\gamma
P^{\pi}\right)^{-1}\|_{\infty}\|\Phi^{\pi}(\theta_{0})\|_{\infty}$
$\displaystyle\overset{(d)}{\leq}z^{i}\,C_{*}\,\frac{1+\gamma}{1-\gamma}\,\|\Phi^{\pi}(\theta_{0})\|_{\infty}\,,$

where Inequalities $(a)$ and $(c)$ follow from the application of the triangle
inequality, Inequality $(b)$ from the convergence properties of the iterative
linear solver in the infinity-norm and Inequality $(d)$ from the upper-bounds
on the infinity-norm of the matrix $I-\gamma P^{\pi}$ and its inverse. From
(72) it follows that
$\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq\alpha\|\Phi^{\pi}(\theta_{0})\|_{\infty}$
if $z^{i}\,C_{*}\,\frac{1+\gamma}{1-\gamma}\leq\alpha$, which leads to
$i\geq\lceil\log_{z}\left(\alpha/C_{*}\cdot\,(1-\gamma)/(1+\gamma)\right)\rceil$.
Therefore we can conclude that, since $(1-\gamma)/(1+\gamma)<1$ and $z<1$, in
both scenarios
$i\geq\lceil\log_{z}\left(\alpha/C_{*}\cdot\,(1-\gamma)/(1+\gamma)\right)\rceil$
is a sufficient number of iterations to guarantee that
$\|\Phi^{\pi}(\theta_{i})\|_{\infty}\leq\alpha\|\Phi^{\pi}(\theta_{0})\|_{\infty}$.

In the context of iPI, as described in Algorithm 3, at the $k$-th iteration we
consider the linear system associated to policy $\pi_{k+1}$ and $\theta_{0}$
is initialized with the current outer iterate $V_{k}$ (step 5 in Algorithm 3).
The approximate solution $\theta_{i}$ generated by the inner solver after
$i>0$ iterations is then used to update the outer iterate (step 11 in
Algorithm 3). The results of Proposition 18 ensure that iPI methods are well-
posed provided an inner solver which is $z$-contractive in some arbitrary
norm, since the stopping condition of the inner loop is verified in at most
$\lceil\log_{z}\frac{\alpha}{C_{*}}\frac{1-\gamma}{1+\gamma}\rceil$ inner
iterations.

###  4.2 Iterative Linear Solvers for Policy Evaluation

Inexact policy iteration methods require at every iteration the approximate
solution of a linear system of the following form

$\left(I-\gamma
P^{\pi}\right)\Delta\theta=-r^{\pi}(\bar{V})\quad\quad\pi\in\Pi\,,$ (73)

where $\bar{V}\in\mathbb{R}^{n}$, $\Delta\theta=\theta-\bar{V}$ and the
approximation level is regulated by the stopping condition in (49). Notice
that, given a $\pi\in\Pi$, (73) is equivalent to $(I-\gamma
P^{\pi})\theta=g^{\pi}$.

In the choice of the solver for a linear system of equations, the structural
properties of the coefficient matrix play a fundamental role. In particular,
for any $\pi\in\Pi$, the coefficient matrix in (73) has eigenvalues contained
in the circle with center $(1,0)$ and radius $\gamma$ in the complex plane
(gargiani_2023, , Lemma 1). Consequently, the coefficient matrix is always
non-singular, while we can not rely on symmetry unless the underlying MDP
enjoys a specific structure. Its symmetric component is

$H^{\pi,\,\gamma}=I-\frac{\gamma}{2}\left(P^{\pi}+{P^{\pi}}^{\top}\right)=I-\gamma
P^{\pi}_{s}\,.$ (74)

An approximate solution of (73) can therefore be provided by any linear and
non-linear iterative method for linear systems with a non-singular coefficient
matrix saad ; iterative_solutions . For a general iPI method, the choice of
the inner solver should though be tailored to the specific structure of (73)
and, as underlined in Proposition 18, it is of crucial importance in order to
obtain fast convergence. In this work we consider Richardson’s method
(iterative_solutions, , Chapter 3.2.1), steepest descent (iterative_solutions,
, Chapter 9.2), the minimal residual method (iterative_solutions, , Chapter
9.4) and GMRES gmres . For each of these methods we analyze the convergence
properties and computational costs when applied to solve (73).

####  4.2.1 Richardson’s Method

The iteration of Richardson’s method for (73) is

$\Delta{\theta}_{i+1}=\left(I-\frac{1}{\nu}\left(I-\gamma
P^{\pi}\right)\right)\Delta\theta_{i}-\frac{1}{\nu}r^{\pi}(\bar{V})\,,$ (75)

with $\nu>0$ and $\Delta{\theta}_{i}={\theta}_{i}-\bar{V}$ for any $i\geq 0$
and $\bar{V}\in\mathbb{R}^{n}$. For compactness, we use
$M_{\text{Rich},\,\nu}=I-\frac{1}{\nu}\left(I-\gamma P^{\pi}\right)$ to denote
its iteration matrix.

######  Proposition 19.

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. VI for policy evaluation is the instance of Richardson’s method
obtained by setting $\nu=1$.

Proof: Consider an arbitrary ${\theta}_{i}\in\mathbb{R}^{n}$ and we denote
with ${\theta}_{i+1}^{\text{VI}}$ and
$\Delta{\theta}_{i+1}^{\text{Rich},\,\nu}$ the VI for policy evaluation and
Richardson’s updates, respectively. Recall that the iteration of VI for policy
evaluation starting from ${\theta}_{i}$ is

${\theta}_{i+1}^{\text{VI}}=T_{\pi}{\theta}_{i}\,.$ (76)

We now consider the Richardson’s iteration starting from
$\Delta{\theta}_{i}={\theta}_{i}-\bar{V}$ with $\nu=1$

$\displaystyle\Delta{\theta}_{i+1}^{\text{Rich},\,1}$ $\displaystyle=\gamma
P^{\pi}\Delta{\theta}_{i}-r^{\pi}(\bar{V})$ (77) $\displaystyle=g^{\pi}+\gamma
P^{\pi}{\theta}_{i}-\bar{V}$ $\displaystyle=T_{\pi}{\theta}_{i}-\bar{V}\,.$

By simplifying the $\bar{V}$-terms on both sides of (77) we obtain that
${\theta}_{i+1}^{\text{Rich},\,1}={\theta}_{i+1}^{\text{VI}}$. Richardson’s
method can be therefore interpreted as a generalization of the value iteration
method for policy evaluation. In particular, the $\nu$-parameter plays a
fundamental role in determining its convergence properties. We now study the
values of $\nu$ for which Richardson’s iterations enjoy monotone convergence
in the infinity-norm (see Section 2.2.5 in iterative_solutions  for a
definition of monotone convergence with respect to a norm).

######  Proposition 20 (monotone convergence in the infinity-norm).

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. Let $\left\\{{\theta}_{i}\right\\}$ be the sequence of iterates
generated by (75) starting from $V_{0}\in\mathbb{R}^{n}$. If
$\nu>\frac{1+\gamma}{2}$, then for any $i\geq 0$

$\|{\theta}_{i+1}-V^{\pi}\|_{\infty}\leq\,h(\nu,\,\gamma)\,\|{\theta}_{i}-V^{\pi}\|_{\infty}\,,$
(78)

where $h(\nu,\gamma)=\frac{|\nu-1|}{\nu}+\frac{\gamma}{\nu}<1$.

Proof: A sufficient condition for monotone convergence in the infinity-norm is
$\|M_{\text{Rich},\,\nu}\|_{\infty}<1\,$ (iterative_solutions, , Th. 2.19). We
start by computing an upper-bound on the infinity-norm of the iteration matrix

$\|M_{\text{Rich},\,\nu}\|_{\infty}\leq\underbrace{\frac{|\nu-1|}{\nu}+\frac{\gamma}{\nu}}_{:=h(\nu,\,\gamma)}\,,$
(79)

where the inequality follows from the application of the triangle-inequality.
We now study the values of $\nu$ for which $h(\nu,\,\gamma)<1$. It is
straightforward to verify that, since $\gamma\in(0,1)$, $h(\nu,\,\gamma)<1$ if
and only if $\nu>\frac{1+\gamma}{2}$. The results of Proposition 20 allow one
to have an upper-bound on the infinity-norm of the error for some fixed value
of $m$. From this analysis, it is though not possible to conclude that there
exists values of $\nu$ which result into a faster monotone convergence than
value iteration for policy evaluation ($\nu=1$). It is also straightforward to
see that

$1=\arg\min_{\nu>\frac{1+\gamma}{2}}h(\nu,\,\gamma)\,,$ (80)

and $h(1,\,\gamma)=\gamma$. A necessary and sufficient condition on
convergence of Richardson’s method for policy evaluation can be formulated by
studying the values of $\nu$ for which
$\rho\left(M_{\text{Rich},\,\nu}\right)<1$ (iterative_solutions, , Th. 2.16).
In the literature, the spectral radius of the iteration matrix of a linear
iteration method is called convergence rate since it determines the asymptotic
convergence rate of the method. Differently from the contraction number
(iterative_solutions, , Section 2.2.5), the convergence rate only allows for
asymptotic statements concerning the convergence of the error. In particular,
$\rho(M_{\text{Rich},\,\nu})$ is a tight measure of the asymptotic convergence
rate of Richardson’s method. In the following theorem, we characterize the
class of MDPs for which it is always possible to select $\nu$ in order to
achieve asymptotic acceleration with respected to value iteration for policy
evaluation, i.e. $\rho(M_{\text{Rich},\,\nu})<\gamma$, and the interval of
$\nu$-values which leads to asymptotic acceleration. These results indicate
that $\rho(M_{\text{Rich}\,,1})=\gamma$ which suggests that for $\nu=1$ the
upper-bound on the contraction number in (79) is tight.

######  Theorem 21 (accelerated convergence rate).

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. If the MDP is regular, there always exists a non-empty interval
$\left(\underline{\nu},\,1\right)$, with

$\underline{\nu}:=\max_{\lambda\in\Lambda(P^{\pi})}\frac{\left(1-\gamma\text{Re}\left\\{\lambda\right\\}\right)-\gamma\sqrt{(\gamma^{2}-1)\text{Im}^{2}\left\\{\lambda\right\\}+\left(1-\gamma\text{Re}\left\\{\lambda\right\\}\right)^{2}}}{1-\gamma^{2}}\,,$
(81)

such that $\rho\left(M_{\text{Rich},\,\nu}\right)<\gamma$ for all
$\nu\in(\underline{\nu},\,1)$.

Proof: We start by characterizing $\rho(M_{\text{Rich},\,\nu})$ as follow

$\displaystyle\rho(M_{\text{Rich},\,\nu})$
$\displaystyle=\rho\left(I-\frac{1}{\nu}\left(I-\gamma P^{\pi}\right)\right)$
(82)
$\displaystyle=\max_{\lambda\in\Lambda(P^{\pi})}\Big{|}1-\frac{1}{\nu}\left(1-\gamma\lambda\right)\Big{|}\,.$

Since $\rho(P^{\pi})=1$, the choice $\nu=1$ leads to
$\rho(M_{\text{Rich}\,,1})=1$. In order to obtain an accelerated rate, there
has to exist $\underline{\nu},\,\bar{\nu}\in\mathbb{R}$ with
$\underline{\nu}<\bar{\nu}$ such that $\rho(M_{\text{Rich},\,\nu})<\gamma$ for
all $\nu\in(\underline{\nu},\,\bar{\nu})$. We therefore study the values of
$\nu$ for which

$\max_{\lambda\in\Lambda(P^{\pi})}\Big{|}1-\frac{1}{\nu}\left(1-\gamma\lambda\right)\Big{|}<\gamma\,.$
(83)

(83) is equivalent to

$\Big{|}1-\frac{1}{\nu}\left(1-\gamma\lambda\right)\Big{|}<\gamma\quad\forall\,\lambda\in\Lambda(P^{\pi})\,,$
(84)

which can be equivalently reformulated as follow

$\left(\frac{\nu-1}{\nu}+\frac{\gamma}{\nu}\text{Re}\left\\{\lambda\right\\}\right)^{2}+\left(\frac{\gamma}{\nu}\text{Im}\left\\{\lambda\right\\}\right)^{2}<\gamma^{2}\quad\forall\,\lambda\in\Lambda(P^{\pi})\,.$
(85)

After some basic algebraic manipulations, (85) becomes

$\nu^{2}(1-\gamma^{2})+2\nu(-1+\gamma\text{Re}\left\\{\lambda\right\\})+1+\gamma^{2}|\lambda|^{2}-2\gamma\text{Re}\left\\{\lambda\right\\}<0\quad\forall\,\lambda\in\Lambda(P^{\pi})\,.$
(86)

To solve the quadratic inequalities in (86), we first compute their
discriminant as a function of $\lambda\in\Lambda(P^{\pi})$

$\Delta(\lambda)=4\gamma^{2}\left((\gamma^{2}-1)\text{Im}^{2}\left\\{\lambda\right\\}+\left(1-\gamma\text{Re}\left\\{\lambda\right\\}\right)^{2}\right)\,.$
(87)

From (87) with some simple algebraic manipulations and since $\gamma\in(0,1)$,
we can conclude that

$\Delta(\lambda)\geq
0\quad\forall\,\lambda\in\Lambda(P^{\pi})\iff(1-\gamma^{2})\text{Im}^{2}\left\\{\lambda\right\\}\leq{(1-\gamma\text{Re}\left\\{\lambda\right\\})^{2}}\quad\forall\,\lambda\in\Lambda(P^{\pi})\,.$
(88)

We know that, since $P^{\pi}$ is row-stochastic, its eigenvalues lie in the
unit circle in the complex plane, i.e.,

$\text{Im}^{2}\left\\{\lambda\right\\}+\text{Re}^{2}\left\\{\lambda\right\\}\leq
1\quad\forall\,\lambda\in\Lambda(P^{\pi})\,.$ (89)

By multiplying left and right hand-sides in (89) by $(1-\gamma^{2})$ and with
some simple algebraic manipulations we obtain that the following condition
holds for all $\lambda\in\Lambda(P^{\pi})$

$(1-\gamma^{2})\text{Im}^{2}\left\\{\lambda\right\\}\leq(1-\gamma^{2})\left(1-\text{Re}^{2}\left\\{\lambda\right\\}\right)\,.$
(90)

Since
$(1-\gamma^{2})\left(1-\text{Re}^{2}\left\\{\lambda\right\\}\right)\leq\left(1-\gamma\text{Re}\left\\{\lambda\right\\}\right)^{2}$,
we can conclude that (90) is more stringent that the right condition in (88) .
Consequently, the right condition in (88) is always verified and therefore
$\Delta(\lambda)\geq 0$ for all $\lambda\in P^{\pi}$. We now solve the
inequalities in (86) and for each $\lambda\in P^{\pi}$ we obtain an interval
$\left(\nu_{1}(\lambda),\,\nu_{2}(\lambda)\right)$ where

$\nu_{1,2}(\lambda)=\frac{(1-\gamma\text{Re}\left\\{\lambda\right\\})\pm\gamma\sqrt{(\gamma^{2}-1)\text{Im}^{2}\left\\{\lambda\right\\}+(1-\gamma\text{Re}\left\\{\lambda\right\\})^{2}}}{(1-\gamma^{2})}\,.$
(91)

The solution of (86) is therefore obtained by intersecting these intervals as
follow

$\nu\in\bigcap_{\lambda\in\Lambda(P^{\pi})}\left(\nu_{1}(\lambda),\,\nu_{2}(\lambda)\right)=\left(\underline{\nu},\,\bar{\nu}\right)\,,$
(92)

where $\underline{\nu}=\max_{\lambda\in\Lambda(P^{\pi})}\nu_{1}(\lambda)$ and
$\bar{\nu}=\min_{\lambda\in\Lambda(P^{\pi})}\nu_{2}(\lambda)$.

We now show that $\bar{\nu}=1$. We first show that $\bar{\nu}\geq 1$ with a
reduction ad absurdum. We assume that there exists
$\hat{\lambda}\in\Lambda(P^{\pi})$ such that $\nu_{2}(\hat{\lambda})<1$. This
would imply that the following inequality holds

$\gamma\sqrt{(\gamma^{2}-1)\text{Im}^{2}\left\\{\hat{\lambda}\right\\}+(1-\gamma\text{Re}\left\\{\hat{\lambda}\right\\})^{2}}<\gamma\left(-\gamma+\text{Re}\left\\{\hat{\lambda}\right\\}\right)\,.$
(93)

We now distinguish two scenarios:
$\text{Re}\left\\{\hat{\lambda}\right\\}\leq\gamma$ and
$\text{Re}\left\\{\hat{\lambda}\right\\}>\gamma$. In the first case, since the
left hand-side in (93) is always non-negative, we can directly conclude that
the original assumption is absurd and therefore $\bar{\nu}\geq 1$. For the
second scenario, we can square both sides in (93) to obtain the following
equivalent condition

$(\gamma^{2}-1)\text{Im}^{2}\left\\{\hat{\lambda}\right\\}<(1-\gamma^{2})\left(\text{Re}^{2}\left\\{\hat{\lambda}\right\\}-1\right)\,.$
(94)

With some algebraic manipulations (94) can be reformulated as follow

$(1-\gamma^{2})(|\hat{\lambda}|^{2}-1)>0\,,$ (95)

which is verified if and only if $|\hat{\lambda}|>1$. Since all eigenvalue of
$P^{\pi}$ have modulus not greater than 1, also for this second scenario we
can conclude that the original assumption is absurd and therefore
$\bar{\nu}\geq 1$. Finally, since row-stochastic matrices have always an
eigenvalue at $(1,0)$ and $\nu_{2}(1)=1$, we can conclude that $\bar{\nu}=1$.
With a similar argument, we can prove that $\underline{\nu}\leq 1$ and that
the equality is achieved if and only if $|\lambda|^{2}=1$ and
$\text{Re}\left\\{\lambda\right\\}\leq\gamma$ with
$\lambda\in\Lambda(P^{\pi})$. While for a general MDP these conditions can be
verified since $P^{\pi}$ could have multiple eigenvalues with modulus 1, for
the case of regular MDPs we can conclude that $\underline{\nu}<1$ since 1 is
the only eigenvalue on the unit circle and all the other eigenvalues lie
strictly inside. We can therefore conclude that for a general MDP it may not
be possible to select $\nu$ to obtain an accelerated rate with respect to VI
since there could a policy $\pi\in\Pi$ for which
$\underline{\nu}=\bar{\nu}=1$, while for regular MDPs
$(\underline{\nu},\,\bar{\nu})$ is never empty since
$\underline{\nu}<\bar{\nu}=1$ and $\rho(M_{\text{Rich},\,\nu})<\gamma$ for all
$\nu\in(\underline{\nu},\,\bar{\nu})$. The results of Theorem 21 show that
value iteration for policy evaluation is not always optimal in the sense of
contraction rate. For the class of regular MDPs for any $\pi\in\Pi$ the
optimal $\nu$ in the sense of contraction rate is strictly less than 1. In
particular, given the policy evaluation task associated to any $\pi\in\Pi$ of
a regular MDP, it is always possible to achieve a faster convergence rate than
$\gamma$ by setting $\nu$ in an appropriate interval which depends on the
spectrum of the transition probability matrix and the discount factor. The
benchmarks in Figure 2 corroborate these results. As displayed in the figure
for a regular MDP with discount factor $\gamma=0.9$, by setting $\nu=0.7$ we
obtain a considerable improvement in terms of convergence rate with respect to
value iteration for policy evaluation, while for $\nu>1$ we observe monotone
convergence but with a slower rate.

Richardson’s iteration can also be used in combination with a non-singular
preconditioning matrix $D\in\mathbb{R}^{n\times n}$. This leads to the
following update

$\Delta{\theta}_{i+1}=\Delta{\theta}_{i}-\frac{1}{\nu}D^{-1}\left(\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}-g^{\pi}\right)\,.$ (96)

Notice that, if $\nu=1$ and $D=\text{diag}\left(I-\gamma P^{\pi}\right)$, then
we recover the Jacobi iteration (iterative_solutions, , Section 3.2.2).
Finally, starting from (96) and with $D=\text{diag}\left(I-\gamma
P^{\pi}\right)$, if the states are not processed all at once, but in a
sequential manner making use of the interim results, then we recover the SOR
iteration (iterative_solutions, , Section 3.2.4), and, if we set $\alpha=1$,
the Gauss-Seidel iteration (iterative_solutions, , Section 3.2.3). The use of
interim results may lead to a faster convergence in terms of number of
iterations, as proved for instance in gargiani_minibatch  and also confirmed
by the empirical results displayed in Figure 2. From a computational
viewpoint, SOR and Gauss-Seidel iterations are not suitable for parallel
systems since they are inherently sequential, while Richardson and Jacobi
iterations are fully parallelizable. The batch-versions of these latter
algorithms constitute a solution to trade-off the two aspects and therefore
may result into more efficient methods when a parallel system is deployed
gargiani_minibatch .

Unfortunately, the evaluation of (81) requires spectral data of the transition
probability matrix whose computation is generally more expensive than the
solution of the original problem. Therefore, practically, in order to enjoy
the acceleration in convergence rate for regular MDPs, it is necessary to
develop heuristic schemes in order to adjust $\nu$ avoiding the computation of
$\Lambda(P^{\pi})$. In the following, we analyze alternative methods which
adopt exact line-search strategies in order to select the step-size.

Figure 2: We consider a regular discounted MDP with $\gamma=0.9$, $n=100$ and
$m=50$ and we solve the policy evaluation task for a random policy $\pi\in\Pi$
with Richardson’s method with different values of $\nu$, Jacobi, Gauss-Seidel
and SOR. The figure displays the distance in infinity-norm of the iterates
from the solution versus number of iterations.

####  4.2.2 Steepest Descent

Instead of solving directly (20), we can equivalently minimize a quadratic
function of the following form

$q^{\pi,\,\gamma}(\Delta\theta\mathchar
24635\relax\;\Xi,\,c)=\frac{1}{2}\langle\Xi\left(\left(I-\gamma
P^{\pi}\right)\Delta\theta+r^{\pi}(\bar{V})\right),\,\left(I-\gamma
P^{\pi}\right)\Delta\theta+r^{\pi}(\bar{V})\rangle+c\,,$ (97)

where $\Xi\succ 0$ and $c$ is an arbitrary constant. The minimization of (97)
yields $V^{\pi}$ since any quadratic function of this form has a unique
minimum which coincides with the solution of (20) (iterative_solutions, ,
Lemma 9.3). In order to minimize (97), we can design an iteration of the
following form

$\Delta{\theta}_{i+1}=\Delta{\theta}_{i}+\eta_{i}\,d_{i}\,,$ (98)

where $d_{i}\in\mathbb{R}^{n}$ is a selected search-direction and $\eta_{i}>0$
is the optimal step-size obtained with exact line-search. These are the main
principles behind the steepest descent method, where $\Xi=I$, $c=0$ and
$d_{i}=-\nabla\,q^{\pi,\,\gamma}(\Delta{\theta}_{i}\mathchar
24635\relax\;I,\,0)$. In particular, with the choice of this quadratic
function and descent-direction, iteration (98) becomes

$\Delta{\theta}_{i+1}=\Delta{\theta}_{i}-\eta_{i}\left(I-\gamma
P^{\pi}\right)^{\top}\left(\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}+r^{\pi}\left(\bar{V}\right)\right)\,,$ (99)

where $\eta_{i}=\frac{\langle-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i},\,-\left(I-\gamma P^{\pi}\right)\nabla
q^{\pi,\,\gamma}(\Delta{\theta}_{i}\mathchar
24635\relax\;I,\,0)\rangle}{\langle\left(I-\gamma P^{\pi}\right)\nabla
q^{\pi,\,\gamma}(\Delta{\theta}_{i}\mathchar
24635\relax\;I,\,0),\,\,\left(I-\gamma P^{\pi}\right)\nabla
q^{\pi,\,\gamma}(\Delta{\theta}_{i}\mathchar 24635\relax\;I,\,0)\rangle}$.

The following proposition characterizes the convergence properties of
iteration (99).

######  Proposition 22.

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. The sequence $\left\\{{\theta}_{i}\right\\}$ generated by
steepest descent (99) with $V_{0}\in\mathbb{R}^{n}$ enjoys global Q-linear
convergence to $\Delta V^{*}=V^{\pi}-\bar{V}$ with contraction number
$\phi(\gamma,\,P^{\pi})=\dfrac{\kappa(\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right))-1}{\kappa(\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right))+1}$ in the
${\left(I-\gamma P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right)}$-norm,
i.e., for any $V_{0}\in\mathbb{R}^{n}$ and $i\geq 0$

$\|{\theta}_{i+1}-V^{\pi}\|_{\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma
P^{\pi}\right)}\,\leq\,\phi(\gamma,\,P^{\pi})\|{\theta}_{i}-V^{\pi}\|_{\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right)}\,.$ (100)

Proof: Each iteration (99) minimizes $q^{\pi,\,\gamma}(\Delta\theta\mathchar
24635\relax\;I,\,0)$ over all the vectors of the form $\Delta\theta-\eta\nabla
q^{\pi,\,\gamma}(\Delta\theta\mathchar 24635\relax\;I,\,0)$. This is also
equivalent to minimizing over the same family of vectors the following
quadratic function

$f(x)=\frac{1}{2}\|\Delta\theta-\Delta{\theta}^{*}\|^{2}_{\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right)}\,,$ (101)

where $\Delta{\theta}^{*}=V^{\pi}-\bar{V}=-\left(I-\gamma
P^{\pi}\right)^{-1}r^{\pi}(\bar{V})$. Since $\left(I-\gamma P^{\pi}\right)$ is
non-singular for any $\pi\in\Pi$, then this implies that $\left(I-\gamma
P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right)\succ 0$. We can therefore
utilize the results of Theorem 5.2 in saad  to conclude the proof, where
$A=\left(I-\gamma P^{\pi}\right)^{\top}\left(I-\gamma P^{\pi}\right)$ and
$b=-\left(I-\gamma P^{\pi}\right)^{\top}r^{\pi}(\bar{V})$.

####  4.2.3 Minimal Residual Method

The minimal residual method minimizes the same quadratic function as steepest
descent, but $d_{i}=-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}$. This results in the following iteration

$\Delta{\theta}_{i+1}=\Delta{\theta}_{i}+\eta_{i}\left(-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}\right)\,,$ (102)

where $\eta_{i}=\frac{\langle\left(I-\gamma
P^{\pi}\right)\left(-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}\right),\,-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}\rangle}{\langle\left(I-\gamma
P^{\pi}\right)\left(-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}\right),\,\left(I-\gamma
P^{\pi}\right)\left(-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{i}\right)\rangle}$. Notice that if we select a
constant step-size $\eta>0$ in iteration (102) we retrieve Richardson’s
iteration where $\nu=\frac{1}{\eta}$ in (75). While in general the selection
of the $\nu$-parameter in Richardson’s method requires knowledge of the
spectrum of the coefficient matrix, the minimal residual method adjusts this
value on the fly based on local information using exact line search.
Convergence results from the literature rely on the assumption that
$H^{\pi,\,\gamma}\succ 0$. In particular, under this assumption, the 2-norm of
the residual vectors of iteration (102) converges to zero with contraction
number at least
$\kappa_{\text{MR}}=\sqrt{1-\frac{\lambda^{2}_{\min}(H^{\pi\,,\gamma})}{\sigma^{2}_{\max}(I-\gamma
P^{\pi})}}$ (saad, , Th. 5.3). As proved by the following lemma, in the
considered problem setting the assumption $H^{\pi\,,\gamma}\succ 0$ is not
necessarily verified for all values of $\gamma\in(0,1)$.

Figure 3: We consider an MDP with $m=500$, $n=40$ and different values of
discount factor. We plot the minimum eigenvalue of $H^{\pi,\,\gamma}$ as a
function of $\gamma$. In line with the results of Lemma 23, for high-values of
the discount factor $\lambda_{\min}(H^{\pi,\,\gamma})<0$. For all tested
values of $\gamma$ in the plot, the minimal residual method has achieved
convergence to the solution. The convergence of some configurations in terms
of 2-norm of the residuals is plotted in Figure 4.

######  Lemma 23.

Consider a general infinite-horizon discounted MDP with finite spaces. For any
$\pi\in\Pi$, $1/\lambda_{\max}\left(P_{s}^{\pi}\right)\in(0,1]$ and
$H^{\pi\,,\gamma}$ as defined in (74) is positive definite if and only if
$\gamma<1/\lambda_{\max}\left(P_{s}^{\pi}\right)$.

Proof: Recall that, since for any $\pi\in\Pi$ $P^{\pi}$ is row-stochastic, its
eigenvalues lie in the circle with radius 1 and center $(0,0)$ in the complex
plane. Consequently, $\text{Re}\left\\{\lambda\right\\}\leq 1$ for all
$\lambda\in\Lambda(P^{\pi})$ and $(1,0)\in\Lambda\left(P^{\pi}\right)$. Since
$H^{\pi\,,\gamma}$ is symmetric by construction, studying its positive
definiteness boils down to studying if the following condition hold

$\lambda_{\min}\left(I-\frac{\gamma}{2}\left(P^{\pi}+{P^{\pi}}^{\top}\right)\right)>0\,.$
(103)

Since

$\lambda_{\min}\left(I-\frac{\gamma}{2}\left(P^{\pi}+{P^{\pi}}^{\top}\right)\right)=1-\gamma\lambda_{\max}\left(P_{s}^{\pi}\right)\,,$
(104)

(103) can also be reformulated as follow

$1-\gamma\lambda_{\max}\left(P_{s}^{\pi}\right)>0\,.$ (105)

We can therefore conclude that $H^{\pi\,,\gamma}\succ 0$ if and only if
$\gamma<\frac{1}{\lambda_{\max}(P^{\pi}_{s})}$. Finally, since
$\text{Re}\left\\{\lambda\right\\}\leq\lambda_{\max}\left(P^{\pi}_{s}\right)$
for all $\lambda\in\Lambda\left(P^{\pi}\right)$ (saad, , Th. 1.20), we can
conclude that $1/\lambda_{\max}(P^{\pi}_{s})\in(0,1]$. The results of Lemma 23
suggest that, in general, $H^{\pi,\,\gamma}\nsucc 0$ for large values of the
discount factor. Empirical evidence though shows that iteration (102) enjoys
convergence also when $H^{\pi\,,\gamma}\nsucc 0$, indicating that the positive
definite assumption can be restrictive. In particular, we conduct a benchmark
using iteration (102) to compute the cost associated with a policy $\pi$ for
an MDP with $n=500$, $m=40$ and different values of discount factor.
Convergence to $V^{\pi}$ is always achieved, even for high values of $\gamma$
for which, as displayed in Figure 3,
$\lambda_{\min}\left(H^{\pi\,,\gamma}\right)<0$. Finally, even when
$H^{\pi\,,\gamma}\succ 0$, the estimate on the contraction rate provided by
the literature tends to be very loose and therefore not representative of the
performance. An example of this is depicted in Figure 4, where we use
$\hat{\kappa}_{\text{MR}}$ to denote the estimated contraction rate.

Figure 4: We consider the same MDP and policy $\pi$ as in Figure 3 and vary
the discount factor. We solve iteratively the resulting policy evaluation
tasks with the minimal residual method. For each benchmark we report in the
legend $\kappa_{\text{MR}}$ and the empirically estimated contraction rate
$\hat{\kappa}_{\text{MR}}$ and we plot the 2-norm of the residuals versus
iteration number.  Figure 5: We consider a discounted MDP with $\gamma=0.7$,
$n=100$ and $m=50$ and we solve the policy evaluation task for a random policy
$\pi\in\Pi$ with Richardson’s method with $\nu=1$, the minimal residual
method, steepest descent and GMRES. The figure displays the distance in
infinity-norm of the iterates from the solution versus number of iterations.
Figure 6: We consider the same discounted MDP as in Figure 5 but with
$\gamma=0.99$ and we solve the policy evaluation task for a random policy
$\pi\in\Pi$ with Richardson’s method with $\nu=1$, the minimal residual
method, steepest descent and GMRES. The figure displays the distance in
infinity-norm of the iterates from the solution versus number of iterations.
Figure 7: We consider a discounted MDP with $\gamma=0.99$, $n=100$ and $m=50$
and we solve the policy evaluation task for a random policy $\pi\in\Pi$ with
Richardson’s method with $\nu=1$, the minimal residual method, steepest
descent and GMRES. The MDP is structurally different from the one deployed for
the benchmarks in Figure 6. The figure displays the distance in infinity-norm
of the iterates from the solution versus number of iterations.

####  4.2.4 GMRES

GMRES gmres  is an iterative method for linear systems with non-singular
coefficient matrices and it belongs to the class of Krylov subspace methods
(saad, , Chapter 6). Starting from an initial guess
$\Delta{\theta}_{0}\in\mathbb{R}^{n}$ with residual
$\Phi^{\pi}_{0}=-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta{\theta}_{0}$, GMRES generates a sequence
$\left\\{\Delta{\theta}_{i}\right\\}$ of approximate solutions to (73) with

$\Delta{\theta}_{i}=\arg\min_{\Delta\theta}\left\\{\|-r^{\pi}(\bar{V})-\left(I-\gamma
P^{\pi}\right)\Delta\theta\|_{2}\,:\,\Delta\theta\in\Delta{\theta}_{0}+\mathcal{K}_{i}\right\\}\,,$
(106)

where $\mathcal{K}_{i}=\text{span}\left\\{\Phi^{\pi}_{0},\,\left(I-\gamma
P^{\pi}\right)\Phi^{\pi}_{0},\,\left(I-\gamma
P^{\pi}\right)^{2}\Phi^{\pi}_{0},\,\dots,\,\left(I-\gamma
P^{\pi}\right)^{i-1}\Phi^{\pi}_{0}\right\\}$ and it is also known as the
$i$-th Krylov subspace (saad, , Chapter 6). In other words, the $i$-th iterate
is the vector that minimizes the 2-norm of the residual of (73) in the affine
space $\Delta{\theta}_{0}+\mathcal{K}_{i}$. Since working directly with
$\left\\{\Phi^{\pi}_{0},\,\left(I-\gamma
P^{\pi}\right)\Phi^{\pi}_{0},\,\left(I-\gamma
P^{\pi}\right)^{2}\Phi^{\pi}_{0},\,\dots,\,\left(I-\gamma
P^{\pi}\right)^{i-1}\Phi^{\pi}_{0}\right\\}$ is not computationally efficient
and assuming that $\|\Phi^{\pi}_{0}\|>0$, the method constructs an orthonormal
basis $\left\\{q_{1},\dots,q_{i}\right\\}$ where
$q_{1}=\Phi^{\pi}_{0}/\|\Phi^{\pi}_{0}\|_{2}$ and
$\left\\{q_{z}\right\\}_{z=2}^{i}$ are generated with the Arnoldi’s method
Arnoldi1951ThePO . The output of this orthonormalization algorithm is the
matrix $Q_{i}\in\mathbb{R}^{n\times i}$, whose columns are the vectors of the
orthonormal basis, and an Hessenberg matrix $H_{i}\in\mathbb{R}^{(i+1)\times
i}$ that satisfies the following relation

$AQ_{i}=Q_{i+1}H_{i}\,.$ (107)

Since any $\Delta\theta\in\Delta{\theta}_{0}+\mathcal{K}_{i}$ can be rewritten
as $\Delta\theta=\Delta{\theta}_{0}+Q_{i}y$ for some $y\in\mathbb{R}^{i}$ and
by exploiting (107), we can rewrite the residual associated to $\Delta\theta$
as a function of $y$ as follows

$\Phi_{i}^{\pi}=\Phi^{\pi}_{0}-\left(I-\gamma
P^{\pi}\right)Q_{i}y=Q_{i+1}\left(\|\Phi_{0}^{\pi}\|_{2}e_{1}-H_{i}y\right)\,,$
(108)

where $e_{1}=[1,0,\dots,0]^{\top}\in\mathbb{R}^{i+1}$. The $i$-th iterate is
selected as the $\Delta\theta\in\Delta{\theta}_{0}+\mathcal{K}_{i}$ which
minimizes the 2-norm of (108), and therefore as
$\Delta{\theta}_{i}=\Delta{\theta}_{0}+Q_{i}\tilde{y}$ where

$\tilde{y}=\arg\min_{y}\|\|\Phi_{0}^{\pi}\|_{2}e_{1}-H_{i}y\|_{2}\,.$ (109)

Unlike for the conjugate gradient method, an orthonormal basis of
$\mathcal{K}_{i}$ can not be computed with a short recurrence. Consequently,
when $i$ increases, the number of stored vectors also increases like $i$ and
the number of multiplications like $0.5i^{2}n$. A practical variant of GMRES
consists in restarting the algorithm after every $i$ iterations.

With exact arithmetic the sequence of iterates generated by GMRES converges to
the solution of (73) in at most $\ell=\sum_{\lambda\in\Lambda(I-\gamma
P^{\pi})}b(\lambda)\leq n$ iterations Campbell1996GMRESAT . For the large-
scale cases where the minimal polynomial has a high degree then GMRES may
require a significant number of iterations to reach convergence. In these
scenarios it becomes particularly important to characterize its convergence
properties in terms of improvement per iteration, as it is often intractable
to run $\ell$-iterations. In this regard, it is known in the literature that
GMRES has particularly favorable convergence properties when the eigenvalues
of the coefficient matrix are clustered in a circle of center $(1,0)$ and
radius $\xi<1$. In Campbell1996GMRESAT  the authors show that in these
scenarios the contraction rate is determined by the radius $\xi$ of the
cluster. These theoretical results are confirmed by extensive numerical
examples, where it is possible to observe that, while for coefficient matrices
with non-clustered eigenvalues the norm of the residuals stagnates and does
not notably decrease up to the very last iteration, for coefficient matrices
with clustered eigenvalues significant and steady progress is observed
starting from the first iteration (gargiani_2023, , Figure 1).

Linear systems arising from policy evaluation tasks in the considered problem
setting have always a coefficient matrix with eigenvalues in the circle with
center $(1,0)$ and radius $\gamma$ (gargiani_2023, , Lemma 1). Therefore,
GMRES operates in the favorable convergence regime. The following proposition
characterizes further the convergence properties of GMRES when deployed for
policy evaluation for general, ergodic and regular MDPs.

######  Proposition 24.

Consider a general infinite-horizon discounted MDP with finite spaces and the
linear system for the evaluation of the cost associated to a policy
$\pi\in\Pi$. The sequence $\left\\{\Phi^{\pi}({\theta}_{i})\right\\}$
generated by GMRES satisfies the following inequality

$\|\Phi^{\pi}({\theta}_{i})\|_{2}\leq
C_{1}\gamma^{i}\|\Phi^{\pi}({\theta}_{0})\|_{2}\quad i=1,\dots\,,$ (110)

where $C_{1}>0$ is a constant and it is independent of $i$. If the MDP is
ergodic and $h_{\pi}<n$, then

$\|\Phi^{\pi}({\theta}_{h_{\pi}+i})\|_{2}\leq
C_{2}|\lambda_{h_{\pi}+1}|^{i}\|\Phi^{\pi}({\theta}_{0})\|_{2}\quad
i=h_{\pi}+1,\dots\,,$ (111)

where $C_{2}>0$ is a constant and it is independent of $i$.

Proof: In the considered setting, the coefficient matrix of the linear system
associated to a policy evaluation task has always eigenvalues contained in the
circle of radius $\gamma$ and center $(1,0)$ (gargiani_2023, , Lemma 1). We
can therefore directly use the results of Proposition 2.2 in
Campbell1996GMRESAT  to obtain Inequality (110). The same results can be used
to prove Inequality (111) for ergodic MDPs, where the first $h_{\pi}$
eigenvalues are treated as outliers. From Proposition 5 and the fact that
$1\leq b(\lambda)\leq a(\lambda)$, we can conclude that $h_{\pi}$ iterations
are needed to process the first $h_{\pi}$ eigenvalues. After the first
$h_{\pi}$ iterations, we can use Inequality (4) in (Campbell1996GMRESAT, ,
Proposition 2.2), where the radius of the cluster is
$|\lambda_{h_{\pi}+1}|<\gamma$. Proposition 24 suggests that for general MDPs
the asymptotic convergence factor of GMRES is given by the radius $\gamma$ of
the circle where the eigenvalues are clustered. For ergodic MDPs, GMRES takes
$h_{\pi}$ iterations to process the eigenvalues of $P^{\pi}$ with modulus 1.
For this phase the improvement per iteration is still described by (110). Then
the asymptotic convergence factor is $|\lambda_{h_{\pi}+1}|<\gamma$. Because
of Proposition 6, the results of Proposition 24 further imply that for the
subclass of regular MDPs only one iteration is required to process the
eigenvalue $(1,0)$ of $P^{\pi}$ and then for $i>1$

$\|\Phi^{\pi}({\theta}_{i})\|_{2}\leq
C_{2}|\lambda_{2}|^{i}\|\Phi^{\pi}({\theta}_{0})\|_{2}\,.$ (112)

As shown in Proposition 24, for general and ergodic MDPs Richardson’s method
can not always achieve an asymptotic convergence rate faster than $\gamma$.
The latter is attained with the choice $\nu=1$, which leads to value iteration
for policy evaluation. These considerations and the results discussed in
Proposition 24 advocate in favor of GMRES for ergodic MDPs in terms of better
asymptotic contraction rate, especially when $h_{\pi}<<n$ and
$|\lambda_{h_{\pi}+1}|<<\gamma$.

##  5 Numerical Evaluation

This section is dedicated to numerically evaluating the performance of iPI
methods. For that we consider an MDP arising from a compartmental model
commonly deployed in epidemiology to mathematically model infectious diseases
hethcote . Details on the model are given in 5.1, while in 5.2 we comment on
the convergence and scalability performances of iPI methods when used to
compute health policies for a population with large size.

###  5.1 Infectious Disease Model

We consider a Susceptible-Infectious-Susceptible (SIS) model to describe the
evolution of an infectious disease with no immunity conferred by the previous
infection hethcote . This model is characterized by two classes: the
susceptible and the infective class. Individuals move from the susceptible
class to the infective class and then back to the susceptible class upon
recovery. As in sis_model , we assume that the infectious period is fixed and
equal to 1; that is, an individual that is infective at time $t$ remains
infectious over the interval $[t,t+1]$, but she/he will recover upon diagnosis
and effective treatment at time $t+1$, re-entering therefore in the
susceptible class. We denote with $s(t)$ and $i(t)$ the number of individuals
in the susceptible and infective class at time $t$, respectively. Since we
consider populations with a fixed-size, then

$s(t)+i(t)=N\quad\forall\,t\geq 0\,.$ (113)

From (113) we obtain that only the information of a single class are needed to
describe the state of the disease and we choose $s(t)$ to be the state of the
model. From now on, we drop the dependency on time, unless needed. We extend
the static model in sis_model  by adding a dynamic component which allows one
to derive health policies. In particular, we design actions, introduce an
action-state dependent stage-cost and make also the distribution of the
driving event in action-dependent. For the actions, we propose
$\mathcal{A}=\mathcal{A}_{1}\times\mathcal{A}_{2}$ as action set, where
$\mathcal{A}_{1}=\left\\{0,1,2,3,4\right\\}$ are the levels of hygiene
measures and $\mathcal{A}_{2}=\left\\{0,1,2,3\right\\}$ are the levels of
social distancing imposed by the public health authorities. Considering
influenza as example, levels of actions in increasing order in the first set
could correspond to no measures, frequent hand washing and disinfection,
mandatory surgical masks, mandatory FFP2 masks and mandatory full body
protection, respectively. Levels of actions in increasing order in the second
set could correspond to no social restrictions, mandatory social distancing,
restaurants and stores closure and full lockdown, respectively. As stage-cost,
we consider the following multi-objective cost function

$g(s,a)=w_{f}c_{f}(a)-w_{q}c_{q}(a)+w_{h}c_{h}(s)\,,$ (114)

where $w_{f},\,w_{q},\,w_{h}\geq 0$, $c_{f}:\mathcal{A}\rightarrow\mathbb{R}$
captures the financial costs and losses determined by the hygiene and social
measures that are put in place, $c_{q}:\mathcal{A}\rightarrow[0,1]$ assigns to
each action a quality of life score, and
$c_{h}:\mathcal{S}\rightarrow\mathbb{R}$ maps the number of infected people to
the medical cost incurred for their treatments.

The driving event from the susceptible class to the infective class is the
random variable $I(t)$, which represents the number of new infections
occurring during the interval $[t,\,t+1]$. As in sis_model , we model the
probability mass function for the driving event as a binomial distribution.
Hence, the probability mass function for the driving event for any
$a\in\mathcal{A}$ and $s\in\mathcal{S}$ is

$P\left(I(t)=i\,|\,s\right)=\begin{cases}\binom{s}{i}q(s,a)^{i}(1-q(s,a))^{s-i},\quad
0\leq i\leq s\\\ \,\,0,\quad\text{else},\end{cases}$ (115)

where $q:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$ with
$q(s,a)=1-\exp\left({-\lambda(a)\beta(s)\psi(a)}\right)$ is the overall
probability that a susceptible person becomes infected. The latter is a
function of $\psi(a)$, $\beta(s)$ and $\lambda(a)$, which are, respectively,
the probability that a susceptible person becomes infected upon contact with
an infectious individual, the probability that the next interaction of a
random susceptible person is with an infectious person, and the contact rate.
These parameters are specific of the considered infectious disease.

Finally, we consider $s(t)=N$ as an absorbing state since no infected
individuals remain in the population.

###  5.2 Benchmarks

We study the empirical convergence properties of iPI methods with different
inner solvers when they are used to solve infinite-horizon discounted MDPs
arising from dynamic SIS-models. For the benchmarks, we consider a dynamic SIS
model with $N=10000$. Fig. 10 displays the sparsity structure of the resulting
transition probability matrices.

We implemented PI and iPI methods in Python, relying on numpy and scipy
libraries for the dense and sparse linear algebra operations, respectively. We
refer to https://gitlab.ethz.ch/gmatilde/SIS_iPI for more details on the
selection of the other parameters of the model as well as for the code. All
benchmarks were run on a Core(TM) i7-10750H CPU @ 2.60GHz.

After the selection of the inner solver, we call the resulting iPI method
inexact inner solver name-policy iteration (iinnersolvername-PI). For
instance, if Richardson’s method is used as inner solver for the inexact
policy evaluation step, the resulting scheme is called iRich-PI. To guarantee
a finite-time execution of the code, in the case the stopping-condition is not
met after $500$ inner iterations, we terminate the inner loop and update the
outer solution. In addition, the overall algorithm is terminated if the
suboptimality falls below a given threshold or in case the execution time
exceeds 500 seconds.

As underlined by our theoretical analysis, the disocunt factor plays a
fundamental role in determining the convergence properties of the inner solver
and, therefore, the overall performance of an iPI scheme. Consequently, we
study two different scenarios: small ($\gamma=0.1$) and large ($\gamma=0.9$)
values of discount factor. In both cases and in line with the theoretical
analysis, PI is the method converging in the smallest number of iterations.
The inexactness of iPI methods results in a higher number of iterations
required to achieve convergence to the solution, where convergence is intended
within a certain tolerance. On the other hand, the CPU-time can be
dramatically reduced thanks to the inexactness and if an efficient inner
solver is deployed. iGMRES-PI is the fastest converging method in terms of
CPU-time and it is $\times$ 1.46 and $\times$ 1.56 faster than PI in the small
and large discount factor scenarios, respectively. iMinRes-PI shows also a
good and stable performance across the two scenarios, while iSD-PI is
competitive only in the small discount factor scenario. In the large discount
factor scenario it indeed fails to achieve convergence within the given time-
budget, and this is attributable to the exacerbating slow convergence rate of
SD as inner solver for this specific scenario. iRich-PI performance is also
not competitive since, in both scenarios, the method takes longer than PI to
achieve convergence in terms of CPU-time. Similarly as for the iSD-PI method,
this is ascribable to the slow convergence rate of the inner solver. In Figure
9(a) we study the impact of the discount factor on convergence in terms CPU-
time of iPI methods with $\alpha=0.1$ and PI when used to solve a dynamic SIS
MDP with $N=10000$. As depicted in the plot, for low-values of discount factor
iMinRes-PI is the fastest converging method, followed by iGMRES-PI, while for
high-values of discount factor iGMRES-PI attains the best performance. As in
the previous benchmarks, iRich-PI is not competitive and for high-values of
discount factor it even fails to achieve convergence within the given time and
inner-iterations limits. A similar consideration holds for iSD-PI, but for a
larger interval of discount factor values. In Figure 9(b) we benchmark iGMRES-
PI wth $\alpha=0.1$ and PI on dynamic SIS MDPs with different population sizes
and $\gamma=0.9$. As to be expected, for small-scale regimes PI tends to be
the fastest converging method, but its performance does not scale well. For
large-scale settings, indeed, iGMRES-PI converges up to more than $\times 3.7$
faster than PI.

From the conducted benchmarks we can conclude that, in line with the
theoretical results, inexactness of iPI methods results in a slower
convergence rate in terms of iterations, but a faster method in terms of CPU-
time, provided an appropriate choice of the inner solver. The selection of the
latter is indeed fundamental to obtain a competitive method.

(a) $\gamma=0.1$

(b) $\gamma=0.9$

Figure 8: We consider a SIS model with actions and population size of 10000.
We solve the associated infinite-horizon discounted MDP with PI and different
iPI methods with $\alpha=0.1$. In particular, as inner solvers, we consider
Richardson with $\nu=1$, steepest descent, minimal residual method and GMRES.
We study the empirical convergence properties of these methods for two
different values of discount factor and plot the distance in infinity-norm of
the iterates from the solution versus number of iterations. In the legend, we
also report the overall CPU-time in seconds.

(a) CPU-time versus discount factor values for a dynamic SIS MDP with
$N=10000$. For the iPI methods we set $\alpha=0.1$.

(b) CPU-time versus population sizes of dynamic SIS MDPs with $\gamma=0.9$.
For the iGMRES-PI we set $\alpha=0.1$.

Figure 9:  Figure 10: Graphical representation of the sparsity structure of
the transition probability matrices for the dynamic SIS model deployed to
generate the benchmarks in Fig. 8. White color is used for the zero entries
and black color for the non-zero entries. The $(i,j)$-plot displays the
transition probability matrix associated to the $(i,j)$-action pair.

##  6 Conclusions

We analyse a novel class of DP methods, inexact policy iteration methods,
which are based on an approximate solution of the policy evaluation. After
discussing their convergence properties and the importance of a good inner
solver selection, we consider different iterative methods for policy
evaluation and characterize their convergence behavior for the specific task
considered. We conclude by showcasing the competitive performance of inexact
policy iteration methods on large-scale MDPs arising from dynamic SIS models.
In line with the derived theoretical results, our numerical evidence shows
that value iteration as inner solver may result into a slow converging overall
scheme. A much better choice is represented by inner solvers that adjust step
and search direction based on locally available information, such as the
minimal residual method and GMRES.

Future extensions of the current work include an improved global convergence
analysis, which should encompass the need for extra assumptions on the
structural properties of the underlying MDP, as well as a distributed
implementation of inexact policy iteration schemes, which would allow to
overcome the memory limitations of a single machine and, consequently, to
scale to potentially gigantic problem sizes making use of modern high-
performance computer clusters.

## Acknowledgement

This work was supported by the European Research Council under the Horizon
2020 Advanced under Grant 787845 (OCAL).

## References

  * [1] Global inexact Newton methods for very large scale nonlinear problems.  IMPACT of Computing in Science and Engineering, 3(4):366–393, 1991.
  * [2] Newton-MR: Inexact Newton method with minimum residual sub-problem solver.  EURO Journal on Computational Optimization, 10:100035, 2022.
  * [3] W. E. Arnoldi.  The principle of minimized iterations in the solution of the matrix eigenvalue problem.  Quarterly of Applied Mathematics, 9:17–29, 1951.
  * [4] R. Bellman.  Dynamic Programming.  Princeton University Press, Princeton, NJ, USA, 1 edition, 1957.
  * [5] A. Berman and R. J. Plemmons.  Nonnegative Matrices in the Mathematical Sciences.  Society for Industrial and Applied Mathematics, 1994.
  * [6] D. P. Bertsekas.  Dynamic Programming and Optimal Control, Vol. II.  Athena Scientific, Belmont, Massachusets, 4th edition, 2012.
  * [7] N. Bäuerle and U. Rieder.  Markov Decision Processes with Applications to Finance.  Springer Berlin Heidelberg, 2011.
  * [8] S. L. Campbell, I. C. F. Ipsen, C. T. Kelley, and C. D. Meyer.  GMRES and the minimal polynomial.  BIT Numerical Mathematics, 36:664–675, 1996.
  * [9] R. S. Dembo, S. C. Eisenstat, and T. Steihaug.  Inexact newton methods.  SIAM Journal on Numerical Analysis, 19(2):400–408, 1982.
  * [10] F. Facchinei and J.-S. Pang.  Finite-Dimensional Variational Inequalities and Complementarity Problems, Vol. I.  Springer New York, NY, 1st edition, 2003.
  * [11] F. Facchinei and J.-S. Pang.  Finite-Dimensional Variational Inequalities and Complementarity Problems, Vol. II.  Springer New York, NY, 1st edition, 2003.
  * [12] A. Flueck and H.-D. Chiang.  Solving the nonlinear power flow equations with an inexact Newton method using GMRES.  IEEE Transactions on Power Systems, 13(2):267–273, 1998.
  * [13] M. Gargiani, D. Liao-McPherson, A. Zanelli, and J. Lygeros.  Inexact GMRES policy iteration for large-scale Markov decision processes.  IFAC-PapersOnLine, 56(2):11249–11254, 2023.  22nd IFAC World Congress.
  * [14] M. Gargiani, A. Martinelli, M. R. Martinez, and J. Lygeros.  Parallel and flexible dynamic programming via the mini-batch Bellman operator.  IEEE Transactions on Automatic Control, 69(1):455–462, 2024.
  * [15] M. Gargiani, A. Zanelli, D. Liao-McPherson, T. Summers, and J. Lygeros.  Dynamic programming through the lens of semismooth Newton-type methods.  IEEE Control Systems Letters, 6:1–1, 01 2022.
  * [16] M. Gargiani, A. Zanelli, D. Liao-McPherson, T. Summers, and J. Lygeros.  Dynamic programming through the lens of semismooth Newton-type methods (extended version).  https://arxiv.org/abs/2203.08678, 2022.
  * [17] S. Gerschgorin.  Uber die Abgrenzung der Eigenwerte einer Matrix.  Izvestija Akademii Nauk SSSR, Serija Matematika, 7(3):749–754, 1931.
  * [18] W. Hackbusch.  Iterative Solution of Large Sparse Systems of Equations.  Springer International Publishing, 2nd edition, 2016.
  * [19] H. W. Hethcote.  The mathematics of infectious diseases.  SIAM Review, 42(4):599–653, 2000.
  * [20] A. F. Izmailov and M. V. Solodov.  Newton-Type Methods for Optimization and Variational Problems.  Springer Cham, 1st edition, 2014.
  * [21] Y. Li.  Approximate methods of optimal control via dynamic programming models.  Doctoral Thesis in Electrical Engineering, KTH Royal Institute of Technology, 2023.
  * [22] J. Martínez and L. Qi.  Inexact Newton methods for solving nonsmooth equations.  Journal of Computational and Applied Mathematics, 60(1):127–145, 1995.  Proceedings of the International Meeting on Linear/Nonlinear Iterative Methods and Verification of Solution.
  * [23] L. Nielsen, E. Jørgensen, and S. Højsgaard.  Embedding a state space model into a Markov decision process.  Annals OR, 190:289–309, 10 2011.
  * [24] M. L. Puterman.  Markov Decision Processes: Discrete Stochastic Dynamic Programming.  John Wiley & Sons, Inc., USA, 1st edition, 1994.
  * [25] A. Rieder.  Inexact Newton regularization using conjugate gradients as inner iteration.  SIAM Journal on Numerical Analysis, 43(2):604–622, 2005.
  * [26] R. T. Rockafellar.  Monotone operators and the proximal point algorithm.  SIAM Journal on Control and Optimization, 14(5):877–898, 1976.
  * [27] Y. Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [28] Y. Saad and M. H. Schultz.  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [29] D. J. White.  A survey of applications of Markov decision processes.  Journal of the Operational Research Society, 44:1073–1096, 1993.
  * [30] R. Yaesoubi and T. Cohen.  Generalized Markov models of infectious disease spread: A novel framework for developing dynamic health policies.  European Journal of Operational Research, 215(3):679–687, 2011.