Visualizing diagrams needed to compute $\operatorname{Tr}(A^3 (A^T)^3)$

Question

I'm looking for help getting Mathematica code to construct diagrammatic expressions like the following, obtained by River Li as a way to compute $\operatorname{Tr}(A^2 (A^T)^2)$ for $d\times d$ matrix $A$ with IID Gaussian entries.

The clearest explanation of such diagrams is given in Terry Tao's RMT book, Section 2.3.4.

He describes diagrammatic procedure for computing $E\operatorname{Tr}A^4$:

Which gives 4 different diagram types, that can be visualized method like above.

---

I'd like a method which takes an expression like $E \operatorname{Tr}(AAAA'A'A')$ where $A$ is a $d\times d$ matrix with IID standard normal entries, and generates a table of corresponding diagrams, like this, hand-derived by River Li.

Background

For the last 4 years I've been slowly organizing Mathematica-related code that helps explain Gaussian expectations (starting with 2020 post cumulants vis , also post on wicks brackets, partition lattice). I'm estimating I'm about halfway done to a nice Wolfram Community post on combinatorics of random matrices.

Cross-posted to Wolfram Community, for more visibility from WRI employees

Answer 1

Explanation why some graphs appear more than ones:

It depends on the number of possible paths - some graphs have only one possibility, the graph below has four.

Functions definitions:

I borrowed removeIsomorphicDoublePaths from @Domen's answer.

removeIsomorphicPaths is variation of removeIsomorphicDoublePaths.

edgesTagged adds tags to edge for proper displaying in Graph.

ef is edge shape function.

removeIsomorphicDoublePaths[doublePaths_] := 
 Module[{doublePathsFlat, hashes, isoHashes}, 
  doublePathsFlat = Flatten /@ doublePaths;
  hashes = ToString@Values[PositionIndex[#]] & /@ doublePathsFlat;
  isoHashes = DeleteDuplicates[hashes];
  Part[doublePaths, First@FirstPosition[hashes, #] & /@ isoHashes]]
removeIsomorphicPaths[doublePaths_] := Module[{hashes, isoHashes},
  hashes = ToString@Values[PositionIndex[#]] & /@ doublePaths;
  isoHashes = DeleteDuplicates[hashes];
  Part[doublePaths, First@FirstPosition[hashes, #] & /@ isoHashes]]
nk[n_, k_] := n!/((n - k)!)
edgesTagged[e_] := Block[{po, v},
  v = {};
  If[(po = 
        Position[Append[Sort[Most@#], Last@#] & /@ v, 
         Append[Sort@Most[#], {#[[-1, 1]], _}]]) != {}, 
     AppendTo[v, Append[Most[#], Last@v[[(Last@po)[[1]]]] + {0, 1}]], 
     AppendTo[v, #]] & /@ 
   MapIndexed[Append[#, {1 + Boole[First@#2 > Length[e]/2], 1}] &, e];
  Join[Most@#, {Append[Last@#, Subtract @@ Most@#]}] & /@ v
  ]
arrow = Graphics[{Line[# + {-1, 0} & /@ {{-1, 1}, {1, 0}, {-1, -1}}]}];
ef = ({Thick, (RGBColor /@ {"#FF6B00", "#0094FF"})[[#2[[3, 
        1]]]], {Arrowheads[{{0.01, 0.5 + 0.01, arrow}}], 
     Arrow@BSplineCurve[{First@#, 
        If[#2[[3, 3]] != 0, 
         Mean@#[[{1, -1}]] + #2[[3, 3]]^2(#2[[3, 1]] - 3/2) {0, 
            0.4 + 0.2 #2[[3, 2]]/Abs[#2[[3, 3]]]}, 
         Splice[Table[Mean@#[[{1, -1}]], 
            3] + {{-0.05 - 0.05 #2[[3, 2]], 0.3} + {0, 
               0.2 #2[[3, 2]]/2}, {0, 0.5} + {0, 
               0.2 #2[[3, 2]]/2}, {0.05 + 0.05 #2[[3, 2]], 0.3} + {0, 
               0.1 #2[[3, 2]]}}Table[{1, (#2[[3, 1]] - 3/2)}, 3]]], 
        Last@#}]}} &);

Code to produce edges for all graphs.

do = 5;
tup = Prepend[#, 1] & /@ Tuples[Range[do], {do - 1}];
gr = GatherBy[#, Last] & /@ {removeIsomorphicPaths[tup], tup};
pairs = Flatten[Tuples /@ Transpose[gr], 1];
jp = Join @@@ Map[Partition[#, 2, 1] &, pairs, {2}];
edges = removeIsomorphicDoublePaths[
   Select[jp, AllTrue[Values@Counts[#], EvenQ] &]];
FullSimplify[
  Total[nk[n, Flatten[#] // Union // Length]*
      Times @@ ((Values@Counts[#] - 1)!!) & /@ edges]] // Apart

---

48 n + 28 n^2 + 28 n^3 + n^5

Code to display graphs:

Graph[DirectedEdge @@@ #, EdgeShapeFunction -> ef, 
    VertexStyle -> Black, VertexSize -> Tiny, 
    VertexCoordinates -> 
     Table[{n, 
       0}, {n, #[[All, {1, 2}]] // Flatten // Union // Length}], 
    ImageSize -> 
     200*(#[[All, {1, 2}]] // Flatten // Union // 
        Length)] & /@ (edgesTagged /@ edges) // Column

---

do=2

---

do=3

---

do=4

---

do=5

image 2

image 3

---

Some interesting graphs for do=6 out of all 559

And another two interesting graphs for do=6 where there is impossibility of having no crossings of edges.

Answer 2

Here is a very very stupid, naive and brute-force approach. Unfortunately, I don't know how to control the position of graph edges, so the graphs are not as pretty as the ones you have. SameStructureQ is by @lericr.

Because this is a brute-force approach, it gets slow and memory-intensive for larger $n$. I have managed to run it up to $n=5$, and the results match those of @Roman.

d = 4;
(* Check whether two lists are isomorphic *)
SameStructureQ[a_, b_] := Values[PositionIndex[a]] == Values[PositionIndex[b]]
generatePaths[n_, d_] := Prepend[1] /@ Tuples[Range[n], {d}]
makeEdges[lst_] := DirectedEdge @@@ Partition[lst, 2, 1]
(* Generate all possible non-isomorphic paths of length d on n$ vertices *)
n$ = d + 1;
paths = generatePaths[n$, d];
(* Generate all pairs of paths and select only those that end in the same vertex )
( For smaller memory consumption, use the following instead: )
( doublePaths = ResourceFunction["SelectTuples"][
                  paths,2,Last[#[[1]]]==Last[#[[2]]]&]; *)
doublePaths = Tuples[paths, 2];
(* Select only those that end in the same vertex *)
doublePaths = Cases[doublePaths, {{, last_}, {, last_}}];
(* Select only those paths for which all edges appear even many times *)
doublePaths = Select[doublePaths, AllTrue[Values@
      Counts[Partition[#[[1]], 2, 1]~Join~Partition[#[[2]], 2, 1]], EvenQ] &];
removeIsomorphicDoublePaths[doublePaths_] := 
  Module[{doublePathsFlat, hashes, isoHashes},
   doublePathsFlat = Flatten /@ doublePaths;
   hashes = ToString@Values[PositionIndex[#]] & /@ doublePathsFlat;
   isoHashes = DeleteDuplicates[hashes];
   Part[doublePaths, First@FirstPosition[hashes, #] & /@ isoHashes]
   ];
(* Select only non-isomorphic double paths *)
doublePathsUnique = removeIsomorphicDoublePaths[doublePaths];
(* Make edges from vertex paths *)
edges = Catenate /@ Map[makeEdges, doublePathsUnique, {2}];
(* Calculate E )
nk[n_, k_] := n!/((n - k)!)
FullSimplify[Total[nk[n, Max[Last /@ #]]Times @@ ((Values@Counts[#] - 1)!!) & /@
     edges]] // Apart
(* Plot graphs )
( Edges connecting the same vertices need to have different tags if 

   you want them to be coloured differently. Hence "". )
annotateEdges[edges_] := Module[{m = Length[edges], edges1, edges2},
  {edges1, edges2} = Partition[edges, m/2];
  MapThread[
    Annotation[Append[#1, #2], 
      EdgeStyle -> ColorData[97][1]] &, {edges1, Range[m/2]}]~Join~
   MapThread[
    Annotation[Append[#1, ToString[#2] <> "*"], 
      EdgeStyle -> Orange] &, {edges2, Range[m/2]}]
  ]
vc = Table[i -> {i, 0}, {i, n$}];
Grid[GatherBy[
  Graph[annotateEdges[#], VertexCoordinates -> vc, 
     VertexStyle -> Black, EdgeStyle -> Thick, 
     EdgeLabels -> "EdgeTag", ImageSize -> 120] & /@ edges, 
  VertexCount], Alignment -> Left]

$$E_d:=\mathbf{E}[ \operatorname{Tr}( A^d (A^{\top})^d)]$$

---

$d=1: \quad E_1 = n^2$

---

$d=2: \quad E_2 = n^3+2 n$

---

$d=3: \quad E_3 = n^4+10 n^2+4 n$

---

$d=4: \quad E_4 = n^5+28 n^3+28 n^2+48 n$

---

$d=5: \quad E_5 = n^6+62 n^4+110 n^3+468 n^2+304 n$

Large image with 559 diagrams ...

---

Edit explanation: In the first version of my answer, I generated tuples only from non-isomorphic paths. However, this misses some allowed combinations (e.g. these four were missing for $d=4$). Hence, I now generate tuples from all possible paths, and remove isomorphic ones afterwards.

Answer 3

This just an extended comment.

When $d=3$ $A$ is a $n \times n$ array of independent unit Gaussian random variables, the trace of $A^d (A^T)^d$ is the sum of products of the form

$$A_{i_1,j_1}A_{i_2,j_2}A_{i_3,j_3}A_{i_4,j_4} A_{i_5,j_5}A_{i_{6},j_{6}}$$

and for general $d$:

$$A_{i_1,j_1}A_{i_2,j_2}\cdots A_{i_{2d},j_{2d}}$$

It seems you want to know the combinations that have positive expectations (as that is what Tao finds for powers of symmetric matrices). And finding out how often those combinations occur will apparently be addressed in a separate question.

From answers to one of your previous questions for any particular $d$ the multiplicity of each distinct element of $A$ has to be even for the expectation to be positive so one can determine the complete set of diagram types with the following:

diagramTypes[d_] := Select[IntegerPartitions[2 d], AllTrue[#, EvenQ] &]
diagramTypes[1]
(* {{2}} )
diagramTypes[2]
( {{4},{2,2}} )
diagramTypes[3]
( {{6},{4,2},{2,2,2}} )
diagramTypes[4]
( {{8},{6,2},{4,4},{4,2,2},{2,2,2,2}} )
diagramTypes[5]
( {{10},{8,2},{6,4},{6,2,2},{4,4,2},{4,2,2,2},{2,2,2,2,2}} *)