Max-norm Projections for Factored MDPs

1 Introduction

Over the last few years, Markov Decision Processes (MDPs) have been used as the basic semantics for optimal planning for decision theoretic agents in stochastic environments. In the MDP framework, the system is modeled via a set of states which evolve stochastically. The key problem with this representation is that, in virtually any real-life domain, the state space is quite large. However, many large MDPs have significant internal structure, and can be modeled compactly if the structure is exploited in the representation.

Factored MDPs [Boutilier et al., 1999] are one approach to representing large, structured MDPs compactly. In this framework, a state is implicitly described by an assignment to some set of state variables. A dynamic Bayesian network (DBN) [Dean and Kanazawa, 1989] can then allow a compact representation of the transition model, by exploiting the fact that the transition of a variable often depends only on a small number of other variables. Furthermore, the momentary rewards can often also be decomposed as a sum of rewards related to individual variables or small clusters of variables.

Even when a large MDP can be represented compactly, e.g., in a factored way, solving it exactly is still intractable: Exact MDP solution algorithms require the manipulation of a value function, whose representation is linear in the number of states, which is exponential in the number of state variables. One approach is to approximate the solution using an approximate value function with a compact representation. A

Ronald Parr

Computer Science Dept. Duke University parr@cs.duke.edu

common choice is the use of linear value functions as an approximation- value functions that are a linear combination of basis functions.

This paper makes a twofold contribution. First, we provide a new approach for approximately solving MDPs using a linear value function. Previous approaches to linear function approximation typically have utilized a least squares (L 2 -norm) approximation to the value function. Least squares approximations are incompatible withmost convergence analyses for MDPs, which are based on max-norm. We provide the first MDP solution algorithms - both value iteration and policy iteration - that use a linear max-norm projection to approximate the value function, thereby directly optimizing the quantity required to obtain the best error bounds.

Second, we show how to exploit the structure of the problem in order to apply this technique to factored MDPs. Our work builds on the ideas of Koller and Parr [1999; 2000], by using factored (linear) value functions, where each basis function is restricted to some small subset of the domain variables. We show that, for a factored MDP and factored value functions, various key operations can be implemented in closed form without enumerating the entire state space. Thus, our max-norm algorithms can be implemented efficiently, even though the size of the state space grows exponentially in the number of variables.

2 Markov decision processes

A Markov Decision Process (MDP) is defined as a 4-tuple (S;;A;;R;; P) where: S is a finite set of jSj states; A is a set of actions; R is a reward function R : S - A 7! IR, such that R(s;; a) represents the reward obtained by the agent in state s after taking action a; and P is a Markovian transition model where P(s 0 j s;; a) represents the probability of going from state s to state s 0 with action a.

We will be assuming that the MDP has an infinite horizon and that future rewards are discounted exponentially with a discount factor 2 [0;; 1). A stationary policy  for an MDP is a mapping  : S 7! A, where (s) is the action the agent takes at state s. The policy is associated with a value function V  2 IR jSj , where V  (s) is the discounted cumulative value that the agent gets if it starts at state s. The value function for a fixed policy is the fixed point of a set of equations that define the value of a state in terms of the value of its possible successor states. More formally, we define:

Definition 2.1 The DP operator, T  , for a fixed stationary policy  is:

T  V(s) = R(s;; (s)) +  s 0 P(s 0 js;; (s))V(s 0 ): V  is the fixed point of T  : V  = T  V  .

The optimal value function V  is also defined by a set of equations. In this case, the value of a state must be the maximal value achievable by any action at that state. More precisely, we define:

For any value function V, we can define the policy obtained by acting greedily relative to V. In other words, at each state, we take the action that maximizes the one-step utility, assuming that V represents our long-termutility achieved at the next state. More precisely, we define

The greedy policy relative to the optimal value function V  is the optimal policy   = Greedy(V  ). There are several algorithms to compute the optimal policy, we will focus on the two most used: value iteration and policy iteration.

Value iteration relies on the fact that the Bellman operator is a contraction - it is guaranteed to reduce the max-norm (L 1 ) distance between any pair of value functions by a factor of at least . This property guarantees that the Bellman operator has a unique fixed point V  [Puterman, 1994]. Value iteration exploits this property, approaching the fixed point through successive applications of the Bellman operator: V (t+1) = T  V (t) . After a finite number of iterations, the greedy policy Greedy(V (t) ) will be the optimal policy.

Policy iteration iterates over policies. Each iteration consists of two phases. Value determination computes, for a policy  (t) , the value function V  (t) , by finding the fixed point of: T  (t) V  (t) = V  (t) . Policy improvement defines the next policy as  (t+1) (s) = Greedy(V  (t) ). It can be shown that this process converges to the optimal policy.

3 SolvingMDPs with max-normprojections

In many domains, the state space is very large, and we need to perform our computations using approximate value functions. A very popular choice is to approximate a value function using linear regression. Here, we define our space of allowable value functions V 2 H  IR jSj via a set of basis functions H = fh 1 ;; : : : ;; h k g. A linear value function over H is a function V that can be written as V(s) = P k j=1 w j h j (s) for some coefficients w = (w 1 ;; : : : ;; w k ) 0 . We define H to be the linear subspace of IR jSj spanned by the basis functions H. It is useful to define an jSj - k matrix A whose columns are the k basis functions, viewed as vectors. Our approximate value function is then represented by Aw.

Linear value functions: The idea of using linear value functions for dynamic programming was proposed, initially, by Bellman et al. [1963] and has been further explored recently [Tsitsiklis and Van Roy, 1996; Koller and Parr, 1999; 2000]. The basic idea is as follows: in the solution algorithms, whether value or policy iteration, we use only value functions within H. Whenever the algorithm takes a step that results in a value function V that is outside this space, we project the result back into the space by finding the value function within the space which is close to V. More precisely:

Definition 3.1 A projection operator  is a mapping  : IR jSj ! H.  is said to be a projection w.r.t. a norm k
k if: V = Aw  such that w  2 argmin w kAw ; Vk.

Unfortunately, these existing algorithms all suffer from a problem that we might call "norm incompatibility." When computing the projection, they all utilize the standard projection operator with respect to L 2 norm or a weighted L 2 norm. On the other hand, most of the convergence and error analyses for MDP algorithms utilize max-norm. This incompatibility has made it difficult to provide error guarantees.

In this section, we propose a new approach that addresses the issue of norm compatibility. Our key idea is the use of a projection operator in L 1 norm. This problem has been studied in the optimization literature as the problem of finding the Chebyshev solution to an overdetermined linear system of equations [Cheney, 1982]. The problem is defined as finding w  such that: w

We will use an algorithm due to Stiefel [1960], that solves this problem by linear programming:

For the solution (w  ;;   ) of this linear program, w  is the solution of Eq. (1) and  is the L 1 projection error. Note that this LP only has k + 1 variables. However, there are 2
jSj constraints, which makes it impractical for large state spaces. In the remainder of this section, we will discuss how this projection addresses the norm incompatibility problem. In Section 4, we will show that, in factored MDPs, all the 2
jSj constraints can be represented efficiently, leading to a tractable algorithm.

Approximate Value iteration: The basic idea of approximate value iteration is quite simple. We define an L 1 projection operator  1 that takes a value function V and finds w that minimizes kAw;Vk 1 . This is an instance of Eq. (1) and can be solved by Eq. (2). The algorithm alternates applications of the Bellman operator T  and projection steps  1 :

:

In standard value iteration, we only need to perform the first step. However, V (t+1) may not be in H, so we need to add the second step, the projection step, to the process. We can analyze this process, bounding the overall error between our approximate value function Aw (t) and the optimal value function V  . This analysis shows that the overall error depends on the single-step max-norm projection errors (t+1) = Aw (t+1) ; V (t+1) 1 . Thus, using the max-norm projection, we can minimize these projection errors directly, we omit the analysis for lack of space. Note that, in approximate value iteration, the error introduced by successive approximations may grow unboundedly. As we will show, this cannot happen in approximate policy iteration.

Approximate Policy iteration: As we discussed, policy iteration is composed of two steps: value determination and policy improvement. Our algorithm performs the policy improvement step exactly. In the value determination step, the value function is approximated through a linear combination of basis functions. Consider the value determination for a policy  (t) . Define R  (t) (s) = R(s;;  (t) (s)), and P  (t) (s 0 j s) = P(s 0 j s;; a =  (t) (s)). We can now rewrite the value determination step in terms of matrices and vectors. If we view V  (t) and R  (t) as jSj-vectors, and P  (t) as an jSj - jSj matrix, we have the equations: V  (t) = R  (t) + P  (t) V  (t) . This is a system of linear equations with one equation for each state, which can only be solved exactly for small jSj. Our goal is to provide an approximate solution, within H. More precisely, we want to find:

This minimization is another instance of an L 1 projection (Eq. (1)), where C = (A ; P  (t) A) and b = R  (t) , and can be solved using the linear program of Eq. (2). Thus, our approximate policy iteration alternates between two steps:

For the analysis, we define the projection error for the policy iteration case, i.e., the error resulting from the approximate value determination step: (t) = Aw (t) ; ; R  (t) + P  (t) Aw (t)
1 . Unlike in approximate value iteration, the error is bounded in this case:

Lemma 3.2 There exists a constant P <1such that P  (t) for all iterations t of the algorithm.

Finally, we define discounted accumulated projection error as (t) = (t) + (t;1) ;; (0) = 0. Lemma 3.2 implies that the accumulated error remains bounded: (t)  P (1; t ) 1; .

We can now bound the error in the value function resulting from approximate policy iteration:

Theorem 3.3 In the approximate policy iteration algorithm, the distance between our approximate value function at iteration t and the optimal value function is bounded by:

In words, the difference between our approximation at iteration t and the optimal value function is bounded by the sum of two terms. The first term is present in standard policy iteration and goes to zero exponentially fast. The second is the discounted accumulated projection error and, as the theorem shows, is bounded. This second term can be minimized by choosing w (t) as the one that minimizes Aw (t) ; ; R  (t) + P  (t) Aw (t)
1 . Therefore, by performing max-norm projections, we can make the bound on the theorem as tight as possible.

4 Solving factoredMDPs

4.1 FactoredMDPs

Our presentation of factored MDPs follows that of [Koller and Parr, 2000]. In a factored MDP, the set of states is described via a set of random variables X = fX 1 ;; : : : ;;X n g, where each X i takes on values in some finite domain Dom(X i ). A state x defines a value x i 2 Dom(X i ) for each variable X i . We define a state transition model  using a dynamic Bayesian network (DBN) [Dean and Kanazawa, 1989]. Let X i denote the variable X i at the current time and X 0 i the variable at the next step. The transition graph of a DBN is a two-layer directed acyclic graph G  whose nodes are fX 1 ;; : : : ;;X n ;;X 0 1 ;; : : : ;;X 0 n g. We denote the parents of X 0 i in the graph by Parents  (X 0 i ). For simplicity of exposition, we assume that Parents  (X 0 i )  X; i.e., all arcs in the DBN are between variables in consecutive time slices. (This assumption can be relaxed, but our algorithmbecomes somewhat more complex.) Each node X 0 i is associated with a conditional probability distribution (CPD) P  (X 0 i j Parents  (X 0 i )). The transition probability P  (x 0 j x) is then defined to be Q i P  (x 0 i j u i ), where u i is the value in x of the variables in Parents  (X 0 i ).

Consider, for example, the problem of optimizing the behavior of a system administrator maintaining a network of n computers. Each machine is connected to some subset of other machines. In one simple network, we might connect the machines in a ring, with machine i connected to machines i + 1 and i ; 1. (In this example, we assume addition and subtraction are performed modulo n). Each machine is associated with a binary random variable F i , representing whether it has failed. The parents of F 0 i are F i , F i;1 , F i+1 . The CPD of F 0 i is such that if F i = true, then F 0 i = true with high probability; i.e., failures tend to persist. If F i = false, then F 0 i is a noisy or of its two other parents; i.e., a failure in either of its neighbors can independently cause machine i to fail.

We can define the transition dynamics of an MDP by defining a separate DBN model  a = hG a ;; P a i for each action a. However, in many cases, different actions have very similar transition dynamics, only differing in their effect on some small set of variables. In particular, in many cases a variable has a default evolution model, which only changes if an action affects it directly [Boutilier et al., 1999]. Therefore, as in [Koller and Parr, 2000], we use the notion of a default transition model  d = hG d ;; P d i. For each action a, we define E ects[a]  X 0 to be the variables in the next state whose local probabilitymodel is different from  d , i.e., those variables X 0 i such that P a (X 0 i j Parents a (X 0 i )) 6= P d (X 0 i j Parents d (X 0 i )).

In our system administrator example, we have an action a i for rebooting each one of the machines, and a default action d for doing nothing. The transition model described above corresponds to the "do nothing" action, which is also the default transitionmodel. The transitionmodel for a i is different from d only in the transition model for the variable F 0 i , which is now F 0 i = true with some small fixed probability, regardless of the status of the neighboringmachines.

Finally, we need to provide a compact representation of the reward function. We assume that the reward function is factored additively into a set of localized reward functions, each of which only depends on a small set of variables.

Definition 4.1 A function f is restricted to a domain C  X if f : Dom(C) 7! IR. If f is restricted to Y and Y  Z, we will use f(z) as shorthand for f(y) where y is the part of the instantiation z that corresponds to variables in Y.

Let R 1 ;; : : : ;;R r be a set of functions, where each R i is restricted to variable cluster W i  fX 1 ;; : : : ;; X n g. The reward function for state x is defined to be P r i=1 R i (x) 2 IR. In our example, we might have a reward function associated with each machine i, which depends on F i .

Factorization may allow us to represent large MDPs very compactly. However, we must still address the problem of solving these MDPs. Solution algorithms rely on the ability to represent value functions and policies, the representations which requires the same number of parameters as the size of the space. One might be tempted to believe that factored transition dynamics and rewards would result in a factored value function, which can thereby be represented compactly. Unfortunately, even in factored MDPs, the value function rarely has any internal structure [Koller and Parr, 1999].

Koller and Parr [1999] suggest that there are many domains where our value function might be "close" to structured, i.e., well-approximated using a linear combination of functions each of which refers only to a small number of variables. More precisely, they define a value function to be a factored (linear) value function if it is a linear function over the basis h 1 ;; : : :;; h k , where each h i is restricted to some subset of variables C i . In our example, we might have basis functions whose domains are pairs of neighboring machines, e.g., one basis function which is an indicator function for each of the four combinations of values for the pair of failure variables.

As shown by Koller and Parr [1999; 2000], factored value functions provide the key to performing efficient computations over the exponential-sized state sets that we have in factored MDPs. The key insight is that restricted domain functions (including our basis functions) allow certain basic operations to be implemented very efficiently. In the remainder of this section, we will show that this key insight also applies in the context of our algorithm.

4.2 FactoredMax-normProjection

The key computational step in both of our algorithms is the solution of Eq. (1) using the linear program in Eq. (2). In our setting, the vectors b and Cw are vectors in IR jSj , where S is our state space. In the case of factored MDPs, our state space is a set of vectors x which are assignments to the state variables X = fX 1 ;; : : : ;; X n g. We can view both Cw and b as functions of these state variables, and hence also their difference. Thus, we can define a function F w (X 1 ;; : : : ;;X n ) such that F w (x i ) = ((Cw) i ; b i ). Note that we have executed a representation shift; we are viewing F w as a function of X, which is parameterized by w.

The size of the state space is exponential in the number of variables. Hence, our goal in this section is to optimize Eq. (1) without explicitly considering each of the exponentially many states. The key is to use the fact that F w has a factored representation. More precisely, Cw has the form P j w j f 0 j (Z j ), where Z j is a subset of X. For example, we might have f 0 1 (X 1 ;;X 2 ) which takes value 1 in states where X 1 = true and X 2 = false and 0 otherwise. Similarly, the vector b in our case is also a sum of restricted domain functions. Thus, we can express F w as a sum P j f w j (Z j ), where f w j may or may not depend on w. In the future, we sometimes drop the superscript w when it is clear from context.

We tackle the problem of a factored solution to the LP in Eq. (2) in two steps.

Maximizing over the state space: First, assume that b and Cw are given, and that our goal is simply to compute max i ((Cw) i ; b i ) = max x F(x), i.e., to find the state x over which F is maximized. Recall that F = P m j=1 f j (Z j ). We can maximize such a function F using non-serial dynamic programming [Bertele and Brioschi, 1972] or cost networks [Dechter, 1999]. The idea is virtually identical to variable elimination in a Bayesian network. We review this construction here, as it is a key component in our solution LP.

Our goal is to compute

j where Z j [x] is the instantiation of the variables in Z j in the assignment x. The key idea is, rather than summing all functions and then doing the maximization, we maximize over variables one at a time. When maximizing over x l , only summands involving x l participate in the maximization. For example, assume

We can first compute the maximum over x 4 ; the functions f 1 and f 2 are irrelevant, so we can push them out. We get

The result of the internal maximization depends on the values of x 2 ;; x 3 ; i.e., we can introduce a new function e 1 (X 2 ;;X 3 ) whose value at the point x 2 ;; x 3 is the value of the internal max expression. Our problem now reduces to computing

having one fewer variable. Next, we eliminate another variable, say X 3 , with the resulting expression reducing to:

The result at this point is a number, which is the desired maximum over x 1 ;; : : : ;; x 4 .

In general, the variable elimination algorithmmaintains a set F of functions, which initially ff 1 ;; : : : ;; f m g. The algorithm then repeats the following steps:

1. Select an uneliminated variable X l ;
2. Take all e 1 ;; : : : ;; e L 2 F whose domain contains X l .
3. Define a new function e = max x l P j e j and introduce it into F. The domain of e is [ L j=1 Dom[e j ] ; fX l g.

The computational cost of this algorithm is linear in the number of new "function values" introduced in the elimination process. More precisely, consider the computation of a new function e whose domain is Z. To compute this function, we need to compute jDom[Z]j different values. The cost of the algorithm is linear in the overall number of these values, introduced throughout the algorithm. As shown in [Dechter, 1999], this cost is exponential in the induced width of the undirected graph defined over the variables X 1 ;; : : : ;; X n , with an edge between X l and X m if they appear together in one of the original functions f j .

Factored LP: Now, consider our original problem of minimizing kCw ; bk 1 = max x jF w (x)j over w. As in Eq. (2), we want to construct a linear programwhich performs this optimization. However, we want a compact LP, that avoids an explicit enumeration of the constraints for the exponentially many states. The first key insight is that we can replace the entire set of constraints -   (Cw) i ; b i for all states i - by the equivalent constraint   max i ((Cw) i ;b i ), or equivalently,   max x F w (x). The second key insight is that this new constraint can be implemented using a construction that follows the structure of variable elimination in cost networks. (An identical construction applies to the complementary constraints:   b i ; (Cw) i .)

Consider any function e used within F (including the original f i 's), and let Z be its domain. For any assignment z to Z, we introduce a variable into the linear program whose value represents u e z . For the initial functions f w i , we include the constraint that u f i z = f w i (z). As f w i is linear in w, this constraint is linear in the LP variables. Now, consider a new function e introduced into F by eliminating a variable X l . Let e 1 ;; : : : ;; e L be the functions extracted from F, and let Z be the domain of the resulting e. We introduce a set of constraints:

where Z j is the domain of e j and (z;; x l )[Z j ] denotes the value of the instantiation (z;; x l ) restricted to Z j . Let e n be the last function generated in the elimination, and recall that its domain is empty. Hence, we have only a single variable u e n . We introduce the additional constraint   u e n .

To understand this construction, consider our simple example above, and assume we want to express the fact that  max x F w (x). We first introduce a set of variables u f 1 x 1 ;;x 2 for every instantiation of values x 1 ;; x 2 to the variables X 1 ;;X 2 . Thus, if X 1 and X 2 are both binary, we would have four such variables. We would then introduce a constraint defining the value of u f 1 x 1 ;;x 2 appropriately. For example, for our f 1 above, we would have u f 1 t;;t = 0 and u f 1 t;;f = w 1 . We would have similar variables and constraints for each f j and each value z in Z j . Note that each of the constraints is a simple equality constraint involving numerical constants and perhaps the weight variables w.

Next, we introduce variables for each of the intermediate expressions. For example, we would have a set of variables u e 1 x 2 ;;x 3 ; for each of them, we would have a set of constraints

one for each value x 4 of X 4 . We would have a similar set of constraint for u e 2 x 1 ;;x 2 in terms of u f 2 x 1 ;;x 3 and u e 1 x 2 ;;x 3 . Note that each constraint is a simple linear inequality.

It is easy to show that minimizing  "drives down" the value of each variable u e z , so that

We can then prove, by induction, that u e n must be equal to max x P j f w j (x). Our constraints on  ensure that it is greater than this value, which is the maximum of P j f w j (x) over the entire state space. The LP, subject to those constraints, will minimize , guaranteeing that we find the vector w that achieves the lowest value for this expression.

Returning to our original formulation, we have that P j f w j is Cw ; b in one set of constraints and b; Cw in the other. Hence our new set of constraints is equivalent to the original set:   Cw ; b and   b ; Cw. Minimizing  finds the w that minimizes the L 1 norm, as required.

4.3 Factored solution algorithms

The factored max-norm projection algorithm described previously is the key to applying our max-norm solution algorithms in the context of factored MDPs.

Approximate value iteration: Let us begin by considering the value iteration algorithm. As described above, the algorithm repeatedly applies two steps. It first applies the Bellman operator to obtain V (t+1) = T  Aw (t) . Let  (t) be the stationary policy  (t) = Greedy(Aw (t) ). Note that  (t) corresponds to the maximizing policy used in the Bellman operator at iteration t, i.e., T  (t) = T  for iteration t. Thus, we can compute V (t+1) = T  Aw (t) by first computing  (t) and then performing a backprojection operation for this fixed policy, i.e., V (t+1) = R  (t) + P  (t) Aw (t) . Assume, for the moment, that P  (t) is a factored transitionmodel; we discuss the computation of the greedy policy and the resulting transition model below. As discussed by Koller and Parr [1999], the backprojection operation can be performed efficiently if the transition model and the value function are both factored appropriately. Furthermore, the resulting function V (t+1) is also factored, although the factors involve larger domains.

To recap their construction briefly, let f be a restricted domain function with domain Y; our goal is to compute P  f. We define the back-projection of Y through  as the set of parents of Y 0 in the transition graph G  ; ;  (Y 0 ) = [ Y 0 2Y 0 Parents  (Y 0 ). It is easy to show that: (P  f)(x) = P y 0 P  (y 0 j z)f(y 0 ), where z is the value of ;  (Y) in x. Thus, we see that (P  f) is a function whose domain is restricted to ;  (Y). Note that the cost of the computation depends linearly on jDom(;  (Y))j, which depends on Y (the domain of f) and on the complexity of the process dynamics.

Therefore, V (t+1) is composed of the sum of two factored functions: the reward function R  (t) which is assumed to be factored, and the backprojected basis functions P  (t) A which are also factored, as we have just shown. Finally, the second step in approximate value iteration is to compute the projection Aw (t+1) =  1 V (t+1) , i.e., find w (t+1) that minimizes Aw (t+1) ;V (t+1) 1 . As both V (t+1) and Aw (t+1) are

factored, we can perform this computation using the factored LP discussed in the previous section.

Approximate policy iteration: Policy iteration also iterates through two steps. The policy improvement step simply computes the greedy policy relative to  (t) . We discuss this step below. The approximate value determination step computes:

:

Again, assuming that P  (t) is factored, we can conclude that C = (A; P  (t) A) is also a matrix whose columns correspond to restricted-domain functions. The target b = R  (t) corresponds to the reward function, which is assumed to be factored. Thus, we can again apply our factored LP.

In our discussion so far, we assumed that we have some mechanism for computing the greedy policy Greedy(Aw (t) ), and that this policy has a compact representation and a factored transitionmodel. As shown in [Koller and Parr, 2000], the greedy policy relative to a factored value function has the form of a decision list. More precisely, the policy can be written in the form ht 1 ;; a 1 i;; ht 2 ;; a 2 i;; : : : ;; ht L ;; a L i, where each t i is an assignment of values to some small subset T i of variables, and each a i is an action. The optimal action to take in state x is the action a j corresponding to the first event t j in the list with which x is consistent.

Koller and Parr show the greedy policy can be represented compactly using such a decision list, and provide an efficient algorithm for computing it. Unfortunately, as they discuss, the resulting transition model is usually not factored. Thus, we cannot simply apply our factored LP construction, as suggested above. However, we can adapt it to deal with this issue.

The basic idea is to introduce cost networks corresponding to each branch in the decision list. Let S i be the set of states x for which t i is the first event in the decision list for which x is consistent. Recall that our LP construction defines a set of constraints that imply that   ((Cw) i ; b i ) for each state i. Instead, we will have a separate set of constraints for the states in each subset S i . For each state in S i , we know that action a i is taken. Hence, we can apply our construction above using P a i - a transition model which is factored by assumption - in place of the non-factored P  (t) .

The only issue is to guarantee that the cost network constraints derived from this transition model are applied only to states in S i . Specifically, we must guarantee that they are applied only to states consistent with t i , but not to states that are consistent with some t j for j < i. To guarantee the first condition, we simply instantiate the variables in T i to take the values specified in t i . That is, our cost network now considers only the variables in fX 1 ;; : : : ;;X n g ;T i , and computes the maximum only over the states consistent with T i = t i . To guarantee the second condition, we ensure that we do not impose any constraints on states associated with previous decisions. This is achieved by adding indicators I j for each previous decision t j , with weight ;1. These will cause the constraints associated with t i to be trivially satisfied by states in S j for j < i. Note that each of these indicators is a restricted domain function of T j and can be handled in the same fashion as all other terms in the factored LP. Thus, for a decision list of size L, our factored LP contains constraints from 2L cost networks.

It is instructive to compare our max-norm policy iteration algorithm with L 2 -projection policy iteration algorithm of Koller and Parr [2000] in terms of computational costs per iteration and implementation complexity. Computing the L 2 projection requires (among other things) a series of dot product operations between basis functions and backprojected basis functions hh i  P  h j i. These expressions are easy to compute if P  refers to the transition model of a particular action a. However, if the policy  is represented as a decision list, as is the result of the policy improvement step, then this step becomes much more complicated. In particular, for every point in the decision list, for every pair of basis functions i and j, and for each assignment to the variables in Dom(h i )[Dom(P a h j ), it requires the solution of a counting problemwhich is ]P-complete in general. Although, Koller and Parr show that this computation can be performed using a Bayesian Network (BN) inference, the algorithm still requires a BN inference for each one of those assignments at each point of the decision list. This makes the algorithmvery difficult to implement efficiently in practice.

The max-norm projection, on the other hand, relies on solving a linear program at every iteration. The size of the linear program depends on the cost networks generated. As we discuss, two cost networks are needed for each point in the decision list. The complexity of each of these cost networks is approximately the same as one of the BN inferences in the counting problem for the L 2 projection. Overall, for each point in the decision list, we have a total of two of these "inferences", as opposed to one for each assignment of Dom(h i ) [Dom(P a h j ) for every pair of basis functions i and j. Thus, the max-norm policy iteration algorithm is substantially less complex computationally than the approach based on L 2 -projection. Furthermore, the use of linear programming allows us to rely on existing LP packages (such as CPLEX), which are very highly optimized.

5 Experimental Results

The factored representation of a value function is most appropriate in certain types of systems: Systems that involvemany variables, but where the strong interactions between the variables are fairly sparse, so that the decoupling of the influence between variables does not induce an unacceptable loss in accuracy. As argued by Herbert Simon [1981] in "Architecture of Complexity", many complex systems have a "nearly decomposable, hierarchic structure", with the subsystems of such systems interacting only weakly between them. We selected to try our algorithm on a problem that we believe characterizes this type of structure.

The problem relates to a system administrator who has to maintain a network of computers; we experimented with various network architectures, shown in Fig. 1. Machines fail randomly, and a faulty machine increases the probability that its neighboring machines will fail. At every time step, the SysAdmin can go to one machine and reboot it, causing it to be working in the next time step with high probability. Each machine receives a reward of 1 when working (except in the ring, where one machine receives a reward of 2, to introduce some asymmetry), a zero reward is given to faulty machines, and the discount factor is = 0:95. Note that the additive structure of the reward function makes it unsuitable for the tree-structured representation used by Boutilier et al. [1999].

The basis functions included independent indicators for each machine, with value 1 if it is working and zero otherwise (each, a restricted domain function of a single variable), and the constant basis, whose value is 1 for all states.

We implemented the factored policy and value iteration algorithms in Matlab, using CPLEX as the LP solver. In our experiments, most of the time was spent inMatlab generating the LP constraints; CPLEX is a remarkably efficient and reliable solver, allowing even large LPs to be solved very quickly. We present only results for policy iteration. The time per iteration was about equal for policy and value iteration, but policy iteration converged in many fewer iterations, only about 4 or 5 iterations in the models we tested.

To evaluate the complexity of the algorithm, tests were performed with increasing number of states, that is, increasing number of machines on the network. Fig. 2 shows the running time for increasing problem sizes, for various architectures. The simplest one is the "Star", where the backprojection of each basis function has domain restricted to two variables and the largest factor in the cost network has domain restricted to two variables. The most difficult one was the "Bidirectional Ring", where factors contain five variables.

Note, the number of states is growing exponentially (indicated by the log scale in Fig. 2), but running times increase only logarithmically in the number of states, i.e., polynomially in the number of variables. This is illustrated by Fig. 3(a), where we fit a 3rd order polynomial to the running times for the "unidirectional ring". Note, the problem size grows quadratically with the number of variables: adding a machine to the network also adds the possible action of fixing that machine. For this problem, the computation cost of our algorithm empirically grows approximately as O ; (jXj
jAj) 1:5
.

For further evaluation, we measured the error in our approximate value function relative to the true optimal value function V  . Note that it is only possible to compute V  for small problems; in our case, we were only able to go up to 10 machines. For comparison, we also evaluated the error in the approximate value function produced by the L 2 -projection algorithm of Koller and Parr. As we discussed above, the L 2 projections in factored MDPs by Koller and Parr [2000] are difficult and time consuming; hence, we were only able to compare the two algorithms for smaller problems, where an equivalent L 2 -projection can be implemented using an explicit state space formulation. Results for both algorithms are presented in Fig. 3(b), showing the relative error of the approximate solutions to the true value function for increasing problem sizes. The results indicate that, for larger problems, the max-norm formulation generates a better approximation of the true optimal value function V  than the L 2 -projection. Here, we used two types of basis functions: the same single variable basis, and pairwise basis, which also includes indicators for neighboring pairs of machines. As expected, pairwise basis generated better approximations.

For these small problems, we can also compare the actual value of the policy generated by our algorithm to the value of the optimal policy. Here, the value of the policy generated by our algorithm is much closer to the value of the optimal policy than the error implied by the difference between our approximate value function and V  . For example, for the "Star" architecture with one server and up to 6 clients, our approximation with single variable basis functions had relative error of 12%, but the policy we generated had the same value as the optimal policy. In this case, the same was true for the policy generated by the L 2 projection. In an "Unidirectional Ring" with 8 machines and pairwise basis, the relative error between our approximation and V  was about 10%, but the resulting policy only had a 6% loss over the optimal policy. For the same problem, the L 2 approximation has a value function error of 12%, and a true policy loss was 9%. In other words, both methods induce policies that have lower errors than the errors in the approximate value function (at least for small problems). However, our algorithm continues to outperform the L 2 algorithm, even relative to actual policy loss.

For the large models, we can no longer compute the correct value function, so we cannot evaluate our results by computing kV  ;Awk 1 . However, the Bellman error, defined as BellmanErr(V) = kT  V ; Vk 1 , can be used to provide a bound: kV  ;Awk 1  BellmanErr(Aw) 1; [Williams and Baird, 1993]. Thus, we use the Bellman error to evaluate our answers for larger models. Fig. 3(c) shows that the Bellman error increases very slowly with the number of states. 1

Finally, we can look at the actual policies generated in our experiments. First, we noted that these tended to be short, the length of the final decision list policy grew approximately linearly with the number of machines. Furthermore, the policy itself is often fairly intuitive. In the "Ring and Star" architecture, for example, the decision list says: If the server is faulty, fix the server; else, if another machine is faulty, fix it.

6 Conclusions

In this paper, we presented new algorithms for approximate value and policy iteration. Unlike previous approaches, our algorithms directly minimize the L 1 error, and therefore have better theoretical performance guarantees than algorithms that optimize the L 2 -norm.

We have shown that the L 1 error can be minimized using linear programming, and provided an approach for representing this LP compactly for factored MDPs, allowing these algorithms to be applied efficiently even for MDPs with exponentially large state spaces. Our algorithms are more efficient and substantially easier to implement than previous algorithms based on the L 2 -projection.

We have presented results on a structured problem, representing a simplified version of a maintenance task. Our results show that our methods scale effectively to very large factored MDPs, and that our approach can exploit problemspecific structure to efficiently generate approximate solutions to complex problems.

The success of our algorithm depends on our ability to capture the most important structure in the value function using a linear factored approximation. This ability, in turn, depends on the choice of the basis functions and on the properties of the domain. The algorithm currently requires the designer to specify the factored basis functions. This is a limitation compared to other algorithms (e.g., [Dearden and Boutilier, 1997]), which are fully automated. However, our experiments indicate that a few simple rules are often quite successful in designing a basis. First, we ensure that the reward function is representable by our basis. A simple basis that, in addition, contains a separate set of indicators for each variable is often successful. We can also add indicators over pairs of variables; most simply, we can choose these according to the DBN transition model, where an indicator is added between variables X i and each one of the variables in Parents(X i ), thus representing one-step influences. This procedure can be extended, adding more basis to represent more influences as required. Thus, the structure of the DBN gives us indications of how to choose the basis functions. Other sources of prior knowledge can also be included for further specifying the basis.

1 Note that computing the Bellman error involves a maximization over the state space. Thus, the complexity of this computation grows exponentially with the number of variables. However, this maximization can be performed by a cost network using the same construction as in our max-norm projection.

The quality of our approximation also depends strongly on the structure of the domain. As we discussed above, our approximation is most successful in systems where variables are tightly coupled to only a small number of other variables. We believe that, for systems of this type, an approximation as a factored linear value function can be very accurate.

There are many possible extensions to this work. In particular, we hope to extend our algorithms to utilize other types of structure in the representation. One interesting direction involves factored action models, where multiple actions are taken simultaneously. Another involves the use of asymmetries (context-specificity) in the value function, as in the work of Dearden and Boutilier [1997], providing a complementary source of structure to the factorization used in our work.

Acknowledgments: We are very grateful to Dirk Ormoneit and Uri Lerner for many useful discussions. This work was supported by the ONR under the MURI program "Decision Making Under Uncertainty" and by the Sloan Foundation.