From Infinite to Finite Programs: Explicit Error Bounds with
Applications to Approximate Dynamic Programming

arXiv:1701.06379v2 [math.OC] 20 Feb 2017

PEYMAN MOHAJERIN ESFAHANI, TOBIAS SUTTER, DANIEL KUHN, AND JOHN LYGEROS

Abstract. We consider linear programming (LP) problems in infinite dimensional spaces that are in general
computationally intractable. Under suitable assumptions, we develop an approximation bridge from the
infinite-dimensional LP to tractable finite convex programs in which the performance of the approximation is
quantified explicitly. To this end, we adopt the recent developments in two areas of randomized optimization
and first order methods, leading to a priori as well as a posteriori performance guarantees. We illustrate
the generality and implications of our theoretical results in the special case of the long-run average cost and
discounted cost optimal control problems for Markov decision processes on Borel spaces. The applicability of
the theoretical results is demonstrated through a constrained linear quadratic optimal control problem and
a fisheries management problem.

Keywords. infinite-dimensional linear programming, Markov decision processes, approximate dynamic
programming, randomized and convex optimization
1. Introduction
Linear programming (LP) problems in infinite dimensional spaces appear in, among other areas, engineering, economics, operations research and probability theory [33, 1, 32]. Infinite LPs offer remarkable modeling
power, subsuming general finite dimensional optimization problems and the generalized moment problem as
special cases. They are, however, often computationally formidable, motivating the study of approximations
schemes.
A particularly rich class of problems that can be modeled as infinite LPs involves Markov decision processes
(MDP) and their optimal control. More often than not, it is impossible to obtain explicit solutions to MDP
problems, making it necessary to resort to approximation techniques. Such approximations are the core of
a methodology known as approximate dynamic programming [8, 6]. Interestingly, a wide range of optimal
control problems involving MDP can be equivalently expressed as static optimization problems over a closed
convex set of measures, more specifically, as infinite LPs [24, 26, 27]. This LP reformulation is particularly
appealing for dealing with unconventional settings involving additional constraints [23, 3], secondary costs
[18], information-theoretic considerations [43], and reachability problems [28, 34]. In addition, the infinite
LP reformulation allows one to leverage the developments in the optimization literature, in particular convex
approximation techniques, to develop approximation schemes for MDP problems. This will also be the
perspective adopted in the present article.
Approximation schemes to tackle infinite LPs have historically been developed for special classes of problems, e.g., the general capacity problem [31], or the generalized moment problem [32]. The literature on
control of MDP with infinite state or action spaces mostly concentrates on approximation schemes with
asymptotic performance guarantees [26, 25], see also the comprehensive book [30] for controlled stochastic
differential equations and [41, 35] for reachability problems in a similar setting. From a practical viewpoint, a
Date: February 22, 2017.
The
authors
are
with
the
Delft
Center
for
Systems
and
Control,
TU
Delft,
The
Netherlands
(P.MohajerinEsfahani@tudelft.nl),
the
Automatic
Control
Laboratory,
ETH
Zurich,
Switzerland,
({sutter,lygeros}@control.ee.ethz.ch), and the Risk Analytics and Optimization Chair, EPFL, Switzerland
(daniel.kuhn@epfl.ch).
1

2

challenge using these schemes is that the convergence analysis is not constructive and does not lead to explicit
error bounds. A wealth of approximation schemes have been proposed in the literature under the names of
approximate dynamic programming [5], neuro-dynamic programming [8], reinforcement learning [29, 46], and
value and/or policy iteration [6]. Most, however, deal with discrete (finite or at most countable) state and
action spaces, while approximation over uncountable spaces remains largely unexplored.
The MDP literature on explicit approximation errors in uncountable settings can, roughly speaking, be
divided to two groups in terms of the performance criteria considered: discounted cost, and average cost. Of
the two, the discounted cost setting has received more attention as the corresponding dynamic programming
operator is a contraction, a useful property to obtain a convergence rate for the approximation error. Examples include the linear programming approach [13, 14], and also a recent series of works [17, 18, 11] on
approximating a probability measure that underlies the random transitions of the dynamics of the system
using different discretization procedures. Long-run average cost problems introduce new challenges due to
loosing the contraction property. The authors in [19] develop approximation schemes leading to finite but nonconvex optimization problems, while [42] investigates the convergence rate of the finite-state approximation
to the original (uncountable) MDP problem.
The approach presented in this article tackles a class of general infinite LPs that, as a special case, cover
both long-run discounted and average cost performance criteria in the optimal control of MDP. The resulting
approximation is based on finite convex programs that are different from the existing schemes. Closest in
spirit to our proposed approximation is the linear programming approach based on constraint sampling in
[13, 14, 45]. Unlike these works, however, we introduce an additional norm constraint that effectively acts
as a regularizer. We study in detail the conditions under which this regularizer can be exploited to bound
the optimizers of the primal and dual programs, and hence provide an explicit approximation error for the
proposed solution.
The proposed approximation scheme involves a restriction of the decision variables from an infinite dimensional space to a finite dimensional subspace, followed by the approximation of the infinite number of
constraints by a finite subset; we develop two complementary methods for performing the latter step. The
structure of the article is illustrates in Figure 1, where the contributions are summarized as follows:
• We introduce a subclass of infinite LPs whose regularized semi-infinite restriction enjoys analytical
bounds for both primal and dual optimizers (Proposition 3.2). The implications for MDP with average
cost (Lemma 3.7) and with discounted cost (Lemma A.2) are also investigated.
• We derive an explicit error bound between the original infinite LP and the regularized semi-infinite
counterpart, providing insights on the impact of the underlying norm structure as well as on how the
choice of basis functions contributes to the approximation error (Theorem 3.3, Corollary 3.5). In the
MDP setting, we recover an existing result as a special case (Corollary 3.9).
• We adopt the recent developments from the randomized optimization literature to propose a finite
convex program whose solution enjoys a priori probabilistic performance bounds (Theorem 4.4).
We extend the existing results to offer also an a posteriori bound under a generic underlying norm
structure. The required conditions and theoretical assertions are validated in the MDP setting (Corollary 4.12).
• In parallel to the randomized approach, we also utilize the recent developments in the structural
convex optimization literature to propose an iterative algorithm for approximating the semi-infinite
program. For this purpose, we extend the setting to incorporate unbounded prox-terms with a certain
growth rate (Theorem 5.3). We illustrate how this extension allows us to deploy the entropy prox-term
in the MDP setting (Lemma 5.10, Corollary 5.8).
Section 2 introduces the main motivation for the work, namely the control of discrete-time MDP and
their LP characterization. Using standard results in the literature we embed these MDP in the more general
framework of infinite LPs. Section 3 studies the link from infinite LPs to semi-infinite programs. Section 4
presents the approximation of semi-infinite programs based on randomization, while Section 5 approaches

3

discrete time
MDP

equivalent

infinite LP
(P): J

Theorem 3.3
Proposition 3.2
robust program
(Pn ): Jn

strong duality

semi-infinite program
(Dn ): Jen

Theorem 4.4

Theorem 5.3

infinite program
semi-infinite programs

scenario program
(Pn,N ): Jn,N

regularized program
(Dn,η ): Jen,η

finite programs

Theorem 6.1
prior & posterior error
Jen,η − Jn,N

Figure 1. Graphical representation of the article structure and its contributions
the same objective using first-order convex optimization methods. Section 6 summarizes the results in the
preceding sections, establishing the approximation error from the original infinite LP to the finite convex
counterparts. Section 7 illustrates the theoretical results through a truncated LQG example and a fisheries
management problem.
Notation. The set R+ denotes the set of non-negative reals and k · kℓp for p ∈ [1, ∞] the standard p-norm in
Rn . Given a function u : S → R, we denote the infinity norm of the function by kuk∞ := sups∈S |u(s)|, and

′
)|
. The space of Lipschitz functions on a set S is
the Lipschitz norm by kukL := sups,s′ ∈S |u(s)|, |u(s)−u(s
ks−s′ kℓ∞
denoted by L (S); define the function 1(s) ≡ 1 for all s ∈ S. We denote the Borel σ-algebra on the (topological) space S by B(S). Measurability is always understood in the sense of Borel. Products of topological
spaces are assumed to be endowed with the product topology and the corresponding product σ-algebra. The
space of finite signed measures (resp. probability measures) on S is denoted by M(S) (resp. P(S)). The
R
Wasserstein norm on the space of signed measures M(S) is defined by kµkW := supkukL ≤1 S u(s)µ(ds) and
can be shown to be the dual of the Lipschitz norm. The set of extreme points of a set A is denoted by E{A}.
Given a bilinear form · , · , the support function of A is defined by σA (y) = supx∈A y, x . The standard
bilinear form in Rn (i.e., the inner product) is denoted by y x.

·

2. Motivation: Control of MDP and LP Characterization

2.1. MDP setting
We briefly recall some standard definitions and refer interested readers to [24, 23, 2] for further details.

Consider a Markov control model S, A, {A(s) : s ∈ S}, Q, ψ , where S (resp. A) is a metric space called the
state space (resp. action space) and for each s ∈ S the measurable set A(s) ⊆ A denotes the set of feasible
actions when the system is in state s ∈ S. The transition law is a stochastic kernel Q on S given the feasible
state-action pairs in K := {(s, a) : s ∈ S, a ∈ A(s)}. A stochastic kernel acts on real valued measurable

4

functions u from the left as
Qu(s, a) :=

Z

S

u(s′ )Q(ds′ |s, a),

and on probability measures µ on K from the right as
Z

Q(B|s, a)µ d(s, a) ,
µQ(B) :=
K

∀(s, a) ∈ K,

∀B ∈ B(S).

Finally ψ : K → R+ denotes a measurable function called the one-stage cost function. The admissible history
spaces are defined recursively as H0 := S and Ht := Ht−1 × K for t ∈ N and the canonical sample space
is defined as Ω := (S × A)∞ . All random variables will be defined on the measurable space (Ω, G) where G
denotes the corresponding product σ-algebra. A generic element ω ∈ Ω is of the form ω = (s0 , a0 , s1 , a1 , . . .),
where si ∈ S are the states and ai ∈ A the action variables. An admissible policy is a sequence π = (πt )t∈N0
of stochastic kernels πt on A given ht ∈ Ht , satisfying the constraints πt (A(st )|ht ) = 1. The set of admissible
policies will be denoted by Π. Given a probability measure ν ∈ P(S) and policy π ∈ Π, by the Ionescu Tulcea
theorem [7, p. 140-141] there exists a unique probability measure Pπν on (Ω, G) such that for all measurable
sets B ⊂ S, C ⊂ A, ht ∈ Ht , and t ∈ N0

Pπν s0 ∈ B = ν(B)

Pπν at ∈ C|ht = πt (C|ht )

Pπν st+1 ∈ B|ht , at = Q(B|st , at ).

The expectation operator with respect to Pπν is denoted by Eπν . The stochastic process Ω, G, Pπν , (st )t∈N0 is
called a discrete-time MDP. For most of the article we consider optimal control problems where the aim is
to minimise a long term average cost (AC) over the set of admissible policies and initial state measures. We
definite the optimal value of the optimal control problem by
" T −1
#
1 π X
AC
:=
inf
lim sup Eν
J
ψ(st , at ) .
(1)
(π,ν)∈Π×P(S) T →∞ T
t=0

We emphasize, however, that the results also apply to other performance objective, including the long-run
discounted cost problem as shown in Appendix A.

2.2. Infinite LP characterization
The problem in (1) admits an alternative LP characterization under some mild assumptions.
Assumption 2.1 (Control model). We stipulate that
(i) the set of feasible state-action pairs is the unit hypercube K = [0, 1]dim(S×A) ;
(ii) the transition law Q is Lipschitz continuous, i.e., there exists LQ > 0 such that for all k, k ′ ∈ K and
all continuous functions u
|Qu(k) − Qu(k ′ )| ≤ LQ kuk∞ kk − k ′ kℓ∞ ;
(iii) the cost function ψ is non-negative and Lipschitz continuous on K with respect to the ℓ∞ -norm.
Assumption 2.1(i) may seem restrictive, however, essentially it simply requires that the state-action set K
is compact. We refer the reader to Example 7.2 where a non-rectangular K is transferred to a hypercube,
and to [26, Chapter 12.3] for further information about the LP characterization in more general settings.
Theorem 2.2 (LP characterization [19, Proposition 2.4]). Under Assumption 2.1,


 inf −ρ
−J AC =

ρ,u

s. t. ρ + u(s) − Qu(s, a) ≤ ψ(s, a),


ρ ∈ R, u ∈ L (S).

∀(s, a) ∈ K

(2)

5


The LP (2) can be expressed in the standard conic form inf x∈X x, c : Ax − b ∈ K by introducing

X = R × L (S)




x
= (ρ, u) ∈ X





 c = (c1 , c2 ) = (−1, 0) ∈ R × M(S)
(3)
b(s, a) = −ψ(s, a)
R



u(s)c
(ds)
x,
c
=
c
ρ
+
2
1

S



Ax(s,
a)
=
−ρ
−
u(s)
+
Qu(s, a)



K = L+ (K),

where M(S) is the set of finite signed measures supported on S, and L+ (K) is the cone of Lipschitz functions
taking non-negative values. It should be noted that the choice of the positive cone K = L+ (K) is justified
since, thanks to Assumption 2.1(ii), the linear operator A maps the elements of X into L (K).
Remark 2.3 (Constrained MDP). The LP characterization of MDP naturally allows us to incorporate constraints in the form of
" T −1
#
1 π X
lim sup Eν
di (st , at ) ≤ ℓi ,
∀i ∈ {1, · · · , I},
T →∞ T
t=0

where the functions di : K → R and constants ℓi reflect our desired specifications. To this end, it suffices to
PI
introduce auxiliary decision variables βi ∈ R+ , and in (2) replace ρ in the objective with ρ − i=1 βi ℓi and
PI
in the constraint with ρ − i=1 βi di , see [23, Theorem 5.2].
Our aim is to derive an approximation scheme for a class of such infinite dimensional LPs, including
problems of the form (2), that comes with an explicit bound on the approximation error.
3. Infinite to Semi-infinite Programs

3.1. Dual pairs of normed vector spaces

The triple X, C, k · k is called a dual pair of normed vector spaces if

• X and C are vector spaces;
• · , · is a bilinear form on X × C that “separates points”, i.e.,
– for each nonzero x ∈ X there is some c ∈ C such that x, c 6= 0,
– for each nonzero c ∈ C there is some x ∈ X such that x, c 6= 0;
• X is equipped with the norm k · k, which together with the bilinear form induces a dual norm in C
defined through kck∗ := supkxk≤1 x, c .

The norm in the vector spaces is used as a means to quantify the performance of the approximation schemes.
In particular, we emphasize that the vector spaces are not necessarily complete with respect to these norms.

Let B, Y, k · k be another dual pair of normed vector spaces. As there is no danger of confusion, we use
the same notation for the potentially different norm and bilinear form for each pair. Let A : X → B be a
linear operator, and K be a convex cone in B. Given the fixed elements c ∈ C and b ∈ B, we define a linear
program, hereafter called the primal program P, as
(
inf
x, c
x∈X
J :=
(P)
s. t. Ax K b

where the conic inequality Ax K b is understood in the sense of Ax − b ∈ K. Throughout this study we
assume that the program P has an optimizer (i.e., the infimum is indeed a minimum), the cone K is closed
and the operator A is continuous where the corresponding topology is the weakest in which the topological
duals of X and B are C and Y, respectively. Let A∗ : Y → C be the adjoint operator of A defined by
Ax, y = x, A∗ y ,

∀x ∈ X,

∀y ∈ Y.

6

Recall that if A is weakly continuous, then the adjoint operator A∗ is well defined as its image is a subset of
C [26, Proposition 12.2.5]. The dual program of P is denoted by D and is given by


 sup b, y

y∈Y
Je :=
(D)
s.
t. A∗ y = c



y ∈ K∗ ,

where K∗ is the dual cone of K defined as K∗ := y ∈ Y : b, y ≥ 0, ∀b ∈ K . It is not hard to see that weak
duality holds, as
e
J = inf sup x, c − Ax − b, y ≥ sup inf x, c − Ax − b, y = J.
x∈X y∈K∗

y∈K∗ x∈X

An interesting question is when the above assertion holds as an equality. This is known as zero duality gap,
also referred to as strong duality particularly when both P and D admit an optimizer [1, p. 52]. Our study
is not directly concerned with conditions under which strong duality between P and D holds; see [1, Section
3.6] for a comprehensive discussion of such conditions. The programs P and D are assumed to be infinite, in
the sense that the dimensions of the decision spaces (X in P, and Y in D) as well as the number of constraints
are both infinite.

3.2. Semi-infinite approximation
Consider a family of linearly independent elements {xn }n∈N ⊂ X, and let Xn be the finite dimensional
subspace generated by the first n elements {xi }i≤n . Without loss of generality, we assume that xi are
normalized, i.e., kxi k = 1. Restricting the decision space X of P to Xn , along with an additional norm
constraint, yields the program

Pn
inf
αi xi , c


 α∈Rn Pi=1
n
s. t.
(4)
Jn :=
i=1 αi Axi K b



kαkR ≤ θP
where k · kR is a given norm on Rn and θP determines the size of the feasible set. In the spirit of dual-paired
normed vector spaces, one can approximate (X, C, k · k) by the finite dimensional counterpart (Rn , Rn , k · kR )
where the bilinear form is the standard inner product. In this view, the linear operator A : X → B may also
be approximated by the linear operator An : Rn → B with the respective adjoint A∗n : Y → Rn defined as
An α :=

n
X
i=1

αi Axi ,

A∗n y :=




Ax1 , y , · · · , Axn , y .

(5)

·

It is straightforward to verify the definitions (5) by noting that An α, y = α A∗n y for all α ∈ Rn and y ∈ Y.
Defining the vector c := [ x1 , c , · · · , xn , c ], we can rewrite the program (4) as


inf α c

 α∈Rn
Jn :=
(Pn )
s. t.
An α K b



kαk ≤ θ .

·

R

P

We call Pn a semi-infinite program, as the decision variable is a finite dimensional vector α ∈ Rn , but the
number of constraints is still in general infinite due to the conic inequality. The additional constraint on the
norm of α in Pn acts as a regularizer and is a key difference between the proposed approximation schemes
and existing schemes in the literature. Methods for choosing the parameter θP will be discussed later.
Dualizing the conic inequality constraint in Pn and using the dual norm definition leads to a dual counterpart

 sup b, y − θP kA∗n y − ckR∗
y∈Y
Jen :=
(Dn )
 s. t. y ∈ K∗ ,

7

wherek · kR∗ denotes the dual norm of k · kR . Note that setting θP = ∞ effectively implies that the second
term of the objective in Dn introduces n hard constraints A∗n y = c (cf. (5)). We study further the connection
between Pn and Dn under the following regularity assumption:
Assumption 3.1 (Semi-infinite regularity). We stipulate that
(i) the program Pn is feasible;

(ii) there exists a positive constant γ such that kA∗n ykR∗ ≥ γkyk∗ for every y ∈ K∗ , and θP is large
enough so that γθP > kbk.

Assumption 3.1(ii) is closely related to the condition
inf∗ sup

y∈K x∈Xn

Ax, y
≥ γ,
kxkkyk∗

that in the literature of numerical algorithms in infinite dimensional spaces, in particular the Galerkin discretization methods for partial differential equations, is often referred to as the “inf-sup” condition, see [20]
for a comprehensive survey. To see this, note that for every x ∈ Xn the definitions in (5) imply that
Ax, y = An α, y = α

·

A∗n y,

x=

n
X

αi xi .

i=1

These conditions are in fact equivalent if the norm k · kR is induced by the original norm on X, i.e., kαkR :=
Pn
k i=1 αi xi k. We note that A∗n maps an infinite dimensional space to a finite dimensional one, and as such
Assumption 3.1(ii) effectively necessitates that the null-space of A∗n intersects the positive cone K∗ only at 0.
In the following we show that this regularity condition leads to a zero duality gap between Pn and Dn , as well
as an upper bound for the dual optimizers. The latter turns out to be a critical quantity for the performance
bounds of this study.
Proposition 3.2 (Duality gap & bounded dual optimizers). Under Assumption 3.1(i), the duality gap between
the programs Pn and Dn is zero, i.e., Jn = Jen . If in addition Assumption 3.1(ii) holds, then for any optimizer
yn⋆ of the program Dn and any lower bound JnLB ≤ Jn we have
kyn⋆ k∗ ≤ θD :=

2θP kckR∗
θP kckR∗ − JnLB
≤
.
γθP − kbk
γθP − kbk

(6)

Proof. Since the elements {xi }i≤n are linearly independent, the feasible set of the decision variable α in
program Pn is a bounded closed subset of a finite dimensional space, and hence compact. Thus, thanks to
the feasibility Assumption 3.1(i) and compactness of the feasible set, the zero duality gap follows because
n
o
n
o
α c + sup b − An α, y = sup
Jn = inf
b, y − α (A∗n y − c) = Jen ,
inf
kαkR ≤θP

·

y∈K∗

y∈K∗ kαkR ≤θP

·

where the first equality holds by the definition of the dual cone K∗ , and the second equality follows from
Sion’s minimax theorem [44, Theorem 4.2]. Thanks to the zero duality gap above, we have
JnLB ≤ Jn = Jen = b, yn⋆ − θP kA∗n yn⋆ − ckR∗ ≤ b, yn⋆ − θP kA∗n yn⋆ kR∗ + θP kckR∗ .

By Assumption 3.1(ii), we then have


Jn ≤ kbkkyn⋆ k∗ − γθP kyn⋆ k∗ + θP kckR∗ = θP kckR∗ − γθP − kbk kyn⋆ k∗ ,

which together with the simple lower bound JnLB := −θP kckR∗ ≤ Jn concludes the proof.



Proposition 3.2 effectively implies that in the program Dn one can add a norm constraint kyk∗ ≤ θD
without changing the optimal value. The parameter θD depends on JnLB , a lower bound for the optimal value
of Jn . A simple choice for such a lower bound is −θP kckR∗ , but in particular problem instances one may be
able to obtain a less conservative bound. We validate the assertions of Proposition 3.2 for long-run average
cost problems in the next section and for long-run discounted cost problems in Appendix A.

8

Program Pn is a restricted version of the original program P (also called an inner approximation [26,
Definition 12.2.13]), and thus J ≤ Jn . However, under Assumption 3.1, we show that the gap Jn − J can
be quantified explicitly. To this end, we consider the projection mapping ΠA (x) := arg minx′ ∈A kx′ − xk, the
operator norm kAk := supkxk≤1 kAxk, and define the set
Bn :=

n
nX
i=1

o
αi xi ∈ Xn : kαkR ≤ θP .

(7)

Theorem 3.3 (Semi-infinite approximation). Let x⋆ and yn⋆ be optimizers for the programs P and Dn ,
respectively, and let rn := x⋆ − ΠBn (x⋆ ) be the projection residual of the optimizer x⋆ onto the set Bn as
defined in (7). Under Assumption 3.1(i), we have 0 ≤ Jn − J ≤ rn , A∗ yn⋆ − c where Jn and J are the
optimal value of the programs Pn and P. In addition, if Assumption 3.1(ii) holds, then

0 ≤ Jn − J ≤ kck∗ + θD kAk krn k,
(8)
where θD is the dual optimizer bound introduced in (6).

Proof. The lower bound 0 ≤ Jn − J is trivial, and we only need to prove the upper bound. Note that since
the optimizer x⋆ ∈ X is a feasible solution of P, then Ax⋆ − b ∈ K. By the definition of the dual cone K∗ , this
implies that Ax⋆ − b, y ≥ 0 for all y ∈ K∗ . Since the dual optimizer yn⋆ belongs to the dual cone K∗ , then
Jn − J ≤ Jn − J + Ax⋆ − b, yn⋆ = Jn − x⋆ , c + Ax⋆ , yn⋆ − b, yn⋆
= Jn + x⋆ , A∗ yn⋆ − c − b, yn⋆

n

= Jn + rn , A∗ yn⋆ − c + ΠBn (x⋆ ), A∗ yn⋆ − c − b, yn⋆ ,

A∗n yn⋆ − c − b, yn⋆ ,
= Jn + rn , A∗ yn⋆ − c + α
e

·

for some α
e ∈ R with norm ke
αkR ≤ θP ; for the last line, see the definition of the operator An in (5) as well
as the vector c in the program Pn . Using the definition of the dual norm and the operators (5), one can
deduce from above that
Jn − J ≤ Jn + rn , A∗ yn⋆ − c + θP kA∗n yn⋆ − ckR∗ − b, yn⋆ = Jn + rn , A∗ yn⋆ − c − Jen ,

which in conjunction with the zero duality gap (Jn = Jen ) establishes the first assertion of the proposition.
The second assertion is simply the consequence of the first part and the norm definitions, i.e.,


rn , A∗ yn⋆ − c = rn , −c + Arn , yn⋆ ≤ krn kkck∗ + kArn kkyn⋆ k∗ ≤ krn k kck∗ + kAkkyn⋆ k∗ .

Invoking the bound on the dual optimizer yn⋆ from Proposition 3.2 completes the proof.



Remark 3.4 (Impact of norms on semi-infinite approximation). We note the following concerning the impact
of the choice of norms on the approximation error:
(i) The only norm that influences the semi-infinite program Pn is k · kR on Rn . When it comes to the
approximation error (8), the norm k · kR may have an impact on the residual rn only if the set Bn in
(7) does not contain ΠXn (x⋆ ), the projection x⋆ on the subspace Xn , where x⋆ is an optimizer of the
infinite program P.

(ii) The norms of the dual pairs of vector spaces only appear in Theorem 3.3 to quantify the approximation error. Note that in (8) the stronger the norm on X, the higher krn k, and the lower kck∗ and
kAk. On the other hand, the stronger the norm on B, the higher kbk and kAk and the lower γ (cf.
Assumption 3.1(ii)).

The error bound (8) can be further improved when X is a Hilbert space. In this case, let Xn denote the
orthogonal complement of Xn . We define the restricted norms by
kck∗n := sup

x∈Xn

x, c
,
kxk

kAkn := sup

x∈Xn

kAxk
.
kxk

It is straightforward to see that by definition kck∗n ≤ kck∗ and kAkn ≤ kAk.

(9)

9

Corollary 3.5 (Hilbert structure). Suppose that X is a Hilbert space and k · k is the norm induced by the
corresponding inner product. Let {xi }i∈N be an orthonormal dense family and k · kR = k · kℓ2 . Let x⋆ be an
optimal solution for P and chose θP ≥ kx⋆ k. Under the assymptions of Theorem 3.3, we have

0 ≤ Jn − J ≤ kckn + θD kAkn ΠXn (x⋆ ) .

Proof. We first note that the ℓ2 -norm on Rn is indeed the norm induced by k · k, since due to the orthonormality of {xi }i∈N we have
v
u n
n
X
uX
αi xi = t
α2 kxi k2 = kαkℓ .
kαkR :=
i

i=1

2

i=1

If θP ≥ kx⋆ k, then ΠBn (x⋆ ) = ΠXn (x⋆ ), i.e., the projection of the optimizer x⋆ on the ball Bn is in fact
the projection onto the subspace Xn . Therefore, thanks to the orthonormality, the projection residual rn =
x⋆ − ΠXn (x⋆ ) belongs to the orthogonal complement Xn . Thus, following the same reasoning as in the proof
of Theorem 3.3, one arrives at a bound similar to (8) but using the restricted norms (9); recall that the norm
in a Hilbert space is self-dual.


3.3. Semi-infinite results in the MDP setting
We now return to the MDP setting in Section 2, and in particular the AC problem (2), to investigate the
application of the proposed approximation scheme. Recall that the AC problem (1) can be recast in an LP
framework in the form of P, see (3). To complete this transition to the dual pairs, we introduce the spaces


 X = R × L (S), C = R × M(S),
(10)
B = L (K),
Y = M(K),

 K = L (K),
∗
K
=
M
(K).
+
+

The bilinear form between each pair (X, C) and (B, Y) is defined in an obvious way (cf. (3)). The linear
operator A : X → B is defined as A(ρ, u)(s, a) := −ρ − u(s) + Qu(s, a), and it can be shown to be weakly
continuous [26, p. 220]. On the pair (X, C) we consider the norms



′
)
 kxk = k(ρ, u)k = max |ρ|, kukL } = max |ρ|, kuk∞ , sups,s′ ∈S u(s)−u(s
ks−s′ kℓ∞ ,
(11a)
R
 kck := sup
x,
c
=
|c
|
+
sup
u(s)c
(ds)
=
|c
|
+
kc
k
.
∗
1
2
1
2
W
kxk≤1
kukL ≤1 S

Recall that k · kL is the Lipschitz norm on L (S) whose dual norm k · kW in M(S) is known as the Wasserstein

norm [47, p. 105]. The adjoint operator A∗ : Y → C is given by A∗ y( · ) := − 1, y , −y( · × A) + yQ( · ) ,
where 1 is the constant function in L (S) with value 1. In the second pair (B, Y), we consider the norms


′
)
 kbk = kbkL := max kbk∞ , supk,k′ ∈K b(k)−b(k
kk−k′ kℓ∞ ,
(11b)
 kyk := sup
∗
kbkL ≤1 b, y = kykW .

A commonly used norm on the set of measures is the total variation whose dual (variational) characterization
is associated with k · k∞ in the space of continuous functions [26, p. 2]. We note that in the positive cone
K∗ = M+ (K) the total variation and Wasserstein norms indeed coincide.
Following the construction in Pn , we consider a collection of n-linearly independent, normalized functions
{ui }i≤n , kui kL = 1, and define the semi-infinite approximation of the AC problem (2) by

inf
−ρ



 (ρ,α)∈R×Rn
n

P
−JnAC =
(12)
αi ui (s) − Qui (s, a) ≤ ψ(s, a), ∀(s, a) ∈ K
s.
t.
ρ
+


i=1


kαkR ≤ θP

Comparing with the program Pn , we note that the finite dimensional subspace Xn ⊂ R×L (S) is the subspace
spanned by the basis elements x0 = (1, 0) and xi = (0, ui ) for all i ∈ {1, · · · , n}, i.e., the subspace Xn is

10

in fact n + 1 dimensional. Moreover, the norm constraint in (12) is only imposed on the second coordinate
of the decision variables (ρ, α) (i.e., kαkR ≤ θP ). The following lemmas address the operator norm and the
respective regularity requirements of Assumption 3.1 for the program (12).
Lemma 3.6 (MDP operator norm). In the AC problem (2) under Assumption 2.1(ii) with the specific norms
defined in (11), the linear operator norm satisfies kI − Qk := supkukL ≤1 ku − QukL ≤ 1 + max{LQ , 1}.
Proof. Using the triangle inequality it is straightforward to see that
kQukL
kQukL
ku − QukL
≤ 1 + sup
≤ 1 + sup
kuk
kuk
L
L
u∈L (S)
u∈L (S) kuk∞
u∈L (S)
n
o
kQuk∞
≤ 1 + max LQ , sup
≤ 1 + max{LQ , 1},
u∈L (S) kuk∞

kI − Qk =

sup

where the second line is an immediate consequence of Assumption 2.1(ii) and the fact that the operator Q is
R
R

a stochastic kernel. Hence, |Qu(s, a)| = | S u(y)Q(dy|s, a)| ≤ kuk∞ ( S Q(dy|s, a)) = kuk∞ .

Lemma 3.7 (MDP semi-infinite regularity). Consider the AC program (2) under Assumption 2.1. Then,
Assumption 3.1 holds for the semi-infinite counterpart in (12) for any positive θP and all sufficiently large γ.
In particular, the dual optimizer bound in Proposition 3.2 simplifies to kyn⋆ kW ≤ θD = 1.

Proof. Since K is compact, for any nonnegative θP , the program (12) is feasible and the optimal value is
bounded; recall that k(Q − I)ui kL ≤ 1 + max{LQ , 1} from Lemma 3.6 and kψk∞ < ∞ thanks to Assumption 2.1(iii). Hence, the optimal value of (12) is bounded and, without loss of generality, one can add a
redundant constraint |ρ| ≤ ω −1 θP , where ω is a sufficiently small positive constant. In this view, the last
constraint kαkR ≤ θP may be replaced with
k(ρ, α)kω := max{ω|ρ|, kαkR } ≤ θP ,

(13)

where k · kω can be cast as the norm on the pair (ρ, α) ∈ R × Rn+1 . Using the ω-norm as defined in (13), we
can now directly translate the program (12) into the semi-infinite framework of Pn . As mentioned above, the
feasibility requirement in Assumption 3.1(i) immediately holds. In addition, observe that for every y ∈ K∗
we have


− 1, y , Qu1 − u1 , y , · · · , Qun − un , y
sup (ρ, α)
kA∗n ykω∗ =
k(ρ,α)kω ≤1

·

= sup −ρ 1, y + sup α
ω|ρ|≤1

≥ω

−1

kαkR ≤1

kykW ,

·



Qu1 − u1 , y , · · · , Qun − un , y



where the third line above follows from the equality 1, y = kykW for every y in the positive cone K∗ ,
and the fact that the second term in the second line is nonnegative. Since ω can be arbitrarily close to 0,
the inf-sup requirement Assumption 3.1(ii) holds for all sufficiently large γ = ω −1 . The second assertion
of the lemma follows from the bound (6) in Proposition 3.2. To show this, recall that in the MDP setting
c = (−1, 0) ∈ R × M(S) (cf. (3)) with the respective vector c = [−1, 0, · · · , 0] ∈ R × Rn (cf. Pn ). Thus,
kckω∗ = supk(ρ,α)kω ≤1 (ρ, α) [−1, 0, · · · , 0] = ω −1 , that helps simplifying the bound (6) to

·

kyn⋆ kW ≤ θD :=

θP kckR∗ − JnLB
θP ω −1 + kψk∞
,
= −1
γθP − kbk
ω θP − kψkL

which delivers the desired assertion when ω tends to 0.



Remark 3.8 (AC dual optimizers bound). As opposed to the general LP in Proposition 3.2, Lemma 3.7
implies that the dual optimizers for the AC problem is not influenced by the primal norm bound θP and is
uniformly bounded by 1. In fact, this result can be strengthened to kyn⋆ kW = 1 due to the special minimax
structure of the AC program (12). This refinement is not needed at this stage and we postpone the discussion to

11

Section 5.2. The feature discussed in this remark does, however, not hold for the class of long-run discounted
cost problems, see Lemma A.2 in Appendix A.
Now we are in a position to translate Theorem 3.3 to the MDP setting for the AC problem (2).
Corollary 3.9 (MDP semi-infinite approximation). Let J AC and u⋆ be the optimal value and an optimizer
for the AC program (2), respectively. Consider the semi-infinite program (12) where θP > kψkL , and let
Pn
Un := { i=1 αi ui : kαkR ≤ θP }. Then, the optimal value of (12) satisfies the inequality

0 ≤ J AC − JnAC ≤ 1 + max{LQ , 1} u⋆ − ΠUn (u⋆ ) L .

Proof. We first note that the existence of the optimizer u⋆ is guaranteed under Assumption 2.1 [26, Theorem 12.4.2]. The proof is a direct application of Theorem 3.3 under the preliminary results in Lemma 3.7

and 3.6. Observe that the projection error is rn := (ρ⋆ , u⋆ ) − ΠUn (ρ⋆ , u⋆ ) = 0, u⋆ − ΠUn (u⋆ ) , resulting in
rn , c = 0. Thanks to this observation Lemma 3.6, the assertion of Theorem 3.3 translates to
0 ≤ J AC − JnAC = Jn − J ≤ rn , A∗ yn⋆ − c = Arn , yn⋆ ≤ kI − Qk krn kL kyn⋆ kW
≤ (1 + max{LQ , 1})ku⋆ − ΠUn (u⋆ )kL .



Observe that if from the beginning we consider the norm k · k∞ on the spaces X and B, it is not difficult
to see that the operator norm in Lemma 3.6 simplifies to 2 (recall that Q is a stochastic kernel). Thus, the
semi-infinite bound reduces to J AC −JnAC ≤ 2ku⋆ −ΠUn (u⋆ )k∞ . One may arrive at this particular observation
through a more straightforward approach: Using the shorthand notation (Q − I)u := Qu − u, we have






J AC − JnAC ≤ min Q − I u⋆ (k) + ψ(k) − min Q − I ΠUn (u⋆ )(k) + ψ(k)
k∈K
k∈K




≤ max Q − I u⋆ − ΠUn (u⋆ ) (k) ≤ Q − I u⋆ − ΠUn (u⋆ ) ∞
k∈K

≤ 2 u⋆ − ΠUn (u⋆ )

∞

.

Theorem 3.3 is a generalization to the above observation in two respects:
• It holds for a general LP that, unlike the AC problem (2), may not necessarily enjoy a min-max
structure.
• The result reflects how the bound on the decision space (i.e., θP in Pn ) influences the dual optimizers
as well as the approximation performance in generic normed spaces.
The latter feature is of particular interest as the boundedness of the decision space is often an a priori
requirement for optimization algorithms, see for instance [38] and the results in Section 5. The approximation
error from the original infinite LP to the semi-infinite version is quantified in terms of the projection residual
of the value function. Clearly, this is where the choice of the finite dimensional ball Un plays a crucial role.
We close this section with a remark on this point.
Remark 3.10 (Projection residual). The residual error u⋆ − ΠUn (u⋆ ) L can be approximated by leveraging
results from the literature on universal function approximation. Prior information about the value function
u⋆ may offer explicit quantitative bounds. For instance, for MDP under Assumption 2.1 we know that u⋆ is
Lipschitz continuous. For appropriate choice of basis functions, we can therefore ensure a convergence rate
of n−1/ dim(S) where dim(S) is the dimension of the state-action set S, see for instance [21] for polynomials
and [40] for the Fourier basis functions.
4. Semi-infinite to Finite Programs: Randomized Approach
We study conditions under which one can provide a finite approximation to the semi-infinite programs
of the form Pn , that are in general known to be computationally intractable — NP-hard [4, p. 16]. We
approach this goal by deploying tools from two areas, leading to different theoretical guarantees for the
proposed solutions. This section focuses on a randomized approach and the next section is dedicated to an

12

iterative gradient-based decent method. The solution of each of these methods comes with a priori as well as
a posteriori performance certificates.

4.1. Randomized approach
We start with a lemma suggesting a simple bound on the norm of the operator An in (5). We will use the
bound to quantify the approximation error of our proposed solutions.
Lemma 4.1 (Semi-infinite operator norm). Consider the operator An : Rn → B as defined in (5). Then,
kAn k := sup

α∈Rn

kAn αk
≤ kAk̺n ,
kαkR

̺n :=

sup kαkℓ1 ,

(14)

kαkR ≤1

where the constant ̺n is the equivalence ratio between the norms k · kR and k · kℓ1 .1
Proof. The proof follows directly from the definition of the operator norm, that is,
kAn αk =
together with the inequality

Pn

i=1

n
X
i=1

αi Axi ≤ kAk

n
X

αi xi ,

i=1

αi xi ≤ kαkℓ1 maxi≤n kxi k = kαkℓ1 , which concludes the proof.



Since K is a closed convex cone, then K∗∗ = K [1, p. 40], and as such the conic constraint in program Pn
can be reformulated as
An α K b

⇐⇒

An α − b, y ≥ 0,

∀y ∈ K := E{y ∈ K∗ : kyk∗ = 1},

(15)

where E{B} denotes the extreme points of the set B, i.e., the set of points that cannot be represented as
a strict convex combination of some other elements of the set. Notice that the norm constraint as well as
the restriction to the extreme points in the definition of K in (15) does not sacrifice any generality, as conic
constraints are homogeneous. These restrictions are introduced to improve the approximation errors. In what
follows, however, one can safely replace the set K with any subset of the cone K∗ whose closure contains K.
This adjustment may be taken into consideration for computational advantages. Let P be a Borel probability
measure supported on K, and {yj }j≤N be independent, identically distributed (i.i.d.) samples generated from
P. Consider the scenario counterpart of the program Pn defined as

minn α c

 α∈R
(Pn,N )
Jn,N :=
s. t.
α A∗n yj ≥ b, yj , j ∈ {1, · · · , N }


kαkR ≤ θP ,

·
·

where the adjoint operator A∗n : B → Rn is introduced in (5). The optimization problem Pn,N is a standard
finite convex program, and thus computationally tractable whenever the norm constraint kαkR ≤ θP is
tractable. Program Pn,N is a relaxation of Pn , i.e., Jn ≥ Jn,N ; note that Jn,N is a random variable, therefore
the relaxation error Jn − Jn,N can only be interpreted in a probabilistic sense.
Definition 4.2 (Tail bound). Given a probability measure P supported on K, we define the function p :
Rn × R+ → [0, 1] as
h
i
p(α, ζ) := P y : σK (−An α + b) < − An α + b, y + ζ ,

where σK ( · ) := supy∈K · , y is the support function of K. We call h : Rn × [0, 1] → R+ a tail bound (TB)
of the program Pn,N , if for all ε ∈ [0, 1] and α we have

h(α, ε) ≥ sup ζ : p(α, ζ) ≤ ε .
1The constant ̺ is indexed by n as it potentially depends on the dimension of α ∈ Rn .
n

13

The TB function in Definition 4.2 can be interpreted as a shifted quantile function of the mapping y 7→
− An α + b, y on K— the “shift” is referred to the maximum value of the mapping which is σK (−An α + b).
TB functions depend on the probability measure P generating the scenarios {yj }j≤N in the program Pn,N , as
well as the properties of the optimization problem. Definition 4.2 is rather abstract and not readily applicable.
The following example suggests a more explicit, but not necessarily optimal, candidate for a TB.
Example 4.3 (TB candidate). Let g : R+ → [0, 1] be a non-decreasing function such that for any κ ∈ K we


have g(γ) ≤ P Bγ (κ) , where Bγ (κ) is the open ball centered at κ with radius γ; note that function g depends
on the choice of the norm on Y. Then, a candidate for a TB function of the program Pn,N is

h(α, ε) := kAn α − bkg −1 (ε) ≤ ̺n kAkkαkR + kbk g −1 (ε),
where the inverse function is understood as g −1 (ε) := sup{γ ∈ R+ : g(γ) ≤ ε}, and ̺n is the constant ratio
defined in (14).
To see this note that according to Definition 4.2 we have
h
i
p(α, ζ) = P y : sup − An α + b, κ − y < ζ
κ∈K
h
i
= inf P y : − An α + b, κ − y < ζ
κ∈K
h
i
≥ inf P y : kAn α − bkky − κk∗ < ζ
κ∈K
h
i
= inf P Bγ(ζ) (κ) ≥ g(γ(ζ)), γ(ζ) := ζkAn α − bk−1 .
κ∈K

Thus, if p(α, ζ) ≤ ε, then g(γ(ζ)) ≤ ε and by construction of the inverse function g −1 we have ζkAn α−bk−1 ≤
g −1 (ε). In view of Definition 4.2, this observation readily suggests that the function h(α, ε) := kAn α−bkg −1(ε)
is indeed a TB candidate, and the suggested upper bound follows readily from Lemma 4.1.
Theorem 4.4 (Randomized approximation error). Consider the programs Pn and Pn,N with the associated
optimum values Jn and Jn,N , respectively. Let Assumption 3.1 hold, α⋆N be the optimizer of the program
Pn,N , and the function h be a TB as in Definition 4.2. Given ε, β in (0, 1), we define
n−1
o
n
X N 
εi (1 − ε)N −i ≤ β .
(16)
N(n, ε, β) := min N ∈ N :
i
i=0
For all positive parameters ε, β and N ≥ N(n, ε, β) we have



PN 0 ≤ Jn − Jn,N ≤ θD h α⋆N , ε ≥ 1 − β,

(17a)

where ̺n is the ratio constant defined in Lemma 4.1. We then have
h
i
PN 0 ≤ Jn − Jn,N ≤ ε ≥ 1 − β .

(17c)

where the constant θD is defined as in (6). In particular, suppose the function h is the TB candidate from
Example 4.3 with corresponding g function, and

−1

(17b)
N ≥ N n, g(zn ε), β ,
zn := θD θP ̺n kAk + kbk

Theorem 4.4 extends the result [36, Theorem 3.6] in two respects:
• The bounds (17) are described in terms of a generic norm and the corresponding dual optimizer
bound.

• Through the optimizer of Pn,N , the bounds involve an a posteriori element (cf. (17a) to (17c)).
Before proceeding with the proof, we first remark on the complexity of the a priori bound of Theorem 4.4,
its implications for an appropriate choice of θP , and its dependence on the dual pair norms.

14

Remark 4.5 (Curse of dimensionality). The TB function h of Example 4.3 may grow exponentially in the
dimension of the support set K (i.e., h(α, ε) ∝ ε− dim(K) ). Since N(n, · , β) admits a linear growth rate,
the a priori bound (17c) effectively leads to an exponential number of samples in the precision level ε, an
observation related to the curse of dimensionality [36, Remark 3.9]. To mitigate this inherent computational
complexity, one may resort to a more elegant sampling approach so that the required number of samples N
has a sublinear rate in the second argument, see for instance [37].
Remark 4.6 (Optimal choice of θP ). In view of the a priori error in Theorem 4.4, the parameter θP may
be chosen so as to minimize the required number of samples. To this end, it suffices to maximize zn defined
in (17b) over all θP > kbkγ −1, see Assumption 3.1(ii), where θD is defined in (6). One can show that the
optimal choice in this respect is analytically available as
s
 kbk
kbk
kbk  kbk
J LB 
⋆
θP :=
,
+
+
− n
γ
γ
̺n kAk
γ
kckR∗
where JnLB is a lower bound on the optimal value of Pn used in (6).

Remark 4.7 (Norm impact on finite approximation). Besides to what has already been highlighted in Remark 3.4, the choice of norms in the dual pairs of normed vector spaces also has an impact on the function
g −1 (ε). More specifically, the stronger the norm in the space B, the larger the balls in the dual space Y, and
thus the smaller the function g −1 .
To prove Theorem 4.4 we need a few preparatory results.
Lemma 4.8 (Perturbation function). Given δ ∈ B, consider the δ-perturbed program of Pn defined as

infn α c

 α∈R
Jn (δ) :=
(Pn (δ))
s. t.
An α K b − δ


kαkR ≤ θP .

·

Under Assumption 3.1, we then have Jn − Jn (δ) ≤ δ, yn⋆ , where yn⋆ is an optimizer of Dn .
Proof. For the proof we first introduce the dual program of Pn (δ):
(
sup b − δ, y − θP kA∗n y − ckR∗
e
y
Jn (δ) :=
s. t. y ∈ K∗ .

(Dn (δ))

We then have

Jn − Jn (δ) = Jen − Jn (δ) = b, yn⋆ − θP kA∗n yn⋆ − ckR∗ − Jn (δ)
= δ, yn⋆ + b − δ, yn⋆ − θP kA∗n yn⋆ − ckR∗ − Jn (δ)

≤ δ, yn⋆ + Jen (δ) − Jn (δ) ≤ δ, yn⋆ ,

where the first line follows from the strong duality (gap-free) between Pn and Dn by Proposition 3.2. The
third line is due to the fact that yn⋆ is a feasible solution of Dn (δ), and the last line follows from weak duality
between Pn (δ) and Dn (δ).

Lemma 4.9 (Perturbation error). Let α⋆N be an optimal solution of Pn,N and assume that δ ∈ B satisfies
the conic inequality An α⋆N K b − δ. Then, under Assumption 3.1, we have 0 ≤ Jn − Jn,N ≤ δ, yn⋆ .
Proof. The lower bound on Jn − Jn,N is trivial since Pn,N is a relaxation of Pn . For the upper bound
the requirement on δ in the program Pn (δ) implies that α⋆N is a feasible solution of Pn (δ). We then have
Jn,N ≥ Jn (δ), and thus 0 ≤ Jn − Jn,N ≤ Jn − Jn (δ). Applying Lemma 4.8 completes the proof.

The following fact follows readily from Definition 4.2.

15



Lemma 4.10 (TB lower bound). If α ∈ Rn satisfies P y : An α − b, y < 0 ≤ ε, then for any TB function
in the sense of Definition 4.2 we have σK (−An α + b) ≤ h(α, ε).
Proof. By the definition of the support function we can equivalently write


p(α, ζ) = P y : An α − b, y < ζ − σK (−An α + b) .

Now setting ζ = σK (−An α + b) in the above relation together with the assumption of Lemma 4.10 yields
p(α, ζ) ≤ ε, which in light of a TB in Definition 4.2 suggests that σK (−An α + b) ≤ h(α, ε).

We follow our discussion with a result from randomized optimization in a convex setting.
Theorem 4.11 (Finite-sample probabilistic feasibility [10, Theorem 1]). Assume that the program Pn,N
admits a unique minimizer α⋆N .2 If N ≥ N(n, ε, β) as defined in (16), then with confidence at least 1 − β


(across multi-scenarios {yj }j≤N ⊂ K) we have P y : An αN − b, y < 0 ≤ ε.
We are now in a position to prove Theorem 4.4.

Proof of Theorem 4.4. By definition of the support function we know that σK (δ) = σconv(K) (δ) where conv(K)
is the convex hull of K. Recall that by definition of the set K in (15), we also have y/kyk∗ ∈ conv(K) for any
y ∈ K∗ . Thus, for any δ ∈ B and y ∈ K∗ we have δ, y ≤ kyk∗ σK (δ). This leads to
0 ≤ Jn − Jn,N ≤ − An α⋆N + b, yn⋆ ≤ kyn⋆ k∗ σK (−An α⋆N + b)
where the second inequality is due to Lemma 4.9 as δ = −An α⋆N + b clearly satisfies the requirements. By
Lemma 4.10 and Theorem 4.11, we know that with probability at least 1 − β we have σK (−An αN + b) ≤
h(αN , ε), which in conjunction with the dual optimizer bound in Proposition 3.2 results in (17a). Now
using the TB candidate in Example 4.3 immediately leads to the first assertion of (17c). Recall that the
solution Pn,N obeys the norm bound kα⋆N kR ≤ θP . Thus, by employing the triangle inequality together with
Lemma 4.1 we arrive at the second assertion (17c).

Theorem 4.4 quantifies the approximation error between programs Pn and Pn,N probabilistically in terms
of the TB functions as introduced in Definition 4.2. The natural question is under what conditions can
the proposed bound be made arbitrarily small. This question is intimately related to the behavior of TB
functions. For the TB candidate proposed in Example 4.3, the question translates to when does the measure
of a ball Bγ (κ) ⊂ K have a lower bound g(γ) uniformly away from 0 with respect to the location of its center:
The answer to this question also depends on the properties of the norm on (B, Y, k · k). A positive answer
to this question requires that the set K can be covered by finitely many balls, indicating that K is indeed
compact with respect to the (dual) norm topology. In the next subsection we study this requirement in more
detail in the MDP setting.

4.2. Randomized results in the MDP setting
We return to the MDP setting and discuss the implication of Theorem 4.4 as the bridge from the semiinfinite program Pn to the finite counterpart Pn,N . Recall the dual pairs of vector spaces setting in (10) with
the assigned norms (11). To construct the finite program Pn,N , we need to sample from the set of extreme
points of P(K), i.e., the set of point measures
 
K := E P(K) = δ(s,a) : (s, a) ∈ K ,

where δ(s,a) denotes a point probability distribution at (s, a) ∈ K. In this view, in order to sample elements
from K it suffices to sample from the state-action feasible pairs (s, a) ∈ K.
2The uniqueness assumption may be relaxed at the expense of solving an auxiliary convex program, see [36, Section 3.3].

16

Corollary 4.12 (MDP finite randomized approximation error). Let {(sj , aj )}j≤N be N i.i.d. samples generated from the uniform distribution on K. Consider the program

inf
−ρ



 (ρ,α)∈Rn+1
n

P
AC
(18)
−Jn,N
=
u
(s
)
−
Qu
(s
,
a
)
≤ ψ(sj , aj ), ∀j ∈ {1, · · · , N }
α
s.
t.
ρ
+
i
j
i
j
j
i


i=1


kαkR ≤ θP .

where the basis functions {ui }i≤n introduced in (12) are normalized (i.e., kui kL = 1). Let LQ be the Lipschitz
constant from Assumption 2.1(ii), and define the constant
−1
,
zn := θP ̺n (max{LQ , 1} + 1) + kψkL

where ̺n is the ratio constant introduced in (14). Then, for all ε, β in (0, 1) and N ≥ N n + 1, (znε)dim(K) , β
defined in (16), we have
h
i
AC
PN 0 ≤ Jn,N
− JnAC ≤ ε ≥ 1 − β.
Proof. Let (ρ⋆N , α⋆N ) be the optimal solution for (18). Observe that in the MDP setting, Assumption 2.1(ii)
implies
kAn α⋆N − bk =

− ρ⋆N +

n
X
i=1

α⋆N (i) (Q − I)ui + ψ

L

≤ (max{LQ , 1} + 1)


≤ (max{LQ , 1} + 1)θP ̺n max kui kL + kψkL ,
i≤n

n
X
i=1

α⋆N (i) ui

L

+ k − ρ⋆N + ψkL
(19)

where the equality k − ρ⋆N + ψkL = kψkL leading to (19) follows from the fact that ψ and ρ⋆ are non-negative
(note that α = 0, ρ = 0 is a trivial feasible solution for (18)). In the second step, we propose a TB candidate
in the sense of Definition 4.2. Note that for any k, k ′ ∈ K, by the definition of the Wasserstein norm we have
kδ{k} − δ{k′ } kW = min{1, kk − k ′ k∞ }. Thus, generating samples uniformly from K leads to




P Bγ (κ) ≥ P Bγ (k) ≥ γ dim(K) ,
∀κ ∈ K, ∀k ∈ K,
(20)

where, with slight abuse of notation, the first ball Bγ (κ) is a subset of the infinite dimensional space Y with
respect to the dual norm k · kW , while the second ball Bγ (k) is a subset of the finite dimensional space K
whose respective norm is k · k∞ . The relation (20) readily suggests a function g : R+ → [0, 1] for Example 4.3,
which together with (19) and the fact that the basis functions are normalized, it yields

h(α, ε) := kAn α − bkg −1 (ε) ≤ θP ̺n (max{LQ , 1} + 1) + kψkL ε1/ dim K .

Recall from Lemma 3.7 that the dual multiplier bound is θD = 1, and feasible solutions α is bounded by θP .
Finally, note that the decision variable of the program (18) is the n + 1 dimensional pair (ρ, α). Given all the
information above, the claim then readily follows from the second result of Theorem 4.4 in (17c).

To select θP , one may minimize the complexity of the a priori bound in Corollary 4.12, which is reflected
through the required number of samples. At the same time, the impact of the bound θP on the approximation
step from infinite to semi-infinite in Corollary 3.9 should also be taken into account. The first factor is
monotonically decreasing with respect to θP , i.e., the smaller the parameter θP , the lower the number of
the required samples. The second factor is presented through the projection residual (cf. Remark 3.10).
Therefore, an acceptable choice of θP is an upper bound for the projection error of the optimal solution onto
the ball Un uniformly in n ∈ N, i.e.,


n
X
α⋆i ui , n ∈ N .
(21a)
θP ≥ sup kα⋆ kR : ΠUn (x⋆ ) =
i=1

17

The above bound may be available in particular cases, e.g., when k · kR = k · kℓ2 it yields the bound
sZ
⋆
u⋆ 2 (s)ds ≤ ku⋆ kL ≤ max{LQ , 1}kψk∞,
kα kℓ2 =

(21b)

S

where LQ is the Lipschitz constant in Assumptions 2.1(ii). We note that the first inequality in (21b) follows
since S is a unit hypercube, and the second inequality follows from [19, Lemma 2.3], see also [19, Section 5]
for further detailed analysis.
5. Semi-infinite to Finite Program: Structural convex optimization
This section approaches the approximation of the semi-infinite program Pn from an alternative perspective
relying on an iterative first order decent method. As opposed to the scenario approach presented in Section 4,
that is probabilistic and starts from the program Pn , the method of this section is deterministic and starts with
the dual counterpart Dn , in particular a regularized version of whose solutions can be computed efficiently. It
turns out that the regularized solution allows one to reconstruct a nearly feasible solution for both programs
Pn and Dn , offering a meaningful performance bound for the approximation step from the semi-infinite
program to a finite program.

5.1. Structural convex optimization
The basis of our approach is the fast gradient method that significantly improves the theoretical and, in
many cases, also the practical convergence speed of the gradient method. The main idea is based on a well
known technique of smoothing nonsmooth functions [39]. To simplify the notation, for a given θP we define
the sets



Y := y ∈ K∗ : kyk∗ ≤ θD ,
A := α ∈ Rn : kαkR ≤ θP ,

where θD is the constant defined in (6). Recall that in the wake of Proposition 3.2 we know that the decision
variables of the dual program Dn may be restricted to the set Y without loss of generality. We modify the
program Dn with a regularization term scaled with the non-negative parameter η and define the regularized
program
n
o
(Dn,η )
Jen,η := sup b, y − θP kA∗n y − ckR∗ − ηd(y) ,
y∈Y

where the regularization function d : Y → R+ , also known as the prox-function, is strongly convex. The
choice of the prox-function depends on the specific problem structure and may have significant impact on the
approximation errors. Given the regularization term η and the parameter α ∈ Rn , we introduce the auxiliary
quantity
n
o
yη⋆ (α) := arg max b − An α, y − ηd(y) .
(22)
y∈Y

It is computationally crucial for the solution method proposed in this part that the prox-function allows us
to have access to the auxiliary variable yη⋆ (α) for each α ∈ Rn . This requirement is formalized as follows.

Assumption 5.1 (Lipschitz gradient). Consider the adjoint operator A∗n in (5) and the optimizer yη⋆ (α) of
the auxiliary quantity (22). We assume that for each α ∈ A the vector A∗n yη⋆ (α) ∈ Rn can be approximated
to an arbitrary precision, and the mapping α 7→ A∗n yη⋆ (α) is Lipschitz continuous with a constant Lη , i.e.,
kA∗n yη⋆ (α) − A∗n yη⋆ (α′ )kR∗ ≤

L
kα − α′ kR ,
η

∀α, α′ ∈ A .

Let ϑ > 0 be the strong convexity parameter of the mapping α 7→ 12 kαk2R with respect to the R-norm. We
then define the operator T : Rn × Rn → Rn as
o
n
1
(23)
kβ − αk2R ,
T(q, α) := arg min q β +
β∈A
2ϑ

·

18

Algorithm 1 Optimal scheme for smooth convex optimization
Choose some w(0) ∈ A
For k ≥ 0 do

Step 1: Define r(k) :=
Step 2: Compute z (k)

Step 3: Set w(k+1) =


c − A∗n yη⋆ (w(k) ) ;
Pk j+1 (j) 
:= T
,0 ,
j=0 2 r
η
L

2
(k)
k+3 z

+

α(k) := T

1 (k)
, w(k)
ϑr

k+1 (k)
.
k+3 α


;

More generally, a different norm can be used in the second term in (23) when ϑ is a different strong convexity
parameter. However, we forgo this additional generality to keep the exposition simple. The operator T is
defined implicitly through a finite convex optimization program whose computational complexity may depend
on the R-norm through the constraint set A . For typical norms in Rn (e.g., k · kℓp ) the pointwise evaluation
of the operator T is computationally tractable. Furthermore, if k · kR = k · kℓ2 , then the definition of (23) has
an explicit analytical description for any pair (q, α) as follows.
Lemma 5.2 (Explicit description of T). Suppose in the definition of the operator (23) the R-norm is the
classical ℓ2 -norm. Then, the operator T admits the analytical description T(q, α) = ξ (α − q) where ξ :=

min 1, θP kq − αk−1
ℓ2 .

Proof. In case of the ℓ2 -norm the strong convexity parameter is ϑ = 1. Now using the classical duality theory,
the objective function of (23) is equal to
n
o
o
n

1
1
2
+ min q β + kβ − αk2ℓ2 + λkβk2ℓ2
q β + kβ − αk2ℓ2 = max − λθP
min
(24)
λ≥0
β
2
2
kβkℓ2 ≤θP
n
λkα − qk2ℓ2 o
kq + 2λαk2ℓ2
q α − kqk2ℓ2
2
+
+
,
= max − λθP
+
λ≥0
1 + 2λ
2(1 + 2λ)2
(1 + 2λ)2

·

·

·

where the second equality in (24) follows by substituting the explicit solution of the unconstrained inner
problem described by

α−q
1
β ⋆ (λ) := arg min q β + kβ − αk2ℓ2 + λkβk2ℓ2 =
.
(25)
β
2
1 + 2λ

·

To find the optimal λ in the right-hand side of (24), it suffices to set the derivative to zero, which yields
λ⋆ := 2θ1P max{kα−qk−θP , 0}. By substituting λ⋆ in (25), we have an optimal solution β ⋆ (λ⋆ ) = ξ(α−q) that
is feasible since kβ ⋆ (λ⋆ )kℓ2 ≤ θP . By virtue of the equality in (24), this concludes the desired assertion. 

Algorithm 1 exploits the information revealed under Assumption 5.1 as well as the operator T to approximate the solution of the program Dn . The following proposition provides explicit error bounds for the
solution provided by Algorithm 1 after k iterations. The result is a slight extension of the classical smoothing
technique in finite dimensional convex optimization [39, Theorem 3] where the prox-function is not necessarily uniformly bounded, a potential difficulty in infinite dimensional spaces. We address this difficulty by
considering a growth rate for the prox-function d evaluated at the optimal solution yη⋆ . We later show how
this extension will help in the MDP setting.
Theorem 5.3 (Smoothing approximation error). Suppose Assumption 5.1 holds with constant L and ϑ is the
strong convexity parameter in the definition of the operator T in (23). Given the regularization term η > 0
and k iterations of Algorithm 1, we define
α
bη := α(k) ,

ybη :=

k
X
j=0

2(j + 1)
y ⋆ (w(j) ).
(k + 1)(k + 2) η

LB
UB
≤ Jn ≤ Jn,η
where
Under Assumption 3.1, the optimal value of the program Pn is bounded by Jn,η
LB
:= b, ybη − θP kA∗n ybη − ckR∗ ,
Jn,η

UB
:= α
bη
Jn,η

sup
· c + y∈Y

b − An α
bη , y

(26)

19

Moreover, suppose there exist positive constants c, C such that



C max log cη −1 , 1 ≥ d yη⋆ (α) ,

∀η > 0,

∀α ∈ A ,

and, given an a priori precision ε > 0, the regularization parameter η and the number of iterations k satisfy
p
CL max{2 log(2cCε−1 ), 1}
ε
√
η≤
,
k
≥
2θ
̺
,
(27)
P
n
2C max{2 log(2cCε−1 ), 1}
ϑε
UB
LB
− Jn,η
≤ ε.
where ̺n is the constant defined in (14). Then, after k iterations of Algorithm 1 we have Jn,η
LB
UB
Proof. Observe that the bounds Jn,η
and Jn,η
in (26) are the values of the programs Dn and Pn evaluated at
ybη and α
bη , respectively. As such, the first assertion follows immediately. Towards the second part, thanks to
the compactness of the set A , the strong duality argument of Sion’s minimax theorem [44] allows to describe
the program Dn,η through
i
h
Jen,η := sup b, y − sup An α, y − α c + ηd(y)
α∈A

y∈Y

h

·

i

sup b − An α, y − ηd(y)
· c + y∈Y
α∈A

= inf α · c + b − An α, yη⋆ (α) − ηd yη⋆ (α) ,
α∈A

= inf α

(28)

where the last equality follows from the definition in (22). Note that the problem (28) belongs to the class of
smooth and strongly convex optimization problems, and can be solved using a fast gradient method developed
by [39]. For this purpose, we define the function

φη (α) := α c + b − An α, yη⋆ (α) − ηd yη⋆ (α) .
(29)

·

Invoking similar techniques to [39, Theorem 1], it can be shown that the mapping α 7→ φη (α) is smooth with
the gradient ∇φη (α) = c − A∗n yη⋆ (α). The gradient ∇φη (α) is Lipschitz continuous by Assumption 5.1 with
constant Lη . Thus, following similar arguments as in the proof of [39, Theorem 3] we have
UB
LB
0 ≤ Jn,η
− Jn,η
≤

 L(θP ̺n )2

Lkα⋆ k2R
+ ηd yη⋆ (α∗ ) ≤
+ Cη max log(cη −1 ), 1 .
2
ϑ(k + 1)(k + 2)η
ϑk η

(30)

Now, it is enough to bound each of the terms in the right-hand side of the above inequality by 12 ε. It should
be noted that this may not lead to an optimal choice of the parameter η, but it is good enough to achieve a
reasonable precision order with respect to ε. To ensure η log(η −1 ) ≤ ε for an ε ∈ (0, 1) , it is not difficult to see
ε
1
1
that it suffices to set η ≤ 2 log(ε
−1 ) . In this observation if we replace η and ε with c η and 2cC ε, respectively,
we deduce that the second term on the right-hand side in (30) bounded by 12 ε. Thus, the desired assertion
follows by equating the first term on the right-hand side in (30) to 12 ε while the parameter η is set as just
suggested.

Remark 5.4 (Computational complexity). Adding the prox-function to the problem Dn ensures thatq
the regu
L −1
larized counterpart Dn,η admits an efficiency estimate (in terms of iteration numbers) of the order O
.
ηε
To construct a smooth ε-approximation for the original problem Dn , the Lipschitz constant Lη can be chosen
of the p
order O(ε−1 log(ε−1 )). Thus, the presented gradient scheme has an efficiency estimate of the order

−1
log(ε−1 ) , see [39] for a more detailed discussion along similar objective.
O ε

Remark 5.5 (Inexact gradient). The error bounds in Theorem 5.3 are introduced based on the availability
of the exact first-order information, i.e., it is assumed that at each iteration the vector r(k) that due to the
bilinear form potentially involves a multi dimensional integration can be computed exactly. In general, the
evaluation of those vectors may only be available approximately. This gives rise to the question of how the
fast gradient method performs in the case of inexact first-order information. We refer the interested reader
to [15] for further details.

20

The a priori bound proposed by Theorem 5.3 involves the positive constants c, C, which are used to
introduce an upper bound for the proxy-term. These constants potentially depend on θD , the size of the dual
feasible set, hence also on θP . Therefore, unlike the randomized approach in Section 4, it is not immediately
clear how θP can be chosen to minimize the complexity of the proposed method, which in this case is the
required number of iterations k suggested in (27) (cf. Remark 4.6). In the next section, we shall discuss how
to address this issue in the MDP setting for particular constants c, C.

5.2. Structural convex optimization results in the MDP setting
To link the approximation method presented in Section 5.1 to the AC program in (12), let us recall the
dual pairs (10) equipped with the norms (11). To simplify the analysis, we refine the assertion in Lemma 3.7
and argue that the dual optimizers are indeed probability measures, i.e.,


Y := y ∈ M+ (K) : kykW = θD = 1 .
(31)
To see this, one can consider the norm k(ρ, α)k := kαkR and follow similar arguments in the proof of
Proposition 3.2. Strictly speaking, this is not a true norm on Rn+1 but it does not affect the technical
argument, in particular strong duality between Pn and Dn . The details are omitted here in the interest of
space. We consider the prox-function as a relative entropy defined by
(

dy
,y y ≪ λ
log dλ
:=
(32)
d(y)
∞
o.w.,
dy
where λ is the uniform measure supported on the set K and dλ
∈ F+ (K) is the Radon-Nikodym derivative
between two measures y and λ. One can inspect that the prox-function (32) is indeed a non-negative function.
The optimizer of the regularized program Dn,η for the AC program (12) is


n
X

dy
αi (Q − I)ui , y − η log dλ
yη⋆ (ρ, α) := arg max
(33)
−ψ+ρ−
,y .
y∈Y

i=1

To see (33), check (22) together with the definitions of the operator An in (5) and the AC problem parameters
in (3). The main reason for such a choice of the regularization term is the fact that the optimizer of the
regularized program (33) admits an analytical expression:
Lemma 5.6 (Entropy maximization [12]). Given a (measurable) function g : K → R and the set Y ⊂
M+ (K) as defined in (31) we have
o exp η −1 g(k)λ(dk)
n
⋆

.
y (dk) := arg max g, y − ηd(y) =
y∈Y
exp η −1 g(k) , λ
Thanks to Lemma 5.6, the analytical description of the dual optimizer in (33) is readily available by setting
g(k) := [b − An α](k) = −ψ(k) + ρ −

n
X
i=1

αi (Q − I)ui (k).

(34)

The last requirement to implement Algorithm 1 is to verify Assumption 5.1, i.e., we need to compute the
Lipschitz constant of the mapping (ρ, α) 7→ A∗n yη⋆ (ρ, α) in which the respective norm is k(ρ, α)k := kαkR . By
definition of the adjoint operator A∗n in (5), it is not difficult to observe that

 

− 1, yη⋆ (ρ, α)
−1

 

 (Q − I)u1 , yη⋆ (ρ, α)   (Q − I)u1 , yη⋆ (ρ, α) 
∗ ⋆



.
An yη (ρ, α) = 
(35)
..
..
=

.
.

 

(Q − I)un , yη⋆ (ρ, α)

(Q − I)un , yη⋆ (ρ, α)

The next lemma addresses the requirement of Assumption 5.1 for the mapping (35).

21

Lemma 5.7 (Lipschitz constant in MDP). Consider the entropy maximizers in Lemma 5.6 with g as defined
in (34) and the adjoint operator in (35). An upper bound for the Lipschitz constant in Assumption 5.1 is
L ≤ 4̺2n where the constant ̺n is the equivalence ratio between the norms k · kR and k · kℓ1 introduced in (14).
Proof. It is straightforward to see that (35) is differentiable with respect to the variable (ρ, α). Hence, it
suffices to bound the norm of the matrix ∇A∗n yη⋆ (ρ, α) ∈ R(n+1)×(n+1) uniformly on (ρ, α). Further, as the
first element of the vector (35) is the constant 1, it only requires ro consider the gradient function with respect
to the variable α ∈ Rn . A direct computation yields

1
1
| ∇α A∗n yη⋆ (ρ, α) ij | =
(Q − I)ui (Q − I)uj , yη⋆ (ρ, α) − (Q − I)ui , yη⋆ (ρ, α) (Q − I)uj , yη⋆ (ρ, α)
η
η
4 (Q − I)ui (Q − I)uj ⋆
(Q − I)ui ⋆
(Q − I)uj ⋆
=
, yη (ρ, α) −
, yη (ρ, α)
, yη (ρ, α)
η
2
2
2
2
4
∀i, j ∈ {1, · · · , n}.
≤ ,
η
where the last inequality is a consequence of the Cauchy-Schwarz inequality and the fact that k(Q−I)uj k∞ ≤ 2
(recall that Q is a stochastic kernel and all the basis functions are normalized). The Lipschitz constant of the
desired mapping can then be upper bounded by
L
≤
η

sup
kαkR ≤ θP
kvkR ≤ 1

∇α A∗n yη⋆ (ρ, α)v

R∗

≤

sup
kΦi kℓ∞ ≤ 1
kvkR ≤ 1

4
(Φ1
η

· v, · · · , Φn · v)

R∗

(36)

·

Recall that by the definition of the dual norm, we have |Φi v| ≤ kΦi kR∗ for all kvkR ≤ 1. Thus, substituting
the scalar variable µi := Φi v in the right-hand side of (36) and eliminating the factor η lead to
L≤

·

sup

sup

kΦi kℓ∞ ≤1 |µi |≤kΦi kR∗

4 (µ1 , · · · , µn )

R∗

= sup 4 (µ1 , · · · , µn )
|µi |≤̺n

R∗

,

where the last statement follows from the definition of the dual norm, and in particular the equality
sup
kΦi kℓ∞ ≤1

kΦi kR∗ =

sup
kΦi kR ≤1

kΦi kℓ1 =: ̺n .

Thus, using the same equality yields
L≤

sup
kµkℓ∞ ≤1

4̺n kµkR∗ = 4̺2n ,

which concludes the first desired assertion.



The performance of Algorithm 1 can now be characterized through the following corollary.
Corollary 5.8 (MDP smoothing approximation error). Consider the operator (5) with the parameters described in (3) for the semi-infinite AC program (12). Given this setting and the Lipschitz constant in
Lemma 5.7, we run Algorithm 1 for k iterations using the entropy function (32) with analytical solution
(33) as the prox-function. We define the constants
r

2 dim(K)
2
C2 := 4θP ̺n
C1 := 2 e ̺n θP (max{LQ , 1} + 1) + kψkL ,
.
ϑ
For every ε ≤ C1 we set the smoothing factor η and the number of iterations k by
p
log(C1 ε−1 )
ε
,
k
≥
C
.
η≤
2
4 dim(K) log(C1 ε−1 )
ε
Then, the outcome of Algorithm 1 as defined in (26) is an ε approximation of the optimal value JnAC in the
sense of Theorem 5.3.
Corollary 5.8 requires one to compute the constants c, C to quantify the a priori bounds. The following
two technical lemmas provide supplementary materials to address this issue.

22

Lemma 5.9. Let K ⊆ [0, 1]m and g : K → R be a Lipschitz continuous function with constant Lg > 0 (with
respect to the ℓ∞ -norm) and the maximum value gmax := maxk∈K g(k). Then, for every η > 0 we have
Z
n mη m o




, 1 exp − min m, Lg η −1 .
exp η −1 g(k) − gmax dk ≥ min
Lg
K

Proof. Let us define the set Z(δ) := {k ∈ K : gmax − g(k) < δ}. Thanks to the Lipschitz continuity of the
function g, we have gmax − g(k) ≤ Lg kk ⋆ − kkℓ∞ where g(k ⋆ ) = gmax . Thus, using this inequality one can
bound the size of the set Z(δ) in the sense of
Z
m
dk ≥ min{(δL−1
∀δ ≥ 0.
g ) , 1},
Z(δ)

By virtue of the above result, one can observe that for every δ > 0
Z
Z




exp η −1 g(k) − gmax dk
exp η −1 g(k) − gmax dk ≥
Z(δ)

K

≥ exp(−η

−1

δ)

Z

Z(δ)

m
dk ≥ exp(−η −1 δ) min{(δL−1
g ) , 1}.

Maximizing the right-hand side of the above inequality over δ suggests to set δ = min{mη, Lg }, which yields
the desired assertion.

In light of Lemma 5.9, we can bound the entropy prox-function (32) evaluated at the optimizer (33).
Lemma 5.10 (Entropy prox-bound). Consider the prox-function (32) and let yη⋆ (ρ, α) be the optimizer of


(33). Then, for every η > 0, ρ, and kαkR ≤ θP , we have d yη⋆ (ρ, α) ≤ C max log(cη −1 ), 1 where

e
θP ̺n (max{LQ , 1} + 1) + kψkL ,
C := dim(K),
c :=
dim(K)

and ̺n is the equivalence ratio between the norms k · kℓ1 and k · kR as defined in (14).

Proof. The result is a direct application of Lemma 5.9. Consider the function g as defined in (34) with
Lipschitz constant Lg ≥ 0; note that the function g, as well as its Lipschitz constant Lg , depends also on the
pair (ρ, α). Observe that



d yη⋆ (ρ, α) = log exp(η −1 g) , yη⋆ (ρ, α) − log exp(η −1 g), λ

= η −1 g, yη⋆ (ρ, α) − log exp(η −1 g), λ

= η −1 g, yη⋆ (ρ, α) − η −1 gmax − log exp(η −1 (g − gmax ), λ

≤ − log exp(η −1 (g − gmax )), λ


n dim(K)η dim(K) o


≤ − log min
(37)
, 1 exp − min dim(K), Lg η −1
Lg
 o
n
 e L 
g
η −1 , 1
≤ dim(K) max log
dim(K)
where the inequality (37) follows from Lemma 5.9. Note also that the Lipschitz constant Lg for the function
g defined in (34) is upper bounded, uniformly in (ρ, α) where kαkR ≤ θP , by
Lg ≤ kg − ρkL ≤

n
X
i=1

αi (Q − I)ui + ψ

L

≤ (max{LQ , 1} + 1)

≤ θP ̺n (max{LQ , 1} + 1) + kψkL .

n
X
i=1

αi ui

L

+ kψkL

We refer to the proof of Corollary 4.12, and in particular the paragraph following (19), for further discussions
regarding Lg . The desired assertion follows from the last two inequalities and the definition of the constant
θD in (31).


23

The proof of Corollary 5.8 follows by replacing the constants in Lemma 5.10 in Theorem 5.3. By contrast
to the randomized approach in Corollary 4.12 where the computational complexity scales exponentially in
dimensional of state-action
space, the complexity of the smoothing technique grows effectively linearly (more
p

precisely O ε−1 log(ε−1 ) , cf. Remark 5.4). The computational difficulty is, however, transferred to Step 1
of Algorithm 1 for computation of A∗n yη⋆ as defined in (35). The following remark elaborates this.
Remark 5.11 (Efficient computation of (35)). When the transition kernel Q and the basis functions ui are
such that the relation (35) involves integration of exponentials of polynomials over simple sets (e.g., box or
a simplex), one may utilize efficient methods that require solving a hierarchy of semidefinite programming
problems to generate upper and lower bounds which asymptotically converge to the true value of integral,
see [32, Section 12.2] and [9]. It is also worth noting that a straightforward computation of (35) for a
small parameter η may be numerically difficult due to the exponential functions. This issue can, however, be
circumvented by a numerically stable technique presented in [39, p. 148].
Regarding the choice of θP , in similar spirit to Section 4, one can target minimizing the complexity of the
a priori bound, in other words the number of iterations k in (27). In the setting of Corollary 5.8, one can
observe that the smaller the parameter θP , the lower the number of the required iterations, leading to the
choice described as in (21).
6. Full Infinite to Finite Programs
The intention in this short section is to combine the two-step process from infinite to semi-infinite programs
in Section 3 and from semi-infinite to finite programs in Section 4 and 5, and hence establish a link from
the original infinite program to finite counterparts. We only present the final result for the general infinite
programs without discussing its implication in the MDP setting, as it is essentially a similar assertion.
Theorem 6.1 (Infinite to finite approximation error). Consider the infinite program P with a solution

{x⋆ , J}, the finite (random) convex program Pn,N with the (random) solution α⋆N , Jn,N , and the output
 LB UB
of Algorithm 1 with values Jn,η , Jn,η . Suppose Assumption 3.1 holds and assume further that there exists
constant d, D so that the projection residual of the optimizer x⋆ onto the finite dimensional ball defined in
Theorem 3.3 is bounded by krn k ≤ Dn−1/d for all n ∈ N. Then, for any number of scenario samples N and
prox-term coefficient η, with probability 1 − β we have


 UB
LB
max Jn,N , Jn,η
− D kck∗ + θD kAk n−1/d ≤ J ≤ min Jn,η
, Jn,N + θD h(α⋆N , ε) .

where θD is as defined in (6) and the function h is a TB in the sense of Definition 4.2. Moreover, given an
a priori precision level ε, if
d


n ≥ D kck∗ + θD kAk ε−1 ,

and the number samples N are chosen as in (17b) or the parameter η together with the number of iterations
of Algorithm 1 is chosen as in (27), then with probability 1 − β we have
n
o
LB
min |J − Jn,N |, |J − Jn,η
| ≤ ε.

The proof follows readily from the link between the infinite program P to the semi-infinite counterpart Pn
in Theorem 3.3, in conjunction with the link between Pn to the finite programs Pn,N and Dn,η in Theorems 4.4
and 5.3, respectively.

The assertion of Theorem 6.1 can be readily translated into the MDP problem by replacing the dual
optimizer bound θD with 1 thanks to Lemma 3.7, and the term (kck∗ + θD kAk) with (1 + max{LQ , 1}) thanks
to Corollary 3.9. In this case, the requirement concerning the projection residual bound krn k ≤ Dn−1/d is
fulfilled due to the Lipschitz continuity of the value function when d = dim(S) and the finite dimensional
approximation is generated by, among others, polynomials [21] or the Fourier basis [40] (cf. Remark 3.10).

24

7. Numerical Examples
We present two numerical examples to illustrate the solution methods and corresponding performance
√
bounds. Throughout this section we consider the norm k · kR = k · kℓ2 , leading to ̺n = n in (14), and we
choose the Fourier basis functions.

7.1. Example 1: truncated LQG
Consider the linear system
st+1 = ϑst + ρat + ξt ,

t ∈ N,

with quadratic stage cost ψ(s, a) = qs2 + ra2 , where q ≥ 0 and r > 0 are given constants. We assume
that S = A = [−L, L] and the parameters ϑ, ρ ∈ R are known. The disturbances {ξt }t∈N are i.i.d. random
variables generated by a truncated normal distribution with known parameters µ and σ, independent of the
initial state s0 . Thus, the process ξt has a distribution density

s−µ
1

σ φ( σ )
, s ∈ [−L, L]
L−µ
Φ( σ )−Φ( −L−µ
)
f (s, µ, σ, L) =
σ

0
o.w.,

where φ is the probability density function of the standard normal distribution, and Φ is its cumulative
R
distribution function. The transition kernel Q has a density function q(y|s, a), i.e., Q(B|s, a) = B q(y|s, a)dy
for all B ∈ B(S), that is given by
q(y|s, a) = f (y − ϑs − ρa, µ, σ, L).

In the special case that L = +∞ the above problem represents the classical LQG problem, whose solution can
be obtained via the algebraic Riccati equation [6, p. 372]. By a simple change of coordinates it can be seen
that the presented system fulfills Assumption 2.1. The following lemma provides the technical parameters
required for the proposed error bounds.
Lemma 7.1 (Truncated LQG properties).pThe error bounds provided by Corollaries 4.12 and 5.8 hold with
the norms kψk∞ = L2 (q + r), kψkL = 4L2 q 2 + r2 , and the Lipschitz constant of the kernel is
LQ =

2L max{ϑ, ρ}

 .
√   L−µ 
− Φ −L−µ
σ 2 2π Φ σ
σ

st
Proof. In regard to Assumption 2.1(i), we consider the change of coordinates s̄t := 2L
+ 21 and āt :=
In the new coordinates, the constants of Lemma 7.1 follow from a standard computation.

at
2L

+ 12 .


Simulation details: For the simulation results we choose the numerical values ϑ = 0.8, ρ = 0.5, σ = 1, µ = 0,
q = 1, r = 0.5, and L = 10. In the first approximation step discussed in Section 3.3, we consider the Fourier


L
L
basis u2k−1 (s) = kπ
cos kπs
sin kπs
and u2k (s) = kπ
L
L .
Randomized approach:

We implement the methodology presented in Section 4.2, resulting in a finite random convex program as
in (18), where the uniform distribution on K = S × A = [−L, L]2 is used to draw the random samples.
Figures 2(a), 2(b), and 2(c) visualize three cases with different number of basis functions n ∈ {2, 10, 100},
respectively. To show the impact of the additional norm constraint, in each case two approximation settings
are examined: the constrained (regularized) one proposed in this article (i.e., θP < ∞), and the unconstrained
one (i.e., θP = ∞). In the former we choose the bound suggested by (21b). In the latter, the resulting
optimization programs of (18) may happen to be unbounded, particularly when the number of samples N is
low; numerically, we capture the behavior of the unbounded θP through a large bound such as θ = 106 . In
each sub-figure, the colored tubes represent the results of 400 independent experiments (shaded areas) as well

25

106

106
θP = kψk∞
θP = ∞

θP = kψk∞
θP = ∞

J AC

J AC

AC
Jn,N

104

AC
Jn,N

104

102

102

100
100

100
100

101

102

103

104

105

101

N

102

103

104

105

N

(a) n = 2 basis functions

(b) n = 10 basis functions

106

100.3

θP = kψk∞
θP = ∞

n=2
n = 10
n = 100

J AC

104

AC
Jn,N

AC
Jn,N

J AC

100.2

102

100
100

100.1
101

102

103
N

(c) n = 100 basis functions

104

105

102

103

104
N

105

106

(d) Zoomed version of the average cost for different n

AC
Figure 2. The objective performance Jn,N
is computed using (18) for Example 7.1. The
AC
red dotted line denoted by J
is the optimal solution approximated by n = 103 and N = 106 .

AC
as the mean value across different experiments (solid and dashed lines) of the objective performance Jn,N
as
a function of the sample size N .

Figure 2(d) depicts a zoomed perspective of the means for the three cases of n. All the results in Figure 2 are
obtained based on 400 independent simulation experiments. It is perhaps not surprising that the optimal value
depicted in red dotted line is very close to the classical LQG example whose exact solution is analytically
available. It can be seen that the randomized approximations asymptotically converge, as suggested by
Theorem 6.1.
The simulation results suggest three interesting features concerning n, the number of basis functions: The
higher the number of basis functions,
(i) the smaller the approximation error (i.e., asymptotic distance for N → ∞ to the red dotted line),
(ii) the lower the variance of approximation with respect to the sampling distribution for each N , and

(iii) the slower the convergence behavior with respect to the sample size N .
The features (i) and (ii) are positive impacts of increasing the number of basis functions. While (i) is predicted
by Corollary 3.9, since the error due to the projection term becomes smaller, it is not entirely clear how to
formally explain (ii). On the contrary, the feature (iii) is indeed a negative impact, as a high number of basis
functions requires a large number of samples N to produce reasonable approximation errors. This phenomena
can be justified through the lens of Corollary 4.12 where the approximation errors grows proportionally to n.

26

3

102

UB − J LB
Jn,η
n,η

101
100
10−1
10−2
101

102
103
104
number of iterations k

UB
Jn,η

posterior approximation

prior & posterior error

ε

LB
Jn,η

J AC

2

1

105

UB − J LB
(a) A priori error ε and a posteriori error Jn,η
n,η

101

102
103
104
number of iterations k

105

UB and lower bound J LB
(b) Upper bound Jn,η
n,η

Figure 3. The results and error bounds are obtained by Algorithm 1 with n = 10 for
Example 7.1. The red dotted line is the optimal solution computed as indicated in Figure 2.
Structural convex optimization:
Algorithm 1 was implemented with the parameters described in Corollary 5.8 leading to deterministic
UB
LB
upper and lower bounds (Jn,η
and Jn,η
, respectively) for the cost function JnAC , see also Theorem 5.3. These
bounds are computationally appealing as they provide a posteriori bounds on the approximation error that
often is significantly smaller than the a priori bounds given by Theorem 5.3. This behavior can be seen in
the simulation results summarized in Figure 3 where the number of basis functions is n = 10. Similar to
Figure 2, the red dotted line is the optimal value of the original infinite program P, which we approximated
by using 103 basis functions and 106 iterations of Algorithm 1; it coincides with the one from the randomized
method.

7.2. Example 2: A fisheries management problem
A natural approximation approach toward dynamic programming problems goes through a discretization
scheme (e.g., discretization the state and/or action spaces). The main objective of this example is to compare
the proposed LP-based approximation of this article with more standard discretization schemes. To this end,
we borrow an example from [24, Section 1.3] and compare our results with the recent discretization method
proposed by [42]. Consider the population growth model, known as Ricker model,
st+1 = ϑ1 at exp(−ϑ2 at + ξt ),

t ∈ N,

where ϑ1 , ϑ2 ∈ R+ , st is the population size in season t, and at is the population to be left for spawning for
the next season, i.e., the difference st − at is the amount of fish captured in season t. The running reward
function, to be maximized is ψ(a, s) = ϕ(s − a), where ϕ is the so-called shifted isoelastic utility function
ϕ(z) := 3(z + 0.5)1/3 − (0.5)1/3 [16, 13, Section 4.1]. The state space is S = [κ, κ], for some κ, κ ∈ R+ .
Since the population left for spawning cannot be greater than the total population, for each s ∈ S, the set
of admissible actions is A(s) = [κ, s]. To fulfill Assumption 2.1(i), following the transformation suggested by
[42], we equivalently reformulate the above problem using the dynamics
st+1 = ϑ1 min(at , st ) exp(−ϑ2 min(at , st ) + ξt ),

t ∈ N,

where the admissible actions set is now the state-independent set A = [κ, κ], and the running reward function
is ψ(a, s) = ϕ(s − a)1{s≥a} . The noise process (ξt )t∈N is a sequence of i.i.d. random variables which have a

27

102

n=2
n = 10
n = 100

AC
Jn,N
+ 10

J AC

101
100

101

102

103

N
AC
Figure 4. The objective performance Jn,N
is computed using (18) for Example 7.2. The
red dotted line is the optimal value approximated by n = 103 and N = 106 , which amounts
to 0 as also reported in [42].

uniform density function g supported on the interval [0, λ]. Thus, the corresponding kernel is
Z 
1

g log ξ − log ϑ1 min(a, s) + ϑ2 min(a, s) dξ, ∀B ∈ B(R) .
Q(B|s, a) =
ξ
B

Note that to make the model consistent, we must have ϑ1 a exp(−ϑ2 a + ξ) ∈ [κ, κ] for all (a, ξ) ∈ [κ, κ] × [0, λ].
By defining an appropriate change of coordinate similar to Lemma 7.1, Assumption 2.1 are fulfilled; we refer
the reader to [42, Section 7.2] for further information and detailed analysis.
Simulation details: The chosen numerical values are λ = 0.5, ϑ1 = 1.1, ϑ2 = 0.1, κ = 7, and κ = 0.005.
Randomized approach:
We implement the methodology presented in Section 4.2, resulting in a finite random convex program
(18), where the uniform distribution on K = S × A = [κ, κ]2 is used to draw the random samples. Figure 4
illustrates three cases with the number of basis functions n ∈ {2, 10, 100} and the bound (21b). The colored
tubes represent the results between [10%, 90%] quantiles (shaded areas) as well as the means (solid lines)
AC
across 400 independent experiments of the objective performance Jn,N
as a function of the sample size N . It
is interesting to note that in this example the optimal solution is captured even with 2 basis functions and
only N = 20 random samples. This becomes even more attractive when we compare the results with a direct
discretization scheme depicted in [42, Figure 2].
Structural convex optimization:
Similar to the LQG example in Section 7.1, we also implement the smoothing methodology for the case of
n = 10. The simulation results are reported in Figure 5.
Appendix A. Infinite-Horizon Discounted-Cost Problems
In the Markov decision process setting, introduced in Section 2.1, let us consider long-run τ -discounted
cost (DC) problems with the discount factor τ ∈ (0, 1) and initial distribution ν ∈ P(X) described as
" n−1
#
X
DC
π
t
J (ν) := inf lim Eν
τ ψ(xt , at ) .
(38)
π∈Π n→∞

t=0

28

3
posterior approximation

prior & posterior error

ε
UB − J LB
Jn,η
n,η

101

10−1

10−3
101

102
103
104
number of iterations k

105

UB − J LB
(a) A priori error ε and a posteriori error Jn,η
n,η

UB
Jn,η
LB
Jn,η

J AC

2

1

0
101

102
103
104
number of iterations k

105

UB and lower bound J LB
(b) Upper bound Jn,η
n,η

Figure 5. The results and error bounds are obtained by Algorithm 1 with n = 10 for
Example 7.2. The red dotted line is the optimal solution computed as indicated in Figure 4.
As in the average cost setting, in Section 2, we assume that the control model satisfies Assumption 2.1. We
refer to [24, Chapter 4] and [26, Chapter 8] for a detailed exposition and required technical assumptions in
more general settings. As for the AC problems, it is well known that the DC problem (38) can be alternatively
characterized by means of infinite LPs (P) and (D) introduced in Section 3.1, where


(X, C) := (C(S), M(S))




(B, Y) := (C(K), M(K))





:= C+ (K)
K


 K∗ := M (K)
+
(39)
 c(B) = −ν(B), B ∈ B(S)




b(s, a) = −ψ(s, a)





A
: X → B, Ax(s, a) := −x(s) + τ Qx(s, a)


 A∗ : Y → C, A∗ y(B) := y(B × A) − τ yQ(B), B ∈ B(S),
Theorem A.1 (LP characterization [24, Theorem 6.3.8]). Under Assumption 2.1, the optimal value J DC of
the DC problem in (38) can be characterized by the LP problem (P) in the setting (39), in the sense that
J = −J DC .

It is known that under similar conditions as in Assumption 2.1 on the control model, the value function
u in the τ -discounted cost optimality equation is Lipschitz continuous [22, Section 2.6] or [17, Theorem 3.1].
We use the norms similar to the AC-setting (10). The next step toward studying the approximation error (8)
for the DC-setting readily follows by Theorem 3.3 combined with the following lemma.
⋆

Lemma A.2 (DC semi-infinite regularity). For the DC-problem (38), characterized by the dual-pair vector
spaces in (39), under Assumption 2.1 we have the operator norm kAk ≤ 1+max{LQ , 1}τ , the inf-sup constant
of Assumption 3.1(ii) γ = 1 − τ , and the dual optimizer norm
ky ⋆ kW ≤ θD =

θP + (1 − τ )−1 kψk∞
.
(1 − τ )θP − kψkL

(40)

Proof. With the norms considered and following a proof similar to Lemma 3.6, the operator norm kAk can
be upper bounded as kAk ≤ 1 + τ. The inf-sup condition, Assumption 3.1(ii), holds with γ = 1 − τ , since
inf sup

y∈K∗ x∈Xn

(1 − τ ) 1, y
Ax, y
≥ inf
= 1 − τ.
kxkkykW y∈K∗
kykW

29

Moreover kνkW = 1 since it is a probability measure. Thus, given the lower bound for the optimal value
JnDC ≥ −(1 − τ )−1 kψk∞ , the assertion of Proposition 3.2 (i.e, the dual optimizers bound in (6)) leads to the
desired assertion (40).

Note that when the norm constraint is neglected, the dual program enforces that any solution yn⋆ in the
program Dn satisfies x, A∗ yn⋆ − c = 0 for all x ∈ Xn (cf. the program D). Assume that a constant function
belongs to the set Xn . Then, the constraint evaluated at the constant function reduces to (1 − τ ) 1, yn⋆ =
(1 − τ )kyn⋆ kW = 1. It is worth noting that this observation can consistently be captured by Lemma A.2 when
θP tends to ∞, in which the bound (40) reduces to kyn⋆ kW ≤ (1 − τ )−1 .
References
[1] E. J. Anderson and P. Nash, Linear Programming in Infinite-dimensional Spaces, John Wiley & Sons, 1987.
[2] A. Arapostathis, V. Borkar, E. Fernndez-Gaucherand, M. Ghosh, and S. Marcus, Discrete-time controlled Markov
processes with average cost criterion: A survey, SIAM Journal on Control and Optimization, 31 (1993), pp. 282–344.
[3] A. Basu and V. S. Borkar, Stochastic control with imperfect models, SIAM Journal on Control and Optimization, 47
(2008), pp. 1274–1300.
[4] A. Ben-Tal, L. Ghaoui, and A. Nemirovski, Robust Optimization, Princeton University Press, 2009.
[5] D. Bertsekas, Convergence of discretization procedures in dynamic programming, IEEE Transactions on Automatic Control, 20 (1975), pp. 415–419.
[6] D. P. Bertsekas, Dynamic Programming and Optimal Control, Vol. II, Athena Scientific, 4th ed., 2012.
[7] D. P. Bertsekas and S. E. Shreve, Stochastic Optimal Control, vol. 139, Academic Press, Inc., 1978.
[8] D. P. Bertsekas and J. N. Tsitsiklis, Neuro-Dynamic Programming, Athena Scientific, 1st ed., 1996.
[9] D. Bertsimas, X. V. Doan, and J. Lasserre, Approximating integrals of multivariate exponentials: A moment approach,
Operations Research Letters, 36 (2008), pp. 205–210.
[10] M. C. Campi and S. Garatti, The exact feasibility of randomized solutions of uncertain convex programs, SIAM Journal
on Optimization, 19 (2008), pp. 1211–1230.
[11] O. L. V. Costa and F. Dufour, A linear programming formulation for constrained discounted continuous control for
piecewise deterministic Markov processes, Journal of Mathematical Analysis and Applications, 424 (2015), pp. 892–914.
[12] I. Csiszár, I-divergence geometry of probability distributions and minimization problems, The Annals of Probability, 3
(1975), pp. 146–158.
[13] D. P. de Farias and B. Van Roy, The linear programming approach to approximate dynamic programming, Operations
Research, 51 (2003), pp. 850–865.
[14] D. P. De Farias and B. Van Roy, On constraint sampling in the linear programming approach to approximate dynamic
programming, Mathematics of Operations Research, 29 (2004), pp. 462–478.
[15] O. Devolder, F. Glineur, and Y. Nesterov, First-order methods of smooth convex optimization with inexact oracle,
Mathematical Programming, (2013), pp. 1–39.
[16] F. Dufour and T. Prieto-Rumeau, Approximation of Markov decision processes with general state space, Journal of
Mathematical Analysis and Applications, 388 (2012), pp. 1254–1267.
[17] F. Dufour and T. Prieto-Rumeau, Finite linear programming approximations of constrained discounted Markov decision
processes, SIAM Journal on Control and Optimization, 51 (2013), pp. 1298–1324.
, Stochastic approximations of constrained discounted Markov decision processes, Journal of Mathematical Analysis
[18]
and Applications, 413 (2014), pp. 856–879.
[19]
, Approximation of average cost Markov decision processes using empirical distributions and concentration inequalities, Stochastics, 87 (2015), pp. 273–307.
[20] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004.
[21] R. T. Farouki, The Bernstein polynomial basis: A centennial retrospective, Computer Aided Geometric Design, 29 (2012),
pp. 379 – 419.
[22] O. Hernández-Lerma, Adaptive Markov Control Processes, vol. 79, Springer, 1989.
[23] O. Hernández-Lerma, J. González-Hernández, and R. López-Martı́nez, Constrained average cost Markov control
processes in Borel spaces, SIAM Journal on Control and Optimization, 42 (2003), pp. 442–468.
[24] O. Hernández-Lerma and J. Lasserre, Discrete-Time Markov Control Processes: Basic Optimality Criteria, Springer,
1996.
[25] O. Hernández-Lerma and J. Lasserre, Approximation schemes for infinite linear programs, SIAM Journal on Optimization, 8 (1998), pp. 973–988.
[26] O. Hernández-Lerma and J. Lasserre, Further Topics on Discrete-Time Markov Control Processes, Springer, 1999.
[27] O. Hernández-Lerma and J. B. Lasserre, The linear programming approach, in Further Topics on Discrete-Time Markov
Control Processes, Springer, 1999, pp. 203–249.

30

[28] N. Kariotoglou, S. Summers, T. Summers, M. Kamgarpour, and J. Lygeros, Approximate dynamic programming for
stochastic reachability, in European Control Conference (ECC), 2013, pp. 584–589.
[29] V. R. Konda and J. N. Tsitsiklis, On actor-critic algorithms, SIAM Journal on Control and Optimization, 42 (2003),
pp. 1143–1166.
[30] H. Kushner and G. G. Yin, Stochastic Approximation and Recursive Algorithms and Applications, vol. 35, Springer, 2003.
[31] H. Lai and S. Wu, Extremal points and optimal solutions for general capacity problems, Mathematical Programming, 54
(1992), pp. 87–113.
[32] J. B. Lasserre, Moments, Positive Polynomials and Their Applications, Imperial College Press, 2009.
[33] D. G. Luenberger, Optimization by Vector Space Methods, John Wiley & Sons Inc., 1969.
[34] P. Mohajerin Esfahani, D. Chatterjee, and J. Lygeros, Motion planning for continuous time stochastic processes: A
dynamic programming approach, IEEE Transactions on Automatic Control, 61 (2016), pp. 2155–2170.
, The stochastic reach-avoid problem and set characterization for diffusions, Automatica, 70 (2016), pp. 43–56.
[35]
[36] P. Mohajerin Esfahani, T. Sutter, and J. Lygeros, Performance bounds for the scenario approach and an extension
to a class of non-convex programs, IEEE Transactions on Automatic Control, 60 (2015), pp. 46–58.
[37] A. Nemirovski and A. Shapiro, Scenario approximations of chance constraints, in Probabilistic and Randomized Methods
for Design under Uncertainty, G. Calafiore and F. Dabbene, eds., Springer, 2006, pp. 3–47.
[38] Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Springer, 2004.
, Smooth minimization of non-smooth functions, Mathematical Programming, 103 (2005), pp. 127–152.
[39]
[40] S. Olver, On the convergence rate of a modified Fourier series, Mathematics of Computation, 78 (2009), pp. 1629–1645.
[41] M. Prandini and J. Hu, Stochastic reachability: Theory and numerical approximation, in Stochastic Hybrid Systems, CRC
Press, 2007, pp. 107–137.
[42] N. Saldi, S. Yüksel, and T. Linder, Asymptotic optimality of finite approximations to Markov decision processes with
general state and action spaces. available at arXiv:1503.02244, 2015.
[43] E. Shafieepoorfard, M. Raginsky, and S. P. Meyn, Rationally inattentive control of Markov processes, SIAM Journal
on Control and Optimization, 54 (2016), pp. 987–1016.
[44] M. Sion, On general minimax theorems, Pacific Journal of Mathematics, 8 (1958), pp. 171–176.
[45] T. Sutter, P. Mohajerin Esfahani, and J. Lygeros, Approximation of constrained average cost Markov control processes,
in 53rd IEEE Conference on Decision and Control (CDC), Dec 2014, pp. 6597–6602.
[46] J. N. Tsitsiklis and B. V. Roy, An analysis of temporal-difference learning with function approximation, IEEE Transactions on Automatic Control, 42 (1997), pp. 674–690.
[47] C. Villani, Topics in Optimal Transportation, American Mathematical Society, 2003.

