Universal Effectiveness of High-Depth Circuits in Variational Eigenproblems

We explore the effectiveness of variational quantum circuits in simulating the ground states of quantum many-body Hamiltonians. We show that generic high-depth circuits, performing a sequence of layer unitaries of the same form, can accurately approximate the desired states. We demonstrate their universal success by using two Hamiltonian systems with very different properties: the transverse field Ising model and the Sachdev-Ye-Kitaev model. The energy landscape of the high-depth circuits has a proper structure for the gradient-based optimization, i.e. the presence of local extrema -- near any random initial points -- reaching the ground level energy. We further test the circuit's capability of replicating random quantum states by minimizing the Euclidean distance.


I. INTRODUCTION
Variational Quantum Eigensolver (VQE) [1,2] is one of the most promising hybrid quantum-classical (HQC) algorithms, which may offer a precise approximation of the ground state of quantum systems.It is based on the iterative application of the following three steps: state preparation, measurement and optimization.Let us briefly describe each step.First, the preparation of a trial state |ψ(θ) is carried out by successive application of unitary quantum gates that depend on variational parameters θ.Second, the measurement step estimates the trial state mean energy, by taking the expectation value of Hamiltonian H of the target system over the trial state.Third, the optimization step adjusts the variational parameters θ of the trial wavefunction to minimize the mean state energy E(θ) by applying a classical optimization algorithm.See [3] and references therein for more details.After a sufficient number of iterations, the variational state |ψ(θ * ) at a convergence point θ = θ * is expected to reproduce well the ground state of the target Hamiltonian H, under the assumption that the state Ansatz |ψ(θ) is parametrically expressible and well-trainable under the gradient based optimization.It is then a crucial question to find such an Ansatz.Ideally, we are looking for a universally effective Ansatz capable of solving the VQE problem associated with an arbitrary target Hamiltonian.Note that it is generically a very challenging task: unless the target Hamiltonian is local, the ground states of non-local interacting systems are relatively close to typical quantum states that constitute most of the Hilbert space.Universality is thus equivalent to demand a circuit Ansatz to approximate any random state |φ at certain values of the circuit parameters.
The existence of a parameter set ϕ, such that |ψ(ϕ) |φ , is therefore a necessary condition for the effectiveness of the variational circuit.On top of that, we must concern if the gradient descent optimization can actually reach the parameters ϕ from randomly initialized ones.Suppose that we use a layered quantum circuit composed of repetitive application of variational layers with the same architecture.Since increasing the number of layers extends the dimension of the parameter space, it can affect the circuit's expressibility only positively.However, there has been a reported tension between increasing depth of layers and trainability of the circuit via minimizing the mean energy, (1), known as the barren plateau phenomenon [4].
When the layered circuit reaches a certain depth such that it evolves to an approximate 2-design, a numerical experiment [4] has shown the exponential decay of the variance of energy gradients ∇ θ E(θ) with respect to the number n of qubits -for random circuit states obtained by uniformly sampling from the parameter space.Combined with another fact that the random energy gradient vanishes on average, Chebyshev's inequality implies that the initial energy gradient can exponentially rarely deviate from zero.Such diminishing gradients hinder the beginning of efficient energy minimization, possibly causing the variational state to be stuck on non-optimal plateaus.
There have been several proposals to overcome the vanishing gradient problem and optimize the variational circuit.The most obvious approach is to incorporate some physics information on the target Hamiltonian [5][6][7][8], e.g. the ground state symmetry, for designing a less generic problem-tailored Ansatz.As an alternative direction towards the universal applicability of the variational algorithms, novel initialization [9,10], architecture [11,12], and optimization [13][14][15] of generic purpose circuits have been developed to enhance the circuit performance.
The main goal of this paper is to demonstrate the inherent capability of variational circuits as universal and accurate eigensolvers for generic Hamiltonians -if they can even approximate close-to-random states.To this end, we will simply put aside the barren plateau problem by sufficiently increasing the classical computation capacity.Specifically, we will consider the high-depth regime of the layered quantum circuit, thus building an exponentially high-dimensional parameter space.In this regime, the variational circuit has sufficiently many layers and consequently many parameters.Hence, the norm of the gradient vector ∇ θ E(θ) can still grow to a finite size, capable of moving the circuit parameter θ from initial points, even though the magnitudes of the individual components |∂ i E(θ)| are exponentially suppressed for the number n of qubits.We will illustrate the effectiveness of the high-depth circuits by solving two concrete VQE problems for the following quantum many-body Hamiltonians: the 1d Ising model in a constant transverse magnetic field, which is a prototypical model of locally interacting spin-chain systems, and the Sachdev-Ye-Kitaev (SYK) model, which is a strongly interacting quantum mechanical system of Majorana fermions [16][17][18].Despite the striking contrast in these two Hamiltonians' ground state properties, we will find that the high-depth circuit achieves a very high fidelity in approximating both ground states.Moreover, we will see that the high-depth circuits can narrow the Euclidean distance from random states up to arbitrary precision, indicating outstanding expressibility and trainability that stems from the highdimensionality of θ.
The remarkable efficiency of the gradient descent optimization in the over-parameterized regime was also reported for approximating random unitary matrices [19] and ground states [8] with certain variational Ansatzes.More generally, over-parameterization is an active research topic that underlies the success of deep learning.It makes large neural networks capable to reach a global minimum during the optimization process, despite the non-convexity of the energy landscape [20,21].While searching the ground state via the energy minimization, we will observe some phenomena similar to what happens during the training of the neural networks with large parameter spaces, summarized as follows: First, the energy landscape looks fairly simple in the local vicinity of randomly chosen points [15].Generically, an initial point is already confined in a certain basin of attraction, such that an emanating optimization trajectory can quickly arrive at a nearby local extremum.Especially if the circuit has enough layers to become an approximate 2-design, almost all uniformly sampled initial states end up with rather homogeneous energy levels.
Second, for high-depth circuits, all the local extrema in the energy landscape, reachable by the VQE optimization from randomly initialized points, are substantially close to the exact ground energy, i.e., the value of the global minimum.These minima are not isolated individ-ual points but develop multiple flat directions.It explains the robust success of the high-depth circuits in solving the VQE problems.
The rest of this paper is organized as follows.Section II introduces the architecture of the variational circuits used throughout this paper.Section III first reviews the occurrence of the barren plateau phenomenon, then examines the optimization of the variational circuits using the VQE example of the 1d Ising model coupled to a uniform transverse magnetic field.Section IV studies the same VQE problem for the SYK model, thus showing the universal effectiveness of the high-depth circuit.Section V addresses the ability of the variational circuits to reproduce typical states by minimizing the Euclidean distance.In particular, it shows that the high-depth Ansatz is highly expressible and trainable, being able to reach any random state.Finally, Section VI concludes with discussions.

II. CIRCUIT ANSATZES
Our focus in this work is to demonstrate the efficiency of high-depth layered circuits in a typical VQE problem, i.e. to approximate the ground state of a given Hamiltonian.To this purpose, here we specify the architecture of the variational circuit used in our numerical experiments.Our circuit state |ψ(θ) is composed of L unitary layers acting sequentially on the initial product state |0 , i.e., where each layer U i (θ i ) has single-qubit y-rotation gates, parameterized by n periodic variables θ i , and controlledz gates operating on all pairs of n qubits.More precisely, once the single-qubit RY gates have acted upon all individual qubits, the entangling CZ gates acting on a-th controlling and b-th targeting qubits are arranged for every integer pair (a, b) satisfying 1 ≤ a < b ≤ n.We have drawn the circuit state (2) with n = 4 qubits in Figure 1.
Note that we have deliberately chosen an unbiased circuit architecture, instead of embedding known properties of the ground states that are under search.Nonetheless, as we will see in Sections III, IV and V, the above circuit can solve the ground states of distinct Hamiltonians by minimizing E(θ) and accurately approximate random quantum states by minimizing the Euclidean distance, as long as the number L of layers is sufficiently large.

III. LOOKING INTO VQE TRAJECTORIES
This section is devoted to the detailed investigation of the VQE optimization procedure, with a particular focus on the effectiveness of the high-depth circuit Ansatz (2).Our exploration will be based on a concrete Hamiltonian system, commonly used in measuring the performance of h u 7 Q A 3 p E T 8 6 9 8 + y 8 O K / z 1 p S z m D l A S 3 D e P g H X z J q L < / l a t e x i t > U 3 (✓ 3 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " U 8 P y / Z 9 1 D T d m c x l T F o 7 q T K r T 9 e r F e 5 6 0 p a z F z g J Z g v X 0 C 3 j y a j w = = < / l a t e x i t > |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 r t M C d F m l j j 1 H variational circuits, i.e. the 1d Ising model in a transverse and uniform magnetic field [8,22].
The 1d transverse field Ising model is defined over a spin lattice of length n, consisting of the spin-spin coupling between nearest neighbors in the z direction, as well as the spin interaction with a background uniform magnetic field along the transverse x direction.Assuming periodic boundary conditions, σ z n+1 ≡ σ z 1 , the ferromagnetic Ising Hamiltonian reads where σ x,y,z i is the Pauli operator acting on the i'th spin, and g denotes the strength of the uniform magnetic field.
The physics of this model has been well-studied [23].When the lattice size scales up to infinity, n → ∞, the system undergoes a quantum phase transition at |g| = 1 between the ordered (|g| < 1) and disordered (|g| > 1) phases.The former phase has the spin-flip Z 2 symmetry that connects the two opposite ferromagnetic ground states, while the latter phase has a unique ground state, with all the spins aligned along the x direction.
We will apply the VQE algorithm to the finite n Ising system, which exhibits some differences from the thermodynamic limit.Specifically, the Z 2 degeneracy in the 0 < |g| < 1 phase is broken at finite n.Our target state will be always non-degenerate and gapped at |g| = 0.All concrete calculations presented in this section have been done for g = 2.

Barren Plateaus and Classical Resolution
It was argued in [4] that optimization of the quantum variational circuits under the random parameter initialization comes with an inherent difficulty, given by the problem of vanishing gradients.For the circuit unitary ensemble that is quantum 2-design, random initial parameters are typically located on a plateau in the energy landscape.It means that the gradient-based optimization cannot even roll out.The odds of having non-vanishing gradients at random points decays exponentially with the system size n, impeding a large-scale application of the variational circuits.Its consequence would be even more detrimental on actual quantum devices, where the sampling noise must be taken into account [4,14].Given its prominence, we begin by reviewing the vanishing gradient problem in the VQE setting.
Consider the ensemble of the energy gradients over the parameter space of the variational circuit in Figure 1.For analyzing the k'th component, ∂ k E(θ), that belongs to the 'th variational layer, it is convenient to group the L layer unitaries (2) into the following two blocks: where As the partial derivative ∂ k acts only on U (θ ), the variance of the k'th gradient component, Var where V k denotes the Pauli operator σ y acting on the k'th spin variable, such that Tr(V k ) = 0 and Tr(V 2 k ) = 2 n .If we assume the quantum 2-design property of U + (θ + ) and/or U − (θ − ), the above integral ( 6) can be replaced with the unitary matrix integral.In that case, the matrix integral can be handled exactly and simplified to 2Tr(H 2 ) 2 2n .
Such simplification (7) happens when both U ± (θ ± ) are 2designs.Instead, if the 2-design condition does not hold for either of U ± (θ ± ), we have the following expression: where ρ = |0 0|, if not U + (θ + ) but only U − (θ − ) is a 2-design.Asymptotically in the system size n, the above expressions ( 7)-( 9) are all bounded as where we find the upper bound by expanding the commutator inside the integral of ( 8) and ( 9), then applying the following trace inequality that holds for two Hermitian matrices A, B [24]: To determine the scaling behavior of the upper bound (10) with respect to the system size n, we examine how Tr(H 2 ) scales with n by extrapolating the values obtained from exact diagonalization of the Hamiltonian (3) up to n ≤ 17.The result is presented in Figure 2. We see that the term Tr(H 2 ) scales as 2 n , such that the exponential factor in the denominator of the upper bound (10) cannot fully be balanced by Tr(H 2 ).Then, inserting the upper bound (10) into Chebyshev's inequality, which states with ∂ k E(θ) being a random variable of zero mean, the probability of having a non-zero and finite derivative is exponentially suppressed with growing n.Therefore, the large-scale VQE problems are expected to suffer from the vanishing gradients.
We have not yet justified the assumption that either of the unitaries U ± (θ ± ) is quantum 2-design.To see if the problem of exponentially vanishing gradients happens in the circuit unitaries of Figure 1, we numerically estimate the variance of initial gradients by collecting 10 3 random energy gradients of the Ising model.Figure 3a clearly exhibits an exponential decay of the partial derivatives with the increasing number n of qubits, where the shaded region displays component-wise fluctuations of the variance for all 1 ≤ k ≤ nL.It shows that the energy landscape for the circuit in Figure 1 indeed contain the barren plateaus, which hinders the circuit optimization towards the ground state.
On the other hand, Figure 3a shows that the variance, at a fixed value of n, converges to a constant by increasing the number L of layers.The variance is independent of L beyond a transition point L 0 , where the circuit unitaries with L ≥ L 0 evolve to approximate 2-design.Due to the saturation, having exponentially many parameters can compensate for the exponential decay of individual components when calculating the gradient norm ∇ θ E(θ) , which appears in the evolution of the circuit energy under the gradient descent with an infinitesimal rate α: Therefore, the vanishing gradient problem inherent to the quantum Hilbert space can be resolved in a classical way, i.e. by using the exponentially high-dimensional θ-space.
In agreement with the reasoning above, Figure 3b exhibits a small initial drop dominated by the transient decrease of Var θ [∂ k E(θ)] for L ≤ L 0 layers, then a steady increase driven by the linear growth in the number of circuit parameters.One can approximate the asymptotic increase rate of the norm ∇ θ E(θ) as which agrees well with the numerically estimated growth rates between log ∇ θ E(θ) and log L in Figure 3b: We remark that the barren plateau phenomenon concerns only the initial steps of the gradient descent update.In almost all physical systems of relevance, the Hamiltonian spectrum is symmetrically distributed around a certain value, which we canonically set to zero.A random state is therefore in a superposition of multiple Hamiltonian eigenstates, whose mean energy is almost certainly zero.On the other hand, the variational mean energy quickly becomes a negative value after a few steps of the VQE optimization, implying that the circuit state is no longer represented by sample statistics of random states.
Having found that the vanishing gradient problem can be trivially avoided at the cost of having exponentially many parameters, we turn to answer if the variational circuit with a sufficiently many layers can actually solve the VQE problems.We will examine the optimization error, the training curve, and the trajectory in the parameter space at different values of the circuit depth L.

Optimizing the Circuit
We now come back to original task of approximating the ground state of the Ising model, Eq. ( 3).
The VQE optimization results of the circuit states (2), with different numbers L of layers and under the random initialization of the circuit parameters θ, can be summarized to the following two features: 1.For small enough L, the minimized energies E(θ * ) that the circuit states can reach are highly variable.
2. For larger L, the E(θ * ) distribution turns concentrated around a mean value which gets smaller.
We illustrate them with the outcomes of the VQE algorithm for the case of a chain having n = 10 sites and for L ∈ {8, 10, 12, 14, 16, 18, 20} layers.To find the optimal parameters θ * that minimize the energy, we use the Adam optimization algorithm [25], which iteratively updates the variational parameters θ by the exponential moving averages of gradients and their squares.It has a clear advantage in convergence speed, being widely used in a variety of deep learning models.The parameter update rule is decided by a choice of three hyperparameters (α, β 1 , β 2 ).We have collected 35 independent VQE runs for the Ising Hamiltonian (3), whose initial parameters are randomly sampled from the uniform distribution U(0, 2π) ⊗nL , after 500 parameter updates with the hyperparameters (α, β 1 , β 2 ) = (0.05, 0.9, 0.999).
The sample distribution of the final VQE energy E(θ * ) is visualized in Figure 4 for different numbers L of layers.One can characterize it as follows.First, the energy distribution at the local extremum θ * clearly exhibits the widespread spectrum for a shallow depth, e.g., for L = 10.
sometimes achieving a relatively good energy level while the gap always persists.Second, by stacking more layers, i.e., L ≥ 12, deeper than the 2-design transition point detected in Figure 3a, the energy distribution starts to concentrate around a value, which is far from the ground energy E 0 .Third, the average value of E(θ * ) continues to decrease for the growing number L of layers, suggesting that the high-depth variational circuits can possibly simulate the ground state using the VQE optimization.Encouraged by the observed shrinkage of the mean and variance of E(θ * ), we also have done the single VQE tryouts for a broader span of (n, L), as summarized in Figure 5.It is clear that being deeper enhances the variational circuit's capability to replicate the ground state |ψ 0 , reaching the ground state energy E 0 and achieving a high level of fidelity with |ψ 0 .The high-depth circuits can achieve a zero error with the following precision: where ∆E denotes the bandwidth of the target Hamiltonian H, defined as the difference between the largest and smallest eigenvalues of H.We note that such precision can be achieved only with an appropriately chosen learning rate α, in order to avoid too-large parameter updates that prevent the fine-level optimization.Here we refrain from the systematic hyperparameter search, which may be more relevant for the case where the average gap between nearby energy levels shrinks, but simply stick to α = 0.05 (L = 4, 6, 8) and α = 0.01 (L = 10).We have found that (15) can be achieved when the circuit depth L passes through the following threshold value: We also remark that deeper circuits do not necessarily lead to better performance with VQE, as displayed by the gentle ramp after passing the threshold point L v .This can be understood as follows: As the space of variational parameters θ has more dimensions, the basin of attractor to local extrema becomes narrower [11], giving a larger value of the estimated inverse volume [26] where the summation is taken over the top-k eigenvalues {λ i (H)} k i=1 of the Hessian matrix For instance, the positive correlation between L and V −1 k (L) can clearly be identified in Figure 7a, drawn for k = 100.As a result, the VQE trajectory is unlikely to land at an exact extremum but wander around nearby points, whose deviation gets larger as the attractor basin becomes narrower and steeper.
Even with the optimal number of layers, such that the circuit state can reach an exact extremum and accurately represent the ground state |ψ 0 of the Ising Hamiltonian (3), the narrowness of the attractor basin still makes the VQE trajectory somewhat unstoppable, passing through the best point θ * and then hopping around in the local neighborhood.Figure 6a shows that the VQE error E(θ τ ) − E 0 slightly increases on average and mildly fluctuates after achieving the minimum error E(θ * ) − E 0 .This residual error can be reduced by making use of popular optimization tricks, such as early stopping or learning rate scheduling, which causes the VQE optimization to stop at the optimal point.For instance, by introducing the exponential decay of the learning rate α, i.e.
at the optimization step τ with a constant value c = 0.3, we could reduce the late time fluctuations as in Figure 6b.So far, we have discussed some important aspects of the VQE optimization.The main observation is that, when supported by the high-dimensional parameter space, randomly initialized variational circuits can approximate the Ising ground state with remarkably high accuracy.The energy gradient neither vanishes nor randomly fluctuates along the optimization trajectory, making it quickly converge to a local minimum that exhibits a small energy gap from the ground energy E 0 and high fidelity with the exact ground state |ψ 0 .We will explore this efficiency of the high-depth circuit in some details by visualizing the VQE trajectory on the energy landscape.

Visualizing the Trajectory
A key characteristic that contributes to the success of the high-depth circuit is the fact that randomly initialized points are likely to be already confined in the basin of attraction of a good attractor, i.e. a local extremum that is close enough to the ground state energy.We illustrate now this point by looking at the actual optimization trajectories under the energy minimization.
For a given initial point θ 0 and the trajectory T (θ 0 ) = {θ τ } thereafter, we identify the optimal parameter θ * as the point in the trajectory T (θ 0 ) having the minimum energy expectation value, E(θ * ) ≤ E(θ τ ) for any τ ≥ 0, (18) which is the best possible representative point of the local extrema that the trajectory T (θ 0 ) converges around.
The associated basin of attractor can be estimated by calculating the Hessian H ij ≡ ∂ i ∂ j E(θ * ) at the optimum θ * .We are interested in the eigenvalue spectrum, {λ i (H)}, from which we distinguish steep and flat directions and calculate the degree of steepness.Similar to the case of deep neural networks [26][27][28], the local extrema in the VQE energy landscape of highdepth circuits are often interconnected in multiple flat directions, whose corresponding Hessian eigenvalues are zero, looking like valleys rather than isolated singular points.Success of the VQE algorithm is not affected by a specific position in flat directions, but only concerned with if a ball, initially at a considerable height, can roll off in steep directions and reach a sufficiently deep gorge.It motivates us to consider the k-dimensional hypersurface S k (θ * ) along the k steepest directions, i.e., spanned by the top-k Hessian eigenvectors.More precisely, we want to examine if the S k (θ * )-projected Euclidean distance between θ * and θ τ decreases along the trajectory T (θ 0 ) as the optimization step τ progresses, until E(θ τ ) = E(θ * ).This will tell us if the entire trajectory T (θ 0 ) is confined in the k-dimensional basin of attraction, while ignoring the movement in the directions of less or zero attraction.
Figure 7b displays the VQE error E(θ τ ) − E 0 and the projection distance θ τ −θ * S100(θ * ) along the actual optimization trajectories, whose step number τ is indicated by color.We have made the following observations: First, when an appropriate value of k is selected, both the dis-tance θ τ − θ * S k (θ * ) and the loss value E(θ τ ) − E 0 continuously decrease on a macroscopic scale.It exhibits that the trajectory T (θ 0 ) converges without escaping from a specific basin of the attractor that encloses a randomly initialized point θ 0 .Second, for a shallow circuit state, the optimization trajectory often makes a slight detour in some orthogonal directions.In contrast, the steady convergence occurs typically for large L. It implies that the vicinity of any randomly initialized parameters is effectively convex.Finally, we find from Figure 7a that the attractor basin in the direction of S k (θ * ) evolves rapidly steeper and narrower as the depth L increases, thereby causing the rapid convergence and substantial late-time fluctuation around θ * .It is noticeable that the initial convergence follows a milder path than later fluctuations in the L = 64 case, while the convergence in the L = 80 case happens along a much steeper route.Taken together, these observations show the quick convergence of high-depth circuits [8,19] under the VQE algorithm.

IV. SOLVING THE SYK MODEL
Our discussions so far have been based on just a particular Hamiltonian, (3), the Ising model in a transverse uniform magnetic field.Since the variational circuit that we use has no features particularly tailored for the Ising model, we expect it to closely replicate the ground states of other Hamiltonians in the high-depth regime.To check the generality of the prior discussions on the efficiency of the high-depth circuits, we will now solve the VQE problem, defined for another Hamiltonian of very distinct nature: the Sachdev-Ye-Kitaev (SYK) model.

The SYK Model
The SYK model [16][17][18] is built out of 2n Majorana fermions in 1d, i.e. the operators γ i , with i = 1, . . ., 2n, satisfying the following anti-commutation relations where δ ij denotes the Kronecker delta.The SYK Hamiltonian is an all-to-all Hamiltonian, which couples all the Majorana fermions together in a fully non-local fashion, consisting of the following q-body interaction terms with q ≥ 2 being an even integer: where the coupling constants J i1...iq are randomly sampled from the Gaussian distribution of mean 0 and variance where J 2 is a constant which we set to be equal to one.The model has recently attracted a widespread attention from different communities, due to some peculiar features it enjoys.It has been shown that when q ≥ 4 the model is highly chaotic [17,18,29,30], although solvable in the large n limit [17,18], thus creating a perfect situation to study relevant questions on quantum chaos which are usually out of reach for other chaotic models.Moreover, the SYK model has intriguing connections with the physics of black holes and quantum gravity, promoting itself as an ideal candidate to address new questions on holography and the AdS/CFT correspondence.
We will focus our attention on the SYK model with q = 4.It has two notable features that can be a source of trouble for any eigensolver algorithms, regardless of being classical or quantum mechanical, whose goal is to reach the ground state.First, the energy spectrum of the SYK model is very dense, especially having a small energy gap between the ground state E 0 and all other excited states.Second, the SYK ground state is much less distinguishable from generic quantum states in the Hilbert space, supporting the volume law scaling of the entanglement entropy.Therefore, the VQE computation with the SYK model must be seen as a highly challenging benchmark [31].
As another characteristic, the SYK model is known to have two-fold degenerate eigenstates, if and only if n = 4k + 2 for any positive integer k [29,30].Denoting the two-fold degenerate ground states by |φ 1 and |φ 2 , that we assume to be normalized and orthogonal to each other, the VQE target states for the SYK model are then given by all the possible linear combinations of the form: with α and β being complex numbers.Therefore, the distance between the circuit state |ψ(θ) and the closest state of the form |φ 0 can be measured by computing with the inequality which is saturated whenever |ψ(θ) takes exactly the form (22).We will numerically show that the high-depth quantum circuit can effectively learn the ground state of the SYK model.It will imply that even complicated systems, involving non-local interactions and a high level of entanglement, can be universally simulated through the variational circuits.

Optimizing the Circuit
We start by discussing if the SYK Hamiltonian causes the vanishing gradient problem for the variational circuit in Figure 1.When at least one of U ± (θ ± ) is a 2-design, the scaling behavior of Tr(H 2 ) will determine if the SYK energy gradient will be exponentially suppressed or not.The Hamiltonian (20) has been exactly diagonalized up to n = 15.Moreover, an approximate analytical formula of the spectral density is at disposal in [32], which we use to extrapolate Tr(H 2 ) for the large system size n.The results are presented in Figure 2. We clearly see, as in the Ising model, that Tr(H 2 ) scales like 2 n , which indicates the abundance of the barren plateaus in the VQE energy landscape of the SYK model.The vanishing gradient problem has also been observed numerically by computing the sample variance of ∂ k E(θ) over a collection of 1000 random parameters, as displayed in Figure 8a.Clearly, the energy gradient ∂ k E(θ) is exponentially suppressed by increasing n.We also see that for the SYK model, contrary to the Ising Hamiltonian, there is almost no transient regime where ∂ k E(θ) is still large, yet decreasing for the growing number L of layers.The lack of the transient regime, which is a consequence of the non-local nature of the SYK interactions [11], is a clear obstacle in using the variational circuit to approximate the SYK ground state in the low-depth regime.
On the other hand, as we can see from Figure 8b, the norm of the gradient vector is again increasing for the growing number L of layers, due to the saturation of the Var θ [∂ k E(θ)] with respect to L. The empirically measured growth rates between log ∇ θ E(θ) and log L are of the SYK Hamiltonian, by repeating the same numerical experiments conducted in Section III 2 under the same choice of hyperparameters (α, β 1 , β 2 ) = (0.05, 0.9, 0.999).First, the sample distribution of the minimized energy E(θ * ) for randomly initialized circuits with 4 ≤ L ≤ 16 layers is illustrated in Figure 9, based on 35 sample VQE runs at each L. A notable observation is the small variance of the final energy, compared to the Ising VQE energy distribution in Figure 4, even at very low depth.We interpret it as another manifestation of the fact that the ground state that we aim to approximate is less distinguishable from generic quantum states.Hence, developing an optimization trajectory demands more classical computing power through the high-dimensional θ-space [31].Another -and perhaps the most important -point to stress is that, the mean value of the minimized circuit energy E(θ * ) decreases by stacking more layers, just like what has happened to the Ising VQE problem.It suggests that the high-depth circuits can reach a very good approximation of the SYK ground states.
Second, we have measured the performance of the layered circuit Ansatz (2) in approximating a ground state of the SYK model (20), as a function of the circuit depth.Specifically, the VQE single-run error E(θ * )−E 0 is drawn We also note from Figure 10b that the fidelity between the optimized and ground states tends to decrease at an intermediate scale of depth, while the VQE error continues to reduce without temporary increase.Such contrasting behavior is due to the dense energy spectrum of the SYK model near the ground energy level E 0 .With an intermediate-depth circuit, the VQE algorithm is going to approximate not the exact ground states, but some low-lying excited states.In this way, the energy continues to decrease but the fidelity does not improve.However, for sufficiently deep circuits, the VQE algorithm can overcome the difficulty and reach an excellent agreement with the ground state.Interestingly, we have seen that the necessary number L v of layers to reach the high precision (15) is roughly in the same order, both for the SYK and Ising models.We will also see that the circuit (2) with L ≥ L v achieves high precision in replicating random states by minimizing the Euclidean distance, as shown in Figure 11.This compatibility highlights the universal effectiveness of the high-depth circuit in approximating any quantum state in the Hilbert space, both generic and non-generic ones.

V. APPROXIMATING RANDOM STATES USING EUCLIDEAN LOSS FUNCTION
Taking one step further, we will look into the highdepth circuit's capability to approximate generic quantum states |φ by minimizing the Euclidean distance.
Recall that the variational circuit |ψ(θ) with L layers depends on nL parameters, which are periodic over the finite interval [0, 2π) up to an overall sign.It is a mapping from the nL-dimensional torus to the Hilbert space of n qubits.As a practical measure for the expressibility and trainability of the circuits in simulating quantum states, we consider whether, for a given quantum state |φ , there exist a point θ * in the parameter space, reachable by the gradient-based optimization such that |ψ(θ * ) |φ .Such measure of expressibility is motivated by [33] but tailored for the hybrid algorithms.Its definition follows: Let us apply the gradient descent to find the minimum distance at the closest point between the circuit and target states, where • is the Euclidean norm of a complex vector.The (in)expressibility of the variational circuit is written as the average of the minimum distances over all quantum states: Notice that the closeness between two quantum states is usually defined by using the trace distance, which is highly sensitive to small parameter changes.Instead, we have adopted the Euclidean distance to improve the convexity of the optimization landscape, such that the closest point θ * can be found as easily as possible.The Euclidean norm defines trivially a convex landscape.Hence, any non-convexity of the optimization landscape is inherited from the variational circuit itself that fixes how to embed the nL-dimensional torus into the Hilbert space.
For an actual estimation of ε, we substitute the Haar unitary integral (25) with the sample mean over m states: Given a target state |φ i , we start with an Ansatz state |ψ(θ) with nL initial parameters θ, randomly sampled from the uniform distribution U(0, 2π) ⊗nL .To reach the optimal parameters θ * that minimize the distance, we use the Adam optimization algorithm with the hyperparameters (α, β 1 , β 2 ) = (0.05, 0.9, 0.999).Note that the overwhelming majority of Hilbert space is filled with generic states with a high degree of entanglement, so the sample states |φ i will almost never be non-generic, lowentangling states similar to the Ising ground state.
Figure 11 shows the sample mean ε m over m = 10 random target states, divided by the number 2 n of state components.The narrow shade displays the fluctuation range of the component-averaged distance |ψ(θ * ) −|φ i • 2 −n across different target states indexed by i = 1, . . ., 10.We have selected the minimum distance among the first τ ≤ 500 optimization steps, which is sufficient since the optimal parameter θ * is empirically always found in an early stage of the optimization.
We find that the (in)expressibility approaches to 0 as the circuit depth L grows.To achieve the following high precision, the depth L has to be greater than the threshold values: Note that the slight ramp-up after the steep fall of ε m=10 is caused by the fluctuation of θ around the optimal point θ * , since the fluctuation size gets amplified for bigger L.
It is a common phenomenon in gradient-based optimization, where a suitable reduction of the learning rate can turn it to the stable convergence θ → θ * .
In Figure 11, we have found the following overall trend: Deeper circuits are not only a superset of shallow circuits but also behave more effectively in the gradient-based optimization.Though the functional form of the variational circuit in Figure 1 is composed of only two types of gates, it can accurately express and reach most quantum states in Hilbert space in the high-depth regime.

VI. DISCUSSION
From the viewpoint of the vanishing gradient problem, adding more parametric gates to the circuit architecture brings both positive and negative effects on trainability.On the negative side, it makes the circuit unitary ensemble closer to quantum 2-design, whose one-and two-point correlators are equivalent to those of Haar random unitaries, so that initial random gradients may decay exponentially with the system size n [4].On the positive side, however, it enlarges the dimensionality of the parameter space and maintains the norm of random gradients to be finite.
Until now, we have explored one limiting regime, where the impact of exponentially many parameters dominates over the barren plateau phenomenon.It turns out that the high-depth circuit has been very effective in solving the ground states of the Ising and SYK Hamiltonians, as well as in simulating generic random states.One may note a certain degree of qualitative similarity between the VQE optimization trajectory of high-depth circuits, demonstrated in Section III 2, and the lazy learning [34] in over-parameterized neural networks.We speculate it as naturally emerging in any systems involving the highdimensional parameter space.It would be interesting to know why local extrema on the energy landscape of the high-depth circuits can reliably reach a zero VQE error, as studied in [35] for the neural networks.
Some interesting phenomena that we found during the optimization of low-and intermediate-depth circuits need to be further understood.For low-depth circuits, it is important to characterize what features of the initial points lead to the observed big difference in their minimized energy levels.We are also curious to know what makes the local minima on the energy landscape of intermediatedepth circuits so homogeneous.More generally, we want to understand how the circuit states evolve along a gradient optimization trajectory on average, in terms of various quantum information measures.
We need to concern two types of errors for the actual use of the variational circuits on the near-term quantum devices [36].Firstly, for accurate estimation of the energy gradient, it is necessary to sample the variational wavefunction repeatedly, ideally exponentially many times in the system size n.Under the assumption that we can collect only a limited number of samples, the energy gradient estimation will be inevitably noisy.Therefore, the expressibility and trainability of the variational circuit needs to be addressed with noisy gradients, or even without direct gradient computation as in [37,38].We expect the simplicity of the energy landscape in the high-depth regime may bring robustness against the sampling noise.Secondly, and more severely, quantum hardware noise restricts our ability to maintain quantum states when a sequence of layer unitary acts on them.A noise-induced mechanism that causes the vanishing gradient problem has also been studied in [39].We would like to investigate if the variational circuit approach can be successful under decoherence of quantum states in the near future.
FIG. 2. Tr(H 2) of the Ising model (3) and the SYK model(20) as a function of the system size n.The dashed lines are the regression lines with the functional form of a • n b 2 n , based on the numerical data denoted as small circles.
FIG. 6. Optimization curves of the circuit state on the Ising model at 8 qubits.The upper and lower plots denote the VQE error E(θτ ) − E0 and the fidelity | ψ(θ * )|φ | 2 between the circuit and true ground states, respectively, as a function of the optimization steps τ .Notice the late-time fluctuation has been alleviated with the learning rate scheduling.

FIG. 11 .
FIG. 11.The parametric (in)expressibility of the layered circuit Ansatz over m = 10 random target samples, εm=10/2 n , normalized by the number 2 n of state components.The narrow shade denotes the fluctuation across m = 10 target states.
Sample average of the Euclidean norm of ∇ θ E(θ) FIG. 3. The barren plateau experiment for the Ising Hamiltonian (3).The shades indicate (a) the variance across gradient components {∂ k E(θ)} nL k=1 , (b) the first and third sample quantiles.
Sample average of the Euclidean norm of ∇ θ E(θ) FIG. 8.The barren plateau experiment for the SYK Hamiltonian (20).The shades indicate (a) the variance across gradient components {∂ k E(θ)} nL k=1 , (b) the first and third sample quantiles.