Central-limit approach to risk-aware Markov decision processes

arXiv:1512.00583v1 [math.OC] 2 Dec 2015

Pengqian Yu∗†

Jia Yuan Yu∗‡

Huan Xu

§

Abstract
Whereas classical Markov decision processes maximize the expected reward, we consider
minimizing the risk. We propose to evaluate the risk associated to a given policy over a longenough time horizon with the help of a central limit theorem. The proposed approach works
whether the transition probabilities are known or not. We also provide a gradient-based policy
improvement algorithm that converges to a local optimum of the risk objective.

1

Introduction

Markov Decision Processes (MDPs) are essential models for stochastic sequential decision-making
problems (e.g., Puterman, 2014; Bertsekas, 1995; Sutton and Barto, 1998). Classical MDPs are
concerned with the expected performance criteria. However, in many practical problems, risk-neutral
objectives may not be appropriate (cf. Ruszczyński, 2010, Example 1).
Risk-aware decision-making is prevalent in the financial mathematics and optimization literature
(e.g., Artzner et al., 2002; Rockafellar, 2007; Markowitz et al., 2000; Von Neumann and Morgenstern,
2007), but limited to single-stage decisions. Risk-awareness has been adopted in multi-stage or
sequential decision problems more recently. A chance-constrained formulation for MDPs with uncertain parameters has been discussed in (Delage and Mannor, 2010). A criteria of maximizing the
probability of achieving a pre-determined target performance has been analyzed in (Xu and Mannor,
2011). In (Ruszczyński, 2010), Markov risk measure is introduced and a corresponding risk-averse
policy iteration method is developed. A mean-variance optimization problem for MDPs is addressed
in (Mannor and Tsitsiklis, 2013). A generalization of percentile optimization with objectives defined
by a measure over percentiles instead of a single percentile is introduced in (Chen and Bowling, 2012).
In terms of computational convenience, actor-critic algorithms for optimizing variance-related risk
measures in both discounted and average reward MDPs have been proposed in (Prashanth and Ghavamzadeh,
2013). Mean and conditional value-at-risk optimization problem in MDPs is solved by policy gradient
and actor-critic algorithms in (Chow and Ghavamzadeh, 2014). A unified policy gradient approach
is proposed in (Tamar et al., 2015) to seek optimal strategy in MDPs with the whole class of coherent
risk measures. Risk-aversion in multiarmed bandit problems is studied in (Yu and Nikolova, 2013),
where algorithms are proposed with PAC guarantees for best arm identification. Robust MDPs (e.g.,
Nilim and El Ghaoui, 2005; Iyengar, 2005) mitigate the sensitivity of optimal policy to ambiguity
∗ This work was supported in part by the EU FP7 project INSIGHT under grant 318225 while J. Yu and P. Yu
were at IBM Research Ireland.
† P. Yu is with the Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive
1, Singapore 117575, Singapore (e-mail: yupengqian@nus.edu.sg).
‡ J. Yu is with Concordia Institute of Information System Engineering, Concordia University, 1455 de Maisonneuve
Blvd. West EV007.635, Montreal, Quebec H3G 1M8, Canada (e-mail: jy@osore.ca).
§ H. Xu is with the Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive
1, Singapore 117575, Singapore (e-mail: mpexuh@nus.edu.sg).

1

in the underlying transition probabilities, and are closely related to the risk-sensitive MDPs. Such
relations are uncovered by Osogami (2012).
In this work, we study the value-at-risk in finite-horizon MDPs. Computing this risk associated with a sequence of decision using the defintion and first principles is intractable. However, we
show that this computation is made tractable by using a central limit theorem for Markov chains
(Jones et al., 2004; Kontoyiannis and Meyn, 2003; Glynn and Meyn, 1996). For a long-enough horizon T and a fixed policy π, we are thus able to evaluate the risk associated with following policy
π over T time steps. Specifically, our first contributions are policy evaluation algorithms whether
the transition probability matrix induced by π is known or not. Under mild conditions, we provide
high-probability error bounds for the evaluation.
For a fixed risk measure ρ, and a space of policies that is parametrized by θ, our second contribution is a policy improvement algorithm that converges in finite iterations, under certain assumptions,
to a locally optimal policy. This approach updates the parameter θ in the direction of gradient of the
risk measure. Compared to the previous work, our proposed method does not explicitly approximate
the value function. Therefore, it does not have approximation error due to the selection of value
function approximator.
Even though we deal mainly with the static value-at-risk risk measure (defined by the cumulative
reward), our results can be easily extended to deal with the conditional- or average-value-at-risk by
using its definition Schied (2004).

1.1

Distinction from related works

It is important to note that our risk measure is very different from the dynamic risk measure analyzed
in (Ruszczyński, 2010), which has the following recursive structure: for T -length horizon MDP and
policy π, the dynamic risk measure ρT is defined as



,
r(xπ0 ) + γρ r(xπ1 ) + γρ r(xπ2 ) + γρ . . . + r(xπT )

where r is a state dependent reward function, γ < 1 is a constant and ρ is a Markov risk measure,
i.e., the evaluation of each static coherent risk measure ρ is not allowed to depend on the whole past.
In contrast, our risk measure emphasizes the statistical properties of the reward accumulated over
multiple time steps: it is similar to that used in (Yu and Nikolova, 2013) in the context of bandit
problems.
From the computational and algorithmic perspective, we evaluate the policy and obtain the risk
value directly. Different from the work (Prashanth and Ghavamzadeh, 2013; Chow and Ghavamzadeh,
2014; Tamar et al., 2015), which indirectly evaluated the risk value by using function approximation,
our proposed method has an explicit policy evaluation step and does not require value function approximation. Therefore, it does not have approximation error due to the selection of value function
approximator and the richness of the features. Moreover, for risk-aware MDPs that are known to be
NP-hard (e.g., Delage and Mannor, 2010; Xu and Mannor, 2011), our method can find approximated
solutions to those problems.
From the conceptual perspective, our work is the first attempt to consider general risk measures with a central-limit approach. Compared with (Tamar et al., 2015), our approach is not
limited to so-called “coherent” risk measures (both static and dynamic), whose defining properties are not satisfied by the value-at-risk. The methods proposed in (Tamar et al., 2015) thus
cannot be applied. Moreover, the time horizon in our setup is allowed to be infinite. Other previous work (Delage and Mannor, 2010; Xu and Mannor, 2011; Mannor and Tsitsiklis, 2013) restrict
themselves to specific risk criteria: Delage and Mannor (2010) considered a set of percentile criteria; Xu and Mannor (2011) seek to find the policy that maximizes the probability of achieving
2

a pre-determined target performance, and Mannor and Tsitsiklis (2013) discussed a mean-variance
optimization problem.
This paper is organized as follows: In Section 2 we provide background and problem statement.
For a fixed policy, we then evaluate its risk value in Section 3. In specific, we discuss the policy
evaluation algorithms in Section 3.1 and 3.2 when the transition kernel is either known or unknown.
The main result is Theorem 2, which bounds the error of true and estimated asymptotic variances.
In Section 4, we provide policy gradient algorithm to improve the fixed policy. A conclusion is
offered in Section 5.

2

Problem formulation

MDP is a tuple hX , A, P, r, T i, where X is a finite set of states, A is a finite set of actions, P is the
transition kernel, T is the (possibly infinite) moderate time horizon, and r(x, a) : X × A → R is
the bounded and measurable reward function. A function F : X → R is called lattice if there are
h > 0 and 0 ≤ d < h, such that F (x)−d
is an integer for x ∈ X . The minimal h for such condition
h
holds is called the span of F . If the function F can be written as a sum F = F0 + Fl , where Fl
is lattice with span h and F0 has zero asymptotic variance (which will be defined later), then F
is called almost-lattice. Otherwise, F is called strongly nonlattice. We assume the functional r is
strongly nonlattice. A policy π is the rule according to which actions are selected at each state.
We parameterize the stochastic policies by {π(·|x; θ), x ∈ X , θ ∈ Θ ⊆ Rk1 }, and denote the set of
policies by Π. Since in this setting a policy π is represented by its k1 -dimensional parameter vector
θ, policy dependent functions can be written as a function of θ in place of π. So, we use π and θ
interchangeably in the paper.
Given a policy π ∈ Π, a Markov reward process X π = {xπt : x ∈ X , t ∈ T } is induced by π. We
write the transition kernel of the induced Markov process and dependent rewards as P π and r(xπt )
for all t ∈ T , respectively. We denote (P π )n as the n-step Markov transition kernel corresponding
to P π . We are interested in the long-term behavior of the sum
STπ =

T
−1
X

r(xπt ),

t=0

T ≥ 1.

If the Markov chain is positive recurrent with invariant probability measure ξ π , we define the mean
of reward function r as
"T −1
#
X
1
ϕπ (r) = lim
E
r(xπt ) ,
(1)
T →∞ T
t=0
where xπ0 is the initial condition and E is the expectation over the Markov process X π . Equivalently,
the mean can be expressed as
X
ϕπ (r) =
ξ π (x)r(x).
(2)
x

Often we can quantify the rate of convergence of (1) by showing that the limit
"T −1
#
X
π
π
π
r̂ (x) = lim E
r(xt ) − T ϕ (r)
T →∞

t=0

exists where, in fact, the function r̂π (x) solves the Poisson equation
P π r̂π = r̂π − r + ϕπ (r)1.
3

(3)

P
Here P π operates on functions r̂π : X → R via P π r̂π (x) = y P (x, y)r̂π (y), r is the vector form of
r(x) and 1 is the vector with all elements being 1.
Let Y denote the set of bounded random variables, and let ρ : Y → R denote a risk measure.
Our risk-aware objective is
min ρ (STπ ) .
(4)
π∈Π

We discuss specific types of risk measures in the coming sections. All proofs are deferred in the
supplementary material due to space constraints.

3

Policy evaluation

In this section, we develop methods to evaluate the risk value ρ(STπ ) for a fixed policy π. In specific,
we consider the Markov chain X π induced by policy π, which is assumed to be geometrically ergodic.
We then apply limit theorems for geometrically ergodic Markov chain in (Kontoyiannis and Meyn,
2003) to obtain the distribution function of cumulative reward when transition kernel is either known
or unknown. Finally, we provide algorithms for computing some examples of risk measures.
Similar to Kontoyiannis and Meyn (2003), we make the following two standard assumptions for
the Markov process X π .
Assumption 1. The Markov process X π is geometrically ergodic with Lyapunov function V (i.e.,
ψ−irreducible, aperiodic and satisfying Condition (V4) of Kontoyiannis and Meyn (2003)) where
V : X → [1, ∞) satisfies ϕπ (V 2 ) < ∞.
The detailed definition of geometric ergodicity can be found in (Kontoyiannis and Meyn, 2003,
Section 2).
Assumption 2. The (measurable)
reward function r : X → [−1, 1] has nontrivial asymptotic vari√
ance σ 2 , limt→∞ varx (Stπ / t) > 0.
Under Assumption 1 and 2, Kontoyiannis and Meyn (2003) showed that the following result
holds.
Theorem 1 (Theorem 5.1, Kontoyiannis and Meyn (2003)). Suppose that X π and the strongly
nonlattice functional r satisfy Assumption 1 and 2, and let GπT (y) denote the distribution function
S π −T ϕπ (r)
:
of the normalized partial sums T σ√T
GπT (y)

,P




STπ − T ϕπ (r)
√
≤y .
σ T

Then, for all xπ0 ∈ X and as T → ∞,

i
1
γ(y) h ̺
(1 − y 2 ) − r̂π (xπ0 ) + o(T − 2 ),
GπT (y) = g(y) + √
2
σ T 6σ

uniformly in y ∈ R, where γ(y) denotes the standard normal density and g(y) √
is the corresponding
distribution function. Here ̺ is a constant related to the third moment of STπ / T . The definitions
of ̺ and term o(T −1/2 ) can be found in the supplementary material.
Unlike (Heyman and Sobel, 2003), which considered finite Markov chains and obtained the moments of the cumulative reward, we further impose geometric ergodicity assumption on the Markov
4

process X π . By doing this, we are able to utilize the explicit form of distribution function given by
(Kontoyiannis and Meyn, 2003, Theorem 5.1) with the time error term o(T −1/2 ). In addition, the
geometric ergodicity property enables us to bound the approximation error in Theorem 2.
Similar to Glynn and Meyn (1996), we define the fundamental kernel Z π associated with fixed π
as
Z π = (H π )−1 ,
where the kernel Ξπ is defined as Ξπ (x, ·) = ξ π (·), and H π , I − P π − Ξπ . If the inverse of H π is
well defined, Glynn and Meyn (1996) stated that the solution to Poisson equation (3) has the form
r̂π = Z π (r − ϕπ (r)1).

(5)

When σ 2 < ∞ and ϕπ (r̂π ) < ∞, Meyn and Tweedie (2012, Equation 17.44) showed that the
asymptotic variance σ 2 can be obtained by
X
E[(Snπ − nϕπ (r))2 ]
=
[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x).
n→∞
n

σ 2 = lim

(6)

x∈X

We remark that the asymptotic variance can also be calculated without the knowledge of r̂π by
Theorem 17.5.3 of Meyn and Tweedie (2012):
σ2 =

X

xπ
0 ∈X

r 0 2 ξ π (xπ0 ) + 2

∞ X
X
k=1

xπ
0 ∈X

r 0 ξ π (xπ0 ) ×

X

xπ
k ∈X

rk P π (x, xπk |x = xπ0 )k ,

where rk = r(xπk ) − ϕπ (r) and P π (x, y|x = x0 )n denotes the n-step transition probability to state y
conditioned on the initial state x0 . It is clear that solving (6) requires the knowledge of transition
kernel P π . In the following subsections, we study the policy evaluation with and without knowing
P π , respectively.

3.1

Transition probability kernel P π is known

In this subsection, we propose an algorithm to obtain the distribution function GπT in Theorem 1,
using which we can evaluate the risk value.
Since transition kernel P π is known, we first calculate the stationary distribution ξ π . Note that
the stationary distribution satisfies
P π ξπ = ξπ .
Thus, the stationary distribution is the eigenvector of matrix P π with eigenvalue 1 (recall that 1 is
always an eigenvalues for matrix P π ), which is unique up to constant multiples under Assumption
1.
After we compute the stationary distribution ξ π , it is straightforward to calculate the mean (2),
solve Poisson equation (3) and get the asymptotic variance (6). Therefore, we are able to obtain the
distribution function GπT (y) defined in Theorem 1 at this stage. Finally, the risk measures that can
be expressed as a function of GπT (y) can be evaluated. We provide two examples below.
Example 3.1 (T -step value-at-risk). Let λ be given and fixed. The T -step value-at-risk is defined
by
ρλ (STπ ) = V aRλ (STπ ) = −q(λ),
where q is the right-continuous quantile function1 of STπ .
1 Formally,

q(λ) = inf{y ∈ R : F (y) > λ}, where F is the distribution function of x.

5

When T is moderate, o(T −1/2 ) vanishes, the distribution function GπT (y) defined in Theorem 1
has the form
i
γ(y) h ̺
(1 − y 2 ) − r̂π (xπ0 ) .
(7)
GπT (y) ≈ g(y) + √
2
σ T 6σ

This allows us to compute the T -step value-at-risk. The procedure is summarized in Algorithm 1.
The next example considers the mean-variance risk measure.
Example 3.2 (Mean-variance risk). Given λ, the mean-variance risk measure is given by
√
ρλ (STπ ) = −µ(STπ ) + λσ 2 (STπ / T ),
√
√
where µ(STπ ) is the mean of STπ , and σ 2 (STπ / T ) is the variance of STπ / T .

The corollary below gives formula to obtain the mean-variance risk. An algorithm can be obtained
similarly for the mean-variance risk simply by replacing line 7 in Algorithm 1 by the corresponding
risk measure. The same procedure can easily be generalized to other risk measures.
Corollary 1. Under Assumption 1 and 2, given λ, the mean-variance risk measure can be computed
by
ρλ (STπ ) = −ϕ(r) + λσ 2 .
where ϕ(r) and σ 2 are defined in (2) and (6), respectively.

Algorithm 1 Policy evaluation for T -step value-at-risk (P π is known)
1: Input: Transition probability matrix P π induced by policy π, initial state xπ
0 , decision horizon
T ∈ N+ , reward function r, and the coefficient λ of value-at-risk.
2: Obtain stationary distribution ξ π by solving the equation P π ξ π = ξ π .
P π
3: Substitute ξ π in (2) for mean ϕπ (r) =
x ξ (x)r(x).
4: Solve the Poisson equation (3) for r̂π by (5).
P
π
2
π π
2 π
5: Compute asymptotic variance in (6): σ 2 =
x∈X [r̂ (x) − (P r̂ (x)) ]ξ (x).
6: Calculate the constant ̺.
7: Output: The T -step value-at-risk, which can be computed by
V]
aRλ (STπ ) = − inf {y ∈ R : GπT (h(y)) > λ} ,
where GπT (y) is defined in (7), and h(y) =

3.2

y−T ϕπ (r)
√
.
σ T

Transition probability kernel P π needs to be estimated

In practice, knowing transition kernel P π is generally hard. In this subsection, we present a policy
evaluation algorithm when transition kernel needs to be estimated, and provide the corresponding
approximation error.
Intuitively, if the estimated transition kernel is accurate enough, the gap between solutions to
estimated and true Poisson equations is small (as shown in Lemma 5), which implies the difference
between estimated and true asymptotic variances is also small (as shown in Theorem 2).
A way of estimating transition kernel is to use the empirical distribution. Denoting the required number of samples by n1 and the estimated transition kernel by Pnπ1 , we make the following
assumption.
6

Assumption 3. There exist ǫ1 (n1 ) > 0 and δ1 (n1 ) > 0 such that ǫ1 (n1 ) → 0 and δ1 (n1 ) → 0 as
n1 → ∞. Moreover, we have the following bound for the error2 :

P ||Pnπ1 − P π ||∞ ≤ ǫ1 (n1 ) ≥ 1 − δ1 (n1 ).

Here the supreme norm || · ||∞ of a matrix A is defined as the maximum absolute row sum of the
matrix:
n
X
|aij |.
||A||∞ , max
1≤i≤m

j=1

We remark Assumption 3 can be easily satisfied if a large collection of observed data is provided.
This could be achievable if one is allowed to simulate the system arbitrarily many times.
Next, we perturb the true transition kernel P π , define the resulting perturbed transition kernel
as P̃ π and denote the stationary distribution associated with P̃ π as ξ˜π . Given P̃ π , we let ξnπ (·) =
P̃ π (x, ·)n be an estimator of ξ˜π (·). Under Assumption 1, we make the following assumption which
establishes the convergence of ||ξnπ − ξ̃ π ||, where || · || is the total variation norm defined as following:
If µ is a signed measure on B(X ), then the total variation norm is defined as
||µ|| = sup |µ(f )| =
|f |≤1

sup µ(A) −

A∈B(X )

inf

A∈B(X )

µ(A).

Assumption 4. Given ǫ > 0 and perturbed transition kernel P̃ π ∈ P . If ||P̃ π − P π ||∞ ≤ ǫ, then
for every n ∈ N+ , the geometrically ergodic Markov chain has a form of convergence
||P̃ π (x, ·)n − ξ˜π (·)|| ≤ M (x)τ n ,

where τ < 1, and M (x) is a nonnegative function.
This assumption essentially guarantees the geometric ergodicity property remains if the original
Markov chain is slightly perturbed.
Assumption 3 and 4 immediately imply that


(8)
P ||Pnπ1 (x, ·)n − ξ˜π (·)|| ≤ M (x)τ n ≥ 1 − δ1 (n1 ),

where ξ˜π is the stationary distribution associated with Pnπ1 . Let τ1 (P π ) be the ergodicity coefficient
of transition kernel P π (Seneta, 1988), which is defined by
τ1 (P π ) = sup ||ν ⊤ P π ||1 .

(9)

||ν||1 =1
ν ⊤ e=0

If the ergodicity coefficient satisfies 0 ≤ τ1 (P π ) < 1, Cho and Meyer (2001, Equation 3.4) shows
||ξ˜π − ξ π ||1 ≤
Therefore, under Assumption 3, we have

P ||ξ˜π − ξ π ||1 ≤

||Pnπ1 − P π ||∞
.
1 − τ1 (P π )

ǫ1 (n1 )
1 − τ1 (P π )



≥ 1 − δ1 (n1 ).

(10)

The ergodicity coefficient τ1 (Pnπ1 ) for the estimated transition kernel Pnπ1 is defined in a same manner.
Matrix H π plays an important role in defining solution to Poisson equation. We impose a
requirement on the smallest singular value σmin (H π ) of H π .
2 We

use ǫ1 (n1 ) and δ1 (n1 ) to show the dependence of the error on n1 .

7

Assumption 5. There exists a constant c > ̟(x)
e
such that the smallest singular value of matrix
H π satisfies σmin (H π ) ≥ c. Here3


p
ǫ1 (n1 )
n2
+ 2M (x)τ
.
̟(x)
e
= |X | ǫ1 (n1 ) +
(11)
1 − τ1 (Pnπ1 ) − ǫ1 (n1 )
Similar to (3), we define the estimated Poisson equation as

Pnπ1 r̂nπ1 ,n2 = r̂nπ1 ,n2 − r + ϕπn2 (r)1,

(12)

where ϕπn2 (r) is the estimated mean
ϕπn2 (r) =

X

ξnπ2 (x)r(x),

(13)

x

with ξnπ2 being the estimation of stationary distribution ξ˜π
ξnπ2 (·) = Pnπ1 (x, ·)n2 .

(14)

We then define the perturbed fundamental kernel as
Znπ1 ,n2 = (Hnπ1 ,n2 )−1 ,
where Hnπ1 ,n2 = I − Pnπ1 − Ξπn2 with Ξπn2 (x, ·) = ξnπ2 (·). The following lemma shows Znπ1 ,n2 is well
defined under Assumption 5.
Lemma 1. Under Assumption 5, Znπ1 ,n2 is well defined.
The solution to estimated Poisson equation (12) can be expressed as
r̂nπ1 ,n2 = Znπ1 ,n2 (r − ϕπn2 (r)1),

(15)

We define the asymptotic variance associated with r̂n1 ,n2 as σn2 1 ,n2 , which can be obtained by
the equation below:
X
σn2 1 ,n2 =
(16)
[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξnπ2 (x).
x∈X

The theorem below provides the condition for n2 such that the difference between estimated and
true asymptotic variances σ 2 and σn2 1 ,n2 defined in (6) and (16) is small. We provide an outline of
the proof below and a complete version can be found in the supplementary material.

Theorem 2. Given policy π, ǫ ≥ 0, λ ∈ (0, 1) and let σn2 1 ,n2 and σ 2 be the estimated and true
asymptotic variances, respectively. Under Assumption 1, 2, 3, 4 and 5, there exists n2 such that the
following two conditions (C1) and (C2) are satisfied
p
|X |(|X | + 1)3 (̟
e − |X |ǫ1 (n1 ))
(C1)
≤ λǫ
(c − ̟)
e 2



p 
α2 (x)|X | (|X | + 1)̟
e p
+ max α1 (x)α2 (x)|X | ≤ (1 − λ)ǫ,
(C2)
e − ǫ1 (n1 ) |X |
+ |X | ̟
x∈X
c
c−̟
e
3 The parameter ̟
e is dependent on the initial state x through nonnegative function M (x). In the following, we
omit x for convenience and emphasize x when such dependence is required to show.

8

where

p

3 

p 
ǫ1 (n1 ) |X |(|X | + 1) |X | 2 ̟(|X
e
| + 1) p
,
e − ǫ1 (n1 ) |X |
+
+ |X | ̟
c−̟
e
c
c−̟
e
|X |(|X | + 1)(2c − ̟)
e
α2 (x) =
.
c(c − ̟)
e

α1 (x) =

Moreover, we have

P(|σn2 1 ,n2 − σ 2 | ≤ ǫ) ≥ 1 − 38δ1 (n1 ).

Proof. By equations (6) and (16), we have

|σn2 1 ,n2 − σ 2 |
X
X
[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξnπ2 (x) −
[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x)
=
x∈X

(1)

=

X

x∈X

+

x∈X

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξnπ2 (x) −

X

x∈X
(2)

X

x∈X

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x) −

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x)

X

x∈X

[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x)

≤ an + b n ,

where
an =

X

x∈X

bn =

X

x∈X

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ](ξnπ2 (x) − ξ π (x)) ,
[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x) −

X

x∈X

[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x) .

P
Equality (1) holds due to subtracting and then adding a same term x∈X [r̂nπ1 ,n2 (x)2 −(Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x).
Inequality (2) holds by triangle inequality. In order to bound the terms an and bn , we need following
lemmas.
First, we derive a bound for |τ1 (Pnπ1 ) − τ1 (P π )| by the lemma below.
Lemma 2. Under Assumption 3, the ergodicity coefficient for the estimated transition kernel Pnπ1
satisfies

P |τ1 (Pnπ1 ) − τ1 (P π )| ≤ ǫ1 (n1 ) ≥ 1 − δ1 (n1 ).

Denoting ||a||2 as the Euclidean norm of vector a, we present a lemma showing that the solution
(5) of Poisson equation (3) is bounded.
Lemma
3. Under Assumption 2 and 5, the solution to Poisson equation r̂π is bounded ||r̂π ||2 ≤
√
|X |(|X |+1)
.
c
Analogously to Lemma 3, we have the following result.
Lemma 4. Under Assumption 1, 2, 3, 4 and 5, the solution (15) to estimated Poisson equation
(12) is bounded such that
!
p
|X |(|X | + 1)
π
≥ 1 − 4δ1 (n1 ).
P ||r̂n1 ,n2 ||2 ≤
c−̟
e
9

The lemma below shows that solutions to the true and estimated Poisson equations are close to
each other with high probability.
Lemma 5. Under Assumption 1, 2, 3, 4 and 5, the difference between solutions to true and estimated
Poisson equation is bounded such that

P ||r̂π − r̂nπ1 ,n2 ||2 ≤ ς˜ ≥ 1 − 7δ1 (n1 ),
where ς˜ =

√

|X |(|X |+1)̟
e
c(c−̟)
e

+

|X |
c



p 
̟
e − ǫ1 (n1 ) |X | .

Recall that for two random variables an and bn , and ǫ ∈ R with λ ∈ (0, 1), by union bound and
monotonicity property of probability, we have
P(an + bn ≤ ǫ) ≥ P(an ≤ λǫ) + P(bn ≤ (1 − λ)ǫ).
Finally, by bounding terms an and bn , we prove the theorem.
Now we are able to compute the T -step value-at-risk, which is summarized as a procedure in
the Algorithm 2. Note that given the asymptotic variance σn1 ,n2 and solution to estimated Poisson
equation r̂n1 ,n2 , when T is moderate, o(T −1/2 ) vanishes, GπT (y) defined in Theorem 1 becomes


γ(y)
̺˜
π
π
2
π
√
GT (y) ≈ g(y) +
(17)
(1 − y ) − r̂n1 ,n2 (x0 ) ,
σn1 ,n2 T 6σn2 1 ,n2
where ̺˜ is defined in a same manner as ̺ by replacing P π and ξ π by Pnπ1 and ξnπ2 , respectively.
In the following, we analysis the mean-variance risk defined in Example 3.2. Corollary 2 gives
formula to obtain the mean-variance risk and provides high-probability bound for the risk measure.
An algorithm can be developed for the mean-variance risk simply by replacing line 9 in Algorithm
2 by the corresponding risk measure.
Corollary 2. Under Assumption 1, 2, 3, 4 and 5, given λ, the estimated mean-variance risk measure
can be computed by
ρ̃λ (STπ ) = −ϕπn2 (r) + λσn2 1 ,n2 .

where ϕπn2 (r) and σn2 1 ,n2 are defined in (13) and (16), respectively. Moreover, we have
P (|ρ̃λ (STπ ) − ρλ (STπ )| ≤ κ) ≥ 1 − 41δ1 (n1 ),

where ρλ (STπ ) is the true mean-variance risk value computed in Corollary 1 and κ = 2M (x)τ n2 +
ǫ1 (n1 )
1−τ1 (P π )−ǫ1 (n1 ) + λǫ.
n1

A similar argument can be developed for T -step value-at-risk defined in Example 3.1. When
̺ = 1, r̂ = 0.5, T = 100, ϕ(r) = 0.1, we plot the distribution functions GT for σ = 1 and 1.1
below. Figure 1 shows the distribution functions are close, which implies the corresponding T -step
value-at-risks are also close to each other (When λ = 0.3, the T -step value-at-risks are 5.134 and
4.632, respectively).

10

1
0.8

GT

0.6
0.4
0.2
0

σ=1
σ=1.1
−20

−10

0

10

y

20

30

40

Figure 1: Distribution functions GT for σ = 1 and 1.1.
Algorithm 2 Policy evaluation for T -step value-at-risk (P π is unknown)
1: Input: Policy π, initial state xπ
0 , reward function r, decision horizon T ∈ N+ , λ ∈ (0, 1) and
estimation parameters ǫ1 (n1 ) and δ1 (n1 ).
2: Obtain the estimated transition kernel Pnπ1 satisfying Assumption 3 with parameter ǫ1 (n1 ) and
δ1 (n1 ).
3: Compute the ergodicity coefficient τ1 (Pnπ1 ) for the estimated transition kernel Pnπ1 .
4: Choose n2 with τ1 (Pnπ1 ) obtained from step 3 such that conditions in Theorem 2 hold.
π
defined by (14).
5: Compute the estimated stationary distribution ξn
2
π
by (15).
6: Calculate the solution of estimated Poisson equation r̂n
1 ,n2
2
7: Obtain the asymptotic variance σn1 ,n2 defined in (16).
8: Calculate the constant ̺
˜.
9: Output: The T -step value-at-risk, which can be computed by
V]
aRλ (STπ ) = − inf {y ∈ R : GπT (h(y)) > λ} ,
where GπT (y) is defined in (17), and h(y) =

4

y−T ϕπ
n2 (r)
√
σn1 ,n2 T

with ϕπn2 (r) defined in (13).

Policy Improvement

In this section, we propose policy gradient methods to improve policy π and solve the optimization
problem (4). Specifically, we compute the gradient of the performance (i.e., the risk) w.r.t. the policy
parameters from the evaluated policies, and improve π by adjusting its parameters in the direction
of the gradient.
In the following, we make a standard assumption in gradient-based MDPs literature (e.g.,
Sutton et al., 1999; Konda and Tsitsiklis, 1999; Bhatnagar et al., 2009; Prashanth and Ghavamzadeh,
2013).
Assumption 6. For any state-action pair (x, a), π(a|x; θ) is continuously differentiable in the parameter θ.
In a recent work, Prashanth and Ghavamzadeh (2013, Section 4) presented actor-critic algorithms for optimizing the risk-sensitive measure that are based on two simultaneous perturbation
11

methods: simultaneous perturbation stochastic approximation (SPSA) (Spall, 1992) and smoothed
functional (SF) (Bhatnagar et al., 2012). We use similar arguments here and further apply those
methods to estimate the gradient of the risk measure with respect to the policy parameter θ, i.e.,
∇θ ρ(STθ ).
The idea in these methods is to estimate the gradient ∇θ ρ(STθ ) using two simulated trajectories
of the system corresponding to policies with parameters θ and θ+ = θ + β△. Here β > 0 is a positive
constant and △ is a perturbation random variable, i.e., a k1 -vector of independent Rademacher (for
SPSA) and Gaussian N (0, 1) (for SF) random variables.
SPSA-based estimate for ∇θ ρ(STθ ) is given by
∂θ(i) ρ(STθ ) ≈

ρ(STθ+β△ ) − ρ(STθ )
,
β△(i)

i = 1, . . . , k1

(18)

where △ is a vector of independent Rademacher random variables and △(i) is its i-th entry. The
advantage of this estimator is that it perturbs all directions at the same time (the numerator is
identical in all k1 components). So, the number of function measurements needed for this estimator
is always two, independent of the dimension k1 .
SF-based estimate for ∇θ ρ(STθ ) is given by
∂θ(i) ρ(STθ ) ≈


△(i) 
ρ(STθ+β△ ) − ρ(STθ ) , i = 1, . . . , k1
β

(19)

where △ is a vector of independent Gaussian N (0, 1) random variables.
The step size of the gradient descent {λ(t)} satisfies the following condition.
Assumption 7. The step size schedule {λ(t)} satisfies
X
λ(t) > 0,
λ(t) = ∞, and
t

X
t

λ(t)2 < ∞.

At each time step t, we update the policy parameter θ as follows: for i = 1, · · · , k1 ,


(i)
(i)
θt+1 = Γi θt − λ(t)∂θ(i) ρ(STθ ) .

(20)

(i)

Note that △t ’s are independent Rademacher and Gaussian N (0, 1) random variables in SPSA and
SF updates, respectively. Γ is an operator that projects a vector θ ∈ Rk1 to the closest point in a
compact and convex set C ⊂ Rk1 . The projection operator is necessary to ensure the convergence
of the algorithms. Specifically, we consider the differential equation
θ̇t = Γ̌(∇θ ρ(STθt )),

(21)

where Γ̌ is defined as follows: for any bounded continuous function f (·),
Γ̌(f (θt )) = lim

τ →0

Γ(θt + τ f (θt )) − θt
.
τ

The projection operator Γ̌(·) ensures that the evolution of θ via the differential equation (21) stays
within the bounded set C ∈ Rk1 . Let Z = {θ ∈ C : Γ̌(∇θ ρ(STθt )) = 0} denote the set of asymptotically stable equilibrium points of the differential equation (21) and Z ǫ = {θ ∈ C : ||θ − θ0 || < ǫ, θ0 ∈
Z} denote the set of points in the ǫ-neighborhood of Z.
The following is our main result, which bounds the error of our policy improvement algorithm.
12

Theorem 3. Under Assumption 1, 6 and 7, for any given ǫ > 0, there exists β0 > 0 such that for
all β ∈ (0, β0 ), θt → θ∗ ∈ Z ǫ almost surely.
The above theorem guarantees the convergence of the SPSA and SF updates to a local optimum of
the risk objective function ρ(STθ ). Finally, we provide the policy improvement method in Algorithm
3.
Algorithm 3 Policy improvement
1: Initialize policy π parameterized by parameter θ1 and stopping criteria ǫ. Let t = 1 and θ0 =
θ1 + 2ǫ1.
2: while ||θt − θt−1 ||∞ > ǫ do
3:
Evaluate the risk value ρ(STθt ) and ρ(STθt +β△t ) associated with policy θt and θt + β△t by
Algorithm 1 or 2. Here △t are independent Rademacher and Gaussian N (0, 1) random variables
in SPSA and SF updates, respectively.
4:
For i = 1, · · · , k1 , update the policy parameter


(i)
(i)
θt+1 = Γi θt − λ(t)∂θ(i) ρ(STθ ) ,
5:
6:
7:

where λ(t) is the step size and ∂θ(i) ρ(STθ ) is defined in (18) for SPSA or (19) for SF.
t = t + 1.
end while
Output: Policy parameter θ.

5

Conclusion

In the context of risk-aware Markov decision processes, we use a central limit theorem efficiently
evaluate the risk associated to a given policy. We also provide a policy improvement algorithm,
which works on a parametrized policy-space in a gradient-descent way. Under mild conditions, it
is guaranteed that the policy evaluation is approximately correct and that the policy improvement
converges to a local optimum of the risk measure. We believe that many important problems that are
usually addressed using standard MDP models should be revisited using our proposed approach to
risk management (e.g., machine replacement, inventory management, portfolio optimization, etc.).

A

Appendix

A.1

Defintion of ̺ in Theorem 1

̺ = ̺1 + ̺3 + ̺3 ,
where
̺1 =

X

r(xπ0 )3 ξ π (xπ0 ),

xπ
0 ∈X

̺2 = 3

∞
X
X

i=−∞
i6=0

xπ
0 ∈X

r(xπ0 )2 ξ π (xπ0 ) ×

X

xπ
i ∈X

r(xπi )P π (x, xπi |x = xπ0 )i ,

13

and
̺3 = 6

∞
X
X

i,j=1

xπ
0 ∈X

r(xπ0 )ξ π (xπ0 ) ×

X

xπ
i ∈X

r(xπi )P π (x, xπi |x = xπ0 )i , ×

X

xπ
j ∈X

r(xπi+j )P π (x, xπj |x = xπi )j .

Here P π (x, y|x = x0 )n denotes the n-step transition probability to state y conditioned on the initial
state x0 .

A.2

Definition of o(T −1/2 ) in Theorem 1

Kontoyiannis and Meyn (2003) showed that the term o(T −1/2 ) can be represented as
o(T

−1/2

1
)=
π

Z

√
A T
√

−A T

MT



iω
√
σ T



− φT (ω)

ν
dω
+√
|ω|
T

where
MT (α) , mT (α) exp(−αE[STπ ]),
with mT (α) = E[exp(αSTπ )], T ≥ 1 and α ∈ C. The characteristic function φT (·) is defined as



−ω 2
̺(iω)3
√
, ω ∈ R.
φT (ω) , exp
1+
2
6σ 3 T


′
E[S π ]
A is chosen large enough so that A > 24(νπ)−1 |ΨT y − σ√TT | for all y ∈ R, T ≥ 1 and arbitrary
ν > 0. Here,
̺
√ (1 − y 2 )γ (y) , y ∈ R.
ΨT (y) = g (y) +
3
6σ T

A.3

Proof of Lemma 1

e π (x, ·) = ξ˜π (·) and
Proof. Let Ξ

△1 = P π − Pnπ1 ,
eπ ,
△2 = Ξπ − Ξ

e π − Ξπ ,
△3 = Ξ
n2

and ||A||2 be the spectral norm of matrix A, which is the largest singular value of matrix A:
||A||2 = σmax (A).
Due to Assumption 3, (8) and (10), we have

p 
P ||△1 ||2 ≤ ǫ1 (n1 ) |X | ≥ 1 − δ1 (n1 ),
!
p
ǫ1 (n1 ) |X |
≥ 1 − δ1 (n1 ),
P ||△2 ||2 ≤
1 − τ1 (P π )


p
P ||△3 ||2 ≤ 2 |X |M (x)τ n2 ≥ 1 − δ1 (n1 ).

Define △ = △1 + △2 + △3 , we have

P (||△||2 ≤ ̟(x)) ≥ 1 − 3δ1 (n1 ),
14

(22)

where4
̟(x) =


p
|X | ǫ1 (n1 ) +

ǫ1 (n1 )
+ 2M (x)τ n2
1 − τ1 (P π )



.

We perform the singular value decomposition for the matrix H π and have
H π = U ΣV ∗ = U Σ1/2 IΣ1/2 V ∗ ,

where U and V ∗ are unitary matrix with dimension |X | × |X |, and Σ is a |X | × |X | diagonal matrix
with non-negative real numbers on the diagonal. Next, by the definition of Znπ1 ,n2 and △, we obtain
Znπ1 ,n2 =(H π + △)−1

=(U Σ1/2 IΣ1/2 V ∗ + U Σ1/2 Σ−1/2 U −1 △(V ∗ )−1 Σ−1/2 Σ1/2 V ∗ )−1
e −1 (U Σ1/2 )−1 ,
=(Σ1/2 V ∗ )−1 (I + △)

e , Σ−1/2 U −1 △(V ∗ )−1 Σ−1/2 . In addition, we have
where △
Since ||Σ−1 ||2 = ||(H π )−1 ||2 =

e 2 = ||Σ−1 ||2 ||△||2 .
||△||

1
σmin (H π ) ,

under Assumption 5, we obtain



e 2 ≤ ̟ ≥ 1 − 3δ1 (n1 ).
P ||△||
c

(23)

e −1 is well defined, which implies Z π
Since ̟ < ̟
e < c, (I + △)
n1 ,n2 is well defined.

A.4

Proof of Lemma 2

Proof. From the definition of τ1 (P π ) and τ1 (Pnπ1 ), we have
|τ1 (Pnπ1 ) − τ1 (P π )|
≤

sup ||ν ⊤ Pnπ1 ||1 − sup ||ν ⊤ P π ||1

||ν||1 =1
ν ⊤ e=0

(1)

≤ sup

||ν||1 =1
ν ⊤ e=0

||ν||1 =1
ν ⊤ e=0

||ν ⊤ Pnπ1 ||1 − ||ν ⊤ P π ||1

(2)

≤ sup ||ν ⊤ (Pnπ1 − P π )||1
||ν||1 =1
ν ⊤ e=0

(3)

≤ sup ||ν||1 ||Pnπ1 − P π ||∞
||ν||1 =1
ν ⊤ e=0
=||Pnπ1 −

P π ||∞

4 The parameter ̟ is dependent on the initial state x through nonnegative function M (x). In the following, we
omit x for convenience and emphasize x when such dependence is required to show. Same argument goes for ̟.
e

15

Here recall that if f, g : A → R are bounded functions, then supA f − supA g ≤ supA |f − g|. For
||ν||1 = 1 and ν ⊤ e = 0, since ||ν ⊤ Pnπ1 ||1 ≤ 1 and ||ν ⊤ P π ||1 ≤ 1, thus inequality (1) holds. Inequality
(2) holds due to reverse triangle property, while Hölder’s inequality is applied to (3).
From Assumption 3, we get

P |τ1 (Pnπ1 ) − τ1 (P π )| ≤ ǫ1 (n1 ) ≥ 1 − δ1 (n1 ).

A.5

Proof of Lemma 3

Proof. Recall the solution can be bounded above by
||r̂π ||2 = ||Z π (r − ϕπ (r)1)||2
≤ ||Z π ||sp ||r − ϕπ (r)1||2
p
(1)
|X |(|X | + 1)
≤
.
c
p
Here (1) is due to ||r − ϕπ (r)1||2 ≤ |X |(|X | + 1) and Assumption 5. Since for x ∈ X ,
|r(x) − ϕπ (r)| = r(x) −

X

ξ π (x′ )r(x′ )

x′ ∈X

≤ |r(x)| +
≤ |X | + 1,

X

ξ π (x′ )|r(x′ )|

x′ ∈X

we have ||r−ϕπ (r)1||∞ ≤ |X |+1, which implies ||r−ϕπ (r)1||2 ≤
to show ||P π r̂π ||2 is also bounded.

p
|X |(|X |+1). It is straightforward

||P π r̂π ||2 ≤ ||P π ||2 ||r̂π ||2
p
≤ |X |||P π ||∞ ||r̂π ||2
≤

A.6

|X |(|X | + 1)
.
c

Proof of Lemma 4

Proof. Recall the solution to estimated Poisson equation can be expressed as
||r̂nπ1 ,n2 ||2 = ||Znπ1 ,n2 (r − ϕπn2 (r)1)||2

≤ ||Znπ1 ,n2 ||2 ||r − ϕπn2 (r)1||2 .

First, note that ||H π − Hnπ1 ,n2 ||2 = ||△||2 . From (22) in the proof of Lemma 1, we have

P ||H π − Hnπ1 ,n2 ||2 ≤ ̟ ≥ 1 − 3δ1 (n1 ).
16

By Weyl’s inequality, we have

which implies


P |σmin (H π ) − σmin (Hnπ1 ,n2 )| ≤ ̟ ≥ 1 − 3δ1 (n1 ).

P ||Znπ1 ,n2 ||2 ≤

Moreover, for x ∈ X ,

1
c−̟



≥ 1 − 3δ1 (n1 ).

X

|r(x) − ϕπn2 (r)| = r(x) −

ξnπ2 (x′ )r(x′ )

x′ ∈X

≤ |r(x)| +
≤ |X | + 1,

X

ξnπ2 (x′ )|r(x′ )|

x′ ∈X

p
we have ||r − ϕπn2 (r)1||∞ ≤ |X | + 1. Thus ||r − ϕπn2 (r)1||2 ≤ |X |(|X | + 1). Therefore, we obtain
!
p
|X |(|X | + 1)
π
≥ 1 − 3δ1 (n1 ).
P ||r̂n1 ,n2 ||2 ≤
c−̟
By Lemma 2, we have
P(̟ ≤ ̟)
e ≥ 1 − δ1 (n1 ).

Applying (24), we have
P

||r̂nπ1 ,n2 ||2

!
p
|X |(|X | + 1)
≤
≥ 1 − 4δ1 (n1 ).
c−̟
e

It is straightforward to show ||Pnπ1 r̂nπ1 ,n2 ||2 is bounded because

Thus, we have

A.7

||Pnπ1 r̂nπ1 ,n2 ||2 ≤||Pnπ1 ||2 ||r̂nπ1 ,n2 ||2
p
≤ |X |||Pnπ1 ||∞ ||r̂nπ1 ,n2 ||2 .



|X |(|X | + 1)
π π
P ||Pn1 r̂n1 ,n2 ||2 ≤
≥ 1 − 4δ1 (n1 ).
c−̟
e

Proof of Lemma 5

Proof. From the proof of Lemma 1, recall that Znπ1 ,n2 can be expressed as
e −1 (U Σ1/2 )−1 .
Znπ1 ,n2 = (Σ1/2 V ∗ )−1 (I + △)

e 2 < 1 under Assumption 5, we have
Since ||△||

e −1 =
(I + △)

∞
X

e n.
(−1)n △

n=0

17

(24)

Furthermore, we have

||Znπ1 ,n2 − Z π ||2
= (Σ1/2 V ∗ )−1

∞
X

!

e n (U Σ1/2 )−1
(−1)n △

n=1

≤||(Σ1/2 V ∗ )−1 ||2 ||(U Σ1/2 )−1 ||2
e 2
||△||
≤
||Σ−1 ||2 .
e 2
1 − ||△||

∞
X

n=1

2

e n
||△||
2

Under Assumption 5 and (23) in the proof of Lemma 1, we have


̟
e
π
π
≥ 1 − 4δ1 (n1 ).
P ||Z − Zn1 ,n2 ||2 ≤
c(c − ̟)
e

The difference of solutions to true and estimated Poisson equation can be expressed as
||r̂π − r̂nπ1 ,n2 ||2

=||Z π (r − ϕπ (r)1) − |Znπ1 ,n2 (r − ϕπn2 (r)1)||2

=||Z π (r − ϕπ (r)1) − Z π (r − ϕπn2 (r)1) + Z π (r − ϕπn2 (r)1) − Znπ1 ,n2 (r − ϕπn2 (r)1)||2

=||(Z π − Znπ1 ,n2 )(r − ϕπn2 (r)1) + Z π (ϕπn2 (r) − ϕπ (r))1||2

≤||(Z π − Znπ1 ,n2 )(r − ϕπn2 (r)1)||2 + ||Z π (ϕπn2 (r) − ϕπ (r))1||2

≤||Z π − Znπ1 ,n2 ||2 ||r − ϕπn2 (r)1||2 + ||Z π ||2 ||(ϕπn2 (r) − ϕπ (r))1||2

≤||Z π − Znπ1 ,n2 ||2 ||r − ϕπn2 (r)1||2 + ||Z π ||2 (||(ϕπ (r) − ϕ̃π (r))1||2 + || ϕ̃π (r) − ϕπn2 (r))1||2
For x ∈ X ,

|r(x) − ϕπn2 (r)| = r(x) −

X

x′ ∈X

≤|r(x)| +
≤|X | + 1,

X

ξnπ2 (x′ )r(x′ )
ξnπ2 (x′ )|r(x′ )|

x′ ∈X

p
we have ||r − ϕπn2 (r)1||∞ ≤ |X | + 1. Thus ||r − ϕπn2 (r)1||2 ≤ |X |(|X | + 1).
(10) implies


ǫ1 (n1 )
π
π
˜
P ||ξ − ξ ||2 ≤
≥ 1 − δ1 (n1 ).
1 − τ1 (P π )
p
Since ||r||2 ≤ |X | and |ϕ̃π (r) − ϕπ (r)| ≤ ||ξ˜π − ξ π ||2 ||r||2 , we have
!
p
ǫ1 (n1 ) |X |
π
π
≥ 1 − δ1 (n1 ),
P |ϕ̃ (r) − ϕ (r)| ≤
1 − τ1 (P π )
which further yields


ǫ1 (n1 )|X |
P ||(ϕ̃π (r) − ϕπ (r))1||2 ≤
≥ 1 − δ1 (n1 )
1 − τ1 (P π )
18



p
due to the fact that ||(ϕ̃π (r) − ϕπ (r))1||2 ≤ |X |||(ϕ̃π (r) − ϕπ (r))1||∞ .
Moreover, (8) implies


P ||ξnπ2 − ξ˜π ||2 ≤ 2M (x)τ n ≥ 1 − δ1 (n1 ),

Similarly, we have


P ||(ϕ̃π (r) − ϕπn2 (r))1||2 ≤ 2M (x)τ n |X | ≥ 1 − δ1 (n1 ).

Therefore, we obtain

where ς =

where ς˜ =

A.8

√

|X |(|X |+1)̟
e
c(c−̟)
e

√

|X |(|X |+1)̟
e
c(c−̟)
e


P ||r̂π − r̂nπ1 ,n2 ||2 ≤ ς ≥ 1 − 6δ1 (n1 )


ǫ1 (n1 )
n
. Finally, due to Lemma 2, we have
+ |Xc | 1−τ
π ) + 2M (x)τ
(P
1

+

|X |
c


P ||r̂π − r̂nπ1 ,n2 ||2 ≤ ς˜ ≥ 1 − 7δ1 (n1 ),


p 
̟
e − ǫ1 (n1 ) |X | and ̟
e is defined in (11).

Proof of Theorem 2

Proof. Recall that the xth element of vector P π r̂π is denoted by P π r̂π (x). Define the vectors
(r̂π )2 and (P π r̂π )2 with their xth elements being (r̂π )2 (x) = r̂π (x)2 and (P π r̂π )2 (x) = (P π r̂π (x))2 ,
respectively. Similarly for the vectors (r̂nπ1 ,n2 )2 and (Pnπ1 r̂nπ1 ,n2 )2 . Lemma 3 implies
|X |(|X | + 1)2
,
c2
|X |2 (|X | + 1)2
,
||(P π r̂π )2 ||∞ ≤
c2

||(r̂π )2 ||∞ ≤

and from Lemma 4 we get


|X |(|X | + 1)2
≥ 1 − 4δ1 (n1 ),
P ||(r̂nπ1 ,n2 )2 ||∞ ≤
(c − ̟)
e 2


|X |2 (|X | + 1)2
P ||(Pnπ1 r̂nπ1 ,n2 )2 ||∞ ≤
≥ 1 − 4δ1 (n1 ).
(c − ̟)
e 2
First, by equations (6) and (16), we have

|σn2 1 ,n2 − σ 2 |
X
X
[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x)
[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξnπ2 (x) −
=
x∈X

x∈X

(1)

=

X

x∈X

+

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξnπ2 (x) −

X

x∈X
(2)

X

x∈X

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x) −

≤ an + b n ,
19

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x)

X

x∈X

[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x)

(25)

where
an =

X

x∈X

bn =

X

x∈X

[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ](ξnπ2 (x) − ξ π (x)) ,
[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x) −

X

x∈X

[r̂π (x)2 − (P π r̂π (x))2 ]ξ π (x) .

P
Equality (1) holds due to subtracting and then adding a same term x∈X [r̂nπ1 ,n2 (x)2 −(Pnπ1 r̂nπ1 ,n2 (x))2 ]ξ π (x).
Inequality (2) holds by triangle inequality. In the following, we bound the terms an and bn individually. Term an can be further expressed as
(1)

an ≤

(2)

X

x∈X

|[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ](ξnπ2 (x) − ξ π (x))|

≤ ||(r̂nπ1 ,n2 )2 − (Pnπ1 r̂nπ1 ,n2 )2 ||∞ ||ξnπ2 − ξ π (·)||1

(3)

≤ ||(r̂nπ1 ,n2 )2 ||∞ + ||(Pnπ1 r̂nπ1 ,n2 )2 ||∞ ||ξnπ2 − ξ˜π ||1 + ||ξ˜π − ξ π ||1 .

Inequalities (1) and (3) hold due to triangle property, while (2) is true because of Hölder’s
inequality. From (8) and (10), we have


(26)
P ||ξnπ2 − ξ̃ π ||1 + ||ξ˜π − ξ π ||1 ≤ ι ≥ 1 − 2δ1 (n1 ).

where ι = 2

p
ǫ1 (n1 )
|X |M (x)τ n2 + 1−τ
π . In addition, under Lemma 2 we obtain
1 (P )


P ||ξnπ2 − ξ̃ π ||1 + ||ξ˜π − ξ π ||1 ≤ ι̃ ≥ 1 − 3δ1 (n1 ),

p
where ι̃ = 2 |X |M (x)τ n2 +

ǫ1 (n1 )
1−τ1 (Pnπ1 )−ǫ1 (n1 ) .

Combining (25) and (26), by the union bound of

probability and definition of ̟,
e we have

!
p
|X |(|X | + 1)3 (̟
e − |X |ǫ1 (n1 ))
P an ≤
≤ 1 − 11δ1 (n1 ).
(c − ̟)
e 2

Next, we bound the term bn .
(1)

bn ≤

(2)

X

x∈X

|{[r̂nπ1 ,n2 (x)2 − (Pnπ1 r̂nπ1 ,n2 (x))2 ] − [r̂π (x)2 − (P π r̂π (x))2 ]}ξ π (x)|

≤ ||(r̂nπ1 ,n2 )2 − (r̂π )2 ||∞ ||ξ π ||1 + ||(Pnπ1 r̂nπ1 ,n2 )2 − (P π r̂)2 ||∞ ||ξ π ||1

(3)

≤ |X |(||(r̂nπ1 ,n2 )2 − (r̂π )2 ||∞ + ||(Pnπ1 r̂nπ1 ,n2 )2 − (P π r̂)2 ||∞ ).

Inequality (1) holds due to triangle property, and (2) is due to Hölder’s inequality, while (3) is true
since ξ π (x) ∈ [0, 1]. Note that
|r̂nπ1 ,n2 (x)2 − r̂π (x)2 |

=|r̂nπ1 ,n2 (x) − r̂π (x)||r̂nπ1 ,n2 (x) + r̂π (x)|

≤|r̂nπ1 ,n2 (x) − r̂π (x)|(|r̂nπ1 ,n2 (x)| + |r̂π (x)|).
20

Furthermore, from Lemma 3, 4 and 5, we have

P ||(r̂nπ1 ,n2 )2 − (r̂π )2 ||∞ ≤ ϑ ≥ 1 − 11δ1 (n1 ).

where
ϑ=

(27)



p
|X |(|X | + 1)(2c − ̟)
e (|X | + 1)̟
e p
|X
|(
̟
e
−
ǫ
(n
)
|X
|)
.
+
1 1
c2 (c − ̟)
e
c−̟
e

In addition, we have

|(Pnπ1 r̂nπ1 ,n2 (x))2 − (P π r̂π (x))2 |

=|Pnπ1 r̂nπ1 ,n2 (x) − P π r̂π (x)| × |Pnπ1 r̂nπ1 ,n2 (x) + P π r̂π (x)|

(1)

≤ (|Pnπ1 r̂nπ1 ,n2 (x) − P π r̂nπ1 ,n2 (x)| + |P π r̂nπ1 ,n2 (x) − P π r̂π (x)|)(|Pnπ1 r̂nπ1 ,n2 (x)| + |P π r̂π (x)|)
!
X
X
′
′
π
′
π
π
′
π ′
π
′
π
=
(Pn1 (x, x ) − P (x, x ))r̂n1 ,n2 (x )| + |
P (x, x )(r̂n1 ,n2 (x ) − r̂ (x ))
×

(2)

≤

(3)

x′ ∈X
(|Pnπ1 r̂nπ1 ,n2 (x)|

X

+ |P r̂ (x)|)

|(Pnπ1 (x, x′ )

x′ ∈X
× (|Pnπ1 r̂nπ1 ,n2 (x)|

x′ ∈X

π π

π

′

− P (x, x

))r̂nπ1 ,n2 (x′ )|

+

X

x′ ∈X

+ |P π r̂π (x)|)

π

′

|P (x, x

)(r̂nπ1 ,n2 (x′ )

π

′

− r̂ (x ))|

!

≤ b1 × (|Pnπ1 r̂nπ1 ,n2 (x)| + |P π r̂π (x)|).

where

b1 = ||Pnπ1 (x, ·) − P π (x, ·)||1 ||r̂nπ1 ,n2 ||∞ + ||P π (x, ·)||1 ||r̂nπ1 ,n2 − r̂π ||∞ .

Inequalities (1) and (2) are due to triangle property, and (3) is because of Hölder’s inequality.
Since ||Pnπ1 (x, ·) − P π (x, ·)||1 ≤ ||Pnπ1 − P π ||∞ , Assumption 3 implies

P ||Pnπ1 (x, ·) − P π (x, ·)||1 ≤ ǫ1 (n1 ) ≥ 1 − δ1 (n1 ).

Since ||P π (x, ·)||1 ≤ |X |, from Lemma 4 and 5, we have

P b1 ≤ b1 ≥ 1 − 12δ1 (n1 ).
where

b1 =

p

3 

p 
ǫ1 (n1 ) |X |(|X | + 1) |X | 2 ̟(|X
e
| + 1) p
e − ǫ1 (n1 ) |X |
.
+
+ |X | ̟
c−̟
e
c
c−̟
e

Moreover, from Lemma 3 and 4, we obtain


P |Pnπ1 r̂nπ1 ,n2 (x)| + |P π r̂π (x)| ≤ η ≥ 1 − 4δ1 (n1 )

|X |(|X | + 1)(2c − ̟)
e
. Note ̟
e is dependent on the initial state x and denote
c(c − ̟)
e
p

3 

p 
ǫ1 (n1 ) |X |(|X | + 1) |X | 2 ̟(|X
e
| + 1) p
,
α1 (x) =
e − ǫ1 (n1 ) |X |
+
+ |X | ̟
c−̟
e
c
c−̟
e
|X |(|X | + 1)(2c − ̟)
e
α2 (x) =
,
c(c − ̟)
e

where η =

21

we conclude



P ||(Pnπ1 r̂nπ1 ,n2 )2 − (P π r̂π )2 ||∞ ≤ max α1 (x)α2 (x) ≥ 1 − 16δ1 (n1 ).
x∈X

(28)

At last, by the union bound of probability, (27) and (28) yield

P bn ≤ bn ≥ 1 − 27δ1 (n1 ),

where

bn =



p
e p
α2 (x)|X | (|X | + 1)̟
+ |X |(̟
e − ǫ1 (n1 ) |X |) + max α1 (x)α2 (x)|X |.
x∈X
c
c−̟
e

Recall that for two random variables an and bn , and ǫ ∈ R with λ ∈ (0, 1), by union bound and
monotonicity property of probability, we have
P(an + bn ≤ ǫ) ≥ P(an ≤ λǫ) + P(bn ≤ (1 − λ)ǫ).
Therefore, if for an

p
|X |(|X | + 1)3 (̟
e − |X |ǫ1 (n1 ))
≤ λǫ,
(c − ̟)
e 2

and for bn



p
α2 (x)|X | (|X | + 1)̟
e p
e − ǫ1 (n1 ) |X |) + max α1 (x)α2 (x)|X | ≤ (1 − λ)ǫ.
+ |X |(̟
x∈X
c
c−̟
e

We can conclude that

P(an + bn ≤ ǫ) ≥ 1 − 38δ1 (n1 ),

which further implies

P(|σn2 1 ,n2 − σ 2 | ≤ ǫ) ≥ 1 − 38δ1 (n1 )

by the monotonicity property and union bound of probability.

A.9

Proof of Corollary 2

Proof. From the definitions of ρ̃λ (STπ ) and ρλ (STπ ), we get
|ρ̃λ (STπ ) − ρλ (STπ )|

=|ϕ(r) − ϕπn2 (r) + λσn2 1 ,n2 − λσ 2 |
(1) X
[ξ π (x) − ξnπ2 (x)]r(x) + λ|σ̃ 2 − σ 2 |
≤
x∈X

(2)

≤

(3)

X

x∈X

[|ξ π (x) − ξ̃ π (x)| + |ξ˜π (x) − ξnπ2 (x)|]|r(x)| + λ|σn2 1 ,n2 − σ 2 |

≤ ||ξ π − ξ̃ π ||1 + ||ξ˜π − ξnπ2 ||1 + λ|σn2 1 ,n2 − σ 2 |.

Here inequalities (1) and (2) are due to triangle property. Inequality (3) is because Hölder’s inequality
and ||r||∞ ≤ 1. Moreover, since Pnπ1 is the ǫ-perturbed transition kernel under Assumption 4, by (8)
and (14), we have


(29)
P ||ξ˜π − ξ π ||1 ≤ 2M (x)τ n2 ≥ 1 − δ1 (n1 ).
n2

22

In addition, from (10) and Lemma 2, we get


ǫ1 (n1 )
P ||ξ˜π − ξ π ||1 ≤
≥ 1 − 2δ1 (n1 ).
1 − τ1 (Pnπ1 ) − ǫ1 (n1 )

(30)

Applying the union bound of probability to (29), (30) and Theorem 2, we have
P (|ρ̃λ (STπ ) − ρλ (STπ )| ≤ κ) ≥ 1 − 41δ1 (n1 ).
where κ = 2M (x)τ n2 +

A.10

ǫ1 (n1 )
1−τ1 (Pnπ1 )−ǫ1 (n1 )

+ λǫ.

Proof of Theorem 3

Proof. Note that, for fixed θ, we are able to evaluate the policy ρ(STθ ). Moreover, (20) can be
considered as a discretization of the differential equation (21).
For SPSA update, Z is an asymptotically stable attractor for the differential equation (21), with
ρ(STθ ) itself serving as a strict Lyapunov function. This can be inferred as follows
dρ(STθ )
= ∇θ ρ(STθ )θ̇ = ∇θ ρ(STθ )Γ̌(−∇θ ρ(STθ )) < 0.
dt
The claim for SPSA update now follows from (Kushner and Clark, 2012, Theorem 5.3.1).
The convergence analysis for SF update follows in a similar manner as SPSA update.

References
Philippe Artzner, Freddy Delbaen, Jean-Marc Eber, and David Heath. Coherent measures of risk.
Risk management: value at risk and beyond, page 145, 2002.
Dimitri P Bertsekas. Dynamic programming and optimal control, volume 1. Athena Scientific Belmont, MA, 1995.
Shalabh Bhatnagar, Richard S Sutton, Mohammad Ghavamzadeh, and Mark Lee. Natural actor–
critic algorithms. Automatica, 45(11):2471–2482, 2009.
Shalabh Bhatnagar, HL Prasad, and LA Prashanth. Stochastic recursive algorithms for optimization:
simultaneous perturbation methods, volume 434. Springer, 2012.
Katherine Chen and Michael Bowling. Tractable objectives for robust policy optimization. In
Advances in Neural Information Processing Systems, pages 2069–2077, 2012.
Grace E Cho and Carl D Meyer. Comparison of perturbation bounds for the stationary distribution
of a Markov chain. Linear Algebra and its Applications, 335(1):137–150, 2001.
Yinlam Chow and Mohammad Ghavamzadeh. Algorithms for CVaR optimization in MDPs. In
Advances in Neural Information Processing Systems, pages 3509–3517, 2014.
Erick Delage and Shie Mannor. Percentile optimization for Markov decision processes with parameter
uncertainty. Operations research, 58(1):203–213, 2010.
Peter W Glynn and Sean P Meyn. A Liapounov bound for solutions of the Poisson equation. The
Annals of Probability, pages 916–931, 1996.
23

Daniel P Heyman and Matthew J Sobel. Stochastic Models in Operations Research: Stochastic
Optimization, volume 2. Courier Corporation, 2003.
Garud N Iyengar. Robust dynamic programming. Mathematics of Operations Research, 30(2):
257–280, 2005.
Galin L Jones et al. On the Markov chain central limit theorem. Probability surveys, 1(299-320):
5–1, 2004.
Vijay R Konda and John N Tsitsiklis. Actor-critic algorithms. In Advances in Neural Information
Processing Systems, volume 13, pages 1008–1014, 1999.
Ioannis Kontoyiannis and Sean P Meyn. Spectral theory and limit theorems for geometrically ergodic
Markov processes. Annals of Applied Probability, pages 304–362, 2003.
Harold Joseph Kushner and Dean S Clark. Stochastic approximation methods for constrained and
unconstrained systems, volume 26. Springer Science & Business Media, 2012.
Shie Mannor and John N Tsitsiklis. Algorithmic aspects of mean–variance optimization in Markov
decision processes. European Journal of Operational Research, 231(3):645–653, 2013.
Harry M Markowitz, G Peter Todd, and William F Sharpe. Mean-variance analysis in portfolio
choice and capital markets, volume 66. John Wiley & Sons, 2000.
Sean P Meyn and Richard L Tweedie. Markov chains and stochastic stability. Springer Science &
Business Media, 2012.
Arnab Nilim and Laurent El Ghaoui. Robust control of Markov decision processes with uncertain
transition matrices. Operations Research, 53(5):780–798, 2005.
Takayuki Osogami. Robustness and risk-sensitivity in Markov decision processes. In Advances in
Neural Information Processing Systems, pages 233–241, 2012.
LA Prashanth and Mohammad Ghavamzadeh. Actor-critic algorithms for risk-sensitive MDPs. In
Advances in Neural Information Processing Systems, pages 252–260, 2013.
Martin L Puterman. Markov decision processes: discrete stochastic dynamic programming. John
Wiley & Sons, 2014.
R Tyrrell Rockafellar. Coherent approaches to risk in optimization under uncertainty. Tutorials in
operations research, 3:38–61, 2007.
Andrzej Ruszczyński. Risk-averse dynamic programming for Markov decision processes. Mathematical programming, 125(2):235–261, 2010.
A. Schied. Risk measures and robust optimization problems, 2004. Lecture notes.
E Seneta. Perturbation of the stationary distribution measured by ergodicity coefficients. Advances
in Applied Probability, pages 228–230, 1988.
James C Spall. Multivariate stochastic approximation using a simultaneous perturbation gradient
approximation. Automatic Control, IEEE Transactions on, 37(3):332–341, 1992.
Richard S Sutton and Andrew G Barto. Reinforcement learning: An introduction, volume 1. MIT
press Cambridge, 1998.
24

Richard S Sutton, David A McAllester, Satinder P Singh, Yishay Mansour, et al. Policy gradient
methods for reinforcement learning with function approximation. In Advances in Neural Information Processing Systems, volume 99, pages 1057–1063. Citeseer, 1999.
Aviv Tamar, Yinlam Chow, Mohammad Ghavamzadeh, and Shie Mannor. Policy gradient for
coherent risk measures. arXiv preprint arXiv:1502.03919, 2015.
John Von Neumann and Oskar Morgenstern. Theory of games and economic behavior. Princeton
university press, 2007.
Huan Xu and Shie Mannor. Probabilistic goal Markov decision processes. In Proceedings of the
Twenty-Second international joint conference on Artificial Intelligence-Volume Volume Three,
pages 2046–2052. AAAI Press, 2011.
Jia Yuan Yu and Evdokia Nikolova. Sample complexity of risk-averse bandit-arm selection. In
Proceedings of the Twenty-Third international joint conference on Artificial Intelligence, pages
2576–2582. AAAI Press, 2013.

25

