# Optim.jl violates inequality constraint

**URL:** https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573
**Category:** Optimization (Mathematical)
**Tags:** optim
**Created:** [February 1, 2022, 12:31pm UTC](https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573 "2022-02-01T12:31:42Z")
**Posts on this page:** 4
**Page:** 1

<div class="post-metadata">

### Author: ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)
#### Post date: [February 1, 2022, 12:31pm UTC](https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573/1 "2022-02-01T12:31:42Z")

</div>

I’m trying to do maximum likelihood estimation with constraints:

```julia-auto
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:

```julia
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

```julia
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: [Optim.jl violates inequality constraint? · GitHub](https://gist.github.com/ForceBru/847d5b214c2102547f056dcb1956cdae)

---

<div class="post-metadata">

### Author: ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)
#### Post date: [February 4, 2022, 4:06pm UTC](https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573/2 "2022-02-04T16:06:18Z")

</div>

Any help? Please? Any pointers? Anything?

- I put `@assert`s everywhere, they all pass.
- I used [autodiff with constraints](https://github.com/JuliaNLSolvers/Optim.jl/issues/971#issuecomment-1029802689) - same error.
- I changed the order of constraints - this didn’t help.
- I changed the inequality constraint from `<= 0` to `>= 0` - it didn’t help
- I replaced broadcasting like `sum(p .* α ./ (1 .- β))` with generator expressions, like `sum(p[k] * α[k] / (1 - β[k]) for k in eachindex(p))`. I know [Symbolics.jl can’t differentiate broadcasted expressions](https://github.com/JuliaSymbolics/Symbolics.jl/issues/498), so maybe they’re the problem here… even though Optim uses ForwardDiff? Yet nothing changed after I did that.
- I tried a finite value instead of `Inf` for the inequality constraint’s upper bound - didn’t help.
- I tried multiple versions of Optim.jl: 1.6.1, 1.6.0, 1.5.0, 1.4.0 and 1.3.0 - same error everywhere

It’s like the inequality constraint isn’t there…

* * *

Here’s the updated code, littered with asserts and without manual derivative computation.

```julia
import Pkg
Pkg.activate(temp=true, io=devnull)
Pkg.add(name="Optim", version="1.6.1", io=devnull)
Pkg.add(["ForwardDiff", "Distributions"], io=devnull)
Pkg.status()

import Random
Random.seed!(0)

using Optim, ForwardDiff, Distributions

const AV = AbstractVector{T} where T;
const AM = AbstractMatrix{T} where T;

# ===== Functions to be optimized =====
function D_t(p::AV, μ::AV, var::AV, b::Real, X::Real)
    K = length(p)
    @assert K == 2

    p1 = sum(@. p^2 * sqrt(π / (b + var)))
    p2 = sum(
        @. p * sqrt(2π / (2b + var)) * exp(-(X - μ)^2 / (4b + 2var))
    )
    p3 = sum(
        p[k] * p[h]
        * sqrt(2π / (2b + var[k] + var[h]))
        * exp(-(μ[k] - μ[h])^2 / (4b + 2var[k] + 2var[h]))
        for k ∈ 1:K, h ∈ 1:K
    )

    sqrt(π/b) + p1 - 2p2 + 2 * (p3 - p1)
end

"E(X^2)"
function var_err(p::AV, μ::AV, ω::AV, α::AV, β::AV)::Real
    @assert length(p) == 2
    @assert length(μ) == length(p)
    @assert length(ω) == length(p)
    @assert length(α) == length(p)
    @assert length(β) == length(p)

    idx = eachindex(p)

    s1 = sum(p[k] * μ[k]^2 for k in idx)
    @assert s1 >= 0
    s2 = sum(p[k] * ω[k] / (1 - β[k]) for k in idx)
    @assert s2 >= 0

    # Exactly the same as `constr[1]`
    # in `constr!(constr, θ::AV)` below
    s3 = 1 - sum(p[k] * α[k] / (1 - β[k]) for k in idx)

    @show ForwardDiff.value(s3)

    var = (s1 + s2) / s3
    if !(var ≥ 0)
        @assert s3 <= 0
        @show ForwardDiff.value(var)
        @show ForwardDiff.value.(p) ForwardDiff.value.(μ) ForwardDiff.value.(ω) ForwardDiff.value.(α) ForwardDiff.value.(β)
        @assert var ≥ 0
    end

    var
end

"E(σ_k^2)"
function expected_var(p::AV, μ::AV, ω::AV, α::AV, β::AV)
    verr = var_err(p, μ, ω, α, β)

    @. (ω + α * verr) / (1 - β)
end

function D(p::AV, μ::AV, ω::AV, α::AV, β::AV, b::Real, X::AV{<:Real})
    var = expected_var(p, μ, ω, α, β)
    @assert size(var) == size(p)

    ret = D_t(p, μ, var, b, X[1])
    for t ∈ 2:length(X)
        @. var = ω + α * X[t-1]^2 + β * var
        ret += D_t(p, μ, var, b, X[t])
    end

    ret
end

"Extract parameters from vector θ"
function get_params(θ::AV)
    @assert length(θ) % 5 == 0
    K = length(θ) ÷ 5 # 5 variables
    @assert K == 2

    p = @view θ[1:K]
    μ = @view θ[K+1:2K]
    ω = @view θ[2K+1:3K]
    α = @view θ[3K+1:4K]
    β = @view θ[4K+1:end]

    @assert length(p) == K
    @assert length(μ) == K
    @assert length(ω) == K
    @assert length(α) == K
    @assert length(β) == K

    @assert all(p .>= 0)
    @assert all(ω .>= 0)
    @assert all(α .>= 0)
    @assert all(β .>= 0)

    (p, μ, ω, α, β)
end

# ===== Constraints =====
function constr!(constraints, θ::AV)
    p, μ, ω, α, β = get_params(θ)
    idx = eachindex(p)

    @assert length(constraints) == 2

    constraints[2] = sum(p) - 1 # == 0
    constraints[1] = 1 - sum(p[k] * α[k] / (1 - β[k]) for k in idx) # >= 0

    constraints
end

constr(θ::AV) = constr!(zeros(eltype(θ), 2), θ)

"Final objective function for Optim"
D(θ::AV, b::Real, X::AV{<:Real}) = D(get_params(θ)..., b, X)

# ===== Try to optimize D =====
@info "Generating random data..."
distr = UnivariateGMM([-1, 1], [0.2, .3], Categorical([.4, .6]))
data = rand(distr, 100)

θ0 = [
    1/2; 1/2; #p
    0; 0; # μ
    rand(2); # ω
    1e-5; 1e-5; # α
    1e-5; 1e-5
]
@info "Initial value" θ0
@show get_params(θ0)

@show constr(θ0)

@info "Setting up optimization problem..."
objective(θ) = D(θ, 1.0, data)

θ_lo = [
    0, 0, # p
    -Inf, -Inf, # μ
    0, 0, # ω
    0, 0, 0, 0 # α, β
]
θ_hi = [
    1, 1,
    Inf, Inf,
    Inf, Inf,
    Inf, Inf, 1, 1 # α, β
]
# Must lie EXACTLY in interior
@assert all(θ0 .> θ_lo)
@assert all(θ0 .< θ_hi)

# 0 <= constr1 < Inf
# constr2 == 0
constr_lo = Float64[0, 0]
constr_hi = Float64[Inf, 0]
@assert all(constr(θ0) .≥ constr_lo)
@assert all(constr(θ0) .≤ constr_hi)

dfc = TwiceDifferentiableConstraints(
    constr!,
    θ_lo, θ_hi, constr_lo, constr_hi,
    :forward # use autodiff!
)

@info "Optimizing..."
ret = optimize(objective, dfc, θ0, IPNewton(), autodiff=:forward)

display(ret)

```

* * *

What am I doing wrong here??

---

<div class="post-metadata">

### Author: ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [February 4, 2022, 4:30pm UTC](https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573/3 "2022-02-04T16:30:29Z")

</div>

What do you mean when you say the “optimization fails”?

And have you tried another package like JuMP or NLopt?

---

<div class="post-metadata">

### Author: ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)
#### Post date: [February 4, 2022, 4:49pm UTC](https://discourse.julialang.org/t/optim-jl-violates-inequality-constraint/75573/4 "2022-02-04T16:49:11Z")

</div>

By “optimization fails” I mean that Optim.jl violates the inequality constraint, as stated in the title.

When I run this code, the output looks like this:

```julia
      Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_8iC0Np/Project.toml`
  [31c24e10] Distributions v0.25.46
  [f6369f11] ForwardDiff v0.10.25
  [429524aa] Optim v1.6.1
[ Info: Generating random data...
┌ Info: Initial value
│ θ0 =
│ 10-element Vector{Float64}:
│ 0.5
│ 0.5
│ 0.0
│ 0.0
│ 0.3076608766072988
│ 0.2803873329461656
│ 1.0e-5
│ 1.0e-5
│ 1.0e-5
└ 1.0e-5
get_params(θ0) = ([0.5, 0.5], [0.0, 0.0], [0.3076608766072988, 0.2803873329461656], [1.0e-5, 1.0e-5], [1.0e-5, 1.0e-5])
constr(θ0) = [0.999989999899999, 0.0]
[ Info: Setting up optimization problem...
constr! = constr!
[ Info: Optimizing...
ForwardDiff.value(s3) = 0.999989999899999
ForwardDiff.value(s3) = Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.999989999899999,-1.000010000100001e-5,-1.000010000100001e-5,-0.0,-0.0,-0.0,-0.0,-0.5000050000500005,-0.5000050000500005,-5.00010000150002e-6,-5.00010000150002e-6)
ForwardDiff.value(s3) = 0.9999809154864348
...
ForwardDiff.value(s3) = 0.9965220053647488
ForwardDiff.value(s3) = Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.9965220053647488,-0.0035224798058174664,-0.003433432142628737,-0.0,-0.0,-0.0,-0.0,-0.5018859397772173,-0.5010145138215388,-0.001767883087688969,-0.0017201993356783809)
ForwardDiff.value(s3) = -0.06428220890828795
ForwardDiff.value(var) = -1.5104538785258217
ForwardDiff.value.(p) = [0.6340530871847483, 0.3659469128152516]
ForwardDiff.value.(μ) = [0.29075393704923885, 0.09769969796228732]
ForwardDiff.value.(ω) = [0.03220713502201594, 0.0002772290095625074]
ForwardDiff.value.(α) = [0.5171137635535503, 0.581144621641674]
ForwardDiff.value.(β) = [0.4868815931762428, 0.49994821923123]
ERROR: LoadError: AssertionError: var ≥ 0

```

So, `ForwardDiff.value(s3)` is positive at first, as it must be (because `s3` is literally the value of the inequality constraint), but then it becomes negative: `ForwardDiff.value(s3) = -0.06428220890828795`. Thus, the model variance becomes negative as well:

```
ForwardDiff.value(var) = -1.5104538785258217

```

This doesn’t make any sense, because that inequality constraint is there specifically to prevent the variance from being negative!

I tried JuMP + Ipopt, and it doesn’t violate the inequality.