# Optim and complex numbers: InexactError thrown

**URL:** <https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338>\
**Category:** Optimization (Mathematical)\
**Tags:** complex-numbers\
**Created:** [November 9, 2018, 11:08am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338 "2018-11-09T11:08:15Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![axelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axelle/32/4488_2.png) [@axelle](https://discourse.julialang.org/u/axelle)\
**Post date:** [November 9, 2018, 11:08am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/1 "2018-11-09T11:08:15Z")

</div>

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.

```julia
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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [November 9, 2018, 11:16am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/2 "2018-11-09T11:16:09Z")

</div>

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

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

```

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

---

<div class="post-metadata">

**Author:** ![axelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axelle/32/4488_2.png) [@axelle](https://discourse.julialang.org/u/axelle)\
**Post date:** [November 9, 2018, 11:17am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/3 "2018-11-09T11:17:23Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 9, 2018, 11:24am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/4 "2018-11-09T11:24:05Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 9, 2018, 11:26am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/5 "2018-11-09T11:26:17Z")

</div>

> [@axelle](#):
>
> what should i use instead to accommodate possibly complex numbers?

Are you sure you should be dealing with complex numbers in this problem (which looks like some economic model)? If you get `DomainError`s 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.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 9, 2018, 11:42am UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/6 "2018-11-09T11:42:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![axelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axelle/32/4488_2.png) [@axelle](https://discourse.julialang.org/u/axelle)\
**Post date:** [November 9, 2018, 12:28pm UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/7 "2018-11-09T12:28:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 9, 2018, 1:13pm UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/8 "2018-11-09T13:13:27Z")

</div>

Since you are already using Optim, may be worth trying this [Optim.jl](http://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 9, 2018, 1:56pm UTC](https://discourse.julialang.org/t/optim-and-complex-numbers-inexacterror-thrown/17338/9 "2018-11-09T13:56:00Z")

</div>

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.