An MCMC Approach to Solving Hybrid Factored MDPs

1 Introduction

Markov decision processes (MDPs) [Bellman, 1957; Puterman, 1994] offer an elegant mathematical framework for solving sequential decision problems in the presence of uncertainty. However, traditional techniques for solving MDPs are computationally infeasible in real-world domains, which are factored and contain both discrete and continuous state and action variables. Recently, approximate linear programming (ALP) [Schweitzer and Seidmann, 1985] has emerged as a promising approach to solving large factored MDPs [de Farias and Roy, 2003; Guestrin et al., 2003]. This work focuses on hybrid ALP (HALP) [Guestrin et al., 2004] and addresses its major computational bottleneck - constraint satisfaction in the domains of discrete and continuous variables.

If the state and action variables are discrete, HALP involves an exponential number of constraints, and if any of the variables are continuous, the number of constraints is infinite. To approximate the constraint space in HALP, two techniques have been proposed: ε-HALP [Guestrin et al., 2004] and Monte Carlo constraint sampling [de Farias and Roy, 2004; Hauskrecht and Kveton, 2004]. The ε-HALP formulation relaxes the continuous portion of the constraint space to an εgrid, which can be compactly satisfied by the methods for discrete-state ALP [Guestrin et al., 2001; Schuurmans and

Milos Hauskrecht

Department of Computer Science University of Pittsburgh milos@cs.pitt.edu

Patrascu, 2002]. However, these methods are exponential in the treewidth of the discretized constraint space, which limits their application to real-world problems. In addition, the εgrid discretization is done blindly and impacts the quality of the approximation. Monte Carlo methods [de Farias and Roy, 2004; Hauskrecht and Kveton, 2004] offer an alternative to the ε-grid discretization and approximate the constraint space in HALP by its finite sample. Unfortunately, the efficiency of Monte Carlo methods is heavily dependent on an appropriate choice of sampling distributions. The ones that yield good approximations and polynomial sample size bounds are closely related to the optimal solutions and rarely known a priori [de Farias and Roy, 2004].

To overcome the limitations of the discussed constraint satisfaction techniques, we propose a novel Markov chain Monte Carlo (MCMC) method for finding the most violated constraint of a relaxed HALP. The method directly operates in the domains of continuous variables, takes into account the structure of factored MDPs, and its space complexity is proportional to the number of variables. Such a separation oracle can be easily embedded into the ellipsoid or cutting plane method for solving linear programs, and therefore constitutes a key step towards solving HALP efficiently.

The paper is structured as follows. First, we introduce hybrid MDPs and HALP [Guestrin et al., 2004], which are our frameworks for modeling and solving large-scale stochastic decision problems. Second, we review existing approaches to solving HALP and discuss their limitations. Third, we compactly represent the constraint space in HALP and formulate an optimization problem for finding the most violated constraint of a relaxed HALP. Fourth, we design a Markov chain to solve this optimization problem and embed it into the cutting plane method. Finally, we test our HALP solver on a set of large control problems and compare its performance to alternative approaches.

2 Hybrid factoredMDPs

Factored MDPs [Boutilier et al., 1995] allow a compact representation of large stochastic planning problems by exploiting their structure. In this section, we review hybrid factored MDPs [Guestrin et al., 2004], which extend this formalism to the domains of discrete and continuous variables.

A hybrid factored MDP with distributed actions (HMDP) [Guestrin et al., 2004] is a 4-tuple M = (X, A, P, R), where X = {X1, . . . ,Xn} is a state space represented by a set of state variables, A = {A1, . . . , Am} is an action space represented by a set of action variables, P(X′ | X,A) is a stochastic transition model of state dynamics conditioned on the preceding state and action choice, and R is a reward model assigning immediate payoffs to state-action configurations 1 .

State variables: State variables are either discrete or continuous. Every discrete variable Xi takes on values from a finite domain Dom(Xi). Following Hauskrecht and Kveton 2004, we assume that every continuous variable is bounded to the [0, 1] subspace. The state is represented by a vector of value assignments x = (xD, xC) which partitions along its discrete and continuous components xD and xC.

Action variables: The action space is distributed and represented by action variables A. The composite action is defined by a vector of individual action choices a = (aD, aC) which partitions along its discrete and continuous components aD and aC .

Transition model: The transition model is given by the conditional probability distribution P(X′ | X,A), where X and X′ denote the state variables at two successive time steps. We assume that the model factors along X′ as P(X′ | X,A) = ∏ni=1 P (X′ i | Par(X′ i )) and can be compactly represented by a dynamic Bayesian network (DBN) [Dean and Kanazawa, 1989]. Usually, the parent set Par(X′ i ) ⊆ X∪A is a small subset of state and action variables which allows for a local parameterization of the model.

Parameterization of transition model: One-step dynamics of every state variable is described by its conditional probability distribution P(X′ i | Par(X′ i )). If X′ i is a continuous variable, its transition function is represented by a mixture of beta distributions [Hauskrecht and Kveton, 2004]:

where πij is the weight assigned to the j-th component of the mixture, and αj = φαij (Par(X′ i )) and βj = φβij (Par(X′ i )) are arbitrary positive functions of the parent set. The mixture of beta distributions provides a very general class of transition functions and yet allows closed-form solutions to the integrals in HALP. If X′ i is a discrete variable, its transition model is parameterized by |Dom(X′ i )| nonnegative discriminant functions θj = φθij (Par(X′ i )) [Guestrin et al., 2004]:

Optimal value function and policy: The quality of a policy is measured by the infinite horizon discounted reward E[∑∞t=0 γ t rt], where γ ∈ [0, 1) is a discount factor and rt is the reward obtained at the time step t. This optimality criterion guarantees that there always exists an optimal policy π∗

Reward model: The reward function is an additive function R(x, a) = ∑j Rj(xj, aj) of local reward functions defined on the subsets of state and action variables Xj and Aj .

1 General state and action space MDP is an alternative name for a hybrid MDP. The term hybrid does not refer to the dynamics of the model, which is discrete-time.

which is stationary and deterministic [Puterman, 1994]. The policy is greedy with respect to the optimal value function V ∗ , which is a fixed point of the Bellman equation [Bellman, 1957; Bertsekas and Tsitsiklis, 1996]:

3 Hybrid ALP

Value iteration, policy iteration, and linear programming are the most fundamental dynamic programming (DP) methods for solving MDPs [Puterman, 1994; Bertsekas and Tsitsiklis, 1996]. However, their computational complexity grows exponentially in the number of used variables, which makes them unsuitable for factored MDPs. Moreover, they are built on the assumption of finite support for the optimal value function and policy, which may not exist if continuous variables are present. Recently, Feng et al. 2004 showed how to solve general state-space MDPs by performing DP backups of piecewise constant and piecewise linear value functions. This approximate method has a lower scale-up potential than HALP, but does not require the design of basis functions.

Linear value function: Value function approximation is a standard approach to solving large factored MDPs. Due to its favorable computational properties, linear value function approximation [Bellman et al., 1963; Roy, 1998]:

has become extremely popular in recent research [Guestrin et al., 2001; Schuurmans and Patrascu, 2002; de Farias and Roy, 2003; Hauskrecht and Kveton, 2004; Guestrin et al., 2004]. This approximation restricts the form of the value function V w to the linear combination of |w| basis functions fi(x), where w is a vector of tunable weights. Every basis function can be defined over the complete state space X, but often is restricted to a subset of state variables Xi.

3.1 HALP formulation

Various methods for fitting of the linear value function approximation have been proposed and analyzed [Bertsekas and Tsitsiklis, 1996]. We adopt a variation on approximate linear programming (ALP) [Schweitzer and Seidmann, 1985], hybrid ALP (HALP) [Guestrin et al., 2004], which extends this framework to the domains of discrete and continuous variables. The HALP formulation is given by:

where αi denotes basis function relevance weight:

ψ(x) is a state relevance density function weighting the quality of the approximation, and Fi(x, a) = fi(x) - γgi(x, a)

is the difference between the basis function fi(x) and its discounted backprojection:

We say that the HALP is relaxed if only a subset of the constraints is satisfied.

The HALP formulation reduces to the discrete-state ALP [Schweitzer and Seidmann, 1985; Schuurmans and Patrascu, 2002; de Farias and Roy, 2003; Guestrin et al., 2003] if the state and action variables are discrete, and to the continuousstate ALP [Hauskrecht and Kveton, 2004] if the state variables are continuous. The quality of this approximation was studied by Guestrin et al. 2004 and bounded with respect to minw ‖V ∗ - V w‖∞,1/L, where ‖·‖∞,1/L is a maxnorm weighted by the reciprocal of the Lyapunov function L(x) = ∑i wL i fi(x). The integrals in the objective function (Equation 4) and constraints (Equation 5) have closed-form solutions if the basis functions and state relevance densities are chosen appropriately [Hauskrecht and Kveton, 2004]. For example, the mixture of beta transition model (Equation 1) yields a closed-form solution to Equation 5 if the basis function fi(x′ ) is a polynomial. Finally, we need an efficient constraint satisfaction procedure to solve HALP.

3.2 Constraint satisfaction in HALP

If the state and action variables are discrete, HALP involves an exponential number of constraints, and if any of the variables are continuous, the number of constraints is infinite. To approximate such a constraint space, two techniques have been proposed recently: ε-HALP [Guestrin et al., 2004] and Monte Carlo constraint sampling [de Farias and Roy, 2004; Hauskrecht and Kveton, 2004].

ε-HALP: The ε-HALP formulation relaxes the continuous portion of the constraint space to an ε-grid by the discretization of continuous variables XC and AC. The new constraint space spans discrete variables only and can be compactly satisfied by the methods for discrete-state ALP [Guestrin et al., 2001; Schuurmans and Patrascu, 2002]. For example, Schuurmans and Patrascu 2002 search for the most violated constraint with respect to the solution w(t) of a relaxed ALP:

and add it to the linear program. If no violated constraint is found, w(t) is an optimal solution to the ALP.

The space complexity of both constraint satisfaction methods [Guestrin et al., 2001; Schuurmans and Patrascu, 2002] is exponential in the treewidth of the constraint space. This is a serious limitation because the cardinality of discretized variables grows with the resolution of the ε-grid. Roughly, if the discretized variables are replaced by binary, the treewidth increases by a multiplicative factor of log2(1/ε + 1), where (1/ε + 1) is the number of discretization points in a single dimension. Therefore, even problems with a relatively small treewidth are intractable for small values of ε. In addition, the ε-grid discretization is done blindly and impacts the quality of the approximation.

Monte Carlo constraint sampling: Monte Carlo methods approximate the constraint space in HALP by its finite sample. De Farias and Van Roy 2004 analyzed constraint sampling for discrete-state ALP and bounded the sample size for achieving good approximations by a polynomial in the number of basis functions and state variables. Hauskrecht and Kveton 2004 applied random constraint sampling to solve continuous-state factored MDPs and later refined their sampler by heuristics [Kveton and Hauskrecht, 2004].

Monte Carlo constraint sampling is easily applied in continuous domains and can be implemented in a space proportional to the number of variables. However, proposing an efficient sampling procedure that guarantees a polynomial bound on the sample size is as hard as knowing the optimal policy itself [de Farias and Roy, 2004]. To lower a high number of sampled constraints in a relaxed HALP, Monte Carlo samplers can be embedded into the cutting plane method. A technique similar to Kveton and Hauskrecht 2004 yields significant speedup with no drop in the quality of the approximation.

4 MCMC constraint sampling

To address the deficiencies of the discussed constraint satisfaction techniques, we propose a novel Markov chain Monte Carlo (MCMC) method for finding the most violated constraint of a relaxed HALP. Before we proceed, we compactly represent the constraint space in HALP and formulate an optimization problem for finding the most violated constraint in this representation.

4.1 Compact representation of constraints

Guestrin et al. 2001 and Schuurmans and Patrascu 2002 showed that the compact representation of constraints is essential in solving ALP efficiently. Following their ideas, we define violation magnitude τ w(x, a):

to be the amount by which the solution w violates the constraints of a relaxed HALP. We represent τ w(x, a) compactly by an influence diagram (ID), where X and A are decision nodes and X′ are random variables. The ID representation is built on the transition model P(X′ | X,A), which is already factored and contains dependencies among the variables X, X′ , and A. We extend the diagram by three types of reward nodes, one for each term in Equation 7: Rj = Rj(xj , aj) for every local reward function, Hi = -wifi(x) for every basis function, and Gi = γwifi(x′ ) for every backprojection. The construction is completed by adding arcs that represent the dependencies of the reward nodes on the variables. Finally, we verify that τ w(x, a) = EP (x′ |x,a)[∑i (Hi+Gi)+∑j Rj ]. Therefore, the decision that maximizes the expected utility in the ID corresponds to the most violated constraint.

We conclude that any algorithm for solving IDs can be used to find the most violated constraint. Moreover, special properties of the ID representation allow its further simplification. If the basis functions are chosen conjugate to the transition model (Section 3.1), we obtain a closed-form solution to EP(x′ |x,a)[Gi] (Equation 5), and thus the random variables X′ can be marginalized out of the diagram. This new representation contains no random variables and is know as a cost network [Guestrin et al., 2001].

4.2 Separation oracle

To find the most violated constraint in the cost network, we use the Metropolis-Hastings (MH) algorithm [Metropolis et al., 1953; Hastings, 1970] and construct a Markov chain whose invariant distribution converges to the vicinity of argmaxz τ w(z), where z = (x, a) and Z = X ∪ A is a joint set of state and action variables. The MetropolisHastings algorithm accepts the transition from a state z to a proposed state z∗ with the acceptance probability A(z, z∗ ) = min {1, p(z∗ )q(z|z∗ ) p(z)q(z∗ |z) }, where q(z∗ | z) is a proposal distribution and p(z) is a target density. Under mild restrictions on p(z) and q(z∗ | z), the chain always converges to the target density p(z) [Andrieu et al., 2003]. In the rest of this section, we discuss the choice of p(z) and q(z∗ | z) to solve our optimization problem.

Target density: The violation magnitude τ w(z) is turned into a density by the transformation p(z) = exp[τ w(z)]. Due to its monotonic character, p(z) retains the same set of global maxima as τ w(z), and thus the search for argmaxz τ w(z) can be performed on p(z). To prove that p(z) is a density, we show that it has a normalizing constant ∑zD ∫zC p(z) dzC, where zD and zC are the discrete and continuous components of the vector z = (zD, zC). As the integrand zC is restricted to the finite space [0, 1]|ZC| , the integral is proper as long as p(z) is bounded, and therefore it is Riemann integrable and finite. To prove that p(z) is bounded, we bound τ w(z). Let Rmax denote the maximum one-step reward in the HMDP. If the basis functions are of unit magnitude, wi can be typically bounded by |wi| ≤ γ(1-γ)-1Rmax, and consequently |τ w(z)| ≤ (|w| γ(1 - γ)-1 + 1)Rmax. Therefore, p(z) is bounded and can be treated as a density function.

To find the mode of p(z), we adopt the simulating annealing approach [Kirkpatrick et al., 1983] and simulate a non-homogeneous Markov chain whose invariant distribution equals to p1/Tt (z), where Tt is a decreasing cooling schedule with limt→∞ Tt = 0. Under weak regularity assumptions on p(z), p∞(z) is a probability density that concentrates on the set of global maxima of p(z) [Andrieu et al., 2003]. If the cooling schedule decreases such that Tt ≥ c/ log2(t + 2), where c is a problem-specific constant independent of t, the chain converges to the vicinity of argmaxz τ w(z) with the probability converging to 1 [Geman and Geman, 1984]. However, this schedule can be too slow in practice, especially for a high initial temperature c. Following the suggestion of Geman and Geman 1984, we overcome this limitation by selecting a smaller value of c than is required by the convergence criterion. As a result, convergence to the global optimum argmaxz τ w(z) is no longer guaranteed.

Proposal distribution: We take advantage of the factored character of Z and adopt the following proposal distribution [Geman and Geman, 1984]:

where z-i and z∗-i are the assignments to all variables but Zi in the original and proposed states. If Zi is a discrete variable, its conditional p(z∗ i | z-i) = p(z1,...,zi-1,z∗ i ,zi+1,...,zn+m) ∑ zi p(z1,...,zi-1,zi,zi+1,...,zn+m) can be derived in a closed form. If Zi is a continuous variable, a closed form of its cumulative density function is not likely to exist. To allow sampling from its conditional, we embed another MH step within the original chain. In the experimental section, we use the Metropolis algorithm with the acceptance probability A(zi, z∗ i ) = min {1, p(z∗ i |z-i) p(zi|z-i) }, where zi and z∗ i correspond to the original and proposed values of Zi. Note that sampling from both conditionals can be performed in the space of τ w(z) and locally.

Finally, we get a non-homogenous Markov chain with the acceptance probability A(z, z∗ ) = min {1, p1/Tt-1 (z∗ i |z-i) p1/Tt-1 (zi |z-i) } that converges to the vicinity of the most violated constraint. A similar chain was derived by Yuan et al. 2004 and applied to find the maximum a posteriori (MAP) configuration of random variables in Bayesian networks.

4.3 Constraint satisfaction

If the MCMC oracle converges to a violated constraint (not necessarily the most violated) in a polynomial time, it guarantees that the ellipsoid method can solve HALP in a polynomial time [Bertsimas and Tsitsiklis, 1997]. However, convergence of our chain within an arbitrary precision requires an exponential number of steps [Geman and Geman, 1984]. Even if this bound is too weak to be of practical interest, it suggests that the time complexity of finding a violated constraint dominates the time complexity of solving HALP. Therefore, the search for violated constraints should be performed efficiently. Convergence speedups that directly apply to our work include hybrid Monte Carlo (HMC) [Duane et al., 1987], slice sampling [Higdon, 1998], and RaoBlackwellization [Casella and Robert, 1996].

To evaluate the MCMC oracle, we embed it into the cutting plane method for solving linear programs, which results in a novel approach to solving HALP. As the oracle is not guaranteed to converge to the most violated constraint, we run the cutting plane method for a fixed number of iterations rather than having the same stopping criterion as Schuurmans and Patrascu 2002 (Section 3.2). Our MCMC solver differs from the Monte Carlo solver in that it samples constraints based on their potential to improve the existing solution, which substitutes for an unknown problem-specific sampling distribution. Comparing to the ε-HALP method, the MCMC oracle directly operates in the domains of continuous variables and its space complexity is linear in the number of variables.

5 Experiments

The performance of the MCMC solver is evaluated against two alternative constraint satisfaction techniques: ε-HALP and uniformMonte Carlo sampling. Due to space limitations, our comparison focuses on two irrigation-network problems [Guestrin et al., 2004], but our conclusions are likely to generalize across a variety of control optimization tasks. The irrigation-network problems are challenging for state-of-theart MDP solvers due to the factored state and action spaces

(Figure 1). The goal of an irrigation network operator is to select discrete water-routing actions AD to optimize continuous water levels XC in multiple interconnected irrigation channels. The transition model is parameterized by beta distributions and represents water flows conditioned on the operation modes of regulation devices. The reward function is additive and given by a mixture of two normal distributions for each channel. The optimal value function is approximated by a linear combination of four univariate piecewise linear basis functions for each channel [Guestrin et al., 2004]. The state relevance density function ψ(x) is uniform. A comprehensive description of the irrigation-network problems can be found in Guestrin et al. 2004.

The ε-HALP solver is implemented with the method of Schuurmans and Patrascu 2002 as described in Section 3.2. The Monte Carlo solver uniformly generates one million constraints and establishes a baseline for the comparison to an uninformatively behaving sampling method. The chain of the MCMC oracle is simulated for 500 steps from the initial temperature c = 0.2, which yields a decreasing cooling schedule from T0 = 0.2 to T500 ≈ 0.02. These parameters were chosen empirically to demonstrate the characteristics of our approach rather than to maximize the performance of the MCMC oracle. All experiments were performed on a Dell Precision 340 workstation with 2GHz Pentium 4 CPU and 1GB RAM. All linear programs were solved by the simplex method in the LP SOLVE package. The results of the experiments are reported in Figure 1.

Based on our results, we draw the following conclusions. First, the MCMC solver (N = 250) achieves the highest objective values on all problems. Higher objective values can be interpreted as closer approximations to the constraint space in HALP since the solvers operate on relaxed versions of HALP. Second, the quality of the MCMC policies (N = 250) surpasses the Monte Carlo ones while both solvers consume approximately the same computation time. This result is due to the informative search for violated constraints in the MCMC solver. Third, the quality of the MCMC policies (N = 250) is close to the ε-HALP ones (ε = 1/16), but does not surpass them significantly. In the irrigation-network problems, the εHALP policies (ε = 1/16) are already close to optimal, and therefore hard to improve. Even if the MCMC solver reaches higher objective values than the ε-HALP solver, the policies may not improve due to the suggestive relationship between our true objective minw ‖V ∗ - V w‖∞,1/L and the objective of HALP minw ‖V ∗ - V w‖1,ψ [Guestrin et al., 2004].

Finally, the computation time of the ε-HALP solver is seriously affected by the topologies of tested networks, which can be explained as follows. For a small ε and large n, the time complexity of formulating cost networks grows approximately by the rates of (1/ε + 1)2 and (1/ε + 1)3 for the ring and ring-of-rings topologies, respectively. The ε-HALP solver spends a significant amount of time by formulating cost networks, which makes its decent time complexity on the ring topology (quadratic in 1/ε + 1) deteriorate on the ring-of-rings topology (cubic in 1/ε + 1). A similar crosstopology comparison of the MCMC solver shows that its computation times differ only by a multiplicative factor of 2. This difference is due to the increased complexity of sampling p(z∗ i | z-i), which is caused by more complex local dependencies, and not the treewidth of the ring-of-rings topology.

6 Conclusions

Development of scalable algorithms for solving large factored MDPs is a challenging task. The MCMC approach presented in this paper is a small but important step in this direction. In particular, our method overcomes the limitations of existing approaches to solving HALP and works directly with continuous variables, generates constraints based on their potential to improve the existing solution, and its space complexity is linear in the number of variables. Moreover, the MCMC

solver seems to be less affected by larger treewidth than the ε-HALP method while delivering substantially better results than uniform Monte Carlo sampling. Empirical results on two large control problems confirm the expected benefits of the approach and its potential to tackle complex real-world optimization problems. The objective of our future research is to eliminate the assumptions placed on the transition model and basis functions, which would make the framework applicable to a broader class of problems.

Acknowledgment

This work was supported in part by National Science Foundation grants CMS-0416754 and ANI-0325353. The first author was supported by an Andrew Mellon Predoctoral Fellowship for the academic year 2004-05. We thank anonymous reviewers for providing insightful comments that led to the improvement of the paper.