# Using ApproxFun for dynamic programming
**URL:** https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754
**Category:** Numerics
**Tags:** question
**Created:** [January 29, 2017, 10:50am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754 "2017-01-29T10:50:50Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [January 29, 2017, 10:50am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/1 "2017-01-29T10:50:50Z")
I am wondering how to use ApproxFun for dynamic programming. I know the methods themselves very well, and I can easily implement a Chebyshev or spline basis and do everything from scratch, but it would be nice to use existing, mature code (this is for teaching). So the question is about how to make use of the framework of `ApproxFun` to the largest extent (or whether I should do something else).
To make things concrete, consider the Bellman equation of a deterministic growth model:
```julia
V(k) = max_c [u(c) + β V(f(k)-c)]
```
where `u` and `f` are increasing concave functions (given) and `β` is a given constant. Usually one approximates `V` in some space (eg Chebyshev, splines), then either solves with
1. Value iteration: for each `k` on some grid, maximize the RHS, approximate it, and iterate. This is slow but guaranteed (if functions are “nice”), since `β<1` it is a contraction. In `ApproxFun`, I can do this by using constructors with specified gridpoints (what the documentation calls _Using ApproxFun for “manual” interpolation_).
2. Policy iteration: solve for `c` as a function of `k` on a grid, then given `c` everything else is linear, so one can use \ as described in the `ApproxFun` manual.
3. Projection: instead of the Bellman equation above in `V`, use the Euler equation
```julia
u'(c(k)) = β f'(k) u'(f(k)-c(k))
```
where `c(k)` is the unknown function, again approximated in some function space (eg Chebyshev, splines, linear interpolation). Decide on an objective (Galerkin, Least Squares, collocation), and minimize the disrepancy between RHS and LHS.
The problems I am running into: it is my impression that `ApproxFun` always tries go give me a very good approximation, with the assumption that the function approximated is relatively cheap. For dynamic programming, the function is relatively expensive, and since all procedures are iterative, it is wasteful to go for extremely good approximations anyway in each iteration step (until we are nearing convergence). So mostly I am using `ApproxFun` as a collection of basis functions.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [January 29, 2017, 11:50am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/2 "2017-01-29T11:50:44Z")
Do you have code that works? Having a profile of what’s too slow would help to answer.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 29, 2017, 2:05pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/3 "2017-01-29T14:05:40Z")
Are you asking for this?
[https://github.com/JuliaApproximation/ApproxFun.jl/issues/459](https://github.com/JuliaApproximation/ApproxFun.jl/issues/459)
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [January 29, 2017, 3:15pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/4 "2017-01-29T15:15:02Z")
A self-contained example for a single step of value iteration is below. `value_iteration` is the “naive” approach, while `value_iteration_grid` just evaluates on a grid. The former is very slow, the latter is OK (I did not optimize anything).
I think the question boils down to this: if I just want to use `ApproxFun` as a way to manage approximation using basis functions, then it seems that I am fighting against the natural (and very nice!) approach it was designed for (construction of basis matrices is tedious, etc). So maybe I should just write my own functions for parallel functionality.
```julia
using ApproxFun
using Parameters
using Optim
# we are solving
# V(k) = max_c u(c) + β V(f(k)-c) where f(k) = A k^α + (1-δ)*k
# with CRRA u(c)
# Euler equation:
# FOC u'(c) = β V'(f(k)-c)
# ENV V'(k) = β V'(f(k)-c) f'(k)
# impose steady state f(k)-c = k
# u'(c) = β u'(c) f'(k) => β f'(k) = 1
# then Aα k^(α-1) - δ =1 <=> k = ((Aα)/(1-δ))^(1/(1-α))
immutable GrowthModel
α::Float64 # capital share
δ::Float64 # depreciation
A::Float64 # productivity
σ::Float64 # CRRA utility coefficient
β::Float64 # discount factor
end
CRRA_utility(c,σ) = σ == 1 ? log(c) : (c^(1-σ))/(1-σ)
function steady_state(model::GrowthModel)
@unpack A, α, δ = model
((A*α)/(1-δ))^(1/(1-α))
end
function optimal_c_rhs(model::GrowthModel, k, V)
@unpack A, α, δ, β, σ = model
y = A*k^α + (1-δ)*k
opt = optimize(c -> CRRA_utility(c, σ) + β * V(y-c),
max(0.0, y-domain(V).b), y-max(0, domain(V).a))
Optim.converged(opt) || error("did not converge")
(Optim.minimizer(opt), Optim.minimum(opt))
end
function value_iteration(model::GrowthModel, V)
Fun(k->optimal_c_rhs(model, k, V)[2], domain(V))
end
function basis_matrix(S, p, M=length(p))
N = length(p)
V = Array{Float64}(N,M)
for k = 1:M
V[:,k] = Fun(S,[zeros(k-1);1]).(p)
end
V
end
function value_iteration_grid(model::GrowthModel, V, N)
d = domain(V)
k_grid = linspace(d.a, d.b, N)
V_grid = map(k->optimal_c_rhs(model, k, V)[2], k_grid)
S = Chebyshev(d)
Fun(S, basis_matrix(S, k_grid) \ V_grid)
end
model = GrowthModel(0.3, 0.05, 1, 1, 0.95)
# initialize
ss = steady_state(model)
V = Fun(identity, (0.5*ss)..(2*ss))
# single step -- this is very slow (127s, 5.27G allocated)
@time value_iteration(model, V)
# just evaluating in a grid, OK
V = value_iteration_grid(model, V, 10)
```
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [January 29, 2017, 3:24pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/5 "2017-01-29T15:24:27Z")
Possibly not, with
```julia
function value_iteration_nonadapt(model::GrowthModel, V, N)
Fun(k->optimal_c_rhs(model, k, V)[2], domain(V), N)
end
```
I get an oscillating behavior, not convergence.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [January 29, 2017, 7:40pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/6 "2017-01-29T19:40:18Z")
The constructor is very unoptimized at the moment. There are a few changes that may help:
1. Reuse the values. At the moment, it’s doubling the grid until convergence is verified, but it’s re-evaluating the function at the same grid over and over. This would be a fairly simple fix.
2. Convergence tolerance. There’s a hard coded `tol` In the constructor, that could easily be changed to a keyword. This is not going to make as much of a difference as you might think, as analytic functions converge exponentially fast.
3. Better convergence verification. Chebfun has a fairly complicated heuristic that could be copied. ApproxFun just naively looks at whether the tail is below the tolerance.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [January 30, 2017, 8:13am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/7 "2017-01-30T08:13:42Z")
Thanks for your help. I realized that my approach is by necessity more manual, and harder to automate with a high-level library like `ApproxFun`. What would really help me is the following interface:
1. representation for a basis with a given number of functions (eg the first `n` Chebyshev polynomials, tent functions for FEM with given nodes),
2. when applicable, its “natural” colocation points (eg Chebysev zeros, nodes of tent functions),
3. a way to construct an approximated function-like object `g` from a set of `(x,f(x))`, especially when `x` is the points in 2. above (which may not require matrix inversion, eg for Chebyshev),
4. quick evaluation of `g(x)`, and possibly `g'(x)` and `g''(x)` (for optimization), but going through AD is fine too.
I could do all this with a thin wrapper around `ApproxFun` constructs, I am trying to decide whether to do this or just write it from scratch.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 10, 2017, 1:58pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/8 "2017-03-10T13:58:17Z")
I have been thinking about this topic for a while (teaching a short course on numerical methods to economists in the meantime), and while I wrote my own code for the course for pedagogical purposes, I still think it would be nice to use `ApproxFun` in general. The usual requirements of numerical problems in economics boil down to the following:
1. Having a function basis `B`, and a set of coefficients `a`, calculate the approximated function `F(B, a)` at arbitrary points (in a grid), fast. Not always the same points, as the calculation involves optimization steps which depend on `a`. The calculation can be used to compute some residual function `g(F)`, which can be used for collocation, Galerkin, etc. `B` can be tensored univariate bases, or something more sophisticated to reduce the size of `a` for problems with many dimensions, eg complete polynomials, Smolyak.
2. Given a set of points on some grid, fit `a`. This is used in value iteration: projection methods are faster, but one needs a good initial guess.
3. Transform domains and function values. Domain transformations, such as `(-1,1)` to `(0,Inf)`, are common. Range transformations are sometimes needed, eg in a savings problem where the fraction of assets saved is in `(0,1)`.
4. Have basis functions which are crude, but preserve monotonicity and convexity, eg linear interpolation (eg linear splines). These have slow convergence, but are essential when the problem is non-smooth.
My understanding is that
1. `ApproxFun` provides a nice DSL for bases and representing functions. But the interface is very clever, trying to hide a lot from the user, and uses adaptive methods, so one would need to build a different API, along the lines of “just calculate, don’t re-approximate unless I ask for it”. But the core is basically there, just a thin layer of syntax is need.
2. Fitting can be done, and is described in the “manual” example of the docs. The API could be made easier though.
3. Transformations are trivial to add.
4. Crude linear bases are missing, the library focuses on spectral methods.
So I am unsure whether to contribute to ApproxFun or start something different. A lot of great code is already implemented in it, and duplicating it seems silly. But OTOH I am not sure PRs with a parallel API and weird bases would be appreciated.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [March 10, 2017, 10:34pm UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/9 "2017-03-10T22:34:44Z")
I’m having trouble understanding the first paragraph but it sounds like you just need function evaluation, which definitely works in ApproxFun. The rest is there too I think for the most part:
1. the domain `Ray()` represents functions on 0…Inf using a map to -1…1
2. `SplineSpace` represents general splines (not sure how much is implemented)
I’d be happy to take PRs, or help you with the code. Though I think some of what you mentioned could live in a separate package (for example, I’m not sure collocation should be built into ApproxFun itself)
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 11, 2017, 8:05am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/10 "2017-03-11T08:05:35Z")
Most of what I need is already in `ApproxFun`, what I need is a way to manage the bookkeeping. Knowing how to do the following (or that there is no support for it and I have to write a thin layer) would be useful:
1. How to denote a family of basis functions with a given number of elements in `ApproxFun`? Eg `Interval(0,1)` would give me a Chebyshev basis appropriately transformed, but is there a way to say “this basis uses the first 10 Chebyshev polynomials”?
2. Return the “natural” collocation points for a basis, eg Chebyshev nodes for Chebyshev polynomials, nodes + extra points for splines, etc.
3. How to provide coefficients for product spaces using the `Fun` constructor?
---
### Author:  [@ZacLN](https://discourse.julialang.org/u/ZacLN)
#### Post date: [March 11, 2017, 9:49am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/11 "2017-03-11T09:49:48Z")
As an alternative to ApproxFun I have a package [SparseGrids](https://github.com/ZacLN/SparseGrids) which I use for solving mid/high dimensional economic models. It uses linear or quadratic basis functions though the latter can be unstable if there are excessive kinks in the approximated functions.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [March 11, 2017, 11:41am UTC](https://discourse.julialang.org/t/using-approxfun-for-dynamic-programming/1754/12 "2017-03-11T11:41:02Z")
> [@Tamas\_Papp](#):
>
> How to denote a family of basis functions with a given number of elements in ApproxFun? Eg Interval(0,1) would give me a Chebyshev basis appropriately transformed, but is there a way to say “this basis uses the first 10 Chebyshev polynomials”?
There’s actually a space `SubSpace(::space,rng)` that can represent finite dimensional spaces. | is overriden as a shorthand:
```julia
julia> Chebyshev()|(1:3)
Chebyshev(【-1.0,1.0】)|1:3
julia> import ApproxFun: ∞
julia> Chebyshev()|(1:2:∞) # odd functions
Chebyshev(【-1.0,1.0】)|1:2:∞
```
Not much functionality is overloaded yet, so definite room for a PR.
> [@Tamas\_Papp](#):
>
> Return the “natural” collocation points for a basis, eg Chebyshev nodes for Chebyshev polynomials, nodes + extra points for splines, etc.
`points(sp,n)` is the right routine:
```julia
julia> points(Chebyshev(),10)
10-element Array{Float64,1}:
0.987688
0.891007
0.707107
0.45399
0.156434
-0.156434
-0.45399
-0.707107
-0.891007
-0.987688
```
> [@Tamas\_Papp](#):
>
> How to provide coefficients for product spaces using the Fun constructor?
Do you mean as a matrix of coefficients? There is support, look at `coefficientmatrix`, `fromtensor` and `totensor`.