# How to properly specify Jacobian & Hessian functions of inequality constraints (Optim)

**URL:** <https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870>\
**Category:** Optimization (Mathematical)\
**Tags:** optim, optimization\
**Created:** [January 24, 2021, 6:12pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870 "2021-01-24T18:12:34Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Post date:** [January 24, 2021, 6:12pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/1 "2021-01-24T18:12:34Z")

</div>

I’m trying to optimize an objective function with 19 variables in Optim with the following inequality constraints:

```julia
0 <= x[1]/3 - x[2] <= 1/3 
5 <= 1/x[3] + 1/x[4] <= 6

```

I’m trying to use either `IPNewton()` or `NewtonTrustRegion`, so I need to supply both a Jacobian and Hessian for the constraints. My question is: what is the correct way to write the Jacobian and Hessian functions?

I believe the constraint function would be

```julia
function con_c!(c,x)
   c[1] = x[1]/3 - x[2]          
   c[2] = 1/x[3] + 1/x[4]
   c
end

```

Would the Jacobian function be

```julia
function con_jacobian!(J,x)
   #first constraint:
   J[1,1] = 1/3
   J[1,2] = -1.0
  #second constraint:
   J[2,3] = -1/(x[3])^2
   J[2,4] = -1/(x[4])^2
   J
end

```

? (I assume all other indices of `J` are automatically set to zero?)

My main question: What would the Hessian function be? This is where I’m most confused. My understanding was that we take Hessians of scalar-valued functions. So do we have to enter _multiple_ Hessians, one for each constraint function (2 in my case)?

I’ve looked at the multiple constraints example given here [https://github.com/JuliaNLSolvers/ConstrainedOptim.jl](https://github.com/JuliaNLSolvers/ConstrainedOptim.jl), but I’m still confused. In the example, it looks like they are adding together two Hessian matrices…? Would greatly appreciate some help.

---

<div class="post-metadata">

**Author:** ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Post date:** [January 26, 2021, 7:36pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/2 "2021-01-26T19:36:27Z")

</div>

Just a note that I’ve posted this question on [Stack Overflow](https://stackoverflow.com/questions/65908057/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-i) because I haven’t received any replies here after 2 days.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [January 26, 2021, 8:33pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/3 "2021-01-26T20:33:11Z")

</div>

Here are the docs: [https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton\_basics/#multiple-constraints](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/#multiple-constraints)

I guess you want something like this, assuming I haven’t made a mistake.

```nohighlight
using Optim

f(x) = sum(x)

function f_grad!(g, x)
    fill!(g, 1.0)
    return g
end

function f_hess!(h, x)
    fill!(h, 0.0)
    return h
end

function c(y, x)
    y[1] = x[1] / 3 - x[2]
    y[2] = 1 / x[3] + 1 / x[4]
    return y
end

function c_jac!(J, x)
    fill!(J, 0.0)
    # c₁'(x) = [1 / 3, -1, 0, 0]
    J[1, 1] = 1 / 3
    J[1, 2] = -1
    # c₂'(x) = [0, 0, -1 / x[3]^2, -1 / x[4]^2]
    J[2, 3] = -1 / x[3]^2
    J[2, 4] = -1 / x[4]^2
    return J
end

function c_hess!(H, x, λ)
    # c₁''(x) = 0
    # c₂''(x) = Diag([0, 0, 2 / x[3]^3, 2 / x[4]^2])
    H[3, 3] += λ[2] * 2 / x[3]^3
    H[4, 4] += λ[2] * 2 / x[4]^3
    return H
end

x0 = [0.5, 0.1, 1.2, 0.2]
df = TwiceDifferentiable(f, f_grad!, f_hess!, x0)
dfc = TwiceDifferentiableConstraints(
    c, 
    c_jac!, 
    c_hess!, 
    [0.0, 0.0, 0.0, 0.0], 
    Float64[], 
    [0.0, 5.0], 
    [1 / 3, 6.0],
)
optimize(df, dfc, x0, IPNewton())

```

You could also write this in JuMP ([Introduction · JuMP](https://jump.dev/JuMP.jl/stable/)), which may be easier:

```nohighlight
using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, x[1:4] >= 0.00001)
@NLconstraint(model, 0 <= x[1] / 3 - x[2] <= 1 / 3)
@NLconstraint(model, 5 <= 1 / x[3] + 1 / x[4] <= 6)
@objective(model, Min, sum(x))

```

---

<div class="post-metadata">

**Author:** ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Post date:** [January 26, 2021, 11:45pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/4 "2021-01-26T23:45:39Z")

</div>

Thanks for your reply! I used your Hessian function code and the algorithm worked! Out of curiousity, I tried replacing `IPNewton()` with `NewtonTrustRegion()`, just to see how the results would compare. But I got this weird error:

```julia
StackOverflowError:

Stacktrace:
 [1] optimize(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{var"#con_c!#89",var"#con_jacobian!#90",var"#con_hess!#91",Float64}, ::Array{Float64,1}, ::NewtonTrustRegion{Float64}, ::Optim.Options{Float64,Nothing}; inplace::Bool, autodiff::Symbol) at C:\Users\Michael\.julia\packages\Optim\CK6Dn\src\multivariate\optimize\interface.jl:148
 [2] optimize(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{var"#con_c!#89",var"#con_jacobian!#90",var"#con_hess!#91",Float64}, ::Array{Float64,1}, ::NewtonTrustRegion{Float64}, ::Optim.Options{Float64,Nothing}) at C:\Users\Michael\.julia\packages\Optim\CK6Dn\src\multivariate\optimize\interface.jl:147
 ... (the last 2 lines are repeated 28980 more times)
 [57963] ModelFit(; ODE_model::Function, ODE_vars::Array{String,1}, vars_to_fit::Array{String,1}, paramIC_dict::Dict{Symbol,Any}, data_dict::Dict{Symbol,String}, f_calc_ICs::typeof(f_ICs), N0::Int64, fit_range::Array{String,1}, date_obs_last::String, forecast_until::String, date_format::String, norm::typeof(mean_abs_rel_error), norm_weights::Array{Int64,1}, norm_scale::Int64, integrator_options::Dict{Symbol,Any}, optimizer_options::Dict{Symbol,Any}) at .\In[60]:317
 [57964] top-level scope at In[61]:15

```

It did work, however, when I removed the inequality constraints. Does `NewtonTrustRegion()` not support inequality constraints? (I checked the docs but they only gave one very short example for this method, and it was without any constraints.)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [January 27, 2021, 12:13am UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/5 "2021-01-27T00:13:48Z")

</div>

I assume it doesn’t support constraints, but it should throw a nicer error. You should open an issue: [https://github.com/JuliaNLSolvers/Optim.jl/issues](https://github.com/JuliaNLSolvers/Optim.jl/issues)

---

<div class="post-metadata">

**Author:** ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Post date:** [January 27, 2021, 9:37pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/6 "2021-01-27T21:37:20Z")

</div>

Thanks, I just opened an issue [here](https://github.com/JuliaNLSolvers/Optim.jl/issues/903).

I also think the documentation could clearer about the capabilities of the different algorithms, particularly regarding the types of constraints that each supports/doesn’t support. It’d be nice to have a table like this: [https://www.mathworks.com/help/optim/ug/problems-handled-by-optimization-toolbox-functions.html](https://www.mathworks.com/help/optim/ug/problems-handled-by-optimization-toolbox-functions.html)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [January 27, 2021, 11:04pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/7 "2021-01-27T23:04:09Z")

</div>

Open an issue if you have suggestions to improve the documentation.

Even better, make a pull request editing the documentation. Here’s the algorithm choice section:  
[https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/user/algochoice.md](https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/user/algochoice.md)  
and here’s the home page:  
[https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/index.md](https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/index.md)

You can click the pencil icon to edit online. Follow the instructions at the bottom of the page to make a pull request.

---

<div class="post-metadata">

**Author:** ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)\
**Post date:** [January 28, 2021, 8:08pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/8 "2021-01-28T20:08:06Z")

</div>

Thanks for the suggestions. For now I’ll open an issue and include some broad suggestions for how I think the documentation could be improved. Since I’m still a Julia newbie, I don’t know enough about the Optim package to be able to actually revise the documentation…

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [January 28, 2021, 11:16pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/9 "2021-01-28T23:16:39Z")

</div>

i did a patch supporting AD on constraints, basically, the “hessian” here, is not a hessian, but an operation over a hessian, adding the lagrange multipliers. a simplified version of what the jacobian and hessian of constraints is doing can be found in this post:

> [@Optim: gradient and Hessian of constraints](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493/2):
>
> No one has implemented AD for those. It should be relatively simple though.

---

<div class="post-metadata">

**Author:** ![luluywang](https://avatars.discourse-cdn.com/v4/letter/l/ee7513/32.png) [@luluywang](https://discourse.julialang.org/u/luluywang)\
**Post date:** [April 24, 2022, 10:34pm UTC](https://discourse.julialang.org/t/how-to-properly-specify-jacobian-hessian-functions-of-inequality-constraints-optim/53870/10 "2022-04-24T22:34:46Z")

</div>

This has actually been implemented but poorly documented. I was tipped off to this by seeing the package author commenting somewhere that they fixed this, but given I can’t seem to find that post on Discourse it’s likely worth copying.

The key is that there is a parameter in the TwiceDifferentiableConstraint constructor that allows for autodiff, but unlike in TwiceDifferentiable it is not a keyword variable.

Here’s a MWE based on the IPNewton example

```julia
using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide
import ForwardDiff as Diff

fun(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
lx = [-0.5, -0.5]; ux = [0.5, 0.5]

function con2_c!(c, x)
    c[1] = x[1]^2 + x[2]^2 ## First constraint
    c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
    return c
end

df = TwiceDifferentiable(fun, x0, autodiff = :forwarddiff)
dfc = TwiceDifferentiableConstraints(con2_c!,
                                     lx, ux, lc, uc, :forwarddiff)
res = optimize(df, dfc, x0, IPNewton())
println(res)
println(res.minimizer)

```