Optim and complex numbers: InexactError thrown

Specific Domains Optimization (Mathematical)
1

Good morning,

I encounter some issues when some complex numbers arise in my optimisation. I think it might be due to my function u(c) = c^gamma where sometimes c turns out negative (c being different expressions defined in code below and not a choice variable: c1_ND, c1_D, c2ND(θ), c2D).

My code is below and I apologize for its length. The reason I am not cutting it is that I am not sure where the error comes from and also that it works for some values of gamma. A working example of it is to limit the values of the parameter gamma can take by doing aversion=2 and riskaversion = collect(range(3.0, length=Para.aversion, stop=30.0)): then it works for those values in the loop on gamma (loop at the end).

My question is: is there a way to deal with complex numbers arising for some values of my parameters gamma with the currently used optim? (Or should i turn to JuMP and define a condition on the expressions of c ( c1_ND, c1_D, c2ND(θ), c2D) - which seems second best to me?)

Many thanks.

using Optim
using QuadGK
using QuantEcon
using Interpolations
using LineSearches

module Para 
    Re=1.3
    Rl=1.1
    E=2.0
    points=2
    lowerB=0.01
    aversion=3
end

    import ..Para

function solve(Gamma::AbstractFloat, centralized::Bool = true)
    
    u(c) = (c+0im).^Gamma
    
    function bk_obj(Sb::AbstractFloat, 
                    c_bar::AbstractFloat,
                    D::AbstractFloat)

        Lf= Para.Rl/(Para.Re- Para.Rl) * D * c_bar
        Sf= Para.E - D - Lf
        
        if centralized
                P_star= (Lf /Sb)
            else 
                P_star= Para.Re / Para.Rl
        end

        θhigh = (Para.Re * (1 - Sb) + Para.Re * (1 - Sb) * (P_star) - c_bar * (P_star) * D) / (c_bar * D * (Para.Re - (Lf / Sb)) )

        c1_ND = D * c_bar
        c1_D = (1 - Sb) + Lf
        c2ND(θ) = (((1 - Sb) - θ* c_bar * D) * Para.Rl + Sb * Para.Re + Lf * Para.Rl + Sf * Para.Re ) / (1 - θ)
        c2D(θ) = (1 - Sb) + Lf + ( (Sb + Sf) * Para.Re ) / (1 - θ)

        val1(θ) = u(θ * c1_ND + (1 - θ)* c2ND(θ)) 
        val2(θ) = u(θ * c1_D + (1 - θ)* c2D(θ)) 

        return -(quadgk(val1, 0.1, θhigh)[1] + quadgk(val2, θhigh, 0.999)[1])

    end #end of obj of bank

        grid = collect(range(Para.lowerB, length=Para.points, stop=Para.E))
        Sb = Vector{Float64}(undef,Para.points)
        cbar = Vector{Float64}(undef,Para.points)
        D_vec = Vector{Float64}(undef,Para.points)
        Res = Array{Float64}
        suboptimizer = GradientDescent(linesearch=LineSearches.Static())

        for (i,D) in enumerate(grid)

            bk(x) = bk_obj(x[1], x[2], D)
            D_vec[i]= D
            lower = [0.001, 0.001]
            upper = [D, 1.0]
            x0 = [0.01, 0.3]

                if i==1
            res = optimize(bk, lower, upper, x0, Fminbox(suboptimizer))
                    else 
            res = optimize(bk, lower, upper, Res, Fminbox(suboptimizer))
        end

            Res= Optim.minimizer(res)
            Sb[i]=Optim.minimizer(res)[1]
            cbar[i]=Optim.minimizer(res)[2]

         end

    D_row = (D_vec,) 

    Sb_int = interpolate(Sb, BSpline(Constant()))
    cbar_int = interpolate(cbar, BSpline(Constant()))
    Sb_func = extrapolate(Sb_int, Line())
    cbar_func = extrapolate(cbar_int, Line())

    function HH_obj(D::AbstractFloat) 

    Lf= Para.Rl/(Para.Re- Para.Rl) * D * cbar_func(D)
    Sf= Para.E - D - Lf
        
        if centralized
            P_star= (Lf /1 - Sb_func(D))
        else 
            P_star= Para.Re / Para.Rl
        end

    θhigh = (Para.Re * (1 - Sb_func(D)) + Para.Re * (1 - Sb_func(D)) * (P_star) - cbar_func(D) * (P_star) * D) / (cbar_func(D) * D * (Para.Re - (Lf / Sb_func(D))) )

    c1_ND = D * cbar_func(D)
    c1_D = (1 - Sb_func(D)) + Lf
    c2ND(θ) = (((1 - Sb_func(D)) - θ* cbar_func(D) * D) * Para.Rl + Sb_func(D) * Para.Re + Lf * Para.Rl + Sf * Para.Re ) / (1 - θ)
    c2D(θ) = (1 - Sb_func(D)) + Lf + ( (Sb_func(D) + Sf) * Para.Re ) / (1 - θ)

    val1(θ) = u(θ * c1_ND + (1 - θ)* c2ND(θ)) 
    val2(θ) = u(θ * c1_D + (1 - θ)* c2D(θ)) 

    return -(quadgk(val1, 0.1, θhigh)[1] + quadgk(val2, θhigh, 0.999)[1])    

    end

    x0_HH = [1.0]
    lowerHH = [Para.lowerB]
    upperHH = [1.9]
    fct_HH(x) = HH_obj(x[1])
    suboptimizer = GradientDescent(linesearch=LineSearches.Static())
    resHH= optimize(fct_HH, lowerHH, upperHH, x0_HH, Fminbox(suboptimizer))
    
    D_star_cent=Optim.minimizer(resHH)[1]
    
     if centralized
            D_star_cent=Optim.minimizer(resHH)[1]
            Sb_star_cent=Sb_func(D_star_cent)
            cbar_star_cent=cbar_func(D_star_cent)
            U_star_cent= -(HH_obj(D_star_cent))
            
            return  D_star_cent, Sb_star_cent, cbar_star_cent, U_star_cent
        
    else 
            D_star_dec=Optim.minimizer(resHH)[1]
            Sb_star_dec=Sb_func(D_star_dec)
            cbar_star_dec=cbar_func(D_star_dec)
            U_star_dec= -(HH_obj(D_star_dec))
            
            return D_star_dec, Sb_star_dec, cbar_star_dec, U_star_dec
    end
    
    
end

    riskaversion = collect(range(3, length=Para.aversion, stop=30))
    Gamma_vec = Vector{Float64}(undef,Para.aversion)
    diffU = Vector{Float64}(undef,Para.aversion)
    diffSb = Vector{Float64}(undef,Para.aversion)
    diffD = Vector{Float64}(undef,Para.aversion)

for (i,Gamma) in enumerate(riskaversion) 
    
    solve(Gamma, true)
    D_star_cent=solve(Gamma, true)[1]
    Sb_star_cent=solve(Gamma, true)[2]
    cbar_star_cent=solve(Gamma, true)[3]
    U_star_cent= solve(Gamma, true)[4]
    
    solve(Gamma, false)
    D_star_dec=solve(Gamma, false)[1]
    Sb_star_dec=solve(Gamma, false)[2]
    cbar_star_dec=solve(Gamma, false)[3]
    U_star_dec= solve(Gamma, false)[4]
    
    Gamma_vec[i]= Gamma
    diffU[i] = U_star_cent - U_star_dec 
    diffSb[i] = Sb_star_cent - Sb_star_dec 
    diffD[i] = D_star_cent - D_star_dec 
        
end

    return diffU, diffSb, diffD, Gamma_vec

and the error:

InexactError: Float64(Float64, 0.0 - 601724.3363405127im)

Stacktrace:
 [1] Type at .\complex.jl:37 [inlined]
 [2] convert(::Type{Float64}, ::Complex{Float64}) at .\number.jl:7
 [3] barrier_combined(::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Array{Float64,1}, ::getfield(Optim, Symbol("##55#61")){OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}},getfield(Optim, Symbol("##54#60")){Array{Float64,1},Array{Float64,1}}}, ::Base.RefValue{Float64}) at C:\Users\arquie\.julia\packages\Optim\wi7Vp\src\multivariate\solvers\constrained\fminbox.jl:65
 [4] (::getfield(Optim, Symbol("##57#63")){getfield(Optim, Symbol("##55#61")){OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}},getfield(Optim, Symbol("##54#60")){Array{Float64,1},Array{Float64,1}}}})(::Nothing, ::Array{Float64,1}) at C:\Users\arquie\.julia\packages\Optim\wi7Vp\src\multivariate\solvers\constrained\fminbox.jl:236
 [5] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Fminbox{GradientDescent{InitialPrevious{Float64},Static,Nothing,getfield(Optim, Symbol("##13#15"))},Float64,getfield(Optim, Symbol("##46#48"))}, ::Optim.Options{Float64,Nothing}) at C:\Users\arquie\.julia\packages\Optim\wi7Vp\src\multivariate\solvers\constrained\fminbox.jl:224
 [6] #optimize#53(::Bool, ::Symbol, ::Function, ::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Fminbox{GradientDescent{InitialPrevious{Float64},Static,Nothing,getfield(Optim, Symbol("##13#15"))},Float64,getfield(Optim, Symbol("##46#48"))}, ::Optim.Options{Float64,Nothing}) at C:\Users\arquie\.julia\packages\Optim\wi7Vp\src\multivariate\solvers\constrained\fminbox.jl:164
 [7] optimize at C:\Users\arquie\.julia\packages\Optim\wi7Vp\src\multivariate\solvers\constrained\fminbox.jl:163 [inlined] (repeats 2 times)
 [8] solve(::Float64, ::Bool) at .\In[16]:71
 [9] top-level scope at .\In[16]:155
2

First off, this is a rather unwieldly example, but I’m quite sure the problem is the mixing of complex and real variables. Why is it that you write

u(c) = (c+0im).^Gamma

? That +0im is surely giving you this error message.

1 Like
3

I had an error message telling me that complex number arise and to use that expression (it was in the error message and i copy pasted it). what should i use instead to accommodate possibly complex numbers?

4

Some general advice:

  1. don’t use a module for parameters. a struct or a NamedTuple should be more idiomatic.

  2. code your objective function separately, so that you can test it for various inputs outside optimization etc.

  3. break up your code to small functions. this should help you locate the error better.

  4. no need to dispatch on ::AbstractFloat.

1 Like
5

Are you sure you should be dealing with complex numbers in this problem (which looks like some economic model)? If you get DomainErrors because of negative bases, try formulating the problem in a way that excludes that, eg restricting the other parameters so that consumption is always (weakly) positive.

1 Like
6

Optim supports complex optimization (as in : optimization of a C^n → R function), but I doubt that’s what you want here. As @Tamas_Papp says, try to ensure the base is always positive.

1 Like
7

Thanks for all your answers. I will take into acocunt also the general advice… Thanks.

It is indeed an economic model. To ensure that the base is always positive, i should then use JuMP and a condition I imagine.

8

Since you are already using Optim, may be worth trying this Optim.jl.

9

Usually I find it nice to parametrize models in terms of a savings rate s \in [0,1]. Then c = sw, where w is beginning of period assets, and R(1-s)w or similar is assets next period. It is easy to extend this to other choices as a fraction of 1.

Then a box-constrained optimizer would work rather well, constraining everything within some [0,1]^n. That said, if it is a parametric problem (eg the solution to a functional equation, like a system of FOCs or value functions), you sould consider trying a very primitive Newton’s method first using the previous starting point, then see if it does not diverge for a couple of iterations, and then only switch to a local optimizer if it does.

If you outline the model in detail, perhaps you can get more specific help — it is harder to reverse-engineer it from the code.