Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JuMP interface for Optimization using Homotopy Continuation #446

Open
freemin7 opened this issue Jan 22, 2021 · 10 comments
Open

JuMP interface for Optimization using Homotopy Continuation #446

freemin7 opened this issue Jan 22, 2021 · 10 comments

Comments

@freemin7
Copy link

freemin7 commented Jan 22, 2021

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
Copy link
Collaborator

PBrdng commented Jan 22, 2021

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
Copy link
Author

freemin7 commented Jan 22, 2021

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
Copy link
Collaborator

PBrdng commented Jan 22, 2021

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
Copy link
Author

freemin7 commented Jan 22, 2021

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
Copy link
Author

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
Copy link
Collaborator

PBrdng commented Jan 25, 2021

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
Copy link
Contributor

blegat commented Jan 25, 2021

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 #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
Copy link
Collaborator

PBrdng commented Feb 1, 2021

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
Copy link
Contributor

blegat commented Feb 1, 2021

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
Copy link
Collaborator

PBrdng commented Feb 5, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants