# Optim: gradient and Hessian of constraints

**URL:** <https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493>\
**Category:** Specific Domains\
**Tags:** optimization\
**Created:** [November 25, 2019, 1:57pm UTC](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493 "2019-11-25T13:57:14Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Alois](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alois/32/22558_2.png) [@Alois](https://discourse.julialang.org/u/Alois)\
**Post date:** [November 25, 2019, 1:57pm UTC](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493/1 "2019-11-25T13:57:14Z")

</div>

Hello, I am using Optim.jl to solve a constrained optimization problem.  
The gradient is not specified, so finite differences are the default.  
This works nicely for the objective, but not for the constraints.  
Thank you!

Here ist my mwe:

```julia
using Optim, Test

function fun(x) # objective
    (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 + (x[3]-x[1])^2
end

function con_c!(c, x)
    c[1]= x[1]^2 + x[2]^2 # 1st constraint
    c[2]= x[2]* sin(x[1])-x[1] # 2nd constraint
    c
end
function con_Jac!(J, x) # Jacobian of the constraint
    J[1,1] = 2*x[1]; J[1,2] = 2*x[2]
    J[2,1] = x[2]* cos(x[1])-1; J[2,2] = sin(x[1])
    J
end
function con_Hess!(h, x, λ) # Hessian of the constraint
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
    h
end

x_initial = [0.3, 0.2, 0.1]
df = TwiceDifferentiable(fun, x_initial) # here, automatic differentiation works
lx = [0, 0., 0]; ux = [1., 1., 1.]
lc = [-Inf, -3.]; uc = [0.5^2, 3.0]
# this works:
dfc= TwiceDifferentiableConstraints(con_c!, con_Jac!, con_Hess!, lx, ux, lc, uc)
# this is not working, but it should:
#dfc= TwiceDifferentiableConstraints(con_c!, lx, ux, lc, uc) # no method matching

@show res = optimize(df, dfc, x_initial, IPNewton())

```

---

<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 25, 2019, 4:15pm UTC](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493/2 "2019-11-25T16:15:33Z")

</div>

> [@Alois](#):
>
> TwiceDifferentiableConstraints

No one has implemented AD for those. It should be relatively simple though.

---

<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:** [November 27, 2019, 3:03am UTC](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493/3 "2019-11-27T03:03:21Z")

</div>

here is a function that implements autodiff on contraints:

```julia
function con_c!(c, x)
    c[1]= x[1]^2 + x[2]^2 # 1st constraint
    c[2]= x[2]* sin(x[1])-x[1] # 2nd constraint
    c
end

function constraint_derivatives(fn!, lx, ux, lc, uc)
nc = length(lc)
nx = length(lx)
cached_c = zeros(eltype(lc),nc)
cached_x = zeros(eltype(lx),nx)
cached_x2 = zeros(eltype(lx),nx*nc,nx)

function con_jac!(jac,x)
    ForwardDiff.jacobian!(jac,fn!,cached_c,x)
    jac
end

function con_jac2!(jac_vec,x)
    jac_mat = reshape(jac_vec,(nc,nx))
    ForwardDiff.jacobian!(jac_mat,fn!,zeros(eltype(x),nc),x)
    jac_vec
end

function con_hess!(hess, x, λ)
    ForwardDiff.jacobian!(cached_x2,con_jac2!,zeros(eltype(x),nc*nx),x)
    h = reshape(cached_x2,(nc,nx,nx))
    for i = 1:nc  
        hi = @view h[:,:,i]
        hess+=λ[i].*hi
    end
    hess
end
return TwiceDifferentiableConstraints(fn!, con_jac!, con_hess!,lx, ux, lc, uc)
end

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

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;
x0 = [0.25, 0.25]
lx = [-0.5, -0.5]; ux = [0.5, 0.5]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
dfc = constraint_derivatives(con_c!,lx,ux,lc,uc)
res = Optim.optimize(df, dfc, x0, IPNewton())

```

a problem is the allocation of arrays of duals that aren’t used when calculating jacobians and hessians, but it works

---

<div class="post-metadata">

**Author:** ![Alois](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alois/32/22558_2.png) [@Alois](https://discourse.julialang.org/u/Alois)\
**Post date:** [December 5, 2019, 3:43pm UTC](https://discourse.julialang.org/t/optim-gradient-and-hessian-of-constraints/31493/4 "2019-12-05T15:43:57Z")

</div>

Dear longemen3000,

thank you very much for this solution, this works well.  
Perhaps you want to contribute your solution to the general package?  
Please let me know if the solution is generally accepted.

Cheers,  
Alois