Faster variational quantum algorithms with quantum kernel-based surrogate models

The two quantum kernels that we will consider in this work are the 'state kernel', which is equivalent to the fidelity between quantum states, and the 'unitary kernel', which quantifies the similarity between two unitary operations in terms of their Hilbert-Schmidt norm. For an ansatz circuit $\vert \psi(\boldsymbol{\theta}) \rangle = U(\boldsymbol{\theta})\vert \boldsymbol{0} \rangle$, with an n-qubit input state $\vert \boldsymbol{0} \rangle = \vert 0 \rangle^{\otimes n}$, the state kernel evaluated between gate angles $\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{^{\prime}}$ is defined as:

$\begin{equation} \begin{split} k_s(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}})&: = \left\vert \langle \psi(\boldsymbol{\theta}^{^{\prime}}) \vert \psi(\boldsymbol{\theta}) \rangle \right\vert^2\; = \left\vert \langle \boldsymbol{0} \right\vert U(\boldsymbol{\theta}^{^{\prime}})U(\boldsymbol{\theta}) \left\vert \boldsymbol{0} \rangle \right\vert^2 \end{split} \end{equation} \tag{ 2 }$

or equivalently by $k_s(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}): = \textrm{Tr}\{\rho(\boldsymbol{\theta}^{^{\prime}})\rho(\boldsymbol{\theta})\}$ (assuming a pure quantum state). As we discuss in section 3.2 and explicitly show in appendix A, the energy function is linear in the feature space induced by the state kernel provided this kernel is based on the same circuit used to calculate $E(\boldsymbol{\theta})$. As a result the state kernel is still expected to be effective when building general surrogate models for $E(\boldsymbol{\theta})$.

The unitary kernel also induces a feature space in which the $E(\boldsymbol{\theta})$ is a linear function. For an ansatz imparting a d × d-dimensional unitary $U(\boldsymbol{\theta})$ the unitary kernel is defined as

$\begin{equation} k_u(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}): = \left\vert \frac{\textrm{Tr}\{U\,^\dagger(\boldsymbol{\theta}^{^{\prime}})U(\boldsymbol{\theta})\}}{d} \right\vert^2, \end{equation} \tag{ 3 }$

which is the (normalized) modulus-squared Hilbert-Schmidt norm between the unitaries $U(\boldsymbol{\theta})$ and $U(\boldsymbol{\theta}^{^{\prime}})$. We show in appendices B and C that while $E(\boldsymbol{\theta})$ is linear in the feature space induced by this kernel this space typically has a much larger dimension than that of the state kernel and so more data is required to build accurate regression models using it.

We will now explicitly construct the feature maps of these two kernels for an ansatz circuit $\vert \psi(\boldsymbol{\theta}) \rangle = U(\boldsymbol{\theta})\vert \boldsymbol{0} \rangle$. Note that while a kernel function and corresponding reproducing kernel Hilbert space are unique, the feature map is not unique and is only defined up to an isometry (which preserves the inner product and thereby the kernel). We consider n-qubit p-parameter ansatzes formed of an initial state $\vert \boldsymbol{0} \rangle = \vert 0 \rangle^{\otimes n}$ acted upon by p parameterized Pauli rotations (PPRs), with n-qubit Pauli operators $\{P_1,\ldots,P_p\}$, which are interleaved with p + 1 arbitrary fixed unitaries $\{R_1,\ldots R_{p+1}\}$. The corresponding unitary $U(\boldsymbol{\theta})$ is given by:

$\begin{equation} U(\boldsymbol{\theta}): = R_{p+1}e^{-iP_p\theta_p/2}R_p\ldots R_2e^{-iP_1\theta_1/2}R_1. \end{equation} \tag{ 4 }$

'Hardware efficient' ansatzes [26] admit this description and are formed of individual PPRs with low-weight Pauli operators; these ansatzes are common in NISQ quantum computing as they often map efficiently onto a device's native gate set [16]. Ansatzes of the form given by (4) are extremely flexible, being able to implement arbitrary parameterized unitaries provided a one uses a suitably many PPRs, a sufficiently sophisticated encoding of the angles $\boldsymbol{\theta}$, and a careful choice of the $\{R_i\}$ (although these may require exponentially deep circuits [61]).

We show in appendix A that the state kernel's feature space can be constructed in terms of an ansatz-dependent $3^p$-element vector of operators on $\mathcal{H}\otimes\mathcal{H}^*$ (where $\mathcal{H}$ is the qubits' Hilbert space) $\boldsymbol{s}$, given by

$\begin{equation} \boldsymbol{s} = \frac{1}{2^p}\left[ \begin{pmatrix} (I \otimes I + P_p \otimes P_p^*)R_p\otimes R_p^*\\ (I\otimes iP_p^* -iP_p\otimes I)R_p\otimes R_p^*\\ (I \otimes I - P_p \otimes P_p^*) R_p\otimes R_p^* \end{pmatrix} \otimes_K\ldots\otimes_K \begin{pmatrix} (I \otimes I + P_1\otimes P_1^*)R_1\otimes R_1^*\\ (I\otimes iP_1^* -iP_1\otimes I)R_1\otimes R_1^*\\ (I\otimes I - P_1 \otimes P_1^*) R_1\otimes R_1^* \end{pmatrix} \right] \end{equation} \tag{ 5 }$

and an angle-dependent vector $\boldsymbol{v}(\boldsymbol{\theta})$, given by

$\begin{equation} \boldsymbol{v}(\boldsymbol{\theta}) = \begin{pmatrix}1\\\sin(\theta_p)\\\cos(\theta_p)\end{pmatrix}\otimes_K\ldots\otimes_K\begin{pmatrix}1\\\sin(\theta_1)\\\cos(\theta_1)\end{pmatrix}. \end{equation} \tag{ 6 }$

Here we have made a distinction between the tensor products ⊗ between the physical and conjugated Hilbert spaces of the qubits and Kronecker products $\otimes_K$ which are used to construct the kernels' feature maps. The $\otimes_K$ act as tensor products between the sub-vectors in $\boldsymbol{s}$ and $\boldsymbol{v}$p imposing matrix multiplication in the former case and ordinary multiplication in the latter. Note that the elements of $\boldsymbol{v}(\boldsymbol{\theta})$ are linearly independent Fourier components and so cannot be decomposed further.

The state kernel is then

$\begin{equation} k_s(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}) = \boldsymbol{v}^T(\boldsymbol{\theta}^{^{\prime}})\boldsymbol{S} \boldsymbol{v}(\boldsymbol{\theta}), \end{equation} \tag{ 7 }$

where $(\boldsymbol{S})_{ij} = \left\langle\!\left\langle\right.\right.\!\!\rho_0 \rvert s_i^\dagger s_j\lvert \rho_0 \left\rangle\!\left\rangle\!\!\right.\right.$ is a Hermitian matrix of inner products between the (unnormalized) states $s_i\lvert \rho_0 \left\rangle\!\left\rangle\!\!\right.\right.$ and $s_j\lvert \rho_0 \left\rangle\!\left\rangle\!\!\right.\right.$ ($s_i$ and $s_j$ are the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ components of $\boldsymbol s$- , respectively), $\rho_0 = \vert \boldsymbol{0} \rangle\langle \boldsymbol{0} \vert$, and $\lvert X \left\rangle\!\left\rangle\!\right.\right.$ denotes the vectorization of X so that $\lvert \rho_0 \left\rangle\!\left\rangle\right.\right. = \vert \boldsymbol{0} \rangle\otimes\vert \boldsymbol{0}^* \rangle$ [62]. As ${\boldsymbol S}$ is a Hermitian matrix of inner products it is a positive-semidefinite Gram matrix and admits a Cholesky decomposition $\boldsymbol{S} = Q^\dagger Q$ (where $Q$ is an upper-triangular matrix). This means that the state kernel can be written as the inner-product of two ${3^p}$-element vectors, $k_s(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}) = \boldsymbol{v}^T(\boldsymbol{\theta}^{^{\prime}})Q^\dagger Q \boldsymbol{v}(\boldsymbol{\theta}) = \boldsymbol{\varphi}_s^\dagger(\boldsymbol{\theta}^{^{\prime}}) \boldsymbol{\varphi}_s(\boldsymbol{\theta})$, and so its feature map $\boldsymbol{\varphi}_s:\mathcal{X}\to\mathcal{F}_s$ is

$\begin{equation} \boldsymbol{\varphi}_s(\boldsymbol{\theta}) = Q\boldsymbol{v}(\boldsymbol{\theta}),\ \textrm{where}\ Q^\dagger Q = \boldsymbol{S}. \end{equation} \tag{ 8 }$

The construction of the feature map for the unitary kernel is similar and involves the same two vectors $\boldsymbol{v}(\boldsymbol{\theta})$ (6) and $\boldsymbol s$ (5). The unitary kernel is given by

$\begin{equation} k_u(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}) = \boldsymbol{v}^T(\boldsymbol{\theta}^{^{\prime}})\boldsymbol{T} \boldsymbol{v}(\boldsymbol{\theta}), \end{equation} \tag{ 9 }$

where, $(\boldsymbol{T})_{ij} = {\textrm{Tr}\{s_i^\dagger s_j\}}/{4^n}$ is a Hermitian Gram matrix of Hilbert-Schmidt inner products between the operators $s_i$ and $s_j$. It too is positive-semidefinite allowing it to be written in Cholesky form $\boldsymbol{T} = B^\dagger B$ (with $B$ upper-triangular), meaning that the unitary kernel is $k_u(\boldsymbol{\theta},\boldsymbol{\theta}^{^{\prime}}) = \boldsymbol{v}^T(\boldsymbol{\theta}^{^{\prime}})B^\dagger B \boldsymbol{v}(\boldsymbol{\theta}) = \boldsymbol{\varphi}_u^\dagger (\boldsymbol{\theta}^{^{\prime}})\boldsymbol{\varphi}_u (\boldsymbol{\theta})$ and its feature map $\boldsymbol{\varphi}_u:\mathcal{X}\to\mathcal{F}_u$ is

$\begin{equation} \boldsymbol{\varphi}_u(\boldsymbol{\theta}) = B\boldsymbol{v}(\boldsymbol{\theta}),\ \textrm{where}\ B^\dagger B = \boldsymbol{T}. \end{equation} \tag{ 10 }$

So far we have assumed that all the gate angles used in the PPRs are independent. In more sophisticated ansatzes whose gate angles are linear combinations of the ansatz parameters one can apply linear transformations to $\boldsymbol{v}(\boldsymbol{\theta})$ such that it contains a set of independent Fourier components and so apply a similar analysis. Provided the couplings between gate parameters are linear (for example if multiple parameterized gates share the same parameter), this will act to reduce the dimension of the induced feature space.

For Bayesian VQE to be effective, the kernel used must result in an GP model which can well-approximate $\tilde{E}(\boldsymbol{\theta})$, the noisy energy function. An interesting question this poses is whether quantum kernels are more ideally suited to this task than typical classical kernels. To justify that these kernels are well suited, we can analyze noiseless energy landscape $E(\boldsymbol{\theta})$ for a Hamiltonian H and an ansatz unitary of the form given in (4). In appendix A we show that this is

$\begin{equation} E(\boldsymbol{\theta}) = \boldsymbol{h}^T\boldsymbol{v}(\boldsymbol{\theta}), \end{equation} \tag{ 11 }$

where $h_i = \left\langle\!\left\langle\right.\right.\!\!R_{p+1}H R_{p+1}^\dagger \rvert s_i\lvert \rho_0 \left\rangle\!\left\rangle\right.\right.$ and s is given in (5). Encouragingly, we see that the angle-dependent part of the quantum kernels' feature spaces $\boldsymbol{v}(\boldsymbol{\theta})$ makes an appearance, mirroring the decomposition of the energy function in [51] in terms of a weighted sum of Fourier components. For a suitable set of energy observations, e.g. at the points $\{(\theta_1,\cdots,\theta_p)\in\{0,\pi/2,\pi\}^{\times p}\}$ or $\{(\theta_1,\cdots,\theta_p)\in\{-2\pi/3,0,2\pi/3\}^{\times p}\}$ it is straightforward to see how the components of h (the weights of the Fourier components in $E(\boldsymbol{\theta})$) can be directly inferred (the latter set of points giving exactly the result in [51]). As the feature spaces of both the state and unitary kernels are related to $\boldsymbol{v}(\theta)$ by linear transformations, the energy is a linear function in both these spaces. More explicitly, $E(\boldsymbol{\boldsymbol{\theta}}) = \boldsymbol{h}^T Q^{-1}\boldsymbol{\varphi}_s(\boldsymbol{\theta}) = \boldsymbol{h}^T B^{-1}\boldsymbol{\varphi}_u(\boldsymbol{\theta})$ (Q and B are defined in (8) and (10), respectively and are assumed to be invertible).

Given a set of training data (pairs of ansatz parameters and corresponding noiseless energy evaluations), the representer theorem [57] guarantees that the minimizer of any regularized empirical loss functional is given by a finite linear sum of weighted kernel evaluations evaluated at the training points (ansatz parameters whose energy values have been observed) [57]. For an l₂-regularized least-squares loss the minimizer is kernel ridge regression [63], which is mathematically equivalent to the posterior mean of an GP (21). Because these kernels have finite dimensional feature spaces, kernel ridge regression or GP models using them should be able to achieve perfect global prediction accuracy (can perfectly predict $E(\boldsymbol{\theta})$ for all θ ) if the size of the training data set is greater than or equal to the feature space dimension (assuming the data is linearly independent in the feature space).

We note that linearity of $E(\boldsymbol{\theta})$ in these kernels' induced feature spaces is achieved without any hyperparameters. This removes the need for any maximum likelihood hyperparameter optimization. It also means that the state or unitary kernel between any two observed points (used to construct the Gram matrix K) only needs to be evaluated a single time.

While both quantum kernels look promising for VQE, one might suspect the state kernel to be more suitable as VQE directly concerns the energy of the state produced by the ansatz, rather than its unitary.

To evaluate the effectiveness of these quantum kernels versus typical classical kernels in kernel-based VQE regression we first consider the problem of attempting to build an GP model for observations of the noiseless $E(\boldsymbol{\theta})$. The Hamiltonian we will use throughout is the anti-ferromagnetic 1D transverse field Ising model with a longitudinal field, given by:

$\begin{equation} H = J\sum_{\langle i,j\rangle}Z_iZ_j +h_x\sum_iX_i+h_z\sum_iZ_i, \end{equation} \tag{ 12 }$

where J is a coupling strength, h_x is the transverse field strength, h_z is the longitudinal field strength, X_i and Z_i are the Pauli X and Z operators acting respectively on the ith qubit, $Z_iZ_j$ denotes simultaneous Pauli Z operators on the ith and jth qubits, and the first sum is taken over nearest-neighbor pairs in a 1D spin chain with periodic boundary conditions. Throughout we use $J = h_z = 0.5$ and $h_x = -0.5$. As this is a real-valued Hamiltonian, its eigenstates can also be taken to be real. This means that for VQE we can limit our ansatz circuit to only yielding real wavefunctions. The ansatz we will use is shown in figure 1 and is designed to be efficiently implementable on a quantum processor with nearest-neighbor connectivity between qubits (which helpfully mirrors the nearest-neighbor terms in the Hamiltonian). It consists of parameterized rotations $RY(\theta_j) = e^{-i\sigma_y\theta_j/2}$ and fixed controlled-X gates to build up entanglement.

Zoom In Zoom Out Reset image size

Let $\mathbf{G}[\mu(\cdot),k(\cdot,\cdot, \boldsymbol{\alpha})]$ be an GP model with covariance/kernel function k, kernel hyperparameters α , and mean function µ. Given a set of observed training and validation energies, $\boldsymbol{y}_{\mathrm{t}} = \{E(\boldsymbol{\theta}_{\mathrm{t}_1}),\ldots,E(\boldsymbol{\theta}_{\mathrm{t}_{\mathrm{N_t}}})\}$ and $\boldsymbol{y}_{\mathrm{v}} = \{E(\boldsymbol{\theta}_{\mathrm{v}_1}),\ldots,E(\boldsymbol{\theta}_{\mathrm{v}_{\mathrm{N_v}}})\}$ at parameters $\mathbf{X}_\mathrm{t} = \{\boldsymbol{\theta}_{\mathrm{t}_1},\ldots,\boldsymbol{\theta}_{\mathrm{t}_{\mathrm{N_t}}}\}$ and $\mathbf{X}_\mathrm{v} = \{\boldsymbol{\theta}_{\mathrm{v}_1},\ldots,\boldsymbol{\theta}_{\mathrm{v}_{\mathrm{N_v}}}\}$ respectively, we define validation score $R^2_\mathrm{v}$ as

$\begin{equation} R_\mathrm{v}^2(\mathbf{X}_\mathrm{v},\boldsymbol{y}_\mathrm{v}|\mathbf{X}_\mathrm{t},\boldsymbol{y}_\mathrm{t}): = 1-{\frac{\sum_{\boldsymbol{\theta}\in \mathbf{X}_\mathrm{v}}[E(\boldsymbol{\theta})-\hat{y}_{\mathbf{G}}(\boldsymbol{\theta})]^2}{\mathrm{Var}[\boldsymbol{\,y}_\mathrm{v}]}}, \end{equation} \tag{ 13 }$

where $\hat{\,y}_{\mathbf{G}}(\boldsymbol{\theta}) = \mathbb{E}[\,y^*|\boldsymbol{\theta},\mathbf{X}_t,\boldsymbol{\,y}_t]$ is the prediction of model G at point θ trained with data $(\mathbf{X}_\mathrm{t},\boldsymbol{y}_{\mathrm{t}})$ (which, if appropriate, includes hyperparameter optimization). A validation score of $R^2_{\mathrm{v}} = 0$ indicates that the GP model is performing as well (with respect to the L₂ loss) as a model that always predicts $\mathbb{E}[\boldsymbol{\,y}_\mathrm{v}]$, while a score of $R^2_{\mathrm{v}} = 1$ shows perfect prediction across the unseen validation set.

Figure 2 shows the improvement in validation score as the size of the training data is increased for GP regression models equipped with different quantum and typical classical kernels. The models are presented with noiseless evaluations of $E(\boldsymbol{\theta})$ for the Hamiltonian given in (12) ($J = h_z = 0.5$ and $h_x = -0.5$). The training and validation sets ($\mathbf{X}_{\mathrm{t}}$ and $\mathbf{X}_{\mathrm{v}}$) are drawn uniformly at random from the interval $\{-\pi,\pi\}$ (for each of the 16 parameters). Both the median $R^2_{\mathrm{v}}$ over 100 repeats (solid line) and the inter-quartile range of the data are shown. The quantum kernel and energy evaluations were performed noiselessly on a classical computer. In all cases, a small diagonal offset was added to the Gram matrices used in GP model prediction $\mathbf{K}\to\mathbf{K}+10^{-10} I$ (equivalent to $\sigma_n^2 = 10^{-10}$ in (21)) to ensure numerical stability.

Zoom In Zoom Out Reset image size

We see that an GP model equipped with the state kernel greatly outperforms the other kernels and reaches the optimal validation score of $R^2_{\mathrm{v}} = 1$ once the size of the training set reaches $N_{\mathrm{t}} = 136$. The validation score of the unitary kernel improves much more slowly, and for $N_{\mathrm{t}} = 136$ achieves an $R^2_{\mathrm{v}}\lt0.1$. However this is still much higher than the scores for the various classical kernels which perform poorly at this task. These classical kernels yield validation scores $R^2_\mathrm{v}\lt0$ and show negligible signs of improvement as the training data set grows.

It is clear, and perhaps unsurprising, that a state kernel built from the same circuit used for a VQE problem provides a much more accurate similarity measure between data points than the other kernels we consider. This is because the energy function is linear in the finite-dimensional feature space induced by the state kernel (as discussed in section 3.2). While $E(\boldsymbol{\theta})$ is also linear in the unitary kernel's feature space, we show in section 6 that this kernel's induced feature space is typically quadratically larger than that of the state kernel. To build a globally accurate GP model, the observations must span a significant portion of this space which explains why k_s vastly outperforms k_u at this task. The classical kernels we use are universal but do not appear particularly well-suited to the complicated periodic and highly oscillatory $E(\boldsymbol{\theta})$. As this is 16-parameter problem, the density of the couple of hundred observations in the input parameter space will be far lower than what is required for these classical kernels to make good use of their universal property.

For BO we are not always interested in creating an accurate global (i.e. $\theta_i\in[-\pi,\pi]$) regression model as in figure 2. While good global accuracy is helpful, in BO only promising regions of parameter space (those that look close to the desired extremum) are explored in detail. On smaller parameter scales the energy landscape will be less dramatically oscillatory and the density of observed data will be higher meaning that the classical kernels may have a better chance of accurately interpolating between observations.

To test this we examine the validation score for GP models when performing local regression of a VQE energy landscape, restricted to a smaller region of the parameter space. In this case, the training set $\mathbf{X}_\mathrm{t}$ and validation set $\mathbf{X}_\mathrm{v}$ are generated by first picking a common anchor point $\boldsymbol{\theta}^{(0)}$ whose elements are chosen uniformly at random in the range $[-\pi,\pi]$. The circuit parameters for $\mathbf{X}_{\mathrm{t}}$ and $\mathbf{X}_{\mathrm{v}}$ are then obtained by sampling parameter vectors θ with elements drawn uniformly at random from the interval $\theta_i\in[\theta^{(0)}_i-\pi/s,\theta^{(0)}_i+\pi/s]$, where s is a scale reduction factor.

Figure 3 illustrates this process, showing the log prediction error $\log_{10}(1-R^2_{\mathrm{v}})$ for a repeat of the simulations in figure 2 on reduced parameter scales with scale factors s = 10 (3(a)) and s = 100 (3(b)). We see that when required to make only local predictions, GP models using the unitary and classical kernels show significant improvement although those using the state kernel are consistently the most accurate. This is to be expected as observations are much denser in the input space. Denser observations allow the universal classical kernel-based models to more effectively interpolate between observations. They also reduce the portion of the state and unitary kernel's induced feature spaces which the observations must span for $k_s$- and $k_u$-based GP models to perform well.

Zoom In Zoom Out Reset image size

The state kernel's advantage is most pronounced for s = 10, which is to be expected as this is closer to global regression. The unitary, RQ, RBF kernels and the Matern kernel with $\nu = 5/2$ all perform similarly well but the Matern kernel with $\nu = 3/2$ struggles when presented with large amounts of data (particularly for s = 10). This is likely because a GP equipped with a Matern kernel is $\left\lceil \nu \right\rceil-1$ times differentiable (once differentiable for $\nu = 3/2$ and twice for $\nu = 5/2$) [63], while those using the unitary, RBF, and RQ kernels are infinitely differentiable. As the noiseless energy function (a finite weighted sum of Fourier components, see 3.2) is smooth and infinitely differentiable this may explain why the singly differentiable $\nu = 3/2$ Matern kernel performs worse than the other, more-times differentiable kernels.

For it to be useful, the predictions from a surrogate model should be significantly easier to compute and optimize than the cost function. Given $m$ observations, making a prediction from an GP surrogate (see section 8.1), as is done when optimizing a GP-based surrogate model in BO, costs $m$ kernel evaluations, while fitting the model requires a Gram-matrix of $\mathcal{O}(m^2)$ kernel evaluations to be calculated.

This raises a potential issue with quantum kernel-based surrogate models, namely that both fitting the surrogate model and using it to make predictions requires repeated executions of a large number of quantum circuits. Unless the number of data points is smaller than the number of circuits required to measure the cost/energy function one could simply measure and optimize the cost directly, avoiding any inaccuracies of the approximate surrogate model. Accordingly, quantum kernels evaluated on-device are unlikely to be useful for producing the kinds of easy-to-optimize surrogate models needed for BO.

It is also worth noting that unless the ansatz used is highly structured, both the state and unitary kernels evaluated for two randomly chosen sets of gate angles will shrink exponentially with the number of qubits used [64]. While this is less of an issue for small NISQ-scale VQE problems, this means that for problems on many qubits one would require exponentially-many shots to accurately resolve these exponentially-small kernel evaluations. This is closely related to the barren-plateau problem and further reinforces our assertion that quantum kernels evaluated on-device are unlikely to be useful for producing globally-accurate surrogate models for VQE.

Device noise also greatly complicates the evaluation of quantum kernels on-device. Although using a kernel based on the noisy operation of a quantum processor may be useful in quantifying the similarity between noisy energy evaluations such kernel functions may require large numbers of samples to evaluate accurately. This would greatly complicate any surrogate model optimization and means that the usual requirements of a kernel, for example positivity, would not necessarily be satisfied. It would also be difficult to ensure that the errors present when evaluating the kernel correspond to those encountered when estimating the cost/energy. For example, one typical implementation of the state kernel requires a circuits of twice the depth used to estimate the energy [65] and so involves more noise and decoherence than would be seen in the cost function. Alternatively, one can implement the kernel using a circuit of the same depth as the cost function but over two sets of qubits [66] which will likely have different noise characteristics.

While GPs naturally accommodate noisy observations it is normally assumed that the kernel function can be computed exactly. Circuit noise and finite sampling errors would mean an on-device quantum kernel evaluation is affected by statistical fluctuations. This would then affect the two key objects in GP modeling; the positive-definite Gram matrix $\mathbf K$ of pair-wise kernel evaluations for the observed points and the vector $\boldsymbol k$ of kernel evaluations between a point of interest and the observed points (see section 8.1 for a full discussion of these). Randomness in these objects can be detrimental in two ways. Firstly, if the Gram matrix ${\mathbf K}$ is composed of noisy kernel evaluations then it may not be positive definite and numerical instabilities can occur when calculating its inverse (needed to make predictions). This can be partially alleviated by increasing the noise strength hyperparameter σ_n (see section 8.1) to ensure an invertible and positive-definite $\mathbf K$ .

Another more serious issue comes when attempting to minimize a surrogate model built using a noisy kernel as part of a BO loop. Gradient descent methods are the standard approach for this however they require repeated accurate evaluations of $\nabla_{\boldsymbol{\theta}} \boldsymbol{k}$, the gradient of kernel evaluations between the query point and the training data. While noise-resilient methods such as SPSA noise-resilient could be used for this optimization, these could be applied directly to the objective function, bypassing the need for a surrogate model. A possible exception would be if the surrogate model is considerably cheaper to evaluate than the cost function, which would be true if the number of observations is smaller than the number of circuits required to evaluate the cost.

Due to these complications, we do not expect significant practical advantages to VQE from GPs which use device-evaluated quantum kernels. Instead we evaluate all the quantum kernels using a classical computer, which can be done tractably provided the number of qubits is relatively small and depth of the quantum circuits involved is relatively low.

For many-qubit VQE problems, in which most kernel evaluations are exponentially suppressed, one could still attempt a highly localized form of Bayesian VQE with classically-evaluated quantum kernels. In such scenarios, the Gram matrix $\mathbf K$ calculated for a set of randomly chosen initial points will be exponentially close to the identity. Similarly, kernel evaluations between a new point of interest $\boldsymbol{\theta}^*$ and these initially chosen points will be close to zero unless $\boldsymbol{\theta}^*$ is in the immediate vicinity of an initial point. Bayesian VQE using an GP with such an exponentially-decaying kernel will be highly localized around whichever initial point gives the lowest energy. This means that at each step in the BO, the maximum of the acquisition function (see section 8.2) will stay close to the most promising initial point (where the kernel is non-negligible). This optimization strategy would still differ from completely local methods like as gradient descent in that all previous observations (however close to each other they may be) are used to decide the next query point rather than just the most recent point seen. However for this to work, one would still need to find a useful starting point for the optimization (i.e. not in a barren-plateau) and have a method for evaluating the kernel. This could either be estimated on-device, would likely be more costly than optimizing the energy directly, or evaluated classically, which would be intractable unless the ansatz has some simple exploitable structure or admits a simplifying approximation (we discuss this in section 5).

We have seen that the state kernel provides a significant advantage over typical classical kernels when performing GP regression of a circuit's energy landscape and now apply these results to VQE using BO. BO is a gradient-free strategy for optimizing noisy expensive-to-evaluate objective functions [56]. As current cloud-based NISQ computers are in high demand, the problem of variationally minimizing the noisy $\tilde{E}(\boldsymbol{\theta})$ is an apt use-case for BO. We describe BO in more detail in section 8.2.

We will quantify the performance of the optimization in terms of the best seen energy error $\mathcal{E}(\boldsymbol{y},\boldsymbol{\theta}_{\mathrm{opt}})$ which we define as the fractional difference between the lowest energy seen at the current stage in the optimization and the minimum achievable noiseless energy for the ansatz. For a set of observed energy values $\boldsymbol y$ and an optimal set of gate angles $\boldsymbol{\theta}_{\mathrm{opt}}$, as defined in (1), this is given by

$\begin{equation} \mathcal{E}(\boldsymbol{y},\boldsymbol{\theta}_{\mathrm{opt}}) : = \frac{\min(\boldsymbol{y}) - E(\boldsymbol{\theta}_\mathrm{opt})}{\vert E(\boldsymbol{\theta}_\mathrm{opt}) \vert}. \end{equation} \tag{ 14 }$

The optimal gate angles $\boldsymbol{\theta}_{\mathrm{opt}}$ for the ansatz shown in figure 1 were found using $10\,000$ attempts of direct gradient-based minimization, giving a minimum energy of $E(\boldsymbol{\theta}_{\mathrm{opt}}) = -2.762194$.

Figure 4 shows the results of noiseless VQE simulations using BO with different kernels. The ansatz circuit used is shown in figure 1 and the Hamiltonian is given in (12) with $J = h_z = 0.5$ and $h_x = -0.5$. At the start of each simulation, 25 initial points were drawn uniformly at random and their energies evaluated to form an initial training data set. The Expected Improvement [67] acquisition function was used with an exploration hyperparameter $\xi$ = 0.01, this value is widely used in EI-based BO [68, 69] and was found to yield the best BO performance across the classical kernels considered. For each kernel, 100 repeats of the optimization were performed. Figure 4(a) shows the median (solid lines), mean (dashed lines), and inter-quartile (filled) range of the energy error. Figure 4(b) gives histograms and Gaussian kernel density estimates (with bandwidth 0.15 [70]) of the final best seen energy errors (after 80 new points have been requested).

Zoom In Zoom Out Reset image size

The state kernel outperforms the other kernels at this VQE task, both in terms of the final energy found and in the speed at which it converges. It frequently reaches an energy less than one part in $10^{-3.5}\approx 0.03\%$ away from the minimum achievable with this ansatz. Its final achieved energies form clusters around a few values, the most noticeable being one at errors between 10⁻⁴ and 10^−3.5 and one that is close to, but still lower than, the majority of energies observed with the unitary and classical kernels. The unitary kernel also performs well, quickly and consistently converging to an error of approximately $4\%$ and exhibits the second lowest mean and median errors overall. However, the data is strongly clustered with this error, implying that achieving a more accurate final energy may be difficult. The classical kernels vary considerably in their performance, generally having much worse mean errors than the quantum kernels. The final energies achieved by the classical kernels form two clusters; one of relatively successful runs with errors ∼6% and one with errors between 30%–60%. The errors in this latter cluster are similar to those seen in the randomly selected initialization data, implying that these optimization runs failed to achieve any significant reduction in the energy. This suggests that while these GP models using these classical kernels sometimes lead to relatively successful VQE, they often immediately and become stuck in an exploitative phase where points close to those in the initial training data are repeatedly queried. It is therefore likely that as well as their advantages in over-all performance, Bayesian VQE using the state or unitary kernel can be more resilient against initialization failures than when using the classical kernels we have considered.

For comparison, figure 5 shows a repeat of this noiseless VQE using the SPSA optimization scheme. We used the Qiskit implementation [71] which closely follows the initial proposal in [72]. We see from the experimental traces 5(a) and histograms in 5(b) that the effectiveness of this optimization strategy can vary but the average final energies typically have errors in line with those obtained for BO using the unitary kernel at around $4\%$. These are also comparable to the best experimental runs for BO using the various classical kernels. The data is strongly clustered around this energy error and only a handful of SPSA runs manage to compete with the low energies frequently achieved with state kernel-based BO. The SPSA runs which cluster around $4\%$ generally require several hundred energy evaluations to achieve this level of accuracy whereas only around 100 evaluations are needed for the various kernel based models. The two SPSA runs with the lowest mean final energies had mean final energy errors on the order of $0.1\%$ which is comparable to the lowest energies achieved with state kernel-based BO. However, these were only seen after 1500 total energy evaluations to be made. This level of final accuracy with SPSA appears to be extremely rare and requires greater than an order of magnitude more energy evaluations than is needed with the state kernel, highlighting the effectiveness and economy of our optimization strategy.

Zoom In Zoom Out Reset image size

The question remains whether a GP surrogate using a noiseless classically-evaluated quantum kernel is useful for Bayesian VQE of noisy quantum circuits. The noisy energy function is unlikely to be a linear function in the feature spaces induced by the noiseless state and unitary kernels. However, provided the device noise is not too great, these noiseless classical evaluations should provide a good approximation to the correct feature space for describing the noisy energy. By using a quantum kernel function based on the ansatz circuit we are able to leverage our prior knowledge of the (noiseless) ansatz circuit whereas standard GP surrogates can only attempt to do this through hyperparameter optimization. Because the observed energy values provide us with information about the noise on the device we are also able to take this noise into account implicitly.

Figure 6 shows the results of a repeat of the Bayesian VQE simulations illustrated in figure 4 with noisy energy evaluations. These were performed with using Qiskit's noisy quantum circuit simulation framework [71] and a noise model derived from the errors present on ibmq_quito in Oct 2021 [73]. The noise includes contributions from gate errors, state preparation and measurement errors (for which Qiskit's readout error mitigation was used), and finite circuit shots (10 000 per circuit per evaluation). The quantum kernel evaluations (where used) were simulated exactly using a tensor network quantum circuit implementation based on the Quimb Python package [74]. As in figure 4, figure 6(a) shows the median (solid lines), mean (dashed lines), and inter-quartile range (filled) of the best seen energy error; the difference between the lowest seen energy and the minimum possible for the ansatz. Figure 6(b) shows histograms and kernel density estimates (with bandwidth 0.15) of the final best seen energy errors (after 80 new points are queried).

Zoom In Zoom Out Reset image size

The noise in the simulations makes it unlikely that the energy can reach the noiseless minimum value. Accordingly, we also show the mean energy error obtained from $10\,000$ repeated noisy evaluations of the energy $\tilde{E}(\boldsymbol{\theta}_{\mathrm{opt}})$, taken at $\boldsymbol{\theta}_{\mathrm{opt}}$, the noiseless optimum gate angles. While $\boldsymbol{\theta}_{\mathrm{opt}}$ may differ slightly from the gate angles that yield the true noisy optimum these evaluations serve as an indicator of the best performance to be expected from a noisy VQE implementation.

We again see the GP models based on the state kernel greatly outperform those using the other kernels. The state kernel frequently manages to achieve energies lower than the mean $\tilde{E}(\boldsymbol{\theta}_{\mathrm{opt}})$ (an error of ∼2%) while also showing the fastest convergence. The unitary kernel also performs well, consistently giving energies lower or at least as low as those found using classical kernels (around $6\%$) but with a significantly lower variation in its performance. When dealing with noisy evaluations, VQE using the classical kernels shows qualitatively similar performance to the noiseless case (although with different final energy errors). Almost all the simulations performed with the RBF kernel and approximately half of those with the other classical kernels give a final errors in the range 30%–60%. As in the noiseless case, these errors are similar to those seen for the randomly-selected initialization points implying that these runs failed almost immediately. The remainder of the classical kernel simulations achieve a much better performance which is close to the ∼6% seen with the unitary kernel.

To test our approximation strategy we compare the predictive accuracy of a GP model built with approximated state kernels of various $\chi_{\mathrm{max}}$ to the accuracy when using the full state kernel. We again consider the Hamiltonian in (12) with $J = h_z = 0.5$ and $h_x = -0.5$. To ensure that the full state kernel is tractable for large numbers of observations we use a brickwork ansatz with the same structure and open boundary conditions as in (1) but with 6 qubits and 20 layers of CX gates on alternating pairs of adjacent qubits separated by RY gates (a total of 50 CX gates and 106 RYs). To engineer a situation in which a low $\chi_{\mathrm{max}}$ MPS approximation is unlikely accurate we choose the parameterized portion of the ansatz, U_B , as the last two layers of RYs (10 of which surround the final CX layer) and pick a fixed set of random gate angles for the remaining gates (which form U_A ). For each simulation we generate two sets of random uniform gate angles $\mathbf{X}_{\mathrm{t}}$ and $\mathbf{X}_{\mathrm{v}}$ and their corresponding noiseless energies $\boldsymbol{y}_{\mathrm{t}}$ and $\boldsymbol{y}_{\mathrm{v}}$ as training and validation data. The size of the training data set is varied while the size of the validation set is fixed to $N_{\mathrm{v}} = 100$.

Figure 7(a) shows the validation scores achieved by GP surrogates equipped with MPS approximated state kernels of different $\chi_{\mathrm{max}}$ compared to one equipped with the full state kernel. The fidelities of the MPS approximations with different maximum bond dimensions to $\vert \psi_A \rangle$ are shown in table 1. To help assess the relative suitability of the different GP models in explaining their training data, figure 7(b) shows the Log Bayes factor between the GP models with different approximated kernels and a GP model using the full state kernel. The Bayes factor, $\mathcal{B}(\mathbf{X},\boldsymbol{y},A,B) = {p(\boldsymbol{y}|\mathbf{X},A)}/{p(\boldsymbol{y}|\mathbf{X},B)}$, [83] between two statistical models A and B is defined as the ratio of the marginal likelihoods of a given set of training data under the two models. Its logarithm is

$\begin{equation} \log \mathcal{B}(\mathbf{X},\boldsymbol{y},A,B) = M_A(\mathbf{X},\boldsymbol{y}) - M_B(\mathbf{X},\boldsymbol{y}), \end{equation} \tag{ 17 }$

where $M_A(\mathbf{X},\boldsymbol{y})$ and $M_B(\mathbf{X},\boldsymbol{y})$ are the log marginal likelihoods of the data $(\mathbf{X},\boldsymbol{y})$ with the models A and B respectively (which for the GP models we use is given by (23)). If the models A and B are assigned equal prior probabilities then the Bayes factor is equivalent to the ratio of the posterior probabilities of the two models $p(A|\boldsymbol{y},\mathbf{X})/p(B|\boldsymbol{y},\mathbf{X})$. This gives the Bayes factor a useful interpretation as how much more plausible one model is at explaining the data than the other [83]. A log Bayes factor $\mathcal{B}(\mathbf{X},\boldsymbol{y},A,B)\gt1$ generally implies that model A is more strongly supported by the data than model B [83].

Zoom In Zoom Out Reset image size

Table 1. Fidelities of MPS approximations with $\vert \psi_A \rangle$ used in simulations. For the maximum bond dimensions allowed in the MPS approximation, the fidelity $\vert \langle \psi_A \vert \tilde{\psi}_A^{(\chi_{\mathrm{max}})} \rangle \vert^2$ with the unapproximated $\vert \psi_A \rangle$ is shown.

$\chi_{\mathrm{max}}$	1	2	3	4	5	6	7
Fidelity	0.0054	0.2904	0.5788	0.8341	0.9382	0.9716	0.9976

We see that MPS approximations to $\vert \psi_A \rangle$ with a higher $\chi_{\mathrm{max}}$ have a larger fidelity with $\vert \psi_A \rangle$ state and their GP models show increasingly similar performance to the GP which uses the full state kernel. Interestingly, the GP models that use approximated state kernels with lower $\chi_{\mathrm{max}}$ often have a higher validation score and Bayes factor (relative to the full state kernel) $\gg$1, provided the size of the training data set is small. When dealing with small numbers of observations these models have both a higher predictive accuracy for unseen energies and are better supported by their training data than the model using $k_s$. This may be because the induced feature spaces of these kernels have a smaller dimension than the full state kernel, being based on a restricted subspace of the n qubit Hilbert space. As a result, less data would be needed to span an appreciable portion of the feature space and so make accurate predictions.

For all the approximated kernels there comes a point where the validation score begins to decrease as more training data is presented. This happens when the size of the training data set becomes close to the dimension of the approximated kernel's feature space. As the energy function will not be completely linear in the approximated kernel's feature space, it becomes increasingly difficult to reconcile these different observations well with a linear model. This feature space saturation will also occur for the GP based on the full state kernel when the validation score reaches $R^2_{\mathrm{v}} = 1$ and the data spans the whole feature space. However, this is not an issue because the energy function is truly linear in the full state kernel's feature space and so energy observations can always be fully reconciled by a linear model in this space. Linear dependence will occur when the number of observations exceeds the feature space dimension at which point the Gram matrix $\mathbf{K}$ also becomes rank deficient; even before this point small eigenvalues will be present in $\mathbf{K}$ causing the GP model's predictions to become numerically unstable.

Typically a small regularizing offset, e.g. $\mathbf{K}\to \mathbf{K}+10^{-10}\mathbf{I}$, is added to Gram matrices to ensure they are invertible by effectively adding additional dimensions to the feature vectors. Provided the offset is small, the regularization has little effect on models with well-suited kernels, like the full state kernel, as the objective function has minimal dependence on the additional components of the feature vectors introduced by the regularization. It can however help smooth out apparent inconsistencies encountered when using a kernel whose feature space is improperly aligned with the objective function by increasing the small eigenvalues present in $\mathbf{K}$ and reducing their disproportionate contribution to the predictions.

Choosing the value for this offset is key to ensuring accurate and numerically stable predictions; it must be sufficiently large to damp the contributions of small eigenvalues of $\mathbf{K}$ but not so large that it dominates the model. We see in figure 7(b) that when the feature space saturation occurs, degrading the validation score, the likelihood of the training data also drops precipitously. As this can be calculated from the training data alone, we can use this as a metric for calibrating the diagonal offset to avoid numerical instabilities. We do this by treating the weighting of the diagonal offset as a kernel hyperparameter $\sigma_n^2$ such that $\mathbf{K}\to\mathbf{K}+\sigma_n^2 I$ and optimizing this to maximize the marginal likelihood of the observations. Adding this hyperparameter is equivalent to assuming that we are dealing with observations corrupted by additive Gaussian noise with zero mean and variance $\sigma_n^2$; discrepancies between the observations and the assumption that they are drawn from a linear function in the kernel's feature space are effectively treated as an additional source of noise. When dealing with real noisy data, adding this noise term to avoid feature space saturation issues introduces negligible additional cost as a this term can be used to simultaneously deal with the observation noise.

Figure 8(a) shows the results of a repeat of the simulations in 7(a) when the approximated kernels are equipped with this additional noise term. For comparison with the classical kernels used in our other simulations, the plot in figure 8(b) shows the validation scores obtained for GP models using these classical kernels as well as the unitary kernel. Because the $U_B$ used in these simulations only has a depth of 1 (a single layer of entangling gates) the unitary kernel can be calculated classically at low cost and so provides a benchmark against which the approximated state kernels should be compared.

Zoom In Zoom Out Reset image size

When given an additional noise term, the approximated state kernels no longer suffer from an eventual drop in predictive accuracy due to feature space saturation. Instead, their performance continues to improve as more data is provided until the validation score eventually begins to saturate. For $\chi_{\mathrm{max}}\gt1$ the approximated state kernels all yield better performance than any of the classical kernels and the unitary kernel. Although the unitary kernel fares better than the classical kernels, its performance increases slowly with the number of training data due to it having a much larger feature space dimension than any of the (approximated) state kernels. Again we see for small amounts of data, the state kernel approximations with low $\chi_{\mathrm{max}}$ have higher validation scores than the full state kernel, particularly with $\chi_{\mathrm{max}} = 2$ for small $N_{\mathrm{t}}$ and $\chi_{\mathrm{max}} = 4$ for intermediate $N_{\mathrm{t}}$. This reinforces the idea that having a simpler model (with a smaller feature space) may provide accuracy benefits in the small-data regime as well as being easier to calculate. This raises the possibility that one could adaptively switch between kernels of varying complexity as more data is obtained, using comparative tools like the Bayes factor to decide which is most suitable at the current optimization stage.

The kernel-based surrogate model used in our Bayesian VQE simulations is an GP. An GP is a collection of random variables with the property that any finite subset has a joint multivariate normal distribution [63]. It is defined in terms of a covariance (a kernel) function $k(\cdot,\cdot)$ and a mean function $\mu(\cdot)$. Given a vector of m observed values $\boldsymbol{y} = (y_1,\ldots,y_m)^T$ of some unknown function at points $\mathbf{X} = (\boldsymbol{x}_1,\ldots, \boldsymbol{x}_m)^T$ (expressed as an m × p matrix, where p in the input space dimension) then the joint distribution with a new point $y^*$ at location $\boldsymbol{x}^*$ is [63]

$\begin{equation} \begin{bmatrix} \ \boldsymbol{y}\ \\ \ y^* \end{bmatrix}\sim\mathcal{N}\left( \begin{bmatrix} \ \boldsymbol{\mu}\ \\ \ \mu(\boldsymbol{x}^*) \end{bmatrix}, \begin{bmatrix} \mathbf{K}&\boldsymbol{k}\\ \boldsymbol{k}^T&k(\boldsymbol{x}^*,\boldsymbol{x}^*) \end{bmatrix} \right), \end{equation} \tag{ 20 }$

where µ is an m element mean vector with $\mu_i = \mu(\boldsymbol{x}_i)$, K is an m × m Gram matrix of pair-wise kernel evaluations between the observed inputs, $K_{ij} = k(\boldsymbol{x}_i,\boldsymbol{x}_j)$, and k is an m element vector of kernel evaluations between the new point $\boldsymbol{x}^*$ and the observed inputs X, $k_i = k(\boldsymbol{x}^*,\boldsymbol{x}_i)$. If noise is present in the observations this is usually approximated to be normally distributed with zero mean. One can show analytically that for observations with normally distributed noise the kernel should be altered to $k(x_i,x_j)\to k(x_i,x_j)+\sigma_n^2\delta_{ij}$, where $\sigma_n^2$ is a noise variance hyperparameter and the indices of δ_ij correspond to indices in the training data meaning that $\mathbf{K}\to \mathbf{K}+\sigma_n^2 I$ but k and $k(\boldsymbol{x}^*,\boldsymbol{x}^*)$ are unchanged. Marginalizing over the observed data, the resulting posterior distribution is itself a normal distribution with a posterior mean (a prediction) $\mathbb{E}[y^*|\boldsymbol{x}^*,\mathbf{X},\boldsymbol{y}] = \hat{y}(\boldsymbol{x}^*)$ given by

$\begin{equation} \hat{y}(\boldsymbol{x}^*) = \mu(\boldsymbol{x}^*) + \boldsymbol{k}^T(\mathbf{K}+\sigma_n^2I)^{-1}(\boldsymbol{y}-\boldsymbol{\mu}) \end{equation} \tag{ 21 }$

and a posterior variance $\mathrm{Var}[\,y^*(\boldsymbol{x}^*)] = {\Delta y(\boldsymbol{x}^*)}^2$ of

$\begin{equation} {\Delta y(\boldsymbol{x}^*)}^2 = k(\boldsymbol{x}^*,\boldsymbol{x}^*) - \boldsymbol{k}^T(\mathbf{K}+\sigma_n^2I)^{-1}\boldsymbol{k}. \end{equation} \tag{ 22 }$

The performance of GP regression depends on the suitability of the kernel function in describing the unknown function. Unless one has significant prior knowledge of the problem at hand, choosing a suitable kernel is often difficult. As a result, problem-agnostic approaches to GP regression often make use of flexible kernels with internal hyperparameters which are varied to fit observed data. A common approach is to find hyperparameter values that maximize the marginal likelihood of the observed data (often the log-likelihood) [63]. Thanks to the normality of the GP's distributions, the marginal log-likelihood $M(\mathbf{X},\boldsymbol{y}) = \log{p(\boldsymbol{y}|\mathbf{X},\alpha)}$ can be calculated analytically for a kernel with hyperparameter α as

$\begin{equation} M(\mathbf{X},\boldsymbol{y}) = -\frac{1}{2}(\boldsymbol{y}-\boldsymbol{\mu})^T(\mathbf{K}+\sigma_n^2\mathbf{I})^{-1}(\boldsymbol{y}-\boldsymbol{\mu})-\frac{1}{2}\log \det[\mathbf{K}+\sigma_n^2\mathbf{I}]-\frac{n}{2}\log 2\pi. \end{equation} \tag{ 23 }$

The gradient of this quantity with respect to a kernel hyperparameter α is given by

$\begin{equation} \begin{split} \frac{\partial M(\mathbf{X},\boldsymbol{y})}{\partial \alpha} = \frac{1}{2}\textrm{Tr}\{(\boldsymbol{\beta}\boldsymbol{\beta}^{T}-\mathbf{K}^{-1}_{\sigma_n})\frac{\partial\mathbf{K}_{\sigma_n}}{\partial\alpha}\},\\ \textrm{where}\ \boldsymbol{\beta} = \mathbf{K}^{-1}_{\sigma_n}\boldsymbol{y},\ \mathbf{K}_{\sigma_n} = \mathbf{K}+\sigma_n^2\mathbf{I}, \end{split} \end{equation} \tag{ 24 }$

and so this quantity can be maximized through gradient ascent at a cost primarily dictated by the calculation of $K_{\sigma_n}^{-1}$ if the kernel evaluations and their gradients are inexpensive to calculate.

The posterior variance (22) depends only on kernel evaluations and not the observations y . If the kernel is fixed, a re-scaling of the training data of the problem $\boldsymbol{y}\to\alpha\boldsymbol{y}$ will leave the posterior variance unchanged whereas it should be a factor of $\alpha^2$ larger. To allow the posterior variance to scale with the size of the observations an external 'signal variance' hyperparameter $\sigma^2$ is often added to the kernel $k(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})\to\sigma^2 k(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})$. We see from (21) and (22) that this does not change posterior mean prediction (up to a re-scaling of $\sigma_n^2$) but the posterior variance is scaled by $\sigma^2$. If the base kernel has no hyperparameters (which is true for the quantum kernels we consider) and the observed data is assumed to be noiseless then it is simple to show the marginal likelihood of the training data is maximized when $\sigma^2_{\mathrm{noiseless}} = (\boldsymbol{y}-\boldsymbol{\mu})\mathbf{K}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})^T/m$ (where K is the Gram matrix for the kernel without $\sigma^2$ and m is the number of observations) – this often serves as a good starting point for maximum likelihood optimization even if noise is present.

The mean function µ is often set to a constant value, most commonly $\mu(\boldsymbol{x}) = 0$ for all $\boldsymbol{x}$, as it is primarily the covariance function which defines a GP's properties; in most cases µ only has a significant impact at points that are far away (with respect to the kernel) from any observed data. However (as shown in section 6) for many ansatz circuits the expected value (taken over the gate angles) of the noiseless energy $\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]$ is given by $\textrm{Tr}\{H\}$ meaning that one should set $\mu(\boldsymbol{x}) = \textrm{Tr}\{H\}$. Alternatively, one can remove this diagonal offset from the Hamiltonian by using $H^{^{\prime}} = H-\textrm{Tr}\{H\}I$ and set µ = 0. We suggest that this should be done generally as calculating or estimating the average energy for a complicated ansatz will be difficult, especially if device noise is present. Removing this offset is good practice in general as a large diagonal offset in H can lead to improper conclusions on the effectiveness of an optimization scheme.

In our simulations we build the GP models using two quantum kernels (described in section 3) and a set of classical kernel functions. These classical kernels are as follows; Matern kernels with smoothness hyperparameter $\nu = 5/2$ and $\nu = 3/2$, the radial basis function (RBF) kernel, and the rational quadratic (RQ) kernel. Their kernel functions are:

$\begin{equation} k_{Matern}(\boldsymbol{x},\boldsymbol{x}^{^{\prime}};\nu = 5/2,l ) = \left(1+\frac{\sqrt{5}}{l}d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})+\frac{5}{3l}d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})^2\right)\exp\left\{-\frac{\sqrt{5}}{l}d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})\right\}, \end{equation} \tag{ 25 }$

$\begin{equation} k_{Matern}(\boldsymbol{x},\boldsymbol{x}^{^{\prime}};\nu = 3/2,l ) = \left(1+\frac{\sqrt{3}}{l}d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})\right)\exp\left\{-\frac{\sqrt{3}}{l}d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})\right\}, \end{equation} \tag{ 26 }$

$\begin{equation} k_{RBF}(\boldsymbol{x},\boldsymbol{x}^{^{\prime}};l ) = \exp\left\{-\frac{d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})^2}{2l^2}\right\},\ \textrm{and} \end{equation} \tag{ 27 }$

$\begin{equation} k_{RQ}(\boldsymbol{x},\boldsymbol{x}^{^{\prime}};l,\alpha) = \left(1+\frac{d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})^2}{2\alpha l^2}\right)^{-\alpha}, \end{equation} \tag{ 28 }$

where $d(\boldsymbol{x},\boldsymbol{x}^{^{\prime}})$ is the squared Euclidean distance between points $\boldsymbol{x}$ and $\boldsymbol{x}^{^{\prime}}$, $l$ are length scale hyperparameters, and $\alpha$ is a scale mixture hyperparameter for the RQ kernel.

These classical kernels require optimization of internal hyperparameters (a length scale $l$ and, for the RQ kernel, a mixing parameter $\alpha$) to ensure that the kernel's infinite-dimensional feature space most-succinctly describes the objective function. By contrast, the state and unitary kernel have finite-dimensional feature spaces in which the (noiseless) energy function is manifestly linear which means no hyperparameter optimization is required for a GP model to make sensible predictions about the energy. Hyperparameter optimization is often the most computationally intensive part of BO as it requires repeated calculation of the Gram matrix K (see (23) and (24)) and its inverse. If we have a kernel with no hyperparameters other than those for the signal variance and noise $k(\boldsymbol{x}_i,\boldsymbol{x}_j) = \sigma^2k_0(\boldsymbol{x}_i,\boldsymbol{x}_j)+\sigma_n^2\delta_{ij}$, where $k_0$ is fixed, then the Gram matrix for the fixed part of the kernel $\mathbf{K}_0$ only needs to be calculated once. This simplifies the hyperparameter optimization as the total Gram matrix $\mathbf{K} = \sigma^2 \mathbf{K}_0 + \sigma_n^2 I$ can be calculated easily for $\sigma^2$ and $\sigma_n^2$ from $\mathbf{K}_0$. Additionally, the eigenvalues/vectors of the full $\mathbf{K}$ can be found analytically in terms of those for $\mathbf{K}_0$, simplifying calculation of the inverse Gram matrices needed to both optimize the hyperparameters and make predictions from the model. This ultimately means that the state and unitary kernels only need to be evaluated once between any two data points within a training data (up to a re-scaling by $\sigma^2$ and a noise offset), in contrast to the classical kernels whose Gram matrices must be entirely recalculated whenever an internal hyperparameter is adjusted.

In BO one attempts to minimize an unknown expensive-to-evaluate objective function by assuming it is randomly drawn from some family of functions [67]. A prior distribution is chosen to encode any initial beliefs about this function, usually that it is sampled from an GP. From observations of the function's values and the prior, a posterior distribution over the chosen family of functions is constructed (the posterior distribution of the GP), quantifying how likely it is that a given function could have produced the observations. An acquisition function is then calculated from the posterior distribution and is maximized/minimized to determine the next point to query. Once the new point has been queried, the posterior distribution is reconstructed with this new data and the process is repeated until convergence to the objective function global minima.

The acquisition function quantifies how 'promising' the querying of a given an unseen point appears for minimizing the objective function. These functions are often designed to balance the exploitation of regions of the parameter space that have already shown good objective function values with the exploration of areas in the parameter space where observations are sparse. By only querying parameter points that appear likely to yield a significant improvement in the objective function, BO can allow global minima to be found with remarkably few iterations. In our simulations we will use the Expected Improvement acquisition function [67] given (for minimization) by $\mathrm{EI}(\boldsymbol{x}) = \mathbb{E}[\mathrm{max}(y_{\mathrm{best}}-\hat{y}(\boldsymbol{x})+\xi,0)|\boldsymbol{x}^*,\mathbf{X},\boldsymbol{y}]$, where the expectation is taken over the surrogate model's posterior distribution, $y_{\mathrm{best}} = \min\{\boldsymbol{y}\}$ is the lowest observation seen so far, and $\xi$ is an exploration hyperparameter. Using of the analytically tractable GP posterior distribution (i.e. normally distributed with mean/variance given by (21) and (22)) we can derive an explicit formula for this [67]:

$\begin{equation} \begin{split} \mathrm{EI}(\boldsymbol{x})&: = (y_\mathrm{best}-\hat{y}(\boldsymbol{x})+\xi)\Phi(Z) + \Delta y(\boldsymbol{x})\phi(Z),\;\ \textrm{where}\ Z = \frac{y_\mathrm{best}-\hat{y}(\boldsymbol{x})+\xi}{\Delta y(\boldsymbol{x})}, \end{split} \end{equation} \tag{ 29 }$

Φ and φ are the cumulative distribution and probability density functions for the standard normal distribution respectively. Expected improvement-based optimization can sometimes be overly greedy and exploitative and so a hyperparameter $\xi$ has been added to control the ratio of exploration to exploitation. The $\xi$ parameter allows more positive observations (up to where $y^*(\boldsymbol{x}) = y_{\mathrm{best}}+\xi$) to count as an improvement, meaning the objective function encourages exploration further away from the point where $y_{\mathrm{best}}$ was observed. In our simulations we used $\xi$ = 0.01 and optimized the surrogate model using the L-BFGS-B gradient-descent algorithm [87].

While it is possible to run BO starting with an observation at a single randomly chosen parameter point, more typically a collection of initial points are used to ensure an initial explorative phase. While there are many designs for these points [88] we use the simplest; random uniform sampling across the parameter space.

Our approximation scheme for the state kernel relies on the parameterized circuit being of the form $\vert \psi(\boldsymbol{\theta}_B) \rangle = U_A(\boldsymbol{\theta}_A)U(\boldsymbol{\theta}_B)U(\boldsymbol{\theta}_C)\vert \boldsymbol{0} \rangle$ where $\boldsymbol{\theta}_A$ and $\boldsymbol{\theta}_C$ are fixed, so that the state kernel can be written in the form (15). We also require that $U_B$ is simple enough that its action on an MPS state of relatively low bond dimension can be simulated classically. If this is the case then by approximating $\vert \psi_A \rangle$ to an MPS $\vert \tilde{\psi}_A^{(\chi_\mathrm{max})} \rangle$, we can achieve a classically tractable approximation to the state kernel.

Figure 9 illustrates our approximation scheme. We begin (9(a)) with a quantum circuit composed of a constant part ($U_A$) and a variational part ($U_B$) whose parameters are to be optimized. For simplicity we omit any $U_C$ that lies in the future light-cone of the variational part of the circuit as this has no effect on the state kernel. The circuit is then converted into a tensor-network representation, which can be done at a cost proportional to the number of gates (an example is shown in 9(b)). The initial state $\vert \boldsymbol{0} \rangle$ is then converted to an MPS with bond dimension χ = 1 and each gate in the circuit is contracted into the MPS to produce $\vert \tilde{\psi}^{(\chi_\mathrm{max})}_A \rangle$ (shown in 9(c)).

Zoom In Zoom Out Reset image size

Figure 10 illustrates the process of contracting a gate G into an MPS $\vert \psi_\mathrm{MPS} \rangle$ (with physical dimension d). As we only consider single- and two-qubit gates (outlined in green and yellow respectively) only two MPS tensors with shared bond dimension $\chi$ are shown (outlined in blue). Vertical open bonds indicate connections to the rest of the MPS state. In 10(a) a single qubit unitary is contracted directly into the MPS tensor at the corresponding site/qubit, altering this tensor but leaving the bond dimensions unchanged. Figure 10(b) shows the contraction of a two-qubit unitary which acts on adjacent sites. The two-qubit unitary is fully contracted into the MPS tensors yielding a joint tensor for the two sites. Singular value decomposition (SVD) is used to split this joint tensor into two new MPS tensors with a new bond dimension $\chi^{^{\prime}}$. The size of $\chi^{^{\prime}}$ depends on the specifics of the unitary and the state to which it is applied but takes a maximum value of $d^2\chi$. Figure 10(c) shows this diagrammatically by giving an alternative scheme for contracting the two-qubit gate. Here the unitary is first decomposed using SVD into two tensors connected by a bond of dimension $d^2$. Each resulting tensor is contracted into their connected MPS tensor yielding an MPS state with a two bonds between the new MPS tensors (of size χ and $d^2$). These bonds are combined to form a single bond of dimension $d^2\chi$. Two qubit gates on non-neighboring sites can be contracted using the same scheme by applying a chain of SWAP gates (as described in [89]) but can significantly increase the bond dimensions.

Zoom In Zoom Out Reset image size

Applying many two-qubit gates will lead to exponential growth in the MPS's bond dimensions. To avoid this we truncate $\chi^{^{\prime}}$ by only retaining the $\chi_\mathrm{max}$ largest singular values obtained from the SVD in the final stage of 10(b). Following truncation the MPS is re-normalized, yielding an approximation $\vert \tilde{\psi}_{MPS} \rangle$ to the transformed state $G\vert \psi_\mathrm{MPS} \rangle$ with a maximum bond dimension of $\chi_{\mathrm{max}}$. If significant singular values of the joint tensor are discarded, the fidelity between $\vert \tilde{\psi}_{MPS} \rangle$ and the target state $G\vert \psi_\mathrm{MPS} \rangle$ will be degraded. We partially mitigate this by maximizing the fidelity $\vert \langle \tilde{\psi}_{\mathrm{MPS}} \vert G\vert \psi_\mathrm{MPS} \rangle \vert^2$ between these two states, optimizing over $\vert \tilde{\psi}_{\mathrm{MPS}} \rangle$. This procedure can be greatly simplified by only optimizing over the two MPS that are altered when producing $\vert \tilde{\psi}_{MPS} \rangle$ (shown as yellow dots with a blue outline in figures 10(b) and (d)). By first contracting over the tensors which are not optimized (creating the tensor shown with red fill figure 10(d) this optimization can be performed extremely efficiently, involving contractions over just 4 tensors. These initial contractions can also be further simplified by first canonicalizing $\vert \psi_\mathrm{MPS} \rangle$ around the two qubits acted upon by the unitary [82] (provided the MPS has open boundaries to allow this). Once all gates in the circuit have been contracted into the MPS using the scheme we have outlined (with the necessary truncation and fidelity optimization) the resulting MPS state $\vert \tilde{\psi}_A^{(\chi_\mathrm{max})} \rangle$ is used as the input state for the approximated state kernel given in (15).