RTNI - A symbolic integrator for Haar-random tensor networks

TL;DR

RTNI is a symbolic computation tool that calculates averages of tensor networks with Haar-random unitaries, aiding quantum information and holography research.

Contribution

It introduces a novel computer algebra package implementing graphical Weingarten calculus for tensor network averages with symbolic dimensions.

Findings

Successfully computes tensor network averages with Haar-random unitaries.

Supports symbolic tensor dimensions and complex tensor contractions.

Applied to entropy calculations in quantum information and holography models.

Abstract

We provide a computer algebra package called Random Tensor Network Integrator (RTNI). It allows to compute averages of tensor networks containing multiple Haar-distributed random unitary matrices and deterministic symbolic tensors. Such tensor networks are represented as multigraphs, with vertices corresponding to tensors or random unitaries and edges corresponding to tensor contractions. Input and output spaces of random unitaries may be subdivided into arbitrary tensor factors, with dimensions treated symbolically. The algorithm implements the graphical Weingarten calculus and produces a weighted sum of tensor networks representing the average over the unitary group. We illustrate the use of this algorithmic tool on some examples from quantum information theory, including entropy calculations for random tensor network states as considered in toy models for holographic duality. Mathematica and Python implementations are supplied.

Equations

$${\int }_{\mathcal{U}(d)}{U}_{i_1{j}_1}\cdots U_{i_p{j}_p}{\bar{U}}_{i_1^{\prime }{j}_1^{\prime }}\cdots {\bar{U}}_{i_p^{\prime }{j}_p^{\prime }}\mathrm{d}U=\sum _{\alpha ,\beta \in {\mathcal{S}}_p}{\delta }_{i_1{i}_{\alpha (1)}^{\prime }}\dots {\delta }_{i_p{i}_{\alpha (p)}^{\prime }}{\delta }_{j_1{j}_{\beta (1)}^{\prime }}\dots {\delta }_{j_p{j}_{\beta (p)}^{\prime }}{\mathrm{Wg}}_d({\alpha }^{-1}\beta ).$$

$${\int }_{\mathcal{U}(d)}{U}_{i_1{j}_1}\cdots U_{i_p{j}_p}{\bar{U}}_{i_1^{\prime }{j}_1^{\prime }}\cdots {\bar{U}}_{i_{p^{\prime }}^{\prime }{j}_{p^{\prime }}^{\prime }}\mathrm{d}U=0.$$

$$(f\ast g)(\pi )$$

$$h_d(\pi )$$

$${\mathrm{Wg}}_d:{\mathcal{S}}_p$$

$$\sigma$$

$$h_d\ast {\mathrm{Wg}}_d\ast h_d$$

$$\sum _{\tau \in {\mathcal{S}}_p}{\mathrm{Wg}}_d(\tau )d^{\#({\tau }^{-1}\pi )}$$

$${\mathrm{Wg}}_d(\sigma )=\frac{1}{(p!{)}^2}\sum _{\lambda ⊢p}\frac{({\chi }^{\lambda }(e){)}^2}{s_{\lambda ,d}(1)}{\chi }^{\lambda }(\sigma )$$

$$s_{\lambda ,d}(1)=\prod _{1\le i<j\le d}\frac{{\lambda }_i-{\lambda }_j+j-i}{i-j}.$$

$${\mathbb{E}}_U(\mathcal{D})=\sum _{\alpha ,\beta }{\mathcal{D}}_{\alpha ,\beta }{\mathrm{Wg}}_d({\alpha }^{-1}\beta ),$$

$$\displaystyle{\color[rgb]{0.0,0.2,0.13}e_{1}}$$

$$\displaystyle{\color[rgb]{0.0,0.2,0.13}e_{2}}$$

$$\displaystyle{\color[rgb]{0.0,0.2,0.13}e_{3}}$$

$$X_1{U}_1{X}_2{U}_3\cdots X_n{U}_n\text{ respectively its trace }\mathrm{Tr}(X_1{U}_1{X}_2{U}_3\cdots X_n{U}_n),$$

$$v=\mathbb{E}\mathrm{Tr}[XUY{U}^{\ast }]={\int }_{U\in \mathcal{U}(d)}XUY{U}^{\ast }\mathrm{d}U.$$

$$v=\frac{(\mathrm{Tr}X)(\mathrm{Tr}Y)}{d}.$$

$$w=\frac{\mathrm{Tr}Y}{d}X.$$

$$[\mathrm{id}\otimes \mathrm{Tr}](UA{U}^{\ast }),$$

$$\mathbb{E}[\mathrm{id}\otimes \mathrm{Tr}](UA{U}^{\ast })=I_n\mathrm{Tr}(A)k{\mathrm{Wg}}_{nk}((1))=\frac{\mathrm{Tr}A}{n}{I}_n.$$

$$\mathcal{T}(X):=\mathbb{E}(U\otimes U)X(U^{\ast }\otimes U^{\ast }),$$

$$Fx\otimes y=y\otimes x\text{ for }x,y\in {\mathbb{C}}^d.$$

$$\mathcal{T}(X)=\frac{1}{d^2-1}\left\{\left[\mathrm{Tr}X-\frac{\mathrm{Tr}(XF)}{d}\right]I+\left[\mathrm{Tr}(XF)-\frac{\mathrm{Tr}X}{d}\right]F\right\}.$$

$$\Phi :M_d(\mathbb{C})$$

$$X$$

$$f:=\mathrm{Tr}[{\omega }_k\cdot [\Phi \otimes \bar{\Phi }]({\omega }_d)]$$

$${\mathbb{C}}^d\otimes {\mathbb{C}}^d\ni {\Omega }_n:=\frac{1}{\sqrt{n}}\sum _{i=1}^n{e}_i\otimes e_i,$$

$$\underset{n\to \infty }{\lim }\mathbb{E}f=t+\frac{1-t}{k^2}.$$

$$({\mathbb{C}}^d{)}^{\otimes (1+{\mathrm{deg}}_x)}\ni |V_x⟩:=U_x|0⟩,$$

$${\Psi }_G:=\left[\underset{e\in E}{⨂}⟨{\Omega }_e|\right]\left[\underset{x\in V}{⨂}|V_x⟩\right],$$

$${\rho }_{G;A}:=[{\mathrm{id}}_A\otimes {\mathrm{Tr}}_{\bar{A}}]\left({\Psi }_G{\Psi }_G^{\ast }\right)\in {\mathcal{M}}_{d^{|A|}}(\mathbb{C}).$$

$${\hat{\rho }}_A:=\frac{{\rho }_A}{\mathrm{Tr}{\rho }_A},$$

Full text

1. 1 Introduction

2. 2 Our contribution

3. 3 Mathematical background: the Weingarten formula

4. 4 Graphical notation for tensors and for the Weingarten formula

5. 5 The RTNI package

6. 5.1 integrateHaarUnitary

7. 5.2 visualizeTN

8. 5.3 MultinomialexpectationvalueHaar

9. 5.4 Weingarten functions and precomputed values

10. 5.5 Differences between the Mathematica and Python packages

11. 6 Examples

12. 6.1 Twirling a matrix

13. 6.2 Several tensor factors

14. 6.3 Bi-partite twirling

15. 6.4 Bell state as input in conjugate random quantum channels

16. 6.5 Random tensor networks for holographic duality

17. 6.6 Example computation of expectation values of matrix products

18. A User interface of Python package

19. B Generating Weingarten functions

20. B.1 How the algorithm works

21. B.2 How to use the Python RTNI package for generating Weingarten functions

1. Introduction

The probabilistic method – made famous by Paul Erdös [Erd47] – now pervades all of discrete mathematics and computer science. In its most well-known form, it is used to demonstrate the existence of certain combinatorial objects by defining a suitable probability space, and showing that randomly chosen elements of that space exhibit the required property with non-zero probability. A prime example of this strategy is the proof of the achievability part of the channel coding theorem by Shannon [Sha48] in classical information theory, established roughly at the same time: here the method is applied to show the existence of capacity-achieving codes for a given communication channel. To mention just one additional famous example beyond the original applications to graph theory, the probabilistic method has been used to estimate sizes of $\varepsilon$-nets (respectively $\varepsilon$-samples) in spaces with limited Vapnik-Chervonenkis-Dimension [VC71].

While the probabilistic method constitutes a powerful mathematical tool, much more direct manifestations of randomness in the real world have been studied throughout the history of physics. Here we are particularly interested in the appearance of random matrices. The consideration of random matrices was initiated by Wigner in order to explain universal features of the eigenvalue spacing of atomic spectra [Wig55], but has since then gone undergone significant developments – see e.g., [GMGW98] for an in-depth review. Significant fields of application include the theory of disorder and associated transport phenomena in condensed matter systems, and the development of random matrix models in quantum field theory.

Random matrices, and more specifically the Haar distribution over a compact Lie group, also play a central role in quantum mechanics because of conservation laws. For example, any physical process respecting certain superselection rules can be described by a quantum channel (a completely positive trace-preserving map) which is covariant under the relevant group action. Thus such a channel can be seen as the result of “supertwirling”, i.e., the consequence of averaging (the inputs and outputs of) a given channel over the unitary group. Because this averaging operation typically leads to a drastic reduction of the number of free parameters, the set of covariant operations (and states) has often been a fruitful testbed and source of examples, e.g., in entanglement theory [Wer89]. Covariant channels (or instruments) also appear e.g., when considering certain operational figures of merit which are related for example to average fidelities (with averages taken over the Haar measure). In certain cases, the associated optimization problem can then be restricted to operations respecting the underlying symmetry. For example, optimal cloning devices where constructed following this approach in [Wer98, KW99]. We note that these are only some of the pioneering contributions to this area.

Interestingly, random matrices are fundamental to quantum information theory even beyond such symmetry considerations. Indeed, it was found that the probabilistic method can be used to establish existence (and it some cases “genericity”) of states and operations with particular properties; here the distributions are typically defined in terms of Haar-random unitaries. Pioneering work in this direction [HLSW04, HLW06] has shown for example the existence and ubiquity of states with a large amount of “locked” classical correlations. The probabilistic method was also used by Hastings [Has09b] to disprove the additivity conjecture for the Holevo-quantity: following a proof strategy successfully applied for Rényi entropies earlier [HW08, CHL*+*08] he demonstrated the existence of a pair of quantum channels for which the Holevo-quantity (which determines the classical capacity) is non-additive. Closer to the spirit of Shannon’s achievability theorem, the probabilistic method has also been used to show the achievability of the quantum capacity (more precisely, the coherent information): here the argument involves picking a random subspace respectively vectors (according to the Haar measure) of a Hilbert space [HHWY08, HSW08]. The corresponding “decoupling approach” is now a standard procedure in quantum information theory.

We refer to [CN16] for a detailed review of applications of random matrix techniques in quantum information theory. In most of these applications, averages are taken over the unitary group on some finite-dimensional space. Typically, one is interested in the moments of some expression (such as a matrix product) involving instances of random unitaries, and these are evaluated in an ad-hoc manner. In recent applications, random tensor network states are considered: here Haar-random unitaries are contracted together with fixed tensors in a network. An example is the consideration of random matrix product states [CGGPG13], and, more recently, the contraction of higher-dimensional random tensor networks in the context of holographic duality [HNQ*+*16] (see Example 6.5 below). In these cases, the evaluation of Haar averages of quantities of interest becomes highly non-trivial, either because higher moments are required, and/or because the average is over a $n$-tuple of independently and identically distributed unitaries.

2. Our contribution

Our work is motivated by the ubiquity of Haar-random averages in a large variety of settings, and the fact that the associated computations are – on a conceptual level – very similar. Indeed, the natural language capturing these problems is that of random tensor networks: these allow to capture a variety of mathematical objects/expressions involving random matrices, including matrix products, tracial expressions, and tensor contractions. Furthermore, expectation values of such random tensor networks are themselves linear combinations of tensor networks (which may represent scalars, matrices or tensors), and – for Haar-random unitaries – these expectation values can be computed algorithmically using the Weingarten calculus.

Even though detailed derivations of the Weingarten calculus can be found in the literature, there appears to be no general implementation relying on symbolic manipulation of tensor networks. Our Random Tensor Network Integrator (RTNI) package111Code available at https://github.com/MotohisaFukuda/RTNI aims to fill this gap by providing routines for computing expectation values of tensor networks involving random unitaries. This is envisioned to facilitate the study and application of objects involving Haar-random unitaries.

Let us briefly mention existing similar software packages dealing with Weingarten integration we are aware of. The IntU package [PM17] is a Mathematica package used for the computation of the averages monomials in the entries of a random unitary matrices, see also https://zksi.iitis.pl/wiki/projects:intu. There are two main differences with our implementation. Most importantly, we allow for symbolic Hilbert space dimensions, where as in IntU these dimensions must be fixed numerical integers. Moreover, our input data is a graph, as opposed to a monomial for IntU; we consider our input format more practical, especially for applications of the Weingarten calculus (and additionally provide subroutines for monomials). The IntHaar package [GK16] is a MAPLE implementation for the Weingarten calculus, see http://sites.math.rutgers.edu/~ag930/Maple%20Packages/IntHaar.txt for the code. Similarly to the IntU package, the program computes averages of polynomials in the entries of the random unitary matrix, but the Weingarten functions are computed symbolically. Moreover, the case of Haar-distributed orthogonal and symplectic matrices is implemented.

Outline

This document describes the RTNI package. In Section 3, we briefly review the underlying Weingarten formula. In Section 4, we recall the graphical Weingarten calculus which is realized by the RTNI package. Section 5 discusses the main routines provided by RTNI and their use. Finally, in Section 6, several explicit example computations are provided, illustrating the use of the package.

3. Mathematical background: the Weingarten formula

Here we briefly review the definition of some relevant mathematical objects, including, in particular, the Weingarten function. We defer the discussion of the graphical Weingarten calculus to Section 4 and instead begin with an “element-wise” expression of the Weingarten formula. This shows how to compute the average of monomials composed of matrix elements of random unitary matrices with respect to the Haar measure. This problem was first considered in the physics literature by Weingarten [Wei78] in the asymptotic limit of large dimension. An explicit expression for fixed matrix dimensions was obtained in [Col03] and [CŚ06] by Collins, respectively Collins and Śniady. The main insight exploited in those works is that the average of tensor product operators of the form $U^{\otimes n}\otimes\bar{U}^{\otimes n}$ is a weighted sum of permutation operators. This is a consequence of Schur-Weyl duality.

Theorem 3.1**.**

Let $d,p$ be positive integers and $i=(i_{1},\ldots,i_{p})$, $i^{\prime}=(i^{\prime}_{1},\ldots,i^{\prime}_{p})$, $j=(j_{1},\ldots,j_{p})$, $j^{\prime}=(j^{\prime}_{1},\ldots,j^{\prime}_{p})$ be $p$-tuples of positive integers from $\{1,2,\ldots,d\}$. Then

[TABLE]

If $p\neq p^{\prime}$ then

[TABLE]

Here we denoted by $\mathcal{U}(d)$ the unitary group acting on an $d$-dimensional Hilbert space, and the integrals are taken with respect to the normalized Haar measure on $\mathcal{U}(d)$. The function $\operatorname{Wg}_{d}$ is called the unitary Weingarten function, see Definition 3.2 below.

We note that software packages such as the IntU package [PM17] mentioned above primarily rely on expression (1). In contrast, our RTNI package exploits the graphical calculus introduced in the next section, including, in particular, expression (11) below.

The unitary Weingarten function $\operatorname{Wg}_{d}(\cdot)$ is a combinatorial object defined as follows. Let $\mathcal{A}(\mathcal{S}_{p})$ be the algebra of complex-valued functions on the symmetric group $\mathcal{S}_{p}$ with product given by the convolution

[TABLE]

The Dirac-Delta function $\delta_{e}(\pi)=\delta_{e,\pi}$ for $\pi\in\mathcal{S}_{p}$ is the identity element in $\mathcal{A}(\mathcal{S}_{p})$. For any $z\in\mathbb{C}$, let $h_{d}\in\mathcal{A}(\mathcal{S}_{p})$ be defined

[TABLE]

where $\#(\pi)$ denotes the number of cycles of the permutation $\pi$. Then $h_{d}$ belongs to the center $\mathcal{Z}(\mathcal{A}(\mathcal{S}_{p}))$.

Definition 3.2**.**

The unitary Weingarten function

[TABLE]

is a function of a dimension parameter $d\in\mathbb{C}$ and of a permutation $\sigma\in\mathcal{S}_{p}$. It is defined as the pseudo-inverse of the element $h_{d}\in\mathcal{A}(\mathcal{S}_{p})$ defined in (4), i.e., the unique element in $\mathcal{Z}(\mathcal{A}(\mathcal{S}_{p}))$ satisfying

[TABLE]

It can be shown that if $d\not\in\{0,\pm 1,\ldots,\pm(p-1)\}$, then $\operatorname{Wg}_{d}$ and $h_{d}$ are inverses of each other

[TABLE]

Furthermore, for integer $d\geq p$, the Weingarten functions can be written as

[TABLE]

Here, $\lambda\vdash p$ means that $\lambda$ is a partition of the integer $p$, $\chi^{\lambda}$ is the character of the irreducible representation of the symmetric group $\mathcal{S}_{p}$ specified by $\lambda$, and $s_{\lambda,d}(1)$ is the Schur polynomial evaluated at the identity (i.e., the dimension of the irreducible representation of $\mathcal{U}(d)$ with highest weight $\lambda$. The latter is given by the formula

[TABLE]

For integer $d<p$, the element $h_{d}\in\mathcal{A}(\mathcal{S}_{p})$ has no inverse, and formula (9) no longer applies, but a separate formula is needed. Indeed, expression (9) has poles when $d<p$. However, in computing averages over the Haar measure as in Theorem 3.1, those poles cancel, yielding a rational function in $d$ on the RHS of (1). This is why one can use the expression of Eq. (9) for the Weingarten function $\operatorname{Wg}_{d}$ for all $d\geq 1$ in the computations as in Theorem 3.1. We refer to [CŚ06] for more details.

4. Graphical notation for tensors and for the Weingarten formula

In this section, we briefly recall the graphical calculus method for computing unitary integrals introduced in [CN10], to which the readers are referred for more details. As discussed below, the RTNI package provides a general implementation of this method.

First, we explain the graphical formalism for tensor calculus, which can be traced back to Penrose [Pen71]; for modern presentations, see [CN10] or [Coe10]. The main object here is a complex vector space of diagrams. An element of this vector space will be represented as a list of pairs of the form $(g_{i},w_{i})$, where $g_{i}$ is a graph (a.k.a. a diagram or tensor network) and $w_{i}$ is a scalar weight; the corresponding element is $\mathcal{D}=\sum_{i}w_{i}\cdot g_{i}$. The diagrams (or graphs) correspond to tensors, and consist of boxes and wires. Boxes in a graph represent tensors. Each box has attached to it symbols of different shapes, where each shape corresponds to a vector space; identical shapes depict isomorphic vector spaces. Those symbols are filled (black) or empty (white), corresponding to primal or dual spaces. Wires connect these symbols, corresponding to tensor contractions $V\times V^{*}\to\mathbb{C}$. This is why each wire must connect symbols of the same shape, i.e. the same dimension. A diagram is a collection of such boxes and wires and corresponds to an element in a tensor product space. If there is more than one connected component in a graph, then the components are multiplied by taking the tensor product.

For example, the tensor (here, a bipartite matrix) $\operatorname{Tr}_{1}[X(Y\otimes Z)]\otimes W$, where $\operatorname{Tr}_{1}$ stands for the partial trace over the first space, is represented in Figure 1. Consider the primal side of the box for the matrix $X$. There are two symbols of different shapes implying that the effective primal space of $X$ is the tensor product of these two spaces.

Next, we describe how to efficiently compute expected values of diagrams which contain boxes of random unitary matrices. This technique in fact can be extracted from Theorem 3.1, once one uses the fact that the averaging process is linear. More precisely, the delta functions in each term in the RHS of (1) simply indicate how one needs to reconnect boxes corresponding to the random unitary operator $U$ and $\bar{U}$. For each pair of permutations $(\alpha,\beta)$ in (1), one eliminates $U$ and $\bar{U}$ boxes and reconnects the wires originally connected to these boxes to get a new diagram. Since these manipulations are made on wires and symbols (boxes), one can regard them as operations on a graph, identifying wires and symbols with edges and vertices, respectively. This process for a fixed pair of permutations is called a removal and the whole process which sums all the new graphs over all permutations is called the graph expansion. Formally, if $\mathcal{D}$ is a diagram containing boxes corresponding to a Haar-distributed random unitary matrix $U$, then its expectation value (with respect to $U$) is

[TABLE]

where $\mathcal{D}_{\alpha,\beta}$ is the diagram obtained from $\mathcal{D}$ by removing the boxes associated with $U$ and $\bar{U}$ and adding extra wires to the resulting diagram, as follows: we connect white decorations of the $i$-th $U$ box with the white decorations of the $\alpha(i)$-th $\bar{U}$ box and, by a similar procedure, the black decorations are paired using the $\beta$ permutation. One can see how this works in examples in Section 6, see in particular Figures 4 and 7.

5. The RTNI package

In this section, we explain how to use the RTNI package, available at https://github.com/MotohisaFukuda/RTNI. Readers are also referred to some examples in Section  6 while reading this section.

The following two functions are the main routines provided by RTNI:

1. (1)

$\textbf{integrateHaarUnitary}(TNList,\ varName,\ inDims,\ outDimns,\ totalDims)$: integrates out the random unitary matrix $varName$ in the tensor networks in $TNList$. 2. (2)

$\textbf{visualizeTN}(TNList,\{EdgeLabeling\rightarrow True\})$: produces a graphical output representing each tensor network (graph) in the list, together with the associated weights.

In Sections 5.1 and 5.2, we provide a detailed explanation of the functionality realized by these routines, as well as the syntax used for inputs and outputs. These two routines are the only routines needed for working with random tensor networks.

We also provide an additional auxiliary routine for matrix products, which covers many simple use-cases. It internally creates a tensor network representing the matrix product (respectively its trace), invokes integrateHaarUnitary, and re-expresses the resulting linear combination of tensor networks in terms of matrix products/tracial expressions.

1. (3)

$\textbf{MultinomialexpectationvalueHaar}(Dim,\ eList,\ variableList,\ usetrace)$: symbolically computes moments of matrix products involving Haar-random unitaries.

The corresponding syntax is explained in Section 5.3.

5.1. integrateHaarUnitary

An input of integrateHaarUnitary specifies the diagrams (tensor networks) containing random matrices to be averaged. At the program level, inputs of the function integrateHaarUnitary divide into five variables as indicated above. In the following, we describe each of them.

1. (1)

TNList is a list of pairs (tnsj,weightj) represeting a weighted linear combination of $n$ (with $n\in\mathbb{N}$ arbitrary) tensor networks, where the $j$-th tensor network tnsj has weight weightj. The syntax is as follows:

{ { tns1, weight1 }, { tns2, weight2 },…, { tnsn, weightn } }

For $n=1$ (i.e., if only one tensor network is considered) and the associated weight is equal to $1$, the alternative expression

tns

may be used. In these expressions, each term “tns” is a list of edges of the form

{ { vertex, vertex }, …, { vertex, vertex } }

and “weight” is a symbolic expression or scalar specifying the weight of the associated tensor network.

Vertices correspond to the different shaped and colored decorations attached to boxes in the graphical calculus from Section 4. Importantly, each term “vertex” is given by a list of four items in the format

{varName, ID, IN/OUT, leg}

where

• •

varName is a string denoting the name of the corresponding box; for boxes associated to Haar random unitary matrices, this string can end in the star symbol “$*$”, in which case this references the adjoint of the unitary operator.

• •

ID is an integer used to distinguish between several copies of the same box; this is in particular used to distinguish different copies of random unitary matrices and their adjoints.

• •

IN/OUT is a string, referring to whether the vertex corresponds to an input (”in”) or an output (”out”) of a matrix.

• •

leg is a positive integer corresponding to which tensor factor the vertex corresponds to.

We display in Figure 2 an example of a tensor network; the three edges appearing are encoded by (from left to right, and top to bottom; the colors are used for clarity only):

[TABLE]

2. (2)

varName is a letter showing over which random matrix the tensor network (respectively linear combination of tensor networks) is to be averaged. 3. (3)

inDims is a list specifying the dimensions of the factors in the tensor product corresponding to the outputs of the random unitary matrix. 4. (4)

outDims is a list specifying the dimensions of the factors in the tensor product corresponding to the inputs of the random unitary matrix. 5. (5)

totalDims is the size of the random unitary matrix to be integrated out.

It is important to note that the elements of the lists outDims, inDimns, as well as totalDims can be symbolic.

An output of integrateHaarUnitary is a list made of pairs of tensor networks and weights:

{ { tns1, weight1 }, { tns2, weight2 },…, { tnsm, weightm } }

Hence the output can be fed into integrateHaarUnitary for iterative calculations. Note that integrateHaarUnitary creates new edges (and new vertices): these correspond to symbols on the random matrix which are not connected to other symbols (which did not appear in TNList). In the Mathematica package the names of such new vertices starts with “dummy”, and in the Python package with “@”. One can also add such extra vertices manually.

5.2. visualizeTN

As the name indicates, this function creates a graphical representation of a single tensor network respectively a weighted list of tensor networks. An input of visualizeTN is the first argument of an input, or an output of integrateHaarUnitary. An example is shown in Fig. 3. The arguments are as follows:

1. (1)

TNList This is either a single tensor network, or a list of (weighted) tensor networks as in the description above. 2. (2)

{EdgeLabeling$\rightarrow$ True} This optional second argument forces the visualization to display edge labels, indicating which input/output is connected to which input/output. This argument can be omitted, in which case edges are not labeled.

5.3. MultinomialexpectationvalueHaar

An input of MultinomialexpectationvalueHaar specifies a “matrix product” expression of the form

[TABLE]

where each $U_{j}\in\{U,U^{*},U^{T},\overline{U}\}$ is either a fixed unitary $U$, its Hermitian adjoint $U^{*}$, its transpose $U^{T}$, or its entriwise complex conjugate $\overline{U}$. The function takes an additional dimension parameter $d$ and computes the Haar averaged expression symbolically. In more detail, the input variables are the following:

1. (1)

Dim is (a symbolical expression giving) the dimension of the matrices that are integrated over. 2. (2)

eList is a list of the form $\{\epsilon_{1},\ldots,\epsilon_{n}\}$ (for any integer $n$). Here $\epsilon_{j}\in\{1,2,3,4\}$ for each $j$. Here the values define the unitary part of the expressions (12) by the following translation table:

[TABLE] 3. (3)

variableList is a list of the form $\{X_{1},\ldots,X_{n}\}$. Here the $j$-term directly specifies the $j$-th variable $X_{j}$ in an expression of the form (12) 4. (4)

usetrace is either $\mathsf{True}$ or $\mathsf{False}$. It determines whether the first or second term in (12) is computed according to the table

[TABLE]

The output of “MultinomialexpectationvalueHaar” is a symbolic expression giving corresponding Haar average. Examples for this function are given in Section 6.6.

5.4. Weingarten functions and precomputed values

A table of precomputed Weingarten function is loaded upon program initialization. These are produced by a program in the Python package, and are stored in a subdirectory called precomputedWG. The latter computes all values of Weingarten functions based on expression (9). At present, all Weingarten functions with $p\leq 20$ are provided, with $d$ being a symbolic parameter. In particular, this means that the Mathematica package can only compute averages of tensor networks containing up to $p=20$ boxes corresponding to a random unitary matrix of some type (and the corresponding number of conjugate/adjoint boxes).

5.5. Differences between the Mathematica and Python packages

There are two implementations for Mathematica and Python, but they work almost identically from the user’s viewpoint, with the replacement of curly brackets (Mathematica) with square brackets (Python).

In the Mathematica package, it is possible to give as an input a single diagram, in which case the function assigns it the weight 1, as specified above. This feature is not implemented in Python.

In the Python package, the input of the routine integrateHaarUnitary is of the form (TNList, randomMatrixList), where

1. (1)

TNList is a (weighted) list of tensor networks as above, and 2. (2)

randomMatrixList is a list of the form

[ [varName, outDims, inDimns, totalDims], …, [varName, outDims, inDimns, totalDims] ]

The routine then computes the average of a tensor network over multiple random unitary matrices at once. In contrast, in Mathematica, the routine needs to be invoked several times.

The Mathematica package loads the (provided) precomputed values of the Weingarten functions. The Python routine first checks if there is a preexisting set of precomputed Weingarten function for $p$ as required. If not, this is computed and saved to the folder named “Weingarten”. To this end, the Python package also creates folders named “SGC” and “SP” which contain the characters of symmetric groups and Schur polynomials. See Appendix B for more details.

6. Examples

Here we present a collection of simple examples of use cases for the package RTNI. We gradually work our way up from basic situations to more involved cases, presenting the code (of the Mathematica implementation) and comparing with analytical computations.

6.1. Twirling a matrix

We start with the simplest possible example, which is to compute the overlap between a fixed matrix $X$ and a “twirled” matrix $UYU^{*}$, where $Y$ is a fixed $d\times d$ matrix and $U$ is a Haar-distributed random unitary matrix (also of size $d$). The result is a scalar

[TABLE]

The diagram corresponding to the integral above, and its Weingarten expansion are depicted in Figure 4; the result is

[TABLE]

In order to implement this in the RTNI package, we first need to input the corresponding tensor network (see Figure 4 top): it is a square, having 4 edges.

⬇

In[1]:= e1 = {{”U”, 1, ”out”, 1}, {”X”, 1, ”in”, 1}};

In[2]:= e2 = {{”Y”, 1, ”out”, 1}, {”U”, 1, ”in”, 1}};

In[3]:= e3 = {{”U*”, 1, ”out”, 1}, {”Y”, 1, ”in”, 1}};

In[4]:= e4 = {{”X”, 1, ”out”, 1}, {”U*”, 1, ”in”, 1}};

In[5]:= g = {e1, e2, e3, e4}

Out[1]= {{{U, 1, out, 1}, {X, 1, in, 1}}, {{Y, 1, out, 1}, {U, 1, in, 1}},

    {{U*, 1, out, 1}, {Y, 1, in, 1}}, {{X, 1, out, 1}, {U*, 1, in, 1}}}

The tensor network described by this code is visualized on the left in Figure 5, using the subroutine visualizeTN; it is essentially identical to Figure 4. The four edges correspond to the four matrix products in the trace. Let us analyze in detail the first edge, $e1$. Its source is the $X$-box having ID 1 (there is only one such box, i.e. $X$ occurs only once in the diagram), and more precisely it is an “input” of that box (hence the ”in” in the third position of the list describing the source). Since there is just one vector space in this example, the fourth parameter is 1. The sink of $e1$ is the output (the ”out” in the third position) of the first (the 1 in the second position) $U$ box; again, the last parameter is 1, indicating that we are using the first (and only) factor of the tensor product here. The tensor network $g$ is the collection of the four edges, and it is displayed using the function visualizeTN. Next, we perform the unitary integration, as follows:

⬇

In[1]:= Eg = integrateHaarUnitary[g, ”U”, {d}, {d}, d]

Out[1]= {{{{{Y, 1, out, 1}, {Y, 1, in, 1}}, {{X, 1, in, 1},

    {X, 1, out, 1}}}, 1/d}}

The resulting weighted sum of tensor networks (in this case having only one summand) is then displayed using the routine visualizeTN, see the right subfigure of Figure 5. There are two single vertex loop graphs corresponding to the two traces and the Weingarten weight.

The command integrateHaarUnitary invoked above is the main routine of our package, and it performs the integral of the tensor network $g$ over the unitaries corresponding to the letter $U$ (these are the boxes $U$ and $U^{*}$), assigning, in order, dimensions to the input spaces of $U$, the output spaces of $U$, and the total dimension of $U$ (here, there is just one input/output space, and all the dimensions are equal to $d$). The output $Eg$ is a (weighted) list of tensor networks (just one in this case), and we recognize in Figure 5 the tensor network with one vertex ($X$) with one loop attached (corresponding to $\operatorname{Tr}X$), the trace of $Y$, and the weight of the tensor network, $1/d$.

Next, we consider the same example but without the trace, i.e., the expression $w=XUYU^{*}$. Now, the result of the computation is no longer a scalar, but equal to

[TABLE]

To compute this, the following code can be used:

⬇

In[1]:= e1 = {{”U”, 1, ”out”, 1}, {”X”, 1, ”in”, 1}};

In[2]:= e2 = {{”Y”, 1, ”out”, 1}, {”U”, 1, ”in”, 1}};

In[3]:= e3 = {{”U*”, 1, ”out”, 1}, {”Y”, 1,”in”, 1}};

In[4]:= g = {e1, e2, e3}

Out[1]= {{{U, 1, out, 1}, {X, 1, in, 1}}, {{Y, 1, out, 1},

    {U, 1,in, 1}}, {{U*, 1, out, 1}, {Y, 1, in, 1}}}

In[5]:= Eg = integrateHaarUnitary[g, ”U”, {d}, {d}, d]

Out[2]= {{{{{Y, 1, out, 1}, {Y, 1, in, 1}},

    {{X, 1, in, 1}, {dummy-U*-IN-1-1, 1, out, 1}}}, 1/d}}

In this case, the RTNI package creates a “dummy” vertex for the input of $U^{*}$, in order to work with a proper tensor network. After the computation, this dummy vertex is connected to the output of $X$, see Figure 6.

We note that both expectation values $\mathbb{E}\operatorname{Tr}[XUYU^{*}]$ and $\mathbb{E}XUYU^{*}$ considered here could alternatively be computed using the subroutine MultinomialexpectationvalueHaar, see Example 6.6 below.

6.2. Several tensor factors

We now consider an example where the random unitary operators are bipartite: they act on a tensor product space. More precisely, we are interested in the partial trace of a twirled operator

[TABLE]

where the identity operator acts on a tensor factor (with label 1) of dimension $n$, while the trace operator acts on a different tensor factor (with label 2) of dimension $k$.

The analytical computation of the Weingarten integral with respect to $U$ is straightforward: the result reads

[TABLE]

The corresponding tensor networks for both the expresssion $[\operatorname{id}\otimes\operatorname{Tr}](UAU^{*})$ and its expectation are illustrated in Fig. 7.

The Mathematica code for the corresponding tensor network and for the computation of the expectation value with respect to the bipartite unitary matrix $U$ is as follows. In this example, one needs to specify that the unitary matrices $U$ and $U^{*}$ act on a tensor product of two vector spaces with respective dimensions $n$ and $k$. The input and output spaces are identical, and the total dimension is $nk$.

⬇

In[1]:= e1 = {{”A”, 1, ”out”, 1}, {”U”, 1, ”in”, 1}};

In[2]:= e2 = {{”A”, 1, ”out”, 2}, {”U”, 1, ”in”, 2}};

In[3]:= e3 = {{”U*”, 1, ”out”, 1}, {”A”, 1, ”in”, 1}};

In[4]:= e4 = {{”U*”, 1, ”out”, 2}, {”A”, 1, ”in”, 2}};

In[5]:= e5 = {{”U”, 1, ”out”, 2}, {”U*”, 1, ”in”, 2}};

In[6]:= g = {e1, e2, e3, e4, e5}

Out[1]= {{{A, 1, out, 1}, {U, 1, in, 1}}, {{A, 1, out, 2},

    {U, 1, in, 2}}, {{U*, 1, out, 1}, {A, 1, in, 1}},

    {{U*, 1, out, 2}, {A, 1, in, 2}}, {{U, 1, out, 2},

    {U*, 1, in, 2}}}

In[7]:= Eg = integrateHaarUnitary[g, ”U”, {n, k}, {n, k}, n k]

Out[2]= {{{{{A, 1, out, 1}, {A, 1, in, 1}}, {{A, 1, out, 2},

    {A, 1, in, 2}}, {{dummy-U-OUT-1-1, 1, in, 1},

    {dummy-U*-IN-1-1, 1, out, 1}}}, 1/n}}

The graphical representations of the input tensor network and of the result of the Weingarten integration are presented in Figure 8. One can easily check that this indeed reproduces the expression (15).

6.3. Bi-partite twirling

Bi- or multi-partite twirling is one of the most useful primitives in quantum information theory, when one has to do with some kind of permutational or rotational symmetry. Introduced by Werner in [Wer89], it has found many applications in the theory of quantum information. We shall be concerned here with the bi-partite case

[TABLE]

where $U\in\mathcal{U}(d)$ is a Haar-distributed random unitary matrix and $X\in\mathcal{M}_{d^{2}}(\mathbb{C})$ is a fixed matrix on which the twirling map $\mathcal{T}$ is acting. The Weingarten computation is more involved in this case, since we are dealing with two copies of the same random unitary matrix (we are in the case $p=2$ of formula (1)).

The analytical computation of the integral in (16) proceeds as follows. First, let $F:\mathbb{C}^{d}\otimes\mathbb{C}^{d}\to\mathbb{C}^{d}\otimes\mathbb{C}^{d}$ denote the flip (or swap) operator, defined by

[TABLE]

Then, the four choices of the permutations $\alpha,\beta\in\mathcal{S}_{2}$ from (1) give four terms, as follows:

[TABLE]

Putting everything together, we obtain

[TABLE]

The Mathematica code used to compute the twirling of an operator is below.

⬇

In[1]:= e1 = {{”X”, 1, ”out”, 1}, {”U”, 1, ”in”, 1}};

In[2]:= e2 = {{”U*”, 1, ”out”, 1}, {”X”, 1, ”in”, 1}};

In[3]:= e3 = {{”X”, 1, ”out”, 2}, {”U”, 2, ”in”, 1}};

In[4]:= e4 = {{”U*”, 2, ”out”, 1}, {”X”, 1, ”in”, 2}};

In[5]:= g = {e1, e2, e3, e4}

Out[1]= {{{X, 1, out, 1}, {U, 1, in, 1}}, {{U*, 1, out, 1},

    {X, 1, in, 1}}, {{X, 1, out, 2}, {U, 2, in, 1}},

    {{U*, 2, out, 1}, {X, 1, in, 2}}}

In[6]:= Eg = integrateHaarUnitary[g, ”U”, {d}, {d}, d]

Out[2]= {{{{{X,1,out,1},{X,1,in,1}},{{X,1,out,2},{X,1,in,2}},

    {{dummy-U-OUT-1-1,1,in,1},{dummy-U*-IN-1-1,1,out,1}},

    {{dummy-U-OUT-2-1,1,in,1},{dummy-U*-IN-2-1,1,out,1}}},1/(-1+d^2)},

    {{{{X,1,out,1},{X,1,in,1}},{{X,1,out,2},{X,1,in,2}},

    {{dummy-U-OUT-1-1,1,in,1},{dummy-U*-IN-2-1,1,out,1}},

    {{dummy-U-OUT-2-1,1,in,1},{dummy-U*-IN-1-1,1,out,1}}},1/(d-d^3)},

    {{{{X,1,out,1},{X,1,in,2}},{{X,1,out,2},{X,1,in,1}},

    {{dummy-U-OUT-1-1,1,in,1},{dummy-U*-IN-1-1,1,out,1}},

    {{dummy-U-OUT-2-1,1,in,1},{dummy-U*-IN-2-1,1,out,1}}},1/(d-d^3)},

    {{{{X,1,out,1},{X,1,in,2}},{{X,1,out,2},{X,1,in,1}},

    {{dummy-U-OUT-1-1,1,in,1},{dummy-U*-IN-2-1,1,out,1}},

    {{dummy-U-OUT-2-1,1,in,1},{dummy-U*-IN-1-1,1,out,1}}},1/(-1+d^2)}}

Note that the output list of (weighted) tensor networks consists of 4 tensor networks, corresponding, respectively, to the entries $(1,1)$, $(2,1)$, $(1,2)$, $(2,2)$ of the table above. The input and output are illustrated in Figure 9.

6.4. Bell state as input in conjugate random quantum channels

In this section, we discuss an application of RTNI to the theory of random quantum channels. Consider a random quantum channel given by a random isometry $V:\mathbb{C}^{d}\to\mathbb{C}^{k}\otimes\mathbb{C}^{n}$,

[TABLE]

Such random quantum channels were first considered in [HW08] in relation to the additivity conjecture in quantum information theory, solved in the negative by Hastings [Has09a]. In such considerations, one bounds the minimum output entropy of the tensor product channel $\Phi\otimes\bar{\Phi}$ by using the overlap

[TABLE]

between the output of the channels acting on a maximally entangled state $\omega_{d}$ with another maximally entangled state (on the output spaces) $\omega_{k}$. Here $\bar{\Phi}$ is the complex conjugate channel. The diagram for the scalar $f$ is given in Figure 10. Recall that the maximally entangled state is defined, in general, by $\omega_{n}:=\Omega_{n}\Omega_{n}^{*}$, with

[TABLE]

where $\{e_{1},\ldots,e_{n}\}$ is a given orthonormal basis of $\mathbb{C}^{n}$. In [HW08], it was shown that, for all (random) isometries $V$, $f\geq d/(kn)$. This bound was improved in [CN10], where it was shown that, in the asymptotic regime where $k$ is fixed, $n\to\infty$, and $d\sim tkn$ for some $t\in(0,1)$,

[TABLE]

We confirm this computation using our symbolic Weingarten integration routines; note that the corresponding edges are also represented in Figure 10.

⬇

In[1]:= e1 = {{”U*”, 2, ”out”, 1}, {”U”, 1, ”in”, 1}};

In[2]:= e2 = {{”U*”, 1, ”out”, 1}, {”U”, 2, ”in”, 1}};

In[3]:= e3 = {{”U”, 1, ”out”, 1}, {”U*”, 1, ”in”, 1}};

In[4]:= e4 = {{”U”, 2, ”out”, 1}, {”U*”, 2, ”in”, 1}};

In[5]:= e5 = {{”U”, 1, ”out”, 2}, {”U*”, 2, ”in”, 2}};

In[6]:= e6 = {{”U”, 2, ”out”, 2}, {”U*”, 1, ”in”, 2}};

In[7]:= g = {e1, e2, e3, e4, e5, e6};

In[8]:= listg = {{g, 1/(d k)}};

In[9]:= Eg = integrateHaarUnitary[listg, ”U”, {d}, {n, k}, n k]

Out[1]= {{{}, -((d^2 k^2 n)/(d k - d k^3 n^2)) - (d k n^2)/(

d k - d k^3 n^2) + (d k^2 n)/(d k^2 n - d k^4 n^3) + (d^2 k n^2)/(

d k^2 n - d k^4 n^3)}}

In[10]:= overlap = Eg[[1, 2]];

In[11]:= overlap = overlap /. {d -> t n k};

In[12]:= Assuming[t > 0 && k > 1, Limit[overlap, n -> Infinity]]

Out[2]= (1 - t)/k^2 + t

Note that we obtain first a formula for $\mathbb{E}f$ at finite $n,d$; we then take the limit $n\to\infty$ to recover (19).

6.5. Random tensor networks for holographic duality

In this section, we illustrate the versatility of the RTNI package on an example taken from [HNQ*+*16]. The authors of [HNQ*+*16] relate the entanglement properties of a model of random tensors to the AdS/CFT correspondence. In our treatment below, we shall consider a simplified version of their model of random tensors (to be precise, we shall not consider bulk states as in [HNQ*+*16, Section 2]) that will be analyzed on some example with the help of the RTNI routines.

To a given unoriented simple graph $G=(V,E)$ on $n=|V|$ vertices, we shall associate an $n$-tensor $\Psi_{G}\in(\mathbb{C}^{d})^{\otimes n}$, as follows. For a vertex $x\in V$ of degree $\mathrm{deg}_{x}$, consider a random pure state (unit vector)

[TABLE]

where $U_{x}$ is a Haar-distributed random unitary matrix and $|0\rangle$ is a fixed vector; moreover, assume that the family $\{U_{x}\}_{x\in V}$ is independent. Among the $\mathrm{deg}_{x}+1$ tensor legs of $|V_{x}\rangle$, the first one is a dangling edge, while the $\mathrm{deg}_{x}$ remaining ones will be contracted with the neighboring vertices to give the state $\Gamma_{G}$ (see Figure 11 for the case of a triangle graph):

[TABLE]

where $|\Omega\rangle\in\mathbb{C}^{d}\otimes\mathbb{C}^{d}$ is the maximally entangled state from (18).

For a fixed (possibly empty) subset $A\subseteq V$, consider the marginal

[TABLE]

It is important to note that neither $\Psi_{G}$ nor $\rho_{G;A}$ are normalized quantum states; this is due to the fact that the tensor contractions implemented by the maximally entangled states $|\Omega_{e}\rangle$ are not isometries. We define the normalized version (we drop the index $G$ to make the notation lighter)

[TABLE]

and we compute, as in [HNQ*+*16], its $2$-Renyi entropy

[TABLE]

where $F(x\otimes y)=y\otimes x$ is the flip operator, $\rho=\Psi\Psi^{*}=\rho_{\emptyset}$ is the (global) density operator, and $F_{AA^{\prime}}=F\otimes\mathrm{id}_{V\backslash A}\otimes\mathrm{id}_{V^{\prime}\backslash A^{\prime}}$ denotes the flip operator acting on the sites associated with $A$ (respectively its copy $A^{\prime}$). In [HNQ*+*16] it is argued that both the numerator and the denominator of the expression above concentrate around their average value (with respect to the randomness in the Haar unitaries $\{U_{x}\}_{x\in V}$), denoted in general by

[TABLE]

In [HNQ*+*16, Section 2.2], the value of averages of the form (20) was shown to be expressible – to lowest order in the limit $d\rightarrow\infty$ – in terms of the minimum energy associated with an Ising model with classical spins on the vertices of the graph $G$, and a magnetic field determined by the sites $A$ of interest. More precisely, the following expression was derived:

[TABLE]

where

[TABLE]

Obviously, in the case where $A=\emptyset$ (which corresponds to the normalizing denominator for the state $\rho_{A}$), the unique optimum is $s\equiv 1$, and we obtain

[TABLE]

Here we use the RTNI package to verify the validity of (21) and (22) for simple example of graphs $G$ 222In [HNQ*+*16], expressions (21) and (22) are subsequently used to derive Ryu-Takayanagi-type formulas for the Rényi entropy. We do not consider this here as this derivation appears to require additional assumptions on connectivity properties of the graph $G$, guaranteeing the minimum-energy considerations have a single domain wall (unlike in triangle graph example considered here). . In the case of the triangle graph $G_{\textrm{triangle}}$ shown in Figure 11, with a single point marginal $A=\{1,2\}$, we obtain again that the optimizing configuration is $s_{x}=1$ for all $x\in V$ and thus we obtain from (22) and (21)

[TABLE]

This result is confirmed by the following RTNI exact, non-asymptotic, computation (see also Figure 12). Note that the pre-factor $d^{-2|E|}$ appears below to take into account the normalization of the maximally entangled states (18), which is not encoded directly into the tensor network.

⬇

In[1]:= g = Graph[{1 <-> 2, 2 <-> 3, 3 <-> 1}];

In[2]:= edgeNormalization = d^(-2 EdgeCount[g]);

In[3]:= marginal = {1, 2};

In[4]:= tn = buildTNfromGraph[g, marginal];

In[5]:= Z0 = edgeNormalization integrateAllUs[tn[[1]], d, d][[1, 2]] //

FullSimplify

Out[1]= (3 + (-2 + d) d)/(d^7 (1 + d) (1 + (-1 + d) d)^3)

In[6]:= Limit[-Log[Z0]/Log[d], d -> Infinity]

Out[2]= 12

In[7]:= ZA = edgeNormalization integrateAllUs[tn[[2]], d, d][[1, 2]] //

FullSimplify

Out[3]= (1 + d^2)/(d^8 (1 + d) (1 + (-1 + d) d)^3)

In[8]:= Limit[-Log[ZA]/Log[d], d -> Infinity]

Out[4]= 13

Here the routine buildTNfromGraph creates two tensor networks starting from a graph $G$ and a subset $A\subset V$ of vertices. The tensor networks represent $\operatorname{Tr}[\rho\otimes\rho]$ and $\operatorname{Tr}[\rho_{A}^{2}]$, respectively (see Fig. 12). The corresponding code and additional examples for e.g., 2D grid graphs are provided in a separate Mathematica notebook.

6.6. Example computation of expectation values of matrix products

Here we illustrate the use of MultinomialexpectationvalueHaar. We use the routine to obtain

[TABLE]

for two $d\times d$-matrices $\alpha,\beta$ and a Haar-random unitary $U$ of size $d$.

⬇

In[1]:= MultinomialexpectationvalueHaar[d,{1,2},{X,Y},True]

Out[2]= (Tr[X] Tr[Y])/d

In[3]:= MultinomialexpectationvalueHaar[d,{1,3},{X,Y},False]

Out[4]= 0

In[5]:= MultinomialexpectationvalueHaar[d,{2,3},{X,Y},False]

Out[6]= X.Transpose[Y]/d

In[7]:= MultinomialexpectationvalueHaar[d,{1,4},{X,Y},True]

Out[8]= Tr[X.Transpose[Y]]/d

For two additional matrices $V,W$, we can compute for example

[TABLE]

using the command

⬇

In[8]:= MultinomialexpectationvalueHaar[d,{1,2,3,4},{V,W,X,Y},False]

Out[9]= -(V.Transpose[Y].X.Transpose[W]/(4-5 d^2+d^4))

    -(V.**Transpose**[Y].X Tr[W])/(4-5 d^2+d^4)

    -(V.X.**Transpose**[W] Tr[Y])/(4-5 d^2+d^4)

    -(V.X Tr[W] Tr[Y])/(4-5 d^2+d^4)

Acknowledgments. I.N.’s research has been supported by the ANR projects StoQ ANR-14-CE25-0003-01 and NEXT ANR-10-LABX-0037-NEXT. I.N. and M.F. acknowledge the hospitality of the TU München, where part of this work was conducted, and are both supported by the PHC Sakura program (project number: 38615VA). R.K. acknowledges support by DFG project no. K05430/1-1 and the Technical University of Munich – Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement no. 291763. M.F. was financially supported by JSPS KAKENHI Grant Number JP16K00005.

Appendix A User interface of Python package

In this section, we show how the Python package works, by using the example in Section 6.2: $[\operatorname{id}\otimes\operatorname{Tr}](UAU^{*})$. First, the following code creates the initial graph and visualizes it.

⬇

In[1]: from IHU_source import *

In[2]: e1 = [[”A”, 1, ”out”, 1], [”U”, 1, ”in”, 1]]

In[3]: e2 = [[”A”, 1, ”out”, 2], [”U”, 1, ”in”, 2]]

In[4]: e3 = [[”U*”, 1, ”out”, 1], [”A”, 1, ”in”, 1]]

In[5]: e4 = [[”U*”, 1, ”out”, 2], [”A”, 1, ”in”, 2]]

In[6]: e5 = [[”U”, 1, ”out”, 2], [”U*”, 1, ”in”, 2]]

In[7]: g = [e1, e2, e3, e4, e5]

In[8]: gw = [g,1]

In[9]: visualizeTN(gw)

Note that the input line 1 shows that source codes are imported from the file “IHU_source.py”, which imports source codes from “WFG_source”.

One difference (compared to the Mathematica code) is in the input line 8, where the weight is explicitly specified, here as $1$. If one wants to use symbolic number for the initial weight, it must be defined as described below using the symbols command. Figure 13 shows the visualization of the initial tensor network before taking average. In this figure, matrices are represented as yellow disks and red arrows carry information about how matrices are connected. For example, consider the arrow from “U*1” to “A1”. The information “[in1:out1][in2:out2]” shows that this is a double edge, where the first bracket come from the edge “e3” and the second the edge “e4”.

Next, the following code averages the input graph over the unitary group and visualize it.

⬇

In[10]: k,n = symbols(’k n’)

In[11]: rm = [”U”,[n,k],[n,k],n*k]

In[12]: Eg = integrateHaarUnitary(gw,rm)

In[13]: print(Eg)

In[14]: visualizeTN(Eg)

Out[1]: [[[[[’@U*’, 1, ’in’, 1], [’@U’, 1, ’out’, 1]],

[[’A’, 1, ’out’, 1], [’A’, 1, ’in’, 1]],

[[’A’, 1, ’out’, 2], [’A’, 1, ’in’, 2]]], 1/n]]

Note that in the input line 2, $k$ and $n$ are declared to be symbols explicitly. One can define other symbols in the same way. The algorithm calculates the average in input line 11 and outputs the result in the next line as shown in output line 1. The output is again a list, where the outer bracket are there in case the output is a sum of several tensor networks with corresponding weights. In this example, since we only have a pair of $U$ and $U^{*}$, the average can be written as a single pair of a tensor network and a weight; the weight is now $1/n$ and the remainder of the output specifies the tensor network. The mark “@” corresponds to “dummy” in the Mathematica package. A visualization of the averaged tensor network is given in Figure 14. Loops are represented by orange disks. In this example, one orange disk is connected to the yellow disk labeled as “A1” by the red arrow with the information “[out1:in1][out2:in2]”. This means $\operatorname{Tr}A$.

Appendix B Generating Weingarten functions

In this section, we provide a brief explanation of how the Python algorithm for generating Weingarten functions works, and show how to use it.

B.1. How the algorithm works

The formula for Weingarten functions (9) consists of two ingredients: characters of symmetric groups and Schur polynomials (i.e., characters of the unitary group) evaluated at the identity. Computation of the latter is straightforward based on the formula (10). The calculation of the former is rather complicated. We use the Murnaghan–Nakayama rule  [Mur71, Nak41a, Nak41b]. There are recursive and non-recursive versions, and we choose to use a recursive one because it uses only Young diagrams, while the non-recursive ones requires manipulation of Young Tableaux. We refer to [Ber04] for a more detailed analysis of the problem of computing character tables of symmetric groups.

To make our paper self-contained, let us explain this recursive formula. The character table of the symmetric group $\mathcal{S}_{p}$ is a square table, with entries $\chi(\alpha,\beta)$ indexed by two ordered partitions $\alpha,\beta$ of $p$. Here $\alpha$ denotes a Young diagram specifying an irreducible representation, and $\chi(\alpha,\beta)=\chi_{\alpha}(\pi)$ is the character evaluated for any permutation $\pi\in\mathcal{S}_{p}$ with cycle type $\beta$. A border strip is a connected set of boxes contained in a Young diagram. Now, fix $\alpha=(\alpha_{1},\ldots,\alpha_{m})$ and $\beta=(\beta_{1},\ldots,\beta_{n})$ and take the Young diagram of the form $\alpha$ and find all border strips within $\alpha$, denoted by $\gamma=(\gamma_{1},\ldots,\gamma_{m})$, such that

1. (1)

$\displaystyle\sum_{i=1}^{m}\gamma_{i}=\beta_{1}$ 2. (2)

removing $\gamma$ from the Young diagram yields another Young diagram 3. (3)

$\gamma$ does not contain $2\times 2$ square.

We call any such border strip valid. From the second condition above it follows that the boxes of a valid border strip occupy the right-hand-side of each row of the Young diagram.

Let us consider some examples. Suppose $\alpha=(5,2,2,1)$. One can find two valid border strips of size $3$ (Figure 15) but none of size $5$ (Figure 16).

Then, calculate the height $h(\gamma)$ of each border strip $\gamma$, which is defined as

[TABLE]

(In the notation introduced in Figure 15, the number of rows is the number of non-zero entries of $\gamma$.) The quantity $\chi(\alpha,\beta)$ then is given by the recursive formula

[TABLE]

where the sum is over all valid border strips $\gamma$ for the diagram $\alpha$. Here the expression $\alpha-\gamma$ refers to the Young diagram obtained by removing the boxes of $\gamma$ from $\alpha$.

Let us illustrate this using the example of the diagram $\alpha=(5,2,2,1)$ considered above. One can compute for example

[TABLE]

where the sume was taken over the valid border strips $(3,0,0,0)$ and $(0,0,2,1)$. By contrast,

[TABLE]

because there is no border strip of size $5$ which can be removed from the Young diagram of $(5,2,2,1)$.

Successive application of the recursive rule (32) leads to a linear combination of $\chi((1),(1))$ and $\chi((),())$, which are both defined to be $1$. However, our program uses the rule (32) just once, and computes the character table of $\mathcal{S}_{p}$ based on all the character tables of $\mathcal{S}_{1},\ldots,\mathcal{S}_{p-1}$. This is why if one computes Weingarten functions for $\mathcal{S}_{p}$, the program automatically generates the character tables of $\mathcal{S}_{1},\ldots,\mathcal{S}_{p}$.

B.2. How to use the Python RTNI package for generating Weingarten functions

One can generate Weingarten functions using the Python package. Example code is given here:

⬇

In[1]: from WFG_source import *

In[2]: k = 3

In[3]: display=’yes’

In[4]: record = ’yes’

In[5]: weigartenFunctionGenerator(k,display,record)

Out[1]: [[[1, 1, 1], (n2 - 2)/(n*(n4 - 5*n**2 + 4))],

[[2, 1], -1/(n4 - 5*n2 + 4)], [[3], 2/(n*(n4 - 5*n2 + 4))]]

where $k$ is the size of the symmetric group. The output gives the Weingarten function for dimension $n$ (i.e., for the unitary group on $\mathbb{C}^{n}$) as a list of pairs $(\alpha,\operatorname{Wg}_{n}(\pi_{\alpha}))$ where $\alpha$ (an ordered partition of $k$) is the cycle type and $\pi_{\alpha}\in\mathcal{S}_{k}$ is any permutation with cycle type $\alpha$.

By choosing ’yes’ or ’no’, one can specify if the algorithm displays the result and if it should be recorded as ”functions$k$” in the folder ”Weingarten”. In case the folder is missing, it is created automatically. In the example, $k=3$ and Weingarten functions for $\mathcal{S}_{3}$ and $\mathcal{U}(n)$ are calculated in the output, being represented as a list. For example, the Weingarten function for the partition (cycle type) $[1,1,1]$ is $(n^{2}-2)/(n(n^{4}-5n^{2}+4))$.

Moreover, one can generate only the character tables of symmetric groups with the following code.

⬇

In[1]: from WFG_source import *

In[2]: k = 3

In[3]: display=’yes’

In[4]: characterTableGenerator(k,display)

Out[1]: {’[1, 1, 1][1, 1, 1]’: 1, ’[1, 1, 1][2, 1]’: -1, ’[1, 1, 1][3]’: 1,

’[2, 1][1, 1, 1]’: 2, ’[2, 1][2, 1]’: 0, ’[2, 1][3]’: -1, ’[3][1, 1, 1]’: 1,

’[3][2, 1]’: 1, ’[3][3]’: 1}

Here, $k$ is the size of the symmetric group, as before, and it is possible to choose if the algorithm shows the result or not. The result will be recorded automatically in the folder “SGC”, which is created in case it is missing. In this example, $k=3$ and each of the output tuples stands for an element of character table. For example, by using the notation $\chi$ in Section B.1, $\chi((2,1),(3))=-1$. All the above data is stored as “.pkl” files, which can be translated to text files by the program called “pkl2text”.