Contributed by: Benoît Legat
Given a Sum-of-Squares constraint on an algebraic set: $$ g_1(x) = , \ldots, g_m(x) = 0 \Rightarrow p(x) \ge 0. $$ We can either use the certificate: $$ p(x) = s(x) + \lambda_1(x) g_1(x) + \cdots + \lambda_m(x) g_m(x), s_0(x) \text{ is SOS}, $$ or $$ p(x) \equiv s(x) \pmod{\la g_1(x), \ldots, g_m(x) \ra}, s_0(x) \text{ is SOS}. $$ the second one leads to a simpler SDP but needs to compute a Gr"obner basis:
@set
macro recognizes variable fix, e.g., x = 1
and provides shortcut.We illustrate this in this example.
using DynamicPolynomials
@polyvar x[1:3]
p = sum(x)^2
using SumOfSquares
S = algebraicset([xi^2 - 1 for xi in x])
Algebraic Set defined by 3 equalities -1//1 + x[1]^2 = 0 -1//1 + x[2]^2 = 0 -1//1 + x[3]^2 = 0
We will now search for the minimum of x
over S
using Sum of Squares Programming.
We first need to pick an SDP solver, see here for a list of the available choices.
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
function min_algebraic(S)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = S)
optimize!(model)
@show termination_status(model)
@show objective_value(model)
end
min_algebraic(S)
termination_status(model) = MathOptInterface.OPTIMAL objective_value(model) = -4.404406006575101e-10
-4.404406006575101e-10
Note that the minimum is in fact 1
.
Indeed, since each variables is odd (it is either -1
or 1
)
and there is an odd number of variables, their sum is odd.
Therefore it cannot be zero!
We can see that the Gröbner basis of S
was computed
@show S.I.gröbner_basis
S.I.algo
S.I.gröbner_basis = false
Buchberger(1.4901161193847656e-8, SemialgebraicSets.presort!, SemialgebraicSets.normal_selection)
The Gröbner basis is simple to compute in this case as the vector
of xi^2 - 1
is already a Gröbner basis.
However, we still need to divide polynomials by the Gröbner basis
which can be simplified in this case.
const MP = MultivariatePolynomials
const SS = SemialgebraicSets
struct HypercubeIdeal{V} <: SS.AbstractPolynomialIdeal
variables::Vector{V}
end
struct HypercubeSet{V} <: SS.AbstractAlgebraicSet
ideal::HypercubeIdeal{V}
end
MP.variables(set::HypercubeSet) = MP.variables(set.ideal)
MP.variables(ideal::HypercubeIdeal) = ideal.variables
Base.similar(set::HypercubeSet, ::Type) = set
SS.ideal(set::HypercubeSet) = set.ideal
function Base.rem(p, set::HypercubeIdeal)
return MP.polynomial(map(MP.terms(p)) do term
mono = MP.monomial(term)
new_mono = one(mono)
for (var, exp) in powers(mono)
if var in set.variables
exp = rem(exp, 2)
end
new_mono *= var^exp
end
MP.coefficient(term) * new_mono
end)
end
H = HypercubeSet(HypercubeIdeal(x))
min_algebraic(H)
termination_status(model) = MathOptInterface.OPTIMAL objective_value(model) = -4.404406006575101e-10
-4.404406006575101e-10
Let's now try to find the correct lower bound:
function min_algebraic_rational(S, d)
model = SOSModel(solver)
@variable(model, q, SOSPoly(MP.monomials(x, 0:d)))
deno = q + 1
@constraint(model, c, deno * p >= deno, domain = S)
optimize!(model)
@show termination_status(model)
end
min_algebraic_rational (generic function with 1 method)
With d = 0
, it's the same as previously
min_algebraic_rational(H, 0)
termination_status(model) = MathOptInterface.INFEASIBLE
INFEASIBLE::TerminationStatusCode = 2
But with d = 1
, we can find the correct lower bound
min_algebraic_rational(H, 1)
termination_status(model) = MathOptInterface.OPTIMAL
OPTIMAL::TerminationStatusCode = 1
This notebook was generated using Literate.jl.