Skip to content

Commit 2b8f0fb

Browse files
blegatblegat
authored
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
1 parent 28f42fa commit 2b8f0fb

File tree

9 files changed

+241
-125
lines changed

9 files changed

+241
-125
lines changed

‎docs/make.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ const _TUTORIAL_SUBDIR = [
1313
"Other Applications",
1414
"Noncommutative and Hermitian",
1515
"Sparsity",
16-
#"Symmetry",
16+
"Symmetry",
1717
"Extension",
1818
]
1919

‎docs/src/tutorials/Symmetry/cyclic.jl‎

Lines changed: 37 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,7 @@ end
5151
# `x_1^3*x_2*x_3^4` into `x_1^4*x_2^3*x_3`.
5252

5353
import MultivariatePolynomials as MP
54+
import MultivariateBases as MB
5455

5556
using SumOfSquares
5657

@@ -83,12 +84,24 @@ Symmetry.SymbolicWedderburn.action(action, g, poly)
8384

8485
# Let's now find the minimum of `p` by exploiting this symmetry.
8586

86-
import CSDP
87-
solver = CSDP.Optimizer
87+
import Clarabel
88+
solver = Clarabel.Optimizer
8889
model = Model(solver)
8990
@variable(model, t)
9091
@objective(model, Max, t)
9192
pattern = Symmetry.Pattern(G, action)
93+
basis = MB.explicit_basis(MB.algebra_element(poly - t))
94+
using SymbolicWedderburn
95+
summands = SymbolicWedderburn.symmetry_adapted_basis(
96+
Float64,
97+
pattern.group,
98+
pattern.action,
99+
basis,
100+
semisimple = true,
101+
)
102+
103+
gram_basis = SumOfSquares.Certificate.Symmetry._gram_basis(pattern, basis, Float64)
104+
92105
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
93106
optimize!(model)
94107
@test termination_status(model) == MOI.OPTIMAL #src
@@ -98,21 +111,19 @@ solution_summary(model)
98111
# Let's look at the symmetry adapted basis used.
99112

100113
gram = gram_matrix(con_ref).blocks #src
101-
@test length(gram) == 3 #src
114+
@test length(gram) == 2 #src
102115
@test gram[1].Q ≈ [0 0; 0 2] #src
103-
@test length(gram[1].basis.polynomials) == 2 #src
104-
@test gram[1].basis.polynomials[1] == 1 #src
105-
@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
116+
polys = gram[1].basis.bases[].elements #src
117+
@test length(polys) == 2 #src
118+
@test polys[1] ≈ 1 #src
119+
@test polys[2] ≈ -sum(x)/√3 #src
106120
@test gram[2].Q ≈ [0.5;;] #src
107-
@test length(gram[2].basis.polynomials) == 1 #src
108-
@test gram[2].basis.polynomials[1] ≈ (x[1] + x[2] - 2x[3])/√6 #src
109-
@test gram[3].Q == gram[2].Q #src
110-
@test length(gram[3].basis.polynomials) == 1 #src
111-
@test gram[3].basis.polynomials[1] ≈ (x[1] - x[2])/√2 #src
112-
for gram in gram_matrix(con_ref).blocks
113-
println(gram.basis.polynomials)
114-
display(gram.Q)
115-
end
121+
@test length(gram[2].basis.bases) == 2 #src
122+
polys = gram[2].basis.bases[1].elements #src
123+
@test polys[] ≈ (x[1] + x[2] - 2x[3])/√6 #src
124+
polys = gram[2].basis.bases[2].elements #src
125+
@test polys[] ≈ (x[1] - x[2])/√2 #src
126+
gram_matrix(con_ref)
116127

117128
# Let's look into more details at the last two elements of the basis.
118129

@@ -137,8 +148,6 @@ b = √3/2
137148

138149
# In fact, these last two basis comes from the real decomposition of a complex one.
139150

140-
import CSDP
141-
solver = CSDP.Optimizer
142151
model = Model(solver)
143152
@variable(model, t)
144153
@objective(model, Max, t)
@@ -153,19 +162,19 @@ solution_summary(model)
153162
gram = gram_matrix(con_ref).blocks #src
154163
@test length(gram) == 3 #src
155164
@test gram[1].Q ≈ [0 0; 0 2] #src
156-
@test length(gram[1].basis.polynomials) == 2 #src
157-
@test gram[1].basis.polynomials[1] == 1 #src
158-
@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
165+
polys = gram[1].basis.bases[].elements #src
166+
@test length(polys) == 2 #src
167+
@test polys[1] ≈ 1 #src
168+
@test polys[2] ≈ -sum(x)/√3 #src
159169
@test gram[2].Q ≈ [0.5;;] rtol = 1e-6 #src
160-
@test length(gram[2].basis.polynomials) == 1 #src
161-
@test gram[2].basis.polynomials[1] ≈ (basis[1] - basis[2] * im) / √2 #src
170+
polys = gram[2].basis.bases[].elements #src
171+
@test length(polys) == 1 #src
172+
@test polys[] ≈ (basis[1] + basis[2] * im) / √2 #src
162173
@test gram[3].Q ≈ [0.5;;] rtol = 1e-6 #src
163-
@test length(gram[3].basis.polynomials) == 1 #src
164-
@test gram[3].basis.polynomials[1] ≈ (basis[1] + basis[2] * im) / √2 #src
165-
for gram in gram_matrix(con_ref).blocks
166-
println(gram.basis.polynomials)
167-
display(gram.Q)
168-
end
174+
polys = gram[3].basis.bases[].elements #src
175+
@test length(polys) == 1 #src
176+
@test polys[] ≈ (basis[1] - basis[2] * im) / √2 #src
177+
gram_matrix(con_ref)
169178

170179
# We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:
171180

‎docs/src/tutorials/Symmetry/dihedral.jl‎

Lines changed: 40 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,6 @@
88
# *Symmetry groups, semidefinite programs, and sums of squares*.
99
# Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128.
1010

11-
1211
using Test #src
1312
import Random #src
1413
Random.seed!(0) #src
@@ -134,9 +133,9 @@ end
134133

135134
# We can exploit this symmetry for reducing the problem using the `SymmetricIdeal` certificate as follows:
136135

137-
import CSDP
136+
import Clarabel
138137
function solve(G)
139-
solver = CSDP.Optimizer
138+
solver = Clarabel.Optimizer
140139
model = Model(solver)
141140
@variable(model, t)
142141
@objective(model, Max, t)
@@ -148,37 +147,44 @@ function solve(G)
148147

149148

150149
g = gram_matrix(con_ref).blocks #src
151-
@test length(g) == 5 #src
152-
@test g[1].basis.polynomials == [y^3, x^2*y, y] #src
153-
@test g[2].basis.polynomials == [-x^3, -x*y^2, -x] #src
154-
for i in 1:2 #src
155-
I = 3:-1:1 #src
156-
Q = g[i].Q[I, I] #src
157-
@test size(Q) == (3, 3) #src
158-
@test Q[2, 2] ≈ 1 rtol=1e-2 #src
159-
@test Q[1, 2] ≈ 5/8 rtol=1e-2 #src
160-
@test Q[2, 3] ≈ -1 rtol=1e-2 #src
161-
@test Q[1, 1] ≈ 25/64 rtol=1e-2 #src
162-
@test Q[1, 3] ≈ -5/8 rtol=1e-2 #src
163-
@test Q[3, 3] ≈ 1 rtol=1e-2 #src
164-
end #src
165-
@test length(g[3].basis.polynomials) == 2 #src
166-
@test g[3].basis.polynomials[1] == 1.0 #src
167-
@test g[3].basis.polynomials[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
168-
@test size(g[3].Q) == (2, 2) #src
169-
@test g[3].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
170-
@test g[3].Q[1, 2] ≈ 0.983 rtol=1e-2 #src
171-
@test g[3].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
172-
@test g[4].basis.polynomials == [x * y] #src
173-
@test size(g[4].Q) == (1, 1) #src
174-
@test g[4].Q[1, 1] ≈ 0 atol=1e-2 #src
175-
@test length(g[5].basis.polynomials) == 1 #src
176-
@test g[5].basis.polynomials[1] ≈ (√2/2)x^2 - (√2/2)y^2 #src
177-
@test size(g[5].Q) == (1, 1) #src
178-
@test g[5].Q[1, 1] ≈ 0 atol=1e-2 #src
179-
for g in gram_matrix(con_ref).blocks
180-
println(g.basis.polynomials)
181-
end
150+
@test length(g) == 4 #src
151+
@test length(g[4].basis.bases) == 2 #src
152+
polys = g[4].basis.bases[1].elements #src
153+
@test length(polys) == 3 #src
154+
@test polys[1] ≈ y^3 #src
155+
@test polys[2] ≈ x^2*y #src
156+
@test polys[3] ≈ y #src
157+
polys = g[4].basis.bases[2].elements #src
158+
@test length(polys) == 3 #src
159+
@test polys[1] ≈ -x^3 #src
160+
@test polys[2] ≈ -x*y^2 #src
161+
@test polys[3] ≈ -x #src
162+
I = 3:-1:1 #src
163+
Q = g[4].Q[I, I] #src
164+
@test size(Q) == (3, 3) #src
165+
@test Q[2, 2] ≈ 1 rtol=1e-2 #src
166+
@test Q[1, 2] ≈ 5/8 rtol=1e-2 #src
167+
@test Q[2, 3] ≈ -1 rtol=1e-2 #src
168+
@test Q[1, 1] ≈ 25/64 rtol=1e-2 #src
169+
@test Q[1, 3] ≈ -5/8 rtol=1e-2 #src
170+
@test Q[3, 3] ≈ 1 rtol=1e-2 #src
171+
polys = g[1].basis.bases[].elements #src
172+
@test length(polys) == 2 #src
173+
@test polys[1] ≈ 1.0 #src
174+
@test polys[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
175+
@test size(g[1].Q) == (2, 2) #src
176+
@test g[1].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
177+
@test g[1].Q[1, 2] ≈ 0.983 rtol=1e-2 #src
178+
@test g[1].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
179+
polys = g[2].basis.bases[].elements #src
180+
@test polys[] ≈ x * y #src
181+
@test size(g[2].Q) == (1, 1) #src
182+
@test g[2].Q[1, 1] ≈ 0 atol=1e-2 #src
183+
polys = g[3].basis.bases[].elements #src
184+
@test polys[] ≈ (√2/2)x^2 - (√2/2)y^2 #src
185+
@test size(g[3].Q) == (1, 1) #src
186+
@test g[3].Q[1, 1] ≈ 0 atol=1e-2 #src
187+
gram_matrix(con_ref)
182188
end
183189
solve(G)
184190

‎docs/src/tutorials/Symmetry/even_reduction.jl‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -31,8 +31,8 @@ G = PermGroup([perm"(1,2)"])
3131

3232
# We can exploit the symmetry as follows:
3333

34-
import CSDP
35-
solver = CSDP.Optimizer
34+
import Clarabel
35+
solver = Clarabel.Optimizer
3636
model = Model(solver)
3737
@variable(model, t)
3838
@objective(model, Max, t)
@@ -45,6 +45,9 @@ value(t)
4545
# We indeed find `-1`, let's verify that symmetry was exploited:
4646

4747
@test length(gram_matrix(con_ref).blocks) == 2 #src
48-
@test gram_matrix(con_ref).blocks[1].basis.polynomials == [1, x^2] #src
49-
@test gram_matrix(con_ref).blocks[2].basis.polynomials == [x] #src
48+
polys = gram_matrix(con_ref).blocks[1].basis.bases[].elements #src
49+
@test polys[1] ≈ 1 #src
50+
@test polys[2] ≈ x^2 #src
51+
polys = gram_matrix(con_ref).blocks[2].basis.bases[].elements #src
52+
@test polys[] ≈ x #src
5053
gram_matrix(con_ref)

‎docs/src/tutorials/Symmetry/permutation_symmetry.jl‎

Lines changed: 25 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -32,8 +32,8 @@ G = PermGroup([perm"(1,2,3,4)"])
3232

3333
# We can use this certificate as follows:
3434

35-
import CSDP
36-
solver = CSDP.Optimizer
35+
import Clarabel
36+
solver = Clarabel.Optimizer
3737
model = Model(solver)
3838
@variable(model, t)
3939
@objective(model, Max, t)
@@ -45,25 +45,28 @@ value(t)
4545

4646
# We indeed find `-1`, let's verify that symmetry was exploited:
4747

48-
g = gram_matrix(con_ref).blocks #src
49-
@test length(g) == 4 #src
50-
@test length(g[1].basis.polynomials) == 2 #src
51-
@test g[1].basis.polynomials[1] == 1.0 #src
52-
@test g[1].basis.polynomials[2] ≈ -0.5 * sum(x) #src
53-
@test size(g[1].Q) == (2, 2) #src
54-
@test g[1].Q[1, 1] ≈ 1.0 atol=1e-6 #src
55-
@test g[1].Q[1, 2] ≈ -1.0 atol=1e-6 #src
56-
@test g[1].Q[2, 2] ≈ 1.0 atol=1e-6 #src
57-
@test length(g[2].basis.polynomials) == 1 #src
58-
@test g[2].basis.polynomials[1] ≈ (x[2] - x[4]) / √2 #src
48+
gram = gram_matrix(con_ref).blocks #src
49+
@test length(gram) == 3 #src
50+
polys = gram[1].basis.bases[].elements #src
51+
@test length(polys) == 2 #src
52+
@test polys[1] ≈ 1 #src
53+
@test polys[2] ≈ -0.5 * sum(x) #src
54+
@test size(gram[1].Q) == (2, 2) #src
55+
@test gram[1].Q[1, 1] ≈ 1.0 atol=1e-6 #src
56+
@test gram[1].Q[1, 2] ≈ -1.0 atol=1e-6 #src
57+
@test gram[1].Q[2, 2] ≈ 1.0 atol=1e-6 #src
58+
@test length(gram[2].basis.bases) == 2 #src
59+
polys = gram[2].basis.bases[1].elements #src
60+
@test length(polys) == 1 #src
61+
@test polys[] ≈ (x[2] - x[4]) / √2 #src
5962
@test size(g[2].Q) == (1, 1) #src
6063
@test g[2].Q[1, 1] ≈ 1.0 atol=1e-6 #src
61-
@test length(g[3].basis.polynomials) == 1 #src
62-
@test g[3].basis.polynomials[1] ≈ (x[1] - x[3]) / √2 #src
63-
@test size(g[3].Q) == (1, 1) #src
64-
@test g[3].Q[1, 1] ≈ 1.0 atol=1e-6 #src
65-
@test length(g[4].basis.polynomials) == 1 #src
66-
@test g[4].basis.polynomials[1] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
67-
@test size(g[4].Q) == (1, 1) #src
68-
@test g[4].Q[1, 1] ≈ 1.0 atol=1e-6 #src
69-
gram_matrix(con_ref).blocks
64+
polys = gram[2].basis.bases[2].elements #src
65+
@test length(polys) == 1 #src
66+
@test polys[1] ≈ (x[1] - x[3]) / √2 #src
67+
polys = gram[3].basis.bases[].elements #src
68+
@test length(polys) == 1 #src
69+
@test polys[] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
70+
@test size(gram[3].Q) == (1, 1) #src
71+
@test gram[3].Q[1, 1] ≈ 1.0 atol=1e-6 #src
72+
gram_matrix(con_ref)

‎src/Certificate/Symmetry/Symmetry.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module Symmetry
33
import LinearAlgebra
44

55
import MutableArithmetics as MA
6+
import StarAlgebras as SA
67
import MultivariatePolynomials as MP
78
import MultivariateBases as MB
89

0 commit comments

Comments
 (0)