A path towards quantum advantage in training deep generative models with quantum annealers

Early numerical studies showed that Boltzmann machines (BMs) with the quasi two-dimensional connectivities typical of current-generation quantum annealers perform poorly [53], even on simple standard datasets such as MNIST [54]. A possible approach to circumvent the limitation above is the following quantum-classical hybrid (QCH) paradigm: 1) employ classical feed-forward NN to encode the data into a compressed representation that is more easily processed by a quantum device and 2) jointly train the NN and the quantum device. This approach was pioneered in reference [55], where the authors introduced a quantum-assisted Helmholtz machine (QAHM), a generative model with latent variables. A QAHM has an inference( encoder), a generation( decoder) network and a generative process that is represented by BM and can thus be simulated by a quantum annealer. QAHMs do not have a well-defined loss function [9], which means that training does not scale well, even to standard machine-learning datasets such as MNIST. This approach was taken a step further in reference [1] by using variational autoencoders (VAE), a class of generative models with latent variables that provide an efficient inference mechanism [56, 57]. VAEs can be trained by minimizing the evidence lower bound (ELBO), a variational lower bound to the exact log-likelihood. The ELBO is a well-defined loss function with fully propagating gradients that can be efficiently optimized via backpropagation. This allows us to achieve competitive and scalable performance on large-scale datasets. Upon discretization of the latent space [58, 59], the generative process can be then realized by CBM or QBM. A QCH variational autoencoder (QVAE) is a VAE with an integrated QBM. Crucially, QVAEs can also be trained by optimizing a lower bound to the exact log-likelihood [1].

Our hybrid approach exploits a large amount of classical computational resources (backpropagation through deep NN performed on GPUs). This is advantageous since it allows us to achieve competitive results on relatively large-scale datasets such as MNIST. However, validating the training with quantum annealers is a major challenge: VAEs can train well even with unstructured priors (e.g., a product of Gaussian or Bernoulli distributions). One thus needs to carefully validate and demonstrate any performance improvement due to the use of structured samples coming from well-trained BMs. To this end, existence of a well-defined loss function, which allows for a quantitative comparison between different models, is critical. In order to single out the performance improvement due to computations that can be off-loaded to a quantum annealer (sampling from the BM), we always focus on comparing the performance of a model trained with a nontrivially connected BM (for example with full or Chimera connectivity D) to a baseline in which we employ a BM with no connectivity; that is, a set of independent Bernoulli variables. After a certain number of gradient updates, training of such models will have used exactly the same amount of classical resources that cannot be off-loaded to the quantum annealer, with the only difference being the computational effort needed to correctly sample from the latent-space BM.

A complete description of our model will be given in section 2, more theoretical details about hybrid quantum-classical VAEs can be found in [1].

We successfully train a convolutional QVAE with 288 latent units and Chimera-structured QBM (a 6-by-6 Chimera graph; see appendix D) using only D-Wave 2000Q quantum annealers to draw the samples required to evaluate the gradients for the QBM parameters. We then use several techniques to validate the trained models. In particular, to perform quantitative evaluation we use auxiliary CBMs trained on the annealer samples to quantitatively estimate the test log-likelihood. Using these auxiliary CBMs is necessary in practice since we can easily evaluate the log-likelihood of reasonably large CBMs, while it would not be feasible to evaluate the log-likelihood of large QBMs. This approach is consistent because we empirically find that the quantum annealer samples from a QBM with a very small effective transverse field and it is thus close to simulating a CBM. With this approach, we show that models trained on Chimera connectivity outperform their trivial ('Bernoulli') baseline and can achieve competitive performance on the MNIST dataset.

The next question is whether our approach offers a path toward obtaining quantum advantage with quantum annealers in machine learning applications. To achieve such a goal, we need to exploit large latent-space BMs that develop complex multimodal probability distributions (i.e., a complex energy landscape) from which sampling is classically inefficient. We give evidence that this is indeed possible within our setup. For example, we demonstrate that the BMs placed in the latent space develop nontrivial modes that are likely to cause Monte Carlo algorithms to have long mixing times. Moreover, we show that training on more complex datasets likely takes advantage of larger BMs to improve performance. In addition, we discuss the role of connectivity, emphasizing its importance even in this hybrid approach and the necessity to develop device-specific classical NN to better exploit physical connectivities such as the Chimera graph.

The structure of the paper is as follows. In section 2 we review VAE and the implementations of discrete latent variables, a necessary step to implement BMs and QBMs in their latent space. Section 3 motivates the use of quantum annealers as samplers to train quantum and classical BMs. In section 4 we report our experiments in training VAEs with D-Wave 2000Q systems. In section 5 we discuss a possible path toward quantum advantage in our setup. Finally, we present our conclusions in section 6.

Training generative models with latent variables via MLE requires the evaluation of the intractable integral of equation (2) to calculate the posterior distribution $p_{{\boldsymbol{\theta}}}(\boldsymbol{\zeta}|\mathbf{x})$. VAEs circumvent this problem by introducing a tractable variational approximation $q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})$ to the true posterior [67], with variational parameters ${\boldsymbol{\phi}}$ (see figure 1). VAEs are then trained by maximizing a variational lower bound $\mathcal L(\mathbf{x},{\boldsymbol{\theta}},{\boldsymbol{\phi}})$ to the log-probabilities $\log p_{{\boldsymbol{\theta}}}(\mathbf{x})$:

$\begin{eqnarray} \mathcal L(\mathbf{x}, {\boldsymbol{\theta}},{\boldsymbol{\phi}}) & = & \log p_{{\boldsymbol{\theta}}}(\mathbf{x}) - \mathbb{E}_{\boldsymbol{\zeta}\sim q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})}\left[ \log \frac{q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})}{p_{{\boldsymbol{\theta}}}(\boldsymbol{\zeta}|\mathbf{x})}\right]\equiv \nonumber \\ &\equiv& \log p_{{\boldsymbol{\theta}}}(\mathbf{x}) - D_{KL}( q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x}) || p_{{\boldsymbol{\theta}}}(\boldsymbol{\zeta}|\mathbf{x}))\,, \end{eqnarray} \tag{ 3 }$

where $D_{KL}( q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x}) || p_{{\boldsymbol{\theta}}}(\boldsymbol{\zeta}|\mathbf{x}))$ is the Kullback-Leibler divergence (KL divergence) between the true and approximating posteriors. Since KL divergences are always non-negative, we have

$\begin{eqnarray} &&\mathcal L(\mathbf{x}, {\boldsymbol{\theta}},{\boldsymbol{\phi}}) \le \log p_{{\boldsymbol{\theta}}}(\mathbf{x}), \end{eqnarray} \tag{ 4 }$

which immediately gives:

$\begin{eqnarray} \mathcal L(\mathbf{X}, {\boldsymbol{\theta}},{\boldsymbol{\phi}}) & \equiv & \mathbb{E}_{\mathbf{x} \sim p_\text{data}}[ \mathcal L(\mathbf{x},{\boldsymbol{\theta}},{\boldsymbol{\phi}})] \nonumber \\ & \le & \mathbb{E}_{\mathbf{x} \sim p_\text{data}}[\log p_{{\boldsymbol{\theta}}}(\mathbf{x})] \equiv L(\mathbf{X}, {\boldsymbol{\theta}})\,, \end{eqnarray} \tag{ 5 }$

where $\mathcal L(\mathbf{X}, {\boldsymbol{\theta}},{\boldsymbol{\phi}})$ is called the ELBO. The ELBO can be written in terms of tractable quantities:

$\begin{eqnarray} \mathcal L(\mathbf{X}, {\boldsymbol{\theta}},{\boldsymbol{\phi}}) & = & \mathbb{E}_{\mathbf{x} \sim p_\text{data}}\Big[ \mathbb{E}_{\boldsymbol{\zeta}\sim q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})}[\log p_{{\boldsymbol{\theta}}}(\mathbf{x}|\boldsymbol{\zeta})] + \nonumber \\ &- & D_{KL}( q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x}) || p_{{\boldsymbol{\theta}}}(\boldsymbol{\zeta}))\Big]\,. \end{eqnarray} \tag{ 6 }$

The marginal $p_{{\boldsymbol{\theta}}}(\mathbf{x}|\boldsymbol{\zeta})$ and approximating posterior $q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})$, also called 'decoder' and 'encoder', respectively. In our set-up, both encoder and decoder are a product of Bernoulli distributions whose parameters (probabilities, or logits) are modeled by the convolutional NN described in section A.

To train VAEs, we need to calculate the derivatives of the objective function (equation (6)) with respect to the generative (${\boldsymbol{\theta}}$) and inference (${\boldsymbol{\phi}}$) parameters. The naive evaluation of $\partial_{{\boldsymbol{\phi}}}$ of terms of the type $\mathbb{E}_{\boldsymbol{\zeta} \sim q_\phi}[f(\boldsymbol{\zeta})]$ is called REINFORCE [68]. With the use of the identity $\partial_{{\boldsymbol{\phi}}} q_\phi = q_\phi \partial_{{\boldsymbol{\phi}}} \log q_\phi$, one has:

$\begin{eqnarray} &&\partial_{{\boldsymbol{\phi}}} \mathbb{E}_{\boldsymbol{\zeta} \sim q_\phi}\left[f(\boldsymbol{\zeta})\right] = \mathbb{E}_{\boldsymbol{\zeta} \sim q_\phi}\left[f(\boldsymbol{\zeta}) \partial_{{\boldsymbol{\phi}}} \log q_{\boldsymbol{\phi}} \right]\,. \end{eqnarray} \tag{ 7 }$

However, the term above has high variance and requires intricate variance-reduction mechanisms to be of practical use [69].

A better approach is to write the random variable $\boldsymbol{\zeta}$ as a deterministic function of the distribution parameters ${\boldsymbol{\phi}}$ and of an additional auxiliary random variable ${\boldsymbol{\rho}}$. The latter is given by a probability distribution $p({\boldsymbol{\rho}})$ that does not depend on ${\boldsymbol{\phi}}$. This reparameterization $\boldsymbol{\zeta}({\boldsymbol{\phi}} , {\boldsymbol{\rho}})$ is appropriately chosen so that one can write $\mathbb{E}_{\boldsymbol{\zeta} \sim q_\phi}[f(\boldsymbol{\zeta})] = \mathbb{E}_{{\boldsymbol{\rho}} \sim p({\boldsymbol{\rho}})}[f( \boldsymbol{\zeta}({\boldsymbol{\phi}}, {\boldsymbol{\rho}}))]$. Therefore, we can move the derivative inside the expectation with no difficulties:

$\begin{eqnarray} &&\partial_{{\boldsymbol{\phi}}} \mathbb{E}_{\boldsymbol{\zeta} \sim q_\phi}\left[f(\boldsymbol{\zeta})\right] = \mathbb{E}_{{\boldsymbol{\rho}} \sim p({\boldsymbol{\rho}})}\left[ \partial_{{\boldsymbol{\phi}}} f( \boldsymbol{\zeta}({\boldsymbol{\phi}}, {\boldsymbol{\rho}})) \right]. \end{eqnarray} \tag{ 8 }$

This is called the reparameterization trick[ 56] and its efficient implementation is responsible for the recent success and proliferation of VAE.

The application of the reparameterization trick as in equation (8) requires that $f( \boldsymbol{\zeta}({\boldsymbol{\phi}}, {\boldsymbol{\rho}}))$ be differentiable, so the latent variables $\boldsymbol{\zeta}$ are continuous. However, discrete latent units can be indispensable to represent the right distributions, such as in attention models, language modeling and reinforcement learning [66, 70, 71]. For example, a latent space composed of discrete variables can learn to disentangle content and style information of images in an unsupervised fashion [72]. Several methods have thus been developed to circumvent the non-differentiability of discrete latent units [69, 73–75]. In the context of VAE, the reparameterization trick has been extended to discrete variables by either relaxation of discrete variables into continuous variables [59, 70, 76] or by introducing smoothing functions [58]. In reference [1], QVAE was introduced based on the implementation of reference [58]. In this work, we follow the implementation of reference [59], which gives biased estimates but provides a much simpler and flexible implementation.

To set up a notation that we keep throughout the paper, we now assume the prior distribution is defined on a set of discrete variables $\mathcal{\boldsymbol{\zeta}} \sim p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}})$, with $\mathcal{\boldsymbol{\zeta}} \in \{0, 1 \}^L$. Given a discrete variable z with mean q and logit $l = \sigma^{-1}(q) = \log(q) - \log(1-q)$ (where $\sigma = 1/[1+\exp(-l)]$ is the sigmoid function), a non-differentiable implementation of equation (8) for discrete variables can be obtained as

$\begin{eqnarray} &&z = \Theta[\rho - (1 - q) ] = \Theta[\sigma^{-1}(\rho) + l )]\,, \end{eqnarray} \tag{ 9 }$

where Θ is the Heaviside function and the random variable ρ ∈ [0, 1] is distributed according to a uniform distribution $\mathcal{U}$. In the second equality, we have used the fact that the inverse sigmoid function is monotonic.

A continuous smoothing (also known as the Gumbel trick [76]), is performed by replacing the Heaviside function with the sigmoid function:

$\begin{eqnarray} &&z = \Theta[\sigma^{-1}(\rho) + l ] \leadsto \boldsymbol{\zeta} = \sigma\left(\frac{\sigma^{-1}(\rho) + l }{\tau} \right)\,, \end{eqnarray} \tag{ 10 }$

where τ is a temperature parameter introduced to control the smoothing. Typically, τ is annealed from large to low values during training. For large values of τ, the bias introduced by substituting $\mathcal{\boldsymbol{\zeta}}$ with $\boldsymbol{\zeta}$ everywhere in the loss function is large, but the gradients propagating through $\boldsymbol{\zeta}$ are also large, facilitating training. Conversely, for low values of τ the bias is reduced but gradients vanish and training stops. Evaluation of trained models is done in the limit $\tau \rightarrow 0$, where $\boldsymbol{\zeta} \rightarrow \mathcal{\boldsymbol{\zeta}}$.

Throughout this paper, we will use BMs to provide powerful and expressive prior distributions defined on discrete variables. For CBMs:

$\begin{eqnarray} p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) & \equiv & {e^{-\mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})}}/Z_{\boldsymbol{\theta}}\,, \quad Z_{\boldsymbol{\theta}} \equiv \sum_{\mathcal{\boldsymbol{\zeta}}} e^{-\mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})} \,, \nonumber \\ \mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) & \equiv & \sum_l z_l b_l + \sum_{l< m} W_{lm} z_l z_m\,. \end{eqnarray} \tag{ 11 }$

To train a VAE with CBM prior, following the prescription of the previous section, we formally replace $p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) \leadsto p_{\boldsymbol{\theta}}(\boldsymbol{\zeta})$. As usual, the gradient of log-probability is given by the difference between a positive and negative phase:

$\begin{equation} \partial \log p_{\boldsymbol{\theta}}({\boldsymbol{\zeta}}_{\boldsymbol{\phi}}) = - \partial \mathcal{H}_{\boldsymbol{\theta}}({\boldsymbol{\zeta}}_{\boldsymbol{\phi}}) + \mathbb{E}_{\bar{\mathcal{\boldsymbol{\zeta}}} \sim p_\theta}[{\partial \mathcal{H}_{\boldsymbol{\theta}}(\bar{\mathcal{\boldsymbol{\zeta}}})}]\,. \end{equation} \tag{ 12 }$

In the equation above, we have highlighted the fact that the smoothed latent samples ${\boldsymbol{\zeta}}_{\boldsymbol{\phi}}$ depend on the variational parameters ${\boldsymbol{\phi}}$. The model samples $\bar{\mathcal{\boldsymbol{\zeta}}}$, however, remain discrete variables sampled from the CBM and are thus not smoothed during training [59].

Once we have a framework to train VAEs with discrete latent variables, we can consider QCH VAEs in which the generative process $\mathcal{\boldsymbol{\zeta}} \sim p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})$ is realized by measuring the computational basis on a given quantum state $\rho_{\boldsymbol{\theta}}$. Such quantum states can be realized by a quantum circuit or via a quantum annealing process controlled by a set of parameters ${\boldsymbol{\theta}}$ we wish to adjust during training of the model.

As introduced in reference [1], QVAE is obtained by assuming the quantum state $\rho_{\boldsymbol{\theta}}$ is a Gibbs state of a transverse field Ising model; i.e., a QBM [52]. The prior $p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})$ distribution is then given by:

$\begin{eqnarray} p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) & \equiv & \text{Tr}[\Lambda_\mathcal{\boldsymbol{\zeta}} {e^{-\mathcal{H}_{\boldsymbol{\theta}}}}]/Z_{\boldsymbol{\theta}}\,, \quad Z_{\boldsymbol{\theta}} \equiv \text{Tr} [e^{-\mathcal{H}_{\boldsymbol{\theta}}}] \,, \nonumber \\ \mathcal{H}_{\boldsymbol{\theta}}& = & \sum_l \sigma^x_l \Gamma_l + \sum_l \sigma^z_l b_l + \sum_{l< m} W_{lm} \sigma^z_l \sigma^z_m, \end{eqnarray} \tag{ 13 }$

where $\boldsymbol{\Gamma}, \mathbf{h}, \mathbf{W} \in \{{\boldsymbol{\theta}} \}$, $\Lambda_\mathcal{\boldsymbol{\zeta}} \equiv | \mathcal{\boldsymbol{\zeta}} \rangle \langle \mathcal{\boldsymbol{\zeta}} |$ is the projector on the classical state $\mathcal{\boldsymbol{\zeta}}$ and $\sigma_l^{x,z}$ are Pauli operators. Unlike for a classical BM, the direct evaluation of the gradients of the $\log p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})$ above is intractable. As discussed in reference [52], a possible workaround is to perform the following substitution:

$\begin{eqnarray} &&p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) = \text{Tr}[\Lambda_\mathcal{\boldsymbol{\zeta}} {e^{-\mathcal{H}_{\boldsymbol{\theta}}}}]/Z_{\boldsymbol{\theta}} \rightarrow \tilde p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) = e^{-\mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})}/Z_{\boldsymbol{\theta}} \end{eqnarray} \tag{ 14 }$

in the ELBO $\mathcal L$, where $\mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})\equiv \langle\mathcal{\boldsymbol{\zeta}}| \mathcal{H}_{\boldsymbol{\theta}} |\mathcal{\boldsymbol{\zeta}}\rangle$, to obtain the so-called quantum ELBO (Q-ELBO) $\tilde{\mathcal L}$. As a consequence of the Golden-Thompson inequality $\text{Tr}[ e^{A} e^{B}] \ge \text{Tr} [e^{A + B}]$, one has:

$\begin{eqnarray} &&p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) \ge \tilde p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) \quad \Rightarrow \quad \mathcal L \ge \tilde{\mathcal L}\,. \end{eqnarray} \tag{ 15 }$

The Q-ELBO $\tilde{\mathcal L}$ is thus a lower bound to the ELBO with tractable gradients that can be used during training. The derivatives of the log-probabilities $\log \tilde p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}})$ can be estimated via sampling from the QBM [52]:

$\begin{equation} \partial \log \tilde p_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) = - \partial \mathcal{H}_{\boldsymbol{\theta}}(\mathcal{\boldsymbol{\zeta}}) + \mathbb{E}_{\bar{\mathcal{\boldsymbol{\zeta}}} \sim p_\theta}[{\partial \mathcal{H}_{\boldsymbol{\theta}}(\bar{\mathcal{\boldsymbol{\zeta}})}}]\,, \end{equation} \tag{ 16 }$

where $\bar{\mathcal{\boldsymbol{\zeta}}}$ are the model samples distributed according to the quantum Boltzmann distribution. The use of the Q-ELBO and its gradients precludes the training of the transverse field $\boldsymbol \Gamma$ [52], which is treated as a constant (hyper-parameter) throughout the training. Training via the Q-ELBO is performed as in the BM case, by smoothing $\mathcal{\boldsymbol{\zeta}} \leadsto \boldsymbol{\zeta}$.

In this section we give more evidence that we have successfully exploited the Chimera-structured QBM prior in the latent space of our convolutional VAE. Validating the training of QCH generative models can be nontrivial and must be assessed carefully, especially when training uses a large amount of classical processing. We trained the deep networks of our model using GTX 1080 Ti GPUs and the quantum annealer is called only to estimate the negative phase. A principled validation strategy is thus to compare a model trained with quantum assistance to a simpler baseline for which the quantum hardware is not required. In our case, a convenient baseline is a model that has the same classical networks but whose prior is trivial. As a trivial prior, we choose a set of independent Bernoulli variables, which is equivalent to a QBM with vanishing weights between latent units. For simplicity, in the following we refer to such prior as Bernoulli prior ⁹ .

For image processing, comparison between different generative models could be done qualitatively by visually inspecting the generated samples [60]. In our research however, visual comparison is inconclusive. In figure 3, for example, we compare samples generated by a trained model with Chimera (figure 3(a)) and Bernoulli (figure 3 b) priors. Both models have been trained by evaluating the negative phase with a D-Wave 2000Q system. The images are also generated using samples coming from the annealers. Given our specific implementation, it is difficult to discern an improvement of visual quality when using a Chimera-structured rather than a Bernoulli prior. To overcome this difficulty, we seek to evaluate quantitatively the generative performance of our model by computing the ELBO defined in equation (3) or by estimating the log-likelihood via an importance sampling technique as described in reference [57]. These quantities are however not accessible when training with analog devices as samplers. In fact, while we assume that samples generated by quantum annealers are distributed according to the required Boltzmann distribution, we must treat quantum annealers as black-box samplers during testing and validation. In other words, the log-probabilities for the quantum generative process $\log\ p^\mathrm{DW}_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}})$ must be assumed to be unknown. Therefore, we validate results by replacing the unknown hardware log-probabilities with those of an auxiliary CBM whose weights are given by the relations in equation (17) ¹⁰ :

$\begin{eqnarray} &&\log p^\mathrm{DW}_{h, J}(\mathcal{\boldsymbol{\zeta}}) \leadsto \log p^\mathrm{CBM}_{b, W}(\mathcal{\boldsymbol{\zeta}})\,. \end{eqnarray} \tag{ 18 }$

This approach can be more rigorously interpreted as validating a fully classical model in which we replace the quantum annealer by the auxiliary CBM defined by the relations above. We can consistently evaluate training using a fully classical auxiliary model since the quantum bias due to the presence of a transverse field is small enough that the estimation of the negative phase of the QBM [80] obtained with the quantum annealer is a good approximation of the negative phase of the auxiliary CBM.

Zoom In Zoom Out Reset image size

In table 1 we report the LL of the auxiliary VAE with 288 latent unit on a Chimera connectivity trained with a D-Wave 2000Q quantum annealer. We compare it to a VAE with a CBM trained end-to-end with PA [1]. We also compare each model with its respective Bernoulli baseline. We have reported the mean and the standard error over 5 independent training runs.

Training on Chimera improves significantly the log-likelihood of the auxiliary model over the Bernoulli baseline. Moreover, the auxiliary models trained with quantum annealers achieved the same log-likelihood as the models trained with PA. Notice that each model and its baseline employed exactly the same amount of classical computational resources. Models trained with structured BMs achieve better performance by requiring thermal sampling, a computational task that can be offloaded to a quantum annealer.

In general we would like to use quantum annealers to sample from the trained generative model. As explained above, treating quantum annealers as black-boxes means we cannot quantitatively evaluate such a model. However, we argue that the log-likelihood of such a QVAE is likely very close to that of the auxiliary VAE. After all, the training assumes the hardware samples are distributed according to a Boltzmann distribution of the auxiliary restricted BM (RBM) and in the previous section we have shown that this assumption is accurate enough to correctly train the auxiliary CBM. We perform an additional visual analysis in figure 4. On the left panel we show a set of digits generated by sampling from a D-Wave 2000Q. On the right panel we use the same trained model but sample from a D-Wave 2000Q after setting its weights to zero. We see that while the annealer still generates plausible digits, it does not generate an ensamble of digits with the correct statistics (in the right panel of figure 4, digits 9 and 4 seem to dominate the scene). This gives evidence that D-Wave 2000Q quantum annealers sampled consistently, such that the classical networks were able to correctly learn the correlations between latent units existing due to the QBM with non-vanishing weights.

Zoom In Zoom Out Reset image size

Exploiting a larger number of latent units to improve the generative performance of a VAE is a popular and active research area. One known obstacle in achieving this is the loss function used for training. We rewrite the ELBO here for convenience by highlighting its two terms:

$\begin{eqnarray} \mathcal L(\mathbf{x}, \mathcal{\boldsymbol{\zeta}},{\boldsymbol{\phi}}) & = & \underbrace{\mathbb{E}_{\mathcal{\boldsymbol{\zeta}}\sim q_{{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}|\mathbf{x})}[\log p_{{\boldsymbol{\theta}}}(\mathbf{x}|\mathcal{\boldsymbol{\zeta}})]}_\mathrm{autoencoding \, term} + \nonumber \\ &- & \underbrace{D_{KL}( q_{{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}|\mathbf{x}) || p_{\mathcal{\boldsymbol{\zeta}}}(\mathcal{\boldsymbol{\zeta}}))}_\mathrm{KL-regularization}\,. \end{eqnarray} \tag{ 19 }$

The first term is sometimes called the 'autoencoding' term and can be thought of as a reconstruction error: an encoded latent configuration $\mathcal{\boldsymbol{\zeta}}eta$ is first sampled from the encoder $q_{{\boldsymbol{\phi}}}(\boldsymbol{\zeta}|\mathbf{x})$ and is subsequently decoded by $p_{{\boldsymbol{\theta}}}(\mathbf{x}|\boldsymbol{\zeta})$. When both approximating posterior and marginal distributions are factorized distributions, this term can be also interpreted as a reconstruction error. Maximizing the ELBO results in maximizing this term, which tends to maximize the number of latent units used to prevent information loss during the encoding-decoding steps. The second term, the KL divergence between the approximating posterior and the prior, has the effect of a regularization term and it is sometimes also called KL regularization. Maximizing the ELBO results in minimizing the KL term, which pushes the approximating posterior close to the prior. This also means the approximating posterior depends less sharply on the inputs $\mathbf{x}$. In the case of factorized distributions, this usually means some latent units are conditionally independent from the input ('inactive units'): $z_{inact} \sim q_{\boldsymbol{\theta}}(z_{inact}|\mathbf{x}) = q_{\boldsymbol{\theta}}(z_{inact})$.

The balance between the autoencoding and KL terms means using VAE can be efficient at lossy data compression by using the right number of latent units. However, the tension between KL and autoencoding terms results in an optimization challenge: during training the model is usually stuck in a local minimum with a suboptimal number of latent units. The main takeaway, for our purposes, is that the number of latent units effectively used is highly dependent on the model, the optimization technique and the training set. To exploit a larger BM, one thus has to work on all these elements. In figure 5 we show the number of active units during a training run of 2000 epochs when the models are trained with either PA or D-Wave 2000Q, with either a Chimera-structured BM or a Bernoulli prior. Lines (bold or dashed) are the means over 5 independent runs while light-color areas delimit the smallest and largest values among the 5 runs. To identify whether a latent unit is active or not, we compute the variance σ of the value of each unit z over the test set and we set a threshold of σ > 0.01 as definition of an active unit.

Zoom In Zoom Out Reset image size

Figure 5 shows several key points that we will expand upon in the next sections. First, it shows the use of a KL annealing technique: the KL term is turned off at the beginning of the training and it is slowly (linearly) turned on within 200 epochs. The figure clearly shows the effect of the KL term in shutting down a large number of active units. Since the number of active units plateaus around a number much lower than 288, using a larger BM usually does not improve performance of our implementation on MNIST. It also shows that connectivity of the BM plays a major role in determining the number of active units, which is much higher with a Chimera-structured BM. Notice also that in the Bernoulli case, both samplers (PA and D-Wave 2000Q) train a model that uses a very similar number of latent units. However, when sampling is nontrivial (as in a Chimera-structured BM) the model trained with the quantum annealer uses a number of units larger than the Bernoulli case, but smaller than the model trained with PA. This is a manifestation (which we will discuss later) of biased sampling with the quantum annealer: sampling quality is good enough to train the model (indeed we obtained a log-likelihood as good as that of the model trained with PA) but exploits a smaller number of latent units.

In the next three sections we discuss three important elements that would allow to build QCH VAEs that can effectively exploit large latent-space BMs. This is a necessary condition to search for quantum advantage in these models: speed-up and scale-up sampling from large BMs using quantum annealers rather than inefficient classical sampling techniques.

5.1.1. Denser connectivities

Connectivity of the BM in the latent space plays an important role in determining both performance of the generative model and number of active latent units. While these two elements are not directly related, we typically observe a correlation between them. In this section, we investigate in more detail the effects of implementing denser connectivities by performing numerical experiments in four different cases: Bernoulli prior and Chimera, Pegasus and fully-connected RBM. Together with Bernoulli and RBM, we pick the connectivities of currently available D-Wave 2000Q quantum annealers and D-Wave's next-generation Pegasus architecture [81] (see also appendix D).

In figure 6 we compare the log-likelihood and the number of active units along training runs obtained using the same classical networks but using a different connectivity for the latent-space RBM prior. At each step during the training run, we employed roughly the same amount of classical resources for back-propagation. Unsurprisingly, using a more capable (dense) BM results in better generative performance (log-likelihood) of the model (left panel)). It also results in using a larger number of latent units (right panel). Notice how Bernoulli and Chimera priors effectively use a number of latent units well below 150, while Pegasus and fully connected use well above 150. Because of this, we could not improve generative performance of Bernoulli and Chimera models by just using larger graphs, at least when training on MNIST. On the other hand, using larger Pegasus and fully connected RBMs would likely have improved the log-likelihood of the model. In figure 6 we show results on the same number of latent units (288) for proper comparison.

Zoom In Zoom Out Reset image size

Working with VAEs allows us to easily take advantage of new connectivities without having to implement a new convolutional VAE. In fact, our model has been fairly optimized to improve performance on Chimera graphs (see the next section). Despite this architecture-specific optimization, just using the denser Pegasus graphs improve performance of the model. Moreover, the flexibility of the VAE hybrid approach allows us to easily adapt the implementation to the slightly different working graphs of different processors with different active/inactive qubits. By using the same implementation, during training the model naturally learns to deactivate latent units corresponding to uncalibrated qubits. We have indeed seamlessly used the same model to train on different D-Wave 2000Q processors, as well as using different groups of qubits within the same processor. In no case was a hard-coded connectivity (which would change for each processor) necessary.

The results of figure 6 show that developing quantum annealers with denser connectivities (such as Pegasus) naturally leads to exploiting larger latent-space BMs, possibly getting us closer to a regime where quantum advantage is possible.

5.1.2. Hardware-specific optimization of classical networks

In our implementation, we have used fairly conventional convolutional NN. As we will discuss in this section, we have implemented only one specific architecture-dependent element in the encoder network that turned out to be very effective at improving performance on the Chimera graph. In general, we believe more elaborate, hardware-specific implementations of both the approximating posterior $q_{{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}|\mathbf{x})$ and the marginal distribution $p_{{\boldsymbol{\theta}}}(\mathbf{x}|\mathcal{\boldsymbol{\zeta}})$ will significantly improve performance of VAE trained with analog devices with quasi two-dimensional connectivities. This is not just a problem of building a model with more capacity. As we discussed in the case of latent-unit use, it is also a problem of optimizing the model hyperparameters. When working with sparsely connected BMs, it is easier for the KL term to push both approximating posterior and BM prior to a local minimum of the loss function in which they are both trivial. Developing hardware-specific hybrid models will thus also aim at reaching local minima during training in which the BM prior is as expressive as possible, exploiting the largest amount of correlations between latent units.

In the context of hybrid VAE models trained with quantum annealers, the considerations above could replace more standard mapping techniques used in the quantum annealing community such as minor embedding [82, 83] and majority voting. The latter techniques typically require a hard-coded specification of the hardware connectivity of each processor, which makes adapting the code to different processors cumbersome. Our perspective in the contest of hybrid generative modeling is to work by adapting classical networks using a high-level specification of the connectivity and letting the model, through stochastic gradient descent, learn the details of the connectivity of each processor (such as the locations of uncalibrated qubits).

We now discuss an example of how the classical networks can be optimized for a given architecture (in our case the Chimera graph). A common technique to implement a more expressive approximating posterior $q_{{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}|\mathbf{x})$ (with a tighter variational bound) is to introduce conditional relationships among latent units, also called hierarchies. In the case of two hierarchies, we first define an approximating posterior for a subset of latent units $\mathcal{\boldsymbol{\zeta}}_1$ and sample from it, then define a second approximating posterior for the remaining latent units $\mathcal{\boldsymbol{\zeta}}_2$ that depends conditionally on both the input data and the sampled values of the first group of latent variables:

$\begin{eqnarray} &&\mathcal{\boldsymbol{\zeta}}_1 \sim q_{1,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_1|\mathbf{x})\,, \quad \mathcal{\boldsymbol{\zeta}}_2 \sim q_{2,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_2,|\boldsymbol{\zeta}_1,\mathbf{x})\,. \end{eqnarray} \tag{ 20 }$

The schematic of our implementation is shown in figure 7(a). We notice that models with a large number of hierarchies (possibly as large as the number of latent units) are possible. Such models, also referred to as autoregressive models [84], are very powerful but have inefficient inference, since sampling must be performed sequentially and cannot be parallelized on modern GPU hardware.

Zoom In Zoom Out Reset image size

The two hierarchical groups $(\mathcal{\boldsymbol{\zeta}}_1, \mathcal{\boldsymbol{\zeta}}_2)$, as well as the physical qubits on the Chimera connectivity, are not equivalent. We can thus build different models by simply choosing the mapping between the two hierarchical groups and an arbitrary bipartition of the physical qubits. Notice that training and deploying each of these models will involve exactly the same amount of classical and quantum computational resources. In our experiments we consider two possible mappings of the two hierarchical groups onto the physical qubits of the Chimera graph, which we call 'Bipartite' and 'Chains'. The Bipartite mapping corresponds to the bipartite structure of the Chimera graph (see figure 7(b). The Chains mapping corresponds to the vertical and horizontal physical layout of qubits of the Chimera architecture (see figure 7(c)). The identification of vertical and horizontal chains of qubits is commonly used to perform a minor embedding of an RBM on the Chimera graph. We stress again, however, that we never perform any minor embedding and we always sample from the native Chimera graph in all cases.

Comparative results of the two mappings, obtained with classical sampling, are shown in figure 8. In the left panel we see that there is a sizeable difference in generative performance between the two mappings, with the Chains mapping performing remarkably better. This better performance is also reflected in the much higher number of latent units exploited by the Chains mapping (see right panel of figure 8). An intuitive explanation of the results above can be given as follows. Let us first write the KL term as:

$\begin{eqnarray} D_{KL}( q_{{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}|\mathbf{x}) || p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}})) & = & D_{KL}( q_{1,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_1|\mathbf{x}) || p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_1)) + \nonumber \\ &+& D_{KL}( q_{2,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1, \mathbf{x}) || p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1))\,. \end{eqnarray} \tag{ 21 }$

In the Bipartite mapping, the conditional $p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1)$ has a simple form that can be computed analytically due to the bipartite structure of the Chimera graph. During training, it is very easy for the model to use the capacity of $q_{2,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1, \mathbf{x})$ to match the simple prior marginal $p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1)$ and be independent from $\mathbf{x}$. As a consequence, a large portion of the representational capacity of the approximating posterior is wasted in representing the simple marginal $p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1)$. In the Chains mapping, in contrast, the marginal $p_{{\boldsymbol{\theta}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1)$ is nontrivial and the approximating posterior $q_{2,{\boldsymbol{\phi}}}(\mathcal{\boldsymbol{\zeta}}_2|\mathcal{\boldsymbol{\zeta}}_1, \mathbf{x})$ has more difficulties in matching it and decoupling $\mathbf{x}$. As a consequence, the model ends up using more efficiently all the variational parameters ${\boldsymbol{\phi}}$. This is another manifestation of the optimization challenge present with VAE models mentioned before, which in this case it is exploited to find better local minima of the loss function.

Zoom In Zoom Out Reset image size

In this section we have shown how a simple architecture-aware modification of the encoder network allows us to train better models and to exploit a given architecture more efficiently. We expect that architecture-aware model-engineering will be crucial to fully exploit large physical connectivities in the latent space of VAE.

5.1.3. Training on larger datasets

Implementing more highly connected BMs and developing classical encoders and decoders tailored to a given connectivity can only go so far in helping to exploit larger latent spaces. Together with other techniques such as KL-term anneal, the ideas mentioned in the previous two sections help reduce the pressure of the KL term to reach suboptimal local minima. In essence, VAE are also efficient lossy encoders. An alternative direction to increase latent space utilization is thus to train on more complex datasets. By doing so, a larger number of latent units is necessary to store enough information such that the reconstruction term (first term in equation (19)) is large enough.

We give numerical evidence of the intuition above by training the same VAE models used in the previous sections on the Fashion MNIST (FMNIST) dataset. A set of images from the FMINST dataset is shown in the left panel of figure 9. FMNIST is the same size as MNIST (60 000, 28 × 28 images) and has the same number of classes. However, its images are more complex with more fine details, including grey-scale features that are important for correct image classification (whereas MNIST digits are substantially black and white). In the right panel of figure 9, we train the same models on MNIST (shaded) and FMNIST (solid) and compare the number of active units during training. We see that, apart from the case with fully connected RBMs, all other models use a substantially larger number of latent units.

Zoom In Zoom Out Reset image size

Exploiting large latent-space BMs is a necessary condition to eventually achieve quantum advantage when sampling with quantum annealers. This condition is however not sufficient. The typical computational bottleneck in training a BM is due to the appearance of well-defined modes. These modes make classical sampling techniques inefficient and slow-mixing, resulting in highly correlated samples used both during training and generation. While making sampling harder, the development of multi-modal distributions is actually an appealing property of BMs, since it allows such models to represent complex and powerful probability distributions. The idea behind searching for quantum advantage in training BM with quantum annealers is, indeed, to exploit quantum resources (such as tunneling) to more efficiently mix between different modes in the landscape defined by the BM.

When trained on visible data, a BM naturally develops complex landscapes to match the complexity of the statistical relationship present in the training data. However, while BMs trained on latent representations can potentially develop well-defined modes, they do not necessarily do so. In fact, one of the capabilities of generative models with latent variables is finding a set of statistically independent latent features [85]. This is typically enforced during model building by using trivial priors such as the product of independent Gaussian (for continuous latent units) or the product of independent Bernoulli (for discrete latent spaces). Even when the prior is potentially complex and trainable, as a BM, the presence of the KL term can push the model during training into local minima in which the trained BM develops a trivial landscape.

In this section we give evidence that BMs trained in the latent space of a VAE model do indeed develop a nontrivial landscape with well-defined modes. As we have shown in section 5.1.1 (see figure 6), BMs with denser connectivities naturally lead to better performing VAEs. This is an indirect indication that we are indeed exploiting the additional capacity and expressivity of more connected BMs. In this section we give more explicit evidence of this. In figure 10, we generate a sequence of images via block Gibbs sampling. The top left image is generated by picking a latent configuration $\mathcal{\boldsymbol{\zeta}}$ out of a uniform distribution over all configurations. This latent sample is then sent through the decoder. Going from left-to-right, top-to-bottom, each subsequent image is obtained after updating all latent units with a sequence of block Gibbs updates (one for Bernoulli, two for the bipartite connectivities Chimera and RBM, four for the quadripartite Pegasus connectivity). As expected, in the Bernoulli case (figure 10(a)), each update results in uncorrelated samples. The Chimera connectivity (figure 10(b)) is able to develop weakly correlated samples, as shown by short sequences of similar images. Correlated samples with well-defined modes are more clearly visible with the Pegasus connectivity ((figure 10(c))). Finally, we confirm that increasing the connectivity up to a fully connected RBM (figure 10(d)) results in long sequences of correlated samples and related to the deep valleys of the RBM energy landscape.

Zoom In Zoom Out Reset image size

The results shown in figure 10 show that BMs trained as priors of generative models with latent variables naturally learn multi-modal, nontrivial probability distributions. These distributions are expressive, making the whole VAE more expressive, while at the same time developing the same types of computational bottlenecks that make classical sampling algorithms inefficient. This paves the way to effectively use quantum annealers as means to more scalable quantum-assisted sampling, enabling us to sample from BMs of sizes and complexity that are infeasible with classical methods.

Using quantum annealers to train large BMs directly on complex data remains challenging. Apart from the unsatisfying performance of using BMs with quasi two-dimensional connectivities on visible data, a major difficulty is biased sampling (and thus inaccurate gradients) obtained with quantum annealers. There are two main sources of bias: control errors and imperfect or incomplete thermalization at the freezing point. While the latter can be improved with appropriate pause-and-ramp annealing schedules, the former can only be improved with technological advancements. Despite these known difficulties, we have shown in the previous sections we have successfully trained large BM (hundreds of units) in the latent space of a VAE solely using samples coming from a D-Wave 2000Q, without using any hard-coded pre or post-processing to the raw data obtained from the annealer.

We interpret our positive results as an indication that training BMs with our setup is relatively robust to noise and control errors. In fact, we can interpret both the encoder and decoder as powerful tools to pre- and post-process data to be sent to the quantum annealer. Using stochastic gradient descent, we train the encoder and decoder to generate a set of latent features that are more easily modeled by the latent-space BM. For example, real images might have strongly correlated, sharp features (such as regions with black or white pixels), which require large weights to be modeled correctly. A precise implementation of such large weights might be challenging for analog devices with finite range such as quantum annealers. Additionally, both encoders and decoders might be able to learn and correct, or at least reduce the effects of, systematically biased sampling.

To investigate the role of noise and control errors in determining sampling quality and performance of the trained models, we perform a set of comparative experiments in which we train the same model on MNIST dataset, using samples coming from three D-Wave 2000Q with different noise profiles. The Baseline and Lower-Noise D-Wave 2000Q are both publicly available on D-Wave's Leap^TM cloud service. We have also included an Interim Lower-Noise processor with an intermediate noise profile that is internally available at D-Wave.

Results are shown in figure 11. In figure 11(a) we report the log-likelihood during training. Models trained on the Lower and Interim Lower-Noise D-Wave 2000Q achieved performance comparable to training with PA. The evaluation of the log-likelihood for the Baseline D-Wave 2000Q diverged due to diverging weights, as can be seen in figure 11(b). The weights of the Baseline processor start diverging after about 100 epochs. The Interim Lower-Noise processor shows an opposite behavior, with weights getting small and plateauing after about 500 epochs. For the Lower-Noise processor, the L1 of the weights plateaus at a value that is closer to that obtained with 'noiseless' PA. Notice that a consistent comparison in figure 11(b) we have reported the weights W rescaled by the effective temperature (see equation (17) and not the 'bare' annealing values J. Figure 11(b) highlights the remarkably different response of three different quantum annealers to our model. Despite such differences, the performance of our hybrid implementation is robust (as shown in figure 11(a)) and does not require any hardware-specific adaptations or fine-tuning. Only while training with the Baseline D-Wave 2000Q, we needed a more aggressive clipping (that is restricting the weights and biases to have narrower range than the maximum allowed) to achieve similar performance and converged log-likelihood estimation.

Zoom In Zoom Out Reset image size

Noise and control errors also manifest in a less efficient use of the latent space, as seen in figure 11(c). All models trained with D-Wave 2000Q use fewer active units than PA. Since the estimate of the log-likelihood of the models trained with the Baseline processor did not converge, we focus on the comparison between the Interim and Lower Noise processor: the latter can exploit a larger number of latent units. We finally show in figure 11(d) the profile for the extracted effective temperature during training, which is remarkably different for the three D-Wave 2000Q. The results shown in figure 11(d) underlines the importance of developing advanced annealing techniques to stabilize temperature during training.

In this section we have demonstrated the robustness to noise of our implementation, as highlighted in figure 11(a). At the same time, figure 11(c) shows an important effect of noise, which is to make it harder for our hybrid model to exploit the optimal number of latent units. We thus anticipate that exploiting large BMs (with thousands of units, eventually) following the directions indicated in the previous sections must be accompanied by continued efforts in reducing sampling bias due to noise and control errors of future-generation quantum annealing devices.