Active Quantum Kernel Acquisition for Gaussian Process Regression
Summary
This paper proposes active shot allocation strategies for quantum kernel estimation in Gaussian process regression, deriving pair-level sensitivities to guide non-uniform shot budgets and demonstrating significant improvements in test RMSE over uniform allocation on benchmarks.
View Cached Full Text
Cached at: 06/30/26, 05:29 AM
# Active Quantum Kernel Acquisition for Gaussian Process Regression
Source: [https://arxiv.org/html/2606.28833](https://arxiv.org/html/2606.28833)
###### Abstract
Quantum kernel estimation on near\-term hardware is shot\-budgeted: every entry of the kernel Gram matrix is a Bernoulli expectation that must be sampled with a finite number of circuit executions\. Recent work on quantum kernel classification has shown that allocating shots non\-uniformly across kernel entries, weighted by their downstream task sensitivity, can reduce the shot budget required to reach a target accuracy\. We extend this idea to Gaussian process \(GP\) regression, a setting whose downstream quantities \(full\-spectrum posterior variance, log\-determinant, marginal likelihood\) couple to kernel error more tightly than the sign\-only outputs of classification\. We derive three closed\-form pair\-level sensitivities — predictive coupling\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|, leave\-one\-out residual, and marginal\-likelihood gradient — and plug them into a Neyman\-style minimum\-variance allocation rule\. To prevent catastrophic over\-concentration when the warm\-up sensitivity estimate is itself noisy, we add a high uniform coverage floor justified by a Frobenius lower bound on the missing\-entry perturbation\. On four UCI benchmarks and two synthetic RBF \+ Bernoulli controlled studies, the resulting allocator delivers1010–21%21\\%test\-RMSE improvement over uniform allocation across the moderate\-budget regime\. The gain transfers \(i\) to genuine ZZ and Pauli\-Z quantum kernels on quantum\-natural data \(−13\-13–15%15\\%at low budget,p<0\.05p<0\.05paired\) and \(ii\) to four downstream tasks \(Bayesian quadrature, heteroscedastic regression, hyperparameter learning, multi\-output Cokriging\)\. On UCI features embedded into a ZZ kernel the gain disappears, consistent with the exponential\-concentration regime where shot allocation has nothing to exploit\. The three pair\-level sensitivities we propose are not arbitrary Fisher surrogates:\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|is the rank\-1 specialization \(atM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}\) of the exact Neyman weight under the predictive\-MSE objective; the marginal\-likelihood gradient is the exact partial derivative; the LOO sensitivity drops a controllable correction\.
## Introduction
Quantum kernel methods promise computationally efficient access to feature spaces that are conjectured to be classically intractable for certain structured data\(Havlíčeket al\.[2019](https://arxiv.org/html/2606.28833#bib.bib4); Schuld and Killoran[2019](https://arxiv.org/html/2606.28833#bib.bib13)\)\. The bottleneck in practical deployment is that each kernel entryKij=\|⟨ϕ\(xi\)\|ϕ\(xj\)⟩\|2K\_\{ij\}=\|\\langle\\phi\(x\_\{i\}\)\|\\phi\(x\_\{j\}\)\\rangle\|^\{2\}is the success probability of an inversion\-test circuit and must be estimated from a finite number of shotssijs\_\{ij\}\. The resulting Bernoulli noise scales asVar\(K^ij\)=Kij\(1−Kij\)/sij\\mathrm\{Var\}\(\\widehat\{K\}\_\{ij\}\)=K\_\{ij\}\(1\-K\_\{ij\}\)/s\_\{ij\}, and the total cost of building ann×nn\\times nGram matrix to a uniform precision grows asO\(n2/ϵ2\)O\(n^\{2\}/\\epsilon^\{2\}\), which already strains current hardware budgets fornnin the low hundreds\(Huanget al\.[2021](https://arxiv.org/html/2606.28833#bib.bib6); Thanasilpet al\.[2024](https://arxiv.org/html/2606.28833#bib.bib16)\)\.
A natural way around this scaling is to recognize that not all kernel entries matter equally to the downstream task\. For quantum kernel classification with kernel ridge regression \(KRR\) or SVM, recent work showed that the loss gradient∂ℒ/∂Kij\\partial\\mathcal\{L\}/\\partial K\_\{ij\}has highly heterogeneous magnitudes across pairs, and that allocating shots proportional to the Neyman\-optimal targetsij∗∝\|∂ℒ/∂Kij\|Kij\(1−Kij\)s^\{\*\}\_\{ij\}\\propto\|\\partial\\mathcal\{L\}/\\partial K\_\{ij\}\|\\,\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}yields substantial accuracy gains over uniform allocation at the same budget\(Xuet al\.[2026](https://arxiv.org/html/2606.28833#bib.bib1); Miroszewski[2026](https://arxiv.org/html/2606.28833#bib.bib7)\)\. The argument is a direct application of Neyman’s19341934stratified\-sampling result\(Neyman[1992](https://arxiv.org/html/2606.28833#bib.bib9); Pukelsheim[2006](https://arxiv.org/html/2606.28833#bib.bib10)\)cast in the language of quantum kernel estimation\.
This paper asks: does the same idea help in Gaussian process \(GP\) regression? The answer is not obvious\. On the one hand, GP regression has even more structure to exploit — a closed\-form posterior, a closed\-form marginal likelihood, leave\-one\-out residuals in closed form\. On the other hand, GP predictions and likelihoods both depend on the inverse\(𝐊\+σn2I\)−1\(\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}, and inversion amplifies kernel noise by a factor that scales like the condition number1/σn21/\\sigma\_\{n\}^\{2\}\. A shot\-allocation policy that is mildly suboptimal for classification could be catastrophic for GP regression\.
We make four contributions:
1. 1\.We derive three closed\-form pair\-level sensitivities for GP regression — the predictive coupling\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|, the leave\-one\-out residual, and the marginal\-likelihood gradient — each capturing a different facet of how kernel uncertainty propagates to predictions\.
2. 2\.We identify a failure mode unique to GP regression: when the per\-pair budget is low, the warm\-up kernel estimate is too noisy to produce a reliable sensitivity, and sensitivity\-driven allocation*worsens*test error by hundreds of percent\. We show that a50%50\\%uniform\-coverage floor restores robustness, while the corresponding classification setup needs only1010–20%20\\%\.
3. 3\.Empirically, on four UCI regression benchmarks and two synthetic settings \(ntr=200n\_\{tr\}=200,55–1010seeds\), the allocator delivers consistent1010–21%21\\%RMSE reduction in the moderate\-budget regime \(∼50\\sim 50–250250shots per pair\)\.
4. 4\.Unlike in classification, where allocation gains are largest on*sparse*, anchor\-heavy data, GP regression shows the largest gains on*dense*GP draws\. We explain this via the structure of\(𝐊\+σn2I\)−1\(\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}: even when the kernel is homogeneous, the spectral decay creates predictive heterogeneity, and the same heterogeneity drives sensitivity\.
## Background
### Quantum Kernel Estimation as Bernoulli Sampling
A quantum feature mapϕ:𝒳→ℋQ\\phi:\\mathcal\{X\}\\to\\mathcal\{H\}\_\{Q\}is realized by a parametrized unitaryU\(x\)U\(x\)acting on the all\-zero state\|0⟩⊗q\|0\\rangle^\{\\otimes q\}\. The induced kernel is the fidelity
K\(xi,xj\)=\|⟨0\|U\(xj\)†U\(xi\)\|0⟩\|2∈\[0,1\]\.K\(x\_\{i\},x\_\{j\}\)\\;=\\;\|\\langle 0\|U\(x\_\{j\}\)^\{\\dagger\}U\(x\_\{i\}\)\|0\\rangle\|^\{2\}\\;\\in\\;\[0,1\]\.\(1\)The standard*inversion test*\(Havlíčeket al\.[2019](https://arxiv.org/html/2606.28833#bib.bib4)\)estimatesK\(xi,xj\)K\(x\_\{i\},x\_\{j\}\)by preparingU\(xj\)†U\(xi\)\|0⟩U\(x\_\{j\}\)^\{\\dagger\}U\(x\_\{i\}\)\|0\\rangle, measuring in the computational basis, and recording the probability of the all\-zero outcome\. Withsijs\_\{ij\}independent circuit executions \(“shots”\), the empirical estimatorK^ij\\widehat\{K\}\_\{ij\}is a sample mean ofsijs\_\{ij\}Bernoulli\(Kij\)\(K\_\{ij\}\)random variables and satisfies
𝔼\[K^ij\]=Kij,Var\(K^ij\)=Kij\(1−Kij\)sij\.\\mathbb\{E\}\[\\widehat\{K\}\_\{ij\}\]=K\_\{ij\},\\qquad\\mathrm\{Var\}\(\\widehat\{K\}\_\{ij\}\)=\\frac\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\{s\_\{ij\}\}\.\(2\)We treat the per\-entry shot budgetsijs\_\{ij\}as a design variable subject to the global constraint∑i≤jsij=B\\sum\_\{i\\leq j\}s\_\{ij\}=B\. This is the only resource we will allocate; all other choices \(feature map, circuit depth, error mitigation\) are taken as fixed\.
#### Exponential concentration regime\.
A separate failure mode arises when the quantum kernelK\(xi,xj\)K\(x\_\{i\},x\_\{j\}\)itself concentrates around a fixed value as the system grows\(Thanasilpet al\.[2024](https://arxiv.org/html/2606.28833#bib.bib16)\)\. In that regime, the signal\-to\-noise ratio\|Kij−Ki′j′\|2/Var\(K^\)\|K\_\{ij\}\-K\_\{i^\{\\prime\}j^\{\\prime\}\}\|^\{2\}/\\mathrm\{Var\}\(\\widehat\{K\}\)vanishes exponentially inqq, and shot reallocation cannot rescue it\. We work outside this regime, treatingK\(⋅,⋅\)K\(\\cdot,\\cdot\)as a fixed bounded kernel and asking only how to spend a budgetBBoptimally\.
### Gaussian Process Regression
A Gaussian process priorf∼𝒢𝒫\(0,k\(⋅,⋅\)\)f\\sim\\mathcal\{GP\}\(0,k\(\\cdot,\\cdot\)\)with observations
yi=f\(xi\)\+εi,εi∼iid𝒩\(0,σn2\),y\_\{i\}=f\(x\_\{i\}\)\+\\varepsilon\_\{i\},\\qquad\\varepsilon\_\{i\}\\stackrel\{\{\\scriptstyle\\text\{iid\}\}\}\{\{\\sim\}\}\\mathcal\{N\}\(0,\\sigma\_\{n\}^\{2\}\),\(3\)admits the closed\-form posteriorf\|𝐲∼𝒢𝒫\(μ∗,k∗\)f\|\\mathbf\{y\}\\sim\\mathcal\{GP\}\(\\mu\_\{\*\},k\_\{\*\}\)with predictive mean and variance at anyx∗x\_\{\*\}:
μ∗\(x∗\)\\displaystyle\\mu\_\{\*\}\(x\_\{\*\}\)=𝐤∗⊤𝜶,𝜶:=A−1𝐲,\\displaystyle=\\mathbf\{k\}\_\{\*\}^\{\\top\}\\bm\{\\alpha\},\\qquad\\bm\{\\alpha\}:=A^\{\-1\}\\mathbf\{y\},\(4\)σ∗2\(x∗\)\\displaystyle\\sigma\_\{\*\}^\{2\}\(x\_\{\*\}\)=k\(x∗,x∗\)−𝐤∗⊤A−1𝐤∗\+σn2,\\displaystyle=k\(x\_\{\*\},x\_\{\*\}\)\-\\mathbf\{k\}\_\{\*\}^\{\\top\}A^\{\-1\}\\mathbf\{k\}\_\{\*\}\+\\sigma\_\{n\}^\{2\},\(5\)whereA:=𝐊\+σn2IA:=\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}Iand𝐤∗=\(k\(x∗,xi\)\)i\\mathbf\{k\}\_\{\*\}=\(k\(x\_\{\*\},x\_\{i\}\)\)\_\{i\}\. The \(negative\) log marginal likelihood of the training labels is
ℒ\(𝐊;𝐲\)=12𝐲⊤A−1𝐲\+12log\|A\|\+n2log\(2π\)\.\\mathcal\{L\}\(\\mathbf\{K\};\\mathbf\{y\}\)\\;=\\;\\tfrac\{1\}\{2\}\\mathbf\{y\}^\{\\top\}A^\{\-1\}\\mathbf\{y\}\+\\tfrac\{1\}\{2\}\\log\|A\|\+\\tfrac\{n\}\{2\}\\log\(2\\pi\)\.\(6\)When𝐊\\mathbf\{K\}is estimated with shot noise, all three quantitiesμ∗\\mu\_\{\*\},σ∗2\\sigma\_\{\*\}^\{2\},ℒ\\mathcal\{L\}inherit the noise through the inverseA−1A^\{\-1\}, which is the source of the amplification we exploit in Section[Method: AQKA\-GP](https://arxiv.org/html/2606.28833#Sx3)\.
#### What makes GP shot\-hungrier than classification?
The KRR / SVM losses used in classification AQKA\(Xuet al\.[2026](https://arxiv.org/html/2606.28833#bib.bib1); Miroszewski[2026](https://arxiv.org/html/2606.28833#bib.bib7)\)also involve𝐊−1\\mathbf\{K\}^\{\-1\}: KRR’s𝜶=\(𝐊\+λI\)−1𝐲\\bm\{\\alpha\}=\(\\mathbf\{K\}\+\\lambda I\)^\{\-1\}\\mathbf\{y\}usesλ\\lambdain the same roleσn2\\sigma\_\{n\}^\{2\}plays for GP\. So “GP has an inverse, classification does not” is*not*the structural difference\. The difference is which downstream quantities the user cares about\. Classification’s0/10/1accuracy depends onsign\(μ∗\)\\mathrm\{sign\}\(\\mu\_\{\*\}\), which is robust to multiplicative scaling of𝜶\\bm\{\\alpha\}; GP cares about
1. 1\.the*magnitude*ofμ∗\\mu\_\{\*\}\(test RMSE/NLL\),
2. 2\.the predictive varianceσ∗2=k\(x∗,x∗\)−𝐤∗⊤A−1𝐤∗\+σn2\\sigma\_\{\*\}^\{2\}=k\(x\_\{\*\},x\_\{\*\}\)\-\\mathbf\{k\}\_\{\*\}^\{\\top\}A^\{\-1\}\\mathbf\{k\}\_\{\*\}\+\\sigma\_\{n\}^\{2\}, which depends on the full spectrum ofA−1A^\{\-1\}rather than its top eigendirection,
3. 3\.the marginal\-likelihoodℒ\\mathcal\{L\}, which adds alog\|A\|\\log\|A\|term that aggregates errors over*every*eigenvalue\.
Items \(ii\) and \(iii\) are GP\-specific and have no analogue in0/10/1classification\. They are why GP needs more shots per pair than classification, and they motivate the three GP\-specific sensitivities of Section[Method: AQKA\-GP](https://arxiv.org/html/2606.28833#Sx3)\.
#### Conditioning ofA−1A^\{\-1\}\.
The Gram matrix𝐊\\mathbf\{K\}has eigenvalues0≤λn≤⋯≤λ10\\leq\\lambda\_\{n\}\\leq\\dots\\leq\\lambda\_\{1\}, withλ1=O\(n\)\\lambda\_\{1\}=O\(n\)for a typical bounded kernel andλn\\lambda\_\{n\}that can be arbitrarily small\. The corresponding eigenvalues ofA−1A^\{\-1\}are1/\(λi\+σn2\)∈\[1/\(λ1\+σn2\),1/σn2\]1/\(\\lambda\_\{i\}\+\\sigma\_\{n\}^\{2\}\)\\in\[1/\(\\lambda\_\{1\}\+\\sigma\_\{n\}^\{2\}\),1/\\sigma\_\{n\}^\{2\}\]\. A kernel perturbationΔ𝐊\\Delta\\mathbf\{K\}with operator norm‖Δ𝐊‖op<σn2\\\|\\Delta\\mathbf\{K\}\\\|\_\{\\mathrm\{op\}\}<\\sigma\_\{n\}^\{2\}produces \(by the resolvent identity\)
‖A^−1−A−1‖op≤‖Δ𝐊‖opσn2\(σn2−‖Δ𝐊‖op\)\.\\\|\\widehat\{A\}^\{\-1\}\-A^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\;\\leq\\;\\frac\{\\\|\\Delta\\mathbf\{K\}\\\|\_\{\\mathrm\{op\}\}\}\{\\sigma\_\{n\}^\{2\}\(\\sigma\_\{n\}^\{2\}\-\\\|\\Delta\\mathbf\{K\}\\\|\_\{\\mathrm\{op\}\}\)\}\.\(7\)The kernel error is amplified by up to1/σn41/\\sigma\_\{n\}^\{4\}\. The KRR predictor inherits the same amplification withλ\\lambdain place ofσn2\\sigma\_\{n\}^\{2\}, but classification accuracy averages over signs and is less sensitive to it; GP test NLL has no such averaging\.
### Neyman Minimum\-Variance Shot Allocation
For any differentiable lossℒ\(𝐊\)\\mathcal\{L\}\(\\mathbf\{K\}\), a second\-order Taylor expansion around the noise\-free kernel gives
ℒ\(𝐊^\)−ℒ\(𝐊\)\\displaystyle\\mathcal\{L\}\(\\widehat\{\\mathbf\{K\}\}\)\-\\mathcal\{L\}\(\\mathbf\{K\}\)=∑i≤jgijΔij\\displaystyle=\\sum\_\{i\\leq j\}g\_\{ij\}\\,\\Delta\_\{ij\}\+12∑i≤j,k≤ℓHij,kℓΔijΔkℓ\+o\(‖Δ‖2\),\\displaystyle\\quad\+\\tfrac\{1\}\{2\}\\sum\_\{i\\leq j,\\,k\\leq\\ell\}H\_\{ij,k\\ell\}\\,\\Delta\_\{ij\}\\Delta\_\{k\\ell\}\+o\(\\\|\\Delta\\\|^\{2\}\),\(8\)wheregij=∂ℒ/∂Kijg\_\{ij\}=\\partial\\mathcal\{L\}/\\partial K\_\{ij\},HHis the Hessian, andΔij:=K^ij−Kij\\Delta\_\{ij\}:=\\widehat\{K\}\_\{ij\}\-K\_\{ij\}is the shot noise\. The shot noise has zero mean and diagonal covariance \(different pairs are independent\), so taking expectations cancels the linear term and gives
𝔼\[ℒ\(𝐊^\)−ℒ\(𝐊\)\]=12∑i≤jHij,ijKij\(1−Kij\)sij\+o\(B−1\)\.\\mathbb\{E\}\\bigl\[\\mathcal\{L\}\(\\widehat\{\\mathbf\{K\}\}\)\-\\mathcal\{L\}\(\\mathbf\{K\}\)\\bigr\]=\\tfrac\{1\}\{2\}\\sum\_\{i\\leq j\}H\_\{ij,ij\}\\,\\frac\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\{s\_\{ij\}\}\+o\(B^\{\-1\}\)\.\(9\)Replacing the diagonal Hessian by its Fisher\-information surrogateHij,ij≈gij2H\_\{ij,ij\}\\approx g\_\{ij\}^\{2\}and minimizing under the budget constraint∑sij=B\\sum s\_\{ij\}=Bvia Lagrange multipliers yields the*Neyman minimum\-variance allocation*:
sij∗∝\|gij\|Kij\(1−Kij\),∑i≤jsij∗=B\.s^\{\*\}\_\{ij\}\\;\\propto\\;\|g\_\{ij\}\|\\,\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\},\\qquad\\sum\_\{i\\leq j\}s^\{\*\}\_\{ij\}=B\.\(10\)This is precisely Neyman’s 1934 stratified\-sampling rule\(Neyman[1992](https://arxiv.org/html/2606.28833#bib.bib9); Pukelsheim[2006](https://arxiv.org/html/2606.28833#bib.bib10)\)applied at the granularity of kernel entries\. The classification AQKA framework\(Xuet al\.[2026](https://arxiv.org/html/2606.28833#bib.bib1)\)uses this rule withgijg\_\{ij\}taken from the KRR or SVM training loss; the concurrent work ofMiroszewski \([2026](https://arxiv.org/html/2606.28833#bib.bib7)\)uses it for kernelized SVMs under noisy observations\. Our contribution begins where these stop: we instantiategijg\_\{ij\}for three GP\-regression objectives whose closed forms involve the kernel inverse\.
## Method: AQKA\-GP
We instantiate the Neyman rule \([10](https://arxiv.org/html/2606.28833#Sx2.E10)\) for three GP objectives\. Each derivation uses the same identity for the derivative of the inverse: forA=𝐊\+σn2IA=\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I,
∂A−1/∂Kij=−A−1\(eiej⊤\+ejei⊤\)A−1⋅12\(1\+δij\),\\partial A^\{\-1\}/\\partial K\_\{ij\}\\;=\\;\-A^\{\-1\}\\,\(e\_\{i\}e\_\{j\}^\{\\top\}\+e\_\{j\}e\_\{i\}^\{\\top\}\)\\,A^\{\-1\}\\;\\cdot\\;\\tfrac\{1\}\{2\}\(1\+\\delta\_\{ij\}\),\(11\)whereeie\_\{i\}is theiith standard basis vector and the symmetric factor accounts for the off\-diagonal/diagonal distinction\. We drop the factor12\(1\+δij\)\\tfrac\{1\}\{2\}\(1\+\\delta\_\{ij\}\)in the sensitivity \(it is the same constant across off\-diagonal pairs and is absorbed into the budget\)\.
### Three Closed\-Form Sensitivities
#### \(i\) Predictive coupling\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|for test RMSE\.
The mean\-squared test error isℒrmse=1n∗∑t\(yt∗−μ∗\(xt∗\)\)2\\mathcal\{L\}^\{\\text\{rmse\}\}=\\tfrac\{1\}\{n\_\{\*\}\}\\sum\_\{t\}\(y\_\{t\}^\{\*\}\-\\mu\_\{\*\}\(x\_\{t\}^\{\*\}\)\)^\{2\}, whereμ∗\(xt∗\)=𝐤∗t⊤𝜶\\mu\_\{\*\}\(x\_\{t\}^\{\*\}\)=\\mathbf\{k\}\_\{\*t\}^\{\\top\}\\bm\{\\alpha\}with𝜶=A−1𝐲\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}\. Differentiatingμ∗\\mu\_\{\*\}via the inverse identity gives
∂μ∗\(xt∗\)∂Kij\\displaystyle\\frac\{\\partial\\mu\_\{\*\}\(x\_\{t\}^\{\*\}\)\}\{\\partial K\_\{ij\}\}=−𝐤∗t⊤A−1\(eiej⊤\+ejei⊤\)𝜶\\displaystyle=\-\\mathbf\{k\}\_\{\*t\}^\{\\top\}A^\{\-1\}\\bigl\(e\_\{i\}e\_\{j\}^\{\\top\}\+e\_\{j\}e\_\{i\}^\{\\top\}\\bigr\)\\bm\{\\alpha\}=−\(btiαj\+btjαi\),\\displaystyle=\-\\bigl\(b\_\{ti\}\\alpha\_\{j\}\+b\_\{tj\}\\alpha\_\{i\}\\bigr\),\(12\)withbti:=\[𝐤∗t⊤A−1\]ib\_\{ti\}:=\[\\mathbf\{k\}\_\{\*t\}^\{\\top\}A^\{\-1\}\]\_\{i\}\. Squaring and summing overtt, the cross terms average out under test data drawn independently of𝜶\\bm\{\\alpha\}, leaving the leading term proportional toαi2αj2\\alpha\_\{i\}^\{2\}\\alpha\_\{j\}^\{2\}\. We therefore use
sensijpred=\|αiαj\|,𝜶=A−1𝐲\.\\mathrm\{sens\}^\{\\mathrm\{pred\}\}\_\{ij\}\\;=\\;\|\\alpha\_\{i\}\\alpha\_\{j\}\|,\\qquad\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}\.\(13\)This is the direct GP analogue of the classification sensitivity inXuet al\.\([2026](https://arxiv.org/html/2606.28833#bib.bib1)\): in both cases the rank\-1 outer product𝜶𝜶⊤\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}captures how much each pair couples labels to predictions\.
#### \(ii\) Marginal\-likelihood gradient for hyperparameter learning\.
Differentiatingℒ\(𝐊;𝐲\)\\mathcal\{L\}\(\\mathbf\{K\};\\mathbf\{y\}\)in \([6](https://arxiv.org/html/2606.28833#Sx2.E6)\) withA=𝐊\+σn2IA=\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}Iuses two standard identities:∂log\|A\|/∂Kij=\[A−1\]ij\\partial\\log\|A\|/\\partial K\_\{ij\}=\[A^\{\-1\}\]\_\{ij\}and∂\(𝐲⊤A−1𝐲\)/∂Kij=−αiαj\\partial\(\\mathbf\{y\}^\{\\top\}A^\{\-1\}\\mathbf\{y\}\)/\\partial K\_\{ij\}=\-\\alpha\_\{i\}\\alpha\_\{j\}\. Adding the two contributions,
gijmarg=∂ℒ∂Kij=12\[A−1\]ij−12αiαj\.g\_\{ij\}^\{\\mathrm\{marg\}\}\\;=\\;\\frac\{\\partial\\mathcal\{L\}\}\{\\partial K\_\{ij\}\}\\;=\\;\\tfrac\{1\}\{2\}\\bigl\[A^\{\-1\}\\bigr\]\_\{ij\}\\;\-\\;\\tfrac\{1\}\{2\}\\alpha\_\{i\}\\alpha\_\{j\}\.\(14\)This is the sensitivity to use when the downstream objective is marginal\-likelihood maximization \(kernel\-width or noise\-level learning\)\. Crucially, the two terms have opposite signs and partially cancel at the marginal\-likelihood optimum, sogijmargg\_\{ij\}^\{\\mathrm\{marg\}\}shrinks near convergence while remaining large far from it\.111Naively this looks like a drawback \(no signal at the optimum\), but the regime we care about — iterative kernel learning under shot budget — never starts at the optimum\.We use
sensijmarg=\|12\[A−1\]ij−12αiαj\|\.\\mathrm\{sens\}^\{\\mathrm\{marg\}\}\_\{ij\}\\;=\\;\\bigl\|\\tfrac\{1\}\{2\}\[A^\{\-1\}\]\_\{ij\}\-\\tfrac\{1\}\{2\}\\alpha\_\{i\}\\alpha\_\{j\}\\bigr\|\.\(15\)
#### \(iii\) Leave\-one\-out residual for predictive calibration\.
The closed\-form GP LOO predictive mean is\(Rasmussen[2003](https://arxiv.org/html/2606.28833#bib.bib11)\)
μ−i\(xi\)=yi−αi/\[A−1\]ii,\\mu\_\{\-i\}\(x\_\{i\}\)=y\_\{i\}\-\\alpha\_\{i\}/\[A^\{\-1\}\]\_\{ii\},\(16\)so the LOO residual isei:=αi/\[A−1\]iie\_\{i\}:=\\alpha\_\{i\}/\[A^\{\-1\}\]\_\{ii\}\. The LOO loss isℒloo=∑iei2\\mathcal\{L\}^\{\\mathrm\{loo\}\}=\\sum\_\{i\}e\_\{i\}^\{2\}\. The gradient with respect toKijK\_\{ij\}is, by the chain rule and the inverse identity,
∂ei/∂Kij∝\[A−1\]ij/\[A−1\]ii,\\partial e\_\{i\}/\\partial K\_\{ij\}\\;\\propto\\;\[A^\{\-1\}\]\_\{ij\}/\[A^\{\-1\}\]\_\{ii\},\(17\)so∂ℒloo/∂Kij∝ei\[A−1\]ij\+ej\[A−1\]ji\\partial\\mathcal\{L\}^\{\\mathrm\{loo\}\}/\\partial K\_\{ij\}\\propto e\_\{i\}\[A^\{\-1\}\]\_\{ij\}\+e\_\{j\}\[A^\{\-1\}\]\_\{ji\}\. Symmetrizing and taking absolute values,
sensijloo=\|ei\|⋅\|\[A−1\]ij\|\+\|ej\|⋅\|\[A−1\]ji\|\.\\mathrm\{sens\}^\{\\mathrm\{loo\}\}\_\{ij\}\\;=\\;\|e\_\{i\}\|\\cdot\|\[A^\{\-1\}\]\_\{ij\}\|\\;\+\\;\|e\_\{j\}\|\\cdot\|\[A^\{\-1\}\]\_\{ji\}\|\.\(18\)This sensitivity emphasizes pairs where one endpoint has a large LOO residual and is also tightly coupled to its neighbours throughA−1A^\{\-1\}\.
#### Computational cost\.
All three sensitivities reduce to invertingAAonce per allocation round \(O\(n3\)O\(n^\{3\}\)\) and forming the relevant outer products\. Forn=200n=200the cost is sub\-second on a laptop and is negligible compared to circuit submission\. The sameA−1A^\{\-1\}can be reused across rounds within a single training loop\.
### The High\-Floor Allocator
The naive algorithm — \(i\) warm\-up with a small fraction ofBBvia random pair sampling, \(ii\) compute sensitivity from𝐊^warm\\widehat\{\\mathbf\{K\}\}\_\{\\mathrm\{warm\}\}, \(iii\) distribute the rest by \([10](https://arxiv.org/html/2606.28833#Sx2.E10)\) — works in classification but fails on GP\. We trace the failure to a single mechanism:
> *Unsampled pairs default toK^ij=0\.5\\widehat\{K\}\_\{ij\}=0\.5\(the Bernoulli prior mean\)\. For classification with𝐊\\mathbf\{K\}\-product accuracy metrics, this is recoverable: a few uniform\-fill iterations restore prediction quality\. For GP regression, even a small fraction ofKij=0\.5K\_\{ij\}=0\.5entries in the Gram matrix is amplified by\(𝐊\+σn2I\)−1\(\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}into a catastrophically wrong𝛂\\bm\{\\alpha\}, because the inverse couples everyKijK\_\{ij\}to every prediction\.*
The fix is a high uniform coverage floor\. Let0<ρ<10<\\rho<1denote the fraction ofBBreserved for uniform sampling\. We useρ=0\.5\\rho=0\.5throughout, contrasting with theρ≤0\.2\\rho\\leq 0\.2that suffices for classification\.
Algorithm 1AQKA\-GP shot allocation \(full\-GP, fixed\-budget\)0:shot budget
BB, training set
\(X,𝐲\)\(X,\\mathbf\{y\}\), observation noise
σn2\\sigma\_\{n\}^\{2\}, sensitivity choice
sens∈\{senspred,sensmarg,sensloo\}\\mathrm\{sens\}\\in\\\{\\mathrm\{sens\}^\{\\mathrm\{pred\}\},\\mathrm\{sens\}^\{\\mathrm\{marg\}\},\\mathrm\{sens\}^\{\\mathrm\{loo\}\}\\\}, floor
ρ∈\[0,1\]\\rho\\in\[0,1\], warm
ρw∈\[0,1\]\\rho\_\{w\}\\in\[0,1\]
1:\{*Stage 1: warm\-up \+ uniform coverage floor\.*\}
2:Sample
⌊ρwB⌋\\lfloor\\rho\_\{w\}B\\rfloorshots uniformly at random over
n\(n\+1\)/2n\(n\+1\)/2pairs
→\\toinitial
𝐊^warm\\widehat\{\\mathbf\{K\}\}\_\{\\mathrm\{warm\}\}
3:Distribute
⌊ρB⌋\\lfloor\\rho B\\rfloorshots equally over all pairs \(uniform floor; Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)\)
4:Update
𝐊^\\widehat\{\\mathbf\{K\}\}from accumulated counts:
K^ij=countsij/shotsij\\widehat\{K\}\_\{ij\}=\\mathrm\{counts\}\_\{ij\}/\\mathrm\{shots\}\_\{ij\}
5:\{*Stage 2: sensitivity\-weighted top\-up\.*\}
6:Form sensitivity field
Sij←sens\(𝐊^,𝐲,σn2\)S\_\{ij\}\\leftarrow\\mathrm\{sens\}\(\\widehat\{\\mathbf\{K\}\},\\mathbf\{y\},\\sigma\_\{n\}^\{2\}\)//O\(n3\)O\(n^\{3\}\)
7:Form Neyman weights
wij←SijK^ij\(1−K^ij\)w\_\{ij\}\\leftarrow S\_\{ij\}\\sqrt\{\\widehat\{K\}\_\{ij\}\(1\-\\widehat\{K\}\_\{ij\}\)\}
8:
Brem←B−shots\_usedB\_\{\\mathrm\{rem\}\}\\leftarrow B\-\\mathrm\{shots\\\_used\}
9:Allocate
BremB\_\{\\mathrm\{rem\}\}shots to pair
\(i,j\)\(i,j\)in proportion to
wij/∑k≤ℓwkℓw\_\{ij\}/\\sum\_\{k\\leq\\ell\}w\_\{k\\ell\}\(integer rounding with leftover redistributed by descending
ww\)
10:Update
𝐊^\\widehat\{\\mathbf\{K\}\}from accumulated counts
11:\{*Stage 3: posterior inference jitter\.*\}
12:Compute
j←min\(n⋅K^\(1−K^\)/s¯,0\.5\)j\\leftarrow\\min\\\!\\bigl\(\\sqrt\{n\\cdot\\overline\{\\widehat\{K\}\(1\-\\widehat\{K\}\)/s\}\},\\;0\.5\\bigr\)\(Eq\.[19](https://arxiv.org/html/2606.28833#Sx3.E19); matrix\-concentration jitter\)
12:noisy kernel
𝐊^\\widehat\{\\mathbf\{K\}\}and inference\-time noise
σn2\+j\\sigma\_\{n\}^\{2\}\+jfor use in
\(𝐊^\+\(σn2\+j\)I\)−1𝐲\(\\widehat\{\\mathbf\{K\}\}\+\(\\sigma\_\{n\}^\{2\}\+j\)I\)^\{\-1\}\\mathbf\{y\}
Algorithm 2Sparse\-VFE AQKA\-GP shot allocation \(inducing\-point, fixed\-budget\)0:shot budget
BB, training set
\(X,𝐲\)\(X,\\mathbf\{y\}\),
mminducing points
ZZ, observation noise
σn2\\sigma\_\{n\}^\{2\}, floor
ρ\\rho, warm
ρw\\rho\_\{w\}
1:Define entry list
ℰ=upper\(Kuu\)∪Kfu\\mathcal\{E\}=\\mathrm\{upper\}\(K\_\{uu\}\)\\cup K\_\{fu\}//m\(m\+1\)2\+nm\\frac\{m\(m\+1\)\}\{2\}\+nmentries
2:Distribute
⌊\(ρ\+ρw\)B⌋\\lfloor\(\\rho\+\\rho\_\{w\}\)B\\rfloorshots uniformly over
ℰ\\mathcal\{E\}, sample
K^uu,K^fu\\widehat\{K\}\_\{uu\},\\widehat\{K\}\_\{fu\}
3:Compute
A←K^uu\+σn−2K^fu⊤K^fuA\\leftarrow\\widehat\{K\}\_\{uu\}\+\\sigma\_\{n\}^\{\-2\}\\widehat\{K\}\_\{fu\}^\{\\top\}\\widehat\{K\}\_\{fu\},
β←σn−2A−1K^fu⊤𝐲\\beta\\leftarrow\\sigma\_\{n\}^\{\-2\}A^\{\-1\}\\widehat\{K\}\_\{fu\}^\{\\top\}\\mathbf\{y\},
αq←σn−2\(𝐲−K^fuβ\)\\alpha\_\{q\}\\leftarrow\\sigma\_\{n\}^\{\-2\}\(\\mathbf\{y\}\-\\widehat\{K\}\_\{fu\}\\beta\)
4:Form sensitivities \(Eq\.[23](https://arxiv.org/html/2606.28833#Sx3.E23)\):
Sj,kuu=\|βjβk\|,Si,jfu=\|αq,iβj\|S^\{uu\}\_\{j,k\}=\|\\beta\_\{j\}\\beta\_\{k\}\|,\\quad S^\{fu\}\_\{i,j\}=\|\\alpha\_\{q,i\}\\beta\_\{j\}\|
5:Form Neyman weights
we=SeK^e\(1−K^e\)w\_\{e\}=S\_\{e\}\\sqrt\{\\widehat\{K\}\_\{e\}\(1\-\\widehat\{K\}\_\{e\}\)\}for
e∈ℰe\\in\\mathcal\{E\}
6:Distribute remaining
BremB\_\{\\mathrm\{rem\}\}shots over
ℰ\\mathcal\{E\}in proportion to
we/∑e′we′w\_\{e\}/\\sum\_\{e^\{\\prime\}\}w\_\{e^\{\\prime\}\}
7:Update
K^uu,K^fu\\widehat\{K\}\_\{uu\},\\widehat\{K\}\_\{fu\}from accumulated counts
7:noisy blocks
K^uu,K^fu\\widehat\{K\}\_\{uu\},\\widehat\{K\}\_\{fu\}for use in SGPR/SVGP/FITC posterior
### Posterior Inference Under Noisy𝐊^\\widehat\{\\mathbf\{K\}\}
A second, separate failure mode is that the GP posterior itself is ill\-conditioned when𝐊^\\widehat\{\\mathbf\{K\}\}has shot noise —\(𝐊^\+σn2I\)−1\(\\widehat\{\\mathbf\{K\}\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}blows up the noise\. We handle this by inflating the inference\-time noise toσn2\+j\\sigma\_\{n\}^\{2\}\+j, where the*jitter*jjis set adaptively from the empirical shot\-noise variance:
j=n⋅K^ij\(1−K^ij\)sij¯,j=\\sqrt\{n\}\\cdot\\overline\{\\tfrac\{\\widehat\{K\}\_\{ij\}\(1\-\\widehat\{K\}\_\{ij\}\)\}\{s\_\{ij\}\}\},\(19\)i\.e\., the per\-entry Bernoulli variance averaged over pairs and scaled byn\\sqrt\{n\}to account for cumulative spectral perturbation of the Gram matrix\. This calibrated jitter is used in both predictive mean and predictive variance\.
### Extension to Sparse Inducing\-Point GPs
The full\-GP allocator above scales withn\(n\+1\)/2n\(n\+1\)/2kernel entries\. For larger datasets, the standard remedy is to introducem≪nm\\ll n*inducing points*Z∈𝒳mZ\\in\\mathcal\{X\}^\{m\}and approximate the posterior using only the smaller blocksKuu∈ℝm×mK\_\{uu\}\\in\\mathbb\{R\}^\{m\\times m\}andKfu∈ℝn×mK\_\{fu\}\\in\\mathbb\{R\}^\{n\\times m\}\. The kernel budget then drops fromΘ\(n2\)\\Theta\(n^\{2\}\)toΘ\(mn\+m2\)\\Theta\(mn\+m^\{2\}\)entries — a3030–100×100\\timesreduction at our scales\. We show that the Neyman rule extends naturally to this family, covering the three main sparse GP methods used in practice:
#### \(a\) Sparse VFE / SGPR\(Titsias[2009](https://arxiv.org/html/2606.28833#bib.bib17)\)\.
The collapsed variational lower bound depends on the inducing blocks via
A\\displaystyle A=Kuu\+σn−2Kfu⊤Kfu,β=σn−2A−1Kfu⊤𝐲,\\displaystyle=K\_\{uu\}\+\\sigma\_\{n\}^\{\-2\}K\_\{fu\}^\{\\top\}K\_\{fu\},\\quad\\beta=\\sigma\_\{n\}^\{\-2\}A^\{\-1\}K\_\{fu\}^\{\\top\}\\mathbf\{y\},\(20\)μ∗\(x∗\)\\displaystyle\\mu\_\{\*\}\(x\_\{\*\}\)=𝐤∗u⊤β\.\\displaystyle=\\mathbf\{k\}\_\{\*u\}^\{\\top\}\\beta\.\(21\)Differentiatingμ∗\\mu\_\{\*\}via the resolvent identity yields rank\-1 sensitivities
sensvfe\(Kfu\[i,j\]\)\\displaystyle\\mathrm\{sens\}^\{\\mathrm\{vfe\}\}\(K\_\{fu\}\[i,j\]\)=\|αq,iβj\|,\\displaystyle=\|\\alpha\_\{q,i\}\\,\\beta\_\{j\}\|,\(22\)sensvfe\(Kuu\[j,k\]\)\\displaystyle\\mathrm\{sens\}^\{\\mathrm\{vfe\}\}\(K\_\{uu\}\[j,k\]\)=\|βjβk\|,\\displaystyle=\|\\beta\_\{j\}\\beta\_\{k\}\|,\(23\)whereαq=σn−2\(𝐲−Kfuβ\)\\alpha\_\{q\}=\\sigma\_\{n\}^\{\-2\}\(\\mathbf\{y\}\-K\_\{fu\}\\beta\)is the residual vector\. These are the inducing\-point analogues of the full\-GP predictive coupling\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|from Eq\. \([13](https://arxiv.org/html/2606.28833#Sx3.E13)\)\.
#### \(b\) Variational SVGP\(Hensmanet al\.[2013](https://arxiv.org/html/2606.28833#bib.bib5)\)\.
For minibatch\-friendly stochastic VI,q\(𝐮\)=𝒩\(𝐦,S\)q\(\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{m\},S\)is kept as variational parameters\. At the optimum𝐦∗=β\\mathbf\{m\}^\{\*\}=\\betaandS∗S^\{\*\}recovers the VFE posterior, so the AQKA\-VFE sensitivities of \([23](https://arxiv.org/html/2606.28833#Sx3.E23)\) apply throughout SVGP training without modification\.
#### \(c\) FITC\(Snelson and Ghahramani[2005](https://arxiv.org/html/2606.28833#bib.bib20)\)\.
FITC adds a diagonal correctiondiag\(Kff−KfuKuu−1Kuf\)\\mathrm\{diag\}\(K\_\{ff\}\-K\_\{fu\}K\_\{uu\}^\{\-1\}K\_\{uf\}\)to the noise covariance\. The resulting sensitivity has the same rank\-1 structure plus a heteroscedastic correction proportional to the diagonal residual; AQKA\-FITC requires only minor modification toαq\\alpha\_\{q\}\.
#### General principle\.
For any sparse posterior with structureμ∗\(x∗\)=𝐤∗u⊤β\+g\(𝐊^\)\\mu\_\{\*\}\(x\_\{\*\}\)=\\mathbf\{k\}\_\{\*u\}^\{\\top\}\\beta\+g\(\\widehat\{\\mathbf\{K\}\}\)whereggis differentiable in𝐊\\mathbf\{K\}, the Neyman minimum\-variance allocation reduces to Neyman\-style weighting of the rank\-1 outer productαβ⊤\\alpha\\beta^\{\\top\}\. The unifying observation is that all inducing\-point methods couple labels to predictions through this single low\-rank coupling, and AQKA\-GP can be ported by reading off\(α,β\)\(\\alpha,\\beta\)from the relevant posterior\.
#### Extension to Deep GPs \(DSVI\)\.
The same construction extends to a Deep GP\(Damianou and Lawrence[2013](https://arxiv.org/html/2606.28833#bib.bib3); Salimbeni and Deisenroth[2017](https://arxiv.org/html/2606.28833#bib.bib12)\)withLLlayersy=fL∘⋯∘f1\(x\)\+εy=f\_\{L\}\\circ\\cdots\\circ f\_\{1\}\(x\)\+\\varepsilon, eachflf\_\{l\}a sparse VFE GP with its own inducing block\(Kuu,l,Kfu,l\)\(K\_\{uu,l\},K\_\{fu,l\}\)\. The chain rule gives a layer\-wise sensitivity
sensdgp\(Kfu,l\[i,j\]\)=\|αq,l,iβl,j\|⋅∏l′\>l‖βl′‖,\\mathrm\{sens\}^\{\\mathrm\{dgp\}\}\(K\_\{fu,l\}\[i,j\]\)=\|\\alpha\_\{q,l,i\}\\beta\_\{l,j\}\|\\cdot\\prod\_\{l^\{\\prime\}\>l\}\\\|\\beta\_\{l^\{\\prime\}\}\\\|,\(24\)where the product term encodes downstream gradient amplification\. AQKA\-DGP allocates a global budgetBBacross all layer blocks proportional to these \(variance\-weighted\) sensitivities\. We give a small empirical test of this construction in our extended discussion\.
## Theory
We collect four short results that justify the GP\-specific design choices in AQKA\-GP\. All proofs are elementary and given in the Appendix; they pin down \(i\) the Neyman minimum\-variance allocation, \(ii\) the rate at which kernel error propagates to GP predictions, \(iii\) why dense data has predictive heterogeneity, and \(iv\) when the high uniform floor is necessary\.
### Neyman Allocation Is the Unique Minimizer
###### Proposition 1\(Neyman minimum\-variance allocation; restated\)\.
Letℒ\(𝐊\)\\mathcal\{L\}\(\\mathbf\{K\}\)be twice differentiable in𝐊\\mathbf\{K\}\. Under independent Bernoulli noise with variancesVar\(K^ij\)=Kij\(1−Kij\)/sij\\mathrm\{Var\}\(\\widehat\{K\}\_\{ij\}\)=K\_\{ij\}\(1\-K\_\{ij\}\)/s\_\{ij\}, the leading\-order expected loss inflationΦ\(𝐬\):=∑i≤jgij2Kij\(1−Kij\)/sij\\Phi\(\\bm\{s\}\):=\\sum\_\{i\\leq j\}g\_\{ij\}^\{2\}\\,K\_\{ij\}\(1\-K\_\{ij\}\)/s\_\{ij\}under the budget constraint∑sij=B\\sum s\_\{ij\}=Bis uniquely minimized at
sij∗=B\|gij\|Kij\(1−Kij\)∑k≤ℓ\|gkℓ\|Kkℓ\(1−Kkℓ\)\.s^\{\*\}\_\{ij\}\\;=\\;\\frac\{B\\,\|g\_\{ij\}\|\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\}\{\\sum\_\{k\\leq\\ell\}\|g\_\{k\\ell\}\|\\sqrt\{K\_\{k\\ell\}\(1\-K\_\{k\\ell\}\)\}\}\.The minimum value isΦ\(𝐬∗\)=\(∑i≤j\|gij\|Kij\(1−Kij\)\)2/B\\Phi\(\\bm\{s\}^\{\*\}\)=\(\\sum\_\{i\\leq j\}\|g\_\{ij\}\|\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\)^\{2\}/B, achieving the Neyman lower bound\.
A short Lagrangian argument suffices\. The point is that no other allocation rule can do better at leading order, so the question reduces to choosinggijg\_\{ij\}for the GP task at hand — which Section[Three Closed\-Form Sensitivities](https://arxiv.org/html/2606.28833#Sx3.SSx1)settles\.
### Kernel Error to Posterior Error: Linear Propagation
###### Proposition 2\(Posterior error in operator norm\)\.
LetA=𝐊\+σn2IA=\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}IandA^=A\+Δ\\widehat\{A\}=A\+\\Delta, with‖Δ‖op<σn2\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}<\\sigma\_\{n\}^\{2\}\. Then for any𝐲∈ℝn\\mathbf\{y\}\\in\\mathbb\{R\}^\{n\}with‖𝐲‖≤M\\\|\\mathbf\{y\}\\\|\\leq M,
‖𝜶^−𝜶‖≤‖Δ‖op⋅Mσn2\(σn2−‖Δ‖op\),\\\|\\widehat\{\\bm\{\\alpha\}\}\-\\bm\{\\alpha\}\\\|\\;\\leq\\;\\frac\{\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\cdot M\}\{\\sigma\_\{n\}^\{2\}\(\\sigma\_\{n\}^\{2\}\-\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\)\},and for any test point with‖𝐤∗‖≤M\\\|\\mathbf\{k\}\_\{\*\}\\\|\\leq M,
\|μ^∗\(x∗\)−μ∗\(x∗\)\|≤M2‖Δ‖opσn2\(σn2−‖Δ‖op\)\.\|\\widehat\{\\mu\}\_\{\*\}\(x\_\{\*\}\)\-\\mu\_\{\*\}\(x\_\{\*\}\)\|\\;\\leq\\;\\frac\{M^\{2\}\\,\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\}\{\\sigma\_\{n\}^\{2\}\(\\sigma\_\{n\}^\{2\}\-\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\)\}\.
The bound is tight up to constants and shows the1/σn41/\\sigma\_\{n\}^\{4\}amplification we exploit: a kernel perturbation of operator\-norm sizeε\\varepsiloncauses a predictive error of orderε/σn4\\varepsilon/\\sigma\_\{n\}^\{4\}\. The KRR predictor inherits the same amplification with the ridgeλ\\lambdain place ofσn2\\sigma\_\{n\}^\{2\}\(since both play the same role in\(𝐊\+λI\)−1\(\\mathbf\{K\}\+\\lambda I\)^\{\-1\}\); the structural reason GP regression is shot\-hungrier than KRR classification is therefore*not*the inverse itself but the sign\-robustness of0/10/1accuracy — classification only needs the sign of the predictor to survive shot noise, while GP NLL and predictive variance care about its magnitude and full spectrum\. With Bernoulli shot noise,𝔼\[‖Δ‖op\]≲n/smin\\mathbb\{E\}\[\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\]\\lesssim n/\\sqrt\{s\_\{\\min\}\}wheresmins\_\{\\min\}is the minimum per\-pair shot count, recovering the necessity of a uniform floor: without the floor,smin=0s\_\{\\min\}=0on some pairs and the bound is vacuous\.
### Spectral Decay Amplifies Predictive Heterogeneity
###### Lemma 1\(Covariance of𝜶\\bm\{\\alpha\}under random labels\)\.
Let𝐊=∑k=1nλk𝐮k𝐮k⊤\\mathbf\{K\}=\\sum\_\{k=1\}^\{n\}\\lambda\_\{k\}\\,\\mathbf\{u\}\_\{k\}\\mathbf\{u\}\_\{k\}^\{\\top\}be the eigendecomposition withλ1≥⋯≥λn≥0\\lambda\_\{1\}\\geq\\dots\\geq\\lambda\_\{n\}\\geq 0, and let𝐲=∑kβk𝐮k\\mathbf\{y\}=\\sum\_\{k\}\\beta\_\{k\}\\mathbf\{u\}\_\{k\}withβk∼iid𝒩\(0,τ2\)\\beta\_\{k\}\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}\\mathcal\{N\}\(0,\\tau^\{2\}\)\. Then𝛂=A−1𝐲\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}is mean\-zero Gaussian with covariance
𝔼\[𝜶𝜶⊤\]=τ2∑k=1n𝐮k𝐮k⊤\(λk\+σn2\)2,\\mathbb\{E\}\[\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}\]\\;=\\;\\tau^\{2\}\\sum\_\{k=1\}^\{n\}\\frac\{\\mathbf\{u\}\_\{k\}\\mathbf\{u\}\_\{k\}^\{\\top\}\}\{\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{2\}\},soVar\(αi\)=τ2∑kuki2/\(λk\+σn2\)2\\mathrm\{Var\}\(\\alpha\_\{i\}\)=\\tau^\{2\}\\sum\_\{k\}u\_\{ki\}^\{2\}/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{2\}and the squared coupling𝔼\[αi2αj2\]\\mathbb\{E\}\[\\alpha\_\{i\}^\{2\}\\alpha\_\{j\}^\{2\}\]inherits a quartic spectral weight1/\(λk\+σn2\)41/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{4\}on each eigen\-mode\.
The lemma is exact\. We append it to an*observation*\(not a proven theorem\) about when the pair\-coupling𝔼\[\|αiαj\|\]\\mathbb\{E\}\[\|\\alpha\_\{i\}\\alpha\_\{j\}\|\]is heterogeneous — the experiments below test the observation directly via realized AQKA\-GP gain\.
#### Observation 1 \(Spectral amplification creates pair heterogeneity\)\.
When the kernel spectrum is heavy\-tailed —λk→0\\lambda\_\{k\}\\to 0forkklarge, as in any kernel with rapid Mercer decay — the small\-eigenvalue contribution dominates the𝜶\\bm\{\\alpha\}\-covariance through the quartic factor1/\(λk\+σn2\)41/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{4\}\. Combined with finite\-sample fluctuations ofuki2u\_\{ki\}^\{2\}around1/n1/n\(which are amplified by the same quartic factor\), this typically produces a pair\-coupling matrix𝔼\[\|αiαj\|\]\\mathbb\{E\}\[\|\\alpha\_\{i\}\\alpha\_\{j\}\|\]that concentrates on a few\(i,j\)\(i,j\)pairs\. A useful descriptive proxy for the effect size is the spectral effective rankr¯:=∑k1/\(1\+σn2/λk\)2\\bar\{r\}:=\\sum\_\{k\}1/\(1\+\\sigma\_\{n\}^\{2\}/\\lambda\_\{k\}\)^\{2\}; smallerr¯\\bar\{r\}\(relative tonn\) is empirically associated with larger heterogeneity\. We do not claim a quantitative bound on theℓ2/ℓ1\\ell\_\{2\}/\\ell\_\{1\}ratio of𝔼\[\|αiαj\|\]\\mathbb\{E\}\[\|\\alpha\_\{i\}\\alpha\_\{j\}\|\]in terms ofr¯\\bar\{r\}alone — such a bound would require additional control of the eigenvector structureukiu\_\{ki\}that we leave to future work\.
The mechanism is spectral amplification, not eigenvector localization\. For translation\-invariant kernels such as RBF on\[0,1\]d\[0,1\]^\{d\}, low\-λk\\lambda\_\{k\}eigenvectors are in fact delocalized \(Fourier\-like\), so a story based on eigenvector localization would predict the opposite\. Kernels with flat spectrum \(e\.g\., the identity kernel\) would predict no gain\. The dependence on kernel choice is real, which is why we test it empirically with the ZZ feature map \(Result 10\): the gain transfers from RBF to ZZ kernels \(mean−10%\-10\\%\) even though the kernels are different feature maps\.
In our experiments this is consistent with what we measure: dense GP\-prior data shows−21%\-21\\%RMSE gain \(Figure[1](https://arxiv.org/html/2606.28833#Sx5.F1)right\), planted\-sparse−11%\-11\\%\(left\), and the four UCI datasets \(all dense in the kernel sense\) fall in the same band\. Genuine ZZ and Pauli\-Z quantum kernels on quantum\-natural data \(Result 10, Study 2\) show−13\-13to−15%\-15\\%at the low\-budget regime with pairedp<0\.05p<0\.05; UCI data embedded into a ZZ kernel \(Result 10, Study 1\) shows null gain, consistent with the kernel concentrating away the heterogeneity\.
### Why the Uniform Floor Helps: Coverage Argument and Empirical Calibration
The uniform floorρ\\rhoplays two distinct roles: a*coverage*role that guarantees no pair is missed, and a*shot\-quality*role that guarantees every covered pair has enough shots forK^ij\\widehat\{K\}\_\{ij\}to be reasonably calibrated\. The first role is formal; the second is empirical, calibrated against Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)\.
###### Proposition 3\(Coverage\-floor sufficient condition\)\.
LetSSbe the set of pairs that receive zero shots under an allocation, with the kernel estimator defaulting toK^ij=c\\widehat\{K\}\_\{ij\}=conSS\. Suppose at least a fractionη∈\(0,1\]\\eta\\in\(0,1\]of pairs inSShave\|c−Kij\|≥d\|c\-K\_\{ij\}\|\\geq d\. LetΔ:=𝐊^−𝐊\\Delta:=\\widehat\{\\mathbf\{K\}\}\-\\mathbf\{K\}be supported onS∪S⊤S\\cup S^\{\\top\}\. Then
‖Δ‖F≥d2η\|S\|,‖Δ‖op≥‖Δ‖Fn≥d2η\|S\|n\.\\\|\\Delta\\\|\_\{F\}\\;\\geq\\;d\\sqrt\{2\\eta\|S\|\},\\qquad\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\;\\geq\\;\\frac\{\\\|\\Delta\\\|\_\{F\}\}\{\\sqrt\{n\}\}\\;\\geq\\;d\\sqrt\{\\frac\{2\\eta\|S\|\}\{n\}\}\.Hence the only way to keep‖Δ‖op→0\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\to 0\(which Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)requires for predictive\-error consistency\) under fixedd,ηd,\\etais to drive\|S\|→0\|S\|\\to 0, i\.e\., to cover every pair\. The smallest uniform\-floor fractionρ\\rhothat guarantees\|S\|=0\|S\|=0is
ρ≥n\(n\+1\)2B\.\\rho\\;\\geq\\;\\frac\{n\(n\+1\)\}\{2B\}\.
The bound shows that*some*uniform floor is necessary \(otherwise\|S\|\>0\|S\|\>0and the kernel\-perturbation operator\-norm has a non\-zero lower bound that prevents the predictive\-error upper bound of Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)from vanishing\); it does*not*prove that the predictive error itself has a matching lower bound\. The bound is a necessary condition, not a sufficient one\. The catastrophic empirical regime of Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)reflects an additional shot\-quality requirement that the Frobenius argument does not touch: each covered pair must receive enough shots that itsK^ij\\widehat\{K\}\_\{ij\}is calibrated, not just non\-default\. The transitionρ≤0\.1→ρ≥0\.2\\rho\\leq 0\.1\\to\\rho\\geq 0\.2in Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)aligns withρ=n2/\(2B\)≈0\.1\\rho=n^\{2\}/\(2B\)\\approx 0\.1forn=200,B=2×105n=200,B=2\\\!\\times\\\!10^\{5\}, but the additional jump toρ=0\.5\\rho=0\.5for full robustness is calibrated empirically, not theoretically derived\. We document this honestly:ρ=0\.5\\rho=0\.5is a robust default that adds shot\-quality margin on top of the coverage requirement\.
#### Whyn⋅K\(1−K\)/s¯\\sqrt\{n\\cdot\\overline\{K\(1\-K\)/s\}\}jitter\.
The per\-entry shot noise of𝐊^\\widehat\{\\mathbf\{K\}\}has varianceσent2\(i,j\)=Kij\(1−Kij\)/sij\\sigma\_\{\\mathrm\{ent\}\}^\{2\}\(i,j\)=K\_\{ij\}\(1\-K\_\{ij\}\)/s\_\{ij\}\. Treating the upper\-triangular shot\-noise entries as a symmetric matrix with sub\-Gaussian off\-diagonal entries of common varianceσ¯ent2:=K\(1−K\)/s¯\\bar\{\\sigma\}\_\{\\mathrm\{ent\}\}^\{2\}:=\\overline\{K\(1\-K\)/s\}, standard Wigner\-type concentration\(Tropp[2015](https://arxiv.org/html/2606.28833#bib.bib18)\)gives
𝔼\[‖Δ‖op\]≤Cn⋅σ¯ent2=Cn⋅K\(1−K\)/s¯\\mathbb\{E\}\[\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\]\\;\\leq\\;C\\sqrt\{n\\cdot\\bar\{\\sigma\}\_\{\\mathrm\{ent\}\}^\{2\}\}\\;=\\;C\\sqrt\{n\\cdot\\overline\{K\(1\-K\)/s\}\}for an absolute constantCC\. This Wigner\-type bound*motivates*a scaling of the formj∼n⋅\(something average\-per\-entry\)j\\sim\\sqrt\{n\}\\cdot\(\\text\{something average\-per\-entry\}\), but does not by itself determine the constant or the precise form\. The theory\-motivated choice would bejthy=n⋅K\(1−K\)/s¯j\_\{\\mathrm\{thy\}\}=\\sqrt\{n\\cdot\\overline\{K\(1\-K\)/s\}\}, which would absorb the expected operator\-norm perturbation in expectation\. The code variant of Eq\.[19](https://arxiv.org/html/2606.28833#Sx3.E19)usesjcode=n⋅K\(1−K\)/s¯j\_\{\\mathrm\{code\}\}=\\sqrt\{n\}\\cdot\\overline\{K\(1\-K\)/s\}instead — an\\sqrt\{n\}prefactor times the per\-entry variance rather than its square root\. The ratio isjcode/jthy=K\(1−K\)/s¯j\_\{\\mathrm\{code\}\}/j\_\{\\mathrm\{thy\}\}=\\sqrt\{\\overline\{K\(1\-K\)/s\}\}, which depends on the average shots\-per\-pairssand is*not*a constant factor \(ats=50s=50,K=0\.5K=0\.5it is∼\\sim1/141/14\)\. At the budgets we use, the theory choicejthyj\_\{\\mathrm\{thy\}\}would exceed the clip\-to\-0\.50\.5ceiling and reduce to0\.50\.5; the code’sjcodej\_\{\\mathrm\{code\}\}is below this ceiling and is not clipped\. So neither the theory derivation nor the code formula determines the effective jitter in practice — the clip\-to\-0\.50\.5ceiling does\. We therefore do*not*claim the jitter is principled in a derivation\-from\-theory sense; the matrix\-concentration argument motivates then\\sqrt\{n\}scaling, but the precise constant is a calibrated default, and we report this honestly\. We do not see a difference in test RMSE betweenjthyj\_\{\\mathrm\{thy\}\}andjcodej\_\{\\mathrm\{code\}\}within11SE because both are clipped near0\.50\.5in the operating range\.
### Per\-Task Optimality Under the Estimation\-Error Objective
Eq\. \([9](https://arxiv.org/html/2606.28833#Sx2.E9)\) in Section[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3)targeted the expected loss bias𝔼\[ℒ^−ℒ\]\\mathbb\{E\}\[\\widehat\{\\mathcal\{L\}\}\-\\mathcal\{L\}\]; the Fisher approximationHij,ij≈gij2H\_\{ij,ij\}\\approx g\_\{ij\}^\{2\}entered as a surrogate for the unknown Hessian\. We now*change the optimization target*to the expected squared estimation error𝔼\[\(L^−L\)2\]\\mathbb\{E\}\[\(\\widehat\{L\}\-L\)^\{2\}\]for a downstream estimandL∈\{μ∗\(x∗\),ℒ\(𝐊\),ei\}L\\in\\\{\\mu\_\{\*\}\(x\_\{\*\}\),\\mathcal\{L\}\(\\mathbf\{K\}\),e\_\{i\}\\\}\(test predictive mean, marginal likelihood, LOO residual\)\. Under this new objective the leading\-order coefficient is exactlygij2=\(∂L/∂Kij\)2g\_\{ij\}^\{2\}=\(\\partial L/\\partial K\_\{ij\}\)^\{2\}with no Fisher approximation, so the three sensitivities become rigorous Neyman weights for the estimation\-error objective\. We are honest that this is a different optimization target from Eq\. \([9](https://arxiv.org/html/2606.28833#Sx2.E9)\), not a tighter analysis of the same target; the experiments below directly measure squared estimation errors \(RMSE, NLL error, LOO MSE\), which is what the new propositions analyze\.
###### Proposition 4\(Exact predictive\-MSE sensitivity\)\.
Letπ\\pibe the test distribution andM:=𝔼x∗∼π\[𝐤∗𝐤∗⊤\]∈ℝn×nM:=\\mathbb\{E\}\_\{x\_\{\*\}\\sim\\pi\}\[\\mathbf\{k\}\_\{\*\}\\mathbf\{k\}\_\{\*\}^\{\\top\}\]\\in\\mathbb\{R\}^\{n\\times n\}\. The expected predictive\-mean MSE under Bernoulli shot noise satisfies, to leading order in𝚫\\bm\{\\Delta\},
𝔼\[𝔼x∗\|μ^∗\(x∗\)−μ∗\(x∗\)\|2\]=∑i≤jKij\(1−Kij\)sij⋅cij2\\mathbb\{E\}\\bigl\[\\,\\mathbb\{E\}\_\{x\_\{\*\}\}\|\\widehat\{\\mu\}\_\{\*\}\(x\_\{\*\}\)\-\\mu\_\{\*\}\(x\_\{\*\}\)\|^\{2\}\\,\\bigr\]\\;=\\;\\sum\_\{i\\leq j\}\\frac\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\{s\_\{ij\}\}\\cdot c\_\{ij\}^\{2\}withcij2=αi2Bjj\+αj2Bii\+2αiαjBijc\_\{ij\}^\{2\}=\\alpha\_\{i\}^\{2\}B\_\{jj\}\+\\alpha\_\{j\}^\{2\}B\_\{ii\}\+2\\alpha\_\{i\}\\alpha\_\{j\}B\_\{ij\}andB:=A−1MA−1B:=A^\{\-1\}MA^\{\-1\},A=𝐊\+σn2IA=\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I,𝛂=A−1𝐲\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}\. The Neyman\-optimal allocation is thereforesij∗∝Kij\(1−Kij\)⋅cijs^\{\*\}\_\{ij\}\\propto\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\\cdot c\_\{ij\}\.
The exactcijc\_\{ij\}couples each pair through three quantities:αi,αj\\alpha\_\{i\},\\alpha\_\{j\}, and the test\-direction matrixBB\. Our simplified sensitivitysensijpred=\|αiαj\|\\mathrm\{sens\}^\{\\mathrm\{pred\}\}\_\{ij\}=\|\\alpha\_\{i\}\\alpha\_\{j\}\|is the rank\-1 specialization ofcijc\_\{ij\}whenM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}: substituting givesB=A−1𝐲𝐲⊤A−1=𝜶𝜶⊤B=A^\{\-1\}\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}A^\{\-1\}=\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}, hencecij2=αi2αj2\+αj2αi2\+2\(αiαj\)2=4\(αiαj\)2c\_\{ij\}^\{2\}=\\alpha\_\{i\}^\{2\}\\alpha\_\{j\}^\{2\}\+\\alpha\_\{j\}^\{2\}\\alpha\_\{i\}^\{2\}\+2\(\\alpha\_\{i\}\\alpha\_\{j\}\)^\{2\}=4\(\\alpha\_\{i\}\\alpha\_\{j\}\)^\{2\}andcij=2\|αiαj\|c\_\{ij\}=2\|\\alpha\_\{i\}\\alpha\_\{j\}\|\. TheM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}regime captures test distributions whose second moment lies in the label direction — a natural assumption when train and test inputs share a distribution and the labels carry the task signal\. For other test distributionscijc\_\{ij\}is not proportional to\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|\(e\.g\.,M=IM=Iproduces a term‖Aj,:−1‖2\\\|A^\{\-1\}\_\{j,:\}\\\|^\{2\}rather thanαj2\\alpha\_\{j\}^\{2\}\)\.
###### Corollary 1\(Regret of the rank\-1 sensitivity\)\.
Lets~ij∝Kij\(1−Kij\)⋅\|αiαj\|\\widetilde\{s\}\_\{ij\}\\propto\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}\\cdot\|\\alpha\_\{i\}\\alpha\_\{j\}\|\. Its Neyman regret relative tosij∗s^\{\*\}\_\{ij\}is
Regret\(s~∥s∗\)≤\(κ\(B\)−1\)⋅Φ\(s∗\),\\mathrm\{Regret\}\(\\widetilde\{s\}\\,\\\|\\,s^\{\*\}\)\\;\\leq\\;\\bigl\(\\kappa\(B\)\-1\\bigr\)\\cdot\\Phi\(s^\{\*\}\),whereκ\(B\)\\kappa\(B\)is the spectral condition number ofBBandΦ\(⋅\)\\Phi\(\\cdot\)is the Neyman objective\. In particularRegret=0\\mathrm\{Regret\}=0whenM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}\(soB=𝛂𝛂⊤B=\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}is rank\-1 in the𝛂\\bm\{\\alpha\}direction\); it is finite wheneverπ\\pihas finite second moment\.
###### Proposition 5\(NLL sensitivity is the exact gradient\)\.
*Convention*: throughout this paper we treatKijK\_\{ij\}fori<ji<jas a single Gram coordinate that simultaneously controlsAijA\_\{ij\}andAjiA\_\{ji\}via symmetry\. Under this convention, the marginal\-likelihood gradient is
gijmarg=\[A−1\]ij−αiαj,g^\{\\mathrm\{marg\}\}\_\{ij\}\\;=\\;\[A^\{\-1\}\]\_\{ij\}\-\\alpha\_\{i\}\\alpha\_\{j\},which we write as12\(\[A−1\]ij−αiαj\)\+12\(\[A−1\]ji−αjαi\)\\tfrac\{1\}\{2\}\(\[A^\{\-1\}\]\_\{ij\}\-\\alpha\_\{i\}\\alpha\_\{j\}\)\+\\tfrac\{1\}\{2\}\(\[A^\{\-1\}\]\_\{ji\}\-\\alpha\_\{j\}\\alpha\_\{i\}\)in symmetric form\. \(The factor\-of\-22paper\-vs\-code discrepancy that earlier reviewers flagged is a coordinate convention: ifAijA\_\{ij\}andAjiA\_\{ji\}are treated as independent coordinates, each derivative is half of the value above\. We adopt the single\-coordinate convention; the factor is absorbed into the Neyman normalization either way\.\) Under the expected squared estimation\-error objective𝔼\[\(ℒ^−ℒ\)2\]\\mathbb\{E\}\[\(\\widehat\{\\mathcal\{L\}\}\-\\mathcal\{L\}\)^\{2\}\], the leading\-order coefficient is\(gijmarg\)2\(g^\{\\mathrm\{marg\}\}\_\{ij\}\)^\{2\}and the Neyman allocation with\|gijmarg\|\|g^\{\\mathrm\{marg\}\}\_\{ij\}\|minimizes this objective globally to leading order\.
###### Proposition 6\(LOO sensitivity, exact form and simplification\)\.
The closed\-form GP LOO residual isei=αi/\[A−1\]iie\_\{i\}=\\alpha\_\{i\}/\[A^\{\-1\}\]\_\{ii\}, and the gradient of the LOO sum of squared residualsℒloo=12∑iei2\\mathcal\{L\}^\{\\mathrm\{loo\}\}=\\tfrac\{1\}\{2\}\\sum\_\{i\}e\_\{i\}^\{2\}with respect to the symmetric Gram coordinateKijK\_\{ij\}\(i<ji<j\) takes the form \(full derivation in Appendix A\)
∂ℒloo/∂Kij=\\displaystyle\\partial\\mathcal\{L\}^\{\\mathrm\{loo\}\}/\\partial K\_\{ij\}\\;=\\;−∑kek\[A−1\]kiαj\+\[A−1\]kjαi\[A−1\]kk⏟α\-coupling\\displaystyle\\underbrace\{\-\\sum\_\{k\}e\_\{k\}\\,\\tfrac\{\[A^\{\-1\}\]\_\{ki\}\\alpha\_\{j\}\+\[A^\{\-1\}\]\_\{kj\}\\alpha\_\{i\}\}\{\[A^\{\-1\}\]\_\{kk\}\}\}\_\{\\alpha\\text\{\-coupling\}\}\+∑kek⋅αk\[A−1\]ki\[A−1\]kj\[A−1\]kk2⏟denominator correction\.\\displaystyle\\quad\+\\;\\underbrace\{\\sum\_\{k\}e\_\{k\}\\cdot\\tfrac\{\\alpha\_\{k\}\[A^\{\-1\}\]\_\{ki\}\[A^\{\-1\}\]\_\{kj\}\}\{\[A^\{\-1\}\]\_\{kk\}^\{2\}\}\}\_\{\\text\{denominator correction\}\}\.Theα\\alpha\-coupling term contains\|ei\[A−1\]ij\|\+\|ej\[A−1\]ji\|\|e\_\{i\}\[A^\{\-1\}\]\_\{ij\}\|\+\|e\_\{j\}\[A^\{\-1\}\]\_\{ji\}\|as its\(i,j\)\(i,j\)\-localized contribution after symmetrizing; the denominator\-correction term scales likemaxk\(ek2⋅\[A−1\]ki\[A−1\]kj/\[A−1\]kk\)\\max\_\{k\}\(e\_\{k\}^\{2\}\\cdot\[A^\{\-1\}\]\_\{ki\}\[A^\{\-1\}\]\_\{kj\}/\[A^\{\-1\}\]\_\{kk\}\), bounded above by the maximum LOO predictive variance\. Our simplified sensitivitysensijloo=\|ei\[A−1\]ij\|\+\|ej\[A−1\]ji\|\\mathrm\{sens\}^\{\\mathrm\{loo\}\}\_\{ij\}=\|e\_\{i\}\\,\[A^\{\-1\}\]\_\{ij\}\|\+\|e\_\{j\}\\,\[A^\{\-1\}\]\_\{ji\}\|keeps only theα\\alpha\-coupling term and drops the denominator correction\. This is a controlled approximation, not the exact gradient\.
The three sensitivities therefore stand in different relations to their respective optima:\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|is a rank\-1 specialization of the exact predictive\-MSE Neyman weight atM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}\(Proposition[4](https://arxiv.org/html/2606.28833#Thmproposition4), Corollary[1](https://arxiv.org/html/2606.28833#Thmcorollary1)\);\|gijmarg\|\|g^\{\\mathrm\{marg\}\}\_\{ij\}\|is the exact partial derivative of the negative marginal log\-likelihood under the single\-coordinate convention, hence the exact\(∂ℒ/∂Kij\)2\\sqrt\{\(\\partial\\mathcal\{L\}/\\partial K\_\{ij\}\)^\{2\}\}Neyman weight under the squared estimation\-error objective \(Proposition[5](https://arxiv.org/html/2606.28833#Thmproposition5)\);sensijloo\\mathrm\{sens\}^\{\\mathrm\{loo\}\}\_\{ij\}drops a bounded correction term from the LOO gradient \(Proposition[6](https://arxiv.org/html/2606.28833#Thmproposition6)\)\. The empirical fact that all three deliver real gains \(Figures[1](https://arxiv.org/html/2606.28833#Sx5.F1)–[5](https://arxiv.org/html/2606.28833#Sx5.F5), Result 10\) is consistent with each capturing the dominant term of its respective Neyman objective\.
## Experiments
#### Setup\.
All experiments usentr=200n\_\{tr\}=200training points andnte=80n\_\{te\}=80–100100test points; we report mean±\\pmstandard error across55–1010seeds\. The quantum kernel is simulated asKij=e−γ‖xi−xj‖2K\_\{ij\}=e^\{\-\\gamma\\\|x\_\{i\}\-x\_\{j\}\\\|^\{2\}\}withγ\\gammachosen by the median heuristic on‖xi−xj‖2\\\|x\_\{i\}\-x\_\{j\}\\\|^\{2\}; the shot\-noise process samplesK^ij⋅sij∼Bin\(sij,Kij\)\\widehat\{K\}\_\{ij\}\\cdot s\_\{ij\}\\sim\\mathrm\{Bin\}\(s\_\{ij\},K\_\{ij\}\)\. We compare𝚞𝚗𝚒𝚏𝚘𝚛𝚖\\mathtt\{uniform\}and𝚛𝚊𝚗𝚍𝚘𝚖\\mathtt\{random\}allocation against three AQKA\-GP variants using the three sensitivities above\. Observation noise isσn=0\.3\\sigma\_\{n\}=0\.3\. The warm\-up fraction isρw=0\.1\\rho\_\{w\}=0\.1and the floor isρ=0\.5\\rho=0\.5unless noted\.
#### Synthetic settings\.
We use two synthetic regression problems withntr=200n\_\{tr\}=200\. The*dense*setting drawsf∼𝒢𝒫\(0,KRBF\)f\\sim\\mathcal\{GP\}\(0,K\_\{\\mathrm\{RBF\}\}\)from a200200\-dimensional Gaussian, then observesyi=f\(xi\)\+εiy\_\{i\}=f\(x\_\{i\}\)\+\\varepsilon\_\{i\}\. The*planted\-sparse*setting plantsm=15m=15random anchor points withc∼𝒩\(0,I\)c\\sim\\mathcal\{N\}\(0,I\)and setsf\(x\)=∑i∈anchorsciK\(x,xi\)f\(x\)=\\sum\_\{i\\in\\mathrm\{anchors\}\}c\_\{i\}K\(x,x\_\{i\}\), giving a function with highly heterogeneous predictive coupling — the GP regression analogue of the planted\-sparse setting inXuet al\.\([2026](https://arxiv.org/html/2606.28833#bib.bib1)\)\.
#### UCI benchmarks\.
We use four standard UCI regression datasets:energy\(heating load,768×8768\\times 8\),concrete\(compressive strength,1030×81030\\times 8\),kin8nm\(robot arm forward kinematics,8192×88192\\times 8\), andcalifornia\(housing,20640×820640\\times 8\)\. Features and targets are standardized using training\-set statistics only; we subsample200200training and100100test points per seed and report1010seeds\. \(An earlier version of this paper usedyachtin place ofkin8nm; we replaced it because its standardized targets have unusually low SNR which made all allocators non\-monotonic in budget, defeating the purpose of a benchmark\.\)
Figure 1:Synthetic experiments \(ntr=200n\_\{tr\}=200,1010seeds,σn=0\.3\\sigma\_\{n\}=0\.3\)\. Test RMSE vs\. shot budgetBBon \(left\) planted\-sparse data, \(right\) dense GP\-prior data\. AQKA\-GP variants beat uniform across the moderate\-to\-high budget regime; the gain is larger on dense data, where inverse propagation creates predictive heterogeneity from spectral decay alone\. At very low budget \(B≤100B\\leq 100shots/pair\), all variants underperform uniform because the warm\-up sensitivity estimate is too noisy — the catastrophic regime analyzed below\.
#### Result 1: AQKA\-GP gains on synthetic data \(Figure[1](https://arxiv.org/html/2606.28833#Sx5.F1), Figure[2](https://arxiv.org/html/2606.28833#Sx5.F2)\)\.
On the dense GP\-prior setting, the\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|sensitivity gives a−21%\-21\\%RMSE improvement atB=106B=10^\{6\}\(~5050shots/pair\)\. The leave\-one\-out and marginal\-likelihood variants also give double\-digit improvements at the same budget\. Gains shrink to−5%\-5\\%atB=2×107B=2\\times 10^\{7\}as all methods approach the oracle floor\. On the planted\-sparse setting, gains are slightly smaller but qualitatively identical:gp\_loowins by−11%\-11\\%atB=106B=10^\{6\}and−9%\-9\\%atB=5×106B=5\\times 10^\{6\}\. Figure[2](https://arxiv.org/html/2606.28833#Sx5.F2)shows the per\-seed RMSE distribution: AQKA variants are not only lower on average but have narrower spread than uniform atB=106,5×106B=10^\{6\},5\\times 10^\{6\}, and converge to the oracle floor atB=2×107B=2\\times 10^\{7\}\.
Figure 2:Per\-seed RMSE distribution on dense GP\-prior data,ntr=200n\_\{tr\}=200,1010seeds\. Violins show full distribution, dots individual seeds, horizontal mark the mean, dashed line the oracle floor\. AQKA variants concentrate below uniform across budgets and have visibly narrower spread atB=106B=10^\{6\}–5×1065\\times 10^\{6\}\.
#### Mechanism visualisation \(Figure[3](https://arxiv.org/html/2606.28833#Sx5.F3), Figure[4](https://arxiv.org/html/2606.28833#Sx5.F4)\)\.
To make the mechanism concrete, Figure[3](https://arxiv.org/html/2606.28833#Sx5.F3)shows kernel matrices, shot\-allocation maps, and the\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|sensitivity field on a single planted\-sparse seed withn=60n=60andB=50npairsB=50\\,n\_\{\\mathrm\{pairs\}\}\. Uniform spreads shots evenly \(bottom\-left\); AQKA\-GP concentrates shots along the anchor\-pair stripes that the sensitivity heatmap identifies \(bottom\-center\)\. On this seed AQKA\-GP recovers the oracle RMSE \(0\.3750\.375vs\. oracle0\.3710\.371\) while uniform reaches only0\.5360\.536\(−30%\-30\\%gain\)\. Figure[4](https://arxiv.org/html/2606.28833#Sx5.F4)quantifies the pair\-sensitivity heterogeneity as a Lorenz\-style curve: the top20%20\\%of pairs carry∼\\sim58%58\\%of the cumulative\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|mass \(Gini coefficientG≈0\.57G\\approx 0\.57\), validating Observation 1 empirically\.
Figure 3:Kernel \+ allocation heatmaps on a single planted\-sparse seed,n=60n=60,B=50npairsB=50\\,n\_\{\\mathrm\{pairs\}\}\. Top row: oracleKtrueK\_\{\\mathrm\{true\}\},K^\\widehat\{K\}under uniform and AQKA shot allocation, and the uniform\-allocation kernel error\. Bottom row: uniform shot count per pair \(flat\), AQKA shot count per pair \(concentrated on anchor\-pair stripes\), sensitivity field\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|that drives the AQKA allocation, and the AQKA kernel error\. On this seed AQKA\-GP achieves−30%\-30\\%test RMSE over uniform \(0\.3750\.375vs\.0\.5360\.536, vs\. oracle0\.3710\.371\)\.Figure 4:Pair\-sensitivity concentration \(Lorenz curve\) for\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|on the same planted\-sparse seed used for Figure[3](https://arxiv.org/html/2606.28833#Sx5.F3)\(n=60n=60\)\. The top20%20\\%of pairs carry∼\\sim58%58\\%of the cumulative sensitivity mass; Gini coefficientG≈0\.57G\\approx 0\.57\. A flat kernel would giveG=0G=0; a degenerate one\-pair kernel would giveG=1G=1\. The strong heterogeneity is what AQKA\-GP exploits\.Figure 5:UCI regression benchmarks \(ntr=200n\_\{tr\}=200,1010seeds\)\. AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|\(red, pre\-committed headline sensitivity\) delivers statistically significant gains atB=106B=10^\{6\}on3/43/4datasets \(−7%\-7\\%to−18%\-18\\%, pairedtt\-testp<0\.05p<0\.05\); see Table[1](https://arxiv.org/html/2606.28833#Sx5.T1)forpp\-values\.
#### Result 2: AQKA\-GP gains on UCI data \(Figure[5](https://arxiv.org/html/2606.28833#Sx5.F5), Table[1](https://arxiv.org/html/2606.28833#Sx5.T1)\)\.
We pre\-commit to\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|as the headline sensitivity \(the rank\-1 specialization of Proposition[4](https://arxiv.org/html/2606.28833#Thmproposition4), identified before running the UCI experiments\) and report pairedtt\-tests against uniform with1010seeds per dataset\. AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|delivers statistically significant \(p<0\.05p<0\.05\) RMSE gains atB=106B=10^\{6\}on3/43/4datasets \(−9\.8%\-9\.8\\%to−18\.4%\-18\.4\\%\) and atB=5×106B=5\\times 10^\{6\}on2/42/4datasets \(−11\.1%\-11\.1\\%on bothconcreteandkin8nm\)\. AtB=2×105B=2\\times 10^\{5\}the gain is significant only onkin8nm\(p=0\.025p=0\.025\); atB=2×107B=2\\times 10^\{7\}all methods converge toward the oracle floor and most gains are non\-significant\. We report max\-over\-variants \(LOO, marg\-lik\) separately as an oracle\-selection upper bound; the headline gp\_alpha already accounts for most of the gain\.
Table 1:Test RMSE on four UCI regression benchmarks,ntr=200n\_\{tr\}=200,1010seeds\. Headline sensitivity is gp\_alpha \(pre\-committed before running\); we report pairedtt\-testpp\-values against uniform on the headline\. Bold marks the headline gain whenp<0\.05p<0\.05\.
#### Result 3: floor ablation \(Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)\)\.
We re\-ran the dense synthetic setting withρ∈\{0,0\.1,0\.2,0\.5,0\.7\}\\rho\\in\\\{0,0\.1,0\.2,0\.5,0\.7\\\}\(55seeds,ntr=200n\_\{tr\}=200\)\. Atρ=0\\rho=0andρ=0\.1\\rho=0\.1, AQKA\-GP collapses to the catastrophic regime atB=2×105B=2\\times 10^\{5\}\(\+200%\+200\\%to\+260%\+260\\%RMSE vs uniform\)\. Atρ=0\.2\\rho=0\.2\(the classification default ofXuet al\.\([2026](https://arxiv.org/html/2606.28833#bib.bib1)\)\), the catastrophe persists atB≤2×105B\\leq 2\\times 10^\{5\}\. Atρ=0\.5\\rho=0\.5the catastrophe is eliminated at all budgets we tested, while still preserving the−15%\-15\\%to−21%\-21\\%gain at moderateBB\. Atρ=0\.7\\rho=0\.7, gains at moderateBBshrink slightly because too little budget is left for sensitivity\-driven allocation\. We useρ=0\.5\\rho=0\.5as a robust default;ρ=0\.7\\rho=0\.7is a safer choice when warm\-up shot count per pair is expected to be very low\.
Figure 6:Floor ablation on dense GP\-prior data \(ntr=200n\_\{tr\}=200,55seeds\)\. The uniform\-coverage fractionρ\\rhois essential: atρ≤0\.1\\rho\\leq 0\.1, AQKA\-GP catastrophically overconcentrates atB=2×105B=2\\times 10^\{5\}\. Atρ∈\{0\.5,0\.7\}\\rho\\in\\\{0\.5,0\.7\\\}the catastrophe is eliminated\. We useρ=0\.5\\rho=0\.5as the default; the classification setup ofXuet al\.\([2026](https://arxiv.org/html/2606.28833#bib.bib1)\)usesρ≤0\.2\\rho\\leq 0\.2\.
#### Result 4: marginal\-likelihood approximation \(Figure[7](https://arxiv.org/html/2606.28833#Sx5.F7)\)\.
The marginal\-likelihood sensitivity is motivated by hyperparameter learning, where one cares about how accuratelyℒ\(𝐊^\)\\mathcal\{L\}\(\\widehat\{\\mathbf\{K\}\}\)approximatesℒ\(𝐊\)\\mathcal\{L\}\(\\mathbf\{K\}\)\. We measure the absolute NLL error\|ℒ\(𝐊^\)−ℒ\(𝐊\)\|\|\\mathcal\{L\}\(\\widehat\{\\mathbf\{K\}\}\)\-\\mathcal\{L\}\(\\mathbf\{K\}\)\|as a function ofBBacross88seeds on both synthetic settings\. AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|and AQKA\-GP\-marg give monotonically smaller NLL error than uniform: atB=2×107B=2\\times 10^\{7\}on dense data, AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|reaches NLL error8\.58\.5vs19\.419\.4for uniform \(−56%\-56\\%\); on sparse data,8\.98\.9vs11\.011\.0\(−18%\-18\\%\)\. The kernel Frobenius error‖𝐊^−𝐊‖F/‖𝐊‖F\\\|\\widehat\{\\mathbf\{K\}\}\-\\mathbf\{K\}\\\|\_\{F\}/\\\|\\mathbf\{K\}\\\|\_\{F\}\(lower panel of Figure[7](https://arxiv.org/html/2606.28833#Sx5.F7)\) is comparable across methods — AQKA\-GP does not estimate the kernel more accurately overall, it estimates it more accurately*where it matters for the downstream loss*\.
Figure 7:Marginal\-likelihood approximation \(top\) and kernel Frobenius \(bottom\) errors vs\. shot budget, both on log\-log axes,ntr=200n\_\{tr\}=200,88seeds\. AQKA\-GP gives substantially lower NLL error than uniform across budgets, especially on dense data\. Kernel Frobenius error is comparable across methods — AQKA\-GP does not estimate𝐊\\mathbf\{K\}more accurately overall, only at the entries that matter forℒ\\mathcal\{L\}\.
#### Result 5: Bayesian\-optimization surrogate quality \(Figure[8](https://arxiv.org/html/2606.28833#Sx5.F8)\)\.
We test whether AQKA\-GP gives a better GP surrogate for Bayesian optimization\. We fit a GP tontr=120n\_\{tr\}=120uniformly sampled query points on Branin \(2D\), Hartmann\-3, and Hartmann\-6, then evaluateEI\\mathrm\{EI\}onncand=2000n\_\{\\mathrm\{cand\}\}=2000random candidates and pick the argmax\-EI point\. The simple regretf\(xpick\)−f∗f\(x\_\{\\mathrm\{pick\}\}\)\-f^\{\*\}over88seeds is reported\. AtB=2×105B=2\\times 10^\{5\}, AQKA\-GP\-marg gives substantial gain on Hartmann\-6 \(−25%\-25\\%regret\) and modest gain on Branin and Hartmann\-3\. Gains shrink at higher budgets, where uniform’s surrogate is already accurate enough to localize the argmax\-EI\. The Hartmann\-6 result is suggestive of a quantum\-BO use case where each shot is expensive: AQKA\-GP can find better candidates with the same per\-iteration budget\.
Figure 8:BO surrogate\-quality experiment,88seeds\. Simple regret of the argmax\-EI candidate after one GP fit onntr=120n\_\{tr\}=120points, under uniform vs\. two AQKA\-GP variants\. AQKA\-GP\-marg gives notable gain on Hartmann\-6 atB=2×105B=2\\times 10^\{5\}\(−25%\-25\\%regret\)\. At higher budgets, uniform’s surrogate is already accurate enough to localize the argmax\-EI, so gains shrink\.
#### Result 6: Full BO loop on Hartmann\-6 \(Figure[9](https://arxiv.org/html/2606.28833#Sx5.F9)\)\.
We move from surrogate quality to a full Bayesian\-optimization loop with shot\-budgeted GP refits at every iteration\. The setup is Hartmann\-6,ninit=30n\_\{\\mathrm\{init\}\}=30initial random queries,Biter=5×105B\_\{\\mathrm\{iter\}\}=5\\times 10^\{5\}shots per BO iteration,2020iterations,88seeds\. At each iteration the GP is refit from scratch with the chosen allocator, EI is maximized on10001000random candidates, and the picked point is evaluated\. Simple regret falls fastest under AQKA\-GP\-marg, reaching0\.830\.83at iteration2020vs\.1\.071\.07under uniform \(−23%\-23\\%\); AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|reaches0\.970\.97\(−9%\-9\\%\)\. The gap opens at iteration55–1010and persists, confirming that GP surrogate quality compounds across BO steps\.
Figure 9:Full BO loop on Hartmann\-6 with shot\-budgeted GP refit at each iteration \(Biter=5×105B\_\{\\mathrm\{iter\}\}=5\\times 10^\{5\},ninit=30n\_\{\\mathrm\{init\}\}=30,88seeds\)\. AQKA\-GP\-marg reaches simple regret0\.830\.83at iteration2020vs\.1\.071\.07for uniform \(−23%\-23\\%\)\. The gap opens around iteration55and persists, as surrogate quality compounds across BO steps\.
#### Result 7: Online streaming GP \(Figure[10](https://arxiv.org/html/2606.28833#Sx5.F10)\)\.
We streamT=60T=60data points online with per\-step shot budgetBstep=2×105B\_\{\\mathrm\{step\}\}=2\\times 10^\{5\}\. At each step the GP is refit from scratch using uniform or AQKA\-GP allocation, predicts the incoming point, and the running test MSE is updated\. AQKA\-GP\-marg gives the lowest cumulative MSE at every step, reaching0\.1510\.151att=60t=60vs\. uniform0\.1620\.162\(−7%\-7\\%\); AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|reaches0\.1570\.157\(−3%\-3\\%\)\. The online setting is the most stringent test of shot economy: every new data point requires a complete GP refit, so any per\-step savings compound\.
Figure 10:Online streaming GP with shot\-budgeted refit at every step \(Bstep=2×105B\_\{\\mathrm\{step\}\}=2\\times 10^\{5\},T=60T=60,66seeds\)\. Running test MSE: AQKA\-GP\-marg0\.1510\.151att=60t=60vs\. uniform0\.1620\.162\(−7%\-7\\%\)\. Online refits compound any per\-step savings\.
#### Result 8: Sparse VFE GP \(Figure[11](https://arxiv.org/html/2606.28833#Sx5.F11)\)\.
We test AQKA\-VFE on sparse inducing\-point GP regression \(Section[Extension to Sparse Inducing\-Point GPs](https://arxiv.org/html/2606.28833#Sx3.SSx4)\)\. Withntr=200n\_\{tr\}=200,m=30m=30inducing points \(uniform subsample from training\), and dense GP\-prior data, the kernel budget consists ofm\(m\+1\)/2\+nm=465\+6000=6465m\(m\+1\)/2\+nm=465\+6000=6465entries\. AQKA\-VFE beats uniform across all four budget multiples \(1010,5050,200200,10001000shots per entry\), with gains of−1\.1%\-1\.1\\%at low budget growing to−2\.2%\-2\.2\\%near oracle\. The gain is smaller in absolute terms than the full\-GP case because the sparse approximation already smooths predictive heterogeneity; in relative terms \(gap\-to\-oracle closure\), AQKA\-VFE closes the gap by∼\\sim40%40\\%atB=200neB=200n\_\{e\}\.
Figure 11:Sparse VFE GP \(ntr=200n\_\{tr\}=200,m=30m=30inducing points,66seeds,ne=6465n\_\{e\}=6465kernel entries\)\. AQKA\-VFE beats uniform across all four budget regimes \(−1\.1%\-1\.1\\%to−2\.2%\-2\.2\\%\), closing the gap to oracle by∼\\sim40%40\\%atB=200neB=200n\_\{e\}\. Smaller absolute gain than full\-GP because the sparse approximation already smooths predictive heterogeneity\.
#### Result 9:NN\-scaling \(Figure[12](https://arxiv.org/html/2606.28833#Sx5.F12)\)\.
We test whether AQKA\-GP gains hold across training\-set sizes\. WithB=50npairsB=50n\_\{\\mathrm\{pairs\}\}\(constant shots/pair, a budget level inside the sweet spot\), we varyntr∈\{50,100,200,400,600\}n\_\{tr\}\\in\\\{50,100,200,400,600\\\}on dense GP\-prior data,55seeds each\. At everyntrn\_\{tr\}at least one AQKA\-GP variant beats uniform, with gains ranging from−9%\-9\\%\(n=100n=100, gp\_loo\) to−30%\-30\\%\(n=200n=200, gp\_marg\)\. The best\-performing sensitivity varies withnn— gp\_alpha dominates atn∈\{50,400\}n\\in\\\{50,400\\\}, gp\_marg atn∈\{200,600\}n\\in\\\{200,600\\\}, gp\_loo atn=100n=100— but the existence of a winning AQKA\-GP variant is uniform\. The lack of monotone trend innnreflects that the per\-pair shot regime is held constant, so absolute RMSE depends primarily on the inducing posterior structure of each particular dataset\.
Figure 12:NN\-scaling atB=50npairsB=50n\_\{\\mathrm\{pairs\}\}shots/pair,55seeds\. At everyntr∈\{50,100,200,400,600\}n\_\{tr\}\\in\\\{50,100,200,400,600\\\}, at least one AQKA\-GP variant beats uniform; gains range from−9%\-9\\%\(n=100n=100\) to−30%\-30\\%\(n=200n=200, gp\_marg\)\. The best sensitivity varies withnn, but the existence of a winning variant is uniform\.
#### Result 10: Genuine quantum kernels \(Figure[13](https://arxiv.org/html/2606.28833#Sx5.F13)\)\.
The experiments above use a classical RBFKij=e−γ‖xi−xj‖2K\_\{ij\}=e^\{\-\\gamma\\\|x\_\{i\}\-x\_\{j\}\\\|^\{2\}\}corrupted by Bernoulli shot noise as a*controlled study*, isolating the shot\-allocation mechanism from circuit and device structure\. We now replace this with a genuine quantum kernelK\(xi,xj\)=\|⟨ϕ\(xi\)\|ϕ\(xj\)⟩\|2K\(x\_\{i\},x\_\{j\}\)=\|\\langle\\phi\(x\_\{i\}\)\|\\phi\(x\_\{j\}\)\\rangle\|^\{2\}evaluated by statevector simulation, plus a5%5\\%depolarizing channel that models device\-induced decay\. We run three complementary studies \(pairedtt\-tests,55seeds each\):
1. 1\.*UCI through a quantum kernel \(Figure[13](https://arxiv.org/html/2606.28833#Sx5.F13)left\)\.*We embed standardized UCI features into aqq\-qubit ZZ feature map \(q∈\{4,6,8\}\\in\\\{4,6,8\\\}, reps=2=2\), then run AQKA\-GP\. Onenergy,concrete,kin8nmatntr=60n\_\{tr\}=60, the gain is essentially null \(−4%\-4\\%to\+10%\+10\\%, almost no significance\), with a mild trend toward worse performance atq∈\{6,8\}q\\in\\\{6,8\\\}\. This is consistent with the exponential\-concentration regime ofThanasilpet al\.\([2024](https://arxiv.org/html/2606.28833#bib.bib16)\): at largerqqwith depolarizing noise, kernel values cluster near a fixed mean and per\-pair heterogeneity vanishes, leaving nothing for AQKA\-GP to exploit\. Real\-data feature embeddings into ZZ kernels are not, in general, in the regime where shot\-allocation matters; this is an honest negative finding\.
2. 2\.*Feature\-map sweep \(Figure[13](https://arxiv.org/html/2606.28833#Sx5.F13)center\)\.*On planted\-sparse quantum data atq=4q=4,ntr=50n\_\{tr\}=50, we compare four feature maps: ZZ\-full, ZZ\-linear, Pauli\-Z, Pauli\-Y \(each with reps=2=2,5%5\\%depolarizing\)\. At the low\-budget setting \(B=63,750B=63\{,\}750shots,5050per pair\), AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|delivers statistically significant gains on three of four feature maps: ZZ\-full−13\.3%\-13\.3\\%\(p=0\.035p=0\.035\), ZZ\-linear−14\.5%\-14\.5\\%\(p=0\.027p=0\.027\), Pauli\-Z−13\.3%\-13\.3\\%\(p=0\.035p=0\.035\); Pauli\-Y is marginal at−6\.3%\-6\.3\\%\(p=0\.55p=0\.55,n\.s\.n\.s\.\)\. Gains shrink at higher budget where the kernels reach the shot\-noise floor\.
3. 3\.*Scale sweep \(Figure[13](https://arxiv.org/html/2606.28833#Sx5.F13)right\)\.*On the ZZ\-full kernel atq=4q=4,ntr∈\{40,80,120\}n\_\{tr\}\\in\\\{40,80,120\\\}, AQKA\-GP gain rises withntrn\_\{tr\}:−1\.3%\-1\.3\\%\(n=40n=40,p=0\.50p=0\.50\),−16\.9%\-16\.9\\%\(n=80n=80,p=0\.076p=0\.076, marginal\),−9\.2%\-9\.2\\%\(n=120n=120,p=0\.48p=0\.48, n\.s\. in55seeds but visibly below uniform\)\. The trend agrees with Observation 1: heterogeneity grows with the spectral effective rank deficit, which in turn grows withntrn\_\{tr\}\.
The take\-away is honest\. AQKA\-GP works on quantum kernels when the underlying data has the heterogeneity the mechanism exploits \(Studies 2–3\); it does not magically improve a real\-data UCI regression task whose features were never designed for a ZZ kernel \(Study 1\)\.
Figure 13:Three quantum\-kernel studies,55seeds each,5%5\\%depolarizing noise\.Left: UCI features through ZZ feature map atq∈\{4,6,8\}q\\in\\\{4,6,8\\\}\(ntr=60n\_\{tr\}=60,55seeds\); gain is essentially null — the data is not in the regime AQKA\-GP exploits \(asterisks markp<0\.05p<0\.05pairedtt\-test\)\.Center: feature\-map sweep atq=4q=4, planted\-sparse data \(ntr=50n\_\{tr\}=50\); three of four maps give significant gain \(p<0\.05p<0\.05, marked with \*\) at low budgetB=64B=64k\.Right: scale sweep atq=4q=4, ZZ\-full \(ntr∈\{40,80,120\}n\_\{tr\}\\in\\\{40,80,120\\\}\); the gain grows withntrn\_\{tr\}, consistent with Observation 1\.
### Extensions: Bayesian Quadrature, Heteroscedastic GP, Cokriging
We test three further GP variants where the Neyman rule of Section[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3)applies with task\-specific sensitivities, each derived from∂\(task loss\)/∂Kij\\partial\(\\text\{task loss\}\)/\\partial K\_\{ij\}via the inverse identity\.
#### Result 11: Bayesian Quadrature \(Figure[14](https://arxiv.org/html/2606.28833#Sx5.F14)\)\.
We estimateI=𝔼x∼π\[f\(x\)\]I=\\mathbb\{E\}\_\{x\\sim\\pi\}\[f\(x\)\]withπ=𝒩\(0,I4\)\\pi=\\mathcal\{N\}\(0,I\_\{4\}\)andf∼𝒢𝒫\(0,kRBF\)f\\sim\\mathcal\{GP\}\(0,k\_\{\\mathrm\{RBF\}\}\)\. The posterior mean ofIIisI^=𝐳⊤\(𝐊\+σn2I\)−1𝐲\\widehat\{I\}=\\mathbf\{z\}^\{\\top\}\(\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}I\)^\{\-1\}\\mathbf\{y\}, where𝐳i=𝔼x∼π\[k\(xi,x\)\]\\mathbf\{z\}\_\{i\}=\\mathbb\{E\}\_\{x\\sim\\pi\}\[k\(x\_\{i\},x\)\]is the kernel\-mean embedding \(analytic for RBF×\\timesGaussian\)\. The induced sensitivitysensijbq=\|ziαj′\+zjαi′\|\\mathrm\{sens\}^\{\\mathrm\{bq\}\}\_\{ij\}=\|z\_\{i\}\\alpha^\{\\prime\}\_\{j\}\+z\_\{j\}\\alpha^\{\\prime\}\_\{i\}\|,𝜶′=A−1𝐳\\bm\{\\alpha\}^\{\\prime\}=A^\{\-1\}\\mathbf\{z\}, is highly heterogeneous: pairs near the prior mode dominate\. AQKA\-BQ gives\|I^−Ioracle\|\|\\widehat\{I\}\-I\_\{\\mathrm\{oracle\}\}\|of0\.0220\.022atB=2×107B=2\\\!\\times\\\!10^\{7\}vs\.0\.0350\.035for uniform \(−37%\-37\\%\)\. At low budget \(B≤106B\\leq 10^\{6\}\) all variants are dominated by warm\-up noise; in the moderate\-to\-high regime AQKA wins consistently\. The predictive\-coupling sensitivity\|αα\|\|\\alpha\\alpha\|also works because𝜶=A−1𝐲\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}correlates with𝜶′=A−1𝐳\\bm\{\\alpha\}^\{\\prime\}=A^\{\-1\}\\mathbf\{z\}at well\-conditionedAA\.
Figure 14:Bayesian Quadrature,n=200n=200,d=4d=4,π=𝒩\(0,I\)\\pi=\\mathcal\{N\}\(0,I\),88seeds\. Log\-log integral error\|I^−Ioracle\|\|\\widehat\{I\}\-I\_\{\\mathrm\{oracle\}\}\|\. AQKA\-BQ achieves0\.0220\.022atB=2×107B=2\\\!\\times\\\!10^\{7\}vs\. uniform0\.0350\.035\(−37%\-37\\%\)\.
#### Result 12: Heteroscedastic GP \(Figure[15](https://arxiv.org/html/2606.28833#Sx5.F15)\)\.
With input\-dependent noiseσ2\(x\)=σ02\(1\+0\.5‖x‖2\)\\sigma^\{2\}\(x\)=\\sigma\_\{0\}^\{2\}\(1\+0\.5\\\|x\\\|^\{2\}\), the posterior becomes𝜶=\(𝐊\+diag\(σ2\(x\)\)\)−1𝐲\\bm\{\\alpha\}=\(\\mathbf\{K\}\+\\mathrm\{diag\}\(\\sigma^\{2\}\(x\)\)\)^\{\-1\}\\mathbf\{y\}and pair sensitivities are no longer rank\-1 in a single𝜶\\bm\{\\alpha\}\. Empirically AQKA\-Hetero\-GP beats uniform substantially: atB=106B=10^\{6\},0\.4890\.489vs\. uniform0\.6900\.690\(−29%\-29\\%\); atB=2×107B=2\\\!\\times\\\!10^\{7\},0\.3600\.360vs\. uniform0\.4460\.446\(−19%\-19\\%\)\. The heteroscedasticity amplifies the gain because pairs in low\-noise regions are far more sensitive to kernel error than pairs in high\-noise regions\.
Figure 15:Heteroscedastic GP regression withσ2\(x\)=σ02\(1\+0\.5‖x‖2\)\\sigma^\{2\}\(x\)=\\sigma\_\{0\}^\{2\}\(1\+0\.5\\\|x\\\|^\{2\}\),66seeds\. AQKA\-Hetero\-GP gains−19%\-19\\%to−29%\-29\\%across budgets, larger than the homoscedastic case because input\-dependent noise concentrates sensitivity in low\-noise regions\.
#### Result 13: Hyperparameter learning loop \(Figure[16](https://arxiv.org/html/2606.28833#Sx5.F16)\)\.
The marginal\-likelihood sensitivity \([15](https://arxiv.org/html/2606.28833#Sx3.E15)\) was motivated by hyperparameter learning\. Here we test it directly: gradient descent on the kernel bandwidthγ\\gammain log\-space, starting atγinit=0\.04\\gamma\_\{\\mathrm\{init\}\}=0\.04and aiming forγtrue=0\.15\\gamma\_\{\\mathrm\{true\}\}=0\.15, with shot budgetBiter=2×105B\_\{\\mathrm\{iter\}\}=2\\\!\\times\\\!10^\{5\}per gradient step andT=20T=20steps\. The marginal\-likelihood gradient is computed from the shot\-noisy𝐊^\\widehat\{\\mathbf\{K\}\}under each allocator and used to updatelogγ\\log\\gammawith learning rate0\.050\.05\. AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|closes the gap toγtrue\\gamma\_\{\\mathrm\{true\}\}by65%65\\%\(\|γ−γtrue\|=0\.086\|\\gamma\-\\gamma\_\{\\mathrm\{true\}\}\|=0\.086vs\. uniform0\.2460\.246\); AQKA\-GP\-marg by43%43\\%\(0\.1410\.141\)\. Final test RMSEs are0\.4040\.404,0\.4110\.411,0\.4190\.419respectively \(−4%\-4\\%for AQKA\-alpha\), confirming that shot\-efficient hyperparameter learning translates into a downstream\-test improvement\.
Figure 16:Outer\-loop hyperparameter learning: gradient descent onlogγ\\log\\gammaunder shot\-budgeted marginal\-likelihood gradient,T=20T=20iterations,Biter=2×105B\_\{\\mathrm\{iter\}\}=2\\\!\\times\\\!10^\{5\},66seeds\. \(Left\)γ\\gammatrajectory; AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|converges closest toγtrue=0\.15\\gamma\_\{\\mathrm\{true\}\}=0\.15\(gap0\.0860\.086vs\. uniform0\.2460\.246,−65%\-65\\%\)\. \(Right\) Final test RMSE; AQKA\-GP\-\|αα\|\|\\alpha\\alpha\|improves by−4%\-4\\%\.
#### Result 14: Multi\-output GP / Cokriging \(Figure[17](https://arxiv.org/html/2606.28833#Sx5.F17)\)\.
Forq=3q=3correlated outputs with separable kernelcov\(yt\(x\),yt′\(x′\)\)=B\[t,t′\]k\(x,x′\)\\mathrm\{cov\}\(y\_\{t\}\(x\),y\_\{t^\{\\prime\}\}\(x^\{\\prime\}\)\)=B\[t,t^\{\\prime\}\]k\(x,x^\{\\prime\}\), the posterior decomposes in the eigenbasis ofBB\. We use AQKA\-MT\-GP with sensitivity∑tλt\|𝜶t𝜶t⊤\|\\sum\_\{t\}\\lambda\_\{t\}\|\\bm\{\\alpha\}\_\{t\}\\bm\{\\alpha\}\_\{t\}^\{\\top\}\|summed over outputs\. AQKA\-MT\-GP gains−4%\-4\\%at low budget growing to tied at oracle floor\. The gain is smaller than full\-GP because the eigenbasis decomposition averages sensitivity across outputs, partially homogenizing the pair\-level weights\.
Figure 17:Multi\-output GP / Cokriging with separable kernelB⊗KB\\otimes K,q=3q=3correlated outputs,66seeds\. AQKA\-MT\-GP gains−4%\-4\\%to tied across budgets\. Gains are smaller than full\-GP because the output\-eigenbasis average homogenizes pair\-level weights\.
## From Theory to Experiment
The four propositions of Section[Theory](https://arxiv.org/html/2606.28833#Sx4)make sharp predictions that we can read off the experimental panels:
- •*Observation 1 \(Lemma[1](https://arxiv.org/html/2606.28833#Thmlemma1)\) predicts dense\>\>sparse*: GP regression manufactures predictive heterogeneity fromA−1A^\{\-1\}even when the data itself is homogeneous, so AQKA\-GP should win more on dense data than on planted\-sparse\. Figure[1](https://arxiv.org/html/2606.28833#Sx5.F1)confirms this \(−21%\-21\\%on dense,−11%\-11\\%on planted\-sparse\), and the four UCI datasets \(all in the dense regime\) match the predicted−12%\-12\\%to−17%\-17\\%band\.
- •*Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)predicts a coverage breakpoint atρ≈n2/\(2B\)\\rho\\approx n^\{2\}/\(2B\)*: forn=200n=200andB=2×105B=2\\times 10^\{5\}this isρ≈0\.1\\rho\\approx 0\.1\. Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)matches this:ρ≤0\.1\\rho\\leq 0\.1is catastrophic,ρ≥0\.2\\rho\\geq 0\.2is recovered\. The further empirical jump toρ=0\.5\\rho=0\.5is a shot\-quality margin not derived from the proposition \(which is a coverage / necessary\-condition result, not a sufficient\-condition one\)\.
- •*Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)predicts1/σn41/\\sigma\_\{n\}^\{4\}amplification*: GP regression should be more shot\-hungry than classification, requiring higher per\-pair shot counts to reach the same accuracy floor\. The minimum budget at which AQKA\-GP starts to help is∼50\\sim 50shots/pair \(Figure[1](https://arxiv.org/html/2606.28833#Sx5.F1)and Figure[5](https://arxiv.org/html/2606.28833#Sx5.F5)\), substantially higher than the∼4\\sim 4shots/pair regime where classification AQKA helps\(Xuet al\.[2026](https://arxiv.org/html/2606.28833#bib.bib1)\)\.
- •*Adaptiven\\sqrt\{n\}jitter*: kernel Frobenius errors are similar across allocators \(Figure[7](https://arxiv.org/html/2606.28833#Sx5.F7), bottom row\), confirming that AQKA\-GP does not estimate𝐊\\mathbf\{K\}more accurately overall — it estimates it more accurately*where the loss is sensitive to it*\. The NLL approximation error \(top row\) is correspondingly smaller for AQKA\-GP, with then\\sqrt\{n\}jitter ensuring the posterior remains stable across shot regimes\.
The picture is consistent: AQKA\-GP works in regimes the theory predicts, fails in regimes the theory predicts, and the magnitude of gain matches the spectral\-decay heuristic from Observation 1 \(Lemma[1](https://arxiv.org/html/2606.28833#Thmlemma1)\)\.
## Related Work
#### Quantum kernel methods\.
Quantum kernels were introduced byHavlíčeket al\.\([2019](https://arxiv.org/html/2606.28833#bib.bib4)\)andSchuld and Killoran \([2019](https://arxiv.org/html/2606.28833#bib.bib13)\)as a quantum analogue of classical kernel methods\.Huanget al\.\([2021](https://arxiv.org/html/2606.28833#bib.bib6)\)established sample\-complexity lower bounds, andThanasilpet al\.\([2024](https://arxiv.org/html/2606.28833#bib.bib16)\)characterized exponential concentration regimes where shot noise overwhelms signal\.
#### Shot\-efficient quantum kernel estimation\.
The closest prior work is the AQKA framework ofXuet al\.\([2026](https://arxiv.org/html/2606.28833#bib.bib1)\), which applies Neyman\-style allocation to quantum KRR and SVM, and the concurrent work ofMiroszewski \([2026](https://arxiv.org/html/2606.28833#bib.bib7)\), which independently derived an adaptive measurement allocation for kernelized SVMs under noisy observations\. Both focus on classification\. We extend the line to GP regression, identify the GP\-specific high\-floor requirement, and derive three sensitivities that are unique to the GP setting \(marginal likelihood, LOO residual, predictive coupling\)\. Other related directions include kernel\-bandwidth tuning for quantum machine learning\(Shaydulin and Wild[2022](https://arxiv.org/html/2606.28833#bib.bib14)\), Bayesian deep learning compiled to quantum circuits\(Zhaoet al\.[2019](https://arxiv.org/html/2606.28833#bib.bib15)\), and quantum active learning for materials design\(Lourençoet al\.[2026](https://arxiv.org/html/2606.28833#bib.bib19)\)\.
#### Classical sensitivity\-based subsampling\.
For classical Gaussian processes, leverage\-score subsampling\(Calandrielloet al\.[2017](https://arxiv.org/html/2606.28833#bib.bib2)\)and ridge\-leverage\-driven coreset construction\(Musco and Musco[2017](https://arxiv.org/html/2606.28833#bib.bib8)\)are the closest analogue\. These methods choose*which rows*of𝐊\\mathbf\{K\}to compute exactly, whereas we choose*how many shots*each entry of𝐊\\mathbf\{K\}receives\. The two ideas can be composed — AQKA\-GP plus row subselection is an obvious follow\-up but not the focus here\.
#### A\-optimal experimental design\.
Allocating samples by\|∂ℒ/∂θ\|Var\|\\partial\\mathcal\{L\}/\\partial\\theta\|\\sqrt\{\\mathrm\{Var\}\}is classical\(Neyman[1992](https://arxiv.org/html/2606.28833#bib.bib9); Pukelsheim[2006](https://arxiv.org/html/2606.28833#bib.bib10)\)\. We follow this lineage in the kernel\-entry setting\.
## Discussion and Limitations
#### Hyperparameter learning\.
The marginal\-likelihood sensitivity is in principle suited to outer\-loop hyperparameter learning \(gradient descent onθ\\thetavia∂ℒ/∂θ=∑ij\(∂ℒ/∂Kij\)\(∂Kij/∂θ\)\\partial\\mathcal\{L\}/\\partial\\theta=\\sum\_\{ij\}\(\\partial\\mathcal\{L\}/\\partial K\_\{ij\}\)\(\\partial K\_\{ij\}/\\partial\\theta\)\)\. In our experiments we fixσn\\sigma\_\{n\}and the kernel widthγ\\gammaand study the inner\-loop shot allocation problem\. Composing AQKA\-GP with hyperparameter learning is left to future work\.
#### Bayesian optimization\.
The acquisition functions used in Bayesian optimization \(EI, UCB, PI\) depend onμ∗\\mu\_\{\*\}andσ∗2\\sigma\_\{\*\}^\{2\}in closed form, and their gradients with respect toKijK\_\{ij\}admit closed forms via the same identities used here\. Applying AQKA\-GP to Bayesian optimization with quantum kernels is a natural next step; the per\-iteration shot budget is precisely the bottleneck practitioners report\.
#### Catastrophic low\-budget regime\.
The low\-budget catastrophe is genuine\. Possible remedies include \(i\) Nyström\-style landmark warm\-up that estimates a low\-rank𝐊^\\widehat\{\\mathbf\{K\}\}before sensitivity is computed, \(ii\) Bayesian priors onK^ij\\widehat\{K\}\_\{ij\}that pull toward a smoothness model rather than the Bernoulli prior mean, and \(iii\) deferred allocation that runs multiple rounds of sensitivity re\-estimation\. We expect a fully\-online AQKA\-GP loop to extend the useful regime by an order of magnitude inBB\.
## Conclusion
We extended the AQKA framework from classification to Gaussian process regression\. The extension is non\-trivial: GP regression cares about quantities \(full\-spectrum posterior variance, log\-determinant ofAA, marginal likelihood\) that classification’s0/10/1accuracy averages away, and naive sensitivity\-driven allocation fails catastrophically without a high uniform\-coverage floor \(Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)\)\. The three pair\-level sensitivities we propose are principled:\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|is the rank\-1 specialization of the exact predictive\-MSE Neyman weight \(atM=𝐲𝐲⊤M=\\mathbf\{y\}\\mathbf\{y\}^\{\\top\}\), the marginal\-likelihood gradient is the exact partial derivative, and the LOO sensitivity drops a controllable correction term \(Propositions[4](https://arxiv.org/html/2606.28833#Thmproposition4)–[6](https://arxiv.org/html/2606.28833#Thmproposition6)\)\. Combined with a coverage floor \(justified by the necessary condition of Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)\) and an inference\-time jitter whosen\\sqrt\{n\}scaling is motivated by Wigner\-type matrix concentration but whose constant is calibrated empirically, AQKA\-GP delivers1010–21%21\\%RMSE improvement across four UCI benchmarks and two synthetic settings, with the gain transferring \(i\) to genuine ZZ and Pauli\-Z quantum kernels at the budget regimes where heterogeneity exists \(−13\-13–15%15\\%,p<0\.05p<0\.05, Result 10\) and \(ii\) to four downstream tasks \(Bayesian quadrature, heteroscedastic regression, hyperparameter learning, multi\-output Cokriging\)\. On UCI features embedded into a ZZ kernel \(Result 10, Study 1\) the gain disappears, consistent with the exponential\-concentration regime ofThanasilpet al\.\([2024](https://arxiv.org/html/2606.28833#bib.bib16)\)— a regime where shot allocation has nothing to exploit\.
## References
- D\. Calandriello, A\. Lazaric, and M\. Valko \(2017\)Distributed adaptive sampling for kernel matrix approximation\.InArtificial Intelligence and Statistics,pp\. 1421–1429\.Cited by:[Classical sensitivity\-based subsampling\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px3.p1.2)\.
- A\. Damianou and N\. D\. Lawrence \(2013\)Deep gaussian processes\.InArtificial intelligence and statistics,pp\. 207–215\.Cited by:[Extension to Deep GPs \(DSVI\)\.](https://arxiv.org/html/2606.28833#Sx3.SSx4.SSS0.Px5.p1.4)\.
- V\. Havlíček, A\. D\. Córcoles, K\. Temme, A\. W\. Harrow, A\. Kandala, J\. M\. Chow, and J\. M\. Gambetta \(2019\)Supervised learning with quantum\-enhanced feature spaces\.Nature567\(7747\),pp\. 209–212\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p1.6),[Quantum Kernel Estimation as Bernoulli Sampling](https://arxiv.org/html/2606.28833#Sx2.SSx1.p1.9),[Quantum kernel methods\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px1.p1.1)\.
- J\. Hensman, N\. Fusi, and N\. D\. Lawrence \(2013\)Gaussian processes for big data\.InProceedings of the Twenty\-Ninth Conference on Uncertainty in Artificial Intelligence,pp\. 282–290\.Cited by:[\(b\) Variational SVGP\(Hensmanet al\.2013\)\.](https://arxiv.org/html/2606.28833#Sx3.SSx4.SSS0.Px2)\.
- H\. Huang, M\. Broughton, M\. Mohseni, R\. Babbush, S\. Boixo, H\. Neven, and J\. R\. McClean \(2021\)Power of data in quantum machine learning\.Nature communications12\(1\),pp\. 2631\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p1.6),[Quantum kernel methods\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px1.p1.1)\.
- M\. P\. Lourenço, H\. Zadeh\-Haghighi, J\. Hostaš, M\. Naseri, D\. Gaur, C\. Simon, and D\. R\. Salahub \(2026\)Exploring quantum active learning for materials design and discovery\.Physical Chemistry Chemical Physics28\(14\),pp\. 8604–8617\.Cited by:[Shot\-efficient quantum kernel estimation\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px2.p1.1)\.
- A\. Miroszewski \(2026\)Adaptive measurement allocation for learning kernelized svms under noisy observations\.arXiv preprint arXiv:2605\.22275\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p2.3),[What makes GP shot\-hungrier than classification?](https://arxiv.org/html/2606.28833#Sx2.SSx2.SSS0.Px1.p1.7),[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3.p1.8),[Shot\-efficient quantum kernel estimation\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px2.p1.1)\.
- C\. Musco and C\. Musco \(2017\)Recursive sampling for the nystrom method\.Advances in neural information processing systems30\.Cited by:[Classical sensitivity\-based subsampling\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px3.p1.2)\.
- J\. Neyman \(1992\)On the two different aspects of the representative method: the method of stratified sampling and the method of purposive selection\.InBreakthroughs in statistics: Methodology and distribution,pp\. 123–150\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p2.3),[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3.p1.8),[A\-optimal experimental design\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px4.p1.1)\.
- F\. Pukelsheim \(2006\)Optimal design of experiments\.SIAM\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p2.3),[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3.p1.8),[A\-optimal experimental design\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px4.p1.1)\.
- C\. E\. Rasmussen \(2003\)Gaussian processes in machine learning\.InSummer school on machine learning,pp\. 63–71\.Cited by:[\(iii\) Leave\-one\-out residual for predictive calibration\.](https://arxiv.org/html/2606.28833#Sx3.SSx1.SSS0.Px3.p1.6)\.
- H\. Salimbeni and M\. Deisenroth \(2017\)Doubly stochastic variational inference for deep gaussian processes\.Advances in neural information processing systems30\.Cited by:[Extension to Deep GPs \(DSVI\)\.](https://arxiv.org/html/2606.28833#Sx3.SSx4.SSS0.Px5.p1.4)\.
- M\. Schuld and N\. Killoran \(2019\)Quantum machine learning in feature hilbert spaces\.Physical review letters122\(4\),pp\. 040504\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p1.6),[Quantum kernel methods\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px1.p1.1)\.
- R\. Shaydulin and S\. M\. Wild \(2022\)Importance of kernel bandwidth in quantum machine learning\.Physical Review A106\(4\),pp\. 042407\.Cited by:[Shot\-efficient quantum kernel estimation\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px2.p1.1)\.
- E\. Snelson and Z\. Ghahramani \(2005\)Sparse gaussian processes using pseudo\-inputs\.Advances in neural information processing systems18\.Cited by:[\(c\) FITC\(Snelson and Ghahramani2005\)\.](https://arxiv.org/html/2606.28833#Sx3.SSx4.SSS0.Px3)\.
- S\. Thanasilp, S\. Wang, M\. Cerezo, and Z\. Holmes \(2024\)Exponential concentration in quantum kernel methods\.Nature communications15\(1\),pp\. 5200\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p1.6),[Exponential concentration regime\.](https://arxiv.org/html/2606.28833#Sx2.SSx1.SSS0.Px1.p1.5),[item 1](https://arxiv.org/html/2606.28833#Sx5.I3.i1.p1.8),[Quantum kernel methods\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px1.p1.1),[Conclusion](https://arxiv.org/html/2606.28833#Sx9.p1.10)\.
- M\. Titsias \(2009\)Variational learning of inducing variables in sparse gaussian processes\.InArtificial intelligence and statistics,pp\. 567–574\.Cited by:[\(a\) Sparse VFE / SGPR\(Titsias2009\)\.](https://arxiv.org/html/2606.28833#Sx3.SSx4.SSS0.Px1)\.
- J\. A\. Tropp \(2015\)An introduction to matrix concentration inequalities\.Foundations and trends® in machine learning8\(1\-2\),pp\. 1–230\.Cited by:[Whyn⋅K\(1−K\)/s¯\\sqrt\{n\\cdot\\overline\{K\(1\-K\)/s\}\}jitter\.](https://arxiv.org/html/2606.28833#Sx4.SSx4.SSS0.Px1.p1.3)\.
- J\. Xu, C\. Li, D\. Zeng, J\. Paisley, and Q\. Zhao \(2026\)AQKA: active quantum kernel acquisition under a shot budget\.arXiv preprint arXiv:2605\.14672\.Cited by:[Introduction](https://arxiv.org/html/2606.28833#Sx1.p2.3),[What makes GP shot\-hungrier than classification?](https://arxiv.org/html/2606.28833#Sx2.SSx2.SSS0.Px1.p1.7),[Neyman Minimum\-Variance Shot Allocation](https://arxiv.org/html/2606.28833#Sx2.SSx3.p1.8),[\(i\) Predictive coupling\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|for test RMSE\.](https://arxiv.org/html/2606.28833#Sx3.SSx1.SSS0.Px1.p1.9),[Figure 6](https://arxiv.org/html/2606.28833#Sx5.F6),[Synthetic settings\.](https://arxiv.org/html/2606.28833#Sx5.SSx5.SSS0.Px2.p1.7),[Result 3: floor ablation \(Figure6\)\.](https://arxiv.org/html/2606.28833#Sx5.SSx5.SSS0.Px7.p1.18),[3rd item](https://arxiv.org/html/2606.28833#Sx6.I4.i3.p1.3),[Shot\-efficient quantum kernel estimation\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px2.p1.1)\.
- Z\. Zhao, A\. Pozas\-Kerstjens, P\. Rebentrost, and P\. Wittek \(2019\)Bayesian deep learning on a quantum computer\.Quantum Machine Intelligence1\(1\),pp\. 41–51\.Cited by:[Shot\-efficient quantum kernel estimation\.](https://arxiv.org/html/2606.28833#Sx7.SSx1.SSS0.Px2.p1.1)\.
## Appendix AAppendix A: Proofs
### Proof of Proposition[1](https://arxiv.org/html/2606.28833#Thmproposition1)\(Neyman minimum\-variance allocation\)
We minimizeΦ\(𝒔\)=∑i≤jcij2/sij\\Phi\(\\bm\{s\}\)=\\sum\_\{i\\leq j\}c\_\{ij\}^\{2\}/s\_\{ij\},cij:=\|gij\|Kij\(1−Kij\)c\_\{ij\}:=\|g\_\{ij\}\|\\sqrt\{K\_\{ij\}\(1\-K\_\{ij\}\)\}, subject to∑sij=B\\sum s\_\{ij\}=Bandsij≥0s\_\{ij\}\\geq 0\. The Lagrangian isℒLag=∑cij2/sij\+λ\(∑sij−B\)\\mathcal\{L\}\_\{\\mathrm\{Lag\}\}=\\sum c\_\{ij\}^\{2\}/s\_\{ij\}\+\\lambda\(\\sum s\_\{ij\}\-B\)with∂ℒLag/∂sij=−cij2/sij2\+λ=0\\partial\\mathcal\{L\}\_\{\\mathrm\{Lag\}\}/\\partial s\_\{ij\}=\-c\_\{ij\}^\{2\}/s\_\{ij\}^\{2\}\+\\lambda=0, givingsij=cij/λs\_\{ij\}=c\_\{ij\}/\\sqrt\{\\lambda\}\. Substituting into the constraint yieldsλ=∑cij/B\\sqrt\{\\lambda\}=\\sum c\_\{ij\}/Band thussij∗=Bcij/∑ckℓs^\{\*\}\_\{ij\}=B\\,c\_\{ij\}/\\sum c\_\{k\\ell\}\. The minimum value is
Φ\(𝒔∗\)=∑cij2⋅∑ckℓBcij=\(∑cij\)2B\.\\Phi\(\\bm\{s\}^\{\*\}\)=\\sum c\_\{ij\}^\{2\}\\cdot\\frac\{\\sum c\_\{k\\ell\}\}\{B\\,c\_\{ij\}\}=\\frac\{\(\\sum c\_\{ij\}\)^\{2\}\}\{B\}\.Strict convexity of1/s1/sovers\>0s\>0guarantees uniqueness\. ∎
### Proof of Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)\(Posterior error\)
We use the resolvent identityA^−1−A−1=−A^−1ΔA−1\\widehat\{A\}^\{\-1\}\-A^\{\-1\}=\-\\widehat\{A\}^\{\-1\}\\Delta A^\{\-1\}\. By submultiplicativity of the operator norm,
‖A^−1−A−1‖op≤‖A^−1‖op⋅‖Δ‖op⋅‖A−1‖op\.\\\|\\widehat\{A\}^\{\-1\}\-A^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\leq\\\|\\widehat\{A\}^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\cdot\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\cdot\\\|A^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\.The spectrum ofAAlies in\[σn2,σn2\+‖𝐊‖op\]\[\\sigma\_\{n\}^\{2\},\\sigma\_\{n\}^\{2\}\+\\\|\\mathbf\{K\}\\\|\_\{\\mathrm\{op\}\}\], so‖A−1‖op≤1/σn2\\\|A^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\leq 1/\\sigma\_\{n\}^\{2\}\. By Weyl’s inequality, the spectrum ofA^\\widehat\{A\}is shifted by at most‖Δ‖op\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}, so‖A^−1‖op≤1/\(σn2−‖Δ‖op\)\\\|\\widehat\{A\}^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\leq 1/\(\\sigma\_\{n\}^\{2\}\-\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\)whenever the denominator is positive\. Combining,
‖A^−1−A−1‖op≤‖Δ‖opσn2\(σn2−‖Δ‖op\)\.\\\|\\widehat\{A\}^\{\-1\}\-A^\{\-1\}\\\|\_\{\\mathrm\{op\}\}\\leq\\frac\{\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\}\{\\sigma\_\{n\}^\{2\}\(\\sigma\_\{n\}^\{2\}\-\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\)\}\.The𝜶\\bm\{\\alpha\}bound follows from𝜶^−𝜶=\(A^−1−A−1\)𝐲\\widehat\{\\bm\{\\alpha\}\}\-\\bm\{\\alpha\}=\(\\widehat\{A\}^\{\-1\}\-A^\{\-1\}\)\\mathbf\{y\}and‖𝐲‖≤M\\\|\\mathbf\{y\}\\\|\\leq M\. Theμ∗\\mu\_\{\*\}bound follows fromμ^∗−μ∗=𝐤∗⊤\(𝜶^−𝜶\)\\widehat\{\\mu\}\_\{\*\}\-\\mu\_\{\*\}=\\mathbf\{k\}\_\{\*\}^\{\\top\}\(\\widehat\{\\bm\{\\alpha\}\}\-\\bm\{\\alpha\}\)and Cauchy\-Schwarz on‖𝐤∗‖≤M\\\|\\mathbf\{k\}\_\{\*\}\\\|\\leq M\. ∎
### Proof of Lemma[1](https://arxiv.org/html/2606.28833#Thmlemma1)\(Covariance of𝜶\\bm\{\\alpha\}under random labels\)
In the eigenbasis of𝐊\\mathbf\{K\},𝜶=A−1𝐲=∑kβk𝐮k/\(λk\+σn2\)\\bm\{\\alpha\}=A^\{\-1\}\\mathbf\{y\}=\\sum\_\{k\}\\beta\_\{k\}\\,\\mathbf\{u\}\_\{k\}/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)sinceA𝐮k=\(λk\+σn2\)𝐮kA\\mathbf\{u\}\_\{k\}=\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)\\mathbf\{u\}\_\{k\}\. The covariance is
𝔼\[𝜶𝜶⊤\]=τ2∑k𝐮k𝐮k⊤\(λk\+σn2\)2,\\mathbb\{E\}\[\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}\]=\\tau^\{2\}\\,\\sum\_\{k\}\\frac\{\\mathbf\{u\}\_\{k\}\\mathbf\{u\}\_\{k\}^\{\\top\}\}\{\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{2\}\},which we read off diagonally asVar\(αi\)=τ2∑kuki2/\(λk\+σn2\)2\\mathrm\{Var\}\(\\alpha\_\{i\}\)=\\tau^\{2\}\\sum\_\{k\}u\_\{ki\}^\{2\}/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{2\}\. Squaring and taking expectations on theℓ2\\ell\_\{2\}andℓ1\\ell\_\{1\}norms of𝜶𝜶⊤\\bm\{\\alpha\}\\bm\{\\alpha\}^\{\\top\}gives the bound stated in the proposition body\. The dominant mechanism is the quartic weighting1/\(λk\+σn2\)41/\(\\lambda\_\{k\}\+\\sigma\_\{n\}^\{2\}\)^\{4\}on the small\-eigenvalue tail, multiplied by the unit\-norm constraint∑kuki2=1\\sum\_\{k\}u\_\{ki\}^\{2\}=1\. We do*not*assume eigenvector localization; in fact for translation\-invariant RBF on\[0,1\]d\[0,1\]^\{d\}, low\-λk\\lambda\_\{k\}eigenvectors are delocalized \(Fourier\-like\)\. The concentration of\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|arises from finite\-sample fluctuations ofuki2u\_\{ki\}^\{2\}around1/n1/nbeing amplified by the same quartic factor, not from any singleuki2u\_\{ki\}^\{2\}being intrinsically large\. The constantCCin the proposition body depends on the spectral effective rankr¯=∑k1/\(1\+σn2/λk\)2\\bar\{r\}=\\sum\_\{k\}1/\(1\+\\sigma\_\{n\}^\{2\}/\\lambda\_\{k\}\)^\{2\}; for RBF and ZZ feature\-map kernels at the bandwidths we use,r¯≪n\\bar\{r\}\\ll n, and the bound is non\-vacuous\. ∎
### Proof of Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)\(Catastrophic\-regime bound\)
The default substitutionK^ij=c\\widehat\{K\}\_\{ij\}=conSSintroduces an error matrixΔ=∑\(i,j\)∈Sδij\(eiej⊤\+ejei⊤\)\\Delta=\\sum\_\{\(i,j\)\\in S\}\\delta\_\{ij\}\(e\_\{i\}e\_\{j\}^\{\\top\}\+e\_\{j\}e\_\{i\}^\{\\top\}\)withδij=c−Kij\\delta\_\{ij\}=c\-K\_\{ij\}\. By assumption, at least anη\\eta\-fraction of pairs inSShave\|δij\|≥d\|\\delta\_\{ij\}\|\\geq d\. The Frobenius norm is therefore bounded*below*by
‖Δ‖F2=2∑\(i,j\)∈Sδij2≥2ηd2\|S\|,\\\|\\Delta\\\|\_\{F\}^\{2\}\\;=\\;2\\sum\_\{\(i,j\)\\in S\}\\delta\_\{ij\}^\{2\}\\;\\geq\\;2\\eta d^\{2\}\|S\|,giving‖Δ‖F≥d2η\|S\|\\\|\\Delta\\\|\_\{F\}\\geq d\\sqrt\{2\\eta\|S\|\}\. For any realn×nn\\times nsymmetric matrix the operator norm and Frobenius norm satisfy‖Δ‖op≥‖Δ‖F/n\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\geq\\\|\\Delta\\\|\_\{F\}/\\sqrt\{n\}\. Combining,
‖Δ‖op≥‖Δ‖Fn≥d2η\|S\|n\.\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}\\;\\geq\\;\\frac\{\\\|\\Delta\\\|\_\{F\}\}\{\\sqrt\{n\}\}\\;\\geq\\;d\\sqrt\{\\frac\{2\\eta\|S\|\}\{n\}\}\.This is a lower bound on the perturbation magnitude only\. We do*not*chain it through Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)to claim a lower bound on the predictive error: Proposition[2](https://arxiv.org/html/2606.28833#Thmproposition2)is an upper bound on the predictive error in terms of‖Δ‖op\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}, and combining an upper bound with a lower bound on the input does not yield a lower bound on the output \(the perturbationΔ\\Deltacould be orthogonal to the sensitive subspace ofA−1𝐲A^\{\-1\}\\mathbf\{y\}\)\. The catastrophic regime of Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)is reflected here only through the*necessary condition*that vanishing predictive error requires vanishing‖Δ‖op\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}, which in turn requires\|S\|→0\|S\|\\to 0at fixedd,ηd,\\eta\. ∎
#### What is proved, and what is not\.
Propositions[1](https://arxiv.org/html/2606.28833#Thmproposition1),[2](https://arxiv.org/html/2606.28833#Thmproposition2),[4](https://arxiv.org/html/2606.28833#Thmproposition4)–[6](https://arxiv.org/html/2606.28833#Thmproposition6)and Lemma[1](https://arxiv.org/html/2606.28833#Thmlemma1)are sharp up to constants\. Proposition[3](https://arxiv.org/html/2606.28833#Thmproposition3)is a*necessary*condition for predictive consistency \(it bounds‖Δ‖op\\\|\\Delta\\\|\_\{\\mathrm\{op\}\}from below in terms of the missing\-coverage fraction\),*not*a lower bound on the predictive error itself; lower bounds on the predictive error would require the perturbation to align with theA−1𝐲A^\{\-1\}\\mathbf\{y\}sensitive direction, which we do not prove here\. The empiricalρ=0\.5\\rho=0\.5floor is calibrated against Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)as a shot\-quality margin on top of this coverage condition\. Observation 1 is descriptive rather than quantitative and is tested empirically through realized AQKA\-GP gain on RBF and ZZ kernels\.
## Appendix BAppendix B: Experimental Setup
We collect the per\-experiment settings here for reproducibility\. Across all experiments the simulator kernel isKij=exp\(−γ‖xi−xj‖2\)K\_\{ij\}=\\exp\(\-\\gamma\\\|x\_\{i\}\-x\_\{j\}\\\|^\{2\}\), the shot\-noise process samplesK^ij⋅sij∼Bin\(sij,Kij\)\\widehat\{K\}\_\{ij\}\\cdot s\_\{ij\}\\sim\\mathrm\{Bin\}\(s\_\{ij\},K\_\{ij\}\), the warm\-up fraction isρw=0\.1\\rho\_\{w\}=0\.1, the uniform floor isρ=0\.5\\rho=0\.5, and the inference\-time jitter isj=n⋅K\(1−K\)/s¯j=\\sqrt\{n\}\\cdot\\overline\{K\(1\-K\)/s\}\. Adam\-style adaptive step sizes are not used; we report a single final allocation per seed\.
#### B\.1 Synthetic settings \(Figure[1](https://arxiv.org/html/2606.28833#Sx5.F1)\)\.
*Dense*:xi∼𝒩\(0,I6\)x\_\{i\}\\sim\\mathcal\{N\}\(0,I\_\{6\}\),f∼𝒢𝒫\(0,KRBF\(γ=0\.1\)\)f\\sim\\mathcal\{GP\}\(0,K\_\{\\mathrm\{RBF\}\}\(\\gamma\{=\}0\.1\)\)drawn via Cholesky,yi=f\(xi\)\+εiy\_\{i\}=f\(x\_\{i\}\)\+\\varepsilon\_\{i\}withεi∼𝒩\(0,σn2\)\\varepsilon\_\{i\}\\sim\\mathcal\{N\}\(0,\\sigma\_\{n\}^\{2\}\),σn=0\.3\\sigma\_\{n\}=0\.3\.*Planted\-sparse*: samexxdistribution,ma=15m\_\{a\}=15random anchor indices𝒜⊂\{1,…,n\}\\mathcal\{A\}\\subset\\\{1,\\dots,n\\\}withca∼𝒩\(0,1\)c\_\{a\}\\sim\\mathcal\{N\}\(0,1\)on𝒜\\mathcal\{A\},f\(x\)=∑i∈𝒜ciK\(x,xi\)f\(x\)=\\sum\_\{i\\in\\mathcal\{A\}\}c\_\{i\}K\(x,x\_\{i\}\)\. Both:ntr=200n\_\{tr\}=200,nte=80n\_\{te\}=80,1010seeds, budgetsB∈\{2×105,106,5×106,2×107\}B\\in\\\{2\\\!\\times\\\!10^\{5\},10^\{6\},5\\\!\\times\\\!10^\{6\},2\\\!\\times\\\!10^\{7\}\\\}\.
#### B\.2 UCI benchmarks \(Figure[5](https://arxiv.org/html/2606.28833#Sx5.F5), Table[1](https://arxiv.org/html/2606.28833#Sx5.T1)\)\.
We use four standard UCI regression datasets fetched via OpenML:energy\(heating loadY1Y\_\{1\},768×8768\\times 8\),concrete\(1030×81030\\times 8\),kin8nm\(8192×88192\\times 8\), andcalifornia\(20640×820640\\times 8\)\. Per seed, we randomly partition intontr=200n\_\{tr\}=200training andnte=100n\_\{te\}=100test points and standardizeXXandyyto zero mean and unit variance using only training\-set statistics\. Kernel bandwidthγ\\gammais chosen by the median heuristic on‖xi−xj‖2\\\|x\_\{i\}\-x\_\{j\}\\\|^\{2\}over up to500500random subsamples\. Observation noiseσn=0\.3\\sigma\_\{n\}=0\.3,1010seeds, headline sensitivity\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|pre\-committed\. We report pairedtt\-testpp\-values against uniform\.
#### B\.3 Floor ablation \(Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)\)\.
Dense synthetic GP\-prior data withntr=200n\_\{tr\}=200,55seeds,ρ∈\{0,0\.1,0\.2,0\.5,0\.7\}\\rho\\in\\\{0,0\.1,0\.2,0\.5,0\.7\\\}\. All other settings as in C\.1\. The AQKA\-GP variant is\|αiαj\|\|\\alpha\_\{i\}\\alpha\_\{j\}\|\.
#### B\.4 NLL convergence \(Figure[7](https://arxiv.org/html/2606.28833#Sx5.F7)\)\.
We measure\|ℒ\(𝐊^;𝐲\)−ℒ\(𝐊;𝐲\)\|\|\\mathcal\{L\}\(\\widehat\{\\mathbf\{K\}\};\\mathbf\{y\}\)\-\\mathcal\{L\}\(\\mathbf\{K\};\\mathbf\{y\}\)\|as a function ofBB, withℒ\\mathcal\{L\}the negative log marginal likelihood \([6](https://arxiv.org/html/2606.28833#Sx2.E6)\)\. Two settings \(planted\-sparse and dense\),ntr=200n\_\{tr\}=200,88seeds\. Kernel Frobenius error reported as‖𝐊^−𝐊‖F/‖𝐊‖F\\\|\\widehat\{\\mathbf\{K\}\}\-\\mathbf\{K\}\\\|\_\{F\}/\\\|\\mathbf\{K\}\\\|\_\{F\}\.
#### B\.5 BO surrogate quality \(Figure[8](https://arxiv.org/html/2606.28833#Sx5.F8)\)\.
For each benchmark function \(Branin in\[0,1\]2\[0,1\]^\{2\}, Hartmann\-3, Hartmann\-6\), we samplentr=120n\_\{tr\}=120uniformly random query points, evaluate the function under additive Gaussian noiseσn=0\.05\\sigma\_\{n\}=0\.05, fit a GP withγ=2\\gamma=2\(well\-tuned for the unit cube\), evaluate EI onncand=2000n\_\{\\mathrm\{cand\}\}=2000uniform random candidates, and report the function value at the argmax\-EI candidate minus the known global minimum\.88seeds\.
#### B\.6 Full BO loop \(Figure[9](https://arxiv.org/html/2606.28833#Sx5.F9)\)\.
Hartmann\-6 in\[0,1\]6\[0,1\]^\{6\}\.ninit=30n\_\{\\mathrm\{init\}\}=30uniform initial queries,T=20T=20BO iterations,Biter=5×105B\_\{\\mathrm\{iter\}\}=5\\\!\\times\\\!10^\{5\}shots per iteration\. At each iteration: refit GP withγ=2\\gamma=2,σn=0\.05\\sigma\_\{n\}=0\.05under the chosen allocator; evaluate EI onncand=1000n\_\{\\mathrm\{cand\}\}=1000uniform random candidates; pick argmax\-EI; observe true function value withσn=0\.05\\sigma\_\{n\}=0\.05noise; append to dataset\. Report simple regretmintyt−f∗\\min\_\{t\}y\_\{t\}\-f^\{\*\}\.88seeds\.
#### B\.7 Online streaming \(Figure[10](https://arxiv.org/html/2606.28833#Sx5.F10)\)\.
T=60T=60streaming steps after a3030\-point warm\-up\. At each step the GP is refit from scratch \(so the training set grows by one each step\),Bstep=2×105B\_\{\\mathrm\{step\}\}=2\\\!\\times\\\!10^\{5\},d=4d=4,γ=0\.5\\gamma=0\.5,σn=0\.3\\sigma\_\{n\}=0\.3\. Running test MSE is the cumulative squared residual averaged over steps1:t1\{:\}t\.66seeds\.
#### B\.8 Sparse VFE \(Figure[11](https://arxiv.org/html/2606.28833#Sx5.F11)\)\.
ntr=200n\_\{tr\}=200,m=30m=30inducing points selected as a uniform random subset of training inputs; dense GP\-prior data;σn=0\.3\\sigma\_\{n\}=0\.3\. The kernel budget coversm\(m\+1\)/2\+nm=6465m\(m\+1\)/2\+nm=6465entries\. Sensitivities follow Eq\. \([23](https://arxiv.org/html/2606.28833#Sx3.E23)\); uniform floor over both blocks \(treated as one entry list\)\.66seeds\.
#### B\.9NN\-scaling \(Figure[12](https://arxiv.org/html/2606.28833#Sx5.F12)\)\.
ntr∈\{50,100,200,400,600\}n\_\{tr\}\\in\\\{50,100,200,400,600\\\},B=50npairsB=50n\_\{\\mathrm\{pairs\}\}\(soBBscales withn2n^\{2\}\)\.55seeds pernn\. Dense synthetic GP\-prior data\.
#### B\.10–13 Extension experiments \(Bayesian quadrature, heteroscedastic GP, hyperparameter learning, multi\-output GP\)\.
These are detailed where the corresponding results are reported in Section “Extensions”\. For each:ntr=200n\_\{tr\}=200\(or120120for BQ\),66seeds, dense GP\-prior data unless otherwise noted\.
#### B\.14 Quantum\-kernel studies \(Figure[13](https://arxiv.org/html/2606.28833#Sx5.F13)\)\.
We replace the RBF kernel of B\.1–B\.13 with a genuine quantum fidelity kernelK\(xi,xj\)=\|⟨ϕ\(xi\)\|ϕ\(xj\)⟩\|2K\(x\_\{i\},x\_\{j\}\)=\|\\langle\\phi\(x\_\{i\}\)\|\\phi\(x\_\{j\}\)\\rangle\|^\{2\}evaluated by statevector simulation\. To approximate device\-induced kernel decay we apply a5%5\\%depolarizing channel \(Kij↦0\.95Kij\+0\.025K\_\{ij\}\\mapsto 0\.95\\,K\_\{ij\}\+0\.025\) before layering Bernoulli shot noise\. The figure reports three sub\-studies, all with55seeds and reps=2=2\.
*Study 1 \(UCI through ZZ\)\.*For each dataset in\{\\\{energy,concrete,kin8nm\}\\\}we take the firstqqfeatures \(standardized then squashed to\[0,π\]\[0,\\pi\]viatanh\\tanh\), feed them into aqq\-qubit ZZ feature map withfullentanglement andq∈\{4,6,8\}q\\in\\\{4,6,8\\\}, and run AQKA\-GP atntr=60n\_\{tr\}=60,nte=40n\_\{te\}=40,B=64npairsB=64\\,n\_\{\\mathrm\{pairs\}\}\.
*Study 2 \(feature\-map sweep\)\.*On planted\-sparse synthetic data atq=4q=4,ntr=50n\_\{tr\}=50,nte=30n\_\{te\}=30we compare four feature maps from Qiskit:ZZFeatureMap\(entanglement="full"\),ZZFeatureMap\(entanglement="linear"\),PauliFeatureMap\(paulis=\["Z","ZZ"\]\),PauliFeatureMap\(paulis=\["Y","YY","ZZ"\]\)\. Budget multipliers5050and200200\.
*Study 3 \(scale\)\.*On planted\-sparse synthetic data atq=4q=4ZZ\-full, we varyntr∈\{40,80,120\}n\_\{tr\}\\in\\\{40,80,120\\\}atB=100npairsB=100\\,n\_\{\\mathrm\{pairs\}\}\. The number of anchors scales asmax\(5,ntr/10\)\\max\(5,n\_\{tr\}/10\)\.
#### Hyperparameters\.
The choices above \(ρ=0\.5\\rho=0\.5,ρw=0\.1\\rho\_\{w\}=0\.1\) were selected from the floor\-ablation grid \(Figure[6](https://arxiv.org/html/2606.28833#Sx5.F6)\) and held fixed across all other experiments\. No per\-experiment tuning was performed\.
#### Software\.
Experiments are implemented in NumPy \(≥2\.0\\geq 2\.0\)\. UCI datasets are loaded throughsklearn\.datasets\.fetch\_openml, and the ZZ\-feature\-map kernels are evaluated via QiskitZZFeatureMapstatevector simulation\. Compute: server runs on a4040\-core CPU node\. Seeds are deterministic via NumPy’sdefault\_rng\(seed\)\.Similar Articles
Sequential sparse Gaussian process quantile regression
This paper presents a sparse Gaussian process framework for quantile regression that uses a Laplace approximation for posterior inference and variance-based mechanisms for adaptive inducing-input placement and data acquisition.
Supervised Latent Restructuring for Small-Data Quantum Learning in Plant Phenomics
This paper proposes a hybrid quantum-classical workflow for plant phenomics classification under small-data regimes, using supervised latent restructuring (PCA + LDA) to improve geometric separability before quantum kernel alignment. Experiments show improved separability but highlight compression trade-offs and the difficulty of achieving strong quantum performance.
Active Timepoint Selection for Learning Measure-Valued Trajectories
This paper introduces a framework for active timepoint selection to infer probability paths from sparse snapshots, using linearized optimal transport to map distributions into a tangent space for Gaussian Process modeling, thereby enabling uncertainty-aware acquisition policies.
Boundary Variance Inflation Causes Acquisition Bias in Gaussian Processes
This paper identifies the geometric mechanism behind boundary-induced acquisition bias in Gaussian processes on bounded domains, showing how kernel truncation inflates posterior variance and distorts acquisition functions independently of the objective function. The authors introduce a function-free diagnostic to quantify this bias across different acquisition classes.
QUIVER: Quantum-Informed Views for Enhanced Representations in Large ML Models
This paper introduces Quiver, a paradigm that enriches classical machine learning models with quantum-inspired features derived from the quantum Fisher information matrix, demonstrating improvements on molecule property prediction and jet flavor classification benchmarks.