DP2019-24

Fast Value Iteration:
An Application ofLegendre-Fenchel Duality to a
Class of Deterministic Dynamic
Programming Problems in
Discrete Time
Ronaldo CARPIO
Takashi KAMIHIGASHI

December 17, 2019

Fast Value Iteration:
An Application of Legendre-Fenchel Duality to a Classof Deterministic Dynamic Programming Problems in
Discrete Time∗
Ronaldo Carpio†

Takashi Kamihigashi‡

November 28, 2019

Abstract
We propose an algorithm, which we call “Fast Value Iteration” (FVI), to computethe value function of a deterministic infinite-horizon dynamic programming problem
in discrete time. FVI is an efficient algorithm applicable to a class of multidimensional dynamic programming problems with concave return (or convex cost) functions
and linear constraints. In this algorithm, a sequence of functions is generated startingfrom the zero function by repeatedly applying a simple algebraic rule involving the
Legendre-Fenchel transform of the return function. The resulting sequence is guaranteed to converge, and the Legendre-Fenchel transform of the limiting function coincides
with the value function.
Keywords: Dynamic programming; Legendre-Fenchel transform; Bellman operator;
convex analysis∗

Part of this research was conducted while Ronaldo Carpio was visiting RIEB, Kobe University, asa Visiting Researcher in 2015. Earlier versions of this paper were presented at the SCE Conference on
Computing in Economics and Finance in Oslo, 2014, and the SAET Conference in Tokyo, 2014. Financial
support from the Japan Society for the Promotion of Science (JSPS KAKENHI Grant Number 15H05729)
is gratefully acknowledged.
†
School of Business and Finance, University of International Business and Economics, Beijing 100029,
China. Email: rncarpio@yahoo.com‡
RIEB and CCSS, Kobe University, 2-1 Rokkodai, Nada, Kobe 657-8501 Japan.
Email:
tkamihig@rieb.kobe-u.ac.jp.

1

Introduction

It has been known since Bellman and Karush (1962a, 1962b, 1963a, 1963b) that LegendreFenchel (LF) duality (Fenchel 1949) can be utilized to solve finite-horizon dynamic programming (DP) problems in discrete time. Although there have been subsequent applications ofLegendre-Fenchel (LF) duality to DP (e.g, Morin and Esogbue 1974, Klein 1990, Esogbue
and Ahn 1990, Klein and Morin 1991, Akian and Gaubert 2003), to our knowledge there has
been no serious attempt to exploit Legendre-Fenchel (LF) duality to develop an algorithm
to solve infinite-horizon DP problems.
In this paper we propose an algorithm, which we call “Fast Value Iteration” (FVI), tocompute the value function of a deterministic infinite-horizon DP problem in discrete time.
FVI is an efficient algorithm applicable to a class of multidimensional DP problems withconcave return functions (or convex cost functions) and linear constraints.
The FVI algorithm is an implementation of what we call the “dual Bellman operator,”
which is a simple algebraic rule involving the Legendre-Fenchel (LF) transform of the return
function. A sequence of functions generated by repeated application of the dual Bellman
operator is guaranteed to converge, and the Legendre-Fenchel (LF) transform of the limiting
function coincides with the value function. Involving no optimization, the dual Bellman
operator offers a dramatic computational advantage over standard computational methods
such as value iteration and policy iteration (e.g., Puterman 2005). We prove that the convergence properties of the iteration of the dual Bellman operator are identical to those ofvalue iteration when applied to a DP problem with a continuous, bounded, concave return
function and a linear constraint.
The rest of the paper is organized as follows. In Section 2 we review some basic concepts
from convex analysis and show some preliminary results. In Section 3 we present the general
DP framework used in our analysis. In Section 4 we apply Legendre-Fenchel (LF) dulaity toa DP problem with a continuous, bounded, concave return function and a linear constraint.
In Section 5 we present our numerical algorithm and compare its performance with that ofmodified policy iteration. In Section 7 we offer some concluding comments.

2

Preliminaries I: Convex Analysis

In this section we review some basic concepts from convex analysis and state some well-known
results. We also establish some preliminary results.
Let N ∈ N. Let R denote the extended real line; i.e., R = [−∞, ∞]. For f : RN → R,
we define f∗ , f ∗ : RN → R byf∗ (p) = inf {p| x − f (x)},∀p ∈ RN ,(1)

f ∗ (p) = sup {p| x − f (x)},∀p ∈ RN ,

(2)

x∈RN
x∈RN

where p and x are N × 1 vectors, and p| is the transpose of p. The functions f∗ and f ∗ arecalled the concave conjugate and convex conjugate of f , respectively.
1

It follows from (1) and (2) that for any functions f, g : RN
+ → R, we havef = −g⇒∀p ∈ RN , f∗ (p) = −g ∗ (−p) = −(−f )∗ (−p).

(3)

This allows us to translate any statement about g and g ∗ to the corresponding statement
about −g and (−g)∗ ; this is useful since most results in convex analysis deal with convex
functions and convex conjugates. In what follows we focus on concave functions and concave
conjugates, and by “conjugate,” we always mean “concave conjugate.” The biconjugate f∗∗of f is defined byf∗∗ = (f∗ )∗ .

(4)

A concave function f : RN → R is called proper if f (x) < ∞ for all x ∈ RN andf (x) > −∞ for at least one x ∈ RN .1 The effective domain of f is defined as
dom f = {x ∈ RN : f (x) > −∞}.

(5)

Let F be the set of proper, concave, upper semicontinuous functions from RN to R. Forf, g : RN → R ∪ {−∞}, the sup-convolution of f and g is defined as(f #g)(x) = sup {f (y) + g(x − y)}.

(6)

y∈RN

The following two lemmas collect the basic properties of conjugates we need later.
Lemma 1 (Rockafellar and Wets 2009, Theorems 11.1, 11.23).
(a) For any f ∈ F , we have f∗ ∈ F and f∗∗ = f .
(b) For any f, g ∈ F , we have (f #g)∗ = f∗ + g∗ .
0(c) Let N 0 ∈ N and u : RN → R. Let L be an N × N 0 matrix. Define (Lu) : RN → R by
(Lu)(x) = sup {u(c) : Lc = x},c∈RN 0∀x ∈ RN .

(7)

Then
(Lu)∗ (p) = u∗ (L| p),∀p ∈ RN .

(8)

Lemma 2 (Hiriart-Urruty, 1986, p. 484). Let f, g ∈ F be such that dom f = dom g = RN
+.
N
∗
∗
N
Suppose that both f and g are bounded on R+ . Then dom f = dom g = R+ , andsup |f (x) − g(x)| = sup |f∗ (p) − g∗ (p)|.
x∈RN
+

(9)

p∈RN
+

The following result is proved in the Appendix.
Lemma 3. Let A be an invertible N × N matrix. Let β ∈ R++ and S = A−1 /β. Letf, v : RN → R be such that f (x) = βv(Ax) for all x ∈ RN . Thenf∗ (p) = βv∗ (S | p).
1(10)This is a different concept from proper maps in topological spaces. For example, the function f : R → R
defined by f (x) = −x−1 for x > 0 and f (x) = −∞ for x ≤ 0 is proper as a concave function but not properin the topological sense.

2

3Preliminaries II: Dynamic Programming

In this section we present the general framework for dynamic programming used in our analysis, and show a standard result based on the contraction mapping theorem. Our exposition
here is based on Stokey and Lucas (1989) and Kamihigashi (2014a, 2014b).
Let N ∈ N and X ⊂ RN . Let x ∈ X. We need the following definitions from Stokey andLucas (1989, pp. 56–57). A correspondence Γ : X → 2X is called lower hemi-continuous at
x if Γ(x) 6= ∅ and if for any y ∈ Γ(x) and sequence {xn }∞n=0 with limn↑∞ xn = x, there exist∞
N ∈ N and a sequence {yn }n=N such that limn↑∞ yn = y and yn ∈ Γ(xn ) for all n ≥ N . A
correspondence Γ : X → 2X is called upper hemi-continuous at x if Γ(x) 6= ∅ and if for any∞
sequence {xn }∞n=0 with limn↑∞ xn = x and sequence {yn }n=0 with yn ∈ Γ(xn ) for all n ∈ Z+ ,∞there exists a subsequence {yni }i=0 of {yn } with limi↑∞ yni ∈ Γ(x).
A correspondence Γ : X → 2X is called continuous at x if it is both lower hemi-continuous
and upper hemi-continuous at x. It is called continuous if it is continuous at each x ∈ X.
Consider the following problem:
sup{xt+1 }∞t=0∞
X

β t r(xt , xt+1 )

(11)

t=0

s.t. x0 ∈ X given,∀t ∈ Z+ , xt+1 ∈ Γ(xt ).

(12)
(13)

In this section we maintain the following assumption.
Assumption 1. (i) β ∈ (0, 1). (ii) Γ : X → 2X is a nonempty, compact-valued, continuous
correspondence.2 (iii) X and gph Γ are convex sets, where gph Γ is the graph of Γ:gph Γ = {(x, y) ∈ RN × RN : y ∈ Γ(x)}.

(14)(iv) r : gph Γ → R is continuous, bounded, and concave.
Let v̂ : X → R be the value function of the problem (11)–(13); i.e., for x0 ∈ X we definev̂(x0 ) =sup{xt+1 }∞t=0∞
X

β t r(xt , xt+1 ) s.t. (13).

(15)

t=0

It is well-known that v̂ satisfies the optimality equation (see Kamihigashi 2008, 2014b):v̂(x) = sup {r(x, y) + βv̂(y)},∀x ∈ X.

(16)y∈Γ(x)

Thus v̂ is a fixed point of the Bellman operator B defined by(Bv)(x) = sup {r(x, y) + βv(y)},∀x ∈ X.

(17)y∈Γ(x)Let C (X) be the space of continuous, bounded, concave functions from X to R equippedwith the sup norm k · k. The following result is proved using standard arguments in the
Appendix. We claim no originality and state it here as a theorem only for reference purposes.
2

See Stokey and Lucas (1989, p. 57) for the definition of a continuous correspondence.

3

Theorem 1. Under Assumption 1, the following statements hold:(a) The Bellman operator B is a contraction on C (X) with modulus β; i.e., B maps C (X)into itself, and for any v, w ∈ C (X) we havekBv − Bwk ≤ βkv − wk.

(18)(b) B has a unique fixed point ṽ in C (X). Furthermore, for any v ∈ C (X),∀i ∈ N,kB i v − ṽk ≤ β i kv − ṽk.

(19)(c) We have ṽ = v̂, where v̂ is defined by (15).
Recent results on convergence and uniqueness that require neither continuity nor concavity can be found in Kamihigashi (2014b) and Kamihigashi et al. (2015).

4

The Dual Bellman Operator

In this section we introduce the “dual Bellman operator,” which traces the iterates of the
Bellman operator in a dual space for a special case of (11)–(13). In particular we consider
the following problem:max ∞{ct ,xt+1 }t=0∞X

β t u(ct )

(20)

t=0

s.t. x0 ∈ RN
+ given,∀t ∈ Z+ , xt+1 = Axt − Dct ,0N
ct ∈ RN+ , xt+1 ∈ R+ ,

(21)
(22)
(23)

where A is an N × N matrix, D is an N × N 0 matrix with N 0 ∈ N, and ct is a N 0 × 1 vector.
Throughout this section we maintain the following assumption.
0

Assumption 2. (i) β ∈ (0, 1). (ii) u : RN
+ → R is continuous, bounded, and concave.
N(iii) A is a nonnegative monotone matrix (i.e., Ax ∈ RN
+ ⇒ x ∈ R+ ). (iv) D is a nonzero
nonnegative matrix.
It is well-known (e.g., Berman and Plemmons 1994, p. 137) that a square matrix is monotone if and only if it is invertible and its inverse is nonnegative. Furthermore, a nonnegativesquare matrix is monotone if and only if it has exactly one nonzero element in each row andin each column (Kestelman 1973). Thus the latter property is equivalent to part (iii) above.
See Subsection 6.1 for a simple economic model satisfying this and other properties.

4

Under Assumption 2, the optimization problem (20)–(23) is a special case of (11)–(13)with
X = RN+,(24)
N0∀x ∈ X,Γ(x) = {y ∈ RN
+ : ∃c ∈ R+ , y = Ax − Dc},r(x, y) = max0 {u(c) : y = Ax − Dc}, ∀(x, y) ∈ gph Γ.

(25)
(26)

c∈RN
+

It is easy to see that Assumption 1 holds under Assumption 2 and (24)–(26).
Before proceeding, we introduce a standard convention for extending a function definedon a subset of Rn to the entire Rn . Given any function g : E → R with E ⊂ Rn , we extend
g to Rn by setting∀x ∈ Rn \ E.

g(x) = −∞,(27)

Note that in general, for any extended real-valued function f defined on a subset of Rn , wehave
sup f (x) = sup f (x),
x∈Rn

(28)

x∈dom fwhere the function f on the left-hand side is the version of f extended according to (27).
Letting L = A−1 D, we can express (22) as
Lct = xt − A−1 xt+1 .

(29)

In view of this and (26), we haver(x, y) = max0 {u(c) : Lc = x − A−1 y},∀(x, y) ∈ gph Γ.

(30)

c∈RN
+

By Assumption 2 and (24)–(27), the Bellman operator B defined by (17) can be writtenas(Bv)(x) = sup {(Lu)(x − z) + βv(Az)},∀x ∈ X,

(31)

z∈RN

where Lu is defined by (7) and (30) with z = A−1 y. The constraint y ∈ Γ(x) in (17) is
implicitly imposed in (31) by the effective domains of u and v, which require, respectively,
0that there exist c ∈ RN
+ with Lc = x − z and that y = Az ∈ X. Following the convention(27), we set(Bv)(x) = −∞,∀x ∈ RN \ X.

(32)

For any f : RN → R with dom f = RN
+ , we write f ∈ C (X) if f is continuous, bounded,
and concave on X. Since u is bounded, (Bv)(x) is well-defined for any x ∈ RN and v : X →R. In particular, for any v ∈ C (X), we have Bv ∈ C (X) by Theorem 1. The following resultshows that the Bellman operator B becomes a simple algebraic rule in the “dual” space ofconjugates.
5

Lemma 4. Let S = A−1 /β. For any v ∈ C (X) we have
(Bv)∗ (p) = u∗ (L| p) + βv∗ (S | p),∀p ∈ RN .

(33)

Proof. Let f (z) = βv(Az) for all z ∈ RN . Let g = Lu. We claim that
Bv = f #g.

(34)

It follows from (31) that (Bv)(x) = (f #g)(x) for all x ∈ X. It remains to show that
(Bv)(x) = (f #g)(x) for all x ∈ RN \ X or, equivalently, (f #g)(x) > −∞ ⇒ x ∈ X.
Let x ∈ RN with (f #g)(x) > −∞. Then there exists z ∈ RN with g(x − z) + f (z) =(Lu)(x − z) + βv(Az) > −∞; i.e.,x − z ∈ dom(Lu),Az ∈ dom v.

(35)

Since A−1 and D are nonnegative by Assumption 2, L is nonnegative; thus dom(Lu) ⊂ RN+.
Hence x − z ≥ 0 and Az ≥ 0 by (35). Since the latter inequality implies that z ≥ 0 bymonotonicity of A, it follows that x ≥ z ≥ 0; i.e, x ∈ X. We have verified (34).
By Lemma 1 we have (Bv)∗ = (f #g)∗ = f∗ + g∗ . Thus for any p ∈ RN+,
(Bv)∗ (p) = f∗ (p) + g∗ (p) = βv∗ (S | p) + u∗ (L| p),(36)

where the second equality uses Lemmas 3 and 1. Now (33) follows.
We call the mapping from v∗ to (Bv)∗ defined by (33) the dual Bellman operator B∗ ;more precisely, for any f : RN → R, we define B∗ f by(B∗ f )(p) = u∗ (L| p) + βf (S | p),∀p ∈ RN .

(37)

Using the dual Bellman operator B∗ , (33) can be written simply as(Bv)∗ = B∗ v∗ ,∀v ∈ C (X).

(38)

According to the convention (27), the domain (rather than the effective domain) of aNfunction f defined on RN+ can always be taken to be the entire R . Although this results in
no ambiguity when we evaluate supx f (x) (recall (28)), it causes ambiguity when we evaluate
kf k. For this reason we specify the definition of k · k as follows:
kf k = sup kf (x)k.

(39)

x∈RN
+

We use this definition of k · k for the rest of the paper. The following result establishes thebasic properties of the dual Bellman operator B∗ .
Theorem 2. For any v ∈ C (X), the following statements hold:6

(a) For any i ∈ N we have
(B i v)∗ = B∗i v∗ ,i

B v = (B∗i v∗ )∗ ,B∗i v∗ ∈ C (X),

(40)
(41)
(42)

where B∗i = (B∗ )i .
(b) The sequence {B∗i v∗ }i∈N converges uniformly to v̂∗ (the conjugate of the value functionv̂); in particular, for any i ∈ N we havekB∗i v∗ − v̂∗ k = kB i v − v̂k ≤ β i kv − v̂k = β i kv∗ − v̂∗ k.

(43)(c) We have v̂ = (v̂∗ )∗ .
Proof. (a) We first note that (40) implies (42) by Lemma 1 and Theorem 1(a). We show by
induction that (40) holds for all i ∈ N.
Note that for i = 1, (40) hods by (38). Suppose that (40) holds for i = n − 1 ∈ N:
(B n−1 v)∗ = B∗n−1 v∗ .

(44)

Let us consider the case i = n. We have B n−1 v ∈ C (X) by Theorem 1. Thus using (38) and(44), we obtain(B n v)∗ = (BB n−1 v)∗ = B∗ (B n−1 v)∗ = B∗ B∗n−1 v∗ = B∗n v∗ .

(45)

Hence (40) holds for i = n. Now by induction, (40) holds for all i ∈ N.
To see (41), let i ∈ N. Since B i v ∈ C (X) by Theorem 1, we have B i v = (B i v)∗∗ byLemma 1. Recalling (4), we have B i v = ((B i v)∗ )∗ = (B∗i v∗ )∗ , where the second equality uses
(40). We have verified (41).
(b) By Lemma 2 and (39), for any i ∈ N we havekB i v − v̂k = k(B i v)∗ − v̂∗ k = kB∗i v∗ − v̂∗ k,(46)

where the second equality uses (40). Thus the first equality in (43) follows; the second
equality follows similiarly. By Theorem 1 the inequality in (43) holds. As a consequence,{B∗i v∗ }i∈N converges uniformly to v̂∗ .
(c) Since v̂ ∈ C (X) by Theorem 1, we have v̂ = v̂∗∗ = (v̂∗ )∗ by Lemma 1. This completes
the proof of Theorem 2.
Figure 1 summarizes the results of Theorem 2. The vertical bidirectional arrows between
Bv and B∗ v∗ , B 2 v and B∗2 v∗ , etc, indicate that any intermediate result obtained by the Bellman operator B can be recovered through conjugacy from the corresponding result obtainedby the dual Bellman operator B∗ , and vice versa. This is formally expressed by statement (a)of Theorem 2. Statement (b) shows that both iterates {B i v} and {B∗i v∗ } converge exactly
the same way. In fact, as shown by Lemma 2, conjugacy preserves the sup norm between
any pair of functions in F whose effective domains are RN+ . The rightmost vertical arrow in
Figure 1 indicates that the value function v̂ can be obtained as the conjugate of the limit of{B∗i v∗ }, as shown in statement (c) of Theorem 2.
7

Figure 1: Bellman operator B and dual Bellman operator B ∗
B

v

5

B

B∗

B∗ v∗

B2v

B···B

B∗

B∗2 v∗

Biv

i↑∞

∗

∗

∗

∗

v∗

Bv

B∗ ···B∗B∗i v∗

v̂
∗

i↑∞

v̂∗

Fast Value Iteration

We exploit the relations expressed in Figure 1 to construct a numerical algorithm. The
upper horizontal arrows in Figure 1 illustrate the standard value iteration algorithm, whichapproximates the value function v̂ by successively computing Bv, B 2 v, B 3 v, · · · until convergence. The same result can be obtained by successively computing B∗ v∗ , B∗2 v∗ , B∗3 v∗ , · · ·until convergence and by computing the conjugate of the last iterate. Theorem 2(b) suggests that this alternative method can achieve convergence in the same number of steps asvalue iteration, but it is considerably faster since each step is a simple algebraic rule without
optimization; recall (37).
Algorithm 1, which we call “Fast Value Iteration,” implements this procedure with afinite number of grid points, using nearest-grid-point interpolation to approximate points
not on the grid. To be precise, we take n grid points p1 , . . . , pn in RN+ as given, and index
them by j ∈ J ≡ {1, . . . , n}. Recall from (42) that it suffices to consider the behavior ofNB∗i v∗ on X = RN
+ . We also take as given a function ρ : R+ → {p1 , . . . , pn } that maps eachNN
point p ∈ R+ to a nearest grid point. We define λ : R+ → J by ρ(p) = pλ(p) ; i.e., λ(p) is theindex of the grid point corresponding to p.
Algorithm 1 requires us to compute the conjugate of the return function u at the beginning
as well as the conjugate of the final iterate at the end. To compute these conjugates, weemploy the linear-time algorithm (linear in the number of grid points) presented in Lucet(1997), which computes the conjugate of a concave function on a box grid. Since the rateof convergence for {B∗i v∗ } is determined by β (as shown in Theorem 2(b)) and the numberof algebraic operations required for each grid point in each iteration of the “while” loop in
Algorithm 1 is independent of the number of grid points, it follows that FVI is a linear-timealgorithm.

8

Algorithm 1: Fast Value IterationN
let n grid points in RN+ be given by p1 , . . . , pn ∈ R+
initialize a, b, w : J → R (i.e., ∀j ∈ J, a(j), b(j), w(j) ∈ R)
initialize g : J → J (i.e., ∀j ∈ J, g(j) ∈ J)
compute u∗ on L| p1 , . . . , L| pn
for j = 1, . . . , nb(j) = 0w(j) = u∗ (L| pj )g(j) = λ(S | pj )fix  > 0
d = 2
while d > 
a=b
for j = 1, . . . , nb(j) = w(j) + βa(g(j))
d = maxj∈J {|a(j) − b(j)|}compute b∗
return b∗

6

Numerical Examples

6.1

The AK Model

As a simple example satisfying all our assumptions, consider the “AK model” of economic
growth used in macroeconomics (e.g. Acemoglu 2009):
max∞{ct }t=0∞X

β t u(ct )

(47)

t=0

s.t. x0 ≥ 0 given,∀t ∈ Z+ , ct + xt+1 = Axt ,ct , xt ≥ 0,

(48)
(49)
(50)

where ct is the consumption in period t, u(ct ) is the utility derived from consumption, xt isthe capital stock at the beginning of period t, and A > 0 is a constant. We assume that theutility function is given byu(c) =cq,qq ∈ (0, 1).

(51)

In Figure 2 we compare the solutions computed by FVI and value iteration using a logspaced grid of 500 points from 10−6 to 1000, and a convergence tolerance of 10−6 . The
parameter values used are β = 0.9, q = 0.2, and A = 0.9. For both algorithms, the initial
9

Figure 2: FVI vs. value iteration for the AK model

guess for the value function is u(x). As expected from our theoretical results, the solutions
produced by the two methods appear to be reasonably close to each other. Our experiencesuggests that accuracy can be improved with a larger number of grid points and a larger
domain for u(c). We make a more serious comparison in the next subsection.

6.2

Numerical Comparison

To illustrate the efficiency of FVI, we compare the performance of FVI with that of modified
policy iteration (MPI), which is a standard method to accelerate value iteration (Puterman
2005, Ch. 6.5). In what follows, we assume the following in (20)–(23):
1 0u(c1 , c2 ) = −(c1 − 10) − (c2 − 10) , β = 0.9, D =
.
0 1
2

2



(52)

Although u above is not bounded, it is bounded on any bounded region that contains
L| p1 , . . . , L| pn ; thus we can treat u as a bounded function for our purposes. Concerning
the matrix A, we consider two cases:

10




1 0
(a) A =
,
0 1




0 1.1
(b) A =
.
1 0

(53)

The grid points for MPI are evenly spread over [0, 20] × [0, 20]. For FVI, the same numberof grid points are evenly spread over a sufficiently large bounding box in the dual space.
We implement both FVI and MPI in Python, using the Scipy 0.13.3 package on a 2.40
GHz i7-3630QM Intel CPU. For MPI, we utilize C++ to find a policy that achieves the
maximum of the right-hand side of the Bellman equation (31) by brute-force grid search.
We use brute-force grid search because a discretized version of a concave function need not be
concave (see Murota 2003); we utilize C++ because brute-force grid search is unacceptably
slow in Python. The resulting policy is used to update the approximate value function 100
times, and the resulting approximate value function is used to find a new policy.
Table 1 shows the number of iterations and total CPU time for FVI and MPI to converge
to a tolerance of 10−5 . For each grid size, the final approximate value functions from FVI
and MPI are compared by computing, at each grid point, the absolute difference divided by
the largest absolute value of the MPI value function; we report the maximum and averagevalues of this difference over all grid points.
Panels (a) and (b) in Figure 3 plot the time to convergence of FVI and MPI against the
number of grid points using the data in Table 1. Panels (a’) and (b’) show the performance of
FVI for an extended range of grid point sizes. These plots indicate that FVI is a linear-timealgorithm, as discussed above. In terms of CPU time, FVI clearly has a dramatic advantage.

11

12
FVI

diff

MPI

diff

(b) MPI(a)

FVI

grid size
iterations
CPU time
iterations
CPU timemax
mean
iterations
CPU time
iterations
CPU time
max
mean

40x40
120
0.045
7
1.636
3.72E-03
1.20E-03
131
0.048
8
1.888
7.95E-03
2.63E-03

80x80
141
0.139
7
7.964
3.36E-03
1.27E-03
193
0.184
8
9.690
5.62E-03
1.19E-03

120x120178
0.299
834.901
2.31E-03
8.77E-04
192
0.311
8
35.887
5.01E-03
1.02E-03

160x160191
0.587
893.305
2.29E-03
6.91E-04
190
0.536
8
96.8484.74E-03
1.40E-03

200x200
209
0.954
8
214.673
1.84E-03
5.44E-04
205
0.943
8
220.004
4.60E-03
1.19E-03240x2402031.398
8
434.626
1.79E-03
7.20E-04
214
1.446
8
448.523
4.52E-03
1.07E-03

280x280219
2.086
8764.269
1.98E-03
6.34E-04
207
1.977
8
796.127
4.46E-03
1.19E-03

Table 1: Number of iterations to convergence, time to convergence in seconds, maximum relative difference, and averagerelative difference for FVI and MPI algorithms. Case (a) assumes (52) and (53)(a), while case (b) assumes (52) and(53)(b).

Figure 3: Time to convergence in seconds vs. number of grid points. Panels (a) & (a’) assume(52) & (53)(a), while panels (b) & (b’) assume (52) & (53)(b).

7

(a)

(b)

(a’)

(b’)

Concluding Comments

In this paper we proposed an algorithm called “Fast Value Iteration” (FVI) to compute thevalue function of a deterministic infinite-horizon dynamic programming problem in discrete
time. FVI is an efficient algorithm that offers a dramatic computational advantage for aclass of problems with concave return (or convex cost) functions and linear constraints.
The algorithm we presented is based on the theoretical results shown for continuous
state problems, but in practice, numerical errors are introduced through discretization andcomputation of conjugates. Although precise error estimates are yet to be established, ournumerical experiments suggest that the difference between the approximate value functions
computed using FVI and MPI, respectively, is rather insignificant.
In practice, one can combine FVI with other numerical methods to achieve a desired
combination of speed and accuracy. For example, to obtain essentially the same MPI valuefunction while economizing on time, one can apply FVI until convergence first and thenswitch to MPI. As in this algorithm, FVI can be used to quickly compute a good approximation of the value function.
In conclusion, we should point out that the theoretical results shown in Section 4 can beextended to problems with more general and nonlinear constraints using a general formulafor the conjugate of a composite function (Hiriart-Urruty 2006). New algorithms based on

13such an extension are left for future research.

A

Appendix

A.1

Proof of Lemma 3

Let p ∈ RN . Note that f∗ (p) = inf x∈RN {p| x − βv(Ax)}. Letting y = Ax and noticing thatx = A−1 y, we havef∗ (p) = inf {p| (A−1 y) − βv(y)}
y∈RN

(54)

= β inf {(p| A−1 /β)y − v(y)}(55)= β inf {(A−1 /β)| p)| y − v(y)}.

(56)

y∈RN

y∈RN

Now (10) follows.

A.2

Proof of Theorem 1

Let C(X) be the space of continuous bounded functions from X to R equipped with the supnorm k · k. Then statement (a) holds with C(X) replacing C (X) by Stokey and Lucas (1989,
Theorem 4.6). Thus if v ∈ C (X) ⊂ C(X), then Bv ∈ C(X); furthermore, Bv is concave by
a standard argument (e.g., Stokey and Lucas, 1989, p. 81). Thus B maps C (X) into itself.
Hence statement (a) holds. It is easy to see that C (X) equipped with the sup norm k · k
is a complete metric space; thus statement (b) follows by the contraction mapping theorem(Stokey and Lucas, 1989, p. 50). Finally, statement (c) holds by Stokey and Lucas (1989,
Theorem 4.3).

References
Acemoglu, D., 2009, Introduction to Modern Economic Growth, Princeton University Press, Princeton.
Akian, M., Gaubert, S., 2003, Spectral theorem for convex monotone homogeneous maps, andergodic control, Nonlinear Analysis 52, 637–679.
Bellman, R., Karush, W., 1962a, On the maximum transform and semigroups of transformations,
Bulletin of the American Mathematical Society 68, 516–518.
Bellman, R., Karush, W., 1962b, Mathematical programming and the maximum transform, SIAM
J. Appl. Math, 10, 550–567.
Bellman, R., Karush, W., 1963a, On the maximum transform, Journal of Mathematical Analysis
and Applications 6, 67–74.

14

Bellman, R., Karush, W., 1963b, Functional equations in the theory of dynamic programming XII:
an application of the maximum transform, Journal of Mathematical Analysis and Applications
6, 155–157.
Berman, A., Plemmons, R.J., 1994, Nonnegative Matrices in the Mathematical Sciences, SIAM,
Philadelphia.
Esogbue, A.O., Ahn, C.W., 1990, Computational experiments with a class of dynamic programmingalgorithms of higher dimensions, Computers Math. Applic. 19, 3–23.
Fenchel, W., 1949, On conjugate convex functions, Canadian Journal of Mathematics 1, 73–77.
Hiriart-Urruty, J.B., 1986, A general formula on the conjugate of the difference of functions, Canadian Mathematical Bulletin 29, 482–485.
Hiriart-Urruty, J.B., 2006, A note on the Legendre-Fenchel transform of convex composite functions,Nonsmooth Mechanics and Analysis, Springer US, 35-46.
Kamihigashi, T., 2008, On the principle of optimality for nonstationary deterministic dynamic
programming, International Journal of Economic Theory 4, 519–525.
Kamihigashi, T., 2014a, An order-theoretic approach to dynamic programming: an exposition,
Economic Theory Bulletin 2, 13–21.
Kamihigashi, T., 2014b, Elementary results on solutions to the Bellman equation of dynamic programming: existence, uniqueness, and convergence, Economic Theory, 1–23.
Kamihigashi, T., Reffett, K., Yao, M., 2015, An application of Kleene’s fixed point theorem todynamic programming: a note, International Journal of Economic Theory 11, 429–434.
Kestelman, H., 1973, Matrices with A ≥ 0 and A−1 ≥ 0, American Mathematical Monthly E2379
1059–1060.
Klein, C.M., 1990, Conjugate duality and its implications in dynamic programming, Math. Comput.
Modelling 14, 151–154.
Klein, C.M., Morin, T.L., 1991, Conjugate duality and the curse of dimensionality, European
Journal of Operational Research 50, 220–228.
Lucet, Y., 1997, Faster than the fast Legendre transform, the linear-time Legendre transform,
Numerical Algorithms 16, 171–185.
Morin, T.L., Esogbue, A.M.O., 1974, The imbedded state space approach to reducing dimensionality in dynamic programs of higher dimensions, Journal of Mathematical Analysis andApplications 48, 801–810.
Murota, K., 2003, Discrete Convex Analysis, SIAM, Philadelphia.
Puterman, M.L., 2005, Markov Decision Processes: Discrete Stochastic Dynamic Programming.
John Wiley & Sons, Hoboken, New Jersey.
Rockafellar, R.T., Wets, R.J.B., 2009, Variational Analysis, Springer, Dordrecht.
Stokey, N., Lucas, R.E. Jr., 1989, Recursive Methods in Economic Dynamics, Harvard University
Press, Cambridge MA.

15