# 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[1] is assets to carry forward, x[2] is labor supply
opt = optimize( x -> -(utility(logit(x[2])*(w[t] + Ī¾[i])+ grid_A[s]*(1+r) - x[1], logit(x[2])) +
Ī²*interp_func(t+1)(x[1]) ),
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[2])
A1dict[i, s, t] = xstar[1]
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[1] is assets to carry forward, x[2] is labor supply
opt = optimize( x -> -(utility(logit(x[2])*(w[t] + Ī¾[i])+ grid_A[s]*(1+r) - x[1], logit(x[2])) +
Ī²*interp_func(t+1)(x[1]) ),
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[2])
A1dict[i, s, t] = xstar[1]
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