Add deterministic property-based regression suite

lib/centralizer.gi

@@ -96,7 +96,19 @@ end );
 
 # Same as DecomposeCentralizerBlocks but converts to full matrices
 InstallGlobalFunction( CentralizerOfRepresentation, function(rho)
-    return List(CentralizerBlocksOfRepresentation(rho), BlockDiagonalMatrix@);
+    local decomposition, basis_change, block_basis;
+
+    decomposition := ComputeUsingMethod@(rho);
+    basis_change := TransposedMat(decomposition.basis);
+    block_basis := List(decomposition.centralizer_basis,
+                        BlockDiagonalMatrix@);
+
+    # The standard block generators commute with the block-diagonalized
+    # representation.  Conjugate them back so that this function returns the
+    # centralizer in the coordinates of the representation supplied by the
+    # caller.
+    return List(block_basis,
+                C -> basis_change * C * basis_change^-1);
 end );
 
 # Given an orthonormal basis for the centralizer (w.r.t the product

lib/groupsum.gi

@@ -1,6 +1,10 @@
 InstallGlobalFunction( GroupSumBSGS, function(G, summand)
     local H, n, chain, groups, m, sum, i, cosets, iso, zero, g, S, right_reps;
 
+    if IsTrivial(G) then
+        return summand(One(G));
+    fi;
+
     # stab chain computation only works on permutation groups
     iso := IsomorphismPermGroup(G);
     H := Image(iso, G);

lib/serre.gi

@@ -155,6 +155,14 @@ InstallGlobalFunction( CanonicalDecomposition, function(rho)
     # The group we are taking representations of
     G := Source(rho);
 
+    # All projector calculations below act on vectors by matrices.  Preserve
+    # the public convenience of accepting permutation representations by
+    # converting once and then consistently using the linear representation.
+    rho := ConvertRhoIfNeeded@(rho);
+    if rho = fail then
+        Error("<rho> must be a finite linear or permutation representation");
+    fi;
+
     # The list of irreps W_i of G over F appearing in rho
     irreps := ValueOption("irreps");
 
@@ -164,8 +172,7 @@ InstallGlobalFunction( CanonicalDecomposition, function(rho)
 
     # We might need to convert here, since this function needs a
     # linear rep
-    irreps := RelevantIrreps@(ConvertRhoIfNeeded@(rho),
-                              irreps);
+    irreps := RelevantIrreps@(rho, irreps);
 
     # if there's only 1 irrep, the canonical summand is just the whole
     # space!

lib/utils.gi

@@ -268,13 +268,14 @@ OrthonormalBasis@ := function(v)
 end;
 
 InstallGlobalFunction( IsOrthonormalSet, function(S, prod)
-    return ForAll(S, v1 -> ForAll(S, function(v2)
-                                     if v1 = v2 then
-                                         return prod(v1, v2) = 1;
-                                     else
-                                         return prod(v1, v2) = 0;
-                                     fi;
-                                 end ));
+    return ForAll([1..Length(S)],
+                  i -> ForAll([1..Length(S)], function(j)
+                      if i = j then
+                          return prod(S[i], S[j]) = 1;
+                      else
+                          return prod(S[i], S[j]) = 0;
+                      fi;
+                  end ));
 end );
 
 KroneckerList@ := function(reps)

tst/README.md

@@ -0,0 +1,39 @@
+# RepnDecomp test strategy
+
+The hand-written tests in this directory complement the examples generated by
+AutoDoc.  They use fixed groups, matrices, and random seeds, and primarily test
+mathematical invariants instead of a particular choice of basis.
+
+The suite checks that:
+
+- returned summands are invariant, irreducible, independent, and exhaustive;
+- collected summands agree exactly with representation isomorphism classes;
+- character inner products recover the expected multiplicities;
+- a returned basis simultaneously block diagonalizes every group generator;
+- the centralizer has dimension `sum(m_i^2)`, is an algebra containing the
+  identity, and commutes with the original representation in its original
+  coordinates;
+- centralizer projection reproduces direct conjugacy-class sums;
+- explicit intertwiners satisfy their defining equations;
+- BSGS group sums equal direct enumeration;
+- tensor, unitarization, LDL, permutation, and orbital operations satisfy their
+  defining identities.
+
+These checks reflect the principal computational use cases and certification
+criteria discussed in:
+
+- Hymabaccus and Pasechnik, *Decomposing Linear Representations of Finite
+  Groups*, <https://arxiv.org/abs/2007.02459>;
+- Vallentin, *Symmetry in semidefinite programs*,
+  <https://arxiv.org/abs/0706.4233>;
+- Rosset, Montealegre-Mora, and Bancal, *RepLAB: a
+  computational/numerical approach to representation theory*,
+  <https://arxiv.org/abs/1911.09154>;
+- Montealegre-Mora et al., *Certifying Numerical Decompositions of Compact
+  Group Representations*, <https://arxiv.org/abs/2101.12244>;
+- Olver, *Representations of the symmetric group are decomposable in
+  polynomial time*, <https://arxiv.org/abs/2211.12592>.
+
+No wall-clock thresholds are asserted.  Performance-sensitive paths are
+covered by small deterministic representatives so the suite remains suitable
+for GAP's `minimal`, `latest`, and `devel` CI jobs.

tst/centralizer-properties.tst

@@ -0,0 +1,53 @@
+# Centralizer and class-sum properties motivated by symmetry-reduced SDPs.
+
+gap> START_TEST("centralizer-properties.tst");
+
+gap> G := DihedralGroup(8);; irreps := IrreducibleRepresentations(G);;
+gap> blockRep := DirectSumOfRepresentations([irreps[4], irreps[4], irreps[5]]);;
+gap> A := REPN_TEST_UpperUnitriangularMat(4);;
+gap> rho := REPN_TEST_ConjugateRepresentation(blockRep, A);;
+gap> centralizer := CentralizerOfRepresentation(rho);;
+gap> Length(centralizer) = 2^2 + 1^2;
+true
+gap> ForAll(centralizer, C -> REPN_TEST_CommutesWithRepresentation(rho, C));
+true
+gap> RankMat(List(centralizer, Flat)) = Length(centralizer);
+true
+gap> centSpace := VectorSpace(Cyclotomics, centralizer);;
+gap> IdentityMat(4) in centSpace;
+true
+gap> ForAll(centralizer, X -> ForAll(centralizer, Y -> X*Y in centSpace));
+true
+
+gap> blockCentralizer := CentralizerBlocksOfRepresentation(rho);;
+gap> Length(blockCentralizer) = 5;
+true
+gap> ForAll(blockCentralizer, blocks -> Length(blocks) = 2);
+true
+gap> List(blockCentralizer[1], Length) = [2,2];
+true
+
+gap> diagRep := BlockDiagonalRepresentation(rho);;
+gap> expandedBlocks := List(blockCentralizer, BlockDiagonalMatrix@RepnDecomp);;
+gap> ForAll(expandedBlocks, C -> REPN_TEST_CommutesWithRepresentation(diagRep, C));
+true
+
+gap> unitaryRep := blockRep;; unitaryCentralizer := CentralizerOfRepresentation(unitaryRep);;
+gap> orthogonalCentralizer := OrthonormalBasis@RepnDecomp(unitaryCentralizer);;
+gap> IsOrthonormalSet(orthogonalCentralizer, InnerProduct@RepnDecomp);
+true
+gap> ForAll(ConjugacyClasses(G), class -> ClassSumCentralizer(unitaryRep, class, orthogonalCentralizer) = Sum(class, g -> Image(unitaryRep,g)));
+true
+gap> ForAll(ConjugacyClasses(G), class -> ClassSumCentralizerNC(unitaryRep, class, orthogonalCentralizer) = Sum(class, g -> Image(unitaryRep,g)));
+true
+
+gap> sizes := [rec(dimension:=1,nblocks:=2), rec(dimension:=3,nblocks:=1)];;
+gap> standardBlocks := SizesToBlocks(sizes);;
+gap> Length(standardBlocks) = 5;
+true
+gap> ForAll(standardBlocks, blocks -> DimensionsMat(BlockDiagonalMatrix@RepnDecomp(blocks)) = [5,5]);
+true
+gap> RankMat(List(standardBlocks, blocks -> Flat(BlockDiagonalMatrix@RepnDecomp(blocks)))) = 5;
+true
+
+gap> STOP_TEST("centralizer-properties.tst", 1);

tst/core-utils.tst

@@ -0,0 +1,69 @@
+# Deterministic tests for small matrix, character, and representation helpers.
+
+gap> START_TEST("core-utils.tst");
+
+gap> Take@RepnDecomp([ 1, 2, 3 ], 0) = [];
+true
+gap> Take@RepnDecomp([ 1, 2, 3 ], 2) = [ 1, 2 ];
+true
+gap> Drop@RepnDecomp([ 1, 2, 3 ], 0) = [ 1, 2, 3 ];
+true
+gap> Drop@RepnDecomp([ 1, 2, 3 ], 3) = [];
+true
+
+gap> original := [ 1, 2, 3 ];; copies := Replicate@RepnDecomp(original, 2);;
+gap> copies[1][1] := 9;; copies[2] = original and original = [ 1, 2, 3 ];
+true
+gap> Replicate@RepnDecomp(7, 3) = [ 7, 7, 7 ];
+true
+
+gap> BlockDiagonalMatrix@RepnDecomp([]) = [];
+true
+gap> BlockDiagonalMatrix@RepnDecomp([ IdentityMat(2), IdentityMat(3) ]) = IdentityMat(5);
+true
+gap> BlockDiagonalMatrix@RepnDecomp([ [[2]], [], [[3,0],[0,4]] ]) = DiagonalMat([2,3,4]);
+true
+
+gap> M := DiagonalMat([ 1, 2, 3, 4 ]);;
+gap> ExtractBlock@RepnDecomp(M, 1, 1, 2) = DiagonalMat([1,2]);
+true
+gap> ExtractBlock@RepnDecomp(M, 1, 2, 2) = NullMat(2,2);
+true
+gap> ExtractBlock@RepnDecomp(M, 2, 2, 2) = DiagonalMat([3,4]);
+true
+
+gap> V := VectorSpace(Rationals, [[1,0],[0,1]], [0,0]);;
+gap> W := MatrixImage@RepnDecomp([[1,0],[0,0]], V);;
+gap> Dimension(W) = 1 and [1,0] in W and not [0,1] in W;
+true
+
+gap> IsOrthonormalSet([[1,0],[0,1]], function(x,y) return x*y; end);
+true
+gap> IsOrthonormalSet([[1,0],[1,0]], function(x,y) return x*y; end);
+false
+
+gap> G := SymmetricGroup(3);; irreps := IrreducibleRepresentations(G);;
+gap> ForAll(irreps, IsFiniteGroupLinearRepresentation);
+true
+gap> ForAll(irreps, rho -> DegreeOfRepresentation(rho) = Length(Image(rho, One(G))));
+true
+gap> chars := List(irreps, CharacterOfRepresentation@RepnDecomp);;
+gap> ForAll([1..Length(chars)], i -> ForAll([1..Length(chars)], function(j) if i=j then return InnerProductOfCharacters@RepnDecomp(chars[i],chars[j],G)=1; else return InnerProductOfCharacters@RepnDecomp(chars[i],chars[j],G)=0; fi; end));
+true
+
+gap> direct := DirectSumOfRepresentations([irreps[1], irreps[3]]);;
+gap> DegreeOfRepresentation(direct) = 3;
+true
+gap> IrrVectorOfRepresentation@RepnDecomp(direct, chars) = [1,0,1];
+true
+
+gap> perm := ActionHomomorphism(G, [1..3]);;
+gap> IsFiniteGroupPermutationRepresentation(perm);
+true
+gap> linear := PermToLinearRep(perm);;
+gap> IsFiniteGroupLinearRepresentation(linear) and DegreeOfRepresentation(linear) = 3;
+true
+gap> ForAll(GeneratorsOfGroup(G), g -> Image(linear,g) = PermutationMat(Image(perm,g),3));
+true
+
+gap> STOP_TEST("core-utils.tst", 1);

tst/decomposition-properties.tst

@@ -0,0 +1,68 @@
+# Property tests for exact decomposition and simultaneous block diagonalization.
+
+gap> START_TEST("decomposition-properties.tst");
+
+gap> G := SymmetricGroup(3);; irreps := IrreducibleRepresentations(G);;
+gap> blockRep := DirectSumOfRepresentations([irreps[1], irreps[3], irreps[3]]);;
+gap> A := REPN_TEST_UpperUnitriangularMat(5);;
+gap> rho := REPN_TEST_ConjugateRepresentation(blockRep, A);;
+gap> info := REPN_ComputeUsingSerre(rho : irreps := irreps);;
+gap> REPN_TEST_IsBlockDiagonalization(rho, info);
+true
+gap> REPN_TEST_IsCollectedDecomposition(rho, Filtered(info.decomposition, x -> Length(x) > 0));
+true
+gap> List(Filtered(info.decomposition, x -> Length(x) > 0), Length) = [1,2];
+true
+gap> Sum(Flat(info.decomposition), summand -> Length(summand.basis)) = 5;
+true
+gap> Length(info.centralizer_basis) = 1^2 + 2^2;
+true
+
+gap> collected := IrreducibleDecompositionCollected(rho);;
+gap> REPN_TEST_IsCollectedDecomposition(rho, collected);
+true
+gap> REPN_TEST_IsInvariantDecomposition(rho, IrreducibleDecomposition(rho));
+true
+
+gap> rhoK := REPN_TEST_ConjugateRepresentation(blockRep, A);;
+gap> infoK := REPN_ComputeUsingSerre(rhoK : irreps := irreps, use_kronecker);;
+gap> REPN_TEST_IsBlockDiagonalization(rhoK, infoK);
+true
+gap> REPN_TEST_IsCollectedDecomposition(rhoK, Filtered(infoK.decomposition, x -> Length(x) > 0));
+true
+
+gap> rhoAlt := REPN_TEST_ConjugateRepresentation(blockRep, A);;
+gap> infoAlt := REPN_ComputeUsingMyMethod(rhoAlt : irreps := irreps);;
+gap> REPN_TEST_IsBlockDiagonalization(rhoAlt, infoAlt);
+true
+gap> REPN_TEST_IsCollectedDecomposition(rhoAlt, infoAlt.decomposition);
+true
+
+gap> canonical := CanonicalDecomposition(REPN_TEST_ConjugateRepresentation(blockRep, A) : irreps := irreps);;
+gap> SortedList(List(canonical, Dimension)) = [1,4];
+true
+gap> ForAll(canonical, V -> IsGInvariant@RepnDecomp(rho, V));
+true
+gap> RankMat(Concatenation(List(canonical, V -> List(Basis(V))))) = 5;
+true
+
+gap> T := Group(());; trivial3 := FuncToHom@RepnDecomp(T, g -> IdentityMat(3));;
+gap> trivialDecomp := IrreducibleDecomposition(trivial3);;
+gap> Length(trivialDecomp) = 3 and REPN_TEST_IsInvariantDecomposition(trivial3, trivialDecomp);
+true
+
+gap> H := SymmetricGroup(4);; hirreps := IrreducibleRepresentations(H);;
+gap> multiplicityFree := DirectSumOfRepresentations([hirreps[1], hirreps[3], hirreps[5]]);;
+gap> mfDecomp := IrreducibleDecomposition(multiplicityFree);;
+gap> Length(mfDecomp) = 3 and REPN_TEST_IsInvariantDecomposition(multiplicityFree, mfDecomp);
+true
+
+gap> perm := ActionHomomorphism(G, [1..3]);;
+gap> permCanonical := CanonicalDecomposition(perm : irreps := irreps);;
+gap> SortedList(List(permCanonical, Dimension)) = [1,2];
+true
+gap> linearPerm := PermToLinearRep(perm);;
+gap> ForAll(permCanonical, V -> IsGInvariant@RepnDecomp(linearPerm, V));
+true
+
+gap> STOP_TEST("decomposition-properties.tst", 1);

tst/isomorphism-properties.tst

@@ -0,0 +1,40 @@
+# Deterministic inputs for all explicit-intertwiner implementations.
+
+gap> START_TEST("isomorphism-properties.tst");
+
+gap> G := SymmetricGroup(4);; irreps := IrreducibleRepresentations(G);;
+gap> rho := DirectSumOfRepresentations([irreps[1], irreps[3], irreps[3]]);;
+gap> B := REPN_TEST_UpperUnitriangularMat(5);;
+gap> tau := REPN_TEST_ConjugateRepresentation(rho, B);;
+gap> AreRepsIsomorphic(rho, tau);
+true
+gap> IsLinearRepresentationIsomorphism(B^-1, rho, tau);
+true
+gap> IsLinearRepresentationIsomorphism(NullMat(5,5), rho, tau);
+false
+
+gap> iso := LinearRepresentationIsomorphism(rho, tau);;
+gap> IsLinearRepresentationIsomorphism(iso, rho, tau);
+true
+gap> isoK := LinearRepresentationIsomorphism(rho, tau : use_kronecker);;
+gap> IsLinearRepresentationIsomorphism(isoK, rho, tau);
+true
+gap> isoOrbit := LinearRepresentationIsomorphism(rho, tau : use_orbit_sum);;
+gap> IsLinearRepresentationIsomorphism(isoOrbit, rho, tau);
+true
+gap> isoSlow := LinearRepresentationIsomorphismSlow(rho, tau);;
+gap> IsLinearRepresentationIsomorphism(isoSlow, rho, tau);
+true
+
+gap> AreRepsIsomorphic(irreps[1], irreps[2]);
+false
+gap> LinearRepresentationIsomorphism(irreps[1], irreps[2]);
+fail
+gap> AreRepsIsomorphic(rho, irreps[1]);
+false
+
+gap> H := CyclicGroup(2);; hrep := IrreducibleRepresentations(H)[1];;
+gap> AreRepsIsomorphic(irreps[1], hrep);
+false
+
+gap> STOP_TEST("isomorphism-properties.tst", 1);