From 6cb02719633bd6a9ddc733cdafe90a5a78533f1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 8 Jul 2024 10:36:55 +0200 Subject: [PATCH] Fix symmetry reduction --- docs/src/tutorials/Symmetry/cyclic.jl | 14 ++++++ src/Certificate/Symmetry/Symmetry.jl | 1 + src/Certificate/Symmetry/wedderburn.jl | 63 +++++++++++++++++++------- src/Certificate/ideal.jl | 2 +- 4 files changed, 63 insertions(+), 17 deletions(-) diff --git a/docs/src/tutorials/Symmetry/cyclic.jl b/docs/src/tutorials/Symmetry/cyclic.jl index f7088f4b0..6cc61c777 100644 --- a/docs/src/tutorials/Symmetry/cyclic.jl +++ b/docs/src/tutorials/Symmetry/cyclic.jl @@ -51,6 +51,7 @@ end # `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 @@ -89,6 +90,19 @@ model = Model(solver) @variable(model, t) @objective(model, Max, t) pattern = Symmetry.Pattern(G, action) +#import MultivariateBases as MB +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) @test termination_status(model) == MOI.OPTIMAL #src diff --git a/src/Certificate/Symmetry/Symmetry.jl b/src/Certificate/Symmetry/Symmetry.jl index 9a49c76d1..c6119c288 100644 --- a/src/Certificate/Symmetry/Symmetry.jl +++ b/src/Certificate/Symmetry/Symmetry.jl @@ -3,6 +3,7 @@ module Symmetry import LinearAlgebra import MutableArithmetics as MA +import StarAlgebras as SA import MultivariatePolynomials as MP import MultivariateBases as MB diff --git a/src/Certificate/Symmetry/wedderburn.jl b/src/Certificate/Symmetry/wedderburn.jl index 3cb73a757..a34d52432 100644 --- a/src/Certificate/Symmetry/wedderburn.jl +++ b/src/Certificate/Symmetry/wedderburn.jl @@ -1,3 +1,18 @@ +function SymbolicWedderburn.decompose( + p::MB.Polynomial, + hom::SymbolicWedderburn.InducedActionHomomorphism, +) + return [hom[p]], [1] +end + +function SymbolicWedderburn.decompose( + p::SA.AlgebraElement, + hom::SymbolicWedderburn.InducedActionHomomorphism, +) + return [hom[k] for k in SA.supp(p)], + [v for (_, v) in SA.nonzero_pairs(SA.coeffs(p))] +end + function SymbolicWedderburn.decompose( k::MP.AbstractPolynomialLike, hom::SymbolicWedderburn.InducedActionHomomorphism, @@ -10,13 +25,13 @@ function SymbolicWedderburn.decompose( return indcs, coeffs end -function SymbolicWedderburn.ExtensionHomomorphism( - action::SymbolicWedderburn.Action, - basis::MB.SubBasis{MB.Monomial}, -) - monos = collect(basis.monomials) - return SymbolicWedderburn.ExtensionHomomorphism(Int, action, monos) -end +#function SymbolicWedderburn.ExtensionHomomorphism( +# action::SymbolicWedderburn.Action, +# basis::MB.SubBasis{MB.Monomial}, +#) +# monos = collect(basis.monomials) +# return SymbolicWedderburn.ExtensionHomomorphism(Int, action, monos) +#end struct VariablePermutation <: SymbolicWedderburn.ByPermutations end _map_idx(f, v::AbstractVector) = map(f, eachindex(v)) @@ -35,6 +50,18 @@ function SymbolicWedderburn.action( end abstract type OnMonomials <: SymbolicWedderburn.ByLinearTransformation end +function SymbolicWedderburn.action( + a::Union{VariablePermutation,OnMonomials}, + el, + p::MB.Polynomial{MB.Monomial}, +) + res = SymbolicWedderburn.action(a, el, p.monomial) + if res isa MP.AbstractMonomial + return MB.Polynomial{MB.Monomial}(res) + else + return MB.algebra_element(res) + end +end function SymbolicWedderburn.action( a::Union{VariablePermutation,OnMonomials}, el, @@ -76,7 +103,7 @@ function SumOfSquares.matrix_cone_type(::Type{<:Ideal{C}}) where {C} return SumOfSquares.matrix_cone_type(C) end function SumOfSquares.Certificate.gram_basis_type(::Type{<:Ideal}) - return Vector{Vector{MB.FixedPolynomialBasis}} + return MB.MultiBasis end SumOfSquares.Certificate.zero_basis_type(::Type{<:Ideal}) = MB.Monomial SumOfSquares.Certificate.zero_basis(::Ideal) = MB.Monomial @@ -108,12 +135,12 @@ function SumOfSquares.Certificate.reduced_basis( end function MA.promote_operation( ::typeof(SumOfSquares.Certificate.reduced_basis), - ::Type{Ideal{S,C}}, + ::Type{<:Ideal{C}}, ::Type{B}, ::Type{D}, ::Type{G}, ::Type{W}, -) where {S,C,B,D,G,W} +) where {C,B,D,G,W} return MA.promote_operation( SumOfSquares.Certificate.reduced_basis, C, @@ -149,6 +176,7 @@ end function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T} # We set `semisimple=true` as we don't support simple yet since it would not give all the simple components but only one of them. + @show T summands = SymbolicWedderburn.symmetry_adapted_basis( T, pattern.group, @@ -223,14 +251,17 @@ function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T} # Moreover, `Q = kron(C, I)` is not block diagonal but we can get a block-diagonal # `Q = kron(I, Q)` by permuting the rows and columns: # `(U[1:d:(1+d*(m-1))] * F)' * basis.monomials`, `(U[2:d:(2+d*(m-1))] * F)' * basis.monomials`, ... - map(1:d) do i - return MB.FixedPolynomialBasis( - (transpose(U[:, i:d:(i+d*(m-1))]) * F) * basis.monomials, - ) - end + return MB.MultiBasis( + map(1:d) do i + return MB.OrthonormalCoefficientsBasis( + transpose(U[:, i:d:(i+d*(m-1))]) * F, + basis, + ) + end, + ) else F = convert(Matrix{T}, R) - [MB.FixedPolynomialBasis(F * basis.monomials)] + MB.MultiBasis([MB.OrthonormalCoefficientsBasis(F, basis)]) end end end diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl index c95b6455e..21a0983e7 100644 --- a/src/Certificate/ideal.jl +++ b/src/Certificate/ideal.jl @@ -16,7 +16,7 @@ Base.:+(::_NonZero, a::_NonZero) = a function _combine_with_gram( basis::MB.SubBasis{B,M}, - gram_bases::AbstractVector{<:MB.SubBasis}, + gram_bases::AbstractVector{<:SA.ExplicitBasis}, weights, ) where {B,M} p = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}()))