Edit on GitHub

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

\[ x^{6} - x^{4}y^{2} - x^{2}y^{4} + y^{6} - x^{4} + 3x^{2}y^{2} - y^{4} - x^{2} - y^{2} + 1 \]

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) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1
SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1

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).sub_gram_matrices
        println(g.basis.polynomials)
    end
end
solve(G)
Iter: 11 Ap: 9.58e-01 Pobj: -1.6713214e-10 Ad: 9.58e-01 Dobj: -4.5283865e-09
Success: SDP solved
Primal objective value: -1.6713214e-10
Dual objective value: -4.5283865e-09
Relative primal infeasibility: 1.16e-14
Relative dual infeasibility: 1.17e-09
Real Relative Gap: -4.36e-09
XZ Relative Gap: 1.21e-09
DIMACS error measures: 2.00e-14 0.00e+00 2.80e-09 0.00e+00 -4.36e-09 1.21e-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.9540693e+00 Ad: 8.42e-01 Dobj:  3.0671681e-01
Iter:  2 Ap: 8.58e-01 Pobj: -1.6732387e+01 Ad: 8.59e-01 Dobj:  7.8448086e+00
Iter:  3 Ap: 8.56e-01 Pobj: -6.5657892e+00 Ad: 7.36e-01 Dobj:  4.4153053e+00
Iter:  4 Ap: 7.93e-01 Pobj: -3.0377721e+00 Ad: 8.60e-01 Dobj:  8.3888624e-01
Iter:  5 Ap: 8.47e-01 Pobj: -1.8654993e+00 Ad: 8.10e-01 Dobj: -2.9577322e-01
Iter:  6 Ap: 7.39e-01 Pobj: -1.2273103e+00 Ad: 8.05e-01 Dobj: -6.7817689e-01
Iter:  7 Ap: 8.93e-01 Pobj: -1.0495527e+00 Ad: 6.76e-01 Dobj: -8.1824309e-01
Iter:  8 Ap: 6.81e-01 Pobj: -9.8438100e-01 Ad: 7.66e-01 Dobj: -8.8707837e-01
Iter:  9 Ap: 8.56e-01 Pobj: -9.5522629e-01 Ad: 7.27e-01 Dobj: -9.1564040e-01
Iter: 10 Ap: 1.00e+00 Pobj: -9.2940978e-01 Ad: 8.95e-01 Dobj: -9.2673695e-01
Iter: 11 Ap: 8.59e-01 Pobj: -9.2935253e-01 Ad: 9.06e-01 Dobj: -9.3052848e-01
Iter: 12 Ap: 1.00e+00 Pobj: -9.2917505e-01 Ad: 1.00e+00 Dobj: -9.3127516e-01
Iter: 13 Ap: 6.96e-01 Pobj: -9.2976072e-01 Ad: 5.87e-01 Dobj: -9.3164776e-01
Iter: 14 Ap: 1.00e+00 Pobj: -9.3001270e-01 Ad: 3.85e-01 Dobj: -9.3189163e-01
Iter: 15 Ap: 3.46e-01 Pobj: -9.3071764e-01 Ad: 1.00e+00 Dobj: -9.3190533e-01
Iter: 16 Ap: 1.00e+00 Pobj: -9.3058083e-01 Ad: 9.82e-01 Dobj: -9.3222213e-01
Iter: 17 Ap: 3.75e-01 Pobj: -9.3121390e-01 Ad: 1.00e+00 Dobj: -9.3246428e-01
Iter: 18 Ap: 4.04e-01 Pobj: -9.3131190e-01 Ad: 5.84e-01 Dobj: -9.3262908e-01
Iter: 19 Ap: 4.39e-01 Pobj: -9.3163773e-01 Ad: 2.42e-01 Dobj: -9.3271332e-01
Iter: 20 Ap: 3.16e-01 Pobj: -9.3173799e-01 Ad: 4.08e-01 Dobj: -9.3277433e-01
Iter: 21 Ap: 6.55e-01 Pobj: -9.3179189e-01 Ad: 6.39e-01 Dobj: -9.3280691e-01
Iter: 22 Ap: 1.00e+00 Pobj: -9.3180218e-01 Ad: 1.00e+00 Dobj: -9.3281932e-01
Iter: 23 Ap: 7.71e-01 Pobj: -9.3182701e-01 Ad: 4.94e-01 Dobj: -9.3282851e-01
Iter: 24 Ap: 1.00e+00 Pobj: -9.3182027e-01 Ad: 9.25e-01 Dobj: -9.3282972e-01
value(t) = -0.9318227936087737
DynamicPolynomials.Polynomial{true, Float64}[xy², x, x³]
DynamicPolynomials.Polynomial{true, Float64}[x²y, y, y³]
DynamicPolynomials.Polynomial{true, Float64}[1.0, -0.7071067811865472x² - 0.7071067811865475y²]
DynamicPolynomials.Polynomial{true, Float64}[xy]
DynamicPolynomials.Polynomial{true, Float64}[-0.7071067811865472x² + 0.7071067811865475y²]

We notice that we indeed find -3825/4096 and that symmetry was exploited.

This page was generated using Literate.jl.