Optim: gradient and Hessian of constraints

Specific Domains
1

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:

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())
1 Like
2

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

1 Like
3

here is a function that implements autodiff on contraints:

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

3 Likes
4

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

1 Like