Fix symmetry reduction

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

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
 

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

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}()))