Fast Approximate Dynamic Programming for Infinite-Horizon Markov Decision Processes
Year: 2021

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.

heading: 1 Introduction
text: Value iteration (VI) is one of the most basic and wide-spread algorithms employed for tackling problems in reinforcement learning (RL) and optimal control [9, 28] formulated as Markov decision processes (MDPs). The VI algorithm simply involves the consecutive applications of the dynamic programming (DP) operator
T J(xt) = min ut
{ C(xt, ut) + γEJ(xt+1) } ,
where C(xt, ut) is the cost of taking the control action ut at the state xt. This fixed point iteration is known to converge to the optimal value function for discount factors γ ∈ (0, 1). However, this algorithm suffers from a high computational cost for large-scale finite state spaces. For problems with a continuous state space, the DP operation becomes an infinite-dimensional 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 a 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-dimensional problems. We refer the reader to [9, 26] for various approximations of VI.
For some problems, however, it is possible to partially address this issue by using duality theory, i.e., approaching the minimization problem in the conjugate domain. In particular, as we will see
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
in Section 3, the minimization in the primal domain in DP can be transformed to a 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 (LLT), has been developed in late 90s [23]. Other than the classical application of LLT (and other fast algorithms for conjugate transform) in solving Hamilton-Jacobi equation [1, 13, 14], these algorithms are used in image processing [24], thermodynamics [12], and optimal transport [17].
The application of conjugate duality for the DP problem is not new and actually goes back to Bellman [4]. Further applications of this idea for reducing the computational complexity were later explored in [15, 18]. However, surprisingly, the application of LLT for solving discrete-time optimal control problems, has been limited. In particular, in [11], the authors propose the “fast value iteration” algorithm (without a rigorous analysis of the complexity and error of the proposed algorithm) for a particular class of infinite-horizon optimal control problems with state-independent stage cost C(x, u) = C(u) and deterministic linear dynamics xt+1 = Axt + But, 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 xt+1 = fs(xt) + But with separable cost C(x, u) = Cs(x) + Ci(u). In particular, we introduced the “discrete conjugate DP” (d-CDP) operator, and provided a detailed analysis of its complexity and error. As we will discuss shortly, the current study is an 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 as a convex function of the “distance” between the current 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., [2, 25]. These works exploit the compatibility of the DP operation with maxplus operations, and approximate the value function as a max-plus linear combination. Recently, in [3, 5], the authors used this idea to propose an approximate VI algorithm for continuous-state, 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-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-time 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 input 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 O(XU) 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 3.11), and use that result to provide specific guidelines on 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 Rn and the corresponding induced 2-norm are denoted by ⟨·, ·⟩ and ∥·∥2, respectively. ∥·∥∞ denotes the infinity norm. We use the superscript d to denote finite (discrete) sets (as in Xd) and discrete functions (as in hd : Xd → R). We use the superscript g to denote grid-like finite sets (as in Xg = Πni=1X g i where X g i ⊂ R). We also use X g sub to denote the sub-grid of Xg derived by omitting the smallest and the largest elements of Xg in each dimension. The cardinality of the finite set Xd or Xg is denoted by X . We use h̃d : Rn → R = R ∪ {∞} to denote a generic extension of a discrete function hd, and hd : Rn → R to denote multilinear interpolation
and extrapolation (LERP) extension of a discrete function with a grid-like domain. Let X,Y be two arbitrary sets in Rn. co(X) is the convex hull of X. We use d(X,Y) := infx∈X,y∈Y ∥x− y∥2 to denote the distance between X and Y. The one-sided Hausdorff distance from X to Y is defined as dH(X,Y) := supx∈X infy∈Y ∥x− y∥2. Let h : Rn → R be an extended real-valued function with a non-empty effective domain dom(h) = X := {x ∈ Rn : h(x) < ∞}. The range of h is denoted by rng(h) := maxx∈X h(x)−minx∈X h(x), and the subdifferential of h at a point x ∈ X is defined as ∂h(x) := { y ∈ Rn : h(x̃) ≥ h(x) + ⟨y, x̃− x⟩ ,∀x̃ ∈ X } . We define L(h) :=
Πni=1 [ L−i (h),L + i (h) ] , where L+i (h) (resp. L − i (h)) is the maximum (resp. minimum) slope of the function h along the i-th dimension. Note that ∂h(x) ⊆ L(h) for all x ∈ 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.

heading: 2 VI in primal domain
text: We are concerned with the infinite-horizon, discounted cost, optimal control problems of the form
J⋆(x) =min Ewt [ ∞∑ t=0 γtC(xt, ut) ∣∣∣∣x0 = x ]
s.t. xt+1 = g(xt, ut, wt), xt ∈ X, ut ∈ U, wt ∼ P(W), ∀t ∈ {0, 1, . . .}, where xt ∈ Rn, ut ∈ Rm, and wt ∈ Rl are the state, input and disturbance variables at time t, respectively; γ ∈ (0, 1) is the discount factor; C : X × U → R is the stage cost; g : Rn × Rm × Rl → Rn describes the dynamics; X ⊂ Rn and U ⊂ Rm describe the state and input constraints, respectively; and, P(·) is the distribution of the disturbance over the support W ⊂ Rl. Assuming the stage cost C is bounded, the optimal value function solves the equation J⋆ = T J⋆, where T is the DP operator (C and J are extended to infinity outside their effective domains) [7, Prop. 1.2.2]
T J(x) := min u
{ C(x, u) + γ · EwJ ( g(x, u, w) )} , ∀x ∈ X. (1)
Indeed, T is γ-contractive in the infinity-norm [7, Prop. 1.2.4]. This property then gives rise to the VI algorithm Jk+1 = T Jk, which converges to J⋆ as k → ∞, for arbitrary initialization J0. 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.3.1]. For numerical implementation of VI, we need to address three issues. First, we need to compute the expectation in (1). 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 ⊂ Rl with a given probability mass function (p.m.f.) p : Wd → [0, 1].
Under the preceding assumption, we have EwJ ( g(x, u, w) ) = ∑ w∈Wd p(w) · J ( g(x, u, w) ) . The second and more important issue is that the optimization problem (1) is infinite-dimensional for the continuous state space X. This renders the exact implementation of VI impossible, except for a few cases with available closed-form solutions. A common solution to this problem is to deploy a sample-based approach, accompanied by a function approximation scheme. To be precise, for a finite subset Xd of X, at each iteration k = 0, 1, . . ., we take the discrete function Jdk : Xd → R as the input, and compute the discrete function Jdk+1 = [ T J̃dk ]d : Xd → R, where J̃dk : X → R is an extension of Jdk . Finally, for each x ∈ Xd, we have to solve the minimization problem in (1) over the control input. Since this minimization problem is often a difficult, non-convex problem, a common approximation involves enumeration over a discretization Ud ⊂ U. Incorporating these approximations, we end up with the approximate VI algorithm Jdk+1 = T dJdk , characterized by the discrete DP (d-DP) operator
T dJd(x) := min u∈Ud
{ C(x, u) + γ · ∑ w∈Wd p(w) · J̃d ( g(x, u, w) )} , ∀x ∈ Xd. (2)
The convergence of approximate VI described above depends on the properties of the extension operation [̃·]. In particular, if the extension operation is non-expansive (in the infinity-norm), then T d is also γ-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 on the time complexity of the d-DP operation. Let the time complexity of a single evaluation of the extension operator [̃·] in (2) be of O(E). Then, the time complexity of the d-DP operation (2) is of O ( XUWE ) . In this regard, note that the scheme described above essentially involves approximating a continuous-state/action MDP with a finite-state/action MDP, and then applying the VI algorithm. This, in turn, implies the lower bound Ω(XU) for the time complexity (corresponding to enumeration over u ∈ Ud for each x ∈ Xd). 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 a particular class of problems, it is possible to exploit the structure of the underlying continuous system in order to achieve a better time complexity in the corresponding discretized problem.

heading: 3 Reducing complexity via conjugate duality
text: 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 finite-horizon, optimal control of deterministic systems. Here, we extend those results for infinite-horizon, 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.

heading: 3.1 VI in conjugate domain
text: 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(x, u, w) = f(x, u) + w = fs(x) + Bu + w, with additive disturbance, where fs : Rn → Rn is a Lipschitz continuous, possibly nonlinear map, and B ∈ Rn×m. (ii) The stage cost C is separable in state and input; that is, C(x, u) = Cs(x) + Ci(u), where the state cost Cs : X → R and the input cost Ci : U → R are Lipschitz continuous. (iii) The constraint sets X ⊂ Rn and U ⊂ Rm are compact. Moreover, for each x ∈ X, the set of admissible inputs U(x) := {u ∈ U : g(x, u, w) ∈ X, ∀w ∈ Wd} is nonempty.
Some remarks are in order regarding the preceding assumptions. We first note that the setting of Assumption 3.1 goes beyond the classical LQR. In particular, it includes nonlinear dynamics, state and input constraints, and non-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 ∈ X. This, in turn, implies that the optimal value in (1) 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 ∈ X and consider the following reformulation of the optimization problem (1)
T J(x) =Cs(x) + min u,z {Ci(u) + γ · EwJ(z + w) : z = f(x, u)} ,
where we used additivity of disturbance and separability of stage cost. The corresponding dual problem then reads as
T̂ J(x) := Cs(x) + max y min u,z {Ci(u) + γ · EwJ(z + w) + ⟨y, f(x, u)− z⟩} , (3)
where y ∈ Rn is the dual variable corresponding to the equality constraint. For the deterministic dynamics of Assumption 3.1-(i), we can obtain the following representation for the dual problem.
Proposition 3.2 (CDP operator). The dual problem (3) equivalently reads as
ϵ(x) := γ · EwJ(x+ w), x ∈ X, (4a) ϕ(y) := C∗i (−B⊤y) + ϵ∗(y), y ∈ Rn, (4b) T̂ J(x) = Cs(x) + ϕ∗ ( fs(x) ) , x ∈ X, (4c)
where [·]∗ 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̂ equivalently reads as
T̂ J(x) = Cs(x) + min u
{ C∗∗i (u) + γ · [EwJ(·+ w)]∗∗ ( f(x, u) )} , (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, it points to a sufficient condition for zero duality gap.
Corollary 3.4 (Equivalence of T and T̂ ). T̂ J = T J if Ci : U → R and J : X → R are convex.
Hence, T̂ has the same properties as T if Ci and J are convex. More importantly, if T and T̂ preserve convexity, then the conjugate VI (ConjVI) algorithm Jk+1 = T̂ Jk, also converges to the optimal value function J⋆, with arbitrary convex initialization J0. For convexity to be preserved, however, we need two more additional assumptions. First, the state cost Cs : X → R needs to be also convex. Then, for T̂ J to be convex, a sufficient condition is the convexity of J ◦ f (jointly in x 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 ⊂ Rn and U ⊂ Rm are convex. (ii) The costs Cs : X → R and Ci : U → R are convex. (iii) The deterministic dynamics f : Rn × Rm → Rn is such that given a convex function J : X → R, the composition J ◦ f is jointly convex in the state and input variables.
We note that the last condition in the preceding assumption usually does not hold for nonlinear dynamics; however, for fs(x) = Ax with A ∈ Rn×n, this is indeed the case for problems satisfying Assumptions 3.1 and 3.5 [6]. Note that, if convexity is 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.

heading: 3.2 ConjVI algorithm
text: The approximate ConjVI algorithm involves consecutive applications of an approximate implementation of the CDP operator (4) until some termination condition is satisfied. Algorithm 1 provides the pseudo-code of this procedure. In particular, we consider solving (4) for a finite set Xd ⊂ X, and terminate the iterations when the difference between two consecutive discrete value functions (in the infinity-norm) is less than a given constant et > 0. Since we are working with a finite subset of the state space, we can restrict the feasibility condition of Assumption 3.1-(iii) to all x ∈ Xd: Assumption 3.6 (Feasibile discretization). We have U(x) ̸= ∅ for all x ∈ Xd.
In what follows, we describe the main steps within the initialization and iterations of Algorithm 1. In particular, the conjugate operations in (4) are handled numerically via the linear-time Legendre transform (LLT) algorithm [23]. LLT is an efficient algorithm for computing the discrete conjugate function over a finite grid-like dual domain. Precisely, to compute the conjugate of the function h : X → R, LLT takes its discretization hd : Xd → R as an input, and outputs hd∗d : Yg → R, for the grid-like dual domain Yg. We refer the reader to [23] for a detailed description of LLT. The main steps of the proposed approximate implementation of the CDP operator (4) are as follows: (i) For the expectation operation in (4a), by Assumption 2.1, we again have EwJ(· + w) =
Algorithm 1 ConjVI: Approximate VI in conjugate domain
Input: dynamics fs : Rn → Rn, B ∈ Rn×m; finite state space Xd ⊂ X; finite input space Ud ⊂ U; state cost function Cds : Xd → R; input cost function Cdi : Ud → R; finite disturbance space Wd and its p.m.f. p : Wd → [0, 1]; discount factor γ; termination bound et. Output: discrete value function Ĵd : Xd → R. initialization:
1: construct the grids Vg, Zg, and Yg; 2: use LLT to compute Cd∗di : Vg → R from Cdi : Ud → R; 3: Jd(x)← 0 and Jd+(x)← Cds (x)−minCdi for x ∈ Xd;
iteration: 4: while ∥∥Jd+ − Jd∥∥∞ ≥ et do 5: Jd ← Jd+; 6: Jd+ ← T̂ dJd : Xd → R according to (6) [Cd∗di in (6c) is already computed in line 2]; 7: end while 8: output Ĵd ← Jd+.
∑ w∈Wd p(w) · J(· + w). Hence, we need to pass the value function Jd : Xd → R through the “scaled expection filter” to obtain εd : Xd → R in (6a) as an approximation of ϵ in (4a). Notice that here we are using an extension J̃d : X → R of Jd (recall that we only have access to the discrete value function Jd). (ii) In order to compute ϕ in (4b), we need access to two conjugate functions. First, for ϵ∗, we use the approximation εd∗d : Yg → R in (6b), by applying LLT to the data points εd : Xd → R for a properly chosen state dual grid Yg ⊂ Rn. We also need the conjugate C∗i of the input cost. If this function is not analytically available, we approximate it as follows: For a properly chosen input dual grid Vg ⊂ Rm, we employ LLT to compute Cd∗di : Vg → R in (6c), using the data points Cdi : Ud → R, where Ud is a finite subset of U. With these conjugate functions at hand, we can now compute φd : Yg → R in (6d), as an approximation of ϕ in (4b). In particular, notice that we use the LERP extension Cd∗di of C d∗d i to approximate C d∗ i at the required point −B⊤y for each y ∈ Yg. (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 ϕ∗ at fs(x) for x ∈ Xd. Here, we use the approximation φd∗d : Zg → R in (6e), by applying LLT to the data points φd : Yg → R for a properly chosen grid Zg ⊂ Rn. Finally, we use the LERP extension φd∗d of φd∗d to approximate φd∗ at the required point fs(x) for each x ∈ Xd, and compute T̂ dJd in (6f) as an approximation of T̂ J in (4c). With these approximations, we introduce the discrete CDP (d-CDP) operator as follows
εd(x) := γ · ∑
w∈Wd p(w) · J̃d(x+ w), x ∈ Xd, (6a)
εd∗d(y) = max x∈Xd
{ ⟨x, y⟩ − εd(x) } , y ∈ Yg, (6b)
Cd∗di (v) = max u∈Ud
{ ⟨u, v⟩ − Cdi (u) } , v ∈ Vg, (6c)
φd(y) := Cd∗di (−B ⊤y) + εd∗d(y), y ∈ Yg, (6d) φd∗d(z) = max y∈Yg { ⟨y, z⟩ − φd(y) } , z ∈ Zg, (6e)
T̂ dJd(x) := Cs(x) + φd∗d ( fs(x) ) , x ∈ Xd. (6f)
The proper construction of the grids Yg, Vg, and Zg will be discussed in Section 3.4.

heading: 3.3 Analysis of ConjVI algorithm
text: 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 1 (consult the Notations in Section 1): (i) The grid Vg is constructed such that co(Vgsub) ⊇ L(Cdi ). (ii) The grid Zg is constructed such that co(Zg) ⊇ fs ( Xd ) . (iii) The construction of Yg, Vg, and Zg requires at most
O(X + U) operations. The cardinality of the grids Yg and Zg (resp. Vg) in each dimension is the same as that of Xd (resp. Ud) in that dimension so that Y,Z = 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 infinity norm; (ii) Given a function J : X → R and its discretization Jd : Xd → R, we have ∥J − J̃d∥∞ ≤ ee for some constant ee ≥ 0.
Our first result concerns the contractiveness of the d-CDP operator.
Theorem 3.9 (Convergence). Let Assumptions 3.7-(ii) and 3.8-(i) hold. Then, the d-CDP operator (6) is γ-contractive in the infinity-norm.
The preceding theorem implies that the approximate ConjVI Algorithm 1 is indeed convergent given that the required conditions are satisfied. In particular, for deterministic dynamics, co(Zg) ⊇ fs ( Xd ) is sufficient for Algorithm 1 to be convergent. We next consider the complexity of our algorithm.
Theorem 3.10 (Complexity). Let Assumption 3.7-(iii) 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 Algorithm 1 are of O(X + U) and Õ(XWE), respectively.
The requirements of Assumption 3.7-(iii) will be discussed in Section 3.4. Recall that each iteration of VI (in primal domain) has a complexity of O(XUWE), where E denotes the complexity of the extension operation used in (2). This observation points to a basic characteristic of the proposed approach: ConjVI reduces the quadratic complexity of VI to a linear one by replacing the minimization operation in the primal domain with a simple addition in the conjugate domain. Hence, for problem class of Assumption 3.1, ConjVI is expected to lead to a reduction in the computational cost. We note that ConjVI, like VI and other approximation schemes that utilize discretization/abstraction of the continuous state and input spaces, still suffers from the so-called “curse of dimensionality.” This is because the sizes X and U of the discretizations increase exponentially with the dimensions n and m of the corresponding spaces. However, for ConjVI, this exponential increase is of rate max{m,n}, compared to the rate m+ n 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⋆ = T J⋆ : X → R and its discretization Jd⋆ : Xd → R, and let Assumption 3.8-(ii) hold for J⋆. Also, let Ĵd : Xd → R be the output of Algorithm 1. Then, ∥Ĵd − Jd⋆ ∥∞ ≤ γ(ee+et)+ed 1−γ , where ed = eu + ev + ex + ey + ez, and
eu = cu · dH(U,Ud), ev = cv · dH ( co(Vg),Vg ) , ex = cx · dH ( X,Xd ) ,
ey = cy ·maxx∈Xd d ( ∂(J⋆ − Cs)(x),Yg ) , ez = cz · dH ( fs(Xd),Zg ) ,
(7)
with constants cu, cv, cx, cy, cz > 0 depending on the problem data.
Let us first note that Assumption 3.5 implies that the DP and CDP operators preserve convexity, and they both have the true optimal value function J⋆ as their fixed point (i.e., the duality gap is zero). Otherwise, the proposed scheme can suffer from large errors due to dualization. Moreover, Assumptions 3.7-(i,ii) on the grids Vg and Zg are required for bounding the error of approximate discrete conjugations using LERP in (6d) and (6f); 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) ee is due to the approximation of the value function using the extension operator [̃·]; (ii) et corresponds to the termination of the algorithm after a finite number of iterations; (iii) ed captures the error due to the discretization of the primal and dual state and input domains.

heading: 3.4 Construction of the grids
text: In this subsection, we provide specific guidelines for the construction of the grids Yg, Vg and Zg. 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 LLT. The presented guidelines aim to minimize the error terms in (7) while taking into account the properties laid out in Assumption 3.7. In particular, the schemes described below satisfy the requirements of Assumption 3.7-(iii).
Construction of Vg. Assumption 3.7-(i) and the error term ev in (7) suggest that we find the smallest input dual grid Vg such that co(Vgsub) ⊇ L(Cdi ). This latter condition essentially means that Vg must “more than cover the range of slope” of the function Cdi ; recall that L(Cdi ) = Πmj=1 [ L−j (C d i ),L − j (C d i ) ] , where L−j (C d i ) (resp. L + j (C d i )) is the minimum (resp. maximum) slope of Cdi along the j-th dimension. Hence, we need to compute/approximate L ± j (C d i ) for j = 1, . . . ,m. A conservative approximation is L−j (Ci) = min ∂Ci/∂uj and L + j (Ci) = max ∂Ci/∂uj , assuming Ci is differentiable. Alternatively, we can directly use the discrete input cost Cdi for computing L±j (C d i ). In particular, if the domain Ud = Ug = Πmj=1U g j of C d i is grid-like and Ci is convex, we can take L−j (C d i ) (resp. L + j (C d i )) to be the minimum first forward difference (resp. maximum last backward difference) of Cdi along the j-th dimension (this scheme requires O(U) operations). Having L±j (C d i ) at our disposal, we can then construct V g sub = Π m j=1V g subj
such that, in each dimension j, Vgsubj is uniform and has the same cardinality as U g j , and co(V g subj ) = [ L−j (C d i ),L + j (C d i ) ] . Finally, we construct Vg by extending Vgsub uniformly in each dimension (by adding a smaller and a larger element to Vgsub in each dimension). Construction of Zg. According to Assumption 3.7-(ii), the grid Zg must be constructed such that co(Zg) ⊇ fs ( Xd ) . This can be simply done by finding the vertices of the smallest box that contains
the set fs ( Xd ) . Those vertices give the diameter of Zg in each dimension. We can then, for example, take Zg to be the uniform grid with the same cardinality as Yg in each dimension (so that Z = Y ). This way, dH ( fs(Xd),Zg ) ≤ dH ( co(Zg),Zg ) , and hence ez in (7) reduces by using finer grids Zg. This construction has a time complexity of O(X). Construction of Yg. Construction of the state dual grid Yg 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 Yg such that Yg ∩ ∂(J⋆ − Cs) ̸= ∅ for all x ∈ Xd so that ey = 0. Even if we had access to the optimal value function J⋆, satisfying such a condition could lead to a dual grid Yg ⊂ Rn of size O(Xn). Such a large size violates Assumption 3.7-(iii) on the size of Yg, and essentially renders the proposed algorithm impractical for dimensions n ≥ 2. A more practical condition is co(Yg)∩ ∂(J⋆ −Cs) ̸= ∅ for all x ∈ Xd so that maxx∈Xd d ( ∂(J⋆ − Cs)(x),Yg ) ≤ dH ( co(Yg),Yg ) , and hence ey reduces by using a finer grid Yg. The latter condition is satisfied if co(Yg) ⊇ L(J⋆ − Cs), i.e., if Yg “covers the range of slops” of (J⋆−Cs). Hence, we need to approximate the range of slopes of (J⋆−Cs). To this end, we first use the fact that J⋆ is the fixed point of DP operator (1) to approximate rng(J⋆−Cs) by R = rng(C d i )+γ·rng(C d s )
1−γ . We then construct the gird Y g = Πni=1Y g i such that, for each dimension i,
we have ±αR/∆iXd ∈ co(Y g i ), where ∆ i Xd denotes the diameter of the projection of X d on the i-th dimension, and α > 0 is a 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 Yg. Alternatively, we can construct Yg dynamically at each iteration in order to minimize the corresponding error in each application of the d-CDP operator given by ey = cy ·maxx∈Xd d ( ∂(T J−Cs)(x),Yg ) ; see Lemma A.6 and Proposition A.8 of the supplementary material. This means that the construction of Yg in Algorithm 1 is moved from line 1 to inside the iterations, after line 5. Similar to the static scheme described above, the aim here is to construct Yg such that co(Yg) ⊇ L(T J − Cs). Since we do not have access to T J (it is the output of the current iteration), we can again use the definition of the DP operator (1) to approximate rng(T J − Cs) by R = rng(Cdi ) + γ · rng(Jd) where Jd is the output of the previous iteration. We then construct the gird Yg = Πni=1Y g i such that, for each dimension i, we again have ±αR/∆iXd ∈ co(Y g i ). This construction has a one-time computational cost of O(U) for computing rng(Cdi ) and a per iteration computational cost of O(X) for computing rng(Jd). Notice, however, that under this dynamic construction, the error bound of Theorem 3.11 does not hold true. More importantly, with a dynamic grid Yg, there is no guarantee for ConjVI to converge.

heading: 4 Numerical simulations
text: 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
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 routine) provided in [23].
Example 1. We use the setup of [20, App. C.2.2] for the optimal control of a 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 ConjVI. We also consider the dynamic scheme for the construction of Yg in ConjVI (hereafter, referred to as ConjVI-d). Moreover, in each implementation of VI and ConjVI(-d), all of the involved discrete sets are 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 1. As shown in Figures 1a, both VI and ConjVI algorithms are convergent with a rate ≤ γ = 0.95. In particular, ConjVI 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 iteration, leads to a significant reduction in the running time of this algorithm compared to VI; see Figure 1b. Since we do 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 Jd computed using these algorithms. Figure 1c 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 dynamic scheme for construction of Yg, we see that ConjVI-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 2. 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, Z = N4 and U, V = N2. We note that we use multi-linear interpolation and extrapolation for the extension operator in (2) for VI. Due to the extrapolation, the extension operator is no longer non-expansive, and hence the convergence of VI is not guaranteed. On the other hand, since the dynamics is deterministic, there is no need for extension in ConjVI (the scaled expectation in (6a) in ConjVI reduces to the simple scaling εd = γ · Jd for deterministic dynamics), and hence the convergence of ConjVI only requires
co(Zg) ⊇ fs ( Xg ) . The results of simulations are shown in the bottom panels of Figure 1. 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 1d, we observe the non-monotone behavior of ConjVI-d. In this regard, recall that when the grid Yg 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-linear extrapolation operation used for extension is expansive.

heading: 5 Final remarks
text: In this paper, we proposed the ConjVI algorithm which reduces the time complexity of the VI algorithm from O(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 to apply the ConjVI algorithm: The dynamics must be of the form x+ = fs(x) +Bu+ w, and the stage cost C(x, u) = Cs(x) + Ci(u) must be separable. Moreover, since ConjVI essentially solves the dual problem, for non-convex problems, it 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 Yg 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 ConjVI 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 infinite-dimensional 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 4 (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 Yg must be controlled by the size X of the discrete state space Xd (Assumption 3.7-(iii)). 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 Yg partially addresses this issue by focusing on the range of slope of Jdk 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 a 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 to a convex function 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 [10, 26]. Moreover, the minimization over u again must be solved for each sample point in each iteration, while application of ConjVI avoids solving this minimization in each iteration. In this regard, let us note that ConjVI actually uses a convex approximation of the value function, which allows for application of a gradient-based algorithm for minimization over u within the ConjVI algorithm. However, such an algorithm has a per iteration complexity of O(XY ) = O(X2) , which is practically inefficient.

heading: Acknowledgments
text: This research is part of a project that has received funding from the European Research Council (ERC) under the grant TRUST-949796. The authors are also grateful to anonymous reviewers for their comments concerning the three remarks in Section 5.

Review summaries:
This paper studies the problem of dynamic programming with continuous state and action spaces.  The 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.  Regarding 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.

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).  The 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.

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.  The 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$. Thus, the slow $\inf$ operation can be replaced with a fast addition operation. The paper provides analyses of the convergence, time complexity, and error of the algorithm. When 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)$. Finally, the paper provides a numerical example implemented in MATLAB.

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.

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.