\documentclass[twoside,11pt]{article}
\usepackage{jair}
\usepackage{theapa-italics} 
\usepackage{codesection}
%\documentstyle[jair,twoside,11pt,theapa,psfig]{article}
\input epsf %% at top of file
\usepackage{dsfont}


\jairheading{19}{2003}{399-468}{1/02}{10/03}
\ShortHeadings{Efficient Solution Algorithms for Factored MDPs}
{Guestrin, Koller, Parr \& Venkataraman}
\firstpageno{399}

%\documentclass{article}
%\usepackage{times}
%\usepackage{ijcai01-mod}
\usepackage{epsfig}
\usepackage{graphics}
%\usepackage{latexsyms}
%\documentclass{amsart}
%\usepackage[active]{srcltx} % SRC Specials: DVI [Inverse] Search
% ----------------------------------------------------------------
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small
\hfuzz2pt % Don't report over-full h-boxes if over-edge is small

% ENVIRONMENTS -------------------------------------------------------
\newsavebox{\savepar}
\newenvironment{boxedtext}
{\noindent
\begin{lrbox}{\savepar}\begin{minipage}{\textwidth}}
{\end{minipage}\end{lrbox}\fbox{\usebox{\savepar}}}

\newenvironment{pseudocode}{\begin{boxedtext} \scshape}{\end{boxedtext}}



\newcommand{\bcsone}{\begin{codesection1}}
\newcommand{\bcstwo}{\begin{codesection2}}
\newcommand{\bcsthree}{\begin{codesection3}}
\newcommand{\bcsfour}{\begin{codesection4}}
\newcommand{\bcsfive}{\begin{codesection5}}
\newcommand{\bcssix}{\begin{codesection6}}
\newcommand{\ecsone}{\end{codesection1}}
\newcommand{\ecstwo}{\end{codesection2}}
\newcommand{\ecsthree}{\end{codesection3}}
\newcommand{\ecsfour}{\end{codesection4}}
\newcommand{\ecsfive}{\end{codesection5}}
\newcommand{\ecssix}{\end{codesection6}}
\newcommand{\codecomment}[1]{{\footnotesize \it //#1}}


% TEXTCOMMANDS -------------------------------------------------------
\newcommand{\tab}{\hspace*{1cm}}

% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
%\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
%\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
%\theoremstyle{assumption}
\newtheorem{assumpt}[thm]{Assumption}
%\numberwithin{equation}{section}
\newtheorem{example}[thm]{Example}

% MATH -----------------------------------------------------------

\newcommand{\bX}{{\bf X}}
\newcommand{\bY}{{\bf Y}}
\newcommand{\bZ}{{\bf Z}}
\newcommand{\bW}{{\bf W}}
\newcommand{\bC}{{\bf C}}
\newcommand{\bB}{{\bf B}}
\newcommand{\bA}{{\bf A}}
\newcommand{\bx}{{\bf x}}
\newcommand{\ba}{{\bf a}}
\newcommand{\bb}{{\bf b}}
\newcommand{\bc}{{\bf c}}
\newcommand{\bv}{{\bf v}}
\newcommand{\bw}{{\bf w}}
\newcommand{\bu}{{\bf u}}
\newcommand{\by}{{\bf y}}
\newcommand{\bz}{{\bf z}}
\newcommand{\bq}{{\bf q}}
\newcommand{\bt}{{\bf t}}
\newcommand{\vecb}{{\bf b}}
\newcommand{\bT}{{\bf T}}
\newcommand{\event}{\bt}
\newcommand{\Event}{\bT}
\newcommand{\vbp}{\bz}
\newcommand{\Scope}{{\rm Scope}}

\newcommand{\proof}{\noindent {\bf Proof: }}
\newcommand{\proofsketch}{\noindent {\bf Proof sketch: }}
\newcommand{\smilyface}[1]{\begin{picture}(10,10)(-.5,0)\setlength{\unitlength}{#1pt}\put(3.4,2){\circle{7.6}}\put(0.2,-1.0){\tiny$\smile$}
\put(1.8,3){\circle{.5}}\put(5.0,3){\circle{.5}}\end{picture}}
%\newcommand{\smilyface}[1]{\begin{picture}(10,10)(-.5,0)\setlength{\unitlength}{#1pt}\put(4,2){\circle{9.5}}\put(0.5,-1.5){\scriptsize$\smile$}
%\put(2.1,4){\circle{.5}}\put(6,4){\circle{.5}}\end{picture}}
%\newcommand{\smilyface}[1]{\begin{picture}(10,10)(-.5,0)\setlength{\unitlength}{#1 pt}\put(4,2){\circle{8.5}}\put(0.1,-1.5){\scriptsize$\smile$} \put(1.1,1.7){\scriptsize$\cdot\hspace{0.5ex}\cdot$}\end{picture}}
%\newcommand{\smilyface}[1]{\scriptsize{\colorbox{black}{\makebox[-0.6ex][l]{\raisebox{-.5ex}{\color{white}$\smile$}}
%\raisebox{.5ex}{\color{white}$\cdot\hspace{0.5ex}\cdot$}}}}
\newcommand{\bbox}{\vrule height7pt width4pt depth1pt}
%\newcommand{\qed}{$\ \ \ \bbox$}
%\newcommand{\qedmath}{\mbox{\qed}}
\def\boxit#1{\vbox{\hrule\hbox{\vrule\kern3pt
   \vbox{\kern3pt#1\kern3pt}\kern3pt\vrule}\hrule}}
\def\Box{\boxit{\null}}
\newcommand{\qedmath}{{\ \ \Box}}
\newcommand{\qed}{{$\qedmath$}}
\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\set}[1]{\left\{#1\right\}}
%\newcommand{\Real}{\makebox[.50ex][l]{I}\makebox{R}}
\newcommand{\Real}{\mathds{R}}
\newcommand{\Natural}{\mathds{N}}
\newcommand{\eps}{\varepsilon}
\newcommand{\epstilde}{\tilde{\varepsilon}}
\newcommand{\epsp}{\eps_P}
\newcommand{\cL}{\mathcal{L}}
\newcommand{\Linf}{\cL_\infty}
\newcommand{\wstar}{w^\ast}
\newcommand{\what}{\widehat{w}}
\newcommand{\phistar}{\phi^\ast}
\newcommand{\tstar}{t^\ast}
\newcommand{\bwstar}{\bw^\ast}
\newcommand{\bwhat}{\widehat{\bw}}
\newcommand{\V}{\mathcal{V}}
\newcommand{\Vstar}{\V^\ast}
\newcommand{\Vbar}{\overline{\V}}
\newcommand{\maxnorm}[1]{\norm{#1}_\infty}
\newcommand{\T}{\mathcal{T}}
\newcommand{\Tbell}{\mathcal{T}^*}
%\newcommand{\Errpi}[1]{Err_\Pi(#1)}
\newcommand{\Piinf}{{\Pi_{\mbox{\tiny{$\infty$}}}}}
\newcommand{\Proj}{{\Pi}}
\newcommand{\phibar}{\overline{\phi}}
\newcommand{\betabar}{\overline{\beta}}
\newcommand{\betap}{\beta_P}
\newcommand{\policy}{\pi}
\newcommand{\policystar}{\policy^\ast}
\newcommand{\policyhat}{\widehat{\policy}}
\newcommand{\basis}{h}
\newcommand{\Basis}{H}
\newcommand{\Basisspace}{\mathcal{H}}
\newcommand{\backbasis}{g}
\newcommand{\Backbasis}{G}
\newcommand{\Basismatrix}{{\rm H}}
\newcommand{\Cluster}{{\bf C}}
\newcommand{\Clusters}{{\cal C}}
\newcommand{\clustval}{\bf c}
\newcommand{\rclust}{{\bf U}}
\newcommand{\discount}{\gamma}
\newcommand{\Rmax}{R_{max}}
\newcommand{\influence}[1]{{\it Influences\/}[{#1}]}
\newcommand{\control}[1]{{\it Control\/}[{#1}]}
\newcommand{\at}[1]{^{(#1)}}
\newcommand{\Parents}{{\emph{Parents}}}
\newcommand{\effects}[1]{{{\it Effects\/}[{#1}]}}
\newcommand{\rewards}[1]{{{\it Rewards\/}[{#1}]}}
\newcommand{\graph}{G}
\newcommand{\trans}{\tau}
\newcommand{\backproj}[2]{\Gamma_{#1}({#2})}
\newcommand{\Val}{{\rm Dom}}
\newcommand{\numagents}{g}
\newcommand{\Bellman}{{\rm BellmanErr}}
\newcommand{\greedy}[1]{{\textsl{Greedy}}({#1})}
\newcommand{\sweights}{{\alpha}}
\newcommand{\relcomp}{I}
\newcommand{\diffun}{\delta}
\newcommand{\Diff}{{\bT}}



\newcommand{\indicator}{\mathcal{I}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\bF}{{\bf F}}
\newcommand{\bS}{{\bf S}}
\newcommand{\bs}{{\bf s}}
\newcommand{\bQ}{{\bf Q}}
\newcommand{\Children}{{\emph{Children}}}
\newcommand{\Parent}{{\emph{Parent}}}
\newcommand{\Clique}{\bQ}
\newcommand{\Cliquevalue}{\bq}
\newcommand{\numcliques}{{\abs{\Clique}}}
\newcommand{\separator}{\bS}
\newcommand{\separatorvalue}{\bs}
\newcommand{\sepmessage}[1]{m_{#1}}

\newcommand{\denselist}{\itemsep 0pt\topsep-8pt\partopsep-8pt}
\newcommand{\bitem}{\begin{itemize}\denselist}
\newcommand{\eitem}{\end{itemize}}
\newcommand{\benum}{\begin{enumerate}\denselist}
\newcommand{\eenum}{\end{enumerate}}
\newcommand{\beqa}{\vspace{0pt}\begin{eqnarray*}}
\newcommand{\eeqa}{\end{eqnarray*}}
\newcommand{\bdm}{\vspace{0pt}\begin{displaymath}}
\newcommand{\edm}{\end{displaymath}}
\newcommand{\beqn}{\vspace{0pt}\begin{equation}}
\newcommand{\eeqn}{\end{equation}}
\newcommand{\tableref}[1]{Table~\ref{#1}}
\newcommand{\figref}[1]{Figure~\ref{#1}}
\newcommand{\eqnref}[1]{Equation~(\ref{#1})}
\newcommand{\secref}[1]{Section~\ref{#1}}
\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\newcommand{\appenref}[1]{Appendix~\ref{#1}}

\newcommand{\true}{\textit{true}}
\newcommand{\false}{\textit{false}}
\newcommand{\ttrue}{\textit{t}}
\newcommand{\ffalse}{\textit{f}}

%%% Trying to save some space:
%\renewenvironment{abstract}{\centerline{\bf Abstract}\vspace{0.2ex}\begin{quote}}{\par\end{quote}\vskip 0.8ex}
%\abovedisplayskip 2pt plus2pt minus5pt%
%\belowdisplayskip \abovedisplayskip
%\abovedisplayshortskip  0pt plus3pt%
%\belowdisplayshortskip  2pt plus3pt minus3pt%

\newcommand{\cparagraph}{\vspace{0pt}\paragraph}
%
\newcommand{\eqn}{\begin{equation}}
\newcommand{\eqnar}{\begin{eqnarray}}
\newcommand{\dotprod}[2]{\langle #1 \bullet #2\rangle}
\newcommand{\captionsmall}[1]{\caption{\footnotesize{#1}}}

\newcommand{\bprojfun}[2]{{\it Backproj_{#1}}({#2})}
\newcommand{\order}{\mathcal{O}}
\newcommand{\constraints}{\Omega}
\newcommand{\dlist}{\Delta}
\newcommand{\dlisthat}{\widehat{\Delta}}
%\newcommand{\indicone}{\makebox[.50ex][l]{1}\makebox{I}}
%\newcommand{\indicone}{{\bf 1}}
\newcommand{\indicone}{\mathds{1}}



\newcommand{\fmbox}[1]{\mbox{\footnotesize #1}}
\newcommand{\apricodd} {\textsc{Apricodd}}
\newcommand{\cplex}{\textsc{Cplex }}
\newcommand{\agents}{a}
\newcommand{\valuerule}[1]{{\left\langle{#1}\right\rangle}}
\newcommand{\rrule}{\rho}
\newcommand{\prule}{\eta}
\newcommand{\lprule}{e}
\newcommand{\context}{\bc}
\newcommand{\Context}{\bC}
\newcommand{\vvalue}{v}
\newcommand{\pvalue}{p}
\newcommand{\lpvalue}{u}
\newcommand{\numrules}{n}
\newcommand{\relevant}[1]{{\it Agents\/}[{#1}]}
\newcommand{\observe}[1]{{\it Obs\/}[{#1}]}
\newcommand{\maxout}[1]{{\it MaxOut\/}\left({#1}\right)}
\newcommand{\rulebproj}[2]{{\it RuleBackproj_{#1}}({#2})}
\newcommand{\rulesplit}[1]{{\it Split\/}\left({#1}\right)}
\newcommand{\residual}[1]{{\it Residual\/}\left({#1}\right)}
\newcommand{\objective}{{\it Objective\/}}
\newcommand{\setP}{\mathcal{P}}
\newcommand{\var}{\it Var\/} 
\newcommand{\union}{{\cup}}
\newcommand{\defref}[1]{Definition~\ref{#1}}

\newcommand{\instin}[2]{{#2}[#1]}
\newcommand{\consist}[2]{#1\sim[#2]}

%CGjairf2 added these italics macros 
\newcommand{\ie}{{\it i.e.\/}}
\newcommand{\eg}{{\it e.g.\/}}
\newcommand{\etal}{\emph{et al.}}


%CGjairf2 moved all proofs to appendix

\newcommand{\commentout}[1]{}
% -------------------------------------------------------------
\begin{document}

\title{Efficient Solution Algorithms for Factored MDPs}

\author{\name Carlos Guestrin \email guestrin@cs.stanford.edu\\
\addr Computer Science Dept., Stanford University
 \AND
\name Daphne Koller \email koller@cs.stanford.edu\\
\addr Computer Science Dept., Stanford University
\AND
\name Ronald Parr \email parr@cs.duke.edu\\
\addr Computer Science Dept., Duke University
\AND
\name Shobha Venkataraman \email shobha@cs.cmu.edu\\
\addr Computer Science Dept., Carnegie Mellon University
}

%\thanks{}%
%\subjclass{}%
%\keywords{}%

%\date{}
%\dedicatory{}%
%\commby{}%
%

% For research notes, remove the comment character in the line below.
% \researchnote
\maketitle

\begin{abstract}
%CGjairf2 shorten a little and rewrote parts
%DK3 - rewrote to shorten considerably
\footnotesize
This paper addresses the problem of planning under uncertainty in large
Markov Decision Processes (MDPs).
%domains.
%RP5:  reworded
%  \emph{Factored} Markov Decision Processes (MDPs), where the
%transition model is represented as a dynamic Bayesian network (or
%influence diagram), provide a compact representation for
% structured
%domains.
%DK4
%Many large MDPs have structure that permits factorization into a compact
%state variable representation, where the transition model is represented as
%a dynamic Bayesian network (or influence diagram).  \emph{Factored} MDPs
%can enjoy an exponential reduction in representation size over an explicit
%representation.
\emph{Factored MDPs} represent a complex state space using state
variables and the transition model using a dynamic Bayesian network.
%DK4
%However, the complexity of exact solution algorithms for
%such MDPs can grow exponentially in the representation size.  
This representation often allows an exponential reduction in the
representation size of structured MDPs, but the complexity of exact
solution algorithms for such MDPs can grow exponentially in the
representation size.  
In this paper, we
present two approximate solution algorithms that exploit structure in
factored MDPs.  Both use an approximate value function represented as a
linear combination of basis functions, where each basis function involves
only a small subset of the domain variables.  
A key contribution of this paper is that it shows how the
basic operations of both algorithms can be performed efficiently in
closed form, by exploiting both \emph{additive} and
\emph{context-specific} structure in a factored MDP.  A central
element of our algorithms is a novel linear program decomposition
technique, analogous to variable elimination in Bayesian networks,
which reduces an exponentially large LP to a provably equivalent, polynomial-sized one. 
One algorithm uses
approximate linear programming, and the second approximate dynamic
programming.  Our dynamic programming algorithm is novel in that it uses
an approximation based on max-norm, a technique that more directly
minimizes the terms that appear in error bounds for approximate MDP
algorithms.  
We provide experimental results on problems
with over $10^{40}$ states, demonstrating a promising indication of the scalability of
our approach, 
%CGjairf4 added:
and compare our algorithm to an existing state-of-the-art approach, showing, in some problems, exponential gains in computation time. 
\end{abstract}


% --------------------------------------------------------------
% --------------------------------------------------------------
\section{Introduction}
%CGjairf2 rewrote introduction

Over the last few years, \emph{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
%RP3:  key -> main
stochastically. The main 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.

\emph{Factored MDPs}~\cite{Boutilier+al:2000} are one approach to
representing large, structured MDPs compactly.  In this framework, a state
is implicitly described by an assignment to some set of \emph{state
variables}. A \emph{dynamic Bayesian network
(DBN)\/}~\cite{Dean+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.

There are two main types of structure that can simultaneously be
exploited in factored MDPs: additive and context-specific
structure. \emph{Additive} structure captures the fact that typical
large-scale systems can often be decomposed into a combination of
locally interacting components. For example, consider the management
of a large factory with many production cells.  Of course, in the long
run, if a cell positioned early in the production line generates
faulty parts, then the whole factory may be affected. However, the
quality of the parts a cell generates depends \emph{directly} only on
the state of this cell and the quality of the parts it receives from
neighboring cells. Such additive structure can also be present in the
reward function. For example, the cost of running the factory depends,
among other things, on the \emph{sum} of the costs of maintaining each
local cell.  

%DK4
%The second type of structure exploited in factored MDPs addresses
%\emph{context-specific} independence. Here we tackle a different type
%of local influence
\emph{Context-specific} structure encodes a different type of locality
of influence: Although a part of a large system may, in general,
be influenced by the state of every other part of this system, at any
given point in time only a small number of parts may influence it
directly. In our factory example, a cell responsible for anodization
may receive parts directly from any other cell in the
factory. However, a work order for a cylindrical part may restrict
this dependency only to cells that have a lathe. Thus, \emph{in the
  context} of producing cylindrical parts, the quality of the anodized
parts depends directly \emph{only} on the state of cells with a lathe.  

Even when a large MDP can be represented compactly, for example, by using
a factored representation,
%RP3:  Exact -> typical exact
solving it exactly may still be intractable: Typical
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 common choice is the use of \emph{linear}
value functions as an approximation --- value functions that are a
linear combination of potentially non-linear basis
functions~\cite{Bellman+al:1963,Sutton:1988,Tsitsiklis+VanRoy:1996a}. 
Our work builds on the ideas of
Koller and Parr~\citeyear{Koller+Parr:1999,Koller+Parr:2000}, by using
\emph{factored (linear) value functions}, where each basis function is
restricted to some small subset of the domain variables.  


%DK4 - rewrote
%This paper presents two new algorithms for computing linear value
%function approximations for factored MDPs. These algorithms can
%exploit \emph{both} additive and context-specific structure in the 
%model, to produce approximate solutions very efficiently.  They differ from
%earlier approaches~\cite{Boutilier+al:2000}, which only exploit
%context-specific structure, and in their use of a linear value
%function, a more sophisticated and expressive function approximation
%method. We provide a more detailed discussion of the differences in
%\secref{sec:Related}.

This paper presents two new algorithms for computing linear value
function approximations for factored MDPs:  one that uses approximate
dynamic programming and another that uses approximate linear programming.   
Both algorithms are based on the use of factored linear value
functions, a highly expressive function approximation method.  This
representation allows the algorithms to take advantage of \emph{both}
additive and context-specific structure, in order to produce
high-quality approximate solutions very efficiently.  The capability
to exploit both types of structure distinguishes these algorithms
differ from earlier approaches~\cite{Boutilier+al:2000}, which only
exploit   context-specific structure.  We provide a more detailed
discussion of the differences in \secref{sec:Related}.

%DK4 - slight rewrite, and merged two paragraphs
We show that, for a factored MDP and factored value functions, various
critical operations for our planning algorithms can be implemented in
closed form without necessarily enumerating the entire state space.
%RPfinal:  Added 'upon'
In particular, both our new algorithms build upon a novel linear
programming decomposition technique. This technique reduces structured
LPs with exponentially many constraints to equivalent,
polynomially-sized ones. This decomposition follows a procedure
analogous to variable elimination that applies both to additively
structured value functions~\cite{Bertele+Brioschi:1972} and to value
functions that also exploit context-specific
structure~\cite{Zhang+Poole:1999}. 
%DK4
%, and uses linear constraints to represent the factors generated in
%the variable elimination procedure. 
Using these basic operations, our planning algorithms can be
implemented efficiently, even though the size of the state space grows 
exponentially in the number of variables.

Our first method is based on the approximate linear programming
algorithm~\cite{Schweitzer+Seidmann:1985}. This algorithm
generates a linear, approximate value function by solving a single
linear program. Unfortunately, the number of constraints in the LP
proposed by Schweitzer and Seidmann grows exponentially in the
number of variables. Using our LP decomposition technique, we exploit
structure in factored MDPs to represent exactly the same optimization 
problem with exponentially fewer constraints. 

In terms of approximate dynamic programming, 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 ($\cL_2$-norm) approximation to the value function.
Least squares approximations are incompatible with most 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 that appears in our
provided error bounds. Second, we show how to exploit the structure of
%RPfinal:  replace ``in order to'' with ``to'' here and elsewhere
%          to comply with recommendations of style guides.
%          Others are not marked.
the problem to apply this technique to factored MDPs, by
again leveraging on our LP decomposition technique.  
%DK4 - no need to belabor this quite so much
%As with approximate linear programming, our max-norm algorithms can be
%implemented efficiently, even though the size of the state space grows
%exponentially in the number of variables.


Although approximate dynamic programming currently possesses stronger
theoretical guarantees, our experimental results suggest that
approximate linear programming is a good alternative. Whereas the
former tends to generate better policies for the same set of basis
functions, due to the simplicity and computational advantages of
approximate linear programming, we can add more basis functions,
obtaining a better policy and still requiring less computation than
the approximate dynamic programming approach.

%CGjairf4 added:
Finally, we present experimental results comparing our approach to the
work of Boutilier \etal~\citeyear{Boutilier+al:2000}, illustrating
some of the tradeoffs between the two methods. In particular, for
problems with significant context-specific structure \emph{in the
  value function}, their approach can be faster due to their efficient
handling of their value function representation. However, there are
cases with significant context-specific structure in the problem,
rather than in the value function, in which their algorithm requires
%RPfinal:  Added changed to say classes of problems
an exponentially large value function representation. In such classes of
problems, we demonstrate that by 
%using a linear value function, \ie,
using a value function that exploits both additive and
context-specific structure, our algorithm can obtain a polynomial-time
near-optimal approximation of 
the true value function. 


\commentout{
The approximate dynamic programming with max-norm projection
algorithms presented here first appeared in Guestrin
\etal~\citeyear{Guestrin+al:2001b}. We first presented the approximate 
linear programming approach, along with extensions to multiagent
systems in Guestrin \etal~\citeyear{Guestrin+al:2001c}.
%RP3: Added
The extension to context specific independence was first presented in
\cite{Guestrin+al:2002d}.  In addition to combining these earlier
results, this paper provides more detailed proofs, analyses and
experimental evaluation.  A detailed description of the relationship to the
work of Boutilier \etal~\citeyear{Boutilier+al:2000} is provided, along
with an empirical comparison.
}

%DK4 - I think we are repeating ourselves far too much.  I found
%myself getting annoyed at the authors ;-)
%I cut this.
%CGjairf2 there will be a paragraph here on the experimental section
%RP5:  killed ``we feel that''
\commentout{
The main contributions of this paper include:
\begin{itemize}
\item A novel and general LP decomposition technique, which can yield an exponential reduction in the number of constraints in the problem.
\item Two new approximate planning algorithms for factored MDPs, which can exploit both additive and context-specific structure:
\begin{itemize}
\item Approximate linear programming for factored MDPs;
\item Approximate dynamic programming with max-norm projection.
\end{itemize}
\item A theoretical analysis of the use of max-norm projection as a tool for general approximate dynamic programming.
\item An efficient method for bounding the quality of the policies generated by factored MDP algorithms based on the Bellman error analysis of Williams and Baird~\citeyear{Williams+Baird:1993}.
%CGjairf2 I hope we can comment out this bullet....
%\item The first empirical comparison of factored linear value function-based algorithms against the tree-based decomposition technique of Boutilier \etal~\citeyear{Boutilier+al:2000}.
\end{itemize}
}

This paper starts with a presentation of factored MDPs and approximate
solution algorithms for MDPs. In \secref{sec: basic ops}, we describe
the basic operations used in our algorithms, including our LP
decomposition technique. 
%DK4 - rewrote 
In \secref{sec:ALP}, we present the first of our two algorithms: the
approximate linear programming algorithm for factored MDPs. The second
algorithm, approximate policy iteration with max-norm projection, is
presented in \secref{sec:API}.  \secref{sec: Bellman error} describes 
an approach for efficiently computing bounds on policy quality based
on the Bellman error.  \secref{sec: CSI} shows how to extend our
methods to deal with context-specific structure. Our paper concludes
with an empirical evaluation in \secref{sec:Experimental results} and
a discussion of related work in \secref{sec:Related}. 

\commentout{
The paper is divided in two main parts: In the first part, we show
the analysis and describe our new algorithms for the general MDP
case. In the second part, we first describe the representation of
a factored MDP. This formulation is well suited for algorithms
that exploit structure to solve large MDPs. In particular, we show
how this factored representation can be exploited to reduce the
complexity of our new algorithms. Finally, we present experimental
results demonstrating the running time of our algorithms and the
quality of our approximations.
}

%RPfinal:  I fixed the citation style here.
This paper is a greatly expanded version of work that was
published before in Guestrin \etal~\citeyear{Guestrin+al:2001b}, and some
of the work presented in Guestrin
\etal~\citeyear{Guestrin+al:2001c,Guestrin+al:2002d}.


% --------------------------------------------------------------
% --------------------------------------------------------------
\section{Factored Markov Decision Processes}
A Markov decision process (MDP) is a mathematical framework for sequential
decision problems in stochastic domains.  It thus provides an underlying
semantics for the task of planning under uncertainty.  We begin with a
concise overview of the MDP framework, and then describe the
representation of \emph{factored MDPs}. 

\subsection{Markov Decision Processes}
\label{sec: MDP}
We briefly review the MDP framework, referring the reader to 
%RP4: changed to citeyear
the books by Bertsekas and Tsitsiklis~\citeyear{Bertsekas+Tsitsiklis:1996}
or Puterman~\citeyear{Puterman:1994} for a more
in-depth review.  A \emph{Markov Decision Process (MDP)\/} is defined as a 
4-tuple $(\bX, A, R, P)$ where: $\bX$ is a finite set of
$\abs{\bX} = N$ states; $A$ is a finite set of actions; $R$ is a {\em
reward function\/} $R: \bX \times A \mapsto \Real$, such that
$R(\bx,a)$ represents the reward obtained by the agent in state
$\bx$ after taking action $a$; and $P$ is a \emph{Markovian
transition model\/} where $P(\bx' \mid \bx,a)$ represents the
probability of going from state $\bx$ to state $\bx'$ with action
$a$. We assume that the rewards are bounded, that is, there
exists $\Rmax$ such that $\Rmax \geq \abs{R(\bx,a)}, \forall
\bx,a$.

\begin{example}
Consider the problem of
optimizing the behavior of a system administrator (SysAdmin)
maintaining a network of $m$ computers.  
%CGjairf2 added
In this network,
each machine is connected
to some subset of the other machines. Various possible network
topologies can be defined in this manner (see \figref{fig:network
topologies} for some examples). 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 $m$.)

Each machine is associated with a binary random variable $X_i$,
representing whether it is working or has failed. At every time
step, the SysAdmin receives a certain amount of money (reward) for
each working machine.  The job of the SysAdmin is to decide which
machine to reboot; thus, there are $m+1$ possible actions at each
time step: reboot one of the $m$ machines or do nothing (only one
machine can be rebooted per time step). If a machine is rebooted,
it will be working with high probability at the next time step.
Every machine has a small probability of failing at each time
step. However, if a neighboring machine fails, this probability
increases dramatically. These failure probabilities define the
transition model $P(\bx' \mid \bx,a)$, where $\bx$ is a particular
assignment describing which machines are working or have failed in
the current time step, $a$ is the SysAdmin's choice of machine to
reboot and $\bx'$ is the resulting state in the next time step. \qed
\end{example}

\begin{figure}[t]
\centerline{\epsfig{figure=Figures/network-topologies.eps,height=5cm}}
\caption{Network topologies tested; the status of a machine is
influence by the status of its parent in the network.}
\label{fig:network topologies}
\end{figure}


We assume that the MDP has an infinite horizon and that future rewards
are discounted exponentially with a discount factor $\discount \in
[0,1)$. A \emph{stationary policy} $\policy$ for an MDP is a mapping
$\policy: \bX \mapsto A$, where $\policy(\bx)$ is the action the agent
takes at state $\bx$. In the computer network problem, for each possible
configuration of working and failing machines, the policy would tell the
SysAdmin which machine to reboot. 
%CGjairf2 The -> Each
Each policy is associated with a {\em
value function\/} $\V_{\policy} \in \Real^{N}$, where $\V_{\policy}(\bx)$
is the discounted cumulative value that the agent gets if it starts at
state $\bx$ and follows policy $\policy$. More precisely, the value
$\V_{\policy}$ of a state $\bx$ under policy $\policy$ is given by:
\bdm
\V_{\policy}(\bx) = E_{\policy}\left[ \left. \sum_{t=0}^{\infty} \discount^t R\left(\bX\at{t}, \policy(\bX\at{t}) \right) \right | \bX\at{0} = \bx  \right],
\edm
where $\bX\at{t}$ is a random variable representing the state of the system after $t$ steps.
In our running example, the value function represents how much money the
SysAdmin expects to collect if she starts acting according to $\policy$
when the network is at state $\bx$.  The value function for a fixed
policy is the fixed point of a set of 
%CGjairf2 added
%RP4:  remove parentheses (linear)
linear
equations that define the value of
a state in terms of the value of its possible successor states.  More
formally, we define: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: DP operator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RP3: Fixed as per reviewer 1's comments
\begin{defn}
%RP1
%The \emph{DP operator}, $\T_{\policy}$, for a fixed stationary
The \emph{DP operator}, $\T_{\policy}$, for a stationary
policy $\policy$ is:
\bdm
    \T_{\policy} \V(\bx) = {R(\bx,\policy(\bx)) + \discount \sum_{\bx'}
    P(\bx'\mid\bx,\policy(\bx))\V(\bx')}.
\end{displaymath}
The value function of policy $\policy$, $\V_{\policy}$, is the fixed point
of the $\T_{\policy}$ operator: 
$\V_{\policy} = \T_{\policy} \V_{\policy}$. \qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The optimal value function $\Vstar$ describes the optimal value the agent can achieve for each starting state. $\Vstar$ is also defined by a set of 
%CGjairf2 added
%RP4:  remove parentheses (non-linear)
non-linear
equations.
%RP3: clarified for reviewer
In this case, the value of a state must be the maximal expected value
achievable by any policy starting at that state.  More precisely, we define:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Bellman operator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
The \emph{Bellman operator}, $\Tbell$, is:
\bdm
    \Tbell \V(\bx) = \max_a[{R(\bx,a) + \discount \sum_{\bx'}
    P(\bx' \mid \bx,a)\V(\bx')}].
\end{displaymath}
The optimal value function $\Vstar$ is the fixed point of $\Tbell$: $ \Vstar = \Tbell \Vstar$.\qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


For any value function $\V$, we can define the policy obtained by
acting greedily relative to $\V$.  In other words, at each state, the agent takes the action that maximizes the one-step utility, assuming that $\V$ represents our long-term utility achieved at the next state. More precisely, we define:
\beqn \label{eq: greedy}
 \greedy{\V}(\bx) = \arg\max_a[R(\bx,a) + \discount \sum_{\bx'} P(\bx' \mid \bx,a) \V(\bx')].
\end{equation}
The greedy policy relative to the optimal value function $\Vstar$ is the
optimal policy $\policystar = \greedy{\Vstar}$. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\subsection{Factored MDPs} \label{sec: Factored MDPs}
\emph{Factored MDPs} are a representation language that allows us to
exploit problem structure to represent exponentially large MDPs very
compactly. The idea of representing a large MDP using a factored model
was first proposed by Boutilier {\it et
al.\/}~\citeyear{Boutilier+al:1995}. 
%Our presentation of factored MDPs follows that of Koller and
%Parr~\citeyear{Koller+Parr:2000}. 

In a factored MDP, the set of states is described via a set of random
variables $\bX = \{X_1,\ldots,X_n\}$, where each $X_i$ takes on values in
some finite domain $\Val(X_i)$.  A state $\bx$ defines a value $x_i \in
\Val(X_i)$ for each variable $X_i$.  In general, we use upper case
letters (\eg, $X$) to denote random variables, and lower case (\eg, $x$) to denote their values.  We use boldface to denote vectors of
variables (\eg, $\bX$) or their values ($\bx$).  For an instantiation
$\by \in \Val(\bY)$ and 
%CGjairf2 set -> subset
a subset of these variables $\bZ \subset \bY$, we use
$\instin{\bZ}{\by}$ to denote the value of the variables $\bZ$ in the
instantiation $\by$.

%CGjairf2 added
In a factored MDP,
we define a state transition model $\trans$ using a \emph{dynamic
Bayesian network (DBN)}~\cite{Dean+Kanazawa:1989}.  Let $X_i$ denote the
variable $X_i$ at the current time and $X_i'$, the same variable at the
next step.  The \emph{transition graph\/} of a DBN is a two-layer
directed acyclic graph $\graph_\trans$ whose nodes are
$\{X_1,\ldots,X_n,X_1',\ldots,X_n'\}$.  We denote the parents of
$X_i'$ in the graph by $\Parents_\trans(X_i')$.  For simplicity of
exposition, we assume that $\Parents_\trans(X_i') \subseteq \bX$;
thus, all arcs in the DBN are between variables in consecutive
%RP1: grammar, and didactic has a moralistic connotation
%time slices. (This assumption is used for didactic purposes only,
time slices. (This assumption is used for expository purposes only;
intra-time slice arcs are handled by a small modification
presented in \secref{sec: one-step lookahead}.) Each node $X_i'$
is associated with a \emph{conditional probability distribution
(CPD)\/} $P_\trans(X_i' \mid \Parents_\trans(X_i'))$.  The
transition probability $P_\trans(\bx' \mid \bx)$ is then defined
to be:
%CGjairf2 made this an equation
\bdm
P_\trans(\bx' \mid \bx)= \prod_i P_\trans(x_i' \mid \bu_i) \ , 
\edm 
where $\bu_i$ is the
value in $\bx$ of the variables in $\Parents_\trans(X_i')$.



\begin{figure}[t]
\begin{center}
\begin{tabular}{ccc}
\multicolumn{2}{c}{
\hspace{-1cm}
\begin{tabular}{cc}
\epsfig{figure=Figures/ring-example.eps,height=4cm} &
\epsfig{figure=Figures/ring-example-dbn.eps,height=5cm}
\end{tabular}}
&
\begin{footnotesize}
\begin{tabular}{|l|c|c|}
\multicolumn{3}{l}{$P(X_i' = \ttrue\mid X_i, X_{i-1},A)$:} \vspace{.2cm}\\
 \hline
 & \multicolumn{2}{|c|}{Action is reboot:}\\
& machine $i$ &  other machine\\
 \hline
\begin{tabular}{c}
 $X_{i-1} = \ffalse \ \wedge$ \\ $X_i = \ffalse$
 \end{tabular}& 1 & 0.0238\\
 \hline
\begin{tabular}{c}
 $X_{i-1} = \ffalse \ \wedge$ \\ $X_i = \ttrue$
 \end{tabular} & 1 & 0.475\\
 \hline
\begin{tabular}{c}
 $X_{i-1} = \ttrue \ \wedge$ \\ $X_i = \ffalse$
 \end{tabular} & 1 & 0.0475\\
 \hline
\begin{tabular}{c}
 $X_{i-1} = \ttrue \ \wedge$ \\ $X_i = \ttrue$
 \end{tabular} & 1 & 0.95\\
 \hline
\end{tabular}
\end{footnotesize}
\\
\tab\hspace{.5cm}{\bf (a)} & \tab\tab{\bf (b)} & \tab{\bf (c)}
\end{tabular}
\caption{Factored MDP example: from a network topology (a) we obtain the factored MDP representation (b) with the CPDs described in (c).
}
\label{fig:ring example}
\end{center}
\end{figure}

\begin{example}
%RPfinal: removed redundant ``for example''
Consider an instance of the {\it SysAdmin\/} problem with
four computers, labelled $M_1,\dots,M_4$, in an unidirectional ring topology
as shown in \figref{fig:ring example}(a).  Our first task in
modeling this problem as a factored MDP is to define the state
space $\bX$.  Each machine is associated with a binary random
variable $X_i$, representing whether it is working or has failed.
Thus, our state space is represented by four random variables: $\{
X_1, X_2, X_3, X_4\}$. The next task is to define the transition
model, represented as a DBN. The parents of the next time step variables
$X_i'$ depend on the network topology. Specifically, the
probability that machine $i$ will fail at the next time step
depends on whether it is working at the current time step and on
the status of its direct neighbors 
%CGjairf2 made this a () 
(parents in the topology) in the network at the current time step. As shown
in \figref{fig:ring example}(b), the parents of $X_i'$ in this
example are $X_i$ and $X_{i-1}$. The CPD of $X_i'$ is such that
if $X_i = \false$, then $X_i' = \false$ with high probability;
that is, failures tend to persist.  If $X_i = \true$, then $X_i'$ is
a \emph{noisy or} of its other parents (in the unidirectional ring
topology $X_i'$ has only one other parent $X_{i-1}$); that is, a
failure in any of its neighbors can independently cause machine
$i$ to fail. \qed
\end{example}

%CGjairf2 changed from here to the end of this subsection to remove the default action model.
%RP4:  Reworded slightly
%We have described how to represent a single transition model as a DBN in
%factored MDPs.
We have described how to represent factored the Markovian transition
dynamics arising from an MDP as a DBN, but we have not directly
addressed the representation of actions.
Generally, we can define the transition dynamics of
an MDP by defining a separate DBN model $\trans_a =
\langle\graph_a,P_a\rangle$ for each action $a$. 

\begin{example}
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. The transition model for $a_i$ is different from $d$ only
in the transition model for the variable $X_i'$, which is now
$X_i'=\true$ with probability one, regardless of the status of the
neighboring machines.  
%CGff1 fixed sentence
\figref{fig:ring example}(c) shows 
the actual CPD for $P(X_i' =
{Working}\mid X_i, X_{i-1},A)$, with one entry for each
assignment to the state variables $X_i$ and $X_{i-1}$, and to the
action $A$. 
\qed
\end{example}

To fully specify an MDP, we also 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. In our example, we might
have a reward function associated with each machine $i$, which
depends on $X_i$. That is, the SysAdmin is paid on a per-machine
basis: at every time step, she receives money for machine $i$ only
if it is working. We can formalize this concept of localized
functions:
\begin{defn}\label{dfn:restricted}
%RP3: scope
%DK3 - changed to emphasize scope
A function $f$ has a \emph{scope} $\Scope[f] = \Cluster
\subseteq \bX$ if $f : \Val(\Cluster) \mapsto \emph{$\Real$}$.
\qed
\end{defn}
If $f$ has scope $\bY$ and $\bY \subset \bZ$,
we use $f(\bz)$ as shorthand for $f(\by)$ where $\by$ is the
part of the instantiation $\bz$ that corresponds to variables in
$\bY$. 

We can now characterize the concept of local rewards. Let
$R_1^a,\ldots,R_r^a$ be a set of functions, where the scope of each $R_i^a$ 
is restricted to variable cluster $\rclust_i^a \subset
%RP3: Removed ``function'' as per review complaint
\{X_1,\ldots,X_n\}$. The reward for taking action $a$ at state $\bx$ is
defined to be $R^a(\bx) = \sum_{i=1}^r R_i^a(\rclust_i^a) \in \Real$. In our
example, we have a reward function $R_i$ associated with each
machine $i$, which depends only $X_i$, and does not depend on the action choice. These local rewards are
represented by the diamonds in \figref{fig:ring example}(b), in
the usual notation for influence
diagrams~\cite{Howard+Matheson:1981}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Approximate Solution Algorithms}
\label{sec:approximate MDP algs}
There are several algorithms to compute the optimal policy in an MDP.
The three most commonly used are value iteration, policy iteration, and
linear programming.  A key component in all three algorithms is the
computation of value functions, as defined in \secref{sec: MDP}.  Recall
that a value function defines a value for each state $\bx$ in the state
%RPfinal:  pluralization
space.  With an explicit representation of the value function as a vector of
values for the different states, the solution algorithms all can be
implemented as a series of simple algebraic steps.  Thus, in this case,
all three can be implemented very efficiently.

Unfortunately, in the case of factored MDPs, the state space is
exponential in the number of variables in the domain.  
In the {\it SysAdmin\/} problem, for example, the state $\bx$ of the system is an 
assignment describing which machines are working or have failed;
that is, a state $\bx$ is an assignment to each random variable
$X_i$. Thus, the number of states is exponential in the number $m$
of machines in the network ($\abs{\bX} = N = 2^m$). Hence, even
representing an explicit value function in problems with more than about
ten machines is infeasible. 
%In general, real-world domains tend
%to have immense state spaces and we need to perform our
%Solution algorithms rely on the ability to represent value
%functions and policies, these representations require 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 trivial factored MDPs,
%RP3: Softened for reviewer #2
% the
%value function rarely has any internal
%structure~\cite{Koller+Parr:1999}.
there is no guarantee that structure in the model is
preserved in the value function~\cite{Koller+Parr:1999}.

In this section, we discuss the use of an \emph{approximate} value
function, that admits a compact representation.  We also describe
approximate versions of these exact algorithms, that use approximate
value functions.  Our description in this section is 
somewhat abstract, and does not specify how the basic operations
required by the algorithms can be performed explicitly.  In later
sections, we elaborate on these issues, and describe the algorithms in
detail.  For brevity, we choose to focus on policy iteration and linear
programming; our techniques easily extend to value iteration.


%RP3:  Fixed as per reviewer #1's comments

\commentout{
\paragraph{Value Iteration:} Iterate over value functions, Relying
on the fact that the Bellman operator
is a \emph{contraction} --- it is guaranteed to reduce the \emph{max-norm}
($\Linf$)\footnote{For any vector $\V$, the max-norm is given by:
$\maxnorm{\V} = \max_{\bx} \abs{\V(\bx)}$.} distance between any pair of
value functions by a factor of at least $\discount$.\footnote{More
formally, for any two value functions $\V_1$ and $\V_2$, the contraction
property of the Bellman operator $\Tbell$ implies that: $\maxnorm{\Tbell
\V_1 - \Tbell \V_2} \leq \discount \maxnorm{\V_1 - \V_2}$.}
The contraction property guarantees that the Bellman operator has a unique
fixed point $\Vstar$~\cite{Puterman:1994}.  Value iteration exploits this
property, starting with an initial estimate $\V\at{0}$, it approaches the
fixed point through successive applications of the Bellman operator:
$\V\at{t+1} = \Tbell \V\at{t}$. After a finite number of iterations, the
greedy policy $\greedy{\V\at{t}}$ will be the optimal policy.
}

\subsection{Linear Value Functions}

A very popular choice for approximating value functions is by using 
\emph{linear regression}, as first proposed by Bellman \etal~\citeyear{Bellman+al:1963}. 
Here, we define our space of allowable
value functions $\V \in \Basisspace \subseteq \Real^N$ via a set of {\em
basis functions\/}:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Linear value function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
A \emph{linear value function} over a set of basis functions
$\Basis = \{\basis_1,\ldots,\basis_k\}$ is a function $\V$ that
can be written as $\V(\bx) = \sum_{j=1}^k w_j \ \basis_j(\bx)$ for
some coefficients $\bw = (w_1,\ldots,w_k)'$. \qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We can now define $\Basisspace$ to be the linear subspace of
$\Real^N$ spanned by the basis functions $\Basis$. It is useful to
define an $N \times k$ matrix $\Basismatrix$ whose columns are the $k$ basis
%RP3: removed comma as per reviewer 3
functions viewed as vectors. In a more compact notation, our
approximate value function is then represented by $\Basismatrix \bw$.



%CGff1 change paragraph
The expressive power of this linear representation is equivalent,
for example, to that of a single layer neural network with
features corresponding to the basis functions defining
$\Basisspace$. Once the features are defined, we must optimize the
coefficients $\bw$ in order to obtain a good approximation for the
true value function. We can view this approach as separating the
problem of defining a reasonable space of features and the induced
space $\Basisspace$, from the problem of searching within the
space. The former problem is typically the purview of domain
experts, while the latter is the focus of analysis and algorithmic
design. Clearly, feature selection is an important issue for
essentially all areas of learning and approximation.  We offer some simple methods
for selecting good features for MDPs in \secref{sec:Conclusion}, but it
is not our goal to address this large and important topic in this paper.

%RP3: reworded for reviewer
%Once we have defined a representation for the approximate value function,
%i.e., a linear combination of basis functions, the problem becomes
Once we have a chosen a linear value function representation and a set of
basis functions, the problem becomes
one of finding values for the weights $\bw$ such that $\Basismatrix \bw$ will yield a
good approximation of the true value function. In this paper, we consider
two such approaches: approximate dynamic programming using policy
iteration and approximate linear programming. In this section, we present these two approaches.  
%CGjairf2 fixed reference
In \secref{sec: basic ops}, we show
how we can exploit problem structure to transform these approaches into
practical algorithms that can deal with exponentially large state spaces.

\subsection{Policy Iteration} 
\label{sec: Approximate Policy Iteration}
\subsubsection{The Exact Algorithm}
%RP4: reworded slightly
The exact policy iteration algorithm iterates over policies, producing
an improved policy at each iteration.  Starting with some
initial policy $\policy\at{0}$, each iteration consists of two phases.
\emph{Value determination} computes, for a policy $\policy\at{t}$, the
value function $\V_{\policy\at{t}}$, by finding the fixed point of the
equation $\T_{\policy\at{t}} \V_{\policy\at{t}} = \V_{\policy\at{t}}$,
that is, the unique solution to the set of linear equations:
\bdm\label{eqn: value determination} \V_{\policy\at{t}}(\bx) =
R(\bx,\policy\at{t}(\bx)) + \discount \sum_{\bx'} P(\bx' \mid
\bx,\policy\at{t}(\bx)) \V_{\policy\at{t}}(\bx'), \forall \bx.  \edm
%CGff1 fixed this paragraph
The \emph{policy improvement} step defines the next policy as
\[
\policy\at{t+1} = \greedy{\V_{\policy\at{t}}}.
\]
It can be shown that this process converges to the optimal policy~\cite{Bertsekas+Tsitsiklis:1996}.  Furthermore, in
practice, the convergence to the optimal policy is often very quick.
%DK3
\commentout{in at
most as many iterations as value iteration~\cite{Puterman:1994}. In 
practice, policy iteration tends to find the optimal policy in many fewer
%RP3: added for reviewer
iterations --- though each iteration is more costly computationally.
}

\subsubsection{Approximate Policy Iteration}
The steps in the policy iteration algorithm require a manipulation of
%RP4:  added often for reviewer sensitivies
both value functions and policies, both of which often cannot be represented
explicitly in large MDPs.  To define a version of the policy iteration
algorithm that uses approximate value functions, we use the following
basic idea:  We restrict the algorithm to using only value functions
within the provided  $\Basisspace$; whenever the algorithm takes a step
that results in a value function $\V$ that is outside this space, we
\emph{project} the result back into the space by finding the value 
%RP3: close -> closest
function within the space which is closest to $\V$.  More precisely:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Projection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
A \emph{projection operator} $\Pi$ is a mapping $\Pi:\Real^N
\rightarrow \Basisspace$.  $\Pi$ is said to be a \emph{projection w.r.t.\ a
norm $\norm{\cdot}$} if \ \ $\Pi \V = \Basismatrix\bwstar \mbox{~such that~}  \bwstar \in
\arg\min_\bw \norm{\Basismatrix \bw - \V}$. \qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
That is, $\Pi \V$ is the linear combination of the basis
functions, that is closest to $\V$ with respect to the chosen
norm.

Our approximate policy iteration algorithm performs the policy
improvement step exactly. In the value determination step, the value
function --- the value of acting according to the current policy
$\policy\at{t}$ --- is approximated through a linear combination of basis 
functions.  

We now consider the problem of value determination for a policy $\policy\at{t}$.
At this point, it is useful to introduce some notation:
Although the rewards are a function of the state and action choice, once
the policy is fixed, the rewards become a function of the state only,
which we denote as $R_{\policy\at{t}}$, where $R_{\policy\at{t}}(\bx) =
R(\bx, \policy\at{t}(\bx))$. Similarly, for the transition model: 
$P_{\policy\at{t}}(\bx' \mid \bx) = P(\bx' \mid \bx,\policy\at{t}(\bx))$.
We can now rewrite the value determination step in terms of matrices and
vectors. If we view $\V_{\policy\at{t}}$ and $R_{\policy\at{t}}$ as
$N$-vectors, and $P_{\policy\at{t}}$ as an $N \times N$ matrix, we have
the equations: 
\bdm
\V_{\policy\at{t}} =
R_{\policy\at{t}} + \discount P_{\policy\at{t}}
\V_{\policy\at{t}}.
\edm
This is a system of linear equations with one equation for each
state, which can only be solved exactly for relatively small $N$.
Our goal is to provide an approximate solution, within $\Basisspace$.
More precisely, we want to find: 
\beqa
  \bw\at{t} & = & \arg\min_\bw \norm{\Basismatrix\bw - \left(R_{\policy\at{t}} + \discount
P_{\policy\at{t}}\Basismatrix\bw\right)};\\
& = & \arg\min_\bw
\norm{\left(\Basismatrix - \discount P_{\policy\at{t}}\Basismatrix\right)\bw\at{t} - R_{\policy\at{t}}}.
\eeqa
Thus, our approximate policy iteration alternates between two steps:
\begin{eqnarray}
\bw\at{t} & = & \arg\min_\bw \norm{\Basismatrix\bw - \left(R_{\policy\at{t}} + \discount
P_{\policy\at{t}}\Basismatrix\bw\right)}; \label{eq:policy iteration projection}\\
  \policy\at{t+1} & = & \greedy{\Basismatrix \bw\at{t}}.
\end{eqnarray}

\subsubsection{Max-norm Projection}
%RP1
%Unfortunately, algorithms previous to the one we proposed
%Unfortunately, algorithms previous to the one we proposed
%RP1
%in~\cite{Guestrin+al:2001b} suffer from a problem that we might
%RPfinal:  ``these'' was ambiguous
An approach along the lines described above
has been used in various papers, with
several recent theoretical and algorithmic results
\cite{Schweitzer+Seidmann:1985,Tsitsiklis+VanRoy:1996a,%
VanRoy:Thesis,Koller+Parr:1999,Koller+Parr:2000}.
%RP3: For reviewer 3
However, these approaches suffer from a problem that we
might call ``norm incompatibility.''  When computing the projection,
%they all utilize the standard projection operator with respect to
they utilize the standard Euclidean projection operator with respect to
the $\cL_2$ norm or a \emph{weighted} $\cL_2$ norm.%
\footnote{Weighted $\cL_2$ norm projections are stable and have
meaningful error bounds when the weights correspond to the stationary
distribution of a fixed policy under evaluation (value determination)
\cite{VanRoy:Thesis}, but they are not stable when combined with
$\Tbell$.  Averagers \cite{Gordon:1995} are stable and non-expansive in
$\Linf$, but require that the mixture weights be determined \emph{a
priori}.  Thus, they do not, in general, minimize $\Linf$ error.} 
On the other hand, most of the convergence and error 
analyses for MDP algorithms utilize max-norm ($\Linf$). This
incompatibility has made it difficult to provide error guarantees.

%RP3: Toned down
%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 $\Linf$ norm.
We can tie the projection operator more closely to the
error bounds through the use of a projection operator in $\Linf$ norm.
The problem of minimizing the $\Linf$ norm has been
studied in the optimization literature as the problem of finding
the Chebyshev solution%
\footnote{The Chebyshev norm is also referred to as max, supremum and
$\Linf$ norms and the minimax solution.} 
to an overdetermined linear system of equations~\cite{Cheney:1982}. The
problem is defined as finding $\bwstar$ such that:
%RP1:  Should the C below be an A?
\beqn\label{eqn:minimax}
  \bwstar \in \arg\min_\bw \maxnorm{C \bw - \vecb}.
\end{equation}

We use an algorithm due to Stiefel~\citeyear{Stiefel:1960}, that solves
this problem by linear programming: 
%RP1:  Again, should these be a_{ij} instead of c_{ij}?
\beqn\label{eqn:max-norm linear-program}
\begin{array}{lll}
%  \mbox{variables: } & w_1,\dots,w_k,\ \phi;
  \mbox{Variables:} & {w_1,\ldots, w_k,\phi} \ ; \\  \mbox{Minimize: } & \phi \ ; \\
  \mbox{Subject to:} &  \phi  \geq  \sum_{j=1}^k c_{ij} w_j - b_i & \mbox{and} \\
                     &  \phi  \geq  b_i - \sum_{j=1}^k c_{ij} w_j,
                                  & i = 1 ... N. \ \
\end{array}
\end{equation}
The constraints in this linear program imply that $\phi
\geq  \abs{\sum_{j=1}^k c_{ij} w_j - b_i}$ for each $i$, or
equivalently, that $\phi \geq  \maxnorm{C \bw - \vecb}$. The
objective of the LP is to minimize $\phi$. Thus, at the solution
$(\bw^*,\phi^*)$ of this linear program, $\bw^*$ is the solution
of \eqnref{eqn:minimax} and $\phi$ is the $\Linf$ projection
error. 

We can use the $\Linf$ projection in the context of the approximate
policy iteration in the obvious way.  When implementing the projection
operation of \eqnref{eq:policy iteration projection}, we can use the
$\Linf$ projection (as in \eqnref{eqn:minimax}), where $C = (\Basismatrix -
\discount P_{\policy\at{t}}\Basismatrix)$ and $\vecb = R_{\policy\at{t}}$.  This
minimization can be solved using the linear program of~(\ref{eqn:max-norm
linear-program}). 

A key point is that this LP only has $k+1$ variables. However, there
are $2N$ constraints, which makes it impractical for large
state spaces. In the {\it SysAdmin\/} problem, for example, the number of
constraints in this LP is exponential in the number of
machines in the network (a total of $2 \cdot 2^m$ constraints for $m$
machines).  
%CGjairf2 fixed ref
In \secref{sec: basic ops}, we show that, in
\emph{factored} 
%RP1: Spacing looked weird; added with factored basis functions
%MDPs}, all the $2 \  N$ constraints can be represented
MDPs with linear value functions, all the $2N$ constraints can be
represented efficiently, leading to a tractable algorithm.

\subsubsection{Error Analysis}
%RP3: Toned down
% In the remainder of this section, we will discuss how
%this projection addresses  the norm incompatibility problem.
We motivated our use of the max-norm projection within the approximate
policy iteration algorithm via its compatibility with standard error
analysis techniques for MDP algorithms.  We now provide a careful
analysis of the impact of the $\Linf$ error introduced by the
projection step.  The analysis provides motivation for the use of a
projection step that directly minimizes this quantity.  We acknowledge,
however, that the main impact of this analysis is motivational.  In
practice, we cannot provide \emph{a priori} guarantees that an $\Linf$
projection will outperform other methods.

%We can also use our $\Linf$ projection in the context of policy
%iteration.  As we discussed in \secref{sec: MDP}, policy iteration is
%composed of two steps: value determination and policy improvement. 

Our goal is to analyze approximate policy iteration in terms of the
amount of error introduced at each step by the
projection operation. If the error is zero, then we are performing 
exact value determination, and no error should accrue.  If the error is
small, we should get an approximation that is accurate.  This result
follows from the analysis below.  More precisely, we define the
projection error as the error resulting from the approximate value
determination step: 
\bdm
\beta\at{t} = \maxnorm{\Basismatrix\bw\at{t} - \left(R_{\policy\at{t}} + 
   \discount P_{\policy\at{t}}\Basismatrix\bw\at{t}\right)}. 
\edm
Note that, by using our max-norm projection, we are finding the set of
weights $\bw\at{t}$ that exactly minimizes the one-step projection error
$\beta\at{t}$. That is, we are choosing the best possible weights with
respect to this error measure. Furthermore, this is exactly the error
measure that is going to appear in the bounds of our theorem. Thus, we can
now make the bounds for each step as tight as possible.  

We first show that the projection error accrued in each step is bounded:
%Unlike in approximate value iteration, the one-step projection error
%$\beta\at{t}$ is bounded for all iterations in the approximate policy
%iteration algorithm: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lemma: value determination error is bounded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace*{-.1cm}
\begin{lem}\label{lem: vd error is bounded}
The value determination error is bounded: There exists a constant
$\betap \leq \Rmax$ such that $\betap \geq \beta\at{t}$ for all
iterations $t$ of the algorithm.
\end{lem}
\proof See \appenref{proof:lem: vd error is bounded}. \qed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Due to the contraction property of the Bellman operator, the overall
accumulated error is a decaying average of the projection error incurred
throughout all iterations:  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: discounted value determination error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
The \emph{discounted value determination error} at iteration $t$
is defined as: $\betabar\at{t} = \beta\at{t} + \discount
\betabar\at{t-1}; \ \ \betabar\at{0} = 0.$\qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\lemref{lem: vd error is bounded} implies that the accumulated error
remains bounded in approximate policy iteration: $\betabar\at{t} \leq
\frac{\betap (1-\discount^t)}{1-\discount}$.
%RP4:  collapsed into one paragraph
We can now bound the loss incurred when acting according to the policy
generated by our approximate policy iteration algorithm, as opposed to
the optimal policy: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theorem: pi error bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thm}\label{thm: pi error bound}
In the approximate policy iteration algorithm, let $\policy\at{t}$ be the
policy generated at iteration $t$. Furthermore, let $\V_{\policy\at{t}}$ be
%RP1: grammar
%the \emph{actual} value of acting according to this policy. Then, the loss
the \emph{actual} value of acting according to this policy. The loss
incurred by using policy $\policy\at{t}$ as opposed to the optimal policy
$\policystar$ with value $\Vstar$ is bounded by:
\beqn \label{eq: pi error bound}
\maxnorm{\Vstar - \V_{\policy\at{t}} }  \leq
  \discount^t \maxnorm{\Vstar - \V_{\policy\at{0}}} +
       \frac{2\discount\betabar\at{t}}{(1-\discount)^2}.
\eeqn
\end{thm}
%RP1:  I wanted to change this to de-emphasize the intial details because
%      they are standard, while emphasizing the latter part of the proof,
%      but I decided that it's not worth the trouble.
\proof See \appenref{proof:thm: pi error bound}.
\qed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%RP3:  Added and rearranged to address reviewer confusion

In words, \eqnref{eq: pi error bound} shows that 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 \lemref{lem: vd error is bounded} shows,
is bounded. This second term can be minimized by choosing
$\bw\at{t}$ as the one that minimizes: \bdm \maxnorm{\Basismatrix\bw\at{t} -
\left(R_{\policy\at{t}} + \discount
P_{\policy\at{t}}\Basismatrix\bw\at{t}\right)}, \edm which is exactly the
computation performed by the max-norm projection. Therefore, this theorem
motivates the use of max-norm projections to minimize the error term that
appears in our bound.

The bounds we have provided so far may seem fairly trivial, as we have
not provided a strong \emph{a priori} bound on $\beta\at{t}$.  
Fortunately, several factors make these bounds interesting despite the lack
of \emph{a priori} guarantees.
If approximate policy iteration converges, as occurred in all of
our experiments, we can obtain a much tighter bound: 
%CGff1 added new line here
If
$\policyhat$ is the policy after convergence, then:
\[\maxnorm{\Vstar - \V_{\policyhat} } \leq
\frac{2\discount\beta_{\policyhat}}{(1-\discount)}, 
\]
where
$\beta_{\policyhat}$ is the one-step max-norm projection error
associated with estimating the value of $\policyhat$.  Since the max-norm
projection operation provides 
%CGjairf2 fixed () -> {}
$\beta_{\policyhat}$, we can easily obtain an
\emph{a posteriori} bound as part of the policy iteration procedure.
More details are provided in \secref{sec: Bellman error}.

%RP3: could could
One could rewrite the bound in \thmref{thm: pi error bound}
in terms of the worst case projection error $\betap$, or the worst
projection error in a cycle of policies, if approximate policy
iteration gets stuck in a cycle.  These formulations would be
closer to the analysis of Bertsekas and
Tsitsiklis~\citeyear[Proposition 6.2,
p.276]{Bertsekas+Tsitsiklis:1996}. However, consider the case
where most policies (or most policies in the final cycle) have a
low projection error, but there are a few policies that cannot be
approximated well using the projection operation, so that they have a
large one-step projection error.  A worst-case bound would be very
loose, because it would be dictated by the error of the most
difficult policy to approximate. On the other hand, using our
discounted accumulated error formulation, errors introduced by
policies that are hard to approximate decay very rapidly. Thus,
the error bound represents an ``average'' case analysis: a
decaying average of the projection errors for policies encountered
at the successive iterations of the algorithm.  As in the convergent case,
this bound can be computed easily as part of the policy iteration procedure
when max-norm projection is used. 

The practical benefit of \emph{a posteriori} bounds is that they can give
meaningful feedback on the impact of the choice of the value function
approximation architecture.  While we are not explicitly addressing the
difficult and general problem of feature selection in this paper, our error
bounds motivate algorithms that aim to minimize the error \emph{given} an
approximation architecture and provide feedback that could be useful in
future efforts to automatically discover or improve approximation
architectures.

\commentout{ We can now specify the complete approximate policy
iteration algorithm: \benum
  \item {\bf Approximate value determination:} for a fixed policy
  $\policy\at{t}$, solve:
  \bdm
  \bw\at{t} = \arg\min_\bw\maxnorm{\left(I - \discount P_{\policy\at{t}}\right)\Basismatrix\bw -
R_{\policy\at{t}}}.
  \end{displaymath}
  \item {\bf Update discounted value determination error:} $\betabar\at{t} =
  \beta\at{t} + \discount \betabar\at{t-1}$.
  \item {\bf Stopping criterion: } if $\Bellman(\Basismatrix\bw\at{t}) \leq \eps$ return $\Basismatrix\bw\at{t}$;
  \item {\bf Policy improvement:} the next policy is obtained by:
  \bdm
    \policy\at{t+1}(s) = \max_a[{R(s,a) + \discount \sum_{s'} P(s'\mids,a)\V_{\policy\at{t}}(s')}]; ~~~ \V_{\policy\at{t}} = \Basismatrix\bw\at{t}.
  \end{displaymath}
\eenum

As in the value iteration case, the stopping criterion guarantees:
$\maxnorm{\Vstar - \Basismatrix \bwstar} \leq \frac{2 \eps}{1 - \discount}$.
In this algorithm, we are also able to maintain a the discounted
value determination error $\betabar\at{t}$, the discounted value
approximation error incurred throughout the execution of the
algorithm. As in the approximate value iteration case, it measures
the amount of error that was introduced by successive projection
operations and can be used to evaluate if our basis is well suited
for the problem. 
} %end commentout
%The policy improvement step can be performed in the usual manner,
%as described in \secref{Policy Iteration}. However, this involves
%storing the new policy, $\policy\at{t+1}$, which could be
%intractable for a large state space.  In the second half of the
%paper, we will present an algorithm that addresses this problem.


%----------------------------------------------
%----------------------------------------------
%----------------------------------------------
\subsection{Approximate Linear Programming}
\label{sec: Approximate Linear Programming}

\subsubsection{The Exact Algorithm}
Linear programming provides an alternative method for solving MDPs.  
It formulates the problem of finding a value function as a linear program
(LP).  Here the LP variables are $V_1,\ldots,V_N$, where $V_i$ represents 
$\V(\bx_i)$: the value of starting at the $i$th state of
the system. The LP is given by: \beqn \label{eq: exact LP for
MDPs}
\begin{array}{lll}
  \mbox{Variables:} & {V_1,\ldots, V_N} \ ; \\
  \mbox{Minimize: } & \sum_{\bx_i} \sweights(\bx_i) \ V_i \ ; \\
  \mbox{Subject to:} & V_i \geq
       [R(\bx_i,a) + \discount \sum_{j} P(\bx_j \mid \bx_i,a) V_j] &
\forall \bx_i \in \bX, a \in A,
\end{array}
\eeqn
where the \emph{state relevance weights} $\sweights$ are positive.
%RP1: They don't need to sum to one
%($\sweights(\bx_i) > 0, \forall \bx_i$) and sum to one
%($\sum_{\bx_i}\sweights(\bx_i)=1$). Note that, in this exact case, the
%($\sweights(\bx_i) > 0, \forall \bx_i$). Note that, in this exact case, the
Note that, in this exact case, the
solution obtained is the same for any positive weight vector.
%RP1: Changed
%$\sweights$. It is interesting to note that policy iteration corresponds to
%a large step application of the Simplex algorithm to this linear
%program. Thus, in practice, policy iteration tends to be faster than this
%linear programming approach~\cite{Puterman:1994}.
%$\sweights$.
It is interesting to note that steps of the simplex algorithm
correspond to policy changes at single states, while steps of policy
iteration can involve policy changes at multiple states.
In practice, policy iteration tends to be faster than the
linear programming approach~\cite{Puterman:1994}.



\commentout{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Bellman error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
The \emph{Bellman error} of a value function $\V$ is defined as:
$\Bellman(\V) = \maxnorm{\Tbell \V - \V}.$\qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%}
If $\Bellman(\V) \leq \eps$, we have that $\maxnorm{\Vstar - \V} \leq
\frac{\eps}{1 - \discount}$~\cite{Williams+Baird:1993}. Thus, the Bellman
error can be used to bound the error between an estimated value function
$\V$ and $\Vstar$.
}

\subsubsection{Approximate Linear Program}

The approximate formulation for the LP approach, first proposed by
Schweitzer and Seidmann~\citeyear{Schweitzer+Seidmann:1985}, restricts
the space of allowable value functions to the linear space spanned by our
basis functions.  In this approximate formulation, the variables are
$w_1,\ldots,w_k$: the weights for our basis functions. The LP is
given by: 
\beqn \label{eq: approx LP for MDPs}
\begin{array}{lll}
  \mbox{Variables:} & {w_1,\ldots,w_k} \ ; \\
  \mbox{Minimize: } & \sum_{\bx} \sweights(\bx) \sum_i w_i \  \basis_i(\bx) \ ; \\
  \mbox{Subject to:} & \sum_i w_i \  \basis_i(\bx) \geq
       [R(\bx,a) + \discount \sum_{\bx'} P(\bx' \mid \bx,a) \sum_i w_i \  \basis_i(\bx')] &
\forall \bx \in \bX, \forall a \in A .
\end{array}
\eeqn 

In other words, this formulation takes the LP in (\ref{eq: exact LP for
MDPs}) and substitutes the explicit state value 
function by a linear value function representation $\sum_i w_i \
\basis_i(\bx)$, or, in our more compact notation, $\V$ is replaced
by $\Basismatrix\bw$. This linear program is guaranteed to be feasible if 
a constant function --- a function with the same constant value
for all states --- is included in the set of basis functions.

In this approximate linear programming formulation, the choice of
state relevance weights, $\sweights$, becomes important.
Intuitively, not all constraints in this LP are binding; that is, the 
constraints are tighter for some states than for others. For
each state $\bx$, the relevance weight $\sweights(\bx)$ indicates
%RP3: fixed as per reviewer complaint
the relative importance of a tight constraint. Therefore, unlike
the exact case, the solution obtained may differ for different
choices of the positive weight vector $\sweights$. Furthermore,
there is, in general, no guarantee as to the quality of the greedy
policy generated from the approximation $\Basismatrix\bw$. However, the
recent work of de Farias and Van Roy~\citeyear{Farias+VanRoy:2001}
provides some analysis of the error relative to that of the best
possible approximation in the subspace, and some guidance as to
selecting $\sweights$ so as to improve the quality of the
approximation. In particular, their analysis shows that this LP
provides the best approximation $\Basismatrix\bwstar$ of the optimal value
function $\Vstar$ in a weighted $\cL_1$ sense subject to the
constraint that $\Basismatrix\bwstar \geq \Tbell \Basismatrix\bwstar$, where the weights
in the $\cL_1$ norm are the state relevance weights $\sweights$.

The transformation from an exact to an approximate problem
formulation has the effect of reducing the number of free
variables in the LP to $k$ (one for each basis function
coefficient), but the number of constraints remains $N \times
|A|$.  In our {\it SysAdmin\/} problem, for example, the number of
constraints in the LP in~(\ref{eq: approx LP for MDPs}) is $(m+1)
\cdot 2^m$, where $m$ is the number of machines in the network.
Thus, the process of generating the constraints and solving the LP
still seems unmanageable for more than a few machines. In the next
section, we discuss how we can use the structure of a factored MDP
to provide for a compact representation and an efficient solution
to this LP.

% -------------------------------------------------------------
% -------------------------------------------------------------
% -------------------------------------------------------------

\section{Factored Value Functions}
\label{sec: basic ops}
The linear value function approach, and the algorithms described in
\secref{sec:approximate MDP algs}, apply to any choice of basis
functions.  In the context of factored MDPs, Koller and
Parr~\citeyear{Koller+Parr:1999} suggest a particular type of basis
function, that is particularly compatible with the structure of a
factored MDP.  They suggest that, although the value function is
typically not structured, there are many cases where it might be
``close'' to structured.  That is, it might be well-approximated using a
linear combination of functions each of which refers only to a small
number of variables.  
%RP4: Collapsed in stray sentence below into previous paragraph 
More precisely, we define:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: factored linear value function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
A \emph{factored (linear) value function} is a linear function
over the basis set $\basis_1,\ldots,\basis_k$, where the scope of 
each $\basis_i$ is restricted to some subset of variables
$\Cluster_i$. \qed 
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Value functions of this type have a long history in the area of
multi-attribute utility theory~\cite{Keeney+Raiffa:1976}. In our
example, we might have a basis function $\basis_i$ for each
machine, indicating whether it is working or not.  Each basis
%RP3: scope
function has scope restricted to $X_i$. These are represented as
diamonds in the next time step in \figref{fig:ring example}(b).

Factored value functions provide the key to performing efficient
computations over the exponential-sized state spaces we have in factored
MDPs.
%~\cite{Koller+Parr:1999,Koller+Parr:2000,Guestrin+al:2001b,Guestrin+al:2001c}.
%RP3: key -> main, scope
The main insight is that restricted scope functions (including our basis
functions) allow for certain basic operations to be implemented very
efficiently.  In the remainder of this section, we 
%CGjairf2 added
show
how structure in
factored MDPs can be exploited to perform two crucial operations very
%RP4:  one-step lookahead value -> one step lookahead
efficiently: one-step lookahead (backprojection), and the
representation of exponentially many constraints in the LPs.
Then, we use these basic building blocks to formulate very efficient
approximation algorithms for factored MDPs, each presented in its own
self-contained section: the approximate linear programming for factored
MDPs in \secref{sec:ALP}, and approximate policy iteration with max-norm  
projection in \secref{sec:API}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{One-step Lookahead}\label{sec: one-step lookahead}

A key step in all of our algorithms is the computation of the
one-step lookahead value of some action $a$. This is necessary,
for example, when computing the greedy policy as in \eqnref{eq:
greedy}. Let's consider the computation of a $Q$ function,
$Q_a(\bx)$, which represents the expected value the agent obtains
after taking action $a$ at the current time step and receiving a
long-term value $\V$ thereafter. This $Q$ function can be computed
by:
\begin{equation}\label{eq: exact q}
Q_a(\bx) = R(\bx,a) + \discount \sum_{\bx'} P(\bx'\mid \bx,a) \V(\bx).
\end{equation}
That is, $Q_a(\bx)$ is given by the current reward plus the
discounted expected future value. Using this notation, we can
express the greedy policy as: $\greedy{\V}(\bx) = \max_a
Q_a(\bx)$.

Recall that we are estimating the long-term value of our policy
using a set of basis functions: $\V(\bx) = \sum_i w_i \
\basis_i(\bx)$. Thus, we can rewrite \eqnref{eq: exact q} as:
\begin{equation}\label{eq: approx q}
Q_a(\bx) = R(\bx,a) + \discount \sum_{\bx'} P(\bx'\mid \bx,a) \sum_i w_i \  \basis_i(\bx).
\end{equation}
The size of the state space is exponential, so that
computing the expectation $\sum_{\bx'} P(\bx'\mid \bx,a) \sum_i
w_i \  \basis_i(\bx)$ seems infeasible. Fortunately, as discussed
by Koller and Parr~\citeyear{Koller+Parr:1999}, this expectation
operation, or backprojection, can be performed efficiently if the
transition model and the value function are both factored
appropriately.
%RP3: Added
The linearity of the value function permits a linear decomposition, where
each summand in the expectation can be viewed as an independent value
function and updated in a manner similar to the value iteration procedure
used by Boutilier \etal~\citeyear{Boutilier+al:2000}.
We now recap the construction briefly,
%RP4:  Collapse into previous paragraphy
by first defining:
\bdm \Backbasis^a(\bx) = \sum_{\bx'} P(\bx'\mid
\bx,a) \sum_i w_i \  \basis_i(\bx') = \sum_i w_i \sum_{\bx'}
P(\bx'\mid \bx,a) \basis_i(\bx').
\edm
Thus, we can compute the
expectation of each basis function separately: \bdm
\backbasis^a_i(\bx) = \sum_{\bx'} P(\bx'\mid \bx,a)
%RP3:  weigh -> weight (reviewer)
\basis_i(\bx'), \edm and then weight them by $w_i$ to obtain the
total expectation $\Backbasis^a(\bx) = \sum_i w_i \
\backbasis^a_i(\bx)$. The intermediate function $\backbasis^a_i$
is called the \emph{backprojection} of the basis function
$\basis_i$ through the transition model $P_a$, which we denote by
$\backbasis^a_i = P_a \basis_i$. Note that, in factored MDPs, the
transition model $P_a$ is factored (represented as a DBN)
%RP3: scope
and the basis functions $\basis_i$ have scope restricted to a
%RP3: key -> important
small set of variables. These two important properties allow us to
compute the backprojections very efficiently.

We now show how some restricted scope function $\basis$
(such as our basis functions) can be backprojected through some
transition model $P_{\trans}$ represented as a DBN $\trans$. Here
$\basis$ has scope restricted to $\bY$; our goal is to compute
$\backbasis = P_\trans \basis$.  We define the
%RP3: scope
\emph{backprojected scope of $\bY$ through $\trans$} as the set of
parents of $\bY'$ in the transition graph $\graph_\trans$;
$\backproj{\trans}{\bY'} = \cup_{Y_i' \in \bY'}
\Parents_{\trans}(Y_i')$. If intra-time slice arcs are included,
so that $\Parents_{\trans}(X_i') \in \{X_1, \dots,
X_n,X_1',\dots,X_n'\}$, then the only change in our algorithm is
%RP3: scope
in the definition of backprojected scope of $\bY$ through $\trans$. The
definition now includes not only direct parents of $Y'$, but
also all variables in $\{X_1, \dots, X_n\}$ that are ancestors of
$Y'$:
\[
  \backproj{\trans}{\bY'} = \{X_j \mid \mbox{there
exist a directed path from $X_j$ to any $X_i' \in \bY'$} \}.
\]
%RP3: scope
Thus, the backprojected scope may become larger, but the functions 
are still factored.

%RP3: scope
We can now show that, if $\basis$ has scope restricted to $\bY$,
%RP3: scope
then its backprojection $\backbasis$ has scope restricted
to the parents of $\bY'$, \ie, $\backproj{\trans}{\bY'}$.
Furthermore, each backprojection can be computed by only
enumerating settings of variables in $\backproj{\trans}{\bY'}$,
rather than settings of all variables $\bX$: \beqa
\backbasis(\bx) & = & (P_\trans \basis)(\bx); \\
     & = & \sum_{\bx'} P_\trans(\bx' \mid \bx) \basis(\bx'); \\
     & = & \sum_{\bx'} P_\trans(\bx' \mid \bx) \basis(\by'); \\
     & = & \sum_{\by'} P_\trans(\by' \mid \bx) \basis(\by')
                 \sum_{\bu' \in (\bx' - \by')} P_\trans(\bu' \mid \bx);\\
     & = & \sum_{\by'} P_\trans(\by' \mid \vbp) \basis(\by');\\
     & = & \backbasis(\vbp);
\eeqa 
%CGjairf2 added lots of ' in \backproj{\trans}{\bY} throughout the paper
where $\vbp$ is the value of $\backproj{\trans}{\bY'}$ in
$\bx$ and the term $\sum_{\bu' \in (\bx' - \by')} P_\trans(\bu'
\mid \bx) = 1$ as it is the sum of a probability distribution over
a complete domain. Therefore, we see that $(P_\trans \basis)$ is a
%RP3:  scope
function whose scope is restricted to $\backproj{\trans}{\bY'}$.
Note that the cost of the computation depends linearly on
$|\Val(\backproj{\trans}{\bY'})|$, which depends on $\bY$ (the
%RP3: scope
scope of $\basis$) and on the complexity of the process dynamics.
%CGjairf1 added
This backprojection procedure is summarized in \figref{code:bproj}.

Returning to our example, consider a basis function $\basis_i$
that is an indicator of variable $X_i$: it takes value $1$ if
the $i$th machine is working and $0$ otherwise. Each $\basis_i$
%RP3: scope
has scope restricted to $X_i'$, thus, its backprojection
%RP3: scope
$\backbasis_i$ has scope restricted to $\Parents_\trans(X_i')$:
$\backproj{\trans}{X_i'} = \{X_{i-1},X_i\}$.


%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Backprojection code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{pseudocode}
$\bprojfun{a}{\basis}$ \ --- \ \ where basis function $\basis$ has scope $\Cluster$.
\bcsone
\item {\bf Define} the scope of the backprojection:  $\backproj{a}{\Cluster'} = \cup_{X_i' \in \Cluster'} \Parents_{a}(X_i')$.
\item {\bf For} each assignment $\by \in \backproj{a}{\Cluster'}$:
\bcstwo
\item $\backbasis^a(\by) = \sum_{\clustval' \in \Cluster'} \prod_{i \mid X_i' \in \Cluster'}  P_{a}(\instin{X_i'}{\clustval'}\mid \by) \basis(\clustval') .$\\
\ecstwo
\item {\bf Return} $\backbasis^a$.
\ecsone
\end{pseudocode}
\caption{Backprojection of basis function $\basis$. }\label{code:bproj}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Representing Exponentially Many Constraints}
\label{sec: Representing exponentially many constraints}

As seen in \secref{sec:approximate MDP algs}, both our approximation
algorithms require the solution of linear programs:  the LP in
(\ref{eqn:max-norm linear-program}) for 
%CGjairf2 removed
%approximate value iteration and 
approximate policy iteration, and the LP in
(\ref{eq: approx LP for MDPs}) for the approximate linear
programming algorithm. These LPs have some common characteristics:
%CGjairf2 rephrased and ()
they have a small number of free variables (for $k$ basis
functions there are $k+1$ free variables in approximate policy iteration and $k$ in approximate linear programming), but the number of constraints is still exponential in the
number of state variables. However, in factored MDPs, these LP
constraints have another very useful property: the functionals in
%RP3: scope
the constraints have restricted scope.
This key observation allows us to
represent these constraints very compactly.
%~\cite{Guestrin+al:2001b}.

First, observe that the constraints in the linear programs are all
of the form:
%%%%%%%
\beqn \label{eq: general constraint form} \phi \geq \sum_i w_i \
c_i(\bx) - b(\bx), \forall \bx, \eeqn
%%%%
where only $\phi$ and $w_1, \dots, w_k$ are free variables in the
LP and $\bx$ ranges over all states. This general form represents
both the type of constraint in the max-norm projection LP in
(\ref{eqn:max-norm linear-program}) and the approximate linear
programming formulation in (\ref{eq: approx LP for
MDPs}).\footnote{The complementary constraints in
(\ref{eqn:max-norm linear-program}), $\phi \geq b(\bx) - \sum_i
w_i \  c_i(\bx)$, can be formulated using an analogous
construction to the one we present in this section by changing the
sign of $c_i(\bx)$ and $b(\bx)$. The approximate linear
programming constraints of (\ref{eq: approx LP for MDPs}) can also
be formulated in this form, as we show in \secref{sec:ALP}.}

%RP3: key -> null
The first insight 
%CGjairf2 added 
in our construction
is that we can replace the entire set of
constraints in \eqnref{eq: general constraint form} by one
equivalent non-linear constraint: 
\beqn \label{eq: max constraint form} 
\phi \geq \max_{\bx}\sum_i w_i \  c_i(\bx) - b(\bx). 
\eeqn
%RP3: key -> null
The second insight is that this new non-linear constraint can be
implemented by a set of linear constraints using a construction
that follows the structure of variable elimination in cost
networks. This insight allows us to exploit structure in
factored MDPs to represent this constraint compactly.

We tackle the problem of representing the constraint in
\eqnref{eq: max constraint form} in two steps: first, computing
the maximum assignment for a fixed set of weights; then,
representing the non-linear constraint by small set of linear
%RP1:
% constraints, using a construct we call \emph{factored LP}.
constraints, using a construction we call the \emph{factored LP}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Maximizing Over the State Space}  \label{sec: Maximizing over the state space}

%CGjairf1 removing overline w's and fixing some notation
The key computation in our algorithms is to represent a non-linear
constraint of the form in \eqnref{eq: max constraint form}
efficiently by a small set of linear constraints. Before
presenting this construction, let's first consider a simpler
problem: Given some \emph{fixed} weights $w_i$, we
would like to compute the maximization: $\phistar =
\max_{\bx}\sum_i w_i \  c_i(\bx) - b(\bx)$, that is,
the state $\bx$, such that the difference between $\sum_i
w_i \ c_i(\bx)$ and $b(\bx)$ is maximal. However, we
cannot explicitly enumerate the exponential number of states and
compute the difference. Fortunately, structure in factored MDPs
allows us to compute this maximum efficiently.


In the case of factored MDPs, our state space is a set of vectors
$\bx$ which are assignments to the state variables $\bX =
\{X_1,\ldots,X_n\}$.  We can view both $C\bw$ and $\vecb$ as
functions of these state variables, and hence also their
difference.  Thus, we can define a function
$F^\bw(X_1,\ldots,X_n)$ such that $F^\bw(\bx) = \sum_i w_i \
c_i(\bx) - b(\bx)$.  Note that we have executed a representation
shift; we are viewing $F^\bw$ as a function of the variables
$\bX$, which is parameterized by $\bw$. Recall that the size of
the state space is exponential in the number of variables. Hence,
our goal in this section is to compute $\max_{\bx} F^\bw(\bx)$
without explicitly considering each of the exponentially many
%RP3:  key -> solution
states.  The solution is to use the fact that $F^\bw$ has a factored
representation. More precisely, $C\bw$ has the form $\sum_i w_i \
c_i(\bZ_i)$, where $\bZ_i$ is a subset  of $\bX$.  For example, we
might have $c_1(X_1,X_2)$ which takes value $1$ in states where
$X_1=\true$ and $X_2=\false$ and $0$ otherwise. Similarly, the
%RP3: scope
vector $\vecb$ in our case is also a sum of restricted scope
functions.  Thus, we can express $F^\bw$ as a sum $\sum_j
f^\bw_j(\bZ_j)$, where $f^\bw_j$ may or may not depend on $\bw$.
In the future, we sometimes drop the superscript $\bw$ when it is
clear from context.

Using our more compact notation, our goal here is simply to
compute $\max_{\bx} \sum_i w_i \ c_i(\bx) - b(\bx) =
\max_{\bx} F^{\bw}(\bx)$, that is, to find the state $\bx$
over which $F^{\bw}$ is maximized.  Recall that $F^{\bw} =
\sum_{j=1}^m f_j(\bZ_j)$.  We can maximize such a function, $F^{\bw}$,
without enumerating every state using \emph{non-serial dynamic
programming}~\cite{Bertele+Brioschi:1972}. The idea is virtually
identical to variable elimination in a Bayesian network. We review
%RP3: key -> central
this construction here, as it is a central component in our solution
LP.

%CGjairf2 fixed a lot of bZ -> \bZ
Our goal is to compute
\[
  \max_{x_1,\ldots,x_n} \sum_j f_j(\instin{\bZ_j}{\bx}).
\]
The main idea is that, 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.  

\begin{example}
Assume
\bdm
  F  =  f_1(x_1,x_2) + f_2(x_1,x_3) +
                         f_3(x_2,x_4) + f_4(x_3,x_4).
\edm
%\beqa
%  F(X_1,X_2,X_3,X_4) & = & f_1(X_1,X_2) + f_2(X_1,X_3) + \\
%                     &   & f_3(X_2,X_4) + f_4(X_3,X_4).
%\eeqa
We therefore wish to compute:
\[
  \max_{x_1,x_2,x_3,x_4} f_1(x_1,x_2) + f_2(x_1,x_3) +
                         f_3(x_2,x_4) + f_4(x_3,x_4).
\]
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
\[
  \max_{x_1,x_2,x_3} f_1(x_1,x_2) + f_2(x_1,x_3) +
       \max_{x_4} [f_3(x_2,x_4) + f_4(x_3,x_4)].
\]
The result of the internal maximization depends on the values of $x_2,x_3$;
thus, 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
\[
  \max_{x_1,x_2,x_3} f_1(x_1,x_2) + f_2(x_1,x_3) + e_1(x_2,x_3),
\]
having one fewer variable.  Next, we eliminate another variable, say
$X_3$, with the resulting expression reducing to:
\[
  \max_{x_1,x_2} f_1(x_1,x_2) + e_2(x_1,x_2) ,
\]
\[
\mbox{where \ \ \ }
 e_2(x_1,x_2) = \max_{x_3} [f_2(x_1,x_3) + e_1(x_2,x_3)].
\]
Finally, we define
\[
  e_3 = \max_{x_1,x_2} f_1(x_1,x_2) + e_2(x_1,x_2).
\]
The result at this point is a number, which is the desired maximum
over $x_1,\ldots,x_4$. While the naive approach of enumerating all
states requires $63$ arithmetic operations if all variables
are binary, using variable elimination we only need to perform
$23$ operations. \qed
\end{example}


%CGjairf1 removed old code description and added new one 
% and rearranged rest of section
\commentout{ In general, the variable elimination algorithm maintains a set
$\F$ of functions,  initially $\{f_1,\ldots,f_m\}$. The algorithm
then repeats the following steps: \benum
\item
Select an uneliminated variable $X_l$;
item
%RP3: take -> let, scope
Let $e_1,\dots ,e_L$ be the functions in  $\F$ whose scope
contains $X_l$.
%RP3:  Added range to summation
\item
Define a new function $e = \max_{x_l} \sum_{j=1}^{L} e_j$
and introduce it into $\F$.
%RP3: scope
The scope of $e$ is $\cup_{j=1}^L \Scope[e_j] - \{X_l\}$.
\eenum
}

The general variable elimination algorithm is described in
\figref{code:VE}. The inputs to the algorithm are the functions to be
maximized $\F = \{f_1,\ldots,f_m\}$ and an elimination ordering $\order$
on the variables, where $\order(i)$ returns the $i$th variable to be
eliminated. As in the example above, for each variable $X_l$ to be
eliminated, we select the relevant functions $e_1,\dots ,e_L$, those
whose scope contains $X_l$.  These functions are removed from the set
$\F$ and we introduce a new function $e = \max_{x_l} \sum_{j=1}^{L}
e_j$. At this point, the scope of the functions in $\F$ no longer depends
on $X_l$, that is, $X_l$ has been `eliminated'. This procedure is
repeated until all variables have been eliminated. The remaining
functions in $\F$ thus have empty scope. The desired maximum
is therefore given by the sum of these remaining functions.  


%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variable elimination code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{pseudocode}
VariableElimination ($\F$, $\order$) \\ 
\tab \codecomment{$\F = \{f_1,\ldots,f_m\}$ is the set of functions to be maximized;}\\
\tab \codecomment{$\order$ stores the elimination order. }
\bcsone
\item {\bf For} $i = 1$ to number of variables:
\bcstwo
\item \codecomment{Select the next variable to be eliminated.}
\item {\bf Let} $l = \order(i)$ ; 
\item \codecomment{Select the relevant functions.}
\item
{\bf Let} $e_1,\dots ,e_L$ be the functions in  $\F$ whose scope
contains $X_l$.
\item \codecomment{Maximize over current variable $X_l$.}
\item
{\bf Define} a new function $e = \max_{x_l} \sum_{j=1}^{L} e_j$ ;
note that $\Scope[e]=\cup_{j=1}^L \Scope[e_j] - \{X_l\}$.
\item \codecomment{Update set of functions.}
\item {\bf Update} the set of functions  $\F = \F \cup \{e\} \setminus \{e_1,\dots ,e_L\}$.
\ecstwo
\item \codecomment{Now, all functions have empty scope and their sum is the maximum value of $f_1+\cdots+f_m$.}
\item {\bf Return} the maximum value $\sum_{e_i \in \F} e_i$.
\ecsone
\end{pseudocode}
\caption{Variable elimination procedure for computing the maximum value $f_1+\cdots+f_m$, where each $f_i$ is a restricted scope function. }\label{code:VE}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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$
%RP3: scope
whose scope is $\bZ$.  To compute this function, we need to
compute $\abs{\Val[\bZ]}$ different values.  The cost of the
algorithm is linear in the overall number of these values,
introduced throughout the execution. As shown
by Dechter~\citeyear{Dechter:1999}, this cost is exponential in the induced
width of the \emph{cost network}, the undirected graph
defined over the variables $X_1,\ldots,X_n$, with an edge between
$X_l$ and $X_m$ if they appear together in one of the original
functions $f_j$.
%RP3:  Changed to address reviewer complaints
%The complexity of this algorithm is, of course,
%dependent on the variable elimination order. Computing the optimal
%elimination order is an NP-hard problem~\cite{Arnborg+al:1987}.
%However, very good heuristics have been developed for this
%problem~\cite{Bertele+Brioschi:1972}.
The complexity of this algorithm is, of course, dependent on the
variable elimination order and the problem structure.  Computing the
optimal elimination order is an NP-hard problem~\cite{Arnborg+al:1987}
and elimination orders yielding low induced tree width do not exist
for some problems.  These issues have been confronted successfully for
a large variety of practical problems in the Bayesian network
community, which has benefited from a large variety of good heuristics
which have been developed for the variable elimination ordering
problem~\cite{Bertele+Brioschi:1972,Kjaerulff:1990,Reed:1992,Becker+Geiger:2001}.
%CGjairf2 added
%RP4:  Remove the e.g, b/c we were using the references as nouns
%(\eg, see~\cite{Bertele+Brioschi:1972,Kjaerulff:1990,Reed:1992,Becker+Geiger:2001}).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Factored LP}\label{sec: Factored LP}


%CGjairf1 umph in the beginning and removed self-citations
In this section, we present the centerpiece of our planning algorithms: a
new, general approach for compactly representing exponentially large sets
of LP constraints in problems with factored structure --- those where the
functions in the constraints can be decomposed as the sum of restricted
scope functions.  
Consider our original problem of representing the non-linear
constraint in \eqnref{eq: max constraint form} compactly. Recall that
we wish to represent the non-linear constraint $\phi \geq
\max_{\bx} \sum_i w_i \  c_i(\bx) - b(\bx)$, or equivalently,
$\phi \geq \max_\bx F^\bw(\bx)$, without generating one constraint
for each state as in \eqnref{eq: general constraint form}. The new, key
insight is that this non-linear constraint can be implemented
using a construction that follows the structure of variable
elimination in cost networks. %~\cite{Guestrin+al:2001b}.

Consider any function $e$ used within $\F$ (including the original
%RP3: scope
$f_i$'s), and let $\bZ$ be its scope.  For any assignment $\bz$
%RP1:
%
%value represents $u^e_\bz$.  For the initial functions $f_i^\bw$,
to $\bZ$, we introduce variable $u^e_\bz$, whose
value represents $e_\bz$, into the linear program.
For the initial functions $f_i^\bw$,
we include the constraint that $u^{f_i}_\bz = f_i^\bw(\bz)$.  As
$f_i^\bw$ is linear in $\bw$, 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,\ldots,e_L$ be the
%RP3: scope
functions extracted from $\F$, and let $\bZ$ be the scope of the
resulting $e$.  We introduce a set of constraints:
\beqn \label{eqn: factored LP max constraints}
  u^e_{\bz} \geq \sum_{j=1}^L u^{e_j}_{\instin{\bZ_j}{(\bz,x_l)}} ~~~~~ \forall x_l .
\eeqn 
%RP3: scope
Let $e_n$ be the last function generated in the
%RP3: scope
elimination, and recall that its scope is empty.  Hence, we have
only a single variable $u^{e_n}$.  We introduce the additional
constraint $\phi \geq u^{e_n}$.

%CGjairf1 added
The complete algorithm, presented in \figref{code:factored LP}, is
divided into three parts: First, we generate equality constraints for
functions that depend on the weights $w_i$ (basis functions). In the
second part, we add the equality constraints for functions that do not
depend on the weights (target functions).  These equality constraints let
us abstract away the differences between these two types of functions and
manage them in a unified fashion in the third part of the algorithm. This
third part follows a procedure similar to variable elimination described
in \figref{code:VE}. However, unlike standard variable elimination where
we would introduce a new function $e$, such that $e = \max_{x_l}
\sum_{j=1}^{L} e_j$, in our factored LP procedure we introduce new LP
variables $u^e_{\bz}$. To enforce the definition of $e$ as the maximum
over $X_l$ of $\sum_{j=1}^{L} e_j$, we introduce the new LP constraints
in \eqnref{eqn: factored LP max constraints}. 

%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Factored LP code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\begin{pseudocode}
FactoredLP ($C$, $\vecb$,$\order$) \\ 
\tab \codecomment{ $C = \{c_1,\ldots,c_k\}$ is the set of basis functions.}\\
\tab \codecomment{ $\vecb = \{b_1,\ldots,b_m\}$ is the set of target functions.}\\
\tab \codecomment{$\order$ stores the elimination order. } \\
\tab \codecomment{Return a (polynomial) set of constraints $\constraints$ equivalent to $\phi \geq \sum_i w_i c_i(\bx) + \sum_j b_j(\bx) , \forall \bx$ . }
\bcsone
\item \codecomment{Data structure for the constraints in factored LP.}
\item {\bf Let} $\constraints = \{\}$ .
\item \codecomment{Data structure for the intermediate functions generated in variable elimination.}
\item {\bf Let} $\F = \{\}$ .
\item \codecomment{Generate equality constraint to abstract away basis functions.}
\item {\bf For} each $c_i \in C$:
\bcstwo
\item {\bf Let} $\bZ = \Scope[c_i]$.
\item {\bf For} each assignment $\bz \in \bZ$, create a new LP variable $u^{f_i}_{\bz}$ and add a constraint to $\constraints$:
\bdm
u^{f_i}_\bz = w_i c_i(\bz).
\edm
%CGjairf2 fixed j -> i
\item {\bf Store} new function $f_i$ to use in variable elimination step: $\F = \F \cup \{f_i\}$.
\ecstwo
\item \codecomment{Generate equality constraint to abstract away target functions.}
\item {\bf For} each $b_j \in \vecb$:
\bcstwo
\item {\bf Let} $\bZ = \Scope[b_j]$.
\item {\bf For} each assignment $\bz \in \bZ$, create a new LP variable $u^{f_j}_{\bz}$ and add a constraint to $\constraints$:
\bdm
u^{f_j}_\bz = b_j(\bz).
\edm
\item {\bf Store} new function $f_j$ to use in variable elimination step: $\F = \F \cup \{f_j\}$.
\ecstwo
\item \codecomment{Now, $\F$ contains all of the functions involved in the LP, our constraints become: $\phi \geq \sum_{e_i \in \F} e_i(\bx)  , \forall \bx$ , which we represent compactly using a variable elimination procedure.}
\item {\bf For} $i = 1$ to number of variables:
\bcstwo
\item \codecomment{Select the next variable to be eliminated.}
\item {\bf Let} $l = \order(i)$ ; 
\item \codecomment{Select the relevant functions.}
\item
{\bf Let} $e_1,\dots ,e_L$ be the functions in  $\F$ whose scope
contains $X_l$, and let $\bZ_j = \Scope[e_j]$.
\item \codecomment{Introduce linear constraints for the maximum over current variable $X_l$.}
\item {\bf Define} a new function $e$ with scope $\bZ = \cup_{j=1}^L \bZ_j - \{X_l\}$
to represent $\max_{x_l} \sum_{j=1}^{L} e_j$.
\item {\bf Add} constraints to $\constraints$ to enforce maximum: for each assignment $\bz \in \bZ$:
\[
  u^e_{\bz} \geq \sum_{j=1}^L u^{e_j}_{\instin{\bZ_j}{(\bz,x_l)}} ~~~~~ \forall x_l .
\]
\item \codecomment{Update set of functions.}
\item {\bf Update} the set of functions $\F = \F \cup \{e\} \setminus \{e_1,\dots ,e_L\}$.
\ecstwo
\item \codecomment{Now, all variables have been eliminated and all functions have empty scope.}
\item {\bf Add} last constraint to $\constraints$:
\bdm
\phi \geq \sum_{e_i \in \F} e_i.
\edm
\item {\bf Return} $\constraints$.
\ecsone
\end{pseudocode}
%CGjairf2 fixed - -> +
\caption{Factored LP algorithm for the compact representation of the exponential set of constraints $\phi \geq \sum_i w_i c_i(\bx) + \sum_j b_j(\bx) , \forall \bx$.}\label{code:factored LP}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{example}
To understand this construction, consider our simple example
above, and assume we want to express the fact that $\phi \geq
\max_\bx F^\bw(\bx)$. 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 have four such variables.  We then
introduce a constraint defining the value of $u^{f_1}_{x_1,x_2}$
appropriately.  For example, for our $f_1$ above, we have
$u^{f_1}_{t,t} = 0$ and $u^{f_1}_{t,f} = w_1$.  We have
similar variables and constraints for each $f_j$ and each value
$\bz$ in $\bZ_j$. Note that each of the constraints is a simple
equality constraint involving numerical constants and perhaps the
weight variables $\bw$.

Next, we introduce variables for each of the intermediate
expressions
%CGjairf2 added and rephrased
generated by variable elimination. For example, when eliminating $X_4$,
we introduce a set of LP variables
$u^{e_1}_{x_2,x_3}$; for each of them, we have a set of
constraints
\[
  u^{e_1}_{x_2,x_3} \geq u^{f_3}_{x_2,x_4} + u^{f_4}_{x_3,x_4}
\]
one for each value $x_4$ of $X_4$.  We 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. \qed
\end{example}

We can now prove that our factored LP construction represents the
same constraint as non-linear constraint in \eqnref{eq: max
constraint form}:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theorem: factored LP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thm}\label{thm: factored LP}
The constraints generated by the factored LP construction are
equivalent to the non-linear constraint in \eqnref{eq: max
constraint form}. That is, an assignment to $(\phi,\bw)$ satisfies
the factored LP constraints if and only if it satisfies the
constraint in \eqnref{eq: max constraint form}.
\end{thm}
\proof See \appenref{proof:thm: factored LP}. \qed


%CGjairf1 changed a little and added umph until end of section
Returning to our original formulation, we have that $\sum_j
f_j^\bw$ is $C\bw-\vecb$ in the original set of constraints. Hence
our new set of constraints is equivalent to the original set:
$\phi \geq \max_{\bx} \sum_i w_i \  c_i(\bx) - b(\bx)$ in
\eqnref{eq: max constraint form}, which in turn is equivalent to
the exponential set of constraints $\phi \geq \sum_i w_i \
c_i(\bx) - b(\bx), \forall \bx$ in \eqnref{eq: general constraint
form}. Thus, we can represent this exponential set of constraints
by a new set of constraints and LP variables. The size of this new set,
as in variable elimination, is exponential only in the induced width of
the cost network, rather than in the total number of variables. 

In this section, we presented a new, general approach for compactly
representing ex\-po\-nen\-tial\-ly-large sets of LP constraints in problems with
factored structure. In the remainder of this paper, we exploit this
construction to design efficient planning algorithms for factored MDPs. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Factored Max-norm Projection} \label{sec: Factored Max-norm}

We can now use our procedure for representing the exponential
number of constraints in \eqnref{eq: general constraint form}
compactly to compute efficient max-norm projections, as in
\eqnref{eqn:minimax}:
\[
 \bwstar \in \arg\min_\bw \maxnorm{C \bw - \vecb}.
\]

The max-norm projection is computed by the linear program
in~(\ref{eqn:max-norm linear-program}). There are two sets of
constraints in this LP: $\phi  \geq  \sum_{j=1}^k c_{ij} w_j -
b_i, \forall i$ and $\phi  \geq  b_i - \sum_{j=1}^k c_{ij} w_j,
\forall i$. Each of these sets is an instance of the constraints
in \eqnref{eq: general constraint form}, which we have just
addressed in the previous section. Thus, if each of the $k$ basis
%RP3: scope
functions in $C$ is a restricted scope function and the target
%RP3: scope
function $\vecb$ is the sum of restricted scope functions, then
we can use our factored LP technique to represent the constraints
in the max-norm projection LP compactly. The correctness of our
algorithm is a corollary of \thmref{thm: factored LP}:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corollary: factored max-norm projection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{cor}\label{cor: factored max-norm projection}
%CGff1 made statement more precise
The solution $(\phistar,\bwstar)$ of a linear program that
minimizes $\phi$ subject to the constraints in
\textsc{FactoredLP}($C$, $-\vecb$,$\order$) and
\textsc{FactoredLP}($-C$, $\vecb$,$\order$), for any elimination
order $\order$ satisfies:
\bdm \bwstar \in \arg\min_\bw \maxnorm{C \bw - \vecb},
\mbox{~~~~and~~~~ } \phistar = \min_\bw \maxnorm{C \bw - \vecb}.
\qedmath
\edm
\commentout{
The solution $(\phistar,\bwstar)$ of our factored max-norm
projection linear program satisfies: \bdm \bwstar \in \arg\min_\bw
\maxnorm{C \bw - \vecb}, \mbox{~~~~and~~~~ } \phistar = \min_\bw
\maxnorm{C \bw - \vecb}. \qedmath\edm}
\end{cor}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The original max-norm projection LP had $k + 1$ variables and two
constraints for each state $\bx$; thus, the number of constraints is
exponential in the number of state variables. On the other hand, our new
factored max-norm projection LP has more variables, but exponentially fewer
constraints. The number of variables and constraints in the new factored LP
is exponential only in the number of state variables in the largest factor
in the cost network, rather than exponential in the total number of state
variables.  As we show in \secref{sec:Experimental results}, this
exponential gain allows us to compute max-norm projections efficiently when
solving very large factored MDPs.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Approximate Linear Programming}\label{sec:ALP}


%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Factored approximate linear program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\begin{pseudocode}
FactoredALP ($P$, $R$, $\discount$, $\Basis$, $\order$, $\sweights$) \\ 
\tab \codecomment{$P$ is the factored transition model.}\\
\tab \codecomment{$R$ is the set of factored reward functions.}\\
\tab \codecomment{$\discount$ is the discount factor.}\\
\tab \codecomment{$\Basis$ is the set of basis functions $\Basis = \{\basis_1, \ldots, \basis_k\}$.}\\
\tab \codecomment{$\order$ stores the elimination order. } \\
\tab \codecomment{$\sweights$ are the state relevance weights. } \\
\tab \codecomment{Return the basis function weights $\bw$ computed by approximate linear programming.}
\bcsone
\item \codecomment{Cache the backprojections of the basis functions.}
\item {\bf For} each basis function $\basis_i \in \Basis$; for each action $a$:
\bcstwo
\item {\bf Let} $\backbasis_i^a = \bprojfun{a}{\basis_i}$.
\ecstwo
\item \codecomment{Compute factored state relevance weights. }
\item {\bf For} each basis function $\basis_i$, compute the factored state relevance weights $\sweights_i$ as in \eqnref{eqn: factored state relevance weights} .
\item \codecomment{Generate approximate linear programming constraints }
\item {\bf Let} $\constraints = \{\}$.
\item {\bf For} each action $a$:
\bcstwo
\item {\bf Let} $\constraints = \constraints \cup  \mbox{FactoredLP}(\{\discount\backbasis_1^{a}-\basis_1,\ldots,\discount\backbasis_k^{a}-\basis_k\}, R^a, \order)$. 
\ecstwo
%RP4:  I was going to change this to add to... to match figure 5, but
%      decided to leave it b/c ``add to'' is ugly
\item \codecomment{So far, our constraints guarantee that $\phi \geq R(\bx,a) + \discount \sum_{\bx'} P(\bx' \mid \bx,a) \sum_i w_i \  \basis_i(\bx') - \sum_i w_i \  \basis_i(\bx)$; to satisfy the approximate linear programming solution in~(\ref{eq: explicit approx LP for MDPs}) we must add a final constraint.}
\item {\bf Let} $\constraints = \constraints \cup \{\phi = 0 \}$.
\item \codecomment{We can now obtain the solution weights by solving an LP.}
\item {\bf Let} $\bw$ be the solution of the linear program: minimize 
%CGjairf2 fixed \basis_i -> w_i
$\sum_i \sweights_i w_i$, subject to the constraints $\constraints$.
\item {\bf Return} $\bw$.
\ecsone
\end{pseudocode}
\caption{Factored approximate linear programming algorithm. }\label{code:factored ALP}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We begin with the simplest of our approximate MDP solution algorithms,
based on the approximate linear programming formulation in
\secref{sec: Approximate Linear Programming}. Using the basic
operations described in \secref{sec: basic ops}, we can formulate an
algorithm that is both simple and efficient.  

\subsection{The Algorithm}
As discussed in \secref{sec: Approximate Linear Programming},
%RP3: capitalizaiton
approximate linear program formulation is based on the linear
programming approach to solving MDPs presented in 
%CGjairf2 fixed ref
\secref{sec: Approximate Linear Programming}. However, in this approximate version, we restrict the space
of value functions to the linear space defined by our basis
functions. More precisely, in this approximate LP formulation, the
variables are $w_1,\ldots,w_k$ --- the weights for our basis
functions. The LP is given by: \beqn \label{eq: explicit approx LP
for MDPs}
\begin{array}{lll}
  \mbox{Variables:} & {w_1,\ldots,w_k} \ ; \\
  \mbox{Minimize: } & \sum_{\bx} \sweights(\bx) \sum_i w_i \  \basis_i(\bx) \ ; \\
  \mbox{Subject to:} & \sum_i w_i \  \basis_i(\bx) \geq
       [R(\bx,a) + \discount \sum_{\bx'} P(\bx' \mid \bx,a) \sum_i w_i \  \basis_i(\bx')] &
\forall \bx \in \bX, \forall a \in A .
\end{array}
\eeqn In other words, this formulation takes the LP in (\ref{eq:
exact LP for MDPs}) and substitutes the explicit state value
function with a linear value function representation $\sum_i w_i \
\basis_i(\bx)$. This transformation from an exact to an
approximate problem formulation has the effect of reducing the
number of free variables in the LP to $k$ (one for each basis
function coefficient), but the number of constraints remains
$\abs{\bX} \times |A|$.  In our {\it SysAdmin\/} problem, for example, the
number of constraints in the LP 
%CGjairf2 fixed ref
in~(\ref{eq: explicit approx LP for MDPs})
is $(m+1) \cdot 2^m$, where $m$ is the number of machines in the
network. However, using our algorithm for representing
exponentially large constraint sets compactly we are able to
compute the solution to this approximate linear programming
algorithm in 
%CGjairf2 emphed
\emph{closed form} with an exponentially smaller LP, as in
\secref{sec: Representing exponentially many constraints}.

First, consider the objective function $\sum_{\bx} \sweights(\bx)
\sum_i w_i \  \basis_i(\bx)$ of the LP (\ref{eq: explicit approx
LP for MDPs}). Naively representing this objective function
requires a summation over a exponentially large state space.
However, we can rewrite the objective and obtain a compact
representation. We first reorder the terms: \bdm \sum_{\bx}
\sweights(\bx) \sum_i w_i \ \basis_i(\bx) = \sum_i w_i \sum_{\bx}
\sweights(\bx) \ \basis_i(\bx). \edm 
Now, consider the state relevance weights $\sweights(\bx)$ as a
distribution over states, so that $\sweights(\bx) > 0$ and $\sum_{\bx}
\sweights(\bx) = 1$. As in backprojections, we can now write: 
\beqn \label{eqn: factored state relevance weights} 
\sweights_i = \sum_{\bx} \sweights(\bx) \ \basis_i(\bx) = \sum_{ \bc_i
\in \Cluster_i} \sweights(\bc_i) \ \basis_i(\bc_i); 
\eeqn 
where $\sweights(\bc_i)$
represents the marginal of the state relevance weights $\sweights$ over
the domain $\Val[\Cluster_i]$ of the basis function $\basis_i$. For example,
if we use uniform state relevance weights as in our experiments ---
$\sweights(\bx) = \frac{1}{\abs{\bX}}$ --- then the marginals
become $\sweights(\bc_i) = \frac{1}{\abs{\Cluster_i}}$. Thus, we
can rewrite the objective function as $\sum_i w_i \  \sweights_i$,
where each basis weight $\sweights_i$ is computed as shown in
\eqnref{eqn: factored state relevance weights}. If the state
relevance weights are represented by marginals, then the cost of
computing each $\sweights_i$ depends exponentially on the size of the
%RP3 scope
scope of $\Cluster_i$ only, rather than exponentially on the number of
state variables. On the other hand, if the state relevance weights
are represented by arbitrary distributions, we need to obtain
the marginals over the $\Cluster_i$'s, which may not be an
%RPfinal:  ``best results are'' -> ``greatest efficiency is''
%          to remove ambiguity between solution quality and efficienty
efficient computation. Thus, greatest efficiency is achieved by
using a compact representation, such as a Bayesian network,
for the state relevance weights.

Second, note that the right side of the constraints in the LP
(\ref{eq: explicit approx LP for MDPs}) correspond to the $Q_a$
functions: \bdm Q_a(\bx) = R^a(\bx) + \discount \sum_{\bx'} P(\bx'
\mid \bx,a) \sum_i w_i \  \basis_i(\bx'). \edm Using the efficient
backprojection operation in factored MDPs described in
\secref{sec: one-step lookahead} we can rewrite the $Q_a$
functions as: \bdm Q_a(\bx) = R^a(\bx) + \discount \sum_i w_i \
\backbasis^a_i(\bx); \edm where $\backbasis^a_i$ is the
backprojection of basis function $\basis_i$ through the transition
%RP3: scope
model $P_a$. As we discussed, if $\basis_i$ has scope restricted
to $\Cluster_i$, then $\backbasis^a_i$ is a restricted
%RP3: scope
scope function of $\backproj{a}{\Cluster_i'}$.

We can precompute the backprojections $\backbasis^a_i$ and
the basis relevance weights $\sweights_i$. The approximate linear
programming LP of (\ref{eq: explicit approx LP for MDPs}) can be
written as: 
\beqn \label{eq: factored approx LP for MDPs}
\begin{array}{lll}
  \mbox{Variables:} & {w_1,\ldots,w_k} \ ; \\
  \mbox{Minimize: } & \sum_i \sweights_i \  w_i \ ; \\
  \mbox{Subject to:} & \sum_i w_i \  \basis_i(\bx) \geq
       [R^a(\bx) + \discount \sum_i w_i \  \backbasis^a_i(\bx)] &
\forall \bx \in \bX, \forall a \in A .
\end{array}
\eeqn

Finally, we can rewrite this LP to use constraints of the same
form as the one in \eqnref{eq: max constraint form}: 
%CGjairf2 reversed \basis_i(\bx) - \discount \backbasis^a_i(\bx)
\beqn
\label{eq: max factored approx LP for MDPs}
\begin{array}{lll}
  \mbox{Variables:} & {w_1,\ldots,w_k} \ ; \\
  \mbox{Minimize: } & \sum_i \sweights_i \  w_i \ ; \\
  \mbox{Subject to:} & 0 \geq \max_{\bx} \left \{ R^a(\bx) + \sum_i w_i \ \left[  \discount \backbasis^a_i(\bx) - \basis_i(\bx)\right]  \right \} &
\forall a \in A.
\end{array}
\eeqn We can now use our factored LP construction in \secref{sec:
Representing exponentially many constraints} to represent these
non-linear constraints compactly. Basically, there is one set
of factored LP constraints for each action $a$. Specifically, we
%RP3: those in
can write the non-linear constraint in the same form as those in
\eqnref{eq: max constraint form} by expressing the functions $C$
as: $c_i(\bx) = \basis_i(\bx) - \discount \backbasis^a_i(\bx)$.
%RP3: scope
Each $c_i(\bx)$ is a restricted scope function; that is, if
%RP3: scope
$\basis_i(\bx)$ has scope restricted to $\Cluster_i$, then
%RP3: scope
$\backbasis^a_i(\bx)$ has scope restricted to
$\backproj{a}{\Cluster_i'}$, which means that $c_i(\bx)$ has
%RP3: scope
scope restricted to $\Cluster_i \cup \backproj{a}{\Cluster_i'}$.
Next, the target function $\vecb$ becomes the reward function
$R^a(\bx)$ which, by assumption, is factored. Finally, in the
constraint in \eqnref{eq: max constraint form}, $\phi$ is a free
variable. On the other hand, in the LP in~(\ref{eq: max factored
approx LP for MDPs}) the maximum in the right hand side must be
less than zero. This final condition can be achieved by adding a
constraint $\phi = 0$. Thus, our algorithm generates a set of
factored LP constraints, one for each action. The total number of
constraints and variables in this new LP is linear in the
number of actions $\abs{A}$ and only exponential in the induced
width of each cost network, rather than in the total number of
variables. The complete factored approximate linear programming algorithm is outlined in \figref{code:factored ALP}.

%CGjairf2 moved this discussion to later
%RP3: the
\commentout{
Note that, in the approximate linear programming algorithm, we never need to compute the decision list policy. The policy is always represented implicitly by the $Q_a$ functions. Thus, this algorithm does not require explicit computation or manipulation of the greedy policy. This difference has two important consequences: one computational and the other in terms of generality. First, not having to compute or consider the decision lists makes this algorithm faster and easier to implement. In terms of representation, we do not require the policies to be compact; thus, we do not need to make the default action assumption.  Therefore, this algorithm can deal with a more general class of problems, where each action can have its own independent DBN transition model.
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{An Example}\label{sec: example LP}
%CGjair1 added whole section

We now present a complete example of the operations required
by the approximate LP algorithm to solve the factored MDP shown in
\figref{fig:ring example}(a). Our presentation follows four steps:
problem representation, basis function selection, backprojections and LP
construction. 

%CGjair2 added \paragraphs in this section
\paragraph{Problem Representation: }
First, we must fully specify the factored MDP model for the problem. The
structure of the DBN is shown in \figref{fig:ring example}(b). This
structure is maintained for all action choices. Next, we must define the
transition probabilities for each action. There are 5 actions in this
problem: do nothing, or reboot one of the 4 machines in the network. The
CPDs for these actions are shown in \figref{fig:ring
example}(c). Finally, we must define the reward function. We decompose
the global reward as the sum of 4 local reward functions, one for each
machine, such that there is a reward if the machine is
working. Specifically, $R_i(X_i = \true) = 1$ and $R_i(X_i = \false) =
0$, breaking symmetry by setting $R_4(X_4 = \true) = 2$. We use a
discount factor of $\discount = 0.9$. 

\paragraph{Basis Function Selection: } 
%The next step is the selection of the basis functions. 
In this simple
example, we use five simple basis functions. First, we include the
constant function $\basis_0 = 1$. Next, we add indicators for each
machine which take value 1 if the machine is working: $\basis_i(X_i =
\true) = 1$ and $\basis_i(X_i = \false) = 0$.

\paragraph{Backprojections: } The first algorithmic step is computing the backprojection of the basis functions, as defined in \secref{sec: one-step lookahead}. The backprojection of the constant basis is simple:
\beqa
\backbasis^a_0 & = & \sum_{\bx'} P_a(\bx' \mid \bx) \basis_0 \ ;\\
& = & \sum_{\bx'} P_a(\bx' \mid \bx) \ 1 \ ;\\
& = & 1 \ .
\eeqa
Next, we must backproject our indicator basis functions $\basis_i$:
\beqa
\backbasis^a_i & = & \sum_{\bx'} P_a(\bx' \mid \bx) \basis_i(x_i') \ ;\\
& = & \sum_{x_1',x_2',x_3',x_4'} \prod_j P_a(x_j' \mid x_{j-1},x_{j}) \basis_i(x_i') \ ;\\
& = & \sum_{x_i'} P_a(x_i' \mid x_{i-1},x_{i}) \basis_i(x_i')
\sum_{\instin{\bX'-\{X_i'\}}{\bx'}} \ \prod_{j\neq i} P_a(x_j' \mid x_{j-1},x_{j}) \ ;\\
& = & \sum_{x_i'} P_a(x_i' \mid x_{i-1},x_{i}) \basis_i(x_i') \ ;\\
& = & P_a(X_i' = \true \mid x_{i-1},x_{i}) \ 1 + P_a(X_i' = \false \mid x_{i-1},x_{i}) \ 0 \ ;\\
& = & P_a(X_i' = \true \mid x_{i-1},x_{i}) \ .
\eeqa
Thus, $\backbasis^a_i$ is a restricted scope function of $\{X_{i-1},X_{i}\}$. We can now use the CPDs in \figref{fig:ring example}(c) to specify $\backbasis^a_i$:
\beqa
\backbasis^{\fmbox{\it reboot $=i$}}_i(X_{i-1},X_i) & = & 
\begin{array}{|c|c|c|}
\hline & X_i = \true & X_i = \false 
\\
\hline
X_{i-1} = \true  & 1 & 1 \\
\hline
X_{i-1} = \false  & 1 & 1 \\
\hline
\end{array} \ ;\\ & \mbox{ } & \\
\backbasis^{\fmbox{\it reboot $\neq i$}}_i(X_{i-1},X_i) & = & 
\begin{array}{|c|c|c|}
\hline & X_i = \true & X_i = \false 
\\
\hline
X_{i-1} = \true  & 0.9 & 0.09 \\
\hline
X_{i-1} = \false  & 0.5 & 0.05 \\
\hline
\end{array} \ .
\eeqa

\paragraph{LP Construction: } To illustrate the factored LPs constructed by our algorithms, we define
the constraints for the approximate linear programming approach presented
above. 
%CGjair2 changed signs of c_i to make them positive in the constraint
First, we define the functions $c_i^a = \discount \backbasis_i^a - \basis_i$, as shown in \eqnref{eq: max factored approx LP
for MDPs}. 
In our example, these functions are $c_0^a = \discount - 1
= -0.1$ for the constant basis, and for the indicator bases: 
\beqa
c^{\fmbox{\it reboot $=i$}}_i(X_{i-1},X_i) & = & 
\begin{array}{|c|c|c|}
\hline & X_i = \true & X_i = \false 
\\
\hline
X_{i-1} = \true  & -0.1 & 0.9 \\
\hline
X_{i-1} = \false  & -0.1 & 0.9 \\
\hline
\end{array} \ ;\\ & \mbox{ } & \\
c^{\fmbox{\it reboot $\neq i$}}_i(X_{i-1},X_i) & = & 
\begin{array}{|c|c|c|}
\hline & X_i = \true & X_i = \false 
\\
\hline
X_{i-1} = \true  & -0.19 & 0.081 \\
\hline
X_{i-1} = \false  & -0.55 & 0.045 \\
\hline
\end{array} \ .
\eeqa
Using this definition of $c_i^a$, the approximate linear programming constraints are given by:
\begin{equation}\label{eq: eg max constraint}
0 \geq \max_{\bx} \sum_i R_i  + \sum_j w_j c_j^a \ \ , \ \forall a \ .
\end{equation}
We present the LP construction for one of the 5 actions:
%CGjairf2 made actions italic in this section
$\mbox{\small \it reboot $=1$}$. Analogous constructions can be made for the other
actions.

In the first set
%CGjairf2 added
of constraints, we abstract away the difference between rewards and
basis functions by introducing LP variables $u$ and equality
constraints.  We begin with the reward functions:
\beqa
u_{x_1}^{R_1}  =  1 \ , \ 
u_{\bar{x}_1}^{R_1}  =  0 \ ; \tab
u_{x_2}^{R_2}  =  1 \ , \ 
u_{\bar{x}_2}^{R_2}  =  0 \ ; \\
u_{x_3}^{R_3}  =  1 \ , \ 
u_{\bar{x}_3}^{R_3}  =  0 \ ; \tab
u_{x_4}^{R_4}  =  2 \ , \ 
u_{\bar{x}_4}^{R_4}  =  0 \ . \\
\eeqa
We now represent the equality constraints for the $c_j^a$ functions for the
$\mbox{\small \it reboot $=1$}$ action. 
%CGjairf2 rephrased
Note that the appropriate basis function
weight from \eqnref{eq: eg max constraint} appears in these constraints:
\bdm
\begin{array}{cccc}
u^{c_0} = - 0.1 \ w_0 \ ; \\
u_{x_1,x_4}^{c_1} = - 0.1 \ w_1 \ , & 
u_{\bar{x}_1,x_4}^{c_1} = 0.9 \ w_1 \ , & 
u_{x_1,\bar{x}_4}^{c_1} = - 0.1 \ w_1 \ , & 
u_{\bar{x}_1,\bar{x}_4}^{c_1} = 0.9 \ w_1 \ ; \\
u_{x_1,x_2}^{c_2} = - 0.19 \ w_2 \ , & 
u_{\bar{x}_1,x_2}^{c_2} = -0.55 \ w_2 \ , &
u_{x_1,\bar{x}_2}^{c_2} = 0.081 \ w_2 \ , &
u_{\bar{x}_1,\bar{x}_2}^{c_2} = 0.045 \ w_2 \ ; \\
u_{x_2,x_3}^{c_3} = - 0.19 \ w_3 \ , & 
u_{\bar{x}_2,x_3}^{c_3} = -0.55 \ w_3 \ , &
u_{x_2,\bar{x}_3}^{c_3} = 0.081 \ w_3 \ , &
u_{\bar{x}_2,\bar{x}_3}^{c_3} = 0.045 \ w_3 \ ; \\
u_{x_3,x_4}^{c_4} = - 0.19 \ w_4 \ , & 
u_{\bar{x}_3,x_4}^{c_4} = -0.55 \ w_4 \ , & 
u_{x_3,\bar{x}_4}^{c_4} = 0.081 \ w_4 \ , &
u_{\bar{x}_3,\bar{x}_4}^{c_4} = 0.045 \ w_4 \ . \\
\end{array}
\edm
Using these new LP variables, our LP constraint from \eqnref{eq: eg max
constraint} for the $\mbox{\small \it reboot $=1$}$ action becomes:
\bdm
0 \geq \max_{X_1,X_2,X_3,X_4} \sum_{i=1}^4 u^{R_i}_{X_i}  + u^{c_0} + \sum_{j=1}^4  u^{c_j}_{X_{j-1},X_j} \ .
\edm
We are now ready for the variable elimination process. We illustrate the elimination of variable $X_4$:
\bdm
0 \geq \max_{X_1,X_2,X_3} \sum_{i=1}^3 u^{R_i}_{X_i}  + u^{c_0} + \sum_{j=2}^3  u^{c_j}_{X_{j-1},X_j} + \max_{X_4}\left[ u^{R_4}_{X_4}  + u^{c_1}_{X_{1},X_4} + u^{c_4}_{X_{3},X_4} \right] \ .
\edm
We can represent the term $\max_{X_4}\left[ u^{R_4}_{X_4}  + u^{c_1}_{X_{1},X_4} + u^{c_4}_{X_{3},X_4} \right]$ by a set of linear constraints, one for each assignment of $X_1$ and $X_3$, using the new LP variables $u^{e_1}_{X_1,X_3}$ to represent this maximum:
\beqa
u^{e_1}_{x_1,x_3} & \geq & u^{R_4}_{x_4}  + u^{c_1}_{x_{1},x_4} + u^{c_4}_{x_{3},x_4} \ ; \\
u^{e_1}_{x_1,x_3} & \geq & u^{R_4}_{\bar{x}_4}  + u^{c_1}_{x_{1},\bar{x}_4} + u^{c_4}_{x_{3},\bar{x}_4} \ ; \\
u^{e_1}_{\bar{x}_1,x_3} & \geq & u^{R_4}_{x_4}  + u^{c_1}_{\bar{x}_{1},x_4} + u^{c_4}_{x_{3},x_4} \ ; \\
u^{e_1}_{\bar{x}_1,x_3} & \geq & u^{R_4}_{\bar{x}_4}  + u^{c_1}_{\bar{x}_{1},\bar{x}_4} + u^{c_4}_{x_{3},\bar{x}_4} \ ; \\
u^{e_1}_{x_1,\bar{x}_3} & \geq & u^{R_4}_{x_4}  + u^{c_1}_{x_{1},x_4} + u^{c_4}_{\bar{x}_{3},x_4} \ ; \\
u^{e_1}_{x_1,\bar{x}_3} & \geq & u^{R_4}_{\bar{x}_4}  + u^{c_1}_{x_{1},\bar{x}_4} + u^{c_4}_{\bar{x}_{3},\bar{x}_4} \ ; \\
u^{e_1}_{\bar{x}_1,\bar{x}_3} & \geq & u^{R_4}_{x_4}  + u^{c_1}_{\bar{x}_{1},x_4} + u^{c_4}_{\bar{x}_{3},x_4} \ ; \\
u^{e_1}_{\bar{x}_1,\bar{x}_3} & \geq & u^{R_4}_{\bar{x}_4}  + u^{c_1}_{\bar{x}_{1},\bar{x}_4} + u^{c_4}_{\bar{x}_{3},\bar{x}_4} \ . \\
\eeqa

\begin{figure}[t]
\centerline
{\epsfig{figure=Figures/explicit-versus-factored-ring.eps,height=5.2cm}}
\caption{Number of constraints in the LP generated by the explicit state representation versus the factored LP construction for the solution of the ring problem with basis functions over single variables and approximate linear programming as the solution algorithm. } \label{fig:explicit versus factored LP constraints}
\end{figure}
We have now eliminated variable $X_4$ and our global non-linear constraint becomes:
\bdm
0 \geq \max_{X_1,X_2,X_3} \sum_{i=1}^3 u^{R_i}_{X_i}  + u^{c_0} + \sum_{j=2}^3  u^{c_j}_{X_{j-1},X_j} + u^{e_1}_{X_1,X_3} \ .
\edm
Next, we eliminate variable $X_3$. The new LP constraints and variables have the form:
\bdm
u^{e_2}_{X_1,X_2} \geq u^{R_3}_{X_3}  + u^{c_3}_{X_{2},X_3} + u^{e_1}_{X_1,X_3} \ , \ \forall \ X_1, X_2, X_3 \ ;
\edm
thus, removing $X_3$ from the global non-linear constraint:
\bdm
0 \geq \max_{X_1,X_2} \sum_{i=1}^2 u^{R_i}_{X_i}  + u^{c_0} + u^{c_2}_{X_{1},X_2} + u^{e_2}_{X_1,X_2} \ .
\edm
We can now eliminate $X_2$, generating the linear constraints:
\bdm
u^{e_3}_{X_1} \geq u^{R_2}_{X_2}  + u^{c_2}_{X_{1},X_2} + u^{e_2}_{X_1,X_2} \ , \ \forall \ X_1, X_2 \ .
\edm
Now, our global non-linear constraint involves only $X_1$:
\bdm
0 \geq \max_{X_1} u^{R_1}_{X_1}  + u^{c_0} + u^{e_3}_{X_1} \ .
\edm
As $X_1$ is the last variable to be eliminated, the scope of the new LP
variable is empty and the linear constraints are given by:
\bdm
u^{e_4} \geq u^{R_1}_{X_1}  + u^{e_3}_{X_1} \ , \ \forall \ X_1 \ .
\edm
All of the state variables have now been eliminated, turning our global
non-linear constraint into a simple linear constraint:
\bdm
0 \geq u^{c_0} + u^{e_4} \ ,
\edm
which completes the LP description for the approximate linear programming solution to the problem in \figref{fig:ring example}.

In this small example with only four state variables, our factored LP
technique generates a total of 89 equality constraints, 115 inequality
constraints and 149 LP variables, while the explicit state representation
in \eqnref{eq: approx LP for MDPs} generates only 80 inequality
constraints and 5 LP variables. However, as the problem size increases, the
number of constraints and LP variables in our factored LP approach grow
as $O(n^2)$,  while the explicit state approach grows exponentially,
at $O(n2^n)$. This scaling effect is illustrated in \figref{fig:explicit
versus factored LP constraints}. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Approximate Policy Iteration with Max-norm Projection}\label{sec:API}

%CGff1 new lead-in

The factored approximate linear programming approach described in the
previous section is both elegant and easy to implement. However,
we cannot, in general, provide strong guarantees about the error
it achieves. 
An alternative is to use the approximate policy iteration described in
\secref{sec: Approximate Policy Iteration}, which does offer certain
bounds on the error.  However, as we shall see, this algorithm is
significantly more complicated, and requires that we place additional
restrictions on the factored MDP.


In particular, approximate policy iteration requires a
representation of the policy at each iteration. In order to obtain
a compact policy representation, we must make an additional
assumption: each action only affects a small number of state
variables. We first state this assumption formally. Then, we show
how to obtain a compact representation of the greedy policy with
respect to a factored value function, under this assumption.
Finally, we describe our factored approximate policy iteration
algorithm using max-norm projections.










%CGjairf2 added subsection here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Default Action Model} \label{sec: Default action model}

In \secref{sec: Factored MDPs}, we presented the factored MDP model, where each action is associated with its own factored transition model represented as a DBN and with its own factored reward function. However, different
actions often 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~\cite{Boutilier+al:2000}. 

This type of structure turns out to be useful for compactly representing
policies, a property which is important in our approximate policy
iteration algorithm.  Thus, in this section of the paper, we restrict
attention to factored MDPs that are defined using a \emph{default
transition model} $\trans_d = \langle\graph_d,P_d\rangle$
\cite{Koller+Parr:2000}.  For each action $a$, we define $\effects{a}
\subseteq \bX'$ to be the variables in the next state 
whose local probability model is different from $\trans_d$, \ie, 
those variables $X_i'$ such that $P_a(X_i' \mid \Parents_a(X_i'))
\neq P_d(X_i' \mid \Parents_d(X_i'))$.

\begin{example}
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 transition
model.  The transition model for $a_i$ is different from $d$ only
in the transition model for the variable $X_i'$, which is now
$X_i'=\true$ with probability one, regardless of the status of the
neighboring machines. Thus, in this example, $\effects{a_i} = X_i'$. \qed
\end{example}

As in the transition dynamics, we can also define the notion of
\emph{default reward model}. In this case, there is a set of reward
functions $\sum_{i=1}^r R_i(\rclust_i)$ associated with the default
action $d$. In addition, each action $a$ can have a reward function
$R^a(\rclust^a)$. Here, the extra reward of action $a$ has scope
restricted to $\rewards{a} = \rclust_i^a \subset
\{X_1,\ldots,X_n\}$. Thus, the total reward associated with action $a$
is given by $R^a + \sum_{i=1}^r R_i$. Note that $R^a$ can also be
factored as a linear combination of smaller terms for an even more compact representation.

We can now build on this additional assumption to define the complete algorithm. 
Recall that the approximate policy iteration algorithm iterates through
two steps: policy improvement and approximate value determination.  We
now discuss each of these steps.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Computing Greedy Policies} \label{sec: Computing greedy policies}

%RP3: of
%RP3: Reworded first sentence
The policy improvement step computes the greedy policy relative
to a value function $\V\at{t-1}$: 
\[
\policy\at{t} = \greedy{\V\at{t-1}}.
\]  
Recall that our value function estimates have the linear form $\Basismatrix \bw$. 
As we described in \secref{sec: one-step lookahead}, the
greedy policy for this type of value function is given by: 
\bdm 
\greedy{\Basismatrix \bw}(\bx) = \arg\max_a Q_a(\bx), 
\edm 
where each $Q_a$ can be represented by: $Q_a(\bx) =  R(\bx,a) + \sum_i w_i \ \backbasis^a_i(\bx). $

If we attempt to represent this policy naively, we are again faced with
the problem of exponentially large state spaces. Fortunately, as shown by
Koller and Parr~\citeyear{Koller+Parr:2000}, the greedy policy relative to a
factored value function has the form of a \emph{decision list}.
More precisely, the policy can be written in the form
%\beqn\label{eq:declist}
$\langle\event_1,a_1\rangle, \langle\event_2,a_2\rangle,\ldots,
     \langle \event_L,a_L\rangle$,
%\eeqn
where each $\event_i$ is an assignment of values to some small
subset $\Event_i$ of variables, and each $a_i$ is an action.  The
greedy action to take in state $\bx$ is the action $a_j$
corresponding to the first event $\event_j$ in the list with which
$\bx$ is consistent. For completeness, we now review the
construction of this decision-list policy.

%RP3:  key -> critical
The critical assumption that allows us to represent the policy as a
compact decision list is the default action assumption described
in \secref{sec: Default action model}. Under this assumption, the $Q_a$
functions can be written as: \bdm Q_a(\bx) =  R^a(\bx) +
\sum_{i=1}^r R_i(\bx) + \sum_i w_i \  \backbasis^a_i(\bx), \edm
%RP3: scope
where $R^a$ has scope restricted to $\rclust^a$. The $Q$ function
%RP3:  being -> is
for the default action $d$ is just: $Q_d(\bx) = \sum_{i=1}^r
R_i(\bx) + \sum_i w_i \ \backbasis^d_i(\bx)$.

We now have a set of linear $Q$-functions which implicitly
describes a policy $\pi$.  It is not immediately obvious that
these $Q$ functions result in a compactly expressible policy.
%RP3: The key -> an important
An important insight is that most of the components in the weighted
combination are identical, so that $\backbasis^a_i$ is
equal to $\backbasis^d_i$ for most $i$. Intuitively, a component
$\backbasis^a_i$ corresponding to the backprojection of basis
function $\basis_i(\Cluster_i)$ is only different if the
action $a$ influences one of the variables in $\Cluster_i$.  More
formally, assume that $\effects{a} \cap \Cluster_i = \emptyset$.
In this case, all of the variables in $\Cluster_i$ have the same
transition model in $\trans_a$ and $\trans_d$. Thus, we have that
$\backbasis^a_i(\bx) = \backbasis^d_i(\bx)$; in other words, the $i$th 
component of the $Q_a$ function is irrelevant when deciding whether
action $a$ is better than the default action $d$. We can define
which components are actually relevant: let $\relcomp_a$ be the
set of indices $i$ such that $\effects{a} \cap \Cluster_i \neq
\emptyset$.  These are the indices of those basis functions whose
backprojection differs in $P_a$ and $P_d$. In our example DBN of
\figref{fig:ring example}, actions and basis functions involve
single variables, so $\relcomp_{a_i} = i$.

Let us now consider the impact of taking action $a$ over the
default action $d$. We can define the impact --- the difference
in value --- as: 
\begin{eqnarray}
 \diffun_a(\bx) & = & Q_a(\bx) - Q_d(\bx); \nonumber\\
                & = & R^a(\bx) + \sum_{i \in \relcomp_a} w_i \  \left[\backbasis^a_i(\bx) - \backbasis^d_i(\bx) \right].
\label{eqn: dlist bonus}
\end{eqnarray}
This analysis shows that $\diffun_a(\bx)$ is a function
%RP3: scope
whose scope is restricted to 
\beqn
\Diff_a = \rclust^a \cup
\left[\cup_{i \in \relcomp_a} \backproj{a}{\Cluster_i'}\right].
\label{eqn:dlist bonus domain}
\eeqn
In our example DBN, $\Diff_{a_2} = \{X_1, X_2\}$.


Intuitively, we now have a situation where we have a ``baseline''
value function $Q_d(\bx)$ which defines a value for each state
$\bx$.  Each action $a$ changes that baseline by adding or
%RP3:  key -> null
subtracting an amount from each state.  The point is that this
amount depends only on $\Diff_a$, so that it is the same for all
states in which the variables in $\Diff_a$ take the same values.

We can now define the 
%CGjairf2 optimal -> greedy
greedy policy relative to our $Q$
functions.  For each action $a$, define a set of
\emph{conditionals} $\langle\event,a,\delta\rangle$, where each
$\event$ is some assignment of values to the variables $\Diff_a$,
and $\delta$ is $\diffun_a(\event)$. Now, sort the
conditionals for all of the actions by order of decreasing
$\delta$: 
\bdm 
\langle\event_1,a_1,\delta_1\rangle,
\langle\event_2,a_2,\delta_2\rangle,\ldots,
     \langle \event_L,a_L,\delta_L\rangle .
\edm 
Consider our optimal action in a state $\bx$.  We would like
to get the largest possible ``bonus'' over the default value.  If
$\bx$ is consistent with $\event_1$, we should clearly take action
$a_1$, as it gives us bonus $\delta_1$.  If not, then we
should try to get $\delta_2$; thus, we should check if $\bx$ is
consistent with $\event_2$, and if so, take $a_2$.  Using this
procedure, we can compute the decision-list policy associated with
our linear estimate of the value function. 
%CGjairf1 added
The complete algorithm for computing the decision list policy is
summarized in \figref{code:dlist}. 

%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% decision list policy code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{pseudocode}
DecisionListPolicy ($Q_a$) \\ 
\tab \codecomment{$Q_a$ is the set of Q-functions, one for each action;}\\
\tab \codecomment{Return the decision list policy $\dlist$.}
\bcsone
\item \codecomment{Initialize decision list.}
\item {\bf Let} $\dlist = \{\}$.
\item \codecomment{Compute the bonus functions.}
\item {\bf For} each action $a$, other than the default action $d$:
\bcstwo
\item {\bf Compute} the bonus for taking action $a$,
\bdm
\diffun_a(\bx) =  Q_a(\bx) - Q_d(\bx) ;
\edm
as in \eqnref{eqn: dlist bonus}. Note that $\diffun_a$ has scope restricted to $\Diff_a$, as in \eqnref{eqn:dlist bonus domain}.
\item \codecomment{Add states with positive bonuses to the (unsorted) decision list.}
\item {\bf For} each assignment $\event \in \Diff_a$:
\bcsthree
\item {\bf If} $\diffun_a(\event) > 0$, add branch to decision list:
\bdm
\dlist = \dlist \cup \{\langle\event,a,\diffun_a(\event)\rangle\} .
\edm 
\ecsthree
\ecstwo
\item \codecomment{Add the default action to the (unsorted) decision list.}
\item {\bf Let} $\dlist = \dlist \cup \{\langle\emptyset,d,0\rangle\}$.
\item \codecomment{Sort decision list to obtain final policy.}
\item {\bf Sort} the decision list $\dlist$ 
%CGjairf2 added
in decreasing order
on the $\diffun$ element of %the tuples 
$\langle\event,a,\diffun\rangle$.
\item {\bf Return} $\dlist$.
\ecsone
\end{pseudocode}
\caption{Method for computing the decision list policy $\dlist$ from the factored representation of the $Q_a$ functions. }\label{code:dlist}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Note that the number of conditionals in the list is $\sum_{a}
|\Val(\Diff_a)|$; $\Diff_a$, in turn, depends on the set of basis
function clusters that intersect with the effects of $a$.  Thus, the
size of the policy depends in a natural way on the interaction between
the structure of our process description and the structure of our
basis functions.
%RP3: Added to address reviewer complaint
In problems where the actions modify a large number of variables, the
policy representation could become unwieldy.  The approximate linear
programming approach in \secref{sec:ALP} is more appropriate in
such cases, as it does not require an explicit representation of the
policy.

\subsection{Value Determination}
%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Factored approximate policy iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\begin{pseudocode}
FactoredAPI ($P$, $R$, $\discount$, $\Basis$, $\order$, $\eps$, $t_{max}$) \\ 
\tab \codecomment{$P$ is the factored transition model.}\\
\tab \codecomment{$R$ is the set of factored reward functions.}\\
\tab \codecomment{$\discount$ is the discount factor.}\\
\tab \codecomment{$\Basis$ is the set of basis functions $\Basis = \{\basis_1, \ldots, \basis_k\}$.}\\
\tab \codecomment{$\order$ stores the elimination order. } \\
\tab \codecomment{$\eps$ Bellman error precision. } \\
\tab \codecomment{$t_{max}$ maximum number of iterations. } \\
\tab \codecomment{Return the basis function weights $\bw$ computed by approximate policy iteration.}
\bcsone
\item \codecomment{Initialize weights}
\item {\bf Let} $\bw\at{0} = {\bf 0}$.
\item \codecomment{Cache the backprojections of the basis functions.}
\item {\bf For} each basis function $\basis_i \in \Basis$; for each action $a$:
\bcstwo
\item {\bf Let} $\backbasis_i^a = \bprojfun{a}{\basis_i}$.
\ecstwo
\item \codecomment{Main approximate policy iteration loop.}
\item {\bf Let} $t = 0$.
\item {\bf Repeat}
\bcstwo
\item \hspace*{-.5cm} \codecomment{\bf \emph{Policy improvement part of the loop.}}
\item \codecomment{Compute decision list policy for iteration $t$ weights.}
\item {\bf Let} $\dlist\at{t} = \mbox{DecisionListPolicy}(R^a + \discount \sum_i w_i\at{t} \backbasis_i^a )$.
\item \hspace*{-.5cm} \codecomment{\bf \emph{Value determination part of the loop.}}
\item \codecomment{Initialize constraints for max-norm projection LP.}
\item {\bf Let} $\constraints^{+} = \{\}$ and $\constraints^{-} = \{\}$.
\item \codecomment{Initialize indicators. }
\item {\bf Let} $\indicator = \{\}$.
\item \codecomment{For every branch of the decision list policy, generate the relevant set of constraints, and update the indicators to constraint the state space for future branches.}
\item {\bf For} each branch $\langle\event_j,a_j\rangle$ in the decision list policy $\dlist\at{t}$:
\bcsthree
\item \codecomment{Instantiate the variables in $\Event_j$ to the assignment given in $\event_j$.}
\item {\bf Instantiate} the set of functions $\{\basis_1-\discount\backbasis_1^{a_j},\ldots,\basis_k-\discount\backbasis_k^{a_j}\}$ with the partial state assignment $\event_j$ and store in $C$.
\item {\bf Instantiate} the target functions $R^{a_j}$ with the partial state assignment $\event_j$ and store in $\vecb$.
\item {\bf Instantiate} the indicator functions $\indicator$ with the partial state assignment $\event_j$ and store in $\indicator'$.
\item \codecomment{Generate the factored LP constraints for the current decision list branch.}
\item {\bf Let} $\constraints^{+} = \constraints^{+} \cup  \mbox{FactoredLP}(C, - \vecb + \indicator', \order)$. 
\item {\bf Let} $\constraints^{-} = \constraints^{-} \cup  \mbox{FactoredLP}(-C, \vecb + \indicator', \order)$. 
\item \codecomment{Update the indicator functions.}
\item {\bf Let} $\indicator_j(\bx) = - \infty \indicone(\bx = \event_j)$ and update the indicators $\indicator = \indicator \cup \indicator_j$. 
\ecsthree
\item \codecomment{We can now obtain the new set of weights by solving an LP, which corresponds to the max-norm projection.}
\item {\bf Let} $\bw\at{t+1}$ be the solution of the linear program: minimize $\phi$, subject to the constraints $\{\constraints^{+} , \constraints^{-}\}$.
\item {\bf Let} $t = t+1$.
\ecstwo
\item {\bf Until} $\Bellman(\Basismatrix \bw\at{t}) \leq \eps$ or $t \geq t_{max}$ or $\bw\at{t-1} = \bw\at{t}$.
\item {\bf Return} $\bw\at{t}$.
\ecsone
\end{pseudocode}
\caption{Factored approximate policy iteration with max-norm projection algorithm. }\label{code:factored API}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In the approximate value determination step our algorithm computes:
\bdm
 \bw\at{t} = \arg\min_\bw \maxnorm{\Basismatrix\bw -
   \left(R_{\policy\at{t}} + \discount P_{\policy\at{t}}\Basismatrix\bw\right)}.
\end{displaymath}
By rearranging the expression, we get: 
\bdm
 \bw\at{t} = \arg\min_\bw \maxnorm{\left(\Basismatrix - \discount
       P_{\policy\at{t}}\Basismatrix\right)\bw - R_{\policy\at{t}}}.
\end{displaymath}
This equation is an instance of the optimization in
\eqnref{eqn:minimax}.  If $P_{\policy\at{t}}$ is factored, we 
can conclude that $C = (\Basismatrix - \discount P_{\policy\at{t}}\Basismatrix)$ is also
a matrix whose columns correspond to restricted-scope functions.
More specifically: 
\bdm 
c_i(\bx) = \basis_i(\bx) - \discount \backbasis^{\policy\at{t}}_i(\bx),
\edm 
where $\backbasis^{\policy\at{t}}_i$ is the backprojection of the basis
function $\basis_i$ through the transition model
$P_{\policy\at{t}}$, as described in \secref{sec: one-step
lookahead}. The target $\vecb = R_{\policy\at{t}}$ corresponds to
the reward function, which for the moment is assumed to be
factored. Thus, we can again apply our factored LP in \secref{sec:
Factored Max-norm} to estimate the value of the policy
$\policy\at{t}$.

Unfortunately, the transition model $P_{\policy\at{t}}$ is not
factored,
%CGjairf2 added 
as a decision list representation for the policy $\policy\at{t}$ will,
in general, induce a transition model $P_{\policy\at{t}}$ which cannot
be represented by a compact DBN.  Nonetheless, we can still generate a
compact LP by exploiting the decision list structure of the policy.
%CGff1 more intuition
The basic idea is to introduce cost
networks corresponding to each branch in the decision list,
ensuring, additionally, that only states consistent with this
branch are considered in the cost network maximization.
Specifically, we have a factored LP construction for each branch
$\langle \event_i,a_i\rangle$. The $i$th cost network only
considers a subset of the states that is consistent with the $i$th
branch of the decision list. 
Let $S_i$ be the set of
states $\bx$ such that $\event_i$ is the first event in the decision
list for which $\bx$ is consistent. That is, for each state $\bx \in
S_i$, $\bx$ is consistent with $\event_i$, but it is \emph{not}
consistent with any $\event_j$ with $j < i$.

Recall that, as in \eqnref{eq: general constraint form}, our LP
construction defines a set of constraints that imply that $\phi
\geq \sum_i w_i \  c_i(\bx) - b(\bx)$ for each state $\bx$.
Instead, we 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_{\policy\at{t}}$.
Similarly, the reward function becomes $R^{a_i}(\bx) +
\sum_{i=1}^r R_i(\bx)$ for this subset of states.

%CGff1 added actual constraints generated by our approach
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 $\event_i$, but not to states that are
consistent with some $\event_j$ for $j < i$.  To guarantee the
first condition, we simply instantiate the variables in $\Event_i$
to take the values specified in $\event_i$.  That is, our cost
network now considers only the variables in
$\{X_1,\ldots,X_n\}-\Event_i$, and computes the maximum only over
the states consistent with $\Event_i=\event_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 $\indicator_j$ for each previous decision
$\event_j$, with weight $-\infty$. More specifically,
$\indicator_j$ is a function that takes value $-\infty$ for states
consistent with $\event_j$ and zero for other all assignments of
$\Event_j$. The constraints for the $i$th branch will be of the
form:
\beqn\label{eq: constraints with decision branch indicators}
\phi \geq R(\bx,a_i)+\sum_l w_l\left(\discount
\backbasis_l(\bx,a_i) - \basis(\bx)\right) + \sum_{j<i}
-\infty\indicone(\bx=\event_j), \tab \forall
\consist{\bx}{\event_i}, 
\eeqn
where $\consist{\bx}{\event_i}$ defines the assignments of $\bX$
consistent with $\event_i$. The introduction of these indicators
causes the constraints associated with $\event_i$ to be trivially
satisfied by states in $S_j$ for $j<i$. Note that each of these
indicators is a restricted-scope function of $\Event_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.  The complete approximate
policy iteration with max-norm projection algorithm is outlined in
\figref{code:factored API}.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Comparisons} \label{sec:PI comparisons}

It is instructive to compare our max-norm policy iteration
algorithm to the $\cL_2$-projection policy iteration algorithm of
Koller and Parr~\citeyear{Koller+Parr:2000} in terms of
computational costs per iteration and implementation complexity.
Computing the $\cL_2$ projection requires (among other things) a
series of dot product operations between basis functions and
backprojected basis functions
$\dotprod{\basis_i}{\backbasis^{\policy}_j}$.  These expressions
are easy to compute if $P_{\pi}$ refers to the transition model of
a particular action $a$.  However, if the policy $\pi$ 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 branch of the decision list, for every pair
of basis functions $i$ and $j$, and for each assignment to the
%RP3: \val -> \scope
variables in $\Scope[\basis_i] \cup \Scope[\backbasis^{a}_j]$, it
requires the solution of a counting problem which is $\sharp
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 branch of the decision list. This
makes the algorithm very 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 only one of the BN inferences in the counting problem for the
$\cL_2$ projection.  Overall, for each branch in the decision
list, we have a total of two of these ``inferences,'' as opposed
%RP3: \val -> \scope
to one for each assignment of $\Scope[\basis_i] \cup
\Scope[\backbasis^a_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 $\cL_2$-projection.  Furthermore, the use of linear programming
allows us to rely on existing LP packages (such as CPLEX), which
are very highly optimized.

%CGjairf2 rephrased and added stuff from here to the end of this section
It is also interesting to compare the approximate policy iteration
algorithm to the approximate linear programming algorithm we presented
in \secref{sec:ALP}. In the approximate linear programming algorithm,
we never need to compute the decision list policy. The policy is
always represented implicitly by the $Q_a$ functions. Thus, this
algorithm does not require explicit computation or manipulation of the
greedy policy. This difference has two important consequences: one
computational and the other in terms of generality.

First, not having to compute or consider the decision lists makes
approximate linear programming faster and easier to implement.  In
this algorithm, we generate a single LP with one cost network for each
action and never need to compute a decision list policy. On the other
hand, in each iteration, approximate policy iteration needs to
generate two LPs for every branch of the decision list of size $L$,
which is usually significantly longer than $\abs{A}$, with a total of
$2L$ cost networks.  In terms of representation, we do not require the
policies to be compact; thus, we do not need to make the default
action assumption.  Therefore, the approximate linear programming
algorithm can deal with a more general class of problems, where each
action can have its own independent DBN transition model.  On the
other hand, as described in \secref{sec: Approximate Policy
Iteration}, approximate policy iteration has stronger guarantees in
terms of error bounds. These differences will be highlighted further
in our experimental results presented in \secref{sec:Experimental
results}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computing Bounds on Policy Quality}\label{sec: Bellman error}

%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Factored Bellman error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\begin{pseudocode}
FactoredBellmanErr ($P$, $R$, $\discount$, $\Basis$, $\order$, $\bwhat$) \\ 
\tab \codecomment{$P$ is the factored transition model.}\\
\tab \codecomment{$R$ is the set of factored reward functions.}\\
\tab \codecomment{$\discount$ is the discount factor.}\\
\tab \codecomment{$\Basis$ is the set of basis functions $\Basis = \{\basis_1, \ldots, \basis_k\}$.}\\
\tab \codecomment{$\order$ stores the elimination order. } \\
\tab \codecomment{$\bwhat$ are the weights for the linear value function. } \\
\tab \codecomment{Return the Bellman error for the value function $\Basismatrix \bwhat$. }
\bcsone
\item \codecomment{Cache the backprojections of the basis functions.}
\item {\bf For} each basis function $\basis_i \in \Basis$; for each action $a$:
\bcstwo
\item {\bf Let} $\backbasis_i^a = \bprojfun{a}{\basis_i}$.
\ecstwo
\item \codecomment{Compute decision list policy for value function $\Basismatrix \bwhat$.}
\item {\bf Let} $\dlisthat = \mbox{DecisionListPolicy}(R^a + \discount \sum_i \what_i \backbasis_i^a )$.
\item \codecomment{Initialize indicators. }
\item {\bf Let} $\indicator = \{\}$.
\item \codecomment{Initialize Bellman error. }
\item {\bf Let} $\eps = 0$.
\item \codecomment{For every branch of the decision list policy, generate the relevant cost networks, solve it with variable elimination, and update the indicators to constraint the state space for future branches.}
\item {\bf For} each branch $\langle\event_j,a_j\rangle$ in the decision list policy $\dlisthat$:
\bcstwo
\item \codecomment{Instantiate the variables in $\Event_j$ to the assignment given in $\event_j$.}
\item {\bf Instantiate} the set of functions $\{\what_1 (\basis_1-\discount\backbasis_1^{a_j}),\ldots,\what_k (\basis_k-\discount\backbasis_k^{a_j})\}$ with the partial state assignment $\event_j$ and store in $C$.
\item {\bf Instantiate} the target functions $R^{a_j}$ with the partial state assignment $\event_j$ and store in $\vecb$.
\item {\bf Instantiate} the indicator functions $\indicator$ with the partial state assignment $\event_j$ and store in $\indicator'$.
\item \codecomment{Use variable elimination to solve first cost network, and update Bellman error, if error for this branch is larger.}
\item {\bf Let} $\eps = \max\left( \eps , \mbox{VariableElimination}(C - \vecb + \indicator', \order)\right)$. 
\item \codecomment{Use variable elimination to solve second cost network, and update Bellman error, if error for this branch is larger.}
\item {\bf Let} $\eps = \max\left( \eps , \mbox{VariableElimination}(- C + \vecb + \indicator', \order)\right)$. 
\item \codecomment{Update the indicator functions.}
\item {\bf Let} $\indicator_j(\bx) = - \infty \indicone(\bx = \event_j)$ and update the indicators $\indicator = \indicator \cup \indicator_j$. 
\ecstwo
\item {\bf Return} $\eps$.
\ecsone
\end{pseudocode}
\caption{Algorithm for computing Bellman error for factored value function $\Basismatrix\bwhat$. }\label{code:factored Bellman Error}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We have presented two algorithms for computing approximate
solutions to factored MDPs. All these algorithms generate linear
value functions which can be denoted by $\Basismatrix \bwhat$, where $\bwhat$
are the resulting basis function weights. In practice, the agent
will define its behavior by acting according to the greedy policy
$\policyhat = \greedy{\Basismatrix \bwhat}$. One issue that remains is how
this policy $\policyhat$ compares to the true optimal policy
$\policystar$; that is, how the \emph{actual} value
$\V_{\policyhat}$ of policy $\policyhat$ compares to $\Vstar$.

In \secref{sec:approximate MDP algs}, we showed some \emph{a priori} 
bounds for the quality of the policy. Another possible procedure is to
compute an \emph{a posteriori} bound. That is, given our resulting
weights $\bwhat$, 
%CGjairf2 rephrased
we compute a bound on the loss of acting
according to the greedy policy $\policyhat$ rather than the
optimal policy. This can be achieved by using the \emph{Bellman
error} analysis of Williams and
Baird~\citeyear{Williams+Baird:1993}.

%RP3: delete extra defined as
The \emph{Bellman error} is defined as  $\Bellman(\V) =
\maxnorm{\Tbell \V - \V}$. Given the greedy policy $\policyhat =
\greedy{\V}$, their analysis provides the bound: \beqn \label{eqn:
Bellman error} \maxnorm{\Vstar - \V_{\policyhat}} \leq
\frac{2\discount\Bellman(\V)}{1 - \discount}. \eeqn Thus, we can
use the Bellman error $\Bellman(\Basismatrix \bwhat)$ to evaluate the quality
of our resulting greedy policy.

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 state variables. Koller and Parr~\citeyear{Koller+Parr:2000} suggested
that structure in the factored MDP can be exploited to compute the
Bellman error efficiently. Here, we show how this error
bound can be computed by a set of cost networks using a similar
construction to the one in our max-norm projection algorithms. 
%RP4: Added
This technique can be used for any $\policyhat$ that can be
represented as a decision list and does not depend on the algorithm used to
determine the policy.  Thus, we can apply this technique to solutions
determined approximate linear programming if the action descriptions permit
a decision list representation of the policy.


%CGjairf2 redid next few paragraphs.
For some set of weights $\bwhat$, the Bellman error is given by:
\beqa
\Bellman(\Basismatrix \bwhat) & = &
\maxnorm{\Tbell \Basismatrix \bwhat - \Basismatrix \bwhat} ; \\ 
& = & \max\left( 
\begin{array}{l}
\max_{\bx} \sum_i w_i \basis_i(\bx) - R_{\policyhat}(\bx) - \discount \sum_{\bx'} P_{\policyhat}(\bx' \mid \bx) \sum_j w_j \basis_j(\bx') \ ,  \\
\max_{\bx} R_{\policyhat}(\bx) + \discount \sum_{\bx'} P_{\policyhat}(\bx' \mid \bx) \sum_j w_j \basis_j(\bx') - \sum_i w_i \basis_i(\bx)
\end{array}
\right).
\eeqa
If the rewards $R_{\policyhat}$ and the transition model $P_{\policyhat}$ are factored appropriately, then we can compute each one of these two maximizations ($\max_{\bx}$) using variable elimination in a cost network as described in \secref{sec: Maximizing over the state space}. However, $\policyhat$ is a decision list policy and it does not induce a factored transition model. Fortunately, as in the approximate policy iteration algorithm in \secref{sec:API}, we can exploit the structure in the decision list to perform such maximization efficiently. 
In particular, as in approximate policy iteration, we will generate two cost networks for each branch in the decision list. To guarantee that our maximization is performed only over states where this branch is relevant, we include the same type of indicator functions, which will force irrelevant states to have a value of $- \infty$, thus guaranteeing that at each point of the decision list policy we obtain the corresponding state with the maximum error. The state with the overall largest Bellman error will be the maximum over the ones generated for each point the in the decision list policy.
The complete factored algorithm for computing the Bellman error is outlined in \figref{code:factored Bellman Error}.

One last interesting note concerns our approximate policy
iteration algorithm with max-norm projection of \secref{sec:API}.
In all our experiments, this algorithm converged, so that $\bw\at{t}
= \bw\at{t+1}$ after some iterations. If such convergence occurs,
then the objective function $\phi\at{t+1}$ 
%CGjairf2 added
of the linear program in
our last iteration
is equal to the Bellman error of the final policy:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lemma: convergent policy iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{lem}\label{lem: convergent policy iteration}
If approximate policy iteration with max-norm projection
converges, so that $\bw\at{t} = \bw\at{t+1}$ for some iteration $t$,
then the max-norm projection error $\phi\at{t+1}$ of the last
iteration is equal to the Bellman error for the final value
function estimate $\Basismatrix\bwhat = \Basismatrix\bw\at{t}$: \bdm \Bellman(\Basismatrix\bwhat) =
\phi\at{t+1}.\edm
\end{lem}
\proof See \appenref{proof:lem: convergent policy iteration}. \qed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Thus, we can bound the loss of acting according to the final
policy $\policy\at{t+1}$ by substituting $\phi\at{t+1}$ into the
Bellman error bound:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corollary: convergent policy iteration bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{cor}\label{cor: convergent policy iteration bound}
If approximate policy iteration with max-norm projection converges
after $t$ iterations to a final value function estimate $\Basismatrix\bwhat$
associated with a greedy policy $\policyhat = \greedy{\Basismatrix \bwhat}$,
then the loss of acting according to $\policyhat$ instead of the
optimal policy $\policystar$ is bounded by: \bdm \maxnorm{\Vstar -
\V_{\policyhat}} \leq \frac{2\discount\phi\at{t+1}}{1 -
\discount}, \edm where $\V_{\policyhat}$ is the \emph{actual}
value of the policy $\policyhat$. \qed
\end{cor}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Therefore, when approximate policy iteration converges we can
obtain a bound on the quality of the resulting policy without
%RP1:
%needing to compute explicitly the Bellman error.
needing to compute the Bellman error explicitly.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Exploiting Context-specific Structure}\label{sec: CSI}
%CGjair1 added whole section

Thus far, we have presented a suite of algorithms which exploit additive
structure in the reward and basis functions and sparse connectivity in
the DBN representing the transition model. However, there exists another
important type of structure that should also be exploited for efficient
decision making: \emph{context-specific independence} (CSI). 
For example, consider an agent responsible for building and maintaining a
house, if the painting task can only be completed after the plumbing and
the electrical wiring have been installed, then the probability that the
painting is done is 0, in all contexts where plumbing or electricity are
not done, \emph{independently} of the agents action.  The representation
we have used so far in this paper would use a table to represent this
type of function.  This table is exponentially large in the number of
variables in the scope of the function, and ignores the context-specific
structure inherent in the problem definition.

Boutilier \emph{et
al.}~\cite{Boutilier+al:1995,Dearden+Boutilier:1997,Boutilier+al:1999,Boutilier+al:2000} have developed a set of algorithms which can exploit CSI in the
transition and reward models to perform efficient (approximate)
planning. Although this approach is often successful in problems where
the \emph{value function} contains sufficient context-specific structure,
the approach is not able to exploit the additive structure which is also
often present in real-world problems. 

In this section, we extend the factored MDP model to include
context-specific structure.  We present a simple, yet effective extension
of our algorithms which can exploit both CSI and additive structure to
obtain efficient approximations for factored MDPs.  We first extend the
factored MDP representation to include context-specific structure and
then show how the basic operations from \secref{sec: basic ops} required
by our algorithms can be performed efficiently in this new representation. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Factored MDPs with Context-specific and Additive Structure}\label{sec: CSI model}

There are several representations for context-specific functions.  The
most common are decision trees~\cite{Boutilier+al:1995}, algebraic
decision diagrams 
%CGjairf2 we had the wrong ref here
%(ADDs)~\cite{Boutilier+al:2000}, and
(ADDs)~\cite{Hoey+al:1999}, and
rules~\cite{Zhang+Poole:1999}.  We choose to use rules as our basic
representation, for two main reasons.  First, the rule-based
representation allows a fairly simple algorithm for variable elimination,
which is a key operation in our framework.  Second, rules are not
required to be mutually exclusive and exhaustive, a requirement that can
be restrictive if we want to exploit additive independence, where
functions can be represented as a linear combination of a
set of non-mutually exclusive functions.  

%CGjairf2 created a definition for consistent and changed the organization of this section a little.

We begin by describing the rule-based representation (along the lines of
Zhang and Poole's presentation~\citeyear{Zhang+Poole:1999}) for the
probabilistic transition model, in particular, the CPDs of our DBN model.  Roughly speaking, each rule
corresponds to some set of CPD entries that are all associated with a
particular probability value.
%CGjairf2 added
These entries with the same value are referred to as \emph{consistent} contexts:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Consistent
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:consistent}
Let $\Context \subseteq \{\bX,\bX'\} $ and $\context \in \Val(\Context)$.
We say that $\context$ is \emph{consistent} with $\bb \in \Val(\bB)$, for $\bB \subseteq \{\bX,\bX'\}$, if
$\context$ and $\bb$ have the same assignment for the variables in $\Context \cap \bB$. \qed
\end{defn}
%CGjairf2 added
The probability of these consistent contexts will be represented by \emph{probability rules}:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Probability rules
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}\label{def: probability rules}
A \emph{probability rule} $\prule = \valuerule{\context:\pvalue}$ is a
function $\prule : \{\bX,\bX'\} \mapsto [0,1]$, where the
\emph{context} $\context \in \Val(\Context)$ for $\Context \subseteq
\{\bX,\bX'\}$ and $\pvalue \in [0,1]$, such that $\prule(\bx,\bx')
= p$ if $(\bx,\bx')$ is consistent with $\context$ and is equal to $1$
otherwise. \qed
\end{defn}

In this case, it is convenient to require that the rules be mutually
exclusive and exhaustive, so that each CPD entry is uniquely defined by
its association with a single rule.  
\begin{defn}\label{def: rule CPD}
A \emph{rule-based conditional probability distribution}
(rule CPD) $P_a$ is a function $P_a: (\{X_i'\} \cup \bX) \mapsto [0,1]$,
composed of a set of probability rules
$\{\prule_1,\prule_2,\dots,\prule_m\}$ whose contexts are mutually
exclusive and exhaustive.  We define:
\bdm
P_a(x_i' \mid \bx) = \prule_j(\bx,\bx'),
\edm
where $\prule_j$ is the unique rule in $P_a$ for which $\context_j$ is
consistent with $(x_i',\bx)$.  We require that, for all $\bx$, 
\bdm
 \sum_{x_i'} P_a(x_i' \mid \bx) = 1. \ \qedmath
\edm
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We can define $\Parents_a(X_i')$ to be the union of the contexts of the
rules in $P_a(X_i' \mid \bX)$.  An example of a CPD represented by a set
of probability rules is shown in \figref{fig:painting rule cpd}.  


\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\hspace{0.2cm}\epsfig{figure=Figures/action-paint-rule-cpd.eps,height=5.2cm}&
\hspace{0.2cm}\epsfig{figure=Figures/action-not-paint-rule-cpd.eps,height=5.2cm}\\
\bf (a)  & \bf (b) 
\vspace{1cm} \\
\begin{tabular}{c}
$\prule_1=\valuerule{\neg\fmbox{Electrical}:0}$\\
$\prule_2=\valuerule{\fmbox{Electrical} \wedge \neg\fmbox{ Plumbing}:0}$\\
$\prule_3=\valuerule{\fmbox{Electrical}\wedge\fmbox{Plumbing}:0.95}$
\end{tabular}  & 
\begin{tabular}{c}
$\prule_4 = \valuerule{\neg\fmbox{Electrical}:0}$\\
$\prule_5 = \valuerule{\fmbox{Electrical} \wedge \neg\fmbox{Plumbing}:0}$\\
$\prule_6 = \valuerule{\fmbox{Electrical} \wedge \fmbox{Plumbing} \wedge
      \neg\fmbox{Painting}:0}$\\
$\prule_7 = \valuerule{\fmbox{Electrical} \wedge \fmbox{Plumbing}\wedge \fmbox{Painting}:0.9}$
\end{tabular}  \\
\bf (c)  & \bf (d)
\end{tabular}
\caption{Example CPDs for variable the $\mbox{\it Painting'\/} = \true$ represented as decision
trees: (a) when the action is paint; (b) when the action is not
paint. The same CPDs can be represented by probability rules as shown in
(c) and (d), respectively.} 
\label{fig:painting rule cpd}
\end{center}
\vspace*{-.5cm}
\end{figure}

Rules can also be used to represent additive functions, such as reward or
basis functions.  We represent such context specific value dependencies
using \emph{value rules}:  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Value rules
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:value rules}
A \emph{value rule} $\rrule = \valuerule{\context:\vvalue}$ is a function
$\rrule: \bX \mapsto \Real$ such that
$\rrule(\bx) = v$ when $\bx$ is consistent with $\context$ and $0$
 otherwise. \qed
\end{defn}
Note that a value rule $\valuerule{\context:\vvalue}$ has a scope
$\Context$.  

It is important to note that value rules are not required to be mutually
exclusive and exhaustive.  Each value rule represents a (weighted)
indicator function, which takes on a value $v$ in states consistent with
some context $\context$, and 0 in all other states.  In any given state,
the values of the zero or more rules consistent with that state are
simply added together.  
\begin{example}
In our construction example, we might have a set of rules: 
\bdm
\begin{array}{c}
\rrule_1 = \valuerule{\mbox{Plumbing} = \mbox{\emph{done}}:100}; \\
\rrule_2 = \valuerule{\mbox{Electricity} = \mbox{\emph{done}}:100}; \\
\rrule_3 = \valuerule{\mbox{Painting} = \mbox{\emph{done}}:100}; \\
\rrule_4 = \valuerule{\mbox{Action} = \mbox{\emph{plumb}}:-10}; \\
\vdots
\end{array}
\edm
which, when summed together, define the reward function $R = \rrule_1 +
\rrule_2 + \rrule_3 + \rrule_4 + \cdots$. \qed
\end{example}

In general, our reward function $R^a$ is represented as a
\emph{rule-based function}:  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: rule function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}\label{def: rule function}
A \emph{rule-based function} $f: \bX \mapsto
\Real$ is composed of a set of rules $\{\rrule_1,\dots,\rrule_\numrules\}$
such that $f(\bx) = \sum_{i=1}^{\numrules} \rrule_i(\bx)$. \ \qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the same manner, each one of our basis functions $\basis_j$ is now
represented as a rule-based function. 

This notion of a rule-based function is related to the tree-structure
functions used by Boutilier \etal~\citeyear{Boutilier+al:2000},
but is substantially more general.  In the tree-structure value
functions, the rules corresponding to the different leaves are mutually
exclusive and exhaustive.  Thus, the total number of different values
represented in the tree is equal to the number of leaves (or rules).  In
the rule-based function representation, the rules are not mutually
exclusive, and their values are added to form the overall function value
for different settings of the variables.  Different rules are added in
different settings, and, in fact, with $k$ rules, one can easily generate
$2^k$ different possible values, as is demonstrated in
\secref{sec:Experimental results}.  Thus, the rule-based functions can
provide a compact representation for a much richer class of value
functions. 

Using this rule-based representation, we can exploit both CSI and
additive independence in the representation of our factored MDP and basis
functions. We now show how the basic operations in \secref{sec: basic
ops} can be adapted to exploit our rule-based representation. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CGjairf2 significantly changed this section
\subsection{Adding, Multiplying and Maximizing Consistent Rules}\label{sec: CSI rule merge}

\newcommand{\split}[2]{Split(#1 \angle #2)}

In our table-based algorithms, we relied on standard sum and product
operators applied to tables. In order to exploit CSI using a
rule-based representation, we must redefine these standard
operations. In particular, the algorithms will need to add or multiply
rules 
%DK4 - this was confusing
%that give different values to the same set of states.
%, \ie, consistent rules. 
that ascribe values to overlapping sets of states.

We will start by defining these operations for rules with the same context:
%%%%%%%%%%%%%%%%%%%
% Defn. rule product
%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:rule product}
%RPfinal:  shortened
Let $\rrule_1 = \valuerule{\context:\vvalue_1}$ and
$\rrule_2 = \valuerule{\context:\vvalue_2}$  be two rules with 
context $\context$. Define the \emph{rule product} as $\rrule_1 \times
\rrule_2 = \valuerule{\context:\vvalue_1 \cdot \vvalue_2}$, and the
\emph{rule sum} as $\rrule_1 + \rrule_2 = \valuerule{\context:\vvalue_1 +
\vvalue_2}$. \qed 
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Note that this definition is restricted to rules with the same context. We will address this issue in a moment. First, we will introduce an additional operation which maximizes a variable from a set of rules, which otherwise share a common context:
%%%%%%%%%%%%%%%%%%%
% Defn. rule maximization
%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:rule max}
Let $Y$ be a variable with $\Val[Y] = \{y_1, \dots, y_k\}$, and let $\rrule_i$, for each $i = 1, \dots, k$, be a rule of the form
$\rrule_i = \valuerule{\context \wedge Y = y_i:\vvalue_i}$. Then for the
rule-based function $f = \rrule_1 + \cdots + \rrule_k$, define the
\emph{rule maximization} over $Y$ as $\max_{Y} f = \valuerule{\context:
\max_i \vvalue_i}$ . \qed 
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
After this operation, $Y$ has been maximized out from the scope of the function $f$.

These three operations we have just described can only be applied 
%DK4
%very restricted settings, \ie, 
to sets of rules that satisfy very stringent conditions. To
make our set of rules amenable to the application of these operations,
we might need to refine some of these rules. We therefore define the
following operation: 
%%%%%%%%%%%%%%%%%%%
% Defn. rule split
%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:rule split}
Let $\rrule = \valuerule{\context:\vvalue}$ be a rule, and $Y$ be a variable. Define the \emph{rule split} $\split{\rrule}{Y}$ of $\rrule$ on a variable $Y$ as follows: If $Y \in \Scope[\Context]$, then $\split{\rrule}{Y} =  \{\rrule\}$; otherwise,
\bdm
\split{\rrule}{Y} = \{\valuerule{\context\wedge Y= y_i :\vvalue} \mid y_i
\in \Val[Y]\} \ . \qedmath 
\edm
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Thus, if we split a rule $\rrule$ on variable $Y$ that is not in the scope of the context of $\rrule$, then we generate a new set of rules, with one for each assignment in the domain of $Y$. 

%RP5: make -> extend
In general, the purpose of rule splitting is to extend the context $\context$
of one rule $\rrule$ coincide with the context $\context'$ of another
consistent rule $\rrule'$.
Naively, we might take all variables in
$\Scope[\Context'] - \Scope[\Context]$ and split $\rrule$ recursively on
each one of them. However, this process creates unnecessarily many rules:
If $Y$ is a variable in $\Scope[\Context'] - \Scope[\Context]$ and we split
$\rrule$ on $Y$, then only one of the $\abs{\Val[Y]}$ new rules generated
will remain consistent with $\rrule'$: the one which has the same
assignment for $Y$ as the one in $\context'$. Thus, only this consistent
rule needs to be split further. We can now define the recursive splitting
procedure that achieves this more parsimonious representation: 
%%%%%%%%%%%%%%%%%%%%%%%
% defn recursive rule splitting
%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn} \label{dfn:recursive splitting}
Let $\rrule = \valuerule{\context:\vvalue}$ be a rule, and $\bb$ be a context
such that $\bb \in \Val[\bB]$. Define the \emph{recursive rule split}
$\split{\rrule}{\bb}$ of $\rrule$ on a context $\bb$ as follows:   
\begin{enumerate}
\item $\{\rrule\}$, if $\context$ is not consistent with $\bb$; else,
\item $\{\rrule\}$, if $\Scope[\bB] \subseteq \Scope[\Context]$; else,
\item $\{\split{\rrule_i}{\bb} \mid \rrule_i \in \split{\rrule}{Y}\}$, for some variable $Y \in \Scope[\bB] - \Scope[\Context]$ . 
\end{enumerate} \qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
In this definition, each variable $Y \in \Scope[\bB] - \Scope[\Context]$ leads to the
generation of $k = |\Val(Y)|$ rules at the step in which it is split.  
However, only one of these $k$ rules is used in the next recursive
step because only one is consistent with $\bb$.
%RP5: added because...
Therefore, the size of the split set is simply $1 + \sum _{Y \in \Scope[\bB] - \Scope[\Context]} (|Dom(Y)| - 1)$.  This size is independent of the order
in which the variables are split within the operation.

Note that only one of the rules in $\split{\rrule}{\bb}$ is consistent
with $\bb$: the one with context $\context\wedge\bb$.
Thus, if we want to add two consistent rules
$\rrule_1 = \valuerule{\context_1:\vvalue_1}$ and
$\rrule_2 = \valuerule{\context_2:\vvalue_2}$, then all we need to 
do is replace these rules by the set: 
\[
\split{\rrule_1}{\context_2} \cup \split{\rrule_2}{\context_1} ,\]
and then simply replace the resulting rules $\valuerule{\context_1\wedge\context_2:\vvalue_1}$ 
and $\valuerule{\context_2\wedge\context_1:\vvalue_2}$ by their sum
$\valuerule{\context_1\wedge\context_2:\vvalue_1+\vvalue_2}$.
%RP5:
Multiplication is performed in an analogous manner.

\begin{example}
Consider adding the following set of consistent rules: 
\[ \begin{array}{l}
\rrule_1 = \valuerule{a \wedge b: 5}, \\
\rrule_2 = \valuerule{a \wedge \neg c \wedge d: 3}.
\end{array} \]
\noindent
In these rules, the context $\context_1$ of $\rrule_1$ is $a \wedge b$, and the
context $\context_2$ of $\rrule_2$ is $a \wedge \neg c \wedge d$. 

Rules $\rrule_1$ and $\rrule_2$ are consistent, therefore, we must split them to perform the addition operation: 
\[ \split{\rrule_1}{\context_2} = \left \{ \begin{array}{l}
\valuerule{a \wedge b \wedge c: 5}, \\
\valuerule{a \wedge b \wedge \neg c \wedge \neg d: 5}, \\
\valuerule{a \wedge b \wedge \neg c \wedge d: 5}.
\end{array}
\right .
\]
Likewise, 
\[ \split{\rrule_2}{\context_1} =  \left\{\begin{array}{l}
\valuerule{a \wedge \neg b \wedge \neg c \wedge d: 3}, \\
\valuerule{a \wedge b \wedge \neg c \wedge d: 3}. 
\end{array} \right.
\]
\noindent
The result of adding rules $\rrule_1$ and $\rrule_2$ is
\[ \begin{array}{l}
\valuerule{a \wedge b \wedge c: 5}, \\
\valuerule{a \wedge b \wedge \neg c \wedge \neg d: 5},\\ 
\valuerule{a \wedge b \wedge \neg c \wedge d: 8}, \\ 
\valuerule{a \wedge \neg b \wedge \neg c \wedge d: 3}. \
\end{array}
\]
\qedmath
\end{example}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Rule-based One-step Lookahead}\label{sec: CSI lookahead}

Using this compact rule-based representation, we are able to compute a one-step lookahead plan efficiently for models with significant context-specific or additive independence.

As in \secref{sec: one-step lookahead} for the table-based case, the rule-based $Q_a$ function can be represented as the sum of the reward function and the discounted expected value of the next state. Due to our linear approximation of the value function, the expectation term is, in turn, represented as the linear combination of the backprojections of our basis functions.
To exploit CSI, we are representing the rewards and basis functions as rule-based functions. To represent $Q_a$ as a rule-based function, it is sufficient for us to show how to represent the backprojection $\backbasis_j$ of the basis function $\basis_j$ as a rule-based function.


Each $\basis_j$ is a rule-based function, which can be written as $\basis_j(\bx)
= \sum_i \rrule_i\at{\basis_j}(\bx)$, where  $\rrule_i\at{\basis_j}$ has
%the form $\valuerule{\rrule_i\at{\basis_j};\context_i\at{\basis_j} :
the form $\valuerule{\context_i\at{\basis_j} :
\vvalue_i\at{\basis_j}}$.  Each rule is a restricted scope function; thus,
we can simplify the backprojection as:
\beqa
\backbasis_j^a(\bx) & = & \sum_{\bx'} P_a(\bx' \mid \bx) \basis_j(\bx') \ ;\\
& = & \sum_{\bx'} P_a(\bx' \mid \bx)\sum_i \rrule_i\at{\basis_j}(\bx'); \\
& = & \sum_i \sum_{\bx'} P_a(\bx' \mid \bx) \rrule_i\at{\basis_j}(\bx'); \\
& = & \sum_i \vvalue_i\at{\basis_j} P_a(\context_i\at{\basis_j}
\mid \bx); 
\eeqa 
where the term $\vvalue_i\at{\basis_j}
P_a(\context_i\at{\basis_j}  \mid \bx)$ can be written as a rule
function. We denote this backprojection operation by
$\rulebproj{a}{\rrule_i\at{\basis_j}}$. 

The backprojection procedure, described in \figref{code:rulebproj},
follows three steps.  First, the relevant rules are selected: In the
CPDs for the variables that appear in the context of $\rrule$, we select
the rules consistent with this context, as these are the only rules that
play a role in the backprojection computation. Second, we multiply all
consistent probability rules to form a local set of mutually-exclusive
rules. This procedure is analogous to the addition procedure described in
\secref{sec: CSI rule merge}. Now that we have represented the
probabilities that can affect $\rrule$ by a mutually-exclusive set, we
can simply represent the backprojection of $\rrule$ by the product of
these probabilities with the value of $\rrule$. That is, the
backprojection of $\rrule$ is a rule-based function with one rule for
each one of the mutually-exclusive probability rules $\prule_i$. The
context of this new value rule is the same as that of $\prule_i$, 
and the value is the product of the probability of $\prule_i$ and
the value of $\rrule$.  

%CGjairf2 removed residual from code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code: Rule-based backprojection 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{pseudocode}
$\rulebproj{a}{\rrule}$ , \ \ where $\rrule$ is given by $\valuerule{\context:\vvalue}$, with $\context \in \Val[\Context]$.
\bcsone
\item {\bf Let} $g = \{\}$.
\item {\bf Select} the set $\setP$ of relevant probability rules:\\
$\setP =$ $\{\mbox{$\prule_j \in P(X_i'\mid\Parents(X_i'))$} \ |\  X_i' \in \Context$ and $\context$ is consistent with $\context_j\}$.
\item {\bf Remove} the $\bX'$ assignments from the context of all rules in $\setP$.
\item \emph{// Multiply consistent rules:}
\item {\bf While} there are two consistent rules
$\prule_1 = \valuerule{\context_1:\pvalue_1}$ and
$\prule_2 = \valuerule{\context_2:\pvalue_2}$: 
    \bcstwo
    \item {\bf If} $\context_1 = \context_2$, replace these two rules by $\valuerule{\context_1:\pvalue_1\pvalue_2}$;
    \item {\bf Else} replace these two rules by the set: $\split{\prule_1}{\context_2} \cup \split{\prule_2}{\context_1}$. 
    \ecstwo
\item \emph{// Generate value rules:}
\item {\bf For} each rule $\prule_i$ in $\setP$:
    \bcstwo
    \item {\bf Update} the backprojection $g = g \cup \{\valuerule{\context_i:\pvalue_i \vvalue} \}.$
    \ecstwo
\item {\bf Return} $g$.
\ecsone
\end{pseudocode}
\caption{Rule-based backprojection. }\label{code:rulebproj}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{example}
For example, consider the backprojection of a simple rule,
$$\rrule = \valuerule{\mbox{ Painting = \emph{done}}:100},$$ through
the CPD in \figref{fig:painting rule cpd}(c) for the paint action: \beqa
\rulebproj{\fmbox{\emph{paint}}}{\rrule} & = & \sum_{\bx'} P_{\fmbox{\emph{paint}}}(\bx' \mid \bx) \rrule(\bx') ;\\
& = & \sum_{\mbox{\tiny Painting}'} P_{\fmbox{\emph{paint}}}(\fmbox{Painting}' \mid
\bx)
\rrule(\fmbox{Painting}') ; \\
& = & 100 \ \prod_{i=1}^3 \prule_i(\mbox{Painting'} = done,\bx)\ .\eeqa
Note that the product of these simple rules is equivalent to the decision tree CPD shown in
\figref{fig:painting rule cpd}(a).  Hence, this product is equal to $0$ in
most contexts, for example, when electricity is not done at time $t$. The product
in non-zero only in one context: in the context associated with rule
$\prule_3$. Thus, we can express the
result of the backprojection operation by a rule-based function with a single rule:  
\bdm
 \rulebproj{\fmbox{\emph{paint}}}{\rrule} = \valuerule{\fmbox{Plumbing}\wedge\fmbox{Electrical} :95}. \edm
Similarly, the backprojection of $\rrule$ when the action is not paint can also be represented by a single rule:
\bdm
 \rulebproj{\neg\fmbox{\emph{paint}}}{\rrule} =
\valuerule{\fmbox{Plumbing}\wedge\fmbox{Electrical}
\wedge\fmbox{Painting} :90}. \ \qedmath
\edm
\end{example}

Using this algorithm, we can now write the \emph{backprojection} of the rule-based basis function $\basis_j$ as:
\begin{equation}
\backbasis_j^a (\bx) = \sum_i \rulebproj{a}{\rrule_i\at{\basis_j}},
\end{equation}
where $\backbasis_j^a$ is a sum of rule-based functions, and therefore also
a rule-based function.
For simplicity of notation, we use $\backbasis_j^a
=\rulebproj{a}{\basis_j}$ to refer to this definition of backprojection. 
Using this notation, we can write $Q_a(\bx) = R^a(\bx) +
\discount\sum_j w_j \backbasis_j^a(\bx)$, which is again a rule-based function.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Rule-based Maximization Over the State Space}\label{sec: CSI VE}

The second key operation required to extend our planning algorithms to exploit CSI is to modify the variable elimination algorithm in \secref{sec: Maximizing over the state space} to handle the rule-based representation.
In \secref{sec: Maximizing over the state space}, we showed that the maximization of a linear combination of table-based functions with restricted scope can be performed efficiently using non-serial
dynamic programming~\cite{Bertele+Brioschi:1972}, or variable
elimination. To exploit structure in rules, we use an
algorithm similar to variable elimination in a Bayesian network with
context-specific independence~\cite{Zhang+Poole:1999}.

Intuitively, the algorithm operates by
selecting the value rules relevant to the variable being maximized in the current iteration. Then, a local maximization is performed over this subset of the rules, generating a new set of rules without the current variable. The procedure is then repeated recursively until all variables have been eliminated. 

More precisely, our algorithm ``eliminates'' variables one by one, where
the elimination process performs a maximization step over the variable's
domain.  Suppose that we are eliminating $X_i$, whose collected value rules
lead to a rule function $f$, and $f$ involves additional variables in some
set $\bB$, so that $f$'s scope is $\bB \cup \{X_i\}$.  We need to compute
the maximum value for $X_i$ for each choice of $\bb \in \Val[\bB]$.  We use
$\maxout{f,X_i}$ to denote a procedure that takes a rule function
$f(\bB,X_i)$ and returns a rule function $g(\bB)$ such that: $g(\bb) =
\max_{x_i} f(\bb,x_i)$.  Such a procedure is an extension of the variable
elimination algorithm of Zhang and Poole~\cite{Zhang+Poole:1999}.

%CGjairf2 removed residual from code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code: Rule-based maxout
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{pseudocode}
$\maxout{f,B}$
\bcsone
\item {\bf Let} $g = \{\}$.
\item {\bf Add} completing rules to $f$: $\valuerule{B=b_i:0}, i = 1,\dots,k$.
\item \emph{// Summing consistent rules:}
\item {\bf While} there are two consistent rules $\rrule_1 = \valuerule{\context_1:\vvalue_1}$ and $\rrule_2 = \valuerule{\context_2:\vvalue_2}$:
    \bcstwo
    \item {\bf If} $\context_1 = \context_2$, then replace these two rules by $\valuerule{\context_1:\vvalue_1 + \vvalue_2}$;
    \item {\bf Else} replace these two rules by the set: $\split{\rrule_1}{\context_2} \cup \split{\rrule_2}{\context_1}$. 
    \ecstwo
\item \emph{// Maximizing out variable $B$:}
\item {\bf Repeat} until $f$ is empty:
    \bcstwo
    \item {\bf If} there are rules $\valuerule{\context \wedge
          B=b_i:\vvalue_i}, \forall b_i \in \Val(B)$ : 
      \bcsthree
\item {\bf Then} remove these rules from $f$ and add rule $\valuerule{\context:\max_i\vvalue_i}$ to $g$;
\ecsthree
    \item {\bf Else} select two rules: $\rrule_i = \valuerule{\context_i \wedge B=b_i
:\vvalue_i}$ and $\rrule_j = \valuerule{\context_j \wedge B=b_j
:\vvalue_j}$ such that $\context_i$ is consistent with
$\context_j$, but not identical, and
%DK2 - formalized
%split $\rrule_i$ on $\context_j$ and $\rrule_j$ on $\context_i$
replace them with $\split{\rrule_i}{\context_j} \cup
    \split{\rrule_j}{\context_i}$ .
    \ecstwo
\item {\bf Return} $g$.
\ecsone
\end{pseudocode}
\caption{Maximizing out variable $B$ from rule function $f$. }\label{code:maxout}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The rule-based variable elimination algorithm maintains a set
$\F$ of value rules,  initially containing the set of rules to be maximized. The algorithm
then repeats the following steps for each variable $X_i$ until all variables have been eliminated: 
\benum
\item Collect all rules which depend on $X_i$ into $f_i$ ---
$f_i = \{\valuerule{\context:\vvalue} \in \F \mid  X_i \in
\Context\}$ --- and remove these rules from $\F$. 
\item Perform the local maximization step over $X_i$:  $g_i = \maxout{f_i,X_i}$;
\item Add the rules in $g_i$ to $\F$; now, $X_i$ has been ``eliminated.'' 
\eenum

The cost of this algorithm is polynomial in the number of new rules
generated in the maximization operation $\maxout{f_i,X_i}$. The number of
rules is never larger and in many cases exponentially smaller than the
complexity bounds on the table-based maximization in \secref{sec: Maximizing over the state space}, which, in
turn, was exponential only in the \emph{induced width} of the
cost network graph~\cite{Dechter:1999}.  However, the computational costs involved in
managing sets of rules usually imply that the computational advantage of
the rule-based approach over the table-based one will only be significant in problems that possess a
fair amount of context-specific structure.

In the remainder of this section, we present the algorithm for
computing the local maximization $\maxout{f_i,X_i}$. In the next section, we show how these ideas can be applied to extending the algorithm in \secref{sec: Factored LP} to exploit CSI in the LP representation for planning in factored MDPs.

The procedure, presented in \figref{code:maxout}, is divided into two
parts: first, all consistent rules are added together as described in
\secref{sec: CSI rule merge}; then, variable $B$ is maximized. This
maximization is performed by generating a set of rules, one for each
assignment of $B$, whose contexts have the same assignment for all
variables except for $B$, as in \defref{dfn:rule max}. This set is then substituted by a single rule
without a $B$ assignment in its context and with value equal to the
maximum of the values of the rules in the original set. 
Note that, to simplify the algorithm, we initially need to
add a set of value rules with $0$ value, which guarantee that our
rule function $f$ is complete (\ie, there is at least one rule consistent with
every context).   

%RP5:
The correctness of this procedure follows directly from the correctness of
the rule-based variable elimination procedure described by Zhang and Poole,
merely by replacing summations with product with $\max$, and products with
products with sums.  We conclude this section with a small example to
illustrate the algorithm:
%RP5:  Added some explanations to this
\begin{example}
Suppose we are maximizing $a$ for the following set of rules: 
\[ \begin{array}{l}
\rrule_1 = \valuerule{\neg a: 1}, \\
\rrule_2 = \valuerule{a \wedge \neg b: 2}, \\
\rrule_3 = \valuerule{a \wedge b \wedge \neg c: 3}, \\
\rrule_4 = \valuerule{\neg a \wedge b: 1}.
\end{array}
\]
When we add completing rules, we get:
\[ \begin{array}{l}
\rrule_5 = \valuerule{\neg a: 0}, \\
\rrule_6 = \valuerule{a : 0}.
\end{array}
\]
In the first part of the algorithm, we need to add consistent rules:
We add $\rrule_5$ to $\rrule_1$ (which remains unchanged), combine $\rrule_1$
with $\rrule_4$, $\rrule_6$ with $\rrule_2$, and then the split
of $\rrule_6$ on the context of $\rrule_3$, to get the following 
inconsistent set of rules: 
\[ \begin{array}{ll}
\rrule_2 = \valuerule{a \wedge \neg b: 2}, \\
\rrule_3 = \valuerule{a \wedge b \wedge \neg c: 3}, \\
\rrule_7 = \valuerule{\neg a \wedge b: 2}, &
\mbox{\footnotesize (from adding $\rrule_4$ to the consistent rule from $\split{\rrule_1}{\bb}$)} \\
\rrule_8 = \valuerule{\neg a \wedge \neg b: 1}, & \mbox{\footnotesize (from
$\split{\rrule_1}{\bb}$)} \\
\rrule_9 = \valuerule{a \wedge b \wedge c: 0}, & \mbox{\footnotesize (from
$\split{\rrule_6}{a \wedge b \wedge \neg c}$)}. \\
\end{array}
\]
Note that several rules with value $0$ are also generated, but not shown
here because they are added to other rules with consistent contexts.
We can move to the second stage (repeat loop) of \textit{MaxOut}. We remove $\rrule_2$, and
$\rrule_8$, and maximize $a$ out of them, to give:
\[ \begin{array}{l}
\rrule_{10} = \valuerule{\neg b: 2}. \\
\end{array} \]
We then select rules $\rrule_3$ and $\rrule_7$ and split $\rrule_7$ on $c$
($\rrule_3$ is split on the empty set and is not changed), 
\[ \begin{array}{l}
\rrule_{11} = \valuerule{\neg a \wedge b \wedge c: 2}, \\
\rrule_{12} = \valuerule{\neg a \wedge b \wedge \neg c: 2}.
\end{array} \]
Maximizing out $a$ from rules $\rrule_{12}$ and $\rrule_3$, we get:
\[ \begin{array}{l}
\rrule_{13} = \valuerule{b \wedge \neg c: 3}. \\
\end{array} \]
We are left with $\rrule_{11}$, which maximized over its counterpart $\rrule_9$
gives
\[ \begin{array}{l}
\rrule_{12} = \valuerule{ b \wedge \neg c: 2}. \\
\end{array} \]
Notice that, throughout this maximization, we have not split on the
variable $C$ when $\neg b \in \context_i$, giving us only 6 distinct
rules in the final result. This is not possible
in a table-based representation, since our functions would then be over
the 3 variables $a$,$b$,$c$, and therefore must have 8 entries.  \qed
\end{example}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Rule-based Factored LP}\label{sec: CSI factored LP}


In \secref{sec: Factored LP}, we showed that the LPs used in our algorithms have exponentially many constraints of the form: 
$\phi \geq \sum_i w_i \
c_i(\bx) - b(\bx), \forall \bx$, which can be substituted by a single,
equivalent,  non-linear constraint: $\phi \geq \max_{\bx}\sum_i w_i \
c_i(\bx) - b(\bx)$. We then showed that, using variable elimination, we can
represent this non-linear constraint by an equivalent set of linear
constraints in a construction we called the factored LP. The number of
constraints in the factored LP is linear in the size of the largest table
generated in the variable elimination procedure. This table-based algorithm
can only exploit additive independence.  We now extend the algorithm in
\secref{sec: Factored LP} to exploit \emph{both} additive and
context-specific structure, by using the rule-based variable elimination
described in the previous section.  

Suppose we wish to enforce the more general constraint $0 \geq \max_{\by}
F^{\bw}(\by)$, where $F^{\bw}(\by) = \sum_j f^{\bw}_j(\by)$ such that each
$f_j$ is a rule.  As in the table-based version, the superscript $\bw$
means that $f_j$ might depend on $\bw$.  Specifically, if $f_j$ comes
from basis function $\basis_i$, it is multiplied by the weight $w_i$; if
$f_j$ is a rule from the reward function, it is not.

In our rule-based factored linear program, we generate LP variables associated with contexts; we call these \emph{LP rules}.  An LP
rule has the form $\valuerule{\context:\lpvalue}$; it is associated
with a context $\context$ and a variable $\lpvalue$ in the linear program.
We begin by transforming all our original rules $f^{\bw}_j$ into LP rules
as follows:  If rule $f_j$ has the form $\valuerule{\context_j:\vvalue_j}$
and comes from basis function $\basis_i$, we introduce an LP rule
$\lprule_j = \valuerule{\context_j:\lpvalue_j}$ and the equality constraint
$\lpvalue_j = w_i \vvalue_j$. If $f_j$ has the same form but comes from a
reward function, we introduce an LP rule of the same form, but the equality
constraint becomes $\lpvalue_j = \vvalue_j$.

Now, we have only LP rules and need to represent the constraint:
$0 \geq max_{\by} \sum_j \lprule_j(\by)$. To represent such a constraint, we
follow an algorithm very similar to the variable elimination procedure in
\secref{sec: CSI VE}. The main difference occurs in the
$\maxout{f,B}$ operation in \figref{code:maxout}. Instead of generating new
value rules, we generate new LP rules, with associated new variables and
new constraints. The simplest case occurs when computing a split or
adding two LP rules. For example, when we add two value rules in the
original algorithm, we instead perform the following operation on their
associated LP rules:  If the LP rules are
$\valuerule{\context:\lpvalue_i}$  and
$\valuerule{\context:\lpvalue_j}$, we replace these by a new rule
$\valuerule{\context:\lpvalue_k}$, associated with a new LP
variable $\lpvalue_k$ 
with context $\context$, whose value should be
$\lpvalue_i + \lpvalue_j$.  To enforce this value constraint, we simply
add an additional constraint to the LP: $\lpvalue_k
= \lpvalue_i + \lpvalue_j$. A similar procedure can be followed when
computing the split.

More interesting constraints are generated when we perform a
maximization. In the rule-based variable elimination algorithm in \figref{code:maxout}, this maximization occurs
when we replace a set of rules:
\bdm
\valuerule{\context \wedge B=b_i:\vvalue_i}, \forall b_i \in \Val(B),
\edm           
by a  new rule 
\bdm \valuerule{\context:\max_i\vvalue_i}. 
\edm Following
the same process as in the LP rule summation above, if we are maximizing
\bdm
\lprule_i = \valuerule{\context \wedge B=b_i:\lpvalue_i}, \forall b_i \in \Val(B),
\edm           
we generate a new LP variable $\lpvalue_k$ associated with
the rule $\lprule_k = \valuerule{\context:\lpvalue_k}$.  
However, we cannot
add the nonlinear constraint $\lpvalue_k = \max_i \lpvalue_i$, but we can
add a set of equivalent linear constraints 
\bdm \lpvalue_k \geq \lpvalue_i, \ \forall i .
\edm


Therefore, using these simple operations, we can exploit structure in the rule functions to represent the nonlinear constraint $\lprule_n \geq max_{\by} \sum_j \lprule_j(\by)$, where $\lprule_n$ is the very last LP rule we generate. A final constraint $\lpvalue_n = \phi$ implies that we are representing exactly the constraints in \eqnref{eq: max
constraint form}, without having to enumerate every state. 

The correctness of our rule-based factored LP construction is a corollary of \thmref{thm: factored LP} and of the correctness of the rule-based variable elimination algorithm~\cite{Zhang+Poole:1999} .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corollary: rule-based factored LP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{cor}\label{cor: rule-based factored LP}
The constraints generated by the \emph{rule-based} factored LP construction are
equivalent to the non-linear constraint in \eqnref{eq: max
constraint form}. That is, an assignment to $(\phi,\bw)$ satisfies
the rule-based factored LP constraints if and only if it satisfies the
constraint in \eqnref{eq: max constraint form}. \qed
\end{cor}
The number of variables and constraints in the rule-based factored LP is
linear in the number of rules generated by the variable elimination
process. In turn, the number of rules is no larger, and often
exponentially smaller, than the number of entries in the table-based
approach. 

To illustrate the generation of LP constraints as just described, we now
present a small example:  
\begin{example}
Let $\lprule_1$, $\lprule_2$, $\lprule_3$, and $\lprule_4$ be the set of
LP rules which depend on the variable $b$ being maximized. Here, rule
$\lprule_i$ is associated with the LP variable $u_i$:
\[ \begin{array}{l}
\lprule_1 = \valuerule{a \wedge b:u_1}, \\
\lprule_2 = \valuerule{a \wedge b \wedge c: u_2}, \\
\lprule_3 = \valuerule{a \wedge \neg b: u_3}, \\
\lprule_4 = \valuerule{a \wedge b \wedge \neg c: u_4}. 
\end{array}
\] 
\noindent
In this set, note that rules $\lprule_1$ and $\lprule_2$ are
consistent. We combine them to generate the following rules: 
\[ \begin{array}{l}
\lprule_5 = \valuerule{a \wedge b \wedge c: u_5}, \\
\lprule_6 = \valuerule{a \wedge b \wedge \neg c: u_1}.  
\end{array}
\]
\noindent
and the constraint $u_1 + u_2 = u_5$. 
Similarly, $\lprule_6$ and $\lprule_4$ may be combined, resulting in: 
\[ \begin{array}{l}
\lprule_7 = \valuerule{a \wedge b \wedge \neg c: u_6}. 
\end{array}
\]
\noindent
with the constraint $u_6 = u_1 + u_4$. 
Now, we have the following three inconsistent rules for the maximization: 
\[ \begin{array}{l}
\lprule_3 = \valuerule{a \wedge \neg b: u_3}, \\
\lprule_5 = \valuerule{a \wedge b \wedge c: u_5}, \\
\lprule_7 = \valuerule{a \wedge b \wedge \neg c: u_6}.
\end{array}
\]
\noindent
Following the maximization procedure, since no pair of rules can be
eliminated right away, we split $\lprule_3$ and $\lprule_5$ to generate
the following rules:  
\[ \begin{array}{l}
\lprule_8 = \valuerule{a \wedge \neg b \wedge c: u_3}, \\
\lprule_9 = \valuerule{a \wedge \neg b \wedge \neg c: u_3}, \\
\lprule_5 = \valuerule{a \wedge b \wedge c: u_5}.
\end{array}
\]
\noindent
We can now maximize $b$ out from $\lprule_8$ and $\lprule_5$, resulting in the following rule and constraints respectively:
\[ \begin{array}{l}
\lprule_{10} = \valuerule{a \wedge c: u_7}, \\
u_7 \geq u_5, \\
u_7 \geq u_3. 
\end{array}
\]
\noindent
Likewise, maximizing $b$ out from $\lprule_9$ and $\lprule_6$, we get:
\[ \begin{array}{l}
\lprule_{11} = \valuerule{a \wedge \neg c: u_8}, \\
u_8 \geq u_3, \\
u_8 \geq u_6;
\end{array}
\]
which completes the elimination of variable $b$ in our rule-based factored LP. \qed
\end{example}

We have presented an algorithm for exploiting both additive and
context-specific structure in the LP construction steps of our planning
algorithms. This rule-based factored LP approach can now be applied
%DK4 - reversed order
directly in our approximate linear programming and
approximate policy iteration algorithms, which were presented in
Sections~\ref{sec:ALP} and~\ref{sec:API}.

The only additional modification required concerns the manipulation of
the decision list policies presented in \secref{sec: Computing greedy
policies}. Although approximate linear programming does not require any
explicit policy representation (or the default action model), approximate
policy iteration require us to represent such policy. Fortunately, no
major modifications are required in the rule-based case. In particular,
the conditionals $\langle\event_i,a_i,\delta_i\rangle$ in the decision
list policies are already context-specific rules. Thus, the policy
representation algorithm in \secref{sec: Computing greedy policies} can
be applied directly with our new rule-based representation.   Therefore, we
now have a complete framework for exploiting both additive and
context-specific structure for efficient planning in factored MDPs. 





% --------------------------------------------------------------
% --------------------------------------------------------------
% --------------------------------------------------------------
\section{Experimental Results}\label{sec:Experimental results}
%CGjairf4 rewrote this section


The factored representation of a value function is most
appropriate in certain types of systems:  Systems that involve
many 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~\citeyear{Simon:1981} in
``Architecture of Complexity,'' many complex systems have a
``nearly decomposable, hierarchical structure,'' with the subsystems
interacting only weakly between themselves. To evaluate our
algorithm, we selected problems that we believe 
%DK4 
%characterizes
exhibit
this type of structure.

In this section, we 
%DK4
%evaluate three aspects of our algorithms. 
perform various experiments intended to explore the performance of our
algorithms.  
First, we compare our factored approximate linear programming (LP) and
approximate policy iteration (PI) algorithms. We
also compare to the $\cL_2$-projection algorithm of Koller and
Parr~\citeyear{Koller+Parr:2000}.  Our second evaluation compares a
table-based implementation to a rule-based implementation that can
exploit CSI. Finally, we present comparisons between our approach and
the algorithms of Boutilier \etal~\citeyear{Boutilier+al:2000}. 


% --------------------------------------------------------------
% --------------------------------------------------------------
\subsection{Approximate LP and Approximate PI}\label{sec:LPxPI}

In order to compare our 
%factored approximate linear programming and approximate policy
%iteration with max-norm projection algorithms, we 
%tested the two approaches 
approximate LP and approximate PI algorithms, we tested both 
on the {\it SysAdmin\/} problem described in
detail in \secref{sec: MDP}.  This problem relates to
a system administrator who has to maintain a network of computers;
we experimented with various network architectures, shown in
\figref{fig:network topologies}. 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. Recall that the state space in this
problem grows exponentially in the number of machines in the network,
that is, a problem with $m$ machines has $2^m$ states. 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 $\discount = 0.95$.
%RP3: Added for reviewer 2
The optimal strategy for rebooting machines will depend upon the topology,
the discount factor, and the status of the machines in the network. If
machine $i$ and machine $j$ are both faulty, the benefit of rebooting $i$
must be weighed against the expected discounted impact of delaying rebooting
$j$ on $j$'s successors.  For topologies such as rings, this policy may be a function
of the status of every single machine in the network.

The basis functions used included independent indicators for each
machine, with value 1 if it is working and zero otherwise (\ie, each one is a
%RP3: scope
restricted scope function of a single variable), and the constant
basis, whose value is 1 for all states. 
We selected straightforward variable elimination orders: for the
``Star'' and ``Three Legs'' topologies, we first eliminated the
variables corresponding to computers in the legs, and the center
computer (server) was eliminated last; for ``Ring,'' we started with
an arbitrary computer and followed the ring order; for ``Ring and
Star,'' the ring machines were eliminated first and then the center
one; finally, for the ``Ring of Rings'' topology, we eliminated the
%DK4: little -> outer, larger -> inner
computers in the outer rings first and then the ones in the inner
ring.  


\commentout{
Note that the additive structure of the reward function makes it
unsuitable for the tree-structured representation used by
%RP3: Changed to 2000 reference
Boutilier \etal~\citeyear{Boutilier+al:2000}
%RP3: added for reviwer #1 Duh...
since the reward function requires space exponential in the  number of
states (a full tree) when expressed in their representation.
}



\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\epsfig{figure=Figures/running-time.eps,width=7.5cm} &
\epsfig{figure=Figures/more-models.eps,width=7.5cm} \\
\textbf{(a)} & \textbf{(b)} \\
~~ \\
\epsfig{figure=Figures/ring-x-biring.eps,width=7.5cm} &
\epsfig{figure=Figures/fitting-polynomial-time.eps,width=7.5cm} \\
\textbf{(c)} & \textbf{(d)}
\end{tabular}
\caption{(a)--(c) Running times for policy iteration with max-norm
projection on variants of the \emph{SysAdmin} problem;
(d) Fitting a polynomial to the running time for the ``Ring'' topology.}
\label{fig:running times} \vspace*{-.4cm}
\end{center}
\end{figure}


\commentout{In our experiments, most of the time was spent in Matlab
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 and linear
%programming. 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.
}

We implemented the factored policy iteration and
linear programming algorithms in Matlab, using CPLEX as the LP
solver. Experiments were performed on a Sun UltraSPARC-II, 359 MHz with 256MB of RAM.
To evaluate the complexity of the approximate policy iteration
with max-norm projection algorithm, tests were performed with
increasing the number of states, that is, increasing number of
machines on the network. \figref{fig:running times} shows the
running time for increasing problem sizes, for various
architectures. The simplest one is the ``Star,'' where the
%RP3: scope
backprojection of each basis function has scope restricted to two
%RP3: scope
variables and the largest factor in the cost network has scope
restricted to two variables. The most difficult one was the
``Bidirectional Ring,'' where factors contain five variables.
%RPfinal:  Globally move punctuation inside of quotations to conform
%          with style guidelines, as per Strunk & White pp. 36-37

Note that the number of states is growing exponentially (indicated by
the log scale in \figref{fig:running times}), but running times
increase only logarithmically in the number of states, or
polynomially in the number of variables.  We illustrate this
behavior in \figref{fig:running times}(d), where we fit a 3rd order
polynomial to the running times for the ``unidirectional ring.''
Note that the size of the problem description 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 factored algorithm empirically grows
approximately as $O\left((n\cdot\abs{A})^{1.5}\right)$, for a problem
with $n$ variables, as opposed to the exponential complexity 
--- $poly\left(2^n, \abs{A}\right)$ --- of the explicit algorithm.

\begin{figure*}
\begin{center}
\begin{tabular}{cc}
\epsfig{figure=Figures/l2-linf-comparison.eps,width=7.5cm} &
\epsfig{figure=Figures/bellman-error.eps,width=7.5cm} \\
{\textbf{(a)}} & {\textbf{(b)} } 
\end{tabular}
\caption{(a) Relative error to optimal value function
$\Vstar$ and comparison to $\cL_2$ projection for ``Ring''; (b) For large
models, measuring Bellman error after convergence.}
\label{fig:errors}
\end{center}
\end{figure*}

For further evaluation, we measured the error in our approximate
value function relative to the true optimal value function
$\Vstar$.  Note that it is only possible to compute $\Vstar$ 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 $\cL_2$-projection
algorithm of Koller and Parr~\citeyear{Koller+Parr:2000}.  As we
discussed in \secref{sec:PI comparisons}, the $\cL_2$ projections in factored MDPs by
Koller and Parr are difficult and time consuming; hence, we were
only able to compare the two algorithms for smaller problems,
where an equivalent $\cL_2$-projection can be implemented using an
explicit state space formulation.  Results for both algorithms are
presented in \figref{fig:errors}(a), 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 $\Vstar$ than the
$\cL_2$-projection. Here, we used two types of basis functions:
the same single variable functions, and pairwise basis functions. The pairwise basis functions contain
indicators for neighboring pairs of machines (\ie,
functions of two variables). As expected, the use of pairwise basis
functions resulted in 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 $\Vstar$. 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 $\cL_2$
projection.  In a ``Unidirectional Ring'' with 8 machines and
pairwise basis, the relative error between our approximation and
$\Vstar$ was about $10\%$, but the resulting policy only had a
$6\%$ loss over the optimal policy.  For the same problem, the
$\cL_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 $\cL_2$ algorithm, even
with respect to actual policy loss.

For large models, we can no longer compute the correct value
function, so we cannot evaluate our results by computing
%RP3:  Removed redundant definition as per reviewer complaint
%$\maxnorm{\Vstar - Aw}$. However, the \emph{Bellman error},
%defined as $\Bellman(\V) = \maxnorm{\Tbell \V - \V}$, can be used
%to provide a bound: $\maxnorm{\Vstar - Aw} \leq
%\frac{\Bellman(\Basismatrix\bw)}{1 - \discount}$~\cite{Williams+Baird:1993}.
$\maxnorm{\Vstar - \Basismatrix\bw}$. Fortunately, as discussed in \secref{sec: Bellman error}, the Bellman error
can be used to provide a bound on the approximation error and can be computed efficiently by exploiting problem-specific structure.
%Thus, we use the Bellman error to evaluate our answers for larger models. 
\figref{fig:errors}(b) shows that the Bellman error increases very
slowly with the number of states.

It is also valuable to look at the actual decision-list policies generated in our
experiments. First, we noted that the lists 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.

\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\epsfig{figure=Figures/lp-versus-pi-runtime.eps,width=7.5cm} &
\epsfig{figure=Figures/lp-versus-pi-value.eps,width=7.5cm} \\
{\bf (a)} &
{\bf(b)}
\end{tabular}
\caption{Approximate LP versus approximate PI on the {\it SysAdmin\/} problem with a ``Ring'' topology: (a) running time; (b) estimated value of policy.}
\label{fig:lp versus pi}
\end{center}
\end{figure}

Thus far, we have presented scaling results for running times and
approximation error for our approximate PI approach.  We now compare
this algorithm to the simpler approximate LP approach of
\secref{sec:ALP}. As shown in \figref{fig:lp versus pi}(a), the
approximate LP algorithm for factored MDPs is significantly faster
than the approximate PI algorithm.  In fact, approximate PI with
single-variable basis functions variables is more costly
computationally than the LP approach using basis functions over
consecutive triples of variables.  As shown in \figref{fig:lp versus
pi}(b), for singleton basis functions, the approximate PI policy
obtains slightly better performance for some problem sizes. However,
as we increase the number of basis functions for the approximate LP
formulation, the value of the resulting policy is much better. Thus,
in this problem, our factored approximate linear programming
formulation allows us to use more basis functions and to obtain a
resulting policy of higher value, while still maintaining a faster
running time. These results, along with the simpler implementation,
suggest that in practice one may first try to apply the approximate
linear programming algorithm before deciding to move to the more
elaborate approximate policy iteration approach.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Comparing Table-based and Rule-based Implementations} \label{sec: tablexcsi results}


\begin{figure}
\begin{center}
\begin{tabular}{c}
\epsfig{figure=Figures/star-process-time.eps,width=8.5cm} \\
{\textbf{(a)}} \\
\epsfig{figure=Figures/ring-process-time.eps,width=8.5cm} \\
{\textbf{(b)} } \\
\epsfig{figure=Figures/r-star-process-time.eps,width=8.5cm} \\
{\textbf{(c)} }
\end{tabular}
\caption{Running time for {\it Process-SysAdmin\/} problem for various
topologies: (a) ``Star''; (b) ``Ring''; (c) ``Reverse star'' 
%DK4 
(with fit function).}
\label{fig: Process-SysAdmin time}
\end{center}
\end{figure}



\begin{figure}[t]
\begin{center}
\begin{tabular}{c}
\epsfig{figure=Figures/cplex-total-time-ring-process.eps,width=8.5cm} \\
\end{tabular}
\caption{Fraction of total running time spent in CPLEX for the {\it Process-SysAdmin\/} problem with a ``Ring'' topology.}
\label{fig: Process-SysAdmin cplex x total time}
\end{center}
\end{figure}

Our next evaluation compares a table-based representation, which exploits only additive independence, to the rule-based representation presented in \secref{sec: CSI}, which can exploit both additive and context-specific independence. For these experiments,
we implemented our factored approximate linear programming algorithm with table-based and rule-based representations in C++, using CPLEX as the LP
solver. Experiments were performed on a Sun UltraSPARC-II, 400 MHz with 1GB of RAM. 

To evaluate and compare the algorithms, we utilized a more complex
extension of the {\it SysAdmin\/} problem. This problem, dubbed the {\it
Process-SysAdmin\/} problem, contains three state variables for each
machine $i$ in the network: $\textit{Load}_i$, $\textit{Status}_i$ and
$Selector_i$. Each computer runs processes and receives rewards when
the processes terminate. These processes are represented by the
$\textit{Load}_i$ variable, which takes values in $\{\textit{Idle},
\textit{Loaded}, \textit{Success}\}$, and the computer receives a
reward when the assignment of $\textit{Load}_i$ is $Success$. The
$\textit{Status}_i$ variable, representing the status of machine $i$,
takes values in \textit{\{Good, Faulty, Dead\}}; if its value is
\textit{Faulty}, then processes have a smaller probability of
terminating and if its value is $\textit{Dead}$, then any running
process is lost and $\textit{Load}_i$ becomes \textit{Idle}. The
status of machine $i$ can become \textit{Faulty} and eventually
\textit{Dead} at random; however, if machine $i$ receives a packet
from a dead machine, then the probability that $\textit{Status}_i$
becomes \textit{Faulty} and then \textit{Dead} increases. The
$\textit{Selector}_i$ variable represents this communication by
selecting one of the neighbors of $i$ uniformly at random at every
time step. The SysAdmin can select at most one computer to reboot at
every time step. If computer $i$ is rebooted, then its status becomes
\textit{Good} with probability $1$, but any running process is lost,
\ie, the $\textit{Load}_i$ variable becomes $\textit{Idle}$. Thus, in
this problem, the SysAdmin must balance several conflicting goals:
rebooting a machine kills processes, but not rebooting a machine may
cause cascading faults in network. Furthermore, the SysAdmin can only
choose one machine to reboot, which imposes the additional tradeoff of
selecting only one of the (potentially many) faulty or dead machines
in the network to reboot.

We experimented with two types of basis functions: ``single+''
includes indicators over all of the joint assignments of
$\textit{Load}_i$, $\textit{Status}_i$ and $Selector_i$, and ``pair''
which, in addition, includes a set of indicators over
$\textit{Status}_i$, $\textit{Status}_j$, and $Selector_i = j$, for
each neighbor $j$ of machine $i$ in the network. The discount factor
was $\discount = 0.95$. The variable elimination order eliminated all
of the $\textit{Load}_i$ variables first, and then followed the same
patterns as in the simple \textit{SysAdmin} problem, eliminating first
$\textit{Status}_i$ and then $\textit{Selector}_i$ when machine $i$ is eliminated. 

\figref{fig: Process-SysAdmin time} compares the running times for the
table-based implementation to the ones for the rule-based
representation for three topologies: ``Star,'' ``Ring,'' and ``Reverse
star.'' The ``Reverse star'' topology reverses the direction of the
influences in the ``Star'': rather than the central machine
influencing all machines in the topology, all machines influence the
central one. These three topologies demonstrate three different levels
of CSI: In the ``Star'' topology, the factors generated by variable
elimination are small. Thus, although the running times are polynomial
in the number of state variables for both methods, the table-based
representation is significantly faster than the rule-based one, due to
the overhead of managing the rules. The ``Ring'' topology illustrates
an intermediate behavior: ``single+'' basis functions induce
relatively small variable elimination factors, thus the table-based
approach is faster. However, with ``pair'' basis the factors are
larger and the rule-based approach starts to demonstrate faster
running times in larger problems. Finally, the ``Reverse star''
topology represents the worst-case scenario for the table-based
approach. Here, the scope of the backprojection of a basis function
for the central machine will involve all computers in the network, as
all machines can potentially influence the central one in the next
time step. Thus, the size of the factors in the table-based variable
elimination approach are exponential in the number of machines in the
network, which is illustrated by the exponential growth in
\figref{fig: Process-SysAdmin time}(c). The rule-based approach can
exploit the CSI in this problem; for example, the status of the central
machine $\textit{Status}_0$ only depends on machine $j$ if the value selector
is $j$, \ie, if $\textit{Selector}_0 = j$. By exploiting CSI, we can solve the
same problem in polynomial time in the number of state variables, as
seen in the second curve in \figref{fig: Process-SysAdmin time}(c).

It is also instructive to compare the portion of the total running
time spent in CPLEX for the table-based as compared to the rule-based
approach. \figref{fig: Process-SysAdmin cplex x total time}
illustrates this comparison. Note that amount of time spent in CPLEX
is significantly higher for the table-based approach. There are two
reasons for this difference: first, due to CSI, the LPs generated by
the rule-based approach are smaller than the table-based ones; second,
rule-based variable elimination is more complex than the table-based
one, due to the overhead introduced by rule management. Interestingly,
the proportion of CPLEX time increases as the problem size increases,
indicating that the asymptotic complexity of the LP solution is higher
than that of variable elimination, thus suggesting that, for larger
%RPn:  Added constraint generation
problems, additional large-scale LP optimization procedures, such as
constraint generation, may be helpful.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Comparison to Apricodd} \label{sec: versus APRICODD} 



\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\epsfig{figure=Figures/linear-fit.eps,width=7.5cm} &
\epsfig{figure=Figures/expon-fit.eps,width=7.5cm} \\
{\textbf{(a)}} & {\textbf{(b)}} 
\end{tabular}
\caption{Comparing Apricodd to rule-based approximate linear programming on the (a) {\it Linear\/} and (b) {\it Expon\/} problems.}
\label{fig: linear, expon time}
\end{center}
\end{figure}


The most closely related work to ours is a line of research that began with
the work of Boutilier \etal~\citeyear{Boutilier+al:1995}. 
In particular, the approximate Apricodd algorithm of Hoey 
\etal~\citeyear{Hoey+al:1999}, which uses  
analytic decision diagrams (ADDs) to represent the value function is a strong alternative approach for solving factored MDPs. As discussed in detail in \secref{sec:Related}, the Apricodd algorithm can successfully exploit context-specific structure in the \emph{value function}, by representing it with the set of mutually-exclusive and exhaustive branches of the ADD. On the other hand, our approach can exploit both additive and context-specific structure in the problem, by using a linear combination of non-mutually-exclusive rules. To better understand this difference, we evaluated both our rule-based approximate linear programming algorithm and Apricodd in two problems, {\it Linear\/} and {\it Expon\/}, designed by Boutilier \etal~\citeyear{Boutilier+al:2000} to illustrate respectively the best-case and the worst-case behavior of their algorithm.
In these experiments, we used the web-distributed version of Apricodd~\cite{Hoey+al:2002}, running it locally on a Linux Pentium III 700MHz with 1GB of RAM.

These two problems involve $n$ binary variables $X_1, \ldots, X_n$ and $n$
deterministic actions $a_1, \ldots, a_n$. The reward is $1$ when all
variables $X_k$ are \true, and is $0$ otherwise. The problem is discounted
by a factor $\discount = 0.99$. The difference between the {\it Linear\/}
and the {\it Expon\/} problems is in the transition probabilities. In the
{\it Linear\/} problem, the action $a_k$ sets the variable $X_k$ to \true\
and makes all \emph{succeeding} variables, $X_i$ for $i>k$, \false.  If the
state space of the {\it Linear\/} problem is seen as a binary number, the
optimal policy is to set repeatedly the largest bit ($X_k$ variable) which
has all preceding bits set to \true. Using an ADD, the optimal value
function for this problem can be represented in linear space, with $n+1$
leaves~\cite{Boutilier+al:2000}. This is the ``best-case'' for Apricodd,
and the algorithm can compute this value function quite
efficiently. \figref{fig: linear, expon time}(a) compares the running time
%RPfinal:  Reworded below.
of Apricodd to that of one of our algorithms with indicator basis functions
between pairs of consecutive variables. Note that both algorithms obtain
the same policy in polynomial time in the number of variables. However, in
such structured problems, the efficient implementation of the ADD package
used in Apricodd makes it faster in this problem.

On the other hand, the {\it Expon\/} problem illustrates the worst-case for Apricodd. 
In this problem, the action $a_k$ sets the variable $X_k$ to \true, if all \emph{preceding} variables, $X_i$ for $i<k$,  are \true, and it makes all preceding variables \false.
If the state space is seen as a binary number, the optimal policy goes
through all binary numbers in sequence, by repeatedly setting the
largest bit ($X_k$ variable) which has all preceding bits set to
\true. Due to discounting, the optimal value function assigns a value
of $\discount^{2^n -j-1}$ to the $j$th binary number, so that the value
function contains exponentially many different values. Using an ADD,
the optimal value function for this problem requires an exponential
number of leaves~\cite{Boutilier+al:2000}, which is illustrated by the  
exponential running time in \figref{fig: linear, expon time}(b).
However, the same value function can be approximated very compactly 
%DK4
%by exploiting additive structure: . Using simple indicators over each 
%variable $X_k$ and along with the constant as basis functions (\ie,
%$n+1$ basis functions), we can represent exponentially many different
%values. 
as a factored linear value function using $n+1$ basis functions: an
indicator over each variable $X_k$ and the constant base.
As shown in \figref{fig: linear, expon time}(b), using this
representation, our factored approximate linear programming algorithm
computes the value function in polynomial time. Furthermore, the
policy obtained by our approach was optimal for this problem. 
%DK4
%Thus indicating that the additive independence exploited by our
%algorithm is key to obtaining a polynomial-time solution to this
%problem. 
Thus, in this problem, the ability to exploit additive independence
allows an efficient polynomial time solution.




\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\epsfig{figure=Figures/ring-rule-appricod-time.eps,width=7.5cm} &
\epsfig{figure=Figures/ring-rule-appricod-value.eps,width=7.5cm} \\
{\textbf{(a)}} & {\textbf{(b)}} \\
\epsfig{figure=Figures/star-rule-appricod-time.eps,width=7.5cm} &
\epsfig{figure=Figures/star-rule-appricod-value.eps,width=7.5cm} \\
{\textbf{(c)}} & {\textbf{(d)}} 
\end{tabular}
\caption{Comparing Apricodd to rule-based approximate linear programming with ``single+'' basis functions on the {\it Process-SysAdmin\/} problem with ``Ring'' topology (a) running time and (b) value of the resulting policy; and with ``Star'' topology (c) running time and (d) value of the resulting policy.}
\label{fig: ring, star, ruleXappricod time,value}
\end{center}
\end{figure}


We have also compared Apricodd to our rule-based approximate linear
programming algorithm on the {\it Process-SysAdmin\/} problem. This
problem has significant additive structure in the reward function and
factorization in the transition model.  Although this type of
structure is not exploited directly by Apricodd, the ADD approximation
steps performed by the algorithm can, in principle, allow Apricodd to
find approximate solutions to the problem. We spent a significant
amount of time attempting to find the best set of parameters for
Apricodd for these problems.%
\footnote{We are very grateful to Jesse Hoey and Robert St-Aubin for their assistance in selecting the
parameters.}
We settled on the ``sift'' method of variable reordering and the
``round'' approximation method with the ``size''  (maximum ADD size)
criteria.  To allow the value function representation to scale with
the problem size, we set the maximum ADD size to $4000 + 400n$ for a
network with $n$ machines.  (We experimented with a variety of
different growth rates for the maximum ADD size; here, as for the
other parameters, we selected the choice that gave the best results
for Apricodd.)  We compared Apricodd with these parameters to our 
rule-based approximate linear programming algorithm with ``single+''
basis functions on a Pentium III 700MHz with 1GB of RAM. These results
are summarized in \figref{fig: ring, star, ruleXappricod time,value}.  

On very small problems (up to 4--5 machines), the performance of the
two algorithms is fairly similar in terms of both the running time and
the  quality of the policies generated.  However, as the problem size
grows, the running time of Apricodd increases rapidly, and becomes
significantly higher than that of our algorithm .  Furthermore, as the
problem size increases, the quality of the policies generated by
Apricodd also deteriorates.  This difference in policy quality is
caused by the different value function representation used by the two
algorithms.  The ADDs used in Apricodd represent $k$ different values
with $k$ leaves; thus, they are forced to agglomerate many different
states and represent them using a single value.  For smaller problems,
such agglomeration can still represent good policies. Unfortunately,
as the problem size increases and the state space grows exponentially,
Apricodd's policy representation becomes inadequate, and the quality
of the policies decreases. On the other hand, our linear value
functions can represent exponentially many values with
only $k$ basis functions, which allows our approach to scale up to
significantly larger problems. 



% -------------------------------------------------------------
% -------------------------------------------------------------
\section{Related Work}\label{sec:Related}

The most closely related work to ours is a line of research that began with
the work of Boutilier \etal~\citeyear{Boutilier+al:1995}.  We 
address this comparison separately below, but we begin this section with
some broader background references.

\subsection{Approximate MDP Solutions}

The field of MDPs, as it is popularly known, was formalized by
Bellman~\citeyear{Bellman:1957} in the 1950's.  The 
importance of value function approximation was recognized at an early stage
by Bellman himself~\citeyear{Bellman+al:1963}.  In the early 1990's the MDP
framework was recognized by AI researchers as a formal framework that could
be used to address the problem of planning under uncertainty
\cite{Dean+al:1993}.

Within the AI community, value function approximation developed
concomitantly with the notion of value function representations for Markov
chains.  Sutton's seminal paper on temporal difference
learning~\citeyear{Sutton:1988}, which addressed the use of value
functions for prediction but not planning, assumed a very general
representation of the value function and noted the connection to general
function approximators such as neural networks.  However, the stability of
this combination was not directly addressed at that time.

Several important developments gave the AI community deeper insight into
the relationship between function approximation and dynamic programming.
Tsitsiklis and Van Roy~\citeyear{Tsitsiklis+VanRoy:1996b} and,
independently,
Gordon~\citeyear{Gordon:1995} 
%DK4 - are we claiming that they were given credit wrongly?!
%are given credit for popularizing the
popularized the analysis of approximate MDP methods via the
contraction properties of the dynamic programming operator and
function approximator.  Tsitsiklis and Van
Roy~\citeyear{Tsitsiklis+VanRoy:1996a} later established a general 
convergence result for linear value function approximators and
$TD(\lambda)$, and Bertsekas and
Tsitsiklis~\citeyear{Bertsekas+Tsitsiklis:1996} unified a large body
of work on approximate dynamic programming under the name of 
\emph{Neuro-dynamic Programming}, also providing many novel
and general error analyses.

Approximate linear programming for MDPs using linear value function
approximation was introduced by Schweitzer and
Seidmann~\citeyear{Schweitzer+Seidmann:1985}, although the approach was
somewhat deprecated until fairly recently due the lack of compelling error
analyses and the lack of an
%DK4
%compelling
effective method for handling the large number of constraints.  Recent
work by de Farias and Van 
Roy~\citeyear{Farias+VanRoy:2001,Farias+VanRoy:2001c} has started to
address these concerns with new error bounds and constraint sampling
methods.  
%In contrast, 
Our approach, rather than sampling constraints, utilizes structure in 
the model and value function to represent all of the constraints
compactly. 
% emphasizes efficient methods of representing \emph{all} constraints
% given a factored model and value function representation.

\subsection{Factored Approaches}

Tatman and Shachter~\citeyear{Tatman+Shachter:1990} considered the
additive decomposition of value nodes in influence diagrams.  A number of
approaches to factoring of general MDPs have been explored in the
literature.  Techniques for exploiting reward functions that decompose
additively were studied by Meuleau \etal~\citeyear{Meuleau+al:1998}, and
by Singh and Cohn~\citeyear{Singh+Cohn:1998}.

%DK4 - rewrote 
The use of factored representations such as dynamic Bayesian networks was
pioneered by Boutilier \etal~\citeyear{Boutilier+al:1995} and has
developed steadily in recent years.   These methods rely on the use of
context-specific structures such as decision trees or analytic decision
diagrams (ADDs)~\cite{Hoey+al:1999} to represent both the transition
dynamics of the DBN and the value function.  The algorithms use 
dynamic programming to partition the state space, representing the
partition using a tree-like structure that branches on state
variables and assigns values at the leaves.  The tree is grown
dynamically as part of the dynamic programming process and the 
algorithm creates new leaves as needed:  A leaf is split by the
application of a DP operator when two states associated with that leaf
turn out to have different values in the backprojected value function.
This process can also be interpreted as a form of model minimization
\cite{Dean+Givan:1997}. 

The number of leaves in a tree used to represent a value function
determines the computational complexity of the algorithm.  It also limits
the number of distinct values that can be assigned to states: since the
leaves represent a partitioning of the state space, every state maps to
exactly one leaf.  
%DK3
%It is easy to construct trivial MDPs which will require exponentially
%large value functions~\cite{Koller+Parr:1999}.  This general difficulty
%was recognized early in the line of research on Factored MDPs. 
However, as was recognized early on, there are trivial MDPs which require
exponentially large value functions.  This observation led to a line of
approximation algorithms aimed at limiting the tree
size~\cite{Boutilier+Dearden:1996} and, later, limiting the ADD 
size~\cite{StAubin+al:2000}.  Kim and Dean~\citeyear{Kim+Dean:2001} also
explored techniques for discovering tree-structured value functions for
factored MDPs.  While these methods permit good approximate solutions to
some large MDPs, their complexity is still determined by the number of
leaves in the representation and the number of distinct values than can be
assigned to states is still limited as well. 

Tadepalli and Ok~\citeyear{Tadepalli+Ok:1996} were the first to apply
linear value function approximation to Factored MDPs.  Linear value
function approximation is a potentially more expressive approximation
method because it can assign unique values to every state in an MDP without
requiring storage space that is exponential in the number of state
variables.  The expressive power of a tree with $k$ leaves can be captured
by a linear function approximator with $k$ basis functions such that basis
function $h_i$ is an indicator function that tests if a state belongs in
the partition of leaf $i$.  Thus, the set of value functions that can be
represented by a tree with $k$ leaves is a subset of the set of value
functions that can be represented by a value function with $k$ basis
functions.
%CGjairf4 added:
Our experimental results in \secref{sec: versus APRICODD} highlight
%DK4 - only one problem...
this difference by showing an example problem that requires
exponentially many leaves in the value function, but that can be
approximated 
%compactly using additive structure.
well using a linear value function.

The main advantage of tree-based value functions is that their structure is
determined dynamically during the solution of the MDP.  In principle, as
the value function representation is derived automatically from the model
description, this approach requires less insight from the user.  In
problems for which the value function can be well approximated by a
relatively small number of values, this approach provides an excellent
solution to the problem.  Our method of linear value function
approximation aims to address what we believe to be the more common case,
where a large range of distinct values is required to achieve a good
approximation.

%CGjairf1 I would like to rephrase this, ask Ron.
%RP3:  Added as per reviwer 1's comments
%RP4:  Rephrased
%DK3: rephrased again
\commentout{
Finally, we note that Schuurmans and
Patrascu~\citeyear{Schuurmans+Patrascu:2001} developed an alternative
adaptation of the cost network from Guestrin \etal~\citeyear{Guestrin+al:2001b} to approximate linear programming using
constraint generation.}
Finally, we note that Schuurmans and
Patrascu~\citeyear{Schuurmans+Patrascu:2001}, based on our earlier work
on max-norm projection using cost networks and linear programs,
independently developed an alternative approach to approximate linear
programming using a cost network.  Our method embeds a cost network
inside a single linear program.  
%DK4
%Their method does not contain a cost network inside the LP, 
By contrast, their method is based on a constraint generation
approach, using a cost network to detect constraint
violations.  When constraint violations are found, a new constraint is
added, repeatedly generating and attempting to solve LPs until a feasible
solution is found.
%CGJairf4 added comment illustrated by our experimental results
Interestingly, as the approach of Schuurmans and
Patrascu uses multiple calls to variable elimination in order to speed
up the LP solution step, it will be most successful when the time spent solving the LP is significantly larger than the time required for variable elimination. As suggested in \secref{sec: tablexcsi results}, the LP solution time is larger for the table-based approach. Thus, Schuurmans and
Patrascu's constraint generation method will probably be more successful in table-based problems than in rule-based ones.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}\label{sec:Conclusion}
%DK4 - this was way too long, and reads more like a rehash of
%everything than a summary; I cut large chunks.  Old version is
%at the end
%CGjairf2 rewrote entire conclusion section

In this paper, we presented new algorithms for approximate linear
programming and approximate dynamic programming (value and policy
iteration) for factored MDPs.  Both of these algorithms leverage on a novel
LP decomposition technique, analogous to variable elimination in cost
networks, which reduces an exponentially large LP to a provably
equivalent, polynomial-sized one.  

Our approximate dynamic programming algorithms are motivated by error
analyses showing the importance of minimizing $\Linf$ error.  These
algorithms are more efficient and substantially easier to implement
than previous algorithms based on the $\cL_2$-projection.  Our
experimental results suggest that they also perform better in
practice.  

Our approximate linear programming algorithm for factored MDPs is simpler,
easier to implement and more general than the dynamic programming
approaches. Unlike our policy iteration algorithm, it does not rely on 
the default action assumption, which states that actions only affect a
small number of state variables.  Although this algorithm does not
have the same theoretical guarantees as max-norm projection
approaches, empirically it seems to be a favorable option.  Our
experiments suggest that approximate policy iteration tends to
generate better policies for the same set of basis
functions.  However, due to the computational advantages, we can
add more basis functions to the approximate linear programming algorithm,
obtaining a better policy and still maintaining a much faster running time
than approximate policy iteration.

Unlike previous approaches, our
algorithms can exploit both additive and context-specific structure in
the factored MDP model. Typical real-world systems possess both of
these types of structure. thus, this feature of our algorithms will
increase the applicability of factored MDPs to more practical
problems.  We demonstrated that exploiting context-specific
independence, by using a rule-based representation instead of the
standard table-based one, can yield exponential improvements in
computational time when the problem has significant amounts of
CSI. However, the overhead of managing sets of rules make it less
well-suited for simpler problems.  We also compared our approach to
the work of Boutilier \etal~\citeyear{Boutilier+al:2000}, which
exploits only context-specific structure.  For problems with significant
context-specific structure in the value function, their approach can
be faster due to their efficient handling of the ADD
representation. However, there are problems with significant
context-specific structure in the problem representation, rather than
in the value function, which require exponentially large ADDs. In some
such problems, we demonstrated that by using a linear value function
our algorithm can obtain a polynomial-time near-optimal approximation
of the true value function.

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 algorithms currently require the designer to specify the
factored basis functions.  This is a limitation compared to the
algorithms of Boutilier \etal~\citeyear{Boutilier+al:2000}, which are fully automated.
However, our experiments suggest that a few simple rules can be
quite successful for designing a basis. First, we ensure that
the reward function is representable by our basis.
A simple basis that, in addition, contained a separate set of
indicators for each variable often did quite well.
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 functions
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.

Nonetheless, a general algorithm for choosing good factored basis
functions still does not exist. However, there are some potential
approaches: First, in problems with CSI, one could apply the
algorithms of Boutilier \etal\ for a few iterations to
%RPn:  Added backprojections (else basis functions would be orthogonal)
generate partial tree-structured solutions. Indicators defined over
the variables in backprojection of the leaves could, in turn, be used
to generate a basis set for such problems. Second, the Bellman error
computation, which can be performed efficiently as shown in
\secref{sec: Bellman error}, does not only provide a bound on the
quality of the policy, but also the actual state where the error is
largest. This knowledge can be used to create a mechanism to
incrementally increase the basis set, adding new basis functions to
tackle states with high Bellman error.

There are many other possible extensions to this work. We have
already pursued extensions to collaborative multiagent systems,
where multiple agents act simultaneously to maximize the global
reward~\cite{Guestrin+al:2001c}, and factored POMDPs, where the
full state is not observed directly, but indirectly through
observation variables~\cite{Guestrin+al:2001d}. However, there are
other settings that remain to be explored. In particular, we hope
to address the problem of learning a factored MDP and planning in
a competitive multiagent system. 

Additionally, in this paper we have tackled problems where the induced
width of the cost network is sufficiently low or that possess sufficient
context-specific structure to allow for the exact solution of our factored
LPs. Unfortunately, some practical problems may have prohibitively large
induced width. We plan to leverage on ideas from loopy belief propagation
algorithms for approximate inference in Bayesian
networks~\cite{Pearl:1988,Yedidia+al:2000} to address this issue.

We believe that the methods described herein significantly further
extend the efficiency, applicability and general usability of factored
models and value functions for the control of practical dynamic systems.

% ------------------------------------------------------------
% ------------------------------------------------------------
\section*{Acknowledgements}
We are very grateful to Craig Boutilier, Dirk Ormoneit and Uri Lerner
 for many useful discussions, and to the anonymous reviewers for their detailed and thorough comments. We also would like to thank Jesse Hoey,
 Robert St-Aubin, Alan Hu, and Craig Boutilier for distributing their
 algorithm and for their very useful assistance in using Apricodd and in selecting its parameters.  This work was supported by the DoD MURI program,
 administered by the Office of Naval Research under Grant
 N00014-00-1-0637, by Air Force contract F30602-00-2-0598 under
 DARPA's TASK program, and by the Sloan Foundation. The first author
 was also supported by a Siebel Scholarship. 






\appendix

\section{Proofs}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Proof of \lemref{lem: vd error is bounded}} \label{proof:lem: vd error is bounded}
There exists a setting to the weights --- the all zero setting 
--- that
yields a bounded max-norm projection error $\betap$ for any policy
($\betap \leq \Rmax$). Our max-norm projection operator chooses
the set of weights that minimizes the projection error
$\beta\at{t}$ for each policy $\policy\at{t}$. Thus, the
projection error $\beta\at{t}$ must be at least as low as the one
given by the zero weights $\betap$ (which is bounded). Thus, the
error remains bounded for all iterations. \qed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
\subsection{Proof of \thmref{thm: pi error bound}}
\label{proof:thm: pi error bound}
First, we need to bound our  approximation of
$\V_{\policy\at{t}}$:
\beqa
\maxnorm{\V_{\policy\at{t}}  - \Basismatrix \bw\at{t}} & \leq &
     \maxnorm{\T_{\policy\at{t}} \Basismatrix \bw\at{t} - \Basismatrix \bw\at{t} } +
     \maxnorm{\V_{\policy\at{t}} - \T_{\policy\at{t}} \Basismatrix \bw\at{t}}
        \mbox{\small~~~;~(triangle inequality;)} \\
     & \leq &   \maxnorm{\T_{\policy\at{t}} \Basismatrix \bw\at{t} - \Basismatrix \bw\at{t}}
                + \discount\maxnorm{\V_{\policy\at{t}} - \Basismatrix \bw\at{t}}
        \mbox{~~~;~($\T_{\policy\at{t}}$ is a contraction.)}
\end{eqnarray*}
Moving the second term to the right hand side and dividing through by
$1-\discount$, we obtain:
\beqn\label{eqn:value-determination-error}
\maxnorm{\V_{\policy\at{t}}  - \Basismatrix \bw\at{t}} \leq \frac{1}{1-\discount}
     \maxnorm{\T_{\policy\at{t}} \Basismatrix \bw\at{t} - \Basismatrix \bw\at{t}}
    = \frac{\beta\at{t}}{1-\discount}.
\end{equation}
For the next part of the proof, we adapt a lemma of Bertsekas and
Tsitsiklis~\citeyear[Lemma 6.2, p.277]{Bertsekas+Tsitsiklis:1996}
to fit into our framework. After some manipulation, this lemma can
be reformulated as: \beqn \maxnorm{\Vstar - \V_{\policy\at{t+1}}}
\leq
    \discount \maxnorm{\Vstar - \V_{\policy\at{t}}} +
    \frac{2\discount}{1-\discount}\maxnorm{\V_{\policy\at{t}}  - \Basismatrix\bw\at{t}}. \label{eqn:Tsitsiklis lemma}
\end{equation}
The proof is concluded by substituting
\eqnref{eqn:value-determination-error} into \eqnref{eqn:Tsitsiklis lemma} and, finally, induction on $t$. \qed


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Proof of \thmref{thm: factored LP}} \label{proof:thm: factored LP}
First, note that the equality constraints represent a
simple change of variable. Thus, we can rewrite \eqnref{eq: max
constraint form} in terms of these new LP variables
$u^{f_i}_{\bz_i}$ as: 
\beqn \label{eq: lp var max constraint form}
\phi \geq  \max_{\bx} \sum_i u^{f_i}_{\bz_i},
\eeqn 
where any assignment to the weights $\bw$ implies an assignment for each 
$u^{f_i}_{\bz_i}$.  After this stage, we only have LP variables.

It remains to show that the
factored LP construction is equivalent to the constraint in
\eqnref{eq: lp var max constraint form}. For a system with $n$
variables $\{X_1,\dots,X_n\}$, we assume, without loss of
generality, that variables are eliminated starting from $X_n$ down
to $X_1$. We now prove the equivalence by induction on the
number of variables.

%CGjairf2 fixed
%The base case is $n=1$, 
The base case is $n=0$, so that the functions $c_i(\bx)$ and $b(\bx)$ in 
\eqnref{eq: max constraint form} all have empty scope.  In this case,
\eqnref{eq: lp var max constraint form} can be written as: 
\beqn\label{eq: base case max} 
\phi \geq  \sum_i u^{e_i}.
\eeqn 
In this case, no transformation is done on the constraint, and
equivalence is immediate.

Now, we assume the result holds for systems with $i-1$ variables
and prove the equivalence for a system with $i$ variables. In such
a system, the maximization can be decomposed into two terms: one
%RP1:
%RP with the factors that \emph{do not} depend on $X_i$, these are
with the factors that \emph{do not} depend on $X_i$, which are
irrelevant to the maximization over $X_i$, and another term with
all the factors that depend on $X_i$. Using this decomposition, we
can write \eqnref{eq: lp var max constraint form} as:
\begin{eqnarray} 
\phi & \geq &  \max_{x_1,\dots,x_i} \sum_j u^{e_j}_{\bz_j}; \nonumber\\
& \geq & \max_{x_1,\dots,x_{i-1}} \left[ \sum_{l\ :\ x_i \not\in
\bz_l} u^{e_l}_{\bz_l} + \max_{x_i} \sum_{j\ :\ x_i \in \bz_j}
u^{e_j}_{\bz_j} \right]. \label{eq: induction step max in}
 \end{eqnarray}
At this point we can define new LP variables $u^e_{\bz}$
corresponding to the second term on the right hand side of the
constraint. These new LP variables must satisfy the following
constraint: 
\beqn \label{eq: induction step max}
  u^e_{\bz} \geq \max_{x_i} \sum_{j=1}^\ell u^{e_j}_{\instin{\bZ_j}{(\bz,x_i)}}.
\eeqn 
This new non-linear constraint is again represented in the
factored LP construction by a set of equivalent linear
constraints: 
\beqn \label{eq: induction step explicit}
  u^e_{\bz} \geq \sum_{j=1}^\ell u^{e_j}_{\instin{\bZ_j}{(\bz,x_i)}}, \forall \bz,x_i.
\eeqn 
%CGjairf2 removed
%As in the base case, 
The equivalence between the non-linear
constraint \eqnref{eq: induction step max} and the set of linear
constraints in \eqnref{eq: induction step explicit} can be shown
by considering binding constraints. For each new LP variable
created $u^e_{\bz}$, there are $\abs{X_i}$ new constraints
created, one for each value $x_i$ of $X_i$. For any assignment to
the LP variables in the right hand side of the constraint in
\eqnref{eq: induction step explicit}, only one of these
$\abs{X_i}$ constraints is relevant. That is, one where
$\sum_{j=1}^\ell u^{e_j}_{\instin{\bZ_j}{(\bz,x_i)}}$ is maximal, which
corresponds to the maximum over $X_i$. Again, if
for each value of $\bz$ more than one assignment to $X_i$ achieves
the maximum, then any of (and only) the constraints corresponding
to those maximizing assignments could be binding. Thus,
\eqnref{eq: induction step max} and \eqnref{eq: induction step
explicit} are equivalent.

Substituting the new LP variables $u^e_{\bz}$ into \eqnref{eq:
induction step max in}, we get: \bdm \phi \geq
\max_{x_1,\dots,x_{i-1}} \sum_{l\ :\ x_i \not\in \bz_l}
u^{e_l}_{\bz_l} + u^e_{\bz}, \edm which does not depend on $X_i$
anymore. Thus, it is equivalent to a system with $i-1$ variables,
concluding the induction step and the proof. \qed


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Proof of \lemref{lem: convergent policy iteration}}
\label{proof:lem: convergent policy iteration}
First note that at iteration $t+1$ the objective function
$\phi\at{t+1}$ of the max-norm projection LP is given by: \bdm
\phi\at{t+1} = \maxnorm{\Basismatrix \bw\at{t+1} - \left(R_{\policy\at{t+1}}
+ \discount P_{\policy\at{t+1}}\Basismatrix\bw\at{t+1}\right) }. \edm
However, by convergence the value function estimates are equal for
both iterations: \bdm \bw\at{t+1} = \bw\at{t}. \edm So we have
that: \bdm \phi\at{t+1} = \maxnorm{\Basismatrix \bw\at{t} -
\left(R_{\policy\at{t+1}} + \discount
P_{\policy\at{t+1}}\Basismatrix\bw\at{t}\right) }. \edm In operator notation,
this term is equivalent to:  \bdm \phi\at{t+1} = \maxnorm{\Basismatrix
\bw\at{t} - \T_{\policy\at{t+1}}\Basismatrix\bw\at{t}}. \edm Note that,
$\policy\at{t+1} = \greedy{\Basismatrix\bw\at{t}}$ by definition. Thus, we
have that: \bdm \T_{\policy\at{t+1}}\Basismatrix\bw\at{t} = \Tbell
\Basismatrix\bw\at{t}. \edm Finally, substituting into the previous
expression, we obtain the result: \bdm \phi\at{t+1} = \maxnorm{\Basismatrix
\bw\at{t} - \Tbell \Basismatrix\bw\at{t}}. \qedmath \edm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% -------------------------------------------------------------
\bibliographystyle{theapa-italics}
%\denselist
\bibliography{mdp}



\end{document}








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{table}[t]
\begin{center}
\begin{tabular}{|l|c|c|}
\hline
Topology & \# vars. in backprojection & \# vars. in largest VE factor\\
\hline
\hline
Star & 2 & 2\\
\hline
Three Legs & & \\
\hline
Ring & & \\
\hline
Bidirectional Ring & & \\
\hline
Ring of Rings & & \\
\hline
\end{tabular}
\end{center}
\caption{Sizes of backprojections and of the largest factor generated by variable elimination for the SysAdmin problem with single basis functions.} \label{table:SysAdmin treewidth}
\end{table}




\section{Approximate Value Iteration with Max-norm Projection}\label{sec:AVI}

%CGjairf1 added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Factored Approximate value iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\begin{pseudocode}
FactoredAVI ($P$, $R$, $\discount$, $\Basis$, $\order$, $\eps$, $t_{max}$) \\ 
\tab \codecomment{$P$ is the factored transition model.}\\
\tab \codecomment{$R$ is the set of factored reward functions.}\\
\tab \codecomment{$\discount$ is the discount factor.}\\
\tab \codecomment{$\Basis$ is the set of basis functions $\Basis = \{\basis_1, \ldots, \basis_k\}$.}\\
\tab \codecomment{$\order$ stores the elimination order. } \\
\tab \codecomment{$\eps$ Bellman error precision. } \\
\tab \codecomment{$t_{max}$ maximum number of iterations. } \\
\tab \codecomment{Return the basis function weights $\bw$ computed by approximate value iteration.}
\bcsone
\item \codecomment{Initialize weights}
\item {\bf Let} $\bw\at{0} = {\bf 0}$.
\item \codecomment{Cache the backprojections of the basis functions.}
\item {\bf For} each basis function $\basis_i \in \Basis$; for each action $a$:
\bcstwo
\item {\bf Let} $\backbasis_i^a = \bprojfun{a}{\basis_i}$.
\ecstwo
\item \codecomment{Main approximate value iteration loop.}
\item {\bf Let} $t = 0$.
\item {\bf Repeat}
\bcstwo
\item \codecomment{Compute decision list policy for iteration $t$ weights.}
\item {\bf Let} $\dlist\at{t} = \mbox{DecisionListPolicy}(R^a + \discount \sum_i w_i\at{t} \backbasis_i^a )$.
\item \codecomment{Initialize constraints for max-norm projection LP.}
\item {\bf Let} $\constraints^{+} = \{\}$ and $\constraints^{-} = \{\}$.
\item \codecomment{Initialize indicators. }
\item {\bf Let} $\indicator = \{\}$.
\item \codecomment{For every branch of the decision list policy, generate the relevant set of constraints, and update the indicators to constraint the state space for future branches.}
\item {\bf For} each branch $\langle\event_j,a_j\rangle$ in the decision list policy $\dlist\at{t}$:
\bcsthree
\item \codecomment{Instantiate the variables in $\Event_j$ to the assignment given in $\event_j$.}
\item {\bf Instantiate} the basis functions $\Basis$ with the partial state assignment $\event_j$ and store in $C$.
\item {\bf Instantiate} the target functions $R^{a_j} + \discount\sum_i w_i\at{t} \backbasis_{a_j}$ with the partial state assignment $\event_j$ and store in $\vecb$.
\item {\bf Instantiate} the indicator functions $\indicator$ with the partial state assignment $\event_j$ and store in $\indicator'$.
\item \codecomment{Generate the factored LP constraints for the current decision list branch.}
\item {\bf Let} $\constraints^{+} = \constraints^{+} \cup  \mbox{FactoredLP}(C, - \vecb + \indicator', \order)$. 
\item {\bf Let} $\constraints^{-} = \constraints^{-} \cup  \mbox{FactoredLP}(-C, \vecb + \indicator', \order)$. 
\item \codecomment{Update the indicator functions.}
\item {\bf Let} $\indicator_j(\bx) = - \infty \indicone(\bx = \event_j)$ and update the indicators $\indicator = \indicator \cup \indicator_j$. 
\ecsthree
\item \codecomment{We can now obtain the new set of weights by solving an LP, which corresponds to the max-norm projection.}
\item {\bf Let} $\bw\at{t+1}$ be the solution of the linear program: minimize $\phi$, subject to the constraints $\{\constraints^{+} , \constraints^{-}\}$.
\item {\bf Let} $t = t+1$.
\ecstwo
\item {\bf Until} $\Bellman(\Basismatrix \bw\at{t}) \leq \eps$ or $t \geq t_{max}$ or $\bw\at{t-1} = \bw\at{t}$.
\item {\bf Return} $\bw\at{t}$.
\ecsone
\end{pseudocode}
\caption{Factored approximate value iteration with max-norm projection algorithm. }\label{code:factored AVI}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The efficient basic operations described in \secref{sec: basic
%RP3: key -> central
ops} are central to applying our max-norm solution algorithms in the
context of factored MDPs. Let us begin by considering the
approximate value iteration algorithm.  As described in
\secref{sec: Approximate Value Iteration}, the algorithm
%RP3: removed commas fore reviewer
repeatedly applies two steps: first a Bellman operator $\Tbell$;
then a max-norm projection $\Piinf$.

First consider the application of the Bellman operator: \beqn
\label{eqn: vbar = bellman} \Vbar\at{t+1} = \Tbell \Basismatrix\bw\at{t}.
\eeqn Let $\policy\at{t}$ be the stationary policy $\policy\at{t}
= \greedy{\Basismatrix\bw\at{t}}$. Note that $\policy\at{t}$ corresponds to
the maximizing policy used in the Bellman operator at iteration
$t$, \ie, $\T_{\policy\at{t}} = \Tbell$ for iteration $t$. Thus,
we can compute $\Vbar\at{t+1} = \Tbell \Basismatrix\bw\at{t}$ by first
computing $\policy\at{t}$ and then performing a
\emph{backprojection} operation for this fixed policy, \ie,
$\Vbar\at{t+1} = R_{\policy\at{t}} + \discount
P_{\policy\at{t}}\Basismatrix\bw\at{t}$.

Assume, for the moment, that $P_{\policy\at{t}}$ is a factored
transition model and that the rewards $R_{\policy\at{t}}$ are also
factored; we discuss the computation of $\Vbar\at{t+1}$ with
the policy $\policy\at{t}$ represented as a decision-list below.
As discussed in \secref{sec: one-step lookahead}, the
backprojection operation can be performed efficiently if the
transition model and the value function are both factored
appropriately.  Furthermore, the resulting function
$\Vbar\at{t+1}$ is also factored, although the factors involve
larger domains than those in the value function estimate $A
\bw\at{t}$. More specifically, denote the backprojection of basis
function $\basis_i$ through $P_{\policy\at{t}}$ as: \bdm
\backbasis^{\policy\at{t}}_i(\bx) = \sum_{\bx'} P(\bx'\mid \bx, a
= \policy\at{t}(\bx)) \basis_i(\bx'). \edm Using this definition,
we can rewrite the first step of the algorithm as: \bdm
\Vbar\at{t+1}(\bx) = R_{\policy\at{t}}(\bx) + \sum_i w_i\at{t} \
\backbasis^{\policy\at{t}}_i(\bx). \edm Therefore, $\Vbar\at{t+1}$
is composed of the sum of two factored functions: the reward
function $R_{\policy\at{t}}$ which is assumed to be factored, and
the backprojected basis functions $\backbasis^{\policy\at{t}}_i$
which are also factored, as we have just shown.

Finally, in the second step in approximate value iteration we
compute the max-norm projection $\Basismatrix\bw\at{t+1} = \Piinf
\Vbar\at{t+1}$, \ie, find $\bw\at{t+1}$ that minimizes
$\maxnorm{\Basismatrix \bw\at{t+1} - \Vbar\at{t+1}}$. This is exactly an
instance of the optimization in \eqnref{eqn:minimax},
where $C = A$ and $\vecb = \Vbar\at{t+1}$. As both $\Vbar\at{t+1}$
and $\Basismatrix\bw\at{t+1}$ are factored, we can perform this computation
using the factored LP discussed in \secref{sec: Factored
Max-norm}.

In our discussion so far, we assumed that the transition model
$P_{\policy\at{t}}$ is factored. As shown in \secref{sec:
Computing greedy policies}, the policy $\policy\at{t}$ can be
represented compactly as a decision list, \ie, it can be written
in the form $\langle\event_1,a_1\rangle,
\langle\event_2,a_2\rangle,\ldots, \langle \event_L,a_L\rangle$,
where each $\event_i$ is an assignment of values to some small
subset $\Event_i$ of variables, and each $a_i$ is an action.  The
optimal action to take in state $\bx$ is the action $a_j$
corresponding to the first event $ \event_j$ in the list with
which $\bx$ is consistent.

Although the policy $\policy\at{t}$ can be represented by a
compact decision list, the resulting transition model is usually
not factored.  Thus, we cannot simply apply our factored LP
construction, as suggested in \secref{sec: Factored LP}.  However,
through a simple change, 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 $\langle \event_i,a_i\rangle$. The
$i$th cost network only considers a subset of the states that
is consistent with the $i$th branch of the decision list. Let
$S_i$ be the set of states $\bx$ such that $\event_i$ is the first
event in the decision list for which $\bx$ is consistent. That is,
for each state $\bx \in S_i$, $\bx$ is consistent with $\event_i$,
but it is \emph{not} consistent with any $\event_j$ with $j < i$.

Recall that, as in \eqnref{eq: general constraint form}, our LP
construction defines a set of constraints that imply that $\phi
\geq \sum_i w_i \  c_i(\bx) - b(\bx)$ for each state $\bx$.
Instead, we 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_{\policy\at{t}}$.
Similarly, the reward function becomes $R^{a_i}(\bx) +
\sum_{i=1}^r R_i(\bx)$ for this subset of states.


%CGff1 more details here
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 $\event_i$, but not to states that are
consistent with some $\event_j$ for $j < i$.  To guarantee the
first condition, we simply instantiate the variables in $\Event_i$
to take the values specified in $\event_i$.  That is, our cost
network now considers only the variables in
$\{X_1,\ldots,X_n\}-\Event_i$, and computes the maximum only over
the states consistent with $\Event_i=\event_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 $\indicator_j$ for each previous decision
$\event_j$, with weight $-\infty$. More specifically,
$\indicator_j$ is a function that takes value $-\infty$ for states
consistent with $\event_j$ and zero for other all assignments of
$\Event_j$. The constraints for the $i$th branch will be of the
form:
\beqn\label{eq: constraints with decision branch indicators}
\phi \geq R(\bx,a_i)+\sum_l w_l\left(\discount
\backbasis_l(\bx,a_i) - \basis(\bx)\right) + \sum_{j<i}
-\infty\indicone(\bx=\event_j), \tab \forall
\consist{\bx}{\event_i},
\eeqn
where $\consist{\bx}{\event_i}$ defines the assignments of $\bX$
consistent with $\event_i$. The introduction of these indicators
causes the constraints associated with $\event_i$ to be trivially
satisfied by states in $S_j$ for $j<i$. Note that each of these
indicators is a restricted-scope function of $\Event_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.  The complete approximate
policy iteration with max-norm projection algorithm is outlined in
\figref{code:factored API}.



\commentout{
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 $\event_i$, but not to states that are
consistent with some $\event_j$ for $j < i$.  To guarantee the
first condition, we simply instantiate the variables in $\Event_i$
to take the values specified in $\event_i$.  That is, our cost
network now considers only the variables in
$\{X_1,\ldots,X_n\}-\Event_i$, and computes the maximum only over
the states consistent with $\Event_i=\event_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 $\indicator_j$ for each previous decision
$\event_j$, with weight $-\infty$. More specifically,
$\indicator_j$ is a function that takes value $-\infty$ for
states consistent with $\event_j$ and zero for other all
assignments of $\Event_j$. The introduction of these indicators
causes the constraints associated with $\event_i$ to be
trivially satisfied by states in $S_j$ for $j<i$. Note that each
%RP3: scope
of these indicators is a restricted scope function of $\Event_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. 
The complete approximate value iteration with max-norm projection algorithm is outlined in \figref{code:factored AVI}.
}

\commentout{
\subsubsection{Approximate Value Iteration}\label{sec: Approximate Value
Iteration}
%RP1:  Changed to start with a more generic description first
\commentout{The basic idea of approximate value iteration is quite simple.  We define an
$\Linf$ projection operator $\Piinf$ that takes a value function $\V$ and
finds $\bw$ that minimizes $\maxnorm{\Basismatrix \bw - \V}$. This is an instance of
\eqnref{eqn:minimax} and can be solved by the linear program in~(\ref{eqn:max-norm linear-program}).
Approximate value iteration algorithm alternates
applications of the Bellman operator $\Tbell$ and projection steps
$\Piinf$:
\beqa
  \Vbar\at{t+1}  & = & \Tbell \Basismatrix\bw\at{t}; \\
  \Basismatrix \bw\at{t+1}  & = & \Piinf \Vbar\at{t+1}.
\end{eqnarray*}
}

The basic idea of approximate value iteration is quite simple.  We define
a value function projection operator, $\Proj$,
that takes a value function $\V$ and
finds $\bw$ that minimizes $\norm{\Basismatrix \bw - \V}$, for some norm
$\norm{\cdot}$.
Approximate value iteration algorithms alternate
applications of the Bellman operator $\Tbell$ and projection steps $\Proj$:
\beqa
  \Vbar\at{t+1}  & = & \Tbell \Basismatrix\bw\at{t}; \\
  \Basismatrix \bw\at{t+1}  & = & \Proj \Vbar\at{t+1}.
\end{eqnarray*}

\commentout{
In standard value iteration, we only need to perform the first
step. However, $\Vbar\at{t+1}$ may not be in $\Basisspace$, 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 $\Basismatrix\bw\at{t}$ and the optimal value function $\Vstar$. This
analysis shows that the overall error depends on the single-step max-norm
projection errors $\phi\at{t+1} = \maxnorm{\Basismatrix\bw\at{t+1} -
\Vbar\at{t+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 show, this cannot happen in
approximate policy iteration. 
%CG3 need space
% , but the max-norm projection makes it as small as possible.
%We omit further details for lack of space.
}

%RP1: Reworked to bring max norm projection back in.
\commentout{
In standard value iteration, we only need to perform the first
step. However, $\Vbar\at{t+1}$ may not be representable by our
basis functions (\ie, not in $\Basisspace$), 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 $\Basismatrix\bw\at{t}$ and the optimal value
function $\Vstar$. This analysis shows that the overall error
depends on the single-step \emph{max-norm} projection errors
$\phi\at{t} = \maxnorm{\Basismatrix\bw\at{t} - \Vbar\at{t}}$. This is exactly
the quantity that is minimized in the max-norm projection. Thus,
using the max-norm projection, we can minimize these projection
errors  directly.
}


In standard value iteration, we only need to perform the first
step. However, $\Vbar\at{t+1}$ may not be representable by our basis
functions (\ie, not in $\Basisspace$), 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 $\Basismatrix\bw\at{t}$ and
the optimal value function $\Vstar$. This analysis shows that the overall
error depends on the single-step \emph{max-norm} errors $\phi\at{t} =
\maxnorm{\Basismatrix\bw\at{t} - \Vbar\at{t}}$. This is exactly the quantity that is
minimized in a max-norm projection, which is an instance of
\eqnref{eqn:minimax} and can be solved by the linear program
in~(\ref{eqn:max-norm linear-program}).  Thus, using the max-norm
projection, we can minimize the terms in our error bounds directly.

Our goal is to analyze approximate value iteration in terms of the
amount of error $\phi\at{t}$ introduced at each step by the
projection operation. If the error is zero, then we are back to
value iteration and the contraction property dictates that the
%result is $\Vstar$. If the error is small, we should get an
final result is $\Vstar$. If the error is small, we should get an
approximation that is accurate.  This result follows from the
analysis below.

Due to the contraction property of the Bellman operator,
as we see in the theorem, the overall accumulated error is a decaying average of the projection error incurred throughout all iterations:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: discounted projection error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{defn}
The \emph{discounted projection error} at iteration $t$ is defined as:
$\phibar\at{t} = \phi\at{t} + \discount \phibar\at{t-1}; \ \
\phibar\at{0} = 0.$
\qed
\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Notice that this error is bounded if $\epsp$ is fixed.
%Intuitively, at every iteration the projection step introduces at
%most $\epsp$ error and every time we apply the Bellman operator
%the error contracts. Then we can say that the total error
%introduced by projections is at most: $\epsp + \discount \epsp +
%\discount^2\epsp + \discount^3\epsp + \dots \leq
%\frac{\epsp}{1-\discount}$. As we see below, we can obtain an
%even tighter bound on the error.
We can now analyze the convergence of our approximate value
iteration procedure using the discounted projection error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theorem: error bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thm}\label{thm: vi error bound}
In the approximate value iteration algorithm, the distance between
our approximate value function at iteration $t$ and the optimal
value function $\Vstar$ is bounded by:
\begin{eqnarray}
\maxnorm{\Vstar - \Basismatrix \bw\at{t}}   & \leq &
   \discount^t \maxnorm{\Vstar - \Basismatrix\bw\at{0}} + \phibar\at{t}; \label{eq: vi error bound}
\end{eqnarray}
If $\epsp = \sup_{i} \phi\at{i}$, we have that
\beqa
\phibar\at{t}  & \leq & \frac{\epsp (1-\discount^t)}{1-\discount}.
\end{eqnarray*}
\end{thm}
%RP1:  Changed b/c this was not a proper induction proof
%\proof By induction on $t$:
\proof First we establish the following recurrence:
\beqa
\maxnorm{\Vstar - \Basismatrix \bw\at{t}}  & \leq &
         \maxnorm{\Vstar - \Tbell \Basismatrix \bw\at{t-1}} +
         \maxnorm{\Tbell \Basismatrix \bw\at{t-1} - \Basismatrix \bw\at{t}}
                    \mbox{\small~~~;~(triangle inequality;)} \\
     & \leq &  \discount \maxnorm{\Vstar - \Basismatrix \bw\at{t-1}} +  \maxnorm{\Tbell \Basismatrix\bw\at{t-1} - \Basismatrix \bw\at{t}}
                    \mbox{\small~~~;~($\Tbell$ is a contraction;)}\\
     & = &  \discount \maxnorm{\Vstar - \Basismatrix \bw\at{t-1}} + \phi\at{t}; \\
&\leq&   \discount \maxnorm{\Vstar - \Basismatrix \bw\at{t-1}} + \epsp. %\mbox{\small~~~bounded $\Linf$ approximation}
\end{eqnarray*}
>From this, the bound is easily established by induction on $t$.\qedmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thus, the distance between our approximation and the optimal value function is bounded by the sum of two terms as shown in \eqnref{eq: vi error bound}. The first, $\discount^t \maxnorm{\Basismatrix \bw\at{0} - \Vstar}$, depends on the error of our initial estimate and goes to zero exponentially quickly. This term is also present in the bounds for the standard value iteration algorithm.

The second term, $\phibar\at{t}$, is the component of the error that was
introduced by the successive projection operations, that is, by our
approximation algorithm.  Note that this second term might not be bounded,
and the maximum projection error $\epsp$ may be infinite.  Indeed, there
are examples where the error $\phi\at{t}$ grows unboundedly at each
iteration. This problem is characteristic of approximate value iteration,
also occurring when $\cL_2$ projections are applied in place of our $\Linf$
projection.%
\footnote{The same two-state MDP example of Bertsekas and
Tsitsiklis~\citeyear{Bertsekas+Tsitsiklis:1996} which causes divergence
for value iteration with unweighted $\cL_2$ projection causes divergence for
$\Linf$ projection.}
%RP3:  Toned down  However, by using our $\Linf$ projection, we minimize
%$\phi\at{t}$ directly, thereby obtaining the best possible
%bounds.
The consolation in using the $\Linf$ projection here is that we appear to
have done the best we could hope to do given our representation and the
knowledge available at each iteration.
%RP3 Furthermore, at every iteration, our algorithm returns the value of
Moreover, we obtain a useful diagnostic since the linear program returns
the $\Linf$ projection error $\phi\at{t}$, allowing us to compute directly
a bound on the amount of error introduced.
%RP3:  Changed because it sounded like this was a special property of
%$\Linf$ with policy iteration,  when it really isn't
% In the next section, we show
%that this unbounded growth in error cannot occur when we apply our $\Linf$
%projection in the policy iteration case.
Since unbounded growth of error cannot occur with approximate policy
iteration algorithms, we next analyze policy iteration with $\Linf$
projection to obtain somewhat more compelling bounds.



\commentout{
%RP1:  Changed 'determine' to 'analyze'
Our goal is to analyze the result of this process in terms of the
amount of error introduced at each step by the projection operation.
%If the error is zero, then the contraction property
%dictates that the result is $\Vstar$. If the error is small, we should
%get an approximation that is accurate.  This result follows from the
%analysis below.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: Linf projection error at t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{defn}
%The \emph{$\Linf$ projection error} at iteration $t$ is: \\
%for projecting $\Tbell \Basismatrix\bw\at{t-1}$ using $\Piinf$ %as:

We define
$\phi\at{t} = \maxnorm{\Piinf \Tbell \Basismatrix\bw\at{t-1} - \Tbell \Basismatrix\bw\at{t-1}}
 = \maxnorm{\Basismatrix \bw\at{t} - \Tbell \Basismatrix \bw\at{t-1}}$.
%\end{eqnarray*}
%\end{defn}
Due to the contraction property of the Bellman operator,
%as we will see in the theorem,
the overall accumulated error is a decaying average of the projection error
incurred throughout all iterations:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition: discounted projection error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{defn}
%The \emph{discounted projection error} at iteration $t$ is defined
%as:
$\phibar\at{t} = \phi\at{t} + \discount \phibar\at{t-1}; \ \
\phibar\at{0} = 0.$
%\qed
%\end{defn}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Notice that this error is bounded if $\epsp$ is fixed.
%Intuitively, at every iteration the projection step introduces at
%most $\epsp$ error and every time we apply the Bellman operator
%the error contracts. Then we can say that the total error
%introduced by projections is at most: $\epsp + \discount \epsp +
%\discount^2\epsp + \discount^3\epsp + \dots \leq
%\frac{\epsp}{1-\discount}$. As we will see below, we can obtain an
%even tighter bound on the error.
%We can now analyze the convergence of our approximate value iteration
%procedure using the discounted projection error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Theorem: error bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thm}\label{thm: error bound}
The distance between our approximate value function at iteration $t$
and the optimal value function $\Vstar$ is bounded by:
\beqa
\maxnorm{\Vstar - \Basismatrix \bw\at{t}}   & \leq &
   \discount^t \maxnorm{\Vstar - \Basismatrix\bw\at{0}} + \phibar\at{t};
\end{eqnarray*}
If $\epsp = \sup_{t} \phi\at{t}$, we have that
\beqa
\phibar\at{t}  & \leq & \frac{\epsp (1-\discount^t)}{1-\discount}.
\end{eqnarray*}
\end{thm}
%\proofsketch By the contraction of $\Tbell$ and induction on
%$t$.\qed

 \commentout{ \proof By induction on $t$:
\beqa
\maxnorm{\Basismatrix \bw\at{t} - \Vstar}  & \leq &
         \maxnorm{\Basismatrix \bw\at{t} - \Tbell \Basismatrix \bw\at{t-1}} +
         \maxnorm{\Tbell \Basismatrix \bw\at{t-1} - \Vstar}
                    \mbox{\small~~~triangle inequality;} \\
     & \leq &  \maxnorm{\Basismatrix \bw\at{t} - \Tbell \Basismatrix\bw\at{t-1}} +
               \discount \maxnorm{\Basismatrix \bw\at{t-1} - \Vstar}
                    \mbox{\small~~~$\Tbell$ is a contraction;}\\
     & = &  \phi\at{t} + \discount \maxnorm{\Basismatrix \bw\at{t-1} - \Vstar} \leq   \epsp + \discount \maxnorm{\Basismatrix \bw\at{t-1} - \Vstar}. %\mbox{\small~~~bounded $\Linf$ approximation}
                         \qedmath
\end{eqnarray*}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thus, the distance between our approximation and the optimal value
function is bounded by the sum of two terms. The first, $\discount^t
\maxnorm{\Basismatrix \bw\at{0} - \Vstar}$, depends on the error of our initial
estimate and goes to zero exponentially quickly. This term is also
present in the bounds for the standard value iteration algorithm
started at $\V_0 = \Basismatrix\bw\at{0}$.  The second term, $\phibar\at{t}$, is
the component of the error that was introduced by the successive
projection operations, that is, by our approximation.  Note that
this second term might not be bounded, and the maximum projection
error $\epsp$ may be infinite.  Indeed, there are examples where the error
$\phi\at{t}$ grows unboundedly at each iteration. This problem is
characteristic of approximate value iteration, also occurring when $\cL_2$
projections are applied in place of our $\Linf$ projection. However, by
using our $\Linf$ projection, we minimize $\phi\at{t}$ directly, thereby
obtaining the best possible bounds.
Furthermore, at every iteration, our
algorithm returns the value of the $\Linf$ projection error $\phi\at{t}$,
allowing us to compute directly the amount of error introduced.

If $\epsp$ is bounded, we can guarantee that our algorithm ``converges'' in
a finite number of iterations, in the sense that it eventually reaches an
$\eps$-sized $\Linf$ ball around the solution and does not leave it.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Corollary: finite iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{cor}\label{cor: finite iterations}
%For any $\eps > \epsp/(1-\discount)$, the approximate iterative fixed
%point algorithm will reach a solution $\bwstar$, with
%$\maxnorm{\Basismatrix \bwstar - \Vstar} \leq \eps$, in a finite number of
%iterations. \qed
There exists an $\eps$ and a $t^*$ such that for all $t > t^*$,
$\maxnorm{\Basismatrix \bw\at{t} - \Vstar} \leq \eps$.
\end{cor}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
} %end commentout

\commentout{
We can now fully specify the approximate value iteration algorithm.  The
algorithm iterates through the following steps:
\benum
  \item {\bf Dynamic programming step:} compute $\Vbar\at{t} = \Tbell
  \Basismatrix\bw\at{t-1}$;
  \item {\bf Stopping criterion: } if $\Bellman(\Basismatrix\bw\at{t-1}) \leq \eps$ return $\Basismatrix\bw\at{t-1}$;
  \item {\bf Projection step:} $\Basismatrix \bw\at{t} = \Piinf
  \Vbar\at{t}$;
  \item {\bf Update discounted projection error:} $\phibar\at{t} =
  \phi\at{t} + \discount \phibar\at{t-1}$.
\eenum A feature of our algorithm is that we are able to
maintain $\phibar\at{t}$, the discounted projection error incurred
throughout the execution of the algorithm. As we saw from
\thmref{thm: error bound}, this quantity can be used to bound the
error introduced by successive projection operations, that is, the
error introduced by approximating our value functions $\V\at{t}$
as a linear combination of basis functions, as opposed to
maintaining them exactly.  This estimate can be used, for example,
to evaluate how closely our chosen basis functions can approximate
the optimal value function, perhaps indicating the need to modify
our basis. It can also be compared to our stopping criterion to
evaluate if our choice of $\eps$ was too optimistic or too
pessimistic.
}
} %end commentout


In this paper, we presented new algorithms for approximate linear
programming and approximate dynamic programming (value and policy
iteration) for factored MDPs.  Both of these algorithms leverage on a novel
LP decomposition technique, analogous to variable elimination in cost
networks, which reduces an exponentially large LP to a provably
equivalent, polynomial-sized one.

%RP5: minor tweaks
Unlike previous approaches, our algorithms can exploit both additive and
context-specific structure in the factored MDP model. Typical real-world
systems possess both of these types of structure. thus, this feature of our
algorithms will increase the applicability of factored MDPs to more
practical problems.

Our approximate linear programming algorithm for factored MDPs is simpler,
easier to implement and more general than the dynamic programming
approaches. Unlike our policy iteration algorithm, it does not make any
assumptions about the structure of the actions because it does need to
represent a policy explicitly. Although this algorithm does not have the
same theoretical guarantees as max-norm projection approaches, empirically
it seems to be a favorable option. In terms of quality of the resulting
policies, our experiments suggest that approximate policy iteration tends
to generate better policies for the same set of basis
functions. Nonetheless, due to the computational advantages, we can add
more basis functions to the approximate linear programming algorithm,
obtaining a better policy and still maintaining a much faster running time
than approximate policy iteration.

Our approximate dynamic programming algorithms are motivated by error
analyses showing the importance of minimizing $\Linf$ error.  We have shown
that the $\Linf$ error can be minimized at each step using linear
programming, and provided an approach for representing this LP compactly
for factored MDPs, allowing these algorithms to be applied efficiently even
%RP5:  exponentially large -> enormous (exponential in what?)
for MDPs with enormous state spaces.  Our algorithms are more
efficient and substantially easier to implement than previous algorithms
based on the $\cL_2$-projection. Unlike our approximate linear programming
algorithm, these approximate dynamic programming algorithms require the
default action assumption, which states that actions only affect a small
number of state variables. If this assumption holds, we can also compute a
bound on the quality of the policies generated by any of our algorithms
using an efficient factored formulation of the Bellman error.


%CGjairf4 added this paragraph
We have presented experimental results on structured problems, 
highlighting three aspects of our approach: First, our factored approximate linear programming approach seems computationally advantageous over our approximate policy iteration algorithm. Although both algorithms run in polynomial time in the number of variables, the simpler approximate LP method is significantly faster. Second, we demonstrated that exploiting context-specific independence, by using a rule-based representation instead of the standard table-based one, can yield exponential improvements in computational time when the problem has significant amounts of CSI. However the overhead of managing sets of rules make it less well-suited for simpler problems. 
Finally, we compared our approach to the work of Boutilier
\etal~\citeyear{Boutilier+al:2000}, illustrating some of the tradeoffs
between the two methods. In particular, for problems with significant
context-specific structure in the value function, their approach can
be faster due to their efficient handling of the ADD
representation. However, there are problems with significant
context-specific structure in the problem representation, rather than
in the value function, which require exponentially large ADDs. In some
such problems, we demonstrated that by using a linear value function
our algorithm can obtain a polynomial-time near-optimal approximation
of the true value function. 

The success of our algorithm depends on our ability to model the problem as
a factored MDP and 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 algorithms currently require the designer to specify the
factored basis functions.  This is a limitation compared to the
algorithms of Boutilier \etal~\citeyear{Boutilier+al:2000}, which are fully automated.
However, our experiments suggest that a few simple rules can be
quite successful for designing a basis. First, we ensure that
the reward function is representable by our basis.
A simple basis that, in addition, contained a separate set of
indicators for each variable often did quite well.
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 functions
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.

Nonetheless, a general algorithm for choosing good factored basis
functions still does not exist. However, there are some potential
approaches: First, in problems with CSI, one could apply the
algorithms of Boutilier \etal for a few iterations to
generate partial tree-structured solutions. These trees could, in
turn, be used to generate a basis set for such problems. Second,
the Bellman error computation, which can be performed efficiently
as shown in \secref{sec: Bellman error}, does not only provide a
bound on the quality of the policy, but also the actual state
where the error is largest. This knowledge can be used to create a
mechanism to incrementally increase the basis set, adding new
basis functions to tackle states with high Bellman error.

There are many other possible extensions to this work. We have
already pursued extensions to collaborative multiagent systems,
where multiple agents act simultaneously to maximize the global
reward~\cite{Guestrin+al:2001c}, and factored POMDPs, where the
full state is not observed directly, but indirectly through
observation variables \cite{Guestrin+al:2001d}. However, there are
other settings that remain to be explored. In particular, we hope
to address the problem of learning a factored MDP and planning in
a competitive multiagent system. 

Additionally, in this paper we have tackled problems where the induced
width of the cost network is sufficiently low or that possess sufficient
context-specific structure to allow for the exact solution of our factored
LPs. Unfortunately, some practical problems may have prohibitively large
induced width. We plan to leverage on ideas from loopy belief propagation
algorithms for approximate inference in Bayesian
networks~\cite{Pearl:1988,Yedidia+al:2000} to address this issue.

We believe that the
methods described herein significantly further extend the
efficiency, applicability and general usability of factored models and value
functions for the control of practical dynamic systems.
