Predicting Ground State Properties: Constant Sample Complexity and Deep Learning Algorithms

Authors: marc  wanner
laura  lewis
chiranjib  bhattacharyya
devdatt  dubhashi
alexandru  gheorghiu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 14:18:44

Abstract: A fundamental problem in quantum many-body physics is that of finding ground states of local
Hamiltonians. A number of recent works gave provably efficient machine learning (ML) algorithms
for learning ground states. Specifically, [Huang et al. Science 2022], introduced an approach for learning
properties of the ground state of an $n$-qubit gapped local Hamiltonian $H$ from only $n^{\mathcal{O}(1)}$ data
points sampled from Hamiltonians in the same phase of matter. This was subsequently improved
by [Lewis et al. Nature Communications 2024], to $\mathcal{O}(\log 𝑛)$ samples when the geometry of the $n$-qubit system is known.
In this work, we introduce two approaches that achieve a constant sample complexity, independent
of system size $n$, for learning ground state properties. Our first algorithm consists of a simple
modification of the ML model used by Lewis et al. and applies to a property of interest known beforehand. Our second algorithm, which applies even if a description of
the property is not known, is a deep neural network model. While empirical results showing the
performance of neural networks have been demonstrated, to our knowledge, this is the first rigorous
sample complexity bound on a neural network model for predicting ground state properties. We also perform numerical experiments that confirm the improved scaling of our approach compared to earlier results.

Main paper: 
 introduction: 
 One of the most important problems in quantum many-body physics is that of finding ground states of quantum systems. This is due to the fact that the ground state describes the behavior of electronic systems (e.g., metals, magnets, etc.) at room temperature well. Thus, understanding This is especially relevant for systems approaching the thermodynamic limit, where the system size can be arbitrarily large. Needing fewer samples also means less work for experimentally preparing ground states of the system. The second question, stated as an open question in [1, 2] is: Question 2. Can we obtain rigorous sample complexity guarantees for neural-network-based ML algorithms for predicting ground state properties?We give positive answers to both questions. We consider the same assumptions as [2] with minimal additional ones that we mention here. First, we show that a simple modification to the approach in [2] allows us to achieve a sample complexity that is independent of the system size. This does, however require knowledge of the property we wish to predict in advance, whereas this is not a requirement in [2] . We view this as a reasonable assumption, since in practice we can imagine preparing ground states of some system in order to measure a specific property of interest. We show the following theorem, stated informally here. The formal statement, including all the assumptions required for proving the result, can be found in Appendix B.Theorem 1 (Informal). Let H(x) be an n-qubit gapped, geometrically local Hamiltonian with ground state ρ(x). Given an observable O, with a known decomposition as a sum of local Pauli operators and given training data {(x ℓ , y ℓ )} N ℓ=1 sampled from an arbitrary distribution, with y ℓ ≈ tr(Oρ(x ℓ )), Figure 1 : A deep network model for predicting ground state properties. Given a vector x ∈ [-1, 1] m that parameterizes a quantum many-body Hamiltonian H(x), the algorithm uses geometric structure to create "local" neural network models fθ P iPi . The ML algorithm then combines the outputs of these local models to predict a property tr(Oρ(x)), where ρ(x) is the ground state of H(x).Here, we decompose O = M i=1 α Pi P i for Pauli operators P i , where the final layer takes a linear combination of the outputs of the local models weighted by some trainable parameters w Pi that intuitively should approximate the Pauli coefficients α Pi .there is an ML algorithm for predicting ground state properties tr(Oρ(x)) to within precision ϵ > 0 using N = O 2 polylog(1/ϵ) training samples.Note that the number of samples N depends only on the desired prediction error ϵ and is independent of the system size. In particular, this means that for a fixed prediction error our algorithm requires only a constant amount of training data. Moreover, the computational complexity of our algorithm improves upon [2] , having O(n) runtime, compared to the previous O(n log n). While removing the log n factor may seem like a small improvement, in practice this can make a significant difference. For instance, for a system of n ∼ 1000 qubits, removing the log n factor would result in a ten-fold reduction in training data and time.Much like in [2] , this result also extends to learning classical representations of ρ(x). In other words, if the algorithm is instead given classical shadows [52] of the ground state as training data, it can then predict a classical representation of ρ(x) for new parameters x. This can mitigate the requirement that the observable is known in Theorem 1, as predicting properties from a classical representation clearly requires knowledge of the observable.Our second result shows the same sample complexity guarantee for a neural network ML algorithm (Figure 1 ) [53] , in which one does not need to know the observable being measured in advance. An additional constraint that we require in this case is that the training data is not sampled according to an arbitrary distribution, but a distribution satisfying some technical assumptions. We note that these assumptions are satisfied for common distributions such as uniform and Gaussian.Theorem 2 (Informal). Let H(x) be an n-qubit gapped, geometrically local Hamiltonian with ground state ρ(x). For any observable O, expressible as a sum of local Pauli operators and given training data {(x ℓ , y ℓ )} N ℓ=1 , sampled from a distribution satisfying certain assumptions with y ℓ ≈ tr(Oρ(x ℓ )), there is a neural network ML algorithm for predicting ground state properties tr(Oρ(x)), for uniform x, to within precision ϵ > 0 using N = O 2 polylog(1/ϵ) training samples under mild assumptions on training.We prove this result by making use of the Koksma-Hlawka inequality from quasi-Monte Carlo theory [54, 55, 56, 57, 58] and combining it with the spectral flow formalism [59, 60, 61] .Similar to Theorem 1, we can also extend this result to learning classical representations of ρ(x) when given classical shadow training data. The formal statement and its proof can be found in Appendix C. We also remark that, much like the setting in [51] , a more favorable scaling with respect to ϵ can be achieved if the number of parameters that the Hamiltonian depends on is constant. In particular, it was shown in [51] that if the number of parameters is constant, the sample complexity scales as N = poly(1/ϵ, log(n)). For our results, this similarly yields N = poly(1/ϵ), preserving the independence on the system size while also achieving a polynomial scaling in 1/ϵ. Furthermore, we perform numerical experiments on system sizes of up to 45 qubits, which support our theoretical findings, and show that, in practice, our deep learning algorithm outperforms previous methods [2] . We describe them in detail in Appendix D, and they are illustrated in Figure 2 . 
problem statement: 
 First, we formally describe the problem setting, which is the same as [2] . We consider a family of n-qubit Hamiltonians H(x) smoothly parameterized by an m-dimensional vector x ∈ [-1, 1] m . We assume that these Hamiltonians are gapped for all choices of parameters x ∈ [-1, 1] m and geometrically local such that they can be written as a sum of local termsEQUATIONwhere the parameter vector x is a concatenation of the constant-dimensional vectors ⃗ x 1 , . . . , ⃗ x L . Each of these constant-dimensional vectors ⃗ x j parameterizes the local interaction term h j (⃗ x j ). Crucially, we assume that each local term h j only depends on a constant number of parameters rather than the entire parameter vector x. We also assume that the geometry of the n-qubit system is known.Throughout this work, we use ρ(x) to denote the ground state of the Hamiltonian H(x) and O to denote an observable that can be written as a sum of geometrically local observables with bounded spectral norm ∥O∥ ∞ ≤ 1. Here, the ground states ρ(x) form a gapped quantum phase of matter. Given samples of quantum states drawn from this phase, we wish to predict expectation values of observables O with respect to other states in the same phase. In other words, we are given training data {(x ℓ , y ℓ )} N ℓ=1 , where y ℓ ≈ tr(Oρ(x ℓ )) approximates the ground state property for a parameter choice x ℓ ∈ [-1, 1] m sampled from some distribution D over the parameter space. We aim to learn a function h * (x) that approximates the ground state property tr(Oρ(x)) for some unseen parameter x while minimizing the amount of training data, or sample complexity, N . How well we learn the ground state property is quantified by the average prediction error (2.2)We describe the precise conditions under which Theorems 1 and 2 hold in the following sections.In this section, we review the previous algorithm from [2] , as our proofs rely on similar ideas. For full details, we refer the reader to [2] and our more detailed presentation in Appendix A.1. The ML algorithm proposed in [2] requires some geometric definitions. Fix a geometrically local Pauli observable P ∈ {I, X, Y, Z} ⊗n .Let δ 1 , B > 0 be efficiently-computable hyperparameters that we define in Appendix A.1. Define the set I P of coordinates c such that x c parameterizes some local term h j(c) that is close to the Pauli P .Here, the distance between two observables d obs is defined as the minimum distance between the qubits that the observables act on, where the distance between qubits is given by the geometry of the system, which we assume to be known. Formally, we define this set of local coordinates asEQUATIONwhere h j(c) is the local term in the Hamiltonian H(x) whose parameters ⃗ x j(c) include the variable x c . The intuition behind this set of coordinates is that it indexes the parameters x c that influence the ground state property tr(P ρ(x)) corresponding to the Pauli P . The algorithm consists of two steps. First, it maps the parameter space [-1, 1] m to a high dimensional space via a nonlinear feature map ϕ. Second, it runs ℓ 1 -regularized linear regression (LASSO) [62, 63, 64] over the feature space.This first step encodes the geometry of the problem. The feature map intuitively projects a given parameter vector x onto the local parameter space {x c : c ∈ I P }. We define this precisely in Appendix A.1. Following the feature mapping, the ML algorithm uses LASSO [62, 63, 64] to learn functions of the form {h(x) = w • ϕ(x) : ∥w∥ 1 ≤ B} for a chosen hyperparameter B > 0. We denote the learned function by h * (x) = w * • ϕ(x). For our purposes, we set B = 2 O(polylog(1/ϵ1)) . This algorithm obtains the following rigorous guarantee.Theorem 3 (Theorem 1 in [2] ). Given n, δ > 0, 1 e > ϵ > 0 and a training data set {(x ℓ , y ℓ )} N ℓ=1 of size N = log(n/δ)2 polylog(1/ϵ) , (2.4) where x ℓ is sampled from an unknown distribution D and |y ℓ -tr(Oρ(x ℓ ))| ≤ ϵ. With a proper choice of the efficiently computable hyperparameters δ 1 , δ 2 , and B, the learned function h A crucial step in the proof is that ground state properties can indeed be approximated by linear functions over the feature space. Along the way, [2] proves that the ground state property can be approximated by a linear combination of "local functions," which are local in that they only depend on parameters with coordinates in the set I P . We relegate further details to Appendix A.1 and [2] .A, B 1 , . . . , B m satisfy A ∪ m j=1 B j = [k] and A ∩ B j = ∅, B i ∩ B j = ∅.3. For each (B l , l) pair, it holds that i > l for all i ∈ B l .This notation gives a compact way of expressing the terms of the mixed derivative and allows us to address each mixed derivative of the terms D êj individually. Each term ‚ l∈A ∂ B l l contains the product of m + 1 matrix-valued functions (including ρ(x), which depend on x. The set A denotes the partial derivatives on the factor ρ(x), which have been differentiated using Equation (C.174) when applying the product rule. We will address the partial derivatives on D êj later in this section, when we derive an upper bound for ‚ l∈A ∂ B l l.The first condition underlines that the first partial derivative on ρ(x) is necessarily computed via Equation (C.174) and thus contained in each term. The second condition reflects that each partial derivative operates on exactly one factor in each summand when applying the product rule. The third condition arises from the order, by which the partial derivatives are computed. For example, when we apply ∂ ∂x ′ j after ∂ ∂xj , the ∂ ∂xj D ê′ j can not occur in any term, since no term contained D ê′ j when the partial derivatives ∂ ∂xj were computed.We can show that the mixed partial derivatives of αP tr(P ρ(x)) can be written in terms of ‚ l∈A ∂ B l l.Lemma 18 (Mixed derivative). Let A k = {A ⊆ [k] : 1 ∈ A} and B A be as in Definition 8. The mixed derivative of the ground state property tr(P ρ(x)) is given by∂ k ∂x 1 . . . ∂x k α P tr(P ρ(x)) = α P A∈A k (B1,...,B |A| )∈B A ‚ l∈A ∂ B l l. (C.183)Proof. We proceed via induction. First, we verify that∂ ∂x k ‚ l∈A ∂ B l l with A ∈ A k-1 and A ∪ m j=1 B j = [k -1]yield summands which fulfill the criteria for summands of the k-th partial derivative stated in Definition 8. Then, we show that each summand of the k-th derivative stems from a unique summand from the (k -1)-th partial derivative.For the first part, let |A| = m. Furthermore, letI s (j) = {j} if s = k ∅ else . (C.184) Then, ∂ ∂x k ‚ l∈A ∂ B l l = m j=1 ‚ lj ∈A ∂ B l j ∪Ij (l) l + ‚ l∈A∪{k} ∂ B l l, (C.185)where each summand fulfills the properties, since k > l for all l ∈ A.Next, let A ′ ∈ A k and let B A ′ be the corresponding collection. Then, it is easy to see that, if k ∈ B lj , it stems from the summand‚ l∈A ′ ∂ B l l with B l1 , . . . , B lj \ {k}, . . . , B lm . If k ∈ A ′ , it stems from ‚ l∈A\{k} ∂ B l l.With this expression in hand, we can move forward to bounding the mixed derivative, which we need in order to bound the Hardy-Krause variation. First, we will upper bound the number of terms in ‚ l∈A ∂ B l l. Then, we derive an upper bound on the individual terms. For the first step, we exploit the conditions on the sets defining the mixed derivative. In the next step, we aim to bound each individual term of the mixed derivative ‚ l∈A ∂ B l l. Therefore, a crucial step is to bound the spectral norm of each factor. We first derive a preliminary result on the mixed derivatives of the factors in D êj (x), which depend on x. This can be done using Duhamel's Formula for the derivative of the exponential map on e H(x) , where we exploit that we only compute the derivative with respect to one parameter at a time, such that we can treat H(x) as a function, which only depends on one parameter. Theorem 18 (Derivative of the exponential map; Theorem 3a in [95] )). Let A(t) : R → C n×n . Then,d dt e A(t) = 1 0 e (1-s)A(t) dA(t) dt e sA(t) ds. (C.187) Lemma 20. Let k ∈ [n], B ⊆ [n] \ {k}, such that ∂ |C| hj ∂x C ∞ ≤ 1 ∀C ⊆ B ∪ {k}. Then ∂ |B| ∂x B e itH(x) ∂h j ∂x k e -itH(x) ∞ ≤ 2 |B|+1 (|B| + 1) |B|+1 (C.188)Proof. By Theorem 18, the mixed derivative equals the sum of terms of the formT = 1 0 • • • 1 0 l f l (s l )ds l . . . ds 1 , (C.189)where f l (s l ) can be any of e (1-s l )iH(x) , e s l iH(x) , ∂ l hj ∂x B l or 1. By our assumption and the Cauchy-Schwartz inequality, each term T satisfies ∥T ∥ ∞ ≤ 1. Furthermore, by the product rule, the number of terms is smaller than |B| j=1 (2j + 1). Since each term of the lth partial derivative (including k) is the product of at most 2l + 1 factors depending on x, such that the (l + 1)th derivative contains at most 2l + 1-times as many factors. Thus, the number of terms is bounded above by|B| j=1 (2j + 1) ≤ |B|+1 j=1 (2j) = 2 n n! ≤ 2 |B|+1 (|B| + 1) |B|+1 , (C.190)as required.Now we can bound the terms ‚ l∈A ∂ B l l.Lemma 21 (Bound components of the derivative). Let ‚ l∈A ∂ B l l be as in Definition 8. Then‚ l∈A ∂ B l l ≤ (2C γ ) |A| |A| s=1 2 |B ls |+1 (|B ls | + 1) |B ls |+1 . (C.191) Proof. Recall that D û(x) = +∞ -∞ W γ (t)e itH(x) ∂H ∂ û (x)e -itH(x) dt, (C.192)where W γ (t), such thatEQUATIONγ|t| log 2 (γ|t|) γ|t| > θ, (C.193)where θ is chosen to be the largest real solution of 35e 2 (γ|t|) 4 exp -2 7 γ|t| log 2 (γ|t|) = 1/2. It is also useful to note that sup t |W γ (t)| = 1/2.By definition of the terms ‚ l∈A ∂ B l l, the Cauchy-Schwartz inequality, and ∥ [A, B]∥ ∞ ≤ 2∥A∥ ∞ ∥B∥ ∞ , we obtain ‚ l∈A ∂ B l l ≤ +∞ -∞ |W γ (t)(2C γ ) |A| |A| s=1 2 |B ls |+1 (|B ls | + 1) |B ls |+1 ≤ (2C γ ) |A| 2 k |A| s=1 k |B ls |+1 ≤ (C 1 k) k , (C.|(C 1 k) k ≤ |α P |C k 1 k 2k = |α P |2 O(k log(k)) . (C.199)Note that when deriving this result, we do not require that the parameters for the mixed derivatives are distinct. Assuming that ∥H(x)∥ W k,∞ ([-1,1] m ) ≤ 1, we can induce an order recover the above bound for any mixed derivative of order k.Corollary 8. If ∥H(x)∥ W k,∞ ([-1,1] m ) ≤ 1, then ∥α P tr(P ρ(x))∥ W k,∞ ≤ |α P |2 O(k log(k)) .Proof. Note that the bound from Lemma 22 is agnostic to the explicit directions êj of the derivatives. Thus, we can choose any mixed derivative λ ∈ N k 0 such that k j=1 λ j = k and fix an order o : dom(λ) → [k]. Then, we can bound the mixed derivative on o(λ) using the same approach as in Lemma 22 to obtain the bound.In this section, we discuss the numerical experiments in detail.As in [2] , we consider the two-dimensional antiferromagnetic Heisenberg model with spin-1/2 particles placed on sites in a two-dimensional lattice. The corresponding Hamiltonian isH = ⟨ij⟩ J ij (X i X j + Y i Y j + Z i Z j ), (D.1)where ⟨ij⟩ denotes all pairs of neighboring sites on the lattice. The coupling terms J ij correspond to the parameters x of the Hamiltonian and are sampled uniformly from [0, 2] (and then mapped to lie in [-1, 1] for our ML algorithm). The goal of the numerical experiment is to predict the two-body correlation functions, i.e., the expectation value ofC ij = 1 3 (X i X j + Y i Y j + Z i Z j ) (D.2)for all neighboring sites ⟨ij⟩.To this end, we generate data similarly to [1, 2] , approximating the ground state and corresponding correlation functions for the Hamiltonian Equation (D.1) of different lattice sizes and choices of coupling parameters J ij . We consider lattice sizes of 4 × 5 = 20 up to 9 × 5 = 45. For each lattice size, we generate two datasets of size 4096, one with uniformly randomly distributed J ij and one where the coupling parameters are distributed as a Sobol sequence. We obtained the data by approximating the ground state using the density-matrix renormalization group (DMRG) [96] based on matrix-product-states (MPS) [97] , as has been done in [1, 2] . The simulations were performed on Nvidia T4 and A40 graphical processing units (GPUs). The former were used for lattice sizes from 4 × 5 up to 7 × 5 while the latter were used for lattice sizes 8 × 5 and 9 × 5. Depending on system size, we required between ≈ 50 and 200 hours on the respective hardware component to simulate one dataset of size 4096.Our deep learning model was also trained on Nvidia T4 and A40 GPUs. We trained the models for all respective correlation terms in parallel, by training a full model f Θ,w ij (we omit the indices for the model's parameters) for each term and minimizing the combined loss function⟨ij⟩ N ℓ=1 |f Θ,w ij (x ℓ ) -(C ij ) ℓ | 2 (D.3)for the sake of time efficiency. For each data point, we trained a combined model for 500 epochs.For the terms of the local models f θ P P , as defined in Definition 6, we used fully connected deep neural networks with five hidden layers of width 200. For training, we used the AdamW optimization algorithm [83] . Depending on the system size and the amount of training data, this took between 0.5 and 20 hours. As a baseline, we compared against the best model from [2] . The code can be found at https://github.com/marcwannerchalmers/learning_ground_states.git.In this section, we discuss the results of the numerical experiments and additional experiments performed that are not mentioned in the main text.First, we perform additional experiments that analyze the scaling of the training/prediction error with respect to various parameters such as system size, local neighborhood size, and training set size (Figures 4 to 6 ). Importantly, in each of these, we see that the training error is small, as required by Theorem 5. Thus, as discussed in the main text, this assumption is satisfied in practice.Moreover, as shown in Figure 2 (Left), the empirical prediction accuracy (RMSE) of the deep learning model is approximately constant with respect to the size of the lattice. Figure 4 (Right) further underlines this statement. The slight increase in prediction error for δ 1 > 0 (size of the local neighborhood in Equation (A.2)) present in Figure 4 (Right) when increasing the system size from 4 × 5 to 5 × 5 may occur due to numerical errors in the data. From system size 5 × 5 onwards, we rather witness random fluctuations in test errors than a systematic increase. Furthermore, we observe that the deep learning model significantly outperforms the regression model with random Fourier features from [2] . On the one hand, we notice that the performance of the latter could be improved, since the hyperparameters considered for hyperparameter tuning were selected for a substantially smaller dataset. This is underlined by the drop in RMSE for the regression model on Figure 5 for δ 1 = 1, whereas a smaller RMSE is possible when choosing δ 1 = 0. On the other hand, we think that the vast body of deep learning research also offers room for practical improvement of our deep learning model.For δ 1 = 0, we believe that our model achieves the best possible prediction error. For training set size larger than 2048, there is little improvement on the prediction error, as opposed to all experiments with δ 1 > 0 (see Figure 2 (Center)). Furthermore, the training error remains relatively large compared to other choices of δ 1 (see Figure 4 ). Hence, we conclude that the error arising from approximating the ground state property via local functions dominates the prediction error.When increasing δ 1 , we witness an increase in prediction error, especially for small training sets. This is consistent with Lemma 10, which states that the bound on the prediction error is a combination of the training error and a term proportional to the star-discrepancy (and thus increases with the dimension of the domain of the local models). Our experimental results underline the balance which must be achieved between the two in order to obtain a small prediction error. This can clearly be observed in Figure 6 . The training error decreases when increasing δ 1 and increases with the size of the training set. Meanwhile, the test error increases when increasing δ 1 and decreases with the size of the training set.Another interesting observation is that the ML algorithm's performance on LDS seem to be almost the same as that of uniformly random points. We believe this is due to the dominance of the local approximation error for small δ 1 and the drastic increase in dimensionality of the local models with increasing δ 1 outweighing the benefit of using LDS in practice. The dominance of approximation error is also a possible explanation for the slight decrease in prediction error with respect to the system size in Figure 2 (Left) and Figure 5 . For our concrete choice of lattice shape and ground state properties, the local approximation error may be decreasing with respect to system size. However, we do not expect this to be the case in general. of terms increases quadratically in system size, the norm of the weights in the final layer can not be bounded by a constant anymore. Furthermore, the properties of the local approximation do not hold true anymore. Hence, the predictive capabilities of a model with δ 1 = 0 may be more limited here than in the geometrically local case. However, the prediction error may also be impacted by possible numerical errors in the training data, as well as the architecture of the local deep neural networks. Overall, these experiments illustrate the necessity of the geometric locality assumption in our theoretical results.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: This research conforms with the code of ethics. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: There are no societal impacts of the work performed. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
main results: 
 In this section, we discuss our rigorous guarantees for predicting ground state properties with constant sample complexity and with neural-network-based ML algorithms.In this section, we show that a simple modification of the algorithm from [2] can achieve a sample complexity that is independent of the system size n, under the additional assumption that the observable O is known. In practice, a scientist often has a specific ground state property in mind that they wish to study, so we view this as a natural assumption. Moreover, this is still an interesting learning problem, as when obtaining the training data via quantum experiments, preparing the ground state ρ(x) in the laboratory for a new choice of parameters x may be difficult experimentally. This in turn means that accurately predicting some property tr(Oρ(x)) for a new choice of x may be challenging, even if the property of interest, O, is known. ML algorithms can allow us to circumvent this issue and generalize from the results of few training data points without needing to prepare the ground state directly.Let O = P ∈{I,X,Y,Z} ⊗n α P P be an observable that can be written as a sum of geometrically local observables. Because O is assumed to be known, we can find this decomposition of O in terms of the Pauli observables P . The overall structure of the algorithm remains the same: perform a nonlinear feature mapping followed by linear regression. However, there are two key differences from the previous algorithm [2] . First, we change the feature mapping of [2] to incorporate the Pauli coefficients α P . We define this new feature mapping φ in Appendix B. The second difference from [2] is that we use ridge regression [65, 66] instead of LASSO [62, 63, 66] . Recall that LASSO learns hypothesis functions of the form {h(x) = w • φ(x) : ∥w∥ 1 ≤ B} for some hyperparameter B > 0. In contrast, ridge regression replaces the ℓ 1 -norm constraint ∥w∥ 1 ≤ B with an ℓ 2 -norm constraint: ∥w∥ 2 ≤ Λ, for some hyperparameter Λ > 0. Namely, for a chosen efficiently-computable hyperparameter Λ > 0, ridge regression finds a vector w * that minimizes the training error subject to the constraint that ∥w∥ 2 ≤ Λ, i.e., minw∈R m ϕ ∥w∥ 2 ≤Λ 1 N N ℓ=1 |w • φ(x ℓ ) -y ℓ | 2 . (3.1)Standard results in machine learning theory give sample complexity upper bounds for ridge regression in terms of Λ and the ℓ 2 -norm of the feature vector φ(x) [65, 66] . The key idea is that with our new feature mapping, we can still approximate the ground state property by a linear function over the feature space, as in [2] , to obtain a low training error. Meanwhile, by incorporating the Pauli coefficients α P into the feature map, we can bound the ℓ 2 -norm of φ(x) by a quantity independent of system size, leveraging bounds on the ℓ 1 -norm of the Pauli coefficients [2, 67] . We note that naively applying ridge regression with the feature map from [2] does not achieve the same guarantees and in fact gives worse scaling than [2] . Similarly, we can also choose a suitable Λ > 0 independent of system size. Thus, we obtain the following guarantee.Theorem 4 (Constant sample complexity). Given n, δ > 0, 1/e > ϵ > 0 and a training data set We compare this result to Theorem 3. For a constant prediction error ϵ = O(1), our proposed algorithm achieves a constant sample complexity N = O(1), compared to the logarithmic sample complexity N = O(log n) of [2] . Moreover, we also improve the computational complexity, achieving a linear-in-n runtime, compared to the previous O(n log n). The scaling with respect to the prediction error ϵ is the same as the previous algorithm [2] . This means that regardless of how large our quantum system is, we need the same amount of samples to predict ground state properties well. This is especially important for settings in which obtaining training data for large systems is difficult.EQUATIONThus far, we have only considered the setting in which we learn a specific ground state property tr(Oρ(x)) for a fixed observable O. Because our training data is given in the form {(x ℓ , y ℓ )} N ℓ=1 , where y ℓ approximates tr(Oρ(x)) for this fixed observable O, if we want to predict a new property for the same ground state ρ(x), we would need to generate new training data. Thus, it may be more useful to learn a ground state representation, from which we could predict tr(Oρ(x)) for many different choices of observables O without requiring new training data. In this case, suppose we are instead given training data {x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 , where σ T (ρ(x ℓ )) is a classical shadow representation [52, 68, 69, 70, 71] of the ground state ρ(x ℓ ). An immediate corollary of Theorem 4 is that we can predict ground state representations with the same sample complexity. This follows from the same proof as Corollary 5 in [2] .Corollary 1 (Learning representations of ground states). Let n, δ > 0, 1/e > ϵ > 0 and δ > 0.Given training data {(x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 of sizeEQUATIONwhere x ℓ is sampled from D and σ T (ρ(x ℓ ) is the classical shadow representation of the ground state ρ(x ℓ ) using T randomized Pauli measurements. For T = Õ(log(n/δ)/ϵ 2 ), with probability at least 1δ, the ML algorithm will produce a ground state representation ρN,T (x) that achievesE x∼D | tr(O ρN,T (x)) -tr(Oρ(x))| 2 ≤ ϵ (3.5)for any observable with ∥O∥ ∞ ≤ 1 that can be written as a sum of geometrically local observables.In this section, we prove the existence of a deep neural network model that can predict ground state properties using a constant number of training samples. In particular, we prove that after training on a constant number of samples from a distribution D on [-1, 1] m satisfying certain technical assumptions, our model can achieve a low prediction error under mild assumptions on training. In this case, for predicting properties tr(Oρ(x)), the observable O need not be known in advance. However, we need to assume that all mixed first order derivatives of the Hamiltonian ∥∂ m H(x)/∂x 1 • • • ∂x m ∥ ∞ ≤ 1 exist and are bounded. This is not much stronger than [2] , which assumes that directional derivatives ∂h j /∂ û are bounded by one for any direction û. Moreover, we also need the training data to be sampled from a distribution D with probability density function g satisfying the following assumptions: g has full support and is continuously differentiable on [-1, 1] m . Also, g is of the form g(x) = L j=1 g j (⃗ x j ). This resembles our assumption on the form of the Hamiltonian H(x). Furthermore, the average prediction error is measured with respect to the same distribution D. We note that these assumptions are satisfied for common distributions such as uniform and Gaussian.As in the previous algorithm [2] , we leverage the geometry of the n-qubit system to approximate the ground state properties by a linear combination of smooth local functions, which only depend on parameters with coordinates in the local coordinate set I P defined in Equation (2.3) . Crucially, the size m ≜ |I P | of the domains of these local functions is independent of the system size.Instead of using a feature map and linear regression to learn the ground state properties, we utilize a deep neural network model defined as follows. Inspired by the local approximation of ground state properties, we define "local models" f θ P P : [-1, 1] m → R, which are neural networks consisting of three layers of affine transformations and applications of a nonlinear activation function. In particular, f θ P P has two hidden layers with the affine transformations given by the trainable weights and biases denoted by θ P . We take hyperbolic tangent, tanh, as the activation function. These local models are then combined into a model f Θ,w : [-1, 1] m → R given by f Θ,w (x) = P ∈S (geo)EQUATIONwhere w P ∈ R are the weights in the last layer and Θ = {θ P } P ∈S (geo) . This model is schematically illustrated in Figure 1 . We refer to Definition 6 in Appendix C for a full description of the model.Consider training data {(x ℓ , y ℓ )} N ℓ=1 , where x ℓ are sampled according to a distribution D satisfying the assumptions described above and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ. The ML algorithm first initializes the weights via standard deep learning initialization procedures, e.g., Xavier initialization [72] . Then, the algorithm performs quasi-Monte Carlo training given the training data, e.g., Adam [73] , to find weights Θ * , w * which minimize the training objective functionEQUATIONwhere λ is some regularization parameter that may depend on ϵ. For this algorithm, we prove the following theorem bounding the average prediction error of our deep neural network model.Theorem 5 (Neural network sample complexity guarantee). Let 1/e > ϵ > 0. Let D be a distribution with probability density function g satisfying the properties stated above. LetEQUATIONwhere the x ℓ 's are sampled from D and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ. Suppose that f Θ * ,w * achieves a value no larger than O(ϵ) on the training objective (Equation(3.7)) with λ(ϵ) = O(ϵ). Additionally, suppose that all parameters Θ * i of f Θ * ,w * satisfy |Θ * i | ≤ W max , for some W max > 0 that is independent of n. Then E x∼D |f Θ * ,w * (x) -tr(Oρ(x))| 2 ≤ ϵ. (3.9)Similar to Theorem 4, for a constant prediction error ϵ = O(1), the deep neural network algorithm achieves constant sample complexity N = O(1). In contrast to Theorem 4, we do not require knowledge about the observable, O. This is a direct consequence of the regularity of w P , which is achieved when the training objective is small. Theorem 6 guarantees, that a model with such regularity can yield a small prediction error.There are, however, some caveats compared to the previous result. First, the training data is restricted to being sampled from a distribution satisfying our technical assumptions stated previously, in contrast to Theorem 4 which holds for data sampled from any arbitrary unknown distribution. Second, in regards to the model, the weights must be bounded by a constant W max . Finally, we cannot guarantee a priori that the network will indeed achieve a low training error. This is due to the fact that our training objective is non-convex and thus, globally optimal weights cannot be found efficiently in general [74] . Even so, we are still able to prove the existence of suitable weights such that the resulting network approximates tr(Oρ(x)) for any x ∈ [-1, 1] m (see Theorem 6 in the next section).However, we view the assumptions made in Theorem 5 as being mild in practice. Small training objectives are commonly achieved in deep learning so we expect our training algorithm to produce a model which fulfills the assumptions of Theorem 5 after O(1) training steps and O(n) runtime. Moreover, it is known that gradient descent provably converges to the global optimum for overparametrized deep neural networks, while the weights remain small, when properly initialized [75] . We verify that these conditions are satisfied in practice through our numerical experiments in Figure 2 and ?? . To our knowledge, Theorem 5 is the first rigorous sample complexity bound on a neural network model for predicting ground state properties.We also note that if the training data is instead sampled according to a low-discrepancy sequence (LDS) [55, 56, 76, 77, 78, 57, 79, 58, 80, 81] , we can obtain better guarantees, but these improvements are hidden in the polylogarithmic factors in the exponential. We discuss learning given data from a LDS in Appendix C. Intuitively, a LDS is a collection of points in the parameter space that covers the space such that there are no large gaps, or discrepancies.Similar to Corollary 1, if we are instead given training data {x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 , where σ T (ρ(x ℓ )) is a classical shadow representation [52, 68, 69, 70, 71] of the ground state ρ(x ℓ ), then we obtain the following immediate corollary of Theorem 5. Corollary 2 (Learning representations of ground states with neural networks). Let n, δ > 0, 1/e > ϵ > 0 and δ > 0. Given training data {(x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 of sizeEQUATIONwhere x ℓ is sampled from a distribution D satisfying the same assumptions as Theorem 5 and σ T (ρ(x ℓ ) is the classical shadow representation of the ground state ρ(x ℓ ) using T randomized Pauli measurements. For T = Õ(log(n/δ)/ϵ 2 ), with probability at least 1δ, the ML algorithm will produce a ground state representation ρN,T (x) that achievesE x∼D | tr(O ρN,T (x)) -tr(Oρ(x))| 2 ≤ ϵ (3.11)for any observable with ∥O∥ ∞ ≤ 1 that can be written as a sum of geometrically local observables.To prove Theorem 5, we first show that our neural network model f Θ,w can approximate the ground state properties well. In particular, we show that there exist weights Θ ′ , w ′ such that f Θ ′ ,w ′ approximates the ground state properties and thus achieves small value for the training objective (Equation (3.7)). Then, we bound the prediction error using tools from deep learning and quasi-Monte Carlo theory [54, 55, 56, 57, 58] . We ensure the existence of f Θ ′ ,w ′ in the following theorem.Theorem 6. For any 1/e > ϵ > 0 and width W , there exist weights Θ ′ , w ′ such that the neural network modelf Θ ′ ,w ′ satisfies |f Θ ′ ,w ′ (x) -tr(Oρ(x))| 2 ≤ ϵ, ∀x ∈ [-1, 1] m . (3.12)Moreover, each parameter Θ i of the network has a magnitude of|Θ i | = 2 O(polylog(1/ϵ)) .This implies that for a suitable choice of regularization parameter λ = O(ϵ), the training objective from Equation (3.7) is also small. We prove this statement by combining results in deep learning regarding tanh neural networks approximating functions [53] with the geometric locality of the system and smoothness of the ground state properties. We note that the weights Θ ′ , w ′ in Theorem 6 are not necessarily the weights Θ * , w * found via the neural network training procedure in Theorem 5.Because the training objective is non-convex, we cannot guarantee convergence to these weights Θ ′ , w ′ . However, assuming that f Θ * ,w * does indeed achieve a low training error (which is often satisfied in practice), we are able to rigorously guarantee that the model will generalize well and achieve a low prediction error in Theorem 5.Notice that the guarantee of Theorem 6 holds for all x and, in particular, does not require our assumptions on the distribution D. The assumption that the network is trained on such data only becomes relevant when bounding the prediction error. While not explicitly stated here, we also note that Theorem 6 gives a bound on the number of trainable parameters |Θ i | that has a similar dependence on ϵ as the model in [2] . Furthermore, the parameters are independent of system size, n. Additional smoothness assumptions on the Hamiltonian H(x) can yield mild improvements on the dependence in terms of ϵ, as briefly discussed in Appendix C.1. Moreover, because of this bound on |Θ i |, applying an additional penalty on the ℓ 2 -norm of the weights Θ can help ensure that the weights remain small. In practice, this is usually satisfied during training when the weights are initialized properly and the inputs are regularized, e.g. [53] . Thus, the condition that |Θ * i | ≤ W max is often satisfied in practice and is not considered a strong assumption in deep learning.To prove the prediction error bound in Theorem 5 assuming that a low training error is achieved, we combine techniques from quasi-Monte Carlo theory applied to deep learning [54] (see Appendix A.2 for a review) along with our knowledge of the geometry of the n-qubit system. In contrast to [54] , we need to characterize the dimension of the input domain in our approach. The reason for doing this is that the approximation error depends on the size m = |I P | (Equation (2.3)) of our local models f θ P P . The central result we use here is the Koksma-Hlawka inequality [56] (see Theorem 10 in Appendix A.2) from quasi-Monte Carlo theory. This produces a bound on the prediction error in terms of the star-discrepancy (see Definition 2 in Appendix A.2) and the Hardy-Krause variation. The star-discrepancy can be controlled by known bounds on the star-discrepancy of random points [82] . We bound the Hardy-Krause variation by explicitly computing the mixed derivatives of the local models f θ P P and the ground state properties tr(Oρ(x)). In particular, we bound the latter using tools from the spectral flow formalism [59, 60, 61] , and this is where the assumption that the mixed first order derivatives of the Hamiltonian are bounded is needed. Putting these steps together, we arrive at the rigorous guarantee in Theorem 5. 
numerical experiments: 
 We conduct numerical experiments to observe the performance of our model in practice. The results demonstrate that our assumptions in Theorem 5 are often satisfied in practice and that our deep learning algorithm outperforms the previous best-known method [2] . Moreover, we generate and utilize significantly more training data than in prior works [1, 2] . The code can be found at https://github.com/marcwannerchalmers/learning_ground_states.git.We consider the classical neural network model discussed in the previous section and defined formally in Definition 6. For each of the local models f θ P P , we use fully connected deep neural networks with five hidden layers of width 200. We train the model with the AdamW optimization algorithm [83] . We measure the training error and prediction error via the root-mean-square error (RMSE). The model is discussed further in Appendix D.As in [2] , we consider the two-dimensional antiferromagnetic random Heisenberg model on between 20 to 45 qubits and predict two-body correlation functions. The corresponding Hamiltonian isEQUATIONwhere ⟨ij⟩ denotes all pairs of neighboring sites on the lattice. The coupling terms J ij correspond to the parameters x of the Hamiltonian and are sampled uniformly from [0, 2].We generate training data similarly to [1, 2] , using the density-matrix renormalization group (DMRG) [9] based on matrix-product-states (MPS) [84] . To assess the performance of our model, we consider both uniformly randomly distributed J ij and coupling parameters, which are distributed as a Sobol sequence. It is easy to see that the distributions and H satisfy the requirements of Theorem 5.In Figure 2 (Left), we see that our deep learning algorithm consistently outperforms the previous best-known ML algorithm from [2] , achieving approximately half the prediction error on the same training data. The prediction error also exhibits a constant scaling with respect to system size, agreeing with our rigorous guarantee in Theorem 5. Another noteworthy observation is that the ML algorithm's performance on LDS is nearly equivalent to its performance on uniformly random points. We discuss a potential reason for this in Appendix D. Figure 2 (Center) illustrates the prediction error scaling with respect to the training set size for various choices of δ 1 (size of the local neighborhood from Equation (2.3)). For δ 1 = 0, the error arising from approximating the ground state property via local functions dominates. For δ 1 > 0, we observe a smaller local approximation error and thus achieve a smaller prediction error for sufficiently large training sets. This is consistent with our theoretical results.Finally, Figure 2 (Right) illustrates that our assumptions in Theorem 5 are mild in practice. Namely, the blue points show that a small training error can be achieved. The red points also demonstrate that the ℓ 1 -norm of the parameters in the last layer and the largest absolute value of the parmeters in the trained neural network remain small. In particular, in Figure 2 , the weights exhibit a scaling independent of system size n. Hence, we find that the assumptions needed to guarantee the prediction error bound in Theorem 5, namely that the training objective is small and the weights of the neural network are small and independent of system size, are fulfilled in our numerical experiments. We provide further details of the numerical experiments in Appendix D. 
discussion: 
 We have shown that we can construct ML models for predicting ground state properties that require only a constant number of training samples, for a fixed prediction error. Specifically, we showed that a simple modification to the linear regression model in [2] only requires 2 polylog(ϵ -1 ) samples in order to achieve a prediction error of ϵ, provided that we know a decomposition of the observable of interest in terms of Pauli operators. We then showed that a neural network model which is trained on 2 polylog(ϵ -1 ) training samples and which achieves O(ϵ) training error on these samples will also have a prediction error of at most ϵ. In this case, knowledge of the observable O is no longer required.Our work leaves open several avenues for future exploration. First, it would be desirable to understand the conditions under which we can prove convergence for the training error. For instance, could the model be changed so as to use a convex objective, thereby avoiding the issues associated with finding a global optimum in a non-convex landscape? Following [49, 50] , we would also like to know whether the results obtained for neural networks can be extended to thermal states or Lindbladian phases of matter. Finally, for both results it would be desirable to improve the scaling with respect to the error ϵ. Currently, the models have quasipolynomial scaling in 1/ϵ and the only case in which we know how to achieve poly(1/ϵ) scaling is when the number of parameters, m, is constant (as in [51] ). These appendices provide the technical details of the ideas discussed in the main text. In Appendix A, we review several important concepts for our proofs such as the algorithm and rigorous guarantee from [2] in Appendix A.1 and background on classical deep learning techniques in Appendix A.2. In Appendix B, we build on [2] to obtain a sample complexity upper bound for predicting ground state properties independent of system size. In Appendix C, we prove our guarantee for predicting ground state properties using neural networks.A.1 Review of previous algorithm and proofIn this section, we review the previous algorithm from [2] along with intermediate results we use throughout our proofs. For full details, we refer the reader to [2] . Throughout this section, let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0. One can think of ϵ 1 as the approximation error caused by the hypothesis of our ML algorithm not exactly capturing the ground state property; ϵ 2 represents the noise in the training data; ϵ 3 corresponds to the generalization error.Recall that we consider a family of n-qubit Hamiltonians H(x) smoothly parameterized by an m-dimensional vector x ∈ [-1, 1] m that satisfies the following assumptions, which we restate from [2] .(a) Physical system: We consider n finite-dimensional quantum systems that are arranged at locations, or sites, in a d-dimensional space, e.g., a spin chain (d = 1), a square lattice (d = 2), or a cubic lattice (d = 3). Unless specified otherwise, our big-O, Ω, Θ notation is with respect to the thermodynamic limit n → ∞. (b) Hamiltonian: H(x) decomposes into a sum of geometrically local terms H(x) = L j=1 h j (⃗ x j ), each of which only acts on an O(1) number of sites in a ball of O(1) radius. Here, ⃗ x j ∈ R q , q = O(1) and x is the concatenation of L vectors ⃗ x 1 , . . . , ⃗ x L with dimension m = Lq = O(n). Individual terms h j (⃗ x j ) obey ∥h j (⃗ x j )∥ ∞ ≤ 1 and also have bounded directional derivative: ∥∂h j /∂ û∥ ∞ ≤ 1, where û is a unit vector in parameter space. (c) Ground-state subspace: We consider the ground state ρ(x) for the Hamiltonian H(x) to be defined as ρ(x) = lim β→∞ e -βH(x) / tr e -βH(x) . This is equivalent to a uniform mixture over the eigenspace of We also assume that the spectral gap of H(x) is bounded from below by some constant γ for all choices of parameters x ∈ [-1, 1] m .The ML algorithm is given a training data set {(x ℓ , y ℓ )} N ℓ=1 , where x ℓ is sampled from some distribution D over the parameter space [-1, 1] m and y ℓ approximates the ground state property:|y ℓ -tr(Oρ(x ℓ ))| ≤ ϵ 2 .The goal is to learn some function h * (x) that achieves a low average prediction errorE x∼D |h * (x) -tr(Oρ(x))| 2 ≤ ϵ. (A.1)The ML algorithm proposed in [2] requires several geometric definitions. We use S (geo) to denote the set of all geometrically local Pauli observables throughout.Let δ 1 , δ 2 , B > 0 be efficiently-computable hyperparameters that we define later. Then, define the set I P of coordinates c that parameterize some local term h j(c) that is close to a Pauli P ∈ {I, X, Y, Z} ⊗n . Here, the distance between two observables d obs is defined as the minimum distance between the qubits that the observables act on, where the distance between qubits is given by the geometry of the system, which we assume to be known. Formally, we define the set of local coordinates as I P ≜ {c ∈ {1, . . . , m} : d obs (h j(c) , P ) ≤ δ 1 }, (A.2) where h j(c) is the local term in the Hamiltonian H(x) whose parameters ⃗ x j(c) include the variable x c . The intuition behind this set of coordinates is that it indexes the parameters x c that influence the ground state property tr(P ρ(x)) corresponding to this Pauli P . Using this intuition, because these parameters x c for c ∈ I P matter most for the property we are trying to learn (as [2] proves and we give the ideas for later), then we can define a new effective parameter space in which all other parameters are set to zero. Moreover, parameters x c for c ∈ I P can be discretized to lie on a lattice. This gives the following set X PX P ≜ x ∈ [-1, 1] m : if c / ∈ I P , x c = 0 if c ∈ I P , x c ∈ {0, ±δ 2 , ±2δ 2 , . . . , ±1} . (A.3)We can also define a set T x,P for each vector x ∈ X P which is the set of parameters x ′ that are close to x for coordinates in I P :T x,P ≜ x ′ ∈ [-1, 1] m : - δ 2 2 < x c -x ′ c ≤ δ 2 2 , ∀c ∈ I P . (A.4)With these definitions in place, we set the hyperparameters as follows. Define δ 1 asδ 1 ≜ max C max log 2 (2C/ϵ 1 ), C 4 , C 5 , max(5900, α, 7(d + 11), θ) b , (A.5)where b, C max , C 4 , C 5 , α, θ, C are all constants. We refer to the supplementary information of [2] for a full description of these constants. Moreover, δ 2 is given byδ 2 ≜ 1 2 √ C ′ |I P | ϵ1 , (A.6)where C ′ is a constant. Finally, we define an additional hyperparameter B > 0 asB ≜ 2 O(polylog(1/ϵ1)) . (A.7)The ML algorithm from [2] utilizes these objects to encode the geometric locality of the system. The algorithm consists of two steps. First, it maps the parameter space[-1, 1] m to a high dimensional space R m ϕ for m ϕ ≜ P ∈S (geo) |X P | = O(n)2 O(polylog(1/ϵ1)) (A.8)via a nonlinear feature map ϕ. Second, it runs ℓ 1 -regularized linear regression (LASSO) over the feature space.This first step encodes the geometry of the problem. In particular, the feature map is defined as follows, where each coordinate of ϕ(x) is indexed by x ′ ∈ X P and P ∈ S (geo)ϕ(x) x ′ ,P ≜ 1[x ∈ T x ′ ,P ]. (A.9)In this way, the feature map ϕ(x) identifies the nearest lattice point to x. The idea is that one can approximate the ground state property well by only approximating it at these representative points and summing up. We make this intuition rigorous in the following section.Following the feature mapping, our ML algorithm uses LASSO [62, 63, 64] to learn functions of the form {h(x) = w • ϕ(x) : ∥w∥ 1 ≤ B}. In particular, for a chosen hyperparameter B > 0, LASSO finds a coefficient vector w * that solves the following optimization problem minimizing the training error subject to the constraint that ∥w∥ 1 ≤ B minw∈R m ϕ ∥w∥ 1 ≤B 1 N N ℓ=1 |w • ϕ(x ℓ ) -y ℓ | 2 . (A.10)We denote the learned function by h * (x) = w * • ϕ(x). Note that the learned function does not need to achieve the minimum training error, but can be some amount ϵ 3 /2 above it. For our purposes, we set B = 2 O(polylog(1/ϵ1)) .This algorithm obtains the following rigorous guarantee. Theorem 7 (Theorem 5 in [2] ).Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0 and δ > 0. Given training data {(x ℓ , y ℓ )} N ℓ=1 of size N = log(n/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) , (A.11) where x ℓ is sampled from D and y ℓ is an estimator of tr(Oρ(x ℓ )) such that |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 , the ML algorithm can produce h * (x) that achieves prediction errorE x∼D |h * (x) -tr(Oρ(x))| 2 ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 (A.12)with probability at least 1δ. The training time for constructing the hypothesis function h and the prediction time for computing h * (x) are upper bounded by O(nN) = n log(n/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) .This rigorous guarantee is proven by first showing that the training errorR(h) = min w 1 N N ℓ=1 |h(x ℓ ) -y ℓ | 2 (A.13) is small.Lemma 1 (Lemma 15 in [2] ). The functiong(x) = P ∈S (geo) x ′ ∈X P f P (x ′ )1[x ∈ T x ′ ,P ] = w ′ • ϕ(x), (A.14) achieves training error R(g) ≤ (ϵ 1 + ϵ 2 ) 2 . (A.15)The proof of this consists of three different steps. First, one can show that tr(Oρ(x)) can be approximated by a sum of smooth local functions, denoted as f (x) = P ∈S (geo) f P (x), where f P (x) = α P tr(P ρ(χ P (x))) for O = P ∈{I,X,Y,Z} ⊗n α P P andχ P (x) c = x c , c ∈ I P 0 c / ∈ I P (A.16)for all c ∈ {1, . . . , m}. In other words, parameters that parameterize local terms h j far away a Pauli P (x c for c / ∈ I P ) do not contribute much to the ground state property, and thus we can simply set them to zero. Formally, this approximation is given in the following lemma. Lemma 2 (Corollary 2 in [2] ). Consider a class of local Hamiltonians {H(x) : x ∈ [-1, 1] m } and an observable O = P ∈{I,X,Y,Z} ⊗n α P P satisfying assumptions (a)-(d). There exists a constant C > 0 such that for any 1/e > ϵ 1 > 0,| tr(Oρ(x)) -f (x)| ≤ Cϵ 1   P ∈S (geo) |α P |   , (A.17)where f (x) = P ∈S (geo) f P (x).Second, one can also show that this sum of local functions f (x) = P ∈S (geo) f P (x) can in turn be approximated by a linear function over the feature space g(x) = w ′ • ϕ(x), where w ′ is a vector with entries indexed by P ∈ S (geo) and x ′ ∈ X P given by w ′ x ′ ,P = f P (x ′ ). Lemma 3 (Corollary 3 in [2] ). For g(x) = w ′ • ϕ(x) and f (x) = P ∈S (geo) f P (x), then writing an observable O = P ∈{I,X,Y,Z} ⊗n α P P , we have|g(x) -f (x)| < ϵ 1   P ∈S (geo) |α P |   (A.18)for any x.This tells us that the hypothesis functions of the ML algorithm indeed approximate the ground state properties well. The final piece needed is a norm inequality bounding the ℓ 1 -norm of the Pauli coefficients. This allows us to bound the terms involving |α P | in Lemma 2 and Lemma 3. In particular, we have the following bound. Theorem 8 (Corollary 4 in [2] ). Let O = P ∈{I,X,Y,Z} ⊗n α P P be an observable that can be written as a sum of geometrically local observables. Then,P |α P | = O(1). (A.19)Given these results, Lemma 1 follows directly by triangle inequality and rescaling ϵ 1 when using Lemma 2 and Lemma 3. Finally, to prove Theorem 7, it remains to bound the generalization error by ϵ 3 . This follows directly from known sample complexity guarantees for the LASSO algorithm [64] , which learns ℓ 1 -regularized linear functions. In order to apply this known result, one needs to provide a regularization parameter, i.e., some B > 0 such that the ML algorithm learns functions of the formh(x) = w • ϕ(x) for ∥w∥ 1 ≤ B.To choose such a B, Lewis et al. bound the ℓ 1 -norm of w ′ , where recall w ′ x ′ ,P = f P (x ′ ). Lemma 4 (Lemma 14 in [2] ). Let w ′ be the vector of coefficients defined by w ′x ′ ,P = f P (x ′ ). Then,∥w ′ ∥ 1 = P ∈S (geo) x ′ ∈X P |f P (x ′ )| = 2 O(polylog(1/ϵ1)) . (A.20) Using B = 2 O(polylog(1/ϵ1)) in the known guarantees for LASSO [64] gives Theorem 7.In this section, we review results in classical deep learning theory for obtaining rigorous guarantees when learning from data sampled according to a low-discrepancy sequence (LDS) [55, 56, 76, 77, 78, 57, 79, 58, 80, 81] . For this discussion, we follow [54, 85] .We consider a neural network or multi-layer perceptron model as a composition of several layers of affine transformations and nonlinear activation functions. Namely, let σ : R → R be a nonlinear activation function. Then, a neural network with L layers is defined as follows. Let d 0 , . . . , d L ∈ N be the dimension (number of neurons or width) of each layer k ∈ {0, . . . , L}. Here, the zeroth layer is the input layer and the Lth layer is the output layer. At each layer k ∈ {0, . . . , L -1} except for the output layer, we define an affine transformationW k : R d k → R d k+1 by W k (x) = A k x + b k for a matrix of weights A k ∈ R d k+1 ×d k and a vector of biases b k ∈ R d k+1 .Then, a neural network is defined asf θ (x) = (W L-1 • σ • • • • • σ • W 0 )(x), (A.21) where σ is applied element-wise and θ = ((A 0 , b 0 ), . . . , (A L-1 , b L-1 )). The hidden layers are the first L -1 layers. Here, θ are the trainable parameters of the neural network, which can be iteratively updated through training on data. A deep neural network is a neural network with at least three layers: L ≥ 3. In this work, we consider the activation functionσ(x) = tanh(x) = e x -e -x e x + e -x . (A.22)We refer to such neural networks with this activation function as tanh neural networks.Suppose a neural network f θ aims to approximate some target function f , given training data{(x ℓ , f (x ℓ )} N ℓ=1 .Then, the training error is defined asR(θ) = 1 N N ℓ=1 |f (x ℓ ) -f θ (x ℓ )| 2 . (A.23)The prediction error is then defined over the whole domain, including unseen data, asR(θ) = E x∼D |f (x) -f θ (x)| 2 , (A.24)where the training data is sampled from some distribution D. A canonical result in deep learning theory [66] is that the generalization error can be bounded by roughlyR(θ) ≲ R(θ) + O 1 √ N , (A.25)where ≲ indicates that we only state this result schematically. Importantly, this means that in order for the neural network f θ to approximate f with high accuracy, many training data points N are needed, which is undesirable. In order to fix this issue, [54] combines ideas from deep learning with tools from quasi-Monte Carlo methods [55, 56, 57, 58] to achieve a generalization error bound ofR(θ) ≲ R(θ) + Õ 1 N , (A.26)where Õ indicates that we are suppressing polylogarithmic factors. The key tool used here is lowdiscrepancy sequences [55, 56, 76, 77, 78, 57, 79, 58, 80, 81] . Intuitively, this is a collection of points that covers domain of the function f in such a way that there are no large gaps, or discrepancies. By filling these gaps, one can ensure that the training data accurately represents the target function, more so than even uniformly random data. We leverage these ideas to obtain our rigorous guarantee on the sample complexity of a deep learning algorithm for predicting ground state properties. In the following, we formally define low-discrepancy sequences and a key inequality in quasi-Monte Carlo theory for obtaining our generalization bound.First, we define the discrepancy of a sequence, which is a measure of uniformity. Definition 1 (Discrepancy [56] ). Let λ be the Lebesgue measure, N ∈ N. Let x = {x ℓ } N ℓ=1 be a sequence of points with x ℓ ∈ [0, 1] d for all ℓ. The discrepancy of the sequence x is defined asD N (d) = sup J∈E |R N (J)|, (A.27)whereR N (J) = 1 N N ℓ=1 1{x ℓ ∈ J} -λ(J) (A.28) for a Lebesgue-measurable set J ⊆ [0, 1] d . Also, E is the set of all rectangular subsets of [0, 1] d , i.e., E = d i=1 [a i , b i ) : 0 ≤ a i < b i ≤ 1 . (A.29)Intuitively, one can consider the discrepancy as a measure of how well the sequence fills rectangular subsets of [0, 1] d . If the discrepancy is small, this means that the sequence fills these subsets well. We can similarly define the star-discrepancy, where the supremum is instead taken over rectangular subsets of [0, 1] d such that one endpoint is 0. Definition 2 (Star-discrepancy [56] ). Let λ be the Lebesgue measure, N ∈ N. Let x = {x ℓ } N ℓ=1 be a sequence of points with x ℓ ∈ [0, 1] d for all ℓ. The star-discrepancy of the sequence x is defined asD * N (d) = sup J∈E * |R N (J)|, (A.30)where R N is defined in Equation (A.28) for a Lebesgue-measurable set J ⊆ [0, 1] d . Also, E * is the set of all rectangular subsets of [0, 1] d , i.e.,E * = d i=1 [0, b i ) : 0 < b i ≤ 1 . (A.31)With these definitions, we can define low-discrepancy sequences. Definition 3 (Low-discrepancy sequence [56] ). A sequence of pointsx = {x ℓ } N ℓ=1 with x ℓ ∈ [0, 1] d for all ℓ is a low-discrepancy sequence if D * N (d) ≤ C (log N ) d N , (A.32)where C is a constant that possibly depends on d but is independent of N .The value of the constant C in this definition depends on the construction of the low-discrepancy sequence. Several constructions of low-discrepancy sequences are known [80, 77, 81, 57] . In this work, we consider Sobol sequences in base 2 [57] . For these sequences, we have the following guarantee Theorem 9 (Theorem 4.17 in [57] ).Let N ∈ N. If x = {x ℓ } N ℓ=1 is a Sobol sequence in base 2 with x ℓ ∈ [0, 1] d for all ℓ, then the star-discrepancy satisfies D * N (d) ≤ C(d) (log N ) d N , (A.33)where C(d) is a constant satisfyingC(d) < 1 d! d log(2d) . (A.34)We state this result without proof and refer to [57] for details on this construction. Another important known discrepancy bound that we will use in Appendix C.3 is the following bound on the stardiscrepancy of uniformly random points. Lemma 5 (Corollary 1 in [82] ). For any d ≥ 1, N ≥ 1 and δ ∈ (0, 1) a (uniformly) randomly generated d-dimensional point set (x 1 , . . . , x N ) satisfiesD * N (d) ≤ 5.7 4.9 + log 1 δ √ d √ N (A.35)with probability at least 1δ.As discussed earlier, low-discrepancy sequences allow us to obtain better sample complexity guarantees for neural networks. The key result in quasi-Monte Carlo theory that enables this is the Koksma-Hlawka inequality [56] . In order to properly state it, we first need to define the Hardy-Krause variation. A full technical definition can be found in, e.g., [54] , but for our purposes, it suffices to consider the following upper bound [86] . Let f be a "sufficiently smooth" function. Then, its Hardy-Krause variation can be upper bounded byV HK (f ) ≤ VHK = [0,1] d ∂ d f (y) ∂y i • • • ∂y d dy + d i=1 VHK (f (i) 1 ), (A.36)wheref (i)1 is the restriction of the function f to the boundary y i = 1. If all of the mixed partial derivatives are continuous, then this inequality is actually an equality [56] . Now, we can state the Koksma-Hlawka inequality. Theorem 10 (Koksma-Hlawka inequality). Let f : [0, 1] d → R be a function whose mixed derivatives are absolutely integrable over its domain with bounded Hardy-Krause variationV HK (f ) < ∞. Let x = {x ℓ } N ℓ=1 be a sequence of N d-dimensional points in [0, 1] d with star-discrepancy D * N (d). Then [0,1] d f (x) dx - 1 N N ℓ=1 f (x ℓ ) ≤ V HK (f )D * N (d). (A.37)This theorem is used in quasi-Monte Carlo methods to estimate the error of approximating an integral of a function f by the empirical average of f evaluated on a sequence of points. Notice that if the sequence x is a low-discrepancy sequence, then by definition, we can upper bound the star-discrepancy D * N . Moreover, recalling the definitions of prediction error and training error (Equations (A.23) and (A.24)), one can see how this relates to our task of bounding the prediction error.To generalize our results to a wider class of distributions, we need to extend these tools for arbitrary measures, rather than just the Lebesgue measure. First, we restate the definition of discrepancy and star-discrepancy [87] . Definition 4 (General Discrepancy [88] ). Let µ be a normalized Borel measure on [0, 1] d . Let x = {x ℓ } N ℓ=1 be a sequence of points with x ℓ ∈ [0, 1] d for all ℓ. The discrepancy with respect to µ of the sequence x is defined asD N (d; µ) = sup J∈E |R N (J; µ)|, (A.38)whereR N (J; µ) = 1 N N ℓ=1 1{x ℓ ∈ J} -µ(J) (A.39)for a Borel-measurable set J ⊆ [0, 1] d . Also, E is the set of all rectangular subsets of [0, 1] d , i.e.,E = d i=1 [a i , b i ) : 0 ≤ a i < b i ≤ 1 . (A.40)Definition 5 (General Star-Discrepancy [87] ). Let µ be a normalized Borel measure on [0, 1] d , and let N ∈ N. Let x = {x ℓ } N ℓ=1 be a sequence of points with x ℓ ∈ [0, 1] d for all ℓ. The star-discrepancy with respect to µ of the sequence x is defined asD * N (d; µ) = sup J∈E * |R N (J; µ)| , (A.41)where R N is defined in Equation (A.39) for a Borel-measurable set J ⊆ [0, 1] d . Also, E * is the set of all rectangular subsets of [0, 1] d , i.e.,E * = d i=1 [0, b i ) : 0 < b i ≤ 1 . (A.42)These definitions coincide with Definition 1 and Definition 2 when µ is the Lebesgue measure λ. Moreover, we can define general low-discrepancy sequences similarly to Definition 3 with respect this general star-discrepancy. There is also a generalized Koksma-Hlawka inequality [87] , which we state below.Theorem 11 (Generalized Koksma-Hlawka inequality; Theorem 1 in [87] ). Let f : [0, 1] d → R be a measurable function whose mixed derivatives are absolutely integrable over its domain with bounded Hardy-Krause variation V HK (f ) < ∞. Let µ be a normalized Borel measure on [0, 1] d , and letx = {x ℓ } N ℓ=1 be a sequence of N d-dimensional points in [0, 1] d with general star-discrepancy D * N (d; µ). Then, 1 N N ℓ=1 f (x ℓ ) - [0,1] d f (x) dµ(x) ≤ V HK (f )D * N (d; µ). (A.43)In this section, we show that with a simple modification of the algorithm from [2] , we can reduce the sample complexity to O(1) for a constant prediction error. We consider all of the same definitions/notation as in Appendix A.1. This section is similar to Section IV in the Supplementary Information of [2] . As in [2] , our algorithm first maps the parameter space [-1, 1] m into a high-dimensional feature space R m ϕ for m ϕ given in Equation (A.8) via a feature map ϕ. Our simple modification is to use the feature map defined byφ(x) x ′ ,P ≜ sign(α P ) |α P |1{x ∈ T x ′ ,P }, (B.1)where each coordinate of ϕ(x) is indexed by P ∈ S (geo) , x ′ ∈ X P . Note that defining the feature map in this way requires knowledge of the observable O = P α P P corresponding to the ground state property to be predicted. However, in practice, this is a natural assumption. The hypothesis class for our proposed ML algorithm consists of linear functions in this feature space, i.e., functions of the form h(x) = w • ϕ(x). Then, our algorithm learns these functions via ridge regression [65, 66] . For a chosen hyperparameter Λ > 0, ridge regression finds a vector w * that solves the following optimization problem minimizing the training error subject to the constraint that ∥w∥ 2 ≤ Λ minw∈R m ϕ ∥w∥ 2 ≤Λ 1 N N ℓ=1 |w • φ(x ℓ ) -y ℓ | 2 , (B.2)where y ℓ approximates tr(Oρ(x ℓ )). We denote the learned function by h * (x) = w * • φ(x). Note that the learned function does not need to achieve the minimum training error, but can be some amount say ϵ 3 /2 above it. For our purposes, we choose the hyperparameter to be Λ = 2 O(polylog(1/ϵ1)) , which we justify in the next section.Note that there are two main differences from the algorithm in [2] . First, recall from Appendix A.1 that the feature map was previously defined as ϕ(x) x ′ ,P = 1{x ∈ T x ′ ,P } for P ∈ S (geo) , x ′ ∈ X P . Second, instead of using LASSO (ℓ 1 -regularized regression), our proposed algorithm uses ridge regression.With this algorithm, we obtain the following guarantee.Theorem 12 (Constant sample complexity; Detailed restatement of Theorem 4).Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0 and δ > 0. Given training data {(x ℓ , y ℓ )} N ℓ=1 of size N = log(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) , (B.3)where x ℓ is sampled from D and y ℓ is an estimator of tr(Oρ(x ℓ )) such that |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 , the ML algorithm can produce h * (x) that achieves prediction errorE x∼D |h * (x) -tr(Oρ(x))| 2 ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 (B.4)with probability at least 1δ. The training time for constructing the hypothesis function h * and the prediction time for computing h * (x) are upper bounded byO(n)polylog(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) . (B.5)Comparing to Theorem 7, notice that our sample complexity guarantee is completely independent of system size n.The theorem in the main text corresponds to ϵ 1 = 0.2ϵ, ϵ 2 = ϵ, and ϵ 3 = 0.4ϵ. In this way,(ϵ 1 + ϵ 2 ) 2 ≤ 1.44ϵ 2 ≤ 0.53ϵ and (ϵ 1 + ϵ 2 ) 2 + ϵ 3 ≤ ϵ.So far, we have only considered the setting in which we learn a specific ground state property tr(Oρ(x)) for a fixed observable O. Because our training data is given in the form {(x ℓ , y ℓ )} N ℓ=1 , where y ℓ approximates tr(Oρ(x)) for this fixed observable O, if we want to predict a new property for the same ground state ρ(x), we would need to generate new training data. Thus, it may be more useful to learn a ground state representation, from which we could predict tr(Oρ(x)) for many different choices of observables O without requiring new training data. In this case, suppose we are instead given training data {x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 , where σ T (ρ(x ℓ )) is a classical shadow representation [52, 68, 69, 70, 71] of the ground state ρ(x ℓ ). An immediate corollary of Theorem 12 is that we can predict ground state representations with the same sample complexity. This follows from the same proof as Corollary 5 in [2] . Corollary 3 (Learning representations of ground states; detailed restatement of Corollary 1). Let1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0 and δ > 0. Given training data {(x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 of size N = log(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) , (B.6)where x ℓ is sampled from D and σ T (ρ(x ℓ ) is the classical shadow representation of the ground state ρ(x ℓ ) using T randomized Pauli measurements. For T = O(log(nN/δ)/ϵ 2 2 ) = Õ(log(n/δ)/ϵ 2 2 ), the ML algorithm can produce a ground state representation ρN,T (x) that achievesE x∼D | tr(O ρN,T (x)) -tr(Oρ(x))| 2 ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 (B.7)with probability at least 1δ, for any observable with eigenvalues between -1 and 1 that can be written as a sum of geometrically local observables.We note that the number of measurements T needed to generate the training data scales as log(n), but the amount of training data still remains constant with respect to system size. We do not consider the number of measurements as contributing to the sample complexity because in our setting, the ML algorithm is given this training data as input and does not generate it itself.To prove Theorem 12, we first derive a bound on the training error. Recall that the training error is defined asR(h) = min w 1 N N ℓ=1 |h(x ℓ ) -y ℓ | 2 . (B.8)Define the vector w with entries indexed by P ∈ S (geo) , x ′ ∈ X P by wx ′ ,P ≜ |α P | tr(P ρ(χ P (x))), (B.9)where χ P (x) is defined in Equation (A.16). Then, notice thatg(x) ≜ w • φ(x) (B.10) = P ∈S (geo) x ′ ∈X P sign(α P )|α P | tr(P ρ(χ P (x)))1{x ∈ T x ′ ,P } (B.11) = P ∈S (geo) x ′ ∈X P α P tr(P ρ(χ P (x)))1{x ∈ T x ′ ,P } (B.12) = w ′ • ϕ(x) (B.13) = g(x), (B.14)where g(x) = w ′ •ϕ(x) with w ′ x ′ ,P = α P tr(P ρ(χ P (x))), ϕ(x) x ′ ,P = 1{x ∈ T x ′ ,P }. By Lemma 1, we know that g(x) approximates the ground state property with low training error, and thus, in turn, g(x) also approximates the ground state property well. The existence of w such that g(x) = w • φ(x) guarantees that the function h * (x) = w * • φ(x) found by performing via ridge regression will also yield a small training error. More formally, we have the following guarantee Lemma 6 (Training error). The functiong(x) = w • φ(x) = P ∈S (geo) x ′ ∈X P α P tr(P ρ(χ P (x)))1{x ∈ T x ′ ,P } (B.15) achieves training error R(g) ≤ (ϵ 1 + ϵ 2 ) 2 . (B.16)Since g(x) = g(x), this follows directly from Lemma 1. Moreover, we can obtain an ℓ 2 -norm bound on w. We can utilize this upper bound to choose the hyperparameter Λ > 0 such that ∥w∥ 2 ≤ Λ. Thus, we have the following lemma, Lemma 7 (ℓ 2 -Norm bound). Let w be the vector of coefficients defined in Equation (B.9). Then, we have ∥ w∥2 2 = P ∈S (geo) x ′ ∈X P |α P || tr(P ρ(χ P (x)))| 2 = 2 O(polylog(1/ϵ1)) . (B.17)Proof. This is a simple consequence of Lemma 4. Explicitly, we have∥ w∥ 2 = P ∈S (geo) x ′ ∈X P |α P || tr(P ρ(χ P (x)))| 2 (B.18) ≤ P ∈S (geo) x ′ ∈X P |α P || tr(P ρ(χ P (x)))| (B.19) = 2 O(polylog(1/ϵ1)) , (B.20)where the second line follows because tr(P ρ(χ P (x))) ≤ 1 and the last line follows by Lemma 4.This justifies our choice of Λ = 2 O(polylog(1/ϵ1)) . Now consider the learned functionh * (x) = w * • φ(x), where w * is found by minimizing the training error subject to the constraint that ∥w∥ 2 ≤ Λ.We do not require the learned function to achieve the minimum training error, but it can be some amount ϵ 3 /2 above it, i.e.,R(h * ) ≤ ϵ 3 2 + min w ∥w∥ 2 ≤Λ 1 N N ℓ=1 |w • φ(x ℓ ) -y ℓ | 2 . (B.21)Since we chose Λ = 2 O(polylog(1/ϵ1)) and we showed in Lemma 7 that ∥ w∥ 2 ≤ Λ, then the minimum training error is at most R(g). We also know that this is bounded by (ϵ 1 + ϵ 2 ) 2 by Lemma 6. This then impliesR(h * ) ≤ ϵ 3 2 + R(g) ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 . (B.22)To prove Theorem 12, it remains to bound the prediction error. We can use a standard result from machine learning theory on the prediction error of ridge regression algorithms [66, 64] . Theorem 13 (Theorem 26.12 in [66] ). Suppose that D is a distribution over X × Y such that with probability 1 we have that ∥x∥ 2 ≤ R. Let H = {x → w • x : ∥w∥ 2 ≤ Λ} and let ℓ : H × Z → R be a loss function of the form ℓ(w, (x, y)) = ϕ(w • x, y) such that for all y ∈ Y, a → ϕ(a, y) is a ρ-Lipschitz function and such that max a∈[-ΛR,ΛR] |ϕ(a, y)| ≤ c. Then, for any δ ∈ (0, 1), with probability of at least 1δ over the choice of an i.i.d. sample of size N , for all h ∈ H,R(h) ≤ RS (h) + 2ρΛR √ N + c 2 log(2/δ) N . (B.23)Here, R(h) denotes the prediction error for the hypothesis h. With this, we can complete the proof of Theorem 12.Proof of Theorem 12. First, let us reframe the theorem in our setting. Consider the input space X to be the parameter space [-1, 1] m and our input variable is x = φ(x). Since the observables we consider have spectral norm at most 1, the output space fulfils Y ⊆ [-1, 1]. The hypothesis set is H = {x → w • φ(x) : ∥w∥ 2 ≤ Λ}, where in the previous section, we set Λ = 2 O(polylog(1/ϵ1)) .It remains to check the conditions of the theorem. We begin by showing that ∥x∥ 2 ≤ R for some R > 0. We have the following computation: = O(1), (B.28) where the second to last line follows because for a given P , x ∈ T x ′ ,P for exactly one x ′ ∈ X P . This is shown in Corollary 3 of [2] . Also, the last line follows by Theorem 8. Thus, we can take R = O(1).∥x∥ 2 = φ(x) • φ(x) = φ(x)Finally, note that ℓ(w, (x, y)) = |w • x -y| = ϕ(w • x, y). Therefore, ϕ(a, y) is a 1-Lipschitz function and fulfils max a∈[-ΛR,ΛR] |ϕ(a, y)| = max a∈[-ΛR,ΛR] |a -y| ≤ ΛR + 1. (B.29) Thus, we can consider ρ = 1 and c = O(1) • 2 O(polylog(1/ϵ1)) + 1.By Equation (B.22), we know that the learned modelh * (x) = w * • φ(x) achieves R(h * ) ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 . (B.30) Plugging in R = O(1), ρ = 1, Λ = 2 O(polylog(1/ϵ1)) and c = O(1) • 2 O(polylog(1/ϵ1)) + 1 into Theorem 13, we have R(h * ) ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 (B.31) + 1 √ N 2O(1) • 2 O(polylog(1/ϵ1)) + O(1) • 2 O(polylog(1/ϵ1)) + 1 2 log(2/δ) (B.32)with probability at least 1δ. In order to bound the prediction error by (ϵ 1 + ϵ 2 ) 2 + ϵ 3 , we need N to be large enough such that1 √ N 2O(1) • 2 O(polylog(1/ϵ1)) + O(1) • 2 O(polylog(1/ϵ1)) + 1 2 log(2/δ) ≤ ϵ 3 2 . (B.33)Solving for N in this inequality and simplifying we haveN = 4 ϵ 2 3 2 O(polylog(1/ϵ1)) (1 + log(1/δ)) 2 = 2 O(log(1/ϵ3)+polylog(1/ϵ1)) log(1/δ). (B.34)Thus, for this N , we can guarantee that R(h * ) ≤ (ϵ 1 + ϵ 2 ) 2 + ϵ 3 , as claimed.On another note, when considering a scenario with a fixed number of parameters m = O(1), much like the setting in [51] , the expression derived from the result in Lemma 10 exhibits polynomial dependence on ϵ. One can incorporate the constant number of parameters by setting m = m. Thus, we recover the exact ground state properties tr(P ρ(x)) in f P and the approximation error resulting from applying Lemma 2 vanishes completely. Furthermore, we can slightly adapt the proof of Lemma 7 and obtain ∥ w∥ where the last step is performed similarly as in the proof of Lemma 4.It remains to analyze the computational time for the ML algorithm's training and prediction.Proof of computational time in Theorem 12. The training time is dominated by the time required for ridge regression over the feature space defined by the feature map ϕ. Recall that the optimization problem under considerations ismin w∈R m ϕ ∥w∥ 2 ≤Λ 1 N N ℓ=1 |w • φ(x ℓ ) -y ℓ | 2 . (B.36)One can show that this is a convex optimization problem so that we can solve its equivalent dual problem instead. This dual optimization problem is given bymax α∈R N -α ⊺ (K + λI)α + 2α • Y, (B.37)where the kernel matrix is K = X ⊺ X, for the feature matrixX ∈ R m ϕ ×N defined by X = ( φ(x 1 ) • • • φ(x N ))and the response vector Y = (y 1 , . . . , y N ) ⊺ . If κ is the maximum time it takes to compute a kernel entry K(x, x ′ ) = φ(x) • φ(x ′ ), then one can show that the time to solve this dual problem is O(κN 2 + N 3 ). Moreover, prediction can be executed in O(κN ). For more details in this analysis, we refer the reader to, e.g., Section 11.3.2 of [64] .In our case, κ = O(m ϕ ) since φ(x) ∈ R m ϕ and the kernel is simply the dot product of two of these vectors. By Equation (A.8), we know thatm ϕ = O(n)2 O(polylog(1/ϵ1)) (B.38) so that κ = O(n)2 O(polylog(1/ϵ1)). Moreover, by Theorem 12,N = log(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) . (B.39)Plugging this into the time required to solve the dual problem for kernel ridge regression, we haveO(κN 2 + N 3 ) = O(n)polylog(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) . (B.40)Moreover, the prediction time is given byO(κN ) = O(n)polylog(1/δ)2 O(log(1/ϵ3)+polylog(1/ϵ1)) . (B.41)In this section, we derive a rigorous guarantee on the sample complexity of a deep-learning based model for predicting ground state properties. Similarly to the previous sections, let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0 throughout. One can think of ϵ 1 as the approximation error caused by our neural network model not exactly capturing the ground state property; ϵ 2 represents the noise in the training data; ϵ 3 corresponds to the generalization error.Recall again the setup, where we consider a family of n-qubit Hamiltonians H(x) = L j=1 h j (⃗ x j ) parameterized by an m-dimensional vector x ∈ [-1, 1] m , which satisfies the assumptions (a)-(c) in Appendix A.1. Let ρ(x) denote the ground state of H(x). We consider the task of predicting ground state properties tr(Oρ(x)) for some observable O that satisfies assumption (d) in Appendix A.1, where we are given training data {(x ℓ , y ℓ )} N ℓ=1 with y ℓ ≈ tr(Oρ(x ℓ )). In particular, suppose |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Furthermore, we also assume that all mixed partial derivatives of order m ≜ |I P | of h j are bounded as∂ m ∂x 1 ∂x 2 . . . ∂x m h j (x) ∞ ≤ 1, (C.1)where I P is the set of local coordinates defined in Equation (A.2). Here, we denote the number of local coordinates by m = |I P | for ease of notation. This is similar in spirit to assumption (b), in which we assume that the local terms have bounded directional derivatives: ∥∂h j /∂ û∥ ∞ ≤ 1,w P f θ P P (x), (C.2)where the "local models"f θ P P : [-1, 1] m → R are given by f θ P P (x) = (W out • tanh • W hidden • tanh • W in • τ -1 )(x), (C.3) with τ -1 (x) = (x + 1)/2 and θ P = [(W in , b in ), (W hidden , b hidden ), (W out , b out )]. Here, W in ∈ R m×W , b in ∈ R W , W hidden ∈ R W ×W , b hidden ∈ R W , W out ∈ R W ×1and b out ∈ R, where W denotes the width of the hidden layers. The weights are given by Θ = {θ P : P ∈ S (geo) } in the local models and w ∈ R in the last layer. Furthermore, we denote the individual parameters by Θ i ∈ R.Using this model, we can establish an objective function that we aim to minimize during the training process. Specifically, this objective function comprises the mean square error along with a lasso penalty applied to the weights w in the final layer. Definition 7 (Training objective). Let f Θ,w be a neural network model as in Definition 6. Let {(x ℓ , y ℓ )} N ℓ=1 be the training data set and λ > 0 be some regularization parameter that may depend on ϵ 1 , ϵ 2 > 0. The training objective is given by1 N N ℓ=1 |f Θ,w (x ℓ ) -y ℓ | 2 + λ∥w∥ 1 .(C.4)Our proposed ML algorithm then operates as in Algorithm 1. Initialize model architecture according to Definition 6 with appropriate hyperparameter δ 1 , width W as in Theorem 16 and weights Θ, w using an appropriate initialization method (e.g., Xavier initialization [72] ); Train with respect to the objective in Definition 7 with appropriate hyperparameter λ > 0 using a quasi-Monte Carlo training algorithm, e.g., Adam [73] until convergence; Obtain locally optimal parameters Θ * , w * ; Result: Classical representation f Θ * ,w * ;After training our model using Algorithm 1, we obtain the following rigorous guarantee. Theorem 14 (Neural network sample complexity guarantee).Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0. Let f Θ * ,w * : [-1, 1] m → R be a neural network model produced from Algorithm 1 trained on data {(x ℓ , y ℓ )} N ℓ=1 of size N = 2 O(polylog(1/ϵ1)+polylog(1/ϵ3)) , (C.5) where the x ℓ 's form a low-discrepancy Sobol sequence and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Suppose that f Θ * ,w * achieves a training error of at most ((ϵ 1 + ϵ 2 ) 2 + ϵ 3 )/2. Additionally, suppose that all parameters Θ * i of f Θ * ,w * satisfy |Θ * i | ≤ W max , for some W max > 0 that is independent of the system size n. Then the neural network f Θ * ,w * achieves prediction errorE x∼U [-1,1] m |f Θ * ,w * (x) -tr(Oρ(x))| 2 ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 , (C.6) where x ∼ U [-1, 1] m denotes sampling x from a uniform distribution over [-1, 1] m .We prove this theorem in the next two sections (Appendices C.1 and C.2). As a corollary of this, we obtain the theorem stated in the main text. We discuss the assumptions that the distribution D must satisfy in depth and prove the corollary in Appendix C.3. The theorem in the main text (Theorem 5) corresponds to ϵ 1 = ϵ 3 = 0.1ϵ and ϵ 2 = ϵ. Hence, 2(ϵ 1 + ϵ 2 ) 2 ≤ 2.44ϵ 2 ≤ 0.9ϵ for 1/e > ϵ > 0 and so 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 ≤ ϵ. Corollary 4 (Neural network sample complexity guarantee; detailed restatement of Theorem 5). Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0, D a distribution with PDF g satisfying the following properties: g has full support and is continuously differentiable on [-1, 1] m . Moreover, g is of the formg(x) = L j=1 g j (⃗ x j ). (C.7) Let f Θ * ,w * : [-1, 1] m → R be a neural network model produced from Algorithm 1 trained on data {(x ℓ , y ℓ )} N ℓ=1 of size N = log(1/δ)2 O(polylog(1/ϵ1)+polylog(1/ϵ3)) , (C.8)where the x ℓ ∼ D and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Suppose that f Θ * ,w * achieves a training error of at most ((ϵ 1 + ϵ 2 ) 2 + ϵ 3 )/2. Additionally, suppose that all parameters Θ * i of f Θ * ,w * satisfy |Θ * i | ≤ W max , for some W max > 0 that is independent of the system size n. Then the neural network f Θ * ,w * achieves prediction errorE x∼D |f Θ * ,w * (x) -tr(Oρ(x))| 2 ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 , (C.9)with probability at least 1δ.Moreover, while we can show the existence of suitable parameters that achieve a low training error, quantified by our training objective in Definition 7, in general, we cannot prove that Algorithm 1 converges to parameters close to this optimum because our training objective is not convex. Thus, to obtain the guarantee in Theorem 14, we need to assume that a low training error is indeed achieved by Algorithm 1. This is commonly satisfied in practice.Similar to Corollary 3, if we are instead given training data {x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 , where σ T (ρ(x ℓ )) is a classical shadow representation [52, 68, 69, 70, 71] of the ground state ρ(x ℓ ), then an immediate corollary of Theorem 14 is that we can predict ground state representations with the same sample complexity. This follows from the same proof as Corollary 5 in [2] . Corollary 5 (Learning representations of ground states with neural networks; detailed restatement of Corollary 2). Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0 and δ > 0. Given training data {(x ℓ , σ T (ρ(x ℓ ))} N ℓ=1 of size N = log(1/δ)2 O(polylog(1/ϵ3)+polylog(1/ϵ1)) , (C.10) where x ℓ is sampled from a distribution D satisfying the same assumptions as Corollary 4 and σ T (ρ(x ℓ ) is the classical shadow representation of the ground state ρ(x ℓ ) using T randomized Pauli measurements. For T = O(log(nN/δ)/ϵ 2 2 ) = Õ(log(n/δ)/ϵ 2 2 ), the ML algorithm can produce a ground state representation ρN,T (x) that achievesE x∼D | tr(O ρN,T (x)) -tr(Oρ(x))| 2 ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 (C.11)with probability at least 1δ, for any observable with eigenvalues between -1 and 1 that can be written as a sum of geometrically local observables.We review low-discrepancy sequences and techniques in quasi-Monte Carlo theory in Appendix A.2, which we use in our proof. To prove Theorem 14, we first show that there exists weights Θ ′ , w ′ such that our proposed neural network f Θ ′ ,w ′ achieves a low training error, i.e., it approximates the ground state properties tr(Oρ(x)) well. We show this using results in classical deep learning theory about approximating arbitrary functions with neural networks [53] . Then, we use the Koksma-Hlawka inequality (Theorem 10) to bound the prediction error of our model, similarly to [54] . As we want to derive a bound which is independent of the size of the physical system, our approach requires some additional steps. Since the dimension of the input domain of our model depends on the size of the physical system, we cannot treat it as constant as in [54] . Therefore, we bound the prediction error with respect to local functions, whose domain size is independent of the system size.Moreover, recall that the Koksma-Hlawka inequality produces a bound in terms of the star-discrepancy (Definition 2) and the Hardy-Krause variation. The star-discrepancy can be bounded by considering low-discrepancy sequences (Definition 3), and the Hardy-Krause variation can be bounded by Equation (A.36). We derive explicit bounds for the Hardy-Krause variation of the ground state properties tr(Oρ(x)), using tools from the spectral flow formalism [59, 60, 61] . To obtain Corollary 4, we follow a similar proof but use results relating the discrepancy with respect to the Lebesgue measure to the discrepancy with respect to arbitrary measures and bounds on the discrepancy of uniformly random points (Appendix A.2).In Appendix C.1, we prove that our model approximates the ground state properties well. Then, in Appendix C.2, we use the Koksma-Hlawka inequality to bound the prediction error of our model. Technical results explicitly bounding the mixed partial derivatives of the ground state properties are found in Appendix C.4. We use these in Appendix C.2 to bound the Hardy-Krause variation. We then generalize this to data sampled from different distributions, as in Corollary 4, in Appendix C.3.In this section, we prove that when choosing the number of parameters and width of the model appropriately, there exists a parameter set for which the deep neural network model approximates the ground state properties well. This shows the existence of a neural network with low training error. The proof is a direct application of the main result from [53] , which proves that tanh neural networks can approximate sufficiently smooth functions, in combination with the bounds on the mixed derivative of tr(P ρ(x)) we derived in Appendix C.4.We consider the local functions defined as in Appendix A.1.2. Namely, define f (x) = P ∈S (geo) α P f P (x), where f P (x) = tr(P ρ(χ P (x))) for O = P ∈{I,X,Y,Z} ⊗n α P P andχ P (x) c = x c , c ∈ I P 0 c / ∈ I P (C.12)for all c ∈ {1, . . . , m}, for I P defined in Equation (A.2). Note that here, we slightly alter the definition from Appendix A.1.2, where we do not include the coefficient α P in the definition of f P (x). Because all parameters with coordinates not in I P are set to 0, we can view f P as a function taking inputs in [-1, 1] m, where recall that we use m = |I P | to denote the number of local coordinates.We show that there exists a neural network that can approximate these local functions f P well.In order to apply the result from [53] to approximate tr(P ρ(χ P (x))), we need to transform its inputs, such that it becomes a map [0, 1] m → R. Therefore, we introduce an appropriate function τ (x) = 2x -1. To avoid confusion when considering the different domains[-1, 1] m versus [0, 1] m , if an input x ∈ [-1, 1] m, we simply denote it by x. If x ∈ [0, 1] m , we denote it by x. We similarly use this notation for other domain dimensions.In the following lemma, we use W k,∞ (Ω) for Ω ⊆ R m to denote the Sobolev spaceW k,∞ (Ω) ≜ f ∈ L ∞ (Ω) : ∂ |α| f ∂x α1 1 • • • ∂x αm m ∈ L ∞ (Ω) for all α ∈ N m 0 with |α| ≤ k (C.13)so that the Sobolev norm is defined as∥f ∥ W k,∞ (Ω) ≜ max |α|=s ∂ |α| f ∂x α1 1 • • • ∂x αm m L ∞ (Ω) . (C.14) for α ∈ N m 0 and the L ∞ -norm is defined by ∥f ∥ L ∞ (Ω) = sup x∈Ω ∥f ∥. (C.15) Lemma 8 (Existence of approximating neural network). Let ϵ 1 > 0, s, M ∈ N. Let ∥H(x)∥ W s,∞ ([-1,1] m ) ≤ 1. Define functions f P : [0, 1] m → R as f P (τ (x)) = tr(P ρ(χ P (τ (x)))),where τ (x) = 2x -1. Then, there exist neural networks f M P , such thatf P • τ -f M P L ∞ ([0,1] m ) ≤ ϵ 1 (C.16)with at most ϵ-m+1 s 12 O( m log( m)) parameters. Furthermore, the weights scale as 2 poly(log(1/ϵ1), m,s) .To prove this, we utilize the main result from [53] , which states that a neural network with tanh activation functions can approximate effectively any function.Theorem 15 (Theorem 5.1 in [53] ). Let d, s ∈ N, R > 0, d > 0 and f ∈ W s,∞ ([0, 1] d ). There exist constants C(d, k, s, f ), N 0 (d) > 0, such that for every N ∈ N with N > N 0 (d) there exists a tanh neural network f N with two hidden layers, one of width at most 3⌈ s 2 ⌉|P s-1,d+1 | + d(N -1) (where|P n,d | = n+d-1 n) and another of width at most 3⌈ d+2 2 ⌉|P d+1,d+1 |N d (or 3⌈ s 2 ⌉ + N -1 and 6N for d = 1), such that,∥f -f N ∥ L ∞ ([0,1] d ) ≤ (1 + δ) C(d, 0, s, f ) N s , (C.17)and for k = 1, . . . , s -1,∥f -f N ∥ W k,∞ ([0,1] d ) ≤ 3 d (1 + δ)(2(k + 1)) 3k max R k , ln k βN s+d+2 C(d, k, s, f ) N s-k , (C.18)where we defineβ = k 3 2 d √ d max{1, ∥f ∥ 1/2 W k,∞ ([0,1] d ) } δ min{1, C(d, k, s, f )} . (C.19) If f ∈ C s ([0, 1] d ), then it holds that C(d, k, s, f ) = max 0≤l≤k 1 (s -l)! 3d 2 s-l |f | W s,∞ ([0,1] d ) , N 0 (d) = 3d 2 , (C.20)and else it holds thatC(d, k, s, f ) = max 0≤l≤k π 1/4 √ s (s -l -1)! 5d 2 s-l |f | W s,∞ ([0,1] d ) , N 0 (d) = 5d 2 . (C.21)In addition, the weights of f N scale as O C -s/2 N d(d+s 2 +k 2 )/2 (s(s + 2)) 3s(s+2) .The proof of Lemma 8 follows by an application of Theorem 15.Proof of Lemma 8. We directly apply Theorem 15, with the input space dimension m, where m is the number of local parameters m = |I P |. By Corollary 8, then we havef P • τ -f M P L ∞ ([0,1] m) ≤ (1 + δ) s! 3 mCs 2 2M s . (C.22)We want to show that this is bounded above by ϵ 1 . By rearranging, we find that this holds when and therefore at mostϵ -1 s 1 (1 + δ) s! 1 s 3 2 mCs 2 ≤ M. (C.(c 1 m) c2 mϵ -m+1 s 1 = 2 O( m(log( m)+log(1/ϵ1)/s)) (C.25)trainable weights in the neural network. The constants c 1 and c 2 are independent of m, but may depend on s. By Theorem 15, the weights scale asO C -s/2 ϵ -1 s 1 ms 2 m( m+s 2 +k 2 )/2 (s(s + 2)) 3s(s+2) = ϵ -m+1 s 1 2 O( m log( m)) . (C.26)Although the dependence on s is not relevant for our final statement, it is important to comment on the effect of the smoothness of H(x). The result states that the dependence of the required parameters with respect to 1/ϵ 1 improves with the highest degree for which all mixed derivatives of H(x) are bounded. When H(x) is analytic, s can be chosen to be very large so that the number of parameters in the neural network almost scales as O( ms log( m) ). The constant scales rather poorly with s; however, this effect is only be visible for very small ϵ 1 .Using Lemma 8, we can show that there exist parameters such that the complete model approximates tr(Oρ(x)) and obtains a small training objective (defined in Definition 7) . The theorem (Theorem 6) in the main text corresponds to ϵ 1 = 0.2ϵ, ϵ 2 = ϵ so that (ϵ 1 + ϵ 2 ) 2 ≤ 1.44ϵ 2 ≤ 0.53ϵ ≤ ϵ. Theorem 16 (Detailed restatement of Theorem 6). For any 1/e > ϵ 1 , ϵ 2 > 0 and appropriate width W , there exist weights Θ ′ , w ′ such that the neural network modelf Θ ′ ,w ′ satisfies |f Θ ′ ,w ′ (x) -tr(Oρ(x))| ≤ ϵ 1 (C.27)for any x ∈ [-1, 1] m . In particular, for any collection of N training data points {(x ℓ , y ℓ )} N ℓ=1 with |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 , we have1 N N ℓ=1 |f Θ ′ ,w ′ (x ℓ ) -y ℓ | 2 + λ∥w ′ ∥ 1 ≤ (ϵ 1 + ϵ 2 ) 2 (C.28)for a suitable choice of regularization parameter λ = O(ϵ 2 1 ). Moreover, each parameter Θ i of the network has a magnitude of|Θ i | = 2 O(polylog(1/ϵ1)) .Proof. Write O = P α P tr(P ρ(x)). By Theorem 8, let D = O(1) be a constant such thatP |α P | ≤ D. (C.29)For every Pauli P , then by Lemma 8, there exist weights θ ′ P such that a neural network f θ ′ P : [0, 1] m → R with two hidden layers as in Definition 6 approximates the local functions f P (τ (x)) = tr(P ρ(χ P (τ (x)))), where τ (x) = 2x -1, up to an error ϵ 1 /(4D), when the width of their hidden layers is chosen as W = ϵ where x ∈ [0, 1] m. Because τ is simply a coordinate transformation, then we also obtainf θ ′ P P (x) -tr(P ρ(χ P (x))) ≤ ϵ 1 4D , (C.31) for x ∈ [-1, 1] m.Furthermore, by Lemma 2 in Appendix A.1.2, the sum of the local functions f (x) = P α P f P (x) approximates the ground state property tr(Oρ(x)) = P α P tr(P ρ(x)) well. In particular, combining Lemma 2 with Theorem 8, we have P α P tr(P ρ(χ P (x))) -P α P tr(P ρ(x)) ≤ ϵ 1 4 . (C.32)This holds when choosing the local radius δ 1 (defined in Equation (A.5)) to be δ 1 = 4C log 2 (1/ϵ 1 ) for some constant C. This implies that m = |I P | = O(polylog(1/ϵ 1 )) (e.g., Equation (S33) of [2] ). Thus, for the model f 1 N N ℓ=1 |f Θ ′ ,w ′ (x ℓ ) -y ℓ | 2 + λ∥w ′ ∥ 1 (C.38) = 1 N N ℓ=1 |f Θ ′ ,w ′ (x ℓ ) -tr(Oρ(x ℓ )) + tr(Oρ(x ℓ )) -y ℓ | 2 + λ∥w ′ ∥ 1 (C.39) ≤ (ϵ 1 /2 + ϵ 2 ) 2 + λ∥w ′ ∥ 1 (C.40) ≤ ϵ 1 2 + ϵ 2 2 + ϵ 2 1 2 (C.41) ≤ (ϵ 1 + ϵ 2 ) 2 . (C.42)Finally, plugging in m = O(polylog(1/ϵ 1 )), by Lemma 8, then we obtain |Θ i | = 2 O(polylog(1/ϵ1)) , as required.In this section, we derive our result on the prediction error to complete the proof of Theorem 14.The central result we use is the Koksma-Hlawka inequality (Theorem 10) from quasi-Monte Carlo theory, which produces a bound in terms of the star-discrepancy (Definition 2) and the Hardy-Krause variation (Equation (A.36)). We review these tools in Appendix A.2. To bound the star-discrepancy, we consider a specific low-discrepancy sequence with guarantees described in Appendix A.2. The Hardy-Krause variation can be bounded by considering the mixed derivatives of our target function tr(Oρ(x)) and our neural network model. We relegate the bounds on the mixed derivatives of tr(Oρ(x)) to Appendix C.4, as the discussion is fairly technical. To bound the mixed derivatives of our model, we consider the following lemma.Lemma 9 (Bound on mixed derivatives of neural network). Let k, d ∈ N. Let f : [-1, 1] d → R be a tanh neural network with two hidden layers of width W ≥ d and maximal weight W max . Then∥ f ∥ W k,∞ ([-1,1] d ) = 2 O(k 2 log(k)+k log(dW Wmax)) . (C.43)Proof. Recall that a tanh deep neural network with two hidden layers is defined as a functionf : [-1, 1] d → R such that f (x) = (W out • tanh • W hidden • tanh • W in )(x), (C.44)where the activation function tanh is applied element-wise. Note that this result holds for any tanh neural network with two hidden layers, where this neural network does not necessarily have to be the same model as Definition 6.Let W L ∈ {W in , W hidden , W out } denote the layers of the neural network that perform an affine transformation for L ∈ {in, hidden, out}. We can also use L ∈ {0, 1, 2}, where 0 corresponds to in, 1 corresponds to hidden, and 2 corresponds to out. Let d L denote the width (number of input neurons) in each layer, where we defined 0 = d in , d 1 = d 2 = W, d 3 = d out . In this way, W L : R d L → R d L+1 for L ∈ {0, 1, 2}. Explicitly, we have W in : R din → R W , W hidden : R W → R W , W out : R W → R dout .Since W L performs an affine transformation, we can write it has W L (x) = (f 1 (x), . . . , f d L+1 (x)), where x ∈ R d L and f i are linear functionsf i (x) = w ⊺ i x + b i for w i ∈ R d L , b i ∈ R.For these linear layers, we observe for any function g : R dg → R d L with input dimension d g and for L ∈ {0, 1, 2}, we havemax 1≤i≤d L+1 ∥(W L • g) i ∥ W k,∞ ([-1,1] d ) = max 1≤i≤d L+1 ∥f i (g(x))∥ W k,∞ ([-1,1] d ) (C.45) ≤ d L+1 i=1 ∥w T i g(x) + b i ∥ W k,∞ ([-1,1] d ) (C.46) ≤ max 1≤j≤dg ∥W L ∥ 1 ∥g(x) j ∥ W k,∞ ([-1,1] d ) , (C.47)where we use the notation∥W L ∥ 1 ≜ d L+1 i=1   |b i | + d L j=1 |w i,j |   , (C.48)where w is a matrix with rows given by the vectors w ⊺ i , w i ∈ R d L . To show this inequality, we used Hölder's inequality in the last step. With this, by factoring out one layer W L at a time, we can bound the Sobolev norm of f . In particular, we have∥ f ∥ W k,∞ ([-1,1] d ) (C.49) = ∥(W out • tanh • W hidden • tanh • W in )(x)∥ W k,∞ ([-1,1] d ) (C.50) ≤ ∥W out ∥ 1 max 1≤j≤W ∥(tanh • W hidden • tanh • W in ) j ∥ W k,∞ ([-1,1] d ) (C.51) ≤ ∥W out ∥ 1 16(e 2 k 4 d 2 ) k (2k) k(k+1) max 1≤j≤W ∥(W hidden • tanh • W in ) j ∥ k W k,∞ ([-1,1] d ) (C.52) ≤ ∥W out ∥ 1 ∥W hidden ∥ k 1 • 16(e 2 k 4 d 2 ) k (2k) k(k+1) max 1≤j≤W ∥(tanh • W in ) j ∥ k W k,∞ ([-1,1] d ) (C.53) ≤ ∥W out ∥ 1 ∥W hidden ∥ k 1 • 16 2 (e 2 k 4 d 2 ) 2k (2k) 2k(k+1) ∥W in ∥ k 1 . (C.54)In the second line, we used Equation (C.47). In the third line, we used the two following inequalities:d m dx m tanh(x) ≤ (2m) m+1 min{exp(-2x), exp(2x)} (C.55)for all x ∈ R, m ∈ N (see Lemma A.4 in [53] ), and [53] ). In the fourth and fifth lines, we used Equation (C.47) and these inequalities again. Furthermore, we used that our inputs are absolutely bounded by 1 in the last step.∥g • f ∥ W n,∞ ≤ 16(e 2 n 4 md 2 ) n ∥g∥ W n,∞ max 1≤i≤m ∥(f ) i ∥ n W n,∞ (C.56) for any functions f ∈ C n (Ω 1 ; Ω 2 ) and g ∈ C n (Ω 2 ; R) defined on Ω 1 ⊂ R d , Ω 2 ⊂ R m with d, m, n ∈ N (see Lemma A.7 inWe can further bound this term using that W max is the maximal weight of f and the width of the hidden layers is lower bounded by W ≥ d.∥ f ∥ W k,∞ ([-1,1] d ) ≤ 16 2 (e 2 k 4 d 2 ) 2k (2k) 2k(k+1) W 2k+1 max W 3k+1 d k = 2 O(k(k log(k)+log(dW Wmax))) . (C.57)Now we have all the necessary tools in order to derive a bound on the generalization error for our neural network model of the form given in Definition 6. In our proof, we first bound the prediction error in terms of functions with 2 m-dimensional domain and on which we can directly apply the Koksma-Hlawka inequality. Then, we use the previous result to obtain an explicit bound. Due to the regularity of the parameters α P and our model parameters w P , this prediction error bound is independent of the system size n.Before stating the formal result bounding the prediction error, we introduce some notation. We define the prediction error of a neural network f Θ,w with weights given by Θ, w asR(Θ) ≜ E x∼U [-1,1] m |f Θ,w (x) -tr(Oρ(x))| 2 , (C.58)where in our case, x ∼ U [-1, 1] m denotes x sampled from a uniform distribution over [-1, 1] m . We suppress w in the notation to avoid cluttering. For a neural network f Θ,w generated from training on some data {(x ℓ , y ℓ )} N ℓ=1 , we can define the training error asR(Θ) ≜ 1 N N ℓ=1 |f Θ,w (x ℓ ) -y ℓ | 2 . (C.59)Moreover, as in our analysis in Appendix C.1, we rely on an approximation of the ground state property tr(Oρ(x)) by a sum of smooth local functions P α P f P (x) (Lemma 2). Namely, combining Lemma 2 and Theorem 8, we have that for ϵ 1 > 0, then choosing δ 1 > 0 as in Equation (A.5), i.e., δ 1 = O(log 2 (1/ϵ 1 )),P α P f P (x) -tr(Oρ(x)) ≤ ϵ 1 2 (C.60)Note that here, again, we slightly alter the definition from Appendix A.1.2, where we do not include the coefficient α P in the definition of f P (x). With these definitions, we have the following guarantee on the prediction error. Lemma 10 (Prediction error bound). Let 1/e > ϵ 1 , ϵ 2 > 0. Consider a tanh neural network f Θ,w : [-1, 1] m → R with architecture defined in Definition 7 with weights Θ i ≤ W max for some W max > 0 independent of the system size n and weights w in the last layer. Suppose we train f Θ,w on data {(x ℓ , y ℓ )} N ℓ=1 of size N , where the x ℓ 's form a low-discrepancy sequence with star-discrepancy D * N and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Then, we haveR(Θ) ≤ R(Θ) + ϵ 2 1 2 + ϵ 2 2 + (∥w∥ 1 + ∥w∥ 2 1 )D * N • 2 O(polylog(W Wmax/ϵ1)) . (C.61)where m = |I P | = O(log 2 (1/ϵ 1 )) for I P defined in Equation (A.2).Proof. Recall the definition of our neural network model in Definition 6. In particular, our model is given by a function f Θ,w : [-1, 1] m → R defined by f Θ,w (x) = P ∈S (geo) w P f θ P P (x), (C.62)where we refer to f θ P P : [-1, 1] m → R as the local models. For f P (x) = tr(P ρ(χ P (x))) as considered in Equation (C.60), we can define the following quantities. Define the training error with respect to this local approximation byRloc (Θ) ≜ 1 N N ℓ=1 f Θ,w (x ℓ ) - P α P f P (x ℓ ) 2 (C.63)Also define the prediction error with respect to the local approximation asR loc (Θ) ≜ E x∼U [-1,1] m f Θ,w (x) - P α P f P (x) 2 , (C.64)where x ∼ U [-1, 1] m means that x is sampled according to the uniform distribution.By Lemma 2, for our choice of δ 1 , we have Equation (C.60):P α P f P (x) -tr(Oρ(x)) ≤ ϵ 1 2 . (C.65)By the triangle inequality, we can bound the prediction error asR(Θ) = E x∼U [-1,1] m f Θ,w (x) - P α P f P (x) + P α P f P (x) -tr(Oρ(x)) 2 ≤ R loc (Θ) + ϵ 2 1 4 .(C.66)By applying the reverse triangle inequality, we can further bound this asR(Θ) ≤ ϵ 2 1 4 + Rloc (Θ) + |R loc (Θ) -Rloc (Θ)| (C.67) = ϵ 2 1 4 + Rloc (Θ) (C.68)+ E x∼U [-1,1] m f Θ,w (x) - P α P f P (x) 2 - 1 N N ℓ=1 f Θ,w (x ℓ ) -α P f P (x ℓ ) 2 (C.69) = ϵ 2 1 4 + Rloc (Θ) (C.70) + E x∼U [-1,1] m P w P f θ P P (x) -α P f P (x) 2 - 1 N N ℓ=1 P w P f θ P P (x ℓ ) -α P f P (x ℓ ) 2 (C.71)We can expand the term in the expectation/sum as followsP w P f θ P P (x) -α P f P (x) 2 (C.72) = P w P f θ P P (x) 2 -2 P w P f θ P P (x) P α P f P (x) + P α P f P (x) 2 (C.73) = P1,P2 w P1 f θ P 1 P1 (x)w P2 f θ P 2 P2 (x) -2w P1 f θ P 1 P1 (x)α P2 f P2 (x) + α P1 f P1 (x)α P2 f P2 (x). (C.74)Plugging this into the absolute value term in Equation (C.71) and upper bounding it with the triangle inequality, we have|R loc (Θ) -Rloc (Θ)| ≤ P1,P2 |w P1 ||w P2 | E x∼U [-1,1] m [f θ P 1 P1 (x)f θ P 2 P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f θ P 2 P2 (x ℓ ) + 2|w P1 |α P2 | E x∼U [-1,1] m [f θ P 1 P1 (x)f P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f P2 (x ℓ ) + |α P1 ||α P2 | E x∼U [-1,1] m [f P1 (x)f P2 (x)] - 1 N N ℓ=1 f P1 (x ℓ )f P2 (x ℓ ) . (C.75)Notice that in the expectation over [-1, 1] m , we can replace this by an expectation over the set of local parameters, i.e., the parameters with coordinates in I P1 ∪ I P2 , which we denote by S P1,P2 . This is because the functions in the expectations are local functions that only depend on these local parameters. The dimension of the domain we integrate over thus becomes independent of the system size n, as |S P1,P2 | ≤ 2 m = 2|I P |.We can now bound this term further using the Koksma-Hlawka inequality (Theorem 10). We apply a simple variable transformation τ (x) = 2x -1 so that the domain of f P • τ becomes [0, 1] m. Furthermore, we denote the domain associated with S P1,P2 by Ω P1,P2 ≜ [0, 1] |S P 1 ,P 2 | . Starting with the first term in Equation (C.75), we obtainE x∼U [0,1] m [f θ P 1 P1 (τ (x))f θ P 2 P2 (τ (x))] - 1 N N ℓ=1 f θ P 1 P1 (τ (x ℓ ))f θ P 2 P2 (τ (x ℓ )) (C.76) = E x∼U (Ω P 1 ,P 2 ) [f θ P 1 P1 (τ (x))f θ P 2 P2 (τ (x))] - 1 N N ℓ=1 f θ P 1 P1 (τ (x ℓ ))f θ P 2 P2 (τ (x ℓ )) (C.77) = S P 1 ,P 2 f θ P 1 P1 (τ (x))f θ P 2 P2 (τ (x))dx - 1 N N ℓ=1 f θ P 1 P1 (τ (x ℓ ))f θ P 2 P2 (τ (x ℓ )) (C.78) ≤ D * N (2 m)V HK (f θ P 1 P1 • f θ P 2 P2 ) • τ , (C.79)where xℓ = τ -1 (x ℓ ), such that Equation (C.76) matches the expression referenced in Equation (C. 75) . Note that we also applied in the last step that the star-discrepancy is increasing with respect to the dimension of the sequence. By application of the chain rule and the Cauchy-Schwartz inequality in the definition of the Hardy-Krause variation (Equation (A.36)), it is easy to see thatV HK (f θ P 1 P1 • f θ P 2 P2 ) • τ ≤ 2 2 mV HK (f θ P 1 P1 • f θ P 2 P2 ). (C.80)For all subsets S ′ ⊆ S P1,P2 , applying the product rule yields∂ |S ′ | ∂x S ′ (f θ P 1 P1 • f θ P 2 P2 ) ≤ A⊆S ′ ∂ |A| ∂x A f θ P 1 P1 ∂ |S ′ \A| ∂x S ′ \A f θ P 2 P2 = 2 O( m log(W Wmax)+ m2 log( m)) , (C.81)where the last equality follows from applying Lemma 9 from Appendix C.4 with d = k = 2 m and |{A : A ⊆ S ′ }| = 2 2 m. Here, we are using the notation ∂ |B| ∂x B to denote the mixed derivative with respect to all parameters x i ∈ B for some set B. Thus, applying Lemma 9 again, we obtainV HK (f θ P 1 P1 • f θ P 2 P2 ) ≤ S ′ ⊆S P 1 ,P 2 ∂ |S ′ | ∂x S ′ (f θ P 1 P1 • f θ P 2 P2 ) = 2 O( m log(W Wmax)+ m2 log( m)) . (C.82)Thus, putting it all together, we see thatE x∼U [-1,1] m [f θ P 1 P1 (x)f θ P 2 P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f θ P 2 P2 (x ℓ ) (C.83) ≤ 2 2 mD * N (2 m)2 O( m log(W Wmax)+ m2 log( m)) . (C.84)The remaining terms in Equation (C.75) can be bounded similarly using Lemma 22 from Appendix C.4. This lemma is applicable to f P because the derivatives are with respect to the local parameters. In this way, we can upper bound Equation (C.75) by|R loc (Θ) -Rloc (Θ)| ≤ P1,P2 (|w P1 ||w P2 | + |w P1 ||α P2 |)2 O( m log(W Wmax)+ m2 log( m)) (C.85) +|α P1 ||α P2 |2 O( m log( m)) D * N (2 m). (C.86)Plugging this back in to Equation (C.67), we haveR(Θ) ≤ ϵ 2 1 4 + Rloc (Θ) + P1,P2 (|w P1 ||w P2 | + |w P1 ||α P2 |)2 O( m log(W Wmax)+ m2 log( m)) (C.87) +|α P1 ||α P2 |2 O( m log( m)) D * N (2 m). (C.88)Lastly, we can bound Rloc (Θ) ≤ ϵ 2 1 /4 + ϵ 2 2 + R(Θ) in the same way as in Equation (C.66): Rloc (Θ) = 1 N N ℓ=1 f Θ,w (x ℓ ) - P α P f P (χ P (x ℓ )) 2 (C.89) ≤ 1 N N ℓ=1 f Θ,w (x ℓ ) -y ℓ 2 + |y ℓ -tr(Oρ(x ℓ ))| 2 + tr(Oρ(x ℓ )) - P α P f P (χ P (x ℓ ))R(Θ) ≤ R(Θ) + ϵ 2 1 2 + ϵ 2 2 + (∥w∥ 1 + ∥w∥ 2 1 )D * N (2 m)2 O(polylog(W Wmax/ϵ1)) . (C.92)Using the previous result and the results from low-discrepancy theory (see Appendix A.2 for a review), we can now show that Algorithm 1 will, under mild assumptions for training, output a model, which yields low prediction error. Thus, using Lemma 10, we can easily prove Theorem 14.Proof of Theorem 14. By Theorem 9, we know that for Sobol sequences in base 2 with points in [0, 1] d , the star-discrepancy is bounded byD * N (d) ≤ C(d) log(N ) d N , (C.93)where C(d) is a constant such thatC(d) < 1 d! d log(2d) . (C.94) Since C(d) = o(1), there exists a constant C, such that C ≥ C(d) for all d > 0. In our case, d = 2 m, so we haveD * N (2 m) ≤ C log(N ) 2 m N . (C.95)Using the assumption that the training objective is not larger than ((ϵ 1 + ϵ 2 ) 2 + ϵ 3 )/2, by Lemma 10, we haveR(Θ * ) = E x∼U [-1,1] m |f Θ * ,w * (x) -tr(Oρ(x))| 2 (C.96) ≤ ϵ 2 1 2 + ϵ 2 2 + (ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 + C ′ log(N ) polylog(1/ϵ1) 2 O(polylog(Wmax/ϵ1)) N (C.97) ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 + C ′ log(N ) polylog(1/ϵ1) 2 O(polylog(Wmax/ϵ1)) N , (C.98)where C ′ is a constant. We also used here that m = |I P | = O(polylog(1/ϵ 1 ). Since the training data has size N = O 2 polylog(1/ϵ1)+polylog(1/ϵ3) , W max can be chosen with respect to ϵ 1 , ϵ 3 and independent of the system size n such thatC ′ log(N ) polylog(1/ϵ1) 2 O(polylog(Wmax/ϵ1)) N ≤ ϵ 3 4 . (C.99)In this way, we obtainR(Θ * ) ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 . (C.100)Since the training objective from Definition 7 is non-convex, we cannot guarantee that our algorithm converges to a neural network with low training error. However, the assumptions made in Theorem 14 are rather mild in practice. Small training errors are a well-known phenomenon in deep learning and usually come at the expense of a larger prediction error, which is referred to as overfitting. Overfitting may arise due to excessive model complexity [89] , i.e. too many trainable parameters. This is reflected by Lemma 10, since the generalization error increases with the width W of the layers. The major challenge in practice lies in finding an appropriate balance between achieving a small training objective and model complexity, rather than only the latter. Furthermore, when the inputs are regularized, the weights usually remain small during training when initialized properly. This was for example observed in [53] .Finally, it is worth noting that in a scenario with a constant number of parameters m = O(1), similar to the setup in [51] , the expression derived from the outcome in Lemma 10 exhibits nearly linear dependence on ϵ. When incorporating the constant number of parameters by setting m = m, we recover the exact ground state properties tr(P ρ(x)) in f P . Thus, ϵ 1 in Lemma 10 becomes 0. Hence, the ability of LDS training to overcome the curse of dimensionality can unfold its full potential, since the domain dimension becomes independent of ϵ and expression Equation (C.88) reduces to a constant multiplied by D * N (2m). By Equation (C.95), we obtain R(Θ) = O(ϵ -(1+δ) ) for any δ > 0 and ϵ small enough, when the conditions of Theorem 5 are fulfilled.In this section, we generalize our results to hold for a wider class of distributions. Recall that our rigorous guarantee proven so far (Theorem 14) holds when the training data is generated according to a low-discrepancy sequence and the prediction error is measured with respect to the uniform distribution. We want to extend this result for different choices of both training and prediction error distributions. Notice that our prediction error bound (Lemma 10) is the only place that requires these assumptions on the distributions. Thus, in this section, we establish bounds on the expected prediction error for a more general family of distributions. We consider the following two cases.1. The training data is generated according to a general low-discrepancy sequence (in the sense of Definitions 4 and 5), and the prediction error is measured with respect to some distribution D.The training data consists of independently and identically distributed (i.i.d.) random samples according to a distribution D, and the prediction error is measured with respect to the same distribution D.There are some conditions on the distributions that we discuss shortly. In Case 1, suppose for example that we want to provide rigorous guarantees on the prediction error when the parameters x ∈ [-1, 1] m are sampled from a standard normal distribution (restricted to [-1, 1] m and normalized appropriately). As normally distributed test samples are more densely populated around the mean and more sparse around the boundary of the input domain, we need to predict more accurately around the mean than close to the boundary. When using a uniform low-discrepancy sequence for training, as in Algorithm 1, the predictive capabilities of our model are not exploited properly. To remedy this, we consider the training data to form a general low-discrepancy sequence, where it is low-discrepancy with respect to a normal distribution. We can relate this general low-discrepancy sequence to an LDS with respect to the Lebesgue measure, which are the sequences considered in Appendix C.2, via the probability integral transform (see, e.g., [90] ). We sometimes refer to LDS with respect to the Lebesgue measure as uniform low-discrepancy sequences. Formally, for any random variable X, which follows some probability distribution P (X ≥ x) ≜ F X (x), the random variable Y = F X (X) follows a uniform distribution. It turns out that the same transformation on LDS produces LDS with respect to other measures than the Lebesgue measure, as illustrated in Figure 3 . Moreover, under some assumptions on the distribution, we can bound the discrepancy with respect to other measures in terms of the discrepancy with respect to the Lebesgue measure, which we know how to bound as in Appendix C.2.In the following, we formalize this argument and adapt it to our problem setting. We refer the reader to Appendix A.2 to review the necessary concepts of generalized (star-)discrepancy, the Koksma-Hlawka inequality, and related results. Then, we demonstrate that a generalization of Lemma 10 and Theorem 14 can be achieved by incorporating these findings with slight adjustments to the proofs.In Case 2, we consider training data sampled i.i.d. from some distribution D and prediction error measured with respect to the same distribution D. To obtain a rigorous guarantee on the prediction error in this case, we leverage a probabilistic bound on the discrepancy of uniformly random points from [82] . Utilizing the previously established framework from Case 1, we can bound the discrepancy of points sampled from D in terms of the discrepancy of uniformly random points. This allows us to establish similar guarantees for Case 2.Before proving each of these cases, we set up our probabilistic framework and define the Borel measure with respect to which our low-discrepancy sequence is defined. Let g ≜ PDF(D) be the probability density function (PDF) of the data distribution and let G ≜ CDF(D) be the corresponding cumulative distribution function (CDF). In the following, assume that the PDF g satisfies the following properties. (c) Component-wise independence: The (random) variables ⃗ x i , ⃗ x j upon which different local terms h i (⃗ x i ), h j (⃗ x j ) of the Hamiltonian depend on, are independent. Hence, the PDF g is of the formg(x) = L j=1 g j (⃗ x j ) (C.101)for PDFs g j .We implicitly assume that g also satisfies all properties of a probability density function. It should be noted that Assumptions (a), (b) could technically be relaxed. We expand more on this later. Notice that if g : [-1, 1] m → R meets these requirements, the same holds for ḡ ≜ g • τ : [0, 1] m → R. Here, we use the notation from the previous section, where a bar denotes that we are working in the domain [0, 1] m as opposed to [-1, 1] m , and τ ( (x)) = 2x -1. Since the available results hold on [0, 1] m , we will mostly work with ḡ and the corresponding CDF Ḡ.We continue to set up the necessary notation to formally state our prediction error bound for Case 1.Let S P1,P2 , Ω P1,P2 be as in the proof of Lemma 10. Namely, let S P1,P2 be the parameters with coordinates in I P1 ∪ I P2 , where I P is defined in Equation (A.2), and let Ω P1,P2 = [0, 1] |S P 1 ,P 2 | . Additionally, define µ P1,P2 ≜ j∈S P 1 ,P 2 Ḡj ( ⃗ xj ) as the probability measure of the marginal for all (random) variables with indices in S P1,P2 . Due to Assumption (c), µ P1,P2 depends on at most 2 m variables. Furthermore, we define µ * ≜ arg maxµ P i ,P j D N (|S Pi,Pj |; µ Pi,Pj ), (C.102)and denote by S * the corresponding coordinate set. S * forms the domain of µ * , and we use d * to denote the dimension of the domain.In both Case 1 and Case 2, the idea is to define a transformation F that maps random variables with an arbitrary distribution to uniformly random variables. Namely, we construct a mapping ϕ such thatE x∼D [u(x)] = E x∼U [-1,1] m [u(ϕ(x))] (C.103)for any function u. In the following, we introduce the transform F ≜ ϕ -1 , as has been introduced in [91, 88, 92] . F can nicely be characterized using ḡ and Ḡ, and assumptions on F are easy to verify for a given data distribution. In fact, if F satisfies a Lipschitz condition, then known results [88] bound the discrepancy with respect to an arbitrary measure in terms of the discrepancy with respect to the Lebesgue measure, i.e., we can directly upper-bound D N (d * ; µ * ) in terms of D N (d * ). Our prediction error bound for more general distributions follows from this result and the results from Appendix C.2.Let g * be defined such that dµ * (x) = g * (x)dx. Also, let A, B ⊆ S * be such that A ∩ B = ∅ and C = S * \ (A ∪ B). Then, we define the conditional marginal PDF asg * (x A |X B = x B ) ≜ [0,1] |C| g * (x)dx C [0,1] |A|+|C| (x B )1 0 • • • (x B ) |B| 0 g * (x)dx (C.104)and the corresponding CDF asG * (X A = x A |X B = x B ) ≜ (x A )1 0 • • • (x A ) |A| 0 g * (x A |x B )dx A . (C.105)For convenience, we refer to the indices of x in S * via x 1 , x 2 , . . . , x d * . We can do this without loss of generality by permuting the order of the parameters. Using these definitions, we can now define the reverse transformation asF : [0, 1] d * → [0, 1] d *, where the indices of F are given byF j (x) ≜ G * (X j = x j |X 1 = x 1 , . . . , X j-1 = x j-1 ). (C.106) If random variables are distributed as X ∼ G * , then F (X) ∼ U [-1, 1] d * (or equivalently U [0, 1] d * under the variable transformation τ (x) = 2x -1), since X 1 , X 2 |X 1 , . . . , X d * |X 1 , . . . , X d * -1 are independent and j F j (X) = G * (X). (C.107)Finally, with this notation set up, we can formally state our result for Case 1. Corollary 6 (Neural network sample complexity guarantee; generalization of Theorem 14). Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0, and let D be a distribution with PDF g fulfilling assumptions (a)-(c) and F according to Equation (C.106). Let f Θ * ,w * : [-1, 1] m → R be a neural network model produced from Algorithm 1 trained on data {(x ℓ , ŷℓ )} N ℓ=1 of sizeN = 2 O(polylog(1/ϵ1)+polylog(1/ϵ3)) , (C.108)where the x ℓ 's form a low-discrepancy Sobol sequence,xℓ = F -1 (x ℓ ) and |ŷ ℓ -tr(Oρ(x ℓ ))| ≤ ϵ 2 .Suppose that f Θ * ,w * achieves a training error of at most((ϵ 1 + ϵ 2 ) 2 + ϵ 3 )/2. Additionally, suppose that all parameters Θ * i of f Θ * ,w * satisfy |Θ * i | ≤ W max ,for some W max > 0 that is independent of the system size n. Then the neural network f Θ * ,w * achieves prediction errorE x∼D |f Θ * ,w * (x) -tr(Oρ(x))| 2 ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 . (C.109)Similarly, we also have a guarantee for Case 2, which is the version we state in the main text and the beginning of this appendix. Corollary 7 (Neural network sample complexity guarantee; generalization of Corollary 6 for random data). Let 1/e > ϵ 1 , ϵ 2 , ϵ 3 > 0, D a distribution with PDF g satisfying assumptions (a)-(c). Letf Θ * ,w * : [-1, 1] m → R be a neural network model produced from Algorithm 1 trained on data {(x ℓ , y ℓ )} N ℓ=1 of size N = log(1/δ)2 O(polylog(1/ϵ1)+polylog(1/ϵ3)) , (C.110)where the x ℓ ∼ D and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Suppose that f Θ * ,w * achieves a training error of at most ((ϵ 1 + ϵ 2 ) 2 + ϵ 3 )/2. Additionally, suppose that all parameters Θ * i of f Θ * ,w * satisfy |Θ * i | ≤ W max , for some W max > 0 that is independent of the system size n. Then the neural network f Θ * ,w * achieves prediction errorE x∼D |f Θ * ,w * (x) -tr(Oρ(x))| 2 ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 , (C.111)with probability at least 1δ.We first prove Corollary 6 similarly to how we proved Theorem 14. In particular, we can prove a generalized version of Lemma 10, where we are given a low-discrepancy sequence with respect to µ * and wish to bound the prediction error with respect to D, as in Case 1. Define the prediction error of a neural network f Θ,w with weights given by Θ, w with respect to a distribution D asR D (Θ) ≜ E x∼D |f Θ,w (x) -tr(Oρ(x))| 2 = [-1,1] m |f Θ,w (x) -tr(Oρ(x))| 2 dG(x), (C.112)where x ∼ D denotes x sampled from the distribution D over [-1, 1] m and dG(x) = g(x)dx. Again, we suppress w in the notation to avoid cluttering. Then, we have the following lemma. Lemma 11 (Generalized prediction error bound). Let 1/e > ϵ 1 , ϵ 2 > 0. Consider a tanh neural network f Θ,w : [-1, 1] m → R with architecture defined in Definition 7 with weights Θ i ≤ W max for some W max > 0 independent of the system size n and weights w in the last layer. Assume that G satisfies assumptions (a)-(c) and |y ℓtr(Oρ(x ℓ ))| ≤ ϵ 2 . Furthermore, suppose we train f Θ,w on data {(x ℓ , y ℓ )} N ℓ=1 of size N , where the τ -1 (x ℓ )'s from a set with star-discrepancy at mostD * N (d; µ * ) in each dimension d. Then, we have R D (Θ) ≤ R(Θ) + ϵ 2 1 2 + ϵ 2 2 + (∥w∥ 1 + ∥w∥ 2 1 )D * N (d * ; µ * ) • 2 O(polylog(W Wmax/ϵ1)) . (C.113) Moreover, if there exist constants b 1 , b 2 such that D * N (d) ≤ b 1 b 2 + log(1/δ) dN with probability at least 1δ, then there exists a constant b1 such thatR D (Θ) ≤ R(Θ)+ ϵ 2 1 2 +ϵ 2 2 +(∥w∥ 1 +∥w∥ 2 1 ) b1 1 + log(1/δ) m N •2 O(polylog(W Wmax/ϵ1)) (C.114)with probability at least 1δ.This lemma can be proven in the same fashion as Lemma 10 with two minor adjustments. One change is the use of a generalization of the Koksma-Hlawka inequality, which we discuss in Theorem 11 in Appendix A.2. To prove the second part of Lemma 11, we need a technical lemma (Lemma 14) to handle the probability of failure in the bound on the star-discrepancy. We relegate this to the end of the section, as it is mainly a technicality.Proof. We proceed in the same way as in Lemma 10, replacing R(Θ) with R D (Θ) and replacingR loc (Θ) with R loc,D (Θ) ≜ E x∼D f Θ,w (x) - P α P f P (x) 2 . (C.115)We follow the proof of Lemma 10 until Equation (C.75). This gives usR D (Θ) ≤ ϵ 2 1 4 + Rloc (Θ) + |R loc,D (Θ) -Rloc (Θ)|, (C.116)where recall thatRloc (Θ) = 1 N N ℓ=1 f Θ,w (x ℓ ) - P α P f P (x ℓ ) 2 , (C.117)as in Equation (C.63). Moreover, we also have the adjusted version of Equation (C.75)|R loc,D (Θ) -Rloc (Θ)| ≤ P1,P2 |w P1 ||w P2 | E x∼D [f θ P 1 P1 (x)f θ P 2 P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f θ P 2 P2 (x ℓ ) + 2|w P1 ||α P2 | E x∼D [f θ P 1 P1 (x)f P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f P2 (x ℓ ) + |α P1 ||α P2 | E x∼D [f P1 (x)f P2 (x)] - 1 N N ℓ=1 f P1 (x ℓ )f P2 (x ℓ ) . (C.118)To bound the first term, we use the generalized Koksma-Hlawka inequality (Theorem 11) to obtainE x∼D [f θ P 1 P1 (x)f θ P 2 P2 (x)] - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f θ P 2 P2 (x ℓ ) (C.119) = [-1,1] m |f Θ,w (x) -tr(Oρ(x))| 2 dG(x) - 1 N N ℓ=1 f θ P 1 P1 (x ℓ )f θ P 2 P2 (x ℓ ) (C.120) = [0,1] m f θ P 1 P1 (τ (x))f θ P 2 P2 (τ (x)) L i=1 ḡi ( ⃗ xi )dx - 1 N N ℓ=1 f θ P 1 P1 (τ (x ℓ ))f θ P 2 P2 (τ (x ℓ )) (C.121) = Ω P 1 ,P 2 f θ P 1 P1 (τ (x))f θ P 2 P2 (τ (x))dµ P1,P2 (x) - 1 N N ℓ=1 f θ P 1 P1 (τ (x ℓ ))f θ P 2 P2 (τ (x ℓ )) (C.122) ≤ D * N (d * ; µ * )V HK (f θ P 1 P1 • f θ P 2 P2 ) • τ . (C.123)Here, in the first equality, we use Assumption (c) on the structure of the PDF g and use the variable transformation τ (x) = 2x -1. In the second equality, similarly to Lemma 10, we notice that because the functions in the expectation are local functions that only depend on parameters in S P1,P2 , we can replace the expectation over the whole domain [0, 1] m with an expectation over just the domain Ω P1,P2 . This step also crucially uses Assumption (c), where the factorization of the PDF g due to independence is needed. The last line uses the generalized Koksma-Hlawka inequality (Theorem 11). The remainder of the proof follows in the same way as Lemma 10.The second part of the statement is a direct consequence of Lemma 14. Specifically, following the proof of Lemma 10, we use Lemma 14 to bound the term in Equation (C.88). This is necessary because the upper bound on the star-discrepancy only holds probabilistically, so we must show that we can still use this upper bound on the sum of several star-discrepancy terms. This is more complicated than a simple union bound, and we relegate the proof and statement of Lemma 14 to the end of this section.In addition to this lemma, we also need a result from [88] , adapted to our definitions above. In the following, we use D N (ω; d) to denote the discrepancy with respect to the Lebesgue measure of a specific sequence ω of length N and dimension d, as in Definition 1. Similarly, we use D N (ω; d; µ) to denote the discrepancy with respect to a measure µ of a specific sequence ω of length N and dimension d, as in Definition 4. Lemma 12 (Theorem 2 in [88] ). Let ω = {x ℓ } N ℓ=1 be an arbitrary sequence on the open ddimensional unit cube with discrepancy D N (ω, d), and let ω = {x ℓ } N ℓ=1 be the sequence defined by F xℓ = x ℓ , where F is defined in Equation (C.106). Moreover, let g be a strictly positive, d-times continuously differentiable PDF, such that g(x) ≥ m > 0 for all x. Let G be the corresponding probability measure (i.e., CDF). Furthermore, let F satisfy∥F (x) -F (y)∥ ≤ K∥x -y∥. (C.124)Then, the discrepancy of ω with respect to G is bounded asD N (ω; d; G) ≤ c (D N (ω; d)) 1 d , (C.125) where c = 2d • 3 d (K + 1) d-1 .The authors in [88] note that the assumption on F in Equation (C.124) is certainly fulfilled when g is continuously differentiable. This is where Assumption (b) is used, where technically, we only require this Lipschitz condition on F . With Lemma 11 and Lemma 12, we are ready to prove Corollary 6.Proof of Corollary 6. We proceed similarly as in the proof of Theorem 14 but this time using Lemma 11 instead of Lemma 10. First, we bound the nonuniform discrepancy of our training inputs ω = {x ℓ } N ℓ=1 . Recall that xℓ = F -1 (x ℓ ) for x ℓ generated according to a low-discrepancy Sobol sequence (i.e., low-discrepancy with respect to the Lebesgue measure). By definition of F , then ω has star-discrepancy D * N (ω; d * ; µ * ). By Assumption (b), we can apply Lemma 12 to obtain D * N (ω; d * ; µ * ) ≤ D N (ω; d * ; µ * ) ≤ c (D N (ω; d * )) 1 d * ≤ 2 d * c (D * N (ω; d * ))R D (Θ * ) = E x∼D |f Θ * ,w * (x) -tr(Oρ(x))| 2 (C.129) ≤ ϵ 2 1 2 + ϵ 2 2 + (ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 + C ′ log(N )2 O(polylog(Wmax/ϵ1)) N 1/d * (C.130) ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 2 + C ′ log(N )2 O(polylog(Wmax/ϵ1)) N 1/polylog(1/ϵ1) , (C.131)where C ′ is a constant. We also used here that d * ≤ 2 m and m = |I P | = O(polylog(1/ϵ 1 ). Since the training data has size N = O 2 polylog(1/ϵ1)+polylog(1/ϵ3) , W max can be chosen with respect to ϵ 1 , ϵ 3 and independent of the system size n such thatC ′′ log(N )2 O(polylog(Wmax/ϵ1)) N 1/polylog(1/ϵ1) ≤ ϵ 3 4 (C.132)for some constant C ′′ . In this way, we obtainR D (Θ * ) ≤ 2(ϵ 1 + ϵ 2 ) 2 + ϵ 3 . (C.133)Finally, we prove Case 2, where the training data is sampled i.i.d. according to a distribution D and the prediction error is also measured with respect to D. Note that we can drop the assumption that F -1 is efficiently computable for this case. The key result we need for this is a bound on the star-discrepancy for uniformly random points (Lemma 5). with probability at least 1δ. The rest of the proof follows in the same way as Corollary 6, but using the above discrepancy bound. Note that because this discrepancy bound only holds probabilistically, we need to use the second part of Lemma 11.Our prediction error bounds for Case 1 and Case 2 (in Corollaries 6 and 7, respectively) look rather similar. However, one can verify that Corollary 6 only requires about square root of the number of samples Corollary 7 uses to achieve a certain risk bound with small enough ϵ, but this advantage is hidden in the polylogarithmic factors in the exponent. Hence, low-discrepancy data yields better theoretical guarantees. However, they did not yield an improvement in our numerical experiments, as discussed in Appendix D. The size N of the training set seems to be very large for low-discrepancy data to have practical effects. In the main text, we present Corollary 7 because it is the more general theoretical statement.In fact, one can also improve the polylogarithmic factors in Corollary 6 by imposing stronger assumptions on the distribution. In particular, the result in Lemma 12 seems surprisingly weak at first glance. Multidimensional transformations do, however, constitute a major challenge, since they generally do not preserve properties such as lines remaining straight or parallel. The boxes over which one optimizes in order to compute the discrepancy can thus change severely in shape, which can strongly alter the discrepancy and makes it difficult to analyze. This can result in a rather poor scaling in terms of the discrepancy with respect to the Lebesgue measure and thus N . However, when g fulfills additional assumptions, we obtain a much better dependence on N by directly applying the Koksma-Hlawka inequality (with respect to the Lebesgue measure) to f • ϕ. Unsurprisingly, this is possible when the mixed derivative of F -1 = ϕ is bounded on [0, 1]. This follows, when g's mixed derivative is bounded [88] . We restate this result below. Lemma 13 (Theorem 1 in [88] ). Let ω = {x ℓ } N ℓ=1 of size N be an arbitrary sequence on the open d-dimensional unit cube with discrepancy D N (ω, d) and ω = {x ℓ } N ℓ=1 the sequence defined by F xℓ = x ℓ , where F is defined in Equation (C.106). Moreover, let g be a strictly positive, d-times continuously differentiable PDF, such that g(x) ≥ m > 0 for all x. Let G be the corresponding probability measure (i.e., CDF). Furthermore, let F satisfy∂ |A| F j ∂x A ≤ M 1 ≤ j ≤ d, A ⊆ {1, . . . , d}. (C.136) Then [0,1 d ] f (x)dG(x) - 1 N N ℓ=1 f (x ℓ ) ≤ d! M m 2d-1 D * N (ω; d)V HK (f ). (C.137)On a high level, the proof works via the observation that the Jacobian of F has g as its determinant, which is strictly positive. Since F • ϕ is the identity, one can write the Jacobian of F • ϕ as a linear system of equations with the derivatives of ϕ as solution. Using Cramer's rule and the assumption on F , one can upper bound the derivative of ϕ. Applying this iteratively, one can show via induction that the mixed derivatives of ϕ are also bounded when the mixed derivatives of F are bounded. Note that (as stated in [88] ) when D is a product of independent distributions, i.e. g(x) = m i=1 g i (x i ) and fulfills assumptions (a)-(c), the conditions for Lemma 13 are also fulfilled. It is important to emphasize that the additional assumption used in Lemma 13 yield a much better dependence of ϵ on N . However, this improvement is hidden in the polylogarithmic factors.We dedicate the last part of this section to the proof the following statement, which we used in the proof of Lemma 11. At a high level, this shows that when we have a probabilistic upper bound on the star-discrepancy, we can still upper bound a sum of star-discrepancies with high probability. (c 1 |α P1 ||α P2 | + c 2 (|w P1 ||α P2 | + |w P1 ||w P2 |))D * N (x P1,P2 ) ≥ t   (C.138) ≤ exp - N t 2 b1 (c 1 ∥α∥ 2 1 + c 2 (∥w∥ 1 ∥α∥ 1 + ∥w∥ 2 1 )) 2|α P2 |D * N (x P1,P2 ) ≥ t   ≤ exp - N t 2 224∥α∥ 2 2 ( m) 2 exp(b 1 )b 2 (C.148)for any P 1 ∈ S (geo) and any t > 0.Proof. Let P 1 ∈ S (geo) . DefineX j ≜ E   P2∈S (geo) |α P2 |D * N (x P1,P2 ) Y j-1 , . . . , Y 0   , (C.149)where we omit the dependence x P1,P2 = x P1,P2 (Y 1 , . . . , Y m ) and Y j ≜ {(x j ) ℓ } N ℓ=1 and x j parameterize h j . We consider all increments, which are not contained in I P1 . Hence, with slight abuse of notation, let index j = 0 refer to all coordinates in I P1 and Y 0 ≜ {{(x j ) :j ∈ I P1 } ℓ } N ℓ=1 . Furthermore, let X 0 ≜ E   P2∈S (geo) |α P2 |D * N (x P1,P2 ) Y 0   (C.150)and Z 1 ≜ X 1 -X 0 . Clearly, X 0 , . . . , X m is a martingale sequence and Z 1 , . |α P2 |D * N (x P1,P2 ) ≥ t   ≤ exp - t 2 28b ′ m i=1 c 2 i (C.160) with b ′ = 2 m exp(b 1 ) and c 2 i = 4b2 m N P2∈S (geo) :i∈I P 2 |α P2 | 2 . Furthermore, m i=1 P2∈S (geo) :i∈I P 2 |α P2 | 2 = P2∈S (geo) |α P2 | 2 m i=1 1{i ∈ I P2 } = m P2∈S (geo) |α P2 | 2 = m∥α∥ 2 2 .(C.161) The result follows from this. Now, we are finally able to prove Lemma 14.Proof of Lemma 14. We need to bound the weighted sum of the star-discrepancies we consider. This requires an extra step, since the star-discrepancy may vary among the sequences in the sum. Note that simply applying the union bound would result in a log(n)-factor. Luckily, only the sum needs to be small, rather than all individual terms needing to be small at once. Recall that we use S P1,P2 to denote the set of parameters with coordinates in I P1 ∪ I P2 . In the following, we use )). For illustrative purposes, we only consider the first term for now. We proceed similarly to the proof of Lemma 17. DefineD * N (xX j ≜ E   P1,P2∈S (geo) |α P1 ||α P2 |D * N (x P1,P2 ) Y j-1 , . . . , Y 1   , (C.162)where the expectation is with respect to the inputs Y j ≜ {(x j ) ℓ } N ℓ=1 and x j parametrize h j . Furthermore, let In the last step, we used the observation that for j to be contained in I P1 ∪ I P2 , it has to be contained in at least one of the two sets. Hence, we can enumerate the admissible coordinate sets by fixing P 1 , such that I P1 contains j and combine it with all I P2 . The factor two arises from doing the same with P 1 when fixing P 2 . Now, for any t > 0, we have X 0 ≜ E   P1,P2∈S (geo) |α P1 ||α P2 |D * N (x S P 1 ,P 2 )   (C.163) and Z 1 ≜ X 1 -X 0 . Clearly, X 0 , . . . ,EQUATIONEQUATIONP1∈S (geo) :j∈I P 1 Pr   P2∈S (geo) |α P2 |D * N (x P1,P2 ) ≥ t 4|α P1 |   (C.168) ≤ 2 P1∈S (geo) :j∈I P 1 exp - N t 2 16 • 224|α P1 | 2 ∥α∥ 2 2 ( m) 2 exp(b 1 )b 2 (C.169) ≤ 2|{P 1 ∈ S (geo) : j ∈ I P1 }| exp - N t 2 16 • 224∥α∥ 2 2 ( m) 2 exp(b 1 )b 2 P1∈S (geo) :j∈I P 1 |α P1 | 2 (C.170) ≤ 2 m exp - N t 2 16 • 224∥α∥ 2 2 ( m)(c 1 |α P1 ||α P2 | + c 2 (|w P1 ||α P2 | + |w P1 ||w P2 |))D * N (x P1,P2 )). (C.173)Using that ∥α∥ 2 2 ≤ ∥α∥ 2 1 and solving for the appropriate δ yields the desired result.Let O = P ∈{I,X,Y,Z} ⊗n α P P be an observable that can be written as a sum of geometrically local observables. In the following, we derive an expression for the mixed partial derivatives of tr(P ρ(x)), using tools from the spectral flow formalism [59, 60, 61] . This allows us to bound the Hardy-Krause variation (Equation (A.36)) in Appendix C.2. Let the spectral gap of H(x) be lower bounded by some constant γ for all choices of parameters x ∈ [-1, 1] m . Then, by the spectral flow formalism [59, 60, 61] , the directional derivative of a ground state of H(x) in the direction defined by the parameter unit vector û is given by∂ ∂ û ρ(x) = û • ∇ x ρ(x) = i[D û(x), ρ(x)] (C.174)where D û(x) = +∞ -∞ W γ (t)e itH(x) ∂H ∂ û (x)e -itH(x)(γ|t|) = 1/2.This allows to us to obtain an expression of the first order derivative of ρ(x) with respect to some parameter x k . Consider the unit vector û = êk ≜ (0, . . . 0, 1, 0, . . . 0) T , where the 1 is in the kth position. Then, the directional derivative in the direction given by e k is∂ ∂ê k ρ(x) = êk • ∇ x ρ(x) = ∂ ∂x k ρ(x) = i[D êk (x), ρ(x)]. (C.177)Hence, we obtain∂ ∂x 1 tr(P ρ(x)) = tr P ∂ ∂x 1 ρ(x) = i tr(P [D ê1 (x), ρ(x)]) = i tr([P, D ê1 (x)]ρ(x)). (C.178)In order to compute the mixed derivative of second order, we now apply the product rule to this expression, which yields∂ 2 ∂x 1 ∂x 2 α P tr(P ρ(x)) = ∂ ∂x 2 iα P tr([P, D ê1 (x)]ρ(x)) (C.179) = iα P tr P, ∂ ∂x 2 D ê1 (x) ρ(x) -tr [P, D ê1 (x)] ∂ ∂x 2 ρ(x) (C.180) = α P tr i P, ∂ ∂x 2 D ê1 (x) ρ(x) -tr([[P, D ê1 (x)], D ê2 (x)]ρ(x)).(C.181)Note that the terms of this expression can be treated similarly as the first partial derivative. For each additional partial derivative, we obtain terms consisting of the product with nested commutators with ρ(x) under the trace. The nested commutators contain D êj (x) or partial derivatives of it, for which we will later derive an explicit form. Hence, we can apply the same scheme until we arrive at the k-th partial derivative. In order to formalize this statement, we need to introduce some additional notation.Throughout the rest of this section, we use the notation ∂ |B| ∂x B to denote the mixed derivative with respect to all parameters x i ∈ B for some set B.Definition 8. Let k ∈ N. Let A ⊆ [k] be a set of size |A| = m. Define an ordering l 1 < l 2 < • • • < l m over the elements l 1 , l 2 , . . . , l m ∈ A of A. Then, we define ‚ l∈A ∂ B l l ≜ i m tr • • • P, ∂ |B l 1 | ∂x B l 1 D êl 1 (x) , ∂ |B l 2 | ∂x B l 2 D êl 2 (x) • • • , ∂ |B lm | ∂x B lm D êlm (x) ρ(x) , (C.182) where B j ⊂ [k].We refer to the nested commutators under the trace as summands. The set A and collection {B l } l∈A ≜ B A satisfy the following conditions:1. Each summand contains D ê1 (x). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: When using data or parts of code from previous papers, this is cited. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.2-local in the sense that P does not act on more than two not necessarily neighboring sites. 
 
 Appendix:  


Decision: Regarding the potential impact of this paper on introducing rigorous machine learning techniques to quantum physics research, I recommend accepting it. Predicting ground state properties is a fundamental task in quantum physics. This paper is the first theoretical work on predicting ground state properties using deep learning and, modulo some assumptions that deserve further study, achieves state-of-the-art sample complexity results. Although there are several gaps between theory and practice in the proposed algorithms, potential concerns are explicitly raised and mostly discussed in the paper.

Below are some additional comments. 
- Some discussion is needed to justify the bounded weights assumption in Theorem 5. 
- The claim that "we cannot guarantee a priori that the network will indeed achieve a low training error" in Section 3.2 deserves further consideration. I agree with the claim, but I believe the reason is the agnostic nature of machine learning. 
- Please consider revising the paper to avoid potential concerns and misunderstandings as those raised by Reviewer u4Mp.

Reviews: review: paper_summary: paper_summary: The authors propose an ML based method to predict properties of ground states of quantum systems which comes with provable guarantees. Improving on recent work by Huang et al and Lewis et al, they give sample complexity bounds which are independent of the number of qubits. This approach is applicable when the observable one is trying to predict is predetermined. The authors also suggest a deep learning based approach for the case where the observable is not known in advance. They support their theoretical wok with numerical experiments.
strength_weakness: strength_weakness: It is unclear to me how the Neural Network generalization result compares to known results in the literature- the setting which the authors study is quite specific and thus it is not easy to relate the result they obtained to those in the theoretical deep learning literature.
questions: questions: How restrivtive is the assumption about the dependence of each local term on a fixed number of parameters? It would be instructive to give physically relevant examples where this does and does not hold.
Presumably the results would not apply to standard Neural Network architectures- what specifically would break in the analysis? 
How does the sample complexity depend on the constant the authors assume the network weights are bounded by?
limitations: limitations: I believe the authors have adequately adressed the limitations of the paper.
paper_summary: paper_summary: In this work, the authors give two algorithms that predict (geometrically local) properties of ground states of gapped geometrically local Hamiltonians. This problem has been introduced by Huang et al. [HKT+22], and the previous best known algorithm is given by Lewis et al. [LHT+24], which uses $\log(n)$ samples, where $n$ is the number of qubits in the Hamiltonian. This paper further improves on the $log(n)$ sample complexity, and gives two algorithms that only use a constant number of samples. The first  algorithm is modified from the algorithm of [LHT+24], changing the regression part of the algorithm from LASSO to kernel ridge regression. The second algorithm uses deep neural network, having the advantage of not needing to know the observables in advance, but requires more restriction on the distribution of the Hamiltonian parameters. The authors complement their theoretical results with numerical simulations.

[HKT+22] Huang, Hsin-Yuan, Richard Kueng, Giacomo Torlai, Victor V. Albert, and John Preskill. "Provably efficient machine learning for quantum many-body problems." Science 377, no. 6613 (2022): eabk3333
[LHT+24] Lewis, Laura, Hsin-Yuan Huang, Viet T. Tran, Sebastian Lehner, Richard Kueng, and John Preskill. "Improved machine learning algorithm for predicting ground state properties." nature communications 15, no. 1 (2024): 895.
strength_weakness: strength_weakness: The part of preliminaries that are restating definition and result of [LHT24+], is not well written, and I believe has led to a critical bug to the first algorithm. In particular, Theorem 8 claims that for every $O\sum_{P} \alpha_P P$ that can be written as sum of geometrically observables, $\sum_{P} |\alpha| =O(1)$. However, the counterpart in  [LHT24+] has extra restrictions: $||O||_{infty}=1$ and $O$ need to be inside a radius of $R=O(1)$.  Therefore, where the authors uses Theorem 8 in equation (B.28) to bound the kernel, the result is incorrect since they do not have $R=O(1)$.

Other minor inconsistencies includes: 

line 642: $S^{geo}$ not defined

line 660: $h_{c(j)}$ not defined
questions: questions: Some typos:

line 121: geometrically [local] observable

line 148: || \omega || -> || w ||
limitations: limitations: The authors adequately addressed the limitations.
paper_summary: paper_summary: This work builds upon the work of Huang et al. and Lewis et al. by introducing two new approaches to get constant sample complexity for predicting properties of a ground state of a many-body local Hamiltonian. The two new approaches are a modified ML model that requires knowledge of the property of interest and a deep neural network model that does not need prior knowledge of the property. In this paper, the authors provide both proves and small experimental evaluations to show that both approaches achieve constant sample complexity, independent of system size.
strength_weakness: strength_weakness: Even though it is a strong paper, the issue addressed here is a specific case that builds upon two other papers. Additionally, the related published works are mostly, if not all, published in physics journals. I do not see why the results shared in this paper are valuable to share with the broader NeurIPS community, especially since the mathematical proofs are very rigorous, I would expect that it is not accessible to the broader audience. Some assumptions and conditions required for the theoretical guarantees may also limit the applicability of the results.
questions: questions: How do the parameters or phases of the random Heisenberg model affect the training performance?
limitations: limitations: Limitations are addressed in the paper.
paper_summary: paper_summary: This paper studies the sample-efficient learnability of properties of grounds states of local Hamiltonians. Ground states of local Hamiltonians are hard to compute, even for quantum computers and to circumvent this hardness, several recent works proposed learning the trace inner product of local observables with the ground state given labeled training data. This setting is exactly PAC learning, i.e. given labeled data from a worst-case distribution, the goal is to get low prediction error wrt the same distribution on future samples. The best sample complexity for this problem is known to be log(n) 2^{polylog(1/eps)}, shown by Lewis et al. 

The main questions addressed in this work are (1) whether the sample complexity can be improved to be independent of the system size and (2) whether there are rigorous guarantees for learning properties of the ground state via neural network based algorithms.
strength_weakness: strength_weakness: I think the paper does not resolve the questions it claims to resolve and does so in a slightly camouflaged way.  

1. Question 1 in the paper asks whether you can get sample complexity that is independent of system size for learning properties of ground states, aka the PAC learning setting for ground states of local Hamiltonians. The answer obtained is yes, under two crucial caveats, the observable is known in advance and the distribution over the training data is not worst-case.  This diverges significantly from the PAC learning model. The same critique holds for Question 2. Further, reference [1] does not state this as an open problem. 

2. The assumptions on the training distribution are not stated upfront and do not appear to be mild. The assumptions include the distribution g to be strictly non-zero on [-1,1], and zero outside; g being continuously differentiable; and component-wise independent. Are there any natural distributions that satisfy all these properties simultaneously? 

3. There is no discussion of why each of these properties is needed, which ones are crucial to the argument and which ones are for technical convenience. Having 4 non-trivial technical assumptions of the pdf is a major weakness, especially since it makes the result incomparable to prior work [1,2], where the setting is truly PAC learning. 

4. In the numerical experiments, I see no discussion of what distribution was used to generate the training data, and how many of the technical conditions this distribution satisfies. It remains unclear to me what the experimental section is trying to convey, since it does not complement the main theorems. 

5. When is it reasonable to expect labeled data for the observable / property you want to learn? What is a real-world scenario where one would expect to obtain such labeled data? 

6. Is there a way to get some non-trivial sample complexity (not necessarily system independent) bound for PAC learning via neural networks, without the extra assumptions on the training distribution?  

7. There are completely unjustified claims such as the neural network achieving low loss and finding a bounded magnitude solution after constant many steps and O(n) time. It is not clear why this should ever be the case.
questions: questions: I included the questions with the weaknesses.
limitations: limitations: Yes.
paper_summary: paper_summary: In this paper, the authors focused on utilizing deep learning methods to predict the ground states. They made an important assumption that brings theoretical improvement to achieve constant sample complexity in the training data. They also made two main alternations to the learning model compared to previous literature, including incorporating Pauli coefficients in feature mapping and utilizing kernel ridge instead of Lasso. Numerical results for up to 45 qubit systems are provided, supporting the theoretical findings.
strength_weakness: strength_weakness: The main concern is that the improvement of this paper against precious works is limited. The main theoretical finding is based on an additional assumption that we know the property we'd like to predict in advance. The proposed learning method has only two minor alternations. These issues prevent me from giving a higher evaluation score, but they do not overshadow the fact that this article is of high quality.
questions: questions: No questions
limitations: limitations: See weakness
paper_summary: paper_summary: The paper presents two novel machine learning algorithms for predicting ground state properties of quantum systems with constant sample complexity, independent of system size. The first algorithm modifies an existing ML model, while the second introduces a deep neural network model, both showing improved scaling in numerical experiments.
strength_weakness: strength_weakness: 1. The training objective for the neural network is non-convex, which poses challenges in finding a global optimum efficiently. The paper does not address how to overcome this issue or guarantee convergence to optimal weights.

2. While the paper claims improved computational complexity, the actual implementation details and computational resources required for the deep learning model are not thoroughly discussed.
questions: questions: 1. The reviewer would appreciate if the authors could elaborate on how the performance of the deep learning model generalizes to Hamiltonians that extend beyond the specific cases examined in the numerical experiments.

2. Can the authors provide more insights into the practical implementation of their algorithms, particularly regarding the initialization and regularization procedures used during training? This will be helpful for readers reproduce the results of the paper.
limitations: limitations: The authors have discussed the limitations.

Comments: comment: Thanks to the authors for their response, and I have no further questions. I'd like to keep the score.
We thank 7tNg for acknowledging that their misunderstanding has been clarified and are happy that they are convinced of the soundness of the results. So, we wonder why the soundness score is still 1. Also, if they have any other issues we could address that leads to the overall score of only 5, we are happy to respond.
Thank you for your reply. I would like to keep my score.
I would like to thank the authors for the response. I have raised my score.
Dear Reviewer KHZF, 

The author-reviewer discussion period is ending soon. Please check if the authors’ response has addressed your concerns and feel free to adjust your score accordingly. If you find the authors’ response unsatisfactory, please explain your reasons and discuss them with the authors immediately.

Best regards,
AC
Dear Reviewer MeLw, 

The author-reviewer discussion period is ending soon. Please check if the authors’ response has addressed your concerns and feel free to adjust your score accordingly. If you find the authors’ response unsatisfactory, please explain your reasons and discuss them with the authors immediately.

Best regards,
AC
Dear Reviewer u4Mp, 

The author-reviewer discussion period is ending soon. Please check if the authors’ response has addressed your concerns and feel free to adjust your score accordingly. If you find the authors’ response unsatisfactory, please explain your reasons and discuss them with the authors immediately.

Best regards,
AC
Dear Reviewer 7tNg,

The author-reviewer discussion period is ending soon. Please check if the authors’ response has addressed your concerns and feel free to adjust your score accordingly. If you find the authors’ response unsatisfactory, please explain your reasons and discuss them with the authors *immediately*.

Best regards,
AC
I thank the authors for the response. I increased my rating for the paper.
I would like to thank the authors for their detailed response. I will leave my score unchanged.