Physics-Informed Bayesian Optimization of Variational Quantum Circuits
Authors: Kim Andrea Nicoli
Christopher J. Anders
Lena Funcke
Tobias Hartung
Karl Jansen
Stefan Kuhn
Klaus Robert Muller
Paolo Stornati
Pan Kessel
Shinichi Nakajima
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian optimization
Expected improvement
Quantum computing
Variational Quantum Eigensolvers
TL;DR: We propose a novel physics-informed kernel and a new acquisition function for the optimization of variational quantum eigensolvers (VQEs).

Abstract: In this paper, we propose a novel and powerful method to harness Bayesian optimization for variational quantum eigensolvers (VQEs) - a hybrid quantum-classical protocol used to approximate the ground state of a quantum Hamiltonian. Specifically, we derive a *VQE-kernel* which incorporates important prior information about quantum circuits: the kernel feature map of the VQE-kernel exactly matches the known functional form of the VQE's objective function and thereby significantly reduces the posterior uncertainty.
Moreover, we propose a novel acquisition function for Bayesian optimization called \emph{Expected Maximum Improvement over Confident Regions} (EMICoRe) which can actively exploit the inductive bias of the VQE-kernel by treating regions with low predictive uncertainty as indirectly "observed". As a result, observations at as few as three points in the search domain are sufficient to determine the complete objective function along an entire one-dimensional subspace of the optimization landscape. 
Our numerical experiments demonstrate that our approach improves over state-of-the-art baselines.

Physics-Informed Bayesian Optimization of Variational Quantum Circuits
Kim A. Nicoli
Christopher J. Anders Lena Funcke
Tobias Hartung Paolo Stornati 10
Karl Jansen
Stefan Kühn" Klaus-Robert Müller
Pan Kessel" Shinichi Nakajima 3,4412
3,4,12
Transddicippliary Research Area (TRA) Matter, University of Bonn, Germany HHelmholtz Institute for Radiation and Nuclear Physics (HISKP)
Berlin Institute for the Foundations of Learning and Data (BIFOL.)
Technische Universität Berlin, Germany, Northeastern University London, UK
6 Khoury College of Computer Sciences, Northeastern University, USA
CQTA, Deutsches Blektronen--Synchorroe (DESY), Zeuthen, Germany
8 Department of Artificial Intelligence, Korea University, Korea MMax Planck Institut für Informatik, Saarbrücken, Germany Institute Photonic Sciences, The Barcelona Institute of Science and Technology (ICFO)
11 Prescient Design, gRED, Roche, Switzerland, 12RIKEN Center for AIP, Japan
Abstract
In this paper. we propose novel and powerful method to harness Bayesian op- timization for Variational Quantum Eigensolvers (VQEs)-a hybrid quantum- classical protocol used to approximate the ground state of quantum Hamiltonian. Specifically, we derive a VQE-kernel which incorporates important prior informa- tion about quantum circuits: the kernel feature map of the VQE-kernel exactly matches the known functional form of the VQE's objective function and thereby significantly reduces the posterior uncertainty. Moreover, propose novel acqui- sition function for Bayesian optimization called Expected Maximum Improvement over Confident Regions (EMICoRe) which can actively exploit the inductive bias of the VQE-kernel by treating regions with low predictive uncertainty as indirectly "observed". As result, observations as few as three points in the search domain are sufficient to determine the complete objective function along an entire one- dimensional subspace of the optimization landscape. Our numerical experiments demonstrate that our approach improves over state-of-the-ee baselines.
Introduction
Quantum computing raises the exciting future prospect to efficiently tackle currently intractable problems in quantum chemistry [], many-body physics [2], combinatorial optimization [3], machine learning 4 and beyond. Rapid progress the last years has led the development of the first noisy intermediate-
quantum (NISQ) devices [3] quantum hardware with several hundreds of qubits that can be harnessed to outperform classical computers
specific tasks (see, e.g., Ref. (6). However, these tasks have been specifically designed to be hard on classical computers and easy on quantum computers, while being of practical use. As have entered the NISQ era, of the grand challenges of quantum computing lies in the development of algorithms for NISQ devices that may exhibit quantum advantage task of practical relevance.
'These authors contributed equally to the work.
Correspondence mnicci@uui-boo de and {anders, makaajiaalttt-bbrllll
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

One promising approach toward using NISQ devices is hybrid quantum classical algorithms, such as VQEs 7 []. which can be used to compute ground states of quantum Hamiltonians. VQEs can be seen as the quantum counterpart of neural networks: while classical neural networks model functions by parametric layers, the parametric quantum circuits used in VQEs represent variational wave functions by parametric quantum gates acting on qubits. During training, we aim to find suitable choice for the variational parameters of the wave function such that the quantum mechanical expectation value of the Hamiltonian is minimized. From the optimization perspective, the quantum mechanical nature of the energy measurement is of little relevance. The training of VQEs can thus be regarded specific, albeit particularly challenging, noisy black-box optimization problem. Namely, we solve
ni022(002x) (*));,
where (a(): is estimated from costly noisy observation, y f (az) + e, on the quantum device.
Efficiently exploring the energy landscape of the VQE is crucial for successful optimization. This requires leveraging strong prior knowledge about the VQE objective function f"().. For instance, gradient-based optimization methods commonly utilize the parameter shift rule, which takes advan- tage of the specific structure of variational quantum circuits to compute gradients efficiently 99 [0]. Another approach is the Nakanishi- -ujtisshi-FFii--- Fuji-Todo (NFT) method [TT] a version of sequential minimal optimization (SMO) [2]. which solves sub-problems sequentially. Nakanishi et al. [1]] derived the explicit functional form of the VQE objective, enabling coordinate-wist global optimization with only two observations per iteration.
Given the highly non-trivial nature of classical optimization in VQE, machine learning represents an appealing tool to leverage the informative underlying physies toward a more effective search for global optimum. Bayesian Optimization (BO) 13 4 has been recently applied to VQE [15]. These methods have the distinct advantage that they can take the inherently noisy nature of the NISQ circuits into account. Unfortunately, these methods are yet to be competitive, especially in the high dimensional regime, because of their poor scalability and
hetptapphssizddeepllrttttfe behavior [16].
To overcome this limitation, this paper introduces novel and powerful method for BO of VQEs, capitalizing on VQE-specific properties as physics-informed prior knowledge. To this end, we propose a VQE-kernel, designed with feature vectors that precisely align with the basis functions of the VQE objective. This incorporation of a strong inductive bias maximizes the statistical efficiency of Gaussian process (GP) regression and guarantees that GP posterior samples reside within the VQE function space. To further harness this powerful inductive bias, we present novel acquisition function named Expected Maximum Improvement over Confident Regions (EMICoRe). EMICoRe operates by initially predicting the posterior variance and treating the points with low posterior variances as "observed" points. Subsequently, these indirectly observed points, alongside the directly observed ones, form the Confident Regions (CoRe) on which safe optimization of the GP mean is conducted to determine the current optimum. EMICoRe evaluates the expected maximum improvement of the best points in CoRe before and after the candidate points are observed. By utilizing EMICoRe, our approach combines the strengths of NFT [11] and BO, complementing each other in synergistic manner. BO enhances the efficiency of NFT by replacing sub-optimal deterministic choices of next observation points, while the NET procedure significantly constrains the exploration space of BO, leading to remarkable improvements in scalability.
Our numerical experiments demonstrate the performance gains of using the VQE kernel, and signifi- cant improvement of our NFT-with- Wwithh-EMICoRe.... approuch in VQE problems with different Hamiltonians, different numbers of qubits, etc. As an additional contribution, we prove that two known important properties of VQE, the parameter shift rule O and the sinusoidal function-fome [1I, are equivalent, implying that they are not two different properties but two expressions of a single property.
Related Work Numerous optimization methods for the VQE protocol have been proposed: gradient- based methods often rely on the parameter shift rule 9 10), while NFT harnesses the specific functional form of the VQE objective [II] to establish SMO. BO has also been applied for VQE minimization [15]. Therein, it was shown that combining periodic kernel [17 and noisy expected improvement acquisition function [18] improves the performance of BO over the standard RBF kernel, and can make BO comparable to the state-of-the-art methods in the regime of small qubits and high observation noise.
2

---

d
22
a)
b) EMICoRe step
Update GP, fit line
d)
2
Og
EMICORe
Cp
máái
#
min
dd
20
2x
mmn
a)
b) NFT step (a
(c) Fit line, find min
d'
NFT
misn
Figure 1: Illustration of EMICoRe (ours) and NFT [11] (baseline) procedures. In each step of NFT (bottom row), a) given the current best point &' and the direction d' to be explored, b) next observation points are chosen in a deterministic fashion. Then, c) the minimum along the line is found by the sinusoidal function fitting to the three points, and d) the next optimization step for the new direction d'+1 starts from the found minimum x x+1 = _ min- The EMICoRe procedure (top row) uses GP regression and BO in the steps highlighted by the light blue box: b) the observation points are chosen by BO with the EMICoRe acquisition funccionn based on the GP trained on the previous observations, and c) minimizing the predictive mean function of the updated GP with the new observations gives the best point 2 aimi.
including hyperparameter tuning of deep neural networks [23]. Most work on BO uses GP regression, and computes an acquisition function whose maximizer is suggested as the next observation point. Many acquisition functions, including lower confidence bound (24 25 probability of improvement (26). Expected Improvement (EI) [27], entropy search 28 and knowledge gradient (30) have been proposed. The most common acquisition function is EI, of which many generalizations exist: noisy EI (NEI) [8] [3] for dealing with observation noise, parallel El [31 for batch sample suggestions, and El per cost [4 for cost sensitivity. Our EMICoRe acquisition function is generalization of noisy EI and parallel EI with the key novelty introducing CoRe, which defines the indirectly observed points. Note the difference between the trust region [16] and CoRe: based on the predictive uncertainty, the former restricts the region to be explored, while the latter expands the observed points.
BO is a versatile tool for black- box optimization in many applications 19 20
2 Preliminaries
2.1 Bayesian Optimization
Let ((() X IR be an unknown (black-box) objective function to be minimized. BO 13 [E]
approximates the objective with a surrogate function f(x:) f ((xx), and suggests promising points to be observed in each step, such that the objective is likely to be improved considerably. A commonly used surrogate is the Gaussian process (GP) regression model [32] with dimennnnala Gaussian likelihood and GP prior,
plylar, Myla,,()) f()) No (y; f(x),
which is trained on the (possibly) noisy observations y= fx) e made in the previous iterations. Here, o² is the variance of observation noise e, and v(-) and k(---) are the prior mean and the kernel (covariance) functions, respectively. Throughout the paper, we set the prior mean function be zero, i.e., v(x) = 0, Vx E X Let X, y) } be N training samples, where X (æI...., ΣN) E V and = ( _ E R Since the GP prior is conjugate for the Gaussian likelihood @ the posterior and thus, for arbitrary M test inputs X' = (æ' 1e'1..... x' Mr) E AM the posterior of the function values E M is the M dimeesniioaal Gaussian with its mean and covariance
of the GP regression model (2) is also a GP, i.e., p([(-)X y) = GP((O); px (xx)(,,) (x),x((:)))))))))))))))))
f' =
"See Chap. of Ref. 33 for the conjugacy in Bayesian inference.
3

---

analytically given
x NN S'x, where
(3)
n'x = KK' (K + + y: = K" (K (K+0 K'.
Here, K = k X) E RNNN K' = kXX E RNNM, and K" = 1(X', E RMXM are the train, train-test, and test kernel matrices, respectively, where k X') denotes the kernel matrix evaluated at each column of X and X' such that (k(X,X'))aa
= kxt,, Lmm). Moreover, IN E TRNN denotes the identity matrix, and the subscript X of posterior means and covariances specifies the input points on which the GP was trained (see Appendix B for details of GP regression)
In the general BO framework 14, M > 1 input points to be observed next are chosen by (approxi- mately) solving the following maximization problem in each iteration: max ax-=1 (X'), where X' denotes the training input points observed until the /-th iteration, and ax (+) is acquisition function computed based on the GP trained on the observations y at X. The acquisition function evaluates the promising-ner ness of the new input points X', and should therefore give high value if ob- serving X' likely improves the current best score considerably. Common choices for the acquisition function are EI 27 34,
x (x')= = f p(/' p(/'( and its variants. Here, f denotes the current best observation, i.c., f (1...N) f(xr). ). and (-)p denotes the expectation value with respect to the distribution p. El covers the case where the observation noise is negligible and only one sample is chosen in each iteration, i.e., o² < 1, M 1, and its analytic solution makes EI handy (see Appendix In the general case where 0² > 0, M 2 1, generalizations of EI should be
considered, such as NE 18 B1,,
NEI (X') (max min(f) min (f (x)))
(5)
which appropriately takes into account the correlation between observations. NEI is estimated by quasi Monte Carlo sampling, and its maximization is approximately performed by sequentially adding a point to X' until M points are collected.
2.2 Variational Quantum Eigensolver (VQE)
The VQE [7] is hybrid quantum-classical computing protocol for estimating the ground-state
energy of a given quantum Hamiltonian for a Q-qubit system. The quantum computer is used to prepare parametric quantum state (2))), which depends on D angular parameters æ E X 10, 2nr 2 This trial state ( is generated by applying D' D) quantum gate operations, G(x) = GD o
G1, to an initial quantum state (10),, i.e., ( = G())) All gates (G |2'- are unitary, and we assume in this paper that xa parametrizes only a single gate Ga (ra), where d' (d) specifies the gate parametrized by Tg. Thus, D of the D' gates are parametrized by each entry of æ exclusively!" IV We consider parametric gates of the form Go (I) (2) = exp (-ixPa /2), where Pa is an arbitrary sequence of the Pauli operators {ap og oy Jq=1 acting on each qubit at most once. This 2
form covers not only the single-qubit gates such as Rx = cxp (-i00rx) ) but also the entangling gates, such as (x) = exp (-ino a. o xxx and Rzz(r) = exp (-ino ixo 2 q o0%) for q * 92.
which are commonly realized trapped-ion quantum hardwares 35 36
The quantum computer evaluates the energy of the resulting quantum state (Wa) by observing
== fx) +E,
where f())) (v (() (vo]G(() (G)))))))(G))))))
(6)
and i denotes the Hermitian conjugate. The observation noise in our numerical experiments will only incorporate shot noise, and we will not consider the hardware-dependee imperfect qubits, gates, and measurements. For each observation, multiple readout shots Nahots are acquired to suppress the variance o (Nooosts of the noise e. Since the observation y is the sum of many random variables, it approximately follows Gaussian distribution, according to the central limit thcorem. The Gaussian likelihood in the GP regression model (2) therefore approximates the observation y well if f(x) R (*) and x² 2 2 With the quantum computer that
errors induced by
Note that the subscript does not necessarily specify dependencies. For example, 2'xx also depends on y. IV In Appendix we discuss how the theory and our method can be extended to the non exclusive parametriza- tion case.

---

provides noisy estimates of (x), a classical computer solves the minimization problem 1 and finds the minimizer ă.
Several approaches, including stochastic gradient descent (SGD) 99 10 SMO 11 and BO [5]
have been proposed for solving the VQE optimization problem ① Among others, sate-of-theeart methods effectively incorporate unique properties of VQE. Let {ea) be the standard basis.
(6) for any parametric quantum circuit G(.), Hermitian operator H, and initial state (wo) satisfies
Proposition 1.
[9] (Parameter shift rule) The VQE objective function (-) in Eq.
0 2 Bag (x) f (x +x+ Rea) - f2 (x- Vo E [0, (m) 1...., D.
(7)
Most SGD approaches rely on Proposition [ which allows accurate gradient estimation from 2 D observations. Another useful property is used for tailoring SMO [2] to VQE:
Proposition 2. II For the VQE objective functiinn (-) in Eq. 6 with any G(-), H, and (oo).
3b E 3 such that
(a) bb vec (c0= (1, cos Fd, sin sin T) Va E [0, 2m)
(8)
where & and vec(.) denote the tensor product and the vectorization operator tensor, respectively,
Proposition provides a strong prior knowledge on the VQE objective function (*(+) we fix the function values at three general points along a coordinate axis, e.g., {x, x + 25 ea for any and d, then the whole function along the axis, i.e., (" + oea), Vo E [0, 2T) is fixed because is a first-order sinusoidal function. Leveraging this property, version of SMO was proposed [T] This optimization strategy, named after its authors Nakanishi-Fujje Fuji-Todo (NFT), finds the global optimum in single direction by sinusoidal fitting from three function values (two observed and one estimated) at each iteration, and showed state-of-the-art performance (see Appendix F. for details of NFT).
3 Proposed Method
In this section, introduce our approach that combines BO and NFT by using novel kernel and a novel acquisition function. After introducing these two ingredients, we propose our approach for VQE optimization. Lastly, prove the equivalence between Propositions I and [ We refer the reader to Appendix F in particular Algorithms to for the pseudo-codes and further algorithmic details complementing the brief introduction to the algorithms presented in the following sections.
3.1
VQE Kernel
We first propose the following VQE kernel, and use for the GP regression model
722 (?²+2con(aa-d)) )
A VOU (z. x')
oî 11
7242
where oZ. ,2 O are kernel parameters. of corresponds to the prior variance, while 72 controls the smoothness of the kernel. For 72 = 1, the VQE kernel is the product of Dirichlet kernels [37, each of which is associated with each direction d. We can show that this kernel (9) exactly gives the finite-dimensional feature space specified by Proposition (the proof is given in Appendix
Theorem 1. The VQE kernel (9)
is decomposed as
V (x, x') = (xx) To(("',, where (xx) (-2 2)0/2 vec (0 V2 cos v2)
do
Let VQQ be the set of all possible VQE objective functions specified by Proposition 2 Theorem guarantees that the support of the GP prior in Eq. (2) matches VVQ if we use the VQE kernel function (9) with prior mean function such that v E FVQE. Thus, the VQE kernel drastically limits the expressivity of GP without model misspecification, which enhances statistical efficiency. A more important consequence of Theorem Tis that Crppositionnhhllee for any sample from the GP posterior with the VQË kernel 9 This implies that if GP is certain about three general points along an axis, it must be certain about the whole 1-D subspace going through those points. Theorem can be generalized to the non-exclusive case, where each entry of x may parametrize multiple gates simultancously (see Appendix E
5

---

3.2 EMICoRe: Expected Maximum Improvement over Confident Regions
In GP regression, the predictive covariance does not depend on the observations y (see Eq.
(4),
implying that, for given training inputs X E X we can compute the predictive variance sx (æ, a) any æ E X before observing the function values. Let us define Confident Regions (CoRe)
Ex = 8x .) S x²)
(10)
which corresponds to the set on which the predicted uncertainty by GP is lower than a threshold k. For sufficiently small K, appropriate kernel (which does not cause model misspecificaation) and sufficiently weak prior, the GP predictions on CoRe are already accurate, and therefore CoRe can
be regarded as "observed points." " This leads the following acquisition function, which we call Expected Maximum Improvement over Confident Regions (EMICoRe):
EMICoRe (X')=
M (max (0, min f(x) min f(x))) ((SO)x)))
(11)
where X XX) X-M+ denotes the augmented training set with the new input points X' E XMM This acquisition function evaluates the expected maximum improvement (per new observation point) when CoRe is expanded from Zx to Xx by observing the objective at the new input points X', EMICoRe can be seen as a generalization of existing methods. If CoRe consists of the training points, it reduces to NEI [18] If we set - oc so that the whole space is in the CoRe, and the random function f(-) Eq. (TT) replaced with its predictive mean of the current (first term) and the updated (second term) GP, EMICoRe reduces to knowledge gradient (KG) [30]. Thus, KG can be seen as a version of EMICoRe that ignores the uncertainty of the updated GP.
3.2.1 NFT-with-EMICoRe
We enhance the state-of-the-an NFT approach [1] with the VQE kernel and EMICoRe. We start from a brief overview of the NFT algorithm (detailed algorithms of NFT, NFT-with-EMIIoRR, and the EMICoRe subroutine are given in Appendix
NFT: First, we initialize with a random point a with 1 = S² (*x")) 8. Then, for each iteration step we proceed as follows:
1. Select an axis d E (1...., D} sequentially or randomly and observe the objective y
(v(, T at deterministically chosen two points X' = (zf, = (a* (at-1 2n 2tr33es. a"
2Tr (3ea) along the axis d.
2. Fit the sinusoidal function F(0) = co c cos 0 + sin @ to the two new observations vi, 1/2 as well as the previous best estimated score g²-1 The optimal shift 9 that minimizes f()) is analytically computed, which is used to get the new optimum s² = _''1_ Ōea.
3. The best score updated g f((").
We stress that if the observation noise is negligible, i.c., y RJ f(x), each step of NFT reaches the global optimum in the one-dimensional subspace along the chosen axis d, and thus performs SMO, see Proposition 2 In this case, the choice of the two new observation points X' is not important, as long as any pair of the three points not exactly the same. However, when the observation noise is significant, the estimated global optimum in the chosen subspace is not necessary accurate, and the accuracy highly depends on the choice of the new observation points X', In addition, errors can be accumulated in the best score ii', and therefore an additional measurement needs to be performed at a' after a certain iteration interval.
NFT-with-EMICoke (ours): We propose to apply BO to NFT by using the VQE kernel and EMICoRe. Specifically, we use the GP regression model with the VQE kernel as a surrogate for BO, and choose new observation points X' by using the EMICoRe acquisition function. NFT-EMICORe starts from TNFT NET iterations until GP gets informative with a sufficient number of observations. After this initial phase, proceed for each iteration follows:
1. Select an axis d E (1...., D} sequentially. and new observation points X' by BO with EMICoRe, based on the previous optimum "-1 the previous training data ( (x111 y''1), and the current CoRe threshold x² (this subroutine for EMICoRe will be explained below).

---

2. We observe y at the new points X' chosen by EMICoRe, and train GP with the updated y (y' "y'T).
training data X' = (x' x')
3. The subspace optimization is performed by fitting sinusoidal function f(0) to the GP posterior means 1 (u(0²=1 2 /3ea),
p=(a="), / p(("=1 + 2 at three points. With the analytic minimum of the sinusoidal function, optimum is computed: a' = _²-1 + dea.
argming f(0), the current
4. For the current best score, we simply use the GP posterior mean já = ==(²") and we set the new CoRe threshold
K 4+1 61- TT-TTvv--
(12)
Tave
Note that the GP posterior mean function lies in the VQE function space FVQE and therefore the fitting is done without error, and the choice of the three points in 3 Note also that the CoRe threshold xd+1 adjusts the required accuracy to the average reduction of the best score over the Tave latest iterations-an the early phase where the energy p decreases steeply, large K encourages crude optimization, while in the converging phase where the energy decreases slowly, a small k enforces accurate optimization. Figure Iillustrates the procedures of NFT and our EMICoRe approach.
affect the accuracy.
The EMICoRe subroutine receives the previous optimum a'-1, the previous training data {x . " and the current CoRe threshold k², and returns the new points X' that maximizes EMICoRe. We fix the number of new samples to M = 2, and perform grid search along the chosen direction d. To this end,
1. We prepare Jsa(Jsa 1) combinations of JsG search grid points ( + (xza))=ii Jsa 2x
where (a1,... ...sa
JsG) T, as candidate set C (X'
Jac
IR Dx2) Jsa
j==
2. For each candidate X E C, we compute the updated GP posterior variance 8x (I., 2), Vo Grid where (x is the posterior covariance function of the GP trained on the
augmented training points X
(X'-1, X), and G (a'=1 + ayea) yog with JoG) are JoG grid points along the axis Joa+1
a = (a1,...a.oo)
2x
3. We obtain discrete approximation to the updated CoRe as 2x = fo E x Grid, sx (xx x) S n²). For simplicity, we approximate the previous CoRe the pre- vious optimum, i.c., 2x-- = =fa-1)) After computing the mean and the covariance of the previous GP posterior, p(f, ff x "²=1), at the previous best point a' and the updated CoRe points (as the test set test xx) which is D-dimensional Gaus- EMICoRe = M (max(0, f
sian for Đ = 1 U 2x1 = 1 + 2 -wwe min min(frest)
we estimate a x
(ftest))) ))) pfff x x y by quasi Monte Carlo sampling.
The subroutine iterates this process for all candidates, and returns the best one,
X'
argmax
EMICoRe ax
(X):
Xec
Parameter Setting Our approach has the crucial advantage that the sensitive parameters can be automatically tuned. The kernel smoothness parameter 7 is optimized by maximizing the marginal likelihood of the GP, and the CoRe threshold k is set to the average energy decrease of the last iterations as explained above 7 The noise variance o² is set by observing f f*(x)) several times at several random points and estimating x² X 2N For the GP prior, the zero mean function ((x) 0, Væ is used, and the prior variance of is roughly set so that the absolute value of the ground-state energy is in the same order as ag. Other parameters, including the number of grid points for search and CoRe discretization, should be set to sufficiently large values. See Appendix F for more details of parameter setting.
In Appendix F.4, we investigate other heuristics for the update, and find that the default update rule performs comparably to the best heuristic.
12
7

---

3.3 Equivalence Between VQE Properties
Our method, properties of VQE, namely the highly constrained objective function form (Proposition 2 One may now wonder if we can further improve NFT-EMICORe by harnessing the other property, e-g., by regularizing GP so that its gradients follow the parameter shift rule (Proposition DT. We give a theory that answers to this question. Although the two properties were separately derived by analyzing the VQE objective ( they are actually mathematically equivalent (the proof is given in Appendix G
adapts BO to the VQE problem by harnessing one of the useful
Theorem 2. For any periodic function f : [0, 2T) D + R. Eq.
Eq. 8
This means that any sample from the GP prior with the VQE kernel (and any prior mean such that v()) E FVQE) already satisfies the parameter shift rule. Theorem implies that the parameter shift rule and the VQE function form are not two different properties, but two expressions of a single property, which is an important result that-to our knowledge-was not known in the literature.
4 Experiments
We numerically demonstrate the performance of our approach for several setups, where the goal is to find the ground state of the quantum Heisenberg Hamiltonian
H (Jxoy of + Jyooo J20f + EJ- + thyo + h2of)] (13)
where for oj o²) represent the Pauli matrices applied on the qubit in site j. The Heisenberg
Hamiltonian, which represents a standard benchmarks of high practical relevance, is commonly used for evaluating the VQE performance (see, e.g.. (8 and its ground -rouud ground state |pas)
along with the ground-stad state energy f = (Wess H (as) which gives a lower-bound of the VQE
objective ( can be analytically computed for small Q. For the variational quantum circuit G(x), we use the T-layered Efficient (2) circuit (see Appendix C with the open boundary, for which the search domain of the angular variables is x E [0, 2x) with D = (2 + (L- Q. We classically simulate the quantum computation with the Qiskit (43] library, and our Python implementation along with detailed tutorials on how to reproduce the results is publicly available on GitHub 44 at
hhtps:/ggihub.comm// cow./micore//emmcrrr
To measure the optimization performance, i.e., the quality of the final solution a² after T iter- ations, we use two metrics: the true achieved lowest energy after T steps, ENG = f²")) (Wol G(a²) HG("")) (w)). which is evaluated by simulating the noiseless observation with Nohuots - Do, and the fidelity FID = (vuslvba)) = (Wasl G(#") (o)) which measures how similar the best solution is to the true ground state. As the cost of observations, we count the total number of observed points, ignoring the cost of classical computation. We test each method 50 times, using the same set of initial points, for each method, for fair comparison. The initial points were randomly drawn from the uniform distribution in [0, 2 for fair comparisons. Details of the VQE circuits and the experimental settings can be found in Appendices C and respectively.
4.1 VQE-kernel Analysis
We first investigate the benefit of using the VQE-kernel for the following optimization problem: the task is minimize the VQE objective (6) for the Ising Hamiltonian, special case of the Heisenberg
Hamiltonian with the coupling parameters set to Jx = -1, Jy = Jz = 0, hix = hy = 0, h2 = 1, with the (L layeeed quantum circuit with (Q- -uuuusss.. Namely, we solve the problem ①
in the (D- = 2(L + 1)Q = 24) -iidennnsio space. We use the standard BO procedure with GP regression, and compare our VQE-kernel with the Gaussian radial basis function (RBF) + and the periodic kernel [17. where the El acquisition function is maximized by L-BFGS 45 Figure
shows the achieved energy (left) and the fidelity (right) with different kernels. In each panel, the left
"pin chain Hamiltonians are widely studied in condensed matter physics (39 and generalized spin chains also emerge from the lattice discretization of field theories in low dimensions (see, e.g., Ref. [40] or Refs. [1] 52 for reviews).
8

---

VQE (ours)
Periodic
1.0
Ground State
-2.75
Energy
Fidelity
-3.00
0.5
-3.25
-3.50
Ground State
0.0
100
150 Trial
100
Trial Density
#Observed Paints
Density
Observed Points
Figure 2: Comparison of our VQE-kernel (red) to the RBF and the periodic kernel benchmarks (blue and orange) in the VQE optimization using the standard BO procedure, for the Ising Hamiltonian with the (L = 3) -layered (Q 3)-qubits quantum circuit. The search domain dimension is D 24, and Nahota = 1024 readout shots are taken for each observation. The energy (left) and the fidelity (right) are plotted, and in each plot, optimization progress is shown with the median (solid) and the 25- and 75-th percentiles (shadows) over 50 trials. The portrait square shows the distribution of the final solution after 150 observations have been performed.
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
Ground
Energy
-5.0
Fidelity
0.5
-5.5
1st Exched Stata
6.0
Graund State
0.0
200
400
600 Trial Density
0
200
400
600 Trial Density
# #Observed Points
Observed Points
1.0
Ground
Energy
Fidelity
0.5
1st Excitd State
0.0
200
400
600 Trial
200
400
600 Trial
# Observed Points
Density
# #Observed Points
Density
Figure 3: Comparison (in the same format as Figure 2 between our EMICoRe (red) and the NFT baselines (green and purple) in the VQE for the Ising (top row) and Heisenberg (bottom row) Hamiltonians with the (L = 3a-leye (Q 5)-qubits quantum eircuit (thus, D 40) and Nshoot = 1024. We confirmed for the Ising Hamiltonian that longer optimization by EMICoRe up to 6000 observed points reaches the ground state with 98% fidelity (see Appendix I.1.
plot shows the progress of optimization (median as solid curve and the 25- and 75-th percentiles as shadows) as a function of the observation cost, i.e., the number of observed points after corresponding iterations. The portrait square on the right shows the distribution (by kernel density estimation [46]) of the best solutions after 150 observations have been performed. We see that the proposed VQE kernel (red), which achieves 0.93 0.05 fidelity, converges more stably than the baseline RBF and periodic kernels, both of which achieve 0.90 + 0.09 fidelity. Therefore, the VQE kernel leads to better statistical efficiency albeit its effect appears rather mild the regime a small number qubits.
NrT-with-EMICoRe Analysis
We now evaluate our EMICoRe approach and show that this improves the state-of-the-an baselines of NFT [T]. Specifically, we compare our method with two versions of NFT: NFT- Sequential updates along each axis sequentially, while NFT-Randon chooses the axis randomly. S the results of VQE for the Ising Hamiltonian (top row) and the Heisenberg Hamiltonian
(bottom row) with the parameters set to Jx Jy Jz = 1 and hx = hy hz 1. For both

---

f (x)
GP
(-O 0,
[X] 1, d=
t=1. 1, [X] 3. d-1
t- 293,
[X] 587, d'=13
(= 294,
[X] - 589,
d'=14
2
x
x
ge
X
5
-
-2
0
0
0
0
-5
5
X2
X
π
Figure Evolution of the GP in NFT-with-EMICoRe The posterior mean (blue solid) with uncertainty (blue shadow) along the direction d to be optimized in the steps / = 0, 1, 293, 294 (columns) are shown before (top row) and after (bottom row) the chosen two new points X" (a4, x4) are observed. The red solid curve shows the true energy f* (x)) On the top, the number of observed samples X until step t and the direction are shown.
Hamiltonians, we use (L 3y-leyered (Q = 5)-bbbsts variational quantum circuits, thus D 40, with each observation performed by Nohots = 1024 readout shots. The figure shows that our NFT. With-EMICoRe method consistently outperforms the baselines in terms of both achieved energy (Ieft) and fidelity (right). Moreover, the variance over the trials is reduced, implying its robustness against the initialization. We conducted additional experiments to answer the important questions of whether our approach converges to the ground state with high fidelity, and how the Hamiltonian, the number of qubits, and the number of readout shots affect the optimization performance. The results reported in Appendix [1 show that our WETT-withEMMCCor with 98% % fidelity after 6000 observed points, see Appendix [.] and consistently outperforms the baseline methods in various settings of Hamiltonians, qubits, and the number of readout shots, see Appendices
method reaches the ground state
I.2 and I.3] Thus, we conclude that our approach of combining NET and BO with the VQE-kemel and the EMICoRe acquisition function is suited for efficient VQE optimization, dealing with different levels of observation noise.
Figure 4 shows the evolution of the GP during optimization by NFT-witthFMMCoRe The blue curve and the shadow show the posterior mean (solid) and the uncertainty (shadow) at four different steps (columns). The top and the bottom rows show the GP before and after the two new points X = (x', xτ)), chosen by EMICoRe, are observed. The red solid curve shows the true energy f"(x). We observe the following: In the early phase (the left two columns), GP is uncertain except at the current optimum é before new points are observed; After observing the chosen two points X', GP is certain on the whole subspace, thanks to the VQE kernel; and in the converging phase (the right two columns), GP is certain over the whole subspace even before new observations, and the fine-tuning of the angular variable æ is performed.
Conclusion
Efficient, noise-resilient algorithms for optimizing hybrid quantum-classica algorithms are pivotal for the NISQ era. In this paper, we propose EMICoRe, a new Bayesian optimization approach, in combination with a novel VQE-kernel which leverages strong inductive biases from physics. Our dnysics-informed search with EMICoRe achieves faster and more stable convergence the minimum energy, thus identifying optimal quantum circuit parameters. Remarkably. the physical inductive bias makes the optimization more resilient observation noise. The insight that more physical information helps to gain better practical understanding of quantum systems is starting point for further algorithmic improvements and further exploration of other possible physical biases. This may ultimately lead to more noise-resilient algorithms for quantum circuits, thus enabling the accessibility of larger-scale experiments in the quantum computing realm. In future work, we plan investigate the effect of herdware-dependent noise, and test our approach on real quantum devices.
10

---

Acknowledgments
The authors thank the referees for their insightful comments, which significantly improved the paper. This work was partially supported by the German Ministry for Education and Research (BMBF) as BIFOLD- - Berlin Institute for the Foundations of Learning and Data (BIFOLD33B), the Einstein Research Unit (ERU) Quantum Project (ERU-2020-607), and the European Union's s HORIZON MSCA Doctoral Networks programme and the AQTIVATE project (101072344). This work is supported with funds from the Ministry of Science, Research and Culture of the State of Brandenburg within the Centre for Quantum Technologies and Applications (CQTA). PS ac- knowledges support from ERC AdG NOQIA; MICIN/AEI (PGC20 10EE2018-01
10001 1033,, LEX2019-000910- 1000900-S110.1033000 13039/55 10001 1033, Plan National FIDEUA PID2019-10690IGB-DDO
FPI). K.-.M.. was partly supported by the Institute of Information & Communications Technology Planning & Evaluation (IITP) grants funded by the government (MSIT) (No. 2019-0-00079, Artificial Intelligence Graduate School Program, Korea University and No. 2022-0- 00984, Development of Artificial Intelligence Technology for Personalized Plug-and-Play Explanation and Verification of
Explanation).
References
[1] Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, et al. Quantum computational chemistry. Reviews of Modern Physics, 92(1), 2020.
[2] Adam Smith, M. S. Kim, Frank Pollmann, et al. Simulating quantum many-body dynamics on a current digital quantum computer. npj Quantum Inf., 5:106, 2019.
[3] Nikolaj Moll, Panagiotis Barkoutsos, Lev S Bishop, et al. Quantum optimization using variational algorithms near-term quantum devices. Quantum Science and Technology, 33):030503,, 2018.
[4] Jacob Biamonte, Peter Wittek, Nicola Pancotti, al. Quantum machine learning. Nature, 549(7671)):99- 202, 2017.
[5] John Preskill. Quantum Computing in the NISQ era and beyond. Quantum, 2:79, 2018
[6] Han-Sen Zhong, Hui Wang, Yu-Hao Deng, et Quantum computational advantage using photons. Science, 370(6523); 1460, 2020.
[7] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, et al. A variational eigenvalue solver on photonic quantum processor, Nature Communications, 5(1):4213, 2014.
[8] Jarrod R McClean, Jonathan Romero, Ryan Babbush, al. The theory of variational hyhrid quantum- classical algorithms. New Journal of Physics, 18(2):023023, 2016.
[9] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii. Quantum circuit learning. Physical Review A, 98(032309), 2018.
[10] Maria Schuld, Ville Bergholm, Christian Gogolin, et al. Evaluating analytic gradients on quantum hardware. Phys. Rev. A, 9933333333,, 2019.
[11] Nakanishi, Ken M. and Fujii, Keisuke and Todo, Synge. Sequential minimal optimization for quantum- classical hybrid algorithms. Phys. Rev. Res., 2:043158, 2020.
[12] John Platt. Sequential minimal optimization A fast algorithm for training support vector machines. Microsoft Research Technical Report, 1998.
[13] Bobak Shahriari, Kevin Swersky, Ziyu Wang, et al. Taking the human out of the loop: review of bayesian optimization. Proceedings of the IEEE, 104(1):148 175, 2016.
[14] P. Frazier. tutorial on bayesian optimization. ArXiv e-prints, 2018.
[15] Giovanni lamnelli and Karl Jansen. Noisy bayesian optimization for variational quantum eigensolvers. ArXiv e-prints, 2021.
[16] David Eriksson, Michael Pearce, Jacob Gardner, et al. Scalable global optimization via local bayesian optimization. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.

---

[17] D. J. MacKay. Introduction to gaussian processes. NATO ASI Series F Computer and Systems Sciences. 1998.
[18] Benjamin Letham, Brian Karrer, Guilherme Ottoni, et al. Constrained bayesian optimization with noisy experiments. ArXiv e-prints, 2017.
[19] Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. J. of Global Optimization, 13(4):455-492, 1998.
[20] Diana M. Negoescu, Peter I. Frazier, and Warren B. Powell. The knowledge-gradicng algorithm for sequencing experiments in drug discovery. INFORMS Journal Computing, 23(3):346-363, 2011.
[21] Peter Frazier and Jialei Wang. Bayesian Optimization for Materials Design, pages 45-75. Springer International Publishing, 2016.
[22] Eric Brochu, Vlad M. Cora, and Nando de Freitas. tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. ArXiv e-prints, 2010.
[23] Jasper Snock, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems, volume 25, Curran Associates, Inc., 2012.
[24] T.L. Lai and Herbert Robbins. Asymptotically efficient adaptive allocation rules. Adv. Appl. Math., 6(1):4-24 22, 1985.
[25] Niranjan Srinivas, Andreas Krause, Sham Kakade, et al. Gaussian process optimization in the bandit setting: No regret and experimental design. In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML'10, page (015-1022,, 2010.
[26] H. J. Kushner. A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, 86(1):97-106, 1964.
[27] Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas. The application of Bayesian methods for seeking the extremum. Towards Global Optimization, 2(117-129):2. 1978.
[28] Philipp Hennig and Christian J. Schuler. Entropy search for information-efficcen global optimization. J. Mach. Learn. Res., 13:1809- 1830-1137,2222 2012.
[29] José Miguel Henrández- Matthew W. Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. In Proceedings of the 27th International Conference on Neural Information Processing Systems Volume 1, NIPS' 14, page 918-926, 2014.
[30] Peter Frazier. Warren Powell, and Savas Dayanik. The knowledge -raddeent policy for correlated normal beliefs. INFORMS Journal Computing, 21(4):599-613, 2009.
[31] David Ginsbourger, Rodolphe Le Riche, and Laurent Carraro. A Multi Parallel Global Optimization based on Gaussian Processes, 2008.
Multi-point points Criterion for Deterministic
[32] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes Machine Learning. MIT Press, Cambridge, MA, USA, 2006.
[33] Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, volume 4. Springer, 2006.
[34] Peter Auer. Using confidence bounds for applotattion-explraataii trade-offs. Mach. Learn. Res., 3:397-422 422, 2003.
[35] S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe. Demonstration of a small programmable quantum computer with atomic qubits. Nature, 536(7614):63 66, 2016.
[36] D. Kielpinski, C. Monroe, and D. Wineland. Architecture for a large-scale ion-trap quantum computer. Nature, 41776890):099-11, 2002.
[37] Walter Rudin. Principles of Mathematical Analysis. McGraw-Hill, 1964.
[38] Manpreet Singh Jattana. Fengping Jin, Hans De Raedt, and Kristel Michielsen. Improved Variational Quantum Eigensolver Via Quasidynamical Evolution. Phys. Rev. Applied, 19(2):024047, 2023.
[39] Subir Sachdev. Quantum phase transitions. Wiley Online Library, 2007.
12

---

[40] Yasar Y Atas, Jinglei Zhang, Randy Lewis, Amin Jahanpour, Jan F Haase, and Christine A Muschik. Su hadrons on quantum computer via variational approach. Nature communications, 12(1):6499, 2021.
[41]
Lena Funcke, Tobias Hartung, Karl Jansen, al. Review quantum computing for lattice field theory. ArXiv e-prints, 2023.
[42] Alberto Di Meglio, Karl Jansen, Ivano Tavemelli, Constantia Alexandrou, Srinivasan Arunachalam, Christian W Bauer, Kenstin Borras, Stefano Carrazza, Arianna Crippa, Vincent Croft, et al. Quantum computing for high-energy physies: State of the and challenges. summary the go4hep working group. ArXiv e-prints, 2023.
[43] H. Abraham et al. Qiskit: An Open-source Framework for Quantum Computing. Zenodo, 2019.
[44] A. Nicoli, C. J. Anders, et al. Emicore: Expected maximum improvement over confident regions.
https://giis gitggttthu...... come wwweereeeceeeeeeeeeeeeeeeoeeeeee eeeemmmccorr.rr core 2023.
[45] Dong Liu and Jorge Nocedal. the limited memory bigs method for large scale optimization. Mathematical programming, 45(1-3):503-528, 1989.
[46] Trevor Hastie, Robert Tibshirani, Jerome Friedman, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
[47] Zhenyu Cai, Ryan Babbush, Simon C. Benjamin, Suguru Endo, William J. Huggins, Ying Li, Jarrod R. McClean, and Thomas E. O'Brien. Quantum Error Mitigation. ArXiv e-prints, 2022.
[48] Mateusz Ostaszewski, Lea M. Trenkwalder, Wojeiech Masarczyk, Fleanor Scerri, and Vedran Dunjko. Reinforcement learning for optimization of variational quantum circuit architectures. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 18182-18194, Curran Associates, Inc., 2021.
[49] Jin-Guo Liu, Liang Mao, Pan Zhang, and Lei Wang. Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits. Machine Learning: Science and Technology, 2(2)022001, 2021.
[50] Hamed Mohammadbagherppor, Patrick Dreher, Mohannad Ibrahim, et al. Exploring Airline Gate- Scheduling Optimization Using Quantum Computers, ArXiv e-prints, 2021.
[51] Yahui Chai, Lena Funcke, Tobias Hartung, et al. Towards Finding an Optimal Flight Gate Assignment on a Digital Quantum Computer, ArXiv e-prints, 2023.
[52] Peter Deglmann, Ansgar Schäfer, and Christian Lennartz. Application quantum calculations the chemical industry an overview. International Journal of Quantum Chemistry, 115(3):107 136, 2015.
[53] Yudong Cao. Jonathan Romero, Jonathan Olson, et al. Quantum chemistry in the age of quantum computing. Chemical reviews, 11911))11855-10015, 2019.
[54] Alexander McCaskey, Zachary Parks, Jacek Jakowski, al. Quantum chemistry as a henchmark for near-term quantum computers. npj Quantum Information, 5(1):99, 2019.
[55] Yudong Cao, Jhonathan Romero, and Alán Aspuru-Guzik. Potential of quantum computing for drug discovery. IBM Journal of Research and Development, 62(6))-11, 2018.
[56] Alexander Heifetz. Quantum mechanics in drug discovery. Springer, 2020.
[57] Nick S Blunt, Joan Camps, Ophelia Crawford, et al. Perspective on the current state-of-the-art quantum computing for drug discovery applications. Journal of Chemical Theory and Compudation, 18(12):7000 7023, 2022.
[58] Joris Kattemölle and Jasper van Wezel. Variational quantum eigensolver for the heisenberg antiferromagnet on the kagome lattice. Phys. Rev. B, 106:214429, 2022.
[59] Vincenzo Lordi and John M. Nichol, Advances and opportunities in materials science for scalable quantum computing. MRS Bulletin, 46(7):589-595, 2021.
[60] Mari Carmen Banuls, Rainer Blatt, Jacopo Catani, et al. Simulating lattice gauge theories within quantum technologics. The European physical journal D, 74:1 2020.
[61] Anthony N Ciavarella and Ivan A Chernyshev. Preparation of the su (3) lattice yang-mills vacuum with variational quantum methods. Phys. Rev. D. 105(7):074504, 2022.
13

---

[62] Jinglei Zhang. Ryan Ferguson. Stefan Kühn, et al. Simulating gauge theories with variational quantum eigensolvers in superconducting microwave cavities. ArXiv e-prints, 2021.
[63] Alain Delgado, Juan Miguel Arrazola, Soran Jahangiri, et al. Variational quantum algorithm for molecular geometry optimization. Phys. Rev. A, 104:052402, 2021.
[64] Bo Peng and Karol Kowalski. Mapping renormalized coupled cluster methods to quantum computers through compact unitary representation of nonunitary operators. Phys. Rev. Res., 4:043172, 2022.
[65] Lena Funcke, Tobias Hartung. Karl Jansen, et al. Dimensional expressivity analysis of parametric quantum circuits. Quamtum, 5:422, 2020.
[66] Eric R. Anschuetz and Bobak T. Kiani. Quantum variational algorithms are swamped with traps. Nature Communications, 13(1):7760, 2022.
[67] David JC MacKay et al. Introduction to gaussian processes. NATO ASI series F computer and systems sciences, 168:133-1666 1998.
[68] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, others, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020.
14

---

Limitations
One significant limitation of the approach presented in this paper, particularly in the context of Variational Quantum Eigensolvers (VQEs), relates to the scalability of Gaussian Processes (GPs). When a large number of points added to the GP training set through additional observations, the computational scalability becomes challenge, especially in scenarios involving large number of observations. However, we consider a potential solution to address this issue by imposing a fixed limit on the training sample size. This approach involves removing previously observed points and replacing them with newer ones. We hypothesize that by leveraging the information from the CoRe, the newly added points would contain significantly more valuable information, making previous observations less informative. Consequently, removing those points from the training set would mitigate the inherent scalability problem associated with GPs. Exploring this idea further is an avenue for future research. In addition, in the current version of our proposed NET--Ri we limited the choice of new observation points to two points along a sequentially chosen axis, which is clearly sub-optimal. We will extend our approach for more flexible choices, e.g., by including sets of points along different directions and different numbers of points to the candidate set, in future work.
Another limitation is related to the current constraints of VQEs and Noisy Intermediate- Scalle Quantum (NISQ) devices. The execution of quantum computations on NISQ devices is currently restricted, in particular by the coherence time of the qubits and the number of operations required to execute the algorithm. Consequently, the mcasurements on a quantum computer are susceptible to errors, which is recognized as one of the most challenging obstacles in the field. Although error mitigation techniques have been proposed [47), developing hybrid are more resilient to the inherent noise of quantum computers remains an open area of research.
algorithms that
Broader Impact
This work is remarkably successful example of the synergistic combination of theories and tech- niques developed in physics and machine learning communities. Namely, the strong physical prior knowledge of the VQE objective function is well incorporated in the kernel machine, leading to our novel VQE kernel, while specialized sequential minimal optimization-the the NFT algorithm- developed in the physies community is adapted to the Bayesian Optimization (BO) framework, mitigating the suboptimality of NFT and the scalability issue (in terms of the search domain dimen- sionality) of BO in a complementary manner. Our novel EMICoRe acquisition function plays an important role: it exploits the correlations between the function values on points, which are not necessarily in the neighborhood of each other, by leveraging a new concept of confident regions. This technique can be applied to general optimization problems where the prior knowledge implies such non-trivial correlations. In addition, the equivalence of the parameter shift rule and the VQE function form, which we have proved without using any physics-specife knowledge, can facilitate theoretical developments both in physics and machine learning communities.
Regarding the societal impact, the authors have thoroughly considered the potential negative conse- quences of the proposed work and have found none.
A Extended Related Work
Since the VQE protocol was first introduced 77. many optimization algorithms have been proposed for minimizing the VQE objective. For grad optimization, the parameter shift rule allows for an efficient gradient computation on NISQ devices 9 10]. Making use of the analytic form of the gradient in typical parametric ansatz circuits, this approach avoids estimating the gradient using finite differences, which would be challenging for NISQ hardware due to the limited accuracy that can be achieved due to noise.
The Nakanishi-Fun Fujih--FFoTooooo (NFT) method [1 harnesses the specific function-form of the VQE objective to establish Sequential Minimal Optimization (SMO) and showed the state-of-the-as performance. The authors focused on the case where the parametric gates are of the form Uxxa)= = exp(-ir, with angular parameters E RD and operators P fulfilling p2 = I and derived an explicit function-form of the VQE objective. The resulting function form implies that, by keeping all parameters (angles) in the circuit fixed except for single one, can identify the complete
15

---

objective function in the corresponding one onedimmesiional subspace by observing only three points, and the global minimum in the subspace can be analytically found. NFT uses this property and performs SMO by choosing one-dimensional subspace sequentially randomly until convergence, providing an efficient and stable algorithm suited for NISQ devices.
BO is a versatile tool for black-box optimization with its applications including engineering system
design [9), drug design [20], material design [2I] and reinforcement learning [22]. Recently,
it started to be extensively used for hyperparameter tuning of deep neural networks (23) Most
work on BO uses the GP regression model, computes an acquisition function, which evaluates the promissing- -nsss of the next candidate points, and suggests its maximizer as the set of next observation points. Many acquisition functions have been proposed. Lower (upper for maximization problems) confidence bound 244 25 optimistically gives high acquisition values at the points with low predictive mean and high uncertainty. Probability of improvement (26] and expected improvement (EI) [27 34
evaluate the probability and the expectation value that the next point can improve the previous optimum. Entropy search [28] [9] searches the point where the entropy at the minimizer is expected to be minimized. Knowledge gradient (30) allows final uncertainty and estimates the improvement of the optimum of the predictions before and after the next sample is included in the training data. The most common acquisition function is EI, and many generalizations have been proposed. Noisy EI (NEI) [18] considers the observation noise and takes the correlations between observations into account, parallel EI [3] considers the case where a batch of new samples are to be suggested, and EI per cost (or per second only the computation time matters) 23 penalizes the acquisition function value based on the (estimated) observation cost.
BO has also been applied to VQE minimization [5]. was shown that combining the periodic kernel 117 and NEI acquisition function [8] significantly improves the performance of BO with the plain RBF kernel and EI, thus making BO comparable to the state-of-the-art methods in the regime of small qubits and high observation noise. Our approach with Expected Maximum Improvement over Confident Regions (EMICoRe) has similarities to existing methods and can be seen as a generalization of them. The key novelty is the introduction of Core Regions (CoRe), which defines the indirectly observed points. Note the difference between the trust region [16] and CoRe: Based on the predictive uncertainty, the former restricts the regions to be explored, while the latter expands the observed points.
For completeness, we note that other machine learning techniques from reinforcement learning 48 and deep generative models [4 have also been applied to improve the classical optimization schemes of VQEs.
B Details of Gaussian Process (GP) Regression and Bayesian Optimization (BO)
In the following, we introduce GP, GP regression, and BO with an uncompressed notation.
B.1 Gaussian Process Regression
A GP 32 is an infinitediimensional generalization of multivariate Gaussian distribution. Let : R be a random function, and denote the density of GP GP(S(); k(:-)) where v(.) and -). are the mean function and the kernel (covariance) function, respectively. Intuitively, stating that a random function ʃ(-) follows GP, i.e. p(f(-)) ()) k(,,))) )), means that the function values f(x) indexed by the continuous input variable x E X' follow the diminie-immesionall version (i.e., process) of the Gaussian distribution. The marginalization property [3] of the Gaussian allows the following definition:
Definition 1. (Gaussian process) GP is the process of a random function such that, if f(-) GP(((-), k(-,)), -)), then any set of input points X = (æ1..... xN) E X it holds that
p(fly, K)= = -ND(FU, K),
(14)
16

---

where
f (f(z1)..... f(x,). E E.,,
RR,,
k(m k()) k=
ERNXN
21) x(
(2 Pei...,,''1! E XM, and let f' & TR,, v' (v(a(),.... ) E TR be the corresponding random function values and the mean function values, respectively, Then, Definition Timpless that the joint distribution of f and f is
Consider another set of input points x'
p(f.f) = NN+ta K)
(15)
where
( TT ERN+M
;=v TTRNM
K' K K"T K") ER(N+M)((N4M0
Here, = ( X) E RNXN K' x(X,') E N and K" = =(X', X') E RMM
where x(X,X') )') denotes the kernel matrix evaluated at each column of X and X' such that (k(X, X") x))),, = k( Tm).
The conditional distribution f' given f can be analytically derived as
p(f f) = Nas (f'; Meand Sooad),
(16)
where
Mondd = + K"T K Kf E RM
S = = KK
E ERM.
In GP regression, X and X' correspond to the training and the test inputs, respectively. The basic idea is to use the joint distribution (T3) the prior distribution on the training and the test points, and transform the likelihood information from f to f' by using the conditional (16).
The GP regression model consists of the Gaussian noise likelihood and GP prior:
plylæ, f()) (y; (z))
p(f()) (G)())))())M)))))
(17)
where x² denotes the observation noise variance. Below, assume that the prior mean function is the constant zero function, i.c., v(x)= 0, Væ. Derivations for the general case can be obtained by re-defining the observation and the random function as t y v(x) and f(x) f(x) v(x). respectively.
Given the training inputs and outputs, X (#1, ........ xN) N and y (y1....,, ,NN E RRN. the posterior of the function values at the training input points f = (()))...,)(( E ERN is given
p(f(,y)
= NN.Ji
(18)
plg ply(X)
where
#1 -= -- (K (K-1 INN BERN
S S II) NxN
Given the test inputs X' = (al n'y)) E X the posterior of the function values at the test f
points
E ERM can be obtained, by using Eqs. 16 and (18) as
p(f'Xx,) 4) /(fIFII(AIX 9)
(19)
where
u = R.T (K + -IIINN DERM, -' K" (K+ IK+ K' E MxM (20)
17

---

The predictive distribution of the output y' - (f(x() + + +EM) E ERM is given as
plu' y) p(u'lf' {, = S'),
(21)
where
11 = 1 E ERMM =' + o E MRM
The marginal distribution of the training outputs is also analytically derived:
p(y]X) = f
Mmarg Smarg):
(22)
where
OEER
S marg = o KERNNN
Marg
The marginal likelihood
22
used for hyperparameter optimization.
B.2 Bayesian Optimization
In BO 14]. a surrogate function, which in most cases is GP regression, equipped with uncertainty estimation, is learned from the currently available observations. new set of points that likely improves the current best score is observed in each iteration. Assume that at the f-th iteration of BO, we have already observed points X E X BO suggests new set of M points X' E X by solving the following problem:
max (X'), X
where (x(+) is an acquisition function computed based on the GP trained on the observations y at X. A popular choice for the acquisition function is Expected Improvement (EI) [27 34
(x') (max(0, f p(PX)))))
which covers the case where the observation noise is negligible and only a single point x' is chosen in each iteration, i.e., 0² 1, 1. Here, f denotes the current best observation, i.e.f = mi...e(.....) (ax,,, ())p denotes the expectation value with respect the distribution p. and ( y) is the posterior distribution (19) of the function value f' at the point x'. El can be analytically computed:
( (
(-_1'x 11x
= (f
+ o
(23)
8'x
where ''x E E and s'x E R are the GP posterior mean and variance indicating the input points on which the GP was trained, and
(20) with their subscripts
Φ(E) = J
Φ(E) 0,
((=)=
(e; 13)
are the probability density function (PDF) and the cumulative distribution function (CDF), respec- tively, the one-dimensional standard Gaussian.
For the general case where o² > M 2 1, Noisy Expected Improvement (NEI) proposed:
18
31 was
NEI (X') = (max (0, min(fnn miin(mm))) X(
(24)
NEI treats the function values f at the already observed points X still random variables, and appropriately takes the correlations between all pairs of old and new points into account. This is beneficial because it can avoid overestimating the expected improvements at, e.g., points close to the current best point and set of two close new points in promising regions. A downside that NEI does not have analytic form, and requires quasi-Monte Carlo sampling for estimation. Moreover, maximization is not straightforward, and an approximate solution is obtained by sequentially adding point to X' until M points are collected.
18

---

Solution
Init. State
12/0)
Var. WF
Updale Gix)
Measure
VQE
Params.
H
x
x -
Find Optim
(shx) B(x)
Task
Figure 5: Illustration of the VQE workflow. In the first step, highly complicated optimization problems can be translated into the Hamiltonian formulation (see, e.g.. (50]51)) The Hamiltonian H and an initial state (do) are plugged into the VQE block (light purple), where the variational quantum circuit is instantiated with random angular parameters x". In the VQE block, the top vertex of the triangle represents the quantum computer, the red arrows and blocks refer to operations running quantum computer, while the green parts refer to classical steps. In the bottom-left green box, the current parameters x are updated with the new best parameters x found during classical optimization routines. Then, the quantum circuit G(x) is updated using the new optimum point, + x, and the energy Exx) for the updated variational wave function |xx) is measured. The VQE block executed for 7 iterations and finally outputs the solution to the task as a variational approximation ((*r) of the ground state and the corresponding optimal parameters x" for the quantum circuit.
Our EMICoRe, proposed in Section 3.2 can be seen generalization of NEL In EMICoRe, the points in the confident regions (CoRe), where the predictive uncertainty after new points would have been observed is lower than threshold searched for over CoRe. If we replace CoRe with the previous and the new training points, EMICoRe reduces NEI.
treated as "indirectly observed", and the best score is
Another related method to our approach is Knowledge Gradient (KG) [3
KG (X') = mmxx (0. min 4x (z") min x(xxx')(x)
((xinn m"EX
(25)
"""x"
which assumes that the minimizer of the GP posterior mean function is used as the best score even the uncertainty the minimizer high-and and estimates the improvement of the best score before and after the new points X' are added to the training data. The second term in Eq. (25) is estimated by simulation: trains the GP on the augmented data (X, X') and (y yy UTT)T. where y' are drawn from the GP posterior p(y'l X, y). and finds the minimizer of the updated GP mean. Iterating this
process provides Monte Carlo samples to estimate the second term in Eq. (25).
In our EMICoRe method, if we set the CoRe threshold x2 + oo so that the entire search domain is in CoRe, and replace the random function f(-) Eq. (TT) with its previous and updated GP means, respectively, EMICoRe reduces to KG. Thus, KG can be seen as version of EMICoRe that ignores the uncertainty of the updated GP.
C Details of Variational Quantum Eigensolvers (VQEs)
The VQE 7 g is hybrid quantum-classica algorithm that uses a classical optimization routine in conjunction with a quantum processor to approximate the ground state of a given Hamiltonian. VQEs are designed to current generation of quantum computers. These devices have a limited number of qubits and high these
Noisy Intermediate-Scaal Quantum (NISQ) devices, which are the
error rates, which makes it challenging to execute complex quantum algorithms faithfully quantum devices. VQEs are a promising approach to tackle this challenge since they use a hybrid dlasicall-uuntunr algorithm, where the quantum device only needs perform relatively simple computations (see Fig. 5 for an illustration the VQE workflow).
VQEs are considered to be potentially relevant for some challenging problems in different scientific domains such as quantum chemistry [52 54 drug discovery (55]
condensed matter
physics (58), materials science [59] and quantum field theories 660 6 62 41 Specifically, finding
19

---

the ground state of a molecule is a very challenging problem growing exponentially in complexity with the number of atoms for classical computing, while for VQE this would instead scale polynomially, Despite being naturally designed to solve problems associated with quantum chemistry and physics, such as calculating molecular energies, optimizing molecular geometries [63]. and simulating strongly correlated systems [64], VQEs have also been applied to other domains including optimization and
combinatorial problems such as the flight gate assignment [50 51].
Given Hamiltonian formulation that can be efficiently measured on quantum device, the variational approach of VQE can be applied to obtain an upper bound for the ground-state energy as well as an approximation for the ground-state wave function. Let (&o) be the initial ansatz for the Q- dimensional (qubit) wave function. Let us assume that we use a parametrized quantum circuit G(z)).
where x E [0, 2n)² represents the angular parameters of quantum gates. The circuit G consists of (2 D) unitary gates:
G(x) GD o.... GI,
(26)
where D of the D' gates depend on one of the angular parameters exclusively, i.e., Fdd parametrizes only a single gate Gar(d) (xz), where d' (d) specifies the gate parametrized by xz We consider the parametric gates of the form
ViIl
J Ga - Ua) = exp { {
(27)
where Pa is an arbitrary sequence of the Pauli operators (o²² og og fq=1 acting on each qubit at most once. This form covers not only single-qubit gates such as Rx(x)) exp (-i0or' 8 ). but also entangling gates such as Rxx(x) exp (-ixo e1 ) and Rzz(2) exp (-iro mn o for q1 # 92. In the matrix representation of quantum mechanics, quantum states are expressed as vectors the computational basis, i.e.,
10) = ()
[1] = (i)
(0) = (1 0),
(0 1),
(28)
Moreover, Pauli operators are expressed as matrices,
xx ( .
0
o² (1 d-(1-) -
acting only on one of the Q qubits non-trivially. The application of any operator quantum state (2)) thus becomes a matrix multiplication in the chosen basis. Single-qubit parametric gates correspond to rotations around the axes in a three-dimensional space representation of a qubit state, known as the Bloch sphere, and are widely used in ansatz circuits for VQEs.
Given the Hamiltonian H, which is an Hermitian operator, the quantum device measures the energy of the resulting quantum state (s>a) G(z)) voo) contaminated with observation noise E, i.e. (29)
f(x) = (xx + where (xx = (Vizl | (())) (a)) HG(=)|Vo).
The observation noise E in our numerical experiments only incorporates the shot noise coming from the intrinsically probabilistic nature of quantum measurements, and the errors from imperfect qubits, gates, and measurements on current NISQ devices are not considered.
The task of the classical computer in VQE is find the optimal angular parameters æ such that f=) is minimal. Given that the ansatz circuit is expressive enough, G(x)) (1)))) then corresponds to the ground state of the Hamiltonian H. Thus, the problem to be solved is noisy black-box optimization:
min @0(00,2m) D
f (z).
The long-term goal of research on VQEs to develop efficient quantum algorithms that can solve problems beyond the capabilities of classical computers. While VQE is a promising approach for exploiting the advantages of quantum computing, they are plagued by some limitations. Specifically, these algorithms require appropriate choices for the quantum circuit [65, [6). To have compact circuits with a moderate number quantum gates is, therefore, essential to favor stable computations and high measurement accuracy.
Vil In Appendix E we discuss how the theorems and our method can be extended to the non- non-exclusive
parametrization case.
20

---

(0)
Ry (xo)
(x(25)
Ry (26)
(x
(212)
Ry (215)
(0)
(x1
(z(
(y)
Rz (210)
Ry (x
(
(0)
Ry (22)
Rx(2s)
O
1
Ry (ax)
Rx(211)
11
O
O
Ry (214)
Rz(217)
Efficient (U) Circuit with default parameters for Q 3 qubits and I. 2 layers (thus a E [0, 20)) Each dashed box indicates a layer of the ansatz circuit. The quantum computation proceeds from left to right with an initial ansatz of the form |0) 8 Each horizontal line corresponds
Figure 6: Illustration of Qiskit's
43
quantum wire representing one qubit. The symbols the wires correspond to gate operations on the respective qubits, starting with two parametrized rotational gates Ry and Rz acting the qubits in the initial 0-th layer. Then, each the following layers is composed of one block of CNOT gates and two rotational gates acting on each qubit.
Efficient SU(2) circuit: One circuit ansatz commonly used is the Efficient SU(2) Circuit from Qiskit [43], which is illustrated in Figure this circuit, two consecutive rotational gates,
(x) exx (-1 I
(-r
x
and Rz(x) exp
(30)
act each qubit in the initial (0-th) layer. The rest of the circuit is composed by stack of L layers, each of which consists of CNOT gates applied to all pairs of qubits, and pair of the rotational gates, parametrized and therefore not updated
Ry and Rz, acting each qubit. The CNOT gates during the optimization process. Therefore, the total number of angular parameters of the circuit is equal the number of Ry and Rz gates in the circuit. total, for setup having Q qubits and layers, the number of angular parameters is
D - 2 Q + (L 2) x Q,
(31)
where the first term counts the rotational gates in the initial layer, while the second term counts those in the latter layers. In experiments, set the initial ansatz (bbo) |oy e which corresponds to the tensor product of Q qubits the (0)-ket state VIII
D Proof of Theorem
Proof. The VQE kernel can be rewritten
VOR (z, x') = o?(?2 2) II (x² + 2cos(Ea - 2'd))
d=1
= = + -
D II (?² + 2 cos cos x'a + 2 sin x. sin 2//)
d=1
0g62 +2 2)
Σ (E(0.12))
( (2 x. cos x'd) x'd) ((t==1) (2sin a sin d'y ((2==))
o<x) ((('),
where 1()) denotes the indicator function (equal to one if the event is true and zero otherwise), and
Φ(x) 00(?2 + 2) vec (8 cos za, 22 xa)" R
which completes the proof.
E
Generalization to Non-Exclusive Parametrization Case
In the VQE, it often beneficial to share the same parameter amongst multiple gates, for example, in the case where the Hamiltonian has certain symmetry, such translation invariance. Doing so,
VIII An equivalent notation often found in the literature is |0) 9 0).
21

---

the variational quantum circuit is guaranteed to generate quantum states that fulfill the symmetry, and thus the number of parameters be optimized is efficiently reduced. To consider this case, we need to assume that some of the entries of the search variable a: are shared parameters, each of which parametrizes multiple gates. The Pakanishi-Fujii- Todo (NFT) algorithm
[1]] can still used in this
case, based on the following generalization of Proposition 2
Proposition 3. (1TV Assume that the d-th entry of the input E [0, (T) parametrizes Va 2 1 gate parameters. Then, for the VQE objective function (+) in Eq. (6)
cos x.
36 [R (1+2V))))) such that [=)x) = b vec
(Vara) sin ad
(32)
sin (Vaza))))
Similarly, our VQE kernel introduced in Theorem can be generalized as follows:
Theorem 3. The (higher-order VQE kernel.
7 72 Va com(u(ar - m'2)
V (x,x') x') = o 11
(33)
72+2Va
is decomposed as V2 (2 (xx))(((x') where
v22 cos xd
Φ(x) 0o (7² + 2Va) DD/2
vec
V2cos(Vazal Visin ad
(34)
22
Proof. Similarly to the exclusive parameterization (or first-order) case, have
72 2 V
cos (v(x, x'd))
VQU (xx x²) = =0 II
72 + 2V
d=1
= on + 2Vd)
2 +
(cos(vra) cos(uz's) + sin(vza) sin(vz'a)}
d=1
= o + 2Va) DD
V4
Σ (63
II
(2 cos(uera) cos(vx'd)) (Eu= 20--) (2sin(war))) sin(wa())) ( (Ea=20)
2V(0....2V.) d=1 ((x) (x)'),
which completes the proof.
Since the VQE objective 32) is the Va-th order sinusoidal function along the rx-axis, NFT can perform the sequential minimal optimization (SMO) by observing 2V, points in each step - which become 2 + points, together with the current optimum point to determine the entire function form the one-dimensional subspace parallel to the ax-axis. Our NFT--wit-EEMICORe approach can similarly be generalized to the non-exclusive cases: for the chosen direction d, the approach observes M = 2Va new points by maximizing the EMICoRe acquisition function, based on the GP regression with the generalized VQE kernel 33
22

---

Algorithm 1: Nakanishi-Fuji-Tod (NFT) method [11] (Baseline)
input
TMI max # of iterations
TRI reset interval
D (+ (9) initialization with x ~ [0,2n] 2x) and 1 = f²(2") +E
output:
i Tha last optimal point
y o T last estimated objective
Tau ensemble of collected observations
begin
for 1 t= to TMr do
Choose direction d E (1, {I...., D) sequentially or randomly;
Find (X') X' (x' x222 ( 2T 2T/ /e, a'-1 2T 23/3 along d:
Observe y/ = y/ (y/ 1/2 at the new points X';
Append(D (X', 1²)) = Dr:
Fit( = [(0)
+ cı cos 0 sin 0 to the three points
{(-2m/3, v1), (0, "), (2T/3. 1/2)}:
FindMin (f) = find analytical minimum θ ar f(0):
Update è with a' _'-1 θeg:
Update (j) = with estimated e = f(0);
10
11
ift mod TRI- then
12
Observe y(="):
/* Perform additional observation */
D' D' u (#', v(*)));
13
14
end
end
end
return DTa a Thi y(
17
F
Algorithm Details
Here, we provide the detailed procedures of the baseline method (Nakanishi-Fuii-Ton (NFT) 11] and our proposed method LNFT-with-EMICCRe and the EMICoRe subroutine).
F. Nakanishi-Fuji-Tood (NFT)
In each step of NFT (Algorithm I axis d (1.....D) is chosen sequentially or randomly, and the next observation points along the axis are set (Step 4 and observed (Step (5). Based on Proposition[] the two new observations together with the previous optimum are fitted to a sinusoidal function, the global optimum along axis analytically computed, establishing sequential
minimal optimization (SMO) [12]. NFT iterates this process from a random initial point and outputs the collected datapoints and the last optimum. Since the optimal point at each step observed, enrors can accumulate over iterations. As remedy, NFT observes the optimal point when the reset interval condition met (Step
directly
F.2 NFT-with-EMICoRe
Our proposed method, NFT--wih-EMMCCRR (Algorithm) 2). replaces the deterministic choice of the next observation points X' in NET with promising choice by BO. After initial (optional) iterations of NET until sufficient training data are collected, we start the EMICoRe subroutine (Algorithm 3) to suggest observation points X' by BO (Step
8 Then, the objective values y' X' are observed 1 updated with the new observed data
(Step 10), and the training data D
{
23

---

Algorithm 2: NFT-wihh-EMICRR
input
TMi of iterations
TNFT # of initial NFT steps
Tave(> TNFT): averaging steps for update.
D (# ( initialization with a. [0, 2x) and 1o f² +
output:
a Thn last optimal point
/ r( F last estimated objective
Divai ensemble of collected observations
TNTT > then
DKFTT (I TNFT, TRI 0, D)
/* see Algorithm 1 Update observations Dzhor DD: /* Collect points from NET step */
*/
end
begin
for 1 TNPT + to Tir do
Choose direction d E {I....,D} D) sequentially or randomly;
Find X = points maximizing the EMICORe($' D d' k*)
/* for EMICoRe sub-routine
Algorithh
Observe y/ = y' = v6) T at the new points X'; Append(D' (X', y²)) - Dr:
10
11
Train GP(D") updated dataset; Compute posterior means μ = (μ(#²-1 2n 23ea), (a (u("==("))))) u(a'') -6'"-11___ 2 T:
13
Fit(f(0)) f(0) + cı cos 0 02 sin 0 to the three points
14
{(-27/3 / p1): P2), (2m/3, wa)}
FindMin (f) = find analytical min 9 = argmingg(o.24 f(0):
15
Update è = with a' _'-1 θeg:
16
17
Evaluate optimal objective ji' = p((*); i 2 TAve then
18
Compute the CoRe threshold for the next iteration: k*+1 2-- Tam-
19
Tave
20
end
end
end
23 return MM, Thi u(t
(X', y'} (Step 11p. The GP updated with the updated training data D = [x", y²), and its mean predictions at three points are fitted by sinusoidal function (Step 14p for finding the new optimal point ø' (Step 15 Before going the next iteration, the CoRe threshold is updated according to the energy decrease in the last iterations (Step 19). In the early stage of the optimization, where the energy ji' decreases steeply, a large encourages crude optimization, while in the converging phase where the energy decreases slowly, a small x enforces accurate optimization steps.
E.3 EMICoRe Subroutine
The EMICoRe subroutine (Algorithm receives the current optimal point å the current training data (X, y). the direction to be explored, and the CoRe threshold k, and returns the suggested observa- tion points. Fixing the number of new observations per step to M 2, we prepare pairs of candidate
points (sampled grid) along the axis d as candidate set C { B Dx2 Jac(Jsa-1) -1) (Step 2).
j==

---

Algorithm 3: EMICoRe subroutine
input
à current optimal point
D (X,y) y current training data
d direction to be explored
K CoRe threshold
params:
M = of suggested points
JsG # of search grid points
Joa # of evaluation grid points
NMC # of Monte Carlo samples
output
X'; suggested observation points
begin
Prepare a candidate set C = (X' E Dx21 DDx))Jaa(Js=_)) f==1
/* with JsG (Jsa - 1) being the # of candidate pairs
for j = to Jsa (JsG 1); 1);
/* denotes one pair of points */
do
Update GP adding the current candidate point GP(X) where X = (X, X')
Compute the posterior variance x (x, x) Yæ in the evaluation grid, along axis d:
/* test GP uncertainty on evaluation points.
Find discrete approximation of the CoRe as Z = fo E X Grid
x (x, - n²);
Use & and X test
2x compute mean and the covariance of GP posterior:
p(f.f plest x. y):
10
Estimate acquisition function by quasi Monte Carlo sampling
EMICoRe (j) ad (max (0, m }}) pt prom xx,)
11
"
12
end
Find pair of observation points that maximizes EMICoRe 13
j argmax, d [EMCCole(j).
end
return Suggested points to observe at next step X' X
For each candidate pair, the predictive variances of the updated GP are computed on æ E t points on test grid, along the direction d (Step This way, a discrete approximation of the CoRe obtained by collecting
grid points where variance is smaller than the threshold x² (Step R).
After computing the mean and the covariance of the current GP at the current optimum & and on the (discrete) CoRe points (Step - which is D-dimensional Gaussian for D l = =1+1 the EMICoRe acquisition Tunction estimated by quasi-Monte Carlo sampling (Step IT. After evaluating the acquisition functions for all candidate pairs of points, the best pair is returned as the suggested observation points (Step T6).
-wich
F4 Parameter Setting
We automatically tune the sensitive parameters: the kernel smoothness parameter y is optimized by maximizing the marginal likelihood in each iteration in the early phase of optimization, and at intervals in the later phase (see Appendix H for the concrete schedule in each experiment). The CoRe threshold & is updated at every step and set to the average energy decrease of the last iterations as in Step 19 of Algorithm 2 which performs comparably to the best heuristic in our investigation
below, see Table |
The noise variance o² is set by observing f (xx) several times at several random
25

---

Table 1: Performance of EMICoRe depending on the choice of hyperparameters Co and C₁ for the CoRe threshold update rule (35). The best results are highlighted in bold, while (1) indicates whether lower (higher) values are better.
Description
Co
C1
Energy +
Fidelity 1
Default (Eq. 12)
0.0
1.0
-5.82 + 0.14
0.85 + 0.16
Extreme (small)
0.1
0.1
-5.22 = 0.11
0.85 = 0.16
High (large)
10.0
10.0
-5.72
0.15
0.82 = 0.16
Extreme (large)
10.0
100.0
-5.70 =: 0.16
0.80 = 0.18
Best
10.0
-5.84 0.09
0.87 + 0.11
points and estimating a² o- 2 (Nshots) before starting the optimization. For the GP prior, the zero mean function (xx) = 0, Va is used, and the prior variance of is roughly set so that the standard deviation do is the same order as the absolute value of the ground-state energy. The parameters MMI JsG, JoG, and NMC should be set to sufficiently large values as long as the (classical) computation is tractable, and the performance is not sensitive TNFT and Tave. The parameter values used in our experiments are given in Appendix
Investigation of CoRe Threshold Update Rules:
In our experiments in Section 4 the CoRe threshold is updated by Eq, (I2). Here, we investigate whether a more fine-tuned update rule can improve the performance. Specifically, test the following protocol:
1+1 max (Co-o, C1
ft TAve f
(35)
T TAve
where Co, Ci 2 0 are the hyperparameters controlling the lower bound and the scaling of the average energy reduction, respectively. We note that a is the standard deviation of the observation noise, and
setting the hyperparameters to Co 0, Ci = 1.0 reduces to the default update rule (12) Table [ shows the achieved energy and fidelity for different values of the hyperparameters C. and C1. after 600 observed points, in the setting of the Ising Hamiltonian at criticality and (L 3) 3)-layered
(Q = quubisss quantum circuit with Nahots = 1024 readout shots. We observe that choosing C. 0.1 and C 10 leads to the best performance; however, we also note that the setting used for the paper 12 achieves a similar performance.
G Proof of Theorem
Proof.
We divide the proof into two steps.
G.1
Eq.
Eq.
8
The parameter shift rule
7 for 1/2 gives
(x)==f- 2 Oua 2-)
) 1"
(x x π 2
Væ E [0,2x) d I,...,,, (36)
which implies the differentiability of fx(x) in the whole domain [0, 2m) For any 1..... D and 2 [0, 2Tr)), consider the one-diment dieeeninnsionl subspace of the domain such that Ad,2 { aee [0, ). and the following restriction of (-) on the subspace:
fa = ='' (a) ( + (xra fa)ea)-
For this restricted function, the parameter shift rule (36) applies as ② Oxx Fia Fā (aa
1 7) fā. (xa "f 7).
Vxa E (0,2 dl D.
(37)
The periodicity of (-) requires that Jā (Ra) can be written as a Fourier series,
Fr.( = (2 + Σ (ca,1, cos (Txz)) + cap,r (zxa) sin (Txa))
(38)
26

---

where (axa) denote the Fourier coefficients, which depend on æ except for dd. Below, we omit the dependence of the Fourier coefficients on ad to avoid cluttering.
Substituting Eq. 38 into the left- and the right-hand sides of Eq. 37 respectively, gives
2
(Td)
2 Σ
T ( cd,1,r sin (TIJ) + cd.2,r cos (Tæa)),
(39)
Ord
==1
06
fà
π Td +
②)
fá.
Ed
Σ
π Cd,1,7 cos (T )) - cos (T (x
2
sin (T T ) sin (T (Ia - π
+ Cd22
Cd.1,r sin sin
cd.2,7 cos (TTAa) sin
7
2
sin
TW'
Cd,,, sin (Tea) + cd.2,r cos (TTA)).
(40)
Since Eq.
37) requires that Eqs. and 440) are equal to each other for any [0, 2T) and
39
D.
must hold that
TT
sin
VT such that ca,1,r # 0 or Cd,2,r # O.
(41)
2
Since Eq. 441 can hold only for = 1, we deduce that
Cd,1,r Cd,2,r 0, VT + 1, I,..., D.
Therefore, the restricted function must be the first-order sinusoidal function: (a,)( = (). + (a.1, As the most general function form that satisfies Eq. 42 for all d I,..., D, we have
cos (Fa) + (f.,,,(f,,.)) sin (ra))
(42)
(x)
Σ be II
1 (S
((S==1) (sin
d=1
b vec (&L (1, cos
va, cas a
Here, E E D takes the value of either 0, or 2, specifying the dependence on consanttnnt, cosine, or sinc-for each entry, and b (be (E 1,2) is the 3D -dimensional coefficient vector indexed by &. With the appropriate bijective mapping : {0, 2 1, ..... 3D consistent with the definition of the vectorization operator vec(-), we defined b R such that bi(e) bg. 1(-) the indicator function, which is equal to one if the event is true and zero otherwise.
G.2 Eq. 8
Eq.
For any (+) in the form of Eq. (8) the left- and right-hand sides of the parameter shift rule (36) can be, respectively, written as
f(x)
(2c--a1, ca. (2 (2) sin (Ir)) + ca,2,1 (10ga)) cos (zra),
(43)
Frd
8r + π Cd) 1r
c (n) { cos (ea + 2 7) cos Id 2
7) }
ed)
2
2
a sin Id 2 }} c..,)1 { sin (Ia + 2
m
2 c..11111 sin (Ia) sin 5) + cd.2, (z) a) cos (Iz) sin
2
2 2- (z) sin (xa) + Ca,2,1 1 (ar cos (Ira)
(44)
which coincide with each other. This completes the proof.
27

---

Table Choice of coupling parameters for the Ising and Heisenberg Hamiltonians for reproducing the experiments in Section and Appendix!!
Ising
Heisenberg
- j couplings (Jx, Jy, Jz)
(-1.0, 0.0, 0.0)
(1.0, 1.0, 1.0) (1.0, 1.0, 1.0)
h couplligg (hx, hy, hz) (o.0, 0.0, -1.0)
Table 3: Additional non-default parameters for reproducing the experiments in Section 4.1
Command
Values
hypperopt
optim=grid, http 1os
1oss=m11
-eeenee- -arrams
sigma 0=0ma_0=1..0 gaammmmmaaaa gamma-2.0
H Experimental Details
Every numerical experiment, unless stated otherwise, consists of 50 independent seeded trials. Every seed (trial) starts with datapoint, DD (æ., ..)), with æ0 being initial point uniformly drawn from [0, 2m, and o being the associated energy at x." evaluated on the quantum computer. Those 50 initial pairs cached and, when optimization trial starts with given seed, the corresponding cached initial pair is loaded. This allows a fair comparison of different optimization methods: all methods start from the same of initialization points.
The Qiskit 43 open-source library used to classically simulate the quantum computer, whereas the rest of the implementation on Intel Xeon Silver 4316 @ 2.30GHz CPUs, and the code with instructions on how to run and reproduce the results is publicly available GitHub 44].
pure Python. All numerical experiments have been performed
We compare the performance of our proposed VQE-kernel for VQE with the Ising Hamiltonian, i.e., the Heisenberg Hamiltonian (13) with the coupling parameters set to
VQE kernel analysis (Section 4.1
Jx = -1, 0,
hx = 0, by 0, hz=
(45)
and open boundary conditions (see Table 2). For the variational quantum circuit G(æ), we use (Q 3)-qubit, (L = 3)-layered Efficient (U) circuit. this case, the search domain E 00, 2T) according to Eq. (31). The number of readout shots is set to Nahota 1024. The baseline kernels are the Gaussian-RBs kernel,
(- | - x 112
RBF (x. x²) = og exp
(46)
2y²
and the Periodic kernel (67].
D
pperood (x, a²) = =0 exp -
Σ d=1
Id
x/d
sin?
(47)
272
2
which are compared with our VQE kernel (9) in terms of the standard BO performance in Figure 2 Each kernel has two hyperparameters, the prior variance of and the smoothness parameter 7, For three kernels, the prior variance is fixed to oă = 1, and the smoothness parameter is automatically tuned by marginal likelihood maximization (grid search) in each iteration in the early stage (t O, O,..., 75), and after every 100 iterations in the later stage (t = 76....).
For the standard BO, we used the El acquisition function, which maximized by L-BFGS 45 In
the code, the SciPy [68] implementation of L-BFGS was used and all experiments were run using the same default parameter set the code. Detailed commands for reproducing the results can found
in Table 3
28

---

Table 4: Standard choice EMICORe hyperparameters for experiments in Section 4.2 and Appendix (unless specified otherwise).
General params
-n n-qbits
[3,5,7}
# of qubits
n-layers
(3,3,5} esu2 False 3100,,00,5000 {100, 300,500> vqe
of circuit layers Circuit name
-- -circuit
--pbc
Open Boundary Conditions # iterations for BO # iterations for BO
-n-redoout (7x)
-n-iter (Tvi)
kernel
Name the kernel
-hhppropp
Hyperparams optimization
optim
grid 20
Grid-search optimization of
max_gamma
Max value for 7y Scheduling for grid-search # steps in grid Loss type
interval
100*1 +00+++00++00+100
steps
120
loss
m11
--acq-params
EMICoRe params
func
func=ei optim=emicore
Base acq. func. type Optimizer type of candidate points # of evaluation points
optim
pairsize (JsG)
20
gridsize (JOG)
100
corethresh (k)
1.0
CoRe threshold
corethresh_widde (Avv
10
# averaging steps update k
coremin_ scale (Co)
0.0
coefficient Co in eq. (35)
corethresh_ scale (C1)
1.0
coefficient C₁ in eq.
(35
samplesize (NMc)
100
# of MC samples
smo steps (TNFT)
# of initial NFT steps
emo-axis
True
Sequential direction choice
EMICoRe analysis (Section 4.2 In this experiment, compare the optimization performance of our NFT-with-EMICoke with the NFT baselines (sequential and random) on VQE with both the Ising Hamiltonian, for which the coupling parameters are given in Eq. 45 and the Heisenberg Hamiltonian, for which the coupling parameters are set
Jx =y
=
hx = by = hz =1
(48)
For the variational quantum circuit G(z), Efficient SU(2) circuit with open boundary conditions, giving (D= dimeennnsao search domain. The number of readout shots is set to N 1024 in the experiment shown in Figure 3
(Q 5)-qubit, (L 4)-layered
In the same format as Figure 3 Figures 17 in Appendix I.3 compare EMICoRe with the baselines for different setups of Q, L, and Nahots-
For EM (Algorithm 2 and and Algorithm where the VQE kernel is used, the prior standard deviation is set to value roughly proportional to the number of qubits Q; specifically for Q = 5 we set dg 6. The smoothness parameter 7 is automatically tuned by marginal likelihood maximization in each iteration in the early stage (t 0, ..... 100), after every 9 iterations in the middle phase (t = 101,... 280), and after every 100 iterations in the last phase (t = 281, The other parameters
set to TMI 300, JsG 20, JoG 100, NMc 100. TNFT 0 and
Tave = 10. All relevant hyperparameters are collected in Table along with the corresponding flags in our code.
The command options that specify the VQE setting, the kernel optimization schedule, and the other parameter settings for EMICoRe, are summarized in Table 2 and Table 4

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
Exctted
1.0
-5.6
Energy
Fidelity
0.5
-5.8
-6.0
State
0.0
2000
4000
6000
4000 #Observed Points
2000
6000 Trial Density
#Observed Points
Density
Figure 7: Energy (left) and fidelity (right) for EMICoRe (ours) and the baselines, NFT-sequential and NFT-random, up to 6000 observations. Results are for the Ising Hamiltonian with a (L = 33-layered (Q = -uubiis quantum circuit and Nshots 1024.
Table 5: Energy and fidelity achieved after 6000 observations in the experiment in Figure 7 results highlighted in bold, while + (?) indicates lower (higher) values are better.
The best
Algorithm
Energy +
Fidelity
EMICoRe (ours)
5.97 0.05
0.98 + 0.04
NFT-random
-5.92
0.08
0.92 1 0.09
NFT-sequential
5.93 + 0.09
0.92 + 0.16
I
Additional Experiments
Here, we report on additional experimental results.
I.1 Convergence to Ground State
The experiments from Section 4.2 compared the performance between our EMICoRe and the NFT baselines up to 600 observed points, where the optimization has not yet converged. Here, perform a longer optimization with up to 6000 observations to confirm the ability of EMICoRe to converge the ground state. In Figure the energy and fidelity plots show the optimization progress for the Ising model with a (L = 3) 3)-laye (Q -)-bbbsss quantum circuit and Nahots 1024 readout shots. The portrait plot on the right shows the distribution over 50 independent trials of the final solutions after 6000 observations. The mean and the standard deviation of the achieved energy and fidelity are summarized Table We observe that EMICoRe achieves an average fidelity above 95% after 1000 observations, and reaches 98% fidelity at 4000 observations. In contrast, the NFT baselines require all 6000 observations in order to achieve a fidelity of 92%, exhibiting a much slower convergence. This result confirms that EMICoRe robustly converges to the ground state, independent of the individual trial's 's initialization.
Note that the GP regression exhibits cubic computational complexity with respect to the number of samples, thus significantly slowing down the optimization process with thousands of observed points. As a remedy, we limit the number of utilized samples by discarding old observations, i.e., we choose inducing points for the GP based on the chronological order of observations. We found in preliminary experiments that choosing the last 100 observations as inducing points is sufficient achieve good results. In the experiments above (see Figure Z we keep the number of inducing points above 100 and below 120, where we discard the 20 oldest points when exceeding a number of 120 observations. This strategy is implemented in our public code [44] through the option --inducer last_slack: retain= sl0cklack=22 where last_slack indicates the criterion to choose the inducing points, and where retain and slack can used to specify the minimum number of points retained and the number of samples that can be observed beyond the minimum, before discarding. Hence, the model discards slack=20 observations when the number of observed points equals to the sum of the two, i.c., slack retained.
30

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
VQE-EI
-5
1.0
Ground
Energy
Fidelity
0.5
-7
Ground State
0.0
100
200
300 Trial Density
100
200
300 Trial Density
# Observed Points
#Observed Points
(a) Ising off-criticality
1st Excled
Ground State
-5.6
Energy
Fidelity
0.5
-5.8
-6.0
Ground State
0.0
200
400
600 Trial
0
400 #Observed Points
200
600 Trial
#Observed Points
Density
Density
(b) Ising at criticality
Figure 8: Energy (left) and fidelity (right) for EMICoRe (ours) compared three different baselines: NFT--sequentia,, NFT-random, and El with VQE kernel. Results for the Ising Hamiltonian off- criticality and criticality are shown in the top and bottom rows respectively.
I.2 Ising Hamiltonian Off Off-Criticality
In Section 4 we focused on the Ising and Heisenberg Hamiltonians with the parameters, Jy, Jz) and h = (hx, by, hz) ) in Equation (3). set to criticality in the thermodynamic limit. Such choices are expected to be most challenging for the VQE optimization because the corresponding ground states tend to be highly entangled due to the quantum phase transition. As an ablation study,
we here conduct experiments for an off-critical setting. Specifically, we evaluate the optimization performance for the Ising Hamiltonian off-criticality [J= (0, 0, -1); h = (1.5, O, 0) Figure R (top) shows the energy (left) and the fidelity (right) achieved by EMICoRe (ours), NFT-sequential. NFT-random, and EI with VQE kernel after 600 iterations for (L 3y-layered (Q = 5) 5)-qubits
Figure quantum circuit with Nahots 1024 readout shots. For comparison, we also show (bottom) the performance for the Ising Hamiltonian at criticality (J = (-1, 0, 0); / (0, 0, -1)]. observe that the off-criticality setting significantly easier, while the plain BO with El (without EMICoRe) falls somewhat short, although closely behind NFT-random at 600 observed points.
I.3 Different Setups for Qubits and Layers
Here, we compare the performance of the baselines (NFT-sequential and NFT-random) our WrTT-with-EMCCoee under different settings of Q. L. and Nahots All figures in this appendix for each figure, the energy (left column)
(Figures g I7 are shown in the same format as Figure and the fidelity (right column) are shown for the Ising (top row) and the Heisenberg (bottom row) Hamiltonians. In each panel. the left plot shows the optimization progress with the median (solid) and the 25- and 75-th percentiles (shadow) over the 50 seeded trials, as described in Appendix H as a function of the observation costs, i.e., the number of observed points. The portrait square on the right shows the distribution of the final solution after 200, 600, and 1000 observations have been performed, respectively, for Q 3, 5, and 7 qubit cases. As mentioned earlier, the prior standard deviation do is set roughly proportional to Q. Specifically we use do 4, 6, 9 for Q 3, 5, 7, respectively.
3-qubit setup: Figures the relatively low-dimensional problem with only = 24, the convergence rate comparable for both baselines and approach. However, the red sharp peak in the density plot of the final solutions
show results for Q = T. 3, and Nahots 256, 512, 1024. Given
31

---

(right portrait square) in each panel implies that the robustness against initialization is improved by our WiFT-with-MMICoR approach, thus highlighting its enhanced stability and voise-resiliency.
5-qubit setup: Figures 12 14 show results for Q = 5, L = 3, and Nshots 256, 512, 1024.
The case for Nohots 1024 identical to Figure 3 in the main text. For all noise levels (Natuots 256, 512, 1024), our EMICoRe (red) consistently achieves lower energy and higher fidelity compared to the baselines, thus demonstrating the superiority of our NFT-with-EMMICRR over NFT [ Remarkably, we observe that in high-noise-level cases, such as for the Heisenberg Hamiltonian for Nshots = 256, 512 (bottom rows of Figure 12 and Figure (3), the achieved energy by the state-of-the- art baselines (purple and green) fail to surpass the energy level of the first excited state for the 600 observed data points, whereas EMICoRe successfully accomplishes this task.
7-qubit setup: Figures 15-17 present results for Q 7, L 5, and Nahots 256, 512, 1024. Again, NFT-wwth-EMICCOM consistently outperforms the baselines all experimental setups. Given the increased complexity associated with Q and D 84, the optimization process becomes more challenging, necessitating greater number of observed points to approach the ground energy. Nonetheless, even with just 1000 observed points. NFT- wirt-with-aa EMICoEE already exhibits significant superiority over NFT [11), particularly in high-noise scenarios such as Nghots 256. We also observed that the optimization process faces difficulties with the Heisenberg Hamiltonian We attribute this behavior to the greater complexity the latter task and defer further analysis of this regime to future studies.
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
1st Excred State
Ground State
-2.5
Energy
Fidelity
-3.00
0.5
33.25
3.50
Ground State
0.0
100
150
200 Trial Density
100
150
200 Trial Density
#Observed Points
#Observed Points
1.0
Ground State
Energy
Fidelity
0.5
Ground State
0.0
50
100
150
200 Trial
50
100
150
200 Trial
#Observed Points
Density
#Observed Points
Density
Figure 9: Comparison (in the same format as Figure 3) between our NFT-withhEMMCoRe (red) and the NFT baselines (green and purple) in VQE for the Ising (top row) and Heisenberg (bottom row) Hamiltonians with the (L = 3yyleyeeed (Q = 3)-qubit quantum circuit (thus, D = 24) and
Nohots 256.
32

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-275
Energy
Fidelity
-3.00
0.5
-3.25
-3.50
0.0
50
100
150
200
Trial Density
0
50
100
150
200 Trial Density
#Observed Points
Observed Points
1sf Lytted State
1.0
Ground State
Energy
-4
Fidelity
0.5
-5-
0.0
150
200 Trial
100 150
Trial Density
# #Observed Points
Density
# #Observed Points
31. 3a-yeeee (Q 3)- quantum circuit (thus, D = 24) and Nahots = 512. The NFT-with-EMMICol (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
Figure 10: Same comparison in Fig. with the (L. =
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-2.75
Energy
Fidelity
33,00
0.5
-.25
-3.50
Grann
0.0
100
150
200 Trial Density
50
100
150
200
Trial Density
#Observed Points
Observed Points
1x Facited State
1.0
Ground State
Energy
-4
Fideltty
0.5
-5
Ground
0.0
150
200 Trial
100 150
Trial Density
# #Observed Points
Density
O'bserved Points
Figure 11: Same comparison in Fig. with the (L= 3) layered (Q 3)-qubit quantum circuit (red) and NFT baselines (green and
(thus, D 24) and Nehots = 1024. The purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
33

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-5.0
Energy
Fidelity
0.5
-5.5
Teeeeee 1st Exclted State
-.0
0.0
200
400
600 Trial Density
0
200
400
600 Trial Density
#Observed Points
Observed Points
1.0
Ground State
Energy
-8-
Fidelity
0.5
1st Excited State
-9
0.0
400 # #Observed Points
200
600 Trial
200
400
Trial Density
Density
#Observed Points
3 3ay-eyed (Q 5)-qubit quantum circuit. (thus, D = 40) and Nahots = 256. The NFT-with-EMMICon (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
Figure 12: Same comparison in Fig. with the (L. =
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
Energy
Fidelity
0.5
-5.5
1st Excited
-6.0
0.0
200
400
600 Trial Density
200
400
600 Trial Density
#Observed Points
Observed Points
1.0
Ground State
Energy
Fidelity
-8-
0.5
10 Excited State
-9
Ground
0.0
400 # Observedd Points
200
600 Trial
400 #Observed Points
Trial Density
Density
= 3) llyeerd (Q 5)-qubit quantum circuit (thus, D = 40) and Nahots = 512. The NFT--with-EMMCoR (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
Figure 13: Same comparison in Fig. with the (L=
34

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-5.0
Energy
Fidelity
0.5
-5.5
Treeeee 1st ExcReo S2302
-6.0
0.0
200
400
600 Trial Density
0
200
400
600 Trial Density
#Observed Points
Observed Points
1.0
Ground State
Energy
Fidelity
-8-
0.5
1st Excited State
-9
0.0
400 # Obserreed Points
200
600 Trial
200
400
Trial Density
Density
#Observed Points
Figure 14: Same comparison in Fig. with the -L= = 3 -ayyeeed (Q 5)-qubit quantum circuit.
(thus, D = 40) and Nshots = 1024. The NFT-with-EMMCCo (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-7
Energy
Fidelity
0.5
1st Facited State
w w
0.0
250
500
750
1000 Trial Density
250
500
750
1000 Trial Density
#Observed Points
# Observed Points
-6
1.0
Ground State
Energy
Fidelity
0.5
-10
-12
Ground
0.0
250
750
1000 Trial
500
750
1000 Trial
#Observed Points
Density
#Observed Points
Density
Figure 15: Same comparison in Fig. with the (L (L= llayeerd (Q 7)-qubit quantum circuit
(thus, D = 84) and Nahots = 256. The NFT--withEMMICol (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
35

---

EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-7
Energy
Fidelity
0.5
1at Excited State 2
0.0
250
500
750
1000 Trial Density
0
250
500
750
1ria
#Observed Points
# Observed Points
Density
1.0
Ground State
Energy
Fidelity
0.5
-10-
-12
Broind
0.0
250
750
1000 Trial
500
750
1000 Trial
#Observed Points
Density
#Observed Points
Density
Figure 16: Same comparison in Fig. with the (L. = 5a-leyed (Q 7)-qubit quantum circuit.
(thus, D = 84) and Nahots = 512. The NFT-withhEMMICol (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
EMICoRe (ours)
NFT-Sequential
NFT-Random
1.0
-7
Energy
Fidelity
0.5
1st Facited State
w
0.0
250
500
750
1000 Trial Density
250
500
750
1000 Trial Density
#Observed Points
# Observed Points
-6
1.0
Ground State
Energy
Fidelity
0.5
-10
-12
Ground
0.0
250
750
1000 Trial
500
750
1000 Trial
#Observed Points
Density
#Observed Points
Density
Figure 17: Same comparison in Fig. with the (L= -ay-eed (Q 7)-qubit quantum circuit
(thus, D 84) and Nehots = 1024. The WFT-with-EMMCoRE (red) and NFT baselines (green and purple) shown for both Ising (top row) and Heisenberg (bottom row) Hamiltonians.
36


Reviews:
[{"review_id": "tLGOfThaUj", "reviewer": "Reviewer_93Tb", "summary": "I have reviewed the rebuttal, and I intend to maintain my decision. I appreciate the authors made the efforts to clarify most of my concerns. However, I believe it is essential to address the issue of hardware noise in the present study, a point that previous research seems to have overlooked. I strongly encourage the authors to delve deeper into this matter, as a robust algorithm should be capable of practical implementation on real quantum computers. I believe the issue of hardware noise is actually an inescapable and important challenge here.\n\nThis manuscript presents an innovative approach known as EMICoRe, which leverages the synergy between Bayesian optimization and prior knowledge of variational quantum eigensolvers (VQE) to enhance the efficiency of the optimization. The authors introduce a novel kernel specifically designed for Bayesian optimization and employ it within the EMICoRe framework. Additionally, the paper conducts numerical experiments to compare the proposed method with the state-of-the-art baselines, demonstrating its superior performance.", "questions": "1. Based on the numerical results presented in the paper, the fidelity between the ground state and the state obtained through the VQE remains consistently below 0.8, even after reaching 400 observed points. It is noteworthy that the VQE circuit being optimized in the study comprises only three layers and five qubits. Therefore, I believe it is crucial to investigate the potential impact of further increasing the number of observed points and determine the minimum number of observed points required to achieve a higher fidelity, such as surpassing 0.95.\n\n2. The manuscript just includes two parameter settings, one for the Ising model and the other for the Heisenberg model. However, to verify the robustness of the proposed method, I suggest that the authors incorporate more scenarios for validation.", "limitations": "I believe the authors have provided thorough discussions on the limitations of their work within the Appendix.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The manuscript introduces a novel optimization algorithm based on Bayesian optimization that incorporates the specific property of the VQE objective function into its design.\n- The paper is well-organized and clearly written. The figures in the paper effectively illustrate the improved performance achieved by the proposed method in comparison to the baselines.\n- The paper establishes the equivalence between two previously proposed VQE properties, providing a foundational basis for the subsequent analysis on VQE objective functions.", "weaknesses": "- The numerical results presented in the paper indicate that the proposed method does not exhibit a clear advantage when the number of observed points is limited, as compared to the baselines.\n- Although the proposed method exhibits improved performance compared to the baseline as the number of observed points increases, it still falls short of achieving the ideal ground state.\n- The proposed method just considers the shot noise but ignores possible noise originating from the quantum circuit itself.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback and appreciation of our manuscript.\n\n## **Weaknesses**\n> The numerical results presented in the paper indicate that the proposed method does not exhibit a clear advantage when the number of observed points is limited, as compared to the baselines.\n\nIndeed, when the number of observed points is small the GP surrogate model is not informative enough to provide an effective acquisition function. When the number of observed points gets larger, and thus the GP surrogate model becomes more informative, our EMICoRe approach starts outperforming standard baselines, i.e., our informed search for new points to observe, based on the CoRe, shows its benefit. \n\n> Although the proposed method exhibits improved performance compared to the baseline as the number of observed points increases, it still falls short of achieving the ideal ground state.\n\nWe would like to point out that the primary objective of our work is to demonstrate that EMICoRe can improve the efficiency of the training process, and not to extract the physics of the Ising model, which is well known in the literature. Thus, we compare the performance of the training process for different methods under a fixed budget of observations. Our method reaches the ground state by simply performing longer optimizations. To demonstrate this, we ran additional experiments with a larger number of observations, see Fig. R2 in the PDF. From this plot, it is evident that taking more observations will make our EMICoRe converge to the ground state as expected.\n\n> The proposed method just considers the shot noise but ignores possible noise originating from the quantum circuit itself.\n\nWhen benchmarking the performance of newly proposed hybrid quantum-classical algorithms, the noise from quantum hardware is usually not included in first-case studies (see e.g., [1,2,3]). A standard procedure in the field is to benchmark first on shot noise only, see the NFT [1] experiments, and consider hardware noise in potential follow-up studies. As stated in the general rebuttal, the main reason is that quantum noise is strongly hardware dependent, e.g., superconducting quantum hardware is affected by fundamentally different types of noise (in particular, CNOT gate noise, decoherence, and measurement noise) compared to, e.g., trapped-ion quantum hardware. Moreover, various error mitigation schemes exist for different types of quantum hardware, which generally require additional error calibration runs on the hardware. An informative study of quantum hardware noise would need to take all these considerations into account and therefore is beyond the scope of the initial proposal of our novel method. \n\n## **Questions**\n> Based on the numerical results presented in the paper, the fidelity between the ground state and the state obtained through the VQE remains consistently below 0.8, even after reaching 400 observed points. It is noteworthy that the VQE circuit being optimized in the study comprises only three layers and five qubits. Therefore, I believe it is crucial to investigate the potential impact of further increasing the number of observed points and determine the minimum number of observed points required to achieve a higher fidelity, such as surpassing 0.95.\n\nWe agree with the referee and we have performed (longer) experiments, see Fig. R2, in the attached PDF. In this analysis, we show that running the optimization longer, i.e., performing more observations, lets the energy converge to the ground state and the fidelity to approach $0.976$.  Looking at the figure, we also note that EMICoRe requires 1000 observations in order to push the average fidelity above $0.95$, while NFT never reaches that value. In Table R2 below, we report the mean and standard deviation for the energy and the fidelity for the EMICoRe and NFT baselines when running for 6000 observations, associated to the experiments mentioned above.\n\n> The manuscript just includes two parameter settings, one for the Ising model and the other for the Heisenberg model. However, to verify the robustness of the proposed method, I suggest that the authors incorporate more scenarios for validation.\n\nOur experiments focus on the Heisenberg model with the Ising model as a special subclass of the Heisenberg one. In both cases, the parameters are chosen such that in the thermodynamic limit, they correspond to critical points, i.e., where a quantum phase transition occurs. This choice of setup is in general among the most challenging scenarios, as the corresponding ground states show the largest amount of entanglement. We thus expect any other choice of parameters to render the problem easier. To complement our results, we added another choice of couplings in the PDF document, see Fig. R1. We plan to further extend our study to other Hamiltonians, e.g., quantum chemistry systems, in future work.\n\n\n### Table R2: Energy and fidelity for longer runs (6000 observations) for the Ising model at criticality. Best results highlighted in bold. For energy, lower is better. For the fidelity, higher is better.\n| kernel             |  Energy  | Fidelity |                      \n|--------------------|----------|----------|\n| **EMICoRe (ours)** | **-5.973762 +- 0.045612**   | **0.975552 +- 0.039288** |\n| NFT-random | -5.929573 +- 0.075909    |  0.918256 +- 0.097190|                         \n| NFT-sequential |  -5.930662 +- 0.097275    |  0.916873 +- 0.155590 | \n\n## **References**\n\n\n- [1] [Nakanishi K. et al., Phys. Rev. Research 2, 043158 (2020).](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043158)\n- [2] [Farhi E. et al., arXiv:1411.4028 (2014)](https://arxiv.org/abs/1411.4028)\n- [3] [Bravo-Prieto C. et al., Quantum 4, 272 (2020).](https://quantum-journal.org/papers/q-2020-05-28-272/) "}, {"review_id": "tQSUI1LZko", "reviewer": "Reviewer_WANE", "summary": "The authors propose a method for Bayesian Optimization for Variational Quantum Eigensolvers, which they call NFT with EMICoRe.  This method uses a novel VQE kernel, which constrains the function space of the Gaussian Process underlying the BO to include only valid VQE objective functions (using the representation in Prop. 2, derived from NFT).  The authors also propose a novel acquisition function for the EMICoRe method (Eq. 11), which optimizes over the expected maximum improvement over confident regions.  In their experimentation, the authors show that their VQE kernel is able to outperform other kernels in a BO setting, and that their NFT-EMICoRe approach is able to outperform other (non-BO) NFT approaches.", "questions": "The threshold parameter \\kappa is introduced in Sec. 3.2, but doesn't appear to be explored in the experimentation.  In particular, did the authors verify that an intermediate value of this parameter is advantageous (hence supporting the use of the EMICoRe acquisition function)?  I may have missed this in the experimentation.\n\nFurther, on lines 133-134, they state that they are not concerned with circuit noise on current NISQ devices, but line 320 states that the experimentation confirms the suitability of their method for such devices - is such noise explicitly included in the simulations?", "limitations": "Yes, limitations are addressed.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper tackles an important problem, the optimization of noisy VQE circuits, and offers a principled solution using BO combined with physical constraints, using state of the art methods (NFT).  The paper is well written, and the experimentation is well chosen to support the method.", "weaknesses": "The experimentation could be expanded.  Particularly, it would be interesting to see how the model performs on an actual quantum implementation.  Also, investigation of a broader ranger of Hamiltonians would be desirable (including ones motivated by practical problems).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for their insightful comments and positive feedback on our work. We take this opportunity to address and respond to the comments below:\n\n## **Weaknesses**\n> The experimentation could be expanded. Particularly, it would be interesting to see how the model performs on an actual quantum implementation. Also, investigation of a broader range of Hamiltonians would be desirable (including ones motivated by practical problems).\n\nWe agree that our method is suited for a broader range of applications as well as Hamiltonians. The Ising and Heisenberg Hamiltonians, however, represent standard benchmarks that are of high practical relevance. They correspond to spin chain Hamiltonians that are widely studied in condensed matter physics, see e.g. Refs. [1,2]. Furthermore, many lattice field theories can be represented as generalized spin chains (see, e.g., Eq. (4)-(6) in [3]). As far as the implementation on actual quantum devices is concerned, the standard benchmarking is usually performed on simulated, noiseless devices which only account for shot noise. We plan to investigate hardware noise in future work. For more details please see the fourth bullet point in the general rebuttal.\n\n## **Questions**\n> The threshold parameter $\\kappa$ is introduced in Sec. 3.2, but doesn't appear to be explored in the experimentation. In particular, did the authors verify that an intermediate value of this parameter is advantageous (hence supporting the use of the EMICoRe acquisition function)? I may have missed this in the experimentation.\n\nThank you for raising this crucial question. We performed an additional investigation for different heuristic for setting the parameter $\\kappa$. This study can be found in the general rebuttal, see Table R1, and will be added to the appendix of the paper in its updated version.\n\n> Further, on lines 133-134, they state that they are not concerned with circuit noise on current NISQ devices, but line 320 states that the experimentation confirms the suitability of their method for such devices - is such noise explicitly included in the simulations? \n\nIn our experiments, we only consider shot noise throughout the paper. We apologize for any confusion our phrasing may have caused. We will revise the main text to make this clearer.\n\n## **References**\n\n- [1] [Funcke L. et al., arXiv preprint arXiv:2302.00467 (2023)](https://arxiv.org/pdf/2302.00467.pdf)\n- [2] [Di Meglio A. et al, arXiv preprint arXiv:2307.03236 (2023).](https://arxiv.org/pdf/2307.03236.pdf)\n- [3] [Atas Y. et al, Nature communications 12, no. 1:6499 (2021).](https://arxiv.org/pdf/2102.08920.pdf)"}, {"review_id": "Qzzb525Q5D", "reviewer": "Reviewer_1i6G", "summary": "In this work, the authors integrate a quantum kernel method with the EMICoRe architecture to further improve the NFT framework of Bayesian Optimization. The simulation results show the advantages of the proposed approach. ", "questions": "(1) What are the quantum advantages of the quantum kernel method? Since the quantum kernel method is combined with another proposed method `EMICoRe', we do not know which part contributes to the improvement of the Bayesian Optimization method. --- resolved", "limitations": "It would be better to reformulate the quantum approach in the authors' proposed Bayesian Optimization framework. ", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "(1) The method of leveraging the quantum method for Bayesian Optimization is interesting. \n\n(2) The investigation of how to incorporate the quantum kernel method is significant. \n\n", "weaknesses": "(1) Since no parametric circuits are implemented in the quantum kernel model shown in Eq. (9), the proposed VQE-kernel is nothing but a quantum kernel learning method, which has been comprehensively studied in previous work in Refs. [1-3]. In particular, Ref. [4] has exhibited the use of quantum kernel learning for improving the performance of Bayesian optimization. --- resolved\n\n(2) The circuit diagram for the quantum kernel learning is not provided such that the experiments cannot be easily reproduced. --- resolved\n\n(3) There are two main contributions to this work: the use of quantum kernel learning and an introduction to EMICoRe. Although the simulation results demonstrate the performance improvement, it is still unknown where the performance gains come from, and the quantum advantages of the quantum kernel are not analyzed at all. --- resolved\n\n[1] Havlíček, Vojtěch, et al. \"Supervised learning with quantum-enhanced feature spaces.\" Nature 567.7747 (2019): 209-212\n\n[2] Wang, Xinbiao, et al. \"Towards understanding the power of quantum kernels in the NISQ era.\" Quantum 5 (2021): 531\n\n[3] Blank, Carsten, et al. \"Quantum classifier with tailored quantum kernel.\" npj Quantum Information 6.1 (2020): 41\n\n[4] Rath, Yannic, Aldo Glielmo, and George H. Booth. \"A Bayesian inference framework for compression and prediction of quantum states.\" The Journal of Chemical Physics 153.12 (2020)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the feedback. However, it seems as if the referee misunderstood the paper, as our work cannot be characterized as quantum kernel learning. In essence, quantum kernel leverages a quantum circuit to calculate a certain kernel, i.e., the calculation of the kernel is non-classical. In our work, we use classical kernels. Specifically, we propose a novel *classical* kernel which is specifically tailored to the VQE setup. More precisely, the corresponding feature vectors precisely correspond to the energy landscape modeled by the VQE. This is an important and novel contribution in its own right. We go beyond this insight and deploy this kernel in Bayesian optimization. Specifically, we propose EMICoRe and demonstrate that it leads to superior performance when compared to the current state-of-the-art method for optimizing VQE, i.e. NFT. In the following, we reply in more detail to the reviewer's comments:\n\n## **Weaknesses**\n\n> Since no parametric circuits are implemented in the quantum kernel model shown in Eq. (9), the proposed VQE-kernel is nothing but a quantum kernel learning method, which has been comprehensively studied in previous work in Refs. [1-3]. In particular, Ref. [4] has exhibited the use of quantum kernel learning for improving the performance of Bayesian optimization.\n\nAs discussed above, our work cannot be characterized as quantum kernel learning, which is the topic of the suggested references. We have added a discussion of the notable differences with respect to quantum kernel learning to the related work section in the updated version of the manuscript.\n\n> The circuit diagram for the quantum kernel learning is not provided such that the experiments cannot be easily reproduced.\n\nAs we do not use a quantum kernel, we cannot provide a circuit diagram for it. In more detail, for quantum kernel learning, one would need two circuits, see Figure 1 of Ref. [2] that the referee has mentioned. In our case, we have a classical kernel. Therefore, we only require one circuit: the Efficient SU(2) circuit, as shown in Figure 6 in the Appendix of our paper for reproducibility. We also would like to emphasize that we provided the code for reproducing all our experiments, including the implementation of the quantum circuit. \n\n> There are two main contributions to this work: the use of quantum kernel learning and an introduction to EMICoRe. Although the simulation results demonstrate the performance improvement, it is still unknown where the performance gains come from, and the quantum advantages of the quantum kernel are not analyzed at all.\n\nAs discussed above, our work does not study quantum kernels, thus, no possible quantum advantage from a quantum kernel can be analyzed. Regarding the quantum advantage of VQE and the performance gain, please refer to the following answer. \n\n## **Questions**\n> What are the quantum advantages of the quantum kernel method? Since the quantum kernel method is combined with another proposed method `EMICoRe`, we do not know which part contributes to the improvement of the Bayesian Optimization method.\n\nAs we pointed out in the **Weakness** section, our work does not consider quantum kernels, so no quantum advantage of a quantum kernel can be studied. Instead, our work proposes a *classical* kernel that is uniquely suited for the VQE setup. Specifically, the corresponding feature vectors precisely align with the energy landscape modeled by the VQE. In particular, it retains the property of VQEs that the energy can be determined along an entire line by only fixing three points. This powerful inductive bias allows us to propose EMICoRe, which can outperform the current state-of-the-art method [1]. For the VQE algorithm, it is very well known that the quantum advantage comes from the efficient implementation of the exponentially large Hilbert space using the qubits of a quantum computer. A quantum computer only needs $Q$ qubits to prepare the quantum state $|\\psi(\\theta)\\rangle$ in the VQE algorithm, while a classical computer would require exponentially large resources, $2^Q$, for the same operation.\n\n## **References**\n- [1] [Nakanishi K. et al., Phys. Rev. Research 2, 043158 (2020).](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043158)\n- [2] [Wang, Xinbiao, et al. \"Towards understanding the power of quantum kernels in the NISQ era.\" Quantum 5 (2021): 531](https://quantum-journal.org/papers/q-2021-08-30-531/)"}, {"review_id": "MUuSKMvp0p", "reviewer": "Reviewer_nt6b", "summary": "The paper introduces a new approach for BO of VQEs. BO can be a good match for this problem since it models the noise (measurement and circuit level) of quantum circuits. The main idea of the paper is to use a kernel that is adapted to the form of the VQE objective assumed when variational parameters are associated to single qubit gates. The paper also introduces EMICoRE and NFT-with-EMICoRE, acquisition function and parameter update strategies. Overall, the idea of using an adapted kernel for this task is well motivated, the paper is well written and presents numerical evaluation of their method. I think however that in the current state, the experimental evaluation is too limited to be convincing of the promise of this method with respect to other optimisation algorithms for VQEs.", "questions": "1. Could you benchmark your method against other methods for VQE optimisation than NFT? And for larger systems?\n2. Why is the ground state energy not reached in figure 3 say? I would expect that it is possible to find a VQE ansatz that can solve this easy problem.\n3. Can you add ablations for the EMICoRE acquisition function? Have you tried different heuristics for updating $\\kappa$?\n4. Can you explain how your method can be used to optimise general parametrised circuits beyond those that have single qubit gates?\n", "limitations": "- limited experiments and benchmarking\n- limited explanation of their choice of parametrised circuits", "rating": 6, "confidence": 4, "soundness": 1, "presentation": 3, "contribution": 2, "strengths": "The paper addresses an important problem in quantum computing. It is well written and accessible for both quantum computing and ML communities I think. It introduces the following technical innovations:\n- Kernel adapted to VQEs: the kernel decomposition based on the form of the cost function of VQEs is a sound contribution.\n- NFT-with-EMICoRE: a new BO algorithm that builds on NFT and a new acquisition function.\n\n", "weaknesses": "1. Specific ansatz: the paper deals only with parametrised single qubit gates. While I believe this is not a theoretical limitation - as one can always decompose any unitary as a product of single qubit gates and CNOTs - it is not clear to me how practical is this parametrisation. The authors could add an explanation of whether this is a limitation for practitioner and what steps one needs to do to apply their algorithm to a generic parametrised quantum circuit\n2. Limited experimental evaluation: the paper only discusses Q=3 and 5 qubit systems, while classically one could simulate easily bigger systems. Also the paper does not compare against other methods in the literature, it only compare against NFT. More benchmarking for larger systems and against other techniques would be required to assess the promise of the method. \n3. Unclear benefit of the kernel: looking at figure 2, it is unclear that the red curve is better than blue curve - in fact on the left figure it looks to me that RBF is doing better.\n4. Missing ablation for the EMICoRE acquisition function. One contribution of the paper is to introduce this new acquisition function, but ablations showing its importance are missing.\n\nMinor:\n- $\\mu_X$ depends on y as well in section 2.1\n- \"the the\" in line 308\n- unclear what conjugacy means in line 97", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for highlighting the strengths of our work. We hereby address the concerns raised in the review:\n\n## **Weaknesses**\n> 1. Specific ansatz: ...\n\nThe set of single-qubit rotation gates and entangling CNOT gates are indeed universal, and allow for synthesizing all unitary operations on a subset of qubits in an efficient manner, see [1]. More generally, this is the essence of the Solovay-Kitaev theorem [2]. Moreover, most commercially available superconducting quantum devices offer single-qubit rotation gates and CNOT entangling gates, such that our ansatz does not pose a restriction in practice. In fact, Qiskit, one of the most prevalent SDKs for developing quantum programs, offers the ansatz structure we use as *\"a heuristic pattern that can be used to prepare trial wave functions for variational quantum algorithms or classification circuit for machine learning\"* (see  [qiskit.org/documentation](https://qiskit.org/documentation/stubs/qiskit.circuit.library.EfficientSU2.html)). Furthermore, our work is not limited to the type of ansätze we studied in the manuscript. The only restriction is that the parametric gates in the ansatz are of the form $R(\\theta) = \\exp(-i\\theta G)$, where $G$ is an arbitrary gate operation that fulfills $G^2 = \\mathbf{1}$, see Eq. (2) in the NFT paper [3]. This comprises any gate where $G$ is a tensor product of Pauli operators, for example parametric entangling gates of the form $\\exp(-i\\theta Z_l\\otimes Z_k)$ with $Z_j$ the Pauli $Z$-matrix acting on qubit $j$.\n\n> 2. Limited experimental evaluation: ...\n\nIn the paper's Appendix G, we show results for Q=7 in Figs. 13-15. The method we compare to, NFT [3], is widely accepted as the current state-of-the-art. Specifically, in [3], the authors showed superior performance compared to many gradient-based and gradient-free methods, i.e., SPSA[4], BFGS[5], N-M[6], Powell[7], CG[8]. We, therefore, refrained from including these inferior baselines to streamline the presentation. Our numerical experiments clearly establish that EMICoRe gains favorable scaling with the system size, i.e., our method can beat the state-of-the-art method by an increasing margin as the number of qubits is increased, see Figs. 7-14 in Appendix G.\n\n> 3. Unclear benefit of the kernel: ...\n\nIn Table R3 below, we show that our VQE kernel reaches, on average, lower energy and higher fidelity, with smaller variance, when compared to RBF and periodic kernels. Furthermore, the VQE kernel is a necessary ingredient in order to perform the EMICoRe strategy which is proven to outperform the state-of-the-art baseline in our experiments. \n\n> 4. Missing ablation ...\n\nIn Fig. 2 in the main paper, we evaluated the performance gain of our VQE kernel without EMICoRe.\nTo directly show the ablation we conducted additional experiments with $Q=5$ and $L=3$ in the rebuttal PDF. We showcase the benefits of EMICoRe by comparing it to the widely adopted Expected Improvement (EI) acquisition function for two types of systems, i.e., the Ising model at and away from criticality.\n\n## **Questions**\n> 1. Benchmark, larger systems\n\nIt is well-known that NFT is the current state-of-the-art algorthm for optimizing VQEs. We therefore refrained from including other inferior baselines.\nFurthermore, as mentioned above, our method can beat the state-of-the-art method by an increasing margin as the system size increases. We refer to our reply in the **Weakness** section above.\n\n> 2. Ground state energy in figure 3\n\nIn Figure 3 from the main paper, our main goal was to demonstrate the faster optimization of EMICoRE compared to NFT at a fixed quantum measurement budget. In Fig. R2 in the attached PDF, we demonstrate that our algorithm can reach the ground state for a larger measurement budget. We will include these experiments in the revised manuscript.\n\n> 3. Ablations\n\nWe performed additional experiments in the direction of an ablation study for EMICoRe. This includes comparing to other acquisition functions, e.g., EI, as shown in Fig. R1 in the rebuttal PDF, and performing an in-depth analysis for setting the CoRe threshold $\\kappa$, shown in the Table R1 from the general rebuttal above.\n\n> 4. General parametrised circuits\n\nAs we pointed out in the **Weakness** section, our method is not limited to the parametric gates we use in our ansatz nor is it limited to single-qubit parametric gates. In general, our approach can be applied to any parametric gate of the form $R(\\theta) = \\exp(-i\\theta G)$ with $G$ a gate that has the property $G^2=\\mathbf{1}$, see Eq. (2) in the NFT paper [3]. This includes a wide range of gates, in particular, $G$ can be an arbitrary Pauli string $\\otimes_{q=1}^QP_q$, where $P_q$ is one of the Pauli matrices {$\\mathbf{1}, X, Y, Z$} acting on qubit $q$. More specifically, this set also includes the entangling gates $R_{XX}=\\exp(-i\\theta X_l\\otimes X_k )$ and $R_{ZZ}=\\exp(-i\\theta Z_l\\otimes Z_k )$, which are more commonly realized in trapped ion setups.\n\n\n### Table R3: Energy and Fidelity for 600 observed points and different kernels, as shown in Fig. 2. Best results highlighted in bold. For energy lower is better. For fidelity, higher is better.\n| Kernel | Energy | Fidelity |\n|-|-|-|\n| **vqe (ours)** | **-3.306890 ± 0.117974** | **0.926119 ± 0.047960** |\n| rbf | -3.269373 ± 0.145235 |  0.895492 ± 0.084795 |\n| periodic | -3.232246 ± 0.282348 | 0.901823 ± 0.095861 |\n\n## **References**\n- [1] A. Barenco A. et al., Phys. Rev. A 52, 3457 (1995).\n- [2] A. Y. Kitaev,  Russ. Math. Surv. 52 _1191_ (1997).\n- [3] K. Nakanishi et al., Phys. Rev. Research 2, 043158 (2020).\n- [4] J. C. Spall, IEEE Transactions on Automatic Control, 37(3):332–341, 1992.\n- [5] D. C. Liu, J. Nocedal, Mathematical Programming volume 45, pages 503–528 (1989).\n- [6] J. A. Nelder and R. Mead, The Computer Journal 7, 308 (1965).\n- [7] M. J. D. Powell, The Computer Journal, Volume 7, Issue 2, 1964, Pages 155–162 (1964).\n- [8] R. Fletcher and C. M. Reeves, Computer Journal 7 (1964)."}]

Author rebuttal:
We thank the four reviewers for their valuable feedback. To streamline our reply, we place the referenced tables at the bottom. 

- Some reviewers suggested additional experiments for other parameter choices of the target Hamiltonians. We stress that the Ising Hamiltonian, which is studied along with the Heisenberg Hamiltonian, is considered at criticality (in our paper) and thus represents a challenging optimization objective. Studying the Hamiltonians for coupling values off criticality can be expected to be less challenging. To confirm this, we have performed additional experiments for different coupling parameters. These are shown in Fig. R1 in the attached PDF.

- Furthermore, some reviewers suggested to explore further possibilities of heuristics to set the CoRe threshold $\kappa$. As a result, we investigated this as part of additional experiments. Specifically, we used the following updating rule for tuning the value of $\kappa$ at each step:
$$
\kappa = \max\left(C_0\cdot\sigma, C_1\cdot\frac{ \hat{\mu}^{t - T_{\mathrm{Ave}}} - \hat{\mu}^{t}}{T_{\mathrm{Ave}}}\right)
$$with hyperparameters $C_0$ and $C_1$. We note that $\sigma$ is the standard deviation of the observation noise, and fixing $C_0=0, C_1=1.0$ reduces to the heuristics used in the original submission. 
We report the results for different choices of the hyperparameters $C_0$ and $C_1$ in Table R1 below. The best combination is highlighted in bold. The experiments were performed for the Ising model at criticality using EMICoRe and a $(L = 3)$-layered $(Q = 5)$-qubits quantum circuit with $N_\textrm{{shots}}=1024$. Each trial ran for 600 observations. Table R1 demonstrates that intermediate values for the hyperparameters $C_0, C_1$ achieve the best result, thus proving the effectiveness of EMICoRe. Also, we found that the combination $C_0=0.1, C_1=10$ improves the result from the original submission. We thank the reviewers for the suggestion, and we will update the results in the paper with the new heuristics. 

- Some reviewers also asked why the ground state is never reached by our method as well as any baseline. This is merely a consequence of taking a limited number of observations, and we performed longer optimizations to confirm this. Specifically, Fig. R2 in the attached PDF demonstrates that running the experiments longer, i.e., observing more points, makes the optimization converge closer to the ground state and raises the fidelity to  $0.976$, see table R2 in the reply to referee 93Tb for details. We explicitly tested that the same holds true for every other choice of parameters and circuit setups investigated in this manuscript, e.g., $(Q,L)=(3,3), (5,3), (7,5)$.  For each of those setups, we achieved similarly high fidelities. We will add these experiments to the revised manuscript. 

- Another question raised by some referees relates to the type of noise investigated in our work. When benchmarking the performance of newly proposed hybrid quantum-classical algorithms, the noise from quantum hardware is usually not included in first-case studies. A standard procedure in the field is to benchmark first on shot noise only and consider hardware noise in potential follow-up studies, see, for example, the NFT paper [1] and the other Refs. [2,3] below. The main reason is that quantum noise is strongly hardware dependent, e.g., superconducting quantum hardware is affected by fundamentally different types of noise (in particular CNOT gate noise, decoherence, and measurement noise) compared to, e.g., trapped-ion quantum hardware. Moreover, various error mitigation schemes exist for different types of quantum hardware, which generally require additional error calibration runs on the hardware. An informative study of quantum hardware noise would need to take all these considerations into account and therefore is beyond the scope of the initial proposal of our novel method. 

Due to some potential misunderstandings by a reviewer, we would like to kindly reiterate the main contributions of our work as follows:

- We propose a novel classical kernel, the VQE kernel, which is uniquely suited for the VQE setup. Specifically, its corresponding feature vectors are the basis of the energy functions that can be modeled by the parameterized quantum circuit.
- Leveraging this powerful inductive bias allows us to propose a novel Bayesian optimization method, EMICoRe, which harnesses the fact that for the parameterized quantum circuit, only three observations allow us to determine the energy landscape along an entire line. 
- We demonstrate in our numerical experiments that our proposed method can outperform the current state-of-the-art method NFT and standard BO schemes.       
----------
### Table R1: Results for a new heuristic of setting $\kappa$ using the hyperparameters $C_0$ and $C_1$. Best results highlighted in bold. For energy, lower is better. For fidelity, higher is better.
| Description     | $C_0$  |   $C_1$     | Energy | Fidelity |
|-----------------|-------------------|---------------------|--------|----------|
| In Manuscript   | 0.0               | 1.0    | -5.817205 ± 0.140868  | 0.857839 ± 0.159668 |
| Extreme (small) | 0.1               | 0.1     | -5.816213 ± 0.113959  | 0.857070 ± 0.163497 |
| High (large)    | 10.0              | 10.0    | -5.719012 ± 0.152881  | 0.824516 ± 0.156895 |
| Extreme (large) | 10.0             | 100.0     | -5.703273 ± 0.156716 | 0.804686 ± 0.179101 |
| **Best**   | **0.1**          | **10.0**   | **-5.842853 ± 0.089015**   | **0.869609 ± 0.112713**|

## References

- [1] [Nakanishi K. et al., Phys. Rev. Research 2, 043158 (2020).](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043158)
- [2] [E. Farhi et al., arXiv:1411.4028 (2014)](https://arxiv.org/abs/1411.4028)
- [3] [Carlos Bravo-Prieto et al., Quantum 4, 272 (2020).](https://quantum-journal.org/papers/q-2020-05-28-272/) 

Meta-review: The authors proposed a novel method, known as EMICoRe, to combine Bayesian optimization and prior knowledge of variational quantum eigensolvers (VQE) to enhance the efficiency of the optimization. During the rebuttal, the authors have addressed all the concerns raised by the reviewers. The authors have also provided additional experiments to demonstrate the performance of the proposed method. All the reviewers appreciate the contributions and vote favorably for acceptance. Therefore, I recommend the acceptance of the paper.
Decision: Accept (poster)