Risk-Averse Approximate Dynamic Programming with
Quantile-Based Risk Measures

arXiv:1509.01920v4 [math.OC] 9 May 2017

Daniel R. Jiang and Warren B. Powell
May 10, 2017

Abstract
In this paper, we consider a finite-horizon Markov decision process (MDP) for which
the objective at each stage is to minimize a quantile-based risk measure (QBRM) of
the sequence of future costs; we call the overall objective a dynamic quantile-based
risk measure (DQBRM). In particular, we consider optimizing dynamic risk measures
where the one-step risk measures are QBRMs, a class of risk measures that includes
the popular value at risk (VaR) and the conditional value at risk (CVaR). Although
there is considerable theoretical development of risk-averse MDPs in the literature, the
computational challenges have not been explored as thoroughly. We propose datadriven and simulation-based approximate dynamic programming (ADP) algorithms to
solve the risk-averse sequential decision problem. We address the issue of inefficient
sampling for risk applications in simulated settings and present a procedure, based on
importance sampling, to direct samples toward the “risky region” as the ADP algorithm
progresses. Finally, we show numerical results of our algorithms in the context of an
application involving risk-averse bidding for energy storage.

1

Introduction

Sequential decision problems, in the form of Markov decision processes (MDPs), are most
often formulated with the objective of minimizing an expected sum of costs or maximizing
an expected sum of rewards [Puterman, 2014, Bertsekas and Tsitsiklis, 1996, Powell, 2011].
However, it is becoming more and more evident that solely considering the expectation is
insufficient as risk-preferences can vary greatly from application to application. Broadly
speaking, the expected value can fail to be useful in settings containing either heavy-tailed
distributions or rare, but high-impact events. For example, heavy-tailed distributions arise
frequently in finance (electricity prices are well-known to possess this feature; see Byström
[2005], Kim and Powell [2011]). In this case, the mean of the distribution itself may not
necessarily be a good representation of the randomness of the problem; instead, it is likely
useful to introduce a measure of risk on the tail of the distribution as well. The rare event
situation is, in a sense, the inverse case of the heavy-tail phenomenon, but it can also benefit
from a risk measure other than the expectation. To illustrate, certain problems in operations
research can be complicated by critical events that happen with small probability, such as
guarding against stock-outs and large back-orders in inventory problems (see Glasserman
and Liu [1996]) or managing the risk of the failure of a high-value asset (see Enders et al.
[2010]). In these circumstances, the merit of a policy might be measured by the number of
times that a bad event happens over some time horizon.

1

One way to introduce risk-aversion into sequential problems is to formulate the objective
using dynamic risk measures [Ruszczyński, 2010]. A rough preview, without formal definitions, of our optimization problem is as follows: we wish to find a policy that minimizes
risk, as assessed by a certain type of dynamic risk measure. The objective can be written
as


min ρα0 C1π + ρα1 C2π + · · · + ραT −1 (CTπ ) · · · ,
π∈Π

where Π is a set of policies, {Ctπ } are costs under policy π, and {ραt } are one-step risk
measures (i.e., components of the overall dynamic risk measure). Precise definitions are
given in the subsequent sections. We focus on the case where the objective at each stage
is to optimize a quantile-based risk measure (QBRM) of future costs; we call the overall
objective a dynamic quantile-based risk measure (DQBRM).
This paper makes the following contributions. First, we describe a new data-driven or
simulation-based ADP algorithm, called Dynamic-QBRM ADP, that is similar in spirit to
established asynchronous algorithms like Q-learning (see Watkins and Dayan [1992]) and
lookup table approximate value iteration (see, e.g., Bertsekas and Tsitsiklis [1996], Powell
[2011]), where one state is updated per iteration. The second contribution of the paper
is a companion sampling procedure to Dynamic-QBRM ADP, which we call risk-directed
sampling (RDS). As we describe above, when dealing with risk, there is a large class of
problems in which we are inherently dealing with rare, but very costly events. Broadly
speaking, the evaluation of a QBRM that is focused on the tail of the distribution (e.g.,
CVaR at, say, the 99% level) depends crucially on efficiently directing the algorithm toward
sampling these “risky” regions. In this part of the paper, we consider the question: is there
a way to learn, as the ADP algorithm progresses, the interesting values of the information
process to sample?
The paper is organized as follows. We first provide a literature review in Section 2.
In Section 3, we give our problem formulation, a brief introduction to dynamic risk measures, and the definition of a class of quantile-based risk measures. Next, we introduce
the algorithm for solving risk-averse MDPs in Section 4 and give a theoretical analysis in
Section 5. In Section 6, we discuss sampling issues and describe the companion sampling
procedure. We show numerical results on an example energy trading application in Section
7 and conclude in Section 8.

2

Literature Review

The theory of dynamic risk measures and the notion of time-consistency (see e.g. Riedel
[2004], Artzner et al. [2006], Cheridito et al. [2006]) is extended to the setting of sequential
optimization problems in Ruszczyński and Shapiro [2006a] and Ruszczyński [2010], in which
it is proved that any time-consistent dynamic risk measure can be written as compositions of
one-step conditional risk measures (these are simply risk measures defined in a conditional
setting, analogous to the conditional expectation for the traditional case). From this, a
Bellman recursion is obtained, becoming a familiar way of characterizing optimal policies.
Building on the theory of dynamic programming, versions of exact value iteration and policy
iteration are also developed in Ruszczyński [2010]. Later, in Çavus and Ruszczyński [2014],
these exact methods are analyzed in the more specific case of undiscounted transient models.
Under the assumption that we use one-step coherent risk measures, as axiomatized in
Artzner et al. [1999], the value functions of a risk-averse Markov decision process with a

2

convex terminal value function can be easily shown to satisfy convexity using the fact that
coherent risk measures are convex and monotone. Therefore, the traditional method of
stochastic dual dynamic programming (SDDP) of Pereira and Pinto [1991] for multistage,
risk-neutral problems, which relies on the convexity of value functions, can be adapted to
the risk-averse case. This idea is successfully explored in Philpott and de Matos [2012],
Shapiro et al. [2013], and Philpott et al. [2013], with applications to the large-scale problem
of hydro-thermal scheduling using one-step mean-CVaR (convex combination of mean and
CVaR) and one-step mean-upper semideviation risk measures. The main drawbacks of riskaverse SDDP are (1) the cost function must be linear in the state, (2) some popular risk
measures, such as value at risk (VaR), are excluded because they are not coherent and
do not imply convex value functions, and (3) the full risk measure (can be recast as an
expectation in certain instances) has to be computed at every iteration. Since no convexity
or linearity assumptions are made in this paper, we take an alternative approach from the
SDDP methods and instead assume the setting of finite state and action spaces, as in Qlearning. At the same time, because the default implementation of our approach does not
take advantage of structure, it is limited to smaller problems. Extensions to the methods
proposed in this paper for exploiting structure can be made by following techniques such as
those discussed in Powell et al. [2004], Nascimento and Powell [2009], and Jiang and Powell
[2015a].
Recursive stochastic approximation methods have been applied to estimating quantiles
in static settings (see Tierney [1983], Bardou and Frikha [2009], and Kan [2011]). Related
to our work, a policy gradient method for optimizing MDPs (with a risk-neutral objective)
under a CVaR constraint is given in Chow and Ghavamzadeh [2014]. All of these methods
are related to ours in the sense that the minimization formula [Rockafellar and Uryasev,
2002, Theorem 10] for CVaR is optimized with gradient techniques. In our multistage setting
with dynamic risk measures, which is also coupled with optimal control, there are some new
interesting complexities, including the fact that every new observation (or data point) is
generated from an imperfect distribution of future costs that is “bootstrapped” from the
previous estimate of the value function. This means that not only are the observations
inherently biased, but the errors compound over time – this was not the case for the static
setting considered in earlier work. Under reasonable assumptions, we analyze both the
almost sure convergence and convergence rates for our proposed algorithms.
Our risk-directed sampling procedure is inspired by adaptive importance sampling strategies from the literature, such as the celebrated cross-entropy method of Rubinstein [1999].
See, e.g., Al-Qaq et al. [1995], Bardou and Frikha [2009], Egloff and Leippold [2010], and
Ryu and Boyd [2015] for other similar approaches. The critical difference in our approach
is that in an ADP setting, we have the added difficulty of not being able to assume perfect
knowledge of the objective function; rather, our observations are noisy and biased. To our
knowledge, this is the first time an adaptive sampling procedure has been combined with
a value function approximation algorithm in the risk-averse MDP literature. The closest
paper is by Kozmı́k and Morton [2014], which considers an importance sampling approach
for policy evaluation.

3

Problem Formulation

In this section, we establish the setting of the paper. In particular, we describe the riskaverse model, introduce the concept of dynamic risk measures, and define a class of quantile-

3

based risk measures.

3.1

Model

We consider an MDP with a finite time-horizon, t = 0, 1, 2, . . . , T , where the last decision is
made at time t = T −1, so that the set of decision epochs is given by T = {0, 1, 2, . . . , T −1}.
Given a probability space (Ω, F, P), we define a discrete-time stochastic process {Wt }Tt=0 ,
with Wt ∈ W for all t, as the exogenous information process in the sequential decision
problem, where Wt is adapted to a filtration {Ft }Tt=0 , with {∅, Ω} = F0 ⊆ F1 ⊆ · · · ⊆
FT ⊆ F. We assume that all sources of randomness in the problem are encapsulated by
the process {Wt } and that it is independent across time. For computational tractability,
we work in the setting of finite state and action spaces. Let the state space be denoted S,
and let the action space be A, where |S| < ∞ and |A| < ∞. The set of feasible actions for
each state s ∈ S, written As , is a subset of A. The set U = {(s, a) ∈ S × A : a ∈ As } is
the set of all feasible state-action pairs. The stochastic process describing the states of the
system is {St }Tt=0 , where St is an Ft -measurable random variable taking values in S, and at
is a feasible action determined by the decision maker using St . Furthermore, let Zt denote
the space of Ft -measurable random variables and Zt,T = Zt × · · · × ZT .
We model the system using a transition function or system model S M : S × A × W → S,
which produces the next state St+1 given a current state St , action at , and an outcome
of the exogenous process Wt+1 : St+1 = S M (St , at , Wt+1 ). The cost for time t is given by
ct (St , at , Wt+1 ), where ct : S × A × W → R is the cost function. A policy is a sequence
of decision functions {Aπ0 , Aπ1 , . . . , AπT −1 } indexed by π ∈ Π, where Π is the index set of
all policies. Each decision function Aπt : S → A is a mapping from a state to a feasible
action, such that Aπt (s) ∈ As for any state s. Let the sequence of costs under a policy π be
represented by the process Ctπ for t = 1, 2, . . . , T , where
π
π
Ctπ = ct−1 (St−1
, Aπt−1 (St−1
), Wt ) ∈ Zt ,

where {Stπ } are the states visited while following policy π. Note that Ctπ refers to the cost
from time t − 1, but the index of t refers to its measurability: Ctπ depends on information
only known at time t.

3.2

Review of Dynamic Risk Measures

In this subsection, we briefly introduce the notion of a dynamic risk measure; for a more detailed treatment, see, e.g., Frittelli and Gianin [2004], Riedel [2004], Pflug and Ruszczyński
[2005], Boda and Filar [2006], Cheridito et al. [2006], and Acciaio and Penner [2011]. Our
presentation closely follows that of Ruszczyński [2010], which develops the theory of dynamic risk measures in the context of MDPs. First, a conditional risk measure is a mapping
ρt,T : Zt,T → Zt that satisfies the following monotonicity requirement: for X, Y ∈ Zt,T and
X ≤ Y (componentwise and almost surely), ρt,T (X) ≤ ρt,T (Y ).
Given a sequence of future costs Ct , . . . , CT , the intuitive meaning of ρt,T (Ct , . . . , CT ) is a
certainty equivalent cost (i.e., at time t, one is indifferent between incurring ρt,T (Ct , . . . , CT )
and the alternative of being subjected to the stream of stochastic future costs). See Rudloff
et al. [2014] for an in-depth discussion regarding the certainty equivalent interpretation in
the context of multistage stochastic models. A dynamic risk measure {ρt,T }Tt=0 is a sequence
of conditional risk measures ρt,T : Zt,T → Zt , which allows us to evaluate the future risk at
any time t using ρt,T . Of paramount importance to the theory of dynamic risk measures is
the notion of time-consistency, which says that if from the perspective of some future time
τ , one sequence of costs is riskier than another and the two sequences of costs are identical
4

from the present until τ , then the first sequence is also riskier from the present perspective
(see Ruszczyński [2010] for the full technical definition).
Other definitions of time-consistency can be found in the literature, e.g., Boda and Filar
[2006], Cheridito and Stadje [2009], and Shapiro [2009]. Though they may differ technically,
these definitions share the same intuitive spirit. Under the conditions:
ρt,T (0, . . . , 0) = 0

and

ρt,T (Ct , Ct+1 , . . . , CT ) = Ct + ρt,T (0, Ct+1 , . . . , CT ),

it is proven in Ruszczyński [2010] that for some one-step conditional risk measures ρt :
Zt+1 → Zt , a time-consistent, dynamic risk measure {ρt,T }Tt=0 can be expressed using the
following nested representation:

ρt,T (Ct , . . . , CT ) = Ct + ρt Ct+1 + ρt+1 (Ct+2 + · · · + ρT −1 (CT ) · · · ) ,

It is thus clear that we can take the reverse approach and define a time-consistent dynamic
risk measure by simply specifying a set of one-step conditional risk measures {ρt }Tt=0 . This
is a common method that has been used in the literature when applying the theory of
dynamic risk measures in practice (see, e.g., Philpott and de Matos [2012], Philpott et al.
[2013], Shapiro et al. [2013], Kozmı́k and Morton [2014], and Rudloff et al. [2014]).

3.3

Quantile-Based Risk Measures

In this paper, we focus on simulation techniques where the one-step conditional risk measure
belongs to a specific class of risk measures called quantile-based risk measures (QBRM).
Although the term quantile-based risk measure has been used in the literature to refer to
risk measures that are similar in spirit to VaR and CVaR (see, e.g., Dowd and Blake [2006],
Neise [2008], Sereda et al. [2010]), it has not been formally defined. First, let us describe
these two popular risk measures, which serve to motivate a more general definition for a
QBRM.
Also known as the quantile risk measure, VaR is a staple of the financial industry (see,
e.g., Duffie and Pan [1997]). Given a real-valued random variable X (representing a loss)
and a risk level α ∈ (0, 1), the VaR or quantile of X is defined to be

VaRα (X) = q α (X) = inf P(X ≤ u) ≥ α .
u

q α (X)

To simplify our notation, we use
in the remainder of this paper. It is well known that
VaR does not satisfy coherency Artzner et al. [1999], specifically the axiom of subadditivity,
an appealing property that encourages diversification. Despite this, several authors have
given arguments in favor of VaR. For example, Danielsson et al. [2005] concludes that in
practical situations, VaR typically exhibits subadditivity. Dhaene et al. [2006] and Ibragimov and Walden [2007] give other points of view on why VaR should not be immediately
dismissed as an effective measure of risk. A nested version of VaR for use in a multistage
setting is proposed in Cheridito and Stadje [2009], though practical implications have not
been explored in the literature.
CVaR is a coherent alternative to VaR and has been both studied and applied extensively
in the literature. Although the precise definitions may slightly differ, CVaR is also known
by names such as expected shortfall, average value at risk, or tail conditional expectation.
Given a general random variable X, the following characterization is given in Rockafellar
and Uryasev [2002]:
n
h

o
+ i
1
1
CVaRα (X) = inf u +
E (X − u)+ = q α (X) +
E X − q α (X)
.
u
1−α
1−α
5

Applications of risk-averse MDPs using dynamic risk measures have largely focused on
combining CVaR with expectation; once again, see Philpott and de Matos [2012], Philpott
et al. [2013], Shapiro et al. [2013], Kozmı́k and Morton [2014], and Rudloff et al. [2014].
For the purposes of this paper, we offer the following general definition of a QBRM
that allows dependence on more than one quantile; the definition includes the above two
examples as special cases.
Definition 1 (Quantile-Based Risk Measure (QBRM)). Let X be a real-valued random
variable. A quantile-based risk measure ρα can be written as the expectation of a function
of X and finitely many of its quantiles. More precisely, ρα takes the form
h
i
ρα (X) = E Φ X, q α1 (X), q α2 (X), . . . , q αm (X) ,
(3.1)
where α ∈ Rm is a vector of m risk levels, α1 , α2 , . . . , αm , and a function Φ : Rm+1 → R,
chosen so that ρα satisfies monotonicity, translation invariance, and positive homogeneity
(see Artzner et al. [1999] for the precise definitions and note that we interpret X as a
random loss or a cost).

Our definition of QBRMs is largely motivated by practical considerations. First, the
definition covers the two most widely used risk measures, VaR and CVaR, as special cases
under a single framework; in addition, the flexibility allows for the specification of more
sophisticated risk measures that may or may not be coherent. As previously mentioned,
there are situations where nonconvex (and thus, not coherent) risk measures are appropriate
[Dhaene et al., 2006]. Another motivation for this definition of a QBRM is that it allows us
to easily construct a risk measure such that VaRα (X) ≤ ρα (X) ≤ CVaRα (X), because, as
Belles-Sampera et al. [2014] points out, one issue with VaR is that it can underestimate large
losses, but at the same time, some practitioners of the financial and insurance industries
find CVaR to be too conservative.
We see that VaR is trivially a QBRM with Φ(X, q) = q. CVaR can also be easily written
1
as a QBRM, using the function Φ(X, q) = q + 1−α
(X − q)+ . Although our approach can
be applied to any risk measure of the form (3.1), we use the CVaR risk measure in the
empirical work of Section 7, due to its popularity in a variety of application areas.

3.4

Dynamic Quantile-Based Risk Measures

Notice that, so far, we have developed QBRMs in a “static” setting (the value of the risk
measure is in R) for simplicity. Given a random variable X ∈ Zt+1 and a risk level α ∈ (0, 1),
the conditional counterpart for the quantile is


qtα (X) = inf P X ≤ U | Ft ≥ α ∈ Zt .
U ∈Zt

Using this new definition, we can similarly extend the definition of a QBRM to the conditional setting by replacing (3.1) with
h
i

ραt (X) = E Φ X, qtα1 (X), qtα2 (X), . . . , qtαm (X) Ft ,

and replacing the required properties of monotonicity, translation invariance, and positive
homogeneity in Definition 1 with their conditional forms given in Ruszczyński [2010] (denoted therein by A2, A3, and A4). For the sake of notational simplicity, let us assume that
all parameters, i.e., m, α1 , . . . , αm , Φ are static over time, but we remark that an extension
to time-dependent (and even state-dependent) versions of the one-step conditional risk measure is possible. Let ρ̃αt be a (conditional) QBRM that measures tail risk. In applications,
6

a weighted combination of a tail risk measure with the traditional expectation ensures that
the resulting policies are not driven completely by the
of the cost distribution;
 tail behavior

we may use QBRMs of the form ραt (X) = (1 − λ) E X | Ft + λ ρ̃αt (X), where λ ∈ [0, 1].
Using one-step conditional risk measures as building blocks, we can define a dynamic
risk measure to be {ραt }Tt=0 , which we refer to as a dynamic quantile-based risk measure
(DQBRM). The dynamic risk measures obtained when ραt = VaRαt and ραt = CVaRαt
(the conditional forms of VaR and CVaR) are precisely the time-consistent risk measures
suggested in Cheridito and Stadje [2009] under the names composed value at risk and
composed conditional value at risk.

3.5

Objective Function

We are interested in finding optimal risk-averse policies under objective functions specified
using a DQBRM. The problem is


α
π
α
π
α
π
min ρ0 C1 + ρ1 C2 + · · · + ρT −1 (CT ) · · · .
(3.2)
π∈Π

The upcoming theorem, proven in Ruszczyński [2010], gives the Bellman-like optimality
equations for a risk-averse model. We state it under the assumption that the current period
contribution is random, differing slightly from the original statement. A point of clarification: the original theorem is proved within the setting where the one-step risk measures
satisfy conditional forms of the axioms of Artzner et al. [1999] for coherent risk measures.
In our setting, however, the QBRM ραt is only assumed to satisfy (conditional forms of)
monotonicity, positive homogeneity, and translation invariance, but not necessarily convexity. The crucial step of the proof given in Ruszczyński [2010] relies only on monotonicity
and an associated interchangeability property (see [Ruszczyński and Shapiro, 2006b, Theorem 7.1], [Ruszczyński, 2010, Theorem 2]). The assumption of convexity is therefore not
necessary for the following theorem.
Theorem 1 (Bellman Recursion for Dynamic Risk Measures, Ruszczyński [2010]). The
sequential decision problem (3.2) has optimal value functions given by

∗
(St+1 ) for all s ∈ S, t ∈ T ,
Vt∗ (s) = min ραt ct (s, at , Wt+1 ) + Vt+1
at ∈As

VT∗ (s)

= 0 for all s ∈ S.

The decision functions of an optimal policy π ∗ are given by

∗
∗
Aπt (s) ∈ arg min ραt ct (s, at , Wt+1 ) + Vt+1
(St+1 ) for all s ∈ S, t ∈ T ,
at ∈As

which map to a minimizing action of the optimality equation.

For computational purposes, we are interested in interchanging the minimization operator and the risk measure ραt and thus appeal to the state-action value function or Q-factor
formulation of the Bellman equation. Define the state-action value function
over the state
∗ (S
action pairs (s, a) ∈ U to be Q∗t (s, a) = ραt ct (s, a, Wt+1 ) + Vt+1
t+1 ) , for t ∈ T and let
Q∗T (s, a) = 0. Thus, the counterpart to the recursion in Theorem 1 is

Q∗t (s, a) = ραt ct (s, a, Wt+1 ) + 0 min Q∗t+1 (St+1 , a0 ) ,
(3.3)
a ∈ASt+1

with the minimization occurring inside of the risk measure.

7

3.6

Some Remarks on Notation

For simplicity, we henceforth refer to Q∗ simply as the optimal value function. Let d = |U|
and D = |U| (T + 1). We consider Q∗ to be a vector in RD with components Q∗t (s, a). We
also frequently use the notation Q∗t ∈ Rd for some t ≤ T , by which we mean Q∗ restricted to
the components Q∗t (s, a) for all (s, a) ∈ U. We adopt this system for any vector in RD (e.g.,
Q̄n , ui,∗ , and ūi,n to be defined later). The norms used in this paper are k · k1 , k · k2 , and
k · k∞ , the l1 -norm, the Euclidean norm, and the maximum norm, respectively. Example
usages of the latter two are
 X
1
2
∗
∗
2
kQt k2 =
Qt (s, a)
and kQ∗t k∞ = max |Q∗t (s, a)|.
(s,a)∈ U

(s,a)∈ U

The following naming convention is used throughout the paper and appendix: stochasi,n
h,n
tic processes denoted using , i.e., q,n
t+1 , t+1 , and t+1 , are conditionally unbiased noise
sequences and represent Monte Carlo sampling error. On the other hand, the processes deq,n
i,n
h,n
, ξt+1
, and ξt+1
, are biased noise and represent approximation error
noted using ξ, i.e., ξt+1
from using a value function approximation. For a vector v, diag(v) is the diagonal matrix
whose entries are the components of v. Lastly, for a nonnegative function f : R → R, its
support is represented by the notation supp(f ) = {x ∈ R : f (x) > 0}.

4

Algorithm

In this section, we introduce the risk-averse ADP algorithm for dynamic quantile-based risk
measures, which aims to approximate the value function Q∗ in order to produce near-optimal
policies.

4.1

Overview of the Main Idea

Like most ADP and reinforcement learning algorithms, the algorithm that we develop in
this paper to solve (3.2) is based on the recursive relationship of (3.3). The basic structure
for the algorithm is a time-dependent version of Q-learning or approximate value iteration
(see [Powell, 2011, Chapter 10] for a discussion). Recall the form of the QBRM:
h
i

ραt (X) = E Φ X, qtα1 (X), qtα2 (X), . . . , qtαm (X) Ft .

The main idea of our approach is to approximate the quantiles qtαi (X) and then combine the
approximations to form an estimate of the risk measure. In essence, every observation of the
exogenous information process (real or simulated data) can be utilized to give an updated
approximation of each of the m quantiles. A second step then takes the observation and the
quantile approximations to generate an refined approximation of the optimal value function
Q∗ . This type of logic is implemented using many concurrent stochastic gradient [Robbins
and Monro, 1951, Kushner and Yin, 2003] steps within a framework that walks through a
single forward trajectory of states and actions on each iteration.
It turns out that there is a convenient characterization of the quantile through the socalled CVaR minimization formula. Given a real-valued, integrable random variable X, a
risk level αi ∈ (0, 1), and u ∈ R, Rockafellar and Uryasev [2000] proves that


1
αi
+
q (X) ∈ arg min E u +
(X − u) .
(4.1)
1 − αi
u∈R

Although the main result of Rockafellar and Uryasev [2000] is that the optimal value of
8

the optimization problem gives the CVaRαi (X), the characterization of the quantile as the
minimizer is particularly useful for our purposes. It suggests the use of stochastic approximation or stochastic gradient descent algorithms [Robbins and Monro, 1951, Kushner and
Yin, 2003] to iteratively optimize (4.1).
With this intuition in mind, let us move back to the context of the MDP and define the
auxiliary variables ui,∗ ∈ RD , for each i ∈ {1, 2, . . . , m}, to be the αi -quantiles of the future
costs (recall that the quantiles αi are defined as an argument to our QBRM in Definition
1). The component at time t and state (s, a) is

αi
ui,∗
ct (s, a, Wt+1 ) + 0 min Q∗t+1 (St+1 , a0 ) ,
(4.2)
t (s, a) = q
a ∈ASt+1

for each i ∈ {1, 2, . . . , m}. Using (3.1), this allows us to take advantage of the equation
h
i
Q∗t (s, a) = E Φ ct (s, a, Wt+1 ) + 0 min Q∗t+1 (St+1 , a0 ), ut1,∗ (s, a), . . . , um,∗
(s,
a)
.
(4.3)
t
a ∈ASt+1

The relationship between Q∗ and ui,∗ is fundamental to our algorithmic approach, which
keeps track of mutually dependent approximations {ūi,n } and {Q̄n } to the optimal values
ui,∗ and Q∗ , respectively.

4.2

The Dynamic-QBRM ADP Algorithm

Before discussing the details, we need some additional notation. Clearly, at each time t, the
random quantity with which we are primarily concerned (and attempt to approximate) is the
future cost given the optimal value function Q∗ . Thus, we explicitly define its distribution
function for every (s, a):
h
i


Ft (x|s, a) = P ct s, a, Wt+1 + 0 min Q∗t+1 St+1 , a0 ≤ x .
a ∈ASt+1

Recall that d is the cardinality of the state-action space. Next, suppose ūit ∈ Rd is an
approximation of ui,∗
t and for each t and i ∈ {1, 2, . . . , m}, define the stochastic gradient
mapping ψti : Rd × Rd × W → Rd to perform the stochastic gradient computation:

ψti ūit , Q̄t+1 ,Wt+1 (s, a)
n
o


1
(4.4)
=1−
1 ct s, a, Wt+1 + 0 min Q̄t+1 St+1 , a0 ≥ ūit (s, a) .
a ∈ASt+1
1 − αi

To avoid confusion, we note that this is the stochastic gradient associated with the minimization formula (4.1) for CVaR, but this step is necessary for any QBRM, even if we are
not utilizing CVaR.
The second piece of notation we need is a specialized, stochastic version of the Bellman
operator to the risk-averse case: for each t, we define the mapping Ht : Rd × · · · × Rd × W →
Rd , with m + 1 arguments in Rd , to represent an approximation of the term within the
expectation of (4.3):

Ht ū1t , . . . ,ūm
t , Q̄t+1 , Wt+1 (s, a)



= Φ ct s, a, Wt+1 + 0 min Q̄t+1 St+1 , a0 , ū1t (s, a), . . . , ūm
t (s, a) .
a ∈ASt+1

Therefore, (4.3) can be rewritten using the stochastic Bellman operator Ht by replacing all
approximate quantities with their true values:
h
i

m,∗
∗
Q∗t (s, a) = E Ht u1,∗
,
.
.
.
,
u
,
Q
,
W
(s,
a)
.
(4.5)
t+1
t+1
t
t
9

The Dynamic-QBRM ADP algorithm that we describe in the next section consists of both
outer and inner iterations: for each outer iteration n, we step through the entire time horizon
of the problem t ∈ T . At time t, iteration n, the relevant quantities for our algorithms are
u,n
q,n
a state-action pair (Stn , ant ) ∈ U and two samples Wt+1
, Wt+1
∈ W from the distribution of
Wt+1 corresponding to the “two steps” of our algorithm, one for approximating the auxiliary
variables ui,∗ and the second for approximating the value function Q∗ . Figure 1 illustrates
the main idea behind the algorithm: we merge the results of m adaptive minimizations of
(4.1), corresponding to estimates of the m quantiles, into an estimate of the optimal value
function, Q̄nt . The estimate Q̄nt is then used to produce estimates of the relevant quantities
for the previous time period. Note that the m objective functions shown in the figure differ
only in their risk levels αi . The arrows on the curves indicate that the minimizations are
achieved via gradient descent steps.

E u+

1
1

↵1

(X

u)+


E u+

1
1

↵2

(X


E u+

u)+

1

1
(X
↵m

u)+

···
ū1,n
t

···

! u1,⇤
t

Q⇤t

ū2,n
t

ūm,n
t

! u2,⇤
t

h
i
Q̄nt ⇡ E Ht ⌅, ⌅, . . . , ⌅, ⌅, Wt+1

! um,⇤
t

Q⇤t+1

Q̄nt+1

Figure 1: Main Algorithm Idea

Now that we are in an algorithmic setting, we consider a new probability space (Ω, G, P),
where G = σ (Stn , ant , Wtu,n , Wtq,n ), n ≥ 0, t ≤ T . In order to describe the history of the
algorithms, we define:

Gtn = σ {(Sτk , akτ , Wτu,k , Wτq,k ), k < n, τ ≤ T } ∪ {(Sτn , anτ , Wτu,n , Wτq,n ), τ ≤ t} ,

for t ∈ T and n ≥ 1, with Gt0 = {∅, Ω} for all t ≤ T . We therefore have a filtration that
n
obeys Gtn ⊆ Gt+1
for t ∈ T and GTn ⊆ G0n+1 , coinciding precisely with the progression of
the algorithm. The random variables (Stn , ant ) are generated according to some sampling
u,n
q,n
policy (to be discussed later) while Wt+1
and Wt+1
are generated from the distribution of
the exogenous process Wt+1 .
Crucial to many ADP algorithms is the stepsize (or learning rate). In our case, we use
n
{γt } and {ηtn } for smoothing new observations with previous estimates, where γtn , ηtn ∈ Rd
for each t and n and are Gtn -measurable. The stepsize γtn is used to update our approximation
n
∗
of ui,∗
t while the stepsize ηt is used to update Qt ; see Algorithm 1. We articulate the
asynchronous nature of our algorithm by imposing the following condition on the stepsizes
(included in Assumption 1 of Section 5):
(s, a) 6= (Stn , ant ) =⇒ γtn (s, a) = ηtn (s, a) = 0,
which causes updates to only happen for states that we actually visit.
Stochastic approximation theory often requires a projection step onto a compact set
(giving bounded iterates) to ensure convergence [Kushner and Yin, 2003]. Hence, for each
10

t and (s, a), let Xtu (s, a) ⊆ R and Xtq (s, a) ⊆ R be compact intervals and let


Xtu = x ∈ Rd : x(s, a) ∈ Xtu (s, a) and Xtq = x ∈ Rd : x(s, a) ∈ Xtq (s, a) ,

be our projection sets at time t. The Euclidean projection operator to a set X is given by
the usual definition:
ΠX (y) = arg min ky − xk22 .
x∈X

These sets may be chosen arbitrarily large in practice and our first theoretical result (almost
sure convergence) will continue to hold. However, there is a tradeoff: if, in addition, we
want our convergence rate results to hold, then these sets also cannot be too large (see
Assumption 3).
The precise steps of Dynamic-QBRM ADP are given in Algorithm 1. A main characteristic of the algorithm is that sequences {ūi,n } and {Q̄n } are intertwined (i.e., depend on each
other). Consequently, there are multiple levels of approximation being used throughout the
steps of the algorithm. The theoretical results of the subsequent sections shed light onto
these issues.
Algorithm 1: Dynamic-QBRM ADP
Input: Initial estimates of the value function Q̄0 ∈ RD and quantiles ūi,0 ∈ RD for
i ∈ {1, 2, . . . , m}.
Stepsize rules γtn and ηtn for all n, t.
Output: Approximations {Q̄n } and {ūi,n }.

Set Q̄nT = 0 for all n.
for n = 1, 2, . . . do
2
Choose an initial state (S0n , an0 ).
for t = 0, 1, . . . , T − 1 do
u,n
q,n
3
Draw samples of the information process Wt+1
, Wt+1
∈ W.
1

4

5

6

7

5

Update auxiliary variable
approximations for i = 1, . . . , m:
n
o
u,n 
i,n
i,n−1
, Q̄n−1
ūt = ΠXtu ūt
− diag(γtn ) ψti ūi,n−1
.
t
t+1 , Wt+1

Compute an estimate of the future cost based on the current approximations:
q,n 
, . . . , ūm,n−1
, Q̄n−1
q̂tn = Ht ū1,n−1
t
t
t+1 , Wt+1 .

Update approximation
n of value function:
o
n
Q̄t = ΠXtq Q̄n−1
− diag(ηtn ) Q̄n−1
− q̂tn .
t
t

n
Choose next state (St+1
, ant+1 ).
end
end

Analysis of Convergence

In this section, we state and prove convergence theorems for Algorithm 1. First, we give an
overview of our analysis and the relationship to existing work.

5.1

A Preview of Results

The two main results of this section are: (1) the almost sure convergence of Dynamic-QBRM
ADP and (2) a convergence rate result under a particular sampling policy called ε-greedy.
11

The proof of almost sure convergence uses techniques from the stochastic approximation
literature [Kushner and Yin, 2003], which were applied to the field of reinforcement learning and Q-learning in Tsitsiklis [1994], Jaakkola et al. [1994] and Bertsekas and Tsitsiklis
[1996]. However, our algorithm differs from risk-neutral Q-learning in that it tracks multi2,n
m,n
ple quantities, ū1,n
, Q̄nt , over a horizon 1, 2, . . . , T . The intuition behind the
t , ūt , . . . , ūt
proof is that multiple “stochastic approximation instances” are pasted together in order to
obtain overall convergence of all relevant quantities. Accordingly, the interdependence of
various approximations means that in several parts of the proof, we require careful analysis of biased noise terms (or approximation error ) in addition to unbiased statistical error.
See, e.g., Kearns and Singh [1999], Even-Dar and Mansour [2004] and Azar et al. [2011],
for convergence rate results for standard Q-learning. The proof technique used to analyze
the high probability convergence rate of risk-neutral Q-learning in Even-Dar and Mansour
[2004] is based on the same types of stochastic approximation results that we utilize in this
paper.
Let us now make a few remarks regarding some simplifying assumptions made in this
paper. As proven
in [Rockafellar and Uryasev, 2002, Theorem 10], the set of minimizers

arg minu∈R E u + (1 − αi )−1 (X − u)+ is a nonempty, closed, and bounded interval for a
general X. We shall for ease of presentation, however, make assumptions (strictly increasing
and continuous cdf, Assumption 2(iii)) to guarantee that ui,∗
t (s, a) is the unique minimizer
0
when X is the optimal future cost ct (s, a, Wt+1 )+mina ∈ASt+1 Q∗t+1 (St+1 , a0 ) and that gradient computations to remain valid. This assumption is sufficient for almost sure convergence
(Theorem 2). To further examine the convergence rate of the algorithm (Theorem 3), we
must additionally have Assumption 3, which states that the density of the future cost exists
and is positive within the constraint sets Xtu — this provides us the technical condition of
strong convexity (discussed more in Section 5.3 below).
Since mina0 ∈ASt+1 Q∗t+1 (St+1 , a0 ) has a discrete distribution, the assumptions hold only
in certain situations: an obvious case is when the current stage cost has a density and is
independent of St+1 . For example, such a property holds when Wt+1 can be written as two
1 , W 2 ) where the current stage cost depends on W 1
independent components (Wt+1
t+1
t+1 and
2 . This model is relevant in a number of applicathe downstream state depends on Wt+1
tions; notable examples include multi-armed bandits [Whittle, 1980], shortest path problems
with random edge costs [Ryzhov and Powell, 2011], trade execution with temporary (and
temporally independent) price impact [Bertimas and Lo, 1998], and energy trading in twosettlement markets [Löhndorf et al., 2013]. Small algorithmic extensions (requiring more
complex notation) to handle the general case are possible, but the fundamental concepts
would remain unchanged. Hence, we will assume the cleaner setting for the purposes of this
paper.

5.2

Almost Sure Convergence

First, we discuss the necessary algorithmic assumptions, many of which are standard to the
field of stochastic approximation.
Assumption 1. For all (s, a) ∈ U and t ∈ T , suppose the following are satisfied:
(i) γtn (s, a) = γ̃tn−1 1{(s,a)=(Stn ,ant )} , for some γ̃tn−1 ∈ R that is Gtn−1 -measurable,

(ii) ηtn (s, a) = η̃tn−1 1{(s,a)=(Stn ,ant )} , for some η̃tn−1 ∈ R that is Gtn−1 -measurable,
(iii)

∞
X

n=1

γtn (s, a) = ∞,

∞
X

n=1

ηtn (s, a) = ∞

a.s.,

12

(iv)

∞
X

n=1

γtn (s, a)2 < ∞,

∞
X

n=1

ηtn (s, a)2 < ∞

a.s.,

(v) ∃ ε > 0, such that state sampling policy satisfies


n
≥ ε and P (S0n , an0 ) = (s, a) GTn−1 ≥ ε,
P (Stn , ant ) = (s, a) Gt−1

u
(vi) the projection sets are chosen large enough so that ui,∗
t ∈ Xt for each i ∈ {1, 2, . . . , m}
q
and Q∗t ∈ Xt .

Assumption 1(i) and (ii) represent the asynchronous nature of the algorithm, sending the
stepsize to zero whenever a state is not visited, while (iii) and (iv) are standard conditions on
the stepsize. Assumption 1(v) is an exploration requirement; by the Extended Borel-Cantelli
Lemma (see Breiman [1992]), sampling with this exploration requirement guarantees that
we will visit every state infinitely often with probability one. In particular, for the case
with an ε-greedy sampling policy (i.e., explore with probability ε, follow current policy
otherwise), this assumption holds. We discuss this policy in greater detail in Section 5.
Part (vi) is a technical assumption. The second group of assumptions that we present are
related to the problem parameters.
Assumption 2. The following hold:
(i) the risk-aversion function Φ : Rm+1 → R (from the QBRM within the one-step conditional risk measure ραt ) is Lipschitz continuous with constant LΦ > 0, i.e., for all
v, w ∈ Rm+1 , |Φ(v) − Φ(w)| ≤ LΦ kv − wk1 ,


(ii) ∃ Cmax > 0 such that E ct (s, a, Wt+1 )2 ≤ Cmax for all (s, a) ∈ U and t ∈ T ,
(iii) the distribution function Ft (x|s, a) is strictly increasing and Lipschitz continuous with
constant LF > 0, i.e.,
|Ft (x|s, a) − Ft (y|s, a)| ≤ LF |x − y|,
for all x, y ∈ R, (s, a) ∈ U, and t ∈ T .
As an example for Assumption 2(i), in the case of the QBRM being CVaR, where
1
1
(X − q)+ , it is easy to see that we can take LΦ = 1 + 1−α
. Assumption
Φ(X, q) = q + 1−α
2(ii) states that the second moment of the cost function is bounded. Assumption 2(iii)
u
and Assumption 1(vi) together imply that ui,∗
t (s, a) is the unique u ∈ Xt (s, a) such that
Ft (u|s, a) = αi .
Lemma 1. Under Assumptions 1 and 2, if Q̄nt+1 → Q∗t+1 almost surely, then the sequences
i,n
i,∗
of iterates ūi,n
t generated in Step 4 of Algorithm 1 satisfy ūt → ut almost surely for each
i ∈ {1, 2, . . . , m}.
Proof. We distinguish between two types of noise sequences (in n) for each t ∈ T , the statisi,n
d
d
tical error and the approximation error, denoted by i,n
t+1 ∈ R and ξt+1 ∈ R , respectively.
The definitions are


u,n 
i i,n−1
i,n
, Q∗t+1 , Wt+1
− E ψti ūi,n−1
, Q∗t+1 , Wt+1 ,
t
t+1 = ψt ūt
i,n
u,n 
u,n 
n−1
ξt+1
= ψti ūti,n−1 , Q̄t+1
, Wt+1
− ψti ūi,n−1
, Q∗t+1 , Wt+1
,
t

and we see that the random variable i,n
t+1 represents the error that the sample gradient
deviates from its mean, computed using the true future cost distribution (i.e., assuming we
i,n
have Q∗t+1 ). On the other hand, ξt+1
is the error between the two evaluations of ψti given
13

u,n
∗
the same sample Wt+1
, due only to the difference between Q̄n−1
t+1 and Qt+1 . Rearranging,
we have
 i i,n−1 ∗

u,n 
i,n
ψti ūti,n−1 , Q̄n−1
, Qt+1 , Wt+1 + i,n
t+1 , Wt+1 = E ψt ūt
t+1 + ξt+1 ,

which implies that the update given in Step 4 of Algorithm 1 can be rewritten as
n
h 
io

i,n−1
i,n
i,n
n
i i,n−1
∗
u ū
ūi,n
=
Π
−
diag(γ
)
E
ψ
ū
,
Q
,
W
+

+
ξ
.
t+1
X
t
t
t+1
t
t
t
t+1
t+1
t

Note that the term in the square brackets is a biased stochastic gradient and observe that
it is bounded (since ψti only takes two finite values). For the present inductive step at time
t, let us fix a state (s, a). It now becomes convenient for us to view ūi,n
t (s, a) as a stochastic
n }
n ⊆ G n+1 ⊆ G n+2 · · · ). It is clear
process in n, adapted to the filtration {Gt+1
(since
G
n≥0
t+1
t+1
t+1
that by the definition of i,n
(s,
a):
t+1


n−1
E i,n
(5.1)
t+1 (s, a) | Gt+1 = 0 a.s.

Therefore, nt+1 (s, a) are unbiased increments that can be referred to as martingale difference
noise. Before continuing, notice the following useful fact:
n−1
n−1
− Q∗t+1
(s, a) ≤ Q̄t+1
min Q∗t+1 (s, a) − min Q̄t+1

a∈As

a∈As

∞

.

(5.2)

The proof follows from min v = − max (−v) and max v − max w ≤ max |v − w|, where the
minimum and maximum
are taken over the components of some vectors v and w. Now, let
u,n 
i,n
n
M
(s, a) and using (5.2), we obtain
St+1 = S s, a, Wt+1 . Expanding the definition of ξt+1
i,n
ξt+1
(s, a)

 n
o

1
u,n 
n
=
1 ct s, a, Wt+1
, a0 ≥ ūi,n−1
(s,
a)
+ 0 min Q∗t+1 St+1
t
a ∈AS n
1 − αi
t+1
o
n

i,n−1
u,n 
n−1
n
0
(s, a)
− 1 ct s, a, Wt+1 + 0 min Q̄t+1 St+1 , a ≥ ūt
a ∈AS n

t+1

 n
o

u,n 
n
1 ct s, a, Wt+1
+ 0 min Q∗t+1 St+1
, a0 ≥ ūi,n−1
(s,
a)
t

1
≤
a ∈AS n
1 − αi
t+1
n


u,n
n−1
n
− Q∗t+1
− 1 ct s, a, Wt+1
, a0 ≥ ūi,n−1
(s, a) + Q̄t+1
+ 0 min Q∗t+1 St+1
t
a ∈AS n

t+1

o
.
∞

Using the shorthand Ft ( · ) = Ft ( ·|s, a) and taking the conditional expectation of both sides,
we get (almost surely)


 i,n


1
i,n−1
n−1
n−1
∗
E ξt+1
(s, a) | Gt+1
≤
Ft ūi,n−1
(s,
a)
+
Q̄
−
Q
(s, a)
t+1 ∞ − Ft ūt
t
t+1
1 − αi
LF
n−1
≤
Q̄t+1
− Q∗t+1 ∞ ,
(5.3)
1 − αi

where the second inequality follows by Assumption 2(iii). Since we assumed in the statement
of the lemma that Q̄nt+1 → Q∗t+1 almost surely, it must be the case that
 i,n

n−1
E ξt+1
(s, a) Gt+1
→ 0 a.s.,
(5.4)
and hence the noise “vanishes asymptotically” in expectation. Now, given the boundedness
of ψti (and hence, finite second moment of ψti ), the unbiasedness property (5.1), the vanishing
noise property (5.4), the stepsize and sampling properties of Assumption 1, and the unique-

14

ness of u∗t (s, a) from Assumption 2(iii), we can apply a classical theorem of stochastic api,∗
proximation, [Kushner and Yin, 2003, Theorem 2.4], to conclude that ūi,n
t (s, a) → ut (s, a)
almost surely for each i. Because we chose an arbitrary (s, a), this convergence holds for all
(s, a) ∈ U.
Theorem 2 (Almost Sure Convergence). Choose initial approximations Q̄0 ∈ RD and
u
ūi,0 ∈ RD for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt for all t ∈ T . Under
Assumptions 1–2, Algorithm 1 generates a sequence of iterates Q̄n that converge almost
surely to the optimal value function Q∗ .
Sketch of Proof: The idea of the proof is to induct backwards on t and repeatedly apply
Lemma 1. At each step, the induction hypothesis is that Q̄nt+1 → Q∗t+1 from which we obtain
i,∗
convergence of ūi,n
t → ut from Lemma 1. After making sure certain technical details are
satisfied, stochastic approximation theory allows us to show that Q̄nt → Q∗t . The full details
are given in Appendix A.

5.3

Convergence Rate

In this section, we discuss the convergence rate (in terms of the expected deviation to Q∗ )
of the procedure described in Algorithm 1. Because we are working with an asynchronous
algorithm where only one state per time period is visited every iteration, it is necessary for us
to specify the assumed state sampling policy for visiting states before deriving convergence
rate results. Due to Assumption 1(i)–(ii), we naturally must consider the sampling policy
and the stepsize sequences jointly, as is done in the upcoming proposition. We employ the
ε-greedy policy, a popular choice that balances exploitation and exploration using a tunable
parameter ε; it is defined as follows. For any iteration n > 0, choose (S0n , an0 ) independently
and uniformly at random. At iteration n, time t > 0, let
n
snt+1 = S M (Stn , ant , Wt+1
)

n−1 n
(st+1 , a)
and ant+1 = arg min Q̄t+1
a∈A

n
Xt+1

and define a Bernoulli random variable
with parameter 1 − εd that is independent
n
from Gt (i.e., the rest of the system). The next state to visit is selected by the rule
(

n = 1,

if Xt+1
snt+1 , ant+1
n
n
St+1 , at+1 =
(5.5)
Choose uniformly over U otherwise.

Note that Assumption 1(v) is clearly satisfied as each state is visited with probability at
least ε. We also choose our stepsize sequences with γ̃tn−1 and η̃tn−1 as deterministic harmonic
sequences (for ease of analysis) so that
γt
ηt
n n
γtn (s, a) =
1
and ηtn (s, a) = 1{(s,a)=(Stn ,ant )} ,
(5.6)
n {(s,a)=(St ,at )}
n

where γt > 0 and ηt > 0 are deterministic, time-dependent constants. Since neither γ̃tn−1 =
γtn /n nor η̃tn = ηtn /n depend on the history of visited states, they are known as centralized
learning rates (see Szepesvari and Littman [1996] for a discussion). The main difficulty
for the centralized stepsizes of (5.6) is that when the frequency of visits to a state decays
quickly enough, then Assumption 1(iii) and (iv) may not hold (for example, consider when
a state is only visited on iteration numbers that are powers of two). In fact, it is often not
immediately obvious when Assumption 1(iii) and (iv) are satisfied. The next proposition
shows that under the ε-greedy sampling policy, the states are visited often enough that the
assumption remains satisfied.
15

Proposition 1. Under the ε-greedy sampling policy given in (5.5) and the deterministic
harmonic stepsizes given in (5.6), Assumption 1(iii)–(iv) is satisfied.
Proof. See Appendix A.
For a technical reason needed to prove the convergence rate results, we need to constrain
the iterates of the algorithm to a region within the support of the distribution of future
costs (stated formally in Assumption 3 below). We shall see that when this assumption is
satisfied, we get a notion of strong convexity. To be more precise, observe that ∂ 2 E[u +
(X − u)+ /(1 − αi )]/∂u2 = fX (u)/(1 − α), where fX is the density of X. If u is constrained
to be within the support of the distribution of X (i.e., where fX > 0), then we are able
to lower bound the second derivative by some constant (that depends on the constraint
set), thereby attaining strong convexity within the region. This is useful for deriving the
convergence rate results (see Lemma 2); unfortunately, this condition is in general difficult
to check in practice.
Assumption 3. The density ft (x|s, a) = dFt (x|s, a)/dx exists and the stochastic approx
imation projection set is within the support of the density: Xtu (s, a) ⊆ supp ft ( ·|s, a) for
all (s, a) ∈ U and t ∈ T . Let

lf = min min x ∈ Xtu (s, a) : ft (x|s, a)
(s,a,t)

be a positive real number that lower bounds the density function over all t and (s, a) ∈ U.

We now provide a few lemmas that will be useful in establishing the final result. The
i,n−1
first lemma relates the error of ūi,n
) and the
t to the error in the last iteration (i.e., of ūt
n−1
).
error of the value function in the next time period (i.e., of Q̄t+1
Lemma 2. Under Assumptions 1(vi)–3, the ε-greedy sampling policy of (5.5), and the deterministic harmonic stepsizes given in (5.6), the sequence of approximations ūi,n
t generated
by Algorithm 1 satisfies, for any κ > 0,
i
h
i,∗ 2
E ūi,n
−
u
t
t
2
h
i
i
i h
γt
γt h n−1
γt2
2
i,∗ 2
∗
≤ 1−
2 ε Clf − κ CLF E ūi,n−1
−
u
+
E
Q̄
−
Q
+
C
αi 2 ,
t+1 2
t
t
t+1
2
n
nκ
n
i2
h
2

L
lf
αi
where Clf = 1−α
, CLF = (1−αFt )2 , and Cαi = max 1, 1−α
.
i
i
Proof. Fix an i. Recall the update equation given in Step 5 of Algorithm 1:
n
o
i,n−1
u,n 
n−1
n
i i,n−1
u
ūi,n
=
Π
ū
−
diag(γ
)
ψ
ū
,
Q̄
,
W
.
Xt
t
t
t
t
t
t+1
t+1

Using the non-expansive property of the projection operator, we have
n
o
 i,∗
i,∗ 2
i,n−1
u,n 
n−1
n
i i,n−1
u ū
u u
ūi,n
=
Π
−
diag(γ
)
ψ
ū
,
Q̄
,
W
−
Π
X
X
t
t
t
t
t
t − ut
t+1
t+1
t
t
2

u,n 2
n
i i,n−1
≤ ūi,n−1
− ui,∗
, Q̄n−1
t
t − diag(γt ) ψt ūt
t+1 , Wt+1 2
γt2
n2
i,n−1
i,∗ T
u,n 
n−1
− 2 (ūt
− ut ) diag(γtn ) ψti ūi,n−1
, Q̄t+1
, Wt+1
.
t

≤ ūi,n−1
− ui,∗
t
t

2
2

+ C αi

Recall from the proof of Lemma 1 that we can write
 i i,n−1 ∗

u,n 
i,n
ψti ūti,n−1 , Q̄n−1
, Qt+1 , Wt+1 + i,n
t+1 , Wt+1 = E ψt ūt
t+1 + ξt+1 ,
16

2
2

(5.7)


T
n
i i,n−1 , Q̄n−1 , W u,n in
from which we see that the cross-term (ūi,n−1
− ui,∗
t
t ) diag(γt ) ψt ūt
t+1
t+1
the chain of inequalities (5.7) can be expanded into three terms. These terms can now
be analyzed separately.
By (4.4) and Assumption 2(iv), we see that the derivative of

E ψti u, Q∗t+1 , Wt+1 in u(s, a) can be expressed as



lf
∂
ft (u(s, a)|s, a)
E ψti u, Q∗t+1 , Wt+1 (s, a) =
≥
.
(5.8)
∂u(s, a)
1 − αi
1 − αi


n−1
Since every state is visited with probability
larger than ε,
we know that E γtn (s, a) Gt+1
≥


εγt
i u∗ , Q∗ , W
.
Combining
this
with
(5.8)
and
E
ψ
=
0,
it
follows
that
(almost
surely)
t+1
t
t
t+1
n
h

 n−1 i
i,∗ T
n
i i,n−1
∗
Gt+1
E (ūi,n−1
−
u
)
diag(γ
)
E
ψ
ū
,
Q
,
W
t+1
t
t
t+1
t
t
t
h
T  i i,n−1 ∗

 i ∗ ∗
i
εγt i,n−1
≥
ūt
− ui,∗
E
ψ
ū
,
Q
,
W
−
E
ψ
u
,
Q
,
W
(5.9)
t+1
t+1
t
t+1
t
t
t+1
t
t
n
εγt lf
2
.
ūi,n−1 − ui,∗
≥
t
2
n 1 − αi t

u,n
Note that γtn (s, a) depends on (Stn , ant ) and i,n
t+1 (s, a) depends on Wt+1 . By independence
u,n
and the unbiased property of i,n
of (Stn , ant ) and Wt+1
t+1 (s, a) of (5.1),
i
h
i,∗ T
n−1
n i,n
(5.10)
G
E (ūi,n−1
−
u
)
diag(γ
)

t
t
t
t+1 = 0.
t+1

Moving on to the third term, using (5.3), the fact that γtn has exactly one nonzero component, and the monotonicity of lp norms, we can deduce
i γ
h
LF
t
i,∗ T
n−1
n i,n
∗
G
E −(ūi,n−1
−
u
)
diag(γ
)
ξ
ūi,n−1
− ui,∗
Q̄n−1
t
t
t
t
t
t+1 ≤
t+1
t+1 − Qt+1 ∞
∞
n
1 − αi
γt LF
n−1
Q̄t+1
− Q∗t+1 2 .
≤
ūi,n−1 − ui,∗
t
2
n 1 − αi t
Using the inequality 2ab ≤ a2 κ + b2 /κ for κ > 0 on the above, we arrive at
i
h
n−1
T
n i,n
)
diag(γ
)
ξ
G
E −2 (ūti,n−1 − ui,∗
t
t
t+1
t+1
≤

γt L2F κ
ūi,n−1 − ui,∗
t
n (1 − αi )2 t

2
2

+

γt
n−1
− Q∗t+1
Q̄t+1
nκ

2
.
2

(5.11)

Finally, the statement of the lemma follows by taking expectations of the inequalities (5.7),
(5.9), (5.10), and (5.11) and combining.
Similarly, the next lemma relates the error of the approximate value function Q̄nt to the
error in the last iteration (i.e., of Q̄n−1
), the error of the value function in the next time
t
period (i.e., of Q̄n−1
),
and
the
error
of
all of the approximate quantiles ūi,n
t . Because the
t+1
analysis is similar to that of Lemma 2, we relegate the proof to Appendix A.
Lemma 3. Under the same conditions as Lemma 2, the sequence of approximations Q̄nt
generated by Algorithm 1 satisfies, for any κ0 , κ1 , . . . , κm > 0,
m
 h
h
i 
i
X
ηt 
2
n
∗ 2
E Q̄t − Qt 2 ≤ 1 −
2 ε − LΦ κ0 − LΦ
κi E Q̄tn−1 − Q∗t 2
n
i=1

ηt LΦ h n−1
+
E Q̄t+1 − Q∗t+1
n κ0

2
2

17

i

+

m
ηt X LΦ h i,n−1
E ūt
− ui,∗
t
n
κi
i=1

2
2

i

+ CH

ηt2
,
n2

where CH bounds the term
i2
h

1,n−1
m,n−1
n−1
(s,
a)
−
H
ū
,
.
.
.
,
ū
,
Q̄
,
W
(s,
a)
≤ CH ,
E Q̄n−1
t
t+1
t
t
t
t+1
for all (s, a) ∈ U.

Proof. See Appendix A.
With these preliminary results in mind, we move on to the theorem that states our
O(1/n) convergence rate and provide a sketch of the proof.
Theorem 3 (Convergence Rate). Choose initial approximations Q̄0 ∈ RD and ūi,0 ∈ RD
u
for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt for all t ∈ T . Under Assumptions
1(vi)–3, the ε-greedy sampling policy of (5.5), and the deterministic harmonic stepsizes given
in (5.6), the sequences of iterates ūi,n for i ∈ {1, . . . , m} and Q̄n generated by Algorithm

2
1 satisfy convergence rates of the form E ūi,n − ui,∗ 2 ≤ O (1/n) for i ∈ {1, . . . , m} and

2
E Q̄n − Q∗ 2 ≤ O (1/n).
Sketch of Proof. The proof is also by backwards induction on t, where the induction hypothesis at each step is that Q̄nt+1 converges at a rate of O(1/n). Applying Lemma 2 along
with some additional analysis, it is possible to show that ūi,n
also converges at a rate of
t
O(1/n). Lemma 3 then completes the proof by showing the desired rate for Q̄nt . The details
are given in Appendix A.

Although both ūi,n and Q̄n converge at a rate of O (1/n), the slower sequence is Q̄n .
To see why this step is slower, one can compare Lemmas 2 and 3. From Lemma 2, we
see that the error of the quantity ūi,n
t depends only on the “last iteration error” (the error
i,n−1
n−1
. In contrast, Lemma 3
of ūt
) and the error of the next stage value function Q̄t+1
n
shows that the error of Q̄t also depends on the error of the current stage quantile ūi,n−1 ,
in addition to its own “last iteration error” and the error of the next stage value function.
The interpretation of these bounds is that accuracy of ūi,n−1 is needed before we obtain
an accurate approximation of Q̄nt , exactly as our intuition would suggest, given Figure 1.
Indeed, this is observed in empirical experiments (see Section 7), motivating the second
contribution of this paper, a procedure aimed toward speeding up the slow step of the ADP
algorithm.

6

Efficient Sampling of the “Risky” Region

Arguably, the biggest practical issue with a Monte Carlo-based algorithm in the setting of
risk-averse decision making is the question of sampling. To illustrate, suppose the one-step
conditional risk measure is CVaR at a level of α = 0.99. Because α is close to 1, the iterates
of the Dynamic-QBRM ADP algorithm are volatile and exhibit poor empirical convergence
rates, as shown in Figure 2. In this section, we discuss a method to address this issue.

6.1

Overview of the Main Idea

As we have mentioned, Dynamic-QBRM ADP can be applied in situations where the distribution of the information process {Wt } is unknown, a common assumption for ADP
algorithms. An example of when such a paradigm can be useful is when one has access to
real data, but no stochastic model from which to simulate. Another example is a black-box
simulator where the user has little to no control of its parameters. In these scenarios, a
18

0

0
Q̄nt
ūnt

-100

Q̄nt
ūnt

-100

-200

-200

-300

-300

-400

-400

-500

-500

-600

-600
-700

-700
0

1000

2000

3000

4000

0

5000

1000

2000

Visits

3000

4000

5000

Visits

Figure 2: Sample Paths of Dynamic-QBRM ADP

good remedy to any sampling issue is to implement an adaptive stepsize rule, similar to
the likes of George and Powell [2006], Schaul et al. [2013], and Ryzhov et al. [2015] (the
third reference develops a stepsize rule specifically in the context of ADP), that can adjust
based on previously observed data points. However, if the distribution of the stochastic
process {Wt }Tt=0 is known, then we can propose a new companion procedure to control the
sampling process as our ADP algorithm progresses. The procedure takes advantage of the
idea of importance sampling (see, e.g., Bucklew [2004]) and is inspired by adaptive sampling
approaches like the cross-entropy method of Rubinstein [1999].

ADP

ūi,n , Q̄n
approximations

RDS

sampling dist.

ADP

ūi,n+1 , Q̄n+1
approximations

...

RDS
sampling dist.

Figure 3: RDS Algorithm Idea

Recall from the results of the previous section that the convergence of Q̄nt , the approximation of the value function, is expected to be slower than that of ūi,n
t , the auxiliary
variable. In this section, we propose a procedure called risk-directed sampling (RDS) to
improve the sampling efficiency for the step of Dynamic-QBRM ADP where Q̄nt is updated,
i.e., Step 6 of Algorithm 1. The main idea, as illustrated in Figure 3, is to use the ADP
approximations to drive the learning of the sampling distribution and simultaneously use
the sampling distribution to generate the samples for updating the ADP approximations.
As the figure shows, the hope is that our procedure updates the distribution in such a way
that we dedicate samples to the regions of high risk (from where we may otherwise not see
many samples).

6.2

Risk-Directed Sampling

Suppose the distribution of the exogenous information Wt+1 has a density pt (w). Notice
that by (4.3), we have for any (s, a) and t,
Z

m,∗
∗
∗
Qt (s, a) = Ht ut1,∗ , u2,∗
(6.1)
t , . . . , ut , Qt+1 , w (s, a) pt (w) dw.
For convenience, we use the shorthand notation


m,∗
∗
Ht∗ (w|s, a) = Ht ut1,∗ , u2,∗
t , . . . , ut , Qt+1 , w (s, a),

to emphasize the variable of integration, w. From the principle of importance sampling
(see, e.g., Bucklew [2004]), it is known that to produce a low-variance estimate of Q∗t (s, a)
19

using Monte Carlo sampling, one should sample from a distribution whose density is nearly
proportional to the absolute value of the integrand of (6.1).
Our approach takes advantage of this proportionality property of the optimal density
and directly constructs an approximation to the absolute value of the integrand of (6.1) by
minimizing a mean squared error. The idea is to capture the “risky regions,” i.e., the parts
of the outcome space where the integrand is large. When the parametric class of sampling
distributions is chosen to be a mixture class (as we do), this introduces a simplification by
allowing us to effectively remove the constraint kθk1 = 1 from the optimization problem.
In addition, we prove that our method converges without exact knowledge of the function
Ht∗ (w|s, a). Instead, the algorithm works in conjunction with Dynamic-QBRM ADP by
using approximations defined by

n−1
Htn (w|s, a) = Ht ū1,n−1
, ū2,n−1
, . . . , ūm,n−1
, Q̄t+1
, w (s, a),
t
t
t

n−1
where ū1,n−1
, ū2,n−1
, . . . , ūm,n−1
, Q̄t+1
are outputs from Dynamic-QBRM ADP. Like the
t
t
t
main ADP method, the procedure is fully adaptive and updates to the sampling distribution
are made at every iteration. To our knowledge, an adaptive importance sampling approach
for constructing risk-averse policies has not been considered in the literature. However,
Kozmı́k and Morton [2014] employs importance sampling from a different perspective: for
the evaluation of risk-averse policies in stochastic programming.
Let {φkt }K
k=1 be the set of densities for K prespecified basis distributions (we also refer
to these as basis functions) for time t, from which we create a mixture distribution used
for sampling. Note that it is often sufficient to specify one set of distributions for all t, as
we do in Section 7. Our goal
Pis to take the traditional regression approach and develop an
approximation of the form k θ̄tk,n (s, a) φkt (w) ≈ Ht∗ (w|s, a) pt (w), where θ̄tk,n (s, a) is the
weight of the k-th basis density at iteration n and state (s, a). For convenience, we also
use the notation θ̄tn (s, a) to be a vector in RK with components θ̄tk,n (s, a) for k = 1 to K.
Choose a (large) compact subset W ⊆ W and define a random variable W u that is uniformly
distributed over W, i.e., assume that its density is given by the function pu (w) = CW 1{w∈W}
for some constant CW . In what follows, we measure the mean squared error using this
distribution, though of course, other reasonable choices may exist and the algorithm can
be adapted accordingly. For a function F : W → R and φ = (φ1 , φ2 , . . . , φK )T , define the
projection operator
h
2 i
Πφ F = arg min E θT φ(W u ) − F (W u ) ,
(6.2)
θ≥0

which maps F to a nonnegative coefficient vector θ corresponding to the best, i.e., minimum
mean squared error, approximation under φ.


Proposition 2. Under the condition that E φ(W u ) φ(W u )T is positive definite, the optimization problem of (6.2) has a unique solution.

Proof. If the positive definiteness condition holds, then the problem is a strictly convex
quadratic program (see, e.g., Boyd and Vandenberghe [2004]).
Exactly in the spirit of importance sampling, we define the optimal sampling density
(within the mixture class) to be the probability density function corresponding to the best fit
of the integrand using the given basis functions, i.e., [Πφ [ |Ht∗ ( · |s, a)| pt (·)]]T φ(w), which is
computed by normalizing the weight vector so that it has unit norm. The sampling density
we use at time t, iteration n, and state (s, a) and the corresponding likelihood ratio are
20

given by

T
p̄n−1
(w|s, a) ∝ θ̄tn−1 (s, a) φ(w)
t

−1
q,n  n−1
q,n
and Lnt+1 (s, a) = pt Wt+1
p̄t
Wt+1
s, a ,

q,n
where Lnt+1 ∈ Rd is defined for each t and n and Wt+1
∈ W is the random sample genern−1
n−1
ated by from p̄tP (w|s, a). Assume that if θ̄t (s, a) = 0 (the zero vector), then we set
p̄n−1
(w|s, a) ∝ k φk (w). Finally, to update our coefficients θ̄tn (s, a) from one iteration to
t
the next, we need to introduce another Gtn -measurable stepsize sequence {βtn } with βtn ∈ Rd
for each t. The description of the new algorithm is given in Algorithm 2 below.
Besides the additional input parameters, there are three points of departure from Algoq,n
rithm 1. In Step 3, the information process sample Wt+1
is drawn according to p̄n−1
(w|s, a)
t
rather than the true distribution, pt (w). In Step 5, we include an additional factor corresponding to the likelihood ratio Lnt+1 in order to apply importance sampling. Step 7 is new
and represents the updating step for the sampling coefficients θ̄tn (s, a) using a stochastic
approximation step. Before moving on to the convergence result, we need another stepsize
assumption and an assumption on the basis functions. In addition, Lemma 5 verifies a
technical condition that we need for the convergence result of Theorem 5.

Algorithm 2: Dynamic-QBRM ADP with Risk-Directed Sampling
Input: Initial estimates of the value function Q̄0 ∈ RD and quantiles ūi,0 ∈ RD for
i ∈ {1, 2, . . . , m}.
Basis distributions φk and initial sampling coefficients θ̄0 (s, a) ∈ RK .

Stepsize rules γtn , ηtn , and βtn for all n, t.
Output: Approximations {Q̄n } and {ūi,n }.

Set Q̄nT = 0 for all n.
for n = 1, 2, . . . do
2
Choose an initial state (S0n , an0 ).
for t = 0, 1, . . . , T − 1 do
u,n
3
Draw a sample of the information process Wt+1
from the distribution of Wt+1 .
q,n
q,n
Draw an IS sample Wt+1 so that Wt+1 (s, a) ∼ p̄n−1
(w|s, a).
t

1

4

5

6

7

Update auxiliary variable approximations for i = 1, . . . m:
n
o
u,n 
i,n−1
, Q̄n−1
ūi,n
− diag(γtn ) ψti ūi,n−1
.
t
t = ΠXtu ūt
t+1 , Wt+1

Compute an estimate of the future cost based on the current approximations:
h
i
q,n 
, . . . , ūtm,n−1 , Q̄n−1
q̂tn = diag(Lnt+1 ) Ht ū1,n−1
t
t+1 , Wt+1 .

Update approximation of value function:
n
o
Q̄nt = ΠXtq Q̄n−1
− diag(ηtn ) Q̄n−1
− q̂tn .
t
t

q,n
Update the samplingh coefficients for eachh state. Let w = Wt+1
(s, a) and

T
n−1
n−1
n
n
θ̄t (s, a) = θ̄t (s, a)−βt (s, a) θ̄t (s, a) φ(w)
i
i+
− Htn (w|s, a) pt (w) φ(w) pu (w) p̄n−1
(w|s, a)−1 ,
t

where [ · ]+ is taken componentwise.

8

n
Choose next state (St+1
, ant+1 ).
end
end

21

Assumption 4. For all (s, a) ∈ U and t ∈ T , suppose βtn is Gtn -measurable and

(i) βtn (s, a) = β̃tn−1 1{(s,a)=(Stn ,ant )} , for some β̃tn−1 ∈ R that is Gtn−1 -measurable.

(ii)

∞
X

n=1

βtn (s, a)

= ∞,

∞
X

n=1

βtn (s, a)2 < ∞

a.s.



(iii) E (β̃tn−1 )2 ≤ O(n−1− ) for some  > 0.

Assumption 5. With regard to the basis distributions φk , k ∈ {1, 2, . . . , K}, the following
hold:


(i) E φ(W u ) φ(W u )T is a positive definite matrix (so that we may apply Proposition 2),
(ii) Let q̂tn ∈ Rd be defined as in Step 5 of Algorithm 2. There exists a constant Cφ such
that the following
for any t and (s, a): given θ ≥ 0 and kθk1 = 1, P
the tails of the
P kholds
k are “heavy” enough to guarantee that if W q,n ∼
k k
distribution
θ
φ
t+1
k θ φ , then
 n
 k
2
E q̂t (s, a) ≤ Cφ .

We remark that due to the compactness of the sets Xtu and Xtq , many of the terms
in the definition of q̂tn are bounded; therefore, the crucial terms that affect
 Assumption
5(ii) are Lnt+1 and the cost function ct . The condition that E βtn (s, a)2 ≤ O(n−1− )
is not particularly difficult to satisfy; for example, our deterministic harmonic stepsizes
satisfy the condition with  = 1. In addition, polynomial rules of the form βtn (s, a) =
n−1/2−/2 1{(s,a)=(Stn ,ant )} work as well. We state the convergence result for Algorithm 2 in
Theorem 5. Since Step 7 of Algorithm 2 does not project to a compact set, we cannot make
use of [Kushner and Yin, 2003, Theorem 2.4] as we did before. Instead, our convergence
result is derived from a theorem of Pflug [1996] for stochastic approximation, which requires
the result of Lemma 5. The rest of the proof is mostly standard and is thus deferred to
Appendix A.
Because the importance sampled stochastic processes are corrected for in expectation
using the likelihood ratio Lnt+1 , we should not expect the rate of convergence to change.
Indeed, we have the following theorem.
Theorem 4 (Convergence Rate). Choose initial approximations Q̄0 ∈ RD and ūi,0 ∈ RD
u
for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt for all t ∈ T . Under Assumptions
1(vi)–5, the ε-greedy sampling policy of (5.5), and the deterministic harmonic stepsizes given
in (5.6), the sequences of iterates ūi,n for i ∈ {1, . . . , m} and Q̄n generated by Algorithm

2
1 satisfy convergence rates of the form E ūi,n − ui,∗ 2 ≤ O (1/n) for i ∈ {1, . . . , m} and

2
E Q̄n − Q∗ 2 ≤ O (1/n).

Proof. The proof is analogous to that of Theorem 3. We need versions of Lemma 2 and
Lemma 3 for Algorithm 2. The difference is that we need to deal with the term Lnt+1 (s, a),
which has expectation equal to 1.

Theorem 5 (Convergence of Risk-Directed Sampling Procedure). Choose initial approxiu
mations Q̄0 ∈ RD and ūi,0 ∈ RD for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt
for all t ∈ T . Under Assumptions 1(vi)–5, the ε-greedy sampling policy of (5.5), and the
deterministic harmonic stepsizes given in (5.6), Algorithm 2 generates a sequence of iterates

22

Q̄n that converge almost surely to the optimal value function Q∗ . Moreover, the sampling
coefficients θ̄tn (s, a) converge to the optimal sampling coefficients under φ:
i
h
a.s.
θ̄tn (s, a) −→ Πφ Ht∗ (·|s, a) pt (·)

for each t and (s, a) ∈ U.
Proof. See Appendix A.

Not surprisingly, the new algorithm retains the same theoretical properties of the standard Dynamic-QBRM ADP without RDS. In addition, we now have a sampling density
that converges to the optimal sampling density as the algorithm progresses (Theorem 5).
What remains for us to explore are the empirical convergence rates of the two approaches.

7

Numerical Results

The recent surge of interest in energy and sustainability has shown that when pertaining to
the question of risk, the problem of optimal control of energy storage assets is an especially
rich application domain. In this section, we illustrate our proposed ADP algorithm by way
of a stylized energy trading and bidding problem in which both heavy tails and extreme
events play a prominent role.

7.1

Model

We consider the problem of using energy storage to trade in the electricity market, i.e.,
energy arbitrage, with the caveat that there is the possibility of some financial penalty
when the amount of stored energy is low. For example, since storage for backup purposes
is rarely in use, one might consider using it to generate a stream of revenue by interacting
with the market. However, this immediately introduces a source of risk: in the rare event
when backup is needed, it is crucial that there is enough energy to cover demand. Failure to
do can cause complications, so it is useful to consider a risk-averse policy to this problem.
See Xi et al. [2014] for a detailed model of the shared storage situation, but solved with
a risk-neutral objective. When no penalty is assessed, we assume that there is a modest
reward. Finally, we introduce a bidding aspect to the problem where one must place bids
prior to the desired transaction time (see Jiang and Powell [2015b] for detailed model of
hour-ahead bidding).
For t = 0 to t = T , let St ∈ S = {0, 1, . . . , Smax } be the amount of energy in storage
and let Pt ∼ log N (µP (t), σP (t)2 ) be the (heavy-tailed) spot price of electricity. Also,
let Ut ∼ N (0, σU2 ) be independent of Pt , and suppose our simple model of storage-based
penalties is as follows. Given two constants 0 < a < b, where a represents the “rate of
reward” and b represents the “rate of penalty,” assume that the reward/penalty assessed at
time t + 1 is given by
h
i
Ft+1 = µS (St ) + Ut+1 b · 1{µS (St )+Ut+1 <0} − a · 1{µS (St )+Ut+1 ≥0} ,

where µS : S → R is a nondecreasing function, signifying that as St increases, the possibility
of penalty decreases (rare event). Our two-dimensional action at each time t is given by
 − +
+
−
+
at = (b−
t , bt ) ∈ A ⊆ (b , b ) : 0 ≤ b ≤ b ≤ bmax ,
where |A| < ∞. There are no constraints, so As = A for all s ∈ S. We call b−
t the buy
+
bid and bt the sell bid : if Pt fluctuates below the buy bid, we are obligated to buy from
23

the market and if Pt rises above the sell bid, we are obligated to sell to the market. In
addition, we are penalized the amount of the spot price if we are to sell but the storage
device is empty, i.e., St = 0. In this problem, our information process Wt is given by the pair
(Pt , Ut ) and is independent of the past. We find it most natural to model this problem in
the sense of maximizing revenues (or “contributions”) rather than minimizing costs; hence,
the contribution function is
h
i
ct (St , at , Wt+1 ) = −Ft+1 + Pt+1 1{b+ <Pt+1 } − 1{b− >Pt+1 } − 1{St =0} 1{b+ <Pt+1 } ,
t

t

t

and the transition function is given by St+1 = [min{St + 1{b− >Pt+1 } − 1{b+ <Pt+1 } , Smax }]+ .
t
t
The objective is to optimize the risk-averse model given in (3.2) with the “min” operator
replaced by a “max” and taking the QBRM to be mean-CVaR at α = 0.99 (this refers to
the lower 0.01 tail of the distribution since we are now in the setting of rewards, not costs).
Various values of λ are considered.

7.2

Parameter Choices

We let T = 12 and the size of the storage device
 be Smax = 6. Next, suppose the values
of µS are chosen so that P µS (St ) + Ut+1 < 0 are 0.1, 0.05, 0.02, 0.01, 0.01, 0.001, and
0.001 for St = 0, 1, . . . , Smax (representing increasingly rare penalties) with σU2 = 1. The
penalty parameters are set to be a = 5 and b = 500, making the penalty events relatively
severe. We denote the (seasonal) mean of the spot price process by m(t) = E(Pt ) =
50 sin(4πt/T ) + 100, with a constant variance v = Var(Pt ) = 3000. By the properties of the
lognormal distribution,
this leads to the parameters µP (t) = log(m(t) (1 + v/m(t)2 )−1/2 )
p
and σP (t) = log(1 + v/m(t)2 ). The maximum bid is chosen to be bmax = 500 and each
dimension (the range from 0 to 500) is discretized into increments of 50, resulting in an
MDP with approximately 7,000 states.
Finally, we discuss our choice of parameters for the risk-directed sampling (RDS) procedure. For φ1 , we simply take the true distribution of Wt+1 and for the remaining φk , we
select a set of bivariate normal distributions placed in a grid, with mean of the first component taking values in {50, 175, 300} and the mean of the second component taking values
in {−3, −1, 1}. The standard deviations are chosen to be 750 and 0.25 for the respective
components. This is a fairly general choice of basis functions that uses very little problem
specific information and/or tuning. Figure 4B shows an example of the shape of the sam+
pling density after 500 iterations of RDS for a fixed state of (St , b−
t , bt ) = (0, 150, 300) at
t = T − 1 with λ = 0.5 (note: the z-axis is not normalized here and the range of the Pt axis
has been decreased to focus on the nonzero areas). As we would expect, the algorithm has
chosen to allocate a relatively large sampling effort toward small values of Ut , presumably
due to the fact that these events are very costly with b = 500 and thus critical to the
estimation of CVaR. The last few observations are shown as red points.

7.3

Evaluating Policies

Here, we describe the procedure for evaluating policies under the objective (3.2). Let
π
{Aπ0 , Aπ1 , . . . , AπT −1 } be a policy; to compute (3.2), let St+1
= S M (s, Aπt (s), Wt+1 ) and we
solve the recursive equations given by

π
π
Vtπ (s) = ραt ct (s, Aπt (s), Wt+1 ) + Vt+1
(St+1
) for all s ∈ S, t ∈ T ,
(7.1)
VTπ (s) = 0 for all s ∈ S.

24

(A) Basis distributions φk (equally weighted)

(B) Shape of sampling density, n = 500

Figure 4: Example Illustration of Risk-Directed Sampling (λ = 0.5)

The value of the policy is then given by V0π (S0 ), but since Wt is continuously distributed,
we cannot, in general, solve these equations exactly. For this reason, we select a large, finite
sample Ω̂ ⊆ Ω, and apply the standard sample average approximation (SAA) technique
of Kleywegt et al. [2002] to the Bellman recursion (7.1) along with the standard linear
programming method suggested in Rockafellar and Uryasev [2000] for computing conditional
value at risk. Next, in order to have an optimality benchmark against which ADP policies
can be compared, we also apply the SAA technique to the Bellman recursion for the optimal
policy specified in Theorem 1 (with ρt = ραt and the “min” replaced by “max”).
We use 50,000 realizations of Wt = (Pt , Ut ) for each t = 1, 2, . . . , T in our simulations
below. For the purposes of our numerical work and Figure 7, we refer to the SAA optimal
∗ (S )) as our benchmark for “100% optimal” and a
policy (which has value function V0,β
0
corresponding SAA myopic policy, denoted by πm , (i.e., the policy which takes the action
to maximize ραt (ct (s, a, Wt+1 )) with zero continuation value to be “0% optimal.” More
precisely, we have
% optimality of policy π (with respect to myopic policy πm ) =

V0π (S0 ) − V0πm (S0 )
.
V0∗ (S0 ) − V0πm (S0 )

Recall that larger V0π (S0 ) is better, as we have switched to the maximization setting for the
numerical application.

7.4

Results

Let us first show a few plots that illustrate the effectiveness of the RDS procedure. Once
+
again, we fix the state to be (St , b−
t , bt ) = (0, 150, 300) at t = T −1 and set λ = 0.5, but now
we focus on the evolution of the approximations Q̄nt and ūnt ; representative sample paths
are shown in Figure 5 (the first plot is from Figure 2, repeated for comparison purposes).
To make the comparison as fair as possible, the random number generator seed is set
so that the sequence of observations used to update ūnt is the same in both cases (i.e., the
gray lines are identical). The true limit points are approximately u∗t (St , at ) ≈ −555 and
Q∗t (St , at ) ≈ −387, and we notice that in the case of Figure 5B, a decent approximation
of Q∗t is obtained around visit 2500, while the approximation in Figure 5A does not settle
until the end. This drastic difference in empirical convergence rate can have a significant
impact on the performance of Dynamic-QBRM ADP as the large errors shown in Figure
5A are propagated backwards in time.
In fact, we can illustrate this impact by examining the value function approximations
back at time t = 0, with and without RDS, compared to the SAA optimal value function.
Figure 6 compares the two methods after N = 5,000,000 iterations by varying the bid
25

0

0
Q̄nt
ūnt

-100

Q̄nt
ūnt

-100

-200

-200

-300

-300

-400

-400

-500

-500

-600

-600

-700

-700
0

1000

2000

3000

4000

5000

0

1000

2000

Visits

3000

4000

5000

Visits

(A) Without RDS

(B) With RDS

Figure 5: Sample Paths of Approximations Generated by Dynamic-QBRM ADP (λ = 0.5)

dimensions of a fixed state S0 = 0. Notice that in Figure 6B, the approximation from the
RDS procedure very closely resembles the optimal value function of Figure 6C even after
errors are propagated through T = 12 steps. The same cannot be said of the approximation
in Figure 6A.
-800

-1200
-1600
500

500

250

b−
0

250
0

0

b+
0

(A) Without RDS

-800

Q∗0 (0, · )

Q̄N
0 (0, · )

Q̄N
0 (0, · )

-800

-1200
-1600
500

500

250

b−
0

250
0

0

b+
0

(B) With RDS

-1200
-1600
500

500

250

b−
0

250
0

0

b+
0

(C) SAA Optimal

Figure 6: Surface Plots of Value Function Approximations at t = 0 (λ = 0.5)

Lastly, we illustrate that in addition to the improved empirical convergence rates observed in Figures 5 and 6, running Dynamic-QBRM ADP with RDS has a noticeable effect
on the resulting policies as well. Using the backward recursive SAA evaluation procedure
described above, we plot optimality percentages of risk-averse policies produced by the
two variants of the algorithm for λ ∈ {0.6, 0.55, 0.5, 0.45, 0.4} in Figure 7. We see that
RDS has the advantage in all cases, especially during the early iterations. Figure 7 uses
a log-scale to display the early progress of the algorithms every 50,000 iterations from
N = 0 to N = 1,000,000 and also the asymptotic progress every 1,000,000 iterations from
N = 1,000,000 to N = 5,000,000. Our simulations also suggest that as λ decreases, the
advantage of using RDS diminishes; this is not surprising because decreasing λ corresponds
to deemphasizing the “risky” events that RDS seeks to sample. In fact, we find that in this
particular example problem, the algorithms are practically indistinguishable for λ ≤ 0.3.
Computation times for the three methods are given in Table 1; notice that we are able to
run the ADP methods for 5 million iterations in roughly half of the time it takes to compute the SAA optimal solution (the large number of scenarios presents a challenge for the
optimization routine implemented in CPLEX). The online companion contains additional
numerical results of evaluating our policies on more “practical” metrics of risk and reward,
rather than the dynamic risk measure objective function. We observe that the policies
generated by our algorithms behave in an intuitively appealing way.

26

100

90

90

80

80

80

70
60
50
RDS
No RDS

40
30

105
106
Number of Iterations (Log Scale)

% Optimality

100

90
% Optimality

% Optimality

100

70
60
50
RDS
No RDS

40
30

100

90

90

80

80

% Optimality

% Optimality

50
RDS
No RDS

40

105
106
Number of Iterations (Log Scale)

(B) λ = 0.55

100

70
60
50
RDS
No RDS

40
30

60

30

105
106
Number of Iterations (Log Scale)

(A) λ = 0.6

70

(C) λ = 0.5

70
60
50
RDS
No RDS

40
30

105
106
Number of Iterations (Log Scale)

(D) λ = 0.45

105
106
Number of Iterations (Log Scale)

(E) λ = 0.4

Figure 7: Comparison of Dynamic-QBRM ADP with and without RDS
Method

λ = 0.6 λ = 0.55 λ = 0.5 λ = 0.45 λ = 0.4

Dynamic-QBRM ADP
Dynamic-QBRM ADP with RDS
SAA Backward Recursion

810.24 819.40
926.25 983.33
2119.74 2058.36

816.40 830.69
989.87 926.85
2163.91 2118.84

811.45
973.50
2072.79

Table 1: CPU Times for N = 5, 000, 000 Iterations of ADP and SAA Optimal Solution (minutes)

8

Conclusion

In this paper, we propose a new ADP algorithm for solving risk-averse MDPs [Ruszczyński,
2010] under a class of quantile-based risk measures. The algorithm uses a two-step updating
procedure to approximate the optimal value function and we prove that it converges almost
surely with a rate of O(1/n). We also show a companion sampling procedure to more
efficiently sample the “risky” region of the outcome space and then prove that the sampling
distribution converges to one that is, in a sense, the best within a parametric class. Next,
we describe an energy storage and bidding application, on which we demonstrate that the
RDS sampling approach provides significant benefits in terms of the empirical convergence
rate. Moreover, we illustrate that the approximate policies produced by Dynamic-QBRM
ADP behave in an intuitive and appealing way in terms of the tradeoff between risk and
reward (see online companion), suggesting that it can be readily applied to a variety of
problem settings.

27

A

Proofs

Theorem 2 (Almost Sure Convergence). Choose initial approximations Q̄0 ∈ RD and
u
ūi,0 ∈ RD for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt for all t ∈ T . Under
n
Assumptions 1–2, Algorithm 1 generates a sequence of iterates Q̄ that converge almost
surely to the optimal value function Q∗ .
Proof. We induct backwards on t, starting with the base case t = T . Since QnT (s, a) =
Q∗T (s, a) = 0 for all n and all (s, a), convergence trivially holds. The induction hypothesis is
that Q̄nt+1 (s, a) → Q∗t+1 (s, a) almost surely, and we aim to show the same statement with t
replacing t + 1. As we did before, we analyze the statistical error and approximation error,
q,n
d
q,n
t+1 , ξt+1 ∈ R , except they are now with respect to the stochastic Bellman operator Ht :
1,∗
m,∗
q,n 
∗
∗
q,n
t+1 = Qt − Ht ut , . . . , ut , Qt+1 , Wt+1 ,
q,n
m,∗
q,n 
1,n−1
q,n 
n−1
∗
ξt+1
= Ht u1,∗
, . . . , ūm,n−1
, Q̄t+1
, Wt+1
.
t , . . . , ut , Qt+1 , Wt+1 − Ht ūt
t

Therefore, the update of Step 6 of Algorithm 1 is equivalent to
io
n
h
q,n
q,n
n−1
n
∗
.
(A.1)
+
ξ
Q̄nt = ΠXtq Q̄n−1
−
diag(η
)
Q̄
−
Q
+

t
t
t
t
t+1
t+1


n−1
It is once again clear that E q,n
t+1 (s, a) Gt+1 = 0 almost surely. Furthermore, we can easily
argue that for some positive constant CH ,
i2
h
q,n 
m,n−1
1,n−1
n−1
,
Q̄
,
W
(s,
a)
≤ CH ,
(A.2)
,
.
.
.
,
ū
E Q̄n−1
(s,
a)
−
H
ū
t t
t
t
t+1
t+1
using Assumption 2(i), Assumption 2(ii), the fact that ραt (0) = Φ(0, . . . , 0) = 0, and the
n
boundedness of the iterates Q̄n−1
t+1 and ūt .
q,n 
n
M
Now, fix an (s, a) and let St+1 = S s, a, Wt+1
. Expanding and using Assumption 2(i)
and the property (5.2) we see that
q,n
|ξt+1
(s, a)| ≤ LΦ

n−1
n
n
(St+1
, a0 )
min Q∗t+1 (St+1
, a0 ) − 0 min
Q̄t+1
n

a0 ∈AS n

+ LΦ
≤ LΦ

a ∈AS

t+1

h

X

ūi,n−1
(s, a)
t

i

∗
Q̄n−1
t+1 − Qt+1

+
∞

−

X
i

t+1

ui,∗
t (s, a)
ūi,n−1
− ui,∗
t
t

i
.
∞

(A.3)

Taking conditional expectation, applying Lemma 1, and using the induction hypothesis
i,∗
n
∗
(which tells us that ūi,n
t → ut and Q̄t+1 → Qt+1 almost surely), we conclude
 q,n

n−1
E ξt+1
(s, a) Gt+1
→ 0 a.s.,

satisfying one of the conditions of [Kushner and Yin, 2003, Theorem 2.4]. In addition, note
the bounded second moment condition of (A.2), the conditional unbiasedness of q,n
t+1 (s, a),
the stepsize, sampling, and truncation properties of Assumption 1. We once again have the
ingredients to apply the stochastic approximation convergence theorem [Kushner and Yin,
2003, Theorem 2.4] (the objective function for applying the theorem is q 7→ kq − Q∗t k22 ) to
the update equation of (A.1) in order to conclude
Q̄nt (s, a) → Q∗t (s, a)

a.s.,

for every (s, a) ∈ U, completing both the inductive step and the proof.

28

Proposition 1. Under the ε-greedy sampling policy given in (5.5) and the deterministic
harmonic stepsizes given in (5.6), Assumption 1(iii)–(iv) is satisfied.
Proof. Let us consider, for a fixed t and (s, a), the stepsize sequence
γt
n n .
γtn (s, a) =
1
n {(s,a)=(St ,at )}
The argument for the case of ηtn (s,P
a) is, of course, exactly the same and is omitted. First,
∞
n
2
surely, is trivial. We thus
notice that the second part, i.e.,
n=1 γt (s, a) < ∞ almost P
n
focus on proving that the second part of the assumption, that ∞
n=1 γt (s, a) = ∞ almost
surely, holds in this case. An iteration n is called a visit to state (s, a) if (s, a) = (Stn , ant )
(so the stepsize is nonzero). Let {∆kt (s, a)}k≥1 be a process that describes the interarrival
times for visits to the state (s, a); in other words, ∆kt (s, a) is the number of iterations that
pass between the (k − 1)-st and k-th visits to the state (s, a). The arrival time sequence
(the iterations for which (s, a) is visited) is thus given by the sum
Ntk (s, a)

k
X

=

0

∆kt (s, a) for k > 0.

k0 =1

To simplify notation, we henceforth drop the dependence of these processes on (s, a) and use
∆kt = ∆kt (s, a) and Ntk = Ntk (s, a). By selecting only the iterations for which the stepsizes
are nonzero, we see that
∞
X

Nk
γt t

∞
X

=

n=1

k=1

Analogous to

Ntk

γtn (s, a).

process, define the deterministic sequence
& k
'
X 0
−2
nkt =
δtk
and δtk =
log k.
log (1 − ε)
k0 =1

Observe that under the ε-greedy sampling policy, the event (s, a) 6= (Stn , ant ) occurs with
probability at most 1 − ε. Hence, by independence, we can show that

k
1
P ∆kt ≥ δtk ≤ (1 − ε)dδt e ≤ 2 .
k

By the Borel-Cantelli Lemma, with probability 1, the events ∆kt ≥ δtk occur finitely
often. Hence, there exists an almost surely finite Kt such that for all k ≥ Kt , it is true that
∆kt ≤ δtk . Clearly, nkt ≤ O(k log k) and
∞
X

nk

γt t =

k=1

PKt −1

∞
X
γt
=∞
k
n
t
k=1

(A.4)

P t −1 i
holds. Let C1 = i=1
and C2 = K
i=1 ∆t be almost surely finite random variables.
Some simple manipulations yield (the following chain of inequalities may be analyzed ω-wise
to obtain the a.s. qualification)
∞
X
k=1

Nk
γt t

γt
Nti

∞
∞
X
X
γt
γt
=
= C1 +
Pk
k
N
C + i=Kt ∆it
k=1 t
k=Kt 2

≥ C1 +

∞
X

k=Kt

C2 +

29

γt
Pk

i
i=Kt δt

≥ C1 +

∞
X

k=Kt

C2 −

where the final equality follows from (A.4).

γt
PKt −1
i=0

δti + nkt

= ∞ a.s.,

Lemma 3. Under the same conditions as Lemma 2, the sequence of approximations Q̄nt
generated by Algorithm 1 satisfies, for any κ0 , κ1 , . . . , κm > 0,
m
i
i 
h
 h
X
ηt 
n−1
∗ 2
n
∗ 2
2 ε − LΦ κ0 − LΦ
κi E Q̄t − Qt 2
E Q̄t − Qt 2 ≤ 1 −
n
i=1

ηt LΦ h n−1
+
E Q̄t+1 − Q∗t+1
n κ0

2
2

i

m
ηt X LΦ h i,n−1
+
− ui,∗
E ūt
t
n
κi

2
2

i=1

i

+ CH

ηt2
,
n2

where CH bounds the term
h
i2

≤ CH ,
E Q̄n−1
(s, a) − Ht ū1,n−1
, . . . , ūm,n−1
, Q̄n−1
t
t
t
t+1 , Wt+1 (s, a)
for all (s, a) ∈ U.

Proof. Note that although we omit its algebraic form, the existence of Cg is guaranteed
by Assumption 2(ii) and the boundedness of Q̄n−1
(s, a). Let us recall the update equation
t
given in Step 7 of Algorithm 1 can be rewritten as
io
n
h
q,n
q,n
n−1
n
∗
.
+
ξ
−
diag(η
)
Q̄
−
Q
+

Q̄nt = ΠXtq Q̄n−1
t
t
t
t
t+1
t+1

Expanding, we have
Q̄nt − Q∗t

2
2

− Q∗t
≤ Q̄n−1
t


q,n  2
− Q∗t + q,n
+ diag(ηtn ) Q̄n−1
t
t+1 + ξt+1 2
(A.5)
T

q,n 
− 2 Q̄tn−1 − Q∗t diag(ηtn ) Q̄tn−1 − Q∗t + q,n
t+1 + ξt+1 .
2
2

We focus on the cross term. First, by the ε-greedy sampling policy, notice that
h
 n−1 i εηt

2
n−1
n−1
∗
∗ T
n
E Q̄t − Qt diag(ηt ) Q̄t − Qt Gt+1 ≥
Q̄tn−1 − Q∗t 2 ,
n
and by the definition of q,n
t+1 ,
h
i

n−1
∗ T
n q,n
E Q̄n−1
−
Q
diag(η
)

G
t
t
t
t+1
t+1 = 0.

(A.6)

(A.7)

Using the bound (A.3), the fact that ηtn contains exactly one nonzero component, and the
monotonicity of the lp norms, we can see that
h
i
h
ηt
n−1
∗ T
n q,n
∗
E − (Q̄n−1
L
Q̄n−1
− Q∗t ∞ Q̄n−1
−
Q
)
diag(η
)
ξ
G
≤
Φ
t
t
t
t
t+1 − Qt+1 ∞
t+1
t+1
n
i
X i,n−1
i,∗
∗
+ Q̄n−1
−
Q
ū
−
u
t
t
t
t
∞
∞
h
ηt
− Q∗t
≤ LΦ Q̄n−1
t
n
+

30

Q̄tn−1

−

i

2

Q∗t 2

n−1
Q̄t+1

X
i

− Q∗t+1
ūi,n−1
t

2

− ui,∗
t

2

i

.

Again applying 2ab ≤ a2 κ + b2 /κ, we see that for any constants κk > 0, k = 0, 1, . . . , m,

m
h

i η
X
2
t
n−1
n−1
∗ T
n φ,n
E −2 (Q̄t − Qt ) diag(ηt ) ξt+1 Gt+1 ≤ LΦ κ0 +
κi Q̄n−1
− Q∗t 2
t
n
i=1


m
X
1
i,n−1
i,∗ 2
ū
− ut 2 .
+
κi t
i=1
(A.8)
The statement of the lemma follows by taking expectations of (A.5), (A.6), (A.7), and (A.8),
using the CH bound, and combining.
1
∗
Q̄n−1
+
t+1 − Qt+1
κ0

2
2

Theorem 3 (Convergence Rate). Choose initial approximations Q̄0 ∈ RD and ūi,0 ∈ RD
u
for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt for all t ∈ T . Under Assumptions
1(vi)–3, the ε-greedy sampling policy of (5.5), and the deterministic harmonic stepsizes given
in (5.6), the sequences of iterates ūi,n for i ∈ {1, . . . , m} and Q̄n generated by Algorithm

2
1 satisfy convergence rates of the form E ūi,n − ui,∗ 2 ≤ O (1/n) for i ∈ {1, . . . , m} and

2
E Q̄n − Q∗ 2 ≤ O (1/n).
Proof. First, we state a special case of Chung’s Lemma (see Chung [1954]), which is useful
for analyzing a specific type of sequence that arises often in recursive optimization algorithms (see, e.g., Juditsky et al. [2009], Rakhlin et al. [2012]).

Lemma 4 (Chung [1954]). Consider a sequence {an }. Suppose that the following recursive
inequality holds for some b > 1 and every n ≥ 1:


b
c
an ≤ 1 −
an−1 + 2 .
n
n
n
o
c
Then, if k ≥ max b−1
, a0 , it follows that an ≤ nk for n ≥ 1.

To prove the theorem, we
an induction
argument for Q̄n . Let us first consider t = T
 use
2
n
∗
as the base case. Since E Q̄T − QT 2 = 0, it is trivially O(1/n), so let us take our

2
induction hypothesis to be E Q̄nt+1 − Q∗t+1 2 ≤ O(1/n) for a particular t < T . Hence,
 n−1
kq
2
q
there exists kt+1
> 0 independent of n such that E Q̄t+1
− Q∗t+1 2 ≤ t+1
n holds for n ≥ 1
(note the n−1 on the left hand side). We first apply Lemma 2 with γt =
so by Lemma 4, we can take
kti = 2 max

(

2 CLF q
4 Cα
i,∗
kt+1 + 2 2i , ūi,0
t − ut
2
2
ε Clf
ε Cl f

2
2

2
ε Clf

and κ =

ε Clf
C LF

,

)



2
2
ki
ki
to satisfy E ūnt − u∗t 2 ≤ 2nt for n ≥ 1. It follows that E ūtn−1 − u∗t 2 ≤ nt for n ≥ 2.
Finally, we apply Lemma 3 with ηt = 2ε , κ0 = 2 Lε Φ , and κi = 2 mεLΦ for all i. Again by
Lemma 4, we can choose
(
)
4 L2Φ q
4 m L2Φ X i 4 CH
q
0
∗ 2
kt = max
k
+
kt + 2 , Q̄t − Qt 2 ,
ε2 t+1
ε2
ε
i
h
i
2
kq
which gives E Q̄nt − Q∗t 2 ≤ nt for n ≥ 2, the O(1/n) rate needed to complete the

31

inductive step. Therefore, we can conclude that


n

E Q̄ −

2
Q∗ 2

=

T
X
t=0


E Q̄nt − Q∗t

2
2

≤ O (1/n) ,


2
as desired. As for the convergence of E ūi,n − ui,∗ 2 , we note that the values kti obtained
through the induction above can be used to deduce the O (1/n) convergence rate.
i,n
q,n
Similar to the sequences ξt+1
and ξt+1
used in the proofs of Lemma 1 and Theorem 2,
h,n
let us introduce another useful process ξt+1 (s, a) ∈ RK for each t and (s, a). As in Step 7
q,n
of Algorithm 2, let w = Wt+1
. We define
h
i
h,n
ξt+1
(s, a) = Ht∗ (w|s, a) − Htn (w|s, a) φ(w) pt (w) pu (w) p̄tn−1 (w|s, a)−1 ,

representing the error from using Htn (w|s, a) to approximate Ht∗ (w|s, a).

Lemma 5. Under the ε-greedy sampling policy of (5.5) and Assumptions 1(vi)–5, it follows
that
∞
i
h
X
h,n
n−1
< ∞ a.s.
(s, a) 2 Gt+1
β̃tn−1 E ξt+1
n=1

for each t and (s, a).

q,n
and P = φ(w) pt (w) pu (w) p̄tn−1 (w|s, a)−1 . We notice
Proof. For simplicity, let w = Wt+1
that
h,n
(s, a)
ξt+1

2

≤ kP k2 Ht∗ (w|s, a) − Htn (w|s, a)

n−1
− Q∗t+1
≤ kP k2 LΦ Q̄t+1

2

+ kP k2 LΦ

X
i

ūi,n−1
− ui,∗
t
t

2

,

where the second inequality follows from the same steps used in deriving (A.3) and then
applying the monotonicity of lp norms. Squaring, applying the inequality 2ab ≤ a2 κ + b2 /κ
for any κ > 0 to the cross terms, and using the fact that each component of P is bounded
(due to pu ), we can write
X
2
2
2
h,n
n−1
− Q∗t+1 2 +
(s, a) 2 ≤ CQ Q̄t+1
Cu,i ūi,n−1
− ui,∗
,
(A.9)
ξt+1
t
t
2
i

for some constants CQ , Cu,1 , . . . , Cu,m ≥ 0. Taking expectation of (A.9) and using the
 h,n

convergence rate result of Theorem 4, we see that E kξt+1
(s, a)k22 ≤ O(1/n). Now, by
Cauchy-Schwarz and Assumption 4(iii),
h
h
i
 h,n
i
h,n
n−1
E β̃tn−1 E kξt+1
(s, a)k2 | Gt+1
= E β̃tn−1 ξt+1
(s, a) 2
q 
  h,n

≤ E (β̃tn−1 )2 E kξt+1
(s, a)k22 ≤ O(n−1−/2 ),
h
 h,n
i
n−1
so it is clear that the terms E β̃tn−1 E kξt+1
(s, a)k2 | Gt+1
are summable. By the monotone
convergence theorem, we conclude
"∞
#
h
i
X
h,n
n−1
E
β̃tn−1 E ξt+1 (s, a) 2 | Gt+1
< ∞,
n=1

from which the statement of the lemma follows (notice that the term within the expectation
must be finite almost surely).
32

Theorem 5 (Convergence of Risk-Directed Sampling Procedure). Choose initial approxiu
mations Q̄0 ∈ RD and ūi,0 ∈ RD for each i ∈ {1, 2, . . . , m} so that Q̄0t ∈ Xtq and ūi,0
t ∈ Xt
for all t ∈ T . Under Assumptions 1(vi)–5, the ε-greedy sampling policy of (5.5), and the
deterministic harmonic stepsizes given in (5.6), Algorithm 2 generates a sequence of iterates
Q̄n that converge almost surely to the optimal value function Q∗ . Moreover, the sampling
coefficients θ̄tn (s, a) converge to the optimal sampling coefficients under φ:
i
h
a.s.
θ̄tn (s, a) −→ Πφ Ht∗ (·|s, a) pt (·)
for each t and (s, a) ∈ U.

Proof. First, notice that the convergence of Q̄n and ūi,n is (mostly) unaffected by the
addition of the sampling procedure. By the principle of importance sampling, the factor of
q,n
diag(Lnt ) corrects, in expectation, for the fact that Wt+1
is sampled from the importance
n−1
distribution p̄t (w|s, a) rather than pt (w). In addition, Assumption 5 implies that the
condition stating that the second moment of the gradient term is finite, which is needed for
[Kushner and Yin, 2003, Theorem 2.4], still holds and therefore the convergence follows.
We now focus on the last part of the theorem and analyze the convergence of the
sampling coefficients. Let
i
h
θt∗ (s, a) = Πφ Ht∗ (·|s, a) pt (·) ,

and we aim to show θ̄tn (s, a) → θt∗ (s, a) almost surely. The proof technique is standard
and uses a supermartingale convergence argument; see, e.g, [Pflug, 1996, Theorem 5.3], but
several aspects need to be adapted for technical reasons in our setting. Throughout this
proof, fix a t and (s, a). To simplify notation, we use the shorthand p̄tn−1 (w) = p̄n−1
(w|s, a),
t
q,n
n (w) = H n (w|s, a), and H ∗ (w) = H ∗ (w|s, a). First, we decompose
w = Wt+1
∼ p̄n−1
(w),
H
t
t
t
t
t
the gradient term. Define a function ht : RK → RK , so that for θ ∈ RK ,
h
i

−1
ht (θ) = E θT φ(w) − |Ht∗ (w)| pt (w) φ(w) pu (w) p̄n−1
.
(w)
t
Note that ht is affine in θ and is the gradient of the (strictly convex) quadratic objective
function in (6.2), meaning that we can find positive constants Ch,1 and Ch,2 such that
ht (θ)

2
2

≤ Ch,1 θ − θt∗ (s, a)

2
2

+ Ch,2 .

(A.10)

Next, the noise term h,n
t+1 (s, a) is given by


T

θ̄tn−1 (s, a) φ(w) − Ht∗ (w) pt (w) φ(w) pu (w) p̄tn−1 (w)−1 − ht θ̄tn−1 (s, a) ,
h,n
t+1 (s, a) =



n−1
n−1
and it follows that E h,n
(s, a) is a detert+1 (s, a) | Gt+1 = 0. Using the fact that ht θ̄t
ministic affine function of θ̄tn−1 (s, a), we observe that h,n
t+1 (s, a) can be written in the form
 n−1

∗
K×K
X θ̄t (s, a) − θt (s, a) + Z where X ∈ R
and Z ∈ RK are random. Since pu (w) = 0
for w outside of a compact set, the entries of both X and Z are bounded. Let k · k2 denote
the operator norm whenever its argument is a matrix and we have

2
2
n−1
h,n
(s, a) − θt∗ (s, a) 2 + kZk2
t+1 (s, a) 2 ≤ kXk2 θ̄t
(A.11)
2
≤ C,1 θ̄tn−1 (s, a) − θt∗ (s, a) 2 + C,2 ,
h,n
for some constants C,1 and C,2 . We also reproduce the definition of ξt+1
(s, a) given previously:
h
i
h,n
ξt+1
(s, a) = Ht∗ (w) − Htn (w) φ(w) pt (w) pu (w) p̄tn−1 (w)−1 .

33

Step 7 of Algorithm 2 can therefore be written
h
h
ii+

h,n
.
θ̄tn (s, a) = θ̄tn−1 (s, a) − βtn (s, a) ht θ̄tn−1 (s, a) + h,n
(s,
a)
+
ξ
(s,
a)
t+1
t+1

For convenience, let us define

Ant = θ̄tn (s, a) − θt∗ (s, a)
Hence,


h,n
and ĥnt = ht θ̄tn−1 (s, a) + h,n
t+1 (s, a) + ξt+1 (s, a).

 n−1
+  ∗
+ 2
θ̄t (s, a) − βtn (s, a) ĥnt
− θt (s, a)
2

n−1 2
n−1 T n
n−1 2
n
≤ kAt k2 − 2 βt (s, a) (At ) ĥt + β̃t
kĥnt k22 .

kAnt k22 ≤

(A.12)

Taking conditional expectation of the cross term, using Assumption 4, noting the inequality
(An−1
)T ht (θ̄tn−1 (s, a)) > 0 (by strict convexity), and using Cauchy-Schwarz, we get
t



n−1
(A.13)
E −2 βtn (s, a) (An−1
)T ĥnt | Gt+1
≤ −2 ε β̃tn−1 (Atn−1 )T ht θ̄tn−1 (s, a)
t
 h,n

n−1
n−1
n−1
+ 2 β̃t kAt k2 E kξt+1 (s, a)k2 Gt+1

Moving onto the third term of (A.12), we have

 h

h,n
n−1
2
2
E kĥnt k22 | Gt+1
≤ kht (θ̄tn−1 (s, a))k22 + E kh,n
t+1 (s, a)k2 + kξt+1 (s, a)k2

h,n
n−1
(s, a))k2 kh,n
+ 2 kh,n
t+1 (s, a)k2
t+1 (s, a)k2 kξt+1 (s, a)k2 + 2 kht (θ̄t
i

h,n
n−1
+ 2 kht (θ̄tn−1 (s, a))k2 kξt+1
(s, a) 2 | Gt+1
i
h 

h,n
n−1
n−1
2
2
2
(s,
a)k
|
G
+
k
h
(
θ̄
(s,
a))k
(s,
a)k
+
kξ
≤ 3 E kh,n
t t
2
2
2
t+1
t+1
t+1
i
h


h,n
n−1
(A.14)
≤ Cĥ,1 kAn−1
k22 + Cĥ,2 + 3 E kξt+1
(s, a)k22 | Gt+1
t

where the second inequality is due to the inequality 2ab ≤ a2 + b2 and the third inequality,
with new constants Cĥ,1 and Cĥ,2 , is due to the bounds (A.10) and (A.11). Combining
(A.12), (A.13), and (A.14),


n−1
E kAnt k22 | Gt+1
≤ kAtn−1 k22 (1 + ζtn−1 ) + µtn−1 − νtn−1 ,

where


2
 h,n
n−1
ζtn−1 = Cĥ,1 β̃tn−1 + 2 β̃tn−1 E kξt+1
(s, a)k2 | Gt+1
,



2
2
h,n
n−1
µn−1
= 3 β̃tn−1 E kξt+1
(s, a)k22 | Gt+1
+ Cĥ,2 β̃tn−1 ,
t

νtn−1 = 2 ε β̃tn−1 (Atn−1 )T ht θ̄tn−1 (s, a) .

P
n−1 2
Note that ∞
< ∞ almost surely by Assumption 4(iii) and the monotone conn=1 β̃t
vergence theorem (we can apply
logic as Lemma 5). This, together with Lemma
P the same
n−1
5, allows us to conclude that ∞
ζ
<
∞ almost surely. Moreover, since β̃tn−1 must
n=1 t 

h,n
n−1
converge to zero and (A.9) implies E kξt+1
(s, a)k22 | Gt+1
converges to zero, we can see
P∞
n−1
that n=1 µt
< ∞ almost surely. The well-known supermartingale
lemma
P∞ convergence
n−1
n k converges and
of Robbins and Siegmund
[1971]
tells
us
that
kA
ν
<
∞
almost
2
t
n=1 t

surely. Let Dδ = limn→∞kAnt k2 > δ . On the event Dδ , using strict convexity, we know
that (Atn−1 )T ht θ̄tn−1 (s, a) is positive and bounded away from zero for any n. Along with
P∞
P
n−1
n−1
= ∞, this implies ∞
= ∞, which in turn shows us that Dδ must occur
n=1 β̃t
n=1 νt
with probability zero for any δ > 0.

34

References
B. Acciaio and I. Penner. Dynamic risk measures. In Advanced Mathematical Methods for
Finance, pages 1–34. Springer Berlin Heidelberg, 2011.
W. A. Al-Qaq, M. Devetsikiotis, and J. K. Townsend. Stochastic gradient optimization of
importance sampling for the efficient simulation of digital communication systems. IEEE
Transactions on Communications, 43(12):2975–2985, 1995.
P. Artzner, F. Delbaen, J. Eber, and D. Heath. Coherent measures of risk. Mathematical
Finance, 9:203–228, 1999.
P. Artzner, F. Delbaen, J. Eber, D. Heath, and H. Ku. Coherent multiperiod risk adjusted
values and Bellman’s principle. Annals of Operations Research, 152(1):5–22, 2006.
M. G. Azar, M. Ghavamzadeh, and H. J. Kappen. Reinforcement learning with a near
optimal rate of convergence. Technical Report, 2011.
O. Bardou and N. Frikha. Computing VaR and CVaR using stochastic approximation and
adaptive unconstrained importance sampling. Monte Carlo Methods and Applications,
15(3):173–210, 2009.
J. Belles-Sampera, M. Guillén, and M. Santolino. Beyond value-at-risk: GlueVaR distortion
risk measures. Risk Analysis, 34(1):121–134, 2014.
D. Bertimas and A. W. Lo. Optimal control of execution costs. Journal of Financial
Markets, 1(1):1–50, 1998.
D. P. Bertsekas and J. N. Tsitsiklis. Neuro-dynamic programming. Athena Scientific, Belmont, MA, 1996.
K. Boda and J. A. Filar. Time consistent dynamic risk measures. Mathematical Methods
of Operations Research, 63(1):169–186, 2006.
S. Boyd and L. Vandenberghe. Convex optimization. Cambridge University Press, 2004.
L. Breiman. Probability. Society of Industrial and Applied Mathematics, Philadelphia, PA,
1992.
J. Bucklew. Introduction to rare event simulation. Springer New York, 2004.
H. Byström. Extreme value theory and extremely large electricity price changes. International Review of Economics & Finance, 14(1):41–55, 2005.
O. Çavus and A. Ruszczyński. Computational methods for risk-averse undiscounted transient markov models. Operations Research, 62(2):401–417, 2014.
P. Cheridito and M. Stadje. Time-inconsistency of VaR and time-consistent alternatives.
Finance Research Letters, 6(1):40–46, 2009.
P. Cheridito, F. Delbaen, and M. Kupper. Dynamic monetary risk measures for bounded
discrete-time processes. Electronic Journal of Probability, 11(3):57–106, 2006.

35

Y. Chow and M. Ghavamzadeh. Algorithms for CVaR optimization in MDPs. In Advances
in Neural Information Processing Systems, pages 3509–3517, 2014.
K. L. Chung. On a stochastic approximation method. The Annals of Mathematical Statistics, 25(3):463–483, 1954.
J. Danielsson, B. Jorgensen, G. Samorodnitsky, M. Sarma, and C. de Vries. Subadditivity
reexamined: The case for Value-at-risk. Preprint, London School of Economics, 2005.
J. Dhaene, R. J. A. Laeven, S. Vanduffel, G. Darkiewicz, and M. J. Goovaerts. Can a
coherent risk measure be too subadditive? Journal of Risk and Insurance, 75(2):365–
386, 2006.
K. Dowd and D. Blake. After VaR: The theory, estimation, and insurance applications of
quantile-based risk measures. Journal of Risk and Insurance, 73(2):193–229, 2006.
D. Duffie and J. Pan. An overview of value at risk. The Journal of Derivatives, 4(3):7–49,
1997.
D. Egloff and M. Leippold. Quantile estimation with adaptive importance sampling. Annals
of Statistics, 38(2):1244–1278, 2010.
J. Enders, W. B. Powell, and D. Egan. A dynamic model for the failure replacement of
aging high-voltage transformers. Energy Systems, 1(1):31–59, 2010.
E. Even-Dar and Y. Mansour. Learning rates for Q-learning. Journal of Machine Learning
Research, 5:1–25, 2004.
M. Frittelli and E. R. Gianin. Dynamic convex risk measures. In Risk measures for the 21st
century, volume 8, pages 227–248. 2004.
A. P. George and W. B. Powell. Adaptive stepsizes for recursive estimation with applications
in approximate dynamic programming. Machine Learning, 65(1):167–198, 2006.
P. Glasserman and T. W. Liu. Rare-event simulation for multistage production-inventory
systems. Management Science, 42(9):1292–1307, 1996.
R. Ibragimov and J. Walden. The limits of diversification when losses may be large. Journal
of Banking and Finance, 31(8):2551–2569, 2007.
T. Jaakkola, M. I. Jordan, and S. P. Singh. On the convergence of stochastic iterative
dynamic programming algorithms. Neural Computation, 6:1185–1201, 1994.
D. R. Jiang and W. B. Powell. An approximate dynamic programming algorithm for monotone value functions. Operations Research, 63(6):1489–1511, 2015a.
D. R. Jiang and W. B. Powell. Optimal hour-ahead bidding in the real-time electricity market with battery storage using approximate dynamic programming. INFORMS Journal
on Computing, 27(3):525–543, 2015b.
A. Juditsky, G. Lan, A. Nemirovski, and A. Shapiro. Robust stochastic approximation
approach to stochastic programming. SIAM Journal on Optimization, 19(4):1574–1609,
2009.
36

Y. S. Kan. On the convergence of a stochastic approximation procedure for estimating the
quantile criterion in the case of a discontinuous distribution function. Automation and
Remote Control, 72(2):283–288, 2011.
M. J. Kearns and S. P. Singh. Finite-sample convergence rates for Q-learning and indirect
algorithms. In Advances in Neural Information Processing Systems 11, volume 12, pages
996–1002, 1999.
J. H. Kim and W. B. Powell. An hour-ahead prediction model for heavy-tailed spot prices.
Energy Economics, 33(6):1252–1266, 2011.
A. J. Kleywegt, A. Shapiro, and T. Homem-de Mello. The sample average approximation
method for stochastic discrete optimization. SIAM Journal on Optimization, 12(2):479–
502, 2002.
V. Kozmı́k and D. P. Morton. Evaluating policies in risk-averse multi-stage stochastic
programming. Mathematical Programming, 2014.
H. J. Kushner and G. G. Yin. Stochastic approximation and recursive algorithms and
applications, volume 35. Springer, 2003.
N. Löhndorf, D. Wozabal, and S. Minner. Optimizing trading decisions for hydro storage
systems using approximate dual dynamic programming. Operations Research, 61(4):810–
823, 2013.
J. M. Nascimento and W. B. Powell. An optimal approximate dynamic programming
algorithm for the lagged asset acquisition problem. Mathematics of Operations Research,
34(1):210–237, 2009.
F. Neise. Risk management in stochastic integer programming. Springer Fachmedien, 2008.
M. V. F. Pereira and L. M. V. G. Pinto. Multi-stage stochastic optimization applied to
energy planning. Mathematical Programming, 52:359–375, 1991.
G. C. Pflug. Optimization of stochastic models. Kluwer Academic Publishers, 1996.
G. C. Pflug and A. Ruszczyński. Measuring risk for income streams. Computational Optimization and Applications, 32(1-2), 2005.
A. B. Philpott and V. L. de Matos. Dynamic sampling algorithms for multi-stage stochastic
programs with risk aversion. European Journal of Operational Research, 218(2):470–483,
2012.
A. B. Philpott, V. L. de Matos, and E. Finardi. On solving multistage stochastic programs
with coherent risk measures. Operations Research, 61(4):957–970, 2013.
W. B. Powell. Approximate dynamic programming: Solving the curses of dimensionality.
Wiley, 2nd edition, 2011.
W. B. Powell, A. Ruszczyński, and H. Topaloglu. Learning algorithms for separable approximations of discrete stochastic optimization problems. Mathematics of Operations
Research, 29(4):814–836, 2004.

37

M. L. Puterman. Markov decision processes: Discrete stochastic dynamic programming.
Wiley, New York, 2014.
A. Rakhlin, O. Shamir, and K. Sridharan. Making gradient descent optimal for strongly
convex stochastic optimization. Proceedings of the 29th International Conference on Machine Learning (ICML-12), pages 449–456, 2012.
F. Riedel. Dynamic coherent risk measures. Stochastic Processes and their Applications,
112(2):185–200, 2004.
H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical
Statistics, pages 400–407, 1951.
H. Robbins and D. Siegmund. A convergence theorem for nonnegative almost supermartingales and some applications. In J. S. Rustagi, editor, Optimizing Methods in Statistics,
pages 233–257. Academic Press, New York, 1971.
R. T. Rockafellar and S. Uryasev. Optimization of conditional value-at-risk. Journal of
Risk, 2:21–41, 2000.
R. T. Rockafellar and S. Uryasev. Conditional value-at-risk for general loss distributions.
Journal of Banking and Finance, 26(7):1443–1471, 2002. ISSN 03784266. doi: 10.1016/
S0378-4266(02)00271-6.
R. Y. Rubinstein. The simulated entropy method for combinatorial and continuous optimization. Methodology and Computing in Applied Probability, 1(2):127–190, 1999.
B. Rudloff, A. Street, and D. M. Valladão. Time consistency and risk averse dynamic
decision models: definition, interpretation and practical consequences. European Journal
of Operational Research, 234(3):743–750, 2014.
A. Ruszczyński. Risk-averse dynamic programming for Markov decision processes. Mathematical Programming, 125(2):235–261, 2010.
A. Ruszczyński and A. Shapiro. Conditional risk mappings. Mathematics of Operations
Research, 31(3):544–561, 2006a.
A. Ruszczyński and A. Shapiro. Optimization of convex risk functions. Mathematics of
Operations Research, 31(3):433–452, 2006b.
E. K. Ryu and S. P. Boyd. Adaptive importance sampling via stochastic convex programming. Working Paper, Stanford University, 2015.
I. O. Ryzhov and W. B. Powell. Information collection on a graph. Operations Research,
59(1):188–201, 2011.
I. O. Ryzhov, P. I. Frazier, and W. B. Powell. A new optimal stepsize for approximate
dynamic programming. IEEE Transactions on Automatic Control, 60(3):743–758, 2015.
T. Schaul, S. Zhang, and Y. LeCun. No more pesky learning rates. Journal of Machine
Learning Research, 28(2):343–351, 2013.

38

E. N. Sereda, E. M. Bronshtein, S. T. Rachev, F. J. Fabozzi, E. W. Sun, and S. V. Stoyanov.
Distortion risk measures in portfolio optimization. In Handbook of Portfolio Construction,
pages 649–673. Springer US, 2010.
A. Shapiro. On a time consistency concept in risk averse multistage stochastic programming.
Operations Research Letters, 37(3):143–147, 2009.
A. Shapiro, W. Tekaya, J. P. da Costa, and M. P. Soares. Risk neutral and risk averse
stochastic dual dynamic programming method. European Journal of Operational Research, 224(2):375–391, 2013.
C. Szepesvari and M. L. Littman. Generalized Markov decision processes: dynamic programming and reinforcement learning algorithms. Technical Report, 1996.
L. Tierney. A space-efficient recursive procedure for estimating a quantile of an unknown
distribution. SIAM Journal on Scientific and Statistical Computing, 4(4):706–711, 1983.
J. N. Tsitsiklis. Asynchronous stochastic approximation and Q-learning. Machine Learning,
16(3):185–202, 1994.
C. J. Watkins and P. Dayan. Q-learning. Machine Learning, 8(3-4):279–292, 1992.
P. Whittle. Multi-armed bandits and the Gittins Index. J. Roy. Statistical Soc. Series B
(Methodological), 42(2):143–149, 1980.
X. Xi, R. Sioshansi, and V. Marano. A stochastic dynamic programming model for cooptimization of distributed energy storage. Energy Systems, 5(3):475–505, 2014.

39

