Nonparametric Approximate Dynamic Programming via the Kernel Method

1. Introduction

Problems of dynamic optimization in the face of uncertainty are frequently posed as Markov decision processes (MDPs). The central computational problem is then reduced to the computation of an optimal “value” or “cost-to-go” function that encodes the value garnered under an optimal policy starting from any given MDP state. MDPs for many problems of practical interest frequently have intractably large state spaces precluding exact computation of the cost-to-go function. Approximate dynamic programming (ADP) is an umbrella term for algorithms designed to produce good approximations to this function. Such approximations then imply a natural “greedy” control policy.

ADP algorithms are often formulated in parametric settings. Here, the user specifies an “approximation architecture” (i.e., a set of basis functions), and the algorithm then produces an approximation in the span of this basis. The strongest theoretical results available for such algorithms typically share the following two features:

• The quality of the approximation produced is comparable with the best possible within the basis (i.e., the parameterization or approximation architecture) specified.

• The sample complexity or computational effort required for these algorithms scales, typically polynomially, with the dimension of the basis.

These results highlight the importance of selecting a “good” approximation architecture and remain somewhat dissatisfying in that additional sampling or computational effort cannot remedy a bad approximation architecture.

In contrast, an ideal nonparametric approach would, in principle, free the user from carefully specifying a suitable low-dimensional approximation architecture. Instead, the user would have the liberty of selecting a very rich architecture (whose specification is potentially implicit, such as via a kernel). The quality of the approximation produced by the algorithm would then improve—gracefully—with the extent of computational or sampling effort expended, ultimately becoming exact. Although they may fall short of the ideal notion of a nonparametric ADP algorithm, there are several existing proposals for nonparametric ADP available:

1.1. Kernelizing Policy Evaluation

These are approaches that are motivated by approximate policy iteration. The key computational step in approximate policy iteration methods is “policy evaluation.” The aim of this step is to find the cost-to-go function for a given policy. This step involves solving the projected Bellman equation, a linear stochastic fixed-point equation. Because of the curse of dimensionality, solving this exactly is infeasible, and a portion of the ADP literature focuses on approximately solving the policy-evaluation step. This usually involves approximating the cost-to-go function for a particular policy with a parametric approximation architecture. Nonparametric approaches of this flavor attempt to make this step nonparametric. Bethke et al. (2023) study the problem of minimizing Bellman errors from the point of view of kernel support vector regression and Gaussian processes. Engel et al. (2003) consider a generative model for the cost-to-go function, which is based on Gaussian processes. Xu et al. (2007) introduce the nonparametric version of the Least-Squares Temporal Difference-$(\lambda )$ method, a numerically stable approach to do policy evaluation. Farahmand et al. (2009a) study regularized versions of some popular policy-evaluation methods and perform the sample-complexity analysis for one policy-evaluation step. It is difficult to provide convergence guarantees for approximate policy iteration schemes in parametric settings, and these difficulties remain in nonparametric variations. It is consequently difficult to characterize the computational effort or sample complexity required to produce a good approximation with such approaches. In practice, these theoretical shortcomings are often bypassed with heuristics, such as stopping after a fixed number of iterations. In contrast, the approach we study in this paper will admit sample complexity and approximation guarantees, but at the expense of assuming access to a certain idealized sampling distribution over states.

1.2. Local Averaging

Another approach to nonparametric ADP has been to use kernel-based local averaging to approximate the solution of an MDP with that of a simpler variation on a sampled state space (e.g., Ormoneit and Glynn 2002, Ormoneit and Sen 2002, and Barreto et al. 2011). Convergence rates for these local averaging methods are exponential in the dimension of the problem state space. As in our setting, Dietterich and Wang (2002) construct kernel-based cost-to-go function approximations. These are subsequently “plugged in” to various ad hoc optimization-based ADP formulations. This plug-in approach is inappropriate because it effectively translates to a nonparametric regression approach with no regularization whatsoever. One cannot expect generalization guarantees for the resulting function class, and, indeed, the authors are unable to provide guarantees. Practically, the absence of regularization is also likely to hurt performance, as is borne out in our experiments. Similarly, Ernst et al. (2005) replace the local averaging procedure used for regression by Ormoneit and Sen (2002) with nonparametric regression procedures, such as tree-based learning methods. This is done again without any theoretical justification. Pazis and Parr (2011) discuss a nonparametric method that explicitly restricts the smoothness of the value function. However, sample complexity results for this method are not provided, and it appears unsuitable for high-dimensional problems (such as, for instance, the queuing problem we consider in our experiments).

1.3. Feature Selection via ${\mathcal{\ell }}_1$-Penalty (Parametric)

Closely related to our work, Kolter and Ng (2009) and Petrik et al. (2010) consider modifying the approximate linear program with an ${\mathcal{\ell }}_1$-regularization term to encourage sparse approximations in the span of a large, but necessarily tractable, set of features. In contrast to this line of work, our approach will allow for approximations in a potentially full-dimensional (i.e., of dimension equal to the cardinality of the state space) approximation architecture, with a constraint on an appropriate ${\mathcal{\ell }}_2$-norm of the weight vector to provide regularization.

In addition, there are still other nonparametric proposals that do not quite fall into the three categories above. For example, Farahmand et al. (2009b) study what is effectively a kernelized version of a specific approximate value iteration algorithm. Again, it is not possible to provide performance guarantees of the type we present in this paper because doing so in the parametric case for approximate value iteration is already hard; see De Farias and Van Roy (2000).

This paper presents what we believe is a practical, nonparametric ADP algorithm that enjoys approximation and sample complexity guarantees that serve as natural counterparts to those available for linear programming (LP) approaches to ADP. In greater detail, we make the following contributions:

• A new mathematical programming formulation. We rigorously develop a kernel-based variation of the “smoothed” approximate LP (SALP) approach to ADP proposed by Desai et al. (2012). The resulting mathematical program, which we dub the regularized smoothed approximate LP (RSALP), is distinct from simply substituting the local averaging approximation above in the SALP formulation. We develop a companion active-set method that is capable of solving this mathematical program rapidly and with limited memory requirements.

• Theoretical guarantees.¹ Our algorithm can be interpreted as solving an approximate linear program in a (potentially infinite-dimensional) Hilbert space. We provide a probabilistic upper bound on the approximation error of the algorithm relative to the best possible approximation one may compute in this space, subject to a regularization constraint. We show that the sample complexity of our algorithm grows polynomially as a function of a regularization parameter. As this regularization parameter is allowed to grow, so does the set of permissible approximations, eventually permitting the best approximation within the architecture, which, for some of the architectures we consider, is exact.

  The sampling requirements for our method are independent of the dimension of the approximation architecture. Instead, they grow with the allowed complexity of the approximation. This result can be seen as the “right” generalization of the prior parametric approximate LP approaches of de Farias and Van Roy (2003, 2004) and Desai et al. (2012), where, in contrast, sample complexity grows with the dimension of the approximating architecture.

• A computational study. To study the efficacy of our approach, we consider an MDP arising from a challenging queueing network scheduling problem. We demonstrate that our method yields significant improvements over tailored heuristics and parametric ADP methods, all while using a generic high-dimensional approximation architecture. In particular, these results suggest the possibility of solving a challenging high-dimensional MDP using an entirely generic approach.

Relative to so-called “model-free” approaches to approximate DP based on simulation (such as temporal difference-learning), we note that the method we propose in the paper relies on the knowledge of the transition probabilities of the MDP. Further, we require that the number of possible next states for any state action pair to tractable. One can rely on sampling when this is not the case, but a rigorous analysis of that approach is beyond the scope of this paper; see Haskell et al. (2016) for a general analysis of so-called “empirical” Bellman operators. These issues are germane to linear programming-based ADP methods broadly.

The organization of the paper is as follows: In Section 2, we formulate an infinite-dimensional LP for our problem and present an effective approximate-solution approach. Section 3 provides theoretical guarantees for the quality of the approximations computed via our nonparametric algorithm. Theorem 1 in that section provides our main guarantee. In Section 4, we provide an active-set method that can be used to efficiently solve the required quadratic optimization problem central to our approach, while respecting practical memory constraints. We also establish the correctness of our active-set approach. Section 5 describes a numerical study for a criss-cross queueing system benchmarking our approach against approximate linear programming approaches and tailor-made heuristics. Section 6 concludes.

2. Formulation

2.1. Preliminaries

Consider a discrete-time Markov decision process with finite state space $\mathcal{S}$ and finite action space $\mathcal{A}$. We denote by x_t and a_t, respectively, the state and action at time t. For notational simplicity, and without loss of generality, we assume that all actions are permissible at any given state. We assume time-homogeneous Markovian dynamics: Conditioned on being at state x and taking action a, the system transitions to state $x^{\prime }$ with probability $p(x,x\prime ,a)$ independent of the past. A policy is a map $\mu :\mathcal{S}\to \mathcal{A}$, so that

represents the expected (discounted, infinite horizon) cost-to-go under policy μ starting at state x, with the discount factor $\alpha \in (0,1)$. Letting Π denote the set of all policies, our goal is to find an optimal policy ${\mu }^{*}$ such that ${\mu }^{*}\in {\text{arg}\hspace{0.17em}\text{max}}_{\mu \in \mathrm{\Pi }}{J}^{\mu }(x)$ for all $x\in \mathcal{S}$ (it is well known that such a policy exists). We denote the optimal cost-to-go function by $J^{*}\triangleq J^{{\mu }^{*}}$. An optimal policy ${\mu }^{*}$ can be recovered as a greedy policy with respect to $J^{*}$,where we define the expectation ${\mathsf{E}}_{x,a}[f(X\prime )]\triangleq {\sum }_{x\prime \in \mathcal{S}}p(x,x\prime ,a)f(x\prime )$, for all functions $f:\mathcal{S}\to \mathbb{R}$ on the state space.

Because in practical applications, the state space $\mathcal{S}$ is often intractably large, exact computation of $J^{*}$ is untenable. ADP algorithms are principally tasked with computing approximations to $J^{*}$ of the form $J^{*}(x)\approx z^{\top }\mathrm{\Phi }(x)\triangleq \tilde{J}(x),$ where $\mathrm{\Phi }:\mathcal{S}\to {\mathbb{R}}^m$ is referred to as an approximation architecture or a basis and must be provided as input to the ADP algorithm. The ADP algorithm computes a “weight” vector $z\in {\mathbb{R}}^m$. One then employs a policy that is greedy with respect to the corresponding approximation $\tilde{J}$.

2.2. Primal Formulation

The approach we propose is based on the LP formulation of dynamic programming. It was observed by Manne (1960) that the optimal cost-to-go $J^{*}$ can be obtained by solving the following linear program:

for any strictly positive state-relevance weight vector $\nu \in {\mathbb{R}}_{+}^{\mathcal{S}}$. Motivated by this, a series of ADP algorithms (Schweitzer and Seidman 1985, de Farias and Van Roy 2003, Desai et al. 2012) have been proposed that compute a weight vector z by solving an appropriate modification of (1). In particular, Desai et al. (2012) propose solving the following optimization problem:

Here, $\kappa >0$ is a penalty parameter, and $\pi \in {\mathbb{R}}_{+}^{\mathcal{S}}$ is a strictly positive distribution on $\mathcal{S}$. In considering the above program, notice that if one insisted that the slack variables s were precisely zero, the Program (2) would be identical to (1), with the additional restriction to value-function approximations of the form $J(x)=z^{\top }\mathrm{\Phi }(x)$. This case is known as the approximate linear program (ALP) and was first proposed by Schweitzer and Seidman (1985). de Farias and Van Roy (2003) provided a pioneering analysis that, stated loosely, showed

for an optimal solution $z^{*}$ to the ALP. Desai et al. (2012) showed that these bounds could be improved upon by smoothing the constraints of the ALP—that is, permitting positive slacks. The resulting Program (2) is called the smoothed approximate linear program. The ALP constraints impose the restriction that $TJ\ge J$, where T is the Bellman operator defined byfor all $x\in \mathcal{S}$ and $J:\mathcal{S}\to \mathbb{R}$. This constraint ensures that the optimal solution lower-bounds the optimal value function at each state. But at the cost of providing a lower bound, the ALP might do a poor job of approximating $J^{*}$ due to presence of constraints associated with certain states that are rarely visited. The SALP prevents this by allowing violations of the Bellman constrains, but weighting the violations in a way that they are strongly penalized at highly visited states. For a more detailed interpretation of SALP, see section 3 of Desai et al. (2012). In both the original ALP and the SALP, one can solve a “sampled” version of the above program for problems with large state spaces.

Now, consider allowing $\mathrm{\Phi }$ to map from $\mathcal{S}$ to a general (potentially infinite-dimensional) Hilbert space $\mathcal{H}$. We use bold letters to denote elements in the Hilbert space $\mathcal{H}$; for example, the weight vector is denoted by $\mathbf{z}\in \mathcal{H}$. We further suppress the dependence on $\mathrm{\Phi }$ and denote the elements of $\mathcal{H}$ corresponding to their counterparts in $\mathcal{S}$ by bold letters. Hence, for example, $\mathbf{x}\triangleq \mathrm{\Phi }(x)$ and $\mathbf{X}\triangleq \mathrm{\Phi }(X)$. Denote the image of the state space under the map $\mathrm{\Phi }$ by $\mathcal{X}\triangleq \mathrm{\Phi }(\mathcal{S})$; $\mathcal{X}\subset \mathcal{H}$. The analogous value-function approximation in this case would be given by

where b is a scalar offset corresponding to a constant basis function.² The following generalization of (2)—which we dub the regularized SALP—then essentially suggests itself:

The only new ingredient in the program above is the fact that we regularize the weight vector $\mathbf{z}$ using the parameter $\mathrm{\Gamma }>0$. Penalizing the objective of (5) according to the square of the norm ${\Vert \mathbf{z}\Vert }_{\mathcal{H}}\triangleq \sqrt{〈\mathbf{z},\mathbf{z}〉}$ anticipates that we will eventually resort to sampling in solving this program. In a sampled setting, regularization is necessary to avoid overfitting and, in particular, to construct an approximation that generalizes well to unsampled states. This regularization, which plays a crucial role both in theory and practice, is easily missed if one directly plugs in a local averaging approximation in place of $z^{\top }\mathrm{\Phi }(x)$, as is the case in the earlier work of Dietterich and Wang (2002), or a more general nonparametric approximation, as in the work of Ernst et al. (2005).

Because the RSALP of (5) can be interpreted as a regularized stochastic optimization problem, one may hope to solve it via its sample average approximation. To this end, define the likelihood ratio $w_x\triangleq {\nu }_x/{\pi }_x$, and let $\widehat{\mathcal{S}}\subset \mathcal{S}$ be a set of N states sampled independently, according to the distribution π. The sample average approximation of (5) is then

We call this program the sampled RSALP. Even if $|\widehat{\mathcal{S}}|$ were small, it is still not clear that this program can be solved effectively: If $\mathcal{H}$ were infinite-dimensional, this is a semi-infinite linear program. We will, in fact, solve the dual to this problem.

2.3. Dual Formulation

We begin by establishing some notation. Let ${\mathcal{N}}_{x,a}\triangleq \{x\}\cup \{x\prime \in \mathcal{S}\hspace{0.17em}:p(x,x\prime ,a)>0\}$ denote the set of states that can be reached starting at a state $x\in \mathcal{S}$ given an action $a\in \mathcal{A}$. For any states $x,x\prime \in \mathcal{S}$ and action $a\in \mathcal{A}$, define $q_{x,x\prime ,a}\triangleq {\mathbf{1}}_{\{x=x\prime \}}-\alpha p(x,x\prime ,a)$. Now, define the symmetric positive semidefinite matrix $Q\in {\mathbb{R}}^{(\widehat{\mathcal{S}}\times \mathcal{A})\times (\widehat{\mathcal{S}}\times \mathcal{A})}$ according to

the vector $R\in {\mathbb{R}}^{\widehat{\mathcal{S}}\times \mathcal{A}}$ according toand the scalar S as

Notice that Q, R, and S depend only on inner products in $\mathcal{X}$ (and other easily computable quantities). The dual to (6) is then given by:

Assuming that Q, R, and S can be easily computed, this quadratic program is tractable—its size is polynomial in the number of sampled states. We may recover a primal solution (i.e., the weight vector ${\mathbf{z}}^{*}$) from an optimal dual solution:

Proposition 1.

Programs (6) and (9) have equal (finite) optimal values. The optimal solution to (9) is attained at some ${\lambda }^{*}$. The optimal solution to (6) is attained at some $(z^{*},s^{*},b^{*})$ with

$${\mathbf{z}}^{*}=\frac{1}{\mathrm{\Gamma }}\left[\frac{1}{N}{\displaystyle \sum _{x\in \widehat{\mathcal{S}}}{w}_x}\mathbf{x}-{\displaystyle \sum _{x\in \widehat{\mathcal{S}},a\in \mathcal{A}}{\lambda }_{x,a}^{*}}(\mathbf{x}-\alpha {\mathsf{E}}_{x,a}[\mathbf{X}\prime ])\right].$$(10)

The proof of Proposition 1 is presented in Appendix A. Having solved this program, we may, using Proposition 1, recover our approximate cost-to-go function $\tilde{J}(x)=〈{\mathbf{z}}^{*},\mathbf{x}〉+b^{*}$ as³

A policy greedy with respect to $\tilde{J}$ is not affected by constant translations; hence, in (11), the value of $b^{*}$ can be set to be zero arbitrarily. Again, note that given ${\lambda }^{*}$, evaluation of $\tilde{J}$ only involves inner products in $\mathcal{X}$.

2.4. Kernels

As pointed out earlier, the sampled RSALP is potentially difficult to work with. Proposition 1 establishes that solving this program (via its dual) is a computation that scales polynomially in N, so that it can be solved efficiently, provided that inner products in $\mathcal{H}$ can be evaluated cheaply. Alternatively, we may have arrived at a similar conclusion by observing that in any optimal solution to (6), we must have that ${\mathbf{z}}^{*}\in \text{span}\hspace{0.17em}\{\mathbf{x}\hspace{0.17em}:\hspace{0.17em}x\in \widehat{\mathcal{S}}\cup \mathcal{N}(\widehat{\mathcal{S}})\}$, where $\mathcal{N}(\widehat{\mathcal{S}})$ denotes the set of states that can be reached from the sampled states of $\widehat{\mathcal{S}}$ in a single transition. Then, one can restrict the feasible region of (6) to this subspace. In either approach, we observe that one need not necessarily have explicitly characterized the feature map $\mathrm{\Phi }(·)$; knowing $〈\mathrm{\Phi }(x),\mathrm{\Phi }(y)〉$ for all $x,y\in \mathcal{S}$ would suffice, and so the algorithm designer can focus on simply specifying these inner products. We will see that when $\mathcal{S}\subset {\mathbb{R}}^n$, these inner products are effectively cheap to compute for practically interesting choices of $\mathrm{\Phi }$ and $\mathcal{H}$. This leads us to what is popularly referred to as the “kernel trick,” which we discuss next without the assumption that $\mathcal{S}$ is necessarily a finite set.

A kernel is a map $K:\mathcal{S}\times \mathcal{S}\to \mathbb{R}$; we will call such a kernel positive definite if, for any finite collection of elements ${\{x_i\}}_{1\le i\le n}$ in $\mathcal{S}$, the Gram matrix $G\in {\mathbb{R}}^{n\times n}$ defined by $G_{ij}\triangleq K(x_i,x_j)$ is symmetric and positive semidefinite. Given such a kernel, we are assured of the existence of a Hilbert space $\mathcal{H}$ and a map $\mathrm{\Phi }:\mathcal{S}\to \mathcal{H}$ such that⁴ $〈\mathrm{\Phi }(x),\mathrm{\Phi }(y)〉=K(x,y)$. The kernel trick then allows us to implicitly work with this space $\mathcal{H}$ by replacing inner products of the form $〈\mathbf{x},\mathbf{y}〉\triangleq 〈\mathrm{\Phi }(x),\mathrm{\Phi }(y)〉$ in (7), (8), and (11) by K(x, y). Of course, the quality of the approximation one may produce depends on $\mathcal{H}$ and, consequently, on the kernel employed.

One case of special interest is where $\mathcal{S}$ is finite and the Gram matrix corresponding to this set is positive definite. A kernel satisfying this condition is known as full-dimensional or full-rank. In this case, we may take $\mathcal{H}$ to be the $|\mathcal{S}|$-dimensional Euclidean space. Working with such a kernel would imply working with an approximation architecture where any function $J:\mathcal{S}\to \mathbb{R}$ can be exactly expressed in the form $J(x)=〈\mathrm{\Phi }(x),\mathbf{z}〉$, for some weight vector $\mathbf{z}\in \mathcal{H}$.

There are a number of kernels one can specify on the state space; we give two examples here for the case where $\mathcal{S}\subset {\mathbb{R}}^n$. The polynomial kernel is defined according to

A corresponding Hilbert space $\mathcal{H}$ is given by the span of all monomials of degree up to d on ${\mathbb{R}}^n$. A Gaussian radial basis kernel is defined according to

The Gaussian kernel is known to be full-dimensional (see, e.g., theorem 2.18 in Scholkopf and Smola 2001), so that employing such a kernel in our setting would correspond to working with an infinite-dimensional approximation architecture. A thorough exposition on this topic, along with many more examples, can be found in the text of Scholkopf and Smola (2001).

2.5. Overall Procedure

Our development thus far suggests the following nonparametric algorithm that requires as input a kernel, K; a state sampling distribution, π; and a distribution that allows us to weight approximation across states, ν. Recall that we defined by w the likelihood ration between ν and π; $w_x={\nu }_x/{\pi }_x$. Lastly, recall that ${\mathcal{N}}_{x,a}$ denotes the set of states reachable with positive probability in one step from state x, having taken action a. Our procedure requires that we set the parameters κ (which we recall penalizes the one sided Bellman error across states) and Γ (which regularizes the set of approximations allowed). Finally, we must specify the number of samples N.

1. Sample N states from a distribution π.

2. Solve (9) with

   $$Q(x,a,x\prime ,a\prime )\triangleq {\displaystyle \sum _{y\in {\mathcal{N}}_{x,a}}{\displaystyle \sum _{y\prime \in {\mathcal{N}}_{x\prime ,a\prime }}{q}_{x,y,a}}}{q}_{x\prime ,y\prime ,a\prime }K(y,y\prime ),$$(12)
   and
   $$R(x,a)\triangleq \left(\mathrm{\Gamma }{g}_{x,a}-\frac{1}{N}{\displaystyle \sum _{y\in \widehat{\mathcal{S}}}{\displaystyle \sum _{x\prime \in {\mathcal{N}}_x,a}{w}_y}}{q}_{x,x\prime ,a}K(y,x\prime )\right).$$(13)

   The value of S may be set to be zero.

3. Let ${\lambda }^{*}$ be the dual optimal. Define an approximate value function by

   $$\tilde{J}(x)\triangleq \frac{1}{\mathrm{\Gamma }}\left[\frac{1}{N}{\displaystyle \sum _{y\in \widehat{\mathcal{S}}}{w}_y}K(x,y)-{\displaystyle \sum _{y\in \widehat{\mathcal{S}},a\in \mathcal{A}}{\lambda }_{y,a}^{*}}(K(x,y)-\alpha {\mathsf{E}}_{y,a}[K(X\prime ,x)])\right].$$(14)

In the next section, we will develop theory to characterize the sample complexity of our procedure and the approximation quality it provides. This theory will highlight the roles of the kernel K and the sampling distributions used.

3. Approximation Guarantees

3.1. Overview

Recall that we are employing an approximation ${\tilde{J}}_{\mathbf{z},b}$ of the form

parameterized by the weight vector $\mathbf{z}$ and the offset parameter b. For the purpose of establishing theoretical guarantees, in this section, we will look at the following variation of the RSALP of (5):

Here, rather than regularizing by penalizing according to the weight vector $\mathbf{z}$ in the objective as in (5), we regularize by restricting the size of the weight vector as a constraint. This regularization constraint is parameterized by the scalar $C\ge 0$. It is easy to see that (15) is equivalent to the original Problem (5), in the following sense: For any $\mathrm{\Gamma }>0$, there exists a $C\ge 0$ so that for all B sufficiently large, (5) and (15) have the same optimal solutions and value.⁵ Let $C(\mathrm{\Gamma })$ be any such value of C, corresponding to a particular Γ.

Now, let

The best approximation to $J^{*}$ within this set has ${\mathcal{\ell }}_{\infty }$-approximation error

Now, observe that if $\text{span}\hspace{0.17em}\{\mathbf{x}\hspace{0.17em}:\hspace{0.17em}x\in \mathcal{S}\}\subset \mathcal{H}$ has dimension $|\mathcal{S}|$ (i.e., if the kernel is full-dimensional), there exists a $\tilde{\mathbf{z}}\in \mathcal{H}$ satisfying ${\tilde{J}}_{\tilde{\mathbf{z}},0}=J^{*}$. Consequently, for $C(\mathrm{\Gamma })$ sufficiently large and any value of $B\ge 0$, we have that $\tilde{\mathbf{z}}\in \mathcal{C}(C(\mathrm{\Gamma }),B)$—the approximation error in (16) reduces to zero. More generally, as $C(\mathrm{\Gamma })$ and B grow large, we expect the approximation error in (16) to monotonically decrease; the precise rate of this decrease will depend, of course, on the feature map $\mathrm{\Phi }$, which, in turn, will depend on the kernel we employ.

Loosely speaking, in this section, we will demonstrate that:

1. For any given $C(\mathrm{\Gamma })$ and B, the exact solution of the RSALP (15) will produce an approximation with ${\mathcal{\ell }}_{\infty }$-error comparable to the optimal ${\mathcal{\ell }}_{\infty }$-error possible for approximations ${\tilde{J}}_{\mathbf{z},b}$ with ${\Vert \mathbf{z}\Vert }_{\mathcal{H}}\le C(\mathrm{\Gamma })$ and $b\le B$.

2. Under a certain set of idealized assumptions, solving a sampled variation of the RSALP (15) with a number of samples $N(C(\mathrm{\Gamma }),B)$ that scales gracefully with $C(\mathrm{\Gamma })$, B, and other natural parameters produces a near-optimal solution to the RSALP with high probability. These sample complexity bounds will not depend on the dimension of the approximation architecture $\mathcal{H}$.

Because $C(\mathrm{\Gamma })$ and B are parameters of the algorithm designer’s choosing, these results will effectively show that with the ability to use a larger number of samples, the designer may choose larger values of $C(\mathrm{\Gamma })$ and B and, consequently, produce approximations of improving quality. Under mild assumptions on the kernel, the approximation error can be made arbitrarily small in this fashion.

3.2. Preliminaries and an Idealized Program

Before proceeding with our analysis, we introduce some preliminary notation and define an idealized sampling distribution. For a vector $x\in {\mathbb{R}}^n$ and a weight vector $v\in {\mathbb{R}}_{+}^n$, we denote by ${\Vert x\Vert }_v\triangleq {\sum }_{i=1}^n{v}_i|x_i|$ the weighted 1-norm of x. We next define an idealized sampling distribution. Let ν be an arbitrary distribution over $\mathcal{S}$ and denote by ${\mu }^{*}$ and $P_{{\mu }^{*}}$ an optimal policy and the transition probability matrix corresponding to this policy, respectively. We define a distribution over $\mathcal{S},\hspace{0.17em}{\pi }_{{\mu }^{*},\nu }$ according to

This idealized distribution will play the role of the sampling distribution π in the sequel.

As before, let $\widehat{\mathcal{S}}$ be a set of N states drawn independently at random from $\mathcal{S}$; here, we pick a specific sampling distribution $\pi ={\pi }_{{\mu }^{*},\nu }$. Given the definition of ${\tilde{J}}_{\mathbf{z},b}$ in (4), the “idealized” sampled program we consider is:

Here, T is the Bellman operator defined in (3).

This program is a sampled variant of (15) and is closely related to the sampled RSALP (6) introduced in the last section. Before proceeding with an analysis of the quality of approximation afforded by solving this program, we discuss its connection with the sampled RSALP (6) of the previous section:

1. The distributions ν and π: We allow for an arbitrary distribution ν, but, given this distribution, require that $\pi ={\pi }_{{\mu }^{*},\nu }$. In particular, the distribution ν might be chosen as the empirical distribution corresponding to N independent draws from $\mathcal{S}$ under some measure; one would then draw N independent samples under $\pi ={\pi }_{{\mu }^{*},\nu }$ to construct the second term in the objective. In a sense, this is the only “serious” idealized assumption we make here. Given the broad nature of the class of problems considered, it is hard to expect meaningful results without an assumption of this sort, and, indeed, much of the antecedent literature considering parametric LP-based approaches makes such an assumption. When the sampling distribution π is a close approximation to ${\pi }_{{\mu }^{*},\nu }$ (say, the likelihood ratio between the two distributions is bounded), then it is possible to provide natural extensions to our results that account for how close the sampling distribution used is to the idealized sampling distribution.

2. Regularization: We regularize the weight vector $\mathbf{z}$ with an explicit constraint on ${\Vert \mathbf{z}\Vert }_{\mathcal{H}}$; this permits a transparent analysis and is equivalent to placing the “soft” regularization term $\frac{\mathrm{\Gamma }}{2}{\Vert \mathbf{z}\Vert }_{\mathcal{H}}^2$ in the objective. In particular, the notation $C(\mathrm{\Gamma })$ makes this equivalence explicit. In addition, we place a separate upper bound on the offset term B. This constraint is not binding if $B>2C{\max }_{x\in \mathcal{S}}{\Vert \mathbf{x}\Vert }_{\mathcal{H}}+{\Vert g\Vert }_{\infty }/(1-\alpha )$, but, again, permits a transparent analysis of the dependence of the sample complexity of solving our idealized program on the offset permitted by the approximation architecture.

3. The choice of κ: The smoothing parameter κ can be interpreted as yet another regularization parameter, here on the (one-sided) Bellman error permitted for our approximation. Our idealized program chooses a specific value for this smoothing parameter, in line with that chosen by Desai et al. (2012).

3.3. The Approximation Guarantee

Let $(\widehat{\mathbf{z}},\widehat{b})$ be an optimal solution to the idealized sampled Program (18) and let $K\triangleq {\max }_{x\in \mathcal{S}}{\Vert \mathbf{x}\Vert }_{\mathcal{H}}$. Further, define

Notice that $\mathrm{\Xi }{(C,B,K,\delta )}^2$ scales as the square of C, K, and B and, further, is $O(\ln 1/\delta )$. The following result will constitute our main approximation guarantee:

Theorem 1.

For any $ϵ,\delta >0$, let $N\ge \mathrm{\Xi }{(C,B,K,\delta )}^2/{ϵ}^2$. If $(\widehat{\mathbf{z}},\widehat{b})$ is an optimal solution to (18), then with probability at least $1-\delta -{\delta }^4$,

$$\begin{array}{ll}{\Vert J^{*}-{\tilde{J}}_{\widehat{\mathbf{z}},\widehat{b}}\Vert }_{1,\nu }\hfill & \le \frac{3+\alpha }{1-\alpha }\underset{{\Vert \mathbf{z}\Vert }_{\mathcal{H}}\le C,\hspace{0.17em}|b|\le B}{\inf }\hspace{0.17em}{\Vert J^{*}-{\tilde{J}}_{\mathbf{z},b}\Vert }_{\infty }+\frac{4ϵ}{1-\alpha }.\hfill \end{array}$$(19)

The result is proved in Appendix B. There, we prove a more refined result, replacing the max-norm in the bound above with an “optimally weighted” variant, which can provide a substantially stronger guarantee (Desai et al. 2012). The remainder of this section is dedicated to parsing and interpreting this guarantee:

1. Optimal approximation error: Ignoring the ϵ-dependent error term in the guarantee, we see that the quality of approximation provided by $(\widehat{\mathbf{z}},\widehat{b})$ is essentially within a constant multiple (at most $(3+\alpha )/(1-\alpha )$) of the optimal (in the sense of ${\mathcal{\ell }}_{\infty }$-error) approximation to $J^{*}$ possible using a weight vector $\mathbf{z}$ and offsets b permitted by the regularization constraints. This is a “structural” error term that will persist, even if one were permitted to draw an arbitrarily large number of samples. It is analogous to the approximation results produced in parametric settings, with the important distinction that one allows comparisons to approximations in potentially full-dimensional sets, which might, as we have argued earlier, be substantially superior.

2. Dimension-independent sampling error: In addition to the structural error above, one incurs an additional additive “sampling” error that scales like $4ϵ/(1-\alpha )$. The result demonstrates that

   $$ϵ=O\left(\frac{(CK+B)\sqrt{\ln 1/\delta }}{\sqrt{N}}\right).$$

   This is an important contrast with existing parametric sample complexity guarantees. In particular, existing guarantees typically depend on the dimension of the space spanned by the basis function architecture $\{\mathbf{x}\hspace{0.17em}:\hspace{0.17em}x\in \mathcal{S}\}$. Here, this space may be full-dimensional, so that such a dependence would translate to a vacuous guarantee. Instead, we see that the dependence on the approximation architecture is through the constants C, K, and B. Of these, K can, for many interesting kernels, be upper-bounded by a constant that is independent of $|\mathcal{S}|$, whereas C and B are user-selected regularization bounds. Put another way, the guarantee allows for arbitrary “simple” (in the sense of ${\Vert \mathbf{z}\Vert }_{\mathcal{H}}$ being small) approximations in a rich feature space, as opposed to restricting us to some a priori fixed, low-dimensional feature space. This yields some intuition for why we expect the approach to perform well, even with a relatively general choice of kernel. Although the dependence on C, K, and B in the bound above is likely loose, the guarantees above show that for a fixed approximation architecture—that is, choice of kernel, C, and B—the sampling error ϵ improves with the number of samples. We see in our experiments that, keeping other parameters fixed, policy performance improves with the number of samples, N.

3. A nonparametric interpretation: As we have argued earlier, in the event that $\text{span}\hspace{0.17em}\{\mathbf{x}\hspace{0.17em}:\hspace{0.17em}x\in \mathcal{S}\}$ is full-dimensional, there exists a choice of C and B for which the optimal approximation error is, in fact, zero. A large number of kernels would guarantee such feature maps. More generally, as C and B grow large, ${\inf }_{{\Vert \mathbf{z}\Vert }_{\mathcal{H}}\le C,|b|\le B}{\Vert J^{*}-{\tilde{J}}_{\mathbf{z},b}\Vert }_{\infty }$ will decrease. In order to maintain the (bound on) sampling error constant, one would then need to increase N at a rate that is roughly $\mathrm{\Omega }({(CK+B)}^2)$. In summary, by increasing the number of samples in the sampled program, we can (by increasing C and B appropriately) hope to compute approximations of increasing quality, approaching an exact approximation. Although the rate at which the structural error term goes to zero with B remains unclear, this is a substantially weaker requirement than choosing a low-dimensional approximation architecture that provides a small approximation error.

4. Numerical Scheme

This section outlines an efficient numerical scheme that we use to solve the regularized SALP. In particular, we would like to solve the sampled dual Problem (9), introduced in Section 2.3, in order to find an approximation to the optimal cost-to-go function. This approach requires solving a quadratic program (QP) with N × A variables, where $N\triangleq |\widehat{\mathcal{S}}|$ is the number of sampled states and $A\triangleq |\mathcal{A}|$ is the number of possible actions. Furthermore, constructing the coefficient matrices Q and R for (9) requires $O(N^2{A}^2{H}^2)$ arithmetic operations, where H is the maximum number of states that can be reached from an arbitrary state-action pair; that is,

These computationally expensive steps may prevent scaling up solution of the QP to a large number of samples. Also, an off-the-shelf QP solver will typically store the matrix Q in memory, so that the memory required to solve our QP with an off-the-shelf solver effectively scales like $O(N^2{A}^2)$.

In this section, we develop an iterative scheme to solve the Program (9) that, by exploiting problem structure, enjoys low computational complexity per iteration and attractive memory requirements. Our scheme is an active-set method in the vein of the approaches used by Osuna et al. (1997) and Joachims (1999), among others, for solving large support vector machine (SVM) problems. The broad steps of the scheme are as follows:

1. At the tth iteration, a subset $\mathcal{B}\subset \widehat{\mathcal{S}}\times \mathcal{A}$ of the decision variables of (9)—the “active set”—is chosen. Only variables in this set may be changed in a given iteration; these changes must preserve feasibility of the new solution that results. Our algorithm will limit the size of the active set to two variables—that is, $|\mathcal{B}|=2$. The methodology for choosing this active set is crucial and will be described in the sequel.

2. Given the subset $\mathcal{B}$, we solve (9) for ${\lambda }^{(t)}$, where all variables except those in $\mathcal{B}$ must remain unchanged. In other words, ${\lambda }_{x,a}^{(t)}\triangleq {\lambda }_{x,a}^{(t-1)}$ for all $(x,a)\notin \mathcal{B}$. This entails the solution of a QP with only $|\mathcal{B}|$ decision variables. In our case, we will be able to solve this problem in closed form.

3. Finally, if the prior step does not result in a decrease in objective value, we conclude that we are at an optimal solution; Proposition 2 establishes that this is, in fact, a correct termination criterion.

In the following section, we will establish an approach for selecting the active set at each iteration and show how Step 2 above can be solved while maintaining feasibility at each iteration. We will establish that Steps 1 and 2 together need no more than $O(\mathit{NA}\log NA)$ arithmetic operations and comparisons and, moreover, that the memory requirement for our procedure scales like O(NA). Finally, we will establish that our termination criterion is correct: In particular, if no descent direction of cardinality at most two exists at a given feasible point, we must be at an optimal solution.

4.1. Subset Selection

The first step in the active-set method is to choose the subset $\mathcal{B}\subset \widehat{\mathcal{S}}\times \mathcal{A}$ of decision variables to optimize over. Given the convex objective in (9), if the prior iteration of the algorithm is at a suboptimal point $\lambda \triangleq {\lambda }^{(t-1)}$, then there exists a direction $d\in {\mathbb{R}}^{\widehat{\mathcal{S}}\times \mathcal{A}}$ such that $\lambda +ϵd$ is feasible with a lower objective value for $ϵ>0$ sufficiently small. To select a subset to optimize over, we look for such a descent direction d of low cardinality ${\Vert d_0\Vert }_0\le q$—that is, a vector d that is zero on all but at most q components. If such a direction can be found, then we can use the set of nonzero indices of d as our set $\mathcal{B}$.

This problem of finding a “sparse” descent direction d can be posed as

Here, $h(\lambda )\triangleq Q\lambda +R$ is the gradient of the objective of (9) at a feasible point λ; thus, the objective $h{(\lambda )}^{\top }d$ seeks to find a direction d of steepest descent. The first three constraints ensure that d is a feasible direction. The constraint ${\Vert d\Vert }_0\le q$ is added to ensure that the direction is of cardinality at most q. Finally, the constraint ${\Vert d\Vert }_{\infty }\le 1$ is added to ensure that the program is bounded and may be viewed as normalizing the scale of the direction d.

The Program (20) is, in general, a challenging mixed-integer program because of the cardinality constraint. Joachims (1999) discusses an algorithm to solve a similar problem of finding a low-cardinality-descent direction in an SVM classification setting. Their problem can be easily solved, provided that the cardinality q is even; however, no such solution seems to exist for our case. However, in our case, when q = 2, there is a tractable way to solve (20). We will restrict attention to this special case—that is, consider only descent directions of cardinality two. In Appendix D, we will establish that this is, in fact, sufficient: If the prior iterate λ is suboptimal, then there will exist a direction of descent of cardinality two.

To begin, define the sets

Consider the following procedure:

1. Sort the set of indices $\widehat{\mathcal{S}}\times \mathcal{A}$ according to their corresponding component values in the gradient vector $h(\lambda )$. Call this sorted list ${\mathcal{L}}_1$.

2. Denote by (x₁, a₁) the largest element of ${\mathcal{L}}_1$ such that $(x_1,a_1)\notin {\mathcal{P}}_1$, and denote by (x₂, a₂) the smallest element of ${\mathcal{L}}_1$ such that $x_2\notin {\mathcal{P}}_2$. Add the tuple $(x_1,a_1,x_2,a_2)$ to the list ${\mathcal{L}}_2$.

3. Consider all $x\in {\mathcal{P}}_2$. For each such x, denote by $(x,a_1)$ the largest element of ${\mathcal{L}}_1$ such that $(x,a_1)\notin {\mathcal{P}}_1$, and denote by $(x,a_2)$ the smallest element of ${\mathcal{L}}_1$. Add each tuple $(x,a_1,x,a_2)$ to ${\mathcal{L}}_2$.

4. Choose the element $(x_1^{*},a_1^{*},x_2^{*},a_2^{*})$ of ${\mathcal{L}}_2$ that optimizes

   $$\underset{(x_1,a_1,x_2,a_2)\in {\mathcal{L}}_2}{\min }\hspace{0.17em}h{(\lambda )}_{x_2,a_2}-h{(\lambda )}_{x_1,a_1}.$$(21)

   Set $d_{x_1^{*},a_1^{*}}=-1,\hspace{0.17em}{d}_{x_2^{*},a_2^{*}}=1$, and all other components of d to zero.

This procedure finds a direction of maximum descent of cardinality two by examining candidate index pairs $(x_1,a_1,x_2,a_2)$, for which $h{(\lambda )}_{x_2,a_1}-h{(\lambda )}_{x_1,a_1}$ is minimal. Instead of considering all $N^2{A}^2$ such pairs, the routine selectively checks only pairs that describe feasible directions with respect to the constraints of (20). Step 2 considers all feasible pairs with different states x, whereas Step 3 considers all pairs with the same state. It is thus easy to see that the output of this procedure is an optimal solution to (20)—that is, a direction of steepest descent of cardinality two. Further, if the value of the Minimal Objective (21) determined by this procedure is nonnegative, then no descent direction of cardinality two exists, and the algorithm terminates.

In terms of computational complexity, this subset selection step requires us to first compute the gradient of the objective function $h(\lambda )\triangleq Q\lambda +R$. If the gradient is known at the first iteration, then we can update it at each step by generating two columns of Q, because λ only changes at two coordinates. Hence, the gradient calculation can be performed in O(NA) time and with O(NA) storage (because it is not necessary to store Q). For Step 1 of the subset-selection procedure, the component indices must be sorted in the order given by the gradient $h(\lambda )$. This operation requires computational effort of the order $O(\mathit{NA}\log NA)$. With the sorted indices, the effort required in the remaining steps to find the steepest direction via the outlined procedure is O(NA). Thus, our subset-selection step requires a total of $O(\mathit{NA}\log NA)$ arithmetic operations and comparisons.

The initialization of the gradient requires $O(N^2{A}^2)$ effort. In the cases where such a computation is prohibitive, one can think of many approaches to approximately evaluate the initial gradient. For example, if we use the Gaussian kernel, the matrix Q will have most of its large entries close to the diagonal. In this case, instead of the expensive evaluation of Q, we can only evaluate the entries $Q(x,a,x\prime ,a\prime )$, where $x=x\prime$, and set the rest of them to zero. This block diagonal approximation might be used to initialize the gradient. Another approach is to sample from the distribution induced by λ to approximately evaluate $Q\lambda$. As the algorithm makes progress, it evaluates new columns of Q. With a bit of bookkeeping, one could get rid of the errors associated with the gradient initialization. The convergence properties of the active-set method with approximate initialization is an issue not tackled in this paper.

4.2. QP Subproblem

Given a prior iterate ${\lambda }^{(t-1)}$ and a subset $\mathcal{B}\triangleq \{(x_1,a_1),(x_2,a_2)\}$ of decision variable components of cardinality two, as computed in Section 4.1, we have the restricted optimization problem

This subproblem has small dimension. In fact, the equality constraint implies that ${\lambda }_{x_1,a_1}+{\lambda }_{x_2,a_2}$ is constant; hence, the problem is, in fact, a one-dimensional QP that can be solved in closed form. Further, to construct this QP, two columns of Q are required to be generated. This requires computation effort of order $O(NA)$.

Note that the subset $\mathcal{B}$ is chosen so that it is guaranteed to contain a descent direction, according to the procedure in Section 4.1. Then, the solution of (20) will produce an iterate ${\lambda }^{(t)}$ that is feasible for the original Problem (22) and has lower objective value than the prior iterate ${\lambda }^{(t-1)}$.

Finally, we state a proposition establishing the correctness of our active-set method: If the prior iterate $\lambda \triangleq {\lambda }^{(t-1)}$ is suboptimal, then there must exist a direction of descent of cardinality two. Our iterative procedure will therefore improve the solution and will only terminate when global optimality is achieved.

Proposition 2.

If $\lambda \in {\mathbb{R}}^{\widehat{\mathcal{S}}\times \mathcal{A}}$ is feasible, but suboptimal, for (9), then there exists a descent direction of cardinality two.

This result is proved in Appendix D.

5. Case Study: A Queuing Network

This section considers the problem of controlling the so-called “Rybko-Stolyar” queuing network illustrated in Figure 1, with the objective of minimizing long-run average delay. There are two “flows” in this network: the first through server 1 followed by server 2 (with buffering at queues 1 and 2, respectively), and the second through server 2 followed by server 1 (with buffering at queues 4 and 3, respectively). Here, all interarrival and service times are exponential with rate parameters summarized in Figure 1.

This specific network has been studied by de Farias and Van Roy (2003), Chen and Meyn (1998), and Kumar and Seidman (1990), for example, and closely related networks have been studied by Harrison and Wein (1989), Kushner and Martins (1996), Martins et al. (1996), and Kumar and Muthuraman (2004). It is widely considered to be a challenging control problem. As such, a lot of thought has been invested in designing scheduling policies with networks of this type in mind. Our goal in this section will be two-fold. First, we will show that the RSALP, used “out-of-the-box” with a generic kernel, can match or surpass the performance of tailor-made heuristics and state-of-the-art parametric ADP methods. Second, we will show that the RSALP can be solved efficiently.

5.1. MDP Formulation

Although the control problem at hand is nominally a continuous time problem, it is routinely converted into a discrete time problem via a standard uniformization device; see Moallemi et al. (2008), for instance, for an explicit example. In the equivalent discrete time problem, at most a single event can occur in a given epoch, corresponding either to the arrival of a job at queues 1 or 4 or the arrival of a service token for one of the four queues with probability proportional to the corresponding rates. The state of the system is described by the number of jobs in each of the four queues, so that $\mathcal{S}\triangleq {\mathbb{Z}}_{+}^4$, whereas the action space $\mathcal{A}$ consists of four potential actions, each corresponding to a matching between servers and queues. We take the single-period cost as the total number of jobs in the system, so that $g_{x,a}={\Vert x\Vert }_1$. Finally, we take $\alpha \triangleq 0.9$ as our discount factor.

5.2. Approaches

The following scheduling approaches were considered for the queueing problem:

5.2.1. RSALP (This Paper).

We solve (9) using the active-set method outlined in Section 4, taking as our kernel the standard Gaussian kernel $K(x,y)\triangleq \exp (-{\Vert x-y\Vert }_2^2/h)$, with the bandwidth parameter⁶ $h\triangleq 100$. Note that this implicitly corresponds to an full-dimensional-basis function architecture. Because the idealized sampling distribution, ${\pi }_{{\mu }^{*},\nu }$, is unavailable to us, we use in its place the geometric distribution

with the sampling parameter ζ set at 0.9, and take $\nu =\pi$. This choice mimics that of de Farias and Van Roy (2003). The regularization parameter Γ was chosen via a line-search;⁷ we report results for $\mathrm{\Gamma }\triangleq {10}^{-6}$. In accordance with the theory, we set the constraint-violation parameter $\kappa \triangleq 2/(1-\alpha )$, as suggested by the analysis of Section 3, as well as by Desai et al. (2012), for the SALP.

5.2.2. SALP (Desai et al. 2012).

The SALP Formulation (2) is, as discussed earlier, the parametric counterpart to the RSALP. It may be viewed as a generalization of the ALP approach proposed by de Farias and Van Roy (2003) and has been demonstrated to provide substantial performance benefits relative to the ALP approach. Our choice of parameters for the SALP mirrors those for the RSALP to the extent possible, so as to allow for an “apples-to-apples” comparison. Thus, as earlier, we solve the sample average approximation of this program using the same sampling distribution π as in (23), and we set $\kappa \triangleq 2/(1-\alpha )$. Being a parametric approach, one needs to specify an appropriate approximation architecture. Approximation architectures that span (at least quadratic) polynomials are known to work well for queueing problems. We use the basis functions used by de Farias and Van Roy (2003) for a similar problem modeled directly in discrete time. In particular, we use all monomials with degree at most 3, which we will call the cubic basis, as our approximation architectures.

5.2.3. Longest Queue First (Generic).

This is a simple heuristic approach: At any given time, a server chooses to work on the longest queue from among those it can service.

5.2.4. Max-Weight (Tassiulas and Ephremides 1992).

Max-Weight is a well-known scheduling heuristic for queueing networks. The policy is obtained as the greedy policy with respect to a value-function approximation of the form

given a parameter $ϵ>0$. This policy has been extensively studied and shown to have a number of good properties—for example, being throughput optimal (Dai and Lin 2005) and offering good performance for critically loaded settings (Stolyar 2004). Via a line-search, we chose to use $ϵ\triangleq 1.5$ as the exponent for our experiments.

5.3. Results

Policies were evaluated by using a common set of arrival-process sample paths. The performance metric we report for each control policy is the long-run average number of jobs in the system under that policy,

where we set $T\triangleq 10,000$. We further average this random quantity over an ensemble of 300 sample paths. The performance measure we report here is the average cost associated with a particular policy; the RSALP is optimizing over the discounted infinite-horizon cost.

Further, in order to generate SALP and RSALP policies, state sampling is required. To understand the effect of the sample size on the resulting policy performance, the different sample sizes listed in Table 2 were used. Because the policies generated involve randomness to the sampled states, we further average performance over 10 sets of sampled states. The results are reported in Tables 1 and 2 and have the following salient features:

1. RSALP outperforms established policies: Approaches such as the Max-Weight or “parametric” ADP with basis-spanning polynomials have been previously shown to work well for the problem of interest. We see that the RSALP with more than 3,000 samples achieves performance that is superior to these extant schemes.

2. Sampling improves performance: This is expected from the theory in Section 3. Ideally, as the sample size is increased, one should relax the regularization. However, for our experiments, we noticed that the performance is quite insensitive to the parameter Γ. Nonetheless, it is clear that larger sample sets yield a significant performance improvement.

3. RSALP is less sensitive to state sampling: We notice from the standard deviation values in Table 2 that our approach gives policies whose performance varies significantly less across different sample sets of the same size.

Table 1. Performance Results in the Queueing Example

Table 2. Performance Results in the Queueing Example

In Figure 2, we see how performance varies for the RSALP method, as a function of the discount factor α and the regularization parameter Γ.

First, consider the choice of regularizer Γ. From Figure 2, we see that performance degrades for large choices of the parameter (presumably because this restricts the richness of the approximation architecture), as well as for very small choices of the parameter (presumably due to poor generalization). Although there is degradation in performance for small values of Γ, this degradation is meaningful, but limited, in our experiments (average cost is $\sim 20\%$ worse for small choice than the optimal choice).

Now, consider the choice of discount factor α. Note that we are concerned with (and report) average cost performance, but solve a discounted problem formulation. Although this may appear to be inconsistent at first glance, note that the average cost criterion can be seen as a weighted average, across states, of the discounted cost incurred starting at those states. The weights in this average correspond to the steady-state distribution induced by the policy in question. This is in keeping with Desai et al. (2012) and de Farias and Van Roy (2003), for example.

Moreover, we view the choice of α as an algorithmic tuning parameter that implicitly regularizes a bias-variance tradeoff in the learning algorithm. Specifically, if α is too low, there is “bias” introduced through the misalignment of discounted and average cost objectives. On the other hand, if α is too large, the sample complexity degrades, and many more samples are needed. This is consistent with the approximation guarantees we have developed and also consistent with the results in Figure 2: Values of α closer to one result in degraded average cost performance. This is in keeping with the view of the discount factor as a regularizer in the broader reinforcement-learning literature (e.g., Jiang et al. 2015, Amit et al. 2020, and Dewanto and Gallagher 2021).

In summary, we view these results as indicative of the possibility that the RSALP may serve as a practical and viable alternative to state-of-the-art parametric ADP techniques.

6. Conclusions

This paper set out to present a practical nonparametric algorithm for approximate dynamic programming, building upon the success of linear programming-based methods that require the user to specify an approximation architecture. We believe that the RSALP, presented and studied here, is such an algorithm. In particular, the theoretical guarantees presented here establish the “nonparametric” nature of the algorithm by showing that increased sampling effort leads to improved approximations. On the empirical front, we have shown that our essentially “generic” procedure was able to match the performance of tailor-made heuristics as well as ADP approaches using prespecified approximation architectures. Nevertheless, several interesting directions for future work are evident at this juncture:

• The choice of kernel: The choice of kernel matters insomuch as the feature map it encodes allows for approximations with small Hilbert norm (i.e., small C). Thus, a good choice of kernel would require fewer samples to achieve a fixed approximation accuracy than a poor choice of kernel. That said, picking a useful kernel is an apparently easier task than picking a low-dimensional architecture—there are many full-dimensional kernels possible, and with sufficient sampling, they will achieve arbitrarily good value-function approximations. Nevertheless, it would be valuable to understand the interplay between the choice of kernel and the corresponding value of C for specific problem classes (asking for anything more general appears rather difficult).

• The active-set method: In future work, we would like to fine-tune/build more extensible software implementing our active-set method for wider use. Given the generic nature of the approach here, we anticipate that this can be a technology for high-dimensional MDPs that can be used out of the box.