Bayesian Parameter Shift Rules in Variational Quantum Eigensolvers
Authors: Samuele Pedrielli
Christopher J. Anders
Lena Funcke
Karl Jansen
Kim Andrea Nicoli
Shinichi Nakajima
Venue: ICLR 2026 Poster (2026)
Keywords: parameter shift rule
variational quantum eigensolver
quantum computing
confidence region
Gaussian process
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Gradient-based optimizer with Gaussian process derivative estimator performs state-of-the-art in variational quantum eigensolvers.

Abstract: Parameter shift rules (PSRs) are key techniques for efficient gradient estimation in variational quantum eigensolvers (VQEs). In this paper, we propose their Bayesian variant, where Gaussian processes with appropriate kernels are used to estimate the gradient of the VQE objective. Our Bayesian PSR offers flexible gradient estimation from observations at arbitrary locations with uncertainty information, and reduces to the generalized PSR in special cases. In stochastic gradient descent (SGD), the flexibility of Bayesian PSR allows reuse of observations in previous steps, which accelerates the optimization process. Furthermore, the accessibility to the posterior uncertainty, along with our proposed notion of gradient confident region (GradCoRe), enables us to minimize the observation costs in each SGD step. Our numerical experiments show that the VQE optimization with Bayesian PSR and GradCoRe significantly accelerates SGD, and outperforms the state-of-the-art methods, including sequential minimal optimization.

Published as a conference paper at ICLR 2026
BAYESIAN PARAMETER SHIFT RULE IN VARIATIONAL QUANTUM EIGENSOLVERS
Samuele Pedrielli 1,2,, Christopher J. Anders', Lena Funcke Karl Jansen 7.8 Kim A. Nicoli, 9.5,6 Shinichi Nakajima 1,2,4+
BIFOLD, Germany, 2 Technische Universität Berlin, Germany 3 Università degli Studi di Padova, Italy, "RIKEN Center for AIP, Japan
5 "Transdisciplinary Research Area (TRA) Matter, University of Bonn, Germany
HHelmholt Institute for Radiation and Nuclear Physies (HISKP), University of Bonn, Germany
Deutsches Elektronen--Snnchrtrro (DESY), Germany
8 *Computation-Be Based Science and Technology Research Center, The Cyprus Institute, Cyprus
9 "Ileendorff Carriers GmbH & Co. KG, Germany
ABSTRACT
Parameter shift rules (PSRs) are key techniques for efficient gradient estimation variational quantum eigensolvers (VQEs). In this paper, we propose their Bayesian variant, where Gaussian processes with appropriate kernels are used to estimate the gradient of the VQE objective. Our Bayesian PSR offers flexible gradient estimation from observations at arbitrary locations with uncertainty information, and reduces to the generalized PSR in special cases. In stochastic gradient descent (SGD), the flexibility of Bayesian PSR allows reuse of observations in previous steps, which accelerates the optimization process. Furthermore, the accessibility to the posterior uncertainty, along with our proposed notion of gradient confident region (GradCoRe), enables to minimize the observation costs in each SGD step. Our numerical experiments show that the VQE optimization with Bayesian PSR and GradCoRe significantly accelerates SGD, and outperforms the state-of-the-art methods, including sequential minimal optimization.
INTRODUCTION
The variational quantum eigensolver (VQE) (Peruzzo et al., 2014 McClean et al. 2016) is hybrid quantum-classien algorithm for approximating the ground state of the Hamiltonian of a given physical system. The quantum part of VQEs uses parameterized quantum circuits to generate trial quantum states and measures the expectation value of the Hamiltonian, i.e., the energy, while the classical part performs energy minimization with noisy observations from the quantum device. Provided that the parameterized quantum circuits can accurately approximate the ground state, the minimized energy gives tight upper bound of the ground state energy of the Hamiltonian.
The observation noise in the quantum device comes from multiple sources. One source of noise is measurement shot noise, which arises from the statistical nature of quantum measurements- outcomes follow the probabilities specified by the quantum state, and finite sampling introduces fluctuations. Since this noise source is random and independent, it can be reduced by increasing the number of measurement shots, to which the variance is inversely proportional. Another source of noise stems from imperfections in the quantum hardware, which have been reduced in recent years by hardware design (Bluvstein et al., 2023), as well as error mitigation (Cai et al., 2023), quantum emror correction (Roffe 2019 Acharya al. 2024), and machine learning (Liao al.) 2024 Nicoli
et al. 2025) techniques. In this paper. we do not consider hardware noise, as is common in papers developing optimization methods (Nakanishi al. 2020 Nicoli et al. 2023b).
Stochastic gradient descent (SGD), sequential minimal optimization (SMO), and Bayesian optimiza- tion (BO) have previously been used to minimize the VQE objective function. Under some mild
"Correspondence {samuel e pedd i e11 1 1e..pdrie11@ccmpppse.... nakaj ma to cooo-berll ie

---

Published as a conference paper at ICLR 2026
Prior
Posterior
f(x)
8
daf(z)
10 2
Figure I: Illustration of our gradient confident region (GradCoRe) approach. Our goal is to minimize [0, )x) where we use GP surrogate f((r) for approximating f (z). Observing f* at points x- and + (green circles) along the d-th direction (solid horizontal line) decreases the uncertainty (dashed curves) not only for predicting ff (x*), but also for predicting 0, ( so that the current optimal point " falls within the GradCoRe (blue square). Our GradCoRe-bar SGD uses the minimum number of measurement shots for achieving required gradient estimation accuracy in each iteration, and thus minimizes the total observation costs over the optimization process.
the true energy / (*xx) over the set of parameters
assumptions (Nakanishi et al. 2020), this objective function is known to have special properties. Based on those properties, SGD methods can use the gradient estimated by so-called parameter shift rules (PSRs) (Mitarai et al. 2018). and specifically designed S.O (Plati 1998) methods, called Nakanishi-Furi Fukanishi-Fuji-Todo Todo (NFT) (Nakamishi et al., 2020), perform one-dimensional subspace optimization with only a few observations in each iteration. lannelli and Jansen 2021 applied BO to solve VQFs as noisy global optimization problems.
Although Gaussian processes (GPs) (Rasmussen and Williams) 2006 have been used in VQEs as common surrogate functions for BO (Frazier 2018), they have also been used to improve SGD-based and SMO-based methods. Nicoli et al. (023a) proposed the VQE kernel-a physics-informed kernel that fully reflects the properties of VQEs and combined SMO and BO with the expected maximumm improvement within confident region (EMICoRe) acquisition function. This allows for identification of the optimal locations to measure on the quantum computer in each SMO iteration. (Tamiya and Yamasaki (2022) combined SGD and BO, and proposed stochastic gradient line BO (SGLBO), which uses BO to identify the optimal step size in each SGD iteration. [Anders et al.) (2024) proposed the subspace in confident region (SubsCoRe) approach, where the observation costs are minimized based on the posterior uncertainty estimation in each SMO iteration.
In this paper, we take a different approach leveraging GPs, and introduce Bayesian parameter shift rule (Bayesian PSR), where the gradient of the VQE objective is estimated using GPs with the VQE kernel. The Bayesian PSR translates into a regularized variant of PSRs if the observations are performed at designated locations. However, our approach offers significant advantages caaantaaes-fexxxiiliiyy and direct access to uncertainty- over existing PSRs (Mitarai et al.) 2018 Wierichs al. (022). More specifically, the Bayesian PSR can observations at any set of locations, which allows the reuse of observations performed previous iterations of SGD. Reusing previous observations along with new observations improves the gradient estimation accuracy, and thus accelerates the optimization process. Furthermore, the uncertainty information can be used to adapt the observation cost in each SGD iteration- in similar spirit to Anders et al. (2044) -which significantly reduces the cost of obtaining new observations, while maintaining required level of accuracy. We implement this adaptive observation cost strategy by introducing novel notion of gradient confidence region (GradCoRe)-the region in which the uncertainty of the gradient estimation is below a specified threshold (see Figure []. Empirical evaluations show that our proposed Bayesian PSR improves the gradient estimator, and SGD equipped with our GradCoRe approach outperforms all previous state-of- of-tthe-rt methods including NET and its variants.
The main contributions are summarized as follows:
We propose Bayesian PSR, a flexible variant of existing PSRs that provides access to uncertainty information.
2

---

Published as a conference paper at ICLR 2026
We theoretically establish the relationship between Bayesian PSR and existing PSRs, reveal- ing the optimality of the shift parameter in first-order PSRs.
We introduce the notion of GradCoRe, and propose an adaptive observation cost strategy for SGD optimization.
We numerically validate our theory and empirically demonstrate the effectiveness of the proposed Bayesian PSR and GradCoRe.
Related work: Finding the optimal set of parameters for a variational quantum circuit a challenging problem, prompting the development of various approaches to improve the optimization in VQEs. Gradient-basee based methods for VQEs often rely on PSRs (Mitarai et al.; 2018) Wierichs et al.) 2022). which enable reasonably accurate gradient estimation of the output of quantum cireuits with respect to their parameters. Nakanishi et (2020) proposed an SMO (Plan) 1998) algorithm, known as NFT, where, each step of SMO, one parameter analytically minimized by performing a few observations. Nicoli et al. (2023a) combined NFT with GP and BO by developing physics-inspired kernel for GP regression and proposing the EMICoRe acquisition function, relying on the concept of confident regions (CoRe). This method improves upon NFT by leveraging the information from observations in previous steps to identify the optimal locations to perform the next observations. Anders al. (2024) leveraged the same notion of CoRe, and proposed SubsCoRe, where, instead of optimizing the observed locations, the minimal number of measurement shots is identified to achieve the required accuracy defined by the CoRe. The resulting algorithm converges to the same energy as NFT with a smaller quantum computation cost, i.e., the total number of measurement shots on a quantum computer. [aamiya and Yamasaki (2022) combined SGD with BO to tackle the excessive cost of standard SGD approaches and used BO to accelerate the convergence by finding the optimal step size. In general context of BO, Müller et al. (202T) proposed gradient information with BO (GIBO) approach, where the uncertainty of the GP-estimated gradient is minimized. Our GradCoRe can be seen as an enhanced version of GIBO, where the theoretically optimal locations are observed with minimum costs based on strong physical information of VQEs.
2 BACKGROUND
Here we briefly introduce Gaussian process (GP) regression and its derivatives, as well as VQEs with their known properties.
2.1 GP REGRESSION AND DERIVATIVE GP
Assume that we aim to learn an unknown function f*() X' + I from the training data X =
(x (II,..., xN)) E N. y (y1,..., BNN E ER,, , (op of R. that fulfills
yn (xx) + En:
M (0,
where No(: n, Σ) denotes the D D-dimensional Gaussian distribution with mean a and covariance Σ. With the Gaussian process (GP) prior p(f(-)) = (P) where 0(.) and k(:,-) are the
prior zero-mean and the kernel (covariance) functions, respectively, the posterior distribution of the function values f = (eof),,.,,(((,,,, X is given as
E RRM at arbitrary test points X' = (a'1,..... x MM E
p( = x) NN ((x.v.o)) : S'x ((x,))), where
(2)
M M(x K' (K + Diag(o)) and S|x,0 = K" - -K" (K + Diag(o)) ()- (3)
K'
are the posterior mean and covariance, respectively (Rasuussen and Williams 2006). Here Diag (v) is the diagonal matrix with specifying the diagonal entries, and K = k(X, X) R N K" ((XX)) E IR N and K" s k(X', X') E R MMM are the train, train-test, and test kernel matrices, respectively, where k(X, X') denotes the kernel matrix evaluated at each column of X and X' such that (k(X, X X') Dm. = k( zmm). We also denote the posterior as p(f()))x, ,y) = GP((O): ((. ()) 8(x.))(( ()))) with the posterior mean [(xx (-) and covariance s[x,o (::)) functions.
Since the derivative operator is linear, the derivative Va. (o.,...., Opʃ) RR of GP samples
0 Since the kernel function corresponds to the also follows GP. Here abbreviate 0₂ Ozs
3

---

Published as a conference paper at ICLR 2026
covariance of GP prior, can straightforwardly handle the derivative outputs by adjusting the kernel so that is consistent with the original kemel defined for non-derivative outputs, Assume that æ is a training or test point with non-derivative output y = (x) + E., and x' and x" are training test points with derivative outputs, y' = 02 (x') E', y² = O. (² + e". Then, the kernel entries involving those three points should be replaced with
k(ar, o²) cov (f))) Ra (xx)) Del
cor(J()))
2.d k(a,, x'),
k( arr""")) cov dx [aa"))
o² DE'
cov( (f(x') [. Oz
192
k( xx"). (5)
The posterior (2) with appropriately replaced kernel entries gives the posterior distribution of deriva- tives at test points. We denote the GP posterior of single component of the derivative as
p p(Ox[())X.)) = GP (aas()) P X v. ()), [x (..))
(6)
with the posterior mean (.) and covariance (:.)) functions for the derivative with respect to zạ. More generally. GP regression can be analytically performed in the case where the training outputs (i.e., observations) and the test outputs (i.e., predictions) contain derivatives with different orders (see Appendix A for more details).
2.2 VARIATIONAL. QUANTUM EIGENSOLVERS AND THEIR PHYSICAL. PROPERTIES
The VQE (Peruzzo et al. 2014; McClean et al. 2016) is a hybrid quantum classicl computing protocol for estimating the ground-star energy of a given quantum Hamiltonian for a Q-qubit system. The quantum computer is used to prepare parametric quantum state |//=), which depends on angular parameters x E X [0, 2T)) This trial state (ea)) is generated by applying D' D) quantum gate operations, G(x))
G o...o G1, to initial quantum state (Wo).. i.e.
[t/a) G(z)) All gates (G D dd'=1 are unitary operators, parameterized by at most one variable {1...., D} be the mapping specifying which one of the variables {xa) parameterizes the d'-th gate. We consider parametric gates of the form G (x) Up (Ta()))
dd- Let d(d')= (1, (1,.... -D'} +
exp (-ixa(a))Pa Pa /2). where Pa is an arbitrary sequence of the Pauli operators {1g. o₂¹ of og fg=1 acting on each qubit most once. This general structure covers both single-qubit gates, such as Rx (x) = exp (-i0o,)) ), and entangling gates acting on multiple qubits simultancously, such as Rxx (z) exp -izoor" X xxx) for qı Q2. commonly realized in trapped-ion quantum hardware setups (Kielpinski et al. 2002 Debnath et al.. 2016).
2
The quantum computer is used to evaluate the energy of the resulting quantum state (Maa by observing
y (x) + E,
where
(*) (v H()) (vol (C))))G(()))o)))))
(7)
and denotes the Hermitian conjugate. For each observation, repeated measurements, called shots, on the quantum computer are performed. Averaging over the number Nshots of shots suppresses the variance x²2 (Nshots) x Nehiot of the observation noise e Since the observation y is the sum of many random variables, it approximately follows the Gaussian distribution, according to the central limit theorem. The Gaussian likelihood ([] therefore approximates the observation y well if o xs x²2 (Nshotss. Using the noisy estimates of fxx obtained from the quantum computer, protocol running on a classical computer is used to solve the following minimization problem:
min (x)),,
thus finding the minimizer ă. i.e., the optimal parameters for the (rotational) quantum gates. Given the high expense of quantum computing resources, the computation cost is primarily driven by quantum operations. As a result, the optimization cost in VQE is typically measured by the total number of measurement shots required during the optimization process We refer to Tilly al. 2022) for further details about VQEs and their challenges.
We do not consider the hardware noise, and therefore, the observation noise e consists only of the measure- ment shot noise.
When the Hamiltonian consists of Nog groups of non-commuting operators, each of which needs to be measured separately, Nahors denotes the number of shots per operator group. Therefore, the number of shots report on the total number of shots per operator group, i.c. the cumulative sum of Nehots over all observations, when evaluating the observation cost.
per observation is Nog x Nehots- In our experiments,
4

---

Published as a conference paper at ICLR 2026
Let Va be the number of gates parameterized by zd, i... Va |(d' E ......'':d d(d") d(d") Mitarai et al. (2018) proved that the VQE objective (7) for Va = 1 satisfies the parameter shift rule (PSR)
(a') f (aos)) - (a' aas)
Va E (0,2m) d= 1..,.... D, a E [0. 2T).
2 sinn a
where {ea)a=1 are the standard basis, and the shift a is typically set to 2. Wierichs et al. 2022
generalized the PSR (9) for arbitrary Va with equidistant observations { x'
2w 2
Ted
2Va-T.
fxf" = 2V -V-a1 (-1)" 5"
(10)
2sin
(2w:13)
4Va
Most gradient- based approaches rely on those PSRs, which allow reasonably accurate gradient estimation from Ed 1 2Va observations. Let
(0) = (7, v2 cos 22 cos 20,..., v2 cos V2O, v2 sin 0, 2 sin 20,... ........ Vin V20)
(1I)
be the (1 (l-dimensional) Va-th order Fourier basis for arbitrary O. Nakanishi et al. (2020) found that the VQE objective function f (.) in Eq. 7 with any G(-), ), H, and (20] can be expressed exactly
fx) -ve
vec (02m) 2b (ma)
(12)
for some E (1+22, where @ and vec(-) denote the tensor product and the vectoriza-
tion operator for tensor, respectively. Based on this property, the Nakanishi-Fuji. Fujiiiii--FiTTooo (NFT)
method (Nakanishi et al. 2020) performs SMO (Platt 1998). where the optimum in a chosen 1D subspace for each iteration is analytically estimated from only + 2V₂ observations (see Appendix B
for the detailed procedure). It was shown that the PSR 9 and the trigonometric polynomial function form (1) are mathematically equivalent (Nicoli et al. 2023a).
Inspired by the function form 12 of the objective, Nicoli et al. (2023a) proposed the VQE kernel
7 > cos( (os(m(na-- (())))aa)))
(2, x²) = og IT
(13)
72+2Va
which is decomposed as (0, x')
with feature maps Φ,, (y2+2V2))73 vec (o))), for GP regression. The kernel parameter 7² controls the smooth-
ness of the function, i.e., suppressing the interaction terms when 72 > 1. When 72 = 1, the Fourier basis (1]) is orthonormal, and the VQE kernel (13) is proportional to the product of Dirichlet kernels (Rudin 1964). The VQE kernel reflects the physical knowledge ([2) of VQE, and thus allows us to perform Bayesian variant of NFT JFT-Bayesian NFT Bayesian SMO-where the 1D subspace optimzation each SMO step is performed with GP (see Appendix [or more details and the perfor- mance comparison between the original NFT and Bayesian NFT). Nicoli et al. (2023a) furthermore enhanced Bayesian NFT with BO, using the notion of confident region (CoRe),
2 _ = E X3; 8|x,o] (z, a²) S x²)
(14)
i... the region in which the uncertainty of the GP prediction is lower than a threshold k. More specifically, they introduced the EMICoRe acquisition function to find the best observation points in each SMO iteration, such that the maximum expected improvement within the CoRe is maximized.
BAYESIAN PARAMETER SHIFT RULES
We propose Bayesian PSR, which estimates the gradient of the VQE objective (7p by the GP posterior
(6) with the VQE kernel (13) along with its derivatives (4) and (5) (which can be explicitly given
as Eqs. 38 and (39 in Appendix C.3). The advantages of Bayesian PSR include: 1) The gradient
estimator has an analytic-form 2) Estimation can be performed using observations at any set of points,
3) Estimation optimal for heteroschedastically noisy observations (from the Bayesian perspective),
as long as the prior with the kernel parameters, 7 and og. is appropriately set, and 4) The posterior
Any cireuit consisting of parametrized rotation gates and non-parametric unitary gates.
5

---

Published as a conference paper at ICLR 2026
V2
Uncersainty dependence ou - -312"22 [z22,"
by O.fie) Observations + General PSR
f(a) of
(xx) 1 Observations + General PSR
a/2
Figure 2: Illustration of the behavior of the Bayesian PSR when Va - (left) and when Va2 (middle). Bayesian PSR prediction (red) coincides with general PSR (green cross) for the designed equidistant observations (magenta crosses). The right plot visualizes the variance (the second equation in Eq. (16)) of the derivative GP prediction at æ', as function of the shift a of observations when Va 1. Intuitively, the minimum uncertainty is achieved with a 2. maximum span (= T) between the two observed points. For all panels, the noise and kernel parameters are set to 0.01, y2 = 9, où 100.
n which corresponds to the
uncertainty can be analytically computed before performing the observations. In Section we propose novel SGD solvers for VQEs that leverage the advantages of Bayesian PSR.
As naturally expected, our Bayesian PSR is a generalization of exisiting PSRs, and reduces to the general PSR (TO) for noiseless and equidistant observations. Let D RR. be the vector with all entries equal to one.
Theorem 3.1. For any x' [0, 2T) and d I..... D. the mean and variance of the derivative GP prediction, given observations
(yo...., 1/2V ,.2v.........- E R2V at 2Va equidistant training points X = (Zo...., 11) E IR Dx2Va for aw x' + 2w+1 Tea with homoschedastic noise 0². 12v, for 0² < 5o, <0n, are
2Va
E3/ (-1)
(d) Pix v.o0 (x')
2
12
+g)g) 0g
((d) (t', 0² (2V3 +0(g)) (15) og
(T2+2V))
2 +V
The proof, the non-asymptotic form of the mean and the variance, and the numerical validation of the theorem are given in Appendix[[] Apparently, the mean prediction (the first equation in Eq. (I5) by Bayesian PSR converges the general PSR (TO) with the uncertainty (the second equation in Eq. (15)) converging to zero in the noiseless limit, i.e., o² -> +0 and hence ye f" (x... In noisy cases, the prior variance of ~ O(a²) suppresses the amplitude of the gradient estimator as regularizer through the first term in the denominator in the first equation of Eq. (15).
Figure illlstrates the behavior of Bayesian PSR when Va (left panel) and when Va = 2 (middle panel). In each panel, given 2Va equidistant observations (magenta crosses), the blue curve shows the (non-deriattive) GP prediction with uncertainty (blue shades), while the red curve shows the derivative GP prediction with uncertainty (red shades). Note the zV. shift of the low uncertainty locations between the GP prediction (blue) and the derivative GP prediction (red). The green cross shows the output of the general PSR (TO) at x' = O, which almost coincides with the Bayesian PSR prediction (red curve) under this setting. Other examples, including cases where the Bayesian regularization is visible, are given in Appendix C
In the simplest first-order case, i.e., where Va 1.Vd = = 1..... D. we can theoretically investigate the optimality of the choice of the shift
Eq. ( (the proof is also given in Appendix C
Theorem 3.2. Assume that Va = 1, Vdl = 1...., D. For any x' [0, 2m)) and d 1, D. the E TR22 at two
mean and variance of the derivative GP prediction, given observations = (v), V2)
training points X (x' a.a, x. + oeaa) E Dx2 with homoschedastic noise (a² o²) are
(dd) P|x "iX.v.o(a'') = (g2 y1 sin
(d)
(at, (x',,²)) = 17722 + 2in Ox
(16)
Again, the mean prediction (the first equation in Eq. (16) is a regularized version of the PSR (9) The uncertainty prediction (the second equation in Eq.. 16 implies that a 7/2 minimizes the

---

Published as a conference paper at ICLR 2026
uncertainty in the noisy case, regardless of o², oy and y (see the right panel in Figure 2 where the variance of the derivative GP prediction at x' is visualized as a function of the shift a of observations for Va 1). This supports most of the use cases of the PSR in the literature (Mitarai et al. 2018). and matches the intuition that the maximum span minimizes the uncertainty.
4 SGD WITH BAYESIAN PSR
In this section, we equip SGD with Bayesian PSR. In the standard implementation of SGD for VQEs, 2Va equidistant points along each direction d = 1, 1..... D are observed for gradient estimation by the general PSR (10) (or by the PSR (9 if Va = 1, Vd) in each SGD iteration.
Bayesian SGD (Bayes-SGD): A straightforward application of Bayesian PSR is replace existing PSRs with Bayesian PSR for gradient estimation, allowing for the reuse of previous observations. We retain R 2Va D latest observations for predetermined R. in our experiments. Reusing previous observations accumulates the gradient information, and thus improves the gradient estimation accuracy, shown in Section 5.2
4.1 GRADIENT CONFIDENT REGION (GRADCORF)
We propose an adaptive observation cost control strategy that leverages the uncertainty information provided by the Bayesian PSR. This strategy adjusts the number of measurement shots for gradient estimation in each SGD iteration so that the variances of the derivative GP prediction at the current optimal point e are below certain thresholds. In a similar fashion the CoRe (14), we define the gradient confident region (GradCoRe)
2|x (1i) = {ar E 3 |x, S N2 Yd
(17)
where = (n]. xb). E RDD are the required accuracy thresholds. Our proposed SGD- based optimizer, named GradCr,,,, measures new equidistant points X = {(x)/) 2w++1 Tea 2V2, D w= d=1 for all directions with the minimum total number of shots such that the current optimal point E is in the GradCoRe (see Figure I
æ + 2V2
Before starting optimization, we evaluate the single-shot observation noise variance o 2 x by collecting measurements at random locations, following Anders et al. (2024). We use this information estimate the observation noise variance as function of the number of shots as 2² (Nahots) Nahotz Let (X',, , o²) be the training data (all previous observations) at the I-th SGD iteration step, and let i E 2V2V be the vector of the numbers of measurement shots at the new equidistant measurement points X for all directions. Before measuring at X the (1 + 1) SGD iteration, solve the following problem:
min 11l s.t. t E (x)(,
(18)
(v) ..... 222D)) T, and xc(t) is the required accuracy dependent on the iteration step t. Informally, we minimize the total measurement budget under the constraint that the posterior gradient variance along each direction di is smaller than the required accuracy threshold. For simplicity, we solve the GradCoRe problem (18) by grid search under the additional constraint that all 2V, D points are measured with an equal number of shots.
where -
We set the required accuracy thresholds to rc() = =x((1)DD, where
(t) max co: D Eb (Ml
(1ix) (r (x (29)
(19)
Namely, k(t) is set proportional to the L2-norm of the estimated gradient at the current optimal point at the t-th SGD iteration, as long it is larger than lower bound. The lower bound co and the slope C1 are hyperparameters to be tuned. This strategy for setting the required accuracy based on the was proposed by Tamiya and Yamasaki (2022).
estimated gradient
In the experiment plots in Section 5 we will refer to SGD-GradCoRe as GradCoRe. Further algorithmic details, including pseudo-code and used hyperparameter values, are given in Appendix D
7

---

Published as a conference paper at ICLR 2026
SGD-128
SGD-256
SGD-512
SGD-1024
Bayes-SGD-125 (ours)
Bayes-SGD-256 (ours)
Bayes-SGD-511 (ours)
Bayes- SSD-1024 (ours)
GradCoRe (ours)
10a
10-1
AEnergy
AFidelity
10a1
10
10a2
0.0
0.2
0.4
0.6
0.8 1. Trial
0.0
0.2
0.4
0.6
0.8 1.0 Trial
Cumulative Nahot x 100 Density
Cumulative Naho x10? Density
Figure 3: Comparison between SGD with PSR (dashed curves) and SGD with Bayesian PSR (solid curves), well GradCoRe (red solid curve), the Ising Hamiltonian with (Q = 5)-qubits [L=3) 33--yeress quantum circuit. The energy (left) and fidelity (right) are plotted as functions of the cumulative Nshots, i.e., the total number of measurement shots. Except GradCoRe equipped with the adaptive shots strategy, the number of shots per observation is set to Nshots 128 (blue). 256 (green), 512 (orange), and 1024 (purple).
SGLBO
Bayes-NFT
EmiCoRe
SubsCoRe
GradCoRe (ours)
10a
10
AEnergy
AFidelity
10-1
10-2.
10-2
0.0
0.2
0.4
0.6
0.8
1.0 Trial
0.2
0.4
0.6
0.8
1.0 Trial
Cumulative Shot x10? Density
Cumulative Na x10 Density
Figure 4: Energy (left) and fidelity (right) achieved within the cumulative number of measurement shots for the Ising Hamiltonian with a (Q = 5)-qubits (L = 3)-layers quantum circuit. The curves correspond to SGLBO (blue), Bayes-NFT (green), EMICoRe (orange), SubsCoRe (purple), and proposed GradCoRe (red).
EXPERIMENTS
5.1 SETUP
We demonstrate the performance of our Bayesian PSR and GradCoRe approaches the same setup used by Nicoli al. (203a)). For all experiments, prepared 100 different random initial points, from which all optimization methods start. Our Python implementation uses Qiskit (Abraham al. 2019) for the classical simulation of quantum hardware. The implementation for reproducing our results is attached as supplemental material.
Hamiltonian and Quantum Circuit: open boundary conditions,
We focus the quantum Heisenberg Hamiltonian with
H = - -ie(Xxx)
+2 h
(20)
8

---

Published as a conference paper at ICLR 2026
where hej the Pauli operators acting on the j-th qubit. For the quantum circuit, a common ansatz, called the L-layered Efficient SU (2) circuit with open boundary conditions, where Va 1, vil (see Nicoli et al. (2023a) for more details).
Evaluation Metries: We compare all methods using two metrics: the best achieved true energy (*),, for (-) defined in Eq. 7. and fidelity ((Gs//ba))* E [0,1]. 1]. The latter is the inner product between the true ground ground-state state wave function (&ss), computed by exact diagonalization of the target Hamiltonian H, and the trial wave function, //ma), corresponding to the quantum state generated by
the circuit using the optimized parameters T. For both metrics, we plot the difference (smaller is better) to the respective target, i.e.,
= f (vos (atcsssH/eaass (bas), AFidelity (Mas (as) (Das (a) = 1 (
A Energy- = (val I (ba) -
(21)
(22)
in log scale. Here, (égs) and (Fes H (as) are the wave function and true energy at the ground respectively, both of which are computed analytically. As a measure of the quantum computation
cost, we consider the total number of measurement shots per operator group (see Footnote ( for all
observations over the whole optimization process.
Baseline Methods: We compare our Bayesian SGD and GradCoRe approaches to the baselines,
including SGD with the PSR (9). Bayesian NFT, SGLBO (Tamiya and Yamasaki [2022) EMICoRe (Nicoli et al. 2023a), and SubsCoRe (Anders et al. 2024). We exclude the original NFT (Nakanishi
et al. 2020) because it is outperformed by Bayesian NET (see Figure 5 in Appendix (B. We also
exclude GIBO (Müleer et al. 202I), which is an even weaker baseline than the original NFT (sce Appendix
All SGD-based methods use the ADAM optimizer with 4, = 0.05, Bs (0.9. 0.999). For the methods not equipped with adaptive cost control (i.e., all methods except SGLBO, SubsCoRe and GradCoRe), we set Nahiot = 1024 for each observation the same setting as Nicoli ct al. -nneess specified explicitly. To avoid error accumulation, all SMO-based methods measure the "center", i.c., the current optimal point without shift, every D + 1 iterations (Nakanishi et al. 2020 Bayes-SGID and GradCoRe estimate the gradient from the R- 2Va D latest observations for R 5, and GradCoRe initially uses the fixed threshold (t) = o²2 /256 before starting the cost adaption after 1) SGD iterations. Further details on the algorithmic and experimental settings are described in Appendix and Appendix respectively.
Algorithm Setting:
5.2 EMPROVEMENT OVER SGD WITH BAYESIAN PSR AND GRADCORE
First, we investigate how our Bayesian PSR and GradCoRe improve SGD. Figure 3 compares SGD with the standard PSR (SGD) and SGD with Bayesian PSR (Bayes-SGD) on the Ising Hamiltonian, i.e., Eq. (20) for (x,
(0, 0, 00, with (Q 5)-qubits (L- 3y-layers quantum circuit. Both for SGD and Bayes-SGD, the optimization performance with Nshots = 128, 256, 512, 1024 measurement shots are shown. The left and right panels plot the difference to the ground-state in true energy (2T) and fidelity
(-1, 0, 0) and (x
achieved by each method as functions of the cumulative Nahots, i..., the total number of measurement shots. To the right of each panel, the trial densiry, i.c., the distribution over the trials computed by kernel-density estimation, after the use of x 10? total measurement shots is depicted. The median, the 25-th and the 75-th percentiles are shown as a solid curve and shades, respectively. We observe that, although Bayesian PSR provides more accurate gradient estimator, as shown in Figure 7 in Appendix F the optimization performance on par with the SGD with the standard PSR. On the other hand, GradCoRe outperforms SGD and Bayes-SGD with different fixed number of shots (Nshota) through the entire optimization process. Note that GradCoRe is built on the Bayesian PSR framework, which provides uncertainty estimation stated in Theorem 3.1 This enables the method to automatically determine the optimal number of measurement shots at each optimization step. The adaptively selected number of shots and the accuracy threshold k(t) for GradCoRe are shown in Appendix F
5.3 COMPARISON WITH STATE-OF-THE-ARI METHODS
Figure 4 compares GradCoRe to the baseline methods, SGLBO, Bayes-NFT, EMICoRe, and Sub- sCoRe. Our GradCoRe, which significantly improves upon SGD as shown in Figure establishes

---

Published as a conference paper at ICLR 2026
itself as the new of exhibiting faster convergence and achieving lower overall energy (see in Appendix [F.T] for statistical significance test results. We also conducted experiments with different Q and L, well as for the Heisenberg Hamiltonian, on which the results are reported in Appendix F.T
CONCLUSION
The physical properties of variational quantum eigensolvers (VQEs) allow us to use specialized optimization methods, i.e., stochastic gradient descent (SGD) with parameter shift rules (PSRs) and a specialized sequential minimal optimization (SMO), called NFT (Nakanishi et al. 2020). Recent research has shown that those properties can be appropriately captured by the physics- informed VQE kernel, with which NFT has been successfully improved through Bayesian machine learning techniques. For instance, observations in previous SMO iterations are used to determinc the optimal measurement points (Nicoli et al. 2023a), and observation costs are minimized based on the uncertainty prediction (Anders et al.] 2024). In this paper, have shown that similar approach can also improve SGD-based methods. Specifically, we proposed Bayesian PSR, where the gradient estimated by derivative Gaussian processes (GPs). Bayesian PSR generalizes existing PSRs to allow for flexible estimation from observations at an arbitrary set of locations. Furthermore, provides uncertainty information, which enables observation cost adaptation through the novel notion of gradient confident region (GradCoRe). Our theoretical analysis revealed the relation between Bayesian PSR and existing PSRs, while our numerical investigation empirically demonstrated the utility of our approaches. We envisage that Bayesian approaches will facilitate further development of more efficient algorithms for VQEs and, more generally, quantum computing. In future work, we aim to explore the optimal combination of existing methods and strategies for selecting the most suitable approaches for specific tasks, i.e., specific Hamiltonians.
ACKNOWLEDGEMENT
The authors thank the reviewers for their constructive comments and discussion for improving the pa- per. The authors also thank Stefan Kühn for valuable inputs and discussions. This work was supported by the German Ministry for Education and Research (BMBF) under the grant BIFOLD25B, the Euro- pean Union's ' HORIZÓN MSCA Doctoral Networks programme project AQTIVATE (101072344) the Deutsche Forschungsgemmensschat (DFG, German Research Foundation) as part of the CRC 1639 NuMeriQS project no. 511713970 and under Germany's 's Excellence Strategy Cluster of Excellence Matter and Light for Quantum Computing (MI4Q) EXC 2004/2 390534769, the Euro- pean Union's 's Horizon Europe Framework Programme (HORIZON) under the ERA Chair scheme with grant agreement no. 101087126, and the Ministry of Science, Research and Culture of the State of Brandenburg within the Centre for Quantum Technologies and Applications (CQTA). C.J.A. is supported by the Bayes duality project, JST CREST Grant Number JPMJCR2112.
REFERENCES
H. Abraham al. Qiskit: An open-source framework for quantum computing. Zenodo, 2019, doi: 10.528 25t2
Rajeev Acharya, Dmitry A. Abanin, et al. Quantum error correction below the surface code threshold. Nature, 2024. doi: 10. 10.10388441586-04--844444
C. J. Anders, K. Nicoli, B. Wu, N. Elosegui, S. Pedrielli, L. Funcke, K. Jansen, S. Kuhn, and S. Nakajima. Adaptive observation cost control for variational quantum eigensolvers. In Proceedings of 41 International Conference Machine Learning (ICML2024), 2024. doi: 0..555/3692000.66221333
Dolev Bluvstein, Simon J Evered, Alexandra Geim, Sophie H Li, Hengyun Zhou, Tom Manovitz, Sepehr Ebadi, Madelyn Cain, Marcin Kalinowski, Dominik Hangleiter, et al. Log- ical quantum processor based on reconfigurable atom arrays. Nature, pages 1-3, 2023, doi:
10. 1 10108/04118888 066--33-66677---.
10

---

Published as a conference paper at ICLR 2026
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. Rev. Mod. Phys., 95:045005, Dec 2023. doi: 10.1 ww..10//ew.Moohy..55.....
Debnath, M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe. Demonstration of small programmable quantum computer with atomic qubits. Nature, 536(7614))63-66, 2016. doi: 10.1038/nature 18644.
Roger Fletcher. Practical methods of optimization. John Wiley & Sons, 2000.
Frazier. A tutorial on Bayesian optimization. ArXiv e-prints, 2018. doi: 10.48550/arXiw.1807, 02811.
Giovanni lannelli and Karl Jansen. Noisy Bayesian optimization for variational quantum eigensolvers. ArXiv e-prints, 2021. doi: 104485504axii.22112000666
D. Kielpinski, C. Monroe, and D. J. Wineland. Architecture for large-scale ion-trap quantum computer. Nature, 4 41766890) 2002. doi: 10. 10.1033/aature00784..
Haoran Liao, Derek S. Wang, Iskandar Sitdikov, Ciro Salcedo, Alireza Seif, and Zlatko K. Minev. Machine learning for practical quantum error mitigation. Nature Machine Intelligence, 6(12): 1478 1486, November 2024. ISSN 2522-5839. doi: 10. 010888444225--004-000077777 URLhttp
7/dx. ord.x.oo.org/1..oo 103.://....00///22//
024 00927-2
Jarrod R McClean, Jonathan Romero, Ryan Babbush, et al. The theory of variational hybrid quantum- classical algorithms. New Journal of Physics, 18(2):023023, 2016. doi: 10.1088/1337-2660/111 023023.
K. Mitarai, M. Negoro, M. Kitagawa, al. Quantum circuit learning. Phys. Rev. 98:032309, 2018. doi: 10. 10.1100//P111 103//0hsseee.AA.... 9. A.9.0033399.
Sarah Müller, Alexander
Rohr, and Sebastian Trimpe. Local policy search with bayesian optimization. In A. Beygelzimer, Y. Dauphin, P. Liang, and Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
Ken M. Nakanishi, Keisuke Fujii, and Synge Todo. Sequential minimal optimization for quantum- classical hybrid algorithms. Phys. Rev. Res., 2:043158, 2020. doi: 10. 10.103//hysReeeesserc..... 043158.
K. A. Nicoli, C. J. Anders, L. Funcke, T. Hartung, K. Jansen, S. Kuhn, K.-R. Müller, P. Stornati, P. Kessel, and S. Nakajima. Physies-informed Bayesian optimization variational quantum circuits. In Advances in Neural Information Processing Systems (NeurIPS2023), 2023a.
K. A. Nicoli, C. J. Anders, et al. EMICoRe: Expected maximum improvement over confident regions. https: /ttp::///tttht com/ emweore/ceeeeeeeee eewie///mmiccsse 2023b.
Kim A. Nicoli, Luca Wagner, and Lena Funcke. Machine-leaming-enhaeee optimization of noise- resilient variational quantum cigensolvers. ArXiv e-prints, 2025, doi: 10.48550%arXi..551...700
Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, al. A variational eigenvalue solver photonic quantum processor. Nature Communications, 5(1):4213, 2014. doi: 10. 00.103///comm55213
John Platt. Sequential minimal optimization A fast algorithm for training support vector machines. Microsoft Research Technical Report, 1998.
C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, USA, 2006. doi: 00.7551//mmppess3320.000111 00t.eeeese/.220.....000.00
Joschka Roffe. Quantum error correction: An introductory guide. Contemporary Physics, 60(3): 226-245, 2019. doi: 10. 10..0000010711..00..66667777
Walter Rudin. Principles of Mathematical Analysis. McGraw-Hill, 1964. doi: 1017/
11

---

Published as a conference paper at ICLR 2026
Shiro Tamiya and Hayata Yamasaki. Stochastic gradient line Bayesian optimization for efficient noise-robust optimization of parameterized quantum circuits. upj Quantum Information, 8(1):90, 2022. doi: 10. 00008//41134-022-00555200
Jules Tilly, Hongxiang Chen, Shuxiang Cao, al. The variational quantum eigensolver: A review of methods and best practices. Physics Reports, 986:1-128 2022. doi: https:/ddoi.or//0011016 10.1/// physrep.2022.08.005
David Wierichs, Josh Izaac, Cody Wang, and Cedric Yen-Yr Yu Lin. General parameter-shift rules for quantum gradients. Quanium, 6:677, March 2022. ISSN 2521-327X. doi: 10.22331/ 9-2022-03-30-677.
12

---

Published as a conference paper at ICLR 2026
A GENERAL GAUSSIAN PROCESSES (GPS) WITH DERIVATIVE OUTPUTS
The derivative GP regression can be straightforwardly extended to the case where both training outputs (i.e., observations), and test outputs (i.e., predictions) contain different orders of derivatives.
Assume that we have a set of input points, and for each input point x E RR, the corresponding output, i.e., observation prediction, is [(x) Oz, (xx), where Ong Far
0 Let us denote the derivative
kernel functions
k(2x Oe', k( x0) O (2r 20)
d = 0, = 0, d = 0, = =,....,
[dd,d') (x,x²)
if d 1, ...., D,d' =0,
Pz x²)
d = I...., D,d = I....,..
For training points X (x(()) I== and test points X' = {m/(m) }M Jm=1, should set the the entries
of the train-train K E R. xN train-test K' E IR xM and test-test K" E IR MMM kernels as
Kn, [(d(ma))d(an)))) (In; Fnt),
(23)
K2 1 k(d(ea) (I Lmm)
(24)
K" mmmm [(d(am)) (Im, Bmm)
(25)
where
{0 if the corresponding output for the input æ is f(x),
d(x) = (d if the corresponding output for the input is fxx).
(q and 3 with the kernel matrices K, K', K" set as Eqs. 23 25
give the posterior GP for
the corresponding test outputs.
For higher-order derivative outputs, we can define the kernels in exactly the same way as above, by applying the same derivative operators to the kernels as the ones applied to the outputs, i.c..
k( x') = ((m)))
)(ro)
gbrbo d/
k(xxx'))
a(r) ... (i..))))))) f(æ) and airi) the corresponding outputs at æ and a' are
g(rb) f(x'), respec-
tively, where of7) denotes the r-th order derivative with respect za. Bal
or
B NAKANISHI-FUI-TTRRE (NFT) ALGORITHM (NAKANISHI ET AL., 2020 AND BAYESIAN NFT
Let (eald be the standard basis. NFT initialized with a random point T0 with first observation 70 = (*( + So. and iterates the following procedure: for each iteration step 4,
1. Select an axis fl..... D) sequentially and observe the objective y E 2V at 2Va points X = (#1...., 22V2) = {@ + ea) 2V2 E IR Dx2V, along the axis d Here a [0, 2 is such that a. + 0, aw # aw for all and w' + w.
2. Apply the ID trigonometric polynomial regression F(() b +b, (0) to the 2V, obser- vations y. together with the previous best estimated score g and analytically compute the new optimum é' x' 1-1 + Ōe., where = argming f()).
3. Update the best score by T² I(0).
Note that if the observation noise is negligible, i.e., y = f*x),, each step of NFT reaches the global optimum in the 1D subspace along the chosen axis d for any choice of a,, and thus performs SMO
"With slight abuse of notation, use the set notation to specify the column vectors of matrix, i.c., ( Øx) { 1a=1-
13

---

Published as a conference paper at ICLR 2026
NFT-128 Bayes-NFT-128
NFT-256
NFT-512
NPT-1004
Bayes-NFT-256
Bayes-NFT-512
Baycs-NFT-1on
10a1
AEnergy
AFidelity
10-2
10"
0.0
0.2
0.4
0.6
0.8
1.0 Trial
0.0
0.2
0.4
0.6
0.8 1.0 Trial
Cumulative Nshot
x 10 Density
Cumulative Nahot
x10? Density
Comparison between NET (Nakanishi et al. 2020) and Bayes-NFT for the Ising Hamilto- nian with a (Q 5)-bbbs (L 3)-layers quantum circuit. The energy (left) and fidelity (right), in the forms of Eqs. (21) and (22). respectively, are plotted as functions of the cumulative Nahots. i.e., the total number of measurement shots. The number of shots per observation is set to Nohnts 128 (blue). 256 (green), 512 (orange), and 1024 (purple).
Figure 5:
exactly. Otherwise, errors can be accumulated in the best score Y., and therefore additional measurement may need to be performed at H' after certain iteration interval.
Bayesian NFT (Bayes-NFT) performs the 1D trigonometric polynomial regression and optimization in Step2 with GP with the VQE kernel (13). where all previous observations are used for training. Using previous observations allows prediction with smaller uncertainty and thus more accurate subspace optimization. Figure comparees the original NFT and Bayesian NFT on the Ising Hamiltonian with
a (Q= = 5) qubits (L = 3y-layers quantum circuit with different number of shots per observation. We observe that using GP generally accelerates the optimization process.
C
PROOFS
Here, we give proofs of theorems in Section 3 and numerically validate them.
C.I PROOF OF THEOREM
We start from a more general theorem than Theorem 3.1 which is proven in Appendix C.3
Theorem C.1. Assume that, for any given point ã [0. 2m) we have observations y (Yo 112 E &RR2Va at 2Va equidistant training points X (arg...., 3224-11 E Dx2va for 2wr+1 a= = 2 2V2 xea with homoschedastic noise o a². 12v, E VV, Then, the mean and variance of the derivative (af(x') prediction at x' = I + for any d = 1...., and a' E [0, (0)2z) are given
1 1(x "xXx),,(('')) (x') = (3*+2Va)o +2V 2V--1 1 (-1) yuu
(dd
con(Vea')
Vasiul (2w+1)e -VV++1/2)a')
4V2 cos Via (22+2Va0="00 (26)
+
2xin
(2w-11 11s
-a(/2)
sin( 4V -2/22)
.V.
(d)
(x', x='====) 1x.o
= o² V2(V=+1)(VV++1)
1V cos (2Va2)))))
3(172222300*/06--2
(6=442Vaj0 (
8V (cun (722a +
(27)
Regardless of the observations, the predictive uncertainty 27 is periodic with respective to a' with the period of T/Va. We can easily get the following corollaries.
14

---

Published as a conference paper at ICLR 2026
1/2 =2, =0.10, 7² =3.00, o3 =10.00
V/a =2, =0.50, ;" =3-3000, =5.00
0, f(a) I an GP
0. f(a) 1 any GP
far) my Amalytic
far) my Amalytic
Observasion
Obsevvatios
General PSR
General PSR
2r
0
2x
C.1 under parameter settings (see above each panel). Figure 6: Numerical validation of Theorem
Given the 2Va equidistant observations (magenta crosses), the derivative GP prediction (blue curve) with uncertainty (blue dashed curves) compared their analytic forms (26) and (27), i.c., the mean function (red curve) and the variance function (red dashed curves), respectively. We observe that our theory perfectly matches the numerical computation. The green cross shows the prediction by the general PSR (10), which almost coincides with Bayesian PSR prediction when o /0/ 0.01 (left panel), while significant difference observed when o² /²/o = 0.1 (right panel).
Corollary C.2. For the test point at x' = Et i.e., a' = O, the mean of the derivative GP prediction is
22 (-1)"
Vala 22vg2=vv)
(d) P v.o (x')
2
4V4
(28)
(T242Va)o² /a0
2V
Corollary C.3. For the test point at x' = 3 + a a'ed, Va' =, Va, 2x Va...., ,2Va 1) /V. variance of the derivative GP prediction is
the
(d) ox Ix
(x', x') o²
Valv=Dirav
(29)
Ignoring high order terms with respect to o /2/oß in Eqs.(2 and 29 gives Theorem 3.1
Figure 6 shows numerical validation of Theorem C.I where the derivative GP prediction (bluc curve) with uncertainty (blue dashed curves) is compared to their analytic forms, i.c., the mean function (66 (red curve) and the variance function (27) (red dashed curves), respectively, under two settings of noise and kernel parameters. We observe that theory perfectly matches the numerical computation. When o² / 0.01 (left panel), the regularization is small enough and the Bayesian PSR prediction (red curve) almost coincides with the general PSR prediction (green cross). On the other hand, when o² /a = 0.1 (right panel), the Bayesian PSR prediction (red) does not match the general PSR prediction (green cross), because of the regularization.
In addition, we can also derive an upper bound of Eq 29). which is useful for setting the grid search range the GGD-GradCoRe algorithm described in Appendix [
Corollary C.4. For the test point x' -c a a'ea, Va' æ VV, 2n /Va. (2Va 1)xz/ //Va the variance of the derivative GP prediction is upper bounded as
(d) ox
(x', x') σ2
22 6
(30)
Proof. Let - (n2 =2 /0 > Then, Eq. 29 can be written as
fxx00 . (x', x²) = o² Va X222 ()(V++)22V++++
D2s 4V
8
(31)
3(E+2V
and its derivative with respect to is given as
83
(x',x') -02
(32)
3(E+2V]"((+446.)
15

---

Published as a conference paper at ICLR 2026
Apparently, (q. is negative for any Va 1, and therefore, the predictive variance s(d) monotonically decreasing with respect to g > Taking the limit 5 0, Eq. reduces to the bound (30).
i[x,o] (a' x²)
Note that the monotonicity of the predictive variance with respect to Ç, and hence with respect to of (when 7 and o² are fixed), matches the intuition that the posterior uncertainty is smaller when the prior variance is smaller.
C.2 MATHEMATICAL PREPARATIONS
Before proving Theorem C.I
we give some mathematical identities on the trigonometric functions.
C.2.1 ROOT OF UNITY
For natural number N E (1,2.... let us define a root of unity PN = 2wi/N such that PN = 1. Then, the following hold:
N-1
PX
1-px
for
k: 1,..., N 1,
(33)
1-03
ΣN- (x+ = PN
PN
ké En=o PN = PN
kb 1- 1-
1-p2
for
k ,,..., N 1,
(34)
also holds for even N that
N
(k- N/2)N
En PN
(n+1/2)ktnN/2
PN
k/2
PN
n(k+N/22
PN
k/2 1
PN
N/2) 0
(35)
1 PN
n=0
for k 1, ,.... N/2-1 1.
C.2.2 PROPERTIES OF DIRICHLET KERNEL The summation in the Dirichlet kernel can be analytically performed as
1 1+ 2 > cos (nx) = 1 1+ 2- n=1
inz
2
- N
=e INx 1- +12N+1)
1-ote
-N441/13) - (N41/22)
-in((N++122z)
(36)
sin(x/2]
Therefore, also holds that
2 > n sin (nx) = EVE (1/Va + cos (nx))
dz
(1+ 11+27 11+2 V22 v= cos (nz)
dz.
(N+1/2) cool (N sin( /2) sin((N+1222z) (1x_1/1)))))))x())2)))))221))))) cos(1r/2)
snn2
N Neoo((N (NN_NN++//2osssn(oo))---aaa
sin (=//)
sin(Nx)
(NN
(37)
2sin (x/2)
sin(a/22)
C.3 PROOF OF THEOREM
For derivative predictions (af('') Oaf(x") the kernels should be modified Eqs. and (5).
For the VQE kernel (13). they are
((x x') = O,; k(xxx')
22 sin (M(mna x(5))
7242 Ye coa(m((r,r (1()),)=((()))
= =g
Hase
(38)
7 72VV
7² +VV
16

---

Published as a conference paper at ICLR 2026
k(0', x²) 022 boy
2 Va v² coo cco(p1t's221)
++2 Vi coofeless- 12.1)
=
Ila
(39)
72=2Va
7=22V
The training kernel matrix for (Ti
E
2
2Va-1 is Toeplitz Ted
2V.
0=
T2V2-2
-2V2-1
71
Tg
7
K o
T2
71
TO
RVVax2V
T2V.-2
TU
71
T2V,-1
71
TD
where
2-22 V/4
cas
TV
2T
(40)
+242V
For test point at x' # + a'e, the test kernel components are
Kg
R1
K2V
K0' on:
where
Yt sin (f 2mmt x-
fww
+2V2
(41)
The first identity
33
for the root of unity implies that
2V
2V
for
w 1,..., 2Va
and therefore
21
2T
2Va
for
w 0, 2V
COS
VW
(42)
2V.
for
W 1,
2V
m=0
2V
2T
sin
VW
0
for
0, ..... ,Va..
(43)
2V
The second identity
34 for the root of unity gives
2
21.
(u+1/2)
IV
Σ
(20 (2++1)
2V
for
w 1,....21 2V
and therefore
21
2V,
for
0.
Σ
cos
(20 + 1)w
-2V4 for
2Vd,
(44)
2V
for
w = 1, 1.... 2Va - 1,
2V
(2v + 1)w π
sin
for
0. 2Va
(45)
2Va
==0
17

---

Published as a conference paper at ICLR 2026
The third identity 35 for the root of unity gives
2V
2V.
e (20+1)
Σ
(-1)" (2 =0
o=0
for w l,..., Va - 1, and therefore
-V111
(-1) cos
(2v 1 2Vd
0 for
w .,..., Vd.
(46)
2V (-1)² sin
2Va for
w = Vd.
(2v 1)w
2V = [ for
w = 0, ..... Va 1.
(47)
1=0
From the symmetry of the trigonometric functions, it holds that
2T
for for
= 1, 2Va
COS
UW
(48)
2V
W = 2, 6, 6......... 2V
Va
2V
2T
2m
sin
vw
Σ
sin
FW
(49)
2V
2V
v= V
Note that the factor -1 in the odd w case in Eq. (48) is because the summand is -1 for Va. while the summands for the other are canceled each other.
By using Eq. (48) Eq. can be written as
Tuo
2
72 + 2
V2
COS
2V.
2nT
7
+
2
for
for
W 0,
1,3, 5, ................................ ,2V 1,
for
2,4,6 2V, 2,
+2
and therefore
o
K
(2V2I2v, + (72-1)10 + +cc²)
7/2 + 28o
o
2V2 + +(1 c) (+-1 1 (1 c) ))
(50)
72 2Vo
where
2V
With the training kernel expression (50,, the matrix inversion lemma gives
(0² + (V2+2Valoo/088 +V2) I2va + (1 2) (08²1 Da a) )²
[K + I2I2v -1= 22+2V2
of
+2V of
1 (o422
(1 <) (>2-1 (3p 9)a 077
I2va + (2242V
7242
%6
(22 ( (1)
9)
12va (24220)0+100 (1
18

---

Published as a conference paper at ICLR 2026
1 (1 c) (
9) a |7}
12+ + (1
e)²
+2V
(22=2Va)=(10tt34 (1 c) (72- 1)
1251
2Va(?2 1) 2V1)) (1 8)
11 + 1o2
+22 %
Iava - (1 c) 72 (201
(672 + +V++VVV00 00 + 2Vn
(02 + 2V2)0 +4V, )0
aa T)7
1 a + b11 + ccc
of
where
72 + 2V
a=
(n2 + 2V o + 2V
72 -1
b
(72 + 2Va) oo 2V
c= (n2 + 2V2)02 + 4V6
For the test kernels
K0
K1
K. 0
12Va-1
2 22 Val x² 72 + +V2
o VV
T²"' = 2
3(7² + 2Va)
with
2 2Σ
Vir
w= usin (v (2m++ 1)s
- -a'))
2V
Kw
22 + 2Va
(51)
have
V2 usin(
Vd 12 aa 2V
1 112 = o/ > 21
w==
F2Va
2V4-1
Va
V/
vo sin (
(2w+1) 2V
a²))
sin (v'
(2m++))
a'))
2V
21. 1
V2 V4
00
(72 + 2V
1==0
(0- (0-21) 2²)
(2w+1)a
a')) cos ((v (o+ 0')
(2w+1)5 - -a)))
cos
2V
2V,
19

---

Published as a conference paper at ICLR 2026
V2
V2
2
0
Σ
vD
(72 + 2Va)²
0=1 v'=1
w=0
cos
(u
0²)
(2w+1)m
- -a/)) -cos
cos
(w++)))
(2w+1)x a'
2Va
0 (72 + 2V
Va
2Va-1
vo
v=1
w=D
(2t0+ (') xr
cos ((v r')a') sin (2w+1)(e-- r'y sin ((v ea)a')
co8
2
2Va
cos
(2wtt) cos ((v: + 1')a') - sin (2w+1) 2Va
sin ((v: + 1')a')
2V
(52)
V4
2
0
(2nr+ 1)(0- cos ((V r')aa)
VU
CO8
(72 + 2V)2
2
(2w+1)(tte)) cos ((V + (o)oo))
cos
2Va
(53)
(73 + 2V32 2(2Vg)
o
|Σ 2
V2
+ V22 cos (2Vaa')
(54)
of (7² + 2V 2(2Vg)
ValVa + 1)(2V) + 1)
+ V2 cos (2Vad')
6
4V3
(Va + 1)(2Va + 1) + Va cos (2Vaa') )
(72 2V,))
Here we used Eqs. and (45) obtain Eqs. and 54
from Eq. (2)
We also have
214
2
v=1 usin (v (2w+1)z a') 2Va
V2
-T
1k 11 K
12v3 = 0
x2 + 2V/
w=
2 2>
Vu
21 1==
sin (v
(2w++) -
-a'))
of
1=
2V
72 + 2Va
2 >> Via
E2V=-1
sin
(2wr+ Djur
2V 008 vat cos
(2m++ Dew
sin va'
2Va
0
72 + 2Vo
0,
and
2V/-1
2 2Σ V/ -==1
sin (v
12ar+1) a')
-T
c = of Σ (-1) w
2Va
72 + +Va
0=
2
V/
mm
10 0
2V
( sin (w
(2w++)1
2Va
- a')
00
72 + 2Va
2 V v -2Vi-1 (-1) (sin (2w+1)ux 2V4 cos vat cos
(2w+1)ex
sin va
2Va
00
72 + 2Va
2V22Va cos Via' 922 2V
o
20

---

Published as a conference paper at ICLR 2026
00 4V2 cos Vaa'
72 + 2V2
Here, we used Eqs. (6) and 47) in the second last equation. Therefore, the mean of the derivative is
(1))
"X ".x..o( (2') =k -/T
(K + +
=k TT
[Ivv. + bl2v. + ccc²) y
0
/T/
cc 4) 12v, L2 + ck TIT bk
TT
k
y +
of
I siu (u( (2utlle -a') (
2V61
4V3 ens Via' 2 (-1)" 72+2V
g +2V2
(2 valn (0( (Lw
2Vo-1 Yhu
-V
cos Vaay))
a
7242Va
+
32+2Va
2V
Σ Vw
Vi Σ
usin (v
(
12 -))) -
4cV cos Vaaa
+2V
w=0
v=1
E
yiu
{2
V2 Vi
usin (o(
(2w lin =a())}
28
-V3(-1)"" van'' (roaaa)ot/ooaa4a
(55)
22=+2V2)0= +2Va
Eq. 37
implies that, for w = ,9,1,...,,,,, - 1, holds that
2
2 sin (v (2w++1)= a') od))
sin (V2
(2
24
-a'))
V/ ens (V++1/22 (2w 10m -a'))
28
#=1
2 sin"
(2m+1) -a 12)
244
sin
-aa) )2) 2V
sin
(2w+
Viaa)
Va coe
(2u
(2wi+ -(V++1/2)a'
2 sin 2
(20++)
24 -aa (2)
sin
(20+11)
-a') -2)) 2V
sin
(-1) Vaae)
Vi2 cns (-1)a
(2u++13n
4V4
-(Va+1/2)a
2 sin
-a') 12) 2V4
(2m+1) -a') 12) 27
1-)
(oc()a 2siu
+
Vasin
(2m 1)- -VV2+112)0'
4V
(2ui+ 1)m
-a 02122
(2m=1)x aa'²"'" 4V.2
4V.
(56)
Substituting Eq. (56) into Eq. 55 gives Eq. (26).
The posterior variance can be computed as
((d) $ [xxr (2', x¹) = 1
k -T
K
+ 0² I2v.)) K
= k
TT
(I2v, + bl2v. L2v, + occ K
= k.
[
-
122 + b(k 12v,)² + c(k' T.T c)))
oVVa(V=+1)22Vat)
3(7"+2V4)
av?
(V++1](V+++1)
+ Va cos (2Vaa') + coo (
(VV cus Via 2
0
(7*+2V+)
72+2Va
o2V +(7++
AV? a.a
(V++
+22+2
6
0,00
4v2 cas( (2V
16V2 cos? Vyn' (22+)V2)=
( (T2+2
o2V2(V++1112Vat1)) 3(72+2V4)
3² +2VA
1V2
(V++1)(2V++1)
-222 2/0))2/0++000 +2V (22+2)))²
2-
4v3 cns( 2V
8V/11+cos 2V (2++)
4224 1 22
21

---

Published as a conference paper at ICLR 2026
0² 12
1V2 coo (2V2aa)
8V (can 2V Va
(57)
which gives Eq. 27
C.4 PROOF OF THEOREM 3.2
In the first order case with Va = 1, vd= = ,,...,,, D. the test VQE kernels for predicting derivatives
Sas(x') (a'),,a[(()) Oxl("")) are
k(, x²) = oî [2sin (Idd - =2/2)
22 2 cos (I x'de)
k(x,, x²) = 0z'y
II
72 +
72 +2
d'
0²
cos (2'a2 x'd))
222 + 2 cos (I'R x'd')
k(0', x²)
k(0', x²") o/
dn Ow'd
7/2 +2 2
7/2 +2
d'
Then, the kernels with the two training points X = (x' wved x' + wee) and the one test point x' are
1
+2oss
m = 20isin a
(-1
20 T²= = 722 2
K= o0
7242
k
+2 201
1
72 +2
7242
With these kernels, the posterior mean is
1 ( ) =k -k= (K + IOoIN) y
-id)
-IT
22in a
1+ 1+ 108
12+2 20
-1
72 +2 (-1
1)
+2
+2 cos 20
1 1 12022/00
7-
+ 0² %o
22+2 cos 20t
2 (-1 1) 7242
1
+2
(2+02/033- (12100202000 2
(2+02/03)3
+2 20
+ 0//0
42
+2
2sin a
1
(-1 l)
72 2 (1+ 02 02/03)
y2+2 cos 2a
7242
1/2 1/1
=2sin a
01+02/08))
+ 2) - (72 + 2 cos 20x
(J2 - y/ sin
(n2/2 0 1002 /03 + 2sn
The posterior variance is
(())
"xx x0) = K"
k'
T
(K + 0² IN) - k
20
Ao sin? 2 a
1 1 02 /og
22 cas 2a
+ 2)22 (-1 1)
3"t
72 +2
+2 cos 2a
1 1+ 120/
1
+2
20
Aog sin? 2
1
(-1 1)
722 +
(72 + 2)²
(1+ (1+02/03)
(22-2cos 2c
1
72
20
4sin² a
72 +2 2
(n/2 + 2)02 103 - 2 cos 20
20
(12 + 2)0² 2//000
72 +
(n2 2)0² 2 + + a
a²
(n2/2+ 1)02 / 0 + 2 sin 2
Thus we obtained 16).
22

---

Published as a conference paper at ICLR 2026
Algorithm 1: (SGD-GradCoRe) Improved SGD algorithm using VQE-derivative kernel GP with the GradCoRe measurement selection subroutine, as described Algorithm The algorithm finds the minimum number of shots required to estimate the gradient wrt. parameter configurations of the quantum circuit to optimize with SGD. The optimization stops when the total number of measurement shots reaches the maximum number of observation shots allowed, ie., Moot-shots To avoid cluttering notation, the algorithm is restricted to the case where Va Generalization arbitrary Va is straightforward.
Input
a:: initial starting point (best point)
Parameters:
Va =
number parameters to optimize, i..., kD"
Niutshoss shoe-soss Total # of shots, i.c., maximum allowed quantum computing budget.
C²2 measurement variance using single shot.
Initial GradCoRe threshold step = O.
Timi Number of steps in beginning to initial GradCoRe threshold No-
co smallest allowed GradCoRe threshold
ci
GradCoRe threshold scaling parameter
Output
x' optimal choice of parameters for the quantum circuit.
n + 0/+ initialize consumed shot budget
f 0/x initialize aptimization step
n² l /* initial to use for Timial steps
x", y", o + O, O, /) initialize empty Gaussian process
while Nurssoon do
/+ choose measurement points & number of shots a.t. a' is in the GradCoRe of x²
X, gradcore_measurenee
(X', vi, a' a x²) (Algorithm
for i (1, ..., } do
+ quantum circuit (parameters-X,, shots-Wd measure chosen points
8: t at2
10
11
end
9. ô (81, ... 8 x (8) (ō1...., )(2)) / concatenate observed values noise
x tt rt (x*, X), (v' g) (a" σ) / add now observations to Gaussian process
13
t - #' p o0th ((1)))))()))) / SGD (or variant) step using
derivative
if 2 Thiis then
(7i *+ l max co: Exn ("x ( ).11)(())))) adapt GradCoRe
threshold
end
18
t 1 update the step
19
+ n + Eg Pa / update the consumed shot budget
end
return x'
D ALGORITHM DETAILS
D. GRADCORE PSEUDO-CODE
Algorithm I describes GradCrrCCoRe in detail. SGD-GradCoRe uses the GradCoRe measurement selection subroutine described in Algorithm 2 which selects measurement points and respective minimum required number of shots estimate the quantum circuit parameter derivative required for
23

---

Published as a conference paper at ICLR 2026
Algorithm 2: (GradCoRe measurement selection subroutine) Select the points to measure and respective minimum number of required shots such that when updating the GP with these new measurements, the GP's derivative uncertainty at the current best point is smaller than the threshold K, i.e., the current point within the GradCoRe.
Input
X, y/ Gaussian process at current step
current best point
K (Ki....,,i (DD GradCoRe thresholds at current step
Parameters:
Va
measurement variance using single shot.
ā shift from best point the previous step (default
Output
points which should be measured and added the GP compute derivative. number shots measured points.
begin
for
E (1, ..., D) do
Xa (# ea, à ea) /* choose points measure along
ã= + Ka /x initialize measurement noise to minimum (most expensive,
& F')
for i [V2ka, do
/* create temporazy GP copies, add points with observation noise
X', y' e (x,X,,) Xa), (y 0,0), (a,
/* find largest observation noise for which a is in the GradCoRe
if x x) S x (0+ > 8) then
84
12
end
13
/* compute shots from variance through single shot
variance x²
14
end
X XI,..., XD) /+ concatenate points to measure
15
16
D - (DI,..., DD) / concatcnate shots to measure per point
17
return X, 7
18
end
the SGD. Note that the grid search range in Line 5 in Algorithm Appendix C
2 is set based on Corollary
C.4 in
D.2 PARAMETER SETTING
Every algorithm used in our benchmarking analysis has several hyperparameters to be set. For transparency and to allow the reproduction of our experiments. we detail the choice of parameters for EMICoRe, SubsCoRe and GradCoRe Table The SGLBO results were obtained using the original code from Tamiya and Yamasaki (2022) and we used the default setting from the original paper. For NFT, Bayes-NFT and Bayes- SGD runs, we used the default parameters specified in Table For algorithmic efficiency, all Mayesian-SMOO methods use the inducer option introduced in Nicoli ct al. (2023a). retaining only the last 2V. D
= 399 measured points once more than 2V. D D 439 points were stored in the GP, where we chose R = 5. Since the discarded points are replaced with a single point predicted from them, the number of the traininig points for the GP is kept constant at R 2V, D 400. On the other hand, Bayesian-SGD methods measure (at most, in the SGD-GradCoRe case) 2V,D = 80 points per SGD step, and we retain R. 2V, D 400 points after more than (R 1) 2V, D 480 points measured. Unlike the Bayesian-SMO
24

---

Published as a conference paper at ICLR 2026
Table 1: Algorithm specific parameter choice for EMICoRe, SubsCoRe and GradCoRe for all experiments (unless specified otherwisc).
Algorithmic specifie parameters
-macq-params
EMICoRe params
as in Nicoli ct al
(003aa)
func
ei
Base acq. func. type Optimizer type # of candidate points evaluation points Gradient strategy for k Order of gradient norm
optim
emicore
pairsize (Jsa)
20
gridsize (Joc)
100
corethresh-strateee
grad
pnerm
corethresh (x)
256
CoRe threshold k # initial steps with fixed Coefficient Co updating k Coefficient C1 for updating g Stabilization interval in SMO steps
corethresh_wice (Tinitial)
40
coremin_scale (Cg)
2048
corethresh_scalee (C1)
1.0
@ebbllizzeinterree
41
samplesize (Nxc)
100
# of MC samples
smo-steps (TNFT)
# initial NFT steps Sequential direction choice
smc-axis
True
--acq-parana
SubsCoRe params
in Anders et al. 2004
optim
subscord center
Optimizer type Alg type SubsCoRe Gradient strategy for k
readout-strategy
corethresh-strategy
grad
pnorm
Order of gradient norm Initial Natos for CoRe initial steps with fixed k Coefficient Co for updating x Coefficient Ci for updating Stabilization interval in SMO steps
corethresh (x)
256
gorethresh_widat (Tiitial))
coremin saanninnscaaleeee (C))
2048
worethresh_scalle (C1)
1.0
stabilize
stabilize_interval interval
coremetric
readout
Metric to CoRe
--acq-params
GradCoRe params
this pape
optim
gradcore grad
Optimizer type Gradient strategy k
corethresh-strategy
рnоrm
2
Order of gradient norm Initial Naws for CoRe # initial steps with fixed Coefficient
corethresh (x)
256
corethresh width (Timnaa)
coremin scale (Ca)
2048
updating k Cocfficient C₁ for updating k
corethresh_scale (C1)
1.4
coremetric ir gooptim
readout 05
Metric to set CoRe learning rate for SGD Optimizer for SGD
adam.
methods, do not add additional inducer based on the prediction from the discarded points, and therefore the number the training points for the GP is kept constant at R. 2V, = 400.
E EXPERIMENTAL DETAILS
As discussed the main text, our experiments focus on the same experimental setup as in Nicoli et al. (203a)) and Anders et al. (2024). Specifically, starting from the quantum Heisenberg Hamiltonian, we reduce it to the special case of the Ising Hamiltonian at the critical point by choosing the suitable couplings, namely
Ising Hamiltonian at criticality: (-1.0, 0.0, 0.0); = (0.0, 0.0, - 1.0).
a...a., "readour" in Nicoli et al. (2023a)
'All hyperparameters specified the table are the default in Nicoli al. 2023a
25

---

Published as a conference paper at ICLR 2026
Table 2: Default choice of circuit parameters and kernel hyperparameters for all experiments (unless specified otherwise).
Deafult params
gbits
qubits
-lllaeers
# of circuit layers
-circuit
esu2
Circuit name
-bbc
False
Open Boundary Conditions max number readouts
n-iter
1x10**7
kernel
vqe
Name of the kernel
- kernel parzel-params params
Bayes-NFT
EmiCoRe
SubsCoRe
GradCore
Bayes-SGD
gamma sigma_D 0
3
3
3
3
10
10
10
10
10
It is important note that due to the finite size of the system at hand, this choice of parameters does not imply criticality but already represents a challenging setup, as discussed in Sec. 1.2 in Nicoli al. (2023a). We stop the optimization when predetermined maximum number of cumulative shots (total measurement budget on the quantum computer) reached; unless specified otherwise, we set this cutoff Nmma 107.
Our implementation of GradCoRe can be found in the supplementary zip file and will be made available on Github upon acceptance. In experiments, the kernel parameters dp and are fixed to the values in Table Furthermore, NFT, Bayes-NFT, Bayes-SGD, SubsCoRe and GradCoRe require fixed shifts for the points to measure at each iteration. our experiments, we always used for SMO based methods (as this makes the uncertainty uniform in the 1D- DD-subspace, discussed for SGD based methods (as this minimizes the uncertainty in the noisy case, discussed in Section 3 unless explicitly stated otherwise.
in Anders et al. (2024) ) and a
Each experiment shown in the paper was repeated 100 times (trials) with differently seeded starting points. We aggregated the statistics from these independent trials and presented them in our plots. We used the same starting point for every algorithm in each trial to ensure fair comparison between all approaches. Note that SGD-based methods do not require measurements at the starting point, but SMO-based methods do. Therefore, each starting point is further paired with a fixed initial measurement.
All experiments were conducted on Intel Xeon Silver 4316 2.30GHz CPUs.
F DETAILED BEHAVIOR OF BAYESIAN PSR AND GRADCORE
Figure shows the gradient esitmation error by PSR (dashed curve) and Bayesian PSR (solid curve) during the SGD optimization. We clearly observe that Bayesian PSR estimates the gradient more accurately than PSR.
Figure 8 shows the behavior of the GradCoRe threshold k(t) (left), and the number v(t) of measure- ment shots (left) that GradCoRe used each SGD iteration.
ADDITIONAL RESULTS
Figures 9 11 show the optimization result for the Ising Hamiltonian with (Q = 7)- (L = 5)- layers circuit, that for the Heisenberg Hamiltonian with (Q 5) -)uuiisss (L 3a-laye circuit, and that for the Heisenberg Hamiltonian with (Q = 7)-qubits (L = -ayles circuit. We observe that the advantages of GradCoRe
prominent for larger qubits.
To confirm the statistical significance that the GradCoRe outperforms the baselines, we applied the Wilcoxon signed rank test. Table 3 shows the p-values for the null hypothesis that there is
26

---

Published as a conference paper at ICLR 2026
PSR
Bayesian- (ours)
0.5
0.4
0.3
0.2
0.1
0
25
50
75
100
SGD Steps
Figure 7: Gradient estimation error by PSR (dashed curve) and Bayesian PSR (solid curve) for
Nshots 1024, evaluated by the 1.2-distance between the estimated gradient jz()) and the true gradient g (.) (computed by the PSR with simulated noiseless measurements).
# Measurement Shots
60000
0.5
0.4
k(
40000
0.3
0.2
20000
0.1
0
200
400
600
800
1000
1200
200
400
600
800
1000
1200
SGD Steps
SGD Steps
The GradCoRe threshold x(t) (left), set according to Eq. (9)), and the number of measurement shots (right) per SGD iteration used by GradCoRe. As expected, the number of shots gradually increases as the GradCoRe threshold decreases, reflecting the flatness of the objective function via the gradient norm estimation.
Figure 8:
difference in performance between GradCoRe and the baselines. In all cases, we observe p-values smaller than 0.05, proving the advantage of GradCoRe.
G PERFORMANCE EVALUATION OF GRADIENT INFORMATION WITH BO (GIBO)
As mentioned in Section 1 the gradient information with BO (GIBO) (Müller et al. 2021 can be considered general BO baseline. Here we compare its performance with other baseline methods. For GIBO, we use the VQE kernel (Nicoli et al., (023a) but do not use the strong prior information about the optimal locations of observed points given by Theorem 3.2 Instead, we minimize the trace of the predictive variance, as acquisition function, following Algorithm in Müller et al. (2021). Specifically, we sequentially optimize 2D points by using the L-BFGS-B (Fletcher (000) optimizer in each SGD iteration. For choosing each point, we start from 10 D initial points drawn from the uniform distribution on the D-dimensional sphere with radius 5 = 0.2 centered at the current best point, and choose the best one (see Table in Appendix A.8 in Müller et al. (202I)). Figure 12 compares the performance of GIBO with the other baseline methods evaluated the main text, i.c., SGD-PSR and NFT with Nohot 1024, well our proposed GradCoRe. We see that GIBO is outperformed by both the the plain SGD-PSR and NFT, even though it uses some physical prior information through the VQE kernel.
The typical solvers for VQEs, including the tfate-of-the-ee art methods, use the strong physical knowl- edge for determining the observed locations they observe the points T + 3 ea for 1,...; : D.
27

---

Published as a conference paper at ICLR 2026
Bayes-NFT
EmiCoRe
SubsCoRc
GradCoRe (ours)
10"
10P
AEnergy
AFidelity
10-1
10-1
10-2
0.0
0.2
0.4
0.6
0.8
1. Trial
0.0
0.2
0.4
0.6
0.8
1.0 Trial
Cumulative Nshot
x x00 Density
Cumulative Nahot
x107 Density
Figure 9: Energy (left) and Fidelity (right) achieved with the cumulative number of measurement
shots for the Ising Hamiltonian with (Q 7)-uubiss (L 5)-layers quantum eircuit.
Bayes-NFT
EmiCoRe
SubsCoRe
GradCoRe (ours)
10-1
AEnergy
AFidelity
10
10-2
0.0
0.2
0.4
1.0 Trial
0.0
0.2
0.4
0.6
0.8
1.0
Trial
Cumulative Nahot
x10? Density
Cumulative Nahot
x10? Density
Figure 10: Energy (left) and Fidelity (right) achieved with the cumulative number of measurement shots for the Heisenberg Hamiltonian with (Q = quubbisss (L. 3y-lrs quantum circuit.
Bayes-NFT
EmiCoRc
SubsCoRe
GradCoRe (ours)
10a
AEnergy
AFidelity
10
0.4
0.6
0.8
Trial
0.0
0.2
0.4
0.6
0.8
1.0 Trial
Cumulative Nohor
x10 Density
Cumulative Nato x100 Density
Figure 11: Energy (left) and Fidelity (right) achieved with the cumulative number of measurement shots for the Heisenberg Hamiltonian with (Q = 7)-qubsss (L = 5)-layers quantum circuit.
In this paper, we proved in Theorem 3.2 that this choice is optimal for minimizing the uncertainty of gradient estimation. Since we use the VQE kernel also for GIBO, its poor performance is due the inaccurate optimization of the acquisition function in D-dimensional space for 2D points, which is known to be challenging (Frazier] 2018). On the contrary, our GradCoRe observes the theoretically optimal locations, based our Theorem 3.2,, instead of tackling the challenging optimization of the multipoint acquisition function. This in principle may be achieved with GIBO equipped with an oracle mulli-points acquisition function optimizer, which does not exist. This further demonstrates that
28

---

Published as a conference paper at ICLR 2026
GIBO
SGD-PSR
NFT
GradCoRe (ours)
10P
AEnergy
AFidelity
10
10
10
10
0.0
0.2
0.4
0.6
0.8
1.0
Trial Density
0.0
0.2
0.4
0.6
0.8
1.0 Trial
Cumulative Nshot
x100
Cumulative Nahot
x107 Density
Figure 12: GradCoRe.
Performance comparison between GIBO (with VQE kernel), SGD-PSR, NFT, and
combining existing techniques developed separately in the machine learning and physics literature can enhance performance.
Table 3: Obtained p-values by performing the Wilcoxon signed rank test. Each column indicates the baseline against which GradCoRe is compared, showing the p-values obtained in different settings. All p-values are smaller than 0.05, implying consistent performance improvement by GradCoRe over all settings.
Ennerg
AFidelity
Experiment
Bayes-NFT
EmiCoRe
SubsCoRe
Bayes-NFT
EmiCoRe
SubsCoRe
Ising (3,5)
9.68e-05
4.23e-09
3.84e-05
1.80e-02
1.42e-07
5.19e-04
Ising (5,7)
1.16e-07
4.02e-10
7,49e-10
1.54e-05
4.06e-09
2.89e-08
Heis. (3,5)
8.57e-11
8.19e-11
2.29e-12
5.30e-06
2.78e-07
3.16e-10 10
Heis. (5,7)
2.37e-16
5,84e-16
2.82e-17
1.31e-02
2.69e-02
8.20e-03
H CLASSICAL COMPUTATION CoST
Following common practices as in previous work on VQE optimization, we focus only on the quantum computation cost when evaluating the optimization performance in Section S Here we discuss that the classical computation cost still negligible even when the gradient is estimated by GP. The cost of GP prediction is dominated by the Cholesky inversion of the kernel matrix, which scales O(N3) in the number of training points N. Therefore, the classical computation cost for estimating all elements of the gradient by Bayesian SGD and GradCoRe is O(D N:3), where the maximum number of
training points is upper-bounded by N R d D 400 (see Appendix D.2 for the choice of N at each step).
Table 4 reports on the classical computation time per step in comparison with the SGD with the standard PSR. Although prediction slows down Bayes-SGD by two orders magnitude, the computation time is still small, and thus it is reasonable to ignore it in the context of VQE optimization.
Algorithm
SGD (with standard PSR) (1.4 0.2) 10
Bayes-SGD
(tstep) + to Otttep
(4.3 + 0.7) 10
Table 4: Wall-clock classical computation time per step for SGD with the standard PSR and Bayes- SGD for the Ising Hamiltonian, simulated using the Efficiens SU (2) ansatz with Q 5 qubits and L = 3 layers. Average and standard deviation are computed over full optimization run, using 107 shots.
29

---

Published as a conference paper at ICLR 2026
I GAUSSIANITY OF OBSERVATION NOISE
Here we empirically demonstrate that the observation noise approximately follows a Gaussian distribution. Figure 13 shows the distribution of energy observations with Nahot 1, 128, 1024 at randomly chosen fixed point in the 40-dimensional parameter space of Efficient SU (2) circuits with Q = 5 qubits and L 3 layers. Note that the center of the distribution is shifted the origin by subtracting the average of the energy observations. Since every single measurement (shot) a sum of multiple terms, the observation noise nearly follows a Gaussian distribution, even for
Nahot
News=1
Neve 88
News 128
New 1024
0.23
Immprtical
Caraddan 2801
0.7
Gazanten (-08.731)
limpitical
Simptrical
Gacamben (o18.177)
Impprtica)
Gauadaan (yn 68)
0.29
9:5
2.5-
0.35
Density
22 4 Dumstyy
3.5
2.30
0.3
0.65
0.5-
0.25
6.1
Deviattin from
Devibtton from moan
Deviation moan
Devattion
mean
Figure 13: Observation noise distribution for Nahot 8, 128, 1024 at a randomly chosen fixed point.
30


Reviews:
[{"review_id": "nM3z2SXS63", "reviewer": "Reviewer_pQ8e", "summary": "This paper introduces the Bayesian Parameter Shift Rule (Bayesian PSR), a method that uses Gaussian Processes (GPs) with a physics-informed VQE kernel to estimate the gradient of the VQE objective. A key feature of this Bayesian estimator is its flexibility, allowing it to reuse observations to improve the accuracy of the gradient estimate. The paper also proposes a novel adaptive cost strategy called the Gradient Confident Region (GradCoRe), which leverages the gradient's posterior uncertainty provided by the GP. At each SGD step, the GradCoRe strategy solves for the minimum number of measurement shots needed to ensure the gradient's uncertainty falls below a dynamically updated accuracy threshold. Numerical experiments show this combined SGD-GradCoRe method accelerates optimization, as measured by the number of shots on the quantum side, and outperforms existing state-of-the-art methods.", "questions": "Figure 3  shows that Bayes-SGD (with data reuse) performs no better than standard SGD. This implies the performance gain of GradCoRe comes from its adaptive shot strategy, not from the Bayesian PSR or data reuse. Do the authors agree? If so, should the paper be re-framed to de-emphasize the Bayesian PSR and focus on GradCoRe as an adaptive shot-counting heuristic that can be applied to any PSR?   \n\nPlease justify the use of a different kernel ($\\gamma=1$) for GradCoRe versus the baselines ($\\gamma=3$). This is a critical confounding variable that is not mentioned in the main text. Does GradCoRe still outperform SubsCoRe if both use $\\gamma=3$? Or if both use $\\gamma=1$?\"\n\nSince the practical contribution of this work is the non-insignificant reduction in number of shots, can you please expand a bit on the computational overhead (on the classical side) introduced by your model. Presumably this should scale as $\\mathcal{O}(N^3)$ and you have chosen N=400 as mentioned in the appendix. The major drawback that VQE suffers from in terms of scalability is the amount of time required for training a model with a substantial number of parameters. A discussion on this computational tradeoff could really help your work.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper's primary contribution is building on recent work by Nicoli et al. and Anders et al. to adapt their concepts to a more general framework. It takes the Confident Region (CoRe) idea, which was previously applied to 1D subspace methods, and adapts it for the D-dimensional gradients used in SGD, creating the Gradient Confident Region (GradCoRe). \n\nThe quality of this adaptation is well-supported on two fronts: theoretically, it provides a solid proof (Theorem 3.1) that its Bayesian PSR estimator is a valid generalization of the standard Parameter Shift Rule, and empirically, it demonstrates superior performance against strong, state-of-the-art baselines, including SubsCoRe from Anders et al that can be seen as a precursor model to the current one. This leads to a reduction of  the total number of required measurement shots by approximately 10-20% to reach the same energy levels in the benchmarks. The paper also maintains a high degree of clarity, using its experiments (Fig. 3) to effectively untangle the two sources of its performance gain. The results in Figure 3 indicate that while data-reuse alone (Bayes-SGD) provides no benefit over standard SGD, the entire performance lift is achieved by the adaptive measurement strategy (GradCoRe), which uses the GP's uncertainty to minimize shot cost.\n\nIts significance comes from successfully broadening the utility of these physics-informed, adaptive-cost Bayesian methods beyond specialized algorithms, making them accessible to the widely-used SGD framework", "weaknesses": "The paper's novelty is somewhat limited, as its main contributions are largely an incremental adaptation of the VQE kernel  and Confident Region (CoRe) concept  developed by Nicoli et al. and Anders et al. While adapting this to the more general SGD framework is a valuable step, it relies heavily on these recently established foundations. \n\nThe paper's primary ablation study (Figure 3) reveals that Bayes-SGD, which incorporates the Bayesian PSR and data reuse, offers no clear performance benefit over standard SGD. This suggests that the paper's performance gain is derived exclusively from the GradCoRe adaptive shot-counting strategy. This undermines the presented value of the Bayesian PSR framework itself.\n\nAdditionally, a limitation in the empirical validation is the inconsistent use of the gamma kernel hyperparameter; Bayes-SGD use gamma=1 while all key baselines are set to gamma=3. The paper lacks an ablation study to quantify the performance sensitivity of all methods to this parameter, making it difficult to fully isolate the architectural gains of GradCoRe and Bayes-SGD from a potentially more favorable hyperparameter setting. \n\nLastly a very minor thing. I would argue that the paper overstates the novelty of Theorem 3.2, which proves the optimality of the π/2 shift. The proof primarily serves as a validation that the Bayesian framework aligns with existing knowledge, rather than being a new discovery.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "9ooeLb3owg", "reviewer": "Reviewer_TZCS", "summary": "This paper presents a Bayesian variant of the well-known parameter shift rules (PSRs). Bayesian PSRs enable flexible gradient estimation from observations at arbitrary locations and unify the standard PSRs. Connecting the Bayesian PSRs with SGD, the paper proposes a new optimizer named gradient confident region (GradeCoRe), which outperforms many STOA methods for VQE training.", "questions": "- Is there an intuitive explanation for why the uncertainty dependence on $\\alpha$ is a \"U\" shape curve (as in the last panel in Figure 2)?\n- Even with Bayesian PSRs, it is a bit counterintuitive that an SGD-based algorithm can outperform Bayesian optimization methods. Is there an explanation for this?\n- Is it possible to establish the optimality argument in the second order? If not, what are the challenges?\n- Does the momentum information (e.g., momentum SGD + Bayesian PSRs) further accelerate the optimizer?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Bayesian PSRs are a natural generalization of the standard PSRs, and they enable more efficient exploitation of previous measurement data by taking uncertainty information into account. \n- The theoretical relationship between Bayesian PSR and existing PSRs is established. Bayesian PSR reduces to standard PSRs in certain cases. Moreover, the optimality of PSRs is established. \n- Leveraging this new gradient Bayesian-based gradient estimation technique, gradient-based optimizers demonstrate improved numerical performance in various VQE problems.", "weaknesses": "- The current Bayesian PSR is derived based on Gaussian processes. Is this an appropriate assumption for VQE? How to justify the GP model?\n- Is there a quantitative way to compare Bayesian PSRs and standard PSRs in terms of the quality of gradient estimation and sample complexity?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ACcBc5o9lU", "reviewer": "Reviewer_9nAn", "summary": "The authors propose a method for optimizing variational quantum eigensolvers. The method aims to reduce the number of measurement shots by estimating gradients using Gaussian Processes (GPs). This Bayesian approach allows the algorithm to reuse previous observations in its gradient estimates. The authors combine this with an adaptive observation cost strategy, \"GradCoRe\", and show that this combined method outperforms comparable approaches for a fixed measurement shot budget.", "questions": "* Is it sensible or standard to assume that the observation noise in the model is Gaussian??\n* The experimental setup (Ising model, $(Q=5)$, $(L=3)$) is taken from Nicoli et al. (2023a). Is this 5-qubit problem considered a sufficient benchmark for \"state-of-the-art\" claims in the field?", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The problem it targets—reducing measurement shots—is high-impact and central to near-term quantum algorithms.\n* They obtain rigorous estimates for the mean and covariance of the gradient GP (Theorem 3.1) and show that their Bayesian PSR converges to the standard PSR in the noiseless limit.\n* In a simplified setting (Theorem 3.2), they use this estimate to demonstrate that the common $\\pi/2$ shift is optimal for minimizing uncertainty, justifying a technique used in the literature.\n* The proposed method is a clever combination of existing tools (SGD, GPs, and physics-informed kernels).", "weaknesses": "* The method appears to be a combination of existing methods (which, as noted, is also a strength), making the novelty of the contribution somewhat unclear.\n* The experiments seem fairly limited, primarily focusing on a $(Q=5)$-qubit system. A comparison on a higher-dimensional problem would be needed to substantiate the \"state-of-the-art\" claim.\n* The graphs comparing techniques (e.g., Figure 4) are hard to parse. I would suggest experimenting with ways to represent the uncertainty without using shaded fills, which tend to obscure each other.\n* The derivation of the derivative GP (Section 2.1) is essential to the method. While the paper states the kernel modifications and alludes to Appendix A, a clearer, self-contained (or better-referenced) derivation in the main text would be beneficial.\n* While the paper focuses on quantum cost (measurement shots), it neglects the *classical* computation cost. The method requires GP regression and solving an optimization problem (Eq. 17) at each step. An empirical comparison of this classical overhead would be helpful.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TdCqMVvZZX", "reviewer": "Reviewer_fLQm", "summary": "This work proposes the Bayesian parameter shift rule (PSR), which uses Gaussian process with VQE kernel, to estimate the gradient of objective in variational quantum eigensolver. Furthermore, guided by the posterior uncertainty, the stochastic gradient decent (SGD) step is completed in the gradient confident area. Numerical results demonstrate that the proposed optimization strategy can achieve faster convergence than the vanilla SGD.", "questions": "Please see the Weakness section", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The manuscript provides a comprehensive theoretical analysis and numerical results, which are convincing.", "weaknesses": "- The proposed method is a natural application of Bayesian optimization with Gaussian process to VQE gradient estimation, which offers limited novelty.\n- The motivation of employing the Bayesian PSR is not clearly clarified. It would be more convincing to discuss the challenges the proposed method tries to address.\n- The manuscript lacks a formal statement or theorem guaranteeing the existence of an optimal sequence for both the shift parameter and shot allocation that ensures the fastest convergence.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)