1

A Linearly Relaxed Approximate Linear
Program for Markov Decision Processes

arXiv:1704.02544v1 [] 9 Apr 2017

Chandrashekar Lakshminarayanan? , Shalabh Bhatnagar? , and Csaba Szepesvari†

Abstract
Approximate linear programming (ALP) and its variants have been widely applied to Markov
Decision Processes (MDPs) with a large number of states. A serious limitation of ALP is that it has an
intractable number of constraints, as a result of which constraint approximations are of interest. In this
paper, we define a linearly relaxed approximation linear program (LRALP) that has a tractable number
of constraints, obtained as positive linear combinations of the original constraints of the ALP. The main
contribution is a novel performance bound for LRALP.

Keywords:

Markov Decision Processes (MDPs), Approximate Linear Programming (ALP),
I. I NTRODUCTION

Markov decision processes (MDPs) have proved to be an indispensable model for sequential
decision making under uncertainty with applications in networking, traffic control, robotics, operations research, business, finance, artificial intelligence, health-care and more (see, e.g., [Whi93;
Rus96a; FS02; HY07; SB10; BR11; Put94; LL12; AA+15; BD17]). In this paper we adopt the
framework of discrete-time, discounted MDPs when a controller steers the stochastically evolving
state of a system while receiving rewards that depends on the states visited and actions chosen.
The goal is to choose the actions so as to maximize the return, defined as the total discounted
expected reward. A controller that uses past state information is called a policy. An optimal
policy is one that maximizes the value no matter where the process is started from [Put94].
In this paper we consider planning problems where the goal is to calculate actions of policies
?

Department of Computer Science and Automation, Indian Institute of Science, Bangalore 560012. E-mail: {chandrul,

shalabh}@csa.iisc.ernet.in
†

Department of Computing Science, University of Alberta, Edmonton, Alberta, Canada T6G 2E8. E-mail:

csaba.szepesvari@ualberta.ca

April 11, 2017

DRAFT

2

that give rise to high values and give new error bounds on the quality of solutions obtained by
solving linear programs of tractable size. To explain the contributions in more details, we start
by describing the computational challenges involved in planning.
The main objective of planning is to compute actions of an optimal policy while interacting
with an MDP model. In finite state-action MDPs, assuming access to individual transition
probabilities and rewards along transitions, various algorithms are available to perform this
computation in time and space that scales polynomially with the number of states and actions.
However, in most practical applications, the MDP is compactly represented and if it is not
infinite, the number of states scale exponentially with the size of the representation of the MDP.
If planners are allowed to perform some fixed amount of calculations for each state encountered,
it is possible to use sampling to make the per-state calculation-cost independent of the size of
the state space [Rus96b; Sze01; KMN02]. Nevertheless, the resulting methods are still quite
limited. In fact, various hardness results show that computing actions of (near-) optimal policies
is intractable in various senses and in various compactly represented MDPs [BT00]. Given these
negative results, it is customary to adopt the modest goal of efficiently computing actions of a
policy that is nearly as good as a policy chosen by a suitable (computationally unbounded, and
well-informed) oracle from a given restricted policy class. Here, within some restrictions (see
below), the policy class can be chosen by the user. The more flexibility the user is given in this
choice, the stronger a planning method is. The problem of planning with limited resources is also
one of the key problem in artificial intelligence (AI). The book of [KM12] gives a relatively fresh,
algorithm-centered summary of existing methods suitable for planning in MDPs. AI research tend
to focus on empirical results through the development of various benchmarks and little if any
effort is devoted to the theoretical understanding of the quality-effort tradeoff exhibited by the
that the various algorithms that are developed in this field.
A popular approach along these lines, which goes back to Schweitzer and Seidmann [SS85],
relies on considering linear approximations to the optimal value function: The idea is that,
similarly to linear regression, a fixed sequence of basis functions are combined linearly. The
user’s task is to use a priori knowledge of the MDP to choose the basis functions so that a good
approximation to the optimal value function will exist in the linear space spanned by the basis
functions. The idea then is to design some algorithm to find the coefficients of the basis functions
that gives a good approximation, while keeping computation cost in check. Finding a good
approximation is sufficient, since at the expense of an extra O(1/ε2 ) randomized computation,

April 11, 2017

DRAFT

3

a uniform O(ε)-approximation to the optimal value function can be used to calculate an action
of an O(ε)-optimal policy at any given state (e.g., follow the ideas in [Sze01; KMN02]; see
also Theorem 3.7 of Kallenberg [Kal17]). Since the number of coefficients can be much smaller
than the number of states, the algorithms that search for the coefficients have the potential to
run efficiently regardless of the number of states.
Following Schweitzer and Seidmann [SS85], most of the literature considers algorithms that
are obtained from restricting exact planning methods to search in the span of the fixed basis
functions when performing computations. In this paper we consider the so-called approximate
linear programming (ALP) approach, which was heavily studied during the last two decades,
e.g., [SP01; Gue+03; FV03; FV04; KH04; PZ09; DFM09; Tay+10; PP11; BFM12; AYBM14].
The basic idea here is to combine a linear program whose solution is the optimal value function
(and thus the number of optimization variables in it scales with the number of states) with a
linear constraint that restricts the optimization variables to lie in the subspace spanned by the
basis functions. As already noted by Schweitzer and Seidmann [SS85], the new LP can still
be kept feasible by just adding one special basis function, while by substituting the “value
function candidates” with their linear expansions, the number of optimization variables becomes
the number of basis functions. As shown by de Farias and Van Roy [FV03], the solution to the
resulting LP is within a constant factor of the best approximation to the optimal value function
within the span of the chosen bases. However, since the number of constraints in the LP is still
proportional to the number of states, it is not obvious whether a solution to the resulting LP can
be found in time independent of the number of states (other computations can be done in time
independent of the number of states, e.g., using sampling, at the price of a controlled increase
of the error, e.g., Theorem 6 of [PZ09]).
Most of the literature is thus devoted to designing methods to select a tractable subset of
the constraints while keeping the approximation guarantees, as well as keeping computations
tractable. Since a linear objective is optimized by a point on the boundary of the feasible region,
knowing the optimizer would be sufficient to eliminate all but as many constraints as the number
of optimization variables. The question is how to find a superset of these, or an approximating set,
without incurring much computational overhead. Schuurmans and Patrascu [SP01] and Guestrin
et al. [Gue+03] propose constraint generation in a setting where the MDP has additional structure
(i.e., factorized transition structure). This additional structure is then exploited in designing
constraint generation methods which are able to efficiently generate violated constraints. A more

April 11, 2017

DRAFT

4

general approach due to de Farias and Van Roy [FV04] is to choose a random subset of the
constraints by choosing states to be included at random from a distribution that reflects the
“importance” of states. While constraint generation can be powerful, it is not known how solution
quality degrades with the budget on the constraints generated (Guestrin et al. note that the number
of constraints generated can be at most exponential in a fundamental quantity, the induced
width of a so-called cost-network, which may be large and is in general hard to control). For
constraint sampling, de Farias and Van Roy [FV04] prove a bound on the suboptimality, but
this bound applies only in the unrealistic scenario when the constraints are sampled from an
idealized distribution, which is related to the stationary distribution of an optimal policy. While
it is possible to extend this result to any sampling distribution, the bound then scales with
the mismatch between the sampling and the idealized distributions, which, in general, will be
uncontrolled. Another weakness of the bound is related to that when constraints are dropped, the
linear program may become unbounded. To prevent this, de Farias and Van Roy [FV04] propose
imposing an extra constraint on the optimization variables. The bound they obtain, however,
scales with the worst approximation error over this constraint set. While in a specific example it
is shown that this error can be controlled, no general results are derived in this direction. Later
works, such as that of Desai, Farias, and Moallemi [DFM09] and Bhat, Farias, and Moallemi
[BFM12], repeat the analysis of de Farias and Van Roy [FV04] in combinations with other ideas.
However, no existing work that we know of addresses the above weaknesses of the result of
de Farias and Van Roy [FV04].
Another interesting approach is to consider the dual linear program, where the optimization
variables are measures over the state-action space and the feasible set is the set of discounted
state-occupation (DSO) measures of all possible policies. By adding an extra linear constraint on
the optimization variables, we arrive at an “approximate dual LP”. Feasibility of the resulting LP
can be ensured by adding basis functions that represent DSO measures of some select policies.
Abbasi-Yadkori, Bartlett, and Malek [AYBM14] considered this approach together and proposed
to use a randomized gradient method to minimize a penalized form of the linear objective to
approximately enforce the constraints. The algorithm computes the parameters of a measure over
the state-action space, from which a policy can be derived by normalization. The main result of
Abbasi-Yadkori, Bartlett, and Malek is a bound on the performance loss of this policy relative

April 11, 2017

DRAFT

5

to the performance of the best DSO measure in the feasible set of the approximate dual LP,1
while iteration cost to obtain an ε-competitive solution is O(1/ε4 ) provided that a number of
conditions hold. On the above complexity bound, the constants hidden are instance dependent,
but do not depend on the number of states or actions. The conditions under which the result is
proven are as follows: (i) the algorithm needs to be able to sample from distributions not too
dissimilar to the idealized distributions q1 , q2 , where q1 is a distribution over state-action pairs, q2
P
is a distribution over states and, e.g., q1 is defined by q1 (s, a) = kφ(s, a)k/ s0 ,a0 kφ(s0 , a0 )k with
φ(s, a) = (φ1 (s, a), . . . , φk (s, a))> and φ1 , . . . , φk being the chosen basis functions specifying the
linear constraints and (s, a) is a state-action pair; (ii) the Markov chains underlying all policies
in the MDP are uniformly fast mixing; (iii) for any state s0 ∈ S and index 1 ≤ i ≤ k, the
P
expression s,a φi (s, a)pa (s, s0 ) can be evaluated in O(1) time, where pa (s, s0 ) is the probability
of transitioning from state s to s0 provided action a is chosen. While the second assumption
limits the scope of MDPs that the result can be applied to, the other two assumptions limit the
choice of the basis functions. Among other things, it is unclear how feasibility can be ensured
while satisfying (i). Nevertheless, Abbasi-Yadkori, Bartlett, and Malek demonstrate promising
empirical results on a queuing problem.
Our main contribution is a new suboptimality bound for the case when the constraint system
is replaced with a smaller, linearly projected constraint system. We also propose a specific
way of adding the extra constraint to keep the resulting LP bounded. Rather than relying on
combinatorial arguments (such as those at the heart of de Farias and Van Roy [FV04]), our
argument uses previously unexploited geometric structure of the linear programs underlying
MDPs. As a result our bound avoids distribution-mismatch terms and we also remove the scaling
with worst approximation error. A specific outcome of our general result is the realization that
it is beneficial to select states so that the “feature vectors” of all states when scaled with a fixed
constant factor are included in the conic hull of the “feature vectors” underlying the selected
states. This suggests to choose the basis functions so that this property can be satisfied by
selecting only a few states. As we will argue, this property holds for several popular choices of
basis functions. A preliminary version of this paper without the theoretical analysis and without
the geometric arguments was published in a short conference communication [LB15].
1

The result shown is more general, allowing to use measures outside of the feasible set. However, for such measures the

performance bound degrades very rapidly and hence the greater generality does not seem to add much to the result.

April 11, 2017

DRAFT

6

II. BACKGROUND
The purpose of this section is to introduce the necessary background before we can present
the problem studied and the main results.
We shall consider finite state-action space, discounted total expected reward MDPs. We note
in passing that the assumption that number of states is finite is mainly made for convenience
and at the expense of a more technical presentation could be lifted. We will comment later
on the assumption concerning the number of actions. Let the set of states, or state space be
S = {1, 2, . . . , S} and let the set of actions be A = {1, 2, . . . , A}. For simplicity, we assume
that all actions are admissible in all states. Given a choice of an action a ∈ A in a state s ∈ S,
the controller incurs a reward (or gain) of ga (s) ∈ [0, 1] and the state moves to a next state
s0 ∈ S with probability pa (s, s0 ). A policy u is a mapping from states to actions.2 When a policy
is followed, the state sequence evolves as a Markov chain with transition probabilities given by
Pu matrix whose (s, s0 )th entry is Pu(s) (s, s0 ). Along the way the rewards generated from gu
.
defined by gu (s) = gu(s) (s). The value of following a policy from a starting state s is denoted
by Ju (s) and is defined as the expected total reward discounted reward. Thus,
∞
. X t t
α (Pu gu )(s) ,
Ju (s) =
t=0

where α ∈ (0, 1) is the so-called discount factor. We call Ju the value function of policy u. The
value function of a policy satisfies the fixed-point equation Ju = Tu Ju where the affine-linear
operator Tu is defined by Tu J = gu + αPu J. An optimal policy, is one that maximizes the
value simultaneously for all initial states. The optimal value function J ∗ is defined by J ∗ (s) =
maxu Ju (s) and is known to be the solution of the fixed-point equation J ∗ = T J ∗ where the
operator T is defined by (T J)(s) = maxu (Tu J)(s), s ∈ S, i.e., the maximization is componentwise. Optimal policies exist and in fact any policy u such that the equation Tu J ∗ = T J ∗ holds
is optimal (e.g., Corollary 3.3 of [Kal17]). A policy u is said to be greedy with respect to (w.r.t.)
J if Tu J = T J ∗ . Thus, any policy that is greedy w.r.t. J ∗ is optimal.
III. T HE L INEARLY R ELAXED ALP
In this section we introduce the computational model used and the “Linearly Relaxed Approximate Linear Program” a relaxation of the ALP.
2

For the scope of this paper, it suffices to restrict our attention to such policies as opposed to considering history dependent

policies. See Chapter 3, and specifically Corollary 3.3 of [Kal17].

April 11, 2017

DRAFT

7

As discussed in the introduction, we are interested in methods that compute a good approximation to the optimal value function. As noted earlier, at the expense of a modest additional
cost, knowing an O(ε) approximation to J ∗ at a few states suffices to compute actions of an
O(ε)-optimal policy. We will take a more general view, and we will consider calculating good
approximations to J ∗ with respect to a weighted 1-norm, where the weights c form a probability
distribution over S. Recall that the weighted 1-norm kJk1,c of a vector J ∈ RS is defined
P
S
as kJk1,c =
s c(s)|J(s)|. Note that here and in what follows we identify elements of R
(functions, mapping S = {1, . . . , S} to the reals) with elements of RS in the obvious way. This
P
allows us to write e.g. c> J, which denotes s c(s)J(s).
To introduce the optimization problem we study, first recall that the optimal value function
J ∗ is the solution of the fixed point equation T J ∗ = J ∗ . It follows from the definition of T
that J ∗ = maxu Tu J ∗ ≥ Tu J ∗ for any u, where ≥ is the componentwise partial ordering of
vectors (≤ is the reverse relation). With some abuse of notation, we also introduce Ta to denote
Tu where u(s) = a for any s ∈ S. It follows that J ∗ ≥ Ta J ∗ for any a ∈ A and also that
T = maxa Ta , where again the maximization is componentwise. We call a vector J that satisfies
J ≥ Ta J for any a ∈ A superharmonic. Note that this is a set of linear inequalities. By our
note on T and (Ta )a , these inequalities can also be written compactly as J ≥ T J. It is not
hard to show then that J ∗ is the smallest superharminoc function (i.e., for any J superharmonic,
.
J ≥ J ∗ ). It also follows that for any c ∈ RS++ = (0, ∞)S , the unique solution to the linear
program min{c> J : J ≥ T J} = min{c> J : J ≥ Ta J, a ∈ A} is J ∗ .
Now, let φ1 , . . . , φk : S → R be k basis functions. The Approximate Linear Program (ALP)
P
of Schweitzer and Seidmann [SS85] is obtained by adding the linear constraints J = ki=1 ri φi
to the above linear program. Eliminating J gives
X
X
X
min{
ri c> φi :
ri φi ≥ ga + α
ri Pa φi , a ∈ A, r = (ri ) ∈ Rk } .
i

i

i

As noted by Schweitzer and Seidmann [SS85], the linear program is feasible as long as 1,
defined as the vector with all components being identically equal to one, is in the span of
{φ1 , . . . , φk }. For the purpose of computations, it is assumed that the values c> φi , i = 1, . . . , k
and the values (Pa φi )(s) and ga (s) can be accessed in constant time. This assumption can be
relaxed to assuming that one can access ga (s) and φi (s) for any (s, a) in constant time, as well
as to that one can efficiently sample from c, from Pa (s, ·) for any (s, a) pair, but the details of
this are the beyond the scope of the present work. As shown by de Farias and Van Roy [FV03],

April 11, 2017

DRAFT

8

.
. P
if rALP denotes the solution to the above ALP then for JALP = i rALP (i)φi = ΦrALP it holds
that kJALP − J ∗ k1,c ≤

2ε
,
1−α

where ε = inf r kJ ∗ − Φrk∞ is the error of approximating the optimal

value with the span of the basis functions φ1 , . . . , φk and kJk∞ = maxs |J(s)| is the maximum
norm and Φ ∈ RS×k is the matrix formed by (φ1 , . . . , φk ). That the error of approximating J ∗
with JALP is O(ε) is significant: The user can focus on finding a good basis, leaving the search
for the “right” coefficients to a linear program solver.
While solving the ALP can be significantly cheaper than solving the LP underlying the MDP
and thus it can be advantageous for moderate-scale MDPs, the number of constraints in the
ALP is SA, hence the ALP is still intractable for huge-scale MDPs. To reduce the number of
constraints, we consider a relaxation of ALP where the constraints are replaced with positive
linear combinations of them. Recalling that the constraints took the form J ≥ ga + αPa J (with
J = Φr), choosing m to be target number of constraints, for 1 ≤ i ≤ m, the ith new constraint
P >
P >
is given by a wi,a
J ≥ a wi,a
(ga + αPa J), where the choice of m and that of the vectors
wi,a ∈ RS+ is left to the user. Note that this results in a linear program with k variables and m
constraints, which can be written as
min c> Φr

r∈Rk

s.t.

X
a

Wa> Φr ≥

X

Wa> (ga + αPa )Φr ,

(1)

a

where Wa = (w1,a , . . . , wm,a ) ∈ RS×m
. Note that the (i, j)th entry of the m × k constraint
+
P >
P >
matrix of the resulting LP is a wi,a φj −α a wi,a
Pa φj and assuming that (wi,a )a has p nonzero
elements, this can be calculated in O(p) time, making the total cost of obtaining the constraint
matrix to be O(mkp) regardless the value of S and A.
We will call the LP in (1) the linearly relaxed approximate linear program (LRALP). Any
LP obtained using any constraint selection/generation process can be represented by choosing
an appropriate binary-valued matrix W > = (W1> , . . . , WA> ) ∈ Rm×SA
. In particular, when the
+
constraints are selected in a random process as suggested by de Farias and Van Roy [FV04], the
matrix W would be a random, binary-valued matrix.
Note that the LRALP may be unbounded. Unboundedness could be avoided by adding an

April 11, 2017

DRAFT

9

extra constraint of the form r ∈ N to the LRALP, for a properly chosen polyhedron N ⊂ Rk .3
However, it seems to us that it is downright misleading to think that guaranteeing a bounded
solution will also lead to reasonable solutions. Thus we will stick to the above simple form,
forcing a discussion of how W should be chosen to get meaningful results.4
Further insight into the choice of W can be gained by considering the Lagrangians of the
ALP and LRALP. To write both LP’s in a similar form let us introduce E = (IS×S , . . . , IS×S )> ,
where IS×S is the S × S identity matrix. Further, let H : RS → RSA be the operator defined by
(HJ)> = ((T1 J)> , . . . , (TA J)> ) .
Note that H, which we call the linear Bellman operator, is a linear operator. Then, the ALP
can be written as
min{c> Φr | EΦr ≥ HΦr} ,

(ALP)

while LRALP takes the form
min{c> Φr | W > EΦr ≥ W > HΦr} .

(LRALP)

Hence, their Lagrangians are LALP (r, λ) = c> Φr + λ> (HΦr − EΦr) LLRALP (r, q) = c> Φr +
q > W > (HΦr − EΦr). Thus, we can view W q as a “linear approximation” to the dual variable
λ ∈ RSA
+ . This suggests that perhaps W should be chosen such that it approximates well the
optimal dual variable. If Φ spans RS , the optimal dual variable λ∗ is known to be the discounted
occupancy measure underlying the optimal policy (Theorem 3.18, [Kal17]), suggesting that the
role of W is very similar to the role of Φ excepts that the subspace spanned by the columns of
W should ideally be close to λ∗ .
3

In particular, to obtain their theoretical result, de Farias and Van Roy [FV04] need the assumption that the set N is bounded

and that it contains rALP . In fact, the error bound derived by de Farias and Van Roy depends on the worst error of approximating
J ∗ with Φr when r ranges over N . Hence, if N is unbounded, their bound is vacuous. In the context of a particular application,
de Farias and Van Roy [FV04] demonstrate that N can be chosen properly to control this term. However, no general construction
is presented to choose N .
4

The only question is whether there is some value in adding constraints beyond choosing W properly. Our position is that

the set N would most likely be chosen based on very little and general information; the useful knowledge is in choosing W ,
not in choosing some general set N . Since randomization does not guarantee bounded solutions, de Farias and Van Roy [FV03]
must use N : In their case, N incorporates all the knowledge that makes the LP bounded.

April 11, 2017

DRAFT

10

IV. M AIN R ESULTS
The purpose of this section is to present our main results. Let rLRA be a solution to the LRALP
given by (1) and let JLRA = ΦrLRA . When multiple solutions exist, we can choose any of them.
For the result, we assume that the LRALP is not unbounded, and hence a solution exist. In fact,
we will assume something much stronger. The discussion of why our assumptions are reasonable
and how to ensure that they hold is postponed to after the presentation of our results. Our main
results bounds the error kJ ∗ − JLRA k1,c .
The bound is given in terms of the approximation error of J ∗ with the basis functions Φ =
∗
∗
(φ1 , . . . , φk ), as well as the deviation between two functions, JALP
, JLRA
: S → R, which we

define next. In particular,
∗
JALP
(s) = min{r> φ(s) | Φr ≥ J ∗ , r ∈ Rk } ,
∗
JLRA
(s) = min{r> φ(s) | W > EΦr ≥ W > EJ ∗ , r ∈ Rk } ,

where s ∈ S. Recall that E : RS → RSA is defined so that (EJ)> = (J > , . . . , J > ), i.e., E
P
∗
stacks its argument A-fold. Hence, W > E = a Wa> . Our strong assumption is that JLRA
is
∗
finite-valued. Note that JALP
≥ J ∗ reflects the error due to using the basis functions (φj )j , and
∗
∗
the magnitude of the deviation JLRA
− JALP
reflects the error introduced due to the relaxed

constraint system.
Following de Farias and Van Roy [FV03; FV04], we will quantify the magnitude of the error
∗
∗
JLRA
− JALP
and also that of the error of approximating J ∗ with the subspace spanned by Φ, in

terms of a weighted maximum norm, kJk∞,ψ = maxs∈S |J(s)|/ψ(s), where ψ : S → R++ is a
positive-valued weighting function.5 As also stressed by de Farias and Van Roy, the appropriate
choice of ψ is crucial for MDPs with huge state-spaces: The problem is that if the range of
values of |J ∗ (s)| in different parts of the state space differ in orders of magnitude, we do not
expect to be able to control the error of approximating J ∗ uniformly over S. By choosing the
weighting function to reflect the magnitude of J ∗ , the weighted maximum norm is controlled as
soon as the relative errors are and this latter goal may be much easier to achieve than controlling
absolute errors.
5

As opposed to de Farias and Van Roy [FV03] and others, our definition uses division and not multiplication with the weights.

We choose this form for mathematical convenience: With this definition, nice duality results hold between weighted 1-norms
and weighted maximum norms.

April 11, 2017

DRAFT

11

Just like de Farias and Van Roy [FV03], we will also require that ψ is a stochastic Lyapunovfunction for the MDP. In particular, we require that the α-discounted stability coefficient
.
βψ = α max kPa ψk∞,ψ
a

(2)

is strictly less than one. This can be seen to imply that H : (RS , k · k∞,ψ ) → (RSA , k · k∞,ψ ) is
a contraction, where for J = (J1> , . . . , JA> )> ∈ RSA we let kJk∞,ψ = maxa kJa k∞,ψ . That H is
a contraction will play a crucial role in our results. Note that the condition βψ < 1 is closely
related to the condition that for any policy u, Pu ψ ≤ ψ, which can be viewed as a stability
condition on the MDP and which appeared in a slightly altered form in studying the stability
of MDPs with infinite state spaces [e.g., CM99]. Note also that one can always choose ψ = 1,
which gives β1 = α < 1. With this, we are ready to state our main result:
Theorem IV.1 (Error Bound for LRALP). Assume that c ∈ RS+ is such that 1> c = 1 and that
W ∈ RSA×m
is nonnegative valued. Let ψ ∈ RS+ be in the column span of Φ and assume that
+
the α-discounted stability coefficient of ψ is βψ < 1. Let ε = inf r∈Rk kJ ∗ − Φrk∞,ψ be the error
of approximation J ∗ using the basis functions in Φ. Then,
kJ ∗ − JLRA k1,c ≤

2c> ψ
∗
∗
(2.5ε + kJALP
− JLRA
k∞,ψ ) .
1 − βψ

Note that the result implicitly assumes that JLRA exists, because if JLRA does not exist then
∗
JLRA
is necessarily unbounded, making the last error term infinite. To ensure that ψ is in the

span of Φ, after choosing ψ, one can add ψ as one of the basis functions. Alternatively, the
bound can also be interpreted to hold for any ψ in the span of Φ with βψ < 1.
As noted earlier, de Farias and Van Roy [FV03] prove a similar error bound for JALP , the
solution of the ALP. In particular, their Theorem 3 states that under identical assumptions as in
our result, kJ ∗ − JALP k1,c ≤

2c> ψε
1−βψ

for ε defined as above (the result we cited previously is a

simplified form of this bound). The larger coefficient of ε is probably an artifact of our analysis.
Note that when W does not reduce the constraints, our bound is only a constant factor larger
∗
∗
than this previous result. The extra term kJALP
− JLRA
k∞,ψ can be seen as the price paid for

relaxing the constraints.
From linear programming theory, it follows that primal boundedness is equivalent to dual
feasibility. Since the dual of min{c> x : Ax ≥ b} is max{y > b : y ≥ 0, c = A> y}, we get
∗
that a necessary and sufficient condition for JLRA
to be finite-valued is that for any s ∈ S,
> >
φ(s) lies in the conic span, {U λ : λ ∈ RSA
+ }, of (the columns) of U = Φ E W . When W

April 11, 2017

DRAFT

12

is such that its constituents W1 , . . . , WA are all identical, the conic span of U is equal to the
conic span of Φ> W1 . It is particularly instructive to consider the case when the common matrix
Wa = (w1,a , . . . , wm,a ) “selects” the m states, i.e., when {w1,a , . . . , wm,a } = {es : s ∈ S0 } for
some S0 ⊂ S, |S0 | ≤ m, where es ∈ {0, 1}S are the s ∈ S vectors in the standard Euclidean
basis. In this case, the condition that φ(s) lies in the conic span of U is equivalent to that φ(s)
.
∗
lies in the conic span of φ(S0 ) = {φ(s0 ) : s0 ∈ S0 }. Thus, to ensure boundedness of JLRA
, the
chosen states should be selected to “conicly cover” all the vectors in φ(S) ⊂ Rk .6
The next theorem shows that magnitudes of the coefficients used in the conic cover control the
∗
∗
0
be the matrix of conic coefficients:
size of kJALP
− JLRA
k∞,ψ . For the theorem we let Λ ∈ RS×S
+
P
For any s ∈ S, φ(s) = s0 ∈S0 Λ(s, s0 )φ(s0 ). After the theorem we give constructions for creating

conic covers.
Theorem IV.2. Assume that W1 = · · · = WA , {w1,a , . . . , wm,a } = {es : s ∈ S0 } and that φ(S)
lies in the conic span of φ(S0 ) with conic coefficients given by Λ. Let ε = inf r kJ ∗ − Φrk∞,ψ .
Then,
∗
∗
∗
kJALP
− JLRA
k∞,ψ ≤ kJALP
− J ∗ k∞,ψ + (1 + kΛψkψ,∞ )ε .

Proof. Let r∗ be such that kJ ∗ − Φr∗ k∞,ψ = ε (this exists by continuity) and let δ = J ∗ − Φr∗ .
∗
Pick any s ∈ S and let rs = argmin{r> φ(s) : W > EΦr ≥ W > EJ ∗ , r ∈ Rk } so that JLRA
(s) =
∗
rs> φ(s). Note that by assumption, for any s0 ∈ S0 , JLRA
(s0 ) = rs> φ(s0 ) ≥ J ∗ (s0 ). Now, notice
∗
∗
∗
∗
that by definition, JLRA
≤ JALP
(the LP defining JLRA
is the relaxation of the LP defining JALP
).

Hence,
∗
∗
∗
∗
0 ≤ JALP
(s) − JLRA
(s) = JALP
(s) − J ∗ (s) + J ∗ (s) − JLRA
(s)
∗
and JLRA
(s) = rs> φ(s) = rs>

P

s0 ∈S0

Λ(s, s0 )φ(s0 ) =

P

s0 ∈S0

∗
Λ(s, s0 )JLRA
(s0 ) ≥

P

s0 ∈S0

Λ(s, s0 )J ∗ (s0 ).

Combining this with the previous inequality we get
∗
J ∗ (s) − J ∗ (s) J ∗ (s) −
J ∗ (s) − JLRA
(s)
≤ ALP
+
0 ≤ ALP
ψ(s)
ψ(s)

6

P

Λ(s, s0 )J ∗ (s0 )
.
ψ(s)

s0 ∈S0

The same implies that, under the same condition on W , boundedness of the LRALP holds if and only if

P

s

c(s)φ(s) is

in the conic span of φ(S0 ). Note that this is easy to fulfill if the support of c has a small cardinality by add all states in the
support of c to S0 .

April 11, 2017

DRAFT

13

Plugging in J ∗ (s) = φ(s)> r∗ +δ(s), using again that φ(s) =

P

s0 ∈S0

Λ(s, s0 )φ(s0 ), and also using

the triangle inequality after taking absolute values, we get
P
P
|J ∗ (s) − s0 ∈S0 Λ(s, s0 )J ∗ (s0 )|
Λ(s, s0 )|δ(s0 )|
|δ(s)|
0
≤
+ s ∈S0
ψ(s)
ψ(s)
ψ(s)
1 X
|δ(s0 )|
|δ(s)|
+
Λ(s, s0 )ψ(s0 )
≤
ψ(s)
ψ(s) s0 ∈S
ψ(s0 )
0
P
Λ(s, s0 )ψ(s0 )
|δ(s)|
0
≤
+ kδk∞,ψ s ∈S0
.
ψ(s)
ψ(s)
Combining this with the previous display and noting that kδk∞,ψ = ε finishes the proof.
Given φ : S → Rk , what is the minimum cardinality set S0 that conicly covers φ(S) and
how to find such a set? Further, how to keep the magnitude of kΛψk∞,ψ small? To control this
latter quantity it seems essential to make sure S0 contains states with high ψ-values. However, if
one is content with a bound that depends on kψk∞ , one can bound kΛψk∞,ψ by kψk∞ ζ where
P
ζ = maxs s0 ∈S0 Λ(s, s0 ), hence, the second term in the previous bound will be bounded by
(1 + kψk∞ ζ)ε.
Let us now return to the problem of finding conic covers. We will proceed by considering
some illustrative examples. As a start, consider the case when the basis functions are binary
valued. In this case, it is sufficient and necessary to choose one state for each binary vector
.
that appears in φ(S) ⊂ {0, 1}k . This gives that m0 = |S0 | ≤ 2k representative states will be
sufficient regardless of the cardinality of S. Further, in this case ζ = 1. For moderate to large
k (e.g., k  20), it will quickly become infeasible to keep 2k constraints. In this case we may
need to restrict what features are considered to guarantee the conic cover condition. Letting
Ai = {s ∈ S : φi (s) = 1}, if for a many pairs i 6= j, Ai and Aj do not overlap then N = |φ(S)|
can be much smaller than 2k . For example, in the commonly used state aggregation procedures
Ai ∩Aj = ∅ for any i 6= j, giving N = k. In the more interesting case of hierarchical aggregation
(when the sets {Ai } form a nested hierarchical partitioning of S), we have m0 ≤ D · k where
D is the depth of the hierarchy.
Another favourable example is the case of separable bases. In this case, the states are assumed
to be factored and the basis functions depend only on a few factors. Let us consider a simple
illustration. By abusing notation (redefining S), let S = S1 ×S2 , let there be k = 2 basis functions
and assume that φi (s) = hi (s) for some hi : Si → R, i = 1, 2. Assume further that 0 ∈ hi (Si )
for both i and specifically let si0 be such that hi (si0 ) = 0. In this case it is not hard to verify that

April 11, 2017

DRAFT

14

if Si0 is such that hi (Si ) is in the conic span of hi (Si0 ) then φ(S) is also in the conic span of
.
S0 = S1 × {s20 } ∪ {s10 } × S2 . The point is that |S0 | ≤ |S1 | + |S2 |, which is a tolerable increase of
growth. This example is not hard to generalize to more general, ANOVA-like basis expansions.
The moral is that as long as their limited order of interaction (which is usually necessary for
information theoretic reasons as well), the number of constraints may grow moderately with the
number of factors (dimensionality) of the state space.
In some cases, finding a conic cover with a small cardinality is not possible. This can already
happen in simple examples such as when S = {1, . . . , S} (as before) and φ(s) = (1, s, s2 )> .
In this case, the only choice is S0 = S. In examples similar to this one one possibility is to
quantize the range of φ, which may loose little on approximation quality, while it creates the
opportunity to construct a small cardinality conic cover.
Note that the bound of [FV04] and our main result can be seen as largely complementary.
Recall that de Farias and Van Roy consider adding an extra constraint r ∈ N , while they
propose to select all A constraints from the ALP corresponding to m states chosen at random
from some distribution µ. Then, with high probability, they show that, provided that rALP ∈ N ,
the extra price paid for relaxing the constraints of the ALP is O(ρεN k/m), where ρ = maxs

µ∗ (s)
,
µ(s)

µ∗ = (1 − α)c> (I − αPu∗ )−1 , u∗ is an optimal policy, and εN = supr∈N kJ ∗ − Φrk∞,ψ .7 The
bound is nontrivial when m ≥ ρεN k. In general, it may be hard to control ρ, or even εN while
ensuring that rALP ∈ N .
V. P ROOF OF T HEOREM IV.1
In this section we present the proof of the main result, Theorem IV.1. The proof uses
contraction-arguments. We will introduce a novel contraction operator, Γ̂ : RS → RS , that
captures the distortion introduced by the extra constraint in ALP and the relaxation in LRALP,
respectively. Then we relate the solution of LRALP to the fixed point of Γ̂.
∗
Note that for the proof it suffices to consider the case when JLRA
is finite-valued because

otherwise the bound is vacuous. Also, recall that it was assumed that ψ lies in the column space
of Φ, while βψ , the α-discounted stability of ψ w.r.t. the MDP (cf. (2)) is strictly below one.
We will let r0 ∈ Rk be such that ψ = Φr0 . We also assumed that the matrix W is nonnegative
P
valued, while c specifies a probability distribution over S: s c(s) = 1 and c ∈ RS+ .
7

The paper presents the results for µ = µ∗ giving ρ = 1, but the analysis easily extends to the general case.

April 11, 2017

DRAFT

15

The operator Γ̂ are defined as follows: For J ∈ RS , s ∈ S,
(Γ̂J)(s) = min{r> φ(s) : W > EΦr ≥ W > EHJ, r ∈ Rk } .
Note that (Γ̂J)(s) mimics the definition of ALP with c = es , except that the constraint J = Φr
is dropped.
Let us now recall some basic results from the theory of contraction maps. First, let us recall
the definition of contractions. Let k · k be a norm on RS and ρ > 0. We say that the map
B : RS → RS is (ρ, k · k)-Lipschitz if for any J, J 0 ∈ RS , kBJ − BJ 0 k ≤ ρkJ − J 0 k. We say
that B is a k · k-contraction with factor ρ if it is (ρ, k · k)-Lipschitz and ρ < 1. It is particularly
easy to check whether a map is a contraction map with respect to a weighted maximum norm
if it is known to be monotone. Here, B is said to be monotone if for any J ≤ J 0 , J, J 0 ∈ RS ,
BJ ≤ BJ 0 also holds, where ≤ is the componentwise partial order between vectors. We start
with the following characterization of monotone contractions with respect to weighted maximum
norms:
Lemma V.1. Let B : RS → RS , ψ : S → R++ , β ∈ (0, 1). The following are equivalent:
(i) B is a monotone contraction map with contraction factor β with respect to k · kψ,∞ .
(ii) For any J, J 0 ∈ RS , t ≥ 0, J ≤ J 0 + tψ implies that BJ ≤ BJ 0 + βtψ.
The proof, which essentially copies that of Lemma 3.1 of [Kal17], is given for completeness:
Proof. Introduce · to denote elementwise products: Thus, (ψ · J)(s) = ψ(s)J(s). We also let
ψ −1 (s) = 1/ψ(s) and we will use the shorthand k · k = k · k∞,ψ .
Let us first prove (i) ⇒ (ii). Thus, assume that B is a monotone contraction map with factor
β. Take any J, J 0 , t > 0, J ≤ J 0 + tψ. We have BJ = B(J + tψ) − BJ 0 + BJ 0 ≤ (ψ −1 · (B(J +
tψ) − BJ 0 )) · ψ + BJ 0 ≤ kB(J + tψ) − BJ 0 kψ + BJ 0 ≤ βtkJ − J 0 kψ + BJ 0 .
For the reverse direction, note that monotonicity follows by taking t = 0. Now, let ε = kJ−J 0 k.
Then, J ≤ J 0 + εψ and J 0 ≤ J + εψ. By monotonicity and the assumed property of B (using
t = ε ≥ 0), −βεψ ≤ BJ − BJ 0 ≤ βεψ, which implies that kBJ − BJ 0 k ≤ β.
Corollary V.2. If B is monotone and there exists some β ∈ [0, 1) such that for any J ∈ RS and
any t > 0,
B(J + tψ) ≤ BJ + βtψ

April 11, 2017

(3)

DRAFT

16

then B is a k · k∞,ψ contraction with factor β.
Proof. Let J, J 0 ∈ RS , t ≥ 0 and assume that J ≤ J 0 + tψ. By monotonicity BJ ≤ B(J 0 + tψ),
while by (3), B(J 0 + tψ) ≤ BJ 0 + βtψ. Hence, BJ ≤ BJ 0 + βtψ. This shows that (ii) of
Lemma V.1 holds. Hence, by this lemma, B is a contraction with factor β with respect to
k · k∞,ψ .
Let us now return to the proof of our main result. Recall that the goal is to bound kJ ∗ −
JLRA k1,c through relating this deviations from the fixed point of Γ̂, which was promised to be a
contraction. Let us thus now prove this. For this, it suffices to show that Γ̂ satisfies the conditions
of Corollary V.2. In fact, we will see this holds with β = βψ .
Proposition V.3. The operator Γ̂ satisfies the conditions of Corollary V.2 with β = βψ , and is
thus a k · k∞,ψ -contraction with coefficient βψ .
Proof. First, note that (as it is well known) H is monotone (all the Pa matrices in the definition
of H are nonnegative valued) and that it satisfies an inequality similar to (3): For any t ≥ 0,
J ∈ RS ,
H(J + tψ) ≤ HJ + βψ t Eψ .

(4)

This follows again because our assumption on ψ implies that for any a ∈ A, αPa ψ ≤ βψ ψ.
.
Let us now prove that Γ̂ is monotone. Given J ∈ RS , let F 0 (J) = { Φr : W > EΦr ≥
W > HJ, r ∈ Rk } . Choose any s ∈ S. Since J1 ≤ J2 , W is nonnegative valued and H is
monotone, we have W > HJ1 ≤ W > HJ2 . Hence, FJ2 ⊂ FJ1 and thus (Γ̂J1 )(s) ≤ (Γ̂J2 )(s).
Since s was arbitrary, monotonicity of Γ̂ follows.
Let us now turn to proving that (3) holds with β = βψ . By definition, for s ∈ S, t ≥ 0,
J ∈ RS , (Γ̂(J + tψ))(s) = min{r> φ(s) : W > EΦr ≥ W > H(J + tψ), r ∈ Rk }. By (4),
H(J + tψ) ≤ HJ + tβψ Eψ and hence W > H(J + tψ) ≤ W > (HJ + tβψ Eψ). Thus, (Γ̂(J +
tψ))(s) ≤ min{r> φ(s) : W > EΦr ≥ W > (HJ + tβψ Eψ), r ∈ Rk }.
To finish, we need the following elementary observation:
Claim V.4. Let A ∈ Ru×v , b ∈ Ru , d ∈ Rv and b0 = Ax0 for some x0 ∈ Rv . Then
min{d> x : Ax ≥ b + b0 , x ∈ Rv }
= min{d> y : Ay ≥ b, y ∈ Rv } + d> x0 .

April 11, 2017

DRAFT

17

Proof of Claim V.4. Set y = x − x0 .
Now, using Claim V.4 with A = W > EΦ, b = W > HJ, d = φ(s), b0 = tβψ W > Eψ and
x0 = tβψ r0 , thanks to Φr0 = ψ we have Ax0 = b0 . Hence the desired statement follows from
the claim.
Let us now return to bounding kJ ∗ − JLRA k1,c . For x ∈ R, let (x)− be the negative part
of x: (x)− = max(−x, 0). Then, |x| = x + 2(x)− . For a vector J ∈ RS , we will write (J)−
to denote the vector obtained by applying the negative part componentwise. We consider the
decomposition
kJLRA − J ∗ k1,c = c>(JLRA − J ∗ ) + 2c>(JLRA − J ∗ )− .

(5)

Let VLRA be the fixed point Γ̂. We know claim the following:
Claim V.5. We have JLRA ≥ VLRA , c> JALP ≥ c> JLRA .
Proof. The inequality c> JALP ≥ c> JLRA follows immediately from the definitions of JALP and
JLRA .
To prove the first part let s ∈ S, c = es and let rs be a solution to LRALP in (1). For s ∈ S,
let V0 (s) = mins0 ∈S rs>0 φ(s).
.
It suffices to show that V1 = Γ̂V0 ≤ V0 ≤ JLRA . Indeed, if this holds then Vn+1 = Γ̂Vn , n ≥ 1,
satisfies Vn+1 ≤ Vn and Vn → VLRA as n → ∞ since Γ̂ is a monotone contraction mapping.
Since rs>0 φ(s) ≥ rs> φ(s) also holds for any s, s0 ∈ S, we have V0 (s) = rs> φ(s). Also, since
JLRA (s) ≥ rs> φ(s), it follows that JLRA ≥ V0 . Now, fix some s ∈ S and define re0 s ,V0 be the
solution to the linear program defining (Γ̂V0 )(s). We need to show that V1 (s) = (Γ̂V0 )(s) =
(re0 s ,V0 )> φ(s) ≤ V0 (s). By the definition of re0 s ,V0 we know that (re0 s ,V0 )> φ(s) ≤ r> φ(s) holds
for any r ∈ Rk such that W > EΦr ≥ W > HV0 . Thus, it suffices to show that rs satisfies
W > EΦrs ≥ W > HV0 . By definition, rs satisfies W > EΦrs ≥ W > HΦrs . Hence, by the monotone
property of H and since W is nonnegative valued, it is sufficient if Φrs ≥ V0 . This however
follows from the definition of V0 .
Thanks to the previous claim, (JLRA − J ∗ )− ≤ (VLRA − J ∗ )− and c> JLRA ≤ c> JALP . Hence,
from (5) we get
kJLRA − J ∗ k1,c ≤ c>(JALP − J ∗ ) + 2c>(VLRA − J ∗ )− .

April 11, 2017

DRAFT

18

>

c ψ
ε, where recall that ε = inf r kJ ∗ −
By Theorem 3 of [FV03], the first term is bounded by 2 1−β
ψ

Φrk∞,ψ . Hence, it remains to bound the second term.
For this, note that for any J ∈ RS , (J)− ≤ |J| and also that kJk1,c ≤ c> ψkJk∞,ψ . Hence,
we switch to bounding kJ ∗ − VLRA k∞,ψ . A standard contraction argument gives
kJ ∗ − VLRA k∞,ψ = kJ ∗ − Γ̂J + Γ̂J ∗ − Γ̂VLRA k∞,ψ
≤ kJ ∗ − Γ̂J ∗ k∞,ψ + kΓ̂J ∗ − VLRA k∞,ψ
≤ kJ ∗ − Γ̂J ∗ k∞,ψ + βψ kΓ̂J ∗ − VLRA k∞,ψ .
Reordering and using another triangle inequality we get
kJ ∗ − VLRA k∞,ψ ≤

∗
∗
kJ ∗ − JALP
k∞,ψ + kJALP
− Γ̂J ∗ k∞,ψ
.
1 − βψ

∗
We bound the term kJ ∗ − JALP
k∞,ψ in the following lemma:
∗
Lemma V.6. We have kJ ∗ − JALP
k∞,ψ ≤ 2ε, where recall that ε = inf r∈Rk kJ ∗ − Φrk∞,ψ .

.
∗
∗
Proof. Let r∗ = argminr∈Rk kJ ∗ − Φrk∞,ψ . First, notice that JALP
≥ J ∗ . Hence, 0 ≤ JALP
− J ∗.
Now let r0 = r∗ + εr0 . Then, Φr0 = Φr∗ + εψ ≥ J ∗ , where the equality follows by the definition
of r0 and the inequality follows by the definition of ε. Hence, r0 is in the feasible set of the
∗
∗
∗
LP defining JALP
and thus JALP
≤ Φr0 . Thus, 0 ≤ JALP
− J ∗ ≤ Φr∗ − J ∗ + εψ. Dividing

componentwise by ψ, taking absolute value and then taking maximum of both sides gives the
result.
∗
The proof of the main result is finished by noting that Γ̂J ∗ = JLRA
and the chaining the

inequalities we derived.
VI. N UMERICAL I LLUSTRATION
In this section, we show via an example in the domain of controlled queues the consequences
of Theorem IV.2, which bounded the error when the constraints are chosen based on selecting
a set of representative states (further preliminary experimental results have been reported in
[LB15]).
Model: We ran the experiments in the context of a queuing model similar to the one in
Section 5.2 of [FV03]. We consider a (simple) small scale model so that we can compare
with the optimal policy. At the same time, we will use a small number of basis functions and
constraints, to “stress-test” the algorithm. The queuing system has a single queue with random

April 11, 2017

DRAFT

19

arrivals and departures. The state of the system is the queue length with the state space given by
S = {0, . . . , S − 1}, where S − 1 is the buffer size of the queue. The action set A = {1, . . . , A}
is related to the service rates. We let st denote the state at time t. The state at time t + 1
when action at ∈ A is chosen is given by st+1 = st + 1 with probability p, st+1 = st − 1 with
probability q(at ) and st+1 = st , with probability (1−p−q(at )). For states st = 0 and st = S −1,
the system dynamics is given by st+1 = st + 1 with probability p when st = 0 and st+1 = st − 1
with probability q(at ) when st = S − 1. The service rates satisfy 0 < q(1) ≤ . . . ≤ q(A) < 1
with q(A) > p so as to ensure ‘stabilizability’ of the queue. The reward associated with action
a ∈ A and state s ∈ S is given by ga (s) = −(s/N + q(a)3 ) (the idea here is to penalize higher
queue lengths and higher service rates).
Parameter Settings: We ran our experiments for S = 1000, A = 4 with q(1) = 0.2, q(2) =
0.4, q(3) = 0.6, q(4) = 0.8, p = 0.4 and α = 1 − S1 . The moderate size of S = 1000 enabled us
to compute the exact value of J ∗ (the most expensive part of the computation). We made use of
polynomial features in Φ (i.e., 1, s, . . . , sk−1 ) since they are known to work reasonably well for
this domain [FV03]. Note hat the conic span conditions will only be met with some lag, unless
all the constraints are selected. Hence, these features allow us to test the limits of the theory.
We chose k = 4, a low number, to counteract that the MDP is small scale.
Experimental Methodology: We compare two different sampling strategies (i) based on the
cone conditions, and (ii) based on constraint sampling. The two strategies are compared via
lookahead policies, wherein, the action at state s is obtained by computing the approximate
value functions of the next states and selecting the action that leads to the larger estimated
value. The details are as follows: Case (i): Except for the corner states i.e., s = 0 and s = 999,
each state 0 < s < S − 1 has two next states namely s0 = s − 1 and s0 = s + 1. We formulate two
separate LRALPs (or just one LRALP for s = 0 and s = S −1) for next states. When formulating
the LRALP for state s0 , we let c = es0 and choose the constraint corresponding to state s0 to
ensure the cone condition to be met for LRALP. We choose 5 more constraints corresponding to
states 1, 200, 400, 600, 800, 999 (uniformly spaced across the state space) and compute Jˆe 0 . The
s

number of constraints is kept very small as a way of emulating that in large-scale problems we
cannot expect a dense covering of the state-space when selecting the constraints. The lookahead
P
policy is formulated as uLRA (s) = argmina∈A ga (s) + s0 ∈S pa (s, s0 )Jˆes0 (s0 ). Case (ii): In a
manner similar to Case (i), we formulate two separate LRALPs for next states. However, as
opposed to the previous case, when formulating the LRALP for state s0 , we sample m = 6

April 11, 2017

DRAFT

20

states (defining the constraints) from a distribution dependent on s0 . We experimented with two
sampling distributions that lead to two the lookahead policies that we denote by uCS−ideal and
uCS , respectively. The sampling distribution that defines uCS−ideal is the sampling distribution
that minimizes the upper bound proved by de Farias and Van Roy [FV04]. In particular, the
−1
sampling distribution used at state s0 is cs0 = e>
s0 (1−α)(I −αPu∗ ) , with es denoting the standard

basis vector which is 1 in the sth co-ordinate and 0 in all the other co-ordinates. This sampling
distribution is used as a baseline; it is unrealistic to assume that one would be able to sample
from this distribution without access to the optimal policy u∗ , which is the quantity of ultimate
0

interest. As a more realistic approach, we also consider sampling from cs0 (s) = κ(1−α)(α)|s −s| ,
where κ > 0 is a normalization factor that ensures that cs0 is a distribution. Again, we sample
m = 6 states. This leads to the policy uCS .
0

4

−200
3

−400
−600

2

−800
1
0

200

400

600

J∗
JCS
JCS−ideal
JLRA

800

0

200

400

600

800 1,000

u∗
uCS

uCS−ideal
uLRA

Fig. 1. Results for a single-queue with polynomial features. On both figures the x axis represents the state space: the length
of the queue. The left-hand-side figure shows the value functions of the various policies computed, alongside with the optimal
value function (higher values are better), while the right-hand side subfigure shows the underlying policies. “CS” and “CS-ideal”
stand for constraint sampling, while LRA stands for choosing the constraints based on geometric principles proposed in the
paper. For further details, see the text.

The results are shown in Fig. 1. The right-hand-side figure shows the policies computed, while
the left-hand-side figure shows their value functions. Since constraint sampling (CS) produces
randomized results, we repeated the simulations 10 times. The results in all cases were quite

April 11, 2017

DRAFT

REFERENCES

21

close, hence we show the plot for a typical run. The plots show that the CS case with the ideal
sampler is slightly worse, which can be attributed to the fact that the in the case of ideal sampler,
the sampling distribution is concentrated near the start state s0 in comparison to the behaviour
0

of the distribution cs0 (s) = (1 − α)(α)|s −s| which distributes the mass more evenly. As can
be seen from the figure, choosing the constraints to (approximately) satisfy the constraint of
the theoretical results reliably produces better results: In fact, the value functions J ∗ and JLRA
are mostly on the top of each other. We expect that in larger domains, differences between
constraints chosen based on the principles discovered in this paper and choosing constraints in
more heuristic ways will lead to similar, or even larger differences. However, the study of this
is left for future work.
VII. C ONCLUSION
In this paper, we introduced and analyzed the linearly relaxed approximate linear program
(LRALP) whose constraints were obtained as positive linear combination of the original constraints of the ALP. The main novel contribution is a theoretical result which gives a geometrically
interpretable bound on the performance loss due to relaxing the constraint sets. Possibilities for
future work include extending the results to other forms of approximate linear programming
in MDPs (e.g., [DFM09]), exploring the idea of approximating dual variables and designing
algorithms that use the newly derived results to actively compute what constraints to select.
R EFERENCES
[AA+15]

M. Abu Alsheikh, D. T. Hoang, D. Niyato, H.-P. Tan, and S. Lin, “Markov decision
processes with applications in wireless sensor networks: A survey,” IEEE Comm.
Surveys & Tutorials, vol. 17, pp. 1239–1267, 2015.

[AYBM14]

Y. Abbasi-Yadkori, P. Bartlett, and A. Malek, “Linear programming for large-scale
Markov decision problems,” in ICML, 2014, pp. 496–504.

[BD17]

R. J. Boucherie and N. M. van Dijk, Eds., Markov decision processes in practice.
Springer, 2017, vol. 248.

[BFM12]

N. Bhat, V. Farias, and C. C. Moallemi, “Non-parametric approximate dynamic
programming via the kernel method,” in NIPS, 2012, pp. 386–394.

[BR11]

N. Bäuerle and U. Rieder, Markov decision processes with applications to finance.
Springer, 2011.

April 11, 2017

DRAFT

REFERENCES

[BT00]

22

V. D. Blondel and J. N. Tsitsiklis, “A survey of computational complexity results
in systems and control,” Automatica, vol. 36, pp. 1249–1274, 2000.

[CM99]

R.-R. Chen and S. P. Meyn, “Value iteration and optimization of multiclass queueing networks,” Queueing Systems, vol. 32, no. 1-3, pp. 65–97, 1999.

[DFM09]

V. V. Desai, V. F. Farias, and C. C. Moallemi, “A smoothed approximate linear
program,” in NIPS, 2009, pp. 459–467.

[FS02]

E. A. Feinberg and A. Shwartz, Handbook of Markov decision processes: Methods
and applications. Kluwer Academic Publishers, 2002.

[FV03]

D. P. de Farias and B. Van Roy, “The linear programming approach to approximate
dynamic programming,” Operations Research, vol. 51, pp. 850–865, 2003.

[FV04]

——, “On constraint sampling in the linear programming approach to approximate
dynamic programming,” Mathematics of Operations Research, vol. 29, pp. 462–
478, 2004.

[Gue+03]

C. Guestrin, D. Koller, R. Parr, and S. Venkataraman, “Efficient solution algorithms
for factored MDPs,” J. of Artificial Intelligence Research, vol. 19, pp. 399–468,
2003.

[HY07]

Q. Hu and W. Yue, Markov decision processes with their applications. Springer,
2007.

[Kal17]

L. Kallenberg. (2017). Markov decision processes: Lecture notes, [Online]. Available: https://goo.gl/yhvrph (visited on 03/16/2017).

[KH04]

B. Kveton and M. Hauskrecht, “Heuristic refinements of approximate linear programming for factored continuous-state Markov decision processes,” in ICAPS,
2004, pp. 306–314.

[KM12]

A. Kolobov and Mausam, “Planning with Markov decision processes: An AI
perspective,” Synthesis Lectures on Artificial Intelligence and Machine Learning,
2012.

[KMN02]

M. Kearns, Y. Mansour, and A. Y. Ng, “A sparse sampling algorithm for nearoptimal planning in large Markov decision processes,” Machine learning, vol. 49,
pp. 193–208, 2002.

[LB15]

C. Lakshminarayanan and S. Bhatnagar, “A generalized reduced linear program
for Markov Decision Processes,” in AAAI, 2015, pp. 2722–2728.

April 11, 2017

DRAFT

23

[LL12]

F. L. Lewis and D. Liu, Eds., Reinforcement learning and approximate dynamic
programming for feedback control. Wiley-IEEE Press, 2012.

[PP11]

J. Pazis and R. Parr, “Non-parametric approximate linear programming for MDPs,”
in AAAI, 2011.

[Put94]

M. L. Puterman, Markov decision processes: Discrete stochastic programming.
New York: John Wiley, 1994.

[PZ09]

M. Petrik and S. Zilberstein, “Constraint relaxation in approximate linear programs,” in ICML, 2009, pp. 809–816.

[Rus96a]

J. Rust, “Numerical dynamic programming in economics,” in Handbook of Computational Economics, vol. 1, Elsevier, North Holland, 1996, pp. 619–729.

[Rus96b]

——, “Using randomization to break the curse of dimensionality,” Econometrica,
vol. 65, pp. 487–516, 1996.

[SB10]

O. Sigaud and O. Buffet, Eds., Markov decision processes in artificial intelligence.
Wiley-ISTE, 2010.

[SP01]

D. Schuurmans and R. Patrascu, “Direct value-approximation for factored MDPs,”
in NIPS, 2001, pp. 1579–1586.

[SS85]

P. J. Schweitzer and A. Seidmann, “Generalized polynomial approximations in
Markovian decision processes,” Journal of Mathematical Analysis and Applications, vol. 110, pp. 568–582, 1985.

[Sze01]

Cs. Szepesvári, “Efficient approximate planning in continuous space Markovian
decision problems,” AI Communications, vol. 13, no. 3, pp. 163–176, 2001.

[Tay+10]

G. Taylor, M. Petrik, R. Parr, and S. Zilberstein, “Feature selection using regularization in approximate linear programs for Markov decision processes,” in ICML,
2010, pp. 871–878.

[Whi93]

D. J. White, “A survey of applications of markov decision processes,” The Journal
of the Operational Research Society, vol. 44, no. 11, pp. 1073–1096, 1993.

April 11, 2017

DRAFT

