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)
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)
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)
Iter: 18 Ap: 9.59e-01 Pobj: -3.5512000e-03 Ad: 9.59e-01 Dobj: -3.5512523e-03 Success: SDP solved Primal objective value: -3.5512000e-03 Dual objective value: -3.5512523e-03 Relative primal infeasibility: 1.17e-13 Relative dual infeasibility: 5.36e-10 Real Relative Gap: -5.20e-08 XZ Relative Gap: 1.22e-09 DIMACS error measures: 1.26e-13 0.00e+00 1.44e-09 0.00e+00 -5.20e-08 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: 1.00e+00 Pobj: -6.8633327e+00 Ad: 9.00e-01 Dobj: 1.9835999e+01 Iter: 2 Ap: 1.00e+00 Pobj: -2.8142041e+00 Ad: 9.00e-01 Dobj: 3.4715380e+00 Iter: 3 Ap: 1.00e+00 Pobj: -1.5335395e-02 Ad: 9.00e-01 Dobj: 4.1414777e-01 Iter: 4 Ap: 1.00e+00 Pobj: -2.4302191e-05 Ad: 1.00e+00 Dobj: 4.4865070e-05 Iter: 5 Ap: 1.00e+00 Pobj: -7.8953240e-07 Ad: 1.00e+00 Dobj: -1.2669850e-06 Iter: 6 Ap: 1.00e+00 Pobj: -5.8902884e-08 Ad: 1.00e+00 Dobj: -5.8810369e-07 Iter: 7 Ap: 1.00e+00 Pobj: -8.1684162e-09 Ad: 1.00e+00 Dobj: 1.9520853e-09
* 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.05095e-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₁]
2×2 SymMatrix{Float64}: 1.07382e-17 0.0 0.0 2.0
1×1 SymMatrix{Float64}: 0.4999999999999981
1×1 SymMatrix{Float64}: 0.4999999999999981
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₁]
2×2 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}: 1.67745e-10+0.0im 0.0+0.0im 0.0+0.0im 2.0+0.0im
1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}: 0.5000000408028554 + 0.0im
1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}: 0.4999999591971749 + 0.0im
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)
And there is a direct correspondance between the representation of the real and complex versions:
(a + b * im) * complex_basis
This notebook was generated using Literate.jl.