# Hypercube

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:

• SemialgebraicSets implements Buchberger's algorithm.
• The @set macro recognizes variable fix, e.g., x = 1 and provides shortcut.
• If you know a \alert{better} way to take modulo, better create your \alert{own} type of algebraic set!

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)
-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
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)
-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)
INFEASIBLE::TerminationStatusCode = 2

But with d = 1, we can find the correct lower bound

min_algebraic_rational(H, 1)
OPTIMAL::TerminationStatusCode = 1