From a1a9b9239668a3b2568a7ec87906e99333933fe8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 20 Aug 2024 13:30:56 +0200 Subject: [PATCH] [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 | 16 +---- docs/src/variables.md | 8 +-- src/Certificate/Symmetry/block_diag.jl | 56 +++++++++++------ test/symmetry.jl | 83 ++++++++++++++++++-------- 4 files changed, 99 insertions(+), 64 deletions(-) diff --git a/docs/src/constraints.md b/docs/src/constraints.md index 354bd1d0e..7e1b87bbe 100644 --- a/docs/src/constraints.md +++ b/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, α) α diff --git a/docs/src/variables.md b/docs/src/variables.md index a74da1032..8ff06a842 100644 --- a/docs/src/variables.md +++ b/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² diff --git a/src/Certificate/Symmetry/block_diag.jl b/src/Certificate/Symmetry/block_diag.jl index 091b9477c..e48663a94 100644 --- a/src/Certificate/Symmetry/block_diag.jl +++ b/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,10 +160,13 @@ 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'` @@ -162,7 +174,7 @@ Given Schur decompositions 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 diff --git a/test/symmetry.jl b/test/symmetry.jl index c267d8524..2a66fee06 100644 --- a/test/symmetry.jl +++ b/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()