Julia 1.6 versus 1.5.4

General Usage
1

I find that Julia 1.5.4 is faster in some use cases than 1.6. Do others experience this? What is the explanation.

Below is a piece of code that illustrates this. It finds the value function for an optimal growth problem—a canonical model in economics.

using Plots
using BenchmarkTools

function optimal_growth_model()
    σ = 1.5
    δ = 0.1
    β = 0.95
    α = 0.30
    ks = ((1-β*(1-δ))/(α*β))^(1/(α-1))
    csy = 1-α*β*δ/(1-β*(1-δ))
    Δ = 0.9
    kmin = (1-Δ)*ks
    kmax = (1+Δ)*ks
    Nk = 1000
    ∂k = (kmax-kmin)/(Nk-1)
    k = kmin:∂k:kmax
    v′ = zeros(Nk)
    i_k′ = zeros(Int64,Nk)
    v′ = (((csy*k.^α)).^(1-σ) .- 1) ./ ((1-σ)*(1-β))
    #v = Array{Float64}(undef,Nk)
    v = zeros(Nk)

    iter, crit, tol = 1,1, 1e-6

    while crit>tol
        for i in eachindex(1:Nk)
            # start by narrowing the interval on which to look for the solution
            i_min = Int64(max(ceil(((1-δ)*k[i] - kmin)/∂k)+1,1))
            i_max = Int64(min(floor((k[i]^α+(1-δ)*k[i]-kmin)/∂k)+1,Nk))

            # find the optimal consumption profile for the much narrow interval
            #note that there is no need for a second-loop here as the operation is vectorized
            c = k[i]^α + (1-δ)*k[i] .- k[i_min:i_max]
            u = (c.^(1-σ) .-1)/(1-σ)

            # for the narrow interval now find the consumption choice that minimizes the consumption profile
            #(v[i],i_tmp) = findmax(u+β*v′[i_min:i_max],dims=1) #this is nK×1 vector. Operate/vary on rows
            res = findmax(u+β*v′[i_min:i_max],dims=1) #this is nK×1 vector. Operate/vary on rows

            # extract both the maximizer and the index that we shall later use to recover the optimal consumption profile
            v[i] = res[1][1]    # we index multiple times to make types consistent
            i_tmp = res[2][1]    # we index multiple times to make types consistent
            i_k′[i] = i_min-1+i_tmp
        end
        crit = maximum(abs.(v-v′))
        v′ .= v
        iter+=1
    end
    k′ = k[i_k′]
    c = k.^α + (1-δ)*k - k′
    u = (c.^(1-σ) .- 1)/(1-σ)
    v = u/(1-β)
    return k′,k,c,u,v
end

function plot_figure()
    k′,k,c,u,v = optimal_growth_model()
    plt1 = plot(k,[k′ k],label=["k′ vs k" "k vs k"],leg=:bottomright,lw=2)
    plt2 = plot(k,c,label="c vs k",leg=:bottomright,lw=2)
    plt3 = plot(k,u,label="u vs k",leg=:bottomright,lw=2)
    plt4 = plot(k,v,label="v vs k",leg=:bottomright,lw=2)
    plot(plt1,plt2,plt3,plt4)
end

@time plot_figure()
1 Like
2

I can confirm this. To be a bit more specific I find

@btime optimal_growth_model()
# 1.6: 2.079 s (1477435 allocations: 2.11 GiB)
# 1.5: 977.718 ms (1477435 allocations: 2.11 GiB)

@time plot_figure()
# 1.6: 4.755011 seconds (first call) then 2.154710 seconds
# 1.5: 3.781330 seconds (first call) then 1.094716 seconds

So plotting isn’t slower, but the computation is.

1 Like
3

If I remove the two (small-union) type instabilities by adjusting the lines

c = k[i]^α .+ (1-δ)*k[i] .- k[i_min:i_max]

and

iter, crit, tol = 1,1.0, 1e-6

I find

julia> @btime optimal_growth_model() # Julia 1.5
  1.025 s (1688435 allocations: 2.53 GiB)

julia> @btime optimal_growth_model() # Julia 1.6
  1.860 s (1688435 allocations: 2.53 GiB)
2 Likes
4

Have you tried to run it with -O2 and -O3 flags? It looks that some O2 optimizations now belong to O3.

1 Like
5

Original code with -O3:

julia> @btime optimal_growth_model() # Julia 1.5
  796.140 ms (1477435 allocations: 2.11 GiB)

julia> @btime optimal_growth_model() # Julia 1.6
  2.040 s (1477435 allocations: 2.11 GiB)

@jovansam Unless someone resolves this here in the next few hours I’d file a bug report over at GitHub.

1 Like
6

I guess this is Performance regression with non-integer powers (v1.6 and newer) · Issue #39976 · JuliaLang/julia · GitHub.

4 Likes
7

So I looked into the non-integer powers issue with a MWE. 1.6.0 is actually competitive or even faster with non-integer powers in my case. So non-integer powers might not be the only cause of regression in this case.

See example code and results below:

using BenchmarkTools
function multiplies(x; α = 0.9)
    y = 0.0
    for i in eachindex(x)
        y = x[i] + y^α
    end
end

@btime multiplies(rand(1000000))

1.5.4
α = 0.9  104.750 ms (2 allocations: 7.63 MiB)
α = 1.0  2.679 ms (2 allocations: 7.63 MiB)
α = 1    1.717 ms (2 allocations: 7.63 MiB)
α = 1.1   9.061 ms (2 allocations: 7.63 MiB)

1.6.0
α = 0.9  46.383 ms (2 allocations: 7.63 MiB)
α = 1.0  2.872 ms (2 allocations: 7.63 MiB)
α = 1    1.462 ms (2 allocations: 7.63 MiB)
α = 1.1  8.651 ms (2 allocations: 7.63 MiB)

PS. Notice that I cannot set α=1 in the original optimal growth model.

8

If I @profile the code it shows that almost all time is spent in pow for me:

julia> using Profile

julia> @profile optimal_growth_model();

julia> Profile.print(noisefloor=2)
Overhead ╎ [+additional indent] Count File:Line; Function
...
2856╎    ╎    ╎    ╎    ╎    2857 @Base/math.jl:899; ^
...

What does it say for you? You can also use ProfileView to get a better view of the profile data. And compare 1.6 vs 1.5.

1 Like
9

With a @profile, I see something similar when comparing 1.5.4 to 1.6.0 indeed.

However, using other methods to solve the same model—but that still involve powers—does not result in this regression. See for example with the code below that uses an interpolation scheme in part of the code to solve the problem.

Here is where the cause of the regression becomes a little confusing for me. It seems more than one issues/thing are combining to create the regression.

using Plots, Dierckx

function optimal_growth_model()
    σ = 1.5
    δ = 0.1
    β = 0.95
    α = 0.30

    ks = ((1-β*(1-δ))/(α*β))^(1/(α-1))
    csy = 1-α*β*δ/(1-β*(1-δ))

    Δ = 0.9
    kmin = (1-Δ)*ks
    kmax = (1+Δ)*ks

    Nk = 20
    Nc = 1000

    ∂k = (kmax-kmin)/(Nk-1)
    k = kmin:∂k:kmax

    v = zeros(Nk)
    v = ((csy.*k.^α).^(1-σ) .- 1)/(1-σ)

    iter, crit, tol = 1,1, 1e-6

    cmin = 0.01
    cmax = kmax^α
    ∂c = (cmax-cmin)/(Nc-1)
    c = cmin:∂c:cmax


    util = (c.^(1-σ) .- 1)/(1-σ)

    k′ = zeros(Nc)
    v′ = zeros(Nc)

    t_v = zeros(Nk)
    i_v = zeros(Int64,Nk)

    while crit>tol
        for i in eachindex(1:Nk)
            k′ = @. k[i]^α + (1-δ)*k[i] - c
            v′ .= Spline1D(k,v; k=3, bc="extrapolate")(k′)
            res = findmax(util + β*v′,dims=1)
            t_v[i] = res[1][1]
            i_v[i] = res[2][1]
        end
    crit = 0
    crit = maximum(abs.(v-t_v))
    v .= t_v
    iter+=1
    end
    @show iter
    k,c[i_v], (k.^α + (1-δ)*k .- c[i_v]), ((c[i_v].^(1-σ) .- 1)/(1-σ)),v

end

function plot_figure()
    @time k,c,k′,u,v = optimal_growth_model()
    #plt1 = plot(k,[k′ k],label=["k′ vs k" "k vs k"],leg=:bottomright,lw=2)
    #plt2 = plot(k,c,label="c vs k",leg=:bottomright,lw=2)
    #plt3 = plot(k,u,label="u vs k",leg=:bottomright,lw=2)
    #plt4 = plot(k,v,label="v vs k",leg=:bottomright,lw=2)
    #plot(plt1,plt2,plt3,plt4)
end

plot_figure()