grlp.tex

%!TEX root =  autocontgrlp.tex
\begin{comment}
\section{Approximate Linear Programming: Successes and Challenges}
When the MDP has a large number of states, it is difficult to solve for $J^*$ using either the linear program \eqref{mdplp} or other full state representation methods such as value iteration or policy iteration \cite{BertB}. A practical solution is to resort to function approximation. Linear function approximation, wherein the solution is searched in the subspace spanned by the column vectors of a given feature matrix $\Phi$.
The approximate linear program (ALP) is obtained by making use of LFA in the LP, i.e., by introducing the new variables $r\in \Re^k$ and adding the extra constraint $J=\Phi r$ in \eqref{mdplp} with $\Phi \in \Re^{n\times k}$ \citep{SchSei85}.
By substitution, this leads to
\begin{align}\label{alp}
\begin{split}
\min_{r\in \Re^k}\, &c^\top \Phi r\,\,\,
\text{s.t.}\mb \Phi r\geq T \Phi r,
\end{split}
\end{align}
where $J\geq TJ$ is a shorthand for the $nd$ constraints in \eqref{mdplp} and $\Phi$ is a feature matrix whose first column is $\one$. Unless specified otherwise we use $\tr$ to denote an arbitrary solution to ALP, and we let $\tj=\Phi \tr$ to denote the corresponding approximate value function and $\tu$ to denote the greedy policy w.r.t. $\tj$.
The following is a preliminary error bound for ALP from \cite{ALP}:
\begin{theorem}[Error Bound for ALP]
\begin{align*}
\norm{J^*-\tj}_{1,c}\leq \frac{2}{1-\alpha}\min_{r}\norm{J^*-\Phi r}_\infty
\end{align*}
\end{theorem}
%For a more detailed treatment of ALP and sophisticated bounds, the reader is referred to \cite{ALP}.
%\subsection{Approximating the Constraints}
The ALP is a linear program in $k$ ($<<n$) variables as opposed to the LP in \eqref{mdplp} which has $n$ variables. Nevertheless, ALP has $nd$ constraints (same as the LP) which is an issue when $n$ is large and calls for constraint approximation/reduction techniques. Most works in literature make use of the underlying structure of the problem to cleverly reduce the number of constraints of ALP. A good example is \cite{gkp}, wherein the structure in factored linear functions is exploited. The use of basis function also helps constraint reduction in \cite{Mor-Kum}. In \cite{ALP-Bor}, the constraints are approximated indirectly by approximating the square of the Lagrange multipliers. In \cite{petrik} the transitional error is reduced ignoring the representational and sampling errors.\par
The most important work in the direction of constraint reduction is constraint sampling \cite{CS} wherein a reduced linear program (RLP) is solved instead of ALP. While the objective of RLP is same as that of ALP, RLP has only $m<<nd$ constraints \emph{sampled} from the original $nd$ constraints of ALP.  The following is a preliminary error bound for RLP from \cite{CS} holds for a special sampling distribution which is dependent on the optimal policy $u^*$ (see \cite{CS} for a detailed presentation):
\begin{theorem}[Error Bound for RLP]
Let $\mu_{u^*}\eqdef(1-\alpha)c^\top (I-\alpha P_{u^*})^{-1}$ be a probability distribution (of discounted number of visits) over the states $S$. Define $\psi_{u^*}$ to be the distribution amongst state-action pairs such that $\psi(s,a)\eqdef \frac{\mu_{u^*}}{d}$ and define $\theta\eqdef=\frac{1+\alpha}{2 c^\top J^*}\underset_{r\in \N}{\sup}\parallelJ^*-\Phi r \parallel_\infty$. Then for a given $\epsilon>0$ and $\delta>0$ it holds that
\begin{align*}
\norm{J^*-\Phi\tilde{r}_{RLP}}_{1,c}\leq \norm{J^*-\Phi\tilde{r}}_{1,c}+\epsilon \norm{J^*}_{1,c}
\end{align*}
for all $m\geq \frac{16d\theta}{(1-\alpha)\epsilon}\big(k\ln\frac{48d\theta}{(1-\alpha)\epsilon}+\ln \frac{2}{\delta}\big)$.
\end{theorem}
%A major gap in the theoretical analysis is that the error bounds are known for only a specific RLP formulated using idealized assumptions, i.e., under knowledge of $u^*$. However, RLP has be found to do well empirically in domains such as Tetris \cite{CST} and controlled queues \cite{CS}.
\subsection{Open Questions}
Interestingly, RLP has nevertheless been found to do well empirically in domain such as Tetris \cite{CST} and controlled queues \cite{CS} even when the constraints were sampled using distribution other than the ideal distribution. This fact indicates a gap in the theoretical analysis and points to the need for a more elaborate theory that addresses the issue of constraint approximation. In this paper, we answer the following questions related to constraint reduction in ALP that have so far remained open. \\
$\bullet$ As a natural generalization of RLP, what happens if we define a generalized reduced linear program (GRLP) whose constraints are positive linear combinations of the original constraints of ALP?\\
$\bullet$ Unlike \cite{CS} which provides error bounds for a specific RLP formulated using an idealized sampling distribution, is it possible to provide error bounds for any GRLP (and hence any RLP)?
In this paper, we address both of the questions above.
\end{comment}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The Linearly Relaxed ALP}\label{sec:grlp}
In this section we introduce the computational model used and the ``Linearly Relaxed Approximate Linear Program''
a relaxation of the ALP.

As discussed in the introduction, we are interested in methods that compute a good approximation to the optimal value function.
As noted earlier, at the expense of a modest additional cost, knowing an $O(\epsilon)$ approximation to $J^*$ at a few states suffices to compute actions of an $O(\epsilon)$-optimal policy. We will take a more general view, and we will consider calculating good approximations to $J^*$ with respect to a weighted $1$-norm, where the weights $c$ form a probability distribution over $\S$. Recall that the weighted $1$-norm $\norm{J}_{1,c}$ of a vector $J\in \Re^S$ is defined as $\norm{J}_{1,c}  = \sum_s c(s) |J(s)|$. Note that here and in what follows we identify elements of $\Re^\S$ (functions, mapping $\S = \{1,\dots,S\}$ to the reals) with elements of $\Re^S$ in the obvious way. This allows us to write e.g. $c^\top J$, which denotes $\sum_s c(s) J(s)$.

To introduce the optimization problem we study, first recall that 
the optimal value function $J^*$ is the solution of the fixed point equation $TJ^* = J^*$. 
It follows from the definition of $T$ that $J^* = \max_u T_u J^* \ge T_u J^*$ for any $u$, 
where $\ge$ is the componentwise partial ordering of vectors ($\le$ is the reverse relation).
With some abuse of notation, we also introduce $T_a$ to denote $T_u$ where $u(s) = a$ for any $s\in \S$. 
It follows that $J^* \ge T_a J^*$ for any $a\in \A$ and also that $T = \max_a T_a$, where again the maximization is componentwise.
We call a vector $J$ that satisfies $J \ge T_a J$ for any $a\in A$ \emph{superharmonic}. Note that this is a set of linear inequalities. 
By our note on $T$ and $(T_a)_a$, these inequalities can also be written compactly as $J \ge T J$.
\if0
Since $T$ is \emph{monotone} (i.e., for any $J_1\le J_2$ it holds that $TJ_1 \le T J_2$) and can be seen to be an $\alpha$-contraction with respect to the maximum norm, if $J$ is superharmonic then $J \ge T J \ge T^2 J \ge \dots \ge J^*$, where the last inequality follows since $T$ is an $\alpha$-contraction and hence $T^n J \to J^*$ in the maximum norm, and hence also componentwise.
Thus, $J^*$ is the ``smallest'' superharmonic function. It follows than that 
for any $c\in \Re_+^S \doteq [0,\infty)^S$ and any superharmonic $J$, $c^\top J^*\le c^\top J$ and since $J^*$ is also superharmonic, $\min \{ c^\top J\,:\, J \ge T J \} = c^\top J^*$ and if $c$ is positive-valued then the minimizer of $c^\top J$ amongst all superharminoc functions is $J^*$.
\fi
It is not hard to show then that $J^*$ is the smallest superharminoc function (i.e., for any $J$ superharmonic, $J\ge J^*$). It also follows that for any  $c\in \Re_{++}^S \doteq (0,\infty)^S$, the unique solution to the linear program $\min\{ c^\top J \,:\, J \ge T J \} =\min \{ c^\top J\,:\, J \ge T_a J, a\in \A \} $ is $J^*$. 

Now, let $\phi_1,\ldots,\phi_k : \S \to \Re$ be $k$ basis functions. 
\todoc{I think we should use $d$ to denote the number of basis functions. Or at least $n$. Low priority.} 
The \emph{Approximate Linear Program} (ALP) of \citet{SchSei85} 
is obtained by adding the linear constraints $J = \sum_{i=1}^k r_i \phi_i$ to the above linear program. Eliminating $J$ then gives $\min\{ \sum_i r_i c^\top \phi_i \,:\, \sum_i r_i \phi_i \ge g_a + \alpha \sum_i r_i P_a \phi_i, a\in \A, r = (r_i)\in \Re^k \}$.
As noted by \citet{SchSei85}, the linear program is feasible as long as $\one$, defined as the vector with all components being identically equal to one, is in the span of $\{\phi_1,\dots,\phi_k\}$. 
\emph{For the purpose of computations, it is assumed that the values $c^\top \phi_i$, $i=1,\dots, k$ and the values $(P_a \phi_i)(s)$ and $g_a(s)$ can be accessed in constant time.} 
This assumption can be relaxed to assuming that one can access $g_a(s)$ and $\phi_i(s)$ for any $(s,a)$ in constant time, as well as to that one can efficiently sample from $c$, from $P_a(s,\cdot)$ for any $(s,a)$ pair, 
but the details of this are the beyond the scope of the present work. As shown by \citet{ALP}, if $\ralp$ denotes the solution to the above ALP then for $\Jalp \doteq \sum_i \ralp(i) \phi_i \doteq \Phi \ralp$ it holds that $\norm{\Jalp - J^*}_{1,c} \le \frac{2 \epsilon}{1-\alpha}$, where $\epsilon = \inf_r \norm{ J^* - \Phi r }_\infty$ is the error of approximating the optimal value with the span of the basis functions $\phi_1,\dots,\phi_k$ and $\norm{J}_\infty = \max_s |J(s)|$ is the maximum norm and $\Phi \in \Re^{S\times k}$ is the matrix formed by $(\phi_1,\dots,\phi_k)$. That the error of approximating $J^*$ with $\Jalp$ is $O(\epsilon)$ is significant: The user can focus on finding a good basis, leaving the search for the ``right'' coefficients to a linear program solver.

While solving the ALP can be significantly cheaper than solving the LP underlying the MDP
and thus it can be advantageous for moderate-scale MDPs,
 \todoc{Give the big-Oh costs!} the number of constraints in the ALP is $SA$, 
hence the ALP is still intractable for huge-scale MDPs.
%While evaluating the objective and checking a single constraint for constraint violations can be done in constant time,
%since the number of constraints in the ALP is equal to the number of state-action pairs, 
%the ALP is still intractable in general.
 To reduce the number of constraints, we consider a relaxation of ALP 
where the constraints are replaced with positive linear 
combinations of them. 
Recalling that the constraints took the form $J \ge g_a + \alpha P_a J$ (with $J = \Phi r$),
choosing $m$ to be target number of constraints, for $1\le i \le m$, the $i$th new constraint is given by
$\sum_a w_{i,a}^\top J \ge \sum_a w_{i,a}^\top(g_a + \alpha P_a J)$, 
where the choice of $m$ and that of the vectors $w_{i,a}\in \Re_+^S$ is left to the user.
Note that this results in a linear program with $k$ variables and $m$ constraints, which can be written as
\begin{align}\label{grlp}
\begin{split}
&\underset{r\in \Re^k}{\min}\, \, c^\top \Phi r\\
&\text{s.t.}\mb  \,\sum_a W_a^\top \Phi r\geq \sum_a W_a^\top (g_a + \alpha P_a \Phi r)\,,
\end{split}
\end{align}
where $W_a = (w_{1,a},\dots,w_{m,a}) \in \Re_+^{S \times m}$. 
Note that the $(i,j)$th entry of the $m\times k$ constraint matrix of the resulting LP is
$\sum_a  w_{i,a}^\top  \phi_j - \alpha \sum_a w_{i,a}^\top P_a \phi_j$ and assuming that $(w_{i,a})_{a}$ has $p$ nonzero
elements, this can be calculated in $O( p )$ time, making the total cost of obtaining the constraint matrix to be $O(mkp)$ regardless the value of $S$ and $A$. \todoc{Somewhere discuss that the cost of obtaining a policy will depend on $A$. 
We can also discuss how $A$ could be folded into states.}

We will call the LP in \eqref{grlp} the \emph{linearly relaxed approximate linear program (LRALP)}. 
Any LP obtained using 
any constraint selection/generation process can be represented by choosing an appropriate binary-valued matrix
$W^\top = (W_1^\top,\dots,W_A^\top)\in \Re_+^{m\times SA}$. In particular, when the constraints are selected
in a random process as suggested by \citet{CS}, the matrix $W$ would be a random, binary-valued matrix.

Note that the LRALP may be unbounded. 
Unboundedness could be avoided by adding an extra constraint of the form $r\in \N$ to the LRALP,
for a properly chosen polyhedron $\N \subset \Re^k$.%
\footnote{
In particular, to obtain their theoretical result, \citet{CS} need the assumption that the set $\N$ 
is bounded and that it contains $\ralp$. In fact, 
the error bound derived by \citeauthor{CS} depends on the \emph{worst} 
error of approximating $J^*$ with $\Phi r$ when
$r$ ranges over $\N$. Hence, if $\N$ is unbounded, their bound is vacuous. In the context of a particular application, 
\citet{CS} demonstrate that $\N$ can be chosen properly to control this term. 
However, no general construction is presented to choose $\N$.}
However, it seems to us that it is downright misleading to think that guaranteeing a bounded solution 
will also lead to reasonable solutions.
Thus we will stick to the above simple form, forcing a discussion of how $W$ should be chosen to get meaningful results.%
\footnote{
The only question is whether there is some value in adding constraints beyond choosing $W$ properly.
Our position is that the set $\N$ would most likely be chosen based on very little and general information;
the useful knowledge is in choosing $W$, not in choosing some general set $\N$. 
Since randomization does not guarantee bounded solutions, \citet{ALP} must use $\N$:
In their case, $\N$ incorporates 
all the knowledge that makes the LP bounded. 
}

Further insight into the choice of $W$ can be gained by
considering the Langrangians of the ALP and LRALP. To write both LP's in a similar form let us introduce $E = (I_{S\times S},\dots,I_{S\times S})^\top$, where $I_{S\times S}$ is the $S\times S$ identity matrix. Further, let $H:\Re^S \to \Re^{SA}$ be the operator defined by $(HJ)^\top = ( (T_1 J)^\top, \dots, (T_A J)^\top )$. Note that $H$ is a linear operator, which we call the \emph{linear Bellman operator}. 
Then, the ALP can be written as $\min\{ c^\top \Phi r \,|\, E \Phi r \ge H \Phi r \}$, while LRALP takes the form  
$\min\{ c^\top \Phi r \,|\, W^\top E \Phi r \ge W^\top H \Phi r \}$. 
Hence, their Langrangians are $\mathcal{L}_{\alp}(r,\lambda) = c^\top \Phi r+\lambda^\top (H\Phi r-E\Phi r)$
$\mathcal{L}_{\lralp}(r,q) = c^\top \Phi r+q^\top W^\top (H\Phi r-E\Phi r)$. Thus, we can view $W q$ as a ``linear approximation''
to the dual variable $\lambda \in \Re_+^{SA}$. 
This suggests that perhaps $W$ should be chosen such that it approximates well the optimal dual variable.
If $\Phi$ spans $\Re^{S}$, the optimal dual variable $\lambda^*$ is known to be the discounted occupancy measure underlying the optimal policy (Theorem 3.18, \cite{Kall17}), suggesting that the role of $W$ is very similar to the role of $\Phi$ excepts that the subspace spanned by the columns of $W$ should ideally be close to $\lambda^*$.

\if0
\textbf{Linear function approximation in Primal and Dual Variables:} Let us look at the Lagrangian of ALP and GRLP in
\eqref{lag} and \eqref{lag2} respectively, i.e.,
\begin{align}\label{lag}
\tilde{L}(r,\lambda)=c^\top \Phi r+\lambda^\top (T\Phi r-\Phi r), \\ \label{lag2}\hat{L}(r,q)=c^\top \Phi r+q^\top W^\top (T\Phi r-\Phi r).
\end{align}
Thus, when $Wq = \lambda$, i.e., when $W$ is a set of basis functions that allow
a low dimensional linear representation of the dual variables $\lambda$,
the two problems are the same.
%Note that $ Wq\approx \lambda$ in \eqref{lag2} is linear function approximation of the Lagrange multipliers.
Hence, while ALP employs LFA in its objective function (i.e., use of $\Phi r$), GRLP employs linear approximation both in the objective function ($\Phi r$) as well as the constraints (use of $W$).
%Further, $W$ can be interpreted as the feature matrix that approximates the Lagrange multipliers as $\lambda\approx Wq$, where $\lambda \in \Re^{nd}, r\in \Re^m$.
To get a sense of how $W$ should be chosen, recall that
the optimal Lagrange multipliers are the discounted number of visits to the ``state-action pairs'' under an optimal policy $u^*$, i.e.,
\begin{align}
\lambda^*(s,u^*(s))&=\big(c^\top(I-\alpha P_{u^*})^{-1}\big)(s)\nn\\
&= \big(c^\top(I+\alpha P_{u^*}+\alpha^2 P_{u^*}^2+\ldots)\big)(s),\nn\\
\lambda^*(s,a)&=0, \qquad \text{for all } a \neq u^*(s),\nn
\end{align}
where $P_{u^*}$ is the probability transition matrix under $u^*$ ($P_{u^*}(s,s') = P_{u^*(s)}(s,s')$, $s,s'\in S$) \cite{dolgov}. Even though we might not have the optimal policy $u^*$ in practice, the fact that $\lambda^*$ is a probability distribution and that it is a linear combination of $\{P_{u^*},P^2_{u^*},\ldots\}$ hints at the kind of features that might be useful for constructing the $W$ matrix.
\fi