JuMP interface for Optimization using Homotopy Continuation

[freemin7]

Homotopy Continuation can be used to solve and optimize highly non convex problems. JuMP at the moment has really limited support for polynomial systems (that my change by JuMP 2.0) however there is an JuMP extension https://github.com/jump-dev/PolyJuMP.jl that supports polynomial systems. Would there be interest that at some future point in time your work is available from the JuMP ecosystem?

Imagine that Homotopy Continuation might also be useful to solve really helpful for very non convex problems by finding minimas of a higher order (3rd/4th) polynomial approximation inside some trust region and do so for different trust regions.

[PBrdng]

Would there be interest that at some future point in time your work is available from the JuMP ecosystem?

Absolutely, yes. It is certainly a good idea to integrate both systems. As a side note: I am also quite interested reasearch-wise to apply these algebraic methods to optimization. So yes, let's do it!

[freemin7]

I would translate a small problem from JuMP in to a polynomial form in your language, so we can have a first idea what is involved in the translation. Do you see a trivial way to implement box constraints like r \in (0, 5)? The only way i can imagine now is: r + s1^2 == 5 r - s2^2 == 0 Where s1, s2 are \in \mathbb{R}. This leads to many additional terms and variables that occur only once.

[PBrdng]

I would not put inequality constraints, but simply check a posteriori, which variables satisfy the constraints.

You should also note that homotopy continuation computes in C, so the variable explosion also causes an explosion in the number of solutions.

[freemin7]

The problem i am modeling is: place 4 circles of arbitrary size on a 2 by 3 unit long in an non overlapping way to maximize the amount of area covered. A JuMP implementation would look like

using JuMP, Ipopt

a = 3
b = 2

xmax2 = [a,a,a,a,b,b,b,b,min(a,b)/2,min(a,b)/2,min(a,b)/2,min(a,b)/2]
xmin2 = eps() * ones(12)

model = Model()

x = @variable(model, x[1:12])

for i in 1:12 set_lower_bound(x[i],xmin2[i]) end
for i in 1:12 set_upper_bound(x[i],xmax2[i]) end

@objective(model, Min, 3*2-π*sum(x[9:12].^2))

@constraint(model,0 .<= x[1:4] .- 2*x[9:12] )
@constraint(model,0 .<= x[5:8] .- 2*x[9:12] )

for i in 2:4 
    for j in 1:(i-1)
            @constraint(model,
0 <= ((x[i]-x[i+8])-(x[j]-x[j+8]))^2 + ((x[i+4]-x[i+8])-(x[j+4]-x[j+8]))^2 - (x[j+8]+x[i+8])^2 )
    end
 end

set_optimizer(model,Ipopt.Optimizer)
optimize!(model)

I am not sure this optimization problem is sensible when we remove the inequality constraints. I am aware that this problem is not "nice" for a Homotopy approach but that's why i choose it.

So i transcribed the problem into your language but i run into a problem with the inequality constrains i think

using HomotopyContinuation

@var r[1:4] x[1:4] y[1:4]
@var s[1:22] λ[1:22]

si = 1;
acc = [];

for i in 1:4
    push!(acc, y[i] + r[i] - 2 + s[si]^2)
    si+=1;
end
for i in 1:4
    push!(acc, x[i] + r[i] - 3 + s[si]^2)
    si+=1;
end
for i in 1:4
    push!(acc, y[i] - r[i] - s[si]^2)
    si+=1;
end
for i in 1:4
    push!(acc, x[i] - r[i] - s[si]^2)
    si+=1;
end
for i in 2:4
    for j in 1:(i-1)
        push!(acc, s[si]^2 + (x[i]-x[j])^2 + (y[i]-y[j])^2 - (r[i]+r[j])^2)
        si += 1;
    end
end

J = (ac -> differentiate(ac, [x;y;r])).(acc)

C = System([- r[1]^2 - r[2]^2 - r[3]^2 - r[4]^2 - sum(J'*λ); acc...], variables=[x;y;r;λ;s],)

res = solve(C)

ERROR: FiniteException: The solution set of the given system has at least dimension 33 > 0. Consider intersecting with an (affine) subspace of codimension 33 to reduce to (possibly) finitely many solutions.

[freemin7]

This "toy problem" has no equality constraints and the unconstrained problem would just blow up in a similar way:

C = System([- r[1]^2 - r[2]^2 - r[3]^2 - r[4]^2], variables=r)

julia> res = solve(C)
ERROR: FiniteException: The solution set of the given system has at least dimension 2 > 0. Consider intersecting with an (affine) subspace of codimension 2 to reduce to (possibly) finitely many solutions.

Any ideas where to go from here?

[PBrdng]

Hi,

for such problems that have inequalities one can use the KTT conditions.

Below, you can find code that writes down those conditions for your toy problem. This system, however, has a lot of structure. Observe that ℓ .* constraints consists of products of variables. Elimination simplifies the problem.

However, in general, for this simple problem, there are 137817 complex solutions. I imagine that more complex problems are nearly infeasible. This is why modelling a problem as zeros of polynomial systems often is not suited to problems involving inequalities.

(I deleted an earlier answer. The reference to KKT makes my answer clearer)

using HomotopyContinuation, LinearAlgebra

@var x[1:4] # x-coordinates of 4 circles
@var y[1:4] # y-coordinates of 4 circles

r = [x[i]^2 + y[i]^2 for i in 1:4]

# x < 2
constraints_x = [2 .- x; x]
# y < 3
constraints_y = [3 .- y; y]
# non-intersection
contraints_r =
[
(x[i] - x[j])^2 + (y[i] - y[j])^2 - r[i] - r[j] for i in 1:4 for j in i+1:4
]

objective = sum(r)
constraints = [constraints_x; constraints_y; contraints_r]
k = length(constraints)
# KKT
@var ℓ[1:k]
∇₁ = differentiate(objective, [x; y])
∇₂ = differentiate(constraints, [x; y])

F = ∇₁ + sum(ℓ[i] .* ∇₂[i,:] for i in 1:k)

# System that we need to solve
f = System([F; ℓ .* constraints])

[blegat]

PolyJuMP models problem where the decision variables are part of the coefficients of the polynomials. We are not trying to find values of the variables of the polynomials, we are rather trying to find values of the decision variables so that the polynomial inequalities and equalities are satisfied for all values of the polynomial variables. In https://github.com/tweisser/PolyModels and https://github.com/lanl-ansi/MomentOpt.jl (cc @tweisser), the aim is to find values of the polynomial variables instead.

@PBrdng Solving optimization problems using solvers for systems of algebraic equations might seem crazy as there may be too many solution but I have also seen it here: https://www.sciencedirect.com/science/article/pii/S1474667015381179 where they present another approach for solving polynomial systems (it's basically Groebner Basis but instead of using the numericaly unstable Gaussian Elimination to in the Groebner basis computation they use SVD).

I would also be interested to write a MOI solver parametrized by a SemialgebraicSets.AbstractAlgebraicSolver such as the ones defined in https://github.com/JuliaAlgebra/SemialgebraicSets.jl/blob/master/src/solve.jl and HomotopyContinuation (since https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/pull/445). One use case is to find algebraic expressions for the solution of optimization problems to be used for writing tests. When I write a tests and I don't know the solution of a problem, it's a bit annoying to have to write @test ... ≈ -0.6180339887498949 but if a solver find that the solution is (1 - √5) / 2 then you can write @test ... ≈ (1 - √5) / 2 :)

[PBrdng]

Solving optimization problems using solvers for systems of algebraic equations might seem crazy as there may be too many solution

It's not crazy for some optimization problems. A rule of thumb is that for few inequalities, solving the KKT equations can be done. Note that HC.jl can solve systems of equations with up to several-ten-thousand solutions easily.

I would also be interested to write a MOI solver parametrized by a SemialgebraicSets.AbstractAlgebraicSolver such as the ones defined in https://github.com/JuliaAlgebra/SemialgebraicSets.jl/blob/master/src/solve.jl and HomotopyContinuation (since #445). One use case is to find algebraic expressions for the solution of optimization problems to be used for writing tests. When I write a tests and I don't know the solution of a problem, it's a bit annoying to have to write @test ... ≈ -0.6180339887498949 but if a solver find that the solution is (1 - √5) / 2 then you can write @test ... ≈ (1 - √5) / 2 :)

You mean something like @test f(x)=0 has a solution with property A?

[blegat]

You mean something like @test f(x)=0 has a solution with property A?

You mean something like @test f(x)=0 has a solution with property A?

Sometimes, you get the closed form expression so you can hardcode the algebraic expression of the solution, see e.g.: https://github.com/jump-dev/MathOptInterface.jl/blob/8d8beb387551d45138857bd472328718ba98ca8e/src/Test/contconic.jl#L4074-L4131. The solution could also be the root of a polynomial of high degree so you don't get a closed form expression but we try to have simple tests so this does not happen too often.

[PBrdng]

I understand, but is it a polynomial (univariate) or a system of polynomial equations?