Quantum Spectral Anomaly Detection Yewei Yuan,1 Michele Minervini,2 Mark M. Wilde,2 and Nana Liu1, 3, 4 1
arXiv:2607.05307v1 [quant-ph] 6 Jul 2026
Global College, Shanghai Jiao Tong University, Shanghai 200240, China 2 School of Electrical and Computer Engineering, Cornell University, Ithaca, New York 14850, United States 3 Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China 4 School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai 200240, China (Dated: July 7, 2026) A core task in quantum anomaly detection is to compute an anomaly score that quantifies how strongly a test quantum state deviates from a given quantum dataset assumed to be normal. Classically, principal component analysis (PCA) for centered data computes the anomaly score by evaluating the test sample relative to the subspace spanned by the selected leading eigenvectors. However, for quantum data that lack a standard centering, explicitly recovering principal eigenvectors, constructing full Gram matrices, or loading quantum-random-access-memory-style data can be more costly than estimating the anomaly score itself. To avoid these costs, we propose Quantum Spectral Anomaly Detection (QSPADE), which computes PCA-like anomaly scores directly from the spectrum of the average state of the normal dataset. By replacing hard PCA rank selection with a smooth, temperature-controlled spectral threshold, QSPADE makes near-threshold spectral components contribute partially to the anomaly score. This makes the score vary continuously rather than jump when a borderline component is included or excluded, and makes it less sensitive to noise or arbitrary hard cutoffs near the threshold. In the zero-temperature limit, QSPADE recovers the hard-projector PCA score. The proposed measurement-based quantum detector can be calibrated with a sample complexity independent of the data dimension. Numerical simulations show that QSPADE behaves like kernel-PCA on encoded classical data and detects changes across a transverse-field Ising transition without predefined order parameters. Consequently, QSPADE gives an efficient framework for both quantum-kernel anomaly detection on encoded classical data and the monitoring of quantum-native systems where diagnostic observables are unknown.
I.
INTRODUCTION
Anomaly detection is a one-class learning problem: given training examples that are assumed to be normal, we want to assign a new test object an anomaly score, with larger values indicating stronger deviation from the normal data. Principal component analysis (PCA) [1] gives a spectral way to construct such scores. After centering the normal data, the sample covariance matrix identifies directions of variation: large eigenvalues correspond to directions strongly represented in the normal data, while small eigenvalues correspond to weakly represented directions. The large-eigenvalue subspace is the principal subspace. Relative to this subspace, PCA produces two anomaly scores. First, it can have a large residual weight outside this subspace, as measured by the residual statistic Q. Second, it can lie inside the subspace but have unusually large weight along a low-variance direction, as measured by Hotelling’s T 2 statistic [2, 3]. For either method, the desired output is therefore a single scalar anomaly score. This distinction is especially important for quantum data. When the data are naturally quantum, such as states produced by a quantum device, or when classical data are already represented as quantum feature states, the task is to estimate this scalar score without requiring full spectral decomposition information. However, the standard PCA workflow is not naturally matched to this goal. Classically, one constructs and diagonalizes a covariance or Gram ma-
trix for centered data, chooses a target rank according to a retained-variance criterion, stores the corresponding eigenvectors, and projects each test point onto them. For quantum data, such tasks may solve a harder problem than anomaly detection itself. These extra costs are also visible in existing quantum anomaly detection methods [4]. Quantum kernel PCA and one-class SVM [5] use quantum kernels to obtain nonlinear feature representations, but rely on a pairwise-kernel representation and quantum-randomaccess-memory-style access assumption. Quantum PCA and quantum linear-algebra algorithms provide direct spectral access to density operators or block-encoded matrices [6, 7]. However, standard qPCA usually aims at extracting or projecting onto principal spectral information, rather than directly calibrating anomaly scores. For quantum native data, reproducing standard PCA also requires treating centering explicitly, since no classical feature vectors or centered Gram matrix are directly available. Variational one-class models and quantum autoencoders provide another route [8–10], but they replace reconstruction by trainable circuits, which introduce optimization loops and trainability concerns such as barren plateaus [11]. Recent works show that full eigenvector recovery is unnecessary when PCA score is the only target output [12, 13]. They reduce the complexity and sensitivity caused by eigenvalue resolution in standard PCA. Following this development, we seek a simple end-to-end quantum anomaly detection method avoiding the extra reconstruction, matrix-access, rank-selection, and opti-
2 mization cost. In this paper, we propose Quantum Spectral Anomaly Detection (QSPADE), a measurement-based approach for computing PCA-like anomaly scores directly from the spectrum of the average state of the normal dataset. Given N ∈ N normal states ρi , the average state is a mixed state X C := ρ̄ = N −1 ρi . i
QSPADE then defines the regularized measurementbased spectral detector C − µI f := Mµ,T f , T where f is a monotone smooth response function taking f values in [0, 1] (so that Mµ,T is a measurement operator), µ is a scalar threshold corresponding to the prescribed acceptance level, and the temperature T > 0 controls how sharply the response changes near the threshold. Specifically, f responds to the comparison results between the eigenvalues of C and the threshold µ. In the limit T → 0, this response becomes a step function, recovering the hard-projector PCA score. Thus, higheigenvalue spectral components of C contribute little to the residual score, low-eigenvalue components contribute strongly, and near-threshold components contribute fracf tionally. For a test state σ, the probability 1−Tr(Mµ,T σ) corresponds to the final residual score. The threshold µ is calibrated by matching the average acceptance probaf bility Tr(Mµ,T C) to a target retained mass, rather than choosing how many eigenvectors to keep. The calibration requires O(η −2 log δ −1 ) sample complexity, independent of the data dimension, where η > 0 and δ ∈ (0, 1) are the target precision and failure probability, respectively. This work contributes an efficient quantum anomaly detection method that focuses on directly evaluating target anomaly scores from multiple copies of given quantum data, bypassing redundant costs. We show how classical Q and Hotelling’s T 2 scores translate to a quantum system. We provide corresponding quantum implementations using the measurement-based spectral detector, including both qumode and qubit circuits for different anomaly scores. The qubit approach is natural for singlethreshold residual score calculation, while the qumode approach excels at resolving how different spectra contribute to the anomaly score. We also identify the distinct calibration resources required by these methods. Finally, we validate the framework numerically: QSPADE successfully recovers the nonlinear decision geometries, and it robustly tracks quantum phase transitions in the transverse-field Ising model without requiring any prior specification of order parameters. The results show that QSPADE is effective for both quantum-embedded classical data and quantum-native data. The rest of this paper is organized as follows. Section II reviews the classical PCA monitoring target. Section III
defines nominal spectral support for quantum data and states the hard spectral detector. Section IV introduces the regularized spectral detector. Section V defines the QSPADE Q and T 2 statistics and gives their sharp-limit correspondence with the classical quantities. Section VI gives the qumode and qubit implementation. Section VII reports numerical experiment results.
II.
CLASSICAL PCA MONITORING BACKGROUND
For d ∈ N, let x1 , . . . , xN ∈ Rd be training vectors assumed to be normal. Classical PCA monitoring first centers the data, N
x̄ :=
1 X xi , N i=1
yi := xi − x̄,
(II.1)
and then forms the sample covariance N
Ccl :=
1 X yi yi⊤ . N i=1
(II.2)
This centering step is part of the model, not merely preprocessing: the resulting statistics detect deviations around the mean vector. A test input x is therefore evaluated through its centered form z = x − x̄. Write the spectral decomposition of the sample covariance as Ccl =
d X
λj uj u⊤ j ,
λ1 ≥ λ2 ≥ · · · ≥ λd ≥ 0, (II.3)
j=1
and let PK :=
K X
uj u⊤ j
(II.4)
j=1
be the rank-K principal projector. The rank K is often chosen by a given retained variance target α ∈ (0, 1): ( ) PK j=1 λj Kα := min K : Pd ≥α . (II.5) j=1 λj For a centered test vector z, the classical residual statistic called Q is defined as X 2 2 QK (z) := ∥(I − PK )z∥ = |⟨uj , z⟩| . (II.6) j>K
It is large when the input has energy outside the retained principal subspace, which is commonly referred to as the nominal support. Hotelling’s statistic instead measures how the test vector is distributed inside the retained principal subspace. It uses the squared components of z along the retained
3 eigenvectors u1 , . . . , uK , and weights each component by the inverse of the corresponding covariance eigenvalue: 2 TK,γ (z) :=
K X |⟨uj , z⟩|2 j=1
λj + γ
,
γ ≥ 0.
(II.7)
The ridge parameter γ is included to make T 2 stable when retained eigenvalues are small or poorly estimated. The two statistics have different roles. In one-class support monitoring, Q is usually the default score because it directly tests whether a sample leaves the principal subspace and does not require inverse spectral weights. Hotelling’s T 2 is useful when anomalies are expected to remain within the nominal support but have abnormal leverage, for example, a large score in a low-variance retained eigenvector. It requires more structure: a reliable choice of which eigenvectors are retained, stable variance denominators, and a calibration rule for the resulting normalized score. QSPADE maintains this distinction. The quantum residual score is defined in Section V A; the quantum Hotelling-type score is defined in Section V B. III.
choices. Therefore, no default centering operation preserves the physical sample space and the Born-probability interpretation in (III.3). QSPADE consequently uses the uncentered average state in (III.1). For classical data, the recommended convention is different. If the goal is to reproduce centered PCA monitoring, one centers in the classical domain first, yi = xi − x̄ and z = x − x̄, and then embeds the centered vectors into quantum states. After this embedding, the resulting objects are treated as quantum data by QSPADE. Let a spectral decomposition of C be as follows: C=
d X
λj |uj ⟩⟨uj |,
The high-eigenvalue eigenspaces of C represent the spectral components most strongly supported by the normal data. The following proposition formalizes this statement. Proposition III.2 (Hard spectral detector). For 1 ≤ K ≤ d,
NOMINAL SUPPORT FOR QUANTUM DATA
We now translate the PCA anomaly detection idea to the quantum setting. In classical PCA, the nominal support is the high-variance spectral region of the covariance matrix. For quantum data, the analogous object is the high-eigenvalue spectral region of the average state of the normal dataset {ρi }N i=1 : 1 X ρi . N i=1
(III.1)
Let H be a d-dimensional Hilbert space, and define the effect set M(H) := {M ∈ Herm(H) : 0 ⪯ M ⪯ I}.
(III.2)
For every measurement effect M , N
Tr(M ρ̄) =
1 X Tr(M ρi ). N i=1
max
M ∈M(H) rank(M )≤K
Tr(M C) =
K X
λj .
(III.5)
j=1
The maximum is attained by any top-K spectral projector PK =
K X |uj ⟩⟨uj |.
(III.6)
j=1
If λK > λK+1 , the top-K projector is unique.
N
C := ρ̄ =
λ1 ≥ λ2 ≥ · · · ≥ λd ≥ 0. (III.4)
j=1
Proof. See Appendix A. Pd Since C is a density operator, j=1 λj = 1. Thus, the quantum analogue of a retained-variance target is a retained-mass target α ∈ (0, 1], for which the hard rank is chosen as K X Kα := min K : λj ≥ α . (III.7) j=1
(III.3)
The operational meaning of (III.3) is simple: an effect M accepts normal data with probability Tr(M C). Thus, learning a nominal support can be phrased as choosing a measurement effect that accepts the average normal state as much as possible under a rank constraint. Remark III.1. Quantum-native states are already normalized objects: pure states are on the unit sphere up to phase, and density operators have trace one and are positive semidefinite. Subtracting a mean state generally produces a matrix that is not a state, and subtracting vector amplitudes is not invariant under global phase
This hard choice accepts the first Kα eigenspaces completely and rejects the remaining eigenspaces completely. Proposition III.2 identifies the standard principal subspace with the optimal hard detector. The next section replaces this discontinuous hard acceptor by a smooth, optimization-derived measurement effect.
IV.
SOFT SPECTRAL DETECTOR
The hard projector PK makes a binary decision on the eigenspaces of C: each eigenspace is either fully accepted or fully rejected. It is sensitive to eigenvalue crossings
4 and becomes unstable when the spectral gap near the cutoff is small. A well-behaved detector should remain a valid effect, increase monotonically with the eigenvalues of C, and admit a principled optimization formulation. We derive such a detector by introducing a regularized optimization problem. We first replace the step response in the binary decision with a smooth monotone response function. Definition IV.1 (Optimization-admissible response). A scalar response f : R → (0, 1) is optimization-admissible if:
Proposition IV.3 (Sharp-resolution limit). For fixed µ, f lim Mµ,T = Π(C > µI) + f (0)Π(C = µI).
T →0+
(IV.5)
f In particular, if λK > µ > λK+1 , then Mµ,T → PK .
Proof. See Appendix C. This parameterized family assigns smooth acceptance weights to the eigencomponents of C at finite resolution, while recovering the hard spectral cutoff in the sharp limit.
i. f ∈ C 1 (R) and f ′ (z) > 0 for every z ∈ R; V.
ii. limz→−∞ f (z) = 0 and limz→+∞ f (z) = 1; iii. the quantile function f −1 belongs to L1 (0, 1); R1 iv. 0 f −1 (m) dm = 0. The associated generator is Z m ϕf (m) := f −1 (u) du,
m ∈ (0, 1),
(IV.1)
0
with continuous extension to [0, 1]. The last condition only fixes the location convention for the response. Logistic, Gaussian, and other smooth thresholds can be placed in this form after centering. Fix T > 0 and µ ∈ R. Define the regularized spectral detector as the optimizer of min {µ Tr(M ) − Tr(M C) + T Tr ϕf (M )} .
0⪯M ⪯I
(IV.2)
Here µ plays the role of a spectral threshold, or equivalently, a Lagrange parameter controlling the detector size, while the temperature parameter T controls the smoothness of the transition. The term − Tr(M C) rewards normal acceptance, µ Tr(M ) penalizes detector size, and T Tr ϕf (M ) smooths the occupation pattern. Theorem IV.2 (Soft Spectral Detector). Let C ⪰ 0, T > 0, and let f be optimization-admissible. Problem (IV.2) has the unique solution C − µI f . (IV.3) Mµ,T = f T
We now turn the soft spectral detector into anomaly scores. We define the two QSPADE scores and show how the scalar threshold µ is empirically calibrated to match a desired retained-mass target (the continuous analogue), and how test states are scored based on this threshold. The residual Q statistic only needs a retained-mass threshold, calibrated directly with the average state C. A Hotelling-type statistic additionally needs spectral windows ranks, calibrated with the maximally mixed probe τmm = I/d. A.
Residual statistic Q and retained-mass threshold
For fixed f and T , define the retained-mass function d X λj − µ f f RC,T (µ) := Tr Mµ,T C = λj f . T j=1
d X λj − µ f |uj ⟩⟨uj |. T j=1
(IV.4)
Proof. See Appendix B. Let Π(C > µI) denote the projection onto the eigenspace of C − µI with strictly positive eigenvalues, and let Π(C = µI) denote the projection onto the zero eigenspace of C − µI.
(V.1)
Since C = ρ̄ is a density operator and f is strictly inf creasing with limits 0 and 1, RC,T (µ) decreases continuously from 1 to 0 as µ increases from −∞ to +∞. Thus, for each retained-mass target α ∈ (0, 1), we choose the unique threshold µα by f RC,T (µα ) = α.
(V.2)
The calibrated detector is defined as Mαf := Mµfα ,T .
Equivalently, f Mµ,T =
QSPADE STATISTICS AND CALIBRATION TARGETS
(V.3)
For a normalized test state σ, define the acceptance score and residual anomaly score under the retained detector Sαf (σ) := Tr(Mαf σ),
Qfα (σ) := 1 − Sαf (σ).
(V.4)
The anomaly decision cutoff for Qfα is fixed in advance. Alternatively, to guarantee a specific false-alarm rate, the decision threshold can be empirically bounded using held-out scores, consistent with standard classical process monitoring.
5 B.
Hotelling’s T 2 from soft spectral windows
A Hotelling-type score requires more information than the residual score: it must resolve how a test state is distributed within the part of the spectrum accepted as normal. For this purpose, we introduce a finite sequence of thresholds and take differences between consecutive soft detectors. These differences define soft spectral windows. Each window collects contributions from eigencomponents lying mainly between two neighboring thresholds, with smooth boundary weights determined by the resolution T . The number of windows, L, is chosen so that the retained windows cover the target retained mass. Choose ordered thresholds +∞ = µ0 > µ1 > · · · > µL ,
(V.5)
and set M0 := 0,
Mℓ := Mµfℓ ,T
(ℓ = 1, . . . , L).
(V.7)
They are positive semidefinite because lowering the threshold increases the spectral response. If the final threshold is −∞, then the windows sum to I. For component-wise Hotelling behavior, these retained window boundaries should be rank-calibrated: Tr(Mµfℓ ,T ) = ℓ,
ℓ = 1, . . . , L
(L ≪ d).
vℓ ≥ α.
(V.10)
ℓ=1
For a test state σ, the test window score is pσ (ℓ) := Tr(Dℓ σ).
(V.11)
For a ridge γ ≥ 0, define the soft leverage score 2,f TL,γ (σ) :=
L X pσ (ℓ) ℓ=1
vℓ + γ
For a pure test state σ = |ψ⟩⟨ψ| with |ψ⟩ = the residual statistic converges to X Qfα (σ) → |aj |2 ,
(V.13) P
j aj |uj ⟩,
(V.14)
j>L
and the soft leverage score converges to 2,f TL,γ (σ) →
L X |aj |2 . λ +γ j=1 j
(V.15)
Thus, when C is obtained from centered classical features, QSPADE recovers the classical PCA residual and Hotelling forms. Proof. See Appendix D.
(V.9)
Since Dℓ is calibrated to represent one spectral component, vℓ plays the role of the corresponding soft eigenvalue. The total number of retained windows L can be chosen as the minimal integer L satisfying L X
Dj → |uj ⟩⟨uj |.
(V.8)
This condition places the threshold µℓ at the ℓ-th PCA rank position in the soft spectrum. Hence the window Dℓ is the soft analogue of the ℓ-th principal component. For each window, define the retained window mass vℓ := Tr(Dℓ C),
Proposition V.2 (Sharp-limit correspondence with classical Q and T 2 ). Suppose that C has nondegenerate positive eigenvalues in the retained range. We consider the limit T → 0+ , with thresholds placed between consecutive eigenvalues. Then, in this limit,
(V.6)
The window effects are defined as Dℓ := Mℓ − Mℓ−1 .
Remark V.1 (Boundary leverage in the soft residual). In standard hard PCA, the Q-statistic is entirely insensitive to variance differences among the retained components, making T 2 strictly necessary for leverage monitoring. In QSPADE, however, the soft residual Qfα applies a graded penalty 1 − f ((λj − µ)/T ) to the boundary modes. Because this fractional penalty naturally tracks the spectral decay in the transition region, the soft Q statistic absorbs some of the variance-weighting functions of a traditional leverage score near the cutoff. Consequently, for many cases, the Q statistic alone is sufficiently informative, relegating T 2 to cases where anomalous leverage deep within the dominant support is explicitly suspected.
.
(V.12)
The key calibration distinction is that Q requires only one retained-mass threshold µα , whereas T 2 requires component-resolved windows in addition, corresponding to thresholds µℓ for L windows.
VI.
QUANTUM MEASUREMENT IMPLEMENTATIONS
f The soft spectral detector Mµ,T can be accessed without explicit eigenvector reconstruction. Table I summarizes the two calibration objects. The qumode implementation records a continuous spectral sample, making it naturally suited for global retained-mass sweeps and window-resolved T 2 leverage. Conversely, the qubit implementation estimates a single target threshold via Hadamard tests for the residual Q detector. Throughout this section, we use C = ρ̄ as a Hamiltonian generator. This access can be supplied either by a Hamiltonian-simulation oracle or block-encoding for C, or, in the sample-access setting, by density matrix exponentiation [6] (see also [14, 15]) using copies of C, which are prepared by sampling normal training states uniformly. Thus the qumode interaction e−ip̂⊗C/T2 in (VI.3) and the qubit evolutions eitω C/T in (VI.19) are both implementations of the same underlying simulation primitive for C.
6 Statistic Default?
Calibrated Quantity
Required Probe
Recommended Method
Sample Complexity
Q
yes
retained-mass threshold µα
C
T2
optional
thresholds µℓ for L windows C, τmm = I/d
qumode (sweep) / qubit (single-target) qumode
O(η −2 log δ −1 ) ((VI.10), (VI.24)) −2 e O(Ldη log δ −1 ) L ≪ d ((VI.14))
for
TABLE I. Summary of QSPADE anomaly scores, required calibration, and sample complexities. The residual score Q is the default one-class support test, whereas the Hotelling-type score T 2 requires additional windows. Here C is the average state, τmm = I/d is the maximally mixed probe (with system dimension d), µℓ are thresholds for the L windows, η is the target additive accuracy, and δ is the failure probability.
A.
Qumode spectral-threshold measurement
We first present a quantum implementation of the f soft detector Mµ,T via an ancillary qumode interacting with C, followed by a classical threshold on the position outcome. Since f ′ > 0 and f (+∞) − f (−∞) = 1, the derivative f ′ is a probability density function. For T1 > 0, define the reflected position density 1 ′ q gf,T1 (q) := f − . (VI.1) T1 T1 Prepare a control qumode in the position wavepacket Z q gf,T1 (q)|q⟩ dq. (VI.2) |ψf,T1 ⟩ := R
T2 > 0,
(VI.3)
measure the position quadrature q̂, and accept if the outcome q > β, where β ∈ R is chosen based on µ and T . Proposition VI.1 (General spectral-threshold measurement). The above procedure implements the binary POVM effect C − βT2 I f := Eβ,T f . (VI.4) 1 ,T2 T1 T2 Equivalently, with T := T1 T2 , f f one has Eβ,T = Mµ,T . 1 ,T2
The shift µI is not simulated in the circuit. Once UC and the control width are fixed, changing µ only changes the classical threshold β applied to the measured position.
B.
Qumode calibration for Q and T 2
We now explain how the qumode-based detector introduced in Section VI A is calibrated for Q and T 2 . For the residual statistic Q, prepare the average state C. With the circuit in Proposition VI.1, define GC (β) := Pr(q > β | C).
(VI.6)
Then
The response f is chosen from the optimizationadmissible class in Definition IV.1. This choice is not only a regularization choice: it also fixes the physical control distribution through f ′ . Two useful examples are the Fermi–Dirac (logistic) response, which connects to entropy-regularized soft spectral measurements [13, 16], and the Gaussian (probit) response, whose control qumode is a standard Gaussian state [17] and connects to the Gaussian error function activation observable from [18]. Their explicit densities and details are given in Appendix E. Apply UC := e−ip̂⊗C/T2 ,
Proof. See Appendix F.
µ := βT2 ,
(VI.5)
f f GC (β) = Tr(Eβ,T C) = Tr(MβT C). 1 ,T2 2 ,T1 T2
(VI.7)
Thus a retained-mass target α is calibrated by GC (βα ) = α,
µα = βα T2 .
(VI.8)
Because the threshold β is applied in post-processing, the same recorded position samples can be reused to calibrate any number of targets α. This global sweep procedure is summarized in Algorithm 1. Algorithm 1 Qumode retained-mass calibration for residual Q Inputs: Average state C, response f , control scales T1 , T2 , successful calibration shots S, target set A ⊂ (0, 1). 1: Prepare |ψf,T1 ⟩ with density (VI.1). 2: for s = 1, . . . , S do 3: Prepare C, apply UC = e−ip̂⊗C/T2 , and record the position outcome qs . 4: end for 5: Form the empirical tail function S X b S (β) = 1 G 1{qs > β}. S s=1
6: for each α ∈ A do b S (βbα ) ≈ α and set µ 7: Choose βbα satisfying G bα = βbα T2 . 8: end for 9: return thresholds {b µα : α ∈ A}.
7 Uniformly over β, the Dvoretzky–Kiefer–Wolfowitz inequality [19] gives ! b S (β) − GC (β)| ≤ η ≥ 1 − δ Pr sup |G (VI.9) β
if S≥
2 1 log . 2 2η δ
(VI.10)
Thus, the retained-mass curve is learned to additive accuracy η with O(η −2 log δ −1 ) samples, independent of d. For T 2 , the same qumode measurement circuit is used to construct component-resolved windows Dℓ . The maximally mixed probe τmm := dI gives f Pr(q > β | τmm ) = Tr(Eβ,T τmm ) 1 ,T2 1 f = Tr(MβT ). 2 ,T1 T2 d
To calibrate the thresholds µℓ for the retained windows to an additive accuracy η for the rank level, the required probability accuracy in the τmm probe is η/d. Because we restrict Hotelling leverage monitoring to L ≪ d principal windows, the relevant binomial success probability for these top modes under the maximally mixed probe is small, p ≈ L/d. Applying Bernstein’s inequality [20] leverages the reduced variance p(1 − p) ≈ L/d, yielding a significantly sharper sample complexity: Ld d 1 + ST 2 = O log . (VI.14) η2 η δ Thus, by exploiting the finite-window truncation L ≪ d, e 2 ) worst-case scaling, reQSPADE breaks the naive O(d ducing the component-resolution rank calibration overe head to O(Ld).
(VI.11) C.
Single-target qubit calibration for residual Q
Thus, the rank-calibrated thresholds are obtained from ℓ Pr(q > βℓ | τmm ) = . d
(VI.12)
For windows Wℓ = (βℓ , βℓ−1 ], the Hotelling denominator is then the retained window mass vℓ = Pr(q ∈ Wℓ | C) = Tr(Dℓ C).
(VI.13)
Algorithm 2 below gives the full Hotelling-type scoring procedure: the maximally mixed probe fixes componentresolved windows, and the average state C supplies the denominators. Algorithm 2 Qumode Hotelling-type score Inputs: Same fixed qumode circuit as Algorithm 1; average state C; test state σ; maximally mixed state τmm = I/d; retained-mass target α; candidate window number Lmax ; ridge γ ≥ 0. bℓ satisfying 1: Use τmm to choose thresholds β c > βbℓ | τmm ) ≈ ℓ , Pr(q d
ℓ = 1, . . . , Lmax .
bℓ , βbℓ−1 ], with βb0 = +∞. 2: Define windows Wℓ = (β 3: Estimate retained window masses c ∈ Wℓ | C), vbℓ = Pr(q 4: Choose
ℓ = 1, . . . , Lmax .
L′ X ′ b = min L : L vbℓ ≥ α . ℓ=1
5: Estimate test-window masses
c ∈ Wℓ | σ), pbσ (ℓ) = Pr(q
b ℓ = 1, . . . , L.
When the retained-mass target α is fixed in advance, it is unnecessary to collect a full qumode distribution. A qubit route can calibrate the single scalar threshold µα by directly evaluating the nonlinear spectral filter Tr[f ((C − µI)/T )ρ] using a hybrid Hadamard-test estimator adapted from Algorithm 2 of [18]. Since a standard Hadamard test [21] evaluates expectation values of the form Tr(eiHt ρ), we bridge the gap to the nonlinear response f by representing f (x) as a bounded superposition of phase functions eitx on the relevant spectral interval (Definition VI.2). By substituting the operator argument into this phase expansion, evaluating the target filter decomposes into estimating a sequence of unitary terms eit(C−µI)/T . Because the scalar shift µI commutes with C, the evolution eit(C−µI)/T trivially factors as eit(C−µI)/T = e−iµt/T eitC/T .
Thus the data register only needs controlled simulations of C; changing the threshold µ only changes a scalar phase on the control qubit. Definition VI.2 (Hadamard-estimable response). Fix a compact interval J ⊂ R. A response f : R → (0, 1) is called (Af , εf )-Hadamard-estimable on J if there exist a real constant bf , a constant Af > 0, a probability space Ωf , real frequencies tω ∈ R, and coefficients χω ∈ C with |χω | ≤ 1 such that fe(x) := bf + Af Eω∼Ωf ℜ χω eitω x (VI.16) satisfies sup |f (x) − fe(x)| ≤ εf .
6: Output 2,f TbL,γ b (σ) =
b L X ℓ=1
pbσ (ℓ) . vbℓ + γ
(VI.15)
(VI.17)
x∈J
The interval J is chosen so that (λ − µ)/T ∈ J for all λ ∈ spec(C) and all thresholds µ considered
8 during calibration. Since C = ρ̄ is a density operator, spec(C) ⊆ [0, 1]; hence for a calibration interval [µmin , µmax ] one may take µmax 1 − µmin , , (VI.18) J := − T T or a smaller interval if tighter spectral bounds are known. By casting f into this phase-superposition form, the evaluation of the nonlinear filter translates directly into random sampling over controlled-unitary evolutions. Representations of the form in (VI.16) arise from Fourier approximations of f or other truncated phase expansions. They naturally accommodate quantum signal processing (QSP) polynomials [7, 22] or linear combinations of unitaries (LCU) [23]. The quantity Af captures the overall sampling prefactor or block-encoding norm of the chosen representation. To construct the single-shot estimator Yµ (ρ) for the target filter Tr[f ((C − µI)/T )ρ], the process is as follows: 1. Sample ω ∼ Ωf and write χω = |χω |eiφω (choosing an arbitrary phase φω if χω = 0).
Proof. See Appendix G. For fixed µ and ρ, averaging m independent copies of Yµ (ρ) gives, by Hoeffding’s inequality [24], additive error at most η + εf with probability at least 1 − δ using ! A2f 1 log (VI.24) m=O η2 δ Hadamard-test samples, multiplied by the average or worst-case cost of simulating eitω C/T over the sampled frequencies. For residual calibration, the probe state is the average state C = ρ̄. Operationally, this means sampling a normal training state uniformly. With ρ = C, Proposition VI.3 gives noisy observations of C − µI f RC,T (µ) = Tr f C . (VI.25) T The threshold µα is then found by Robbins–Monro stochastic approximation [25] applied to the monotone equation f RC,T (µ) = α.
2. Prepare the data register in state ρ and a control qubit in |0⟩. Apply a Hadamard gate to the control, apply the controlled unitary Uω := eitω C/T
(VI.19)
on the data register, apply the phase θω (µ) := φω −
µtω T
(VI.20)
to the control qubit.
Algorithm 3 Single-target qubit calibration for residual Q Inputs: Hadamard-estimable response f , resolution T , target retained mass α ∈ (0, 1), interval [µmin , µmax ] containing the solution, initial µ0 ∈ [µmin , µmax ], step sizes ηn , average state C = ρ̄, mini-batch sizes mn ≥ 1, number of iterations N . 1: for n = 0, 1, . . . , N − 1 do 2: Draw independent samples
3. Measure the control qubit in the X basis to obtain a signed outcome Zω (ρ) ∈ {−1, +1}. Conditioned on the sampled frequency ω, the expected value of this measurement outcome exactly recovers the required real part: E[Zω (ρ) | ω] = ℜ[eiθω (µ) Tr(Uω ρ)].
Yµ (ρ) := bf + Af |χω |Zω (ρ).
n) (C) (C), . . . , Yµ(m Yµ(1) n n
3:
using (VI.22). Form mn X bn (µn ) = 1 R Yµ(s) (C). mn s=1 n
(VI.21)
Scaling and translating this outcome yields the final single-shot estimator:
(VI.26)
4:
Update bn (µn ) − α) . µn+1 = Π[µmin ,µmax ] µn + ηn (R
(VI.22)
Proposition VI.3 (Control-shifted Hadamard estimator). Suppose that f is (Af , εf )-Hadamard-estimable on J and that C can be simulated for the times tω /T . For every density operator ρ and every µ satisfying (spec(C) − µ) /T ⊆ J , the single-shot estimator Yµ (ρ) defined in (VI.22) satisfies C − µI E Yµ (ρ) − Tr f ρ ≤ εf . (VI.23) T Moreover, the estimator is strictly bounded, Yµ (ρ) ∈ [bf − Af , bf + Af ].
5: end for 6: return µ bα = µN .
Algorithm 3 implements this stochastic calibration. f The sign in the update is chosen because RC,T is decreasing: if the observed retained mass is larger than α, then the threshold is too low and must be increased. Under bounded conditional variance and the standard Robbins– Monro step-size conditions ∞ X n=1
ηn = ∞,
∞ X n=1
ηn2 < ∞
(VI.27)
9
D.
Deployment
Regardless of the calibration route, once the calibrated threshold µ bα is obtained, the residual detector is fixed as Mµbfα ,T . Deployment on a test state σ estimates the acceptance score and the residual anomaly score, Sαf (σ) = Tr(Mµbfα ,T σ),
Qfα (σ) = 1 − Sαf (σ). (VI.28)
The same measurement primitive is used with the calibration input replaced by the test state σ. For a binary qumode threshold measurement, additive score accuracy η requires O(η −2 log δ −1 ) shots. For the qubit Hadamard route, the corresponding sample count is O(A2f η −2 log δ −1 ), with the per-shot cost determined by the simulation of eitω C/T in the chosen representation. If Hotelling leverage T 2 is requested, the qumode implementation is the natural route: its continuous position output q can be binned into all calibrated windows from the same set of measurement records. The deployment on a test state σ is included in Algorithm 2. VII.
NUMERICAL EXPERIMENTS
We report two simulations that test the two intended uses of QSPADE. The first uses classical data presented as quantum feature states and is a consistency check: the QSPADE residual and leverage scores reproduce kernelPCA-like nonlinear spectral geometry. The second is the genuinely quantum-native setting, in which P the normal object is the mixed state C = ρ̄ = N −1 i |ψi ⟩⟨ψi | of a transverse-field Ising model. A.
Classical data: the kernel-PCA behavior
The normal data in this numerical experiment are a one-class curved point cloud (a single “moon”) in [0, 1]2 . Each point x is mapped to a normalized feature state |ϕ(x)⟩ by a fixed feature map after classical centering. Because the QSPADE scores depend on the data only through the inner products ⟨ϕ(x) | ϕ(x′ )⟩ = k(x, x′ ), QSPADE on classical data is a kernel method, and in the
Q residual
1.0
T2 in Q-accepted
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
combined max-percentile
Classical RBF kernel-PCA
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 1.0 1.0 1.0 0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 1.0 1.0 1.0 0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
1.0
0.8
normal-calibrated percentile
Classical PCA
1.0
QSAD (fourier)
satisfied, for instance, by an algebraically decaying sequence ηn = c/nν with ν ∈ (0.5, 1], the recursion converges almost surely to µα as N → ∞ when the estimator is unbiased. A deterministic representation bias εf moves the limiting point to an O(εf /κ) neighborhood of f µα whenever the local slope satisfies |(RC,T )′ (µ)| ≥ κ > 0 around the solution. In practical quantum implementations subject to finite-shot binomial noise, the raw trajectory µn may exhibit oscillations. Applying Polyak– PN Ruppert averaging [26, 27] µ̄N = N1 n=1 µn over the iterations provides a robust estimate of the true threshold.
0.6
0.4
0.2
0.0
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
FIG. 1. Normal-calibrated percentile maps over [0, 1]2 for the moon dataset; the pale scattered points are the normal training data. Rows: linear classical PCA, RBF kernel-PCA, and QSPADE with a Fourier feature map. Columns: residual Q, within-support T 2 (shown inside the Q-accepted region), and the combined max-percentile. Linear PCA gives a rigid straight band; the kernel baseline and QSPADE both produce curved acceptance regions that follow the data. The cyan contour marks the calibrated decision boundary.
sharp-resolution limit, the residual Qfα coincides with the kernel-PCA reconstruction error. This experiment makes the correspondence explicit and quantitative. We compare three monitors on identical data: (i) linear PCA in the raw two-dimensional space (centered, single retained component); (ii) standard centered radial basis function (RBF) kernel-PCA, the canonical kernel oneclass baseline, with the median-heuristic bandwidth; and (iii) QSPADE with a fixed representative feature map, a four-qubit tensor-product Fourier map. The feature maps, kernels, and all simulation parameters are specified in Appendix H 1. For each monitor, in Figure 1 we report the residual Q, the within-support Hotelling T 2 , and their combined percentile, calibrated against the empirical distribution of normal scores. Linear PCA produces rigid, straight acceptance bands that cannot follow the curvature of the data, so normal points near the extremes of the arc are flagged as anomalous. The RBF kernel-PCA and QSPADE both produce curved acceptance regions that track the manifold; the boundaries are similar though not identical, since they realize different kernels. In Table II, we quantify detection accuracy by the area under the receiver-operating-characteristic curve (ROCAUC), a threshold-independent score equal to 1 for perfect normal/anomaly separation and 1/2 for chance. Across a suite of geometries that spans an elongated
10 1.0
0.93 0.87 0.82
0.93 0.88 0.88
TABLE II. ROC-AUC for one-class detection on three classical geometries: an elongated Gaussian blob, the moon, and a ring. Linear PCA degrades on curved or non-convex dataset, while the RBF kernel-PCA baseline and QSPADE with a Fourier feature map remain comparable. QSPADE supplies a PCA-style one-class anomaly score for quantum-kernel methods: whenever a quantum feature map is advantageous for classical data, the same access model can be used for calibrated anomaly detection.
Gaussian blob (the linear-correlation case), the moon, and a ring, linear PCA degrades from ROC-AUC ≈ 0.89 on the blob to ≈ 0.57 on the ring, whereas both the kernel baseline and QSPADE remain in the 0.82–0.93 range (see also Figure 7 in Appendix H). The agreement is not specific to one kernel: a comparison across further classical kernels (Laplacian and polynomial) and further feature maps shows the same data-dependent behavior on both sides (Figure 6 in Appendix H).
B.
Quantum-native data: support across a phase transition
We now apply QSPADE to genuinely quantum-native data. This experiment tests whether a soft retained-mass support detector can monitor a quantum phase change when the normal state has a graded spectrum and no clear PCA rank. We use the transverse-field Ising Hamiltonian HTFIM := −J
n−1 X
Zi Zi+1 − h
n X
hard cutoff K(α) = 7
0.8
10−2
sampling resolution 1/N
0.6
10−3 10−4 10−5
0.4 λj(C) f((λj − μα)/T) (colored by T)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
10−2
0.2 0.0
10−3
spectral index j
FIG. 2. Spectrum and soft occupations of the average normal state C. Black points show the eigenvalues λj (C) on a log scale, and the dotted line marks the sampling resolution 1/N for N = 600. The vertical dashed line marks the hard retained-mass rank K(α) = 7 for α = 0.99. Colored curves, read on the right axis, show the soft occupations f ((λj − µα )/T ) for different resolutions T . The figure shows why a hard rank is ambiguous on this graded spectrum, while QSPADE assigns fractional weights to near-cutoff modes.
inside the ladder, and no spectral gap singles out a correct number of components. Anomalies are instead low-energy states of the same chain tuned to the paramagnetic phase at hB = 1.2, drawn from a 0.7/0.2/0.1 mixture of its three lowest states; we call states from this phase the anomaly class. All reported scalar scores are averages over fresh samples from the corresponding source: the ordered-phase lowenergy mixture for normal states, the paramagnetic mixture for the anomaly class, and the TFIM ground state at the displayed field values for the sweep experiment. The evaluation protocol is detailed in Appendix H 1. We compare three residual anomaly scores:
(VII.1)
• Raw overlap residual: Qraw (ψ) = 1 − ⟨ψ|C|ψ⟩, which is density-weighted.
on n = 10 spins with J = 1, which has a quantum phase transition at the critical field hc = J between an ordered (ferromagnetic) phase with long-range ZZ correlations and a disordered (paramagnetic) phase polarized along the field. The normal source, in our setting, populates the lowenergy ladder of the ordered chain (hA = 0.4) with geometrically decreasing probabilities: each sample is an eigenstate |mj ⟩ of HTFIM (hA ) drawn with probability wj ∝ 2−j over the M = 12 lowest modes, P with a small random perturbation, so that C = N −1 i |ψi ⟩⟨ψi |, built from N = 600 such samples, is a genuinely mixed state with a graded spectrum. This thermal-like occupation makes the retained rank ambiguous by construction: the empirical eigenvalues of C (Fig. 2) decay geometrically from 0.52 through the sampling resolution 1/N ≈ 1.7 × 10−3 into the noise floor, the expected permode training counts wj N , falling from 300 to below one
• Hard top-K residual: Qhard (ψ) = 1 − ⟨ψ|PK |ψ⟩, with the rank K to be chosen.
i=1
Xi ,
10−1
temperature T
0.89 0.66 0.57
soft occupation f((λj − μα)/T)
blob moon ring
eigenvalue λj
geometry linear PCA RBF kernel-PCA QSPADE (Fourier)
i=1
• QSPADE residual: Qfα (ψ) = 1 − ⟨ψ|Mµfα ,T |ψ⟩, fixed by the single retained-mass target α = 0.99. On a graded spectrum the rank choice that hard PCA requires is ill-posed, in the sense that the data single out no particular K: the retained-mass rule maps the reasonable targets α ∈ {0.95, 0.98, 0.99, 0.995} to four different ranks K = 5, 6, 7, 8, and beyond K = 7 the cutoff reaches the sampling resolution, where the trailing empirical eigenvectors are no longer reliable. The retained-mass calibration instead places the threshold µα inside the graded tail, between λ7 and λ8 , and the soft response retains the boundary modes fractionally, a weighting that no projector can represent (Fig. 2). The resolution is chosen to resolve the spectral scale of the tail: the scalar scores
11
1.0
0.6 0.4 0.2 0.0
1 2 3 4 5 6 7 8 9 10 11 12 anomaly
0.8 0.6 0.4 0.2 0.0
10−3
normal mode index j
quoted below use T = 3 × 10−3 , and the figures display a family of resolutions around it. Figure 3 shows the consequence of this rank ambiguity mode by mode. Each hard detector is binary: for example, the mode m6 , a perfectly valid sector carrying w6 ≈ 1.6% of the average state (nine expected training samples), is scored fully anomalous at K = 4 and fully normal at K = 6, so the choice of K flips entire sectors between the two verdicts and the score carries no information about how marginal a sector is. Once the cutoff reaches the sampling resolution (K = 8, where λ8 ≈ 1/N ) the hard score mislabels individual modes erratically, because the trailing eigenvectors mix directions whose eigenvalues are degenerate at the noise level, so increasing K is not a remedy. Instead, the soft residual, with one α and no rank choice, grades the same ladder monotonically: well-sampled modes score near zero, boundary modes take intermediate values that track their rarity, and only the modes at or below the sampling resolution, which no data-driven detector could certify, score as anomalous. The rare-but-valid sector, the 1.6% mode m6 above, makes the comparison concrete (Fig. 4). The densityweighted Qraw assigns it a mean score of 0.99, above the genuine near-critical anomalies at 0.86; the hard projector at K = 4 commits the same inversion (1.00 against 0.45), while at K = 6, whose cutoff retains exactly the first six modes, it absorbs the sector silently, with no trace of its rarity. The calibrated Qfα keeps the rare sector low (0.02 at T = 3 × 10−3 ) and clearly below the anomaly (at 0.40), provided the resolution resolves the tail: for T much larger than the spacing of the tail eigenvalues the detector over-retains the tail and this contrast
dominant normal (50%)
rare normal (1.6%)
anomaly (paramagnetic)
10−3
FIG. 4. Mean residual anomaly score for a dominant normal mode (50%), a rare normal mode (1.6%), and the paramagnetic anomaly class. The density-weighted Qraw and the hard K = 4 projector score the rare valid mode above the anomaly, while the hard K = 6 projector gives it zero score and hides its rarity. QSPADE assigns the rare mode an intermediate, temperature-dependent score; when T resolves the spectral tail, the rare valid mode remains below the anomaly. This shows soft scoring in QSPADE separates rarity within the normal support from genuine support mismatch.
0.8
anomaly score Q
FIG. 3. Mean residual anomaly score for each normal mode j, with the expected number of training samples wj N shown on the top axis, and for the anomaly class in the rightmost column. Grey dashed curves are hard top-K projectors; colored curves are QSPADE residuals at different resolutions T . Hard ranks flip whole normal modes between zero and unit score, and become unstable when the cutoff reaches the samplingresolution regime (K = 8). QSPADE replaces this all-ornothing behavior by a graded transition from well-sampled modes to unresolved tail modes.
10−2
QQSPADE (colored by T) Qhard (T → 0 limit) Qraw (naive) critical point hc
0.6 0.4 0.2 0.0
0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00
10−2
temperature T
10−2
mean anomaly score
0.8
Qhard, K = 4 Qhard, K = 6 Qhard, K = 8 QQSPADE (colored by T)
temperature T
mean anomaly score
1.0
Qraw Qhard (K = 4) Qhard (K = 6) QQSPADE (colored by T)
temperature T
expected samples per normal mode
300 150 75 38 19 9 5 2 1 1 0 0
10−3
transverse field h/J
FIG. 5. QSPADE residual score of the TFIM ground state as the transverse field h is swept across the transition. The score is near zero at the normal training field hA = 0.4 (grey line), stays low in the ordered phase, and rises through the critical point hc = 1 (dashed) into the paramagnetic phase. Colored curves show different resolutions T ; as T → 0 they approach the hard retained-rank projector K(α) (black markers), while larger T gives a softer monitor. The density-weighted Qraw (dotted) is offset upward even at the normal field.
degrades. This finite resolution and the smooth occupation profile are precisely what the spectral-threshold measurement of Section VI A produces, fixed by one retainedmass target rather than a discrete rank ladder, and the hard projector at the retained-mass rank is recovered as its T → 0 endpoint. Sweeping the transverse field across the critical point
12 displays the whole soft family at once (Fig. 5). The QSPADE residual is smallest at the normal training field hA = 0.4, stays low throughout the ordered phase (h < hc ), and rises through the transition into the paramagnetic phase (h > hc ), so anomaly detection tracks the quantum phase boundary. The family interpolates smoothly between a permissive large-T monitor and the hard projector recovered as T → 0.
VIII.
CONCLUSION
We introduced QSPADE, a measurement-based framework for PCA-style anomaly detection with quantum data. The main contribution is to replace the reconstruction-first PCA workflow, including forming a covariance or Gram matrix, recovering principal eigenvectors, and selecting a hard rank, with directly getting an anomaly score from the average normal state C. For the anomaly score, the residual Q score is the basic one-class detector and requires only retained-mass calibration, while the Hotelling-type T 2 score can be added when component-resolved leverage information is needed. The soft spectral threshold is what makes this formulation different from standard PCA. Instead of accepting or rejecting whole eigenspaces by a discrete rank cutoff, QSPADE lets boundary spectral components contribute fractionally, with the hard-projector PCA score recovered in the zero-temperature limit T → 0. This soft weighting reduces the sensitivity to eigenvalue crossings or small spectral gaps near the cutoff, since near-threshold modes are not forced into a binary rank decision. It also avoids the need to resolve adjacent near-cutoff eigenvalues solely to decide a hard rank. As a result, the default residual detector has dimension-independent measurement-shot complexity, while the optional leverage calibration is con-
[1] I. T. Jolliffe, Principal Component Analysis, 2nd ed., Springer Series in Statistics (Springer, 2002). [2] H. Hotelling, The generalization of Student’s ratio, The Annals of Mathematical Statistics 2, 360 (1931). [3] J. E. Jackson and G. S. Mudholkar, Control procedures for residuals associated with principal component analysis, Technometrics 21, 341 (1979). [4] S. Corli, L. Moro, D. Dragoni, M. Dispenza, and E. Prati, Quantum machine learning algorithms for anomaly detection: A review, Future Generation Computer Systems 166, 107632 (2025). [5] N. Liu and P. Rebentrost, Quantum machine learning for quantum anomaly detection, Physical Review A 97, 042315 (2018). [6] S. Lloyd, M. Mohseni, and P. Rebentrost, Quantum principal component analysis, Nature Physics 10, 631 (2014). [7] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics, in Proceedings of the 51st Annual ACM SIGACT Symposium
e fined to the finite-window overhead O(Ld). The numerical experiments demonstrate the two intended roles of the framework. On encoded classical data, QSPADE reproduces kernel-PCA-like nonlinear decision regions. On quantum-native transverse-field Ising data, it detects the change across the phase transition without specifying an order parameter. These results support QSPADE as both a quantum-kernel anomaly detector and a data-driven monitor for quantum native systems. Our results also raise concrete implementation questions. A natural next step is a small-scale hardware demonstration of QSPADE. The qumode construction is well matched to Gaussian responses and spectral sweeps, but the required qumode–qubit coupling may make the Hadamard-test route the more immediate option on qubit devices. At the circuit level, density matrix exponentiation gives a standard sample-access primitive, but tailored block-encodings or compiled simulations of C may reduce the overhead, which deserves further exploration. Finally, the same calibrated support viewpoint could be extended to multi-class quantum state discrimination, with separate learned spectral supports for different classes.
ACKNOWLEDGEMENTS
MMW acknowledges support from the National Science Foundation under grant no. 2611810. NL acknowledges funding from the Science and Technology Commission of Shanghai Municipality (STCSM) grant no. 24LZ1401200 (21JC1402900), NSFC grants no. 12471411 and no. 12341104, the Shanghai Jiao Tong University 2030 Initiative, the Shanghai Pilot Program for Basic Research, and the Fundamental Research Funds for the Central Universities.
on Theory of Computing (STOC) (2019) pp. 193–204. [8] G. Park, J. Huh, and D. K. Park, Variational quantum one-class classifier, Machine Learning: Science and Technology 4, 015006 (2023). [9] H. Oh and D. K. Park, Quantum support vector data description for anomaly detection, Machine Learning: Science and Technology 5, 035052 (2024). [10] K. Kottmann, F. Metz, J. Fraxanet, and N. Baldelli, Variational quantum anomaly detection: Unsupervised mapping of phase diagrams on a physical quantum computer, Physical Review Research 3, 043184 (2021). [11] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018). [12] S. M. Hossain and S. Bhattacharjee, Filtered spectral projection for quantum principal component analysis (2026), arXiv:2603.13441 [stat.ML]. [13] Y. Yuan, M. Minervini, M. M. Wilde, and N. Liu, Quantum principal component analysis without eigenvector re-
13 covery (2026), arXiv:2605.27942 [quant-ph]. [14] S. Kimmel, C. Y.-Y. Lin, G. H. Low, M. Ozols, and T. J. Yoder, Hamiltonian simulation with optimal sample complexity, npj Quantum Information 3, 13 (2017). [15] B. Go, H. Kwon, S. Park, D. Patel, and M. M. Wilde, Sample-based Hamiltonian and Lindbladian simulation: Non-asymptotic analysis of sample complexity, Quantum Science and Technology 10, 045058 (2025). [16] N. Liu and M. M. Wilde, Fermi-dirac thermal measurements: A framework for quantum hypothesis testing and semidefinite optimization (2026), arXiv:2603.04061 [quant-ph]. [17] C. Weedbrook, S. Pirandola, R. Garcı́a-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Reviews of Modern Physics 84, 621 (2012). [18] A. He, N. Liu, and M. M. Wilde, Fermi–Dirac machines as quantizations of neurons (2026), arXiv:2605.24386 [quant-ph]. [19] A. Dvoretzky, J. Kiefer, and J. Wolfowitz, Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator, The Annals of Mathematical Statistics 27, 642 (1956). [20] S. Boucheron, G. Lugosi, and P. Massart, Concentration Inequalities: A Nonasymptotic Theory of Independence (Oxford University Press, 2013). [21] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Quantum algorithms revisited, Proceedings of the Royal Society A 454, 339 (1998). [22] G. H. Low and I. L. Chuang, Optimal Hamiltonian simulation by quantum signal processing, Physical Review Letters 118, 010501 (2017). [23] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating Hamiltonian dynamics with a truncated Taylor series, Physical Review Letters 114, 090502 (2015). [24] W. Hoeffding, Probability inequalities for sums of bounded random variables, Journal of the American Statistical Association 58, 13 (1963). [25] H. Robbins and S. Monro, A stochastic approximation method, The Annals of Mathematical Statistics 22, 400 (1951). [26] D. Ruppert, Efficient estimations from a slowly convergent Robbins–Monro process, Tech. Rep. Technical Report No. 781 (School of Operations Research and Industrial Engineering, Cornell University, Ithaca, NY, 1988). [27] B. T. Polyak and A. B. Juditsky, Acceleration of stochastic approximation by averaging, SIAM Journal on Control and Optimization 30, 838 (1992). [28] R. Bhatia, Matrix Analysis, Graduate Texts in Mathematics, Vol. 169 (Springer, New York, 1997).
Appendix A: Proof of Proposition III.2
Proof. For any effect 0 ⪯ M ⪯ I, let QM be the orthogonal projector onto the range of M . Since 0 ⪯ M ⪯ QM and rank(QM ) = rank(M ) ≤ K, Tr(M C) ≤ Tr(CQM ).
(A.1)
Ky Fan’s maximum principle [28] gives Tr(CQM ) ≤
K X
λj .
(A.2)
j=1
Equality is achieved by M = PK . If λK > λK+1 , the top-K spectral subspace is unique.
Appendix B: Proof of Theorem IV.2
Proof. First, we establish various properties of the generator ϕf . By the fundamental theorem of calculus and the inverse-function theorem, its derivatives are ϕ′f (m) = f −1 (m),
ϕ′′f (m) =
1 f ′ (f −1 (m))
> 0. (B.1)
Since f ′ > 0, the inverse f −1 is strictly increasing, making ϕf strictly convex on (0, 1). Integrability gives a finite continuous extension at the endpoints. The definition yields ϕf (0) = 0, and the normalization condition in R1 Definition IV.1 ensures ϕf (1) = 0 f −1 (u) du = 0. Consequently, the objective function in (IV.2) is strictly convex on the operator interval because M 7→ Tr ϕf (M ) is strictly convex and the remaining terms are linear. Since f maps R onto (0, 1), the derivative ϕ′f (m) = f −1 (m) tends to −∞ as m ↓ 0 and to +∞ as m ↑ 1. Thus, the optimum must lie in the open interval 0 ≺ M ≺ I. The Fréchet first-order condition for the objective is µI − C + T ϕ′f (M ) = 0.
(B.2)
Rearranging the terms yields ϕ′f (M ) =
C − µI . T
(B.3)
Applying the inverse derivative (ϕ′f )−1 = f via functional calculus gives (IV.3). Strict convexity guarantees that this optimal solution is unique.
Appendix C: Proof of Proposition IV.3
Proof. For each eigenvalue λj , the scalar argument (λj − µ) /T tends to +∞, −∞, or 0 depending on whether λj > µ, λj < µ, or λj = µ, respectively. The result follows from the boundary conditions in Definition IV.1.
Appendix D: Proof of Proposition V.2
Proof. By Proposition IV.3, each tail detector Mℓ converges to the hard projector onto the modes above threshold µℓ . Rank-calibrated thresholds isolate consecutive
14 nondegenerate eigenmodes in the sharp limit, so that Dj = Mj − Mj−1 → |uj ⟩⟨uj |. Substitution into (V.4) and (V.12) gives the residual and leverage limits, because pσ (j) → |aj |2 and vj = Tr(Dj C) → λj .
Appendix E: Response examples for qumode spectral-threshold measurement
This appendix records two concrete optimizationadmissible responses used in the qumode construction of Section VI A. In both cases the symmetry f (−x) = R1 1 − f (x) implies 0 f −1 (m) dm = 0, so the centering convention in Definition IV.1 is satisfied.
1.
Fermi–Dirac response
The increasing Fermi–Dirac, or logistic, response is fFD (x) :=
1 . 1 + e−x
(E.1)
Its quantile and generator are m , 1−m ϕFD (m) = m log m + (1 − m) log(1 − m). −1 fFD (m) = log
(E.2) (E.3)
where erf(·) denotes the error function. Its derivative is the standard normal density 2 1 ′ fG (x) = √ e−x /2 . 2π
(E.9)
The quantile and generator are −1 fG (m) = Φ−1 (m),
(E.10)
2 1 1 −1 ϕG (m) = − √ exp − Φ (m) . 2 2π
(E.11)
The reflected position density is the Gaussian density q2 1 exp − 2 , (E.12) gG,T1 (q) = √ 2T1 2πT1 so the control wavepacket is the Gaussian qumode state Z q2 1 |ψG,T1 ⟩ = exp − |q⟩ dq. (E.13) 2 1/4 4T12 R (2πT1 ) Again T = T1 T2 and µ = βT2 , and the implemented effect is C − µI G (E.14) Eβ,T = Φ 1 ,T2 T 1 C − µI = . (E.15) I + erf √ 2 2T
Its derivative is the logistic probability density ′ fFD (x) =
x e−x 1 = sech2 . −x 2 (1 + e ) 4 2
Appendix F: Proof of Proposition VI.1
(E.4)
Therefore the reflected position density in (VI.1) is 1 q 2 gFD,T1 (q) = sech , (E.5) 4T1 2T1 and the corresponding control wavepacket is Z q 1 √ sech |ψFD,T1 ⟩ = |q⟩ dq. 2T1 R 2 T1
(E.6)
With T = T1 T2 and µ = βT2 , Proposition VI.1 gives the effect −1 C − µI FD Eβ,T = I + exp − . (E.7) 1 ,T2 T
P Proof. Let C = j λj |uj ⟩⟨uj |. Under e−ip̂λj /T2 , the position distribution is shifted by λj /T2 . The density conditioned on eigenmode |uj ⟩ is given by λj gf,T1 q − . (F.1) T2 The tail probability is Z ∞ q − λj /T2 1 ′ f − dq T1 T1 β Z (λj /T2 −β)/T1 = f ′ (z) dz −∞ λj − βT2 =f . T1 T2
(F.2) (F.3)
Functional calculus gives (VI.4). 2.
Gaussian response Appendix G: Proof of Proposition VI.3
The Gaussian, or probit, response is the standard normal cumulative distribution function 1 x := √ fG (x) = Φ(x) 1 + erf (E.8) 2 2
Proof. For a sampled ω, write χω = |χω |eiφω . The identity shift factors as eitω (C−µI)/T = e−iµtω /T eitω C/T .
(G.1)
15 kernel-PCA (polynomial) AUC = 0.799 1.0
0.8
0.8
0.8
θω (µ) = φω − µtω /T
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
Set Yµ (ρ) = bf + Af |χω |Zω . Taking the expectation over the measurement outcome from the Hadamard test and over ω gives EYµ (ρ) h i = bf + Af Eω ℜ χω e−iµtω /T Tr eitω C/T ρ C − µI . = Tr ρ fe T
kernel-PCA (Laplacian) AUC = 0.891 1.0
0.8
QSAD (degree-two polynomial) QSAD (random features) AUC = 0.764 AUC = 0.688 1.0 1.0 1.0 0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
QSAD (Fourier) AUC = 0.879
0.6
0.4
0.2
0.0
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
(G.4) (G.5)
The uniform scalar approximation error implies an operator-norm error at most εf , which proves (VI.23). The range bound follows from |χω | ≤ 1 and |Zω | = 1. Appendix H: Supplementary numerical results
This appendix collects supporting material for the experiments of Section VII: the feature maps and simulation parameters, the resolution–shot trade-off, and the kernel and geometry comparisons.
1.
1.0
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
(G.3)
x2
outputs a sign variable Zω ∈ {−1, +1} satisfying h i E[Zω | ω] = ℜ eiθω (µ) Tr eitω C/T ρ .
(G.2)
kernel-PCA (RBF) AUC = 0.873
Q percentile
1.0
x2
A Hadamard test with controlled unitary eitω C/T and control phase
Experiments setup
In the classical-data experiment of Section VII A, each classical point x = (x1 , x2 ) ∈ [0, 1]2 is mapped to a feature state |ϕ(x)⟩ by one of the maps below; we write x̃i = 2xi −1. All maps produce a proper pure state of unit Euclidean norm: the proportionality sign in ϕ ∝ (· · · ) denotes this normalization (the amplitude vector is divided by its norm). The Fourier map is used in the main text as a fixed representative; the comparison across feature maps and classical kernels appears in Fig. 6. • Degree-two polynomial : ϕ ∝ (1, x̃1 , x̃2 , x̃21 , x̃1 x̃2 , x̃22 ). • Fourier (four qubits, 16-dimensional): the tensor product over the two coordinates of the peraxis features (cos πt, cos 2πt, sin πt, sin 2πt), t ∈ {x1 , x2 }. • Gaussian random features (four qubits): ϕ ∝ cos(W x + b) with Wij ∼ N (0, 2ξ) and bi ∼ Unif[0, 2π), a 16-dimensional Monte-Carlo approx′ 2 imation of the RBF kernel e−ξ∥x−x ∥ . The classical kernel-PCA baselines use the RBF kernel ′ 2 e−ξ∥x−x ∥ (the baseline of the main text), the Laplacian
x1
x1
x1
FIG. 6. Residual percentile on the moon dataset for three classical kernel-PCA baselines (top row) and QSPADE with three feature maps (bottom row), with the one-class detection AUC (against a uniform background) in each panel title.
−1 ′ kernel e−ξ1 ∥x−x ∥1 with ξ1 = mediani,j ∥xi − xj ∥1 , and the cubic polynomial kernel (1 + x · x′ )3 . For the classical-data experiment, the normal clouds are a single moon (and, for the geometry sweep, an elongated Gaussian blob and a ring), each with N = 450–500 points rescaled into [0, 1]2 . The RBF kernel-PCA bandwidth is set by the median heuristic ξ = 2 mediani,j ∥xi − −1 xj ∥ 2 and its retained rank by the same explainedvariance target as QSPADE’s retained mass; the QSPADE detector uses α = 0.88, T = 0.1 in the threemethod grid of Fig. 1, α = 0.88, T = 0.035 in the kernel comparison of Fig. 6, and α = 0.90, T = 0.05 in the geometry sweep underlying Table II. For the quantum-native experiment, the normal ladder and anomalies are as in Section VII B (M = 12 modes with weights ∝ 2−j , N = 600 training samples, perturbation amplitude 0.05), with α = 0.99 and, unless a family of resolutions T ∈ [10−3 , 3×10−2 ] is shown, T = 3×10−3 . The expected per-mode training counts wj N fall from 300 at j = 1 to nine at j = 6, about one at j = 9, and below one for the last modes. Sector and per-mode scores are means over 300 states sampled as in training (a single normal mode per sector, the 0.7/0.2/0.1 paramagnetic mixture for the anomaly class); the detection AUC compares 600 normal-mixture samples against 300 anomalous ones; and the field sweep of Fig. 5 tracks the single ground state at each field h. All scores are converted to normal-calibrated percentiles through the empirical cumulative distribution of the normal scores; detection AUC uses uniformbackground anomalies (classical) or paramagnetic-phase states (quantum-native). Centering of the training and test data follows each method’s standard convention: linear PCA centers the inputs with the training-set mean, the kernel-PCA baselines center in feature space through
16 blob
1.0
moon
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
ring 1.0
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
0.4
0.4
0.4
Q percentile
Classical KPCA (RBF)
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 1.0 1.0 1.0
0.2
0.2
0.2
0.4
QSAD (fourier)
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 1.0 1.0 1.0 0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
blob, moon, and ring geometries, that underlies Table II: linear PCA cannot follow the curved or non-convex supports, most strikingly the ring, where its straight acceptance band runs through the empty centre, whereas the RBF kernel-PCA baseline and QSPADE both track the manifold. 1.0
0.9
detection AUC
Classical PCA
1.0
0.2
0.8 0.7 0.6 0.5
∞ shots 200 shots 50 shots
10−3 0.0
0.0 0.0 0.0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
FIG. 7. Residual-percentile decision regions for linear PCA (top), RBF kernel-PCA (middle), and QSPADE (bottom) on the blob, moon, and ring. Linear PCA cannot follow the curved or non-convex supports; both nonlinear monitors do.
Figure 6 places the classical and the quantum side of the kernel-PCA correspondence of Section VII A side by side on the moon: classical kernel-PCA with polynomial, RBF, and Laplacian kernels in the top row, and QSPADE with the degree-two polynomial, random-feature, and Fourier maps in the bottom row, under the same calibration and with the one-class AUC in each panel. The columns pair related kernels. The two polynomial methods produce similarly broad acceptance regions, the weakest separation among the classical baselines; the random-feature map is a sixteen-dimensional approximation of the exact RBF kernel above it, visibly coarser at this encoding size and the weakest panel on the QSPADE side; and the Laplacian and Fourier panels give the tightest boundaries on each side. Because the classical-data experiment is a consistency check rather than an accuracy claim, the encoding is not tuned for performance; the Fourier map serves as a fixed representative in the main text, with the comparison across kernels and maps reported here. Figure 7 shows the three-method comparison across the three synthetic classical datasets, the
10−2
QSPADE temperature T
10−1
FIG. 8. Detection ROC-AUC versus resolution T at infinite and finite measurement-shot budgets for the quantumnative experiment (n = 10): normal-mixture samples from the graded ladder at hA = 0.4 versus paramagnetic anomalies at hB = 1.2. Smaller T is markedly more robust to shot noise; the dotted line is the density-weighted Qraw baseline.
2.
the standard double-centering of the Gram matrix, and the feature maps act on the mean-centered data.
20 shots Qraw (naive)
Resolution–shot trade-off (quantum-native experiment)
The QSPADE score is read off finite measurement shots: the acceptance ⟨ψ|Mµfα ,T |ψ⟩ is estimated from m binary accept/reject outcomes, so the empirical residual carries binomial sampling noise of order m−1/2 . The data and calibration are those of the quantum-native experiment of Section VII B: the detector is trained on the graded normal ladder of the ordered TFIM at hA = 0.4 and calibrated at α = 0.99, and the AUC separates fresh normal-mixture samples from paramagnetic anomalies at hB = 1.2. Figure 8 plots the detection ROC-AUC against T for several shot budgets, the ∞-shot curve being the exact score. A sharper detector (small T ) places the normal and anomalous acceptances near the well-separated extremes 0 and 1, where the per-state variance is small, and stays robust down to a few tens of shots; a softer detector (large T ) compresses the acceptances, so the normal–anomaly gap shrinks relative to the shot noise and the AUC degrades. This is the cost side of the resolution knob: the larger T that yields a more graded, information-rich score also demands more measurement shots to resolve.