# An Analysis of Quantile Temporal-Difference Learning

Mark Rowland<EMAIL_ADDRESS>
DeepMindRémi Munos
DeepMind Mohammad Gheshlaghi Azar
DeepMind Yunhao Tang
DeepMind Georg Ostrovski
DeepMind Anna Harutyunyan
DeepMind Karl Tuyls
DeepMind Marc G. Bellemare
Google Research, Brain Team Will Dabney
DeepMind

###### Abstract

We analyse quantile temporal-difference learning (QTD), a distributional
reinforcement learning algorithm that has proven to be a key component in
several successful large-scale applications of reinforcement learning. Despite
these empirical successes, a theoretical understanding of QTD has proven
elusive until now. Unlike classical TD learning, which can be analysed with
standard stochastic approximation tools, QTD updates do not approximate
contraction mappings, are highly non-linear, and may have multiple fixed
points. The core result of this paper is a proof of convergence to the fixed
points of a related family of dynamic programming procedures with probability
1, putting QTD on firm theoretical footing. The proof establishes connections
between QTD and non-linear differential inclusions through stochastic
approximation theory and non-smooth analysis.

##  1 Introduction

In distributional reinforcement learning, an agent aims to predict the full
probability distribution over future returns it will encounter (Morimura et
al., 2010b, a; Bellemare et al., 2017, 2023), in contrast to predicting just
the mean return, as in classical reinforcement learning (Sutton and Barto,
2018). A widely-used family of algorithms for distributional reinforcement
learning is based on the notion of learning _quantiles_ of the return
distribution, an approach that originated with Dabney et al. (2018b), who
introduced the quantile temporal-difference (QTD) learning algorithm. This
approach has been particularly successful in combination with deep
reinforcement learning, and has been a central component in several recent
real-world applications of reinforcement learning, including sim-to-real
stratospheric balloon navigation (Bellemare et al., 2020), robotic
manipulation (Bodnar et al., 2020), and algorithm discovery (Fawzi et al.,
2022), as well as on benchmark simulated domains such as the Arcade Learning
Environment (Bellemare et al., 2013; Machado et al., 2018; Dabney et al.,
2018b, a; Yang et al., 2019) and racing simulation (Wurman et al., 2022).

Despite these empirical successes of QTD, little is known about its behaviour
from a theoretical viewpoint. In particular, questions regarding the
asymptotic behaviour of the algorithm – Do its predictions converge? Under
what conditions? What is the qualitative character of the predictions when
they do converge? – were left open. A core reason for this is that unlike
classical TD, and other distributional reinforcement learning algorithms such
as categorical temporal-difference learning (Rowland et al., 2018; Bellemare
et al., 2023), the updates of QTD rely on asymmetric $L^{1}$ losses. As a
result, these updates do not approximate the application of a contraction
mapping, are highly non-linear (even in the tabular setting), and also may
have multiple fixed points (depending on the exact structure of the reward
distributions of the environment, as we will see), and their analysis requires
a distinct set of tools to those typically used to analyse temporal-difference
learning algorithms.

In this paper, we prove the convergence of QTD – notably under weaker
assumptions than are required in typical proofs of convergence for classical
TD learning – establishing it as a sound algorithm with theoretical
convergence guarantees, and paving the way for further analysis and
investigation. The more general conditions stem from the structure of the QTD
updates (namely, their boundedness), and the proof is obtained through the use
of stochastic approximation theory with differential inclusions, which have
seldom been used in the context of RL theory previously. The remainder of the
paper is structured as follows.

We begin by providing background on Markov decision processes, classical TD
learning, and quantile regression in Section 2. After motivating the QTD
algorithm in Section 3, we describe the related family of quantile dynamic
programming (QDP) algorithms, and provide a convergence analysis of these
algorithms in Section 4. We then present the main result, a convergence
analysis of QTD, in Section 5. The proof relies on the stochastic
approximation framework set out by Benaïm et al. (2005), arguing that the QTD
algorithm approximates a continuous-time differential inclusion, and then
constructing a Lyapunov function to demonstrate that the limiting behaviour of
trajectories of the differential inclusion matches that of the QDP algorithms
introduced earlier. Finally, in Section 6, we analyse the limit points of QTD,
bounding their approximation error to the true return distributions of
interest, and investigating the kinds of approximation artefacts that arise
empirically.

##  2 Background

We first introduce background concepts and notation.

###  2.1 Markov decision processes

We consider a Markov decision process specified by finite state and action
spaces $\mathcal{X}$ and $\mathcal{A}$, transition kernel
$P_{\mathcal{X}}:\mathcal{X}\times\mathcal{A}\rightarrow\mathscr{P}(\mathcal{X})$,
reward distribution function
$P_{\mathcal{R}}:\mathcal{X}\times\mathcal{A}\rightarrow\mathscr{P}_{1}(\mathbb{R})$,
and discount factor $\gamma\in[0,1)$. Here, $\mathscr{P}(\mathcal{X})$ is the
set of probability distributions over the finite set $\mathcal{X}$, and
$\mathscr{P}_{1}(\mathbb{R})$ is the set of probability distributions over
$\mathbb{R}$ (with its usual Borel $\sigma$-algebra) with finite mean.

Given a policy $\pi:\mathcal{X}\rightarrow\mathscr{P}(\mathcal{A})$ and an
initial state $x_{0}\in\mathcal{X}$, an agent interacting with the environment
using the policy $\pi$ generates a sequence of states, actions and rewards
$(X_{t},A_{t},R_{t})_{t=0}^{\infty}$, called a _trajectory_ , the joint
distribution of which is determined by the transition dynamics and reward
distributions of the environment, and the policy of the agent. More precisely,
we have

  * •

$X_{0}=x_{0}$, and for each $t\geq 0$:

  * •

$A_{t}\mid X_{0:t},A_{0:t-1},R_{0:t-1}\sim\pi(\cdot|X_{t})$;

  * •

$R_{t}\mid X_{0:t},A_{0:t},R_{0:t-1}\sim P_{\mathcal{R}}(\cdot|X_{t},A_{t})$;

  * •

$X_{t+1}\mid X_{0:t},A_{0:t},R_{0:t}\sim P_{\mathcal{X}}(\cdot|X_{t},A_{t})$.

The distribution of the trajectory is thus parametrised by the initial state
$x_{0}$, and the policy $\pi$. To illustrate this dependency, we use the
notation $\mathbb{P}^{\pi}_{x_{0}}$ and $\mathbb{E}^{\pi}_{x_{0}}$ to denote
the probability distribution and expectation operator corresponding to this
distribution, and will write $P^{\pi}(\cdot|x)$ for the joint distribution
over a reward–next-state pair when the current state is $x$.

###  2.2 Predicting expected returns and the return distribution

The quality of the agent’s performance on the trajectory is quantified by the
_discounted return_ , or simply the _return_ , given by

$\displaystyle\sum_{t=0}^{\infty}\gamma^{t}R_{t}\,.$ (1)

The return is a random variable, whose sources of randomness are the random
selections of actions made according to $\pi$, the randomness in state
transitions, and the randomness in rewards observed. Typically in
reinforcement learning, a single scalar summary of performance is given by the
expectation of this return over all these sources of randomness. For a given
policy, this is summarised across each possible starting state via the _value
function_ $V^{\pi}:\mathcal{X}\rightarrow\mathbb{R}$, defined by

$\displaystyle
V^{\pi}(x)=\mathbb{E}^{\pi}_{x}\Bigg{[}\sum_{t=0}^{\infty}\gamma^{t}R_{t}\Bigg{]}\,.$
(2)

Learning the value function of a policy $\pi$ from sampled trajectories
generated through interaction with the environment is a central problem in
reinforcement learning, referred to as the _policy evaluation task_.

Each expected return is a scalar summary of a much more rich, complex object:
the probability distributions of the random return in Equation (1) itself.
_Distributional reinforcement learning_ (Bellemare et al., 2023) is concerned
with the problem of learning to predict the _probability distribution_ over
returns, in contrast to just their expected value. Mathematically, the goal is
to learn the return-distribution function
$\eta^{\pi}:\mathcal{X}\rightarrow\mathscr{P}(\mathbb{R})$; for each state
$x\in\mathcal{X}$, $\eta^{\pi}(x)$ is the probability distribution of the
random return in Expression (1) when the trajectory begins at state $x$, and
the agent acts using policy $\pi$. Mathematically, we have

$\displaystyle\eta^{\pi}(x)=\mathcal{D}^{\pi}_{x}\Bigg{(}\sum_{t\geq
0}\gamma^{t}R_{t}\Bigg{)}\,,$

where $\mathcal{D}^{\pi}_{x}$ extract the probability distribution of a random
variable under $\mathbb{P}^{\pi}_{x}$.

There are several distinct motivations for aiming to learn these more complex
objects. First, the richness of the distribution provides an abundance of
signal for an agent to learn from, in contrast to a single scalar expectation.
The strong performance of deep reinforcement learning agents that incorporate
distributional predictions is hypothesised to be related to this fact (Dabney
et al., 2018b; Barth-Maron et al., 2018; Dabney et al., 2018a; Yang et al.,
2019). Second, learning about the full probability distribution of returns
makes possible the use of _risk-sensitive_ performance criteria; one may be
interested in not only the expected return under a policy, but also the
variance of the return, or the probability of the return being under a certain
threshold.

Unlike the value function $V^{\pi}$, which is an element of
$\mathbb{R}^{\mathcal{X}}$, and can therefore be straightforwardly represented
on a computer (up to floating-point precision), the return-distribution
function $\eta^{\pi}$ is not representable. Each object $\eta^{\pi}(x)$ is a
probability distribution over the real numbers, and, informally speaking,
probability distributions have infinitely many degrees of freedom.
Distributional reinforcement learning algorithms therefore typically work with
a subset of distributions that _are_ amenable to parametrisation on a computer
(Bellemare et al., 2023). Common choices of subsets include categorical
distributions (Bellemare et al., 2017), exponential families (Morimura et al.,
2010b), and mixtures of Gaussian distributions (Barth-Maron et al., 2018).
Quantile temporal-difference learning, the core algorithm of study in this
paper, aims to learn a particular set of quantiles of the return distribution,
as described in Section 3.

###  2.3 Monte Carlo and temporal-difference learning

To foreshadow our description and motivation of quantile temporal-difference
learning, we recall a line of thinking that interprets the classical TD
learning update rule as an approximation to Monte Carlo learning; this
material is common to many introductory texts on reinforcement learning
(Sutton and Barto, 2018), and we present it here to make a direct analogy with
QTD. First, we may observe that, under the condition that all reward
distributions have finite variance, $V^{\pi}(x)$ is the unique minimiser of
the following loss function over $v$:

$\displaystyle\mathcal{L}^{\pi}_{x}(v)=\frac{1}{2}\mathbb{E}^{\pi}_{x}\Big{[}\Big{(}v-\sum_{t=0}^{\infty}\gamma^{t}R_{t}\Big{)}^{2}\Big{]}\,.$

This well-known characterisation of the expectation of a random variable is
readily verified by, for example, observing that the loss is convex and
differentiable in $v$, and solving the equation
$\partial_{v}\mathcal{L}^{\pi}_{x}(v)=0$. This motivates an approach to
learning $V^{\pi}(x)$ based on stochastic gradient descent on the loss
function $\mathcal{L}^{\pi}_{x}$. We maintain an estimate
$V\in\mathbb{R}^{\mathcal{X}}$ of the value function, and each time a
trajectory $(X_{t},A_{t},R_{t})_{t\geq 0}$ beginning at state $x$ is observed,
we can obtain an unbiased estimator of the negative gradient of
$\mathcal{L}^{\pi}_{x}(V(x))$ as

$\displaystyle\sum_{t=0}^{\infty}\gamma^{t}R_{t}-V(x)\,,$

and update $V(x)$ by taking a step in the direction of this negative gradient,
with some step size $\alpha$:

$\displaystyle V(x)\leftarrow
V(x)+\alpha\Big{(}\sum_{t=0}^{\infty}\gamma^{t}R_{t}-V(x)\Big{)}\,.$ (3)

This is a _Monte Carlo_ algorithm, so called because it uses Monte Carlo
samples of the random return to update the estimate $V$.

A popular alternative to this Monte Carlo algorithm is temporal-difference
learning, which replaces samples from the random return with a _bootstrapped_
approximation to the return, obtained from a transition $(x,A,R,X^{\prime})$
by combining the immediate reward $R$ with the current estimate of the
expected return obtained at $X^{\prime}$, resulting in the return estimate

$\displaystyle R+\gamma V(X^{\prime})\,,$ (4)

and the corresponding update rule

$\displaystyle V(x)\leftarrow V(x)+\alpha(R+\gamma V(X^{\prime})-V(x))\,.$ (5)

While the mean-return estimator in Expression (4) is generally _biased_ ,
since $V(X^{\prime})$ is not generally equal to the true expected return
$V^{\pi}(X^{\prime})$, it is often a lower-variance estimate, since we are
replacing the _random_ return from $X^{\prime}$ with an estimate of its
expectation (Sutton, 1988; Sutton and Barto, 2018; Kearns and Singh, 2000).

This motivates the TD learning rule in Expression (5) based on the Monte Carlo
update rule in Expression (3), with the understanding that this algorithm can
be applied more generally, with access only to sampled transitions (rather
than full trajectories), and may result in more accurate estimates of the
value function, due to lower-variance updates, and the propensity of TD
algorithms to “share information” across states. Note however that this does
not _prove_ anything about the behaviour of temporal-difference learning, and
a fully rigorous theory of the asymptotic behaviour emerged several years
after TD methods were formally introduced (Sutton, 1984, 1988; Watkins, 1989;
Watkins and Dayan, 1992; Dayan, 1992; Dayan and Sejnowski, 1994; Jaakkola et
al., 1994; Tsitsiklis, 1994).

##  3 Quantile temporal-difference learning and quantile dynamic programming

We now present the main algorithms of study in this paper.

###  3.1 Quantile regression

To motivate QTD, we begin by considering how we might adapt a Monte Carlo
algorithm such as that in Expression (3) to learn about the distribution of
returns, rather than just their expected value. We cannot learn the return
distribution in its entirety with a finite collection of parameters; the space
of return distributions is infinite-dimensional, so we must instead be
satisfied with learning an approximation of the return distribution by
selecting a _probability distribution representation_ (Bellemare et al., 2023,
Chapter 5) – a subset of probability distributions parametrised by a finite-
dimensional set of parameters. The approach of quantile temporal-difference
learning is to learn an approximation of the form

$\displaystyle\eta(x)=\sum_{i=1}^{m}\frac{1}{m}\delta_{\theta(x,i)}\,;$ (6)

an equally-weighted mixture of Dirac deltas, for each state $x\in\mathcal{X}$.
The quantile-based approach to distributional reinforcement learning aims to
have the particle locations $(\theta(x,i))_{i=1}^{m}$ approximate certain
_quantiles_ of $\eta^{\pi}(x)$.

######  Definition 3.1.

For a probability distribution $\nu\in\mathscr{P}(\mathbb{R})$ and parameter
$\tau\in(0,1)$, the set of $\tau$-quantiles of $\nu$ is given by the set

$\displaystyle\\{z\in\mathbb{R}:F_{\nu}(z)=\tau\\}\cup\inf\\{y\in\mathbb{R}:F_{\nu}(y)>\tau\\}\,,$

where $F_{\nu}:\mathbb{R}\rightarrow[0,1]$ is the CDF of $\nu$, defined by
$F_{\nu}(t)=\mathbb{P}_{Z\sim\nu}(Z\leq t)$ for all $t\in\mathbb{R}$.

Expanding on this definition, if the set
$\\{z\in\mathbb{R}:F_{\nu}(z)=\tau\\}$ is non-empty, then the $\tau$-quantiles
are precisely the values $z$ such that $\mathbb{P}_{Z\sim\nu}(Z\leq z)=\tau$.
If however this set is empty (which may arise when $F_{\nu}$ has points of
discontinuity), then the quantile is the smallest value $y$ such that
$\mathbb{P}_{Z\sim\nu}(Z\leq y)>\tau$. Note also that if $F_{\nu}$ is strictly
increasing, this guarantees uniqueness of each $\tau$-quantile for
$\tau\in(0,1)$; this is often a useful property in the analysis we consider
later. These different cases are illustrated in Figure 1. The generalised
inverse CDF of $\nu$, $F^{-1}_{\nu}:(0,1)\rightarrow\mathbb{R}$, is defined by

$\displaystyle F^{-1}_{\nu}(\tau)=\inf\\{y:F_{\nu}(y)\geq\tau\\}\,,$

and provides a way of uniquely specifying a quantile for each level $\tau$. In
cases where there is not a unique $\tau$-quantile (see Figure 1),
$F^{-1}_{\nu}(\tau)$ corresponds to the _left-most_ or _least_ valid
$\tau$-quantile. We also introduce the notation

$\displaystyle\bar{F}^{-1}_{\nu}(\tau)=\inf\\{y:F_{\nu}(y)>\tau\\}\,,$

which corresponds to the _right-most_ or _greatest_ $\tau$-quantile; notice
the strict inequality that appears in the definition, in contrast to that of
$F^{-1}_{\nu}(\tau)$. If $F^{-1}_{\nu}$ is continuous at $\tau$, then
$F^{-1}_{\nu}(\tau)=\bar{F}^{-1}_{\nu}(\tau)$, as is the case for
$\tau=\tau_{1}$ and $\tau=\tau_{3}$ in Figure 1. However, if $F_{\nu}$ has a
flat region for the value $\tau$ (as is the case for $\tau=\tau_{2}$ in Figure
1), then $F^{-1}_{\nu}(\tau)$ and $\bar{F}^{-1}_{\nu}(\tau)$ are distinct, and
correspond to the boundary points of this flat region.

Figure 1: The three distinct scenarios that arise in defining quantiles.
Firstly, there is a value $z_{1}$ for which $F_{\nu}(z_{1})=\tau_{1}$ and at
which $F_{\nu}$ is strictly increasing. Therefore $z_{1}$ is the unique
$\tau_{1}$-quantile of $\nu$. Next, there is an interval
$[z_{2},z_{2}^{\prime}]$ on which $F_{\nu}$ equals $\tau_{2}$, therefore all
elements in this interval are $\tau_{2}$-quantiles of $\nu$. Finally, there is
no value $z$ such that $F_{\nu}(z)=\tau_{3}$, and the unique
$\tau_{3}$-quantile is therefore defined by the infimum part of the
definition.

Algorithmically, we aim for $\theta(x,i)$ to approximate a $\tau_{i}$-quantile
of $\eta^{\pi}(x)$, where $\tau_{i}=\nicefrac{{2i-1}}{{2m}}$. To build
learning algorithms that achieve this, we require an incremental algorithm
that updates $\theta(x,i)$ in response to samples from the target distribution
$\eta^{\pi}(x)$, which converges to a $\nicefrac{{2i-1}}{{2m}}$-quantile of
$\eta^{\pi}(x)$.

Such an approach is available by using the quantile regression loss. We define
the quantile regression loss associated with distribution
$\nu\in\mathscr{P}(\mathbb{R})$ and quantile level $\tau\in(0,1)$ as a
function of $v$ by

$\displaystyle\mathbb{E}_{Z\sim\nu}[(\tau\mathbbm{1}\\{Z>v\\}+(1-\tau)\mathbbm{1}\\{Z<v\\})|Z-v|]\,.$
(7)

This loss is the expectation of an asymmetric absolute value loss, in which
positive and negative errors are weighted according to the parameters $\tau$
and $1-\tau$ respectively. Just as the expected squared loss encountered above
encodes the mean as its unique minimiser, the quantile regression loss encodes
the $\tau$-quantiles of $\nu$ as the unique minimisers; see e.g. Koenker
(2005) for further background. Thus, applying the quantile regression loss to
the problem of estimating $\tau$-quantiles of the return distribution, we
arrive at the loss

$\displaystyle\mathcal{L}^{\tau,\pi}_{x}(v)=\mathbb{E}^{\pi}_{x}\Big{[}(\tau\mathbbm{1}\\{\Delta>0\\}+(1-\tau)\mathbbm{1}\\{\Delta<0\\})|\Delta|\Big{]}\,,\quad\text{where
}\Delta=\sum_{t=0}^{\infty}\gamma^{t}R_{t}-v\,.$

Given an observed return $\sum_{t\geq 0}\gamma^{t}R_{t}$ from the state $x$,
we therefore have that an unbiased estimator of the negative
gradient111Technically speaking, we are assuming that differentiation and
expectation can be interchanged here. Further, under certain circumstances the
loss is only _sub_ -differentiable. As our principal goal in this section is
to provide intuition for QTD, we do not comment further on these technical
details. The convergence results later in the paper deal with these issues
carefully. of this loss is

$\displaystyle\tau\mathbbm{1}\Big{\\{}\sum_{t=0}^{\infty}\gamma^{t}R_{t}>v\Big{\\}}-(1-\tau)\mathbbm{1}\Big{\\{}\sum_{t=0}^{\infty}\gamma^{t}R_{t}<v\Big{\\}}\,,$

which motivates an update rule of the form

$\displaystyle\theta(x,i)\leftarrow\theta(x,i)+\alpha\Big{(}\tau_{i}\mathbbm{1}\Big{\\{}\sum_{t=0}^{\infty}\gamma^{t}R_{t}>\theta(x,i)\Big{\\}}-(1-\tau_{i})\mathbbm{1}\Big{\\{}\sum_{t=0}^{\infty}\gamma^{t}R_{t}<\theta(x,i)\Big{\\}}\Big{)}\,.$
(8)

This can be rewritten as

$\displaystyle\theta(x,i)\leftarrow\theta(x,i)+\alpha\Big{(}\tau_{i}-\mathbbm{1}\Big{\\{}\sum_{t=0}^{\infty}\gamma^{t}R_{t}<\theta(x,i)\Big{\\}}\Big{)}\,.$
(9)

This is essentially the application of the stochastic gradient descent method
for quantile regression to learning quantiles of the return distribution.

###  3.2 Quantile temporal-difference learning

We can motivate and describe the quantile temporal-difference learning
algorithm (Dabney et al., 2018b; Bellemare et al., 2023) by modifying the
Monte Carlo algorithm in Expression (8) in a similar manner to the
modification that led to the TD algorithm in Expression (5). We replace the
Monte Carlo return

$\displaystyle\sum_{t=0}^{\infty}\gamma^{t}R_{t}$

based on a full trajectory, with an approximate sample from the return
distribution derived from an observed transition $(x,R,X^{\prime})$, and the
estimate $\eta(X^{\prime})$ of the return distribution at state $X^{\prime}$.
If the return distribution estimate $\eta(X^{\prime})$ takes the form given in
Equation (6), as is the case for the probability distribution representation
considered here, then such a sample return is obtained as

$\displaystyle R+\gamma\theta(X^{\prime},J)\,,$

with $J$ sampled uniformly from $\\{1,\ldots,m\\}$. This yields the update
rule

$\displaystyle\theta(x,i)\leftarrow\theta(x,i)+\alpha\Big{(}\tau_{i}-\mathbbm{1}\Big{\\{}R+\gamma\theta(X^{\prime},J)<\theta(x,i)\Big{\\}}\Big{)}\,.$

We can consider also a variance-reduced version of this update, in which we
average over updates performed under different realisations of $J$, leading to
the update

$\displaystyle\theta(x,i)\leftarrow\theta(x,i)+\frac{\alpha}{m}\sum_{j=1}^{m}\Big{(}\tau_{i}-\mathbbm{1}\Big{\\{}R+\gamma\theta(X^{\prime},j)<\theta(x,i)\Big{\\}}\Big{)}\,.$
(10)

This is precisely the quantile temporal-difference learning update, presented
in Algorithm 1 below, which underlies many recent successful applications of
reinforcement learning at scale (Dabney et al., 2018b, a; Yang et al., 2019;
Bellemare et al., 2020; Wurman et al., 2022; Fawzi et al., 2022). Similar to
other temporal-difference learning algorithms, QTD updates its parameters
$((\theta(x,i))_{i=1}^{m}:x\in\mathcal{X})$ on the basis of sample transitions
$(x,r,x^{\prime})$ generated through interaction with the environment via the
policy $\pi$, comprising a state, reward, and next state.

0: Quantile estimates $\theta\in\mathbb{R}^{\mathcal{X}\times[m]}$ , Observed
transition $(x,r,x^{\prime})$ , Learning rate $\alpha$.

1: Set $\tau_{i}=\tfrac{2i-1}{2m}$ for each $i=1,\ldots,m$.

2: for $i=1,\ldots,m$ do

3: Set
$\theta^{\prime}(x,i)\leftarrow\theta(x,i)+\alpha\frac{1}{m}\sum_{j=1}^{m}\left[\tau_{i}-\mathbbm{1}\Big{\\{}r+\gamma\theta(x^{\prime},j)-\theta(x,i)<0\Big{\\}}\right]$

4: end for

5: for $i=1,\ldots,m$ do

6: Set $\theta(x,i)\leftarrow\theta^{\prime}(x,i)$

7: end for

8: return $((\theta^{\prime}(x,i))_{i=1}^{m}:x\in\mathcal{X})$

Algorithm 1 QTD update

Whilst the QTD update makes use of temporal-difference errors
$r+\gamma\theta(x^{\prime},j)-\theta(x,i)$, there are two key differences to
the use of analogous quantities in classical TD learning. First, the TD errors
influence the update only through their sign, not their magnitude. Second, the
predictions at each state $(\theta(x,i))_{i=1}^{m}$ are indexed by $i$, and
each update includes a distinct term $\tau_{i}$ (equal to
$\nicefrac{{2i-1}}{{2m}}$). The presence of these terms causes the learnt
parameters to make distinct predictions, as described in Section 3.1.
Practical implementations of QTD use these precise values for $\tau_{i}$,
equally spaced out on $[0,1]$, as proposed by Dabney et al. (2018b). Much of
the analysis in this paper goes through straightforwardly for other values of
$\tau_{i}$, though we will see in Section 6 that this choice is well motivated
in that it provides the best bounds on distribution approximation.

The discussion above provides _motivation_ for the form of the QTD update
given in Algorithm 1, and intuition as to why this algorithm might perform
reasonably, and learn a sensible approximation to the return distribution.
However, it stops short of providing an explanation of how the algorithm
should be expected to behave, or providing any theoretical guarantees as to
what the algorithm will in fact converge to. A core goal of the sections that
follow is to answer these questions, and put QTD on firm theoretical footing.

###  3.3 Motivating examples

Before undertaking an analysis of QTD, we pause to provide several numerical
examples of its behaviour in small environments. These examples provide
further intuition for the characteristics of the algorithm, illustrate the
breadth of qualitative behaviours it can exhibit, and provide motivation for
the kinds of theoretical questions we might hope to answer.

######  Example 3.2.

Consider the chain MDP illustrated at the top of Figure 2. The random return
at each state is a sum of independent Gaussian random variables, and hence the
return distribution at each state is Gaussian. The centre plot in Figure 2
illustrates the evolution of $m=5$ quantile estimates learnt by QTD, using a
constant learning rate of 0.01, and updating all states at each update. The
estimates eventually settle after around 5,000 updates, with small
oscillations around this point. The bottom of Figure 2 compares the true
return distribution at each state (in blue), with the approximation learnt by
QTD (in black), and the approximation obtained with the true value of the five
quantiles of interest (grey). The behaviour of QTD in this case raises several
questions: Can it be shown that QTD is guaranteed to stabilise/converge around
a certain point? Can a guarantee be given on the quality of the approximate
distributions learnt by QTD?

Figure 2:  Top: A chain MDP with four states. Each transition yields a reward
distributed according to $N(0,1)$; from $x_{3}$, the episode ends. The
discount factor is $\gamma=0.9$. Middle: The progress of QTD, run with $m=5$
quantiles, over the course of 10,000 updates. Bottom: True CDF of the return
distribution (blue) at each state, along with the final estimate produced by
QTD (black), and the approximation produced by the quantiles of the return
distribution (grey).

######  Example 3.3.

For a different perspective on the behaviour of QTD, consider a two-state MDP
with transition dynamics as illustrated in the top-left of Figure 3, and
discount factor $\gamma=0.5$. The reward obtained when transitioning from
state $x_{1}$ is distributed as N(2, 1), and the reward obtained when
transitioning from state $x_{2}$ is distributed as N(-1, 1); here, we write
$N(\mu,\sigma^{2})$ for the normal distribution with mean $\mu$ and variance
$\sigma^{2}$. We consider the case of learning a single quantile (the median)
at each of these two states, taking $m=1$; this allows us to plot the full
phase space of the QTD algorithm in a two-dimensional plot.

The top-right of Figure 3 shows a path taken by QTD under this MDP. In
addition, the streamplot illustrates the direction of the _expected_ update
that QTD undertakes at each point in phase space. We empirically observe
convergence of the algorithm to a point. Additionally, the expected update
direction changes smoothly; the result is a vector field that appears to point
towards the point of convergence from all directions.

The bottom-left of Figure 3 shows a path taken by QTD under a modified version
of the MDP, in which the reward distributions N(2, 1) and N(-1, 1) are
replaced with $\delta_{2}$ and $\delta_{-1}$, respectively. We observe that
the algorithm still converges to a point, although the vector field of
expected update directions is now piecewise constant, with discontinuities
along several lines. This behaviour is typical of QTD; the less ‘smooth’ the
reward distributions in the MDP, the more abrupt the changes in behaviour we
typically observe with QTD.

Finally, we consider a modified version of the MDP in which all transition
probabilities are $\nicefrac{{1}}{{2}}$, rewards from state $x_{1}$ are always
$2$, and rewards from state $x_{2}$ are always $-1$. In this case, QTD no
longer appears to converge to a point, but instead converges to the set
bounded by the four grey lines appearing in the bottom-right of Figure 3, and
subsequently performing a random walk over this set. This collection of
examples illustrates that QTD can exhibit a fairly wide family of behaviours
depending on the characteristics of the environment.

Figure 3:  Top left: The example Markov decision process described in Example
3.3. Top right: Example dynamics of QTD with $m=1$ in this environment, when
reward distributions are Gaussian. Also included are the directions of
expected update, in blue. Bottom left: Example dynamics and expected update
directions when reward distributions are Dirac deltas. Bottom right: Example
dynamics and expected updates with modified environment transition
probabilities.

###  3.4 Quantile dynamic programming

Recall the QTD update given in Equation (10). As described in Section 3.2,
this update serves, on average, to move $\theta(x,i)$ in the direction of the
$\tau_{i}$-quantiles of the distribution of the random variable
$R+\theta(X^{\prime},J)$, where $(x,R,X^{\prime})$ is a random transition
generated by interacting with the environment using $\pi$, and
$J\sim\text{Unif}(\\{1,\ldots,m\\})$.

Suppose we were able to update $\theta(x,i)$ not just with a single gradient
step in this direction, but instead were able to update it to take on exactly
this quantile value. This motivates a dynamic programming alternative to QTD,
_quantile dynamic programming_ (QDP), which directly calculates these
quantiles iteratively, in a similar manner to iterative policy evaluation in
classical reinforcement learning (Bertsekas and Tsitsiklis, 1996).

The mathematical structure of such an algorithm is given in Algorithm 2. This
stops short of being an implementable algorithm, since we do not describe in
what format the transition probabilities and reward distributions are
available, which are required to evaluate the inverse CDFs that arise in the
algorithm. However, for MDPs in which transition probabilities and reward
distributions are available, QDP is an algorithmic framework of interest in
its own right, and to this end we provide several concrete implementations in
Appendix B.

The QDP template in Algorithm 2 is parametrised by the interpolation
parameters $\lambda\in[0,1]^{\mathcal{X}\times[m]}$. These parameters control
exactly which quantile is chosen when the desired quantile level $\tau_{i}$
corresponds to a flat region of the CDF for the distribution $\nu$ (the second
case in Figure 1). QDP was originally presented by Bellemare et al. (2023) in
the case $\lambda(x,i)\equiv 0$; the presentation here generalises QDP to a
family of algorithms, parametrised by $\lambda$.

Our interest in QDP stems from the fact that QTD can be viewed as
approximating the behaviour of the QDP algorithms, without requiring access to
the transition structure and reward distributions of the environment. In
particular, we will show that under appropriate conditions, the asymptotic
behaviour of QTD and QDP are equivalent: they both converge to the same
limiting points. Figure 4 illustrates the behaviour of the QDP algorithm in
the environment described in Example 3.3; since the reward distributions in
this example have strictly increasing CDFs, QDP behaves identically for all
choices of interpolation parameters $\lambda$. QTD and QDP appear to have the
same asymptotic behaviour, converging to the same limiting point. In cases
where QTD appears to converge to a set, such as in the bottom-right plot of
Figure 3.3, the relationship is slightly more complicated, and there is a
correspondence between the asymptotic behaviour of QTD and the family of
dynamic programming algorithms parametrised by $\lambda$, as illustrated at
the bottom of Figure 4. Thus, to understand the asymptotic behaviour of QTD,
we begin by analysing the asymptotic behaviour of QDP.

Figure 4:  Top left: Illustration of QDP (dashed purple) and QTD (solid red)
on the first MDP from Example 3.3, with Gaussian rewards. Top right:
Illustration of QDP and QTD on the first MDP from Example 3.3, with
deterministic rewards. Bottom: Values of $\lambda$ and corresponding fixed
points of QDP in the final MDP from Example 3.3.

0: Quantile estimates $((\theta(x,i))_{i=1}^{m}:x\in\mathcal{X})$,
Interpolation parameters $\lambda\in[0,1]^{\mathcal{X}\times[m]}$.

1: for $x\in\mathcal{X}$ do

2: Let $(x,R,X^{\prime})$ be a random transition under $\pi$, and
$J\sim\text{Unif}(\\{1,\ldots,m\\})$.

3: Set $\nu$ to be the distribution of $R+\gamma\theta(X^{\prime},J)$.

4: for $i=1,\ldots,m$ do

5: Set
$\theta(x,i)\leftarrow(1-\lambda(x,i))F^{-1}_{\nu}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\nu}(\tau_{i})$.

6: end for

7: end for

8: return $((\theta(x,i))_{i=1}^{m}:x\in\mathcal{X})$

Algorithm 2 Quantile dynamic programming

##  4 Convergence of quantile dynamic programming

We can decompose the update QDP performs into the composition of several
operators. Algorithm 2 manipulates tables of the form
$((\theta(x,i))_{i=1}^{m}:x\in\mathcal{X})$. For a given state $x$, the vector
$(\theta(x,i))_{i=1}^{m}$ represents the estimated
$\nicefrac{{2i-1}}{{2m}}$-quantiles of the return distribution at state $x$,
for $i=1,\ldots,m$. In mathematically analysing the algorithm, it is useful to
be able to refer to the distribution encoded by these quantiles:

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\delta_{\theta(x,i)}\,,$ (11)

and reason about the transformations undertaken by Algorithm 2 directly in
terms of distributions. To this end, if we write
$\eta(x)\in\mathscr{P}(\mathbb{R})$ for the probability distribution
associated with the quantile estimates $(\theta(x,i))_{i=1}^{m}$, we can
interpret the transformation performed by Algorithm 2 as comprising two parts,
which we now describe in turn.

First, the variable $\eta(x)$ is assigned the distribution of $R+\gamma
G(X^{\prime})$, where $R,X^{\prime}$ are the random reward and next-state
encountered from the initial state $x$ with policy $\pi$, and
$(G(y):y\in\mathcal{X})$ is an independent collection of random variables,
with each $G(y)$ distributed according to $\eta(y)$.

We write
$\mathcal{T}^{\pi}:\mathscr{P}(\mathbb{R})^{\mathcal{X}}\rightarrow\mathscr{P}(\mathbb{R})^{\mathcal{X}}$
for this transformation. The function $\mathcal{T}^{\pi}$ is known as the
_distributional Bellman operator_ (Bellemare et al., 2017; Rowland et al.,
2018; Bellemare et al., 2023). In terms of the above definition via
distributions of random variables, $\mathcal{T}^{\pi}$ can be written

$\displaystyle(\mathcal{T}^{\pi}\eta)(x)=\mathcal{D}_{\pi}(R+\gamma
G(X^{\prime}))\,,$

where $(x,R,X^{\prime})$ is a random environment transition beginning at $x$,
independent of $(G(y):y\in\mathcal{X})$, and $\mathcal{D}_{\pi}$ extracts the
distribution of its argument when $(x,R,X^{\prime})$ is generated by sampling
an action from $\pi$. See Bellemare et al. (2023) for further background on
the distributional Bellman operator.

In general, $\mathcal{T}^{\pi}\eta$ may comprise much more complicated
distributions than $\eta$ itself, with many more atoms, or possibly infinite
support, if reward distributions are infinitely-supported. Algorithm 2 does
not return these full transformed distributions, but rather approximations, or
_projections_ , of these distributions, obtained by keeping only information
about certain quantiles (in the inner for-loop of Algorithm 2); this is the
second distribution transformation the algorithm undertakes. Each choice of
interpolation parameters $\lambda$ corresponds to a different projection
operator, denoted
$\Pi^{\lambda}:\mathscr{P}(\mathbb{R})^{\mathcal{X}}\rightarrow\mathscr{P}(\mathbb{R})^{\mathcal{X}}$,
and defined by

$\displaystyle(\Pi^{\lambda}\eta)(x)=\frac{1}{m}\sum_{i=1}^{m}\delta_{(1-\lambda(x,i))F^{-1}_{\eta(x)}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\eta(x)}(\tau_{i})}\,.$
(12)

Thus, the composition $\Pi^{\lambda}\mathcal{T}^{\pi}$ – the projected
distributional Bellman operator – is a transformation on the space of return-
distribution functions $\mathscr{P}(\mathbb{R})^{\mathcal{X}}$. We will also
find it useful to abuse notation slightly and consider
$\Pi^{\lambda}\mathcal{T}^{\pi}$ as an operator on the space
$\mathbb{R}^{\mathcal{X}\times[m]}$ of parameters that QDP and QTD operate
over. The understanding is that an input
$\theta\in\mathbb{R}^{\mathcal{X}\times[m]}$ is first re-interpreted as a
collection of distributions as in Expression (11), with
$\Pi^{\lambda}\mathcal{T}^{\pi}$ applied as defined above to this collection
of probability distributions, and then finally extracting the support of the
output distributions, which take the form

$\displaystyle\sum_{i=1}^{m}\frac{1}{m}\delta_{z_{i}}\,,$

to return an element of $\mathbb{R}^{\mathcal{X}\times[m]}$. We will also
write $\mathcal{T}^{\pi}\theta$ for the element of
$\mathscr{P}(\mathbb{R})^{\mathcal{X}}$ obtained by applying
$\mathcal{T}^{\pi}$ to the distributions $(\eta(x):x\in\mathcal{X})$ defined
by

$\displaystyle\eta(x)=\sum_{i=1}^{m}\frac{1}{m}\delta_{\theta(x,i)}\,.$

######  Remark 4.1.

This convention highlights that there are two complementary views of
distributional reinforcement learning algorithms, through finite-dimensional
sets of parameters, and through probability distributions. The view in terms
of probability distributions is often useful in contraction analysis, and in
measuring approximation error, while we will see that the parameter view is
key to the stochastic approximation analysis that follows, and is ultimately
the way in which these algorithms are implemented.

With this convention, $\Pi^{\lambda}\mathcal{T}^{\pi}\theta$ is precisely the
table $\theta^{\prime}$ output by Algorithm 2 on input $\theta$, and so the
QDP algorithm is mathematically equivalent to repeated application of the
operator $\Pi^{\lambda}\mathcal{T}^{\pi}$ to an initial collection of quantile
estimates. To understand the long-term behaviour of QDP, we can therefore seek
to understand this projected operator $\Pi^{\lambda}\mathcal{T}^{\pi}$.

###  4.1 Convergence analysis

We will show that $\Pi^{\lambda}\mathcal{T}^{\pi}$ is a contraction mapping
with respect to an appropriate metric over return-distribution functions.
Building on the analysis in the case of $\lambda\equiv 0$ carried out by
Dabney et al. (2018b) and Bellemare et al. (2023), we use the
Wasserstein-$\infty$ metric
$w_{\infty}:\mathscr{P}(\mathbb{R})\times\mathscr{P}(\mathbb{R})\rightarrow[0,\infty]$,
defined by

$\displaystyle
w_{\infty}(\nu,\nu^{\prime})=\sup_{t\in(0,1)}|F^{-1}_{\nu}(t)-F^{-1}_{\nu^{\prime}}(t)|\,,$

and its extension to return-distribution functions,
$\bar{w}_{\infty}:\mathscr{P}(\mathbb{R})^{\mathcal{X}}\times\mathscr{P}(\mathbb{R})^{\mathcal{X}}\rightarrow[0,\infty]$,
given by

$\displaystyle\bar{w}_{\infty}(\eta,\eta^{\prime})=\max_{x\in\mathcal{X}}\sup_{t\in(0,1)}|F^{-1}_{\eta(x)}(t)-F^{-1}_{\eta^{\prime}(x)}(t)|\,.$

Both $w_{\infty}$ and $\bar{w}_{\infty}$ fulfil all the requirements of a
metric, except that they may assign infinite distances (Villani, 2009; see
also Bellemare et al. (2023) for a detailed discussion specifically in the
context of reinforcement learning). We must therefore take some care as to
when distances are finite. The following is established by Bellemare et al.
(2023, Proposition 4.15).

######  Proposition 4.2.

The distributional Bellman operator
$\mathcal{T}^{\pi}:\mathscr{P}(\mathbb{R})^{\mathcal{X}}\rightarrow\mathscr{P}(\mathbb{R})^{\mathcal{X}}$
is a $\gamma$-contraction with respect to $\bar{w}_{\infty}$. That is,

$\displaystyle\bar{w}_{\infty}(\mathcal{T}^{\pi}\eta,\mathcal{T}^{\pi}\eta^{\prime})\leq\gamma\bar{w}_{\infty}(\eta,\eta^{\prime})\,,$

for all $\eta,\eta^{\prime}\in\mathscr{P}(\mathbb{R})^{\mathcal{X}}$.

Next, we show that the projection operator $\Pi^{\lambda}$ cannot expand
distances as measured by $\bar{w}_{\infty}$, generalising the proof given by
Bellemare et al. (2023) in the case $\lambda\equiv 0$; the proof is given in
Appendix A.1.

######  Proposition 4.3.

The projection operator
$\Pi^{\lambda}:\mathscr{P}(\mathbb{R})^{\mathcal{X}}\rightarrow\mathscr{P}(\mathbb{R})^{\mathcal{X}}$
is a non-expansion with respect to $\bar{w}_{\infty}$. That is, for any
$\eta,\eta^{\prime}\in\mathscr{P}(\mathbb{R})^{\mathcal{X}}$, we have

$\displaystyle\bar{w}_{\infty}(\Pi^{\lambda}\eta,\Pi^{\lambda}\eta^{\prime})\leq\bar{w}_{\infty}(\eta,\eta^{\prime})\,.$

Finally, we put these two results together to obtain our desired conclusion.
In stating this result, it is useful here to introduce the notation

$\displaystyle\mathscr{F}_{\text{Q},m}=\Big{\\{}\sum_{i=1}^{m}\frac{1}{m}\delta_{z_{i}}:z_{i}\in\mathbb{R}\text{
for }i=1,\ldots,m\Big{\\}}\,,$

for the set of probability distributions representable with $m$ quantile
locations.

######  Proposition 4.4.

The projected operator
$\Pi^{\lambda}\mathcal{T}^{\pi}:\mathscr{F}_{\text{Q},m}^{\mathcal{X}}\rightarrow\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$
is a $\gamma$-contraction with respect to $\bar{w}_{\infty}$. Hence,
$\Pi^{\lambda}\mathcal{T}^{\pi}$ has a unique fixed point in
$\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$, which we denote
$\hat{\eta}^{\pi}_{\lambda}$. Further, given any initial
$\eta_{0}\in\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$, the sequence
$(\eta_{k})_{k=0}^{\infty}$ defined iteratively by
$\eta_{k+1}=\Pi^{\lambda}\mathcal{T}^{\pi}\eta_{k}$ for $k\geq 0$ satisfies
$\bar{w}_{\infty}(\eta_{k},\hat{\eta}^{\pi}_{\lambda})\leq\gamma^{k}\bar{w}_{\infty}(\eta_{0},\hat{\eta}^{\pi}_{\lambda})\rightarrow
0$.

###### Proof.

That
$\Pi^{\lambda}\mathcal{T}^{\pi}:\mathscr{F}_{\text{Q},m}^{\mathcal{X}}\rightarrow\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$
is a $\gamma$-contraction with respect to $\bar{w}_{\infty}$ follows directly
from Propositions 4.2 and 4.3:

$\displaystyle\bar{w}_{\infty}(\Pi^{\lambda}\mathcal{T}^{\pi}\eta,\Pi^{\lambda}\mathcal{T}^{\pi}\eta^{\prime})\leq\bar{w}_{\infty}(\mathcal{T}^{\pi}\eta,\mathcal{T}^{\pi}\eta^{\prime})\leq\gamma\bar{w}_{\infty}(\eta,\eta^{\prime})\,.$

Next, observe that $\bar{w}_{\infty}$ assigns finite distance to all pairs of
return-distribution functions in $\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$, and
further, this set is complete with respect to $\bar{w}_{\infty}$. Hence, we
may apply Banach’s fixed point theorem to obtain the existence of the unique
fixed point $\hat{\eta}^{\pi}_{\lambda}$ in
$\mathscr{F}_{\text{Q},m}^{\mathcal{X}}$. The final claim follows by
induction, and the contraction property established for
$\Pi^{\lambda}\mathcal{T}^{\pi}$. ∎

Note that the fixed point $\hat{\eta}^{\pi}_{\lambda}$ depends on $\lambda$,
and therefore implicitly on $m$. We also introduce the notation
$\hat{\theta}^{\pi}_{\lambda}\in\mathbb{R}^{\mathcal{X}\times[m]}$ for the
parameters of this collection of distributions, which is what the QDP
algorithm really operates over, so that we have

$\displaystyle\hat{\eta}^{\pi}_{\lambda}(x)=\sum_{i=1}^{m}\frac{1}{m}\delta_{\hat{\theta}^{\pi}_{\lambda}(x,i)}\,.$

Note that the convergence result of Proposition 4.4 also implies convergence
of the estimated quantile locations to $\hat{\theta}^{\pi}_{\lambda}$. In
Section 6, we will analyse the fixed point $\hat{\eta}^{\pi}_{\lambda}$, and
understand how closely it approximates the true return-distribution function
$\eta^{\pi}$. For now, having established convergence of QDP through
contraction mapping theory, we can return to QTD and demonstrate its own
convergence to the same fixed points.

##  5 Convergence of quantile temporal-difference learning

We now present the convergence analysis of QTD. We will consider a
_synchronous_ version of QTD, in which all states are updated using
independent transitions at each algorithm step, given by:

$\displaystyle\theta_{k+1}(x,i)=\theta_{k}(x,i)+\alpha_{k}\frac{1}{m}\sum_{j=1}^{m}(\tau_{i}-\mathbbm{1}\\{R_{k}(x)+\gamma\theta_{k}(X^{\prime}_{k}(x),j)<\theta_{k}(x,i)\\})\,,$
(13)

where given $x$ and $k$, we have $(R_{k}(x),X^{\prime}_{k}(x))\sim
P^{\pi}(\cdot|x)$, independently of the transitions used at all other
states/time steps, and $(\alpha_{k})_{k=0}^{\infty}$ is a sequence of step
sizes. The assumption of synchronous updates makes the analysis easier to
present, though it is possible to extend the analysis to the asynchronous
case, where a single state is updated at each algorithm time step (as would be
the case in fully online QTD, or an implementation using a replay buffer); see
Section 5.7. We now state the main convergence result of the paper.

######  Theorem 5.1.

Consider the sequence $(\theta_{k})_{k=0}^{\infty}$ defined by an initial
point $\theta_{0}\in\mathbb{R}^{\mathcal{X}\times[m]}$, the iterative update
in Equation (13), and non-negative step sizes satisfying the condition

$\displaystyle\sum_{t=0}^{\infty}\alpha_{k}=\infty\,,\quad\alpha_{k}=o(1/\log
k)\,.$ (14)

Then $(\theta_{k})_{k=0}^{\infty}$ converges almost surely to the set of fixed
points of the projected distributional Bellman operators
$\\{\Pi^{\lambda}\mathcal{T}^{\pi}:\lambda\in[0,1]^{\mathcal{X}\times[m]}\\}$;
that is,

$\displaystyle\inf_{\lambda\in[0,1]^{\mathcal{X}\times[m]}}\|\theta_{k}-\hat{\theta}^{\pi}_{\lambda}\|_{\infty}\rightarrow
0$

with probability 1.

Of particular note is the generality of this result. It does not require
finite-variance conditions on rewards (as is typically the case with
convergence results for classical TD); it holds for any collection of reward
distributions with the finite mean property set out at the beginning of the
paper. Some intuition as to why this is the case is that the finite-variance
conditions typically encountered are to ensure that the updates performed in
classical TD learning cannot grow in magnitude too rapidly. Since the updates
performed in QTD are bounded, this is not a concern, meaning that the proof
does not rely on such conditions. We note also that the step size conditions
are weaker than the typical Robbins-Monro conditions used in classical TD
analyses (see, e.g., Bertsekas and Tsitsiklis, 1996), which enforce square-
summability, also to avoid the possibility of divergence due to unbounded
noise in the classical TD learning.

The proof is based on the ODE method for stochastic approximation; in
particular we use the framework set out by Benaïm (1999) and Benaïm et al.
(2005). This involves interpreting the QTD update as a noisy Euler
discretisation of a differential equation (or more generally, a differential
inclusion). The broad steps are then to argue that the trajectories of the
differential equation/inclusion converge to some set of fixed points in a
suitable way (i.e., in a way that is robust to small perturbations), and that
the asymptotic behaviour of QTD, forming a noisy Euler discretisation, matches
the asymptotic behaviour of the true trajectories. This then allows us to
deduce that the QTD iterates converge to the same set of fixed points as the
true trajectories. We begin by elucidating the connection to differential
equations and differential inclusions.

###  5.1 The QTD differential equation

Taking the expectation over the random variables $R_{k}(x)$ and
$X^{\prime}_{k}(x)$ in Equation (13) conditional on the algorithm history up
to time $k$ yields an expected increment of

$\displaystyle\alpha_{k}\left(\tau_{i}-\mathbb{P}^{\pi}_{x}(R+\theta_{k}(X^{\prime},J)<\theta_{k}(x,i))\right)\,.$
(15)

We now briefly introduce an assumption on the MDP reward structure that
simplifies the analysis that follows. This assumption guarantees that the two
“difficult” cases of flat and vertical regions of CDFs (see Figure 1) do not
arise; we will lift this assumption later.

######  Assumption 5.2.

For each state $x\in\mathcal{X}$, the reward distribution at $x$ has a CDF
which is strictly increasing, and Lipschitz continuous.

As described in Section 4, the distribution of $R+\theta_{k}(X^{\prime},J)$
given the initial state $x$ is in fact equal to the application of the
distributional Bellman operator $\mathcal{T}^{\pi}$ applied to the return-
distribution function $\eta_{k}\in\mathscr{P}(\mathbb{R})^{\mathcal{X}}$ given
by

$\displaystyle\eta_{k}(x)=\frac{1}{m}\sum_{i=1}^{m}\delta_{\theta_{k}(x,i)}\,.$

Under Assumption 5.2, and in particular the assumption of continuous reward
CDFs, this yields a concise rewriting of the increment as

$\displaystyle\alpha_{k}\left(\tau_{i}-F_{(\mathcal{T}^{\pi}\theta_{k})(x)}(\theta_{k}(x,i))\right)\,.$

We may therefore intuitively interpret Equation (13) as a noisy discretisation
of the differential equation

$\displaystyle\partial_{t}\vartheta_{t}(x,i)=\tau_{i}-F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i))\,,$
(16)

which we refer to as the QTD differential equation (or QTD ODE). Note also
that Assumption 5.2 guarantees the global existence and uniqueness of
solutions to this differential equation.

######  Remark 5.3.

Calling back to Figure 3, the trajectories of the QTD ODE are obtained
precisely by integrating the vector fields that appear in these plots. In
contrast to the ODE that emerges when analysing classical TD learning (both in
tabular and linear function approximation settings) (Tsitsiklis and Van Roy,
1997), the right-hand side of Equation (16) is non-linear in the parameters
$\vartheta_{t}$, meaning that we are outside the domain of linear stochastic
approximation methods.

###  5.2 The QTD differential inclusion

In lifting Assumption 5.2, a few complications arise. Firstly, if
$F_{(\mathcal{T}^{\pi}\theta)(x)}$ is not continuous at $\theta(x,i)$, then
the right-hand side of the QTD ODE in Equation (16) is modified to

$\displaystyle\tau_{i}-\mathbb{P}_{Z\sim(\mathcal{T}^{\pi}\vartheta_{t})(x)}(Z<\vartheta_{t}(x,i))\,;$

the difference is the strict inequality. Now the right-hand side of the
differential equation itself is not continuous; in general, solutions may not
even exist for this differential equation. The situation is illustrated in the
bottom-left panel of Figure 3; the lines in this plot illustrate points of
discontinuity of the vector field to be integrated, and there are instances
where the vector field either side of such a line of discontinuity “pushes”
back into the discontinuity. In such cases, the differential equation has no
solution in the usual sense. This phenomenon is known as sliding, or sticking,
from cases when it arises in the modelling of physical systems with
potentially discontinuous forces (such as static friction models in
mechanics).

Filippov (1960) proposed a method to deal with such non-existence issues, by
relaxing the definition of the dynamics at points of discontinuity.
Technically, Filippov’s proposal is to allow the derivative to take on any
value in the convex hull of possible limiting values as we approach the point
of discontinuity. In our case, we consider redefining the dynamics at points
of discontinuity as follows:

$\displaystyle\partial_{t}\vartheta_{t}(x,i)\in[\tau_{i}-F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i)),\tau_{i}-F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i)-)]\,,$
(17)

where $F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i)-)$ denotes
$\lim_{s\uparrow\vartheta_{t}(x,i)}F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(s)$.
This refines the dynamics so that for each coordinate $(x,i)$, the derivative
may take on either the left or right limit around $\vartheta_{t}(x,i)$, or any
value in between; this is a looser relaxation than Filippov’s proposal, and is
easier to work with in our analysis.

Equation (17) is a _differential inclusion_ , as opposed to a differential
equation; the derivative is constrained to a set at each instant, rather than
constrained to a single value. We refer to Equation (17) specifically as the
QTD differential inclusion (or QTD DI). Note that if
$F_{(\mathcal{T}^{\pi}\theta)(x)}$ is continuous at $\theta(x,i)$, then the
right-hand side of Equation (17) reduces to the singleton
$\\{\tau_{i}-F_{(\mathcal{T}^{\pi}\theta)(x)}(\theta(x,i))\\}$, and we thus
obtain the ODE dynamics considered previously.

###  5.3 Solutions of differential inclusions

We briefly recall some key concepts regarding solutions of differential
inclusions; a full review of the theory of differential inclusions is beyond
the scope of this article, and we refer the reader to the standard references
by Aubin and Cellina (1984), Clarke et al. (1998), and Smirnov (2002).

######  Definition 5.4.

Let $H:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ be a set-valued map. The
path $(z_{t})_{t\geq 0}$ is a solution to the differential inclusion
$\partial_{t}z_{t}\in H(z_{t})$ if there exists an integrable function
$g:[0,\infty)\rightarrow\mathbb{R}^{n}$ such that

$\displaystyle z_{t}=\int_{0}^{t}g_{s}\mathrm{d}s$ (18)

for all $t\geq 0$, and $g_{t}\in H(z_{t})$ for almost all $t\geq 0$.

Note that Definition 5.4 does not require that $z_{t}$ is _differentiable_
with derivative $g_{t}$, but only the weaker integration condition in Equation
(18). We then have the following existence result (see e.g. Smirnov (2002) for
a proof).

######  Proposition 5.5.

Consider a set-valued map $H:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$,
and suppose that $H$ is a _Marchaud map_ : that is,

  * •

the set $\\{(z,h):z\in\mathbb{R}^{n},h\in H(z)\\}$ is closed.

  * •

For all $z\in\mathbb{R}^{n}$, $H(z)$ is non-empty, compact, and convex.

  * •

There exists a constant $C>0$ such that for all $z\in\mathbb{R}^{n}$,

$\displaystyle\max_{h\in H(z)}\|h\|\leq C(1+\|z\|)\,.$

Then the differential inclusion $\partial_{t}z_{t}\in H(z_{t})$ has a global
solution, for any initial condition.

It is readily verified that the QTD DI satisfies the requirements of this
result, and we are therefore guaranteed global solutions to this differential
inclusion, under any initial conditions.

###  5.4 Asymptotic behaviour of differential inclusion trajectories

Recall that our goal is to show that the trajectories of the QTD differential
inclusion must approach the fixed points of QDP. A key tool in doing so is the
notion of a Lyapunov function; the following definition is based on Benaïm et
al. (2005).

######  Definition 5.6.

Consider a Marchaud map $H:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$, and
a subset $\Lambda\subseteq\mathbb{R}^{n}$. A continuous function
$L:\mathbb{R}^{n}\rightarrow[0,\infty)$ is said to be a _Lyapunov function_
for the differential inclusion $\partial_{t}z_{t}\in H(z_{t})$ and subset
$\Lambda$ if for any solution $(z_{t})_{t\geq 0}$ of the differential
inclusion and $0\leq s<t$, we have $L(z_{t})<L(z_{s})$ for all
$z_{s}\not\in\Lambda$ and $L(z)=0$ for all $z\in\Lambda$.

Intuitively, $L$ is a Lyapunov function if it decreases along trajectories of
the differential inclusion, and is minimal precisely on $\Lambda$. Lyapunov
functions are a central tool in dynamical systems for demonstrating
convergence, and in the sections that follow, we will consider the QTD
differential inclusion, and take $\Lambda$ to be the set of fixed points of
the family of QDP algorithms.

###  5.5 QTD as a stochastic approximation to the QTD differential inclusion

We can now give the proof of our core result, Theorem 5.1. The abstract
stochastic approximation result at the heart of the convergence proof of QTD
is presented below. It is a special case of the general framework described by
Benaïm et al. (2005), the proof of which is given in Appendix A.2.

######  Theorem 5.7.

Consider a Marchaud map $H:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$, and
the corresponding differential inclusion $\partial_{t}z_{t}\in H(z_{t})$.
Suppose there exists a Lyapunov function $L$ for this differential inclusion
and a subset $\Lambda\subseteq\mathbb{R}^{n}$. Suppose also that we have a
sequence $(\theta_{k})_{k\geq 0}$ satisfying

$\displaystyle\theta_{k+1}=\theta_{k}+\alpha_{k}(g_{k}+w_{k})\,,$

where:

  * •

$(\alpha_{k})_{k=0}^{\infty}$ satisfy the conditions
$\sum_{k=0}^{\infty}\alpha_{k}=\infty$, $\alpha_{k}=o(1/\log(k))$;

  * •

$g_{k}\in H(\theta_{k})$ for all $k\geq 0$;

  * •

$(w_{k})_{k=0}^{\infty}$ is a bounded martingale difference sequence with
respect to the natural filtration generated by $(\theta_{k})_{k=0}^{\infty}$;
that is, there is an absolute constant $C$ such that $\|w_{k}\|_{\infty}<C$
almost surely, and $\mathbb{E}[w_{k}|\theta_{0},\ldots,\theta_{k}]=0$.

If further $(\theta_{k})_{k=0}^{\infty}$ is bounded almost surely (that is,
$\sup_{k\geq 0}\|\theta_{k}\|_{\infty}<\infty$ almost surely), then
$\theta_{k}\rightarrow\Lambda$ almost surely.

The intuition behind the conditions of the theorem are as follows. The
Marchaud map condition ensures the differential inclusion of interest has
global solutions. The existence of the Lyapunov function guarantees that
trajectories of the differential inclusion converge in a suitably stable sense
to $\Lambda$. The step size conditions, martingale difference condition, and
boundedness conditions mean that the iterates $(\theta_{k})_{k=0}^{\infty}$
will closely track the differential inclusion trajectories, and hence exhibit
the same asymptotic behaviour. We can now give the proof of Theorem 5.1, first
requiring the following proposition, which is proven in Appendix A.3.

######  Proposition 5.8.

Under the conditions of Theorem 5.1, the iterates
$(\theta_{k})_{k=0}^{\infty}$ are bounded almost surely.

###### Proof of Theorem 5.1.

We see that for the QTD sequence $(\theta_{k})_{k=0}^{\infty}$ and the QTD DI
and QDP invariant set
$\Lambda=\\{\hat{\theta}^{\pi}_{\lambda}:\lambda\in[0,1]^{\mathcal{X}\times[m]}\\}$,
the conditions of Theorem 5.7 are satisfied, except perhaps for the
boundedness of $(\theta_{k})_{k=0}^{\infty}$, and the existence of the
Lyapunov function. The fact that the sequence $(\theta_{k})_{k=0}^{\infty}$ is
bounded almost surely is Proposition 5.8; its proof is somewhat technical, and
given in the appendix. The construction of a valid Lyapunov function is given
in Proposition 5.11 below, which completes the proof. ∎

######  Remark 5.9.

What makes the relaxation to the differential inclusion work in this analysis?
We have already seen that some kind of relaxation of the dynamics is required
in order to define a valid continuous-time dynamical system; the original ODE
may not have solutions in general. If we relax the dynamics too much (an
extreme example would be the differential inclusion
$\vartheta_{t}(x,i)\in\mathbb{R}$), what goes wrong? The answer is that there
are too many resulting solutions, which do not exhibit the desired asymptotic
behaviour. Thus, the differential inclusion in Equation (17) is in some sense
just the right level of relaxation of the differential equation we started
with, since trajectories of the QTD DI are still guaranteed to converge to the
QDP fixed points.

###  5.6 A Lyapunov function for the QDP fixed points

In this section, we prove the existence of a Lyapunov function required in
order to use Theorem 5.7 to prove Theorem 5.1. We treat the case when
Assumption 5.2 holds separately as the proof is instructive, and considerably
simpler than the general case. Under this assumption, note that all
projections $\Pi^{\lambda}$ behave identically on the image of
$\mathcal{T}^{\pi}$, since all resulting CDFs are strictly increasing. We
therefore introduce the notation $\Pi$ to refer to any such projection in this
case, and the notation $\hat{\theta}^{\pi}_{m}$ to refer to the unique fixed
point of $\Pi\mathcal{T}^{\pi}$.

######  Proposition 5.10.

Consider the ODE in Equation (16), and suppose Assumption 5.2 holds. A
Lyapunov function for the equilibrium point $\hat{\theta}^{\pi}_{m}$ is given
by

$\displaystyle
L(\theta)=\max_{x\in\mathcal{X}}\max_{i=1,\ldots,m}|\theta(x,i)-\hat{\theta}^{\pi}_{m}(x,i)|\,.$

###### Proof.

We immediately observe that $L$ is continuous, non-negative, and takes on the
value $0$ only at $\hat{\theta}^{\pi}_{m}$. To show that $L(\vartheta_{t})$ is
decreasing, where $(\vartheta_{t})_{t\geq 0}$ is an ODE trajectory, suppose
$(x,i)$ is a state-index pair attaining the maximum in $L(\vartheta_{t})$. It
is sufficient to show that $\vartheta_{t}(x,i)$ is moving towards
$\hat{\theta}^{\pi}_{m}(x,i)$, or expressed mathematically,

$\displaystyle\partial_{t}\vartheta_{t}(x,i)\overset{\text{S}}{=}\hat{\theta}^{\pi}_{m}(x,i)-\vartheta_{t}(x,i)\,,$

where we use $a\overset{\text{S}}{=}b$ as shorthand for _equality of signs_
$\text{sign}(a)=\text{sign}(b)$, where

$\displaystyle\text{sign}(z)=\begin{cases}1&\text{ if }z>0\,,\\\ 0&\text{ if
}z=0\,,\\\ -1&\text{ if }z<0\,.\\\ \end{cases}$

Now note that

$\displaystyle\partial_{t}\vartheta_{t}(x,i)$
$\displaystyle=\tau_{i}-F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i))$
$\displaystyle\overset{\text{S}}{=}F^{-1}_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\tau_{i})-F^{-1}_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}(\vartheta_{t}(x,i)))$
$\displaystyle=(\Pi\mathcal{T}^{\pi}\vartheta_{t})(x,i)-\vartheta_{t}(x,i)\,,$

where the sign equality follows from Assumption 5.2; since all reward CDFs are
strictly increasing, so too is $F_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}$, and
so $F^{-1}_{(\mathcal{T}^{\pi}\vartheta_{t})(x)}$ is strictly monotonic.
Additionally, from the contractivity of $\Pi\mathcal{T}^{\pi}$ with respect to
$\overline{w}_{\infty}$ (see Proposition 4.4), we have

$\displaystyle|(\Pi\mathcal{T}^{\pi}\vartheta_{t})(x,i)-\hat{\theta}^{\pi}_{m}(x,i)|$
$\displaystyle\leq\overline{w}_{\infty}(\Pi\mathcal{T}^{\pi}\vartheta_{t},\Pi\mathcal{T}^{\pi}\hat{\theta}^{\pi}_{m})$
$\displaystyle\leq\gamma\max_{x\in\mathcal{X}}\max_{i=1,\ldots,m}|\vartheta_{t}(x,i)-\hat{\theta}^{\pi}_{m}(x,i)|$
$\displaystyle=\gamma|\vartheta_{t}(x,i)-\hat{\theta}^{\pi}_{m}(x,i)|\,.$

From this, we deduce

$\displaystyle(\Pi\mathcal{T}^{\pi}\vartheta_{t})(x,i)-\vartheta_{t}(x,i)\overset{\text{S}}{=}\hat{\theta}^{\pi}_{m}(x,i)-\vartheta_{t}(x,i)\,,$

and so

$\displaystyle\partial_{t}\vartheta_{t}(x,i)\overset{\text{S}}{=}\hat{\theta}^{\pi}_{m}(x,i)-\vartheta_{t}(x,i)\,.$

We therefore have that $L(\vartheta_{t})$ is decreasing at $t$, strictly so if
$\vartheta_{t}\not=\hat{\theta}^{\pi}_{m}$, as required to establish the
result. ∎

The proof of Proposition 5.10 also sheds further light on the mechanisms
underlying the QTD algorithm. A key step in the argument is to show that for
the state-index pairs $(x,i)$ such that $\vartheta_{t}(x,i)$ is maximally
distant from the fixed point $\theta^{\pi}_{m}(x,i)$, the expected update
under QTD moves this coordinate of the estimate in the same direction as
gradient descent on a squared loss from the fixed point. However, the fact
that it is only the _sign_ of the update that has this property, and not its
magnitude, means that the empirical rate of convergence and stability of QTD
can be expected to be somewhat different from methods based on an $L^{2}$
loss, such as classical TD.

We now state the Lyapunov result in the general case; the proof is somewhat
more involved, and is given in Appendix A.4.

######  Proposition 5.11.

The function

$\displaystyle
L(\theta)=\min_{\lambda\in[0,1]^{\mathcal{X}\times[m]}}\max_{(x,i)}|\theta(x,i)-\hat{\theta}^{\pi}_{\lambda}(x,i)|$
(19)

is a Lyapunov function for the differential inclusion in Equation (17) and the
set of fixed points
$\\{\hat{\theta}^{\pi}_{\lambda}:\lambda\in[0,1]^{\mathcal{X}\times[m]}\\}$.

###  5.7 Extension to asynchronous QTD

Our convergence results have focused on the synchronous case of QTD. However,
in practice, it is often of interest to implement _asynchronous_ versions of
TD algorithms, in which only a single state is updated at a time. More
formally, an asynchronous version of QTD computes the sequence
$(\theta_{k})_{k\geq 0}$ defined by an initial estimate
$\theta\in\mathbb{R}^{\mathcal{X}\times[m]}$, a sequence of transitions
$(X_{k},R_{k},X^{\prime}_{k})_{k\geq 0}$, and the update rule

$\displaystyle\theta_{k+1}(x,i)=\theta_{k}(x,i)+\beta_{x,k}\frac{1}{m}\sum_{j=1}^{m}(\tau_{i}-\mathbbm{1}\\{R_{k}+\gamma\theta_{k}(X^{\prime}_{k},j)<\theta_{k}(x,i)\\})$

for $x=X_{k}$, and $\theta_{k+1}(x,i)=\theta_{k}(x,i)$ otherwise. Here, the
step size $\beta_{x,k}$ depends on both $x$ and $k$, and is typically selected
so that each state individually makes use of a fixed step size sequence
$(\alpha_{k})_{k=0}^{\infty}$, by taking
$\beta_{x,k}=\alpha_{\sum_{l=0}^{k}\mathbbm{1}\\{X_{l}=x\\}}$. This models the
online situation where a stream of experience $(X_{k},R_{k})_{k\geq 0}$ is
generated by interacting with the environment using the policy $\pi$, and
updates are performed setting $X^{\prime}_{k}=X_{k+1}$, and also the setting
in which the tuples $(X_{k},R_{k},X^{\prime}_{k})_{k\geq 0}$ are sampled
i.i.d. from a replay buffer, among others.

Convergence of QTD in such asynchronous settings can also be proven; Perkins
and Leslie (2013) extend the analysis of Benaïm (1999) and Benaïm et al.
(2005), incorporating the approach of Borkar (1998), to obtain convergence
guarantees for asynchronous stochastic approximation algorithms approximating
differential inclusions. In the interest of space, we do not provide the full
details of the proof here, but instead sketch the key differences that arise
in the analysis in Appendix C.

##  6 Analysis of the QTD limit points

In general, the limiting points $\hat{\eta}^{\pi}_{\lambda}$ for QTD/QDP will
not be the same as the true return-distribution function $\eta^{\pi}$. On the
one hand, this is clear; each return-distribution function
$\hat{\eta}^{\pi}_{\lambda}$ is in the image of the projection
$\Pi^{\lambda}$, so each constituent probability distribution must be of the
form $\tfrac{1}{m}\sum_{i=1}^{m}\delta_{z_{i}}$, whereas the true return
distributions need not take on this form. In addition, the magnitude of this
approximation error is not immediately clear. Each application of the
projection $\Pi^{\lambda}$ in the dynamic programming process causes some loss
of information, and the quality of the fixed point
$\hat{\eta}^{\pi}_{\lambda}$ is affected by the build up of these
approximations over time.

Measuring approximation error in $\bar{w}_{\infty}$ typically turns out to be
uninformative, as $\bar{w}_{\infty}$ is a particularly strict notion of
distance between probability distributions, as discussed in the context of
distributional RL by Rowland et al. (2019) and Bellemare et al. (2023). In
particular, fixed points $\hat{\eta}^{\pi}_{\lambda}$ that intuitively provide
a good approximation to $\eta^{\pi}$ may have high
$\bar{w}_{\infty}$-distance, and the $\bar{w}_{\infty}$-distance generally
does not decrease with $m$ (Bellemare et al., 2023). Instead, we use the
Wasserstein-1 metric, and its extension to return-distribution functions,
defined by

$\displaystyle
w_{1}(\nu,\nu^{\prime})=\int_{0}^{1}|F^{-1}_{\nu}(t)-F^{-1}_{\nu^{\prime}}(t)|\;\mathrm{d}t\,,\quad\bar{w}_{1}(\eta,\eta^{\prime})=\max_{x\in\mathcal{X}}w_{1}(\eta(x),\eta^{\prime}(x))\,,$

for all $\nu,\nu^{\prime}\in\mathscr{P}(\mathbb{R})$, and
$\eta,\eta^{\prime}\in\mathscr{P}(\mathbb{R})^{\mathcal{X}}$. The following
result improves on the analysis given by Bellemare et al. (2023) for the case
of $\lambda\equiv 0$, establishing an upper bound on the $\overline{w}_{1}$
distance between $\hat{\eta}^{\pi}_{\lambda}$ and $\eta^{\pi}$ for any
$\lambda$, essentially by showing that the errors accumulated in dynamic
programming can be made arbitrarily small by increasing $m$, which controls
the richness of the distribution representation.

######  Proposition 6.1.

For any $\lambda\in[0,1]^{\mathcal{X}\times[m]}$, if all reward distributions
are supported on $[R_{\textsc{min}},R_{\textsc{max}}]$, then we have

$\displaystyle\bar{w}_{1}(\hat{\eta}^{\pi}_{\lambda},\eta^{\pi})\leq\frac{V_{\textsc{max}}-V_{\textsc{min}}}{2m(1-\gamma)}\,,$

where $V_{\textsc{max}}=R_{\textsc{max}}/(1-\gamma)$, and similarly
$V_{\textsc{min}}=R_{\textsc{min}}/(1-\gamma)$.

######  Remark 6.2.

This bound also provides motivation for the specific values of
$(\tau_{i})_{i=1}^{m}$ that QTD uses. A similar convergence analysis and
fixed-point analysis can be straightforwardly carried out for a version of the
QTD algorithm with other values for $(\tau_{i})_{i=1}^{m}$; by tracing through
the proof of Proposition 6.1, it can be seen that the bound is proportional to
$\max(\tau_{1},\max((\tau_{i+1}-\tau_{i})/2:i=2,\ldots,m-1),1-\tau_{m})$,
which is minimised precisely by the choice of $(\tau_{i})_{i=1}^{m}$ used by
QTD.

###  6.1 Instance-dependent bounds

The result above implicitly assumes the worst-case projection error is
incurred at all states with each application of the Bellman operator. In
environments where this is not the case, the fixed point can be shown to be of
considerably better quality. We describe an example of an instance-dependent
quality bound here.

######  Proposition 6.3.

Consider an MDP such that for any trajectory, after $k$ time steps all
encountered transition distributions and reward distributions are Dirac
deltas. If all reward distributions in the MDP are supported on
$[R_{\textsc{min}},R_{\textsc{max}}]$, then for any
$\lambda\in[0,1]^{\mathcal{X}\times[m]}$, we have

$\displaystyle\bar{w}_{1}(\hat{\eta}^{\pi}_{\lambda},\eta^{\pi})\leq\frac{(V_{\textsc{max}}-V_{\textsc{min}})(1-\gamma^{k})}{2m(1-\gamma)}\,.$

######  Remark 6.4.

One particular upshot of this bound for practitioners is that for agents in
near-deterministic environments using near-deterministic policies, it may be
possible to use $m=o((1-\gamma)^{-1})$ quantiles and still obtain accurate
approximations to the return-distribution function via QTD and/or QDP. It is
interesting to contrast this result for quantile-based distributional
reinforcement learning against the case when using categorical distribution
representations (Bellemare et al., 2017; Rowland et al., 2018; Bellemare et
al., 2023). In this latter case, fixed point error continues to be accumulated
even when the environment has solely deterministic transitions and rewards,
due to the well-documented phenomenon of the approximate distribution
‘spreading its mass out’ under the Cramér projection (Bellemare et al., 2017;
Rowland et al., 2018; Bellemare et al., 2023). Our observation here leads to
immediate practical advice for practitioners (in environments with mostly
deterministic transitions, a quantile representation may be preferred to a
categorical representation, leading to less approximation error), and raises a
general question that warrants further study: how can we use prior knowledge
about the structure of the environment to select a good distribution
representation?

We conclude this section by noting that many variants of Proposition 6.3 are
possible; one can for example modify the assumption that rewards are
deterministic to an assumption that rewards distributions are supported on a
‘small’ interval, and still obtain a fixed-point bound that improves over the
instance-independent bound of Proposition 6.1. There are a wide variety of
such modifications that could be imagined, and we believe this to be an
interesting direction for future research and applications.

###  6.2 Qualitative analysis of QDP fixed points

The analysis in the previous section establishes quantitative upper bounds on
the quality of the fixed point learnt by QTD, and guarantees that with enough
atoms an arbitrarily accurate approximation of the return-distribution
function (as measured by $w_{1}$) can be learnt. We now take a closer look at
the way in which approximation errors may manifest in QTD and QDP.

######  Example 6.5.

Consider the two-state Markov decision process (with a single action) whose
transition probabilities are specified by the left-hand side of Figure 5, such
that a deterministic reward of $2$ is obtained in state $x_{1}$, and $-1$ in
state $x_{2}$; further, let us take a discount factor $\gamma=0.9$. The centre
panel of this figure shows various estimates of the CDF for the return
distribution at state $x_{1}$. The ground truth estimate in black is obtained
from Monte Carlo sampling. The CDFs in purple, blue, green, and orange are the
points of convergence for QDP with $m=2,5,10,100$, respectively. For $m=100$,
a very close fit to the true return distribution is obtained. However, for
small $m$ in particular, the distribution is heavily skewed to the right. In
the case of $m=2$, half of the probability mass is placed on the greatest
possible return in this MDP – namely 20 – even though with probability 1 the
true return is less than this value. What is the cause of this behaviour from
QDP? This question is answered by investigating the dynamic programming update
itself in more detail.

In this MDP, the result of the QDP operator applied to the fixed point
$\theta$ is to update each particle location with a ‘backed-up’ particle
appearing in the distributions $\mathcal{T}^{\pi}\theta$. When such settings
arise, tracking which backed-up particles are allocated to which other
particles helps us to understand the behaviour of QDP, and the nature of the
approximation incurred. We also gain intuition about the situation, since the
QDP operator is behaving like a an affine policy evaluation operator on
$\mathcal{X}\times[m]$ _locally_ around the fixed point. We can visualise
which particles are assigned to one another by a QDP operator application
through _local quantile back-up diagrams_ ; the right-hand side of Figure 5
shows the local quantile back-up diagram for particular MDP. We observe that
$\theta(x_{1},2)$ backs up from itself, and hence learns a value that
corresponds to observing a self-transition at every state, with a reward of 2;
under the discount factor of $0.9$, this corresponds to a return of 20. This
is the source of the drastic over-estimation of returns in the approximation
obtained with $m=2$, and the fact that all other state-quantile pairs
implicitly bootstrap from this estimate leads to the over-estimation leaking
out into all quantiles estimated in this case. As $m$ increases, we observe
from the CDF plot that there is always one particle that learns this maximal
return of 20, but that this has less effect on the other quantiles; indeed in
the orange curve, we obtain a very good approximation (in $w_{1}$) to the true
return distribution despite this particle with a maximal value of 20 remaining
present. We can interpret the increase in $m$ as preventing pathological self-
loops/small cycles in the quantile backup diagram from “leaking out” and
degrading the quality of other quantile estimates; this provides a
complementary perspective on the approximation artefacts that occur in QDP/QTD
fixed points to the quantitative upper bounds in the previous section.

We expect the local quantile back-up diagram introduced in Example 6.5 to be a
useful tool for developing intuition, as well as further analysis, of QDP and
QTD. As described in the example itself, being able to define the local back-
up diagram depends on the structure of the MDP being such that the QDP
operator obtains each new coordinate value from a single backed-up particle
location. It is an interesting question as to how the definition of such local
back-up diagrams could be generalised to apply in situations where this does
not hold, such as with absolutely continuous reward distributions.

Figure 5:  Left: An example MDP. Centre: The fixed point return distribution
estimates for state $x_{1}$ obtained by QDP for $m=2,5,20,100$ (solid purple,
dotted blue, dashed green, and dash-dotted orange, respectively) compared to
ground truth in solid black. Right: The corresponding local quantile backup
diagram at the fixed point for $m=2$, illustrating potential approximation
artefacts in QDP fixed points.

##  7 Related work

Stochastic approximation theory with differential inclusions. The ODE method
was introduced by Ljung (1977) as a means of analysing stochastic
approximation algorithms, and was subsequently extended and refined by Kushner
and Clark (1978); standard references on the subject include Kushner and Yin
(2003); Borkar (2008); Benveniste et al. (2012); see also Meyn (2022) for an
overview in the context of reinforcement learning. The framework we follow in
this paper is set out by Benaïm (1999), and was extended by Benaïm et al.
(2005) to allow for differential inclusions. Perkins and Leslie (2013) later
extended this analysis further to allow for asynchronous algorithms, building
on the approach introduced by Borkar (1998), and extended, with particular
application to reinforcement learning, by Borkar and Meyn (2000).

Differential inclusion theory. Differential inclusions have found application
across a wide variety of fields, including control theory (Wazewski, 1961),
economics (Aubin, 1991) differential game theory (Krasovskii and Subbotin,
1988), and mechanics (Monteiro Marques, 2013). The approach to modelling
differential equations with discontinuous right-hand sides via differential
inclusions was introduced by Filippov (1960). Standard references on the
theory of differential inclusions include Aubin and Cellina (1984); Clarke et
al. (1998); Smirnov (2002); see also Bernardo et al. (2008) on the related
field of piecewise-smooth dynamical systems. Joseph and Bhatnagar (2019) also
use tools combining stochastic approximation and differential inclusions from
Benaïm et al. (2005) to analyse (sub-)gradient descent as a means of
estimating quantiles of fixed distributions. Within reinforcement learning and
related fields more specifically, differential inclusions have played a key
role in the analysis of game-theoretic algorithms based on fictitious play
(Brown, 1951; Robinson, 1951); see Benaïm et al. (2006); Leslie and Collins
(2006); Benaïm and Faure (2013) for examples. More recently, Gopalan and
Thoppe (2022) used differential inclusions to analyse TD algorithms for
control with linear function approximation.

Quantile regression. Quantile regression as a methodology for statistical
inference was introduced by Koenker and Bassett (1978). Koenker (2005) and
Koenker et al. (2017) provide detailed surveys of the field. Quantile
temporal-difference learning may be viewed as fusing quantile regression with
the bootstrapping approach (_learning a guess from a guess_ , as Sutton and
Barto (2018) express it) that is core to much of the reinforcement learning
methodology.

Quantiles in reinforcement learning. The approach to distributional
reinforcement learning based on quantiles was introduced by Dabney et al.
(2018b). A variety of modifications and extensions were then considered in the
deep reinforcement learning setting (Dabney et al., 2018a; Yang et al., 2019;
Zhou et al., 2020; Luo et al., 2021), as well as further developments on the
theoretical side (Lhéritier and Bondoux, 2022). A summary of the appraoch is
presented by Bellemare et al. (2023). Gilbert and Weng (2016) study the
problem of optimising quantile criteria in end-state MDPs. Li et al. (2022)
consider the risk-sensitive control problem of optimising particular quantiles
of the return distribution, and derive a dynamic programming algorithm that
maintains a value function over state and the “target quantile-to-go”.

##  8 Conclusion

We have provided the first convergence analysis for QTD, a popular and
effective distributional reinforcement learning algorithm. In contrast to the
analysis of many classical temporal-difference learning algorithms, this has
required the use of tools from the field of differential inclusions and
branches of stochastic approximation theory that deal with the associated
dynamical systems. Due to the structure of the QTD algorithm, such as its
bounded-magnitude updates, these convergence guarantees hold under weaker
conditions than are generally used in the analysis of TD algorithms. These
results establish the soundness of QTD, representing an important step towards
understanding its efficacy and practical successes, and the theoretical tools
used here are likely to be useful in further analyses of (distributional)
reinforcement learning algorithms. We believe further research into the
theory, practice and applications of QTD, at a variety of scales, are
important directions for foundational reinforcement learning research. Some
natural directions following on from this work are to understand the
interaction of QTD with function approximation, and adaptive optimisation.

Acknowledgements

We thank Tor Lattimore for detailed comments on an earlier draft, and David
Abel, Bernardo Avila Pires, Diana Borsa, Yash Chandak, Daniel Guo, Clare Lyle,
and Shantanu Thakoor for helpful discussions. The simulations in this paper
were generated using the Python 3 language, and made use of the NumPy (Harris
et al., 2020), SciPy (Virtanen et al., 2020), and Matplotlib (Hunter, 2007)
libraries.

##  Appendix A Proofs

In this section, we provide proofs for results which are not proven in the
main text.

###  A.1 Proof of Proposition 4.3

Let $\eta,\eta^{\prime}\in\mathscr{P}(\mathbb{R})$. We have

$\displaystyle\Big{|}\left((1-\lambda(x,i))F^{-1}_{\eta(x)}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\eta(x)}(\tau_{i})\right)\\!-\\!\left((1-\lambda(x,i))F^{-1}_{\eta^{\prime}(x)}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\eta^{\prime}(x)}(\tau_{i})\right)\Big{|}$
$\displaystyle\qquad\qquad\leq(1-\lambda(x,i))|F^{-1}_{\eta(x)}(\tau_{i})-F^{-1}_{\eta^{\prime}(x)}(\tau_{i})|+\lambda(x,i)|\bar{F}^{-1}_{\eta(x)}(\tau_{i})-\bar{F}^{-1}_{\eta^{\prime}(x)}(\tau_{i})|\,.$

Clearly

$\displaystyle|F^{-1}_{\eta(x)}(\tau_{i})-F^{-1}_{\eta^{\prime}(x)}(\tau_{i})|\leq\bar{w}_{\infty}(\eta,\eta^{\prime})\,.$

Additionally, we have

$\displaystyle|\bar{F}^{-1}_{\eta(x)}(\tau_{i})-\bar{F}^{-1}_{\eta^{\prime}(x)}(\tau_{i})|$
$\displaystyle=|\lim_{s\downarrow\tau_{i}}F^{-1}_{\eta(x)}(s)-\lim_{s\downarrow\tau_{i}}F^{-1}_{\eta^{\prime}(x)}(s)|$
$\displaystyle=\lim_{s\downarrow\tau_{i}}|F^{-1}_{\eta(x)}(s)-F^{-1}_{\eta^{\prime}(x)}(s)|$
$\displaystyle\leq\bar{w}_{\infty}(\eta,\eta^{\prime})\,.$

Putting this together, we obtain

$\displaystyle\Big{|}\left((1-\lambda(x,i))F^{-1}_{\eta(x)}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\eta(x)}(\tau_{i})\right)\\!-\\!\left((1-\lambda(x,i))F^{-1}_{\eta^{\prime}(x)}(\tau_{i})+\lambda(x,i)\bar{F}^{-1}_{\eta^{\prime}(x)}(\tau_{i})\right)\Big{|}$
$\displaystyle\qquad\leq\bar{w}_{\infty}(\eta,\eta^{\prime})\,,$

as required.

###  A.2 Proof of Theorem 5.7

Theorem 5.7 is essentially a special case of the general results presented in
Benaïm et al. (2005), in the form needed for the proof of convergence of QTD.
To explain how to obtain Theorem 5.7 from the results of Benaïm et al. (2005),
first, we associate a continuous-time path $(\bar{\theta}(t))_{t\geq 0}$ with
the iterates $(\theta_{k})_{k=0}^{\infty}$ by linear interpolation, in
particular defining $\bar{\theta}(\sum_{k=0}^{s}\alpha_{k})=\theta_{s}$, and
linearly interpolating in between. The continuous-time path
$(\bar{\theta}(t))_{t\geq 0}$ satisfies the definition of a _perturbed
solution_ of the Marchaud differential inclusion with probability 1, as
defined by Definition II of Benaïm et al. (2005), since: (i) $\bar{\theta}$ is
piecewise linear, hence absolutely continuous; (ii) the difference
$\|\theta_{k+1}-\theta_{k}\|_{\infty}$ is $O(\alpha_{k})$, due to the growth
condition on $H$ and since $\bar{\theta}$ is bounded by assumption; and (iii)
the lim-sup condition holds with probability 1 thanks to the boundedness of
the martingale difference sequence $(w_{k})_{k=0}^{\infty}$ and Proposition
1.4 of Benaïm et al. (2005); see also Theorem 5.3.3 of Kushner and Yin (2003).

Next, since we assume $\bar{\theta}$ is bounded, Theorem 4.2 of Benaïm et al.
(2005) applies so that we deduce that it is an asymptotic pseudotrajectory of
the differential inclusion (w.p.1). We then have that
$(\bar{\theta}(t))_{t\geq 0}$ is a bounded asymptotic pseudotrajectory
(w.p.1), so Theorem 4.3 of Benaïm et al. (2005) applies, and we deduce that
the set of limit points of $(\bar{\theta}(t))_{t\geq 0}$ is internally chain
transitive (w.p.1). But now by Proposition 3.27 of Benaïm et al. (2005)
applied to the Lyapunov function $L$ and the set $\Lambda$, all internally
chain transitive sets are contained within $\Lambda$. Since
$(\bar{\theta}(t))_{t\geq 0}$ is bounded, we deduce that it converges to
$\Lambda$ (w.p.1). It therefore follows that the discrete sequence
$(\theta_{k})_{k=0}^{\infty}$ converges to $\Lambda$ with probability 1, as
required.

###  A.3 Proof of Proposition 5.8

Roughly, the intuition of the proof is that the structure of the QTD
differential inclusion means that when $\|\theta_{k}\|_{\infty}$ is
sufficiently large, the coordinates of $\theta_{k}$ furthest from the origin
are moved back towards the origin by the differential inclusion. We then argue
that the martingale noise cannot cause divergence, which completes the
argument.

Differential inclusion update direction. To begin with the analysis of the
differential inclusion, fix $\delta>0$ such that $1-\delta>\gamma$, and let
$M>0$ be such that for all $(x,a)\in\mathcal{X}\times\mathcal{A}$, we have

$\displaystyle F_{P_{\mathcal{R}}(x,a)}((1-\delta-\gamma)M)>1-1/(4m)\,,\quad
F_{P_{\mathcal{R}}(x,a)}(-(1-\delta-\gamma)M)<1/(4m)\,.$

We then introduce the events

$\displaystyle I^{+}_{k}(x,i)$ $\displaystyle=\\{\|\theta_{k}\|_{\infty}>M\,,\
\ \theta_{k}(x,i)>(1-\delta)\|\theta_{k}\|_{\infty}\\}\,,$ $\displaystyle
I^{-}_{k}(x,i)$ $\displaystyle=\\{\|\theta_{k}\|_{\infty}>M\,,\ \
\theta_{k}(x,i)<-(1-\delta)\|\theta_{k}\|_{\infty}\\}\,.$

which, roughly speaking, hold when $\theta_{k}$ has at least one large
coordinate (in absolute value), and $\theta_{k}(x,i)$ is a positive
(respectively, negative) coordinate close to the maximum value. When
$I^{+}_{k}(x,i)$ holds, we have

$\displaystyle\phantom{=}\tau_{i}-F_{(\mathcal{T}^{\pi}\theta_{k})(x)}(\theta_{k}(x,i)-)$
(20)
$\displaystyle=\tau_{i}-\sum_{x^{\prime}\in\mathcal{X}}P(x^{\prime}|x,a)\pi(a|x)\frac{1}{m}\sum_{j=1}^{m}\lim_{s\uparrow
0}F_{\mathcal{R}(x,a)}(\theta_{k}(x,i)-\gamma\theta_{k}(x^{\prime},j)+s)$
$\displaystyle\overset{(a)}{\leq}\tau_{i}-\sum_{x^{\prime}\in\mathcal{X}}P(x^{\prime}|x,a)\pi(a|x)\frac{1}{m}\sum_{j=1}^{m}(1-1/(4m))$
$\displaystyle\leq(1-1/(2m))-(1-(1/(4m))$ $\displaystyle\leq-1/(4m)\,,$ (21)

and hence the differential inclusion moves $\theta_{k}(x,i)$ towards the
origin. Inequality (a) follows since on $I^{+}_{k}(x,i)$, we have

$\displaystyle\theta_{k}(x,i)-\gamma\theta_{k}(x^{\prime},j)\geq(1-\delta)\|\theta_{k}\|_{\infty}-\gamma\|\theta_{k}\|_{\infty}\geq(1-\delta-\gamma)M\,.$

Analogously, we conclude that on $I^{-}_{k}(x,i)$, we have

$\displaystyle\tau_{i}-F_{(\mathcal{T}^{\pi}\theta_{k})(x)}(\theta_{k}(x,i))\geq
1/(4m)\,,$

and so the differential inclusion moves $\theta_{k}(x,i)$ towards the origin
in this case too.

Chaining updates and reasoning about noise. To describe the relationship
between successive iterates in the sequence $(\theta_{k})_{k\geq 0}$, we
introduce the notation
$\theta_{k+1}=\theta_{k}+\alpha_{k}g_{k}+\alpha_{k}w_{k}$, where $w_{j}$ is
martingale difference noise, and hence $g_{j}$ is an expected update
direction, from the right-hand side of the QTD differential inclusion. By
boundedness of the update noise and the step size assumptions, we have from
Proposition 1.4 of Benaïm et al. (2005) (see also Theorem 5.3.3 of Kushner and
Yin (2003)) that

$\displaystyle\lim_{k}\sup\Big{\\{}\big{\|}\sum_{j=k}^{k+l}\alpha_{j}w_{j}\big{\|}_{\infty}:l\geq
0\text{ and }\sum_{j=k}^{k+l}\alpha_{j}\leq 8m+1\Big{\\}}=0\,,$

almost surely. In particular, letting $\varepsilon\in(0,1)$, there almost-
surely exists $K$ (which depends on the realisation of the martingale noise)
such that

$\displaystyle\sup\Big{\\{}\big{\|}\sum_{j=k}^{k+l}\alpha_{j}w_{j}\big{\|}_{\infty}:l\geq
0\text{ and }\sum_{j=k}^{k+l}\alpha_{j}\leq 8m+1\Big{\\}}<\varepsilon$

for all $k\geq K$, and further such that $\alpha_{k}<1$ for all $k\geq K$.

Let us additionally take $\bar{M}\geq M$ such that $\delta\bar{M}\geq
4(8m+1)$. Suppose that for some $k\geq K$,
$\|\theta_{k}\|_{\infty}\geq\bar{M}+(8m+1)$. Let $l$ be minimal such that
$\sum_{j=k}^{k+l}\alpha_{j}>8m$. Then we have
$\|\theta_{k+j}-\theta_{k}\|_{\infty}\leq 8m+1$ for all $0\leq j\leq l$, and
so $\|\theta_{k+j}\|_{\infty}\geq\bar{M}$ for all $0\leq j\leq l$. Further, if
$\theta_{k}(x,i)$ satisfies
$\theta_{k}(x,i)>\|\theta_{k}\|_{\infty}(1-\delta)+2(8m+1)$, then we have

$\displaystyle\phantom{=}\theta_{k+j}(x,i)$
$\displaystyle\geq\theta_{k}(x,i)-(8m+1)$
$\displaystyle\geq(1-\delta)\|\theta_{k}\|_{\infty}+(8m+1)$
$\displaystyle\geq(1-\delta)(\|\theta_{k+j}\|_{\infty}-(8m+1))+(8m+1)$
$\displaystyle=(1-\delta)\|\theta_{k+j}\|_{\infty}\,,$

so $I^{+}_{k+j}(x,i)$ holds for all $0\leq j\leq l$, and hence

$\displaystyle|\theta_{k+l+1}(x,i)|\leq\|\theta_{k}\|_{\infty}-\sum_{j=k}^{k+l}\alpha_{j}\times
1/(4m)+\|\sum_{j=k}^{k+l}\alpha_{j}w_{j}\|_{\infty}<\|\theta_{k}\|_{\infty}-2+\varepsilon<\|\theta_{k}\|-1\,.$
(22)

Similarly, if $\theta_{k}(x,i)<-\|\theta_{k}\|_{\infty}(1-\delta)-2(8m+1)$,
$I^{-}_{k+j}(x,i)$ holds for all $0\leq j\leq l$, and we reach the same
conclusion as in Equation (22). Finally, if
$|\theta_{k}(x,i)|\leq\|\theta_{k}\|_{\infty}(1-\delta)+2(8m+1)$, then since
$\delta\|\theta_{k}\|_{\infty}>\delta\bar{M}$, we have
$|\theta_{k}(x,i)|\leq\|\theta_{k}\|_{\infty}-2(8m+1)$, and hence
$|\theta_{k+l+1}(x,i)|\leq\|\theta_{k}\|_{\infty}-(8m+1)$. Putting these
components together, we have

$\displaystyle\|\theta_{k+l+1}\|_{\infty}<\|\theta_{k}\|_{\infty}+1\,,\text{
and }\max_{0\leq j\leq
l}\|\theta_{k+j}\|_{\infty}\leq\|\theta_{k}\|_{\infty}+(8m+1)\,,$

as required to establish boundedness.

###  A.4 Proof of Proposition 5.11

We first state and prove a useful lemma that allows us to compare QDP fixed
points for different values of $\lambda$. Throughout this section, we will
adopt the shorthand $\theta^{\lambda}$ for $\hat{\theta}^{\pi}_{\lambda}$.

######  Lemma A.1.

Let $\lambda,\lambda^{\prime}\in[0,1]^{\mathcal{X}\times[m]}$. Then we have

$\displaystyle\|\theta^{\lambda}-\theta^{\lambda^{\prime}}\|_{\infty}\leq
C\|\lambda-\lambda^{\prime}\|_{\infty}\,,$

where $C$ is a constant depending only on the reward distributions of the MDP
and $\gamma$.

###### Proof.

By the triangle inequality, we have

$\displaystyle\|\theta^{\lambda}-\theta^{\lambda^{\prime}}\|_{\infty}$
$\displaystyle\leq\|\theta^{\lambda}-\Pi^{\lambda^{\prime}}\mathcal{T}^{\pi}\theta^{\lambda}\|_{\infty}+\|\Pi^{\lambda^{\prime}}\mathcal{T}^{\pi}\theta^{\lambda}-\theta^{\lambda^{\prime}}\|_{\infty}$
$\displaystyle=\|\Pi^{\lambda}\mathcal{T}^{\pi}\theta^{\lambda}-\Pi^{\lambda^{\prime}}\mathcal{T}^{\pi}\theta^{\lambda}\|_{\infty}+\|\Pi^{\lambda^{\prime}}\mathcal{T}^{\pi}\theta^{\lambda}-\Pi^{\lambda^{\prime}}\mathcal{T}^{\pi}\theta^{\lambda^{\prime}}\|_{\infty}$
$\displaystyle\leq\|(\Pi^{\lambda}-\Pi^{\lambda^{\prime}})\mathcal{T}^{\pi}\theta^{\lambda}\|_{\infty}+\gamma\|\theta^{\lambda}-\theta^{\lambda^{\prime}}\|_{\infty}$
$\displaystyle\implies\|\theta^{\lambda}-\theta^{\lambda^{\prime}}\|_{\infty}$
$\displaystyle\leq\frac{1}{1-\gamma}\|(\Pi^{\lambda}-\Pi^{\lambda^{\prime}})\mathcal{T}^{\pi}\theta^{\lambda}\|_{\infty}\,.$

Now we aim to bound $\|\theta^{\lambda}\|_{\infty}$, and hence the term on the
right-hand side above. Note that in general for a mixture distribution
$\nu=\sum_{i=1}^{n}p_{i}\nu_{i}$, we have
$F^{-1}_{\nu}(\tau)\geq\min\\{F^{-1}_{\nu_{i}}(\tau):i=1,\ldots,n\\}$, since

$\displaystyle\mathbb{P}_{Z\sim\nu}(Z\leq\min\\{F^{-1}_{\nu_{i}}(\tau_{n}):i=1,\ldots,n\\})$
$\displaystyle=\sum_{i=1}^{n}p_{i}\mathbb{P}_{Z_{i}\sim\nu_{i}}(Z_{i}\leq\min\\{F^{-1}_{\nu_{i}}(\tau):i=1,\ldots,n\\})$
$\displaystyle\leq\sum_{i=1}^{n}p_{i}\mathbb{P}_{Z_{i}\sim\nu_{i}}(Z_{i}\leq
F^{-1}_{\nu_{i}}(\tau))$ $\displaystyle\leq\tau\,.$

Thus, it follows that the $\nicefrac{{1}}{{2m}}$ quantile of
$\mathcal{T}^{\pi}\theta^{\lambda}(x)$ is at least as great as

$\displaystyle\min_{x}F_{\mathcal{R}^{\pi}(x)}^{-1}(\nicefrac{{1}}{{2m}})-\gamma\|\theta^{\lambda}\|_{\infty}\,.$

By analogous reasoning, we obtain that
$\bar{F}^{-1}_{(\mathcal{T}^{\pi}\theta^{\lambda})(x)}(\nicefrac{{2m-1}}{{2m}})$
is no greater than

$\displaystyle\max_{x}\bar{F}_{\mathcal{R}^{\pi}(x)}^{-1}(\nicefrac{{2m-1}}{{2m}})+\gamma\|\theta^{\lambda}\|_{\infty}\,.$

From these facts, it follows that

$\displaystyle\|\theta^{\lambda}\|_{\infty}\leq\frac{1}{1-\gamma}\max\left(|\min_{x}F_{\mathcal{R}^{\pi}(x)}^{-1}(\nicefrac{{1}}{{2m}})|,|\max_{x}\bar{F}_{\mathcal{R}^{\pi}(x)}^{-1}(\nicefrac{{2m-1}}{{2m}})|\right)\,,$

and hence

$\displaystyle\|(\Pi^{\lambda}-\Pi^{\lambda^{\prime}})\mathcal{T}^{\pi}\theta^{\lambda}\|_{\infty}\leq
C\|\lambda-\lambda^{\prime}\|_{\infty}\,,$

as required for the statement of the result. ∎

We now turn to the proof of Proposition 5.11. First, we observe that the
infimum over $\lambda$ in Equation (19) is attained, since Lemma A.1
establishes that $\lambda\mapsto\theta^{\lambda}$ is continuous (in fact
Lipschitz), and $[0,1]^{\mathcal{X}\times[m]}$ is compact. We therefore have
that $L$ is continuous, non-negative, and takes on the value $0$ only on the
set of fixed points
$\\{\theta^{\lambda}:\lambda\in[0,1]^{\mathcal{X}\times[m]}\\}$.

For the decreasing property, let $(\vartheta_{t})_{t\geq 0}$ be a solution to
the differential inclusion in Equation (17), and as in Definition 5.4, let
$g:[0,\infty)\rightarrow\mathbb{R}^{\mathcal{X}\times[m]}$ satisfy

$\displaystyle\vartheta_{t}=\int_{0}^{t}g_{s}\;\mathrm{d}s\,,$ (23)

with $g_{t}(x,i)\in H^{\pi}_{x,i}(\vartheta_{t})$ for all $(x,i)$, and for
almost all $t\geq 0$, where we have introduced the notation

$\displaystyle
H^{\pi}_{x,i}(\theta)=[\tau_{i}-F_{(\mathcal{T}^{\pi}\theta)(x)}(\theta(x,i)),\tau_{i}-F_{(\mathcal{T}^{\pi}\theta)(x)}(\theta(x,i)-)]\,.$

As in the proof of Proposition 5.10, we will show that $L(\vartheta_{t})$ is
locally decreasing outside of the fixed point set, which is enough for the
global decreasing property. Further, by continuity of $L(\vartheta_{t})$, it
is enough to show this property for almost all $t\geq 0$. We will therefore
consider a value of $t\geq 0$ at which the above inclusion for $g_{t}$ holds.

Let $\overline{\lambda}$ attain the minimum in the definition of
$L(\vartheta_{t})$. Write $\theta^{\overline{\lambda}}$ for the corresponding
fixed point for conciseness, and let $(x,i)$ be a $\overline{\lambda}$-_argmax
index with respect to $\vartheta_{t}$_; a state-particle pair achieving the
maximum in the definition of the norm
$\|\vartheta_{t}-\theta^{\overline{\lambda}}\|_{\infty}$. First, we consider
the cases where $H^{\pi}_{x,i}(\vartheta_{t})$ is _not_ a singleton. Now, if
$0\in H^{\pi}_{x,i}(\vartheta_{t})$, then we have
$(\Pi^{\lambda}\mathcal{T}^{\pi}\vartheta_{t})(x,i)=\vartheta_{t}(x,i)$, and
with the same logic as above, we have
$\vartheta_{t}=\theta^{\overline{\lambda}}$, and hence $\vartheta_{t}$ is in
the fixed point set, and $L(\vartheta_{t})$ is constant. If $0\not\in
H^{\pi}_{x,i}(\vartheta_{t})$, then as in the proof of Proposition 5.10, it
can be shown that any element of $H^{\pi}_{x,i}(\vartheta_{t})$ has the same
sign as

$\displaystyle(\Pi^{\lambda}\mathcal{T}^{\pi}\vartheta_{t})(x,i)-\vartheta_{t}(x,i)\,.$
(24)

In the case of Proposition 5.10, continuity of the derivative then allowed us
to deduce that $|\vartheta_{t}(x,i)-\theta^{\overline{\lambda}}(x,i)|$ is
locally decreasing. Here, we require a related concept of continuity for the
set-valued map $\theta\mapsto H^{\pi}_{x,i}(\theta)$, namely that it is upper
semicontinuous (see, e.g., Smirnov, 2002); for a given
$\theta\in\mathbb{R}^{\mathcal{X}\times[m]}$ and any given $\varepsilon>0$,
there exists $\delta>0$ such that if
$\|\theta^{\prime}-\theta\|_{\infty}<\delta$, then
$H^{\pi}_{x,i}(\theta^{\prime})\subseteq\\{h+v:h\in
H^{\pi}_{x,i}(\theta)\,,|v|<\varepsilon\\}$. From this, it follows that any
element of $H^{\pi}_{x,i}(\vartheta_{t+s})$, for sufficiently small positive
$s$, has the same sign as the expression in Equation (24), and so from
Equation (23), we have that
$|\vartheta_{t}(x,i)-\theta^{\overline{\lambda}}(x,i)|$ is locally decreasing,
as required.

Now, when $H^{\pi}_{x,i}(\vartheta_{t})$ _is_ a singleton, if it is non-zero,
then by the same argument as in the proof of Proposition 5.10, the
corresponding element has the same sign as the expression in Equation (24),
and so as above, we conclude that
$|\vartheta_{t}(x,i)-\theta^{\overline{\lambda}}(x,i)|$ is locally decreasing.

Finally, the case where there exists an argmax index $(x,i)$ with
$H^{\pi}_{x,i}(\vartheta_{t})=\\{0\\}$ requires more care, and we will need to
reason about the effects of perturbing $\lambda$ to show that the Lyapunov
function is decreasing. For some intuition as to what the problem is, if
$H^{\pi}_{x,i}(\vartheta_{t+s})=\\{0\\}$ for small positive $s$, then the
coordinate $\vartheta_{t+s}(x,i)$ is static, as it lies on the flat region of
the CDF $F_{(\mathcal{T}^{\pi}\vartheta_{t+s})(x)}$ at level $\tau_{i}$, and
so the distance $|\vartheta_{t+s}(x,i)-\theta^{\overline{\lambda}}(x,i)|$ is
not decreasing. We explain how to deal with this case below.

####  A.4.1 Perturbative argument

We introduce the notation $J_{0}\subseteq\mathcal{X}\times[m]$ for the set of
$\overline{\lambda}$-argmax indices with respect to $\vartheta_{t}$. Assuming
that $\|\vartheta_{t}-\theta^{\overline{\lambda}}\|_{\infty}$ is not locally
decreasing, it must be locally constant (it cannot increase, by the arguments
above). Now consider $s>0$ sufficiently small so that (i) no coordinates not
in $J_{0}$ can be a $\overline{\lambda}$-argmax index with respect to
$\vartheta_{t+s}$, so that $J$, the set of $\overline{\lambda}$-argmax indices
with respect to $\vartheta_{t+s}$, satisfies $J\subseteq J_{0}$, (ii) all
indices $(x,i)\in J$ satisfy $H^{\pi}_{x,i}(\vartheta_{t+u})=\\{0\\}$ for all
$u\in[0,2s]$.

We will now demonstrate the existence of a parameter
$\lambda^{\prime}\in[0,1]^{\mathcal{X}\times[m]}$ such that
$\|\vartheta_{t+s}-\theta^{\lambda^{\prime}}\|_{\infty}<\|\vartheta_{t}-\theta^{\overline{\lambda}}\|_{\infty}$,
which establishes the locally decreasing property of the Lyapnuov function, as
required. To do so, we introduce a modification of the fixed point map
$\lambda\mapsto\theta^{\lambda}$. Letting $\mu\in\mathbb{R}^{J}$, and defining
$\overline{\lambda}[\mu]\in\mathbb{R}^{\mathcal{X}\times[m]}$ to be the
replacement of the $J$ coordinates of $\overline{\lambda}$ with the
corresponding coordinates of $\mu$, we consider the map

$\displaystyle
h_{\overline{\lambda}}:[0,1]^{J}\rightarrow\mathbb{R}^{J}\,,\quad
h(\mu)=\mathrm{P}_{J}\theta^{\overline{\lambda}[\mu]}\,,$

where
$\mathrm{P}_{J}:\mathbb{R}^{\mathcal{X}\times[m]}\rightarrow\mathbb{R}^{J}$
extracts the $J$ coordinates. At an intuitive level, this map allows us to
study the effect of perturbing the $J$ coordinates of $\overline{\lambda}$ on
the corresponding coordinates of the fixed point.

####  A.4.2 Case 1: $\overline{\lambda}_{J}$ is in the interior of $[0,1]^{J}$

We now first consider the case where $\overline{\lambda}_{J}$, the $J$
coordinates of $\overline{\lambda}$, lies in the interior of $[0,1]^{J}$, that
is $(0,1)^{J}$. By Lemma A.1, $h_{\overline{\lambda}}$ is continuous, since it
is the composition of the continuous maps $\mu\mapsto\overline{\lambda}[\mu]$,
$\lambda\mapsto\theta^{\lambda}$, and $\theta\mapsto\mathrm{P}_{J}\theta$. It
is also injective in a neighbourhood of $\overline{\lambda}_{J}$. This can be
seen by noting first that the fixed points $\theta^{\overline{\lambda}[\mu]}$
are distinct for distinct values of $\mu$ sufficiently close to
$\overline{\lambda}_{J}$; if $\mu\not=\mu^{\prime}$ are each sufficiently
close to $\overline{\lambda}_{J}$, then we have

$\displaystyle\Pi^{\overline{\lambda}[\mu^{\prime}]}\mathcal{T}^{\pi}\theta^{\overline{\lambda}[\mu]}\not=\Pi^{\overline{\lambda}[\mu]}\mathcal{T}^{\pi}\theta^{\overline{\lambda}[\mu]}=\theta^{\overline{\lambda}[\mu]}\,,$

where the inequality follows from the fact that since
$\theta^{\overline{\lambda}[\mu]}$ is continuous in $\mu$, for $\mu$
sufficiently close to $J$ there is a flat region of
$F_{(\mathcal{T}^{\pi}\theta^{\overline{\lambda}[\mu]})(x)}$ at level
$\tau_{i}$, for any $(x,i)\in J$. To complete the injectivity argument, we
cannot have
$\mathrm{P}_{J}\theta^{\overline{\lambda}[\mu^{\prime}]}=\mathrm{P}_{J}\theta^{\overline{\lambda}[\mu]}$
if
$\theta^{\overline{\lambda}[\mu^{\prime}]}\not=\theta^{\overline{\lambda}[\mu]}$,
as the contraction maps $\Pi^{\overline{\lambda}[\mu]}\mathcal{T}^{\pi}$ and
$\Pi^{\overline{\lambda}[\mu]}\mathcal{T}^{\pi}$ are equal on coordinates not
in $J$, and these two maps would therefore have the same fixed point, a
contradiction.

We may now appeal to the invariance of domain theorem (Brouwer, 1912) to
deduce that since $h_{\overline{\lambda}}$ is a continuous injective map
between an open subset of $[0,1]^{J}$ containing $\overline{\lambda}_{J}$
(here we are using the assumption that $\overline{\lambda}_{J}$ lies in the
interior of $[0,1]^{J}$) and the Euclidean space $\mathbb{R}^{J}$ of equal
dimension, it is an open map on this domain; that is, it maps open sets to
open sets. Hence, we can perturb $\theta^{\overline{\lambda}}$ in the $J$
coordinates in any direction we want by locally modifying the $J$ coordinates
of $\overline{\lambda}$. In particular, we can move all $J$ coordinates of
$\theta^{\overline{\lambda}}$ closer to those of
$(\vartheta_{t+s}(x,i):(x,i)\in J)$. Let
$\lambda^{\prime}\in(0,1)^{\mathcal{X}\times[m]}$ be such a modification of
$\overline{\lambda}$, taken to be close enough to $\overline{\lambda}$ so that
all coordinates outside $J$ have sufficiently small perturbations so that they
cannot be $\lambda^{\prime}$-argmax indices with respect to $\vartheta_{t+s}$.
We then have that
$\|\vartheta_{t+s}-\theta^{\lambda^{\prime}}\|_{\infty}<\|\vartheta_{t}-\theta^{\lambda}\|_{\infty}$,
as required.

####  A.4.3 Case 2: $\overline{\lambda}_{J}$ is on the boundary of $[0,1]^{J}$

In the more general case when $\overline{\lambda}_{J}$ may lie on the boundary
of $[0,1]^{J}$, we can apply the same argument to an extension of the function
$h_{\overline{\lambda}}$, by increasing its domain from $[0,1]^{J}$ to an open
neighbourhood of this domain in $\mathbb{R}^{J}$. We define this extension
simply by extending the definition of $\Pi^{\lambda}$ in Equation (12) to
allow coordinates of $\lambda$ to lie outside the range $[0,1]$. We lose the
non-expansiveness of $\Pi^{\lambda}$ (in $L^{\infty}$) under this extension,
but if $\lambda_{\textrm{min}},\lambda_{\textrm{max}}$ are the minimum and
maximum coordinates of $\lambda$, respectively, it is easy verified (by
modifying the proof of Proposition 4.3) that $\Pi^{\lambda}$ is
$\max(1-\lambda_{\textrm{min}},\lambda_{\textrm{max}})$-Lipschitz, and so if
we extend the function to a domain where
$\lambda_{\textrm{max}},1-\lambda_{\textrm{min}}\leq\gamma^{-1/2}$, the
composition $\Pi^{\lambda}\mathcal{T}^{\pi}$ is a $\gamma^{1/2}$-contraction
in $L^{\infty}$, and hence has a unique fixed point $\theta^{\lambda}$.

By the same arguments as above, the extended map $h_{\overline{\lambda}}$ is
continuous and injective in a neighbourhood of $\overline{\lambda}_{J}$ on
this extended domain, and hence we may again apply the invariance of domain
theorem to obtain that $h_{\overline{\lambda}_{\setminus J}}$ is locally
surjective around $\overline{\lambda}_{J}$. However, since
$\overline{\lambda}_{J}$ lies on the boundary of the original domain, we must
additionally check that we can perturb $\overline{\lambda}_{J}$ to obtain
$\mu$ in such a way that we obtain the desired perturbation of
$\theta^{\overline{\lambda}}$, without the parameters $\mu$ leaving the set
$[0,1]^{J}$. To do this, we first rule out $\overline{\lambda}_{J}$ lying on
certain parts of the boundary.

######  Lemma A.2.

If $(x,i)\in J$ and $\vartheta_{t+s}(x,i)<\theta^{\overline{\lambda}}(x,i)$,
then $\overline{\lambda}(x,i)>0$. Similarly, if
$\vartheta_{t+s}(x,i)>\theta^{\overline{\lambda}}(x,i)$, then
$\overline{\lambda}(x,i)<1$.

###### Proof.

We prove the claim when
$\vartheta_{t+s}(x,i)<\theta^{\overline{\lambda}}(x,i)$; the other case
follows analogously. If $\overline{\lambda}(x,i)=0$, then since
$\vartheta_{t+s}(x,i)$ corresponds to the flat region at level $\tau_{i}$ of
the CDF $F_{(\mathcal{T}^{\pi}\vartheta_{t+s})(x)}$, we must have
$(\Pi^{\overline{\lambda}}\mathcal{T}^{\pi}\vartheta_{t+s})(x,i)\leq\vartheta_{t+s}(x,i)$
since $\overline{\lambda}(x,i)=0$, and so the chosen quantile at level
$\tau_{i}$ by the projection $\Pi^{\overline{\lambda}}$ is the left-most point
of this flat region. We therefore have

$\displaystyle\overline{w}_{\infty}(\Pi^{\overline{\lambda}}\mathcal{T}^{\pi}\vartheta_{t},\theta^{\overline{\lambda}})\geq|(\Pi^{\overline{\lambda}}\mathcal{T}^{\pi}\vartheta_{t})(x,i)-\theta^{\overline{\lambda}}(x,i)|\geq|\vartheta_{t}(x,i)-\theta^{\overline{\lambda}}(x,i)|=\overline{w}_{\infty}(\vartheta_{t},\theta^{\overline{\lambda}})\,,$

contradicting contractivity of $\Pi^{\overline{\lambda}}\mathcal{T}^{\pi}$
around $\theta^{\overline{\lambda}}$. ∎

We write
$v=\text{sign}((\vartheta_{t+s})_{J}-\theta^{\overline{\lambda}}_{J})\in\mathbb{R}^{J}$,
where the sign mapping is applied elementwise, and introduce the notation
$\mathrm{N}(v)=\\{\alpha\odot v:\alpha\in\mathbb{R}^{n}_{>0}\\}$ for the
(open) orthant containing the vector $v$. We are therefore seeking a
perturbation $\mu$ of $\overline{\lambda}_{J}$ such that
$\theta^{\overline{\lambda}[\mu]}_{J}$ lies in a direction in $\mathrm{N}(v)$
from $\theta^{\overline{\lambda}}_{J}$, and further such that the perturbation
to $\theta^{\overline{\lambda}}$ is sufficiently small that no index that was
not an argmax in $\|\vartheta_{t+s}-\theta^{\overline{\lambda}}\|_{\infty}$
can become one in
$\|\vartheta_{t+s}-\theta^{\overline{\lambda}[\mu]}\|_{\infty}$; under these
conditions, we have
$\|\vartheta_{t+s}-\theta^{\overline{\lambda}[\mu]}\|_{\infty}<\|\vartheta_{t+s}-\theta^{\overline{\lambda}}\|_{\infty}$,
as required. Lemma A.2 then guarantees that a (sufficiently small)
perturbation of $\overline{\lambda}_{J}$ in any direction in $\mathrm{N}(v)$
remains within $[0,1]^{J}$, so it is sufficient to show that a perturbation in
such a direction achieves the desired perturbation of
$\theta^{\overline{\lambda}}$.

Differentiability. Now, if the extended map $\lambda\mapsto\theta^{\lambda}$
is differentiable at $\overline{\lambda}$, then differentiating through the
fixed-point equation $\theta^{\lambda}=G(\lambda,\theta^{\lambda})$ (where we
write $G(\lambda,\theta)=\Pi^{\lambda}\mathcal{T}^{\pi}\theta$ for
conciseness) yields

$\displaystyle\nabla_{\lambda}\theta^{\lambda}=\partial_{\lambda}G(\lambda,\theta^{\lambda})+\partial_{\theta}G(\lambda,\theta^{\lambda})\nabla_{\lambda}\theta^{\lambda}\,;$

differentiability of $G$ in $\theta$ results from differentiability of the map
$\lambda\mapsto G(\lambda,\theta^{\lambda})$, and continuous differentiability
of $G$ in $\lambda$. Since $\theta\mapsto G(\lambda,\theta)$ is contractive in
$L^{\infty}$ with factor $\gamma^{1/2}$ (on the extended domain), and by
coordinatewise monotonicity of $\theta\mapsto G(\lambda,\theta)$, it follows
that $\partial_{\theta}G(\lambda,\theta^{\lambda})$ is non-negative and
strictly substochastic, with row $L^{1}$ norms bounded by $\gamma^{1/2}$, the
contraction factor for the extended set of contraction mappings. We remark as
a point of independent interest that this is a kind of Bellman equation for
$\nabla_{\lambda}\theta^{\lambda}$, with
$\partial_{\theta}G(\lambda,\theta^{\lambda})$ taking the role of the
transition matrix, and $\partial_{\lambda}G(\lambda,\theta^{\lambda})$ taking
the role of a collection of cumulants; in fact, the structure of
$\partial_{\theta}G(\lambda,\theta^{\lambda})$ coincides with the local
quantile back-up diagrams described in Example 6.5. We therefore have

$\displaystyle\nabla_{\lambda}\theta^{\lambda}=(I-\partial_{\theta}G(\lambda,\theta^{\lambda}))^{-1}\partial_{\lambda}G(\lambda,\theta^{\lambda})\,.$

By extracting the principal submatrix on the $J$ coordinates, we obtain a
derivative for $h_{\overline{\lambda}}(\overline{\lambda}_{J})$. The following
lemma is useful in reasoning about the structure of this principal submatrix.

######  Lemma A.3.

Let $Q_{1}\in\mathbb{R}^{n\times n}$ be strictly substochastic, and let
$K\subseteq[n]$. Then the principal submatrix on the $K$ coordinates of
$(I-Q_{1})^{-1}$ can be expressed as $(I-Q_{2})^{-1}$, with
$Q_{2}\in\mathbb{R}^{K\times K}$ strictly substochastic.

###### Proof.

We interpret $Q_{1}$ as the transition matrix of a Markov chain
$(Z_{t})_{t\geq 0}$ that includes a non-zero probability of termination at
each state. Each row of the matrix $(I-Q_{1})^{-1}$ is then the pre-
termination visitation measure associated with a particular initial state in
the Markov chain. Now let $Q_{2}$ be the strictly substochastic matrix defined
by

$\displaystyle Q_{2}(z_{1},z_{2})=\mathbb{P}((Z_{t})_{t\geq 0}$
$\displaystyle\text{ does not terminate before returning to }K,$
$\displaystyle\text{ first state on return is }z_{2}\mid Z_{0}=z_{1})\,.$

By construction, the pre-termination visitation distribution $(I-Q_{2})^{-1}$
is identical to the principal submatrix of $(I-Q_{1})^{-1}$ on the $K$
coordinates, as required. ∎

From Lemma A.3, we therefore obtain that $\nabla
h_{\overline{\lambda}}(\overline{\lambda}_{J})$ has the form

$\displaystyle\nabla
h_{\overline{\lambda}}(\overline{\lambda}_{J})=(I-Q)^{-1}D\,,$

with $D\in\mathbb{R}^{J\times J}$ diagonal, with positive elements on the
diagonal (from monotonicity of $\lambda\mapsto G(\lambda,\theta)$), with
$Q\in\mathbb{R}^{J\times J}$ strictly substochastic. The derivative is
therefore invertible, and we obtain the derivative of the inverse of the form

$\displaystyle\nabla
h^{-1}_{\overline{\lambda}}(\theta^{\lambda}_{J})=D^{-1}(I-Q)\,.$

From strict substochasticity of $Q$, and since $v\in\\{\pm 1\\}^{J}$, it
follows that for the desired perturbation direction $v$, we have

$\displaystyle\nabla
h^{-1}_{\overline{\lambda}}(\theta^{\lambda}_{J})v\overset{\text{S}}{=}v\,,$

and so $\nabla
h^{-1}_{\overline{\lambda}}(\theta^{\lambda}_{J})v\in\mathrm{N}(v)$, where the
equality of signs applies elementwise. Therefore, a perturbation of
$\overline{\lambda}_{J}$ in a direction in $\mathrm{N}(v)$ is achieved by a
sufficiently small perturbation of $\overline{\lambda}_{J}$ in a direction in
$\mathrm{N}(v)$, as required.

Non-differentiability. If $\lambda\mapsto\theta^{\lambda}$ is not
differentiable at $\overline{\lambda}$, we instead use techniques from non-
smooth analysis to complete the argument. First, since
$\lambda\mapsto\theta^{\lambda}$ is Lipschitz (by Lemma A.1), it is
differentiable almost everywhere by Rademacher’s theorem (Rademacher, 1919).
By adapting the argument made by Clarke (1976, Lemma 3), by Fubini’s theorem,
for almost all $\lambda_{\setminus
J}\in\mathbb{R}^{\mathcal{X}\times[m]\setminus J}$, the map
$\lambda\mapsto\theta^{\lambda}$ is differentiable at $(\lambda_{\setminus
J},\mu)$ for almost all $\mu$ with $(\lambda_{\setminus J},\mu)$ in the
extended domain. The map $(\lambda_{\setminus
J},\mu)\mapsto(\lambda_{\setminus J},h_{\lambda_{\setminus J}}(\mu))$ is
Lipschitz and locally injective around $\overline{\lambda}$, and hence maps
sufficiently small open neighbourhoods of $\overline{\lambda}$ to open
neighbourhoods of $(\overline{\lambda}_{\setminus
J},h_{\overline{\lambda}_{\setminus J}}(\overline{\lambda}_{J}))$. Further,
since each $h_{\lambda_{\setminus J}}$ is Lipschitz, and so absolutely
continuous, the inverse map $h^{-1}_{\lambda_{\setminus J}}$ is almost-
everywhere differentiable within such a neighbourhood. Following the analysis
of the differentiable case, we therefore deduce

$\displaystyle\nabla h^{-1}_{\lambda_{\setminus
J}}(\theta)v\overset{\text{S}}{=}v$

for almost all $\lambda_{\setminus J}$ in a ball $B$ around
$\overline{\lambda}_{\setminus J}$, and (for each such $\lambda_{\setminus
J}$) for almost all $\theta$ in the $L^{\infty}$ ball $B^{\prime}$ with centre
$h_{\overline{\lambda}_{\setminus J}}(\overline{\lambda}_{J})$ and radius
$\rho$, for some radius $\rho>0$. We further take $B$ and $B^{\prime}$ to be
of small enough radii so that this directional derivative is bounded on this
set, so that $h_{\lambda}^{-1}$ is locally Lipschitz on $B^{\prime}$ for each
$\lambda_{\setminus J}\in B$ (and hence absolutely continuous), and so that
for any $\theta\in B^{\prime}$, we have
$\text{sign}(\theta-(\vartheta_{t+s})_{J})=\text{sign}(\theta^{\overline{\lambda}}_{J}-(\vartheta_{t+s})_{J})$,
and so that for no $\mu$ in the preimage of $B^{\prime}$ under
$h_{\lambda_{\setminus J}}$ can have that $\|\theta^{(\lambda_{\setminus
J},\mu)}-\vartheta_{t+s}\|_{\infty}$ has new argmax coordinates outside of
$J$.

Let us consider $\tilde{\lambda}\in B$ at which the almost-everywhere
differentiability condition holds. By applying the same argument with Fubini’s
theorem, for almost all $\bar{\theta}$ in
$B(h_{\tilde{\lambda}}(\overline{\lambda}_{J}),\rho/4)$, the inverse
$h_{\tilde{\lambda}}^{-1}$ is differentiable almost everywhere on
$\\{\bar{\theta}+uv:u\in[0,\rho/2]\\}$.

Now, defining $\mu_{\tau}=h^{-1}_{\tilde{\lambda}}(\bar{\theta}+\tau v)$ for
$\tau\in[0,\rho/2]$, we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\tau}\mu_{\tau}=\nabla
h^{-1}_{\tilde{\lambda}}(\bar{\theta}+\tau v)v$

for almost all $\tau$, and by absolute continuity of
$h^{-1}_{\tilde{\lambda}}$, it follows that

$\displaystyle\mu_{\rho/2}=\mu_{0}+\int_{0}^{\rho/2}\frac{\mathrm{d}}{\mathrm{d}\tau}\mu_{\tau}\;\mathrm{d}\tau\,.$

Hence, $\mu_{\varepsilon}-\mu_{0}\in\mathrm{N}(v)$, and by construction
$h_{\tilde{\lambda}}(\mu_{\rho/2})=\bar{\theta}+\rho v/2$. By continuity of
$\lambda\mapsto\theta^{\lambda}$ and its inverse, and since $\tilde{\lambda}$
and $\bar{\theta}$ can be chosen above to be arbitrarily close to
$\overline{\lambda}_{\setminus J}$ and $h_{\overline{\lambda}_{\setminus
J}}(\overline{\lambda}_{J})$ respectively, we may consider a sequence of these
parameters converging to $\overline{\lambda}_{\setminus J}$ and
$h_{\overline{\lambda}_{\setminus J}}(\overline{\lambda}_{J})$, such that the
values of $\mu_{\rho/2}$ as constructed above also converge (by compactness),
and thus conclude the existence of $\bar{\mu}_{\rho/2}$ such that
$\bar{\mu}_{\rho/2}-\overline{\lambda}_{J}\in\overline{\mathrm{N}(v)}$, and
$h_{\overline{\lambda}}(\bar{\mu}_{\rho/2})=h_{\overline{\lambda}}(\overline{\lambda}_{J})+\rho
v/2$, as required.

###  A.5 Proof of Proposition 6.1

We begin with the observation that for any return-distribution function
$\eta\in\mathscr{P}([V_{\textsc{min}},V_{\textsc{max}}])$, for the projection
$\Pi^{\lambda}$ onto $\mathscr{F}_{\text{Q},m}$ (for any
$\lambda\in[0,1]^{\mathcal{X}\times[m]}$), we have

$\displaystyle\bar{w}_{1}(\Pi^{\lambda}\eta,\eta)\leq\frac{V_{\textsc{max}}-V_{\textsc{min}}}{2m}\,.$

Using this observation, we have

$\displaystyle\bar{w}_{1}(\hat{\eta}^{\pi}_{\lambda},\eta^{\pi})$
$\displaystyle\overset{(a)}{\leq}\bar{w}_{1}(\hat{\eta}^{\pi}_{\lambda},\mathcal{T}^{\pi}\hat{\eta}^{\pi}_{\lambda})+\bar{w}_{1}(\mathcal{T}^{\pi}\hat{\eta}^{\pi}_{\lambda},\eta^{\pi})$
$\displaystyle\overset{(b)}{=}\bar{w}_{1}(\Pi^{\lambda}\mathcal{T}^{\pi}\hat{\eta}^{\pi}_{\lambda},\mathcal{T}^{\pi}\hat{\eta}^{\pi}_{\lambda})+\bar{w}_{1}(\mathcal{T}^{\pi}\eta^{\pi},\mathcal{T}^{\pi}\eta^{\pi})$
$\displaystyle\overset{(c)}{\leq}\frac{V_{\textsc{max}}-V_{\textsc{min}}}{2m}+\gamma\bar{w}_{1}(\hat{\eta}^{\pi}_{\lambda},\eta^{\pi})\,.$

Here, (a) follows from the triangle inequality, (b) follows as
$\hat{\eta}^{\pi}_{\lambda}$, $\eta^{\pi}$ are fixed points of
$\Pi^{\lambda}\mathcal{T}^{\pi}$, $\mathcal{T}^{\pi}$, respectively, and (c)
follows from the application of the inequality at the beginning of the proof
and contractivity of $\mathcal{T}^{\pi}$. Rearranging then gives the desired
result.

###  A.6 Proof of Proposition 6.3

From the assumptions of the proposition, we have
$\hat{\eta}^{\pi}_{\lambda}=(\Pi^{\lambda}\mathcal{T}^{\pi})^{k}\eta^{\pi}$.
Then observe that following the argument for the proof of Proposition 6.1, we
have, for any $l\in\\{1,\ldots,k\\}$,

$\displaystyle\bar{w}_{1}((\Pi^{\lambda}\mathcal{T}^{\pi})^{l}\eta^{\pi},\eta^{\pi})\leq\gamma\bar{w}_{1}((\Pi^{\lambda}\mathcal{T}^{\pi})^{l-1}\eta^{\pi},\eta^{\pi})+\frac{V_{\textsc{max}}-V_{\textsc{min}}}{2m}\,.$

Chaining these inequalities yields the required statement.

##  Appendix B Implementations of quantile dynamic programming

Here, we describe two concrete implementations of QDP, which may be of
independent interest to the reader. Algorithm 3 (Bellemare et al., 2023)
describes an implementation when the reward distributions are available as
input to the algorithm as a list of outcomes and probabilities.

0: Quantile estimates $((\theta(x,i)_{i=1}^{m}:x\in\mathcal{X})$, Transition
and reward probabilities $(P^{\pi}(x^{\prime},r\mid
x):x,x^{\prime}\in\mathcal{X})$, Interpolation parameters
$\lambda\in[0,1]^{\mathcal{X}\times[m]}$.

1: for $x\in\mathcal{X}$ do

2: Set Targets as empty list {List of outcome/probability pairs}

3: for $x^{\prime}\in\mathcal{X}$ do

4: for $r\in\mathscr{R}$ do

5: for $j=1,\ldots,m$ do

6: Append $(r+\gamma\theta(x^{\prime},j),P^{\pi}(x^{\prime},r|x)/m)$ to
Targets

7: end for

8: end for

9: end for

10: Sort Targets ascending according to outcomes.

11: for $i=1,\ldots,m$ do

12: Find minimal outcome $q^{\prime}$ such that cumulative probability is
$\geq\nicefrac{{2i-1}}{{2m}}$.

13: Set $\theta^{\prime}(x,i)\leftarrow q^{\prime}$.

14: end for

15: end for

16: return $((\theta^{\prime}(x,i)_{i=1}^{m}:x\in\mathcal{X})$

Algorithm 3 Quantile dynamic programming (finitely-supported rewards)

Algorithm 4 makes use of a root-finding subroutine (such as
scipy.optimize.root_scalar), and can be used when the CDFs of the reward
distributions are available as input, and can be queried at individual points.
A common use case for this implementation is the case of Gaussian rewards.
Note that the root-finding subroutine is called on a monotonic scalar
function, and therefore strong guarantees can be given on the approximate
solution returned when the reward CDFs of the MDP are continuous.
Nevertheless, note that Algorithm 4 does not exactly implement the operator
$\Pi^{\lambda}\mathcal{T}^{\pi}$ due to this root-finding approximation error.

0: Quantile estimates $((\theta(x,i)_{i=1}^{m}:x\in\mathcal{X})$, Transition
probabilities $(P^{\pi}(x^{\prime}\mid x):x,x^{\prime}\in\mathcal{X})$, Reward
CDFs $(F_{\mathcal{R}^{\pi}(x)}:x\in\mathcal{X})$.

1: for $x\in\mathcal{X}$ do

2: Construct function
$\displaystyle\phi_{x}:t\mapsto\sum_{x^{\prime}\in\mathcal{X}}P^{\pi}(x^{\prime}|x)\sum_{j=1}^{m}F_{\mathcal{R}^{\pi}(x)}(t-\gamma\theta(x^{\prime},j))$

3: for $i=1,\ldots,m$ do

4: Use a scalar root-finding subroutine to find $\theta^{\prime}(x,i)$
approximately satisfying
$\displaystyle\phi_{x}(\theta^{\prime}(x,i))=\tau_{i}$

5: end for

6: end for

7: return $((\theta^{\prime}(x,i)_{i=1}^{m}:x\in\mathcal{X})$

Algorithm 4 Quantile dynamic programming (reward CDFs)

##  Appendix C Convergence of asynchronous QTD updates

Here, we describe the key considerations in extending our analysis to a proof
of convergence for asynchronous versions of QTD; our discussion follows the
approach of Perkins and Leslie (2013).

Step size restrictions. Typically, more restrictive assumptions on step sizes,
beyond the Robbins-Monro conditions, are required for asynchronous convergence
guarantees. See, for example, Assumption A2 of Perkins and Leslie (2013); note
that the typical Robbins-Monro step size schedule of $\alpha_{k}\propto
1/k^{\rho}$ for $\rho\in(1/2,1]$ satisfies these requirements.

Conditions on the sequence of states $(X_{k})_{k\geq 0}$ to be updated.
Additionally, different states are required to be updated “comparably often”;
assuming that $(X_{k})_{k\geq 0}$ forms an aperiodic irreducible time-
homogeneous Markov chain is sufficient, and this conditions holds when either
(i) $\pi$ generates such a Markov chain over the state space of the MDP of
interest, or (ii) when the states to be updated are sampled i.i.d. from a
fixed distribution supported on the entirety of the state space, amongst other
settings. See Assumption A4 of Perkins and Leslie (2013) for further details.

Modified differential inclusion. The QTD differential inclusion in Equation
(17) must be broadened to account for the possibility of different states
being updated with different frequencies, leading to a differential inclusion
of the form

$\displaystyle\vartheta_{t}(x,i)\in\\{\omega h:\omega\in(\delta,1]\,,h\in
H^{\pi}_{x,i}(\vartheta_{t})\\}\,,$

where $\delta$ represents a minimum relative update frequency for the state
$x$, derived from the conditions on $(X_{k})_{k\geq 0}$ described above.
Because of the structure of the Lyapunov function for the QTD DI in Equation
(19), it is readily verified that this remains a valid Lyapunov function for
this broader differential inclusion, for the same invariant set of QDP fixed
points.

## References

  * Aubin (1991) Jean-Pierre Aubin.  _Viability theory_.  Springer Birkhauser, 1991.
  * Aubin and Cellina (1984) Jean-Pierre Aubin and Arrigo Cellina.  _Differential inclusions: Set-valued maps and viability theory_.  Springer Science & Business Media, 1984.
  * Barth-Maron et al. (2018) Gabriel Barth-Maron, Matthew W. Hoffman, David Budden, Will Dabney, Dan Horgan, Dhruva TB, Alistair Muldal, Nicolas Heess, and Timothy Lillicrap.  Distributed distributional deterministic policy gradients.  In _Proceedings of the International Conference on Learning Representations_ , 2018.
  * Bellemare et al. (2013) Marc G. Bellemare, Yavar Naddaf, Joel Veness, and Michael Bowling.  The Arcade Learning Environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ , 47:253–279, June 2013.
  * Bellemare et al. (2017) Marc G. Bellemare, Will Dabney, and Rémi Munos.  A distributional perspective on reinforcement learning.  In _Proceedings of the International Conference on Machine Learning_ , 2017.
  * Bellemare et al. (2020) Marc G. Bellemare, Salvatore Candido, Pablo Samuel Castro, Jun Gong, Marlos C. Machado, Subhodeep Moitra, Sameera S. Ponda, and Ziyu Wang.  Autonomous navigation of stratospheric balloons using reinforcement learning.  _Nature_ , 588(7836):77–82, 2020.
  * Bellemare et al. (2023) Marc G. Bellemare, Will Dabney, and Mark Rowland.  _Distributional reinforcement learning_.  MIT Press, 2023.
  * Benaïm (1999) Michel Benaïm.  Dynamics of stochastic approximation algorithms.  In _Seminaire de probabilites XXXIII_ , pages 1–68. Springer, 1999\.
  * Benaïm and Faure (2013) Michel Benaïm and Mathieu Faure.  Consistency of vanishingly smooth fictitious play.  _Mathematics of Operations Research_ , 38(3):437–450, 2013.
  * Benaïm et al. (2005) Michel Benaïm, Josef Hofbauer, and Sylvain Sorin.  Stochastic approximations and differential inclusions.  _SIAM Journal on Control and Optimization_ , 44(1):328–348, 2005.
  * Benaïm et al. (2006) Michel Benaïm, Josef Hofbauer, and Sylvain Sorin.  Stochastic approximations and differential inclusions, Part II: Applications.  _Mathematics of Operations Research_ , 31(4):673–695, 2006.
  * Benveniste et al. (2012) Albert Benveniste, Michel Métivier, and Pierre Priouret.  _Adaptive algorithms and stochastic approximations_.  Springer Science & Business Media, 2012.
  * Bernardo et al. (2008) Mario Bernardo, Chris Budd, Alan Richard Champneys, and Piotr Kowalczyk.  _Piecewise-smooth dynamical systems: Theory and applications_.  Springer Science & Business Media, 2008.
  * Bertsekas and Tsitsiklis (1996) Dimitri P. Bertsekas and John N. Tsitsiklis.  _Neuro-dynamic programming_.  Athena Scientific, 1996.
  * Bodnar et al. (2020) Cristian Bodnar, Adrian Li, Karol Hausman, Peter Pastor, and Mrinal Kalakrishnan.  Quantile QT-Opt for risk-aware vision-based robotic grasping.  In _Robotics: Science and Systems_ , 2020.
  * Borkar (1998) Vivek S. Borkar.  Asynchronous stochastic approximations.  _SIAM Journal on Control and Optimization_ , 36(3):840–851, 1998.
  * Borkar (2008) Vivek S. Borkar.  _Stochastic approximation: A dynamical systems viewpoint_.  Springer, 2008.
  * Borkar and Meyn (2000) Vivek S. Borkar and Sean P. Meyn.  The ODE method for convergence of stochastic approximation and reinforcement learning.  _SIAM Journal on Control and Optimization_ , 38(2):447–469, 2000.
  * Brouwer (1912) Luitzen E. J. Brouwer.  Beweis der invarianz des $n$-imensionalen gebiets.  _Mathematische Annalen_ , 71(3):305–313, 1912\.
  * Brown (1951) George W. Brown.  Iterative solution of games by fictitious play.  _Act. Anal. Prod Allocation_ , 13(1):374, 1951\.
  * Clarke (1976) Francis H. Clarke.  On the inverse function theorem.  _Pacific Journal of Mathematics_ , 64(1):97–102, 1976.
  * Clarke et al. (1998) Francis H. Clarke, Yuri S. Ledyaev, Ronald J. Stern, and Peter R. Wolenski.  _Nonsmooth analysis and control theory_.  Springer Science & Business Media, 1998.
  * Dabney et al. (2018a) Will Dabney, Georg Ostrovski, David Silver, and Rémi Munos.  Implicit quantile networks for distributional reinforcement learning.  In _Proceedings of the International Conference on Machine Learning_ , 2018a.
  * Dabney et al. (2018b) Will Dabney, Mark Rowland, Marc G. Bellemare, and Rémi Munos.  Distributional reinforcement learning with quantile regression.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2018b.
  * Dayan (1992) Peter Dayan.  The convergence of TD($\lambda$) for general $\lambda$.  _Machine learning_ , 8(3-4):341–362, 1992.
  * Dayan and Sejnowski (1994) Peter Dayan and Terrence J. Sejnowski.  TD($\lambda$) converges with probability 1.  _Machine Learning_ , 14(3):295–301, 1994.
  * Fawzi et al. (2022) Alhussein Fawzi, Matej Balog, Aja Huang, Thomas Hubert, Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Francisco J. R. Ruiz, Julian Schrittwieser, Grzegorz Swirszcz, David Silver, Demis Hassabis, and Pushmeet Kohli.  Discovering faster matrix multiplication algorithms with reinforcement learning.  _Nature_ , 610(7930):47–53, 2022.
  * Filippov (1960) A. F. Filippov.  Differential equations with discontinuous right-hand side.  _Mat. Sb. (N.S.)_ , 51(93):99–128, 1960.
  * Gilbert and Weng (2016) Hugo Gilbert and Paul Weng.  Quantile reinforcement learning.  In _Proceedings of the Asian Workshop on Reinforcement Leanring_ , 2016.