I’m trying to do maximum likelihood estimation with constraints:
p[1] + p[2] + p[3] + ... == 1
p[1] * a[1] / (1 - b[1]) + p[2] * a[2] / (1 - b[2]) + ... < 1
However, Optim seems to ignore the inequality constraint.
Constraints code
I set up both constraints like this:
function constr1(θ::AV)
par = get_params(θ)
sum(par[:p]) - 1 # == 0
end
function constr2(θ::AV)
par = get_params(θ)
p, α, β = par[:p], par[:α], par[:β]
sum(p .* α ./ (1 .- β)) - 1 # <= 0
end
constr(θ::AV) = [constr1(θ), constr2(θ)]
constr!(c::AV, θ::AV)::AV = (c .= constr(θ); c)
constr_jac!(J::AM, θ::AV)::AM = ForwardDiff.jacobian!(J, constr, θ)
function constr_hess!(H::AM, θ::AV, λ::AV)::AM
myH = similar(H)
ForwardDiff.hessian!(myH, constr1, θ)
H .+= λ[1] .* myH
ForwardDiff.hessian!(myH, constr2, θ)
H .+= λ[2] .* myH
H
end
Here, get_params extracts parameter vectors p, μ, ω, α, β from the parameter vector θ.
Setting up optimization problem
K = 2
# Initial guess
θ0 = [
ones(K) ./ K; #p
zeros(K); # μ
rand(K); # ω
zeros(K) .+ 1e-5; # α
zeros(K) .+ 1e-5
]
objective(θ) = D(θ, 1.0, data)
# Lower bounds
θ_lo = [
zeros(K); # p
fill(-Inf, K); # μ
zeros(K); # ω
zeros(K); zeros(K) # α, β
]
# Upper bounds
θ_hi = [
ones(K);
fill(Inf, K);
fill(Inf, K);
fill(Inf, K); fill(1., K); # α, β
]
# constr1 == 0
# -Inf <= constr2 <= 0
constr_lo = [0, -Inf]
constr_hi = [0, 0.]
dfc = TwiceDifferentiableConstraints(
constr!, constr_jac!, constr_hess!,
θ_lo, θ_hi, constr_lo, constr_hi
)
optimize(objective, dfc, θ0, IPNewton(), autodiff=:forward)
As far as I understand, constr(θ) must return a vector of two elements (because there are two constraints), where:
- the first element is exactly zero, since the constraint is
== 0 - the second element must be negative, since the constraint is
... - 1 < 0
Indeed, constr(θ0) == [0.0, -0.999989999899999], so the initial guess seems fine.
Problem
However, optimization fails after some steps. If I @show constr(ForwardDiff.value.(θ)) inside the objective function, I get this:
constr(ForwardDiff.value.(θ)) = [4.440892098500626e-16, 0.06428113096769317]
…where the second element is not negative, so the constraint is not satisfied.
What am I doing wrong with the inequality constraint?
A reproducible example can be found here: https://gist.github.com/ForceBru/847d5b214c2102547f056dcb1956cdae