Cautious optimism for deep parameterized quantum circuits

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Cautious optimism for deep parameterized quantum circuits

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Cautious optimism for deep parameterized quantum circuits Marie Kempkes,1, 2, ∗ Elies Gil-Fuster,3, 4 Carlos Bravo-Prieto,3 Aroosa Ijaz,5, 6 Alissa Wilms,3, 7 Jens Eisert,3 Evert van Nieuwenburg,1 and Vedran Dunjko1 1 Leiden University, Niels Bohrweg 1, 2333 CA Leiden, Netherlands Volkswagen Group Innovation, Berliner Ring 2, 38440 Wolfsburg, Germany 3 Dahlem Center for Complex Quantum Systems, Freie Universität Berlin, 14195 Berlin, Germany 4 Fraunhofer Heinrich Hertz Institute, 10587 Berlin, Germany 5 Department of Physics and Astronomy, University of Waterloo, ON N2L 3G1, Canada 6 Vector Institute, Toronto, ON M5G 0C6, Canada 7 Porsche Digital GmbH, 71636 Ludwigsburg, Germany

arXiv:2607.21409v1 [quant-ph] 23 Jul 2026

A central challenge in quantum machine learning is understanding the scaling behavior of parameterized quantum circuits (PQCs). In particular, it remains unclear how their performance on unseen data changes as the number of trainable parameters increases. Prior works have derived formal generalization guarantees for quantum models, but it is well-known that many such results do not fully characterize generalization behavior in practice. In this work, we show that gradient-based PQCs can exhibit improved performance on unseen data as model size increases, displaying the phenomenon of double descent. This contrasts with the traditional view that larger models lead to degraded generalization. We provide analytical results rigorously underpinning this behavior by leveraging add-one-in perturbation techniques and spectral properties of random matrices. We support these results with numerical experiments on re-uploading PQCs across several data sets and training set sizes, consistently observing the predicted double descent behavior. While other obstacles on the path toward practical quantum machine learning remain, our finding that deeper parameterized quantum circuits do not necessarily exhibit degraded performance provides reasons for cautious optimism.

Quantum machine learning (QML) has emerged as a promising paradigm at the intersection of quantum computing and data-driven learning, with the potential to outperform classical methods for specific learning tasks. Indeed, a growing body of work has identified concrete settings in which quantum models can offer computational advantages over classical counterparts [1–11]. Despite this substantial progress, many fundamental aspects of QML remain poorly understood. In particular, it is still unclear whether, and under what precise conditions, quantum advantage can be expected, especially as models are scaled to larger numbers of trainable parameters. While several limitations to the trainability of parameterized quantum circuits (PQCs) are now well characterized [12–16], considerable effort has also been devoted to understanding their generalization behavior through bounds on the difference between empirical and expected risk. These bounds depend, for example, on the training set size, the number of trainable gates, and the data-encoding strategy [17–24]. However, they do not directly address how the expected risk attained by a trained model changes as the number of parameters crosses the so-called interpolation threshold, namely, the critical point at which the number of parameters becomes comparable to the number of training points. In classical machine learning, studying the behavior of models across this threshold has led to substantial insights through the so-called double descent phenomenon [25]. Contrary to traditional expectations from statistical learning theory, double descent describes a non-monotonic relationship between test error and model size. As the number of parameters increases, the test error initially follows a classical Ushaped curve, peaking near the interpolation threshold. Remarkably, further overparameterization often leads to a second descent phase, during which the test error decreases again

and eventually saturates. This behavior marks a transition from the underparameterized to the overparameterized regime and has been observed across a wide range of classical models and learning tasks [25–28]. In the realm of quantum machine learning, however, analogous scaling behavior has so far been much less explored. Ref. [29] has shown that quantum kernel methods can, in principle, exhibit a similar scaling behavior in the overparameterized regime, assuming trainability can be maintained. More recently, Ref. [30] has reported non-monotonic test-error dynamics during the training of overparameterized PQCs. This naturally raises the question of whether PQCs trained with gradient-based methods also exhibit double descent as a function of the number of trainable parameters. To date, some evidence has pointed toward a standard U-shaped risk curve for PQCs, suggesting that overparameterization may not yield a second descent [31]. This conclusion is obtained for circuits in which concentration effects dominate, a regime often associated with poor trainability. In this work, we instead focus on PQCs operating in a trainable regime and study their generalization properties as the number of trainable parameters increases. Building on an existing lower bound on the expected risk, we show that, under reasonable assumptions, the bound attains a maximum at the interpolation threshold. We further derive a corresponding upper bound and show that it likewise attains a maximum at interpolation. Together, these results establish a double descent behavior in the risk bounds of PQCs analogous to that observed in classical machine learning. We further support our theoretical findings with numerical experiments on a variety of datasets. Overall, our results provide new insight into the scaling behavior of gradient-based quantum models and suggest that overparameterization may play a role in PQCs

2 similar to that observed in classical machine learning, offering cautious optimism for the prospects of quantum machine learning at scale.

PRELIMINARIES

We consider the ubiquitous setting of supervised learning with quantum learning models based on parametrized quantum circuits (PQCs). Refer to Appendix A for a fully-formal presentation of our framework. We consider a learning task on d-dimensional inputs x and K-dimensional labels y, and assume we are given a training set S = {(xi , yi )}N i=1 of N input-output pairs, sampled i.i.d. from the underlying distribution that defines the problem. Consider data re-uploading PQCs [32], where datadependent layers Uℓ (x) and trainable layers Vℓ (ϑ) are alternated, for ℓ ∈ {1, . . . , L}. Here, ϑ is a p-dimensional vector of trainable parameters. Starting from an initial state, the data re-uploading PQC prepares a parametrized state ρ(x; ϑ) that depends on both the data x and the parameters ϑ. To fully specify our model, we fix observables Ok for k ∈ {1, . . . , K}, from which we define the hypothesis functions. Specifically, the k th output of the hypothesis function is given by the expectation value of Ok with respect to the parametrized state ρ(x; ϑ) [fϑ (x)]k = Tr{ρ(x; ϑ)Ok } .

(1)

Given a loss function ℓ((x, y); ϑ), our goal is to minimize the expected risk, defined as the average loss over the problem distribution L(ϑ) = E[ℓ((x, y); ϑ)].

(2)

In practice, we assume the underlying distribution is unknown, and we explicitly optimize the empirical risk – the average loss over the training set N

L̂S (ϑ) =

1 X ℓ((xi , yi ); ϑ), N i=1

(3)

where (xi , yi ) are the elements of the training set S. In this work, we focus on the mean-squared error loss ℓ((x, y); ϑ) = 1 2 2 ∥fϑ (x) − y∥ . We next follow the formalism in Ref. [33], with further technical details being deferred to Appendix B. Our study is specialized to the final parameters reached during training, which we denote by ϑ̂S . We make the explicit assumption that ϑ̂S is a local minimum of the empirical risk, as is common in gradient-based QML. Under this assumption, ϑ̂S satisfies ∇ϑ L̂S (ϑ̂S ) = 0

and

ĤS (ϑ̂S ) ⪰ 0,

(4)

where ĤS (ϑ) := ∇2ϑ L̂S (ϑ) is the Hessian matrix of the empirical risk.

We now derive an approximation of the expected risk via a so-called add-one-in analysis. Specifically, we consider inserting an additional pair (x′ , y ′ ) to the training set S. Denote by ϑ̂S ′ the final training parameters, where S ′ is the union of S and the new datum. The value of the loss ℓ((x′ , y ′ ); ϑ̂S ′ ) will depend on the new datum. With these, define the so-called add-one-in loss L′ (S) as h i L′ (S) := E ℓ (x′ , y ′ ); ϑ̂S ′ , (5) where the average is over (x′ , y ′ ) being drawn from the problem distribution. While not immediately obvious, the quantity L′ (S) plays a central role in deriving an approximation of the expected risk. The key idea is that the transition from S to S ′ constitutes a small perturbation of the training set, allowing the parameter displacement ϑ̂S ′ − ϑ̂S and the associated change in loss at (x′ , y ′ ) to be approximated via so-called influence functions [34]. Averaging over the choice of the additional sample (x′ , y ′ ) then yields an approximation of the expected risk that is accurate up to corrections of order O(N −2 ). To state the resulting expression, we introduce the gradients of the loss function Jℓ ((x, y); ϑ) = ∇ϑ ℓ((x, y); ϑ) and the gradients of the hypothesis functions Jf ((x, y); ϑ) = ∇ϑ fϑ (x). Denote by C(ϑ) := E [Jℓ ((x, y); ϑ)Jℓ ((x, y); ϑ)⊺ ]

(6)

the average of the outer product of the loss gradients over the problem distribution. Analogously, Cf (ϑ) is given by the average of the outer product of the hypothesis gradients Cf (ϑ) := E [Jf ((x, y); ϑ)Jf ((x, y); ϑ)⊺ ] .

(7)

As a technical comment, we condition the average in Cf (ϑ) to be only over samples (x, y) for which ∥Jf ((x, y); ϑ)∥2 is greater than zero. Note that both C(ϑ) and Cf (ϑ) are uncentered covariance matrices. With these quantities, Theorem 3 in Ref. [33] provides an error decomposition for the expected risk h i−1 1 ′ C(ϑ̂S ) (8) Tr ĤS (ϑ̂S ) L(ϑ̂S ) = L (S)+ N +1 evaluated at the local minimum ϑ̂S , where the inverse of the Hessian is understood as the Moore–Penrose pseudoinverse. The decomposition holds up to second- and higher-order corrections O N −2 . Building on this decomposition, Theorem 4 in Ref. [33] derives a lower bound on the expected risk in the scalar-output case (K = 1) by rewriting the trace term and applying standard eigenvalue inequalities, yielding L(ϑ̂S ) ≥ L′ (S) +

2 1 α ξmin λmin (Cf (ϑ̂S )) . N + 1 λmin (ĤS (ϑ̂S ))

(9)

2 Here, ξmin is the smallest non-zero value of ∥Jf ((x, y); ϑ)∥2 , and the theorem assumes that data can be found for which

3 ∥Jf ((x, y); ϑ)∥2 > 0 with probability α over the problem distribution. Further, λmin denotes the smallest non-zero eigenvalue of a matrix. We note that the lower bound scales inversely with the smallest non-zero eigenvalue of the Hessian matrix of the empirical risk at local minima, λmin (ĤS (ϑ̂S )). The authors of Ref. [33] then carefully investigate the smallest non-zero eigenvalues of the matrices arising in classical neural networks. At this point, we depart from their formalism, and rather specialize to usual gradient-based PQCs.

DOUBLE DESCENT BEHAVIOR IN PQCS

Our goal is to characterize the behavior of the lower bound on the expected risk in Eq. (9) as we vary the number of trainable parameters in our model p, as well as the number of training data N . In particular, we consider the asymptotic limit where p, N → ∞ with fixed ratio p/N → γ. The analysis relies on a series of formal assumptions, listed in Appendix C together with all further technical details. These assumptions are meant to capture the behavior of QML models in practice, when undergoing gradient-based training, and they mostly concern statements about the spectral properties of matrices related to the Hessian and the covariance matrices introduced above. We expect our assumptions to be largely ansatzagnostic and to hold for a broad class of trainable PQCs, including data re-uploading architectures. We further test a selection of these assumptions numerically in Appendix F. We can now state our main result, which holds in the simpler case of learning tasks with a single output dimension K = 1. Theorem 1 (Double descent peak at interpolation – informal). With the definitions above, and under reasonable assumptions, for Gaussian-distributed input data x ∼ N (0, Id ), outputs y ∈ R, and considering the mean-squared error loss, the lower bound on L(ϑ̂S ) attains a maximum at p = N with high probability in the limit p, N → ∞ with fixed ratio p/N → γ. The formal version of this result is Theorem C.3, which can be found in Appendix C. The proof can be found in Appendix D. The main idea for the proof is to carefully control the spectral properties of several random matrices, which we achieve by combining our assumptions with seminal results in random matrix theory. Ultimately, we argue that the smallest non-zero eigenvalue of the dominant contribution to the Hessian matrix follows the behavior predicted by the so-called Marčenko-Pastur law in the stated limit. Our analysis is inspired by the framework of Ref. [33], with suitable modifications to the underlying assumptions. The assumption of i.i.d. Gaussian inputs is standard in highdimensional random matrix analyses and has been extensively studied in classical literature [35, 36]. The earlier work on double descent in QML [29] already argues that such assumptions often extend beyond strictly Gaussian data. Notably, the datasets considered in our numerical experiments in Sec-

tion are non-Gaussian, suggesting that the predicted behavior may extend beyond the assumptions of the theorem. Theorem C.3 is restricted to the scalar-output case K = 1. Nevertheless, the interpolation threshold in the proof arises from the spectral properties of a matrix of size p × N K in the general case, suggesting that the relevant transition should more generally occur when the number of parameters matches the number of scalar training constraints N K. We therefore conjecture that, under the assumptions of Theorem C.3, the lower bound on the expected risk exhibits a double descent peak at p = N K also for K > 1. Our numerical experiments provide evidence supporting this conjecture. Finally, in Appendix E we provide further results concerning upper bounds to the expected risk. Specifically, in Theorem E.1 we give the upper bound L(ϑ̂S ) ≤ L′ (S) +

2 R λmax (Cf (ϑ̂S )) 1 α ξmax . (10) N +1 λmin (ĤS (ϑ̂S ))

Similarly as in Eq. (9), this bound holds up to second-order corrections O(N −2 ), and is dominated by spectral properties of the Hessian matrix ĤS (ϑ̂S ) and the covariance matrix 2 Cf (ϑ̂S ). Here ξmax is defined analogously as before, and R is the rank of the Hessian matrix, i.e., the number of non-zero eigenvalues. Under a slightly different set of assumptions, we show in Theorem E.2 that this upper bound also peaks at interpolation p = N , with high probability in the same limit as above.

EMPIRICAL EVIDENCE OF DOUBLE DESCENT IN PQCS

The theoretical results presented in the previous section hold in the asymptotic limit and rely on assumptions about the local geometry of the loss landscape at local minima as well as the data being Gaussian. In this section, we provide numerical evidence that gradient-based data re-uploading PQCs display the predicted double descent behavior. In particular, we show that, when training converges, the test loss, used here as a proxy for the expected risk, exhibits a pronounced peak near the interpolation threshold p = N K, where p is the number of trainable parameters, N is the training set size, and K is the output dimension. We consider re-uploading PQCs where each input x is encoded through angle embeddings on n = 8 qubits, followed by trainable single-qubit rotations and nearest-neighbor CZ gates. The model output is obtained by measuring Pauli-X expectation values on K qubits, giving fθ (x) ∈ RK . Before evaluating the loss, we rescale these expectation-value outputs by a fixed factor c = 150. Empirically, we find that such a rescaling substantially improves optimization and makes convergence to local minima significantly more reliable. We perform numerical experiments on three datasets, shown in Fig. 1. The first is the MNIST-1D dataset [37], used as a multiclass classification task with K = 8 classes. The second is the Fashion MNIST dataset [38], also restricted to

4 (a)

MNIST-1D

(b)

(c)

Regression

8 Test loss

Fashion MNIST

6 4

N=21 N=30 N=39 N=48

2 200 400 Number of parameters

600

20 200 400 Number of parameters

600

200 400 Number of parameters

600

Figure 1. Test loss as a function of the number of parameters p for (a) MNIST-1D and (b) Fashion MNIST classification, and (c) a multidimensional synthetic regression task, for different training set sizes N . Vertical dashed lines indicate the predicted interpolation thresholds p = N K. In all three tasks, the test loss peaks close to interpolation and decreases again in the overparameterized regime, displaying a characteristic double descent profile. The shaded areas correspond to the standard deviation for ten independent experiment repetitions, each using independently sampled training data.

K = 8 classes. Finally, the third dataset is a synthetic multidimensional regression task with input dimension d = 8 and output dimension K = 8, generated from a noisy linear target function. For each dataset, we vary the number of training samples N ∈ {21, 30, 39, 48} and vary the circuit depth, thereby scanning across the underparameterized and overparameterized regimes. Further implementation details are provided in Appendix G. The vertical dashed lines in Fig. 1 mark the corresponding interpolation thresholds p = N K. These results provide empirical evidence that double descent appears in trainable PQCs optimized by gradient descent. Importantly, the second descent occurs in regimes where the number of trainable parameters exceeds the number of training samples, showing that increasing the size of the PQC further does not necessarily degrade performance on unseen data, contrary to the traditional expectation that larger models generalize worse. Instead, once the model passes the interpolation threshold and training succeeds, additional parameters can improve generalization. Note, however, that this does not imply that overparameterized models outperform underparameterized ones. Rather, it supports the main conclusion of our analysis: overparameterization can play a constructive role in gradient-based quantum machine learning. We provide further details in Appendix F and some perspectives in Appendix H. In particular, we show numerical results systematically testing the formal assumptions we introduce for Theorem C.3.

DISCUSSION

In this work, we investigate the scalability of gradientbased parameterized quantum circuits (PQCs) by analyzing how their expected risk evolves across different parameterization regimes. By combining an add-one-in perturbation analysis with spectral arguments from random matrix theory, and building on an existing lower bound on the expected risk, we show that the bound attains a maximum at the interpolation

threshold. We further derive a corresponding upper bound and show that it exhibits the same qualitative dependence on model size. Our analysis therefore establishes that, under mild assumptions, the risk bounds of a broad class of PQCs exhibit double descent behavior. We confirm these findings through comprehensive numerical experiments with data re-uploading PQCs on both classification and regression tasks. Our work offers a complementary perspective to the growing literature on scalable QML under the slogan “train classical, deploy quantum” [39–41] and extends previous findings of double descent in quantum kernel methods [29]. The presented results indicate that overparameterized regimes may be compatible with favorable learning behavior [25, 42, 43], providing grounds for cautious optimism about the prospects of deep PQCs. The reason why we advocate for cautious optimism is twofold. First, our analysis is restricted to trainable PQCs. Ref. [31] already pointed out that, in concentration-dominated regimes where trainability breaks down, larger models may not exhibit favorable generalization behavior. Our results should therefore be understood as a statement about what may be possible if trainability can be maintained as the number of parameters increases. We regard trainability as a prerequisite for practically relevant quantum machine learning and therefore focus on the generalization properties of models operating in such a regime. In this sense, we do not solve wellknown trainability problems of QML [12, 15, 44, 45]. At the same time, we expect the assumptions underlying our analysis to hold for a broad class of trainable PQCs. Consequently, the qualitative phenomenon identified here may remain relevant for future trainable quantum models of practical relevance. The second reason for caution is that our results do not imply that overparameterized PQCs outperform underparameterized ones. Rather, they show that increasing the number of parameters beyond the interpolation threshold need not be detrimental to generalization, contrary to the traditional expectation from statistical learning theory that increasing model size ultimately leads to poorer generalization [46–49]. Determin-

5 ing when overparameterization leads to overall improvements in expected risk, and under which conditions such improvements can be realized in practice, remain important research directions. Another noteworthy observation from our numerical experiments is that rescaling the PQC outputs by a large constant improves convergence to local minima. One possible explanation is that expectation values are typically concentrated around zero, so that rescaling increases the density of lowtraining-loss solutions in parameter space by shifting the target labels away from the boundaries of the output interval. Understanding this effect theoretically is likely of independent interest and is left for future work. Overall, this work provides a framework for understanding how the generalization performance of parameterized quantum circuits evolves under scaling. We expect these results to inform future studies of learning, generalization, and overparameterization in increasingly large quantum models, and ultimately provide cautious optimism for deep parameterized quantum circuits. Code and data availability. The code and data used in this study are available in GitHub [50]. Acknowledgments. The authors would like to thank Evan Peters for insightful comments on an earlier version of this manuscript. E.G.-F., C.B.-P. and J.E. acknowledge the BMFTR (MUNIQC-Atoms, HYBRID, QuSol, Hybrid++), the DFG (CRC 183 and SPP 2514), Bifold, the Quantum Flagship (Millenion, PasQuanS2), the Munich Quantum Valley, Berlin Quantum and the European Research Council (DebuQC) for financial support. E.G.-F. further acknowledges support from a 2023 Google PhD Fellowship. V.D. acknowledges support from the European Union’s Horizon Europe program through the ERC CoG BeMAIQuantum (Grant No. 101124342). V.D. and E.vN acknowledge the Dutch National Growth Fund (NGF) as part of the Quantum Delta NL programme. GPT-5.6 was used to assist with code development and mathematical consistency checks during the preparation of this manuscript. All scientific ideas, analyses, and conclusions are those of the authors, and all AI-assisted code and text were carefully reviewed, manually checked, and validated prior to inclusion. Disclaimer. The results, opinions, and conclusions expressed in this publication are not necessarily those of Volkswagen Aktiengesellschaft or those of the European Union or the European Research Council Executive Agency. None of the parties involved can be held responsible for them.

∗

[email protected] [1] V. Havlı́ček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019). [2] R. Sweke, J.-P. Seifert, D. Hangleiter, and J. Eisert, On the quantum versus classical learnability of discrete distributions, Quantum 5, 417 (2021).

[3] N. Pirnay, R. Sweke, J. Eisert, and J.-P. Seifert, A superpolynomial quantum-classical separation for density modelling, Physical Review A 107, 042416 (2023). [4] R. Molteni, S. C. Marshall, and V. Dunjko, Quantum machine learning advantages beyond hardness of evaluation, arXiv:2504.15964 (2025). [5] R. Molteni, C. Gyurik, and V. Dunjko, Exponential quantum advantages in learning quantum observables from classical data, npj Quantum Information 12, 19 (2026). [6] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Reviews of Modern Physics 91, 045002 (2019). [7] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017). [8] S. Masot-Llima, E. Gil-Fuster, C. Bravo-Prieto, J. Eisert, and T. Guaita, Prospects for quantum advantage in machine learning from the representability of functions, arXiv:2512.15661 (2025). [9] A. Barthe, M. Y. Rad, M. Grossi, and V. Dunjko, Quantum advantage in learning quantum dynamics via fourier coefficient extraction, arXiv:2506.17089 (2025). [10] R. Bandyopadhyay, R. Molteni, J. Eisert, V. Dunjko, and S. Jerbi, Provable learning separation for predicting timeevolution of quantum many-body systems, arXiv:2607.06472 (2026). [11] J. Eisert and J. Preskill, Mind the gaps: The fraught road to quantum advantage, arXiv:2510.19928 (2025). [12] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018). [13] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 1791 (2021). [14] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise-induced barren plateaus in variational quantum algorithms, Nature Communications 12, 6961 (2021). [15] E. R. Anschuetz and B. T. Kiani, Quantum variational algorithms are swamped with traps, Nature Communications 13, 7760 (2022). [16] A. A. Mele, A. Angrisani, S. Ghosh, S. Khatri, J. Eisert, D. Stilck França, and Y. Quek, Noise-induced shallow circuits and the absence of barren plateaus, Nature Physics 22, 751 (2026). [17] M. C. Caro, E. Gil-Fuster, J. J. Meyer, J. Eisert, and R. Sweke, Encoding-dependent generalization bounds for parametrized quantum circuits, Quantum 5, 582 (2021). [18] A. Abbas, D. Sutter, C. Zoufal, A. Lucchi, A. Figalli, and S. Woerner, The power of quantum neural networks, Nature Computational Science , 403 (2021). [19] L. Banchi, J. Pereira, and S. Pirandola, Generalization in quantum machine learning: A quantum information standpoint, PRX Quantum 2, 040321 (2021). [20] M. C. Caro, H.-Y. Huang, M. Cerezo, K. Sharma, A. Sornborger, L. Cincio, and P. J. Coles, Generalization in quantum machine learning from few training data, Nature Communications 13, 4919 (2022). [21] Y. Du, Z. Tu, X. Yuan, and D. Tao, Efficient measure for the expressivity of variational quantum algorithms, Physical Review Letters 128, 080506 (2022). [22] E. Gil-Fuster, J. Eisert, and C. Bravo-Prieto, Understanding quantum machine learning also requires rethinking generalization, Nature Communications 15, 2277 (2024).

6 [23] P. Rodriguez-Grasa, M. C. Caro, J. Eisert, E. Gil-Fuster, F. J. Schreiber, and C. Bravo-Prieto, A PAC-Bayesian approach to generalization for quantum models, arXiv:2603.22964 (2026). [24] E. Peters and M. Schuld, Generalization despite overfitting in quantum machine learning models, Quantum 7, 1210 (2023). [25] M. Belkin, D. Hsu, S. Ma, and S. Mandal, Reconciling modern machine-learning practice and the classical bias-variance trade-off, Proceedings of the National Academy of Sciences 116, 15849 (2019). [26] P. Nakkiran, G. Kaplun, Y. Bansal, T. Yang, B. Barak, and I. Sutskever, Deep double descent: where bigger models and more data hurt, Journal of Statistical Mechanics: Theory and Experiment 2021, 124003 (2021). [27] T. Hastie, A. Montanari, S. Rosset, and R. J. Tibshirani, Surprises in high-dimensional ridgeless least squares interpolation, Annals of Statistics 50, 949 (2022). [28] S. Mei and A. Montanari, The generalization error of random features regression: Precise asymptotics and double descent curve, Communications on Pure and Applied Mathematics 75, 667 (2022). [29] M. Kempkes, A. Ijaz, E. Gil-Fuster, C. Bravo-Prieto, J. Spiegelberg, E. van Nieuwenburg, and V. Dunjko, Double descent in quantum kernel methods, PRX Quantum 7, 010312 (2026). [30] D. Pranjić, M. Roth, and C. Tutschku, Grokking and epoch-wise double descent in quantum neural networks, arXiv:2607.08350 (2026). [31] Y. Du, Y. Yang, D. Tao, and M.-H. Hsieh, Problem-dependent power of quantum neural networks on multiclass classification, Physical Review Letters 131, 140601 (2023). [32] A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, and J. I. Latorre, Data re-uploading for a universal quantum classifier, Quantum 4, 226 (2020). [33] S. P. Singh, A. Lucchi, T. Hofmann, and B. Schölkopf, Phenomenology of double descent in finite-width neural networks, in International Conference on Learning Representations (ICLR 2022) (2022). [34] F. R. Hampel, R. E. M, R. P. J, and S. W. A, Robust Statistics: The Approach Based on Influence Functions, Vol. 196 (1986). [35] M. E. A. Seddik, C. Louart, M. Tamaazousti, and R. Couillet, Random matrix theory proves that deep learning representations of GAN-data behave as Gaussian mixtures, in International Conference on Machine Learning (PMLR, 2020). [36] R. Couillet and Z. Liao, Random Matrix Methods for Machine Learning (Cambridge University Press, 2022). [37] S. Greydanus and D. Kobak, Scaling down deep learning with MNIST-1D, in Proceedings of the 41st International Conference on Machine Learning, ICML (2024). [38] H. Xiao, K. Rasul, and R. Vollgraf, Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms, arXiv:1708.07747 (2017). [39] A. Huang, W. Maxwell, V. Belis, E. Peters, J. Pye, S. Jahangiri, and J. Bowles, Spectral Born machines: classically trainable quantum generative models for discrete data, arXiv:2607.06675 (2026). [40] E. Recio-Armengol, S. Ahmed, and J. Bowles, Train on classical, deploy on quantum: scaling generative quantum machine learning to a thousand qubits, arXiv:2503.02934 (2026). [41] A. Kurkin, K. Shen, S. Pielawa, H. Wang, and V. Dunjko, Universality and kernel-adaptive training for classically trained, quantum-deployed generative models, arXiv:2510.08476

(2025). [42] S. D’Ascoli, M. Refinetti, G. Biroli, and F. Krzakala, Double trouble in double descent: Bias and variance(s) in the lazy regime, in Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 119, edited by H. D. III and A. Singh (PMLR, 2020) pp. 2280–2290. [43] Y. Dar, V. Muthukumar, and R. G. Baraniuk, A farewell to the bias-variance tradeoff? An overview of the theory of overparameterized machine learning, arXiv:2109.02355 (2021). [44] E. Anschuetz, A unified theory of quantum neural network loss landscapes, in International Conference on Learning Representations, Vol. 2025 (2025) pp. 97859–97918. [45] S. Thanasilp, S. Wang, M. Cerezo, and Z. Holmes, Exponential concentration in quantum kernel methods, Nature Communications 15 (2024). [46] V. N. Vapnik and A. Y. Chervonenkis, On the uniform convergence of relative frequencies of events to their probabilities, in Measures of complexity: festschrift for alexey chervonenkis (Springer, 2015) pp. 11–30. [47] V. Vapnik, An overview of statistical learning theory, IEEE Transactions on Neural Networks 10, 988 (1999). [48] L. G. Valiant, A theory of the learnable, Commun. ACM 27, 1134–1142 (1984). [49] S. Shalev-Shwartz and S. Ben-David, Understanding machine learning (Cambridge University Press, Cambridge, England, 2014). [50] Cautious optimism for deep parameterized quantum circuits, GitHub repository (2026). [51] Y. Fang, K. Loparo, and X. Feng, Inequalities for the trace of matrix product, IEEE Transactions on Automatic Control 39, 2489 (1994). [52] P. Yaskov, Necessary and sufficient conditions for the Marchenko-Pastur theorem, Electronic Communications in Probability 21, 1 (2016). [53] C. Louart and R. Couillet, Concentration of measure and large random matrices with an application to sample covariance matrices, arXiv:1805.08295 (2021). [54] R. Mathias, Singular values and singular value inequalities (CRC Press, 2013). [55] K. Karhadkar, M. Murray, H. Tseran, and G. Montúfar, Mildly overparameterized ReLU networks have a favorable loss landscape, arXiv:2305.19510 (2024). [56] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, and N. Killoran, Pennylane: Automatic differentiation of hybrid quantumclassical computations, arXiv:1811.04968 (2018). [57] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Scikit-learn: Machine learning in python, The Journal of Machine Learning Research 12, 2825 (2011). [58] N. E. Karoui, The spectrum of kernel random matrices, Annals of Statistics 38, 1 (2010). [59] X. Cheng and A. Singer, The spectrum of random inner-product kernel matrices, Random Matrices: Theory and Applications 2, 1350010 (2013). [60] Z. Fan and A. Montanari, The spectral norm of random innerproduct kernel matrices, Probability Theory and Related Fields 173, 27 (2019).

Supplementary Material for “Cautious optimism for deep parameterized quantum circuits” Appendix A: Framework

Throughout this work, we consider supervised learning tasks on a data space Z := X × Y, where X ⊆ Rd denotes the input domain and Y ⊆ RK the output domain. Samples z = (x, y) are assumed to be drawn i.i.d. from an unknown but fixed distribution D on Z. Given a training set S = {zi }N i=1 , the goal is to learn a hypothesis fϑ ∈ F, where F = {fϑ : X → Y | ϑ ∈ Θ} denotes the hypothesis class and Θ ⊆ Rp the parameter space. The hypothesis class is realized by parameterized quantum circuits (PQCs). For a classical input x, the model prepares a quantum state ρ(x, ϑ) = U (x, ϑ)ρ0 U † (x, ϑ),

(A1)

where ρ0 is the initial state and U (x, ϑ) is the unitary implemented by the PQC. A common architecture consists of alternating data-dependent and parameter-dependent layers, U (x, ϑ) =

L Y

Vl (ϑ)Ul (x),

(A2)

l=1

commonly referred to as data re-uploading [32]. We consider vector-valued models fϑ : X → Y ⊆ RK with outputs fϑ (x) = [fϑ (x)]1 , . . . [fϑ (x)]K ,

(A3)

where each component is obtained as the expectation value [fϑ (x)]k = Tr{ρ(x, ϑ)Ok } ,

k = 1, . . . , K,

(A4)

of a Hermitian observable. For classification tasks with K > 1, the model output fϑ (x) is interpreted as a vector of class scores, while the predicted label is typically given by argmaxk [fϑ (x)]k . For K = 1, a setting commonly used in binary classification, the predicted label can instead be obtained by thresholding the output, e.g., via sign(fϑ (x)). Let a learning task be specified by a data distribution D(Z), which we assume to be fixed and unknown. We denote by N S = {(xi , yi )}N i=1 a training set of size N composed of i.i.d. samples zi := (xi , yi ) from the problem distribution S ∼ D . Next, let ℓ : Z × Θ → R≥0 be a twice-differentiable loss function, which given z ∈ Z and a parameter specification ϑ ∈ Θ acts as a measure of distance between the prediction fϑ (x) and the correct label y. From the loss function ℓ and the problem distribution D we define the expected risk functional L : Θ → R≥0 as L(ϑ) := E [ℓ(z, ϑ)] . z∼D

(A5)

The goal in supervised learning is to find a set of parameters that achieve a small expected risk, though this task cannot be tackled head-on due to the assumption that D is unknown. Rather, the common approach is to use a proxy quantity relative to a given training set S, which we denote empirical risk L̂S : Θ → R≥0 , defined as N

L̂S (ϑ) :=

1 X ℓ(zi , ϑ). N i=1

(A6)

This way, the empirical risk is an unbiased estimator for the expected risk. In the following, we abuse language in that we identify a function fϑ (also called hypothesis) with its parameters ϑ. Definition A.1 (Local minimum). Given the model F, a training set S, and the empirical risk functional L̂S , we call a parameter vector a local minimum, denoted by ϑ̂S , if it fulfills ∇ϑ L̂S (ϑ̂S ) = 0,

(A7)

ĤS (ϑ̂S ) ⪰ 0.

(A8)

To avoid confusion with other differential operators that appear below, ∇ϑ represents the gradient operator ⊺ ∂ ∂ · · · ∂ϑ ∇ϑ := ∂ϑ 1 p

(A9)

with respect to the parameters. Specifically, for any ϑ ∈ Θ, the gradient of the loss with respect to the parameters is a pdimensional column vector ∇ϑ L̂S (ϑ) ∈ Rp .

8 Remark (Local minima as output). From here on, we make the explicit assumption that standard gradient-based training techniques produce local minima as output. We denote by ĤS (ϑ) the Hessian matrix of the empirical risk ĤS (ϑ) = ∇ϑ ∇⊺ϑ L̂S (ϑ) =

∂ 2 L̂S (ϑ) ∂ϑi ∂ϑj

!p .

(A10)

i,j=1

The inverse of the Hessian. Our results below are stated in terms of the inverse of the Hessian matrix, which does not exist unless the matrix is full rank. Throughout this work, we therefore interpret [ĤS (ϑ)]−1 as the Moore-Penrose pseudoinverse of ĤS (ϑ), obtained by inverting only the non-zero eigenvalues. Appendix B: Expected risk at local minima

Since double descent is fundamentally a phenomenon of the expected risk, we first derive an explicit expression for the expected risk attained attained at local minima. To this end, we analyze the effect of adding a single sample to the training set. Following Ref. [33], we consider an add-one-in perturbation of the empirical risk. Under this perturbation, the empirical risk changes only slightly, and the corresponding local minimum is therefore expected to remain close to that of the original training set. This motivates a local expansion around the original local minimum. Definition B.1 (Add-one-in perturbation). Given a training set S and a datum z ′ ∈ Z not already-contained in S, we refer to S ′ = S ∪ {z ′ } as the add-one-in perturbation of S. We do not include the z ′ dependence of S ′ explicitly because below we only make statements that hold in expectation over z ∼ D. The distinction we keep in mind regarding notation is the size of the sets: for |S| = N we have |S ′ | = N + 1. Note that if we call DS the uniform distribution over S, and similarly for S ′ , then DS ′ = (1 − ε)DS + εδz′ . Here δz′ is a delta distribution centered at z ′ , and ′

ε=

1 . N +1

(B1)

In the regime of large N , we have that ε becomes small, which justifies calling DS ′ a “perturbation” of DS . Again, we abuse language and we identify the sets S and S ′ with their uniform distributions. Lemma B.2 (Influence function for add-one-in perturbation, adapted from Ref. [33, Prop. 2]). Let F be a hypothesis family, S a training set of size N . Let S ′ be the add-one-in perturbation of S, and let ϑ̂S , ϑ̂S ′ be local minima with respect to S and S ′ respectively. Then, the first-order contribution from adding samples to the training set is as h i−1 ℓ(z ′ , ϑ̂S ′ ) − ℓ(z ′ , ϑ̂S ) 1 = −∇ϑ ℓ(z ′ , ϑ̂S )⊺ ĤS (ϑ̂S ) . (B2) ∇ϑ ℓ(z ′ , ϑ̂S ) + O ε N Definition B.3 (Add-one-in loss). Suppose we are given the model F, a training set S and its corresponding local minimum ϑ̂S . Further, for any add-one-in perturbation S ′ = S ∪ {z ′ }, let ϑ̂S ′ be its local minimum. Then, we denote the add-one-in loss L′ (S) as h i L′ (S) := ′E ℓ(z ′ , ϑ̂S ′ ) . (B3) z ∼D

Definition B.4 (Uncentered covariance of loss gradients). For any parameter vector ϑ ∈ Rp , we define the uncentered covariance matrix of the loss gradients as C(ϑ) := E [∇ϑ ℓ(z, ϑ)∇ϑ ℓ(z, ϑ)⊺ ] . z∼D

(B4)

Lemma B.5 (Error decomposition, Theorem 3 in Ref. [33]). Let F be the hypothesis family of a learning model, let S be a training set of size N , let ℓ be a twice-differentiable loss function, let ϑ̂S be a local minimum, let L′ (S) be the add-one-in loss, and let ĤS (ϑ̂S ) be the Hessian matrix of the empirical loss at the local minimum. Then, the expected risk L(ϑ̂S ) fulfills h i−1 1 1 ′ Tr ĤS (ϑ̂S ) C(ϑ̂S ) + O . (B5) L(ϑ̂S ) = L (S) + N +1 N2 We refer to the summand with the trace as the complexity term.

9 Proof. We prove this statement directly, starting from Lemma B.2 and taking an expectation over D. We start by taking the equation in the Lemma and substituting ε=

1 , N +1

(B6)

in order to get h i−1 ℓ(z ′ , ϑ̂S ′ ) − ℓ(z ′ , ϑ̂S ) 1 = −∇ϑ ℓ(z ′ , ϑ̂S )⊺ ĤS (ϑ̂S ) ∇ϑ ℓ(z ′ , ϑ̂S ) + O ε N h i−1 1 ℓ(z ′ , ϑ̂S ′ ) − ℓ(z ′ , ϑ̂S ) = −ε∇ϑ ℓ(z ′ , ϑ̂S )⊺ ĤS (ϑ̂S ) ∇ϑ ℓ(z ′ , ϑ̂S ) + O N2 h i −1 1 1 ′ ⊺ ′ =− . ∇ϑ ℓ(z , ϑ̂S ) ĤS (ϑ̂S ) ∇ϑ ℓ(z , ϑ̂S ) + O N +1 N2 Next we perform a common step to turn a vector-matrix-vector product into the trace of a matrix-matrix product as h i−1 1 1 ′ ′ ′ ⊺ ′ ′ ℓ(z , ϑ̂S ) − ℓ(z , ϑ̂S ) = − ∇ϑ ℓ(z , ϑ̂S ) + O ∇ϑ ℓ(z , ϑ̂S ) ĤS (ϑ̂S ) N +1 N2 i h −1 1 1 ′ ⊺ ′ =− Tr ∇ϑ ℓ(z , ϑ̂S ) ĤS (ϑ̂S ) ∇ϑ ℓ(z , ϑ̂S ) + O N +1 N2 h i −1 1 1 Tr ĤS (ϑ̂S ) ∇ϑ ℓ(z ′ , ϑ̂S )∇ϑ ℓ(z ′ , ϑ̂S )⊺ + O . =− N +1 N2

(B7) (B8) (B9)

(B10) (B11) (B12)

We now take the expectation over z ′ ∼ D on both sides and simplify the terms that do not depend on z ′ , in order to get h i−1 h i 1 1 ′ ′ ′ ′ ⊺ ′ Ĥ ( ϑ̂ ) E ℓ(z , ϑ̂ ) − ℓ(z , ϑ̂ ) = E − Tr ∇ ℓ(z , ϑ̂ )∇ ℓ(z , ϑ̂ ) + O , (B13) S S S S ϑ S ϑ S z ′ ∼D z ′ ∼D N +1 N2 h h i i−1 h i h i 1 1 ′ ′ ⊺ ′ ′ ′ E ∇ ℓ(z , ϑ̂ )∇ ℓ(z , ϑ̂ ) + O Ĥ ( ϑ̂ ) Tr . (B14) E ℓ(z , ϑ̂ ) − E ℓ(z , ϑ̂ ) = − ϑ S ϑ S S S S S z ′ ∼D z ′ ∼D z ′ ∼D N +1 N2 Notice that the two expectations on the left-hand side are precisely the add-one-in loss and the expected risk, and the expectation value on the right-hand side corresponds to the uncentered covariance of the gradients. From here we need only shift terms around to complete the proof, h i−1 1 1 ′ C(ϑ̂S ) + O L (S) − L(ϑ̂S ) = − Tr ĤS (ϑ̂S ) (B15) N +1 N2 h i−1 1 1 C(ϑ̂S ) + O Tr ĤS (ϑ̂S ) . (B16) L(ϑ̂S ) = L′ (S) + N +1 N2

We next reproduce the lower bound for this error decomposition, following the formalism in Ref. [33]. To this end, we restrict our attention to the scalar-output case K = 1 and introduce a mild regularity assumption on the loss gradients. We start by introducing two auxiliary lemmas that will be used in the proof. Lemma B.6. Let A, B ∈ Rm×m be symmetric matrices, with A, B ≥ 0 further being positive semi-definite. Then λmin (A) Tr{B} ≤ Tr{AB} ≤ λmax (A) Tr{B}.

(B17)

Here, λmin (A) corresponds to the smallest eigenvalue of A, and similarly for λmax . The proof of Lemma B.6 can be found in Ref. [51]. Lemma B.7. Under Assumption 1 and in the case K = 1, the complexity term from Lemma B.5 fulfills h h i−1 i−1 2 ⊺ Tr ĤS (ϑ̂S ) C(ϑ̂S ) = α E ξz Tr ĤS (ϑ̂S ) ∇ϑ fϑ̂S (z)∇ϑ fϑ̂S (z) . z∼D ξz2 >0

(B18)

10 Proof. We prove this statement directly. We first use the chain rule to relate ∇ϑ ℓ(z, ϑ) and ∇ϑ fϑ (z), and then we separate the expectation value into the two disjoint conditions ξz2 > 0 and ξz2 = 0, !p p ∂ℓ(z, ϑ) ∂fϑ (z) ∂ℓ(z, ϑ̂s ) ∂ℓ(z, ϑ) = ∇ϑ fϑ (x). (B19) ∇ϑ ℓ(z, ϑ) := = ∂ϑm ∂fϑ (x) ∂ϑm m=1 ∂fϑ (x) m=1

p In the case K = 1, ∂ℓ(z,ϑ) ∂fϑ (x) ∈ R is a scalar and ∇ϑ fϑ (x) ∈ R is a p-dimensional column vector. We now substitute the chain

rule into the definition of C(ϑ̂S ), to obtain   2 i h ′ ∂ℓ(z , ϑ̂ ) S ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺  . C(ϑ̂S ) = ′E ∇ϑ ℓ(z ′ , ϑ̂S )∇ϑ ℓ(z ′ , ϑ̂S )⊺ = ′E  z ∼D z ∼D ∂fϑ (x)

(B20)

The squared term inside the expectation value is exactly ∥∇f ℓ(z ′ , ϑ̂S )∥2 in the case of K = 1, which we defined as ξz2′ in Assumption 1,   2 h i ′  ∂ℓ(z , ϑ̂S ) ∇ϑ f (x′ )∇ϑ f (x′ )⊺  = E ξz2′ ∇ϑ f (x′ )∇ϑ f (x′ )⊺ . (B21) E ϑ̂S ϑ̂S ϑ̂S ϑ̂S z ′ ∼D z ′ ∼D ∂fϑ (x) Now, we split the expectation value over two disjoint conditions Ez′ ∼D [ · ] = P(ξz2′ > 0) E z′ ∼D [ · ] + P(ξz2′ = 0) E z′ ∼D [ · ]. ξz2′ >0 2 Assumption 1 can be restated as: there exists an α > 0 such that Pz∼D (ξz > 0) = α. With this, we obtain

h E ′

z ∼D

i h i ξz2′ ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺ = P(ξz2′ > 0) ′E ξz2′ ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺

ξz2′ =0

(B22)

z ∼D ξz2′ >0

+ P(ξz2′ = 0) ′E

z ∼D ξz2′ =0

= α ′E

z ∼D ξz2′ >0

ξz2′ ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺

h i ξz2′ ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺ .

(B23) (B24)

The second summand vanishes due to the ξz2′ = 0 condition. Directly substituting these in the complexity term finishes the proof,    h h  i−1 h i−1 i α ′E ξz2′ ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺ ĤS (ϑ̂S ) (B25) Tr ĤS (ϑ̂S ) C(ϑ̂S ) = Tr   z ∼D   2 ξz′ >0 h i−1 = α ′E ξz2′ Tr ĤS (ϑ̂S ) ∇ϑ fϑ̂S (x′ )∇ϑ fϑ̂S (x′ )⊺ . (B26) z ∼D ξz2′ >0

Assumption 1 (Regularity assumption). There exists a sample z ∼ D with non-zero probability α such that ξz2 := ∥∇f ℓ(z, ϑ̂S )∥2 > 0, where ∇f denotes the gradient with respect to the model output. For the mean-squared error loss considered in this work, Assumption 1 merely requires that predictions of the model with parameters ϑ̂S differ from the true targets on a subset of the data distribution with non-zero probability. Definition B.8 (Uncentered covariance of function gradients). For any parameter vector ϑ ∈ Rp , we define the uncentered covariance matrix of the function gradients as Cf (ϑ) := E [∇ϑ fϑ (x)∇ϑ fϑ (x)⊺ ] . z∼D ξz2 >0

(B27)

Combining the error decomposition from Lemma B.5 with Assumption 1 yields the following lower bound on the expected risk in terms of the spectrum of Cf (ϑ̂S ) and the Hessian matrix of the empirical risk ĤS (ϑ̂S ) at local minima.

11 Theorem B.9 (Lower bound on the expected risk, Theorem 4 Ref. [33]). Under the setting of Lemma B.5, under Assumption 1 and for K = 1, the expected risk L(ϑ̂S ) at the local minimum ϑ̂S fulfills the lower bound L(ϑ̂S ) ≥ L′ (S) +

2 λmin (Cf (ϑ̂S )) 1 α ξmin . N + 1 λmin (ĤS (ϑ̂S ))

(B28)

2 Here, we define ξmin = min z∼D {ξz2 }, and λmin (·) the smallest non-zero eigenvalue of its matrix argument. ξz2 >0

Proof. The proof follows directly from Lemmas B.5, B.6, and B.7. 2 We first provide a lower bound for the complexity term in Lemma B.7 using the definition of ξmin h i−1 h i−1 2 ⊺ Tr ĤS (ϑ̂S ) C(ϑ̂S ) = α E ξz Tr ĤS (ϑ̂S ) ∇ϑ fϑ̂S (z)∇ϑ fϑ̂S (z) z∼D ξz2 >0

 o i−1 h 2 ≥ α ξmin Tr Ĥ (ϑ̂ ) E ∇ϑ fϑ̂S (z)∇ϑ fϑ̂S (z)⊺ z∼D   S S ξz2 >0 h i−1 2 = α ξmin Tr ĤS (ϑ̂S ) Cf (ϑ̂S ) .  h

Lemma B.6 allows us to lower-bound this quantity in terms of the eigenvalues of Cf (ϑ̂S ) as h h i−1 i−1 Tr ĤS (ϑ̂S ) Cf (ϑ̂S ) ≥ λmin (Cf (ϑ̂S )) Tr ĤS (ϑ̂S ) .

(B29)

(B30)

(B31)

(B32)

Now, under the assumption that ĤS (ϑ̂S ) is full-rank, we can immediately relate the trace of its inverse to the eigenvalues of the Hessian. Namely: the trace of a matrix is the sum of its eigenvalues, and for Hermitian matrices, the eigenvalues of the inverse are the inverse of the eigenvalues λi ([ĤS (ϑ̂S )]−1 ) = (λi (ĤS (ϑ̂S )))−1 , h h p i−1 X i−1 Tr ĤS (ϑ̂S ) = λi ĤS (ϑ̂S ) (B33) i=1

R X i=1 λi

≥

ĤS (ϑ̂S )

1 λmin (ĤS (ϑ̂S ))

(B34)

(B35)

We call R = rank(ĤS (ϑ̂S )), namely the number of non-zero eigenvalues of the Hessian. From the first to the second line, we replace a sum over p eigenvalues (dimension of the matrix) with a sum over the R-many inverses of non-zero eigenvalues. In the last step, we lower-bounded the value of a sum of positive numbers by the value of the largest element in the sum. We can finally insert these bounds in the error decomposition (Eq. (B5) from Lemma B.5), to get h i−1 1 1 L(ϑ̂S ) = L′ (S) + Tr ĤS (ϑ̂S ) C(ϑ̂S ) + O (B36) N +1 N2 h i−1 1 1 2 α ξmin Tr ĤS (ϑ̂S ) Cf (ϑ̂S ) +O (B37) ≥ L′ (S) + N +1 N2 h i−1 2 α ξmin 1 λmin (Cf (ϑ̂S )) Tr ĤS (ϑ̂S ) +O (B38) ≥ L′ (S) + N +1 N2 ! 2 α ξmin λmin (Cf (ϑ̂S )) 1 1 ′ ≥ L (S) + +O (B39) N +1 N2 λmin (ĤS (ϑ̂S )) 2 λmin (Cf (ϑ̂S )) 1 α ξmin 1 . (B40) = L′ (S) + +O N + 1 λmin (ĤS (ϑ̂S )) N2 This completes the proof. In particular, the lower bound is governed by the smallest non-zero eigenvalue of the Hessian matrix of the empirical risk at local minima, λmin (ĤS (ϑ̂S )), a quantity that will play a central role in our subsequent analysis. A complementary upper bound, obtained using analogous techniques, is presented in Theorem E.1 in Appendix E.

12 Appendix C: Double descent behavior in PQCs

We now establish that, as the number of trainable parameters p is varied, the lower bound on the expected risk exhibits a double descent peak. While this does not directly imply that the expected risk itself exhibits double descent, numerical results presented later provide evidence for a corresponding phenomenon. Restricting to the scalar-output case K = 1, we show that the peak occurs at p = N , corresponding to the interpolation threshold. Numerical results presented later indicate the same qualitative behavior for K > 1. We therefore keep the general notation in terms of K throughout, although the rigorous results derived below are restricted to K = 1. Our analysis is inspired by the framework of Ref. [33], with suitable modifications to the underlying assumptions. A central ingredient in the derivation is a decomposition of the Hessian of the empirical risk into the so-called outer-product Hessian ĤSo (ϑ) and the functional Hessian ĤSf (ϑ), ĤS (ϑ) = ĤSo (ϑ) + ĤSf (ϑ),

(C1)

where ĤSo (ϑ) =

N h i 1 X ∇ϑ fϑ (xi ) ∇f ∇⊺f ℓ(zi , ϑ) ∇ϑ fϑ (xi )⊺ , N i=1

ĤSf (ϑ) =

1 X ⊺ ∇ ℓ(zi , ϑ)∇ϑ ∇⊺ϑ fϑ (xi ). N i=1 f

(C2)

(C3)

To prove our main result, Theorem C.3, we require several assumptions on the quantities appearing in the Hessian decomposition and the ensuing spectral analysis. This section culminates with the statement of Theorem C.3, and its proof (including a derivation for the decomposition of the Hessian) is provided in Appendix D in a self-contained way. We first introduce two quantities that will play a central role in formulating these assumptions. p Definition C.1 (Sample Jacobian of the function). For a training set S = {(xi , yi )}N i=1 and parameter vector ϑ ∈ R , the sample Jacobian of the function is defined as ZS (ϑ) := ∇ϑ fϑ (x1 ), . . . , ∇ϑ fϑ (xN ) ∈ Rp×N K . (C4)

Definition C.2 (Uncentered sample covariance of function gradients). For a training set S = {(xi , yi )}N i=1 and parameter vector ϑ ∈ Rp , the uncentered sample covariance of the function gradients is defined as ĈSf (ϑ) :=

1 ZS (ϑ)ZS (ϑ)⊺ ∈ Rp×p . N

(C5)

The following assumptions serve different purposes in the analysis. Assumption 2 concerns the behavior of the functional Hessian near interpolation, while the subsequent assumptions concern properties of the outer-product Hessian and the covariance of function gradients. Assumption 2. In the underparameterized regime p < N K, the largest eigenvalue λmax ĤSf (ϑ̂S ) of the functional Hessian at local minima decreases with the number of trainable parameters p in a neighborhood of the interpolation threshold p = N K. For the mean-squared error loss, the functional Hessian depends explicitly on the training residuals. Assumption 2 is therefore consistent with the expectation that these residuals decrease as the number of trainable parameters increases. We investigate this assumption numerically later. The next assumption concerns the second term in the Hessian decomposition, namely the outer-product Hessian ĤSo (ϑ). For the scalar-output mean-squared error loss, this term reduces to the uncentered sample covariance of the function gradients, ĈSf (ϑ). This follows immediately from the fact that the first derivative of the loss is ∇f ℓ(zi , ϑ) = (fϑ (xi ) − yi ) while its second derivative ∇f (fϑ (xi ) − yi )) = 1 . We can, therefore, write ĤSo (ϑ) = ĈSf (ϑ) =

1 ZS (ϑ)ZS (ϑ)⊺ . N

(C6)

Our next assumption addresses the rank of the sample Jacobian of the function in the overparameterized regime. Assumption 3. In the overparameterized regime p ≥ N K, the sample Jacobian of the function at local minima ZS (ϑ̂S ) ∈ Rp×N K has full rank, rank ZS (ϑ̂S ) = N K. (C7)

13 Intuitively, Assumption 3 means that the trainable parameters affect the N K training outputs in sufficiently many independent ways. One may expect this condition to hold in sufficiently expressive PQCs, and we provide numerical evidence that it is satisfied in the settings considered here. We further assume that the spectrum of the uncentered (sample) covariance of function gradients remains stable during training. Assumption 4. The smallest eigenvalue of the uncentered sample covariance of the function gradients remains comparable to its value at initialization uniformly in the number of trainable parameters p: λmin ĈSf (ϑ̂S ) ≍ λmin ĈSf (ϑ0 ) . (C8) More precisely, there exist constants A, B > 0 independent of p such that A λmin ĈSf (ϑ0 ) ≤ λmin ĈSf (ϑ̂S ) ≤ B λmin ĈSf (ϑ0 ) .

(C9)

We further assume that an analogous relation holds for the uncentered covariance of function gradients Cf , with constants A′ , B ′ > 0. These relations hold with high probability with respect to the random draw of the training set S, the initialization ϑ0 , and the local minimum ϑ̂S returned by the training procedure. Assumption 4 formalizes a form of spectral stability during training. We provide empirical evidence that this assumption holds in the PQCs considered here. The next Assumption 5 is a regularity condition on the covariance of function gradients at initialization. It ensures that the spectrum remains bounded and non-degenerate as the number of trainable parameters varies, thereby enabling random matrix theoretic analysis, as we discuss in more detail in Appendix D. Assumption 5 (Spectrum of the uncentered covariance of function gradients at initialization). The uncentered covariance of function gradients at initialization has bounded spectrum. That is, there exist constants c, c′ > 0, independent of the number of trainable parameters p, such that c ≤ λmin (Cf (ϑ0 )) ≤ λmax (Cf (ϑ0 )) ≤ c′ .

(C10)

We expect this assumption to hold for circuit architectures with non-redundant parameterizations and generic initializations. Our next assumption relates the smallest non-zero eigenvalue of the empirical-risk Hessian to the spectra of its outer-product and functional contributions. Assumption 6 (Weyl-like inequality). The smallest non-zero eigenvalue of the empirical-risk Hessian at ϑ̂S satisfies λmin (ĤSo (ϑ̂S )) + λmin (ĤSf (ϑ̂S )) ≤ λmin (ĤS (ϑ̂S )) ≤ λmin (ĤSo (ϑ̂S )) + λmax (ĤSf (ϑ̂S )).

(C11)

Note that if rank(ĤSo (ϑ̂S )) = rank(ĤS (ϑ̂S )), then Weyl’s inequality directly implies Assumption 6. One may therefore interpret this assumption as requiring that the functional Hessian does not alter the rank structure induced by the outer-product Hessian. The final assumption concerns the add-one-in loss term appearing in the expected-risk decomposition. Assumption 7 (Add-one-in loss). The add-one-in loss L′ (S) does not dominate the lower bound in Theorem B.9 in a neighborhood of the interpolation threshold. Since L′ (S) corresponds to the loss obtained after perturbing the training set by a single sample, one may expect it to remain close to the empirical risk of the trained model, which is typically small around the interpolation threshold. We are now ready to present the main result of this work. The following theorem shows that, under the above assumptions, the lower bound on the expected risk from Theorem B.9 attains a maximum at the interpolation threshold, leading to a double descent peak as the number of parameters increases. Theorem C.3 (Double descent peak at interpolation). Let S = {(xi , yi )}N i=1 be a dataset with d-dimensional Gaussian normal inputs xi ∼ N (0, Id ) and outputs yi ∈ R. Let fϑ : Rd → R be a parameterized quantum circuit as in Eq. (A4) for K = 1, consisting of p parameters. Under Assumptions 1–7, and considering the mean-squared error loss function, the lower bound on L(ϑ̂S ) at a local minimum from Theorem B.9 attains a maximum at p = N with high probability in the limit N, p → ∞ with fixed ratio p/N → γ.

14 Proof sketch. Starting from the lower bound of Theorem B.9, it suffices to analyze the dependence on p of the Hessian-dependent contribution appearing in the bound. To this end, we use the decomposition from Eq. (C1). Using the upper bound from Assumption 6 together with Assumption 4, this term can be further bounded by the smallest non-zero eigenvalue of the uncentered sample covariance of function gradients at initialization. Consequently, we need to characterize the behavior of λmin (ĈSf (ϑ0 )) as a function of the parameter dimension p. Under Assumption 5 and using standard singular value inequalities, central results in random matrix theory dictate that λmin (ĈSf (ϑ0 )) approaches zero at the interpolation threshold p = N , while remaining bounded away from zero on either side of this point with high probability. It therefore remains to control the second term in the denominator involving λmax (ĤSf (ϑ̂S )). In the underparameterized regime, Assumption 2 implies that λmax (ĤSf (ϑ̂S )) decreases with increasing p, causing the lower bound to increase. In the overparameterized regime, we can show that Assumption 3 implies ĤSf (ϑ̂S ) = 0 at local minima. Combining these observations with the fact that λmin (ĈSf (ϑ0 )) is minimized at the interpolation threshold, we conclude that the lower bound attains its maximum at p = N . A corresponding result on the double descent peak holds for the upper bound derived in Appendix E; see Theorem E.2, which shows that the upper bound likewise has a maximum at the interpolation threshold p = N . Theorem C.3 is restricted to the scalar-output case K = 1. Nevertheless, the interpolation threshold in the proof arises from the spectral properties of the sample Jacobian ZS (ϑ) ∈ Rp×N K , suggesting that the relevant transition should more generally occur when the number of parameters matches the number of scalar training constraints N K. This motivates the following conjecture. Conjecture C.4 (Interpolation threshold for K > 1). Under the setting of Theorem C.3, except that K > 1, the lower bound on the expected risk exhibits a double descent peak at p = N K. Our numerical experiments provide evidence supporting this conjecture. Appendix D: Proof of Theorem C.3

In this section, we prove Theorem C.3, which establishes that the lower bound on the expected risk from Theorem B.9 attains its maximum at the interpolation threshold. We first collect several auxiliary results concerning the spectrum of the Hessian of the empirical risk and its decomposition into outer-product and functional contributions. These ingredients are then combined with random matrix arguments to characterize the behavior of the lower bound as a function of the parameter dimension. Lemma D.1 (Decomposition of the Hessian). For a given training set S = {zi := (xi , yi )}N i=1 , we consider the empirical risk PN L̂S (ϑ) = N1 i=1 ℓ(zi , ϑ). We define the Hessian of the empirical risk ĤS (ϑ) := ∇ϑ ∇⊺ϑ L̂S (ϑ). Then, the decomposition ĤS (ϑ) = ĤSo (ϑ) + ĤSf (ϑ), outer-product Hessian: ĤSo (ϑ) :=

(D1)

N X

h i 1 ∇ϑ fϑ (xi ) ∇f ∇⊺f ℓ(zi , ϑ) ∇⊺ϑ fϑ (xi ), N i=1

(D2)

functional Hessian: ĤSf (ϑ) :=

1 X ⊺ ∇ ℓ(zi , ϑ)∇ϑ ∇⊺ϑ fϑ (xi ) N i=1 f

(D3)

holds. Proof. The proof follows the direct application of the chain rule, being careful about the sizes of the algebraic objects involved. (⊺) To be precise, we treat ∇ϑ and ∇⊺ϑ as operator-valued column and row vectors, respectively; and similarly for ∇f . ⊺ ⊺ For any x, we have fϑ (x) ∈ RK , ∇ϑ fϑ (x) ∈ RK×p , and ∇ϑ ∇ϑ fϑ (x) ∈ RK×(p×p) . The second ∇ϑ acts along a different direction than the original dimension of the output. For any z = (x, y), when instantiating the chain rule, we treat ℓ(z, ϑ) as ℓ((f, y); ϑ), with f = f (x). Then, for any z, we write ∇⊺ϑ ℓ(z, ϑ) = ∇⊺f ℓ(z, ϑ) ∇⊺ϑ fϑ (x) . | {z } | {z } | {z } 1×p

1×K

(D4)

K×p

Then, for the second derivative, we have 1×K

∇ϑ |

{

∇⊺f ℓ(z, ϑ) {z

p×K

= ∇ϑ fϑ (x) ∇f ∇⊺f ℓ(z, ϑ) . {z } } | {z } | p×K

K×K

(D5)

15 We use this for the chain rule of the second derivative, to get ∇ϑ ∇⊺ϑ ℓ(z, ϑ) = ∇ϑ ∇⊺f ℓ(z, ϑ)∇⊺ϑ fϑ (x) | {z } p×p = ∇ϑ ∇⊺f ℓ(z, ϑ) ∇⊺ϑ fϑ (x) + ∇⊺f ℓ(z, ϑ) (∇ϑ ∇⊺ϑ fϑ (x)) = ∇ϑ fϑ (x) ∇f ∇⊺f ℓ(z, ϑ) ∇⊺ϑ fϑ (x) + ∇⊺f ℓ(z, ϑ) ∇ϑ ∇⊺ϑ fϑ (x) . | {z } | {z } {z } | {z } | {z } | p×K

K×K

K×p

1×K

(D6) (D7) (D8)

K×(p×p)

Having proven this for arbitrary z = (x, y), the linearity of the differential operators guarantees that it holds for each summand in ĤS (ϑ), and hence for the whole sum. For completeness, we also give expressions for the entries of the two summands, avoiding the matrix-matrix product notation. Written out in terms of scalars, the (j, j ′ )th entry of each of the summands in terms of the k th entry of the function [fϑ (x)]k takes the form h i ĤSo (ϑ)

1 X X ∂[fϑ (xi )]k ∂ 2 ℓ(zi , ϑ) ∂ [fϑ (xi )]k′ , = N i=1 ′ ∂ϑj ∂fk ∂fk′ ∂ϑj ′ j,j ′

(D9)

k,k =1

ĤSf (ϑ)

i j,j

= ′

1 X X ∂ℓ(zi , ϑ) ∂ 2 [fϑ (xi )]k . N i=1 ∂fk ∂ϑj ∂ϑj ′

(D10)

k=1

Lemma D.2 (Theorem 2.1 in Ref. [52]). Let X be a matrix with N columns xp ∈ Rp . Assume xp are isotropic (i.e., E[xp x⊺p ] = Ip ) random vectors and so-called concentration of quadratic forms " −1 −1 # 1 ⊺ 1 1 ⊺ ⊺ x XX + εIp XX + εIp xp − Tr →0 (D11) p p N N is satisfied, as N → ∞ with p = p(N ) and p/N → γ. Then, the empirical spectral distribution of the sample covariance 1 ⊺ N XX converges weakly to the Marčenko-Pastur (MP) distribution µMP in the asymptotic limit. Its non-zero eigenvalue density is given by p (λ+ − λ)(λ − λ− ) µMP (dλ) = dλ, (D12) 2πγλ √ λ± = (1 ± γ)2 , (D13) where λ denotes an eigenvalue of the sample covariance matrix N1 XX ⊺ , and λ− and λ+ denote the lower and upper edges of the support of the non-zero spectrum, respectively. Definition D.3 (Exponentially-concentrated random vectors [35, Def. 2.1]). Let x ∈ Rp a random vector x ∼ P. We say x is an exponentially-concentrated random vector if, for any 1-Lipschitz function h : Rp → R, the following holds: h i 2 P f (x) − E [f (x)] > t ≤ Ae−Bt , (D14) x∼P

x∼P

for constants A, B > 0. Lemma D.4 (Theorem 2.40 in Ref. [53]). Let x ∈ Rp be an exponentially-concentrated random vector x ∼ P, and let A ∈ Rp×p be a matrix with bounded operator norm. Then, x⊺ Ax − Tr{A} is sharply concentrated. Definition D.5 (Data whitening). Let us call gi = ∇ϑ fϑ0 (xi ) ∈ Rp for each i ∈ {1, . . . , N }, and recall that Cf (ϑ0 ) := Ezi ∼D, ξz2 >0 [gi gi⊺ ]. Under Assumption 5, we refer to the following as whitened random vectors: g̃i = Cf (ϑ0 )−1/2 gi . At the matrix level, we define Z̃S (ϑ0 ) = (g̃1 , . . . , g̃N ) ∈ Rp×N , which corresponds to Z̃S (ϑ0 ) = Cf (ϑ0 )−1/2 ZS (ϑ0 ). Anal˜ ogously, we define C̃Sf (ϑ0 ) = Z̃S (ϑ0 )Z̃S (ϑ0 )⊺ /N (we drop the hat ĈSf to avoid cluttering notation), which corresponds to C̃Sf (ϑ0 ) = Cf (ϑ0 )−1/2 ĈSf (ϑ0 )Cf (ϑ0 )−1/2 . We call C̃Sf (ϑ0 ) the whitened uncentered sample covariance of the function gradient at initialization.

16 Lemma D.6 (Gaussian vectors [35, Prop. 2.2]). Let x ∼ N (0, Id ). Then x is an exponentially-concentrated random vector. Lemma D.7 (Lipschitz stability of concentrated random vectors [35, Prop. 2.3]). Let x ∈ Rd be an exponentially-concentrated random vector, and let G : Rd → Rp be an L-Lipschitz function, where L may depend on p. Then the random vector G(x) is also exponentially concentrated. Lemma D.8 (Properties of whitened vectors). Under Assumption 5, let fϑ be a parameterized quantum circuit as in Eq. (A4) with K = 1, and let g̃i = Cf (ϑ0 )−1/2 gi with gi = ∇ϑ fϑ0 (xi ) ∈ Rp as in Def. D.5. For Gaussian normal inputs xi ∼ N (0, Id ), the vectors g̃i are independent, isotropic, and fulfill concentration of quadratic forms. Proof. Consider the notation from Def. D.5. Then, ZS (ϑ0 ) = (g1 , . . . , gN ) ∈ Rp×N , and still ĈSf (ϑ0 ) = ZS (ϑ0 )ZS (ϑ0 )⊺ /N . Since the training data xi are sampled i.i.d. and ϑ0 is fixed at initialization, it follows that the random vectors gi are independent among themselves. This implies that g̃i are also independent. We immediately have that the whitened random vectors are isotropic Ezi ∼D [g̃i g̃i⊺ ] = Ip , since i h E [g̃i g̃i⊺ ] = E Cf (ϑ0 )−1/2 gi gi⊺ Cf (ϑ0 )−1/2 = Cf (ϑ0 )−1/2 E [gi gi⊺ ] Cf (ϑ0 )−1/2 (D15) zi ∼D

zi ∼D

−1/2

= Cf (ϑ0 )

Cf (ϑ0 )Cf (ϑ0 )

−1/2

= Ip .

(D16)

⊺ Here, we have used that Cf (ϑ0 )−1/2 = Cf (ϑ0 )−1/2 is symmetric and does not depend on zi . Since further xi ∼ N (0, Id ), we know by Lemma D.6 that they are concentrated random vectors. For standard PQCs, the map x 7→ ∇ϑ fϑ0 (x) is Lipschitz continuous, since the parameter-shift rule expresses gradient components as differences of Lipschitz-continuous function evaluations. It therefore follows from Lemma D.7 that the function gradients gi = ∇ϑ fϑ0 (xi ) are also concentrated random vectors. Moreover, by construction, g̃i = Cf (ϑ0 )−1/2 gi . Since Cf (ϑ0 )−1/2 is a linear map with bounded operator norm by Assumption 5, a second application of Lemma D.7 shows that the whitened vectors g̃i are concentrated random vectors as well. Therefore, the concentration of quadratic forms follows from Lemma D.4. Lemma D.9 (Whitened-vector matrix converges to MP law). In the setting of Lemma D.8, the empirical spectral distribution of the whitened uncentered sample covariance of the function gradient at initialization C̃Sf (ϑ0 ) converges weakly almost surely to the Marčenko-Pastur law in the limit where N, p → ∞ and p/N → γ. Proof. Our strategy is to invoke Lemma D.2, which requires two conditions: independent, isotropic columns and concentration of quadratic forms, which we confirm using Lemma D.8. Corollary D.10 (Smallest eigenvalue convergence for whitened vectors). The smallest non-zero eigenvalue of pthe whitened uncentered sample covariance of the function gradient at initialization C̃Sf (ϑ0 ) converges almost surely to (1 − p/N )2 . Proof. Follows directly from Lemma D.9. Lemma D.11 (Smallest eigenvalue convergence for original vectors). In the setting of Lemma D.8, the smallest non-zero eigenvalue of the (original) sample covariance of the function gradient at initialization ĈSf (ϑ0 ) achieves a minimum at N = p, in the limit where N, p → ∞ and p/N → γ. Proof. We relate the eigenvalues of both sample covariance matrices: the whitened one and the original one. In both cases, the smallest non-zero eigenvalues of the sample covariance matrix λmin (C̃Sf (ϑ0 )) correspond to the smallest non-zero singular 2 values of the matrix of function gradients σmin (Z̃S (ϑ0 )) (and analogously for ĈSf (ϑ0 ) and ZS (ϑ0 )). Using the identity Z̃S (ϑ0 ) = Cf (ϑ0 )−1/2 ZS (ϑ0 ), standard singular value inequalities (see, e.g., Ref. [54]) imply σmin Z̃S (ϑ0 ) σmin Cf (ϑ0 )1/2 ≤ σmin (ZS (ϑ0 )) ≤ σmin Z̃S (ϑ0 ) σmax Cf (ϑ0 )1/2 . (D17) p By Corollary D.10, σmin Z̃S (ϑ0 ) converges to 1 − p/N with high probability, which is σmin Z̃S (ϑ0 ) → 0 in the stated limit of p = N . From Assumption 5 it then follows thatσmin (ZS (ϑ0 )) → 0 in the same limit, thus achieving a minimum. Lemma D.12 (Adapted from Ref. [55]). Let K = 1, p ≥ N , and ϑ̂S a local minimum as in Definition A.1. Under Assumption 3, and considering the mean-squared error loss function, ϑ̂S is a global minimum of the empirical risk h i ϑ̂S ∈ arg min L̂S (ϑ) . (D18) ϑ

17 Proof. We start with the definition of local minima: ∇ϑ L̂S (ϑ) = 0. From the definition of the empirical risk, this becomes PN i=1 ∇ϑ ℓ(zi , ϑ̂S ) = 0. We first use the chain rule ∇ϑ ℓ(zi , ϑ̂S ) = ∇ϑ fϑ̂S (xi )∇f ℓ(zi , ϑ̂S ), and we note that, for the meansquared error loss ℓ(z, ϑ) = (fϑ (x) − y)2 /2, we have ∇f ℓ(z, ϑ) = fϑ (x) − y. We introduce the vector of errors ÊS = (fϑ̂S (xi ) − yi )N i=1 , and we re-write the sum over xi ∈ S as a matrix-vector multiplication, in order to get N X

∇ϑ ℓ(zi , ϑ̂S ) =

i=1

N X

∇ϑ fϑ̂S (xi )∇f ℓ(zi , ϑ̂S ) =

i=1

N X

∇ϑ fϑ̂S (xi )(fϑ (xi ) − yi ) = ZS (ϑ̂S )ÊS .

(D19)

i=1

With this, the condition that ϑ̂S is a local minimum becomes ZS (ϑ̂S )ÊS = 0. Assumption 3 dictates that, in the overparameterized regime p ≥ N , the sample Jacobian matrix fulfills rank(ZS (ϑ̂S )) = N (is full rank). From this, it follows that the only possibility for ZS (ϑ̂S )ÊS = 0 to hold is that ÊS = 0. Since ÊS is the vector of errors, it being zero means that the model makes no errors on the training set, which by definition implies the empirical risk achieves its lowest possible value 0. This confirms that ϑ̂S is a global minimum of L̂S . Theorem C.3 (Double descent peak at interpolation). Let S = {(xi , yi )}N i=1 be a dataset with d-dimensional Gaussian normal inputs xi ∼ N (0, Id ) and outputs yi ∈ R. Let fϑ : Rd → R be a parameterized quantum circuit as in Eq. (A4) for K = 1, consisting of p parameters. Under Assumptions 1–7, and considering the mean-squared error loss function, the lower bound on L(ϑ̂S ) at a local minimum from Theorem B.9 attains a maximum at p = N with high probability in the limit N, p → ∞ with fixed ratio p/N → γ. Proof. Let us start from the lower bound in Theorem B.9, given by L(ϑ̂S ) ≥ L′ (S) +

2 λmin (Cf (ϑ̂S )) 1 αξmin + O N −2 . N +1 λ Ĥ (ϑ̂ ) min

(D20)

f o We first turn our attention to the decomposition of ĤS (ϑ̂S) from Lemma D.1: ĤS (ϑ̂S ) = ĤS (ϑ̂S ) + ĤS (ϑ̂S ). Then, we

apply the upper bound from Assumption 6: λmin ĤS (ϑ̂S ) ≤ λmin ĤSo (ϑ̂S ) + λmax ĤSf (ϑ̂S ) . Also, specializing to the mean-squared error loss simplifies the expression for the outer-product Hessian, since the first derivative yields ∇f ℓ(zi , ϑ) = (fϑ (xi ) − yi ) and the second derivative ∇f (fϑ (xi ) − yi )) = 1. With these, we have ĤSo (ϑ) :=

N N h i 1 X 1 X ∇ϑ fϑ (xi ) ∇f ∇⊺f ℓ(zi , ϑ) ∇ϑ fϑ (xi )⊺ = ∇ϑ fϑ (xi )∇ϑ fϑ (xi )⊺ =: ĈSf (ϑ). N i=1 N {z } | i=1

(D21)

The outer-product Hessian becomes the uncentered sample covariance of the function gradient. Substituting both this and the inequality from Assumption 6 into the lower bound, we obtain L(ϑ̂S ) ≥ L′ (S) +

1 N +1λ

min

2 αξmin λmin (Cf (ϑ̂S )) + O N −2 . ĈSf (ϑ̂S ) + λmax ĤSf (ϑ̂S )

(D22)

With Assumption 4, we can further replace the ϑ̂S dependence by the value of the parameters at initialization ϑ0 , for some constants A, B > 0, to get L(ϑ̂S ) ≥ L′ (S) +

1 N + 1 Bλ

αξ 2 Aλ (C (ϑ )) min min f 0 + O N −2 . f f min ĈS (ϑ0 ) + λmax ĤS (ϑ̂S )

(D23)

h i First of all, recall that the first term in both regimes, L′ (S) = Ez′ ∼D ℓ(z ′ , ϑ̂S ′ ) , denotes the add-one-in loss. Since it depends on a training set perturbed by only a single sample, we do not expect its variation with p to dominate the Hessian-dependent complexity term near interpolation. Let us now consider the overparameterized regime p ≥ N . From Assumption 3 and Lemma D.12, for our choice of mean f squared error loss, we have that ĤS (ϑ̂S ) = 0. The denominator in the lower bound then becomes simply Bλmin ĈSf (ϑ̂S ) . Wielding Lemma D.11 we note that λmin ĈSf (ϑ̂S ) achieves a minimum at p = N , and hence it is an increasing function of p for p ≥ N . With this, we confirm that the lower bound of L(ϑ̂S ) is decreasing in the overparameterized regime.

18 In the underparameterized regime p < N , the local minimum ϑ̂S is not necessarily a global minimum of the empirical risk, and hence both summands in the denominator of the lower bound remain. In this case, we combine Assumption 2 and f Lemma D.11 to argue that, with p < N , both summands are decreasing functions of p. For λmax ĤS (ϑ̂S ) , Assumption 2 explicitly states this fact; and for λmin ĈSf (ϑ0 ) , it follows λmin ĈSf (ϑ0 ) having a minimum at p = N . We thus confirm that the lower bound of L(ϑ̂S ) is increasing in the underparameterized regime near interpolation. The last two paragraphs combined complete the proof: the lower bound of L(ϑ̂S ) has a maximum at p = N . Appendix E: Upper bound

In this section we complement Theorems B.9 and C.3 with analogous results for an upper bound on the expected risk. We use the same notation and language as in Appendices B and C. Theorem E.1 (Upper bound on the expected risk). Under the setting of Lemma B.5, under Assumption 1 and for K = 1, the expected risk L(ϑ̂S ) at the local minimum ϑ̂S fulfills the following upper bound: 2 R λmax (Cf (ϑ̂S )) 1 1 α ξmax ′ . (E1) +O L(ϑ̂S ) ≤ L (S) + N +1 N2 λmin (ĤS (ϑ̂S )) 2 = max z∼D {ξz2 }, R = rank(ĤS (ϑ̂S )), λmin (·) the smallest non-zero eigenvalue of its matrix argument, Here, we define ξmax ξz2 >0

and analogously for the largest eigenvalue λmax (·). Proof. The proof follows the steps of the proof of Theorem B.9, using the opposite direction of Lemma B.6, h i−1 1 1 L(ϑ̂S ) = L′ (S) + C(ϑ̂S ) + O Tr ĤS (ϑ̂S ) N +1 N2 h i−1 1 1 2 Cf (ϑ̂S ) +O α ξmax Tr ĤS (ϑ̂S ) ≤ L′ (S) + N +1 N2 i−1 h 2 1 α ξmax +O λmax (Cf (ϑ̂S )) Tr ĤS (ϑ̂S ) ≤ L′ (S) + N +1 N2 ! 2 1 α ξmax λmax (Cf (ϑ̂S )) R +O ≤ L′ (S) + N +1 N2 λmin (ĤS (ϑ̂S )) 2 1 R λmax (Cf (ϑ̂S )) 1 α ξmax +O = L′ (S) + . N +1 N2 λmin (ĤS (ϑ̂S ))

(E2) (E3) (E4) (E5) (E6)

The only qualitative difference in this expression is in the trace of the inverse of the Hessian. This time, we upper bound a sum of positive terms by the largest of the terms (λmin (ĤS (ϑ̂S )))−1 times the number of terms R. This concludes the proof. Assumption 2’. In the underparametrized regime p < N K, the smallest eigenvalue λmin (ĤSf (ϑ̂S )) of the functional Hessian at local minima decreases with the number of trainable parameters p in a neighborhood of the interpolation threshold p = N K. Assumption 4’. The largest eigenvalue of the uncentered covariance of the function gradients remains comparable to its value at initialization, uniformly in the number of parameters p, λmax (Cf (ϑ̂S )) ≍ λmax (Cf (ϑ0 )). This relation holds with high probability with respect to the random draw of the training set S, the initialization ϑ̂0 , and the local minimum ϑ̂S returned by the training procedure. Theorem E.2 (Upper bound peak at interpolation). Let S = {(xi , yi )}N i=1 be a dataset with d-dimensional Gaussian normal inputs xi ∼ N (0, Id ) and outputs yi ∈ R. Let fθ : Rd → R be a parameterized quantum circuit as in Eq. (A4) for K = 1, consisting of p parameters. Under Assumptions 1, 2’, 3, 4, 4’, and 5–7, and considering the mean-squared error loss function, the upper bound on L(ϑ̂S ) at a local minimum from Theorem E.1 attains a maximum at p = N with high probability in the limit N, p → ∞ with fixed ratio p/N → γ. Proof. The proof follows the steps of the proof of Theorem C.3. We start from the upper bound in Theorem E.1, which is 2 1 1 α ξmax R λmax (Cf (ϑ̂S )) +O . (E7) L(ϑ̂S ) ≤ L′ (S) + N +1 N2 λmin (ĤS (ϑ̂S ))

λmax (Hf )

N=21 N=30 N=39 N=48

0.3 0.2

(b) 400

0.4

Jacobian rank

(a)

300

200

100

0.1 100

200 300 Number of parameters

400

200 400 Number of parameters

After training Initialization

10−1 10−4 10−7

600

200 400 Number of parameters

600

Figure F.1. Numerical verification of the assumptions underlying the double descent analysis for the MNIST-1D dataset, for different training set sizes N . Vertical dashed lines indicate the predicted interpolation thresholds p = N K. (a) Largest eigenvalue of the functional Hessian ĤSf (ϑ̂S ) after training. In the underparameterized regime p < N K, λmax (ĤSf (ϑ̂S )) decreases as p approaches interpolation, supporting Assumption 2. (b) Rank of the sample Jacobian of the function ZS (ϑ̂S ) after training. For p ≥ N K, the rank saturates at N K, supporting Assumption 3. (c) Smallest non-zero eigenvalue of the uncentered sample covariance matrix of function gradients ĈSf at initialization and after training, supporting Assumption 4. The shaded areas correspond to the standard deviation for ten independent experiment repetitions, each using independently sampled training data.

We first turn our attention to the decomposition of ĤS (ϑ̂S ) from Lemma D.1: ĤS (ϑ̂S ) = ĤSo (ϑ̂S ) + ĤSf (ϑ̂S ). Then, we apply the other direction of the inequality in Assumption 6, which is λmin ĤS (ϑ̂S ) ≥ λmin ĤSo (ϑ̂S ) + λmin ĤSf (ϑ̂S ) . (E8) Again, specializing to the mean-squared error loss simplifies the expression for the outer-product Hessian, we have ĤSo (ϑ) = ĈSf (ϑ). Substituting both this and the inequality from Assumption 6 into the upper bound, we obtain L(ϑ̂S ) ≤ L′ (S) +

1 N +1λ

αξ 2 Rλ (C (ϑ̂ )) + O N −2 . max max f S f f min ĈS (ϑ̂S ) + λmin ĤS (ϑ̂S )

(E9)

With Assumptions 4 and 4’, we can further replace the ϑ̂S dependence by the value of the parameters at initialization ϑ0 , for some constants A, B > 0, L(ϑ̂S ) ≤ L′ (S) +

1 N + 1 Aλ

2 RBλmax (Cf (ϑ0 )) αξmax + O N −2 . f f min ĈS (ϑ0 ) + λmin ĤS (ϑ̂S )

(E10)

Let us now consider the overparameterized regime p ≥ N . From Assumption 3 and Lemma D.12, for our choice of mean f squared error loss, we have that ĤS (ϑ̂S ) = 0. The denominator in the lower bound then becomes simply Aλmin ĈSf (ϑ̂S ) . Wielding Lemma D.11 we note that λmin ĈSf (ϑ̂S ) achieves a minimum at p = N , and hence it is an increasing function of p for p ≥ N . With this, we confirm that the upper bound of L(ϑ̂S ) is decreasing in the overparameterized regime. In the underparameterized regime p < N , the local minimum ϑ̂S is not necessarily a global minimum of the empirical risk, and hence both summands in the denominator of the upper bound remain. In this case, we combineAssumption 2’ and Lemma D.11 to argue that, with p < N , both summands are decreasing functions of p. For λmin ĤSf (ϑ̂S ) , Assumption 2’ explicitly states this fact; and for λmin ĈSf (ϑ0 ) , it follows from λmin ĈSf (ϑ0 ) having a minimum at p = N . We thus confirm that the upper bound of L(ϑ̂S ) is increasing in the underparameterized regime near interpolation. The last two paragraphs combined complete the proof: the upper bound of L(ϑ̂S ) has a maximum at p = N . Appendix F: Further empirical evidence of double descent in PQCs

Complimenting our main numerical findings reported in Fig. 1, we now turn our attention to a selection of the assumptions we introduced for Theorem C.3 in Appendix C. We present our results in Fig. F.1, where we probe further the mechanism

20 behind the double descent peak. Here, we focus on the MNIST-1D dataset. By systematically testing for different assumptions, we strengthen the evidence that the peak at interpolation in Fig. 1 is accurately captured by our main analytical contribution, Theorem C.3. • In Fig. F.1(a) we test Assumption 2: in the underparameterized regime p < N K, the largest eigenvalue of the functional Hessian, λmax (ĤSf (ϑ̂S )), decreases as the number of parameters approaches the interpolation threshold. This indicates that the contribution of the functional Hessian becomes less dominant near interpolation, as required in the proof of the double descent peak. • In Fig. F.1(b) we test Assumption 3: we show that, once p ≥ N K, the rank of the sample Jacobian of the function ZS (ϑ̂S ) saturates at N K for all training-set sizes considered. Thus, in the overparameterized regime, the Jacobian has full rank with respect to the N K scalar training constraints. • In Fig. F.1(c) we test Assumption 4: we compare the smallest non-zero eigenvalue of the uncentered sample covariance matrix of function gradients ĈSf at initialization and after training, displaying a dip near the interpolation threshold. Appendix G: Further details on the numerical experiments

In this section, we provide further implementation details for the numerical experiments presented in Figs. 1 and F.1. All quantum circuit simulations were performed with the PennyLane software library [56], and all experiments use n = 8 qubits and K = 8 output dimensions. The circuit consists of an initial angle embedding of the classical input, followed by L data re-uploading layers. Each layer contains trainable single-qubit rotations RX , RZ , and RY on every qubit, followed by nearestneighbor CZ gates arranged in a ring. The input is then re-uploaded before the next trainable layer. The model output is obtained by measuring Pauli-X expectation values on the first K qubits and multiplying the resulting vector by the fixed scale factor c = 150. Since each layer contains three trainable parameters per qubit, the total number of trainable parameters is p = 3nL. In the experiments, the depth is varied as L = 2, . . . , 25, corresponding to p = 48, . . . , 600 trainable parameters. All models are trained with the mean squared error loss using the Adam optimizer, and the training is performed for 2500 epochs. For the MNIST-1D and Fashion MNIST classification tasks, we use 8 classes. The labels are encoded as one-hot vectors and subsequently mapped from {0, 1} to {−1, 1}, so that the correct class has value +1 and all other classes have value −1. This allows the classification tasks to be trained using the same mean squared error objective as in the theoretical analysis. Classical inputs are reduced to dimension d = n = 8 using principal component analysis, implemented with scikit-learn [57], and then rescaled component-wise to the interval [−π/2, π/2] before being passed to the angle embedding. For the synthetic regression task, inputs are sampled uniformly from [−π/2, π/2]8 , and targets are generated from a multidimensional linear model y = XW + ζ, where W is the all-ones matrix in R8×8 and ζ is Gaussian noise with standard deviation 0.5. To scan across the interpolation threshold, we vary the number of training samples as N ∈ {21, 30, 39, 48}. Since K = 8, the corresponding interpolation thresholds are p = N K ∈ {168, 240, 312, 384}. For each value of N and each circuit depth, the test loss is evaluated on 1000 test samples. These values are chosen such that the interpolation thresholds coincide exactly with parameter counts attainable by varying the number of data re-uploading layers. The results are averaged over 10 independent repetitions, and shaded regions in the figures indicate the standard deviation over these repetitions. The code is publicly available in Ref. [50].

Appendix H: Universality classes of double descent

The proof of Theorem C.3 establishes the double descent peak by relating the smallest non-zero eigenvalue of the empirical covariance matrix of function gradients to the soft edge of the Marčenko–Pastur law. This naturally raises the question of how much of the phenomenon depends on the specific Wishart ensemble arising from independent Gaussian data, and how much is instead a manifestation of a broader random-matrix universality principle. Viewed abstractly, the proof relies on only two ingredients. First, the smallest non-zero singular value of the empirical Jacobian develops a soft edge at the interpolation threshold. Second, the associated singular vectors are sufficiently delocalized so that generic data directions overlap with them in a statistically uniform way. The specific Gaussian assumptions merely provide one convenient realization of these properties. A natural question is to what extent the present analysis depends on the specific Wishart ensemble underlying the Marčenko– Pastur law. From the perspective of the proof, the essential ingredient is not Gaussianity itself but rather the emergence of a soft spectral edge governing the smallest non-zero singular value of the empirical covariance matrix. This suggests that analogous double descent behavior may persist for considerably broader classes of covariance ensembles. Prominent examples include

21 correlated Wishart ensembles, random kernel matrices, and deformed covariance ensembles [58–60]. More generally, nonlinear kernel models whose limiting spectral measures differ from the classical Marčenko–Pastur law may provide a mathematically natural framework for extending the present analysis beyond linearized parameterizations [28, 59]. These observations suggest that the mechanism underlying Theorem C.3 should be understood as a manifestation of a broader universality class rather than of a single random matrix model. Conjecture H.1 (Universality of double descent). Consider the setting of Theorem C.3, replacing the Gaussian Wishart covariance model by a random covariance ensemble whose smallest non-zero singular value exhibits universal soft-edge behavior and whose singular vectors remain asymptotically delocalized. Then the lower bound on the expected risk exhibits an interpolation peak at the corresponding critical aspect ratio, giving rise to double descent behavior. Establishing the precise universality class governing double descent in parameterized quantum circuits, and identifying which structural properties of the underlying random matrix ensemble are genuinely required, constitute interesting open problems.

Record · ID 394409 · SHA-256 9cae1ba2d3cb4674
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.