# Method Error and Domain Error when using NLsolve and Optim

**URL:** <https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418>\
**Category:** New to Julia\
**Created:** [July 16, 2018, 1:07pm UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418 "2018-07-16T13:07:40Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lippi](https://avatars.discourse-cdn.com/v4/letter/l/a3d4f5/32.png) [@Lippi](https://discourse.julialang.org/u/Lippi)\
**Post date:** [July 16, 2018, 1:07pm UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/1 "2018-07-16T13:07:40Z")

</div>

Hey all,

I am trying to solve a simple function with two variables (a profit function). I tried two ways: using Optim and hence function, gradient and hessian, as well as NLsolve, and hence only gradient and hessian. However, both either give me “DomainError: Exponentiation yielding a complex result requires a complex argument.” but mostly “MethodError: no method matching”. My first idea was that there is a type conversion within the iteration that let’s the solvers crash. Do you have any ideas???

```julia
function f!(F::AbstractArray, J::AbstractArray, x::AbstractArray;
             pa=1.0,
             pn=1.0,
             r=0.05,
             q=1.0,
             s=0.1341,
             Aa=4.04,
             γ=0.54,
             α1= 2.22*10.0^(-14),
             α2=0.75*10.0^(-2),
             β1= 0.32^(-2),
             β2= 0.17)
   
    F[1] = pa*((1 + β1*x[2]^β2)/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2))*(Aa*s)^(1-γ)*γ*x[1]^(γ-1) - (1+r)*q + q
    F[2] = pa*((β2*β1*x[2]^(β2-1)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) - (β2*β1^(2)*x[2]^(2*β2 -1) + β2*β1*x[2]^(β2-1)))/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2)^2)*(Aa*s)^(1-γ)*x[1]^γ - (1+r)*pn
    
    gprime = (β2-1)*β2*β1*x[2]^(β2-2)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) + (β1*β2)^(2)*x[2]^(2*β2 -2) - (2*β2 -1)*β2*β1^(2)*x[2]^(2*β2-2) - (β2-1)*β2*β1*x[2]^(β2-2)
    hprime = 2*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2)*β2*β1*x[2]^(β2-1)
    g = β2*β1*x[2]^(β2-1)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) - (β2*β1^(2)*x[2]^(2*β2 -1) + β2*β1*x[2]^(β2-1))
    h = (1 + β1*x[2]^β2 + α1*μ + α2*μ^2)^2
        
    J[1, 1] = pa*((1 + β1*x[2]^β2)/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2))*(Aa*s)^(1-γ)*γ*(γ-1)*x[1]^(γ-2)
    J[1, 2] = pa*(g/h)*(Aa*s)^(1-γ)*γ*x[1]^(γ-1)
    J[2, 1] = pa*(g/h)*(Aa*s)^(1-γ)*γ*x[1]^(γ-1)
    J[2, 2] = (gprime*h - g*hprime)/(h^2)*pa*(Aa*s)^(1-γ)*x[1]^γ
end

x0 = [0.0, 0.0]
solution = nlsolve(f!, x0)

```

Thanks a lot!!!

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [July 18, 2018, 9:12am UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/2 "2018-07-18T09:12:10Z")

</div>

Could you post a MWE with your error? The above gives me `μ not defined`, so you might have defined a global somewhere earlier in your code?

Your errors seem to indicate two problems:

(1) DomainError: this likely comes from exponentiating negative numbers with exponents between -1 and 1 (e.g. in `x[2]^β2`); you could get around this by constraining your search to positive `x`’s (if you’re looking at a profit function I’d assume these are factor inputs and therefore positive in any case?)

(2) MethodError: this is harder to diagnose without seeing all your actual function calls, but given that your only function seems to be `f!` this likely comes from writing the function definition in a way that is incompatible with the way the function is called in the solver. `Optim.optimize` for example would expect a function which takes one input array (in your case `x`), and returns the function value to minimise; depending on the method of `optimize` you’re calling you might also need to include the gradient (as a separate function!). Your function seems to return the `[2,2]` element of the Jacobian matrix, which is probably not what you’re actually trying to minimise?

---

<div class="post-metadata">

**Author:** ![Lippi](https://avatars.discourse-cdn.com/v4/letter/l/a3d4f5/32.png) [@Lippi](https://discourse.julialang.org/u/Lippi)\
**Post date:** [July 18, 2018, 10:22am UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/3 "2018-07-18T10:22:36Z")

</div>

> [@nilshg](#):
>
> MW

Hi,

first of all, one mistake I figured out is that the gradient and the hessian are not defined for zero values since I divide by the inputs and have negative exponents in gradient and hessian. This indeed would be solved by restricting the input values \>0.  
Let me give a different way I coded it using optim (so that I premultiplied function, gradient and hessian with -1 since I want to maximize)

```julia
const pa=1.0
const pn=1.0
const r=0.05
const q=1.0
const s=0.1341
const Aa=4.04
const γ=0.54
const α1= 2.22*10.0e-3
const α2=0.75*10.0e-2
const β1= 0.32e-2
const β2= 0.17
const μ = 0.99

function f(x)
    -1(pa*((1 + β1*x[2]^β2)/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2))*(Aa*s)^(1-γ)*x[1]^γ - (1+r)*(q*x[1] + pn*x[2]) + q*x[1])
end

function g!(G, x)
    
    G[1] = -1*(pa*((1 + β1*x[2]^β2)/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2))*(Aa*s)^(1-γ)*γ*x[1]^(γ-1) - (1+r)*q + q)
    G[2] = -1*(pa*((β2*β1*x[2]^(β2-1)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) - (β2*β1^(2)*x[2]^(2*β2 -1) + β2*β1*x[2]^(β2-1)))/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2)^2)*(Aa*s)^(1-γ)*x[1]^γ - (1+r)*pn)
    
end

function h!(H,x)
    
    gprime = (β2-1)*β2*β1*x[2]^(β2-2)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) + (β1*β2)^(2)*x[2]^(2*β2 -2) - (2*β2 -1)*β2*β1^(2)*x[2]^(2*β2-2) - (β2-1)*β2*β1*x[2]^(β2-2)
    hprime = 2*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2)*β2*β1*x[2]^(β2-1)
    g = β2*β1*x[2]^(β2-1)*(1 + β1*x[2]^β2 + α1*μ + α2*μ^2) - (β2*β1^(2)*x[2]^(2*β2 -1) + β2*β1*x[2]^(β2-1))
    h = (1 + β1*x[2]^β2 + α1*μ + α2*μ^2)^2
        
    H[1, 1] = -1*(pa*((1 + β1*x[2]^β2)/(1 + β1*x[2]^β2 + α1*μ + α2*μ^2))*(Aa*s)^(1-γ)*γ*(γ-1)*x[1]^(γ-2))
    H[1, 2] = -1*(pa*(g/h)*(Aa*s)^(1-γ)*γ*x[1]^(γ-1))
    H[2, 1] = -1*(pa*(g/h)*(Aa*s)^(1-γ)*γ*x[1]^(γ-1))
    H[2, 2] = -1*((gprime*h - g*hprime)/(h^2)*pa*(Aa*s)^(1-γ)*x[1]^γ)
    H
end

initial_x = [0.5, 0.5]
optimize(f, g!, h!, initial_x)

```

Thanks a lot

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [July 18, 2018, 10:59am UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/4 "2018-07-18T10:59:44Z")

</div>

Thanks, that makes more sense. What you probably want to do then is use box constrained optimization to set lower bounds, described in the [Optim.jl manual here](http://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/#box-constrianed-optimization). The function call would then be:

```julia
lower = [1e-5, 1e-5]
upper = [Inf, Inf]
initial_x = [0.5, 0.5]
inner_optimizer = GradientDescent()
results = optimize(f, g!, lower, upper, initial_x, Fminbox(inner_optimizer))

```

The problem with this is that for some reasons when I tried it just now, `Fminbox(inner_optimizer))` throws an error, which seems unrelated to your code (the error also shows up when simply running the example from the Optim.jl docs - I’ll file an issue)

One other option that I’ve used in my own thesis is the `NLopt` package which wraps the widely used [NLopt](https://nlopt.readthedocs.io/en/latest/) library. For this you would do something like:

```julia
const pa=1.0; const pn=1.0; const r=0.05; const q=1.0; const s=0.1341; const Aa=4.04
const γ=0.54; const α₁= 2.22*10.0e-3; const α₂=0.75*10.0e-2; const β₁= 0.32e-2;
const β₂= 0.17; const μ = 0.99

function obj(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = -1*(pa*((1 + β₁*x[2]^β₂)/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s)^(1-γ)*γ*x[1]^(γ-1) - (1+r)*q + q)
        grad[2] = -1*(pa*((β₂*β₁*x[2]^(β₂-1)*(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2) - (β₂*β₁^(2)*x[2]^(2*β₂ -1) + β₂*β₁*x[2]^(β₂-1)))/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2)^2)*(Aa*s)^(1-γ)*x[1]^γ - (1+r)*pn)
    end

    -(pa*((1 + β₁*x[2]^β₂)/(1 + β₁*x[2]^β₂ + α₁*μ + α₂*μ^2))*(Aa*s)^(1-γ)*x[1]^γ - (1+r)*(q*x[1] + pn*x[2]) + q*x[1])
end

using NLopt
opt = Opt(:LD_MMA, 2)
lower_bounds!(opt, [1e-5, 1e-5])
min_objective!(opt, obj)
(minf, minx, ret) = NLopt.optimize(opt, [0.5, 0.5])

```

(Note that I haven’t thought at all about your problem so the choice of bounds and/or optimisation algorithm here is likely suboptimal but it should hopefully get you going!)

Hope that helps!

P.S. send my regards to Sevi - he was one of the examiners in my viva (small world)!

---

<div class="post-metadata">

**Author:** ![Lippi](https://avatars.discourse-cdn.com/v4/letter/l/a3d4f5/32.png) [@Lippi](https://discourse.julialang.org/u/Lippi)\
**Post date:** [July 18, 2018, 11:13am UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/5 "2018-07-18T11:13:56Z")

</div>

Thanks a lot,

I will give it a try.

Yeah, what a small world 🙂 I will send him your regards!

---

<div class="post-metadata">

**Author:** ![Lippi](https://avatars.discourse-cdn.com/v4/letter/l/a3d4f5/32.png) [@Lippi](https://discourse.julialang.org/u/Lippi)\
**Post date:** [July 18, 2018, 1:42pm UTC](https://discourse.julialang.org/t/method-error-and-domain-error-when-using-nlsolve-and-optim/12418/6 "2018-07-18T13:42:14Z")

</div>

haaaa it worked out!

Thank you so much!