# Help speeding up finite horizon DP problem

Crossposted in the Julia discourse, but trying here to get feedback on my technique.

I am struggling to speed up the computation of solutions for a finite horizon dynamic programming problem. I am using interpolation to avoid computing the value function outside of points in the grid. Right now it takes my computer around 7 minutes to solve this 5-period model.

Please have a look and let me know what I can improve. Any suggestions you may have will help me!

Thanks!

``````# VFI finite horizon
using Optim
using Distributions
using Random
using Interpolations

# transformation functions
# open interval (a, b)
function logit(x; a = 0, b = 1)
(b - a) * (exp(x)/(1 + exp(x))) + a
end

function utility(c::R, L::R) where R <: Real
if c <= 0 || L <= 0
return -Inf
else
return log(c) + log(1 - L)
end
end

const T = 5 # terminal period
const β = 0.95 # discount factor
const r = 0.05 # interest rate
const dist = LogNormal(0, 1) # distribution of ξ
const n = 1 # number of points of support of ξ

Random.seed!(1234)
const w = Vector{Float64}(undef, T) # exogenous wages

w .= (900 .+ 20.0 .* (1:T) .- 0.5 .* (1:T).^2)

const grid_A = -1000:10:8_000

const ξ = quantile.(dist, 0:1/n:prevfloat(1.0))

const Vdict = Array{Float64}(undef, (n, length(grid_A), T))
const Cdict = Array{Float64}(undef, (n, length(grid_A), T))
const Ldict = Array{Float64}(undef, (n, length(grid_A), T))
const A1dict = Array{Float64}(undef, (n, length(grid_A), T))
const convdict = Array{Bool}(undef, (n, length(grid_A), T))

# solve last period
for s in 1:length(grid_A)
for i in 1:n
opt = optimize( x -> -utility(w[T]*x + grid_A[s]*(1+r), x), 0.0, 1.0 )
Ldict[i, s, T] = Optim.minimizer(opt)
Cdict[i, s, T] = (w[T] + ξ[i])*Ldict[i, s, T] + grid_A[s]*(1+r)
A1dict[i, s, T] = 0.0
Vdict[i, s, T] = -Optim.minimum(opt)
end
end

interp_func(t) = LinearInterpolation(grid_A, sum(Vdict[i, :, t] for i in 1:n) ./ n, extrapolation_bc = Line())

# solve rest of periods
@time for t in T-1:-1:1
for s in 1:length(grid_A)
for i in 1:n
initial_x = [A1dict[i, s, t+1], 0.0]
# x is assets to carry forward, x is labor supply
opt = optimize( x -> -(utility(logit(x)*(w[t] + ξ[i])+ grid_A[s]*(1+r) - x, logit(x)) +
β*interp_func(t+1)(x) ),
initial_x, method = NewtonTrustRegion(), autodiff = :forward, iterations = 100_000)
Vdict[i, s, t] = -Optim.minimum(opt)
xstar = Optim.minimizer(opt)
Ldict[i, s, t] = logit(xstar)
A1dict[i, s, t] = xstar
Cdict[i, s, t] = Ldict[i, s, t] * (w[t] + ξ[i]) + grid_A[s] * (1+r) - A1dict[i, s, t]
convdict[i, s, t] = Optim.converged(opt)
end
end
println("period ", t, " finished")
end
``````

HI @amrods, I have to keep my answer short but the first thing I suggest you do is break your code up into functions. In particular, those loops where most of the execution time is occurring should be inside functions. The reason is that Julia does type inference at the function boundaries, and this is easier for the compiler when the functions are relatively small. They also need to be type stable, so one function can call another with prior knowledge of the type that will be returned.

1 Like

Thanks @john.stachurski, that was also suggested in the Julia discourse forum. I revised the program, see below.
However, I posted here as well to gather feedback about any “meta” tricks that economists employ to solve these models. For example, one can use linear algebra to solve a particular class of DP problems, which is something that programmers (like in the Julia forum) would overlook.

``````using Optim
using Interpolations

struct DP{Ty <: Real}
u::Function
T::Int
β::Ty
r::Ty
n::Int
w::Vector{Ty}
grid_A::Vector{Ty}
ξ::Vector{Ty}
end

function solvelast!(dp::DP{Ty}, Ldict::Array{Ty}, Cdict::Array{Ty}, A1dict::Array{Ty}, Vdict::Array{Ty}) where Ty <: Real
utility = dp.u
grid_A = dp.grid_A
n = dp.n
w = dp.w
ξ = dp.ξ
r = dp.r
T = dp.T
for s in 1:length(grid_A)
for i in 1:n
opt = optimize( x -> -utility(w[T]*x + grid_A[s]*(1+r), x), 0.0, 1.0 )
Ldict[i, s, T] = Optim.minimizer(opt)
Cdict[i, s, T] = (w[T] + ξ[i])*Ldict[i, s, T] + grid_A[s]*(1+r)
A1dict[i, s, T] = 0.0
Vdict[i, s, T] = -Optim.minimum(opt)
end
end
return Ldict, Cdict, A1dict, Vdict
end

# transformation functions
# open interval (a, b)
function logit(x; a = 0, b = 1)
(b - a) * (exp(x)/(1 + exp(x))) + a
end

# interpolation
interp_func(t) = LinearInterpolation(grid_A, sum(Vdict[i, :, t] for i in 1:n) ./ n, extrapolation_bc = Line())

function solverest!(dp::DP{Ty}, Ldict::Array{Ty}, Cdict::Array{Ty}, A1dict::Array{Ty}, Vdict::Array{Ty}, convdict::Array{Bool}; t0::Int=1) where Ty <: Real
utility = dp.u
grid_A = dp.grid_A
n = dp.n
w = dp.w
ξ = dp.ξ
r = dp.r
T = dp.T

for t in T-1:-1:t0
for s in 1:length(grid_A)
for i in 1:n
initial_x = [A1dict[i, s, t+1], 0.0]
# x is assets to carry forward, x is labor supply
opt = optimize( x -> -(utility(logit(x)*(w[t] + ξ[i])+ grid_A[s]*(1+r) - x, logit(x)) +
β*interp_func(t+1)(x) ),
initial_x, method = NewtonTrustRegion(), autodiff = :forward, iterations = 100_000)
Vdict[i, s, t] = -Optim.minimum(opt)
xstar = Optim.minimizer(opt)
Ldict[i, s, t] = logit(xstar)
A1dict[i, s, t] = xstar
Cdict[i, s, t] = Ldict[i, s, t] * (w[t] + ξ[i]) + grid_A[s] * (1+r) - A1dict[i, s, t]
convdict[i, s, t] = Optim.converged(opt)
end
end
println("period ", t, " finished")
end
return Ldict, Cdict, A1dict, Vdict, convdict
end
``````

I then instantiate the problem and solve it:

``````using Distributions

function u(c::R, L::R) where R <: Real
if c <= 0 || L <= 0
return -Inf
else
return log(c) + log(1 - L)
end
end

T = 5 # terminal period
β = 0.95 # discount factor
r = 0.05 # interest rate
dist = LogNormal(0, 1) # distribution of ξ
n = 1 # number of points of support of ξ

w = (900 .+ 20.0 .* (1:T) .- 0.5 .* (1:T).^2)

grid_A = collect(-1000:10.0:8_000)

ξ = quantile.(dist, 0:1/n:prevfloat(1.0))

Vdict = Array{Float64}(undef, (n, length(grid_A), T))
Cdict = Array{Float64}(undef, (n, length(grid_A), T))
Ldict = Array{Float64}(undef, (n, length(grid_A), T))
A1dict = Array{Float64}(undef, (n, length(grid_A), T))
convdict = Array{Bool}(undef, (n, length(grid_A), T))

dp = DP(u, T, β, r, n, w, grid_A, ξ)
solvelast!(dp, Ldict, Cdict, A1dict, Vdict)
solverest!(dp, Ldict, Cdict, A1dict, Vdict, convdict)``````

Sorry, I just noticed this post. The biggest source of performance issues is always making mistakes with the creation of `struct` where you end up with abstract types. Your code

Is almost generic, but you have the `u::Function` there, which is not a concrete type. You could try adding in another `F <: Function` generic type with `u::F``. It is tricky to tell if this will cause all your problems (i.e. it depends how it is passed across function boundaries), but it is usually the source of all woes.

This is why the notes recommend using `NamedTuples` instead of structs for holding parameter until you are an expert. See https://lectures.quantecon.org/jl/introduction_to_types.html#Syntax-for-Creating-Concrete-Types and https://lectures.quantecon.org/jl/fundamental_types.html#Tuples-and-Named-Tuples Even experts make concrete vs. abstract type problems when declaring structs.

I think that the best style for this sort of code (and the one with the maximum performance, when done right) is to never declare a single type in your code unless necessary. Let the compiler do its job.

• Unless you are dispatching to different algorithms based on the types, there is never a performance advantage of using a `struct` compared to a named tuple. Furthermore, since you never pick the type, named tuples are idiot-proof when it comes to the types.
• The performance of your `function solvelast!(dp::DP{Ty}, Ldict::Array{Ty}, Cdict::Array{Ty}, A1dict::Array{Ty}, Vdict::Array{Ty}, convdict::Array{Bool}; t0::Int=1) where Ty <: Real` is identical if you just go
``````function solvelast!(dp, Ldict, Cdict, A1dict, Vdict, convdict; t0 = 1)
``````
• There is an added advantage in never declaring any types in your functions or types: you can use generic types you hadn’t thought of (e.g. `StaticArrays`, `GPUArrays`, `BigFloat`, etc.) seamlessly, and you use auto-differentiation where appropriate. If you start declaring types all over the place, it gets very difficult to write generic code.
• See https://lectures.quantecon.org/jl/introduction_to_types.html#Tips-and-Tricks-for-Writing-Generic-Functions but beyond the “use named tuples for sets of parameters” and “don’t declare types for functions”, the biggest issue is how to create new vectors and values inside of your code. If you get into the habit of using `one`, `zero`, and `similar`, then it gets easy. `ones` and `zeros` require a little more work with an `eltype` because they have hidden floating points if you aren’t careful.
1 Like