# Interior-Point Newton tries points outside of the constraint set

**URL:** <https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215>\
**Category:** New to Julia\
**Tags:** optim, newton, constraint\
**Created:** [March 21, 2025, 10:12am UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215 "2025-03-21T10:12:23Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [March 21, 2025, 10:12am UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/1 "2025-03-21T10:12:23Z")

</div>

Hi, I’m using Interior-Point Newton to solve a constrained optimization problem where my objective function is not (always) defined outside the constraint set. I found that Interior-Point Newton sometimes tries points outside of the constraint set. Does that mean that I should use a different solution method, or am I doing something wrong? See the below MWE in which I tried to replicate the behaviour.

```julia
using Optim, ForwardDiff, Random, Distributions
Random.seed!(42)

f_unforced(x) = -prod(x)
f_forced(x) = sum(x .^ 2) < 1 ? -prod(x) : error("The value $([xi.value for xi in x]) is not allowed")

norm_c!(c, x) = (c[1] = sum(y -> y^2, x); c)

function check_optims(d)
    con = TwiceDifferentiableConstraints(norm_c!, -ones(d), ones(d), [-Inf], [1.])

    for i in 1:1000
        x_init = rand(d)
        while sum(x_init .^ 2) > 1
            x_init = rand(d)
        end
        x_init_copy = copy(x_init)
        # If we use a function that does not force the constraint, there is no problem
        opt2 = Optim.optimize(f_unforced, con, x_init, IPNewton(); autodiff = :forward)
        try
            # If we use a function that forces the constraint, we may get an error
            opt1 = Optim.optimize(f_forced, con, x_init, IPNewton(); autodiff = :forward)
        catch
            display("Failed at the i'th attempt: $i")
            if sum(x_init_copy .^ 2) > 1
                display("Initial value was $x_init_copy, which is not allowed")
            else
                display("Initial value was $x_init_copy, which is allowed")
            end
            error()
        end
    end
    return nothing
end

check_optims(2)

```

gives

```julia
"Failed at the i'th attempt: 1"
"Initial value was [0.6293451231426089, 0.4503389405961936], which is allowed"
ERROR: 
Stacktrace:
 [1] error()
   @ Base ./error.jl:44
 [2] check_optims(d::Int64)
   @ Main ~/Documents/ConOptStop_Julia_v0.1/test_optim.jl:30
 [3] top-level scope
   @ ~/Documents/ConOptStop_Julia_v0.1/test_optim.jl:36

caused by: The value [0.9466761003329859, 0.45177899946620614] is not allowed
# A long stacktrace

```

If I should use a different solver, I’m curious about your suggestions, so please let me outline my set-up: I typically have box-constraints + a norm constraint in dimensions 1-3, and the objective is (roughly) polynomial in shape (although I do not have a nice closed form available!). It needs to be computed many times (10.000s of times). These repetitions are independent, so I’m currently multithreading the optimizations; I’d like to keep that (or a similar optimization) in there! Of course I can “fix” `f_forced` by setting it to be `f_unforced` if the constraint is met, and infinity otherwise, but I’m not sure if that interferes with how interior-point Newton works?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 21, 2025, 10:31am UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/2 "2025-03-21T10:31:21Z")

</div>

I’m not sure how Optim’s interior-point method works, but you might find this post informative

> [@OptimizationMOI Ipopt violating inequality constraint](https://discourse.julialang.org/t/optimizationmoi-ipopt-violating-inequality-constraint/92608/2):
>
> That is a frequent misconception. If a problem has a constraint g(x) \leq 0, for an interior-point method to ensure satisfaction of the constraint throughout, it would have to penalize g(x) in the logarithmic barrier. However, that’s not what IPOPT does. Firstly, IPOPT adds slack variables s to transform the constraint to g(x) + s = 0 and s \geq 0. Secondly, the inequality on s is penalized via the log barrier, i.e., s \> 0 (strictly) is maintained throughout. BUT g(x) + s = 0 is NOT maintained…

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [March 21, 2025, 10:42am UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/3 "2025-03-21T10:42:44Z")

</div>

Thank you, that explains a lot. Based on that thread and my previous experience with these packages, I think the time-efficient way forward for me is to reduce to a (purely) box-constraint problem through an explicit coordinate transformation for now!

edit: the coordinate transformations of course introduce a non-differentiability; I’ll instead just have to accept that I’ll need to make sure my functions are defined on a slightly bigger set for now (which fortunately isn’t impossible).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2025, 11:19am UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/4 "2025-03-21T11:19:34Z")

</div>

> [@JADekker](#):
>
> edit: the coordinate transformations of course introduce a non-differentiability; I’ll instead just have to accept that I’ll need to make sure my functions are defined on a slightly bigger set for now (which fortunately isn’t impossible).

It sounds like you want to minimize f(x) subject to \Vert x \Vert \le 1, i.e. over the interior of a unit ball. One formulation that uses box constraints without introducing a non-differentiability is

\min\_{r \in \mathbb{R}, x \in \mathbb{R}^n} f(r x / \Vert x \Vert) \\ \text{s.t. } r \le 1

which introduces one additional scalar variable r with a simple box constraint (which many optimization algorithms can strictly enforce). (Even though you’ve “added a dimension”, the optimization shouldn’t really get harder, because you’ve effectively subtracted a dimension from x: the objective function is insensitive to the length of x.)

It’s differentiable except at x=0, but that doesn’t matter — the gradient is \perp x, so if you initialize x to a random point on the unit sphere, the length of x should stay \approx 1 during optimization.

Alternatively, if you know that your constraint is _active_ at the optimum, i.e. that your optimum is _on_ the unit sphere, you can simply optimize f(x / \Vert x \Vert).

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [March 21, 2025, 12:27pm UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/6 "2025-03-21T12:27:44Z")

</div>

> [@stevengj](#):
>
> > [@JADekker](#):
> >
> > edit: the coordinate transformations of course introduce a non-differentiability; I’ll instead just have to accept that I’ll need to make sure my functions are defined on a slightly bigger set for now (which fortunately isn’t impossible).
> 
> It sounds like you want to minimize f(x) subject to \Vert x \Vert \le 1, i.e. over the interior of a unit ball. One formulation that uses box constraints without introducing a non-differentiability is
> 
> \min\_{r \in \mathbb{R}, x \in \mathbb{R}^n} f(r x / \Vert x \Vert) \\ \text{s.t. } r \le 1
> 
> which introduces one additional scalar variable r with a simple box constraint (which many optimization algorithms can strictly enforce). (Even though you’ve “added a dimension”, the optimization shouldn’t really get harder, because you’ve effectively subtracted a dimension from x: the objective function is insensitive to the length of x.)
> 
> It’s differentiable except at x=0, but that doesn’t matter — the gradient is \perp x, so if you initialize x to a random point on the unit sphere, the length of x should stay \approx 1 during optimization.
> 
> Alternatively, if you know that your constraint is _active_ at the optimum, i.e. that your optimum is _on_ the unit sphere, you can simply optimize f(x / \Vert x \Vert).

Thank you, I hadn’t heard about that method before, that seems exactly what I need!

---

<div class="post-metadata">

**Author:** ![JADekker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jadekker/32/210281_2.png) [@JADekker](https://discourse.julialang.org/u/JADekker)\
**Post date:** [March 21, 2025, 2:12pm UTC](https://discourse.julialang.org/t/interior-point-newton-tries-points-outside-of-the-constraint-set/127215/7 "2025-03-21T14:12:59Z")

</div>

I’ve tried a very simple implementation of this, but my optimizer is having some trouble finding the optimum. Am I overlooking something?

```julia
using Optim, ForwardDiff, Random, Distributions, BenchmarkTools
Random.seed!(42)

function sphere_trick(a::Vector{T}) where {T<:Real}
    a_coords = view(a, 1:length(a)-1)
    if all(x -> x == 0., a_coords)
        return a_coords
    else
        return a_coords * a[end] / (sqrt(sum(y -> y^2, a_coords)))
    end
end

f(x) = -prod(x)

norm_c!(c, x) = (c[1] = sum(y -> y^2, x); c)

function check_optims(d)
    con = TwiceDifferentiableConstraints(norm_c!, -ones(d), ones(d), [-Inf], [1.])
    con_simple = TwiceDifferentiableConstraints(-vcat(ones(d), [0.]), vcat(ones(d), [1.]))

    x_init = rand(d)
    while sum(x_init .^ 2) > 1
        x_init = rand(d)
    end
    # If we use a function that does not force the constraint, there is no problem
    opt_con = Optim.optimize(f, con, x_init, IPNewton(); autodiff = :forward)

    # If we use a function that forces the constraint, we may get an error
    opt_r = Optim.optimize(x -> f(sphere_trick(x)), con_simple, vcat(x_init, [0.5]), IPNewton(); autodiff = :forward)

    x_con = opt_con.minimizer
    x_r = opt_r.minimizer
    display(x_con)
    display(x_r)
    display(f(x_con))
    display(f(x_r[1:end-1]))
    return nothing
end

check_optims(2)

```

gives

```julia
# Original optimum
2-element Vector{Float64}:
 0.7071067811865474
 0.7071067811865475

# Found by the new method
3-element Vector{Float64}:
 8.584942598615005e-13
 8.584943133736025e-13
 0.9999999999999599

# Original optimal value
-0.49999999999999983

# New found value (which is suboptimal)
-7.370124401549779e-25

```