Fast Approximate Dynamic Programming for Infinite-Horizon Markov Decision Processes
Authors: Mohamad Amin Sharifi Kolarijani
Gyula Felix Max
Peyman Mohajerin Esfahani
Venue: NeurIPS 2021 Poster (2021)
Keywords: Stochastic optimal control
value iteration
input-affine systems
Fenchel duality
computational complexity.

Abstract: In this study, we consider the infinite-horizon, discounted cost, optimal control of stochastic nonlinear systems with separable cost and constraints in the state and input variables. Using the linear-time Legendre transform, we propose a novel numerical scheme for implementation of the corresponding value iteration (VI) algorithm in the conjugate domain. Detailed analyses of the convergence, time complexity, and error of the proposed algorithm are provided. In particular, with a discretization of size $X$ and $U$ for the state and input spaces, respectively, the proposed approach reduces the time complexity of each iteration in the VI algorithm from $O(XU)$ to $O(X+U)$, by replacing the minimization operation in the primal domain with a simple addition in the conjugate domain.

Fast Approximate Dynamic Programming for Infinite-Horizon Markov Decision Processes
M. S. Kolarijani
G. Max
Delft Center for Systems and Control Delft University of Technology
Delft Center for Systems and Control Delft University of Technology
The Netherlands
The Netherlands
M. A. A.A.Sharifikole. @awaannnttueffff...
G.F. F....Mxx Maax@tudelee
P. Mohajerin Esfahani
Delft Center for Systems and Control Delft University of Technology
The Netherlands
P. je fahan wtw.ee.fff....
Abstract
this study, we consider the infinite-horizon, discounted cost, optimal control of stochastic nonlinear systems with separable cost and constraints in the state and input variables. Using the linear-time Legendre transform, we propose a novel numerical scheme for implementation of the corresponding value iteration (VT) algorithm in the conjugate domain. Detailed analyses of the convergence, time complexity, and error of the proposed algorithm are provided. In particular, with a discretization of size X and U for the state and input spaces, respectively, the proposed approach reduces the time complexity of each iteration in the VI algorithm from (XU) to O(X + U), by replacing the minimization operation in the primal domain with simple addition in the conjugate domain.
Introduction
Value iteration (VT) is one of the most basic and wide-spread algorithms employed for tackling problems reinforcement learning (RL) and optimal control 99 28 formulated as Markov decision processes (MDPs). The VI algorithm simply involves the consecutive applications of the dynamic programming (DP) operator
TJ(x)) = min (C(x,, ue) + (E(4)))) ar
where C( (z) the cost of taking the control action u, the state 24. This fixed point iteration known to converge to the optimal value function for discount factors E (0, 1). However, this algorithm suffers from high computational cost for large-scale finite state spaces. For problems with a continuous state space, the DP operation becomes an infinite-dimensionall optimization problem, rendering the exact implementation of VI impossible in most cases. A common solution is to incorporate function approximation techniques and compute the output of the DP operator for finite sample (i.e., a discretization) of the underlying continuous state space. This approximation again suffers from a high computational cost for fine discretizations of the state space, particularly in high- aigh-dimensional problems. We refer the reader
9 26 for various approximations of VI.
For some problems, however, is possible to partially address this issue by using duality theory, ie., approaching the minimization problem in the conjugate domain. In partieular, will see
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

in Section 3 the minimization in the primal domain in DP can be transformed to simple addition in the dual domain, at the expense of three conjugate transforms. However. proper application of this transformation relies on efficient numerical algorithms for conjugation. Fortunately, such an
algorithm, known as linear-time Legendre transform (LLTT), has been developed in late 90s [23]. Other than the classical application of LLT (and other fast algorithms for conjugate transform) in solving Hamilton-Jace Jacobi equation [3 [4], these algorithms are used in image processing [24).
thermodynamics [12) and optimal transport [17].
The application of conjugate duality for the DP problem not new and actually goes back to Bellman 4 Further applications of this idea for reducing the computational complexity were later explored in [T3] [8]. However, surprisingly, the application of LLT for solving discrete-time optimal control problems, has been limited. In particular, in [. the authors propose the "fast value iteration" algorithm (without rigorous analysis of the complexity and error of the proposed algorithm) for a particular class of infinite-horizz stage cost C(a,, u) = C(u) and deterministic linear dynamics It+1 = Ar, + Bug, where A is a non-negative, monotone, invertible matrix. More recently, in [20]. we also considered the application of LLT for solving the DP operation in finite-horizon, optimal control of input-affine dynamics Tz-1 1 fax))) + Bu, with separable cost C(x, u) (x + Cu)). In particular, we introduced the "discrete conjugate DP" (d-CDP) operator, and provided a detailed analysis of its complexity and eиror. As we will discuss shortly, the current study is extension of the corresponding d-CDP algorithm that, among other things, considers infinite horizon, discounted cost problems. We note that the algorithms developed in [16] 24 for distance transform can also potentially tackle the optimal control problems similar to the ones of interest in the current study. In particular, these algorithms require the stage cost to be reformulated and next states. While this property might arise naturally, it can generally be restrictive, as it is in the problem class considered in this study, Another line of work that is closely related to ours invloves utilizing max-plus algebra in solving deterministic, continuous-state continuous-time, optimal control problems; see, e.g..
optimal control problems with state-independent
convex function of the "distance" between the current
2 25 These works exploit the compatibility of the DP operation with max- plus operations, and approximate the value function as max-plus linear combination. Recently, in S.. the authors used this idea to propose an approximate VI algorithm for continuous-ste deterministic MDPs. In this regard, we note that the proposed approach in the current study also involves approximating the value function as a max-plus linear combination, namely, the maximum of affine functions. The key difference is however that by choosing a grid-like (factorized) set of slopes for the linear terms (i.e., the basis of the max-plus linear combination). we take advantage of linear time complexity of LLT in computing the constant terms (i.e., the coefficients of the max-p max-plus linear combination).
Main contribution. In this study, we focus on an approximate implementation of VI involving discretization of the state and input spaces for solving the optimal control problem of discrete- systems, with continuous state-input space. Building upon our earlier work (20), we employ conjugate duality to speed-up VI for problems with separable stage cost (in state and input) and input-affine
dynamics. We propose the conjugate VI (ConjVI) algorithm based on a modified version of the d-CDP operator introduced in 20), and extend the existing results in three directions: We consider infinite-horizon, discounted cost problems with stochastic dynamics, while incorporating a numerical scheme for approximation of the conjugate of imput cost. The main contributions of this paper are then as follows: (i) we provide sufficient conditions for the convergence of ConjVI (Theorem 3.9); (ii) we show that ConjVI can achieve a linear time complexity of O(X +U) in each iteration (Theorem 3. 10) compared to the quadratic time complexity of OXU)) of the standard VI, where X and U are the cardinalities of the discrete state and input spaces, respectively; (iii) we analyze the error of ConjVI (Theorem B.11). and use that result to provide specific guidelines the construction of the discrete dual domain (Section 3.4); (iv) we provide a MATLAB package for the implementation of the proposed ConjVI algorithm.
Notations. The standard inner product in R2 and the corresponding induced 2-norm are denoted by (---) and 11-1|g, respectively. I|-llsc denotes the infinity norm. We use the superscript d to denote
finite (discrete) sets (as in x and discrete functions (as in hd x + IR). We use the superscript g to denote grid-like finite sets (as in X mT,, X where X c RR.. We also use X sub, to denote the sub-grid of g derived by omitting the smallest and the largest elements of X in each dimension. The
cardinality of the finite set x or X is denoted by X. We use hd
R" + k = RU {co] to denote a
generic extension of a discrete function hd, and hd R + R to denote multilinear interpolation
2

---

and extrapolation (LERP) extension of a discrete function with grid-like domain. Let X, Y be two arbitrary sets in R co(X) is the convex hull of X. We use d(X, Y) == inf infzexx.ev || |-- 8|12 to denote the distance between X and Y. The one-sided Hausdorff distance from X to V is defined as (X, Y) == > yer ||2 - y||2- Let h R - R be an extended real-valued function with a non-empty effective domain dom(h)- = X == {ar E R" h(x)) < oc). The range of h is denoted by rug(h) = max h(z)- - minzex ((x), and the subdifferential of h at a point æ E X is defined as 0h(x)
{v E R" h(#) 2 h(z) + (y, E - z). Vē E x}. We define L(h) == 1 [L, (h), L (h)]. where L (h) (resp. L, (h)) is the maximum (resp. minimum) slope of the function along the i-th dimension. Note that ah(x) C L(h) for all x E X. We report the complexities using the standard big-O notations O and O where the latter hides the logarithmic factors. In this study, we are mainly concerned with the dependence of the computational complexities on the size of the finite sets involved (discretization of the primal and dual domains). In particular, we ignore the possible dependence of the computational complexities on the dimension of the variables, unless they appear in the power of the size of those discrete sets.
We note that the extended version of this article, including the technical proofs, is available in the supplementary material.
2 VI in primal domain
We are concerned with the infinite-horizon, discounted cost, optimal control problems of the form
a
J+(x)= min E Σ CIx, 24) 10= = x
S.L. Fr+-1 = g(ze, U2 we); Te E X, u. U, we ~ P(W), W 10.
where 2, E R", Ur RT, and w, E R are the state, input and disturbance variables at time I, respectively; 7 E (0,1) is the discount factor; C² X x U - R is the stage cost; Rn x Ih" x R + Rd describes the dynamics; X C RR" and U c Rm describe the state and input constraints, respectively; and, IP(.) is the distribution of the disturbance over the support W R Assuming the stage cost C' is bounded, the optimal value function solves the equation J. = -T. where T is the DP operator (C and J are extended to infinity outside their effective domains) 77 Prop. 1.2.2]
TJ(x) = min (C(r, u) + J (g(ar, u, m,)) Va E X.
(1)
u
Indeed, T is contraaciee the infinity- norm 7 Prop. 1.2.4]. This property then gives rise the VI algorithm J++1 = Ja, which converges to J. as k + co, for arbitrary initialization Jo. Moreover, assuming that the composition J g (for each w) and the cost C' are jointly convex in the state and input variables, T also preserves convexity [8] Prop. 3.
For numerical implementation of VI, we need address three issues. First, we need to compute the expectation in ( In order to simplify the exposition and include the computational cost of this operation explicitly, we consider disturbances with finite support in this study:
Assumption 2.1 (Disturbance with finite support). The disturbance w has a finite support Wd c R with given probability mass function (p.m.f.) p : Wp - [0, 1).
Under the preceding assumption, we have E.. [g(x, u, w)) = = Eweyaa p(w) J(((,, u, w)). The for
second and more important issue is that the optimization problem ①1 is infinite-dimension the continuous state space X. This renders the exact implementation of VI impossible, except for few cases with available closed-for solutions. A common solution to this problem is to deploy a sample-basee based approach, accompanied by a function approximation scheme. To be precise, for a finite subset x of X, at each iteration & 0, L..., take the discrete function J : x IR
as the input, and compute the discrete function J [T]]]d x R. where J : X - RR an extension of Jg. Finally, for each x E x we have solve the minimization problem in (T] over the control input. Since this minimization problem is often a difficult, non-convex problem, a common approximation involves enumeration over discretization Ud c U. Incorporating these approximations, we end up with the approximate VI algorithm Jé = Jd discrete DP (d-DP) operator
JdL characterized by the
(²) = min (((x, u) + Ewewe p(w). (o(e, Vx E xd.
(2)
3

---

The convergence of approximate VI described above depends on the properties of the extension operation H. In partieular, if the extension operation non-expansive (in the infinity-norm), then T is also 7-contractive. The error of this approximation also depends on the extension operation and its representative power. We refer the interested reader to 7 10, 26 for detailed discussions on the convergence and error of different approximation schemes for VI.
The d-DP operator and the corresponding VI algorithm will be our benchmark for evaluating the performance of the alternative algorithm developed in this study. To this end, we finish this section with some remarks the time complexity of the d-DP operation. Let the time complexity of single evaluation of the extension operator [:] in 2 be of (EE). Then, the time complexity of the d-DP operation 2 of O (XUWE). In this regard, note that the scheme described above essentially involves approximating
MDP with finitesstate/action MDP, and then applying the VI algorithm. This, in turn, implies the lower bound S2 for the time complexity (corresponding to enumeration over u E Ud for each æ E X This lower bound is also compatible with the best existing time complexities in the literature for VI for finite MDPs; see, e.g., 3 27]. However, as we will see in the next section, for particular class of problems, it is possible exploit the structure of the underlying continuous system in order to achieve a better time complexity in the corresponding discretized problem.
mattiuouss-sttteacctio
3 Reducing complexity via conjugate duality
In this section, we present the class of problems that allows us to employ conjugate duality and propose an alternative path for solving the corresponding DP operator. We also present the numerical scheme for implementing the proposed alternative path, and analyze its convergence, complexity, and error. We note that the proposed algorithm and its analysis are based on the d-CDP algorithm presented in [20 Sec. 5] for horie-horizon, optimal control of deterministic systems. Here, we extend those results for horizo,, discounted cost, optimal control of stochastic systems. Moreover, unlike [20], our analysis includes the case where the conjugate of input cost is not analytically available and has to be computed numerically; see 20 Assump. 5.1] for more details.
3.1 VI in conjugate domain
Throughout this section, we assume that the problem data satisfy the following conditions.
Assumption 3.1 (Problem class). The problem data has the following properties: (i) The dynamics is of the form g(z, u. w) f(z, u) + w fa(x) + Bu + w, with additive disturbance, where fa R= + R" is Lipschitz continuous, possibly nonlinear map, and B E Raxm (ii) The stage cost C is separable in state and input; that is, C(z, u) = Cz(x) + Ci(u), where the state cost C2 X -> P. and the input cost C, U - R are Lipschitz continuous. (iii) The constraint sets X C R and U C I are compact. Moreover, for each x E X, the set of admissible inputs U(x)= = {u U g(x, u, w) E X, Vu E Wi) is nonempty.
Some remarks are in order regarding the preceding assumptions. We first note that the setting of Assumption 3. gooes beyond the classical LQR. In particular, it includes nonlinear dynamics, state and input constraints, and nou-quadratic stage costs. Second, the properties laid out in Assumption 3.1 imply that the set of admissible inputs U(x) is a compact set for each x E This, in turn, implies that the optimal value in (T) is achieved if J X + R is also assumed to be lower semi-continuous.
For the problem class of Assumption (3.1) we can use duality theory to present an alternative path for computing the output of the DP operator. This path forms the basis for the algorithm proposed in this study. Fix X and consider the following reformulation of the optimization problem ①T)
TJ(x) (x) + min (u) + En + w) : f(, u))
u,z
where we used additivity of disturbance and separability of stage cost. The corresponding dual problem then reads as
Fí(x) = Cx) + max min (Ci(w) + Em + w) + f(z u) - 2)}. u.²
(3)
where y E R is the dual variable corresponding to the equality constraint. For the deterministic dynamics of Assumption (1), we can obtain the following representation for the dual problem.

---

Proposition 3.2 (CDP operator). The dual problem
3 equivalently reads as
<(z) = 7
I EEX,
(4a)
d(y) = C (v),, VER,
(4b)
C.(x) + (Ia(t)) xEX,
(4c)
where [1]* denotes the conjugate operation.
Following (20), we call the operator T in (4) the conjugate DP (CDP) operator. We next provide an alternative representation of the CDP operator that captures the essence of this operation.
Proposition 3.3 (CDP reformulation). The CDP operator T equivalentyy reads as
T(x) = Cx(x) + min (Cr + [. + w)] (f(z,
(5)
where [-]** denotes the biconjugate operation.
The preceding result implies that the indirect path through the conjugate domain essentially involves substituting the input cost and (expectation of the) value function by their biconjugates. In particular, points to a sufficient condition for zero duality gap.
Corollary 3.4 (Equivalence of 'T and 7). FJ = TJ Tf U R and J X - IR are convex.
Hence, T has the same properties as T if C, and J are convex. More importantly, if T and 7 preserve convexity, then the conjugate VI (ConjVI VI) algorithm Jn+1 = T Jk. also converges to the optimal value function J.. with arbitrary convex initialization Jo. For convexity to be preserved, however, need two more additional assumptions. First, the state cost C, : X - R needs to be also convex. Then, for TJ be convex, sufficient condition the convexity of J f (jointly in æ and u), given that J/ is convex. The following assumption summarizes the sufficient conditions for equivalence of VI and ConjVI algorithms.
Assumption 3.5 (Convexity). Consider the following properties for the constraints, costs, and dynamics: (i) The sets X C RR" and C In' are convex. (ii) The costs Cg X - R and C U R are convex. (iii) The deterministic dynamics
R Rm Rh is such that given convex function J % + R, the composition J o is jointly convex in the state and inpul variables.
We note that the last condition in the preceding assumption usually does not hold for nonlinear dynamics; however, for fz(z) = Ax with E RRax,, this is indeed the case for problems satisfying Assumptions 3.1 and 3.5 6 Note that, if convexity not preserved, then the alternative path suffers from duality gap in the sense that in each iteration it uses the convex envelop of (the expectation of) the output of the previous iteration.
3.2 ConjVI algorithm
The approximate ConjVI algorithm involves consecutive applications of an approximate implementa- tion of the CDP operator (4) until some termination condition is satisfied. Algorithm T provides the pseudo-code of this procedure. In particular, we consider solving (4) for a finite set Xd C X, and terminate the iterations when the difference between infinity-norm) is less than a given constant e, > O. Since we are working with a finite subset of the state space, we can restrict the feasibility condition of Assumption 3.1 3.1
consecutive discrete value functions (in the
all Xd.
Assumption 3.6 (Feasibile discretization). We have U(z) + 0 for all E X..
In what follows, we describe the main steps within the initialization and iterations of Algorithm I In particular, the conjugate operations in 4 are handled numerically via the time Legendre transform (LLT)) algorithm (23] LLT an efficient algorithm for computing the discrete conjugate function over finite grid-like dual domain. Precisely, to compute the conjugate of the function h + R, LLT takes its discretization hd xd + R as an input, and outputs Aded Y - R, for the grid-like dual domain Ys, We refer the reader 23 for a detailed description of LLT.
The main steps of the proposed approximate implementation of the CDP operator (4) are as fol- lows: (i) For the expectation operation in (4a), by Assumption 2.1 we again have E (- + w)
5

---

Algorithm 1 ConjVI: Approximate VI in conjugate domain
Input: dynamies fa R= R", E I=*m; finite state space x X; finite input space U C U; state cost function C x - R; imput cost function C -td -> R; finite disturbance space Wd and its p.m.f. p: 87a + [0, 1]; discount factor y termination bound ex.
Output: discrete value function J² x +
initialization:
1: construct the grids Vr 2", and Y6;
2: use LLT to compute Ci y8 R from CP Dd - R:
3: (z) + and (x) t (x) min Cf for X
iteration:
4: while | - J 2 et do
5: + J4:
6:
J + :xa + R according to G [C in (6c) is already computed in line 2):
7:
end while
8:
output Jd < J4.
Emern p(w). J. + w). Hence, we need to pass the value function Jd_ x4 - I through the
"scaled expection filter" to obtain ea x + R in (6a) as an approximation of in 44a Notice that
here we are using an extension Jd= + E of J (recall that we only have access to the discrete value function ) (ii) In order to compute oi in 4b)), we need access to two conjugate functions. First, for E", we use the approximation eded Y
in (6b), by applying LLT to the data points ed X + R for a properly chosen state dual grid Ys C ER".. We also need the conjugate Cj of the input cost. If this function is not analytically available, we approximate it as follows: For a properly chosen input dual grid V8 C I we employ LLT to compute Ciadd 8/8 - R in (6c). using the data points G Ud TR, where [J is finite subset of U. With these conjugate functions at hand, (4b). In particular, notice that we use the LERP extension Cilod of Giood to approximate Cp at the required point - y for each y E Ws, (iii) To be able to compute the output according to (4c). we need to perform another conjugate transform. In particular, we need the value of o* at fz( for æ E X Here, we use the approximation aá Z + IR in (be). by applying LLT to the data points edd Y - R for a properly chosen grid z C R". Finally, we use the LERP extension yodted of yolend to approximate dd at the required point fz( for each x E X and compute Ad in (61) as an approximation of TJ in (c).. With these approximations, we introduce the discrete CDP (d-CDP) operator as follows
can now compute pd Y + R in (6d), an approximation of Φ in
(x))= = 7.
Σ p(u) J + w): xEX. wEW
(6a)
ed)d()) = max { (1, 3) 5 (z)) VE Ys, =E
(6b)
(ol) = max { (u, u) (u))
v vews,
(6c)
Nelid
(g)) /ilnd BB y) + (y))
ye Ys, ZE x E
(6d)
w gddd (z) = max (y,, 2) - (1))) yere
(6c)
()) = Cz(z) + (dd (f(
(6f)
The proper construction of the grids Ys, vs, and 2 will be discussed in Section 3.4
3.3 Analysis of ConjVI algorithm
We now provide our main theoretical results concerning the convergence, complexity, and error of the proposed algorithm. Let us begin with presenting the assumptions to be called in this subsection.
Assumption 3.7 (Grids). Consider the following properties for the grids in Algorithm (consult the
Notations in Section Ty: (i) The grid V/ is constructed such that co (9 %ub) 2 (L())). (i) The grid z5
is constructed such that (Z3) 2 fo (x)). (iii) The construction of V VE, and 23 requires most

---

O(X + U) operations. The cardinality of the grids Y and 2 (resp. 4%) in each dimension is the same as that of (resp. Ua) in that dimension so that Y, Y = X and V = U.
Assumption 3.8 (Extension operator). Consider the following properties for the operator [-] in (6a (i) The extension operator is non-expansive in the infimity norm; (ii) Given a function % - IR and
its discretization Jd x - R. we have || - Ja |a I for some constant Ce 2 0.
Our first result concerns the contractiveness of the d-CDP operator.
Theorem 3,9 (Convergence). Let Assumptions 3.7 (i) and and 3.8 hold. Then, the d-CDP operator (6)
7-contructive in the infinity-norm.
The preceding theorem implies that the approximate ConjVI Algorithm Tis indeed convergent given that the required conditions are satisfied. In particular, for deterministic dynamics, co(2²) fa (x') is sufficient for Algorithm to be convergent. We next consider the complexity of our algorithm.
Theorem 3.10 (Complexity). Let Assumption 3.7 (iäi) hold. Also assume that each evaluation of the extension operator [-] in (6a) requires O(E) operations. Then, the time complexities of initialization and each iteration in Algoritim 1 are of O(X + U) and (XW E), respectively.
The requirements of Assumption (ii)
(iii) will be discussed in Section 3.4 Recall that each iteration of VI (in primal domain) has complexity of O (XUW E), where E denotes the complexity of the extension operation used in ( This observation points to a basic characteristic of the proposed approach: ConjVI reduces the quadratic complexity of VI to linear one by replacing the minimization operation in the primal domain with simple addition in the conjugate domain. Hence, for problem class of Assumption 3.1] ConjVI expected lead reduction in the computational cost. We note that ConjVe VI, like VI and other approximation schemes that utilize discretization/aasttattioo of the continuous state and input spaces, still suffers from the so-called "curse of dimensionality." This because the sizes X and U of the discretizations increase exponentially with the dimensions and m of the corresponding spaces. However, for ConjVI, this exponential increase of rate max (m, n}, compared to the rate m + for VI.
Theorem 3.11 (Error). Let Assumptions 3.5
3.7 (i.ii). and 3.8 (i) hold. Consider the true optimal
value function J. = TJ X
+ R and its discretization J + R. and let Assumption 3.8
hold for J.. Also, let Jd x - TR. be the output of Algorithm Then, [| J (c.tt. F
where eu + ev ex + By + Ez and
eu = cu du (U, U By oy max (Ø(J. ())),))))))))) Ex Ca du (fa(Xd) 2z),
cv du ( (co(V)), VR), ex Cx du (X, xx),
with constants Guy Cv. Cx, Cy. C22 0 depending on the problem data.
Let first note that Assumption 3.5 implies that the DP and CDP operators preserve convexity, and they both have the true optimal value functiin J. as their fixed point (i.c., the duality gap is zero). Otherwise, the proposed scheme can suffer from large errors due to dualization. Moreover, Assumptions ii.ii on the grids V and 2h are required for bounding the error of approximate discrete conjugations using LERP in (6d) and (6t); see the proof of Lemmas A.5 and A.7 of the supplementary material. The remaining sources of error in the proposed approximate implementation of ConjVI are captured by the three main error terms: (i) e, is due to the approximation of the value function using the extension operator [:]: (ii) e, corresponds to the termination of the algorithm after finite number of iterations; (iii) ea captures the error due to the discretization of the primal and dual state and input domains.
3.4 Construction of the grids
In this subsection, we provide specific guidelines for the construction of the grids YR, V and 2 We note that these discrete sets must be grid-like since they form the dual grid for the three conjugate transforms that are handled using L.IT. The presented guidelines aim to minimize the error terms in 7 while taking into account the properties laid out Assumption 3.7 In particular, the schemes described below satisfy the requirements of Assumption 3.7 (ii).
7

---

Construction of Ve. Assumption 3.7 (i) and the error term Ev in 7 suggest that we find the smallest input dual grid V such that co(Vs sub) 2 ((C)). This latter condition essen- tially means that V must "more than cover the range of slope" of the function Ct recall that L(C - [L, (C,),, (G])], where L, (CT)) (resp. I (CI)) is the minimum (resp. maxi- mum) slope of Of along the j-th dimension. Hence, we need to compute/approximate 1 (Ci)) for j ,,..,,,,,, m. A conservative approximation is 7 (C) min aC1/ou; and I. (C) = max 0C//aus. assuming Ci is differentiable. Alternatively, we can directly use the discrete input cost CP for com- puting L (C)).. In particular, if the domain Ud = U I
U3 of Cj is grid-like and Cj is convex, we can take L (Ci') (resp. L (Cf)) to be the minimum first forward difference (resp. maximum last backward difference) of Ci along the j-th dimension (this scheme requires (U) operations). Having I (CF) at our disposal, we can then construct V sub y= v suby such that, in each dimension j.
V Yuuby is uniform and has the same cardinality as Ug. and co co V {uby) |L, (CP),, L (c))]- Finally,
we construct y by extending sub uniformly in each dimension (by adding a smaller and a larger element to v sub in each dimension).
Construction of 25, According to Assumption 3.7 the grid 2 must be constructed such that cO(2R) 2 fx.)). This can be simply done by finding the vertices of the smallest box that contains the set fs (x*). Those vertices give the diameter of 2 in each dimension. We can then, for example, take 2 to be the uniform grid with the same cardinality as Y in each dimension (so that Z = Y). This way, dи (fa (x=),, 23) S du co(2)), 2z), and hence e₂ in 7 reduces by using finer grids 2 This construction has time complexity of (()).
Construction of Ys. Construction of the state dual grid V8 is more involved. According to Theorem 3.11 we need to choose a grid that minimizes ey in (7). This can be done by choosing YE such that VE r a(J. C.) + 0 for all x E X so that ey = 0. Even if we had access to the optimal value function J., satisfying such a condition could lead Such a large size violates Assumption 3.7 on the size of Ys5 and essentially renders the proposed algorithm impractical for dimensions n 2 2. A more practical condition is c0 (J.. C.) + 0 for all x E X so that max d (Ø(J.- - (1)( Ys) S du (co(Ys), Y5), and hence ey reduces by using a finer grid Ys, The latter condition is satisfied if co(2 2 L(J. C.), i... if Y "covers the range of slops" of (J. CC). Hence, we need to approximate the range of slopes of (J.-Ca) To this end, we first use the fact that J, is the fixed point of DP operator [T] to approximate rng(J,- CC.
dual grid YS C RR" of size O(x)).
by R = rmgi (nggg((())))))) We then construct the gird Y = IT WH such that, for each dimension i, 1-7 we have = E co())) where AzA denotes the diameter of the projection of x on the i-th dimension, and > 0 is scaling factor mainly depending on the dimension of the state space. This construction has a one-time computational cost of O(X- + U).
Dynamic construction of Yz. Alternatively, we can construct YE dynamically at each iteration in order to minimize the corresponding error in each application of the d-CDP operator given by ey = cy = mxx d (B(TJ-C)) Y see Lemma A.6 and Proposition A.8 of the supplementary material. This means that the construction of YS Algorithm i moved from line I to inside the iterations, after line 5 Similar to the static scheme described above, the aim here is to construct Ys such that co 2 L[T C.). Since we do not have access to J (it the output of the current iteration), we can again use the definition of the DP operator (1) to approximate Ing(T J Cs)) by R = rng(C,) + + .rng(Jd) where Jd is the output of the previous iteration. We then construct the gird Y = IT Y such that, for each dimension i, we again have =/// E co(Y)). This construction has a one-time computational cost of (U) for computing rng(C,)) and a per iteration computational cost of O(X) for computing rng(Ja). Notice, however, that under this dynamic construction, the error bound of Theorem 3. does not hold true. More importantly, with a dynamic grid Y there is no guarantee for ConjVI to converge.
4 Numerical simulations
We now showcase the application of the ConjVI algorithm through two numerical examples, and compare its performance with the VI algorithm. The details of these simulations (and another numerical example) are provided in Section 4 of the supplementary material. We also provide the
8

---

10
V7
10
CVI
CVIdd
M-dill
10"
Tizse (pwc)
Averua coet
60
10a
50
10
30
10
109
10
10
10f
10²
10
10
(a) Convergence (N 41).
(b) Running time
(c) Average cost (100 runs)
10
CVI CV14
UNI
CVI CVhd
DVIdd
1--1
10
10
Aranuse cost
Tlor
30
10
10"
20
10"4
10 Grd stow X D
10"
10"
10"
10
freention,
Grid size, X
(d) Convergence (N 25).
(c) Running time
(1) Average cost (100 runs)
Figure I: vs. ConjVI (CVI): optimal control inverted pendulum (top) and batch reactor (bottom). The black dotted line (a) corresponds to exponential convergence with rate y = 0.95. CVI-d corresponds to dynamic construction of the dual grid Y in the ConjVI algorithm.
ConjVI MATLAB package [21] for the implementation of the proposed algorithm. The package also includes the numerical simulations of this section. We note that multiple routines in the developed
package are borrowed from the d-CDP MATLAB package 22] Also, for the discrete conjugation
(LLT), we used the MATLAB package (in particular, the LLTd d routine) provided in [23]
We the setup of (20) App. C.2.2] for the optimal control of noisy inverted pendulum with two state variables and one input channel. In particular. we use nearest neighbor extension (which is non-expansive) for the extension operators in (2) for VI and in (6a) for ConjVL. We also consider the dynamic scheme for the construction of Y in ConjVI (hereafter, referred to as Co-jV Moreover, in each implementation of VI and ConjVI(-d) all of the involved discrete sets uniform grids with the same size N in each dimension, i.e., X, Y, Z N2 and U, V = N. We are particularly interested in the performance of these algorithms, as the size of the discretizations increases. The results of our numerical simulations are shown in the top panels of Figure T As shown in Figures Ta both VI and ConjVI algorithms are convergent with rate < 7 = 0.95. particular, Conj VI terminates in 57 iterations, compared to 101 iterations required for VI to terminate. As expected, this faster convergence, combined with the lower time complexity of ConjVI in each itcration, leads to a significant reduction in the running time of this algorithm compared to VI; see Figure Since we not have access to the true optimal value function, in order to evaluate the performance of the outputs of the these algorithms, we consider the performance of the greedy policy with respect to the discrete value function Jdl computed using these algorithms. Figure Ic reports the average cost of 100 random instances over T = 100 steps. As shown, the reduction in the running time in ConjVI comes with an increase in the cost of the controlled trajectories, particularly for coarse discretizations of the state-action space. On the other hand, using the dymamic scheme for construction of Ys, we see that Conj VI-d, with a small increase in the running time compared to ConjVI, achieves almost the same performance as VI concerning the quality of the greedy actions.
Example 1.
We now consider the optimal control of an unstable deterministic batch reactor with four states and two inputs. The setup is borrowed from 19 Sec. 6]. Once again, all of the involved discrete sets are uniform grids with the same size N in each dimension, i.e., X, Y.Y.Z Z No and U, V N We note that we use multi- -inffr- interpolation and extrapolation for the extension operator in 2 for VI. Due to the extrapolation, the extension operator no longer non-expansive, and hence the convergence of VI is not guaranteed. On the other hand, since the dynamics is deterministic, there no need for extension in ConjVI (the scaled expectation (6a) in ConjVI reduces to the simple scaling dd y ld for deterministic dynamics), and hence the convergence of ConjVI only requires
Example 2.

---

co(2*) 2 fx++). The results of simulations are shown in the bottom panels of Figure Once again, we see the trade-off between the time complexity and the greedy control performance in VI and ConjVI. On the other hand, ConjVI-d has the same control performance as VI with an insignificant increase in running time compared to ConjVI. Moreover, in Figure ld we observe the non-monotone behavior of ConjVI-d. In this regard, recall that when the grid V is constructed dynamically and varies at each iteration, the d-CDP operator is not necessarily contractive, which is here the case for the first six iterations. The VI algorithm is also showing a non-monotone behavior, where for the first nine iterations the d-DP operation is actually expansive. As we noted earlier, this is because the multi-lineal extrapolation operation used for extension is expansive.
5 Final remarks
In this paper, we proposed the ConjVI algorithm which reduces the time complexity of the VI algorithm from (XU) to O(X + U). This better time complexity however comes at the expense of restricting the class of problem. In particular, there are two main conditions that must be satisfied in order to be able apply the ConjVI algorithm: The dynamics must be of the form x+ [((x)) + Bu + w. and the stage cost C( u) = C.(x) + Ci()) must be separable. Moreover, since ConjVI essentially solves the dual problem, for non-convex problems, suffers from a non-zero duality gap. Based on our simulation results, we also notice a trade-off between computational complexity and control action quality: While ConjVI has a lower computational cost, VI generates better control actions. However, the dynamic scheme for the construction of state dual grid Y allows us to achieve almost the same performance as VI when it comes to the quality of control actions, with a small extra computational burden. In what follows, we provide our final remarks on the limitations of the proposed ConjVr VI algorithm and its relation to existing approximate VI algorithms.
Relation to existing approximate VI algorithms. The basic idea for complexity reduction introduced in this study can be potentially combined with and further improve the existing sample-based VI algorithms. These sample- based algorithms solely focus on transforming the optimization in DP problems into computationally tractable ones, and in general, they have a time complexity of O(XU), depending on the product of the cardinalities of the discrete state and action spaces. The proposed ConjVI algorithm, on the other hand, focuses on reducing this time complexity to O(X + U). by avoiding the minimization over input in each iteration. Take, for example, the aggregation technique in (26 Sec. 8.1] that leads to a piece-wise constant approximation of the value function. It is straightforward to combine ConjVI with this type of state space aggregation. Indeed, the first numerical example of Section f (optimal control of inverted pendulum) essentially uses such an aggregation by approximating the value function via nearest neighbor extension.
Cost functions with a large Lipschitz constant. Recall that for the proposed ConjVI algorithm to be computationally efficient, the size Y of the state dual grid Y5 must be controlled by the size X of the discrete state space x (Assumption (ii)). Then, as the range of slope of the value function J. increases, the corresponding error ey in 7 due to discretization of the dual state space increases. The proposed dynamic approach for construction of YE partially addresses this issue by focusing on the range of slope of J in each iteration in order to minimize the discretization error of the same iteration k. However, when the cost function has a large Lipschitz constant, even this latter approach can fail to provide a good approximation of the value function. (See the supplementary material for numerical example).
Gradient-based algorithms for solving the minimization over input. Let us first note that the minimization over u in sample-based VI algorithms usually involves solving a difficult non-convex problem. This is particularly due to that fact that the extension operation employed in these algorithms for approximating the value function using the sample points does not lead in u (e.g., take kernel-based approximations or neural networks). This is why in MDP and RL literature, it is actually quite common to consider a finite action space in the first place [I0 26). Moreover, the minimization over again must be solved for each sample point in each iteration, while application of ConjVI avoids solving this minimization each iteration. In this regard, let us note that ConjVI actually uses a convex approximation of the value function, which allows for application of gradient-based algorithm for minimization over u within the Conj algorithm. However, such an algorithm has per iteration complexity of O(XY) o(x²). which is practically inefficient.
convex function
10

---

Acknowledgments
This research is part of a project that has received funding from the European Research Council (ERC) under the grant TRUST-94979. The authors are also grateful to anonymous reviewers for their comments concerning the three remarks in Section
References
[1] Y. Achdou, F. Camilli, and L. Coriias. On numerical approximation of the Hamilton-Jacon Jacobi-transpat system arising in high frequency approximations. Discrete & Continuous Dynamical Systems-Series B, 19(3), 2014.
[2] M. Akian, S. Gaubert, and A. Lakhoua. The max-plus finite element method for solving deterministic optimal control problems: Basic properties and convergence analysis. SIAM Journal Control and Optimization, 17-2888 2008.
[3] F. Bach. Max-plus matching pursuit deterministic Markov decision processes. arXiv preprint arXiv: Seiv:19066 08524, 2019.
[4] R. Bellman W. Karush. Mathematical programming and maximum transform. Journal of the Society for Industrial and Applied Mathematics, 10(3):550-567, 1962.
[5] E. Berthier and Bach. Max-plus linear approximations for deterministic continuous-state markov decision processes. IEEE Control Systems Letters, pages 1-1, 2020.
[6] D. Bertsekas, I Linear convex stochastic control problems over an infinite horizon. JEEF Transactions on Automatic Control, 18(3):11-311 1973.
[7] D. P. Bertsekas. Dynamic Programming and Optimal Control, Vol. II. Athena Scientific, Belmont, MA, 3rd edition, 2007.
[8] D. Bertsekas. Convex Optimization Theory. Athena Scientific, Belmont, MA, 2009.
[9] D. P. Bertsekas. Reinforcement Learning and Optimal Control. Athena Scientific, Belmont, MA, 2019.
[10] L. Busoniu, R. Babuska, B. De Schutter, and D. Ernst. Reinforcement learning and dynamic programming using function approximators. CRC press, 2017.
[11] R. Carpio and T. Kamihigashi. Fast value iteration: application of Segendre-Fencher duality to class of deterministic dynamic programming problems discrete time, Journal of Difference Equations and Applications, 26(2):209-22) 2020.
[12] L. Contento, A. Ern, and R. Vermiglio. linear- time approximate convex envelope algorithm using double Legendre- Fenchel transform with application to phase separation. Computational Optimization and Applications, 60(1):231-28 2015.
2015.
[13] L. Corrias. Fast Legendre- Fenchell transform and applications to Hamilton-Jacee equations and conserva- tion laws. SIAM Journal Numerical Analysis, 33(4):1534-1555 1996.
[14] G. Costeseque and J.-P. Lebacque. A variational formulation for higher order macroscopic traffic flow models: Numerical investigation. Transportation Research Part B: Methodological,, 70:112 133, 2014.
[15] A. O. Esogbue and C. W. Ahn. Computational experiments with class of dynamic programming algorithms of higher dimensions. Computers & Mathematics with Applications, 19(11):3 -33, 1990.
[16] P. Felzenszwalb and D. Huttenlocher. Distance transforms of sampled functions. Theory of computing. 8(1):415-425 2012.
[17] Jacobs and Léger. A fast approach to optimal transport: The back-and aic--nd--oo forth method. Numerische Mathematik, 146(3):513 544, 2020.
[18] C. M. Klein and T. L.. Morin. Conjugate duality and the curse dimensionality. European Journal of Operational Revearch, 50(2):220- 228, 1991.
[19] A. S. Kolarijani, S. C. Bregman, P. Mohajerin Esfahani, and T. Keviczky. A decentralized event- approach for robust model predictive control. IEEE Transactions on Automatic Control, 65(8):3517-333) 2020.

---

[20] M. A. S. Kolarijani and P. Mohajerin Esfahani. Fast approximate dynamic programming for inpul-affine dynamics. preprint arXiv:2008. 10362, 2021.
[21] M. A. S. Kolarijani and P. Mohajerin Esfahani. Conjugate value iteration (ConjVI) MATLAB package. Licensed under the MIT License, available online at /ttttp:///itbub http:/u....oommmiinnnaaaaaatat jani/ ConjVi 2021.
[22] M. A. Kolarijani and P. Mohajerin Esfahani. Discrete conjugate dynamic programming (d-CDP) MATLAB package. Licensed under the MIT License, available online attttt::///gthub com/
AminKolari aaankollrraann/CCCDD CDP 2021.
[23] Y. Lucet. Faster than Algorithms, 16(2):171-185, 1997.
Legendre transform, the linear-time Legendre transform. Numerical
[24] Y. Lucet. New sequential exact Euclidean distance transform algorithms based on convex analysis. Image and Vision Computing, 27(1):37 44, 2009.
[25] M. McEneaney, Max-plus eigenvector representations solution nonlinear oo problems: basic concepts. IEEE Trunsactions on Automatic Control, 48(7):1150-1163, 2003.
[26] B. Powell. Approximate Dynamic Programming: Solving the Curses of Dimensionality. John Wiley Sons, Hoboken, NJ, 2nd edition. 2011.
[27] Sidford, M. Wang, Wu, and Ye. Variance reduced value iteration and faster algorithms for solving Markov decision processes. In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 770-787. SIAM, 2018.
[28] R. S. Sutton and A. G. Barto. Reinforcement Learning: An Introduction. MIT Press, 2018
12


Reviews:
[{"review_id": "sjF23J5O47", "reviewer": "Reviewer_mb6w", "summary": "This paper proposed a new value iteration algorithm (CVI) for optimal control of stochastic system with continuous state space. The proposed algorithm discretizes the state space and performs the value iteration step in the \"conjugate domain\" to reduce one-step computational complexity. It was shown the proposed algorithm converges under conditions on the stochastic system, discretization, and etc. The proposed algorithm CVI, its variant CVI-d are compared to classical value iteration (VI) over a noisy inverted pendulum experiment.", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well written and the language is precise.\n\nThe main benefit of the propose algorithm is the reduced per-step complexity (Thm 3.10). However, no guarantees on how many iterations will be needed for the algorithm to stop (see Alg.1 line 4) given the error tolerance parameter e_t in Algo 1. Therefore, Thm 3.10 along with Thm 3.11 cannot really justify the claim of \"fast approximate dynamic programming\". More on this point in the experiment part below.\n\nI understand, in the inverted pendulum experiment, the CVI algorithm needs less iteration than VI, but we need more evidence to be convinced this is a fast value iteration algorithm. The experiment also brings another question, Fig.1 c shows the resulting greedy policy learned from CVI is worse than classic VI in all state grid settings; it is much worse with coarse grids. The variant of CVI, CIV-d (with dynamic grid in dual space), was able to find better policies, but there is no guarantee for that.\n\nFor originality, I am not very familiar with literature in this area (sample based value iteration). It will help a lot if the authors clarify the technical contributions of the current algorithm compared to d-CDP from [18, reference from the paper]. Since the authors claim that the current algorithm extends d-CDP from infinite-horizon, discounted cost to stochastic dynamics.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "See the main comments above.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "hpK-RgShPSM", "reviewer": "Reviewer_53cU", "summary": "The paper develops an approximate value iteration method for an infinite-horizon, discounted-cost Markov Decision Processes (MDPs) that satisfy a given set of regulatory assumptions. The idea is to work on a dual space that replaces the DP value function by a reformulation written in terms of biconjugate operations. The authors present an algorithm that applies the conjugate operators iteratively, and present convergence and other structural results associated with the resulting errors. Numerical results compare the proposed approach against a traditional value iteration and a variant that generated dynamic discretization grids.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Originality**\n\nThe paper is an extension the previous work by Kolarijani and Esfahani (2020), which lays out the fundamental ideas underlying conjugate dual operators in MDPs and has similar error-type analysis. The key contribution of this paper is that it extends those ideas in different ways, incorporating stochastic dynamics (more briefly discussed in that paper) and involves a more general class of MDPs (e.g., infinite-horizon and resulting numerical approximations). \n\nWhile such extensions are no doubt a challenging and admirable endeavor, their contribution is much more technical and builds upon the novelty of that previous work. Assuming the previous work is a separate publication and the current one does not encompass that one, my concern is that the originality of this paper is more restricted to an MDP audience, given the (usual) wide array of assumptions that must be in place for this approach to be applicable. \n\n**Quality** and **Clarity**\n\nThe paper is excellently well written and rigorous. It is, however, quite dense to read: some concepts I was only able to grasp after reading the \"extended manuscript\" in the supplemental material, which includes more examples (I appreciate the authors for that). Nonetheless, it is overall a very pleasant and fluid read. There are also six pages of technical proofs, some involved but most using standard contraction and similar operations - I did my best to verify them and they seem correct. \n\n\n**Significance**\n\nAs mentioned before, my major concern is that the paper seems to be more suitable to a more specialized audience of NeurIPS due to its strong technical aspects. The numerical simulations are very limited, and in my view are not very effective in demonstrating the practical relevance of their new approach.\n\nHowever, the paper reads very well and the conjugate reformulation is a relevant (and impactful) contribution. It does reduce the VI complexity significantly and can be of potential use in different dynamic programming models. \n\n**Questions and Minor Points** \n\n- How restrictive is assumption 3.1-(ii)?\n- In Assumption 3.5: \"feasibile\" --> \"feasible\"", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper indeed lacks a discussion of the limitations of their methodology. In particular, authors could be clearer with respect to trade-offs of the conjugate formulation and the traditional VI algorithm.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "F6dKJTnccij", "reviewer": "Reviewer_cP9r", "summary": "The authors consider value iteration [TJ](x) = min_u [c(x,u) + gamma E J(x’)] for problems where the transition function and cost functions are additively separable: c(x,u)=c1(x)+c2(u) and x’=f1(x)+f2(u) + w (Assumptions 3.1). Furthermore, f1(x) is Lipschitz continuous, f2(u) is linear, c1(x) is convex and Lipschitz continuous, and c2(u) is linear in u. The disturbance w has a finite support (Assumption 2.1). The state space for x and decision space for u are continuous. Thus, the authors consider problems such that [TJ](x) = c1(x) + min_u [c2(u) + gamma sum_i  p_i J(f1(x)+Bu+w_i)]. A problem that fits this description is the stabilization of an inverted pendulum with discretized disturbances, studied in Section 4. With this choice of assumptions, and assuming the value function is convex (which holds true if f1(x) is linear or linearized, that is, if x’=Ax+Bu+w), convex duality can be applied to the minimization problem over (u,x”) given x and subject to x”=Ax+Bu, which gives another representation of the optimization problem, based on the convex conjugate of c1 and the convex conjugate of e(x”)=gamma E J(x”+w).\n\nThe authors propose and study an algorithm that relies on computing these convex conjugate functions numerically in order to provide an alternative way to the computation of the value function J. They establish error bounds (Theorem 3.13) in terms of theoretical quantities of the problem data.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The results appear to be new. While dualization is a classical concept in optimization, and there is preliminary work on applying these ideas to deterministic systems, the extension to stochastic systems and the results on error bounds appear to be new.\nThe work is technically sound given the assumptions. Matlab codes are provided. In the opinion of this reviewer, the limitations of the work are\n(i) the scope of problems that can be treated by the proposed approach, compared to classical VI. For instance, the authors state that the curse of dimensionality from classical VI has exponential rate m+n, while their method has exponential rate max(m,n). But classical VI does not rely on a convex value function assumption that restrict the admissible class of transition and cost functions. The test problem presented in the numerical section does not itself satisfy all the assumptions assumed by the theory (as pointed out by line 463 of the extended manuscript).\n(ii) the practicality of the bounds, compared to classical VI. The error bounds of the proposed approach do not seem to easily lead to a stopping criterion as they rely non-universal constants and quantities that seem difficult to evaluate or bound. \n(iii) it is not clear if the assumption that the transition function is linear to ensure convexity of the value function and fully justify a dualization approach is needed for the approach to be sound. The numerical examples provided actually use a system where control techniques using a linearization of the transition function are known to work well. \n(iv) the fundamental reason that would explain that solving in the dual space is easier than solving in the primal space is not fully apparent. The test problem has a unidimensional decision variable, so it is difficult to draw strong conclusions on scaling from there. \nThe submission is clearly written. The results are significant, although they seem to be most relevant for certain classes of control problems that are already well solved (it is hard from the current manuscript to identify problems that would significantly gain from the new approach).\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The contributions are essentially mathematical. The work does not appear to have potential negative societal impact.", "review_text": "", "extra_scores": {"ethical_concerns": "NA", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "6Ue4w9Y_SyH", "reviewer": "Reviewer_9wHP", "summary": "This paper presents a fast method for solving certain classes of optimal control problems, by implementing a value iteration (VI) algorithm that takes advantage of convex duality. \nThe basic idea is to take advantage of the fact that for two functions $f_1, f_2$, we have $(f_1 \\square f_2)^* = f_1^* + f_2^*$, where $f^*$ denotes the Legendre-Fenchel transform of $f$, and $(f_1 \\square f_2)(x) = \\inf \\\\{f_1(x_1) + f_2(x_2) : x_1 + x_2 = x \\\\} $ denotes the infimal convolution of $f_1, f_2$.\nThus, the slow $\\inf$ operation can be replaced with a fast addition operation.\nThe paper provides analyses of the convergence, time complexity, and error of the algorithm.\nWhen compared to a \"naive\" VI algorithm in which each iteration takes $O(X U)$ (where $X$, $U$ are the grid sizes of the state and input spaces, respectively), the algorithm presented in this paper can achieve $O(X+U)$.\nFinally, the paper provides a numerical example implemented in MATLAB.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper presents an fast algorithm that is probably novel to most researchers and practitioners who need to solve optimal control problems.\nThe ideas first appeared in Bellman & Karush (1962), and the algorithm was most recently presented and analyzed in Kolarijani and Esfahani (2020). \nThe latter paper applied the algorithm to finite-horizon optimal control problems with deterministic input-affine dynamics, and provided detailed analyses of the algorithm's time complexity and error properties.\nThis current paper goes beyond that paper in applying the algorithm to discounted, infinite-horizon problems with stochastic dynamics.\nThe marginal contribution of this paper is incremental; however, the novelty and performance of the algorithm may be of interest to researchers who need to solve problems that fall within the class handled by this algorithm.\n\nThe paper is well-written, clear and concise. The mathematical details seem correct (though I did not check them closely). I have addressed some limitations of the paper in the \"Limitations\" section below.\n\nOverall, the paper is interesting and presents a novel algorithm with the potential for large performance gains in solving optimal control problems. The applicability of the algorithm to more complicated, real-world problems could be an interesting and useful research question.\n\nBellman and Karush (1962). Mathematical programming and the maximum transform. Journal\nof the Society for Industrial and Applied Mathematics, 10(3):550-567.  \n\nKolarijani and Esfahani (2020). Fast approximate dynamic programming for input-affine\n409 dynamics. preprint arXiv:2008.10362, 2020.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations\n\n1) The paper clearly states the restriction on the class of problems that the algorithm applies to: principally, the stage cost function $C(x, u)$ must be jointly convex.\nThe baseline comparison algorithm is a naive version of VI. \nA single iteration of VI must solve the Bellman equation\n$$\n  J(x) = \\min_u \\\\{ C(x, u) + \\gamma E\\[ J(g(x, u, w) \\] \\\\}\n$$\nIf $X, U$ are the grid sizes of the state space and input space, respectively, then a brute-force implementation will search through every possible $u$ for every possible grid point value of $x$, giving $O(X U)$ time complexity.\nHowever, since it is assumed that $C(x,u)$ is convex, it is not hard to show that the value function $J$ must also be convex, and thus faster optimization algorithms that use gradient information can be used.\nComparing the performance of CVI against a gradient-aware algorithm would provide a more realistic comparison for the speed advantage of CVI.\n\n\n2) The numerical example provided is a two-dimensional noisy inverted pendulum, which may be too simplistic to provide an indication of the potential strengths and weaknesses of the algorithm.\nIt would be useful to see how the algorithm handles more complicated problems, and how some practical issues in implementing the algorithm are addressed.\nFor example, when constructing the approximate Legendre-Fenchel transform $f^*$, the grid space of $f^*$ must be large enough to encompass the full range of slopes of the primal function $f$.\nFor certain functions, this range may be very large, extending to infinity (e.g. $\\log(x)$ or $x^{-\\alpha}$).\nEspecially in higher dimensions, it is non-trivial to construct an approximation for $f^*$ that encompasses such a large range while maintaining accuracy.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "5hDxCxspxa", "reviewer": "Reviewer_Hfvo", "summary": "This paper studies the problem of dynamic programming with continuous state and action spaces. \nThe authors first propose the VI which naively discretizes the state and action space and performs the classical value iteration. The authors use this method as a baseline for comparison with CVI developed in this paper. \nRegarding the CVI, this algorithm performs VI in the the conjugate domain. In particular, by taking the dual of the objective function, we take the problem to the dual space. By doing so, instead of taking minimum over all the actions, the problem is reduced to a summation. Hence, the complexity of the CVI is better than the VI in the order of the cardinality of the action space.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "- Since the J function can be bounded above by O(1/(1-gamma)), I do not understand the purpose of the Theorem 3.11. Isn't it result in a trivial bound?\n- Comparison in the experimental result is not sufficient. For instance, can the authors compare the performance of their algorithm with the ones in [8] and [22]? Specifically section 8.1 in [22] that talks about aggregation techniques to deal with the curse of dimensionality, and references in section 8.1 in the book. Also look at table 4.4 in [22] for a variety of the references to compare with. \n- There is no Assumption 3.7-(ii), there is no Fig 1-d.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "- Lack of comparison with the related works.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper suggests a fast method for solving a class of optimal control problems by introducing a variant of the Value Iteration (VI) algorithm that benefits from the convex conjugate transform. This improves the computational complexity of solving each step of VI from O(|X||U|) to O(|X| + |U|).

Most reviewers are positive about this work. One of the reviewers (mb6w) is on the negative side, but I believe their concerns and questions are adequately answered (the reviewer did not acknowledge the authors' response). Therefore, I would recommend acceptance of this paper.

I would encourage the authors to incorporate reviewers' comments and suggestions in order to improve their paper. Some of improvements can be along the followings:

- Emphasizing more how this work can be relevant to NeurIPS community, as this paper is different from a typical paper on MDP or RL appearing at NeurIPS (no focus on learning; focusing on a special class of problems, etc.)
- Clarifying and emphasizing the technical contributions of this work compared to [18].
- Including larger-scale experiments.

I would also note that the acronym of CVI has already been used in the context of Value Iteration:
*Conservative Value Iteration* (Kozuno, Uchibe, Donya, "Theoretical Analysis of Efficiency and Robustness of Softmax and Gap-Increasing Operators in Reinforcement Learning," AISTATS, 2019) and *Characteristic Value Iteration* (Farahmand, "Value Function in Frequency Domain and the Characteristic Value Iteration Algorithm," NeurIPS, 2019), so the authors may want to use another acronym to reduce possible confusion. There is no need to cite these work, as they are not related enough.
Decision: Accept (Poster)