Adaptive directional gradients for parameterised quantum circuits

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Adaptive directional gradients for parameterised quantum circuits

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Adaptive directional gradients for parameterised quantum circuits Brian Coyle,1, 2 Snehal Raj,3, 4 Virag Umathe,1 El Amine Cherrat,5 and Elham Kashefi1, 3

arXiv:2606.09734v1 [quant-ph] 8 Jun 2026

School of Informatics, University of Edinburgh, Edinburgh, UK 2 Fujitsu Research of Europe Ltd., United Kingdom 3 LIP6, CNRS, Sorbonne Université, Paris, France 4 QC Ware, Palo Alto, USA and Paris, France 5 Quantum Signals, Paris, France

Training parameterised quantum circuits (PQCs) on quantum hardware is bottlenecked by the measurement cost of gradient estimation, which under the parameter-shift rule scales linearly in the number of trainable parameters and dominates the total shot budget of training at scale. In this work, we propose a framework of forward gradient estimators for PQCs, based on the forward mode of automatic differentiation, that yields an unbiased estimator of the gradient by averaging a freely tunable number of random directional derivatives and recovers SPSA, random coordinate descent, and the parameter-shift rule as limiting cases, with no ancilla qubits or controlled-gate overhead. We prove that stochastic quantum forward gradient descent converges under standard assumptions, with an explicit second-moment expansion that interpolates between the single-direction extreme of SPSA and the full-gradient extreme of parameter-shift. Within this framework we derive QUIVER (Quantum Iterative V-adaptive Estimator Rule), an adaptive optimiser for parameterised circuits whose update rule follows from a closed-form minimum measurement-cost allocation. We show numerically that forward gradients train Hamming-weight-preserving orthogonal quantum neural networks with up to 60 qubits and 1770 parameters on the ECG5000 and MNIST datasets orders of magnitude more efficiently than the parameter-shift rule. We also demonstrate that our proposed QUIVER optimiser can outperform iCANS and gCANS measurement-frugal optimisers on optimisation problems using the quantum approximate optimisation algorithm and quantum simulation with the variational quantum eigensolver.

INTRODUCTION

Gradient-based training is the dominant paradigm for parameterised quantum circuits (PQCs) [1, 2], from the variational quantum eigensolver (VQE) in chemistry to quantum neural networks (QNNs) in quantum machine learning (QML). As problem sizes grow, so does the number of trainable parameters, a trend compounded by recent evidence that quantum models, like their classical counterparts, benefit from overparameterisation [3– 5]. The cost of extracting gradient information from such models is itself a function of the ansatz expressibility [6], the measurement cost per gradient step therefore grows as O(N ) under the parameter-shift rule, placing gradientbased PQC training under increasing pressure as models scale [7]. Classically, the efficiency of the backpropagation algorithm [8] and its generalisation into automatic differentiation (AD) pipelines [9] ensures that gradient computation scales modestly with the number of parameters. No such general efficiency exists for quantum models. All known approaches to efficient quantum gradient extraction require either fault-tolerant circuits with deep ancilla overhead [10], architectures with restricted expressibility [11–13], or approximate gradient estimators that remove the explicit O(N ) parameter dependence from the gradient rule itself. The third class is attractive for near-term devices, but the parameter dependence is not eliminated. In SPSA [14], a single random perturbation replaces N parameter-shift evaluations, yet the N dependence re-emerges in the estimator variance. Sim-

ilarly, random coordinate descent (RCD) [15] computes one exact gradient component per step, but convergence slows by O(N ). In this work, we show that these and other approximate estimators are special cases of a single framework: forward gradients [16]1 . Forward gradients reconstruct the full gradient from V random directional derivatives, where V is a free parameter independent of N . By varying V and the direction distribution, one recovers SPSA (V = 1, Rademacher), RCD (V = 1, basis vectors), and the parameter-shift rule (V = N , basis vectors) as limiting cases. For classical neural networks, forward gradients have shown backpropagation-competitive performance without the backward lock [19] or memory overhead [20, 21]. Quantum circuits admit no known analogue of backpropagation, making the forward mode a natural gradient strategy for this setting. Casting SPSA, RCD, and the parameter-shift rule as limiting cases of this framework provides a unified treatment of estimators that have until now been analysed separately, and exposes V as an explicit lever through which one can interpolate between the cheapest and most expensive gradient strategies. We make three contributions. First, we establish a convergence bound for stochastic gradient descent with the

1 Closely related ideas appear under different names in the classi-

cal literature, including directional gradient descent for reinforcement learning [17] and sketched gradient methods for large-scale optimisation [18]; these differ in scope but share the idea of estimating gradients from a small number of random projections.

2 central-difference forward gradient estimator on PQCs, incorporating both the O(ε2 ) finite-difference bias and quantum measurement shot noise (Section IV), which to our knowledge is the first such result for this estimator class on parameterised quantum circuits. The bound reveals an exact no-free-lunch result on convex losses: the V -fold saving in per-step cost is exactly compensated by a V -fold increase in the number of steps required, so the total shot cost is V -independent. The bound also identifies the finite-difference step size ε as the dominant hyperparameter: shot noise is amplified by 1/ε2 , and the bias–variance trade-off induced by ε is empirically more consequential than the choice of V . Second, we derive quiver (Quantum Iterative Vadaptive Estimator Refinement), an adaptive optimiser built on the forward gradient estimator (Section VII). Extending the CABS framework [22] and the quantum iCANS/gCANS family [23, 24] to the directionalderivative setting reveals a structural property absent from the parameter-shift rule: measurement noise concentrates uniformly across random Rademacher directions, neutralising the per-direction shot-allocation lever that makes iCANS effective. Maximising expected gain per measurement shot over the remaining lever (the number of directions V ) yields a closed-form joint (V, M ) rule that is variance-optimal (Rademacher directions uniquely minimise the estimator second moment at fixed V ) and saturates the Cramér–Rao lower bound for gradient recovery from a quantum shot oracle up to a constant that vanishes as N → ∞. We demonstrate numerically that quiver outperforms iCANS and gCANS on VQE ground-state optimisation and QAOA MaxCut benchmarks. Third, we validate the forward gradient estimator empirically at parameter counts N ∼ 102 –103 , training orthogonal quantum neural networks [12, 25, 26] with up to 60 qubits (N = 1770 parameters) on ECG5000 time-series and MNIST image classification (Section V). Across four problem classes, a fixed V ≪ N estimator reaches the accuracy of the parameter-shift rule at a fraction of its total shot budget, with the cost saving growing with N . We also show heuristically that a decaying V -schedule outperforms both fixed-V endpoints at matched budget (Section VI); we present this as a motivating observation rather than a systematic result, and as the empirical precursor to the quiver treatment of Section VII. The remainder of this paper is organised as follows. Section II reviews gradient descent, automatic differentiation, and quantum gradient extraction including the parameter-shift rule and resource-frugal optimisers. Section III–Section III C propose forward gradients for quantum circuits and develop finite-difference directional derivatives. Section V and Section VI validate the estimator on classification, ground-state estimation, and combinatorial optimisation benchmarks, and Section VII develops the gain-theoretic analysis, noise concentration result, and the quiver optimiser.

II. A.

BACKGROUND

Gradient descent and automatic differentiation

The goal of gradient-based optimisation is to find parameters θ ∗ = arg minθ L(θ) that minimise a loss function L(θ) computed from the output of some model f (θ). In gradient descent, one moves iteratively in the direction of steepest descent with step size η: θ t+1 = θ t − η∇L(θ t )

(1)

where θ ∈ RN and [∇L(θ)]j := ∂L(θ)/∂θj . For vectorvalued losses the gradient generalises to the Jacobian J ij := ∂Li/∂θj . In practice we can only compute ese(θ)2 of the true gradient, ideally unbiased : timators g E[e g (θ)] = ∇L(θ). In classical ML the stochasticity arises from data subsampling; for quantum models, measurement shot noise introduces a second source. The most widely used variant, Adam [27], replaces the raw gradient with exponentially weighted first and second moment estimates: θ t+1 = θ t − η p

c m

, bt + ε v t t−1 c := β1 m c m + (1 − β1 )e g (θ t ), t

t−1

b := β2 v b v

+ (1 − β2 )e g (θ )

(2) (3) (4)

computed elementwise. In practice, learning rate scheduling η → η(t) is also applied. The efficiency of gradient descent depends on how the gradient is computed. In modern machine learning, automatic differentiation (AD) [9] is the standard approach, implemented natively in JAX [28], PyTorch [29], and TensorFlow [30]. AD operates in two modes: the forward mode computes Jacobian-vector products J v (yielding one column of J per pass), while the reverse mode computes vector-Jacobian products v ⊤ J (one row per pass)3 . Since most ML losses map RN → RM with N ≫ M , the reverse mode is overwhelmingly preferred: a single pass computes ∂Li /∂θj for all j simultaneously. This is the backpropagation algorithm [8, 31, 32], though it requires storing intermediate activations and is inherently sequential. We will return to the forward mode in Section III A, where it plays a central role.

Computing quantum gradients

For parameterised quantum circuits, the situation is fundamentally different. A trainable unitary U(θ) :=

2 We use g e interchangeably throughout. e(θ) and ∇L

3 In JAX, these correspond to .jvp() and .vjp() respectively.

(a) Variational/QML algorithm training e gradient feedback ∇L

digit = 7

|0⟩ |0⟩

ECG = ‘anomaly’

Machine learning

|0⟩

Uenc (x)

U (θ)

L(θ)

|0⟩

loss function

|0⟩

e (b) Gradient estimators ∇L

(i) Parameter-shift

(ii) SPSA

cost 2N M

cost 2M

Optimisation

Simulation

(d) QUIVER adapts V and M

Vt = 8

(iii) RCD

(iv) Forward V =2

cost 2M

cost 4M

Loss

lower loss at fixed budget target accuracy

SPSA RCD Parameter-shift Forward (this work)

N/V fewer shots

Total measurement shots

Vt = 6

explore

8 Vt 64 2

Vt = 3

exploit

training step t

80 60 Mv 40 20

FIG. 1. Overview of the forward gradient framework and the quiver optimiser. (a) Variational/QML algorithm training loop: a parametric quantum model U (θ) (preceded by a data-encoding unitary Uenc (x)) is measured to produce a loss L(θ) suited to a range of tasks (quantum simulation, combinatorial optimisation, machine learning). A gradient estimator e feeds back into U (θ) to update its parameters. (b) Gradient estimators ∇L, e ∇L shown as random directions probed in the N -dimensional parameter space at fixed per-direction shot count M . The parameter-shift rule (PS) evaluates 2N circuits per step covering all N basis directions (cost 2N M ). SPSA draws a single random direction (2M shots). RCD draws a single basis-aligned direction (2M shots). The proposed forward gradient estimator (this work) draws V random directions (2V M shots), with V as a free parameter interpolating between SPSA (V = 1) and PS (V = N ). (c) Practical training advantage (cartoon): forward gradients with V ≪ N reach the same loss as PS in approximately N/V times fewer total shots, and achieve lower loss than SPSA and RCD at fixed budget. (d) The quiver optimiser adapts the number of directions Vt and shots-per-direction Mv over training, where t indexes the training step and v a particular sampled direction: early in training, many noisy directions are sampled at low shot count. Later, Vt is pruned to a few directions with high shot count, concentrating the budget along the gradient direction ∇L.

U (θ1 ) · · · U (θd ) is applied to an initial state |ψ⟩4 and measured with an observable O, yielding output functions f (θ) = ⟨ψ| U † (θ)O U (θ) |ψ⟩. Unlike classical models, we cannot directly access intermediate states to apply the chain rule, so dedicated gradient extraction methods are required. The Hadamard test [10, 34, 35] is the most general approach, computing gradients via controlled unitaries and an ancillary qubit, but the controlled operations are expensive on near-term hardware. The parameter-shift rule [36] avoids this overhead entirely, computing gradi-

4 Typically |0⟩⊗n , but for QML applications one may encode data via |ψ⟩ = Uenc (x) |0⟩⊗n , possibly with data reuploading [33].

ents using the same circuits as function evaluation with no ancillae or controlled gates. The most general form handles arbitrary gate generators [37, 38]; for the common case of gates U (θj ) = eiθj Gj where Gj has two distinct eigenvalues (e.g. Pauli matrices), the rule simplifies to ∂f (θ) 1 = f θ + π2 ej − f θ − π2 ej ∂θj 2

(5)

requiring just two circuit evaluations per parameter. More generally, the gradient can be written as a linear combination of R shifted function evaluations, ∂f/∂θj = PR j j r=1 yr fr (θ), where R depends on the spectrum of the generator (e.g. R = 4 for beam-splitter gates [12, 39]). Crucially, parameter-shift rules yield exact derivatives

4 with non-infinitesimal shifts, distinguishing them from finite-difference approximations. Quantum backpropagation. Both the Hadamard test and parameter-shift rule compute each gradient component independently, giving an O(N ) lower bound on gradient evaluation for N -parameter models5 . Ref. [10] developed a quantum backpropagation algorithm achieving O(polylog(N )) gradient scaling, but the construction requires a fault-tolerant device with coherent quantum random-access memory and O(poly(n)) ancilla overhead, neither of which is available on current hardware. In a complementary direction, Refs. [11–13] identify model families with efficient gradient scaling by restricting the gate set to subgroups that admit classical back-substitution, at the cost of reduced representational capacity. C.

1 X om (θ), M m=1

E[oM ] = ⟨ψ| O |ψ⟩θ

(6)

which is an unbiased estimator of the expectation value. For bounded observables ∥O∥ ≤ 1 (e.g. a Pauli string, or a normalised Hamiltonian term), each single-shot outcome has variance at most one, so Var[oM ] ≤ 1/M , and 2 shot noise is fully controlled by M . This is the σshot /M term that will reappear in the forward gradient convergence analysis of Section IV B. Each shifted function evaluation in eq. (5) is estimated by such a sample mean, so the gradient estimator for the j-th component becomes e (θ)]j := [∇f

R X

yrj oj,M (θ) r

Approximate gradient estimators: SPSA and RCD

Two well-known alternatives avoid the explicit O(N ) per-step cost by estimating a single random gradient component or direction rather than the full gradient vector. SPSA. The SPSA algorithm [14, 42, 43] estimates the gradient via a single random perturbation with Rademacher vector u (ui ∈ {−1, +1}): eSPSA (θ) = g

(7)

r=1

where oj,M is the M -shot sample mean for the r-th r shifted circuit of parameter j, and unbiasedness follows by linearity. In practice one typically sets Mrj := M uniformly, though this can be relaxed [41]. By standard Hoeffding arguments, estimating eq. (6) to within precision ξ with probability 1 − δ requires M = O(ξ −2 log(1/δ)) shots. For the full N -parameter gradient eq. (7), each of N components requires R shifted circuits, giving a total

5 More

Estimators for quantum gradients

Unlike classical models where gradients can often be computed exactly, each circuit evaluation in the parameter-shift rule eq. (5) must be estimated from finite measurement samples. Following [40], let om (θ) denote the outcome of a single measurement trial, and define the M -shot sample mean oM (θ) :=

cost of O(RN ξ −2 log(N/δ)) shots. This O(N ) scaling is the fundamental bottleneck of parameter-shift gradient estimation.

precisely, O(N ) applies per parameter; the full gradient costs O(N 2 ) if each shifted circuit is a separate execution. See [10, 11] for precise statements.

f (θ + εu) − f (θ − εu) u 2ε

(8)

requiring only two circuit evaluations per step, independent of N . SPSA has been applied to VQE [44], quantum natural gradients [45], combinatorial optimisation [46], and quantum control with very few shots per evaluation [47]. However, the N -dependence removed from the per-step cost reappears in the estimator variance, which scales linearly in N [10]. The total cost to convergence is therefore comparable to the parameter-shift rule, but in practice SPSA performs well at the scales of variational quantum algorithms studied in the literature thus far [44, 48]. Random coordinate descent. RCD [15, 49, 50] takes a complementary approach: it selects a single parameter index j uniformly at random and computes the exact gradient in that component via the parameter-shift rule: eRCD (θ) = [∇f ]j ej , g

j ∼ U ({1, . . . , N })

(9)

As with SPSA, the per-step saving is offset by O(N ) slower convergence [15]. Both methods highlight a fundamental tension: removing the explicit N -dependence from the gradient rule shifts it to the variance or convergence rate.

Adaptive shot allocation: the xCANS family

A complementary approach to reducing measurement cost is to keep the O(N ) per-step structure of the parameter-shift rule but allocate shots adaptively across gradient components, spending more on noisy or informative directions and less on well-determined ones. The xCANS family of optimisers [23, 24] formalises this via the expected gain of a gradient step.

5 1.

Expected gain

e For gradient descent with an unbiased estimator ∇L on an L-Lipschitz loss, the improvement in the loss from one step to the next can be lower bounded [22]. We define the gain G as this lower bound: t e G := η∇L(θ t )⊤ ∇L(θ )−

Lη 2 e ∥∇L(θ t )∥2 2

(10)

satisfying L(θ t ) − L(θ t+1 ) ≥ G. Taking the expectation (over data subsampling classically, or measurement shots quantumly) and using unbiasedness, the expected gain becomes Lη 2 Lη 2 E [G] = η − ∥∇L∥2 − Tr(Σ(θ t )) (11) 2 2M where Σ(θ) is the covariance of the gradient estimator and M the number of measurement shots. The gain is positive when η < 2/L and the noise term Tr(Σ)/M is sufficiently small relative to ∥∇L∥2 .

Hamiltonian terms and parameter-shift terms, achieving unbiased gradient estimates from a single measurement shot per step. Refoqus [41] generalises these ideas to data-dependent QML losses (e.g. MSE), where shot budgets must also be distributed across data points. Closely related is the metric-aware shot-budgeting analysis of [52], which derives an optimal allocation between matrix and vector entries of the quantum natural gradient. A complementary line of work reduces the percomponent shot cost not by reallocating shots but by reusing measurement data: classical-shadow methods such as CoVaR [53] and adaptive informationally complete POVMs [54] can estimate many observables from a single measurement record, attacking the same O(N ) bottleneck from a different direction. All of these optimisers, however, assume per-parameter gradient computation via the parameter-shift rule. Table I (Section III) summarises the methods reviewed in this section together with the forward-gradient framework and the quiver optimiser developed in this work. III.

Shot allocation rules

iCANS [23] maximises eq. (11) per gradient component, yielding Mi =

ei 2Lη σ 2 − Lη gei2

(12)

where gei and σ ei are empirical estimates of the i-th gradient component and its variance. Components with large signal-to-noise ratio receive fewer shots; noisy components receive more. The global variant gCANS [24] couples the allocation across all parameters by replacing the per-component denominator with the full gradient norm: PN ei n=1 σ en 2Lη σ Mi = 2 2 − Lη ∥e g∥

(13)

In practice, both algorithms use bias-corrected exponential moving averages for gei and σ ei . The critical feature that makes these allocations effective is that measurement variances σ ei differ across parameters in the parameter-shift setting, since each shifted circuit perturbs a different gate.

The gradient methods reviewed above each pay an O(N ) cost in some form: the parameter-shift rule explicitly, in N separate component evaluations per step; SPSA implicitly, through an O(N ) amplification of the estimator variance from packing the full gradient into a single random perturbation; and RCD through an O(N ) slow-down in convergence from updating one coordinate at a time. In this section we propose forward gradients, an estimator class that unifies these methods as special cases of a single random-directional-derivative estimator with a tunable parameter V (the number of random directions sampled per step). The O(N ) trade-off persists, but is consolidated into a single parameter V through which one interpolates between these trade-off extremes, as the remainder of the paper develops. Table I maps the landscape: the top block places parameter-shift, SPSA, and RCD as fixed-(V, direction) instances of the estimator; the middle block lists existing shot-adaptive optimisers, all built on the parameter-shift estimator; and the bottom block summarises the quiver family introduced in this work, which adapts V (and M ) on top of the forward-gradient estimator. A.

FORWARD GRADIENTS

Forward gradients in machine learning

Extensions

Several works extend the xCANS framework by introducing additional sources of adaptivity. Random operator sampling [40, 51] distributes shots P across terms in a Hamiltonian decomposition O = k ck Ok , weighting each term by |ck |; when combined with iCANS, this yields the Rosalin optimiser [51]. The doubly-stochastic parameter-shift rule [40] goes further, sampling both

The forward mode of AD computes Jacobian-vector products J v = ∇v f for a given direction v. As discussed in Section II A, this is less efficient than backpropagation for the typical N ≫ M regime, since each forward pass yields only one directional derivative rather than the full gradient. Nevertheless, a productive line of work [16– 18, 56] has shown that a useful gradient estimator can be constructed from a small number of such directional

6 TABLE I. Unified taxonomy of gradient estimators and shot-adaptive optimisers. Gradient estimators and shotadaptive optimisers for variational quantum algorithms, framed as instances of the forward-gradient estimator eq. (15). Top block: gradient estimators, parameterised by (V, direction distribution), with the directional derivative computed via the parameter-shift rule (PSR) for basis directions and central finite difference (FD) for Rademacher directions. Middle block: existing shot-adaptive optimisers, all built on the parameter-shift estimator (V = N , basis, PSR). Bottom block: the quiver family introduced in this work, built on the forward-gradient estimator (tunable V , Rademacher, FD). Loss type: “Lin.” = linear in circuit expectation values; “Any” = general. ∗ SSD [55] computes each directional derivative exactly via an inner product circuit with ⌈log2 N ⌉ ancilla qubits and controlled gates, rather than by central FD. See Section III C for a full comparison. Method V Directions v Gradient estimators Parameter-shift [36] N basis ej SPSA [14] 1 Rademacher RCD [15] 1 N ej , j uniform SSD∗ [55] 1 µ=0, σ=1 Forward gradient (this work) tunable µ=0, σ=1 Shot-adaptive optimisers, parameter-shift base Op. sampling [51] N basis ej Double stoch. [40] N basis ej iCANS [23] N basis ej gCANS [24] N basis ej Rosalin [51] N basis ej Refoqus [41] N basis ej Shot-adaptive optimisers, forward-gradient base (this work) quiver fixed-M tunable µ=0, σ=1 quiver joint (V, M ) tunable µ=0, σ=1

derivatives, with multi-tangent aggregation [56] explicitly studying V > 1 via orthogonal projection of the random tangents: e F f (θ) := (∇v f ) v ∇ v

(14)

The full gradient estimator is obtained by averaging over V sampled directions V := {v ℓ }Vℓ=1 : e F f (θ) = 1 ∇ V

V X (∇vℓ f ) v ℓ

(15)

ℓ=1

When the components of each v ℓ are drawn independently with zero mean and unit variance, eq. (15) is an unbiased estimator of ∇f (θ) [16]. In the classical setting, the estimator’s variance scales as O(N ) [21, 57], so more directions are needed as the number of parameters grows; this dimension penalty can be partially circumvented by repeated sampling of each direction [58], the classical analogue of the per-direction shot count M in the quantum setting, and the multi-tangent extension that averages over V > 1 random directions with explicit orthogonal projection [56] is the closest classical antecedent of the V -interpolation unification of Section III B. We establish the corresponding bound for the quantum estimator (with the additional contribution from measurement shot noise) in Section IV, which to our knowledge is the first such result for the central-difference forward gradient estimator on parameterised quantum circuits. A separate line of work addresses variance reduction by exploiting the layered structure of classical neural networks: activity perturbation [21, 59], local losses with auxiliary per-layer networks [20, 21], and Jacobian ma-

Adapts

Loss

Derivation principle

− − − − −

Any Any Any Lin. Any

− − − − −

Mej per op Mej + ops Mej per param Mej per param (global) Mej + ops Mej + batch

Lin. Lin. Lin. Lin. Lin. Any

min-variance op weights random-coord. per shift gain per shot gain per shot gain per shot gain per shot

V V , Mv

Lin. Lin.

gain per shot min-cost at MSE target (CRB)

nipulation that reduces both bias and variance by exploiting low-dimensional gradient structure [60]. These mitigations do not transfer to quantum circuits, where the loss is computed by a single non-decomposable circuit evaluation [61]: classical forward-mode AD makes a single directional derivative O(1) at the cost of one function evaluation, whereas no known quantum subroutine recovers this scaling, so the quantum saving comes instead from the V ≪ N random-projection reduction we develop below. Two classical extensions are natural follow-ups but lie outside the scope of this work: the second-order forward-mode AD of [62] adds directional curvature information at the cost of one hyperdual-number forward pass per direction, and the randomsubspace exact-Generalized-Gauss–Newton construction of SOFO [63] would require a quantum analogue of the insubspace Hessian inverse. Quantum-side counterparts to these second-order constructions already exist (quantum natural gradient [64], simultaneous-perturbation natural gradient (QN-SPSA) [45], and parameter-shift Hessian estimators [65]), and a head-to-head comparison of a directional-Hessian forward-gradient construction against these baselines is a natural follow-up that we do not pursue here.

Existing gradient estimators as special cases

The forward gradient framework eq. (15) recovers SPSA, RCD, and the parameter-shift rule as special cases under specific choices of V and the direction distribution. SPSA eq. (8) corresponds to V = 1 with a Rademacher direction and central finite-difference directional deriva-

7 tive. RCD eq. (9) corresponds to V = 1 with v = ej drawn uniformly from the N basis vectors, rescaled by a factor of N to give an unbiased estimate, and the parameter-shift rule for the directional derivative. The full parameter-shift rule eq. (5) is then recovered in the deterministic limit where all N basis vectors are enumerated in a single step rather than sampled; this is the V → N deterministic limit of the framework, in which direction randomness is replaced by full enumeration. As a consistency check, the embedding reproduces the known behaviour of its limiting cases: SPSA’s O(N ) variance amplification and RCD’s O(N ) slow-down relative to parameter-shift both emerge directly from the second-moment expansion developed below. The new ingredient the common framing supplies is an explicit V that interpolates between the cheapest single-direction estimator (V = 1, highest variance) and the exact full gradient (V = N , zero direction variance, highest cost), and it is this lever that the rest of the paper explores, both as a fixed hyperparameter in Section V and as an adaptive quantity in Section VII. C.

Quantum directional derivatives

Applying eq. (15) to quantum circuits requires an efficient method for computing directional derivatives ∇v f . Computing the exact directional derivative via the parameter-shift rule requires evaluating the full gradient first (at cost O(RN )) and then projecting onto v, which negates any saving from the forward gradient framework. However, Ref. [66] observed that for gates of the form eiθG , loss functions expressible as linear measurements of observables are trigonometric polynomials in the parame√ ters, with maximum angular frequency bounded by N . This band-limitation means that finite-difference approximations with large step sizes ε can achieve high accuracy using only Q circuit evaluations per direction, with Q ≪ RN . In particular, defining the finite-difference estimator: M

shots. The total measurement cost is 2V M shots per step, far below the parameter-shift cost of 2RN M when V ≪ N. The choice of finite-difference step ε governs a bias–variance trade-off in the central-difference estimator eq. (16): shot noise is amplified by 1/ε2 , while the leading bias scales as O(ε2 ). We find ε = 0.1 to be the best-performing value on a representative VQE example, and use it for all subsequent experiments. App. B derives the corresponding closed-form optimum and shows that it agrees with this empirical choice. A recent complementary approach [55] (stochastic shadows descent) computes ∇f · v as a single blockencoded subroutine using N controlled copies of the ansatz and log2 N ancillas, removing the finite-difference bias. The practical advantage of forward gradients comes from a structural classical fact: forward-mode automatic differentiation computes a directional derivative in one pass at the cost of a single function evaluation, so reconstructing the gradient from V ≪ N random projections gives an N/V saving. No known quantum subroutine recovers this scaling; every exact directionalderivative method on a PQC pays O(N ), either as N parameter-shift queries or as the depth-N block encoding of Ref. [55]. The circuit of Ref. [55] produces one exact DD per query, but a V -projection reconstruction has total gate count V N ·depth(U ), a factor V worse than plain parameter-shift, with equality only at V = 1. The bias removal therefore moves the N cost from query count to query depth without saving anything over PS, and pushes the construction past the NISQ regime. The quantum approximate forward gradient estimator proposed here accepts an O(ε2 ) bias on the directional derivative instead, recovering the classical V ≪ N saving with O(depth(U )) per query; the bias is analytically controlled via the optimal step size ε⋆ that balances finite-difference bias against shot noise (App. B).

IV. CONVERGENCE OF QUANTUM FORWARD GRADIENT DESCENT

e εv f := f (θ + εv) − f (θ − εv) ∇ (16) 2ε PM where f M (θ ′ ) := 1/M m=1 om (θ ′ ) is the M -shot estimate of the circuit output at parameters θ ′ . Ref. [66] validated this approach numerically, showing that Q = 2 (central difference) with step sizes ε ∈ {0.05, 0.1, 0.2, 0.3} suffices for accurate nuclear gradient estimation in VQE. Higher-order stencils (Q = 4, 6) can improve accuracy at the cost of additional circuit evaluations; in this work we use the central difference throughout. Combining eq. (16) with eq. (15), the quantum forward gradient estimator becomes V X e f= 1 e ε ℓ f ) vℓ ∇ (∇ v V F

(17)

ℓ=1

requiring 2V circuit evaluations per gradient step (two per direction for the central difference), each with M

The reduced per-step cost of the forward gradient estimator eq. (17) naturally raises the question of whether this saving translates to a reduction in total measurement cost to reach a fixed accuracy, or whether the increased variance of the estimator compensates. We analyse this in stages, first assuming access to exact directional derivatives, then incorporating the finite-difference approximation and quantum measurement noise.

Convergence with exact directional derivatives

Consider the idealised setting where each directional derivative ∇vℓ f = v ℓ ·∇f is computed exactly (no finitedifference error, no shot noise). The forward gradient estimator eq. (15) is then an unbiased estimator of ∇f . Its

8 second moment is characterised by the following lemma (proved in App. C). Lemma 1 (Second-moment expansion). Let {v ℓ }Vℓ=1 be drawn i.i.d. with E[viℓ ] = 0, E[viℓ viℓ′ ] = δii′ , and let κ := E[(viℓ )4 ] be the fourth moment of one component (κ = 1 for Rademacher, κ = 3 for standard Gaussian). Then i h F e f ∥2 = N + V + κ − 2 ∥∇f ∥2 . Ev ∥∇ V

(18)

Throughout the remainder of this paper we use Rademacher directions (κ = 1) unless stated otherwise, so the prefactor reduces to (N + V − 1)/V . The convergence rates for forward gradient descent follow from Lemma 1 and are stated as a proposition (proved in App. D). Proposition 1 (Forward gradient convergence). Let f ∈ C 1 (RN ) be L-smooth and satisfy the PL inequality with constant µ > 0, and let directions be drawn from the Rademacher distribution (κ = 1). (i) Convergence of forward gradient descent. With learning rate η = V /(L(N + V − 1)), forward gradient descent satisfies E[f (θ (T ) )] − f ∗ ≤ 1−

µV L(N +V −1)

(f (θ (0) ) − f ∗ )

(19)

(ii) Convergence of quantum forward gradient descent. 2 With per-direction shot noise variance σshot and learning rate η ∈ [0, 1/(2µ)], the expected suboptimality satisfies T

E[f (θ (T ) )] − f ∗ ≤ (1 − 2µη) (f0 − f ∗ ) +

2 Lη (N +V −1) σshot · 4µ VM

(20)

Proof sketch. The forward gradient estimator is unbiased (Proposition 3). Lemma 1 bounds its second moment by β 2 ∥∇f ∥2 with β 2 = (N + V − 1)/V (using κ = 1). Substituting into the standard PL-SGD bound [40, 67] yields (i). For (ii), adding the constant shot-noise contri2 bution σshot /M to the second moment gives a boundedvariance estimator; the bounded-variance form of the same SGD bound then yields the asymptotic floor directly. The formal derivation is in App. D. To reach f (θ (T ) ) − f ∗ ≤ δ under (i), we therefore require N L 1 T =O · log (21) V µ δ steps, where the last equality holds for V ≪ N . Since each step costs O(V ) circuit evaluations, the total circuit cost is LN 1 T ×V =O log (22) µ δ

which is independent of V . This is the no-free-lunch result for forward gradients: the V -fold saving in per-step cost is exactly compensated by a V -fold increase in the number of steps required, and the total circuit evaluation budget is the same as for the full gradient (V = N ). In this exact setting, the choice of V affects only the granularity of progress (many cheap steps vs few expensive ones), not the total cost.

Convergence under quantum measurement noise

In practice, each directional derivative is estimated from M measurement shots, contributing a per-direction 2 2 variance σshot /M to the estimator, where σshot := e Em [Varv [∇v fm ]]. Proposition 1(ii) captures this: the first term of eq. (20) is the optimisation error, which decreases geometrically and is independent of how the shot budget is allocated, while the second term is an asymptotic error floor (a residual suboptimality set by the shot noise that cannot be averaged away within a single gradient step). The floor is proportional to (N + V + κ − 2)/(V M ). To compare methods at equal cost, we impose a fixed total per-step shot budget B and allocate it as M = B/(2V ) shots per direction, so that each of the V directions contributes two circuit evaluations each estimated with M shots. Substituting into the floor gives 2 Lη 2(N +V +κ−2) σshot · 4µ B

(23)

For V ≪ N the +V term in the numerator is negligible, 2 so the floor reduces to approximately Lη/4µ·2N σshot/B with a correction of order V /N ; crucially, this expression is independent of V , since increasing V reduces the directionaveraging variance by 1/V while the corresponding decrease in M = B/(2V ) raises the per-direction shot noise by the same factor, and the two effects cancel exactly. This is the no-free-lunch result: at fixed per-step shot budget, no choice of V lowers the asymptotic shot-noise residual. The experimental savings reported in Section V and Section VI therefore cannot come from reducing this residual; they arise before it is reached, where the nonconvex landscape still provides strong gradient signal and a larger V finds better descent directions early in training.

FORWARD GRADIENTS AT SCALE

To validate these practical savings, we train Hammingweight preserving orthogonal circuits [12, 25, 26] on ECG5000 time-series classification (n up to 60, N up to 1770 parameters) and MNIST image classification (n up to 50, N up to 1225), comparing forward gradient estimators against the parameter-shift rule, SPSA, and RCD.

(a) QAOA MaxCut, n=16

(b) VQE TFIM, n=10

(d) MNIST, n=50

100 Exact: 100%

Exact: 76%

Exact: 64%

Test accuracy (%)

100

|E − E0 |

Approx. ratio (%)

85 80

50 40 40 30

30 20

N =99

102

105

10−1

Total shots

N =160

102

N =1770

105

107

Total shots Forward (best V )

109

Total shots PS

SPSA

N =1225

109

1011

Total shots

RCD

FIG. 2. Forward gradients match or exceed all baselines across four problem domains: quantum optimisation, quantum simulation, time-series classification, and image classification. Each panel compares the best forward gradient configuration against the parameter-shift rule (PS), SPSA (V = 1) and RCD on the total measurement shots axis. All stochastic methods are run at a matched per-step shot budget B = 1000; PS uses M = 10 shots per parameter and its per-step cost therefore grows as 2N M . (a) QAOA MaxCut, n = 16, depth 3, N = 99, best forward V = 10. (b) VQE ground-state estimation (transverse-field Ising model, n = 10, d = 8), N = 160, best forward V = 10. (c) ECG5000 time-series classification, n = 60, N = 1770 parameters, best forward V = 25. (d) MNIST image classification, n = 50, N = 1225, best forward V = 10. Forward gradients reach near-optimal performance on all four domains, at a fraction of the PS shot cost on the three panels (a, b, c) for which a PS baseline is shown, while SPSA and RCD degrade at large N . The MNIST panel (d) reports forward, SPSA, and RCD curves; PS is omitted from this panel because its per-step cost 2N M is computationally prohibitive at these system sizes, and the PS-vs-forward comparison on MNIST is reported separately at n ∈ {10, 20, 40} on a 10×-subsampled training set in Fig. 5. Shaded bands: one standard deviation over three random seeds (all four panels). The QAOA panel is shown on a random instance of a weighted Erdős–Rényi graph. All methods use the Adam optimiser at identical hyperparameters; only the gradient estimator differs.

All stochastic methods (forward, SPSA, RCD) are run at a fixed per-step shot budget of B = 1000 shots, so that comparisons between them are at matched measurement cost; the parameter-shift rule is run with M = 10 shots per parameter, so its per-step cost grows as 2N M and is not matched. Results are mean accuracy over three random seeds, reported against total shots used. Fig. 2 summarises the results, confirming empirically the cartoon behaviour of Fig. 1(c). Full hyperparameters are in App. L.

Fig. 4 compares forward methods against the parameter-shift rule. Despite using a fixed per-step budget (independent of N ), forward gradients reach comparable or superior accuracy to PS at all system sizes. The accuracy gap between forward and exact gradients remains small (1–2 percentage points) up to n = 60, while the total shot cost is a fraction of that required by PS. The best-performing forward configuration uses V = 10 at most system sizes, with V = 25 preferred at n = 60.

ECG5000

Fig. 3 compares forward gradient methods against SPSA and RCD at a fixed shot budget per step. Forward gradients with V ≥ 10 outperform single-direction methods (V = 1), with the gap widening at larger N . At n = 60 (N = 1770), SPSA collapses to 63.8% while forward methods with V = 25 match the exact gradient baseline at 76.0%.

MNIST. On MNIST image classification the same ordering holds (Fig. 5, App. L): across system sizes n ∈ {10, 20, 40, 50} (N up to 1225), forward gradients with V = 10 track the exact-gradient baseline at the matched B = 1000-shot per-step budget, while SPSA and RCD degrade with n and the parameter-shift-versus-forward shot ratio scales linearly in N .

n = 10, N = 45

n = 20, N = 190

Test accuracy (%)

n = 10, N = 45

n = 20, N = 190

0.0

2.5

5.0

n = 30, N = 435

20 0

Total shots ×109

Exact gradient Fwd V =2, M =250 Fwd V =10, M =50

Fwd V =25, M =20 Fwd V =50, M =10

Total shots ×109

SPSA M =500 RCD M =500

FIG. 3. Forward gradients consistently outperform SPSA and RCD at matched shot budget. ECG5000: forward gradient methods (V > 1) vs SPSA (V = 1) and RCD. Each panel shows a different system size n. Forward gradients with V = 10 consistently outperform singledirection methods.

0.5

1.0

1.5 ×1010

n = 60, N = 1770

Test accuracy (%)

n = 50, N = 1225

0.0 80

7.5 ×109

n = 40, N = 780

Test accuracy (%)

×109

n = 50, N = 1225

20 0

Total shots ×1010

Exact gradient PS M =10 PS M =50

n = 40, N = 780

3 ×1010

n = 60, N = 1770

Fwd V =2, M =250 Fwd V =10, M =50

Total shots ×1010

Fwd V =25, M =20 Fwd V =50, M =10

FIG. 4. Forward gradients reach parameter-shift accuracy at a fixed, N -independent shot cost. ECG5000: forward gradient methods vs parameter-shift rule (PS) with M shots per parameter. Forward methods use a fixed per-step budget independent of N , yet reach comparable accuracy. We do not include the larger-shot M = 50 baseline at n = 50, 60 due to the extensive runtimes required (2N M shots per step for N ∈ {1225, 1770}).

Ground-state estimation and combinatorial optimisation

Beyond classification, we validate forward gradients on two further problem classes for which the optimisation landscape and the measurement noise are qualitatively different. The ground-state estimation benchmark is a transverse-field Ising model (TFIM) on n = 10 qubits with d = 8 layers (N = 160 parameters), trained to minimise the energy expectation value of a hardware-efficient ansatz. The combinatorial optimisation benchmark is a depth-3 QAOA for MaxCut on 16-vertex weighted Erdős–Rényi graphs with N = 99 parameters, trained to maximise the approximation ratio. Both are shown in panels (c) and (d) of Fig. 2, under the same matched per-step budget B = 1000 convention as the classification

benchmarks. On VQE, forward gradients with V = 10 reach a mean best-over-training energy error |Ebest − Eexact | = 0.115 ± 0.010 at a 5 × 106 -shot budget (mean and standard deviation over three random seeds; energy readout from a 104 -shots/term noisy oracle excluded from the training budget), compared with 0.203 ± 0.021 for Adamoptimised parameter-shift at M = 10 shots per parameter over three seeds and 0.215 for SPSA at the same budget (single-seed reference configuration). On QAOA, forward gradients with V = 10 reach an approximation ratio of 0.988 at a 5 × 106 -shot budget, within ∼ 0.01 of the exact-gradient baseline and well ahead of SPSA at the same budget. At this system size

Test accuracy (%)

n = 20, N = 190

10 107

Test accuracy (%)

109

1010

n = 40, N = 780

108 50 40

10 107

Test accuracy (%)

108

2)/V in the gain expression eq. (18) decreases with V , suggesting that fewer directions yield more precise gradient estimates per shot. However, the number of steps to convergence grows as N/V (eq. (21)), favouring more directions. The optimal balance between these competing effects shifts during training: early on, when the gradient is large and the loss landscape has many viable descent directions, a large V provides broad coverage at low cost per direction. As training progresses and the gradient shrinks, precision matters more than coverage, and concentrating shots on fewer directions becomes preferable.

n = 20, N = 190

108

109

1010

n = 50, N = 1225

108 50

1010

n = 40, N = 780

This suggests a simple scheduling strategy: decay V from a large initial value Vmax to a small final value Vmin over the course of training, while keeping the total shot budget per step B = 2V M fixed, so that as V decreases, the shots-per-direction M = B/(2V ) increases automatically. We consider three monotone schedule shapes, parameterised by training step t ∈ [0, T ]:

1010

n = 50, N = 1225

10 107

108

109

1010

Total shots Exact gradient Fwd V =2, M =250

108

1010

1012

Total shots

Fwd V =10, M =50 Fwd V =25, M =20

Fwd V =50, M =10 SPSA M =500

linear:

Vt = Vmax − (Vmax − Vmin ) Tt ,

cosine:

Vt = Vmin + 21 (Vmax − Vmin )

RCD M =500 PS M =10

exponential: FIG. 5. MNIST image classification reproduces the ECG5000 ordering across system sizes. MNIST 10class classification reproduces the ECG5000 ordering on image data. Rows: system sizes n. Left column: forward gradient family vs SPSA and RCD; right column: forward gradients vs PS (M = 10 shots per parameter), trained on a 10×-subsampled training set so the larger PS comparison is tractable. Shaded bands: ±1σ across three random seeds. Same hyperparameters and matched B = 1000-shot per-step budget as Fig. 2.

(N = 99), the parameter-shift rule remains a viable comparison point, but its per-step cost 2N M grows rapidly with system size and becomes the dominant factor at the classification scales of Section V A. The four panels of Fig. 2 span ground-state estimation, combinatorial optimisation, and supervised learning, and the matchedbudget forward gradient estimator is competitive with or better than the best available baseline on every one.

VI.

ADAPTIVE DIRECTION SCHEDULING A.

(24) 1 + cos πt T

Motivation

The results of the previous section show that forward gradients with a fixed number of directions V can train quantum circuits at scale with practical shot savings. A question this raises is whether V itself can be adapted during training. The convergence analysis of Section IV provides the motivation. The noise amplification factor (N + V + κ −

Vt = Vmin

Vmax Vmin

(T −t)/T

, (25) (26)

In each case Vt is rounded to the nearest integer and clamped to [Vmin , Vmax ]. The per-step shot cost B = 2Vt Mt is unchanged; only the allocation between the number of directions and shots per direction varies.

Results

We illustrate this across four benchmarks (Fig. 6). In each case, a decreasing V -schedule outperforms both fixed-V endpoints at matched shot budget: exponential decay on VQE TFIM (n = 10, d = 8, N = 160), linear decay on QAOA MaxCut (n = 16, N = 99), and linear or cosine decay on ECG5000 and MNIST (n = 50, N = 1225). The primary driver of the improvement is the direction of decay rather than the precise functional form, though the optimal shape and range of Vmax → Vmin are problem-dependent and the same schedule does not achieve the same margin on every benchmark. These results are presented as a heuristic illustration rather than a systematic study: the observed improvement is problemand configuration-dependent, and we make no claim that a decaying V -schedule will consistently outperform fixedV endpoints in general. Rather, the observation that the training signal itself drives a natural preference for fewer directions as training progresses motivates the quiver optimiser of Section VII, which determines V automatically from the measurement data.

12 (a) QAOA n=16 (N =99)

(b) ECG5000 (n=50, N =1225)

(d) VQE (TFIM n=10, N =160)

85 80 75 70

Fixed V =2 Fixed V =30 Linear V : 30→2

104

105

Total shots

106

40 Fixed V =50 Fixed V =2 Linear V : 50→2

30 Fixed V =150 Fixed V =5 Cosine V : 150→5

20 0

Total shots ×109

Total shots×1010

Energy error (E − E0 )

Test accuracy (%)

Approx. ratio (%)

100

10−1

Fixed V =20 Fixed V =2 Exp. V : 20→2

104

106

Total shots

FIG. 6. Decreasing V -schedules outperform fixed-V endpoints across all problem domains. Each panel shows a scheduled V (blue) against the two fixed-V endpoints (grey solid = Vmax , grey dashed = Vmin ). The scheduled curve reaches higher performance at matched shot budgets than either endpoint on all four benchmarks: (a) QAOA MaxCut (n = 16, N = 99), linear V : 30 → 2; (b) ECG5000 (n = 50, N = 1225), linear V : 50 → 2; (c) MNIST (n = 50, N = 1225), cosine V : 150 → 5; (d) VQE TFIM (n = 10, N = 160), exponential V : 20 → 2. The direction of decay is the primary driver of the improvement; the optimal functional form and V -range are problem-dependent.

VII.

QUANTUM-NATIVE OPTIMISATION FOR FORWARD GRADIENTS

The forward gradient results of Section V and Section VI treat the direction vectors as purely random: drawn from an isotropic Rademacher distribution, with V either fixed or scheduled by a classical rule that ignores the measurement outcomes themselves. A natural question, in the spirit of the CABS [22] and iCANS/gCANS [23, 24] families of shot-adaptive optimisers for the parameter-shift rule, is whether the quantum measurement statistics generated during training can themselves be used to guide the estimator. In this section we answer that question by constructing quiver (Quantum Iterative V-adaptive Estimator Refinement), an adaptive optimiser for forward gradients. The derivation that follows assumes a loss linear in circuit expectation values, as is the case for the VQE and QAOA benchmarks of Section V; an extension to non-linear losses (such as the MSE and cross-entropy classification losses), in the spirit of Refoqus [41], is left to future work. The construction proceeds in five steps. In Section VII A we derive the expected gain of a forward gradient step, extending the per-component gain decomposition of iCANS to the directional-derivative setting. The resulting expression exposes two adaptive levers: the number of shots M allocated to each direction, and the number of directions V itself. In Section VII B we show that the per-direction shot lever (the natural forward-gradient analogue of iCANS) is neutralised by a concentration-of-measure phenomenon that is structural to random-direction estimators and absent from the parameter-shift setting. In Section VII C we turn to the remaining lever and find that the natural “fixed-M , optimise V for gain-per-shot” construction yields a closedform rule with a signal-independent floor: it picks a rea-

sonable constant V automatically but does not adapt within a run. We then decouple V and M into two independent targets, a per-direction signal budget and an absolute estimator-variance budget, and derive a joint closed-form for (Vt⋆ , Mt⋆ ) that genuinely adapts. Section VII D re-derives this rule from a single optimisation principle (minimum shot cost at a target estimator accuracy with a per-direction shot floor) and shows that the resulting allocation (Theorem 1) uses Rademacher directions that minimise variance among isotropic distributions (Proposition 2) and matches the Cramér–Rao lower bound on gradient recovery from a shot-noise oracle up to a vanishing constant (Corollary 1). Finally Section VII E validates the joint rule on VQE TFIM and compares it head-to-head against fixed parameter-shift, iCANS, gCANS, and the fixed-M variant at matched total shot budget.

Gain expression for forward gradients

Our starting point is a per-direction decomposition of the expected gain eq. (10), which makes the contribution of each random direction explicit as a competition between a signal term and a noise term. This is the forward-gradient analogue of the per-component gain decomposition used by CABS [22] and iCANS [23] for classical and parameter-shift gradients respectively.

Lemma 2 (Per-direction gain decomposition). Substituting Lemma 1 into eq. (10) yields:

13 (i) Em Ev [G F ] decomposes as

ℓ=1

Mℓ∗ =

h i e vℓ Lm 2 Varm ∇ C∥∇L∥2 − (∇vℓ L)2

optimum

Mℓ∗

, 2V C := Lη(N +V +κ−2) .

(28)

Each direction ℓ contributes a signal term (∇vℓ L)2 and e vℓ Lm ] to the gain penalty, ama noise term 1/M Varm [∇ plified by (N +V +κ−2)/V . In eq. (28), directions better aligned with the gradient receive more shots (smaller denominator); in practice the unknown true quantities are replaced by empirical estimates (see App. F for positivity conditions). The allocation eq. (28) is the forwardgradient analogue of iCANS: shots are distributed across directions in proportion to the ratio of per-direction variance to a signal-dependent denominator. The next subsection shows why this per-direction allocation is structurally inert for forward gradients, at which point we return to the remaining lever (the number of directions V ) in Section VII C. B.

Why iCANS-style allocation fails for forward gradients

The allocation eq. (28) has the same structure as the iCANS allocation [23]: shots are distributed according to the ratio of per-direction measurement variance (numerator) to a signal-dependent denominator. For iCANS, both quantities vary across parameters, because each canonical direction ej shifts a single gate while all others remain fixed, producing quantum states that depend on which parameter is being differentiated. For forward gradients, each direction v ℓ perturbs all N parameters simultaneously with i.i.d. components, and the measurement variance of the resulting directional derivative is, under smoothness assumptions on the parameterised expectation value, a Lipschitz function of v ℓ . Classical concentration-of-measure results for Lipschitz functions of independent bounded random variables then yield the following heuristic, which we treat as a working assumption and verify empirically below. Assumption 1 (Noise concentration across random directions). Let v ℓ have i.i.d. components (Rademacher or

QAOA MaxCut n=8 n = 10

1.15

n=8 n = 10

1.10 2 / 2 ,

V F 1 X η∥∇L∥2 Em Ev G ≈ V ℓ=1 h i Lη 2 N + V + κ − 2 1 2 e − (∇vℓ L) + Varm ∇vℓ Lm 2 V M V 1 X =: γvℓ . (27) V

(ii) The gain-per-shot arg maxMℓ [γvℓ /Mℓ ] is

VQE TFIM

1.20

1.05 1.00 0.95 0.90 0.85 0.80

102

Number of parameters N

103

102

Number of parameters N

103

FIG. 7. Directional-derivative variance concentrates uniformly across random directions. For each configuration (Vprobe = 32 Rademacher directions; ε = 0.1, M = 50, 2 R = 5000), the 32 individual per-direction variances σ∇,ℓ are 2 plotted normalised by their configuration mean σ̄∇ . The solid line at y = 1 is the configuration mean; all points lie within ±15% of unity, and the spread is typically below ±10%. Left: VQE TFIM (n = 8, 10). Right: QAOA MaxCut (n = 8, 10). Each series sweeps N up to ∼ 103 by increasing depth d. The tight clustering around unity, across all N and both problem types, confirms Assumption 1.

2 e vℓ Lm ] denote the perGaussian) and let σ∇,ℓ := Varm [∇ 2 direction measurement variance. If σ∇,ℓ is Lipschitz in ℓ v with Lipschitz constant independent of N , then 2 2 σ∇,ℓ ≈ σ̄∇

for all ℓ,

(29)

2 2 ]. The Lipschitz condition ensures := Ev [σ∇,ℓ with σ̄∇ 2 2 with fluctuations O(Lσ ) inσ∇,ℓ concentrates around σ̄∇ dependent of N . We prove this for local Hamiltonians and bounded-depth ansätze in App. H and verify it empirically below.

Because each v ℓ averages over N i.i.d. perturbations, the per-gate variance structure that iCANS exploits is washed out. The numerator of eq. (28) becomes approx2 imately σ̄∇ for all ℓ, leaving only the denominator as a source of variation. But the denominator variation comes from the signal gℓ2 = (∇vℓ L)2 , which for random directions is a random projection of the gradient rather than a structured per-parameter quantity. The allocation re2 duces to Mℓ∗ ≈ σ̄∇ /(const − gℓ2 ), which provides only weak, stochastic variation compared to the structured numerator-driven allocation of iCANS. We verify this empirically in Fig. 7. For each configuration we draw Vprobe = 32 independent random directions e vℓ Lm ] for at the initial parameters and estimate Varm [∇ each by repeating the noisy central-difference estimator R = 5000 times at M = 50 shots per evaluation, with ε = 0.1. Series cover VQE TFIM and QAOA MaxCut at n ∈ {8, 10} each, with layer depth d increasing along each series to sweep N up to ∼ 103 . For each configuration 2 we plot the 32 individual values σ∇,ℓ normalised by their 2 mean σ̄∇ ; the points cluster tightly around unity (within ±10% typically) for all N and all four series, confirming 2 2 σ∇,ℓ ≈ σ̄∇ regardless of problem type or qubit count. Substituting eq. (29) into the per-direction gain makes

14 the structure explicit: Lη 2 (N +V +κ−2) . 2V (30) The noise penalty (third term) is identical for all ℓ, while the signal cost (second term) varies through gℓ2 . Since the only adaptive lever in the allocation (the numerator) has been removed by concentration, a different strategy is needed: rather than allocating shots non-uniformly across directions, we adapt the number of directions itself. This is the remaining degree of freedom in the forward gradient estimator, and the one that the V scheduling heuristic of Section VI already exploits empirically. In the next subsection we derive a principled rule for this adaptation.

γvℓ = η∥∇L∥2 − A gℓ2 −

2 A σ̄∇ , Mℓ

A :=

The quiver optimiser

The per-direction allocation of eq. (28) having been neutralised by the noise concentration, the remaining adaptive lever is the number of directions V itself, and optionally the scalar per-direction shot count M . Each epoch, the optimiser maintains two exponential moving averages at no extra measurement cost: gbt2 of the squared e F L∥2 and σ forward gradient norm ∥∇ bt2 of the mean per2 direction measurement variance σ̄∇ . We use these EMAs to update (Vt , Mt ) in two stages: a gain-per-shot closedform that yields a reasonable but static allocation, and a two-target refinement that recovers genuine within-run adaptation. Gain-per-shot optimum at fixed M . With M held fixed, the gain expression eq. (30) is a function of V alone. Lemma 3 (Fixed-M optimal V ). Taking the expectation of γvℓ over isotropic directions under Assumption 1 and maximising E[G F ]/(2V M ) over V , the first-order condition yields V⋆ =

and does not move thereafter, because at fixed M the optimum eq. (31) has a high-SNR floor

2 2Lη (M ∥∇L∥2 + σ̄∇ )(N + κ − 2) , 2 2 M ∥∇L∥ (2 − Lη) − Lη σ̄∇

(31)

clamped to [Vmin , Vmax ] = [2, N ], with V ⋆ = Vmax when the denominator is non-positive (proved in App. E 4). V ⋆ scales with N + κ − 2, reflecting the O(N ) variance amplification of forward gradients; both inputs are quantities already returned by the forward gradient step; and M is the single hyperparameter, with per-step cost Bt = 2Vt M varying naturally with Vt . The fixed-M optimum does not adapt. On VQE TFIM, per-epoch logging of Vt from eq. (31) shows that Vt converges to a fixed value on the first update step

⋆ Vfloor =

2Lη(N + κ − 2) 2 − Lη

(32)

which depends only on (L, η, N ), not on gbt2 or σ bt2 . The 2 2 formula leaves this floor only when σ b ≳ M gb , a condition that never holds during VQE training at the conservative learning rate η = 3 × 10−3 required for N ≥ 80 convergence: across 500 logged training epochs on VQE TFIM n = 8, d = 10 the signal-to-noise ratio M gb2 /b σ2 decreases from ≈ 4300 to ≈ 1400, an order-of-magnitude drop but still well above the threshold for adaptation. At Lη = 0.045, eq. (32) evaluates to V ≈ 15 at N = 320, which is precisely the steady-state value observed in the run data. The fixed-M rule therefore sets a constant V calibrated to (L, η, N ) with no additional tuning, and, as we show in the large-N panel of Fig. 8, this alone is enough to outperform iCANS/gCANS at the same shot budget when Lη ≪ 1. It does not, however, respond to training dynamics: within a run, Vt is flat. The joint (V, M ) rule. When the gradient signal dominates the measurement noise (b σ 2 ≪ M gb2 ), holding M fixed collapses the gain-per-shot objective to a function of (L, η, N ) alone, so V never moves in response to training dynamics. Allowing M to vary jointly with V restores this dependence: the pair (Vt⋆ , Mt⋆ ) is then the solution of a minimum-cost allocation problem whose inputs are the current signal and noise estimates, derived in Section VII D (Theorem 1), Vt⋆ =

(N − 1 + α) gbt2 , τ2

Mt⋆ =

Nσ bt2 α gbt2

(33)

with two hyperparameters: a dimensionless ratio α > 0 that controls the per-direction shot count and a target absolute variance τ 2 > 0 for the reconstructed gradient estimate. In this form the qualitative behaviour of the rule is immediate: M ⋆ grows as training reduces the signal (since gbt2 decreases), while V ⋆ shrinks with gbt2 . Noise dependence is preserved in M ⋆ through σ bt2 , and ⋆ 2 V tracks gbt directly. Both are computed once per epoch from the EMAs already maintained by the optimiser. We take eq. (33) as the defining update of quiver; the fixedM rule eq. (31) is recovered as a special case when M is held constant, and as we show in Section VII D, both rules are instances of the same minimum-cost allocation principle with different constraint choices. The closed form Vt⋆ ∝ (N − 1 + α)b gt2 /τ 2 inherits its N -dependence from both the explicit N prefactor and the gradient-norm scaling of the deployed problem; substituting the empirically observed gb2 ∝ N −1/2 scaling for our trained PQCs √ recovers the empirical V ⋆ ∝ N trend reported in our V -sweep experiments without further tuning. Because gbt2 and σ bt2 are noisy early in training, unclamped updates to Vt and Mt can swing by an order

15 of magnitude in a single epoch; we therefore apply a multiplicative clamp on each update and hold (Vt , Mt ) fixed during a short initial warmup while the estimates stabilise. Neither modification changes the fixed point of the rule. The full update is in Algorithm 1.

a constant that vanishes in N (Corollary 1). Parametershift saturates the same CRB on the same oracle, so both estimators are CRB-optimal and the practical advantage of Section V is a constant-factor gap within that class, not an information-theoretic separation.

Algorithm 1 quiver joint (V, M ) optimiser

bV for Rademacher is variance-optimal. Write g the forward gradient estimator eq. (15) with exact directional derivatives and g := ∇L(θ). The squared error E∥b g V − g∥2 depends on the choice of direction distribution; the following proposition identifies which distribution minimises it at fixed V .

Input: learning rate η, number of parameters N , targets (α, τ 2 ), warmup Tw , EMA decay β, rate-limit bounds (r↓ , r↑ ) = (0.7, 1.5), clamps [Vmin , Vmax ], [Mmin , Mmax ], initial (V0 , M0 ), finite-difference step ε 1: Initialise θ 0 , Adam state ϕ0 ; V ← V0 , M ← M0 2: g b2 ← None; σ b2 ← None 3: for t = 0, 1, 2, . . . do 4: Sample V random directions {v (ℓ) }Vℓ=1 (Rademacher) 5: Measure directional derivatives dℓ ← LeM (θ t + εv (ℓ) ) − LeM (θ t − εv (ℓ) ) /(2ε) e ← 1 PV dℓ v (ℓ) ▷ reconstructed gradient 6: ∇L ℓ=1 V 2 2 2 e 7: gt ← ∥∇L∥ ; σt ← Varℓ (dℓ ) 8: if gb2 = None then 9: gb2 ← gt2 ; σ b2 ← σt2 10: else 11: gb2 ← β gb2 + (1 − β) gt2 12: σ b2 ← β σ b2 + (1 − β) σt2 13: end if e η) 14: (θ t+1 , ϕt+1 ) ← Adam(θ t , ϕt , ∇L; 15: if t ≥ Tw then ▷ update (V, M ) after warmup 16: M⋆ ← Nσ b2 /(α gb2 ) 17: V ⋆ ← (N − 1 + α) gb2 /τ 2 18: V ← clip V ⋆ , max(r↓ V, Vmin ), min(r↑ V, Vmax ) ⋆ 19: M ← clip M , max(r↓ M, Mmin ), min(r↑ M, Mmax ) 20: end if 21: end for

Optimal allocation

The joint rule eq. (33) was stated in the previous subsection as the defining update of quiver, but not derived. We now derive it from a single optimisation principle: among all (V, M ) pairs that achieve a prescribed gradient-estimation accuracy with at least a minimum shot count per direction, choose the one that uses the fewest total measurements. For any zero-mean unitvariance direction distribution the minimisation yields a rule of the same form, with the component kurtosis κ appearing in place of 1; substituting Rademacher directions (κ = 1, Proposition 2) gives exactly eq. (33). The derivation has three parts. First, among isotropic independentcomponent direction distributions, Rademacher uniquely minimises the estimator second moment at fixed V (Proposition 2). Second, solving the cost-minimisation problem over (V, M ) when the accuracy constraint is tight gives the joint rule eq. (33) in closed form (Theorem 1); the two hyperparameters (α, τ 2 ) correspond to a reparameterisation of the minimum-shot floor and the accuracy target respectively. Third, the resulting shot budget matches the Cramér–Rao lower bound for unbiased gradient estimation from a shot-noise oracle, up to

Proposition 2 (Minimum fourth-moment direction distribution). Let p be an isotropic distribution on RN with independent components of zero mean and unit variance, and write κ := E[vi4 ] for the component kurtosis. Then E∥b g V − g∥2 =

∥g∥2 (κ + N − 2) . V

(34)

By E[vi4 ] ≥ 1 (Jensen, with equality iff |vi | is almostsurely constant), Rademacher directions (vi ∈ {±1}) uniquely minimise the right-hand side at fixed V , giving E∥b g V − g∥2 ≥ ∥g∥2 (N − 1)/V (proved in App. I). The result is stated within the class of isotropic independent-component distributions with unit second moment, which contains Rademacher and standard Gaussian and excludes unit-norm distributions such as directions drawn uniformly from the sphere S N −1 (the sphere distribution has dependent components). The classical forward-gradient literature [16, 21] uses both families interchangeably; in the independent-component class Rademacher is optimal and this is the class we work with throughout. Minimum-cost allocation at target accuracy. We now allow the directional derivatives to be estimated with M shots each, so the estimator picks up an addi2 tional shot-noise contribution of σ̄∇ N/(V M ) (derived in the same way as the diagonal term above with Varm [dℓ ] = 2 σ̄∇ ). Using Rademacher directions, the mean-squared error of the full estimator is 2 ∥g∥2 (N − 1) σ̄∇ N + . V VM (35) We seek the allocation that minimises the total perstep measurement cost 2V M subject to two practical requirements: the estimator MSE must not exceed a userspecified target τ 2 , and the per-direction shot count must not fall below a minimum Mmin . Formally,

MSE(V, M ) := E∥b g V − g∥2 =

min

V,M ∈R>0

2V M

s.t. MSE(V, M ) ≤ τ 2 , M ≥ Mmin .

(36)

16 Theorem 1 (Optimal (V, M ) allocation). For any τ 2 > 0, Mmin ≥ 1, and ∥g∥2 > 0, the continuous relaxation of eq. (36) has the unique minimiser 2 (N − 1)∥g∥2 + N σ̄∇ /Mmin , τ2 (37) with both constraints tight, and minimum cost 2 2 (N − 1)∥g∥2 Mmin + N σ̄∇ ⋆ ⋆ C = 2V Mmin = . (38) τ2

M ⋆ = Mmin ,

V⋆ =

bound is below 1.5% across all parameter counts at which we report experiments and below 0.5% for N ≥ 240. Remark. Corollary 1 is not a separation result: the same Fisher-information calculation applied to 2 parameter-shift gives BPS ≥ 4N σm /(π 2 τ 2 ), matching eq. (40) up to a constant, so parameter-shift is equally CRB-optimal on this oracle and both estimators lie on the same information-theoretically optimal surface. Together, Theorem 1, Proposition 2, and Corollary 1 give quiver a constructive derivation from a measurementcost optimisation problem.

2 Writing α := N σ̄∇ /(Mmin ∥g∥2 ), eq. (37) becomes

V⋆ =

(N − 1 + α) ∥g∥2 , τ2

M⋆ =

2 N σ̄∇ , α ∥g∥2

(39)

which is exactly the joint update of eq. (33). Thus (α, τ 2 ) is a reparameterisation of (Mmin , τ 2 ), and the joint rule is the constructive minimum-cost allocation (proved in App. I 1). Three remarks. First, τ 2 in eq. (39) is treated as a hyperparameter, tuned once per problem as part of a small (α, τ 2 ) grid. Second, (V ⋆ , M ⋆ ) is the unique minimiser: every other allocation satisfying the two constraints has shot cost strictly greater than C ⋆ . Third, Theorem 1 2 , ∥g∥2 ) as fixed inputs, but both change over treats (σ̄∇ training. The algorithm of Section VII C re-solves the optimum each epoch with the current EMA estimates 2 , ∥g∥2 ). At every epoch (Vt , Mt ) (b σt2 , gbt2 ) in place of (σ̄∇ is the optimum for the EMA estimates of that epoch. Matching the Cramér–Rao lower bound. The optimal cost eq. (38) is a construction: it tells us the minimum shot budget achievable by a forward-gradient estimator meeting the two practical requirements. A natural complementary question is the information-theoretic lower bound on shot cost for any unbiased estimator of g that queries a shot-noise oracle. Comparing the two tells us whether the joint rule is near-optimal among all possible estimators, not just within its own class. The Cramér–Rao inequality provides the lower bound.

Validation

We validate the joint rule eq. (33) in two steps. First, we compare quiver head-to-head against iCANS, gCANS, and plain parameter-shift (with Adam optimiser) at two (N, η) operating points: (N = 80, η = 0.05) where iCANS and gCANS can fully adapt, and (N = 160, η = 0.003) where they cannot (Fig. 8); quiver outperforms the full PS-based family in both regimes. Second, we inspect the trajectory of (Vt , Mt ) during training: both adapt over the run, and the rule is stable across a smooth hyperparameter surface (Fig. 9). Comparison against iCANS and gCANS. The iCANS/gCANS shot rule eq. (12) contains an Lη prefactor that collapses the allocation to smin whenever Lη ≪ 1; quiver’s joint (V, M ) rule eq. (33) has no such prefactor. In Fig. 8a, at (N =80, η=0.05) with Lη ≈ 0.95, iCANS and gCANS allocate above smin and adapt, but quiver still reaches a lower final energy error. In Fig. 8b, at (N =160, η=0.003) with Lη ≈ 0.06 (η is the largest stable value at this depth, and N =160 is the minimum at which exact gradient descent reaches the target; see Fig. 12), iCANS and gCANS collapse onto smin and reduce to plain parameter-shift, while quiver maintains a clear margin. The same ordering holds across a broader N sweep on n=8 and n=10 TFIM (App. J, Fig. 12) and on QAOA MaxCut (Fig. 10).

In the bias-neutral limit, the forward-gradient estimator at the optimal allocation of Theorem 1 achieves B ⋆ /BCRB = (N − 1 + α)/N → 1 as N → ∞ (proved in App. I 2).

Trajectory of (Vt , Mt ). Fig. 9 shows the joint rule on the large-N VQE benchmark (TFIM n = 8, d = 20, N = 320) at the best configuration from the (α, τ 2 ) grid (α = 0.2, τ 2 = 6400). After the 50-epoch warmup, Vt responds to the decay of gbt2 and Mt responds to the ratio σ bt2 /b gt2 ; both are shown in Fig. 9a,b. The training loss descends monotonically on both the epoch and shot axes, confirming that the warmup period and multiplicative clamp absorb the per-step EMA noise without suppressing adaptation. We selected (α, τ 2 ) from a 3 × 3 hyperparameter sweep (App. J).

The saturation rate B ⋆ /BCRB = (N − 1 + α)/N is tight at finite N : at α = 0.2 it evaluates to 0.985, 0.9925, 0.9950, 0.9963 at N = 80, 160, 240, 320 respectively, so the gap to the information-theoretic lower

Robustness to iCANS hyperparameter tuning. For QAOA on weighted MaxCut the iCANS Lipschitz constant L is upper-bounded by the cost-Hamiltonian opP erator norm Lop := ij |wij |, the safe choice of [23] used

Corollary 1 (CRB-level optimality). Let a quantum shot oracle return i.i.d. outcomes with mean f (θ +s) and 2 variance σm for shifts with ∥s∥2 ≤ ε2 N . Any unbiased estimator of g := ∇f (θ) from B queries satisfies E∥b g − g∥2 ≥

2 N σm . B ε2

(40)

(a) n = 10, d = 4, η = 0.05, 3 seeds

(b) n = 10, d = 8, η = 0.003, 3 seeds

100

N = 80 100

101

budget 5×106

100

iCANS gCANS Parameter-shift QUIVER joint (V, M)

Energy error |E − E0 |

budget 5×106

iCANS gCANS Parameter-shift QUIVER joint (V, M)

N = 160 102

103

104

Total shots

105

106

103

104

105

106

Total shots

FIG. 8. quiver outperforms iCANS and gCANS across the (N, η) regime crossover. VQE TFIM n = 10, matched 5 × 106 -shot budget, three random seeds (±1σ bands); the red vertical line marks the training cutoff. (a) N = 80, η = 0.05 (Lη ≈ 0.95): iCANS/gCANS allocate above smin and genuinely adapt; quiver leads all methods. (b) N = 160, η = 0.003 (Lη ≈ 0.06): the Lη-dependent prefactor in the iCANS/gCANS shot rule pins the allocation at smin , so both reduce to plain parameter-shift with identical trajectories; quiver leads by a margin well outside the seed-to-seed spread.

(a) Vt , Mt vs. epoch

40 20

102

2000

Epoch

Energy error

(16, 3)

10 2

(d) loss vs. shots

Energy error 0

(14, 3)

10 1

(12, 3)

101

100

(n, d)

(10, 3)

gb2t σ 2t

103

EMA

(8, 2)

Residual gap 1 − AR

Vt Mt

Vt , M t

(b) signal / noise EMAs

iCANS (best of LR grid) QUIVER joint (best of α grid) 34

100

10 1

100

103

Total shots

106

FIG. 9. quiver trajectory on the large-N VQE benchmark. TFIM n = 8, d = 20 (N = 320), three random seeds (faded lines) and their mean (bold). (a) Both Vt and Mt adapt after the 50-epoch warmup, spanning an order-of-magnitude dynamic range over the run. (b) The signal EMA gbt2 decays while the noise EMA σ bt2 is nearly flat, driving the joint update. (c,d) Training loss on epoch and shot axes: monotone descent across all seeds.

Number of parameters N

FIG. 10. quiver on QAOA MaxCut. Residual gap 1 − AR (log scale, ±1σ) vs. parameter count N at matched 106 shot budget, averaged over five weighted Erdős–Rényi graphs and three seeds. quiver leads iCANS at every N tested and maintains tight across-sample bands throughout, confirming the result is not specific to VQE or the TFIM Hamiltonian.

in Fig. 10; tuning L below this bound voids the convergence guarantee.6 We swept iCANS at n = 16, d = 3 over η ∈ {0.01, 0.05, 0.1, 0.3} and L/Lop ∈ {1, 1/3, 1/10}, and

6 We use the iCANS shot-allocation rule with an Adam parameter

update rather than Algorithm 2’s per-parameter learning-rate cap ηi ≤ χ2i /(L(χ2i + ζi /si )); at L = Lop , η = 0.1 this cap would force ηi ≲ 0.03, so the deviation favours iCANS.

Gradient estimation by the parameter-shift rule has been a structural bottleneck in training parameterised quantum circuits: its per-step measurement cost scales linearly in the number of parameters and dominates the total shot budget at any practical scale. In this work we replaced that fixed O(N ) overhead with a continuously tunable lever. Random-direction forward gradients reduce the per-step cost from O(N ) to O(V ) for any V between 1 (SPSA) and N (full parameter-shift), and the quiver optimiser sets that lever automatically from the measurement statistics already produced during training. The result is a single estimator family that recovers SPSA, RCD, and parameter-shift as limiting cases of the same construction, and a single optimiser that interpolates between them in response to the running signal-to-noise ratio. Empirically, a fixed V ≪ N forward estimator reaches the accuracy of parametershift at a fraction of the shot cost across all four benchmarks tested (ECG5000 and MNIST classification, VQE TFIM, QAOA MaxCut), with the saving growing as N increases. quiver’s joint (V, M ) rule then outperforms iCANS, gCANS, and Adam-optimised parameter-shift at

matched shot budget on VQE TFIM and QAOA MaxCut, including the conservative Lη ≪ 1 regime where the iCANS/gCANS allocation collapses onto smin and reduces to fixed-shot parameter-shift. On the hardware side, forward-gradient estimators evaluate the circuit only at shifted parameters θ±εv ℓ and require no ancillas, mid-circuit measurements, or modified gate sets beyond what parameter-shift already needs. The V directional derivatives per step are independent circuit jobs and run in parallel across QPUs, so wall-clock cost is set by perstep shots rather than by V . The framework is compatible with existing PQC training pipelines and transfers directly to any variational algorithm whose loss is linear in circuit expectation values. Open questions remain in three directions. The convergence analysis of Section IV shows that the per-step saving cannot be converted into an asymptotic speedup on convex losses, yet the empirical savings on non-convex quantum learning landscapes grow with N , and a sharper characterisation of when and why this happens would close that gap. The noiseconcentration result underpinning quiver is proved here only for local Hamiltonians on bounded-depth ansätze (App. H), leaving open the regimes of chemistry Hamiltonians in the molecular-orbital basis, long-range spin systems, and non-translation-invariant problems. Finally, generalising the quiver gain decomposition to non-linear losses (such as negative log-likelihood, in the style of Refoqus [41]), and a head-to-head comparison between the joint optimum (V ⋆ , M ⋆ ) and the classical repeatedsampling optimum ℓ⋆ of [58] at matched V × M budgets, are both natural extensions.

[1] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nat Rev Phys 3, 625 (2021). [2] K. Bharti et al., Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022). [3] M. Larocca, N. Ju, D. García-Martín, P. J. Coles, and M. Cerezo, Theory of overparametrization in quantum neural networks, Nat Comput Sci 3, 542 (2023). [4] A. Delgado, F. Rios, and K. E. Hamilton, Identifying overparameterization in Quantum Circuit Born Machines (2023), arXiv:2307.03292. [5] D. García-Martín, M. Larocca, and M. Cerezo, Effects of noise on the overparametrization of quantum neural networks, Phys. Rev. Res. 6, 013295 (2024). [6] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, Connecting ansatz expressibility to gradient magnitudes and barren plateaus, PRX Quantum 3, 010313 (2022). [7] M. Schuld and N. Killoran, Is quantum advantage the right goal for quantum machine learning?, PRX Quantum 3, 030101 (2022). [8] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning representations by back-propagating errors, Nature 323, 533 (1986).

[9] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic Differentiation in Machine Learning: a Survey, Journal of Machine Learning Research 18, 1 (2018). [10] A. Abbas, R. King, H.-Y. Huang, W. J. Huggins, R. Movassagh, D. Gilboa, and J. R. McClean, On quantum backpropagation, information reuse, and cheating measurement collapse, in Advances in Neural Information Processing Systems, Vol. 36 (2023) arXiv:2305.13362. [11] J. Bowles, D. Wierichs, and C.-Y. Park, Backpropagation scaling in parameterised quantum circuits, Quantum 9, 1873 (2025). [12] B. Coyle, S. Raj, N. Mathur, E. A. Cherrat, N. Jain, S. Kazdaghli, and I. Kerenidis, Training-efficient density quantum machine learning, npj Quantum Inf. 11, 172 (2025), arXiv:2405.20237. [13] K. Chinzei, S. Yamano, Q. H. Tran, Y. Endo, and H. Oshima, Trade-off between Gradient Measurement Efficiency and Expressivity in Deep Quantum Neural Networks, npj Quantum Inf. 11, 79 (2025). [14] J. Spall, Multivariate stochastic approximation using a simultaneous perturbation gradient approximation, IEEE Transactions on Automatic Control 37, 332 (1992).

raised the per-parameter shot ceiling from smax = 200 to 1000. No configuration beat the safe baseline: across the grid, the shot count saturates at smax within a few epochs and iCANS reduces to fixed-shot parameter-shift.

VIII.

DISCUSSION AND OUTLOOK

19 [15] Z. Ding, T. Ko, J. Yao, L. Lin, and X. Li, Random coordinate descent: A simple alternative for optimizing parameterized quantum circuits, Phys. Rev. Res. 6, 033029 (2024). [16] A. G. Baydin, B. A. Pearlmutter, D. Syme, F. Wood, and P. Torr, Gradients without Backpropagation (2022), arXiv:2202.08587. [17] D. Silver, A. Goyal, I. Danihelka, M. Hessel, and H. v. Hasselt, Learning by Directional Gradient Descent, in International Conference on Learning Representations (2022). [18] F. Hanzely, K. Mishchenko, and P. Richtarik, SEGA: Variance Reduction via Gradient Sketching, in Advances in Neural Information Processing Systems, Vol. 31 (2018) arXiv:1809.03054. [19] G. Hinton, The Forward-Forward Algorithm: Some Preliminary Investigations (2022), arXiv:2212.13345. [20] L. Fournier, S. Rivaud, E. Belilovsky, M. Eickenberg, and E. Oyallon, Can Forward Gradient Match Backpropagation?, in Fortieth International Conference on Machine Learning (2023) arXiv:2306.06968. [21] M. Ren, S. Kornblith, R. Liao, and G. Hinton, Scaling Forward Gradient With Local Losses, in International Conference on Learning Representations (2023) arXiv:2210.03310. [22] L. Balles, J. Romero, and P. Hennig, Coupling Adaptive Batch Sizes with Learning Rates, in Uncertainty in Artificial Intelligence (2017) arXiv:1612.05086. [23] J. M. Kübler, A. Arrasmith, L. Cincio, and P. J. Coles, An Adaptive Optimizer for Measurement-Frugal Variational Algorithms, Quantum 4, 263 (2020). [24] A. Gu, A. Lowe, P. A. Dub, P. J. Coles, and A. Arrasmith, Adaptive shot allocation for fast convergence in variational quantum algorithms (2021), arXiv:2108.10434. [25] J. Landman, N. Mathur, Y. Y. Li, M. Strahm, S. Kazdaghli, A. Prakash, and I. Kerenidis, Quantum Methods for Neural Networks and Application to Medical Image Classification, Quantum 6, 881 (2022). [26] L. Monbroussou, J. Landman, A. B. Grilo, R. Kukla, and E. Kashefi, Trainability and Expressivity of HammingWeight Preserving Quantum Circuits for Machine Learning, Quantum 9, 1745 (2025). [27] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, in International Conference on Learning Representations (2015) arXiv:1412.6980. [28] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018). [29] A. Paszke et al., PyTorch: An Imperative Style, High-Performance Deep Learning Library (2019), arXiv:1912.01703. [30] Martín Abadi et al., TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems (2015), software available from tensorflow.org. [31] A. Griewank, K. Kulshreshtha, and A. Walther, On the numerical stability of algorithmic differentiation, Computing 94, 125 (2012). [32] J. Schmidhuber, Deep learning in neural networks: An overview, Neural Networks 61, 85 (2015). [33] A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, and J. I. Latorre, Data re-uploading for a universal quantum

classifier, Quantum 4, 226 (2020). [34] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz, Quantum Sci. Technol. 4, 014008 (2018). [35] E. Farhi and H. Neven, Classification with Quantum Neural Networks on Near Term Processors (2018), arXiv:1802.06002. [36] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Phys. Rev. A 98, 032309 (2018). [37] D. Wierichs, J. Izaac, C. Wang, and C. Y.-Y. Lin, General parameter-shift rules for quantum gradients, Quantum 6, 677 (2022). [38] O. Kyriienko and V. E. Elfving, Generalized quantum circuit differentiation rules, Phys. Rev. A 104, 052417 (2021). [39] G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, and R. M. Parrish, Local, expressive, quantum-number-preserving VQE ansätze for fermionic systems, New J. Phys. 23, 113010 (2021). [40] R. Sweke, F. Wilde, J. Meyer, M. Schuld, P. K. Faehrmann, B. Meynard-Piganeau, and J. Eisert, Stochastic gradient descent for hybrid quantum-classical optimization, Quantum 4, 314 (2020). [41] C. Moussa, M. H. Gordon, M. Baczyk, M. Cerezo, L. Cincio, and P. J. Coles, Resource frugal optimizer for quantum machine learning, Quantum Sci. Technol. 8, 045019 (2023). [42] J. C. Spall, A Stochastic Approximation Technique for Generating Maximum Likelihood Parameter Estimates, in 1987 American Control Conference (1987) pp. 1161– 1167. [43] S. Bhatnagar, H. Prasad, and L. Prashanth, Stochastic Approximation Algorithms, in Stochastic Recursive Algorithms for Optimization (Springer, 2013) pp. 17–28. [44] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Strategies for solving the Fermi-Hubbard model on nearterm quantum computers, Phys. Rev. B 102, 235122 (2020). [45] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner, Simultaneous Perturbation Stochastic Approximation of the Quantum Fisher Information, Quantum 5, 567 (2021). [46] N. Jain, B. Coyle, E. Kashefi, and N. Kumar, Graph neural network initialisation of quantum approximate optimisation, Quantum 6, 861 (2022). [47] F. Sauvage and F. Mintert, Optimal quantum control with poor statistics, PRX Quantum 1, 020322 (2020). [48] X. Bonet-Monroig, H. Wang, D. Vermetten, B. Senjean, C. Moussa, T. Bäck, V. Dunjko, and T. E. O’Brien, Performance comparison of optimization methods on variational quantum algorithms, Physical Review A 107, 032407 (2023), arXiv:2111.13454 [quant-ph]. [49] Y. Nesterov, Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems, SIAM J. Optim. 22, 341 (2012). [50] P. Richtárik and M. Takáč, Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function, Math. Program. 144, 1 (2014). [51] A. Arrasmith, L. Cincio, R. D. Somma, and P. J. Coles, Operator Sampling for Shot-frugal Optimization in Variational Algorithms (2020), arXiv:2004.06252. [52] B. van Straaten and B. Koczor, Measurement cost of metric-aware variational quantum algorithms, PRX Quantum 2, 030324 (2021).

20 [53] G. Boyd and B. Koczor, Training variational quantum circuits with CoVaR: Covariance root finding with classical shadows, Phys. Rev. X 12, 041022 (2022). [54] G. García-Pérez, M. A. C. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, Learning to measure: Adaptive informationally complete generalized measurements for quantum algorithms, PRX Quantum 2, 040342 (2021). [55] S. Pramanik and M. G. Chandra, Stochastic Shadow Descent: Training Parametrized Quantum Circuits with Shadows of Gradients (2025), arXiv:2511.12168. [56] K. Flügel, D. Coquelin, M. Götz, and C. Debus, Beyond Backpropagation: Optimization with Multi-Tangent Forward Gradients (2024), arXiv:2410.17764. [57] T. Bos and J. Schmidt-Hieber, Convergence guarantees for forward gradient descent in the linear regression model, Journal of Statistical Planning and Inference 233, 106174 (2024). [58] N. Dexheimer and J. Schmidt-Hieber, Improving the Convergence Rates of Forward Gradient Descent with Repeated Sampling (2024), arXiv:2411.17567. [59] U. Singhal, B. Cheung, K. Chandra, J. Ragan-Kelley, J. B. Tenenbaum, T. A. Poggio, and S. X. Yu, How to guess a gradient (2023), arXiv:2312.04709. [60] Z. Wang, S. Markou, and A. Campbell, Towards Scalable Backpropagation-Free Gradient Estimation (2025), arXiv:2511.03110. [61] K. Panchal, S. Choudhary, Y. Brun, and H. Guan, The Cost of Avoiding Backpropagation (2025), arXiv:2506.21833. [62] A. D. Cobb, A. G. Baydin, B. A. Pearlmutter, and S. Jha, Second-Order Forward-Mode Automatic Differentiation for Optimization, in International Conference on Learning Representations (2025) arXiv:2408.10419.

[63] Y. Yu, R. Xia, Q. Ma, M. Lengyel, and G. Hennequin, Second-Order Forward-Mode Optimization of Recurrent Neural Networks for Neuroscience, in Advances in Neural Information Processing Systems, Vol. 37 (2024). [64] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, Quantum Natural Gradient, Quantum 4, 269 (2020). [65] A. Mari, T. R. Bromley, and N. Killoran, Estimating the gradient and higher-order derivatives on quantum hardware, Physical Review A 103, 012405 (2021). [66] R. M. Parrish, G.-L. R. Anselmetti, and C. Gogolin, Analytical Ground- and Excited-State Gradients for Molecular Electronic Structure Theory from Hybrid Quantum/Classical Methods (2021), arXiv:2110.05040. [67] M. M. Wolf, Mathematical Foundations of Supervised Learning (Lecture notes, Technical University of Munich, 2023). [68] M. Talagrand, Concentration of measure and isoperimetric inequalities in product spaces, Publications Mathématiques de l’IHÉS 81, 73 (1995). [69] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 1791 (2021). [70] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardwareefficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017). [71] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014), arXiv:1411.4028 [quant-ph]. [72] R. Herrman, P. C. Lotshaw, J. Ostrowski, T. S. Humble, and G. Siopsis, Multi-angle quantum approximate optimization algorithm, Scientific Reports 12, 6781 (2022), arXiv:2109.11455.

Appendix A: Unbiasedness of the forward gradient estimator

We prove that the V -direction, M -shot forward gradient estimator eq. (A2) is unbiased in the ε → 0 limit, adapting the classical argument of [16] to the quantum finite-difference setting with shot noise. Proposition 3 (Unbiasedness in the ε → 0 limit). Let f (θ) = ⟨ψ(θ)|O|ψ(θ)⟩ be a circuit expectation value, and write õm (θ ′ ) for the outcome of the m-th single-shot measurement of O on |ψ(θ ′ )⟩, satisfying Em [õm (θ ′ )] = f (θ ′ ). Let {v ℓ }Vℓ=1 be drawn i.i.d. from a distribution p on RN satisfying E[v] = 0 and E[vv ⊤ ] = IN . Suppose the central finite-difference stencil is consistent, in the sense that for every measurement realisation m, lim

ε→0

Q X

γqv,ε õm θ εv,q = ∇v fm (θ) = v ⊤ ∇fm (θ),

(A1)

q=1

where fm denotes the single-shot estimator. Then the V -direction, M -shot forward gradient estimator Q V M 1 X 1 X X vℓ ,ε eF (θ) = γq õm θ εvℓ ,q v ℓ g V M m=1 q=1

(A2)

ℓ=1

F e (θ) = ∇f (θ), with expectation taken jointly over directions and measurement noise. satisfies limε→0 E g Proof. By the tower property, E[e g F ] = Ev Em [e g F | v] . The inner expectation over the unbiased single-shot quantum estimator Em [õm (θ)] = ⟨ψ(θ)|O|ψ(θ)⟩ = f (θ), combined with stencil consistency eq. (A1), gives V F 1 X ℓ⊤ e |v = v ∇f (θ) v ℓ . lim Em g ε→0 V ℓ=1

(A3)

21 Taking the expectation over directions and using the i.i.d. assumption together with the second-moment identity E[vv ⊤ ] = IN , Ev (v ⊤ ∇f ) v = Ev [vv ⊤ ] ∇f = ∇f (θ), (A4) which is independent of V . Hence limε→0 E[e g F ] = ∇f (θ). The proposition extends to losses linear in circuit observables (such as the MSE loss used in the main text) by PD linearity of expectation: the gradient of L(θ) = 1/D d=1 ℓd (θ) inherits unbiasedness from each term ℓd by linearity. The finite-ε bias is O(ε2 ) for the central stencil [40]; the bias–variance balance leading to the ε⋆ ∝ M −1/6 rule is derived in App. B.

Appendix B: Bias–variance trade-off and the optimal step size

We derive the closed-form optimum ε⋆ that balances finite-difference bias against shot-noise amplification, and verify it against the empirical U-curve of Fig. 11. 2 /(2M ε2 ) (derived Variance. The central-difference estimator divides by 2ε, so shot noise enters the variance as σm 2 in full below as eq. (B2)), amplified as 1/ε , and taking ε large suppresses noise but introduces O(ε2 ∥∇3 f ∥) Taylorremainder bias. Because quantum expectation values are band-limited trigonometric polynomials [66], the bias term remains small up to surprisingly large ε, and the optimal choice is set almost entirely by the noise amplification.

Closed-form ε⋆ . expectation

Expanding f (θ ± εv) to fourth order around θ, the central difference estimator eq. (16) has h i 2 e εv f = (∇f · v) + ε ∂v3 f + O(ε4 ), Em ∇ 6

(B1)

P where ∂v3 f := i,j,k vi vj vk ∂i ∂j ∂k f is the third directional derivative. The leading bias is therefore ε2 C3 /6 where C3 := |∂v3 f | depends on v and θ but not on ε, while the variance is i h 2 e εv f = σm , Varm ∇ 2M ε2

(B2)

2 where σm is the per-shot variance of the observable at θ±εv (upper-bounded by the operator norm of the Hamiltonian). The mean-squared error of the estimator is the sum of squared bias and variance,

MSE(ε) =

ε4 C32 σ2 + m2. 36 2M ε

(B3)

Minimising over ε, the first-order condition d/dε MSE = 0 yields the closed-form optimum ε⋆ =

2 9 σm M C32

1/6

(B4)

2 with minimum MSE MSE(ε⋆ ) = 1/2 [3−1/3 + 32/3 ] (C3 σm /M )2/3 .

Scaling properties. Three properties follow directly from eq. (B4): (i) ε⋆ ∝ M −1/6 , so doubling M decreases the optimal step by only ≈ 11% and transferring ε⋆ between runs of different M is safe over the factor-of-ten range 1/3 typical in practice; (ii) ε⋆ ∝ σm , so a noisier observable calls for a larger step since more shot noise can be absorbed −1/3 ⋆ into a larger denominator; , so a flatter (more band-limited) landscape tolerates a larger step. Band√ (iii) ε ∝ C3 limitation to frequency N [66] means C3 = O(N 3/2 ∥H∥) in the worst case, so ε⋆ is weakly N -dependent, decreasing as N −1/2 .

VQE TFIM n = 8, d = 20 (N = 320), matched 5M-shot budget Energy error |Ebest E0|

Calibration and transferability. Fig. 11 shows the empirical U-curve on a representative VQE problem: at fixed V = 10 and matched shot budget, the best energy error is achieved at ε = 0.1, with a near-3× degradation at ε = 0.01 (noise-dominated) and a 5× degradation at ε = 0.5 (bias-dominated). The closed form eq. (B4) is most usefully read as a transferability tool rather than as a calibration-free predictor: the constant C3 is problemspecific and is recovered from a single U-curve fit (here, 2 C3 ∼ 103 at M = 50 and σm ∼ ∥H∥2 ∼ 10, recovering ⋆ the empirical ε ≈ 0.1 used throughout the main text), but once C3 is calibrated at any one operating point the scal1/3 ing laws ε⋆ ∝ M −1/6 , σm , N −1/2 make the optimal step transferable across (M, σm , N ) without re-running the full sweep.

Forward V = 10, M = 50 Optimum: = 0.1, err = 0.109

1.0 0.8 noise 1/ 2

bias

0.6

0.4 0.2 10 2

10 1

Finite-difference step

FIG. 11. ε = 0.1 sits at the bias–variance optimum. VQE TFIM n=8, d=20, V =10, M =50, matched 5M-shot budget, seed 0. Noise amplification (∝ 1/ε2 ) dominates at small ε; finite-difference bias (∝ ε2 ) dominates at large ε.

Appendix C: Second-moment expansion of the forward gradient estimator

We prove Lemma 1, which bounds the second moment of the V -direction forward gradient estimator and underlies the convergence rates of Section IV. Lemma 4. Let {v ℓ }Vℓ=1 be drawn i.i.d. with E[viℓ ] = 0, E[viℓ viℓ′ ] = δii′ , and let κ := E[(viℓ )4 ] be the fourth moment of e F L(θ) = 1/V PV (∇ e vℓ LM ) v ℓ be one component (so κ = 1 for Rademacher and κ = 3 for standard Gaussian). Let ∇ ℓ=1 the V -direction M -shot forward gradient estimator. Then h Ev

i e M 2, e F L(θ) 2 = N + V + κ − 2 ∇L ∇ V

(C1)

e M is the (full-coordinate) M -shot estimator of the true gradient. In the noiseless limit M → ∞ this reduces where ∇L F e L 2 → (N +V +κ−2)/V · ∥∇L∥2 , which is the bound used in part 2 of Lemma 5 below. to ∇ e F L]k = 1/V M PM PV PN v ℓ (∂j fm ) v ℓ , where ∂j fm is the (single-shot) partial derivative Proof. Write Zk := [∇ k m=1 ℓ=1 j=1 j of the loss in coordinate j. Squaring and summing, Ev ∥Z∥2 =

N M 1 X X V 2M 2 ′

V N X X

h i ′ ′ (∂j fm )(∂j ′ fm′ ) E vjℓ vkℓ vjℓ′ vkℓ .

k=1 m,m =1 ℓ,ℓ′ =1 j,j ′ =1

′

The fourth-moment object E[vjℓ vkℓ vjℓ′ vkℓ ] depends on whether ℓ = ℓ′ and on the index pattern in (j, j ′ , k): ′

′

• ℓ ̸= ℓ′ . By independence, E[vjℓ vkℓ ]E[vjℓ′ vkℓ ] = δjk δj ′ k . • ℓ = ℓ′ . The expectation is E[vjℓ vjℓ′ (vkℓ )2 ], with three non-vanishing index patterns: – j = j ′ = k: E[(vkℓ )4 ] = κ; – j = j ′ ̸= k: E[(vjℓ )2 (vkℓ )2 ] = 1; – j ̸= j ′ , or any pattern where a lone v component appears: 0 by zero-mean.

(C2)

23 Substituting these cases,

Ev ∥Z∥

1 = 2 2 V M

M X

X N N X V κ (∂k fm )(∂k fm′ ) + (N − 1) (∂j fm )(∂j fm′ )

m,m′ =1

j=1

k=1

+ (V − 1)

N X

(∂k fm )(∂k fm′ )

(C3)

k=1

N M N +V +κ−2 1 X X (∂k fm )(∂k fm′ ) V M2 ′

(C4)

N +V +κ−2 e M 2 , ∇L V

(C5)

k=1 m,m =1

e M by the true gradient ∇L. which is eq. (C1). Taking M → ∞ replaces ∇L For Rademacher directions κ = 1 the prefactor reduces to (N + V − 1)/V ; for standard Gaussian κ = 3 it becomes (N + V + 1)/V . Both choices are bounded above by (N + V + 1)/V , so the convergence rate eq. (19) obtained by substituting β 2 = (N + V + κ − 2)/V into Lemma 5 (App. D) is robust across reasonable direction distributions.

Appendix D: Polyak–Łojasiewicz inequality and stochastic gradient descent

We state the Polyak–Łojasiewicz (PL) inequality and the PL-SGD convergence theorem [40, 67], then prove Proposition 1 by substituting the second-moment bound of App. C into the SGD bound. The PL inequality. A function f ∈ C 1 (RN ) that attains its global minimum at θ ⋆ is said to satisfy the (global) PL inequality with constant µ > 0 if ⋆ 2 1 ∀θ ∈ RN . (D1) 2 ∥∇f (θ)∥ ≥ µ f (θ) − f (θ ) The PL condition is in general not satisfied globally for parameterised quantum circuits [15], but it is often a reasonable assumption locally. Writing X = {θ ⋆ } for the set of global minima, the local PL condition states that there exist δf , µ > 0 such that for all θ in the sublevel neighbourhood N (X ) := f −1 ([f (θ ⋆ ), δf ]), 2 1 2 ∥∇f (θ)∥ ≥ µ

f (θ) − f (θ ⋆ )

∀θ ∈ N (X ).

(D2)

The local condition is what is invoked implicitly throughout Section IV: convergence is to a neighbourhood of a minimum, not to the global minimum. Proposition 4 (Proposition 1 restated). Let f ∈ C 1 (RN ) be L-smooth and satisfy the PL inequality with constant µ > 0, and let directions be drawn from the Rademacher distribution (κ = 1). (i) Convergence of forward gradient descent. With learning rate η = V /(L(N + V − 1)), forward gradient descent satisfies E[f (θ (T ) )] − f ⋆ ≤

1−

µV L(N +V −1)

f (θ (0) ) − f ⋆ .

(D3)

2 (ii) Convergence under shot noise. With per-direction shot noise variance σshot and η ∈ [0, 1/(2µ)],

E[f (θ (T ) )] − f ⋆ ≤ (1 − 2µη)

2 Lη(N +V −1) σshot f (θ (0) ) − f ⋆ + . 4µV M

(D4)

The proof uses the following cited lemma. Lemma 5 (SGD under PL and bounded second moment, [40, 67]). Let f ∈ C 1 (RN ) be L-smooth, satisfy the PL inequality eq. (D1) with constant µ > 0, and attain its global minimum at θ ⋆ . Let g (1) , . . . , g (T ) be unbiased gradient estimators and consider the iterates θ (t+1) = θ (t) − η g (t) (θ (t) ).

24 1. If E[∥g (t) (θ)∥2 ] ≤ γ 2 for all θ, t and η ∈ [0, 1/(2µ)], then Lη γ 2 . E[f (θ (T ) )] − f (θ ⋆ ) ≤ (1 − 2µη)T f (θ (0) ) − f (θ ⋆ ) + 4µ

(D5)

2. If E[∥g (t) (θ)∥2 ] ≤ β 2 ∥∇f (θ)∥2 for all θ, t and η = 1/(Lβ 2 ), then µ T f (θ (0) ) − f (θ ⋆ ) . E[f (θ (T ) )] − f (θ ⋆ ) ≤ 1 − 2 Lβ

(D6)

Proof of Proposition 4. Part (i). By App. A, E[e g F (θ)] = ∇f (θ), so the estimator is unbiased. By Lemma 4 with κ = 1 (Rademacher), F 2 N + V − 1 E ∥e g ∥ = ∥∇f ∥2 =: β 2 ∥∇f ∥2 . V This is the bounded relative second moment condition of part 2 of Lemma 5. Setting η = 1/(Lβ 2 ) = V /(L(N +V −1)) and substituting β 2 = (N + V − 1)/V , T µ T µV (0) ⋆ (f (θ (0) ) − f ⋆ ), E[f (θ (T ) )] − f ⋆ ≤ 1 − (f (θ ) − f ) = 1 − Lβ 2 L(N + V − 1) which is part (i).

e F LM adds a shot-noise contribution. By the variance Part (ii). With M shots per direction the M -shot estimator ∇ decomposition (proof of Lemma 4), the full second moment is 2 F 2 N + V − 1 N + V − 1 σshot N +V −1 E ∥e g ∥ = ∥∇f ∥2 + · ≤ ∥∇f ∥2 + γ 2 , V V M V

2 where γ 2 := (N + V − 1)σshot /(V M ) is the shot-noise residual. This satisfies the bounded absolute second moment condition of part 1 of Lemma 5 with this γ 2 . Substituting into part 1 gives

E[f (θ (T ) )] − f ⋆ ≤ (1 − 2µη)T (f (θ (0) ) − f ⋆ ) +

2 Lηγ 2 Lη(N + V − 1)σshot = (1 − 2µη)T (f (θ (0) ) − f ⋆ ) + , 4µ 4µV M

which is part (ii). Appendix E: Expected gain decomposition and per-direction shot allocation

We derive the per-direction gain decomposition eq. (30) for forward gradient descent, then derive the optimal per-direction shot allocation Mℓ⋆ that maximises gain per measurement. Proof strategy.

The argument proceeds in four steps.

1. Simplify the gain via unbiasedness. Starting from the one-step expected loss decrease eq. (10), the cross-term e F L] reduces to ∥∇L∥2 by unbiasedness, leaving the second moment E[∥∇ e F L∥2 ] as the only quantity to ∇L⊤ E[∇ bound. e F L∥2 ] into a sum over 2. Decompose the second moment per direction. Lemma 6 (proved below) splits E[∥∇ e vℓ Lm ]/M . directions, each contributing a signal term (∇vℓ L)2 and a shot-noise term Varm [∇ 3. Derive the per-direction gain formula. Substituting Steps 1–2 into the gain and grouping by direction yields eq. (E3), from which the learning-rate criterion follows directly. 4. Optimise the shot allocation. Treating the per-direction gain as a function of Mℓ alone and solving the firstorder condition gives the closed-form Mℓ⋆ eq. (E5). Under noise concentration (Assumption 1) the same approach applied to V at fixed M gives the fixed-M optimal V ⋆ eq. (31).

25 Notation. Throughout this appendix we work with two sources of randomness: the random forward directions {v ℓ }Vℓ=1 (drawn i.i.d. with zero mean, unit variance, kurtosis κ) and the per-shot quantum measurement noise. We write Ev , Em for the corresponding expectations, and E = Ev Em for the joint expectation. For a directional derivative e v LM = 1/M P ∇ e we write ∇v L for the true value and ∇ m v Lm for the M -shot Monte Carlo estimator. Setup. Starting from eq. (10), the unbiasedness of the forward gradient estimator (Proposition 3) implies e F L] = ∥∇L∥2 , so the gain expression simplifies to ∇L⊤ E[∇ E[G F ] = η ∥∇L∥2 −

Lη 2 h e F 2 i E ∇ L . 2

(E1)

The second-moment term is exactly the object bounded by Lemma 1 of App. C.

Variance decomposition over measurements

To express the gain in a form where each random direction contributes a separate signal and noise term we decompose the measurement-side expectation of the directional-derivative variance. e vℓ Lm ] = ∇vℓ L Lemma 6 (Variance-with-measurement decomposition). With unbiased single-shot estimators Em [∇ and i.i.d. measurement trials, V h i h i 2 1 X 1 M e e Em Varv ∇v L ≈ ∇v ℓ L + Varm ∇vℓ Lm . V M

(E2)

ℓ=1

e v LM = Ev [(∇ e v LM )2 ]. Proof. Since the directional derivatives have zero mean over directions (by isotropy), Varv ∇ Approximate this expectation by the empirical average over the V sampled directions: V X e vℓ LM 2 . e v LM )2 ≈ 1 ∇ Ev (∇ V ℓ=1

e vℓ LM = 1/M Expand each squared term using ∇

m=1 ∇v ℓ Lm :

e v ℓ LM 2 = 1 ∇ M2

M X

e vℓ Lm ∇ e vℓ Lm′ . ∇

m,m′ =1

Taking Em of each summand splits into two cases. Off-diagonal (m ̸= m′ ). The M (M − 1) such terms factorise by independence of the two shots: h i e v ℓ Lm · ∇ e vℓ Lm′ = ∇vℓ L 2 . Em ∇ Diagonal (m = m′ ). The M such terms use Em [X 2 ] = (Em X)2 + Varm X: h Em

h i i e vℓ Lm 2 = ∇vℓ L 2 + Varm ∇ e vℓ Lm . ∇

Combining diagonal and off-diagonal contributions and dividing by M 2 : h Em

i h i e vℓ Lm ] + M (M − 1)(∇vℓ L)2 e vℓ LM 2 = 1 M (∇vℓ L)2 + Varm [∇ ∇ M2 h i 2 1 e v ℓ Lm . = ∇v ℓ L + Varm ∇ M

Averaging over ℓ yields eq. (E2).

26 2.

Per-direction gain and learning-rate criterion

Substituting Lemma 1 and Lemma 6 into eq. (E1): E[G F ] = η ∥∇L∥2 −

Lη 2 h e F 2 i E ∥∇ L∥ 2

V e v ℓ Lm ] Lη 2 N + V + κ − 2 1 X Varm [∇ · · (∇vℓ L)2 + 2 V V M ℓ=1 V e vℓ Lm ] 1 X Varm [∇ Lη 2 N + V + κ − 2 , = η ∥∇L∥2 − (∇vℓ L)2 + V 2 V M ℓ=1 | {z } =: γvℓ

≈ η ∥∇L∥2 −

where the second line uses Lemma 1 for the second-moment term and Lemma 6 for the per-direction variance, and the third line folds the constant signal term η∥∇L∥2 inside the sum. This gives the per-direction decomposition V

E[G F ] ≈

1 X γv ℓ , V

γvℓ := η ∥∇L∥2 −

ℓ=1

Lη 2 N + V + κ − 2 2 V

i h 2 1 e vℓ Lm , ∇v ℓ L + Varm ∇ M

(E3)

which is eq. (30) of the main text. Each direction contributes a (constant) signal term η ∥∇L∥2 and a penalty composed of a directional-derivative term (∇vℓ L)2 (the part of the gradient projected onto v ℓ ) and a measurement-noise term e vℓ Lm ]/M . Varm [∇ Requiring E[G F ] > 0 (i.e. the loss decreases in expectation) gives the learning-rate criterion η<

2V ∥∇L∥2 , e v LM L(N + V + κ − 2) Em Varv ∇

(E4)

the forward-gradient analogue of the CABS [22] and iCANS [23] learning-rate criteria. 3.

Optimal per-direction shot allocation

Allowing the number of shots to depend on the direction, M → Mℓ , the per-direction gain γvℓ from eq. (E3) becomes a function of Mℓ alone. Maximising the gain-per-shot γvℓ /Mℓ over Mℓ and rearranging yields the optimal per-direction allocation referenced from Section VII A. Lemma 7 (Optimal per-direction shot allocation). Let L have L-Lipschitz gradients and let {v ℓ }Vℓ=1 have zero mean, unit variance, and kurtosis κ. Then the per-direction shot count that maximises the per-shot gain is h i e vℓ Lm 2 Varm ∇ (E5) Mℓ⋆ = 2 . 2V ∥∇L∥2 − ∇vℓ L Lη(N + V + κ − 2) 2 e vℓ Lm ], gℓ := (∇vℓ L)2 , and c := Lη2/2 · N +V +κ−2/V . Then the per-shot gain is Since γvℓ Proof. Write σ∇,ℓ := Varm [∇ contains a 1/Mℓ term from the noise, the gain-per-shot is

i c σ2 γv ℓ 1 h ∇,ℓ = η∥∇L∥2 − c gℓ − . Mℓ Mℓ Mℓ2 2 Differentiating with respect to Mℓ (treating ∥∇L∥2 , gℓ , σ∇,ℓ as Mℓ -independent): 2 2c σ∇,ℓ d γv ℓ η∥∇L∥2 − c gℓ =− + = 0. dMℓ Mℓ Mℓ2 Mℓ3

Multiplying through by Mℓ3 > 0: 2 2c σ∇,ℓ = Mℓ η∥∇L∥2 − c gℓ .

27 Solving for Mℓ and re-expanding c = Lη /2 · N +V +κ−2/V : 2

Mℓ⋆ =

2 2 e vℓ Lm ] 2 σ∇,ℓ 2c σ∇,ℓ 2 Varm [∇ , = = 2 2 2V η∥∇L∥ − c gℓ η∥∇L∥ 2 − (∇ ℓ L)2 ∥∇L∥ − gℓ v Lη(N + V + κ − 2) c

which is eq. (E5). The denominator of eq. (E5) can be negative when the directional derivative is well-aligned with the true gradient; the regime analysis of when this happens, and the practical mitigation, is the subject of App. F.

Fixed-M optimal V

e vℓ Lm ] ≈ σ̄ 2 for all ℓ. For Under Assumption 1 the per-direction measurement variance concentrates, Varm [∇ ∇ 2 2 isotropic zero-mean unit-variance directions, Ev [(∇vℓ L) ] = ∥∇L∥ . Taking this expectation in the per-direction gain eq. (E3) and averaging over the V directions: E[G F ] ≈ η ∥∇L∥2 −

Lη 2 N + V + κ − 2 σ̄ 2 ∥∇L∥2 + M∇ . 2 V

(E6)

With M fixed, the gain-per-shot objective is f (V ) := E[G F ]/(2V M ). Define the shorthand 2 2 B := Lη2 ∥∇L∥2 + σ̄∇ /M ,

A := η∥∇L∥2 ,

so that eq. (E6) reads E[G F ] ≈ A − B(N + V + κ − 2)/V = (A − B) − B(N + κ − 2)/V . Dividing by 2V M : f (V ) =

B(N + κ − 2) A−B − . 2V M 2V 2 M

Differentiating with respect to V : f ′ (V ) = −

B(N + κ − 2) A−B + . 2V 2 M V 3M

Setting f ′ (V ⋆ ) = 0 gives (A − B)/(2(V ⋆ )2 M ) = B(N + κ − 2)/(V ⋆ )3 M , so V ⋆ (A − B) = 2B(N + κ − 2) and V⋆ =

2B(N + κ − 2) . A−B

Re-substituting A and B and cancelling Lη 2 /2 from numerator and denominator: V⋆ =

2 Lη 2 (∥∇L∥2 + σ̄∇ /M )(N + κ − 2) 2

2 /M ) η∥∇L∥2 − Lη2 (∥∇L∥2 + σ̄∇

which after dividing numerator and denominator by η yields eq. (31). This is a maximum (not minimum) since f (V ) → −∞ as V → 0+ and f (V ) → 0− as V → ∞. The denominator A − B > 0 is precisely the learning-rate criterion eq. (E4); when it fails the gain-per-shot is decreasing in V for all V > 0, and V ⋆ is set to Vmax = N .

Appendix F: Regime analysis of per-direction shot allocation

The per-direction allocation Mℓ⋆ of Lemma 7 is well-defined only when its denominator is positive. We analyse when this fails and what it means operationally.

28 bℓ with ∥b Alignment decomposition. Decompose an arbitrary direction as v ℓ = ∥v ℓ ∥ v v ℓ ∥ = 1, and let ϕℓ be the ℓ b . Then angle between ∇L and v 2 2 ∇vℓ L = ∇L · v ℓ = ∥∇L∥2 ∥v ℓ ∥2 cos2 ϕℓ . Defining the dimensionless constant C := 2V /(Lη(N + V + κ − 2)), the denominator of eq. (E5) becomes Dℓ = ∥∇L∥2 C − ∥v ℓ ∥2 cos2 ϕℓ .

(F1)

(F2)

The sign of Dℓ depends on the product ∥v ℓ ∥2 cos2 ϕℓ : even at moderate misalignment, a large vector norm ∥v ℓ ∥ can drive Dℓ < 0 purely from scale. Restricting to unit-norm directions (∥b v ℓ ∥ = 1) decouples scale from alignment, (F3)

Dℓ = ∥∇L∥2 (C − cos2 ϕℓ ), so the sign depends purely on alignment and on (V, L, η, N, κ). Three regimes.

For unit-norm directions:

• C > 1: Dℓ > 0 for every v ℓ , so Mℓ⋆ > 0 and grows with alignment (∂Mℓ⋆ /∂ cos2 ϕℓ > 0). Better-aligned directions receive more shots, as one would expect intuitively. • C = 1: Dℓ ≥ 0, with equality only at perfect alignment (cos2 ϕℓ = 1). Allocation diverges at exact alignment but is finite elsewhere. • C < 1: cones around ±∇L where cos2 ϕℓ > C have Dℓ < 0, formally requiring negative shot counts. The third regime occurs when an aggressive choice of η, L, or small V drives the per-direction gain-per-shot negative, signalling that the squared projection (∇L · v ℓ )2 already dominates the gain expression and additional shots in that direction provide no benefit. The remedy is to push C ≥ 1 (e.g. by tightening the learning-rate criterion eq. (E4)) or to enforce a minimum allocation Mℓ ≥ Mmin . For a representative configuration (κ = 3, η = 0.1, N = 100, V = 20, L = 1) we have C ≈ 3.31, comfortably in the first regime. b = v/∥v∥ rather than raw i.i.d. vectors Normalised random directions. Working with normalised directions v PN 2 changes the per-coordinate moments. Since i=1 vbi ≡ 1, the i.i.d. assumption gives E[b vi2 ] =

1 , N

Var(b vi ) =

2 1 − E[b vi ] , N

(F4)

and zero-mean for normalised coordinates is not alone. A sufficient extra condition is √ guaranteed by zero-mean of vP symmetry under v 7→ −v, which makes t 7→ t/ t2 + S an odd map (with S = k̸=i vk2 ) and yields E[b vi ] = 0. Under F b = N (∇f · v b) v b recovers this assumption the second moment is 1/N , so a normalised forward gradient estimator g unbiasedness with the rescaling factor N . Appendix G: Gradient-Lipschitz constant of the VQE loss

The gradient-Lipschitz constant L of the VQE loss bounds the step size for convergence and enters the iCANS/gCANS shot-allocation multiplier. We derive it here for TFIM with J = h = 1 and open boundary conditions so that it is available to the noise-concentration proof that follows. The VQE loss is L(θ) := ⟨ψ(θ)|H|ψ(θ)⟩ with |ψ(θ)⟩ := U (θ)|0⟩⊗n . The ansatz is hardware-efficient with each layer a product of single-qubit rotations e−iθj Pj /2 for Pj ∈ {X, Y } followed by a fixed CZ ring. All generators satisfy G2j = I and ∥Gj ∥ = 1. For two parameter vectors θ 1 , θ 2 , ∥∇L(θ 1 ) − ∇L(θ 2 )∥ ≤ L ∥θ 1 − θ 2 ∥,

L ≤ ∥H∥,

which follows from standard parameter-shift arguments: the j-th gradient component is ∂j L = i/2⟨ψ(θ)|[U † HU , G ]|ψ ⟩ with |ψ ⟩ the state at parameter j, so |∂ ∂ L| ≤ ∥[[H, G ], G ]∥ ≤ 4∥H∥ ∥G ∥ ∥G ∥ = >j j <j <j i j i j i j >j 4∥H∥, and the operator-norm bound on the Hessian gives L ≤ ∥H∥ (with the factor of 4 absorbed by the half-angle rotation convention; see [23, 40] for the same convention).

29 For TFIM, HTFIM = −J

n−1 X i=1

Zi Zi+1 − h

n X

Xi ,

i=1

the operator norm is ∥H∥ ≤ (n − 1)|J| + n|h| by the triangle inequality. At n = 8, J = h = 1, OBC, this gives ∥H∥ ≤ 15, and so L ≤ 15. The bound is tight within an O(1) factor. Appendix H: Proof of the noise-concentration assumption for local Hamiltonians

We first restate the assumption from Section VII B, then prove it for the circuits and Hamiltonians used in this paper via the following lemma. 2 e v ℓ Lm ] Assumption 2 (Assumption 1 restated). Let v ℓ have i.i.d. Rademacher components and let σ∇,ℓ := Varm [∇ 2 2 2 ℓ denote the per-direction measurement variance, with mean σ̄∇ := Ev [σ∇,ℓ ]. The function σ∇,ℓ (v ) is Lipschitz in v ℓ 2 2 with Lipschitz constant Lσ independent of N , so that σ∇,ℓ ≈ σ̄∇ for all ℓ with fluctuations O(Lσ ). This is verified empirically for the benchmarks of this paper in Fig. 7.

The following lemma gives a concrete, quantitative version of this condition. PJ Lemma 8 (Noise concentration for local Hamiltonians). Let H = j=1 hj be a k-local Hamiltonian with ∥hj ∥ ≤ hmax , and let U (θ) be a 1D brick-layer hardware-efficient ansatz of depth d with N parameters and Pauli generators Gi with 2 e vℓ Lm ] for Rademacher v ℓ ∈ {−1, +1}N , and let σ̄ 2 := Ev [σ 2 ]. Then for every ∥Gi ∥ = 1. Let σ∇,ℓ (v) := Varm [∇ ∇ ∇,ℓ δ > 0, p ξ(k, d) h2max δ2 2 2 , Lσ = Pr |σ∇,ℓ − σ̄∇ | > δ ≤ 2 exp − 2 , 2Lσ εM where ξ(k, d) = kd + d(d + 1) is the parameter count in the backward light cone of a k-local term at depth d, and Lσ is independent of N at fixed (k, d). In particular, Assumption 1 holds with Lipschitz constant Lσ . Setting. Throughout the proof, O is a Pauli operator entering one of the hj terms, U (θ) = UN · · · U1 with Ui = e−iθi Gi /2 , and ⟨O⟩(θ) := ⟨ψ(θ)|O|ψ(θ)⟩. 2 2 Variance formula for σ∇,ℓ . Write σ± (v) := Varm [om (θ ± εv)] for the per-shot variance at the two shifted parameter vectors. The variance of the central-difference estimator is then M 2 (v) σ 2 (v) + σ− f (θ + εv) − f M (θ − εv) 2 σ∇,ℓ := Var = + . (H1) 2ε 4ε2 M

Each measurement of O returns an eigenvalue, which by definition lies in [λmin , λmax ]. Popoviciu’s inequality (Var[X] ≤ 2 (b − a)2 /4 for X ∈ [a, b]) then gives σ± ≤ 1/4(λmax − λmin )2 , so 2 σ∇,ℓ (v) =

2 2 σ+ (v) + σ− (v) (λmax − λmin )2 ≤ . 2 4ε M 8ε2 M

(H2)

2 Proof of Lemma 8. Since v ℓ has N i.i.d. bounded (±1) components and σ∇,ℓ is a smooth function of them, the bound 2 follows from Talagrand’s inequality once we establish that σ∇,ℓ is Lipschitz with an N -independent constant.

Lemma 9 (Talagrand’s concentration inequality for Lipschitz functions, [68]). Let X1 , . . . , XN be independent random variables with |Xi | ≤ 1, and let f : RN → R be L-Lipschitz with respect to the Euclidean norm. Then for any t > 0, t2 Pr |f (X1 , . . . , XN ) − E[f ]| > t ≤ 2 exp − 2 . 2L 2 We apply this to f = σ∇,ℓ as a function of the N i.i.d. ±1 Rademacher components of v ℓ . The lemma requires 2 σ∇,ℓ to be Lipschitz; the Lipschitz constant of a differentiable function is its gradient-norm supremum, so we need an 2 N -independent bound on ∥∇v σ∇,ℓ ∥, which we derive in three steps.

30 1. Light-cone bound [69]. Write U (θ) = UN UN −1 · · · U1 where each Ui = e−iθi Gi /2 is a single-qubit rotation. The expectation value decomposes as ⟨O⟩(θ) = ⟨0|⊗n U † (θ) O U (θ)|0⟩⊗n . Gate Ui acts on qubit qi . By the product rule, ∂⟨O⟩/∂θi involves commuting Gi through the subsequent gates to the measurement. If qubit qi is never in the support of any gate between layer i and the support of O, then Gi acts as the identity on the relevant subspace and ∂⟨O⟩/∂θi = 0 exactly. Following Cerezo et al. [69], the backward light cone LB (O) is the set of gates whose output qubits are causally connected to the support of O; all parameters outside LB (O) have zero gradient. For a k-local term hj in a 1D brick-layer ansatz of depth d, each layer extends LB (hj ) by one qubit in each direction, so at layer l (counted back from the measurement) the cone has width k + 2l. Summing over all d layers gives the total parameter count: ξ(k, d) =

d X

(H3)

(k + 2l) = kd + d(d + 1) = O (k + d) d ,

l=1

which is independent of N at fixed (k, d). 2 2 2. Chain rule for ∇v σ± . Since σ± (v) = Varm [om (θ ± εv)] and Varm [om (θ ′ )] = ⟨O2 ⟩(θ ′ ) − ⟨O⟩2 (θ ′ ), differentiating in vi gives 2 ∂σ± ∂ 2 ⟨O ⟩ − ⟨O⟩2 θ±εv = ∓2ε ⟨O⟩(θ ± εv) ∂θi ⟨O⟩(θ ± εv). = ±ε ∂vi ∂θi

(H4)

2 By the light-cone bound of step (i), ∂vi σ± = 0 for all θi ∈ / LB (hj ), leaving at most ξ(k, d) eq. (H3) non-zero components. For the in-cone indices, we bound each factor separately:

• |∂θi ⟨O⟩| ≤ ∥Gi ∥ ∥O∥ ≤ hmax , by the parameter-shift Lipschitz bound of App. G with ∥Gi ∥ = 1; here hmax := maxj ∥hj ∥ is the maximum local-term operator norm defined in the Setting above; • |⟨O⟩(θ ′ )| ≤ ∥O∥ ≤ hmax for any θ ′ . 2 Hence |∂vi σ± | ≤ 2ε · hmax · hmax = 2εh2max . 2 3. Bounding ∥∇v σ∇,ℓ ∥. Differentiating eq. (H1) with respect to vi gives 2 ∂vi σ∇,ℓ =

2 2 (v) (v) + ∂vi σ− ∂vi σ+ . 4ε2 M

(H5)

2 Step (ii) showed |∂vi σ± | ≤ 2εh2max via eq. (H4). Substituting into eq. (H5): 2 |∂vi σ∇,ℓ |≤

2εh2max + 2εh2max h2max = . 4ε2 M εM

2 Crucially, each vi appears in the light cone of at most O(1) Hamiltonian terms by k-locality, so ∇v σ∇,ℓ has at most O(ξ(k, d)) non-zero components independently of J (the number of Hamiltonian terms) and hence of N . Summing the squared components over the non-zero indices:

2 ∥∇v σ∇,ℓ ∥2 =

X i : θi ∈LB

2 |∂vi σ∇,ℓ |2 ≤ ξ(k, d) ·

2 2 hmax . εM

Taking the square root defines p 2 Lσ := ∥∇v σ∇,ℓ ∥≤

which is independent of N at fixed (k, d).

ξ(k, d) h2max =O εM

(k + d)d h2max εM

! ,

(H6)

31 2 4. Concentration via Lemma 9. With Lσ from eq. (H6) established, σ∇,ℓ (v ℓ ) is Lσ -Lipschitz in the N i.i.d. ±1 2 2 2 Rademacher components of v ℓ , and Ev [σ∇,ℓ ] = σ̄∇ by definition. Applying Lemma 9 with f = σ∇,ℓ and L = Lσ : for every δ > 0, δ2 2 2 Pr |σ∇,ℓ − σ̄∇ | > δ ≤ 2 exp − 2 . (H7) 2Lσ p Equivalently, for any η ∈ (0, 1), choosing δ = Lσ 2 log(2/η): 2 2 Pr |σ∇,ℓ − σ̄∇ | ≤ Lσ

(H8)

p 2 log(2/η) ≥ 1 − η.

2 2 So σ∇,ℓ lies within O(Lσ ) of σ̄∇ with probability at least 1 − η, where Lσ = is independent of N .

ξ(k, d) h2max /(εM ) from eq. (H6)

Appendix I: Proof of Proposition 2: Rademacher minimises estimator variance

Proposition 2 claims that among all isotropic independent-component direction distributions, Rademacher uniquely minimises E∥b g V − g∥2 . To prove this we first derive a closed-form expression for the MSE as a function of the bV is unbiased (E[b component kurtosis κ, then minimise over κ. Since g g V ] = g), the MSE equals the variance: E∥b g V − g∥2 = E∥b g V ∥2 − ∥g∥2 . By i.i.d. sampling of the V directions,   " # X X ′ ′ 1 1 E∥b g V ∥2 = 2 E (g · v ℓ )2 ∥v ℓ ∥2 + 2 E  (g · v ℓ )(g · v ℓ ) (v ℓ · v ℓ ) . V V ′

Diagonal and off-diagonal decomposition.

ℓ

(I1)

ℓ̸=ℓ

′

For the off-diagonal terms (ℓ ̸= ℓ′ ), independence of the directions and isotropy (E[vv ⊤ ] = IN ) give E[(g·v ℓ )(g·v ℓ )(v ℓ · ′ v ℓ )] = (E[(g · v)v])2 = ∥g∥2 component-wise, so the off-diagonal sum contributes V (V − 1)∥g∥2 /V 2 = (V − 1)∥g∥2 /V . Diagonal term and kurtosis. For the diagonal terms, expand ∥v∥2 = E[(g · v)2 ∥v∥2 ] =

2 j vj and use component independence:

gi gj E[vi vj ∥v∥2 ]

i,j

gi2 E[vi2 vj2 ]

i,j

gi2 E[vi4 ] +

E[vi2 ]E[vj2 ]

j̸=i

gi2 (κ + N − 1) = ∥g∥2 (κ + N − 1),

(I2)

where the cross terms i ̸= j in the second line vanish because E[vi vj ] = 0 for i ̸= j (isotropy), and we used E[vi2 ] = 1 and E[vi4 ] = κ for i = j. The diagonal contribution to E∥b g V ∥2 is therefore ∥g∥2 (κ + N − 1)/V . Combining.

Adding diagonal and off-diagonal contributions and subtracting ∥g∥2 : ∥g∥2 (κ + N − 1) (V − 1)∥g∥2 + − ∥g∥2 V V ∥g∥2 (κ + N − 2) = , V

E∥b g V − g∥2 =

which is eq. (34).

(I3)

32 Minimisation over the direction distribution. Since κ = E[vi4 ] and E[vi2 ] = 1, Jensen’s inequality gives κ = E[vi4 ] ≥ (E[vi2 ])2 = 1, with equality iff vi2 is almost surely constant. Combined with E[vi2 ] = 1, this forces |vi | = 1 a.s., and combined with E[vi ] = 0 it forces vi ∈ {+1, −1} with equal probability (the Rademacher distribution). This uniquely minimises eq. (34) over all isotropic independent-component distributions with unit variance, giving the bound E∥b g V − g∥2 ≥ ∥g∥2 (N − 1)/V . 1.

Proof of Theorem 1: optimal allocation

Setup. The MSE of the Rademacher forward gradient estimator with V directions and M shots per direction (eq. (35)) is 2 (N − 1)∥g∥2 N σ̄∇ + , V VM

(I4)

MSE(V, M ) ≤ τ 2 ,

(I5)

MSE(V, M ) = and the cost-minimisation problem is min 2V M

V, M >0

s.t.

M ≥ Mmin .

The proof has five steps. Proof.

1. Eliminate V . The MSE constraint in eq. (I5) is equivalent to V ≥ ϕ(M ) :=

2 (N − 1)∥g∥2 + N σ̄∇ /M . 2 τ

(I6)

Since 2V M is strictly increasing in V at fixed M > 0, the infimum is attained at V = ϕ(M ). Substituting reduces eq. (I5) to 2 2 (N − 1)∥g∥2 M + N σ̄∇ min 2M ϕ(M ) = min . (I7) M ≥Mmin M ≥Mmin τ2 2. Minimise over M . The reduced objective eq. (I7) is affine in M with slope 2(N − 1)∥g∥2 /τ 2 > 0 (N ≥ 2 and ∥g∥2 > 0 by assumption), so its minimum on [Mmin , ∞) is (I8)

M ⋆ = Mmin . 3. Recover V ⋆ and C ⋆ . Substituting eq. (I8) into eq. (I6) gives V ⋆ = ϕ(Mmin ) =

2 /Mmin (N − 1)∥g∥2 + N σ̄∇ , τ2

(I9)

2 2 (N − 1)∥g∥2 Mmin + N σ̄∇ , τ2

(I10)

which is eq. (37), and the minimum cost is C ⋆ = 2V ⋆ Mmin = which is eq. (38). 4. Uniqueness. Both constraints in eq. (I5) are tight at (V ⋆ , M ⋆ ). Any feasible (V, M ) with M > Mmin has strictly larger cost (Step 2), and any with V > ϕ(M ) at the same M also has strictly larger cost (Step 1). Hence (V ⋆ , M ⋆ ) is the unique minimiser. 2 2 5. Reparameterisation. Set α := N σ̄∇ /(Mmin ∥g∥2 ) > 0, so Mmin = N σ̄∇ /(α∥g∥2 ). Substituting into eq. (I9):

V⋆ =

(N − 1)∥g∥2 + α∥g∥2 (N − 1 + α)∥g∥2 = , τ2 τ2

2 and M ⋆ = N σ̄∇ /(α∥g∥2 ). These are eq. (39).

(I11)

33 2.

Proof of Corollary 1: CRB-level optimality

The proof has two parts. First we establish the Cramér–Rao lower bound eq. (40) on the MSE of any unbiased estimator of g that queries the shot-noise oracle a total of B times. Second we show that the forward-gradient estimator at the optimal allocation of Theorem 1 attains this bound up to a constant that vanishes as N → ∞. Lower bound. We treat the oracle output at each query as a noisy linear measurement of the gradient. A query at shift sk returns a single shot with mean f (θ +sk ) ≈ f (θ)+s⊤ k g (local linearisation of f around θ, valid for the small 2 shifts ∥sk ∥ = O(ε) we use) and variance σm . Since the mean depends linearly on g and the variance is independent of g, the Fisher information of this query for the gradient parameter g is Ik =

sk s⊤ k . 2 σm

After B i.i.d. queries the total Fisher information is B

IB =

1 X sk s⊤ k, 2 σm k=1

b by the trace of the and the multi-parameter Cramér–Rao inequality bounds the MSE of any unbiased estimator g inverse Fisher information matrix, −1 E∥b g − g∥2 ≥ tr(IB ). 2 For Rademacher shifts sk ∈ {±ε}N drawn i.i.d. with E[sk s⊤ k ] = ε IN , the expected Fisher information is the diagonal 2 2 2 2 2 /(ε2 τ 2 ) queries are matrix E[IB ] = (Bε /σm ) IN , whose inverse trace is N σm /(Bε ). This is eq. (40): BCRB = N σm required to drive the MSE to a target τ 2 . 2 Achievability. Theorem 1 expresses the minimum cost C ⋆ in terms of the per-direction shot variance σ̄∇ . In the bias-neutral limit (small ε, finite-difference bias negligible), this variance reduces to 2 σ̄∇ =

2 σm 2M ε2

2 (the central-difference numerator has variance 2σm from two independent M -shot means, divided by (2ε)2 ). Substituting this expression into eq. (38) gives the QUIVER shot count

B⋆ =

2 2 2(N − 1 + α) N σ̄∇ (N − 1 + α) N σm = . 2 2 2 ατ ε ατ

Comparing to BCRB above, N −1+α B⋆ → 1 as N → ∞, = BCRB N for any fixed hyperparameter α > 0. quiver therefore saturates the Cramér–Rao lower bound up to a multiplicative factor that vanishes as 1/N .

Appendix J: iCANS/gCANS baselines and N -scaling of the VQE regime crossover

iCANS [23] and gCANS [24] are the closest parameter-shift-based shot-adaptive baselines to quiver: both adapt the per-parameter shot count via a gain-per-shot argument on the parameter-shift estimator. This appendix details the mechanism by which they collapse onto the floor smin at the conservative learning rates required at large N , reports the hyperparameter sweeps used to set the iCANS/gCANS baselines in Fig. 8 and Fig. 12, and extends the regime-crossover comparison of Fig. 8b across a broader N sweep.

(a) VQE TFIM, n = 8 5

Circuit depth d 20

(b) VQE TFIM, n = 10

Circuit depth d 10

200

250

300

Best energy error |E − E0 |

4 × 10 1 3 × 10 1 2 × 10 1

10 1

6 × 10 2 100 Exact GD

200

300

400

500

Number of parameters N iCANS (best of grid)

600

100

gCANS (best of grid)

150

Number of parameters N QUIVER fixed-M

QUIVER joint (V, M)

FIG. 12. quiver leads at all depths and budgets on the VQE N -scaling benchmark. Best energy error as a function of parameter count N for VQE TFIM at matched 5 × 106 -shot budget. (a) n = 8, N ∈ [80, 480]. (b) n = 10, N ∈ [80, 300]. For each method we plot the best error reached across a small hyperparameter grid at each N . Black: exact gradient descent. Orange: best iCANS. Green: best gCANS. Blue: best fixed-M quiver. Purple: best quiver joint (V, M ). On (b), the dotted vertical line marks the smallest depth at which exact gradient descent meets the target energy error δE = 0.1 on n = 10 TFIM (d = 8, N = 160); the next-smaller depth d = 5 (N = 100) fails to reach the target.

Mechanism.

The iCANS allocation for parameter i is MiiCANS =

σ b2 2Lη · 2i , 2 − Lη gbi

where σ bi2 , gbi2 are EMA estimates of the per-coordinate variance and squared gradient. For our VQE TFIM benchmark L = 15 (App. G); at η = 0.01 the Lη prefactor is 2Lη/(2 − Lη) ≈ 0.16, so to allocate above smin a parameter must have σ bi2 /b gi2 ≳ 12. Across all N ∈ {80, 160, 240, 320} we tested, this is rarely the case and the allocation is dominated by smin for almost every parameter. gCANS uses a global-norm variant of the same rule (one shot count per step rather than one per parameter), but inherits the same Lη prefactor and converges to a similar final error. Only at small N and aggressive η (e.g. η = 0.05 at N = 80) does the prefactor exceed unity and both methods adapt above smin ; in that regime training is unstable and the final error is worse. Sweep configurations. For the N -scaling comparison of Fig. 12 we ran iCANS at each N over four (η, sinit ) pairs: (0.003, 5), (0.003, 10), (0.01, 10), (0.05, 50), with smin = 2 throughout. For the regime-crossover figure (Fig. 8) we additionally swept sinit ∈ {2, 5, 10} at (N = 80, η = 0.05) for panel (a); panel (b) at (N = 160, η = 0.003) uses the single representative configuration (0.003, 5). In both figures we plot the best-performing iCANS and gCANS configuration at each N .

Appendix K: Circuit descriptions

The three benchmark circuit families used in this work are described below. Figure 13 shows representative n = 4 examples for each.

35 (a) |0⟩

|0⟩ |0⟩ |0⟩

(b) RBS

RBS RBS

RBS

|0⟩ |0⟩ |0⟩ |0⟩

y1 y2 y3

(c)

(d)

|0⟩

RZ •

• RY

RZ •

•

|0⟩

RZ • •

|0⟩

••

|0⟩

• • RY

••

|0⟩

bins:

V (x)

U (θ)

y0 y1 y2 y3

ℓ0 = y0 + y1 , ℓ1 = y2 + y3

RX (β0 ) (1)

RX (β1 ) RX (β2 )

RX (β0 ) (2)

RX (β3 )

RX (β1 ) RX (β2 ) RX (β3 )

FIG. 13. Circuit diagrams for the three benchmark families (n=4). (a) OQNN [12, 25, 26]: the loader (dashed) P encodes x into the unary state |ϕ(x)⟩ = i xi |ei ⟩; the variational pyramid (solid, N =n(n−1)/2 RBS gates) yields amplitudes yi = ⟨ei |ψ⟩. (b) Compact Pview (C=2): loader V (x) and variational U (θ) produce {yi }, which are partitioned into C contiguous bins giving logits ℓc = j∈binc yj . The original orthogonal-QNN proposal uses truncation (ℓc = yn−C+c ); binning uses the full distribution and reduces variance at large n. (c) HEA for VQE [70]: per-qubit RY RZ followed by CZ entanglers on adjacent pairs and the closing long-range CZ(0, n−1); N =2nd. (d) Multi-angle QAOA [71, 72]: H ⊗n initialises the uniform (l) superposition; each layer applies cost block UC (one global ZZ angle and n per-qubit Z angles) and per-qubit RX (βj,l ) mixers; N =(2n+1)d.

Orthogonal quantum neural network (OQNN)

The classification benchmarks (ECG5000 and MNIST) use a Hamming-weight-preserving orthogonal quantum neural network following [12, 25, 26]. An n-qubit circuit acts entirely within the single-excitation (Hamming-weight-1) subspace, spanned by the computational basis states {|ei ⟩}ni=1 where |ei ⟩ has a single 1 in position i. Within this subspace the reconfigurable beam-splitter (RBS) gate RBS(θ) on qubits (i, j) acts as a 2D Givens rotation, RBS(θ) |ei ⟩ = cos θ |ei ⟩ + sin θ |ej ⟩,

RBS(θ) |ej ⟩ = − sin θ |ei ⟩ + cos θ |ej ⟩,

and the full circuit of N RBS gates arranged in a triangular (staircase) layout implements an arbitrary n × n specialorthogonal matrix O(θ) ∈ SO(n) on the amplitude vector, with N = n(n − 1)/2 parameters. Classical input x ∈ Rn P P (unit-normalised) is encoded as |ϕ(x)⟩ = i xi |ei ⟩; the circuit maps this to |ψ(θ)⟩ = i yi (θ)|ei ⟩ where y = O(θ)x. The output n-vector y is the logit vector after the readout step below. Readout: binning vs. truncation. The original proposal for orthogonal QNNs [12, 25, 26] reads class scores by truncation: ℓc = yn−C+c , discarding the first n − C components of y. At the system sizes used here (n up to 60, C ∈ {5, 10}) this discards up to 98% of the output. We instead use a binning readout: partition the n outputs into C contiguous bins of size b = ⌊n/C⌋ (discarding at most C − 1 remainder components) and sum within each bin, ℓc =

cb X

yj ,

c = 1, . . . , C.

j=(c−1)b+1

Binning uses the full output vector and improves convergence at large n at no change in parameter count or gradient cost. Hardware-efficient ansatz for VQE

The VQE ground-state benchmarks use a hardware-efficient ansatz (HEA) [70] composed of d identical layers, each acting on n qubits. In layer k each qubit q ∈ {0, . . . , n − 1} receives a pair of single-qubit rotations RY (θq,k ) RZ (ϕq,k ), followed by a closed CZ ring: CZ gates on pairs (0, 1), (1, 2), . . . , (n − 2, n − 1) and a closing CZ(0, n − 1). The initial state is |0⟩⊗n . Total parameter count is N = 2nd. The CZ ring provides all-to-nearest-neighbour entanglement within

36 each layer at no additional parameter cost; for the 1D Ising Hamiltonian used here the ring topology matches the interaction structure of the Hamiltonian.

Multi-angle QAOA

The combinatorial optimisation benchmarks use a multi-angle QAOA (MA-QAOA) [71, 72] variant with d layers. Standard QAOA uses a single global cost angle γl and a single global mixer angle βl per layer l. MA-QAOA assigns independent parameters to each qubit and, in the cost layer, additionally includes per-qubit local-field terms, yielding greater expressivity at shallower depth than the standard QAOA ansatz: (l) UC = exp −iγzz,l

n−1 Y wij Zi Zj · exp −iγz,j,l Zj ,

(i,j)∈E (l)

UM =

n−1 Y

j=0

exp −iβj,l Xj ,

j=0

where the cost layer carries one global ZZ angle γzz,l shared across all edges, n per-qubit Z angles {γz,j,l }, and the mixer carries n independent RX angles {βj,l } per qubit. Total parameter count is N = (2n + 1)d; for n = 16 and ⊗n d = 3 this gives N = 99. The circuit is initialised P in the uniform superposition |+⟩ and the energy measured under 1 the weighted MaxCut Hamiltonian HC = − /2 (i,j)∈E wij (I − Zi Zj ).

Appendix L: Experimental hyperparameters

All ECG5000 and MNIST experiments use the following setup unless otherwise stated. Training uses Adam [27] with learning rate η = 0.001 and learning rate scheduling. Directional derivatives are estimated via central finite difference with step size ε = 0.1. Results are averaged over 3 random seeds. Shot budget. Forward gradient methods use a fixed budget of B = 1000 shots per gradient step, divided as B = 2V M where M is the number of shots per direction. For example, V = 10 gives M = 50. The parametershift rule uses M = 10 shots per parameter, giving B = 2N M total shots per step (which grows with N ). SPSA corresponds to V = 1, M = 500; RCD uses M = 500 for a single coordinate. ECG5000. 5-class classification (1 normal rhythm and 4 arrhythmia classes) on the UCR ECG5000 dataset. System sizes n ∈ {10, 20, 30, 40, 50, 60} with N = n(n − 1)/2 parameters. Forward methods train for 10,000 epochs; exact gradient and PS for 2,000 epochs. Shot budget termination at 5 × 108 total shots. MNIST. 10-class classification on downsampled MNIST. System sizes n ∈ {10, 20, 40, 50}. Forward and SPSA/RCD methods train for 5,000 epochs; exact gradient and PS for 1,000 epochs. Same shot budget and Adam learning rate as ECG5000. VQE. TFIM Hamiltonian with J = h = 1, open boundary conditions. Hardware-efficient ansatz with nL layers, N = 2nnL parameters. The hero and regime-crossover figures use n = 10; the N -scaling, trajectory, and ε-sweep figures use n = 8. Forward gradient learning rate η = 0.003; PS learning rate η = 0.01. Parameter initialisation scale 0.1. Energy error is reported as |Ebest − Eexact | where Eexact is computed by exact diagonalisation of the Hamiltonian. QAOA MaxCut. Weighted Erdős–Rényi graphs on n = 16 vertices with edge probability p = 0.5 and weights drawn uniformly from [0.1, 1.0]. Depth d = 3 ansatz with N = (n+1)d + nd = (2n+1)d = 99 parameters. Forward gradient learning rate η = 0.1; ε = 0.1. Ground energy Eopt computed by brute force. Approximation ratio reported as r = Ebest /Eopt .

37 Adaptive V ⋆ and quiver parameters. The quiver optimiser of Section VII C uses exponential moving 2 averages with decay β = 0.9 for both gbt2 and σ bt2 . The sample estimator for σ̄∇ at each epoch is the empirical variance M V e vℓ L } of the V directional derivative estimates {∇ around their mean, which under the noise-concentration result ℓ=1 2 of Section VII B converges to σ̄∇ /M . The per-direction shot count M is the only hyperparameter of the optimiser; clamp bounds are Vmin = 2 and Vmax = N unless stated otherwise. All runs with quiver use the same Adam learning rate as the corresponding fixed-V baseline in Fig. 2.

Record · ID 267649 · SHA-256 0ff8339fc82eab06
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.