arXiv CS open access

Supervised Learning with Quantum Hamiltonian Embeddings

Riccardo Molteni et al. · arxiv_cs

arXiv CS · Papers

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Supervised Learning with Quantum Hamiltonian Embeddings Riccardo Molteni,1, 2, ∗ Minjun Jeon,1, 3, ∗ Paolo Braccia,1 Abhinav Deshpande,1 Vedran Dunjko,2 Liang Jiang,1, 4 Su-un Lee,1, 4 Swarnadeep Majumder,1 Javier Robledo Moreno,1 Parthasarathy Suryanarayanan,1 and Kunal Sharma1, † 1 IBM Research 2 Applied Quantum Algorithms ⟨aQa L ⟩, Leiden University, The Netherlands

3 Department of Materials, University of Oxford, Parks Road, Oxford OX1 3PH, United Kingdom 4 The University of Chicago

arXiv:2610.26833v1 [quant-ph] 19 Oct 2026

(Dated: October 20, 2026)

Machine learning is widely regarded as a promising application in which quantum computing could offer a computational advantage. However, general quantum machine learning models with provable guarantees on performance, trainability, and efficient inference remain limited. In this work, we introduce quantum Hamiltonian embeddings, a supervised learning framework in which a classical training dataset is encoded into a single effective, parameterized Hamiltonian. Training optimizes the parameters of this data-dependent Hamiltonian, whose approximate ground state defines the classification rule. Once the ground state has been approximated using a quantum algorithm, inference can be performed either classically using a classical representation of the learned state, or with a quantum representation of the state. We establish efficient methods for gradient estimation and prove favorable generalization guarantees. Furthermore, under standard cryptographic assumptions, we show that Hamiltonian embeddings can realize hypothesis functions beyond the reach of efficient classical algorithms. We then demonstrate the framework at scale on IBM Heron processors using a standard benchmark dataset. Together, these results establish Hamiltonian embeddings as a new framework for quantum machine learning, combining theoretical guarantees with a path for experimental investigations at scale on quantum hardware.

CONTENTS

I. Introduction

II. Hamiltonian embedding model A. Construction of the model 1. Choice of the observable 2. Choice of the unitary embedding 3. Choice of the loss function 4. Ground-state approximation B. Efficient gradient estimation and absence of barren plateaus

4 5 5 7 8 10 10

III. Properties of the Hamiltonian embedding model A. Generalization bounds B. Classical inference and surrogate models 1. Fourier surrogates 2. Pauli surrogates 3. Deep Clifford encodings 4. Dequantizability of the model C. Encoding a classically hard problem in the Hamiltonian embedding model

12 12 14 14 14 15 16 17

IV. Experiments A. The model instance used in the experiments B. Results

18 18 19

V. Discussion

Acknowledgments ∗ The first two authors contributed equally to this work. † kunals@ibm.com

2 References

A. Energy-based generalization bounds 1. Lipschitz assumption

23 27

B. Gradient evaluation and the Hellmann–Feynman relation 1. Gradients with approximate ground states

28 29

C. Absence of barren plateaus in the Hamiltonian embedding model

D. Universality of the Hamiltonian Embedding model

E. Solving a classically hard problem with Hamiltonian embedding 1. The discrete cube root problem 2. Definition of the learning problem 3. Classical hardness of the learning problem 4. Quantum solution with Hamiltonian embedding

35 35 36 37 37

F. A construction with highly entangled ground states 1. Entangled graph states as ground states 2. Entanglement properties of the ground state 3. Numerical demonstration

41 41 42 42

G. A provable separation between Pauli surrogates and Hamiltonian embedding 1. Idealized embedding Hamiltonian and its ground state: perfect classification 2. Embedding Hamiltonian with finite training data a. Empirical accuracy of Hamiltonian embedding b. Accuracy of classical surrogates 3. Required bond dimension for good classification with Hamiltonian embedding

43 45 46 47 48 50

H. Numerical setup 1. Encoding circuit and embedding Hamiltonian 2. Classical iterative solver 3. Optimization and model selection

51 51 52 53

I. Image classification benchmark 1. Protocol 2. Experimental Baselines

53 53 53

J. Hardware experiments 1. The hardware solver 2. Model training insights

54 54 54

K. A numerical counterexample to Pauli classical surrogate models 1. Construction of the learning task 2. Numerical experiment

55 56 58

3 I.

INTRODUCTION

The prospect of quantum computational advantage in learning tasks has driven significant interest in quantum machine learning (QML) [BWP+ 17, CC25]. For instance, in supervised learning, multiple architectures have been proposed [HCT+ 19, SP21, PCGL20, RML14, JGM+ 24] and provable advantages have been established in particularly constructed settings [LAT21, AG26, MGD26, ZZN+ 26, GGPD25]. A broader goal in QML is to develop general models that are efficiently trainable, admit generalization guarantees, cannot be efficiently simulated classically in general, while also being implementable at scale on quantum processors. However, the difficulty of simultaneously satisfying these requirements is particularly evident for parameterized quantum circuits [BLSF19, CAB+ 21, SBSW20]. Shallow-depth PQCs are well suited for current quantum processors but can often be simulated efficiently on classical computers [CLG+ 25, GGPD25], whereas deeper circuits may be difficult to optimize because of barren plateaus [MBS+ 18, CSV+ 21, LTW+ 25]. Similarly, beyond a few examples [GGC+ 24, HBK+ 25], identifying effective quantum kernels for more general datasets remains challenging. In this work, we introduce Hamiltonian embedding, a classification algorithm that can be efficiently implemented on a quantum processor, in which the classifier is constructed by solving a physically motivated energy-minimization problem. The method encodes the entire training dataset into a single parameterized Hamiltonian. Each training example contributes a term obtained by dressing a fixed, label-dependent operator, namely by conjugating it with a parameterized unitary that encodes the corresponding input. Then training optimizes the encoding parameters to lower the ground-state energy of this data-dependent Hamiltonian, and an approximate ground state of the optimized Hamiltonian defines the classifier. In particular, for a new input, the trained encoding specifies an observable whose expectation value with respect to the learned ground state determines the predicted label. We identify tasks for which the physical constraints imposed by the learned quantum state provide a useful inductive bias, yielding a provable reduction in sample complexity relative to classically regularized surrogate models. We then show that the Hamiltonian embedding framework has several properties desirable for QML. During optimization, a quantum algorithm prepares approximations to the required ground states. In the case where the encoding unitaries allow the expectation values of the dressed observables to be efficiently computed from a classical representation of these ground states, both the loss and its gradients with respect to the trainable parameters can be evaluated classically. In particular, we show that the gradients can be expressed as expectation values and evaluated using a parameter-shift rule. Shallow local encoding unitaries provide one such setting, enabling classical inference after training. For a general family of these encodings, we also show the model avoids barren plateaus for a relevant subset of parameters. Importantly, in this regime the model is not known to be classically simulable, as it involves approximating ground states of local Hamiltonians. Although the shallow encoding may allow for classical surrogates, we exhibit a task on which a specific model instance provably achieves substantially higher test accuracy than a corresponding natural surrogate, at the same training-set size. Moreover, even with shallow encoding unitaries, we show that the Hamiltonian embedding framework can implement universal approximators for measurable binary target functions. Specifically, for any such function and any desired accuracy, we provide a construction using finitely many qubits such that the hypothesis function implemented by the model achieves an arbitrarily small average approximation error under the uniform input distribution. Irrespective of the encoding depth, we also show the model admits favorable generalization guarantees for a wide class of encoding unitaries. If training yields a sufficiently low-energy approximate ground state, matrix Hoeffding inequalities and covering number arguments establish generalization to unseen inputs. We show that for fixed generalization error, the required dataset size scales linearly, up to logarithmic factors, with the number of qubits and variational parameters. Finally, the model can realize hypothesis functions that lie beyond the reach of efficient classical algorithms. Although preparing ground states of local Hamiltonians is hard in the worst case, even for quantum computers, we identify structured instances that encode classically hard learning problems for which the relevant ground states can nevertheless be efficiently approximated on a quantum computer. Under standard cryptographic assumptions, these instances yield explicit learning problems solvable on a quantum processor within the Hamiltonian embedding framework, while remaining classically hard. We note that restricting the encoding unitaries to those for which the relevant dressed-observable expectation values can be evaluated classically places Hamiltonian embedding within the shadow model framework introduced in [JGM+ 24]. A shadow QML model uses a quantum computer during training to obtain a quantum advice state, whose classical representation then supports entirely classical inference. Whereas Ref. [JGM+ 24] introduces shadow models as a general framework, our construction provides an explicit method to obtain the advice state, namely as the ground state of the embedding Hamiltonian. Ground states of data-dependent Hamiltonians have also been used for classification in [YGH+ 25]. There, however, the encoding is fixed and the algorithm consists only of preparing a low-energy state of the resulting data-

4 encoding Hamiltonian, without optimizing over a trainable family of Hamiltonians. The Hamiltonian embedding framework introduced here instead learns both the advice state and the Hamiltonian whose ground state realizes it. To assess whether the model can be trained to achieve low validation and test error, we study it in numerical experiments on a standard ML task and implement it on IBM Heron processors. The results show successful optimization and competitive predictive performance relative to the quantum and classical baselines considered. The manuscript is organized as follows. • In Section II, we introduce the Hamiltonian embedding model and provide a complete description of the training and prediction procedures. We also discuss the choice of the embedding unitaries and observables, derive an efficient gradient-estimation procedure, and establish the absence of barren plateaus for a relevant subset of parameters. • The main theoretical results are in Section III. We derive analytical generalization bounds and study classical surrogates of the trained model. We also identify learning tasks for which the physical constraints imposed by the learned quantum ground state act as a useful inductive bias, yielding improved sample efficiency over classically regularized surrogate models. Finally, we exhibit a learning task that is efficiently solvable within the Hamiltonian embedding framework but classically hard under a standard cryptographic assumption. • In Section IV, we present numerical results for the binary Fashion-MNIST classification task. We benchmark the Hamiltonian embedding framework on systems up to 50 qubits against quantum kernels studied in Ref. [HBM+ 21] and classical machine-learning baselines. Across this range, the model remains trainable and achieves competitive test performance. We also implement the training procedure on IBM quantum processors, using hardware-generated samples in the Sample-based Krylov Quantum Diagonalization ground-state solver [YMI+ 25, MN25, SKI+ 25], and demonstrate successful training on a 50-qubit instance.

II.

HAMILTONIAN EMBEDDING MODEL

Although our framework can be readily extended to multiclass classification, here we address the problem of N , where each input x is sampled from binary classification. In particular, we are given a training set T = {( xℓ , yℓ )}ℓ= ℓ 1 d a distribution D over Ω ⊆ R , and the corresponding label is generated by an unknown function f : Ω → {0, 1}, such that yℓ = f ( xℓ ). The goal of the learning algorithm is to infer an approximate map f˜ that, with high probability, agrees with the true mapping f on previously unseen test data sampled from the same distribution D . Given a dataset T , the model constructs a family of dataset-dependent hypothesis functions as follows. On an n-qubit Hilbert space, we define an encoding unitary U ( x, α) parameterized by both the input x and trainable parameters α ∈ R M , together with a family of label-dependent operators {O(y)}y , assigning a measurement operator O(y) to each class label y ∈ {0, 1}. Then, for each data point z = ( x, y), with y = f ( x ), we define a corresponding Hamiltonian term as Hz (α) = U † ( x, α)O(y)U ( x, α). The data Hamiltonian on the whole dataset is defined as 1 HT (α) = N

∑

ℓ=1

I − U † ( xℓ , α)O(yℓ )U ( xℓ , α)

! (1)

= I−

1 U † ( xℓ , α)O(yℓ )U ( xℓ , α) N∑ ℓ

(2)

= I−

1 Hzℓ (α) . N∑ ℓ

(3)

For each α, we define |ψT ,α ⟩ to be the ground state of HT (α). The state |ψT ,α ⟩ is then used to construct the corresponding hypothesis function implemented by the model. In particular, given a new input point x ′ , we apply the encoding unitary U ( x ′ , α) to |ψT ,α ⟩ and measure it with respect to O(y), for y ∈ {0, 1}. As we show in Section II A 1, we can choose O(0) + O(1) = I, so that the expectation values of these positive operators on the encoded state define a normalized probability distribution over the labels. Specifically, the probability PψT ,α ,α ( x ′ , y) ∈ [0, 1] assigned to label y for input x ′ is given by the following expression: h i PψT ,α ,α ( x ′ , y) = Tr O(y) U ( x ′ , α)|ψT ,α ⟩⟨ψT ,α |U † ( x ′ , α) .

(4)

5 The corresponding hypothesis function assigns a label ỹ ∈ {0, 1} to the input x ′ as follows: f˜T ,α ( x ′ ) = arg max PψT ,α ,α ( x ′ , y).

(5)

y∈{0,1}

For a given training set T , the resulting family of candidate hypotheses is therefore HT = { f˜T ,α : α ∈ R M }. We note that the family of hypothesis functions depends on the training data through the ground states of HT (α). Moreover, HT is a subset of the broader, data-independent family of hypothesis functions obtained by allowing arbitrary normalized n-qubit states in the same classifier construction. The model selects a hypothesis from HT by choosing α according to the following training procedure. Let us consider the case when 0 ⪯ O(y) ⪯ I for every y ∈ {0, 1}, then 0 ⪯ N1 ∑ℓ Hzℓ (α) ⪯ I. This implies that the energy of HT (α) is always lower bounded by zero for any α. In particular, the energy can reach the minimum value zero if the ground state |ψT ,α ⟩ of HT (α) is such that ⟨ψT ,α | U † ( x, α)O(y)U ( x, α) |ψT ,α ⟩ = 1 for every training point ( xℓ , yℓ ) in T . The optimal parameter vector α⋆ is then chosen to minimize the ground-state energy over the family of Hamiltonians { HT (α)}α . As we see below, this guarantees the correct classification of the training data. Moreover, in Section III A, we show that our model also satisfies favorable generalization bounds. In particular, let us define the empirical energy cost function as: bT (ψ, α) ≡ Ez∈T ℓenergy (ψ, Hz (α)) = ⟨ψ| HT (α)|ψ⟩, (6) C where ℓenergy is the energy loss, defined as:

ℓenergy (ψ, Hz (α)) = 1 − ⟨ψ| Hz (α)|ψ⟩.

(7)

Then the training objective is to find the parameter value α⋆ and the state ψT ,α⋆ that minimize the cost function bT (ψ, α). We note that the algorithm may not reach the exact ground state of the training Hamiltonian HT (α⋆ ). In C other words, the learned state ψT ,α⋆ may only be an approximate ground state. However, as we show in Section III A, it is sufficient for the learned state to achieve a low energy with respect to the training Hamiltonian: a small empirical energy cost already guarantees good performance on unseen input points. In practice, this optimization is carried out through an alternating nested procedure. Starting from an initial value of α, the algorithm first computes the corresponding ground state of HT (α) and then updates α to decrease the cost. Repeating this procedure yields a sequence of parameter values and associated ground states until convergence to a final pair (α⋆ , ψT ,α⋆ ). We detail the whole learning procedure in Algorithm 1. Although reaching zero energy loss guarantees correct classification bT as α is updated does not necessarily lead to a on the training dataset, a decrease in the empirical energy cost C simultaneous increase in classification accuracy. We discuss this mismatch in Section II A 3, where we also consider alternative loss functions that more directly target the classification error. Once training is complete, the learned parameters α⋆ and ground state |ψT ,α⋆ ⟩ define the hypothesis used to predict labels for unseen inputs. An important feature of the Hamiltonian embedding framework is that, for particular classes of encoding unitaries U ( x, α), these predictions can be evaluated entirely classically. As discussed in Section III B, this is achieved by constructing a classical surrogate σ (ψT ,α⋆ ) of the learned quantum model from a finite number of measurements on the quantum device. The corresponding classical inference procedure is summarized in Algorithm 2. A.

Construction of the model

The Hamiltonian embedding framework provides substantial freedom in how the encoding Hamiltonian is constructed, allowing problem-specific choices to be incorporated. In the following, we outline practical guidelines for selecting the label-dependent observable O(y), the encoding unitary U ( x, α), the choice of the loss function and the ground-state optimization method. 1.

Choice of the observable

The label-dependent observable O(y) is measured at inference time to determine which label y ∈ {0, 1} should be assigned to a given input x. Recall that the energy loss in Eq. (7) is minimized when, for every training point ( xℓ , yℓ ) with yℓ = f ( xℓ ), the learned ground state |ψT ,α⋆ ⟩ satisfies Tr[O(yℓ )U ( xℓ , α⋆ )|ψT ,α⋆ ⟩⟨ψT ,α⋆ |U † ( xℓ , α⋆ )] = 1.

(8)

6 Algorithm 1 Supervised Training via Hamiltonian Embedding Input:

N , x ∈ Rd and y ∈ {0, 1}. 1. Training set T = {( xℓ , yℓ )}ℓ= ℓ 1 ℓ 2. Encoding unitary U ( x, α). 3. Desired training energy ε ≥ 0, learning rate η > 0, maximum number of iterations Tmax . Output: Trained parameters α⋆ and final state |ψT ,α⋆ ⟩.

1: O(y) ← 21 [ I + (−1)y A]

▷ Construct the label observable ▷ Initialize the variational parameters

2: α0 ← Unif([−1, 1] M ) 3: t ← 0 4: while t < Tmax do 5: for ℓ = 1 to N do 6: zℓ ← ( xℓ , yℓ ) 7: Hzℓ (αt ) ← U † ( xℓ , αt )O(yℓ )U ( xℓ , αt ) 8: Store a classical representation of Hzℓ (αt ) 9: end for N 10: HT (αt ) ← I − N1 ∑ℓ= 1 Hzℓ ( αt ) 11: |ψT ,αt ⟩ ≈ arg min ⟨ψ| HT (αt ) |ψ⟩ 12: 13: 14: 15: 16: 17: 18: 19: 20: 21:

▷ Embed the training point ▷ Construct the current embedding Hamiltonian ▷ Find an approximate ground state

|ψ⟩

▷ Evaluate the current training energy

Et ← ⟨ψT ,αt | HT (αt ) |ψT ,αt ⟩ if Et ≤ ε then break end if for ℓ = 1 to N do mαt ,ψT ,αt ( xℓ , yℓ ) ← ⟨ψT ,αt | Hzℓ (αt ) |ψT ,αt ⟩ end for bT (ψT ,α , αt ) according to Eq. (6) Evaluate C t for j = 1 to M do ∂ b C (ψT ,αt , α) g j (t) ← ∂α j T

▷ Evaluate the correct label score

▷

α=αt

22: end for 23: g ( t ) ← g1 ( t ) , . . . , g M ( t ) 24: αt+1 ← αt − ηg(t) 25: t ← t+1 26: end while 27: α⋆ ← αt N 1 ⋆ 28: HT (α⋆ ) ← I − N ∑ℓ= 1 Hzℓ ( α ) 29: |ψT ,α⋆ ⟩ ≈ arg min ⟨ψ| HT (α⋆ ) |ψ⟩

Estimate the gradient via parameter shift, keeping the state fixed

▷ Update the variational parameters

|ψ⟩

30: return (α⋆ , |ψT ,α⋆ ⟩)

Algorithm 2 Classical inference via Hamiltonian embedding Input:

1. New input point x ′ . 2. Trained parameters α⋆ . 3. Classical representation σ (ψT ,α⋆ ) of the learned state |ψT ,α⋆ ⟩. 4. Label-dependent observable O(y) for y ∈ {0, 1}. Output: Predicted label ỹ ∈ {0, 1}. 1: for y ∈ {0, 1} do

e ( x ′ , y) ← U † ( x ′ , α⋆ )O(y)U ( x ′ , α⋆ ) O e ( x ′ , y) ← ∑i ci ( x ′ , y, α⋆ ) Pi 3: O 4: py ( x ′ ) ← ∑i ci ( x ′ , y, α⋆ ) Tr[ Pi σ (ψT ,α⋆ )] 5: end for 6: ỹ ← arg maxy∈{0,1} py ( x ′ ) 7: return ỹ 2:

▷ Construct observable in Eq. (39) for x ′ associated with label y ▷ Decompose O′ into local Pauli operators ▷ Estimate the label score using the classical representation σ(ψT ,α⋆ ) ▷ Output the label with largest score

7 This occurs when, for every training example xℓ , the encoded state U ( xℓ , α⋆ ) |ψT ,α⋆ ⟩ lies in the eigenvalue-1 eigenspace of the observable O(y) associated with the correct label y ∈ {0, 1}. We note that Eq. (8) provides a sufficient, but not necessary, condition for correct classification. Since prediction in Eq. (5) is performed by selecting the label associated with the largest expectation value, it is enough that the score assigned to the correct label exceeds that of the other label. To give these observables a direct probabilistic interpretation, we can impose the natural condition O(0) + O(1) = I, so that their expectation values in Eq. (4) define a normalized probability distribution over the two classes. This condition is equivalent to parameterizing the two label observables in terms of a single Hermitian operator A as O(y) =

I + (−1)y A , 2

(9)

where y ∈ {0, 1} labels the two classes. We take A to be normalized such that its spectrum lies in [−1, 1], with extremal eigenvalues λmin ( A) = −1 and λmax ( A) = 1. We note that an operator A with this normalization can always be obtained whenever the extremal eigenvalues λmin and λmax of the original Hermitian operator are known, 2A−(λmax +λmin ) I by applying the affine transformation A 7→ . The choice of A plays a central role in the construction of λmax −λmin the model, as it determines how the two classes are represented in Hilbert space and can therefore be tailored to the learning task under consideration. We note however that the condition O(0) + O(1) = I is not necessary to define the model as labels can be assigned by comparing the expectation values of the two operators, without interpreting these values as probabilities. As we show in Appendix F and Appendix G, the choice of A, and O(y) in general, can influence properties of the learned ground state, such as its entanglement structure, and can thereby impact the classical simulability of the resulting model. From a practical point of view, the choice of the operator A and the corresponding observable O(y) affects the form of the final embedding Hamiltonian. Since conjugation by a constant-depth unitary U ( x, α) increases the support of a local observable only by a constant factor, choosing A to be local ensures that the resulting embedding Hamiltonian remains local and with polynomially many terms. In our numerical experiments, we choose the operator A to be a simple diagonal Hamiltonian of the form A=

1 n Zi . n i∑ =1

(10)

With this convention, the model perfectly classifies the training set whenever, for each training example, the learned ground state |ψT ,α⋆ ⟩ is mapped by the data-encoding unitary to a state in the eigenspace of O(y) associated with the correct label. In the particular case of the choice in Eq. (10), a zero energy cost is obtained when, for each training point ( xℓ , yℓ ) ∈ T ( |0n ⟩ yℓ = 0, ⋆ (11) U ( xℓ , α ) |ψT ,α⋆ ⟩ = |1n ⟩ yℓ = 1. As discussed above, the choice of A can influence the properties of the ground states learned by the model. For example, in Appendices F and G, we present alternative choices of O(y) that lead to highly entangled learned ground states. 2.

Choice of the unitary embedding

As discussed above, the embedded Hamiltonian HT (α) is determined by the choice of the parameterized encoding unitary U ( x, α). This choice influences the expressiveness of the model, and different encoding strategies may be better suited to different learning tasks. Since training requires the ground state of the embedding Hamiltonian to be prepared at each optimization step, the corresponding quantum ground-state preparation procedure should be efficient. A necessary condition is that the embedding Hamiltonian admit an efficient classical description suitable as input to the quantum ground-state preparation algorithm. For example, when O(y) can be expressed as a sum of local Pauli operators, choosing U ( x, α) to have constant depth ensures that each data-dependent term H( x,y) (α) = U † ( x, α)O(y)U ( x, α) remains of bounded locality O(1). In Appendix D, we show that the model achieves universality even with constant-depth encodings U and local O(y), using a finite number of qubits whose scaling however need not be efficient. We stress, however, that a constant-depth encoding unitary is not strictly required by our framework. For instance, ground-state solvers based on quantum phase estimation only require the ability to efficiently implement controlled time evolution under the target Hamiltonian. Therefore, the encoding unitaries

8 U ( x, α) may have depth growing with the system size, provided that they admit an efficient circuit description from which the corresponding Hamiltonian evolution can be implemented efficiently. In particular, this description need not be restricted to a decomposition in terms of Pauli operators. Another important consideration in the choice of the parameterized unitary U ( x, α) is that the cost function should admit efficient differentiation through the parameter-shift rule, as will be discussed in Section II B. For this to hold, the trainable parameters should enter the circuit through gates generated by operators with a finite, typically twopoint, spectrum, such as Pauli rotations of the form e−iθP/2 with P2 = I. Keeping all of this in consideration, in the numerical experiments in this work we choose the following encoding. The unitary U ( x, α) is taken to be a constant-depth circuit composed of layers of bounded-locality Pauli rotations, whose angles depend on both the input x ∈ Rd and the trainable parameters α ∈ R M . Concretely, we consider circuits of the form: K

i U ( x, α) = ∏ ∏ exp − 2 l =1 g =1

αlg,0 +

∑

r =1

αlg,r xr

Pgl

# ,

(12)

where K = O(1) denotes the number of layers, G = O(n) the number of parameterized gates per layer, and Pgℓ are bounded-locality Pauli operators. The products are taken according to the prescribed ordering of the gates in the circuit. In the numerical experiments, we instantiate Eq. (12) with K = 3 with each layer applying, on every qubit j, an independently parameterized Zj − Yj − Zj sequence of single-qubit rotations, followed by a nearest-neighbor Zj Zj+1 rotation. This architecture satisfies the locality requirement by construction: conjugating a one-local observable through a constant-depth local circuit only expands its support to a constant-size light cone. At the same time, because every trainable gate is a Pauli rotation, derivatives of the cost can be evaluated using the parameter-shift rule. Finally, we note that the choice of the encoding unitary U ( x, α), together with the observables O(y), determines whether inference on new input points can be performed entirely classically. This is the case, for instance, when U ( x, α) has constant depth and O(y) is local. Efficient classical inference is not restricted to this setting, and certain classes of deeper encoding unitaries also admit efficient classical surrogates. We discuss these different regimes in Section III B. Furthermore, in the classically hard task discussed in Section III C, neither the embedded Hamiltonian e (y) is local. Nonetheless, the construction still admits an efficient quantum procedure nor the dressed observable O for preparing the relevant ground state, as well as a problem-specific classical procedure for representing and classically evaluating the corresponding hypothesis. 3.

Choice of the loss function

The Hamiltonian embedding framework naturally defines the energy loss introduced in Eq. (7). For fixed training parameters α and state ψT ,α , its empirical average over the training set is precisely the energy of the corresponding embedding Hamiltonian HT (α). This choice is particularly convenient, since minimizing the loss is directly connected to finding low-energy states of the embedded Hamiltonian HT (α). Moreover, reaching zero energy loss on a training example guarantees that the corresponding point is classified correctly. However, the energy loss is not identical to the classification error and, in particular, a decrease in energy is not necessarily accompanied by a simultaneous increase in classification accuracy. To see this, let us consider the prediction rule in Eq. (5). In order for a point x to be correctly classified, it is sufficient that the model finds a state ψT ,α⋆ such that PψT ,α⋆ ,α⋆ ( x, y) > 0.5 if y is the correct label. This happens whenever the margin mα⋆ ,ψT ,α⋆ ( x, y) = ⟨ψT ,α⋆ |U † ( x, α⋆ )O(y)U ( x, α⋆ )|ψT ,α⋆ ⟩ is larger than 0.5. The training classification error is then defined as N b T (ψ, α) = 1 ∑ 1 mα,ψ ( xℓ , yℓ ) ≤ 0.5 , R N ℓ=1

(13)

which only depends on whether each example lies on the correct side of the decision boundary. We recall that, at each optimization step t, the state ψT ,αt is held fixed while the encoding parameters α are updated. For a fixed state ψ, minimizing the energy loss ℓenergy in Eq. (7) favors parameters α that increase the margins mα,ψ ( x, y) of the training examples ( x, y) ∈ T . However, since the energy loss ℓenergy is linear in mα,ψ ( x, y), it continues to reward increases in the margin even after a sample has already been classified correctly, i.e. when mα,ψ ( x, y) > 0.5. Consequently, the model may continue to lower the energy by increasing the confidence assigned to points

9 that are already correctly classified, without this necessarily leading to a corresponding reduction in the empirical classification error of Eq. (13). Directly minimizing Eq. (13) is not well suited to gradient-based optimization, since the error R̂T changes only when the score mα,ψ ( xℓ , yℓ ) of a training example crosses the classification threshold, and therefore provides no useful gradient for most parameter values. As in classical supervised learning, one can instead optimize a continuous surrogate of the 0–1 loss that more directly reflects the classification objective. We introduce two alternative losses as possible extensions aimed more directly at reducing classification error, although the resulting parameter updates need not decrease the Hamiltonian energy. For example, let us define a centered margin m̃α,ψ ( x, y) = 2mα,ψ ( x, y) − 1 so that m̃α,ψ ∈ [−1, 1]. We then consider the logistic loss ℓlog (m̃) = log 1 + e− βm̃ , (14) where β > 0 controls the sharpness of the loss around the decision boundary, and the hinge loss (γ)

ℓhinge (m̃) = max{0, γ − m̃},

(15)

where γ > 0 specifies a desired classification margin. In contrast to the energy loss, the hinge loss assigns no additional penalty to examples whose margin already exceeds γ, thereby focusing the optimization on misclassified points and on correctly classified points lying close to the decision boundary. Similarly, the logistic loss progressively suppresses the contribution of examples with large positive margins while assigning greater weight to misclassified examples and to points close to the decision boundary. More generally, given a surrogate loss ℓcls , we define the corresponding empirical classification-oriented cost as N bcls (ψ, α) = 1 ∑ ℓcls m̃ψ,α ( xℓ , yℓ ) . C T N ℓ=1

(16)

Importantly, introducing this alternative cost does not modify the embedding Hamiltonian or the ground-state preparation step. At optimization iteration t, the state is still obtained as an approximate ground state of the embedding Hamiltonian,

|ψT ,αt ⟩ ≈ arg min ⟨ψ| HT (αt ) |ψ⟩ .

(17)

|ψ⟩

The resulting state is then kept fixed while the encoding parameters are updated using the gradient of the chosen loss, bcls (ψT ,α , α) α t +1 = α t − η ∇ α C t T

(18)

α=αt

We note that the Hamiltonian embedding algorithm in Algorithm 1 is recovered by choosing ℓcls = ℓenergy . Furthermore, one may also interpolate between the two objectives through b(λ) = (1 − λ)C bT + λC bcls , C T T

λ ∈ [0, 1],

(19)

which allows the parameter optimization to balance the search for low-energy Hamiltonians with a loss that more directly targets classification performance. The choice of the outer-loop loss therefore provides an additional degree of freedom in the Hamiltonian embedding framework: the embedding Hamiltonian determines the state used by the model, while the loss function determines how the parameters are optimized within the family of embedding Hamiltonians. Unless stated otherwise, the numerical results presented in Section IV are obtained using the canonical energy-based objective in Eq. (7), which we find to be sufficient to achieve good classification performance on the tasks considered. We also note that replacing the energy cost in Eq. (6) with either the hinge or logistic cost does not increase the complexity of gradient estimation. Indeed, for a generic differentiable surrogate loss ℓcls (m), the derivative with respect to a variational parameter α j can be written using the chain rule as ∂m e α,ψ ( x, y) ∂ ′ e α,ψ ( x, y) = ℓcls e α,ψ ( x, y) ℓcls m . m ∂α j ∂α j

(20)

10 e α,ψ ( x, y) = 2⟨ψ|U † ( x, α)O(y)U ( x, α)|ψ⟩ − 1, we have Since m

e α,ψ ( x,y) ∂m = 2 ∂α∂ ⟨ψ|U † ( x, α)O(y)U ( x, α)|ψ⟩. ∂α j j

The first factor is a simple classical function of the margin, while the second corresponds to the derivative of the expectation value mα,ψ ( x, y) = ⟨ψ|U † ( x, α)O(y)U ( x, α)|ψ⟩. As discussed in Section II B, this derivative can be efficiently estimated through the parameter-shift rule. Therefore, the same quantum measurements used to estimate gradients of the energy cost can also be used to evaluate the gradients of the alternative surrogate losses, followed by a simple classical post-processing step. 4.

Ground-state approximation

At the heart of the Hamiltonian embedding approach is the use of a quantum computer to prepare the ground state of the embedding Hamiltonian HT (α) for varying parameter values α. Several ground-state preparation methods can be employed, and the most suitable choice depends on the structure of HT (α) and on the capabilities of the available quantum device. In the general setting of local Hamiltonians, a standard approach to ground-state preparation is quantum phase estimation (QPE) [CEMM98, AL99]. Given an efficient implementation of the time evolution generated by HT (α) and an initial state with non-negligible overlap with the ground state, QPE resolves the energy eigenvalues and projects the state onto the corresponding eigenspaces, producing the ground state with probability proportional to the initial overlap. When the Hamiltonian is gapped and the overlap is sufficiently large, QPE provides a systematic procedure for preparing the ground state with controllable precision. The Sample-based Krylov Quantum Diagonalization (SKQD) algorithm offers an alternative that avoids phase estimation and also has rigorous guarantees in certain regimes [YMI+ 25]. In the SKQD algorithm, time-evolution circuits are applied at multiple time steps and measured in a chosen basis, such as the computational basis, producing a set of sampled bitstrings that defines a reduced subspace onto which the Hamiltonian is projected and classically diagonalized. This procedure is compatible with current quantum hardware and is robust against device noise, provided the sampled subspace retains sufficient support on the low-energy sector. Both methods require an initial reference state with non-negligible overlap with the true ground state, ideally at least inverse-polynomial in the system size. In practice, the applicability of these ground-state solvers therefore depends on the embedding Hamiltonians remaining in a region of parameter space where reference states with sufficiently large ground-state overlap can be efficiently obtained. For structured Hamiltonian families, such reference states may be constructed by mean-field or tensor-network methods such as the density matrix renormalization group (DMRG) [Whi92, Sch05]. Identifying families of embedding Hamiltonians that simultaneously admit efficient quantum ground-state preparation while remaining difficult to simulate classically remains an important open question. Finally, the ground state of the data Hamiltonian could also be approximated through adiabatic preparation. Assuming that a ground state of HT (α0 ) is available, one could obtain a ground state of HT (α1 ) entirely quantumly by following an adiabatic path in parameter space. Efficient preparation requires, in particular, that the spectral gap along this path remain bounded below by an inverse polynomial in the system size. Exploring choices of the embedding U and observables O for which such gapped paths can be guaranteed is an interesting direction for future work.

Efficient gradient estimation and absence of barren plateaus

One of the main bottlenecks in training variational quantum ansätze, for example parameterized quantum circuits, is the presence of barren plateaus in the optimization landscape [MBS+ 18, QWZ+ 23]. In this regime, gradients become exponentially small with the system size. Since gradients are estimated from measurement outcomes on a quantum computer, identifying the correct descent direction may then require exponentially many measurement shots. In the Hamiltonian embedding model, gradients can be efficiently evaluated through the parameter-shift rule. Furthermore, for a class of Pauli-rotation encodings satisfying the conditions stated below, we can ensure that a relevant subset of trainable parameters avoids barren plateaus. We also note that, in contrast to standard PQC settings, the use of a shallow encoding circuit does not by itself imply efficient classical simulability of the framework, as finding the ground state of the corresponding local data Hamiltonian may remain hard. The gradient expression can be derived directly from the cost function in Eq. (6). As described in Algorithm 1, the optimization proceeds through a nested procedure. Given an initial parameter value α0 , the algorithm first computes the ground state of the encoding Hamiltonian HT (α0 ), outputting ψT ,α0 . The state is then kept fixed while b(ψT ,α , α). In particular, the algorithm searches for a the parameters are updated by minimizing the cost function C 0 new parameter value α1 for which the energy of HT (α1 ) evaluated on ψT ,α0 is minimized. As derived in Appendix B,

11 b(ψT ,α , α) in Eq. (6) with respect to parameter α j . Recalling the this update is guided by the gradient of the cost C 0 structure of HT (α) from Eq. (1), we obtain: ∂Hzℓ (α) 1 N ∂ b C (ψT ,α0 , α) = − ∑ ⟨ψT ,α0 | |ψT ,α0 ⟩ , ∂α j N ℓ=1 ∂α j

(21)

where N = |T | is the size of the full training set. In the encoding unitary U ( x, α) in Eq. (12), each parameterized gate has a rotation angle that depends linearly on the input features. For notational simplicity, we omit the constant bias term in the following derivation. The bias parameters are treated identically by regarding them as coefficients of an additional constant feature x0 = 1. For an input xi = ( xi1 , . . . , xid ), the angle of the j-th parameterized gate is d

θ j ( xi , α) = ∑ α j,k xik

(22)

k =1

where α j,k is the trainable coefficient associated with the k-th feature in the j-th gate. We can therefore decompose the full encoding unitary around this gate as U ( x i , α ) = U L j ( x i , α ) U j θ j ( x i , α ) UR j ( x i , α ), (23) where UR j ( xi , α) contains the gates applied before the j-th parameterized gate, while UL j ( xi , α) contains those applied after it. Since Uj is a Pauli rotation gate, its dependence on θ j , and consequently on each parameter α j,k , can be differentiated using the standard parameter-shift rule. Since differentiating the gate angle with respect to α j,k contributes a factor xik , and assuming that Uj is a Pauli rotation gate, the parameter-shift rule gives b ∂C 1 N k h b( j,+) b( j,−) i =− xi Ci − Ci , ∂α j,k 2N i∑ =1

(24)

where b( j,+) = ⟨ψT ,α | U † U † C Rj j 0 i b( j,−) = ⟨ψT ,α | U † U † C Rj j 0 i

∑

π α j,k xik +

∑

π α j,k xik −

2 2

UL† j O(yi )UL j Uj UL† j O(yi )UL j Uj

∑

π α j,k xik +

∑

π α j,k xik −

2 2

UR j |ψT ,α0 ⟩ ,

(25)

UR j |ψT ,α0 ⟩ .

(26)

where O(yi ) is the same labeled local observable as in Eq. (9). After the parameter-update step, the algorithm returns a new parameter value αt+1 . It then constructs the updated Hamiltonian HT (αt+1 ) and computes a new bT (ψT ,α , αt+1 ), and the optimization loop ground state ψT ,αt+1 . This state is used to evaluate the updated cost C t +1 continues. A key point is that, during each parameter update, the current ground state is kept fixed. Thus, when computing the gradient with respect to α, we do not differentiate through the dependence of ψT ,α on α. At iteration t, the algorithm first finds ψT ,αt for the Hamiltonian HT (αt ) and then updates the parameters by minimizing bT (ψT ,α , α). This update rule is variationally well motivated. If the new parameters satthe frozen-state cost C t b bT (ψT ,α , αt ), then, after recomputing the ground state ψT ,α , the energy cannot increase, isfy CT (ψT ,αt , αt+1 ) ≤ C t t +1 bT (ψT ,α , αt+1 ) ≤ C bT (ψT ,α , αt ). Moreover, although the ground state ψT ,α does depend implicitly on α, the i.e., C t t +1 Hellmann-Feynman theorem implies that, at each iteration, the derivative of the ground-state energy at the current parameter value can be evaluated while keeping the corresponding ground state fixed, as in Eq. (21). Further details are provided in Appendix B, where we also analyze the case of approximate ground states. As in the inference step, the gradient can also be estimated classically once an approximate ground state has been obtained from the quantum device. In particular, a classical-shadow representation of ψT ,αt can be used to estimate the expectation values of the two shifted observables ! ! π π ± † † k † k Oi,j = UR j Uj ∑ α j,k xi ± UL j O(yi )UL j Uj ∑ α j,k xi ± UR j . (27) 2 2 k k For constant depth U ( x, α) and local O(y), even the dressed observables remain local, or at most low-weight, so their expectation values can be efficiently inferred from the local Pauli classical shadows. Classical shadows can

12 estimate expectation values of local observables to inverse-polynomial precision with only polynomially many measurements. In order to train the model efficiently, the cost-function gradients must remain sufficiently large to provide a useful optimization signal, and in particular should not vanish exponentially with the system size as in barren-plateau regimes. For a general class of shallow encodings, we show in Appendix C that, once the energy cost is strictly below 1/2 by at least an inverse-polynomial amount, a situation that we observe in our numerical experiments and that corresponds to a non-trivial energy advantage, the gradient variance of the frozen-state cost admits an inversepolynomial lower bound for at least an inverse-polynomial fraction of parameters. In particular, in case of shallow encodings with linearly many parameters per layer, a constant energy advantage is sufficient to have a lower bound for a constant fraction of trainable parameters. This guarantees a useful training direction whenever the assumptions of the theorem are satisfied. Theorem 1 (Absence of barren plateaus (informal)). Consider a Hamiltonian embedding model defined on n qubits with an observable O(y) as in Eq. (9) with A = n−1 ∑in=1 Zi , M trainable parameters α ∈ R M , and an encoding U ( x, α) of depth K whose final layer applies one independently trainable X- or Y-rotation to each qubit, with rotation angles independent of the input x. At iteration t, keep the prepared state |ψT ,αt ⟩ fixed and define bT (ψT ,α , α) = ⟨ψT ,α | HT (α)|ψT ,α ⟩. C t t t bT (ψT ,α , αt ) ≤ 1 − γ, for some γ > 0, then at least 2γ2 n final-layer parameters satisfy If C t 2 " # bT (ψT ,α , α) ∂C γ2 t Varα j ≥ 2. ∂α j 4n

(28)

(29)

Here, the variance is taken over the corresponding rotation angle uniformly distributed in [0, 2π ], while the dataset, prepared state ψT ,αt , and all other parameters remain fixed at their iteration-t values. Consequently, for encodings U of K = Θ(1) layer, M = Kn parameters and γ = Θ(1), a constant fraction of trainable parameters has inverse-polynomial gradient variance. The precise assumptions on the encoding U and the complete proof of the theorem are given in Appendix C. III.

PROPERTIES OF THE HAMILTONIAN EMBEDDING MODEL

The Hamiltonian embedding model exhibits several desirable properties for a quantum machine learning model. First, we show that once trained, the model generalizes well to unseen inputs. Moreover, for certain families of encoding unitaries and label observables, the model can be evaluated entirely classically once the optimal ground state has been obtained during training. Finally, we identify a machine learning task believed to be classically hard, yet efficiently solvable on a quantum device using the Hamiltonian embedding framework. A.

Generalization bounds

Let |ψT ,α⋆ ⟩ be the ground state found by the algorithm after the training phase. The label for a new point x ′ is assigned by the hypothesis function feT ,α⋆ in Eq. (5). We define the total classification risk as the probability the model misclassifies a point x sampled from the input distribution D over the whole input space Ω ⊆ Rd : R(ψ, α) = Pr [ feT ,α ( x ) ̸= f ( x )] x ∼D

(30)

where we recall that f is the target labeling function and feT ,α is the model implemented function for the Hamiltonian parameters α ∈ Θ ⊆ R M and the corresponding reached ground state ψT ,α . We show that the total classification risk R(ψ, α) can be controlled through the expected energy cost over the full input distribution, computed as C (ψ, α) = ⟨ψ| HD (α) |ψ⟩, where HD (α) := Ez∼D [ I − Hz (α)] is the population Hamiltonian at parameter α. Specifically, HD represents the average encoding Hamiltonian associated with a random labeled input drawn from the full data distribution. In order to connect the total risk to the training procedure, we first bound the difference between the total cost energy and the empirical cost in Eq. (6): b(ψT ,α⋆ , α⋆ ) = ⟨ψT ,α⋆ | HD (α⋆ ) − HT (α⋆ ) |ψT ,α⋆ ⟩ . C (ψT ,α⋆ , α⋆ ) − C

(31)

13 b and This difference is the energy generalization gap. The training algorithm minimizes the empirical energy cost C ⋆ obtains Hamiltonian parameters α and an associated ground state |ψT ,α⋆ ⟩. We then show that a small empirical cost together with a small generalization gap implies a small classification error on unseen input points. In particular, we find that the probability of observing a large deviation ∥ HD (α⋆ ) − HT (α⋆ )∥∞ can be bounded with a dependence that scales linearly with the number of qubits n and the number of trainable parameters M, up to logarithmic factors. The proof relies on a matrix Hoeffding inequality, which controls the energy deviation for a fixed parameter value, together with a covering-number argument, which extends the bound uniformly over the full parameter space. For the purpose of deriving a uniform generalization bound, we restrict the parameter search to a compact domain Θ = [−1, 1] M . This restriction is specific to the theoretical result below and is not required by the Hamiltonian embedding framework itself. In particular, we obtain the following theorem. Theorem 2 (Uniform energy-based generalization bound). Assume a Hamiltonian embedding model such that 0 ⪯ Hz (α) ⪯ I and that there exists a constant L > 0 such that ∥ Hz (α) − Hz ( β)∥∞ ≤ L∥α − β∥∞ for all points z = ( x, y) ∈ Ω × {0, 1} and all α, β ∈ [−1, 1] M . Then, for every ε ∈ (0, 4L], with probability at least 1 − δ over the draw of N-size training datasets, s sup α∈[−1,1] M

∥ HT (α) − HD (α)∥∞ ≤

2D 24L ε 8 log + M log + . N δ ε 4

(32)

where D = 2n is the Hilbert space dimension. The derived energy generalization gap bound can be used to control the population risk of the trained model. In b(ψT ,α⋆ , α⋆ ). Then, the number of particular, suppose that the training procedure reaches an empirical energy cost C training samples required to ensure that the population risk is upper bounded by twice the empirical energy cost, up to a constant additive error, scales linearly with the number of qubits and the number of trainable parameters, up to logarithmic factors. This is stated in the following corollary of Theorem 2. Corollary 3 (Bound on the total risk). Assume the setting of Theorem 2. Let (ψT ,α⋆ , α⋆ ) be the output of the trainb(ψT ,α⋆ , α⋆ ) = ⟨ψT ,α⋆ | HT (α⋆ ) |ψT ,α⋆ ⟩ be the empirical energy reached by the trained model, and let ing procedure. Let C i h C (ψT ,α⋆ , α⋆ ) = ⟨ψT ,α⋆ | HD (α⋆ ) |ψT ,α⋆ ⟩ be its population energy. Moreover, let R(ψT ,α⋆ , α⋆ ) = Prx∼D f˜ψ ⋆ ,α⋆ ( x ) ̸= f ( x ) T ,α

be the classification risk of the trained model. Then, with probability at least 1 − δ over the draw of the training set T , the classification risk of the trained model can be bounded by b(ψT ,α⋆ , α⋆ ) + ε, R(ψT ,α⋆ , α⋆ ) ≤ 2C

(33)

provided that the number of training samples satisfies 2D 24L 128 + M log . N ≥ 2 log δ ε ε

(34)

For simplicity, we stated the bound for variational parameters α ∈ [−1, 1] M . We note that the same argument applies to any bounded parameter domain Θ = [− R, R] M , with the corresponding modification of the covering number. Finally, in Appendix A we note that many natural encoding unitaries U ( x, α), including the one we propose in Section II A, satisfy the Lipschitz assumption on the terms Hz (α) required by Theorem 2. Thus, a sample complexity linear in n and M, for fixed ϵ, δ, up to logarithmic factors is enough to guarantee a classification error less than ε in the case of perfect training, i.e. when α⋆ is such that HT (α⋆ ) has zero ground energy and the true ground state of HT (α⋆ ) is found during the training. We note that, if the ground states of the data Hamiltonian admit a polynomial-size circuit representation, generalization guarantees could also be obtained from existing QML bounds [CHC+ 22], whose sample complexity is controlled by the number of parameterized gates in the circuit. However, our bound does not depend on the circuit complexity required to prepare the ground state, but only on the number M of parameters entering the embeddingHamiltonian family and on the number of qubits n. In particular we notice that, even when the ground state is efficiently preparable, the number of parameterized gates required by such a preparation circuit may be substantially larger than M.

14 B.

Classical inference and surrogate models

An interesting feature of the Hamiltonian embedding model is that, for specific choices of the encoding unitaries U ( x, α), inference can be carried out entirely classically, thereby removing the need to access a quantum device at deployment. In these settings, the role of the quantum computer is confined to the training stage, where it is used to prepare a suitable advice state by minimizing the energy of the data-embedding Hamiltonian. Quantum models that allow for classical inference after quantum training were introduced in [JGM+ 24] under the name shadow models. In a shadow model, the quantum computer provides a quantum advice state obtained during training. We can then see the Hamiltonian embedding model as a special case of shadow model in which the advice state is obtained as the ground state of the data-dependent Hamiltonian. Once the trained ground state has been obtained, the structure of the encoding unitary U ( x, α⋆ ) can be exploited to evaluate the prediction rule in Eq. (4) without further access to the quantum device. In this section, we present different classes of encoding unitaries U ( x, α) for which inference can be performed efficiently classically, together with the corresponding construction of classical surrogate models. Finally, in Section III B 4 we also identify settings in which the learned ground state acts as a useful inductive bias, enabling quantum Hamiltonian embedding to solve the task with fewer training examples than a classically regularized surrogate. 1.

Fourier surrogates

Let us consider the case of a constant depth encoding unitary U ( x, α). A first example of a classical surrogate for evaluating the trained model is based on the Fourier decomposition of the final circuit appearing in Eq. (4). Indeed, parameterized quantum circuits are known to admit a representation as generalized trigonometric polynomials [SSM21], whose set of frequencies is determined solely by the data-encoding gates. Specifically, let x ∈ Rd be the input points, then any parameterized quantum circuit of the kind f α ( x ) = ⟨ψ| U † ( x, α) MU ( x, α) |ψ⟩

(35)

can be equivalently expressed in the form f α (x) =

∑ cω (α)eiω·x ,

(36)

ω∈Ωα

where Ωα denotes the set of accessible frequencies and cω (α) ∈ C are the corresponding Fourier coefficients. The set of accessible frequencies Ωα only depends on the gates U ( x, α) used to encode the inputs x. In our model, where U ( x, α) is shallow and the measured observables are local, the number of distinct frequencies in Ωα can be shown to scale polynomially with the input size. In the case where the frequency spectrum Ωα is polynomially large, recent results have shown that a classical representation of the circuit can be learned efficiently. In particular, we can use the results in [SEM23] to learn an approximation for the coefficients cω (α) in Eq. (36) and construct a classical surrogate of the form: gα ( x ) =

∑ c̃ω (α)eiω·x ,

(37)

ω∈Ωα

with the guarantee that

Pr sup | f α ( x ) − gα ( x )| ≤ ϵ ≥ 1 − δ, x

(38)

using polynomially many queries in 1/δ and 1/ϵ of the original quantum circuit. 2.

Pauli surrogates

For constant depth encoding U ( x, α), we can define a different classical surrogate model by exploiting the locality of the observables considered at inference time. In particular, rather than learning the Fourier coefficients of the final circuit, we construct a classical representation of the trained model by expanding the dressed observable in a basis of local Pauli operators and estimating their expectation values on the learned ground state using classical shadows. When A in Eq. (9) is a sum of local Pauli operators, these estimates provide a classical description sufficient

15 to evaluate the inference rule entirely through classical-shadow post-processing [HKP20]. In particular, once the optimized state |ψT ,α⋆ ⟩ is obtained during training, the probability PψT ,α⋆ ,α⋆ ( x ′ , y) in Eq. (4) for a new input x ′ can be estimated from a classical-shadow representation of the learned ground state. Indeed, let us consider the dressed observable constructed by conjugating the original observable in Eq. (9) with the trained encoding unitary for a new input point x ′ : Õ( x ′ , y) =

I + (−1)y Ã , 2

Ã( x ) = U † ( x ′ , α⋆ ) AU ( x ′ , α⋆ ).

(39)

Since we choose A, and therefore O(y), to be local, conjugating it by a constant-depth local encoding unitary e ( x′ , y) = U ( x, α) only spreads its support within a constant-size light cone. Therefore, the dressed observable O U † ( x ′ , α⋆ )O(y)U ( x ′ , α⋆ ) can then be efficiently decomposed into local Pauli operators. Indeed, if U is a constantdepth circuit composed of local gates, the backward light cone of each measured term contains at most k = O(1) qubits. Each such term can therefore be expanded into at most 4k Pauli operators. Since the observable contains O(n) local terms, we can write K

e( x ′ ) = ∑ ci ( x ′ , α⋆ ) Pi , A

(40)

i =1

with K = O(n4k ), where the Pi are Pauli observables of weight at most k. Moreover, the coefficients ci ( x ′ , α⋆ ) can be computed efficiently classically by simulating the constant-size backward light cones. Once the approximate ground state |ψT ,α⋆ ⟩ has been prepared quantumly, a classical shadow built from ! 3k maxx ∥c( x, α⋆ )∥21 K log (41) S=O δ ϵ2 measurements suffices to estimate PψT ,α⋆ ,α⋆ ( x ′ , y) to additive error ϵ, uniformly over all x ′ ∈ Ω, with probability at least 1 − δ [HKP20]. 3.

Deep Clifford encodings

A constant depth U ( x, α) is not the only case in which the Hamiltonian embedding model allows for an entirely classical inference after training. As an example, let us consider an encoding unitary of the form U ( x, α) = U1 ( x, α)VU2 (α), where U2 (α) is independent of the input x, V is a deep Clifford circuit, and U1 ( x, α) is a shallowdepth unitary. Similarly to the Pauli-surrogate case, we can define the dressed observable obtained by conjugating O(y) with the shallow part of the encoding, K

e 1 ( x, y) = U † ( x, α⋆ )O(y)U1 ( x, α⋆ ) = ∑ ci ( x, α⋆ ) Pi , O 1

(42)

i =1

where K = poly(n) and the Pi are local Pauli observables. Since V is a Clifford circuit, conjugation by V maps every Pauli observable to another Pauli observable, V † Pi V = Qi .

(43)

Although the resulting operators Qi may have support on a number of qubits growing with the system size, the number of terms in the decomposition remains polynomial. The prediction probability rule of Eq. (4) can therefore be written as K

PψT ,α⋆ ,α⋆ ( x, y) = ∑ ci ( x, α⋆ )mi , i =1

mi = ⟨ψT ,α⋆ | U2† (α⋆ ) Qi U2 (α⋆ ) |ψT ,α⋆ ⟩ .

(44)

Crucially, since U2 is independent of x, the coefficients mi need to be estimated only once after training. Since the operators Qi may have large Pauli weight, product-Pauli classical shadows may not be efficient. However, as only polynomially many such observables appear, their expectation values can be measured directly on the quantum device and stored classically. Predictions for new inputs can then be obtained entirely classically by computing the coefficients ci ( x, α⋆ ) and evaluating the above sum.

16 4.

Dequantizability of the model

The existence of efficient classical surrogates for the trained quantum model raises the natural question of whether the full learning procedure can itself be reproduced classically, potentially eliminating the need for quantum computation altogether. For example, let us consider the Pauli surrogate. Suppose that the dressed observable Õ admits a decomposition into polynomially many Pauli operators, as in Eq. (40). The prediction function of the trained quantum model can then be written as K

f˜α⋆ ( x ) = ∑ ci ( x, α⋆ )mi (α⋆ ),

(45)

where mi (α⋆ ) = ⟨ψT ,α⋆ | Pi |ψT ,α⋆ ⟩ are the expectation values of the Pauli operators Pi on the learned ground state, while ci ( x, α⋆ ) are the coefficients appearing in the Pauli decomposition of the dressed observable as in Eq (40). We can then define a fully classical surrogate of the form K

gα ( x ) = ∑ ci ( x, α)m̃i (α),

(46)

where the coefficients m̃i (α) ∈ R are treated as trainable classical parameters and are not required to correspond to expectation values of Pauli operators on any quantum state. In the quantum model of Eq. (45), the moments mi (α⋆ ) are not independent parameters, but must be jointly compatible with a single physical quantum state. These compatibility constraints can in principle act as a useful inductive bias, although they do not by themselves imply either easier optimization or improved generalization relative to classical alternatives. To investigate this question, in Appendix G we construct a learning problem for which the benefit of these physical compatibility constraints can be established analytically, yielding a provable separation in the number of training examples required by Hamiltonian embedding and a corresponding classically regularized surrogate. We complement this result in Appendix K with a finite-size numerical example exhibiting the same qualitative behavior. Importantly, both tasks are classically tractable, and the purpose of these examples is to isolate the generalization advantage of the physical compatibility constraints imposed by the learned ground state rather than to demonstrate a computational quantum advantage. A more constrained surrogate can be obtained when each dressed readout term Ã j = U ( x, α)† Pj U ( x, α) in Eq. (40) is supported on a fixed region L j of the n qubits, independent of x and α, with | L j | = O(log n). Writing the observable A in Eq. (9) as A = ∑Kj a j Pj , we define a classical surrogate as K gσ,α ( x ) = ∑ a j Tr σj Ã j ,

(47)

where each σj is a density matrix on the fixed subsystem L j , with no global compatibility constraints, and Ã j ( x, α) bT (σ, α) = has support on L j . At fixed α, from the energy-loss definition in Eq. (6), the empirical energy becomes C N 1 K 1 −1 ∑ℓ=1 sℓ Ã j ( xℓ , α) and sℓ = (−1)yℓ . Since the matrices σj can be opti2 − 2 ∑ j=1 a j Tr[ σj B j ( α )], where B j ( α ) = N mized independently, energy minimization reduces to K eigenproblems by choosing each σj as a projector onto a

top eigenvector of a j Bj (α). These matrices have dimension at most 2| L j | = poly(n), so the inner optimization runs in polynomial-time when the dressed terms can be computed efficiently and achieves an energy lower or equal to the ground-state energy of the original embedding Hamiltonian HT (α). Under ∑Kj=1 | a j | ≤ 1 and uniform operatornorm Lipschitz bounds on the individual dressed terms Ã( x, α), the same proof strategy of Theorem 2 extends to this relaxation guaranteeing also the same generalization bounds. However, we note that this surrogate construction may not remain efficient when the encoding unitary U ( x, α) is such that the supports of each dressed term Ã j ( x, α) depend on x. In this case, the support of their union over all training points may involve all n qubits, making the dimension 2| L j | of σj scale exponentially. For instance, the construction in Appendix G also exhibits such inputdependent supports, providing an example in which the efficiency guarantee for this surrogate no longer applies. Finally, we note that the polynomial-time training statements of the surrogates do not cover the outer optimization over α, which is not convex and could benefit from the inductive bias induced by the ground state of the embedding Hamiltonian. In particular, a separate question is whether the full Hamiltonian embedding learning procedure can itself be efficiently reproduced classically. The existence of efficient classical surrogates for inference does not, in general, imply that the corresponding model can also be learned efficiently classically. In particular, one natural route beyond

17 the shallow-surrogate regime is to use encodings whose dressed observables have growing light cones and do not admit a polynomial-sized, efficiently evaluable feature representation. Increasing the depth of U can therefore invalidate the particular Pauli- and Fourier-surrogate arguments used above. Depth alone, however, is not sufficient, since structured deep circuits may remain classically tractable and low-bond-dimension learned states may still permit efficient tensor-network evaluation. The relevant question is whether the task-dependent moments defining the classifier can be learned or computed classically, rather than circuit depth or entanglement alone. The construction in the following section provides an explicit setting in which the trained hypothesis is efficiently evaluable but learning it classically is nevertheless conjectured to be hard under cryptographic assumptions. C.

Encoding a classically hard problem in the Hamiltonian embedding model

When introducing a new quantum machine learning method, it is natural to ask whether it can be applied to problems that lie beyond the reach of classical algorithms. In the Hamiltonian embedding framework, the training procedure requires finding the ground state of an encoding Hamiltonian. Since ground-state preparation is, in general, hard even for quantum computers, it is not immediately clear whether this approach can provide a rigorous quantum speedup. In this section we present a learning task that is believed to be classically intractable, yet can be expressed within our Hamiltonian-embedding framework and solved using a quantum computer. The learning problem is the binary classification task described in [JGM+ 24], where input x ∈ Z N = {0, 1, ..., N − 1} are labeled in two classes y = {0, 1} by the following labeling function:  N−1 1, if √ 3 x mod N ∈ s, s + , 2 (48) y( x ) = f N,s ( x ) =  0, otherwise. where N = pq is a large n-bit integer product of the two primes p and q and the interval in Eq. (48) is interpreted cyclically in Z N . In particular, the labeling function partitions inputs into two classes according to whether their discrete cube root (DCR) modulo N lies in a specified half-interval, determined by s. Under the DCR assumption, computing discrete cube roots modulo large semiprimes N is believed to be classically intractable [KV94], and the associated learning problem is also believed to be classically hard [JGM+ 24]. Notably, the DCR function is a trapdoor function: there exists a key d ∈ Z N that allows one to evaluate the function efficiently. More precisely, for the √ DCR function, there exists a d such that 3 x mod N = x d mod N. Since modular exponentiation can be computed efficiently classically, evaluating the DCR function on different inputs x becomes easy once d is known. However, determining the key d requires factoring the integer N, which is believed to be computationally hard for classical algorithms for large integers N. The Hamiltonian embedding framework can reproduce the DCR labeling function in Eq. (48) as its hypothesis function. More precisely, one can construct an encoding Hamiltonian whose ground state induces, through the prediction rule in Eq. (E25), a hypothesis fe that well approximates the DCR labeling function. In particular, we consider the following encoding unitary U ( x ) defined over the Hilbert space H = span{|d, s⟩ : d, s ∈ Z N } ⊗ C2 and the observable O(y): U ( x ) |d, s⟩|b⟩ = |d, s⟩ |b ⊕ gd,s ( x )⟩

& O(y) = I ⊗ |y⟩⟨y|,

x, d, s ∈ Z N , y ∈ {0, 1}.

(49)

where gd,s is defined as: gd,s ( x ) =

 1, 

N−1 mod N ∈ s, s + , 2 otherwise

if x d

(50)

and it is classically computable as modular exponentiation can be computed efficiently. We then construct the embedding Hamiltonian HDCR as HDCR = ∑ H ( xℓ , yℓ ).

(51)

ℓ

We note that the encoding Hamiltonian HDCR admits an efficient classical description, since the encoding unitary U in Eq. (49) can be implemented by a polynomial-size classical reversible circuit, so that each Hamiltonian term associated with a training point ( xℓ , yℓ ) ∈ T can be efficiently specified as: H ( xℓ , yℓ ) = I − U ( xℓ )† O(yℓ )U ( xℓ ) + I ⊗ |1⟩⟨1|.

(52)

18 The ground space of HDCR contains states |d, s, 0⟩ associated with hypotheses gd,s that are consistent with the training data. As shown in Appendix E, a quantum algorithm can efficiently find a ground state |d⋆ , s′ , 0⟩, whose corresponding hypothesis gd⋆ ,s′ has small population error on the DCR learning task. In particular, following the construction in Section II, for a new point x the predicted label is the value of y that maximizes Tr O(y) U ( x ) |d⋆ , s′ , 0⟩⟨d⋆ , s′ , 0| U ( x )† . Therefore, we can prove the following theorem: Theorem 4 (Encoding of a hard problem (informal)). There exists a classification problem for which the Hamiltonian embedding framework yields a classically efficiently describable data-embedding Hamiltonian whose ground states induce, through the prediction rule in Eq. (5), hypotheses with arbitrarily small classification error. Moreover, a quantum computer can efficiently find a corresponding ground state and construct the associated hypothesis. In contrast, assuming the classical hardness of computing discrete cube roots, no efficient classical learner can achieve the same guarantee with non-negligible probability. The full specification of the result, together with the formal proof, can be found in Appendix E. Finally, we note a distinction between Hamiltonian embedding and some earlier examples of quantum models that efficiently solve learning tasks believed to be classically intractable [LAT21, AG26, MGD26, ZZN+ 26, GGPD25], concerning where the hard computation takes place. In these examples, the quantum models are explicitly constructed to incorporate the correct hard computation required to solve the specific task into their circuit architectures. In the discrete-logarithm-based constructions of [LAT21, GGPD25], for example, the encoding circuits incorporate discretelogarithm computations, while the remaining training optimization then becomes easy. In our construction, the hard computation instead takes place during data-dependent quantum optimization through ground-state preparation, while the encoding circuits themselves perform efficiently computable classical operations. Here ground-state preparation is treated as a subroutine. The encoding circuits therefore need not be designed a priori to perform a specific hard computation and, in fact, may also be classically tractable, as is the case in the classically hard learning task presented here. IV.

EXPERIMENTS

We next test whether the optimization in Algorithm 1 remains feasible at nontrivial system sizes and whether the resulting low-energy states achieve low classification error on a standard benchmark dataset. In particular, we evaluate the Hamiltonian embedding algorithm on the binary classification task of Fashion-MNIST dress and shirt images, with inputs reduced by PCA to n features for an n-qubit model. We first train the model using a classical iterative solver as a proxy [KPM+ 26] to approximate the ground state at each iteration. This benchmark uses 100 random seeds for each of the 25 even system sizes n = 2, 4, . . . , 50, with each seed determining both the data split and the initial model parameters. We compare test accuracy against classical machine-learning baselines and the fidelity and projected quantum kernels in [HBM+ 21]. We then implement the full hybrid training procedure, using sample-based Krylov quantum diagonalization on IBM Heron processors for ground-state approximation. The selection criteria for the hardware runs are given below and in Appendix J. A.

The model instance used in the experiments

All Hamiltonian embedding experiments use the same model architecture and training objective introduced in Section II, which we summarize here and discuss in detail in Appendix H. The encoding unitary is the constantdepth circuit in Eq. (12) with K = 3 layers. We choose K = 3 to explore nontrivial Hamiltonians while keeping the number of Pauli terms computationally tractable. Each layer applies three single-qubit rotations, generated by Z, Y, and Z, to every qubit, followed by a ring of nearest-neighbor two-qubit ZZ rotations. The circuit therefore contains 12n parameterized gates, whose trainable coefficients are α ∈ R12n×(d+1) . For an input x ∈ Rd , the angle of gate g is θ g ( x, α) = ∑rd=1 α g,r xr + α g,0 . This is equivalent to the encoding in Eq. (12) with a constant feature appended to the input. We use the one-local readout A = n1 ∑in=1 Zi from Section II A. The outer optimization minimizes the energy loss defined in Eq. (7), which is sufficient to achieve good classification performance in our experiments. Two features make the experiments feasible at these system sizes. First, the light-cone argument in Section II A shows that each data point contributes a linear combination of the same O(n) bounded-weight Pauli observables to the embedding Hamiltonian HT (α). We obtain these strings by propagating the readout observable through the circuit in the Heisenberg picture and combining identical terms. Second, as discussed in Section II B, the state is held fixed during each parameter update, so the gradient is given by the expectation value of ∂α HT (α) in the current state.

19 a

1.0

1.0 0.932

0.920 0.896

0.9

0.8

Best classical ML Classical iterative solver Projected quantum kernel (E1) Quantum kernel (E1) Projected quantum kernel (E2, w = 4)

0.7

test accuracy

0.813

0.8

0.698

0.7

0.6

0.536

0.5 0.5 0

qubits (PCA dimension, n)

SVC-rbf (best classical)

HE trained on QPU

Classical iterative solver

Projected quantum kernel (E1)

Projected quantum kernel (E2, w=4)

Quantum kernel

Figure 1. a) Test accuracy versus system size for various methods on a binary classification problem from the Fashion-MNIST dataset. The methods are: best-performing classical baseline (gray), classical iterative solver (orange) used in HE framework, quantum kernel based on a product-state embedding (blue), 1-reduced density matrix (1-RDM) projected quantum kernel with the same embedding [HBM+ 21] (green), and the 1-RDM projected quantum kernel based on the E2 embedding of Ref. [HBM+ 21] (pink). b) Held-out test accuracy of the five runs with the highest validation accuracy at n = 50. In this plot, we additionally show the accuracy obtained by the HE model trained on hardware (yellow).

We optimize the parameters using Adam [KB15] together with an ℓ2 penalty. At each training epoch, we evaluate the model on the validation set. After training, we keep the variational parameters and associated learned quantum state from the epoch with the highest validation accuracy, and evaluate the resulting model once on the held-out test set. We first study Hamiltonian embedding framework using a classical method for approximating ground states that we call “classical iterative solver” [KPM+ 26]. In these simulations, we approximate the ground state at each training iteration using an iterative solver based on repeated applications of HT (α), deterministically retaining the bitstrings with the highest probabilities, akin to a power iteration with truncation. We reserve the term “HE” for the hardware runs using the real-time SKQD protocol. For the hardware runs, we approximate the ground state at each training iteration using SKQD [YMI+ 25]. Starting from a reference state, we generate a sequence of Krylov states by Trotterized real-time evolution under the current training Hamiltonian and measure them in the computational basis. The sampled bitstrings span a reduced subspace onto which we project HT (α) and diagonalize it classically. The classical iterative solver and SKQD may select different basis states and therefore produce different reduced subspaces. Both methods nevertheless project the full embedding Hamiltonian onto the selected subspace, so the lowest projected energy remains a variational upper bound on the true ground-state energy. The classical solver provides a useful reference for the optimization and classification performance of Hamiltonian embedding, although its performance need not predict that of SKQD.

Results

We compare Hamiltonian embedding with the fidelity and projected quantum kernels studied in Ref. [HBM+ 21], as well as with classical machine-learning methods. We use the same image-classification benchmark studied in Ref. [HBM+ 21], which is a binary classification of Fashion-MNIST images [XRV17] into the classes dress and shirt. We first use principal component analysis (PCA) to reduce the input dimension to the number of qubits n. For each run, a pool of 800 samples is divided into 448 training, 112 validation, and 240 held-out test examples. We first split the data 70/30 into training and test sets, and then split the training set 80/20 into the final training and validation sets. At each n, all methods use the same random seeds and therefore the same data splits. Figure 1(a) reports the mean test accuracy across the 100 runs at each system size, with error bars showing one standard deviation. We note that the best classical ML baseline, which corresponds to the support vector classifier

20 with a radial basis function kernel, performs the best, reaching an accuracy of approximately 0.946. The classical iterative solver, which may be viewed as a proxy for the Hamiltonian embedding framework, achieves a higher mean test accuracy than the best reported quantum kernel baseline, with margins ranging from 0.018 to 0.067. The mean accuracy increases from 0.880 at n = 2 to approximately 0.924 at intermediate system sizes before gradually decreasing to 0.886 at n = 50. The fidelity kernel based on the product-state embedding E1 from [HBM+ 21] decreases from 0.798 to 0.515, whereas the projected 1-RDM kernel based on the same embedding decreases more gradually from a maximum of 0.887 to 0.831. We also study a projected quantum kernel constructed from the 1-RDMs of feature states generated using a simplified version of the E2 embedding from Refs. [HCT+ 19, HBM+ 21]. These feature states are prepared using circuits inspired by the instantaneous quantum polynomial-time (IQP) model, which we call IQP-like circuits. To make their classical simulation tractable, we restrict the interactions to bands of neighboring qubits of width w = 4, allowing the circuits and resulting kernels to be evaluated using matrix product states. The test accuracy of this kernel decreases rapidly with n and levels off near 0.65, above that of the fidelity kernel at large n. We next turn to the hardware results. In these runs, we prepare the SKQD Krylov states described above on IBM Heron processors and measure them in the computational basis. We then project HT (αt ) in the sampled bitstring subspace and diagonalize it classically. We evaluate the gradients and perform the Adam updates classically, following Algorithm 1. We select five hardware runs based on the results obtained using the classical iterative solver to select data splits and initial model parameters α that lead to high validation accuracy. As shown in Fig. 1(b), Hamiltonian embedding is competitive with the best classical baseline and outperforms the other quantum methods considered. All n = 50 runs were performed on ibm pittsburgh, a 156-qubit IBM Heron r2 processor. The 50 data qubits were mapped to a linear chain chosen greedily to minimize two-qubit gate errors. Moreover, we apply dynamical decoupling and measurement twirling but no additional error mitigation. These experiments demonstrate that Hamiltonian embedding models can be studied at scale on current quantum processors.

DISCUSSION

In this work, we introduced Hamiltonian embedding, a QML framework for classification in which the entire training set is encoded in a single parameterized Hamiltonian. Training minimizes the ground-state energy or a more general loss function over the encoding parameters, with the loss function estimated by a quantum algorithm at each iteration. The ground state of the optimized Hamiltonian then defines a classifier for unseen inputs. For shallow encoding circuits, measurements of the prepared states yield classical representations that allow the loss and its gradients to be evaluated efficiently. Moreover, the representation of the final state supports inference without further quantum computation. Therefore, the role of quantum processors is restricted to approximating ground states during training. We further showed, analytically and numerically, that the physical compatibility constraints imposed by the learned state can provide a useful inductive bias for certain tasks. Finally, under standard cryptographic assumptions, Hamiltonian embedding can encode certain learning tasks that are intractable for efficient classical algorithms. For generic learning tasks, however, it remains unclear when the inductive bias of Hamiltonian embedding offers an advantage over alternative classical or quantum methods. This question is particularly difficult for classical datasets, where symmetry-based arguments may not apply directly [GGC+ 24, RBN+ 22], and will likely require heuristic model design together with systematic, task-specific benchmarking. Practical use also requires the model to remain trainable and to generalize as the system size grows. For fixed generalization error, our analytical bounds show that the required number of training samples scales linearly with the number of qubits and trainable parameters, up to logarithmic factors. Our numerical results provide complementary evidence for trainability. On binary Fashion-MNIST, Hamiltonian embedding can be trained at scale, including implementations on IBM quantum processors, while maintaining competitive test performance as the system size increases. In the same benchmarks, the generalization performance of the quantum-kernel method of Ref. [HBM+ 21] deteriorates with system size, whereas Hamiltonian embedding also achieves higher accuracy than the projected quantum kernels designed to mitigate this behavior. Hamiltonian embedding therefore offers a promising alternative to existing QML models. Several aspects of the framework remain unexplored. In particular, the freedom to choose both the encoding and the readout observable may enable constructions with improved trainability or expressivity. A central open question is which learning problems are naturally suited to Hamiltonian embedding and whether it can offer advantages over classical methods beyond cryptographically motivated constructions.

21 ACKNOWLEDGMENTS

We acknowledge discussions with Aram Harrow, Will Kirby, Sarah Sheldon, Kevin Smith, and Minh Tran. R.M. and V.D. acknowledge the support of the Dutch Research Council (NWO) through the NWO/NWA project Divide and Quantum (Project No. 1389.20.241). V.D. also acknowledges the support of the European Union’s Horizon Europe program through the ERC CoG BeMAIQuantum (Grant No. 101124342). This work was also supported by the Dutch National Growth Fund (NGF) as part of the QDNL program. LLM models were used to assist with numerical simulations and parts of the analytical work, including selected steps in the proofs. The authors critically reviewed and verified all incorporated suggestions and take full responsibility for the accuracy and integrity of this work.

[ACGS88] Werner Alexi, Benny Chor, Oded Goldreich, and Claus P. Schnorr. RSA and Rabin Functions: Certain Parts are as Hard as the Whole. SIAM Journal on Computing, 17(2):194–209, 1988. doi:10.1137/0217013. [AG26] Eric R. Anschuetz and Xun Gao. Arbitrary Polynomial Separations in Trainable Quantum Machine Learning. Quantum, 10:1976, 2026. doi:10.22331/q-2026-01-20-1976. [AL99] Daniel S. Abrams and Seth Lloyd. Quantum Algorithm Providing Exponential Speed Increase for Finding Eigenvalues and Eigenvectors. Physical Review Letters, 83(24):5162–5165, 1999. doi:10.1103/PhysRevLett.83.5162. [BLSF19] Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini. Parameterized quantum circuits as machine learning models. Quantum Science and Technology, 4(4):043001, 2019. doi:10.1088/2058-9565/ab4eb5. [BWP+ 17] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd. Quantum machine learning. Nature, 549(7671):195–202, 2017. doi:10.1038/nature23474. [CAB+ 21] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles. Variational quantum algorithms. Nature Reviews Physics, 3(9):625–644, 2021. doi:10.1038/s42254-021-00348-9. [CC25] Su Yeon Chang and M. Cerezo. A Primer on Quantum Machine Learning. (arXiv:2511.15969), 2025. arXiv:2511.15969. [CEMM98] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca. Quantum algorithms revisited. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 454(1969):339–354, 1998. doi:10.1098/rspa.1998.0164. [CHC+ 22] Matthias C. Caro, Hsin-Yuan Huang, M. Cerezo, Kunal Sharma, Andrew Sornborger, Lukasz Cincio, and Patrick J. Coles. Generalization in quantum machine learning from few training data. Nature Communications, 13(1):4919, 2022. doi: 10.1038/s41467-022-32550-3. [CLG+ 25] M. Cerezo, Martin Larocca, Diego Garcı́a-Martı́n, N. L. Diaz, Paolo Braccia, Enrico Fontana, Manuel S. Rudolph, Pablo Bermejo, Aroosa Ijaz, Supanut Thanasilp, Eric R. Anschuetz, and Zoë Holmes. Does provable absence of barren plateaus imply classical simulability? Nature Communications, 16(1):7907, 2025. doi:10.1038/s41467-025-63099-6. [CSV+ 21] M. Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J. Coles. Cost Function Dependent Barren Plateaus in Shallow Parametrized Quantum Circuits. Nature Communications, 12(1):1791, 2021. doi:10.1038/s41467-021-21728-w. [GGC+ 24] Jennifer R. Glick, Tanvi P. Gujarati, Antonio D. Córcoles, Youngseok Kim, Abhinav Kandala, Jay M. Gambetta, and Kristan Temme. Covariant quantum kernels for data with group structure. Nature Physics, 20(3):479–483, 2024. doi:10. 1038/s41567-023-02340-9. [GGPD25] Elies Gil-Fuster, Casper Gyurik, Adrian Pérez-Salinas, and Vedran Dunjko. On the relation between trainability and dequantization of variational quantum learning models. In The Thirteenth International Conference on Learning Representations, 2025. URL: https://openreview.net/forum?id=TdqaZbQvdi. [HBK+ 25] Laura J. Henderson, Kerstin Beer, Salini Karuvade, Riddhi Gupta, Angela White, and Sally Shrapnel. Quantum advantage without exponential concentration: Trainable kernels for symmetry-structured data. 2025. arXiv:2509.14337. [HBM+ 21] Hsin-Yuan Huang, Michael Broughton, Masoud Mohseni, Ryan Babbush, Sergio Boixo, Hartmut Neven, and Jarrod R. McClean. Power of data in quantum machine learning. Nature Communications, 12(1):2631, 2021. doi:10.1038/ s41467-021-22539-9. [HCT+ 19] Vojtech Havlicek, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow, Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta. Supervised learning with quantum enhanced feature spaces. Nature, 567(7747):209–212, 2019. doi: 10.1038/s41586-019-0980-2. [HKP20] Hsin-Yuan Huang, Richard Kueng, and John Preskill. Predicting many properties of a quantum system from very few measurements. Nature Physics, 16(10):1050–1057, 2020. doi:10.1038/s41567-020-0932-7. [JGM+ 24] Sofiene Jerbi, Casper Gyurik, Simon C. Marshall, Riccardo Molteni, and Vedran Dunjko. Shadows of quantum machine learning. Nature Communications, 15(1):5676, 2024. doi:10.1038/s41467-024-49877-8. [KB15] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. [KPM+ 26] William Kirby, Bibek Pokharel, Javier Robledo Moreno, Kevin C. Smith, Sergey Bravyi, Abhinav Deshpande, Constantinos Evangelinos, Bryce Fuller, James R. Garrison, Ben Jaderberg, Caleb Johnson, Petar Jurcevic, Su-un Lee, Simon Martiel, Mario Motta, Seetharami Seelam, Oles Shtanko, Kevin J. Sung, Minh Tran, Vinay Tripathi, Kazuhiro Seki, Kazuya Shinjo, Han Xu, Lukas Broers, Tomonori Shirakawa, Seiji Yunoki, Kunal Sharma, and Antonio Mezzacapo. Observation

22 of Improved Accuracy over Classical Sparse Ground-State Solvers using a Quantum Computer. (arXiv:2603.03496), 2026. arXiv:2603.03496. [KV94] Michael J. Kearns and Umesh Vazirani. An Introduction to Computational Learning Theory. The MIT Press, 1994. doi: 10.7551/mitpress/3897.001.0001. [LAT21] Yunchao Liu, Srinivasan Arunachalam, and Kristan Temme. A rigorous and robust quantum speed-up in supervised machine learning. Nature Physics, 17(9):1013–1017, 2021. doi:10.1038/s41567-021-01287-z. [LTW+ 25] Martı́n Larocca, Supanut Thanasilp, Samson Wang, Kunal Sharma, Jacob Biamonte, Patrick J. Coles, Lukasz Cincio, Jarrod R. McClean, Zoë Holmes, and M. Cerezo. Barren plateaus in variational quantum computing. Nature Reviews Physics, 7(4):174–189, 2025. doi:10.1038/s42254-025-00813-9. [MBS+ 18] Jarrod R. McClean, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and Hartmut Neven. Barren plateaus in quantum neural network training landscapes. Nature Communications, 9(1):4812, 2018. doi:10.1038/s41467-018-07090-4. [MGD26] Riccardo Molteni, Casper Gyurik, and Vedran Dunjko. Exponential quantum advantages in learning quantum observables from classical data. npj Quantum Information, 12(1):19, 2026. doi:10.1038/s41534-025-01162-2. [MN25] Mathias Mikkelsen and Yuya O. Nakagawa. Quantum-selected configuration interaction with time-evolved state. Physical Review Research, 7(4):043043, 2025. doi:10.1103/75pv-hbrx. [PCGL20] Adrián Pérez-Salinas, Alba Cervera-Lierta, Elies Gil-Fuster, and José I. Latorre. Data re-uploading for a universal quantum classifier. Quantum, 4:226, 2020. doi:10.22331/q-2020-02-06-226. [QWZ+ 23] Han Qi, Lei Wang, Hongsheng Zhu, Abdullah Gani, and Changqing Gong. The barren plateaus of quantum neural networks: Review, taxonomy and trends. Quantum Information Processing, 22(12):435, 2023. doi:10.1007/ s11128-023-04188-7. [RBN+ 22] Michael Ragone, Paolo Braccia, Quynh T. Nguyen, Louis Schatzki, Patrick J. Coles, Frederic Sauvage, Martin Larocca, and M. Cerezo. Representation Theory for Geometric Quantum Machine Learning. 2022. arXiv:2210.07980. [RML14] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd. Quantum Support Vector Machine for Big Data Classification. Physical Review Letters, 113(13):130503, 2014. doi:10.1103/PhysRevLett.113.130503. [SBSW20] Maria Schuld, Alex Bocharov, Krysta M. Svore, and Nathan Wiebe. Circuit-centric quantum classifiers. Physical Review A, 101(3):032308, 2020. doi:10.1103/PhysRevA.101.032308. [Sch05] U. Schollwöck. The density-matrix renormalization group. Reviews of Modern Physics, 77(1):259–315, 2005. doi:10.1103/ RevModPhys.77.259. [SEM23] Franz J. Schreiber, Jens Eisert, and Johannes Jakob Meyer. Classical Surrogates for Quantum Learning Models. Physical Review Letters, 131(10):100803, 2023. doi:10.1103/PhysRevLett.131.100803. [Sho97] Peter W. Shor. Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer. SIAM Journal on Computing, 26(5):1484–1509, 1997. doi:10.1137/S0097539795293172. [SKI+ 25] Kenji Sugisaki, Shu Kanno, Toshinari Itoko, Rei Sakuma, and Naoki Yamamoto. Hamiltonian simulation-based quantum-selected configuration interaction for large-scale electronic structure calculations with a quantum computer. Physical Chemistry Chemical Physics, 27(38):20869–20884, 2025. doi:10.1039/d5cp02202a. [SP21] Maria Schuld and Francesco Petruccione. Quantum Models as Kernel Methods, pages 217–245. Springer International Publishing, Cham, 2021. doi:10.1007/978-3-030-83098-4_6. [SSM21] Maria Schuld, Ryan Sweke, and Johannes Jakob Meyer. The effect of data encoding on the expressive power of variational quantum machine learning models. Physical Review A, 103(3):032430, 2021. doi:10.1103/PhysRevA.103.032430. [Tro12] Joel A. Tropp. User-Friendly Tail Bounds for Sums of Random Matrices. Foundations of Computational Mathematics, 12(4):389–434, 2012. doi:10.1007/s10208-011-9099-z. [Ver26] Roman Vershynin. High-Dimensional Probability. Cambridge University Press, New York, 1 edition, 2026. [Whi92] Steven R. White. Density matrix formulation for quantum renormalization groups. Physical Review Letters, 69(19):2863– 2866, 1992. doi:10.1103/PhysRevLett.69.2863. [XRV17] Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: A Novel Image Dataset for Benchmarking Machine Learning Algorithms. 2017. arXiv:1708.07747. [YGH+ 25] Qi Ye, Shuangyue Geng, Zizhao Han, Weikang Li, L.-M. Duan, and Dong-Ling Deng. Quantum automated learning with provable and explainable trainability. (arXiv:2502.05264), 2025. arXiv:2502.05264. [YMI+ 25] Jeffery Yu, Javier Robledo Moreno, Joseph T. Iosue, Luke Bertels, Daniel Claudino, Bryce Fuller, Peter Groszkowski, Travis S. Humble, Petar Jurcevic, William Kirby, Thomas A. Maier, Mario Motta, Bibek Pokharel, Alireza Seif, Amir Shehata, Kevin J. Sung, Minh C. Tran, Vinay Tripathi, Antonio Mezzacapo, and Kunal Sharma. Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization. (arXiv:2501.09702), 2025. arXiv:2501.09702. [ZZN+ 26] Haimeng Zhao, Alexander Zlokapa, Hartmut Neven, Ryan Babbush, John Preskill, Jarrod R. McClean, and Hsin-Yuan Huang. Exponential quantum advantage in processing massive classical data. (arXiv:2604.07639), 2026. arXiv:2604.07639.

23 Appendix A: Energy-based generalization bounds

The bound on the energy generalization gap stated in Theorem 2 follows from two main technical ingredients. First, for any fixed value of α, we bound the probability that the deviation ∥ HD (α) − HT (α)∥ is large over the random draw of the training set T . Second, we extend this fixed-parameter concentration bound uniformly over the entire parameter space, thereby obtaining control of the deviation simultaneously for all α in the parameter space. N where x is sampled from an underlying Consider a classification problem with training points T = { xℓ , yℓ }ℓ= 1 distribution D over an input space Ω and yℓ ∈ {0, 1}. Let HT (α) denote the training Hamiltonian defined in Eq. (1), constructed from an arbitrary encoding unitary U ( x, α) parameterized by α ∈ Θ = [−1, 1] M . For any labeled input z = ( x, y), we define Hz′ (α) = I − Hz (α), where Hz (α) = U † ( x, α)O(y)U ( x, α). We then define the population Hamiltonian as HD (α) = Ez∼D [ Hz′ (α)], namely, the average encoding Hamiltonian associated with a random labeled input drawn from the data distribution. The key idea for the proof of Theorem 2 relies on controlling how much the empirical energy cost optimized during training can differ from the expected energy of the resulting trained state when averaged over the underlying data distribution. This amounts to controlling the deviation between the training Hamiltonian HT (α) and the population Hamiltonian HD (α). Specifically, we first show that for any fixed α, the following lemma holds. Lemma 5. Let HT = I − N1 ∑ℓ Hzℓ (α) be a Hamiltonian embedding a training set T as in Eq. (1), such that 0 ⪯ Hz (α) ⪯ I for all points z = ( x, y) ∈ Ω × {0, 1} and all α ∈ Θ, and let HD (α) = Ez∼D [ Hz′ (α)] be the population Hamiltonian, with Hz′ (α) = I − Hz (α). Then, for every α ∈ Θ, with probability at least 1 − δ over the draw of N-size training datasets N , T = {( xℓ , yℓ )}ℓ= 1 s 8 2D log . (A1) H ( α ) − H ( α )∥ ≤ ∥ T D ∞ N δ where D = 2n is the Hilbert space dimension. The proof of Lemma 5 relies on the matrix Hoeffding inequality. We report here the relative theorem, the proof can be found in [Tro12]. Theorem 6 (Matrix Hoeffding (Theorem 1.3 in [Tro12])). Consider a finite sequence { Xk } of independent, random, selfadjoint matrices with dimension D, and let { Ak } be a sequence of fixed self-adjoint matrices. Assume that each random matrix satisfies Ek [ X k ] = 0

Xk2 ⪯ A2k almost surely

and

(A2)

Then, for all t ≥ 0, ( Pr

λmax

∑ Xk k

)

≥t

≤ D · e−t /8σ , where σ2 :=

∑ A2k k

(A3)

∞

In our proof of Lemma 5 we will actually need to bound the operator norm of the sum of random matrices, not just the largest eigenvalue. However, we can obtain a bound on the operator norm of a sum of random matrices ∑k Xk by applying Theorem 6 twice: one to ∑k Xk and one to − ∑k Xk as explained in the following corollary.

Corollary 7 (Matrix Hoeffding for operator norm (Corollary of Theorem 1.3 in [Tro12])). Consider a finite sequence { Xk } of independent, random, self-adjoint matrices with dimension D, and let { Ak } be a sequence of fixed self-adjoint matrices. Assume that each random matrix satisfies Ek [ X k ] = 0

Xk2 ⪯ A2k almost surely

and

(A4)

Then, for all t ≥ 0, ( Pr

∑ Xk k

) ∞

≥t

≤ 2D · e−t /8σ , where σ2 :=

∑ A2k k

(A5)

∞

Proof. For a self-adjoint matrix A, the operator norm satisfies ∥ A∥∞ = max {λmax ( A), λmax (− A)} . Therefore, ∥ A∥∞ ≥ t is equivalent to λmax ( A) ≥ t

λmax (− A) ≥ t.

(A6)

24 Applying Tropp’s one-sided matrix Hoeffding bound to ∑k Xk gives " ! #

t2 ≥ t ≤ Dexp − 2 8σ

∑ Xk

Pr λmax

Applying the same bound to − ∑k Xk gives " Pr λmax

− ∑ Xk

Hence, by the union bound, " Pr

∑ Xk k

# ∞

t2 ≥ t ≤ Dexp − 2 8σ #

≥ t ≤ Pr λmax

∑ Xk

≤ 2Dexp −

k t2

8σ2

(A7)

(A8)

− ∑ Xk

≥ t + Pr λmax

≥t

(A9)

We are now ready to give the complete proof of our Lemma 5 Proof of Lemma 5. We define Xz (α) = Hz′ (α) − HD (α). Note that HT (α) − HD (α) = N1 ∑z∈T Xz (α) and that Ez∼D [ Xz (α)] = Ez∼D [ Hz′ (α)] − HD (α) = HD (α) − HD (α) = 0.

(A10)

Moreover, − I ⪯ Xz ⪯ I since we assumed 0 ⪯ Hz′ ⪯ I for any point z and therefore also 0 ⪯ Hd ⪯ I. Hence, ∥ Xz ∥∞ ≤ 1 for any data point z. We can then bound the difference in norm between HT (α) and HD (α) by

∥ HT (α) − HD (α)∥∞ ≤

1 Xz ( α ) N z∑ ∈T

In order to prove the lemma is then sufficient to bound the quantity

(A11) ∞

1 N ∑z∈T Xz ( α ) ∞ . Since the samples z ∈ T are

drawn independently, the matrices Xz (α) are also independent for every fixed α. Moreover, under the assumption 0 ⪯ Hz (α) ⪯ I, they satisfy ∥ Xz (α)∥∞ ≤ 1 for all z. We can therefore apply the matrix Hoeffding inequality of Corollary 7. In particular, in our case Xk2 ⪯ I and therefore σ2 = ∑kN=1 I " Pr T

1 N Xi N i∑ =1

= N to obtain:

Nε2 ≥ ε ≤ 2Dexp − 8 #

∞

where D = 2n is the Hilbert space dimension. Finally, after setting the right hand to be δ, we get s " # 1 N 8 2D PrT Xi ≤ log ≥ 1−δ. N i∑ N δ =1

(A12)

(A13)

∞

which, combined with Eq. (A11) proves the bound stated in Lemma 5. Lemma 5 bounds the probability of having a large difference in operator norm between the training and population Hamiltonian for a fixed α. Importantly, this probability is taken over different random datasets T . In general, the set of trained parameter α⋆ will depend on the particular training set used in the training, therefore when averaging over different T this dependence must be taken into account. When the Hamiltonian embedding is Lipschitz continuous in α, the fixed-α bound can be extended uniformly over the whole parameter space. This is done by first applying the bound to all points of a finite covering net and then using the Lipschitz property to control the deviation at parameters outside the net. This is shown in the full proof of Theorem 2, that we restate here for convenience.

25 Theorem 2 (Uniform energy-based generalization bound). Assume a Hamiltonian embedding model such that 0 ⪯ Hz (α) ⪯ I and that there exists a constant L > 0 such that ∥ Hz (α) − Hz ( β)∥∞ ≤ L∥α − β∥∞ for all points z = ( x, y) ∈ Ω × {0, 1} and all α, β ∈ [−1, 1] M . Then, for every ε ∈ (0, 4L], with probability at least 1 − δ over the draw of N-size training datasets, s 8 2D 24L ε log + M log + . (32) sup ∥ HT (α) − HD (α)∥∞ ≤ N δ ε 4 α∈[−1,1] M where D = 2n is the Hilbert space dimension. Proof. Let the parameters take values in Θ = [−1, 1] M . Equip the parameter space with the ℓ∞ norm ∥α − β∥∞ := max j=1,...,M |α j − β j |. Then [−1, 1] M is exactly the radius-one ℓ∞ ball in R M . Let η ∈ (0, 1]. An η-net Aη ⊂ Θ is a finite set such that for every α ∈ Θ, there exists β ∈ Aη satisfying ∥α − β∥∞ ≤ η. It is known [Ver26] that the covering numbers of a ball BR of radius R > 0 around a fixed point in R M , for 0 ≤ η ≤ R is:

N ( BR , ||.||, η ) = |Aη | ≤

2R η

≤

M 3 . η

(A14)

Where in the last step we consider R = 1, as we are considering a ball of radius R = 1 around 0 (that corresponds

to Θ = [−1, 1] M ). For each fixed β ∈ Aη , Eq. (A1) gives ProbT [∥ HT ( β) − HD ( β)∥ ≥ t] ≤ 2Dexp − Nt8 bound the probability that there is a large deviation for at least a point in the net. Define the bad event Eβ := {∥ HT ( β) − HD ( β)∥ ≥ t} . The event that at least one net point is bad is independent. It gives " ProbT

β∈Aη Eβ .

. We now (A15)

The union bound does not require the events Eβ to be #

sup ∥ HT ( β) − HD ( β)∥ ≥ t ≤

β∈Aη

∑ ProbT (Eβ )

(A16)

β∈Aη

Nt2 ≤ 2D |Aη |exp − 8 M 3 Nt2 ≤ 2D exp − . η 8

To make the right-hand side at most δ, it is enough to choose s 8 2D 3 t= log + M log . N δ η

(A17) (A18)

(A19)

Hence, with probability at least 1 − δ, s sup ∥ HT ( β) − HD ( β)∥ ≤

β∈Aη

8 2D 3 log + M log . N δ η

(A20)

The finite-net bound is not yet enough: the final trained parameter α⋆ (T ) need not lie exactly on the grid of points in Aη . However, we can extend the bound to a supremum of all the parameter points in α ∈ [−1, 1] M by assuming regularity of the Hamiltonian family. In particular, Lipschitz continuity implies the existence of a L ≥ 0 such that for every data point z and every α, β ∈ [−1, 1] M it holds that ∥ Hz (α) − Hz ( β)∥∞ ≤ L|α − β|∞ . We notice that both the training Hamiltonian HT (α) and the full population Hamiltonian inherit the same Lipschitz constant:

∥ HT (α) − HT ( β)∥ ≤

1 N ∥ Hzi (α) − Hzi ( β)∥ N i∑ =1

≤ L∥α − β∥∞ .

(A21) (A22)

26 and, similarly,

∥ HD (α) − HD ( β)∥ ≤ Ez [∥ Hz (α) − Hz ( β)∥] ≤ L∥α − β∥∞ .

(A23) (A24)

Define ∆ T (α) := HT (α) − HD (α). Combining Eqs. (A22) and (A24),

∥∆ T (α) − ∆ T ( β)∥ ≤ ∥ HT (α) − HT ( β)∥ + ∥ HD (α) − HD ( β)∥ ≤ 2L∥α − β∥∞ .

(A25) (A26)

Now take an arbitrary α ∈ [−1, 1] M . Choose β ∈ Aη such that ∥α − β∥∞ ≤ η. Then

∥∆ T (α)∥ ≤ ∥∆ T ( β)∥ + ∥∆ T (α) − ∆ T ( β)∥ ≤ ∥∆ T ( β)∥ + 2Lη.

(A27) (A28)

Taking the supremum over all α gives sup α∈[−1,1] M

∥∆ T (α)∥ ≤ sup ∥∆ T ( β)∥ + 2Lη.

(A29)

β∈Aη

We can now conclude the proof of the Theorem by inserting Eq. (A20) in Eq. (A29) in order to obtain: s sup α∈[−1,1] M

∥ HT (α) − HD (α)∥ ≤

2D 24L ε 8 log + M log + . N δ ε 4

(A30)

ε . This concludes the proof of Theorem 2. where we chose η = 8L

Theorem 2 provides a bound on the energy generalization gap defined in Eq. (31), it holds for any α and therefore even for the α⋆ obtained after the training. The energy generalization gap can be directly used to bound the total risk R defined in Eq. (30) as we show in Corollary 3. Corollary 3 (Bound on the total risk). Assume the setting of Theorem 2. Let (ψT ,α⋆ , α⋆ ) be the output of the trainb(ψT ,α⋆ , α⋆ ) = ⟨ψT ,α⋆ | HT (α⋆ ) |ψT ,α⋆ ⟩ be the empirical energy reached by the trained model, and let ing procedure. Let C h i ⋆ C (ψT ,α⋆ , α ) = ⟨ψT ,α⋆ | HD (α⋆ ) |ψT ,α⋆ ⟩ be its population energy. Moreover, let R(ψT ,α⋆ , α⋆ ) = Prx∼D f˜ψ ⋆ ,α⋆ ( x ) ̸= f ( x ) T ,α

(33)

(34)

Proof. Theorem 2 bounds the energy generalization gap for any parameter α ∈ Θ = [−1, 1] M and, in particular, for the α⋆ obtained after the model is trained. In particular, for any ϵ > 0, with probability larger than 1 − δ over different training data T , it holds that s sup α∈[−1,1] M

∥ HT (α) − HD (α)∥ ≤

8 2D 24L ε log + M log + . N δ ε 4

Suppose now that the number of training samples satisfies 128 2D 24L N ≥ 2 log + M log , δ ε ε

(A31)

(A32)

27 then Theorem 2 implies that, with probability 1 − δ, supα∈[−1,1] M ∥ HT (α) − HD (α)∥ ≤ 2ε . Recalling the definition of total cost, we obtain: C (ψT ,α⋆ , α⋆ ) = ⟨ψT ,α⋆ | HD (α⋆ ) |ψT ,α⋆ ⟩ ≤ ⟨ψT ,α⋆ | HT (α⋆ ) |ψT ,α⋆ ⟩ + ∥ HD (α⋆ ) − HT (α⋆ )∥ ε ≤ Cb(ψT ,α⋆ , α⋆ ) + . 2

(A33) (A34) (A35)

It now remains to relate the population energy to the total classification error R. Recall that, for a data point z = ( x, y) embedded through the Hamiltonian Hz (α), the energy loss evaluated on the state ψ is defined in Eq. (7) as:

ℓenergy (ψ, Hz (α)) = 1 − ⟨ψ| Hz (α)|ψ⟩.

(A36)

Moreover, the model prediction rule for a point x is given by Eq. (5), which we also state here for convenience: feψT ,α⋆ ,α⋆ ( x ′ ) = argmax PψT ,α⋆ ,α⋆ ( x ′ , ỹ).

(A37)

ỹ∈{0,1}

where PψT ,α ,α ( x, y) = Tr O(y)U ( x, α)|ψT ,α ⟩⟨ψT ,α |U † ( x, α) . We notice that if ( x, y) is misclassified, then the model assigns at least as much probability to the wrong label as to the true label, in other words 1 (A38) ℓenergy (ψ, Hz (α)) = 1 − ⟨ψT ,α | Hz (α) |ψT ,α ⟩ ≥ . 2 for misclassified point z. This directly implies that 1{ f˜α,ψT ,α ( x ) ̸= f ( x )} ≤ 2ℓenergy (ψ, Hz (α)) for any point z = ( x, y). Therefore, taking the expectation over different point z for x ∼ D and remembering that R(ψ, α) = Prx∼D [ feψT ,α ,α ( x ) ̸= f ( x )], we obtain: R(ψT ,α⋆ , α⋆ ) ≤ 2C (ψT ,α⋆ , α⋆ ) ≤ 2Cb(ψT ,α⋆ , α⋆ ) + ε.

(A39) (A40)

where in the second equation we used the bound on the population cost in Eq. (A33) given by the training cost plus the energy generalization gap assuming the sample size given in Eq. (A32).

Lipschitz assumption

In Theorem 2, we assumed that the parameterized Hamiltonians used to encode the data satisfy a Lipschitz continuity condition. Namely, we required that there exists a constant L > 0 such that ∥ Hz (α) − Hz ( β)∥ ≤ L∥α − β∥∞ for any α, β ∈ [−1, 1] M and all the labeled input z = ( x, y). We now show that this is a natural assumption, and that it holds for common choices of encoding unitaries U ( x, α) used in the construction of Hz (α) = U † ( x, α)O(y)U ( x, α). Suppose that the data are encoded through rotation angles given by linear combinations of the input features, for example U (α, x ) = e−iθα ( x)G ,

θα ( x ) = ∑ α j x j .

(A41)

j =1

If Hz (α) = U ( x, α)† O(y)U ( x, α) with ∥O(y)∥ ≤ 1, then differentiating with respect to α j gives a commutator-type

expression. Since

∂θα ( x ) ∂α j = x j , one obtains

∂Hz (α) ∂α j

≤ 2| x j |∥ G ∥.

(A42)

≤ 2Bx ∥ G ∥.

(A43)

If the features are uniformly bounded, | x j | ≤ Bx , then ∂Hz (α) ∂α j

28 Therefore, moving coordinate by coordinate from α to β, M

∥ Hz (α) − Hz ( β)∥ ≤ 2Bx ∥ G ∥ ∑ |α j − β j |.

(A44)

∥ Hz (α) − Hz ( β)∥ ≤ 2MBx ∥ G ∥∥α − β∥∞ .

(A45)

j =1

In particular,

Thus, for Pauli generators with ∥ G ∥ = 1, one may take L ≤ 2MBx , up to convention-dependent factors such as the presence of a 1/2 in the rotation generator. Appendix B: Gradient evaluation and the Hellmann–Feynman relation

Here we justify why, when evaluating the gradient of the cost function in Eq. (6), the ground state can be treated as fixed at the value corresponding to the current parameters, despite its implicit dependence on α. We formalize the result in Theorem 8. N Theorem 8 (Energy cost gradient). Assume that the data Hamiltonian HT (α) = I − N1 ∑ℓ= 1 Hzℓ ( α ), as defined in Eq. (1), is continuously differentiable with a nondegenerate ground state. Define the energy cost as in Eq. (6)

bT (ψT ,α , α) = ⟨ψT ,α | HT (α)|ψT ,α ⟩. C

(B1)

Then, for every parameter α j , ∂ b C (ψT ,α , α) ∂α j T

α=αt

=−

∂Hzℓ (α) 1 N ⟨ψT ,αt | ∑ N ℓ=1 ∂α j

α=αt

|ψT ,αt ⟩.

(B2)

Proof. Remember that the cost function as a function of the parameter α is defined as: bT (ψT ,α , α) = ⟨ψT ,α | HT (α) |ψT ,α ⟩ , C

(B3)

where |ψT ,α ⟩ is a normalized ground state of HT (α), satisfying HT (α) |ψT ,α ⟩ = E0 (α) |ψT ,α ⟩ with E0 (α) the ground energy of the data Hamiltonian HT (α) at the point α. bT (ψT ,α , α) Assuming that the ground state is non-degenerate and differentiable with respect to α, differentiating C with respect to α j gives bT (ψT ,α , α) ∂C ∂HT (α) = ∂ j ψT ,α HT (α) |ψT ,α ⟩ + ⟨ψT ,α | |ψT ,α ⟩ ∂α j ∂α j

+ ⟨ψT ,α | HT (α) ∂ j ψT ,α ,

(B4)

where ∂ j ψT ,α ≡ ∂ |ψT ,α ⟩ /∂α j . Using the eigenvalue equation, the first and third terms combine to give E0 (α)

∂ j ψT ,α ψT ,α + ψT ,α ∂ j ψT ,α

(B5)

Since |ψT ,α ⟩ is normalized, differentiating ⟨ψT ,α |ψT ,α ⟩ = 1 yields ∂ j ψT ,α ψT ,α + ψT ,α ∂ j ψT ,α = 0. Therefore, bT (ψT ,α , α) ∂C ∂HT (α) = ⟨ψT ,α | |ψT ,α ⟩ , ∂α j ∂α j

(B6)

which is the Hellmann-Feynman theorem. In particular, at iteration t of Algorithm 1, after computing the ground state |ψT ,αt ⟩ of HT (αt ), Eq. (B6) gives bT (ψT ,α , α) ∂C ∂α j

α=αt

= ⟨ψT ,αt |

∂HT (α) ∂α j

α=αt

|ψT ,αt ⟩ ,

(B7)

29 Hence, although the ground state depends implicitly on the parameters α, its derivative does not need to be evaluated when computing the gradient of the ground-state energy at the current parameter value. Equivalently, the bT (ψT ,α , α) has the same gradient as the ground-state energy E0 (α) at α = αt . frozen-state cost C t N Finally, using HT (α) = I − N1 ∑ℓ= 1 Hzℓ ( α ), we obtain bT (ψT ,α , α) ∂C ∂α j

α=αt

=−

∂Hz (α) 1 N ∑ ⟨ψT ,αt | ∂αℓ j N ℓ= 1

α=αt

|ψT ,αt ⟩ ,

(B8)

which gives the gradient used in Eq. (21) in the main text. 1.

Gradients with approximate ground states

In practice, the ground states in the optimization loop are obtained through quantum algorithms like QPE or eα of HT (α). Defining the corresponding approximate SKQD. These may return only an approximate ground state ψ energy as e(α) = ψ eα , eα HT (α) ψ E

(B9)

its exact derivative with respect to α j is e(α) ∂HT (α) e ∂E ψα = ∂ j ψeα HT (α) ψeα + ψeα ∂α j ∂α j

+ ψeα HT (α) ∂ j ψeα ,

(B10)

eα is an exact eigenstate, the Hellmann-Feynman cancellation does eα /∂α j . Therefore, unless ψ eα ≡ ∂ ψ where ∂ j ψ e(α). not hold exactly for the derivative of E Nevertheless, the gradient obtained by keeping the approximate ground state fixed provides a controlled approximation to the exact ground-state gradient. Let |ψT ,α ⟩ denote the exact ground state of HT (α) and define its fidelity eα |2 . Using Hölder’s inequality together with the pure-state tracewith the approximate state as F (α) = | ψT ,α ψ p eα ψ eα − |ψT ,α ⟩⟨ψT ,α | = 2 1 − F (α), we obtain distance relation ψ 1 eα ∂HT (α) ψ eα − ⟨ψT ,α | ∂HT (α) |ψT ,α ⟩ ≤ 2 ∂HT (α) ψ ∂α j ∂α j ∂α j

1 − F ( α ).

Using the Hellmann–Feynman relation in Eq. (B6), this gives eα ψ

∂HT (α) ∂E (α) ∂HT (α) eα − 0 ψ ≤2 ∂α j ∂α j ∂α j

1 − F ( α ).

(B11)

Hence, if the fidelity of the ground state found by the quantum solver√is F (α) ≥ 1 − ϵ, the error in the gradient ϵ. Moreover, as discussed in Appendix C, obtained from the approximate ground state is at most 2 ∂HT (α)/∂α j for the standard Pauli encoding introduced in the main text, the operator norm ∂HT (α)/∂α j is bounded. Appendix C: Absence of barren plateaus in the Hamiltonian embedding model

One appealing feature of Hamiltonian embedding is that, for local encoding Hamiltonians, once the datasetencoding ground state has been obtained, inference reduces to estimating a local observable on the learned state. Together with enabling classical-shadow-based inference, this locality structure also allows us to establish trainability guarantees for the model. In particular, for a particular class of Pauli-rotation encoding U ( x, α), we can theoretically prove that a nontrivial subset of the trainable parameters avoids barren plateaus. We present the proof for the simple observable A = 1/n ∑i Zi and a family of shallow Hamiltonian embedding architectures. The same argument extends to observables that are sums of bounded-locality Pauli terms, provided that each Pauli term entering the bound anticommutes with at least one final-layer rotation generator acting on its support.

30 Theorem 9 (Cost-level absence of barren plateaus). Let U ( x, α) be a depth-K layered unitary of the form U ( x, α) = K ( x ) · · · V 1 ( α1 )U 1 ( x ), where α = ( α1 , . . . , αK ) contains M trainable parameters. Assume that A = 1 n Z V K (αK )Uenc enc n ∑ i =1 i and the final trainable layer contains one independent X- or Y-rotation on each readout qubit, parameterized by the angles α1K , . . . , αnK with αiK ∈ [0, 2π ]. Let h i b(ψT ,α , α) = E( x,y)∈T ℓeng ψT ,α , H( x,y) (α) Ct (α) := C (C1) t t be the empirical energy cost at optimization step t, where the state ψT ,αt is kept fixed when differentiating with respect to α. Assume that, for some γ > 0, Ct (αt ) ≤

1 − γ. 2

(C2)

Then the coordinate-averaged variance of the gradients of the cost function satisfies 1 M Varαµ M µ∑ =1

∂Ct (α) ∂αµ

≥

γ2 , 2Mn

(C3)

where Varαµ denotes the coordinate-wise variance with respect to αµ ∼ Unif([0, 2π ]) while keeping all the other parameters fixed at their optimization-step values. Proof. At optimization step t, we keep the state |ψT ,αt ⟩ fixed and optimize only the circuit parameters appearing in the data-dependent Hamiltonian. For a single training example zℓ = ( xℓ , yℓ ), the fixed-state energy loss is defined as in Eq. (7):

ℓeng (ψT ,αt , Hzℓ (α)) = 1 − ⟨ψT ,αt | Hzℓ (α) |ψT ,αt ⟩ ,

(C4)

where the state ψT ,αt is kept fixed when differentiating with respect to α. Let ρt = |ψT ,αt ⟩ ⟨ψT ,αt |. Recalling the I +(−1)y A and A = n1 ∑nj=1 Zj , we define expression for HT (α), with O(y) = 2 h i h j ( x; α) := Tr ρt U ( x, α)† Zj U ( x, α) . Let us define s( x ) = (−1)y( x) . Then the fixed-state single-sample energy loss is ℓeng (ψT ,αt , H( x,y) (α)) = s( x ) n 2n ∑ j=1 h j ( x; α ), while the frozen total energy cost at optimization step t is

Ct (α) =

1 1 n − Ex∈T s( x )h j ( x; α) . ∑ 2 2n j=1

(C5) 1 2 −

(C6)

The average is taken over the training set T . By the assumption in Eq. (C2), Ct (αt ) ≤ 12 − γ. Equivalently, 1 n Ex∈T s( x )h j ( x; αt ) ≥ 2γ. n j∑ =1

(C7)

We now isolate the dependence of the encoding unitary U ( x, α) on the final trainable layer. After fixing all previous layers, including both the data-encoding gates and the trainable gates in layers 1, . . . , K − 1, we write U ( x, α) = VK (αK )U<K ( x, α<K ), where U<K ( x, α<K ) denotes the product of all layers before the final one, and VK (αK ) is the K last trainable layer. At optimization step t, the parameters α<K are fixed at α< t , so we absorb them into the frozen < K unitary Wt ( x ) := U<K ( x, αt ). Thus, for the purpose of analyzing the final layer at optimization step t, we write U ( x, αK ) = VK (αK )Wt ( x ). In particular, we take the final layer to consist of one trainable single-qubit rotation on each readout qubit, αK = (α1K , . . . , αnK ). By assumption, the rotation acting on qubit i is generated by a Pauli operator that does not commute with the measured observable Zi . For example, if the final gate on qubit i is an Xi -rotation, then R Xi (αiK )† Zi R Xi (αiK ) = Zi cos αiK + Yi sin αiK . Equivalently, after freezing ρt , the dataset, all earlier parameters, and all final-layer parameters except αiK , there exist real coefficients ai ( x ) and bi ( x ) such that we can rewrite hi in Eq. (C5) as hi ( x; α) = ai ( x ) cos αiK + bi ( x ) sin αiK .

(C8)

31 Note that the coefficients ai ( x ) and bi ( x ) are independent of αiK . We now pass to the empirical cost by averaging ℓeng over the input x. Since αiK is independent of the data point x, we can define the labeled empirical coefficients Ai := Ex∈T [s( x ) ai ( x )] and Bi := Ex∈T [s( x )bi ( x )] and obtain that: (C9) Ex∈T [s( x )hi ( x; α)] = Ai cos αiK + Bi sin αiK . q 2 2 Define the empirical amplitude Ri := Ai + Bi . By the amplitude-phase form of a sinusoid, there exists a phase φi such that Ex∈T [s( x )hi ( x; α)] = Ri cos αiK + φi . (C10) This implies that, for every value of αiK ,

|Ex∈T [s( x )hi ( x; α)]| ≤ Ri .

(C11)

Ex∈T [s( x )hi ( x; αt )] ≤ |Ex∈T [s( x )hi ( x; αt )]| ≤ Ri .

(C12)

Evaluating Eq. (C11) at α = αt we obtain:

Therefore, n1 ∑in=1 Ri ≥ 2γ and, by Jensen’s inequality, 1 n 2 Ri ≥ n i∑ =1

1 n Ri n i∑ =1

≥ 4γ2 .

(C13)

Equivalently, ∑in=1 R2i ≥ 4nγ2 . Using Eq. (C9), the frozen empirical cost can be written as Ct (αK ) =

1 1 n Ai cos αiK + Bi sin αiK , − ∑ 2 2n i=1

(C14)

K where all earlier-layer parameters are fixed at α< t . The derivative of the frozen empirical cost with respect to the final-layer parameter αiK can be computed using Eq. (C14)

∂Ct (αK ) 1 K K A sin α − B cos α . = i i i i 2n ∂αiK

(C15)

Following the standard barren-plateau analysis, we evaluate the gradient variance with respect to independently and uniformly distributed final-layer angles αiK ∼ Unif([0, 2π ]), while keeping all other parameters fixed. Since " # ∂Ct Eα K = 0, (C16) i ∂αiK we obtain VarαK i

∂Ct ∂αiK

!2  ∂Ct  = Eα K  i ∂αiK 2 1 K K = 2 EαK Ai sin αi − Bi cos αi i 4n 

A2i + Bi2 R2 = i2 . 2 8n 8n

Summing over the n final-layer parameters and using Eq. (C13), we therefore find ! n ∂Ct 1 n γ2 ∑ VarαiK ∂αK = 8n2 ∑ R2i ≥ 2n . i =1 i =1 i

(C17)

(C18)

32 Finally, let M denote the total number of trainable parameters in the circuit. Since the variance associated with every parameter direction is nonnegative, the average variance over all trainable parameters is lower bounded by the contribution of the n final-layer parameters alone: ! ∂Ct ∂Ct 1 n 1 M Varαµ VarαK ≥ i M µ∑ ∂αµ M i∑ ∂αiK =1 =1

≥

γ2 . 2Mn

(C19)

This proves Eq. (C3). If the ansatz has depth K and O(n) trainable parameters per layer, then M = O(Kn). Hence, there exists a constant c > 0 such that 1 M ∂Ct cγ2 Var . (C20) ≥ αµ ∑ M µ =1 ∂αµ Kn2

We now show that the results in Theorem 9 are enough to guarantee that a non-trivial fraction of parameters in the model have gradient variance that remains inverse-polynomial in the system size, and therefore avoid barren plateaus. Corollary 10 (Non-negligible fraction of parameters avoids barren plateaus.). Under the assumptions of Theorem 9, let ( ! ) γ2 ∂Ct Sγ := i ∈ {1, . . . , n} : VarαK ≥ 2 . (C21) i 4n ∂αiK Then

| Sγ | ≥

2γ2 n. 1 − 2γ2

(C22)

Hence, if the total number of trainable parameters satisfies M = O(Kn), an Ω(γ2 /K ) fraction of all trainable parameters has gradient variance at least γ2 /(4n2 ). In particular, for K = O(1) number of layers and γ = Θ(1), a constant fraction of the trainable parameters has at least inverse-polynomial gradient variance. Proof. From Eq. (C13) in Theorem 9, we have ∑in=1 R2i ≥ 4nγ2 . Moreover, Ri ≤ 1 for every i. Indeed, by Eq. (C10), Ri = max |Ex∈T [s( x )hi ( x; α)]| , αiK

(C23)

and since s( x ) ∈ {−1, +1} and hi ( x; α) = Tr[ρt U ( x, α)† Zi U ( x, α)] ∈ [−1, 1], the absolute value of the expectation is at most one. Define the set of indices i whose corresponding R2i is greater than 2γ2 n o Seγ := i ∈ {1, . . . , n} : R2i ≥ 2γ2 ,

(C24)

and let m := |Seγ |. For the indices in Seγ we use the bound R2i ≤ 1 derived above, while for the remaining n − m indices we have R2i < 2γ2 . Therefore, n

∑ R2i ≤ m + 2γ2 (n − m).

(C25)

i =1

Combining this with Eq. (C13) gives 4nγ2 ≤ m + 2γ2 (n − m), and hence m(1 − 2γ2 ) ≥ 2nγ2 .

(C26)

33 Since the assumption Ct (αt ) ≤ 1/2 − γ implies 0 < γ ≤ 1/2, we have 1 − 2γ2 > 0, and thus the number of indices i whose R2i is greater than 2γ2 is m≥

2γ2 n. 1 − 2γ2

(C27)

For every i ∈ Seγ , the same argument as in Theorem 9 allows us to bound the variance of the gradient for the parameter αiK using the quantity Ri ! R2i ∂Ct γ2 VarαK = ≥ . (C28) i 8n2 4n2 ∂αiK Therefore Seγ ⊆ Sγ , and hence

|Sγ | ≥ |Seγ | ≥

2γ2 n. 1 − 2γ2

(C29)

Finally, if M = O(Kn), then 2 | Sγ | γ =Ω . M K

(C30)

Thus, for constant depth K = O(1) and γ ≥ 1/ poly(n), an inverse-polynomial fraction of all trainable parameters has inverse-polynomial gradient variance. Our result concerns the cost evaluated for different α with the prepared ground state held fixed, as in line 21 of Algorithm 1. It shows that varying a non-negligible fraction of the parameters can produce cost differences of at least inverse-polynomial size. We note that, in general, the absence of barren plateaus alone does not guarantee successful training, as optimization may still encounter local minima or other stationary points. Appendix D: Universality of the Hamiltonian Embedding model

In this Appendix, we show that the Hamiltonian embedding model allows for universal construction. In particular, given a measurable target function f : Ω ⊂ Rd → {0, 1} with inputs from Ω = [−πL, πL]d , where L > 0 is an arbitrary parameter determining the size of the input domain, we construct an HE model with a fixed product encoding U ( x ) whose ground-state classifier approximates f up to ε error under the uniform input distribution on Ω. Theorem 11 (Universality with constant-depth encodings). Let Ω = [−πL, πL]d be the input space, µ be the uniform distribution on Ω and let f : Ω → {0, 1} be a measurable function. For every ε, δ ∈ (0, 1), there exist finite n, N ∈ N and a N , with x ∼ µ the following holds. Let constant depth encoding U ( x ) such that given a training dataset T = {( xℓ , f ( xℓ ))}ℓ= 1 |ψT ⟩ be a ground state of HT = I −

1 N ∑ U (xℓ )† O( f (xℓ ))U (xℓ ), N ℓ= 1

(D1)

where O(y) = ( I + (−1)y A)/2 with A = n−1 ∑in=1 Zi . Then, the HE classifier f˜T defined by the prediction rule in Eq. (5), using the fixed encoding U ( x ) and the ground state of HT , satisfies with at least 1 − δ probability over the training dataset, Pr [ f˜T ( x ) ̸= f ( x )] ≤ ε.

x ∼µ

(D2)

Proof. First, we write the target function as s( x ) = (−1) f ( x) and compute its Fourier coefficients b sk = Eµ [s( x )e−ik· x/L ] = iϕ |b sk |e k , where µ is the uniform distribution over Ω. Since s is bounded and measurable, it belongs to L2 (µ). Completeness of the Fourier basis therefore ensures that, for every ε ∈ (0, 1), a finite real Fourier sum p exists: p( x ) = ∑ ak cos(k · x/L + ϕk ), k∈K

(D3)

34 where ak = 2|b sk | for k ̸= 0, a0 = |b s0 | for the constant term, and ϕk = arg b sk . Parseval’s identity then guarantees that with a finite |K | it holds that ∥s − p∥2L2 (µ) ≤ ε/4. Here K is the set of frequencies retained from the original Fourier decomposition of s, excluding those with zero Fourier coefficients and keeping one frequency from each nonzero pair {k, −k }. For example, we can take K to contain frequencies k ∈ Zd whose entries have absolute value at most a cutoff M, chosen such that ∥s − p∥2L2 (µ) ≤ ε/4. We now construct the encoding U such that the HE model approximates the function p. Set W = ∑k∈K ak . For each frequency k ∈ K, we choose integer multiplicities mk such that ∑k∈K mk = n and |mk /n − ak /W | ≤ 1/n, where n is the number of qubits we use to encode the HE model. We then define the following HE model U (x) =

O k∈K

RY (k · x/L)⊗mk ,

1 n Zi , n i∑ =1

O(y) =

I + (−1)y A . 2

(D4)

Thus, each frequency k ∈ K appearing in the truncated Fourier series p specifies a single-qubit rotation, and its Fourier magnitude determines the fraction of qubits encoding that frequency. The phases ϕk will be assigned during training in the following way. Let us consider a draw of independent inputs xℓ ∼ µ and labeled by yℓ = f ( xℓ ). Define the empirical correlations Ck =

1 N ∑ (−1)yℓ e−ik·xℓ /L . N ℓ= 1

(D5)

Notice that Ck approaches the Fourier coefficient b sk of the target s as N → ∞. The real and imaginary part of Ck are given by Euler’s formula: Re(Ck ) =

1 N ∑ (−1)yℓ cos(k · xℓ /L), N ℓ= 1

Im(Ck ) = −

1 N ∑ (−1)yℓ sin(k · xℓ /L). N ℓ= 1

(D6)

Since the encoding in Eq. (D4) acts independently on each qubit, the identity RY (θ )† ZRY (θ ) = cos θ Z − sin θ X gives for each qubit i: U ( x )† Zi U ( x ) = cos(k i · x/L) Zi − sin(k i · x/L) Xi ,

(D7)

where k i is the frequency assigned to qubit i. We can now write the embedded Hamiltonian in terms of Ck . Substituting these identities and the readout O(y) into the definition of HT in Eq. (1) yields HT = I −

1 N 1 1 n † U ( x ) O ( y ) U ( x ) = I − Re(Cki ) Zi + Im(Cki ) Xi . ∑ ∑ ℓ ℓ ℓ N ℓ=1 2 2n i=1

(D8)

By writing Ck in polar form, Ck = |Ck |eiθk , it is easy to show that the ground state of the embedded Hamiltonian in Eq. (D8) is given by product state of one-qubit states RY (arg Cki )|0⟩, for the qubits encoding the frequency k i . In fact, each one-qubit Hamiltonian term is real and symmetric in the computational basis, so a ground state can be chosen with real amplitudes. Since every normalized real qubit state can be written as RY ( β)|0⟩, it suffices to minimize the energy over this one-parameter family. Its Pauli expectations are ⟨ Z ⟩ = cos β and ⟨ X ⟩ = sin β, so the corresponding energy term has expectation Ek ( β) =

1 1 [1 − |Ck | (cos θk cos β + sin θk sin β)] = [1 − |Ck | cos( β − θk )]. 2n 2n

(D9)

This reaches the lowest eigenvalue, Ek = 1/(2n)[1 − |Ck |], when β = θk . Therefore, the one-qubit ground state is RY (arg Ck )|0⟩. Since the Hamiltonian terms act on separate qubits, their ground states can be combined into a product ground state. It is important to note that the empirical phase θk is distinct from the target Fourier phase ϕk = arg b sk . However, as Ck → b sk , the learned phase θk approaches ϕk modulo 2π. We now show that the HE classifier obtained from the learned ground state approximates the target labels encoded N by s. The full ground state is |ψT ⟩ = k∈K ( RY (arg Ck )|0⟩)⊗mk , which implements the following prediction score function: m (D10) gT ( x ) := ⟨ψT |U ( x )† AU ( x )|ψT ⟩ = ∑ k cos(k · x/L + arg Ck ) . n k∈K

35 Indeed, the encoding rotation and the ground-state rotation combine into RY (k · x/L + arg Ck ) on each qubit, whose Z expectation is the corresponding cosine. The model predicts label 0 when gT ( x ) ≥ 0 and label 1 otherwise. For fixed n, define the limiting score for N → ∞ mk cos(k · x/L + ϕk ). n k∈K

g∞ ( x ) = ∑

(D11)

Since each cosine has magnitude at most one, Wmk W |K | − ak ≤ . n n k∈K

(D12)

√

(D13)

∥Wg∞ − p∥∞ ≤ ∑

√ We can then choose n such that this bound is at most ε/4. Secondly, for a finite dataset, arg Ck will only be an approximation of ϕk = arg b sk . However, from the fact that Ck → b sk as N → ∞, we have that for any δ > 0, we can √ first choose n such that W |K |/n ≤ ε/4, and then a finite N such that, with probability at least 1 − δ, ∥WgT − Wg∞ ∥∞ ≤

ε/4.

On this event, we can bound the average approximation error of the HE model with respect the target s as:

∥s − WgT ∥ L2 (µ) ≤ ∥s − p∥ L2 (µ) + ∥ p − Wg∞ ∥∞ + ∥Wg∞ − WgT ∥∞ √ √ √ √ ε ε ε ≤ + + = ε. 2 4 4

(D14) (D15)

Since W > 0, a mistaken prediction implies |s( x ) − WgT ( x )| ≥ 1. Therefore, Prx∼µ [ f˜T ( x ) ̸= f ( x )] ≤ ∥s − WgT ∥2L2 (µ) ≤ ε. The construction guarantees universality with finite, target-dependent number of qubits n and training points N. We note that a polynomial scaling of n and N is not guaranteed for arbitrary target functions.

Appendix E: Solving a classically hard problem with Hamiltonian embedding 1.

The discrete cube root problem

We first recall the discrete cube root problem and the corresponding classical hardness assumption. Let p and q be two distinct large prime numbers of the form 3k + 2 and 3k′ + 2, respectively, and assume that they are represented by approximately the same number of bits. Let N = pq be their product, which we take to be an n-bit integer, and denote by Z N = {0, . . . , N − 1} the integers modulo N. We consider the discrete cube function f N : ZN → ZN ,

f N ( y ) = y3

mod N.

(E1)

For primes of the form specified above, this map is a bijection on Z N , and therefore its inverse is well defined. We denote this inverse by

√ −1 gN (x) = f N (x) = 3 x

mod N.

(E2)

Indeed, since p = 3k + 2 and q = 3k′ + 2, one has gcd(3, p − 1) = gcd(3, q − 1) = 1. Hence cubing is invertible modulo both p and q, and therefore also modulo N = pq by the Chinese remainder theorem. Definition 12 (Discrete cube root problem). Let p and q be two distinct primes of the form 3k + 2 and 3k′ + 2, and let N = pq be an n-bit integer. Given as input N and x ∈ Z N , the discrete cube root problem asks to output y ∈ Z N such that y3 = x

mod N.

(E3)

36 The forward function f N can be computed efficiently on a classical computer using modular exponentiation. By contrast, computing the inverse function g N from only N and x is believed to be classically hard when the factorization of N is unknown. This is formalized by the following assumption, which is a special case of the RSA assumption with exponent e = 3. Definition 13 (Discrete cube root assumption). For any polynomial P(·), there is no classical algorithm A running in time P(n) such that, given as input an n-bit integer N = pq, where p and q are random primes of the form 3k + 2 and 3k′ + 2, and a random x ∈ Z N , outputs a value y ∈ Z N satisfying y3 = x

mod N

(E4)

with probability at least 1/P(n). The probability is taken over the random choice of the primes p, q, the input x ∈ Z N , and the internal randomness of A. The discrete cube function is also a trapdoor function. If the factors p and q of N are known, then one can compute Euler’s totient φ( N ) = ( p − 1)(q − 1). Since p = 3k + 2 and q = 3k′ + 2, both p − 1 and q − 1 are congruent to 1 modulo 3, and therefore gcd(3, φ( N )) = 1. It follows that 3 has a multiplicative inverse modulo φ( N ). Hence, by Euclid’s algorithm, one can efficiently find an integer d such that 3d = 1 mod φ( N ), or equivalently 3d = rφ( N ) + 1 for some integer r. We now show why this value of d gives a secret key for inverting the discrete cube function. Let x = f N (y) = y3 mod N for some y ∈ Z N . Then x d = (y3 )d = y3d = yrφ( N )+1

mod N.

(E5)

Using φ( N ) = ( p − 1)(q − 1) together with Fermat’s little theorem modulo p and modulo q, one obtains yrφ( N )+1 = y mod p and yrφ( N )+1 = y mod q. Since p and q are distinct primes, congruence modulo both p and q implies congruence modulo N = pq. Therefore, x d = y mod N. Thus, for every x ∈ Z N , the discrete cube root can be computed efficiently as g N ( x ) = x d mod N. The integer d therefore acts as a trapdoor key: it can be computed efficiently from the factorization of N, but is believed to be inaccessible to efficient classical algorithms given only N. Hence, knowledge of the factorization makes the inverse efficiently computable, whereas without this trapdoor the problem is assumed to be classically intractable.

Definition of the learning problem

Based on the DCR function in Eq. (E2), we define the following function that classifies the input x ∈ Z N of our learning problem. For a fixed shift parameter s ∈ Z N , define g N,s : Z N → {0, 1}, by  1, if g ( x ) ∈ s, s + N − 1 N 2 (E6) y( x ) = g N,s ( x ) =  0, otherwise where, as in Eq. 48, the interval is interpreted cyclically in Z N . Equivalently, using the key d as in Eq. (E5), N−1 . gdN,s ( x ) = 1 x d mod N ∈ s, s + 2

(E7)

Based on the labeling function above, we define the following learning problem. Let N be a large semiprime of the form N = pq as above, we formally define the concept class C N as the family of all concepts function g N,s in Eq. (E6):

C N = { g N,s | s ∈ Z N } .

(E8)

Then, the learner is given training data labeled by one of the concept in C N ,

T = {( xℓ , yℓ )}m ℓ=1 ,

xℓ ∈ Z N ,

yℓ = gdN,s⋆ ( xℓ ) ∈ {0, 1}, ⋆

(E9)

37 for some hidden target parameters (d⋆ , s⋆ ), where gdN,s denotes the representation of g N,s associated with the secret key d. Importantly, the learner does not know the true (d⋆ , s⋆ ), it can only retrieve information about them through the training data. We also assume that the learner has knowledge of N from the specification of the concept class C N or from the training data (for example we could trivially append it to all the training points xℓ ∈ Z N as ( xℓ , N, yℓ )). The goal of the learning task is to implement a learning algorithm A that meets the following learning condition. Definition 14 (Efficient PAC learning). Let C N be the DCR concept class. We say that C N is efficiently classically learnable if there exists a classical learning algorithm A with the following property. For every target concept g N,s ∈ C N , every probability distribution D over Z N , and every ε, δ > 0, given a training set T = {( xℓ , g N,s ( xℓ ))}m ℓ=1 with xℓ ∼ D drawn independently, the algorithm outputs a hypothesis h : Z N → {0, 1} such that Pr [h( x ) ̸= g N,s ( x )] ≤ ε

x ∼D

(E10)

with probability at least 1 − δ over the random draw of the training set and the internal randomness of A. Moreover, A runs in time polynomial in n, 1/ε, and 1/δ. 3.

Classical hardness of the learning problem

The classical hardness result for the DCR learning task relies on the result of [ACGS88]. In particular, they show that even a weak predictor for suitable bits of the discrete cube root is sufficient to recover the full inverse function. Based on this, [JGM+ 24] showed that no efficient classical learner can achieve accuracy significantly better than random guessing on the concept class defined from the discrete cube root function. We first recall the notion of an ϵ(n)-oracle. Definition 15 (ϵ(n)-oracle). Let O N,s be a probabilistic oracle such that, on input x ∈ Z N , it outputs the label g N,s ( x ) correctly with probability 1/2 + ϵ(n) over the random choice of x and the internal randomness of the oracle. We say that O N,s is an ϵ(n)-oracle for g N,s . The result of Alexi et al. [ACGS88] implies that a weak oracle with inverse-polynomial advantage over random guessing is already sufficient to break the discrete cube root assumption. Lemma 16 (Corollary (a) to Theorem 1 in [ACGS88]). For any concept g N,s in the DCR concept class, given access to a 1/ poly(n)-oracle √ O N,s for g N,s , there exists a polynomial-time classical algorithm that uses O N,s to compute the discrete cube root g N ( x ) = 3 x mod N for every x ∈ Z N with constant success probability.

Using this result, the authors in [JGM+ 24] obtain the following classical hardness statement for the corresponding learning problem. Lemma 17 (Classical hardness of the DCR learning task). Assume the discrete cube root assumption. Then no polynomialtime classical learning algorithm can, with constant success probability over the random generation of the training data and its internal randomness, output a hypothesis h whose expected error on the DCR concept g N,s is smaller than random guessing by an inverse-polynomial amount. More precisely, for the uniform distribution over Z N , no efficient classical learner can achieve Pr

x ∼U ( Z N )

[h( x ) ̸= g N,s ( x )] ≤

1 1 − . 2 poly(n)

(E11)

Indeed, if a classical learner satisfying the learning condition in Definition 14 existed, then the hypothesis it outputs would define a 1/ poly(n)-oracle for the label function g N,s . By Lemma 16, this oracle could then be used to compute the full discrete cube root function g N ( x ) in polynomial time with constant success probability. This would contradict the discrete cube root assumption. Therefore, under this assumption, the DCR learning task cannot be solved efficiently by classical learning algorithms with non-negligible advantage over random guessing. 4.

Quantum solution with Hamiltonian embedding

We now show how this problem can be encoded within the Hamiltonian-embedding framework. In particular, we construct an embedding Hamiltonian of the form HT = ∑ℓ H ( xℓ ) whose ground state can be used to correctly classify new inputs, in the same way as described in Section II. This allows us to show in the following lemma the quantum learnability of the DCR task within the Hamiltonian embedding model. Together with Lemma 17, this formally establishes the result anticipated in the main text as Theorem 4.

38 Lemma 18. Let CN be the DCR concept class, and let T = {( xℓ , yℓ )}m ℓ=1 be the training data associated to a concept ⋆ gdN,s⋆ ∈ CN with xℓ sampled independently from the uniform distribution U (Z N ). For every ϵ, δ > 0, there exists m = poly(n, 1/ϵ, log(1/δ)), an encoding unitary U ( x ), and a local label observable O(y), such that, with probability at least 1 − δ over the draw of T , the ground state of the corresponding embedding Hamiltonian HT induces a hypothesis f˜ : Z N → {0, 1} as in Eq. (5) satisfying i h ⋆ (E12) Pr f˜( x ) ̸= gdN,s⋆ ( x ) ≤ ϵ. x ∼U (Z N )

Moreover, such a ground state can be prepared efficiently by a quantum algorithm, and once obtained, the induced hypothesis f˜ can be efficiently evaluated classically. Proof. The embedded Hamiltonian is defined on the Hilbert space H = span{|d, s⟩ : d, s ∈ Z N } ⊗ C2 , where the last qubit serves as an ancilla with computational basis {|0⟩, |1⟩}. For each input x ∈ Z N , define the unitary U ( x ) |d, s⟩|b⟩ = |d, s⟩ |b ⊕ gdN,s ( x )⟩.

(E13)

We note that gdN,s ( x ) can be evaluated efficiently classically, since its computation only requires modular exponentiation and a classical comparison. Hence, there exists a polynomial-size classical circuit computing gdN,s ( x ) from x, d, and s. By reversibly implementing this computation, with at most polynomial overhead, one obtains a polynomialsize quantum circuit realizing the unitary U ( x ). Define now the observable O(y) = I ⊗ |y⟩⟨y|,

y ∈ {0, 1}.

(E14)

We can now define the embedding Hamiltonian HT for the training set T . For one training example ( xℓ , yℓ ), define the Hamiltonian term H ( xℓ , yℓ ) = I − U ( xℓ )† O(yℓ )U ( xℓ ) + I ⊗ |1⟩⟨1|.

(E15)

The dataset Hamiltonian HT is then: m

HT = ∑ H ( xℓ , yℓ ).

(E16)

ℓ=1

We now show that the ground state of HT can be used to correctly predict the labels of unseen inputs. In particular, in analogy to Eq. (5), for a new point xthe predicted label is the value of y that maximizes Tr O(y) U ( x ) |d⋆ , s⋆ , 0⟩⟨d⋆ , s⋆ , 0| U ( x )† , with |d⋆ , s⋆ , 0⟩ ground state of HT . We first show that the energy of each term H ( xℓ , yℓ ) is minimized by states of the form |d, s, 0⟩ with (d, s) such that ⋆ d g N,s ( xℓ ) = yℓ . Then, since the training labels are generated by the target concept yℓ = gdN,s⋆ ( xℓ ), the state |d⋆ , s⋆ , 0⟩ associated with the true parameters d⋆ and s⋆ minimizes every term in HT and thus has zero energy. Let us consider each term H ( xℓ , yℓ ) in Eq. (E15). We can expand U ( xℓ )† O(yℓ )U ( xℓ ) as follows: U ( x ℓ ) † O ( y ℓ )U ( x ℓ ) =

∑

(d,s): gdN,s ( xℓ )=yℓ

|d, s⟩⟨d, s| ⊗ |0⟩⟨0| +

∑

(d,s): gdN,s ( xℓ )̸=yℓ

|d, s⟩⟨d, s| ⊗ |1⟩⟨1|.

(E17)

Eq. (E17) follows from computing the action of U ( xℓ )† O(yℓ )U ( xℓ ) on generic states |d, s, b⟩ using the definitions of U ( x ) and O(y) in Eqs. (E13) and (E14): U † ( xℓ )O(yℓ )U ( xℓ ) |d, s, b⟩ = δy ,b⊕ gd ( x ) |d, s, b⟩ ℓ

N,s

ℓ

(E18)

Therefore, by rewriting the Hamiltonian H ( xℓ , yℓ ) in Eq. (E15) using Eq. (E17) we obtain: H ( xℓ , yℓ ) = I −

∑

(d,s): gdN,s ( xℓ )=yℓ

|d, s⟩⟨d, s| ⊗ |0⟩⟨0| −

∑

(d,s): gdN,s ( xℓ )̸=yℓ

|d, s⟩⟨d, s| ⊗ |1⟩⟨1| + I ⊗ |1⟩⟨1|

(E19)

We now evaluate the energy of a generic basis state |d, s, b⟩ to determine which states minimize the Hamiltonian terms H ( xℓ , yℓ ).

39 • If b = 1, i.e. |ψ⟩ = |d, s, 1⟩, then (

⟨ψ| H ( xℓ , yℓ )|ψ⟩ =

2, if gdN,s ( xℓ ) = yℓ , 1, if gdN,s ( xℓ ) ̸= yℓ .

(E20)

0, if gdN,s ( xℓ ) = yℓ , 1, if gdN,s ( xℓ ) ̸= yℓ .

(E21)

• If b = 0, i.e. |ψ⟩ = |d, s, 0⟩, then (

⟨ψ| H ( xℓ , yℓ )|ψ⟩ =

Hence, each term H ( xℓ , yℓ ) in the dataset Hamiltonian HT is minimized precisely by states of the form |d, s, 0⟩ satisfying gdN,s ( xℓ ) = yℓ . Since the total Hamiltonian is given by HT = ∑m ℓ=1 H ( xℓ , yℓ ), then for a generic basis state | d, s, 0⟩ one has i m h ⟨d, s, 0| HT |d, s, 0⟩ = ∑ 1 gdN,s ( xℓ ) ̸= yℓ .

(E22)

ℓ=1

Therefore, the energy of a state |d, s, 0⟩ is exactly the number of training examples misclassified by the hypothesis specified by (d, s):

⟨d, s, 0| HT |d, s, 0⟩ = {training points misclassified by gdN,s } .

(E23)

It follows that a basis state |d, s, 0⟩ has zero energy if and only if the corresponding hypothesis gdN,s correctly classifies ⋆ every training example. Since the target hypothesis gdN,s⋆ is consistent with the training data, the ground energy of HT is zero, and its ground space is spanned by the states |d, s, 0⟩ corresponding to hypotheses that are consistent with the training set. We note that the exponent d⋆ can be recovered efficiently by a quantum algorithm by factoring N [Sho97]. Once ⋆ d is known, the remaining parameter s can be learned efficiently from the training set by first computing the DCR √ ⋆ of each input: xℓ 7−→ xℓd mod N = 3 xℓ mod N = g N ( xℓ ), with g N the DCR function defined in Eq. (E2). The transformed dataset Te = {( g N ( xℓ ), yℓ )}m ℓ=1 then reduces the problem to learning the shift of a one-dimensional interval classifier. Since the data are realizable, a shift s′ that is consistent with all the training examples can be found efficiently. While several shifts may be consistent with a polynomial-size training set, sufficiently many training samples guarantee that a suitable consistent shift s′ has small population error under the uniform distribution U (Z N ). In particular, [JGM+ 24] shows that, with O(log(1/δ)/ϵ′ ) uniformly sampled examples, one can identify a shift s′ satisfying |s⋆ − s′ | ≤ ϵ′ N with probability at least 1 − δ, where |s⋆ − s′ | denotes the cyclic distance between the two shifts in Z N . Let us fix ϵ′ = ϵ/2. This in turn guarantees Pr g N,s′ ( x ) ̸= g N,s⋆ ( x ) ≤ ϵ. (E24) x ∼U (Z N )

Therefore, the quantum learning algorithm finds the state |d⋆ , s′ , 0⟩ which has zero energy with respect to the data Hamiltonian HT . Once a ground state |d⋆ , s′ , 0⟩ is obtained, a new point x is classified by the same labeling rule as in Eq. (5). In particular, the probability associated to the label y for an input x is given by i h Pd⋆ ,s′ ( x, y) = Tr O(y) U ( x ) |d⋆ , s′ , 0⟩⟨d⋆ , s′ , 0| U ( x )† . (E25) Since, by the construction of U ( x ) in Eq. (E13), U ( x )|d⋆ , s′ , 0⟩ = |d⋆ , s′ , gdN,s′ ( x )⟩, it follows that ( ⋆ 1, if gdN,s′ ( x ) = y, Pd⋆ ,s′ ( x, y) = ⋆ 0, if gdN,s′ ( x ) ̸= y. ⋆

(E26)

Hence, the prediction rule induced by the ground state |d⋆ , s′ , 0⟩ is precisely the hypothesis gdN,s′ . By Eq. (E24), this hypothesis has population error at most ϵ under U (Z N ), which proves the lemma. Finally, we note that, once the state |d⋆ , s′ , 0⟩ has been obtained quantumly, the prediction probability Pd⋆ ,s′ ( x, y) in Eq. (E25) can be evaluated entirely classically, since U ( x ) acting on |d⋆ , s′ , 0⟩ is efficiently computable by directly ⋆ evaluating gdN,s′ ( x ). ⋆

40 In Lemma 18, we constructed an embedding Hamiltonian HT whose ground space contains a state |d⋆ , s′ , 0⟩ as⋆ sociated with a hypothesis gdN,s′ that solves the DCR learning task, and we showed that such a state can be found efficiently by a quantum algorithm. In the following lemma, we show that, more generally, every ground state of HT induces a hypothesis with small population error on the DCR task. n o Lemma 19. Let H N = gdN,s : d, s ∈ Z N , be the hypothesis class of functions gdN,s as in Eq. (E7). Given the training set

d T = {( xℓ , yℓ )}m ℓ=1 , with yℓ = g N,s⋆ ( xℓ ), generated by drawing the inputs independently from the uniform distribution over Z N , let HT be the embedding Hamiltonian as in Eq. (E16). For any ϵ, δ > 0, if 1 1 2 log N + log , m≥ ϵ δ ⋆

then, with probability at least 1 − δ over the draw of the training set, every zero-energy basis state |d, s, 0⟩ of HT induces a hypothesis gdN,s satisfying h

x ∼U ( Z N )

i ⋆ gdN,s ( x ) ̸= gdN,s⋆ ( x ) ≤ ϵ.

In particular, the conclusion holds independently of whether d = d⋆ . Proof. By Eq. (E23), for every basis state |d, s, 0⟩, h i m ⟨d, s, 0| HT |d, s, 0⟩ = ∑ 1 gdN,s ( xℓ ) ̸= yℓ .

(E27)

ℓ=1

Since the target hypothesis gdN,s⋆ is consistent with the training data, the ground energy of HT is zero. Therefore, any zero-energy basis state |d, s, 0⟩ corresponds to a hypothesis gdN,s that is consistent with all training examples. Now consider a fixed hypothesis h ∈ H N whose population error with respect to the target hypothesis is larger than ϵ, namely h i ⋆ h( x ) ̸= gdN,s⋆ ( x ) > ϵ. (E28) R(h) := Pr ⋆

x ∼U ( Z N )

The probability that h nevertheless agrees with all m independently drawn training examples is

(1 − R(h))m ≤ (1 − ϵ)m ≤ e−mϵ . Since d, s ∈ Z N , the hypothesis class satisfies |H N | ≤ N 2 . Hence, by a union bound, h i b T (h) = 0 ≤ N 2 e−mϵ . Pr ∃h ∈ H N : R(h) > ϵ and R T

(E29)

(E30)

If 1 m≥ ϵ

1 2 log N + log δ

(E31)

then the right-hand side is at most δ. Therefore, with probability at least 1 − δ, every hypothesis in H N that is consistent with the training set has population error at most ϵ. Since every zero-energy basis state of HT corresponds to such a consistent hypothesis, the claimed result follows. We note that the same argument can be extended to arbitrary ground states that are superpositions of zero-energy basis states. Indeed, for a ground state |ψ⟩ = ∑d,s cd,s |d, s, 0⟩, the prediction probabilities are convex combinations of those associated with the individual hypotheses gdN,s , since the states |d, s⟩ are mutually orthogonal. Consequently, if every zero-energy basis hypothesis has population error at most ϵ/2, the classifier induced by |ψ⟩ has population

41 error at most ϵ, since a misclassification can occur only when at least half of the total weight is assigned to basis hypotheses that err on that input. We note that this example establishes the existence of a hard task that can be solved within the Hamiltonian embedding framework, for which a quantum algorithm can efficiently prepare the ground state of the embedded Hamiltonian, while doing so classically is believed to be hard. This is sufficient to demonstrate an advantage for that specific task. For more general learning problems, however, we do not expect to be able to design a problem-specific algorithm for preparing the ground state of the embedded Hamiltonian. Instead, one would likely rely on more generic quantum ground-state preparation methods, such as QPE or SKQD. Finally, for the ground state |d⋆ , s′ , 0⟩ constructed by the quantum learning algorithm, inference can be performed ⋆ entirely classically. Indeed, once the parameters d⋆ and s′ are known, the corresponding hypothesis gdN,s′ ( x ) can be evaluated efficiently using classical modular exponentiation and a classical comparison. This differs from the settings considered in Section III B, where the classical surrogates are constructed using task-agnostic procedures. Appendix F: A construction with highly entangled ground states

In Section II A 1, we discussed that the Hamiltonian embedding framework allows for different choices of the Hermitian operator A entering the observable O(y) in Eq. (9). The choice of the label observable O(y) can significantly affect the performance of the model, and different operators may be more suitable depending on the learning task under consideration. Although identifying which choices of O(y) are best suited to a given learning problem is itself an interesting question, in this Appendix and Appendix G we show that the choice of O(y) can influence properties of the ground state ψT ,α⋆ found by the Hamiltonian embedding model, such as its degree of entanglement. In this Appendix, we provide a simple example of a label dependent observable O(y) for which the ground states obtained within the Hamiltonian embedding framework exhibit a high degree of entanglement. We note that, in this particular construction, the resulting states remain efficiently classically representable, as they are derived from Clifford states. More elaborate choices of O(y) could in principle give rise to ground states with more complex structure. The purpose of this example is therefore not to establish classical hardness results as we did in Section III C, but rather to illustrate how the choice of O(y) can directly influence the properties of the ground state found by the model. 1.

Entangled graph states as ground states

In order to obtain entangled ground states, we start from the parent Hamiltonian of a graph state, which can exhibit a large amount of bipartite entanglement. Specifically, for n qubits, we consider a graph composed of two cycles of n/2 vertices each, connected by a perfect matching between corresponding vertices. To each vertex a we associate the graph-state stabilizer Ka = Xa

∏ Zq ,

(F1)

q∈ N ( a)

where N ( a) denotes the set of vertices adjacent to a. The corresponding parent Hamiltonian is constructed by introducing one local check for each qubit h a = 12 ( I − Ka ) , and defining H = n1 ∑na=1 h a . Since Ka has eigenvalues ±1, each check contributes energy 0 when Ka = +1 and energy 1 when Ka = −1. The ground state is therefore the unique simultaneous +1 eigenstate of all the stabilizers Ka , namely the graph state associated with the chosen graph. In order to parameterize this construction with respect to the label y, we associate the two classes y = 0 and y = 1 with two different stabilizer syndromes. In particular, we choose the two syndromes β(0) = (0, . . . , 0), β(1) = (1, 0, 1, 0, . . . , 1, 0) with both β(0) , β(1) ∈ {0, 1}n and define the label-dependent checks as (y) 1 βa h a (y) = I − (−1) Ka , 2

(F2)

We can define then the label-dependent observable as the parent Hamiltonian of the graph states with the signature

β(0) or β(1) according to the label y:

" # (y) 1 n 1 1 n βa O(y) = ∑ h a (y) = I − ∑ (−1) Ka . n a =1 2 n a =1

(F3)

42 (1)

For y = 0, the ground state satisfies Ka = +1 for every qubit a, whereas for y = 1 it satisfies Ka = (−1) β a , corresponding to the alternating stabilizer syndrome. To complete the construction of the Hamiltonian embedding model, we now dress the observable O(y) with a particular choice of encoding unitary U ( x, α). In particular, for n qubits we choose the simple unitary U ( x, α) =

n O a =1

Y RZ αZ a x a RY α a x a .

(F4)

n 2n denoting the set of trainable parameters. We notice that there are no entangling gates with α = {αYa , α Z a } a =1 ∈ R in U ( x, α). This, however, does not preclude the embedding Hamiltonian from having entangled ground states, as we demonstrate below. For a datapoint ( x, y) the encoding Hamiltonian is then (y)

ea ( x, α) 1 n I − (−1) β a K , Hx,y (α) = ∑ n a =1 2

(F5)

ea ( x, α) = U ( x, α)† Ka U ( x, α). For a training set T = {( xi , yi )} N , the total Hamiltonian is HT (α) = where K i =1 N 1 H ( α ) . ∑ x ,y i i N i =1 2.

Entanglement properties of the ground state

Our numerical experiments show that, for this construction, the ground states learned by the model exhibit a high degree of entanglement. This can be understood directly from the structure of the label Hamiltonians. The E ground state Gβ(y)

of the labeled observable O(y) in Eq. (F3) is the simultaneous eigenstate satisfying Ka Gβ =

(−1) β a Gβ , ∀ a. For β(0) , this is the graph state | G0 ⟩ =

∏ CZij |+⟩⊗n ,

(F6)

(i,j)∈ E

while for β(1) we have 



| G1 ⟩ = 

∏

(1) a:β a =1

Za  | G0 ⟩ .

(F7)

The embedding Hamiltonian associated with the full training set is HT (α) =

1 H . |T | (x ,y∑)∈T (xℓ ,yℓ ) ℓ

(F8)

ℓ

Its ground state ψT⋆ ,α is in general determined by the competition between the different training terms and therefore does not coincide with the ground state of any individual term. The case of zero energy loss is achieved when the learned parameters α⋆ and ground state |ψT ,α⋆ ⟩ satisfy U ( xℓ , α⋆ ) |ψT ,α⋆ ⟩ = Gyℓ ,

(F9)

for every training example ( xℓ , yℓ ). Since U ( x, α) is a product of single-qubit rotations, it is a local unitary and therefore preserves the bipartite entanglement spectrum. Hence, we can expect that the learned ground state |ψT ,α⋆ ⟩ follows the same bipartite entanglement as the corresponding graph states | G0 ⟩ and | G1 ⟩. 3.

Numerical demonstration

To study the entanglement properties of the ground states reached by the constructed HE model, we consider the following constructed classification task. As input data x, we use images from the Fashion-MNIST dataset and apply PCA to reduce their dimensionality to n, matching the number of qubits in the model. Rather than using the

43 original Fashion-MNIST class labels, we assign labels through a planted teacher model. Specifically, we fix a highly entangled graph-state ground state |ψT ⟩ together with a set of teacher parameters αT , and we label each input x according to the class y ∈ {0, 1} whose corresponding Hamiltonian term Hx,y (α T ) in Eq. (F5) has the lower energy, i.e. y( x ) = arg min ⟨ψT | Hx,y (α T ) |ψT ⟩ y∈{0,1}

(F10)

We note that, in this Appendix, O(y) is defined as the label-dependent parent Hamiltonian, so that the state associated with the correct label minimizes its expectation value and prediction is performed via an arg miny rule. This differs only by convention from the main text, where O(y) denotes the complementary score observable and prediction is therefore performed via an arg maxy rule. Moreover, the operators used here do not satisfy O(0) + O(1) = I. Prediction selects the label with the lower expectation value, without interpreting these values as probabilities. The resulting labeled dataset is then used to train a Hamiltonian-embedding student model of the same form, without providing it with either the teacher state or the teacher parameters. We notice that the purpose of this experiment is to study the entanglement properties of ground states that can be learned within the Hamiltonian-embedding model, rather than the ability of the outer optimization to locate such solutions from arbitrary initializations. We therefore initialize the student parameters in the neighborhood of the planted teacher, α(0) = α T + ϵ, with ϵ j ∼ N (0, 0.252 ). In particular, we choose the teacher parameters α T by sampling each entry independently from a uniform distribution over [−π, π ], and take as the teacher state |ψT ⟩ a highly entangled graph-syndrome state,

|ψT ⟩ = Z β T | G0 ⟩ ,

(F11)

where | G0 ⟩ is the graph state introduced in Eq. (F6) and β T ∈ {0, 1}n is a fixed nontrivial stabilizer syndrome for each system size. We use N = 800 training samples, together with 200 validation samples and 400 test samples for each run. For the outer parameter optimization, we use the energy-based analogue of the hinge loss introduced in Eq. (15). Since classification is performed by selecting the label with the lowest energy, we define the class-energy margin as (γ)

∆α ( x, y) = E1−y ( x ) − Ey ( x ), and optimize ℓhinge ( x, y) = max{0, γ − ∆α ( x, y)}, with γ = 0.05. At each optimization step, the ground state of the corresponding training Hamiltonian is recomputed exactly and then kept fixed while updating the encoding parameters. After training, we study how accurately the final exact learned ground state can be represented by an MPS of bounded bond dimension. We fix the qubit ordering used in the simulation and convert the exact state into an MPS, truncating its Schmidt decompositions so that the maximum bond dimension is bounded by χ. For each value of χ, we evaluate both the fidelity of the resulting truncated MPS with the exact learned state and the test accuracy obtained when this MPS is used in place of the exact state. As shown in Fig. 2, low-bond-dimension MPS approximations in this fixed ordering have low fidelity with the learned ground state, and the corresponding classification accuracy is substantially reduced. Increasing χ progressively recovers both the exact state and the classification performance. Appendix G: A provable separation between Pauli surrogates and Hamiltonian embedding

In Section III B, we showed that, for certain choices of the encoding unitary U ( x, α), the trained Hamiltonian embedding model can be evaluated classically once a suitable classical surrogate has been constructed from the final state ψT ,α⋆ . This naturally raises the question of whether such a surrogate could instead be trained directly, thereby removing the need to obtain the quantum state altogether. In this Appendix, we present an example in which we analytically show that the physical constraints imposed by the Hamiltonian embedding model lead to better predictions on unseen data than the corresponding Pauli classical surrogate. Moreover, the efficiency guarantee for the classical surrogate based on reduced density matrices introduced in Eq. (47) does not extend to this construction. In order to prove a separation between the two models, we use properties of stabilizer codes and stabilizer states. Similar to Appendix F, we consider a stabilizer code with stabilizers Ka , a ∈ [n]. Concretely, in our example, we have a graph G with stabilizers of the form Ka = Xa ∏q∈ NG (a) Za , where NG ( a) is the set of neighboring vertices of a in the graph. This graph and the set of stabilizers is known and specifies the “concept class” of the labeling functions. In particular, each concept is associated with a distinct assignment specifying whether each Ka acts as a stabilizer or an antistabilizer of the ground state. Thus, there are 2n possible functions for a given graph, indexed by the syndrome vector β whose coordinates β a are +1 (-1) depending on whether Ka is a stabilizer (antistabilizer). We now describe the function more explicitly. Let the data points have coordinates xi = | x |i ζ i with | x |i ∈ [0, 1] and ζ i = sign( xi ). For an input vector x, we form a query by considering which coordinates i have a magnitude above

44 (b) 1.0

0.8

0.9

Test accuracy

Fidelity

(a) 1.0

0.6

0.4

0.8

0.7

0.6 0.2 0.5 0.0 21

Bond cap χ n=6

n=8

Bond cap χ n = 10

n = 12

n = 14

n = 16

n = 18

n = 20

Figure 2. Left: fidelity between the final exact learned ground state and its MPS approximation obtained by Schmidt truncation at maximum bond dimension χ, using the fixed qubit ordering of the simulation. Right: test accuracy obtained by replacing the exact learned state with the corresponding truncated MPS state.

some cutoff ξ. Define this to be the “active set” S = {i : | xi | ≥ ξ }. We imagine the active set S to be a query to the function since the function depends only on the signs {ζ i , i ∈ S} (and not on | xi |). For each such input, construct an “error vector” Wζ := ∏i∈S:ζ i =−1 Zi , supported on the collection of places where the query has a negative sign. Now define η (ζ, a) to be +1 if the stabilizer Ka commutes with the error Wζ and −1 otherwise. The function to be learned is: y β ( x ) = ∏ y β( x),a ,

where

(G1)

a∈S

y β( x),a = β a η (ζ, a). In our case where the stabilizers are of the form Ka = Xa ∏q∈ NG (a) Zq , whether an error vector Wζ commutes with Ka or not is decided by whether ζ a = ±1. In other words, we have Wζ† Ka Wζ = ζ a Ka ,

(G2)

so in this case η (ζ, a) = ζ a . Therefore, the function is yβ (x) = 

=

∏ βa

a∈S

∏

a:| x | a ≥ξ

∏ ζa

a∈S

! (G3)



β a · sign( x a ) ,

(G4)

where we have made the dependence on x explicit by writing out S and ζ a in terms of x. We define some more operators for notational convenience. Let L a = −iZa Ka = Ya ∏q∈ NG (a) Zq . This operator anticommutes both with Za and Ka . Also, define LS = ∏ a∈S L a and similarly ZS and KS . Note that there is no ambiguity with respect to ordering in these definitions because all the different L a operators (respectively, Za and Ka ) commute with each other. For a given active set S, let RS be a Clifford operator that applies the map R†S Z1 RS = LS .

(G5)

Now consider the unitary U ( x ) = RS ∏ a∈S exp[−iα a x a Za /2]. The form of U ( x ) matches the definition in this work of a unitary that is efficiently constructible from the input x. Moreover, the gates that depend on the magnitudes of x

45 have the same form as those considered inhthe numerical i sections of this work, namely Pauli rotations whose angles depend linearly on the input features, exp − 2i (⃗α · ⃗x ) P , where P denotes the corresponding Pauli generator. Now U ( x )† Z1 U ( x ) = ∏ exp[+iα a x a Za /2] R†S Z1 RS ∏ exp[−iαb xb Zb /2] a∈S

b∈S

= ∏ exp[+iα a x a Za /2] LS ∏ exp[−iαb xb Zb /2]. a∈S

(G6) (G7)

b∈S

Since Za commutes with all Lb except L a , we get αa xa αa xa αa xa αa xa U ( x )† Z1 U ( x ) = ∏ cos I + i sin Za L a cos I − i sin Za 2 2 2 2 a∈S

(G8)

= ∏ (cos (α a x a ) L a + i sin (α a x a ) Za L a )

(G9)

= ∏ (cos (α a x a ) L a + sin (α a x a )Ka ).

(G10)

a∈S

Recall x a = | x a |ζ a , ζ a = ±1. Then

U ( x )† Z1 U ( x ) = ∏ (cos(α a | x a |) L a + ζ a sin(α a | x a |)Ka ) .

(G11)

a∈S

We note that this construction also illustrates the limitation of the reduced-density-matrix surrogate in Eq. (47) in Section III B 4. Here, the readout is A = Z1 , and the support of the dressed observable U ( x, α)† Z1 U ( x, α) depends on the active set S( x ) and contains its qubits. Consequently, a fixed region supporting all training observables must SN contain ℓ= 1 S ( xℓ ), the set of qubits active in at least one training example. Whenever the sampled active sets cover all n qubits, the fixed supporting region for all the training observables U ( x )† Z1 U ( x ) must contain the entire register. Therefore, the corresponding density matrix σ1 has dimension 2n . Thus, the polynomial-time argument based on explicitly diagonalizing matrices on small fixed regions does not apply. We remark that the accuracy separation proved below concerns the Pauli surrogate introduced in Eq. (46). 1.

Idealized embedding Hamiltonian and its ground state: perfect classification

The Hamiltonian term corresponding to this data point is I − yU † Z1 U 2 I − y ∏ a (cos (α a | x a |) L a + ζ a sin (α a | x a |)Ka ) = 2 I − ∏ a β a ζ a (cos (α a | x a |) L a + ζ a sin (α a | x a |)Ka ) = . 2

H ( x, y) =

Rewrite this in terms of binary variables z ∈ {0, 1}|S| as  1  H ( x, y) = × I − 2

∏ βa

a∈S

(G12) (G13) (G14)



∑

Pz  ,

(G15)

z∈{0,1}|S|

where Pz = ∏ a:za =0 ζ a cos (α a | x a |) L a × ∏ a:za =1 sin (α a | x a |)Ka . If we average this over all data points from the distribution, then considering the different signs of x a at the same magnitude, i.e. averaging over ζ a = ±1 while keeping | x a | fixed: we get Pz + Pzā = 0

if z a = 0,

(G16)

where z ā is the binary vector z with the bit at position a flipped: z āa = 1 − z a . The average over signs gives

∑

z∈{0,1}n

Pz = ∏ sin (α a | x a |)Ka . a∈S

(G17)

46 This means the embedding Hamiltonian averaged over different signs for the same | x a | gives " ! # 1 HD (| x |) = × I − ∏ β a ∏ sin (α a | x a |)Ka 2 a∈S a∈S I − ∏ a∈S ( β a sin (α a | x a |)Ka ) . 2

(G18) (G19)

For the full embedding Hamiltonian, we should also average over the distribution of x for all coordinates, denoted g ( x ): HD =

dn xg( x )

I − ∏ a∈S β a (sin (α a | x a |)Ka ) . 2

(G20)

This is a commuting Hamiltonian whose terms are products of β a Ka . The number of products that appear is precisely |S|, which is the sparsity of the input vector x above the cutoff ξ. If the distribution over the input data is supported on 1-sparse inputs x, then these stabilizers β a Ka are represented in the Hamiltonian. Since these stabilizers commute with each other, the full embedding Hamiltonian is frustration-free. Its ground space is the +1 eigenspace of the operators ∪S {∏ a∈S β a Ka }. Note that if the distribution that the input data x is drawn from is not supported on all coordinates, the group generated by the stabilizer generators ∪S {∏ a∈S β a Ka } may not be full rank. Now consider classification. Assume that the ground state |ψ⟩ of the embedding Hamiltonian is found by the ground state solver (which we treat as a black box in this Appendix). The expectation value of U ( x )† Z1 U ( x ) in this state is:

⟨ψ| U ( x )† Z1 U ( x ) |ψ⟩ = ⟨ψ| ∏ (cos (α a | x a |) L a + ζ a sin (α a | x a |)Ka ) |ψ⟩ .

(G21)

a∈S

We have ⟨ψ| L a |ψ⟩ = 0 because L a anticommutes with β a Ka , which is the stabilizer of |ψ⟩. Therefore,

⟨ψ| U ( x )† Z1 U ( x ) |ψ⟩ = ⟨ψ| ∏ (cos (α a | x a |) L a + ζ a sin (α a | x a |)Ka ) |ψ⟩ = ∏ β a ζ a sin (α a | x a |). a∈S

(G22)

a∈S

The sign of this is ∏ a∈S β a ζ a sign(α a ), which matches the function y β ( x ) whenever sign(α a ) = 1 for all a ∈ S. Therefore, classification is perfect as long as α a > 0. 2.

Embedding Hamiltonian with finite training data

In this subsection, we specify the distribution over inputs and prove that there exists a regime of number of training points such that the empirical Hamiltonian embedding model classifies the data nearly perfectly; in this same regime we can also prove an upper bound on the accuracy of the classical Pauli surrogates described in Section III B 2. This separation between the classical surrogate and the Hamiltonian embedding model proves that this model generalizes better when the model is suited to the task. We choose the input training data to be drawn from the following distribution. First, we choose the size of the active set according to the distribution Pr(|S| = 2) =

4 , 7

Pr(|S| = 3) =

2 , 7

Pr(|S| = 4) =

1 . 7

(G23)

Then, for a given size |S|, we pick the active coordinates uniformly. Let a denote a location that is active, that is a ∈ S. The input vector x at coordinate a has magnitude | x | a = r a , where r a is a random number in [ξ, 1]. This number r a is fixed for the duration of the experiment, that is, any other training point with a an active location will have | x | a = r a . Now consider the empirical embedding Hamiltonian as opposed to the idealized one in Eq. (G20). The difference is that in practice we may not have queried every possible input, which is precisely the case where we expect the Hamiltonian embedding model to generalize better by virtue of its inductive bias. The empirical embedding Hamiltonian is   ! 1 1   HT (α) = H ( xℓ , yℓ ) = ∑ ×  I − ∏ βa (G24) ∑ ∑ Pzℓ (α) , N ℓ 2N ℓ ℓ a∈S Sℓ | | z∈{0,1}

47 where Pzℓ (α) = ∏ a:zℓ =0 ζ aℓ cos (α a x ℓa ) L a × ∏ a:zℓ =1 sin (α a x ℓa )Ka and ℓ = 1, 2, . . . N denotes the different input traina

ing data points. Since, as mentioned before, for different data points we have the same magnitude x ℓa = r a ∀ℓ, we can rewrite this as: Pzℓ (α) =

∏ ζ aℓ cos (αa ra ) La × ∏ sin (αa ra )Ka .

a:zℓa =0

(G25)

a:zℓa =1

Notice that in this case, since we do not assume that the different signs ζ a queried at the same magnitude r a cancel out, and thus the cosine term survives. This can potentially be problematic because L a anticommutes with Ka and would make the Hamiltonian frustrated. Nevertheless, the empirical embedding Hamiltonian can have the same ground state |ψ⟩ as the full embedding Hamiltonian if the angles α can be trained. In this case, setting α∗a = 2rπa ensures that the cosine term gives 0 and thus only the all-sine term zℓ = 1 . . . 1 survives: Pzℓ (α) = ∏ a∈Sℓ Ka , giving " !# 1 ∗ × I − ∏ β a Ka . (G26) HT (α ) = ∑ 2N ℓ ℓ a∈S

As we have seen before, the Hamiltonian has as ground space the +1 eigenspace of the stabilizer group ∪ℓ {∏ a∈Sℓ β a Ka }, with energy ⟨ HT ⟩ = 0. As we have restricted ourselves to positive semidefinite Hamiltonians, this proves that α∗ is an optimum. a.

Empirical accuracy of Hamiltonian embedding

Consider the binary matrix A, whose ℓth row records the active set Sℓ . We are interested in the event that A has rank n over F2 . This is equivalent to whether the sampled stabilizer products generate the full stabilizer group ⟨ β a Ka : a ∈ [n]⟩. In this case, the empirical embedding Hamiltonian would have the unique ground state |ψ⟩. Note that merely requiring every coordinate a to appear in some sampled query is not sufficient for this spanning condition. On the other hand, if N is large enough that every possible input in the support of the distribution is observed, the classical surrogate can also classify perfectly by memorization. Therefore, we now investigate the case where N is bounded. Specifically, we choose N = 4n log n. Our analysis of accuracy depends on the probability that the stabilizer group ∪ℓ {∏ a∈Sℓ β a Ka } has full rank. Consider the matrix A N ×n whose rows correspond to training data. For each data point ℓ, we have Aℓ,a = 1 if a is in the active set Sℓ , and 0 otherwise. Now A is not full rank if and only if there exists a nonzero binary vector v such that A · v = 0. Denote by T the set of locations i where vi = 1 and let | T | = t. Then Aℓ · v = 0 ⇔ Sℓ ∩ T = 0 mod 2.

(G27)

Suppose we fix a v and consider whether for a random query ℓ from the distribution v is “detected” by ℓ, i.e. whether |S ∩ T | = 1 mod 2. The probability of this event, conditioned on |S| = 2, is: t × (n − t) Pr |S ∩ T | = 1 |S| = 2 = . (G28) (n2 ) Similarly, conditioning on |S| = 3, we get the two cases |S ∩ T | = 1 or 3, giving t t × (n− ( t ) × (n−t) 2 ) Pr |S ∩ T | = 1 mod 2 |S| = 3 = + 3 n 0 . n (3) (3) Lastly, when |S| = 4, we again have the possibilities |S ∩ T | = 1 or 3, with the probability t t × (n− ( t ) × (n−t) 3 ) Pr |S ∩ T | = 1 mod 2 |S| = 4 = + 3 n 1 . n (4) (4)

(G29)

(G30)

Denote the combined probability 4 2 Pr |S ∩ T | = 1 |S| = 2 + Pr |S ∩ T | = 1 mod 2 |S| = 3 7 7 1 + Pr |S ∩ T | = 1 mod 2 |S| = 4 , 7

rt = Pr [|S ∩ T | = 1 mod 2] =

(G31)

48 where we have made use of the distribution over the active set sizes. The probability that no row detects v is given by (1 − rt ) N ≤ e− Nrt . By a union bound over all nonzero v ∈ F2n , n n − Nrt Pr[rank( A) < n] ≤ ∑ . (G32) e t t =1 We set N = ⌈4n log n⌉ and divide the sum into two regimes. First, suppose 1 ≤ t ≤ n/2. Retaining only the contribution from 2-sparse queries gives rt ≥

= Consequently,

4 t(n − t) 7 (n2 )

(G33)

8t(n − t) 4t ≥ . 7n(n − 1) 7n

(G34)

16 t log n. 7

(G35)

Nrt ≥

Using (nt) ≤ (en/t)t , the contribution of this regime to Eq. (G32) is at most ⌊n/2⌋

∑

t =1

⌊n/2⌋ n − Nrt en t −16t/7 n e ≤ ∑ t t t =1 e t −9/7 ≤∑ = O n . 9/7 t ≥1 n

(G36) (G37)

Next, suppose n/2 < t ≤ n. A 3-sparse query lying entirely inside T has odd intersection with T, and hence t

rt ≥

1 2 (3) ≥ 7 (n3 ) 56

(G38)

for all sufficiently large n. The contribution of this regime is therefore at most n log n n − N/56 = o (1). 2 e ≤ exp n log 2 − 14

(G39)

Combining the two regimes gives Pr[rank( A) < n] = o (1),

(G40)

so A has rank n with probability 1 − o (1). This means that the accuracy of the Hamiltonian embedding model, assuming that the ground-state solver is perfect, would be 1 − o (1). b.

Accuracy of classical surrogates

For an active set S, define PS,z := cS,z ( x, α) :=

∏

cos(α a | x a |)

a∈S:z a =0 a∈S:z a =0

∏

a∈S:z a =1

Ka ,

(G41)

∏

a∈S:z a =1

ζ a sin(α a | x a |),

(G42)

where z ∈ {0, 1}|S| . The inference observable and its Pauli surrogate are then U ( x, α)† Z1 U ( x, α) = gS ( x, α) =

∑

cS,z ( x, α) PS,z ,

(G43)

∑

cS,z ( x, α)mS,z .

(G44)

z∈{0,1}|S| z∈{0,1}|S|

49 Here, the expectation value of PS,z , namely ⟨ϕ| PS,z |ϕ⟩, is replaced by the coefficient mS,z , which is unconstrained. For the graph stabilizers considered here, PS,z has an X or Y on precisely the qubits in S. Pauli observables belonging to different active sets are therefore distinct, so the free surrogate has a separate coefficient block for each possible query S. For fixed α, the loss function used in ridge regression is 2



J (m; α) =

1 N  ∑ yℓ − 2N ℓ= 1

λ  cSℓ ,z ( x ℓ , α)mSℓ ,z  + ∑ ∑ m2S,z , 2 S∈Q z∈{0,1}|S| Sℓ z∈{0,1}| |

∑

(G45)

where the second sum S ∈ Q is over all supports Q = {S ⊆ [n]}. If a support S is absent from the training set, its coefficient block does not enter the data-fitting term. The unique ridge solution is therefore m⋆S,z (α) = 0

for all z,

S∈ / QT ,

(G46)

where QT is the set of observed supports. This holds for every α, and hence also when α and m are optimized jointly for the classical surrogate. We can evaluate the coefficients explicitly when α is fixed to the Hamiltonian-embedding optimum constructed above, α⋆a r a = π2 . Writing β S := ∏ a∈S β a , we then have cS,1S ( x, α⋆ ) = ∏ ζ a ,

cS,z ( x, α⋆ ) = 0

a∈S

( z ̸ = 1S ),

(G47)

and the label is y β ( x ) = β S ∏ a∈S ζ a . If NS denotes the number of training examples with support S, the ridge solution is found by setting ∂J (m; α) =0 mS,z

(G48)

=−

NS N β δz,1 + S mS,z δz,1S + λmS,z , N S S N

(G49)

giving m⋆S,z (α⋆ ) =

 β 

NS , NS + λN

NS > 0 and z = 1S , otherwise.

(G50)

Consequently, every observed support is classified correctly: y β ( x ) gS ( x, α⋆ ) =

NS > 0, NS + λN

NS > 0,

(G51)

whereas an unobserved support has zero score. Let Mk (T ) be the number of distinct k-element supports observed during training. A new query picked at random according to the distribution in Eq. (G23), will have been seen previously with probability 4

qseen (T ) = ∑ pk k =2

Mk (T ) . (nk)

(G52)

On a previously unobserved support the accuracy is 1/2. Therefore, for any fixed or independently optimized α, the ridge classifier’s accuracy is Aridge (T ; α) =

1 − qseen (T ) + qseen (T ) Aseen (T ; α), 2

(G53)

where Aseen is the classifier’s accuracy on previously seen points. This gives Aridge (T ; α) −

1 1 ≤ qseen (T ). 2 2

(G54)

50 When α is the optimal Hamiltonian-embedding solution, Eq. (G50) gives Aseen = 1, and hence Aridge (T ; α⋆ ) =

1 1 + qseen (T ). 2 2

(G55)

Thus even if the surrogate optimizes α independently its accuracy remains bounded; moreover, fixing α = α⋆ shows that the bound is achievable. Lastly, we bound qseen . For a particular k-element support, the probability of appearing at least once in N independent training queries is pk N 1− 1− n . (G56) (k) Averaging Eq. (G53) over the training set gives " # 1 1 4 pk N ET [ Aridge ] ≤ + ∑ pk 1 − 1 − n 2 2 k =2 (k) 1 1 4 N p2k + n 2 2 k∑ =2 ( k ) 4 1 1 N 16 + + , ≤ + 2 98 (n2 ) (n3 ) (n4 )

≤

where the first inequality is an equality at α⋆ . In particular, for N = Θ(n log n), 1 log n Aridge = + O . 2 n

(G57) (G58) (G59)

(G60)

The difference between the surrogate and the Hamiltonian embedding model is that the latter combines the observed stabilizer products: once their incidence vectors (i.e. the rows in A) have full rank, the underlying generator signs β a are fixed and every support can be classified. The fact that a single state is consistent with all Pauli expectation values and also informs the expectation values of unseen Paulis by virtue of the Pauli algebra is what provides the inductive bias in this case, allowing Hamiltonian embedding to generalize better than the classical surrogate. 3.

Required bond dimension for good classification with Hamiltonian embedding

In this subsection, we empirically study this previous example by instantiating a family of graphs and their corresponding codes. We randomly generate a simple, connected 3-regular graph, which fixes the stabilizer code family, at each system size. As mentioned earlier, the stabilizers are β i Xi ∏ j∈ N (i) Zj , with the signs β ∈ ±1 that describe the function also drawn uniformly at random at every vertex. We set the cutoff ξ = 0.6 and draw the radius of every active coordinate | xi | ∼ Unif[0.6, 1] for every realization and fix them thereafter. The signs of the input coordinates are again chosen uniformly at random. This fixes the distribution over input vectors ∈ Rn . We take the same regime of N = 4n log n, with the same distribution over queries as Eq. (G23). In order to train the Hamiltonian embedding model, we adopt a different numerical technique from the one considered in Section IV. Specifically, we apply imaginary time evolution on a classical MPS representation of the state. The MPS has a fixed ordering and no attempt was made at optimizing the ordering in order to better represent the ideal stabilizer ground state, since this strategy does not scale in general. After each application of imaginary time evolution, we compress the MPS up to a maximum bond dimension χ, with a numerical cutoff of 10−8 for the singular value decomposition (SVD). For the numerical optimization, we parameterize the rotation angles as α a = αc + ∆α u a , where αc =

π (1 + ξ −1 ), 4

u a ∈ [−1, 1], ∆α =

π −1 ( ξ − 1). 4

(G61)

(G62)

π Thus α a ∈ [ π2 , 2ξ ], an interval containing the analytic optimum α⋆a = 2rπa for r a ∈ [ξ, 1]. We initialize u a ∼ Unif[−0.75, 0.75] and optimize u using Adam, projecting it onto [−1, 1]n after each update. During this procedure, we optimize the energy of the parameterized Hamiltonian HT (α), keeping the state fixed.

51 (b) 1.0

0.9

0.8

Test A0.01

Test A0.1

(a) 1.0

0.7 0.6

0.5

0.5 21

210

211

Bond cap χ

210

211

210

211

Bond cap χ

(c)

(d)

0.5

1.0

Mean signed test margin

Training energy

0.4 0.3 0.2 0.1 0.0

0.8 0.6 0.4 0.2 0.0

210

211

Bond cap χ

n = 12

n = 24

n = 36

n = 48

n = 60

Figure 3. Test accuracy Aη as a function of bond dimension at system sizes n ∈ {12, 24, 36, 48, 60}, averaged over 10 trials, with (a) η = 0.1, (b) η = 0.01. (c) Minimum energy seen during training as a function of bond dimension. (d) The signed margin y( x ) U ( x )† Z1 U ( x ) .

We illustrate the accuracy of the trained model as a function of the maximum bond dimension cap χ in Fig. 3. Writing m( x ) = y( x )⟨U ( x )† Z1 U ( x )⟩, each test example receives score 1 when m( x ) ≥ η, 0 when m( x ) ≤ −η, and 1/2 when |m( x )| < η. The accuracy Aη is the average of these scores. We empirically find that increasing bond dimension is necessary in order to achieve the same accuracy at larger system sizes. Similarly, with increasing bond dimension, the energy and the observable’s margins also improve. In Fig. 4, we also plot the empirical accuracy of the classical surrogates and compare with the one from Hamiltonian embedding. For the points labeled as coming from HE, we show the highest test score observed over all bond dimensions. We notice that the classical surrogates trained on the ridge loss show a gradual decay in accuracy with respect to system size n.

Appendix H: Numerical setup

In this Appendix, we discuss the Hamiltonian embedding model used in all the experiments of Section IV, together with the classical simulation of the training loop. The protocols of the individual tasks are given in Appendices I and J.

Encoding circuit and embedding Hamiltonian

We use the encoding unitary of Eq. (12) with K = 3 layers. Layer ℓ applies to every qubit j the sequence of single-qubit Pauli rotations e e

ℓ Z Z −iθ j,4 j j+1 /2

ℓ Z /2 −iθ ℓ Y /2 −iθ ℓ Z /2 −iθ j,3 j j,2 j j,1 j

, followed by a ring of nearest-neighbor two-qubit rotations

with periodic boundary conditions. All 12n angles are independently trainable, and each is a linear

52 1.0

Test A0.01

0.9

0.8

0.7

0.6

0.5 12

System size n HE, best over χ ≤ 2048 (each trial) Ridge, HE-trained α (mean ± 95% CI) Ridge, independently optimized α (mean ± 95% CI)

Figure 4. Accuracy of surrogate models. The green curve shows the surrogate that optimizes the mS,z from Eq. (G44), keeping cS,z (α) fixed by taking α = α∗ obtained from the Hamiltonian embedding model. The orange curve optimizes over both mS,z and α (which determines cS,z (α)). Also shown in orange dots for comparison are the best (over bond dimensions χ) accuracies for the Hamiltonian embedding model for each trial.

function of the input feature and of the trainable weights, d

θ g ( x, α) = ∑ α g,r xr + α g,0 ,

g = 1, . . . , 12n,

(H1)

r =1

so the trainable parameters are α ∈ R12n×(d+1) and they are a total of M = 12n(d + 1). Every trainable gate is generated by a Pauli operator with P2 = I, so the cost admits the parameter-shift rule of Section II B. The two-qubit rotations are nearest-neighbor, so the circuit remains of constant depth with a light cone of constant size. Lastly, the label observable is set to A = n1 ∑in=1 Zi throughout. We consider a constant-depth encoding unitary U (α, x ) such that the observable A can be evolved under U (α, x ) classically. In particular, using backpropagation, we store a list of which Pauli terms appear in the embedding Hamiltonian together with the functional dependence of their coefficients from the features and trainable weights.

Classical iterative solver

We describe the training procedure when using the classical iterative solver. At the first training step, we start from the reference state |ref⟩ = |0n ⟩ and construct truncated approximations to the six vectors { HTk (α0 ) |ref⟩}5k=0 obtained by repeatedly applying the data Hamiltonian evaluated at the initial parameters α0 . After each application, we keep up to 8,000 computational-basis bitstrings with the largest amplitudes. We then diagonalize HT (α0 ) in the subspace spanned by the collected computational-basis states, obtaining an approximate ground state. Let us refer to this procedure of applying the Hamiltonian, collecting bitstrings, and diagonalizing as one “round”. We perform three rounds in total at fixed α0 , where the ground state from each round is used as the reference state for the next, and the diagonalization in each round is in the subspace spanned by all collected bitstrings from previous rounds. This gives us the approximate ground state for HT (α0 ). At subsequent training steps, we use an analogous procedure, taking the approximate ground state obtained at step t − 1 as the reference state. Here, each step involves a single round with only one application of HT (αt ), yielding truncated approximations to the two Krylov states { HTk (αt ) |ref⟩}1k=0 . Furthermore, for optimization steps t > 0, the truncation limit is reduced from 8,000 to 3, s000 computational-basis components per state, retaining those with the largest amplitudes.

53 3.

Optimization and model selection

The trainable parameters are initialized by drawing from independent Gaussian distributions α g,r ∼ N (0, 0.32 ), using a fixed seed for each run to ensure reproducibility. To limit overfitting, we add a squared-ℓ2 penalty to the training objective, breg (ψ, α) = C bT (ψ, α) + λ ∥α∥2 , C 2 T 2

(H2)

where α denotes the full parameter vector and λ = 0.005. At iteration t, we hold the state |ψt ⟩ used for gradient evaluation fixed and compute gt = ⟨ψt | ∇α HT (α) |ψt ⟩|α=αt + λαt .

(H3)

For the parameter update step, we use Adam with learning rate 0.05 and constants β 1 = 0.9, β 2 = 0.999, and ϵ = 10−8 . Training is performed for 80 iterations for the image-classification benchmark described in the next section. Model selection after training is performed by scoring the performance at each training epoch against a validation set that is split off the training data. In practice we hold out 20% of the latter, and the final reported model corresponds to the epoch with the highest validation accuracy. Classification uses the score f ( x ) = ⟨ψ| U † ( x, α) A U ( x, α) |ψ⟩ with label 0 assigned when f ( x ) ≥ 0 and label 1 otherwise. Appendix I: Image classification benchmark 1.

Protocol

The task considered is the binary classification of the Fashion-MNIST dataset [XRV17], restricted to the classes dress (label 3) and shirt (label 6), which is the setting used in Appendix L of Ref. [HBM+ 21]. Following that protocol, the two classes are pooled (14,000 images of 784 pixels each), the pool is reduced by principal component analysis to n components, each component is standardized to zero mean and unit variance, and Npool = 800 points are subsampled without replacement using the seed of the run. A stratified 70/30 train–test split of the 800 points is then drawn from the same seed, and a further stratified 80/20 split of the training part provides the inner validation set used for model selection. For every pair (qubit count, seed) the subsample and the splits are drawn once and shared by every method, so all comparisons are paired. 2.

Experimental Baselines

The quantum baselines are those of Ref. [HBM+ 21], evaluated on the same splits and with the same hyperparameters. Three data embeddings x 7→ |ϕ( x )⟩ are used: a product of single-qubit rotations, one per feature (E1); an IQP-type circuit in which two layers of Hadamards sandwich diagonal phases with linear terms x j Zj and all-toall quadratic terms x j x j′ Zj Zj′ (E2); and a Trotterized one-dimensional Heisenberg evolution for time t = n/3, with couplings set by the features, acting on a fixed Haar-random product state (E3). For each embedding we compute the fidelity quantum kernel Kij = | ⟨ϕ( xi )| ϕ( x j )⟩|2 and the projected quantum kernel built from the single-qubit reduced density matrices, Kij = exp[−γ ∑q ∥ρq ( xi ) − ρq ( x j )∥2F ], together with its generalizations in which the reduced density matrices of all pairs (2-RDM) or all triples (3-RDM) of qubits enter the exponent. Each kernel is fed to a support vector machine with a precomputed Gram matrix. A penalty term C is scanned over an 18-point grid and the bandwidth over γ ∈ {0.25, 0.5, 1, 2, 3, 4, 5, 20}/(n Var( X )) (see Ref. [HBM+ 21] for more details). Finally to compare against the performance of hardware runs (in Fig. 1 (b)), the best kernel results are selected based on their validation accuracy. Simulating these baselines faithfully at 50 qubits requires embedding-specific methods, and not every embedding reaches every n. The E1 states are products of single-qubit states, so both kernels admit exact closed-form solutions. The E2 circuit generates volume-law entanglement through its all-to-all diagonal phases and admits no efficient tensor-network contraction, so it is simulated by exact statevector simulations, which bounds its reach by memory. The fidelity kernel is computed up to n = 16, and the projected kernels up to n = 24 (n = 20 at 2-RDM order). The E3 evolution also generates extensive entanglement, and hence its exact matrix-product-state bond dimension grows as 2⌈n/2⌉ . We simulate it by time-evolving-block-decimation with 20 Trotter steps for the total time t = n/3, capping

54 the bond dimension as χ = min(256, 2⌈n/2⌉ ). With this choice, the MPS can represent exactly any quantum state up to n ≤ 16, and the only approximation error comes from the finite Trotter step, which we measured at n = 16 as a projected-kernel feature error of 1.6 × 10−2 . We do not evaluate E3 beyond n = 16, where χ = 256 would cease to be exact. We note that the classical hardness of simulating the E2 and E3 embeddings at large n does not, by itself, confer a learning advantage, as the corresponding kernels are precisely the ones that concentrate and fail on this task. For the classical baseline, we consider support vector machines with radial basis function and linear kernels, random forests, gradient boosting, adaptive boosting, and a multilayer perceptron. We use the same n-dimensional features for all methods. The hyperparameters for these methods follow the same grids of [HBM+ 21]. For each model class, hyperparameters are tuned by three-fold cross-validation on the training split. The selected model is then evaluated once on the test split.

Appendix J: Hardware experiments 1.

The hardware solver

The hardware runs replace only the sampling stage of the solver of Appendix H 2 by a quantum device, whereas the construction of HT , the subspace union, the projection, the diagonalization, the gradient and the Adam update are performed classically on a laptop. At iteration t the current HT (αt ) is converted into a Pauli operator and the circuits k e |ref⟩ 7−→ e−i δt HT (αt ) |ref⟩ ,

k = 0, . . . , kmax ,

(J1)

are built by first-order Lie–Trotter product formulas, with the step δt = π/(2 ∑i |ci |) of Appendix H 2, and measured e T is a pruned version of HT keeping only its n P Pauli terms of largest coefficient in the computational basis. Here H e T is determined as follows. First, Pauli terms with magnitude (identity excluded). The pruned Hamiltonian H − 3 coefficients |ci | < 10 max j |c j | are dropped from the set of operators in the expansion. Next, we keep the n P = 400 Pauli terms of largest magnitude and assess whether the deepest transpiled circuit has at most 4,000 two-qubit gates. We iteratively reduce n P by factors of 0.6 until this is the case. The size of this set of Pauli terms n P is kept fixed for the remaining training epochs. The resulting value of n P is 240 for the 50-qubit hardware run. For comparison, the full HT contains over 20,000 terms. Note that the pruning applies only to the evolution operator, while the projection uses the full HT . At the first iteration the reference is |0n ⟩ and the thirteen circuits k ∈ {0, . . . , 12} are measured with 25,000 shots each. At later iterations the references are the four highest-weight bitstrings of the previous ground state, the circuits k ∈ {1, 2, 3} are measured for each reference (12 circuits at 8,192 shots), and the entire previous subspace is added to the union, so that the subspace can follow the training trajectory of the Hamiltonian. Circuits are mapped to a chain of n physical qubits chosen from the calibration data at run time, to minimize the summed two-qubit-gate and readout errors along the chain by a greedy search over the coupling graph. These circuits are then executed with dynamical decoupling and measurement twirling enabled; no other error mitigation is used.

Model training insights

For every hardware run we compare the quantumly trained HE model to the classical iterative solver applied to the same data cell (Fig. 5). In all five runs, the accuracy of the hardware runs is greater, by +0.004 to +0.050 (mean +0.023). Finally, to test whether the seed-to-seed variation reflects genuine differences in learned models or merely finitesample noise in the 240-example test set, we perform a cross-seed transfer analysis in simulation with the classical iterative solver. We evaluate twenty independently trained cells on all twenty held-out test splits (Fig. 6). This allows us to compute the seed-to-seed variation for the model initialization (indexed by rows) and the data splits (indexed by columns). The near-uniformity of the map shows that the accuracy barely changes when a model is applied to a different data sample. Both the native (diagonal) and transfer (off-diagonal) accuracies average ≈ 0.885. Thus, a model trained on one data sample performs just as well on an independent sample, confirming that HE learns transferable, sample-independent structure from the data.

55 0.95 classical iterative solver HE trained on QPU

0.94 0.93

+0.050

+0.013

+0.021

test accuracy

0.92 0.91

+0.004

+0.029

0.90 0.89 0.88 0.87 0.86

trial-1

trial-2

trial-3

trial-4

trial-5

1.0

1000

.875

.892

.846

.875

.887

.863

.900

.879

.867

.875

.892

.912

.875

.863

.846

.875

.904

.892

1001

.887

.858

.850

.887

.871

.875

.863

.867

.892

.858

.879

.854

.879

.858

.908

.875

.896

.863

.871

.892

1002

.875

.887

.863

.867

.887

.863

.887

.867

.912

.858

.883

.896

.863

.896

.883

.908

.917

.896

.900

1003

.921

.904

.892

.900

.908

.912

.900

.912

.892

.917

.929

.917

.879

.925

.896

.933

.925

.917

1004

.863

.854

.867

.842

.858

.867

.883

.838

.908

.875

.883

.850

.858

.867

.879

.858

.854

.871

1005

.900

.917

.900

.917

.904

.917

.887

.917

.896

.917

.904

.917

.879

.921

.904

.929

.921

.946

.925

1006

.875

.879

.867

.896

.867

.858

.887

.875

.879

.850

.896

.887

.875

.854

.863

.904

.883

.892

.887

1007

.908

.883

.863

.875

.892

.900

.879

.904

.867

.896

.892

.879

.850

.912

.887

.900

.912

.917

1008

.925

.904

.896

.921

.900

.908

.904

.921

.887

.912

.908

.912

.887

.925

.912

.929

.925

.938

.929

1009

.883

.879

.854

.883

.875

.846

.875

.854

.887

.858

.896

.875

.838

.883

.887

.879

.863

.900

1010

.896

.867

.883

.892

.875

.879

.883

.904

.863

.908

.883

.846

.867

.875

.904

.887

.896

1011

.912

.904

.896

.900

.912

.883

.921

.879

.917

.908

.875

.921

.912

.929

.925

.917

.912

1012

.883

.900

.846

.896

.887

.879

.904

.863

.912

.867

.896

.904

.892

.863

.900

.892

.908

.912

.892

.883

1013

.904

.896

.879

.904

.887

.875

.904

.887

.912

.883

.912

.879

.900

.875

.900

.896

.908

.912

.904

1014

.875

.842

.887

.908

.867

.887

.871

.892

.863

.879

.904

.896

.863

.858

.892

.900

.892

.871

.925

1015

.875

.850

.842

.871

.863

.854

.875

.829

.875

.863

.867

.892

.896

.838

.887

.875

.867

.883

.887

.879

1016

.883

.875

.854

.867

.879

.863

.858

.850

.871

.858

.883

.900

.871

.842

.863

.871

.900

.875

.900

1017

.883

.863

.871

.887

.871

.896

.854

.896

.858

.904

.900

.883

.858

.900

.879

.892

.887

.892

1018

.875

.879

.854

.879

.850

.854

.879

.833

.879

.867

.854

.842

.875

.850

.883

.846

.871

1019

.879

.875

.858

.879

.900

.867

.871

.875

.896

.883

.900

.904

.892

.838

.887

.883

.917

.896

.912

.896

0.9

test accuracy

0.8

0.7

1019

1018

1017

1016

1015

1014

1013

1012

1011

1010

1009

1008

1007

1006

1005

1004

1003

1002

1001

0.6

1000

model trained on seed

Figure 5. A comparison of the test accuracy from the Hamiltonian embedding model, when used in conjunction with a classical iterative solver (blue circles) and the SKQD method on hardware (orange squares) for the dataset of Fig. 1(b).

0.5

evaluated on seed’s test split (data sample)

Figure 6. Cross-seed transfer at n = 50. Each cell is the test accuracy of a model trained on one seed (row) and evaluated on another seed’s held-out data split (column); the red boxes mark the diagonal, where a model is tested on its own split.

Appendix K: A numerical counterexample to Pauli classical surrogate models

In Appendix G, we presented a learning task for which the Hamiltonian embedding model provably requires fewer training examples than the corresponding classically regularized surrogate. Here, we present a different example as a proof of principle, numerically showing that retaining the quantum state during learning yields a substantial performance advantage over the corresponding Pauli-based classical surrogate. Specifically, we compare the performance of a particular Hamiltonian embedding model, defined by a fixed choice of encoding unitary U ( x, α) and observable O(y), with that of the corresponding Pauli classical surrogate. The latter takes the form e i ( α ), gα ( x ) = ∑ ϕi ( x, α) m i

(K1)

56 where the coefficients ϕi ( x, α) are defined through the Pauli decomposition K

U † ( x, α) AU ( x, α) = ∑ ϕi ( x, α) Pi ,

(K2)

i =1

I +(−1)y A

. Here, Pi are the K = poly(n) Pauli operators whose expectation values are measured by the where O(y) = 2 e i (α) ∈ R are unconstrained classical parameHamiltonian embedding model as explained in Section III B 2, while m ters replacing the physical quantum moments mi (α) = ⟨ψT ,α | Pi |ψT ,α ⟩ for the quantum model in Eq. (45). Thus, the two models use the same input-dependent parameterized Pauli feature map, but differ in that the coefficients of the Hamiltonian embedding model must arise as expectation values of a single physical quantum state, whereas those of the classical surrogate are optimized freely. Schematically, f HE ( x ) = ∑ ϕi ( x, α) mi (α),

mi (α) = ⟨ψT ,α | Pi |ψT ,α ⟩ ,

(K3)

gridge ( x ) = ∑ ϕi ( x, α) m̃i (α),

m̃i (α) ∈ R .

(K4)

For a fixed value of the trainable parameter vector α, the coefficients of the classical surrogate are obtained by solving the ridge regression problem " # 2 λ 1 L ⋆ 2 yℓ T m̃ (α) = arg min (K5) ∑ (−1) − ϕ(xℓ , α) m̃ + 2 ∥m̃∥2 . m̃ 2L ℓ= 1 Since the feature map ϕ( x, α) itself depends on α, the optimal vector m̃⋆ depends implicitly on the choice of the representation parameters. The corresponding loss is then minimized with respect to the trainable parameters α. 1.

Construction of the learning task

The task is constructed such that the requirement that the Pauli moments mi (α) be jointly compatible with a physical quantum state provides a useful inductive bias in the low-data regime. To emphasize this effect, we consider a teacher–student setting in which the inputs x are labeled according to y=

1 − s( x ) 2

h i s( x ) = sign ⟨ψ⋆ | U † ( x, α⋆ ) AU ( x, α⋆ ) |ψ⋆ ⟩ ,

(K6)

for a fixed teacher parameter α⋆ and teacher state |ψ⋆ ⟩. We then use the same parameterized unitary U ( x, α) in the dataset Hamiltonian, with α initialized randomly. The goal of the Hamiltonian embedding model is therefore to learn parameters α and the corresponding ground state |ψT ,α ⟩ such that the prediction rule in Eq. (5) reproduces the teacher labels. In particular, we choose the following setting. We consider a system of n = 2M qubits, grouped into M two-qubit blocks (2j − 1, 2j) for j ∈ {1, .., M }. The teacher state |ψ⋆ ⟩ is the product of Bell states on each block

|ψ⋆ ⟩ = Φ+

⊗M

Φ+ =

|00⟩ + |11⟩ √ . 2

(K7)

The inputs x are chosen to be a tuple of four different variables x = ( j, θ, β, q) with different roles: • j selects the active two-qubit block and it is sampled uniformly at random in j ∈ {1, ..., M }; • θ is sampled from a two-pole Gaussian mixture centered at 0 and π, θ∼

1 1 N (0, σθ2 ) + N (π, σθ2 ), 2 2

σθ = 0.12,

(K8)

where N (µ, σ2 ) denotes a normal distribution with mean µ and variance σ2 .

• β ∼ Unif[−π, π ];

• q = (q0 , q1 , q2 , q3 ) is a random unit vector such that q20 + q21 + q22 + q23 = 1. It parametrizes a random single-qubit unitary V (q) ∈ SU (2), which is applied to the active two-qubit block through V (q) ⊗ V (q)∗ .

57 The label-dependent observable is defined as O(y, j) =

I + (−1)y A j , 2

2 1 X ⊗ X − Y ⊗ Y. 3 3

(K9)

and A j denotes the operator A acting on the block (2j − 1, 2j) and identity elsewhere. The encoding unitary acts nontrivially only on the two-qubit block selected by the input index j. Specifically, U ( x = ( j, θ, β, q), α) = I ⊗2( j−1) ⊗ Uloc (θ, β, q; α j ) ⊗ I ⊗2( M− j) ,

(K10)

h i Uloc (θ, β, q; α j ) = C (θ ) Rn(θ ) β + α j h(θ ) ⊗ I S(q),

(K11)

C (θ ) = RY (θ/2) ⊗ RY (−θ/2), n(θ ) = (sin(θ/2), 0, cos(θ/2)) ,

(K12) (K13)

h(θ ) = sin2 (θ/2), S(q) = V (q) ⊗ V (q)∗ ,

(K14) (K15)

where

with

and V ( q ) = q0 I − i ( q1 X + q2 Y + q3 Z ) .

(K16)

We note that the unitary S(q) = V (q) ⊗ V (q)∗ has the Bell state |Φ+ ⟩ as eigenstate. In particular (V ⊗ V ∗ ) |Φ+ ⟩ = |Φ+ ⟩ for any single-qubit unitary V. Therefore q can rotate the Pauli representation seen by the learner while leaving the Bell teacher prediction unchanged. Finally, we also note that the local encoding Uloc (θ, β, q; α j ) can be implemented entirely using standard single-qubit rotation gates. In particular, Rn(θ ) ( φ) can be decomposed as RY (θ/2) R Z ( φ) RY (−θ/2), while V (q) is itself a single-qubit rotation parameterized by the quaternion q. We then label the data by the rule in Eq. (K6), with the choice of α⋆j = 0 for all j ∈ {1, ..., M }. The teacher corresponding expectation value is then Φ+ Uloc (θ, β, q; α⋆ )† A j Uloc (θ, β, q; α⋆ ) Φ+ =

1 + 2 cos θ cos β. 3

(K17)

We note that, although the observable A j explicitly depends on the input index j, the same construction can be written with a fixed observable by absorbing this dependence into the encoding unitary. Choosing A1 as the fixed e ( x, α) = observable and letting Pj swap the pair (2j − 1, 2j) with (1, 2), so that Pj† A1 Pj = A j , we can define U (1,2)

(1,2)†

(1,2)

e † ( x, α) A1 U e ( x, α) = P† U Uloc (θ, β, q; α j ) Pj , which gives U A1 Uloc Pj and therefore reproduces the original jj loc dependent construction. From the labeling rule in Eq. (K6) and Eq. (K17), the labels are then given by Eq. (K6) with s = sign[(1 + 2 cos θ ) cos β] .

(K18)

The nuisance variable q and the block index j do not appear in the labeling rule, but they make the learning problem for the classical surrogate more demanding. In particular, the index j determines which two-qubit block is probed by e j in the classical surrogate. On each example, with each block associated with its own vector of Pauli coefficients m the other hand, q changes the Pauli decomposition of the encoded observable. Indeed, conjugation by the random unitary V (q) rotates Pauli operators into different linear combinations of X, Y, and Z, so that different values of q make different examples probe different combinations of the Pauli coefficients. In the low-data regime, the training e which can fit the observed set may therefore be compatible with many choices of the free surrogate coefficients m, data while generalizing differently to unseen inputs. In contrast, the Hamiltonian embedding model restricts the coefficients m to be jointly compatible with the Pauli moments of a physical quantum state. This provides an additional inductive bias that favors coefficient vectors consistent with the teacher construction and, in our numerical experiment, leads to better generalization before the classical surrogate has received enough data to identify a wellgeneralizing solution.

700

Training samples L90

600 500 400 300 200 100 0

Qubits n Hamiltonian embedding

Classical surrogate (λ = 0.1)

Figure 7. Training examples needed to reach 90% test accuracy for Hamiltonian embedding and classical surrogate with λ = 0.1. The sweep uses 20 independent seeds. Central values are the first crossings of the mean test-accuracy curve across the 20 seeds, while the error bars are the standard deviations of the seed-wise first-crossing values. The bars therefore show seed-to-seed variability rather than a conventional confidence interval.

Table I. First tested training size L for which the mean test accuracy over all seeds reaches 90%. The λ sweep uses 20 independent seeds and 5000 test points per run. Model / λ HE 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 1

n=4 56 160 160 160 160 160 160 128 128 144

6 72 216 216 216 216 216 216 216 216 240

8 96 320 320 320 320 320 320 288 256 256

10 140 400 400 400 400 400 400 360 320 360

12 144 432 432 432 432 432 432 432 384 432

14 168 560 560 560 560 560 560 448 448 448

16 192 640 640 640 640 640 640 576 512 640

18 216 648 648 648 648 648 648 648 576 576

20 240 800 800 800 800 800 800 720 640 640

Numerical experiment

We compare the Hamiltonian embedding model and the Pauli classical surrogate for systems ranging from 4 to 20 qubits, with n ∈ {4, 6, 8, 10, 12, 14, 16, 18, 20}, using 20 independent random seeds and 5000 fresh test examples per run. For each system size, we vary the number of training examples supplied to both models. Specifically, for M = n/2 two-qubit pairs, the total training-set size is chosen as L = Ms, where s is varied over a fixed grid of samples per pair. For each (n, L, seed), we also vary the strength λ of the regularization term in the ridge regression of Eq. (K5). We then evaluate the two models on the same training/test data and with the same initial random α. As a common metric, for each model and system size we define L90 as the smallest training-set size L for which the mean test accuracy over the 20 random seeds reaches 90%. This quantity is plotted in Fig. 7. Because L is evaluated on a discrete grid, the reported values are first grid crossings. To quantify seed-to-seed variability, for each seed r we also define the first seed-wise crossing. The error bars used in Fig. 7 are the standard deviation of these 20 seed-wise first-crossing values. The plotted central point is the smallest training-set size L for which the mean accuracy over the random seeds reaches 90%. Therefore the bars should be interpreted as a diagnostic of run-to-run variability, rather than as a confidence interval around the central point. The complete results are shown in Table I. The first row is the Hamiltonian embedding model. The remaining

59 Table II. Hamiltonian embedding compared with the best ridge threshold observed over the tested λ grid. The final column gives the ridge-to-HE ratio averaged over system sizes. LHE 90 ridge,best L90

4 56 128

6 72 216

8 96 256

10 140 320

12 144 384

14 168 448

16 192 512

18 216 576

20 240 640

mean ratio – 2.62×

rows give the classical surrogate L90 for each ridge regularization value. The weakly regularized regime, from 10−8 through 10−3 , gives the same first-crossing thresholds on the tested grid. Moderate regularization improves ridge regression. The strongest improvement is typically around λ = 10−1 , although the optimum is not identical at every system size. To make the comparison maximally favorable to the classical surrogate, we also report the comparison between the L90 value for the Hamiltonian embedding and the classical surrogate at the best regularization strength λ, for ridge,best ridge each system size, i.e. L90 (n) = minλ∈Λ L90 (n, λ), where Λ is the tested grid. As the best performance for the classical surrogate is found at λ = 0.1 for almost all system sizes, we plot the comparison with the Hamiltonian embedding model in Fig. 7. Finally, we emphasize that the labeling rule in Eq. (K18) is deliberately simple, and a classical learner provided with additional structural information could easily solve the learning task. Nevertheless, the example serves to demonstrate that there exist settings in which the physical bias imposed by requiring the Pauli moments mi to arise from a valid quantum state can improve the performance of the Hamiltonian embedding model relative to the corresponding classical surrogate.

Supervised Learning with Quantum Hamiltonian Embeddings

Analyze & explore related research in Conceptio

Related documents