Deep quantum neural networks form Gaussian processes
It is well known that artificial neural networks initialized from independent
and identically distributed priors converge to Gaussian processes in the limit
of large number of neurons per hidden layer. In this work we prove an analogous
result for Quantum Neural Networks (QNNs). Namely, we show that the outputs of
certain models based on Haar random unitary or orthogonal deep QNNs converge to
Gaussian processes in the limit of large Hilbert space dimension $d$. The
derivation of this result is more nuanced than in the classical case due to the
role played by the input states, the measurement observable, and the fact that
the entries of unitary matrices are not independent. An important consequence
of our analysis is that the ensuing Gaussian processes cannot be used to
efficiently predict the outputs of the QNN via Bayesian statistics.
Furthermore, our theorems imply that the concentration of measure phenomenon in
Haar random QNNs is worse than previously thought, as we prove that expectation
values and gradients concentrate as $\mathcal{O}\left(\frac{1}{e^d
\sqrt{d}}\right)$. Finally, we discuss how our results improve our
understanding of concentration in $t$-designs.
Diego García-Martín, Martin Larocca, M. Cerezo
# Deep quantum neural networks form Gaussian processes

###### Abstract

It is well known that artificial neural networks initialized from independent and identically distributed priors converge to Gaussian processes in the limit of large number of neurons per hidden layer. In this work we prove an analogous result for Quantum Neural Networks (QNNs). Namely, we show that the outputs of certain models based on Haar random unitary or orthogonal deep QNNs converge to Gaussian processes in the limit of large Hilbert space dimension \(d\). The derivation of this result is more nuanced than in the classical case due the role played by the input states, the measurement observable, and the fact that the entries of unitary matrices are not independent. An important consequence of our analysis is that the ensuing Gaussian processes cannot be used to efficiently predict the outputs of the QNN via Bayesian statistics. Furthermore, our theorems imply that the concentration of measure phenomenon in Haar random QNNs is much worse than previously thought, as we prove that expectation values and gradients concentrate as \(\mathcal{O}\left(\frac{1}{e^{d}\sqrt{d}}\right)\) - exponentially in the Hilbert space dimension. Finally, we discuss how our results improve our understanding of concentration in \(t\)-designs.

Neural Networks (NNs) have revolutionized the fields of Machine Learning (ML) and artificial intelligence. Their tremendous success across many fields of research in a wide variety of applications [1; 2; 3] is certainly astonishing. While much of this success has come through heuristics, the past few decades have witnessed a significant increase in our theoretical understanding of their inner workings. One of the most interesting results regarding NNs is that fully-connected models with a single hidden layer converge to Gaussian Processes (GPs) in the limit of large number of hidden neurons, when the parameters are initialized from independent and identically distributed (i.i.d.) priors [4]. More recently, it has been shown that i.i.d.-initialized, fully-connected, multi-layer NNs also converge to GPs in the infinite-width limit [5]. Furthermore, other architectures, such as convolutional NNs [6], transformers [7] or recurrent NNs [8] are also GPs under certain assumptions. More than just a mathematical curiosity, the correspondence between NNs and GPs opened up the possibility of performing exact Bayesian inference for regression and learning tasks using wide NNs [9; 4].

With the advent of quantum computers, there has been an enormous interest in merging quantum computing with ML, leading to the thriving field of Quantum Machine Learning (QML) [10; 11; 12; 13; 14]. Rapid progress has been made in this field, largely fueled by the hope that QML may provide a quantum advantage in the near-term for some practically-relevant problems. While the prospects for such a practical quantum advantage remain unclear [15], a number of promising analytical results have already been put forward [16; 17; 18; 19]. Still, much remains to be known about QML models.

In this work, we contribute to the QML body of knowledge by proving that under certain conditions, the outputs of deep Quantum Neural Networks (QNNs) - i.e., parametrized quantum circuits acting on input states drawn from a training set- converge to GPs in the limit of large Hilbert space dimension (see Fig. 1). Our results are derived for QNNs that are Haar random over the unitary and orthogonal groups. Unlike in the classical case, where the proof of the emergence of GPs stems from the central limit theorem, the situation becomes more intricate in the quantum setting as the entries of the QNN are not independent - the rows and columns of a unitary matrix are constrained to be mutually orthonormal. Hence, our proof strategy boils down to showing that each moment of the QNN's output distribution converges to that of a multivariate Gaussian. In addition, we show that in contrast to classical NNs, the Bayesian distribution of the QNN is inefficient for predicting the model's outputs. We then use our results to provide a precise characterization of the concentration of measure phenomenon in deep random quantum circuits [20; 21; 22; 23; 24; 25]. Here, our theorems indicate that the expectation values, as well as the gradients, of Haar random processes concentrate exponentially faster than reported in previous barren plateau studies [20; 21]. Finally, we discuss how our results can be leveraged to study QNNs that are not fully Haar random but instead form \(t\)-designs, which constitutes a much more practical assumption [26; 27; 28].
Gaussian processes and classical machine learning

We begin by introducing GPs.

**Definition 1** (Gaussian process).: _A collection of random variables \(\{X_{1},X_{2},\ldots\}\) is a GP if and only if, for every finite set of indices \(\{1,2,\ldots,m\}\), the vector \((X_{1},X_{2},\ldots,X_{m})\) follows a multivariate Gaussian distribution, which we denote as \(\mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})\). Said otherwise, every linear combination of \(\{X_{1},X_{2},\ldots,X_{m}\}\) follows a univariate Gaussian distribution._

In particular, \(\mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})\) is determined by its \(m\)-dimensional mean vector \(\mathbf{\mu}=(\mathbb{E}[X_{1}],\ldots,\mathbb{E}[X_{m}])\), and its \(m\times m\) dimensional covariance matrix with entries \((\mathbf{\Sigma})_{\alpha\beta}=\mathrm{Cov}[X_{\alpha},X_{\beta}]\).

GPs are extremely important in ML since they can be used as a form of kernel method to solve learning tasks [4, 9]. For instance, consider a regression problem where the data domain is \(\mathscr{X}=\mathbb{R}\) and the label domain is \(\mathscr{Y}=\mathbb{R}\). Instead of finding a single function \(f:\mathscr{X}\rightarrow\mathscr{Y}\) which solves the regression task, a GP instead assigns probabilities to a set of possible \(f(x)\), such that the probabilities are higher for the "more likely" functions. Following a Bayesian inference approach, one then selects the functions that best agree with some set of empirical observations [9, 14].

Under this framework, the output over the distribution of functions \(f(x)\), for \(x\in\mathscr{X}\), is a random variable. Then, given a set of training samples \(x_{1},\ldots,x_{m}\), and some covariance function \(\kappa(x,x^{\prime})\), Definition 1 implies that if one has a GP, the outputs \(f(x_{1}),\ldots,f(x_{m})\) are random variables sampled from some multivariate Gaussian distribution \(\mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})\). From here, the GP is used to make predictions about the output \(f(x_{m+1})\) (for some new data instance \(x_{m+1}\)), given the previous observations \(f(x_{1}),\ldots,f(x_{m})\). Explicitly, one constructs the joint distribution \(P(f(x_{1}),\ldots,f(x_{m}),f(x_{m+1}))\) from the averages and the covariance function \(\kappa\), and obtains the sought-after "predictive distribution" \(P(f(x_{m+1})|f(x_{1}),\ldots,f(x_{m}))\) via marginalization. The power of the GP relies on the fact that this distribution usually contains less uncertainty than \(P(f(x_{m+1}))=\mathcal{N}(\mathbb{E}[f(x_{m+1})],\kappa(x_{m+1},x_{m+1}))\) (see the Methods).

## II Haar random deep QNNs form GPs

In what follows we consider a setting where one is given repeated access to a dataset \(\mathscr{D}\) containing pure quantum states \(\{\rho_{i}\}_{i}\) on a \(d\)-dimensional Hilbert space. We will make no assumptions regarding the origin of these states, as they can correspond to classical data encoded in quantum states [29, 30], or quantum data obtained from some quantum mechanical process [31, 32]. Then, we assume that the states are sent through a deep QNN, denoted \(U\). While in general \(U\) can be parametrized by some set of trainable parameters \(\mathbf{\theta}\), we leave such dependence implicit for the ease of notation. At the output of the circuit one measures the expectation value of a traceless Hermitian operator taken from a set \(\mathscr{O}=\{O_{j}\}_{j}\) such that \(\mathrm{Tr}[O_{j}O_{j^{\prime}}]=d\delta_{j,j^{\prime}}\) and \(O_{j}^{2}=\mathds{1}\), for all \(j,j^{\prime}\) (e.g., Pauli strings). We denote the QNN outputs as

\[C_{j}(\rho_{i})=\mathrm{Tr}\big{[}U\rho_{i}U^{\dagger}O_{j}\big{]}\,. \tag{1}\]

Then, we collect these quantities over some set of states from \(\mathscr{D}\) and some set of measurements from \(\mathscr{O}\) in a vector

\[\mathscr{C}=(C_{j}(\rho_{i}),\ldots,C_{j^{\prime}}(\rho_{i^{\prime}}),\ldots). \tag{2}\]

As we will show below, in the large-\(d\) limit \(\mathscr{C}\) converges to a GP when the QNN unitaries \(U\) are sampled according to the Haar measure on the degree-\(d\) unitary \(\mathbb{U}(d)\) or orthogonal \(\mathbb{O}(d)\) groups (see Fig. 1). We will henceforth use the notation \(\mathbb{E}_{\mathbb{U}(d)}\) and \(\mathbb{E}_{\mathbb{O}(d)}\) to respectively denote Haar averages over \(\mathbb{U}(d)\) and \(\mathbb{O}(d)\). Moreover, we assume that when the circuit is sampled from \(\mathbb{O}(d)\), the states in \(\mathscr{D}\) and the measurement operators in \(\mathscr{O}\) are real valued.

Figure 1: **Schematic of our main results.** It is well known that certain classical NNs with \(N_{h}\) neurons per hidden layer become GPs when \(N_{h}\rightarrow\infty\). That is, given inputs \(x_{1}\) and \(x_{2}\), and corresponding outputs \(y_{1}\) and \(y_{2}\), then the joint probability \(P(y_{1},y_{2})\) is a multivariate Gaussian \(\mathcal{N}(\mathbf{0},\mathbf{\Sigma})\). In this work, we show that a similar result holds under certain conditions for deep QNNs in the limit of large Hilbert space dimension, \(d\rightarrow\infty\). Now, given quantum states \(\rho_{1}\) and \(\rho_{2}\), \(C(\rho)=\mathrm{Tr}[U\rho U^{\dagger}O]\) is such that \(P(C(\rho_{1}),C(\rho_{2}))=\mathcal{N}(\mathbf{0},\mathbf{\Sigma})\).
### Moment computation in the large-\(d\) limit

As we discuss in the Methods, we cannot rely on simple central-limit-theorem arguments to show that \(\mathscr{C}\) forms a GP. Hence, our proof strategy is based on computing all the moments of the vector \(\mathscr{C}\) and showing that they asymptotically match those of a multivariate Gaussian distribution. To conclude the proof we show that these moments unequivocally determine the distribution, for which we can use Carleman's condition [33, 34]. We refer the reader to the Supplemental Information (SI) for the detailed proofs of the results in this manuscript.

First, we present the following lemma.

**Lemma 1**.: _Let \(C_{j}(\rho_{i})\) be the expectation value of a Haar random QNN as in Eq. (1). Then for any \(\rho_{i}\in\mathscr{D}\), \(O_{j}\in\mathscr{O}\), \(\mathbb{E}_{\mathbb{U}(d)}[C_{j}(\rho_{i})]=\mathbb{E}_{\mathbb{O}(d)}[C_{j} (\rho_{i})]=0\,.\) (3) Moreover, for any pair of states \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\) and operators \(O_{j},O_{j^{\prime}}\in\mathscr{O}\) we have_

\[\operatorname{Cov}_{\mathbb{U}(d)}[C_{j}(\rho_{i})C_{j^{\prime}}(\rho_{i^{ \prime}})]=\operatorname{Cov}_{\mathbb{O}(d)}[C_{j}(\rho_{i})C_{j^{\prime}}( \rho_{i^{\prime}})]=0\,,\]

_if \(j\neq j^{\prime}\) and_

\[\mathbf{\Sigma}_{i,i^{\prime}}^{\mathbb{U}} =\frac{d}{d^{2}-1}\left(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime }}]-\frac{1}{d}\right)\,, \tag{4}\] \[\mathbf{\Sigma}_{i,i^{\prime}}^{\mathbb{O}} =\frac{2(d+1)}{(d+2)(d-1)}\left(\operatorname{Tr}[\rho_{i}\rho_{ i^{\prime}}]\left(1-\frac{1}{d+1}\right)-\frac{1}{d+1}\right)\,, \tag{5}\]

_if \(j=j^{\prime}\). Here, we have defined \(\mathbf{\Sigma}_{i,i^{\prime}}^{G}=\operatorname{Cov}_{G}[C_{j}(\rho_{i})C_{j}( \rho_{i^{\prime}})]\), where \(G=\mathbb{U}(d),\mathbb{O}(d)\)._

Lemma 1 shows that the expectation value of the QNN outputs is always zero. More notably, it indicates that the covariance between the outputs is null if we measure different observables (even if we use the same input state and the same circuit). This implies that the distributions \(C_{j}(\rho_{i})\) and \(C_{j^{\prime}}(\rho_{i^{\prime}})\) are independent if \(j\neq j^{\prime}\). That is, knowledge of the measurement outcomes for one observable and different input states does not provide any information about the outcomes of other measurements, at these or any other input states. Therefore, in what follows we will focus on the case where \(\mathscr{C}\) contains expectation values for different states, but the same operator. In this case, Lemma 1 shows that the covariances will be positive, zero, or negative depending on whether \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]\) is larger, equal, or smaller than \(\frac{1}{d}\), respectively.

We now state a useful result.

**Lemma 2**.: _Let \(\mathscr{C}\) be a vector of \(k\) expectation values of a Haar random QNN as in Eq. (2), where one measures the same operator \(O_{j}\) over a set of \(k\) states \(\rho_{1},\dots,\rho_{k}\in\mathscr{D}\). In the large-\(d\) limit, if \(k\) is odd then \(\mathbb{E}_{\mathbb{U}(d)}\left[C_{j}(\rho_{1})\cdots C_{j}(\rho_{k})\right]= \mathbb{E}_{\mathbb{O}(d)}\left[C_{j}(\rho_{1})\cdots C_{j}(\rho_{k})\right]=0\). Moreover, if \(k\) is even and if a) \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{ \operatorname{poly}(\log(d))}\right)\) for all \(i,i^{\prime}\), we have_

\[\mathbb{E}_{\mathbb{U}(d)}\left[C_{j}(\rho_{1})\cdots C_{j}(\rho _{k})\right] =\frac{1}{d^{k/2}}\sum_{\sigma\in T_{k}}\prod_{\{t,t^{\prime}\} \in\sigma}\operatorname{Tr}[\rho_{t}\rho_{t^{\prime}}] \tag{6}\] \[=\frac{\mathbb{E}_{\mathbb{O}(d)}\left[C_{j}(\rho_{1})\cdots C_{ j}(\rho_{k})\right]}{2^{k/2}}\,,\]

_where the summation runs over all the possible disjoint pairing of indexes in the set \(\{1,2,\dots,k\}\), \(T_{k}\), and the product is over the different pairs in each pairing; while if b) \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]=0\) for all \(i,i^{\prime}\), we have_

\[\mathbb{E}_{\mathbb{U}(d)}\left[C_{j}(\rho_{1})\cdots C_{j}(\rho _{k})\right] =\frac{k!}{2^{k/2}(k/2)!}\frac{1}{d^{k}} \tag{7}\] \[=\frac{\mathbb{E}_{\mathbb{O}(d)}\left[C_{j}(\rho_{1})\cdots C_{ j}(\rho_{k})\right]}{2^{k/2}}\,.\]

Using Lemma 2 as our main tool, we will be able to prove that deep QNNs form GPs for different types of datasets. In Table 1 we present a summary of our main results.

### Positively correlated GPs

We begin by studying the case when the states in the dataset satisfy \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{ \operatorname{poly}(\log(d))}\right)\) for all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\). According to Lemma 1, this implies that the variables are positively correlated. In the large \(d\) limit, we can derive the following theorem.

**Theorem 1**.: _Under the same conditions for which Lemma 2(a) holds, the vector \(\mathscr{C}\) forms a GP with mean vector \(\mathbf{\mu}=\mathbf{0}\) and covariance matrix given by \(\mathbf{\Sigma}_{i,i^{\prime}}^{\mathbb{U}}=\frac{\mathbf{\Sigma}_{i,i^{\prime}}^{ \mathbb{O}}}{2}=\frac{\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]}{d}\)._

\begin{table}
\begin{tabular}{|l|c|c|c|} \hline Dataset. For all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\): & GP & Correlation & Statement \\ \hline \hline \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{ \operatorname{poly}(\log(d))}\right)\) & Yes & Positive & Theorem 1 \\ \hline \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]=\frac{1}{2}\) & Yes & Null & Theorem 2 \\ \hline \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]=0\) & Yes & Negative & Theorem 3 \\ \hline \end{tabular}
\end{table}
Table 1: **Summary of main results.** In the first column we present conditions for the states in the dataset under which the deep QNN’s outputs form GPs. In the remaining columns we report the correlation in the GP variables and the associated theorem where the main result is stated. In all cases we assume that we measure the same operator \(O_{j}\) for all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\). In Theorem 5, we extend some of these results to the cases where the conditions are only met on average when sampling states over \(\mathscr{D}\).
Theorem 1 indicates that the covariances for the orthogonal group are twice as large as those arising from the unitary group. In Fig. 2, we present results obtained by numerically simulating a unitary Haar random QNN for a system of \(n=18\) qubits. The circuits were sampled using known results for the asymptotics of the entries of unitary matrices [33]. In the left panels of Fig. 2, we show the corresponding two-dimensional GP obtained for two initial states that satisfy \(\mathrm{Tr}[\rho_{i},\rho_{i^{\prime}}]\in\Omega(1)\). Here, we see that the variables are positively correlated in accordance with the prediction in Theorem 1.

The fact that the outputs of deep QNNs form GPs unravels a deep connection between QNNs and quantum kernel methods. While it has already been pointed out that QNN-based QML constitutes a form of kernel-based learning [35], our results materialize this connection for the case of Haar random circuits. Notably, we can recognize that the kernel arising in the GP covariance matrix is proportional to the Fidelity kernel, i.e., to the Hilbert-Schmidt inner product between the data states [36, 37, 35]. Moreover, since the predictive distribution of a GP can be expressed as a function of the covariance matrix (see the Methods), and thus of the kernel entries, our results further cement the fact that quantum models such as those in Eq. (1) are functions in the reproducing kernel Hilbert space [35].

### Uncorrelated GPs

We now consider the case when \(\mathrm{Tr}[\rho_{i}\rho_{i^{\prime}}]=\frac{1}{d}\) for all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\). We find the following result.

**Theorem 2**.: _Let \(\mathscr{C}\) be a vector of \(k\) expectation values of operators in \(\mathscr{C}\) over a set of \(k\) states \(\rho_{1},\dots,\rho_{k}\in\mathscr{D}\). If \(\mathrm{Tr}[\rho_{i}\rho_{i^{\prime}}]=\frac{1}{d}\) for all \(i,i^{\prime}\), then in the large \(d\)-limit \(\mathscr{C}\) forms a GP with mean vector \(\mathbf{\mu}=\mathbf{0}\) and diagonal covariance matrix_

\[\mathbf{\Sigma}_{i,i^{\prime}}^{\mathrm{U}}=\frac{\mathbf{\Sigma}_{i,i^{\prime}}^{ \mathrm{O}}}{2}=\begin{cases}\frac{1}{d}&\text{if }i=i^{\prime}\\ 0&\text{if }i\neq i^{\prime}\end{cases}\quad. \tag{8}\]

In the right panel of Fig. 2, we plot the GP corresponding to two initial states such that \(\mathrm{Tr}[\rho_{i},\rho_{i^{\prime}}]=\frac{1}{d}\). In this case, the variables appear uncorrelated as predicted by Theorem 2.

### Negatively correlated GPs

Here we study the case when \(\mathrm{Tr}[\rho_{i}\rho_{i^{\prime}}]=0\) for all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\). We prove the following theorem.

**Theorem 3**.: _Under the same conditions for which Lemma 2(b) holds, the vector \(\mathscr{C}\) forms a GP with mean vector \(\mathbf{\mu}=\mathbf{0}\) and covariance matrix_

\[\mathbf{\Sigma}_{i,i^{\prime}}^{\mathrm{U}(d)}=\begin{cases}\frac{1}{d+1}\text{ if }i=i^{\prime}\\ -\frac{1}{(d^{2}-1)}\text{ if }i\neq i^{\prime}\end{cases}\quad, \tag{9}\]

_and_

\[\mathbf{\Sigma}_{i,i^{\prime}}^{\mathrm{O}(d)}=\begin{cases}\frac{2}{d+1}\text{ if }i=i^{\prime}\\ -\frac{1}{(d+2)(d-1)}\text{ if }i\neq i^{\prime}\end{cases}\quad. \tag{10}\]

Note that since we are working in the large \(d\) limit, we could have expressed the entries of the covariance matrices of Theorem 3 as \(\mathbf{\Sigma}_{i,i}^{\mathrm{U}}=\frac{\mathbf{\Sigma}_{i,i}^{\mathrm{O}}}{2}=\frac {1}{d}\), and \(\mathbf{\Sigma}_{i,i^{\prime}}^{\mathrm{U}}=\frac{\mathbf{\Sigma}_{i,i}^{\mathrm{O}}} {2}=-\frac{1}{d^{2}}\) for \(i\neq i^{\prime}\). However, we find it convenient to present their full form as it will be important below.

Figure 2: **Two-dimensional GPs.** We plot the joint probability density function, as well as its scaled marginals, for the measurement outcomes at the output of a unitary Haar random QNN acting on \(n=18\) qubits. The measured observable is \(O_{j}=Z_{1}\), where \(Z_{1}\) denotes the Pauli \(z\) operator on the first qubit. Moreover, the input states are: \(\rho_{1}=|0\rangle\langle 0|^{\otimes n}\) and \(\rho_{2}=|\mathrm{GHZ}\rangle\langle\mathrm{GHZ}|\) with \(|\mathrm{GHZ}\rangle=\frac{1}{\sqrt{2}}(|0\rangle^{\otimes n}+|1\rangle^{ \otimes n})\), for the left panel; \(\rho_{1}\) and \(\rho_{3}=|\Psi\rangle\langle\Psi|\) with \(|\Psi\rangle=\frac{1}{\sqrt{d}}\,|0\rangle^{\otimes n}+\sqrt{1-\frac{1}{d}}\,| 1\rangle^{\otimes n}\) for the right panel. In both cases we took \(10^{4}\) samples.
### Deep QNN outcomes, and their linear combination

In this section, and the following ones, we will study the implications of Theorems 1, 2 and 3. Unless stated otherwise, the corollaries we present can be applied to all considered datasets (see Table 1).

First, we study the univariate probability distribution \(P(C_{j}(\rho_{i}))\).

**Corollary 1**.: _Let \(C_{j}(\rho_{i})\) be the expectation value of a Haar random QNN as in Eq. (1). Then, for any \(\rho_{i}\in\mathscr{D}\) and \(O_{j}\in\mathscr{O}\), we have_

\[P(C_{j}(\rho_{i}))=\mathcal{N}(0,\sigma^{2})\,, \tag{11}\]

_where \(\sigma^{2}=\frac{1}{d},\frac{2}{d}\) when \(U\) is Haar random over \(\mathbb{U}(d)\) and \(\mathbb{O}(d)\), respectively._

Corollary 1 shows that when a single state from \(\mathscr{D}\) is sent through the QNN, and a single operator from \(\mathscr{O}\) is measured, the outcomes follow a Gaussian distribution with a variance that vanishes inversely proportional with the Hilbert space dimension. This means that for large problems sizes, we can expect the results to be extremely concentrated around their mean (see below for more details). In Fig. 3 we compare the predictions from Corollary 1 to numerical simulations. We find that the simulations match very closely our theoretical results for both the unitary and the orthogonal groups. Moreover, we can observe that the standard deviation for orthogonal Haar random QNNs is larger than that for unitary ones. In Fig. 3 we also plot the quotient \(\frac{\mathbb{E}[C_{j}(\rho_{i})^{k}]}{\mathbb{E}[C_{j}(\rho_{i})^{2}]^{k/2}}\) obtained from our numerics, and we verify that it follows the value \(\frac{k!}{2^{k/2}(k/2)!}\) of a Gaussian distribution.

At this point, it is worth making an important remark. According to Definition 1, if \(\mathscr{C}\) forms a GP, then any linear combination of its entries will follow a univariate Gaussian distribution. In particular, if \(\{C_{j}(\rho_{1}),C_{j}(\rho_{2}),\ldots,C_{j}(\rho_{m})\}\subseteq\mathscr{C}\), then \(P(C_{j}(\widetilde{\rho}))\) with \(\widetilde{\rho}=\sum_{i=1}^{m}c_{i}\rho_{i}\) will be equal to \(\mathcal{N}(0,\widetilde{\sigma}^{2})\) for some \(\widetilde{\sigma}\). Note that the real-valued coefficients \(\{c_{i}\}_{i=1}^{m}\) need not be a probability distribution, meaning that \(\widetilde{\rho}\) is not necessarily a quantum state. The previous then raises an important question: What happens if \(\widetilde{\rho}\propto\openone\)? A direct calculation shows that \(C_{j}(\widetilde{\rho})=\sum_{i=1}^{d}C_{j}(c_{i}\rho_{i})\propto\text{Tr} \big{[}U\openone U^{\dagger}O_{j}\big{]}=\text{Tr}[O_{j}]=0\). How can we then unify these two perspectives? On the one hand \(C_{j}(\widetilde{\rho})\) should be normally distributed, but on the other hand we know that it is always constant. To solve this issue, we note that the only dataset we considered where the identity can be constructed is the one where \(\text{Tr}[\rho_{i}\rho_{i^{\prime}}]=0\) for all \(i,i^{\prime}\)1. In that case, we can leverage Theorem 3 along with the identity \(\widetilde{\sigma}^{2}=\text{Var}_{G}[\sum_{i=1}^{d}C_{j}(\rho_{i})]=\sum_{i, i^{\prime}}\text{Cov}_{G}[C_{j}(\rho_{i}),C_{j}(\rho_{i^{\prime}})]\) to explicitly prove that \(\text{Var}_{G}[\sum_{i=1}^{d}C_{j}(\rho_{i})]=0\) (for \(G=\mathbb{U}(d),\mathbb{O}(d)\)). Hence, we find a zero-variance Gaussian distribution, i.e., a delta distribution in the QNN's outcomes (as expected).

Footnote 1: This follows from the fact that if \(\mathscr{D}\) contains a complete basis then for any \(\widetilde{\rho}\in\mathscr{D}^{\perp}\), one has that if \(\text{Tr}[\widetilde{\rho}_{i}]=0\) for all \(\rho_{i}\in\mathscr{D}\), then \(\widetilde{\rho}=0\). Here, \(\mathscr{D}^{\perp}\) denotes the kernel of the projector onto the subspace spanned by the vectors in \(\mathscr{D}\).

### Predictive power of the deep QNN's GP

Let us now study the predictive distribution of the QNN's GP. We consider a scenario where we send \(k\) states \(\rho_{1},\ldots,\rho_{k}\in\mathscr{D}\) to the QNN and measure the same operator \(O_{j}\) at its output. Moreover, we assume that there exists some (statistical) noise in the measurement process, so that

Figure 3: **Probability density function for \(C_{j}(\rho_{i})\), for Haar random QNNs and different problem sizes.** We consider unitary and orthogonal QNNs with \(n\)-qubits, and we take \(\rho_{i}=|0\rangle\langle 0|^{\otimes n}\), and \(O_{j}=Z_{1}\). The colored histograms are built from \(10^{4}\) samples in each case, and the solid black lines represent the corresponding Gaussian distributions \(\mathcal{N}\left(0,\sigma^{2}\right)\), where \(\sigma^{2}\) is given in Corollary 1. The insets show the numerical versus predicted value of \(\mathbb{E}[C_{j}(\rho_{i})^{k}]/\mathbb{E}[C_{j}(\rho_{i})^{2}]^{k/2}\). For a Gaussian distribution with zero mean, such quotient is \(\frac{k!}{2^{k/2}(k/2)!}\) (solid black line).
we actually estimate the quantities \(y(\rho_{i})=C_{j}(\rho_{i})+\varepsilon_{i}\), where the noise terms \(\varepsilon_{i}\) are assumed to be independently drawn from the same distribution \(P(\varepsilon_{i})=\mathcal{N}(0,\sigma_{N}^{2})\). For simplicity, we assume that the noise is given by finite sampling and that \(\sigma_{N}^{2}=\frac{1}{N}\), where \(N\in\mathcal{O}(\mathrm{poly}(\log(d)))\) is the number of shots used to estimate each \(y(\rho_{i})\). We then prove the following result.

**Theorem 4**.: _Consider a GP obtained from a Haar random QNN. Given the set of observations \((y(\rho_{1}),\ldots,y(\rho_{m}))\) obtained from \(N\in\mathcal{O}(\mathrm{poly}(\log(d)))\) measurements, then the predictive distribution of the GP is trivial:_

\[P(C_{j}(\rho_{m+1})|C_{j}(\rho_{1}),\ldots,C_{j}(\rho_{m}))=P(C_{j}(\rho_{m+1} ))=\mathcal{N}(0,\sigma^{2})\,,\]

_where \(\sigma^{2}\) is given by Corollary 1._

Theorem 4 shows that by spending only a polylogarithmic-in-\(d\) number of measurements, one cannot use Bayesian statistical theory to learn any information about new outcomes given previous ones.

### Concentration of measure

In this section we show that Corollary 1 provides a more precise characterization of the concentration of measure and the barren plateau phenomena for Haar random circuits than that found in the literature [20; 21; 22; 23; 24; 25]. First, it implies that deep orthogonal QNNs will exhibit barren plateaus, a result not previously known. Second, we recall that in standard barren plateau analyses, one only looks at the first two moments of the distribution of cost values \(C_{j}(\rho_{i})\) (or, similarly, of gradient values \(\partial_{\theta}C_{j}(\rho_{i})\)). Then one uses Chebyshev's inequality, which states that for any \(c>0\), the probability \(P(|X|\geqslant c)\leqslant\frac{\mathrm{Var}[X]}{c^{2}}\), to prove that \(P(|C_{j}(\rho_{i})|\geqslant c)\) and \(P(|\partial_{\theta}C_{j}(\rho_{i})|\geqslant c)\) are in \(\mathcal{O}(\frac{1}{d})\)[25; 21]. However, having a full characterization of \(P(C_{j}(\rho_{i}))\) allows us to compute tail probabilities and obtain a much tighter bound. For instance, for \(U\) being Haar random over \(\mathbb{U}(d)\), we find

**Corollary 2**.: _Let \(C_{j}(\rho_{i})\) be the expectation value of a Haar random QNN as in Eq. (1). Assuming that there exists a parametrized gate in \(U\) of the form \(e^{-i\theta H}\) for some Pauli operator \(H\), then_

\[P(|C_{j}(\rho_{i})|\geqslant c),\,P(|\partial_{\theta}C_{j}(\rho_{i})| \geqslant c)\in\mathcal{O}\left(\frac{1}{ce^{dc^{2}}\sqrt{d}}\right)\,.\]

Corollary 2 indicates that the QNN outputs, and their gradients, actually concentrate with a probability which vanishes exponentially with \(d\). In an \(n\)-qubit system, where \(d=2^{n}\), then \(P(|C_{j}(\rho_{i})|\geqslant c)\) and \(P(|\partial_{\theta}C_{j}(\rho_{i})|\geqslant c)\) are doubly exponentially vanishing with \(n\). The tightness of our bound arises from the fact that Chebyshev's inequality is loose for highly narrow Gaussian distributions. Corollary 2 also implies that the narrow gorge region of the landscape [25], i.e., the fraction of non-concentrated \(C_{j}(\rho_{i})\) values, also decreases exponentially with \(d\).

In the Methods we furthermore show how our results can be used to study the concentration of functions of QNN outcomes, e.g., standard loss functions used in the literature, like the mean-squared error.

### Implications for \(t\)-designs

We now note that our results allow us to characterize the output distribution for QNN's that form \(t\)-designs, i.e., for QNNs whose unitary distributions have the same properties up to the first \(t\) moments as sampling random unitaries from \(\mathbb{U}(d)\) with respect to the Haar measure. With this in mind, one can readily see that the following corollary holds.

**Corollary 3**.: _Let \(U\) be drawn from a \(t\)-design. Then, under the same conditions for which Theorems 1, 2 and 3 hold, the vector \(\mathscr{C}\) matches the first \(t\) moments of a GP._

Corollary 3 extends our results beyond the strict condition of the QNN being Haar random to being a \(t\)-design, which is a more realistic assumption [26; 27; 28]. In particular, we can study the concentration phenomenon in \(t\)-designs: using an extension of Chebyshev's inequality to higher order moments leads to \(P(|C_{j}(\rho_{i})|\geqslant c)\), \(P(|\partial_{\theta}C_{j}(\rho_{i})|\geqslant c)\in\mathcal{O}\left(\frac{ \left(2\lfloor\frac{t}{2}\rfloor\right)!}{\left\lfloor\frac{t}{2}\right\rfloor( dc^{2})\left\lfloor\frac{t}{2}\right\rfloor(\lfloor\frac{t}{2}\rfloor)!}\right)\) (see the SI for a proof). Note that for \(t=2\) we recover the known concentration result for barren plateaus, but for \(t\geqslant 4\) we obtain new polynomial-in-\(d\)-tigter bounds.

### Generalized datasets

Up to this point we have derived our theorems by imposing strict conditions on the overlaps between every pair of states in the dataset. However, we can extend these results to the cases where the conditions are only met on average when sampling states over \(\mathscr{D}\).

**Theorem 5**.: _The results of Theorems 1 and 2 will hold, on average, if \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\operatorname{Tr}[\rho_ {i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{\mathrm{poly}(\log(d))}\right)\) and \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\operatorname{Tr}[\rho_ {i}\rho_{i^{\prime}}]=\frac{1}{d}\), respectively._

As discussed in the Methods section these extensions make our results more practical as \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\operatorname{Tr}[\rho_ {i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{\mathrm{poly}(\log(d))}\right)\).
\(\Omega\left(\frac{1}{\mathrm{poly}(\log(d))}\right)\) holds on standard multi-class classification settings [29; 32], while \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\mathrm{Tr}[\rho_{i}\rho_ {i^{\prime}}]=\frac{1}{d}\) holds when the dataset is composed of Haar random states.

## III Discussion and Outlook

In this manuscript we have shown that under certain conditions, the output distribution of deep Haar random QNNs converges to a Gaussian process in the limit of large Hilbert space dimension. While this result had been conjectured in [13], a formal proof was still lacking. We remark that although our result mirrors its classical counterpart -that certain classical NNs form GPs-, there exist nuances that differentiate our findings from the classical case. For instance, we need to make assumptions on the states processed by the QNN, as well as on the measurement operator. Moreover, some of these assumptions are unavoidable, as Haar random QNNs will not necessarily always converge to a GP. As an example, we have that if \(O_{i}\) is a projector onto a computational basis state, then one recovers a Porter-Thomas distribution [38]. Ultimately, these subtleties arise because the entries of unitary matrices are not independent. In contrast, classical NNs are not subject to this constraint.

It is worth noting that our theorems have further implications beyond those discussed here. We envision that our methods and results will be useful in more general settings where Haar random unitaries / \(t\)-designs are considered, such as quantum information scramblers and black holes [24; 39; 40], many-body physics [41], quantum decoulpers and quantum error correction [42]. Finally, we leave for future work several potential generalizations of our results. For instance, one could envision proving a general result that combines our Theorems 1, 2, and 3 into a single setting. Moreover, it could be interesting to study if GPs arise in other architectures such as quantum convolutional neural networks [43], or re-uploading circuits [30; 44], among others.

## IV Methods

### Infinitely-wide neural networks as Gaussian processes

Here we will briefly review the seminal work of Ref. [4], which proved that artificial NNs with a single infinitely-wide hidden layer form GPs. Our main motivation for reviewing this result is that, as we will see below, the simple technique used in its derivation cannot be directly applied to the quantum case.

For simplicity let us consider a network consisting of a single input neuron, \(N_{h}\) hidden neurons, and a single output neuron (see Fig. 1). The input of the network is \(x\in\mathbb{R}\), and the output is given by

\[f(x)=b+\sum_{l=1}^{N_{h}}v_{l}h_{l}(x)\,, \tag{12}\]

where \(h_{l}(x)=\phi(a_{l}+u_{l}x)\) models the action of each neuron in the hidden layer. Specifically, \(u_{l}\) is the weight between the input neuron and the \(l\)-th hidden neuron, \(a_{l}\) is the respective bias and \(\phi\) is some (non-linear) activation function such as the hyperbolic tangent or the sigmoid function. Similarly, \(v_{l}\) is the weight connecting the \(l\)-th hidden neuron to the output neuron, and \(b\) is the output bias. From Eq. (12) we can see that the output of the NN is a weighted sum of the hidden neurons' outputs plus some bias.

Next, let us assume that the \(v_{l}\) and \(b\) are taken i.i.d. from a Gaussian distribution with zero mean and standard deviations \(\sigma_{v}/\sqrt{N_{h}}\) and \(\sigma_{b}\), respectively. Likewise, one can assume that the hidden neuron weights and biases are taken i.i.d. from some Gaussian distributions. Then, in the limit of \(N_{h}\to\infty\), one can conclude via the central limit theorem that, since the NN output is a sum of infinitely many i.i.d. random variables, then it will converge to a Gaussian distribution with zero mean and variance \(\sigma_{b}^{2}+\sigma_{v}^{2}\mathbb{E}[h_{l}(x)^{2}]\). Similarly, it can be shown that in the case of multiple inputs \(x_{1},\ldots,x_{m}\) one gets a multivariate Gaussian distribution for \(f(x_{1}),\ldots,f(x_{m})\), i.e., a GP [4].

Naively, one could try to mimic the technique in Ref. [4] to prove our main results. In particular, we could start by noting that \(C_{j}(\rho_{i})\) can always be expressed as

\[C_{j}(\rho_{i})=\sum_{k,k^{\prime},r,r^{\prime}=1}^{d}u_{kk^{\prime}}\rho_{k^ {\prime}r}u_{r^{\prime}r}^{*}o_{r^{\prime}k}\,, \tag{13}\]

where \(u_{kk^{\prime}}\), \(u_{r^{\prime}r}^{*}\), \(\rho_{k^{\prime}r}\) and \(o_{r^{\prime}l}\) are the matrix entries of \(U\) and \(U^{\dagger}\), \(\rho\) and \(O\), respectively. Although Eq. (13) is a summation over a large number of random variables, we cannot apply the central limit theorem (or its variants) here, since the matrix entries \(u_{kk^{\prime}}\) and \(u_{rr^{\prime}}^{*}\) are not independent [33].

In fact, the correlation between the entries in the same row, or column, of a Haar random unitary are of order \(\frac{1}{d}\), while those in different rows, or columns, are of order \(\frac{1}{d^{2}}\). This small, albeit critical, difference makes it such that we cannot simply use the central limit theorem to prove that \(\mathscr{C}\) converges to a GP. Instead, we need to rely on the techniques described in the main text.
### Learning with the Gaussian process

In this section we will review the basic formalism for learning with Gaussian processes. Let \(\mathbf{C}\) be a Gaussian process. Then, by definition, given a collection of inputs \(\{x_{i}\}_{i=1}^{m}\), \(\mathbf{C}\) is determined by its \(m\)-dimensional mean vector \(\mathbf{\mu}\), and its \(m\times m\)-dimensional covariance matrix \(\mathbf{\Sigma}\). In what follows we will assume that the mean of \(\mathbf{C}\) is zero, and that the entries of its covariance matrix are expressed as \(\kappa(x_{i},x_{i^{\prime}})\). That is,

\[P\!\left(\!\!\left(\!\!\begin{array}{c}C(x_{1})\\ \vdots\\ C(x_{m})\end{array}\!\!\right)\!\!\right)=\mathcal{N}\!\left(\!\mathbf{\mu}=\left( \begin{array}{c}0\\ \vdots\\ 0\end{array}\!\!\right),\mathbf{\Sigma}=\left(\!\!\begin{array}{ccc}\kappa(x_{ 1},x_{1})&\cdots&\kappa(x_{1},x_{m})\\ \vdots&\vdots\\ \kappa(x_{m},x_{1})&\cdots&\kappa(x_{m},x_{m})\end{array}\!\!\right)\!\!\right)\!.\]

The previous allows us to know that, _a priori_, the distribution of values for any \(f(x_{i})\) will take the form

\[P(C(x_{i}))=\mathcal{N}(0,\sigma_{i}^{2})\,, \tag{14}\]

with \(\sigma_{i}^{2}=\kappa(x_{i},x_{i})\).

Now, let us consider the task of using \(m\) observations, which we will collect in a vector \(\mathbf{y}\), to predict the value at \(x_{m+1}\). First, if the observations are noiseless, then \(\mathbf{y}=(y(x_{1}),\cdots,y(x_{m}))\) is equal to \(\mathbf{C}=(C(x_{1}),\cdots,C(x_{m}))\). That is, \(\mathbf{C}=\mathbf{y}\). Here, we can use the fact that \(C\) forms a Gaussian process to find [9; 45]

\[P(C(x_{m+1})|\mathbf{C}) =P(C(x_{m+1})|C(x_{1}),C(x_{2}),\ldots,C(x_{m}))\] \[=\mathcal{N}\left(\mu(C(x_{m+1})),\sigma^{2}(C(x_{m+1}))\right)\,, \tag{15}\]

where \(\mu(C(x_{m+1}))\) and \(\sigma^{2}(C(x_{m+1}))\) respectively denote the mean and variance of the associated Gaussian probability distribution, and which are given by

\[\mu(C(x_{m+1})) =\mathbf{m}^{T}\cdot\mathbf{\Sigma}^{-1}\cdot\mathbf{C} \tag{16}\] \[\sigma^{2}(C(x_{m+1})) =\sigma^{2}_{m+1}-\mathbf{m}^{T}\cdot\mathbf{\Sigma}^{-1}\cdot\mathbf{m}\,. \tag{17}\]

The vector \(\mathbf{m}\) has entries \(\mathbf{m}_{i}=\kappa(x_{m+1},x_{i})\). We can compare Eqs. (14) and (15) to see that using Bayesian statistics to obtain the predictive distribution of \(P(C(x_{m+1})|\mathbf{C})\) shifts the mean from zero to \(\mathbf{m}^{T}\cdot\mathbf{\Sigma}^{-1}\cdot\mathbf{C}\) and the variance is decreased from \(\sigma^{2}_{m+1}\) by a quantity \(\mathbf{m}^{T}\cdot\mathbf{\Sigma}^{-1}\cdot\mathbf{m}\). The decrease in variance follows from the fact that we are incorporating knowledge about the observations, and thus decreasing the uncertainty.

In a realistic scenario, we can expect that noise will occur during our observation procedure. For simplicity we model this noise as Gaussian noise, so that \(y(x_{i})=C(x_{i})+\varepsilon_{i}\), where the noise terms \(\varepsilon_{i}\) are assumed to be independently drawn from the same distribution \(P(\varepsilon_{i})=\mathcal{N}(0,\sigma_{N}^{2})\). Now, since we have assumed that the noise is drawn independently, we know that the likelihood of obtaining a set of observations \(\mathbf{y}\) given the model values \(\mathbf{C}\) is given by \(P(\mathbf{y}|\mathbf{C})=\mathcal{N}(\mathbf{C},\sigma_{N}^{2}\mathds{1})\). In this case, we can find the probability distribution [9; 45]

\[P(C(x_{m+1})|\mathbf{C}) =\int d\mathbf{C}P(x_{m+1}|\mathbf{C})P(\mathbf{C}|\mathbf{y})\] \[=\int d\mathbf{C}P(C(x_{m+1})|\mathbf{C})P(\mathbf{y}|\mathbf{C})P(\mathbf{C})/P(\mathbf{y})\] \[=\mathcal{N}\left(\widetilde{\mu}(C(x_{m+1})),\widetilde{\sigma}^ {2}(C(x_{m+1}))\right)\,, \tag{18}\]

where now we have

\[\widetilde{\mu}(C(x_{m+1})) =\mathbf{m}^{T}\cdot(\mathbf{\Sigma}+\sigma_{N}^{2}\mathds{1})^{-1}\cdot \mathbf{C} \tag{19}\] \[\widetilde{\sigma}^{2}(C(x_{m+1})) =\sigma_{m+1}^{2}-\mathbf{m}^{T}\cdot(\mathbf{\Sigma}+\sigma_{N}^{2} \mathds{1})^{-1}\cdot\mathbf{m}\,. \tag{20}\]

In the first and the second equality we have used the explicit decomposition of the probability, along with Bayes and marginalization rules. We can see that the probability is still governed by a Gaussian distribution but where the inverse of \(\mathbf{\Sigma}\) has been replaced by the inverse of \(\mathbf{\Sigma}+\sigma_{N}^{2}\mathds{1}\).

### Concentration of functions of QNN outcomes

In the main text we have evaluated the distribution of QNN outcomes and their linear combinations. However, in many cases one is also interested in evaluating a function of the elements of \(\mathscr{C}\). For instance, in a standard QML setting the QNN outcomes are used to compute some loss function \(\mathcal{L}(\mathscr{C})\) which one wishes to optimize [10; 11; 12; 13; 14]. While we do not aim here at exploring all possible relevant functions \(\mathcal{L}\), we will present two simple examples that showcase how our results can be used to study the distribution of \(\mathcal{L}(\mathscr{C})\), as well as its concentration.

First, let us consider the case when \(\mathcal{L}(C_{j}(\rho_{i}))=C_{j}(\rho_{i})^{2}\). It is well known that given a random variable with a Gaussian distribution \(\mathcal{N}(0,\sigma^{2})\), then its square follows a Gamma distribution \(\Gamma(\frac{1}{2},2\sigma^{2})\). Hence, we know that \(P(\mathcal{L}(C_{j}(\rho_{i})))=\Gamma(\frac{1}{2},2\sigma^{2})\). Next, let us consider the case when \(\mathcal{L}(C_{j}(\rho_{i}))=(C_{j}(\rho_{i})-y_{i})^{2}\) for \(y_{i}\in[-1,1]\). This case is relevant for supervised learning as the mean-squared error loss function is composed of a linear combination of such terms. Here, \(y_{i}\) corresponds to the label associated to the state \(\rho_{i}\). We can exactly compute all the moments of \(\mathcal{L}(C_{j}(\rho_{i}))\) as

\[\mathbb{E}_{G}[\mathcal{L}(C_{j}(\rho_{i}))^{k}]=\sum_{r=0}^{2k} \binom{2k}{r}\mathbb{E}_{G}[C_{j}(\rho_{i})^{r}](-y_{i})^{2k-r}\,, \tag{21}\]

for \(G=\mathbb{U}(d),\mathbb{O}(d)\). We can then use Lemma 2 to obtain

\[\mathbb{E}_{\mathbb{U}(d)}[C_{j}(\rho_{i})^{r}]=\frac{r!}{d^{r/2}2^{r/2}(r/2)!} =\frac{\mathbb{E}_{\mathbb{O}(d)}[C_{j}(\rho_{i})^{r}]}{2^{r/2}}\,,\]
if \(r\) is even, and \(\mathbb{E}_{\mathbb{U}(d)}[C_{j}(\rho_{i})^{r}]=\mathbb{E}_{\mathbb{O}(d)}[C_{j}( \rho_{i})^{r}]=0\) if \(r\) is odd. We obtain

\[\mathbb{E}_{\mathbb{U}(d)}[\mathcal{L}(C_{j}(\rho_{i}))^{k}]=\frac{2^{k}}{(-d) ^{k}}M\left(-k,\frac{1}{2},-\frac{dy^{2}}{2}\right)\,, \tag{22}\]

with \(M\) the Kummer's confluent hypergeometric function.

Furthermore, we can also study the concentration of \(\mathcal{L}(C_{j}(\rho_{i}))\) and show that \(P\left(|\mathcal{L}(C_{j}(\rho_{i}))-\mathbb{E}_{\mathbb{U}(d)}(\mathcal{L}(C _{j}(\rho_{i})))|\geqslant c\right)\), where the average \(\mathbb{E}_{\mathbb{U}(d)}(\mathcal{L}(C_{j}(\rho_{i})))=y_{i}^{2}+\frac{1}{d}\), is in \(\mathcal{O}\left(\frac{1}{|\sqrt{c}+y_{i}|e^{d|\sqrt{c}+y_{i}|^{2}}\sqrt{d}}\right)\).

### Motivation for the generalized datasets

In Theorem 5 we generalized the results of Theorems 1 and 2 to hold on average when a) \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\operatorname{Tr}[\rho_ {i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{\operatorname{poly}(\log(d))}\right)\) and b) \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}\operatorname{Tr}[\rho _{i}\rho_{i^{\prime}}]=\frac{1}{d}\), respectively. Interestingly, these two cases have practical relevance. Let us start with Case a). Consider a multiclass classification problem, where each state \(\rho_{i}\) in \(\mathscr{D}\) belongs to one of \(Y\) classes, with \(Y\in\mathcal{O}(1)\), and where the dataset is composed of an (approximately) equal number of states from each class. That is, for each \(\rho_{i}\) we can assign a label \(y_{i}=1,\ldots,Y\). Then, we assume that the classes are well separated in the Hilbert feature space, a standard and sufficient assumption for the model to be able to solve the learning task [29; 32]. By well separated we mean that

\[\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}] \in\Omega\left(\frac{1}{\operatorname{poly}(\log(d))}\right)\,, \quad\text{if}\quad y_{i}=y_{i^{\prime}}\,, \tag{23}\] \[\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}] \in\mathcal{O}\left(\frac{1}{2^{n}}\right)\,,\quad\text{if}\quad y _{i}\neq y_{i^{\prime}}\,. \tag{24}\]

In this case, it can be verified that for any pair of states \(\rho_{i}\) and \(\rho_{i^{\prime}}\) sampled from \(\mathscr{D}\), one has \(\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}[\operatorname{Tr}[ \rho_{i}\rho_{i^{\prime}}]]\in\Omega\left(\frac{1}{\operatorname{poly}(\log(d) )}\right)\).

Next, let us evaluate Case b). Such situation arises precisely if the sates in \(\mathscr{D}\) are Haar random states. Indeed, we can readily show that

\[\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\mathscr{D}}[ \operatorname{Tr}[\rho_{i}\rho^{\prime}]] =\mathbb{E}_{\rho_{i},\rho_{i^{\prime}}\sim\operatorname{Haar}}[ \operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]]\] \[=\int_{\mathbb{U}(d)}d\mu(U)\operatorname{Tr}\!\left[U\rho_{0}U^ {\dagger}V\rho_{0}^{\prime}V^{\dagger}\right]\] \[=\int_{\mathbb{U}(d)}d\mu(U)\operatorname{Tr}\!\left[U\rho_{0}U^ {\dagger}\rho_{0}^{\prime}\right]\] \[=\frac{\operatorname{Tr}[\rho_{0}]\operatorname{Tr}[\rho_{0}^{ \prime}]}{d}\] \[=\frac{1}{d}\,. \tag{25}\]

Here, in the first equality we have used that sampling Haar random pure states \(\rho_{i}\) and \(\rho_{i^{\prime}}\) from the Haar measure is equivalent to taking two reference states \(\rho_{0}\) and \(\rho_{0}^{\prime}\) and evolving them with Haar random unitaries. In the second equality we have used the left-invariance of the Haar measure, and in the third equality we have explicitly performed the integration (see the SI).

### Sketch of the proof of our main results

Since our main results are mostly based on Lemmas 1 and 2, we will here outline the main steps used to prove these Lemmas. In particular, to prove them we need to calculate, in the large \(d\) limit, quantities of the form

\[\mathbb{E}_{G}\left[\operatorname{Tr}\left[U^{\otimes k}\Lambda(U^{\dagger})^{ \otimes k}O^{\otimes k}\right]\right]\,, \tag{26}\]

for arbitrary \(k\), and for \(G=\mathbb{U}(d),\mathbb{O}(d)\). Here, the operator \(\Lambda\) is defined as \(\Lambda=\rho_{1}\otimes\cdots\otimes\rho_{k}\), where the pure states \(\rho_{i}\) belong to \(\mathscr{D}\), and where \(O\) is an operator in \(\mathscr{O}\). The first moment (\(k=1\)), \(\mathbf{\mu}\), and the second moments (\(k=2\)), \(\mathbf{\Sigma}_{i,\nu}^{G}\) can be directly computed using standard formulas for integration over the unitary and orthogonal groups (see the SI). This readily recovers the results in Lemma 1. However, for larger \(k\) a direct computation quickly becomes intractable, and we need to resort to asymptotic Weingarten calculations. More concretely, let us exemplify our calculations for the unitary group and for the case when the states in the dataset are such that \(\operatorname{Tr}\!\left[\rho_{i}\rho_{i_{i^{\prime}}}\right]\in\Omega\left( \frac{1}{\operatorname{poly}(\log(d))}\right)\) for all \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\). As shown in the SI, we can prove the following lemma.

**Lemma 3**.: _Let \(X\) be an operator in \(\mathcal{B}(\mathcal{H}^{\otimes k})\), the set of bounded linear operators acting on the \(k\)-fold tensor product of a \(d\)-dimensional Hilbert space \(\mathcal{H}\). Let \(S_{k}\) be the symmetric group on \(k\) items, and let \(P_{d}\) be the subsystem permuting representation of \(S_{k}\) in \(\mathcal{H}^{\otimes k}\). Then, for large Hilbert space dimension (\(d\to\infty\)), the twirl of \(X\) over \(\mathbb{U}(d)\) is_

\[\mathbb{E}_{\mathbb{U}(d)}[U^{\otimes k}X(U^{\dagger})^{\otimes k}]= \frac{1}{d^{k}}\sum_{\sigma\in S_{k}}\operatorname{Tr}[XP_{d}( \sigma)]P_{d}(\sigma^{-1})\] \[+\frac{1}{d^{k}}\sum_{\sigma,\Pi\in S_{k}}c_{\sigma,\pi} \operatorname{Tr}[XP_{d}(\sigma)]P_{d}(\pi)\,,\]

_where the constants \(c_{\sigma,\Pi}\) are in \(\mathcal{O}(1/d)\)._

We recall that the subsystem permuting representation of a permutation \(\sigma\in S_{k}\) is

\[P_{d}(\sigma)=\sum_{i_{1},\ldots,i_{k}=0}^{d-1}|i_{\sigma^{-1}(1)},\ldots,i_{ \sigma^{-1}(k)}\rangle\langle i_{1},\ldots,i_{k}|\,. \tag{27}\]
Lemma 3 implies that (26) is equal to

\[\begin{split}\mathbb{E}_{\mathbb{U}(d)}&\left[\operatorname{ Tr}\left[U^{\otimes k}\Lambda(U^{\dagger})^{\otimes k}O^{\otimes k}\right] \right]\\ &=\frac{1}{d^{k}}\sum_{\sigma\in S_{k}}\operatorname{Tr}[\Lambda P _{d}(\sigma)]\operatorname{Tr}\bigl{[}P_{d}(\sigma^{-1})O^{\otimes k}\bigr{]} \\ &+\frac{1}{d^{k}}\sum_{\sigma,\Pi\in S_{k}}c_{\sigma,\pi} \operatorname{Tr}[\Lambda P_{d}(\sigma)]\operatorname{Tr}\bigl{[}P_{d}(\pi)O^ {\otimes k}\bigr{]}\,.\end{split} \tag{28}\]

We now note that, by definition, since \(O\) is traceless and such that \(O^{2}=\openone\), then \(\operatorname{Tr}\bigl{[}P_{d}(\sigma)O^{\otimes k}\bigr{]}=0\) for odd \(k\) (and for all \(\sigma\)). This result implies that all the odd moments are exactly zero, and also that the non-zero contributions in Eq. (28) for the even moments come from permutations consisting of cycles of even length. We remark that as a direct consequence, the first moment, \(\mathbb{E}_{\mathbb{U}(d)}\left[\operatorname{Tr}\bigl{[}U\rho_{i}U^{\dagger} O\bigr{]}\right]\), is zero for any \(\rho_{i}\in\mathscr{D}\), and thus we have \(\boldsymbol{\mu}=\boldsymbol{0}\). To compute higher moments, we show that \(\operatorname{Tr}\bigl{[}P_{d}(\sigma)O^{\otimes k}\bigr{]}=d^{r}\) if \(k\) is even and \(\sigma\) is a product of \(r\) disjoint cycles of even length. The maximum of \(\operatorname{Tr}\bigl{[}P_{d}(\sigma)O^{\otimes k}\bigr{]}\) is therefore achieved when \(r\) is maximal, i.e., when \(\sigma\) is a product of \(k/2\) disjoint transpositions (cycles of length two), leading to \(\operatorname{Tr}\bigl{[}P_{d}(\sigma)O^{\otimes k}\bigr{]}=d^{k/2}\). Then, we look at the factors \(\operatorname{Tr}[\Lambda P_{d}(\sigma)]\) and include them in the analysis. We have that for all \(\pi\) and \(\sigma\) in \(S_{k}\),

\[\frac{1}{d^{k}}\Bigl{|}(c_{\sigma,\pi}\operatorname{Tr}[\Lambda P _{d}(\sigma)]\operatorname{Tr}\bigl{[}P_{d}(\pi)O^{\otimes k}\bigr{]} \tag{29}\] \[+c_{\sigma^{-1},\pi}\operatorname{Tr}\bigl{[}\Lambda P_{d}(\sigma ^{-1})\bigr{]}\bigr{)}\operatorname{Tr}\bigl{[}P_{d}(\pi)O^{\otimes k}\bigr{]} \Bigr{|}\in\mathcal{O}\left(\frac{1}{d^{\frac{k+2}{2}}}\right)\,.\]

Moreover, since \(\operatorname{Tr}[\rho_{i}\rho_{i^{\prime}}]\in\Omega\left(\frac{1}{\operatorname {poly}(\log(\log))}\right)\) for all pair of states \(\rho_{i},\rho_{i^{\prime}}\in\mathscr{D}\), it holds that

\[\frac{1}{d^{k}}\operatorname{Tr}[\Lambda P_{d}(\sigma)]\operatorname{Tr} \bigl{[}P_{d}(\sigma^{-1})O^{\otimes k}\bigr{]}\in\Omega\left(\frac{1}{d^{k/2 }}\right)\,, \tag{30}\]

if \(\sigma\) is a product of \(k/2\) disjoint transpositions, and

\[\begin{split}\frac{1}{d^{k}}\Bigl{|}& \operatorname{Tr}[\Lambda P_{d}(\sigma)]\operatorname{Tr}\bigl{[}P_{d}(\sigma^ {-1})O^{\otimes k}\bigr{]}+\\ &\operatorname{Tr}\bigl{[}\Lambda P_{d}(\sigma^{-1})\bigr{]} \operatorname{Tr}\bigl{[}P_{d}(\sigma)O^{\otimes k}\bigr{]}\Bigr{|}\in \mathcal{O}\left(\frac{1}{d^{\frac{k+2}{2}}}\right)\,,\end{split} \tag{31}\]

for any other \(\sigma\). We remark that if \(\sigma\) consist only of transpositions, then it is its own inverse, that is, \(\sigma=\sigma^{-1}\).

It immediately follows that for fixed \(k\) and \(d\to\infty\), the second sum in Eq. (28) is suppressed at least inversely proportional to the dimension of the Hilbert space with respect to the first one (i.e. exponentially in the number of qubits for QNNs made out of qubits). Likewise, the contributions in the first sum in (28) coming from permutations that are not the product of \(k/2\) disjoint transpositions are also suppressed inversely proportional to the Hilbert space dimension. Therefore, in the large \(d\) limit we arrive at

\[\mathbb{E}_{\mathbb{U}(d)}\left[\operatorname{Tr}\left[U^{\otimes k}\Lambda(U ^{\dagger})^{\otimes k}O^{\otimes k}\right]\right]=\frac{1}{d^{k/2}}\sum_{ \sigma\in T_{k}}\prod_{\{t,t^{\prime}\}\in\sigma}\operatorname{Tr}[\rho_{t} \rho_{t^{\prime}}]\,, \tag{32}\]

where we have defined as \(T_{k}\subseteq S_{k}\) the set of permutations which are exactly given by a product of \(k/2\) disjoint transpositions. Note that this is precisely the statement in Lemma 2.

From here we can easily see that if every state in \(\Lambda\) is the same, i.e., if \(\rho_{i}=\rho\) for \(i=1,\ldots,k\), then \(\operatorname{Tr}[\rho_{t}\rho_{t^{\prime}}]=1\) for all \(t,t^{\prime}\), and we need to count how many terms are there in Eq. (32). Specifically, we need to count how many different ways there exist to split \(k\) elements into pairs (with \(k\) even). A straightforward calculation shows that

\[\sum_{\sigma\in T_{k}}\prod_{\{t,t^{\prime}\}\in\sigma}1=\frac{1}{(k/2)!}{k\choose 2,2,\ldots,2}=\frac{k!}{2^{k/2}(k/2)!}\,. \tag{33}\]

Thus, we arrive at

\[\mathbb{E}_{\mathbb{U}(d)}\left[\operatorname{Tr}\left[U^{\otimes k}\Lambda(U^{ \dagger})^{\otimes k}O^{\otimes k}\right]\right]=\frac{1}{d^{k/2}}\frac{k!}{2^ {k/2}(k/2)!}\,. \tag{34}\]

Identifying \(\sigma^{2}=\frac{1}{d}\) implies that the moments \(\mathbb{E}_{\mathbb{U}(d)}\left[\operatorname{Tr}\left[U\rho U^{\dagger}O \right]^{k}\right]\) exactly match those of a Gaussian distribution \(\mathcal{N}(0,\sigma^{2})\).

To prove that these moments unequivocally determine the distribution of \(\mathscr{C}\), we use Carleman's condition.

**Lemma 4** (Carleman's condition, Hamburger case [34]).: _Let \(\gamma_{k}\) be the (finite) moments of the distribution of a random variable \(X\) that can take values on the real line \(\mathbb{R}\). These moments determine uniquely the distribution of \(X\) if_

\[\sum_{k=1}^{\infty}\gamma_{2k}^{-1/2k}=\infty. \tag{35}\]

Explicitly, we have

\[\sum_{k=1}^{\infty}\left(\frac{1}{d^{k}}\frac{(2k)!}{2^{k}k!} \right)^{-1/2k} =\sqrt{2d}\sum_{k=1}^{\infty}\left((2k)\cdots(k+1)\right)^{-1/2k}\] \[\geqslant\sum_{k=1}^{\infty}\left((2k)^{k}\right)^{-1/2k}\] \[=\sum_{k=1}^{\infty}\frac{1}{\sqrt{2k}}=\infty\,. \tag{36}\]
Hence, according to Lemma 4, Carleman's condition is satisfied, and \(P(C_{j}(\rho_{i})\) is distributed following a Gaussian distribution.

A similar argument can be given to show that the moments of \(\mathscr{C}\) match those of a GP. Here, we need to compare Eq. (32) with the \(k\)-th order moments of a GP, which are provided by Isserlis theorem [46]. Specifically, if we want to compute a \(k\)-th order moment of a GP, then we have that \(\mathbb{E}[X_{1}X_{2}\cdots X_{k}]=0\) if \(k\) is odd, and

\[\mathbb{E}[X_{1}X_{2}\cdots X_{k}]=\sum_{\sigma\in T_{k}}\prod_{\{t,t^{ \prime}\}\in\sigma}\mathrm{Cov}[X_{t},X_{t^{\prime}}]\,, \tag{37}\]

if \(k\) is even. Clearly, Eq. (32) matches Eq. (37) by identifying \(\mathrm{Cov}[X_{t},X_{t^{\prime}}]=\frac{\mathrm{Tr}[\rho_{t}\rho_{t^{\prime}} ]}{d}\). We can again prove that these moments uniquely determine the distribution of \(\mathscr{C}\) from the fact that since its marginal distributions are determinate via Carleman's condition (see above), then so is the distribution of \(\mathscr{C}\)[34]. Hence, \(\mathscr{C}\) forms a GP.

## Acknowledgements

We acknowledge Francesco Caravelli, Frederic Sauvage, Lorenzo Leone, and Cinthia Huerta for useful conversations. D.G-M. was supported by the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory (LANL) under project number 20230049DR. M.L. acknowledges support by the Center for Nonlinear Studies at Los Alamos National Laboratory (LANL). M.C. acknowledges support by the LDRD program of LANL under project number 20230527ECR. This work was also supported by LANL ASC Beyond Moore's Law project.
