Interior-Point Newton tries points outside of the constraint set

New to Julia
, ,
1

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.

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

"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?

2

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

1 Like
3

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).

1 Like
4

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).

1 Like
6

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

7

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

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

# 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