Fix edge case for block diagonalization

src/Certificate/Symmetry/block_diag.jl

@@ -180,7 +180,7 @@
 
 function ordered_block_diag(As, d)
     U = block_diag(As, d)
-    U === nothing && return nothing
+    isnothing(U) && return nothing
     iU = U'
     @assert iU ≈ inv(U)
     Bs = [iU * A * U for A in As]
@@ -199,7 +199,8 @@
         # should work, this trick is similar to [CGT97].
         #
         # [CGT97] Corless, R. M.; Gianni, P. M. & Trager, B. M.
-        # A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots Proceedings of the 1997 international symposium on Symbolic and algebraic computation,
+        # A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots
+        # Proceedings of the 1997 international symposium on Symbolic and algebraic computation,
         # 1997, 133-140
         R = sum(λ .* refs)
         C = sum(λ .* Cs)
@@ -247,13 +248,14 @@
             offset = 0
             for v in values(blocks)
                 @assert iszero(length(v) % d)
-                if length(v) == d
-                    V = Z[:, v]
-                else
+                V = Z[:, v]
+                if length(v) != d
                     Cs = [B[v, v] for B in Bs]
-                    V = block_diag(Cs, d)
-                    V === nothing && break
-                    V *= transpose(Z[:, v])
+                    _V = block_diag(Cs, d)
+                    if isnothing(_V)
+                        break
+                    end
+                    V *= _V
                 end
                 U[:, offset.+eachindex(v)] = V
                 offset += length(v)

src/Certificate/Symmetry/wedderburn.jl

@@ -185,14 +185,14 @@
         if d > 1
             if m > 1
                 U = ordered_block_diag(S, d)
+                if isnothing(U)
+                    error(
+                        "Could not simultaneously block-diagonalize into $m identical $dx$d blocks",
+                    )
+                end
             else
                 U = Matrix{T}(LinearAlgebra.I, N, N)
             end
-            if U === nothing
-                error(
-                    "Could not simultaneously block-diagonalize into $m identical $dx$d blocks",
-                )
-            end
             # From Example 1.7.3 of
             # Sagan, The symmetric group, Springer Science & Business Media, 2001
             # we know that there exists `C` such that `Q = kron(C, I)` if we use

test/symmetry.jl

@@ -184,6 +184,29 @@ 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)
+    α = 0.75
+    A = [
+        Matrix{Int}(I, 6, 6),
+        [
+            1 0 0 0 0 0
+            0 -1 0 0 0 0
+            0 0 0 0 1 0
+            0 0 0 0 0 -1
+            0 0 1 0 0 0
+            0 0 0 -1 0 0
+        ],
+        [
+            -0.5 α 0 0 0 0
+            -α -0.5 0 0 0 0
+            0 0 -0.5 α 0 0
+            0 0 -α -0.5 0 0
+            0 0 0 0 -0.5 α
+            0 0 0 0 -α -0.5
+        ],
+    ]
+    d = 2
+    U = SumOfSquares.Certificate.Symmetry.ordered_block_diag(A, d)
+    @test SumOfSquares.Certificate.Symmetry.ordered_block_check(U, A, d)
 end
 
 function runtests()