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
 
@@ -83,12 +84,24 @@ Symmetry.SymbolicWedderburn.action(action, g, poly)
 
 # Let's now find the minimum of `p` by exploiting this symmetry.
 
-import CSDP
-solver = CSDP.Optimizer
+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)
 @test termination_status(model) == MOI.OPTIMAL #src
@@ -98,7 +111,7 @@ solution_summary(model)
 # Let's look at the symmetry adapted basis used.
 
 gram = gram_matrix(con_ref).blocks #src
-@test length(gram) == 3                       #src
+@test length(gram) == 2                       #src
 @test gram[1].Q ≈ [0 0; 0 2] #src
 @test length(gram[1].basis.polynomials) == 2 #src
 @test gram[1].basis.polynomials[1] == 1 #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,
@@ -75,17 +102,47 @@ end
 function SumOfSquares.matrix_cone_type(::Type{<:Ideal{C}}) where {C}
     return SumOfSquares.matrix_cone_type(C)
 end
+function _multi_basis_type(::Type{<:MB.SubBasis{B,M}}) where {B,M}
+    T = Float64
+    SC = SA.SparseCoefficients{M,T,Vector{M},Vector{T}}
+    AE = SA.AlgebraElement{
+        MB.Algebra{MB.FullBasis{B,M},B,M},
+        T,
+        SC,
+    }
+    return Vector{MB.SemisimpleBasis{
+        AE,
+        Int,
+        MB.FixedPolynomialBasis{
+            B,
+            M,
+            T,
+            SC,
+        },
+    }}
+end
 function MA.promote_operation(
     ::typeof(SumOfSquares.Certificate.gram_basis),
-    ::Type{<:Ideal},
-)
-    return Vector{Vector{MB.FixedPolynomialBasis}}
+    ::Type{<:Ideal{C}},
+) where {C}
+    return _multi_basis_type(MA.promote_operation(SumOfSquares.Certificate.gram_basis, C))
 end
 function MA.promote_operation(
     ::typeof(SumOfSquares.Certificate.zero_basis),
-    ::Type{<:Ideal},
-)
-    return MB.Monomial
+    ::Type{<:Ideal{C}},
+    ::Type{B},
+    ::Type{D},
+    ::Type{G},
+    ::Type{W},
+) where {C,B,D,G,W}
+    return MA.promote_operation(
+        SumOfSquares.Certificate.zero_basis,
+        C,
+        B,
+        D,
+        G,
+        W,
+    )
 end
 SumOfSquares.Certificate.zero_basis(::Ideal) = MB.Monomial
 function SumOfSquares.Certificate.reduced_polynomial(
@@ -155,8 +212,16 @@ function SumOfSquares.Certificate.gram_basis(cert::Ideal, poly)
     return _gram_basis(cert.pattern, basis, T)
 end
 
+function _fixed_basis(F, basis)
+    return MB.FixedPolynomialBasis([
+        MB.implicit(MB.algebra_element(row, basis))
+        for row in eachrow(F)
+    ])
+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,
@@ -231,14 +296,16 @@ 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.SemisimpleBasis(
+                map(1:d) do i
+                    return _fixed_basis(
+                        transpose(U[:, i:d:(i+d*(m-1))]) * F,
+                        basis,
+                    )
+                end,
+            )
         else
-            F = convert(Matrix{T}, R)
-            [MB.FixedPolynomialBasis(F * basis.monomials)]
+            return MB.SemisimpleBasis([_fixed_basis(convert(Matrix{T}, R), 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}()))

src/gram_matrix.jl

@@ -137,28 +137,32 @@ function GramMatrix(Q::AbstractMatrix, monos::AbstractVector)
     return GramMatrix(Q, MB.SubBasis{MB.Monomial}(sorted_monos), σ)
 end
 
-function _term_element(α, p::MB.Polynomial{B,M}) where {B,M}
-    return MB.algebra_element(
-        MB.sparse_coefficients(MP.term(α, p.monomial)),
-        MB.FullBasis{B,M}(),
+function _gram_operate!(op, p, α, row, col, args::Vararg{Any,N}) where {N}
+    return MA.operate!(
+        op,
+        p,
+        true * row, # TODO `*` shouldn't be defined for group elements so this is a hack
+        α * col,
+        args...,
     )
 end
 
+function _gram_operate!(op, p, α, row::MB.MultiPoly, col::MB.MultiPoly, args::Vararg{Any,N}) where {N}
+    for (r, c) in zip(row.polynomials, col.polynomials)
+        _gram_operate!(op, p, α, r, c, args...)
+    end
+end
+
 function MA.operate!(
     op::SA.UnsafeAddMul{typeof(*)},
     p::SA.AlgebraElement,
     g::GramMatrix,
     args::Vararg{Any,N},
 ) where {N}
-    for col in eachindex(g.basis)
-        for row in eachindex(g.basis)
-            MA.operate!(
-                op,
-                p,
-                _term_element(true, SA.star(g.basis[row])),
-                _term_element(g.Q[row, col], g.basis[col]),
-                args...,
-            )
+    for row in eachindex(g.basis)
+        row_star = SA.star(g.basis[row])
+        for col in eachindex(g.basis)
+            _gram_operate!(op, p, g.Q[row, col], row_star, g.basis[col], args...)
         end
     end
     return p