# Dihedral symmetry of the Robinson form

Adapted from: Example 5.4 of [GP04]

[GP04] Gatermann, Karin and Parrilo, Pablo A. Symmetry groups, semidefinite programs, and sums of squares. Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128.

We start by defining the Dihedral group of order 8. This group is isomorphic to the following permutation group:

using PermutationGroups
d = perm"(1, 2, 3, 4)"
c = perm"(1, 3)"
G = PermGroup([c, d])
Permutation group on 2 generators generated by
(1,3)
(1,2,3,4)

We could rely on this isomorphism to define this group. However, in order to illustrate how to do symmetry reduction with a custom group, we show in this example what should be implemented to define a new group.

import GroupsCore

struct DihedralGroup <: GroupsCore.Group
n::Int
end

struct DihedralElement <: GroupsCore.GroupElement
n::Int
reflection::Bool
id::Int
end

Implementing GroupsCore API:

Base.one(G::DihedralGroup) = DihedralElement(G.n, false, 0)

Base.eltype(::Type{DihedralGroup}) = DihedralElement
Base.IteratorSize(::Type{DihedralGroup}) = Base.HasLength()

function Base.iterate(G::DihedralGroup, prev::DihedralElement=DihedralElement(G.n, false, -1))
if prev.id + 1 >= G.n
if prev.reflection
return nothing
else
next = DihedralElement(G.n, true, 0)
end
else
next = DihedralElement(G.n, prev.reflection, prev.id + 1)
end
return next, next
end

GroupsCore.order(::Type{T}, G::DihedralGroup) where {T} = convert(T, 2G.n)
GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, false, 1), DihedralElement(G.n, true, 0)]

Base.rand not needed for our purposes here

Base.parent(g::DihedralElement) = DihedralGroup(g.n)
function Base.:(==)(g::DihedralElement, h::DihedralElement)
return g.n == h.n && g.reflection == h.reflection && g.id == h.id
end

function Base.inv(el::DihedralElement)
if el.reflection || iszero(el.id)
return el
else
return DihedralElement(el.n, false, el.n - el.id)
end
end
function Base.:*(a::DihedralElement, b::DihedralElement)
a.n == b.n || error("Cannot multiply elements from different Dihedral groups")
id = mod(a.reflection ? a.id - b.id : a.id + b.id, a.n)
return DihedralElement(a.n, a.reflection != b.reflection, id)
end

Base.copy(a::DihedralElement) = DihedralElement(a.n, a.reflection, a.id)

optional functions:

function GroupsCore.order(el::DihedralElement)
if el.reflection
return 2
else
if iszero(el.id)
return 1
else
return div(el.n, gcd(el.n, el.id))
end
end
end

The Robinson form is invariant under the following action of the Dihedral group on monomials: The action of each element of the groups is to map the variables x, y to:

idrotationreflection
0x, yy, x
1-y, x-x, y
2-x, -y-y, -x
3y, -xx, -y
using SumOfSquares
using DynamicPolynomials
@polyvar x y
struct DihedralAction <: Symmetry.OnMonomials end
import SymbolicWedderburn
SymbolicWedderburn.coeff_type(::DihedralAction) = Float64
function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial)
if iseven(el.reflection + el.id)
var_x, var_y = x, y
else
var_x, var_y = y, x
end
sign_x = 1 <= el.id <= 2 ? -1 : 1
sign_y = 2 <= el.id ? -1 : 1
return mono([x, y] => [sign_x * var_x, sign_y * var_y])
end

poly = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1
$$$1 - y^{2} - x^{2} - y^{4} + 3x^{2}y^{2} - x^{4} + y^{6} - x^{2}y^{4} - x^{4}y^{2} + x^{6}$$$

We can verify that poly is indeed invariant under the action of each element of the group as follows.

G = DihedralGroup(4)
for g in G
@show SymbolicWedderburn.action(DihedralAction(), g, poly)
end
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶
SymbolicWedderburn.action(DihedralAction(), g, poly) = 1 - y² - x² - y⁴ + 3x²y² - x⁴ + y⁶ - x²y⁴ - x⁴y² + x⁶

We can exploit this symmetry for reducing the problem using the SymmetricIdeal certificate as follows:

import CSDP
function solve(G)
solver = CSDP.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, DihedralAction())
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
optimize!(model)
@show value(t)

for g in gram_matrix(con_ref).blocks
println(g.basis.polynomials)
end
end
solve(G)
Iter: 11 Ap: 9.58e-01 Pobj: -1.6774514e-10 Ad: 9.58e-01 Dobj: -4.5554076e-09
Success: SDP solved
Primal objective value: -1.6774514e-10
Dual objective value: -4.5554076e-09
Relative primal infeasibility: 1.10e-14
Relative dual infeasibility: 1.18e-09
Real Relative Gap: -4.39e-09
XZ Relative Gap: 1.22e-09
DIMACS error measures: 1.89e-14 0.00e+00 2.82e-09 0.00e+00 -4.39e-09 1.22e-09
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 7.31e-01 Pobj: -5.9540730e+00 Ad: 8.42e-01 Dobj:  3.0671681e-01
Iter:  2 Ap: 8.67e-01 Pobj: -1.6845198e+01 Ad: 8.68e-01 Dobj:  7.9208056e+00
Iter:  3 Ap: 8.41e-01 Pobj: -6.5167134e+00 Ad: 6.91e-01 Dobj:  4.6340995e+00
Iter:  4 Ap: 7.89e-01 Pobj: -3.0718725e+00 Ad: 8.67e-01 Dobj:  8.7844427e-01
Iter:  5 Ap: 8.59e-01 Pobj: -1.8420731e+00 Ad: 8.47e-01 Dobj: -3.3091073e-01
Iter:  6 Ap: 7.40e-01 Pobj: -1.2180543e+00 Ad: 8.00e-01 Dobj: -6.8726801e-01
Iter:  7 Ap: 8.82e-01 Pobj: -1.0523759e+00 Ad: 6.48e-01 Dobj: -8.1732010e-01
Iter:  8 Ap: 6.75e-01 Pobj: -9.8406190e-01 Ad: 7.70e-01 Dobj: -8.8525802e-01
Iter:  9 Ap: 8.70e-01 Pobj: -9.5552035e-01 Ad: 7.26e-01 Dobj: -9.1517421e-01
Iter: 10 Ap: 1.00e+00 Pobj: -9.2933412e-01 Ad: 8.92e-01 Dobj: -9.2629286e-01
Iter: 11 Ap: 8.61e-01 Pobj: -9.2914325e-01 Ad: 9.11e-01 Dobj: -9.3036068e-01
Iter: 12 Ap: 1.00e+00 Pobj: -9.2898054e-01 Ad: 1.00e+00 Dobj: -9.3115787e-01
Iter: 13 Ap: 7.19e-01 Pobj: -9.2961473e-01 Ad: 5.76e-01 Dobj: -9.3155210e-01
Iter: 14 Ap: 1.00e+00 Pobj: -9.2967228e-01 Ad: 3.55e-01 Dobj: -9.3180771e-01
Iter: 15 Ap: 4.26e-01 Pobj: -9.3065963e-01 Ad: 1.00e+00 Dobj: -9.3175040e-01
Iter: 16 Ap: 1.00e+00 Pobj: -9.3042297e-01 Ad: 1.00e+00 Dobj: -9.3211446e-01
Iter: 17 Ap: 5.06e-01 Pobj: -9.3106101e-01 Ad: 6.25e-01 Dobj: -9.3240546e-01
Iter: 18 Ap: 7.18e-01 Pobj: -9.3146998e-01 Ad: 4.44e-01 Dobj: -9.3257853e-01
Iter: 19 Ap: 1.00e+00 Pobj: -9.3128112e-01 Ad: 3.48e-01 Dobj: -9.3271014e-01
Iter: 20 Ap: 4.85e-01 Pobj: -9.3206037e-01 Ad: 1.00e+00 Dobj: -9.3257740e-01
Iter: 21 Ap: 1.00e+00 Pobj: -9.3178232e-01 Ad: 1.00e+00 Dobj: -9.3276737e-01
Iter: 22 Ap: 7.26e-01 Pobj: -9.3197870e-01 Ad: 7.39e-01 Dobj: -9.3290369e-01
Iter: 23 Ap: 4.04e-01 Pobj: -9.3208730e-01 Ad: 4.06e-01 Dobj: -9.3295435e-01
Iter: 24 Ap: 1.00e+00 Pobj: -9.3209505e-01 Ad: 1.00e+00 Dobj: -9.3296372e-01
Iter: 25 Ap: 7.03e-01 Pobj: -9.3212572e-01 Ad: 7.08e-01 Dobj: -9.3298098e-01
Iter: 26 Ap: 9.98e-01 Pobj: -9.3212657e-01 Ad: 9.95e-01 Dobj: -9.3298156e-01
value(t) = -0.9321266534577521
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-0.7071067811865472y² + 0.7071067811865475x²]
We notice that we indeed find -3825/4096 and that symmetry was exploited.