# Optim and complex numbers: InexactError thrown

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 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(θ))

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}

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(θ))

end

x0_HH = [1.0]
lowerHH = [Para.lowerB]
upperHH = [1.9]
fct_HH(x) = HH_obj(x[1])
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


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

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?

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.

4. no need to dispatch on ::AbstractFloat.

1 Like

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

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