Skip to content

Commit d28eb76

Browse files
authored
[Symmetry] Fix sign scaling in ordered_block_diag (#384)
* Bugs in ordered_block_diag * Fix * Remove test * Fix format * Fix doc
1 parent 4af4b3e commit d28eb76

File tree

4 files changed

+99
-64
lines changed

4 files changed

+99
-64
lines changed

docs/src/constraints.md

Lines changed: 2 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -32,13 +32,7 @@ julia> X = monomials(x, 0:2)
3232
3333
julia> using SumOfSquares
3434
35-
julia> model = Model()
36-
A JuMP Model
37-
Feasibility problem with:
38-
Variables: 0
39-
Model mode: AUTOMATIC
40-
CachingOptimizer state: NO_OPTIMIZER
41-
Solver name: No optimizer attached.
35+
julia> model = Model();
4236
4337
julia> @variable(model, p, Poly(X))
4438
(_[1])·1 + (_[2])·x₃ + (_[3])·x₂ + (_[4])·x₁ + (_[5])·x₃² + (_[6])·x₂x₃ + (_[7])·x₂² + (_[8])·x₁x₃ + (_[9])·x₁x₂ + (_[10])·x₁²
@@ -132,13 +126,7 @@ julia> @polyvar x y
132126
133127
julia> using SumOfSquares
134128
135-
julia> model = SOSModel()
136-
A JuMP Model
137-
Feasibility problem with:
138-
Variables: 0
139-
Model mode: AUTOMATIC
140-
CachingOptimizer state: NO_OPTIMIZER
141-
Solver name: No optimizer attached.
129+
julia> model = SOSModel();
142130
143131
julia> @variable(model, α)
144132
α

docs/src/variables.md

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -33,13 +33,7 @@ We can now create our polynomial variable `p` as follows:
3333
```jldoctest variables
3434
julia> using SumOfSquares
3535
36-
julia> model = Model()
37-
A JuMP Model
38-
Feasibility problem with:
39-
Variables: 0
40-
Model mode: AUTOMATIC
41-
CachingOptimizer state: NO_OPTIMIZER
42-
Solver name: No optimizer attached.
36+
julia> model = Model();
4337
4438
julia> @variable(model, p, Poly(X))
4539
(_[1])·1 + (_[2])·y + (_[3])·x + (_[4])·y² + (_[5])·xy + (_[6])·x²

src/Certificate/Symmetry/block_diag.jl

Lines changed: 38 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -53,39 +53,48 @@ end
5353
# We can multiply by `Diagonal(d)` if `d[i] * conj(d[i]) = 1`.
5454
# So in the real case, `d = ±1` but in the complex case, we have more freedom.
5555
function _sign_diag(
56-
A::AbstractMatrix{T},
57-
B::AbstractMatrix{T};
56+
As::Vector{<:AbstractMatrix{T}},
57+
Bs::Vector{<:AbstractMatrix{T}};
5858
tol = Base.rtoldefault(real(T)),
5959
) where {T}
60-
n = LinearAlgebra.checksquare(A)
60+
n = LinearAlgebra.checksquare(As[1])
6161
d = ones(T, n)
6262
for j in 2:n
6363
if T <: Real
6464
minus = zero(real(T))
6565
not_minus = zero(real(T))
6666
for i in 1:(j-1)
67-
a = A[i, j]
68-
b = B[i, j]
69-
minus = max(minus, abs(a + b))
70-
not_minus = max(not_minus, abs(a - b))
67+
for (Ai, Bi) in zip(As, Bs)
68+
a = Ai[i, j]
69+
b = Bi[i, j]
70+
minus = max(minus, abs(a + b))
71+
not_minus = max(not_minus, abs(a - b))
72+
end
7173
end
7274
if minus < not_minus
7375
d[j] = -one(T)
74-
B[:, j] = -B[:, j]
75-
B[j, :] = -B[j, :]
76+
for B in Bs
77+
B[:, j] = -B[:, j]
78+
B[j, :] = -B[j, :]
79+
end
7680
end
7781
else
78-
i = argmax(abs.(B[1:(j-1), j]))
79-
if abs(B[i, j]) <= tol
82+
k = argmax(eachindex(Bs)) do k
83+
return maximum(abs.(Bs[k][1:(j-1), j]))
84+
end
85+
i = argmax(abs.(Bs[k][1:(j-1), j]))
86+
if abs(Bs[k][i, j]) <= tol
8087
continue
8188
end
82-
rot = A[i, j] / B[i, j]
89+
rot = As[k][i, j] / Bs[k][i, j]
8390
# It should be unitary but there might be small numerical errors
8491
# so let's normalize
8592
rot /= abs(rot)
8693
d[j] = rot
87-
B[:, j] *= rot
88-
B[j, :] *= conj(rot)
94+
for B in Bs
95+
B[:, j] *= rot
96+
B[j, :] *= conj(rot)
97+
end
8998
end
9099
end
91100
return d
@@ -151,18 +160,21 @@ function _rotate_complex(
151160
end
152161

153162
"""
154-
orthogonal_transformation_to(A, B)
163+
orthogonal_transformation_to(A, B, Ais, Bis)
155164
156165
Return an orthogonal transformation `U` such that
157166
`A = U' * B * U`
167+
and
168+
`Ai[k] = U' * Bi[k] * U`
169+
for all `(Ai, Bi) in zip(Ais, Bis)`
158170
159171
Given Schur decompositions
160172
`A = Z_A * S_A * Z_A'`
161173
`B = Z_B * S_B * Z_B'`
162174
Since `P' * S_A * P = D' * S_B * D`, we have
163175
`A = Z_A * P * Z_B' * B * Z_B * P' * Z_A'`
164176
"""
165-
function orthogonal_transformation_to(A, B)
177+
function orthogonal_transformation_to(A, B, Ais, Bis)
166178
As = LinearAlgebra.schur(A)
167179
_reorder!(As)
168180
T_A = As.Schur
@@ -173,7 +185,15 @@ function orthogonal_transformation_to(A, B)
173185
Z_B = Bs.vectors
174186
P = _rotate_complex(T_A, T_B)
175187
T_A = P' * T_A * P
176-
d = _sign_diag(T_A, T_B)
188+
# If `T_A` and `T_B` are diagonal, they don't help
189+
# determing the diagonal `d`.
190+
# We provide the other matrices of `Ais` and `Bis`
191+
# in order to help determine the sign in these cases.
192+
Ais = [P' * Z_A' * A * Z_A * P for A in Ais]
193+
push!(Ais, T_A)
194+
Bis = [Z_B' * B * Z_B for B in Bis]
195+
push!(Bis, T_B)
196+
d = _sign_diag(Ais, Bis)
177197
D = LinearAlgebra.Diagonal(d)
178198
return _try_integer!(Z_B * D * P' * Z_A')
179199
end
@@ -204,7 +224,7 @@ function ordered_block_diag(As, d)
204224
# 1997, 133-140
205225
R = sum.* refs)
206226
C = sum.* Cs)
207-
V = orthogonal_transformation_to(R, C)
227+
V = orthogonal_transformation_to(R, C, refs, Cs)
208228
@assert R V' * C * V
209229
for i in eachindex(refs)
210230
@assert refs[i] V' * Cs[i] * V

test/symmetry.jl

Lines changed: 58 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -39,11 +39,29 @@ function test_linsolve()
3939
end
4040
end
4141

42-
function _test_orthogonal_transformation_to(A, B)
43-
U = SumOfSquares.Certificate.Symmetry.orthogonal_transformation_to(A, B)
42+
function _test_orthogonal_transformation_to(A, B, As, Bs)
43+
U = SumOfSquares.Certificate.Symmetry.orthogonal_transformation_to(
44+
A,
45+
B,
46+
As,
47+
Bs,
48+
)
4449
@test A U' * B * U
4550
end
4651

52+
function _test_orthogonal_transformation_to(λ, As, Bs)
53+
A = sum(λ[i] * As[i] for i in eachindex(As))
54+
B = sum(λ[i] * Bs[i] for i in eachindex(Bs))
55+
_test_orthogonal_transformation_to(A, B, As, Bs)
56+
return
57+
end
58+
59+
function _test_orthogonal_transformation_to(A, B)
60+
_test_orthogonal_transformation_to(A, B, typeof(A)[], typeof(B)[])
61+
_test_orthogonal_transformation_to(A, B, [A], [B])
62+
return
63+
end
64+
4765
function _test_orthogonal_transformation_to(T::Type)
4866
A1 = T[
4967
0 -1
@@ -152,6 +170,10 @@ function _test_orthogonal_transformation_to(T::Type)
152170
1 1 0
153171
]
154172
_test_orthogonal_transformation_to(A1, A2)
173+
A1 = T[1 0; 0 -1]
174+
A2 = T[0 1; 1 0]
175+
_test_orthogonal_transformation_to([1, 1], [A1, A2], [-A1, A2])
176+
_test_orthogonal_transformation_to([2, 1], [A1, A2], [-A1, A2])
155177
return
156178
end
157179

@@ -161,29 +183,40 @@ function test_orthogonal_transformation_to()
161183
end
162184
end
163185

164-
function test_block_diag()
165-
# From `dihedral.jl` example
166-
A = [
167-
[
168-
0 -1 0 0 0 0
169-
1 0 0 0 0 0
170-
0 0 0 0 0 -1
171-
0 0 0 0 1 0
172-
0 0 0 -1 0 0
173-
0 0 1 0 0 0
174-
],
175-
[
176-
0 1 0 0 0 0
177-
1 0 0 0 0 0
178-
0 0 0 0 0 1
179-
0 0 0 0 1 0
180-
0 0 0 1 0 0
181-
0 0 1 0 0 0
182-
],
183-
]
184-
d = 2
186+
function _test_block_diag(A, d)
185187
U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
186188
@test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
189+
return
190+
end
191+
192+
function test_block_diag_dihedral()
193+
# From `dihedral.jl` example
194+
d = 2
195+
A2 = [
196+
0 1 0 0 0 0
197+
1 0 0 0 0 0
198+
0 0 0 0 0 1
199+
0 0 0 0 1 0
200+
0 0 0 1 0 0
201+
0 0 1 0 0 0
202+
]
203+
A1 = [
204+
0 -1 0 0 0 0
205+
1 0 0 0 0 0
206+
0 0 0 0 0 -1
207+
0 0 0 0 1 0
208+
0 0 0 -1 0 0
209+
0 0 1 0 0 0
210+
]
211+
_test_block_diag([A1, A2], d)
212+
# Using `GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, true, 1), DihedralElement(G.n, true, 0)]`, we get
213+
# see https://github.com/jump-dev/SumOfSquares.jl/issues/381#issuecomment-2296711306
214+
A1 = Matrix(Diagonal([1, -1, 1, -1, 1, -1]))
215+
_test_block_diag([A1, A2], d)
216+
return
217+
end
218+
219+
function test_block_diag_alpha()
187220
α = 0.75
188221
A = [
189222
Matrix{Int}(I, 6, 6),
@@ -205,8 +238,8 @@ function test_block_diag()
205238
],
206239
]
207240
d = 2
208-
U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
209-
@test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
241+
_test_block_diag(A, d)
242+
return
210243
end
211244

212245
function runtests()

0 commit comments

Comments
 (0)