Help speeding up finite horizon DP problem

amrods · April 5, 2019, 12:35pm

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

john.stachurski · April 5, 2019, 6:17pm

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.

Please read this lecture: https://lectures.quantecon.org/jl/need_for_speed.html

amrods · April 5, 2019, 6:38pm

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)

jlperla · June 12, 2019, 5:01pm

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.

Topic		Replies	Views
Efficient value function iteration in Julia Lectures	9	2737	January 24, 2022
Example for finite horizon DP problem	3	1491	February 6, 2019
More than 3 dimensions in Discrete Dynamic Programming	3	1404	May 21, 2019
Clarification on Bellman and the need for Linear Interp Lectures	3	864	November 2, 2018
Dynamic Programming Squared - Value iteration Lectures	3	1418	November 14, 2019

Help speeding up finite horizon DP problem

Related Topics