Fix symmetry reduction (#375)

* Fix symmetry reduction * Fix * Fixes * Fixes * Fixes * Reenable Symmetry * Fixes * Fix * Remove commented code * Fix format * Updates * Update * Fix * Fix format
jump-dev · Dec 15, 2024 · 2b8f0fb · 2b8f0fb
1 parent 28f42fa
commit 2b8f0fb
Show file tree

Hide file tree

Showing 9 changed files with 241 additions and 125 deletions.
diff --git a/docs/make.jl b/docs/make.jl
@@ -13,7 +13,7 @@ const _TUTORIAL_SUBDIR = [
     "Other Applications",
     "Noncommutative and Hermitian",
     "Sparsity",
-    #"Symmetry",
+    "Symmetry",
     "Extension",
 ]
 

diff --git 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
 
@@ -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,21 +111,19 @@ 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
-@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
+polys = gram[1].basis.bases[].elements #src
+@test length(polys) == 2 #src
+@test polys[1] ≈ 1 #src
+@test polys[2] ≈ -sum(x)/√3 #src
 @test gram[2].Q ≈ [0.5;;] #src
-@test length(gram[2].basis.polynomials) == 1 #src
-@test gram[2].basis.polynomials[1] ≈ (x[1] + x[2] - 2x[3])/√6 #src
-@test gram[3].Q == gram[2].Q #src
-@test length(gram[3].basis.polynomials) == 1 #src
-@test gram[3].basis.polynomials[1] ≈ (x[1] - x[2])/√2 #src
-for gram in gram_matrix(con_ref).blocks
-    println(gram.basis.polynomials)
-    display(gram.Q)
-end
+@test length(gram[2].basis.bases) == 2 #src
+polys = gram[2].basis.bases[1].elements #src
+@test polys[] ≈ (x[1] + x[2] - 2x[3])/√6 #src
+polys = gram[2].basis.bases[2].elements #src
+@test polys[] ≈ (x[1] - x[2])/√2 #src
+gram_matrix(con_ref)
 
 # Let's look into more details at the last two elements of the basis.
 
@@ -137,8 +148,6 @@ b = √3/2
 
 # In fact, these last two basis comes from the real decomposition of a complex one.
 
-import CSDP
-solver = CSDP.Optimizer
 model = Model(solver)
 @variable(model, t)
 @objective(model, Max, t)
@@ -153,19 +162,19 @@ solution_summary(model)
 gram = gram_matrix(con_ref).blocks #src
 @test length(gram) == 3                       #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
-@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
+polys = gram[1].basis.bases[].elements #src
+@test length(polys) == 2 #src
+@test polys[1] ≈ 1 #src
+@test polys[2] ≈ -sum(x)/√3 #src
 @test gram[2].Q ≈ [0.5;;] rtol = 1e-6 #src
-@test length(gram[2].basis.polynomials) == 1 #src
-@test gram[2].basis.polynomials[1] ≈ (basis[1] - basis[2] * im) / √2  #src
+polys = gram[2].basis.bases[].elements #src
+@test length(polys) == 1 #src
+@test polys[] ≈ (basis[1] + basis[2] * im) / √2  #src
 @test gram[3].Q ≈ [0.5;;] rtol = 1e-6 #src
-@test length(gram[3].basis.polynomials) == 1 #src
-@test gram[3].basis.polynomials[1] ≈ (basis[1] + basis[2] * im) / √2 #src
-for gram in gram_matrix(con_ref).blocks
-    println(gram.basis.polynomials)
-    display(gram.Q)
-end
+polys = gram[3].basis.bases[].elements #src
+@test length(polys) == 1 #src
+@test polys[] ≈ (basis[1] - basis[2] * im) / √2 #src
+gram_matrix(con_ref)
 
 # We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:
 

diff --git a/docs/src/tutorials/Symmetry/dihedral.jl b/docs/src/tutorials/Symmetry/dihedral.jl
@@ -8,7 +8,6 @@
 # *Symmetry groups, semidefinite programs, and sums of squares*.
 # Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128.
 
-
 using Test #src
 import Random #src
 Random.seed!(0) #src
@@ -134,9 +133,9 @@ end
 
 # We can exploit this symmetry for reducing the problem using the `SymmetricIdeal` certificate as follows:
 
-import CSDP
+import Clarabel
 function solve(G)
-    solver = CSDP.Optimizer
+    solver = Clarabel.Optimizer
     model = Model(solver)
     @variable(model, t)
     @objective(model, Max, t)
@@ -148,37 +147,44 @@ function solve(G)
 
 
     g = gram_matrix(con_ref).blocks #src
-    @test length(g) == 5                       #src
-    @test g[1].basis.polynomials == [y^3, x^2*y, y] #src
-    @test g[2].basis.polynomials == [-x^3, -x*y^2, -x] #src
-    for i in 1:2                        #src
-        I = 3:-1:1                      #src
-        Q = g[i].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
-    end #src
-    @test length(g[3].basis.polynomials) == 2 #src
-    @test g[3].basis.polynomials[1] == 1.0 #src
-    @test g[3].basis.polynomials[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
-    @test size(g[3].Q) == (2, 2)             #src
-    @test g[3].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
-    @test g[3].Q[1, 2] ≈ 0.983 rtol=1e-2  #src
-    @test g[3].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
-    @test g[4].basis.polynomials == [x * y] #src
-    @test size(g[4].Q) == (1, 1)       #src
-    @test g[4].Q[1, 1] ≈ 0   atol=1e-2 #src
-    @test length(g[5].basis.polynomials) == 1 #src
-    @test g[5].basis.polynomials[1] ≈ (√2/2)x^2 - (√2/2)y^2 #src
-    @test size(g[5].Q) == (1, 1)       #src
-    @test g[5].Q[1, 1] ≈ 0   atol=1e-2 #src
-    for g in gram_matrix(con_ref).blocks
-        println(g.basis.polynomials)
-    end
+    @test length(g) == 4            #src
+    @test length(g[4].basis.bases) == 2 #src
+    polys = g[4].basis.bases[1].elements #src
+    @test length(polys) == 3 #src
+    @test polys[1] ≈ y^3 #src
+    @test polys[2] ≈ x^2*y #src
+    @test polys[3] ≈ y #src
+    polys = g[4].basis.bases[2].elements #src
+    @test length(polys) == 3 #src
+    @test polys[1] ≈ -x^3 #src
+    @test polys[2] ≈ -x*y^2 #src
+    @test polys[3] ≈ -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
+    polys = g[1].basis.bases[].elements #src
+    @test length(polys) == 2 #src
+    @test polys[1] ≈ 1.0 #src
+    @test polys[2] ≈ -(√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
+    polys = g[2].basis.bases[].elements #src
+    @test polys[] ≈ x * y #src
+    @test size(g[2].Q) == (1, 1)       #src
+    @test g[2].Q[1, 1] ≈ 0   atol=1e-2 #src
+    polys = g[3].basis.bases[].elements #src
+    @test polys[] ≈ (√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
 solve(G)
 

diff --git a/docs/src/tutorials/Symmetry/even_reduction.jl b/docs/src/tutorials/Symmetry/even_reduction.jl
@@ -31,8 +31,8 @@ G = PermGroup([perm"(1,2)"])
 
 # We can exploit the symmetry as follows:
 
-import CSDP
-solver = CSDP.Optimizer
+import Clarabel
+solver = Clarabel.Optimizer
 model = Model(solver)
 @variable(model, t)
 @objective(model, Max, t)
@@ -45,6 +45,9 @@ 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.polynomials == [1, x^2] #src
-@test gram_matrix(con_ref).blocks[2].basis.polynomials == [x] #src
+polys = gram_matrix(con_ref).blocks[1].basis.bases[].elements #src
+@test polys[1] ≈ 1 #src
+@test polys[2] ≈ x^2 #src
+polys = gram_matrix(con_ref).blocks[2].basis.bases[].elements #src
+@test polys[] ≈ x #src
 gram_matrix(con_ref)
diff --git a/docs/src/tutorials/Symmetry/permutation_symmetry.jl b/docs/src/tutorials/Symmetry/permutation_symmetry.jl
@@ -32,8 +32,8 @@ G = PermGroup([perm"(1,2,3,4)"])
 
 # We can use this certificate as follows:
 
-import CSDP
-solver = CSDP.Optimizer
+import Clarabel
+solver = Clarabel.Optimizer
 model = Model(solver)
 @variable(model, t)
 @objective(model, Max, t)
@@ -45,25 +45,28 @@ value(t)
 
 # We indeed find `-1`, let's verify that symmetry was exploited:
 
-g = gram_matrix(con_ref).blocks    #src
-@test length(g) == 4                          #src
-@test length(g[1].basis.polynomials) == 2     #src
-@test g[1].basis.polynomials[1] == 1.0        #src
-@test g[1].basis.polynomials[2] ≈ -0.5 * sum(x) #src
-@test size(g[1].Q) == (2, 2)                  #src
-@test g[1].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-@test g[1].Q[1, 2] ≈ -1.0 atol=1e-6            #src
-@test g[1].Q[2, 2] ≈ 1.0 atol=1e-6            #src
-@test length(g[2].basis.polynomials) == 1     #src
-@test g[2].basis.polynomials[1] ≈ (x[2] - x[4]) / √2 #src
+gram = gram_matrix(con_ref).blocks    #src
+@test length(gram) == 3                          #src
+polys = gram[1].basis.bases[].elements #src
+@test length(polys) == 2     #src
+@test polys[1] ≈ 1        #src
+@test polys[2] ≈ -0.5 * sum(x) #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.bases) == 2 #src
+polys = gram[2].basis.bases[1].elements #src
+@test length(polys) == 1     #src
+@test polys[] ≈ (x[2] - x[4]) / √2 #src
 @test size(g[2].Q) == (1, 1)                  #src
 @test g[2].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-@test length(g[3].basis.polynomials) == 1     #src
-@test g[3].basis.polynomials[1] ≈ (x[1] - x[3]) / √2 #src
-@test size(g[3].Q) == (1, 1)                  #src
-@test g[3].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-@test length(g[4].basis.polynomials) == 1     #src
-@test g[4].basis.polynomials[1] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
-@test size(g[4].Q) == (1, 1)                  #src
-@test g[4].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-gram_matrix(con_ref).blocks
+polys = gram[2].basis.bases[2].elements #src
+@test length(polys) == 1     #src
+@test polys[1] ≈ (x[1] - x[3]) / √2 #src
+polys = gram[3].basis.bases[].elements #src
+@test length(polys) == 1     #src
+@test polys[] ≈ (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)
diff --git 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