Year: 2022

Abstract:
One of the most important properties of neural networks is the clustering of local minima of the loss function near the global minimum, enabling efficient training. Though generative models implemented on quantum computers are known to be more expressive than their traditional counterparts, it has empirically been observed that these models experience a transition in the quality of their local minima. Namely, below some critical number of parameters, all local minima are far from the global minimum in function value; above this critical parameter count, all local minima are good approximators of the global minimum. Furthermore, for a certain class of quantum generative models, this transition has empirically been observed to occur at parameter counts exponentially large in the problem size, meaning practical training of these models is out of reach. Here, we give the first proof of this transition in trainability, specializing to this latter class of quantum generative models. We use techniques inspired by those used to study the loss landscapes of classical neural networks. We also verify that our analytic results hold experimentally even at modest model sizes.

heading: 1 INTRODUCTION

heading: 1.1 MOTIVATION
text: One of the great successes of neural networks is the efficiency at which they are trained via gradientbased methods. Though training algorithms often involve the optimization of complicated, nonconvex, high-dimensional functions, training via gradient descent in many contexts manages to converge to local minima that are good approximations of the global minimum in loss function value. This phenomenon has begun to be understood in the context of random matrix theory, particularly when applied to dense classifiers (Choromanska et al., 2015; Chaudhari & Soatto, 2017).
Quantum generative models hold great promise in their ability to sample from probability distributions out of the reach of classical models (Arute et al., 2019; Gao et al., 2021). Though deep quantum generative models are believed to be difficult to train due to vanishing gradients (McClean et al., 2018; Cerezo et al., 2021; Marrero et al., 2020), the situation for shallow models is murkier. Some shallow quantum models have empirically shown great promise in being trainable (Wiersema et al., 2020; Kim et al., 2020; Kim & Oz, 2021), while others have empirically been shown to suffer from poor distributions of local minima (Kiani et al., 2020; Campos et al., 2021). Numerically, all of these models have been seen to experience a phase transition in trainability: below some critical depth, local minima are poor approximators of the global minimum. Above this critical depth, they are good approximators. This transition has been poorly understood analytically, as typically the distribution of local minima monotonically improves as the size of the model increases (Choromanska et al., 2015; Chaudhari & Soatto, 2017).

heading: 1.2 OUR CONTRIBUTIONS
text: In this work, we are the first to analytically show the presence of a computational phase transition in the training of a certain class of quantum generative models. To achieve this, we first show that this class of randomized quantum models is approximated in distribution by a Wishart random field
on the hypertorus. We are then able to use techniques from Morse theory to exactly calculate the distribution of local minima (and general critical points) of this random field. Finally, we analyze this distribution in the limit of large model size, and analytically show the presence of this trainability phase transition. Roughly, we show that in this limit the expected density of local minima for a model with p parameters and Hilbert space dimension ∼ m (exponential in the problem size) at loss value 0 ≤ E ≤ 12 follows a generalization of the beta distribution (Gordy, 1998):
E [Crt0 (E)] ∼ e−mEEm− p 2 (1− 2E)p . (1)
This distribution experiences a transition in behavior at p = 2m: when p < 2m, local minima are exponentially concentrated (i.e. with width m−1) far away from the global minimum E = 0, implying poor optimization performance in this regime. When p ≥ 2m, this distribution is exponentially concentrated at E = 0, implying good optimization performance. We also verify our results numerically, demonstrating this concentration of minima even at small problem sizes.
For the class of quantum generative models we consider, our results mirror the empirical results of Kiani et al. (2020); Campos et al. (2021) in that only unreasonably overparameterized quantum models have good local minima. Though these results are pessimistic, we emphasize here that our results only apply to a certain class of quantum generative models. We are also able to give a heuristic explanation based on our proof techniques as to how one may be able to construct models of a reasonable size that are still trainable at the expense of computational overhead in implementing the model, as seen empirically in Wiersema et al. (2020); Kim et al. (2020); Kim & Oz (2021).

heading: 2 PRELIMINARIES

heading: 2.1 QUANTUM GENERATIVE MODELS
text: A quantum system of size n is naturally represented by a quantum state, which is a normalized vector |ψ〉 ∈ C2n . Here, we use the typical physics notation |·〉 to denote a vector, instead of (say)ψ, when we are describing a quantum state. A quantum state in C2n can be considered a generalization of probability distributions over 2n states (i.e. over states described by n bits), where the norm squared of entries of |ψ〉 give the distribution. The task in quantum generative modeling is to prepare a state |ψ〉 that is close under some metric to a given target state |ψtarget〉. As operations that map probability distributions to probability distributions are naturally described by stochastic matrices, operations that map quantum states to quantum states are naturally described by unitary matrices; equivalently, they are described by the matrix exponentials of 2n × 2n skewHermitian matrices. Multiple quantum operations can then be described by the sequential matrix multiplication of various matrix exponentials. It is then apparent that one can construct a quantum generative model by parameterizing these matrix exponentials, giving rise to the layered structure:
|θ〉 ≡ q∏ i=1 Ui (θ) |ψ0〉 ≡ q∏ i=1 e−iθiQi |ψ0〉 , (2)
where here q is the depth of the model and Qi are fixed Hermitian matrices. Typically, q is polynomial in n, i.e. logarithmic in the dimension of the initial state vector |ψ0〉. In most quantum generative models, various θi are completely dependent, e.g. θi = θi+5 for all i (Romero et al., 2018). For simplicity, we assume throughout this work that each independent parameter appears a constant number r times in the model, and the total number of independent parameters is given by p = q/r. Though the linear nature of the model may be surprising, this model structure (for large enough depth q) is known to be a universal approximator of all (pure) quantum states, even for a fixed number of allowed Qi (Solovay, 1995; Kitaev, 1997).
For completely general (i.e. dense) Qi in equation 2, the model is not efficient to implement on a quantum computer. Thus, due to their efficiency in physical implementation, these Qi are typically taken to be members of the n-qubit Pauli group Pn, which is a generating subset of the additive group of all n × n Hermitian matrices. The Pauli group is also convenient to study analytically, as is the normalizer of the group (called the Clifford group). The assumption that each Qi is a Pauli operator also allows us to use a single parameter θi for each layer of the model without loss of generality; in principle, more parameters can describe each layer by parameterizing sparse Qi. However, as the
Pauli group generates Hermitian matrices, this is a special case of the class of models we consider here (at the expense of larger q and new dependencies among the θi).
In quantum generative modeling, we are often not given the target state |ψtarget〉 directly. Instead, we are in many cases given a 2n × 2n Hermitian matrix H (called the problem Hamiltonian) where |ψtarget〉 is the eigenvector associated with the smallest eigenvalue of H (called the ground state). In this formulation, optimization proceeds via the minimization of:
FVQA (θ) = 〈θ|H |θ〉 . (3)
In physics notation, 〈θ| is the conjugate transpose of the complex vector |θ〉. Assuming no degeneracies in the eigenspectrum ofH and a sufficiently expressive model, the minimizer |θ∗〉 of equation 3 is the ground state of H , up to an overall phase due to the quadratic nature of equation 3. Assuming H is efficiently expressible as the weighted sum of O (poly (n)) Pauli matrices, equation 3 and its gradients can be efficiently measured on a quantum computer (Peruzzo et al., 2014; Romero et al., 2018). This and similar formulations of quantum generative modeling with equation 3 as the loss function are called variational quantum algorithms (VQAs) (Peruzzo et al., 2014). Generally, the goal of these algorithms is to find the state |θ〉 that optimizes equation 3, up to potentially some constant additive error in loss. Though there are other formulations of quantum generative modeling, we here focus on VQAs as they do not require coherent access to data |ψtarget〉, which is generally believed to be difficult (Aaronson, 2015).
Typically, models in VQAs come in one of two flavors: Hamiltonian agnostic models, and Hamiltonian informed models. Hamiltonian agnostic models are constructed such that the Qi present in the model definition are independent of H , and are generally more efficient to implement. This is most analogous to the case in classical generative modeling, where the model structure is usually independent from the specific choice of data distribution. We prove the existence of a computational phase transition for a class of Hamiltonian agnostic models in this work. We also give some heuristic and numerical evidence that a similar transition exists for Hamiltonian informed ansatzes, but that the trainable phase is more easily reached.

heading: 2.2 RANDOM FIELDS ON MANIFOLDS
text: As in previous results on the loss landscapes of machine learning models (Choromanska et al., 2015; Chaudhari & Soatto, 2017), we will map the distribution of a randomized class of quantum generative models to a random field on a manifold. This will then enable us to use standard mathematical techniques to study the distribution of critical points of the model.
Though they can be expressed in many ways, here we will be interested in random fields of the form:
FRF (σ) ∝ Λ∑
i1,...,ir,i′1,...,i ′ r=1
σi1 . . . σirJi1,...,ir,i′1,...,i′rσi′1 . . . σi′r . (4)
Here, σ ∈M is some point on a manifold, and J is a random variable. In the context of most studies of machine learning loss landscapes, M is typically the hypersphere and J a symmetric matrix of i.i.d. Gaussian random variables, i.e. a Gaussian orthogonal ensemble (GOE) matrix.
We will instead find that VQAs are naturally described as Wishart hypertoroidal random fields (WHRFs). For these models, the manifold M is a tensor product embedding of the hypertorus into exponentially large Euclidean space; that is, points on this embedding are described by the Kronecker product:
σ = ⊗ i ( cos (θi) sin (θi) ) (5)
for angles −π ≤ θi ≤ π. Furthermore, in these models, J is drawn from a normalized complex Wishart distribution. The complex Wishart distribution is a natural multivariate generalization of the gamma distribution, and is given by the distribution of the square of a complex Gaussian random matrix. Specifically, forX ∈ Cn×m a matrix with i.i.d. complex Gaussian columns with covariance matrix Σ, the matrix
W = 1
m X ·X† (6)
is normalized complex Wishart distributed with scale matrix Σ and m degrees of freedom. We will find that the degrees of freedom m will greatly affect the distribution of local minima of the WHRF, and thus also of the class of quantum generative models that we consider.

heading: 3 QUANTUM GENERATIVE MODELS AS RANDOM FIELDS
text: We first show that a certain randomized class of Hamiltonian agnostic VQAs can be expressed as WHRFs. This will allow us to more easily study the critical points of the model using techniques from random matrix theory. Though we leave the full statement and proof for Appendix A, we give an informal statement and discussion here.
Theorem 1 (VQAs as WHRFs, informal). Consider the class of models
|θ〉 ≡ q∏ i=1 Ui (θ) |ψ0〉 ≡ q∏ i=1 e−iθiQi |ψ0〉 , (7)
where eachQi is drawn uniformly from the Pauli group Pn and |ψ0〉 is the first column of a uniformly random member of the Clifford group. Let p be the number of distinct θi, and let r = q/p. Under reasonable assumptions on the eigenvalues ofH (with minimum eigenvalue λ1 and mean eigenvalue λ), the random loss function
FH (θ) = 〈θ|H |θ〉 − λ1
λ− λ1 (8)
converges in distribution to the random field
FWHRF (θ) = 2p∑ i1,...,ir,i′1,...,i ′ r=1 wi1 . . . wirJi1,...,ir,i′1,...,i′rwi′1 . . . wi′r , (9)
where w are points on the hypertorus ( S1 )×p
parameterized by θ and J is a complex Wishart random matrix normalized by its number of degrees of freedom m.
Note that for convenience, we have shifted and scaled the typical VQA loss such that it is always greater than zero and independent from overall scalings of the problem Hamiltonian.
In the course of this mapping, we find that the degrees of freedom of the Wishart matrix J (formally a real number) is given by the ratio:
m ≡ ‖H − λ1‖2∗∥∥H − λ∥∥2F . (10)
Here, ‖·‖∗ denotes the nuclear norm, and ‖·‖F the Frobenius norm. Generally, this ratio is exponential in n, particularly when modeling the class of ground states typically represented by VQAs (Peruzzo et al., 2014; Farhi et al., 2014; Romero et al., 2018). Though we are unable to prove Theorem 1 for Hamiltonian informed models, there are heuristic reasons to believe that they are described by a similar random field with m = O (poly (n)), as opposed to equation 10 (see Appendix A.3 and empirical evidence in Sec. 5). We will later find that a number of independent model parameters p that is twice the degrees of freedom m of the matrix J marks the transition from the underparameterized to the overparameterized regime of FWHRF, where the quality of local minima improves.
The general idea for showing this equivalence relies on the path integral expansion of the VQA loss function. Effectively, this is just a Taylor expansion of the unitary matrices composing the model, which is exact even at a finite number of terms. One can then show that terms in this expansion can be assumed independent with negligible error in distribution, and then show that the resulting random process is asymptotically a WHRF. The reasonable assumptions on the eigenvalues of H are essentially just a requirement that the eigenvalues of H are not “unnaturally” spread out; for the quantum states VQAs typically model, this is never the case. We give a full description of these requirements with the full proof in Appendix A.

heading: 4 THE LOSS LANDSCAPE OF WISHART HYPERTOROIDAL RANDOM FIELDS

heading: 4.1 EXACT RESULTS
text: Having shown that VQAs can be described as WHRFs, we now focus discussion entirely on WHRFs. Our strategy for showing the distribution of critical points of this random field will be similar to that in Auffinger et al. (2013), where similar results were shown for Gaussian spherical random fields. Namely, we will lean heavily on the Kac–Rice formula, which gives the expected number of critical points of a certain index at a given range of function values for random fields on manifolds. We give an informal description of the Kac–Rice formula here, with the formal version given in Appendix B. Lemma 1 (Kac–Rice formula (Adler & Taylor, 2009), informal). Let M be a compact, oriented manifold. Assume a random field F (σ) on M is sufficiently nice. Then, the number of critical points of index at most k with F (σ) ∈ B for an open set B ⊂ R is
E [Crtk (B)] = ∫ M E [∣∣det (∇2F (σ))∣∣1 {F (σ) ∈ B}1{ι (∇2F (σ)) ≤ k} |∇F (σ) = 0] × pσ (∇F (σ) = 0) dσ ,
(11) where ∇· is the covariant gradient, ι (·) is the index of ·, pσ is the probability density of ∇F (σ) at σ, and dσ is the volume element on M .
From Lemma 1, we see that when the joint distribution of ∇2F , ∇F , and F is known, then the expected number of critical points with function values in an open set B can be calculated. Perhaps surprisingly, as in the Gaussian case, the joint distribution of these derivatives for WHRFs is fairly simple. Once again leaving the full proof for Appendix C, we show the distribution of the Hessian conditioned to be at a critical point of function value x can be described by the shifted sum of a Wishart matrix with an independent GOE matrix. Similarly, the distribution of the gradient conditioned on the function value being x is given by a normal distribution. Lemma 2 (Hessian and gradient distributions, informal). The scaled Hessian m∂i∂jFWHRF (w) conditioned on FWHRF (w) = x and ∂kFWHRF (w) = 0 is distributed as
mC̃ij (x) = −2rmxδij + rWij + r √ 2mxNij , (12)
where W is Wishart distributed with 2m degrees of freedom, N GOE distributed, and they are independent. Furthermore, the scaled gradient m∂kFWHRF (w) conditioned on FWHRF (w) = x is distributed as
mG̃k (x) = √
2mrxNk, (13) where Nk are i.i.d. standard normal distributions independent from all Wij and Nij .
With all of the pieces in place, we are able to explicitly calculate the expected distribution of local minima in WHRFs via the Kac–Rice formula (with full calculations left for Appendix C). In Sec. 5 we find empirical evidence that these results hold not only in expectation, but in distribution; we leave further analytic investigation of this to future work. Theorem 2 (Distribution of critical points in WHRFs, informal). Let
µC(x) = 1
p p∑ i=1 δ ( λ C(x) i ) (14)
be the empirical spectral measure of the random matrix
C (x) = r
m
( W + √ 2mxN ) , (15)
where W is Wishart distributed with 2m degrees of freedom, N GOE distributed, and they are independent. λCi (x) is the ith smallest eigenvalue of C (x). Then, the distribution of the expected number of critical points of index k of a WHRF at a function value E > 0 is given by
E [Crtk (E)] = (π r ) p 2 Γ (m) −1 m(1+γ)mEC(E) [ ep ∫ ln(|λ−2rE|)dµC(E)1 { λ C(E) k+1 ≥ 2rE }] E(1−γ)m−1e−mE ,
(16)
where γ = p
2m . (17)
We call the parameter γ the overparameterization factor. It describes the ratio between the number of parameters of the model p and twice the degrees of freedom m of the model. We will later find that γ governs the phase transition between an underparameterized phase of the model—where local minima are far from the global minimum—and an overparameterized phase—where local minima are good approximators of the global minimum.
4.2 ASYMPTOTIC RESULTS AS p→∞
Though Theorem 2 gives the exact distribution of critical points, it is difficult to use in practice. This difficulty comes from the expectation over eigenvalues of the sum of independent Wishart and Gaussian matrices. Surprisingly, however, the eigenvalues of both Wishart and Gaussian orthogonal matrices converge to fixed distributions. Essentially, asymptotically in the size of the matrix, the eigenvalue distribution of all normalized Wishart matrices are the same (given by the Marchenko– Pastur distribution) and the eigenvalue distribution of all Gaussian orthogonal matrices are the same (given by the Wigner semicircle distribution).
Luckily, we can characterize the asymptotic distribution of eigenvalues of the sums of these matrices using the tools of free probability theory. Roughly, free probability theory is the probability theory of noncommutative random variables (e.g. random matrices). As the distribution of the sum of two random variables in commutative probability theory can be described by the convolution of the distributions of the two independent random variables, so can the free convolution of the distributions of two freely independent noncommutative random variables. Using the asymptotic free independence of Wishart and Gaussian orthogonal random variables, we are able to show that asymptotically the eigenvalue distribution of their sum weakly converges to the free convolution of a Marchenko–Pastur distribution with a semicircle distribution.
However, weak convergence is not enough; due to the exponential factor in the expectation in Theorem 2, any large deviations from the asymptotic convergence—even if they occur with exponentially vanishing probability—can potentially cause large deviations from the naive application of free probability theory. Thus, our results rely on using large deviations theory to show that to (logarithmic) leading order these deviations do not contribute to the final result. This is due to the contribution to the expectation from the deviations being dominated by what is predicted by free probability theory. These results can be summarized via the following theorem (proved in Appendix D): Theorem 3 (Logarithmic asymptotics of the local minima distribution, informal). Let dµ∗E be the free convolution of a scaled Marchenko–Pastur and scaled Wigner semicircle distribution, with λ∗E,1 the infimum of its support. Let p,m 1 with p2m = γ = O (1). Then, the expected distribution of local minima of a WHRF at a fixed function value E > 0 is given by
1 p ln (E [Crt0 (E)]) = 1 2 ln
( πq
2γ
) + 1
2γ (1− E) + 1 2
( γ−1 − 1 ) ln (E)
+ ∫ ln (∣∣∣∣λr − 2E ∣∣∣∣1{λ∗E,1r ≥ 2E }) dµ∗E + o (1) .
(18)
Note that, though we only prove the asymptotic distribution of local minima in Theorem 3, we expect similar theorems to also hold for critical points of constant index k (taking λ∗E,1 7→ λ∗E,k in the integrand). The only difference in the derivation is the exact form of the large deviations of the kth smallest eigenvalue of C (x). This is similar to the case in Gaussian hyperspherical random fields, which are often used to model neural network loss functions (Auffinger et al., 2013; Choromanska et al., 2015; Chaudhari & Soatto, 2017).

heading: 4.3 DISCUSSION OF THE CRITICAL POINT DISTRIBUTION
text: Let us now discuss the implications of Theorem 3. Note first that the rescaled logarithmic number of critical points diverges when q → ∞. Following the derivation closely, one finds that this is
due to an exponentially suppressed (when m is exponential in n) gradient. We believe that this is a manifestation of the “barren plateau” phenomenon, where for many deep VQA models it can be shown that there is an exponentially vanishing variance of the gradient over the loss landscape (McClean et al., 2018; Cerezo et al., 2021; Marrero et al., 2020). This interpretation suggests that these barren plateau regions are filled with many small “bumps” that are exponentially shallow. Furthermore, note that this class of random fields exhibits banded behavior in the eigenvalues. That is, local minima only exist in the band 0 ≤ E ≤ E0, where E0 is the solution to
λ∗E0,1 = 2rE0. (19)
This banded behavior is similar to that in the Gaussian spherical case. We will see, however, that this does not give necessarily good guarantees on the distribution of local minima. This is due to E0 being generally far from 0 when γ < 1 as p,m→∞. To illustrate this, we focus now on two cases: p ≥ 2m (the overparameterized regime) and p m (the underparameterized regime).
4.3.1 p ≥ 2m
First, let us consider when p ≥ 2m, i.e. γ ≥ 1. In this limit, the Wishart term of C is low-rank, and µ∗E has support on eigenvalues ≤ 0 for all E ≥ 0. Therefore, the condition λ∗E,1 ≥ 2Er is never satisfied, and to leading order in p there are no local minima at a function value E > 0. That is, all local minima are global minima in the p → ∞ limit when γ ≥ 1. Though the choice of model is slightly different, we suspect that a related phenomenon may be what gives rise to the phase transition in training numerically observed in Kiani et al. (2020); Wiersema et al. (2020); Kim et al. (2020); Kim & Oz (2021); Campos et al. (2021).
4.3.2 p m
When the number of distinct parameters p is poly (n) and considering a physically relevant problem Hamiltonian such that the number of degrees of freedom m is exp (n), we have that p m (i.e. γ 1) for large n. In this limit, the spectral distribution µ∗E is dominated by the Wishart term of C, as its eigenvalues are O (1) while the eigenvalues of the GOE term are O (√ γ ) . Furthermore,
the Marchenko–Pastur distribution in this limit only has support at λ = 1 + O (√ γ ) . Therefore, the expected number of local minima at a function value E will be proportional to
E [Crt0 (E)] ∝ e−mE+o(p)Em− p 2 (1− 2E + O (√γ))p 1 { 0 ≤ E + O (√γ) ≤ 1
2
} . (20)
In particular, up to shifts on the order of √ γ, the distribution of local minima is roughly that of a compound confluent hypergeometric (CCH) distribution (Gordy, 1998). The CCH distribution can be considered a generalization of the beta distribution, and for our parameters has mean on the order of 12−γ and standard deviation on the order ofm
−1. Restoring the overall scaling in equation 8, this implies that in this limit the local minima of the variational loss function exponentially concentrate (in expectation) near half the mean eigenvalue of H − λ1 instead of the smallest eigenvalue. Even worse, the CCH distributed form of the local minima implies that, even when beginning at an initial function value well below half of the mean eigenvalue of H − λ1, the found loss will only improve by a fraction of the initial function value before the optimization algorithm finds a local minimum. This is insufficient to find the optimal loss to constant additive error when beginning training at a random point, as is often the goal in VQAs (Peruzzo et al., 2014). Empirically, we find that this occurs not just in expectation but also for individual model instances in Sec. 5.

heading: 5 NUMERICAL EXPERIMENTS
text: We now test our analytic predictions using numerical simulations. First, we investigate the empirical performance of the class of randomized models we study theoretically, and give numerical evidence of things we were unable to prove. Then, we give numerical evidence that, for models dependent on the objective Hamiltonian, the effective degrees of freedom parameter m can be much smaller than predicted. In all cases, we numerically test the predictions of our results by modeling the ground state of the 1D n site spinless Fermi–Hubbard Hamiltonian (Negele & Orland, 1998) at half filling. Here, we take units such that the mean eigenvalue of the considered Hamiltonian (minus its smallest eigenvalue) is E = 1. We give further details of our numerical simulations in Appendix E.

heading: 5.1 EMPIRICAL PERFORMANCE OF RANDOM ANSATZES
text: First, we analyzed the performance of a VQA on this loss function via the random model construction procedure defined in Theorem 1. Previous numerical results on related Hamiltonian agnostic ansatzes have already shown the concentration of local minima far away in loss value from the global minimum below some degrees-of-freedom transition, and concentration at the global minimum above this transition (Kiani et al., 2020; Campos et al., 2021). Here, we tracked where our analysis predicts the local minima to lie as a function of γ for γ 1, up to deviations on the order of √ γ that arise from numerically considering the problem at finite size (as discussed in Sec. 4.3.2).
Concretely, for a given training instance and depth q = p, we generated an ansatz |θ〉 composed of p layers of Pauli rotations, where each Pauli rotation was chosen uniformly from all nonidentity Pauli matrices on n qubits. The numbers of model layers we consider are typical of current physical implementations of Hamiltonian agnostic VQAs (Kandala et al., 2017). A summary of the normalized distribution of found local minima for the randomized model with model dimension 2n = 64, 256 is given in Figure 1, along with the predicted region in which all local minima should lie in the p→∞ limit as discussed in Theorem 3. See Appendix E for details on how this distribution was generated.
We see that almost all found local minima lie within the predicted region, even at small p, n. In particular, for small γ, the distribution of local minima is almost entirely localized within √ γ of the
predicted 12 − γ (in units of the mean eigenvalue of H − λ1). Finally, we numerically observe that the distribution of local minima are qualitatively similar in expectation and for a single choice of random model in Figure 2.

heading: 5.2 EMPIRICAL PERFORMANCE OF A HAMILTONIAN INFORMED MODEL
text: Previous numerical results (Wiersema et al., 2020) on VQAs have shown that only a moderate number of model parameters suffices for efficient training when using a Hamiltonian informed model. As discussed in Sec. 3, we believe this is due to this class of models effectively limiting the degrees of freedom m of the associated WHRF model; to test this, we performed more numerical experiments using a Hamiltonian informed ansatz. We once again tracked where our analysis predicts the local minima to lie as a function of γ for γ 1, up to deviations on the order of√γ.
We show the empirical distribution of local minima in Fig. 3 for 2n = 256, along with the predicted region local minima should lie as discussed in Sec. 4.3.2. The predicted local minima distribution is overly pessimistic (particularly at larger p). We suspect this is due to the fact that the ansatz is constructed in a way that minimizes the effective degrees of freedom of the WHRF m such that γ is close to 1 for smaller p than is predicted analytically.

heading: 6 CONCLUSION
text: Though variational quantum algorithms are perhaps the most promising way to use the error-prone quantum devices of today for practical computational tasks, there are many caveats with regard to their trainability. In particular, previous work has shown that utilizing deep quantum models that are independent of the problem Hamiltonian can introduce a vanishing gradient phenomenon where, though the model is expressive enough to capture the ground state of interest, in practice optimizing the loss function is infeasible (McClean et al., 2018; Cerezo et al., 2021; Marrero et al., 2020). We extended these results by showing a particular class of random models independent of the problem instance not only can exhibit these vanishing gradients at large depth, but also has a concentration of local minima near the mean eigenvalue of the objective Hamiltonian. This is in contrast to the case in traditional neural networks, where even generic model structure tends to lead to a concentration of local minima in a band near the global minimum of the loss function.
Though our results may not seem encouraging for quantum generative models, we emphasize that we expect our analytic results to hold only when the model is independent of the problem Hamiltonian. Indeed, we found empirically (and heuristically) good performance for a particular Hamiltonian informed ansatz, where our analytic results seem much too pessimistic. In principle, this new way of thinking about variational quantum algorithms may inform future quantum generative model design; we leave for future work the study of how various model choices may impact the distribution of critical points of the loss function positively, and how practical considerations such as noisy model implementations may play a role.
REPRODUCIBILITY STATEMENT
The full statement of Theorem 1, along with all of its assumptions and its proof, is given in Appendix A. Similarly, Lemma 2 and Theorem 2 are fully stated and proved in Appendix C. These results rely on the Kac–Rice formula, informally stated as Lemma 1; for completeness, we give the full assumptions of this formula in Appendix B, and there also show the assumptions are met for the model we consider. Furthermore, the full assumptions and proof of Theorem 3 is given in Appendix D, along with the proofs of supplementary lemmas needed to prove the statement. Finally, full details of our numerical simulations and how to reproduce the results are given in Appendix E.

heading: ACKNOWLEDGMENTS
text: We are grateful to Tongyang Li, John Napp, and Aram Harrow for insightful discussion and helpful remarks regarding a draft of this work. E.R.A. is supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. 4000063445, and a Lester Wolfe Fellowship and the Henry W. Kendall Fellowship Fund from M.I.T.

heading: B THE KAC–RICE FORMULA AND ITS ASSUMPTIONS
text: For completeness, we state the formal version of Lemma 1—with all assumptions—here. We borrow heavily from (Adler & Taylor, 2009). By ∇f , we mean the covariant gradient of f . Lemma 3 (Kac–Rice formula (Adler & Taylor, 2009)). Let M be a compact, oriented, N - dimensional C1 manifold with C1 Riemannian metric g. Let B ⊂ RK be an open set such that ∂B has dimension K − 1. Let f : M → RK be a random field on M , and let ι (·) denote the index of ·. Furthermore, assume that:
1. All components of f , ∇f , and ∇2f are almost surely continuous and have finite variances over M .
2. The marginal density pt (∇f (t)) of ∇f at t ∈M is continuous at ∇f = 0. 3. The conditional densities pt ( ∇f (t) | f (t) ,∇2f (t) ) are bounded above and continuous
at ∇f = 0, uniformly in t ∈M . 4. The conditional densities pt ( det ( ∇2f (t) ) |∇f (t) = 0 ) are continuous in the neighbor-
hood of det ( ∇2f ) = 0 and ∇f (t) = 0, uniformly in t ∈M .
5. The conditional densities pt (f (t) |∇f (t) = 0) are continuous for all f and for all ∇f in a neighborhood of 0, uniformly in t ∈M .
6. The Hessian moments are bounded, i.e.
sup t∈M max i,j
E [∣∣∣(∇2f (t))
i,j ∣∣∣N] <∞. (61) 7. The moduli of continuity with respect to (the canonical metric induced by) g of each com-
ponent of f , ∇f , and ∇2f all satisfy P [ω (η) > ] = o ( ηN )
(62)
for all > 0 as η → 0+. Then,
E [ Crtfk (B) ] = ∫ M pσ (∇f (σ) = 0)
× E [∣∣det (∇2f (σ))∣∣1 {f (σ) ∈ B}1{ι (∇2f (σ)) ≤ k} |∇f (σ) = 0]dσ , (63)
where dσ is the volume element induced by g on M .
It is obvious by Lemma 4 that conditions 2-6 are satisfied by WHRFs given B = (0, u). Furthermore, as F is a polynomial in {cos (θi) , sin (θi)}, F and its derivatives are continuous for any value of the components of m−1J , and all have finite variance. Similarly, it is easy to see that the modulus of continuity of f and its gradients go as Jηr as η → 0+, where J is the largest component of m−1J . As the distributions of the components of a Wishart matrix have exponential tails, the probability that J = Ω (η−r) is indeed o ( ηN ) and therefore all conditions are satisfied by WHRFs.

heading: C THE LOSS LANDSCAPE OF WISHART HYPERTOROIDAL RANDOM FIELDS
text: C.1 THE JOINT DISTRIBUTION OF FWHRF AND ITS DERIVATIVES
In order to utilize the Kac–Rice formula (Lemma 3), we must calculate the joint distribution of the random field
FWHRF (θ) = m −1 2p∑ i1,...,ir,i′1,...,i ′ r=1 wi1 . . . wirJi1,...,ir,i′1,...,i′rwi′1 . . . wi′r , (64)
with its derivatives. In the course of proving Theorem 5, we already have shown that the function value is gamma distributed (see equation 57). Here, we explicitly calculate the distribution of the Hessian when given the function value and that the covariant gradient is zero, and also calculate the distribution of the gradient given the function value. We will heavily lean on the rotational invariance property of the distribution discussed in the proof of Theorem 5, with n once again the fixed point with all θi = 0. Note that for the given embedding of the hypertorus into R2 p
, the Christoffel symbols are zero (i.e. we are considering the Euclidean hypertorus) and thus for the most part we can ignore the distinction between covariant and normal derivatives. Here, we choose local coordinates θ such that:
wi = ( cos (θi) sin (θi) ) . (65)
Perhaps surprisingly, we will find that conditioned on being at a critical point at a specified energy, the Hessian takes the simple form of a normalized and shifted Wishart matrix summed with a normalized GOE matrix. The gradient conditioned on the function value is similarly simple, given by independent Gaussian variables. Lemma 4 (Hessian and gradient distributions). The scaled Hessian m∂i∂jFWHRF (w) conditioned on FWHRF (w) = x and ∂kFWHRF (w) = 0 is distributed as
mC̃ij (x) = −2rmxδij + rWij + r √ 2mxNij , (66)
where W ∼ Wp (2m, Ip) and N ∼ GOEp are independent. Furthermore, the scaled gradient m∂kFWHRF (w) conditioned on FWHRF (w) = x is distributed as
mG̃k (x) = √ 2mrxNk, (67)
where Nk are i.i.d. standard normally distributed random variables independent from all Wij and Nij .
Proof. Without loss of generality we take w = n. Let i ∈ {1, 2}×p be the vector with the ith component equal to 2 and all others equal to 1, (i, j) similar with both the ith and jth component, and b the vector with all components equal to 1. Taking derivatives explicitly yields
m∂iFWHRF (n) = 2 Re { J(i,b,...,b),(b,...,b) } + . . .+ 2 Re { J(b,...,i),(b,...,b) } (68)
and m∂i∂jFWHRF (n) = −2rδijJ(b,...,b),(b,...,b) + 2 Re { J(i,b,...,b),(j,b,...,b) } + 2 Re { J(i,b,...,b),(b,j,...,b) } + . . .+ 2 Re { J(b,...,b,i),(b,...,b,j)
} + 2 Re { J((i,j),b,...,b),(b,...,b) } + 2 Re { J(i,j,...,b),(b,...,b) } + . . .+ 2 Re { J(b,...,b,(i,j)),(b,...,b) } .
(69) As J is a Wishart matrix with identity scale matrix, it can be written as X ·X† for X a 2q ×m matrix with i.i.d. standard complex normal entries. By performing an LQ decomposition of X , one can then by inspection determine the distributions of the entries of J (Srivastava, 2003; Yu et al., 2014). For ease of notation, we let τ : {1, 2}×q → {1, . . . , 2q} be a mapping between representations of the indices of J , with the convention τ ((b, . . . , b)) = 1. We then find (taking τ ((i, . . . , 2)) < τ ((j, . . . , 2)) WLOG) that
2J(b,...,b),(b,...,b) = 2mFWHRF (n) , (70) 2 Re { J(i,...,b),(b,...,b) } = √ 2mFWHRF (n)M(b,...,b),(i,...,b), (71)
2 Re { J(i,j,...,b),(b,...,b) } = √ 2mFWHRF (n)M(b,...,b),(i,j,...,b); (72)
and, for τ ((i, . . . , b)) ≤ m,
2 Re { J(i,...,b),(j,...,b) } = √ 2Γ(i,...,b)M(i,...,b),(j,...,b)
+ τ((i,...,b))−1∑ µ=1 Mτ−1(µ),(i,...,b)Mτ−1(µ),(j,...,b) + τ((i,...,b))−1∑ µ=1 M̃τ−1(µ),(i,...,b)M̃τ−1(µ),(j,...,b)
(73) and otherwise
2 Re { J(i,...,b),(j,...,b) } = m∑ µ=1 Mτ−1(µ),(i,...,b)Mτ−1(µ),(j,...,b)
+ m∑ µ=1 M̃τ−1(µ),(i,...,b)M̃τ−1(µ),(j,...,b).
(74)
Here, M and M̃ are symmetric with off-diagonal entries i.i.d. drawn from the standard normal distribution, and Γτ−1(µ) ≡ M2τ−1(µ),τ−1(µ) has entries i.i.d. drawn from Γ (m− µ+ 1, 1). Note that each √ 2Γ is chi-square distributed with 2 (m− µ+ 1) degrees of freedom; therefore, equa-
tion 73 and equation 74 can be considered as elements of a real Wishart matrix W̃ with 2m degrees
of freedom. Also, note that equation 73 and equation 74 are independent of ∂kFWHRF (n) when conditioned on FWHRF (n) ≡ x = 0. If x 6= 0, the condition ∂kFWHRF (n) = 0 is equivalent to taking each sum over the elements ofM from µ = 2 instead of µ = 1, which is equivalent to taking the convention τ ((b, . . . , b)) = 2q and shifting the indices of M and M̃ . Therefore, the (scaled) Hessian conditioned on FWHRF (n) = x and ∂kFWHRF (n) = 0 is distributed as
mC̃ij (x) = −2rmxδij + ( O · W̃ ·Oᵀ ) ij + r √ 2mxNij ; (75)
here, N ∼ GOEp (with the convention that diagonal entries ∼ N (0, 2) and off-diagonal entries ∼ N (0, 1)), and O is a matrix such that Oiµ = 1 if and only if τ−1 (µ) is of the form (b, . . . , i, . . . , b), and is otherwise equal to 0. The invariance of the Wishart distribution under orthogonal transformations and partitioning (Eaton, 2007; Srivastava, 2003) leads to the final result.
C.2 THE EXACT DISTRIBUTION OF CRITICAL POINTS
Given the joint distribution of FWHRF, its gradient, and its Hessian, we are now equipped to calculate the expected number of critical points of a given index k using the Kac–Rice formula (Lemma 3). Theorem 6 (Distribution of critical points in WHRFs). Let
µC(x) = 1
p p∑ i=1 δ ( λ C(x) i ) (76)
be the empirical spectral measure of the random matrix
C (x) = r
m
( W + √ 2mxN ) , (77)
where W ∼ Wp (2m, Ip) and N ∼ GOEp are independent and λCi (x) is the ith smallest eigenvalue of C (x). Then, the distribution of the expected number of critical points of index k at an energy E > 0 of FWHRF is given by E [Crtk (E)]
= (π r ) p 2 Γ (m) −1 m(1+γ)mEC(E) [ ep ∫ ln(|λ−2rE|)dµC(E)1 { λ C(E) k+1 ≥ 2rE }] E(1−γ)m−1e−mE ,
(78) where
γ = p
2m . (79)
Proof. As discussed in Appendix B, the assumptions of the Kac–Rice formula (i.e. Lemma 3) are satisfied. Furthermore, due to the invariance of the Wishart distribution with respect to rotations on the hypertorus (Eaton, 2007; Srivastava, 2003), we can integrate out the volume element independently; the volume of ( S1 )×p
is ∫ (S1)×p dw = (2π) p . (80)
Additionally, we have from Lemma 4 that the probability density of the gradient vector being zero at any w conditioned on HWHRF (w) = x is
pw (∇HWHRF (w) = 0 | HWHRF (w) = x) = ( 4πrx
m
)− p2 . (81)
Taking the expectation over x via equation 57 and using the Hessian distrribution from Lemma 4, we have from Lemma 1 that
E [Crtk (B = (0, E))] = (π r ) p 2 Γ (m) −1 m(1+γ)m
× E∫
0
EC(x) [ ep ∫ ln(|λ−2rx|)dµC(x)1 { λ C(x) k+1 ≥ 2rx }] x(1−γ)m−1e−mx dx .
(82)
Taking the derivative of this cumulative distribution with respect to E yields the final result.

heading: D LOGARITHMIC ASYMPTOTICS VIA FREE PROBABILITY THEORY
text: Though equation 78 is exact, it is difficult to use in practice. Luckily, we are able to use a surprising fact about the eigenvalue distributions of Wishart and GOE matrices; asymptotically, the empirical spectral distributions of these matrices weakly converge to fixed distributions. Concretely, in the limit p → ∞ where γ = p2m is held constant, the eigenvalue distribution of W /2m where W ∼ Wp (2m, Ip) weakly converges to the Marchenko–Pastur distribution (Marčenko & Pastur, 1967):
dµM.P. = ( 1− γ−1 ) 1 {γ > 1} δ (λ) dλ+ 1
2πγλ
√( (1 + √ γ) 2 − λ )( λ− (1−√γ)2 ) dλ . (83)
Similarly, the eigenvalue distribution ofN/ √ p whereN ∼ GOEp weakly converges to the Wigner semicircle distribution (Wigner, 1958):
dµs.c. = 1
2π
√ 4− λ2 dλ . (84)
Furthermore, by using free probability theory one can find the asymptotic distribution of eigenvalues for a weighted sum of these matrices, given their eigenbases are in “generic position” with respect to each other. We now give a brief review of free probability theory—at least in the context of random matrix theory—here. Later, we will also briefly review large deviations theory, which we use to bound the probability of large deviations from the weak convergence of the eigenvalue distributions of Wishart and GOE matrices. Note that, as we are unable to control large deviations in Theorem 4, in principle large deviations in the weak convergence of VQAs to WHRFs could dominate the large deviations in WHRFs; however, as discussed in Appendix A.3, this provably does not occur at shallow enough depths with respect to n, and there are reasons to believe it does not occur even at large depths (which we additionally give numerical evidence for in Sec. 5).
We begin by reviewing the techniques in free probability theory and large deviations theory that we use in studying the asymptotic behavior of equation 78.
D.1 FREE PROBABILITY THEORY
Free probability theory is the study of noncommutative random variables. Specializing to random matrix theory on N ×N matrices, we define the unital linear functional
φ (X) ≡ 1 N E [tr (X)] (85)
as the free analogue of the expectation. Note that the eigenvalues of a matrix A are completely constrained by the trace of powers Ak—therefore, one can study the average distribution of the eigenvalues of a random matrix A via the moments φ ( Ak ) . Free independence (or freeness) is a generalization of the notion of independence in commutative probability theory to free probability theory. In the context of random matrix theory, two N ×N random matrices A and B are said to be freely independent if the mixed moments are identically zero; that is,
φ ((Am1 − φ (Am1)) (Bn1 − φ (Bn1)) . . . (Amk − φ (Amk)) (Bnk − φ (Bnk))) = 0 (86) for all ni,mi ∈ N. Roughly, the free independence of two random matrices means that their eigenbases are in “generic position” from one another.
Taking the analogy with commutative probability theory further, the analogue of the momentgenerating function associated with the distribution of a random variable is the Stieltjes transform of the measure µ:
Gµ (z) =
∫ dµ (t)
z − t , (87)
which can be inverted via the Stieltjes inversion formula:
dµ (t) = − 1 π lim →0+ Im {Gµ (t+ i )} dt . (88)
Similarly, the free analogue of the cumulant-generating function is the R-transform, which can be defined via the Stieltjes transform as the solution to the implicit equation:
Rµ (Gµ (z)) + 1
Gµ (z) = z. (89)
The R-transform is important in that, if two random variables A and B are freely independent with probability measures µA and µB respectively, the probability measure µA+B of A+B satisfies
RµA+B = RµA +RµB . (90) This can be interpreted as the free analog of the additivity of cumulants for commutative random variables. The probability measure µA+B is called the free convolution of µA and µB , and is denoted using the notation µA+B = µA µB . (91) Thus, given the probability distributions of two free random variablesA andB, there is a prescription for determining the probability distribution of their sum by taking their free convolution, just as the convolution in commutative probability theory describes the distribution of the sum of random variables.
D.2 LARGE DEVIATIONS THEORY
In order to bound the probability of large deviations from the weak convergence of the eigenvalue distribution of C to its asymptotic limit we will use results from large deviations theory, which we briefly review here. A sequence of measures {µn} is said to satisfy a large deviation principle in the limit n → ∞ with speed s (n) and lower semicontinuous rate function I with codomain [0,∞] if and only if (Dembo & Zeitouni, 2010)
− inf x∈Γ◦ I (x) ≤ lim inf n→∞
1
s (n) ln (µn (Γ )) ≤ lim sup
n→∞
1
s (n) ln (µn (Γ )) ≤ − inf x∈Γ I (x) (92)
for all Borel measurable sets Γ that all µn are defined on. Here, Γ denotes the closure of Γ and Γ ◦ the interior of Γ . The rate function I is said to be good if all level sets of I are compact. Large deviations theory will be useful for us to bound the probabilities of large deviations of the empirical spectral distribution of µC(x) as p→∞, and show that they do not contribute to leading order in the (logarithmic) asymptotic distribution of critical points. We do this using Varadhan’s lemma, which we state now. Lemma 5 (Varadhan’s lemma (Dembo & Zeitouni, 2010)). Suppose {µn} satisfies a large deviation principle with speed s (n) and good rate function I and let φ be a real-valued continuous function. Further assume either the tail condition
lim M→∞ lim sup n→∞
1 s (n) ln ( EXn∼µn [ es(n)φ(Xn)1 {φ (Xn) ≥M} ]) = −∞, (93)
or the moment condition for some γ > 1
lim sup n→∞
1 s (n) ln ( EXn∼µn [ eγs(n)φ(Xn) ]) <∞. (94)
Then,
lim n→∞
1 s (n) ln ( EXn∼µn [ es(n)φ(Xn) ]) = sup x (φ (x)− I (x)) . (95)
D.3 LOGARITHMIC ASYMPTOTICS OF THE DISTRIBUTION OF CRITICAL POINTS
Equipped with these mathematical tools, we prove our first result on the asymptotic behavior of µC(x), which is present in the expectation of equation 78. Lemma 6 (Asymptotic behavior of µC(x)). Define G∗x (z) as the implicit solution of the equation
8r3γ2xG∗x (z) 3 − 2rγ (z + 2rx)G∗x (z) 2 + (z − 2r (1− γ))G∗x (z)− 1 = 0 (96)
with the smallest imaginary part. Define
dµ∗x ≡ − 1
π Im {G∗x} dλ . (97)
Let p,m → ∞ as γ = p2m is held constant. Then, the empirical spectral measure µC(x) satisfies a large deviation principle as p→∞ with speed p2 with good rate function uniquely minimized at µ∗x with a value of 0.
Proof. The empirical spectral measure of the random matrixN/ √ p satisfies a large deviation principle at a scale p2, with good rate function minimized by Wigner’s semicircle law (Arous & Guionnet, 1997). Similarly, the empirical spectral measure of the random matrix W /2m satisfies a large deviation principle at a scale p2, with good rate function minimized by the Marchenko–Pastur distribution (Hiai & Petz, 1998). As the R-transform of the empirical spectral distribution of A satisfies the scaling property RaA (z) = aRA (az) , (98) the R-transform of the empirical spectral distribution of the weighted GOE term ofC is of the form
RGOE (z) = 4r2γxz (99)
and the R-transform of the weighted Wishart term is of the form
RWishart (z) = 2r
1− 2rγz . (100)
By the asymptotic freeness of independent GOE and Wishart matrices (Hiai & Petz, 2006), µC(x) converges weakly to the fixed measure µ∗ with R-transform
Rx (z) = RWishart (z) +RGOE (z) . (101)
Equation 96 and equation 97 now follow from inverting the R-transform Rx via equation 89 and equation 88, respectively.
We now consider large deviations in the weak convergence µC(x) µ∗x. Conditioning on the empirical spectral distribution of W /2m and using the “strongest growth wins” principle (Dembo & Zeitouni, 2010), we have that µC(x) satisfies a large deviation principle with speed p2 with rate function given by
I (µ) = inf µW/2m
( JµW/2m (µ) +K ( µW /2m )) ; (102)
here, K is the rate function governing convergence of the empirical spectral distribution of the Wishart ensemble (Hiai & Petz, 1998) and JµW/2m is the rate function governing convergence of the empirical spectral distribution of a fixed matrix with asymptotic eigenvalue distribution µW /2m summed with a GOE matrix (Guionnet & Zeitouni, 2002). This sum is obviously uniquely minimized by µ = µ∗, when I (µ∗) = 0.
Now, we examine the asymptotic behavior of the smallest eigenvalue λC(x)1 of C (x). Unlike the empirical spectral measure µC(x) which satisfies a large deviation principle at a speed p2, we will see that this eigenvalue satisfies a large deviation principle at a speed p, with deviations at this speed to the left of the asymptotic value λ∗x,1.
Lemma 7 (Asymptotic behavior of λC(x)1 ). Let λ∗x,1 be the infimum of the support of µ∗x as defined in equation 97. Then, the smallest eigenvalue λC(x)1 of C (x) satisfies a large deviation principle with speed p with good rate function that is infinite at y > λ∗x,1 and is uniquely minimized at y = λ ∗ x,1 with a value of 0.
Proof. The limiting smallest eigenvalues λW1 , λ N 1 of W /2m and N/
√ p both satisfy large devia-
tion principles with speed p that are infinite for λ1 in the bulk of their respective limiting empirical spectral distributions (Johansson, 2000; Arous et al., 2001). As in the proof of Lemma 6, we condition on large deviations of these eigenvalues (Dembo & Zeitouni, 2010) and therefore have that the rate function governing λC(x)1 is
I (y) = inf λW1 ,λ N 1
( JλW1 ,λN1 (y) +K ( λW1 ) + L ( λN1 )) ; (103)
here, K is the rate function governing the convergence of λW1 , L that of λ N 1 , and J that of the smallest eigenvalue C (x) conditioned on the eigenvalue distributions of W and N . Using known results on the large deviations of the smallest eigenvalue of the sum of two matrices with fixed eigenvalues (i.e. JλW1 ,λN1 ) (Guionnet & Maı̈da, 2020), we see that I (y) is infinite for y > λ ∗ x,1 and is uniquely minimized at y = λ∗x,1 with a value of 0.
Using Lemmas 6 and 7, we can prove the following logarithmic asymptotics on the expectation term in equation 78. We will find that neither the large deviations in the convergence µC(x) or λ C(x) 1 will contribute to leading order in the logarithmic asymptotics of Crtk (E), as at a speed p the only large deviations are λC(x)1 ≤ λ∗x,1 which are dominated by λ C(x) 1 = λ ∗ x,1 in the expectation. Lemma 8 (Logarithmic asymptotics of the determinant). Let dµ∗E be the spectral measure given in equation 97, with λ∗E,1 the infimum of its support. Let p,m 1 with p 2m = γ = O (1). Then,
1 p ln ( EC(E) [ ep ∫ ln(|λ−2rE|)dµC(E)1 { λ C(E) k+1 ≥ 2rE }]) = ∫ ln ( 1 { λ∗E,1 ≥ 2rE } |λ− 2rE| ) dµ∗E + o (1) .
(104)
Proof. As µC(E) satisfies a large deviation principle with speed p2 with rate function minimized at µ∗E by Lemma 6, and as 1 { λ C(E) 1 ≥ 2rE } ≤ 1, we have that the tail condition of Varadhan’s
lemma at speed p is satisfied (Dembo & Zeitouni, 2010) and therefore
lim p→∞
1 p ln ( EC(E) [ ep ∫ ln(|λ−2rE|)dµC(E)1 { λ C(E) 1 ≥ 2rE }]) = sup λ∈R (∫ ln (1 {λ ≥ 2rE} |λ− 2rE|) dµ∗E − I (λ) ) .
(105)
Here, I is as in equation 103. The supremum over λ is obviously achieved when λ = λ∗E,1 by the properties of I discussed in Lemma 7, giving the leading order term in equation 104. The result being exact in the p→∞ limit gives the subleading o (1).
Using Lemma 8, we can therefore finally calculate the logarithmic asymptotic distribution of local minima of a WHRF.
Theorem 7 (Logarithmic asymptotics of the local minima distribution). Let dµ∗E be the spectral measure given in equation 97, with λ∗E,1 the infimum of its support. Let p,m 1 with p 2m = γ = O (1). Then, the expected distribution of local minima of FWHRF at a fixed energy E > 0 is given by
1 p ln (E [Crt0 (E)]) = 1 2 ln
( πq
2γ
) + 1
2γ (1− E) + 1 2
( γ−1 − 1 ) ln (E)
+ ∫ ln (∣∣∣∣λr − 2E ∣∣∣∣1{λ∗E,1r ≥ 2E }) dµ∗E + o (1) .
(106)
Proof. The result follows directly from applying Lemma 8 to Theorem 6.
Note that, though we only prove the asymptotic distribution of local minima in Theorem 7, we expect similar theorems to also hold for critical points of constant index k (taking λ∗E,1 7→ λ∗E,k in the integrand). The only difference in the derivation is the exact form of the large deviations of the kth smallest eigenvalue of C (x). This is similar to the case in Gaussian hyperspherical random fields (Auffinger et al., 2013).

heading: E DETAILS OF THE NUMERICAL SIMULATIONS
text: We now give further details on the numerical simulations performed in Sec. 5. We performed all simulations via Qiskit (Abraham et al., 2019), and used standard gradient descent (via the method of finite differences) to optimize the VQA loss function
F (θ) = 〈θ|HT ,U |θ〉 (107)
until convergence. HT,U is the 1D n site spinless Fermi–Hubbard Hamiltonian (Negele & Orland, 1998)
HT ,U = − n−1∑ i=1 Ti ( c†i ci+1 + c † i+1ci ) + n−1∑ i=1 Uic † i cic † i+1ci+1, (108)
where c is the fermionic annihilation operator. Ti and Ui are i.i.d. normally distributed in order to break translational invariance. In our simulations, these random variables were centered at T = 1 and U = 2, respectively, and each had a variance of 10−2.
Our implementation of gradient descent used a learning rate of 0.05 and a momentum of 0.9, and halted when either the function value improved by no more than 10−5 or after 106 iterations, whichever came first. We initialized each instance at a uniformly random point in parameter space, with each parameter initialized within [−2π, 2π]. To estimate the empirical distribution of local minima for the studied instances of the varitional quantum eigensolver (VQE) (Peruzzo et al., 2014), we repeated this procedure 52 times, using a new ansatz and uniformly random starting point for each training instance. We also verified numerically that m as defined in equation 26 is at least on the order of 2n for H1,2, though we directly used equation 26 when computing γ. In all plotted instances, we normalize the energy scale by a factor of cVQA, where cVQA = λ− λ1; (109) this is just the overall factor of equation 8. These units are such that the mean eigenvalue of H − λ1 in the subspace of interest is at E = 1.
In Sec. 5.2, we tested our analytic results against a Hamiltonian informed ansatz. Specifically, we used the Hamiltonian variational ansatz (HVA) (Wecker et al., 2015). For the Fermi–Hubbard Hamiltonian of equation 108, each HVA layer is of the form
UT ,Ui (θi) = e −iθi,2HT,U,odde−iθi,2HT,U,evene−iθi,1HT,U,Coulomb . (110)
Here, HT ,U ,Coulomb is composed of the terms proportional to Ui in HT ,U , HT ,U ,even the hopping terms on even links, and HT ,U ,odd the hopping terms on odd links. We took the starting state |ψ0〉 to be the computational basis state |1〉 on the first n2 qubits and |0〉 on the other n 2 qubits. To observe the effects of scaling the number of independent parameters p, we overparameterize our ansatz at a fixed overall depth by fixing the total number of ansatz layers UT ,Ui to be 6, but introduce extra parameters to govern each evolution. For instance, for a multiplicative factor f = 2, we double the number of parameters by splitting into a sum of two terms
HT ,U ,Coulomb = H (1) T ,U ,Coulomb +H (2) T ,U ,Coulomb, (111)
HT ,U ,even = H (1) T ,U ,even +H (2) T ,U ,even, (112)
HT ,U ,odd = H (1) T ,U ,odd +H (2) T ,U ,odd, (113)
and parameterize the evolution under each term separately. For f = 1, this ansatz preserves the fermion number of the initial state; thus, for these simulations we calculate m in this n2 -fermion subspace. For large f , this parameterization breaks the fermion number conservation of the ansatz, but still preserves the parity of the fermion number. In practice, then, the γ we compute should be considered an upper bound on the true γ, strengthening our empirical results.

Review summaries:
The work attempts to connect the optimization of quantum generative models with random fields on manifolds. Some theoretical and numerical analyses are provided.   The major idea is to apply random matrix theory for variational quantum circuits. When the authors provide some numerical experiments, overall the validation is not very strong and could be improved.   The connection between applying random fields and quantum circuit training is not very direct based on the current writing. Some related references in the field are missing. (e.g., [1,2]).  ***  **Reference**  1. Random matrix spectral form factor of dual-unitary quantum circuits." Communications in Mathematical Physics 387.1 (2021): 597-620.  2. Abrupt transitions in variational quantum circuit training." Physical Review A 103.3 (2021): 032607.

The paper analyzes a general class of variational quantum circuits using random matrix theory, and gives a theoretical characterization of the distribution of high-quality local minima of the associated Hamiltonian loss function. A type of phase transition is predicted, where the model exhibits an abundance of high-quality local minima only at very large parameter counts, and a wide gap between global and typical local minima at more typical parameter counts. These predictions are compared with small numerical experiments, as well as previous works on this topic.

This paper obtains rather impressive results about the problem of estimating a quantum ground state energy using a variational method which could be implemented on a quantum computer.  Given a Hamiltonian H, the method generates quantum states by applying a succession of n-qubit Pauli operators to a random initial state (both chosen without knowledge of H), and measure the expectation of H. Previous work suggested that there is a critical (and exponentially large) number of parameters required to get a good estimate.  Here a theorem is formulated and proven that there is a phase transition of this type, with good estimates only above a threshold in the number of parameters, and with an explicit expression for the threshold.  The proof uses a combination of physics methods (series expansions), random matrix theory and the Kac-Rice method to count critical points of a random function.