Finding invariant matrix given group elements

Question

I have a question about speeding up / optimising the following calculation. I have a feeling there is a way to rewrite it but I can't quite see what to do.

In $d=3$ dimensions, given a set of coordinates $\vec{z} = (z_0,z_1,z_2)$, I can construct a list of monomials of the coordinates at some degree $k$. For example at degree 2 I would have $\vec{s}_{(2)}=(z_0^2,z_0 z_1, z_1^2, \ldots)$. Given the action of a $3\times3$ matrix $g$ on the coordinates, I can compute the induced action on the vector of monomials: $$ \vec{z} \to g\cdot\vec{z} \qquad \implies \qquad \vec{s}_{(2)} \to \tilde{g}\cdot \vec{s}_{(2)}.$$

I can also define an inner product on the list of monomials via $$ s_{(2)}^T H s_{(2)}^* $$ where $H$ is some hermitian matrix.

My question is how to solve for the allowed form of $H$ such that the inner product is invariant under the induced action of $g$ (or a list of such matrices). In other words, what is the allowed form of $H$ such that the resulting polynomial $ s_{(2)}^T H s_{(2)}^* $ will be invariant under $\vec{z} \to g\cdot\vec{z}$? In terms of the matrix acting on the vector of monomials, the condition on $H$ should be $$H = \tilde{g}^T H \tilde{g}^* .$$

In the example code that follows, I take $g$ to permute the coordinates or multiply by phases, but it can be more complicated in general.

The following is some code that (seems to) solves this problem, however it becomes extremely slow when the degree of the monomials (and the size of the $H$ matrix) increases. Also, I'm interested in cases where there are more than 3 coordinates. The size of the monomial basis then grows very quickly, as does, for example, the number of elements of the permutation group. Thanks to this, I've been trying to come up with a better way to find the invariant form of $H$. It seems most of the time is spent in SolveAlways.

Any suggestions for optimising this would be greatly appreciated! Thanks in advance!

Edit: I've updated the code so you can also change the number of coordinates by changing d at the start.

Also, here are some timings on MMA 12.1. For just the Zn set of matrices, at $d=3$ and $k=2,4,6,8,10$ the timings are 0.2s, 1.3s, 7.5s, 37s, 180s. Ideally I'd like to be able to do this for $d=5$ and $k=10$ in a reasonable amount of time.

(* define vectors of coordinates and conjugates *)
d=3;
zv=Array[Subscript[z,#]&,d,0]
zbv=Array[Subscript[zb,#]&,d,0];

(* function to generate list of monomials at degree m *)
monomialList[vars_List,m_Integer?NonNegative]:=Flatten[Map[Apply[Times,vars^#]&,Table[FrobeniusSolve[ConstantArray[1,Length[vars]],k],{k,m,m}],{2}]]

(* pick degree of monomials *)
k=2;

(* output is list of monomials at degree k in Subscript[z,i] *)
s=monomialList[zv,k]
sb=monomialList[zbv,k];


(* Define action of roots of unity by matrices *)
(* imagining these are elements that leave invariant the equation \!\(
\*SubscriptBox[\(\[Sum]\), \(i\)]
\*SuperscriptBox[
SubscriptBox[\(z\), \(i\)], \(deg\)]\)=0 *)
deg=4;
Zn=Flatten[Table[DiagonalMatrix[ReplacePart[ConstantArray[1,Length[zv]],l->Exp[2I \[Pi] k/deg]]],{l,1,Length[zv]},{k,0,deg-1}],{{1,2},{3},{4}}];

(* Define action of symmetric group by matrices *)
pp=Permutations[zv];
ress=Table[Normal[SparseArray[#->1&/@Table[{j,Position[pp[[k]],zv[[j]]][[1,1]]},{j,1,Length[zv]}]]],{k,1,Length[pp]}];
Sn=Table[ress[[kk]],{kk,1,Length[ress]}];


(* function that takes a list of coordinates,list of monomials and a list of matrices that act on the coordinates *)
(* output is a matrix H which defines an inner product invariant under the action of the matrices *)
computeInvariantH[z$list_,mono$list_,g$list_]:=Module[{g$subs,num,A,H,H$inv},

(* convert action of group elements on Subscript[z, i] to a list of substitution rules *)
g$subs=Thread[z$list->#]&/@(g$list.z$list);

(* how many monomials? *)
num=mono$list//Length;

(* define generic matrices with correct dimensions to act on list of monomials *)
(* H will be the invariant matrix in the end *)
(* A will be the matrix that describes the group action on the monomial basis *)
A=Array[a,{num,num}];
H=Array[h,{num,num}];

(* convert the list of group elements which act on Subscript[z,i] to a list of num x num matrices that act on monomial basis *)
(* we known how a group element acts on coordinates and converted this to a substitution rule above *)
(* we can then use this rule to see how it acts on the monomials *)
(* we then solve for the form of the matrix A that would give the same action on the monomials *)
A$list=Table[A/.(SolveAlways[A.mono$list==(mono$list/.g$subs[[n]]),zv]//First),{n,1,g$subs//Length}];

(* compute the allowed form of H such that s^T.H.s^* is invariant under s\[Rule]A.s *)
(* implies that H must satisfy H=A^T.H.A^* *)
(* we get one of these equations for each group element and then solve then simultaneously for the components of H *)
H$inv=H/.(Solve[Outer[Times,ConstantArray[1,g$subs//Length],H]==MapThread[Dot,{Transpose[A$list,2<->3].H,Conjugate[A$list]}]]//First//FullSimplify);

(* define real and imaginary parts of matrix *)
H$inv$re=H$inv/.{h->r};
H$inv$im=H$inv/.{h->i};

(* impose real is symmetric and imaginary is antisymmetric *)
re$subs=Solve[H$inv$re==Transpose[H$inv$re]]//First;
im$subs=Solve[H$inv$im==-Transpose[H$inv$im]]//First;

(* sub in constraints *)
H$inv$re=H$inv$re/.re$subs;
H$inv$im=H$inv$im/.im$subs;

(* form of invariant hermitian matrix *)
H$inv=H$inv$re+I H$inv$im;

H$inv]


(* compute invariant H for Zn permutations *)
H$inv$Zn=computeInvariantH[zv,s,Zn]

(* explicitly check form of inner product left invariant by all Zn elements *)
MapThread[Dot,{(A$list.s).H$inv$Zn,(Conjugate[A$list].sb)}]==ConstantArray[1,Zn//Length] s.H$inv$Zn.sb//Simplify


(* compute invariant H for Sn permutations *)
H$inv$Sn=computeInvariantH[zv,s,Sn]

(* explicitly check form of inner product left invariant by all Sn elements *)
MapThread[Dot,{(A$list.s).H$inv$Sn,(Conjugate[A$list].sb)}]==ConstantArray[1,Sn//Length] s.H$inv$Sn.sb//Simplify

Edit: The other thing I've been thinking about is whether one can construct the general form of the invariant polynomial that would be given by $$ s_{(2)}^T H s_{(2)}^* $$ directly, possibly using something like this but for more than one variable. I could then read off $H$ by computing the derivatives with respect to $s$ and $s^*$. An issue with this that jumps to mind is that I would usually compute the invariant polynomials by acting with the Reynold's operator which requires an average over all of the group elements (and not just the generators that I give in the code). One can compute the group elements using the answer from here. However, the number of group elements grows very quickly and this seems like the wrong way to do it. Anyway, some code to do this is

(* function to symmetrize a polynomial over a set of group elements *)
symmetrize$poly[f_,m_]:=With[{n=Length[m]},1/n Simplify@Total@Map[f/.Thread[zv-> (#.zv)]/.Thread[zbv-> (Conjugate[#].zbv)]&,m]];

(* all monomials constructed from s and s^* *)
monos=Outer[Times,s,sb]//Flatten;

(* function to find all elements of group *)
(* above are only the generators of the group *)
(* Reynold's operator is an average over all group elements, not the generators *)
groupElements[G_]:=Module[{count=0},NestWhile[(Print[count++];
Union[#~Join~Flatten[Outer[Dot,Join[G],#,1],1]])&,{IdentityMatrix[d]},Length[#2]!=Length[#1]&,2,99]]

(* only averaging over all of the group elements gives polynomials that are honestly invariant *)
inv$monos=DeleteCases[symmetrize$poly[monos,groupElements[Zn]]//DeleteDuplicates,0]

(* find form of H matrix *)
coeffs=Array[b,inv$monos//Length];
H$inv$poly=D[D[coeffs.inv$monos/.Thread[s->Array[a[#]&,s//Length]]/.Thread[sb->Array[ab[#]&,s//Length]],{Array[a[#]&,s//Length]}],{Array[ab[#]&,s//Length]}];

(* define real and imaginary parts of matrix *)
H$inv$poly$re=H$inv$poly/.{b->br};
H$inv$poly$im=H$inv$poly/.{b->bi};

(* impose real is symmetric and imaginary is antisymmetric *)
re$subs$poly=Solve[H$inv$poly$re==Transpose[H$inv$poly$re]]//First;
im$subs$poly=Solve[H$inv$poly$im==-Transpose[H$inv$poly$im]]//First;

(* sub in constraints *)
H$inv$poly$re=H$inv$poly$re/.re$subs$poly;
H$inv$poly$im=H$inv$poly$im/.im$subs$poly;

H$inv$poly=H$inv$poly$re+I H$inv$poly$im;

(* Check this matches the form given by the previous method *)
Length[Variables[H$inv$poly]]==Length[Variables[H$inv$Zn]]
Solve[H$inv$poly==H$inv$Zn]

Edit: Following the answer by Daniel, I've updated the code as below. Repeating the timings above we have for just the Zn set of matrices, at $d=3$ and $k=2,4,6,8,10$: 0.03s, 0.2s, 0.65s, 2.3s, 7.5s. That's a speed-up of 24x! A huge improvement. Now most of the time is spent solving for the actual form of H rather than coming up with the matrices for the generators acting on the monomials. I can't think of a better way than just solving the linear equations for H with Solve.

(* define vectors of coordinates and conjugates *)
d=3;
zv=Array[Subscript[z,#]&,d,0]
zbv=Array[Subscript[zb,#]&,d,0];

(* function to generate list of monomials at degree m *)
monomialList[vars_List,m_Integer?NonNegative]:=Flatten[Map[Apply[Times,vars^#]&,Table[FrobeniusSolve[ConstantArray[1,Length[vars]],k],{k,m,m}],{2}]]

(* pick degree of monomials *)
k=2;

(* output is list of monomials at degree k in Subscript[z,i] *)
s=monomialList[zv,k]
sb=monomialList[zbv,k];


(* Define action of roots of unity by matrices *)
(* imagining these are elements that leave invariant the equation Total[zv^deg] *)
deg=4;
Zn=Flatten[Table[DiagonalMatrix[ReplacePart[ConstantArray[1,Length[zv]],l->Exp[2I \[Pi] k/deg]]],{l,1,Length[zv]},{k,0,deg-1}],{{1,2},{3},{4}}];

(* Define action of symmetric group by matrices *)
pp=Permutations[zv];
ress=Table[Normal[SparseArray[#->1&/@Table[{j,Position[pp[[k]],zv[[j]]][[1,1]]},{j,1,Length[zv]}]]],{k,1,Length[pp]}];
Sn=Table[ress[[kk]],{kk,1,Length[ress]}];

(* functions suggested by Daniel *)
action[g_,vars_,monoms_List]:=Map[action[g,vars,#]&,monoms]
action[g_,vars_,monom_]:=monom/.Thread[vars->g.vars]

actionVector[g_,vars_,monoms_,monom_]:=Module[{avec=action[g,vars,monom]},Map[Coefficient[avec,#]&,monoms]]
actionMatrix[g_,vars_,monoms_List]:=Map[actionVector[g,vars,monoms,#]&,monoms]



(* function that takes a list of coordinates,list of monomials and a list of matrices that act on the coordinates *)
(* output is a matrix H which defines an inner product invariant under the action of the matrices *)
computeInvariantH[z$list_,mono$list_,g$list_]:=Module[{g$list$nodups,num,H,H$inv,H$inv$re,H$inv$im,re$subs,im$subs,h},

    (* remove duplicates from list of group generators *)
    g$list$nodups=g$list//DeleteDuplicates;

(* how many monomials? *)
num=mono$list//Length;

(* define generic matrices with correct dimensions to act on list of monomials *)
(* H will be the invariant matrix in the end *)
(* A will be the matrix that describes the group action on the monomial basis *)
H=Array[h,{num,num}];

(* convert the list of group elements which act on Subscript[z,i] to a list of num x num matrices that act on monomial basis *)
(* we known how a group element acts on coordinates and converted this to a substitution rule above *)
(* we can then use this rule to see how it acts on the monomials *)
(* we then solve for the form of the matrix A that would give the same action on the monomials *)
A$list=actionMatrix[#,zv,s]&/@g$list$nodups;

(* compute the allowed form of H such that s^T.H.s^* is invariant under s\[Rule]A.s *)
(* implies that H must satisfy H=A^T.H.A^* *)
(* we get one of these equations for each group element and then solve then simultaneously for the components of H *)
H$inv=H/.(Solve[Transpose[Inverse/@A$list,2<->3].H==Transpose[(ConjugateTranspose/@A$list).Transpose[H],2<->3]]//First//Simplify);

(* define real and imaginary parts of matrix *)
H$inv$re=H$inv/.{h->r};
H$inv$im=H$inv/.{h->i};

(* impose real is symmetric and imaginary is antisymmetric *)
re$subs=Solve[H$inv$re==Transpose[H$inv$re]]//First;
im$subs=Solve[H$inv$im==-Transpose[H$inv$im]]//First;

(* sub in constraints *)
H$inv$re=H$inv$re/.re$subs;
H$inv$im=H$inv$im/.im$subs;

(* form of invariant hermitian matrix *)
H$inv=H$inv$re+I H$inv$im;

H$inv]

(* compute invariant H for Zn permutations *)
H$inv$Zn=computeInvariantH[zv,s,Zn];//AbsoluteTiming

(* explicitly check form of inner product left invariant by all Zn elements *)
MapThread[Dot,{(A$list.s).H$inv$Zn,(Conjugate[A$list].sb)}]==ConstantArray[1,Zn// DeleteDuplicates//Length] s.H$inv$Zn.sb//Simplify


(* compute invariant H for Sn permutations *)
H$inv$Sn=computeInvariantH[zv,s,Sn];//AbsoluteTiming

(* explicitly check form of inner product left invariant by all Sn elements *)
MapThread[Dot,{(A$list.s).H$inv$Sn,(Conjugate[A$list].sb)}]==ConstantArray[1,Sn// DeleteDuplicates//Length] s.H$inv$Sn.sb//Simplify

Answer 1

This is not intended as a complete answer but might get you past some sticky points in converting an action to a matrix representation.

First set up your variables and group actions.

zv = Array[z, d, 0];

deg = 4;
Zn = Flatten[
   Table[DiagonalMatrix[
     ReplacePart[ConstantArray[1, Length[zv]], 
      l -> Exp[2 I \[Pi] k/deg]]], {l, 1, Length[zv]}, {k, 0, 
     deg - 1}], {{1, 2}, {3}, {4}}];

I define the monomial producer a bit differently.

monomialList[vars_, deg_] := 
 Union@Flatten@Outer[Times, Apply[Sequence, ConstantArray[vars, deg]]]

monoms2 = monomialList[zv, 2]

(* Out[514]= {z[2]^2, z[1] z[2], z[1]^2, z[0] z[2], z[0] z[1], z[0]^2} *)

Now we create a way to perform the group action on a monomial or list thereof.

action[g_, vars_, monoms_List] := Map[action[g, vars, #] &, monoms]
action[g_, vars_, monom_] := monom /. Thread[vars -> g . vars]

Example:

In[516]:= action[Zn[[4]], zv, monoms2]

(* Out[516]= {z[2]^2, z[1] z[2], 
 z[1]^2, -I z[0] z[2], -I z[0] z[1], -z[0]^2} *)

Finally we convert an action to a matrix byt examining the effect on each monomial (which gives a vector).

actionVector[g_, vars_, monoms_, monom_] := Module[
  {avec = action[g, vars, monom]},
  Map[Coefficient[avec, #] &, monoms]
  ]
actionMatrix[g_, vars_, monoms_List] := 
 Map[actionVector[g, vars, monoms, #] &, monoms]

Again an example using one of the group elements.

In[521]:= actionMatrix[Zn[[4]], zv, monoms2]

(* Out[521]= {{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 
  0}, {0, 0, 0, -I, 0, 0}, {0, 0, 0, 0, -I, 0}, {0, 0, 0, 0, 0, -1}} *)

So this is a start.