Deep Quantum Neural Networks are Gaussian Process
The overparameterization of variational quantum circuits, as a model of
Quantum Neural Networks (QNN), not only improves their trainability but also
serves as a method for evaluating the property of a given ansatz by
investigating their kernel behavior in this regime. In this study, we shift our
perspective from the traditional viewpoint of training in parameter space into
function space by employing the Bayesian inference in the Reproducing Kernel
Hilbert Space (RKHS). We observe the influence of initializing parameters using
random Haar distribution results in the QNN behaving similarly to a Gaussian
Process (QNN-GP) at wide width or, empirically, at a deep depth. This outcome
aligns with the behaviors observed in classical neural networks under similar
circumstances with Gaussian initialization. Moreover, we present a framework to
examine the impact of finite width in the closed-form relationship using a $
1/d$ expansion, where $d$ represents the dimension of the circuit's Hilbert
space. The deviation from Gaussian output can be monitored by introducing new
quantum meta-kernels. Furthermore, we elucidate the relationship between GP and
its parameter space equivalent, characterized by the Quantum Neural Tangent
Kernels (QNTK). This study offers a systematic way to study QNN behavior in
over- and under-parameterized scenarios, based on the perturbation method, and
addresses the limitations of tracking the gradient descent methods for
higher-order corrections like dQNTK and ddQNTK. Additionally, this
probabilistic viewpoint lends itself naturally to accommodating noise within
our model.
Ali Rad
# Deep Quantum Neural Networks are Gaussian Process

###### Abstract

The overparameterization of variational quantum circuits, as a model of Quantum Neural Networks (QNN), not only improves their trainability but also serves as a method for evaluating the property of a given ansatz by investigating their kernel behavior in this regime. In this study, we shift our perspective from the traditional viewpoint of training in parameter space into function space by employing the Bayesian inference in the Reproducing Kernel Hilbert Space (RKHS). We observe the influence of initializing parameters using random Haar distribution results in the QNN behaving similarly to a Gaussian Process (QNN-GP) at wide width or, empirically, at a deep depth. This outcome aligns with the behaviors observed in classical neural networks under similar circumstances with Gaussian initialization. Moreover, we present a framework to examine the impact of finite width in the closed-form relationship using a \(1/d\) expansion, where \(d\) represents the dimension of the circuit's Hilbert space. The deviation from Gaussian output can be monitored by introducing new quantum meta-kernels. Furthermore, we elucidate the relationship between GP and its parameter space equivalent, characterized by the Quantum Neural Tangent Kernels (QNTK). This study offers a systematic way to study QNN behavior in over- and under-parameterized scenarios, based on the perturbation method, and addresses the limitations of tracking the gradient descent methods for higher-order corrections like dQNTK and ddQNTK. Additionally, this probabilistic viewpoint lends itself naturally to accommodating noise within our model.

## I Introduction

The capacity for learning and trainability in quantum circuits is an intriguing subject that merits further exploration and discovery. The issue of how trainable variational quantum circuits and quantum neural networks (QNNs) poses a significant challenge that needs to be tackled. The standard approach to training QNN is within the context of _parameter space_. The aim is to find the optimum point in the loss landscape using local gradient-based methods. However, these methods might encounter the Barren Plateau problem [14], where the magnitude of the gradient decreases exponentially as the size of the Hilbert space increases.

Alternatively, we can shift the training process into _functions space_, viewing QNNs as linear combinations of Kernels that exist within the Reproducing Kernel Hilbert Space (RKHS)[20; 21]. The learning task in the parametric space of gate's parameters can be transformed into the task of identifying the correct coefficients of kernel regression parameters, which are of the order of the input dataset. This is typically much less than the dimension of the Hilbert space. The process of identifying these coefficients is usually more convex than implementing the original gradient descent method in the parameter space.

Although the quantum kernel's behavior in different quantum circuit regimes is still a largely uncharted area, recent findings suggest that kernel-based methods can be more successful than gradient methods in mitigating the Barren Plateau problem in under-parameterized regimes[2; 18]. With this evidence, we can now turn our attention in this work to the other part of the spectrum: the over-parameterized regime. We hope that this will enhance our understanding and lay the groundwork for more generic and practical applications in the NISQ era. In this regime, both the functional and parametric methods present more straightforward and comprehensible approaches. Similar to classical Neural Networks (CNNs)[7], the output of the QNN can be determined using the kernel regression method. The kernel employed in this lazy training method is referred to as the Quantum Neural Tangent Kernel (QNTK)[1; 10; 12; 13; 22].

An interesting observation has been made regarding CNN when they are over-parameterized through wide layers or deep networks: their output tends to converge towards a Gaussian distribution and Gaussian Process (GP) [8; 9; 15; 16]. Given this observation, it's feasible to replace the prior distribution of parameters in classical neural networks with a Gaussian process prior. This process shifts the focus from parameter spaces to the space of functions. A Gaussian process is a probabilistic model that defines a distribution over functions. the idea is to use a prior probability on functions, which is specified as a Gaussian process that determines the shape of the function space, and then update this with the observed data to get a posterior distribution over functions. In contrast to kernel regression, in which one typically needs to make decisions about the kernel hyperparameters, while in Gaussian process regression, these are typically learned from the data. Furthermore, while kernel regression provides a point estimate, Gaussian process regression provides a full probabilistic model, allowing for uncertainty estimation in predictions. we can think of Gaussian process regression as a probabilistic, kernel-based method, and it's a generalization of kernel regression that also provides uncertainty estimates.
### Quantum Meta Kernels

We begin by study the standard structure of variational quantum circuits, integral to quantum neural network models. The structure is given by:

\[U(\mathbf{\theta})=\prod_{\ell=1}^{L}U_{i}(\theta_{i})W_{i}, \tag{1}\]

where, \(U_{i}(\theta_{i})\) represents a unitary gate with a variational parameter \(\theta_{i}\). \(W\) stands for a fixed and unparameterized section of the circuit. In a common scenario, \(U_{i}(\theta_{i})\) can be represented as \(U_{i}(\theta)=e^{i\theta_{i}X_{i}}\), where \(X_{i}\) is a Hermitian operator. The layer depth is denoted by \(L\), which is equivalent to the number of parameters \(\mathbf{\theta}=\{\theta_{1},\cdots,\theta_{L}\}\). By defining \(\tilde{U}_{i}(\theta_{i})(\rho):=U_{i}(\theta_{i})\rho U_{i}^{\dagger}(\theta _{i})\) ans \(\tilde{W}_{i}:=W_{i}\rho W_{i}^{\dagger}\). We can express the whole quantum circuit as a quantum model that for a given input \(\rho_{\alpha}\), it return an output (array or scalar):

\[f_{i,\alpha}=\text{Tr}(\mathcal{U}(\mathbf{\theta})(\rho_{\alpha})O_{i}) \tag{2}\]

such that

\[\mathcal{U}(\mathbf{\theta})(\cdot)=\bigcirc_{\ell=1}^{L}(\tilde{U}_{\ell}(\theta _{\ell})\circ\tilde{W}_{\ell})(\cdot) \tag{3}\]

In order to analyze the probabilistic behavior of the outputs from our quantum models, it's necessary to examine the k-point correlation function, given by

\[\mathbb{E}_{\mathbf{\theta}}[f_{i_{1};\alpha_{1}}(\mathbf{\theta})\cdots f_{i_{k}; \alpha_{k}}(\mathbf{\theta})]=\mathbb{E}_{\mathbf{\theta}}[\text{Tr}\Big{(}\mathcal{ U}^{(k)}(\mathbf{\theta})(\rho_{\alpha}^{\otimes k})O^{\otimes k}\Big{)}] \tag{4}\]

Where \(\mathcal{U}^{(k)}\) represents an extension of Eq.3, which is defined as follows:

\[\mathcal{U}^{(k)}(\mathbf{\theta})(\rho_{\alpha}^{\otimes k}):=\bigcirc_{\ell=1}^ {L}(\tilde{U}_{\ell}^{\otimes k}(\theta_{\ell})\circ\tilde{W}_{\ell}^{\otimes k })(\rho_{\alpha}^{\otimes k}) \tag{5}\]

If the gates in a quantum model satisfy the t-design condition

\[\frac{1}{|S|}\sum_{s=1}^{s=|S|}P_{t,t}(U)=\int_{\tilde{U}}d\mu(U)UP_{t,t}(U) \tag{6}\]

then, it is possible to analytically evaluate the expectation values. This holds true for overparameterized and deep circuits, where evidence suggests that the \(t\)-design condition can be met. In such cases, we can derive analytical and closed-form expressions for these expectations for a random Haar measure ensembles:

\[\int_{\tilde{U}(d)}U_{i_{1},j_{1}}\cdots U_{i_{p},j_{p}}U_{i^{ \prime},j^{\prime}}^{\dagger}\cdots U_{i^{\prime}_{p^{\prime}},j^{\prime}_{p}} ^{\dagger}\mu(U)\] \[= \sum_{\alpha,\beta\in S_{p}}\delta_{i_{1}i^{\prime}_{\alpha(1)}} \cdots\delta_{i_{p}i^{\prime}_{\alpha(p)}}\delta_{j_{1}j^{\prime}_{\beta(1)}} \cdots\delta_{j_{p}j^{\prime}_{\beta(p)}}W_{g,d}(\alpha^{-1}\beta), \tag{7}\]

where The \(W_{g,d}\) is Weingarten function defined on the symmetric group \(S_{p}\)[3; 17]

To gain a deeper understanding of the non-Gaussian characteristics of the distribution, it is often beneficial to study the connected \(n\)-point correlator. This can be achieved by removing the contributions from the _vacuum

Figure 1: Encoding the classical or quantum data can be seen in two way. (top): Encoding in Feature space which is Hilbert space, \(|\phi(x)\rangle\in F\). The quantum circuit evolves the feature vector using a parameter \(\theta\) to predict the label. Training aims to find the optimal \(\theta\) that achieves the highest criteria, such as minimizing the loss function, for accurate label prediction. (button):Encoding data as a quantum kernel \(K(\mathcal{X},\cdot)\). One example of quantum kernel is Fidelity kernel defines as: \(K^{F}(x,x^{\prime})=|\langle\phi(x),\phi(x^{\prime})\rangle|^{2}=\text{Tr}[ \rho(x)\rho(x^{\prime})]\). The Reproducing Kernel Hilbert Space (RKHS) is the feature space that has canonical feature space if for \(K(\cdot,x)\in\mathcal{H}\) for all \(x\in\mathcal{X}\) it has a reproducing property \(f(x)=\langle f,K(\cdot,x)\rangle_{\mathcal{H}}\). In this method, a training method for quantum circuits that seeks to find the optimal coefficients \(\{c_{i}\}\), which effectively describe the output function \(f(x)=\sum_{i}c_{i}K(x,x_{i})\) based on the observed data points \(\{x_{i}\}\). This technique is commonly referred to as kernel regression. This method can be extended to Bayesian inference which probabilistic kernel regression and predict the \(p(f)\) instead of \(f\).
state in our observable, according to the Wick's theorem:

\[\mathbb{E}_{\mathbf{\varrho}}[f_{1}f_{2}\cdots f_{n}]=\mathbb{E}_{\mathbf{ \varrho}}[f_{1}f_{2}\cdots f_{n}]|_{\text{connected}}\\ +\sum_{\text{all pairing}}\mathbb{E}_{\mathbf{\varrho}}[f_{a_{1,1}} \cdots f_{\alpha_{m_{1},1}}]\cdots\mathbb{E}_{\mathbf{\varrho}}[f_{\alpha_{m,1}} \cdots f_{\alpha_{m,m}}], \tag{8}\]

where the term "all-pairing" denotes all feasible subgroups with a collective count of \(m\). By adjusting the \(f\to f-\mathbb{E}[f]\) we can observe that \(\mathbb{E}[f]|_{\text{connected}}=\mathbb{E}[-f]_{\text{connected}}\) This assumption implies that the odd order of correlators will vanish.

Through careful calculation, as described in the appendix, we derive the following result for two-point correlator :

\[\mathbb{E}[f_{i_{1},\alpha_{1}}f_{i_{2},\alpha_{2}}]|_{\text{ connected}}\\ =\frac{\operatorname{Tr}\left(O_{i_{2}}O_{i_{1}}\right)}{d-d^{3} }+\frac{\operatorname{Tr}\left(O_{i_{2}}O_{i_{1}}\right)\operatorname{Tr} \left(\rho_{\alpha_{1}}\rho_{\alpha_{2}}\right)}{d^{2}-1}\\ +\frac{\operatorname{Tr}\left(O_{i_{1}}\right)\operatorname{Tr} \left(O_{i_{2}}\right)}{d^{2}-1}+\frac{\operatorname{Tr}\left(O_{i_{1}}\right) \operatorname{Tr}\left(O_{i_{2}}\right)\operatorname{Tr}\left(\rho_{\alpha_{1} }\right)\operatorname{Tr}\left(\rho_{\alpha_{2}}\right)}{d^{2}} \tag{9}\]

and similarly, for four-point correlator we obtain:

\[\mathbb{E}[f_{i_{1},\alpha_{1}}f_{i_{2},\alpha_{2}}f_{i_{3}, \alpha_{3}}f_{i_{4},\alpha_{4}}]|_{\text{connected}}\\ =\frac{d^{4}-8d^{2}+6}{d^{2}(d^{6}-14d^{4}+49d^{2}-36)}\\ [\operatorname{Tr}\left(\rho_{\alpha_{1}}\rho_{\alpha_{2}}\rho_{ \alpha_{3}}\rho_{\alpha_{4}}\right)\times\mathcal{V}_{4}(\{O_{i}\})\\ +\operatorname{Tr}\left(\{\rho_{\alpha_{i}}\rho_{\alpha_{j}}\rho _{\alpha_{k}}\}\right)\times\mathcal{V}_{3}(\{O_{i}\}))\\ +\mathcal{O}(\frac{1}{d^{5}}). \tag{10}\]

Here, \(\mathcal{V}_{k}\) represents a linear combination of functions that map subsets of \(\{O_{i}\},\{\rho_{i}\},i\in[1,\cdots k]\) to scalar values using the \(\operatorname{Tr}(\cdot)\) operator.

For large values of \(d\), we can approximate by expansion in terms of Hilbert space dimension :

\[\mathbb{Q}:=\mathbb{E}[f_{1}f_{2}]=\frac{Q^{[2]}}{d^{2}}+\frac{Q^{[3]}}{d^{3}} +\cdots:=\hat{\mathbb{Q}}+\hat{\mathbb{Q}}+\mathcal{O}(\frac{1}{d^{4}}) \tag{11}\]

and

\[\mathbb{V}:=\mathbb{E}[f_{1}f_{2}f_{3}f_{4}]=\frac{V^{[4]}}{d^{4}}+\frac{V^{[5 ]}}{d^{5}}+\cdots:=\hat{\mathbb{V}}+\hat{\mathbb{V}}+\mathcal{O}(\frac{1}{d^{5 }}) \tag{12}\]

such that for example:

\[\mathbb{E}[f_{i_{1},\alpha_{1}}f_{i_{2},\alpha_{2}}]|_{\text{conn.}}=\frac{1}{ d^{2}}[\operatorname{Tr}(\rho_{\alpha_{1}}\rho_{\alpha_{2}})\operatorname{Tr}(O_{i_{1 }}O_{i_{2}})]+\mathcal{O}(\frac{1}{d^{3}}) \tag{13}\]

To examine the impact of the number of layers, denoted as \(L\), we can shift our perspective back to the

Figure 2: In variational quantum circuits, as a model for quantum neural networks (QNNs), the number of qubits, denoted as \(n\) and defined as \(n=\log(d)\), can be interpreted as the width of the model. To investigate the impact of the width (\(n\)) and depth (\(L\)) on the output distribution, we conducted simulations of a quantum circuit following the described structure. Each layer of the circuit consists of random two-qubit gates followed by entanglement gates. The observables used in the simulations were constructed from the Pauli \(Z\) matrix, forming a complete basis. We generated \(10^{5}\) random Haar instances for each circuit with given values of \(n\) and \(L\). The joint distribution, \(P(f_{i},f_{j})\), of the output was plotted for two fixed indices \(i\) and \(j\). The input of the circuit was constructed from the MNIST dataset, reduced size with the Principal Component Analysis(PAC), and then mapped to the feature space using the ZZZFeaturemap[5]. As observed, as the quantum neural network (QNN) becomes deeper, the output distribution tends to converge more towards a Gaussian distribution.
parametric space. In the classical context, it has been observed that \(\mathbb{E}[f_{1}f_{2}]\) can be approximated by averaging over \(\mathbb{E}[\langle\nabla\mathbf{\theta}f,\nabla\mathbf{\theta}f\rangle]\), based on Bochner's theorem. This concept has been extended to the quantum scenario in [11; 13]. For the structure like Eq.1, if we define:

\[[\mathbb{H}_{ij,\alpha\beta}]_{\mu\nu}:=\frac{\partial f_{i,\alpha}}{\partial \theta_{\mu}}\frac{\partial f_{j,\beta}}{\partial\theta_{\nu}} \tag{14}\]

then the evolution of the output can be described by the following differential equation:

\[df_{i,\alpha}=-\sum_{\mu,\nu}\eta^{\mu\nu}[\mathbb{H}_{ij,\alpha}]_{\mu\nu} \nabla_{f_{j,\alpha}}\mathcal{L}+\mathcal{O}(\eta^{2}). \tag{15}\]

By defining the left and right operators as \(U_{L,\mu}=\prod_{\ell=1}^{\mu-1}W_{\ell}U_{\ell}\), \(U_{R,\mu}=\prod_{\ell=\mu+1}^{L}W_{\ell}U_{\ell}\), \(V_{L,\mu}=U_{L,\mu}W_{\mu}U_{\mu}\), and \(V_{R,\mu}=U_{R,\mu}\), the derivative of the output with respect to a specific parameter (or layer), expressed in a product form, can be represented as follows:

\[[\mathbb{H}_{ij,\alpha\beta}]_{\mu\nu} =\operatorname{Tr}\Bigl{(}U_{R,\mu}^{\dagger}[X_{\mu},U_{\mu}^{ \dagger}W_{\mu}^{\dagger}U_{L,\mu}^{\dagger}O_{i}U_{L,\mu}W_{\mu}U_{\mu}]U_{R,\mu}\rho_{\alpha}\Bigr{)} \tag{16}\] \[\times\operatorname{Tr}\Bigl{(}U_{R,\nu}^{\dagger}[X_{\nu},U_{ \nu}^{\dagger}W_{\nu}^{\dagger}U_{L,\nu}^{\dagger}O_{j}U_{L,\nu}W_{\nu}U_{\nu }]U_{R,\nu}\rho_{\beta}\Bigr{)}\]

Now, by taking the Haar average over all random unitaries, utilizing Equation 7, We obtain

\[\mathbb{E}[\mathbb{H}_{ij,\alpha\beta}]_{\mu\nu}=\frac{2d}{(d-1) (d+1)(d^{2}+d)}\times\] \[\times\operatorname{Tr}(X_{\mu}X_{\nu})=\frac{c}{d^{2}} \operatorname{Tr}(O_{i}O_{j})\operatorname{Tr}(\rho_{\alpha}\rho_{\beta})+ \mathcal{O}(\frac{1}{d^{3}})\]

where we have defined \(c:=\frac{\operatorname{Tr}(X_{\mu}X_{\nu})}{d}\), which is typically of order one (\(\mathcal{O}(1)\)). Based on certain considerations [11], we can assume that the kernel can be approximated by a frozen one and that the variables \(t\) are independent in the later stages of training. For simplicity, let's assume \(\eta_{\mu\nu}=\delta_{\mu\nu}\). Then the sum over all parameter sizes \(L\) arises as follows:

\[\sum_{\mu\nu}\eta^{\mu\nu}[H_{ij,\alpha\beta}]_{\mu\nu}\sim\frac{L}{d^{2}} \operatorname{Tr}(O_{i}O_{j})\operatorname{Tr}(\rho_{\alpha}\rho_{\beta})+ \mathcal{O}(\frac{1}{d^{3}}) \tag{18}\]

This result motivates us to seek a similar property for the quantum neural network as well, akin to what is observed in classical neural networks:

\[\mathbb{Q}(\rho_{x},\rho_{x^{\prime}}) =\mathbb{E}_{\theta\sim p(\theta)}[f_{\theta}(\rho_{x})\cdot f_{ \theta}(\rho_{x^{\prime}})] \tag{19}\] \[\sim\mathbb{E}_{\theta\sim p(\theta)}[(\nabla_{\theta}f_{\theta} (\rho_{x}),\nabla_{\theta}f_{\theta}(\rho_{x^{\prime}}))]=\mathbb{Q}_{\rm NTK}\]

Thus, as we will explore in the subsequent section, the impact of depth, denoted as \(L\), becomes evident through the faster convergence to a Gaussian distribution:

\[f(t)-f(0)\propto(e^{-L\hat{Q}})^{t\eta} \tag{20}\]

## II Bayesian Learning or How to Get rid of Parameters

For a given quantum model, denoted as \(\mathcal{M}\), the circuit architecture associated with it will possess a collection of parameters that may either remain fixed or variable. These variable parameters, referred to as \(\mathbf{\theta}=\{\theta_{i},i\in[1,\cdots,L]\}\) need to be tuned such that the quantum cir

Figure 3: The phase transition of the energy landscape of a quantum neural network (QNN) occurs when transitioning from the underparameterized regime to the overparameterized regime. In the underparameterized regime, where the number of parameters in the QNN is low compared to the dimension of the Hilbert space (or the degrees of freedom), the QNN’s highly non-convex loss landscape contains unfavorable local minima with eigenvalues of similar magnitudes as the global minimum. However, as the number of parameters increases, the landscape becomes simpler and exhibits numerous favorable local minima. In this scenario, training and converging towards the actual global minimum become much easier. b) The effect of overparameterization in quantum neural networks (QNNs) is the emergence of Gaussian-like output distributions. As the QNN becomes overparameterized, meaning that the number of parameters exceeds the necessary amount to fit the data, the output distribution of the QNN tends to exhibit Gaussian characteristics. This observation suggests that overparameterization promotes a smoother and more concentrated output distribution, resembling a Gaussian distribution.
cuit performs the desired output. Since the variational parameters are assumed to be independent of each other, the prior distribution of these parameters is merely the product of the prior distribution of each individual parameter such that \(p(\mathbf{\theta}|\mathcal{M})=\prod_{\ell=1}^{L}p(\theta_{\ell}|\mathcal{M})\).

Recent investigations have delved into how the prior distribution of the \(p(\mathbf{\theta}|\mathcal{M})\) as a method of initialization strategies of quantum neural networks, such as Gaussian initialization and reduced-domain parameter initialization [23]. These studies suggest that such strategies can enhance convergence speed and alleviate the problem of barren plateaus. However, it should be noted that these assumptions do not guarantee that the global minimum can be reached within a reasonable number of iterations using these initialization methods, and there is a risk of getting trapped in undesirable local minima.

### Marginalization of Parameters

When considering a specific model \(\mathcal{M}\), our prior hypotheses regarding the behavior of the observables we seek to learn are taken into account. The probability distribution describing our prior belief and prediction of the outcome of a quantum circuit, denoted as \(p(\mathbb{P}|\mathbf{\theta},\mathcal{M})\), can be expressed as follows:

\[p(f_{\mathbb{P}}|\mathcal{M})=\int\prod_{i}^{L}d\theta_{i}\,p(f_{\mathbb{P}}| \theta_{i},\mathcal{M})\,p(\theta_{i}|\mathcal{M}) \tag{21}\]

Upon observing the data points of interest \(\mathbb{O}\), we can update our beliefs and make predictions for new data points \(\mathbb{P}\) using the Bayes' rule:

\[p(f_{\mathbb{P}}|f_{\mathbb{O}},\mathcal{M})=\int\prod_{i}^{L}d\theta_{i}\,p(f _{\mathbb{O}}|\theta,\mathcal{M})\,p(\theta|f_{\mathbb{O}},\mathcal{M}) \tag{22}\]

However, this equation can be simplified by integrating out the parameters, resulting in the following equation [19]:

\[p(f_{\mathbb{P}}|f_{\mathbb{O}},\mathcal{M})=\frac{p(f_{\mathbb{P}},f_{ \mathbb{O}}|\mathcal{M})}{p(f_{\mathbb{O}}|\mathcal{M})} \tag{23}\]

This formula allows us to update our beliefs and make predictions based on the observed data points without explicitly involving any parameters. In fact, we have bypassed the intermediate step of updating our posterior on the parameters, which can be obtained using different schemes such as Maximum A Posteriori (MAP) or Maximum Likelihood Estimation (MLE). For example, The optimal parameter values with MAP can be found by maximizing the posterior probability \(\mathbf{\theta}^{*}=\operatorname*{arg\,max}_{\theta}p(\theta|\mathbb{O},\mathcal{ M})\).

## III Probabilistic Model of QNN

In this section, we present a proposition for a probabilistic approach to quantum neural networks. In this novel viewpoint of quantum neural networks, the output of the circuit is represented as a probability distribution, contrasting the traditional approach where the circuit's output is a delta function across all potential outcomes. The key advantage of this perspective lies in its ability to incorporate uncertainty within the framework, enabling increased flexibility in noise analysis. Furthermore, owing to the inherent probabilistic nature of quantum circuits, the variance from the expectation becomes more pronounced, underscoring its significance.

In general, assuming that the output distribution is symmetric with respect to the origin and neglecting odd-numbered correlation functions, we can describe the output distribution using Equation 3:

In a general framework, assuming symmetry around the origin in the output distribution and excluding correlations of odd order, we can characterize the output distribution as [19]

\[p(f)=\frac{e^{-S(f)}}{\mathcal{Z}}, \tag{24}\]

where the action \(S(\cdot)\) is defined as:

\[\begin{split} S(f)&=\frac{1}{2}\sum_{\mu,\nu=1}^{n} K^{\mu\nu}f_{\mu}f_{\nu}\\ &+\sum_{m=2}^{\Lambda}\frac{1}{(2m)!}\sum_{\mu_{1},\cdots\mu_{2m}= 1}C^{\mu_{1}\cdots\mu_{2m}}f_{\mu_{1}}\cdots f_{\mu_{2m}}.\end{split} \tag{25}\]

In this formulation, \(C^{\mu_{1}\cdots\mu_{2m}}\) represents the \(m\)-th order couplings, and \(\Lambda\) serves as the truncation threshold for limiting the series expansion up to the \(\Lambda\)-th order. Furthermore, the partition function \(\mathcal{Z}\) in the denominator of Equation 3 is defined as:

\[\mathcal{Z}=\int d^{n}fe^{-S(f)} \tag{26}\]

Here, \(C^{\mu_{1}\cdots\mu_{2m}}\) represents the \(m\)-th order couplings defined as

\[\mathbb{E}[f_{1}f_{2}\cdots f_{m}]=\int df^{m}p(f)f_{1}f_{2}\cdots f_{m}, \tag{27}\]

and \(\Lambda\) is the limit for truncating the series to the \(\Lambda\)-th order. The partition function, \(\mathcal{Z}\) in the denominator of Equation 3, is given by:

\[\mathcal{Z}=\int d^{n}fe^{-\mathcal{S}(f)} \tag{28}\]

As the number of qubits, denoted by \(n\), increases, the Hilbert space dimension expands exponentially as \(d=2^{n}\). In the asymptotic regime where \(d\) approaches infinity, higher-order corrections within the correlation
relation become negligible. Specifically, we observe that \(\mathbb{E}[f_{1}\cdots f_{m}]\sim\mathcal{O}(1/d^{m^{\prime}+2})\) for \(m^{\prime}>0\). This implies that the action can be perturbed with respect to the dimension of the Hilbert space:

\[S=S^{(0)}+\frac{S^{(2)}}{d^{2}}+\frac{S^{(4)}}{d^{4}}+\mathcal{O}\left(\frac{1}{ d^{6}}\right). \tag{29}\]

This perturbation approach bears resemblance to the \(\frac{1}{d}\) expansion technique employed in the domains of statistical mechanics and quantum field theory [19], where \(S^{(0)}\) characterizes the vacuum state of the theory and \(S^{(2)}\) represents the first-order interaction and so on.

### Infinite Width QNN

In this section, we want to study the extreme over-parameterized case \(d\rightarrow\infty\) where all higher order terms than quadratic interaction can be neglected: \(S\approx S^{(0)}+S^{(2)}/d^{2}\). In this regime, it becomes possible to express the output distribution of the quantum neural network as a Gaussian distribution.

Our learning task, aimed at determining the cacheetic of this distribution, entails observing samples from ensembles of data points and subsequently predicting the output of the quantum neural network for new, unseen data points. We represent the observed datapoints and their corresponding sets as \(\mathbb{O}=\{\rho_{\beta},f_{i,\beta}:=y_{\beta,i};\beta\in[1,\cdots|\mathbb{ O}|],i\in[1,\cdots,n]\}\).Here, \(|\mathbb{O}|\) denotes the size of the set of observable operators. In this representation, \(\rho_{\beta}\) denotes the density matrix of the \(\beta\)-th data point from the observed ensemble, while \(y_{i,\beta}:=f_{i,\beta}\) represents the true label of this particular data point. Similarly, we define the prediction sets as \(\mathbb{P}=\{\rho_{\alpha},=f_{j,\alpha};\alpha\in[1,\cdots|\mathbb{P}|],j\in[ 1,\cdots,n]\}\). With these two sets, we can represent our updated posterior as follows:

\[p(f_{\mathbb{O}},f_{\mathbb{P}})=\frac{1}{Z}\exp(-\frac{1}{2}\sum_{i,j=1}^{n} \sum_{\mu_{1},\mu_{2}\in\mathbb{P}\cup\mathbb{O}}K^{ij,\mu_{1}\mu_{2}}f_{i, \mu_{1}}f_{j,\mu_{2}}) \tag{30}\]

where \(Z\) is the normalization factor, given by \(|2\pi K|^{n^{2}}\).

Utilizing Bayes' rule, as presented in Equation 23, we can gain deeper insights into the QNN's output distribution for a particular data point, denoted as \(f_{\mathbb{P}}:=f\), based on previously observed data points, \(f_{\mathbb{O}}\). Similar computations to those used in classical neural networks, as detailed in [19] and outlined in Appendix A, allow us to express this distribution as a Gaussian distribution, or potentially as a Gaussian process \(p(f)=\frac{1}{Z}\exp(S)\) with Gaussian action:

\[p(f) =\frac{1}{Z}\exp(-S),\] \[S =\frac{1}{2}\sum_{\begin{subarray}{c}i,j=1\\ \beta_{1},\beta_{2}\in\mathbb{P}\end{subarray}}^{n}\mathbb{K}^{ij,\beta_{1} \beta_{2}}[f_{i,\beta_{1}}-m_{i,\beta_{1}}^{\infty}][f_{j,\beta_{2}}-m_{j, \beta_{2}}^{\infty}], \tag{31}\]

which is characterized by the following mean and covariance:

\[[m]_{i,\mu} =\sum_{\lambda,\sigma\in\mathbb{O}}K_{ij,\mu\lambda}K^{ij,\lambda \sigma}y_{i,\sigma}:=K(\mathbb{P},\mathbb{O})K^{-1}(\mathbb{O},\mathbb{O}) \mathbf{y},\] \[\mathbb{K}]_{ij,\mu\nu} =K_{ij,\mu\nu}-\sum_{\lambda,\sigma\in\mathbb{O}}K_{ij,\mu\lambda }K^{ij,\lambda\sigma}K_{ij,\sigma\nu}\] \[=K(\mathbb{P},\mathbb{O})-K(\mathbb{P},\mathbb{P})K^{-1}(\mathbb{ O},\mathbb{O})K(\mathbb{O},\mathbb{P}). \tag{32}\]

Here, \(K\) is a two-point correlation function representing the kernel function. Therefore, we can establish the explicit relationship between the mean and covariance in terms of the inputs, denoted by \(\{\rho_{\beta}\}_{\beta\in\mathbb{O}}\) and observables \(\{O\}\):

\[[m]_{i,\beta} =\sum_{\alpha_{1},\alpha_{2}\in\mathbb{O}}\mathbb{E}[f_{i}f_{j}]_{ \beta,\alpha_{1}}\mathbb{E}[f_{i}f_{j}]_{\alpha_{1},\alpha_{2}}^{-1}y_{i,\alpha _{2}}\] \[=\sum_{\alpha_{1},\alpha_{2}\in\mathbb{O}}[\text{Tr}(\rho_{\beta }\rho_{\alpha_{1}})][\text{Tr}(\rho_{\mathbf{z}_{1}}\rho_{\mathbf{z}_{2}})]_{\alpha_{1 }\alpha_{2}}^{-1}y_{i,\alpha_{2}}+\mathcal{O}(\frac{1}{d^{2}}) \tag{33}\]

Furthermore, the covariance can be calculated using the following expression:

\[[\mathbb{K}]_{ij,\beta_{1}\beta_{2}}=\frac{\text{Tr}(O_{i}O_{j}) }{d^{2}}[\text{Tr}(\rho_{\beta_{1}}\rho_{\beta_{2}})\] \[-\sum_{\alpha_{1},\alpha_{2}\in\mathbb{O}}\text{Tr}(\rho_{\beta_{1 }}\rho_{\beta_{2}})[\text{Tr}(\rho_{\mathbf{z}_{1}}\rho_{\mathbf{z}_{2}})]_{\beta_{2} \alpha_{2}}^{-1}\,\text{Tr}(\rho_{\alpha_{2}}\rho_{\beta_{2}})]+\mathcal{O}( \frac{1}{d^{4}}). \tag{34}\]

In literature, the term \(\text{Tr}(\rho_{x}\rho_{x})\) is occasionally denoted as the Quantum fidelity kernel \(K^{\mathbb{F}}(x,x^{\prime})\)[6]. Using this notation, we can make the previous equation more understandable:

\[\mathbb{K}_{ij,\mathbb{P}}=\frac{\text{Tr}(O_{i}O_{j})}{d^{2}}[K^{\mathbb{P}}_{ \mathbb{P},\mathbb{O}}-K^{\mathbb{P}}_{\mathbb{P},\mathbb{P}}(K^{\mathbb{P}})_{ \circ,\circ}^{-1}K^{\mathbb{P}}_{\mathbb{O},\mathbb{P}}]+\mathcal{O}(\frac{1}{d ^{4}}) \tag{35}\]

### Parameter space representation

Within the context of parameter space, our learning endeavor can be seen as the creation of a parametric model \(f_{\text{QNN}}(\rho;\mathbf{\theta})\), accompanied by a distribution \(p(\mathbf{\theta})\). The aim is to locate the best value of \(\mathbf{\theta}=\mathbf{\theta}^{*}\), ensuring minimal error considering available computational resources (like the number of qubits, gates, and observed datapoints). The goal is to have our model closely approximate the actual distribution and align accurately with the real datapoints \(f_{\text{QNN}}(x;\mathbf{\theta}^{*})\approx f(x)\). The aspiration is that, with proper initialization near the actual global minimum, we can locate this global minimum using gradient descent methods. As quantum neural networks lean towards overparameterization, this assumption becomes more credible. Given this assumption, the
output of the Quantum Neural Network (QNN) can be approximated using a Taylor series:

\[\begin{split}& f_{\alpha}^{i,(t)}=f_{\alpha}^{i,(0)}-\eta\sum_{ij, \beta}\mathbb{H}^{ij,\beta\alpha}\epsilon_{j,\beta}\leftrightarrow\mathbb{Q}_{ \text{NTK}}\\ &+\frac{\eta^{2}}{2}\sum_{ij_{1}j_{2},\beta_{1}\beta_{2}}d\mathbb{ H}^{ij_{1}j_{2},\alpha_{\beta}\beta_{2}}\epsilon_{j_{1},\beta_{1}}\epsilon_{j_{2}, \beta_{2}}\leftrightarrow\text{d}\mathbb{Q}_{\text{NTK}}\\ &-\frac{\eta^{3}}{6}\sum_{ij_{1}\cdots\alpha_{3}}d\mathbb{H}^{ij _{1}\cdots\alpha_{3}}\epsilon_{j_{1},\beta_{1}}\epsilon_{j_{2},\beta_{2}} \epsilon_{j_{3}\beta_{3}}\leftrightarrow\text{d}\mathbb{Q}_{\text{NTK}}\\ &+\mathcal{O}(\frac{1}{d^{2}})\end{split} \tag{36}\]

where \(\epsilon=\nabla_{f}\mathcal{L}\) denotes the derivative of the loss function with respect to the output function and \(\mathbb{H}\) is defined as:

\[\mathbb{H}_{ij,\alpha\beta}=\sum_{\mu\nu}\eta_{\mu\nu}\frac{\partial f_{i, \alpha}}{\partial\theta_{\mu}}\frac{\partial f_{j,\beta}}{\partial\theta_{ \nu}}. \tag{37}\]

As the width of the Quantum Neural Network (QNN) approaches infinity, \(\mathbb{H}\) can be approximated by a frozen kernel, referred to as the Quantum Neural Tangent Kernel (QNTK):

\[\mathbb{H}_{ij,\alpha\beta}=\mathbb{Q}_{ij,\alpha\beta}+\mathcal{O}(\frac{1}{d ^{2}}) \tag{38}\]

By comparing the expansion of the parameter space with the functional method outlined in Sec.III, we can draw the following analogy:

\[\begin{split}\mathbb{Q}_{\text{NTK}}:&\mathbb{E}[ \frac{\partial f_{i_{1},\alpha_{1}}}{\partial\theta_{\mu}}\frac{\partial f_{j,\alpha_{2}}}{\partial\theta_{\nu}}]\sim\frac{1}{d^{2}}\operatorname{Tr}( \rho_{\alpha_{1}}\rho_{\alpha_{2}})\operatorname{Tr}(O_{i}O_{j})\\ &\sim\mathbb{E}[f_{i_{1},\alpha_{1}}f_{j,\alpha_{2}}]=K_{ij, \alpha_{1}\alpha_{2}}\end{split} \tag{39}\]

Likewise, for the other terms in the Eq.36 expansion, these are known as quantum meta kernels, which correspond to the third and fourth-order correlations:

\[\begin{split}\text{d}\mathbb{Q}_{\text{NTK}}:&\mathbb{ E}[\frac{\partial^{2}f_{i_{1},\alpha_{1}}}{\partial\theta_{\mu}\partial\theta_{ \nu}}\frac{\partial f_{i_{2},\alpha_{2}}}{\partial\theta_{\lambda}}\frac{ \partial f_{i_{3},\alpha_{3}}}{\partial\theta_{\sigma}}]=0\\ &\sim\mathbb{E}[f_{i_{1},\alpha_{1}}f_{i_{2},\alpha_{2}}f_{i_{3}, \alpha_{3}}],\end{split} \tag{40}\]

and

\[\begin{split}\text{d}\mathbb{Q}_{\text{NTK}}:&\mathbb{ E}[\frac{\partial^{3}f}{\partial\theta_{1}\partial\theta_{2}\partial\theta_{3}} \frac{\partial f}{\partial\theta_{4}}\frac{\partial f}{\partial\theta_{5}} \frac{\partial f}{\partial\theta_{6}}]\\ &+\mathbb{E}[\frac{\partial^{2}f}{\partial\theta_{1}\partial \theta_{2}}\frac{\partial^{2}f}{\partial\theta_{3}\partial\theta_{4}}\frac{ \partial f}{\partial\theta_{5}}\frac{\partial f}{\partial\theta_{6}}]\\ &\leftrightarrow\mathbb{E}[f_{i_{1},\alpha_{1}}f_{i_{2},\alpha_{2} }f_{i_{3},\alpha_{3}}f_{i_{4},\alpha_{4}}]=\mathbb{V}_{i_{1}i_{2}i_{3}i_{4}, \alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}\end{split} \tag{41}\]

Solving Eq.36 entails dealing with a complex set of differential equations, which generally do not yield a closed-form solution [19]. However, in the limit of \(d\rightarrow\infty\), it is acceptable to only retain the linear term. The solution to the differential equation, in terms of steps (represented here by the continuous variable \(t\)), can be expressed as:

\[\begin{split} f_{i,\alpha}^{(t)}&\simeq m_{i,\alpha} ^{(t)}+\Sigma_{i,\alpha}^{(t)}=\sum_{\beta_{1},\beta_{2}\in\mathbb{O}}\mathbb{Q }_{\alpha\beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}(1-e^{-\eta\mathbb{Q}t})y_{i,\beta_{2}}\\ &+f_{i,\alpha}^{(0)}-\sum_{\beta_{1},\beta_{2}\in\mathbb{O}} \mathbb{Q}_{\alpha\beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}(1-e^{-\eta \mathbb{Q}t})f_{i,\beta_{2}}^{(0)}.\end{split} \tag{42}\]

In the limit of large iteration times, we can determine the output of the circuit (without actually performing any training) as a linear kernel regression:

\[f_{\text{QNN}}(\rho_{x})\approx f_{\mathbb{Q}_{\text{NTK}}}=\mathbb{Q}_{\text {NTK}}(\rho_{x},\rho_{X})^{T}\cdot\mathbb{Q}_{\text{NTK}}(\rho_{X},\rho_{X})^{ -1}\cdot Y, \tag{43}\]

and the optimized parameters at the end of training would be determined as

\[\mathbf{\theta}^{*}=\mathbf{\theta}^{(0)}-\eta\nabla_{\mathbf{\theta}}f\cdot\mathbb{Q}_{ \text{NTK}}^{-1}(\mathbf{f}-\mathbf{y}) \tag{44}\]

To see the connection with the Gaussian process, let's consider the expected output for a particular input \(\rho_{\alpha}\), which can be given as:

\[\mathbb{E}[f_{i,\alpha}^{(t)}]=\sum_{\beta_{1},\beta_{2}\in\mathbb{O}}\mathbb{ Q}_{\alpha\beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}(1-e^{-\eta\mathbb{Q}t})y_{i,\beta_{2}}, \tag{45}\]

which, in the limit of infinite iterations, converges to \(m_{i,\alpha}^{(\infty)}:=\sum_{\beta_{1},\beta_{2}\in\mathbb{O}}\mathbb{Q}_{ \alpha\beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}y_{i,\beta_{2}}\).

Likewise, the covariance of the circuit's output will take the following form:

\[\begin{split}&\mathbb{E}[f_{i,\alpha_{1}}^{(t)}f_{j,\alpha_{2}}^{(t)}] |_{\text{connected}}=\mathbb{Q}_{\alpha_{1}\alpha_{2}}\\ &-\sum_{\beta_{1},\beta_{2}\in\mathbb{P}}\mathbb{Q}_{\alpha_{2} \beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}(1-e^{-\eta\mathbb{Q}t})\mathbb{Q}_{ \alpha_{1}\beta_{2}}\\ &-\sum_{\beta_{1},\beta_{2}\in\mathbb{P}}\mathbb{Q}_{\alpha_{1} \beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}(1-e^{-\eta\mathbb{Q}t})\mathbb{Q}_{ \alpha_{2}\beta_{2}}\\ &+\sum_{\beta_{1},\cdots,\beta_{4}\in\mathbb{P}}\mathbb{Q}_{\alpha_{1} \beta_{1}}\mathbb{Q}^{\beta_{1}\beta_{2}}[1-e^{-\eta\mathbb{Q}t}]\mathbb{Q}_{ \alpha_{2}\beta_{3}}[1-e^{-\eta\mathbb{Q}t}]\\ &\times\mathbb{Q}^{\beta_{3}\beta_{4}}\mathbb{Q}_{\beta_{4}\beta_{2}}: =\Sigma_{\alpha_{1}\alpha_{2}}^{ij,t}.\end{split} \tag{46}\]

The higher order moments are of the order \(\mathcal{O}(1/d^{2})\) and can be disregarded. Based on this observation, we can infer that the most likely distribution to describe the statistics of the output is a Gaussian process. Upon comparing with the functional description of the output distribution, we observe that we obtain the same distribution as predicted by the Gaussian process, particularly as \(t\rightarrow\infty\), and substituting \(\mathbb{Q}\) with \(\mathbb{K}\):

\[f^{(t)}\sim\mathcal{GP}(m^{(t)},\Sigma^{(t)})\overset{\text{lim}}{\to }\mathcal{GP}(m^{(\infty)},\mathbb{K}) \tag{47}\]

comparing the linear kernel regression method with the Gaussian process,

\[\begin{split}& f\sim f_{\mathbb{Q}_{\text{NTK}}}\leftrightarrow\text{ Kernel Regression}\\ & p(f)\sim\mathcal{N}(\mu,\Sigma)\leftrightarrow\text{Gaussian Process}\end{split} \tag{48}\]
### Finite Width QNN

When working within the finite width regime of the parameter space, it becomes necessary to consider higher-order terms, such as dQNTK and ddQNTK, which are based on the higher-order derivative terms dH and dHH. Using traditional methods to obtain higher-order corrections in the output distribution, such as for dQNTK, can be unfeasible and not provide a closed-form solution. However, in this section, we demonstrate that by utilizing Bayesian learning, we can obtain higher-order corrections of the output distribution in a closed form, at least for the case of dQNTK.

To enhance the model's ability to capture the finite width regime, it becomes necessary to go beyond the Gaussian distribution and utilize representation learning techniques. By incorporating additional correlation functions into the model, it becomes better equipped to accurately represent the underlying reality.

To obtain a nearly-Gaussian (NG) distribution, comparable to ddQNTK, we need to incorporate the next leading order in finite width, which is of \(O(1/d)\). This is because in the nearly-Gaussian distribution, the odd terms are zero or equivalently, dQNTK is zero.

\[S_{\text{NG}}:=\frac{1}{2}\sum_{\mu,\nu}^{n}K^{\mu\nu}f_{\mu}f_{\nu}-\frac{ \lambda}{4!}\sum_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}=1}^{n}V^{\nu_{1}\nu_{2}\nu_{3} \nu_{4}}f_{\nu_{1}}f_{\nu_{2}}f_{\nu_{3}}f_{\nu_{4}}. \tag{49}\]

such that

\[\begin{split} K^{\alpha_{1}\alpha_{2}}&=\hat{ \mathbb{Q}}^{\alpha_{1}\alpha_{2}}+\mathcal{O}(\frac{1}{d^{3}})\\ V^{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}&= \hat{\mathbb{V}}^{(\alpha_{1}\alpha_{2})(\alpha_{3}\alpha_{4})}+\mathcal{O}( \frac{1}{d^{5}})\end{split} \tag{50}\]

By employing the same calculations as we did for the infinite width scenario, we can determine the contributions of finite width structures to the mean and covariance:

\[\begin{split}&\mathbb{E}[f_{i,\beta}]=m_{i,\beta}^{\infty}+\frac{ 1}{3!}\sum_{j_{1},\beta\in\mathbb{P}}K_{ij_{1},\beta\beta_{1}}\\ &\times[\sum_{\begin{subarray}{c}j_{2}j3j_{4}\\ \alpha_{2}\alpha_{3}\alpha_{4}\in\mathbb{O}\end{subarray}}V^{ji_{1}j_{2}j_{3} j_{4},\beta_{2}\alpha_{3}\alpha_{4}}y_{j_{2},\alpha_{2}}y_{j_{3},\alpha_{3}}y_{j_{4}, \alpha_{4}}]\\ &\sum_{\begin{subarray}{c}j_{1}j_{2}j_{3}\\ \beta_{1}\beta_{2}\beta_{3}\in\mathbb{P}\end{subarray}}[K_{ij_{1},\beta_{ \beta}\beta_{1}}K_{j_{2}j_{3},\beta_{2}\beta_{3}}+K_{ij_{2},\beta\beta_{2}}K_{ ij_{3},\beta_{1}\beta_{3}}\\ &+K_{ij_{3},\beta\beta_{3}}K_{j_{1}j_{2},\beta_{1}\beta_{2}}] \sum_{j_{4},\alpha\in\mathbb{O}}V^{ji_{1}j_{2}j_{3}j_{4},\beta_{1}\beta_{2} \beta_{3}\alpha}y_{j_{4},\alpha}\\ &+\mathcal{O}(\frac{1}{d^{4}}),\end{split} \tag{51}\]

and

\[\begin{split}&\Sigma^{\prime}=\mathbb{E}[f_{i_{1},\mu_{1}}f_{i_{2}, \mu_{2}}]|_{\text{conn.}}=K_{i_{1}i_{2},\mu_{1}\mu_{2}}\\ &-\frac{1}{2}\sum_{\begin{subarray}{c}j_{1}j_{2}j_{3}j_{4}=1\\ j_{1}\nu_{2}\nu_{3}\nu_{4}\end{subarray}}V^{ji_{1}j_{2}j_{3}j_{4},\nu_{1}\nu _{2}\nu_{3}\nu_{4}}K_{i_{1}j_{1},\mu_{1}\nu_{1}}K_{i_{2}j_{2},\mu_{2}\nu_{2}},\\ &\times K_{j_{3}j_{4},\nu_{3}\nu_{4}}\end{split} \tag{52}\]

where \(m_{i,\alpha}^{(\infty)}=\sum_{\beta_{1},\beta_{2}}K_{\alpha_{1}\beta_{1}}K^{ \beta_{1}\beta_{2}}y_{i,\beta_{2}}\) thus e can refine the kernel regression method by incorporating finite-size corrections:

\[m^{\prime}=\bar{f}=m^{(\infty)}+m^{(1/d)}+\mathcal{O}(\frac{1}{d^{2}}) \tag{53}\]

If we aim to maintain the Gaussian approximation, the corrected covariance, which arises from higher orders, and the effects of wiring with others, can be described as

\[f\approx\mathcal{GP}(m^{\prime},\Sigma^{\prime}). \tag{54}\]

However, and in the large \(d\) limit and the output distribution can be approximate better by forth order correlator as :

\[p(f)=p(f|y_{\alpha},\hat{\mathbb{Q}}_{\mu\nu},\hat{\mathbb{V}}_{\mu\nu\delta \sigma}). \tag{55}\]

such that \(\hat{\mathbb{Q}}=\Sigma^{\prime}\) and \(\hat{\mathbb{V}}\) given by (see appendix for more details)

\[\begin{split}&\mathbb{V}_{i_{1}i_{2}i_{3}i_{4},\mu_{1}\mu_{2}\mu_{3} \mu_{4}}=\sum_{j_{1}j_{2}j_{3}j_{4}}\sum_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}V^{ji_{1 }j_{2}j_{3}j_{4},\nu_{1}\nu_{2}\nu_{3}\nu_{4}}\\ &\times K_{i_{1}j_{1},\mu_{1}\nu_{1}}K_{i_{2}j_{2},\mu_{2}\nu_{2}} K_{i_{3}j_{3},\mu_{3}\nu_{3}}K_{i_{4}j_{4},\mu_{4}\nu_{4}}+\mathcal{O}(\frac{1}{d^{5}}) \end{split} \tag{56}\]

## IV Conclusion

In this study, we studied the variational quantum circuits, which serve as a model for quantum neural networks from the sense of function space, in stead of parmart space. We explored the impact of over-parameterization through the expansion of width and depth. We observed theoretically and by simulation that this led to a shift in the output distribution of the quantum models towards a Gaussian process. We also examined how the finite size of the width can skew the distribution towards non-Gaussianity. We demonstrated how this can be achieved using a \(1/d\) expansion in closed form relation, setting the stage for systematic consideration of higher-order corrections.

Note added: This work was presented informally to the scientific community during the APS March Meeting 2023, and toward the completion of this work, we became aware of work with a similar topic on arXiv [4].