5
$\begingroup$

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
$\endgroup$
5
  • $\begingroup$ What version do you run on and what is a typical time you got for some deg and number of coordinates? $\endgroup$ – Alex Trounev Jan 27 at 11:04
  • $\begingroup$ The action of a given group element on monomials of a given degree can be represented by a square matrix of dimension equal to the monomial count (as you indicate). Call this matrix G. Then you have G.H.ConjugateTranspose[G]==H. Since the group action is invertible, you can get linear equations in the variables of H as H.ConjugateTranspose[G]==Inverse[G].H. If your current code not setting this up as a set of simultaneous linear systems (one per group generator), that would be an avenue to explore. $\endgroup$ – Daniel Lichtblau Jan 27 at 15:15
  • $\begingroup$ @AlexTrounev I've updated the code so that it runs for varying numbers of coordinates and added some timings in the post. I also renamed the degree of the monomials to k and left deg as a parameter in the Zn group elements (these two are not the same in general). $\endgroup$ – ala10 Jan 27 at 17:23
  • $\begingroup$ @DanielLichtblau I think that is what I am doing at the moment though I might not have set it up in the most straightforward way. The line H$inv=H/.(Solve... solves for the components of H as a simultaneous system. I think the bigger problem is finding the matrix G that acts on the monomials given the action on the coordinates - it seems my code spends the most time in this bit. At the moment I'm taking $g\cdot\vec z$, converting it to a list of a substitution rules, using those rules on the monomials and then using SolveAlways to find the G that would have given the transformation. $\endgroup$ – ala10 Jan 27 at 17:27
  • 1
    $\begingroup$ (1) Also you can add conditions to the effect that H is Hermitian. Separating real and imaginary parts in the equations gives rise to more equations but they might be simpler. (2) The solver code is a bit complicated for my taste. But I suspect the bottleneck might be in doing exact linear algebra. Might want to experiment with this at machine precision to see if that is in fact where it becomes slow. $\endgroup$ – Daniel Lichtblau Jan 27 at 23:12
5
+100
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Amazing! This is a huge speed-up! The SolveAlways that I was using to find the form of the matrices took 120s for the $d=3$, $k=10$ example, whereas it takes 0.17s with your method! Now I'm wondering if there is a way to speed-up solving for the allowed components of H? (the H$inv=H/.(Solve... part). This step alone takes 7s for the same example and represents 97% of the total running time! I guess I need to go back to your comment and think about how to best represent this step using linear algebra instead of just a brute-force Solve. $\endgroup$ – ala10 Jan 27 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.