[Symmetry] Fix sign scaling in ordered_block_diag (#384)

* Bugs in ordered_block_diag

* Fix

* Remove test

* Fix format

* Fix doc

docs/src/constraints.md

@@ -32,13 +32,7 @@ julia> X = monomials(x, 0:2)
 
 julia> using SumOfSquares
 
-julia> model = Model()
-A JuMP Model
-Feasibility problem with:
-Variables: 0
-Model mode: AUTOMATIC
-CachingOptimizer state: NO_OPTIMIZER
-Solver name: No optimizer attached.
+julia> model = Model();
 
 julia> @variable(model, p, Poly(X))
 (_[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
 
 julia> using SumOfSquares
 
-julia> model = SOSModel()
-A JuMP Model
-Feasibility problem with:
-Variables: 0
-Model mode: AUTOMATIC
-CachingOptimizer state: NO_OPTIMIZER
-Solver name: No optimizer attached.
+julia> model = SOSModel();
 
 julia> @variable(model, α)
 α

docs/src/variables.md

@@ -33,13 +33,7 @@ We can now create our polynomial variable `p` as follows:
 ```jldoctest variables
 julia> using SumOfSquares
 
-julia> model = Model()
-A JuMP Model
-Feasibility problem with:
-Variables: 0
-Model mode: AUTOMATIC
-CachingOptimizer state: NO_OPTIMIZER
-Solver name: No optimizer attached.
+julia> model = Model();
 
 julia> @variable(model, p, Poly(X))
 (_[1]) + (_[2])y + (_[3])x + (_[4])y² + (_[5])xy + (_[6])x²

src/Certificate/Symmetry/block_diag.jl

@@ -53,39 +53,48 @@ end
 # We can multiply by `Diagonal(d)` if `d[i] * conj(d[i]) = 1`.
 # So in the real case, `d = ±1` but in the complex case, we have more freedom.
 function _sign_diag(
-    A::AbstractMatrix{T},
-    B::AbstractMatrix{T};
+    As::Vector{<:AbstractMatrix{T}},
+    Bs::Vector{<:AbstractMatrix{T}};
     tol = Base.rtoldefault(real(T)),
 ) where {T}
-    n = LinearAlgebra.checksquare(A)
+    n = LinearAlgebra.checksquare(As[1])
     d = ones(T, n)
     for j in 2:n
         if T <: Real
             minus = zero(real(T))
             not_minus = zero(real(T))
             for i in 1:(j-1)
-                a = A[i, j]
-                b = B[i, j]
-                minus = max(minus, abs(a + b))
-                not_minus = max(not_minus, abs(a - b))
+                for (Ai, Bi) in zip(As, Bs)
+                    a = Ai[i, j]
+                    b = Bi[i, j]
+                    minus = max(minus, abs(a + b))
+                    not_minus = max(not_minus, abs(a - b))
+                end
             end
             if minus < not_minus
                 d[j] = -one(T)
-                B[:, j] = -B[:, j]
-                B[j, :] = -B[j, :]
+                for B in Bs
+                    B[:, j] = -B[:, j]
+                    B[j, :] = -B[j, :]
+                end
             end
         else
-            i = argmax(abs.(B[1:(j-1), j]))
-            if abs(B[i, j]) <= tol
+            k = argmax(eachindex(Bs)) do k
+                return maximum(abs.(Bs[k][1:(j-1), j]))
+            end
+            i = argmax(abs.(Bs[k][1:(j-1), j]))
+            if abs(Bs[k][i, j]) <= tol
                 continue
             end
-            rot = A[i, j] / B[i, j]
+            rot = As[k][i, j] / Bs[k][i, j]
             # It should be unitary but there might be small numerical errors
             # so let's normalize
             rot /= abs(rot)
             d[j] = rot
-            B[:, j] *= rot
-            B[j, :] *= conj(rot)
+            for B in Bs
+                B[:, j] *= rot
+                B[j, :] *= conj(rot)
+            end
         end
     end
     return d
@@ -151,18 +160,21 @@ function _rotate_complex(
 end
 
 """
-    orthogonal_transformation_to(A, B)
+    orthogonal_transformation_to(A, B, Ais, Bis)
 
 Return an orthogonal transformation `U` such that
 `A = U' * B * U`
+and
+`Ai[k] = U' * Bi[k] * U`
+for all `(Ai, Bi) in zip(Ais, Bis)`
 
 Given Schur decompositions
 `A = Z_A * S_A * Z_A'`
 `B = Z_B * S_B * Z_B'`
 Since `P' * S_A * P = D' * S_B * D`, we have
 `A = Z_A * P * Z_B' * B * Z_B * P' * Z_A'`
 """
-function orthogonal_transformation_to(A, B)
+function orthogonal_transformation_to(A, B, Ais, Bis)
     As = LinearAlgebra.schur(A)
     _reorder!(As)
     T_A = As.Schur
@@ -173,7 +185,15 @@ function orthogonal_transformation_to(A, B)
     Z_B = Bs.vectors
     P = _rotate_complex(T_A, T_B)
     T_A = P' * T_A * P
-    d = _sign_diag(T_A, T_B)
+    # If `T_A` and `T_B` are diagonal, they don't help
+    # determing the diagonal `d`.
+    # We provide the other matrices of `Ais` and `Bis`
+    # in order to help determine the sign in these cases.
+    Ais = [P' * Z_A' * A * Z_A * P for A in Ais]
+    push!(Ais, T_A)
+    Bis = [Z_B' * B * Z_B for B in Bis]
+    push!(Bis, T_B)
+    d = _sign_diag(Ais, Bis)
     D = LinearAlgebra.Diagonal(d)
     return _try_integer!(Z_B * D * P' * Z_A')
 end
@@ -204,7 +224,7 @@ function ordered_block_diag(As, d)
         # 1997, 133-140
         R = sum(λ .* refs)
         C = sum(λ .* Cs)
-        V = orthogonal_transformation_to(R, C)
+        V = orthogonal_transformation_to(R, C, refs, Cs)
         @assert R ≈ V' * C * V
         for i in eachindex(refs)
             @assert refs[i] ≈ V' * Cs[i] * V

test/symmetry.jl

@@ -39,11 +39,29 @@ function test_linsolve()
     end
 end
 
-function _test_orthogonal_transformation_to(A, B)
-    U = SumOfSquares.Certificate.Symmetry.orthogonal_transformation_to(A, B)
+function _test_orthogonal_transformation_to(A, B, As, Bs)
+    U = SumOfSquares.Certificate.Symmetry.orthogonal_transformation_to(
+        A,
+        B,
+        As,
+        Bs,
+    )
     @test A ≈ U' * B * U
 end
 
+function _test_orthogonal_transformation_to(λ, As, Bs)
+    A = sum(λ[i] * As[i] for i in eachindex(As))
+    B = sum(λ[i] * Bs[i] for i in eachindex(Bs))
+    _test_orthogonal_transformation_to(A, B, As, Bs)
+    return
+end
+
+function _test_orthogonal_transformation_to(A, B)
+    _test_orthogonal_transformation_to(A, B, typeof(A)[], typeof(B)[])
+    _test_orthogonal_transformation_to(A, B, [A], [B])
+    return
+end
+
 function _test_orthogonal_transformation_to(T::Type)
     A1 = T[
         0 -1
@@ -152,6 +170,10 @@ function _test_orthogonal_transformation_to(T::Type)
         1 1 0
     ]
     _test_orthogonal_transformation_to(A1, A2)
+    A1 = T[1 0; 0 -1]
+    A2 = T[0 1; 1 0]
+    _test_orthogonal_transformation_to([1, 1], [A1, A2], [-A1, A2])
+    _test_orthogonal_transformation_to([2, 1], [A1, A2], [-A1, A2])
     return
 end
 
@@ -161,29 +183,40 @@ function test_orthogonal_transformation_to()
     end
 end
 
-function test_block_diag()
-    # From `dihedral.jl` example
-    A = [
-        [
-            0 -1 0 0 0 0
-            1 0 0 0 0 0
-            0 0 0 0 0 -1
-            0 0 0 0 1 0
-            0 0 0 -1 0 0
-            0 0 1 0 0 0
-        ],
-        [
-            0 1 0 0 0 0
-            1 0 0 0 0 0
-            0 0 0 0 0 1
-            0 0 0 0 1 0
-            0 0 0 1 0 0
-            0 0 1 0 0 0
-        ],
-    ]
-    d = 2
+function _test_block_diag(A, d)
     U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
     @test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
+    return
+end
+
+function test_block_diag_dihedral()
+    # From `dihedral.jl` example
+    d = 2
+    A2 = [
+        0 1 0 0 0 0
+        1 0 0 0 0 0
+        0 0 0 0 0 1
+        0 0 0 0 1 0
+        0 0 0 1 0 0
+        0 0 1 0 0 0
+    ]
+    A1 = [
+        0 -1 0 0 0 0
+        1 0 0 0 0 0
+        0 0 0 0 0 -1
+        0 0 0 0 1 0
+        0 0 0 -1 0 0
+        0 0 1 0 0 0
+    ]
+    _test_block_diag([A1, A2], d)
+    # Using `GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, true, 1), DihedralElement(G.n, true, 0)]`, we get
+    # see https://github.com/jump-dev/SumOfSquares.jl/issues/381#issuecomment-2296711306
+    A1 = Matrix(Diagonal([1, -1, 1, -1, 1, -1]))
+    _test_block_diag([A1, A2], d)
+    return
+end
+
+function test_block_diag_alpha()
     α = 0.75
     A = [
         Matrix{Int}(I, 6, 6),
@@ -205,8 +238,8 @@ function test_block_diag()
         ],
     ]
     d = 2
-    U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
-    @test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
+    _test_block_diag(A, d)
+    return
 end
 
 function runtests()