Universal Effectiveness of High-Depth Circuits in Variational Eigenproblems

It was argued in McClean et al. (2018) 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 McClean et al. (2018); Skolik et al. (2020). 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, $\partial_{k}E(\bm{\theta})$, that belongs to the $\ell$’th variational layer, it is convenient to group the $L$ layer unitaries (2) into the following two blocks:

where $\bm{\theta}_{q}=\{\theta_{p}\,|\,\lfloor\frac{p}{n}\rfloor=q\}$, $\bm{\theta}_{+}=\bigcup_{a\leq\ell}\bm{\theta}_{a}$, $\bm{\theta}_{-}=\bigcup_{a>\ell}\bm{\theta}_{a}$,

As the partial derivative $\partial_{k}$ acts only on $U_{\ell}(\bm{\theta}_{\ell})$, the variance of the $k$’th gradient component, $\text{Var}_{\bm{\theta}}[\partial_{k}E(\bm{\theta})]$, is

where $V_{k}$ denotes the Pauli operator $\sigma^{y}$ acting on the $k$’th spin variable, such that $\text{Tr}(V_{k})=0$ and $\text{Tr}(V_{k}^{2})=2^{n}$.

If we assume the quantum 2-design property of $U_{+}(\bm{\theta}_{+})$ and/or $U_{-}(\bm{\theta}_{-})$, 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

Such simplification (7) happens when both $U_{\pm}(\bm{\theta}_{\pm})$ are 2-designs. Instead, if the 2-design condition does not hold for either of $U_{\pm}(\bm{\theta}_{\pm})$, we have the following expression:

if only $U_{+}(\bm{\theta}_{+})$ is a 2-design while $U_{-}(\bm{\theta}_{-})$ is not, or

where $\rho=|0\rangle\langle 0|$, if not $U_{+}(\bm{\theta}_{+})$ but only $U_{-}(\bm{\theta}_{-})$ 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$ Yang and Feng (2002):

To determine the scaling behavior of the upper bound (10) with respect to the system size $n$, we examine how $\text{Tr}(\mathcal{H}^{2})$ scales with $n$ by extrapolating the values obtained from exact diagonalization of the Hamiltonian (3) up to $n\leq 17$. The result is presented in Figure 2. We see that the term $\text{Tr}(\mathcal{H}^{2})$ scales as $2^{n}$, such that the exponential factor in the denominator of the upper bound (10) cannot fully be balanced by $\text{Tr}(\mathcal{H}^{2})$. Then, inserting the upper bound (10) into Chebyshev’s inequality, which states

with $\partial_{k}E(\bm{\theta})$ 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_{\pm}(\bm{\theta}_{\pm})$ 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\leq k\leq 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\geq 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 $\|\nabla_{\bm{\theta}}E(\bm{\theta})\|$, which appears in the evolution of the circuit energy under the gradient descent with an infinitesimal rate $\alpha$:

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 $\bm{\theta}$-space.

In agreement with the reasoning above, Figure 3b exhibits a small initial drop dominated by the transient decrease of $\text{Var}_{\bm{\theta}}[\partial_{k}E(\bm{\theta})]$ for $L\leq 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 $\|\nabla_{\bm{\theta}}E(\bm{\theta})\|$ as

which agrees well with the numerically estimated growth rates between $\log\|\nabla_{\bm{\theta}}E(\bm{\theta})\|$ 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$.

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 $\bm{\theta}$, can be summarized to the following two features:

1. 1.

   For small enough $L$, the minimized energies $E(\bm{\theta^{*}})$ that the circuit states can reach are highly variable.

2. 2.

   For larger $L$, the $E(\bm{\theta^{*}})$ 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\in\{8,10,12,14,16,18,20\}$ layers. To find the optimal parameters $\bm{\theta}^{*}$ that minimize the energy, we use the Adam optimization algorithm Kingma and Ba (2017), which iteratively updates the variational parameters $\bm{\theta}$ 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 $(\alpha,\beta_{1},\beta_{2})$. We have collected 35 independent VQE runs for the Ising Hamiltonian (3), whose initial parameters are randomly sampled from the uniform distribution $\mathcal{U}(0,2\pi)^{\otimes nL}$, after 500 parameter updates with the hyperparameters $(\alpha,\beta_{1},\beta_{2})=(0.05,0.9,0.999)$.

The sample distribution of the final VQE energy $E(\bm{\theta^{*}})$ 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 $\bm{\theta^{*}}$ clearly exhibits the widespread spectrum for a shallow depth, e.g.,

sometimes achieving a relatively good energy level while the gap always persists. Second, by stacking more layers, i.e., $L\geq 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(\bm{\theta^{*}})$ 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(\bm{\theta^{*}})$, 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 $|\psi_{0}\rangle$, reaching the ground state energy $E_{0}$ and achieving a high level of fidelity with $|\psi_{0}\rangle$. The high-depth circuits can achieve a zero error with the following precision:

where $\Delta E$ denotes the bandwidth of the target Hamiltonian $\mathcal{H}$, defined as the difference between the largest and smallest eigenvalues of $\mathcal{H}$. We note that such precision can be achieved only with an appropriately chosen learning rate $\alpha$, 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 $\alpha=0.05$ ($L=4,6,8$) and $\alpha=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 $\bm{\theta}$ has more dimensions, the basin of attractor to local extrema becomes narrower Cerezo et al. (2020), giving a larger value of the estimated inverse volume Wu et al. (2017)

where the summation is taken over the top-$k$ eigenvalues $\{{\lambda_{i}(H)}\}_{i=1}^{k}$ of the Hessian matrix $H_{ij}\equiv\partial_{i}\partial_{j}E({\bm{\theta})}$. For instance, the positive correlation between $L$ and $V_{k}^{-1}(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 $|\psi_{0}\rangle$ of the Ising Hamiltonian (3), the narrowness of the attractor basin still makes the VQE trajectory somewhat unstoppable, passing through the best point $\bm{\theta^{*}}$ and then hopping around in the local neighborhood. Figure 6a shows that the VQE error $E({\bm{\theta}_{\tau}})-E_{0}$ slightly increases on average and mildly fluctuates after achieving the minimum error $E(\bm{\theta^{*}})-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 $\alpha$, i.e.

at the optimization step $\tau$ 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 $|\psi_{0}\rangle$. We will explore this efficiency of the high-depth circuit in some details by visualizing the VQE trajectory on the energy landscape.

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 $\bm{\theta}_{0}$ and the trajectory $\mathcal{T}(\bm{\theta}_{0})=\{\bm{\theta}_{\tau}\}$ thereafter, we identify the optimal parameter $\bm{\theta^{*}}$ as the point in the trajectory $\mathcal{T}(\bm{\theta}_{0})$ having the minimum energy expectation value,

which is the best possible representative point of the local extrema that the trajectory $\mathcal{T}(\bm{\theta}_{0})$ converges around. The associated basin of attractor can be estimated by calculating the Hessian $H_{ij}\equiv\partial_{i}\partial_{j}E(\bm{\theta^{*}})$ at the optimum $\bm{\theta^{*}}$. We are interested in the eigenvalue spectrum, $\{\lambda_{i}(H)\}$, from which we distinguish steep and flat directions and calculate the degree of steepness.

Similar to the case of deep neural networks Sagun et al. (2017); Chaudhari et al. (2017); Wu et al. (2017), the local extrema in the VQE energy landscape of high-depth 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 $\mathcal{S}_{k}(\bm{\theta^{*}})$ along the $k$ steepest directions, i.e., spanned by the top-$k$ Hessian eigenvectors. More precisely, we want to examine if the $\mathcal{S}_{k}(\bm{\theta^{*}})$-projected Euclidean distance between $\bm{\theta^{*}}$ and $\bm{\theta}_{\tau}$ decreases along the trajectory $\mathcal{T}(\bm{\theta}_{0})$ as the optimization step $\tau$ progresses, until $E(\bm{\theta}_{\tau})=E(\bm{\theta^{*}})$. This will tell us if the entire trajectory $\mathcal{T}(\bm{\theta}_{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(\bm{\theta}_{\tau})-E_{0}$ and the projection distance $\|\bm{\theta}_{\tau}-\bm{\theta^{*}}\|_{\mathcal{S}_{100}(\bm{\theta^{*}})}$ along the actual optimization trajectories, whose step number $\tau$ is indicated by color. We have made the following observations: First, when an appropriate value of $k$ is selected, both the distance $\|\bm{\theta}_{\tau}-\bm{\theta^{*}}\|_{\mathcal{S}_{k}(\bm{\theta^{*}})}$ and the loss value $E(\bm{\theta}_{\tau})-E_{0}$ continuously decrease on a macroscopic scale. It exhibits that the trajectory $\mathcal{T}(\bm{\theta}_{0})$ converges without escaping from a specific basin of the attractor that encloses a randomly initialized point $\bm{\theta}_{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 $\mathcal{S}_{k}(\bm{\theta^{*}})$ evolves rapidly steeper and narrower as the depth $L$ increases, thereby causing the rapid convergence and substantial late-time fluctuation around $\bm{\theta^{*}}$. 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 Kiani et al. (2020); Wiersema et al. (2020) under the VQE algorithm.

The SYK model Sachdev and Ye (1993); Kitaev (2015); Maldacena and Stanford (2016) is built out of $2n$ Majorana fermions in $1d$, i.e. the operators $\gamma_{i}$, with $i=1,\,\dots,2n$, satisfying the following anti-commutation relations

where $\delta_{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\geq 2$ being an even integer:

where the coupling constants $J_{i_{1}\dots i_{q}}$ 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\geq 4$ the model is highly chaotic Kitaev (2015); Maldacena and Stanford (2016); García-García and Verbaarschot (2016); Cotler et al. (2017), although solvable in the large $n$ limit Kitaev (2015); Maldacena and Stanford (2016), 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 Kim and Oz (2021).

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$ García-García and Verbaarschot (2016); Cotler et al. (2017). Denoting the two-fold degenerate ground states by $|\phi_{1}\rangle$ and $|\phi_{2}\rangle$, 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 $\alpha$ and $\beta$ being complex numbers. Therefore, the distance between the circuit state $|\psi(\bm{\theta})\rangle$ and the closest state of the form $|\phi_{0}\rangle$ can be measured by computing

with the inequality which is saturated whenever $|\psi(\bm{\theta})\rangle$ 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.

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_{\pm}(\bm{\theta}_{\pm})$ is a $2$-design, the scaling behavior of $\text{Tr}(\mathcal{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 García-García and Verbaarschot (2017), which we use to extrapolate $\text{Tr}(\mathcal{H}^{2})$ for the large system size $n$. The results are presented in Figure 2. We clearly see, as in the Ising model, that $\mathrm{Tr}(\mathcal{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 $\partial_{k}E(\bm{\theta})$ over a collection of $1000$ random parameters, as displayed in Figure 8a. Clearly, the energy gradient $\partial_{k}E(\bm{\theta})$ 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 $\partial_{k}E(\bm{\theta})$ 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 Cerezo et al. (2020), 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 $\text{Var}_{\bm{\theta}}[\partial_{k}E(\bm{\theta})]$ with respect to $L$. The empirically measured growth rates between $\log\|\nabla_{\bm{\theta}}E(\bm{\theta})\|$ and $\log L$ are

matching the simple estimation formula (14). It tells that the gradient-based optimization can at least launch for high-depth quantum circuits, which bypass the vanishing gradient problem using an exponentially high-dimensional parameter space.

As the next step, we examine the performance of the variational circuit Ansatz (2) in reaching the ground state of the SYK Hamiltonian, by repeating the same numerical experiments conducted in Section 2 under the same choice of hyperparameters $(\alpha,\beta_{1},\beta_{2})=(0.05,0.9,0.999)$.

First, the sample distribution of the minimized energy $E(\bm{\theta}^{*})$ for randomly initialized circuits with $4\leq L\leq 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 $\bm{\theta}$-space Kim and Oz (2021). Another – and perhaps the most important – point to stress is that, the mean value of the minimized circuit energy $E(\bm{\theta}^{*})$ 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({\bm{\theta}^{*}})-E_{0}$ is drawn in Figure 10. The high-depth circuit performs very well, reaching a zero error with the following accuracy,

when the depth $L$ arrives at the following values.

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\geq 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.