...
 
Commits (4)
......@@ -4,6 +4,12 @@
# Reading the declaration part of the package.
#
# stuff from cohcfg package
ReadPackage( "RepnDecomp", "lib/cohcfg/misc.g");
ReadPackage( "RepnDecomp", "lib/cohcfg/classmatr.g");
ReadPackage( "RepnDecomp", "lib/cohcfg/SubgroupCC.g");
ReadPackage( "RepnDecomp", "lib/cohcfg/classsum.g");
ReadPackage( "RepnDecomp", "lib/utils.gd" );
ReadPackage( "RepnDecomp", "lib/serre.gd" );
ReadPackage( "RepnDecomp", "lib/centralizer.gd" );
......
......@@ -54,23 +54,25 @@ BasisChangeMatrixSimilar@ := function(X, Y)
end;
InstallGlobalFunction( BlockDiagonalRepresentationFast, function(rho, args...)
local G, char_rho_basis, irreps, isomorphic_collected, summands, new_img, g, basis_change, basis, full_space_list, current_space_list, chars, new_rho, irrep_list, r, F, ret_basis, all_sizes, centralizer_blocks;
local G, char_rho_basis, irreps, isomorphic_collected, summands, new_img, g, basis_change, basis, full_space_list, current_space_list, chars, new_rho, irrep_list, r, ret_basis, all_sizes, centralizer_blocks, rho_cent_basis, new_rho_cent_basis;
G := Source(rho);
# If we are given a list of irreps, we use them and assume that it
# is complete
irreps := [];
if Length(args) > 0 then
if Length(args) >= 1 then
irreps := args[1];
else
irreps := IrreducibleRepresentationsDixon(G);
fi;
F := Cyclotomics;
# We can also be given a field to work over
if Length(args) > 1 then
F := args[2];
# we might also be given a basis for the centraliser of rho, which
# we can use to speed up some calculations
if Length(args) >= 2 then
rho_cent_basis := args[2];
else
rho_cent_basis := fail;
fi;
# We could just use Irr(G) here, but the ordering of Irr and
......@@ -94,6 +96,16 @@ InstallGlobalFunction( BlockDiagonalRepresentationFast, function(rho, args...)
summands := List(Flat(isomorphic_collected), r -> r.rep);
# We can also compute the basis for the centralizer ring, since we
# know the block sizes and relevant dimensions
all_sizes := List([1..Size(chars)], i -> rec(dimension := chars[i][1],
nblocks := char_rho_basis[i]));
# Don't use the blocks that don't appear
centralizer_blocks := SizesToBlocks@(Filtered(all_sizes, r -> r.nblocks > 0));
new_rho_cent_basis := List(centralizer_blocks, BlockDiagonalMatrix);
# This is the block diagonal rep with minimal size blocks
new_rho := DirectSumRepList(summands);
......@@ -105,7 +117,14 @@ InstallGlobalFunction( BlockDiagonalRepresentationFast, function(rho, args...)
# representation isomorphism.
# Note: this is where the heavy lifting of the function is
basis_change := LinearRepresentationIsomorphismSlow(new_rho, rho);
# If we were given a basis for centraliser of rho, do it fast
if rho_cent_basis <> fail then
basis_change := LinearRepresentationIsomorphism(new_rho, rho);
else
basis_change := LinearRepresentationIsomorphism(new_rho, rho,
new_rho_cent_basis, rho_cent_basis);
fi;
basis := TransposedMat(basis_change);
......@@ -118,19 +137,12 @@ InstallGlobalFunction( BlockDiagonalRepresentationFast, function(rho, args...)
for irrep_list in isomorphic_collected do
current_space_list := [];
for r in irrep_list do
Add(current_space_list, VectorSpace(F, Take@(basis, r.dim)));
Add(current_space_list, VectorSpace(Cyclotomics, Take@(basis, r.dim)));
basis := Drop@(basis, r.dim);
od;
Add(full_space_list, current_space_list);
od;
# We can also compute the basis for the centralizer ring, since we
# know the block sizes and relevant dimensions
all_sizes := List([1..Size(chars)], i -> rec(dimension := chars[i][1],
nblocks := char_rho_basis[i]));
# Don't use the blocks that don't appear
centralizer_blocks := SizesToBlocks@(Filtered(all_sizes, r -> r.nblocks > 0));
return rec(basis := ret_basis,
diagonal_rep := new_rho,
......
......@@ -145,3 +145,23 @@ InstallGlobalFunction( ClassSumCentralizer, function(rho, class, cent_basis)
coeffs := ClassSumCentralizerCoeffs@(rho, class, cent_basis);
return Sum([1..Length(coeffs)], i -> coeffs[i] * cent_basis[i]);
end );
# Does the group sum (just for convenience, nothing clever)
GroupSumWithCentralizer@ := function(rho, cent_basis)
local G, classes;
G := Source(rho);
classes := ConjugacyClasses(G);
return Sum(classes, class -> ClassSumCentralizer(rho, class, cent_basis));
end;
# assumes rho is faithful permutation representation, calculates
RepresentationCentralizerPermRep@ := function(rho)
local H, T, two_orbit_reps;
H := Range(rho); # is perm group
T := CohCfgFromPermGroup(H); # computes conjugacy classes and orbitals
two_orbit_reps := CCTransversal(T); # list of reps of 2-orbits (pairs)
# the orbital matrices themselves, this gives a basis for the
# centraliser
return List(two_orbit_reps, rep -> OrbitalMatrix@(H, rep));
end;
#ReadPackage("cohcfg", "lib/classmatr.g");
# this file is currently undocumented
BindGlobal("GetSubCC", function (CC, H) # CC a CohCfg object, H a subgroup
local SubCC, reps,
i, x;
# Sanity check
if not IsCohCfg(CC) then
Error("GetSubCC: CC must be an IsCohCfg object.\n");
elif not IsSubgroup(CC.group, H) then
Error("GetSubCC: H must be a subgroup of the underlying group of CC.\n");
fi;
# Create sub-coherent configuration
SubCC := CohCfgFromPermGroup(H, CC.Omega);
SubCC.isSubCC := true;
# Determine which G-orbital each H-orbital is contained in
reps := CCTransversal(SubCC);
SubCC.parentOrbitalMapping := List([1..SubCC.dimension], x -> CCElementContainingPair(CC, reps[x]));
SubCC.childOrbitalMapping := List([1..CC.dimension], x -> Positions(SubCC.parentOrbitalMapping, x));
SubCC.parentCC := CC;
return SubCC;
# SubCC contains all members of a CohCfg object record, as well as the following:
# .isSubCC - identifies the record as a SubCohCfg object
# .parentCC - a handle for the parent CohCfg object from which SubCC was created
# .parentOrbitalMapping - a list mapping SubCC's own orbitals to those of the parent CC
# .childOrbitalMapping - a list mapping orbitals of the parent CC to lists of orbitals of SubCC
end);
BindGlobal("IsSubCC", function (CC) # CC an object
return IsCohCfg and IsBound(CC.isSubCC) and CC.isSubCC;
end);
BindGlobal("GetParentCC", function (CC) # CC a SubCC object
# Sanity check
if not IsSubCC(CC) then
Error("GetParentCC: CC must be a sub- coherent configuration object.\n");
fi;
# Get parent CC object
return CC.parentCC;
end);
This diff is collapsed.
####################################################################
# let g in G, and A_i's the orbitals of G. Computing
# g^G as a linear combination of A_i's.
####################################################################
# ClassSum(group, conjugacyclass)
# ClassSum(group, element, sizeofelementsconjugacyclass)
BindGlobal("ClassSum", function(arg)
local arcs, oinf, i, g, G, n, len, v, x, y, fps, onums, sonums, m, orbs, T;
T := arg[1];
G := CCGroup(T);
if not IsPerm(arg[2]) then
if G <> ActingDomain(arg[2]) then
Error("ClassSum: The second argument must be a conjugacy class of the underlying group of the first argument.\n");
fi;
g := Representative(arg[2]);
len := Size(arg[2]);
else # no checks are made; g is a representative, len the size of g^G
g := arg[2];
len := arg[3];
fi;
n := CCDegree(T);
v := 0*[1..T.dimension]; # coefficients of the decomposition
# if 1=0 then # skip this!!!
# fps := Filtered([1..n], x -> x^g = x); # fixed points
# # (essentially, character computation for each G-orbit)
# onums := List(fps, x -> T.orbitMapping[x]);
# sonums := Set(onums);
# m := List(sonums, x -> Size(Filtered(onums, y -> y = x)));
# orbs := List(sonums, x -> CCElementContainingPair(T, [T.orbitNumberRepresentatives[x],
# T.orbitNumberRepresentatives[x]
# ]));
#
# # g hits i-th diagonal 2-orbit orbs[i] m[i] times;
# # there are len elements in g^G, so we get len*m[i] hits;
# # by symmetry, the contribution of i-th diagonal element is len*m[i]/(length of i-th orbit)
#
# for i in [1..Length(m)] do
# v[orbs[i]] := len*m[i]/CCFiberSizes(T)[i];
# od;
# fi; # end of "skip this!!!"
# cycles of the permutation
arcs := List(Cycles(g, [1..n]), x -> rec(
#o := CCElementContainingPair(T, [x[1], x[2]]),
o := CCElementContainingPair(T, [x[Length(x)], x[1]]),
l := Length(x))
);
# g hits i-th 2-orbit arcs[i].o by a cycle of length arcs[i].l
# there are len elements in g^G, so we get len hits by a arcs[i].l-cycle;
# by symmetry, the contribution of this 2-orbit is len*arcs[i].l/(# of 1s in its nxn-matrix)
for i in [1..Length(arcs)] do
oinf := CCLocalIndex(T, arcs[i].o);
x := oinf.orbitNumber;
v[arcs[i].o] := v[arcs[i].o] + len*arcs[i].l / (
CCFiberSizes(T)[x] * T.twoOrbitNumbers[x].subdegree[oinf.twoOrbitNumber]
);
od;
return v;
end);
####################################################################
# computing all the classsums and storing them in the corr. record
####################################################################
BindGlobal("ClassSums", function(arg)
local T, x, sums, g_i, l_i, cc;
T:=arg[1];
if not IsBound(T.classSums) then
if not IsBound(arg[2]) then
cc := ConjugacyClasses(CCGroup(T));
g_i := List(cc, Representative);
l_i := List(cc, Size);
else
g_i := arg[2];
l_i := arg[3];
fi;
sums := List([1 .. Length(g_i)], x -> ClassSum(T, g_i[x], l_i[x]));
T.classSums := rec(sums := sums, reps := g_i, lengths := l_i);
fi;
return T.classSums;
end);
####################################################################
# Computing the projection of the regular representation
# onto a representation with the character \chi.
#
# We want to find
# (dim(\chi) / |G|) * \sum_C \chi(C)^* \bar C
# Here the index C in "\sum_C" ranges over conjugacy classes of G,
# and \bar C is the sum of the elements in C, after they have been
# taken through the representation from G to the basis algebra of
# the coherent configuration and from there through the isomorphism
# from the basis algebra to the intersection algebra of the CC.
#
# Note that
# \sum_C \chi(C)^* \bar C
# = \sum_C (\chi(C)^* \sum_{i=1}^d ClassSum(C)_i P_i)
# = \sum_{i=1}^d (\sum_C ClassSum(C)_i \chi(C)^*) P_i)
####################################################################
BindGlobal("ProjComp", function(T, chi)
local i, v, M, d, s_r, sums, s, c;
s_r := ClassSums(T);
return (chi[1] / Size(CCGroup(T))) * Sum([1 .. CCDimension(T)], i ->
Sum([1 .. Length(s_r.reps)], C ->
s_r.sums[C][i] * ComplexConjugate(chi[C])
) * CCIntersectionMat(T, i)
);
end);
####################################################################
# naive diagonalisation (block diagonalization of the regular representation)
####################################################################
BindGlobal("ProjIrr", function(arg)
local ct, T, x, p, G, irrs, projs, linindeprows;
linindeprows := function(V)
local M, r, len, t;
M := [];
r := RankMat(V);
len := 0;
while len < r do
t := PositionProperty(V, y -> RankMat(Concatenation(M, [y])) > len);
Append(M, [V[t]]);
V := V{[t+1..Length(V)]};
len := len+1;
od;
return M;
end;
T := arg[1];
CCPopulateCoeffs(T);
ClassSums(T);
ct := [];
if not IsBound(arg[2]) then
G := CCGroup(T);
ct := CharacterTable(G);
irrs := Irr(ct);
else
irrs := arg[2]; # TODO instead of accepting irreducibles as a list, accept a character table and connect it to the group
fi;
# G := CCGroup(T);
# if not IsBound(arg[2]) then
# ct := CharacterTable(G);
# else
# ct := arg[2];
# ConnectGroupAndCharacterTable(G, ct);
# fi;
# irrs := Irr(ct);
projs := List(irrs, x -> rec(p := ProjComp(T, x), dim := x[1]));
projs := Filtered(projs, x -> x.p <> 0*projs[1].p);
projs := List(projs, x -> rec(p := x.p, rank := RankMat(x.p), dim := x.dim));
if Sum(List(projs, G -> G.rank)) <> T.dimension then
Error("ProjIrr: wrong projectors? ");
fi;
# "stack" the individual X (which trim matrix to its block component)
return rec( # TODO this record is ugly, make it sensible (xref BlocksOfMat)
projmat := Concatenation(List(projs, G -> linindeprows(G.p))),
blocksizes := List(projs, G -> G.rank),
dims := List(projs, G -> G.dim),
ct := ct
);
end);
####################################################################
# obtaining blocks of a matrix M in the regular representation
####################################################################
BindGlobal("BlocksOfMat", function(P, M) # TODO coordinate with format returned by ProjIrr
local u, v, b, X, i, bsizes, p, V, k;
p := P.projmat*M*P.projmat^-1;
bsizes := P.blocksizes;
V := [];
i := 0;
k := 0;
for b in bsizes do
k := k + 1;
if b > 0 then
X := NullMat(b,b);
for u in [1..b] do
for v in [1..b] do
X[u][v] := p[i + u][i + v];
od;
od;
Add(V, rec(blksz := X, char := P.dims[k]));
i:=i + b;
fi;
od;
return V;
end);
# Replace and extend GRAPE's OrbitNumbers function, using native GAP functions
# G is a group
# Omega is a set upon which G acts (please make sure it is closed under action by G, as this will not be checked!)
BindGlobal("OrbitsInfo", function (G, Omega)
local orbits, mapping, representatives,
i, j;
if IsPosInt(Omega) then Omega := [1..Omega]; # backwards-compatibility with GRAPE's OrbitNumbers() syntax
elif not IsHomogeneousList(Omega) then
Error("OrbitsInfo: malformed Omega");
fi;
orbits := OrbitsDomain(G, Omega); # this returns a list of the orbits; WARNING: undefined behavior when Omega is not closed!
mapping := EmptyPlist(Length(Omega));
for i in [1..Length(orbits)] do
for j in orbits[i] do
mapping[j] := i; # element j in Omega is in the i-th orbit; WARNING: mapping domain may be too big when Omega is not closed!
od;
od;
MakeImmutable(mapping);
representatives := List(orbits, x -> x[1]); # retrieves the first element of each orbit
MakeImmutable(representatives);
return rec(
orbits := orbits, # a list of the orbits of G acting on Omega
mapping := mapping, # a mapping from Omega to orbits, i.e. mapping[i] = j iff i in Omega is in orbit j
orbitNumbers := ~.mapping, # backwards-compatibility with GRAPE's OrbitNumbers() syntax
representatives := representatives # a list whose i-th entry is a representative element from the i-th orbit of Omega
);
end);
......@@ -21,7 +21,84 @@ MatrixBasis@ := function(n)
return List(coords, make_mat);
end;
InstallGlobalFunction( LinearRepresentationIsomorphism, function(rho, tau)
# Calculates the isomorphism using the cool fact about the product
# (see below)
InstallGlobalFunction( LinearRepresentationIsomorphism, function(rho, tau, args...)
local G, n, matrix_basis, vector_basis, alpha, triv, fixed_space, A, A_vec, rho_cent_basis, tau_cent_basis, triv_proj;
# if set, these bases for the centralisers will be used to avoid
# summing over G
rho_cent_basis := fail;
tau_cent_basis := fail;
if Length(args) >= 2 then
rho_cent_basis := args[1];
tau_cent_basis := args[2];
fi;
if not AreRepsIsomorphic(rho, tau) then
return fail;
fi;
G := Source(rho);
n := DegreeOfRepresentation(rho);
# We want to find a matrix A s.t. tau(g)*A = A*rho(g) for all g in
# G. We do this by finding fixed points of the linear maps A ->
# tau(g)*A*rho(g^-1). This is done by considering the
# representation alpha: g -> (A -> tau(g)*A*rho(g^-1)) and finding a
# vector in the canonical summand corresponding to the trivial
# irrep. i.e. a vector which is fixed by all g, which is exactly
# what we want.
# The trick we use to avoid summing over G is to notice that alpha
# is actually \tau \otimes \rho^*, i.e. g -> tau(g) \otimes
# \rho(g^-1)^T.
# Standard basis for M_n(C)
matrix_basis := MatrixBasis@(n);
# the representation alpha : G -> GL(V) (V is the space of matrices)
alpha := FuncToHom@(G, g -> KroneckerProduct(Image(tau, g), TransposedMat(Image(rho, g^-1))));
# the projection of V onto V_triv, the trivial canonical summand,
# is just given by the sum over whole group of alpha(g)
#
# we can get a projection into the same space by doing the sum
# over g in G, h in G of tau(g) \otimes rho^*(h), and we can
# calculate these group sums using the centraliser bases
# if there was no basis, just sum over the whole group
if rho_cent_basis = fail or tau_cent_basis = fail then
triv_proj := Sum(G, g -> Image(alpha, g));
else
# we can just do (sum_{g in G} tau(g)) \otimes (sum_{g in G} rho^*(g))
# which still gives a projection
# The group sum for rho^* is the same as for rho, but
# transposed (the relabelling g -> g^-1 is just a bijection and ^T is linear)
triv_proj := KroneckerProduct(GroupSumWithCentralizer@(tau, tau_cent_basis),
TransposedMat(GroupSumWithCentralizer@(rho, rho_cent_basis)));
fi;
A := NullMat(n, n);
# Keep picking matrices until we get an invertible one. This would
# happen with probability 1 if we really picked uniformly random
# vectors over a ball in C^n^2.
repeat
# we pick a "random vector" and project it to get a fixed one
A_vec := triv_proj * Flat(RandomInvertibleMat(n));
A := WrapMatrix@(A_vec, n);
until RankMat(A) = n;
return A;
end );
# Calculate the iso without using the cool fact about the conjugation
# action being given by \tau \otimes \rho^*
LinearRepresentationIsomorphismNoProduct@ := function(rho, tau)
local G, n, matrix_basis, vector_basis, alpha, triv, fixed_space, A, A_vec, alpha_f;
if not AreRepsIsomorphic(rho, tau) then
......@@ -83,7 +160,7 @@ InstallGlobalFunction( LinearRepresentationIsomorphism, function(rho, tau)
until RankMat(A) = n;
return A;
end );
end;
# checks if it is the case that for all g in G, tau(g)*A = A*rho(g)
......
# These are the implementations of Serre's formulas from his book
# Linear Representations of Finite Groups.
# This is used for speed in some special cases
LoadPackage("cohcfg");
MatrixImage@ := function(p, V)
local F;
......