Cyclic symmetry for Sums of Hermitian Squares

We start by defining the Cyclic group.

using GroupsCore
import PermutationGroups

struct CyclicElem <: GroupElement
    n::Int
    id::Int
end
Base.:(==)(a::CyclicElem, b::CyclicElem) = a.n == b.n && a.id == b.id
Base.inv(el::CyclicElem) = CyclicElem(el.n, (el.n - el.id) % el.n)

function Base.:*(a::CyclicElem, b::CyclicElem)
    return CyclicElem(a.n, (a.id + b.id) % a.n)
end
Base.:^(el::CyclicElem, k::Integer) = CyclicElem(el.n, (el.id * k) % el.n)

Base.conj(a::CyclicElem, b::CyclicElem) = inv(b) * a * b
Base.:^(a::CyclicElem, b::CyclicElem) = conj(a, b)

function PermutationGroups.order(el::CyclicElem)
    return div(el.n, gcd(el.n, el.id))
end

struct CyclicGroup <: Group
    n::Int
end
Base.eltype(::CyclicGroup) = CyclicElem
Base.one(c::Union{CyclicGroup, CyclicElem}) = CyclicElem(c.n, 0)
PermutationGroups.gens(c::CyclicGroup) = [CyclicElem(c.n, 1)]
PermutationGroups.order(::Type{T}, c::CyclicGroup) where {T} = convert(T, c.n)
function Base.iterate(c::CyclicGroup, prev::CyclicElem=CyclicElem(c.n, -1))
    id = prev.id + 1
    if id >= c.n
        return nothing
    else
        next = CyclicElem(c.n, id)
        return next, next
    end
end

Now we define that the cyclic group acts on monomial by permuting variables cyclically. So for instance, CyclicElem(3, 1) would transform x_1^3*x_2*x_3^4 into x_1^4*x_2^3*x_3.

import MultivariatePolynomials as MP

using SumOfSquares

struct Action{V<:MP.AbstractVariable} <: Symmetry.OnMonomials
    variables::Vector{V}
end
Symmetry.SymbolicWedderburn.coeff_type(::Action) = Float64
function Symmetry.SymbolicWedderburn.action(a::Action, el::CyclicElem, mono::MP.AbstractMonomial)
    return prod(MP.powers(mono), init=MP.constant_monomial(mono)) do (var, exp)
        index = findfirst(isequal(var), a.variables)
        new_index = mod1(index + el.id, el.n)
        return a.variables[new_index]^exp
    end
end

using DynamicPolynomials
@polyvar x[1:3]
action = Action(x)
g = CyclicElem(3, 1)
Symmetry.SymbolicWedderburn.action(action, g, x[1]^3 * x[2] * x[3]^4)

\[ x_{1}^{4}x_{2}^{3}x_{3} \]

The following polynomial poly is invariant under the action of the group G.

N = 3
G = CyclicGroup(N)
poly = sum(x[i] * x[mod1(i + 1, N)] for i in 1:N) + sum(x.^2)
Symmetry.SymbolicWedderburn.action(action, g, poly)

\[ x_{3}^{2} + x_{2}x_{3} + x_{2}^{2} + x_{1}x_{3} + x_{1}x_{2} + x_{1}^{2} \]

Let's now find the minimum of p by exploiting this symmetry.

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

* Solver : CSDP

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Problem solved to optimality."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 2.22045e-16
  Dual objective value : -2.20229e-10

* Work counters
  Solve time (sec)   : 1.49012e-03

Let's look at the symmetry adapted basis used.

for gram in gram_matrix(con_ref).blocks
    println(gram.basis.polynomials)
    display(gram.Q)
end

DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[1.0, -0.5773502691896257x₃ - 0.5773502691896257x₂ - 0.5773502691896257x₁]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-0.816496580927726x₃ + 0.4082482904638631x₂ + 0.4082482904638631x₁]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-0.7071067811865477x₂ + 0.7071067811865477x₁]

Let's look into more details at the last two elements of the basis.

basis = [(x[1] + x[2] - 2x[3])/√6, (x[1] - x[2])/√2]

2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:
 -0.8164965809277261x₃ + 0.4082482904638631x₂ + 0.4082482904638631x₁
 -0.7071067811865475x₂ + 0.7071067811865475x₁

This actually constitutes the basis for an invariant subspace corresponding to a group character of degree 2 and multiplicity 1. This means that it decomposes the semidefinite matrix into 2 blocks of size 1-by-1 that are equal. Indeed, we see above that gram.Q is identically equal for both. As the group is generated by one element g, we can just verify it by verifying its invariance under g. The image of each element under the basis is:

image = [Symmetry.SymbolicWedderburn.action(action, g, p) for p in basis]

2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:
 0.4082482904638631x₃ + 0.4082482904638631x₂ - 0.8164965809277261x₁
 -0.7071067811865475x₃ + 0.7071067811865475x₂

We can see that they are both still in the same 2-dimensional subspace.

a = -1/2
b = √3/2
[a -b; b a] * basis

2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:
 0.4082482904638631x₃ + 0.4082482904638629x₂ - 0.816496580927726x₁
 -0.7071067811865476x₃ + 0.7071067811865475x₂ + 5.551115123125783e-17x₁

In fact, these last two basis comes from the real decomposition of a complex one.

import CSDP
solver = CSDP.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, action)
cone = SumOfSquares.NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}()
con_ref = @constraint(model, poly - t in cone, symmetry = pattern)
optimize!(model)
solution_summary(model)

for gram in gram_matrix(con_ref).blocks
    println(gram.basis.polynomials)
    display(gram.Q)
end

Iter:  8 Ap: 1.00e+00 Pobj:  2.2204460e-16 Ad: 9.00e-01 Dobj: -2.2022945e-10
Success: SDP solved
Primal objective value: 2.2204460e-16
Dual objective value: -2.2022945e-10
Relative primal infeasibility: 8.86e-16
Relative dual infeasibility: 3.82e-10
Real Relative Gap: -2.20e-10
XZ Relative Gap: 2.29e-09
DIMACS error measures: 1.53e-15 0.00e+00 7.89e-10 0.00e+00 -2.20e-10 2.29e-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: 1.00e+00 Pobj: -8.9892401e+00 Ad: 9.07e-01 Dobj:  4.8627163e+01
Iter:  2 Ap: 1.00e+00 Pobj: -5.4165963e+00 Ad: 9.18e-01 Dobj:  1.0591774e+01
Iter:  3 Ap: 1.00e+00 Pobj: -6.3166907e-01 Ad: 8.98e-01 Dobj:  3.0288603e+00
Iter:  4 Ap: 1.00e+00 Pobj: -4.2564141e-02 Ad: 9.09e-01 Dobj:  3.9201832e-01
Iter:  5 Ap: 1.00e+00 Pobj: -4.3454962e-03 Ad: 9.28e-01 Dobj:  4.0461107e-02
Iter:  6 Ap: 1.00e+00 Pobj: -4.2815474e-04 Ad: 9.34e-01 Dobj:  3.8633036e-03
Iter:  7 Ap: 1.00e+00 Pobj: -4.3920196e-05 Ad: 9.88e-01 Dobj:  1.7392285e-04
Iter:  8 Ap: 1.00e+00 Pobj: -2.3698636e-06 Ad: 9.99e-01 Dobj:  4.5447245e-06
Iter:  9 Ap: 1.00e+00 Pobj: -9.4705116e-08 Ad: 1.00e+00 Dobj: -2.3966141e-06
Iter: 10 Ap: 1.00e+00 Pobj: -4.0174413e-09 Ad: 9.69e-01 Dobj: -9.9538646e-08
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(1.0 - 0.0im), (-0.5773502691896257 - 0.0im)x₃ + (-0.5773502691896257 - 0.0im)x₂ + (-0.5773502691896257 - 0.0im)x₁]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(-0.577350269189626 - 0.0im)x₃ + (0.2886751345948128 + 0.5000000000000001im)x₂ + (0.28867513459481314 - 0.49999999999999994im)x₁]
DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(-0.577350269189626 - 0.0im)x₃ + (0.28867513459481314 - 0.49999999999999994im)x₂ + (0.2886751345948128 + 0.5000000000000001im)x₁]

We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:

complex_basis = basis[1] + im * basis[2]
image = Symmetry.SymbolicWedderburn.action(action, g, complex_basis)

\[ (0.4082482904638631 - 0.7071067811865475im)x_{3} + (0.4082482904638631 + 0.7071067811865475im)x_{2} + (-0.8164965809277261 + 0.0im)x_{1} \]

And there is a direct correspondance between the representation of the real and complex versions:

(a + b * im) * complex_basis

\[ (0.4082482904638631 - 0.7071067811865476im)x_{3} + (0.4082482904638629 + 0.7071067811865475im)x_{2} + (-0.816496580927726 + 5.551115123125783e-17im)x_{1} \]

This page was generated using Literate.jl.