[GAP Forum] A question about matrices
Jean Michel
jmichel at math.jussieu.fr
Fri May 12 10:22:23 BST 2006
On Thu, May 11, 2006 at 08:33:38PM -0700, D N wrote:
> Suppose A and B are two square matrices of the same dimension (say 1000). Is there an efficient way to check whether A = PBP^t
> for some permutation matrix P (P^t = transpose of P).
I have had a need to solve exactly that problem. 1000x1000 is very large,
especially if the matrices have a large stabilizer in the permutation group,
but I have had success with the heuristic routine below (the main new
idea in PermMatMat is DistHelpedRepresentativeOperation)
#############################################################################
##
#A matrix.g CHEVIE library Jean Michel
##
#Y Copyright (C) 2002 - 2006 University Paris VII.
##
## This file contains routines which augment the capabilities of GAP3
## in mainpulating matrices developed during the writing of Chevie.
## These routines should be in some other part of the GAP library, as
## they are not specific to Chevie
##
#############################################################################
##
#F DecomposedMat( <M> ) . . . . . . . . . Find if the square matrix M
#F admits a block decomposition.
##
## Define a graph G with vertices [1..Length(M)] and with an edge
## between i and j if either M[i][j] or M[j][i] is non-zero.
## DecomposedMat return a list of lists l such that l[1],l[2], etc..
## are the vertices in each connected component of G. In other words,
## M{l[1]}{l[1]},M{l[2]}{l[2]},etc... are blocks of M.
##
DecomposedMat:=function(M)
local l, i, j, k, cl, ci, cj, erg;
l:=Length(M);
# cl contains numbers for the equivalence classes, first
# each element forms a class:
cl:=[1..l];
for i in [1..l] do
for j in [i+1..l] do
# if new relation i~j and classes are different then join classes:
if (M[i][j]<>0 or M[j][i]<>0) and cl[i]<>cl[j] then
ci:=cl[i]; cj:=cl[j];
for k in [1..l] do
if cl[k]=cj then
cl[k]:=ci;
fi;
od;
fi;
od;
od;
# extract the classes from information in cl:
erg:=[];
for i in [1..l] do
if IsBound(erg[cl[i]]) then
AddSet(erg[cl[i]],i);
else
erg[cl[i]]:=[i];
fi;
od;
return Set(erg);
end;
#############################################################################
##
#F IsLowerTriangularMat( <M> ) . . . . . true iff <M> is lower triangular
##
IsLowerTriangularMat:=function(M)local n;n:=Length(M);
return ForAll([1..n-1],i->ForAll([i+1..n],j->M[i][j]=0*M[i][j]));
end;
#############################################################################
##
#F IsDiagonalMat( <mat> ) . . . . . . . . . . . . true if <mat> is diagonal
##
IsDiagonalMat:=function(M)local i, j;
for i in [1..Length(M)] do
for j in [1..Length(M[1])] do
if i<>j and M[i][j]<>0*M[i][j] then return false;fi;
od;
od;
return true;
end;
#############################################################################
##
#F DiagonalMat( <m1>,...,<mn> ) diagonal block-matrix
##
## returns the block diagonal matrix with diagonal entries
## <m1> ... <mn>
## m_i may be rectangular or empty; blocks of size 1x1 may be given as scalars
## DiagonalMat( <v> ) where v is a list of scalars is like DiagonalMat(v1,..vn)
##
DiagonalMat:=function(arg)local res,fr,fc,r,c,m;
if Length(arg)=1 and IsList(arg[1]) and Length(arg[1])>0 and not
IsList(arg[1][1]) then arg:=arg[1];fi;
fr:=function(m)if IsMat(m) then return Length(m);else return 1;fi;end;
fc:=function(m)if not IsMat(m) then return 1;
elif Length(m)=0 then return 0; else return Length(m[1]);fi;end;
res:=NullMat(Sum(arg,fr),Sum(arg,fc));
r:=0;c:=0;
for m in arg do
if IsMat(m) then res{r+[1..fr(m)]}{c+[1..fc(m)]}:=m;
else res[r+1][c+1]:=m;fi;
r:=r+fr(m);c:=c+fc(m);
od;
return res;
end;
#############################################################################
##
#F IsNormalizing( <lst>, <mat> ) . . . . . . . true if matrix <mat> lets
#F set <lst> of vectors invariant
##
IsNormalizing:=function(l,M)
return Set(l*M)=Set(l);
end;
############################################################################
# EigenvaluesMat: eigenvalues of a square matrix over the cyclotomics
# whose eigenvalues are roots of unity or 0. false is returned for other
# matrices.
#
EigenvaluesMat:=function(mat)local p;
p:=CycPol(CharacteristicPolynomial(mat));
if not IsCyc(p.coeff) then return false;fi;
return Concatenation(List(p.vcyc,x->List([1..x[2]],i->
E(Denominator(x[1]))^Numerator(x[1]))));
end;
############################################################################
#
# ExteriorPower(A,m) . . returns the m-th exterior power of square matrix m
#
ExteriorPower:=function(A,m) local comb;
comb:=Combinations([1..Length(A)],m);
return List(comb,i->List(comb,j->DeterminantMat(A{i}{j})));
end;
#############################################################################
##
#F KroneckerMatrix( <M1>[, <M2>[,...]]) . . . . . . like 'KroneckerProduct',
#F but works with arbitrary number of arguments
##
KroneckerMatrix:=function(arg)
local res, i;
res:=arg[1];
for i in [2..Length(arg)] do
res:=MatricesOps.KroneckerProduct(res,arg[i]);
od;
return res;
end;
#############################################################################
##
#F ProportionalityCoefficient( <v>, <w>) . . . . . . 'v/w'
#F v and w are vectors. Returns scalar l such that v=l*w if such
## exists, false otherwise.
##
ProportionalityCoefficient:=function(v,w)local i;
i:=PositionProperty(w,x->x<>0);
if i=false then if w=v then return 0;else return false;fi;fi;
i:=v[i]/w[i];
if v<>w*i then return false;fi;
return i;
end;
#############################################################################
##
#F RepresentativeDiagonalConjugations(M,N) . . tests if square matrices
# M and N are conjugate by a diagonal matrix
#
# find a diagonal matrix D=DiagonalMat(d=[1,d2,..,dn]) such that N=M^D
# thus N[i][j]=M[i][j]*d[j]/d[i]
#
RepresentativeDiagonalConjugation:=function(M,N)local d,n,i,j,c;
d:=M[1]*0;d[1]:=1;n:=Length(M);
for i in [1..n] do
for j in [i+1..n] do
# Print("i=",i," j=",j," d=",d,"\n");
if M[i][j]<>0 then
if N[i][j]=0 then return false;fi;
if d[i]<>0 then
c:=d[i]*N[i][j]/M[i][j];
if d[j]<>0 then if c<>d[j] then return false;fi;
else d[j]:=c;
fi;
fi;
fi;
od;
od;
if 0 in d then return false;fi;
if N<>M^DiagonalMat(d) then return false;fi;
return d;
end;
#############################################################################
##
#F OnMatrices(M,p) . . . Simultaneaous action on rows and columns of a
# permutation p on the square matrix M
#
OnMatrices:=function(M,p)return Permuted(List(M,y->Permuted(y,p)),p);end;
#############################################################################
##
# MatStab([g,]M[,extra]) . . . . . stabilizer of square matrix M
#
# If <g> given returns stabilizer in subgroup g of Symmetricgroup(Length(M))
# If vector <extra> given returns group stabilizing also <extra>
#
MatStab:=function(arg)local stab,M,g,r,l,I,p,s,k,n,i,j,e,inds;
if IsGroup(arg[1]) then g:=arg[1];arg:=arg{[2..Length(arg)]};fi;
M:=arg[1];k:=Length(M);
inds:=function(I)local iM;
iM:=List(I,function(i)local res;
res:=[Collected(M[i]{I}),Collected(M{I}[i]),M[i][i]];
if Length(arg)=2 then Add(res,arg[2][i]);fi;
return res;end);
return List(Set(iM),x->I{Filtered([1..Length(I)],j->iM[j]=x)});
end;
stab:=function(I)local ind,g,p,iM;
ind:=inds(I);
if Length(ind)>1 then return Concatenation(List(ind,J->stab(J)));
elif Length(I)>1 then
if Length(I)>7 then InfoChevie("#Large Block:",I,"\n");fi;
g:=MatStab(CoxeterGroupSymmetricGroup(Length(I)),M{I}{I});
p:=MappingPermListList([1..Length(I)],I);
return [rec(gens:=List(g.generators,x->x^p),ind:=I)];
else return [];
fi;
end;
if IsBound(g) then
for r in inds([1..k]) do g:=Stabilizer(g,Set(r),OnSets);od;
s:=Concatenation(List([1..k],i->List([1..k],j->[M[i][j],k*(i-1)+j])));
Sort(s);
l:=[];j:=0;
for i in [1..Length(s)] do
if i=1 or s[i][1]<>s[i-1][1] then j:=j+1;l[j]:=[];fi;
Add(l[j],s[i][2]);
od;
n:=Cartesian([1..k],[1..k]);
e:=Group(List(g.generators,p->PermListList(n,List(n,x->OnTuples(x,p)))),());
for s in l do e:=Stabilizer(e,s,OnSets);od;
return Group(List(e.generators,p->PermList(List([1..k],i->n[i^p][2]))),());
fi;
g:=Group(());I:=[];
for r in stab([1..k]) do
Append(I,r.ind);p:=MappingPermListList(I,[1..Length(I)]);
g:=Group(Concatenation(g.generators,List(r.gens,x->x^p)),());
g:=MatStab(g,M{I}{I});
od;
return Group(List(g.generators,x->x^(p^-1)),());
end;
#############################################################################
##
#F DistHelpedRepresentativeOperation(G,x,y,opr,dist) . . Heuristic version of
## 'RepresentativeOperation' when a distance is known between x and y.
#
# In the case where x and y live in a space on which we have a distance,
# we can try to get from x to y in G by multiplying by the generator
# which brings us closer. Coupled with the trick of applying a random
# perturbation if we fall in a hole which is not our goal, this often
# works surprisingly well. 'dist(x,y)' should return rational numbers.
#
DistHelpedRepresentativeOperation:=function(G,x,y,opr,dist)local p,d,prev,cv,x1;
InfoChevie("# group:",Size(G)," too big - trying random walk\n");
# best generator in G towards x=y
cv:=function(x)local minimum,mp,i,nn;
mp:=();
minimum:=dist(x,y);
InfoChevie(minimum," \c");
for i in G.generators do
nn:=dist(opr(x,i),y);
if nn<minimum then minimum:=nn; mp:=i;fi;
od;
if mp<>() then InfoChevie(Position(G.generators,mp),"->",minimum," \c");fi;
return mp;
end;
p:=();
while true do
x1:=opr(x,p);
prev:=dist(x1,y);
if prev=0 then InfoChevie("\n"); return p;fi;
p:=p*cv(x1);
d:=dist(opr(x,p),y);
if d=prev then
# Print(Format(opr(x,p)-y),"\n");
InfoChevie("\n#stalled -- restarting at a random element of G\n");
p:=p*Random(G);
fi;
od;
end;
#############################################################################
##
# PermMatMat(M, N[, extra1 ,extra2]) . . . . p such that OnMatrices(M,p)=N
#
# An analogue of PermListList for matrices. If in addition the vectors
# extra1 and extra2 are given, find p such that Permuted(extra1,p)=extra2
#
PermMatMat:=function(arg)local ind,l,I,J,r,p,e,g,n,h,sg,M,N;
M:=arg[1];N:=arg[2];
sg:=n->Group(Concatenation(List([1..n-1],i->List([i+1..n],j->(i,j)))),());
ind:=function(I,J)local iM,iN,p;
iM:=List(I,function(i)local res;
res:=[Collected(M[i]{I}),Collected(M{I}[i]),M[i][i]];
if Length(arg)>2 then Add(res,arg[3][i]);fi;
return res;end);
iN:=List(J,function(i)local res;
res:=[Collected(N[i]{J}),Collected(N{J}[i]),N[i][i]];
if Length(arg)>2 then Add(res,arg[4][i]);fi;
return res;end);
if Set(iM)<>Set(iN) then return false;fi;
iM:=List(Set(iM),x->I{Filtered([1..Length(I)],j->iM[j]=x)});
iN:=List(Set(iN),x->J{Filtered([1..Length(J)],j->iN[j]=x)});
if List(iM,Length)<>List(iN,Length) then return false;fi;
if Length(iM)=1 then
if Length(I)>8 then InfoChevie("large block:",Length(I),"\n");
p:=DistHelpedRepresentativeOperation(sg(Length(I)),M{I}{I},N{J}{J},
OnMatrices,function(M,N)return Sum(M-N,x->Number(x,y->y<>0));end);
else p:=RepresentativeOperation(sg(Length(I)),M{I}{I},N{J}{J},OnMatrices);
fi;
if p=false then return false;fi;
return [[Permuted(I,p),J]];
else p:=Zip(iM,iN,ind);
if false in p then return false;
else return Concatenation(p);
fi;
fi;
end;
l:=ind([1..Length(M)],[1..Length(N)]);
if l=false then return false;fi;
I:=[];J:=[];g:=Group(());
for r in l do
Append(I,r[1]);Append(J,r[2]);
n:=Length(r[1]);
p:=MappingPermListList([1..n],[1..n]+Length(I)-n);
h:=MatStab(M{r[1]}{r[1]});
h:=Group(List(h.generators,x->x^p),());
g:=Group(Concatenation(g.generators,h.generators),());
if n>1 then
InfoChevie("# I=",Length(I)," newstab=",Size(h)," all=",Size(g),"\c");fi;
e:=RepresentativeOperation(g,M{I}{I},N{J}{J},OnMatrices);
if e=false then return false;
else I:=Permuted(I,e);
fi;
g:=MatStab(g,M{I}{I});
if n>1 then InfoChevie(" stab=",Size(g),"\n");fi;
od;
return MappingPermListList(I,J);
end;
------------------------------------------------------------------------
Jean MICHEL, Equipe des groupes finis, Institut de Mathematiques UMR7586
Bureau 9D17 tel.(33)144278119, 175, rue du Chevaleret 75013 Paris
More information about the Forum
mailing list