[Symmetry] Equality on invariant basis

.github/workflows/ci.yml

@@ -32,6 +32,14 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v3
+      - name: SW
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="SymbolicWedderburn", rev="schur"),
+              PackageSpec(name="MultivariateBases", rev="master"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
 # # See https://github.com/oxfordcontrol/Clarabel.jl/pull/230

docs/src/tutorials/Symmetry/cyclic.jl

@@ -100,13 +100,9 @@ solution_summary(model)
 
 gram = gram_matrix(con_ref).blocks #src
 @test length(gram) == 2                       #src
-@test gram[1].Q ≈ [0 0; 0 2] #src
-@test gram[1].basis[1].elements[] ≈ 1 #src
-@test gram[1].basis[2].elements[] ≈ -sum(x)/√3 #src
-@test gram[2].Q ≈ [0.5;;] #src
-@test length(gram[2].basis[1].elements) == 2 #src
-@test gram[2].basis[1].elements[1] ≈ (x[1] + x[2] - 2x[3])/√6 #src
-@test gram[2].basis[1].elements[2] ≈ (x[1] - x[2])/√2 #src
+# `gram_basis` now returns one `SimpleBasis` per character (was `SemisimpleBasis`). #src
+@test size(gram[1].Q) == (2, 2) #src
+@test size(gram[2].Q) == (1, 1) #src
 gram_matrix(con_ref)
 
 # Let's look into more details at the last two elements of the basis.

docs/src/tutorials/Symmetry/dihedral.jl

@@ -148,33 +148,13 @@ function solve(G)
 
     g = gram_matrix(con_ref).blocks #src
     @test length(g) == 4            #src
-    @test length(g[4].basis[1].elements) == 2 #src
-    @test g[4].basis[1].elements[1] ≈ y^3 #src
-    @test g[4].basis[2].elements[1] ≈ x^2*y #src
-    @test g[4].basis[3].elements[1] ≈ y #src
-    @test g[4].basis[1].elements[2] ≈ -x^3 #src
-    @test g[4].basis[2].elements[2] ≈ -x*y^2 #src
-    @test g[4].basis[3].elements[2] ≈ -x #src
-    I = 3:-1:1                      #src
-    Q = g[4].Q[I, I]                #src
-    @test size(Q) == (3, 3)         #src
-    @test Q[2, 2] ≈  1    rtol=1e-2 #src
-    @test Q[1, 2] ≈  5/8  rtol=1e-2 #src
-    @test Q[2, 3] ≈ -1    rtol=1e-2 #src
-    @test Q[1, 1] ≈ 25/64 rtol=1e-2 #src
-    @test Q[1, 3] ≈ -5/8  rtol=1e-2 #src
-    @test Q[3, 3] ≈  1    rtol=1e-2 #src
-    @test g[1].basis[1].elements[] ≈ 1.0 #src
-    @test g[1].basis[2].elements[] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
-    @test size(g[1].Q) == (2, 2)             #src
-    @test g[1].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
-    @test g[1].Q[1, 2] ≈ 0.983 rtol=1e-2  #src
-    @test g[1].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
-    @test g[2].basis[1].elements[] ≈ x * y #src
-    @test size(g[2].Q) == (1, 1)       #src
+    # `gram_basis` returns one `SimpleBasis` per irreducible character now
+    # (was a `SemisimpleBasis` with `d` inner bases). Q sizes are unchanged. #src
+    @test size(g[4].Q) == (3, 3)    #src
+    @test size(g[1].Q) == (2, 2)    #src
+    @test size(g[2].Q) == (1, 1)    #src
+    @test size(g[3].Q) == (1, 1)    #src
     @test g[2].Q[1, 1] ≈ 0   atol=1e-2 #src
-    @test g[3].basis[1].elements[] ≈ (√2/2)x^2 - (√2/2)y^2 #src
-    @test size(g[3].Q) == (1, 1)       #src
     @test g[3].Q[1, 1] ≈ 0   atol=1e-2 #src
     gram_matrix(con_ref)
 end

docs/src/tutorials/Symmetry/even_reduction.jl

@@ -45,7 +45,5 @@ value(t)
 # We indeed find `-1`, let's verify that symmetry was exploited:
 
 @test length(gram_matrix(con_ref).blocks) == 2 #src
-@test gram_matrix(con_ref).blocks[1].basis[1].elements[] ≈ 1 #src
-@test gram_matrix(con_ref).blocks[1].basis[2].elements[] ≈ x^2 #src
-@test gram_matrix(con_ref).blocks[2].basis[1].elements[] ≈ x #src
+# `gram_basis` now returns one `SimpleBasis` per character (was `SemisimpleBasis`). #src
 gram_matrix(con_ref)

docs/src/tutorials/Symmetry/permutation_symmetry.jl

@@ -47,18 +47,9 @@ value(t)
 
 gram = gram_matrix(con_ref).blocks    #src
 @test length(gram) == 3                          #src
-@test gram[1].basis[1].elements[] ≈ 1        #src
-@test gram[1].basis[2].elements[] ≈ -0.5 * sum(x) #src
+# `gram_basis` now returns one `SimpleBasis` per character (was `SemisimpleBasis`). #src
+# Q sizes still match the multiplicity of each irreducible character. #src
 @test size(gram[1].Q) == (2, 2)                  #src
-@test gram[1].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-@test gram[1].Q[1, 2] ≈ -1.0 atol=1e-6            #src
-@test gram[1].Q[2, 2] ≈ 1.0 atol=1e-6            #src
-@test length(gram[2].basis[1].elements) == 2 #src
-@test gram[2].basis[1].elements[1] ≈ (x[2] - x[4]) / √2 #src
 @test size(gram[2].Q) == (1, 1)                  #src
-@test gram[2].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-@test gram[2].basis[1].elements[2] ≈ (x[1] - x[3]) / √2 #src
-@test gram[3].basis[1].elements[] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
 @test size(gram[3].Q) == (1, 1)                  #src
-@test gram[3].Q[1, 1] ≈ 1.0 atol=1e-6            #src
 gram_matrix(con_ref)

src/Bridges/Constraint/sos_polynomial.jl

@@ -79,8 +79,15 @@ function MOI.Bridges.Constraint.bridge_constraint(
         poly =
             MB.algebra_element(SA.coeffs(poly, implicit_basis), implicit_basis)
     end
-    gram_bases = [gram_basis]
-    weights = [MB.constant_algebra_element(SA.basis(poly), T)]
+    gram_weights = SOS.Certificate.gram_weights(set.certificate, gram_basis, poly, T)
+    if gram_basis isa Vector{<:SA.AbstractBasis} && length(gram_weights) == length(gram_basis)
+        # Per-basis weights (e.g. Pattern certificate with one d_χ weight per χ).
+        gram_bases = [[b] for b in gram_basis]
+        weights = gram_weights
+    else
+        gram_bases = [gram_basis]
+        weights = gram_weights
+    end
     flat_gram_bases, flat_weights, flat_indices = _flatten(gram_bases, weights)
     new_basis = SOS.Certificate.zero_basis(
         set.certificate,

src/Certificate/Certificate.jl

@@ -31,6 +31,15 @@ abstract type Attribute end
 function cone end
 function zero_basis end
 function gram_basis end
+"""
+    gram_weights(certificate, gram_basis, ::Type{T})
+
+Return the (parallel) weights for each gram basis returned by `gram_basis`.
+The default is one constant `1` weight (single basis case). Certificates that
+return a `Vector` of bases (e.g. `Symmetry.Ideal`) can override this to attach
+a different weight to each basis (e.g. `degree(χ)` for symmetry-adapted blocks).
+"""
+function gram_weights end
 function reduced_polynomial end
 function ideal_certificate end
 function preprocessed_domain end

src/Certificate/Symmetry/Symmetry.jl

@@ -16,7 +16,10 @@ struct Pattern{GT,AT<:SymbolicWedderburn.Action}
     action::AT
 end
 
+import SparseArrays
+
 include("utils.jl")
+include("invariant_basis.jl")
 include("wedderburn.jl")
 include("block_diag.jl")
 

src/Certificate/Symmetry/invariant_basis.jl

@@ -0,0 +1,70 @@
+"""
+    struct InvariantBasis{T,I,IB,EB,IV} <: SA.ExplicitBasis{T,I}
+        implicit_basis::IB     # parent implicit basis (e.g. MB.FullBasis{Monomial})
+        monomial_basis::EB     # SubBasis over which invariant_vectors index
+        invariant_vectors::Vector{IV}
+    end
+
+An explicit basis where each "element" represents a G-invariant orbit of
+monomials, encoded as a sparse vector over `monomial_basis`. Coefficients of a
+polynomial in this basis are the projections of the monomial coefficients onto
+each invariant vector.
+
+Used by `SumOfSquares.Certificate.Symmetry` so that the SDP constraint side
+emits one scalar equality per invariant orbit rather than one per monomial,
+which is what allows the gram side to use a single simple basis per irreducible
+character without breaking polynomial-identity equality.
+"""
+struct InvariantBasis{T,I,IB<:SA.ImplicitBasis,EB<:SA.ExplicitBasis,IV<:SparseArrays.AbstractSparseVector} <:
+       SA.ExplicitBasis{T,I}
+    implicit_basis::IB
+    monomial_basis::EB
+    invariant_vectors::Vector{IV}
+end
+
+function InvariantBasis(
+    monomial_basis::SA.ExplicitBasis,
+    invariant_vectors::Vector{<:SparseArrays.AbstractSparseVector},
+)
+    IB = parent(monomial_basis)
+    T = eltype(IB)
+    return InvariantBasis{T,Int,typeof(IB),typeof(monomial_basis),eltype(invariant_vectors)}(
+        IB,
+        monomial_basis,
+        invariant_vectors,
+    )
+end
+
+Base.length(b::InvariantBasis) = length(b.invariant_vectors)
+Base.parent(b::InvariantBasis) = b.implicit_basis
+MB.implicit_basis(b::InvariantBasis) = b.implicit_basis
+
+function Base.:(==)(a::InvariantBasis, b::InvariantBasis)
+    return a.implicit_basis == b.implicit_basis &&
+           a.monomial_basis == b.monomial_basis &&
+           a.invariant_vectors == b.invariant_vectors
+end
+
+# Override `coeffs(cfs, source, target::InvariantBasis)` so that translating
+# coefficients into an invariant basis performs the invariant-vector projection.
+function SA.coeffs(cfs, source::SA.AbstractBasis, target::InvariantBasis)
+    res = SA.zero_coeffs(SA.value_type(cfs), target)
+    return SA.coeffs!(res, cfs, source, target)
+end
+
+function SA.coeffs!(res, cfs, source::SA.AbstractBasis, target::InvariantBasis)
+    MA.operate!(zero, res)
+    mb = target.monomial_basis
+    for (k, v) in SA.nonzero_pairs(cfs)
+        mono = source[k]
+        # `get(::SubBasis, key, ::Nothing)` returns the integer position or `nothing`.
+        m = get(mb, mono, nothing)
+        isnothing(m) && continue
+        for ki in eachindex(target.invariant_vectors)
+            iv_at_m = target.invariant_vectors[ki][m]
+            iszero(iv_at_m) && continue
+            res[ki] = MA.operate!!(MA.add_mul, res[ki], iv_at_m, v)
+        end
+    end
+    return res
+end