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 Clarabel
function solve(G)
solver = Clarabel.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)
gram_matrix(con_ref)
end
solve(G)BlockDiagonalGramMatrix with 4 blocks:
[1] Block with row/column basis:
Simple basis:
FixedBasis([1.0·1 + 0.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³, 0.0·1 + 0.0·y + 0.0·x - 0.7071067811865472·y² + 0.0·xy - 0.7071067811865475·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³])
And entries in a 2×2 SymMatrix{Float64}:
1.933837593145518 0.9833204153958562
0.9833204153958562 0.5000001619295955
[2] Block with row/column basis:
Simple basis:
FixedBasis([0.0·1 + 0.0·y + 0.0·x + 0.0·y² + 1.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³])
And entries in a 1×1 SymMatrix{Float64}:
1.8739974090671922e-8
[3] Block with row/column basis:
Simple basis:
FixedBasis([0.0·1 + 0.0·y + 0.0·x - 0.7071067811865472·y² + 0.0·xy + 0.7071067811865475·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³])
And entries in a 1×1 SymMatrix{Float64}:
3.230737160780864e-8
[4] Block with row/column basis:
Semisimple basis with 2 simple sub-bases:
FixedBasis([0.0·1 + 0.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 1.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³, 0.0·1 + 0.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 1.0·x²y + 0.0·x³, 0.0·1 + 1.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³])
FixedBasis([0.0·1 + 0.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y - 1.0·x³, 0.0·1 + 0.0·y + 0.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ - 1.0·xy² + 0.0·x²y + 0.0·x³, 0.0·1 + 0.0·y - 1.0·x + 0.0·y² + 0.0·xy + 0.0·x² + 0.0·y³ + 0.0·xy² + 0.0·x²y + 0.0·x³])
And entries in a 3×3 SymMatrix{Float64}:
1.0000000008594137 -0.9999998695393238 -0.6250000485614923
-0.9999998695393238 0.9999997407910516 0.62499996290667
-0.6250000485614923 0.62499996290667 0.390625067611156We notice that we indeed find -3825/4096 and that symmetry was exploited.
This page was generated using Literate.jl.