Optim.jl violates inequality constraint?

Output.txt

$ julia-1.7 violation.jl
      Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_KPFxWE/Project.toml`
  [31c24e10] Distributions v0.25.45
  [f6369f11] ForwardDiff v0.10.25
  [429524aa] Optim v1.6.0
[ Info: Generating random data...
constr(θ0) = [0.0, -0.999989999899999]
[ Info: Setting up optimization problem...
[ Info: Optimizing...
constr(ForwardDiff.value.(θ)) = [0.0, -0.999989999899999]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), 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)]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9999809154864348]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.9999809154864348,1.8783647734079602e-5,1.938538407109344e-5,0.0,0.0,0.0,0.0,0.5000138795191937,0.5000048002693726,9.392084575039054e-6,9.692785090612151e-6)]
constr(ForwardDiff.value.(θ)) = [0.0, -0.6571770023852826]
constr(ForwardDiff.value.(θ)) = [0.0, -0.8437648213800245]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9251922334621798]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9633653518010368]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9818618990495095]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9909678837284107]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9954859309381525]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9977362951686276]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9988593219630577]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9994202977148056]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9997006513270554]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9998407945857204]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9999108578304923]
constr(ForwardDiff.value.(θ)) = [0.0, -0.999945887357026]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.999945887357026,5.458440262554897e-5,5.364087019852074e-5,0.0,0.0,0.0,0.0,0.5000265306197792,0.5000119711853303,2.729364947080642e-5,2.6821077244058798e-5)]
constr(ForwardDiff.value.(θ)) = [0.0, -0.49888590866184923]
constr(ForwardDiff.value.(θ)) = [0.0, -0.7975662139092156]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9076118770138573]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9557128606433749]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9782829392521926]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9892243314447801]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9946121702573992]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9972857426163824]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9986174870703866]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9992821044420798]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9996141001077974]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9997800197721305]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.9997800197721305,0.0002108412647806276,0.0002291204872091176,0.0,0.0,0.0,0.0,0.500131225511487,0.5000671190035185,0.00010544830014312723,0.00011457562194334594)]
constr(ForwardDiff.value.(θ)) = [0.0, -0.6255602507702114]
constr(ForwardDiff.value.(θ)) = [0.0, -0.8353706705490167]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9222815961709998]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9621097699038457]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9812036025007097]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9905551744607627]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9951832776349353]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9974855490308131]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9986337570226074]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9992071312459745]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9994936361829256]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.9994936361829256,0.0005231302337052975,0.0004895931457580833,0.0,0.0,0.0,0.0,0.5002675683695269,0.5001365873928526,0.0002617050899563316,0.00024486344513037923)]
constr(ForwardDiff.value.(θ)) = [2.220446049250313e-16, -0.6020530279487354]
constr(ForwardDiff.value.(θ)) = [0.0, -0.830062550700005]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9205893609432059]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9613494501839233]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9807324974552256]
constr(ForwardDiff.value.(θ)) = [2.220446049250313e-16, -0.9901889150365905]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9948600080543417]
constr(ForwardDiff.value.(θ)) = [0.0, -0.997181476933656]
constr(ForwardDiff.value.(θ)) = [2.220446049250313e-16, -0.9983387167499669]
constr(ForwardDiff.value.(θ)) = [0.0, -0.998916466076902]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.998916466076902,0.0010780912359542793,0.001088979658870051,0.0,0.0,0.0,0.0,0.500589071665551,0.5003220442275863,0.0005396806909771191,0.0005448405290481236)]
constr(ForwardDiff.value.(θ)) = [2.220446049250313e-16, -0.588207034569716]
constr(ForwardDiff.value.(θ)) = [0.0, -0.82358998422475]
constr(ForwardDiff.value.(θ)) = [2.220446049250313e-16, -0.9172301862417158]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9594202925369444]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9794887943850324]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9892808077031642]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9941179471996311]
constr(ForwardDiff.value.(θ)) = [0.0, -0.9965220053663375]
constr(ForwardDiff.value.(θ)) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}[Dual{ForwardDiff.Tag{typeof(objective), Float64}}(0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0), Dual{ForwardDiff.Tag{typeof(objective), Float64}}(-0.9965220053663375,0.003522479803975012,0.0034334321412939853,0.0,0.0,0.0,0.0,0.5018859398908682,0.5010145139321996,0.0017678830871645996,0.0017201993353895973)]
constr(ForwardDiff.value.(θ)) = [4.440892098500626e-16, 0.06428113096769317]
ForwardDiff.value(var) = -1.510475486007552
ForwardDiff.value.(p) = [0.6340529090454623, 0.36594709095453803]
ForwardDiff.value.(μ) = [0.2907535090341941, 0.09769954139102985]
ForwardDiff.value.(ω) = [0.03220710544918265, 0.0002772290095625074]
ForwardDiff.value.(α) = [0.5171133614862504, 0.5811441322125269]
ForwardDiff.value.(β) = [0.48688147475346827, 0.4999480981633673]
ERROR: LoadError: AssertionError: var ≥ 0
Stacktrace:
  [1] var_err(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}}, μ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}}, ω::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}}, α::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}}, β::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(objective), Float64}, Float64, 10}})
    @ Main ~/violation.jl:47

violation.jl

import Pkg
Pkg.activate(temp=true, io=devnull)
Pkg.add(["Optim", "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)

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

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

"E(X^2)"
function var_err(p::AV, μ::AV, ω::AV, α::AV, β::AV)
    s1 = sum(@. p * μ^2)
    s2 = sum(@. p * ω / (1 - β))
    s3 = 1 - sum(@. p * α / (1 - β))

    var = (s1 + s2) / s3
    if !(var ≥ 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, μ, ω, α, β)

    ret = @. (ω + α * verr) / (1 - β)
    @assert all(ret .≥ 0)
    ret
end

function D(p::AV, μ::AV, ω::AV, α::AV, β::AV, b::Real, X::AV{<:Real})
    var = expected_var(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 paremeters from vector θ"
function get_params(θ::AV)
    @assert length(θ) % 5 == 0
    K = length(θ) ÷ 5 # 5 variables

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

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

# ===== Constraints =====
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(θ)]

"Final objective function for Optim"
function D(θ::AV, b::Real, X::AV{<:Real})
    @show constr(ForwardDiff.value.(θ))
    D(get_params(θ)..., b, X)
end

# ===== Constraints' Jac and Hess =====
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

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

K = 2
θ0 = [
    ones(K) ./ K; #p
    zeros(K); # μ
    rand(K); # ω
    zeros(K) .+ 1e-5; # α
    zeros(K) .+ 1e-5
]

@show constr(θ0)

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

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

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

dfc = TwiceDifferentiableConstraints(
    constr!, constr_jac!, constr_hess!,
    θ_lo, θ_hi, constr_lo, constr_hi
)

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

display(ret)