Skip to content

Speed up the linear algebra part of character table calculations for (large) groups #4773

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
hulpke merged 6 commits into from
Mar 2, 2022
Merged

Speed up the linear algebra part of character table calculations for (large) groups #4773

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b1bd112
PositionNonZero for vector objects
hulpke Feb 16, 2022
f65e003
ENHANCE Methods for matrix objects
hulpke Feb 16, 2022
ce880cf
ENHANCE: Matrix objects for ZmodnZ
hulpke Feb 16, 2022
5eaa5a4
ENHANCE: Use matrix objects in Dixon-Schneider
hulpke Feb 16, 2022
eec16ca
Test and changes in test output and examples
hulpke Feb 16, 2022
e689844
HACK: Rank method higher than in package
hulpke Feb 17, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 0 additions & 33 deletions doc/ref/ctbl.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -515,39 +515,6 @@ gap> Length( ro.irreducibles );
gap> DxIncludeIrreducibles( d, ro.irreducibles );
]]></Example>
<P/>
The tensor products of the nonlinear characters among each other are reduced
with the irreducible characters.
The result is split according to the spaces found, which yields characters
of smaller norms, but no new irreducibles.
<P/>
<Example><![CDATA[
gap> nlc:= Filtered( d.irreducibles, i -> i[1] > 1 );;
gap> t:= Tensored( nlc, nlc );;
gap> ro:= ReducedCharacters( c, d.irreducibles, t );; ro.irreducibles;
[ ]
gap> List( ro.remainders, i -> ScalarProduct( c, i, i) );
[ 2, 2, 4, 4, 4, 4, 13, 13, 18, 18, 19, 21, 21, 36, 36, 29, 34, 34,
42, 34, 48, 54, 62, 68, 68, 78, 84, 84, 88, 90, 159, 169, 169, 172,
172, 266, 271, 271, 268, 274, 274, 280, 328, 373, 373, 456, 532,
576, 679, 683, 683, 754, 768, 768, 890, 912, 962, 1453, 1453, 1601,
1601, 1728, 1739, 1739, 1802, 2058, 2379, 2414, 2543, 2744, 2744,
2920, 3078, 3078, 4275, 4275, 4494, 4760, 5112, 5115, 5115, 5414,
6080, 6318, 7100, 7369, 7369, 7798, 8644, 10392, 12373, 12922,
14122, 14122, 18948, 21886, 24641, 24641, 25056, 38942, 44950,
78778 ]
gap> t:= SplitCharacters( d, ro.remainders );;
gap> List( t, i -> ScalarProduct( c, i, i ) );
[ 2, 2, 4, 2, 2, 4, 4, 3, 6, 5, 5, 9, 9, 4, 12, 13, 18, 18, 20, 18,
20, 24, 26, 32, 32, 16, 42, 59, 69, 69, 72, 72, 36, 72, 78, 78, 84,
122, 117, 127, 117, 127, 64, 132, 100, 144, 196, 256, 456, 532,
576, 679, 683, 683, 754, 768, 768, 890, 912, 962, 1453, 1453, 1601,
1601, 1728, 1739, 1739, 1802, 2058, 2379, 2414, 2543, 2744, 2744,
2920, 3078, 3078, 4275, 4275, 4494, 4760, 5112, 5115, 5115, 5414,
6080, 6318, 7100, 7369, 7369, 7798, 8644, 10392, 12373, 12922,
14122, 14122, 18948, 21886, 24641, 24641, 25056, 38942, 44950,
78778 ]
]]></Example>
<P/>
Finally we calculate the characters induced from all cyclic subgroups and
obtain the missing irreducibles by applying the LLL-algorithm to them.
<P/>
Expand Down
129 changes: 56 additions & 73 deletions lib/ctblgrp.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ local D,C,cl,pl;
SortParallel(cl,pl,ClassComparison);
D.perm:=PermList(pl);
D.permlist:=pl;
D.invpermlist:=ListPerm(D.perm^-1,Length(D.permlist));
D.currentInverseClassNo:=0;

D.characterTable:=C;
Expand All @@ -104,32 +105,6 @@ local D,C,cl,pl;
return D;
end);

# #############################################################################
# ##
# #F DxPowerClass(<D>,<cl>,<pow>) . . . . . . . . . . . . . evaluate powermap
# ##
# DxPowerClass := function(D,nu,power)
# local p,primes,cl;
# cl:=nu;
# power:=power mod D.characterTable.orders[cl];
# if power=0 then
# return 1;
# elif power=1 then
# return cl;
# else
# primes:=Factors(power);
# for p in primes do
# if not IsBound(D.characterTable.powermap[p]) then
# return D.ClassElement(D,
# Representative(D.classes[nu])^power);
# else
# cl:=D.characterTable.powermap[p][cl];
# fi;
# od;
# return cl;
# fi;
# end;

#############################################################################
##
#F DxCalcAllPowerMaps(<D>) . . . . . . . calculate power maps for char.table
Expand Down Expand Up @@ -270,13 +245,11 @@ end;

#############################################################################
##
#F DxNiceBasis(v,space)
#F DxIsInSpace(v,space)
##
DxIsInSpace := function(v,r)
if not IsBound(r.space) then
r.space:=VectorSpace(Field(r.base[1]),r.base);
fi;
return v in r.space;
local nr,nc,m,i,j;
return SolutionMat(Matrix(BaseDomain(r.base[1]),r.base),v)<>fail;
Comment thread
hulpke marked this conversation as resolved.
end;

#############################################################################
Expand All @@ -286,9 +259,9 @@ end;
DxNiceBasis := function(d,r)
local b;
if not IsBound(r.niceBasis) then
b:=List(r.base,i->i{d.permlist}); # copied and permuted according to order
b:=Matrix(BaseDomain(r.base[1]),List(r.base,i->i{d.permlist})); # copied and permuted according to order
TriangulizeMat(b);
b:=List(b,i->Permuted(i,d.perm)); # permuted back
b:=List(b,i->i{d.invpermlist}); # permuted back
r.niceBasis:=Immutable(b);
fi;
Assert(1,Length(r.niceBasis)=Length(r.base));
Expand Down Expand Up @@ -336,6 +309,7 @@ DxRegisterModularChar := function(D,c)
local d,p;
# it may happen,that an irreducible character will be registered twice:
# 2-dim space,1 Orbit,combinatoric split. Avoid this!
if IsPlistRep(c) then c:=Vector(c); fi;
Copy link
Copy Markdown
Member

@fingolfin fingolfin Mar 1, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are the entries of D.modulars modified later on, or will they be immutable? In the latter case, you could also try

Suggested change
if IsPlistRep(c) then c:=Vector(c); fi;
if IsPlistRep(c) then c:=ImmutableVector(SOME_RING, c); fi;

provided SOME_RING can be derived from D or can be passed in some other way to this function.

Even if it is mutable, it might be a good idea to pass the basering and/or filter: The problem with Vector(c) is that if one passes a list-of-GF(4) elements to it, then sometimes one will get a GF(4) vector yet sometimes a GF(2) vector, which usually is not desirable

if not(c in D.modulars) then
Add(D.modulars,c);
d:=Int(c[1]);
Expand Down Expand Up @@ -412,25 +386,28 @@ local i,pcomp,m,r,D,neue,tm,news,opt;
od;

if opt<>true then
pcomp:=NullspaceMat(D.projectionMat*TransposedMat(D.modulars));
for i in [1..Length(D.raeume)] do
r:=D.raeume[i];
if r.dim = Length(r.base[1]) then
# trivial case: Intersection with full space in the beginning
r:=rec(base:=pcomp);
else
r:=rec(base:=SumIntersectionMat(pcomp,r.base)[2]);
fi;
r.dim:=Length(r.base);
# note stabilizer
if IsBound(D.raeume[i].stabilizer) then
r.stabilizer:=D.raeume[i].stabilizer;
fi;
if r.dim>0 then
DxActiveCols(D,r);
fi;
D.raeume[i]:=r;
od;
pcomp:=NullspaceMat(D.projectionMat*TransposedMatImmutable(Matrix(D.field,D.modulars)));
if Length(pcomp)>0 then
for i in [1..Length(D.raeume)] do
r:=D.raeume[i];
if r.dim = Length(r.base[1]) then
# trivial case: Intersection with full space in the beginning
r:=rec(base:=pcomp);
else
r:=rec(base:=SumIntersectionMat(pcomp,
Matrix(BaseDomain(r.base[1]),r.base))[2]);
Copy link
Copy Markdown
Contributor

@ThomasBreuer ThomasBreuer Feb 23, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it perhaps make sense to make sure that the base component is a MatrixObj from the beginning, then calling Matrix here would not be necessary?

Copy link
Copy Markdown
Contributor Author

@hulpke hulpke Feb 23, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It could be. In the old setup a matrix is just a list of vectors, and there was no need to distinguish the concepts. Here, one has to, and it seemed naively that a basis should be a list of vectors rather than a matrix.

fi;
r.dim:=Length(r.base);
# note stabilizer
if IsBound(D.raeume[i].stabilizer) then
r.stabilizer:=D.raeume[i].stabilizer;
fi;
if r.dim>0 then
DxActiveCols(D,r);
fi;
D.raeume[i]:=r;
od;
fi;
fi;
D.raeume:=Filtered(D.raeume,i->i.dim>0);
end );
Expand Down Expand Up @@ -627,7 +604,7 @@ local ml,nu,ret,r,p,v,alo,ofs,orb,i,j,inv,b;
alo:=Length(b);
od;
od;
inv:=b^(-1);
inv:=Matrix(BaseDomain(b[1]),b)^(-1);
for i in ml do
v:=i*inv;
for r in [1..Length(D.raeume)] do
Expand All @@ -636,7 +613,7 @@ local ml,nu,ret,r,p,v,alo,ofs,orb,i,j,inv,b;
p[j]:=v[j];
od;
p:=p*b;
if p<>nu then
if not IsZero(p) then
AddSet(ret,DxLiftCharacter(D,p));
fi;
od;
Expand All @@ -653,7 +630,7 @@ DxEigenbase := function(M,f)
local dim,i,k,eigenvalues,base,minpol,bases;
k:=Length(M);

minpol:=MinimalPolynomial(f,M);
minpol:=MinimalPolynomial(BaseDomain(M),M);

Assert(2,IsDuplicateFree(RootsOfUPol(minpol)));
eigenvalues:=Set(RootsOfUPol(minpol));
Expand All @@ -664,7 +641,6 @@ DxEigenbase := function(M,f)
if base=[] then
Error("This can`t happen: Wrong Eigenvalue ???");
else
#TriangulizeMat(base);
dim:=dim+Length(base);
Add(bases,base);
fi;
Expand Down Expand Up @@ -705,7 +681,7 @@ InstallGlobalFunction(SplitStep,function(D,bestMat)
if ForAny(raeume,i->i.dim>1) then
bestMatCol:=D.requiredCols[bestMat];
bestMatSplit:=D.splitBases[bestMat];
M:= NullMat( k, k, 0 );
M:= ZeroMatrix(Integers,k,k);
Info(InfoCharacterTable,1,"Matrix ",bestMat,",Representative of Order ",
Order(D.classreps[bestMat]),
",Centralizer: ",D.centralizers[bestMat]);
Expand All @@ -729,15 +705,21 @@ InstallGlobalFunction(SplitStep,function(D,bestMat)
dim:=raum.dim;
# cut out the 'active part' for computation of an eigenbase
activeCols:=DxActiveCols(D,raum);
N:= NullMat( dim, dim, o );
N:= ZeroMatrix(f,dim,dim);
for row in [1..dim] do
Row:=base[row]*M;
#Row:=base[row]*M;
Row:=CompatibleVector(M);
for col in [1..Length(Row)] do
Row[col]:=base[row,col];
od;
Row:=Row*M;
Comment on lines +711 to +715
Copy link
Copy Markdown
Member

@fingolfin fingolfin Mar 1, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note to self: I think we need a function doing this; i.e. which takes a matrix M, and a list of co/a vector v, and returns a new vector in "compatible" representation.

Though in this case, perhaps base could be turned into a list of vectors of the right type from the start, i.e. as part of its construction?

Copy link
Copy Markdown
Contributor

@ThomasBreuer ThomasBreuer Mar 1, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@fingolfin What about the syntax Vector( <M>, <v> )?

for col in [1..dim] do
N[row][col]:=Row[activeCols[col]];
N[row,col]:=Row[activeCols[col]];
od;
od;
eigen:=DxEigenbase(N,f);
# Base umrechnen
base:=Matrix(BaseDomain(base[1]),base);
eigenbase:=List(eigen.base,i->List(i,j->j*base));
#eigenvalues:=List(eigen.values,i->i/D.classiz[bestMat]);

Expand Down Expand Up @@ -1092,7 +1074,7 @@ local i,newRaeume,raum,neuer,j,ch,irrs,mods,incirrs,incmods,nb,rt,neuc;
Info(InfoCharacterTable,1,"TwoDimSpace image");
nb:=D.asCharacterMorphism(raum.base[1],j);
neuc:=List(irrs,i->D.asCharacterMorphism(i,j));
if not ForAny([i+1..Length(D.raeume)],i->nb in D.raeume[i].space) then
if not ForAny([i+1..Length(D.raeume)],i->DxIsInSpace(nb,D.raeume[i])) then
incirrs:=Union(incirrs,neuc);
incmods:=Union(incmods,List(mods,i->D.asCharacterMorphism(i,j)));
else
Expand Down Expand Up @@ -1747,7 +1729,7 @@ end;
StandardClassMatrixColumn := function(D,M,r,t)
local c,gt,s,z,i,T,w,e,j,p,orb;
if t=1 then
M[D.inversemap[r]][t]:=D.classiz[r];
M[D.inversemap[r],t]:=D.classiz[r];
else
orb:=DxGaloisOrbits(D,r);
z:=D.classreps[t];
Expand All @@ -1756,9 +1738,9 @@ StandardClassMatrixColumn := function(D,M,r,t)
p:=RepresentativeAction(Stabilizer(orb.group,r),c,t);
if p<>fail then
# was the first column of the galois class active?
if ForAny(M,i->i[c]>0) then
if ForAny([1..NrRows(M)],i->M[i,c]>0) then
for i in D.classrange do
M[i^p][t]:=M[i][c];
M[i^p,t]:=M[i,c];
od;
Info(InfoCharacterTable,2,"by GaloisImage");
return;
Expand Down Expand Up @@ -1787,19 +1769,19 @@ StandardClassMatrixColumn := function(D,M,r,t)
else # only strong test possible
s:=D.ClassElement(D,e);
fi;
M[s][t]:=M[s][t]+T[2][i];
M[s,t]:=M[s,t]+T[2][i];
od;
if w then # weak discrimination possible ?
gt:=Set(Filtered(orb.orbits,i->Length(i)>1));
for i in gt do
if i[1] in orb.identifees then
# were these classes detected weakly ?
e:=M[i[1]][t];
e:=M[i[1],t];
if e>0 then
Info(InfoCharacterTable,2,"GaloisIdentification ",i,": ",e);
fi;
for j in i do
M[j][t]:=e/Length(i);
M[j,t]:=e/Length(i);
od;
fi;
od;
Expand All @@ -1808,7 +1790,7 @@ StandardClassMatrixColumn := function(D,M,r,t)
for i in [1..Length(T[1])] do
s:=D.ClassElement(D,T[1][i] * z);
Unbind(T[1][i]);
M[s][t]:=M[s][t]+T[2][i];
M[s,t]:=M[s,t]+T[2][i];
od;
fi;
fi;
Expand Down Expand Up @@ -1997,7 +1979,7 @@ local G, # group
D.field:=f;
D.one:=One(f);
D.z:=z;
r:=rec(base:=Immutable( IdentityMat(k,D.one) ),dim:=k);
r:=rec(base:=List(IdentityMat(k,D.one),Vector),dim:=k);
D.raeume:=[r];

# Galois group operating on the columns
Expand All @@ -2011,11 +1993,12 @@ local G, # group
D.galOp:=[];
D.irreducibles:=[];

M:= IdentityMat(k);
fk:=D.one/(Size(G) mod prime);
M:= ZeroMatrix(D.field,k,k);
for i in [1..k] do
M[i][i]:=D.classiz[i] mod prime;
M[i,i]:=(D.classiz[i] mod prime)*D.one*fk;
od;
D.projectionMat:=M*(D.one/(Size(G) mod prime));
D.projectionMat:=M;

#if (USECTPGROUP or Size(G)<2000 or k*10>=Size(G))
# and IsBound(G.isAgGroup) and G.isAgGroup
Expand Down Expand Up @@ -2043,7 +2026,7 @@ local G, # group
# CharacterMorphisms.
D.raeume[1].stabilizer:=CharacterMorphismGroup(D);
m:=First(D.classes,i->Size(i)>1);
if Size(m)>8 then
if m<>fail and Size(m)>8 then
D.maycent:=true;
fi;
fi;
Expand Down
10 changes: 10 additions & 0 deletions lib/matobj.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,8 @@ function( basedomain )
return IsGF2VectorRep;
elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then
return Is8BitVectorRep;
elif IsFinite(basedomain) and IsZmodnZObj(Zero(basedomain)) then
return IsZmodnZVectorRep;
fi;
return IsPlistVectorRep;
end);
Expand All @@ -156,6 +158,8 @@ function( basedomain )
return IsGF2MatrixRep;
elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then
return Is8BitMatrixRep;
elif IsFinite(basedomain) and IsZmodnZObj(Zero(basedomain)) then
return IsZmodnZMatrixRep;
fi;
return IsPlistMatrixRep;
end);
Expand Down Expand Up @@ -293,6 +297,10 @@ InstallMethod( Matrix,

InstallMethod( Matrix,
[ IsSemiring, IsList ],
# rank higher than a method in the semigroups package, which otherwise jumps
# in and causes an error when testing
# line 318 of semigroups-3.4.0/gap/elements/semiringmat.gi
20,
Comment on lines +300 to +303
Copy link
Copy Markdown
Member

@fingolfin fingolfin Feb 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we'll have to sort it out together with the Semigroups folks, though, don't want to break that package (paging @james-d-mitchell).

function( R, list )
if Length(list) = 0 then
Error( "list must be not empty" );
Expand Down Expand Up @@ -1345,6 +1353,8 @@ InstallMethod( Characteristic,
[ IsMatrixOrMatrixObj ],
M -> Characteristic( BaseDomain( M ) ) );

InstallOtherMethod(DefaultFieldOfMatrix,[IsMatrixOrMatrixObj],BaseDomain);
Copy link
Copy Markdown
Member

@fingolfin fingolfin Feb 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The documentation for DefaultFieldOfMatrix explicitly states that either a field or fail is returned; this may return other things.

Is it really necessary to make this overload?

Copy link
Copy Markdown
Contributor Author

@hulpke hulpke Feb 18, 2022 •
edited by fingolfin
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have been using DefaultFiledOfMatrix in my code to get the definition field of an (ordinary) matrix. How would we get it in a uniform way, regardless whether an object is IsMatrix or IsMatrixObj. Will BaseDomain work for the "classical" matrices?

As for "necessary", it is used in the standard matrix routines.

Copy link
Copy Markdown
Member

@fingolfin fingolfin Mar 1, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes BaseDomain works for "classical" matrices (but there was a bug: it only looked .

But since it has different semantics, it is not so clear what the "correct" answer should be for e.g. an identity matrix with integer entries: should it be Integers or Rationals? This is of course clear with DefaultFiledOfMatrix. But BaseDomain currently returns Integers. We could of course change that, but in the end, it'll always be a guessing game; which is after all one of the primary reasons why we want MatrixObj: so that we can stop guessing and start specifying what we want.



#############################################################################
##
Expand Down
19 changes: 19 additions & 0 deletions lib/matobjnz.gd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##

# represent vectors/matrices over Z/nZ by nonnegative integer lists
# in the range [0..n-1], but reduce after
# arithmetic. This way avoid always wrapping all entries separately

DeclareRepresentation( "IsZmodnZVectorRep",
IsVectorObj and IsPositionalObjectRep and HasBaseDomain, [] );

DeclareRepresentation( "IsZmodnZMatrixRep",
IsRowListMatrix and IsPositionalObjectRep and HasBaseDomain, [] );
Loading