Variational optimization for quantum problems using deep generative networks

The VGON model

The architecture of VGON, shown in Fig. 1, consists of two deep feed-forward neural networks, the encoder E_ω and the decoder D_ϕ are connected via a latent layer \({{{\mathcal{Z}}}}\) containing a normal distribution \({{{\mathcal{N}}}}({{{\boldsymbol{\mu }}}}({{{\boldsymbol{z}}}}),{{{{\boldsymbol{\sigma }}}}}^{{{{\boldsymbol{2}}}}}({{{\boldsymbol{z}}}}))\), where the mean μ and the standard deviation σ are provided by E_ω. During the training stage, input data x₀ is sampled from a distribution P(x₀), which in all our tests is the uniform distribution over the parameter space. It is then mapped to the latent distribution \({{{\mathcal{N}}}}({{{\boldsymbol{\mu }}}}({{{\boldsymbol{z}}}}),{{{{\boldsymbol{\sigma }}}}}^{{{{\boldsymbol{2}}}}}({{{\boldsymbol{z}}}}))\) by the encoder network E_ω. Next the decoder network D_ϕ(z) maps data z sampled from the latent distribution \({{{\mathcal{N}}}}({{{\boldsymbol{\mu }}}}({{{\boldsymbol{z}}}}),{{{{\boldsymbol{\sigma }}}}}^{{{{\boldsymbol{2}}}}}({{{\boldsymbol{z}}}}))\) to a distribution minimizing the objective function h(x). This minimization is achieved by iteratively updating the parameters ω and ϕ in E_ω and D_ϕ, respectively. Due to the existence of a stochastic latent layer, the gradients cannot be propagated backwards in the network. We solve this issue by using the reparameterization trick²⁴.

The key difference between VGON and VAE lies in the objective function (also called the loss function in the machine learning literature): instead of asking the output data distribution to approximate the input distribution by maximizing a given similarity measure, VGON simply requires the output data to minimize any objective function that is appropriate for the target problem. In addition, the Kullback-Leibler (KL) divergence between the latent distribution \({{{\mathcal{N}}}}({{{\boldsymbol{\mu }}}}({{{\boldsymbol{z}}}}),{{{{\boldsymbol{\sigma }}}}}^{{{{\boldsymbol{2}}}}}({{{\boldsymbol{z}}}}))\) and the standard normal distribution \({{{\mathcal{N}}}}({{{\boldsymbol{0}}}},{{{\boldsymbol{I}}}})\) is also minimized during training, as part of the objective function.

After the objective function converges to within a given tolerance, the training stage is complete. To utilize the trained model, the encoder network is disabled, and data sampled from a standard normal distribution \({{{\mathcal{N}}}}({{{\boldsymbol{0}}}},{{{\boldsymbol{I}}}})\) are fed into the decoder network. Depending on the characteristics of the objective function, the corresponding output distribution can be tightly centered around one or multiple optimum values. Moreover, it is worth mentioning that requiring the latent layer to follow a normal distribution not only facilitates efficient optimization of the objective function but also simplifies the sampling process, since the KL divergence between two normal distributions can be analytically evaluated and sampling from a normal distribution is computationally efficient.

The goal of VGON can be seen as finding a way which maps a simple distribution defined over the latent space, i.e., the Gaussian, to a complicated one, i.e., a distribution whose samples can minimize the objective function with high probability²⁹. This shares the spirit of transport theory, where given two probability measures μ and ν on spaces X and Y, we call a map T: X → Y a transport map if T_*(μ) = ν, where T_*(μ) is the pushforward of μ by T, representing the process of transferring (or “pushing forward”) the measure μ from X to Y via the measurable function T. In an optimal transport problem³⁰, one is interested in finding a map T that minimizes the transport cost ∫ c(x, T_*(x))dμ(x), subject to the constraint that the pushforward measure satisfies T_*(μ) = ν³¹. In addition to optimal transport problems, there are problems that do not have an explicitly defined target distribution, where the task is evaluating the loss function directly on the generated samples. In these problems, an optimal T can either be found analytically, such as in inverse transform sampling, where both distributions are one-dimensional³², or T can be learned/optimized from a parameterized T_θ on training data, such as generative model like normalizing flows^22,33, where both distributions are high-dimensional but T is invertible and constructed using neural networks. In VGON, the latent space usually has a much lower dimension than the output layer, making T surjective, which means every point in the target space can be reached by applying the decoder to some latent input³⁴. This surjectivity relaxes the requirement for invertibility and enables VGON to easily cover the complex target distributions. In our experience, the best optimization results come from training the encoder-latent space-decoder triple as a whole, even though it is possible to achieve good results without including the encoder in the training process.

To show that VGON can work well, we first use it to solve a variational optimization problem with a known unique optimal solution: finding the minimum ground state energy density among a class of quantum many-body models that matches the lower bound certified by an SDP relaxation^3,35. More specifically, we consider a class of infinite 1D translation-invariant (TI) models with fixed couplings³⁶, and the optimization variables are the local observables. The ground state energy density of this class of models has a lower bound certified by a variant of the NPA hierarchy³⁶. However, there is no guarantee that any infinite TI quantum many-body Hamiltonian, if couplings are fixed but the local observables can be arbitrary, can achieve this bound. Meanwhile, by optimizing 3-dimensional local observables with SGD and computing the ground state with uniform matrix product state algorithms, a Hamiltonian whose ground state energy density matches the above lower bound to seven significant digits has been found. We replace SGD with VGON to conduct the same optimization, and find that the converged model can (almost) deterministically generate Hamiltonians whose ground state energy density matches the NPA lower bound to eight significant digits, reaching the precision limit of commercial SDP solvers (see Supplementary Note 1 for more technical details). Below we apply VGON to several more complicated problems.

Finding the optimal state for entanglement detection

Entanglement detection plays a central role in quantum tasks such as secure communication and distributed computing, where entanglement serves as a fundamental resource. Suppose two players, Alice and Bob, receive a bipartite quantum state ρ from a source, then they want to determine whether ρ is entangled, with the smallest statistical error.

They can either perform the experiment independently in their respective laboratories and subsequently communicate the outcomes from Alice to Bob, or choose to forgo communication entirely. In the first scenario they are implementing a local operations and one-way classical communication (1-LOCC) protocol while in the second scenario they are implementing a local operations (LO) one. Implementing the 1-LOCC protocol experimentally requires fast real-time switching of Bob’s measurement settings and a quantum memory to store Bob’s half of the state while Alice performs her measurement and communicates the result. Do these extra experimental complexities yield tangible advantages such as reduced statistical error probabilities? In fact, it has been shown that for some simple states, such an advantage does exist, but it is too small to be useful³⁷. In fact, the advantage is highly dependent on the choice of target states and is hard to estimate analytically. The goal of this task is to identify high-dimensional quantum states for which this advantage, defined as the gap between the minimum statistical error probabilities in LO and 1-LOCC protocols, is as large as possible. This gap quantifies the practical advantage of allowing one-way communication between the parties in entanglement detection.

Specifically, given a quantum state ρ, the advantage is defined to be the gap between the minimum probabilities p₂ of committing false-negative errors (a.k.a. type-II errors, defined as a source distributes an entangled state, but Alice and Bob conclude the state they received is separable) when Alice and Bob employ LO and 1-LOCC protocols, and the two protocols have the same probability of making false-positive errors, i.e., type-I error, denoted p₁ and defined as the scenario in which they conclude that they have received an entangled state, while the source actually distributes a separable one. The desired quantity can be computed by solving two SDPs that share the same objective but differ in constraint structure. Specifically, the following SDP is solved twice, once under LO protocol \(P\in {{{{\mathcal{P}}}}}^{LO}\) and once under 1-LOCC protocol \(P\in {{{{\mathcal{P}}}}}^{1{\mbox{-}}LOCC}\). The final result is obtained by subtracting the respective optimal values of p₂³⁸:

Here \({{{{\mathcal{S}}}}}^{* }\) denotes the dual of the separable set \({{{\mathcal{S}}}}\). The positive operator-valued measure (POVM) operators M_Y(N)(P) can be constructed as \({M}_{Y(N)}(P)={\sum }_{x,y,a,b}P(x,y,Y(N)| a,b)\,{A}_{x}^{a}\otimes {B}_{y}^{b}\), where \(\{{A}_{x}^{a}\}\) (\(\{{B}_{y}^{b}\}\)) are predetermined measurements performed by Alice (Bob) with x (y) being measurement labels and a (b) being outcomes, and P is the shorthand for the distribution P(x, y, Y(N)∣a, b), which specifies the detection strategy by assigning probabilities to particular combinations of settings, outcomes, and decisions. Here Y and N denote the decisions corresponding to the presence or absence of entanglement, respectively.

For a random quantum state ρ, it turns out that the gap calculated above is usually very small, as shown by the green dots in Fig. 2a. In order to observe the gap under noisy experimental conditions, We focus on a linear optical setup that generates bipartite qutrit states, which also allows us to parameterize the state space. We first employ SGD to maximize the gap by starting from random pure bipartite qutrit states. The results are shown in Fig. 2a. The SGD algorithm gets trapped easily in local maxima and needs to compute the gradient by solving dozens of SDPs. Optimizing the gap for 79,663 random pure states is computationally very costly (see Supplementary Note 2.1 for more details). Most of these states exhibit gaps around 0.0036 before optimization, while the largest gap afterwards is 0.083722.

The results of using VGON to maximize the gap are depicted in Fig. 2b and summarized in Table 1. Based on 3000 sets of initial parameters produced by uniform sampling, the model converges after less than 2 h of training. After that, we use it to generate 10,000 output states. We find that over 98% of them manifest gaps over 0.08, while over 50% of them have gaps larger than 0.0835. A similar performance has been observed even when choosing the initial quantum states from a variational submanifold of the space of all mixed states, where out of 10,000 states generated by the converged VGON model, 83 have gaps larger than 0.0837, with an average purity of 0.99999. For comparison, we also apply seven other global optimization algorithms and a Multilayer Perceptron (MLP)—a neural network with multiple fully connected layers—to this task, and find that VGON consistently outperforms all the baseline methods (see Table 2 and Supplementary Note 2 for more details). Importantly, using VGON to solve an experiment-relevant instance of this problem allows us to experimentally demonstrate the advantage of 1-LOCC protocols over LO ones in detecting entanglement, where we can observe an error probability that is impossible to achieve with local operations alone³⁸.

Alleviating the effect of barren plateaux in variational quantum algorithms

On problems with a moderate size of optimization variables, VGON has shown its ability to quickly converge to the (nearly) optimal output distribution and generate high quality solutions with high probability. In near-term hybrid quantum-classical algorithms such as the variational quantum eigensolver (VQE)⁵, however, the number of classical parameters can quickly reach thousands or tens of thousands. The performance of such a hybrid algorithm can be hard to predict. On the classical part, the problems of vanishing gradients and having multiple minima are often present^{39,40,41,42,43}. On the quantum part, the choice of ansätze greatly affects the expressivity of the quantum circuit, making the certification of global optimality difficult^{44,45,46,47,48}.

For example, in a typical VQE algorithm, a parameterized variational circuit U(θ) is used to approximately generate the ground state of a target Hamiltonian H. The circuit structure usually loosely follows the target Hamiltonian and is often called an ansatz. Then by setting the energy of the output state \(\left\vert \psi ({{{\boldsymbol{\theta }}}})\right\rangle =U({{{\boldsymbol{\theta }}}})\left\vert 00\cdots 0\right\rangle\) with respect to H, i.e., 〈ψ(θ)∣H∣ψ(θ)〉, as the objective function, the algorithm iteratively updates the parameters in the quantum circuit by applying gradient-based methods on a classical computer. When the algorithm converges, the output quantum state will likely be very close to the ground state of H.

However, when the size of quantum systems increases, gradients vanish exponentially. This is primarily because the random initializations of parameterized unitaries conform to the statistics of a unitary 2-design^41,49, making the working of gradient-based optimization difficult. To overcome the BP problem, several strategies have been proposed and investigated, with the small-angle initialization (VQE-SA) method being identified as an effective technique^50,51,52. It initializes parameters θ to be close to the zero vector, which differs from the statistics of the parameters from a 2-design and thus may alleviate the BP problem.

The advantage of VGON over VQE-SA in alleviating BPs can be seen when we use them both, with the same parameterized quantum circuit, to compute the ground state energy of the Heisenberg XXZ model. Its Hamiltonian with periodic boundary conditions is given by:

where \({\sigma }_{x,y,z}^{i}\) denote the Pauli operators at site i. The ansatz for the parameterized quantum circuit is inspired by the matrix product state encoding⁵³. It consists of sequential blocks of nearest-neighbor unitary gates, each of which is made of 11 layers of single qubit rotations and CNOT gates (see Supplementary Note 3.2 for more details).

By choosing N = 18, the circuit contains 816 blocks and 12,240 variational parameters. The average ground state energy, computed using exact diagonalization (ED), is −1.7828. We use both VQE-SA and VGON to compute the same quantity, with each method repeated 10 times. The results are shown in Table 3 and Fig. 3, where the mean values and the 95% confidence intervals of these methods are visualized. The dark-blue and the green lines represent the average energy for VQE-SA and VGON, whose mean values at the last iteration are −1.7613 and −1.7802, respectively. Furthermore, to compare the performance of the two methods in a more fine-grained manner, we also calculate the fidelity between the states produced by the quantum circuit and the exact ground state. As the purple line depicted, VGON can achieve a 99% fidelity by around 880 iterations, while the VQE-SA method can only achieve 78.25% fidelity within the same number of iterations. We would like to remark that since the quantum computational resource consumed in VGON is similar to that in VQE-SA and classical computational resource is relatively cheap, the wall clock time cost per iteration for VGON can be comparable to that of VQE-SA. Moreover, the batch training of VGON can lead to a more stable convergence. For practical on-hardware quantum optimization, batch evaluation also enhances noise robustness by averaging over multiple circuit executions, thereby mitigating stochastic fluctuations, suppressing outliers, and accelerating convergence. For another comparison, VQE with uniformly random initial parameters can barely provide meaningful results due to the presence of barren plateaux, where the mean value of the average energy is −0.1367 after 1000 iterations across ten repetitions, as illustrated by the light-blue line in Fig. 3.

Identifying degenerate ground state space of quantum models

Deterministic gradient-based optimization methods are predisposed to follow a single path, therefore hampering their ability to efficiently detect multiple optima. A unique advantage of generative models is the ability to produce diverse samples of objects, all of which may minimize the objective function. In optimization, this leads to the possibility of finding multiple optimal solutions with a single stage of training. We now show that when appropriately trained, VGON exemplifies such an effective capability for generating multiple (nearly) optimal solutions simultaneously. This capability can be largely ascribed to its integration of randomness and the adoption of batch training. The former facilitates broader exploration within the variational manifold, while the latter, which involves processing subsets of data samples concurrently, supports the collective identification of multiple optimal solutions.

A natural multi-optima problem in quantum many-body physics is the exploration of degenerate ground spaces of quantum many-body Hamiltonians. We apply VGON to two Hamiltonians with known degenerate ground states: the Majumdar-Ghosh (MG)^54,55 model in Eq. (2), and a Heisenberg-like model in Eq. (3) coming from one of the contextuality witnesses presented in ref. ³⁷:

where \({{{{\boldsymbol{\sigma }}}}}^{i}=({\sigma }_{x}^{i},{\sigma }_{y}^{i},{\sigma }_{z}^{i})\) are Pauli operators at site i. We take N = 10 for H_MG, and N = 11 for H₂₃₂, making their ground state spaces five- and two-fold degenerate, respectively. An orthonormal basis for their respective degenerate ground state spaces is computed by the ED method, which outputs five vectors \(\left\vert {v}_{1}\right\rangle \ldots \left\vert {v}_{5}\right\rangle\) spanning the ground state space of H_MG, and two vectors \(\left\vert {u}_{1}\right\rangle\) and \(\left\vert {u}_{2}\right\rangle\) spanning that of H₂₃₂.

The overall objective of this problem is similar to the previous one: finding the ground states of H_MG and H₂₃₂ with variational quantum circuits. We maintain the same circuit layout as in the previous problem, and use 36 and 60 blocks of unitary gates for each Hamiltonian respectively. Profiting from the use of mini-batches to estimate gradients, a common technique in training neural networks, VGON can effectively evaluate many different circuits simultaneously. Meanwhile, to enhance intra-batch diversity, a penalty term consisting of the mean cosine similarity among all pairs of sets of circuit parameters in the same batch is added to the objective function. This penalty term, together with the mean energy of the states in the batch, ensures a balance between maintaining the diversity of the generated outputs and minimizing the energy. Further details can be found in Supplementary Note 3.1.

As a result, unlike VQE-based algorithms aiming to generate multiple energy eigenstates, the objective function of VGON is model-agnostic. In other words, with no prior knowledge of the degeneracy of the ground space, VGON is capable of generating orthogonal or linearly independent states in it. In comparison, to achieve a diversity of outputs with VQE-based algorithms^56,57, it is essential to provide diverse inputs for the VQE model. However, attaining this diversity can result in barren plateaux within the optimization landscape. Though employing VQE-SA may address this problem, it could significantly diminish the diversity of inputs, as it tends to constrain inputs to values near zero.

We generate 1000 output states for each Hamiltonian using a VGON model trained with the above objective function. We find that the vast majority of these states have energy low enough to be treated as ground states. Figure 4 shows the overlap between the generated states and the basis of their ground state space. In Fig. 4a, the generated states for H₂₃₂ fall into two orthogonal classes, which form an orthonormal basis of the ground state space. For H_MG, Fig. 4b shows that most of them are linearly independent and span the same space as \(\left\vert {v}_{1}\right\rangle \ldots \left\vert {v}_{5}\right\rangle\).