Training quantum Boltzmann machines with the β-variational quantum eigensolver

The QBM [6, 7] is defined as the n-qubit Gibbs state (mixed density matrix)

$$\begin{align} \sigma_w = \frac{e^{-\beta H_w}}{Z}, \end{align} \tag{ 1 }$$

with $Z = \textrm{Tr}(e^{-\beta H_w})$ the partition function, β > 0 the fixed inverse temperature and H_w a parameterized Hamiltonian. In this work, we consider Hamiltonians of the form

$$\begin{align} H_w = \sum_{r} w_r H_r, \end{align} \tag{ 2 }$$

which consists of a sum of a polynomial number of non-commuting qubit-operators H_r with weights w_r . The specific form of the operators H_r does not matter at this point and will be given in the result section below. The quantum Hamiltonian H_w generalizes the role of the energy function in a classical Boltzmann machine [20].

In order to find the optimal weights $w^*$, we aim to minimize the quantum relative entropy

$$\begin{align} S\left(\eta \| \sigma_w\right) = \textrm{Tr}\left(\eta \log{\eta}\right) - \textrm{Tr}\left(\eta \log{\sigma_w}\right) \end{align} \tag{ 3 }$$

from σ_w to the target density matrix η [8, 9]. This information-theoretic measure generalizes the classical Kullback–Leibler divergence to density matrices and has nice properties for the QBM model ansatz, e.g. the logarithm cancels the exponential in equation (1). The target η is either some quantum system with unknown Hamiltonian H_t at finite temperature, i.e. $\eta = \frac{e^{-\beta_t H_t}}{Z_t}$, or a classical data set embedded into a quantum state.

The optimization of $S(\eta \|\sigma_w)$ can be done with gradient descent. By noticing that only the second term in equation (3) depends on w_r , and using Duhamel's formula for the derivative of a matrix exponential [21], we obtain the partial derivatives

$$\begin{align} \frac{\partial S}{\partial w_r} = \textrm{tr}\left(H_r \eta\right) - \textrm{tr}\left(H_r \sigma_w\right) . \end{align} \tag{ 4 }$$

Similar to the classical case [20], the gradient is given by the difference between statistics under the data density matrix η and the model density matrix σ_w . The data statistics do not change during training and need to be obtained only once. However, the QBM statistics change after every iteration of gradient descent and must be recomputed. The QBM density matrix σ_w in equation (1) is given by the matrix exponential of the Hamiltonian, whose dimension scales exponentially in the number of qubits. This means that computing the full density matrix becomes computationally intractable for larger system sizes.

While the classical BM has the same exponential scaling of the number of states, there exist several approximation algorithms that can effectively compute the necessary statistics, often based on Markov-chain Monte Carlo [22]. However, for quantum systems, no general MCMC methods exist due to the negative sign problem for non-stoquastic models [23, 24]. Since the QBM is not restricted to the stoquastic domain and can include spin-glass models, quantum Monte Carlo is not an option in general. In the following section, we review the β-VQE method that with help of a quantum computer can be used to potentially circumvent this issue.

The QBM can be trained on classical and quantum data. Since classical data is not naturally represented in density matrices, the data is embedded in a quantum system. There are various ways to do this. We use the pure-state embedding presented in [9]. From a classical data set $\mathcal{D}$, the empirical probability distribution is constructed counting the relative occurrence of each binary vector in the data, $q(\boldsymbol{s}^{i}) = \frac{1}{|\mathcal{D}|} \sum_{\boldsymbol{s}^j \in \mathcal{D}} \delta_{\boldsymbol{s}^{i},\boldsymbol{s}^j}$. We create a pure quantum state representation from the empirical distribution by mapping the n-dimensional binary vector s ⁱ to a spin-$\frac{1}{2}$ system with n particles. Each 1 in s ⁱ corresponds with a particle in the spin-up state $\vert 1\rangle$ and each 0 to a particle in spin-down state $\vert 0\rangle$ in the σ^z basis. Let:

$$\begin{align} \vert \boldsymbol{s}^i\rangle = \bigotimes_{j = 0}^{n} \vert s^i_j\rangle. \end{align} \tag{ 5 }$$

The quantum state is obtained by taking the superposition of all spin states corresponding to the binary vectors with coefficients equal to $\sqrt{q(\boldsymbol{s})}$, i.e. $\vert \psi\rangle = \sum_{\boldsymbol{s} \in \mathcal{D}} \sqrt{q(\textbf{s})} \vert \boldsymbol{s}\rangle$. The density matrix embedding is then constructed by taking the outer product of this state,

$$\begin{align} \eta = |{\psi}\rangle\langle{\psi}| = \sum_{\boldsymbol{s}, \boldsymbol{s}^{\prime}} \sqrt{q\left(\mathbf{\boldsymbol{s}}\right)q\left(\mathbf{\boldsymbol{s}^{\prime}}\right)}|{\boldsymbol{s}}\rangle\langle{\boldsymbol{s}^{\prime}}|.\end{align} \tag{ 6 }$$

Recently, Liu et al [12] proposed a variational method to represent a mixed density matrix called β-VQE. It combines a classical network with a quantum circuit (figure 1) into a parameterized density matrix ansatz given by

$$\begin{align} \rho_{\theta, \phi} = \sum_{\mathbf{s}}p_{\phi}\left(\mathbf{s}\right) U_{\mathbf{\theta}} |{\mathbf{s}}\rangle\langle{\mathbf{s}}|U^{\dagger}_{\mathbf{\theta}}. \end{align} \tag{ 7 }$$

Here $p_{\phi}(\mathbf{s})$ is a parameterized probability distribution over binary states implemented by a classical network, e.g. a parameterized Bernoulli distribution or an autoregressive network [25]. The binary states are then used as inputs for a parameterized quantum circuit U_θ , which transforms them into quantum states. The combined system $\rho_{\theta, \phi}$ can accurately represent finite-temperature Gibbs states of quantum spin systems, e.g. see [12]. Here we aim to see if it can be used for the wide range of Hamiltonians needed to train the QBM on quantum and classical data.

The optimal parameters of $\rho_{\theta, \phi}$ for a fixed QBM Hamiltonian H_w are found by minimizing the variational free energy

$$\begin{align} \mathcal{F}\left(\sigma_w, \rho_{\theta, \phi}\right) = \textrm{Tr}\left(\rho_{\theta, \phi} \log{\rho_{\theta, \phi}}\right) + \textrm{Tr}\left(\rho_{\theta, \phi} H_w\right). \end{align} \tag{ 8 }$$

This is equivalent to minimizing the quantum relative entropy from σ_w to $\rho_{\theta, \phi}$. The minimization can be done with gradient descent on the classical network and quantum circuit parameters simultaneously. The required gradients [12] are given by

$$\begin{align} \nabla_{\theta} \mathcal{F} = \sum_{\mathbf{s}} p_{\phi}\left(\mathbf{s}\right) \left[ \nabla_{\theta} \langle \mathbf{s}|{U^{\dagger}_{\theta} H_{w} U_{\theta}}|{\mathbf{s}} \rangle \right], \end{align} \tag{ 9 }$$

for the circuit, and

$$\begin{align} \nabla_{\phi} \mathcal{F} & = \sum_{\mathbf{s}} p_{\phi}\left(\mathbf{s}\right) \left[ \left(f\left(\mathbf{s}\right) - b\right) \nabla_{\phi} \log{p_{\phi}\left(\mathbf{s}\right)} \right], \end{align} \tag{ 10 }$$

for the classical network. Here we used a control variate to reduce the variance as in [12]: $f(\mathbf{s}) = \ln{p_{\phi}(\mathbf{s})} + \langle\mathrm{s}|{U_\theta^\dagger H_w U_\theta}|{\mathbf{s}}\rangle$ and $b = \sum_{\mathbf{s}} p_{\phi}(\mathbf{s}) f(\mathbf{s})_{p_{\phi}}$. For quantum gates of the form $R_G(\theta_j) = e^{- i \theta_j G /2}$ where G is a Pauli operator, the gradient w.r.t. the circuit parameters, equation (9), reduces to a parameter shift rule: a weighted sum of energy expectation values [26]. Suppose the circuit is a product of L such gates, i.e. $U_{\theta} = R_{G_L}(\theta_L) R_{G_{(L-1)}}(\theta_{L-1}) \cdots R_{G_1}(\theta_1)$. Then

$$\begin{align} \frac{\partial}{\partial \theta_i} \langle\mathrm{s}|{U^{\dagger}_{\theta} H_{w} U_{\theta}}|{\mathbf{s}}\rangle & = \frac{1}{2} \left[ \langle{\psi_{i-1}}|{R^{\dagger}_{G_i}\left(\theta_i + \tfrac{\pi}{2}\right) \tilde{H}_w R_{G_i}\left(\theta_i + \tfrac{\pi}{2}\right)}|{\psi_{i-1}}\rangle\right.\nonumber\\ & \quad \left. - \langle{\psi_{i-1}}|{R^{\dagger}_{G_i}\left(\theta_i + \tfrac{\pi}{2}\right) \tilde{H}_w R_{G_i}\left(\theta_i + \tfrac{\pi}{2}\right)}|{\psi_{i-1}}\rangle \right] , \end{align} \tag{ 11 }$$

where for convenience, $\vert \psi_{i-1}\rangle = \prod_{j = 1}^{i-1} R_{G_j}(\theta_j) \vert \mathbf{s}\rangle$ and $\tilde{H}_w = \prod_{j = (i+1)}^{L} R^{\dagger}_{G_j}(\theta_j) H_w \prod_{j = L}^{i+1} R_{G_j}(\theta_j)$.

Computation of the gradient is tractable on a quantum computer if the sum is approximated by sampling a polynomial number of bitstrings from p_φ (s), and if there is no barren plateau in the optimization landscape [19]. The gradient w.r.t. the parameters of the classical network, equation (10), has the form of a standard weighted entropy gradient, which for classical networks can be efficiently approximated from sampling [27].

We propose a slight modification of the sampling approach for the gradients originally presented in [12]. Instead of randomly sampling from all computational basis states, we choose the R states with the highest probability and renormalize p_φ on those states. This results in a β-VQE density matrix of rank R

$$\begin{align} \rho_{\theta, \phi} = \frac{1}{\sum_{k = i_1}^{i_R} p_{\phi}\left(\boldsymbol{s}_k\right)} \sum_{j = i_1}^{i_{R}} p_{\phi}\left(\boldsymbol{s}_{j}\right) U_\theta| \boldsymbol{s}_j\rangle\langle{\boldsymbol{s}_j}| U_\theta^\dagger. \end{align} \tag{ 12 }$$

In the following, we refer to this as truncated-rank β-VQE. As an immediate consequence, the gradients for the variational free energy equations (9) and (10) are also truncated for this model. This means that we can heuristically choose a small R so that the optimization can be performed at a reduced computational cost.

We now describe the nested-loop algorithm to train a QBM with a truncated rank β-VQE. Recall that the aim is to train the QBM so that it resembles some target data embedded in η. This is schematically shown in figure 1. The algorithm starts with a simple ansatz for the QBM Hamiltonian H_w , e.g. a Heisenberg XXZ model or a random spin-glass model. In the inner loop, the β-VQE ansatz $\rho_{\theta, \phi}$ is trained to represent the QBM σ_w by minimizing equation (8). In the outer loop, $\rho_{\theta, \phi}$ is used to compute approximate QBM statistics $\textrm{tr}( H_r \sigma_w) \approx \textrm{tr}(H_r \rho_{\theta^*,\phi^*})$. These are used for training the QBM by gradient descent with equation (4).

Since H_w changes only slightly between two QBM iterations, we use a warm-start strategy for the parameters $\theta, \phi$ of the β-VQE ansatz $\rho_{\theta, \phi}$ between the steps of the outer loop. We reuse the converged parameters of $\rho_{\theta, \phi}$ in one β-VQE inner loop as the initial parameters for the next β-VQE inner loop. Hence, while the β-VQE has to run for a relatively long time in the first inner loop, we find that it converges much faster in subsequent iterations.

The total number of statistic estimations of the nested-loop algorithm is of order

$$\begin{align} \mathcal{O}\left( n_{\text{QBM}} \times n_{\beta\text{-VQE}} \times R \times n_{\theta} \right). \end{align} \tag{ 13 }$$

Here n_QBM is the number of QBM iterations for convergence to some fixed accuracy, $n_{\beta\text{-VQE}}$ the number of β-VQE iterations, R the rank of the β-VQE, and n_θ the number of circuit parameters.

The classical data we consider in this paper is introduced in [30]. A patch of 160 connected ganglion cells from the salamander retina is connected to a multi-electrode array, where the activity of multiple neurons is recorded simultaneously. The ganglion cells are then exposed to 297 repeats of a movie. Every timeframe of 20 ms, the activity for each neuron is recorded and denoted by a 1 (spike) or a 0 (no spike), resulting in a binary vector of size 160 for each time frame. Each iteration of the movie consists of 953 timeframes.

We selected subsets of the ganglion cells with high mutual information (MI). This selection procedure is described in detail in appendix A. To obtain a representative training set, the data is divided into iterations of the movie. The first 30 iterations are selected for the training set $\mathcal{D}_{\text{train}}$, resulting in a set of $28\,590$ binary vectors. Other iterations of the movie are grouped per 10 to form test sets, resulting in 26 separate test sets. Consequently, all timeframes occur in an equal ratio within the training set and the test sets. Since many test sets are available, this provides an excellent method for testing how well the QBM generalizes to unseen data.

In figure 3, we show the convergence of the quantum relative entropy during training. We observe that $S(\eta\|\sigma_w)$ decreases monotonically and saturates at S = 0.12. This is close to the value obtained from learning with the exact gradient (dashed line). To compare these results to classical models, the QBM probability distribution over classical states can be obtained through projection, i.e. $p_{\text{QBM}}(\boldsymbol{s}^i) = \textrm{Tr}(|\boldsymbol{s}^i\rangle\langle{\boldsymbol{s}^i}| \sigma_w)$. We can then directly compare the classical relative entropy, or KL-divergence, to the classical BM. The QBM achieved a KL-divergence of 0.010 on the training data. In contrast, the classical BM only achieves a KL-divergence of 0.022 when fully converged. These results generalize well to unseen data. The QBM outperforms the classical BM, achieving a lower KL-divergence on all test sets. A histogram of the results is shown in appendix A.

Zoom In Zoom Out Reset image size

Note that the norm of the gradient, shown in the inset of figure 3, is $\mathcal{O}(10^{-4})$, which is still significantly nonzero. In contrast, the classical BM can be trained until the norm of the gradient reaches machine precision. The convergence of the QBM was numerically studied in [9] as a function of the rank of the target state η where it was shown that convergences deteriorates with lower rank. It has been shown that later training iterations tend to increase the magnitude of the weights, akin to lowering the effective temperature of the QBM. During this phase, the ground state of the QBM, $\psi_{w}^{0}$, already accurately represents the target state despite the nonzero gradient, achieving better results than the classical BM. This is also reflected by the high fidelity of F = 0.998 between the ground state and the target.

We investigate some important aspects of the learning algorithm, starting with the rank used in the truncated-rank β-VQE inner loop. In figure 4 we show the fidelity for different states obtained at the end of our algorithm. The fidelity between the QBM and η, indicated by the blue crosses, increases substantially with increased β-VQE rank, even though the target is a pure state. This might seem counter-intuitive but is a direct result of the low-rank approximation of the gradient. For a very low rank β-VQE, the approximation to the gradient becomes vanishingly small once the ground state accurately approximates the target pure state. During that stage, the exact σ_w is still very mixed, but the gradient is approximated as if it were close to a pure state. Conversely, when a higher rank approximation is used, the gradient is still nonzero even if the ground state of σ_w closely resembles the target pure state, allowing the QBM to converge further towards a pure state. The fidelity between the target and the ground state of the QBM is depicted in figure 4 by the yellow triangles. High fidelity of F > 0.99 is achieved with as few as two samples.

Zoom In Zoom Out Reset image size

Another important aspect is the scaling of the performance of the algorithm with the number of features (qubits) n. In order to investigate this, we train the QBM on subsets of the salamander retina data of different sizes. We set the circuit depth equal to the number of features, which was heuristically found to be a good rule of thumb, see appendix B. Motivated by the previous sampling analysis, we also scale the rank of β-VQE linearly with n. In figure 5 we show the fidelity obtained with these heuristic guidelines. The QBM trained by the nested-loop algorithm is able to obtain a fidelity of F > 0.98 for all system sizes studied here. Note that the fidelity decreases for larger systems. This is because of a model mismatch between the QBM and the target, as this drop is observed also in the exact QBM. A more expressive Hamiltonian ansatz, e.g. with the inclusion of third-order interaction terms, might yield more accurate results for larger system sizes.

Zoom In Zoom Out Reset image size

In our second example, we train a QBM to model quantum data. As target state η we take the finite temperature Gibbs state $\eta = e^{-\beta H_{\mathrm{XXZ}}}/Z$ with

$$\begin{align} H_{\mathrm{XXZ}} = \sum_{i = 1}^{n-1}J\left(\sigma^{x}_i\sigma^{x}_{i+1}+\sigma^y_i\sigma^y_{i+1}\right)+\Delta \sigma_i^z\sigma^z_{i+1}, \end{align} \tag{ 16 }$$

with $J = -1$, $\Delta = -\frac{1}{2}$, and β = 1. Quantum data taken from systems at finite temperatures is usually characterized by a high rank. However, we show that the β-VQE can be used to get an efficient low-rank approximation to the target state. We train the QBM with the nested-loop algorithm, again analyzing the effect of the rank. The results for the quantum data are shown in figure 6. In contrast to the results for classical data shown in figure 4, we require a higher rank approximation to obtain a high fidelity. By increasing the rank, the β-VQE is more capable of representing the mixed QBM state. In the truncated-rank β-VQE, one has the freedom to increase the rank freely during training, at the cost of a higher computational demand. For mixed target states, it can be efficient to start out with a low-rank approximation for early iterations of the outer loop, when only a rough approximation of the QBM gradient is sufficient for training, and increase the rank as the gradient approaches zero.

Zoom In Zoom Out Reset image size

All the simulations in previous sections assume access to the exact quantum state. This is possible via classical statevector simulations which are intractable in general. To overcome this, one can execute the required circuits on a quantum computer. This introduces two new sources of error. First, the quantum computer is affected by various sources of hardware noise [31, 32]. While little information is available on the noise for a specific quantum chip, we believe this noise is biased in general. Thus, this is a source of error that cannot be completely mitigated by taking many measurements. Second, the output statistics are approximated from a finite number of shots (measurements). This is important, for example, when estimating the gradient of β-VQE in equation (9). When hardware noise is absent, the error due to a finite number of shots (M) scales as $\epsilon \propto \frac{1}{\sqrt{M}}$. In this section, we take a closer look at these sources of errors in the context of QBM training.

We use Qiskit [33] to simulate both the noise-free and noisy hardware. For the noisy hardware simulations, we use the fake provider module from Qiskit to mimic the behavior of IBM's quantum devices. We use a subset of the salamander retina data of size n = 4 and the Hamiltonian ansatz given in equation (14). For β-VQE we truncate the rank to R = 2 and use the circuit in figure 2 with a depth of d = 2.

Figure 7 compares the QBM training in three different settings in terms of estimating the expectation values: using exact values, using a noise-free quantum simulator with 5000 shots, and finally using a noisy quantum hardware simulator with 5000 shots. The FakeKolkata backend from Qiskit, which mimics the behavior of IBM's ibmq_kolkata device, is used to produce the noisy hardware results.

Zoom In Zoom Out Reset image size

As expected, the learning stops earlier in the presence of noise. An increase in the noise causes the errors introduced by the β-VQE algorithm to reach the size of the QBM gradient faster. The size of the error compared to the size of the gradient is shown in figure 8 for the noisy and noise-free quantum simulators, when both are using the same number of shots to approximate the expectation values. As expected, a higher (lower) quantum relative entropy (fidelity) is achieved using a noisy hardware compared to a noise-free one.

Zoom In Zoom Out Reset image size

Finally, we deploy the QBM model trained using the FakeKolkata backend on the corresponding real quantum hardware device, ibmq_kolkata. We sample the thermal state in the computational basis and compare its KL-divergence to the target distribution. The target distribution is resampled a number of times in order to obtain different test sets. The results are shown in figure 9. As expected, the ibmq_kolkata (red bars) yields KL-divergences higher than those of the noise-free quantum simulator (blue bar). However, deploying the same QBM model on another IBM device, ibmq_perth (yellow bars), results in a KL-divergence about twice as high as using ibmq_kolkata, demonstrating the recent improvements made to noise reduction in quantum hardware. We expect that further advancements in noise reduction will contribute significantly to improving these results. Note that QBM models trained in the noisy setting are to some extent bound to the specific hardware [7] (data not shown). Further investigation is needed to examine how hardware noise affects QBM training by increasing the number of qubits.

Zoom In Zoom Out Reset image size