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
endImplementing 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
endThe 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:
| id | rotation | reflection |
|---|---|---|
| 0 | x, y | y, x |
| 1 | -y, x | -x, y |
| 2 | -x, -y | -y, -x |
| 3 | y, -x | x, -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)
endSymbolicWedderburn.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}[y³, x²y, y]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-x³, -xy², -x]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[1.0, -0.7071067811865472y² - 0.7071067811865475x²]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[xy]
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.
This page was generated using Literate.jl.