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
endNow 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
import MultivariateBases as MB
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 Clarabel
solver = Clarabel.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, action)
basis = MB.explicit_basis(MB.algebra_element(poly - t))
using SymbolicWedderburn
summands = SymbolicWedderburn.symmetry_adapted_basis(
Float64,
pattern.group,
pattern.action,
basis,
semisimple = true,
)
gram_basis = SumOfSquares.Certificate.Symmetry._gram_basis(pattern, basis, Float64)
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
optimize!(model)
solution_summary(model)* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -7.30048e-11
Dual objective value : 2.25773e-10
* Work counters
Solve time (sec) : 6.65317e-04
Barrier iterations : 5
Let's look at the symmetry adapted basis used.
gram_matrix(con_ref)BlockDiagonalGramMatrix with 2 blocks:
[1] Block with row/column basis:
Simple basis:
FixedBasis([1.0·1 + 0.0·x₃ + 0.0·x₂ + 0.0·x₁, 0.0·1 - 0.5773502691896257·x₃ - 0.5773502691896257·x₂ - 0.5773502691896257·x₁])
And entries in a 2×2 SymMatrix{Float64}:
7.300482440797396e-11 0.0
0.0 2.0
[2] Block with row/column basis:
Semisimple basis with 2 simple sub-bases:
FixedBasis([0.0·1 - 0.816496580927726·x₃ + 0.4082482904638631·x₂ + 0.4082482904638631·x₁])
FixedBasis([0.0·1 + 0.0·x₃ - 0.7071067811865477·x₂ + 0.7071067811865477·x₁])
And entries in a 1×1 SymMatrix{Float64}:
0.4999999999999991Let'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] * basis2-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.
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)
gram_matrix(con_ref)BlockDiagonalGramMatrix with 3 blocks:
[1] Block with row/column basis:
Simple basis:
FixedBasis([(1.0 - 0.0im)·1 + (0.0 + 0.0im)·x₃ + (0.0 + 0.0im)·x₂ + (0.0 + 0.0im)·x₁, (0.0 + 0.0im)·1 + (-0.5773502691896257 - 0.0im)·x₃ + (-0.5773502691896257 - 0.0im)·x₂ + (-0.5773502691896257 - 0.0im)·x₁])
And entries in a 2×2 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:
8.099519427648694e-10 + 0.0im 0.0 + 0.0im
0.0 + 0.0im 2.0000000000000013 + 0.0im
[2] Block with row/column basis:
Simple basis:
FixedBasis([(0.0 + 0.0im)·1 + (-0.577350269189626 - 0.0im)·x₃ + (0.28867513459481314 - 0.49999999999999994im)·x₂ + (0.2886751345948128 + 0.5000000000000001im)·x₁])
And entries in a 1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:
0.4999999999999958 + 0.0im
[3] Block with row/column basis:
Simple basis:
FixedBasis([(0.0 + 0.0im)·1 + (-0.577350269189626 - 0.0im)·x₃ + (0.2886751345948128 + 0.5000000000000001im)·x₂ + (0.28867513459481314 - 0.49999999999999994im)·x₁])
And entries in a 1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:
0.5000000000000036 + 0.0imWe 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.